JetEnergyShift.cc

Go to the documentation of this file.

#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/global/EDProducer.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"

#include "DataFormats/PatCandidates/interface/Jet.h"
#include "DataFormats/PatCandidates/interface/MET.h"

class JetEnergyShift : public edm::global::EDProducer<> {
public:
  explicit JetEnergyShift(const edm::ParameterSet&);

private:
  void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;

private:
  edm::EDGetTokenT<std::vector<pat::Jet>> inputJetsToken_;
  edm::EDGetTokenT<std::vector<pat::MET>> inputMETsToken_;
  std::string outputJets_;
  std::string outputMETs_;
  double scaleFactor_;
  double jetPTThresholdForMET_;
  double jetEMLimitForMET_;
};

JetEnergyShift::JetEnergyShift(const edm::ParameterSet& cfg)
    : inputJetsToken_(consumes<std::vector<pat::Jet>>(cfg.getParameter<edm::InputTag>("inputJets"))),
      inputMETsToken_(consumes<std::vector<pat::MET>>(cfg.getParameter<edm::InputTag>("inputMETs"))),
      scaleFactor_(cfg.getParameter<double>("scaleFactor")),
      jetPTThresholdForMET_(cfg.getParameter<double>("jetPTThresholdForMET")),
      jetEMLimitForMET_(cfg.getParameter<double>("jetEMLimitForMET")) {
  // use label of input to create label for output
  outputJets_ = cfg.getParameter<edm::InputTag>("inputJets").label();
  outputMETs_ = cfg.getParameter<edm::InputTag>("inputMETs").label();
  // register products
  produces<std::vector<pat::Jet>>(outputJets_);
  produces<std::vector<pat::MET>>(outputMETs_);
}

void JetEnergyShift::produce(edm::StreamID, edm::Event& event, const edm::EventSetup& setup) const {
  edm::Handle<std::vector<pat::Jet>> jets;
  event.getByToken(inputJetsToken_, jets);

  edm::Handle<std::vector<pat::MET>> mets;
  event.getByToken(inputMETsToken_, mets);

  auto pJets = std::make_unique<std::vector<pat::Jet>>();
  auto pMETs = std::make_unique<std::vector<pat::MET>>();

  double dPx = 0.;
  double dPy = 0.;
  double dSumEt = 0.;

  for (std::vector<pat::Jet>::const_iterator jet = jets->begin(); jet != jets->end(); ++jet) {
    pat::Jet scaledJet = *jet;
    scaledJet.scaleEnergy(scaleFactor_);
    pJets->push_back(scaledJet);
    // consider jet scale shift only if the raw jet pt and emf
    // is above the thresholds given in the module definition
    if (jet->correctedJet("raw").pt() > jetPTThresholdForMET_ && jet->emEnergyFraction() < jetEMLimitForMET_) {
      dPx += scaledJet.px() - jet->px();
      dPy += scaledJet.py() - jet->py();
      dSumEt += scaledJet.et() - jet->et();
    }
  }

  // scale MET accordingly
  pat::MET met = *(mets->begin());
  double scaledMETPx = met.px() - dPx;
  double scaledMETPy = met.py() - dPy;
  pat::MET scaledMET(
      reco::MET(met.sumEt() + dSumEt,
                reco::MET::LorentzVector(
                    scaledMETPx, scaledMETPy, 0, sqrt(scaledMETPx * scaledMETPx + scaledMETPy * scaledMETPy)),
                reco::MET::Point(0, 0, 0)));
  pMETs->push_back(scaledMET);
  event.put(std::move(pJets), outputJets_);
  event.put(std::move(pMETs), outputMETs_);
}

#include "FWCore/Framework/interface/MakerMacros.h"
DEFINE_FWK_MODULE(JetEnergyShift);