RecoAlgos/plugins/JetDeltaRValueMapProducer.cc

0001 /* \class JetDeltaRValueMapProducer
0002  *
0003  * Associates jets using delta-R matching, and writes out
0004  * a valuemap of single float variables based on a StringObjectFunction.
0005  * This is used for miniAOD.
0006  *
0007  * \author: Sal Rappoccio
0008  *
0009  *
0010  * for more details about the cut syntax, see the documentation
0011  * page below:
0012  *
0013  *   https://twiki.cern.ch/twiki/bin/view/CMS/SWGuidePhysicsCutParser
0014  *
0015  *
0016  */
0017
0018 #include "FWCore/Framework/interface/stream/EDProducer.h"
0019 #include "DataFormats/PatCandidates/interface/Jet.h"
0020 #include "DataFormats/JetReco/interface/Jet.h"
0021 #include "DataFormats/JetReco/interface/PFJet.h"
0022 #include "DataFormats/JetReco/interface/CaloJet.h"
0023
0024 #include "FWCore/Framework/interface/MakerMacros.h"
0025 #include "CommonTools/Utils/interface/StringObjectFunction.h"
0026 #include "FWCore/Framework/interface/Event.h"
0027 #include "DataFormats/Common/interface/ValueMap.h"
0028 #include "DataFormats/Math/interface/deltaR.h"
0029
0030 template <class T, class C = T>
0031 class JetDeltaRValueMapProducer : public edm::stream::EDProducer<> {
0032 public:
0033   typedef edm::ValueMap<float> JetValueMap;
0034
0035   JetDeltaRValueMapProducer(edm::ParameterSet const& params)
0036       : srcToken_(consumes<typename edm::View<T> >(params.getParameter<edm::InputTag>("src"))),
0037         matchedToken_(consumes<typename edm::View<C> >(params.getParameter<edm::InputTag>("matched"))),
0038         distMax_(params.getParameter<double>("distMax")),
0039         value_(params.existsAs<std::string>("value") ? params.getParameter<std::string>("value") : ""),
0040         values_(params.existsAs<std::vector<std::string> >("values")
0041                     ? params.getParameter<std::vector<std::string> >("values")
0042                     : std::vector<std::string>()),
0043         valueLabels_(params.existsAs<std::vector<std::string> >("valueLabels")
0044                          ? params.getParameter<std::vector<std::string> >("valueLabels")
0045                          : std::vector<std::string>()),
0046         lazyParser_(params.existsAs<bool>("lazyParser") ? params.getParameter<bool>("lazyParser") : false),
0047         multiValue_(false) {
0048     if (!value_.empty()) {
0049       evaluationMap_.insert(std::make_pair(
0050           value_, std::unique_ptr<StringObjectFunction<C> >(new StringObjectFunction<C>(value_, lazyParser_))));
0051       produces<JetValueMap>();
0052     }
0053
0054     if (!valueLabels_.empty() || !values_.empty()) {
0055       if (valueLabels_.size() == values_.size()) {
0056         multiValue_ = true;
0057         for (size_t i = 0; i < valueLabels_.size(); ++i) {
0058           evaluationMap_.insert(std::make_pair(
0059               valueLabels_[i],
0060               std::unique_ptr<StringObjectFunction<C> >(new StringObjectFunction<C>(values_[i], lazyParser_))));
0061           produces<JetValueMap>(valueLabels_[i]);
0062         }
0063       } else
0064         edm::LogWarning("ValueLabelMismatch")
0065             << "The number of value labels does not match the number of values. Values will not be evaluated.";
0066     }
0067   }
0068
0069   ~JetDeltaRValueMapProducer() override {}
0070
0071 private:
0072   void produce(edm::Event& iEvent, const edm::EventSetup& iSetup) override {
0073     edm::Handle<typename edm::View<T> > h_jets1;
0074     iEvent.getByToken(srcToken_, h_jets1);
0075     edm::Handle<typename edm::View<C> > h_jets2;
0076     iEvent.getByToken(matchedToken_, h_jets2);
0077
0078     std::vector<float> values(h_jets1->size(), -99999);
0079     std::map<std::string, std::vector<float> > valuesMap;
0080     if (multiValue_) {
0081       for (size_t i = 0; i < valueLabels_.size(); ++i)
0082         valuesMap.insert(std::make_pair(valueLabels_[i], std::vector<float>(h_jets1->size(), -99999)));
0083     }
0084     std::vector<bool> jets1_locks(h_jets1->size(), false);
0085
0086     for (typename edm::View<C>::const_iterator ibegin = h_jets2->begin(), iend = h_jets2->end(), ijet = ibegin;
0087          ijet != iend;
0088          ++ijet) {
0089       float matched_dR2 = 1e9;
0090       int matched_index = -1;
0091
0092       for (typename edm::View<T>::const_iterator jbegin = h_jets1->begin(), jend = h_jets1->end(), jjet = jbegin;
0093            jjet != jend;
0094            ++jjet) {
0095         int index = jjet - jbegin;
0096
0097         if (jets1_locks.at(index))
0098           continue;  // skip jets that have already been matched
0099
0100         float temp_dR2 = reco::deltaR2(ijet->eta(), ijet->phi(), jjet->eta(), jjet->phi());
0101         if (temp_dR2 < matched_dR2) {
0102           matched_dR2 = temp_dR2;
0103           matched_index = index;
0104         }
0105       }  // end loop over src jets
0106
0107       if (matched_index >= 0) {
0108         if (matched_dR2 > distMax_ * distMax_)
0109           edm::LogInfo("MatchedJetsFarApart") << "Matched jets separated by dR greater than distMax=" << distMax_;
0110         else {
0111           jets1_locks.at(matched_index) = true;
0112           if (!value_.empty())
0113             values.at(matched_index) = (*(evaluationMap_.at(value_)))(*ijet);
0114           if (multiValue_) {
0115             for (size_t i = 0; i < valueLabels_.size(); ++i)
0116               valuesMap.at(valueLabels_[i]).at(matched_index) = (*(evaluationMap_.at(valueLabels_[i])))(*ijet);
0117           }
0118         }
0119       }
0120     }  // end loop over matched jets
0121
0122     if (!value_.empty()) {
0123       std::unique_ptr<JetValueMap> jetValueMap(new JetValueMap());
0124
0125       JetValueMap::Filler filler(*jetValueMap);
0126       filler.insert(h_jets1, values.begin(), values.end());
0127       filler.fill();
0128
0129       // put in Event
0130       iEvent.put(std::move(jetValueMap));
0131     }
0132     if (multiValue_) {
0133       for (size_t i = 0; i < valueLabels_.size(); ++i) {
0134         std::unique_ptr<JetValueMap> jetValueMap(new JetValueMap());
0135
0136         JetValueMap::Filler filler(*jetValueMap);
0137         filler.insert(h_jets1, valuesMap.at(valueLabels_[i]).begin(), valuesMap.at(valueLabels_[i]).end());
0138         filler.fill();
0139
0140         // put in Event
0141         iEvent.put(std::move(jetValueMap), valueLabels_[i]);
0142       }
0143     }
0144   }
0145
0146   const edm::EDGetTokenT<typename edm::View<T> > srcToken_;
0147   const edm::EDGetTokenT<typename edm::View<C> > matchedToken_;
0148   const double distMax_;
0149   const std::string value_;
0150   const std::vector<std::string> values_;
0151   const std::vector<std::string> valueLabels_;
0152   const bool lazyParser_;
0153   bool multiValue_;
0154   std::map<std::string, std::unique_ptr<const StringObjectFunction<C> > > evaluationMap_;
0155 };
0156
0157 typedef JetDeltaRValueMapProducer<reco::Jet> RecoJetDeltaRValueMapProducer;
0158 typedef JetDeltaRValueMapProducer<reco::Jet, pat::Jet> RecoJetToPatJetDeltaRValueMapProducer;
0159 typedef JetDeltaRValueMapProducer<pat::Jet> PatJetDeltaRValueMapProducer;
0160
0161 DEFINE_FWK_MODULE(RecoJetDeltaRValueMapProducer);
0162 DEFINE_FWK_MODULE(RecoJetToPatJetDeltaRValueMapProducer);
0163 DEFINE_FWK_MODULE(PatJetDeltaRValueMapProducer);