RecoTau/plugins/RecoTauGenericJetRegionProducer.cc

0001 /*
0002  * RecoTauJetRegionProducer
0003  *
0004  * Given a set of Jets, make new jets with the same p4 but collect all the
0005  * particle-flow candidates from a cone of a given size into the constituents.
0006  *
0007  * Author: Evan K. Friis, UC Davis
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   // Build Ptrs for all the PFCandidates
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   // Get the jets
0110   edm::Handle<reco::CandidateView> jetView;
0111   evt.getByToken(Jets_token, jetView);
0112   // Convert to a vector of JetRefs
0113   edm::RefVector<std::vector<JetType> > jets = reco::tau::castView<edm::RefVector<std::vector<JetType> > >(jetView);
0114   size_t nJets = jets.size();
0115
0116   // Get the association map matching jets to Candidates
0117   // (needed for reconstruction of boosted taus)
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       // Get a ref to jet
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   // Get the original product, so we can match against it - otherwise the
0134   // indices don't match up.
0135   edm::ProductID originalId = jets.id();
0136   edm::Handle<std::vector<JetType> > originalJets;
0137   size_t nOriginalJets = 0;
0138   // We have to make sure that we have some selected jets, otherwise we don't
0139   // actually have a valid product ID to the original jets.
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   // Keep track of the indices of the current jet and the old (original) jet
0154   // -1 indicates no match.
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     // Get a ref to jet
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     // Make an initial copy.
0166     newJets->emplace_back(*jetRef);
0167     JetType& newJet = newJets->back();
0168     // Clear out all the constituents
0169     newJet.clearDaughters();
0170     // Loop over all the PFCands
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     // Match the index of the jet we just made to the index into the original
0199     // collection.
0200     //matchInfo[jetRef.key()] = ijet;
0201     matchInfo[jetRef.key()] = nNewJets;
0202     nNewJets++;
0203   }
0204
0205   // Put our new jets into the event
0206   edm::OrphanHandle<std::vector<JetType> > newJetsInEvent = evt.put(std::move(newJets), "jets");
0207
0208   // Create a matching between original jets -> extra collection
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   // RecoTauGenericJetRegionProducer
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   // RecoTauGenericJetRegionProducer
0239   RecoTauGenericJetRegionProducer::fillDescriptionsBase(descriptions, "RecoTauJetRegionProducer");
0240 }
0241
0242 template <>
0243 void RecoTauGenericJetRegionProducer<pat::Jet, pat::PackedCandidate>::fillDescriptions(
0244     edm::ConfigurationDescriptions& descriptions) {
0245   // RecoTauGenericJetRegionProducer
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);