saturFamos.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 <XANADOO/XANAElectronCandidate/interface/XANAElectronTrack.h>
#include <XANADOO/XANAElectronCandidate/interface/XANAElectronSeed.h>
#include <XANADOO/XANAMuons/interface/XANAMuonCandidate.h>
#include <XANADOO/XANATriggerInfo/interface/XANATriggerInfo.h>
#include <XANADOO/XANAJets/interface/XANAJet.h>

#include <XANADOO/XANAGeometry/interface/XANAEcalGeometry.h>

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

#include <string>


using namespace std;

int main(int argc, char** argv)
{
// version TChain
   TChain *T = new TChain("ESD");
   T->Add("/data/scratch/collard/Production_xana151/single/famos/xSingleE4000.root");
   /*
   T->Add("/data/scratch/collard/Production_xana151/single/famos/xSinglePt50.root");
   T->Add("/data/scratch/collard/Production_xana151/single/famos/xSinglePt150.root");
   T->Add("/data/scratch/collard/Production_xana151/single/famos/xSinglePt500.root");
   T->Add("/data/scratch/collard/Production_xana151/single/famos/xSinglePt1500.root");
   T->Add("/data/scratch/collard/Production_xana151/single/famos/xSinglePt2500.root");
   */
 
// example of how to suppress reading of branches
// has to be put before setting the branch addresses!
// T->SetBranchStatus("GeantEvent*",0);

// Reconstructed event
  XANAEsdEvent *event = new XANAEsdEvent();  
  T->SetBranchAddress("Event.",&event);
  // MC information (generator and GEANT event)
  XANAGeneratorEvent *generatorEvent = new XANAGeneratorEvent();
  T->SetBranchAddress("GeneratorEvent.",&generatorEvent);
  XANAGeantEvent *geantEvent = new XANAGeantEvent();
  T->SetBranchAddress("GeantEvent.",&geantEvent);
  
  // Definitions
  Short_t nc = 0, nh = 0, nhc = 0, ns = 0, npc = 0,  nt = 0, nv = 0, nvt = 0;
  Short_t ne = 0, nm = 0, nelt = 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;
  TObjArray *electronTracks = 0;
  TObjArray *muons = 0;

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

  // Histo
  //
  //
  TH1F *histo1001 = new TH1F("histo1001","E rec / E true           ",100,0.8,1.05);

  TH2F *histo1002 = new TH2F("histo1002","S1 vs Red5  (all)",100,0.,350.,150,0.,4500.);
  TH2F *histo1003 = new TH2F("histo1003","S1 vs Red5 (S1<1.7TeV)",100,0.,175.,150,0.,1700.);
  TH1F *histo1004 = new TH1F("histo1004","E_SC corr / E_SC         ",100,0.5,1.5);

  TH1F *histo1011 = new TH1F("histo1011","S1/Red5 bin Red5:0-20",50,0.,24.0);
  TH1F *histo1012 = new TH1F("histo1012","S1/Red5 bin Red5:20-30",50,0.,24.0);
  TH1F *histo1013 = new TH1F("histo1013","S1/Red5 bin Red5:30-40",50,0.,24.0);
  TH1F *histo1014 = new TH1F("histo1014","S1/Red5 bin Red5:40-50",50,0.,24.0);
  TH1F *histo1015 = new TH1F("histo1015","S1/Red5 bin Red5:50-60",50,0.,24.0);
  TH1F *histo1016 = new TH1F("histo1016","S1/Red5 bin Red5:60-70",50,0.,24.0);
  TH1F *histo1017 = new TH1F("histo1017","S1/Red5 bin Red5:70-80",50,0.,24.0);
  TH1F *histo1018 = new TH1F("histo1018","S1/Red5 bin Red5:80-90",50,0.,24.0);
  TH1F *histo1019 = new TH1F("histo1019","S1/Red5 bin Red5:90-100",50,0.,24.0);
  TH1F *histo1020 = new TH1F("histo1020","S1/Red5 bin Red5:100-110",50,0.,24.0);
  TH1F *histo1021 = new TH1F("histo1021","S1/Red5 bin Red5:110-120",50,0.,24.0);

  TH1F *histo1111 = new TH1F("histo1111","S1 bin Red5:0-20",50,0.,1700.0);
  TH1F *histo1112 = new TH1F("histo1112","S1 bin Red5:20-30",50,0.,1700.0);
  TH1F *histo1113 = new TH1F("histo1113","S1 bin Red5:30-40",50,0.,1700.0);
  TH1F *histo1114 = new TH1F("histo1114","S1 bin Red5:40-50",50,0.,1700.0);
  TH1F *histo1115 = new TH1F("histo1115","S1 bin Red5:50-60",50,0.,1700.0);
  TH1F *histo1116 = new TH1F("histo1116","S1 bin Red5:60-70",50,0.,1700.0);
  TH1F *histo1117 = new TH1F("histo1117","S1 bin Red5:70-80",50,0.,1700.0);
  TH1F *histo1118 = new TH1F("histo1118","S1 bin Red5:80-90",50,0.,1700.0);
  TH1F *histo1119 = new TH1F("histo1119","S1 bin Red5:90-100",50,0.,1700.0);
  TH1F *histo1120 = new TH1F("histo1120","S1 bin Red5:100-110",50,0.,1700.0);
  TH1F *histo1121 = new TH1F("histo1121","S1 bin Red5:110-120",50,0.,1700.0);



  //   Start main loop on all events
  Int_t nb = 0; 
  Int_t nevent = (Int_t)T->GetEntries();

//  nevent =2;
  
  cout <<  "Nombre d entrees a traiter : " << nevent << endl; 
  for (Int_t i=0; i<nevent; i++) {
    nb += T->GetEntry(i); 


          // generator event
    mc_np   = generatorEvent->getNumberOfParticles();
    mc_nsp  = generatorEvent->getNumberOfSignalParticles();
    mc_nssp = generatorEvent->getNumberOfStableSignalParticles();
    mc_npp  = generatorEvent->getNumberOfPileupParticles();
    mc_nspp = generatorEvent->getNumberOfStablePileupParticles();
    TObjArray *genpart = 0;      genpart=generatorEvent->getParticles();
    int ilink=0; int first =1;
    for (int gp=0;gp<mc_np;gp++) {
     XANAGeneratorParticle *gpart = (XANAGeneratorParticle *)((*genpart)[gp]);
     if (gpart->getPartId()==11 && first==1)
     {
             ilink=gp;
             first=2;
     }
    }
    XANAGeneratorParticle *gElec = (XANAGeneratorParticle *)((*genpart)[ilink]);
    

    // Reconstructed event
    
    nc  = event->getNumberOfEmClusters();      clusters = event->getEmClusters(); 
    nh  = event->getNumberOfTrackHits();          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(); 
    nm  = event->getNumberOfMuonCandidates(); muons = event->getMuonCandidates(); 
    nelt    = event->getNumberOfElectronTracks();
    electronTracks = event->getElectronTracks();


    // electrons
    Int_t i_ze_el=-1;
    Float_t etmax=0.;

    for (Int_t ie=0; ie<ne; ie++) {
      XANAElectronCandidate *elec = (XANAElectronCandidate*)((*electrons)[ie]);
      if ( (fabs( elec->getSuperCluster()->getPosition().eta() ) <  1.4442 ||
            fabs( elec->getSuperCluster()->getPosition().eta() ) >  1.566 )
            && elec->getPt()/elec->getSuperCluster()->getEnergyScaleFactor() > etmax)
      {
            i_ze_el=ie;
            etmax=elec->getPt()/elec->getSuperCluster()->getEnergyScaleFactor();
      }
    }
    
    if (i_ze_el>-1)
    {
               XANAElectronCandidate *rElec = (XANAElectronCandidate*)((*electrons)[i_ze_el]);


               // E_rec/E_true
               float E_true = gElec->getMomentum().e();
               float E_rec = rElec->getSuperClusterEnergy()/rElec->getSuperCluster()->getEnergyScaleFactor();

               histo1001->Fill(E_rec/E_true);

               // S1 vs Red5
               float S1 = rElec->getSuperCluster()->getSum1();
               float Red5= rElec->getSuperCluster()->getSum9()-rElec->getSuperCluster()->getSum4();


               histo1002->Fill(Red5,S1);
               if (S1<1700)
               {

                histo1003->Fill(Red5,S1);
                if        (0<Red5 && Red5<=20)  { histo1011->Fill(S1/Red5); histo1111->Fill(S1); }
                else if  (20<Red5 && Red5<=30)  { histo1012->Fill(S1/Red5); histo1112->Fill(S1); }
                else if  (30<Red5 && Red5<=40)  { histo1013->Fill(S1/Red5); histo1113->Fill(S1); }
                else if  (40<Red5 && Red5<=50)  { histo1014->Fill(S1/Red5); histo1114->Fill(S1); }
                else if  (50<Red5 && Red5<=60)  { histo1015->Fill(S1/Red5); histo1115->Fill(S1); }
                else if  (60<Red5 && Red5<=70)  { histo1016->Fill(S1/Red5); histo1116->Fill(S1); }
                else if  (70<Red5 && Red5<=80)  { histo1017->Fill(S1/Red5); histo1117->Fill(S1); }
                else if  (80<Red5 && Red5<=90)  { histo1018->Fill(S1/Red5); histo1118->Fill(S1); }
                else if  (90<Red5 && Red5<=100) { histo1019->Fill(S1/Red5); histo1119->Fill(S1); }
                else if (100<Red5 && Red5<=110) { histo1020->Fill(S1/Red5); histo1120->Fill(S1); }
                else if (110<Red5 && Red5<=120) { histo1011->Fill(S1/Red5); histo1121->Fill(S1); }
                }
                else 
                {
                  // verif de la saturation
                  float NewS1= 10.9*Red5;
//                float NewS1= 11.0*Red5;
                  float NewErec= E_rec -S1 +NewS1;
                  histo1004->Fill(NewErec/E_rec);
                }


               
    }


    event->clear(); 
    generatorEvent->clear();
    geantEvent->clear();
  }
  
//  f.Close();
  //
     TFile *filehisto = new TFile("Saturation.root","RECREATE");

     histo1001->Write();
     histo1002->Write();
     histo1003->Write();
     histo1004->Write();

     histo1011->Write();
     histo1012->Write();
     histo1013->Write();
     histo1014->Write();
     histo1015->Write();
     histo1016->Write();
     histo1017->Write();
     histo1018->Write();
     histo1019->Write();
     histo1020->Write();
     histo1021->Write();
     histo1111->Write();
     histo1112->Write();
     histo1113->Write();
     histo1114->Write();
     histo1115->Write();
     histo1116->Write();
     histo1117->Write();
     histo1118->Write();
     histo1119->Write();
     histo1120->Write();
     histo1121->Write();


     filehisto->Close();

  
}

---