source: JEM-EUSO/esaf_cc_at_lal/packages/simulation/lightsources/src/NaganoFluoCalculator.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: 12.6 KB
Line 
1// ESAF : Euso Simulation and Analysis Framework
2// contact person : Naoto SAKAKI <n-sakaki@riken.jp>
3// $Id: NaganoFluoCalculator.cc 2776 2006-11-23 16:53:19Z moreggia $
4
5#include <stdexcept>
6#include "NaganoFluoCalculator.hh"
7#include "EsafSpectrum.hh"
8#include "Atmosphere.hh"
9#include "EConst.hh"
10
11#include <TMath.h>
12#include <TProfile.h>
13#include "TGraph.h"
14
15using namespace TMath;
16using namespace sou;
17using namespace EConst;
18
19ClassImp(NaganoFluoCalculator)
20
21// number of fluorescence bins
22const Int_t NaganoFluoCalculator::fNumWavelengths = 15;
23const Double_t hPa = 100.*pascal;
24
25// 391nm and 428nm are 1N bands and others are 2P bands
26const Double_t NaganoFluoCalculator::fWavelength[]  = {
27  316*nm, 329*nm, 337*nm, 354*nm, 358*nm,
28  376*nm, 381*nm, 391*nm, 394*nm, 400*nm,
29  406*nm, 414*nm, 420*nm, 427*nm, 428*nm};
30
31const Double_t NaganoFluoCalculator::fPhi0[]  = {
32  4.80e-4, 0.880e-4, 10.01e-4, 0.769e-4, 7.82e-4,
33  1.20e-4, 2.46e-4,   9.60e-4, 0.42e-4,  0.847e-4,
34  1.49e-4, 0.327e-4,  0.86e-4, 0.069e-4, 4.57e-4};
35
36const Double_t NaganoFluoCalculator::fRefPressure[] = {
37    23.0*hPa, 40.2*hPa, 19.2*hPa, 30.6*hPa, 18.1*hPa,
38    34.1*hPa, 19.4*hPa,  5.02*hPa,24.2*hPa, 24.2*hPa,
39    12.3*hPa, 19.3*hPa,  7.3*hPa, 72.*hPa, 3.86*hPa};
40
41
42//______________________________________________________________________________
43NaganoFluoCalculator::NaganoFluoCalculator() {
44    //
45    // Constructor
46    //
47    fName = "Nagano";
48    Msg(EsafMsg::Info) << "NaganoFluoCalculator built" << MsgDispatch;
49
50    TString formula = "gaus(0)";
51    for (Int_t i(1); i<fNumWavelengths; i++)
52        formula += Form("+gaus(%d)",i*3);
53     
54    fNaganoFluo = new TFormula("fluo - nagano", formula ); 
55}
56
57//______________________________________________________________________________
58NaganoFluoCalculator::~NaganoFluoCalculator() {
59    //
60    // Destructor
61    //
62
63    SafeDelete(fNaganoFluo);
64}
65
66
67//______________________________________________________________________________
68Double_t NaganoFluoCalculator::GetFluoYield(const Double_t NF_alt,
69                   const Double_t NF_energy, EsafSpectrum* FluoSpectrum) const {
70    //
71    // Get the wavelength spectrum of the fluorescence emission as a function
72    // of electron energy
73    //
74    // Input:
75    //     Double_t NF_alt         the altitude where the electron is moving
76    //     Double_t NF_energy      the energy of the electron
77    // Output:
78    //     fTotalYield             total fluorescence yield between 300 and 430 nm
79    //                             per electron per unit length
80    //     EsafSpectrum* FluoSpectrum fluorescence light spectrum per electron
81    //                             per unit length
82    //
83
84    // [1] Nagano et al. Astroparticle Physics, 20 (2003) 293.(astro-ph/0303193)
85    // [2] Nagano et al. astro-ph/0406474 (2004)
86    // measurements of yield for 0.85MeV electrons
87
88#ifdef DEBUG
89    Msg(EsafMsg::Debug) << "NF "<<NF_alt/km << " " <<NF_energy/MeV << MsgDispatch;
90#endif
91    // Get Atmosphere
92    const Atmosphere* atmo = Atmosphere::Get();
93
94    Double_t density     = atmo->Air_Density(NF_alt);
95    Double_t temperature = atmo->Temperature(NF_alt);
96    if(temperature<=0)
97        FatalError("Invalid Temperature in NaganoFluoCalculator");
98
99    const Double_t dedx_ref = GetdEdX( 0.85*MeV );
100
101    // density and temperature dependance from Nagano et al.
102    const static Double_t A[] = {
103        20.5*m2/kg,   3.91*m2/kg, 45.6*m2/kg, 3.68*m2/kg, 37.8*m2/kg,
104        6.07*m2/kg, 12.7*m2/kg,  50.8*m2/kg, 2.25*m2/kg,  4.58*m2/kg,
105        8.18*m2/kg,  1.83*m2/kg,  4.9*m2/kg, 0.40*m2/kg, 26.5*m2/kg
106    };
107    const static Double_t B[] = {
108        2.14*m3/kg, 1.22*m3/kg, 2.56*m3/kg, 1.60*m3/kg,  2.72*m3/kg,
109        1.44*m3/kg, 2.53*m3/kg, 9.80*m3/kg, 2.03*m3/kg,  2.03*m3/kg,
110        3.99*m3/kg, 2.55*m3/kg, 6.8 *m3/kg, 0.68*m3/kg, 12.7*m3/kg
111    };
112
113    Double_t fTotalYield = 0.0;
114    Double_t Yield[fNumWavelengths];
115
116    for (Int_t i=0; i<fNumWavelengths; i++) {
117        Yield[i] = GetdEdX( NF_energy ) / dedx_ref
118            *density*A[i]/(1+density*B[i]*sqrt(temperature));
119        fTotalYield += Yield[i];
120
121        fNaganoFluo->SetParameter(i*3,Yield[i]);
122        fNaganoFluo->SetParameter(i*3+1,fWavelength[i]);
123        fNaganoFluo->SetParameter(i*3+2,0.05*nm);
124    }
125    if ( FluoSpectrum !=0 )
126        FluoSpectrum->Reset(fNaganoFluo, 130, 300*nm, 430*nm);
127
128#ifdef DEBUG
129    Msg(EsafMsg::Debug) << "NF TotalYield/m "<<fTotalYield*m  << MsgDispatch;
130#endif
131
132    return fTotalYield;
133}
134
135//______________________________________________________________________________
136Double_t NaganoFluoCalculator::GetFluoYield_dE(const Double_t NF_alt, 
137                       const Double_t NF_dE, EsafSpectrum* FluoSpectrum) const {
138    //
139    // Get the wavelength spectrum of the fluorescence emission as a function
140    // of energy deposit
141    //
142    // Input:
143    //     Double_t NF_alt         the altitude where the electron(s) is(are) moving
144    //     Double_t NF_dE         the energy deposit of particle(s)
145    // Output:
146    //     fTotalYield             total fluorescence yield between 300 and 430 nm
147    //                             per electron per unit length
148    //     EsafSpectrum* FluoSpectrum fluorescence light spectrum per electron
149    //     per unit length
150    //
151
152
153    // Get Atmosphere
154    const Atmosphere* atmo = Atmosphere::Get();
155    Double_t temperature = atmo->Temperature(NF_alt);
156    Double_t pressure    = atmo->Pressure(NF_alt);
157    if(temperature<=0) FatalError("Invalid Temperature in NaganoFluoCalculator");
158   
159    // const Int_t nbin = 15;  // number of fluorescence bins
160    Double_t Yield[fNumWavelengths];
161    Double_t fTotalYield = 0.0;
162
163    for (Int_t nwl=0; nwl<fNumWavelengths; nwl++) {
164        Double_t RefPressure=fRefPressure[nwl]*sqrt(temperature/(293.0*kelvin));
165        Double_t Phi=fPhi0[nwl]/(1+pressure/RefPressure);
166        Yield[nwl]=NF_dE*Phi/(Hplanck()*Clight()/fWavelength[nwl]);
167        fTotalYield+=Yield[nwl];
168
169        fNaganoFluo->SetParameter(nwl*3,Yield[nwl]);
170        fNaganoFluo->SetParameter(nwl*3+1,fWavelength[nwl]);
171        fNaganoFluo->SetParameter(nwl*3+2,0.05*nm);
172    }
173    if ( FluoSpectrum !=0 )
174        FluoSpectrum->Reset(fNaganoFluo, 130, 300*nm, 430*nm);
175
176#ifdef DEBUG
177    Msg(EsafMsg::Debug) << "NF TotalYield/m "<<fTotalYield*m  << MsgDispatch;
178#endif
179
180    return fTotalYield;
181}
182
183/* to be completed
184//______________________________________________________________________________
185Double_t NaganoFluoCalculator::GetdE(const Double_t NF_alt,
186                                     const EsafSpectrum* FluoSpectrum) const {
187    //
188    // To be completed
189    //
190
191    // Get Atmosphere
192    const Atmosphere* atmo = Atmosphere::Get();
193    Double_t temperature = atmo->Temperature(NF_alt);
194    Double_t pressure    = atmo->Pressure(NF_alt);
195
196    Double_t fEnergyDeposit = 0.0;
197    Int_t SpectNbin = FluoSpectrum->GetAxis().GetNbins();
198    for(Int_t i=0; i<SpectNbin; i++) {
199        for(Int_t nwl=0; nwl<fNumWavelengths; nwl++) {
200            if(FluoSpectrum->GetAxis().GetBinLowEdge(i) <= fWavelength[nwl] &&
201                    fWavelength[nwl] < FluoSpectrum->GetAxis().GetBinLowEdge(i+1)  ) {
202                Double_t RefPressure=fRefPressure[nwl]*sqrt(temperature/(293.0*kelvin));
203                Double_t Phi=fPhi0[nwl]/(1+pressure/RefPressure);
204                fEnergyDeposit += FluoSpectrum->GetWeight(i) * h_Planck * c_light / fWavelength[nwl] / Phi;
205            }
206        }
207    }
208    return fEnergyDeposit;
209}
210*/
211
212//______________________________________________________________________________
213Double_t  NaganoFluoCalculator::GetFluoYield(const Double_t NF_alt, TF12* EnergyDistribution,
214                                             EsafSpectrum* FluoSpectrum) const {
215    //
216    //
217    //
218    Double_t Emax = EnergyDistribution->GetXmax();  // unit to be fixed in shower (MeV)
219    Double_t Emin = EnergyDistribution->GetXmin();
220    Double_t E=Emin;
221    Double_t dE = (Emax-Emin)/2000.;
222    Double_t fTotalYield =0.;
223
224    // get the wavelenght spectrum and total yield for 1.4 MeV
225    Double_t eneref = 1.4*MeV;
226    Double_t RefYield = GetFluoYield(NF_alt,eneref,FluoSpectrum);
227    Double_t ref_dEdX = GetdEdX( eneref );
228    Double_t y1(0.), y2(0.), value1(0.), value2(0.), integral(0.);
229   
230    // integral over the energy spectrum
231    value2 = EnergyDistribution->Eval(E);
232    y2 = GetdEdX (E*MeV) * value2;
233    for(Int_t i=0; i<2000; i++) {
234        y1 = y2;
235        E += dE;
236        value1 = value2;
237        value2 = EnergyDistribution->Eval(E);
238        y2 = GetdEdX (E*MeV) * value2;
239        fTotalYield += 0.5 * (y1 + y2);
240        integral += 0.5 * (value1 + value2);
241    }
242    fTotalYield *= RefYield / ref_dEdX / integral ;  // + energy spectrum normalization
243
244    return fTotalYield;
245}
246
247//_____________________________________________________________________________
248Double_t  NaganoFluoCalculator::GetFluoYieldHisto(const Double_t KF_alt, const TH1F* ehisto, EsafSpectrum* FluoSpectrum) const {
249    Double_t fTotalYield =0.;
250
251    // get the wavelenght spectrum and total yield for 1.4 MeV
252    Double_t eneref = 1.4*MeV;
253    Double_t RefYield = GetFluoYield(KF_alt,eneref,FluoSpectrum);
254    Double_t Int = ehisto->Integral();
255
256    if (Int == 0) return 0; 
257    Float_t dEdX_reference = GetdEdX( eneref ); 
258    for (Int_t i(0); i < ehisto->GetNbinsX(); i++) {
259        Double_t E  = ehisto->GetBinCenter(i);
260        Double_t dE = ehisto->GetBinWidth(i); 
261        fTotalYield += RefYield *  GetdEdX (E*MeV)/dEdX_reference * (ehisto->GetBinContent(i))/Int * dE;
262    }
263    return fTotalYield;
264}
265
266
267//______________________________________________________________________________
268void NaganoFluoCalculator::PlotYield() {
269    //
270    // Plot the fluorescence yield along Nadir
271    //
272/*
273    TProfile* TotalYield = (TProfile*)gROOT->FindObject("TotalYield_N");
274    if ( !TotalYield ) TotalYield = new TProfile("TotalYield_N","Total Yield along Nadir",100,0,50,0,6);
275    TProfile* Yield_337 = (TProfile*)gROOT->FindObject("YieldN_337");
276    if ( !Yield_337 ) Yield_337 = new TProfile("YieldN_337","337nm Yield along Nadir",100,0,50,0,6);
277    TProfile* Yield_357 = (TProfile*)gROOT->FindObject("YieldN_357");
278    if ( !Yield_357 ) Yield_357 = new TProfile("YieldN_357","357nm Yield along Nadir",100,0,50,0,6);
279    TProfile* Yield_391 = (TProfile*)gROOT->FindObject("YieldN_391");
280    if ( !Yield_391 ) Yield_391 = new TProfile("YieldN_391","391nm Yield along Nadir",100,0,50,0,6);
281   
282    EsafSpectrum* spectrum = new EsafSpectrum(357*nm);   
283   
284    Double_t alt,total,y337,y357,y391;
285
286    for (Double_t u=0; u<100;u++) {
287        alt=u/2*km;
288        Double_t energy = 80.*MeV;
289        total=GetFluoYield(alt,energy,spectrum)*m;
290        TotalYield->Fill(u/2,total);
291        Int_t nb = spectrum->GetAxis().GetNbins();
292        y337 = 0.;
293        y357 = 0.;
294        y391 = 0.;
295        for ( Int_t i=0; i<nb; i++) {
296          if ( 334*nm<=spectrum->GetAxis().GetBinLowEdge(i) && spectrum->GetAxis().GetBinLowEdge(i+1)<=338*nm )
297          y337 = y337 + spectrum->GetWeight(i);     
298          if ( 354*nm<=spectrum->GetAxis().GetBinLowEdge(i) && spectrum->GetAxis().GetBinLowEdge(i+1)<=358*nm )
299          y357 = y357 + spectrum->GetWeight(i);     
300          if ( 388*nm<=spectrum->GetAxis().GetBinLowEdge(i) && spectrum->GetAxis().GetBinLowEdge(i+1)<=392*nm )
301          y391 = y391 + spectrum->GetWeight(i);     
302        } 
303        Yield_337->Fill(u/2 , y337);       
304        Yield_357->Fill(u/2 , y357);
305        Yield_391->Fill(u/2 , y391);
306    }
307*/
308}
309
310//______________________________________________________________________________
311Double_t NaganoFluoCalculator::GetdEdX(const Double_t NF_energy) const {
312    //
313    // Get the dE/dX for electrons of energy KF_energy
314    // same function as in KakimotoFluoCalculator
315    //
316    // !! WARNING !! needed to be checked for high energy incident electrons (>1MeV) :
317    //               all the energy lost in dX is not totally deposited in dX (effect of delta rays)
318    //               a part of energy can be brought away by energetic electrons created by ionisation
319    //              (see Belz et al. Astropart. 25 (2006) and others)
320    //
321   Double_t energy[33]={.01,.02,.03,.04,.06,.08,.1,.2,.3,.4,.5,.6,.7,.8,.9,1.,
322                        1.25,1.5,2.,4.,6.,8.,10.,20.,40.,60.,80.,100.,200.,400.,
323                        600.,800.,1000.};
324   Double_t dEdX[33]={19.95,11.68,8.564,6.904,5.15,4.229,3.66,2.486,2.097,1.914,
325                      1.813,1.753,1.716,1.693,1.679,1.67,1.665,1.67,1.693,1.799,
326                      1.879,1.94,1.988,2.144,2.29,2.355,2.395,2.424,2.51,2.59,
327                      2.633,2.66,2.681};
328
329   Double_t dEdXE = 0;
330   for (Int_t i=0; i<33; i++) {
331       if ( energy[i]*MeV<=NF_energy && NF_energy<= energy[i+1]*MeV ) 
332      dEdXE = dEdX[i] + (dEdX[i+1]-dEdX[i])*(NF_energy-energy[i]*MeV)/(energy[i+1]*MeV-energy[i]*MeV); 
333   }
334   return dEdXE*MeV/(g/cm2);
335}
Note: See TracBrowser for help on using the repository browser.