source: JEM-EUSO/esaf_lal/tags/v1_r0/esaf/packages/simulation/lightsources/src/SimpleCrkCalculator.cc @ 117

// ESAF : Euso Simulation and Analysis Framework
// $Id: SimpleCrkCalculator.cc 2992 2011-09-28 10:22:44Z fenu $
// Anne Stutz created Apr, 27 2004

#include "SimpleCrkCalculator.hh"
#include "EsafSpectrum.hh"
#include "Atmosphere.hh"
#include "Config.hh"
#include "EConst.hh"

using namespace sou;
using namespace TMath;
using EConst::FineStructureConst;
using EConst::ElectronMassC2;

ClassImp(SimpleCrkCalculator)

//____________________________________________________________________________________________
  SimpleCrkCalculator::SimpleCrkCalculator() : CrkCalculator() {
    //
    // ctor
    //
    ConfigFileParser* pConf = Config::Get()->GetCF("LightSource", "ShowerLightSource");        
    fLambdaMin = (Int_t) pConf->GetNum("ShowerLightSource.fLambdaMin")*nm;
    fLambdaMax = (Int_t) pConf->GetNum("ShowerLightSource.fLambdaMax")*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, Bool_t thr) const {
    //
    // get the yield of the Cerenkov emission for a fixed energy
    // if thr == false, energy threshold of cerenkov production is not taken into account
    //
    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)) && thr) {
        CrkYield = TwoPi() * FineStructureConst() * (2*delta - pow(ElectronMassC2()/SC_energy,2));
        // wavelenght normalisation
    }
    else if(!thr) CrkYield = TwoPi() * FineStructureConst() * 2*delta * (1 - pow(GetEnergyThreshold(0.)/SC_energy,2));
    else CrkYield=0.;
    CrkYield = CrkYield * ( 1/fLambdaMin - 1/fLambdaMax );

    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 = EnergyDistribution->Integral(Emin,Emax);

    while(true) {
        CrkYield += GetCrkYield(SC_alt,E*MeV) * EnergyDistribution->Eval(E) * dE ; // unit to be fixed in shower
        E += dE;
        if(( Emax-Emin )<( E-Emin )) break;
    }
    CrkYield /= IntSpec;

    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 ElectronMassC2()/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;
}