NtupleWriter14/src/JetTrackSelectorTool.cxx

//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// 17.07.2007, AUTHOR: OLIVER KORTNER
// Modified: 29.07.2007 by O. Kortner, options for a special treatment of
//                                     multivertex events.
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
//:: IMPLEMENTATION OF METHODS DEFINED IN THE CLASS JetTrackSelectorTool ::
//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

//::::::::::::::::::
//:: HEADER FILES ::
//::::::::::::::::::

// standard C++ //
#include <iostream>
#include <fstream>

// Gaudi //
#include "GaudiKernel/Service.h"
#include "GaudiKernel/MsgStream.h"

// track particles //
#include "Particle/TrackParticleContainer.h"
#include "egammaEvent/ElectronContainer.h"
#include "egammaEvent/egammaParamDefs.h"
#include "muonEvent/MuonContainer.h"

// jet reconstruction //
#include "JetUtils/JetCollectionHelper.h"
#include "JetEvent/Jet.h"
#include "JetEvent/JetCollection.h"
#include "JetRec/JetAlgToolBase.h"

//#include "JetRecTools/JetTrackSelectorTool.h"
#include "NtupleWriter14/JetTrackSelectorTool.h"
//#include "TrkEventPrimitives/FitQuality.h"

//::::::::::::::::::::::::
//:: NAMESPACE SETTINGS ::
//::::::::::::::::::::::::

using namespace std;
using namespace Rec;

//*****************************************************************************

//:::::::::::::::::
//:: CONSTRUCTOR ::
//:::::::::::::::::

JetTrackSelectorTool::JetTrackSelectorTool(const std::string & name,
        const std::string & type,
        const IInterface * parent) : JetAlgToolBase(name, type, parent) {

// JOB OPTIONS //

        m_pt_min = 0.0; // no cut on the transverse momentum of a track by
                        // default
        declareProperty("pTMin", m_pt_min);

        m_output_collection_name = string("TrackParticlesForJets");
        declareProperty("outputCollectionName", m_output_collection_name);

        m_track_particle_container_name = string("TrackParticleCandidate");
        declareProperty("TrackParticleContainer",
                                        m_track_particle_container_name);

        m_ignore_multivertices = false;
        declareProperty("ignoreMultivertices", m_ignore_multivertices);

        m_primary_vertex_container = string("VxPrimaryCandidate");
        declareProperty("vertexContainer", m_primary_vertex_container);

        m_vert_coll_svc_name = string("MPIHiggsAnalysis::VertexCollectionSvc");
        declareProperty("VertexCollectionSvc", m_vert_coll_svc_name);

        m_electron_container = string("ElectronCollection");
        declareProperty("electronContainer", m_electron_container);

        m_muon_container = string("StacoMuonCollection");
        declareProperty("muonContainer", m_muon_container);

        m_remove_isol_leptons = false;
        declareProperty("removeIsolLeptons", m_remove_isol_leptons);

        m_electron_cone   = 0.;
        m_electron_et_cut = 0.;
        m_muon_cone       = 0.;
        m_muon_et_cut     = 0.;
        m_inner_cone      = 0.01;

        declareProperty("electronCone",  m_electron_cone);
        declareProperty("electronEtCut", m_electron_et_cut);
        declareProperty("muonCone",      m_muon_cone);
        declareProperty("muonEtCut",     m_muon_et_cut);
        declareProperty("innerCone",     m_inner_cone);
        

// RESET POINTERS //

        m_storeGate = 0;
        m_log = 0;
        m_vert_coll = 0;

}

//*****************************************************************************

//::::::::::::::::
//:: DESTRUCTOR ::
//::::::::::::::::

JetTrackSelectorTool::~JetTrackSelectorTool(void) {}

//*****************************************************************************

//:::::::::::::::::::::::
//:: METHOD initialize ::
//:::::::::::::::::::::::

