File indexing completed on 2024-04-06 12:13:42
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028 #include <memory>
0029 #include <iostream>
0030
0031
0032 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0033
0034 #include "FWCore/Framework/interface/global/EDFilter.h"
0035 #include "FWCore/Framework/interface/Event.h"
0036 #include "FWCore/Framework/interface/Frameworkfwd.h"
0037 #include "FWCore/Framework/interface/MakerMacros.h"
0038 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0039
0040 #include "GeneratorInterface/GenFilters/plugins/MCFilterZboostHelper.h"
0041
0042 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0043
0044 using namespace edm;
0045 using namespace std;
0046
0047 namespace edm {
0048 class HepMCProduct;
0049 }
0050
0051 class PythiaFilterMultiAncestor : public edm::global::EDFilter<> {
0052 public:
0053 explicit PythiaFilterMultiAncestor(const edm::ParameterSet&);
0054
0055 bool filter(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0056
0057 private:
0058 bool isAncestor(HepMC::GenParticle* particle, int IDtoMatch) const;
0059
0060 const edm::EDGetTokenT<edm::HepMCProduct> token_;
0061 const int particleID;
0062 const double minpcut;
0063 const double maxpcut;
0064 const double minptcut;
0065 const double maxptcut;
0066 const double minetacut;
0067 const double maxetacut;
0068 const double minrapcut;
0069 const double maxrapcut;
0070 const double minphicut;
0071 const double maxphicut;
0072
0073 const int status;
0074 const std::vector<int> motherIDs;
0075 const std::vector<int> daughterIDs;
0076 const std::vector<double> daughterMinPts;
0077 const std::vector<double> daughterMaxPts;
0078 const std::vector<double> daughterMinEtas;
0079 const std::vector<double> daughterMaxEtas;
0080
0081 const int processID;
0082
0083 const double betaBoost;
0084 };
0085
0086 PythiaFilterMultiAncestor::PythiaFilterMultiAncestor(const edm::ParameterSet& iConfig)
0087 : token_(consumes<edm::HepMCProduct>(
0088 edm::InputTag(iConfig.getUntrackedParameter("moduleLabel", std::string("generator")), "unsmeared"))),
0089 particleID(iConfig.getUntrackedParameter("ParticleID", 0)),
0090 minpcut(iConfig.getUntrackedParameter("MinP", 0.)),
0091 maxpcut(iConfig.getUntrackedParameter("MaxP", 10000.)),
0092 minptcut(iConfig.getUntrackedParameter("MinPt", 0.)),
0093 maxptcut(iConfig.getUntrackedParameter("MaxPt", 10000.)),
0094 minetacut(iConfig.getUntrackedParameter("MinEta", -10.)),
0095 maxetacut(iConfig.getUntrackedParameter("MaxEta", 10.)),
0096 minrapcut(iConfig.getUntrackedParameter("MinRapidity", -20.)),
0097 maxrapcut(iConfig.getUntrackedParameter("MaxRapidity", 20.)),
0098 minphicut(iConfig.getUntrackedParameter("MinPhi", -3.5)),
0099 maxphicut(iConfig.getUntrackedParameter("MaxPhi", 3.5)),
0100 status(iConfig.getUntrackedParameter("Status", 0)),
0101 motherIDs(iConfig.getUntrackedParameter("MotherIDs", std::vector<int>{0})),
0102 daughterIDs(iConfig.getUntrackedParameter("DaughterIDs", std::vector<int>{0})),
0103 daughterMinPts(iConfig.getUntrackedParameter("DaughterMinPts", std::vector<double>{0.})),
0104 daughterMaxPts(iConfig.getUntrackedParameter("DaughterMaxPts", std::vector<double>{10000.})),
0105 daughterMinEtas(iConfig.getUntrackedParameter("DaughterMinEtas", std::vector<double>{-10.})),
0106 daughterMaxEtas(iConfig.getUntrackedParameter("DaughterMaxEtas", std::vector<double>{10.})),
0107 processID(iConfig.getUntrackedParameter("ProcessID", 0)),
0108 betaBoost(iConfig.getUntrackedParameter("BetaBoost", 0.)) {
0109
0110 }
0111
0112
0113 bool PythiaFilterMultiAncestor::isAncestor(HepMC::GenParticle* particle, int IDtoMatch) const {
0114 for (HepMC::GenVertex::particle_iterator ancestor = particle->production_vertex()->particles_begin(HepMC::ancestors);
0115 ancestor != particle->production_vertex()->particles_end(HepMC::ancestors);
0116 ++ancestor) {
0117
0118
0119
0120
0121 if (abs((*ancestor)->pdg_id()) == abs(IDtoMatch)) {
0122
0123 return true;
0124 }
0125 }
0126
0127
0128 return false;
0129 }
0130
0131
0132 bool PythiaFilterMultiAncestor::filter(edm::StreamID, edm::Event& iEvent, const edm::EventSetup&) const {
0133 using namespace edm;
0134 bool accepted = false;
0135 Handle<HepMCProduct> evt;
0136 iEvent.getByToken(token_, evt);
0137
0138 const HepMC::GenEvent* myGenEvent = evt->GetEvent();
0139
0140 if (processID == 0 || processID == myGenEvent->signal_process_id()) {
0141 for (HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end();
0142 ++p) {
0143 HepMC::FourVector mom = MCFilterZboostHelper::zboost((*p)->momentum(), betaBoost);
0144 double rapidity = 0.5 * log((mom.e() + mom.pz()) / (mom.e() - mom.pz()));
0145
0146 if (abs((*p)->pdg_id()) == particleID && mom.rho() > minpcut && mom.rho() < maxpcut &&
0147 (*p)->momentum().perp() > minptcut && (*p)->momentum().perp() < maxptcut && mom.eta() > minetacut &&
0148 mom.eta() < maxetacut && rapidity > minrapcut && rapidity < maxrapcut && (*p)->momentum().phi() > minphicut &&
0149 (*p)->momentum().phi() < maxphicut) {
0150
0151 for (std::vector<int>::const_iterator motherID = motherIDs.begin(); motherID != motherIDs.end(); ++motherID) {
0152
0153 if (status == 0 && *motherID == 0) {
0154 accepted = true;
0155 }
0156 if (status != 0 && *motherID == 0) {
0157 if ((*p)->status() == status)
0158 accepted = true;
0159 }
0160
0161
0162 if (status == 0 && *motherID != 0) {
0163
0164 if (isAncestor(*p, *motherID)) {
0165 accepted = true;
0166 }
0167 }
0168 if (status != 0 && *motherID != 0) {
0169
0170 if ((*p)->status() == status && isAncestor(*p, *motherID)) {
0171 accepted = true;
0172 }
0173 }
0174 }
0175
0176
0177 if (accepted & (!daughterIDs.empty())) {
0178
0179
0180
0181
0182 uint good_dau = 0;
0183 uint good_dau_cc = 0;
0184 for (HepMC::GenVertex::particle_iterator dau = (*p)->end_vertex()->particles_begin(HepMC::children);
0185 dau != (*p)->end_vertex()->particles_end(HepMC::children);
0186 ++dau) {
0187 for (unsigned int i = 0; i < daughterIDs.size(); ++i) {
0188
0189
0190 if ((*dau)->pdg_id() == daughterIDs[i]) {
0191 if ((*dau)->momentum().perp() < daughterMinPts[i])
0192 continue;
0193 if ((*dau)->momentum().perp() > daughterMaxPts[i])
0194 continue;
0195 if ((*dau)->momentum().eta() < daughterMinEtas[i])
0196 continue;
0197 if ((*dau)->momentum().eta() > daughterMaxEtas[i])
0198 continue;
0199 ++good_dau;
0200 }
0201
0202 if (-(*dau)->pdg_id() == daughterIDs[i]) {
0203 if ((*dau)->momentum().perp() < daughterMinPts[i])
0204 continue;
0205 if ((*dau)->momentum().perp() > daughterMaxPts[i])
0206 continue;
0207 if ((*dau)->momentum().eta() < daughterMinEtas[i])
0208 continue;
0209 if ((*dau)->momentum().eta() > daughterMaxEtas[i])
0210 continue;
0211 ++good_dau_cc;
0212 }
0213 }
0214 }
0215 if (good_dau < daughterIDs.size() && good_dau_cc < daughterIDs.size())
0216 accepted = false;
0217 }
0218 }
0219
0220 if (accepted)
0221 break;
0222 }
0223
0224 } else {
0225 accepted = true;
0226 }
0227
0228 if (accepted) {
0229 return true;
0230 } else {
0231 return false;
0232 }
0233 }
0234
0235 DEFINE_FWK_MODULE(PythiaFilterMultiAncestor);