source: JEM-EUSO/esaf_cc_at_lal/packages/simulation/lightsources/src/KakimotoFluoCalculator.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: 9.5 KB
Line 
1// ESAF : Euso Simulation and Analysis Framework
2// $Id: KakimotoFluoCalculator.cc 2776 2006-11-23 16:53:19Z moreggia $
3// Anne Stutz created Mar, 23 2004
4
5#include <stdexcept>
6#include "KakimotoFluoCalculator.hh"
7#include "EsafSpectrum.hh"
8#include "Atmosphere.hh"
9#include <TROOT.h>
10#include <TProfile.h>
11#include "TGraph.h"
12
13using namespace sou;
14
15ClassImp(KakimotoFluoCalculator)
16
17//________________________________________________________________________________________________________________________
18KakimotoFluoCalculator::KakimotoFluoCalculator() : FluoCalculator() {
19    //
20    // ctor
21    //
22   
23    fName = "Kakimoto";
24    Msg(EsafMsg::Info) << "KakimotoluoCalculator built" << MsgDispatch;
25    //    PlotYield();
26}
27
28//________________________________________________________________________________________________________________________
29KakimotoFluoCalculator::~KakimotoFluoCalculator() {
30    //
31    // dtor
32    //
33}
34
35//________________________________________________________________________________________________________________________
36Double_t KakimotoFluoCalculator::GetFluoYield(const Double_t KF_alt, const Double_t KF_energy, EsafSpectrum* FluoSpectrum) const {
37    //
38    // get the wavelenght spectrum of the fluorescence emission
39    // Use Kakimoto measurement of 3 main peaks, supplemented with Bunner lines for the remaining yield between 300-400 nm
40    // See Auger GapNote 2002-067, Bruce Dawson 2004 : Model used by HiRes
41    //
42    // NB : Kakimoto used Davidson spectrum to reconstruct total yield within [300-400], present model uses Bunner..
43    //
44
45    // Get Atmosphere
46    const Atmosphere* atmo = Atmosphere::Get();
47
48    // 1. DATA INPUTS
49    //    Kakimoto et al. measurements of yield at 760 mmHg, 15°C for 1.4 MeV electrons,
50    //    supplemented with Bunner spectrum
51    Int_t nbin = 12;  // number of fluorescence bins
52    Double_t Wavelenght[] = {298,   316,   325,  334,   343,   352,   361,  370,   379,   388,   397,   406};
53    Double_t RefYield[]   = {0.017, 0.118, 0.01, 1.109, 0.005, 0.058, 1.19, 0.017, 0.152, 0.495, 0.045, 0.035};
54    for(Int_t j=0; j < nbin; j++) {
55        Wavelenght[j] *= nm;
56        RefYield[j]   /= m;
57    }
58    Double_t total2P = 0.;
59    Double_t total1N = RefYield[9];
60    for(Int_t j=0; j < nbin; j++) total2P += RefYield[j];
61    total2P -= total1N;
62   
63    // 1.1 density and temperature dependance from Kakimoto et al. NIMA372(1996)527-533
64    Double_t A1 = 89.0*m2/kg;   // +/- 1.7 m(2)/kg
65    Double_t A2 = 55.0*m2/kg;   //
66    Double_t B1 = 1.85*m3/kg;   // +/- 0.04 m(3)/kg.K(-1/2)
67    Double_t B2 = 6.50*m3/kg;   //
68
69 
70    // 2. get atmosphere parameters
71    Double_t density = atmo->Air_Density(KF_alt);
72    Double_t temp = atmo->Temperature(KF_alt);
73    if(temp == 0.) Msg(EsafMsg::Panic) << "Invalid Temperature in KakimotoFluoCalculator"<<MsgDispatch;
74   
75   
76    // 3. energy dependence
77    Double_t energy = KF_energy;
78    Double_t dedx_ref = GetdEdX(1.4*MeV);
79    for(Int_t i=0; i < nbin; i++) RefYield[i] *= GetdEdX(energy) / dedx_ref;
80
81
82    // 4. build spectrum
83    Double_t Yield[nbin];
84    Double_t fTotalYield = 0.;
85    string KF_string = "gaus(0)";
86    char gaussian[10];
87    for(Int_t i=1; i < nbin; i++) {
88        sprintf(gaussian,"+gaus(%d)",i*3);
89        KF_string += gaussian;
90    }
91    const char* KF_char = KF_string.c_str();
92    TFormula* kakimoto = new TFormula("fluo - kakimoto",KF_char);
93
94    for(Int_t i=0; i < nbin; i++) {
95        if(i != 9) Yield[i] = RefYield[i]/total2P * density*A1/(1 + density*B1*sqrt(temp)); // normalized cause A1 already contains the total 2P yield
96        else       Yield[i] = RefYield[i]/total1N * density*A2/(1 + density*B2*sqrt(temp)); // normalized cause A2 already contains the total 1N yield
97
98        fTotalYield += Yield[i];
99
100        kakimoto->SetParameter(i*3,Yield[i]);
101        kakimoto->SetParameter(i*3+1,Wavelenght[i]);
102        kakimoto->SetParameter(i*3+2,0.05*nm);
103    }
104    if ( FluoSpectrum !=0 ) FluoSpectrum->Reset(kakimoto, 130, 280*nm, 410*nm);
105    SafeDelete(kakimoto);
106   
107    return fTotalYield;
108}
109
110//________________________________________________________________________________________________________________________
111Double_t  KakimotoFluoCalculator::GetFluoYield(const Double_t KF_alt, TF12* EnergyDistribution, EsafSpectrum* FluoSpectrum) const {
112    //
113    // integration over energy spectrum
114    //
115    Double_t Emax = EnergyDistribution->GetXmax();  // unit to be fixed in shower (MeV)
116    Double_t Emin = EnergyDistribution->GetXmin();
117    Double_t E=Emin;
118    Double_t dE = (Emax-Emin)/2000.;
119    Double_t fTotalYield =0.;
120
121    // get the wavelenght spectrum and total yield for 1.4 MeV
122    Double_t eneref = 1.4*MeV;
123    Double_t RefYield = GetFluoYield(KF_alt,eneref,FluoSpectrum);
124    Double_t ref_dEdX = GetdEdX( eneref );
125    Double_t y1(0.), y2(0.), value1(0.), value2(0.), integral(0.);
126   
127    // integral over the energy spectrum
128    value2 = EnergyDistribution->Eval(E);
129    y2 = GetdEdX (E*MeV) * value2;
130    for(Int_t i=0; i<2000; i++) {
131        y1 = y2;
132        E += dE;
133        value1 = value2;
134        value2 = EnergyDistribution->Eval(E);
135        y2 = GetdEdX (E*MeV) * value2;
136        fTotalYield += 0.5 * (y1 + y2);
137        integral += 0.5 * (value1 + value2);
138    }
139    // energy spectrum normalization : dE does not counted in integration cause present in both integral and fTotalYield
140    fTotalYield *= RefYield / ref_dEdX / integral ; 
141
142
143    return fTotalYield;
144}
145
146//________________________________________________________________________________________________________________________
147Double_t  KakimotoFluoCalculator::GetFluoYieldHisto(const Double_t KF_alt, const TH1F* ehisto, EsafSpectrum* FluoSpectrum) const {
148    //
149    // integration over energy spectrum
150    //
151    Double_t fTotalYield =0.;
152    Double_t Int=ehisto->Integral();
153    if (Int==0) return 0;
154
155    // get the wavelenght spectrum and total yield for 1.4 MeV
156    Double_t eneref = 1.4*MeV;
157    Double_t RefYield = GetFluoYield(KF_alt,eneref,FluoSpectrum);
158
159    for (Int_t i(0); i < ehisto->GetNbinsX(); i++) {
160        Double_t E  = ehisto->GetBinCenter(i);
161        Double_t dE = ehisto->GetBinWidth(i); 
162        fTotalYield += RefYield *  GetdEdX (E*MeV)/ GetdEdX( eneref ) * ehisto->GetBinContent(i) * dE;
163    }
164    fTotalYield /= Int;   // energy spectrum normalization
165
166    return fTotalYield;
167}
168
169//________________________________________________________________________________________________________________________
170Double_t KakimotoFluoCalculator::GetdEdX(const Double_t KF_energy) const {
171    //
172    // get the dE/dX for electrons of energy KF_energy
173    //
174    // !! WARNING !! needed to be checked for high energy incident electrons (>1MeV) :
175    //               all the energy lost in dX is not totally deposited in dX (effect of delta rays)
176    //               a part of energy can be brought away by energetic electrons created by ionisation
177    //              (see Belz et al. Astropart. 25 (2006) and others)
178    //
179    Double_t energy[33]=
180 {.01  , .02  , .03  , .04  , .06, .08   , .1  ,.2   , .3    , .4   , .5   , .6   , .7   , .8   , .9   , 1.  , 1.25 , 1.5 , 2.   , 4.   , 6.   , 8.  , 10.  , 20.  , 40. , 60.  , 80.  , 100. ,200.,400.,600.,800.,1000.};
181    Double_t dEdX[33]=
182 {19.95, 11.68, 8.564, 6.904, 5.15, 4.229, 3.66, 2.486, 2.097, 1.914, 1.813, 1.753, 1.716, 1.693, 1.679, 1.67, 1.665, 1.67, 1.693, 1.799, 1.879, 1.94, 1.988, 2.144, 2.29, 2.355, 2.395, 2.424,2.51,2.59,2.633,2.66,2.681};
183
184    Double_t dEdXE = 0;
185    for (Int_t i=0; i<33; i++) {
186        if ( energy[i]*MeV<=KF_energy && KF_energy<= energy[i+1]*MeV ) 
187           dEdXE = dEdX[i] + (dEdX[i+1]-dEdX[i])*(KF_energy-energy[i]*MeV)/(energy[i+1]*MeV-energy[i]*MeV); 
188    }
189    return dEdXE;
190}
191
192//________________________________________________________________________________________________________________________
193void KakimotoFluoCalculator::PlotYield() {
194    //
195    // plot the fluorescence yield along Nadir
196    //
197/*
198    TProfile* TotalYield = (TProfile*)gROOT->FindObject("TotalYield");
199    if ( !TotalYield ) TotalYield = new TProfile("TotalYield","Total Yield along Nadir",100,0,50,0,6);
200    TProfile* Yield_337 = (TProfile*)gROOT->FindObject("Yield_337");
201    if ( !Yield_337 ) Yield_337 = new TProfile("Yield_337","337nm Yield along Nadir",100,0,50,0,6);
202    TProfile* Yield_357 = (TProfile*)gROOT->FindObject("Yield_357");
203    if ( !Yield_357 ) Yield_357 = new TProfile("Yield_357","357nm Yield along Nadir",100,0,50,0,6);
204    TProfile* Yield_391 = (TProfile*)gROOT->FindObject("Yield_391");
205    if ( !Yield_391 ) Yield_391 = new TProfile("Yield_391","391nm Yield along Nadir",100,0,50,0,6);
206   
207    EsafSpectrum* spectrum = new EsafSpectrum(357*nm);   
208   
209    Double_t alt,total,y337,y357,y391;
210
211    for (Double_t u=0; u<100;u++) {
212        alt=u/2*km;
213        Double_t energy = 80.*MeV;
214        total=GetFluoYield(alt,energy,spectrum)*m;
215        TotalYield->Fill(u/2,total);
216        Int_t nb = spectrum->GetAxis().GetNbins();
217        y337 = 0.;
218        y357 = 0.;
219        y391 = 0.;
220        for ( Int_t i=0; i<nb; i++) {
221          if ( 334*nm<=spectrum->GetAxis().GetBinLowEdge(i) && spectrum->GetAxis().GetBinLowEdge(i+1)<=338*nm )
222              y337 = y337 + spectrum->GetWeight(i);     
223          if ( 354*nm<=spectrum->GetAxis().GetBinLowEdge(i) && spectrum->GetAxis().GetBinLowEdge(i+1)<=358*nm )
224              y357 = y357 + spectrum->GetWeight(i);       
225          if ( 388*nm<=spectrum->GetAxis().GetBinLowEdge(i) && spectrum->GetAxis().GetBinLowEdge(i+1)<=392*nm )
226              y391 = y391 + spectrum->GetWeight(i) ;     
227        } 
228        Yield_337->Fill(u/2 , y337);       
229        Yield_357->Fill(u/2 , y357);
230        Yield_391->Fill(u/2 , y391);
231    }
232*/
233}
Note: See TracBrowser for help on using the repository browser.