GenFilters/plugins/MCSmartSingleParticleFilter.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    MCSmartSingleParticleFilter
0004 // Class:      MCSmartSingleParticleFilter
0005 //
0006 /*
0007
0008  Description: filter events based on the Pythia particleID, the Pt and the production vertex
0009
0010  Implementation: inherits from generic EDFilter
0011
0012 */
0013 //         Created:  J. Alcaraz, 04/07/2008
0014 //
0015
0016 #include "DataFormats/Common/interface/Handle.h"
0017 #include "FWCore/Framework/interface/Event.h"
0018 #include "FWCore/Framework/interface/Frameworkfwd.h"
0019 #include "FWCore/Framework/interface/global/EDFilter.h"
0020 #include "FWCore/Framework/interface/MakerMacros.h"
0021 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0022 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0023 #include "FWCore/Utilities/interface/EDGetToken.h"
0024 #include "FWCore/Utilities/interface/InputTag.h"
0025 #include "GeneratorInterface/GenFilters/plugins/MCFilterZboostHelper.h"
0026 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0027
0028 #include <cmath>
0029 #include <cstdlib>
0030 #include <iostream>
0031 #include <string>
0032 #include <vector>
0033
0034 class MCSmartSingleParticleFilter : public edm::global::EDFilter<> {
0035 public:
0036   explicit MCSmartSingleParticleFilter(const edm::ParameterSet&);
0037
0038   bool filter(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0039
0040 private:
0041   const edm::EDGetTokenT<edm::HepMCProduct> token_;
0042   std::vector<int> particleID;
0043   std::vector<double> pMin;
0044   std::vector<double> ptMin;
0045   std::vector<double> etaMin;
0046   std::vector<double> etaMax;
0047   std::vector<int> status;
0048   std::vector<double> decayRadiusMin;
0049   std::vector<double> decayRadiusMax;
0050   std::vector<double> decayZMin;
0051   std::vector<double> decayZMax;
0052   const double betaBoost;
0053 };
0054
0055 using namespace std;
0056
0057 MCSmartSingleParticleFilter::MCSmartSingleParticleFilter(const edm::ParameterSet& iConfig)
0058     : token_(consumes<edm::HepMCProduct>(
0059           iConfig.getUntrackedParameter<edm::InputTag>("moduleLabel", edm::InputTag("generator", "unsmeared")))),
0060       betaBoost(iConfig.getUntrackedParameter("BetaBoost", 0.)) {
0061   vector<int> defpid;
0062   defpid.push_back(0);
0063   particleID = iConfig.getUntrackedParameter<vector<int> >("ParticleID", defpid);
0064   vector<double> defpmin;
0065   defpmin.push_back(0.);
0066   pMin = iConfig.getUntrackedParameter<vector<double> >("MinP", defpmin);
0067
0068   vector<double> defptmin;
0069   defptmin.push_back(0.);
0070   ptMin = iConfig.getUntrackedParameter<vector<double> >("MinPt", defptmin);
0071
0072   vector<double> defetamin;
0073   defetamin.push_back(-10.);
0074   etaMin = iConfig.getUntrackedParameter<vector<double> >("MinEta", defetamin);
0075   vector<double> defetamax;
0076   defetamax.push_back(10.);
0077   etaMax = iConfig.getUntrackedParameter<vector<double> >("MaxEta", defetamax);
0078   vector<int> defstat;
0079   defstat.push_back(0);
0080   status = iConfig.getUntrackedParameter<vector<int> >("Status", defstat);
0081
0082   vector<double> defDecayRadiusmin;
0083   defDecayRadiusmin.push_back(-1.);
0084   decayRadiusMin = iConfig.getUntrackedParameter<vector<double> >("MinDecayRadius", defDecayRadiusmin);
0085
0086   vector<double> defDecayRadiusmax;
0087   defDecayRadiusmax.push_back(1.e5);
0088   decayRadiusMax = iConfig.getUntrackedParameter<vector<double> >("MaxDecayRadius", defDecayRadiusmax);
0089
0090   vector<double> defDecayZmin;
0091   defDecayZmin.push_back(-1.e5);
0092   decayZMin = iConfig.getUntrackedParameter<vector<double> >("MinDecayZ", defDecayZmin);
0093
0094   vector<double> defDecayZmax;
0095   defDecayZmax.push_back(1.e5);
0096   decayZMax = iConfig.getUntrackedParameter<vector<double> >("MaxDecayZ", defDecayZmax);
0097
0098   // check for same size
0099   if ((pMin.size() > 1 && particleID.size() != pMin.size()) ||
0100       (ptMin.size() > 1 && particleID.size() != ptMin.size()) ||
0101       (etaMin.size() > 1 && particleID.size() != etaMin.size()) ||
0102       (etaMax.size() > 1 && particleID.size() != etaMax.size()) ||
0103       (status.size() > 1 && particleID.size() != status.size()) ||
0104       (decayRadiusMin.size() > 1 && particleID.size() != decayRadiusMin.size()) ||
0105       (decayRadiusMax.size() > 1 && particleID.size() != decayRadiusMax.size()) ||
0106       (decayZMin.size() > 1 && particleID.size() != decayZMin.size()) ||
0107       (decayZMax.size() > 1 && particleID.size() != decayZMax.size())) {
0108     edm::LogError("Configuration")
0109         << "WARNING: MCPROCESSFILTER : size of MinPthat and/or MaxPthat not matching with ProcessID size!!";
0110   }
0111
0112   // if pMin size smaller than particleID , fill up further with defaults
0113   if (particleID.size() > pMin.size()) {
0114     for (unsigned int i = pMin.size(); i < particleID.size(); i++) {
0115       pMin.push_back(0.);
0116     }
0117   }
0118   // if ptMin size smaller than particleID , fill up further with defaults
0119   if (particleID.size() > ptMin.size()) {
0120     for (unsigned int i = ptMin.size(); i < particleID.size(); i++) {
0121       ptMin.push_back(0.);
0122     }
0123   }
0124   // if etaMin size smaller than particleID , fill up further with defaults
0125   if (particleID.size() > etaMin.size()) {
0126     for (unsigned int i = etaMin.size(); i < particleID.size(); i++) {
0127       etaMin.push_back(-10.);
0128     }
0129   }
0130   // if etaMax size smaller than particleID , fill up further with defaults
0131   if (particleID.size() > etaMax.size()) {
0132     for (unsigned int i = etaMax.size(); i < particleID.size(); i++) {
0133       etaMax.push_back(10.);
0134     }
0135   }
0136   // if status size smaller than particleID , fill up further with defaults
0137   if (particleID.size() > status.size()) {
0138     for (unsigned int i = status.size(); i < particleID.size(); i++) {
0139       status.push_back(0);
0140     }
0141   }
0142
0143   // if decayRadiusMin size smaller than particleID , fill up further with defaults
0144   if (particleID.size() > decayRadiusMin.size()) {
0145     for (unsigned int i = decayRadiusMin.size(); i < particleID.size(); i++) {
0146       decayRadiusMin.push_back(-10.);
0147     }
0148   }
0149   // if decayRadiusMax size smaller than particleID , fill up further with defaults
0150   if (particleID.size() > decayRadiusMax.size()) {
0151     for (unsigned int i = decayRadiusMax.size(); i < particleID.size(); i++) {
0152       decayRadiusMax.push_back(1.e5);
0153     }
0154   }
0155
0156   // if decayZMin size smaller than particleID , fill up further with defaults
0157   if (particleID.size() > decayZMin.size()) {
0158     for (unsigned int i = decayZMin.size(); i < particleID.size(); i++) {
0159       decayZMin.push_back(-1.e5);
0160     }
0161   }
0162   // if decayZMax size smaller than particleID , fill up further with defaults
0163   if (particleID.size() > decayZMax.size()) {
0164     for (unsigned int i = decayZMax.size(); i < particleID.size(); i++) {
0165       decayZMax.push_back(1.e5);
0166     }
0167   }
0168
0169   // check if beta is smaller than 1
0170   if (std::abs(betaBoost) >= 1) {
0171     edm::LogError("MCSmartSingleParticleFilter") << "Input beta boost is >= 1 !";
0172   }
0173 }
0174
0175 bool MCSmartSingleParticleFilter::filter(edm::StreamID, edm::Event& iEvent, const edm::EventSetup&) const {
0176   bool accepted = false;
0177   edm::Handle<edm::HepMCProduct> evt;
0178   iEvent.getByToken(token_, evt);
0179
0180   const HepMC::GenEvent* myGenEvent = evt->GetEvent();
0181
0182   for (HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end();
0183        ++p) {
0184     for (unsigned int i = 0; i < particleID.size(); i++) {
0185       if (particleID[i] == (*p)->pdg_id() || particleID[i] == 0) {
0186         if ((*p)->momentum().perp() > ptMin[i] && ((*p)->status() == status[i] || status[i] == 0)) {
0187           HepMC::FourVector mom = MCFilterZboostHelper::zboost((*p)->momentum(), betaBoost);
0188           if (mom.rho() > pMin[i] && mom.eta() > etaMin[i] && mom.eta() < etaMax[i]) {
0189             if (!((*p)->production_vertex()))
0190               continue;
0191
0192             double decx = (*p)->production_vertex()->position().x();
0193             double decy = (*p)->production_vertex()->position().y();
0194             double decrad = sqrt(decx * decx + decy * decy);
0195             if (decrad < decayRadiusMin[i])
0196               continue;
0197             if (decrad > decayRadiusMax[i])
0198               continue;
0199
0200             double decz = (*p)->production_vertex()->position().z();
0201             if (decz < decayZMin[i])
0202               continue;
0203             if (decz > decayZMax[i])
0204               continue;
0205
0206             accepted = true;
0207           }
0208         }
0209       }
0210     }
0211   }
0212   return accepted;
0213 }
0214
0215 DEFINE_FWK_MODULE(MCSmartSingleParticleFilter);