source: svn/trunk/Examples/src/Analysis_Ex.cc@ 820

/***********************************************************************
**                                                                    **
**          /----------------------------------------------\          **
**         |  Delphes, a framework for the fast simulation  |         **
**         |       of a  generic collider experiment        |         **
**          \-------------  arXiv:0903.2225v1  ------------/          **
**                                                                    **
**                                                                    **
**   This package uses:                                               **
**   ------------------                                               **
**   ROOT: Nucl. Inst. & Meth. in Phys. Res. A389 (1997) 81-86        **
**   FastJet algorithm: Phys. Lett. B641 (2006) [hep-ph/0512210]      **
**   Hector: JINST 2:P09005 (2007) [physics.acc-ph:0707.1198v2]       **
**   FROG: [hep-ex/0901.2718v1]                                       **
**   HepMC: Comput. Phys. Commun.134 (2001) 41                        **
**                                                                    **
** ------------------------------------------------------------------ **
**                                                                    **
**   Main authors:                                                    **
**   -------------                                                    **
**                                                                    **
**                Severine Ovyn                Xavier Rouby           **
**          severine.ovyn@uclouvain.be      xavier.rouby@cern         **
**                                                                    **
**         Center for Particle Physics and Phenomenology (CP3)        **
**                Universite catholique de Louvain (UCL)              **
**                     Louvain-la-Neuve, Belgium                      **
**                                                                    **
**                      Copyright (C) 2008-2009,                      **
**                        All rights reserved.                        **
**                                                                    **
***********************************************************************/

#include "Examples/interface/Analysis_Ex.h"
#include <iostream>
#include <sstream>
#include <fstream>
#include <iomanip>

using namespace std;

//******************************Debut de l'analyse****************************************
//*****************************************************************************************

Analysis_Ex::Analysis_Ex(string CardWithCuts,string LogName)
{
  string temp_string;
  istringstream curstring;
  
  ifstream fichier_a_lire(CardWithCuts.c_str());
  if(!fichier_a_lire.good()) {
    cout << "DataCardname " << CardWithCuts << " not found, use default values" << endl;
    return;
  }
  
  while (getline(fichier_a_lire,temp_string)) {
    curstring.clear(); // needed when using several times istringstream::str(string)
    curstring.str(temp_string);
    string varname;
    float value;
    
    if(strstr(temp_string.c_str(),"#")) { }//remove comments
    else if(strstr(temp_string.c_str(),"PT_ELEC")){curstring >> varname >> value; PT_ELEC = value;}
    else if(strstr(temp_string.c_str(),"PT_MUON")){curstring >> varname >> value; PT_MUON = value;}
    else if(strstr(temp_string.c_str(),"INV_MASS_LL")){curstring >> varname >> value; INV_MASS_LL = value;}
  }
  
  ofstream f_out(LogName.c_str(),ofstream::app);
  
  f_out<<"*******************************************************************"<<endl;
  f_out << left << setw(30) <<"Cut values used in the analysis: "<<""
        << right << setw(37) <<"-------------------------------------"<<"\n";
  f_out << left <<setw(50) << "Invariant mass of the leptons: "<<""
        << right <<setw(17) << INV_MASS_LL <<"\n";
  f_out<<"*******************************************************************"<<endl;
  
}

Analysis_Ex::~Analysis_Ex()
{
}

