Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

ShowerLightSource - source file

// ESAF : Euso Simulation and Analysis Framework
// $Id: ShowerLightSource.cc,v 1.62 2005/10/27 13:54:35 moreggia Exp $
// Alessandro Thea created Nov, 24 2003

#include "ShowerLightSource.hh"
#include "Config.hh"

#include "ListPhotonsInAtmosphere.hh"
#include "FluoCalculator.hh"
#include "CrkCalculator.hh"
#include "LightSourceFactory.hh"

#include "Atmosphere.hh"

#include "EAtmosphere.hh"
#include "EAtmosphereBunchAdder.hh"
#include "utils.hh"
#include <TH1F.h>

#include <vector> // test

using namespace sou;

ClassImp(ShowerLightSource)

//____________________________________________________________________________________________
   ShowerLightSource::ShowerLightSource() : LightSource("SHOWER"), EsafMsgSource(),
                                           fFluocalcul(0),fCrkcalcul(0), fShowerNemax(0) {
    //
    // ctor
    //
    fPh_in_atmo = new ListPhotonsInAtmosphere;
    if(!fPh_in_atmo) Msg(EsafMsg::Panic) << "Pb for memory allocation of fPh_in_atmo"<<MsgDispatch;
    Configure();
    Msg(EsafMsg::Info) << "ShowerLightSource built" << MsgDispatch;
}

//_____________________________________________________________________________________________
 ShowerLightSource::~ShowerLightSource() {
    //
    // dtor
    //
    SafeDelete(fPh_in_atmo);
    SafeDelete(fFluocalcul);
    SafeDelete(fCrkcalcul);
}

//______________________________________________________________________________________________
 void ShowerLightSource::Configure() {
    //
    // configure fluorescence and cerenkov calculators 
    //

    // Get fluorescence calculator
    string fluoname = Conf()->GetStr("ShowerLightSource.FluoCalculator");
    fFluocalcul = LightSourceFactory::Get()->GetFluoCalculator( fluoname );
    Msg(EsafMsg::Info) << "Fluo calculator name " <<fFluocalcul->GetName()<< MsgDispatch;

    // Get Cerenkov calculator
    string crkname = Conf()->GetStr("ShowerLightSource.CrkCalculator");
    fCrkcalcul = LightSourceFactory::Get()->GetCrkCalculator( crkname );
    Msg(EsafMsg::Info) << "Cerenkov calculator name " <<fCrkcalcul->GetName()<< MsgDispatch;

    // Get Energy Distribution type and name
    fEnergyDistributionType = Conf()->GetStr("ShowerLightSource.EnergyDistributionType");
    fEnergyDistributionName = Conf()->GetStr("ShowerLightSource.EnergyDistributionName");
    if(fEnergyDistributionType == "single") {
        if(!(fEnergyDistributionName == "'80MeV'" || fEnergyDistributionName == "'24MeV'"))
	    Msg(EsafMsg::Panic) << "<Configure> Wrong type of electrons SINGLE energy = "<<fEnergyDistributionName<< MsgDispatch;
    }
    else if(fEnergyDistributionType == "parametrized") {
        if(!(fEnergyDistributionName == "hillas" || fEnergyDistributionName == "giller" || fEnergyDistributionName == "nerling"))
	    Msg(EsafMsg::Panic) << "<Configure> Wrong type of PARAMETRIZED electrons energy distribution = "<<fEnergyDistributionName<< MsgDispatch;
    }

    // Get Lateral Distribution Name
    fLateralDistributionName = Conf()->GetStr("ShowerLightSource.LateralDistributionName");

    // Get Angular Distribution Name
    fAngularDistributionName = Conf()->GetStr("ShowerLightSource.AngularDistributionName");
    
    // to know whether electrons angular deviation used for Yield calculation
    string usedev = Conf()->GetStr("ShowerLightSource.fUseAngDev");
    if(usedev == "yes") {
        if(fAngularDistributionName == "NULL") {
	    Msg(EsafMsg::Warning) <<"<Configure()> fUseAngDev cannot be YES if no angular distribution -> has been set to NO "<<MsgDispatch;
	    fUseAngDev = false;
	}
	else fUseAngDev = true;
    }
    else if(usedev == "no") fUseAngDev = false;
    else Msg(EsafMsg::Panic) <<"<Configure()> Wrong argument fot fUseAngDev -> "<<usedev<<MsgDispatch;

}

