// ESAF : Euso Simulation and Analysis Framework
// $Id: RootInputModule.cc,v 1.52 2005/11/08 10:40:59 naumov Exp $
// Marco Pallavicini created Oct, 16 2003
#include "RootInputModule.hh"
#include "EEvent.hh"
#include "EHeader.hh"
#include "ETruth.hh"
#include "EDetector.hh"
#include "EDetPhoton.hh"
#include "EFee.hh"
#include "EAFee.hh"
#include "EMacroCell.hh"
#include "EMacroCellHit.hh"
#include "ERunParameters.hh"
#include "EAtmosphere.hh"
#include "ESinglePhoton.hh"
#include "EBunchPhotons.hh"
#include "RecoEvent.hh"
#include "RecoPhotonOnPupilData.hh"
#include "RecoCellInfo.hh"
#include "RecoPixelData.hh"
#include "RecoPhotoElectronData.hh"
#include "RecoAnalogData.hh"
#include "RecoRootEvent.hh"
#include <string>
#include <iomanip>
#include <TFile.h>
#include <TBranch.h>
#include <TTree.h>
#include <map>
#include <vector>
using namespace sou;
ClassImp(RootInputModule)
//_____________________________________________________________________________
RootInputModule::RootInputModule() : InputModule( string("Root") ) {
//
// ctor
//
fEventCounter = 0;
SetRecoEvent((RecoEvent*)0);
}
//_____________________________________________________________________________
RootInputModule::~RootInputModule() {
//
// dtor
//
if ( fRootFile )
fRootFile->Close();
}
//_____________________________________________________________________________
void RootInputModule::DestroyEvent() {
//
// Destroy event
//
if ( GetRecoEvent() )
delete GetRecoEvent();
SetRecoEvent((RecoEvent*)0);
}
//_____________________________________________________________________________
Bool_t RootInputModule::Init() {
//
// init: open file, tree, branches
//
// get file name and open root file
string fn = Conf()->GetStr("RootInputModule.FileName");
fInputFileName = fn;
if ( !( fRootFile = (TFile*)gROOT->GetListOfFiles()->FindObject( fn.c_str())))
fRootFile = new TFile(fn.c_str());
else
Msg(EsafMsg::Panic) << "TFile already existing in RootInputModule" << MsgDispatch;
// check file
if ( fRootFile->IsZombie() ) {
Msg(EsafMsg::Panic) << "Invalid root file in RootInputModule" << MsgDispatch;
}
// get tree
fNumEvents = 0;
fTree = (TTree*)fRootFile->Get("etree");
if ( fTree ) {
fNumEvents = (Int_t)fTree->GetEntries();
if (fNumEvents == 1)
Msg(EsafMsg::Info) << "Root file contains " << " 1 event" << MsgDispatch;
else
Msg(EsafMsg::Info) << "Root file contains " << fNumEvents << " events" << MsgDispatch;
}
else {
Msg(EsafMsg::Panic) << "Invalid root file " << fn << MsgDispatch;
}
fFirstEvent = (Int_t)Conf()->GetNum("RootInputModule.FirstEvent");
fLastEvent = (Int_t)Conf()->GetNum("RootInputModule.LastEvent");
Int_t maxEvents = fNumEvents;
if ( fFirstEvent < 0 ) {
Msg(EsafMsg::Warning) << "FirstEvent<0 in cfg file. Setting it to 0" << MsgDispatch;
fFirstEvent = 0;
} else if (fFirstEvent > maxEvents - 1) {
Msg(EsafMsg::Warning) << "FirstEvent> NumEvents in file. Setting it to the last event" << MsgDispatch;
fFirstEvent = fNumEvents - 1;
}
if ( Conf()->GetBool("RootInputModule.TruncateNumberOfEvents") ) {
maxEvents = (Int_t)Conf()->GetNum("RootInputModule.MaxEvents");
} else {
maxEvents = fNumEvents - fFirstEvent;
}
if ( maxEvents > (fNumEvents-fFirstEvent) ) maxEvents = fNumEvents - fFirstEvent;
if (fLastEvent < 0 || fLastEvent >= (maxEvents + fFirstEvent) ) fLastEvent = maxEvents + fFirstEvent - 1;
if ( fNumEvents == 1 || fLastEvent < fFirstEvent ) fLastEvent = fFirstEvent;
if ( fLastEvent < fFirstEvent )
Msg(EsafMsg::Warning) << "LastEvent < FirstEvent. Processing only FirstEvent." << MsgDispatch;
fEv = new EEvent();
// let the event branch on the tree
fEv->SetBranches( fTree );
// // get run parameters tree
// fRunTree = (TTree*)fRootFile->Get("runtree");
// if ( !fRunTree ) {
// throw runtime_error("Invalid root file "+fn);
// }
//
// fRunPars = new ERunParameters();
// fRunTree->SetBranchAddress("runpars", &fRunPars);
// fRunBranch = fRunTree->GetBranch("runpars");
// fRunTree->GetEntry(0);
fRunPars = fEv->GetRunPars();
if(fRunPars) {
Msg(EsafMsg::Debug) << "Detector Configuration:" << MsgDispatch;
Msg(EsafMsg::Debug) << "Number of pixels = " << fRunPars->GetPixelMap().GetNumPixels() << MsgDispatch;
Msg(EsafMsg::Debug) << "Number of pmts = " << fRunPars->GetNumPmts() << MsgDispatch;
Msg(EsafMsg::Debug) << "Number of pixels per pmt = " << fRunPars->GetPmtData().GetNumPads() << MsgDispatch;
}
// start from the beginning
fEventCounter = fFirstEvent;
return kTRUE;
}
//_____________________________________________________________________________
Bool_t RootInputModule::Done() {
//
// End of run
//
return kTRUE;
}
//_____________________________________________________________________________
RecoEvent *RootInputModule::GetEvent() {
//
// Return on event reading from root file
//
// get data
fTree->GetEntry(fEventCounter);
if ( fEventCounter++ > fLastEvent ) {
return (RecoEvent*)0;
}
Msg(EsafMsg::Info) << "=========================> Processing Event " << setw(4) <<
fEventCounter - 1 << " <===========================" << MsgDispatch;
EHeader *header = fEv->GetHeader();
ETruth *truth = fEv->GetTruth();
EDetector *detector = fEv->GetDetector();
EAtmosphere *atmosphere = fEv->GetAtmosphere();
//FIXME
//if ( detector == 0 ) throw an error
// destroy previous event and create a new one
DestroyEvent();
RecoEvent *ev = new RecoEvent();
// fill reco event
ev->fHeader.fNum = header->GetNum();
ev->fHeader.fRun = header->GetRun();
// fill photons on pupil data
ev->fHeader.fTruePhotonsOnPupil = 0;
ev->fHeader.fTrueFluoOnPupil = 0;
ev->fHeader.fTrueBackCerOnPupil = 0;
ev->fHeader.fGtuLength = 0;
if ( atmosphere ) {
Int_t ns = atmosphere->GetNumSingles();
Int_t nbunch = atmosphere->GetNumBunches();
Msg(EsafMsg::Debug)<<" Number of single photons "<< ns <<MsgDispatch;
Msg(EsafMsg::Debug)<<" Number of bunches "<< nbunch <<MsgDispatch;
for (Int_t i=0; i<ns; i++) {
ESinglePhoton *sph = atmosphere->GetSingle(i);
if (! ( sph->IsAbsorbed()) ) {
ev->fHeader.fTruePhotonsOnPupil++;
if ( sph->GetType() == 0 ) ev->fHeader.fTrueFluoOnPupil++;
if ( sph->GetType() == 1 && sph->GetHistory()>=2 ) ev->fHeader.fTrueBackCerOnPupil++;
ev->PutRecoPhotonOnPupilData( new RecoPhotonOnPupilData( sph->GetType(), sph->GetHistory(),
sph->GetWl(),sph->GetTof()+sph->GetDate() ));
}
}
// get 3D position of fluo profile maximum
Double_t Weightmax = 0.;
Int_t indexofmax = 0;
EBunchPhotons *bunch = 0;
for(Int_t i=0; i<nbunch; i++) {
bunch = atmosphere->GetBunch(i);
if(bunch->GetType() == 0) {
if(bunch->GetWeight() > Weightmax) {
Weightmax = bunch->GetWeight();
indexofmax = i;
}
}
}
if ( nbunch>0 ) {
bunch = atmosphere->GetBunch(indexofmax);
ev->fHeader.fTrueFluoMaxPos = 0.5*(bunch->GetShowerPosi() + bunch->GetShowerPosf());
}
}
else Msg(EsafMsg::Info) <<"Atmosphere object is NULL : no photons on pupil data"<<MsgDispatch;
Msg(EsafMsg::Debug)<<" Number of photons on pupil "<< ev->fHeader.fTruePhotonsOnPupil<<MsgDispatch;
Msg(EsafMsg::Debug)<<" Number of fluorescence photons on pupil "<< ev->fHeader.fTrueFluoOnPupil<<MsgDispatch;
Msg(EsafMsg::Debug)<<" Number of backscattered cerenkov photons on pupil "<< ev->fHeader.fTrueBackCerOnPupil<<MsgDispatch;
// GTU length interval
if(detector) ev->fHeader.fGtuLength = detector->GetElectronics()->GetGtuLength();
// fill photoelectrons info
ev->fHeader.fTruePhotoElectrons = 0;
ev->fHeader.fTrueNumFastOR = 0;
ev->fHeader.fRunPars = fRunPars;
ev->fHeader.fDetector = detector;
Int_t fNumPhotons = 0;
if(detector) {
fNumPhotons = detector->GetNumPhotons();
for(Int_t i=0; i < fNumPhotons; i++) {
EDetPhoton *ph = detector->GetPhoton(i);
if ( ph->MadeCount() ) {
ev->fHeader.fTruePhotoElectrons++;
ev->PutRecoPhotoElectronData(
new RecoPhotoElectronData(ph->GetFe(), ph->GetPixelUID(),
ph->GetGtu(), ph->GetMacroCell(), (Int_t)ph->MadeFastOR(),
ph->GetXCell(), ph->GetYCell(), ph->GetPos()[0],
ph->GetPos()[1]) );
if ( ph->MadeFastOR() )
ev->fHeader.fTrueNumFastOR++;
}
}
Msg(EsafMsg::Debug) << "Number of Photons = " << fNumPhotons << MsgDispatch;
if ( fNumPhotons ) {
Msg(EsafMsg::Debug) << "Number of detected photoelectrons = " << ev->fHeader.fTruePhotoElectrons << MsgDispatch;
Msg(EsafMsg::Debug) << "Number of photelectrons counted in macrocells = " << ev->fHeader.fTrueNumFastOR << MsgDispatch;
}
// fill macrocell data
ev->fHeader.fNumCells = detector->GetNumCells();
ev->fHeader.fNumCellHits = detector->GetNumCellHits();
Msg(EsafMsg::Debug) << "Number of macrocells = " << detector->GetNumCells() << MsgDispatch;
Msg(EsafMsg::Debug) << "Number of macrocell hits = " << detector->GetNumCellHits() << MsgDispatch;
// loop on macrocells
for(Int_t i=0; i<detector->GetNumCells(); i++) {
EMacroCell* mc = detector->GetMacroCell(i);
Int_t cell_id = mc->GetMCId();
RecoCellInfo *info = new RecoCellInfo( cell_id );
info->fTriggered = mc->HasTriggered();
info->fTriggerWord = mc->GetTriggerWord();
ev->PutRecoCellInfo(cell_id, info);
}
// fill macrocell hits (if any)
for(Int_t i=0; i<detector->GetNumCellHits(); i++) {
EMacroCellHit* mch = detector->GetMacroCellHit(i);
Int_t cell_id = mch->GetCellId();
RecoCellInfo *info = ev->GetRecoCellInfo(cell_id);
if ( info == NULL) {
Msg(EsafMsg::Warning) << "Found a EMacroCellHit not corresponding to any EMacroCell" << MsgDispatch;
}else {
info->AddHit( mch->GetChUId(),mch->GetCounts(),mch->GetGtu() );
ev->PutRecoCellInfo(cell_id, info);
}
}
// fill pixel front end data
ev->fHeader.fNumFee = detector->GetNumFee();
Msg(EsafMsg::Debug) << "Number of front end electronics = " << detector->GetNumFee() << MsgDispatch;
// loop on efee hits
ev->fHeader.fNumActiveFee = 0;
for (Int_t i=0; i<detector->GetNumFee(); i++) {
EFee* fee = detector->GetFee(i);
Int_t hits = fee->GetNumHits();
if ( hits != 0) {
Int_t pid = fee->GetChUId();
if ( fRunPars->GetPixelMap().GetThetaFOV(pid) != -1. ) {
ev->fHeader.fNumActiveFee++;
RecoPixelData *pd = new RecoPixelData( pid, fee->GetNumHits(), fee->GetGtu() );
const EAnglePixelMap& map = fRunPars->GetPixelMap();
pd->fTheta = map.GetThetaFOV(pid);
pd->fThetaSigma = map.GetSigmaThetaFOV(pid);
pd->fPhi = map.GetPhiFOV(pid);
pd->fPhiSigma = map.GetSigmaPhiFOV(pid);
Int_t pmtid = 1 + (pid-1)/(fRunPars->GetPmtData().GetNumPads());
pd->fMacroCellId = fRunPars->GetPmtGeo(pmtid)->GetMC();
pd->SetCountsSignal(fee->GetNumSignals());
ev->PutRecoPixelData( pd );
}
}
}
// add MC information on the photons type giving the signal
if(detector) {
// save map of <PixelUID, vector <Photons entry in detector->GetPhoton()>>
map <Int_t, vector<Int_t> > photons;
fNumPhotons = detector->GetNumPhotons();
for(Int_t j=0; j < fNumPhotons; j++) {
EDetPhoton *ph = detector->GetPhoton(j);
if ( ph->MadeCount() ) photons[ph->GetPixelUID()].push_back(j);
}
vector<RecoPixelData*> fPixels = (vector<RecoPixelData*>)ev->GetRecoPixels();
vector<RecoPixelData*>::iterator itPixels;
for (itPixels = fPixels.begin(); itPixels != fPixels.end(); itPixels++) {
Int_t uid = (*itPixels)->GetPixelId();
vector<Int_t> pe = photons[uid];
for (UInt_t i=0; i<pe.size(); i++) {
EDetPhoton *ph = detector->GetPhoton(pe[i]);
if (ph->GetType() == kFluorescence && ph->GetGtu() == (*itPixels)->GetGtu()) (*itPixels)->AddCountFluo();
if (ph->GetType() == kDirectCherenkov && ph->GetGtu() == (*itPixels)->GetGtu()) (*itPixels)->AddCountCherDirect();
if (ph->GetType() == kBackCherenkov && ph->GetGtu() == (*itPixels)->GetGtu()) (*itPixels)->AddCountCherScatt();
}
// calculate number of noise p.e. as the difference of observed and accounted to the signal
Int_t noise = (*itPixels)->GetCounts() -
(*itPixels)->GetCountsFluo() -
(*itPixels)->GetCountsCherDirect() -
(*itPixels)->GetCountsCherScatt();
(*itPixels)->SetCountsNoise(noise);
}
}
Msg(EsafMsg::Debug) << "Number of active front end electronics = " << ev->fHeader.fNumActiveFee << MsgDispatch;
// fill analog front end data
ev->fHeader.fNumAFee = detector->GetNumAFee();
Msg(EsafMsg::Debug) << "Number of analog front end electronics = " << detector->GetNumAFee()<< MsgDispatch;
// loop on eafee hits
for (Int_t i=0; i<detector->GetNumAFee(); i++) {
EAFee* afee = detector->GetAFee(i);
ev->PutRecoAnalogData( new RecoAnalogData( afee->GetGtu(), afee->GetFEId() ) );
}
}
else {
Msg(EsafMsg::Info) << "No Detector found in Datan"
<<"GTUlength drawn from file Profile.cfg" << MsgDispatch;
}
if(ev->fHeader.fGtuLength == 0) {
Float_t gtulength = Conf()->GetNum("RootInputModule.GTUlength");
ev->fHeader.fGtuLength = gtulength;
}
// add MC truth to the header and
ev->fHeader.fTrueTheta = truth->GetTrueTheta();
ev->fHeader.fTruePhi = truth->GetTruePhi();
ev->fHeader.fTrueEnergy = truth->GetTrueEnergy();
ev->fHeader.fTrueInitPos = truth->GetTrueInitPos();
ev->fHeader.fTrueX1 = truth->GetTrueX1()*cm2/g;
ev->fHeader.fTrueEarthImpact = truth->GetTrueEarthImpact();
ev->fHeader.fTrueEarthAge = truth->GetTrueEarthAge();
ev->fHeader.fTrueShowerMaxPos = truth->GetTrueShowerMaxPos();
ev->fHeader.fTrueShowerXMax = truth->GetTrueShowerXMax()*cm2/g;
ev->fHeader.fTrueHclouds = truth->GetTrueHclouds()/km;
// save pointer to reco event and return it
SetRecoEvent( ev );
return GetRecoEvent();
}
//_____________________________________________________________________________
Bool_t RootInputModule::SaveRootData(RecoRootEvent* fRecoRootEvent) {
ETruth *truth = fEv->GetTruth();
// save truth data into the root file
fRecoRootEvent->GetTruth().SetEnergy(truth->GetTrueEnergy());
fRecoRootEvent->GetTruth().SetTheta(truth->GetTrueTheta());
fRecoRootEvent->GetTruth().SetPhi(truth->GetTruePhi());
fRecoRootEvent->GetTruth().SetX1(truth->GetTrueX1()*cm2/g);
fRecoRootEvent->GetTruth().SetInitPos(truth->GetTrueInitPos());
fRecoRootEvent->GetTruth().SetEarthImpact(truth->GetTrueEarthImpact());
fRecoRootEvent->GetTruth().SetEarthAge(truth->GetTrueEarthAge());
fRecoRootEvent->GetTruth().SetXMax(truth->GetTrueShowerXMax()*cm2/g);
fRecoRootEvent->GetTruth().SetMaxPos(truth->GetTrueShowerMaxPos());
fRecoRootEvent->GetTruth().SetHclouds(truth->GetTrueHclouds()/km);
// Fill Header info about the event
fRecoRootEvent->GetRecoHeader().Set(GetRecoEvent()->GetHeader().GetRun(),
GetRecoEvent()->GetHeader().GetNum());
return kTRUE;
}