// ESAF : Euso Simulation and Analysis Framework
// contact person : Naoto SAKAKI <n-sakaki@riken.jp>
// $Id: NaganoFluoCalculator.cc,v 1.15 2005/04/25 16:28:55 moreggia Exp $
#include <stdexcept>
#include "NaganoFluoCalculator.hh"
#include "EsafSpectrum.hh"
#include "Atmosphere.hh"
#include <TROOT.h>
#include <TProfile.h>
#include "TGraph.h"
ClassImp(NaganoFluoCalculator)
// number of fluorescence bins
const Int_t NaganoFluoCalculator::fNumWavelengths = 15;
const Double_t hPa = 100.*pascal;
// 391nm and 428nm are 1N bands and others are 2P bands
const Double_t NaganoFluoCalculator::fWavelength[] = {
316*nm, 329*nm, 337*nm, 354*nm, 358*nm,
376*nm, 381*nm, 391*nm, 394*nm, 400*nm,
406*nm, 414*nm, 420*nm, 427*nm, 428*nm};
const Double_t NaganoFluoCalculator::fPhi0[] = {
4.80e-4, 0.880e-4, 10.01e-4, 0.769e-4, 7.82e-4,
1.20e-4, 2.46e-4, 9.60e-4, 0.42e-4, 0.847e-4,
1.49e-4, 0.327e-4, 0.86e-4, 0.069e-4, 4.57e-4};
const Double_t NaganoFluoCalculator::fRefPressure[] = {
23.0*hPa, 40.2*hPa, 19.2*hPa, 30.6*hPa, 18.1*hPa,
34.1*hPa, 19.4*hPa, 5.02*hPa,24.2*hPa, 24.2*hPa,
12.3*hPa, 19.3*hPa, 7.3*hPa, 72.*hPa, 3.86*hPa};
//______________________________________________________________________________
NaganoFluoCalculator::NaganoFluoCalculator() {
//
// Constructor
//
fName = "Nagano";
Msg(EsafMsg::Info) << "NaganoFluoCalculator built" << MsgDispatch;
TString formula = "gaus(0)";
for (Int_t i(1); i<fNumWavelengths; i++)
formula += Form("+gaus(%d)",i*3);
fNaganoFluo = new TFormula("fluo - nagano", formula );
}
//______________________________________________________________________________
NaganoFluoCalculator::~NaganoFluoCalculator() {
//
// Destructor
//
SafeDelete(fNaganoFluo);
}
//______________________________________________________________________________
Double_t NaganoFluoCalculator::GetFluoYield(const Double_t NF_alt,
const Double_t NF_energy, EsafSpectrum* FluoSpectrum) const {
//
// Get the wavelength spectrum of the fluorescence emission as a function
// of electron energy
//
// Input:
// Double_t NF_alt the altitude where the electron is moving
// Double_t NF_energy the energy of the electron
// Output:
// fTotalYield total fluorescence yield between 300 and 430 nm
// per electron per unit length
// EsafSpectrum* FluoSpectrum fluorescence light spectrum per electron
// per unit length
//
// [1] Nagano et al. Astroparticle Physics, 20 (2003) 293.(astro-ph/0303193)
// [2] Nagano et al. astro-ph/0406474 (2004)
// measurements of yield for 0.85MeV electrons
#ifdef DEBUG
Msg(EsafMsg::Debug) << "NF "<<NF_alt/km << " " <<NF_energy/MeV << MsgDispatch;
#endif
// Get Atmosphere
const Atmosphere* atmo = Atmosphere::Get();
Double_t density = atmo->Air_Density(NF_alt);
Double_t temperature = atmo->Temperature(NF_alt);
if(temperature<=0)
FatalError("Invalid Temperature in NaganoFluoCalculator");
const Double_t dedx_ref = GetdEdX( 0.85*MeV );
// density and temperature dependance from Nagano et al.
const static Double_t A[] = {
20.5*m2/kg, 3.91*m2/kg, 45.6*m2/kg, 3.68*m2/kg, 37.8*m2/kg,
6.07*m2/kg, 12.7*m2/kg, 50.8*m2/kg, 2.25*m2/kg, 4.58*m2/kg,
8.18*m2/kg, 1.83*m2/kg, 4.9*m2/kg, 0.40*m2/kg, 26.5*m2/kg
};
const static Double_t B[] = {
2.14*m3/kg, 1.22*m3/kg, 2.56*m3/kg, 1.60*m3/kg, 2.72*m3/kg,
1.44*m3/kg, 2.53*m3/kg, 9.80*m3/kg, 2.03*m3/kg, 2.03*m3/kg,
3.99*m3/kg, 2.55*m3/kg, 6.8 *m3/kg, 0.68*m3/kg, 12.7*m3/kg
};
Double_t fTotalYield = 0.0;
Double_t Yield[fNumWavelengths];
for (Int_t i=0; i<fNumWavelengths; i++) {
Yield[i] = GetdEdX( NF_energy ) / dedx_ref
*density*A[i]/(1+density*B[i]*sqrt(temperature));
fTotalYield += Yield[i];
fNaganoFluo->SetParameter(i*3,Yield[i]);
fNaganoFluo->SetParameter(i*3+1,fWavelength[i]);
fNaganoFluo->SetParameter(i*3+2,0.05*nm);
}
if ( FluoSpectrum !=0 )
FluoSpectrum->Reset(fNaganoFluo, 130, 300*nm, 430*nm);
#ifdef DEBUG
Msg(EsafMsg::Debug) << "NF TotalYield/m "<<fTotalYield*m << MsgDispatch;
#endif
return fTotalYield;
}
//______________________________________________________________________________
Double_t NaganoFluoCalculator::GetFluoYield_dE(const Double_t NF_alt,
const Double_t NF_dE, EsafSpectrum* FluoSpectrum) const {
//
// Get the wavelength spectrum of the fluorescence emission as a function
// of energy deposit
//
// Input:
// Double_t NF_alt the altitude where the electron(s) is(are) moving
// Double_t NF_dE the energy deposit of particle(s)
// Output:
// fTotalYield total fluorescence yield between 300 and 430 nm
// per electron per unit length
// EsafSpectrum* FluoSpectrum fluorescence light spectrum per electron
// per unit length
//
// Get Atmosphere
const Atmosphere* atmo = Atmosphere::Get();
Double_t temperature = atmo->Temperature(NF_alt);
Double_t pressure = atmo->Pressure(NF_alt);
if(temperature<=0) FatalError("Invalid Temperature in NaganoFluoCalculator");
// const Int_t nbin = 15; // number of fluorescence bins
Double_t Yield[fNumWavelengths];
Double_t fTotalYield = 0.0;
for (Int_t nwl=0; nwl<fNumWavelengths; nwl++) {
Double_t RefPressure=fRefPressure[nwl]*sqrt(temperature/(293.0*kelvin));
Double_t Phi=fPhi0[nwl]/(1+pressure/RefPressure);
Yield[nwl]=NF_dE*Phi/(h_Planck*c_light/fWavelength[nwl]);
fTotalYield+=Yield[nwl];
fNaganoFluo->SetParameter(nwl*3,Yield[nwl]);
fNaganoFluo->SetParameter(nwl*3+1,fWavelength[nwl]);
fNaganoFluo->SetParameter(nwl*3+2,0.05*nm);
}
if ( FluoSpectrum !=0 )
FluoSpectrum->Reset(fNaganoFluo, 130, 300*nm, 430*nm);
#ifdef DEBUG
Msg(EsafMsg::Debug) << "NF TotalYield/m "<<fTotalYield*m << MsgDispatch;
#endif
return fTotalYield;
}
/* to be completed
//______________________________________________________________________________
Double_t NaganoFluoCalculator::GetdE(const Double_t NF_alt,
const EsafSpectrum* FluoSpectrum) const {
//
// To be completed
//
// Get Atmosphere
const Atmosphere* atmo = Atmosphere::Get();
Double_t temperature = atmo->Temperature(NF_alt);
Double_t pressure = atmo->Pressure(NF_alt);
Double_t fEnergyDeposit = 0.0;
Int_t SpectNbin = FluoSpectrum->GetAxis().GetNbins();
for(Int_t i=0; i<SpectNbin; i++) {
for(Int_t nwl=0; nwl<fNumWavelengths; nwl++) {
if(FluoSpectrum->GetAxis().GetBinLowEdge(i) <= fWavelength[nwl] &&
fWavelength[nwl] < FluoSpectrum->GetAxis().GetBinLowEdge(i+1) ) {
Double_t RefPressure=fRefPressure[nwl]*sqrt(temperature/(293.0*kelvin));
Double_t Phi=fPhi0[nwl]/(1+pressure/RefPressure);
fEnergyDeposit += FluoSpectrum->GetWeight(i) * h_Planck * c_light / fWavelength[nwl] / Phi;
}
}
}
return fEnergyDeposit;
}
*/
//______________________________________________________________________________
Double_t NaganoFluoCalculator::GetFluoYield(const Double_t NF_alt, TF12* EnergyDistribution,
EsafSpectrum* FluoSpectrum) const {
//
//
//
Double_t Emax = EnergyDistribution->GetXmax(); // unit to be fixed in shower (MeV)
Double_t Emin = EnergyDistribution->GetXmin();
Double_t E=Emin;
Double_t dE = 0.5*MeV;
Double_t fTotalYield =0.;
// get the wavelength spectrum and total yield for 1.4 MeV
Double_t eneref = 1.4*MeV;
Double_t RefYield = GetFluoYield(NF_alt,eneref,FluoSpectrum);
Double_t Integral = EnergyDistribution->Integral(Emin,Emax);
if (Integral == 0) return 0;
Float_t dEdX_reference = GetdEdX( eneref );
// integral over the energy spectrum
while(kTRUE) {
fTotalYield += GetdEdX (E*MeV)*EnergyDistribution->Eval(E);
E += dE;
if(( Emax-Emin )<( E-Emin )) break;
}
fTotalYield *= RefYield/dEdX_reference*dE/Integral;
return fTotalYield;
}
//_____________________________________________________________________________
Double_t NaganoFluoCalculator::GetFluoYieldHisto(const Double_t KF_alt, const TH1F* ehisto, EsafSpectrum* FluoSpectrum) const {
Double_t fTotalYield =0.;
// get the wavelenght spectrum and total yield for 1.4 MeV
Double_t eneref = 1.4*MeV;
Double_t RefYield = GetFluoYield(KF_alt,eneref,FluoSpectrum);
Double_t Int = ehisto->Integral();
if (Int == 0) return 0;
Float_t dEdX_reference = GetdEdX( eneref );
for (Int_t i(0); i < ehisto->GetNbinsX(); i++) {
Double_t E = ehisto->GetBinCenter(i);
Double_t dE = ehisto->GetBinWidth(i);
fTotalYield += RefYield * GetdEdX (E*MeV)/dEdX_reference * (ehisto->GetBinContent(i))/Int * dE;
}
return fTotalYield;
}
//______________________________________________________________________________
void NaganoFluoCalculator::PlotYield() {
//
// Plot the fluorescence yield along Nadir
//
/*
TProfile* TotalYield = (TProfile*)gROOT->FindObject("TotalYield_N");
if ( !TotalYield ) TotalYield = new TProfile("TotalYield_N","Total Yield along Nadir",100,0,50,0,6);
TProfile* Yield_337 = (TProfile*)gROOT->FindObject("YieldN_337");
if ( !Yield_337 ) Yield_337 = new TProfile("YieldN_337","337nm Yield along Nadir",100,0,50,0,6);
TProfile* Yield_357 = (TProfile*)gROOT->FindObject("YieldN_357");
if ( !Yield_357 ) Yield_357 = new TProfile("YieldN_357","357nm Yield along Nadir",100,0,50,0,6);
TProfile* Yield_391 = (TProfile*)gROOT->FindObject("YieldN_391");
if ( !Yield_391 ) Yield_391 = new TProfile("YieldN_391","391nm Yield along Nadir",100,0,50,0,6);
EsafSpectrum* spectrum = new EsafSpectrum(357*nm);
Double_t alt,total,y337,y357,y391;
for (Double_t u=0; u<100;u++) {
alt=u/2*km;
Double_t energy = 80.*MeV;
total=GetFluoYield(alt,energy,spectrum)*m;
TotalYield->Fill(u/2,total);
Int_t nb = spectrum->GetAxis().GetNbins();
y337 = 0.;
y357 = 0.;
y391 = 0.;
for ( Int_t i=0; i<nb; i++) {
if ( 334*nm<=spectrum->GetAxis().GetBinLowEdge(i) && spectrum->GetAxis().GetBinLowEdge(i+1)<=338*nm )
y337 = y337 + spectrum->GetWeight(i);
if ( 354*nm<=spectrum->GetAxis().GetBinLowEdge(i) && spectrum->GetAxis().GetBinLowEdge(i+1)<=358*nm )
y357 = y357 + spectrum->GetWeight(i);
if ( 388*nm<=spectrum->GetAxis().GetBinLowEdge(i) && spectrum->GetAxis().GetBinLowEdge(i+1)<=392*nm )
y391 = y391 + spectrum->GetWeight(i);
}
Yield_337->Fill(u/2 , y337);
Yield_357->Fill(u/2 , y357);
Yield_391->Fill(u/2 , y391);
}
*/
}
//______________________________________________________________________________
Double_t NaganoFluoCalculator::GetdEdX(const Double_t NF_energy) const {
//
// Get the dE/dX for electrons of energy KF_energy
// same function as in KakimotoFluoCalculator
//
Double_t energy[33]={.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.};
Double_t dEdX[33]={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};
Double_t dEdXE;
for (Int_t i=0; i<33; i++) {
if ( energy[i]*MeV<=NF_energy && NF_energy<= energy[i+1]*MeV )
dEdXE = dEdX[i] + (dEdX[i+1]-dEdX[i])*(NF_energy-energy[i]*MeV)/(energy[i+1]*MeV-energy[i]*MeV);
}
return dEdXE*MeV/(g/cm2);
}