source: JEM-EUSO/esaf_lal/tags/v1_r0/esaf/packages/simulation/lightsources/src/TestLightSource.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: 39.3 KB
Line 
1// ESAF : Euso Simulation and Analysis Framework
2// $Id: TestLightSource.cc 2681 2006-05-03 10:58:33Z moreggia $
3// Anne Stutz created Dec,  2 2003
4//
5// This class brings test photons in atmosphere to Euso
6// The type of events supported are the following:
7//     SPOT  light from a single point at the same time
8//     RANDSPOT light from a random point at the same time
9//     TRACK light making a track in atmosphere
10//     RANDTRACK
11//
12// The parameters needed are the following :
13//     TestLightSource.Option         see above can be SPOT or TRACK
14//     TestLightSource.Photons      number of photons generated per event (if<1000 SinglePhoton else BunchOfPhotons)
15//     TestLightSource.Theta1       lower zenith angle                                 
16//     TestLightSource.Theta2       higher zenith angle                                 
17//     TestLightSource.Phi1         lower azimuth angle
18//     TestLightSource.Phi2         higher azimuth angle
19//     TestLightSource.FirstPointX  X first point source in km
20//     TestLightSource.FirstPointY  Y first point source in km
21//     TestLightSource.FirstPointZ  Z first point source in km
22//     TestLightSource.LongExtension = 0          longitudinal extension of the bunch in g/cm2 can be 0
23//     TestLightSource.LatDistribution = NKG2      lateral distribution of the bunch can be NKG NULL
24//     TestLightSource.AngDistribution = baltru   angular distribution of the bunch can be baltru NULL
25//     TestLightSource.SpectrumType   can be FLUO,MONO,CERENKOV
26//     TestLightSource.Lambda         Wavelenght if MONO
27//     TestLightSource.DirectionType  can be EUSO,ISOTROPIC,MONO
28
29#include "TestLightSource.hh"
30#include "ListPhotonsInAtmosphere.hh"
31#include "EsafRandom.hh"
32#include "Config.hh"
33#include "EsafSpectrum.hh"
34#include "FluoCalculator.hh"
35#include "LightSourceFactory.hh"
36#include "TFormula.h"
37#include "TMath.h"
38#include "EConst.hh"
39#include "EarthVector.hh"
40#include "RadiativeFactory.hh"
41#include "Ground.hh"
42#include "Atmosphere.hh"
43
44ClassImp(TestLightSource)
45
46using namespace TMath;
47using namespace sou;
48using namespace EConst;
49
50 /* dNe/dr derived from 'standard' NGK formula for Lateral distribution. r is the distance to shower axis
51 Rm=moliere radius is a parameter . The integration of NGK2 over r from 0. to infinity is normalized to 1 . The integral
52 from radius r1 to radius r2 give the fraction of total number of electrons which lie inside such interval.*/
53Double_t TestLightSourceNKG2(Double_t *x, Double_t *par) {
54  Double_t R = x[0], age = x[1], Rm=par[0];
55  Double_t e1=2.0;
56  Double_t e2=4.5;
57  Double_t D=R/Rm;
58  return Gamma(e2-age)/Gamma(age)/Gamma(e2-2.*age)*Power(D,age-(e1-1.0))
59        *Power((1.0+D),age-e2)/Rm;
60}
61
62Double_t TestLightSourceBaltru(Double_t *x, Double_t *par) {
63    //
64    // dNe/dtheta(theta,Et) from Baltrusaitis et al. J.Phys.G:Nucl. Phys. 13 (1987)
65    // where theta is the angle between the electrons and the shower axis and Et
66    // (MeV) is the energy thr for electrons considered. The integration of
67    // dNe/dtheta(theta,Et) over dtheta from 0 to Pi() is normalized to 1.  The
68    // integral from theta1 to theta2 give the fraction of total number of
69    // electrons which lie inside such angular interval (distribution in phi is
70    // supposed uniform).
71    //
72   
73    Double_t theta = x[0], Et = x[1];
74    Double_t a = 0.85;
75    Double_t b = 0.66;
76    Double_t theta0 = a*Power(Et,-b);
77    Double_t pigreco = Pi();
78   
79    return Exp(-theta/theta0) / theta0 / (1.- Exp(-pigreco/theta0));
80}
81
82//_________________________________________________________________________________________________
83TestLightSource::TestLightSource() : LightSource("TEST"), EsafMsgSource(), fFluocalcul(0) {
84    //
85    // ctor
86    //
87
88    Msg(EsafMsg::Info) << "Start Building TestLightSource" << MsgDispatch;
89
90    fPh_in_atmo = new ListPhotonsInAtmosphere;
91    if(!fPh_in_atmo) Msg(EsafMsg::Panic) << "Pb for memory allocation of fPh_in_atmo" << MsgDispatch;
92 
93    Configure();
94    Msg( EsafMsg::Info) << "TestLightSource built"<<MsgDispatch;
95}
96
97//_________________________________________________________________________________________________
98TestLightSource::~TestLightSource() {
99    //
100    // dtor
101    //
102   
103    SafeDelete(fPh_in_atmo);
104    SafeDelete(fFluocalcul);
105    SafeDelete(fLateralDistribution);
106    SafeDelete(fAngularDistribution);
107}
108
109//________________________________________________________________________
110void TestLightSource::Configure() {
111    //
112    // Configure TestLightSource
113    //
114
115    // Get detector position
116    ConfigFileParser* pConf = Config::Get()->GetCF("General","Euso");
117    fEUSO.SetZ(pConf->GetNum("Euso.fAltitude")*km);
118   
119    // Get fluorescence calculator if required
120    string fluoname = Conf()->GetStr("TestLightSource.SpectrumType");
121    if(fluoname == "FLUO_kakimoto") {
122        fluoname = "kakimoto";
123        fFluocalcul = LightSourceFactory::Get()->GetFluoCalculator( fluoname );
124        Msg(EsafMsg::Info) << "Fluo calculator name " <<fFluocalcul->GetName()<< MsgDispatch;
125    }
126
127    // Lateral and Angular distributions
128    fLateralDistribution = NULL;
129    fAngularDistribution = NULL;
130    string LateralDistributionName = Conf()->GetStr("TestLightSource.LatDistribution");
131    string AngularDistributionName = Conf()->GetStr("TestLightSource.AngDistribution");
132 
133    if (LateralDistributionName == "NKG2" ) 
134        fLateralDistribution = new TF2("LatDist",TestLightSourceNKG2,0.001,5000,0,2,1);
135    else if(LateralDistributionName == "NULL" )
136        fLateralDistribution = NULL;
137    else Msg(EsafMsg::Panic) << "Lateral distribution name is not valid : "<<LateralDistributionName<< MsgDispatch;
138   
139    if (AngularDistributionName == "baltru" )
140        fAngularDistribution = new TF2("AngDist",TestLightSourceBaltru,0.,Pi(),.5,1000.);
141    else if(AngularDistributionName == "NULL" )
142        fAngularDistribution = NULL;
143    else Msg(EsafMsg::Panic) << "Angular distribution name is not valid : "<<AngularDistributionName<< MsgDispatch;
144}
145
146//_________________________________________________________________________________________________
147void TestLightSource::Reset() {
148    //
149    // reset internal list of photons
150    //
151
152    if(fPh_in_atmo) fPh_in_atmo->Reset();
153    if(fFluocalcul) fFluocalcul->Reset();
154}
155
156//______________________________________________________________________________
157Bool_t TestLightSource::ClearMemory() {
158    //
159    // release all the mem hold in the buffers
160    //
161   
162    Reset();
163    return fPh_in_atmo->ClearMemory();
164   
165}
166
167//_________________________________________________________________________________________________
168PhotonsInAtmosphere *TestLightSource::Get( const PhysicsData *dummy) {
169    //
170    // return the list of photons in atmosphere
171    //
172
173    Reset();
174    string option = Conf()->GetStr("TestLightSource.Option");
175    Double_t nbPhotons = Conf()->GetNum("TestLightSource.Nbph");
176    Double_t fFirstPointX = (Conf()->GetNum("TestLightSource.FirstPointX"))*km;
177    Double_t fFirstPointY = (Conf()->GetNum("TestLightSource.FirstPointY"))*km;
178    Double_t fFirstPointZ = (Conf()->GetNum("TestLightSource.FirstPointZ"))*km;
179
180    string impactMode = Conf()->GetStr("TestLightSource.ImpactMode");
181    Double_t ImpactXmin = Conf()->GetNum("TestLightSource.ImpactXmin")*km;
182    Double_t ImpactXmax = Conf()->GetNum("TestLightSource.ImpactXmax")*km;
183    Double_t ImpactYmin = Conf()->GetNum("TestLightSource.ImpactYmin")*km;
184    Double_t ImpactYmax = Conf()->GetNum("TestLightSource.ImpactYmax")*km;
185
186
187    TRandom* rndm = EsafRandom::Get();
188
189    EarthVector Posi;
190    if (fFirstPointZ > 100*km ) {
191        Msg(EsafMsg::Warning) << "In TestLightSource FirstPointZ set at 100 km" << MsgDispatch; 
192        fFirstPointZ = 100*km;
193    }
194    Posi.SetXYZ(fFirstPointX,fFirstPointY,fFirstPointZ);
195
196    // spot in fixed position
197    if ( option == "SPOT" ) MakeSpot(Posi,nbPhotons);
198
199    // spot mimicking shower maximum, depends on theta
200    else if ( option == "HmaxSPOT") {
201        Double_t x = rndm->Rndm()*(ImpactXmax - ImpactXmin) + ImpactXmin;
202        Double_t y = rndm->Rndm()*(ImpactYmax - ImpactYmin) + ImpactYmin;
203        Double_t z = - (x*x + y*y)/(2*EarthRadius());   // approched value
204        MakeHmaxSpot(EarthVector(x,y,z),nbPhotons);
205    }
206
207    // spot mimicking shower maximum, depends on theta
208    else if ( option == "ShowerLike") {
209        Double_t x = rndm->Rndm()*(ImpactXmax - ImpactXmin) + ImpactXmin;
210        Double_t y = rndm->Rndm()*(ImpactYmax - ImpactYmin) + ImpactYmin;
211        Double_t z = - (x*x + y*y)/(2*EarthRadius());   // approched value
212        MakeShowerLike(EarthVector(x,y,z),nbPhotons);
213    }
214
215    // For bgnd sutdies : photons generated at (0,0,50km) isotropically
216    else if ( option == "Bgnd") {
217        MakeBgnd(nbPhotons);
218    }
219
220    // track in fixed position
221    else if ( option == "TRACK" ) {   
222        // random direction between theta1 and theta2, and between phi1 and phi2
223        Double_t theta1 = DegToRad() *Conf()->GetNum("TestLightSource.Theta1");
224        Double_t theta2 = DegToRad() *Conf()->GetNum("TestLightSource.Theta2");
225        Double_t phi1 = DegToRad() *Conf()->GetNum("TestLightSource.Phi1");
226        Double_t phi2 = DegToRad() *Conf()->GetNum("TestLightSource.Phi2");
227
228        Double_t thmax = theta2;
229        Double_t thmin = theta1;
230        if (theta1>theta2) {
231            thmax = thmin;
232            thmin = theta2;
233        }
234        Double_t theta = thmin + rndm->Rndm()*(thmax-thmin);
235
236        Double_t phmax = phi2;
237        Double_t phmin = phi1;
238        if (phi1>phi2) {
239            phmax = phmin;
240            phmin = phi2;
241        }
242        Double_t phi = phmin + rndm->Rndm()*(phmax-phmin);
243
244        MakeTrack(theta,phi,Posi,nbPhotons); 
245    }
246
247    else {
248        Msg(EsafMsg::Panic)<<"Invalid parameter TestLightSource.Option = " << option <<MsgDispatch;
249        return (PhotonsInAtmosphere*)0;
250    }
251
252    return fPh_in_atmo;
253}
254
255//___________________________________________________________________________________________
256PhotonsInAtmosphere *TestLightSource::MakeSpot(const EarthVector& Posi, Double_t nbPhotons) {
257    //
258    // Produce light in a single spot
259    //
260
261    Msg(EsafMsg::Info)<<"TestLightSource SPOT CALLED at (km) "<<Posi.X()/km<<" "
262        <<Posi.Y()/km<<" "<<Posi.Z()/km<<MsgDispatch;
263
264    string directiontype = Conf()->GetStr("TestLightSource.DirectionType");
265    string photontype = Conf()->GetStr("TestLightSource.PhotonType");
266    string photonDescription = Conf()->GetStr("TestLightSource.Descrip");
267    Double_t theta1 = DegToRad() *Conf()->GetNum("TestLightSource.Theta1");
268    Double_t theta2 = DegToRad() *Conf()->GetNum("TestLightSource.Theta2");
269    Double_t phi1 = DegToRad() *Conf()->GetNum("TestLightSource.Phi1");
270    Double_t phi2 = DegToRad() *Conf()->GetNum("TestLightSource.Phi2");
271    PhotonType phtype = Fluo;
272    if(photontype == "cerenkov") phtype = Cerenkov;
273
274    TRandom* rndm = EsafRandom::Get();
275
276    string thetaMode = Conf()->GetStr("TestLightSource.ThetaMode");
277    Double_t thmax = theta2;
278    Double_t thmin = theta1;
279    Double_t phmax = phi2;
280    Double_t phmin = phi1;
281    if (theta1>theta2) {
282        thmax = thmin;
283        thmin = theta2;
284    }
285    if (phi1>phi2) {
286        phmax = phmin;
287        phmin = phi2;
288    }
289    Double_t theta = thmin + rndm->Rndm()*(thmax-thmin);
290    Double_t phi   = phmin + rndm->Rndm()*(phmax-phmin);
291    EarthVector Omega(1);
292    if(thetaMode == "local") {
293        if(theta > PiOver2()) Msg(EsafMsg::Panic)<<"theta > 90° not possible in mode \"local\" " <<MsgDispatch;
294        EarthVector v1 = Posi.Unit();
295        EarthVector vrot;
296        Omega.SetXYZ(1,0,0);
297        EarthVector Uz(0,0,1);
298        vrot = Uz.Cross(v1);
299        if(vrot.Mag() > 0) {
300            Omega.Rotate(v1.Theta(),vrot);
301            Omega.Rotate(phi,v1);
302            vrot = v1.Cross(Omega);
303            Omega.Rotate(Pi()/2+theta,vrot);
304        }
305        else Omega.SetMagThetaPhi(1.,Pi() - theta,phi + Pi());
306        theta = Omega.Theta();
307        phi   = Omega.Phi();
308    }
309    else if(thetaMode == "mes") {
310        theta = Pi() - theta;
311        phi += Pi();
312    }
313    else Msg(EsafMsg::Panic)<<"Invalid parameter TestLightSource.ThetaMode = " << thetaMode <<MsgDispatch;
314
315
316    if (directiontype == "UNIQUE") {
317        Msg(EsafMsg::Info)<<"                with theta between "<<thmin/deg<<" and "<<thmax/deg<<MsgDispatch;
318        Msg(EsafMsg::Info)<<"                with phi between   "<<phmin/deg<<" and "<<phmax/deg<<MsgDispatch;
319    }
320    else if (directiontype == "ISOTROPIC") {
321        Msg(EsafMsg::Info)<<"                with isotropic direction" << MsgDispatch;
322    }
323    else if (directiontype == "EUSO") {
324        Msg(EsafMsg::Info)<<"                direct to EUSO" << MsgDispatch;
325    }
326    else Msg(EsafMsg::Panic)<<"Invalid parameter TestLightSource.DirectionType = " << directiontype <<MsgDispatch;
327
328    // wavelenght spectrum
329    EsafSpectrum spectrum(357*nm);
330    string spectrumtype = Conf()->GetStr("TestLightSource.SpectrumType");
331    Double_t TotalYield = MakeSpectrum(Posi,&spectrum,spectrumtype,photontype);
332
333    //generation of SinglePhoton or BunchOfPhotons
334    if(photonDescription == "single") {
335        // only single photons direct to Euso
336        if (directiontype == "ISOTROPIC" || directiontype == "UNIQUE") 
337            Msg(EsafMsg::Info)<<"In TestLightSource SinglePhoton are only emitted direct to Euso"<< MsgDispatch;
338        for (int i=0; i<(int)nbPhotons; i++) {         
339            Double_t wl = spectrum.GetLambda();
340            EarthVector dir(0,0,0);
341            dir = EUSO() - Posi;
342            SinglePhoton *p = new SinglePhoton(0,wl,Posi,dir,phtype,Direct);
343            fPh_in_atmo->Add(p);
344        }
345    } 
346
347    // one bunch of nbPhotons
348    else if(photonDescription == "bunch"){ 
349        EarthVector dir(1);
350        dir.SetTheta(theta);
351        dir.SetPhi(phi);
352        EarthVector Posf = GetLongitudinalExtension(Posi,dir);
353        Double_t tf = (Posf-Posi).Mag()/Clight();
354        BunchOfPhotons* b = new BunchOfPhotons(nbPhotons,TotalYield,Posi,Posf,0,tf,spectrum,dir,phtype);
355        if( fLateralDistribution ) b->SetParentLateral(GetLateralDistribution(Posi));
356        if( fAngularDistribution && (phtype==Cerenkov) ) b->SetParentAngular(GetAngularDistribution(Posi.Zv()));
357
358        fPh_in_atmo->Add(b);
359
360        // track needed for ckov AND fluo (for last transfer until detector)
361        BuildLightTrack(Posi,dir);
362        //DELETE
363        // for CERENKOV RadiativeTransfer handling, when AlongTrack_CSPropagator is used
364        //if(phtype == Cerenkov) BuildLightTrack(Posi,dir);  //DELETE
365    }
366    else Msg(EsafMsg::Panic)<<"<MakeSpot> Invalid photonDescription parameter = " << photonDescription<<MsgDispatch;
367
368    return fPh_in_atmo; 
369}
370
371//___________________________________________________________________________________________
372PhotonsInAtmosphere *TestLightSource::MakeBgnd(Double_t nbPhotons) {
373    //
374    // produce background photons at (0,0,98km), sin(2theta) distribution, in the MES, between [300nm,MaxWl] (MaxWl configured)
375    //
376
377    Msg(EsafMsg::Info)<<"TestLightSource Background mode called" <<MsgDispatch;
378    Msg(EsafMsg::Warning)<<"Be carefull about TOA altitude set in PureMCRT (should be > point source altitude)" <<MsgDispatch;
379
380    PhotonType phtype = Cerenkov;
381    TRandom* rndm = EsafRandom::Get();
382    Double_t Zstart = (Conf()->GetNum("TestLightSource.FirstPointZ"))*km;
383
384    // init
385    Double_t costh1 = 1.;  // cos(2*0)
386    Double_t costh2 = -1.; // cos(2*halfpi)
387    Double_t theta = 0.;
388    Double_t phi   = 0.;
389    EarthVector Posi(0,0,Zstart);
390    EarthVector dir(1);
391   
392    // wavelenght spectrum
393    EsafSpectrum spectrum(357*nm);
394    string spectrumtype = Conf()->GetStr("TestLightSource.SpectrumType");
395    if(spectrumtype != "FLAT") Msg(EsafMsg::Panic)<<"Bgnd option must use FLAT spectrum"<<MsgDispatch;
396    MakeSpectrum(Posi,&spectrum,spectrumtype,"dummy");
397
398   
399    // create photons
400    for (Int_t i=0; i<(Int_t)nbPhotons; i++) {         
401        Double_t wl = spectrum.GetLambda();
402        theta = 0.5*ACos(costh1 - rndm->Rndm()*(costh1 - costh2));
403        theta = Pi() - theta;
404        phi = rndm->Rndm()*TwoPi();
405        dir.SetMagThetaPhi(1.,theta,phi);
406        SinglePhoton *p = new SinglePhoton(0,wl,Posi,dir,phtype,Direct);
407        fPh_in_atmo->Add(p);
408    }
409
410    return fPh_in_atmo; 
411}
412
413//___________________________________________________________________________________________
414PhotonsInAtmosphere *TestLightSource::MakeHmaxSpot(const EarthVector& impact, Double_t nbPhotons) {
415    //
416    // Produce light in a single spot
417    //
418
419    Msg(EsafMsg::Info)<<"TestLightSource HmaxSPOT with impact at (km) "<<impact.X()/km<<" "
420        <<impact.Y()/km<<" "<<impact.Z()/km<<MsgDispatch;
421
422    string directiontype = Conf()->GetStr("TestLightSource.DirectionType");
423    string photontype = Conf()->GetStr("TestLightSource.PhotonType");
424    string photonDescription = Conf()->GetStr("TestLightSource.Descrip");
425    Double_t theta1 = DegToRad() *Conf()->GetNum("TestLightSource.Theta1");
426    Double_t theta2 = DegToRad() *Conf()->GetNum("TestLightSource.Theta2");
427    Double_t phi1 = DegToRad() *Conf()->GetNum("TestLightSource.Phi1");
428    Double_t phi2 = DegToRad() *Conf()->GetNum("TestLightSource.Phi2");
429    PhotonType phtype = Fluo;
430    if(photontype == "cerenkov") phtype = Cerenkov;
431
432    TRandom* rndm = EsafRandom::Get();
433
434    // bunch direction (for cerenkov bunches only)
435    string thetaMode = Conf()->GetStr("TestLightSource.ThetaMode");
436    Double_t thmax = theta2;
437    Double_t thmin = theta1;
438    Double_t phmax = phi2;
439    Double_t phmin = phi1;
440    if (theta1>theta2) {
441        thmax = thmin;
442        thmin = theta2;
443    }
444    if (phi1>phi2) {
445        phmax = phmin;
446        phmin = phi2;
447    }
448    if(thmax > 80*DegToRad()) {
449        thmax = 80*DegToRad();
450        Msg(EsafMsg::Warning)<<"With HmaxSPOT, theta limited to 80 deg -> theta=80 applied"<<MsgDispatch;
451    }
452    Double_t theta_true = thmin + rndm->Rndm()*(thmax-thmin);
453    Double_t phi = phmin + rndm->Rndm()*(phmax-phmin);
454    Double_t theta = theta_true;
455
456    EarthVector Omega(1);
457    // if theta given is local -> theta in MES is calculated
458    if(thetaMode == "local") {
459        EarthVector v1 = impact.Unit();
460        EarthVector vrot;
461        Omega.SetXYZ(1,0,0);
462        EarthVector Uz(0,0,1);
463        if(vrot.Mag() > 0) {
464            Omega.Rotate(v1.Theta(),vrot);
465            Omega.Rotate(phi,v1);
466            vrot = v1.Cross(Omega);
467            Omega.Rotate(Pi()/2+theta,vrot);
468        }
469        else Omega.SetMagThetaPhi(1.,Pi() - theta,phi + Pi());
470        theta = Omega.Theta();
471        phi   = Omega.Phi();
472    }
473    else if(thetaMode == "mes") Msg(EsafMsg::Panic)<<"mes ThetaMode NOT possible with HmaxSPOT option"<<MsgDispatch;
474    else Msg(EsafMsg::Panic)<<"Invalid parameter TestLightSource.ThetaMode = " << thetaMode <<MsgDispatch;
475
476    // get Hmax from theta_true
477    EarthVector Posi = GetHmaxPos(impact,Omega,theta_true);
478    Msg(EsafMsg::Info)<<"TestLightSource HmaxSPOT with INITIAL pos at (km) "<<Posi.X()/km<<" " <<Posi.Y()/km<<" "<<Posi.Z()/km<<MsgDispatch;
479
480    if (directiontype == "UNIQUE") {
481        Msg(EsafMsg::Info)<<"                with theta between "<<thmin/deg<<" and "<<thmax/deg<<MsgDispatch;
482        Msg(EsafMsg::Info)<<"                with phi between   "<<phmin/deg<<" and "<<phmax/deg<<MsgDispatch;
483    }
484    else if (directiontype == "ISOTROPIC") {
485        Msg(EsafMsg::Info)<<"                with isotropic direction" << MsgDispatch;
486    }
487    else if (directiontype == "EUSO") {
488        Msg(EsafMsg::Info)<<"                direct to EUSO" << MsgDispatch;
489    }
490    else Msg(EsafMsg::Panic)<<"Invalid parameter TestLightSource.DirectionType = " << directiontype <<MsgDispatch;
491
492    // wavelenght spectrum
493    EsafSpectrum spectrum(357*nm);
494    string spectrumtype = Conf()->GetStr("TestLightSource.SpectrumType");
495    Double_t TotalYield = MakeSpectrum(Posi,&spectrum,spectrumtype,photontype);
496
497    // only single photons direct to Euso
498    if(photonDescription == "single") {
499        if (directiontype == "ISOTROPIC" || directiontype == "UNIQUE") 
500            Msg(EsafMsg::Info)<<"In TestLightSource SinglePhoton are only emitted direct to Euso"<< MsgDispatch;
501
502        for (int i=0; i<(int)nbPhotons; i++) {         
503            Double_t wl = spectrum.GetLambda();
504            EarthVector dir(0,0,0);
505            dir = EUSO() - Posi;
506            SinglePhoton *p = new SinglePhoton(0,wl,Posi,dir,phtype,Direct);
507            fPh_in_atmo->Add(p);
508        }
509    } 
510
511    // one bunch of nbPhotons
512    else if(photonDescription == "bunch"){ 
513        EarthVector dir(1);
514        dir.SetTheta(Omega.Theta());
515        dir.SetPhi(Omega.Phi());
516        EarthVector Posf = GetLongitudinalExtension(Posi,dir);
517        Double_t tf = (Posf - Posi).Mag()/Clight();
518        BunchOfPhotons* b = new BunchOfPhotons(nbPhotons,TotalYield,Posi,Posf,0,tf,spectrum,dir,phtype);
519        if( fLateralDistribution ) b->SetParentLateral(GetLateralDistribution(Posi));
520        if( fAngularDistribution && (phtype==Cerenkov) ) b->SetParentAngular(GetAngularDistribution(Posi.Zv()));
521
522        fPh_in_atmo->Add(b);
523
524        // track needed for ckov AND fluo (for last transfer until detector)
525        BuildLightTrack(Posi,dir);
526        //DELETE
527        // for CERENKOV RadiativeTransfer handling, when AlongTrack_CSPropagator is used
528        //if(phtype == Cerenkov) BuildLightTrack(Posi,dir);  //DELETE
529    }
530    else Msg(EsafMsg::Panic)<<"<MakeHmaxSpot> Invalid photonDescription parameter = " << photonDescription<<MsgDispatch;
531
532    return fPh_in_atmo; 
533}
534
535//___________________________________________________________________________________________
536PhotonsInAtmosphere *TestLightSource::MakeShowerLike(const EarthVector& impact, Double_t nbPhotons) {
537    //
538    // Produce light mimiking a point-like shower : one bunch fluo, one bunch ckov
539    // So far, not usable in BunchRadiativeTransfer mode
540    //
541
542    Msg(EsafMsg::Info)<<"TestLightSource ShowerLike with impact at (km) "<<impact.X()/km<<" "
543        <<impact.Y()/km<<" "<<impact.Z()/km<<MsgDispatch;
544
545    string photonDescription = Conf()->GetStr("TestLightSource.Descrip");
546    if(photonDescription == "single") Msg(EsafMsg::Warning)<<"<MakeShowerLike> single photons generation not possible in this mode"<<MsgDispatch;
547    Double_t theta1 = DegToRad() *Conf()->GetNum("TestLightSource.Theta1");
548    Double_t theta2 = DegToRad() *Conf()->GetNum("TestLightSource.Theta2");
549    Double_t phi1 = DegToRad() *Conf()->GetNum("TestLightSource.Phi1");
550    Double_t phi2 = DegToRad() *Conf()->GetNum("TestLightSource.Phi2");
551
552    TRandom* rndm = EsafRandom::Get();
553
554    // bunch direction (for cerenkov bunches only)
555    string thetaMode = Conf()->GetStr("TestLightSource.ThetaMode");
556    Double_t thmax = theta2;
557    Double_t thmin = theta1;
558    Double_t phmax = phi2;
559    Double_t phmin = phi1;
560    if (theta1>theta2) {
561        thmax = thmin;
562        thmin = theta2;
563    }
564    if (phi1>phi2) {
565        phmax = phmin;
566        phmin = phi2;
567    }
568    if(thmax > 80*DegToRad()) {
569        thmax = 80*DegToRad();
570        Msg(EsafMsg::Warning)<<"With ShowerLike, theta limited to 80 deg -> theta=80 applied"<<MsgDispatch;
571    }
572    Double_t theta_true = thmin + rndm->Rndm()*(thmax-thmin);
573    Double_t phi = phmin + rndm->Rndm()*(phmax-phmin);
574    Double_t theta = theta_true;
575
576    EarthVector Omega(1);
577    // if theta given is local -> theta in MES is calculated
578    if(thetaMode == "local") {
579        EarthVector v1 = impact.Unit();
580        EarthVector vrot;
581        Omega.SetXYZ(1,0,0);
582        EarthVector Uz(0,0,1);
583        if(vrot.Mag() > 0) {
584            Omega.Rotate(v1.Theta(),vrot);
585            Omega.Rotate(phi,v1);
586            vrot = v1.Cross(Omega);
587            Omega.Rotate(Pi()/2+theta,vrot);
588        }
589        else Omega.SetMagThetaPhi(1.,Pi() - theta,phi + Pi());
590        theta = Omega.Theta();
591        phi   = Omega.Phi();
592    }
593    else if(thetaMode == "mes") Msg(EsafMsg::Panic)<<"mes ThetaMode NOT possible with ShowerLike option"<<MsgDispatch;
594    else Msg(EsafMsg::Panic)<<"Invalid parameter TestLightSource.ThetaMode = " << thetaMode <<MsgDispatch;
595
596    // get Hmax from theta_true
597    EarthVector Posi = GetHmaxPos(impact,Omega,theta_true);
598    Msg(EsafMsg::Info)<<"TestLightSource ShowerLike with INITIAL pos at (km) "<<Posi.X()/km<<" " <<Posi.Y()/km<<" "<<Posi.Z()/km<<MsgDispatch;
599    Msg(EsafMsg::Info)<<"TestLightSource ShowerLike with INITIAL (theta, phi) = ("<<Posi.Theta()*RadToDeg()<<", " <<Posi.Phi()*RadToDeg() <<") "<<MsgDispatch;
600
601    // wavelenght spectrum
602    EsafSpectrum spectrumfluo(357*nm);
603    EsafSpectrum spectrumckov(357*nm);
604    string spectrumtype = Conf()->GetStr("TestLightSource.SpectrumType");
605    Double_t TotalYieldfluo = MakeSpectrum(Posi,&spectrumfluo,spectrumtype,"fluo");
606    Double_t TotalYieldckov = MakeSpectrum(Posi,&spectrumckov,"LAMBDA2","cerenkov");
607
608    // one fluo bunch
609    EarthVector dir(1);
610    dir.SetTheta(Omega.Theta());
611    dir.SetPhi(Omega.Phi());
612    EarthVector Posf = GetLongitudinalExtension(Posi,dir);
613    Double_t tf = (Posf - Posi).Mag()/Clight();
614    BunchOfPhotons* bfluo = new BunchOfPhotons(nbPhotons,TotalYieldfluo,Posi,Posf,0,tf,spectrumfluo,dir,Fluo);
615    if( fLateralDistribution ) bfluo->SetParentLateral(GetLateralDistribution(Posi));
616
617    fPh_in_atmo->Add(bfluo);
618
619    // one ckov bunch
620    BunchOfPhotons* bckov = new BunchOfPhotons(nbPhotons,TotalYieldckov,Posi,Posf,0,tf,spectrumckov,dir,Cerenkov);
621    if( fLateralDistribution ) bckov->SetParentLateral(GetLateralDistribution(Posi));
622    if( fAngularDistribution ) bckov->SetParentAngular(GetAngularDistribution(Posi.Zv()));
623
624    fPh_in_atmo->Add(bckov);
625
626    // track needed for ckov AND fluo (for last transfer until detector)
627    //BuildLightTrack(Posi,dir);
628    //DELETE
629    // for CERENKOV RadiativeTransfer handling, when AlongTrack_CSPropagator is used
630    //if(phtype == Cerenkov) BuildLightTrack(Posi,dir);  //DELETE
631
632    return fPh_in_atmo; 
633}
634
635//_____________________________________________________________________________________________________
636PhotonsInAtmosphere *TestLightSource::MakeTrack(Double_t theta, Double_t phi, const EarthVector& Posi, Double_t nbPhotons) {
637    //
638    // Produce light along a track
639    //
640
641    Msg(EsafMsg::Info)<<"TestLightSource TRACK CALLED theta="<<theta<<" phi ="<<phi<<MsgDispatch;
642    Msg(EsafMsg::Info)<<"               Starting point at (km) "<<Posi.X()/km<<" "<<Posi.Y()/km<<" "<<Posi.Z()/km<<MsgDispatch;
643    string directiontype = Conf()->GetStr("TestLightSource.DirectionType");
644    string photontype = Conf()->GetStr("TestLightSource.PhotonType");
645    string photonDescription = Conf()->GetStr("TestLightSource.Descrip");
646    PhotonType phtype = Fluo;
647    if(photontype == "cerenkov") phtype = Cerenkov;
648
649    string thetaMode = Conf()->GetStr("TestLightSource.ThetaMode");
650    EarthVector Omega(1);
651    if(thetaMode == "local") {
652        EarthVector v1 = Posi.Unit();
653        EarthVector vrot;
654        Omega.SetXYZ(1,0,0);
655        EarthVector Uz(0,0,1);
656        if(vrot.Mag() > 0) {
657            Omega.Rotate(v1.Theta(),vrot);
658            Omega.Rotate(phi,v1);
659            vrot = v1.Cross(Omega);
660            Omega.Rotate(Pi()/2+theta,vrot);
661        }
662        else Omega.SetMagThetaPhi(1.,Pi() - theta,phi + Pi());
663        theta = Omega.Theta();
664        phi   = Omega.Phi();
665    }
666    else if(thetaMode == "mes") {
667        theta = Pi() - theta;
668        phi += Pi();
669    }
670    else Msg(EsafMsg::Panic)<<"Invalid parameter TestLightSource.ThetaMode = " << thetaMode <<MsgDispatch;   
671    EarthVector DirTrack(1);
672    DirTrack.SetTheta(theta);
673    DirTrack.SetPhi(phi);
674
675    EsafSpectrum spectrum(357*nm);
676    Double_t TotalYield;
677
678    //generation of SinglePhoton or BunchOfPhotons
679    // impact on ground
680    Ground* ground = RadiativeFactory::Get()->GetGround();
681    EarthVector impact = ground->GetImpact(Posi,DirTrack);
682    Msg(EsafMsg::Info)<<"               Impact on ground at "<<impact.X()/km<<" "<<impact.Y()/km<<" "<<impact.Z()/km<< MsgDispatch;
683    Double_t magmax = (Posi - (ground->GetImpact(Posi,DirTrack))).Mag();
684
685    // only single photons
686    if(photonDescription == "single") {
687        Msg(EsafMsg::Info)<<"                direct to EUSO" << MsgDispatch;
688        if (directiontype == "ISOTROPIC" || directiontype == "UNIQUE") 
689            Msg(EsafMsg::Info)<<"In TestLightSource SinglePhoton are only emitted direct to Euso"<< MsgDispatch;
690        for (int i=0; i<(int)nbPhotons; i++) {
691            Double_t mag = magmax*(i+1)/nbPhotons;
692            EarthVector pos = Posi + mag*DirTrack;
693            Double_t t = mag/Clight();
694
695            string spectrumtype = Conf()->GetStr("TestLightSource.SpectrumType");
696            TotalYield = MakeSpectrum(pos,&spectrum,spectrumtype,photontype);
697
698            Double_t wl = spectrum.GetLambda();
699            EarthVector dir = EUSO()-pos;
700            SinglePhoton *p = new SinglePhoton(t,wl,pos,dir,phtype,Direct); 
701
702            fPh_in_atmo->Add(p);
703        }
704    }
705    // 100 bunchs of photons
706    else if(photonDescription == "bunch"){ 
707        if (directiontype == "EUSO") Msg(EsafMsg::Info)<<"In TestLightSource BunchPhoton are not emitted directly to Euso"<< MsgDispatch;
708
709        Double_t nbPhotonsInBunch = nbPhotons/100.;
710        for (Float_t i=0; i<100; i++) {
711            Double_t mag = magmax*(i+1)/100;
712            EarthVector pos = Posi + mag*DirTrack;
713            Double_t t = mag/Clight();
714
715            string spectrumtype = Conf()->GetStr("TestLightSource.SpectrumType");
716            TotalYield = MakeSpectrum(pos,&spectrum,spectrumtype,photontype);
717
718            EarthVector posf = GetLongitudinalExtension(pos,DirTrack);
719            if ( posf == pos ) break;
720            Double_t tf = t + (posf-pos).Mag()/Clight();
721
722            BunchOfPhotons*  b = new BunchOfPhotons(nbPhotonsInBunch,TotalYield,pos,posf,t,tf,spectrum,DirTrack,phtype);
723
724            if( fLateralDistribution ) b->SetParentLateral(GetLateralDistribution(pos));
725            if( fAngularDistribution && (phtype==Cerenkov) ) b->SetParentAngular(GetAngularDistribution(Posi.Zv()));
726
727            fPh_in_atmo->Add(b);
728        }
729
730        // for CERENKOV RadiativeTransfer handling, when AlongTrack_CSPropagator is used
731        if(phtype == Cerenkov) BuildLightTrack(Posi,DirTrack);
732    }
733    else Msg(EsafMsg::Panic)<<"<MakeTrack> Invalid photonDescription parameter = " << photonDescription<<MsgDispatch;
734
735    return fPh_in_atmo; 
736}
737
738//_________________________________________________________________________________________________
739Double_t TestLightSource::MakeSpectrum(const EarthVector& pos,EsafSpectrum* spectrum, string spectrumtype, string photontype) {
740    //
741    // Calculation of the wavelenght spectrum
742    //
743
744    Double_t TotalYield=0;
745
746    if (spectrumtype == "MONO") {
747        Double_t lambda = (Conf()->GetNum("TestLightSource.Lambda"))*nm;
748        if ( spectrum!=0) spectrum -> Reset(lambda);
749        TotalYield = 1.;
750        fLambdaMin = lambda - 1*nm;
751        fLambdaMax = lambda + 1*nm;
752    }
753    else if (spectrumtype == "FLUO_kakimoto" && photontype == "fluo") {
754        Double_t energy=80.*MeV;
755        TotalYield = fFluocalcul->GetFluoYield(pos.Zv(),energy,spectrum);
756        fLambdaMin = spectrum->GetLambdaMin();
757        fLambdaMax = spectrum->GetLambdaMax();
758    }
759    else if (spectrumtype == "LAMBDA2" && photontype == "cerenkov") {
760         TFormula cerenkov("cerenkov","1 /(x*x)");
761         if ( spectrum !=0 ) spectrum -> Reset(&cerenkov,100,300*nm,450*nm);
762         TotalYield = 1.;
763         fLambdaMin = spectrum->GetLambdaMin();
764         fLambdaMax = spectrum->GetLambdaMax();
765    }
766    else if (spectrumtype == "FLAT") {
767         Double_t lambda = (Conf()->GetNum("TestLightSource.Lambda"))*nm;
768         TFormula flat("flat","1");
769         if ( spectrum !=0 ) spectrum -> Reset(&flat,100,300*nm,lambda);
770         TotalYield = 1.;
771         fLambdaMin = spectrum->GetLambdaMin();
772         fLambdaMax = spectrum->GetLambdaMax();
773    }
774    else Msg(EsafMsg::Panic)<<"Invalid parameter TestLightSource.SpectrumType or PhotonType = " << spectrumtype<<", "<<photontype <<MsgDispatch;
775
776    if(!spectrum) Msg(EsafMsg::Panic)<<"Pb of memory allocation in TestLightSource"<<MsgDispatch;
777
778    // check if lowtran FORTRAN CODE WLrange setting matches
779    ConfigFileParser* pConfig = Config::Get()->GetCF("Atmosphere","LowtranAtmosphere");
780    Double_t low_wlmin = pConfig->GetNum("LowtranAtmosphere.Linf")*nm;
781    Double_t low_wlmax = pConfig->GetNum("LowtranAtmosphere.Lsup")*nm;
782    if(fLambdaMin < low_wlmin || fLambdaMax > low_wlmax) {
783        Msg(EsafMsg::Warning)<< "<MakeSpectrum> if Lowtran FORTRAN CODE is used : ";
784        Msg(EsafMsg::Warning)<< "Lambda bounds SHOULD match light generation ones" << MsgDispatch;
785    }
786   
787    return TotalYield;
788}
789
790
791//_________________________________________________________________________________________________
792const TF12 TestLightSource::GetLateralDistribution(const EarthVector& pos) {
793    //
794    // return TF12 the lateral distribution at pos for age = 1
795    //   
796
797    TF12 LatDist;
798    Double_t age = 1.;
799
800    // if NKG in meter gives moliere radius as parameter in meter
801    if ( fLateralDistribution) {
802        if ( fLateralDistribution->GetNpar() == 1 ) { 
803           // Get Atmosphere for Density
804           const Atmosphere* atmo = Atmosphere::Get();
805           Double_t Rm = 9.6 * gram/cm2 / atmo->Air_Density( pos.Zv() ) / m; // in meters
806           Msg(EsafMsg::Debug)<<"rayon moliere =  " << Rm <<MsgDispatch;
807           fLateralDistribution->SetParameters(Rm,0);
808        }   
809        else Msg(EsafMsg::Panic)<<"Lateral Distribution should have 1 parameter" <<MsgDispatch;
810        LatDist = TF12("LDS_name",fLateralDistribution,age,"X"); 
811    }
812    else Msg(EsafMsg::Panic)<<"No Lateral Distribution" <<MsgDispatch;
813
814    return LatDist;
815}
816
817//_________________________________________________________________________________________________
818const TF12 TestLightSource::GetAngularDistribution(Double_t alt) {
819    //
820    // return Pointer to the angular distribution
821    //
822
823   TF12 AngDist;
824
825   Double_t EthCer = GetEnergyThreshold(alt);
826   if ( fAngularDistribution ) AngDist = TF12("ADS_name",fAngularDistribution,EthCer,"X");
827   else Msg(EsafMsg::Panic)<<"No Angular Distribution "<<MsgDispatch;
828
829   return AngDist;
830}
831
832//__________________________________________________________________________________________________
833EarthVector TestLightSource::GetLongitudinalExtension(const EarthVector& Pos,const EarthVector& Dir) {
834    //
835    // return last point of a bunch with first point Pos and mean direction Dir
836    //
837    Double_t LgExt = Conf()->GetNum("TestLightSource.LongExtension")*g/cm2;
838    EarthVector Posf = Pos;
839    /*
840    Double_t L,h(0),depth(0);
841
842    L = LgExt / Atmosphere::Get()->Air_Density(h)/100.;
843    Int_t cycle=0;
844    while(1) {
845        Posf += Dir*L;
846        if ( Posf.IsUnderSeaLevel() ) return Pos;
847        depth += Atmosphere::Get()->Air_Density(Posf.Zv())*L;
848        if ( depth >= LgExt ) break;
849        if ( cycle>100000 ) {
850          MsgForm( EsafMsg::Debug,"Next point not found : cycle > 100000 with Lgext %f",LgExt/g*cm2 );
851          return Pos;
852        }
853    }
854    */
855    Int_t status = Atmosphere::Get()->InvertGrammage(Pos,Dir,LgExt,Posf);
856    if(status != 0) {
857        Msg(EsafMsg::Warning) << "<GetLongitudinalExtension> : CANNOT FIND FINAL POSITION, InvertGrammage status = "<<status<<MsgDispatch;
858        return Pos;
859    }
860   
861    return Posf;
862}
863
864//__________________________________________________________________________________________________
865EarthVector TestLightSource::GetHmaxPos(const EarthVector& impact,const EarthVector& Direc,Double_t theta) const {
866    //
867    // get Hmax 3D-position -- for HmaxSPOT option only
868    // Here theta MUST be LOCAL zenith angle (at impact)
869    //
870   
871    EarthVector rtn(0,0,0);
872    EarthVector Dir(Direc);
873    if(Dir.Mag() != 1) Dir = Dir.Unit();
874    Double_t tolerance = 0.5*m;
875    if(theta > 80) Msg(EsafMsg::Warning) << "theta > 80deg NOT foreseen in <TestLightSource::GetHmaxPos>"<<MsgDispatch;
876   
877    // get Hmax value from LOCAL theta
878    Double_t Hmax = (1.89 - 7.59*Log(Cos(theta))) * km;
879   
880    // get 3D-position
881    EarthVector first = impact;
882    EarthVector last = impact - (30*km/Cos(Dir.Theta() - Pi()))*Dir;
883    EarthVector middle(1.);
884
885    while((last.Zv() - first.Zv()) > tolerance) {
886        middle = 0.5*(last+first);
887        if(fabs(middle.Zv() - Hmax) < tolerance) {rtn = middle; break;}
888        if(Hmax < middle.Zv()) last = middle;
889        else first = middle;
890    }
891    if(rtn.Mag() == 0) rtn = first;
892   
893#ifdef DEBUG
894cout<<" IN GETHMAXPOS :"<<endl;
895cout<<"impact position = "<<impact<<endl;
896cout<<"Hmax foreseen = "<<Hmax<<endl;
897cout<<"Hmax implemented = "<<rtn.Zv()<<endl;
898#endif
899   
900    return rtn;
901}
902
903//__________________________________________________________________________________________________
904void TestLightSource::BuildLightTrack(const EarthVector& posinit, const EarthVector& dir) {
905    //
906    // for CERENKOV RadiativeTransfer handling, when AlongTrack_CSPropagator is used
907    //
908   
909    // init
910    ConfigFileParser* pConfig = Config::Get()->GetCF("RadiativeTransfer","BunchRadiativeTransfer");
911    Double_t depthstep = pConfig->GetNum("BunchRadiativeTransfer.DepthStep")*g/cm2;
912    fPh_in_atmo->ClearTrack();
913    const Atmosphere* atmo = Atmosphere::Get();
914    EarthVector presentpos(posinit);
915    EarthVector nextpos(posinit);
916    EarthVector posi(posinit);
917    if(posi.IsUnderSeaLevel()) {
918        Msg(EsafMsg::Warning) << "<BuildLightTrack> Starting pos is under sea level -> set to Nadir"<<MsgDispatch;
919        posi.SetMag(0.);
920    }
921    Int_t status = -100;  // for Atmosphere::InvertGrammage() status
922   
923    // look at impact and total depth
924    EarthVector impact = atmo->ImpactASL(posi,dir);
925    if(impact.Z() == HUGE) impact = atmo->ImpactAtTOA(posi,dir);
926    if(impact.Z() == HUGE) Msg(EsafMsg::Warning) << "<BuildLightTrack> : No GROUND impact nor TOA impact -> SHOULD NOT"<<MsgDispatch;
927    Double_t FinalDepth = atmo->Grammage(posi,impact);
928   
929    // get nb of TrackLight steps
930    UInt_t nb = UInt_t(floor(FinalDepth/depthstep)) + 2;  // +1 for posi,  +1 for tuning last step exit
931    if(nb > 1000) {
932        Msg(EsafMsg::Warning) << "<BuildLightTrack> : Tracklength is huge             = " << (impact - posi).Mag()/km<<" km"<<MsgDispatch;
933        Msg(EsafMsg::Warning) << "<BuildLightTrack> : Corresponding depth is huge too = " << FinalDepth*cm2/g <<" g/cm2"<<MsgDispatch;
934        Msg(EsafMsg::Warning) << "<BuildLightTrack> : Thus nb of TRACKLIGHT steps set to 1000"<<MsgDispatch;
935        nb = 1000;
936        depthstep = FinalDepth / (nb - 2); // -1 for posi,  -1 for tuning last step exit
937    }
938    fPh_in_atmo->SetNbTrackSteps(nb);
939   
940    // fill the track until last step entry (last step exit done by hand below
941    for(UInt_t i=0; i<nb-1; i++) {
942        presentpos = nextpos;
943        fPh_in_atmo->FillTrack(presentpos);
944        if(i == (nb-2)) break;  // to avoid line below for last step of the loop
945        status = atmo->InvertGrammage(presentpos,dir,depthstep,nextpos);
946        if(status != 0) Msg(EsafMsg::Warning) << "<BuildLightTrack> InvertGrammage() status is problematic -> status = " << status <<MsgDispatch;
947    }
948   
949    // fine tune last step
950    fPh_in_atmo->FillTrack(impact);
951
952
953#ifdef DEBUG
954cout<<"impact position = "<<impact<<endl;
955    FinalDepth = atmo->Grammage(presentpos,impact);
956cout<<"finaldepth w.r.t. depthstep = "<<FinalDepth/depthstep<<endl;
957    if(FinalDepth/depthstep > 1) Msg(EsafMsg::Debug) << "<BuildLightTrack> : Last TRACKLIGHT step tuning failed, it counts as "<<FinalDepth/depthstep<<" times a normal step"<<MsgDispatch;   
958#endif
959
960}
961
962//____________________________________________________________________________________________
963Double_t TestLightSource::GetEnergyThreshold(Double_t SC_alt) const {
964    //
965    // get the energy threshold for cerenkov emission
966    //
967
968    // Get Atmosphere Index
969    const Atmosphere* atmo = Atmosphere::Get();
970    Double_t delta = atmo->Index_Minus1(SC_alt);
971    // When Index is not implemented in the Atmosphere use the Corsika formula
972    if ( delta == 0) delta = 0.000283 * atmo->Air_Density( SC_alt )/atmo->Air_Density(0);
973   
974    return ElectronMassC2()/sqrt(2*delta);
975}
976
977
Note: See TracBrowser for help on using the repository browser.