// ESAF : Euso Simulation and Analysis Framework
// $Id: KakimotoFluoCalculator.cc,v 1.16 2005/05/09 17:34:46 moreggia Exp $
// Anne Stutz created Mar, 23 2004
#include <stdexcept>
#include "KakimotoFluoCalculator.hh"
#include "EsafSpectrum.hh"
#include "Atmosphere.hh"
#include <TROOT.h>
#include <TProfile.h>
#include "TGraph.h"
ClassImp(KakimotoFluoCalculator)
// ctor
KakimotoFluoCalculator::KakimotoFluoCalculator() : FluoCalculator() {
fName = "Kakimoto";
Msg(EsafMsg::Info) << "KakimotoluoCalculator built" << MsgDispatch;
// PlotYield();
}
// dtor
KakimotoFluoCalculator::~KakimotoFluoCalculator() {
}
//________________________________________________________________________________________________________________________
Double_t KakimotoFluoCalculator::GetFluoYield(const Double_t KF_alt, const Double_t KF_energy, EsafSpectrum* FluoSpectrum) const {
//
// get the wavelenght spectrum of the fluorescence emission
//
// Get Atmosphere
const Atmosphere* atmo = Atmosphere::Get();
// Kakimoto et al. measurements of yield at 760 mmHg, 15°C for 1.4 MeV electrons
Int_t nbin = 3; // number of fluorescence bins
Double_t Wavelenght[] = {337*nm,357*nm,391*nm};
Double_t RefYield[]= {1.109/m,1.19/m,0.495/m};
Double_t total2P = RefYield[0]+RefYield[1];
Double_t total1N = RefYield[2];
Double_t density = atmo->Air_Density(KF_alt);
Double_t temp = atmo->Temperature(KF_alt);
if(temp==0)
Msg(EsafMsg::Panic) << "Invalid Temperature in KakimotoFluoCalculator"<<MsgDispatch;
// energy dependance
Double_t eneref = 1.4*MeV;
Double_t energy = KF_energy;
for( Int_t i=0; i<nbin; i++) RefYield[i] = GetdEdX (energy )/ GetdEdX( eneref ) * RefYield[i];
// density and temperature dependance from Kakimoto et al. NIMA372(1996)527-533
Double_t A1 = 89.0*m2/kg; // +/- 1.7 m(2)/kg
Double_t A2 = 55.0*m2/kg; //
Double_t B1 = 1.85*m3/kg; // +/- 0.04 m(3)/kg.K(-1/2)
Double_t B2 = 6.50*m3/kg; //
Double_t Yield[nbin];
Double_t fTotalYield = 0.0;
string KF_string = "gaus(0)";
char* gaussian = 0;
for (Int_t i=1; i<nbin; i++) {
gaussian = new char[8];
sprintf(gaussian,"+gaus(%d)",i*3);
KF_string+=gaussian;
delete[] gaussian;
gaussian = 0;
}
const char* KF_char = KF_string.c_str();
TFormula* kakimoto = new TFormula("fluo - kakimoto",KF_char);
for (Int_t i=0; i<nbin; i++) {
if ( Wavelenght[i] != 391*nm )
Yield[i] =RefYield[i]/total2P
*density*A1/(1+density*B1*sqrt(temp));
else
Yield[i] = RefYield[i]/total1N
*density*A2/(1+density*B2*sqrt(temp));
fTotalYield = fTotalYield + Yield[i];
kakimoto->SetParameter(i*3,Yield[i]);
kakimoto->SetParameter(i*3+1,Wavelenght[i]);
kakimoto->SetParameter(i*3+2,0.05*nm);
}
if ( FluoSpectrum !=0 ) FluoSpectrum->Reset(kakimoto, 100, 300*nm, 400*nm);
// cout << "KF TotalYield/m "<<fTotalYield*m << endl;
SafeDelete(kakimoto);
return fTotalYield;
}
//________________________________________________________________________________________________________________________
Double_t KakimotoFluoCalculator::GetFluoYield(const Double_t KF_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 = (Emax-Emin)/1000.;
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);
// integral over the energy spectrum
while(true) {
fTotalYield += RefYield * GetdEdX (E*MeV)/ GetdEdX( eneref ) * EnergyDistribution->Eval(E) * dE;
E += dE;
if(( Emax-Emin )<( E-Emin )) break;
}
return fTotalYield;
}
//________________________________________________________________________________________________________________________
Double_t KakimotoFluoCalculator::GetFluoYieldHisto(const Double_t KF_alt, const TH1F* ehisto, EsafSpectrum* FluoSpectrum) const {
Double_t fTotalYield =0.;
Double_t Int=ehisto->Integral();
if (Int==0) return 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);
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)/ GetdEdX( eneref ) * (ehisto->GetBinContent(i))/Int * dE;
}
return fTotalYield;
}
//________________________________________________________________________________________________________________________
void KakimotoFluoCalculator::PlotYield() {
//
// plot the fluorescence yield along Nadir
//
/*
TProfile* TotalYield = (TProfile*)gROOT->FindObject("TotalYield");
if ( !TotalYield ) TotalYield = new TProfile("TotalYield","Total Yield along Nadir",100,0,50,0,6);
TProfile* Yield_337 = (TProfile*)gROOT->FindObject("Yield_337");
if ( !Yield_337 ) Yield_337 = new TProfile("Yield_337","337nm Yield along Nadir",100,0,50,0,6);
TProfile* Yield_357 = (TProfile*)gROOT->FindObject("Yield_357");
if ( !Yield_357 ) Yield_357 = new TProfile("Yield_357","357nm Yield along Nadir",100,0,50,0,6);
TProfile* Yield_391 = (TProfile*)gROOT->FindObject("Yield_391");
if ( !Yield_391 ) Yield_391 = new TProfile("Yield_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 KakimotoFluoCalculator::GetdEdX(const Double_t KF_energy) const {
//
// get the dE/dX for electrons of energy KF_energy
//
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<=KF_energy && KF_energy<= energy[i+1]*MeV )
dEdXE = dEdX[i] + (dEdX[i+1]-dEdX[i])*(KF_energy-energy[i]*MeV)/(energy[i+1]*MeV-energy[i]*MeV);
}
return dEdXE;
}