source: JEM-EUSO/esaf_cc_at_lal/packages/simulation/generators/showers/src/ConexFileShowerGenerator.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: 15.9 KB
Line 
1// $Id$
2// Author: kenjikry   2009/03/31
3
4/*****************************************************************************
5 * ESAF: Euso Simulation and Analysis Framework                              *
6 *                                                                           *
7 *  Id: ConexFileGenerator                                                   *
8 *  Package: <packagename>                                                   *
9 *  Coordinator: <coordinator>                                               *
10 *                                                                           *
11 *****************************************************************************/
12
13//_____________________________________________________________________________
14//
15// ConexFileGenerator
16//
17// <extensive class description>
18//
19//   Config file parameters
20//   ======================
21//
22//   <parameter name>: <parameter description>
23//   -Valid options: <available options>
24//
25
26#include "ConexFileShowerGenerator.hh"
27#include <TFile.h>
28#include <TTree.h>
29#include "MCTruth.hh"
30#include "ShowerTrack.hh"
31#include "Atmosphere.hh"
32#include <TMath.h>
33#include <TSystem.h>
34#include <TLorentzVector.h>
35#include <iostream>
36#include "EConst.hh"
37using namespace sou;
38using namespace std;
39using TMath::RadToDeg;
40
41ClassImp(ConexFileShowerGenerator)
42
43//_____________________________________________________________________________
44ConexFileShowerGenerator::ConexFileShowerGenerator(): EventGenerator("conex"), EsafMsgSource(){
45  //
46  // Constructor
47  //
48  fFirstEvent = 0;
49  fFileName = "";
50  fTruth = NULL;
51  fTrack = NULL; 
52  fConexHeader = NULL;
53  fConexShower = NULL;
54  fConexFirstInteraction = NULL;
55  pMom = NULL;
56}
57
58//_____________________________________________________________________________
59ConexFileShowerGenerator::~ConexFileShowerGenerator() {
60  //
61  // Destructor
62  //
63}
64
65
66
67//_____________________________________________________________________________
68void ConexFileShowerGenerator::OpenFile() {
69  //
70  // open the input file
71  //
72  if ( Conf()->IsDefined( "ConexFileShowerGenerator.FileName" ) ) {
73    fFileName = Conf()->GetStrPath("ConexFileShowerGenerator.FileName");
74    fFirstEvent = (Int_t)Conf()->GetNum("ConexFileShowerGenerator.FirstEvent") ;
75  }
76  else {
77    Msg(EsafMsg::Panic) << "Unknown input file name in ConexFileShowerGenerator" << MsgDispatch;
78  }
79 
80  fForceTheta = Conf()->GetBool("ConexFileShowerGenerator.fForceTheta");
81  fForcePhi = Conf()->GetBool("ConexFileShowerGenerator.fForcePhi");
82  fTheta = Conf()->GetNum("ConexFileShowerGenerator.fTheta");
83  fPhi =Conf()->GetNum("ConexFileShowerGenerator.fPhi");
84
85  fCoreX = Conf()->GetNum("ConexFileShowerGenerator.fCoreX");
86  fCoreY = Conf()->GetNum("ConexFileShowerGenerator.fCoreY");
87  fCoreH = Conf()->GetNum("ConexFileShowerGenerator.fCoreH");
88
89  fInputFile = TFile::Open(fFileName);
90  fConexHeader = (TTree*)fInputFile->Get("Header");
91  fConexShower = (TTree*)fInputFile->Get("Shower");
92  fCurrentEvent=fFirstEvent;
93  InitTrees();
94}
95
96//_____________________________________________________________________________
97void ConexFileShowerGenerator::CloseFile() {
98  //
99  // close the input file
100  //
101  fInputFile->Close();
102  fConexHeader = NULL;
103  fConexShower = NULL;
104  fConexFirstInteraction = NULL;
105}
106
107//_____________________________________________________________________________
108void ConexFileShowerGenerator::InitTrees() {
109  //
110  // initialize the trees in conex file
111  //
112        float OutputVersion;
113        fConexHeader->SetBranchAddress("OutputVersion",&OutputVersion);
114        fConexHeader->GetEntry(0);
115
116        if(OutputVersion<2.0){
117                Msg(EsafMsg::Info) << "Reading old-style CONEX file version 1.x" << MsgDispatch;
118
119          fConexFirstInteraction = (TTree*)fInputFile->Get("FirstInteraction");
120
121                fConexShower->SetBranchAddress("nX",&nX); // number of points in slant depth
122                fConexShower->SetBranchAddress("N",N);     // particles(X)
123                fConexShower->SetBranchAddress("dEdX",dEdX); // Energy deposit(X)
124                fConexShower->SetBranchAddress("Electrons",Electrons);
125                fConexHeader->SetBranchAddress("X",X);       // slant depth
126                fConexHeader->SetBranchAddress("H",H);       // Height
127                fConexHeader->SetBranchAddress("D",D);       // Distance
128                fConexHeader->SetBranchAddress("delX",&delX);// depth step
129
130                fConexFirstInteraction->SetBranchAddress("depth",depth);
131
132                fConexShower->SetBranchAddress("lgE",&fitpars[0]);
133                fConexShower->SetBranchAddress("zenith",&fitpars[1]);
134                fConexShower->SetBranchAddress("azimuth",&fitpars[10]);
135                fConexShower->SetBranchAddress("Xmax",&fitpars[9]);
136
137                fConexFirstInteraction->SetBranchAddress("height",height);
138                fConexFirstInteraction->SetBranchAddress("idInt",IdInt1);
139                fConexFirstInteraction->SetBranchAddress("nPart",&nPart);
140                fConexFirstInteraction->SetBranchAddress("nInt",&nInt);
141        } else if(OutputVersion>2.0 && OutputVersion<3.0){
142                Msg(EsafMsg::Info) << "Reading new-style CONEX file version 2.x " ;
143
144                // Somewhere before 2.4 the FirstInteraction tree was renamed to
145                // LeadingInteractions by Ralf Ulrich. The version bump to 2.4 happend
146                // later. So there might be rare case of 2.3 outputs with wrongly labled
147                // trees.
148                if (OutputVersion < 2.4){
149                  Msg(EsafMsg::Info) << "with FirstInteraction tree" << MsgDispatch;
150                        fConexFirstInteraction = (TTree*)fInputFile->Get("FirstInteraction");
151                } else {
152                  Msg(EsafMsg::Info) << "with LeadingInteractions tree" << MsgDispatch;
153                        fConexFirstInteraction = (TTree*)fInputFile->Get("LeadingInteractions");
154                }
155
156                fConexShower->SetBranchAddress("nX",&nX); // number of points in slant depth
157                fConexShower->SetBranchAddress("N",N);     // particles(X)
158                fConexShower->SetBranchAddress("dEdX",dEdX); // Energy deposit(X)
159                fConexShower->SetBranchAddress("Electrons",Electrons);
160                fConexShower->SetBranchAddress("X",X);       // slant depth
161                fConexShower->SetBranchAddress("H",H);       // Height
162                fConexShower->SetBranchAddress("D",D);       // Distance
163
164                // depth step is not defined any more, we need to derive it from data
165
166
167
168                fConexShower->SetBranchAddress("lgE",&fitpars[0]);
169                fConexShower->SetBranchAddress("zenith",&fitpars[1]);
170                fConexShower->SetBranchAddress("azimuth",&fitpars[10]);
171                fConexShower->SetBranchAddress("Xmax",&fitpars[9]);
172
173                fConexFirstInteraction->SetBranchAddress("depth",depth);
174                fConexFirstInteraction->SetBranchAddress("height",height);
175                fConexFirstInteraction->SetBranchAddress("idInt",IdInt1);
176                fConexFirstInteraction->SetBranchAddress("nPart",&nPart);
177                fConexFirstInteraction->SetBranchAddress("nInt",&nInt);
178
179                // FIXME: Get X binning for each shower, find nicer way to get binning.
180                fConexShower->GetEntry(0);
181                delX = X[1] - X[0];
182
183                Msg(EsafMsg::Info) << "CONEX file X binning is " << delX << " g/cm^2" << MsgDispatch;
184        } else {
185                Msg(EsafMsg::Panic) << "CONEX file version not known!" << MsgDispatch;
186        }
187
188
189
190}
191
192//_____________________________________________________________________________
193Bool_t ConexFileShowerGenerator::LoadTrack() {
194  //
195  // load the shower track. Convert it from Conex format to ESAF
196  //
197
198  Msg(EsafMsg::Info) << "Reading shower id " << fCurrentEvent << MsgDispatch;
199
200  fConexHeader->GetEntry(0);
201  fConexShower->GetEntry(fCurrentEvent);
202  fConexFirstInteraction->GetEntry(fCurrentEvent);
203 
204  if(!fTrack) fTrack = new ShowerTrack();
205  else fTrack->Reset();
206 
207  if (!pMom) pMom = new TLorentzVector;
208  else       pMom->SetXYZT(0,0,0,0);
209 
210  fTrack->fEnergy = pow(10.,fitpars[0]) *eV; 
211
212  if (fForceTheta){
213    fTrack->fTheta = fTheta;
214  } else {
215    fTrack->fTheta = fitpars[1] * TMath::DegToRad();
216  }
217
218  if (fForcePhi){
219    fTrack->fPhi = fPhi;
220  } else {
221    fTrack->fPhi = fitpars[10] * TMath::DegToRad();  // Definition phi needs to be checked.
222  }
223  fTrack->fX1 = depth[0]/(cm2/g);  // Need to be checked
224
225
226// KS (090331)
227//      Temporarily theta and phi angles to behand-in inputted
228
229//  fTrack->fTheta  =  45. * TMath::DegToRad();  // Commented out so that ONEX Theta is not overwritten
230  fTrack->fPhi    =  23. * TMath::DegToRad();
231
232//      Temporarily theta and phi angles to behand-in inputted
233  core_x = fCoreX * km;
234  core_y = fCoreY * km;
235  core_z = TMath::Sqrt(((6371. + fCoreH)*km)*((6371. + fCoreH)*km)-core_x*core_x-core_y*core_y)-6371.*km;
236
237
238
239// Note that the Earth's radii in CONEX is 6378.14 km.     
240
241
242  // Shower axis with respect to core
243  cnx_l = TMath::Sin(fTrack->fTheta)*TMath::Cos(fTrack->fPhi);
244  cnx_m = TMath::Sin(fTrack->fTheta)*TMath::Sin(fTrack->fPhi);
245  cnx_n = TMath::Cos(fTrack->fTheta);
246
247    //
248
249  cout << core_x << " " << core_y << " " << core_z << " km" << endl;
250 
251
252
253  // Vector of zenith direction with respect to the core location
254  TVector3 core_zenith(core_x/6371./km,core_y/6371./km,(core_z+6371.*km)/6371./km);
255 
256  cout << core_x/6371./km << " " << core_y/6371./km << " " << (core_z+6371.*km)/6371./km << "  corepos_vector " << endl;
257 
258  Float_t rotate_x,rotate_y; 
259 
260  rotate_y = TMath::ATan2(core_x,core_z+6371.*km);
261  rotate_x = TMath::ATan2(core_y,core_z+6371.*km);
262 
263  //  tmp_lat = TMath::ASin(Core_y/6371./km);
264
265  cout << rotate_x * TMath::RadToDeg() << " " << rotate_y * TMath::RadToDeg() <<  " rotation " << endl;
266 
267 
268  // Primary particle trajectory vector (skyward)
269  geo_l = (cnx_l * TMath::Cos(rotate_y) + cnx_n * TMath::Sin(rotate_y));
270  geo_m = (cnx_m) * TMath::Cos(rotate_x) - (- cnx_l * TMath::Sin(rotate_y) + cnx_n * TMath::Cos(rotate_y)) * TMath::Sin(rotate_x);
271  geo_n = (cnx_m) * TMath::Sin(rotate_x) + (- cnx_l * TMath::Sin(rotate_y) + cnx_n * TMath::Cos(rotate_y)) * TMath::Cos(rotate_x);
272 
273#ifndef DEBUG
274  cout << cnx_l << " " << cnx_m << " " << cnx_n << " " << cnx_l*cnx_l+cnx_m*cnx_m+cnx_n*cnx_n << " " << fTrack->fTheta*TMath::RadToDeg()  << " deg " << fTrack->fPhi*TMath::RadToDeg() << " deg " << endl;
275 
276 
277  cout << geo_l << " " << geo_m << " " << geo_n << " " << geo_l*geo_l+geo_m*geo_m+geo_n*geo_n << " " << TMath::ACos(geo_n)*TMath::RadToDeg() << " " << TMath::ATan2(geo_m,geo_l)*TMath::RadToDeg()  << " final_vector " << endl;
278 
279  cout << TMath::ACos((geo_l*core_x+geo_m*core_y+geo_n*(core_z+6371.*km))/(6371.*km))*TMath::RadToDeg() << " " << (geo_l*core_x+geo_m*core_y+geo_n*(core_z+6371.*km))/(6371.*km) << " inner prod" << endl;
280 
281  //exit(-1);
282  #endif
283
284 
285  TVector3 dir(geo_l,geo_m,geo_n);
286 
287  fTrack->fDirVers = -dir;
288   
289  fTrack->fEthrEl = 0;
290 
291
292//  pMom->SetXYZT(-fTrack->fEnergy*TMath::Sin(fTrack->fTheta)*TMath::Cos(fTrack->fPhi),-fTrack->fEnergy*TMath::Sin(fTrack->fTheta)*TMath::Sin(fTrack->fPhi),-fTrack->fEnergy*TMath::Cos(fTrack->fTheta),fTrack->fEnergy);
293 
294  pMom->SetXYZT(-geo_l*fTrack->fEnergy,-geo_m*fTrack->fEnergy,-geo_n*fTrack->fEnergy,fTrack->fEnergy);
295 
296
297
298
299
300  //  cout << fitpars[9] << " " << fitpars[1] << " " << fitpars[10] << " " << fitpars[0] << endl;
301  cout << geo_l << " " <<  geo_m << " " << geo_n <<  " " <<  " ve" << endl;
302  cout << core_x << " " <<  core_y << " " << core_z <<  " " <<  " ve" << endl;
303
304
305   
306
307  // Age parameter setting (not yet fully working)
308  for (int i = 0 ; i < nX ; i++ ) {
309    ShowerStep s = GetShowerStep(i);
310    fTrack->Add(s);
311  }
312  return kTRUE;
313}
314
315//_____________________________________________________________________________
316const ShowerStep& ConexFileShowerGenerator::GetShowerStep(Int_t i) {
317  fStep.Clear();
318 
319  if (i==0) {
320
321#ifdef ZERO_ASSUMPTION   
322    prevdepth = X[0];
323   
324    Double_t L = D[0]*m; // from CONEX data (distance from the core on Earth along shower axis)
325   
326#endif
327
328    prevdepth = depth[0];
329
330    cout << depth[0] << endl;
331
332    //re-compute 3D initial position
333    Double_t ct = - TMath::Cos(fTrack->GetTheta());
334    Double_t y1 = height[0]/6378140.; // height/EarthRadius
335    Double_t L = (ct + TMath::Sqrt(ct*ct + 2*y1 + y1*y1))*6378140.*m; // track length
336   
337//    Double_t x = L*TMath::Sin(fTrack->GetTheta())*TMath::Cos(fTrack->GetPhi());
338//    Double_t y = L*TMath::Sin(fTrack->GetTheta())*TMath::Sin(fTrack->GetPhi());
339//    Double_t z = L*TMath::Cos(fTrack->GetTheta());
340
341
342    Double_t x = core_x + L * geo_l;
343    Double_t y = core_y + L * geo_m;
344    Double_t z = core_z + L * geo_n;
345   
346    cout << x << " " << y << " " << z << endl;
347    cout << geo_l << " " <<  geo_m << " " << geo_n <<  " " <<  L << endl;
348
349    prevpoint.SetXYZ(x,y,z);
350    fTrack->fInitPos.SetXYZ(x,y,z);
351   
352    // cerr << "showerstep " << y1 << " " << prevdepth << " " <<  height[nInt-1] << " " << i << " " <<currentpoint.Zv()/km << "  " << X[i-1] << height[i] << " " << L <<  " " << ct << endl;
353  }
354  else{
355    prevdepth = X[i-1];
356  }
357   
358  currentpoint = prevpoint;
359 
360//  if (currentpoint.Zv()/km < 0) return fStep;   // In case core is not at (0,0)
361  TVector3 dir = pMom->Vect().Unit(); 
362 
363  Double_t depth_step=0;
364  while (kTRUE) {
365    Double_t step = 10*m;
366    currentpoint += step*dir;
367    if (currentpoint.Zv()/km<0) break;
368    depth_step += Atmosphere::Get()->Air_Density(currentpoint)*step;
369    if (depth_step>=delX*g/cm2) break;
370  }
371 
372  fStep.fTimei= (fTrack->fInitPos - prevpoint).Mag() / EConst::Clight();
373  fStep.fTimef= (fTrack->fInitPos - currentpoint).Mag() / EConst::Clight();
374  fStep.fXi = prevdepth;
375  fStep.fXf = X[i];
376  fStep.fXYZi = prevpoint;
377  fStep.fXYZf = currentpoint;
378  fStep.fNcharged = N[i];
379  fStep.fNelectrons = Electrons[i];
380 
381  //  cout << X[i] << endl;
382 
383  //  fStep.fEnergyLoss = dEdX[i] * Xdel * GeV;
384 
385  // Temporary age parameter approximation
386 
387  if ( fStep.fXi > depth[0] ) fStep.fAgei =  2. / (1 + (fitpars[9]-depth[0]) /(fStep.fXi-depth[0]) );
388  else fStep.fAgei = 0.;
389 
390  if ( fStep.fXf > depth[0] ) fStep.fAgef =  2. / (1 + (fitpars[9]-depth[0]) /(fStep.fXf-depth[0]) );
391  else fStep.fAgef = 0.;
392 
393  prevpoint = currentpoint;
394  return fStep;
395}
396
397//_____________________________________________________________________________
398void ConexFileShowerGenerator::Reset() {
399 
400}
401
402//______________________________________________________________________________
403PhysicsData* ConexFileShowerGenerator::Get() {
404
405  //
406  // Get a new event
407  //
408 
409 
410  // build and return a track
411 
412  if ( !fInputFile ) OpenFile(); 
413 
414  while (1) {
415   
416    if ( fTrack ) {
417      delete fTrack;
418      fTrack=NULL;
419    }
420    if ( fTruth ) {
421      delete fTruth;
422      fTruth=NULL;
423    }
424   
425    if (fCurrentEvent >= fFirstEvent && fCurrentEvent<100) {
426      if ( LoadTrack() ){
427       
428        MsgForm(EsafMsg::Info,"ConexFile event %d", fCurrentEvent);
429        MsgForm(EsafMsg::Info,
430                "ShowerTrack produced : Energy = %.2E MeV, Theta= %.2f deg, Phi= %.2f deg",
431                fTrack->GetEnergy(), fTrack->GetTheta()*RadToDeg(),fTrack->GetPhi()*RadToDeg());
432        MsgForm(EsafMsg::Info,
433                "         First interaction X1 = %.3f at : (%.2f,%.2f,%.2f) km gives a shower with %d steps",
434                fTrack->GetX1()*cm2/g,fTrack->GetInitPos().X()/km,fTrack->GetInitPos().Y()/km,
435                fTrack->GetInitPos().Z()/km, fTrack->GetNumStep());
436      } else {
437        Msg(EsafMsg::Warning) << "ConexFileShowerGenerator: End of file "
438                              << fFileName << " reached? " << MsgDispatch;
439        fTrack = (ShowerTrack*)NULL;
440      }
441      fCurrentEvent++;
442      break;
443    } 
444   
445    fCurrentEvent++;
446   
447  }
448  fTrack->CheckTrack();
449  return fTrack;
450 
451}
452
453//KS
454//______________________________________________________________________________
455MCTruth* ConexFileShowerGenerator::GetTruth() {
456  //
457  // Method returning the Truth object
458  //
459 
460  // init
461  if(!fTruth) fTruth = new MCTruth();
462 
463  // can be filled even if showertrack has no steps
464  fTruth->SetEnergy(fTrack->GetEnergy());
465  fTruth->SetThetaPhi(fTrack->GetTheta(),fTrack->GetPhi());
466 
467// temporarily
468  Int_t particode = Int_t(floor(1+0.5)); // proton (default in CONEX)
469  fTruth->SetParticle(particode);
470 
471  //    fTruth->SetTOAImpact(fTOAImpact);
472  //    fTruth->SetEarthImpact(fEarthImpact);
473 
474  // filled only if steps exist   
475  if (fTrack->Size() > 0) {
476    const ShowerStep& s = fTrack->GetLastStep();
477    fTruth->SetEarthAge(s.GetAgef());
478    // get the shower maximum
479    Float_t MaxElectrons(0);
480    Int_t MaxIndex(-1);
481    for (Int_t i(0); i < Int_t(fTrack->Size()) ; i++) {
482      const ShowerStep& s = fTrack->GetStep(i);
483      if (MaxElectrons < s.GetNelectrons()) {
484        MaxIndex = i;
485        MaxElectrons = s.GetNelectrons();
486      }
487    }
488    if (MaxIndex != -1) {
489      const ShowerStep& s = fTrack->GetStep(MaxIndex);
490      fTruth->SetShowerMax(0.5*(s.GetXYZi() + s.GetXYZf()),0.5*(s.GetXi() + s.GetXf())*(g/cm2));
491    }
492  }
493  return fTruth;
494}
495
496
Note: See TracBrowser for help on using the repository browser.