JetAnalyzers/src/JetIdSelector.cc

0001 /** \class CaloJetIdSelector
0002  *
0003  * Select a subset of the CaloJet collection based on the quality stored in a ValueMap
0004  *
0005  * \author: Kalanand Mishra, Fermilab
0006  *
0007  * usage:
0008  *
0009  * from RecoJets.JetProducers.ak5JetID_cfi import *
0010  * module CaloJetsLooseId = cms.EDProducer("CaloJetIdSelector",
0011  *    src     = cms.InputTag( "ak5CaloJets" ),
0012  *    idLevel = cms.string("LOOSE"),
0013  *    jetIDMap = cms.untracked.InputTag("ak5JetID")
0014  *               ### must provide jet ID value map for CaloJets
0015  * )
0016  *
0017  *
0018  * module PFJetsLooseId = cms.EDProducer("PFJetIdSelector",
0019  *    src     = cms.InputTag( "ak5PFJets" ),
0020  *    idLevel = cms.string("LOOSE"),
0021  * )
0022  *
0023  *
0024  *
0025  */
0026
0027 #include "PhysicsTools/SelectorUtils/interface/JetIDSelectionFunctor.h"
0028 #include "FWCore/Framework/interface/Frameworkfwd.h"
0029 #include "FWCore/Framework/interface/one/EDProducer.h"
0030 #include "FWCore/Framework/interface/Event.h"
0031 #include "FWCore/Framework/interface/EventSetup.h"
0032 #include "FWCore/Framework/interface/MakerMacros.h"
0033
0034 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0035 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0036 #include "FWCore/ServiceRegistry/interface/Service.h"
0037
0038 #include "DataFormats/Common/interface/Handle.h"
0039 #include "DataFormats/JetReco/interface/JetCollection.h"
0040 #include "DataFormats/JetReco/interface/CaloJet.h"
0041 #include "DataFormats/JetReco/interface/CaloJetCollection.h"
0042 #include "DataFormats/JetReco/interface/GenJet.h"
0043 #include "DataFormats/JetReco/interface/GenJetCollection.h"
0044 #include "DataFormats/JetReco/interface/PFJet.h"
0045 #include "DataFormats/JetReco/interface/PFJetCollection.h"
0046 #include "DataFormats/JetReco/interface/JPTJet.h"
0047 #include "DataFormats/JetReco/interface/JPTJetCollection.h"
0048 #include "DataFormats/JetReco/interface/GenJet.h"
0049 #include "DataFormats/JetReco/interface/GenJetCollection.h"
0050
0051 #include <memory>
0052 #include <vector>
0053 #include <sstream>
0054
0055 ////////////////////////////////////////////////////////////////////////////////
0056 // class definition
0057 ////////////////////////////////////////////////////////////////////////////////
0058 template <typename T>
0059 class JetIdSelector : public edm::one::EDProducer<> {
0060 public:
0061   typedef std::vector<T> JetCollection;
0062   // construction/destruction
0063   explicit JetIdSelector(const edm::ParameterSet& iConfig);
0064   ~JetIdSelector() override;
0065
0066   // member functions
0067   void produce(edm::Event& iEvent, const edm::EventSetup& iSetup) override;
0068   void endJob() override;
0069
0070 private:
0071   // member data
0072   edm::InputTag src_;
0073   std::string qualityStr;
0074   edm::InputTag jetIDMap_;
0075   std::string moduleLabel_;
0076
0077   unsigned int nJetsTot_;
0078   unsigned int nJetsPassed_;
0079   JetIDSelectionFunctor* jetIDFunctor;
0080   bool use_pfloose;
0081   bool use_pfmedium;
0082   bool use_pftight;
0083 };
0084
0085 using namespace std;
0086
0087 ////////////////////////////////////////////////////////////////////////////////
0088 // construction/destruction
0089 ////////////////////////////////////////////////////////////////////////////////
0090 //______________________________________________________________________________
0091 template <typename T>
0092 JetIdSelector<T>::JetIdSelector(const edm::ParameterSet& iConfig)
0093     : src_(iConfig.getParameter<edm::InputTag>("src")),
0094       qualityStr(iConfig.getParameter<string>("idLevel")),
0095       jetIDMap_(iConfig.getUntrackedParameter<edm::InputTag>("jetIDMap", edm::InputTag("ak5JetID"))),
0096       moduleLabel_(iConfig.getParameter<string>("@module_label")),
0097       nJetsTot_(0),
0098       nJetsPassed_(0) {
0099   produces<JetCollection>();
0100
0101   use_pfloose = false;
0102   use_pfmedium = false;
0103   use_pftight = false;
0104
0105   if (qualityStr == "MINIMAL") {
0106     jetIDFunctor = new JetIDSelectionFunctor(JetIDSelectionFunctor::PURE09, JetIDSelectionFunctor::MINIMAL);
0107     use_pfloose = true;
0108   } else if (qualityStr == "LOOSE_AOD") {
0109     jetIDFunctor = new JetIDSelectionFunctor(JetIDSelectionFunctor::PURE09, JetIDSelectionFunctor::LOOSE_AOD);
0110     use_pfloose = true;
0111   } else if (qualityStr == "LOOSE") {
0112     jetIDFunctor = new JetIDSelectionFunctor(JetIDSelectionFunctor::PURE09, JetIDSelectionFunctor::LOOSE);
0113     use_pfloose = true;
0114   } else if (qualityStr == "MEDIUM") {
0115     jetIDFunctor = new JetIDSelectionFunctor(JetIDSelectionFunctor::PURE09, JetIDSelectionFunctor::LOOSE);
0116     // There is no medium quality for CaloJet !!
0117     use_pfmedium = true;
0118   } else if (qualityStr == "TIGHT") {
0119     jetIDFunctor = new JetIDSelectionFunctor(JetIDSelectionFunctor::PURE09, JetIDSelectionFunctor::TIGHT);
0120     use_pftight = true;
0121   } else
0122     throw cms::Exception("InvalidInput") << "Expect quality to be one of MINIMAL, LOOSE_AOD, LOOSE, MEDIUM, TIGHT"
0123                                          << std::endl;
0124 }
0125
0126 //______________________________________________________________________________
0127 template <typename T>
0128 JetIdSelector<T>::~JetIdSelector() {
0129   if (jetIDFunctor)
0130     delete jetIDFunctor;
0131 }
0132
0133 ////////////////////////////////////////////////////////////////////////////////
0134 // implementation of member functions
0135 ////////////////////////////////////////////////////////////////////////////////
0136
0137 //______________________________________________________________________________
0138 template <typename T>
0139 void JetIdSelector<T>::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0140   auto selectedJets = std::make_unique<JetCollection>();
0141   edm::Handle<reco::JetView> hjets;  // uncorrected jets!
0142   iEvent.getByLabel(src_, hjets);
0143   auto const& jets = *hjets;
0144
0145   // handle to the jet ID variables, filled in the loop
0146   edm::Handle<reco::JetIDValueMap> hJetIDMap;
0147   unsigned int idx = 0;
0148   bool passed = false;
0149
0150   for (edm::View<reco::Jet>::const_iterator ibegin = jets.begin(), iend = jets.end(), iJet = ibegin; iJet != iend;
0151        ++iJet) {
0152     // initialize the boolean flag to false
0153     passed = false;
0154
0155     //calculate the Calo jetID
0156     auto const& jet = *iJet;
0157     const std::type_info& type = typeid(jet);
0158     if (type == typeid(reco::CaloJet)) {
0159       const reco::CaloJet& calojet = static_cast<const reco::CaloJet&>(jet);
0160       edm::RefToBase<reco::Jet> jetRef = jets.refAt(idx);
0161       if (not hJetIDMap.isValid()) {
0162         iEvent.getByLabel(jetIDMap_, hJetIDMap);
0163       }
0164       reco::JetID const& jetId = (*hJetIDMap)[jetRef];
0165       passed = (*jetIDFunctor)(calojet, jetId);
0166     }
0167
0168     //calculate the PF jetID
0169     if (type == typeid(reco::PFJet)) {
0170       const reco::PFJet& pfjet = static_cast<const reco::PFJet&>(jet);
0171       bool passingLoose = false;
0172       bool passingMedium = false;
0173       bool passingTight = false;
0174       bool ThisIsClean = true;
0175       //apply following only if |eta|<2.4: CHF>0, CEMF<0.99, chargedMultiplicity>0
0176       if ((pfjet.chargedHadronEnergy() / pfjet.energy()) <= 0.0 && fabs(pfjet.eta()) < 2.4)
0177         ThisIsClean = false;
0178       if ((pfjet.chargedEmEnergy() / pfjet.energy()) >= 0.99 && fabs(pfjet.eta()) < 2.4)
0179         ThisIsClean = false;
0180       if (pfjet.chargedMultiplicity() <= 0 && fabs(pfjet.eta()) < 2.4)
0181         ThisIsClean = false;
0182
0183       // always require #Constituents > 1
0184       if (pfjet.nConstituents() <= 1)
0185         ThisIsClean = false;
0186
0187       if (ThisIsClean && (pfjet.neutralHadronEnergy() / pfjet.energy()) < 0.99 &&
0188           (pfjet.neutralEmEnergy() / pfjet.energy()) < 0.99)
0189         passingLoose = true;
0190
0191       if (ThisIsClean && (pfjet.neutralHadronEnergy() / pfjet.energy()) < 0.95 &&
0192           (pfjet.neutralEmEnergy() / pfjet.energy()) < 0.95)
0193         passingMedium = true;
0194
0195       if (ThisIsClean && (pfjet.neutralHadronEnergy() / pfjet.energy()) < 0.90 &&
0196           (pfjet.neutralEmEnergy() / pfjet.energy()) < 0.90)
0197         passingTight = true;
0198
0199       if (use_pfloose && passingLoose)
0200         passed = true;
0201       if (use_pfmedium && passingMedium)
0202         passed = true;
0203       if (use_pftight && passingTight)
0204         passed = true;
0205     }
0206
0207     if (type == typeid(reco::GenJet) || type == typeid(reco::JPTJet)) {
0208       edm::LogWarning("JetId") << "Criteria for jets other than CaloJets and PFJets are not yet implemented";
0209       passed = true;
0210     }  // close GenJet, JPT jet
0211
0212     const T& goodJet = static_cast<const T&>(jet);
0213     if (passed)
0214       selectedJets->push_back(goodJet);
0215
0216     idx++;
0217   }  // close jet iterator
0218
0219   nJetsTot_ += jets.size();
0220   nJetsPassed_ += selectedJets->size();
0221   iEvent.put(std::move(selectedJets));
0222 }
0223
0224 //______________________________________________________________________________
0225 template <typename T>
0226 void JetIdSelector<T>::endJob() {
0227   stringstream ss;
0228   ss << "nJetsTot=" << nJetsTot_ << " nJetsPassed=" << nJetsPassed_
0229      << " fJetsPassed=" << 100. * (nJetsPassed_ / (double)nJetsTot_) << "%\n";
0230   cout << "++++++++++++++++++++++++++++++++++++++++++++++++++"
0231        << "\n"
0232        << moduleLabel_ << "(JetIdSelector) SUMMARY:\n"
0233        << ss.str() << "++++++++++++++++++++++++++++++++++++++++++++++++++" << endl;
0234 }
0235
0236 ////////////////////////////////////////////////////////////////////////////////
0237 // plugin definition
0238 ////////////////////////////////////////////////////////////////////////////////
0239 typedef JetIdSelector<reco::CaloJet> CaloJetIdSelector;
0240 typedef JetIdSelector<reco::PFJet> PFJetIdSelector;
0241 typedef JetIdSelector<reco::JPTJet> JPTJetIdSelector;
0242 typedef JetIdSelector<reco::GenJet> GenJetIdSelector;
0243
0244 DEFINE_FWK_MODULE(CaloJetIdSelector);
0245 DEFINE_FWK_MODULE(PFJetIdSelector);
0246 DEFINE_FWK_MODULE(JPTJetIdSelector);
0247 DEFINE_FWK_MODULE(GenJetIdSelector);