/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_5_3_14/src/PhysicsTools/SelectorUtils/interface/PFJetIDSelectionFunctor.h

Go to the documentation of this file.

#ifndef PhysicsTools_SelectorUtils_interface_PFJetIDSelectionFunctor_h
#define PhysicsTools_SelectorUtils_interface_PFJetIDSelectionFunctor_h

#include "DataFormats/PatCandidates/interface/Jet.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include "PhysicsTools/SelectorUtils/interface/Selector.h"

#include <TMath.h>
class PFJetIDSelectionFunctor : public Selector<pat::Jet>  {

 public: // interface

  enum Version_t { FIRSTDATA, N_VERSIONS };
  enum Quality_t { LOOSE, TIGHT, N_QUALITY};

  PFJetIDSelectionFunctor() {}
  
 PFJetIDSelectionFunctor( edm::ParameterSet const & params ) 
 {
   std::string versionStr = params.getParameter<std::string>("version");
   std::string qualityStr = params.getParameter<std::string>("quality");

   if ( versionStr == "FIRSTDATA" ) 
     version_ = FIRSTDATA;
   else
     version_ = FIRSTDATA;  

   if      ( qualityStr == "LOOSE") quality_ = LOOSE;
   else if ( qualityStr == "TIGHT") quality_ = TIGHT;
   else quality_ = LOOSE;

    push_back("CHF" );
    push_back("NHF" );
    push_back("CEF" );
    push_back("NEF" );
    push_back("NCH" );
    push_back("nConstituents");


    // Set some default cuts for LOOSE, TIGHT
    if ( quality_ == LOOSE ) {
      set("CHF", 0.0);
      set("NHF", 0.99);
      set("CEF", 0.99);
      set("NEF", 0.99);
      set("NCH", 0);
      set("nConstituents", 1);
    } else if ( quality_ == TIGHT ) {
      set("CHF", 0.0);
      set("NHF", 0.9);
      set("CEF", 0.99);
      set("NEF", 0.9);
      set("NCH", 0);
      set("nConstituents", 1);      
    }


    // Now check the configuration to see if the user changed anything
    if ( params.exists("CHF") ) set("CHF", params.getParameter<double>("CHF") );
    if ( params.exists("NHF") ) set("NHF", params.getParameter<double>("NHF") );
    if ( params.exists("CEF") ) set("CEF", params.getParameter<double>("CEF") );
    if ( params.exists("NEF") ) set("NEF", params.getParameter<double>("NEF") );
    if ( params.exists("NCH") ) set("NCH", params.getParameter<int>   ("NCH") );
    if ( params.exists("nConstuents") ) set("nConstituents", params.getParameter<int> ("nConstituents") );

    if ( params.exists("cutsToIgnore") )
      setIgnoredCuts( params.getParameter<std::vector<std::string> >("cutsToIgnore") );


    indexNConstituents_ = index_type (&bits_, "nConstituents");
    indexNEF_ = index_type (&bits_, "NEF");
    indexNHF_ = index_type (&bits_, "NHF");
    indexCEF_ = index_type (&bits_, "CEF");
    indexCHF_ = index_type (&bits_, "CHF");
    indexNCH_ = index_type (&bits_, "NCH");

    retInternal_ = getBitTemplate();
    
  }


 PFJetIDSelectionFunctor( Version_t version,
                          Quality_t quality ) :
  version_(version), quality_(quality)
 {

    push_back("CHF" );
    push_back("NHF" );
    push_back("CEF" );
    push_back("NEF" );
    push_back("NCH" );
    push_back("nConstituents");


    // Set some default cuts for LOOSE, TIGHT
    if ( quality_ == LOOSE ) {
      set("CHF", 0.0);
      set("NHF", 0.99);
      set("CEF", 0.99);
      set("NEF", 0.99);
      set("NCH", 0);
      set("nConstituents", 1);
    } else if ( quality_ == TIGHT ) {
      set("CHF", 0.0);
      set("NHF", 0.9);
      set("CEF", 0.99);
      set("NEF", 0.9);
      set("NCH", 0);
      set("nConstituents", 1);      
    }


    indexNConstituents_ = index_type (&bits_, "nConstituents");
    indexNEF_ = index_type (&bits_, "NEF");
    indexNHF_ = index_type (&bits_, "NHF");
    indexCEF_ = index_type (&bits_, "CEF");
    indexCHF_ = index_type (&bits_, "CHF");
    indexNCH_ = index_type (&bits_, "NCH");

    retInternal_ = getBitTemplate();   
 }
                           

