PatUtils/plugins/ShiftedPFCandidateProducerForPFNoPUMEt.cc

0001 /** \class ShiftedPFCandidateProducerForPFNoPUMEt
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/PFNoPUMETProducer.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/stream/EDProducer.h"
0029 #include "FWCore/Framework/interface/ESHandle.h"
0030 #include "FWCore/Framework/interface/Event.h"
0031 #include "FWCore/Framework/interface/EventSetup.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 ShiftedPFCandidateProducerForPFNoPUMEt : public edm::stream::EDProducer<> {
0041 public:
0042   explicit ShiftedPFCandidateProducerForPFNoPUMEt(const edm::ParameterSet&);
0043
0044 private:
0045   void produce(edm::Event&, const edm::EventSetup&) override;
0046
0047   std::string moduleLabel_;
0048
0049   edm::EDGetTokenT<reco::PFCandidateCollection> srcPFCandidatesToken_;
0050   edm::EDGetTokenT<reco::PFJetCollection> srcJetsToken_;
0051
0052   edm::FileInPath jetCorrInputFileName_;
0053   std::string jetCorrPayloadName_;
0054   edm::ESGetToken<JetCorrectorParametersCollection, JetCorrectionsRecord> jetCorrPayloadToken_;
0055   std::string jetCorrUncertaintyTag_;
0056   std::unique_ptr<JetCorrectorParameters> jetCorrParameters_;
0057   std::unique_ptr<JetCorrectionUncertainty> jecUncertainty_;
0058
0059   bool jecValidFileName_;
0060
0061   double minJetPt_;
0062
0063   double shiftBy_;
0064
0065   double unclEnUncertainty_;
0066 };
0067
0068 namespace {
0069   constexpr double dR2Match = 0.01 * 0.01;
0070 }
0071
0072 ShiftedPFCandidateProducerForPFNoPUMEt::ShiftedPFCandidateProducerForPFNoPUMEt(const edm::ParameterSet& cfg)
0073     : srcPFCandidatesToken_(consumes<reco::PFCandidateCollection>(cfg.getParameter<edm::InputTag>("srcPFCandidates"))),
0074       srcJetsToken_(consumes<reco::PFJetCollection>(cfg.getParameter<edm::InputTag>("srcJets"))) {
0075   jetCorrUncertaintyTag_ = cfg.getParameter<std::string>("jetCorrUncertaintyTag");
0076
0077   jecValidFileName_ = cfg.exists("jetCorrInputFileName");
0078   if (jecValidFileName_) {
0079     jetCorrInputFileName_ = cfg.getParameter<edm::FileInPath>("jetCorrInputFileName");
0080     if (jetCorrInputFileName_.location() == edm::FileInPath::Unknown)
0081       throw cms::Exception("ShiftedJetProducerT")
0082           << " Failed to find JEC parameter file = " << jetCorrInputFileName_ << " !!\n";
0083     edm::LogInfo("ShiftedPFCandidateProducerForPFNoPUMEt")
0084         << "Reading JEC parameters = " << jetCorrUncertaintyTag_ << " from file = " << jetCorrInputFileName_.fullPath()
0085         << "." << std::endl;
0086     jetCorrParameters_ =
0087         std::make_unique<JetCorrectorParameters>(jetCorrInputFileName_.fullPath(), jetCorrUncertaintyTag_);
0088     jecUncertainty_ = std::make_unique<JetCorrectionUncertainty>(*jetCorrParameters_);
0089   } else {
0090     edm::LogInfo("ShiftedPFCandidateProducerForPFNoPUMEt")
0091         << "Reading JEC parameters = " << jetCorrUncertaintyTag_ << " from DB/SQLlite file." << std::endl;
0092     jetCorrPayloadName_ = cfg.getParameter<std::string>("jetCorrPayloadName");
0093     jetCorrPayloadToken_ = esConsumes(edm::ESInputTag("", jetCorrPayloadName_));
0094   }
0095
0096   minJetPt_ = cfg.getParameter<double>("minJetPt");
0097
0098   shiftBy_ = cfg.getParameter<double>("shiftBy");
0099
0100   unclEnUncertainty_ = cfg.getParameter<double>("unclEnUncertainty");
0101
0102   produces<reco::PFCandidateCollection>();
0103 }
0104
0105 void ShiftedPFCandidateProducerForPFNoPUMEt::produce(edm::Event& evt, const edm::EventSetup& es) {
0106   edm::Handle<reco::PFCandidateCollection> originalPFCandidates;
0107   evt.getByToken(srcPFCandidatesToken_, originalPFCandidates);
0108
0109   edm::Handle<reco::PFJetCollection> jets;
0110   evt.getByToken(srcJetsToken_, jets);
0111
0112   std::vector<const reco::PFJet*> selectedJets;
0113   for (reco::PFJetCollection::const_iterator jet = jets->begin(); jet != jets->end(); ++jet) {
0114     if (jet->pt() > minJetPt_)
0115       selectedJets.push_back(&(*jet));
0116   }
0117
0118   if (!jetCorrPayloadName_.empty()) {
0119     JetCorrectorParametersCollection const& jetCorrParameterSet = es.getData(jetCorrPayloadToken_);
0120     const JetCorrectorParameters& jetCorrParameters = (jetCorrParameterSet)[jetCorrUncertaintyTag_];
0121     jecUncertainty_ = std::make_unique<JetCorrectionUncertainty>(jetCorrParameters);
0122   }
0123
0124   auto shiftedPFCandidates = std::make_unique<reco::PFCandidateCollection>();
0125
0126   for (reco::PFCandidateCollection::const_iterator originalPFCandidate = originalPFCandidates->begin();
0127        originalPFCandidate != originalPFCandidates->end();
0128        ++originalPFCandidate) {
0129     const reco::PFJet* jet_matched = nullptr;
0130     for (auto jet : selectedJets) {
0131       for (const auto& jetc : jet->getPFConstituents()) {
0132         if (deltaR2(originalPFCandidate->p4(), jetc->p4()) < dR2Match) {
0133           jet_matched = jet;
0134           break;
0135         }
0136       }
0137       if (jet_matched)
0138         break;
0139     }
0140
0141     double shift = 0.;
0142     if (jet_matched != nullptr) {
0143       jecUncertainty_->setJetEta(jet_matched->eta());
0144       jecUncertainty_->setJetPt(jet_matched->pt());
0145
0146       shift = jecUncertainty_->getUncertainty(true);
0147     } else {
0148       shift = unclEnUncertainty_;
0149     }
0150
0151     shift *= shiftBy_;
0152
0153     reco::Candidate::LorentzVector shiftedPFCandidateP4 = originalPFCandidate->p4();
0154     shiftedPFCandidateP4 *= (1. + shift);
0155
0156     reco::PFCandidate shiftedPFCandidate(*originalPFCandidate);
0157     shiftedPFCandidate.setP4(shiftedPFCandidateP4);
0158
0159     shiftedPFCandidates->push_back(shiftedPFCandidate);
0160   }
0161
0162   evt.put(std::move(shiftedPFCandidates));
0163 }
0164
0165 #include "FWCore/Framework/interface/MakerMacros.h"
0166
0167 DEFINE_FWK_MODULE(ShiftedPFCandidateProducerForPFNoPUMEt);