source: trunk/Mars/msim/MPhotonData.cc

Last change on this file was 19366, checked in by tbretz, 6 months ago
These two lines are obsolete.
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.