JetProducers/plugins/InputGenJetsParticleSelector.cc

0001 /* \class GenJetInputParticleSelector
0002 *
0003 *  Selects particles that are used as input for the GenJet collection.
0004 *  Logic: select all stable particles, except for particles specified in
0005 *  the config file that come from
0006 *  W,Z and H decays, and except for a special list, which can be used for
0007 *  unvisible BSM-particles.
0008 *  It is also possible to only selected the partonic final state,
0009 *  which means all particles before the hadronization step.
0010 *
0011 *  The algorithm is based on code of Christophe Saout.
0012 *
0013 *  Usage: [example for no resonance from nu an mu, and deselect invisible BSM
0014 *         particles ]
0015 *
0016 *  module genJetParticles = InputGenJetsParticleSelector {
0017 *                InputTag src = "genParticles"
0018 *                bool partonicFinalState = false
0019 *                bool excludeResonances = true
0020 *                vuint32 excludeFromResonancePids = {13,12,14,16}
0021 *                bool tausAsJets = false
0022 *                vuint32 ignoreParticleIDs = {   1000022, 2000012, 2000014,
0023 *                                                2000016, 1000039, 5000039,
0024 *                                                4000012, 9900012, 9900014,
0025 *                                                9900016, 39}
0026 *        }
0027 *
0028 *
0029 * \author: Christophe Saout, Andreas Oehler
0030 *
0031 * Modifications:
0032 *
0033 *    04.08.2014: Dinko Ferencek
0034 *                Added support for Pythia8 (status=22 for intermediate resonances)
0035 *    23.09.2014: Dinko Ferencek
0036 *                Generalized code to work with miniAOD (except for the partonicFinalState which requires AOD)
0037 *
0038 */
0039
0040 #include "InputGenJetsParticleSelector.h"
0041 #include "FWCore/Framework/interface/MakerMacros.h"
0042 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0043 //#include <iostream>
0044 #include <memory>
0045 #include "CommonTools/CandUtils/interface/pdgIdUtils.h"
0046
0047 using namespace std;
0048
0049 InputGenJetsParticleSelector::InputGenJetsParticleSelector(const edm::ParameterSet &params)
0050     : inTag(params.getParameter<edm::InputTag>("src")),
0051       prunedInTag(params.exists("prunedGenParticles") ? params.getParameter<edm::InputTag>("prunedGenParticles")
0052                                                       : edm::InputTag("prunedGenParticles")),
0053       partonicFinalState(params.getParameter<bool>("partonicFinalState")),
0054       excludeResonances(params.getParameter<bool>("excludeResonances")),
0055       tausAsJets(params.getParameter<bool>("tausAsJets")),
0056       ptMin(0.0) {
0057   if (params.exists("ignoreParticleIDs"))
0058     setIgnoredParticles(params.getParameter<std::vector<unsigned int> >("ignoreParticleIDs"));
0059   setExcludeFromResonancePids(params.getParameter<std::vector<unsigned int> >("excludeFromResonancePids"));
0060   isMiniAOD =
0061       (params.exists("isMiniAOD") ? params.getParameter<bool>("isMiniAOD") : (inTag.label() == "packedGenParticles"));
0062
0063   if (isMiniAOD && partonicFinalState) {
0064     edm::LogError("PartonicFinalStateFromMiniAOD")
0065         << "Partonic final state not supported for MiniAOD. Falling back to the stable particle selection.";
0066     partonicFinalState = false;
0067   }
0068
0069   produces<reco::CandidatePtrVector>();
0070
0071   input_genpartcoll_token_ = consumes<reco::CandidateView>(inTag);
0072   if (isMiniAOD)
0073     input_prunedgenpartcoll_token_ = consumes<reco::CandidateView>(prunedInTag);
0074 }
0075
0076 InputGenJetsParticleSelector::~InputGenJetsParticleSelector() {}
0077
0078 void InputGenJetsParticleSelector::setIgnoredParticles(const std::vector<unsigned int> &particleIDs) {
0079   ignoreParticleIDs = particleIDs;
0080   std::sort(ignoreParticleIDs.begin(), ignoreParticleIDs.end());
0081 }
0082
0083 void InputGenJetsParticleSelector::setExcludeFromResonancePids(const std::vector<unsigned int> &particleIDs) {
0084   excludeFromResonancePids = particleIDs;
0085   std::sort(excludeFromResonancePids.begin(), excludeFromResonancePids.end());
0086 }
0087
0088 bool InputGenJetsParticleSelector::isParton(int pdgId) const {
0089   pdgId = (pdgId > 0 ? pdgId : -pdgId) % 10000;
0090   return (pdgId > 0 && pdgId < 6) || pdgId == 9 || (tausAsJets && pdgId == 15) || pdgId == 21;
0091   // tops are not considered "regular" partons
0092   // but taus eventually are (since they may hadronize later)
0093 }
0094
0095 bool InputGenJetsParticleSelector::isHadron(int pdgId) {
0096   pdgId = (pdgId > 0 ? pdgId : -pdgId) % 10000;
0097   return (pdgId > 100 && pdgId < 900) || (pdgId > 1000 && pdgId < 9000);
0098 }
0099
0100 bool InputGenJetsParticleSelector::isResonance(int pdgId) {
0101   // gauge bosons and tops
0102   pdgId = (pdgId > 0 ? pdgId : -pdgId) % 10000;
0103   return (pdgId > 21 && pdgId <= 42) || pdgId == 6 || pdgId == 7 || pdgId == 8;  //BUG! was 21. 22=gamma..
0104 }
0105
0106 bool InputGenJetsParticleSelector::isIgnored(int pdgId) const {
0107   pdgId = pdgId > 0 ? pdgId : -pdgId;
0108   std::vector<unsigned int>::const_iterator pos =
0109       std::lower_bound(ignoreParticleIDs.begin(), ignoreParticleIDs.end(), (unsigned int)pdgId);
0110   return pos != ignoreParticleIDs.end() && *pos == (unsigned int)pdgId;
0111 }
0112
0113 bool InputGenJetsParticleSelector::isExcludedFromResonance(int pdgId) const {
0114   pdgId = pdgId > 0 ? pdgId : -pdgId;
0115   std::vector<unsigned int>::const_iterator pos =
0116       std::lower_bound(excludeFromResonancePids.begin(), excludeFromResonancePids.end(), (unsigned int)pdgId);
0117   return pos != excludeFromResonancePids.end() && *pos == (unsigned int)pdgId;
0118 }
0119
0120 static unsigned int partIdx(const InputGenJetsParticleSelector::ParticleVector &p, const reco::Candidate *particle) {
0121   InputGenJetsParticleSelector::ParticleVector::const_iterator pos = std::lower_bound(p.begin(), p.end(), particle);
0122   if (pos == p.end() || *pos != particle)
0123     throw cms::Exception("CorruptedData") << "reco::GenEvent corrupted: Unlisted particles"
0124                                              " in decay tree."
0125                                           << std::endl;
0126
0127   return pos - p.begin();
0128 }
0129
0130 static void invalidateTree(InputGenJetsParticleSelector::ParticleBitmap &invalid,
0131                            const InputGenJetsParticleSelector::ParticleVector &p,
0132                            const reco::Candidate *particle) {
0133   unsigned int npart = particle->numberOfDaughters();
0134   if (!npart)
0135     return;
0136
0137   for (unsigned int i = 0; i < npart; ++i) {
0138     unsigned int idx = partIdx(p, particle->daughter(i));
0139     if (invalid[idx])
0140       continue;
0141     invalid[idx] = true;
0142     //cout<<"Invalidated: ["<<setw(4)<<idx<<"] With pt:"<<particle->daughter(i)->pt()<<endl;
0143     invalidateTree(invalid, p, particle->daughter(i));
0144   }
0145 }
0146
0147 int InputGenJetsParticleSelector::testPartonChildren(InputGenJetsParticleSelector::ParticleBitmap &invalid,
0148                                                      const InputGenJetsParticleSelector::ParticleVector &p,
0149                                                      const reco::Candidate *particle) const {
0150   unsigned int npart = particle->numberOfDaughters();
0151   if (!npart) {
0152     return 0;
0153   }
0154
0155   for (unsigned int i = 0; i < npart; ++i) {
0156     unsigned int idx = partIdx(p, particle->daughter(i));
0157     if (invalid[idx])
0158       continue;
0159     if (isParton((particle->daughter(i)->pdgId()))) {
0160       return 1;
0161     }
0162     if (isHadron((particle->daughter(i)->pdgId()))) {
0163       return -1;
0164     }
0165     int result = testPartonChildren(invalid, p, particle->daughter(i));
0166     if (result)
0167       return result;
0168   }
0169   return 0;
0170 }
0171
0172 InputGenJetsParticleSelector::ResonanceState InputGenJetsParticleSelector::fromResonance(
0173     ParticleBitmap &invalid, const ParticleVector &p, const reco::Candidate *particle) const {
0174   unsigned int idx = partIdx(p, particle);
0175   int id = particle->pdgId();
0176
0177   if (invalid[idx])
0178     return kIndirect;
0179
0180   if (isResonance(id) && (particle->status() == 3 || particle->status() == 22)) {
0181     return kDirect;
0182   }
0183
0184   if (!isIgnored(id) && (isParton(id)))
0185     return kNo;
0186
0187   unsigned int nMo = particle->numberOfMothers();
0188   if (!nMo)
0189     return kNo;
0190
0191   for (unsigned int i = 0; i < nMo; ++i) {
0192     ResonanceState result = fromResonance(invalid, p, particle->mother(i));
0193     switch (result) {
0194       case kNo:
0195         break;
0196       case kDirect:
0197         if (particle->mother(i)->pdgId() == id || isResonance(id))
0198           return kDirect;
0199         if (!isExcludedFromResonance(id))
0200           break;
0201         [[fallthrough]];
0202       case kIndirect:
0203         return kIndirect;
0204     }
0205   }
0206   return kNo;
0207 }
0208
0209 bool InputGenJetsParticleSelector::hasPartonChildren(ParticleBitmap &invalid,
0210                                                      const ParticleVector &p,
0211                                                      const reco::Candidate *particle) const {
0212   return testPartonChildren(invalid, p, particle) > 0;
0213 }
0214
0215 //######################################################
0216 //function NEEDED and called per EVENT by FRAMEWORK:
0217 void InputGenJetsParticleSelector::produce(edm::StreamID, edm::Event &evt, const edm::EventSetup &evtSetup) const {
0218   auto selected_ = std::make_unique<reco::CandidatePtrVector>();
0219
0220   ParticleVector particles;
0221
0222   edm::Handle<reco::CandidateView> prunedGenParticles;
0223   if (isMiniAOD) {
0224     evt.getByToken(input_prunedgenpartcoll_token_, prunedGenParticles);
0225
0226     for (edm::View<reco::Candidate>::const_iterator iter = prunedGenParticles->begin();
0227          iter != prunedGenParticles->end();
0228          ++iter) {
0229       if (iter->status() !=
0230           1)  // to avoid double-counting, skipping stable particles already contained in the collection of PackedGenParticles
0231         particles.push_back(&*iter);
0232     }
0233   }
0234
0235   edm::Handle<reco::CandidateView> genParticles;
0236   evt.getByToken(input_genpartcoll_token_, genParticles);
0237
0238   std::map<const reco::Candidate *, size_t> particlePtrIdxMap;
0239
0240   for (edm::View<reco::Candidate>::const_iterator iter = genParticles->begin(); iter != genParticles->end(); ++iter) {
0241     particles.push_back(&*iter);
0242     particlePtrIdxMap[&*iter] = (iter - genParticles->begin());
0243   }
0244
0245   std::sort(particles.begin(), particles.end());
0246   unsigned int size = particles.size();
0247
0248   ParticleBitmap selected(size, false);
0249   ParticleBitmap invalid(size, false);
0250
0251   for (unsigned int i = 0; i < size; i++) {
0252     const reco::Candidate *particle = particles[i];
0253     if (invalid[i])
0254       continue;
0255     if (particle->status() == 1)  // selecting stable particles
0256       selected[i] = true;
0257     if (partonicFinalState && isParton(particle->pdgId())) {
0258       if (particle->numberOfDaughters() == 0 && particle->status() != 1) {
0259         // some brokenness in event...
0260         invalid[i] = true;
0261       } else if (!hasPartonChildren(invalid, particles, particle)) {
0262         selected[i] = true;
0263         invalidateTree(invalid, particles, particle);  //this?!?
0264       }
0265     }
0266   }
0267
0268   for (size_t idx = 0; idx < size; ++idx) {
0269     const reco::Candidate *particle = particles[idx];
0270     if (!selected[idx] || invalid[idx]) {
0271       continue;
0272     }
0273
0274     if (excludeResonances && fromResonance(invalid, particles, particle)) {
0275       invalid[idx] = true;
0276       //cout<<"[INPUTSELECTOR] Invalidates FROM RESONANCE!: ["<<setw(4)<<idx<<"] "<<particle->pdgId()<<" "<<particle->pt()<<endl;
0277       continue;
0278     }
0279
0280     if (isIgnored(particle->pdgId())) {
0281       continue;
0282     }
0283
0284     if (particle->pt() >= ptMin) {
0285       selected_->push_back(genParticles->ptrAt(particlePtrIdxMap[particle]));
0286       //cout<<"Finally we have: ["<<setw(4)<<idx<<"] "<<setw(4)<<particle->pdgId()<<" "<<particle->pt()<<endl;
0287     }
0288   }
0289   evt.put(std::move(selected_));
0290 }
0291
0292 //define this as a plug-in
0293 DEFINE_FWK_MODULE(InputGenJetsParticleSelector);