File indexing completed on 2024-04-06 12:13:42
0001
0002 #include "GeneratorInterface/GenFilters/plugins/PythiaProbeFilter.h"
0003
0004 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0005 #include <iostream>
0006 #include <memory>
0007
0008 using namespace edm;
0009 using namespace std;
0010 using namespace Pythia8;
0011
0012 PythiaProbeFilter::PythiaProbeFilter(const edm::ParameterSet& iConfig)
0013 : token_(consumes<edm::HepMCProduct>(
0014 edm::InputTag(iConfig.getUntrackedParameter("moduleLabel", std::string("generator")), "unsmeared"))),
0015 particleID(iConfig.getUntrackedParameter("ParticleID", 0)),
0016 MomID(iConfig.getUntrackedParameter("MomID", 0)),
0017 GrandMomID(iConfig.getUntrackedParameter("GrandMomID", 0)),
0018 chargeconju(iConfig.getUntrackedParameter("ChargeConjugation", true)),
0019 nsisters(iConfig.getUntrackedParameter("NumberOfSisters", 0)),
0020 naunts(iConfig.getUntrackedParameter("NumberOfAunts", 0)),
0021 minptcut(iConfig.getUntrackedParameter("MinPt", 0.)),
0022 maxptcut(iConfig.getUntrackedParameter("MaxPt", 14000.)),
0023 minetacut(iConfig.getUntrackedParameter("MinEta", -10.)),
0024 maxetacut(iConfig.getUntrackedParameter("MaxEta", 10.)),
0025 countQEDCorPhotons(iConfig.getUntrackedParameter("countQEDCorPhotons", false)) {
0026
0027 vector<int> defID;
0028 defID.push_back(0);
0029 exclsisIDs = iConfig.getUntrackedParameter<vector<int> >("SisterIDs", defID);
0030 exclauntIDs = iConfig.getUntrackedParameter<vector<int> >("AuntIDs", defID);
0031 identicalParticle = false;
0032 for (unsigned int ilist = 0; ilist < exclsisIDs.size(); ++ilist) {
0033 if (fabs(exclsisIDs[ilist]) == fabs(particleID))
0034 identicalParticle = true;
0035 }
0036
0037 edm::LogInfo("PythiaProbeFilter::PythiaProbeFilter") << "Creating pythia8 instance for particle properties" << endl;
0038 if (!fLookupGen.get())
0039 fLookupGen = std::make_unique<Pythia>();
0040 }
0041
0042 PythiaProbeFilter::~PythiaProbeFilter() {
0043
0044
0045 }
0046
0047 bool PythiaProbeFilter::AlreadyExcludedCheck(std::vector<unsigned int> excludedList, unsigned int current_part) const {
0048 bool result = false;
0049 for (unsigned int checkNow : excludedList) {
0050 if (current_part != checkNow)
0051 continue;
0052 result = true;
0053 break;
0054 }
0055 return result;
0056 }
0057
0058
0059
0060
0061 bool PythiaProbeFilter::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 const HepMC::GenEvent* myGenEvent = evt->GetEvent();
0068
0069
0070 for (HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end();
0071 ++p) {
0072
0073 if (fabs((*p)->pdg_id()) == fabs(particleID))
0074 if ((*p)->pdg_id() != particleID && !chargeconju)
0075 continue;
0076
0077 if (fabs((*p)->pdg_id()) != fabs(particleID) && chargeconju)
0078 continue;
0079
0080 if ((*p)->momentum().perp() < minptcut || (*p)->momentum().perp() > maxptcut)
0081 continue;
0082 if ((*p)->momentum().eta() < minetacut || (*p)->momentum().eta() > maxetacut)
0083 continue;
0084
0085 bool excludeTagParticle = false;
0086 if (naunts == 0) {
0087 if ((*p)->production_vertex()) {
0088 for (HepMC::GenVertex::particle_iterator anc = (*p)->production_vertex()->particles_begin(HepMC::parents);
0089 anc != (*p)->production_vertex()->particles_end(HepMC::parents);
0090 ++anc) {
0091 if (fabs((*anc)->pdg_id()) != fabs(MomID) && chargeconju)
0092 continue;
0093 else if ((*anc)->pdg_id() != MomID && !chargeconju)
0094 continue;
0095 int nsis = 0;
0096 int exclsis = 0;
0097 std::vector<unsigned int> checklistSis;
0098 if ((*anc)->end_vertex()) {
0099 for (HepMC::GenVertex::particle_iterator sis = (*anc)->end_vertex()->particles_begin(HepMC::children);
0100 sis != (*anc)->end_vertex()->particles_end(HepMC::children);
0101 ++sis) {
0102
0103 if ((*p)->pdg_id() == (*sis)->pdg_id() && (identicalParticle || !chargeconju))
0104 continue;
0105 if (fabs((*p)->pdg_id()) == fabs((*sis)->pdg_id()) && !identicalParticle && chargeconju)
0106 continue;
0107
0108 if ((*sis)->pdg_id() == 22 && !countQEDCorPhotons)
0109 continue;
0110 nsis++;
0111
0112 for (unsigned int ilist = 0; ilist < exclsisIDs.size(); ++ilist) {
0113 if (AlreadyExcludedCheck(checklistSis, ilist)) {
0114 continue;
0115 }
0116 if (fabs(exclsisIDs[ilist]) == fabs((*sis)->pdg_id()) && chargeconju) {
0117 exclsis++;
0118 checklistSis.push_back(ilist);
0119 }
0120 if (exclsisIDs[ilist] == (*sis)->pdg_id() && !chargeconju) {
0121 exclsis++;
0122 checklistSis.push_back(ilist);
0123 }
0124 }
0125 }
0126 }
0127 if (nsis == exclsis && nsis == nsisters) {
0128 excludeTagParticle = true;
0129 break;
0130 }
0131 }
0132 }
0133 } else if (naunts > 0) {
0134
0135 if ((*p)->production_vertex()) {
0136 for (HepMC::GenVertex::particle_iterator anc = (*p)->production_vertex()->particles_begin(HepMC::parents);
0137 anc != (*p)->production_vertex()->particles_end(HepMC::parents);
0138 ++anc) {
0139 if (fabs((*anc)->pdg_id()) != fabs(MomID) && chargeconju)
0140 continue;
0141 else if ((*anc)->pdg_id() != MomID && !chargeconju)
0142 continue;
0143 int nsis = 0;
0144 int exclsis = 0;
0145 std::vector<unsigned int> checklistSis;
0146 int naunt = 0;
0147 int exclaunt = 0;
0148 std::vector<unsigned int> checklistAunt;
0149 if ((*anc)->end_vertex()) {
0150 for (HepMC::GenVertex::particle_iterator sis = (*anc)->end_vertex()->particles_begin(HepMC::children);
0151 sis != (*anc)->end_vertex()->particles_end(HepMC::children);
0152 ++sis) {
0153
0154 if ((*p)->pdg_id() == (*sis)->pdg_id() && (identicalParticle || !chargeconju))
0155 continue;
0156 if (fabs((*p)->pdg_id()) == fabs((*sis)->pdg_id()) && !identicalParticle && chargeconju)
0157 continue;
0158
0159 if ((*sis)->pdg_id() == 22 && !countQEDCorPhotons)
0160 continue;
0161 nsis++;
0162 for (unsigned int ilist = 0; ilist < exclsisIDs.size(); ++ilist) {
0163 if (AlreadyExcludedCheck(checklistSis, ilist))
0164 continue;
0165 if (fabs(exclsisIDs[ilist]) == fabs((*sis)->pdg_id()) && chargeconju) {
0166 exclsis++;
0167 checklistSis.push_back(ilist);
0168 }
0169 if (exclsisIDs[ilist] == (*sis)->pdg_id() && !chargeconju) {
0170 exclsis++;
0171 checklistSis.push_back(ilist);
0172 }
0173 }
0174 }
0175 }
0176
0177 if (nsis != exclsis || nsis != nsisters)
0178 break;
0179 if ((*anc)->production_vertex()) {
0180 for (HepMC::GenVertex::particle_iterator granc =
0181 (*anc)->production_vertex()->particles_begin(HepMC::parents);
0182 granc != (*anc)->production_vertex()->particles_end(HepMC::parents);
0183 ++granc) {
0184 if (fabs((*granc)->pdg_id()) != fabs(GrandMomID) && chargeconju)
0185 continue;
0186 else if ((*granc)->pdg_id() != GrandMomID && !chargeconju)
0187 continue;
0188 for (HepMC::GenVertex::particle_iterator aunt = (*granc)->end_vertex()->particles_begin(HepMC::children);
0189 aunt != (*granc)->end_vertex()->particles_end(HepMC::children);
0190 ++aunt) {
0191 if ((*aunt)->pdg_id() == (*anc)->pdg_id())
0192 continue;
0193 if ((*aunt)->pdg_id() == 22 && !countQEDCorPhotons)
0194 continue;
0195 naunt++;
0196 for (unsigned int ilist = 0; ilist < exclauntIDs.size(); ++ilist) {
0197 if (AlreadyExcludedCheck(checklistAunt, ilist))
0198 continue;
0199 if (fabs(exclauntIDs[ilist]) == fabs((*aunt)->pdg_id()) && chargeconju) {
0200 exclaunt++;
0201 checklistAunt.push_back(ilist);
0202 }
0203 if (exclauntIDs[ilist] == (*aunt)->pdg_id() && !chargeconju) {
0204 exclaunt++;
0205 checklistAunt.push_back(ilist);
0206 }
0207 }
0208 }
0209 }
0210 }
0211
0212 if (naunt == exclaunt && naunt == naunts) {
0213 excludeTagParticle = true;
0214 break;
0215 }
0216 }
0217 }
0218 }
0219 if (excludeTagParticle)
0220 continue;
0221 accepted = true;
0222 break;
0223 }
0224
0225 if (accepted) {
0226 return true;
0227 } else {
0228 return false;
0229 }
0230 }