File indexing completed on 2024-04-06 12:13:41
0001 #include "GeneratorInterface/GenFilters/plugins/PythiaDauVFilter.h"
0002 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0003 #include <iostream>
0004 #include <memory>
0005
0006 #include <vector>
0007
0008 using namespace edm;
0009 using namespace std;
0010 using namespace Pythia8;
0011
0012 PythiaDauVFilter::PythiaDauVFilter(const edm::ParameterSet& iConfig)
0013 : fVerbose(iConfig.getUntrackedParameter("verbose", 0)),
0014 token_(consumes<edm::HepMCProduct>(
0015 edm::InputTag(iConfig.getUntrackedParameter("moduleLabel", std::string("generator")), "unsmeared"))),
0016 particleID(iConfig.getUntrackedParameter("ParticleID", 0)),
0017 motherID(iConfig.getUntrackedParameter("MotherID", 0)),
0018 chargeconju(iConfig.getUntrackedParameter("ChargeConjugation", true)),
0019 ndaughters(iConfig.getUntrackedParameter("NumberDaughters", 0)),
0020 maxptcut(iConfig.getUntrackedParameter("MaxPt", 14000.)) {
0021
0022 vector<int> defdauID;
0023 defdauID.push_back(0);
0024 dauIDs = iConfig.getUntrackedParameter<vector<int> >("DaughterIDs", defdauID);
0025 vector<double> defminptcut;
0026 defminptcut.push_back(0.);
0027 minptcut = iConfig.getUntrackedParameter<vector<double> >("MinPt", defminptcut);
0028 vector<double> defminetacut;
0029 defminetacut.push_back(-10.);
0030 minetacut = iConfig.getUntrackedParameter<vector<double> >("MinEta", defminetacut);
0031 vector<double> defmaxetacut;
0032 defmaxetacut.push_back(10.);
0033 maxetacut = iConfig.getUntrackedParameter<vector<double> >("MaxEta", defmaxetacut);
0034
0035 edm::LogInfo("PythiaDauVFilter") << "----------------------------------------------------------------------" << endl;
0036 edm::LogInfo("PythiaDauVFilter") << "--- PythiaDauVFilter" << endl;
0037 for (unsigned int i = 0; i < dauIDs.size(); ++i) {
0038 edm::LogInfo("PythiaDauVFilter") << "ID: " << dauIDs[i] << " pT > " << minptcut[i] << " " << minetacut[i]
0039 << " eta < " << maxetacut[i] << endl;
0040 }
0041 edm::LogInfo("PythiaDauVFilter") << "maxptcut = " << maxptcut << endl;
0042 edm::LogInfo("PythiaDauVFilter") << "particleID = " << particleID << endl;
0043 edm::LogInfo("PythiaDauVFilter") << "motherID = " << motherID << endl;
0044
0045
0046 edm::LogInfo("PythiaDauVFilter") << "Creating pythia8 instance for particle properties" << endl;
0047 if (!fLookupGen.get())
0048 fLookupGen = std::make_unique<Pythia>();
0049 }
0050
0051 PythiaDauVFilter::~PythiaDauVFilter() {
0052
0053
0054 }
0055
0056
0057
0058
0059
0060
0061 bool PythiaDauVFilter::filter(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
0062 using namespace edm;
0063 bool accepted = false;
0064 Handle<HepMCProduct> evt;
0065 iEvent.getByToken(token_, evt);
0066
0067 int OK(1);
0068 vector<int> vparticles;
0069
0070 HepMC::GenEvent* myGenEvent = new HepMC::GenEvent(*(evt->GetEvent()));
0071
0072 if (fVerbose > 5) {
0073 edm::LogInfo("PythiaDauVFilter") << "looking for " << particleID << endl;
0074 }
0075
0076 for (HepMC::GenEvent::particle_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end(); ++p) {
0077 if ((*p)->pdg_id() != particleID)
0078 continue;
0079
0080
0081 if (0 != motherID) {
0082 OK = 0;
0083 for (HepMC::GenVertex::particles_in_const_iterator des = (*p)->production_vertex()->particles_in_const_begin();
0084 des != (*p)->production_vertex()->particles_in_const_end();
0085 ++des) {
0086 if (fVerbose > 10) {
0087 edm::LogInfo("PythiaDauVFilter") << "mother: " << (*des)->pdg_id() << " pT: " << (*des)->momentum().perp()
0088 << " eta: " << (*des)->momentum().eta() << endl;
0089 }
0090 if (abs(motherID) == abs((*des)->pdg_id())) {
0091 OK = 1;
0092 break;
0093 }
0094 }
0095 }
0096 if (0 == OK)
0097 continue;
0098
0099
0100 int ndauac = 0;
0101 int ndau = 0;
0102 if (fVerbose > 5) {
0103 edm::LogInfo("PythiaDauVFilter") << "found ID: " << (*p)->pdg_id() << " pT: " << (*p)->momentum().perp()
0104 << " eta: " << (*p)->momentum().eta() << endl;
0105 }
0106 vparticles.push_back((*p)->pdg_id());
0107 if ((*p)->end_vertex()) {
0108 for (HepMC::GenVertex::particle_iterator des = (*p)->end_vertex()->particles_begin(HepMC::children);
0109 des != (*p)->end_vertex()->particles_end(HepMC::children);
0110 ++des) {
0111 ++ndau;
0112 if (fVerbose > 5) {
0113 edm::LogInfo("PythiaDauVFilter") << "ID: " << (*des)->pdg_id() << " pT: " << (*des)->momentum().perp()
0114 << " eta: " << (*des)->momentum().eta() << endl;
0115 }
0116 for (unsigned int i = 0; i < dauIDs.size(); ++i) {
0117 if ((*des)->pdg_id() != dauIDs[i])
0118 continue;
0119 if (fVerbose > 5) {
0120 edm::LogInfo("PythiaDauVFilter") << "i = " << i << " pT = " << (*des)->momentum().perp()
0121 << " eta = " << (*des)->momentum().eta() << endl;
0122 }
0123 if ((*des)->momentum().perp() > minptcut[i] && (*des)->momentum().perp() < maxptcut &&
0124 (*des)->momentum().eta() > minetacut[i] && (*des)->momentum().eta() < maxetacut[i]) {
0125 ++ndauac;
0126 vparticles.push_back((*des)->pdg_id());
0127 if (fVerbose > 2) {
0128 edm::LogInfo("PythiaDauVFilter")
0129 << " accepted this particle " << (*des)->pdg_id() << " pT = " << (*des)->momentum().perp()
0130 << " eta = " << (*des)->momentum().eta() << endl;
0131 }
0132 break;
0133 }
0134 }
0135 }
0136 }
0137
0138
0139 if (ndau >= ndaughters && ndauac == ndaughters) {
0140 accepted = true;
0141 if (fVerbose > 0) {
0142 edm::LogInfo("PythiaDauVFilter") << " accepted this decay: ";
0143 for (unsigned int iv = 0; iv < vparticles.size(); ++iv)
0144 edm::LogInfo("PythiaDauVFilter") << vparticles[iv] << " ";
0145 edm::LogInfo("PythiaDauVFilter") << " from mother = " << motherID << endl;
0146 }
0147 break;
0148 }
0149 }
0150
0151 if (!accepted && chargeconju) {
0152 OK = 1;
0153
0154 for (HepMC::GenEvent::particle_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end(); ++p) {
0155 if ((*p)->pdg_id() != -particleID)
0156 continue;
0157
0158
0159 if (0 != motherID) {
0160 OK = 0;
0161 for (HepMC::GenVertex::particles_in_const_iterator des = (*p)->production_vertex()->particles_in_const_begin();
0162 des != (*p)->production_vertex()->particles_in_const_end();
0163 ++des) {
0164 if (fVerbose > 10) {
0165 edm::LogInfo("PythiaDauVFilter") << "mother: " << (*des)->pdg_id() << " pT: " << (*des)->momentum().perp()
0166 << " eta: " << (*des)->momentum().eta() << endl;
0167 }
0168 if (abs(motherID) == abs((*des)->pdg_id())) {
0169 OK = 1;
0170 break;
0171 }
0172 }
0173 }
0174 if (0 == OK)
0175 continue;
0176
0177 if (fVerbose > 5) {
0178 edm::LogInfo("PythiaDauVFilter") << "found ID: " << (*p)->pdg_id() << " pT: " << (*p)->momentum().perp()
0179 << " eta: " << (*p)->momentum().eta() << endl;
0180 }
0181 vparticles.push_back((*p)->pdg_id());
0182 int ndauac = 0;
0183 int ndau = 0;
0184 if ((*p)->end_vertex()) {
0185 for (HepMC::GenVertex::particle_iterator des = (*p)->end_vertex()->particles_begin(HepMC::children);
0186 des != (*p)->end_vertex()->particles_end(HepMC::children);
0187 ++des) {
0188 ++ndau;
0189 if (fVerbose > 5) {
0190 edm::LogInfo("PythiaDauVFilter") << "ID: " << (*des)->pdg_id() << " pT: " << (*des)->momentum().perp()
0191 << " eta: " << (*des)->momentum().eta() << endl;
0192 }
0193 for (unsigned int i = 0; i < dauIDs.size(); ++i) {
0194 int IDanti = -dauIDs[i];
0195 if (!(fLookupGen->particleData.isParticle(IDanti)))
0196 IDanti = dauIDs[i];
0197 if ((*des)->pdg_id() != IDanti)
0198 continue;
0199 if (fVerbose > 5) {
0200 edm::LogInfo("PythiaDauVFilter") << "i = " << i << " pT = " << (*des)->momentum().perp()
0201 << " eta = " << (*des)->momentum().eta() << endl;
0202 }
0203 if ((*des)->momentum().perp() > minptcut[i] && (*des)->momentum().perp() < maxptcut &&
0204 (*des)->momentum().eta() > minetacut[i] && (*des)->momentum().eta() < maxetacut[i]) {
0205 ++ndauac;
0206 vparticles.push_back((*des)->pdg_id());
0207 if (fVerbose > 2) {
0208 edm::LogInfo("PythiaDauVFilter")
0209 << " accepted this particle " << (*des)->pdg_id() << " pT = " << (*des)->momentum().perp()
0210 << " eta = " << (*des)->momentum().eta() << endl;
0211 }
0212 break;
0213 }
0214 }
0215 }
0216 }
0217 if (ndau >= ndaughters && ndauac == ndaughters) {
0218 accepted = true;
0219 if (fVerbose > 0) {
0220 edm::LogInfo("PythiaDauVFilter") << " accepted this anti-decay: ";
0221 for (unsigned int iv = 0; iv < vparticles.size(); ++iv)
0222 edm::LogInfo("PythiaDauVFilter") << vparticles[iv] << " ";
0223 edm::LogInfo("PythiaDauVFilter") << " from mother = " << motherID << endl;
0224 }
0225 break;
0226 }
0227 }
0228 }
0229
0230 delete myGenEvent;
0231 return accepted;
0232 }
0233
0234
0235 DEFINE_FWK_MODULE(PythiaDauVFilter);