StatusCode JetTrackSelectorTool::initialize(void) {

// VARIABLES //

        m_log = new MsgStream(msgSvc(), name()); // message service

// MESSAGE //

        *m_log << MSG::INFO << "Initializing tool..." << endreq;

// SET POINTERS //

// store gate //
        StatusCode sc(service("StoreGateSvc", m_storeGate));
        if (!sc.isSuccess()) {
                *m_log << MSG::FATAL << "Cannot retrieve StoreGateSvc!"
                        << endreq;
                return sc;
        }

// VertexCollectionSvc //
        sc = service(m_vert_coll_svc_name, m_vert_coll);
        if (!sc.isSuccess()) {
                *m_log << MSG::FATAL << "Cannot retrieve "
                        << m_vert_coll_svc_name << "!" << endreq;
                return sc;
        }

        if(m_remove_isol_leptons){
            *m_log << MSG::INFO 
                   << "Ignore tracks wich overlap with an isolated lepton" 
                   << endreq
                   << "Isolation: " << endreq
                   << "  electron: etcone(" <<  m_electron_cone << ") < " 
                   <<  m_electron_et_cut << endreq
                   << "  muon:     etcone(" <<  m_muon_cone << ") < " 
                   <<  m_muon_et_cut << endreq;
        }

        return StatusCode::SUCCESS;

}

//*****************************************************************************

//::::::::::::::::::::
//:: METHOD execute ::
//::::::::::::::::::::

