// ESAF : Euso Simulation and Analysis Framework
// $Id: SimpleCrkCalculator.cc,v 1.14 2005/02/11 15:30:47 moreggia Exp $
// Anne Stutz created Apr, 27 2004
#include "SimpleCrkCalculator.hh"
#include "EsafSpectrum.hh"
#include "Atmosphere.hh"
ClassImp(SimpleCrkCalculator)
//____________________________________________________________________________________________
SimpleCrkCalculator::SimpleCrkCalculator() : CrkCalculator() {
//
// ctor
//
fLambdaMin = 300*nm;
fLambdaMax = 450*nm;
fName = "simple";
Msg(EsafMsg::Info)<< "SimpleCrkCalculator built" << MsgDispatch;
}
//____________________________________________________________________________________________
SimpleCrkCalculator::~SimpleCrkCalculator() {
//
// dtor
//
}
//____________________________________________________________________________________________
Double_t SimpleCrkCalculator::GetCrkYield(const Double_t SC_alt, const Double_t SC_energy) const {
//
// get the yield of the Cerenkov emission for a fixed energy
//
Double_t CrkYield;
// Get Atmosphere for Index
const Atmosphere* atmo = Atmosphere::Get();
Double_t delta = atmo->Index_Minus1(SC_alt);
// When Index is not implemented in the Atmosphere use the Corsika formula
if ( delta == 0)
delta = 0.000283 * atmo->Air_Density( SC_alt )/atmo->Air_Density(0);
if(delta==0) throw runtime_error("Invalid Index in SimpleCrkCalculator");
// Cerenkov Yield
if ( SC_energy > GetEnergyThreshold(SC_alt)) {
CrkYield = twopi * fine_structure_const * (2*delta - pow(electron_mass_c2/SC_energy,2));
// wavelenght normalisation
CrkYield = CrkYield * ( 1/fLambdaMin - 1/fLambdaMax );
}
else CrkYield=0.;
return CrkYield;
}
//____________________________________________________________________________________________
Double_t SimpleCrkCalculator::GetCrkYield(const Double_t SC_alt, TF12* EnergyDistribution) const {
//
// get the Cerenkov yield integrated on an energy distribution
//
Double_t CrkYield = 0.;
Double_t Emax = EnergyDistribution->GetXmax(); // unit to be fixed in shower
Double_t Emin = EnergyDistribution->GetXmin();
Double_t E = GetEnergyThreshold(SC_alt)/MeV;
Double_t dE = (Emax-E)/1000.;
Double_t IntSpec=0;
while(true) {
CrkYield += GetCrkYield(SC_alt,E*MeV) * EnergyDistribution->Eval(E) * dE ; // unit to be fixed in shower
IntSpec += EnergyDistribution->Eval(E) * dE;
E += dE;
if(( Emax-Emin )<( E-Emin )) break;
}
return CrkYield;
}
//____________________________________________________________________________________________
Double_t SimpleCrkCalculator::GetCrkYield(const Double_t SC_alt, const TH1F* ehisto) const {
//
// get the Cerenkov yield integrated on an energy histogram
//
Double_t CrkYield = 0.;
Double_t Int = ehisto->Integral();
if (Int == 0) return 0;
for (Int_t i(0); i < ehisto->GetNbinsX(); i++) {
Double_t E = ehisto->GetBinCenter(i);
Double_t dE = ehisto->GetBinWidth(i);
CrkYield +=GetCrkYield(SC_alt,E*MeV) * (ehisto->GetBinContent(i))/Int *dE ;
}
return CrkYield;
}
//____________________________________________________________________________________________
Double_t SimpleCrkCalculator::GetEnergyThreshold(const Double_t SC_alt) const {
//
// get the energy threshold for the cerenkov emission
//
// Get Atmosphere for Index
const Atmosphere* atmo = Atmosphere::Get();
Double_t delta = atmo->Index_Minus1(SC_alt);
// When Index is not implemented in the Atmosphere use the Corsika formula
if ( delta == 0)
delta = 0.000283 * atmo->Air_Density( SC_alt )/atmo->Air_Density(0);
return electron_mass_c2/sqrt(2*delta);
}
//_________________________________________________________________________________
EsafSpectrum* SimpleCrkCalculator::GetCrkSpectrum() const{
//
// get the wavelenght spectrum of the cerenkov emission
//
TFormula* cerenkov = new TFormula("cerenkov","1/(x*x)");
EsafSpectrum* CrkSpectrum = new EsafSpectrum(cerenkov,100,fLambdaMin,fLambdaMax);
delete cerenkov;
return CrkSpectrum;
}