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, typename TN = float>
0031 class JetDeltaRValueMapProducer : public edm::stream::EDProducer<> {
0032 public:
0033   typedef TN value_type;
0034   typedef edm::ValueMap<value_type> JetValueMap;
0035   typedef std::vector<value_type> ValueVector;
0036
0037   JetDeltaRValueMapProducer(edm::ParameterSet const& params)
0038       : srcToken_(consumes<typename edm::View<T> >(params.getParameter<edm::InputTag>("src"))),
0039         matchedToken_(consumes<typename edm::View<C> >(params.getParameter<edm::InputTag>("matched"))),
0040         distMax_(params.getParameter<double>("distMax")),
0041         value_(params.existsAs<std::string>("value") ? params.getParameter<std::string>("value") : ""),
0042         values_(params.existsAs<std::vector<std::string> >("values")
0043                     ? params.getParameter<std::vector<std::string> >("values")
0044                     : std::vector<std::string>()),
0045         valueLabels_(params.existsAs<std::vector<std::string> >("valueLabels")
0046                          ? params.getParameter<std::vector<std::string> >("valueLabels")
0047                          : std::vector<std::string>()),
0048         lazyParser_(params.existsAs<bool>("lazyParser") ? params.getParameter<bool>("lazyParser") : false),
0049         multiValue_(false) {
0050     if (!value_.empty()) {
0051       if (value_ != "index") {
0052         evaluationMap_.insert(std::make_pair(
0053             value_, std::unique_ptr<StringObjectFunction<C> >(new StringObjectFunction<C>(value_, lazyParser_))));
0054       }
0055       produces<JetValueMap>();
0056     }
0057
0058     if (!valueLabels_.empty() || !values_.empty()) {
0059       if (valueLabels_.size() == values_.size()) {
0060         multiValue_ = true;
0061         for (size_t i = 0; i < valueLabels_.size(); ++i) {
0062           evaluationMap_.insert(std::make_pair(
0063               valueLabels_[i],
0064               std::unique_ptr<StringObjectFunction<C> >(new StringObjectFunction<C>(values_[i], lazyParser_))));
0065           produces<JetValueMap>(valueLabels_[i]);
0066         }
0067       } else
0068         edm::LogWarning("ValueLabelMismatch")
0069             << "The number of value labels does not match the number of values. Values will not be evaluated.";
0070     }
0071   }
0072
0073   ~JetDeltaRValueMapProducer() override {}
0074
0075 private:
0076   void produce(edm::Event& iEvent, const edm::EventSetup& iSetup) override {
0077     edm::Handle<typename edm::View<T> > h_jets1;
0078     iEvent.getByToken(srcToken_, h_jets1);
0079     edm::Handle<typename edm::View<C> > h_jets2;
0080     iEvent.getByToken(matchedToken_, h_jets2);
0081
0082     ValueVector values(h_jets1->size(), -99999);
0083
0084     std::map<std::string, std::vector<float> > valuesMap;
0085     if (multiValue_) {
0086       for (size_t i = 0; i < valueLabels_.size(); ++i)
0087         valuesMap.insert(std::make_pair(valueLabels_[i], std::vector<float>(h_jets1->size(), -99999)));
0088     }
0089     std::vector<bool> jets1_locks(h_jets1->size(), false);
0090
0091     for (typename edm::View<C>::const_iterator ibegin = h_jets2->begin(), iend = h_jets2->end(), ijet = ibegin;
0092          ijet != iend;
0093          ++ijet) {
0094       float matched_dR2 = 1e9;
0095       int matched_index = -1;
0096       int iindex = ijet - ibegin;
0097
0098       for (typename edm::View<T>::const_iterator jbegin = h_jets1->begin(), jend = h_jets1->end(), jjet = jbegin;
0099            jjet != jend;
0100            ++jjet) {
0101         int jindex = jjet - jbegin;
0102
0103         if (jets1_locks.at(jindex))
0104           continue;  // skip jets that have already been matched
0105
0106         float temp_dR2 = reco::deltaR2(ijet->eta(), ijet->phi(), jjet->eta(), jjet->phi());
0107         if (temp_dR2 < matched_dR2) {
0108           matched_dR2 = temp_dR2;
0109           matched_index = jindex;
0110         }
0111       }  // end loop over src jets
0112
0113       if (matched_index >= 0) {
0114         if (matched_dR2 > distMax_ * distMax_)
0115           edm::LogInfo("MatchedJetsFarApart") << "Matched jets separated by dR greater than distMax=" << distMax_;
0116         else {
0117           jets1_locks.at(matched_index) = true;
0118           if (!value_.empty()) {
0119             if (value_ == "index") {
0120               values.at(matched_index) = iindex;
0121             } else {
0122               values.at(matched_index) = (*(evaluationMap_.at(value_)))(*ijet);
0123             }
0124           }
0125           if (multiValue_) {
0126             for (size_t i = 0; i < valueLabels_.size(); ++i)
0127               valuesMap.at(valueLabels_[i]).at(matched_index) = (*(evaluationMap_.at(valueLabels_[i])))(*ijet);
0128           }
0129         }
0130       }
0131     }  // end loop over matched jets
0132
0133     if (!value_.empty()) {
0134       std::unique_ptr<JetValueMap> jetValueMap(new JetValueMap());
0135
0136       typename JetValueMap::Filler filler(*jetValueMap);
0137       filler.insert(h_jets1, values.begin(), values.end());
0138       filler.fill();
0139
0140       // put in Event
0141       iEvent.put(std::move(jetValueMap));
0142     }
0143     if (multiValue_) {
0144       for (size_t i = 0; i < valueLabels_.size(); ++i) {
0145         std::unique_ptr<JetValueMap> jetValueMap(new JetValueMap());
0146
0147         typename JetValueMap::Filler filler(*jetValueMap);
0148         filler.insert(h_jets1, valuesMap.at(valueLabels_[i]).begin(), valuesMap.at(valueLabels_[i]).end());
0149         filler.fill();
0150
0151         // put in Event
0152         iEvent.put(std::move(jetValueMap), valueLabels_[i]);
0153       }
0154     }
0155   }
0156
0157   const edm::EDGetTokenT<typename edm::View<T> > srcToken_;
0158   const edm::EDGetTokenT<typename edm::View<C> > matchedToken_;
0159   const double distMax_;
0160   const std::string value_;
0161   const std::vector<std::string> values_;
0162   const std::vector<std::string> valueLabels_;
0163   const bool lazyParser_;
0164   bool multiValue_;
0165   std::map<std::string, std::unique_ptr<const StringObjectFunction<C> > > evaluationMap_;
0166 };
0167
0168 typedef JetDeltaRValueMapProducer<reco::Jet> RecoJetDeltaRValueMapProducer;
0169 typedef JetDeltaRValueMapProducer<reco::Jet, pat::Jet> RecoJetToPatJetDeltaRValueMapProducer;
0170 typedef JetDeltaRValueMapProducer<pat::Jet> PatJetDeltaRValueMapProducer;
0171 typedef JetDeltaRValueMapProducer<reco::Jet, reco::GenJet, int> RecoJetToGenJetDeltaRValueMapProducer;
0172
0173 DEFINE_FWK_MODULE(RecoJetDeltaRValueMapProducer);
0174 DEFINE_FWK_MODULE(RecoJetToPatJetDeltaRValueMapProducer);
0175 DEFINE_FWK_MODULE(PatJetDeltaRValueMapProducer);
0176 DEFINE_FWK_MODULE(RecoJetToGenJetDeltaRValueMapProducer);