source: JEM-EUSO/esaf_cc_at_lal/packages/simulation/lightsources/src/ShowerLightSource.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: 21.4 KB
Line 
1// ESAF : Euso Simulation and Analysis Framework
2// $Id: ShowerLightSource.cc 2927 2011-06-16 18:02:41Z mabl $
3// Alessandro Thea created Nov, 24 2003
4
5#include "ShowerLightSource.hh"
6#include "Config.hh"
7
8#include "ListPhotonsInAtmosphere.hh"
9#include "FluoCalculator.hh"
10#include "CrkCalculator.hh"
11#include "LightSourceFactory.hh"
12
13#include "Atmosphere.hh"
14
15#include "EAtmosphere.hh"
16#include "EAtmosphereBunchAdder.hh"
17#include "utils.hh"
18#include <TH1F.h>
19
20#include <vector> // test
21
22// NB : Light generation only for shower part where " Ne > 0.001*Nemax "
23
24using namespace sou;
25
26ClassImp(ShowerLightSource)
27
28//____________________________________________________________________________________________
29  ShowerLightSource::ShowerLightSource() : LightSource("SHOWER"), EsafMsgSource(),
30                                           fFluocalcul(0),fCrkcalcul(0), fShowerNemax(0) {
31    //
32    // ctor
33    //
34    fPh_in_atmo = new ListPhotonsInAtmosphere;
35    if(!fPh_in_atmo) Msg(EsafMsg::Panic) << "Pb for memory allocation of fPh_in_atmo"<<MsgDispatch;
36
37    // set wl bounds for light generation
38    ConfigFileParser* pConfig = Config::Get()->GetCF("LightSource","ShowerLightSource");
39    fLambdaMin = pConfig->GetNum("ShowerLightSource.fLambdaMin")*nm;
40    fLambdaMax = pConfig->GetNum("ShowerLightSource.fLambdaMax")*nm;
41    // check if lowtran FORTRAN CODE setting matches
42    pConfig = Config::Get()->GetCF("Atmosphere","LowtranAtmosphere");
43    Double_t low_wlmin = pConfig->GetNum("LowtranAtmosphere.Linf")*nm;
44    Double_t low_wlmax = pConfig->GetNum("LowtranAtmosphere.Lsup")*nm;
45    if(fLambdaMin < low_wlmin || fLambdaMax > low_wlmax) {
46        Msg(EsafMsg::Warning)<< "<ShowerLightSource ctor> if Lowtran FORTRAN CODE is used : ";
47        Msg(EsafMsg::Warning)<< "Lambda bounds SHOULD match light generation ones" << MsgDispatch;
48    }
49   
50    Configure();
51    Msg(EsafMsg::Info) << "ShowerLightSource built. Will produce light within ["<<fLambdaMin/nm<<","<<fLambdaMax/nm<<"] nm" << MsgDispatch;
52
53}
54
55//_____________________________________________________________________________________________
56ShowerLightSource::~ShowerLightSource() {
57    //
58    // dtor
59    //
60    SafeDelete(fPh_in_atmo);
61    SafeDelete(fFluocalcul);
62    SafeDelete(fCrkcalcul);
63}
64
65//______________________________________________________________________________________________
66void ShowerLightSource::Configure() {
67    //
68    // configure fluorescence and cerenkov calculators
69    //
70
71    // Get fluorescence calculator
72    string fluoname = Conf()->GetStr("ShowerLightSource.FluoCalculator");
73    fFluocalcul = LightSourceFactory::Get()->GetFluoCalculator( fluoname );
74    Msg(EsafMsg::Info) << "Fluo calculator name " <<fFluocalcul->GetName()<< MsgDispatch;
75
76    // Get Cerenkov calculator
77    string crkname = Conf()->GetStr("ShowerLightSource.CrkCalculator");
78    fCrkcalcul = LightSourceFactory::Get()->GetCrkCalculator( crkname );
79    fCrkcalcul->SetWlRange(fLambdaMin,fLambdaMax);
80    Msg(EsafMsg::Info) << "Cerenkov calculator name " <<fCrkcalcul->GetName()<< MsgDispatch;
81
82    // Get Energy Distribution type and name
83    fEnergyDistributionType = Conf()->GetStr("ShowerLightSource.EnergyDistributionType");
84    fEnergyDistributionName = Conf()->GetStr("ShowerLightSource.EnergyDistributionName");
85    if(fEnergyDistributionType == "single") {
86        if(!(fEnergyDistributionName == "'80MeV'" || fEnergyDistributionName == "'24MeV'"))
87            Msg(EsafMsg::Panic) << "<Configure> Wrong type of electrons SINGLE energy = "<<fEnergyDistributionName<< MsgDispatch;
88    }
89    else if(fEnergyDistributionType == "parametrized") {
90        if(!(fEnergyDistributionName == "hillas" || fEnergyDistributionName == "giller" || fEnergyDistributionName == "nerling"))
91            Msg(EsafMsg::Panic) << "<Configure> Wrong type of PARAMETRIZED electrons energy distribution = "<<fEnergyDistributionName<< MsgDispatch;
92    }
93
94    // Get Lateral Distribution Name
95    fLateralDistributionName = Conf()->GetStr("ShowerLightSource.LateralDistributionName");
96
97    // Get Angular Distribution Name
98    fAngularDistributionName = Conf()->GetStr("ShowerLightSource.AngularDistributionName");
99   
100    // to know whether electrons angular deviation used for Yield calculation
101    string usedev = Conf()->GetStr("ShowerLightSource.fUseAngDev");
102    if(usedev == "yes") {
103        if(fAngularDistributionName == "NULL") {
104            Msg(EsafMsg::Warning) <<"<Configure()> fUseAngDev cannot be YES if no angular distribution -> has been set to NO "<<MsgDispatch;
105            fUseAngDev = false;
106        }
107        else fUseAngDev = true;
108    }
109    else if(usedev == "no") fUseAngDev = false;
110    else Msg(EsafMsg::Panic) <<"<Configure()> Wrong argument fot fUseAngDev -> "<<usedev<<MsgDispatch;
111
112}
113
114//________________________________________________________________________________________
115void ShowerLightSource::Reset() {
116    //
117    // reset internal list of photons
118    //
119    if(fPh_in_atmo) fPh_in_atmo->Reset();
120    if(fFluocalcul) fFluocalcul->Reset();
121    if(fCrkcalcul) fCrkcalcul->Reset();
122    fShowerNemax = 0;
123}
124
125//______________________________________________________________________________
126Bool_t ShowerLightSource::ClearMemory() {
127    //
128    // release all the mem hold in the buffers
129    //
130
131    Reset();
132    return fPh_in_atmo->ClearMemory();
133}
134
135//_________________________________________________________________________________________
136PhotonsInAtmosphere* ShowerLightSource::Get( const PhysicsData* data ) {
137    //
138    // generate photons in atmosphere from shower
139    //
140
141    Reset();
142
143    // get track
144    if ( data->Type() != "shower" )
145        Msg(EsafMsg::Panic) << "Wrong PhysicsData in ShowerLightSource. Must be shower"<<MsgDispatch;
146    ShowerTrack *track = (ShowerTrack*)data;
147
148    UInt_t nbshowstep = track->GetNumStep();
149    if(!nbshowstep) return (PhotonsInAtmosphere*)0;
150
151    // energy spectrum of charged particles
152    const char* ED_char = fEnergyDistributionName.c_str();
153    TString EDname;
154    EDname.Append(ED_char);
155#ifdef DEBUG
156    Msg(EsafMsg::Debug)<< "Energy Distribution Name "<< EDname << MsgDispatch;
157#endif
158    TF2* EnergyDistribution = NULL;
159    if(fEnergyDistributionType == "parametrized")
160        EnergyDistribution = track->GetEnergyDistribution(EDname);
161
162
163    // electron angular distribution
164    const char* AD_char = fAngularDistributionName.c_str();
165    TString ADname;
166    ADname.Append(AD_char);
167#ifdef DEBUG
168    Msg(EsafMsg::Debug)<< "Angular Distribution Name "<< ADname << MsgDispatch;
169#endif
170    TF2* AngularDistribution = NULL;
171    if(!((fAngularDistributionName == "NULL")||(fAngularDistributionName == "histos")) )
172        AngularDistribution = track->GetAngularDistribution(ADname);
173
174
175    // electron lateral distribution
176    const char* LD_char = fLateralDistributionName.c_str();
177    TString LDname;
178    LDname.Append(LD_char);
179#ifdef DEBUG
180    Msg(EsafMsg::Debug)<< "Lateral Distribution Name "<< LDname << MsgDispatch;
181#endif
182    TF2* LateralDistribution = NULL;
183    if(!((fLateralDistributionName == "NULL")||(fLateralDistributionName == "histos")) )
184        LateralDistribution = track->GetLateralDistribution(LDname);
185
186    // if BunchRadiativeTransfer is gonna to be used :
187    // initialization of the photons track -- assess the mean depth of ShowerSteps (assumed cst)
188    ConfigFileParser* pConfig = Config::Get()->GetCF("LightToEuso","StandardLightToEuso");
189    string RT = pConfig->GetStr("StandardLightToEuso.fRadiativeTransfer");
190    string generator_type = pConfig->GetStr("StandardLightToEuso.fGenerator");
191    UInt_t stepjump(0);
192    Double_t depthstep(0.);
193    if(RT == "bunch") {
194        pConfig = Config::Get()->GetCF("RadiativeTransfer","BunchRadiativeTransfer");
195        depthstep = pConfig->GetNum("BunchRadiativeTransfer.DepthStep")*g/cm2;
196        Double_t totalgram = Atmosphere::Get()->Grammage(track->FirstPos(),track->LastPos());
197        Double_t meanshowdepthstep =  totalgram / nbshowstep;
198        if(meanshowdepthstep > 0) stepjump = (UInt_t) floor(depthstep/meanshowdepthstep);
199        else Msg(EsafMsg::Warning)<<"PB : Total grammage travelled by showertrack = " <<totalgram*cm2/g<<" g/cm2" << MsgDispatch;
200        if(stepjump == 0) stepjump = 1;
201        fPh_in_atmo->ClearTrack();
202        if(!(nbshowstep % stepjump)) fPh_in_atmo->SetNbTrackSteps(nbshowstep/stepjump + 1);
203        else fPh_in_atmo->SetNbTrackSteps(nbshowstep/stepjump + 2);
204    }
205
206    // Estimate the total number of electrons
207    Float_t nTotal(0),nTracked(0);
208    for (UInt_t i=0;i<nbshowstep;i++) {
209        const ShowerStep& s = (*track)[i];
210#ifdef DEBUG
211        EarthVector tesst(s.GetXYZi());
212        EarthVector tesst2(s.GetXYZi() - (*track)[0].GetXYZi());
213        Msg(EsafMsg::Debug)<<"POSITION STEP IN UNISIM = " << tesst << MsgDispatch;
214        Msg(EsafMsg::Debug)<<"POSITION STEP IN UNISIM = " << tesst2.Mag()/km << MsgDispatch;
215        Msg(EsafMsg::Debug)<<"STEP ALTITUDE IN UNISIM = " << tesst.Zv() << MsgDispatch;
216        Msg(EsafMsg::Debug)<<"STEP length IN UNISIM = " << (s.GetXYZf() - s.GetXYZi()).Mag()/km << MsgDispatch;
217#endif
218        if (fShowerNemax < s.GetNelectrons()) fShowerNemax = s.GetNelectrons();
219        if(RT == "bunch")
220            if(!(i % stepjump)) fPh_in_atmo->FillTrack(s.GetXYZi());
221        nTotal += s.GetNelectrons();
222    }
223    if(RT == "bunch") fPh_in_atmo->FillTrack((*track)[nbshowstep-1].GetXYZf());
224
225    // in case UnisimGenerator has created ShowerTrack, fTrack is supplemented until its impact on ground or TOA
226    if(generator_type == "Unisimfile" || RT == "bunch") {
227        Int_t invert_status = 0;
228        EarthVector lastpos((*track)[nbshowstep-1].GetXYZf());
229        if(lastpos.IsUnderSeaLevel()) Msg(EsafMsg::Warning)<<"Unisim track seems to end "<<lastpos.Zv()/km<<"km under sea level "<< MsgDispatch;
230        else if(lastpos.Zv() > 10*m) {
231            EarthVector trackDir(lastpos - (*track)[0].GetXYZi());
232            trackDir = trackDir.Unit();
233            EarthVector nextpos(lastpos);
234            Int_t count_debug = 0;
235            // loop until ground or TOA reached
236            while(!(invert_status == 2 || invert_status == 3)) {
237                if(count_debug++ > 10000000) {
238                    Msg(EsafMsg::Warning)<<"<Get, build tracklight for Unisim event> Infinite loop broken by hand"<< MsgDispatch;
239                    break;
240                }
241                invert_status = Atmosphere::Get()->InvertGrammage(lastpos,trackDir,depthstep,nextpos);
242                if(invert_status == -1 || invert_status == 1) {
243                    Msg(EsafMsg::Warning)<<"<Get, build tracklight for Unisim event> WRONG invert status = "<<invert_status <<"SHOULD NOT occur"<< MsgDispatch;
244                    break;
245                }
246                fPh_in_atmo->FillTrack(nextpos,true);
247                lastpos = nextpos;
248            }
249        }
250        Msg(EsafMsg::Info)<<"NB of light track steps = " << fPh_in_atmo->GetNbTrackSteps() << MsgDispatch;
251        if(invert_status == 3 && !track->GetHitGround())
252            Msg(EsafMsg::Warning)<<"<Get, build tracklight for Unisim event> InvertGrammage tells impact on ground BUT NOT Unisim"<< MsgDispatch;
253    }
254
255#ifdef DEBUG
256    /*
257       EarthVector tesst((*track)[nbshowstep-1].GetXYZf());
258       EarthVector tesst2((*track)[nbshowstep-1].GetXYZf() - (*track)[0].GetXYZi());
259       Msg(EsafMsg::Debug)<<"END POSITION of last STEP IN UNISIM = " << tesst << MsgDispatch;
260       Msg(EsafMsg::Debug)<<"END POSITION of last STEP IN UNISIM = " << tesst2.Mag()/km << MsgDispatch;
261       Msg(EsafMsg::Debug)<<"END POSITION of last STEP, ALTITUDE IN UNISIM = " << tesst.Zv() << MsgDispatch;
262       for (UInt_t i=0;i<fPh_in_atmo->GetNbTrackSteps();i++) {
263       EarthVector tesst(fPh_in_atmo->GetTrackStep(i));
264       EarthVector tesst2(fPh_in_atmo->GetTrackStep(0) - tesst);
265       Msg(EsafMsg::Debug)<<"POSITION STEP IN PhotonsList = " << tesst << MsgDispatch;
266       Msg(EsafMsg::Debug)<<"LAST POSITION STEP IN PhotonsList = " << tesst2.Mag()/km << MsgDispatch;
267       Msg(EsafMsg::Debug)<<"LAST STEP ALTITUDE IN PhotonsList = " << tesst.Zv() << MsgDispatch;
268       if(i>0) Msg(EsafMsg::Debug)<<"STEP (before) length IN PhotonsList = " << (fPh_in_atmo->GetTrackStep(i-1) - fPh_in_atmo->GetTrackStep(i)).Mag()/km << MsgDispatch;
269       }
270     */
271#endif
272
273
274    // generate light
275    Int_t progress = 0;
276
277    for (UInt_t i=0;i<nbshowstep;i++) {
278        const ShowerStep& s = (*track)[i];
279        // generate fluorescence bunch
280        BunchOfPhotons* bfluo = MakeFluoStep(s,EnergyDistribution,LateralDistribution,AngularDistribution);
281        fPh_in_atmo->Add(bfluo);
282        // generate cerenkov bunch
283        BunchOfPhotons* bcer = MakeCerenkovStep(s,EnergyDistribution,LateralDistribution,AngularDistribution);
284        fPh_in_atmo->Add(bcer);
285
286        nTracked += s.GetNelectrons();
287
288        if (100.*(i+1)/nbshowstep >= progress) {
289            Msg(EsafMsg::Info).SetProgress(progress);
290            Msg(EsafMsg::Info) << "Processing:" << MsgCount;
291            progress+=10;
292        }
293    }
294
295    return fPh_in_atmo;
296
297}
298//_________________________________________________________________________________________
299MCTruth* ShowerLightSource::Truth() {
300    return NULL;
301}
302
303//_________________________________________________________________________________________
304BunchOfPhotons* ShowerLightSource::MakeFluoStep(const ShowerStep& step,TF2* EnergyDistribution,
305                                                TF2* LateralDistribution, TF2* AngularDistribution) {
306    //
307    // generate a fluorescence bunch of photons for this shower step
308    //
309
310    //mean position in space
311    EarthVector posi = step.GetXYZi();
312    EarthVector posf = step.GetXYZf();
313#ifdef DEBUG
314    //Msg(EsafMsg::Debug)<<" X "<< posi.X()/km <<" Y "<<posi.Y()/km<<" Z "<<posi.Z()/km<<MsgDispatch;
315#endif
316    EarthVector pos=0.5*(posi+posf);
317
318
319    // fluorescence yield for the energy spectrum of charged particles
320    EsafSpectrum spectrum(357*nm);  //initialisation
321    Double_t age = 0.5*(step.GetAgef()+step.GetAgei());
322    Double_t TotalYield = 0;
323
324    if( step.GetNelectrons() > 0.001*fShowerNemax) {
325
326        if ( fEnergyDistributionType == "histos" ) {       
327            const TH1F* ehisto = step.GetHistoEnergy();
328            TotalYield = fFluocalcul->GetFluoYieldHisto(pos.Zv(),ehisto,&spectrum); 
329    #ifdef DEBUG
330           // Msg(EsafMsg::Debug)<<"integration over TH1F - Yield Fluo " << TotalYield*m << MsgDispatch;
331    #endif
332        }
333        else if ( EnergyDistribution && fEnergyDistributionType == "parametrized") {
334            TF12 EnergyDistributionStep("EDS_name",EnergyDistribution,age,"X");
335            TotalYield = fFluocalcul->GetFluoYield(pos.Zv(),&EnergyDistributionStep,&spectrum);
336    #ifdef DEBUG
337           // Msg(EsafMsg::Debug)<<"integration over TF12 - Yield Fluo " << TotalYield*m << MsgDispatch;
338    #endif
339        }
340        else if(fEnergyDistributionType == "single") {
341            if(fEnergyDistributionName == "'80MeV'") TotalYield = fFluocalcul->GetFluoYield(pos.Zv(),80*MeV,&spectrum);
342            else if(fEnergyDistributionName == "'24MeV'") TotalYield = fFluocalcul->GetFluoYield(pos.Zv(),24*MeV,&spectrum);
343    #ifdef DEBUG
344           // Msg(EsafMsg::Debug)<<" E= 80 MeV - Yield Fluo " << TotalYield*m << MsgDispatch;
345    #endif
346        }
347        else Msg(EsafMsg::Panic)<<"<MakeFluoStep> Wrong type of EnergyDistribution = "<<fEnergyDistributionType<< MsgDispatch;
348    }
349    else TotalYield = 0.;
350
351    // projection of angular distribution along X (fixed value of Y)
352    // energy threshold fixed to minimum value to get the larger spectrum
353    // presently REMOVE because of electrons Energy cut issue
354    //TF12* AngularDistributionStep = 0;
355    //if(AngularDistribution) AngularDistributionStep = new TF12("ADS_name",AngularDistribution,0.5,"X");
356
357    // total weight of the bunch of photons
358    Double_t dl = (step.GetXYZf()-step.GetXYZi()).Mag(); 
359    //if(fUseAngDev) dl /= cos(AngularDistributionStep->Mean(AngularDistributionStep->GetXmin(),AngularDistributionStep->GetXmax()));
360    Double_t nbPhotonsInBunch = step.GetNelectrons() * TotalYield * dl;
361#ifdef DEBUG
362   // Msg(EsafMsg::Debug)<< "nelec "<<step.GetNelectrons()<<" nbpib fluo "<<nbPhotonsInBunch << " dl= " << dl/m << MsgDispatch;
363#endif
364    // Create a new bunch of Photons with its ParentBunch
365    BunchOfPhotons* b = new BunchOfPhotons(nbPhotonsInBunch,TotalYield,posi,posf,step.GetTimei(),
366                                           step.GetTimef(),spectrum,step.GetParentTrack()->GetDirVers(),Fluo,age);
367
368   
369    // Set the lateral distribution to the bunch if it exists
370    if ( fLateralDistributionName == "histos" ) {
371        const TH1F* lhisto = step.GetHistoLateral();
372        if (lhisto) b->SetParentLateral(*lhisto);
373    }
374    else if ( LateralDistribution ) {
375        if ( LateralDistribution->GetNpar() == 1 ) {
376            // gives moliere radius as parameter in meter
377            // Get Atmosphere for Density
378            const Atmosphere* atmo = Atmosphere::Get();
379            ///// WARNING : if standard NKG used, Rm should be calculated at 2X0 back in shower dvpt
380            /////           but if NKG for hadrons used, this correction must not be made (as it comes from a fit at the present depth)
381            Double_t Rm = 9.6 * gram/cm2 / atmo->Air_Density( pos.Zv() ) / m; // in meters
382            LateralDistribution->SetParameters(Rm,0);
383       }
384        TF12 LateralDistributionStep("LDS_name",LateralDistribution,age,"X");
385        b->SetParentLateral(LateralDistributionStep);
386    }
387   
388    //SafeDelete(AngularDistributionStep);
389   
390    return b;
391}
392
393
394//_________________________________________________________________________________________
395BunchOfPhotons* ShowerLightSource::MakeCerenkovStep(const ShowerStep& step, TF2* EnergyDistribution,
396                                                    TF2* LateralDistribution, TF2* AngularDistribution) {
397    //
398    // Generate a Cerenkov Bunch Of Photons for this shower step
399    //
400
401    //mean position in space
402    EarthVector posi = step.GetXYZi();
403    EarthVector posf = step.GetXYZf();
404    EarthVector pos=0.5*(posi+posf);
405
406    // Cerenkov yield for the energy spectrum of charged particles
407    Double_t age = 0.5*(step.GetAgef()+step.GetAgei());
408    Double_t TotalYield = 0;
409   
410    if( step.GetNelectrons() > 0.001*fShowerNemax) {
411       if ( fEnergyDistributionType == "histos" ) {
412           const TH1F* ehisto = step.GetHistoEnergy();
413           TotalYield = fCrkcalcul->GetCrkYield(pos.Zv(),ehisto);
414   #ifdef DEBUG
415        //   Msg(EsafMsg::Debug)<<"integration over TH1F - Yield Cerenkov " << TotalYield*m << MsgDispatch;
416   #endif
417       }
418       else if ( EnergyDistribution && fEnergyDistributionType == "parametrized" ) {
419           TF12 EnergyDistributionStep("EDS_name",EnergyDistribution,age,"X");
420           TotalYield = fCrkcalcul->GetCrkYield(pos.Zv(),&EnergyDistributionStep);
421   #ifdef DEBUG
422         //  Msg(EsafMsg::Debug)<<"integration over TF12 - Yield Cerenkov " << TotalYield*m << MsgDispatch;
423   #endif
424       }
425       else if(fEnergyDistributionType == "single") {
426           if(fEnergyDistributionName == "'80MeV'") TotalYield = fCrkcalcul->GetCrkYield(pos.Zv(),80*MeV);
427           else if(fEnergyDistributionName == "'24MeV'") TotalYield = fCrkcalcul->GetCrkYield(pos.Zv(),24*MeV,false); // false to remove energy threshold for ckov prod
428   #ifdef DEBUG
429        //   Msg(EsafMsg::Debug)<<"E=80MeV - Yield Cerenkov " << TotalYield*m << MsgDispatch;
430   #endif
431       }
432       else Msg(EsafMsg::Panic)<<"<MakeCerenkovStep> Wrong type of EnergyDistribution = "<<fEnergyDistributionType<< MsgDispatch;
433    }
434    else TotalYield = 0;
435
436    // Cerenkov wavelenght spectrum
437    EsafSpectrum *spectrum = fCrkcalcul->GetCrkSpectrum();
438
439    // projection  angular distribution ofalong X (fixed value of Y)
440    TF12* AngularDistributionStep = 0;
441    if(AngularDistribution) AngularDistributionStep = new TF12("ADS_name",AngularDistribution,fCrkcalcul->GetEnergyThreshold(pos.Zv())/MeV,"X");
442
443    // weight of the bunch of photons
444    Double_t dl = (step.GetXYZf() - step.GetXYZi()).Mag();
445    if(fUseAngDev) dl /= cos(AngularDistributionStep->Mean(AngularDistributionStep->GetXmin(), AngularDistributionStep->GetXmax()));
446    Double_t nelec = step.GetNelectrons();
447    Double_t nbPhotonsInBunch = nelec*TotalYield*dl;
448#ifdef DEBUG
449   // Msg(EsafMsg::Debug)<< "nelec "<<nelec<<" nbpib cer "<<nbPhotonsInBunch << " dl= " << dl/m << MsgDispatch;
450#endif
451
452    // Create a new bunch of Photons with its ParentBunch
453    BunchOfPhotons* b = new BunchOfPhotons(nbPhotonsInBunch,TotalYield,posi,posf,step.GetTimei(),
454                        step.GetTimef(),*spectrum,step.GetParentTrack()->GetDirVers(),Cerenkov,age);
455
456
457    // Set the lateral distribution to the parent bunch
458    if ( fLateralDistributionName == "histos" ) {
459        const TH1F* lhisto = step.GetHistoLateral();
460        if (lhisto) b->SetParentLateral(*lhisto);
461        Msg(EsafMsg::Warning)<<"Histogram for Angular distribution not handled here, it SHOULD"<< MsgDispatch; 
462    }
463    else if ( LateralDistribution ) {
464        if ( LateralDistribution->GetNpar() == 1 ) {
465            // gives moliere radius as parameter in meter
466            // Get Atmosphere for Density
467            const Atmosphere* atmo = Atmosphere::Get();
468            ///// WARNING : if standard NKG used, Rm should be calculated at 2X0 back in shower dvpt
469            /////           but if NKG for hadrons used, this correction must not be made (as it comes from a fit at the present depth)
470            Double_t Rm = 9.6 * gram/cm2 / atmo->Air_Density( pos.Zv() ) / m; // in meters
471            LateralDistribution->SetParameters(Rm,0);
472        }
473        TF12 LateralDistributionStep("LDS_name",LateralDistribution,age,"X");
474        b->SetParentLateral(LateralDistributionStep);
475    }
476   
477    // Set the angular distribution to the parent bunch
478    if ( AngularDistribution ) b->SetParentAngular(*AngularDistributionStep);
479   
480    SafeDelete(spectrum);
481    SafeDelete(AngularDistributionStep);
482   
483    return b;
484}
485
486
Note: See TracBrowser for help on using the repository browser.