JetProducers/plugins/BoostedTauSeedsProducer.cc

0001
0002 /* =============================================================================
0003  *       Filename:  BoostedTauSeedsProducer.cc
0004  *
0005  *    Description:  Take the two subjets found by CMSBoostedTauSeedingAlgorithm
0006  *                  and add the data-formats for
0007  *                 o seeding the reconstruction of hadronic taus (PFJets, collection of PFCandidates not within subjet)
0008  *                 o computation of electron and muon isolation Pt-sums, excluding the particles in the other subjet (collection of PFCandidates within subjet, used to define "Vetos")
0009  *
0010  *        Created:  10/22/2013 16:05:00
0011  *
0012  *         Authors:  Christian Veelken (LLR)
0013  *
0014  * =============================================================================
0015  */
0016
0017 #include "FWCore/Framework/interface/Frameworkfwd.h"
0018 #include "FWCore/Framework/interface/EventSetup.h"
0019 #include "FWCore/Framework/interface/Event.h"
0020 #include "FWCore/Framework/interface/stream/EDProducer.h"
0021 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0022 #include "FWCore/Utilities/interface/InputTag.h"
0023 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0024
0025 #include "DataFormats/JetReco/interface/Jet.h"
0026 #include "DataFormats/Common/interface/View.h"
0027 #include "DataFormats/JetReco/interface/PFJet.h"
0028 #include "DataFormats/JetReco/interface/PFJetCollection.h"
0029 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0030 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidateFwd.h"
0031 #include "DataFormats/PatCandidates/interface/Jet.h"
0032 #include "DataFormats/PatCandidates/interface/PackedCandidate.h"
0033 #include "DataFormats/Common/interface/AssociationMap.h"
0034 #include "DataFormats/Common/interface/Ref.h"
0035 #include "DataFormats/Common/interface/RefProd.h"
0036 #include "DataFormats/Common/interface/Handle.h"
0037 #include "DataFormats/Math/interface/deltaR.h"
0038
0039 #include <TMath.h>
0040
0041 #include <string>
0042 #include <iostream>
0043 #include <iomanip>
0044
0045 template <class JetType, class CandType>
0046 class GenericBoostedTauSeedsProducer : public edm::stream::EDProducer<> {
0047 public:
0048   explicit GenericBoostedTauSeedsProducer(const edm::ParameterSet&);
0049   ~GenericBoostedTauSeedsProducer() override {}
0050   void produce(edm::Event&, const edm::EventSetup&) override;
0051
0052 private:
0053   typedef edm::AssociationMap<edm::OneToMany<std::vector<JetType>, std::vector<CandType>, unsigned int>>
0054       JetToPFCandidateAssociation;
0055   typedef std::vector<JetType> JetTypeCollection;
0056   typedef std::vector<CandType> CandTypeCollection;
0057
0058   std::string moduleLabel_;
0059
0060   typedef edm::View<reco::Jet> JetView;
0061   edm::EDGetTokenT<JetView> srcSubjets_;
0062   edm::EDGetTokenT<CandTypeCollection> srcPFCandidates_;
0063
0064   int verbosity_;
0065 };
0066
0067 template <class JetType, class CandType>
0068 GenericBoostedTauSeedsProducer<JetType, CandType>::GenericBoostedTauSeedsProducer(const edm::ParameterSet& cfg)
0069     : moduleLabel_(cfg.getParameter<std::string>("@module_label")) {
0070   srcSubjets_ = consumes<JetView>(cfg.getParameter<edm::InputTag>("subjetSrc"));
0071   srcPFCandidates_ = consumes<CandTypeCollection>(cfg.getParameter<edm::InputTag>("pfCandidateSrc"));
0072
0073   verbosity_ = (cfg.exists("verbosity")) ? cfg.getParameter<int>("verbosity") : 0;
0074
0075   produces<JetTypeCollection>();
0076   produces<JetToPFCandidateAssociation>("pfCandAssocMapForIsolation");
0077   //produces<JetToPFCandidateAssociation>("pfCandAssocMapForIsoDepositVetos");
0078 }
0079
0080 namespace {
0081   typedef std::vector<std::unordered_set<uint32_t>> JetToConstitMap;
0082
0083   template <class JetType, class CandType>
0084   JetType convertToPFJet(const reco::Jet& jet, const reco::Jet::Constituents& jetConstituents) {
0085     // CV: code for filling pfJetSpecific objects taken from
0086     //        RecoParticleFlow/PFRootEvent/src/JetMaker.cc
0087     double chargedHadronEnergy = 0.;
0088     double neutralHadronEnergy = 0.;
0089     double chargedEmEnergy = 0.;
0090     double neutralEmEnergy = 0.;
0091     double chargedMuEnergy = 0.;
0092     int chargedMultiplicity = 0;
0093     int neutralMultiplicity = 0;
0094     int muonMultiplicity = 0;
0095     for (auto const& jetConstituent : jetConstituents) {
0096       const CandType* pfCandidate = dynamic_cast<const CandType*>(jetConstituent.get());
0097       if (pfCandidate) {
0098         switch (std::abs(pfCandidate->pdgId())) {
0099           case 211:  // charged hadron
0100             chargedHadronEnergy += pfCandidate->energy();
0101             ++chargedMultiplicity;
0102             break;
0103           case 11:  // electron
0104             chargedEmEnergy += pfCandidate->energy();
0105             ++chargedMultiplicity;
0106             break;
0107           case 13:  // muon
0108             chargedMuEnergy += pfCandidate->energy();
0109             ++chargedMultiplicity;
0110             ++muonMultiplicity;
0111             break;
0112           case 22:  // photon
0113           case 2:   // electromagnetic in HF
0114             neutralEmEnergy += pfCandidate->energy();
0115             ++neutralMultiplicity;
0116             break;
0117           case 130:  // neutral hadron
0118           case 1:    // hadron in HF
0119             neutralHadronEnergy += pfCandidate->energy();
0120             ++neutralMultiplicity;
0121             break;
0122           default:
0123             edm::LogWarning("convertToPFJet") << "PFCandidate: Pt = " << pfCandidate->pt()
0124                                               << ", eta = " << pfCandidate->eta() << ", phi = " << pfCandidate->phi()
0125                                               << " has invalid pdgID = " << pfCandidate->pdgId() << " !!" << std::endl;
0126             break;
0127         }
0128       } else {
0129         edm::LogWarning("convertToPFJet")
0130             << "Jet constituent: Pt = " << jetConstituent->pt() << ", eta = " << jetConstituent->eta()
0131             << ", phi = " << jetConstituent->phi() << " is not of type PFCandidate !!" << std::endl;
0132       }
0133     }
0134     reco::PFJet::Specific pfJetSpecific;
0135     pfJetSpecific.mChargedHadronEnergy = chargedHadronEnergy;
0136     pfJetSpecific.mNeutralHadronEnergy = neutralHadronEnergy;
0137     pfJetSpecific.mChargedEmEnergy = chargedEmEnergy;
0138     pfJetSpecific.mChargedMuEnergy = chargedMuEnergy;
0139     pfJetSpecific.mNeutralEmEnergy = neutralEmEnergy;
0140     pfJetSpecific.mChargedMultiplicity = chargedMultiplicity;
0141     pfJetSpecific.mNeutralMultiplicity = neutralMultiplicity;
0142     pfJetSpecific.mMuonMultiplicity = muonMultiplicity;
0143
0144     reco::PFJet pfJet(jet.p4(), jet.vertex(), pfJetSpecific, jetConstituents);
0145     pfJet.setJetArea(jet.jetArea());
0146
0147     return JetType(pfJet);
0148   }
0149
0150   void getJetConstituents(const reco::Jet& jet, reco::Jet::Constituents& jet_and_subjetConstituents) {
0151     reco::Jet::Constituents jetConstituents = jet.getJetConstituents();
0152     for (auto const& jetConstituent : jetConstituents) {
0153       const reco::Jet* subjet = dynamic_cast<const reco::Jet*>(jetConstituent.get());
0154       if (subjet) {
0155         getJetConstituents(*subjet, jet_and_subjetConstituents);
0156       } else {
0157         jet_and_subjetConstituents.push_back(jetConstituent);
0158       }
0159     }
0160   }
0161
0162   template <class CandType>
0163   std::vector<edm::Ref<std::vector<CandType>>> getPFCandidates_exclJetConstituents(
0164       const reco::Jet& jet,
0165       const edm::Handle<std::vector<CandType>>& pfCandidates,
0166       const JetToConstitMap::value_type& constitmap,
0167       bool invert) {
0168     auto const& collection_cand = (*pfCandidates);
0169     std::vector<edm::Ref<std::vector<CandType>>> pfCandidates_exclJetConstituents;
0170     size_t numPFCandidates = pfCandidates->size();
0171     for (size_t pfCandidateIdx = 0; pfCandidateIdx < numPFCandidates; ++pfCandidateIdx) {
0172       if (!(deltaR2(collection_cand[pfCandidateIdx], jet) < 1.0))
0173         continue;
0174       bool isJetConstituent = constitmap.count(pfCandidateIdx);
0175       if (!(isJetConstituent ^ invert)) {
0176         edm::Ref<std::vector<CandType>> pfCandidate(pfCandidates, pfCandidateIdx);
0177         pfCandidates_exclJetConstituents.push_back(pfCandidate);
0178       }
0179     }
0180     return pfCandidates_exclJetConstituents;
0181   }
0182
0183   void printJetConstituents(const std::string& label, const reco::Jet::Constituents& jetConstituents) {
0184     std::cout << "#" << label << " = " << jetConstituents.size() << ":" << std::endl;
0185     unsigned idx = 0;
0186     for (auto const& jetConstituent : jetConstituents) {
0187       std::cout << " jetConstituent #" << idx << ": Pt = " << (*jetConstituent).pt()
0188                 << ", eta = " << (*jetConstituent).eta() << ", phi = " << (*jetConstituent).phi() << std::endl;
0189       ++idx;
0190     }
0191   }
0192 }  // namespace
0193
0194 template <class JetType, class CandType>
0195 void GenericBoostedTauSeedsProducer<JetType, CandType>::produce(edm::Event& evt, const edm::EventSetup& es) {
0196   if (verbosity_ >= 1) {
0197     std::cout << "<BoostedTauSeedsProducer::produce (moduleLabel = " << moduleLabel_ << ")>:" << std::endl;
0198   }
0199
0200   edm::Handle<JetView> subjets;
0201   evt.getByToken(srcSubjets_, subjets);
0202   if (verbosity_ >= 1) {
0203     std::cout << "#subjets = " << subjets->size() << std::endl;
0204   }
0205   assert((subjets->size() % 2) == 0);  // CV: ensure that subjets come in pairs
0206
0207   edm::Handle<CandTypeCollection> pfCandidates;
0208   evt.getByToken(srcPFCandidates_, pfCandidates);
0209   if (verbosity_ >= 1) {
0210     std::cout << "#pfCandidates = " << pfCandidates->size() << std::endl;
0211   }
0212
0213   auto selectedSubjets = std::make_unique<JetTypeCollection>();
0214   edm::RefProd<JetTypeCollection> selectedSubjetRefProd = evt.getRefBeforePut<JetTypeCollection>();
0215
0216   auto selectedSubjetPFCandidateAssociationForIsolation =
0217       std::make_unique<JetToPFCandidateAssociation>(&evt.productGetter());
0218   //auto selectedSubjetPFCandidateAssociationForIsoDepositVetos = std::make_unique<JetToPFCandidateAssociation>(&evt.productGetter());
0219
0220   // cache for jet->pfcandidate
0221   JetToConstitMap constitmap(subjets->size());
0222
0223   // fill constituents map
0224   const auto& thesubjets = *subjets;
0225   for (unsigned i = 0; i < thesubjets.size(); ++i) {
0226     for (unsigned j = 0; j < thesubjets[i].numberOfDaughters(); ++j) {
0227       constitmap[i].emplace(thesubjets[i].daughterPtr(j).key());
0228     }
0229   }
0230
0231   for (size_t idx = 0; idx < (subjets->size() / 2); ++idx) {
0232     const reco::Jet* subjet1 = &subjets->at(2 * idx);
0233     const reco::Jet* subjet2 = &subjets->at(2 * idx + 1);
0234     assert(subjet1 && subjet2);
0235     if (verbosity_ >= 1) {
0236       std::cout << "processing jet #" << idx << ":" << std::endl;
0237       std::cout << " subjet1: Pt = " << subjet1->pt() << ", eta = " << subjet1->eta() << ", phi = " << subjet1->phi()
0238                 << ", mass = " << subjet1->mass() << " (#constituents = " << subjet1->nConstituents()
0239                 << ", area = " << subjet1->jetArea() << ")" << std::endl;
0240       std::cout << " subjet2: Pt = " << subjet2->pt() << ", eta = " << subjet2->eta() << ", phi = " << subjet2->phi()
0241                 << ", mass = " << subjet2->mass() << " (#constituents = " << subjet2->nConstituents()
0242                 << ", area = " << subjet2->jetArea() << ")" << std::endl;
0243     }
0244
0245     if (!(subjet1->nConstituents() >= 1 && subjet1->pt() > 1. && subjet2->nConstituents() >= 1 && subjet2->pt() > 1.))
0246       continue;  // CV: skip pathological cases
0247
0248     // find PFCandidate constituents of each subjet
0249     reco::Jet::Constituents subjetConstituents1;
0250     getJetConstituents(*subjet1, subjetConstituents1);
0251     reco::Jet::Constituents subjetConstituents2;
0252     getJetConstituents(*subjet2, subjetConstituents2);
0253     if (verbosity_ >= 1) {
0254       printJetConstituents("subjetConstituents1", subjetConstituents1);
0255       printJetConstituents("subjetConstituents2", subjetConstituents2);
0256     }
0257
0258     selectedSubjets->push_back(convertToPFJet<JetType, CandType>(*subjet1, subjetConstituents1));
0259     edm::Ref<JetTypeCollection> subjetRef1(selectedSubjetRefProd, selectedSubjets->size() - 1);
0260     selectedSubjets->push_back(convertToPFJet<JetType, CandType>(*subjet2, subjetConstituents2));
0261     edm::Ref<JetTypeCollection> subjetRef2(selectedSubjetRefProd, selectedSubjets->size() - 1);
0262
0263     // find all PFCandidates that are not constituents of the **other** subjet
0264     std::vector<edm::Ref<std::vector<CandType>>> pfCandidatesNotInSubjet1 =
0265         getPFCandidates_exclJetConstituents<CandType>(*subjet1, pfCandidates, constitmap[2 * idx + 1], false);
0266     std::vector<edm::Ref<std::vector<CandType>>> pfCandidatesNotInSubjet2 =
0267         getPFCandidates_exclJetConstituents<CandType>(*subjet2, pfCandidates, constitmap[2 * idx], false);
0268     if (verbosity_ >= 1) {
0269       std::cout << "#pfCandidatesNotInSubjet1 = " << pfCandidatesNotInSubjet1.size() << std::endl;
0270       std::cout << "#pfCandidatesNotInSubjet2 = " << pfCandidatesNotInSubjet2.size() << std::endl;
0271     }
0272
0273     // build JetToPFCandidateAssociation
0274     // (key = subjet, value = collection of PFCandidates that are not constituents of subjet)
0275     for (auto const& pfCandidate : pfCandidatesNotInSubjet1) {
0276       selectedSubjetPFCandidateAssociationForIsolation->insert(subjetRef1, pfCandidate);
0277     }
0278     for (auto const& pfCandidate : pfCandidatesNotInSubjet2) {
0279       selectedSubjetPFCandidateAssociationForIsolation->insert(subjetRef2, pfCandidate);
0280     }
0281   }
0282
0283   evt.put(std::move(selectedSubjets));
0284   evt.put(std::move(selectedSubjetPFCandidateAssociationForIsolation), "pfCandAssocMapForIsolation");
0285 }
0286
0287 typedef GenericBoostedTauSeedsProducer<reco::PFJet, reco::PFCandidate> BoostedTauSeedsProducer;
0288 typedef GenericBoostedTauSeedsProducer<pat::Jet, pat::PackedCandidate> PATBoostedTauSeedsProducer;
0289
0290 #include "FWCore/Framework/interface/MakerMacros.h"
0291 DEFINE_FWK_MODULE(BoostedTauSeedsProducer);
0292 DEFINE_FWK_MODULE(PATBoostedTauSeedsProducer);