// ESAF : Euso Simulation and Analysis Framework
// $Id: SlastShowerGenerator.cc,v 1.53 2005/05/12 09:15:11 stutz Exp $
// Author: Dmitry V.Naumov 03/03/2004
/*****************************************************************************
* ESAF: Euso Simulation and Analysis Framework *
* *
* Id: SlastShowerGenerator *
* Package: Shower *
* Coordinator: Sergio.Bottai, Dmitry Naumov *
* *
*****************************************************************************/
//_____________________________________________________________________________
//
// SlastShowerGenerator
//
// Slast is short of "Shower Light Attenuated to the Space Telescope".
// Originally written by Dmitry V.Naumov as a collection of F77 routines doing
// the full chain:
// Shower generation, light production (fluorescence and cherenkov), light
// attenuation in the atmosphere.
// See SlastLightToEuso C++ / Fortran interface which is doing the full chain.
//
// SlastShowerGenerator is the simple generator written in pure C++ which is
// NOT doing the full simulation chain. Instead it returnes only the
// ShowerTrack object which is taken later by ShowerLightSource and
// RadiativeTransfer parts of the ESAF package.
//
// ESAF prodides some testing macros which can help to understand how to use
// SlastShowerGenerator. See macros/CheckDistrs.C as an example.
//
// A simple way to use SlastShowerGenerator is directly in root:
//
// root[] gSystem->Load("libatmosphere.so");
// root[] gSystem->Load("libgenbase.so");
// root[] gSystem->Load("libshowers.so");
// root[] SlastShowerGenerator *ShowerGenerator = new SlastShowerGenerator;
// root[] ShowerTrack *track = ShowerGenerator->Get();
// root[] track->DrawXYZ()
//
// as an example.
//
// Now having the track object in hand one can use whatever functions avalaible
// for it. See ShowerTrack class for documentation of ShowerTrack object
#include "SlastShowerGenerator.hh"
#include "MCTruth.hh"
#include "ShowerTrack.hh"
#include "Atmosphere.hh"
#include "Config.hh"
#include "EConst.hh"
#include <TMath.h>
#include <TRandom2.h>
#include <TProfile.h>
ClassImp(SlastShowerGenerator)
using EConst::EarthRadius;
using namespace TMath;
//_________________________________________________________________________________________
SlastShowerGenerator::SlastShowerGenerator() : EventGenerator("slast++"), EsafMsgSource(),
fTrack(0), fTruth(0), fInCome(0), fOutCome(0), fQuiet(kFALSE) {
// This is ctor
if(!Init()) Msg(EsafMsg::Panic)<<"SlastShowerGenerator can not be initialized"<<MsgDispatch;
}
//_________________________________________________________________________________________
Bool_t SlastShowerGenerator::Init()
{
// SlastShowerGenerator is initialized reading from config/LightToEuso/GeneratorLightToEuso.cfg file
// Edit that file to setup desired parameters of the showers such as energy, theta, phi intervals
// or initial position in case of unique parameters.
// Init() is called in SlastShowerGenerator ctor
Msg(EsafMsg::Info) << "SlastShowerGenerator Enabled" << MsgDispatch;
// Enabling the atmosphere
Atmosphere::Get()->GetType();
Msg(EsafMsg::Info) << "SlastShowerGenerator Enables Atmosphere " << Atmosphere::Get()->GetType() << MsgDispatch;
// EUSO detector initializations
ConfigFileParser *pConfig = Config::Get()->GetCF("General","Euso");
fEusoAltitude= pConfig->GetNum("Euso.fAltitude")*km;
fEusoVector.SetXYZ(0,0,fEusoAltitude);
// General Shower Generator initializations
pConfig = Config::Get()->GetCF("LightToEuso","GeneratorLightToEuso");
fFoV = pConfig->GetNum("GeneratorLightToEuso.FoV")*deg;
fEnergyMin = pConfig->GetNum("GeneratorLightToEuso.EnergyRangeMin");
fEnergyMax = pConfig->GetNum("GeneratorLightToEuso.EnergyRangeMax");
fEnergySlope = pConfig->GetNum("GeneratorLightToEuso.EnergySlope");
fThetaMin = pConfig->GetNum("GeneratorLightToEuso.ThetaRangeMin")*deg;
fThetaMax = pConfig->GetNum("GeneratorLightToEuso.ThetaRangeMax")*deg;
fPhiMin = pConfig->GetNum("GeneratorLightToEuso.PhiRangeMin")*deg;
fPhiMax = pConfig->GetNum("GeneratorLightToEuso.PhiRangeMax")*deg;
fType = pConfig->GetStr("GeneratorLightToEuso.UhecrType");
fDepthStep = pConfig->GetNum("GeneratorLightToEuso.DepthStep")*g/cm2;
fFirstPointX = pConfig->GetNum("GeneratorLightToEuso.InteractionVectorX")*km;
fFirstPointY = pConfig->GetNum("GeneratorLightToEuso.InteractionVectorY")*km;
fFirstPointZ = pConfig->GetNum("GeneratorLightToEuso.InteractionVectorZ")*km;
fFirstType = pConfig->GetStr("GeneratorLightToEuso.InteractionType");
if (!fQuiet) {
MsgForm(EsafMsg::Info,"Generate %s event(s) in:",fType.c_str());
MsgForm(EsafMsg::Info,"Energy interval [%.4E,%.4E]",fEnergyMin,fEnergyMax);
MsgForm(EsafMsg::Info,"Theta interval [%.4f,%.4f]",fThetaMin,fThetaMax);
MsgForm(EsafMsg::Info,"Phi interval [%.4f,%.4f]",fPhiMin,fPhiMax);
}
fAtomicMass = EConst::AtomicMass(fType);
// Initialize SLAST specific options
pConfig = Config::Get()->GetCF("LightToEuso","SlastLightToEuso");
fEnergyDistribution = pConfig->GetStr("SlastLightToEuso.EnergyDistributionParametrization");
fShowerParametrization = pConfig->GetStr("SlastLightToEuso.ShowerParametrization");
return kTRUE;
}
//_________________________________________________________________________________________
SlastShowerGenerator::~SlastShowerGenerator() {
// This is dtor
SafeDelete(fTrack);
SafeDelete(fTruth);
}
//_________________________________________________________________________________________
void SlastShowerGenerator::SetShower(Double_t e,Double_t t,Double_t p,EarthVector pos) {
// This is a special method to setup desired showertrack parameters from the code (not from the config file)
// This can be a usefull function in the reconstruction code
SetEnergy(fEnergyMin = fEnergyMax = e);
SetTheta(fThetaMin = fThetaMax = t);
SetPhi(fPhiMin = fPhiMax = p);
fFirstPointX = pos.X();
fFirstPointY = pos.Y();
fFirstPointZ = pos.Z();
}
//_________________________________________________________________________________________
void SlastShowerGenerator::GetX1() {
// Returns atmosphere depth before the first interaction.
// Attention: units are ESAF internal. In order to get in human way one has to multiple by cm2/g
if ( fFirstType == "POS" ) {
fX1 = Atmosphere::Get()->Grammage(fInitPoint,-fOmega,"dir");
// cout << " fInitPoint " << fInitPoint.X() <<" "<<fInitPoint.Y()<<" "<<fInitPoint.Z()<<endl;
// cout << " fX1 POS " << fX1*cm2/g << endl;;
}
else if (fFirstType == "X1" ) {
ConfigFileParser *pConfig = Config::Get()->GetCF("LightToEuso","GeneratorLightToEuso");
fX1 = pConfig->GetNum("GeneratorLightToEuso.InteractionX1")*g/cm2;
// cout << " fX1 X1 " << fX1*cm2/g << endl;;
}
else
fX1 = -Log(gRandom->Rndm())*HadronInteractionLength(fAtomicMass,fEnergy)*g/cm2;
// cout <<" fX1 std " << fX1*cm2/g<<endl;
}
//_________________________________________________________________________________________
Double_t SlastShowerGenerator::HadronInteractionLength(Double_t A, Double_t E) {
// SlastShowerGenerator internal function which computes the hadronic interactio length.
// Used by GetX1() method
Float_t R0 = 1.287e-13, Aair = 14.483, q = 0.93, AN = 14.00674, AO = 15.9994;
Float_t AA = 39.948, AC = 12.011, CN = 0.7847, CO = 0.2105, CA = 0.0047, CC = 0.0001;
Float_t S_A = Pi()*Power(R0,2);
Float_t F_A1 = Power(S_A*Na(),-1);
Float_t F_A2 = Aair*F_A1;
Float_t X = Power(A,1./3);
Float_t Y = 1./X;
Float_t B = Power(AN,1./3);
Float_t S = CN*Power(B+X-q*(1./B+Y),2); // Nitrogen (N2);
B = Power(AO,1./3);
S += CO*Power(B+X-q*(1./B+Y),2); // Oxigen (O2)
B = Power(AA,1./3);
S += CA*Power(B+X-q*(1./B+Y),2); // Argon (Ar)
B = Power(AC,1./3);
S += CC*Power(B+X-q*(1./B+Y),2); // Carbon (in CO2)
return F_A2/S*NucleonAirCrossSection(1.e11/A)/NucleonAirCrossSection(E/A);
}
//_________________________________________________________________________________________
Double_t SlastShowerGenerator::NucleonAirCrossSection(Double_t E) {
/*
Computes NUCLEN + AIR CROSS SECTION:
... USES AN EMPIRICAL PARAMETRIZATION OF TOTAL INELASTIC
... NUCLEON AIR CROSS SECTION (in mbarn)
... INPUT: NUCLEON ENERGY in eV
... OUTPUT: NUCLEON AIR CROSS SECTION (in mbarn)
... REFERENCES:
... 1. Mielke H.H., Foller, M, Knapp J.
... // J.Phys.G: Nucl.Part.Phys. 1994. V 20, P637
... 2. V.A. Naumov, T.S. Sinegovskaya (eq (17))
*/
return 290 - 8.7*Log(E/1.e9) + 1.14*Power(Log(E/1.e9),2);
}
//_________________________________________________________________________________________
const ShowerStep& SlastShowerGenerator::GetShowerStep() {
// Filling the ShowerStep object with relevant information
fStep.Clear();
fStep.fXi = fXcurrent;
fStep.fXf = fXnext;
fStep.fXYZi = fCurrentPoint;
fStep.fXYZf = fNextPoint;
fStep.fTimei = fTimecurrent;
fStep.fTimef = fTimenext;
GetShowerParametrization(fXcurrent);
fStep.fAgei = fAge;
GetShowerParametrization(fXnext);
fStep.fAgef = fAge;
GetShowerParametrization((fXnext+fXcurrent)/2);
fStep.fNelectrons = fNe;
return fStep;
}
//_________________________________________________________________________________________
void SlastShowerGenerator::GetShowerParametrization(Double_t x) {
// Assigning of fNe - number of electrons in the current step according to the Shower Parameterization.
// GIL stands for Greizen-Ilina-Linsley which is the QGSJET model
// parameterization.
// GFA stands for Gaussian Function in Age
// GHF stands for Gaisser Hillas Function
if (fShowerParametrization == "GIL"){
GIL(x);}
else if (fShowerParametrization == "GFA"){
GFA(x);}
else if (fShowerParametrization == "GHF"){
GHF(x);}
else
throw invalid_argument("SlastShowerGenerator does not know about ShowerParametrization "+
fShowerParametrization);
}
//_________________________________________________________________________________________
Bool_t SlastShowerGenerator::GetFarFirstPoint() {
// This method generates the first ``far'' point which is ensured to be visible by EUSO FoV.
// 1. Generate point uniformly on the Earth surface underneath of EUSO
Double_t x,y,z, Zmin, H = fEusoVector.Z() + EarthRadius();
Zmin = H*Power(Sin(fFoV),2) +
Cos(fFoV)*Sqrt( Power(EarthRadius(),2) -
Power(H,2)*Power(sin(fFoV),2) );
while (kTRUE) {
gRandom->Sphere(x,y,z,EConst::EarthRadius());
if (z > Zmin) break;
}
fInitPoint.SetXYZ(x,y,z);
fEarthImpact = fInitPoint;
fEarthImpact.SetZ(fInitPoint.Z()-EarthRadius());
// 2. Generate random Theta and Phi
fPhi = (fPhiMax - fPhiMin)*gRandom->Rndm() + fPhiMin;
Double_t r1 = Cos(2*fThetaMin);
Double_t r2 = Cos(2*fThetaMax);
fTheta = 0.5*ACos(r1 - gRandom->Rndm()*(r1-r2));
// fOmega.SetXYZ(-Sin(fTheta)*Cos(fPhi),-Sin(fTheta)*Sin(fPhi),-Cos(fTheta));
EarthVector v1 = fInitPoint.Unit();
EarthVector vrot;
fOmega.SetXYZ(1,0,0);
EarthVector Uz(0,0,1);
vrot = Uz.Cross(v1);
fOmega.Rotate(v1.Theta(),vrot);
fOmega.Rotate(fPhi,v1);
vrot = v1.Cross(fOmega);
fOmega.Rotate(TMath::Pi()/2-fTheta,vrot);
cout << "fOmega1 " <<-Sin(fTheta)*Cos(fPhi)<<" "<<-Sin(fTheta)*Sin(fPhi)<<" "<<-Cos(fTheta)<<endl;
cout << "fOmega2 " <<fOmega.X()<<" "<<fOmega.Y()<<" "<<fOmega.Z()<<endl;
// 3. Jump to the higher altitude along the generated direction to the sphere of Radius rEarth + h
Float_t h = 80*km;
// 4. Find first need track length
Double_t L = fOmega*fInitPoint + Sqrt(h*(2*fInitPoint.Mag() + h) + Power(-fOmega*fInitPoint,2));
fInitPoint -= fOmega*L;
fInitPoint.SetZ(fInitPoint.Z()-EarthRadius());
fHitGround = 1;
return kTRUE;
}
//_________________________________________________________________________________________
void SlastShowerGenerator::GetFirstPoint(){
// Once ``far'' first point is found this method find the real first point according to the generated X1
Double_t L = 0.1*km, depth(0), h(0);
depth = Atmosphere::Get()->Grammage(fInitPoint,-fOmega,"dir");
while(true) {
if ( fFirstType == "POS" ) break;
fInitPoint += fOmega*L;
if (!fInitPoint.IsUnderSeaLevel())
h = fInitPoint.Zv();
else
break;
depth += Atmosphere::Get()->Air_Density(h)*L;
if(depth>=fX1) break;
}
fXcurrent = depth-fX1;
fXnext = fXcurrent;
fCurrentPoint = fInitPoint;
fNextPoint = fCurrentPoint;
fTimecurrent = 0;
fTimenext = 0;
}
//_________________________________________________________________________________________
Bool_t SlastShowerGenerator::GetNextStep() {
// This method checks if next point is visible, not outside of FoV and makes the next step if
// it is OK.
if ((fEusoVector-fCurrentPoint).Theta()<=fFoV)
fInCome = kTRUE;
if (fInCome && ((fEusoVector-fCurrentPoint).Theta()>fFoV || fCurrentPoint.IsUnderSeaLevel()) )
fOutCome = kTRUE;
if (fOutCome) {
Msg(EsafMsg::Debug) << " fOutCome true " << MsgDispatch;
return kFALSE;
}
if (fCurrentPoint.IsUnderSeaLevel()) return kFALSE;
fCurrentPoint = fNextPoint;
fNextPoint = fCurrentPoint;
Double_t L, h(0), depth(0);
Int_t cycling(0);
while(1) {
cycling++;
L = fDepthStep/Atmosphere::Get()->Air_Density(h)/10;
fNextPoint += fOmega*L;
if (!fNextPoint.IsUnderSeaLevel())
h = fNextPoint.Zv();
else
return kFALSE;
depth += Atmosphere::Get()->Air_Density(h)*L;
if (depth>=fDepthStep) break;
if (cycling>100000) {
MsgForm(EsafMsg::Debug,"NextPoint not found : cycling > 10000 with depth %f ",depth/g*cm2);
return kFALSE;
}
}
if (h<0){
Msg(EsafMsg::Debug) << " Shower reaches the ground " << MsgDispatch;
return kFALSE;
}
fXcurrent = fXnext;
fXnext = fXcurrent + depth;
fTimecurrent = (fCurrentPoint-fInitPoint).Mag()/EConst::C();
fTimenext = (fNextPoint-fInitPoint).Mag()/EConst::C();
return kTRUE;
}
//_________________________________________________________________________________________
void SlastShowerGenerator::GIL(Double_t x)
{
// Assigning of fNe - number of electrons in the current step according to GIL parameterization.
// (GIL stands for Greizen-Ilina-Linsley which is the QGSJET model parameterization.
Float_t t0 = 37.15;
Float_t t = x*cm2/g/t0;
Float_t t_max = 1.7 + 0.76*(Log(fEnergy/0.81e8) - Log(fAtomicMass));
fAge = 2*t/(t+t_max);
fNe = 0;
if(fAge<=0) return;
fNe = fEnergy/1.45e9*Exp(t-t_max-2*t*Log(fAge));
}
//_________________________________________________________________________________________
void SlastShowerGenerator::GFA(Double_t x)
{
// Assigning of fNe - number of electrons in the current step according to Gaussian Function in Age parameterization.
// Formulae from C. Song, Astropart. Physics 22(2004)151
// Xmax, Nmax, sigma values fitted and extrapolated
Float_t Xmax,Nmax;
Float_t sigma,p0,p1;
if (fAtomicMass==1) {
sigma = 0.3206 - 0.006597*Log10(fEnergy);
p0 = -8.9475;
p1 = 0.9874;
}
else {
sigma = 0.470-0.0135*Log10(fEnergy);
p0 = -9.043;
p1 = 0.9923;
}
// Xmax=(-3426.*sigma+1410.)*g/cm2;
Xmax = ((200*Log10(fEnergy)-7026)*sigma+1410.)*g/cm2;
Nmax = pow(10,p1*Log10(fEnergy)+p0);
fAge = 3*x/(x+2.*Xmax);
if(fAge<=0) return;
fNe = Nmax*Exp((-1./(2.*sigma*sigma))*Power(fAge-1.,2));
return;
}
//_________________________________________________________________________________________
void SlastShowerGenerator::GHF(Double_t x)
{
// Assigning of fNe - number of electrons in the current step according to
// Gaisser Hillas parametrisation
// Formulae from C. Song, Astropart. Physics 22(2004)151
// Xmax, Nmax, lambda, X1, X0 values fitted and extrapolated
Float_t Xmax,Nmax;
Float_t X,X1,X0,lambda;
if (fAtomicMass==1) {
X1 = 106.2 - 3.15*Log10(fEnergy);
X0 = 401.9 - 25.79*Log10(fEnergy);
Nmax = Exp(-20.95 + 2.292*Log10(fEnergy));
Xmax = -258.1 + 54.81*Log10(fEnergy);
lambda= 98.6 -1.922*Log10(fEnergy);
}
else {
X1 = 29.72 - 1.031*Log10(fEnergy);
X0 = 471.1 - 29.49*Log10(fEnergy);
Nmax = Exp(-21.68 + 2.329*Log10(fEnergy));
Xmax = -430.7 + 59.08*Log10(fEnergy);
lambda= 182.4 -6.1*Log10(fEnergy);
}
Xmax = Xmax*g/cm2;
X1 = X1*g/cm2;
X0 = X0*g/cm2;
// X=x+fX1;
X=x;
lambda = lambda*g/cm2;
fNe = Nmax * Power((X-X0)/(Xmax-X0),(Xmax-X0)/lambda) * Exp((Xmax-X)/lambda);
fAge = 3*x/(x+2.*Xmax); //should be defined ?
return;
}
//_________________________________________________________________________________________
PhysicsData* SlastShowerGenerator::Get() {
// method returning the ShowerTrack object
// generate random Energy, Theta, Phi and the interaction point
// energy slope is for differential spectrum
if ( fEnergySlope == -1.) { // flat differential spectrum in log scale
fEnergy = fEnergyMin * Exp(gRandom->Rndm() * Log(fEnergyMax/fEnergyMin) );
}
else {
Double_t e1 = Power(fEnergyMin,fEnergySlope + 1);
Double_t e2 = Power(fEnergyMax,fEnergySlope + 1);
fEnergy = Power( e1 + (e2-e1)*gRandom->Rndm(),1/(fEnergySlope + 1) );
}
// Check if only a specific point should be generated
if( fFirstType == "POS") {
fInitPoint.SetXYZ(fFirstPointX,fFirstPointY,fFirstPointZ);
fPhi = (fPhiMax - fPhiMin)*gRandom->Rndm() + fPhiMin;
Double_t r1 = Cos(2*fThetaMin);
Double_t r2 = Cos(2*fThetaMax);
fTheta = 0.5*ACos(r1 - gRandom->Rndm()*(r1-r2));
fOmega.SetXYZ(-Sin(fTheta)*Cos(fPhi),-Sin(fTheta)*Sin(fPhi),-Cos(fTheta));
DevelopShower();
//if (DevelopShower())
// fTrack->fInitPos.SetXYZ(fInitPoint.X(),fInitPoint.Y(),fInitPoint.Z()); // Reset Initial Position to required
}
else {
Int_t cycle(0);
while (kTRUE) {
if ( cycle++>1000 ) break;
if ( GetFarFirstPoint() && DevelopShower()) break;
}
}
if (!fQuiet) {
MsgForm(EsafMsg::Info,"SlastShowerGenerator produced a ShowerTrack: Energy = %.2E MeV Theta = %.2f deg Phi = %.2f deg",
fTrack->GetEnergy(), fTrack->GetTheta()*RadToDeg(),fTrack->GetPhi()*RadToDeg());
MsgForm(EsafMsg::Info," originated at: (%.2f,%.2f,%.2f) km with %d steps",
fTrack->GetInitPos().X()/km,fTrack->GetInitPos().Y()/km,fTrack->GetInitPos().Z()/km, fTrack->GetNumStep());
}
return fTrack;
}
//_________________________________________________________________________________________
void SlastShowerGenerator::Reset() {
// Reset ShowerTrack and SlastShowerGenerator internal parameters
fHitGround = 0;
fInCome = kFALSE;
fOutCome = kFALSE;
fInitPoint.SetXYZ(0,0,0);
fOmega.SetXYZ(0,0,0);
fCurrentPoint.SetXYZ(0,0,0);
fNextPoint.SetXYZ(0,0,0);
fEarthImpact.SetXYZ(0,0,0);
}
//_________________________________________________________________________________________
void SlastShowerGenerator::GetShowerInfo() {
// Fills ShowerTrack parameters
fTrack->fEnergy = fEnergy*eV;
fTrack->fTheta = fTheta;
fTrack->fPhi = fPhi;
fTrack->fX1 = fX1;
fTrack->fEthrEl = 0;
fTrack->fDirVers = fOmega;
fTrack->fInitPos = fInitPoint;
fTrack->fHitGround = fHitGround;
fTrack->fEarthImpact = fEarthImpact;
}
//_________________________________________________________________________________________
Bool_t SlastShowerGenerator::DevelopShower() {
// Method developing the shower
GetX1();
GetFirstPoint();
// Develop a shower of shower steps
if( !fTrack)
fTrack = new ShowerTrack();
else
fTrack->Reset();
Int_t cycle(0);
while (kTRUE) {
if(cycle++>100000)
return kFALSE;
if(!GetNextStep())
break;
if( fInCome) {
ShowerStep s = GetShowerStep();
fTrack->Add(s);
}
}
GetShowerInfo();
return kTRUE;
}
//_________________________________________________________________________________________
MCTruth* SlastShowerGenerator::GetTruth() {
// Method returning the Truth object
if (!fTrack->Size()) return NULL;
if (!fTruth)
fTruth = new MCTruth();
fTruth->SetEnergy(fTrack->GetEnergy());
fTruth->SetThetaPhi(fTrack->GetTheta(),fTrack->GetPhi());
fTruth->SetFirstInt(fInitPoint,fX1);
const ShowerStep& s = fTrack->GetLastStep();
fTruth->SetEarthImpact(fEarthImpact,s.GetAgef());
// get the shower maximum
Float_t MaxElectrons(0);
Int_t MaxIndex(-1);
for (UInt_t i(0); i<fTrack->Size(); i++) {
const ShowerStep& s = fTrack->GetStep(i);
if (MaxElectrons < s.GetNelectrons()) {
MaxIndex = i;
MaxElectrons = s.GetNelectrons();
}
}
if (MaxIndex != -1) {
const ShowerStep& s = fTrack->GetStep(MaxIndex);
fTruth->SetShowerMax(0.5*(s.GetXYZi() + s.GetXYZf()),0.5*(s.GetXi() + s.GetXf()));
}
return fTruth;
}
