#include "SlastShowerGenerator.hh"
#include "EsafRandom.hh"
#include "MCTruth.hh"
#include "ShowerTrack.hh"
#include "Atmosphere.hh"
#include "Config.hh"
#include "EConst.hh"
#include <TMath.h>
#include <TRandom.h>
#include <TProfile.h>
#include <TH1F.h>
#include <TGraph.h>
#include "NumbersFileParser.hh"
ClassImp(SlastShowerGenerator)
using namespace //root.cern.ch/root/html/TMath.html">TMath;
using namespace sou;
using namespace EConst;
SlastShowerGenerator::SlastShowerGenerator(Bool_t quiet) : EventGenerator("slast++"), EsafMsgSource(),
fTrack(0), fTruth(0), fInCome(0), fOutCome(0), fQuiet(kFALSE) {
fQuiet = quiet;
fSpectrumRdnm = 0;
fEnergy = -1;
if(!Init()) Msg(EsafMsg::Panic)<<"SlastShowerGenerator can not be initialized"<<MsgDispatch;
}
Bool_t SlastShowerGenerator::Init() {
ConfigFileParser *pConfig = Config::Get()->GetCF("LightToEuso","GeneratorLightToEuso");
fFoV = pConfig->GetNum("GeneratorLightToEuso.FoV")*deg;
fSpectrumType = pConfig->GetStr("GeneratorLightToEuso.SpectrumType");
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");
fEusoAltitude = pConfig->GetNum("GeneratorLightToEuso.altitude")*km;
fEusoVector.SetXYZ(0,0,fEusoAltitude);
fRejectFakeEvents = pConfig->GetBool("GeneratorLightToEuso.RejectFakeEvents");
fRejectNoXmax = pConfig->GetBool("GeneratorLightToEuso.RejectNoXmax");
if(!fRejectFakeEvents && fRejectNoXmax)
Msg(EsafMsg::Panic) << "<Init> if RejectFakeEvents not required, RejectNoXmax cannot be required" << MsgDispatch;
if((fThetaMax - fThetaMin) < kTolerance && fFirstType != "POS") {
Msg(EsafMsg::Warning) << "<Init> theta bounds width cannot be zero : theta max increased by 1 degree" << MsgDispatch;
fThetaMax += 1*DegToRad();
}
if((fPhiMax - fPhiMin) < kTolerance && fFirstType != "POS") {
Msg(EsafMsg::Warning) << "<Init> phi bounds width cannot be zero : phi max increased by 1 degree" << MsgDispatch;
fPhiMax += 1*DegToRad();
}
if (!fQuiet) {
Msg(EsafMsg::Info) << "SlastShowerGenerator Enabled" << MsgDispatch;
Msg(EsafMsg::Info) << "SlastShowerGenerator Enables Atmosphere " << Atmosphere::Get()->GetType() << MsgDispatch;
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);
}
else {
MsgForm(EsafMsg::Debug,"Generate %s event(s) in:",fType.c_str());
MsgForm(EsafMsg::Debug,"Energy interval [%.4E,%.4E]",fEnergyMin,fEnergyMax);
MsgForm(EsafMsg::Debug,"Theta interval [%.4f,%.4f]",fThetaMin,fThetaMax);
MsgForm(EsafMsg::Debug,"Phi interval [%.4f,%.4f]",fPhiMin,fPhiMax);
}
fAtomicMass = AtomicMass(fType);
if(fSpectrumType == "GZKHiRes2005") {
if(fEnergyMin < 1e19 || fEnergyMax > 1e21)
Msg(EsafMsg::Panic) << "<Init> GZKHiRes2005 spectrum stands within [1e19 - 1e21] only !" << MsgDispatch;
NumbersFileParser nf("config/Tools/GZKHiRes2005.dat",2);
vector<Double_t> log10E_v = nf.GetCol(0);
vector<Double_t> fluxE3_v = nf.GetCol(1);
//root.cern.ch/root/html/TArrayD.html">TArrayD log10E(log10E_v.size()), fluxE3(log10E_v.size());
for(size_t i=0; i < log10E_v.size(); i++) {
log10E[i] = log10E_v[i];
fluxE3[i] = fluxE3_v[i];
}
//root.cern.ch/root/html/TGraph.html">TGraph* gr = new TGraph(log10E_v.size(),log10E.GetArray(),fluxE3.GetArray());
fSpectrumRdnm = new //root.cern.ch/root/html/TH1F.html">TH1F("GZK","GZK",10000,Log10(fEnergyMin),Log10(fEnergyMax));
for(Int_t i=0; i < 10000; i++) {
fSpectrumRdnm->//root.cern.ch/root/html/TH1F.html#TH1F:SetBinContent">SetBinContent(i+1,gr->Eval(fSpectrumRdnm->GetBinCenter(i+1),0,"") / pow(10,fSpectrumRdnm->GetBinCenter(i+1)*3. - 24.));
}
}
pConfig = Config::Get()->GetCF("LightToEuso","SlastLightToEuso");
fEnergyDistribution = pConfig->GetStr("SlastLightToEuso.EnergyDistributionParametrization");
fShowerParametrization = pConfig->GetStr("SlastLightToEuso.ShowerParametrization");
return kTRUE;
}
SlastShowerGenerator::~SlastShowerGenerator() {
SafeDelete(fTrack);
SafeDelete(fTruth);
SafeDelete(fSpectrumRdnm);
}
void SlastShowerGenerator::SetShower(Double_t e,Double_t t,Double_t p,EarthVector pos) {
SetEnergy(fEnergyMin = fEnergyMax = e);
SetTheta(fThetaMin = fThetaMax = t);
SetPhi(fPhiMin = fPhiMax = p);
fFirstPointX = pos.X();
fFirstPointY = pos.Y();
fFirstPointZ = pos.Z();
fFirstType = "POS";
}
void SlastShowerGenerator::GetX1() {
if ( fFirstType == "POS" ) {
fX1 = Atmosphere::Get()->Grammage(fInitPoint,-fOmega,"dir");
}
else if (fFirstType == "X1" ) {
ConfigFileParser *pConfig = Config::Get()->GetCF("LightToEuso","GeneratorLightToEuso");
fX1 = pConfig->GetNum("GeneratorLightToEuso.InteractionX1")*g/cm2;
}
else {
//root.cern.ch/root/html/TRandom.html">TRandom *rndm = EsafRandom::Get();
fX1 = -Log(rndm->Rndm())*HadronInteractionLength(fAtomicMass,fEnergy)*g/cm2;
}
}
Double_t SlastShowerGenerator::HadronInteractionLength(Double_t A, Double_t E) {
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);
B = Power(AO,1./3);
S += CO*Power(B+X-q*(1./B+Y),2);
B = Power(AA,1./3);
S += CA*Power(B+X-q*(1./B+Y),2);
B = Power(AC,1./3);
S += CC*Power(B+X-q*(1./B+Y),2);
return F_A2/S*NucleonAirCrossSection(1.e11/A)/NucleonAirCrossSection(E/A);
}
Double_t SlastShowerGenerator::NucleonAirCrossSection(Double_t E) {
... 2. V.A. Naumov, T.S. Sinegovskaya (eq (17))
*/
return 290 - 8.7*Log(E*1.e-9) + 1.14*Power(Log(E*1.e-9),2);
}
const ShowerStep& SlastShowerGenerator::GetShowerStep() {
fStep.Clear();
fStep.fXi = fXcurrent + fX1;
fStep.fXf = fXnext + fX1;
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) {
if (fShowerParametrization == "GIL") GIL(x);
else Msg(EsafMsg::Panic) << "<GetShowerParametrization> Only GIL formula available so far" <<MsgDispatch;
}
Bool_t SlastShowerGenerator::FindImpactOnTOA() {
//root.cern.ch/root/html/TRandom.html">TRandom *rndm = EsafRandom::Get();
ConfigFileParser* pConfig = Config::Get()->GetCF("LightToEuso","GeneratorLightToEuso");
string ImpactMode = pConfig->GetStr("GeneratorLightToEuso.ImpactMode");
string thetaMode = pConfig->GetStr("GeneratorLightToEuso.ThetaMode");
const Atmosphere* atmo = Atmosphere::Get();
#ifdef DEBUG
Int_t counterr = 0;
#endif
if(ImpactMode == "TOA") {
Double_t TOA_alt = atmo->GetTOAAltitude();
Double_t Euso_alt = fEusoVector.Z();
Double_t FoV_radius = Euso_alt * Tan(fFoV + 1*deg);
Double_t L1 = Sqrt( (EarthRadius() + TOA_alt)*(EarthRadius() + TOA_alt) - EarthRadius()*EarthRadius() );
Double_t Zmin,x,y,z;
Zmin = EarthRadius()*Sqrt(1 + 2*TOA_alt/EarthRadius() + (TOA_alt*TOA_alt - FoV_radius*FoV_radius)/(EarthRadius()*EarthRadius()));
Zmin -= 2*L1 * Sqrt(1 - EarthRadius()*EarthRadius() / (EarthRadius() + TOA_alt) / (EarthRadius() + TOA_alt));
#ifdef DEBUG
#endif
Bool_t eventisfound = kFALSE;
Int_t infloop = 0;
while(!eventisfound) {
#ifdef DEBUG
counterr++;
#endif
if(infloop++ > 100000) {
Msg(EsafMsg::Warning) << "<FindImpactOnTOA> TOA impact mode : Inifinite loop broken by hand" << MsgDispatch;
return kFALSE;
}
z = (EarthRadius()+TOA_alt - Zmin)*rndm->Rndm() + Zmin;
Double_t phi = rndm->Rndm()*TwoPi();
x = EarthRadius()*Sin(ACos(z / (EarthRadius()+TOA_alt) ))*Cos(phi);
y = EarthRadius()*Sin(ACos(z / (EarthRadius()+TOA_alt) ))*Sin(phi);
fPhi_local = (TwoPi() - 0)*rndm->Rndm() + 0;
if(thetaMode == "sinus") {
Double_t r1 = 1.;
Double_t r2 = -1.;
fTheta_local = 0.5*ACos(r1 - rndm->Rndm()*(r1-r2));
}
else if(thetaMode == "linear") fTheta_local = (PiOver2() - 0)*rndm->Rndm() + 0;
else if(thetaMode == "FlatHmax") fTheta_local = ACos( Exp((1.89 - 10)*rndm->Rndm() / 7.59));
else Msg(EsafMsg::Panic) << "Wrong argument for theta Mode -> " <<thetaMode << MsgDispatch;
EarthVector TOA_impact_try(x,y,z);
EarthVector v1 = TOA_impact_try.Unit();
EarthVector vrot;
fOmega.SetXYZ(1,0,0);
EarthVector Uz(0,0,1);
vrot = Uz.Cross(v1);
if(vrot.Mag() > kTolerance) {
fOmega.Rotate(v1.Theta(),vrot);
fOmega.Rotate(fPhi_local,v1);
vrot = v1.Cross(fOmega);
fOmega.Rotate(PiOver2()+fTheta_local,vrot);
}
else fOmega.SetMagThetaPhi(1.,Pi() - fTheta_local,fPhi_local + Pi());
fTheta_mes = Pi() - fOmega.Theta();
if(fTheta_mes < fThetaMin || fTheta_mes > fThetaMax) continue;
if(fOmega.Phi() < 0) fPhi_mes = (fOmega.Phi() + TwoPi()) - Pi();
else fPhi_mes = fOmega.Phi() + Pi();
if(fPhi_mes < fPhiMin || fPhi_mes > fPhiMax) continue;
TOA_impact_try.SetZ(TOA_impact_try.Z() - EarthRadius());
fOutOfFoV = kFALSE;
EarthVector intersection, impactASL;
if(IsInFoV(TOA_impact_try)) eventisfound = kTRUE;
else if(!fRejectFakeEvents) {
EarthVector intersection = FoVintersection(TOA_impact_try,fOmega);
EarthVector impactASL = atmo->ImpactASL(TOA_impact_try,fOmega);
intersection = FoVintersection(TOA_impact_try,fOmega);
impactASL = atmo->ImpactASL(TOA_impact_try,fOmega);
if(intersection.Z() == -HUGE || intersection.Z() == HUGE || intersection.IsUnderSeaLevel()) fOutOfFoV = kTRUE;
else if(intersection.Zv() > TOA_alt) fOutOfFoV = kTRUE;
else if((impactASL - TOA_impact_try).Mag() < (intersection - TOA_impact_try).Mag()) fOutOfFoV = kTRUE;
eventisfound = kTRUE;
}
else {
intersection = FoVintersection(TOA_impact_try,fOmega);
impactASL = atmo->ImpactASL(TOA_impact_try,fOmega);
if(intersection.Z() == -HUGE || intersection.Z() == HUGE || intersection.IsUnderSeaLevel()) continue;
else if(intersection.Zv() > TOA_alt) continue;
else if((impactASL - TOA_impact_try).Mag() < (intersection - TOA_impact_try).Mag()) continue;
else eventisfound = kTRUE;
fOutOfFoV = kFALSE;
}
}
fTOAImpact.SetXYZ(x,y,z - EarthRadius());
#ifdef DEBUG
Msg(EsafMsg::Debug) <<"********** Nb of trials = "<<counterr<<MsgDispatch;
#endif
fEarthImpact = atmo->ImpactASL(fTOAImpact,fOmega);
fHitGround = 0;
if(fEarthImpact.Z() == -HUGE) Msg(EsafMsg::Warning) << "<FindImpactOnTOA> ImpactASL says fTOAImpact is under sea level" << MsgDispatch;
else if(fEarthImpact.Z() == HUGE) Msg(EsafMsg::Debug) << "<FindImpactOnTOA> No ImpactASL found" << MsgDispatch;
else fHitGround = 1;
if(fHitGround == 0) {
fEarthImpact = atmo->ImpactAtTOA(fTOAImpact,fOmega);
if(fEarthImpact.Z() == -HUGE) Msg(EsafMsg::Warning) << "<FindImpactOnTOA> ImpactAtTOA says fTOAImpact is under sea level" << MsgDispatch;
else if(fEarthImpact.Z() == HUGE) Msg(EsafMsg::Warning) << "<FindImpactOnTOA> No ImpactAtTOA found. IT SHOULD NOT OCCUR HERE" << MsgDispatch;
}
}
else if(ImpactMode == "ASL") {
Msg(EsafMsg::Warning) << "<FindImpactOnTOA> ASL impact mode : Cosmic Ray LOCAL direction is taken at sea level" << MsgDispatch;
Msg(EsafMsg::Warning) << "<FindImpactOnTOA> ASL impact mode : DOES NOT CHECK IF TRACK CUTS THE FOV" << MsgDispatch;
Int_t infloop = 0;
while(kTRUE) {
if(infloop++ > 100000) {
Msg(EsafMsg::Warning) << "<FindImpactOnTOA> ASL impact mode : Inifinite loop broken by hand" << MsgDispatch;
return kFALSE;
}
Double_t x,y,z;
Double_t Xmin = pConfig->GetNum("GeneratorLightToEuso.ImpactXmin")*km;
Double_t Ymin = pConfig->GetNum("GeneratorLightToEuso.ImpactYmin")*km;
Double_t Xmax = pConfig->GetNum("GeneratorLightToEuso.ImpactXmax")*km;
Double_t Ymax = pConfig->GetNum("GeneratorLightToEuso.ImpactYmax")*km;
x = rndm->Rndm()*(Xmax - Xmin) + Xmin;
y = rndm->Rndm()*(Ymax - Ymin) + Ymin;
z = EarthRadius() - (x*x + y*y)/(2*EarthRadius());
fEarthImpact.SetXYZ(x,y,z);
fEarthImpact.SetZ(fEarthImpact.Z() - EarthRadius());
fHitGround = 1;
fPhi_local = (TwoPi() - 0)*rndm->Rndm() + 0;
if(thetaMode == "sinus") {
Double_t r1 = 1.;
Double_t r2 = -1.;
fTheta_local = 0.5*ACos(r1 - rndm->Rndm()*(r1-r2));
}
else if(thetaMode == "linear") fTheta_local = (PiOver2() - 0)*rndm->Rndm() + 0;
else if(thetaMode == "FlatHmax") fTheta_local = ACos( Exp((1.89 - 10)*rndm->Rndm() / 7.59));
else Msg(EsafMsg::Panic) << "Wrong argument for theta Mode -> " <<thetaMode << MsgDispatch;
EarthVector v1(x,y,z);
v1 = v1.Unit();
EarthVector vrot;
fOmega.SetXYZ(1,0,0);
EarthVector Uz(0,0,1);
vrot = Uz.Cross(v1);
if(vrot.Mag() > kTolerance) {
fOmega.Rotate(v1.Theta(),vrot);
fOmega.Rotate(fPhi_local,v1);
vrot = v1.Cross(fOmega);
fOmega.Rotate(Pi()/2+fTheta_local,vrot);
}
else fOmega.SetMagThetaPhi(1.,Pi() - fTheta_local,fPhi_local + Pi());
fTheta_mes = Pi() - fOmega.Theta();
fOutOfFoV = kFALSE;
if(fOmega.Phi() < 0) fPhi_mes = (fOmega.Phi() + TwoPi()) - Pi();
else fPhi_mes = fOmega.Phi() + Pi();
if(fPhi_mes < fPhiMin || fPhi_mes > fPhiMax || fTheta_mes < fThetaMin || fTheta_mes > fThetaMax) continue;
else break;
}
fTOAImpact = atmo->ImpactAtTOA(fEarthImpact,EarthVector(-fOmega));
if(fTOAImpact.Z() == -HUGE) Msg(EsafMsg::Warning) << "<FindImpactOnTOA> ImpactAtTOA says fEarthImpact is under sea level" << MsgDispatch;
if(fTOAImpact.Z() == HUGE) Msg(EsafMsg::Warning) << "<FindImpactOnTOA> No ImpactAtTOA found : it Should not occur here" << MsgDispatch;
fEarthImpact = atmo->ImpactASL(fTOAImpact,fOmega);
fHitGround = 0;
if(fEarthImpact.Z() == -HUGE) Msg(EsafMsg::Warning) << "<FindImpactOnTOA> ImpactASL says fTOAImpact is under sea level" << MsgDispatch;
else if(fEarthImpact.Z() == HUGE) Msg(EsafMsg::Debug) << "<FindImpactOnTOA> No ImpactASL found" << MsgDispatch;
else fHitGround = 1;
if(fHitGround == 0) {
Msg(EsafMsg::Warning) << "<FindImpactASL> No earth impact found. IT SHOULD NOT OCCUR IN ASL mode" << MsgDispatch;
fEarthImpact = atmo->ImpactAtTOA(fTOAImpact,fOmega);
if(fEarthImpact.Z() == -HUGE) Msg(EsafMsg::Warning) << "<FindImpactOnTOA> ImpactAtTOA says fTOAImpact is under sea level" << MsgDispatch;
else if(fEarthImpact.Z() == HUGE) Msg(EsafMsg::Warning) << "<FindImpactOnTOA> No ImpactAtTOA found. IT SHOULD NOT OCCUR HERE" << MsgDispatch;
}
}
else Msg(EsafMsg::Panic) << "Wrong argument for GeneratorLightToEuso.ImpactMode config parameter" <<MsgDispatch;
return kTRUE;
}
Bool_t SlastShowerGenerator::IsInFoV(const EarthVector& pos) const {
EarthVector Ux(1,0,0);
EarthVector pos_e = pos;
pos_e -= fEusoVector;
pos_e.Rotate(Pi(),Ux);
if(pos_e.Theta() < fFoV) return kTRUE;
else return kFALSE;
}
EarthVector SlastShowerGenerator::FoVintersection(const EarthVector& pos, const EarthVector& dir) const {
Double_t w = Cos(fFoV)*Cos(fFoV);
EarthVector direc = dir.Unit();
EarthVector rtn;
if( (pos.Zv() + kAltitudeTolerance) < 0 ) {
Msg(EsafMsg::Warning) << "<FoVintersection> Starting position is underground, SHOULD NOT OCCUR" << MsgDispatch;
rtn.SetXYZ(0,0,-HUGE);
return rtn;
}
EarthVector Ux(1,0,0);
EarthVector pos_e = pos;
EarthVector direc_e = direc;
pos_e -= fEusoVector;
direc_e.Rotate(Pi(),Ux);
pos_e.Rotate(Pi(),Ux);
Double_t mag(0);
Double_t a = direc_e.Z()*direc_e.Z() - w;
Double_t b = 2*pos_e.Z()*direc_e.Z() - 2*w*pos_e*direc_e;
Double_t c = pos_e.Z()*pos_e.Z() - pos_e.Mag2()*w;
pair<Int_t,Double_t*>& p = findRoots(a,b,c);
if(p.first == 0) {
rtn.SetXYZ(0,0,HUGE);
return rtn;
}
if(p.first == 1) mag = p.second[0];
else if(p.first == 2) {
if((p.second[0]+kAltitudeTolerance) * (p.second[1]+kAltitudeTolerance) < 0) mag = Max(p.second[0],p.second[1]);
else mag = Min(p.second[0],p.second[1]);
}
if(mag < -kAltitudeTolerance) {
rtn.SetXYZ(0,0,HUGE);
return rtn;
}
rtn = pos_e + mag*direc_e;
rtn.Rotate(-Pi(),Ux);
rtn += fEusoVector;
return rtn;
}
Bool_t SlastShowerGenerator::GetFirstPoint(){
Bool_t rtn = kTRUE;
if(fFirstType == "POS") {
fXcurrent = 0;
fXnext = fXcurrent;
fCurrentPoint = fInitPoint;
fNextPoint = fCurrentPoint;
fTimecurrent = 0;
fTimenext = 0;
return rtn;
}
Int_t status = Atmosphere::Get()->InvertGrammage(fTOAImpact,fOmega,fX1,fInitPoint);
if(status == -1) Msg(EsafMsg::Info) << "<GetFirstPoint> InvertGrammage() status is -1, it Should not" << MsgDispatch;
else if(status == 0) {
fXcurrent = 0;
fXnext = fXcurrent;
fCurrentPoint = fInitPoint;
fNextPoint = fCurrentPoint;
fTimecurrent = 0;
fTimenext = 0;
}
else if(status == 1)
Msg(EsafMsg::Warning) << "<GetFirstPoint> InvertGrammage() status is 1, it SHOULD NOT" << MsgDispatch;
else if(status == 2) {
if(Atmosphere::Get()->Grammage(fTOAImpact,Atmosphere::Get()->ImpactAtTOA(fTOAImpact,fOmega)) < fX1+kTolerance)
Msg(EsafMsg::Info) << "<GetFirstPoint> TOA reached : Primary has not interacted in the atmosphere" << MsgDispatch;
else Msg(EsafMsg::Warning) << "<GetFirstPoint> PROBLEM : TOA reached, BUT primary SHOULD have interacted in the atmosphere" << MsgDispatch;
rtn = kFALSE;
}
else if(status == 3) {
if(Atmosphere::Get()->Grammage(fTOAImpact,Atmosphere::Get()->ImpactASL(fTOAImpact,fOmega)) < fX1+kTolerance)
Msg(EsafMsg::Info) << "<GetFirstPoint> Primary has hit the ground without interacting in the atmosphere " << MsgDispatch;
else Msg(EsafMsg::Warning) << "<GetFirstPoint> PROBLEM : ground reached, BUT primary SHOULD have interacted before" << MsgDispatch;
rtn = kFALSE;
}
else Msg(EsafMsg::Warning) << "<GetFirstPoint> InvertGrammage() status unknown : "<< status << MsgDispatch;
return rtn;
}
Bool_t SlastShowerGenerator::IsXmaxInFoV(){
Bool_t rtn = kTRUE;
if (fShowerParametrization != "GIL")
Msg(EsafMsg::Panic) << "<GetShowerParametrization> Only GIL formula available so far" <<MsgDispatch;
Float_t x0 = 37.15*g/cm2;
Float_t Xmax = (1.7 + 0.76*(Log(fEnergy/0.81e8) - Log(fAtomicMass))) * x0 + fX1;
EarthVector Maxpos(1);
Int_t status = Atmosphere::Get()->InvertGrammage(fTOAImpact,fOmega,Xmax,Maxpos);
if(status == -1) Msg(EsafMsg::Warning) << "<GetFirstPoint> InvertGrammage() status is -1, it Should not" << MsgDispatch;
else if(status == 1) Msg(EsafMsg::Warning) << "<GetFirstPoint> InvertGrammage() status is 1, it SHOULD NOT" << MsgDispatch;
else if(status == 2 || status == 3) rtn = kFALSE;
else if(status == 0) rtn = IsInFoV(Maxpos);
else Msg(EsafMsg::Panic) << "<GetFirstPoint> InvertGrammage() status unknown : "<< status << MsgDispatch;
return rtn;
}
Bool_t SlastShowerGenerator::GetNextStep() {
fInCome = kTRUE;
if(fCurrentPoint.IsUnderSeaLevel()) return kFALSE;
fCurrentPoint = fNextPoint;
Int_t status = Atmosphere::Get()->InvertGrammage(fCurrentPoint,fOmega,fDepthStep,fNextPoint);
if(status == -1) {
Msg(EsafMsg::Warning) << "<GetNextStep> InvertGrammage() status is -1, it Should not" << MsgDispatch;
return kFALSE;
}
else if(status == 0) {
fXcurrent = fXnext;
fXnext = fXcurrent + fDepthStep;
fTimecurrent = (fCurrentPoint-fInitPoint).Mag()/Clight();
fTimenext = (fNextPoint-fInitPoint).Mag()/Clight();
}
else if(status == 1) {
Msg(EsafMsg::Warning) << "<GetNextStep> InvertGrammage() status is 1, it Should not" << MsgDispatch;
return kFALSE;
}
else if(status == 2) {
Msg(EsafMsg::Info) << "<GetNextStep> Shower goes out atmosphere : TOA reached " << MsgDispatch;
return kFALSE;
}
else if(status == 3) {
Msg(EsafMsg::Info) << "<GetNextStep> Shower reaches the ground " << MsgDispatch;
return kFALSE;
}
else Msg(EsafMsg::Warning) << "<GetNextStep> InvertGrammage() status unknown : "<< status << MsgDispatch;
return kTRUE;
}
void SlastShowerGenerator::GIL(Double_t x) {
Float_t x0 = 37.15;
Float_t t = x*cm2/g/x0;
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) {
Msg(EsafMsg::Warning)<<"<GIL> Asking for depth <= 0. Zero returned"<<MsgDispatch;
fNe = 0.;
return;
}
fNe = fEnergy/1.45e9*Exp(t-t_max-2*t*Log(fAge));
}
void SlastShowerGenerator::GFA(Double_t x) {
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;
Nmax = Exp(-20.95 + 2.292*Log10(fEnergy));
Xmax = -258.1 + 54.81*Log10(fEnergy);
}
else if(fAtomicMass == 56) {
sigma = 0.470-0.0135*Log10(fEnergy);
p0 = -9.043;
p1 = 0.9923;
Nmax = Exp(-21.68 + 2.329*Log10(fEnergy));
Xmax = -430.7 + 59.08*Log10(fEnergy);
}
else Msg(EsafMsg::Panic)<<"Gaussian parametrization in age not available for primaries A = "<<fAtomicMass<<MsgDispatch;
Xmax = Xmax*g/cm2;
fAge = 3*x/(x+2.*Xmax);
if(fAge<=0) {
Msg(EsafMsg::Warning)<<"<GFA> Asking for depth <= 0. Zero returned"<<MsgDispatch;
fNe = 0.;
return;
}
fNe = 1.11*Nmax * Exp( -(fAge-1)*(fAge-1) / (2.*sigma*sigma) );
return;
}
void SlastShowerGenerator::GHF(Double_t x) {
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 if(fAtomicMass == 56) {
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);
}
else Msg(EsafMsg::Panic)<<"Gaisser-Hillas parametrization not available for primaries A = "<<fAtomicMass<<MsgDispatch;
Xmax = Xmax*g/cm2;
X1 = X1*g/cm2;
X0 = X0*g/cm2;
X=x;
lambda = lambda*g/cm2;
fNe = 1.11*Nmax * Power((X-X0)/(Xmax-X0),(Xmax-X0)/lambda) * Exp((Xmax-X)/lambda);
fAge = 3*x/(x+2.*Xmax);
if(fAge<=0) {
Msg(EsafMsg::Warning)<<"<GHF> Asking for depth <= 0. Zero returned"<<MsgDispatch;
fNe = 0.;
return;
}
return;
}
PhysicsData* SlastShowerGenerator::Get() {
Reset();
//root.cern.ch/root/html/TRandom.html">TRandom *rndm = EsafRandom::Get();
if(fSpectrumType == "powerlaw") {
if ( fEnergySlope == -1.) {
fEnergy = fEnergyMin * Exp(rndm->Rndm() * Log(fEnergyMax/fEnergyMin) );
}
else {
Double_t e1 = Power(fEnergyMin,fEnergySlope + 1);
Double_t e2 = Power(fEnergyMax,fEnergySlope + 1);
fEnergy = Power( e1 + (e2-e1)*rndm->Rndm(),1/(fEnergySlope + 1) );
}
}
else if(fSpectrumType == "GZKHiRes2005") {
fEnergy = -HUGE;
while((fEnergy < fEnergyMin) || (fEnergy > fEnergyMax)) {
Double_t log10E = fSpectrumRdnm->//root.cern.ch/root/html/TH1F.html#TH1F:GetRandom">GetRandom();
fEnergy = Power(10.,log10E);
}
}
else Msg(EsafMsg::Panic)<<"<Get> Wrong argument for spectrum type == "<< fSpectrumType <<MsgDispatch;
if( fFirstType == "POS") {
fInitPoint.SetXYZ(fFirstPointX,fFirstPointY,fFirstPointZ);
fPhi_mes = (fPhiMax - fPhiMin)*rndm->Rndm() + fPhiMin;
fTheta_mes = (fThetaMax - fThetaMin)*rndm->Rndm() + fThetaMin;
fOmega.SetXYZ(-Sin(fTheta_mes)*Cos(fPhi_mes),-Sin(fTheta_mes)*Sin(fPhi_mes),-Cos(fTheta_mes));
Msg(EsafMsg::Info) << "<Get()> Pos option, LOCAL theta and phi set equal to MES values" << MsgDispatch;
fTheta_local = fTheta_mes;
fPhi_local = fPhi_mes;
fTOAImpact = Atmosphere::Get()->ImpactAtTOA(fInitPoint,EarthVector(-fOmega));
fEarthImpact = Atmosphere::Get()->ImpactASL(fInitPoint,fOmega);
fHitGround = 0;
if(fEarthImpact.Z() == -HUGE) Msg(EsafMsg::Warning) << "<Get(), POS option> ImpactASL says fInitPoint is under sea level" << MsgDispatch;
else if(fEarthImpact.Z() == HUGE) Msg(EsafMsg::Debug) << "<Get(), POS option> Shower does not reach the ground" << MsgDispatch;
else fHitGround = 1;
if(fHitGround == 0) {
fEarthImpact = Atmosphere::Get()->ImpactAtTOA(fInitPoint,fOmega);
if(fEarthImpact.Z() == -HUGE) Msg(EsafMsg::Warning) << "<Get(), POS option> ImpactAtTOA says fInitPoint is under sea level" << MsgDispatch;
else if(fEarthImpact.Z() == HUGE) Msg(EsafMsg::Warning) << "<Get(), POS option> No ImpactAtTOA found. IT SHOULD NOT OCCUR HERE" << MsgDispatch;
}
fRejectFakeEvents = kFALSE;
fOutOfFoV = kFALSE;
DevelopShower();
}
else {
Int_t cycle(0);
while (kTRUE) {
if ( cycle++>1000000 ) break;
if ( FindImpactOnTOA() && DevelopShower()) break;
}
}
if (!fQuiet) {
MsgForm(EsafMsg::Info,"ShowerTrack produced : Energy = %.2E MeV, Theta_local = %.2f deg, Phi_local = %.2f deg",
fTrack->GetEnergy(), fTheta_local*RadToDeg(),fPhi_local*RadToDeg());
MsgForm(EsafMsg::Info," Impact on TOA : (%.2f,%.2f,%.2f) km with MES angles : Theta_mes = %.2f, Phi_mes = %.2f",
fTOAImpact.X()/km,fTOAImpact.Y()/km,fTOAImpact.Z()/km,fTheta_mes*RadToDeg(),fPhi_mes*RadToDeg());
MsgForm(EsafMsg::Info," First interaction X1 = %.3f at : (%.2f,%.2f,%.2f) km gives a shower with %d steps",
fTrack->fX1*cm2/g,fTrack->GetInitPos().X()/km,fTrack->GetInitPos().Y()/km,fTrack->GetInitPos().Z()/km, fTrack->GetNumStep());
}
else {
MsgForm(EsafMsg::Debug,"SlastShowerGenerator produced a ShowerTrack: Energy = %.2E MeV Theta = %.2f deg Phi = %.2f deg",
fTrack->GetEnergy(), fTrack->GetTheta()*RadToDeg(),fTrack->GetPhi()*RadToDeg());
MsgForm(EsafMsg::Debug," originated with X1 = %.3f at: (%.2f,%.2f,%.2f) km with %d steps",
fTrack->fX1*cm2/g,fTrack->GetInitPos().X()/km,fTrack->GetInitPos().Y()/km,fTrack->GetInitPos().Z()/km, fTrack->GetNumStep());
}
return fTrack;
}
void SlastShowerGenerator::Reset() {
fHitGround = 0;
fInCome = kFALSE;
fOutCome = kFALSE;
fInitPoint.SetXYZ(0,0,HUGE);
fOmega.SetXYZ(0,0,HUGE);
fCurrentPoint.SetXYZ(0,0,HUGE);
fNextPoint.SetXYZ(0,0,HUGE);
fEarthImpact.SetXYZ(0,0,HUGE);
fTOAImpact.SetXYZ(0,0,HUGE);
fTheta_local = -HUGE;
fPhi_local = -HUGE;
fTheta_mes = -HUGE;
fPhi_mes = -HUGE;
fX1 = -HUGE;
fEnergy = -HUGE;
}
void SlastShowerGenerator::GetShowerInfo() {
fTrack->fEnergy = fEnergy*eV;
fTrack->fTheta = fTheta_mes;
fTrack->fPhi = fPhi_mes;
fTrack->fX1 = fX1;
fTrack->fEthrEl = 0;
fTrack->fDirVers = fOmega;
fTrack->fInitPos = fInitPoint;
fTrack->fHitGround = fHitGround;
fTrack->fEarthImpact = fEarthImpact;
}
Bool_t SlastShowerGenerator::DevelopShower() {
GetX1();
Bool_t has_interacted = GetFirstPoint();
if(!fTrack) fTrack = new ShowerTrack();
else fTrack->Reset();
if(fRejectFakeEvents && (!has_interacted || fOutOfFoV)) return kFALSE;
if(fRejectNoXmax) {
if(!IsXmaxInFoV()) {
Msg(EsafMsg::Info)<<"Shower Xmax is not in FoV --> retry for another event"<<MsgDispatch;
return kFALSE;
}
}
if(!has_interacted) {
fX1 = -HUGE;
fInitPoint.SetXYZ(0,0,HUGE);
GetShowerInfo();
Msg(EsafMsg::Info)<<"Cosmic ray has not interacted in the atmosphere"<<MsgDispatch;
return kTRUE;
}
if(!fRejectFakeEvents && fOutOfFoV) {
GetShowerInfo();
Msg(EsafMsg::Info)<<"ShowerTrack does NOT cut the FoV"<<MsgDispatch;
return kTRUE;
}
Int_t cycle(0);
while(kTRUE) {
if(cycle++>100000) return kFALSE;
if(!GetNextStep()) break;
if( fInCome) {
ShowerStep s = GetShowerStep();
fTrack->Add(s);
}
}
GetShowerInfo();
Msg(EsafMsg::Info)<<"ShowerTrack cut the FoV"<<MsgDispatch;
return kTRUE;
}
MCTruth* SlastShowerGenerator::GetTruth() {
if(!fTruth) fTruth = new MCTruth();
fTruth->SetEnergy(fTrack->GetEnergy());
fTruth->SetThetaPhi(fTrack->GetTheta(),fTrack->GetPhi());
fTruth->SetFirstInt(fInitPoint,fX1);
Int_t particode = Int_t(floor(fAtomicMass+0.5));
fTruth->SetParticle(particode);
fTruth->SetTOAImpact(fTOAImpact);
fTruth->SetEarthImpact(fEarthImpact);
if (fTrack->Size() > 0) {
const ShowerStep& s = fTrack->GetLastStep();
fTruth->SetEarthAge(s.GetAgef());
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;
}
Double_t SlastShowerGenerator::GetXmax() {
Float_t x0 = 37.15*g/cm2;
Float_t fEnergyScale;
if (fEnergy<0) fEnergyScale = 1.e20;
else fEnergyScale = fEnergy;
Float_t t_max = 1.7 + 0.76*(Log(fEnergyScale/0.81e8) - Log(fAtomicMass));
return t_max*x0;
}
ROOT page - Class index - Class Hierarchy - Top of the page
This page has been automatically generated. If you have any comments or suggestions about the page layout send a mail to ROOT support, or contact the developers with any questions or problems regarding ROOT.