void Analysis_Ex::Run(ExRootTreeReader *treeReaderGen, ExRootTreeReader *treeReaderRec, ExRootTreeReader *treeReaderTrig, ExRootTreeWriter *treeWriter)
{ 
  total=0;//initialisation of total number of events
  cut_trig=0;
  cut_1=0;//initialisation of counter for cut 1
  cut_2=0;
  //access the branches************************
  //to get the generator level information
  const TClonesArray *GEN = treeReaderGen->UseBranch("Particle");
  
  //to get the reconstructed level information
  const TClonesArray *JET    = treeReaderRec->UseBranch("Jet");
  const TClonesArray *TAUJET = treeReaderRec->UseBranch("TauJet");
  const TClonesArray *PHOTO  = treeReaderRec->UseBranch("Photon");
  const TClonesArray *ELEC   = treeReaderRec->UseBranch("Electron");
  const TClonesArray *MUON   = treeReaderRec->UseBranch("Muon");
  const TClonesArray *TRACKS = treeReaderRec->UseBranch("Tracks");
  const TClonesArray *CALO   = treeReaderRec->UseBranch("CaloTower");
  
  //to get the VFD reconstructed level information
  const TClonesArray *ZDC    = treeReaderRec->UseBranch("ZDChits");
  const TClonesArray *RP220 = treeReaderRec->UseBranch("RP220hits");
  const TClonesArray *FP420 = treeReaderRec->UseBranch("FP420hits");
  
  //to get the trigger information
  const TClonesArray *TRIGGER  = treeReaderTrig->UseBranch("TrigResult");

  //Define the branches that will be filled during the analysis
  ExRootTreeBranch *INVMASS = treeWriter->NewBranch("INVMass", TRootInvm::Class());
  TRootInvm *inv_mass; 
  //*******************************************
  
  //run on the events
  Long64_t entry, allEntries = treeReaderRec->GetEntries();
  cout << "** Chain contains " << allEntries << " events" << endl;
  total=allEntries;

  //general information
  float E,Px,Py,Pz;
  float PT,Eta,Phi;

  //lepton information
  bool IsolFlag;
 
  //bjet information
  bool Btag;

  //Particle level information
  int PID, Status, M1,M2,D1,D2;
  float Charge, T, X, Y, Z, M;

  //VFD information
  float S,q2,Tx,Ty;  
  int side;
  bool hadronic_hit; 
        
  for(entry = 0; entry < allEntries; ++entry)
    {
      treeReaderGen->ReadEntry(entry);//access information of generated information
      treeReaderRec->ReadEntry(entry);//access information of reconstructed information
      treeReaderTrig->ReadEntry(entry);//access information of Trigger information
      
      //*****************************************************
      //Example how to run on the generator level information
      //*****************************************************
      TIter itGen((TCollection*)GEN);
      TRootC::GenParticle *gen;
      itGen.Reset();
      while( (gen = (TRootC::GenParticle*) itGen.Next()) )
        {
          PID      = gen->PID;      // particle HEP ID number 
          Status   = gen->Status;   // particle status 
          M1       = gen->M1;       // particle 1st mother 
          M2       = gen->M2;       // particle 2nd mother 
          D1       = gen->D1;       // particle 1st daughter 
          D2       = gen->D2;       // particle 2nd daughter 
          Charge = gen->Charge;   // electrical charge
          
          T      = gen->T;        // particle vertex position (t component)
          X      = gen->X;        // particle vertex position (x component) 
          Y      = gen->Y;        // particle vertex position (y component) 
          Z      = gen->Z;        // particle vertex position (z component)
          M      = gen->M;        // particle mass
        }
      
      
      //***********************************************
      //Example how to run on the reconstructed objects
      //***********************************************
      
      //access the Electron branch
      TIter itElec((TCollection*)ELEC);
      TRootElectron *elec;
      itElec.Reset();
      while( (elec = (TRootElectron*) itElec.Next()) )
        {
          E        = elec->E;        // particle energy in GeV
          Px       = elec->Px;       // particle momentum vector (x component) in GeV
          Py       = elec->Py;       // particle momentum vector (y component) in GeV
          Pz       = elec->Pz;       // particle momentum vector (z component) in GeV
          
          PT       = elec->PT;       // particle transverse momentum in GeV
          Eta      = elec->Eta;      // particle pseudorapidity
          Phi      = elec->Phi;      // particle azimuthal angle in rad
          IsolFlag = elec->IsolFlag;  // is the particule isolated?
        }
      //Running on the muon branch is identical:
      TIter itMuon((TCollection*)MUON);
      TRootMuon *muon;
      itMuon.Reset();
      while( (muon = (TRootMuon*) itMuon.Next()) ){}
      
      //access the Photon branch
      TIter itGam((TCollection*)PHOTO);
      TRootPhoton *gam;
      itGam.Reset();
      while( (gam = (TRootPhoton*) itGam.Next()) )
        {
          E        = gam->E;        // particle energy in GeV
          Px       = gam->Px;       // particle momentum vector (x component) in GeV
          Py       = gam->Py;       // particle momentum vector (y component) in GeV
          Pz       = gam->Pz;       // particle momentum vector (z component) in GeV
          
          PT       = gam->PT;       // particle transverse momentum in GeV
          Eta      = gam->Eta;      // particle pseudorapidity
          Phi      = gam->Phi;      // particle azimuthal angle in rad
        }
      
      //access the jet branch
      TIter itJet((TCollection*)JET);
      TRootJet *jet;
      itJet.Reset();
      while( (jet = (TRootJet*) itJet.Next()) )
        {
          E        = jet->E;        // particle energy in GeV
          Px       = jet->Px;       // particle momentum vector (x component) in GeV
          Py       = jet->Py;       // particle momentum vector (y component) in GeV
          Pz       = jet->Pz;       // particle momentum vector (z component) in GeV
          
          PT       = jet->PT;       // particle transverse momentum in GeV
          Eta      = jet->Eta;      // particle pseudorapidity
          Phi      = jet->Phi;      // particle azimuthal angle in rad
          Btag     = jet->Btag;     // is the jet BTagged 
        }
      //Running on the tau-jet branch is identical:
      TIter itTaujet((TCollection*)TAUJET);
      TRootTauJet *taujet;
      itTaujet.Reset();
      while( (taujet = (TRootTauJet*) itTaujet.Next()) ){}
      
      //access the track branch
      TIter itTrack((TCollection*)TRACKS);
      TRootTracks *tracks;
      itTrack.Reset();
      while( (tracks = (TRootTracks*) itTrack.Next()) )
        {
          E        = tracks->E;        // particle energy in GeV
          Px       = tracks->Px;       // particle momentum vector (x component) in GeV
          Py       = tracks->Py;       // particle momentum vector (y component) in GeV
          Pz       = tracks->Pz;       // particle momentum vector (z component) in GeV
          
          PT       = tracks->PT;       // particle transverse momentum in GeV
          Eta      = tracks->Eta;      // particle pseudorapidity
          Phi      = tracks->Phi;      // particle azimuthal angle in rad
        }
      
      //Running on the calo branch is identical:
      TIter itCalo((TCollection*)CALO);
      TRootCalo *calo;
      itCalo.Reset();
      while( (calo = (TRootCalo*) itCalo.Next()) ){}
      
      //***************************************************
      //Example how to run on the VFD reconstructed objects
      //***************************************************
      
      //access the ZDC branch
      TIter itZdc((TCollection*)ZDC);
      TRootZdcHits *zdc;
      itZdc.Reset();
      while( (zdc = (TRootZdcHits*) itZdc.Next()) )
        {
          E   = zdc->E;     // particle energy in GeV
          T   = zdc->T;     // time of flight [s]
        /*
          Px  = zdc->Px;    // particle momentum vector (x component) in GeV
          Py  = zdc->Py;    // particle momentum vector (y component) in GeV
          Pz  = zdc->Pz;    // particle momentum vector (z component) in GeV
          
          PT  = zdc->PT;    // particle transverse momentum in GeV
          Eta = zdc->Eta;   // particle pseudorapidity
          Phi = zdc->Phi;   // particle azimuthal angle in rad
        */
          side  = zdc->side;  // -1 or +1
          //hadronic_hit = zdc->hadronic_hit; // true if neutron, false if photon
        }
      
      //access the RP220 branch
      TIter itRp220((TCollection*)RP220);
      TRootRomanPotHits *rp220;
      itRp220.Reset();
      
      while( (rp220 = (TRootRomanPotHits*) itRp220.Next()) )
        {
          //T   = rp220->T;     // time of flight to the detector [s]
          S   = rp220->S;     // distance to the IP [m]
          E   = rp220->E;     // reconstructed energy [GeV]
          q2  = rp220->q2;    // reconstructed squared momentum transfer [GeV^2]
          
          X   = rp220->X;     // horizontal distance to the beam [um]
          Y   = rp220->Y;     // vertical distance to the beam [um]
          
          Tx  = rp220->Tx;    // angle of the momentum in the horizontal (x,z) plane [urad]
          Ty  = rp220->Ty;    // angle of the momentum in the verical (y,z) plane [urad]
         
          T   = rp220->T;     // time of arrival of the particle in the detector [s] 
          side  = rp220->side;  // -1 or 1
        }
      //running on FP420 branch is identical
      TIter  itFp420((TCollection*)FP420);
      TRootRomanPotHits *fp420;
      itFp420.Reset();
      while( (fp420 = (TRootRomanPotHits*) itFp420.Next()) ){}
      
      //*********************************************
      //Example how to run on the trigger information
      //*********************************************
      
      TRootTrigger *trig;
      int NumTrigBit = TRIGGER->GetEntries();
      //get the global response of the trigger
      
      bool GlobalResponse=false;
      if(NumTrigBit!=0)GlobalResponse=true;
      //cout<<"GlobalResponse "<<GlobalResponse<<endl;
      for(int i=0; i < NumTrigBit-1; i++){
        trig = (TRootTrigger*)TRIGGER->At(i);
        cout<<"The event has been accepted by the trigger number: "<<trig->Accepted<<endl;
      }
      
     //********************************
     //Example of a very small analysis
     //********************************

     TLorentzVector Lept[2];

     if(NumTrigBit==0)continue;    //event not accepted by the trigger
     cut_trig++;//event accepted
        
     TSimpleArray<TRootElectron> el=SubArrayEl(ELEC,PT_ELEC);//the central isolated electrons, pt > PT_ELEC GeV 
     TSimpleArray<TRootMuon> mu=SubArrayMu(MUON,PT_MUON);//the central isolated electrons, pt > PT_MUON GeV 

     Int_t numElec=el.GetEntries();

     if(el.GetEntries()+mu.GetEntries()!=2)continue;//Exactly 2 isolated leptons are needed
     cut_1++;//event accepted
     for(Int_t i=0;i < numElec; i++)Lept[i].SetPxPyPzE(el[i]->Px,el[i]->Py,el[i]->Pz,el[i]->E);
     for(Int_t k = numElec; k < (numElec+mu.GetEntries()); k++)Lept[k].SetPxPyPzE(mu[k-numElec]->Px,mu[k-numElec]->Py,mu[k-numElec]->Pz,mu[k-numElec]->E);
    cout<<"normalement il y a quelque chose... "<<endl;
     //Example how to white a branch in the output file
     inv_mass=(TRootInvm*) INVMASS->NewEntry();
     inv_mass->M=(Lept[0]+Lept[1]).M();
 
     if((Lept[0]+Lept[1]).M() > INV_MASS_LL )continue;// the invariant mass should be < INV_MASS_LL
     cut_2++;//event accepted

     treeWriter->Fill(); 
   }
 treeWriter->Write();

}

