Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

RootInputModule - source file

// 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;
}
About Us | EUSO Official Website | Web pages created by Roberto Pesce and Alessandro Thea - Last Update 14-May-2005 21:31