Modules/interface/CorrectedJetProducer.h

0001 #ifndef CorrectedJetProducer_h
0002 #define CorrectedJetProducer_h
0003
0004 #include <sstream>
0005 #include <string>
0006 #include <vector>
0007
0008 #include "CommonTools/Utils/interface/PtComparator.h"
0009 #include "DataFormats/Common/interface/Handle.h"
0010 #include "DataFormats/Common/interface/Ref.h"
0011 #include "DataFormats/Common/interface/RefToBase.h"
0012 #include "FWCore/Framework/interface/global/EDProducer.h"
0013 #include "FWCore/Framework/interface/Event.h"
0014 #include "FWCore/Framework/interface/EventSetup.h"
0015 #include "FWCore/Utilities/interface/InputTag.h"
0016 #include "FWCore/Utilities/interface/transform.h"
0017 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0018 #include "JetMETCorrections/JetCorrector/interface/JetCorrector.h"
0019
0020 namespace edm {
0021   class ParameterSet;
0022 }
0023
0024 namespace reco {
0025   template <class T>
0026   class CorrectedJetProducer : public edm::global::EDProducer<> {
0027   public:
0028     typedef std::vector<T> JetCollection;
0029     explicit CorrectedJetProducer(const edm::ParameterSet& fParameters);
0030     ~CorrectedJetProducer() override {}
0031     void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0032
0033   private:
0034     const edm::EDGetTokenT<JetCollection> mInput;
0035     const std::vector<edm::EDGetTokenT<reco::JetCorrector> > mCorrectorTokens;
0036     const bool mVerbose;
0037   };
0038 }  // namespace reco
0039
0040 // ---------- implementation ----------
0041
0042 namespace reco {
0043
0044   template <class T>
0045   CorrectedJetProducer<T>::CorrectedJetProducer(const edm::ParameterSet& fConfig)
0046       : mInput(consumes<JetCollection>(fConfig.getParameter<edm::InputTag>("src"))),
0047         mCorrectorTokens(
0048             edm::vector_transform(fConfig.getParameter<std::vector<edm::InputTag> >("correctors"),
0049                                   [this](edm::InputTag const& tag) { return consumes<reco::JetCorrector>(tag); })),
0050         mVerbose(fConfig.getUntrackedParameter<bool>("verbose", false)) {
0051     std::string alias = fConfig.getUntrackedParameter<std::string>("alias", "");
0052     if (alias.empty())
0053       produces<JetCollection>();
0054     else
0055       produces<JetCollection>().setBranchAlias(alias);
0056   }
0057
0058   template <class T>
0059   void CorrectedJetProducer<T>::produce(edm::StreamID, edm::Event& fEvent, const edm::EventSetup& fSetup) const {
0060     // FIXME - use something more efficient instead of an std::vector
0061     std::vector<reco::JetCorrector const*> correctors(mCorrectorTokens.size(), nullptr);
0062
0063     // look for correctors
0064     for (unsigned i = 0; i < mCorrectorTokens.size(); i++) {
0065       edm::Handle<reco::JetCorrector> handle;
0066       fEvent.getByToken(mCorrectorTokens[i], handle);
0067       correctors[i] = handle.product();
0068     }
0069     edm::Handle<JetCollection> jets;                           //Define Inputs
0070     fEvent.getByToken(mInput, jets);                           //Get Inputs
0071     std::unique_ptr<JetCollection> result(new JetCollection);  //Corrected jets
0072     typename JetCollection::const_iterator jet;
0073     for (jet = jets->begin(); jet != jets->end(); jet++) {
0074       const T* referenceJet = &*jet;
0075       int index = jet - jets->begin();
0076       edm::RefToBase<reco::Jet> jetRef(edm::Ref<JetCollection>(jets, index));
0077       T correctedJet = *jet;  //copy original jet
0078       if (mVerbose)
0079         std::cout << "CorrectedJetProducer::produce-> original jet: " << jet->print() << std::endl;
0080       for (unsigned i = 0; i < mCorrectorTokens.size(); ++i) {
0081         if (!(correctors[i]->vectorialCorrection())) {
0082           // Scalar correction
0083           double scale = 1.;
0084           if (!(correctors[i]->refRequired()))
0085             scale = correctors[i]->correction(*referenceJet);
0086           else
0087             scale = correctors[i]->correction(*referenceJet, jetRef);
0088           if (mVerbose)
0089             std::cout << "CorrectedJetProducer::produce-> Corrector # " << i << ", correction factor: " << scale
0090                       << std::endl;
0091           correctedJet.scaleEnergy(scale);  // apply scalar correction
0092           referenceJet = &correctedJet;
0093         } else {
0094           // Vectorial correction
0095           reco::JetCorrector::LorentzVector corrected;
0096           double scale = correctors[i]->correction(*referenceJet, jetRef, corrected);
0097           if (mVerbose)
0098             std::cout << "CorrectedJetProducer::produce-> Corrector # " << i << ", correction factor: " << scale
0099                       << std::endl;
0100           correctedJet.setP4(corrected);  // apply vectorial correction
0101           referenceJet = &correctedJet;
0102         }
0103       }
0104       if (mVerbose)
0105         std::cout << "CorrectedJetProducer::produce-> corrected jet: " << correctedJet.print() << std::endl;
0106       result->push_back(correctedJet);
0107     }
0108     NumericSafeGreaterByPt<T> compJets;
0109     // reorder corrected jets
0110     std::sort(result->begin(), result->end(), compJets);
0111     // put corrected jet collection into event
0112     fEvent.put(std::move(result));
0113   }
0114
0115 }  // namespace reco
0116
0117 #endif