PatUtils/plugins/ShiftedPFCandidateProducerForNoPileUpPFMEt.cc

0001 /** \class ShiftedPFCandidateProducerForNoPileUpPFMEt
0002  *
0003  * Vary energy of PFCandidates which are (are not) within jets of Pt > 10 GeV
0004  * by jet energy uncertainty (by 10% "unclustered" energy uncertainty)
0005  *
0006  * NOTE: Auxiliary class specific to estimating systematic uncertainty
0007  *       on PFMET reconstructed by no-PU MET reconstruction algorithm
0008  *      (implemented in JetMETCorrections/Type1MET/src/NoPileUpPFMETProducer.cc)
0009  *
0010  *       In case all PFCandidates not within jets of Pt > 30 GeV would be varied
0011  *       by the 10% "unclustered" energy uncertainty, the systematic uncertainty
0012  *       on the reconstructed no-PU MET would be overestimated significantly !!
0013  *
0014  * \author Christian Veelken, LLR
0015  *
0016  *
0017  *
0018  */
0019
0020 #include "CondFormats/JetMETObjects/interface/JetCorrectionUncertainty.h"
0021 #include "CondFormats/JetMETObjects/interface/JetCorrectorParameters.h"
0022 #include "DataFormats/Candidate/interface/Candidate.h"
0023 #include "DataFormats/JetReco/interface/PFJet.h"
0024 #include "DataFormats/JetReco/interface/PFJetCollection.h"
0025 #include "DataFormats/Math/interface/deltaR.h"
0026 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0027 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidateFwd.h"
0028 #include "FWCore/Framework/interface/ESHandle.h"
0029 #include "FWCore/Framework/interface/Event.h"
0030 #include "FWCore/Framework/interface/EventSetup.h"
0031 #include "FWCore/Framework/interface/stream/EDProducer.h"
0032 #include "FWCore/ParameterSet/interface/FileInPath.h"
0033 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0034 #include "FWCore/Utilities/interface/InputTag.h"
0035 #include "JetMETCorrections/Objects/interface/JetCorrectionsRecord.h"
0036
0037 #include <string>
0038 #include <vector>
0039
0040 class ShiftedPFCandidateProducerForNoPileUpPFMEt : public edm::stream::EDProducer<> {
0041 public:
0042   explicit ShiftedPFCandidateProducerForNoPileUpPFMEt(const edm::ParameterSet&);
0043
0044 private:
0045   void produce(edm::Event&, const edm::EventSetup&) override;
0046
0047   edm::EDGetTokenT<reco::PFCandidateCollection> srcPFCandidatesToken_;
0048   edm::EDGetTokenT<reco::PFJetCollection> srcJetsToken_;
0049
0050   edm::FileInPath jetCorrInputFileName_;
0051   std::string jetCorrPayloadName_;
0052   edm::ESGetToken<JetCorrectorParametersCollection, JetCorrectionsRecord> jetCorrPayloadToken_;
0053   std::string jetCorrUncertaintyTag_;
0054   std::unique_ptr<JetCorrectorParameters> jetCorrParameters_;
0055   std::unique_ptr<JetCorrectionUncertainty> jecUncertainty_;
0056
0057   double minJetPt_;
0058
0059   double shiftBy_;
0060
0061   double unclEnUncertainty_;
0062 };
0063
0064 ShiftedPFCandidateProducerForNoPileUpPFMEt::ShiftedPFCandidateProducerForNoPileUpPFMEt(const edm::ParameterSet& cfg)
0065     : srcPFCandidatesToken_(consumes<reco::PFCandidateCollection>(cfg.getParameter<edm::InputTag>("srcPFCandidates"))),
0066       srcJetsToken_(consumes<reco::PFJetCollection>(cfg.getParameter<edm::InputTag>("srcJets"))) {
0067   jetCorrUncertaintyTag_ = cfg.getParameter<std::string>("jetCorrUncertaintyTag");
0068   if (cfg.exists("jetCorrInputFileName")) {
0069     jetCorrInputFileName_ = cfg.getParameter<edm::FileInPath>("jetCorrInputFileName");
0070     if (jetCorrInputFileName_.location() == edm::FileInPath::Unknown)
0071       throw cms::Exception("ShiftedJetProducerT")
0072           << " Failed to find JEC parameter file = " << jetCorrInputFileName_ << " !!\n";
0073     jetCorrParameters_ =
0074         std::make_unique<JetCorrectorParameters>(jetCorrInputFileName_.fullPath(), jetCorrUncertaintyTag_);
0075     jecUncertainty_ = std::make_unique<JetCorrectionUncertainty>(*jetCorrParameters_);
0076   } else {
0077     jetCorrPayloadName_ = cfg.getParameter<std::string>("jetCorrPayloadName");
0078     jetCorrPayloadToken_ = esConsumes(edm::ESInputTag("", jetCorrPayloadName_));
0079   }
0080
0081   minJetPt_ = cfg.getParameter<double>("minJetPt");
0082
0083   shiftBy_ = cfg.getParameter<double>("shiftBy");
0084
0085   unclEnUncertainty_ = cfg.getParameter<double>("unclEnUncertainty");
0086
0087   produces<reco::PFCandidateCollection>();
0088 }
0089
0090 void ShiftedPFCandidateProducerForNoPileUpPFMEt::produce(edm::Event& evt, const edm::EventSetup& es) {
0091   edm::Handle<reco::PFCandidateCollection> originalPFCandidates;
0092   evt.getByToken(srcPFCandidatesToken_, originalPFCandidates);
0093
0094   edm::Handle<reco::PFJetCollection> jets;
0095   evt.getByToken(srcJetsToken_, jets);
0096
0097   std::vector<const reco::PFJet*> selectedJets;
0098   for (reco::PFJetCollection::const_iterator jet = jets->begin(); jet != jets->end(); ++jet) {
0099     if (jet->pt() > minJetPt_)
0100       selectedJets.push_back(&(*jet));
0101   }
0102
0103   if (!jetCorrPayloadName_.empty()) {
0104     const JetCorrectorParametersCollection& jetCorrParameterSet = es.getData(jetCorrPayloadToken_);
0105     const JetCorrectorParameters& jetCorrParameters = (jetCorrParameterSet)[jetCorrUncertaintyTag_];
0106     jecUncertainty_ = std::make_unique<JetCorrectionUncertainty>(jetCorrParameters);
0107   }
0108
0109   auto shiftedPFCandidates = std::make_unique<reco::PFCandidateCollection>();
0110   for (reco::PFCandidateCollection::const_iterator originalPFCandidate = originalPFCandidates->begin();
0111        originalPFCandidate != originalPFCandidates->end();
0112        ++originalPFCandidate) {
0113     const reco::PFJet* jet_matched = nullptr;
0114     for (std::vector<const reco::PFJet*>::iterator jet = selectedJets.begin(); jet != selectedJets.end(); ++jet) {
0115       std::vector<reco::PFCandidatePtr> jetConstituents = (*jet)->getPFConstituents();
0116       for (std::vector<reco::PFCandidatePtr>::const_iterator jetConstituent = jetConstituents.begin();
0117            jetConstituent != jetConstituents.end() && !jet_matched;
0118            ++jetConstituent) {
0119         if (deltaR(originalPFCandidate->p4(), (*jetConstituent)->p4()) < 1.e-2)
0120           jet_matched = (*jet);
0121       }
0122     }
0123
0124     double shift = 0.;
0125     if (jet_matched) {
0126       jecUncertainty_->setJetEta(jet_matched->eta());
0127       jecUncertainty_->setJetPt(jet_matched->pt());
0128
0129       shift = jecUncertainty_->getUncertainty(true);
0130     } else {
0131       shift = unclEnUncertainty_;
0132     }
0133
0134     shift *= shiftBy_;
0135
0136     reco::Candidate::LorentzVector shiftedPFCandidateP4 = originalPFCandidate->p4();
0137     shiftedPFCandidateP4 *= (1. + shift);
0138
0139     reco::PFCandidate shiftedPFCandidate(*originalPFCandidate);
0140     shiftedPFCandidate.setP4(shiftedPFCandidateP4);
0141
0142     shiftedPFCandidates->push_back(shiftedPFCandidate);
0143   }
0144
0145   evt.put(std::move(shiftedPFCandidates));
0146 }
0147
0148 #include "FWCore/Framework/interface/MakerMacros.h"
0149
0150 DEFINE_FWK_MODULE(ShiftedPFCandidateProducerForNoPileUpPFMEt);