void Analysis_Ex::WriteOutput(string LogName)
{
  ofstream f_out(LogName.c_str(),ofstream::app);

  f_out<<"*******************************************************************"<<endl;
  f_out << left  << setw(20) << "Numer of Events "<<""
        << right << setw(15) << total <<"\n";  
  f_out << left  << setw(17) << " Accepted by the trigger "<<""
        << right << setw(20) << cut_trig <<"\n";
  f_out << left  << setw(17) <<" 2 leptons "<< ""
        << right << setw(20) << cut_1 << ""        
        << right << setw(15) << cut_1/total << "\n";        
  f_out << left  << setw(17) <<" Invariant mass "<< ""
        << right << setw(20) << cut_2 << ""
        << right << setw(15) << cut_2/total << "\n";
  f_out<<"*******************************************************************"<<endl;
  f_out<<"                                                                   "<<endl;

}

TSimpleArray<TRootElectron> Analysis_Ex::SubArrayEl(const TClonesArray *ELEC,float pt)
{
    TIter itElec((TCollection*)ELEC);
    TRootElectron *elec;
    itElec.Reset();
    TSimpleArray<TRootElectron> array;
      while( (elec = (TRootElectron*) itElec.Next()) )
        {
         if(elec->PT<pt)continue;
         array.Add(elec);
        }
    return array;
}

TSimpleArray<TRootMuon> Analysis_Ex::SubArrayMu(const TClonesArray *MUON,float pt)
{
   TIter itMuon((TCollection*)MUON);
   TRootMuon *muon;
   itMuon.Reset();
   TSimpleArray<TRootMuon> array;
   while( (muon = (TRootMuon*) itMuon.Next()) )
     {
        if(muon->PT<pt)continue;
        array.Add(muon);
     }
   return array;
}