---

readESDHist.cpp

Go to the documentation of this file.

#include <XANADOO/domain/interface/XANAEsdEvent.h>
#include <XANADOO/XANAClusters/interface/XANACluster.h>
#include <XANADOO/XANAClusters/interface/XANACaloRecHit.h>
#include <XANADOO/XANAClusters/interface/XANAEmCluster.h>
#include <XANADOO/XANAClusters/interface/XANASuperCluster.h>
#include <XANADOO/XANAClusters/interface/XANAHadCluster.h>
#include <XANADOO/XANATracks/interface/XANATrack.h>
#include <XANADOO/XANAMcInfo/interface/XANAGeneratorEvent.h>
#include <XANADOO/XANAMcInfo/interface/XANAGeantEvent.h>
#include <XANADOO/XANAMcInfo/interface/XANAGeneratorParticle.h>
#include <XANADOO/XANAElectronCandidate/interface/XANAElectronCandidate.h>

#include "TTree.h"
#include "TFile.h"
#include "TBranch.h"
#include "TClonesArray.h"
#include "TObjArray.h"
#include "TH1.h"
#include "TH2.h"
#include "TSystem.h"

#include <string>

// Comment this line if not single electron
#define SINGLE_ELECTRON


using namespace std;

int main(int argc, char** argv)
{
  string inputFileName = "h300eemm.root";
  string histoFileName = "hist_"+inputFileName;
 
  //   Connect to file generated
  TFile f(inputFileName.c_str());

 //   Booking the histograms
  
  // Clusters
  TH1F hClustersNumber("hClustersNumber", "Number of clusters", 20, 0., 20.);
  TH1F hClustersEnergy("hClustersEnergy", "Energy of cluster", 100, 0., 500.);
  TH1F hSuperClustersNumber("hSuperClustersNumber", "Number of super clusters", 50, 0., 50.);
  TH1F hSuperClustersEnergy("hSuperClustersEnergy", "Energy of super cluster", 100, 0., 500.);
  TH1F hBremsRecovered("hBremsRecovered","Energy of recovered bremsstrahlung photons",100,0.,200.);
  TH2F hBremsDetaDphi("hBremsDetaDphi","Delta phi vs Delta eta for brems clusters",100,-0.3,0.3,100,-0.3,0.3);
  
  TH1F hHadClustersNumber("hHadClustersNumber", "Number of hadron clusters", 20, 0., 20.);
  TH1F hHadClustersEnergy("hHadClustersEnergy", "Energy of hadron cluster", 100, 0., 500.);
  
  TH1F hPreshClustersNumber("hPreshClustersNumber", "Number of preshower clusters", 20, 0., 20.);

  // Hits
  TH1F hHitsNumber("hHitsNumber", "Number of hits", 100, 0., 100.);
  TH1F hHitsEnergy("hHitsEnergy", "Energy of hit", 200, 0., 200.);
  
  // Tracks
  TH1F hTracksNumber("hTracksNumber", "Number of tracks", 10, 0., 10.);
  TH1F hTracksMomentumAtVertex("hTracksMomentumAtVertex", "Momentum of track at vertex", 100, 0., 500.);
  TH1F hTracksMomentumAtFirstPoint("hTracksMomentumAtFirstPoint", "Momentum of track at first point", 100, 0., 500.);
  TH1F hTracksMomentumAtLastPoint("hTracksMomentumAtLastPoint", "Momentum of track at last point", 100, 0., 500.);
 
  //Vertex 
  TH1F hVertexNumber("hVertexNumber", "Number of vertices", 10, 0., 10.);
  TH1F hVertexTracksNumber("hVertexTracksNumber", "Number of vertex tracks", 10, 0., 10.);
  
  // Monte Carlo information 
  
  // generator event    
  TH1F hGenNumberOfParticles("hGenNumberOfParticles","Generator info: number of particles", 100, 0., 100.);
  
#ifdef SINGLE_ELECTRON
  TH1F hGenElectronMomentum("hGenElectronMomentum","Generator info: Electron momentum", 100, 0., 500.);
  TH1F hGenElectronPt("hGenElectronPt","MGeneratorC info: Electron pT", 200, 0., 200.);
  TH1F hGenElectronEta("hGenElectronEta","Generator info: Electron eta", 100, -3.,3.);
  TH1F hGenElectronPhi("hGenElectronPhi","Generator info: Electron phi", 100, -3.2,3.2);
  XANAGeneratorParticle *tmpGenElectron = 0;
#endif
  
  // GEANT event
  TH1F hGeantNumberOfTracks("hGeantNumberOfTracks","GEANT info: Number of tracks", 50, 0., 50.);
  TH1F hGeantNumberOfVertices("hGeantNumberOfVertices","GEANT info: Number of vertices", 50, 0., 50.);

  // electron candidates
  TH1F hElectronNumber("hElectronNumber", "Number of electrons", 10, 0., 10.);    
  TH1F hElectronEscOverP("hElectronEscOverP", "Electron: E supercluster over P", 100, 0.5, 1.5);
  TH1F hElectronDetaSc("hElectronDetaSc", "Electron: Delta eta super cluster", 100, -0.2, 0.2);   
  TH1F hElectronHoE("hElectronHoE", "Electron: Hadronic over Em", 100, 0.5, 1.5);
  
  //   Read Tree named "T" in memory. Tree pointer is assigned the same name
  TTree *T = (TTree*)f.Get("ESD");
   
// Reconstructed event
  XANAEsdEvent *event = 0;  
  T->SetBranchAddress("Event",&event);
  
  Short_t nc = 0, nh = 0, nhc = 0, ns = 0, npc = 0,  nt = 0, nv = 0, nvt = 0;  
  Short_t nbrems = 0, ne = 0;
  
  TObjArray *clusters = 0;
  TObjArray *hadClusters = 0;
  TObjArray *superClusters = 0;
  TObjArray *preshClusters = 0;
  TRefArray *hits = 0;
  TObjArray *tracks = 0;
  TObjArray *vertices = 0;
  TObjArray *electrons = 0;

  XANAEmCluster *tmpEmCluster = 0;
  XANASuperCluster *tmpSuperCluster = 0;
  XANAEmCluster *tmpSeedCluster = 0;
  XANAEmCluster *tmpBremCluster = 0; 
  XANAHadCluster *tmpHadCluster = 0;
  XANAElectronCandidate *tmpElectronCandidate = 0;

  // MC information (generator and GEANT event)
  XANAGeneratorEvent *generatorEvent = new XANAGeneratorEvent();
  TBranch *branchGeneratorEvent = T->GetBranch("GeneratorEvent");
  branchGeneratorEvent->SetAddress(&generatorEvent);
    
  XANAGeantEvent *geantEvent = new XANAGeantEvent();
  TBranch *branchGeantEvent = T->GetBranch("GeantEvent");
  branchGeantEvent->SetAddress(&geantEvent);

  Short_t mc_np = 0, mc_nsp = 0, mc_nssp = 0, mc_npp = 0, mc_nspp = 0, mc_ngt = 0, mc_ngv = 0;

  //   Start main loop on all events
  Int_t nb = 0; 
  Int_t nevent = (Int_t)T->GetEntries();
  for (Int_t i=0; i<nevent; i++) {
    nb += T->GetEntry(i); 

    cout << "================ Event: " << i << " ================" << endl;    
    
    // Reconstructed event
    
    cout << "Reconstructed:" ;
    
    nc  = event->getNumberOfEmClusters();      clusters = event->getEmClusters(); 
    nh  = event->getNumberOfEmClusterHits();          hits = event->getEmClusterHits();   
    nhc = event->getNumberOfHadClusters();   hadClusters = event->getHadClusters();   
    ns  = event->getNumberOfSuperClusters(); superClusters = event->getSuperClusters(); 
    npc = event->getNumberOfPreshClusters(); preshClusters = event->getPreshClusters(); 
    nt  = event->getNumberOfTracks();        tracks = event->getTracks();
    nv  = event->getNumberOfVertices();      vertices = event->getVertices();   
    nvt = event->getNumberOfVertexTracks();  
    ne  = event->getNumberOfElectronCandidates(); electrons = event->getElectronCandidates();
 
    cout << "  clust, hits, hadClu, supers, presho, tracks, vertex, vTrack, nElec" << endl;
    printf("  \t \t %d\t%d\t%d\t%d\t%d\t%d\t%d\t%d\t%d\n", nc, nh, nhc, ns, npc, nt, nv, nvt, ne);

    // Monte Carlo information
   
    // generator event
    mc_np   = generatorEvent->getNumberOfParticles();
    mc_nsp  = generatorEvent->getNumberOfSignalParticles();
    mc_nssp = generatorEvent->getNumberOfStableSignalParticles();
    mc_npp  = generatorEvent->getNumberOfPileupParticles();
    mc_nspp = generatorEvent->getNumberOfStablePileupParticles();
    
    // GEANT event
    mc_ngt = geantEvent->getNumberOfGeantTracks();
    mc_ngv = geantEvent->getNumberOfGeantVertices();
    
    cout << "MC Info: ";
    cout << "parti, signl, stSig, pileup, stPil, geantTr, geantVe" << endl; 
    printf("\t   %d\t %d\t%d\t%d\t%d\t%d\t%d\n", mc_np, mc_nsp, mc_nssp, mc_npp, mc_nspp, mc_ngt, mc_ngv);  

    
    // Filling the histograms
    
    // Clusters
    
    hClustersNumber.Fill((Axis_t) nc); 
    for (int ic = 0; ic < nc; ic++)
    {  
      tmpEmCluster = (XANAEmCluster *) clusters->At(ic);
      hClustersEnergy.Fill((Axis_t) tmpEmCluster->getEnergy());   
    }
    
    hSuperClustersNumber.Fill((Axis_t) ns);   
    for (int is = 0; is < ns; is++)
    {
      tmpSuperCluster = (XANASuperCluster *) superClusters->At(is);
      hSuperClustersEnergy.Fill((Axis_t) tmpSuperCluster->getEnergy());
      tmpSeedCluster = (XANAEmCluster *) tmpSuperCluster->getSeedCluster();
      hBremsRecovered.Fill((Axis_t) tmpSuperCluster->getEnergy() - tmpSeedCluster->getEnergy());
      nbrems = tmpSuperCluster->getNumberOfBrems();
      for (int ib = 0; ib < nbrems; ib++)
      {
        tmpBremCluster = (XANAEmCluster *)  tmpSuperCluster->getBremClusters()->At(ib);
        hBremsDetaDphi.Fill(tmpSeedCluster->getPosition().getEta()-tmpBremCluster->getPosition().getEta(),
                            tmpSeedCluster->getPosition().getPhi()-tmpBremCluster->getPosition().getPhi());
      }
    }
    
    hHadClustersNumber.Fill((Axis_t) nh);
    for (int ihc = 0; ihc < nhc; ihc++)
    {  
      tmpHadCluster = (XANAHadCluster *) hadClusters->At(ihc);
      hHadClustersEnergy.Fill((Axis_t) tmpHadCluster->getEnergy());   
    }
    
    hPreshClustersNumber.Fill((Axis_t) npc);  
    
    // Hits
    XANACaloRecHit *tmpCaloRecHit;
    
    hHitsNumber.Fill((Axis_t) nh); 
    for (int ih = 0; ih < nh; ih++)
    {  
      tmpCaloRecHit = (XANACaloRecHit *) hits->At(ih);
      hHitsEnergy.Fill((Axis_t) tmpCaloRecHit->getEnergy());   
    }
     
    // Tracks
   XANATrack *tmpTrack;
    
    hTracksNumber.Fill((Axis_t) nt); 
    for (int it = 0; it < nt; it++)
    {  
      tmpTrack = (XANATrack *) tracks->At(it);
      hTracksMomentumAtVertex.Fill((Axis_t) tmpTrack->getMomentumAtVertex().mag());
      hTracksMomentumAtFirstPoint.Fill((Axis_t) tmpTrack->getMomentumAtFirstPoint().mag()); 
      hTracksMomentumAtLastPoint.Fill((Axis_t) tmpTrack->getMomentumAtLastPoint().mag());  
    }    
      
    // Vertex 
    hVertexNumber.Fill((Axis_t) nv); 
    hVertexTracksNumber.Fill((Axis_t) nvt); 
     
    // Monte Carlo information
  
    // generator event
    hGenNumberOfParticles.Fill((Axis_t) mc_np);
     
#ifdef SINGLE_ELECTRON
    tmpGenElectron = (XANAGeneratorParticle *) generatorEvent->getSignalParticles()->At(0);
    hGenElectronMomentum.Fill((Axis_t) tmpGenElectron->getMomentum().vect().mag());
    hGenElectronPt.Fill((Axis_t) tmpGenElectron->getMomentum().perp());
    hGenElectronEta.Fill((Axis_t) tmpGenElectron->getMomentum().pseudoRapidity());
    hGenElectronPhi.Fill((Axis_t) tmpGenElectron->getMomentum().phi());
#endif    
   
    // GEANT event
    hGeantNumberOfTracks.Fill((Axis_t) mc_ngt);
    hGeantNumberOfVertices.Fill((Axis_t) mc_ngv);
      

     // electron candidates
    hElectronNumber.Fill((Axis_t) ne);
    if (ne != 0) {
      cout << "Electrons: Esc/p,   DetaSC,     HoE" << endl;   
      for (int ie = 0; ie < ne ; ie++) {
        tmpElectronCandidate = (XANAElectronCandidate *) electrons->At(ie);
        
        hElectronEscOverP.Fill((Axis_t) tmpElectronCandidate->getESuperClusterOverP());
        hElectronDetaSc.Fill((Axis_t) tmpElectronCandidate->getDeltaEtaSuperClusterAtVtx());
        hElectronHoE.Fill((Axis_t) tmpElectronCandidate->getHadronicOverEm());
        
        printf("  %d, \t%f, %f, %f \n", ie+1,
                tmpElectronCandidate->getESuperClusterOverP(),
                tmpElectronCandidate->getDeltaEtaSuperClusterAtVtx(),
                tmpElectronCandidate->getHadronicOverEm());  
      }
    }

    
    event->clear();
    generatorEvent->clear();
    geantEvent->clear(); 
  }
  
  cout << endl;
  cout << nevent << " events and " << nb << " bytes read " << endl;  
  cout << endl;
  
  f.Close();


// Exporting histograms to file
  TFile hFile(histoFileName.c_str(),"RECREATE");
  
  hClustersNumber.Write();
  hClustersEnergy.Write();
  hSuperClustersNumber.Write();
  hSuperClustersEnergy.Write();
  hHadClustersNumber.Write();
  hHadClustersEnergy.Write();
  hPreshClustersNumber.Write();
  hHitsNumber.Write();
  hHitsEnergy.Write();
  hTracksNumber.Write();
  hTracksMomentumAtVertex.Write();
  hTracksMomentumAtFirstPoint.Write();
  hTracksMomentumAtLastPoint.Write();
  hBremsRecovered.Write();
  hBremsDetaDphi.Write();
  hVertexNumber.Write();
  hVertexTracksNumber.Write();
  hGenNumberOfParticles.Write(); 
#ifdef SINGLE_ELECTRON
  hGenElectronMomentum.Write();
  hGenElectronPt.Write();
  hGenElectronEta.Write();
  hGenElectronPhi.Write();
#endif 
  hGeantNumberOfTracks.Write();
  hGeantNumberOfVertices.Write();
  hElectronNumber.Write();
  hElectronEscOverP.Write();
  hElectronDetaSc.Write();
  hElectronHoE.Write();  
    
  hFile.Close();
  
 
}

---