zprimAnal.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 <XANADOO/XANAnalysisRS/interface/RSElecSelector.h>
#include <XANADOO/XANAnalysisRS/interface/RSCorrectionFunc.h>
#include <XANADOO/XANAnalysisRS/interface/RSSaturation.h>
#include <XANADOO/XANAnalysisRS/interface/RSPhotonAssoc.h>
#include <XANADOO/XANAnalysisRS/interface/ClassRSElec.h>
#include <XANADOO/XANAnalysisRS/interface/RSVertexCorr.h>
#include <XANADOO/XANAnalysisRS/interface/RSDistribAng.h>
#include <XANADOO/XANAnalysisRS/interface/RSMassEstim.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;

Short_t ne_candidatelectron=0,ns_supercluster=0;
TObjArray *genpart = 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;


  
int main(int argc, char** argv)
{
  //   Connect to file generated
  TFile *f;
    //  if (argc==1) f = new TFile("SE04_E25_Pall_PT30.root");
  if (argc==1) f = new TFile("/data/scratch/collard/Production_xana151/zprim/famos/xZprimALRM_m1tev.root");

  else f = new TFile(*++argv);
 
    //   Read Tree named "T" in memory. Tree pointer is assigned the same name
  TTree *T = (TTree*)f->Get("ESD");

   
// version TChain
//   TChain *T = new TChain("ESD");
//   T->Add("file1.root");
//   T->Add("file2.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);
  


 
  // Histo
  // sans contrainte de masse
  TH1F *histo1000 = new TH1F("histo1000","M Zprime gen ee ",280,0.,7.);
  TH1F *histo1001 = new TH1F("histo1001","M Zprime gen",280,0.,7.);
  TH1F *histo1002 = new TH1F("histo1002","M Zprime rec",280,0.,7.);
  TH1F *histo1003 = new TH1F("histo1003","Delta M     ",100,-50.,50.);
  TH1F *histo1004 = new TH1F("histo1004","CosTheta gen",20,-1.,1.);
  TH1F *histo1005 = new TH1F("histo1005","CosTheta rec",20,-1.,1.);
  TH1F *histo1006 = new TH1F("histo1006","CosTheta gen",10,-1.,1.);
  TH1F *histo1007 = new TH1F("histo1007","CosTheta rec",10,-1.,1.);
  TH1F *histo1008 = new TH1F("histo1008","Delta Cos th",100,-0.5,0.5);
  TH1F *histo1009 = new TH1F("histo1009","Delta Cos th",100,-0.2,0.2);
  
  // 1 TeV 
  TH1F *histo2000 = new TH1F("histo2000","M Zprime gen ee ",24,0.4,1.6);
  TH1F *histo2001 = new TH1F("histo2001","M Zprime gen ",24,0.4,1.6);
  TH1F *histo2002 = new TH1F("histo2002","M Zprime rec ",24,0.4,1.6);
  TH1F *histo2003 = new TH1F("histo2003","M Zprime gen",32,0.5,1.3);
  TH1F *histo2004 = new TH1F("histo2004","M Zprime rec",32,0.5,1.3);
  TH1F *histo2005 = new TH1F("histo2005","M Zprime gen",80,0.5,1.3);
  TH1F *histo2006 = new TH1F("histo2006","M Zprime rec",80,0.5,1.3);
  TH1F *histo2007 = new TH1F("histo2007","M Zprime gen",160,0.5,1.3);
  TH1F *histo2008 = new TH1F("histo2008","M Zprime rec",160,0.5,1.3);
  TH1F *histo2009 = new TH1F("histo2009","M Zprime gen",800,0.5,1.3);
  TH1F *histo2010 = new TH1F("histo2010","M Zprime rec",800,0.5,1.3);
  TH1F *histo2011 = new TH1F("histo2011","CosTheta gen",20,-1.,1.);
  TH1F *histo2012 = new TH1F("histo2012","CosTheta rec",20,-1.,1.);
  TH1F *histo2013 = new TH1F("histo2013","CosTheta gen",10,-1.,1.);
  TH1F *histo2014 = new TH1F("histo2014","CosTheta rec",10,-1.,1.);

  // 3 TeV 
  TH1F *histo3000 = new TH1F("histo3000","M Zprime gen ee ",30,1.5,4.5);
  TH1F *histo3001 = new TH1F("histo3001","M Zprime gen ",30,1.5,4.5);
  TH1F *histo3002 = new TH1F("histo3002","M Zprime rec ",30,1.5,4.5);
  TH1F *histo3003 = new TH1F("histo3003","M Zprime gen",40,1.8,3.8);
  TH1F *histo3004 = new TH1F("histo3004","M Zprime rec",40,1.8,3.8);
  TH1F *histo3005 = new TH1F("histo3005","M Zprime gen",100,1.8,3.8);
  TH1F *histo3006 = new TH1F("histo3006","M Zprime rec",100,1.8,3.8);
  TH1F *histo3007 = new TH1F("histo3007","M Zprime gen",200,1.8,3.8);
  TH1F *histo3008 = new TH1F("histo3008","M Zprime rec",200,1.8,3.8);
  TH1F *histo3009 = new TH1F("histo3009","M Zprime gen",2000,1.8,3.8);
  TH1F *histo3010 = new TH1F("histo3010","M Zprime rec",2000,1.8,3.8);
  TH1F *histo3011 = new TH1F("histo3011","CosTheta gen",20,-1.,1.);
  TH1F *histo3012 = new TH1F("histo3012","CosTheta rec",20,-1.,1.);
  TH1F *histo3013 = new TH1F("histo3013","CosTheta gen",10,-1.,1.);
  TH1F *histo3014 = new TH1F("histo3014","CosTheta rec",10,-1.,1.);

  // 5 TeV 
  TH1F *histo4000 = new TH1F("histo4000","M Zprime gen ee ",40,3.0,7.0);
  TH1F *histo4001 = new TH1F("histo4001","M Zprime gen ",40,3.0,7.0);
  TH1F *histo4002 = new TH1F("histo4002","M Zprime rec ",40,3.0,7.0);
  TH1F *histo4003 = new TH1F("histo4003","M Zprime gen",50,3.5,6.0);
  TH1F *histo4004 = new TH1F("histo4004","M Zprime rec",50,3.5,6.0);
  TH1F *histo4005 = new TH1F("histo4005","M Zprime gen",125,3.5,6.0);
  TH1F *histo4006 = new TH1F("histo4006","M Zprime rec",125,3.5,6.0);
  TH1F *histo4007 = new TH1F("histo4007","M Zprime gen",250,3.5,6.0);
  TH1F *histo4008 = new TH1F("histo4008","M Zprime rec",250,3.5,6.0);
  TH1F *histo4009 = new TH1F("histo4009","M Zprime gen",2500,3.5,6.0);
  TH1F *histo4010 = new TH1F("histo4010","M Zprime rec",2500,3.5,6.0);
  TH1F *histo4011 = new TH1F("histo4011","CosTheta gen",20,-1.,1.);
  TH1F *histo4012 = new TH1F("histo4012","CosTheta rec",20,-1.,1.);
  TH1F *histo4013 = new TH1F("histo4013","CosTheta gen",10,-1.,1.);
  TH1F *histo4014 = new TH1F("histo4014","CosTheta rec",10,-1.,1.);

  // 5 TeV  sans saturation
  TH1F *histo5000 = new TH1F("histo5000","M Zprime gen ee ",40,3.0,7.0);
  TH1F *histo5001 = new TH1F("histo5001","M Zprime gen ",40,3.0,7.0);
  TH1F *histo5002 = new TH1F("histo5002","M Zprime rec ",40,3.0,7.0);
  TH1F *histo5003 = new TH1F("histo5003","M Zprime gen",50,3.5,6.0);
  TH1F *histo5004 = new TH1F("histo5004","M Zprime rec",50,3.5,6.0);
  TH1F *histo5005 = new TH1F("histo5005","M Zprime gen",125,3.5,6.0);
  TH1F *histo5006 = new TH1F("histo5006","M Zprime rec",125,3.5,6.0);
  TH1F *histo5007 = new TH1F("histo5007","M Zprime gen",250,3.5,6.0);
  TH1F *histo5008 = new TH1F("histo5008","M Zprime rec",250,3.5,6.0);
  TH1F *histo5009 = new TH1F("histo5009","M Zprime gen",2500,3.5,6.0);
  TH1F *histo5010 = new TH1F("histo5010","M Zprime rec",2500,3.5,6.0);
  TH1F *histo5011 = new TH1F("histo5011","CosTheta gen",20,-1.,1.);
  TH1F *histo5012 = new TH1F("histo5012","CosTheta rec",20,-1.,1.);
  TH1F *histo5013 = new TH1F("histo5013","CosTheta gen",10,-1.,1.);
  TH1F *histo5014 = new TH1F("histo5014","CosTheta rec",10,-1.,1.);


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

//  nevent =10;

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


       // Monte Carlo information
   
      // generator event
    Short_t mc_np = 0, mc_nsp = 0, mc_nssp = 0, mc_npp = 0, mc_nspp = 0 ;
    mc_np   = generatorEvent->getNumberOfParticles();
    mc_nsp  = generatorEvent->getNumberOfSignalParticles();
    mc_nssp = generatorEvent->getNumberOfStableSignalParticles();
    mc_npp  = generatorEvent->getNumberOfPileupParticles();
    mc_nspp = generatorEvent->getNumberOfStablePileupParticles();
    genpart=generatorEvent->getParticles();
    int npart=0,izprime=11,ilink1=0,ilink2=0;
    for (int gp=0;gp<mc_np;gp++) {
         XANAGeneratorParticle *gpart = (XANAGeneratorParticle *)((*genpart)[gp]);
         if ((gpart->getPartId()==11 || gpart->getPartId()==-11) &&  gpart->getMo1()==12 )
         { 
               if (npart==0)
               {
                    // elec 1
                    npart++;
                    ilink1=gp;
               }
               else if (npart==1)
               {
                    //elec 2
                   npart++;
                   ilink2=gp;
               }
               else 
               {
                   cout << " pas possible d'avoir 3 electrons rattaches au Zprime " << endl;
               }
         } 
    }
    XANAGeneratorParticle *gZprim = (XANAGeneratorParticle *)((*genpart)[izprime]);
    XANAGeneratorParticle *gElec1 = (XANAGeneratorParticle *)((*genpart)[ilink1]);
    XANAGeneratorParticle *gElec2 = (XANAGeneratorParticle *)((*genpart)[ilink2]);

    float ZMassGen=gZprim->getMass();
    float eeMassGen=RSMassGen(gElec1,gElec2);
    float costhetaGen = RSdistribAngGen(gElec1,gElec2);

    histo1000->Fill( eeMassGen/1000.);
    histo2000->Fill( eeMassGen/1000.);
    histo3000->Fill( eeMassGen/1000.);
    histo4000->Fill( eeMassGen/1000.);
    histo5000->Fill( eeMassGen/1000.);

    histo1001->Fill( ZMassGen/1000.);

    histo2001->Fill( ZMassGen/1000.);
    histo2003->Fill( ZMassGen/1000.);
    histo2005->Fill( ZMassGen/1000.);
    histo2007->Fill( ZMassGen/1000.);
    histo2009->Fill( ZMassGen/1000.);
    
    histo3001->Fill( ZMassGen/1000.);
    histo3003->Fill( ZMassGen/1000.);
    histo3005->Fill( ZMassGen/1000.);
    histo3007->Fill( ZMassGen/1000.);
    histo3009->Fill( ZMassGen/1000.);

    histo4001->Fill( ZMassGen/1000.);
    histo4003->Fill( ZMassGen/1000.);
    histo4005->Fill( ZMassGen/1000.);
    histo4007->Fill( ZMassGen/1000.);
    histo4009->Fill( ZMassGen/1000.);
   
    histo5001->Fill( ZMassGen/1000.);
    histo5003->Fill( ZMassGen/1000.);
    histo5005->Fill( ZMassGen/1000.);
    histo5007->Fill( ZMassGen/1000.);
    histo5009->Fill( ZMassGen/1000.);

    histo1004->Fill( costhetaGen );
    histo1006->Fill( costhetaGen );

    // Reconstructed event
    
    // Definitions
    Short_t nc = 0, nh = 0, nhc = 0, npc = 0,  nt = 0, nv = 0, nvt = 0;
    Short_t nm = 0, nelt = 0 ;    
  

    nc  = event->getNumberOfEmClusters();         clusters = event->getEmClusters(); 
    nh  = event->getNumberOfTrackHits();          hits = event->getEmClusterHits();   
    nhc = event->getNumberOfHadClusters();        hadClusters = event->getHadClusters();   
    ns_supercluster  = 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_candidatelectron  = event->getNumberOfElectronCandidates(); 
                                                  electrons = event->getElectronCandidates(); 
    nm  = event->getNumberOfMuonCandidates();     muons = event->getMuonCandidates(); 
    nelt    = event->getNumberOfElectronTracks(); electronTracks = event->getElectronTracks();


    
    // electrons
    Int_t indexel1=-1, indexel2=-1;
    
    // Selection des 2 candidats electrons les plus energetiques

    RSSelection(indexel1, indexel2);

    // On ne travaille qu'avec les 2 Candidats Electrons les + energetiques
    
    if (indexel1>-1 && indexel2>-1)
    {
               XANAElectronCandidate *elec1 = (XANAElectronCandidate*)((*electrons)[indexel1]);
               XANAElectronCandidate *elec2 = (XANAElectronCandidate*)((*electrons)[indexel2]);
              
               float energy1 =elec1->getSuperClusterEnergy()/elec1->getSuperCluster()->getEnergyScaleFactor();
               float eta1    =elec1->getSuperCluster()->getPosition().eta();
               float phi1    =elec1->getSuperCluster()->getPosition().phi();
               
               float energy2 =elec2->getSuperClusterEnergy()/elec2->getSuperCluster()->getEnergyScaleFactor();
               float eta2    =elec2->getSuperCluster()->getPosition().eta();
               float phi2    =elec2->getSuperCluster()->getPosition().phi();
              
               // Traitement de la saturation
               int iflag_satur=0;
               RSSaturation(elec1, energy1, iflag_satur);
               RSSaturation(elec2, energy2, iflag_satur);              
               
               // Association des photons
               float EtIsol1=0.;
               float EtIsol2=0.;
               RSAssocPhoton(elec1, energy1, EtIsol1);
               RSAssocPhoton(elec2, energy2, EtIsol2);
               
               // Correction de l'energie pour eta
               
               bool FAMOS=true;
               float xcor1=RSCorrectionFunc(energy1,  eta1,  FAMOS);
               float xcor2=RSCorrectionFunc(energy2,  eta2,  FAMOS);
               xcor1=1./xcor1;
               xcor2=1./xcor2;
               energy1*=xcor1;
               energy2*=xcor2;
               
               // Candidat RSElec
               
               RSElec *cElec1 = new RSElec(elec1,energy1,eta1,phi1,EtIsol1);
               RSElec *cElec2 = new RSElec(elec2,energy2,eta2,phi2,EtIsol2);

                // Masse Reconstruite
                float eeMassRec=RSMassRec(cElec1,cElec2);
                
                // Cut en iso (H/E pas encore dispo dans FAMOS 1_2_0)
                if ( cElec1->getElecIsol()< 0.02* cElec1->getElecEt() &&
                     cElec2->getElecIsol()< 0.02* cElec2->getElecEt() )
                    {
                        float costhetaRec = RSdistribAngRec(cElec1,cElec2);

                        histo1002->Fill( eeMassRec /1000. );
                        histo1003->Fill( ZMassGen-eeMassRec);
                        
                        histo1005->Fill(costhetaRec);
                        histo1007->Fill(costhetaRec);
                        histo1008->Fill( costhetaGen-costhetaRec);
                        histo1009->Fill( costhetaGen-costhetaRec);

                        histo2002->Fill( eeMassRec /1000. );
                        histo2004->Fill( eeMassRec /1000. );
                        histo2006->Fill( eeMassRec /1000. );
                        histo2008->Fill( eeMassRec /1000. );
                        histo2010->Fill( eeMassRec /1000. );
                        
                        histo3002->Fill( eeMassRec /1000. );
                        histo3004->Fill( eeMassRec /1000. );
                        histo3006->Fill( eeMassRec /1000. );
                        histo3008->Fill( eeMassRec /1000. );
                        histo3010->Fill( eeMassRec /1000. );

                        histo4002->Fill( eeMassRec /1000. );
                        histo4004->Fill( eeMassRec /1000. );
                        histo4006->Fill( eeMassRec /1000. );
                        histo4008->Fill( eeMassRec /1000. );
                        histo4010->Fill( eeMassRec /1000. );

                        if (iflag_satur==0)
                        {
                          histo5002->Fill( eeMassRec /1000. );
                          histo5004->Fill( eeMassRec /1000. );
                          histo5006->Fill( eeMassRec /1000. );
                          histo5008->Fill( eeMassRec /1000. );
                          histo5010->Fill( eeMassRec /1000. );
                        }

                        if (900<eeMassRec && eeMassRec<1100)
                        {
                          histo2011->Fill( costhetaGen );
                          histo2012->Fill( costhetaRec );
                          histo2013->Fill( costhetaGen );
                          histo2014->Fill( costhetaRec );
                        }
                        else if (2750<eeMassRec && eeMassRec<3250)
                        {
                          histo3011->Fill( costhetaGen );
                          histo3012->Fill( costhetaRec );
                          histo3013->Fill( costhetaGen );
                          histo3014->Fill( costhetaRec );
                        }
                        else if (4500<eeMassRec && eeMassRec<5500)
                        {
                          histo4011->Fill( costhetaGen );
                          histo4012->Fill( costhetaRec );
                          histo4013->Fill( costhetaGen );
                          histo4014->Fill( costhetaRec );
                          if (iflag_satur==0)
                          {
                            histo5011->Fill( costhetaGen );
                            histo5012->Fill( costhetaRec );
                            histo5013->Fill( costhetaGen );
                            histo5014->Fill( costhetaRec );
                          }
                        }



                    }



    }


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

     histo1000->Write();
     histo1001->Write();
     histo1002->Write();
     histo1003->Write();
     histo1004->Write();
     histo1005->Write();
     histo1006->Write();
     histo1007->Write();
     histo1008->Write();
     histo1009->Write();


     histo2000->Write();
     histo2001->Write();
     histo2002->Write();
     histo2003->Write();
     histo2004->Write();
     histo2005->Write();
     histo2006->Write();
     histo2007->Write();
     histo2008->Write();
     histo2009->Write();
     histo2010->Write();
     histo2011->Write();
     histo2012->Write();
     histo2013->Write();
     histo2014->Write();

     histo3000->Write();
     histo3001->Write();
     histo3002->Write();
     histo3003->Write();
     histo3004->Write();
     histo3005->Write();
     histo3006->Write();
     histo3007->Write();
     histo3008->Write();
     histo3009->Write();
     histo3010->Write();
     histo3011->Write();
     histo3012->Write();
     histo3013->Write();
     histo3014->Write();

     histo4000->Write();
     histo4001->Write();
     histo4002->Write();
     histo4003->Write();
     histo4004->Write();
     histo4005->Write();
     histo4006->Write();
     histo4007->Write();
     histo4008->Write();
     histo4009->Write();
     histo4010->Write();
     histo4011->Write();
     histo4012->Write();
     histo4013->Write();
     histo4014->Write();

     histo5000->Write();
     histo5001->Write();
     histo5002->Write();
     histo5003->Write();
     histo5004->Write();
     histo5005->Write();
     histo5006->Write();
     histo5007->Write();
     histo5008->Write();
     histo5009->Write();
     histo5010->Write();
     histo5011->Write();
     histo5012->Write();
     histo5013->Write();
     histo5014->Write();

     
     filehisto->Close();

  
}

---