// ESAF : Euso Simulation and Analysis Framework
// $Id: UnisimFileShowerGenerator.cc,v 1.12 2005/05/08 22:28:21 thea Exp $
// sergio bottai created Jan, 22 2004
#include "UnisimFileShowerGenerator.hh"
#include "SystemOfUnits.hh"
#include "MCTruth.hh"
#include "ShowerTrack.hh"
#include <zlib.h>
ClassImp(UnisimFileShowerGenerator)
//______________________________________________________________________________
UnisimFileShowerGenerator::UnisimFileShowerGenerator() :
EventGenerator("Unisimfile"), fTrack(0),fFile(NULL), fTruth(0) {
//
// Construtor
//
fFileName = "";
fFirstEvent = 0;
fCurrentEvent = 0;
}
//______________________________________________________________________________
UnisimFileShowerGenerator::~UnisimFileShowerGenerator() {
//
// Destructor
//
if ( fTrack ) {
delete fTrack;
fTrack=NULL;
}
if ( fTruth ) {
delete fTruth;
fTruth=NULL;
}
}
//______________________________________________________________________________
void UnisimFileShowerGenerator::ReplaceInputFile(const char* fn) {
//
//
//
if ( fFile )
throw runtime_error("UnisimFileShowerGenerator: Cannot change name when file is already open!");
if (!fn) {
cerr << "UnisimFileShowerGenerator: NULL file name ignoredn";
return;
}
Conf()->ReplaceStr("UnisimFileShowerGenerator.FileName",fn);
cout << "Input file name changed to " << fn << endl;
}
//______________________________________________________________________________
void UnisimFileShowerGenerator::Open() {
//
// Open file
//
if ( Conf()->IsDefined( "UnisimFileShowerGenerator.FileName" ) ) {
fFileName = Conf()->GetStr("UnisimFileShowerGenerator.FileName");
fFirstEvent = (Int_t)Conf()->GetNum("UnisimFileShowerGenerator.FirstEvent") ;
}
else
throw runtime_error("Unknown input file name in UnisimFileShowerGenerator");
fFile = (FILE*)gzopen( fFileName.c_str(), "rb" );
if( fFile == NULL )
throw runtime_error("Can't open file "+fFileName);
else
cout << "File: " << fFileName << " successfully opened " << endl;
}
//______________________________________________________________________________
void UnisimFileShowerGenerator::Close() {
//
// Close file
//
if( fFile ) gzclose( (void*)fFile ) ;
}
//______________________________________________________________________________
PhysicsData* UnisimFileShowerGenerator::Get() {
//
// Get a new event
//
// build and return a track
if ( !fFile ) Open();
Int_t event=0;
while (1) {
if ( fTrack ) {
delete fTrack;
fTrack=NULL;
}
if ( fTruth ) {
delete fTruth;
fTruth=NULL;
}
fCurrentEvent++;
if (fCurrentEvent >= fFirstEvent)
{
if ( LoadTrack(event) ){
return fTrack;
}
else {
cout << "UnisimFileShowerGenerator: End of file " << fFileName << " reachedn ? ";
return (ShowerTrack*)0;
}
}
else
{
if ( !LoadTrack(event) ){
cout << "UnisimFileShowerGenerator: End of file " << fFileName << " reachedn ? ";
}
}
}
return fTrack;
}
//______________________________________________________________________________
MCTruth* UnisimFileShowerGenerator::GetTruth() {
//
// get pointer to MonteCarlo Truth
//
if ( fTrack ) {
return fTruth;
}
else { cout << "Event to be loaded with LoadTrack before GetTruth" << endl;
return (MCTruth*)0;
}
}
//______________________________________________________________________________
bool UnisimFileShowerGenerator::LoadTrack( Int_t event ) {
//
// Load one event from file
//
if ( !LoadHeader() ) return false;
char s[5010];
fTruth= new MCTruth();
LoadTruth();
fTrack= new ShowerTrack();
fTrack->fEthrEl=0.;
fTrack->fDirVers[0]=fHeader["ELDIRCOS"];
fTrack->fDirVers[1]=fHeader["EMDIRCOS"];
fTrack->fDirVers[2]=fHeader["ENDIRCOS"];
fTrack->fInitPos[0]=fHeader["EX0COMES"]*m;
fTrack->fInitPos[1]=fHeader["EY0COMES"]*m;
fTrack->fInitPos[2]=fHeader["EZ0COMES"]*m;
if(fHeader["HITGROUN"]==0) fTrack->fHitGround=false ;
else fTrack->fHitGround=true ;
fTrack->fEnergy = fHeader["EPRMENRG"]*GeV;
fTrack->fTheta = fHeader["EPOLANGD"]*deg;
fTrack->fPhi = fHeader["EAZIANGD"]*deg;
fTrack->fX1 = fHeader["EPRMDPTH"]*g/cm2;
fTrack->fEarthImpact.SetXYZ(fHeader["EXDMPPOS"]*m,fHeader["EYDMPPOS"]*m,fHeader["EZDMPPOS"]*m);
while( 1 ) {
char* ans = gzgets((void*)fFile,s,5000);
if ( ans == Z_NULL ) {
cout << "Error reading gzgets in UnisimFileShowerGenerator::Load()n";
delete fTrack;
fTrack=NULL;
return false;
}
if (strncmp(s," BEVTBEND EEND",14) == 0) return true;
if (strncmp(s," STEPNUMB",9) == 0) {
char key;
int n;
if ( sscanf(s,"%s %d %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf",
&key,&n,&fStep.fXi,&fStep.fXf,&fStep.fXYZi[0],&fStep.fXYZi[1],&fStep.fXYZi[2],
&fStep.fXYZf[0],&fStep.fXYZf[1],&fStep.fXYZf[2],
&fStep.fTimei,&fStep.fTimef,&fStep.fAgei,&fStep.fAgef,&fStep.fNelectrons,
&fStep.fNcharged,&fStep.fEloss,&fStep.fNcherenkov) != 18 ) {
cout << "Data format error in UnisimFileShowerGenerator::LoadHeader() n";
cout << s << endl;
return false;
}
fStep.fXYZi[0]=fStep.fXYZi[0]*m;
fStep.fXYZi[1]=fStep.fXYZi[1]*m;
fStep.fXYZi[2]=fStep.fXYZi[2]*m;
fStep.fXYZf[0]=fStep.fXYZf[0]*m;
fStep.fXYZf[1]=fStep.fXYZf[1]*m;
fStep.fXYZf[2]=fStep.fXYZf[2]*m;
fStep.fTimei=fStep.fTimei*microsecond;
fStep.fTimef=fStep.fTimef*microsecond;
fStep.fXi=fStep.fXi*gram/cm2;
fStep.fXf=fStep.fXf*gram/cm2;
fTrack->Add(fStep);
}
}
return true;
}
//______________________________________________________________________________
bool UnisimFileShowerGenerator::LoadHeader() {
//
// Load header and store it into a map (key,value)
//
fHeader.clear();
char s[5010];
while(1) {
char* ans = gzgets((void*)fFile,s,5000);
if ( ans == Z_NULL ) {
cout << "Error reading gzgets in UnisimFileShowerGenerator::LoadHeader() (2) n";
return false;
}
if (strncmp(s," BEVTBUFF EVTB",14) != 0) {
char key[200];
double val=0.;
if ((strncmp(s," EVTHEADR EVTH",14) != 0) && (strncmp(s," RUNHEADR RUNH",14) != 0)) {
if ( sscanf(s,"%s %lf",key,&val) != 2 ) {
cout << "Data format error in UnisimFileShowerGenerator::LoadHeader() n";
cout << s << endl;
return false;
}
fHeader[key]=val;
}
} else {
return true;
}
}
return true;
}
//______________________________________________________________________________
void UnisimFileShowerGenerator::LoadTruth() {
//
//
//
fTruth->SetEnergy(fHeader["EPRMENRG"]*GeV); // in ESAF units (MeV)
fTruth->SetThetaPhi(fHeader["EPOLANGD"]*deg, fHeader["EAZIANGD"]*deg);
EVector v;
v[0] = fHeader["EX0COMES"]*m;
v[1] = fHeader["EY0COMES"]*m;
v[2] = fHeader["EZ0COMES"]*m;
fTruth->SetFirstInt(v,fHeader["EPRMDPTH"]);
v[0] = fHeader["EXDMPPOS"]*m;
v[1] = fHeader["EYDMPPOS"]*m;
v[2] = fHeader["EZDMPPOS"]*m;
// FIXME, Age on the Earth is not stored in files. Puting arbitrary number 2.
fTruth->SetEarthImpact(v,2);
// In case no value for ShowerMax is present Puting arbitrary number 0.
v[0] = 0;
v[1] = 0;
v[2] = 0;
map<string,double>::iterator p;
p = fHeader.find("EXMAXPOS") ;
if ( p != fHeader.end() ) v[0] = fHeader["EXMAXPOS"]*m;
p = fHeader.find("EYMAXPOS") ;
if ( p != fHeader.end() ) v[1] = fHeader["EYMAXPOS"]*m;
p = fHeader.find("EZMAXPOS") ;
if ( p != fHeader.end() ) v[2] = fHeader["EZMAXPOS"]*m;
fTruth->SetShowerMax(v,fHeader["ESHOWMAX"]*g/cm2);
fTruth-> SetHclouds(0.);
p = fHeader.find("ESTAHEIG") ;
if ( p != fHeader.end() ) fTruth-> SetHclouds(fHeader["ESTAHEIG"]*km); // in ESAF units
}