Type1MET/plugins/Type1PFMET.cc

0001 /**\class Type1PFMET
0002 \brief Computes the Type-1 corrections for pfMET. A specific version of the Type1MET class from the JetMETCorrections/Type1MET package.
0003
0004 \todo Unify with the Type1MET class from the JetMETCorrections/Type1MET package
0005
0006 \author Michal Bluj
0007 \date   February 2009
0008 */
0009
0010 #include "DataFormats/METReco/interface/MET.h"
0011 #include "DataFormats/METReco/interface/CaloMET.h"
0012 #include "DataFormats/JetReco/interface/PFJet.h"
0013 #include "FWCore/Framework/interface/Frameworkfwd.h"
0014 #include "FWCore/Framework/interface/global/EDProducer.h"
0015 #include "FWCore/Framework/interface/Event.h"
0016 #include "FWCore/Framework/interface/MakerMacros.h"
0017 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0018 #include "DataFormats/Common/interface/Handle.h"
0019 #include "DataFormats/JetReco/interface/PFJetCollection.h"
0020 #include "DataFormats/METReco/interface/METFwd.h"
0021 #include "JetMETCorrections/JetCorrector/interface/JetCorrector.h"
0022 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0023 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0024
0025 #include <memory>
0026 #include <cstring>
0027
0028 class Type1PFMET : public edm::global::EDProducer<> {
0029 public:
0030   explicit Type1PFMET(const edm::ParameterSet&);
0031   explicit Type1PFMET();
0032   ~Type1PFMET() override;
0033   void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0034
0035   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0036
0037 private:
0038   edm::EDGetTokenT<reco::METCollection> tokenUncorMet;
0039   edm::EDGetTokenT<reco::PFJetCollection> tokenUncorJets;
0040   edm::EDGetTokenT<reco::JetCorrector> correctorToken;
0041   double jetPTthreshold;
0042   double jetEMfracLimit;
0043   double jetMufracLimit;
0044   void run(const reco::METCollection& uncorMET,
0045            const reco::JetCorrector& corrector,
0046            const reco::PFJetCollection& uncorJet,
0047            double jetPTthreshold,
0048            double jetEMfracLimit,
0049            double jetMufracLimit,
0050            reco::METCollection* corMET) const;
0051 };
0052
0053 #include "FWCore/Framework/interface/MakerMacros.h"
0054 DEFINE_FWK_MODULE(Type1PFMET);
0055
0056 void Type1PFMET::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0057   // pfType1MET
0058   // Type-1 met corrections (AK4PFJets)
0059   // remember about including ES producer definition e.g. JetMETCorrections.Configuration.L2L3Corrections_Summer08Redigi_cff
0060   edm::ParameterSetDescription desc;
0061   desc.add<edm::InputTag>("inputUncorJetsTag", edm::InputTag("ak4PFJets"));
0062   desc.add<double>("jetEMfracLimit", 0.95);  // to remove electron which give rise to jets
0063   desc.add<double>("jetMufracLimit", 0.95);  // to remove muon which give rise to jets
0064   desc.add<std::string>("metType", "PFMET");
0065   desc.add<double>("jetPTthreshold", 20.0);
0066   // pfMET should be not corrected for HF 0.7
0067   desc.add<edm::InputTag>("inputUncorMetLabel", edm::InputTag("pfMET"));
0068   desc.add<edm::InputTag>("corrector", edm::InputTag("ak4PFL2L3Corrector"));
0069   descriptions.add("pfType1MET", desc);
0070 }
0071
0072 using namespace reco;
0073
0074 // PRODUCER CONSTRUCTORS ------------------------------------------
0075 Type1PFMET::Type1PFMET(const edm::ParameterSet& iConfig) {
0076   tokenUncorMet = consumes<METCollection>(iConfig.getParameter<edm::InputTag>("inputUncorMetLabel"));
0077   tokenUncorJets = consumes<PFJetCollection>(iConfig.getParameter<edm::InputTag>("inputUncorJetsTag"));
0078   correctorToken = consumes<JetCorrector>(iConfig.getParameter<edm::InputTag>("corrector"));
0079   jetPTthreshold = iConfig.getParameter<double>("jetPTthreshold");
0080   jetEMfracLimit = iConfig.getParameter<double>("jetEMfracLimit");
0081   jetMufracLimit = iConfig.getParameter<double>("jetMufracLimit");
0082   produces<METCollection>();
0083 }
0084 Type1PFMET::Type1PFMET() {}
0085
0086 // PRODUCER DESTRUCTORS -------------------------------------------
0087 Type1PFMET::~Type1PFMET() {}
0088
0089 // PRODUCER METHODS -----------------------------------------------
0090 void Type1PFMET::produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
0091   using namespace edm;
0092   Handle<PFJetCollection> inputUncorJets;
0093   iEvent.getByToken(tokenUncorJets, inputUncorJets);
0094   Handle<JetCorrector> corrector;
0095   iEvent.getByToken(correctorToken, corrector);
0096   Handle<METCollection> inputUncorMet;                         //Define Inputs
0097   iEvent.getByToken(tokenUncorMet, inputUncorMet);             //Get Inputs
0098   std::unique_ptr<METCollection> output(new METCollection());  //Create empty output
0099   run(*(inputUncorMet.product()),
0100       *(corrector.product()),
0101       *(inputUncorJets.product()),
0102       jetPTthreshold,
0103       jetEMfracLimit,
0104       jetMufracLimit,
0105       &*output);                  //Invoke the algorithm
0106   iEvent.put(std::move(output));  //Put output into Event
0107 }
0108
0109 void Type1PFMET::run(const METCollection& uncorMET,
0110                      const reco::JetCorrector& corrector,
0111                      const PFJetCollection& uncorJet,
0112                      double jetPTthreshold,
0113                      double jetEMfracLimit,
0114                      double jetMufracLimit,
0115                      METCollection* corMET) const {
0116   if (!corMET) {
0117     std::cerr << "Type1METAlgo_run-> undefined output MET collection. Stop. " << std::endl;
0118     return;
0119   }
0120
0121   double DeltaPx = 0.0;
0122   double DeltaPy = 0.0;
0123   double DeltaSumET = 0.0;
0124   // ---------------- Calculate jet corrections, but only for those uncorrected jets
0125   // ---------------- which are above the given threshold.  This requires that the
0126   // ---------------- uncorrected jets be matched with the corrected jets.
0127   for (PFJetCollection::const_iterator jet = uncorJet.begin(); jet != uncorJet.end(); ++jet) {
0128     if (jet->pt() > jetPTthreshold) {
0129       double emEFrac = jet->chargedEmEnergyFraction() + jet->neutralEmEnergyFraction();
0130       double muEFrac = jet->chargedMuEnergyFraction();
0131       if (emEFrac < jetEMfracLimit && muEFrac < jetMufracLimit) {
0132         double corr = corrector.correction(*jet) - 1.;  // correction itself
0133         DeltaPx += jet->px() * corr;
0134         DeltaPy += jet->py() * corr;
0135         DeltaSumET += jet->et() * corr;
0136       }
0137     }
0138   }
0139   //----------------- Calculate and set deltas for new MET correction
0140   CorrMETData delta;
0141   delta.mex = -DeltaPx;  //correction to MET (from Jets) is negative,
0142   delta.mey = -DeltaPy;  //since MET points in direction opposite of jets
0143   delta.sumet = DeltaSumET;
0144   //----------------- Fill holder with corrected MET (= uncorrected + delta) values
0145   const MET* u = &(uncorMET.front());
0146   double corrMetPx = u->px() + delta.mex;
0147   double corrMetPy = u->py() + delta.mey;
0148   MET::LorentzVector correctedMET4vector(corrMetPx, corrMetPy, 0., sqrt(corrMetPx * corrMetPx + corrMetPy * corrMetPy));
0149   //----------------- get previous corrections and push into new corrections
0150   std::vector<CorrMETData> corrections = u->mEtCorr();
0151   corrections.push_back(delta);
0152   //----------------- Push onto MET Collection
0153   MET result = MET(u->sumEt() + delta.sumet, corrections, correctedMET4vector, u->vertex(), u->isWeighted());
0154   corMET->push_back(result);
0155
0156   return;
0157 }