//________________________________________________________________________________________
 void ShowerLightSource::Reset() {
    //
    // reset internal list of photons
    //
    if(fPh_in_atmo) fPh_in_atmo->Reset();
    if(fFluocalcul) fFluocalcul->Reset();
    if(fCrkcalcul) fCrkcalcul->Reset();
    fShowerNemax = 0;
}

//_________________________________________________________________________________________ 
 PhotonsInAtmosphere* ShowerLightSource::Get( const PhysicsData* data ) {
    //
    // generate photons in atmosphere from shower
    //
    
    Reset();

    // get track
    if ( data->Type() != "shower" )
        Msg(EsafMsg::Panic) << "Wrong PhysicsData in ShowerLightSource. Must be shower"<<MsgDispatch;
    ShowerTrack *track = (ShowerTrack*)data;
    
    UInt_t nbshowstep = track->GetNumStep();
    if(!nbshowstep) return (PhotonsInAtmosphere*)0;

    // energy spectrum of charged particles 
    const char* ED_char = fEnergyDistributionName.c_str();
    TString EDname;
    EDname.Append(ED_char);
#ifdef DEBUG
    Msg(EsafMsg::Debug)<< "Energy Distribution Name "<< EDname << MsgDispatch;
#endif
    TF2* EnergyDistribution = NULL;
    if(fEnergyDistributionType == "parametrized")
        EnergyDistribution = track->GetEnergyDistribution(EDname);

    
    // electron angular distribution  
    const char* AD_char = fAngularDistributionName.c_str();
    TString ADname;
    ADname.Append(AD_char);
#ifdef DEBUG
    Msg(EsafMsg::Debug)<< "Angular Distribution Name "<< ADname << MsgDispatch;
#endif
    TF2* AngularDistribution = NULL;
    if(!((fAngularDistributionName == "NULL")||(fAngularDistributionName == "histos")) ) 
        AngularDistribution = track->GetAngularDistribution(ADname);


    // electron lateral distribution  
    const char* LD_char = fLateralDistributionName.c_str();
    TString LDname;
    LDname.Append(LD_char);
#ifdef DEBUG
    Msg(EsafMsg::Debug)<< "Lateral Distribution Name "<< LDname << MsgDispatch;
#endif
    TF2* LateralDistribution = NULL;
    if(!((fLateralDistributionName == "NULL")||(fLateralDistributionName == "histos")) ) 
        LateralDistribution = track->GetLateralDistribution(LDname);

    // if BunchRadiativeTransfer is gonna to be used : 
    // initialization of the photons track -- assess the mean depth of ShowerSteps (assumed cst)
    ConfigFileParser* pConfig = Config::Get()->GetCF("LightToEuso","StandardLightToEuso");
    string RT = pConfig->GetStr("StandardLightToEuso.fRadiativeTransfer");
    UInt_t stepjump(0);
    if(RT == "bunch") {
	pConfig = Config::Get()->GetCF("RadiativeTransfer","BunchRadiativeTransfer");
	Double_t depthstep = pConfig->GetNum("BunchRadiativeTransfer.DepthStep")*g/cm2;
	Double_t meanshowdepthstep = Atmosphere::Get()->Grammage(track->FirstPos(),track->LastPos()) / nbshowstep;
	stepjump = (UInt_t) floor(depthstep/meanshowdepthstep);
	if(!stepjump) stepjump = 1;
	fPh_in_atmo->ClearTrack();
	if(!(nbshowstep % stepjump)) fPh_in_atmo->SetNbTrackSteps(nbshowstep/stepjump + 1);
	else fPh_in_atmo->SetNbTrackSteps(nbshowstep/stepjump + 2);
    }

    // Estimate the total number of electrons
    Float_t nTotal(0),nTracked(0); 
    for (UInt_t i=0;i<nbshowstep;i++) {
        const ShowerStep& s = (*track)[i];
        if (fShowerNemax < s.GetNelectrons()) fShowerNemax = s.GetNelectrons();
        if(RT == "bunch")
	    if(!(i % stepjump)) fPh_in_atmo->FillTrack(s.GetXYZi());
        nTotal += s.GetNelectrons();
    }
    if(RT == "bunch") fPh_in_atmo->FillTrack((*track)[nbshowstep-1].GetXYZf());

    // generate light
    Int_t progress = 0;

    for (UInt_t i=0;i<nbshowstep;i++) {
        const ShowerStep& s = (*track)[i];
	// generate fluorescence bunch
        BunchOfPhotons* bfluo = MakeFluoStep(s,EnergyDistribution,LateralDistribution,AngularDistribution);
        fPh_in_atmo->Add(bfluo);
        // generate cerenkov bunch
        BunchOfPhotons* bcer = MakeCerenkovStep(s,EnergyDistribution,LateralDistribution,AngularDistribution);
        fPh_in_atmo->Add(bcer);

        nTracked += s.GetNelectrons();

	if (100.*(i+1)/nbshowstep >= progress) {
	   Msg(EsafMsg::Info).SetProgress(progress);
	   Msg(EsafMsg::Info) << "Processing:" << MsgCount;
	   progress+=10;
	 }
    }

    return fPh_in_atmo;
   
}
//_________________________________________________________________________________________ 
 MCTruth* ShowerLightSource::Truth() {
    return NULL;
}

