// ESAF : Euso Simulation and Analysis Framework
// $Id: TestLightSource.cc,v 1.31 2005/04/14 12:44:17 moreggia Exp $
// Anne Stutz created Dec, 2 2003
//
// This class brings test photons in atmosphere to Euso
// The type of events supported are the following:
// SPOT light from a single point at the same time
// RANDSPOT light from a random point at the same time
// TRACK light making a track in atmosphere
// RANDTRACK
//
// The parameters needed are the following :
// TestLightSource.Type see above can be SPOT or TRACK
// TestLightSource.Photons number of photons generated per event (if<1000 SinglePhoton else BunchOfPhotons)
// TestLightSource.Theta1 lower zenith angle
// TestLightSource.Theta2 higher zenith angle
// TestLightSource.Phi1 lower azimuth angle
// TestLightSource.Phi2 higher azimuth angle
// TestLightSource.FirstPointX X first point source in km
// TestLightSource.FirstPointY Y first point source in km
// TestLightSource.FirstPointZ Z first point source in km
// TestLightSource.LongExtension = 0 longitudinal extension of the bunch in g/cm2 can be 0
// TestLightSource.LatDistribution = NKG2 lateral distribution of the bunch can be NKG NULL
// TestLightSource.AngDistribution = baltru angular distribution of the bunch can be baltru NULL
// TestLightSource.SpectrumType can be FLUO,MONO,CERENKOV
// TestLightSource.Lambda Wavelenght if MONO
// TestLightSource.DirectionType can be EUSO,ISOTROPIC,MONO
#include "TestLightSource.hh"
#include "ListPhotonsInAtmosphere.hh"
#include "EsafRandom.hh"
#include "Config.hh"
#include "EsafSpectrum.hh"
#include "FluoCalculator.hh"
#include "LightSourceFactory.hh"
#include "TFormula.h"
#include "TMath.h"
#include "EConst.hh"
#include "EarthVector.hh"
#include "RadiativeFactory.hh"
#include "Ground.hh"
#include "Atmosphere.hh"
ClassImp(TestLightSource)
using EConst::C;
/* dNe/dr derived from 'standard' NGK formula for Lateral distribution. r is the distance to shower axis
Rm=moliere radius is a parameter . The integration of NGK2 over r from 0. to infinity is normalized to 1 . The integral
from radius r1 to radius r2 give the fraction of total number of electrons which lie inside such interval.*/
Double_t TestLightSourceNKG2(Double_t *x, Double_t *par) {
Double_t R = x[0], age = x[1], Rm=par[0];
Double_t e1=2.0;
Double_t e2=4.5;
Double_t D=R/Rm;
return TMath::Gamma(e2-age)/TMath::Gamma(age)/TMath::Gamma(e2-2.*age)*TMath::Power(D,age-(e1-1.0))
*TMath::Power((1.0+D),age-e2)/Rm;
}
/* dNe/dtheta(theta,Et) from Baltrusaitis et al. J.Phys.G:Nucl. Phys. 13 (1987) where theta is the
angle between the electrons and the shower axis and Et (MeV) is the energy thr for electrons considered.
The integration of dNe/dtheta(theta,Et) over dtheta from 0 to pi is normalized to 1.
The integral from theta1 to theta2 give the fraction of total number of electrons which lie inside
such angular interval (distribution in phi is supposed uniform)
.*/
Double_t TestLightSourceBaltru(Double_t *x, Double_t *par) {
Double_t theta = x[0], Et = x[1];
Double_t a=0.85;
Double_t b=-0.66;
Double_t theta0=a*TMath::Power(Et,-b);
return TMath::Exp(-theta/theta0)/theta0/(1.-TMath::Exp(-180/theta0) );
}
//_________________________________________________________________________________________________
TestLightSource::TestLightSource() : LightSource("TEST"), EsafMsgSource(), fFluocalcul(0) {
// ctor
Msg(EsafMsg::Info) << "Start Building TestLightSource" << MsgDispatch;
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) << "TestLightSource built"<<MsgDispatch;
}
//_________________________________________________________________________________________________
TestLightSource::~TestLightSource() {
//
// dtor
//
SafeDelete(fPh_in_atmo);
SafeDelete(fFluocalcul);
SafeDelete(fLateralDistribution);
SafeDelete(fAngularDistribution);
}
//________________________________________________________________________
void TestLightSource::Configure() {
//
// Configure TestLightSource
//
// Get fluorescence calculator
string fluoname = Conf()->GetStr("TestLightSource.FluoCalculator");
fFluocalcul = LightSourceFactory::Get()->GetFluoCalculator( fluoname );
Msg(EsafMsg::Info) << "Fluo calculator name " <<fFluocalcul->GetName()<< MsgDispatch;
// Lateral and Angular distributions
fLateralDistribution = NULL;
fAngularDistribution = NULL;
string LateralDistributionName = Conf()->GetStr("TestLightSource.LatDistribution");
string AngularDistributionName = Conf()->GetStr("TestLightSource.AngDistribution");
if (LateralDistributionName == "NKG2" )
fLateralDistribution = new TF2("LatDist",TestLightSourceNKG2,0.001,5000,0,2,1);
if (AngularDistributionName == "baltru" )
fAngularDistribution = new TF2("AngDist",TestLightSourceBaltru,0.,TMath::Pi(),.5,1000.);
}
//_________________________________________________________________________________________________
void TestLightSource::Reset() {
//
// reset internal list of photons
//
if(fPh_in_atmo) fPh_in_atmo->Reset();
if(fFluocalcul) fFluocalcul->Reset();
}
//_________________________________________________________________________________________________
PhotonsInAtmosphere *TestLightSource::Get( const PhysicsData *dummy) {
//
// return the list of photons in atmosphere
//
Reset();
string type = Conf()->GetStr("TestLightSource.Type");
Double_t nbPhotons = Conf()->GetNum("TestLightSource.Photons");
Double_t fFirstPointX = (Conf()->GetNum("TestLightSource.FirstPointX"))*km;
Double_t fFirstPointY = (Conf()->GetNum("TestLightSource.FirstPointY"))*km;
Double_t fFirstPointZ = (Conf()->GetNum("TestLightSource.FirstPointZ"))*km;
TRandom* rndm = EsafRandom::Get();
EarthVector Posi;
if (fFirstPointZ > 100*km ) {
Msg(EsafMsg::Warning) << "In TestLightSource FirstPointZ set at 100 km" << MsgDispatch;
fFirstPointZ = 100*km;
}
Posi.SetXYZ(fFirstPointX,fFirstPointY,fFirstPointZ);
// spot in fixed position
if ( type == "SPOT" ) {
MakeSpot(Posi,nbPhotons);
}
else if ( type == "RANDSPOT" ) { // spot in random position // FIXME not uniform
return NULL;
}
// track in fixed position
else if ( type == "TRACK" ) {
// random direction between theta1 and theta2, and between phi1 and phi2
Double_t theta1 = TMath::DegToRad() *Conf()->GetNum("TestLightSource.Theta1");
Double_t theta2 = TMath::DegToRad() *Conf()->GetNum("TestLightSource.Theta2");
Double_t phi1 = TMath::DegToRad() *Conf()->GetNum("TestLightSource.Phi1");
Double_t phi2 = TMath::DegToRad() *Conf()->GetNum("TestLightSource.Phi2");
Double_t thmax = theta2;
Double_t thmin = theta1;
if (theta1>theta2) {
thmax = thmin;
thmin = theta2;
}
Double_t theta = thmin + rndm->Rndm()*(thmax-thmin);
Double_t phmax = phi2;
Double_t phmin = phi1;
if (phi1>phi2) {
phmax = phmin;
phmin = phi2;
}
Double_t phi = phmin + rndm->Rndm()*(phmax-phmin);
MakeTrack(theta,phi,Posi,nbPhotons);
}
else if ( type == "RANDTRACK" ) { // track in random position
return NULL;
}
else {
Msg(EsafMsg::Panic)<<"Invalid parameter TestLightSource.Type = " << type <<MsgDispatch;
return (PhotonsInAtmosphere*)0;
}
return fPh_in_atmo;
}
//___________________________________________________________________________________________
PhotonsInAtmosphere *TestLightSource::MakeSpot(const EarthVector& Posi, Double_t nbPhotons) {
// Produce light in a single spot
Msg(EsafMsg::Info)<<"TestLightSource SPOT CALLED at (km) "<<Posi.X()/km<<" "
<<Posi.Y()/km<<" "<<Posi.Z()/km<<MsgDispatch;
string directiontype = Conf()->GetStr("TestLightSource.DirectionType");
Double_t theta1 = TMath::DegToRad() *Conf()->GetNum("TestLightSource.Theta1");
Double_t theta2 = TMath::DegToRad() *Conf()->GetNum("TestLightSource.Theta2");
Double_t phi1 = TMath::DegToRad() *Conf()->GetNum("TestLightSource.Phi1");
Double_t phi2 = TMath::DegToRad() *Conf()->GetNum("TestLightSource.Phi2");
Double_t thmax = theta2;
Double_t thmin = theta1;
Double_t phmax = phi2;
Double_t phmin = phi1;
if (theta1>theta2) {
thmax = thmin;
thmin = theta2;
}
if (phi1>phi2) {
phmax = phmin;
phmin = phi2;
}
if (directiontype == "MONO") {
Msg(EsafMsg::Info)<<" with theta between "<<thmin<<" and "<<thmax<<MsgDispatch;
Msg(EsafMsg::Info)<<" with phi between "<<phmin<<" and "<<phmax<<MsgDispatch;
}
else if (directiontype == "ISOTROPIC") {
Msg(EsafMsg::Info)<<" with isotropic direction" << MsgDispatch;
}
else if (directiontype == "EUSO") {
Msg(EsafMsg::Info)<<" direct to EUSO" << MsgDispatch;
}
else Msg(EsafMsg::Panic)<<"Invalid parameter TestLightSource.DirectionType = " << directiontype <<MsgDispatch;
// wavelenght spectrum
EsafSpectrum spectrum(357*nm);
Double_t TotalYield = MakeSpectrum(Posi,&spectrum);
//generation of SinglePhoton or BunchOfPhotons
TRandom* rndm = EsafRandom::Get();
if (nbPhotons<1000) {
// only single photons direct to Euso
if (directiontype == "ISOTROPIC" || directiontype == "MONO")
Msg(EsafMsg::Info)<<"In TestLightSource SinglePhoton are only emitted direct to Euso"<< MsgDispatch;
for (int i=0; i<(int)nbPhotons; i++) {
Double_t wl = spectrum.GetLambda();
EarthVector dir(0,0,0);
dir = EUSO()-Posi;
SinglePhoton *p = new SinglePhoton(0,wl,Posi,dir,Fluo,Direct);
fPh_in_atmo->Add(p);
}
}
else {
// one bunch of nbPhotons
EarthVector dir(1);
dir.SetTheta(thmin + rndm->Rndm()*(thmax-thmin));
dir.SetPhi(phmin + rndm->Rndm()*(phmax-phmin));
EarthVector Posf = GetLongitudinalExtension(Posi,dir);
Double_t tf = (Posf-Posi).Mag()/C();
PhotonType type=Fluo;
if (directiontype == "MONO") type = Cerenkov;
BunchOfPhotons* b = new BunchOfPhotons(nbPhotons,TotalYield,Posi,Posf,0,tf,spectrum,dir,type);
if( fLateralDistribution ) b->SetParentLateral(GetLateralDistribution(Posi));
if( fAngularDistribution ) b->SetParentAngular(GetAngularDistribution());
fPh_in_atmo->Add(b);
}
return fPh_in_atmo;
}
//_____________________________________________________________________________________________________
PhotonsInAtmosphere *TestLightSource::MakeTrack(const Double_t& theta, const Double_t& phi, const EarthVector& Posi, Double_t nbPhotons) {
//
// Produce light along a track
//
Msg(EsafMsg::Info)<<"TestLightSource TRACK CALLED theta="<<theta<<" phi ="<<phi<<MsgDispatch;
Msg(EsafMsg::Info)<<" Starting point at (km) "<<Posi.X()/km<<" "<<Posi.Y()/km<<" "<<Posi.Z()/km<<MsgDispatch;
string directiontype = Conf()->GetStr("TestLightSource.DirectionType");
EarthVector DirTrack(1);
DirTrack.SetTheta(theta);
DirTrack.SetPhi(phi);
EsafSpectrum spectrum(357*nm);
Double_t TotalYield;
//generation of SinglePhoton or BunchOfPhotons
// impact on ground
Ground* ground = RadiativeFactory::Get()->GetGround();
EarthVector impact = ground->GetImpact(Posi,DirTrack);
Msg(EsafMsg::Info)<<" Impact on ground at "<<impact.X()/km<<" "<<impact.Y()/km<<" "<<impact.Z()/km<< MsgDispatch;
Double_t magmax = (Posi - (ground->GetImpact(Posi,DirTrack))).Mag();
// only single photons
if (nbPhotons<1000) {
Msg(EsafMsg::Info)<<" direct to EUSO" << MsgDispatch;
if (directiontype == "ISOTROPIC" || directiontype == "MONO")
Msg(EsafMsg::Info)<<"In TestLightSource SinglePhoton are only emitted direct to Euso"<< MsgDispatch;
for (int i=0; i<(int)nbPhotons; i++) {
Double_t mag = magmax*(i+1)/nbPhotons;
EarthVector pos = Posi + mag*DirTrack;
Double_t t = mag/C();
TotalYield = MakeSpectrum(pos,&spectrum);
Double_t wl = spectrum.GetLambda();
EarthVector dir = EUSO()-pos;
SinglePhoton *p = new SinglePhoton(t,wl,pos,dir,Fluo,Direct);
fPh_in_atmo->Add(p);
}
}
// 100 bunchs of photons
else {
if (directiontype == "EUSO")
Msg(EsafMsg::Info)<<"In TestLightSource BunchPhoton are not emitted directly to Euso but isotropically"<< MsgDispatch;
Double_t nbPhotonsInBunch = nbPhotons/100.;
PhotonType type=Fluo;
if (directiontype == "MONO") type = Cerenkov;
for (Float_t i=0; i<100; i++) {
Double_t mag = magmax*(i+1)/100;
EarthVector pos = Posi + mag*DirTrack;
Double_t t = mag/C();
TotalYield = MakeSpectrum(pos,&spectrum);
EarthVector posf = GetLongitudinalExtension(pos,DirTrack);
if ( posf == pos ) break;
Double_t tf = t + (posf-pos).Mag()/C();
BunchOfPhotons* b = new BunchOfPhotons(nbPhotonsInBunch,TotalYield,pos,posf,t,tf,spectrum,DirTrack,type);
if( fLateralDistribution ) b->SetParentLateral(GetLateralDistribution(Posi));
if( fAngularDistribution ) b->SetParentAngular(GetAngularDistribution());
fPh_in_atmo->Add(b);
}
}
return fPh_in_atmo;
}
//_________________________________________________________________________________________________
Double_t TestLightSource::MakeSpectrum(const EarthVector& pos,EsafSpectrum* spectrum) {
// Calculation of the wavelenght spectrum
Double_t TotalYield=0;
string spectrumtype = Conf()->GetStr("TestLightSource.SpectrumType");
if (spectrumtype == "MONO") {
double lambda = (Conf()->GetNum("TestLightSource.Lambda"))*nm;
if ( spectrum!=0) spectrum -> Reset(lambda);
TotalYield = 1.;
}
else if (spectrumtype == "FLUO") {
Double_t energy=80.*MeV;
TotalYield = fFluocalcul->GetFluoYield(pos.Zv(),energy,spectrum);
}
else if (spectrumtype == "CERENKOV") {
TFormula cerenkov("cerenkov","1 /(x*x)");
if ( spectrum !=0 ) spectrum -> Reset(&cerenkov,150,300*nm,450*nm);
TotalYield = 1.;
}
else Msg(EsafMsg::Panic)<<"Invalid parameter TestLightSource.SpectrumType = " << spectrumtype <<MsgDispatch;
if(!spectrum) Msg(EsafMsg::Panic)<<"Pb of memory allocation in TestLightSource"<<MsgDispatch;
return TotalYield;
}
//_________________________________________________________________________________________________
const TF12 TestLightSource::GetLateralDistribution(const EarthVector& pos) {
//
// return TF12 the lateral distribution at pos for age = 1
//
TF12 LatDist;
Double_t age = 1.;
// if NKG in meter gives moliere radius as parameter in meter
if ( fLateralDistribution) {
if ( fLateralDistribution->GetNpar() == 1 ) {
// Get Atmosphere for Density
const Atmosphere* atmo = Atmosphere::Get();
Double_t Rm = 9.6 * gram/cm2 / atmo->Air_Density( pos.Zv() ) / m; // in meters
Msg(EsafMsg::Debug)<<"rayon moliere = " << Rm <<MsgDispatch;
fLateralDistribution->SetParameters(Rm,0);
}
else Msg(EsafMsg::Panic)<<"Lateral Distribution should have 1 parameter" <<MsgDispatch;
LatDist = TF12("LDS_name",fLateralDistribution,age,"X");
}
else Msg(EsafMsg::Panic)<<"No Lateral Distribution" <<MsgDispatch;
return LatDist;
}
//_________________________________________________________________________________________________
const TF12 TestLightSource::GetAngularDistribution() {
//
// return Pointer to the angular distribution
//
TF12 AngDist;
Double_t EthCer = 20.;
if ( fAngularDistribution ) AngDist = TF12("ADS_name",fAngularDistribution,EthCer,"X");
else Msg(EsafMsg::Panic)<<"No Angular Distribution "<<MsgDispatch;
return AngDist;
}
//__________________________________________________________________________________________________
EarthVector TestLightSource::GetLongitudinalExtension(const EarthVector& Pos,const EarthVector& Dir) {
//
// return last point of a bunch with first point Pos and mean direction Dir
//
Double_t LgExt = (Conf()->GetNum("TestLightSource.LongExtension"))*g/cm2;
Double_t L,h(0),depth(0);
EarthVector Posf = Pos;
L = LgExt / Atmosphere::Get()->Air_Density(h)/100.;
Int_t cycle=0;
while(1) {
Posf += Dir*L;
if ( Posf.IsUnderSeaLevel() ) return Pos;
depth += Atmosphere::Get()->Air_Density(Posf.Zv())*L;
if ( depth >= LgExt ) break;
if ( cycle>100000 ) {
MsgForm( EsafMsg::Debug,"Next point not found : cycle > 100000 with Lgext %f",LgExt/g*cm2 );
return Pos;
}
}
return Posf;
}
//__________________________________________________________________________________________________
EarthVector TestLightSource::EUSO() const {
//
// TOOL FUNCTION : Returns the EUSO position
//
static Int_t counter = 0;
static EarthVector pos(0,0,0);
if(!counter) {
ConfigFileParser* pConf = Config::Get()->GetCF("General","Euso");
pos.SetZ(pConf->GetNum("Euso.fAltitude")*km);
}
if(!counter) counter++;
return pos;
}