NanoAOD/plugins/CandMCMatchTableProducer.cc

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, "  // continues below
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), "  // continues below
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       //std::cout << "cand #" << i << ": pT = " << cands->ptrAt(i)->pt() << ", eta = " << cands->ptrAt(i)->eta() << ", phi = " << cands->ptrAt(i)->phi() << std::endl;
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;  // go ahead with electrons, as those may be matched to a dressed lepton
0125
0126       switch (type_) {
0127         case MMuon:
0128           if (match->isPromptFinalState())
0129             flav[i] = 1;  // prompt
0130           else if (match->isDirectPromptTauDecayProductFinalState())
0131             flav[i] = 15;  // tau
0132           else
0133             flav[i] = getParentHadronFlag(match);  // 3 = light, 4 = charm, 5 = b
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);  // prompt electron or photon
0166           else if (match->isDirectPromptTauDecayProductFinalState())
0167             flav[i] = 15;  // tau
0168           else
0169             flav[i] = getParentHadronFlag(match);  // 3 = light, 4 = charm, 5 = b
0170           break;
0171         case MPhoton:
0172           if (match->isPromptFinalState() && match->pdgId() == 22)
0173             flav[i] = 1;  // prompt photon
0174           else if ((match->isPromptFinalState() || match->isDirectPromptTauDecayProductFinalState()) &&
0175                    abs(match->pdgId()) == 11)
0176             flav[i] = 11;  // prompt electron
0177           break;
0178         case MTau:
0179           // CV: assignment of status codes according to https://twiki.cern.ch/twiki/bin/viewauth/CMS/HiggsToTauTauWorking2016#MC_Matching
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());  // sanity
0209       if (mom.key() >= match.key())
0210         continue;  // prevent circular refs
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);