JetMCAlgos/plugins/HadronAndPartonSelector.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    HadronAndPartonSelector
0004 // Class:      HadronAndPartonSelector
0005 //
0006 /**\class HadronAndPartonSelector HadronAndPartonSelector.cc PhysicsTools/JetMCAlgos/plugins/HadronAndPartonSelector.cc
0007  * \brief Selects hadrons and partons from a collection of GenParticles
0008  *
0009  * This producer selects hadrons, partons, and leptons from a collection of GenParticles and stores vectors of EDM references
0010  * to these particles in the event. The following hadrons are selected:
0011  *
0012  * - b hadrons that do not have other b hadrons as daughters
0013  *
0014  * - c hadrons that do not have other c hadrons as daughters
0015  *
0016  * Older Fortran Monte Carlo generators (Pythia6 and Herwig6) follow the HEPEVT [1] particle status code convention while
0017  * newer C++ Monte Carlo generators (Pythia8, Herwig++, and Sherpa) follow the HepMC [2] particle status code convention.
0018  * However, both conventions give considerable freedom in defining the status codes of intermediate particle states. Hence,
0019  * the parton selection is generator-dependent and is described in each of the parton selectors separately.
0020  *
0021  * Using the provenance information of the GenEventInfoProduct, the producer attempts to automatically determine what generator
0022  * was used to hadronize events and based on that information decides what parton selection mode to use. It is also possible
0023  * to enforce any of the supported parton selection modes.
0024  *
0025  * The selected hadrons and partons are finally used by the JetFlavourClustering producer to determine the jet flavour.
0026  *
0027  * The following leptons are selected:
0028  *
0029  * - status==1 electrons and muons
0030  *
0031  * - status==2 taus
0032  *
0033  *
0034  * [1] http://cepa.fnal.gov/psm/stdhep/
0035  * [2] http://lcgapp.cern.ch/project/simu/HepMC/
0036  */
0037 //
0038 // Original Author:  Dinko Ferencek
0039 //         Created:  Tue Nov  5 22:43:43 CET 2013
0040 //
0041 //
0042
0043 // system include files
0044 #include <memory>
0045
0046 // user include files
0047 #include "FWCore/Framework/interface/Frameworkfwd.h"
0048 #include "FWCore/Framework/interface/stream/EDProducer.h"
0049
0050 #include "FWCore/Framework/interface/Event.h"
0051 #include "FWCore/Framework/interface/MakerMacros.h"
0052
0053 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0054
0055 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
0056 #include "FWCore/Common/interface/Provenance.h"
0057 #include "DataFormats/Candidate/interface/Candidate.h"
0058 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0059 #include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
0060 #include "PhysicsTools/JetMCUtils/interface/CandMCTag.h"
0061 #include "CommonTools/CandUtils/interface/pdgIdUtils.h"
0062 #include "PhysicsTools/JetMCAlgos/interface/BasePartonSelector.h"
0063 #include "PhysicsTools/JetMCAlgos/interface/Pythia6PartonSelector.h"
0064 #include "PhysicsTools/JetMCAlgos/interface/Pythia8PartonSelector.h"
0065 #include "PhysicsTools/JetMCAlgos/interface/HerwigppPartonSelector.h"
0066 #include "PhysicsTools/JetMCAlgos/interface/SherpaPartonSelector.h"
0067
0068 //
0069 // constants, enums and typedefs
0070 //
0071 typedef std::shared_ptr<BasePartonSelector> PartonSelectorPtr;
0072
0073 //
0074 // class declaration
0075 //
0076
0077 class HadronAndPartonSelector : public edm::stream::EDProducer<> {
0078 public:
0079   explicit HadronAndPartonSelector(const edm::ParameterSet&);
0080   ~HadronAndPartonSelector() override;
0081
0082   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0083
0084 private:
0085   void produce(edm::Event&, const edm::EventSetup&) override;
0086
0087   // ----------member data ---------------------------
0088   const edm::EDGetTokenT<GenEventInfoProduct> srcToken_;                // To get handronizer module type
0089   const edm::EDGetTokenT<reco::GenParticleCollection> particlesToken_;  // Input GenParticle collection
0090
0091   std::string partonMode_;  // Parton selection mode
0092   bool fullChainPhysPartons_;
0093   bool partonSelectorSet_;
0094   PartonSelectorPtr partonSelector_;
0095 };
0096
0097 //
0098 // static data member definitions
0099 //
0100
0101 //
0102 // constructors and destructor
0103 //
0104 HadronAndPartonSelector::HadronAndPartonSelector(const edm::ParameterSet& iConfig)
0105     :
0106
0107       srcToken_(mayConsume<GenEventInfoProduct>(iConfig.getParameter<edm::InputTag>("src"))),
0108       particlesToken_(consumes<reco::GenParticleCollection>(iConfig.getParameter<edm::InputTag>("particles"))),
0109       partonMode_(iConfig.getParameter<std::string>("partonMode")),
0110       fullChainPhysPartons_(iConfig.getParameter<bool>("fullChainPhysPartons"))
0111
0112 {
0113   //register your products
0114   produces<reco::GenParticleRefVector>("bHadrons");
0115   produces<reco::GenParticleRefVector>("cHadrons");
0116   produces<reco::GenParticleRefVector>("algorithmicPartons");
0117   produces<reco::GenParticleRefVector>("physicsPartons");
0118   produces<reco::GenParticleRefVector>("leptons");
0119
0120   partonSelectorSet_ = false;
0121   partonSelector_ = nullptr;
0122 }
0123
0124 HadronAndPartonSelector::~HadronAndPartonSelector() {
0125   // do anything here that needs to be done at desctruction time
0126   // (e.g. close files, deallocate resources etc.)
0127 }
0128
0129 //
0130 // member functions
0131 //
0132
0133 // ------------ method called to produce the data  ------------
0134 void HadronAndPartonSelector::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0135   // determine hadronizer type (done only once per job)
0136   if (partonMode_ == "Auto") {
0137     edm::Handle<GenEventInfoProduct> genEvtInfoProduct;
0138     iEvent.getByToken(srcToken_, genEvtInfoProduct);
0139
0140     std::string moduleName = "";
0141     if (genEvtInfoProduct.isValid()) {
0142       const edm::StableProvenance& prov = iEvent.getStableProvenance(genEvtInfoProduct.id());
0143       moduleName = edm::moduleName(prov, iEvent.processHistory());
0144       if (moduleName == "ExternalGeneratorFilter") {
0145         moduleName = edm::parameterSet(prov, iEvent.processHistory()).getParameter<std::string>("@external_type");
0146         edm::LogInfo("SpecialModule") << "GEN events are produced by ExternalGeneratorFilter, "
0147                                       << "which is a wrapper of the original module: " << moduleName;
0148       }
0149     }
0150
0151     if (moduleName.find("Pythia6") != std::string::npos)
0152       partonMode_ = "Pythia6";
0153     else if (moduleName.find("Pythia8") != std::string::npos)
0154       partonMode_ = "Pythia8";
0155     else if (moduleName.find("ThePEG") != std::string::npos)
0156       partonMode_ = "Herwig++";
0157     else if (moduleName.find("Herwig7") != std::string::npos)
0158       partonMode_ = "Herwig++";
0159     else if (moduleName.find("Sherpa") != std::string::npos)
0160       partonMode_ = "Sherpa";
0161     else
0162       partonMode_ = "Undefined";
0163   }
0164
0165   // set the parton selection mode (done only once per job)
0166   if (!partonSelectorSet_) {
0167     if (partonMode_ == "Undefined")
0168       edm::LogWarning("UndefinedPartonMode")
0169           << "Could not automatically determine the hadronizer type and set the correct parton selection mode. "
0170              "Parton-based jet flavour will not be defined.";
0171     else if (partonMode_ == "Pythia6") {
0172       partonSelector_ = PartonSelectorPtr(new Pythia6PartonSelector());
0173       edm::LogInfo("PartonModeDefined") << "Using Pythia6 parton selection mode.";
0174     } else if (partonMode_ == "Pythia8") {
0175       partonSelector_ = PartonSelectorPtr(new Pythia8PartonSelector());
0176       edm::LogInfo("PartonModeDefined") << "Using Pythia8 parton selection mode.";
0177     } else if (partonMode_ == "Herwig++") {
0178       partonSelector_ = PartonSelectorPtr(new HerwigppPartonSelector());
0179       edm::LogInfo("PartonModeDefined") << "Using Herwig++ parton selection mode.";
0180     } else if (partonMode_ == "Sherpa") {
0181       partonSelector_ = PartonSelectorPtr(new SherpaPartonSelector());
0182       edm::LogInfo("PartonModeDefined") << "Using Sherpa parton selection mode.";
0183     } else
0184       throw cms::Exception("InvalidPartonMode") << "Parton selection mode is invalid: " << partonMode_
0185                                                 << ", use Auto | Pythia6 | Pythia8 | Herwig++ | Sherpa" << std::endl;
0186
0187     partonSelectorSet_ = true;
0188   }
0189
0190   edm::Handle<reco::GenParticleCollection> particles;
0191   iEvent.getByToken(particlesToken_, particles);
0192
0193   auto bHadrons = std::make_unique<reco::GenParticleRefVector>();
0194   auto cHadrons = std::make_unique<reco::GenParticleRefVector>();
0195   auto partons = std::make_unique<reco::GenParticleRefVector>();
0196   auto physicsPartons = std::make_unique<reco::GenParticleRefVector>();
0197   auto leptons = std::make_unique<reco::GenParticleRefVector>();
0198
0199   // loop over particles and select b and c hadrons and leptons
0200   for (reco::GenParticleCollection::const_iterator it = particles->begin(); it != particles->end(); ++it) {
0201     // if b hadron
0202     if (CandMCTagUtils::hasBottom(*it)) {
0203       // check if any of the daughters is also a b hadron
0204       bool hasbHadronDaughter = false;
0205       for (size_t i = 0; i < it->numberOfDaughters(); ++i) {
0206         if (CandMCTagUtils::hasBottom(*(it->daughter(i)))) {
0207           hasbHadronDaughter = true;
0208           break;
0209         }
0210       }
0211       if (hasbHadronDaughter)
0212         continue;  // skip excited b hadrons that have other b hadrons as daughters
0213
0214       bHadrons->push_back(reco::GenParticleRef(particles, it - particles->begin()));
0215     }
0216
0217     // if c hadron
0218     if (CandMCTagUtils::hasCharm(*it)) {
0219       // check if any of the daughters is also a c hadron
0220       bool hascHadronDaughter = false;
0221       for (size_t i = 0; i < it->numberOfDaughters(); ++i) {
0222         if (CandMCTagUtils::hasCharm(*(it->daughter(i)))) {
0223           hascHadronDaughter = true;
0224           break;
0225         }
0226       }
0227       if (hascHadronDaughter)
0228         continue;  // skip excited c hadrons that have other c hadrons as daughters
0229
0230       cHadrons->push_back(reco::GenParticleRef(particles, it - particles->begin()));
0231     }
0232
0233     // status==1 electrons and muons
0234     if ((reco::isElectron(*it) || reco::isMuon(*it)) && it->status() == 1)
0235       leptons->push_back(reco::GenParticleRef(particles, it - particles->begin()));
0236
0237     // status==2 taus
0238     if (reco::isTau(*it) && it->status() == 2)
0239       leptons->push_back(reco::GenParticleRef(particles, it - particles->begin()));
0240   }
0241
0242   // select algorithmic partons
0243   if (partonMode_ != "Undefined") {
0244     partonSelector_->run(particles, partons);
0245   }
0246
0247   // select physics partons
0248   for (reco::GenParticleCollection::const_iterator it = particles->begin(); it != particles->end(); ++it) {
0249     if (!fullChainPhysPartons_) {
0250       if (!(it->status() == 3 || ((partonMode_ == "Pythia8") && (it->status() == 23))))
0251         continue;
0252     }
0253     if (!CandMCTagUtils::isParton(*it))
0254       continue;  // skip particle if not a parton
0255     physicsPartons->push_back(reco::GenParticleRef(particles, it - particles->begin()));
0256   }
0257
0258   iEvent.put(std::move(bHadrons), "bHadrons");
0259   iEvent.put(std::move(cHadrons), "cHadrons");
0260   iEvent.put(std::move(partons), "algorithmicPartons");
0261   iEvent.put(std::move(physicsPartons), "physicsPartons");
0262   iEvent.put(std::move(leptons), "leptons");
0263 }
0264
0265 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0266 void HadronAndPartonSelector::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0267   //The following says we do not know what parameters are allowed so do no validation
0268   // Please change this to state exactly what you do use, even if it is no parameters
0269   edm::ParameterSetDescription desc;
0270   desc.setUnknown();
0271   descriptions.addDefault(desc);
0272 }
0273
0274 //define this as a plug-in
0275 DEFINE_FWK_MODULE(HadronAndPartonSelector);