//_________________________________________________________________________________________ 
 BunchOfPhotons* ShowerLightSource::MakeFluoStep(const ShowerStep& step,TF2* EnergyDistribution,
                                                TF2* LateralDistribution, TF2* AngularDistribution) {
    //
    // generate a fluorescence bunch of photons for this shower step 
    //

    //mean position in space 
    EarthVector posi = step.GetXYZi();
    EarthVector posf = step.GetXYZf();
#ifdef DEBUG
    //Msg(EsafMsg::Debug)<<" X "<< posi.X()/km <<" Y "<<posi.Y()/km<<" Z "<<posi.Z()/km<<MsgDispatch;
#endif
    EarthVector pos=0.5*(posi+posf);


    // fluorescence yield for the energy spectrum of charged particles
    EsafSpectrum spectrum(357*nm);  //initialisation
    Double_t age = 0.5*(step.GetAgef()+step.GetAgei());
    Double_t TotalYield = 0;

    if( step.GetNelectrons() > 0.001*fShowerNemax) {

	if ( fEnergyDistributionType == "histos" ) {        
            const TH1F* ehisto = step.GetHistoEnergy();
            TotalYield = fFluocalcul->GetFluoYieldHisto(pos.Zv(),ehisto,&spectrum);  
    #ifdef DEBUG
	   // Msg(EsafMsg::Debug)<<"integration over TH1F - Yield Fluo " << TotalYield*m << MsgDispatch;
    #endif
	}
	else if ( EnergyDistribution && fEnergyDistributionType == "parametrized") {
            TF12 EnergyDistributionStep("EDS_name",EnergyDistribution,age,"X");
            TotalYield = fFluocalcul->GetFluoYield(pos.Zv(),&EnergyDistributionStep,&spectrum);
    #ifdef DEBUG
	   // Msg(EsafMsg::Debug)<<"integration over TF12 - Yield Fluo " << TotalYield*m << MsgDispatch;
    #endif
	}
	else if(fEnergyDistributionType == "single") {
            if(fEnergyDistributionName == "'80MeV'") TotalYield = fFluocalcul->GetFluoYield(pos.Zv(),80*MeV,&spectrum);
            else if(fEnergyDistributionName == "'24MeV'") TotalYield = fFluocalcul->GetFluoYield(pos.Zv(),24*MeV,&spectrum);
    #ifdef DEBUG
	   // Msg(EsafMsg::Debug)<<" E= 80 MeV - Yield Fluo " << TotalYield*m << MsgDispatch;
    #endif
	}
	else Msg(EsafMsg::Panic)<<"<MakeFluoStep> Wrong type of EnergyDistribution = "<<fEnergyDistributionType<< MsgDispatch;
    }
    else TotalYield = 0.;

    // projection of angular distribution along X (fixed value of Y)
    // energy threshold fixed to minimum value to get the larger spectrum
    // presently REMOVE because of electrons Energy cut issue
    //TF12* AngularDistributionStep = 0;
    //if(AngularDistribution) AngularDistributionStep = new TF12("ADS_name",AngularDistribution,0.5,"X");

    // total weight of the bunch of photons
    Double_t dl = (step.GetXYZf()-step.GetXYZi()).Mag(); 
    //if(fUseAngDev) dl /= cos(AngularDistributionStep->Mean(AngularDistributionStep->GetXmin(),AngularDistributionStep->GetXmax()));
    Double_t nbPhotonsInBunch = step.GetNelectrons() * TotalYield * dl;
#ifdef DEBUG
   // Msg(EsafMsg::Debug)<< "nelec "<<step.GetNelectrons()<<" nbpib fluo "<<nbPhotonsInBunch << " dl= " << dl/m << MsgDispatch;
#endif
    // Create a new bunch of Photons with its ParentBunch
    BunchOfPhotons* b = new BunchOfPhotons(nbPhotonsInBunch,TotalYield,posi,posf,step.GetTimei(),
                                           step.GetTimef(),spectrum,step.GetParentTrack()->GetDirVers(),Fluo);

   
    // Set the lateral distribution to the bunch if it exists
    if ( fLateralDistributionName == "histos" ) {
        const TH1F* lhisto = step.GetHistoLateral();
        if (lhisto) b->SetParentLateral(*lhisto);
    }
    else if ( LateralDistribution ) {
        if ( LateralDistribution->GetNpar() == 1 ) {
            // gives moliere radius as parameter in meter
            // Get Atmosphere for Density
            const Atmosphere* atmo = Atmosphere::Get();
            Double_t Rm = 9.6 * gram/cm2 / atmo->Air_Density( pos.Zv() ) / m; // in meters
            LateralDistribution->SetParameters(Rm,0);
        }
        TF12 LateralDistributionStep("LDS_name",LateralDistribution,age,"X");
        b->SetParentLateral(LateralDistributionStep);
    }
    
    //SafeDelete(AngularDistributionStep);
    
    return b;
}


