File indexing completed on 2023-03-17 11:21:54
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011 #include "DataFormats/JetReco/interface/PFJet.h"
0012 #include "DataFormats/PatCandidates/interface/Jet.h"
0013 #include "DataFormats/Common/interface/Association.h"
0014 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0015 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidateFwd.h"
0016 #include "DataFormats/PatCandidates/interface/PackedCandidate.h"
0017 #include "DataFormats/Common/interface/AssociationMap.h"
0018
0019 #include "RecoTauTag/RecoTau/interface/ConeTools.h"
0020 #include "RecoTauTag/RecoTau/interface/RecoTauCommonUtilities.h"
0021
0022 #include "FWCore/Framework/interface/stream/EDProducer.h"
0023 #include "FWCore/Framework/interface/EventSetup.h"
0024 #include "FWCore/Framework/interface/Event.h"
0025 #include "FWCore/Framework/interface/makeRefToBaseProdFrom.h"
0026 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0027 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0028
0029 #include <FWCore/ParameterSet/interface/ConfigurationDescriptions.h>
0030 #include <FWCore/ParameterSet/interface/ParameterSetDescription.h>
0031
0032 #include <string>
0033 #include <iostream>
0034
0035 template <class JetType, class CandType>
0036 class RecoTauGenericJetRegionProducer : public edm::stream::EDProducer<> {
0037 public:
0038 typedef edm::AssociationMap<edm::OneToOne<reco::JetView, reco::JetView> > JetMatchMap;
0039 typedef edm::AssociationMap<edm::OneToMany<std::vector<JetType>, std::vector<CandType>, unsigned int> >
0040 JetToCandidateAssociation;
0041 explicit RecoTauGenericJetRegionProducer(const edm::ParameterSet& pset);
0042 ~RecoTauGenericJetRegionProducer() override {}
0043
0044 void produce(edm::Event& evt, const edm::EventSetup& es) override;
0045
0046 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0047 static void fillDescriptionsBase(edm::ConfigurationDescriptions& descriptions, const std::string& name);
0048
0049 private:
0050 std::string moduleLabel_;
0051
0052 edm::InputTag inputJets_;
0053 edm::InputTag pfCandSrc_;
0054 edm::InputTag pfCandAssocMapSrc_;
0055
0056 edm::EDGetTokenT<std::vector<CandType> > pf_token;
0057 edm::EDGetTokenT<reco::CandidateView> Jets_token;
0058 edm::EDGetTokenT<JetToCandidateAssociation> pfCandAssocMap_token;
0059
0060 double minJetPt_;
0061 double maxJetAbsEta_;
0062 double deltaR2_;
0063
0064 int verbosity_;
0065 };
0066
0067 template <class JetType, class CandType>
0068 RecoTauGenericJetRegionProducer<JetType, CandType>::RecoTauGenericJetRegionProducer(const edm::ParameterSet& cfg)
0069 : moduleLabel_(cfg.getParameter<std::string>("@module_label")) {
0070 inputJets_ = cfg.getParameter<edm::InputTag>("src");
0071 pfCandSrc_ = cfg.getParameter<edm::InputTag>("pfCandSrc");
0072 pfCandAssocMapSrc_ = cfg.getParameter<edm::InputTag>("pfCandAssocMapSrc");
0073
0074 pf_token = consumes<std::vector<CandType> >(pfCandSrc_);
0075 Jets_token = consumes<reco::CandidateView>(inputJets_);
0076 pfCandAssocMap_token = consumes<JetToCandidateAssociation>(pfCandAssocMapSrc_);
0077
0078 double deltaR = cfg.getParameter<double>("deltaR");
0079 deltaR2_ = deltaR * deltaR;
0080 minJetPt_ = cfg.getParameter<double>("minJetPt");
0081 maxJetAbsEta_ = cfg.getParameter<double>("maxJetAbsEta");
0082
0083 verbosity_ = cfg.getParameter<int>("verbosity");
0084
0085 produces<std::vector<JetType> >("jets");
0086 produces<JetMatchMap>();
0087 }
0088
0089 template <class JetType, class CandType>
0090 void RecoTauGenericJetRegionProducer<JetType, CandType>::produce(edm::Event& evt, const edm::EventSetup& es) {
0091 if (verbosity_) {
0092 std::cout << "<RecoTauJetRegionProducer::produce (moduleLabel = " << moduleLabel_ << ")>:" << std::endl;
0093 std::cout << " inputJets = " << inputJets_ << std::endl;
0094 std::cout << " pfCandSrc = " << pfCandSrc_ << std::endl;
0095 std::cout << " pfCandAssocMapSrc_ = " << pfCandAssocMapSrc_ << std::endl;
0096 }
0097
0098 edm::Handle<std::vector<CandType> > pfCandsHandle;
0099 evt.getByToken(pf_token, pfCandsHandle);
0100
0101
0102 typedef edm::Ptr<CandType> CandPtr;
0103 std::vector<CandPtr> pfCands;
0104 pfCands.reserve(pfCandsHandle->size());
0105 for (size_t icand = 0; icand < pfCandsHandle->size(); ++icand) {
0106 pfCands.push_back(CandPtr(pfCandsHandle, icand));
0107 }
0108
0109
0110 edm::Handle<reco::CandidateView> jetView;
0111 evt.getByToken(Jets_token, jetView);
0112
0113 edm::RefVector<std::vector<JetType> > jets = reco::tau::castView<edm::RefVector<std::vector<JetType> > >(jetView);
0114 size_t nJets = jets.size();
0115
0116
0117
0118 edm::Handle<JetToCandidateAssociation> jetToPFCandMap;
0119 std::vector<std::unordered_set<unsigned> > fastJetToPFCandMap;
0120 if (!pfCandAssocMapSrc_.label().empty()) {
0121 evt.getByToken(pfCandAssocMap_token, jetToPFCandMap);
0122 fastJetToPFCandMap.resize(nJets);
0123 for (size_t ijet = 0; ijet < nJets; ++ijet) {
0124
0125 const edm::Ref<std::vector<JetType> >& jetRef = jets[ijet];
0126 const auto& pfCandsMappedToJet = (*jetToPFCandMap)[jetRef];
0127 for (const auto& pfCandMappedToJet : pfCandsMappedToJet) {
0128 fastJetToPFCandMap[ijet].emplace(pfCandMappedToJet.key());
0129 }
0130 }
0131 }
0132
0133
0134
0135 edm::ProductID originalId = jets.id();
0136 edm::Handle<std::vector<JetType> > originalJets;
0137 size_t nOriginalJets = 0;
0138
0139
0140 if (nJets) {
0141 try {
0142 evt.get(originalId, originalJets);
0143 } catch (const cms::Exception& e) {
0144 edm::LogError("MissingOriginalCollection") << "Can't get the original jets that made: " << inputJets_
0145 << " that have product ID: " << originalId << " from the event!!";
0146 throw e;
0147 }
0148 nOriginalJets = originalJets->size();
0149 }
0150
0151 auto newJets = std::make_unique<std::vector<JetType> >();
0152
0153
0154
0155 std::vector<int> matchInfo(nOriginalJets, -1);
0156 newJets->reserve(nJets);
0157 size_t nNewJets = 0;
0158 for (size_t ijet = 0; ijet < nJets; ++ijet) {
0159
0160 const edm::Ref<std::vector<JetType> >& jetRef = jets[ijet];
0161 if (jetRef->pt() - minJetPt_ < 1e-5)
0162 continue;
0163 if (std::abs(jetRef->eta()) - maxJetAbsEta_ > -1e-5)
0164 continue;
0165
0166 newJets->emplace_back(*jetRef);
0167 JetType& newJet = newJets->back();
0168
0169 newJet.clearDaughters();
0170
0171 for (const auto& pfCand : pfCands) {
0172 bool isMappedToJet = false;
0173 if (jetToPFCandMap.isValid()) {
0174 auto temp = jetToPFCandMap->find(jetRef);
0175 if (temp == jetToPFCandMap->end()) {
0176 edm::LogWarning("WeirdCandidateMap") << "Candidate map for jet " << jetRef.key() << " is empty!";
0177 continue;
0178 }
0179 isMappedToJet = fastJetToPFCandMap[ijet].count(pfCand.key());
0180 } else {
0181 isMappedToJet = true;
0182 }
0183 if (reco::deltaR2(*jetRef, *pfCand) < deltaR2_ && isMappedToJet)
0184 newJet.addDaughter(pfCand);
0185 }
0186 if (verbosity_) {
0187 std::cout << "jet #" << ijet << ": Pt = " << jetRef->pt() << ", eta = " << jetRef->eta()
0188 << ", phi = " << jetRef->eta() << ","
0189 << " mass = " << jetRef->mass() << ", area = " << jetRef->jetArea() << std::endl;
0190 auto jetConstituents = newJet.daughterPtrVector();
0191 int idx = 0;
0192 for (const auto& jetConstituent : jetConstituents) {
0193 std::cout << " constituent #" << idx << ": Pt = " << jetConstituent->pt() << ", eta = " << jetConstituent->eta()
0194 << ", phi = " << jetConstituent->phi() << std::endl;
0195 ++idx;
0196 }
0197 }
0198
0199
0200
0201 matchInfo[jetRef.key()] = nNewJets;
0202 nNewJets++;
0203 }
0204
0205
0206 edm::OrphanHandle<std::vector<JetType> > newJetsInEvent = evt.put(std::move(newJets), "jets");
0207
0208
0209 auto matching = (nJets != 0)
0210 ? std::make_unique<JetMatchMap>(
0211 edm::makeRefToBaseProdFrom(edm::RefToBase<reco::Jet>(jets[0]), evt), newJetsInEvent)
0212 : std::make_unique<JetMatchMap>();
0213 for (size_t ijet = 0; ijet < nJets; ++ijet) {
0214 matching->insert(edm::RefToBase<reco::Jet>(jets[ijet]),
0215 edm::RefToBase<reco::Jet>(edm::Ref<std::vector<JetType> >(newJetsInEvent, matchInfo[ijet])));
0216 }
0217 evt.put(std::move(matching));
0218 }
0219
0220 template <class JetType, class CandType>
0221 void RecoTauGenericJetRegionProducer<JetType, CandType>::fillDescriptionsBase(
0222 edm::ConfigurationDescriptions& descriptions, const std::string& name) {
0223
0224 edm::ParameterSetDescription desc;
0225 desc.add<edm::InputTag>("src", edm::InputTag("ak4PFJets"));
0226 desc.add<double>("deltaR", 0.8);
0227 desc.add<edm::InputTag>("pfCandAssocMapSrc", edm::InputTag(""));
0228 desc.add<int>("verbosity", 0);
0229 desc.add<double>("maxJetAbsEta", 2.5);
0230 desc.add<double>("minJetPt", 14.0);
0231 desc.add<edm::InputTag>("pfCandSrc", edm::InputTag("particleFlow"));
0232 descriptions.add(name, desc);
0233 }
0234
0235 template <>
0236 void RecoTauGenericJetRegionProducer<reco::PFJet, reco::PFCandidate>::fillDescriptions(
0237 edm::ConfigurationDescriptions& descriptions) {
0238
0239 RecoTauGenericJetRegionProducer::fillDescriptionsBase(descriptions, "RecoTauJetRegionProducer");
0240 }
0241
0242 template <>
0243 void RecoTauGenericJetRegionProducer<pat::Jet, pat::PackedCandidate>::fillDescriptions(
0244 edm::ConfigurationDescriptions& descriptions) {
0245
0246 RecoTauGenericJetRegionProducer::fillDescriptionsBase(descriptions, "RecoTauPatJetRegionProducer");
0247 }
0248
0249 typedef RecoTauGenericJetRegionProducer<reco::PFJet, reco::PFCandidate> RecoTauJetRegionProducer;
0250 typedef RecoTauGenericJetRegionProducer<pat::Jet, pat::PackedCandidate> RecoTauPatJetRegionProducer;
0251
0252 #include "FWCore/Framework/interface/MakerMacros.h"
0253 DEFINE_FWK_MODULE(RecoTauJetRegionProducer);
0254 DEFINE_FWK_MODULE(RecoTauPatJetRegionProducer);