Merge pull request cms-sw#23 from jhgoh/unifiedJetProducer

Unified jet producer

Configuration/python/ntuple_template_cff.py

@@ -38,7 +38,6 @@
 ntupleSequenceElEl = cms.Sequence(
     pileupWeight + pdfWeight
   + nEventsNtupleElEl
-  + jetUnc
   + goodMuons + goodElectrons * goodJets
   + jpsiToMuMu
   * ElEl
@@ -47,7 +46,6 @@
 ntupleSequenceMuMu = cms.Sequence(
     pileupWeight + pdfWeight
   + nEventsNtupleMuMu
-  + jetUnc
   + goodMuons + goodElectrons * goodJets
   + jpsiToMuMu
   * MuMu
@@ -56,7 +54,6 @@
 ntupleSequenceMuEl = cms.Sequence(
     pileupWeight + pdfWeight
   + nEventsNtupleMuEl
-  + jetUnc
   + goodMuons + goodElectrons * goodJets
   + jpsiToMuMu
   * MuEl
@@ -65,7 +62,6 @@
 ntupleSequenceMuJets = cms.Sequence(
     pileupWeight + pdfWeight
   + nEventsNtupleMuJets
-  + jetUnc
   + goodMuons + goodElectrons * goodJets
   + jpsiToMuMu
   * MuJets
@@ -74,7 +70,6 @@
 ntupleSequenceElJets = cms.Sequence(
     pileupWeight + pdfWeight
   + nEventsNtupleElJets
-  + jetUnc
   + goodMuons + goodElectrons * goodJets
   + jpsiToMuMu
   * ElJets

GenericNtuple/plugins/KGenericNtupleMaker.cc

@@ -370,7 +370,7 @@ void KGenericNtupleMaker::analyze(const edm::Event& event, const edm::EventSetup
   }
 
   // Do jets
-  edm::Handle<JetRefs> jetHandle;
+  edm::Handle<Jets> jetHandle;
   edm::Handle<pat::JetToValue> fJESUpHandle, fJESDnHandle;
   edm::Handle<pat::JetToValue> fJERHandle;
   edm::Handle<pat::JetToValue> fJERUpHandle, fJERDnHandle;
@@ -387,7 +387,7 @@ void KGenericNtupleMaker::analyze(const edm::Event& event, const edm::EventSetup
   double fJER = 1, fJERUp = 1, fJERDn = 1;
   for ( size_t i=0, n=jetHandle->size(); i<n; ++i )
   {
-    edm::Ref<Jets> jetRef = jetHandle->at(i);
+    edm::Ref<Jets> jetRef(jetHandle, i);
     const pat::Jet& jet = *jetRef;
     const double jetEta = jet.eta();
     if ( std::abs(jetEta) > 2.5 ) continue;

GenericNtuple/python/genericNtupleMaker_cfi.py

@@ -27,7 +27,7 @@
     jetMET = cms.PSet(
         jet = cms.InputTag("goodJets"),
         met = cms.InputTag("patMETsPFlow"),
-        unc = cms.InputTag("jetUnc"),
+        unc = cms.InputTag("goodJets"),
         minNumber = cms.uint32(0),
         leptonDeltaR = cms.double(0.5),
         bTagType = cms.string("combinedSecondaryVertexBJetTags"),

LeptonAnalysis/test/electronTnP/ntuple_MC_cfg.py

@@ -239,7 +239,7 @@
   + process.patSequenceComplete
   + process.hltHighLevel
   + process.leptonSelectionSequence
-  + process.jetUnc * process.goodJets
+  + process.goodJets
   * process.tpPairs * process.tp
 )
 

LeptonAnalysis/test/electronTnP/ntuple_RD_cfg.py

@@ -237,7 +237,7 @@
   + process.patSequenceComplete
   + process.hltHighLevel
   + process.leptonSelectionSequence
-  + process.jetUnc * process.goodJets
+  + process.goodJets
   * process.tpPairs * process.tp
 )
 

LeptonAnalysis/test/muonTnP/ntuple_RD_cfg.py

@@ -264,7 +264,7 @@
   + process.patSequenceComplete
   + process.hltHighLevel
   + process.muonSelectionSequence
-#  * process.goodElectrons + process.jetUnc * process.goodJets
+#  * process.goodElectrons + process.goodJets
 #  * process.tpPairs * process.muonEffs
 )
 

RecoSelectorTools/plugins/KJetUncProducer.cc → RecoSelectorTools/plugins/KJetMetUncSelector.cc

@@ -1,6 +1,6 @@
 #include "FWCore/Framework/interface/Frameworkfwd.h"
 #include "FWCore/Framework/interface/MakerMacros.h"
-#include "FWCore/Framework/interface/EDProducer.h"
+#include "FWCore/Framework/interface/EDFilter.h"
 
 #include "FWCore/Framework/interface/Event.h"
 #include "FWCore/ParameterSet/interface/ParameterSet.h"
@@ -13,6 +13,10 @@
 #include "CondFormats/JetMETObjects/interface/FactorizedJetCorrector.h"
 #include "CondFormats/JetMETObjects/interface/JetCorrectorParameters.h"
 #include "CondFormats/JetMETObjects/interface/JetCorrectionUncertainty.h"
+#include "PhysicsTools/SelectorUtils/interface/PFJetIDSelectionFunctor.h"
+//#include "CommonTools/UtilAlgos/interface/StringCutObjectSelector.h"
+#include "CommonTools/Utils/interface/PtComparator.h"
+#include "DataFormats/Math/interface/deltaR.h"
 
 #include "KrAFT/RecoSelectorTools/interface/Types.h"
 
@@ -21,39 +25,59 @@
 using namespace edm;
 using namespace std;
 
-class KJetUncProducer : public edm::EDProducer
+class KJetMetUncSelector : public edm::EDFilter
 {
 public:
-  KJetUncProducer(const edm::ParameterSet& pset);
-  ~KJetUncProducer() {};
+  KJetMetUncSelector(const edm::ParameterSet& pset);
+  ~KJetMetUncSelector() {};
+
   typedef std::vector<pat::Jet> Jets;
   typedef std::vector<pat::MET> METs;
 
 private:
   void beginJob() {};
-  void produce(edm::Event& event, const edm::EventSetup& eventSetup);
+  bool filter(edm::Event& event, const edm::EventSetup& eventSetup);
   void endJob() {};
 
-  edm::InputTag jetLabel_, metLabel_;
+  bool doClean_;
+  double overlapDeltaR_;
+  unsigned int minNumber_;
+  unsigned int maxNumber_;
+
+  edm::InputTag jetLabel_;
+  edm::InputTag metLabel_;
+  std::vector<edm::InputTag> overlapCandLabels_;
 
-  JetCorrectionUncertainty *jecUncCalculator_;
+  PFJetIDSelectionFunctor* isGoodJet_;
+  JetCorrectionUncertainty* jecUncCalculator_;
+  double minPt_, maxEta_;
 
   bool isMC_;
 
 private:
-  void getJERFactor(const double jetEta, double& cJER, double& cJERUp, double& cJERDn)
+  void loadJEC(const double jetPt, const double jetEta, double& fJECUp, double& fJECDn)
+  {
+    jecUncCalculator_->setJetPt(jetPt);
+    jecUncCalculator_->setJetEta(jetEta);
+    fJECUp = 1+jecUncCalculator_->getUncertainty(true);
+    jecUncCalculator_->setJetPt(jetPt);
+    jecUncCalculator_->setJetEta(jetEta);
+    fJECDn = 1-jecUncCalculator_->getUncertainty(false);
+  };
+
+  void loadJER(const double jetEta, double& cJER, double& cJERUp, double& cJERDn) const
   {
     if      ( jetEta < 0.5 ) { cJER = 1.052; cJERUp = 1.115; cJERDn = 0.990; }
     else if ( jetEta < 1.1 ) { cJER = 1.057; cJERUp = 1.114; cJERDn = 1.001; }
     else if ( jetEta < 1.7 ) { cJER = 1.096; cJERUp = 1.161; cJERDn = 1.032; }
     else if ( jetEta < 2.3 ) { cJER = 1.134; cJERUp = 1.228; cJERDn = 1.042; }
     else if ( jetEta < 5.0 ) { cJER = 1.288; cJERUp = 1.488; cJERDn = 1.089; }
     else { cJER = cJERUp = cJERDn = 1; }
-  }
+  };
 
 };
 
-KJetUncProducer::KJetUncProducer(const edm::ParameterSet& pset)
+KJetMetUncSelector::KJetMetUncSelector(const edm::ParameterSet& pset)
 {
   isMC_ = pset.getParameter<bool>("isMC");
 
@@ -65,6 +89,26 @@ KJetUncProducer::KJetUncProducer(const edm::ParameterSet& pset)
   const std::string jecLevel = pset.getParameter<string>("jecLevel");
   jecUncCalculator_ = new JetCorrectionUncertainty(JetCorrectorParameters(jecFilePath.fullPath(), jecLevel));
 
+  // Selection cuts
+  edm::ParameterSet selectionPSet = pset.getParameter<edm::ParameterSet>("selection");
+  isGoodJet_ = new PFJetIDSelectionFunctor(selectionPSet.getParameter<edm::ParameterSet>("jetId"));
+  minPt_ = selectionPSet.getParameter<double>("minPt");
+  maxEta_ = selectionPSet.getParameter<double>("maxEta");
+
+  // Cleaning
+  edm::ParameterSet cleanPSet = pset.getParameter<edm::ParameterSet>("cleaning");
+  doClean_ = cleanPSet.getParameter<bool>("doClean");
+  if ( doClean_ )
+  {
+    overlapDeltaR_ = cleanPSet.getParameter<double>("overlapDeltaR");
+    overlapCandLabels_ = cleanPSet.getParameter<std::vector<edm::InputTag> >("overlapCands");
+    if ( overlapCandLabels_.empty() ) doClean_ = false;
+  }
+
+  minNumber_ = pset.getParameter<unsigned int>("minNumber");
+  maxNumber_ = pset.getParameter<unsigned int>("maxNumber");
+
+  produces<std::vector<pat::Jet> >();
   produces<pat::JetToValue>("up");
   produces<pat::JetToValue>("dn");
   produces<METs>("up");
@@ -78,44 +122,65 @@ KJetUncProducer::KJetUncProducer(const edm::ParameterSet& pset)
     produces<METs>("resUp");
     produces<METs>("resDn");
   }
+
 }
 
-void KJetUncProducer::produce(edm::Event& event, const edm::EventSetup& eventSetup)
+bool KJetMetUncSelector::filter(edm::Event& event, const edm::EventSetup& eventSetup)
 {
   edm::Handle<Jets> jetHandle;
   event.getByLabel(jetLabel_, jetHandle);
 
   edm::Handle<METs> metHandle;
   event.getByLabel(metLabel_, metHandle);
 
+  // Prepare JEC and JER factors
   std::vector<double> fJECsUp;
   std::vector<double> fJECsDn;
   std::vector<double> fJERs;
   std::vector<double> fJERsUp;
   std::vector<double> fJERsDn;
 
-  pat::MET met = metHandle->at(0);
+  // Prepare MET
+  const pat::MET& met = metHandle->at(0);
   const double metX = met.px(), metY = met.py();
   double metUpX = metX, metUpY = metY;
   double metDnX = metX, metDnY = metY;
   double metResX   = metX, metResY   = metY;
   double metResUpX = metX, metResUpY = metY;
   double metResDnX = metX, metResDnY = metY;
 
+  // Collect candidates for overlap removal
+  std::vector<const reco::Candidate*> overlapCands;
+  if ( doClean_ )
+  {
+    for ( int iLabel=0, nLabel=overlapCandLabels_.size(); iLabel<nLabel; ++iLabel )
+    {
+      edm::Handle<edm::View<reco::Candidate> > overlapCandHandle;
+      event.getByLabel(overlapCandLabels_.at(iLabel), overlapCandHandle);
+
+      //if ( !overlapCandHandle.isValid() ) continue;
+
+      for ( int i=0, n=overlapCandHandle->size(); i<n; ++i )
+      {
+        overlapCands.push_back(&(overlapCandHandle->at(i)));
+      }
+    }
+  }
+
+  // Now start to build jet collection
+  std::auto_ptr<std::vector<pat::Jet> > cleanJets(new std::vector<pat::Jet>());
+
   for ( int i=0, n=jetHandle->size(); i<n; ++i )
   {
     pat::Jet jet = jetHandle->at(i);
-    reco::Candidate::LorentzVector jetP4 = jet.p4();
+    if ( !(*isGoodJet_)(jet) ) continue;
+    const reco::Candidate::LorentzVector jetP4 = jet.p4();
     const double jetPt = jetP4.pt();
+    const double jetEta = jetP4.eta();
 
-    // JEC and uncertainties
-    jecUncCalculator_->setJetPt(jetP4.pt());
-    jecUncCalculator_->setJetEta(jetP4.eta());
-    const double fJECUp = 1+jecUncCalculator_->getUncertainty(true);
-    jecUncCalculator_->setJetPt(jetP4.pt());
-    jecUncCalculator_->setJetEta(jetP4.eta());
-    const double fJECDn = 1-jecUncCalculator_->getUncertainty(false);
-
+    // Calculate JEC uncertanties
+    double fJECUp, fJECDn;
+    loadJEC(jetPt, jetEta, fJECUp, fJECDn);
     math::XYZTLorentzVector jetUpP4 = jetP4*fJECUp;
     math::XYZTLorentzVector jetDnP4 = jetP4*fJECDn;
 
@@ -131,17 +196,16 @@ void KJetUncProducer::produce(edm::Event& event, const edm::EventSetup& eventSet
       const reco::GenJet* genJet = jet.genJet();
       if ( genJet and genJet->pt() > 10 )
       {
+        double cJER, cJERUp, cJERDn;
+        loadJER(jetEta, cJER, cJERUp, cJERDn);
+
         const math::XYZTLorentzVector& rawJetP4 = jet.correctedP4(0);
         const double rawPx = rawJetP4.px();
         const double rawPy = rawJetP4.py();
 
         const double genJetPt = genJet->pt();
         const double dPt = jetPt-genJetPt;
 
-        const double jetEta = std::abs(jet.eta());
-        double cJER, cJERUp, cJERDn;
-        getJERFactor(jetEta, cJER, cJERUp, cJERDn);
-
         fJER   = max(0., (genJetPt+dPt*cJER  )/jetPt);
         fJERUp = max(0., (genJetPt+dPt*cJERUp)/jetPt);
         fJERDn = max(0., (genJetPt+dPt*cJERDn)/jetPt);
@@ -163,8 +227,25 @@ void KJetUncProducer::produce(edm::Event& event, const edm::EventSetup& eventSet
       }
     }
 
-    // Put JES,JER factors
-    //edm::Ref<Jets> jetRef(jetHandle, i);
+    // Check overlap
+    if ( doClean_ )
+    {
+      bool isOverlap = false;
+      for ( int j=0, m=overlapCands.size(); j<m; ++j )
+      {
+        if ( deltaR(jet.p4(), overlapCands.at(j)->p4()) < overlapDeltaR_ )
+        {
+          isOverlap = true;
+        }
+      }
+      if ( isOverlap ) continue;
+    }
+
+    // Check acceptance
+    if ( std::abs(jetEta) <= maxEta_ ) continue;
+    if ( jetPt*fJECDn*fJERDn < minPt_ ) continue;
+
+    cleanJets->push_back(jet);
     fJECsUp.push_back(fJECUp);
     fJECsDn.push_back(fJECDn);
     if ( isMC_ )
@@ -175,21 +256,23 @@ void KJetUncProducer::produce(edm::Event& event, const edm::EventSetup& eventSet
     }
   }
 
-  pat::MET metUp, metDn;
-  metUp.setP4(reco::Candidate::LorentzVector(metUpX, metUpY, 0, hypot(metUpX, metUpY)));
-  metDn.setP4(reco::Candidate::LorentzVector(metDnX, metDnY, 0, hypot(metDnX, metDnY)));
+  const unsigned int nCleanJet = cleanJets->size();
+  edm::OrphanHandle<pat::JetCollection> outColl = event.put(cleanJets);
 
   std::auto_ptr<pat::JetToValue> fJECMapUp(new pat::JetToValue);
   std::auto_ptr<pat::JetToValue> fJECMapDn(new pat::JetToValue);
   pat::JetToValue::Filler fillJECUp(*fJECMapUp);
   pat::JetToValue::Filler fillJECDn(*fJECMapDn);
-  fillJECUp.insert(jetHandle, fJECsUp.begin(), fJECsUp.end());
-  fillJECDn.insert(jetHandle, fJECsDn.begin(), fJECsDn.end());
+  fillJECUp.insert(outColl, fJECsUp.begin(), fJECsUp.end());
+  fillJECDn.insert(outColl, fJECsDn.begin(), fJECsDn.end());
   fillJECUp.fill();
   fillJECDn.fill();
   event.put(fJECMapUp, "up");
   event.put(fJECMapDn, "dn");
 
+  pat::MET metUp, metDn;
+  metUp.setP4(reco::Candidate::LorentzVector(metUpX, metUpY, 0, hypot(metUpX, metUpY)));
+  metDn.setP4(reco::Candidate::LorentzVector(metDnX, metDnY, 0, hypot(metDnX, metDnY)));
   std::auto_ptr<METs> metsUp(new METs);
   std::auto_ptr<METs> metsDn(new METs);
   metsUp->push_back(metUp);
@@ -205,9 +288,9 @@ void KJetUncProducer::produce(edm::Event& event, const edm::EventSetup& eventSet
     pat::JetToValue::Filler fillJER(*fJERMap);
     pat::JetToValue::Filler fillJERUp(*fJERMapUp);
     pat::JetToValue::Filler fillJERDn(*fJERMapDn);
-    fillJER.insert(jetHandle, fJERs.begin(), fJERs.end());
-    fillJERUp.insert(jetHandle, fJERsUp.begin(), fJERsUp.end());
-    fillJERDn.insert(jetHandle, fJERsDn.begin(), fJERsDn.end());
+    fillJER.insert(outColl, fJERs.begin(), fJERs.end());
+    fillJERUp.insert(outColl, fJERsUp.begin(), fJERsUp.end());
+    fillJERDn.insert(outColl, fJERsDn.begin(), fJERsDn.end());
     fillJER.fill();
     fillJERUp.fill();
     fillJERDn.fill();
@@ -231,7 +314,10 @@ void KJetUncProducer::produce(edm::Event& event, const edm::EventSetup& eventSet
     event.put(metsResDn, "resDn");
   }
 
+  if ( nCleanJet < minNumber_ or nCleanJet > maxNumber_ ) return false;
+
+  return true;
 }
 
-DEFINE_FWK_MODULE(KJetUncProducer);
+DEFINE_FWK_MODULE(KJetMetUncSelector);