File indexing completed on 2024-04-06 12:01:09
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
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;
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 }
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 }
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
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
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);