Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

RootInputModule - source file

// 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;
}
About Us | EUSO Official Website | Web pages created by Roberto Pesce and Alessandro Thea - Last Update Wed Nov 16 16:57:39 2005 Wed Nov 16 16:29:22 2005