// ESAF : Euso Simulation and Analysis Framework
// $Id: RootInputModule.cc,v 1.36 2005/03/30 09:19:52 naumov Exp $
// Marco Pallavicini created Oct, 16 2003
#include "RootInputModule.hh"
#include "EEvent.hh"
#include "EDetector.hh"
#include "EDetPhoton.hh"
#include "EFee.hh"
#include "EAFee.hh"
#include "EMacroCellHit.hh"
#include "ERunParameters.hh"
#include "RecoEvent.hh"
#include "RecoPhotoElectronData.hh"
#include "RecoCellInfo.hh"
#include "RecoPixelData.hh"
#include "RecoPhotoElectronData.hh"
#include "RecoAnalogData.hh"
#include "RecoRootEvent.hh"
#include <string>
#include <TFile.h>
#include <TBranch.h>
#include <TTree.h>
ClassImp(RootInputModule)
//_____________________________________________________________________________
// ctor
RootInputModule::RootInputModule() : InputModule( string("Root") ) {
fEventCounter = 0;
SetRecoEvent((RecoEvent*)0);
}
//_____________________________________________________________________________
// dtor
RootInputModule::~RootInputModule() {
if ( fRootFile )
fRootFile->Close();
}
//_____________________________________________________________________________
// destroy event
void RootInputModule::DestroyEvent() {
if ( GetRecoEvent() )
delete GetRecoEvent();
SetRecoEvent((RecoEvent*)0);
}
//_____________________________________________________________________________
// init: open file, tree, branches
Bool_t RootInputModule::Init() {
// get file name and open root file
string fn = Conf()->GetStr("RootInputModule.FileName");
if ( !( fRootFile = (TFile*)gROOT->GetListOfFiles()->FindObject( fn.c_str())))
fRootFile = new TFile(fn.c_str());
else
throw runtime_error("TFile already existing in RootInputModule");
// check file
if ( fRootFile->IsZombie() ) {
throw runtime_error("Invalid root file in RootInputModule");
}
// 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 {
throw runtime_error("Invalid root file "+fn);
}
fFirstEvent = (Int_t)Conf()->GetNum("RootInputModule.FirstEvent");
fLastEvent = (Int_t)Conf()->GetNum("RootInputModule.LastEvent");
if (fLastEvent < 0) {
if ( Conf()->GetStr("RootInputModule.TruncateNumberOfEvents")=="yes" ) {
Int_t maxEvents = (Int_t)Conf()->GetNum("RootInputModule.MaxEvents");
if ( fNumEvents > maxEvents ) {
Msg(EsafMsg::Warning) << "Reco will process only the first " << maxEvents << MsgDispatch;
fLastEvent = maxEvents - 1;
} else {
fLastEvent = fNumEvents -1;
}
}
else if ( Conf()->GetStr("RootInputModule.TruncateNumberOfEvents")=="no" ) {
fLastEvent = fNumEvents -1;
}
else {
Msg(EsafMsg::Panic) << "Wrong value in RootInputModule.TruncateNumberOfEvents." << MsgDispatch;
throw runtime_error("RootInputModule: wrong value in .cfg file.");
}
}
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);
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;
}
//_____________________________________________________________________________
// end of run
Bool_t RootInputModule::Done() {
return kTRUE;
}
//_____________________________________________________________________________
// return on event reading from root file
RecoEvent *RootInputModule::GetEvent() {
// get data
fTree->GetEntry(fEventCounter);
if ( fEventCounter++ > fLastEvent ) {
return (RecoEvent*)0;
}
Msg(EsafMsg::Info) << "Reading event " << fEventCounter - 1 << MsgDispatch;
EHeader *header = fEv->GetHeader();
ETruth *truth = fEv->GetTruth();
EDetector *detector = fEv->GetDetector();
//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();
// GTU length interval
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 = 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.fNumCellHits = detector->GetNumCellHits();
Msg(EsafMsg::Debug) << "Number of macrocell hits = " << detector->GetNumCellHits() << MsgDispatch;
// loop on macrocell hits
for(Int_t i=0; i<detector->GetNumCellHits(); i++) {
EMacroCellHit* mch = detector->GetMacroCellHit(i);
Int_t cell = mch->GetCellId();
RecoCellInfo *info = ev->GetRecoCellInfo( cell );
if ( info == NULL)
info = new RecoCellInfo(cell);
info->AddHit( mch->GetChUId(),mch->GetCounts(),mch->GetGtu() );
ev->PutRecoCellInfo(cell, 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
Bool_t gSign = (Bool_t)Conf()->GetNum("RootInputModule.ProcessOnlySignal");
ev->fHeader.fNumActiveFee = 0;
for (Int_t i=0; i<detector->GetNumFee(); i++) {
EFee* fee = detector->GetFee(i);
Int_t hits;
if (!gSign) hits = fee->GetNumHits();
else hits = fee->GetNumSignals();
if ( hits != 0) {
Int_t pid = fee->GetChUId();
if ( fRunPars->GetPixelMap().GetThetaFOV(pid) != -1. ) {
ev->fHeader.fNumActiveFee++;
RecoPixelData *pd = new RecoPixelData( pid, hits, 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();
ev->PutRecoPixelData( pd );
}
}
}
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() ) );
}
// 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;
// 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());
// Fill Header info about the event
fRecoRootEvent->GetRecoHeader().Set(GetRecoEvent()->GetHeader().GetRun(),
GetRecoEvent()->GetHeader().GetNum());
return kTRUE;
}