// 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;
}