// $Id: CorsikaFileShowerGenerator.cc,v 1.7 2005/01/27 10:23:57 dattoli Exp $
// Author: M.Dattoli 2004/08/31
/*****************************************************************************
* ESAF: Euso Simulation and Analysis Framework *
* *
* Id: CorsikaFileShowerGenerator *
* Package: ShowerGenerators *
* Coordinator: Sergio.Bottai, Dmitry.Naumov *
* *
*****************************************************************************/
//_____________________________________________________________________________
// Corsika File Shower Generator
// =============================
//
// <extensive description>
//
// Config file parameters
// ======================
//
// FlagThinning [string]: flag for thinning option;
// can be t (thinning enabled)
// or f (thinning disabled)
//
// HowManyShowers [int] : number of showers simulated in corsika run
//
// ParticleFileName [string]: name of corsika normal particle output file
//
// LongFileName [string]: name of corsika longitudinal output file
//
#include "CorsikaFileShowerGenerator.hh"
#include "MCTruth.hh"
#include <ctime>
#include <cstdlib>
#include <cmath>
#include "PhysicalConstants.hh"
#include "SystemOfUnits.hh"
ClassImp(CorsikaFileShowerGenerator)
//______________________________________________________________________________
CorsikaFileShowerGenerator::CorsikaFileShowerGenerator():
EventGenerator("Corsikafile"),EsafMsgSource(), fTrack(0),
fParticleFile(NULL), fLongFile(NULL), fTruth(0)
//
// Constructor
//
{
fParticleFileName = "";
fLongFileName = "";
}
//______________________________________________________________________________
CorsikaFileShowerGenerator::~CorsikaFileShowerGenerator()
//
// Destructor
//
{
if ( fTrack ) delete fTrack;
if ( fTruth ) delete fTruth;
fTrack=NULL;
fTruth=NULL;
}
//______________________________________________________________________________
void CorsikaFileShowerGenerator::ReplaceInputFile(const char* pfn,const char* lfn)
//
// Replace input files
//
{
// Particle file
if ( fParticleFile )
Msg(EsafMsg::Panic)<<"Can not change name when particle file is already open!"<<MsgDispatch;
if (!pfn)
{
Msg(EsafMsg::Panic)<<" Particle file name ignored!" <<MsgDispatch;
return;
}
Conf()->ReplaceStr("CorsikaFileShowerGenerator.ParticleFileName",pfn);
Msg(EsafMsg::Info)<<"Input file name changed to" << pfn <<MsgDispatch;
//Longitudinal file: info about energy loss, n electrons, ncharged, ncherenkov, depth
if ( fLongFile )
Msg(EsafMsg::Panic)<<"Can not change name when particle file is already open!"<<MsgDispatch;
if (!lfn)
{
Msg(EsafMsg::Panic)<<" Particle file name ignored!" <<MsgDispatch;
return;
}
Conf()->ReplaceStr("CorsikaFileShowerGenerator.LongFileName",lfn);
Msg(EsafMsg::Info)<<"Input file name changed to" << lfn <<MsgDispatch;
}
//______________________________________________________________________________
void CorsikaFileShowerGenerator::Open()
//
// Open file
//
{
//particle file
if ( Conf()->IsDefined( "CorsikaFileShowerGenerator.ParticleFileName" ) )
fParticleFileName = Conf()->GetStr("CorsikaFileShowerGenerator.ParticleFileName");
else
Msg(EsafMsg::Panic)<<"Unknown input particle file name in CorsikaFileShowerGenerator"<<MsgDispatch;
fParticleFile = fopen( fParticleFileName.c_str(), "rb" );
if( fParticleFile == NULL )
Msg(EsafMsg::Panic)<<"Can't open file "+fParticleFileName << MsgDispatch;
else
Msg(EsafMsg::Info) << "File: " << fParticleFileName << " successfully opened " << MsgDispatch;
//longitudinal file
if ( Conf()->IsDefined( "CorsikaFileShowerGenerator.LongFileName" ) )
fLongFileName = Conf()->GetStr("CorsikaFileShowerGenerator.LongFileName");
else
Msg(EsafMsg::Panic)<<"Unknown input longitudinal ditribution file name in CorsikaFileShowerGenerator"<<MsgDispatch;
fLongFile = fopen( fLongFileName.c_str(), "rb" );
if( fLongFile == NULL )
Msg(EsafMsg::Panic)<<"Can't open file "+fLongFileName << MsgDispatch;
else
Msg(EsafMsg::Info) << "File: " << fLongFileName << " successfully opened " << MsgDispatch;
}
//______________________________________________________________________________
void CorsikaFileShowerGenerator::Close()
//
// Close file
//
{
if( fParticleFile ) fclose(fParticleFile ) ;
if( fLongFile ) fclose(fLongFile ) ;
}
//______________________________________________________________________________
PhysicsData* CorsikaFileShowerGenerator::Get()
{
//
// Get a new event
//
// build and return a track
if ( !fParticleFile || !fLongFile) Open();
if ( Conf()->IsDefined( "CorsikaFileShowerGenerator.FlagThinning" ) )
thinning = Conf()->GetStr("CorsikaFileShowerGenerator.FlagThinning");
else Msg(EsafMsg::Panic) <<"thinning flag missing in config file" << MsgDispatch;
if (thinning=="f") {
n1 = 272;
nmax = 7;
Msg(EsafMsg::Info)<<"Thinning disabled"<<MsgDispatch;
}
else if (thinning=="t") {
n1 = 311;
nmax = 8;
Msg(EsafMsg::Info)<<"Thinning enabled"<<MsgDispatch;
}
n2 = n1 + 1;
ConfigFileParser *pConfig = Config::Get()->GetCF("Generators","CorsikaFileShowerGenerator");
howmany = pConfig->GetNum("CorsikaFileShowerGenerator.HowManyShowers");
//random choice of event to be processed among simulated events
srand(time(0));
event =(rand()%(Int_t)(howmany))+1;
Msg(EsafMsg::Debug) << "Processing event n : " << event << MsgDispatch;
if ( LoadTrack(event) ) {
Msg(EsafMsg::Debug) << "Load Track Event" << MsgDispatch;
return fTrack;
}
else {
Msg(EsafMsg::Panic) << "End of file " << fParticleFileName << " reached? " << MsgDispatch;
return (ShowerTrack*)0;
}
return fTrack;
}
//______________________________________________________________________________
void CorsikaFileShowerGenerator::ReadParticleHeader()
{
//
// Read particle file header.
//
Msg(EsafMsg::Debug) << "Reading particle header.." << MsgDispatch;
while(!feof(fParticleFile))
{
//Reading BlockHeader (273,21) or (312,21)
nin = fread(&r,4,1,fParticleFile);
for (k=1; k<=21;k++)
{
pos = ftell(fParticleFile);
//Reading SubBlockHeader
nin = fread(&header,sizeof(header),1,fParticleFile);
//Reading data in SubBlock; Storing them into map
if (header[0]=='R'&&header[3]=='H')
{
float data_runh[n1];
for (i=0; i<n1; i++) {data_runh[i]=0;}
//Reading RUNH Block (once for run)
nin = fread (&data_runh, sizeof(data_runh), 1, fParticleFile);
for (i=0; i<=4; i++) {
aatm[i]=data_runh[253+i];
batm[i]=data_runh[258+i];
catm[i]=data_runh[263+i];
}
} // end of if "runh"
else if (header[0]=='R'&&header[3]=='E')
{
float data_rune[n1];
for (i=0; i<n1; i++) {data_rune[i]=0;}
//reading RUNE Block (once per run)
nin = fread (&data_rune, sizeof(data_rune), 1,fParticleFile );
} // end of if "rune"
else if (header[0]=='E'&&header[3]=='H')
{
float data_evth[n1];
for (i=0; i<n1; i++) {data_evth[i]=0;}
//Reading EVTH Block (1 for event)
nin = fread (&data_evth, sizeof(data_evth), 1,fParticleFile );
ev_num= (Int_t)(data_evth[0]);
if (ev_num==event)
{
Msg(EsafMsg::Debug) << "Reading event n " << ev_num << MsgDispatch;
z_firstint = data_evth[5]; // cm
Msg(EsafMsg::Debug)<< "z_firstint: " << z_firstint <<MsgDispatch;
thick_firstint = GetThickness(abs(z_firstint)); //g/cm^2
Msg(EsafMsg::Debug)<< "thick_firstint: " << thick_firstint <<MsgDispatch;
epolangd = data_evth[9]*180/pi;
eaziangd = data_evth[10]*180/pi;
eldircos = sin(data_evth[9])*cos(data_evth[10]);
emdircos = sin(data_evth[9])*sin(data_evth[10]);
endircos = cos(data_evth[9]);
energy = data_evth[2];
Msg (EsafMsg::Debug) << "epolangd: " << epolangd
<< "n eaziangd: " << eaziangd
<< "n eldircos: " << eldircos
<< "n emdircos: " << emdircos
<< "n endircos: " << endircos
<< "n energy: " << energy << MsgDispatch;
} //end of if (ev_num==event)
} //end of if "evth"
else if (header[0]=='E'&&header[3]=='E')
{
float data_evte[n1];
for (i=0; i<n1; i++) {data_evte[i]=0;}
//Reading EVTE Block (1 for event)
nin = fread (&data_evte, sizeof(data_evte), 1, fParticleFile);
} // end of if "evte"
else if (header[0]=='L'&&header[3]=='G')
{
float data_long[n1];
for (i=0; i<n1; i++) {data_long[i]=0;}
//Reading LONG Block
nin = fread (&data_long, sizeof(data_long), 1, fParticleFile);
} // end of if "long"
else
{
fseek(fParticleFile,pos,SEEK_SET);
float data_block[n2];
for (i=0; i<n2; i++) {data_block[i]=0;}
//Reading DATA BLOCK
nin = fread (&data_block, sizeof(data_block), 1,fParticleFile );
} // end of else... DATA BLOCK
} // end of for (k=1; k<=21;k++)
k=1;
//Reading BlockEnd (273,21) o (312,21)
nin = fread(&r,4,1,fParticleFile);
}
}
//______________________________________________________________________________
void CorsikaFileShowerGenerator::ReadLong()
//
// Read long file
//
{
Msg(EsafMsg::Debug) << "Reading Long File.." << MsgDispatch;
for ( i=0; i<10; i++)
{ en[i]=0;
part[i]=0;
}
nshow=0;
pos = ftell(fLongFile);
while (!feof(fLongFile) && nshow<=event)
{
fscanf(fLongFile,"%*s%*s%*s%f%*s%*s%*s%f%*s%*s%*s%i",&nstep,&step_long,&nshow);
fgets(longitudinal,sizeof(longitudinal),fLongFile);
pos=ftell(fLongFile);
Msg(EsafMsg::Debug)<< "nstep, step_long: " << nstep << "; t"<< step_long << MsgDispatch;
fgets(longitudinal,sizeof(longitudinal),fLongFile);
if (nshow!=event)
{
for (j=1; j<= 2*nstep+12; j++ ){
fgets(longitudinal,sizeof(longitudinal),fLongFile);
pos = ftell(fLongFile);
fseek(fLongFile,pos,SEEK_SET);
}
fseek(fLongFile,pos,SEEK_SET);
nshow = nshow+1;
}
else
{
for ( j=1; j <= nstep; j++)
{
for ( i=0; i<10; i++)
{
fscanf(fLongFile,"%f",&part[i]);
pos=ftell(fLongFile);
fseek(fLongFile,pos,SEEK_SET);
} // reading 1 raw
particle.insert(pair<int,float>(j,part[0]));
particle.insert(pair<int,float>(j,part[0] - step_long));
particle.insert(pair<int,float>(j,part[3] + part[2]));
particle.insert(pair<int,float>(j,part[7]));
particle.insert(pair<int,float>(j,part[9]));
} // reading block
fgets(longitudinal,sizeof(longitudinal),fLongFile);
fgets(longitudinal,sizeof(longitudinal),fLongFile);
fgets(longitudinal,sizeof(longitudinal),fLongFile);
for ( j=1; j <=nstep; j++)
{
for ( i=0; i<10; i++)
{ fscanf(fLongFile,"%f",&en[i]);
} // reading 1 raw
particle.insert(pair<int,float>(j,en[2]));
} // reading block
fgets(longitudinal,sizeof(longitudinal),fLongFile);
fgets(longitudinal,sizeof(longitudinal),fLongFile);
fgets(longitudinal,sizeof(longitudinal),fLongFile);
pos=ftell(fLongFile);
fseek(fLongFile,pos,SEEK_SET);
fscanf(fLongFile,"%*s%*s%f%f%f%*f%*f%*f",&Nemax,&to,&tmax);
for (int i = 1; i<=8; i++) {
fgets(longitudinal,sizeof(longitudinal),fLongFile);
}
nshow=nshow+1;
} // end of else
pos=ftell(fLongFile);
}
}
//______________________________________________________________________________
bool CorsikaFileShowerGenerator::LoadHeader()
//
// Fill header object
//
{
fHeader.clear();
ReadParticleHeader();
x0comes = xi = 0;
y0comes = yi = 0;
if (z_firstint>=0) {
z0comes = z_firstint;}
else {
z0comes = - z_firstint;}
zmax = GetAltitude(tmax);
vertical_distance = (abs(z_firstint) - zmax);
coeff_f = GetCoeff(vertical_distance);
xmax = GetXstep(coeff_f);
ymax = GetYstep(coeff_f);
float coeff = GetCoeff(z0comes);
fHeader.insert(pair<string,float>("EPRMENRG",energy));
fHeader.insert(pair<string,float>("EPRMDPTH",thick_firstint));
fHeader.insert(pair<string,float>("EX0COMES",x0comes));
fHeader.insert(pair<string,float>("EY0COMES",y0comes));
fHeader.insert(pair<string,float>("EZ0COMES",z0comes));
fHeader.insert(pair<string,float>("EPOLANGD",epolangd));
fHeader.insert(pair<string,float>("EAZIANGD",eaziangd));
fHeader.insert(pair<string,float>("ELDIRCOS",eldircos));
fHeader.insert(pair<string,float>("EMDIRCOS",emdircos));
fHeader.insert(pair<string,float>("ENDIRCOS",endircos));
fHeader.insert(pair<string,float>("EXDMPPOS",GetXstep(coeff)));
fHeader.insert(pair<string,float>("EYDMPPOS",GetYstep(coeff)));
fHeader.insert(pair<string,float>("EZDMPPOS",0.));
fHeader.insert(pair<string,float>("EXMAXPOS",xmax));
fHeader.insert(pair<string,float>("EYMAXPOS",ymax));
fHeader.insert(pair<string,float>("EZMAXPOS",zmax));
return true;
}
//______________________________________________________________________________
bool CorsikaFileShowerGenerator::LoadTrack(Int_t event)
//
// Load one event from file
//
{
ReadLong();
if (!LoadHeader()) return false;
fTrack = new ShowerTrack();
fTrack->fEthrEl = 0.;
fTrack->fEnergy = fHeader["EPRMENRG"]*GeV; //in eV
fTrack->fTheta = fHeader["EPOLANGD"]*deg;
fTrack->fPhi = fHeader["EAZIANGD"]*deg;
fTrack->fX1 = fHeader["EPRMDPTH"]*g/cm2;
fTrack->fDirVers[0] = fHeader["ELDIRCOS"];
fTrack->fDirVers[1] = fHeader["EMDIRCOS"];
fTrack->fDirVers[2] = fHeader["ENDIRCOS"];
fTrack->fInitPos[0] = fHeader["EX0COMES"]*cm;
fTrack->fInitPos[1] = fHeader["EY0COMES"]*cm;
fTrack->fInitPos[2] = fHeader["EZ0COMES"]*cm;
if(fHeader["HITGROUN"]==0)
{
fTrack->fHitGround = false;
}
else
{
fTrack->fHitGround = true;
}
coeff_i = 0;
xi=0;
yi=0;
zi = abs(z_firstint);
timei = 0;
for(j=1; j<=nstep; j++)
{
int count=0;
for (pnum=particle.lower_bound(j); pnum!=particle.upper_bound(j); ++pnum)
{
step[count] = pnum->second;
count++;
}
// step_bound (g/cm^2)
fStep.fXi = (step[1]+ thick_firstint)*gram/cm2;
fStep.fXf = (step[1]+ thick_firstint+step_long)*gram/cm2;
// step_bound (cartesian coordinates)
thickness = step_long*endircos;
zf = GetAltitude(step[1]+ thick_firstint+thickness);
vertical_distance = abs(zi-zf);
coeff_f = GetCoeff(vertical_distance);
xf = GetXstep(coeff_f);
yf = GetYstep(coeff_f);
if (zi>0) {
fStep.fXYZi[2] = zi*cm;}
else {
fStep.fXYZi[2] = 0;}
if (zf>0) {
fStep.fXYZf[2] = zf*cm;}
else {
fStep.fXYZf[2] = 0;}
if (fStep.fXYZi[2]!=0 || fStep.fXYZf[2]!=0 )
{
// fStep.fXYZf[2] = zf;
fStep.fXYZi[0] = xi*cm;
fStep.fXYZf[0] = xf*cm;
fStep.fXYZi[1] = yi*cm;
fStep.fXYZf[1] = yf*cm;
// time
distance = sqrt(abs(pow(xi,2) + pow(yi,2) + pow(zi,2) - (pow(xf,2) + pow(yf,2) + pow(zf,2))));
fStep.fTimei = timei*nanosecond;
timef = GetStepTime(distance) + timei;
fStep.fTimef = timef*nanosecond;
fStep.fNelectrons = step[2];
fStep.fNcharged = step[3];
fStep.fNcherenkov = step[4];
fStep.fEloss = step[5];
// age
fStep.fAgei = GetAge(step[1] + thick_firstint);
fStep.fAgef = GetAge(step[1] + thick_firstint+step_long);
fTrack->Add(fStep);
coeff_i = coeff_f;
xi = xf;
yi = yf;
zi = zf;
timei = timef;
}
}
return true;
}
//______________________________________________________________________________
MCTruth* CorsikaFileShowerGenerator::GetTruth()
//
// Get pointer to MonteCarlo Truth
//
{
if (!fTrack->Size()) return NULL;
if (!fTruth) {
fTruth = new MCTruth();}
fTruth->SetEnergy(fHeader["EPRMENRG"]*GeV); //in eV
fTruth->SetThetaPhi(fHeader["EPOLANGD"]*deg,fHeader["EAZIANGD"]*deg);
EVector v;
v[0] = fHeader["EX0COMES"]*cm;
v[1] = fHeader["EY0COMES"]*cm;
v[2] = fHeader["EZ0COMES"]*cm;
fTruth->SetFirstInt(v, fHeader["EPRMDPTH"]);
v[0] = fHeader["EXDMPPOS"]*cm;
v[1] = fHeader["EYDMPPOS"]*cm;
v[2] = fHeader["EZDMPPOS"]*cm;
// in case no value for ShowerMax is present, put arbitrary number 0
v[0] = 0;
v[1] = 0;
v[2] = 0;
map<string,double>::iterator showmax;
showmax = fHeader.find(string("EXMAXPOS"));
if (showmax != fHeader.end()) v[0] = fHeader["EXMAXPOS"]*cm;
showmax = fHeader.find(string("EYMAXPOS"));
if (showmax != fHeader.end()) v[1] = fHeader["EYMAXPOS"]*cm;
showmax = fHeader.find(string("EZMAXPOS"));
if (showmax != fHeader.end()) v[2] = fHeader["EZMAXPOS"]*cm;
fTruth->SetShowerMax(v,fHeader["ESHOWMAX"]);
return fTruth;
}
//______________________________________________________________________________
float CorsikaFileShowerGenerator::GetAltitude(float thick)
{
float h;
if (thick >= GetThickness(4.e+05)) {
h = -catm[0]*log((thick - aatm[0])/batm[0]);}
if (thick >= GetThickness(10.e+05) && thick < GetThickness(4.e+05)) {
h = -catm[1]*log((thick - aatm[1])/batm[1]);}
if (thick >= GetThickness(40.e+05) && thick< GetThickness(10.e+05)) {
h = -catm[2]*log((thick - aatm[2])/batm[2]);}
if (thick >= GetThickness(100.e+05) && thick< GetThickness(40.e+05)) {
h = -catm[3]*log((thick - aatm[3])/batm[3]);}
if (thick <= GetThickness(100.e+05) ) {
h = (aatm[4]-thick)*catm[4]/batm[4];}
return h;
}
//______________________________________________________________________________
float CorsikaFileShowerGenerator::GetThickness(float z1)
{
float thick_f;
float z = abs(z1);
if (z<= 4.e+05) {
thick_f = aatm[0]+batm[0]*exp(-z/catm[0]);}
else if (z> 4.e+05 && z<= 10.e+05) {
thick_f = aatm[1]+batm[1]*exp(-z/catm[1]);}
else if (z> 10.e+05 && z<= 40.e+05) {
thick_f = aatm[2]+batm[2]*exp(-z/catm[2]);}
else if (z> 40.e+05 && z<= 100.e+05) {
thick_f = aatm[3]+batm[3]*exp(-z/catm[3]);}
else if (z> 100.e+05) {
thick_f = aatm[4]-batm[4]*z/catm[4];}
return thick_f;
}
//______________________________________________________________________________
float CorsikaFileShowerGenerator::GetCoeff(float h)
{
float coeff = h/endircos;
return coeff;
}
//______________________________________________________________________________
float CorsikaFileShowerGenerator::GetXstep(float t)
{
float x = xi + t*eldircos;
return x;
}
//______________________________________________________________________________
float CorsikaFileShowerGenerator::GetYstep(float t)
{
float y = yi + t*emdircos;
return y;
}
//______________________________________________________________________________
float CorsikaFileShowerGenerator::GetStepTime(float h)
{
float time = 0.1*h/(c_light);
return time;
}
//______________________________________________________________________________
float CorsikaFileShowerGenerator::GetAge(float h)
{
float age;
float t = h-thick_firstint;
age = 3*t/(t + 2*tmax);
return age;
}