StatusCode JetTrackSelectorTool::execute(
                                JetAlgToolBase::jetcollection_t * theJets) {

// VARIABLES //

    // static std::ofstream out_e("et_incone_trkjets_e.txt");
    //static std::ofstream out_mu("et_incone_trkjets_mu.txt");

        const TrackParticleContainer *track_particle_container(0);
                                        // pointer to the track particles

        const ElectronContainer *electron_container(0); 
                                        // container with reconstructed electrons

        const Analysis::MuonContainer *muon_container(0); 
                                        // container with reconstructed muons

        jetcollection_t* writeCollection(0); // output jet collection
        Hep3Vector primary_vertex(0.0, 0.0, 0.0); // primary vertex of highest
                                                  // momentum


// GET THE PRIMARY VERTEX OF HIGHEST TRANSVERSE MOMENTUM, IF REQUESTED //

        if (!m_ignore_multivertices) {
                primary_vertex = highest_pt_vertex();
        }


// GET THE POINTER TO THE PARTICLE CONTAINER //

        StatusCode sc(m_storeGate->retrieve(track_particle_container,
                                        m_track_particle_container_name));
        if (sc.isFailure() || track_particle_container==0) {
                *m_log << MSG::WARNING
                        << "No track particle container found in store gate!"
                        << endreq;
                return sc;
        }
        *m_log << MSG::DEBUG
                << "Track particle container retrieved successfully!"
                << endreq;


        sc = m_storeGate->retrieve(electron_container, m_electron_container);
        if(sc.isFailure() || electron_container==0) {
            *m_log << MSG::WARNING
                   << "No electron container found in the store gate!"
                   << endreq;
            return false;
        }
        *m_log << MSG::DEBUG << "Electron container retrieved successfully!"
               << endreq;


        sc = m_storeGate->retrieve(muon_container, m_muon_container);
        if(sc.isFailure() || muon_container==0) {
            *m_log << MSG::WARNING
                   << "No muon container found in the store gate!"
                   << endreq;
            return false;
        }
        *m_log << MSG::DEBUG << "Muon container retrieved successfully!"
               << endreq;

    
// CREATE A NEW JET COLLECTION (TO BE USED AS THE INPUT OF THE JET RECO.) //

        writeCollection = new jetcollection_t();
                               
        // build vectors containing the isolated leptons                         //
        // isolation is done by cutting on the et sum of the tracks in the cone  //

        //ElectronContainer *isol_electron_container = new ElectronContainer(); 
        vector<const Analysis::Electron*>       isol_electron_vector(0);
        vector<const Analysis::Muon*>           isol_muon_vector(0);

        if ( m_remove_isol_leptons ){

            //electrons
            for (unsigned int iE=0; iE<electron_container->size(); iE++) {

                if ((*electron_container)[iE]
                    ->trackParticle()->reconstructedVertex() == NULL){
                    continue;
                }

                if (((*electron_container)[iE]->isem(egammaPID::ElectronLoose))!=0 ){
                  continue;
                } 

                //if (((*electron_container)[iE]->isem() & 0x7)!=0 ){
                //continue;
                //} 
 
                // if ((*electron_container)[iE]->isEM()!=0 ){
//                   continue;
//                 } 


                if ((*electron_container)[iE]->author()!=
                    egammaParameters::AuthorElectron &&
                    (*electron_container)[iE]->author()!=
                    egammaParameters::AuthorSofte) {
                    continue;
                }

                
                if ((*electron_container)[iE]->trackParticle()==0) {  
                    continue;
                }

                if (!m_ignore_multivertices){
                
                    if ((*electron_container)[iE]
                        ->trackParticle()->reconstructedVertex()->recVertex().position() != primary_vertex){
                    
                        continue; 
                    }
                }

                HepLorentzVector p((*electron_container)[iE]->hlv());
 
                double etcone = get_energy_in_cone( m_electron_cone, p,
                                                    (*electron_container)[iE]->trackParticle()
                                                    ->reconstructedVertex()->recVertex().position(),
                                                    track_particle_container
                                                    );
                //out_e << etcone << endl;
                //Electron *electron = const_cast<Electron*>((*electron_container)[iE]);
            
                if( etcone < m_electron_et_cut ){
                    //isol_electron_vector.push_back(electron);
                    isol_electron_vector.push_back((*electron_container)[iE]);
                }
                
            }

            //muons
            for (unsigned int iM=0; iM<muon_container->size(); iM++) {

                // if (!(*muon_container)[iM]->hasCombinedMuonTrackParticle() || !(*muon_container)[iM]->isHighPt() 
//                     || !(*muon_container)[iM]->bestMatch()) {
//                     continue;
//                 }

                if( (*muon_container)[iM]->track()->reconstructedVertex()==NULL ){
                    continue;
                } 

                if (!m_ignore_multivertices){
                
                  //   if( (*muon_container)[iM]->track()->reconstructedVertex()==NULL ){
//                         continue;
//                     } 
                    
                    if ((*muon_container)[iM]
                        ->track()->reconstructedVertex()->recVertex().position() != primary_vertex){
                    
                        continue; 
                    }
                }

                HepLorentzVector p((*muon_container)[iM]->hlv());
 
                double etcone = get_energy_in_cone( m_muon_cone, p,
                                                    (*muon_container)[iM]->track()
                                                    ->reconstructedVertex()->recVertex().position(),
                                                    track_particle_container
                                                    );
                
                //out_mu << etcone << endl;
                //Analysis::Muon *muon = const_cast<Analysis::Muon*>((*muon_container)[iM]);
            
                if( etcone < m_muon_et_cut ){
                    //isol_muon_vector.push_back(muon);
                    isol_muon_vector.push_back((*muon_container)[iM]);
                }
                
            }

        }

        // loop over track particles and add them to the jet collection //
        // if they pass all the cuts                                       //
       
        //static std::ofstream out("eta_trkjets.txt");
        for (unsigned int k=0; k<track_particle_container->size(); k++) {
            
            if((*track_particle_container)[k]->reconstructedVertex() == NULL){
                    continue;
            }

            HepLorentzVector p((*track_particle_container)[k]->hlv());
            //const Trk::FitQuality *fq = (*track_particle_container)[k]->fitQuality();
            
            //out << p.eta() << "\t\t" <<  p.et() << endl;
                //<< fq->chiSquared() << "\t\t" << fq->numberDoF()  << endl;
            
            
            //check pt
            if (p.perp()<m_pt_min) {
                continue;
            }
              
            //check vertex
            if (!m_ignore_multivertices){
                
                if ((*track_particle_container)[k]
                    ->reconstructedVertex()->recVertex().position() != primary_vertex){
                
                    continue;
                }
            } 

           
            //check overlap with isol leptons
            if ( m_remove_isol_leptons ){
                
                bool overlap(false);
                
                //electrons
                for (unsigned int iE=0; iE<isol_electron_vector.size(); iE++) {
                 
                    if( (*track_particle_container)[k] 
                        == isol_electron_vector[iE]->trackParticle() ){
                       
                        overlap = true;
                        break;
                    }
                } 
                
                if(overlap){
                    continue;
                }

                overlap = false;
                
                //muons
                for (unsigned int iM=0; iM<isol_muon_vector.size(); iM++) {
                 
                    if( (*track_particle_container)[k] 
                        == isol_muon_vector[iM]->inDetTrackParticle() ){
                        
                        overlap = true;
                        break;
                    }
                } 

                if(overlap){
                    continue;
                }

            } //m_remove_isol_leptons
            
            
            writeCollection->push_back(new Jet(p));

        } //for track_particle_container 


// record the collection in the store gate //
        if (!JetCollectionHelper::record_jets(m_storeGate, writeCollection,
                                                m_output_collection_name)) {
                *m_log << MSG::ERROR
                        << "Cannot convert initial list to JetCollection!"
                        << endreq;
                return StatusCode::FAILURE;
        }

// free memory //
        JetCollectionHelper::remove_object(writeCollection,
                                                writeCollection->begin(),
                                                writeCollection->end());
        delete writeCollection;

// retrieve collection from StoreGate to satisfy Jet constructor //
        const INavigable4MomentumCollection *readCollection;
        sc = m_storeGate->retrieve(readCollection, m_output_collection_name);
        if (sc.isFailure()) {
                *m_log << MSG::ERROR
                        << "Cannot read back initial jet collection!"
                        << endreq;
                return sc;
        }

// build intermediate jet list //
        for (Jet::index_type theIndex=0; theIndex < readCollection->size();
                                                                theIndex++) {
                theJets->push_back(new Jet(readCollection, theIndex));
        }


        return StatusCode::SUCCESS;

}