//_________________________________________________________________________________________ 
 BunchOfPhotons* ShowerLightSource::MakeCerenkovStep(const ShowerStep& step, TF2* EnergyDistribution,
                                                    TF2* LateralDistribution, TF2* AngularDistribution) {
    //
    // Generate a Cerenkov Bunch Of Photons for this shower step
    //

    //mean position in space 
    EarthVector posi = step.GetXYZi();
    EarthVector posf = step.GetXYZf();
    EarthVector pos=0.5*(posi+posf);

    // Cerenkov yield for the energy spectrum of charged particles
    Double_t age = 0.5*(step.GetAgef()+step.GetAgei());
    Double_t TotalYield;
    
    if( step.GetNelectrons() > 0.001*fShowerNemax) {
       if ( fEnergyDistributionType == "histos" ) {
           const TH1F* ehisto = step.GetHistoEnergy();
           TotalYield = fCrkcalcul->GetCrkYield(pos.Zv(),ehisto);
   #ifdef DEBUG
	//   Msg(EsafMsg::Debug)<<"integration over TH1F - Yield Cerenkov " << TotalYield*m << MsgDispatch; 
   #endif
       }
       else if ( EnergyDistribution && fEnergyDistributionType == "parametrized" ) {
           TF12 EnergyDistributionStep("EDS_name",EnergyDistribution,age,"X");
           TotalYield = fCrkcalcul->GetCrkYield(pos.Zv(),&EnergyDistributionStep);
   #ifdef DEBUG
	 //  Msg(EsafMsg::Debug)<<"integration over TF12 - Yield Cerenkov " << TotalYield*m << MsgDispatch;
   #endif
       }
       else if(fEnergyDistributionType == "single") {
           if(fEnergyDistributionName == "'80MeV'") TotalYield = fCrkcalcul->GetCrkYield(pos.Zv(),80*MeV);
           else if(fEnergyDistributionName == "'24MeV'") TotalYield = fCrkcalcul->GetCrkYield(pos.Zv(),24*MeV,false); // false to remove energy threshold for ckov prod
   #ifdef DEBUG
	//   Msg(EsafMsg::Debug)<<"E=80MeV - Yield Cerenkov " << TotalYield*m << MsgDispatch;
   #endif
       }
       else Msg(EsafMsg::Panic)<<"<MakeCerenkovStep> Wrong type of EnergyDistribution = "<<fEnergyDistributionType<< MsgDispatch;
    }
    else TotalYield = 0;

    // Cerenkov wavelenght spectrum 
    EsafSpectrum *spectrum = fCrkcalcul->GetCrkSpectrum();

    // projection  angular distribution ofalong X (fixed value of Y)
    TF12* AngularDistributionStep = 0;
    if(AngularDistribution) AngularDistributionStep = new TF12("ADS_name",AngularDistribution,fCrkcalcul->GetEnergyThreshold(pos.Zv())/MeV,"X");

    // weight of the bunch of photons
    Double_t dl = (step.GetXYZf() - step.GetXYZi()).Mag();
    if(fUseAngDev) dl /= cos(AngularDistributionStep->Mean(AngularDistributionStep->GetXmin(), AngularDistributionStep->GetXmax()));
    Double_t nelec = step.GetNelectrons();
    Double_t nbPhotonsInBunch = nelec*TotalYield*dl;
#ifdef DEBUG
   // Msg(EsafMsg::Debug)<< "nelec "<<nelec<<" nbpib cer "<<nbPhotonsInBunch << " dl= " << dl/m << MsgDispatch;
#endif

    // Create a new bunch of Photons with its ParentBunch
    BunchOfPhotons* b = new BunchOfPhotons(nbPhotonsInBunch,TotalYield,posi,posf,step.GetTimei(),
                        step.GetTimef(),*spectrum,step.GetParentTrack()->GetDirVers(),Cerenkov);


    // Set the lateral distribution to the parent bunch
    if ( fLateralDistributionName == "histos" ) {
        const TH1F* lhisto = step.GetHistoLateral();
        if (lhisto) b->SetParentLateral(*lhisto);
        Msg(EsafMsg::Warning)<<"Histogram for Angular distribution not handled here, it SHOULD"<< MsgDispatch;	
    }
    else if ( LateralDistribution ) {
        if ( LateralDistribution->GetNpar() == 1 ) {
            // gives moliere radius as parameter in meter
            // Get Atmosphere for Density
            const Atmosphere* atmo = Atmosphere::Get();
            Double_t Rm = 9.6 * gram/cm2 / atmo->Air_Density( pos.Zv() ) / m; // in meters
            LateralDistribution->SetParameters(Rm,0);
        }
        TF12 LateralDistributionStep("LDS_name",LateralDistribution,age,"X");
        b->SetParentLateral(LateralDistributionStep);
    }
    
    // Set the angular distribution to the parent bunch
    if ( AngularDistribution ) b->SetParentAngular(*AngularDistributionStep);
    
    SafeDelete(spectrum);
    SafeDelete(AngularDistributionStep);
    
    return b;
}


About Us | EUSO Official Website | Web pages created by Roberto Pesce and Alessandro Thea - Last Update Wed Nov 16 16:57:39 2005 Wed Nov 16 16:29:22 2005