source: JEM-EUSO/esaf_lal/tags/v1_r0/esaf/packages/simulation/generators/showers/src/ShowerStep.cc @ 117

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

ESAF version compilable on mac OS
File size: 9.6 KB
Line 
1// $Id: ShowerStep.cc 2972 2011-08-01 08:34:15Z mabl $
2// Author: Alessandro Thea, Sergio Bottai, Dmitry Naumov   2004/12/08
3
4/*****************************************************************************
5 * ESAF: Euso Simulation and Analysis Framework                              *
6 *                                                                           *
7 *  Id: ShowerStep                                                           *
8 *  Package: Shower                                                          *
9 *  Coordinator: Sergio.Bottai, Dmitry.Naumov                                *
10 *                                                                           *
11 *****************************************************************************/
12
13//_____________________________________________________________________________
14//
15// ShowerStep is the universal container of Atmospheric Air Shower
16// produced by Ultra High Energy Cosmic Ray entering the atmosphere.
17// ShowerStep is a vector element of the ShowerTrack object which is
18// produced by interfaces to specific generators like CORSICA, AIRES, UNISIM, SLAST, etc
19// ShowerStep contains relevant information of EAS development at the given
20// position in space
21//
22//
23
24#include "ShowerStep.hh"
25#include "ShowerTrack.hh"
26#include "EConst.hh"
27#include <TF12.h>
28
29ClassImp(ShowerStep)
30using namespace sou;
31using namespace EConst;
32
33//______________________________________________________________________________
34ShowerStep::ShowerStep() :  EsafMsgSource(), fXi(0.),fXf(0.),fXYZi(0.,0.,0.),fXYZf(0.,0.,0.),
35    fTimei(0.),fTimef(0.),fAgei(0.),fAgef(0.),fNelectrons(0.),fStepID(-1) {
36    //
37    // Constructor
38    //
39    fNcharged         = 0;
40    fEloss            = 0;
41    fNcherenkov       = 0;
42    fHistoEnergy      = NULL;
43    fHistoLateral     = NULL;
44    fHistoEneAng      = NULL;
45    fHistoRadPhiEle   = NULL;
46    fHistoRadDTimeEle = NULL;
47    fHistoRadPhiEloss = NULL;
48    fHistoAngCher     = NULL;
49    fParentTrack      = NULL;
50}
51
52//______________________________________________________________________________
53ShowerStep::ShowerStep(const ShowerStep &s) {
54    //
55    // Copy Constructor
56    //
57    fXi               = s.fXi;
58    fXf               = s.fXf;
59    fXYZi             = s.fXYZi;
60    fXYZf             = s.fXYZf;
61    fTimei            = s.fTimei;
62    fTimef            = s.fTimef;
63    fAgei             = s.fAgei;
64    fAgef             = s.fAgef;
65    fNelectrons       = s.fNelectrons;
66    fNcharged         = s.fNcharged;
67    fEloss            = s.fEloss;
68    fNcherenkov       = s.fNcherenkov;
69    fParentTrack      = s.fParentTrack;
70    fStepID           = s.fStepID;
71
72    fHistoEnergy      = s.fHistoEnergy  ? (TH1F*)s.fHistoEnergy->Clone()  : 0;
73    fHistoLateral     = s.fHistoLateral ? (TH1F*)s.fHistoLateral->Clone() : 0;
74    fHistoEneAng      = s.fHistoEneAng  ? (TH2F*)s.fHistoEneAng->Clone()  : 0;
75    fHistoRadPhiEle   = s.fHistoRadPhiEle ? (TH2F*)s.fHistoRadPhiEle->Clone() : 0;
76    fHistoRadDTimeEle = s.fHistoRadDTimeEle ? (TH2F*)s.fHistoRadDTimeEle->Clone() : 0;
77    fHistoRadPhiEloss = s.fHistoRadPhiEloss ? (TH2F*)s.fHistoRadPhiEloss->Clone() : 0;
78    fHistoAngCher     = s.fHistoAngCher ? (TH1F*)s.fHistoAngCher->Clone() : 0;
79}
80
81//______________________________________________________________________________
82ShowerStep::~ShowerStep() {
83    //
84    // Destructor
85    //
86    Clear();
87}
88
89//______________________________________________________________________________
90void ShowerStep::Clear() {
91    //
92    // Clear the step
93    //
94
95
96    fXi = 0.;
97    fXf = 0.;
98    fXYZi.SetXYZ(0.,0.,0.);
99    fXYZf.SetXYZ(0.,0.,0.);
100
101    fTimei            = 0.;
102    fTimef            = 0.;
103    fAgei             = 0.;
104    fAgef             = 0.;
105    fNelectrons       = 0;
106    fStepID           =-1;
107
108    fNcharged         = 0;
109    fEloss            = 0;
110    fNcherenkov       = 0;
111
112    fParentTrack      = NULL;
113
114    SafeDelete(fHistoEnergy);
115    SafeDelete(fHistoLateral);
116    SafeDelete(fHistoEneAng);
117    SafeDelete(fHistoRadPhiEle);
118    SafeDelete(fHistoRadDTimeEle);
119    SafeDelete(fHistoRadPhiEloss);
120    SafeDelete(fHistoAngCher);
121}
122
123//______________________________________________________________________________
124const TH1F* ShowerStep::GetHistoEnergy() {
125    //
126    // This method returns pointer to the histogram with energy distribution of
127    // electrons at the given step filled by the Shower Generator. If Shower
128    // Generator does not provide this information the histogram is filled
129    // automatically according to configured analitical distribution.
130    //
131    if (fHistoEnergy) 
132        return  fHistoEnergy;
133    else {
134        // no parent, no histo
135        if ( !fParentTrack ) return 0;
136        // check if the histo for the given age interval already exists. In
137        // this case use it otherwise create the new one
138        Float_t age = (fAgei + fAgef)/2;
139        Int_t Age_index = -1;
140        Int_t fNumEnergyHistos = fParentTrack->GetNumEnergyHistos();
141        for (Int_t i(0); i<fNumEnergyHistos; i++) {
142            if (i*2./fNumEnergyHistos < age && age <= (i+1)*2./fNumEnergyHistos) Age_index = i;
143        }
144        // define histo unique name
145        char name[100], title[100];
146        sprintf(name,"ShowerStep.HistoEnergy.AgeIndex%d",Age_index);
147        sprintf(title,"Energy distribution");
148
149        Int_t Nx(1000);
150        Double_t Xbins[Nx+1]; 
151        Double_t MaxElectrons(1.e4);
152        Xbins[0] = 0.; 
153        for (Int_t i(1);i<=Nx;i++) {
154            if (i<=100)
155                Xbins[i] = i*0.01;                          // step size of 10 keV for [0,1] MeV
156            else if (i>=100 && i<= 500)
157                Xbins[i] = Xbins[i-1] + 99./400.;             // step size of 0.25 MeV for [1,100] MeV
158            else   
159                Xbins[i] = Xbins[i-1] + 1.;                 // step size of 1 MeV for [100, infty] MeV
160        }
161        fHistoEnergy = new TH1F(name,title,Nx,&Xbins[0]);
162        fHistoEnergy->SetDirectory(0);
163        TF12 EnergyDistribution("EnergyDistr",fParentTrack->GetEnergyDistribution("default"),
164                (fAgei + fAgef)/2,"X");
165
166        if (fNelectrons<MaxElectrons) {
167            fHistoEnergy->FillRandom("EnergyDistr",(Long64_t)fNelectrons); 
168            if (fNelectrons) fHistoEnergy->Scale(1/fNelectrons);
169        }
170        else {
171            fHistoEnergy->FillRandom("EnergyDistr",(Long64_t)MaxElectrons);
172            if (MaxElectrons) fHistoEnergy->Scale(1/MaxElectrons); 
173        } 
174        return fHistoEnergy;
175    }
176}
177
178//______________________________________________________________________________
179const TH1F* ShowerStep::GetHistoEnergy() const {
180    //
181    // const version of GetHistoEnergy
182    //
183
184    return ((ShowerStep*)this)->GetHistoEnergy();
185}
186
187//______________________________________________________________________________
188const TH1F* ShowerStep::GetHistoLateral() {
189    //
190    // This method returns pointer to the histogram with lateral distribution
191    // of electrons at the given step filled by the Shower Generator. If Shower
192    // Generator does not provide this information the histogram is filled
193    // automatically according to configured analitical distribution.
194    //
195    if (fHistoLateral) 
196        return  fHistoLateral;
197    else {
198        // no parent, no histo
199        if ( !fParentTrack ) return 0;
200        // check if the histo for the given age interval already exists. In
201        // this case use it otherwise create the new one
202        Float_t age = (fAgei + fAgef)/2;
203        Int_t Age_index = -1;
204        Int_t fNumLateralHistos = fParentTrack->GetNumLateralHistos();
205        for (Int_t i(0); i<fNumLateralHistos; i++) {
206            if (i*2./fNumLateralHistos < age && age <= (i+1)*2./fNumLateralHistos) Age_index = i;
207        }
208        // define histo unique name
209        char name[100], title[100];
210        sprintf(name,"ShowerStep.HistoLateral.AgeIndex%d",Age_index);
211        sprintf(title,"Lateral distribution");
212       
213        Int_t Nx(1000);
214        Double_t MaxElectrons(1.e4);
215        fHistoLateral = new TH1F(name,title,Nx,0,100);
216        fHistoLateral->SetDirectory(0);
217        TF12 LateralDistribution("LateralDistr",fParentTrack->GetLateralDistribution("default"),(fAgei + fAgef)/2,"X");
218
219        if (fNelectrons<MaxElectrons) {
220            fHistoLateral->FillRandom("LateralDistr",(Long64_t)fNelectrons);
221            if (fNelectrons) fHistoLateral->Scale(1/fNelectrons);
222        }
223        else {
224            fHistoLateral->FillRandom("LateralDistr",(Long64_t)MaxElectrons);
225            if (MaxElectrons) fHistoLateral->Scale(1/MaxElectrons); 
226        } 
227        return fHistoLateral;
228    }           
229}
230
231//______________________________________________________________________________
232const TH1F* ShowerStep::GetHistoLateral() const {
233    //
234    // const version of GetHistoLateral
235    //
236
237    return ((ShowerStep*)this)->GetHistoLateral();
238}
239
240//______________________________________________________________________________
241Bool_t ShowerStep::Dump() const {
242    //
243    // dump information about this step
244    //
245    Msg(EsafMsg::Info) << Form("BEGIN      Dumping step N %d",fStepID) << MsgDispatch;
246    Msg(EsafMsg::Info) << Form("POSITION B : X %f Y %f Z %f",fXYZi.X()/km,fXYZi.Y()/km,fXYZi.Z()/km) << MsgDispatch;
247    Msg(EsafMsg::Info) << Form("POSITION E : X %f Y %f Z %f",fXYZf.X()/km,fXYZf.Y()/km,fXYZf.Z()/km) << MsgDispatch;
248    Msg(EsafMsg::Info) << Form("TIME     B : %f",fTimei/microsecond) << MsgDispatch;
249    Msg(EsafMsg::Info) << Form("TIME     E : %f",fTimef/microsecond) << MsgDispatch;
250    Msg(EsafMsg::Info) << Form("AGE        : %f",fAgei) << MsgDispatch;
251    Msg(EsafMsg::Info) << Form("NELECTRONS : %f",fNelectrons) << MsgDispatch;
252    Msg(EsafMsg::Info) << Form("END        Dumping step N %d\n",fStepID) << MsgDispatch;
253    return kTRUE;
254
255}
Note: See TracBrowser for help on using the repository browser.