// ESAF : Euso Simulation and Analysis Framework
// $Id: TestLightSource.cc,v 1.40 2005/10/27 13:54:35 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.Option 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 namespace TMath;
using namespace sou;
using namespace EConst;
/* 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 Gamma(e2-age)/Gamma(age)/Gamma(e2-2.*age)*Power(D,age-(e1-1.0))
*Power((1.0+D),age-e2)/Rm;
}
Double_t TestLightSourceBaltru(Double_t *x, Double_t *par) {
//
// 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 theta = x[0], Et = x[1];
Double_t a = 0.85;
Double_t b = 0.66;
Double_t theta0 = a*Power(Et,-b);
Double_t pigreco = Pi();
return Exp(-theta/theta0) / theta0 / (1.- Exp(-pigreco/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 detector position
ConfigFileParser* pConf = Config::Get()->GetCF("General","Euso");
fEUSO.SetZ(pConf->GetNum("Euso.fAltitude")*km);
// 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.,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 option = Conf()->GetStr("TestLightSource.Option");
Double_t nbPhotons = Conf()->GetNum("TestLightSource.Nbph");
Double_t fFirstPointX = (Conf()->GetNum("TestLightSource.FirstPointX"))*km;
Double_t fFirstPointY = (Conf()->GetNum("TestLightSource.FirstPointY"))*km;
Double_t fFirstPointZ = (Conf()->GetNum("TestLightSource.FirstPointZ"))*km;
string impactMode = Conf()->GetStr("TestLightSource.ImpactMode");
Double_t ImpactXmin = Conf()->GetNum("TestLightSource.ImpactXmin")*km;
Double_t ImpactXmax = Conf()->GetNum("TestLightSource.ImpactXmax")*km;
Double_t ImpactYmin = Conf()->GetNum("TestLightSource.ImpactYmin")*km;
Double_t ImpactYmax = Conf()->GetNum("TestLightSource.ImpactYmax")*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 ( option == "SPOT" ) MakeSpot(Posi,nbPhotons);
// spot mimicking shower maximum, depends on theta
else if ( option == "HmaxSPOT" ) {
Double_t x = rndm->Rndm()*(ImpactXmax - ImpactXmin) + ImpactXmin;
Double_t y = rndm->Rndm()*(ImpactYmax - ImpactYmin) + ImpactYmin;
Double_t z = - (x*x + y*y)/(2*EarthRadius()); // approched value
MakeHmaxSpot(EarthVector(x,y,z),nbPhotons);
}
// track in fixed position
else if ( option == "TRACK" ) {
// random direction between theta1 and theta2, and between phi1 and phi2
Double_t theta1 = DegToRad() *Conf()->GetNum("TestLightSource.Theta1");
Double_t theta2 = DegToRad() *Conf()->GetNum("TestLightSource.Theta2");
Double_t phi1 = DegToRad() *Conf()->GetNum("TestLightSource.Phi1");
Double_t phi2 = 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 {
Msg(EsafMsg::Panic)<<"Invalid parameter TestLightSource.Option = " << option <<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");
string photontype = Conf()->GetStr("TestLightSource.PhotonType");
string photonDescription = Conf()->GetStr("TestLightSource.Descrip");
Double_t theta1 = DegToRad() *Conf()->GetNum("TestLightSource.Theta1");
Double_t theta2 = DegToRad() *Conf()->GetNum("TestLightSource.Theta2");
Double_t phi1 = DegToRad() *Conf()->GetNum("TestLightSource.Phi1");
Double_t phi2 = DegToRad() *Conf()->GetNum("TestLightSource.Phi2");
PhotonType phtype = Fluo;
if(photontype == "cerenkov") phtype = Cerenkov;
TRandom* rndm = EsafRandom::Get();
string thetaMode = Conf()->GetStr("TestLightSource.ThetaMode");
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;
}
Double_t theta = thmin + rndm->Rndm()*(thmax-thmin);
Double_t phi = phmin + rndm->Rndm()*(phmax-phmin);
EarthVector Omega(1);
if(thetaMode == "local") {
EarthVector v1 = Posi.Unit();
EarthVector vrot;
Omega.SetXYZ(1,0,0);
EarthVector Uz(0,0,1);
vrot = Uz.Cross(v1);
if(vrot.Mag() > 0) {
Omega.Rotate(v1.Theta(),vrot);
Omega.Rotate(phi,v1);
vrot = v1.Cross(Omega);
Omega.Rotate(Pi()/2+theta,vrot);
}
else Omega.SetMagThetaPhi(1.,Pi() - theta,phi + Pi());
theta = Omega.Theta();
phi = Omega.Phi();
}
else if(thetaMode == "mes") {
theta = Pi() - theta;
phi += Pi();
}
else Msg(EsafMsg::Panic)<<"Invalid parameter TestLightSource.ThetaMode = " << thetaMode <<MsgDispatch;
if (directiontype == "UNIQUE") {
Msg(EsafMsg::Info)<<" with theta between "<<thmin/deg<<" and "<<thmax/deg<<MsgDispatch;
Msg(EsafMsg::Info)<<" with phi between "<<phmin/deg<<" and "<<phmax/deg<<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,photontype);
//generation of SinglePhoton or BunchOfPhotons
if(photonDescription == "single") {
// only single photons direct to Euso
if (directiontype == "ISOTROPIC" || directiontype == "UNIQUE")
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,phtype,Direct);
fPh_in_atmo->Add(p);
}
}
// one bunch of nbPhotons
else if(photonDescription == "bunch"){
EarthVector dir(1);
dir.SetTheta(theta);
dir.SetPhi(phi);
EarthVector Posf = GetLongitudinalExtension(Posi,dir);
Double_t tf = (Posf-Posi).Mag()/Clight();
BunchOfPhotons* b = new BunchOfPhotons(nbPhotons,TotalYield,Posi,Posf,0,tf,spectrum,dir,phtype);
if( fLateralDistribution ) b->SetParentLateral(GetLateralDistribution(Posi));
if( fAngularDistribution && (phtype==Cerenkov) ) b->SetParentAngular(GetAngularDistribution(Posi.Zv()));
fPh_in_atmo->Add(b);
// track needed for ckov AND fluo (for last transfer until detector)
BuildLightTrack(Posi,dir);
//DELETE
// for CERENKOV RadiativeTransfer handling, when AlongTrack_CSPropagator is used
//if(phtype == Cerenkov) BuildLightTrack(Posi,dir); //DELETE
}
else Msg(EsafMsg::Panic)<<"<MakeSpot> Invalid photonDescription parameter = " << photonDescription<<MsgDispatch;
return fPh_in_atmo;
}
//___________________________________________________________________________________________
PhotonsInAtmosphere *TestLightSource::MakeHmaxSpot(const EarthVector& impact, Double_t nbPhotons) {
//
// Produce light in a single spot
//
Msg(EsafMsg::Info)<<"TestLightSource HmaxSPOT with impact at (km) "<<impact.X()/km<<" "
<<impact.Y()/km<<" "<<impact.Z()/km<<MsgDispatch;
string directiontype = Conf()->GetStr("TestLightSource.DirectionType");
string photontype = Conf()->GetStr("TestLightSource.PhotonType");
string photonDescription = Conf()->GetStr("TestLightSource.Descrip");
Double_t theta1 = DegToRad() *Conf()->GetNum("TestLightSource.Theta1");
Double_t theta2 = DegToRad() *Conf()->GetNum("TestLightSource.Theta2");
Double_t phi1 = DegToRad() *Conf()->GetNum("TestLightSource.Phi1");
Double_t phi2 = DegToRad() *Conf()->GetNum("TestLightSource.Phi2");
PhotonType phtype = Fluo;
if(photontype == "cerenkov") phtype = Cerenkov;
TRandom* rndm = EsafRandom::Get();
// bunch direction (for cerenkov bunches only)
string thetaMode = Conf()->GetStr("TestLightSource.ThetaMode");
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(thmax > 80*DegToRad()) {
thmax = 80*DegToRad();
Msg(EsafMsg::Warning)<<"With HmaxSPOT, theta limited to 80 deg -> theta=80 applied"<<MsgDispatch;
}
Double_t theta_true = thmin + rndm->Rndm()*(thmax-thmin);
Double_t phi = phmin + rndm->Rndm()*(phmax-phmin);
Double_t theta = theta_true;
EarthVector Omega(1);
// if theta given is local -> theta in MES is calculated
if(thetaMode == "local") {
EarthVector v1 = impact.Unit();
EarthVector vrot;
Omega.SetXYZ(1,0,0);
EarthVector Uz(0,0,1);
if(vrot.Mag() > 0) {
Omega.Rotate(v1.Theta(),vrot);
Omega.Rotate(phi,v1);
vrot = v1.Cross(Omega);
Omega.Rotate(Pi()/2+theta,vrot);
}
else Omega.SetMagThetaPhi(1.,Pi() - theta,phi + Pi());
theta = Omega.Theta();
phi = Omega.Phi();
}
else if(thetaMode == "mes") Msg(EsafMsg::Panic)<<"mes ThetaMode NOT possible with HmaxSPOT option"<<MsgDispatch;
else Msg(EsafMsg::Panic)<<"Invalid parameter TestLightSource.ThetaMode = " << thetaMode <<MsgDispatch;
// get Hmax from theta_true
EarthVector Posi = GetHmaxPos(impact,Omega,theta_true);
Msg(EsafMsg::Info)<<"TestLightSource HmaxSPOT with INITIAL pos at (km) "<<Posi.X()/km<<" " <<Posi.Y()/km<<" "<<Posi.Z()/km<<MsgDispatch;
if (directiontype == "UNIQUE") {
Msg(EsafMsg::Info)<<" with theta between "<<thmin/deg<<" and "<<thmax/deg<<MsgDispatch;
Msg(EsafMsg::Info)<<" with phi between "<<phmin/deg<<" and "<<phmax/deg<<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,photontype);
// only single photons direct to Euso
if(photonDescription == "single") {
if (directiontype == "ISOTROPIC" || directiontype == "UNIQUE")
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,phtype,Direct);
fPh_in_atmo->Add(p);
}
}
// one bunch of nbPhotons
else if(photonDescription == "bunch"){
EarthVector dir(1);
dir.SetTheta(Omega.Theta());
dir.SetPhi(Omega.Phi());
EarthVector Posf = GetLongitudinalExtension(Posi,dir);
Double_t tf = (Posf - Posi).Mag()/Clight();
BunchOfPhotons* b = new BunchOfPhotons(nbPhotons,TotalYield,Posi,Posf,0,tf,spectrum,dir,phtype);
if( fLateralDistribution ) b->SetParentLateral(GetLateralDistribution(Posi));
if( fAngularDistribution && (phtype==Cerenkov) ) b->SetParentAngular(GetAngularDistribution(Posi.Zv()));
fPh_in_atmo->Add(b);
// track needed for ckov AND fluo (for last transfer until detector)
BuildLightTrack(Posi,dir);
//DELETE
// for CERENKOV RadiativeTransfer handling, when AlongTrack_CSPropagator is used
//if(phtype == Cerenkov) BuildLightTrack(Posi,dir); //DELETE
}
else Msg(EsafMsg::Panic)<<"<MakeHmaxSpot> Invalid photonDescription parameter = " << photonDescription<<MsgDispatch;
return fPh_in_atmo;
}
//_____________________________________________________________________________________________________
PhotonsInAtmosphere *TestLightSource::MakeTrack(Double_t theta, 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");
string photontype = Conf()->GetStr("TestLightSource.PhotonType");
string photonDescription = Conf()->GetStr("TestLightSource.Descrip");
PhotonType phtype = Fluo;
if(photontype == "cerenkov") phtype = Cerenkov;
string thetaMode = Conf()->GetStr("TestLightSource.ThetaMode");
EarthVector Omega(1);
if(thetaMode == "local") {
EarthVector v1 = Posi.Unit();
EarthVector vrot;
Omega.SetXYZ(1,0,0);
EarthVector Uz(0,0,1);
if(vrot.Mag() > 0) {
Omega.Rotate(v1.Theta(),vrot);
Omega.Rotate(phi,v1);
vrot = v1.Cross(Omega);
Omega.Rotate(Pi()/2+theta,vrot);
}
else Omega.SetMagThetaPhi(1.,Pi() - theta,phi + Pi());
theta = Omega.Theta();
phi = Omega.Phi();
}
else if(thetaMode == "mes") {
theta = Pi() - theta;
phi += Pi();
}
else Msg(EsafMsg::Panic)<<"Invalid parameter TestLightSource.ThetaMode = " << thetaMode <<MsgDispatch;
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(photonDescription == "single") {
Msg(EsafMsg::Info)<<" direct to EUSO" << MsgDispatch;
if (directiontype == "ISOTROPIC" || directiontype == "UNIQUE")
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/Clight();
TotalYield = MakeSpectrum(pos,&spectrum,photontype);
Double_t wl = spectrum.GetLambda();
EarthVector dir = EUSO()-pos;
SinglePhoton *p = new SinglePhoton(t,wl,pos,dir,phtype,Direct);
fPh_in_atmo->Add(p);
}
}
// 100 bunchs of photons
else if(photonDescription == "bunch"){
if (directiontype == "EUSO") Msg(EsafMsg::Info)<<"In TestLightSource BunchPhoton are not emitted directly to Euso"<< MsgDispatch;
Double_t nbPhotonsInBunch = nbPhotons/100.;
for (Float_t i=0; i<100; i++) {
Double_t mag = magmax*(i+1)/100;
EarthVector pos = Posi + mag*DirTrack;
Double_t t = mag/Clight();
TotalYield = MakeSpectrum(pos,&spectrum,photontype);
EarthVector posf = GetLongitudinalExtension(pos,DirTrack);
if ( posf == pos ) break;
Double_t tf = t + (posf-pos).Mag()/Clight();
BunchOfPhotons* b = new BunchOfPhotons(nbPhotonsInBunch,TotalYield,pos,posf,t,tf,spectrum,DirTrack,phtype);
if( fLateralDistribution ) b->SetParentLateral(GetLateralDistribution(pos));
if( fAngularDistribution && (phtype==Cerenkov) ) b->SetParentAngular(GetAngularDistribution(Posi.Zv()));
fPh_in_atmo->Add(b);
}
// for CERENKOV RadiativeTransfer handling, when AlongTrack_CSPropagator is used
if(phtype == Cerenkov) BuildLightTrack(Posi,DirTrack);
}
else Msg(EsafMsg::Panic)<<"<MakeTrack> Invalid photonDescription parameter = " << photonDescription<<MsgDispatch;
return fPh_in_atmo;
}
//_________________________________________________________________________________________________
Double_t TestLightSource::MakeSpectrum(const EarthVector& pos,EsafSpectrum* spectrum,string photontype) {
//
// 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 == "COMPLETE" && photontype == "fluo") {
Double_t energy=80.*MeV;
TotalYield = fFluocalcul->GetFluoYield(pos.Zv(),energy,spectrum);
}
else if (spectrumtype == "COMPLETE" && photontype == "cerenkov") {
TFormula cerenkov("cerenkov","1 /(x*x)");
if ( spectrum !=0 ) spectrum -> Reset(&cerenkov,100,300*nm,450*nm);
TotalYield = 1.;
}
else Msg(EsafMsg::Panic)<<"Invalid parameter TestLightSource.SpectrumType or PhotonType = " << spectrumtype<<", "<<photontype <<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(Double_t alt) {
//
// return Pointer to the angular distribution
//
TF12 AngDist;
Double_t EthCer = GetEnergyThreshold(alt);
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;
EarthVector Posf = Pos;
/*
Double_t L,h(0),depth(0);
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;
}
}
*/
Int_t status = Atmosphere::Get()->InvertGrammage(Pos,Dir,LgExt,Posf);
if(status != 0) {
Msg(EsafMsg::Warning) << "<GetLongitudinalExtension> : CANNOT FIND FINAL POSITION, InvertGrammage status = "<<status<<MsgDispatch;
return Pos;
}
return Posf;
}
//__________________________________________________________________________________________________
EarthVector TestLightSource::GetHmaxPos(const EarthVector& impact,const EarthVector& Direc,Double_t theta) const {
//
// get Hmax 3D-position -- for HmaxSPOT option only
// Here theta MUST be LOCAL zenith angle (at impact)
//
EarthVector rtn(0,0,0);
EarthVector Dir(Direc);
if(Dir.Mag() != 1) Dir = Dir.Unit();
Double_t tolerance = 0.5*m;
if(theta > 80) Msg(EsafMsg::Warning) << "theta > 80deg NOT foreseen in <TestLightSource::GetHmaxPos>"<<MsgDispatch;
// get Hmax value from LOCAL theta
Double_t Hmax = (1.89 - 7.59*Log(Cos(theta))) * km;
// get 3D-position
EarthVector first = impact;
EarthVector last = impact - (30*km/Cos(Dir.Theta() - Pi()))*Dir;
EarthVector middle(1.);
while((last.Zv() - first.Zv()) > tolerance) {
middle = 0.5*(last+first);
if(fabs(middle.Zv() - Hmax) < tolerance) {rtn = middle; break;}
if(Hmax < middle.Zv()) last = middle;
else first = middle;
}
if(rtn.Mag() == 0) rtn = first;
#ifdef DEBUG
cout<<" IN GETHMAXPOS :"<<endl;
cout<<"impact position = "<<impact<<endl;
cout<<"Hmax foreseen = "<<Hmax<<endl;
cout<<"Hmax implemented = "<<rtn.Zv()<<endl;
#endif
return rtn;
}
//__________________________________________________________________________________________________
void TestLightSource::BuildLightTrack(const EarthVector& posinit, const EarthVector& dir) {
//
// for CERENKOV RadiativeTransfer handling, when AlongTrack_CSPropagator is used
//
// init
ConfigFileParser* pConfig = Config::Get()->GetCF("RadiativeTransfer","BunchRadiativeTransfer");
Double_t depthstep = pConfig->GetNum("BunchRadiativeTransfer.DepthStep")*g/cm2;
fPh_in_atmo->ClearTrack();
const Atmosphere* atmo = Atmosphere::Get();
EarthVector presentpos(posinit);
EarthVector nextpos(posinit);
EarthVector posi(posinit);
if(posi.IsUnderSeaLevel()) {
Msg(EsafMsg::Warning) << "<BuildLightTrack> Starting pos is under sea level -> set to Nadir"<<MsgDispatch;
posi.SetMag(0.);
}
Int_t status = -100; // for Atmosphere::InvertGrammage() status
// look at impact and total depth
EarthVector impact = atmo->ImpactASL(posi,dir);
if(impact.Z() == HUGE) impact = atmo->ImpactAtTOA(posi,dir);
if(impact.Z() == HUGE) Msg(EsafMsg::Warning) << "<BuildLightTrack> : No GROUND impact nor TOA impact -> SHOULD NOT"<<MsgDispatch;
Double_t FinalDepth = atmo->Grammage(posi,impact);
// get nb of TrackLight steps
UInt_t nb = UInt_t(floor(FinalDepth/depthstep)) + 2; // +1 for posi, +1 for tuning last step exit
if(nb > 1000) {
Msg(EsafMsg::Warning) << "<BuildLightTrack> : Tracklength is huge = " << (impact - posi).Mag()/km<<" km"<<MsgDispatch;
Msg(EsafMsg::Warning) << "<BuildLightTrack> : Corresponding depth is huge too = " << FinalDepth*cm2/g <<" g/cm2"<<MsgDispatch;
Msg(EsafMsg::Warning) << "<BuildLightTrack> : Thus nb of TRACKLIGHT steps set to 1000"<<MsgDispatch;
nb = 1000;
depthstep = FinalDepth / (nb - 2); // -1 for posi, -1 for tuning last step exit
}
fPh_in_atmo->SetNbTrackSteps(nb);
// fill the track until last step entry (last step exit done by hand below
for(UInt_t i=0; i<nb-1; i++) {
presentpos = nextpos;
fPh_in_atmo->FillTrack(presentpos);
if(i == (nb-2)) break; // to avoid line below for last step of the loop
status = atmo->InvertGrammage(presentpos,dir,depthstep,nextpos);
if(status != 0) Msg(EsafMsg::Warning) << "<BuildLightTrack> InvertGrammage() status is problematic -> status = " << status <<MsgDispatch;
}
// fine tune last step
fPh_in_atmo->FillTrack(impact);
#ifdef DEBUG
cout<<"impact position = "<<impact<<endl;
FinalDepth = atmo->Grammage(presentpos,impact);
cout<<"finaldepth w.r.t. depthstep = "<<FinalDepth/depthstep<<endl;
if(FinalDepth/depthstep > 1) Msg(EsafMsg::Debug) << "<BuildLightTrack> : Last TRACKLIGHT step tuning failed, it counts as "<<FinalDepth/depthstep<<" times a normal step"<<MsgDispatch;
#endif
}
//____________________________________________________________________________________________
Double_t TestLightSource::GetEnergyThreshold(Double_t SC_alt) const {
//
// get the energy threshold for cerenkov emission
//
// Get Atmosphere Index
const Atmosphere* atmo = Atmosphere::Get();
Double_t delta = atmo->Index_Minus1(SC_alt);
// When Index is not implemented in the Atmosphere use the Corsika formula
if ( delta == 0) delta = 0.000283 * atmo->Air_Density( SC_alt )/atmo->Air_Density(0);
return ElectronMassC2()/sqrt(2*delta);
}