File indexing completed on 2023-10-25 09:48:27
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
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
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
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
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
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
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
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
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
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
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
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
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);