source: JEM-EUSO/esaf_cc_at_lal/packages/simulation/generators/showers/src/CorsikaFileShowerGenerator.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: 22.8 KB
Line 
1// $Id: CorsikaFileShowerGenerator.cc 2930 2011-06-16 20:58:10Z mabl $
2// Author: M.Dattoli   2004/08/31
3
4/*****************************************************************************
5 * ESAF: Euso Simulation and Analysis Framework                              *
6 *                                                                           *
7 *  Id: CorsikaFileShowerGenerator                                           *
8 *  Package: ShowerGenerators                                                *
9 *  Coordinator: Sergio.Bottai, Dmitry.Naumov                                *
10 *                                                                           *
11 *****************************************************************************/
12
13//_____________________________________________________________________________
14//   Corsika File Shower Generator
15//   =============================
16//
17//   <extensive description>
18//
19//   Config file parameters
20//   ======================
21//   
22//   FlagThinning [string]:  flag for thinning option;
23//                           can be t (thinning enabled)
24//                           or f (thinning disabled)
25//
26//   HowManyShowers [int] : number of showers simulated in corsika run
27//
28//   ParticleFileName [string]: name of corsika normal particle output file
29//
30//   LongFileName [string]: name of corsika longitudinal output file
31//
32
33#include "CorsikaFileShowerGenerator.hh"
34#include "MCTruth.hh"
35#include <ctime>
36#include <cstdlib>
37#include <cmath>
38 
39using namespace sou;
40using namespace TMath;
41using namespace EConst;
42
43ClassImp(CorsikaFileShowerGenerator)
44
45//______________________________________________________________________________
46CorsikaFileShowerGenerator::CorsikaFileShowerGenerator(): 
47           EventGenerator("Corsikafile"),EsafMsgSource(), fTrack(0), 
48           fParticleFile(NULL), fLongFile(NULL), fTruth(0) 
49    //
50    // Constructor
51    //
52{
53    fParticleFileName = "";
54    fLongFileName = "";
55}
56
57
58//______________________________________________________________________________
59CorsikaFileShowerGenerator::~CorsikaFileShowerGenerator() 
60    //
61    // Destructor
62    //
63{
64    if ( fTrack )  delete fTrack;
65    if ( fTruth )  delete fTruth;
66        fTrack=NULL;
67        fTruth=NULL;
68}
69
70
71
72//______________________________________________________________________________
73void CorsikaFileShowerGenerator::ReplaceInputFile(const char* pfn,const char* lfn) 
74    //
75    // Replace input files
76    //
77{
78
79    // Particle file
80    if ( fParticleFile )
81         Msg(EsafMsg::Panic)<<"Can not change name when particle file is already open!"<<MsgDispatch;
82    if (!pfn) 
83    {
84        Msg(EsafMsg::Panic)<<" Particle file name ignored!" <<MsgDispatch;
85        return;
86    }
87    Conf()->ReplaceStr("CorsikaFileShowerGenerator.ParticleFileName",pfn);
88    Msg(EsafMsg::Info)<<"Input file name changed to" << pfn <<MsgDispatch;
89
90    //Longitudinal file: info about energy loss, n electrons, ncharged, ncherenkov, depth
91    if ( fLongFile )
92        Msg(EsafMsg::Panic)<<"Can not change name when particle file is already open!"<<MsgDispatch;
93    if (!lfn) 
94    { 
95        Msg(EsafMsg::Panic)<<" Particle file name ignored!" <<MsgDispatch;
96        return;
97    }
98    Conf()->ReplaceStr("CorsikaFileShowerGenerator.LongFileName",lfn);
99    Msg(EsafMsg::Info)<<"Input file name changed to" << lfn <<MsgDispatch;
100
101}
102
103
104
105//______________________________________________________________________________
106void CorsikaFileShowerGenerator::Open() 
107    //
108    // Open file
109    //
110{
111
112  //particle file
113    if ( Conf()->IsDefined( "CorsikaFileShowerGenerator.ParticleFileName" ) )
114        fParticleFileName = Conf()->GetStr("CorsikaFileShowerGenerator.ParticleFileName");
115
116    else
117        Msg(EsafMsg::Panic)<<"Unknown input particle file name in CorsikaFileShowerGenerator"<<MsgDispatch;
118
119    fParticleFile = fopen( fParticleFileName.c_str(), "rb" );
120
121    if( fParticleFile == NULL ) 
122        Msg(EsafMsg::Panic)<<"Can't open file "+fParticleFileName << MsgDispatch;
123    else
124        Msg(EsafMsg::Info) << "File: " << fParticleFileName << " successfully opened " << MsgDispatch;
125
126  //longitudinal file
127    if ( Conf()->IsDefined( "CorsikaFileShowerGenerator.LongFileName" ) )
128        fLongFileName = Conf()->GetStr("CorsikaFileShowerGenerator.LongFileName");
129
130    else
131        Msg(EsafMsg::Panic)<<"Unknown input longitudinal ditribution file name in CorsikaFileShowerGenerator"<<MsgDispatch;
132
133    fLongFile = fopen( fLongFileName.c_str(), "rb" );
134
135    if( fLongFile == NULL ) 
136        Msg(EsafMsg::Panic)<<"Can't open file "+fLongFileName << MsgDispatch;
137    else
138        Msg(EsafMsg::Info) << "File: " << fLongFileName << " successfully opened " << MsgDispatch;
139
140
141}
142
143
144//______________________________________________________________________________
145void CorsikaFileShowerGenerator::Close() 
146    //
147    // Close file
148    //
149{
150
151    if( fParticleFile ) fclose(fParticleFile ) ;
152    if( fLongFile  ) fclose(fLongFile ) ;
153
154   
155}
156
157
158//______________________________________________________________________________
159PhysicsData* CorsikaFileShowerGenerator::Get() 
160{
161    //
162    // Get a new event
163    //
164
165    // build and return a track
166    if ( !fParticleFile || !fLongFile) Open();
167 
168    if ( Conf()->IsDefined( "CorsikaFileShowerGenerator.FlagThinning" ) )
169        thinning = Conf()->GetStr("CorsikaFileShowerGenerator.FlagThinning");
170    else Msg(EsafMsg::Panic) <<"thinning flag missing in config file" << MsgDispatch;
171    if (thinning=="f") {
172        n1 = 272;
173        nmax = 7;
174        Msg(EsafMsg::Info)<<"Thinning disabled"<<MsgDispatch;
175        }
176
177        else if (thinning=="t") {
178        n1 = 311;
179        nmax = 8;
180        Msg(EsafMsg::Info)<<"Thinning enabled"<<MsgDispatch;
181        }
182
183        n2 = n1 + 1;
184 
185        ConfigFileParser *pConfig = Config::Get()->GetCF("Generators","CorsikaFileShowerGenerator");
186        howmany = pConfig->GetNum("CorsikaFileShowerGenerator.HowManyShowers");
187
188        //random choice of event to be processed among simulated events
189        srand(time(0));
190        event =(rand()%(Int_t)(howmany))+1;
191        Msg(EsafMsg::Debug) << "Processing event n : " << event << MsgDispatch;
192 
193
194       if ( LoadTrack(event) ) {
195        Msg(EsafMsg::Debug) << "Load Track Event" << MsgDispatch;   
196        return fTrack;
197        }
198       else {
199        Msg(EsafMsg::Panic) << "End of file " << fParticleFileName << " reached? " << MsgDispatch;
200        return (ShowerTrack*)0;
201     }
202 
203       return fTrack;   
204       }
205
206//______________________________________________________________________________
207void CorsikaFileShowerGenerator::ReadParticleHeader() 
208{
209    //
210    // Read particle file header.
211    //
212
213  Msg(EsafMsg::Debug) << "Reading particle header.." << MsgDispatch;
214    while(!feof(fParticleFile))
215    {
216        //Reading BlockHeader (273,21) or (312,21)
217
218        nin = fread(&r,4,1,fParticleFile);
219        for (k=1; k<=21;k++)
220        {
221            pos = ftell(fParticleFile); 
222
223            //Reading SubBlockHeader   
224            nin = fread(&header,sizeof(header),1,fParticleFile);
225
226            //Reading data in SubBlock; Storing them into map
227            if (header[0]=='R'&&header[3]=='H') 
228            {
229
230                float data_runh[n1];
231                for (i=0; i<n1; i++) {data_runh[i]=0;}
232                //Reading RUNH Block (once for run)
233                nin = fread (&data_runh, sizeof(data_runh), 1, fParticleFile);
234                for (i=0; i<=4; i++) {
235                aatm[i]=data_runh[253+i];
236                batm[i]=data_runh[258+i];
237                catm[i]=data_runh[263+i];
238                }
239            }  // end of if "runh"
240
241            else if (header[0]=='R'&&header[3]=='E') 
242            {
243                float data_rune[n1];
244
245                for (i=0; i<n1; i++) {data_rune[i]=0;}
246                //reading RUNE Block (once per run)
247                nin = fread (&data_rune, sizeof(data_rune), 1,fParticleFile );
248            } // end of if "rune"
249
250            else if (header[0]=='E'&&header[3]=='H') 
251            {
252                float data_evth[n1];
253
254                for (i=0; i<n1; i++) {data_evth[i]=0;}
255
256                //Reading EVTH Block (1 for event)
257                nin = fread (&data_evth, sizeof(data_evth), 1,fParticleFile );
258                ev_num= (Int_t)(data_evth[0]);
259
260cout <<"LOOK "<< ev_num << " " << event << endl;
261//                 if (ev_num==event)
262                if (1) 
263                {
264                Msg(EsafMsg::Debug) << "Reading event n " << ev_num << MsgDispatch;
265                    z_firstint = data_evth[5];  // cm
266                    Msg(EsafMsg::Debug)<< "z_firstint: " << z_firstint <<MsgDispatch; 
267                    thick_firstint = GetThickness(abs(z_firstint)); //g/cm^2
268                    Msg(EsafMsg::Debug)<< "thick_firstint: " << thick_firstint <<MsgDispatch; 
269                    epolangd = data_evth[9]*180/Pi();
270                    eaziangd = data_evth[10]*180/Pi();
271                    eldircos = Sin(data_evth[9])*Cos(data_evth[10]);
272                    emdircos = Sin(data_evth[9])*Sin(data_evth[10]);
273                    endircos = Cos(data_evth[9]);
274                    energy = data_evth[2];
275                    Msg (EsafMsg::Debug) << "epolangd: " << epolangd
276                                        << "\n eaziangd: " << eaziangd
277                                        << "\n eldircos: " << eldircos
278                                        << "\n emdircos: " << emdircos
279                                        << "\n endircos: " << endircos
280                                        << "\n energy: " << energy << MsgDispatch;
281                   
282                }  //end of if (ev_num==event)
283            } //end of if "evth"
284
285            else if (header[0]=='E'&&header[3]=='E') 
286            {
287                float data_evte[n1];
288
289                for (i=0; i<n1; i++) {data_evte[i]=0;}
290                //Reading EVTE Block (1 for event)
291                nin = fread (&data_evte, sizeof(data_evte), 1, fParticleFile);
292
293            } // end of if "evte"
294
295            else if (header[0]=='L'&&header[3]=='G') 
296            {
297                float data_long[n1];
298                for (i=0; i<n1; i++) {data_long[i]=0;}
299
300                //Reading LONG Block
301                nin = fread (&data_long, sizeof(data_long), 1, fParticleFile);
302            } // end of if "long"
303
304            else 
305            {
306                fseek(fParticleFile,pos,SEEK_SET);
307                float data_block[n2];
308                for (i=0; i<n2; i++) {data_block[i]=0;}
309
310                //Reading DATA BLOCK
311                nin = fread (&data_block, sizeof(data_block), 1,fParticleFile );
312
313            } // end of else... DATA BLOCK
314
315        } // end of for (k=1; k<=21;k++)
316
317        k=1;
318
319        //Reading BlockEnd  (273,21) o (312,21)
320        nin = fread(&r,4,1,fParticleFile);
321    }
322
323}
324
325//______________________________________________________________________________
326void CorsikaFileShowerGenerator::ReadLong()
327    //
328    // Read long file
329    //
330{
331
332    // These dummies are just for silencing gcc warning for unused fscanf/fgets
333    // results.
334    // FIXME: Only temporal, but better than 18 warnings
335
336    char *fgetDummy;
337    int fscanfDummy;
338
339Msg(EsafMsg::Debug) << "Reading Long File.." << MsgDispatch;
340                for ( i=0; i<10; i++) 
341                { en[i]=0;
342                  part[i]=0;
343                }
344    nshow=0;
345    pos = ftell(fLongFile); 
346    while (!feof(fLongFile) && nshow<=event)
347    {
348       fscanfDummy = fscanf(fLongFile,"%*s%*s%*s%f%*s%*s%*s%f%*s%*s%*s%i",&nstep,&step_long,&nshow);
349       fgetDummy = fgets(longitudinal,sizeof(longitudinal),fLongFile);
350       pos=ftell(fLongFile);
351       Msg(EsafMsg::Debug)<< "nstep, step_long:  " << nstep << "; \t"<< step_long << MsgDispatch;
352       fgetDummy = fgets(longitudinal,sizeof(longitudinal),fLongFile);
353                                                                                                                           
354       if (nshow!=event)
355         {
356           for (j=1; j<= 2*nstep+12; j++ ){
357               fgetDummy = fgets(longitudinal,sizeof(longitudinal),fLongFile);
358             pos = ftell(fLongFile);
359             fseek(fLongFile,pos,SEEK_SET);
360           }
361           fseek(fLongFile,pos,SEEK_SET);
362           nshow = nshow+1;
363         }
364                                                                                                                           
365       else
366         {
367           for ( j=1; j <= nstep; j++)
368          {
369             for ( i=0; i<10; i++)
370              {
371               fscanfDummy = fscanf(fLongFile,"%f",&part[i]);
372               pos=ftell(fLongFile);
373               fseek(fLongFile,pos,SEEK_SET);
374              } // reading 1 raw
375             particle.insert(pair<int,float>(j,part[0]));
376             particle.insert(pair<int,float>(j,part[0] - step_long));
377             particle.insert(pair<int,float>(j,part[3] + part[2]));
378             particle.insert(pair<int,float>(j,part[7]));
379             particle.insert(pair<int,float>(j,part[9]));
380         } // reading block
381                                                                                                                           
382             fgetDummy = fgets(longitudinal,sizeof(longitudinal),fLongFile);
383             fgetDummy = fgets(longitudinal,sizeof(longitudinal),fLongFile);
384             fgetDummy = fgets(longitudinal,sizeof(longitudinal),fLongFile);
385                                                                                                                           
386           for ( j=1; j <=nstep; j++)
387            {
388              for ( i=0; i<10; i++)
389               { fscanfDummy = fscanf(fLongFile,"%f",&en[i]);
390                } // reading 1 raw
391              particle.insert(pair<int,float>(j,en[2]));
392            } // reading block
393           
394              fgetDummy = fgets(longitudinal,sizeof(longitudinal),fLongFile);
395              fgetDummy = fgets(longitudinal,sizeof(longitudinal),fLongFile);
396              fgetDummy = fgets(longitudinal,sizeof(longitudinal),fLongFile);
397
398              pos=ftell(fLongFile);
399              fseek(fLongFile,pos,SEEK_SET);
400              fscanfDummy = fscanf(fLongFile,"%*s%*s%f%f%f%*f%*f%*f",&Nemax,&to,&tmax);
401              for (int i = 1; i<=8; i++) {
402                  fgetDummy = fgets(longitudinal,sizeof(longitudinal),fLongFile);
403                }
404
405             nshow=nshow+1;
406         } // end of else
407                                                                                                                           
408       pos=ftell(fLongFile);
409     }
410                                                                                                                           
411
412}
413
414
415//______________________________________________________________________________
416bool CorsikaFileShowerGenerator::LoadHeader()
417    //
418    // Fill header object
419    //
420{
421    fHeader.clear();
422    ReadParticleHeader();
423
424    x0comes = xi = 0;
425    y0comes = yi = 0;
426    if (z_firstint>=0) {
427            z0comes = z_firstint;}
428    else {
429        z0comes = - z_firstint;}
430    zmax = GetAltitude(tmax);
431
432    vertical_distance = (abs(z_firstint) - zmax);
433    coeff_f = GetCoeff(vertical_distance);
434   
435    xmax = GetXstep(coeff_f);
436    ymax = GetYstep(coeff_f);
437    float coeff = GetCoeff(z0comes);
438
439    fHeader.insert(pair<string,float>("EPRMENRG",energy));
440    fHeader.insert(pair<string,float>("EPRMDPTH",thick_firstint));
441    fHeader.insert(pair<string,float>("EX0COMES",x0comes)); 
442    fHeader.insert(pair<string,float>("EY0COMES",y0comes)); 
443    fHeader.insert(pair<string,float>("EZ0COMES",z0comes)); 
444    fHeader.insert(pair<string,float>("EPOLANGD",epolangd)); 
445    fHeader.insert(pair<string,float>("EAZIANGD",eaziangd));
446    fHeader.insert(pair<string,float>("ELDIRCOS",eldircos)); 
447    fHeader.insert(pair<string,float>("EMDIRCOS",emdircos));
448    fHeader.insert(pair<string,float>("ENDIRCOS",endircos));
449    fHeader.insert(pair<string,float>("EXDMPPOS",GetXstep(coeff)));
450    fHeader.insert(pair<string,float>("EYDMPPOS",GetYstep(coeff)));
451    fHeader.insert(pair<string,float>("EZDMPPOS",0.));
452    fHeader.insert(pair<string,float>("EXMAXPOS",xmax)); 
453    fHeader.insert(pair<string,float>("EYMAXPOS",ymax)); 
454    fHeader.insert(pair<string,float>("EZMAXPOS",zmax)); 
455    return true;
456
457}
458
459//______________________________________________________________________________
460bool CorsikaFileShowerGenerator::LoadTrack(Int_t event)
461    //
462    // Load one event from file
463    //
464{
465
466    ReadLong();
467    if (!LoadHeader()) return false;
468
469
470    fTrack = new ShowerTrack();
471
472    fTrack->fEthrEl = 0.;
473    fTrack->fEnergy = fHeader["EPRMENRG"]*GeV;   //in eV
474    fTrack->fTheta  = fHeader["EPOLANGD"]*deg;
475    fTrack->fPhi    = fHeader["EAZIANGD"]*deg;
476    fTrack->fX1     = fHeader["EPRMDPTH"]*g/cm2;
477   
478    fTrack->fDirVers[0] = fHeader["ELDIRCOS"];
479    fTrack->fDirVers[1] = fHeader["EMDIRCOS"];
480    fTrack->fDirVers[2] = fHeader["ENDIRCOS"];
481   
482    fTrack->fInitPos[0] = fHeader["EX0COMES"]*cm;
483    fTrack->fInitPos[1] = fHeader["EY0COMES"]*cm;
484    fTrack->fInitPos[2] = fHeader["EZ0COMES"]*cm;
485cout << "first point " << fTrack->fInitPos[2] << endl;
486
487    if(fHeader["HITGROUN"]==0) 
488    { 
489        fTrack->fHitGround = false;
490    }
491    else 
492    {
493        fTrack->fHitGround = true;
494    }
495
496
497
498    coeff_i = 0;
499    xi=0;
500    yi=0;
501    zi = abs(z_firstint);
502    timei = 0;
503    for(j=1; j<=nstep; j++)
504    {
505
506
507        int count=0;
508        for (pnum=particle.lower_bound(j); pnum!=particle.upper_bound(j); ++pnum)
509        {
510            step[count] = pnum->second;
511            count++;
512        }
513
514
515// step_bound (g/cm^2)
516        fStep.fXi = (step[1]+ thick_firstint)*gram/cm2;
517        fStep.fXf = (step[1]+ thick_firstint+step_long)*gram/cm2;
518
519// step_bound (cartesian coordinates)
520       
521        thickness = step_long*endircos;
522        zf = GetAltitude(step[1]+ thick_firstint+thickness);
523        cout << zf << " " << step[1] << " " << thick_firstint << " " << thickness << " " << step_long << " " << endircos << endl;                                                                       
524        vertical_distance = abs(zi-zf);
525        coeff_f = GetCoeff(vertical_distance);
526        xf = GetXstep(coeff_f);
527        yf = GetYstep(coeff_f);
528
529       if (zi>0) {
530        fStep.fXYZi[2] = zi*cm;}
531        else {
532        fStep.fXYZi[2] = 0;}
533
534        if (zf>0) {
535        fStep.fXYZf[2] = zf*cm;}
536        else {
537        fStep.fXYZf[2] = 0;}
538 if (fStep.fXYZi[2]!=0 || fStep.fXYZf[2]!=0 ) 
539{
540//        fStep.fXYZf[2] = zf;
541        fStep.fXYZi[0] = xi*cm;
542        fStep.fXYZf[0] = xf*cm;
543        fStep.fXYZi[1] = yi*cm;
544        fStep.fXYZf[1] = yf*cm;
545// time
546        distance = sqrt(abs(pow(xi,2) + pow(yi,2) + pow(zi,2) - (pow(xf,2) + pow(yf,2) + pow(zf,2))));
547
548        fStep.fTimei = timei*nanosecond; 
549        timef = GetStepTime(distance) + timei;
550        fStep.fTimef = timef*nanosecond;
551        fStep.fNelectrons = step[2];
552        fStep.fNcharged = step[3];
553        fStep.fNcherenkov = step[4];
554        fStep.fEloss = step[5];
555
556// age
557
558        fStep.fAgei = GetAge(step[1] + thick_firstint);
559        fStep.fAgef = GetAge(step[1] + thick_firstint+step_long);
560
561        fTrack->Add(fStep);
562        coeff_i = coeff_f;
563        xi = xf;
564        yi = yf;
565        zi = zf;
566        timei = timef;
567} 
568
569
570   }
571
572    return true;
573
574}
575
576//______________________________________________________________________________
577MCTruth* CorsikaFileShowerGenerator::GetTruth() 
578    //
579    // Get pointer to MonteCarlo Truth
580    //
581{
582
583
584    if (!fTrack->Size()) return NULL;
585    if (!fTruth) {
586         fTruth = new MCTruth();}
587
588    fTruth->SetEnergy(fHeader["EPRMENRG"]*GeV);   //in eV
589    fTruth->SetThetaPhi(fHeader["EPOLANGD"]*deg,fHeader["EAZIANGD"]*deg);
590    EVector v;
591    v[0] = fHeader["EX0COMES"]*cm;
592    v[1] = fHeader["EY0COMES"]*cm;
593    v[2] = fHeader["EZ0COMES"]*cm;
594    fTruth->SetFirstInt(v, fHeader["EPRMDPTH"]);
595    v[0] = fHeader["EXDMPPOS"]*cm;
596    v[1] = fHeader["EYDMPPOS"]*cm;
597    v[2] = fHeader["EZDMPPOS"]*cm;
598                                                                                                                             
599                                                                                                                             
600    // in case no value for ShowerMax is present, put arbitrary number 0
601                                                                                                                             
602    v[0] = 0;
603    v[1] = 0;
604    v[2] = 0;
605    map<string,double>::iterator showmax;
606    showmax = fHeader.find(string("EXMAXPOS"));
607    if (showmax != fHeader.end()) v[0] = fHeader["EXMAXPOS"]*cm;
608    showmax = fHeader.find(string("EYMAXPOS"));
609    if (showmax != fHeader.end()) v[1] = fHeader["EYMAXPOS"]*cm;
610    showmax = fHeader.find(string("EZMAXPOS"));
611    if (showmax != fHeader.end()) v[2] = fHeader["EZMAXPOS"]*cm;
612                                                                                                                             
613    fTruth->SetShowerMax(v,fHeader["ESHOWMAX"]);
614
615  return fTruth;
616                                                                                                                           
617}
618
619//______________________________________________________________________________
620float CorsikaFileShowerGenerator::GetAltitude(float thick)
621{
622 //FIXME: Better throw an error here if thickness is not valid.
623 float h=-HUGE;
624
625 if (thick >= GetThickness(4.e+05)) {
626        h = -catm[0]*log((thick - aatm[0])/batm[0]);}
627 if (thick >= GetThickness(10.e+05) && thick < GetThickness(4.e+05)) {
628        h = -catm[1]*log((thick - aatm[1])/batm[1]);}
629 if (thick >= GetThickness(40.e+05) && thick< GetThickness(10.e+05)) {
630        h = -catm[2]*log((thick - aatm[2])/batm[2]);}
631 if (thick >= GetThickness(100.e+05) && thick< GetThickness(40.e+05)) {
632        h = -catm[3]*log((thick - aatm[3])/batm[3]);}
633 if (thick <= GetThickness(100.e+05) ) {
634        h = (aatm[4]-thick)*catm[4]/batm[4];}
635 //cout << h << " " << thick << endl;
636 return h;
637}
638
639//______________________________________________________________________________
640float CorsikaFileShowerGenerator::GetThickness(float z1)
641{
642 float thick_f;
643 float z = abs(z1);
644 if (z<= 4.e+05) { 
645        thick_f = aatm[0]+batm[0]*exp(-z/catm[0]);}
646 else if (z> 4.e+05 && z<= 10.e+05) {
647        thick_f = aatm[1]+batm[1]*exp(-z/catm[1]);}
648 else if (z> 10.e+05 && z<= 40.e+05) {
649        thick_f = aatm[2]+batm[2]*exp(-z/catm[2]);}
650 else if (z> 40.e+05 && z<= 100.e+05) {
651        thick_f = aatm[3]+batm[3]*exp(-z/catm[3]);}
652 else if (z> 100.e+05) {
653        thick_f = aatm[4]-batm[4]*z/catm[4];}
654 else
655     throw invalid_argument("Out of parameterization for given hight!");
656
657 return thick_f;
658}
659
660//______________________________________________________________________________
661float CorsikaFileShowerGenerator::GetCoeff(float h)
662{
663 float coeff = h/endircos;
664 return coeff;
665}
666
667
668
669//______________________________________________________________________________
670float CorsikaFileShowerGenerator::GetXstep(float t)
671{
672 float x = xi + t*eldircos; 
673 return x;
674}
675
676
677//______________________________________________________________________________
678float CorsikaFileShowerGenerator::GetYstep(float t)
679{
680 float y = yi + t*emdircos; 
681 return y;
682}
683
684//______________________________________________________________________________
685float CorsikaFileShowerGenerator::GetStepTime(float h)
686{
687 float time = 0.1*h/Clight();
688 return time;
689}
690
691
692//______________________________________________________________________________
693float CorsikaFileShowerGenerator::GetAge(float h)
694{
695 float age;
696    float t = h-thick_firstint;
697    age = 3*t/(t + 2*tmax);
698 return age;
699}
700                                                                                                                           
701
Note: See TracBrowser for help on using the repository browser.