  // 
  // Accessor from PAT jets
  // 
  bool operator()( const pat::Jet & jet, pat::strbitset & ret )  
  {
    if ( version_ == FIRSTDATA ) {
      if ( jet.currentJECLevel() == "Uncorrected" || !jet.jecSetsAvailable() ) 
        return firstDataCuts( jet, ret );
      else 
        return firstDataCuts( jet.correctedJet("Uncorrected"), ret );
    }
    else {
      return false;
    }
  }
  using Selector<pat::Jet>::operator();

  // 
  // Accessor from *CORRECTED* 4-vector, EMF, and Jet ID. 
  // This can be used with reco quantities. 
  // 
  bool operator()( const reco::PFJet & jet, pat::strbitset & ret )  
  {
    if ( version_ == FIRSTDATA ) return firstDataCuts( jet, ret );
    else {
      return false;
    }
  }

  bool operator()( const reco::PFJet & jet )
  {
    retInternal_.set(false);
    operator()(jet, retInternal_);
    setIgnored(retInternal_);
    return (bool)retInternal_;
  }
  
  // 
  // cuts based on craft 08 analysis. 
  // 
  bool firstDataCuts( reco::Jet const & jet,
                      pat::strbitset & ret) 
  {    
    ret.set(false);

    // cache some variables
    double chf = 0.0;
    double nhf = 0.0;
    double cef = 0.0;
    double nef = 0.0;
    int    nch = 0;
    int    nconstituents = 0;

    // Have to do this because pat::Jet inherits from reco::Jet but not reco::PFJet
    reco::PFJet const * pfJet = dynamic_cast<reco::PFJet const *>(&jet);
    pat::Jet const * patJet = dynamic_cast<pat::Jet const *>(&jet);
    reco::BasicJet const * basicJet = dynamic_cast<reco::BasicJet const *>(&jet);

    if ( patJet != 0 ) {
      if ( patJet->isPFJet() ) {
        chf = patJet->chargedHadronEnergyFraction();
        nhf = ( patJet->neutralHadronEnergy() + patJet->HFHadronEnergy() ) / patJet->energy();
        cef = patJet->chargedEmEnergyFraction();
        nef = patJet->neutralEmEnergyFraction();
        nch = patJet->chargedMultiplicity();
        nconstituents = patJet->numberOfDaughters();
      } 
      // Handle the special case where this is a composed jet for
      // subjet analyses
      else if ( patJet->isBasicJet() ) {
        double e_chf = 0.0;
        double e_nhf = 0.0;
        double e_cef = 0.0;
        double e_nef = 0.0;
        nch = 0;
        nconstituents = 0;

        for ( reco::Jet::const_iterator ibegin = patJet->begin(),
                iend = patJet->end(), isub = ibegin;
              isub != iend; ++isub ) {
          reco::PFJet const * pfsub = dynamic_cast<reco::PFJet const *>( &*isub );
          e_chf += pfsub->chargedHadronEnergy();
          e_nhf += (pfsub->neutralHadronEnergy() + pfsub->HFHadronEnergy());
          e_cef += pfsub->chargedEmEnergy();
          e_nef += pfsub->neutralEmEnergy();
          nch += pfsub->chargedMultiplicity();
          nconstituents += pfsub->numberOfDaughters();
        }
        double e = patJet->energy();
        if ( e > 0.000001 ) {
          chf = e_chf / e;
          nhf = e_nhf / e;
          cef = e_cef / e;
          nef = e_nef / e;
        } else {
          chf = nhf = cef = nef = 0.0;
        }
      }
    } // end if pat jet
    else if ( pfJet != 0 ) {
      // CV: need to compute energy fractions in a way that works for corrected as well as for uncorrected PFJets
      double jetEnergyUncorrected = 
        pfJet->chargedHadronEnergy() 
       + pfJet->neutralHadronEnergy()
       + pfJet->photonEnergy()
       + pfJet->electronEnergy()
       + pfJet->muonEnergy()
       + pfJet->HFHadronEnergy()
       + pfJet->HFEMEnergy();
      if ( jetEnergyUncorrected > 0. ) {
        chf = pfJet->chargedHadronEnergy() / jetEnergyUncorrected;
        nhf = ( pfJet->neutralHadronEnergy() + pfJet->HFHadronEnergy() ) / jetEnergyUncorrected;
        cef = pfJet->chargedEmEnergy() / jetEnergyUncorrected;
        nef = pfJet->neutralEmEnergy() / jetEnergyUncorrected;
      }
      nch = pfJet->chargedMultiplicity();
      nconstituents = pfJet->numberOfDaughters();
    } // end if PF jet
    // Handle the special case where this is a composed jet for
    // subjet analyses
    else if ( basicJet != 0 ) {
      double e_chf = 0.0;
      double e_nhf = 0.0;
      double e_cef = 0.0;
      double e_nef = 0.0;
      nch = 0;
      nconstituents = 0;
      
      for ( reco::Jet::const_iterator ibegin = basicJet->begin(),
              iend = patJet->end(), isub = ibegin;
            isub != iend; ++isub ) {
        reco::PFJet const * pfsub = dynamic_cast<reco::PFJet const *>( &*isub );
        e_chf += pfsub->chargedHadronEnergy();
        e_nhf += (pfsub->neutralHadronEnergy() + pfsub->HFHadronEnergy());
        e_cef += pfsub->chargedEmEnergy();
        e_nef += pfsub->neutralEmEnergy();
        nch += pfsub->chargedMultiplicity();
        nconstituents += pfsub->numberOfDaughters();
      }
      double e = basicJet->energy();
      if ( e > 0.000001 ) {
        chf = e_chf / e;
        nhf = e_nhf / e;
        cef = e_cef / e;
        nef = e_nef / e;
      }
    } // end if basic jet


    // Cuts for all |eta|:
    if ( ignoreCut(indexNConstituents_) || nconstituents > cut(indexNConstituents_, int() ) ) passCut( ret, indexNConstituents_);
    if ( ignoreCut(indexNEF_)           || ( nef < cut(indexNEF_, double()) ) ) passCut( ret, indexNEF_);
    if ( ignoreCut(indexNHF_)           || ( nhf < cut(indexNHF_, double()) ) ) passCut( ret, indexNHF_);    
    // Cuts for |eta| < 2.4:
    if ( ignoreCut(indexCEF_)           || ( cef < cut(indexCEF_, double()) || std::abs(jet.eta()) > 2.4 ) ) passCut( ret, indexCEF_);
    if ( ignoreCut(indexCHF_)           || ( chf > cut(indexCHF_, double()) || std::abs(jet.eta()) > 2.4 ) ) passCut( ret, indexCHF_);
    if ( ignoreCut(indexNCH_)           || ( nch > cut(indexNCH_, int())    || std::abs(jet.eta()) > 2.4 ) ) passCut( ret, indexNCH_);    

    //std::cout << "<PFJetIDSelectionFunctor::firstDataCuts>:" << std::endl;
    //std::cout << " jet: Pt = " << jet.pt() << ", eta = " << jet.eta() << ", phi = " << jet.phi() << std::endl;
    //ret.print(std::cout);

    setIgnored( ret );
    return (bool)ret;
  }
  
 private: // member variables
  
  Version_t version_;
  Quality_t quality_;

  index_type indexNConstituents_;
  index_type indexNEF_;
  index_type indexNHF_;
  index_type indexCEF_;
  index_type indexCHF_;
  index_type indexNCH_;
  
};

#endif