File indexing completed on 2023-03-17 11:10:49
0001
0002
0003
0004
0005
0006
0007
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
0058
0059
0060 edm::ParameterSetDescription desc;
0061 desc.add<edm::InputTag>("inputUncorJetsTag", edm::InputTag("ak4PFJets"));
0062 desc.add<double>("jetEMfracLimit", 0.95);
0063 desc.add<double>("jetMufracLimit", 0.95);
0064 desc.add<std::string>("metType", "PFMET");
0065 desc.add<double>("jetPTthreshold", 20.0);
0066
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
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
0087 Type1PFMET::~Type1PFMET() {}
0088
0089
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;
0097 iEvent.getByToken(tokenUncorMet, inputUncorMet);
0098 std::unique_ptr<METCollection> output(new METCollection());
0099 run(*(inputUncorMet.product()),
0100 *(corrector.product()),
0101 *(inputUncorJets.product()),
0102 jetPTthreshold,
0103 jetEMfracLimit,
0104 jetMufracLimit,
0105 &*output);
0106 iEvent.put(std::move(output));
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
0125
0126
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.;
0133 DeltaPx += jet->px() * corr;
0134 DeltaPy += jet->py() * corr;
0135 DeltaSumET += jet->et() * corr;
0136 }
0137 }
0138 }
0139
0140 CorrMETData delta;
0141 delta.mex = -DeltaPx;
0142 delta.mey = -DeltaPy;
0143 delta.sumet = DeltaSumET;
0144
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
0150 std::vector<CorrMETData> corrections = u->mEtCorr();
0151 corrections.push_back(delta);
0152
0153 MET result = MET(u->sumEt() + delta.sumet, corrections, correctedMET4vector, u->vertex(), u->isWeighted());
0154 corMET->push_back(result);
0155
0156 return;
0157 }