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