test_print.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("rsG1tev5_100.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);
  
  //   Start main loop on all events
  Int_t nb = 0; 
  Int_t nevent = (Int_t)T->GetEntries();

  nevent =5;

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

    cout << " Graviton for event " << i << endl;

       // 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,iRSgrav=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 RSgrav " << endl;
               }
         } 
    }
    XANAGeneratorParticle *gRSgrav = (XANAGeneratorParticle *)((*genpart)[iRSgrav]);
    XANAGeneratorParticle *gElec1 = (XANAGeneratorParticle *)((*genpart)[ilink1]);
    XANAGeneratorParticle *gElec2 = (XANAGeneratorParticle *)((*genpart)[ilink2]);

    float trueMassGen=gRSgrav->getMass();
    float eeMassGen=RSMassGen(gElec1,gElec2);
    float costhetaGen = RSdistribAngGen(gElec1,gElec2);

    float eeMassRec=-1.;
    float costhetaRec = -2.;


    // 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();

               cout << "RSSelection: elec1 " << energy1 << " Eta " << eta1 << " Phi " << phi1 << endl;
               cout << "RSSelection: elec2 " << energy2 << " Eta " << eta2 << " Phi " << phi2 << endl;
              
               // Traitement de la saturation
               int iflag_satur=0;
               RSSaturation(elec1, energy1, iflag_satur);
               cout << "RSSaturation: elec1 " << energy1 << " Sat? " << iflag_satur << endl;

               RSSaturation(elec2, energy2, iflag_satur);              
               cout << "RSSaturation: elec2 " << energy2 << " Sat? " << iflag_satur << endl;
               
               // Association des photons
               float EtIsol1=0.;
               float EtIsol2=0.;
               RSAssocPhoton(elec1, energy1, EtIsol1);
               RSAssocPhoton(elec2, energy2, EtIsol2);

               cout << "RSAssocPhoton: elec1 " << energy1 << " Iso " << EtIsol1 << endl; 
               cout << "RSAssocPhoton: elec2 " << energy2 << " Iso " << EtIsol2 << endl;
               
               // 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;

               cout << "RSCorrectionFunc: elec1 " << energy1 << " xcor1 " << xcor1 << endl; 
               cout << "RSCorrectionFunc: elec2 " << energy2 << " xcor2 " << xcor2 << endl;
               
               // Candidat RSElec
               
               RSElec *cElec1 = new RSElec(elec1,energy1,eta1,phi1,EtIsol1);
               RSElec *cElec2 = new RSElec(elec2,energy2,eta2,phi2,EtIsol2);

               cout << " RSElec: cElec1 " << cElec1->getElecEnergy() << " Eta" << cElec1->getElecEta() <<
                       " Phi " << cElec1->getElecPhi() << " Iso " << cElec1->getElecIsol() << " Px " <<  cElec1->getElecPx() <<
                       " Py " << cElec1->getElecPy() << " Pz " << cElec1->getElecPz() << " Et " <<  cElec1->getElecEt() << endl;
               cout << " RSElec: cElec2 " << cElec2->getElecEnergy() << " Eta" << cElec2->getElecEta() <<
                       " Phi " << cElec2->getElecPhi() << " Iso " << cElec2->getElecIsol() << " Px " <<  cElec2->getElecPx() <<
                       " Py " << cElec2->getElecPy() << " Pz " << cElec2->getElecPz() << " Et " <<  cElec2->getElecEt() << endl;

                // Masse Reconstruite
                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() )
                    {
                        costhetaRec = RSdistribAngRec(cElec1,cElec2);

                    }


    cout << " Graviton Mass " << trueMassGen << " ee-gen " << eeMassGen << " rec " << eeMassRec << endl;
    cout << " Graviton costheta : ee-gen " << costhetaGen << " rec " << costhetaRec << endl; 

    cout << endl;


    }


    event->clear(); 
    generatorEvent->clear();
    geantEvent->clear();
  }
  
}

---