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

Last change on this file since 18677 was 18677, checked in by tbretz, 3 years ago
Fixed a numerical problem which could happen when the values were store with lower precision
File size: 11.0 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(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(Float_t f[7], Int_t i)
223{
224    const UInt_t n = TMath::Nint(f[0]);
225
226    if (n==0)
227        // FIXME: Do we need to decode the rest anyway?
228        return kCONTINUE;
229
230    // Check reuse
231    if (i >=0)
232      {
233       const Int_t reuse = (n/1000)%100; // Force this to be 1!
234       if (reuse!=i)
235           return kCONTINUE;
236      }
237
238    // This seems to be special to mmcs
239    fWavelength = n%1000;
240    fPrimary    = MMcEvtBasic::ParticleId_t(n/100000);
241
242    // x=north, y=west
243    //fPosX = f[1];  // [cm]
244    //fPosY = f[2];  // [cm]
245    //fCosU = f[3];  // cos to x
246    //fCosV = f[4];  // cos to y
247    // x=west, y=south
248    fPosX =  f[2];  // [cm]
249    fPosY = -f[1];  // [cm]
250
251    fCosU =  f[4];  // cos to x
252    fCosV = -f[3];  // cos to y
253
254    fTime =  f[5];  // [ns]
255
256    fProductionHeight = f[6]; // [cm]
257
258    // Now reset all data members which are not in the stream
259    fTag    = -1;
260    fWeight =  1;
261
262    return kTRUE;
263}
264
265// --------------------------------------------------------------------------
266//
267// Set the data member according to the 8 shorts read from a eventio-file.
268// This function MUST reset all data-members, no matter whether these are
269// contained in the input stream.
270//
271// Currently we exchange x and y and set y=-y to convert Corsikas coordinate
272// system into our own.
273//
274Int_t MPhotonData::FillEventIO(Short_t f[8])
275{
276    // From 5.5 compact_bunch:
277    // https://www.mpi-hd.mpg.de/hfm/~bernlohr/iact-atmo/iact_refman.pdf
278
279    // photons in this bunch f[6]/100.
280
281    fPosY             = -f[0]/10.;              // ypos relative to telescope [cm]
282    fPosX             =  f[1]/10.;              // xpos relative to telescope [cm]
283    fCosV             = -f[2]/30000.;           // cos to y
284    fCosU             =  f[3]/30000.;           // cos to x
285
286    fTime             =  f[4]/10.;              // a relative arival time [ns]
287    fProductionHeight =  pow(10, f[5]/1000.);   // altitude of emission a.s.l. [cm]
288    fWavelength       =  f[7];                  // wavelength [nm]: 0 undetermined, <0 already in p.e.
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 1;
296}
297
298// --------------------------------------------------------------------------
299//
300// Set the data member according to the 8 floats read from a eventio-file.
301// This function MUST reset all data-members, no matter whether these are
302// contained in the input stream.
303//
304// Currently we exchange x and y and set y=-y to convert Corsikas coordinate
305// system into our own.
306//
307Int_t MPhotonData::FillEventIO(Float_t f[8])
308{
309    // photons in this bunch
310    const UInt_t n = TMath::Nint(f[6]);
311    if (n==0)
312        return 0;
313
314    fPosX             =  f[1];              // xpos relative to telescope [cm]
315    fPosY             = -f[0];              // ypos relative to telescope [cm]
316    fCosU             =  f[3];              // cos to x
317    fCosV             = -f[2];              // cos to y
318    //fTime             =  f[4];              // a relative arival time [ns]
319    //fProductionHeight =  f[5];              // altitude of emission [cm]
320    fWavelength       =  0;                 // so far always zeor = unspec. [nm]
321
322    // Now reset all data members which are not in the stream
323    fPrimary    = MMcEvtBasic::kUNDEFINED;
324    fTag        = -1;
325    fWeight     =  1;
326
327    return n-1;
328}
329
330/*
331// --------------------------------------------------------------------------
332//
333// Read seven floats from the stream and call FillCorsika for them.
334//
335Int_t MPhotonData::ReadCorsikaEvt(istream &fin)
336{
337    Float_t f[7];
338    fin.read((char*)&f, 7*4);
339
340    const Int_t rc = FillCorsika(f);
341
342    return rc==kTRUE ? !fin.eof() : rc;
343}
344
345// --------------------------------------------------------------------------
346//
347// Read eight floats from the stream and call FillRfl for them.
348//
349Int_t MPhotonData::ReadRflEvt(istream &fin)
350{
351    Float_t f[8];
352    fin.read((char*)&f, 8*4);
353
354    const Int_t rc = FillRfl(f);
355
356    return rc==kTRUE ? !fin.eof() : rc;
357}
358*/
359
360// --------------------------------------------------------------------------
361//
362// Print contents. The tag and Weight are only printed if they are different
363// from the default.
364//
365void MPhotonData::Print(Option_t *) const
366{
367    gLog << inf << endl;
368//    gLog << "Num Photons:      " << fNumPhotons << " from " << MMcEvtBasic::GetParticleName(fPrimary) << endl;
369    gLog << "Origin:           " << MMcEvtBasic::GetParticleName(fPrimary) << endl;
370    gLog << "Wavelength:       " << fWavelength << "nm" << endl;
371    gLog << "Pos X/Y  Cos U/V: " << fPosX << "/" << fPosY << "   " << fCosU << "/" << fCosV << endl;
372    gLog << "Time/Prod.Height: " << fTime << "ns/" << fProductionHeight/100 << "m" << endl;
373    if (fTag>=0)
374        gLog << "Tag:              " << fTag << endl;
375    if (fWeight!=1)
376        gLog << "Weight:           " << fWeight << endl;
377}
Note: See TracBrowser for help on using the repository browser.