source: trunk/Mars/msim/MPhotonData.cc

Last change on this file was 19551, checked in by tbretz, 6 weeks ago
When the bunch size is 1 then the photons still carry their original probability (that is my interpretation) which is stored as 'bunch size'. Therefore, all bunch sizes <= 1 are ok.
File size: 16.1 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    // My understanding is that each photon internally gets a weight (e.g. lambda^-2)
263    // and according to this weight a dice is thrown. The weights are still written
264    // to the output file but only for the surviving photons
265    if (f[0]>1)
266        cout << "MPhotonData::FillCorsika: WARNING - Bunch size > 1 (" << f[0] << ")" <<endl;
267
268#ifdef __MMCS__
269    // Check reuse
270    if (i >=0)
271    {
272        const Int_t reuse = (n/1000)%100; // Force this to be 1!
273        if (reuse!=i)
274        {
275            cout << "REUSE " << reuse << " " << i << " " << n << endl;
276            return kCONTINUE;
277        }
278    }
279
280    // This seems to be special to mmcs
281    fWavelength = n%1000;
282    fPrimary    = MMcEvtBasic::ParticleId_t(n/100000);
283#else
284    fWavelength = 0;
285    fPrimary    = MMcEvtBasic::kUNDEFINED;
286#endif
287
288    // x=north, y=west
289    //fPosX = f[1];  // [cm]
290    //fPosY = f[2];  // [cm]
291    //fCosU = f[3];  // cos to x
292    //fCosV = f[4];  // cos to y
293    // x=west, y=south
294    fPosX =  f[2];  // [cm]
295    fPosY = -f[1];  // [cm]
296
297    fCosU =  f[4];  // cos to x
298    fCosV = -f[3];  // cos to y
299
300    fTime =  f[5];  // [ns]
301
302    fProductionHeight = f[6]; // [cm]
303
304    // Now reset all data members which are not in the stream
305    fTag    = -1;
306    fWeight =  1;
307
308    return kTRUE;
309}
310
311Int_t MPhotonData::FillCorsikaThin(const Float_t f[8], Int_t i)
312{
313    // DATAB2(LHCER(IBUF)+1,IBUF) = PHOTCM*WTCER/MAX(1.D-10,PROBTH)
314
315    // For the THIN option the photon bunch size is multiplied with the
316    // thinning weight. See Sect. 4.89 page 99.
317
318    /*
319     #if __THIN__
320     C      (CONTAINING UP TO 39 BUNCHES, 8 WORDS EACH)
321     C
322     C    8*(N-1)+1  NUMBER OF PHOTONS IN BUNCH
323     C                  (FOR STANDARD PARTICLE OUTPUT FILE:
324     C                   99.E5 + NINT(NUMBER OF PHOTONS IN BUNCH)*10 + 1
325     C    8*(N-1)+2  X- COORDINATE IN CM
326     C    8*(N-1)+3  Y- COORDINATE IN CM
327     C    8*(N-1)+4  DIRECTION COSINUS TO X AXIS
328     C    8*(N-1)+5  DIRECTION COSINUS TO Y AXIS
329     C    8*(N-1)+6  T TIME SINCE FIRST INTERACTION (OR ENTRANCE INTO
330     C               ATMOSPHERE) IN NSEC
331     C    8*(N-1)+7  PRODUCTION HEIGHT OF BUNCH IN CM
332     C    8*(N-1)+8  WEIGHT OF BUNCH
333     #else
334     */
335
336    // From the Corsika manual:
337    //
338    // f[0] : n Number of Cherenkov photons in bunch
339    //          (For THIN option multiplied with thinning weight)
340    // f[1] : x [cm]
341    // f[2] : y [cm]
342    // f[3] : u direction cosine (to x axis)  [u = sin(theta)cos(phi)]
343    // f[4] : v direction cosine (to y axis)  [v = sin(theta)sin(phi)]
344    // f[5] : t [ns] time to first interaction or since start of atmosphere (see TSTART)
345    // f[6] : h [ch] bunch production height (except MCERFI=3)
346    // f[7] : w weight of bunch [only with THIN option]
347
348    // f[0] = MCERFI==1/2/3 && !THIN ? PHOTCM : PHOTCM*WTCER/MAX(1e-10, PROBTH)
349    // f[1] = XCER
350    // f[2] = YXCER
351    // f[3] = UEMIS
352    // f[4] = VEMIS
353    // f[5] = CARTIM
354    // f[6] = MCERFI<3 ? ZEMIS : CERDIST
355    // #if __THIN__
356    // f[7] = MCERFI==1 ? WTCER/MAX(1e-10, PROBTH) : (CEFFIC || CERWLEN ? WL*WLFLAG : WLFLAG);
357    // #endif
358    //
359    // WLFLAG = [CEFFIC=-1] [CERWLEN=1] [ELSE=0]
360    //
361    // WL = 1./(1./WAVLGU+(1./WAVLGL-1./WAVLGU)*RDM(IRDM))
362    //
363    // WL-MIN = 1./(1./WAVLGU+(1./WAVLGL-1./WAVLGU)*0) = WAVLGU
364    // WL-MAX = 1./(1./WAVLGU+(1./WAVLGL-1./WAVLGU)*1) = WAVLGL
365    //
366
367    //const UInt_t n /*fWeight*/ = TMath::Nint(f[0]);
368
369    // Empty bunch (this happend at the end of events when
370    // the remaining block is filled with zeroes)
371    if (f[0]==0)
372        return kCONTINUE;
373
374    // My understanding is that each photon internally gets a weight (e.g. lambda^-2)
375    // and according to this weight a dice is thrown. The weights are still written
376    // to the output file but only for the surviving photons
377    if (f[0]>1)
378        cout << "MPhotonData::FillCorsikaThin: WARNING - Bunch size > 1 (" << f[0] << ")" <<endl;
379
380    // x=north, y=west
381    //fPosX = f[1];  // [cm]
382    //fPosY = f[2];  // [cm]
383    //fCosU = f[3];  // cos to x
384    //fCosV = f[4];  // cos to y
385    // x=west, y=south
386    fPosX =  f[2];  // [cm]
387    fPosY = -f[1];  // [cm]
388
389    fCosU =  f[4];  // cos to x
390    fCosV = -f[3];  // cos to y
391
392    fTime =  f[5];  // [ns]
393
394    fProductionHeight = f[6]; // [cm]
395    // f[7]<0: Photoelectron bunches of specific wavelength (if __CEFFIC__)
396    // f[7]>0: __CERWLEN__ (but __CEFFIC takes priority)
397
398    fWavelength = TMath::Abs(f[7]);
399
400    // Now reset all data members which are not in the stream
401    fPrimary = MMcEvtBasic::kUNDEFINED;
402    fTag     = -1;
403    fWeight  =  1;
404
405    return kTRUE;
406}
407
408// --------------------------------------------------------------------------
409//
410// Set the data member according to the 8 shorts read from a eventio-file.
411// This function MUST reset all data-members, no matter whether these are
412// contained in the input stream.
413//
414// Currently we exchange x and y and set y=-y to convert Corsikas coordinate
415// system into our own.
416//
417Int_t MPhotonData::FillEventIO(const Short_t f[8])
418{
419    // From 5.5 compact_bunch:
420    // https://www.mpi-hd.mpg.de/hfm/~bernlohr/iact-atmo/iact_refman.pdf
421
422    // photons in this bunch f[6]/100.
423
424    fPosY             = -f[0]/10.;              // ypos relative to telescope [cm]
425    fPosX             =  f[1]/10.;              // xpos relative to telescope [cm]
426    fCosV             = -f[2]/30000.;           // cos to y
427    fCosU             =  f[3]/30000.;           // cos to x
428
429    fTime             =  f[4]/10.;              // a relative arival time [ns]
430    fProductionHeight =  pow(10, f[5]/1000.);   // altitude of emission a.s.l. [cm]
431    fWavelength       =  TMath::Abs(f[7]);      // wavelength [nm]: 0 undetermined, <0 already in p.e.
432
433    // Now reset all data members which are not in the stream
434    fPrimary    = MMcEvtBasic::kUNDEFINED;
435    fTag        = -1;
436    fWeight     =  1;
437
438    return 1;
439}
440
441// --------------------------------------------------------------------------
442//
443// Set the data member according to the 8 floats read from a eventio-file.
444// This function MUST reset all data-members, no matter whether these are
445// contained in the input stream.
446//
447// Currently we exchange x and y and set y=-y to convert Corsikas coordinate
448// system into our own.
449//
450Int_t MPhotonData::FillEventIO(const Float_t f[8])
451{
452    // photons in this bunch
453    const UInt_t n = TMath::Nint(f[6]);
454    if (n==0)
455        return 0;
456
457    fPosX             =  f[1];              // xpos relative to telescope [cm]
458    fPosY             = -f[0];              // ypos relative to telescope [cm]
459    fCosU             =  f[3];              // cos to x
460    fCosV             = -f[2];              // cos to y
461    //fTime             =  f[4];              // a relative arival time [ns]
462    //fProductionHeight =  f[5];              // altitude of emission [cm]
463    fWavelength       =  0;                 // so far always zeor = unspec. [nm]
464
465    // Now reset all data members which are not in the stream
466    fPrimary    = MMcEvtBasic::kUNDEFINED;
467    fTag        = -1;
468    fWeight     =  1;
469
470    return n-1;
471}
472
473/*
474// --------------------------------------------------------------------------
475//
476// Read seven floats from the stream and call FillCorsika for them.
477//
478Int_t MPhotonData::ReadCorsikaEvt(istream &fin)
479{
480    Float_t f[7];
481    fin.read((char*)&f, 7*4);
482
483    const Int_t rc = FillCorsika(f);
484
485    return rc==kTRUE ? !fin.eof() : rc;
486}
487
488// --------------------------------------------------------------------------
489//
490// Read eight floats from the stream and call FillRfl for them.
491//
492Int_t MPhotonData::ReadRflEvt(istream &fin)
493{
494    Float_t f[8];
495    fin.read((char*)&f, 8*4);
496
497    const Int_t rc = FillRfl(f);
498
499    return rc==kTRUE ? !fin.eof() : rc;
500}
501*/
502
503// --------------------------------------------------------------------------
504//
505// Print contents. The tag and Weight are only printed if they are different
506// from the default.
507//
508void MPhotonData::Print(Option_t *) const
509{
510    gLog << inf << endl;
511    //    gLog << "Num Photons:      " << fNumPhotons << " from " << MMcEvtBasic::GetParticleName(fPrimary) << endl;
512    if (fPrimary!=MMcEvtBasic::kUNDEFINED)
513        gLog << "Origin:           " << MMcEvtBasic::GetParticleName(fPrimary) << endl;
514    gLog << "Wavelength:       " << fWavelength << "nm" << endl;
515    gLog << "Pos X/Y  Cos U/V: " << fPosX << "/" << fPosY << "   " << fCosU << "/" << fCosV << endl;
516    gLog << "Time/Prod.Height: " << fTime << "ns/" << fProductionHeight/100 << "m" << endl;
517    if (fTag>=0)
518        gLog << "Tag:              " << fTag << endl;
519    if (fWeight!=1)
520        gLog << "Weight:           " << fWeight << endl;
521}
Note: See TracBrowser for help on using the repository browser.