Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-05-08 02:18:42 

0001 // original author: RK18 team
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/NanoAOD/interface/FlatTable.h"
0010 #include "FWCore/Framework/interface/Event.h"
0011 #include "FWCore/Framework/interface/global/EDProducer.h"
0012 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0013 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0014 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0015 
0016 class MCFinalStateSelector : public edm::global::EDProducer<> {
0017 public:
0018   MCFinalStateSelector(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 == "ProbeTracks")
0035       type_ = MTrack;
0036     else if (type == "Kshort")
0037       type_ = MKshort;
0038     else if (type == "Other")
0039       type_ = MOther;
0040     else
0041       throw cms::Exception("Configuration", "Unsupported objType '" + type + "'\n");
0042 
0043     switch (type_) {
0044       case MMuon:
0045         flavDoc_ =
0046             "1 = prompt muon (including gamma*->mu mu), 15 = muon from prompt "
0047             "tau, "  // continues below
0048             "511 = from B0, 521 = from B+/-, 0 = unknown or unmatched";
0049         break;
0050 
0051       case MElectron:
0052         flavDoc_ =
0053             "1 = prompt electron (including gamma*->mu mu), 15 = electron from "
0054             "prompt tau, "                              // continues below
0055             "22 = prompt photon (likely conversion), "  // continues below
0056             "511 = from B0, 521 = from B+/-, 0 = unknown or unmatched";
0057         break;
0058       case MPhoton:
0059         flavDoc_ = "1 = prompt photon, 13 = 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 = unknown or unmatched";
0065         break;
0066       case MTrack:
0067         flavDoc_ =
0068             "1 = prompt, 511 = from B0, 521 = from B+/-, 0 = unknown or "
0069             "unmatched";
0070         break;
0071 
0072       case MOther:
0073         flavDoc_ = "1 = from hard scatter, 0 = unknown or unmatched";
0074         break;
0075       case MKshort:
0076         flavDoc_ = "1 = from hard scatter, 0 = unknown or unmatched";
0077         break;
0078     }
0079 
0080     if (type_ == MTau) {
0081       candMapVisTau_ =
0082           consumes<edm::Association<reco::GenParticleCollection>>(params.getParameter<edm::InputTag>("mcMapVisTau"));
0083     }
0084   }
0085 
0086   ~MCFinalStateSelector() override {}
0087 
0088   void produce(edm::StreamID id, edm::Event& iEvent, const edm::EventSetup& iSetup) const override {
0089     edm::Handle<reco::CandidateView> cands;
0090     iEvent.getByToken(src_, cands);
0091     unsigned int ncand = cands->size();
0092 
0093     auto tab = std::make_unique<nanoaod::FlatTable>(ncand, objName_, false, true);
0094 
0095     edm::Handle<edm::Association<reco::GenParticleCollection>> map;
0096     iEvent.getByToken(candMap_, map);
0097 
0098     edm::Handle<edm::Association<reco::GenParticleCollection>> mapVisTau;
0099     if (type_ == MTau) {
0100       iEvent.getByToken(candMapVisTau_, mapVisTau);
0101     }
0102 
0103     std::vector<int> key(ncand, -1), flav(ncand, 0);
0104     for (unsigned int i = 0; i < ncand; ++i) {
0105       // std::cout << "cand #" << i << ": pT = " << cands->ptrAt(i)->pt() << ",
0106       // eta = " << cands->ptrAt(i)->eta() << ", phi = " <<
0107       // cands->ptrAt(i)->phi() << std::endl;
0108       reco::GenParticleRef match = (*map)[cands->ptrAt(i)];
0109       reco::GenParticleRef matchVisTau;
0110       if (type_ == MTau) {
0111         matchVisTau = (*mapVisTau)[cands->ptrAt(i)];
0112       }
0113       if (match.isNonnull())
0114         key[i] = match.key();
0115       else if (matchVisTau.isNonnull())
0116         key[i] = matchVisTau.key();
0117       else
0118         continue;
0119       switch (type_) {
0120         case MMuon:
0121           if (match->isPromptFinalState())
0122             flav[i] = 1;  // prompt
0123           else
0124             flav[i] = getParentHadronFlag(match);  // pdgId of mother
0125           break;
0126         case MElectron:
0127           if (match->isPromptFinalState())
0128             flav[i] = (match->pdgId() == 22 ? 22 : 1);  // prompt electron or photon
0129           else
0130             flav[i] = getParentHadronFlag(match);  // pdgId of mother
0131           break;
0132         case MPhoton:
0133           if (match->isPromptFinalState())
0134             flav[i] = (match->pdgId() == 22 ? 1 : 13);  // prompt electron or photon
0135           break;
0136         case MTau:
0137           // CV: assignment of status codes according to
0138           // https://twiki.cern.ch/twiki/bin/viewauth/CMS/HiggsToTauTauWorking2016#MC_Matching
0139           if (match.isNonnull() && match->statusFlags().isPrompt() && abs(match->pdgId()) == 11)
0140             flav[i] = 1;
0141           else if (match.isNonnull() && match->statusFlags().isPrompt() && abs(match->pdgId()) == 13)
0142             flav[i] = 2;
0143           else if (match.isNonnull() && match->isDirectPromptTauDecayProductFinalState() && abs(match->pdgId()) == 11)
0144             flav[i] = 3;
0145           else if (match.isNonnull() && match->isDirectPromptTauDecayProductFinalState() && abs(match->pdgId()) == 13)
0146             flav[i] = 4;
0147           else if (matchVisTau.isNonnull())
0148             flav[i] = 5;
0149           break;
0150         case MTrack:
0151           if (match->isPromptFinalState())
0152             flav[i] = 1;  // prompt
0153           else
0154             flav[i] = getParentHadronFlag(match);  // pdgId of mother
0155           break;
0156         case MKshort:
0157           if (match->isPromptFinalState())
0158             flav[i] = 1;  // prompt
0159           else
0160             flav[i] = getParentHadronFlag(match);  // pdgId of mother
0161           break;
0162         default:
0163           flav[i] = match->statusFlags().fromHardProcess();
0164       };
0165     }
0166 
0167     tab->addColumn<int>(branchName_ + "Idx", key, "Index into genParticle list for " + doc_);
0168     // for(auto ij : flav) std::cout << " flav = " << ij << " " << (uint8_t)ij
0169     // << std::endl; tab->addColumn<uint8_t>(branchName_+"Flav", flav, "Flavour
0170     // of genParticle for "+doc_+": "+flavDoc_,
0171     // nanoaod::FlatTable::UInt8Column);
0172     tab->addColumn<int>(branchName_ + "Flav", flav, "Flavour of genParticle for " + doc_ + ": " + flavDoc_);
0173 
0174     iEvent.put(std::move(tab));
0175   }
0176 
0177   static int getParentHadronFlag(const reco::GenParticleRef match) {
0178     for (unsigned int im = 0, nm = match->numberOfMothers(); im < nm; ++im) {
0179       reco::GenParticleRef mom = match->motherRef(im);
0180       assert(mom.isNonnull() && mom.isAvailable());  // sanity
0181       if (mom.key() >= match.key())
0182         continue;  // prevent circular refs
0183       int id = std::abs(mom->pdgId());
0184       return id;
0185     }
0186     return 0;
0187   }
0188 
0189 protected:
0190   const std::string objName_, branchName_, doc_;
0191   const edm::EDGetTokenT<reco::CandidateView> src_;
0192   const edm::EDGetTokenT<edm::Association<reco::GenParticleCollection>> candMap_;
0193   edm::EDGetTokenT<edm::Association<reco::GenParticleCollection>> candMapVisTau_;
0194   enum MatchType { MMuon, MElectron, MTau, MPhoton, MTrack, MOther, MKshort } type_;
0195   std::string flavDoc_;
0196 };
0197 
0198 #include "FWCore/Framework/interface/MakerMacros.h"
0199 DEFINE_FWK_MODULE(MCFinalStateSelector);
[ 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!