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