File indexing completed on 2023-03-17 11:04:34
0001 #include "GeneratorInterface/GenFilters/plugins/PythiaDauVFilterMatchID.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 PythiaDauVFilterMatchID::PythiaDauVFilterMatchID(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("PythiaDauVFilterMatchID")
0036 << "----------------------------------------------------------------------" << endl;
0037 edm::LogInfo("PythiaDauVFilterMatchID") << "--- PythiaDauVFilterMatchID" << endl;
0038 for (unsigned int i = 0; i < dauIDs.size(); ++i) {
0039 edm::LogInfo("PythiaDauVFilterMatchID")
0040 << "ID: " << dauIDs[i] << " pT > " << minptcut[i] << " " << minetacut[i] << " eta < " << maxetacut[i] << endl;
0041 }
0042 edm::LogInfo("PythiaDauVFilterMatchID") << "maxptcut = " << maxptcut << endl;
0043 edm::LogInfo("PythiaDauVFilterMatchID") << "particleID = " << particleID << endl;
0044 edm::LogInfo("PythiaDauVFilterMatchID") << "motherID = " << motherID << endl;
0045
0046
0047 edm::LogInfo("PythiaDauVFilterMatchID") << "Creating pythia8 instance for particle properties" << endl;
0048 if (!fLookupGen.get())
0049 fLookupGen = std::make_unique<Pythia>();
0050 }
0051
0052 PythiaDauVFilterMatchID::~PythiaDauVFilterMatchID() {
0053
0054
0055 }
0056
0057
0058
0059
0060
0061
0062 bool PythiaDauVFilterMatchID::filter(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
0063 using namespace edm;
0064 bool accepted = false;
0065 Handle<HepMCProduct> evt;
0066 iEvent.getByToken(token_, evt);
0067
0068 int OK(1);
0069 vector<int> vparticles;
0070
0071 HepMC::GenEvent* myGenEvent = new HepMC::GenEvent(*(evt->GetEvent()));
0072
0073 if (fVerbose > 5) {
0074 edm::LogInfo("PythiaDauVFilterMatchID") << "looking for " << particleID << endl;
0075 }
0076
0077 for (HepMC::GenEvent::particle_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end(); ++p) {
0078 if ((*p)->pdg_id() != particleID)
0079 continue;
0080
0081
0082 if (0 != motherID) {
0083 OK = 0;
0084 for (HepMC::GenVertex::particles_in_const_iterator des = (*p)->production_vertex()->particles_in_const_begin();
0085 des != (*p)->production_vertex()->particles_in_const_end();
0086 ++des) {
0087 if (fVerbose > 10) {
0088 edm::LogInfo("PythiaDauVFilterMatchID")
0089 << "mother: " << (*des)->pdg_id() << " pT: " << (*des)->momentum().perp()
0090 << " eta: " << (*des)->momentum().eta() << endl;
0091 }
0092 if (abs(motherID) == abs((*des)->pdg_id())) {
0093 OK = 1;
0094 break;
0095 }
0096 }
0097 }
0098 if (0 == OK)
0099 continue;
0100
0101
0102 std::vector<decayTarget> targets;
0103 for (unsigned int i = 0; i < dauIDs.size(); i++) {
0104 decayTarget target;
0105 target.pdgID = dauIDs[i];
0106 target.minPt = minptcut[i];
0107 target.maxPt = maxptcut;
0108 target.minEta = minetacut[i];
0109 target.maxEta = maxetacut[i];
0110 targets.push_back(target);
0111 }
0112 if (fVerbose > 10) {
0113 edm::LogInfo("PythiaDauVFilterMatchID") << "created target: ";
0114 for (unsigned int i = 0; i < targets.size(); i++) {
0115 edm::LogInfo("PythiaDauVFilterMatchID") << targets[i].pdgID << " ";
0116 }
0117 edm::LogInfo("PythiaDauVFilterMatchID") << endl;
0118 }
0119
0120
0121 int ndau = 0;
0122 if (fVerbose > 5) {
0123 edm::LogInfo("PythiaDauVFilterMatchID") << "found ID: " << (*p)->pdg_id() << " pT: " << (*p)->momentum().perp()
0124 << " eta: " << (*p)->momentum().eta() << endl;
0125 }
0126 vparticles.push_back((*p)->pdg_id());
0127 if ((*p)->end_vertex()) {
0128 for (HepMC::GenVertex::particle_iterator des = (*p)->end_vertex()->particles_begin(HepMC::children);
0129 des != (*p)->end_vertex()->particles_end(HepMC::children);
0130 ++des) {
0131 if (TMath::Abs((*des)->pdg_id()) == 22) {
0132 continue;
0133 }
0134 ++ndau;
0135 if (fVerbose > 5) {
0136 edm::LogInfo("PythiaDauVFilterMatchID") << "ID: " << (*des)->pdg_id() << " pT: " << (*des)->momentum().perp()
0137 << " eta: " << (*des)->momentum().eta() << endl;
0138 }
0139 {
0140 int matchedIdx = -1;
0141 for (unsigned int i = 0; i < targets.size(); i++) {
0142 if ((*des)->pdg_id() != targets[i].pdgID) {
0143 continue;
0144 }
0145 if (fVerbose > 5) {
0146 edm::LogInfo("PythiaDauVFilterMatchID") << "i = " << i << " pT = " << (*des)->momentum().perp()
0147 << " eta = " << (*des)->momentum().eta() << endl;
0148 }
0149
0150 if ((*des)->momentum().perp() > targets[i].minPt && (*des)->momentum().perp() < targets[i].maxPt &&
0151 (*des)->momentum().eta() > targets[i].minEta && (*des)->momentum().eta() < targets[i].maxEta) {
0152 vparticles.push_back((*des)->pdg_id());
0153 if (fVerbose > 2) {
0154 edm::LogInfo("PythiaDauVFilterMatchID")
0155 << " accepted this particle " << (*des)->pdg_id() << " pT = " << (*des)->momentum().perp()
0156 << " eta = " << (*des)->momentum().eta() << endl;
0157 }
0158 matchedIdx = i;
0159 break;
0160 }
0161 }
0162 if (matchedIdx > -1) {
0163 targets.erase(targets.begin() + matchedIdx);
0164 }
0165 if (fVerbose > 10) {
0166 edm::LogInfo("PythiaDauVFilterMatchID") << "remaining targets: ";
0167 for (unsigned int i = 0; i < targets.size(); i++) {
0168 edm::LogInfo("PythiaDauVFilterMatchID") << targets[i].pdgID << " ";
0169 }
0170 edm::LogInfo("PythiaDauVFilterMatchID") << endl;
0171 }
0172 }
0173 }
0174 }
0175
0176 if (ndau == ndaughters && targets.empty()) {
0177 accepted = true;
0178 if (fVerbose > 0) {
0179 edm::LogInfo("PythiaDauVFilterMatchID") << " accepted this decay: ";
0180 for (unsigned int iv = 0; iv < vparticles.size(); ++iv)
0181 edm::LogInfo("PythiaDauVFilterMatchID") << vparticles[iv] << " ";
0182 edm::LogInfo("PythiaDauVFilterMatchID") << " from mother = " << motherID << endl;
0183 }
0184 break;
0185 }
0186 }
0187
0188 int anti_particleID = -particleID;
0189 if (!(fLookupGen->particleData.isParticle(anti_particleID))) {
0190 anti_particleID = 0;
0191 if (fVerbose > 5)
0192 edm::LogInfo("PythiaDauVFilterMatchID")
0193 << "Particle " << particleID << " is its own anti-particle, skipping further testing " << endl;
0194 }
0195 if (!accepted && chargeconju && anti_particleID) {
0196 OK = 1;
0197
0198 for (HepMC::GenEvent::particle_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end(); ++p) {
0199 if ((*p)->pdg_id() != anti_particleID)
0200 continue;
0201
0202
0203 if (0 != motherID) {
0204 OK = 0;
0205 for (HepMC::GenVertex::particles_in_const_iterator des = (*p)->production_vertex()->particles_in_const_begin();
0206 des != (*p)->production_vertex()->particles_in_const_end();
0207 ++des) {
0208 if (fVerbose > 10) {
0209 edm::LogInfo("PythiaDauVFilterMatchID")
0210 << "mother: " << (*des)->pdg_id() << " pT: " << (*des)->momentum().perp()
0211 << " eta: " << (*des)->momentum().eta() << endl;
0212 }
0213 if (abs(motherID) == abs((*des)->pdg_id())) {
0214 OK = 1;
0215 break;
0216 }
0217 }
0218 }
0219 if (0 == OK)
0220 continue;
0221
0222
0223 std::vector<decayTarget> targets;
0224 for (unsigned int i = 0; i < dauIDs.size(); i++) {
0225 decayTarget target;
0226 int IDanti = -dauIDs[i];
0227 if (!(fLookupGen->particleData.isParticle(IDanti)))
0228 IDanti = dauIDs[i];
0229 target.pdgID = IDanti;
0230 target.minPt = minptcut[i];
0231 target.maxPt = maxptcut;
0232 target.minEta = minetacut[i];
0233 target.maxEta = maxetacut[i];
0234 targets.push_back(target);
0235 }
0236 if (fVerbose > 10) {
0237 edm::LogInfo("PythiaDauVFilterMatchID") << "created anti target: ";
0238 for (unsigned int i = 0; i < targets.size(); i++) {
0239 edm::LogInfo("PythiaDauVFilterMatchID") << targets[i].pdgID << " ";
0240 }
0241 edm::LogInfo("PythiaDauVFilterMatchID") << endl;
0242 }
0243
0244 if (fVerbose > 5) {
0245 edm::LogInfo("PythiaDauVFilterMatchID") << "found ID: " << (*p)->pdg_id() << " pT: " << (*p)->momentum().perp()
0246 << " eta: " << (*p)->momentum().eta() << endl;
0247 }
0248 vparticles.push_back((*p)->pdg_id());
0249 int ndau = 0;
0250 if ((*p)->end_vertex()) {
0251 for (HepMC::GenVertex::particle_iterator des = (*p)->end_vertex()->particles_begin(HepMC::children);
0252 des != (*p)->end_vertex()->particles_end(HepMC::children);
0253 ++des) {
0254 if (TMath::Abs((*des)->pdg_id()) == 22) {
0255 continue;
0256 }
0257 ++ndau;
0258 if (fVerbose > 5) {
0259 edm::LogInfo("PythiaDauVFilterMatchID")
0260 << "ID: " << (*des)->pdg_id() << " pT: " << (*des)->momentum().perp()
0261 << " eta: " << (*des)->momentum().eta() << endl;
0262 }
0263 {
0264 int matchedIdx = -1;
0265 for (unsigned int i = 0; i < targets.size(); i++) {
0266 if ((*des)->pdg_id() != targets[i].pdgID) {
0267 continue;
0268 }
0269 if (fVerbose > 5) {
0270 edm::LogInfo("PythiaDauVFilterMatchID") << "i = " << i << " pT = " << (*des)->momentum().perp()
0271 << " eta = " << (*des)->momentum().eta() << endl;
0272 }
0273
0274 if ((*des)->momentum().perp() > targets[i].minPt && (*des)->momentum().perp() < targets[i].maxPt &&
0275 (*des)->momentum().eta() > targets[i].minEta && (*des)->momentum().eta() < targets[i].maxEta) {
0276 vparticles.push_back((*des)->pdg_id());
0277 if (fVerbose > 2) {
0278 edm::LogInfo("PythiaDauVFilterMatchID")
0279 << " accepted this particle " << (*des)->pdg_id() << " pT = " << (*des)->momentum().perp()
0280 << " eta = " << (*des)->momentum().eta() << endl;
0281 }
0282 matchedIdx = i;
0283 break;
0284 }
0285 }
0286 if (matchedIdx > -1) {
0287 targets.erase(targets.begin() + matchedIdx);
0288 }
0289 if (fVerbose > 10) {
0290 edm::LogInfo("PythiaDauVFilterMatchID") << "remaining targets: ";
0291 for (unsigned int i = 0; i < targets.size(); i++) {
0292 edm::LogInfo("PythiaDauVFilterMatchID") << targets[i].pdgID << " ";
0293 }
0294 edm::LogInfo("PythiaDauVFilterMatchID") << endl;
0295 }
0296 }
0297 }
0298 }
0299 if (ndau == ndaughters && targets.empty()) {
0300 accepted = true;
0301 if (fVerbose > 0) {
0302 edm::LogInfo("PythiaDauVFilterMatchID") << " accepted this decay: ";
0303 for (unsigned int iv = 0; iv < vparticles.size(); ++iv)
0304 edm::LogInfo("PythiaDauVFilterMatchID") << vparticles[iv] << " ";
0305 edm::LogInfo("PythiaDauVFilterMatchID") << " from mother = " << motherID << endl;
0306 }
0307 break;
0308 }
0309 }
0310 }
0311
0312 delete myGenEvent;
0313 return accepted;
0314 }