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

Last change on this file since 19339 was 19339, checked in by tbretz, 6 months ago
Added the 8-float format as FillEvtCorsikaThin, calculate absolute wavelength (right now, I don't know what to do with the information from the negative WL, print the origin only if it is defined, added some docu from the corsika docu, removed MMCS specific stuff by a pre-compiler directive -- do we still have MMCS data?
File size: 14.6 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), fMirrorTag(-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    const Double_t sinw2 = fCosU*fCosU + fCosV*fCosV;
120    return sinw2>1 ? 1 : sinw2;
121}
122
123// --------------------------------------------------------------------------
124//
125// return the cosine of the Theta-angle == sqrt(1-CosU^2-CosV^2)
126//
127Double_t MPhotonData::GetCosW() const
128{
129    return TMath::Sqrt(GetCosW2());
130}
131
132// --------------------------------------------------------------------------
133//
134// return the sine of the Theta-angle == sqrt(CosU^2+CosV^2)
135//
136Double_t MPhotonData::GetSinW() const
137{
138    return TMath::Sqrt(GetSinW2());
139}
140
141// --------------------------------------------------------------------------
142//
143// Return the theta angle in radians
144//
145Double_t MPhotonData::GetTheta() const
146{
147    return TMath::ASin(GetSinW());
148}
149
150// --------------------------------------------------------------------------
151//
152// Return a TQuaternion with the first three components x, y, and z
153// and the fourth component the time.
154//
155TQuaternion MPhotonData::GetPosQ() const
156{
157    return TQuaternion(GetPos3(), fTime);
158}
159
160// --------------------------------------------------------------------------
161//
162// return a TQuaternion with the first three components the direction
163// moving in space (GetDir3()) and the fourth component is the
164// one devided by the speed of light (converted to cm/ns)
165//
166// FIXME: v in air!
167//
168TQuaternion MPhotonData::GetDirQ() const
169{
170    return TQuaternion(GetDir3(), 1./(TMath::C()*100/1e9));
171}
172
173// --------------------------------------------------------------------------
174//
175// Set the wavelength to a random lambda^-2 distributed value
176// between wmin and wmax.
177//
178void MPhotonData::SimWavelength(Float_t wmin, Float_t wmax)
179{
180    const Double_t w = gRandom->Uniform(wmin, wmax);
181
182    fWavelength = TMath::Nint(wmin*wmax / w);
183}
184
185
186// --------------------------------------------------------------------------
187//
188// Set the data member according to the 8 floats read from a reflector-file.
189// This function MUST reset all data-members, no matter whether these are
190// contained in the input stream.
191//
192Int_t MPhotonData::FillRfl(const Float_t f[8])
193{
194    // Check coordinate system!!!!
195    fWavelength = TMath::Nint(f[0]);
196    fPosX = f[1];  // [cm]
197    fPosY = f[2];  // [cm]
198    fCosU = f[3];  // cos to x
199    fCosV = f[4];  // cos to y
200    fTime = f[5];  // [ns]
201    fProductionHeight = f[6];
202
203    // f[7]: Camera inclination angle
204
205    fPrimary    = MMcEvtBasic::kUNDEFINED;
206//    fNumPhotons =  1;
207    fTag        = -1;
208    fWeight     =  1;
209
210    return kTRUE;
211}
212
213// --------------------------------------------------------------------------
214//
215// Set the data member according to the 7 floats read from a corsika-file.
216// This function MUST reset all data-members, no matter whether these are
217// contained in the input stream.
218//
219// Currently we exchange x and y and set y=-y to convert Corsikas coordinate
220// system intpo our own.
221//
222Int_t MPhotonData::FillCorsika(const Float_t f[7], Int_t i)
223{
224    // From the Corsika manual:
225    //
226    // f[0] : n Number of Cherenkov photons in bunch
227    //          (For THIN option multiplied with thinning weight)
228    // f[1] : x [cm]
229    // f[2] : y [cm]
230    // f[3] : u direction cosine (to x axis)  [u = sin(theta)cos(phi)]
231    // f[4] : v direction cosine (to y axis)  [v = sin(theta)sin(phi)]
232    // f[5] : t [ns] time to first interaction or since start of atmosphere (see TSTART)
233    // f[6] : h [ch] bunch production height (except MCERFI=3)
234    // f[7] : w weight of bunch [only with THIN option]
235
236    // f[0] = MCERFI==1/2/3 && !THIN ? PHOTCM : PHOTCM*WTCER/MAX(1e-10, PROBTH)
237    // f[1] = XCER
238    // f[2] = YXCER
239    // f[3] = UEMIS
240    // f[4] = VEMIS
241    // f[5] = CARTIM
242    // f[6] = MCERFI<3 ? ZEMIS : CERDIST
243    // #if __THIN__
244    // f[7] = MCERFI==1 ? WTCER/MAX(1e-10, PROBTH) : (CEFFIC || CERWLEN ? WL*WLFLAG : WLFLAG);
245    // #endif
246    //
247    // WLFLAG = [CEFFIC=-1] [CERWLEN=1] [ELSE=0]
248    //
249    // WL = 1./(1./WAVLGU+(1./WAVLGL-1./WAVLGU)*RDM(IRDM))
250    //
251    // WL-MIN = 1./(1./WAVLGU+(1./WAVLGL-1./WAVLGU)*0) = WAVLGU
252    // WL-MAX = 1./(1./WAVLGU+(1./WAVLGL-1./WAVLGU)*1) = WAVLGL
253    //
254
255    const UInt_t n /*fWeight*/ = TMath::Nint(f[0]);
256
257    // Empty bunch (this happend at the end of events when
258    // the remaining block is filled with zeroes)
259    if (f[0]==0)
260        return kCONTINUE;
261
262    if (f[0]!=1)
263        cout << "WARNING - Bunch size !=1 (" << f[0] << ")" <<endl;
264
265#ifdef __MMCS__
266    // Check reuse
267    if (i >=0)
268    {
269        const Int_t reuse = (n/1000)%100; // Force this to be 1!
270        if (reuse!=i)
271        {
272            cout << "REUSE " << reuse << " " << i << " " << n << endl;
273            return kCONTINUE;
274        }
275    }
276
277    // This seems to be special to mmcs
278    fWavelength = n%1000;
279    fPrimary    = MMcEvtBasic::ParticleId_t(n/100000);
280#else
281    fWavelength = 0;
282    fPrimary    = MMcEvtBasic::kUNDEFINED;
283#endif
284
285    // x=north, y=west
286    //fPosX = f[1];  // [cm]
287    //fPosY = f[2];  // [cm]
288    //fCosU = f[3];  // cos to x
289    //fCosV = f[4];  // cos to y
290    // x=west, y=south
291    fPosX =  f[2];  // [cm]
292    fPosY = -f[1];  // [cm]
293
294    fCosU =  f[4];  // cos to x
295    fCosV = -f[3];  // cos to y
296
297    fTime =  f[5];  // [ns]
298
299    fProductionHeight = f[6]; // [cm]
300
301    // Now reset all data members which are not in the stream
302    fTag    = -1;
303    fWeight =  1;
304
305    return kTRUE;
306}
307
308Int_t MPhotonData::FillCorsikaThin(const Float_t f[8], Int_t i)
309{
310    // From the Corsika manual:
311    //
312    // f[0] : n Number of Cherenkov photons in bunch
313    //          (For THIN option multiplied with thinning weight)
314    // f[1] : x [cm]
315    // f[2] : y [cm]
316    // f[3] : u direction cosine (to x axis)  [u = sin(theta)cos(phi)]
317    // f[4] : v direction cosine (to y axis)  [v = sin(theta)sin(phi)]
318    // f[5] : t [ns] time to first interaction or since start of atmosphere (see TSTART)
319    // f[6] : h [ch] bunch production height (except MCERFI=3)
320    // f[7] : w weight of bunch [only with THIN option]
321
322    // f[0] = MCERFI==1/2/3 && !THIN ? PHOTCM : PHOTCM*WTCER/MAX(1e-10, PROBTH)
323    // f[1] = XCER
324    // f[2] = YXCER
325    // f[3] = UEMIS
326    // f[4] = VEMIS
327    // f[5] = CARTIM
328    // f[6] = MCERFI<3 ? ZEMIS : CERDIST
329    // #if __THIN__
330    // f[7] = MCERFI==1 ? WTCER/MAX(1e-10, PROBTH) : (CEFFIC || CERWLEN ? WL*WLFLAG : WLFLAG);
331    // #endif
332    //
333    // WLFLAG = [CEFFIC=-1] [CERWLEN=1] [ELSE=0]
334    //
335    // WL = 1./(1./WAVLGU+(1./WAVLGL-1./WAVLGU)*RDM(IRDM))
336    //
337    // WL-MIN = 1./(1./WAVLGU+(1./WAVLGL-1./WAVLGU)*0) = WAVLGU
338    // WL-MAX = 1./(1./WAVLGU+(1./WAVLGL-1./WAVLGU)*1) = WAVLGL
339    //
340
341    const UInt_t n /*fWeight*/ = TMath::Nint(f[0]);
342
343    // Empty bunch (this happend at the end of events when
344    // the remaining block is filled with zeroes)
345    if (f[0]==0)
346        return kCONTINUE;
347
348    if (f[0]!=1)
349        cout << "WARNING - Bunch size !=1 (" << f[0] << ")" <<endl;
350
351
352    // x=north, y=west
353    //fPosX = f[1];  // [cm]
354    //fPosY = f[2];  // [cm]
355    //fCosU = f[3];  // cos to x
356    //fCosV = f[4];  // cos to y
357    // x=west, y=south
358    fPosX =  f[2];  // [cm]
359    fPosY = -f[1];  // [cm]
360
361    fCosU =  f[4];  // cos to x
362    fCosV = -f[3];  // cos to y
363
364    fTime =  f[5];  // [ns]
365
366    fProductionHeight = f[6]; // [cm]
367
368    fWavelength = TMath::Abs(f[7]);
369
370    // Now reset all data members which are not in the stream
371    fPrimary = MMcEvtBasic::kUNDEFINED;
372    fTag     = -1;
373    fWeight  =  1;
374
375    return kTRUE;
376}
377
378// --------------------------------------------------------------------------
379//
380// Set the data member according to the 8 shorts read from a eventio-file.
381// This function MUST reset all data-members, no matter whether these are
382// contained in the input stream.
383//
384// Currently we exchange x and y and set y=-y to convert Corsikas coordinate
385// system into our own.
386//
387Int_t MPhotonData::FillEventIO(const Short_t f[8])
388{
389    // From 5.5 compact_bunch:
390    // https://www.mpi-hd.mpg.de/hfm/~bernlohr/iact-atmo/iact_refman.pdf
391
392    // photons in this bunch f[6]/100.
393
394    fPosY             = -f[0]/10.;              // ypos relative to telescope [cm]
395    fPosX             =  f[1]/10.;              // xpos relative to telescope [cm]
396    fCosV             = -f[2]/30000.;           // cos to y
397    fCosU             =  f[3]/30000.;           // cos to x
398
399    fTime             =  f[4]/10.;              // a relative arival time [ns]
400    fProductionHeight =  pow(10, f[5]/1000.);   // altitude of emission a.s.l. [cm]
401    fWavelength       =  TMath::Abs(f[7]);      // wavelength [nm]: 0 undetermined, <0 already in p.e.
402
403    // Now reset all data members which are not in the stream
404    fPrimary    = MMcEvtBasic::kUNDEFINED;
405    fTag        = -1;
406    fWeight     =  1;
407
408    return 1;
409}
410
411// --------------------------------------------------------------------------
412//
413// Set the data member according to the 8 floats read from a eventio-file.
414// This function MUST reset all data-members, no matter whether these are
415// contained in the input stream.
416//
417// Currently we exchange x and y and set y=-y to convert Corsikas coordinate
418// system into our own.
419//
420Int_t MPhotonData::FillEventIO(const Float_t f[8])
421{
422    // photons in this bunch
423    const UInt_t n = TMath::Nint(f[6]);
424    if (n==0)
425        return 0;
426
427    fPosX             =  f[1];              // xpos relative to telescope [cm]
428    fPosY             = -f[0];              // ypos relative to telescope [cm]
429    fCosU             =  f[3];              // cos to x
430    fCosV             = -f[2];              // cos to y
431    //fTime             =  f[4];              // a relative arival time [ns]
432    //fProductionHeight =  f[5];              // altitude of emission [cm]
433    fWavelength       =  0;                 // so far always zeor = unspec. [nm]
434
435    // Now reset all data members which are not in the stream
436    fPrimary    = MMcEvtBasic::kUNDEFINED;
437    fTag        = -1;
438    fWeight     =  1;
439
440    return n-1;
441}
442
443/*
444// --------------------------------------------------------------------------
445//
446// Read seven floats from the stream and call FillCorsika for them.
447//
448Int_t MPhotonData::ReadCorsikaEvt(istream &fin)
449{
450    Float_t f[7];
451    fin.read((char*)&f, 7*4);
452
453    const Int_t rc = FillCorsika(f);
454
455    return rc==kTRUE ? !fin.eof() : rc;
456}
457
458// --------------------------------------------------------------------------
459//
460// Read eight floats from the stream and call FillRfl for them.
461//
462Int_t MPhotonData::ReadRflEvt(istream &fin)
463{
464    Float_t f[8];
465    fin.read((char*)&f, 8*4);
466
467    const Int_t rc = FillRfl(f);
468
469    return rc==kTRUE ? !fin.eof() : rc;
470}
471*/
472
473// --------------------------------------------------------------------------
474//
475// Print contents. The tag and Weight are only printed if they are different
476// from the default.
477//
478void MPhotonData::Print(Option_t *) const
479{
480    gLog << inf << endl;
481    //    gLog << "Num Photons:      " << fNumPhotons << " from " << MMcEvtBasic::GetParticleName(fPrimary) << endl;
482    if (fPrimary!=MMcEvtBasic::kUNDEFINED)
483        gLog << "Origin:           " << MMcEvtBasic::GetParticleName(fPrimary) << endl;
484    gLog << "Wavelength:       " << fWavelength << "nm" << endl;
485    gLog << "Pos X/Y  Cos U/V: " << fPosX << "/" << fPosY << "   " << fCosU << "/" << fCosV << endl;
486    gLog << "Time/Prod.Height: " << fTime << "ns/" << fProductionHeight/100 << "m" << endl;
487    if (fTag>=0)
488        gLog << "Tag:              " << fTag << endl;
489    if (fWeight!=1)
490        gLog << "Weight:           " << fWeight << endl;
491}
Note: See TracBrowser for help on using the repository browser.