// ESAF : Euso Simulation and Analysis Framework
// $Id: ShowerLightSource.cc,v 1.53 2005/05/10 16:54:10 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
ClassImp(ShowerLightSource)
//____________________________________________________________________________________________
ShowerLightSource::ShowerLightSource() : LightSource("SHOWER"), EsafMsgSource(),
fFluocalcul(0),fCrkcalcul(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 Name
fEnergyDistributionName = Conf()->GetStr("ShowerLightSource.EnergyDistributionName");
// Get Lateral Distribution Name
fLateralDistributionName = Conf()->GetStr("ShowerLightSource.LateralDistributionName");
// Get Angular Distribution Name
fAngularDistributionName = Conf()->GetStr("ShowerLightSource.AngularDistributionName");
}
//________________________________________________________________________________________
void ShowerLightSource::Reset() {
//
// reset internal list of photons
//
if(fPh_in_atmo) fPh_in_atmo->Reset();
if(fFluocalcul) fFluocalcul->Reset();
if(fCrkcalcul) fCrkcalcul->Reset();
}
//_________________________________________________________________________________________
PhotonsInAtmosphere* ShowerLightSource::Get( const PhysicsData* data ) {
//
// generate photons in atmosphere from shower
//
Reset();
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 == 0) 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( !((fEnergyDistributionName == "NULL")||(fEnergyDistributionName == "histos")) )
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);
// initialization of the photons track -- assess the mean depth of ShowerSteps (assumed cst)
ConfigFileParser* 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;
UInt_t stepjump = (UInt_t) floor(depthstep/meanshowdepthstep);
if(stepjump == 0) stepjump = 1;
fPh_in_atmo->ClearTrack();
if(nbshowstep % stepjump == 0) 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(i % stepjump == 0) fPh_in_atmo->FillTrack(s.GetXYZi());
nTotal += s.GetNelectrons();
}
fPh_in_atmo->FillTrack((*track)[nbshowstep-1].GetXYZf());
Bool_t next(0),subnext(0);
Float_t left(0),right(10),subleft(0),subright(2.5);
Msg(EsafMsg::Info) << "Processing: 0%" << MsgFlush;
for (UInt_t i=0;i<nbshowstep;i++) {
const ShowerStep& s = (*track)[i];
// generate fluorescence bunch
BunchOfPhotons* bfluo = MakeFluoStep(s,EnergyDistribution,LateralDistribution);
fPh_in_atmo->Add(bfluo);
// generate cerenkov bunch
BunchOfPhotons* bcer = MakeCerenkovStep(s,EnergyDistribution,LateralDistribution,AngularDistribution);
fPh_in_atmo->Add(bcer);
nTracked += s.GetNelectrons();
Float_t fraction = 100*i/nbshowstep;
if (left<fraction && fraction <=right)
next = kFALSE;
else if (fraction>right) {
next = kTRUE;
left = right;
right += 10;
}
if (subleft<fraction && fraction <=subright) {
subnext = kFALSE;
}
if (fraction > subright) {
subnext = kTRUE;
subleft = subright;
subright +=2;
}
if (next) Msg(EsafMsg::Info) << left << "%" << MsgFlush;
if (subnext) Msg(EsafMsg::Info)<< "." << MsgFlush;
}
Msg(EsafMsg::Info) << " -done. " << MsgDispatch;
return fPh_in_atmo;
}
//_________________________________________________________________________________________
MCTruth* ShowerLightSource::Truth() {
return NULL;
}
//_________________________________________________________________________________________
BunchOfPhotons* ShowerLightSource::MakeFluoStep(const ShowerStep& step,TF2* EnergyDistribution,
TF2* LateralDistribution) {
//
// 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 ( fEnergyDistributionName == "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 ) {
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 {
TotalYield = fFluocalcul->GetFluoYield(pos.Zv(),80*MeV,&spectrum);
#ifdef DEBUG
// Msg(EsafMsg::Debug)<<" E= 80 MeV - Yield Fluo " << TotalYield*m << MsgDispatch;
#endif
}
// total weight of the bunch of photons
Double_t dl = (step.GetXYZf()-step.GetXYZi()).Mag();
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);
}
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 ( fEnergyDistributionName == "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 ) {
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 {
TotalYield = fCrkcalcul->GetCrkYield(pos.Zv(),80*MeV);
#ifdef DEBUG
// Msg(EsafMsg::Debug)<<"E=80MeV - Yield Cerenkov " << TotalYield*m << MsgDispatch;
#endif
}
// Cerenkov wavelenght spectrum
EsafSpectrum *spectrum = fCrkcalcul->GetCrkSpectrum();
// weight of the bunch of photons
Double_t dl = (step.GetXYZf()-step.GetXYZi()).Mag();
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);
}
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 ) {
TF12 AngularDistributionStep("ADS_name",AngularDistribution,
fCrkcalcul->GetEnergyThreshold(pos.Zv())/MeV,"X");
b->SetParentAngular(AngularDistributionStep);
}
SafeDelete(spectrum);
return b;
}
