JetMCAlgos/plugins/PartonSelector.cc

0001 //
0002 // Select the partons status 2 and 3 for MC Jet Flavour
0003 // Author: Attilio
0004 // Date: 10.10.2007
0005 //
0006
0007 //=======================================================================
0008
0009 // user include files
0010 #include "FWCore/Framework/interface/global/EDProducer.h"
0011 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0012 #include "FWCore/Utilities/interface/InputTag.h"
0013
0014 #include "FWCore/Framework/interface/Event.h"
0015 #include "FWCore/Framework/interface/EventSetup.h"
0016 #include "FWCore/Framework/interface/MakerMacros.h"
0017 #include "FWCore/Framework/interface/ESHandle.h"
0018 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0019
0020 #include "DataFormats/Common/interface/Handle.h"
0021 #include "DataFormats/Common/interface/Ref.h"
0022 //#include "DataFormats/Candidate/interface/CandidateFwd.h"
0023 //#include "DataFormats/HepMCCandidate/interface/GenParticleCandidate.h"
0024 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0025 #include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
0026
0027 #include <memory>
0028 #include <string>
0029 #include <iostream>
0030 #include <vector>
0031
0032 using namespace std;
0033 using namespace reco;
0034 using namespace edm;
0035
0036 class PartonSelector : public edm::global::EDProducer<> {
0037 public:
0038   PartonSelector(const edm::ParameterSet&);
0039   ~PartonSelector() override;
0040
0041 private:
0042   void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0043   bool withLeptons;         // Optionally specify leptons
0044   bool withTop;             // Optionally include top quarks in the list
0045   bool acceptNoDaughters;   // Parton with zero daugthers are not considered by default, make it configurable
0046   unsigned int skipFirstN;  // Default skips first 6 particles, make it configurable
0047   edm::EDGetTokenT<reco::GenParticleCollection> tokenGenParticles_;  // input collection
0048 };
0049 //=========================================================================
0050
0051 PartonSelector::PartonSelector(const edm::ParameterSet& iConfig) {
0052   produces<reco::GenParticleRefVector>();
0053   withLeptons = iConfig.getParameter<bool>("withLeptons");
0054   tokenGenParticles_ = consumes<reco::GenParticleCollection>(iConfig.getParameter<edm::InputTag>("src"));
0055   if (iConfig.exists("acceptNoDaughters")) {
0056     acceptNoDaughters = iConfig.getParameter<bool>("acceptNoDaughters");
0057   } else {
0058     acceptNoDaughters = false;
0059   }
0060   if (iConfig.exists("skipFirstN")) {
0061     skipFirstN = iConfig.getParameter<unsigned int>("skipFirstN");
0062   } else {
0063     skipFirstN = 6;
0064   }
0065   if (iConfig.exists("withTop")) {
0066     withTop = iConfig.getParameter<bool>("withTop");
0067   } else {
0068     withTop = false;
0069   }
0070 }
0071
0072 //=========================================================================
0073
0074 PartonSelector::~PartonSelector() {}
0075
0076 // ------------ method called to produce the data  ------------
0077
0078 void PartonSelector::produce(StreamID, Event& iEvent, const EventSetup& iEs) const {
0079   //edm::Handle <reco::CandidateView> particles;
0080   edm::Handle<reco::GenParticleCollection> particles;
0081   iEvent.getByToken(tokenGenParticles_, particles);
0082   edm::LogVerbatim("PartonSelector") << "=== GenParticle size:" << particles->size();
0083   int nPart = 0;
0084
0085   auto thePartons = std::make_unique<GenParticleRefVector>();
0086
0087   for (size_t m = 0; m < particles->size(); m++) {
0088     // Don't take into account first 6 particles in generator list
0089     if (m < skipFirstN)
0090       continue;
0091
0092     const GenParticle& aParticle = (*particles)[m];
0093
0094     bool isAParton = false;
0095     bool isALepton = false;
0096     int flavour = abs(aParticle.pdgId());
0097     if (flavour == 1 || flavour == 2 || flavour == 3 || flavour == 4 || flavour == 5 || (flavour == 6 && withTop) ||
0098         flavour == 21)
0099       isAParton = true;
0100     if (flavour == 11 || flavour == 12 || flavour == 13 || flavour == 14 || flavour == 15 || flavour == 16)
0101       isALepton = true;
0102
0103     //Add Partons status 3
0104     if (aParticle.status() == 3 && isAParton) {
0105       thePartons->push_back(GenParticleRef(particles, m));
0106       nPart++;
0107     }
0108
0109     //Add Partons status 2
0110     int nparton_daughters = 0;
0111     if ((aParticle.numberOfDaughters() > 0 || acceptNoDaughters) && isAParton) {
0112       for (unsigned int i = 0; i < aParticle.numberOfDaughters(); i++) {
0113         int daughterFlavour = abs(aParticle.daughter(i)->pdgId());
0114         if ((daughterFlavour == 1 || daughterFlavour == 2 || daughterFlavour == 3 || daughterFlavour == 4 ||
0115              daughterFlavour == 5 || daughterFlavour == 6 || daughterFlavour == 21)) {
0116           nparton_daughters++;
0117         }
0118       }
0119       if (nparton_daughters == 0) {
0120         nPart++;
0121         thePartons->push_back(GenParticleRef(particles, m));
0122       }
0123     }
0124
0125     //Add Leptons
0126     // Here you have to decide what to do with taus....
0127     // Now all leptons, including e and mu from leptonic tau decays, are added
0128     if (withLeptons && aParticle.status() == 3 && isALepton) {
0129       thePartons->push_back(GenParticleRef(particles, m));
0130       nPart++;
0131     }
0132   }
0133
0134   edm::LogVerbatim("PartonSelector") << "=== GenParticle selected:" << nPart;
0135   iEvent.put(std::move(thePartons));
0136 }
0137
0138 DEFINE_FWK_MODULE(PartonSelector);