PatExamples/plugins/JetEnergyShift.cc

0001 #include "FWCore/Framework/interface/Event.h"
0002 #include "FWCore/Framework/interface/global/EDProducer.h"
0003 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0004
0005 #include "DataFormats/PatCandidates/interface/Jet.h"
0006 #include "DataFormats/PatCandidates/interface/MET.h"
0007
0008 /**
0009    \class   JetEnergyShift JetEnergyShift.h "PhysicsTools/PatExamples/plugins/JetEnergyShift.h"
0010
0011    \brief   Plugin to shift the jet energy scale and recalculate the MET accordingly
0012
0013    Plugin to shift the jet energy scale and recalculate the MET accordingly. The module
0014    mimics the assumption that the jet energy scale (JES) has been estimated wrong by a
0015    factor of _scaleFactor_, corresponding to a L2L3 corrected jet. The p4 of the patJet
0016    is beeing rescaled. All other patJet properties stay the same. The MET is recalculated
0017    taking the shifted JES into account for the Type1 MET correction. For the patMET the
0018    rescaled sumET and the p4 are stored. The different correction levels are lost for
0019    the new collection. The module has the following parameters:
0020
0021   inputJets            : input collection for  MET (expecting patMET).
0022   inputMETs            : input collection for jets (expecting patJets).
0023   scaleFactor          : scale factor to which to shift the JES.
0024   jetPTThresholdForMET : pt threshold for (uncorrected!) jets considered for Type1 MET
0025                          corrections.
0026   jetEMLimitForMET     : limit in em fraction for Type1 MET correction.
0027
0028   For expected parameters for _jetPTThresholdForMET_ and _jetEMLimitForMET_ have a look
0029   at: JetMETCorrections/Type1MET/python/MetType1Corrections_cff.py. Two output collections
0030   are written to file with instance label corresponding to the input label of the jet
0031   and met input collections.
0032 */
0033
0034 class JetEnergyShift : public edm::global::EDProducer<> {
0035 public:
0036   /// default constructor
0037   explicit JetEnergyShift(const edm::ParameterSet&);
0038
0039 private:
0040   /// rescale jet energy and recalculated MET
0041   void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0042
0043 private:
0044   /// jet input collection
0045   edm::EDGetTokenT<std::vector<pat::Jet>> inputJetsToken_;
0046   /// met input collection
0047   edm::EDGetTokenT<std::vector<pat::MET>> inputMETsToken_;
0048   /// jet output collection
0049   std::string outputJets_;
0050   /// MET output collection
0051   std::string outputMETs_;
0052   /// scale factor for the rescaling
0053   double scaleFactor_;
0054   /// threshold on (raw!) jet pt for Type1 MET corrections
0055   double jetPTThresholdForMET_;
0056   /// limit on the emf of the jet for Type1 MET corrections
0057   double jetEMLimitForMET_;
0058 };
0059
0060 JetEnergyShift::JetEnergyShift(const edm::ParameterSet& cfg)
0061     : inputJetsToken_(consumes<std::vector<pat::Jet>>(cfg.getParameter<edm::InputTag>("inputJets"))),
0062       inputMETsToken_(consumes<std::vector<pat::MET>>(cfg.getParameter<edm::InputTag>("inputMETs"))),
0063       scaleFactor_(cfg.getParameter<double>("scaleFactor")),
0064       jetPTThresholdForMET_(cfg.getParameter<double>("jetPTThresholdForMET")),
0065       jetEMLimitForMET_(cfg.getParameter<double>("jetEMLimitForMET")) {
0066   // use label of input to create label for output
0067   outputJets_ = cfg.getParameter<edm::InputTag>("inputJets").label();
0068   outputMETs_ = cfg.getParameter<edm::InputTag>("inputMETs").label();
0069   // register products
0070   produces<std::vector<pat::Jet>>(outputJets_);
0071   produces<std::vector<pat::MET>>(outputMETs_);
0072 }
0073
0074 void JetEnergyShift::produce(edm::StreamID, edm::Event& event, const edm::EventSetup& setup) const {
0075   edm::Handle<std::vector<pat::Jet>> jets;
0076   event.getByToken(inputJetsToken_, jets);
0077
0078   edm::Handle<std::vector<pat::MET>> mets;
0079   event.getByToken(inputMETsToken_, mets);
0080
0081   auto pJets = std::make_unique<std::vector<pat::Jet>>();
0082   auto pMETs = std::make_unique<std::vector<pat::MET>>();
0083
0084   double dPx = 0.;
0085   double dPy = 0.;
0086   double dSumEt = 0.;
0087
0088   for (std::vector<pat::Jet>::const_iterator jet = jets->begin(); jet != jets->end(); ++jet) {
0089     pat::Jet scaledJet = *jet;
0090     scaledJet.scaleEnergy(scaleFactor_);
0091     pJets->push_back(scaledJet);
0092     // consider jet scale shift only if the raw jet pt and emf
0093     // is above the thresholds given in the module definition
0094     if (jet->correctedJet("raw").pt() > jetPTThresholdForMET_ && jet->emEnergyFraction() < jetEMLimitForMET_) {
0095       dPx += scaledJet.px() - jet->px();
0096       dPy += scaledJet.py() - jet->py();
0097       dSumEt += scaledJet.et() - jet->et();
0098     }
0099   }
0100
0101   // scale MET accordingly
0102   pat::MET met = *(mets->begin());
0103   double scaledMETPx = met.px() - dPx;
0104   double scaledMETPy = met.py() - dPy;
0105   pat::MET scaledMET(
0106       reco::MET(met.sumEt() + dSumEt,
0107                 reco::MET::LorentzVector(
0108                     scaledMETPx, scaledMETPy, 0, sqrt(scaledMETPx * scaledMETPx + scaledMETPy * scaledMETPy)),
0109                 reco::MET::Point(0, 0, 0)));
0110   pMETs->push_back(scaledMET);
0111   event.put(std::move(pJets), outputJets_);
0112   event.put(std::move(pMETs), outputMETs_);
0113 }
0114
0115 #include "FWCore/Framework/interface/MakerMacros.h"
0116 DEFINE_FWK_MODULE(JetEnergyShift);