//*****************************************************************************

//::::::::::::::::::::::::::::::
//:: METHOD highest_pt_vertex ::
//::::::::::::::::::::::::::::::

Hep3Vector JetTrackSelectorTool::highest_pt_vertex(void) const {

        return m_vert_coll->getVertices()[
                                m_vert_coll->getVertices().size()-1].position();

// ///////////////
// // VARIABLES //
// ///////////////
// 
//      const Rec::TrackParticleContainer *track_particle_container(0);
//                                      // pointer to the track particles
//      const VxContainer *vertex_container(0); // pointer to the vertices
//      double max_Et; // maximum transerve vertex momentum found so far
//      Hep3Vector primary_vertex(0.0, 0.0, 0.0); // auxiliary primary vertex
//      vector<MPIHiggsAnalysis::Vertex> vertices(m_vert_coll->getVertices());
//                                                              // vertices
//      static ofstream out("debug_jets.txt");
// 
// ///////////////////////////////////////////
// // GET A POINTER TO THE VERTEX CONTAINER //
// ///////////////////////////////////////////
// 
//      StatusCode sc(m_storeGate->retrieve(vertex_container,
//                                              m_primary_vertex_container));
//      if (sc.isFailure() || vertex_container==0) {
//              *m_log << MSG::ERROR
//                      << "No vertex container "
//                      << m_primary_vertex_container
//                      << " found in the store gate!"
//                      << endreq;
//              return primary_vertex;
//      }
// 
//      *m_log << MSG::DEBUG
//              << "Vertex container "
//              << m_primary_vertex_container
//              << " retrieved successfully!"
//              << endreq;
// 
//      if (vertex_container->size()==0) {
//              return primary_vertex;
//      }
// 
// /////////////////////////////////////////////////////////////////////
// // CALCULATE SUM OF TRANSVERSE ENERGIES EMERGING FROM THE VERTICES //
// /////////////////////////////////////////////////////////////////////
// 
// // loop over the vertices //
//      sc = (m_storeGate->retrieve(track_particle_container,
//                                      m_track_particle_container_name));
//      if (sc.isFailure() || track_particle_container==0) {
//              *m_log << MSG::WARNING
//                      << "No track particle container found in store gate!"
//                      << endreq;
//              return primary_vertex;
//      }
//      *m_log << MSG::DEBUG
//              << "Track particle container retrieved successfully!"
//              << endreq;
// 
//      vector<double> Et_vert(vertex_container->size(), 0.0); // Et(vertex)
//      for (unsigned int k=0; k<track_particle_container->size(); k++) {
//              for (unsigned int l=0; l<vertex_container->size(); l++) {
//                      Trk::RecVertex vert((*vertex_container)[l]->recVertex());
//                      if ((*track_particle_container
//                              )[k]->reconstructedVertex()->position()==
//                              vert.position()) {
//                              Et_vert[l] = Et_vert[l]+
//                                      (*track_particle_container)[k]->hlv(
//                                                              ).et();
//                      }
//              }
//      }
// 
// /////////////////////////////////////////////////////////////////
// // GET THE PRIMARY VERTEX WITH THE HIGHEST TRANSVERSE MOMENTUM //
// /////////////////////////////////////////////////////////////////
//      
//      if (Et_vert.size()>0) {
//              max_Et = Et_vert[0];
//      } else {
//              max_Et = 0.0;
//      }
//      primary_vertex = ((*vertex_container)[0]->recVertex()).position();
// 
// // loop over the vertices //
//      for (unsigned int k=1; k<vertex_container->size(); k++) {
//              if (Et_vert[k]>max_Et) {
//                      max_Et = Et_vert[k];
//                      primary_vertex = ((*vertex_container)[k]->recVertex()
//                                                              ).position();
//              }
//      }
// 
//      out << vertex_container->size()<< "\t"
//              << m_vert_coll->getVertices().size() << "\t"
//              << max_Et << "\t"
//              << m_vert_coll->getVertices()[
//                              m_vert_coll->getVertices().size()-1].Et()
//              << "\t"
//              << primary_vertex << "\t"
//              << m_vert_coll->getVertices()[
//                              m_vert_coll->getVertices().size()-1].position()
//              << endl;
// 
//      return primary_vertex;

}

//*****************************************************************************

//:::::::::::::::::::::::::::::::
//:: METHOD get_energy_in_cone ::
//:::::::::::::::::::::::::::::::

double 
JetTrackSelectorTool::get_energy_in_cone(const double & cone,
                                         const HepLorentzVector & p, 
                                         const Hep3Vector & vertex,
                                         const TrackParticleContainer *track_particle_container) {
//-------------
// VARIABLES --
//-------------

    double ET(0.);   // transverse energy in the cone

//-------------------------------------------------------------
// LOOP OVER THE TRACKS AND CALCULATE THE ENERGY IN THE CONE --
//-------------------------------------------------------------

    for (unsigned int k=0; k<track_particle_container->size(); k++) {

        if (!m_ignore_multivertices){
            
            if ( (*track_particle_container)[k]
                 ->reconstructedVertex()->recVertex().position()!=vertex) {
                
                continue; 
            }
        }
    
        if ( p.deltaR((*track_particle_container)[k]->hlv()) < cone && 
             p.deltaR((*track_particle_container)[k]->hlv()) > m_inner_cone) {
    
            ET = ET+(*track_particle_container)[k]->hlv().et();
        }
        
    }
    

    return ET;
    
}

---