Back to home page

Project CMSSW displayed by LXR
 
 

    

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, "  // continues below
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), "  // continues below
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       // std::cout << "cand #" << i << ": pT = " << cands->ptrAt(i)->pt() << ",
0118       // eta = " << cands->ptrAt(i)->eta() << ", phi = " <<
0119       // cands->ptrAt(i)->phi() << std::endl;
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;  // go ahead with electrons, as those may be matched to a
0140                    // dressed lepton
0141 
0142       switch (type_) {
0143         case MMuon:
0144           if (match->isPromptFinalState())
0145             flav[i] = 1;  // prompt
0146           else if (match->isDirectPromptTauDecayProductFinalState())
0147             flav[i] = 15;  // tau
0148           else
0149             flav[i] = getParentHadronFlag(match);  // 3 = light, 4 = charm, 5 = b
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);  // prompt electron or photon
0182           else if (match->isDirectPromptTauDecayProductFinalState())
0183             flav[i] = 15;  // tau
0184           else
0185             flav[i] = getParentHadronFlag(match);  // 3 = light, 4 = charm, 5 = b
0186           break;
0187         case MPhoton:
0188           if (match->isPromptFinalState() && match->pdgId() == 22)
0189             flav[i] = 1;  // prompt photon
0190           else if ((match->isPromptFinalState() || match->isDirectPromptTauDecayProductFinalState()) &&
0191                    abs(match->pdgId()) == 11)
0192             flav[i] = 11;  // prompt electron
0193           break;
0194         case MTau:
0195           // CV: assignment of status codes according to
0196           // https://twiki.cern.ch/twiki/bin/viewauth/CMS/HiggsToTauTauWorking2016#MC_Matching
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;  // prompt
0211           else
0212             flav[i] = getParentHadronFlag(match);  // pdgId of mother
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());  // sanity
0232       if (mom.key() >= match.key())
0233         continue;  // prevent circular refs
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);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!