source: trunk/Mars/msim/MPhotonData.cc@ 15658

Last change on this file since 15658 was 10060, checked in by rohlfs, 14 years ago
two new command line arguments of readcorsika: -A=arrayNo and -T=telescopeNo. New design of program flow in MCorsikaRead: It is now determined by the order of the data blocks in the input file.
File size: 9.7 KB
Line 
1/* ======================================================================== *\
2!
3! *
4! * This file is part of CheObs, the Modular Analysis and Reconstruction
5! * Software. It is distributed to you in the hope that it can be a useful
6! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
7! * It is distributed WITHOUT ANY WARRANTY.
8! *
9! * Permission to use, copy, modify and distribute this software and its
10! * documentation for any purpose is hereby granted without fee,
11! * provided that the above copyright notice appears in all copies and
12! * that both that copyright notice and this permission notice appear
13! * in supporting documentation. It is provided "as is" without express
14! * or implied warranty.
15! *
16!
17!
18! Author(s): Thomas Bretz, 12/2000 <mailto:thomas.bretz@epfl.ch>
19! Author(s): Qi Zhe, 06/2007 <mailto:qizhe@astro.uni-wuerzburg.de>
20!
21! Copyright: CheObs Software Development, 2000-2010
22!
23!
24\* ======================================================================== */
25
26/////////////////////////////////////////////////////////////////////////////
27//
28// MPhotonData
29//
30// Storage container to store Corsika events
31//
32// For details on the coordinate systems see our Wiki.
33//
34// Version 1:
35// ----------
36// * First implementation
37//
38// Version 2:
39// ----------
40// - fNumPhotons
41//
42/////////////////////////////////////////////////////////////////////////////
43#include "MPhotonData.h"
44
45#include <fstream>
46#include <iostream>
47
48#include <TMath.h>
49#include <TRandom.h>
50
51#include "MLog.h"
52#include "MLogManip.h"
53
54ClassImp(MPhotonData);
55
56using namespace std;
57
58// --------------------------------------------------------------------------
59//
60// Default constructor.
61//
62MPhotonData::MPhotonData(/*const char *name, const char *title*/)
63 : fPosX(0), fPosY(0), fCosU(0), fCosV(0), fTime(0), fWavelength(0),
64 /*fNumPhotons(1),*/ fProductionHeight(0), fPrimary(MMcEvtBasic::kUNDEFINED),
65 fTag(-1), fWeight(1)
66{
67 // fName = name ? name : "MPhotonData";
68 // fTitle = title ? title : "Corsika Event Data Information";
69}
70
71/*
72MPhotonData::MPhotonData(const MPhotonData &ph)
73: fPosX(ph.fPosX), fPosY(ph.fPosY), fCosU(ph.fCosU), fCosV(ph.fCosV),
74fTime(ph.fTime), fWavelength(ph.fWavelength), fNumPhotons(ph.fNumPhotons),
75fProductionHeight(ph.fProductionHeight), fPrimary(ph.fPrimary),
76fTag(ph.fTag), fWeight(ph.fWeight)
77{
78}
79*/
80
81// --------------------------------------------------------------------------
82//
83// Copy function. Copy all data members into obj.
84//
85void MPhotonData::Copy(TObject &obj) const
86{
87 MPhotonData &d = static_cast<MPhotonData&>(obj);
88
89// d.fNumPhotons = fNumPhotons;
90 d.fPosX = fPosX;
91 d.fPosY = fPosY;
92 d.fCosU = fCosU;
93 d.fCosV = fCosV;
94 d.fWavelength = fWavelength;
95 d.fPrimary = fPrimary;
96 d.fTime = fTime;
97 d.fTag = fTag;
98 d.fWeight = fWeight;
99 d.fProductionHeight = fProductionHeight;
100
101 TObject::Copy(obj);
102}
103
104// --------------------------------------------------------------------------
105//
106// Return the square cosine of the Theta-angle == 1-CosU^2-CosV^2
107//
108Double_t MPhotonData::GetCosW2() const
109{
110 return 1 - GetSinW2();
111}
112
113// --------------------------------------------------------------------------
114//
115// Return the square sine of the Theta-angle == CosU^2+CosV^2
116//
117Double_t MPhotonData::GetSinW2() const
118{
119 return fCosU*fCosU + fCosV*fCosV;
120}
121
122// --------------------------------------------------------------------------
123//
124// return the cosine of the Theta-angle == sqrt(1-CosU^2-CosV^2)
125//
126Double_t MPhotonData::GetCosW() const
127{
128 return TMath::Sqrt(GetCosW2());
129}
130
131// --------------------------------------------------------------------------
132//
133// return the sine of the Theta-angle == sqrt(CosU^2+CosV^2)
134//
135Double_t MPhotonData::GetSinW() const
136{
137 return TMath::Sqrt(GetSinW2());
138}
139
140// --------------------------------------------------------------------------
141//
142// Return the theta angle in radians
143//
144Double_t MPhotonData::GetTheta() const
145{
146 return TMath::ASin(GetSinW());
147}
148
149// --------------------------------------------------------------------------
150//
151// Return a TQuaternion with the first three components x, y, and z
152// and the fourth component the time.
153//
154TQuaternion MPhotonData::GetPosQ() const
155{
156 return TQuaternion(GetPos3(), fTime);
157}
158
159// --------------------------------------------------------------------------
160//
161// return a TQuaternion with the first three components the direction
162// moving in space (GetDir3()) and the fourth component is the
163// one devided by the speed of light (converted to cm/ns)
164//
165// FIXME: v in air!
166//
167TQuaternion MPhotonData::GetDirQ() const
168{
169 return TQuaternion(GetDir3(), 1./(TMath::C()*100/1e9));
170}
171
172// --------------------------------------------------------------------------
173//
174// Set the wavelength to a random lambda^-2 distributed value
175// between wmin and wmax.
176//
177void MPhotonData::SimWavelength(Float_t wmin, Float_t wmax)
178{
179 const Double_t w = gRandom->Uniform(wmin, wmax);
180
181 fWavelength = TMath::Nint(wmin*wmax / w);
182}
183
184
185// --------------------------------------------------------------------------
186//
187// Set the data member according to the 8 floats read from a reflector-file.
188// This function MUST reset all data-members, no matter whether these are
189// contained in the input stream.
190//
191Int_t MPhotonData::FillRfl(Float_t f[8])
192{
193 // Check coordinate system!!!!
194 fWavelength = TMath::Nint(f[0]);
195 fPosX = f[1]; // [cm]
196 fPosY = f[2]; // [cm]
197 fCosU = f[3]; // cos to x
198 fCosV = f[4]; // cos to y
199 fTime = f[5]; // [ns]
200 fProductionHeight = f[6];
201
202 // f[7]: Camera inclination angle
203
204 fPrimary = MMcEvtBasic::kUNDEFINED;
205// fNumPhotons = 1;
206 fTag = -1;
207 fWeight = 1;
208
209 return kTRUE;
210}
211
212// --------------------------------------------------------------------------
213//
214// Set the data member according to the 7 floats read from a corsika-file.
215// This function MUST reset all data-members, no matter whether these are
216// contained in the input stream.
217//
218// Currently we exchange x and y and set y=-y to convert Corsikas coordinate
219// system intpo our own.
220//
221Int_t MPhotonData::FillCorsika(Float_t f[7], Int_t i)
222{
223 const UInt_t n = TMath::Nint(f[0]);
224
225 if (n==0)
226 // FIXME: Do we need to decode the rest anyway?
227 return kCONTINUE;
228
229 // Check reuse
230 if (i >=0)
231 {
232 const Int_t reuse = (n/1000)%100; // Force this to be 1!
233 if (reuse!=i)
234 return kCONTINUE;
235 }
236
237 // This seems to be special to mmcs
238 fWavelength = n%1000;
239 fPrimary = MMcEvtBasic::ParticleId_t(n/100000);
240
241 // x=north, y=west
242 //fPosX = f[1]; // [cm]
243 //fPosY = f[2]; // [cm]
244 //fCosU = f[3]; // cos to x
245 //fCosV = f[4]; // cos to y
246 // x=west, y=south
247 fPosX = f[2]; // [cm]
248 fPosY = -f[1]; // [cm]
249
250 fCosU = f[4]; // cos to x
251 fCosV = -f[3]; // cos to y
252
253 fTime = f[5]; // [ns]
254
255 fProductionHeight = f[6]; // [cm]
256
257 // Now reset all data members which are not in the stream
258 fTag = -1;
259 fWeight = 1;
260
261 return kTRUE;
262}
263
264// --------------------------------------------------------------------------
265//
266// Set the data member according to the 8 floats read from a eventio-file.
267// This function MUST reset all data-members, no matter whether these are
268// contained in the input stream.
269//
270// Currently we exchange x and y and set y=-y to convert Corsikas coordinate
271// system into our own.
272//
273Int_t MPhotonData::FillEventIO(Float_t f[8])
274{
275 // photons in this bunch
276 const UInt_t n = TMath::Nint(f[6]);
277 if (n==0)
278 return 0;
279
280 f[6] = n-1;
281
282 fPosX = f[1]; // xpos relative to telescope [cm]
283 fPosY = -f[0]; // ypos relative to telescope [cm]
284 fCosU = f[3]; // cos to x
285 fCosV = -f[2]; // cos to y
286 fTime = f[4]; // a relative arival time [ns]
287 fProductionHeight = f[5]; // altitude of emission [cm]
288 fWavelength = TMath::Nint(f[7]); // so far always zeor = unspec. [nm]
289
290 // Now reset all data members which are not in the stream
291 fPrimary = MMcEvtBasic::kUNDEFINED;
292 fTag = -1;
293 fWeight = 1;
294
295 return n;
296}
297
298/*
299// --------------------------------------------------------------------------
300//
301// Read seven floats from the stream and call FillCorsika for them.
302//
303Int_t MPhotonData::ReadCorsikaEvt(istream &fin)
304{
305 Float_t f[7];
306 fin.read((char*)&f, 7*4);
307
308 const Int_t rc = FillCorsika(f);
309
310 return rc==kTRUE ? !fin.eof() : rc;
311}
312
313// --------------------------------------------------------------------------
314//
315// Read eight floats from the stream and call FillRfl for them.
316//
317Int_t MPhotonData::ReadRflEvt(istream &fin)
318{
319 Float_t f[8];
320 fin.read((char*)&f, 8*4);
321
322 const Int_t rc = FillRfl(f);
323
324 return rc==kTRUE ? !fin.eof() : rc;
325}
326*/
327
328// --------------------------------------------------------------------------
329//
330// Print contents. The tag and Weight are only printed if they are different
331// from the default.
332//
333void MPhotonData::Print(Option_t *) const
334{
335 gLog << inf << endl;
336// gLog << "Num Photons: " << fNumPhotons << " from " << MMcEvtBasic::GetParticleName(fPrimary) << endl;
337 gLog << "Origin: " << MMcEvtBasic::GetParticleName(fPrimary) << endl;
338 gLog << "Wavelength: " << fWavelength << "nm" << endl;
339 gLog << "Pos X/Y Cos U/V: " << fPosX << "/" << fPosY << " " << fCosU << "/" << fCosV << endl;
340 gLog << "Time/Prod.Height: " << fTime << "ns/" << fProductionHeight << "cm" << endl;
341 if (fTag>=0)
342 gLog << "Tag: " << fTag << endl;
343 if (fWeight!=1)
344 gLog << "Weight: " << fWeight << endl;
345}
Note: See TracBrowser for help on using the repository browser.