HiggsAnalysis/HiggsAnalysis-00-01-00/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 "ElectronPhotonIDEvent/ElectronContainer.h"
#include "MuonIDEvent/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 "HiggsAnalysis/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;

        m_cut_start = 0;
        m_cut_p = 0;
        m_cut_vertex = 0;
        m_cut_nhits = 0;
        m_cut_d0 = 0;
        m_cut_z0 = 0;


}

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

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

JetTrackSelectorTool::~JetTrackSelectorTool(void) {

    cout << "m_cut_start: " << m_cut_start << endl;
    cout << "m_cut_p: " << m_cut_p << endl;
    cout << "m_cut_vertex: " <<m_cut_vertex << endl;
    cout << "m_cut_nhits: " << m_cut_nhits<< endl;
    cout << "m_cut_d0: " << m_cut_d0<< endl;
    cout << "m_cut_z0: " << m_cut_z0<< endl;
    
}

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

//:::::::::::::::::::::::
//:: 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_dR_trkjets_e.txt");
    //static std::ofstream out_mu("et_dR_trkjets_mu.txt");

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

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

        const 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 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]
                    ->track()->reconstructedVertex() == NULL){
                    continue;
                }

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



                if((*electron_container)[iE]->author()!=
                   ElectronParameters::AuthorEgamma &&
                   (*electron_container)[iE]->author()!=
                   ElectronParameters::AuthorEgammaSofte){
                    continue;
                }
                
                if (!(*electron_container)[iE]->hasTrack()) {  
                    continue;
                }

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

                HepLorentzVector p((*electron_container)[iE]->hlv());
 
                double etcone = get_energy_in_cone( m_electron_cone, p,
                                                    (*electron_container)[iE]->track()
                                                    ->reconstructedVertex()->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()->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()->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;
            const Trk::TrackSummary* trksum = (*track_particle_container)[k]->trackSummary();
            int trk_num_SCTHit  = trksum->get(Trk::numberOfSCTHits);
            int trk_num_PHit    = trksum->get(Trk::numberOfPixelHits);   

            //const Trk::MeasuredPerigee* perigee = (*track_particle_container)[k]->measuredPerigee();
            //double d0 = perigee->parameters()[Trk::d0]; 
            //double z0 = perigee->parameters()[Trk::z0];
            //double sinTh = perigee->sinTheta();
            
            m_cut_start++;
            
            //check pt
            if (p.perp()<m_pt_min) {
                continue;
            }
            m_cut_p++;
 
            //check vertex
            if (!m_ignore_multivertices){
                
                if ((*track_particle_container)[k]
                    ->reconstructedVertex()->position() != primary_vertex){
                
                    continue;
                }
            } 
            m_cut_vertex++;

            if((trk_num_SCTHit+trk_num_PHit)<7){
                continue;
            }
            m_cut_nhits++;

            //if(fabs(d0)>2){
            //  continue;
            //}
            //m_cut_d0++;
            
            //if(fabs(z0*sinTh)>10){
            //  continue;
            //}
            //m_cut_z0++;
            
            //cout << (*track_particle_container)[k]->reconstructedVertex()->position() << endl;
            //cout << trk_num_SCTHit << ", " << trk_num_PHit << ", " 
            //   << fabs(d0) << ", " << fabs(z0*sinTh) << endl;
                       
            //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]->track() ){
                       
                        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 {

    
    //for (unsigned int k=0; k<m_vert_coll->getVertices().size(); k++){
    //  cout << m_vert_coll->getVertices()[k].Et() << endl;
    //}


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


}

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

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

---