source: JEM-EUSO/esaf_cc_at_lal/packages/simulation/generators/showers/src/UnisimFileShowerGenerator.cc @ 114

Last change on this file since 114 was 114, checked in by moretto, 11 years ago

actual version of ESAF at CCin2p3
File size: 10.4 KB
Line 
1// ESAF : Euso Simulation and Analysis Framework
2// $Id: UnisimFileShowerGenerator.cc 2729 2006-06-08 14:54:45Z naumov $
3// sergio bottai created Jan, 22 2004
4
5#include "UnisimFileShowerGenerator.hh"
6#include "MCTruth.hh"
7#include "ShowerTrack.hh"
8#include <zlib.h>
9#include "EarthVector.hh"
10#include "Atmosphere.hh"
11#include <TMath.h>
12#include <TSystem.h>
13
14using namespace sou;
15using TMath::RadToDeg;
16
17ClassImp(UnisimFileShowerGenerator)
18
19//______________________________________________________________________________
20UnisimFileShowerGenerator::UnisimFileShowerGenerator() : 
21    EventGenerator("Unisimfile"), fTrack(0),fFile(NULL), fTruth(0) {
22    //
23    // Construtor
24    //
25   
26    fFileName = "";
27    fFirstEvent = 0;
28    fCurrentEvent = 0;
29}
30
31//______________________________________________________________________________
32UnisimFileShowerGenerator::~UnisimFileShowerGenerator() {
33    //
34    // Destructor
35    //
36
37    if ( fTrack ) {
38        delete fTrack;
39        fTrack=NULL; 
40    }
41    if ( fTruth ) {
42        delete fTruth; 
43        fTruth=NULL;
44    }
45   
46   
47}
48
49//______________________________________________________________________________
50void UnisimFileShowerGenerator::ReplaceInputFile(const char* fn) {
51    //
52    //
53    //
54   
55    if ( fFile )
56        throw runtime_error("UnisimFileShowerGenerator: Cannot change name when file is already open!");
57    if (!fn) {
58        Msg(EsafMsg::Error) << "UnisimFileShowerGenerator: NULL file name ignored" << MsgDispatch;
59        return;
60    }
61    Conf()->ReplaceStr("UnisimFileShowerGenerator.FileName",fn);
62    Msg(EsafMsg::Info) << "Input file name changed to " << fn << MsgDispatch;
63}
64
65
66
67//______________________________________________________________________________
68void UnisimFileShowerGenerator::Open() {
69    //
70    // Open file
71    //
72
73    if ( Conf()->IsDefined( "UnisimFileShowerGenerator.FileName" ) ) {
74        fFileName = Conf()->GetStrPath("UnisimFileShowerGenerator.FileName");
75        fFirstEvent = (Int_t)Conf()->GetNum("UnisimFileShowerGenerator.FirstEvent") ;
76        if (fFirstEvent < 0 ) {
77            Msg(EsafMsg::Warning) << "First event must be >= 0. Now forced fFirstEvent=0" << MsgDispatch;
78            fFirstEvent = 0;
79        }
80    }
81    else
82        Msg(EsafMsg::Panic) << "Unknown input file name in UnisimFileShowerGenerator" << MsgDispatch;
83
84    fFile = (FILE*)gzopen( fFileName.c_str(), "rb" );
85
86    if( fFile == NULL ) 
87        Msg(EsafMsg::Panic) << "Unknown input file name in UnisimFileShowerGenerator" << MsgDispatch;
88    else
89        Msg(EsafMsg::Info) << "File: " << fFileName << " successfully opened " << MsgDispatch;
90
91}
92
93//______________________________________________________________________________
94void UnisimFileShowerGenerator::Close() {
95    //
96    // Close file
97    //
98
99    if( fFile  ) gzclose( (void*)fFile ) ;
100   
101}
102
103
104//______________________________________________________________________________
105PhysicsData* UnisimFileShowerGenerator::Get() {
106    //
107    // Get a new event
108    //
109
110
111    // build and return a track
112
113    if ( !fFile ) Open(); 
114    Int_t event=0;
115
116    while (1) {
117
118        if ( fTrack ) {
119            delete fTrack;
120            fTrack=NULL;
121        }
122        if ( fTruth ) {
123            delete fTruth;
124            fTruth=NULL;
125        }
126
127        if (fCurrentEvent >= fFirstEvent) {
128            if ( LoadTrack(event) ){
129
130                MsgForm(EsafMsg::Info,"UnisimFile event %d, header id %d", fCurrentEvent, (Int_t)fHeader["EVTNUMBR"]);
131                MsgForm(EsafMsg::Info,
132                        "ShowerTrack produced : Energy = %.2E MeV, Theta= %.2f deg, Phi= %.2f deg",
133                        fTrack->GetEnergy(), fTrack->GetTheta()*RadToDeg(),fTrack->GetPhi()*RadToDeg());
134                MsgForm(EsafMsg::Info,
135                        "         First interaction X1 = %.3f at : (%.2f,%.2f,%.2f) km gives a shower with %d steps",
136                        fTrack->fX1*cm2/g,fTrack->GetInitPos().X()/km,fTrack->GetInitPos().Y()/km,
137                        fTrack->GetInitPos().Z()/km, fTrack->GetNumStep());
138            } else {
139                Msg(EsafMsg::Warning) << "UnisimFileShowerGenerator: End of file "
140                    << fFileName << " reached? " << MsgDispatch;
141                fTrack = (ShowerTrack*)NULL;
142            }
143            fCurrentEvent++;
144            break;
145        } else {
146            if ( !LoadTrack(event) ){
147                Msg(EsafMsg::Warning) << "UnisimFileShowerGenerator: End of file "
148                    << fFileName << " reached? " << MsgDispatch;
149                fTrack = (ShowerTrack*)NULL;
150                fCurrentEvent++;
151                break;
152            }
153        }
154
155        fCurrentEvent++;
156
157    }
158    fTrack->CheckTrack();
159    return fTrack;
160
161}
162
163//______________________________________________________________________________
164MCTruth* UnisimFileShowerGenerator::GetTruth() {
165    //
166    // get pointer to MonteCarlo Truth
167    //
168
169    if ( fTrack ) {
170        return fTruth;
171    }
172        else { Msg(EsafMsg::Warning) << "Event to be loaded with LoadTrack before GetTruth" << MsgDispatch;
173        return (MCTruth*)0;
174    }
175   
176
177}
178
179//______________________________________________________________________________
180bool UnisimFileShowerGenerator::LoadTrack( Int_t event ) {
181    //
182    // Load one event from file
183    //
184   
185    if ( !LoadHeader() ) return false;
186
187    char s[5010]; 
188
189    fTruth= new MCTruth();
190     
191    LoadTruth();
192
193    fTrack= new ShowerTrack();
194    fTrack->fEthrEl=0.;
195   
196    fTrack->fDirVers[0]=fHeader["ELDIRCOS"];
197    fTrack->fDirVers[1]=fHeader["EMDIRCOS"];
198    fTrack->fDirVers[2]=fHeader["ENDIRCOS"];
199   
200    fTrack->fInitPos[0]=fHeader["EX0COMES"]*m;
201    fTrack->fInitPos[1]=fHeader["EY0COMES"]*m;
202    fTrack->fInitPos[2]=fHeader["EZ0COMES"]*m;
203   
204    if(fHeader["HITGROUN"]==0) fTrack->fHitGround=false ;
205    else fTrack->fHitGround=true ;
206   
207    fTrack->fEnergy    = fHeader["EPRMENRG"]*GeV;
208    fTrack->fTheta     = fHeader["EPOLANGD"]*deg; 
209    fTrack->fPhi       = fHeader["EAZIANGD"]*deg;
210    fTrack->fX1        = fHeader["EPRMDPTH"]*g/cm2;
211    fTrack->fEarthImpact.SetXYZ(fHeader["EXDMPPOS"]*m,fHeader["EYDMPPOS"]*m,fHeader["EZDMPPOS"]*m);
212   
213    while( 1 ) {
214        char* ans = gzgets((void*)fFile,s,5000);
215
216        if ( ans == Z_NULL ) {
217            Msg(EsafMsg::Error) << "Error reading gzgets in UnisimFileShowerGenerator::Load()\n" << MsgDispatch;
218            delete fTrack;
219            fTrack=NULL;
220            return false;
221        }
222        if (strncmp(s," BEVTBEND EEND",14) == 0) return true;
223        if (strncmp(s," STEPNUMB",9) == 0) {
224
225            char key;
226            int n;
227            if ( sscanf(s,"%s %d %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf",
228                        &key,&n,&fStep.fXi,&fStep.fXf,&fStep.fXYZi[0],&fStep.fXYZi[1],&fStep.fXYZi[2],
229                        &fStep.fXYZf[0],&fStep.fXYZf[1],&fStep.fXYZf[2],
230                        &fStep.fTimei,&fStep.fTimef,&fStep.fAgei,&fStep.fAgef,&fStep.fNelectrons,
231                        &fStep.fNcharged,&fStep.fEloss,&fStep.fNcherenkov) != 18 ) {
232                Msg(EsafMsg::Error) << "Data format error in UnisimFileShowerGenerator::LoadHeader()" << MsgDispatch;
233                Msg(EsafMsg::Error) << s << MsgDispatch;
234                return false;
235            }
236
237            fStep.fXYZi[0]=fStep.fXYZi[0]*m;
238            fStep.fXYZi[1]=fStep.fXYZi[1]*m;
239            fStep.fXYZi[2]=fStep.fXYZi[2]*m;
240
241            fStep.fXYZf[0]=fStep.fXYZf[0]*m;
242            fStep.fXYZf[1]=fStep.fXYZf[1]*m;
243            fStep.fXYZf[2]=fStep.fXYZf[2]*m;
244
245            fStep.fTimei=fStep.fTimei*microsecond;     
246            fStep.fTimef=fStep.fTimef*microsecond;     
247
248            fStep.fXi=fStep.fXi*gram/cm2;
249            fStep.fXf=fStep.fXf*gram/cm2;
250
251            fTrack->Add(fStep);
252
253        }
254    }
255    return true;
256}
257
258//______________________________________________________________________________
259bool UnisimFileShowerGenerator::LoadHeader() {
260    //
261    // Load header and store it into a map (key,value)
262    //
263
264
265    fHeader.clear();
266
267    char s[5010];
268    while(1) {
269        char* ans = gzgets((void*)fFile,s,5000);
270
271        if ( ans == Z_NULL ) {
272            Msg(EsafMsg::Error) << "Error reading gzgets in UnisimFileShowerGenerator::LoadHeader() (2) " << MsgDispatch;
273            return false;
274        }
275        if (strncmp(s," BEVTBUFF EVTB",14) == 0) break;
276
277        char key[200];
278        double val=0.;
279        if ((strncmp(s," EVTHEADR EVTH",14) != 0) && (strncmp(s," RUNHEADR RUNH",14) != 0)) {
280            if ( sscanf(s,"%s %lf",key,&val) != 2 ) {
281                Msg(EsafMsg::Error) << "Data format error in UnisimFileShowerGenerator::LoadHeader()" << MsgDispatch;
282                Msg(EsafMsg::Error) << s << MsgDispatch;
283                return false;
284            }
285            fHeader[key]=val;
286        }
287    }
288    return true;
289}
290
291
292//______________________________________________________________________________
293void UnisimFileShowerGenerator::LoadTruth() {
294    //
295    //
296    //
297
298    fTruth->SetEnergy(fHeader["EPRMENRG"]*GeV); // in ESAF units (MeV)
299    fTruth->SetThetaPhi(fHeader["EPOLANGD"]*deg, fHeader["EAZIANGD"]*deg);
300    EVector v;
301    v[0] = fHeader["EX0COMES"]*m;
302    v[1] = fHeader["EY0COMES"]*m;
303    v[2] = fHeader["EZ0COMES"]*m;
304    fTruth->SetFirstInt(v,fHeader["EPRMDPTH"]*g/cm2);
305    v[0] = fHeader["EXDMPPOS"]*m;
306    v[1] = fHeader["EYDMPPOS"]*m;
307    v[2] = fHeader["EZDMPPOS"]*m;
308    // FIXME, Age on the Earth is not stored in files. Puting arbitrary number 2.
309    fTruth->SetEarthImpact(v,2);     
310   
311    // recalculate TOAimpact
312    EarthVector first_inter(fHeader["EX0COMES"]*m,fHeader["EY0COMES"]*m,fHeader["EZ0COMES"]*m);
313    EarthVector Omega(1);
314    EarthVector toaimpact(1);
315    Omega.SetMagThetaPhi(1.,fHeader["EPOLANGD"]*deg,fHeader["EAZIANGD"]*deg);
316    toaimpact = Atmosphere::Get()->ImpactAtTOA(first_inter,EarthVector(-Omega));
317    fTruth->SetTOAImpact(toaimpact);
318           
319    // In case no value for ShowerMax is present  Puting arbitrary number 0.
320    v[0] = 0;
321    v[1] = 0;
322    v[2] = 0;
323    map<string,double>::iterator p;
324    p = fHeader.find("EXMAXPOS") ;
325    if ( p != fHeader.end() ) v[0] = fHeader["EXMAXPOS"]*m;
326    p = fHeader.find("EYMAXPOS") ;
327    if ( p != fHeader.end() ) v[1] = fHeader["EYMAXPOS"]*m;
328    p = fHeader.find("EZMAXPOS") ;
329    if ( p != fHeader.end() ) v[2] = fHeader["EZMAXPOS"]*m;         
330
331
332    fTruth->SetShowerMax(v,fHeader["ESHOWMAX"]*g/cm2);
333   
334    fTruth-> SetHclouds(0.);
335    p = fHeader.find("ESTAHEIG") ;
336    if ( p != fHeader.end() ) fTruth-> SetHclouds(fHeader["ESTAHEIG"]*km); // in ESAF units
337
338}
Note: See TracBrowser for help on using the repository browser.