// ESAF : Euso Simulation and Analysis Framework
// $Id: SlastLightToEuso.cc,v 1.3 2005/05/03 07:26:35 naumov Exp $
// created Oct, 28 2002
#include "SlastLightToEuso.hh"
#include "EEventTruthAdder.hh"
#include "Config.hh"
#include "ShowerTrack.hh"
#include "EShower.hh"
#include "EShowerFiller.hh"
#include "Atmosphere.hh"
#include "EConst.hh"
ClassImp(SlastLightToEuso)
// interface to fortran functions
extern "C" int slastinit_();
extern "C" int slastsimulation_();
extern "C" {
struct slast2esaf_header {
Float_t Energy;
Float_t Theta;
Float_t Phi;
Float_t X1;
Float_t Rint[3];
Float_t AgeEarth;
Float_t Rimpact[3];
Float_t Rmax[3];
Float_t Xmax;
Float_t ShowerHitGround;
};
struct slast2esaf_shower {
int sn;
Float_t ShowerXi[10000];
Float_t ShowerXf[10000];
Float_t ShowerRxi[10000];
Float_t ShowerRyi[10000];
Float_t ShowerRzi[10000];
Float_t ShowerRxf[10000];
Float_t ShowerRyf[10000];
Float_t ShowerRzf[10000];
Float_t ShowerTimei[10000];
Float_t ShowerTimef[10000];
Float_t ShowerAgei[10000];
Float_t ShowerAgef[10000];
Float_t ShowerNe[10000];
};
struct slast2esaf_fluoresence {
int fn;
Float_t fxshower[500000];
Float_t fyshower[500000];
Float_t fzshower[500000];
Float_t fxfocal[500000];
Float_t fyfocal[500000];
Float_t fzfocal[500000];
Float_t ftheta1[500000];
Float_t fphi1[500000];
Float_t flambda[500000];
Float_t ftime[500000];
};
struct slast2esaf_cherenkov {
int cn;
Float_t cxshower[500000];
Float_t cyshower[500000];
Float_t czshower[500000];
Float_t cxfocal[500000];
Float_t cyfocal[500000];
Float_t czfocal[500000];
Float_t ctheta1[500000];
Float_t cphi1[500000];
Float_t clambda[500000];
Float_t ctime[500000];
};
extern struct slast2esaf_header esafheader_;
extern struct slast2esaf_shower esafshower_;
extern struct slast2esaf_fluoresence esaffluor_;
extern struct slast2esaf_cherenkov esafcherenkov_;
}
// constructor
//_________________________________________________________________________________________
SlastLightToEuso::SlastLightToEuso() : LightToEuso("SLAST"),fPhotons(0),fTruth(0) {
Init();
}
//_________________________________________________________________________________________
void SlastLightToEuso::Init() {
// Enabling the atmosphere
Msg(EsafMsg::Info) << "Atmosphere enabled " << Atmosphere::Get()->GetType() << MsgDispatch;
PrepareDataCards();
slastinit_();
fPhotons = new ListPhotonsOnPupil(new vector<ParentPhoton*>);
}
// destructor
//_________________________________________________________________________________________
SlastLightToEuso::~SlastLightToEuso() {
if ( fTruth )
delete fTruth;
}
// returns the appropriate PhotonsOnPupil object for one event
//_________________________________________________________________________________________
PhotonsOnPupil *SlastLightToEuso::Get(const DetectorGeometry* dg) {
// Slast banner
Msg(EsafMsg::Info) << "SLAST.f77 called" << MsgDispatch;
// call SLAST fortran code
slastsimulation_();
Msg(EsafMsg::Info) << "SLAST.f77 completed succesfully" << MsgDispatch;
// convert data to ListPhotonsOnPupil and return it
fPhotons->Clear();
// loop on photons coming from slast
double x[3];
double y[3];
double th=0.,ph=0.,w=0.,t=0.;
int fn = esaffluor_.fn;
if( esaffluor_.fn == 0 && esafcherenkov_.cn == 0) return fPhotons;
if (fn > 500000) fn = 500000;
// Filling fluorescence information from SLAST
for(int i=0; i<fn; i++)
{
x[0] = esaffluor_.fxshower[i]*m; // in m
x[1] = esaffluor_.fyshower[i]*m; // in m
x[2] = esaffluor_.fzshower[i]*m; // in m
y[0] = esaffluor_.fxfocal[i]*m;
y[1] = esaffluor_.fyfocal[i]*m;
y[2] = esaffluor_.fzfocal[i]*m;
th = esaffluor_.ftheta1[i]*deg;
ph = esaffluor_.fphi1[i]*deg;
w = esaffluor_.flambda[i]*nm ;
t = esaffluor_.ftime[i]*microsecond;
ParentPhoton *p = new ParentPhoton(i,x,y,th,ph,w,t,(ParentPhotonType)1);
fPhotons->Add( p );
}
int cn = esafcherenkov_.cn;
if (cn > 500000) cn = 500000;
// Filling cherenkov information from SLAST
for(int i=0; i<cn; i++)
{
x[0] = esafcherenkov_.cxshower[i]*m;
x[1] = esafcherenkov_.cyshower[i]*m;
x[2] = esafcherenkov_.czshower[i]*m;
y[0] = esafcherenkov_.cxfocal[i]*m;
y[1] = esafcherenkov_.cyfocal[i]*m;
y[2] = esafcherenkov_.czfocal[i]*m;
th = esafcherenkov_.ctheta1[i]*deg;
ph = esafcherenkov_.cphi1[i]*deg;
w = esafcherenkov_.clambda[i]*nm ;
t = esafcherenkov_.ctime[i]*microsecond;
ParentPhoton *p = new ParentPhoton(i+fn,x,y,th,ph,w,t,(ParentPhotonType)2);
fPhotons->Add( p );
}
MsgForm(EsafMsg::Info,"Fluorescent %d Cherenkov %d Total %d photons generated",fn,cn,fn+cn);
// set smontecarlo truth for the shower (clearing previous event first)
SafeDelete(fTruth);
EEvent* event = EEvent::GetCurrent();
if ( event ) {
EEventTruthAdder a( GetTruth());
event->Fill( a );
// Fill the Shower Track Object
ShowerTrack *track = new ShowerTrack();
ShowerStep s;
int sn = esafshower_.sn;
if (sn > 10000) sn = 10000;
for(int i=0;i<sn;i++) {
s.fXi = esafshower_.ShowerXi[i]*g/cm2;
s.fXf = esafshower_.ShowerXf[i]*g/cm2;
s.fXYZi.SetXYZ(esafshower_.ShowerRxi[i]*m,esafshower_.ShowerRyi[i]*m,esafshower_.ShowerRzi[i]*m);
s.fXYZf.SetXYZ(esafshower_.ShowerRxf[i]*m,esafshower_.ShowerRyf[i]*m,esafshower_.ShowerRzf[i]*m);
s.fTimei = esafshower_.ShowerTimei[i]*microsecond;
s.fTimef = esafshower_.ShowerTimef[i]*microsecond;
s.fAgei = esafshower_.ShowerAgei[i];
s.fAgef = esafshower_.ShowerAgef[i];
s.fNelectrons = esafshower_.ShowerNe[i];
track->Add(s);
}
track->fEnergy = esafheader_.Energy;
track->fTheta = esafheader_.Theta;
track->fPhi = esafheader_.Phi;
track->fX1 = esafheader_.X1*g/cm2;
track->fEthrEl = 0;
track->fDirVers.SetXYZ(TMath::Sin(esafheader_.Theta)*TMath::Cos(esafheader_.Phi),
TMath::Sin(esafheader_.Theta)*TMath::Sin(esafheader_.Phi),
TMath::Cos(esafheader_.Theta));
track->fInitPos.SetXYZ(esafheader_.Rint[0]*m,esafheader_.Rint[1]*m,esafheader_.Rint[2]*m);
track->fHitGround = (esafheader_.ShowerHitGround == 1 ? kTRUE: kFALSE);
EShowerFiller f( track );
event->Fill( f );
delete track;
}
system("rm -rf TAPE*");
return fPhotons;
}
// re-configure SLAST
//_________________________________________________________________________________________
void SlastLightToEuso::Configure() {
}
// returns Montecarlo shower truth informations for this event
//_________________________________________________________________________________________
MCTruth* SlastLightToEuso::GetTruth() {
if (!fTruth) {
fTruth = new MCTruth();
fTruth->SetLightToEuso( this );
fTruth->SetEnergy( esafheader_.Energy*eV);
fTruth->SetThetaPhi(esafheader_.Theta,esafheader_.Phi);
EVector v;
v[0] = esafheader_.Rint[0]*m;
v[1] = esafheader_.Rint[1]*m;
v[2] = esafheader_.Rint[2]*m;
fTruth->SetFirstInt(v,esafheader_.X1*g/cm2);
v[0] = esafheader_.Rimpact[0]*m;
v[1] = esafheader_.Rimpact[1]*m;
v[2] = esafheader_.Rimpact[2]*m;
fTruth->SetEarthImpact(v,esafheader_.AgeEarth);
v[0] = esafheader_.Rmax[0]*m;
v[1] = esafheader_.Rmax[1]*m;
v[2] = esafheader_.Rmax[2]*m;
fTruth->SetShowerMax(v,esafheader_.Xmax*g/cm2);
}
return fTruth;
}
// returns simulated shower additional infos for this event
//_________________________________________________________________________________________
PhysicsData* SlastLightToEuso::GetPhysics() {
return NULL;
}
// interface to fortran code
//_________________________________________________________________________________________
void SlastLightToEuso::PrepareDataCards() {
//
// Prepare data card to read by slast77 engine
//
ConfigFileParser *pConfig = Config::Get()->GetCF("General","Euso");
Float_t EarthRadius = EConst::EarthRadius()/km;
Float_t ISSVectorX = 0;
Float_t ISSVectorY = 0;
Float_t ISSVectorZ = pConfig->GetNum("Euso.fAltitude");
Float_t EntrancePupilDiameter = 2*pConfig->GetNum("Euso.fRadius")*mm/m;
// reading common to slast77 and slast++ options from GeneratorLightToEuso file
pConfig = Config::Get()->GetCF("LightToEuso","GeneratorLightToEuso");
Float_t FOV = pConfig->GetNum("GeneratorLightToEuso.FoV");
Float_t InteractionVectorX = pConfig->GetNum("GeneratorLightToEuso.InteractionVectorX");
Float_t InteractionVectorY = pConfig->GetNum("GeneratorLightToEuso.InteractionVectorY");
Float_t InteractionVectorZ = pConfig->GetNum("GeneratorLightToEuso.InteractionVectorZ");
Float_t ThetaRangeMin = pConfig->GetNum("GeneratorLightToEuso.ThetaRangeMin");
Float_t ThetaRangeMax = pConfig->GetNum("GeneratorLightToEuso.ThetaRangeMax");
Float_t PhiRangeMin = pConfig->GetNum("GeneratorLightToEuso.PhiRangeMin");
Float_t PhiRangeMax = pConfig->GetNum("GeneratorLightToEuso.PhiRangeMax");
Float_t EnergyRangeMin = pConfig->GetNum("GeneratorLightToEuso.EnergyRangeMin");
Float_t EnergyRangeMax = pConfig->GetNum("GeneratorLightToEuso.EnergyRangeMax");
Float_t EnergySlope = pConfig->GetNum("GeneratorLightToEuso.EnergySlope");
string fType = pConfig->GetStr("GeneratorLightToEuso.UhecrType");
Float_t UhecrType = EConst::AtomicMass(fType);
// reading specific to slast77 options from its own config
Float_t WaveRangeMin = Conf()->GetNum("SlastLightToEuso.WaveRangeMin");
Float_t WaveRangeMax = Conf()->GetNum("SlastLightToEuso.WaveRangeMax");
string DoCherenkov = Conf()->GetStr("SlastLightToEuso.DoCherenkov");
string DoFluorescence = Conf()->GetStr("SlastLightToEuso.DoFluorescence");
string AtmType = Conf()->GetStr("SlastLightToEuso.AtmosphericType");
Float_t AtmTemperature = Conf()->GetNum("SlastLightToEuso.AtmTemperature");
Float_t Albedo = Conf()->GetNum("SlastLightToEuso.Albedo");
Float_t GTU = Conf()->GetNum("SlastLightToEuso.GTU");
string AtmCurvature = Conf()->GetStr("SlastLightToEuso.AtmCurvature");
string EnergyDistributionParametrization = Conf()->GetStr("SlastLightToEuso.EnergyDistributionParametrization");
string ShowerParametrization = Conf()->GetStr("SlastLightToEuso.ShowerParametrization");
FILE *datacard = fopen("auxilar/global.datacard","w");
int doch=0,dofl=0,atm_type=2,curvature=2,energydistribution=1,showerparametrization=1;
fprintf(datacard,"LISTn");
fprintf(datacard,"%s%fn","ERAD ",EarthRadius);
fprintf(datacard,"%s%f %f %fn","ISSF ",ISSVectorX,ISSVectorY,ISSVectorZ);
fprintf(datacard,"%s%fn","FOV ",FOV);
fprintf(datacard,"%s%fn","SDIA ",EntrancePupilDiameter);
fprintf(datacard,"%s%f %fn","WAVE ",WaveRangeMin,WaveRangeMax);
if(DoCherenkov == "yes") doch=1;
if(DoFluorescence == "yes") dofl=1;
fprintf(datacard,"%s%dn","DOCH ",doch);
fprintf(datacard,"%s%dn","DOFL ",dofl);
if(AtmType == "Isothermic") atm_type=1;
if(AtmType == "USStandard") atm_type=2;
fprintf(datacard,"%s%dn","ATMO ",atm_type);
fprintf(datacard,"%s%fn","TEMP ",AtmTemperature);
fprintf(datacard,"%s%fn","ALBE ",Albedo);
fprintf(datacard,"%s%fn","GTU ",GTU);
if(AtmCurvature == "Curved") curvature = 1;
if(AtmCurvature == "Planar") curvature = 2;
fprintf(datacard,"%s%dn","CURV ",curvature);
fprintf(datacard,"%s%f %f %fn","RINT ",InteractionVectorX,InteractionVectorY,InteractionVectorZ);
fprintf(datacard,"%s%f %fn","THET ",ThetaRangeMin,ThetaRangeMax);
fprintf(datacard,"%s%f %fn","PHI ",PhiRangeMin,PhiRangeMax);
fprintf(datacard,"%s%e %en","EINT ",EnergyRangeMin,EnergyRangeMax);
fprintf(datacard,"%s%dn","ERAN ",EnergySlope);
fprintf(datacard,"%s%fn","TYPE ",UhecrType);
if(EnergyDistributionParametrization == "Hillas") energydistribution=1;
fprintf(datacard,"%s%dn","DIST ",energydistribution);
if(ShowerParametrization == "GIL") showerparametrization=1;
fprintf(datacard,"%s%dn","SHOW ",showerparametrization);
fclose(datacard);
}
