// ESAF : Euso Simulation and Analysis Framework
// $Id: UnisimPhotonsOnPupil.cc,v 1.12 2005/05/08 22:28:22 thea Exp $
// M. Pallavicini created Jul, 2 2002
#include "UnisimPhotonsOnPupil.hh"
#include "utils.hh"
#include "EVector.hh"
#include "Config.hh"
#include "SystemOfUnits.hh"
#include <zlib.h>
ClassImp(UnisimPhotonsOnPupil)
//_________________________________________________________________________________________
UnisimPhotonsOnPupil::UnisimPhotonsOnPupil() : PhotonsOnPupil(), fTheTruth(0){
//
// Constructor
//
fPhotons = 0;
// local buffer for photons
fBufferSize = 200000;
NewPhotonBuffer(fBufferSize);
/* Unisim photons want to know where they must be put.
Now such parameters are set in FileUnisimLightToEuso.
In the future they should be retrieved from a detector
configuration set of informations
*/
ConfigFileParser *pippo=Config::Get()->GetCF("LightToEuso","FileUnisimLightToEuso") ;
fZdiskEuso = pippo->GetNum("FileUnisimLightToEuso.zdiskEuso") ;
fRmaxDiskEuso = pippo->GetNum("FileUnisimLightToEuso.rmaxdiskEuso") ;
fEvCounter = 0;
}
//_________________________________________________________________________________________
UnisimPhotonsOnPupil::~UnisimPhotonsOnPupil() {
//
// Destructor
//
DeletePhotonBuffer();
if ( fTheTruth )
delete fTheTruth;
}
//_________________________________________________________________________________________
void UnisimPhotonsOnPupil::NewPhotonBuffer(Int_t s) {
//
// Create a larger buffer and copy what is already in the old one
//
ParentPhoton** newbuff = (ParentPhoton**) new char[4*s];
for(Int_t i=0; i<s; i++)
newbuff[i]=0;
if ( fPhotons ) {
for(Int_t i=0; i<s; i++) {
if (i<fBufferSize)
newbuff[i] = fPhotons[i];
}
delete[] fPhotons;
}
fPhotons = newbuff;
}
//_________________________________________________________________________________________
void UnisimPhotonsOnPupil::DeletePhotonBuffer() {
//
// Clean the buffer and release memory
//
if ( fPhotons ) {
for(Int_t i=0; i<fBufferSize; i++) {
if (fPhotons[i])
delete fPhotons[i];
fPhotons[i]=0;
}
delete[] fPhotons;
}
}
//_________________________________________________________________________________________
Bool_t UnisimPhotonsOnPupil::SkipEvent( FILE* fp ) {
//
// Skip one event
//
Clear();
char s[5010];
fEvCounter++;
Bool_t firstLine = kTRUE;
while( kTRUE ) {
char* ans = gzgets((void*)fp,s,5000);
if ( ans == Z_NULL ) {
MsgForm(EsafMsg::Warning,"Error reading gzgets in UnisimPhotonsOnPupil::Load()");
return kFALSE;
}
if (!firstLine && strncmp(s," RUNHEADR RUNH",14) == 0)
return kTRUE;
if ( firstLine ) firstLine = kFALSE;
}
}
//_________________________________________________________________________________________
Bool_t UnisimPhotonsOnPupil::Load( FILE* fp ) {
//
// Load one event from file
//
Clear();
fEvCounter++;
Int_t nstart=0;
if ( !LoadHeader( fp ) )
return false;
char s[5010];
while( 1 ) {
char* ans = gzgets((void*)fp,s,5000);
if ( ans == Z_NULL ) {
MsgForm(EsafMsg::Warning,"Error reading gzgets in UnisimPhotonsOnPupil::Load()");
return false;
}
if (strncmp(s," RUNHEADR RUNH",14) == 0)
return true;
Int_t n;
double x[3],tm,w,th,ph,r,ph2;
Int_t t;
if ( sscanf(s,"%d %lf %lf %lf %lf %lf %lf %lf %lf %lf %d",
&n,&x[0],&x[1],&x[2],&tm,&w,&th,&ph,&r,&ph2,&t) != 11 ) {
MsgForm(EsafMsg::Warning,"Error reading file in UnisimPhotonsOnPupil");
return false;
}
nstart++;
if ( n != nstart )
MsgForm(EsafMsg::Warning,"WARNING: Photon counter not consequtive");
x[0]*=m;
x[1]*=m;
x[2]*=m;
tm *= microsecond;
// angles are all degrees in the file
th *= deg;
ph *= deg;
ph2 *= deg;
// Wavelength is forced to 357 nm later
w *= nm;
// compute position on the pupil
/* unisim photon positions (r, ph2) are given on a disk normal to photon direction
in local frame. Now we have to put them on a disk normal to nadir direction
*/
TVector3 phdir ;
phdir.SetX( sin(th)*sin(ph) ) ;
phdir.SetY( sin(th)*cos(ph) ) ;
phdir.SetZ( cos(th) ) ;
phdir = phdir.Unit() ;
TVector3 rdisk ;
r *= m;
// rdisk in local frame
rdisk.SetX( r*cos(ph2) ) ;
rdisk.SetY( r*sin(ph2) ) ;
rdisk.SetZ( 0 ) ;
// now transform rdsik to global frame
rdisk.RotateUz(phdir) ;
// calculate the impact point of the photon on a disk perpendicular to nadir
Double_t tprop ;
tprop = - rdisk.Z()/phdir.Z() ;
TVector3 ydiskunisim ;
ydiskunisim = rdisk + tprop*phdir ;
// put the photons on the EUSO-ESAF entrance disk at Z coordinate = zdiskEuso
Double_t rdiskactual ;
rdiskactual = sqrt( ydiskunisim.X()*ydiskunisim.X() + ydiskunisim.Y()*ydiskunisim.Y() ) ;
if ( fRmaxDiskEuso >= rdiskactual ) {
double y[3];
y[0] = ydiskunisim.X() ;
y[1] = ydiskunisim.Y() ;
y[2] = fZdiskEuso ;
if (!fPhotons[fNphotons])
fPhotons[fNphotons] = new ParentPhoton();
fPhotons[fNphotons++]->Build(n,x,y,th,ph,w,tm,(ParentPhotonType)t);
}
}
return true;
}
//_________________________________________________________________________________________
Photon* UnisimPhotonsOnPupil::Get() {
//
// returns a Photon object
//
//if ( (( fCurrent % 200 ) == 0) && (fCurrent) )
// cout << "Event: " << event << " Photon: " << fCurrent << "n";
if ( fCurrent == fNphotons )
return (Photon*)0;
return &(fPhotons[fCurrent++]->GetPhoton());
}
//_________________________________________________________________________________________
void UnisimPhotonsOnPupil::Clear() {
//
// Reset
//
fNphotons = 0;
fCurrent = 0;
fHeader.clear();
fTheTruth = 0;
}
//_________________________________________________________________________________________
Bool_t UnisimPhotonsOnPupil::LoadHeader( FILE* fp) {
//
// Load header and store it into a map (key,value)
//
char s[5010];
while(1) {
char* ans = gzgets((void*)fp,s,5000);
if ( ans == Z_NULL ) {
MsgForm(EsafMsg::Warning, "Error reading gzgets in UnisimPhotonsOnPupil::LoadHeader() (2)");
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 ) {
MsgForm(EsafMsg::Warning,"Data format error in UnisimPhotonsOnPupil::LoadHeader()");
Msg(EsafMsg::Warning) << s << MsgDispatch;
return false;
}
fHeader[key]=val;
}
} else {
for(Int_t i=0;i<11;i++) // skip 11 lines
gzgets((void*)fp,s,5000);
break;
}
}
return true;
}
//_________________________________________________________________________________________
MCTruth* UnisimPhotonsOnPupil::Truth() {
if (!fTheTruth) {
fTheTruth = new MCTruth();
fTheTruth->SetEnergy(fHeader["EPRMENRG"]*GeV); // now in MeV
fTheTruth->SetThetaPhi(fHeader["EPOLANGD"]*deg, fHeader["EAZIANGD"]*deg);
EVector v;
v[0] = fHeader["EX0COMES"]*m;
v[1] = fHeader["EY0COMES"]*m;
v[2] = fHeader["EZ0COMES"]*m;
fTheTruth->SetFirstInt(v,fHeader["EPRMDPTH"]*g/cm2);
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.
fTheTruth->SetEarthImpact(v,2);
// In case no value for ShowerMax is present Puting arbitrary number -100.
v[0] = -100;
v[1] = -100;
v[2] = -100;
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;
fTheTruth->SetShowerMax(v,fHeader["ESHOWMAX"]*g/cm2);
;
fTheTruth->SetHclouds(0.);
p = fHeader.find("ESTAHEIG") ;
if ( p != fHeader.end() ) fTheTruth-> SetHclouds(fHeader["ESTAHEIG"]*km); // in ESAF units
}
return fTheTruth;
}