File indexing completed on 2025-05-08 02:18:43
0001 #include <DataFormats/Math/interface/deltaR.h>
0002
0003 #include <iostream>
0004 #include <vector>
0005
0006 #include "DataFormats/Candidate/interface/Candidate.h"
0007 #include "DataFormats/Common/interface/View.h"
0008 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0009 #include "DataFormats/JetReco/interface/GenJetCollection.h"
0010 #include "DataFormats/NanoAOD/interface/FlatTable.h"
0011 #include "FWCore/Framework/interface/Event.h"
0012 #include "FWCore/Framework/interface/global/EDProducer.h"
0013 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0014 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0015 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0016
0017 class CandMCMatchTableProducer : public edm::global::EDProducer<> {
0018 public:
0019 CandMCMatchTableProducer(edm::ParameterSet const& params)
0020 : objName_(params.getParameter<std::string>("objName")),
0021 branchName_(params.getParameter<std::string>("branchName")),
0022 doc_(params.getParameter<std::string>("docString")),
0023 src_(consumes<reco::CandidateView>(params.getParameter<edm::InputTag>("src"))),
0024 candMap_(consumes<edm::Association<reco::GenParticleCollection>>(params.getParameter<edm::InputTag>("mcMap"))) {
0025 produces<nanoaod::FlatTable>();
0026 const std::string& type = params.getParameter<std::string>("objType");
0027 if (type == "Muon")
0028 type_ = MMuon;
0029 else if (type == "Electron")
0030 type_ = MElectron;
0031 else if (type == "Tau")
0032 type_ = MTau;
0033 else if (type == "Photon")
0034 type_ = MPhoton;
0035 else if (type == "Track")
0036 type_ = MTrack;
0037 else if (type == "Other")
0038 type_ = MOther;
0039 else
0040 throw cms::Exception("Configuration", "Unsupported objType '" + type + "'\n");
0041
0042 switch (type_) {
0043 case MMuon:
0044 flavDoc_ =
0045 "1 = prompt muon (including gamma*->mu mu), 15 = muon from prompt "
0046 "tau, "
0047 "5 = muon from b, 4 = muon from c, 3 = muon from light or unknown, "
0048 "0 = unmatched";
0049 break;
0050 case MElectron:
0051 flavDoc_ =
0052 "1 = prompt electron (including gamma*->mu mu), 15 = electron from "
0053 "prompt tau, 22 = prompt photon (likely "
0054 "conversion), "
0055 "5 = electron from b, 4 = electron from c, 3 = electron from light "
0056 "or unknown, 0 = unmatched";
0057 break;
0058 case MPhoton:
0059 flavDoc_ = "1 = prompt photon, 11 = prompt electron, 0 = unknown or unmatched";
0060 break;
0061 case MTau:
0062 flavDoc_ =
0063 "1 = prompt electron, 2 = prompt muon, 3 = tau->e decay, 4 = "
0064 "tau->mu decay, 5 = hadronic tau decay, 0 = "
0065 "unknown or unmatched";
0066 break;
0067 case MTrack:
0068 flavDoc_ =
0069 "1 = prompt, 511 = from B0, 521 = from B+/-, 0 = unknown or "
0070 "unmatched";
0071 break;
0072 case MOther:
0073 flavDoc_ = "1 = from hard scatter, 0 = unknown or unmatched";
0074 break;
0075 }
0076
0077 if (type_ == MTau) {
0078 candMapVisTau_ =
0079 consumes<edm::Association<reco::GenParticleCollection>>(params.getParameter<edm::InputTag>("mcMapVisTau"));
0080 }
0081
0082 if (type_ == MElectron) {
0083 candMapDressedLep_ =
0084 consumes<edm::Association<reco::GenJetCollection>>(params.getParameter<edm::InputTag>("mcMapDressedLep"));
0085 mapTauAnc_ = consumes<edm::ValueMap<bool>>(params.getParameter<edm::InputTag>("mapTauAnc"));
0086 genPartsToken_ = consumes<reco::GenParticleCollection>(params.getParameter<edm::InputTag>("genparticles"));
0087 }
0088 }
0089
0090 ~CandMCMatchTableProducer() override {}
0091
0092 void produce(edm::StreamID id, edm::Event& iEvent, const edm::EventSetup& iSetup) const override {
0093 const auto& candProd = iEvent.get(src_);
0094 auto ncand = candProd.size();
0095
0096 auto tab = std::make_unique<nanoaod::FlatTable>(ncand, objName_, false, true);
0097
0098 const auto& mapProd = iEvent.get(candMap_);
0099
0100 edm::Handle<edm::Association<reco::GenParticleCollection>> mapVisTau;
0101 if (type_ == MTau) {
0102 iEvent.getByToken(candMapVisTau_, mapVisTau);
0103 }
0104
0105 edm::Handle<edm::Association<reco::GenJetCollection>> mapDressedLep;
0106 edm::Handle<edm::ValueMap<bool>> mapTauAnc;
0107 edm::Handle<reco::GenParticleCollection> genParts;
0108 if (type_ == MElectron) {
0109 iEvent.getByToken(candMapDressedLep_, mapDressedLep);
0110 iEvent.getByToken(mapTauAnc_, mapTauAnc);
0111 iEvent.getByToken(genPartsToken_, genParts);
0112 }
0113
0114 std::vector<int16_t> key(ncand, -1);
0115 std::vector<uint8_t> flav(ncand, 0);
0116 for (unsigned int i = 0; i < ncand; ++i) {
0117
0118
0119
0120 const auto& cand = candProd.ptrAt(i);
0121 reco::GenParticleRef match = mapProd[cand];
0122 reco::GenParticleRef matchVisTau;
0123 reco::GenJetRef matchDressedLep;
0124 bool hasTauAnc = false;
0125 if (type_ == MTau) {
0126 matchVisTau = (*mapVisTau)[cand];
0127 }
0128 if (type_ == MElectron) {
0129 matchDressedLep = (*mapDressedLep)[cand];
0130 if (matchDressedLep.isNonnull()) {
0131 hasTauAnc = (*mapTauAnc)[matchDressedLep];
0132 }
0133 }
0134 if (match.isNonnull())
0135 key[i] = match.key();
0136 else if (matchVisTau.isNonnull())
0137 key[i] = matchVisTau.key();
0138 else if (type_ != MElectron)
0139 continue;
0140
0141
0142 switch (type_) {
0143 case MMuon:
0144 if (match->isPromptFinalState())
0145 flav[i] = 1;
0146 else if (match->isDirectPromptTauDecayProductFinalState())
0147 flav[i] = 15;
0148 else
0149 flav[i] = getParentHadronFlag(match);
0150 break;
0151 case MElectron:
0152 if (matchDressedLep.isNonnull()) {
0153 if (matchDressedLep->pdgId() == 22)
0154 flav[i] = 22;
0155 else
0156 flav[i] = (hasTauAnc) ? 15 : 1;
0157
0158 float minpt = 0;
0159 const reco::GenParticle* highestPtConstituent = nullptr;
0160 for (auto& consti : matchDressedLep->getGenConstituents()) {
0161 if (abs(consti->pdgId()) != 11)
0162 continue;
0163 if (consti->pt() < minpt)
0164 continue;
0165 minpt = consti->pt();
0166 highestPtConstituent = consti;
0167 }
0168 if (highestPtConstituent) {
0169 auto iter =
0170 std::find_if(genParts->begin(), genParts->end(), [highestPtConstituent](reco::GenParticle genp) {
0171 return (abs(genp.pdgId()) == 11) && (deltaR(genp, *highestPtConstituent) < 0.01) &&
0172 (abs(genp.pt() - highestPtConstituent->pt()) / highestPtConstituent->pt() < 0.01);
0173 });
0174 if (iter != genParts->end()) {
0175 key[i] = iter - genParts->begin();
0176 }
0177 }
0178 } else if (!match.isNonnull())
0179 flav[i] = 0;
0180 else if (match->isPromptFinalState())
0181 flav[i] = (match->pdgId() == 22 ? 22 : 1);
0182 else if (match->isDirectPromptTauDecayProductFinalState())
0183 flav[i] = 15;
0184 else
0185 flav[i] = getParentHadronFlag(match);
0186 break;
0187 case MPhoton:
0188 if (match->isPromptFinalState() && match->pdgId() == 22)
0189 flav[i] = 1;
0190 else if ((match->isPromptFinalState() || match->isDirectPromptTauDecayProductFinalState()) &&
0191 abs(match->pdgId()) == 11)
0192 flav[i] = 11;
0193 break;
0194 case MTau:
0195
0196
0197 if (match.isNonnull() && match->statusFlags().isPrompt() && abs(match->pdgId()) == 11)
0198 flav[i] = 1;
0199 else if (match.isNonnull() && match->statusFlags().isPrompt() && abs(match->pdgId()) == 13)
0200 flav[i] = 2;
0201 else if (match.isNonnull() && match->isDirectPromptTauDecayProductFinalState() && abs(match->pdgId()) == 11)
0202 flav[i] = 3;
0203 else if (match.isNonnull() && match->isDirectPromptTauDecayProductFinalState() && abs(match->pdgId()) == 13)
0204 flav[i] = 4;
0205 else if (matchVisTau.isNonnull())
0206 flav[i] = 5;
0207 break;
0208 case MTrack:
0209 if (match->isPromptFinalState())
0210 flav[i] = 1;
0211 else
0212 flav[i] = getParentHadronFlag(match);
0213 break;
0214 default:
0215 flav[i] = match->statusFlags().fromHardProcess();
0216 };
0217 }
0218
0219 tab->addColumn<int16_t>(branchName_ + "Idx", key, "Index into genParticle list for " + doc_);
0220 tab->addColumn<uint8_t>(branchName_ + "Flav",
0221 flav,
0222 "Flavour of genParticle (DressedLeptons for electrons) for " + doc_ + ": " + flavDoc_);
0223
0224 iEvent.put(std::move(tab));
0225 }
0226
0227 static int getParentHadronFlag(const reco::GenParticleRef match) {
0228 bool has4 = false;
0229 for (unsigned int im = 0, nm = match->numberOfMothers(); im < nm; ++im) {
0230 reco::GenParticleRef mom = match->motherRef(im);
0231 assert(mom.isNonnull() && mom.isAvailable());
0232 if (mom.key() >= match.key())
0233 continue;
0234 int id = std::abs(mom->pdgId());
0235 if (id / 1000 == 5 || id / 100 == 5 || id == 5)
0236 return 5;
0237 if (id / 1000 == 4 || id / 100 == 4 || id == 4)
0238 has4 = true;
0239 if (mom->status() == 2) {
0240 id = getParentHadronFlag(mom);
0241 if (id == 5)
0242 return 5;
0243 else if (id == 4)
0244 has4 = true;
0245 }
0246 }
0247 return has4 ? 4 : 3;
0248 }
0249
0250 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0251 edm::ParameterSetDescription desc;
0252 desc.add<std::string>("objName")->setComment("name of the nanoaod::FlatTable to extend with this table");
0253 desc.add<std::string>("branchName")
0254 ->setComment(
0255 "name of the column to write (the final branch in the nanoaod will "
0256 "be <objName>_<branchName>Idx and "
0257 "<objName>_<branchName>Flav");
0258 desc.add<std::string>("docString")->setComment("documentation to forward to the output");
0259 desc.add<edm::InputTag>("src")->setComment(
0260 "physics object collection for the reconstructed objects (e.g. "
0261 "leptons)");
0262 desc.add<edm::InputTag>("mcMap")->setComment(
0263 "tag to an edm::Association<GenParticleCollection> mapping src to gen, "
0264 "such as the one produced by MCMatcher");
0265 desc.add<std::string>("objType")->setComment(
0266 "type of object to match (Muon, Electron, Tau, Photon, Other), taylors "
0267 "what's in t Flav branch");
0268 desc.addOptional<edm::InputTag>("mcMapVisTau")
0269 ->setComment(
0270 "as mcMap, but pointing to the visible gen taus (only if objType "
0271 "== Tau)");
0272 desc.addOptional<edm::InputTag>("mcMapDressedLep")
0273 ->setComment(
0274 "as mcMap, but pointing to gen dressed leptons (only if objType == "
0275 "Electrons)");
0276 desc.addOptional<edm::InputTag>("mapTauAnc")
0277 ->setComment(
0278 "Value map of matched gen electrons containing info on the tau "
0279 "ancestry");
0280 desc.addOptional<edm::InputTag>("genparticles")->setComment("Collection of genParticles to be stored.");
0281 descriptions.add("candMcMatchTable", desc);
0282 }
0283
0284 protected:
0285 const std::string objName_, branchName_, doc_;
0286 const edm::EDGetTokenT<reco::CandidateView> src_;
0287 const edm::EDGetTokenT<edm::Association<reco::GenParticleCollection>> candMap_;
0288 edm::EDGetTokenT<edm::Association<reco::GenParticleCollection>> candMapVisTau_;
0289 edm::EDGetTokenT<edm::Association<reco::GenJetCollection>> candMapDressedLep_;
0290 edm::EDGetTokenT<edm::ValueMap<bool>> mapTauAnc_;
0291 edm::EDGetTokenT<reco::GenParticleCollection> genPartsToken_;
0292 enum MatchType { MMuon, MElectron, MTau, MPhoton, MTrack, MOther } type_;
0293 std::string flavDoc_;
0294 };
0295
0296 #include "FWCore/Framework/interface/MakerMacros.h"
0297 DEFINE_FWK_MODULE(CandMCMatchTableProducer);