GenFilters/plugins/MCMultiParticleFilter.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    MCMultiParticleFilter
0004 // Class:      MCMultiParticleFilter
0005 //
0006 /*
0007
0008  Description: Filter to select events with an arbitrary number of given particle(s).
0009
0010  Implementation: derived from MCSingleParticleFilter
0011
0012 */
0013 //
0014 // Original Author:  Paul Lujan
0015 //         Created:  Wed Feb 29 04:22:16 CST 2012
0016 //
0017 //
0018
0019 #include "DataFormats/Common/interface/Handle.h"
0020 #include "FWCore/Framework/interface/global/EDFilter.h"
0021 #include "FWCore/Framework/interface/Event.h"
0022 #include "FWCore/Framework/interface/Frameworkfwd.h"
0023 #include "FWCore/Framework/interface/MakerMacros.h"
0024 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0025 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0026 #include "FWCore/Utilities/interface/EDGetToken.h"
0027 #include "FWCore/Utilities/interface/InputTag.h"
0028 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0029
0030 #include <cmath>
0031 #include <cstdlib>
0032 #include <vector>
0033
0034 //
0035 // class declaration
0036 //
0037
0038 class MCMultiParticleFilter : public edm::global::EDFilter<> {
0039 public:
0040   explicit MCMultiParticleFilter(const edm::ParameterSet&);
0041
0042 private:
0043   bool filter(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0044
0045   // ----------member data ---------------------------
0046
0047   const edm::EDGetTokenT<edm::HepMCProduct> token_;
0048   const int numRequired_;              // number of particles required to pass filter
0049   const bool acceptMore_;              // if true (default), accept numRequired or more.
0050                                        // if false, accept events with exactly equal to numRequired.
0051   const std::vector<int> particleID_;  // vector of particle IDs to look for
0052   // the four next variables can either be a vector of length 1 (in which case the same
0053   // value is used for all particle IDs) or of length equal to the length of ParticleID (in which
0054   // case the corresponding value is used for each).
0055   std::vector<int> motherID_;   // mother ID of particles (optional)
0056   std::vector<double> ptMin_;   // minimum Pt of particles
0057   std::vector<double> etaMax_;  // maximum fabs(eta) of particles
0058   std::vector<int> status_;     // status of particles
0059   std::vector<double> decayRadiusMin;
0060   std::vector<double> decayRadiusMax;
0061   std::vector<double> decayZMin;
0062   std::vector<double> decayZMax;
0063 };
0064
0065 MCMultiParticleFilter::MCMultiParticleFilter(const edm::ParameterSet& iConfig)
0066     : token_(consumes<edm::HepMCProduct>(
0067           iConfig.getUntrackedParameter<edm::InputTag>("src", edm::InputTag("generator", "unsmeared")))),
0068       numRequired_(iConfig.getParameter<int>("NumRequired")),
0069       acceptMore_(iConfig.getParameter<bool>("AcceptMore")),
0070       particleID_(iConfig.getParameter<std::vector<int> >("ParticleID")),
0071       ptMin_(iConfig.getParameter<std::vector<double> >("PtMin")),
0072       etaMax_(iConfig.getParameter<std::vector<double> >("EtaMax")),
0073       status_(iConfig.getParameter<std::vector<int> >("Status")) {
0074   //here do whatever other initialization is needed
0075
0076   // default pt, eta, status cuts to "don't care"
0077   std::vector<double> defptmin(1, 0);
0078   std::vector<double> defetamax(1, 999.0);
0079   std::vector<int> defstat(1, 0);
0080   std::vector<int> defmother;
0081   defmother.push_back(0);
0082   motherID_ = iConfig.getUntrackedParameter<std::vector<int> >("MotherID", defmother);
0083
0084   std::vector<double> defDecayRadiusmin;
0085   defDecayRadiusmin.push_back(-1.);
0086   decayRadiusMin = iConfig.getUntrackedParameter<std::vector<double> >("MinDecayRadius", defDecayRadiusmin);
0087
0088   std::vector<double> defDecayRadiusmax;
0089   defDecayRadiusmax.push_back(1.e5);
0090   decayRadiusMax = iConfig.getUntrackedParameter<std::vector<double> >("MaxDecayRadius", defDecayRadiusmax);
0091
0092   std::vector<double> defDecayZmin;
0093   defDecayZmin.push_back(-1.e5);
0094   decayZMin = iConfig.getUntrackedParameter<std::vector<double> >("MinDecayZ", defDecayZmin);
0095
0096   std::vector<double> defDecayZmax;
0097   defDecayZmax.push_back(1.e5);
0098   decayZMax = iConfig.getUntrackedParameter<std::vector<double> >("MaxDecayZ", defDecayZmax);
0099
0100   // check for same size
0101   if ((ptMin_.size() > 1 && particleID_.size() != ptMin_.size()) ||
0102       (etaMax_.size() > 1 && particleID_.size() != etaMax_.size()) ||
0103       (status_.size() > 1 && particleID_.size() != status_.size()) ||
0104       (motherID_.size() > 1 && particleID_.size() != motherID_.size()) ||
0105       (decayRadiusMin.size() > 1 && particleID_.size() != decayRadiusMin.size()) ||
0106       (decayRadiusMax.size() > 1 && particleID_.size() != decayRadiusMax.size()) ||
0107       (decayZMin.size() > 1 && particleID_.size() != decayZMin.size()) ||
0108       (decayZMax.size() > 1 && particleID_.size() != decayZMax.size())) {
0109     edm::LogWarning("MCMultiParticleFilter") << "WARNING: MCMultiParticleFilter: size of PtMin, EtaMax, motherID, "
0110                                                 "decayRadiusMin, decayRadiusMax, decayZMin, decayZMax"
0111                                                 "and/or Status does not match ParticleID size!"
0112                                              << std::endl;
0113   }
0114
0115   // Fill arrays with defaults if necessary
0116   while (ptMin_.size() < particleID_.size())
0117     ptMin_.push_back(defptmin[0]);
0118   while (etaMax_.size() < particleID_.size())
0119     etaMax_.push_back(defetamax[0]);
0120   while (status_.size() < particleID_.size())
0121     status_.push_back(defstat[0]);
0122   while (motherID_.size() < particleID_.size())
0123     motherID_.push_back(defmother[0]);
0124
0125   // if decayRadiusMin size smaller than particleID , fill up further with defaults
0126   if (particleID_.size() > decayRadiusMin.size()) {
0127     for (unsigned int i = decayRadiusMin.size(); i < particleID_.size(); i++) {
0128       decayRadiusMin.push_back(-10.);
0129     }
0130   }
0131   // if decayRadiusMax size smaller than particleID , fill up further with defaults
0132   if (particleID_.size() > decayRadiusMax.size()) {
0133     for (unsigned int i = decayRadiusMax.size(); i < particleID_.size(); i++) {
0134       decayRadiusMax.push_back(1.e5);
0135     }
0136   }
0137
0138   // if decayZMin size smaller than particleID , fill up further with defaults
0139   if (particleID_.size() > decayZMin.size()) {
0140     for (unsigned int i = decayZMin.size(); i < particleID_.size(); i++) {
0141       decayZMin.push_back(-1.e5);
0142     }
0143   }
0144   // if decayZMax size smaller than particleID , fill up further with defaults
0145   if (particleID_.size() > decayZMax.size()) {
0146     for (unsigned int i = decayZMax.size(); i < particleID_.size(); i++) {
0147       decayZMax.push_back(1.e5);
0148     }
0149   }
0150 }
0151
0152 // ------------ method called to skim the data  ------------
0153 bool MCMultiParticleFilter::filter(edm::StreamID, edm::Event& iEvent, const edm::EventSetup&) const {
0154   edm::Handle<edm::HepMCProduct> evt;
0155   iEvent.getByToken(token_, evt);
0156
0157   int nFound = 0;
0158
0159   const HepMC::GenEvent* myGenEvent = evt->GetEvent();
0160
0161   for (HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end();
0162        ++p) {
0163     for (unsigned int i = 0; i < particleID_.size(); ++i) {
0164       if ((particleID_[i] == 0 || std::abs(particleID_[i]) == std::abs((*p)->pdg_id())) &&
0165           (*p)->momentum().perp() > ptMin_[i] && std::fabs((*p)->momentum().eta()) < etaMax_[i] &&
0166           (status_[i] == 0 || (*p)->status() == status_[i])) {
0167         if (!((*p)->production_vertex()))
0168           continue;
0169
0170         double decx = (*p)->production_vertex()->position().x();
0171         double decy = (*p)->production_vertex()->position().y();
0172         double decrad = sqrt(decx * decx + decy * decy);
0173         if (decrad < decayRadiusMin[i])
0174           continue;
0175         if (decrad > decayRadiusMax[i])
0176           continue;
0177
0178         double decz = (*p)->production_vertex()->position().z();
0179         if (decz < decayZMin[i])
0180           continue;
0181         if (decz > decayZMax[i])
0182           continue;
0183
0184         if (motherID_[i] == 0) {  // do not check for mother ID if not sepcified
0185           nFound++;
0186           break;  // only match a given particle once!
0187         } else {
0188           bool hascorrectmother = false;
0189           for (HepMC::GenVertex::particles_in_const_iterator mo = (*p)->production_vertex()->particles_in_const_begin();
0190                mo != (*p)->production_vertex()->particles_in_const_end();
0191                ++mo) {
0192             if ((*mo)->pdg_id() == motherID_[i]) {
0193               hascorrectmother = true;
0194               break;
0195             }
0196           }
0197           if (hascorrectmother) {
0198             nFound++;
0199             break;  // only match a given particle once!
0200           }
0201         }
0202       }
0203     }  // loop over targets
0204
0205     if (acceptMore_ && nFound == numRequired_)
0206       break;  // stop looking if we don't mind having more
0207   }           // loop over particles
0208
0209   if (nFound == numRequired_) {
0210     return true;
0211   } else {
0212     return false;
0213   }
0214 }
0215
0216 //define this as a plug-in
0217 DEFINE_FWK_MODULE(MCMultiParticleFilter);