From 5d91b5b59d43d535e336e82e565fb76a08fe4515 Mon Sep 17 00:00:00 2001 From: Ka Wa Ho Date: Fri, 15 May 2026 16:46:59 -0500 Subject: [PATCH 1/8] minimal changes to make it compile --- .../interface/myEvtRandomEngine.h | 5 + .../plugins/EvtGen/EvtGenInterface.cc | 16 +- .../plugins/EvtGen/EvtGenInterface.h | 3 +- .../plugins/myEvtRandomEngine.cc | 4 + .../interface/Pythia8EvtGenPatched.h | 624 ++++++++++++++++++ .../Pythia8Interface/plugins/HMC3Py8EGun.cc | 2 +- .../plugins/Pythia8Hadronizer.cc | 2 +- .../plugins/Pythia8HepMC3Hadronizer.cc | 2 +- .../Pythia8Interface/src/Py8GunBase.cc | 2 +- 9 files changed, 642 insertions(+), 18 deletions(-) create mode 100644 GeneratorInterface/Pythia8Interface/interface/Pythia8EvtGenPatched.h diff --git a/GeneratorInterface/EvtGenInterface/interface/myEvtRandomEngine.h b/GeneratorInterface/EvtGenInterface/interface/myEvtRandomEngine.h index be9ab5345067f..abdf86ff981a1 100644 --- a/GeneratorInterface/EvtGenInterface/interface/myEvtRandomEngine.h +++ b/GeneratorInterface/EvtGenInterface/interface/myEvtRandomEngine.h @@ -32,6 +32,10 @@ class myEvtRandomEngine : public EvtRandomEngine { double random() override; + void setSeed(unsigned long int seed) override; + + unsigned long int lastSeed() const override; + void setRandomEngine(CLHEP::HepRandomEngine* v) { the_engine = v; } CLHEP::HepRandomEngine* engine() const { return the_engine; } @@ -40,5 +44,6 @@ class myEvtRandomEngine : public EvtRandomEngine { void throwNullPtr() const; CLHEP::HepRandomEngine* the_engine; + unsigned long int m_lastSeed = 0; }; #endif diff --git a/GeneratorInterface/EvtGenInterface/plugins/EvtGen/EvtGenInterface.cc b/GeneratorInterface/EvtGenInterface/plugins/EvtGen/EvtGenInterface.cc index c71fe59e637a5..d7b6deaa82d7b 100644 --- a/GeneratorInterface/EvtGenInterface/plugins/EvtGen/EvtGenInterface.cc +++ b/GeneratorInterface/EvtGenInterface/plugins/EvtGen/EvtGenInterface.cc @@ -36,6 +36,7 @@ #include "EvtGenBase/EvtDecayBase.hh" #include "EvtGenBase/EvtId.hh" #include "EvtGenBase/EvtDecayTable.hh" +#include "EvtGenBase/EvtPDL.hh" #include "EvtGenBase/EvtParticle.hh" #include "EvtGenBase/EvtParticleFactory.hh" #include "EvtGenBase/EvtHepMCEvent.hh" @@ -50,8 +51,6 @@ using namespace gen; using namespace edm; -CLHEP::HepRandomEngine* EvtGenInterface::fRandomEngine; - EvtGenInterface::EvtGenInterface(const ParameterSet& pset) { fPSet = new ParameterSet(pset); the_engine = new myEvtRandomEngine(nullptr); @@ -496,8 +495,8 @@ HepMC::GenEvent* EvtGenInterface::decay(HepMC::GenEvent* evt) { } EvtId idEvt = EvtPDL::evtIdFromStdHep(idHep); int ipart = idEvt.getId(); - EvtDecayTable* evtDecayTable = EvtDecayTable::getInstance(); - if (!isforced && isDefaultEvtGen && ipart != -1 && evtDecayTable->getNMode(ipart) != 0) { + EvtDecayTable& evtDecayTable = EvtDecayTable::getInstance(); + if (!isforced && isDefaultEvtGen && ipart != -1 && evtDecayTable.getNMode(ipart) != 0) { addToHepMC(*p, idEvt, evt, true); // generate decay and remove daugther if they are forced } } @@ -616,17 +615,10 @@ void EvtGenInterface::update_particles(HepMC::GenParticle* partHep, HepMC::GenEv void EvtGenInterface::setRandomEngine(CLHEP::HepRandomEngine* v) { the_engine->setRandomEngine(v); - fRandomEngine = v; } double EvtGenInterface::flat() { - if (!fRandomEngine) { - throw cms::Exception("LogicError") - << "EvtGenInterface::flat: Attempt to generate random number when engine pointer is null\n" - << "This might mean that the code was modified to generate a random number outside the\n" - << "event and beginLuminosityBlock methods, which is not allowed.\n"; - } - return fRandomEngine->flat(); + return the_engine->random(); } bool EvtGenInterface::findLastinChain(HepMC::GenParticle*& p) { diff --git a/GeneratorInterface/EvtGenInterface/plugins/EvtGen/EvtGenInterface.h b/GeneratorInterface/EvtGenInterface/plugins/EvtGen/EvtGenInterface.h index 0dad34911846c..f3ae6fa5153bc 100644 --- a/GeneratorInterface/EvtGenInterface/plugins/EvtGen/EvtGenInterface.h +++ b/GeneratorInterface/EvtGenInterface/plugins/EvtGen/EvtGenInterface.h @@ -43,7 +43,7 @@ namespace gen { const std::vector& operatesOnParticles() override { return m_PDGs; } HepMC::GenEvent* decay(HepMC::GenEvent*) override; void setRandomEngine(CLHEP::HepRandomEngine* v) override; - static double flat(); + double flat(); private: bool addToHepMC(HepMC::GenParticle* partHep, const EvtId& idEvt, HepMC::GenEvent* theEvent, bool del_daug); @@ -67,7 +67,6 @@ namespace gen { int BmixingOption = 1; edm::ParameterSet* fPSet; - static CLHEP::HepRandomEngine* fRandomEngine; myEvtRandomEngine* the_engine; }; } // namespace gen diff --git a/GeneratorInterface/EvtGenInterface/plugins/myEvtRandomEngine.cc b/GeneratorInterface/EvtGenInterface/plugins/myEvtRandomEngine.cc index 43506f9a9b132..5d39dd28369a3 100644 --- a/GeneratorInterface/EvtGenInterface/plugins/myEvtRandomEngine.cc +++ b/GeneratorInterface/EvtGenInterface/plugins/myEvtRandomEngine.cc @@ -30,6 +30,10 @@ double myEvtRandomEngine::random() { return the_engine->flat(); } +void myEvtRandomEngine::setSeed(unsigned long int seed) { m_lastSeed = seed; } + +unsigned long int myEvtRandomEngine::lastSeed() const { return m_lastSeed; } + void myEvtRandomEngine::throwNullPtr() const { throw edm::Exception(edm::errors::LogicError) << "The EvtGen code attempted to a generate random number while\n" << "the engine pointer was null. This might mean that the code\n" diff --git a/GeneratorInterface/Pythia8Interface/interface/Pythia8EvtGenPatched.h b/GeneratorInterface/Pythia8Interface/interface/Pythia8EvtGenPatched.h new file mode 100644 index 0000000000000..c0cec08a0fe9a --- /dev/null +++ b/GeneratorInterface/Pythia8Interface/interface/Pythia8EvtGenPatched.h @@ -0,0 +1,624 @@ +// Local patched copy of Pythia8Plugins/EvtGen.h adapted to EvtGen 03.x API. +// Diverges from the upstream pythia8 header only in: +// - EvtGenRandom now overrides setSeed(unsigned long)/lastSeed() pure virtuals +// - EvtDecayTable::getInstance() returns EvtDecayTable& (not pointer) + +#ifndef CMSSW_Pythia8_EvtGen_Patched_H +#define CMSSW_Pythia8_EvtGen_Patched_H + +#include "Pythia8/Pythia.h" +#include "EvtGen/EvtGen.hh" +#include "EvtGenBase/EvtRandomEngine.hh" +#include "EvtGenBase/EvtParticle.hh" +#include "EvtGenBase/EvtParticleFactory.hh" +#include "EvtGenBase/EvtPDL.hh" +#include "EvtGenBase/EvtDecayTable.hh" +#include "EvtGenBase/EvtParticleDecayList.hh" +#include "EvtGenBase/EvtDecayBase.hh" +#include "EvtGenExternal/EvtExternalGenList.hh" + +namespace Pythia8 { + +//========================================================================== + +// A class to wrap the Pythia random number generator for use by EvtGen. + +class EvtGenRandom : public EvtRandomEngine { + +public: + + // Constructor. + EvtGenRandom(Rndm *rndmPtrIn) {rndmPtr = rndmPtrIn;} + + // Return a random number. + double random() override {if (rndmPtr) return rndmPtr->flat(); else return -1.0;} + + // EvtGen 03 pure virtuals. Record-only: Pythia owns the actual RNG state. + void setSeed(unsigned long int seed) override {lastSeedVal = seed;} + unsigned long int lastSeed() const override {return lastSeedVal;} + + // The random number pointer. + Rndm *rndmPtr; + +private: + unsigned long int lastSeedVal{0}; + +}; + +//========================================================================== + +// A class to perform decays via the external EvtGen decay program, +// see http://evtgen.warwick.ac.uk/, the program manual provided with +// the EvtGen distribution, and D. J. Lange, +// Nucl. Instrum. Meth. A462, 152 (2001) for details. + +// EvtGen performs a series of decays from some initial particle +// decay, rather than just a single decay, and so EvtGen cannot be +// interfaced through the standard external DecayHandler class without +// considerable complication. Consequently, EvtGen is called on the +// complete event record after all steps of Pythia are completed. + +// Oftentimes a specific "signal" decay is needed to occur once in an +// event, and all other decays performed normally. This is possible +// via reading in a user decay file (with readDecayFile) and creating +// aliased particles with names ending with signalSuffix. By default, +// this is "_SIGNAL". When decay() is called, all particles in the +// Pythia event record that are of the same types as the signal +// particles are collected. One is selected at random and decayed via +// the channel(s) defined for that aliased signal particle. All other +// particles are decayed normally. The weight for the event is +// calculated and returned. + +// It is also possible to specify a status needed to consider a +// particle as a signal candidate. This can be done by modifying the +// signals map, e.g. if the tau- is a signal candidate, then +// EvtGenDecays.signals[15].status = 201 +// will only only select as candidates any tau- with this status. This +// allows the event record to be changed before decays, so only +// certain particles are selected as possible signal candidates +// (e.g. passing kinematic requirements). + +// Please note that particles produced from a signal candidate decay +// are not searched for additional signal candidates. This means that +// if B0 and tau- have been designated as signal, then a tau- from a +// W- decay would be a signal candidate, while a tau- from a B0 decay +// would not. This restriction arises from the additional complexity +// of allowing recursive signal decays. The following statuses are +// used: 93 for particles decayed with EvtGen, 94 for particles +// oscillated with EvtGen, 95 for signal particles, and 96 for signal +// particles from an oscillation. + +class EvtGenDecays { + +public: + + // Constructor. + EvtGenDecays(Pythia *pythiaPtrIn, string decayFile, string particleDataFile, + EvtExternalGenList *extPtrIn = 0, EvtAbsRadCorr *fsrPtrIn = 0, + int mixing = 1, bool xml = false, bool limit = true, + bool extUse = true, bool fsrUse = true); + + // Destructor. + ~EvtGenDecays() { + if (evtgen) delete evtgen; + if (extOwner && extPtr) delete extPtr; + if (fsrOwner && fsrPtr) delete fsrPtr; + } + + // Perform all decays and return the event weight. + double decay(); + + // Stop EvtGen decaying a particle. + void exclude(int id) {excIds.insert(id);} + + // Update the Pythia particle database from EvtGen. + void updatePythia(); + + // Update the EvtGen particle database from Pythia. + void updateEvtGen(); + + // Read an EvtGen user decay file. + void readDecayFile(string decayFile, bool xml = false) { + evtgen->readUDecay(decayFile.c_str(), xml); updateData();} + + // External model pointer and FSR model pointer. + bool extOwner, fsrOwner; + EvtExternalGenList *extPtr; + EvtAbsRadCorr *fsrPtr; + std::list models; + + // Map of signal particle info. + struct Signal {int status; EvtId egId; vector bfs; vector map; + EvtParticleDecayList modes;}; + map signals; + + // The suffix indicating an EvtGen particle or alias is signal. + string signalSuffix; + +protected: + + // Update the particles to decay with EvtGen, and the selected signals. + void updateData(bool final = false); + + // Update the Pythia event record with an EvtGen decay tree. + void updateEvent(Particle *pyPro, EvtParticle *egPro, + vector *pySigs = 0, vector *egSigs = 0, + vector *bfs = 0, double *wgt = 0); + + // Check if a particle can decay. + bool checkVertex(Particle *pyPro); + + // Check if a particle is signal. + bool checkSignal(Particle *pyPro); + + // Check if an EvtGen particle has oscillated. + bool checkOsc(EvtParticle *egPro); + + // Number of times to try a decay sampling (constant). + static const int NTRYDECAY = 1000; + + // The pointer to the associated Pythia object. + Pythia *pythiaPtr; + + // Random number wrapper for EvtGen. + EvtGenRandom rndm; + + // The EvtGen object. + EvtGen *evtgen; + + // Set of particle IDs to include and exclude decays with EvtGen. + set incIds, excIds; + + // Flag whether the final particle update has been performed. + bool updated; + + // The selected signal iterator. + map::iterator signal; + + // Parameters used to check if a particle should decay (as set via Pythia). + double tau0Max, tauMax, rMax, xyMax, zMax; + bool limitTau0, limitTau, limitRadius, limitCylinder, limitDecay; + +}; + +//-------------------------------------------------------------------------- + +// The constructor. + +// The EvtGenDecays object is associated with a single Pythia +// instance. This is to ensure a consistent random number generator +// across the two, as well as any updates to particle data, etc. Note +// that if multiple EvtGenDecays objects exist, that they will modify +// one anothers particle databases due to the design of EvtGen. + +// This constructor also sets all particles to be decayed by EvtGen as +// stable within Pythia. The parameters within Pythia used to check if +// a particle should be decayed, as described in the "Particle Decays" +// section of the Pythia manual, are set. Note that if the variable +// "limit" is set to "false", then no check will be made before +// decaying a particle with EvtGen. + +// The constructor is designed to have the exact same form as the +// EvtGen constructor except for these five differences. +// (1) The first variable is the pointer to the Pythia object. +// (2) The third last argument is a flag to limit decays based on the +// Pythia criteria (based on the particle decay vertex). +// (3) The second last argument is a flag if external models should be +// passed to EvtGen (default is true). +// (4) The last argument is a flag if an FSR model should be passed +// to EvtGen (default is true). +// (5) No random engine pointer is passed, as this is obtained from +// Pythia. + +// pythiaPtrIn: the pointer to the associated Pythia object. +// decayFile: the name of the decay file to pass to EvtGen. +// particleDataFile: the name of the particle data file to pass to EvtGen. +// extPtrIn: the optional EvtExternalGenList pointer, this must be +// be provided if fsrPtrIn is provided to avoid double +// initializations. +// fsrPtrIn: the EvtAbsRadCorr pointer to pass to EvtGen. +// mixing: the mixing type to pass to EvtGen. +// xml: flag to use XML files to pass to EvtGen. +// limit: flag to limit particle decays based on Pythia criteria. +// extUse: flag to use external models with EvtGen. +// fsrUse: flag to use radiative correction engine with EvtGen. + +EvtGenDecays::EvtGenDecays(Pythia *pythiaPtrIn, string decayFile, + string particleDataFile, EvtExternalGenList *extPtrIn, + EvtAbsRadCorr *fsrPtrIn, int mixing, bool xml, bool limit, + bool extUse, bool fsrUse) : extPtr(extPtrIn), fsrPtr(fsrPtrIn), + signalSuffix("_SIGNAL"), pythiaPtr(pythiaPtrIn), rndm(&pythiaPtr->rndm), + updated(false) { + + // Initialize EvtGen. + if (!extPtr && fsrPtr) { + cout << "Error in EvtGenDecays::" + << "EvtGenDecays: extPtr is null but fsrPtr is provided\n"; + return; + } + if (extPtr) extOwner = false; + else {extOwner = true; extPtr = new EvtExternalGenList();} + if (fsrPtr) fsrOwner = false; + else {fsrOwner = true; fsrPtr = extPtr->getPhotosModel();} + models = extPtr->getListOfModels(); + evtgen = new EvtGen(decayFile.c_str(), particleDataFile.c_str(), + &rndm, fsrUse ? fsrPtr : 0, extUse ? &models : 0, mixing, xml); + + // Get the Pythia decay limits. + if (!pythiaPtr) return; + limitTau0 = pythiaPtr->settings.flag("ParticleDecays:limitTau0"); + tau0Max = pythiaPtr->settings.parm("ParticleDecays:tau0Max"); + limitTau = pythiaPtr->settings.flag("ParticleDecays:limitTau"); + tauMax = pythiaPtr->settings.parm("ParticleDecays:tauMax"); + limitRadius = pythiaPtr->settings.flag("ParticleDecays:limitRadius"); + rMax = pythiaPtr->settings.parm("ParticleDecays:rMax"); + limitCylinder = pythiaPtr->settings.flag("ParticleDecays:limitCylinder"); + xyMax = pythiaPtr->settings.parm("ParticleDecays:xyMax"); + zMax = pythiaPtr->settings.parm("ParticleDecays:zMax"); + limitDecay = limit && (limitTau0 || limitTau || + limitRadius || limitCylinder); + +} + +//-------------------------------------------------------------------------- + +// Perform all decays and return the event weight. + +// All particles in the event record that can be decayed by EvtGen are +// decayed. If a particle is a signal particle, then this is stored in +// a vector of signal particles. A signal particle is only stored if +// its status is the same as the status provided in the signals map. A +// negative status in the signal map indicates that all statuses +// should be accepted. After all signal particles are identified, one +// is randomly chosen and decayed as signal. The remainder are decayed +// normally. + +// Forcing a signal decay changes the weight of an event from unity, +// and so the relative event weight is returned, given the forced +// signal decay. A weight of 0 indicates no signal in the event, while +// a weight of -1 indicates something is wrong, e.g. either the Pythia +// or EvtGen pointers are not available or the number of tries has +// been exceeded. For the event weight to be valid, one should not +// change the absolute branching fractions in the signal and inclusive +// definitions, but rather just remove the unwanted decay channels +// from the signal decay definition. + +double EvtGenDecays::decay() { + + // Reset the signal and signal counters. + if (!pythiaPtr || !evtgen) return -1; + if (!updated) updateData(true); + + // Loop over all particles in the Pythia event. + Event &event = pythiaPtr->event; + vector pySigs; vector egSigs, egPrts; + vector bfs; double wgt(1.); + for (int iPro = 0; iPro < event.size(); ++iPro) { + + // Check particle is final and can be decayed by EvtGen. + Particle *pyPro = &event[iPro]; + if (!pyPro->isFinal()) continue; + if (incIds.find(pyPro->id()) == incIds.end()) continue; + if (pyPro->status() == 93 || pyPro->status() == 94) continue; + + // Decay the progenitor with EvtGen. + EvtParticle *egPro = EvtParticleFactory::particleFactory + (EvtPDL::evtIdFromStdHep(pyPro->id()), + EvtVector4R(pyPro->e(), pyPro->px(), pyPro->py(), pyPro->pz())); + egPrts.push_back(egPro); + egPro->setDiagonalSpinDensity(); + evtgen->generateDecay(egPro); + pyPro->tau(egPro->getLifetime()); + if (!checkVertex(pyPro)) continue; + + // Add oscillations to event record. + if (checkOsc(egPro)) + updateEvent(pyPro, egPro, &pySigs, &egSigs, &bfs, &wgt); + + // Undo decay if signal (duplicate to stop oscillations). + else if (checkSignal(pyPro)) { + pySigs.push_back(pyPro->index()); + egSigs.push_back(egPro); + bfs.push_back(signal->second.bfs[0]); + wgt *= 1 - bfs.back(); + egPro->deleteDaughters(); + EvtParticle *egDau = EvtParticleFactory::particleFactory + (EvtPDL::evtIdFromStdHep(pyPro->id()), + EvtVector4R(pyPro->e(), pyPro->px(), pyPro->py(), pyPro->pz())); + egDau->addDaug(egPro); + egDau->setDiagonalSpinDensity(); + + // If not signal, add to event record. + } else updateEvent(pyPro, egPro, &pySigs, &egSigs, &bfs, &wgt); + } + if (pySigs.size() == 0) { + for (int iPrt = 0; iPrt < (int)egPrts.size(); ++iPrt) + egPrts[iPrt]->deleteTree(); + return 0; + } + + // Determine the decays of the signal particles (signal or background). + vector modes; int force(-1), n(0); + for (int iTry = 1; iTry <= NTRYDECAY; ++iTry) { + modes.clear(); force = pythiaPtr->rndm.pick(bfs); n = 0; + for (int iSig = 0; iSig < (int)pySigs.size(); ++iSig) { + if (iSig == force) modes.push_back(0); + else modes.push_back(pythiaPtr->rndm.flat() > bfs[iSig]); + if (modes.back() == 0) ++n; + } + if (pythiaPtr->rndm.flat() <= 1.0 / n) break; + if (iTry == NTRYDECAY) { + wgt = 2.; + cout << "Warning in EvtGenDecays::decay: maximum " + << "number of signal decay attempts exceeded"; + } + } + + // Decay the signal particles and mark forced decay. + for (int iSig = 0; iSig < (int)pySigs.size(); ++iSig) { + EvtParticle *egSig = egSigs[iSig]; + Particle *pySig = &event[pySigs[iSig]]; + signal = signals.find(pySig->id()); + if (egSig->getNDaug() > 0) egSig = egSig->getDaug(0); + if (modes[iSig] == 0) egSig->setId(signal->second.egId); + else { + signal->second.modes.getDecayModel(egSig); + egSig->setChannel(signal->second.map[egSig->getChannel()]); + } + if (iSig == force) + pySig->status(pySig->status() == 92 || pySig->status() == 94 ? 96 : 95); + evtgen->generateDecay(egSig); + updateEvent(pySig, egSigs[iSig]); + } + + // Delete all EvtGen particles and return weight. + for (int iPrt = 0; iPrt < (int)egPrts.size(); ++iPrt) + egPrts[iPrt]->deleteTree(); + return 1. - wgt; + +} + +//-------------------------------------------------------------------------- + +// Update the Pythia particle database from EvtGen. + +// Note that only the particle spin type, charge type, nominal mass, +// width, minimum mass, maximum mass, and nominal lifetime are +// set. The name string is not set. + +void EvtGenDecays::updatePythia() { + if (!pythiaPtr || !evtgen) return; + for (int entry = 0; entry < (int)EvtPDL::entries(); ++entry) { + EvtId egId = EvtPDL::getEntry(entry); + int pyId = EvtPDL::getStdHep(egId); + pythiaPtr->particleData.spinType (pyId, EvtPDL::getSpinType(egId)); + pythiaPtr->particleData.chargeType(pyId, EvtPDL::chg3(egId)); + pythiaPtr->particleData.m0 (pyId, EvtPDL::getMass(egId)); + pythiaPtr->particleData.mWidth (pyId, EvtPDL::getWidth(egId)); + pythiaPtr->particleData.mMin (pyId, EvtPDL::getMinMass(egId)); + pythiaPtr->particleData.mMax (pyId, EvtPDL::getMaxMass(egId)); + pythiaPtr->particleData.tau0 (pyId, EvtPDL::getctau(egId)); + } +} + +//-------------------------------------------------------------------------- + +// Update the EvtGen particle database from Pythia. + +// The particle update database between Pythia and EvtGen is not +// symmetric. Particularly, it is not possible to update the spin +// type, charge, or nominal lifetime in the EvtGen particle database. + +void EvtGenDecays::updateEvtGen() { + if (!pythiaPtr || !evtgen) return; + int pyId = pythiaPtr->particleData.nextId(1); + while (pyId != 0) { + EvtId egId = EvtPDL::evtIdFromStdHep(pyId); + EvtPDL::reSetMass (egId, pythiaPtr->particleData.m0(pyId)); + EvtPDL::reSetWidth (egId, pythiaPtr->particleData.mWidth(pyId)); + EvtPDL::reSetMassMin(egId, pythiaPtr->particleData.mMin(pyId)); + EvtPDL::reSetMassMax(egId, pythiaPtr->particleData.mMax(pyId)); + pyId = pythiaPtr->particleData.nextId(pyId); + } +} + +//-------------------------------------------------------------------------- + +// Update the particles to decay with EvtGen, and the selected signals. + +// If final is false, then only signals are initialized in the signal +// map. Any particle or alias that ends with signalSuffix is taken as +// a signal particle. If final is true all particle entries in EvtGen +// are checked to see if they should be set stable in Pythia. If an +// EvtGen particle has no decay modes, then Pythia is still allowed to +// decay the particle. Additionally, the signal decay channels are +// turned off for the non-aliased signal particle. + +void EvtGenDecays::updateData(bool final) { + + // Loop over the EvtGen entries. + if (!pythiaPtr) return; + EvtDecayTable &egTable = EvtDecayTable::getInstance(); + for (int iEntry = 0; iEntry < (int)EvtPDL::entries(); ++iEntry) { + EvtId egId = EvtPDL::getEntry(iEntry); + int pyId = EvtPDL::getStdHep(egId); + if (egTable.getNModes(egId) == 0) continue; + if (excIds.find(pyId) != excIds.end()) continue; + + // Stop Pythia from decaying the particle and include in decay set. + if (final) { + incIds.insert(pyId); + pythiaPtr->particleData.mayDecay(pyId, false); + } + + // Check for signal. + string egName = EvtPDL::name(egId); + if (egName.size() <= signalSuffix.size() || egName.substr + (egName.size() - signalSuffix.size()) != signalSuffix) continue; + signal = signals.find(pyId); + if (signal == signals.end()) { + signals[pyId].status = -1; + signal = signals.find(pyId); + } + signal->second.egId = egId; + signal->second.bfs = vector(2, 0); + if (!final) continue; + + // Get the signal and background decay modes. + vector egList = egTable.getDecayTable(); + int sigIdx = egId.getAlias(); + int bkgIdx = EvtPDL::evtIdFromStdHep(pyId).getAlias(); + if (sigIdx > (int)egList.size() || bkgIdx > (int)egList.size()) continue; + EvtParticleDecayList sigModes = egList[sigIdx]; + EvtParticleDecayList bkgModes = egList[bkgIdx]; + EvtParticleDecayList allModes = egList[bkgIdx]; + double sum(0); + + // Sum signal branching fractions. + for (int iMode = 0; iMode < sigModes.getNMode(); ++iMode) { + EvtDecayBase *mode = sigModes.getDecayModel(iMode); + if (!mode) continue; + signal->second.bfs[0] += mode->getBranchingFraction(); + sum += mode->getBranchingFraction(); + } + + // Sum remaining background branching fractions. + for (int iMode = 0; iMode < allModes.getNMode(); ++iMode) { + EvtDecayBase *mode = allModes.getDecayModel(iMode); + if (!mode) continue; + bool match(false); + for (int jMode = 0; jMode < sigModes.getNMode(); ++jMode) + if (mode->matchingDecay(*(sigModes.getDecayModel(jMode)))) { + match = true; break;} + if (match) bkgModes.removeMode(mode); + else { + signal->second.map.push_back(iMode); + signal->second.bfs[1] += mode->getBranchingFraction(); + sum += mode->getBranchingFraction(); + } + } + signal->second.modes = bkgModes; + for (int iBf = 0; iBf < (int)signal->second.bfs.size(); ++iBf) + signal->second.bfs[iBf] /= sum; + } + if (final) updated = true; + +} + +//-------------------------------------------------------------------------- + +// Update the Pythia event record with an EvtGen decay tree. + +// The production vertex of each particle (which can also be obtained +// in EvtGen via EvtParticle::get4Pos()) is set by the decay vertex of +// its mother, which in turn is calculated from the mother's +// lifetime. The status code 93 is used to indicate an external decay, +// while the status code 94 is used to indicate an oscillated external +// decay. If the progenitor has a single daughter with the same ID, +// this daughter is used as the progenitor. This is used to prevent +// double oscillations. + +// If the arguments after egPro are no NULL and a particle in the +// decay tree is a signal particle, the decay for this particle is +// removed and the particle is stored as a signal candidate in the +// pySigs and egSigs vectors, to be decayed later. However, if any of +// these arguments is NULL then the entire tree is written. + +void EvtGenDecays::updateEvent(Particle *pyPro, EvtParticle *egPro, + vector *pySigs, vector *egSigs, + vector *bfs, double *wgt) { + + // Set up the mother vector. + if (!pythiaPtr) return; + EvtParticle* egMom = egPro; + if (egPro->getNDaug() == 1 && egPro->getPDGId() == + egPro->getDaug(0)->getPDGId()) egMom = egPro->getDaug(0); + Event &event = pythiaPtr->event; + vector< pair > + moms(1, pair(egMom, pyPro->index())); + + // Loop over the mothers. + while (moms.size() != 0) { + + // Check if particle can decay. + egMom = moms.back().first; + int iMom = moms.back().second; + Particle* pyMom = &event[iMom]; + moms.pop_back(); + if (!checkVertex(pyMom)) continue; + bool osc(checkOsc(egMom)); + + // Set the children of the mother. + pyMom->daughters(event.size(), event.size() + egMom->getNDaug() - 1); + pyMom->statusNeg(); + Vec4 vProd = pyMom->vDec(); + for (int iDtr = 0 ; iDtr < (int)egMom->getNDaug(); ++iDtr) { + EvtParticle *egDtr = egMom->getDaug(iDtr); + int id = egDtr->getPDGId(); + EvtVector4R p = egDtr->getP4Lab(); + int idx = event.append(id, 93, iMom, 0, 0, 0, 0, 0, p.get(1), + p.get(2), p.get(3), p.get(0), egDtr->mass()); + Particle *pyDtr = &event.back(); + pyDtr->vProd(vProd); + pyDtr->tau(egDtr->getLifetime()); + if (pySigs && egSigs && bfs && wgt && checkSignal(pyDtr)) { + pySigs->push_back(pyDtr->index()); + egSigs->push_back(egDtr); + bfs->push_back(signal->second.bfs[0]); + (*wgt) *= 1 - bfs->back(); + egDtr->deleteDaughters(); + } + if (osc) pyDtr->status(94); + if (egDtr->getNDaug() > 0) + moms.push_back(pair(egDtr, idx)); + } + } + +} + +//-------------------------------------------------------------------------- + +// Check if a particle can decay. + +// Modified slightly from ParticleDecays::checkVertex. + +bool EvtGenDecays::checkVertex(Particle *pyPro) { + if (!limitDecay) return true; + if (limitTau0 && pyPro->tau0() > tau0Max) return false; + if (limitTau && pyPro->tau() > tauMax) return false; + if (limitRadius && pow2(pyPro->xDec()) + pow2(pyPro->yDec()) + + pow2(pyPro->zDec()) > pow2(rMax)) return false; + if (limitCylinder && (pow2(pyPro->xDec()) + pow2(pyPro->yDec()) + > pow2(xyMax) || abs(pyPro->zDec()) > zMax) ) return false; + return true; +} + +//-------------------------------------------------------------------------- + +// Check if a particle is signal. + +bool EvtGenDecays::checkSignal(Particle *pyPro) { + signal = signals.find(pyPro->id()); + if (signal != signals.end() && (signal->second.status < 0 || + signal->second.status == pyPro->status())) return true; + else return false; +} + +//-------------------------------------------------------------------------- + +// Check if an EvtGen particle has oscillated. + +// The criteria defined here for oscillation is a single daughter but +// with a different ID from the mother. + +bool EvtGenDecays::checkOsc(EvtParticle *egPro) { + if (egPro && egPro->getNDaug() == 1 && + egPro->getPDGId() != egPro->getDaug(0)->getPDGId()) return true; + else return false; +} + +//========================================================================== + +} // end namespace Pythia8 + +#endif // CMSSW_Pythia8_EvtGen_Patched_H diff --git a/GeneratorInterface/Pythia8Interface/plugins/HMC3Py8EGun.cc b/GeneratorInterface/Pythia8Interface/plugins/HMC3Py8EGun.cc index 1863b84327e9b..c2a284cd68f6a 100644 --- a/GeneratorInterface/Pythia8Interface/plugins/HMC3Py8EGun.cc +++ b/GeneratorInterface/Pythia8Interface/plugins/HMC3Py8EGun.cc @@ -8,7 +8,7 @@ // EvtGen plugin // -#include "Pythia8Plugins/EvtGen.h" +#include "GeneratorInterface/Pythia8Interface/interface/Pythia8EvtGenPatched.h" namespace gen { diff --git a/GeneratorInterface/Pythia8Interface/plugins/Pythia8Hadronizer.cc b/GeneratorInterface/Pythia8Interface/plugins/Pythia8Hadronizer.cc index 4dfad6481224c..6494f327cdbf8 100644 --- a/GeneratorInterface/Pythia8Interface/plugins/Pythia8Hadronizer.cc +++ b/GeneratorInterface/Pythia8Interface/plugins/Pythia8Hadronizer.cc @@ -48,7 +48,7 @@ using namespace Pythia8; // EvtGen plugin // -#include "Pythia8Plugins/EvtGen.h" +#include "GeneratorInterface/Pythia8Interface/interface/Pythia8EvtGenPatched.h" #include "FWCore/AbstractServices/interface/RandomNumberGenerator.h" #include "FWCore/Concurrency/interface/SharedResourceNames.h" diff --git a/GeneratorInterface/Pythia8Interface/plugins/Pythia8HepMC3Hadronizer.cc b/GeneratorInterface/Pythia8Interface/plugins/Pythia8HepMC3Hadronizer.cc index ab5683d735dae..607ff0889aa55 100644 --- a/GeneratorInterface/Pythia8Interface/plugins/Pythia8HepMC3Hadronizer.cc +++ b/GeneratorInterface/Pythia8Interface/plugins/Pythia8HepMC3Hadronizer.cc @@ -47,7 +47,7 @@ using namespace Pythia8; // EvtGen plugin // -#include "Pythia8Plugins/EvtGen.h" +#include "GeneratorInterface/Pythia8Interface/interface/Pythia8EvtGenPatched.h" #include "FWCore/AbstractServices/interface/RandomNumberGenerator.h" #include "FWCore/Concurrency/interface/SharedResourceNames.h" diff --git a/GeneratorInterface/Pythia8Interface/src/Py8GunBase.cc b/GeneratorInterface/Pythia8Interface/src/Py8GunBase.cc index b560d45205be8..2957c936e2fee 100644 --- a/GeneratorInterface/Pythia8Interface/src/Py8GunBase.cc +++ b/GeneratorInterface/Pythia8Interface/src/Py8GunBase.cc @@ -6,7 +6,7 @@ // EvtGen plugin // -#include "Pythia8Plugins/EvtGen.h" +#include "GeneratorInterface/Pythia8Interface/interface/Pythia8EvtGenPatched.h" using namespace Pythia8; From 8a349144b560ea3031d52504054da02d229e1351 Mon Sep 17 00:00:00 2001 From: Ka Wa Ho Date: Sat, 16 May 2026 16:19:07 -0500 Subject: [PATCH 2/8] move to thread_local --- .../plugins/EvtGen/EvtGenInterface.cc | 64 ++++++++++++------- .../plugins/EvtGen/EvtGenInterface.h | 8 ++- .../interface/Pythia8EvtGenPatched.h | 13 +++- 3 files changed, 57 insertions(+), 28 deletions(-) diff --git a/GeneratorInterface/EvtGenInterface/plugins/EvtGen/EvtGenInterface.cc b/GeneratorInterface/EvtGenInterface/plugins/EvtGen/EvtGenInterface.cc index d7b6deaa82d7b..227ca0e13ff84 100644 --- a/GeneratorInterface/EvtGenInterface/plugins/EvtGen/EvtGenInterface.cc +++ b/GeneratorInterface/EvtGenInterface/plugins/EvtGen/EvtGenInterface.cc @@ -37,6 +37,7 @@ #include "EvtGenBase/EvtId.hh" #include "EvtGenBase/EvtDecayTable.hh" #include "EvtGenBase/EvtPDL.hh" +#include "EvtGenBase/EvtRandom.hh" #include "EvtGenBase/EvtParticle.hh" #include "EvtGenBase/EvtParticleFactory.hh" #include "EvtGenBase/EvtHepMCEvent.hh" @@ -48,12 +49,15 @@ #include #include +thread_local std::unique_ptr gen::EvtGenInterface::m_EvtGen{}; + using namespace gen; using namespace edm; + EvtGenInterface::EvtGenInterface(const ParameterSet& pset) { - fPSet = new ParameterSet(pset); - the_engine = new myEvtRandomEngine(nullptr); + fPSet = std::make_unique(pset); + the_engine = std::make_unique(nullptr); } void EvtGenInterface::SetDefault_m_PDGs() { @@ -283,12 +287,14 @@ void EvtGenInterface::SetDefault_m_PDGs() { } } +// Destructor is a no-op: m_EvtGen is a thread_local std::unique_ptr that cleans up at +// thread exit. EvtGenInterface::~EvtGenInterface() {} -void EvtGenInterface::init() { - // flags for pythia8 - fSpecialSettings.push_back("Pythia8:ParticleDecays:mixB = off"); - // +void EvtGenInterface::ensureEvtGenOnThread() { + // This thread already has its EvtGen + if (m_EvtGen) + return; edm::FileInPath decay_table(fPSet->getParameter("decay_table")); edm::FileInPath pdt(fPSet->getParameter("particle_property_file")); @@ -304,8 +310,8 @@ void EvtGenInterface::init() { std::getenv("PYTHIA8DATA"); // Specify the pythia xml data directory to use the default PYTHIA8DATA location if (tmp == nullptr) { - edm::LogError("EvtGenInterface::~EvtGenInterface") - << "EvtGenInterface::init() PYTHIA8DATA not defined. Terminating program "; + edm::LogError("EvtGenInterface::ensureEvtGenOnThread") + << "PYTHIA8DATA not defined. Terminating program "; exit(0); } @@ -338,21 +344,12 @@ void EvtGenInterface::init() { std::list::iterator it = userModels.begin(); std::advance(it, i); TString name = (*it)->getName(); - edm::LogInfo("EvtGenInterface::~EvtGenInterface") << "Adding user model: " << name; + edm::LogInfo("EvtGenInterface::ensureEvtGenOnThread") << "Adding user model: " << name; myExtraModels.push_back(*it); } - // Set up the incoherent (1) or coherent (0) B mixing option - BmixingOption = fPSet->getUntrackedParameter("B_Mixing", 1); - if (BmixingOption != 0 && BmixingOption != 1) { - throw cms::Exception("Configuration") << "EvtGenProducer requires B_Mixing to be 0 (coherent) or 1 (incoherent) \n" - "Please fix this in your configuration."; - } - - ////////////////////////////////////////////////////////////////////////////////////////////////////////// - // Create the EvtGen generator object, passing the external generators - m_EvtGen = new EvtGen( - decay_table.fullPath().c_str(), pdt.fullPath().c_str(), the_engine, radCorrEngine, &myExtraModels, BmixingOption); + m_EvtGen = std::make_unique( + decay_table.fullPath().c_str(), pdt.fullPath().c_str(), the_engine.get(), radCorrEngine, &myExtraModels, BmixingOption); // Add additional user information if (fPSet->exists("user_decay_file")) { @@ -369,9 +366,8 @@ void EvtGenInterface::init() { int tmp_creation = mkstemp(user_decay_tmp); FILE* tmpf = std::fopen(user_decay_tmp, "w"); if (!tmpf || (tmp_creation == -1)) { - edm::LogError("EvtGenInterface::~EvtGenInterface") - << "EvtGenInterface::init() fails when trying to open a temporary file for embedded user.dec. Terminating " - "program "; + edm::LogError("EvtGenInterface::ensureEvtGenOnThread") + << "fails when trying to open a temporary file for embedded user.dec. Terminating program "; exit(0); } for (unsigned int i = 0; i < user_decay_lines.size(); i++) { @@ -381,6 +377,23 @@ void EvtGenInterface::init() { std::fclose(tmpf); m_EvtGen->readUDecay(user_decay_tmp); } +} + +void EvtGenInterface::init() { + // flags for pythia8 + fSpecialSettings.push_back("Pythia8:ParticleDecays:mixB = off"); + + // Set up the incoherent (1) or coherent (0) B mixing option. Must be set before + // ensureEvtGenOnThread() because the EvtGen ctor consumes it. + BmixingOption = fPSet->getUntrackedParameter("B_Mixing", 1); + if (BmixingOption != 0 && BmixingOption != 1) { + throw cms::Exception("Configuration") << "EvtGenProducer requires B_Mixing to be 0 (coherent) or 1 (incoherent) \n" + "Please fix this in your configuration."; + } + + // Construct the EvtGen instance for this thread (no-op if already constructed). Streams + // running decay() later on a different thread will lazy-construct on that thread instead. + ensureEvtGenOnThread(); // setup pdgid which the generator/hadronizer should not decay if (fPSet->exists("operates_on_particles")) { @@ -451,6 +464,10 @@ void EvtGenInterface::init() { } HepMC::GenEvent* EvtGenInterface::decay(HepMC::GenEvent* evt) { + // If the framework scheduled this stream onto a thread that has not yet built an EvtGen, + // construct one now. No-op on threads that already initialised in init(). + ensureEvtGenOnThread(); + if (the_engine->engine() == nullptr) { throw edm::Exception(edm::errors::LogicError) << "The EvtGen code attempted to use a random number engine while\n" @@ -615,6 +632,7 @@ void EvtGenInterface::update_particles(HepMC::GenParticle* partHep, HepMC::GenEv void EvtGenInterface::setRandomEngine(CLHEP::HepRandomEngine* v) { the_engine->setRandomEngine(v); + EvtRandom::setRandomEngine(the_engine.get()); } double EvtGenInterface::flat() { diff --git a/GeneratorInterface/EvtGenInterface/plugins/EvtGen/EvtGenInterface.h b/GeneratorInterface/EvtGenInterface/plugins/EvtGen/EvtGenInterface.h index f3ae6fa5153bc..bb5f810784254 100644 --- a/GeneratorInterface/EvtGenInterface/plugins/EvtGen/EvtGenInterface.h +++ b/GeneratorInterface/EvtGenInterface/plugins/EvtGen/EvtGenInterface.h @@ -52,8 +52,10 @@ namespace gen { bool findLastinChain(HepMC::GenParticle*& p); bool hasnoDaughter(HepMC::GenParticle* p); void go_through_daughters(EvtParticle* part); + void ensureEvtGenOnThread(); - EvtGen* m_EvtGen; // EvtGen main object + // One EvtGen per thread + static thread_local std::unique_ptr m_EvtGen; std::vector forced_id; // EvtGen Id's of particles which are to be forced by EvtGen std::vector forced_pdgids; // PDG Id's of particles which are to be forced by EvtGen @@ -65,9 +67,9 @@ namespace gen { std::vector polarize_pol; std::map polarizations; int BmixingOption = 1; - edm::ParameterSet* fPSet; + std::unique_ptr fPSet; - myEvtRandomEngine* the_engine; + std::unique_ptr the_engine; }; } // namespace gen #endif diff --git a/GeneratorInterface/Pythia8Interface/interface/Pythia8EvtGenPatched.h b/GeneratorInterface/Pythia8Interface/interface/Pythia8EvtGenPatched.h index c0cec08a0fe9a..b95ef93b43def 100644 --- a/GeneratorInterface/Pythia8Interface/interface/Pythia8EvtGenPatched.h +++ b/GeneratorInterface/Pythia8Interface/interface/Pythia8EvtGenPatched.h @@ -1,7 +1,16 @@ -// Local patched copy of Pythia8Plugins/EvtGen.h adapted to EvtGen 03.x API. +// EvtGen.h is a part of the PYTHIA event generator. +// Copyright (C) 2026 Torbjorn Sjostrand. +// PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details. +// Please respect the MCnet Guidelines, see GUIDELINES for details. +// Author: Philip Ilten. + +// Local patched copy of Pythia8Plugins/EvtGen.h adapted to EvtGen 03.x.x API. // Diverges from the upstream pythia8 header only in: // - EvtGenRandom now overrides setSeed(unsigned long)/lastSeed() pure virtuals // - EvtDecayTable::getInstance() returns EvtDecayTable& (not pointer) +// EvtGen.h is a part of the PYTHIA event generator. + +// This file contains an EvtGen interface. HepMC and EvtGen must be enabled. #ifndef CMSSW_Pythia8_EvtGen_Patched_H #define CMSSW_Pythia8_EvtGen_Patched_H @@ -33,7 +42,7 @@ class EvtGenRandom : public EvtRandomEngine { // Return a random number. double random() override {if (rndmPtr) return rndmPtr->flat(); else return -1.0;} - // EvtGen 03 pure virtuals. Record-only: Pythia owns the actual RNG state. + // EvtGen 03.x.x pure virtuals. Record-only: Pythia owns the actual RNG state. void setSeed(unsigned long int seed) override {lastSeedVal = seed;} unsigned long int lastSeed() const override {return lastSeedVal;} From 72ed0d79874786b240ca5633a0f7403233330237 Mon Sep 17 00:00:00 2001 From: Ka Wa Ho Date: Mon, 18 May 2026 13:50:56 -0500 Subject: [PATCH 3/8] Apply ThreadHandoff to ConcurrentExternalDecayDriver --- .../ExternalDecays/BuildFile.xml | 1 + .../interface/ConcurrentExternalDecayDriver.h | 4 +- .../interface/EvtGenThreadOwner.h | 34 +++++ .../src/ConcurrentExternalDecayDriver.cc | 118 ++++-------------- 4 files changed, 63 insertions(+), 94 deletions(-) create mode 100644 GeneratorInterface/ExternalDecays/interface/EvtGenThreadOwner.h diff --git a/GeneratorInterface/ExternalDecays/BuildFile.xml b/GeneratorInterface/ExternalDecays/BuildFile.xml index 64b095a3fc707..887efbc72d59e 100644 --- a/GeneratorInterface/ExternalDecays/BuildFile.xml +++ b/GeneratorInterface/ExternalDecays/BuildFile.xml @@ -1,6 +1,7 @@ + diff --git a/GeneratorInterface/ExternalDecays/interface/ConcurrentExternalDecayDriver.h b/GeneratorInterface/ExternalDecays/interface/ConcurrentExternalDecayDriver.h index 80f45590dfa09..244d95cee8199 100644 --- a/GeneratorInterface/ExternalDecays/interface/ConcurrentExternalDecayDriver.h +++ b/GeneratorInterface/ExternalDecays/interface/ConcurrentExternalDecayDriver.h @@ -21,6 +21,7 @@ namespace lhef { namespace gen { + class EvtGenThreadOwner; class EvtGenInterfaceBase; class TauolaInterfaceBase; class PhotosInterfaceBase; @@ -45,8 +46,9 @@ namespace gen { private: bool fIsInitialized; + std::unique_ptr fThreadOwner; + std::unique_ptr fEvtGenInterface; //std::unique_ptr fTauolaInterface; - //std::unique_ptr fEvtGenInterface; //std::unique_ptr fPhotosInterface; std::vector fPDGs; std::vector fSpecialSettings; diff --git a/GeneratorInterface/ExternalDecays/interface/EvtGenThreadOwner.h b/GeneratorInterface/ExternalDecays/interface/EvtGenThreadOwner.h new file mode 100644 index 0000000000000..bb0baee995e1d --- /dev/null +++ b/GeneratorInterface/ExternalDecays/interface/EvtGenThreadOwner.h @@ -0,0 +1,34 @@ +#ifndef gen_EvtGenThreadOwner_h +#define gen_EvtGenThreadOwner_h + +// Pins one dedicated pthread per stream and marshals every EvtGen call onto it. +// Mirrors SimG4Core's OscarMTProducer + omt::ThreadHandoff pattern + +#include "SimG4Core/Application/interface/ThreadHandoff.h" +#include "FWCore/ServiceRegistry/interface/ServiceRegistry.h" + +namespace gen { + + class EvtGenThreadOwner { + public: + explicit EvtGenThreadOwner(int stackSize) : m_handoff{stackSize} {} + + EvtGenThreadOwner(const EvtGenThreadOwner&) = delete; + EvtGenThreadOwner& operator=(const EvtGenThreadOwner&) = delete; + + template + void run(F&& f) { + auto token = edm::ServiceRegistry::instance().presentToken(); + m_handoff.runAndWait([token, &f]() { + edm::ServiceRegistry::Operate guard{token}; + f(); + }); + } + + private: + omt::ThreadHandoff m_handoff; + }; + +} + +#endif diff --git a/GeneratorInterface/ExternalDecays/src/ConcurrentExternalDecayDriver.cc b/GeneratorInterface/ExternalDecays/src/ConcurrentExternalDecayDriver.cc index 57ccbe49710af..c86db8ecc2279 100644 --- a/GeneratorInterface/ExternalDecays/src/ConcurrentExternalDecayDriver.cc +++ b/GeneratorInterface/ExternalDecays/src/ConcurrentExternalDecayDriver.cc @@ -1,4 +1,5 @@ #include "GeneratorInterface/ExternalDecays/interface/ConcurrentExternalDecayDriver.h" +#include "GeneratorInterface/ExternalDecays/interface/EvtGenThreadOwner.h" #include "GeneratorInterface/Core/interface/FortranInstance.h" #include "GeneratorInterface/EvtGenInterface/interface/EvtGenFactory.h" @@ -20,100 +21,55 @@ using namespace gen; using namespace edm; -ConcurrentExternalDecayDriver::ConcurrentExternalDecayDriver(const ParameterSet& pset) : fIsInitialized(false) { +// Matches SimG4Core's OscarMTProducer default (`workerThreadStackSize`). +constexpr int kEvtGenThreadStackSize = 10 * 1024 * 1024; + +ConcurrentExternalDecayDriver::ConcurrentExternalDecayDriver(const ParameterSet& pset) + : fIsInitialized(false), fThreadOwner(std::make_unique(kEvtGenThreadStackSize)) { std::vector extGenNames = pset.getParameter >("parameterSets"); for (unsigned int ip = 0; ip < extGenNames.size(); ++ip) { const std::string& curSet = extGenNames[ip]; - throw cms::Exception("ThreadUnsafeDecayer") << "The decayer " << curSet << " is not thread-friendly."; - /* - if (curSet == "EvtGen") { - fEvtGenInterface = std::unique_ptr( - EvtGenFactory::get()->create("EvtGen", pset.getUntrackedParameter(curSet))); - exSharedResources.emplace_back(edm::SharedResourceNames::kEvtGen); - exSharedResources.emplace_back(edm::SharedResourceNames::kPythia6); - exSharedResources.emplace_back(gen::FortranInstance::kFortranInstance); - } else if (curSet == "EvtGen1" || curSet == "EvtGen130") { - fEvtGenInterface = std::unique_ptr( - EvtGenFactory::get()->create("EvtGen130", pset.getUntrackedParameter(curSet))); - exSharedResources.emplace_back(edm::SharedResourceNames::kEvtGen); - exSharedResources.emplace_back(edm::SharedResourceNames::kPythia8); - exSharedResources.emplace_back(edm::SharedResourceNames::kTauola); - exSharedResources.emplace_back(edm::SharedResourceNames::kPhotos); - exSharedResources.emplace_back(gen::FortranInstance::kFortranInstance); - } else if (curSet == "Tauola" || curSet == "Tauolapp" || curSet == "Tauolapp114") { - fTauolaInterface = std::unique_ptr( - TauolaFactory::get()->create("Tauolapp114", pset.getUntrackedParameter(curSet))); - fPhotosInterface = std::unique_ptr( - PhotosFactory::get()->create("Photos2155", pset.getUntrackedParameter(curSet))); - fPhotosInterface->configureOnlyFor(15); - fPhotosInterface->avoidTauLeptonicDecays(); - exSharedResources.emplace_back(edm::SharedResourceNames::kTauola); - exSharedResources.emplace_back(edm::SharedResourceNames::kPhotos); - } else if (curSet == "Photos" || curSet == "Photos2155") { - if (!fPhotosInterface) { - fPhotosInterface = std::unique_ptr( - PhotosFactory::get()->create("Photos2155", pset.getUntrackedParameter(curSet))); - exSharedResources.emplace_back(edm::SharedResourceNames::kPhotos); - } - } else if (curSet == "Photospp" || curSet == "Photospp356") { - if (!fPhotosInterface) { - fPhotosInterface = std::unique_ptr( - PhotosFactory::get()->create("Photospp356", pset.getUntrackedParameter(curSet))); - exSharedResources.emplace_back(edm::SharedResourceNames::kPhotos); - } + if (curSet == "EvtGen1" || curSet == "EvtGen130") { + fThreadOwner->run([&]() { + fEvtGenInterface = std::unique_ptr( + EvtGenFactory::get()->create("EvtGen130", pset.getUntrackedParameter(curSet))); + }); + } else { + throw cms::Exception("ThreadUnsafeDecayer") << "The decayer " << curSet << " is not thread-friendly."; } - */ } } -ConcurrentExternalDecayDriver::~ConcurrentExternalDecayDriver() = default; +ConcurrentExternalDecayDriver::~ConcurrentExternalDecayDriver() { + if (fEvtGenInterface) { + fThreadOwner->run([&]() { fEvtGenInterface.reset(); }); + } +} HepMC::GenEvent* ConcurrentExternalDecayDriver::decay(HepMC::GenEvent* evt, lhef::LHEEvent* lheEvent) { - /* if (fTauolaInterface) - fTauolaInterface->SetLHE(lheEvent); */ return decay(evt); } HepMC::GenEvent* ConcurrentExternalDecayDriver::decay(HepMC::GenEvent* evt) { if (!fIsInitialized) return evt; - /* - if (fEvtGenInterface) { - evt = fEvtGenInterface->decay(evt); - if (!evt) - return nullptr; - } - if (fTauolaInterface) { - evt = fTauolaInterface->decay(evt); + if (fEvtGenInterface) { + fThreadOwner->run([&]() { evt = fEvtGenInterface->decay(evt); }); if (!evt) return nullptr; } - if (fPhotosInterface) { - evt = fPhotosInterface->apply(evt); - if (!evt) - return nullptr; - } - */ return evt; } void ConcurrentExternalDecayDriver::init(const edm::EventSetup& es) { if (fIsInitialized) return; - /* - if (fTauolaInterface) { - fTauolaInterface->init(es); - for (std::vector::const_iterator i = fTauolaInterface->operatesOnParticles().begin(); - i != fTauolaInterface->operatesOnParticles().end(); - i++) - fPDGs.push_back(*i); - } if (fEvtGenInterface) { - fEvtGenInterface->init(); + fThreadOwner->run([&]() { fEvtGenInterface->init(); }); for (std::vector::const_iterator i = fEvtGenInterface->operatesOnParticles().begin(); i != fEvtGenInterface->operatesOnParticles().end(); i++) @@ -123,39 +79,15 @@ void ConcurrentExternalDecayDriver::init(const edm::EventSetup& es) { } } - if (fPhotosInterface) { - fPhotosInterface->init(); - // for tauola++ - if (fPhotosInterface) { - for (unsigned int iss = 0; iss < fPhotosInterface->specialSettings().size(); iss++) { - fSpecialSettings.push_back(fPhotosInterface->specialSettings()[iss]); - } - } - } - */ - fIsInitialized = true; return; } -void ConcurrentExternalDecayDriver::statistics() const { - /* if (fTauolaInterface) - fTauolaInterface->statistics(); - if (fPhotosInterface) - fPhotosInterface->statistics(); - */ - // similar for EvtGen if needed - return; -} +void ConcurrentExternalDecayDriver::statistics() const { return; } void ConcurrentExternalDecayDriver::setRandomEngine(CLHEP::HepRandomEngine* v) { - /* - if (fTauolaInterface) - fTauolaInterface->setRandomEngine(v); - if (fEvtGenInterface) - fEvtGenInterface->setRandomEngine(v); - if (fPhotosInterface) - fPhotosInterface->setRandomEngine(v); - */ + if (fEvtGenInterface) { + fThreadOwner->run([&]() { fEvtGenInterface->setRandomEngine(v); }); + } } From a6dfc854f144bca819033911cc112a814ba79437 Mon Sep 17 00:00:00 2001 From: Ka Wa Ho Date: Mon, 18 May 2026 15:42:02 -0500 Subject: [PATCH 4/8] move pythia patch to cmsdist --- .../interface/Pythia8EvtGenPatched.h | 633 ------------------ .../Pythia8Interface/plugins/HMC3Py8EGun.cc | 2 +- .../plugins/Pythia8Hadronizer.cc | 2 +- .../plugins/Pythia8HepMC3Hadronizer.cc | 2 +- .../Pythia8Interface/src/Py8GunBase.cc | 2 +- 5 files changed, 4 insertions(+), 637 deletions(-) delete mode 100644 GeneratorInterface/Pythia8Interface/interface/Pythia8EvtGenPatched.h diff --git a/GeneratorInterface/Pythia8Interface/interface/Pythia8EvtGenPatched.h b/GeneratorInterface/Pythia8Interface/interface/Pythia8EvtGenPatched.h deleted file mode 100644 index b95ef93b43def..0000000000000 --- a/GeneratorInterface/Pythia8Interface/interface/Pythia8EvtGenPatched.h +++ /dev/null @@ -1,633 +0,0 @@ -// EvtGen.h is a part of the PYTHIA event generator. -// Copyright (C) 2026 Torbjorn Sjostrand. -// PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details. -// Please respect the MCnet Guidelines, see GUIDELINES for details. -// Author: Philip Ilten. - -// Local patched copy of Pythia8Plugins/EvtGen.h adapted to EvtGen 03.x.x API. -// Diverges from the upstream pythia8 header only in: -// - EvtGenRandom now overrides setSeed(unsigned long)/lastSeed() pure virtuals -// - EvtDecayTable::getInstance() returns EvtDecayTable& (not pointer) -// EvtGen.h is a part of the PYTHIA event generator. - -// This file contains an EvtGen interface. HepMC and EvtGen must be enabled. - -#ifndef CMSSW_Pythia8_EvtGen_Patched_H -#define CMSSW_Pythia8_EvtGen_Patched_H - -#include "Pythia8/Pythia.h" -#include "EvtGen/EvtGen.hh" -#include "EvtGenBase/EvtRandomEngine.hh" -#include "EvtGenBase/EvtParticle.hh" -#include "EvtGenBase/EvtParticleFactory.hh" -#include "EvtGenBase/EvtPDL.hh" -#include "EvtGenBase/EvtDecayTable.hh" -#include "EvtGenBase/EvtParticleDecayList.hh" -#include "EvtGenBase/EvtDecayBase.hh" -#include "EvtGenExternal/EvtExternalGenList.hh" - -namespace Pythia8 { - -//========================================================================== - -// A class to wrap the Pythia random number generator for use by EvtGen. - -class EvtGenRandom : public EvtRandomEngine { - -public: - - // Constructor. - EvtGenRandom(Rndm *rndmPtrIn) {rndmPtr = rndmPtrIn;} - - // Return a random number. - double random() override {if (rndmPtr) return rndmPtr->flat(); else return -1.0;} - - // EvtGen 03.x.x pure virtuals. Record-only: Pythia owns the actual RNG state. - void setSeed(unsigned long int seed) override {lastSeedVal = seed;} - unsigned long int lastSeed() const override {return lastSeedVal;} - - // The random number pointer. - Rndm *rndmPtr; - -private: - unsigned long int lastSeedVal{0}; - -}; - -//========================================================================== - -// A class to perform decays via the external EvtGen decay program, -// see http://evtgen.warwick.ac.uk/, the program manual provided with -// the EvtGen distribution, and D. J. Lange, -// Nucl. Instrum. Meth. A462, 152 (2001) for details. - -// EvtGen performs a series of decays from some initial particle -// decay, rather than just a single decay, and so EvtGen cannot be -// interfaced through the standard external DecayHandler class without -// considerable complication. Consequently, EvtGen is called on the -// complete event record after all steps of Pythia are completed. - -// Oftentimes a specific "signal" decay is needed to occur once in an -// event, and all other decays performed normally. This is possible -// via reading in a user decay file (with readDecayFile) and creating -// aliased particles with names ending with signalSuffix. By default, -// this is "_SIGNAL". When decay() is called, all particles in the -// Pythia event record that are of the same types as the signal -// particles are collected. One is selected at random and decayed via -// the channel(s) defined for that aliased signal particle. All other -// particles are decayed normally. The weight for the event is -// calculated and returned. - -// It is also possible to specify a status needed to consider a -// particle as a signal candidate. This can be done by modifying the -// signals map, e.g. if the tau- is a signal candidate, then -// EvtGenDecays.signals[15].status = 201 -// will only only select as candidates any tau- with this status. This -// allows the event record to be changed before decays, so only -// certain particles are selected as possible signal candidates -// (e.g. passing kinematic requirements). - -// Please note that particles produced from a signal candidate decay -// are not searched for additional signal candidates. This means that -// if B0 and tau- have been designated as signal, then a tau- from a -// W- decay would be a signal candidate, while a tau- from a B0 decay -// would not. This restriction arises from the additional complexity -// of allowing recursive signal decays. The following statuses are -// used: 93 for particles decayed with EvtGen, 94 for particles -// oscillated with EvtGen, 95 for signal particles, and 96 for signal -// particles from an oscillation. - -class EvtGenDecays { - -public: - - // Constructor. - EvtGenDecays(Pythia *pythiaPtrIn, string decayFile, string particleDataFile, - EvtExternalGenList *extPtrIn = 0, EvtAbsRadCorr *fsrPtrIn = 0, - int mixing = 1, bool xml = false, bool limit = true, - bool extUse = true, bool fsrUse = true); - - // Destructor. - ~EvtGenDecays() { - if (evtgen) delete evtgen; - if (extOwner && extPtr) delete extPtr; - if (fsrOwner && fsrPtr) delete fsrPtr; - } - - // Perform all decays and return the event weight. - double decay(); - - // Stop EvtGen decaying a particle. - void exclude(int id) {excIds.insert(id);} - - // Update the Pythia particle database from EvtGen. - void updatePythia(); - - // Update the EvtGen particle database from Pythia. - void updateEvtGen(); - - // Read an EvtGen user decay file. - void readDecayFile(string decayFile, bool xml = false) { - evtgen->readUDecay(decayFile.c_str(), xml); updateData();} - - // External model pointer and FSR model pointer. - bool extOwner, fsrOwner; - EvtExternalGenList *extPtr; - EvtAbsRadCorr *fsrPtr; - std::list models; - - // Map of signal particle info. - struct Signal {int status; EvtId egId; vector bfs; vector map; - EvtParticleDecayList modes;}; - map signals; - - // The suffix indicating an EvtGen particle or alias is signal. - string signalSuffix; - -protected: - - // Update the particles to decay with EvtGen, and the selected signals. - void updateData(bool final = false); - - // Update the Pythia event record with an EvtGen decay tree. - void updateEvent(Particle *pyPro, EvtParticle *egPro, - vector *pySigs = 0, vector *egSigs = 0, - vector *bfs = 0, double *wgt = 0); - - // Check if a particle can decay. - bool checkVertex(Particle *pyPro); - - // Check if a particle is signal. - bool checkSignal(Particle *pyPro); - - // Check if an EvtGen particle has oscillated. - bool checkOsc(EvtParticle *egPro); - - // Number of times to try a decay sampling (constant). - static const int NTRYDECAY = 1000; - - // The pointer to the associated Pythia object. - Pythia *pythiaPtr; - - // Random number wrapper for EvtGen. - EvtGenRandom rndm; - - // The EvtGen object. - EvtGen *evtgen; - - // Set of particle IDs to include and exclude decays with EvtGen. - set incIds, excIds; - - // Flag whether the final particle update has been performed. - bool updated; - - // The selected signal iterator. - map::iterator signal; - - // Parameters used to check if a particle should decay (as set via Pythia). - double tau0Max, tauMax, rMax, xyMax, zMax; - bool limitTau0, limitTau, limitRadius, limitCylinder, limitDecay; - -}; - -//-------------------------------------------------------------------------- - -// The constructor. - -// The EvtGenDecays object is associated with a single Pythia -// instance. This is to ensure a consistent random number generator -// across the two, as well as any updates to particle data, etc. Note -// that if multiple EvtGenDecays objects exist, that they will modify -// one anothers particle databases due to the design of EvtGen. - -// This constructor also sets all particles to be decayed by EvtGen as -// stable within Pythia. The parameters within Pythia used to check if -// a particle should be decayed, as described in the "Particle Decays" -// section of the Pythia manual, are set. Note that if the variable -// "limit" is set to "false", then no check will be made before -// decaying a particle with EvtGen. - -// The constructor is designed to have the exact same form as the -// EvtGen constructor except for these five differences. -// (1) The first variable is the pointer to the Pythia object. -// (2) The third last argument is a flag to limit decays based on the -// Pythia criteria (based on the particle decay vertex). -// (3) The second last argument is a flag if external models should be -// passed to EvtGen (default is true). -// (4) The last argument is a flag if an FSR model should be passed -// to EvtGen (default is true). -// (5) No random engine pointer is passed, as this is obtained from -// Pythia. - -// pythiaPtrIn: the pointer to the associated Pythia object. -// decayFile: the name of the decay file to pass to EvtGen. -// particleDataFile: the name of the particle data file to pass to EvtGen. -// extPtrIn: the optional EvtExternalGenList pointer, this must be -// be provided if fsrPtrIn is provided to avoid double -// initializations. -// fsrPtrIn: the EvtAbsRadCorr pointer to pass to EvtGen. -// mixing: the mixing type to pass to EvtGen. -// xml: flag to use XML files to pass to EvtGen. -// limit: flag to limit particle decays based on Pythia criteria. -// extUse: flag to use external models with EvtGen. -// fsrUse: flag to use radiative correction engine with EvtGen. - -EvtGenDecays::EvtGenDecays(Pythia *pythiaPtrIn, string decayFile, - string particleDataFile, EvtExternalGenList *extPtrIn, - EvtAbsRadCorr *fsrPtrIn, int mixing, bool xml, bool limit, - bool extUse, bool fsrUse) : extPtr(extPtrIn), fsrPtr(fsrPtrIn), - signalSuffix("_SIGNAL"), pythiaPtr(pythiaPtrIn), rndm(&pythiaPtr->rndm), - updated(false) { - - // Initialize EvtGen. - if (!extPtr && fsrPtr) { - cout << "Error in EvtGenDecays::" - << "EvtGenDecays: extPtr is null but fsrPtr is provided\n"; - return; - } - if (extPtr) extOwner = false; - else {extOwner = true; extPtr = new EvtExternalGenList();} - if (fsrPtr) fsrOwner = false; - else {fsrOwner = true; fsrPtr = extPtr->getPhotosModel();} - models = extPtr->getListOfModels(); - evtgen = new EvtGen(decayFile.c_str(), particleDataFile.c_str(), - &rndm, fsrUse ? fsrPtr : 0, extUse ? &models : 0, mixing, xml); - - // Get the Pythia decay limits. - if (!pythiaPtr) return; - limitTau0 = pythiaPtr->settings.flag("ParticleDecays:limitTau0"); - tau0Max = pythiaPtr->settings.parm("ParticleDecays:tau0Max"); - limitTau = pythiaPtr->settings.flag("ParticleDecays:limitTau"); - tauMax = pythiaPtr->settings.parm("ParticleDecays:tauMax"); - limitRadius = pythiaPtr->settings.flag("ParticleDecays:limitRadius"); - rMax = pythiaPtr->settings.parm("ParticleDecays:rMax"); - limitCylinder = pythiaPtr->settings.flag("ParticleDecays:limitCylinder"); - xyMax = pythiaPtr->settings.parm("ParticleDecays:xyMax"); - zMax = pythiaPtr->settings.parm("ParticleDecays:zMax"); - limitDecay = limit && (limitTau0 || limitTau || - limitRadius || limitCylinder); - -} - -//-------------------------------------------------------------------------- - -// Perform all decays and return the event weight. - -// All particles in the event record that can be decayed by EvtGen are -// decayed. If a particle is a signal particle, then this is stored in -// a vector of signal particles. A signal particle is only stored if -// its status is the same as the status provided in the signals map. A -// negative status in the signal map indicates that all statuses -// should be accepted. After all signal particles are identified, one -// is randomly chosen and decayed as signal. The remainder are decayed -// normally. - -// Forcing a signal decay changes the weight of an event from unity, -// and so the relative event weight is returned, given the forced -// signal decay. A weight of 0 indicates no signal in the event, while -// a weight of -1 indicates something is wrong, e.g. either the Pythia -// or EvtGen pointers are not available or the number of tries has -// been exceeded. For the event weight to be valid, one should not -// change the absolute branching fractions in the signal and inclusive -// definitions, but rather just remove the unwanted decay channels -// from the signal decay definition. - -double EvtGenDecays::decay() { - - // Reset the signal and signal counters. - if (!pythiaPtr || !evtgen) return -1; - if (!updated) updateData(true); - - // Loop over all particles in the Pythia event. - Event &event = pythiaPtr->event; - vector pySigs; vector egSigs, egPrts; - vector bfs; double wgt(1.); - for (int iPro = 0; iPro < event.size(); ++iPro) { - - // Check particle is final and can be decayed by EvtGen. - Particle *pyPro = &event[iPro]; - if (!pyPro->isFinal()) continue; - if (incIds.find(pyPro->id()) == incIds.end()) continue; - if (pyPro->status() == 93 || pyPro->status() == 94) continue; - - // Decay the progenitor with EvtGen. - EvtParticle *egPro = EvtParticleFactory::particleFactory - (EvtPDL::evtIdFromStdHep(pyPro->id()), - EvtVector4R(pyPro->e(), pyPro->px(), pyPro->py(), pyPro->pz())); - egPrts.push_back(egPro); - egPro->setDiagonalSpinDensity(); - evtgen->generateDecay(egPro); - pyPro->tau(egPro->getLifetime()); - if (!checkVertex(pyPro)) continue; - - // Add oscillations to event record. - if (checkOsc(egPro)) - updateEvent(pyPro, egPro, &pySigs, &egSigs, &bfs, &wgt); - - // Undo decay if signal (duplicate to stop oscillations). - else if (checkSignal(pyPro)) { - pySigs.push_back(pyPro->index()); - egSigs.push_back(egPro); - bfs.push_back(signal->second.bfs[0]); - wgt *= 1 - bfs.back(); - egPro->deleteDaughters(); - EvtParticle *egDau = EvtParticleFactory::particleFactory - (EvtPDL::evtIdFromStdHep(pyPro->id()), - EvtVector4R(pyPro->e(), pyPro->px(), pyPro->py(), pyPro->pz())); - egDau->addDaug(egPro); - egDau->setDiagonalSpinDensity(); - - // If not signal, add to event record. - } else updateEvent(pyPro, egPro, &pySigs, &egSigs, &bfs, &wgt); - } - if (pySigs.size() == 0) { - for (int iPrt = 0; iPrt < (int)egPrts.size(); ++iPrt) - egPrts[iPrt]->deleteTree(); - return 0; - } - - // Determine the decays of the signal particles (signal or background). - vector modes; int force(-1), n(0); - for (int iTry = 1; iTry <= NTRYDECAY; ++iTry) { - modes.clear(); force = pythiaPtr->rndm.pick(bfs); n = 0; - for (int iSig = 0; iSig < (int)pySigs.size(); ++iSig) { - if (iSig == force) modes.push_back(0); - else modes.push_back(pythiaPtr->rndm.flat() > bfs[iSig]); - if (modes.back() == 0) ++n; - } - if (pythiaPtr->rndm.flat() <= 1.0 / n) break; - if (iTry == NTRYDECAY) { - wgt = 2.; - cout << "Warning in EvtGenDecays::decay: maximum " - << "number of signal decay attempts exceeded"; - } - } - - // Decay the signal particles and mark forced decay. - for (int iSig = 0; iSig < (int)pySigs.size(); ++iSig) { - EvtParticle *egSig = egSigs[iSig]; - Particle *pySig = &event[pySigs[iSig]]; - signal = signals.find(pySig->id()); - if (egSig->getNDaug() > 0) egSig = egSig->getDaug(0); - if (modes[iSig] == 0) egSig->setId(signal->second.egId); - else { - signal->second.modes.getDecayModel(egSig); - egSig->setChannel(signal->second.map[egSig->getChannel()]); - } - if (iSig == force) - pySig->status(pySig->status() == 92 || pySig->status() == 94 ? 96 : 95); - evtgen->generateDecay(egSig); - updateEvent(pySig, egSigs[iSig]); - } - - // Delete all EvtGen particles and return weight. - for (int iPrt = 0; iPrt < (int)egPrts.size(); ++iPrt) - egPrts[iPrt]->deleteTree(); - return 1. - wgt; - -} - -//-------------------------------------------------------------------------- - -// Update the Pythia particle database from EvtGen. - -// Note that only the particle spin type, charge type, nominal mass, -// width, minimum mass, maximum mass, and nominal lifetime are -// set. The name string is not set. - -void EvtGenDecays::updatePythia() { - if (!pythiaPtr || !evtgen) return; - for (int entry = 0; entry < (int)EvtPDL::entries(); ++entry) { - EvtId egId = EvtPDL::getEntry(entry); - int pyId = EvtPDL::getStdHep(egId); - pythiaPtr->particleData.spinType (pyId, EvtPDL::getSpinType(egId)); - pythiaPtr->particleData.chargeType(pyId, EvtPDL::chg3(egId)); - pythiaPtr->particleData.m0 (pyId, EvtPDL::getMass(egId)); - pythiaPtr->particleData.mWidth (pyId, EvtPDL::getWidth(egId)); - pythiaPtr->particleData.mMin (pyId, EvtPDL::getMinMass(egId)); - pythiaPtr->particleData.mMax (pyId, EvtPDL::getMaxMass(egId)); - pythiaPtr->particleData.tau0 (pyId, EvtPDL::getctau(egId)); - } -} - -//-------------------------------------------------------------------------- - -// Update the EvtGen particle database from Pythia. - -// The particle update database between Pythia and EvtGen is not -// symmetric. Particularly, it is not possible to update the spin -// type, charge, or nominal lifetime in the EvtGen particle database. - -void EvtGenDecays::updateEvtGen() { - if (!pythiaPtr || !evtgen) return; - int pyId = pythiaPtr->particleData.nextId(1); - while (pyId != 0) { - EvtId egId = EvtPDL::evtIdFromStdHep(pyId); - EvtPDL::reSetMass (egId, pythiaPtr->particleData.m0(pyId)); - EvtPDL::reSetWidth (egId, pythiaPtr->particleData.mWidth(pyId)); - EvtPDL::reSetMassMin(egId, pythiaPtr->particleData.mMin(pyId)); - EvtPDL::reSetMassMax(egId, pythiaPtr->particleData.mMax(pyId)); - pyId = pythiaPtr->particleData.nextId(pyId); - } -} - -//-------------------------------------------------------------------------- - -// Update the particles to decay with EvtGen, and the selected signals. - -// If final is false, then only signals are initialized in the signal -// map. Any particle or alias that ends with signalSuffix is taken as -// a signal particle. If final is true all particle entries in EvtGen -// are checked to see if they should be set stable in Pythia. If an -// EvtGen particle has no decay modes, then Pythia is still allowed to -// decay the particle. Additionally, the signal decay channels are -// turned off for the non-aliased signal particle. - -void EvtGenDecays::updateData(bool final) { - - // Loop over the EvtGen entries. - if (!pythiaPtr) return; - EvtDecayTable &egTable = EvtDecayTable::getInstance(); - for (int iEntry = 0; iEntry < (int)EvtPDL::entries(); ++iEntry) { - EvtId egId = EvtPDL::getEntry(iEntry); - int pyId = EvtPDL::getStdHep(egId); - if (egTable.getNModes(egId) == 0) continue; - if (excIds.find(pyId) != excIds.end()) continue; - - // Stop Pythia from decaying the particle and include in decay set. - if (final) { - incIds.insert(pyId); - pythiaPtr->particleData.mayDecay(pyId, false); - } - - // Check for signal. - string egName = EvtPDL::name(egId); - if (egName.size() <= signalSuffix.size() || egName.substr - (egName.size() - signalSuffix.size()) != signalSuffix) continue; - signal = signals.find(pyId); - if (signal == signals.end()) { - signals[pyId].status = -1; - signal = signals.find(pyId); - } - signal->second.egId = egId; - signal->second.bfs = vector(2, 0); - if (!final) continue; - - // Get the signal and background decay modes. - vector egList = egTable.getDecayTable(); - int sigIdx = egId.getAlias(); - int bkgIdx = EvtPDL::evtIdFromStdHep(pyId).getAlias(); - if (sigIdx > (int)egList.size() || bkgIdx > (int)egList.size()) continue; - EvtParticleDecayList sigModes = egList[sigIdx]; - EvtParticleDecayList bkgModes = egList[bkgIdx]; - EvtParticleDecayList allModes = egList[bkgIdx]; - double sum(0); - - // Sum signal branching fractions. - for (int iMode = 0; iMode < sigModes.getNMode(); ++iMode) { - EvtDecayBase *mode = sigModes.getDecayModel(iMode); - if (!mode) continue; - signal->second.bfs[0] += mode->getBranchingFraction(); - sum += mode->getBranchingFraction(); - } - - // Sum remaining background branching fractions. - for (int iMode = 0; iMode < allModes.getNMode(); ++iMode) { - EvtDecayBase *mode = allModes.getDecayModel(iMode); - if (!mode) continue; - bool match(false); - for (int jMode = 0; jMode < sigModes.getNMode(); ++jMode) - if (mode->matchingDecay(*(sigModes.getDecayModel(jMode)))) { - match = true; break;} - if (match) bkgModes.removeMode(mode); - else { - signal->second.map.push_back(iMode); - signal->second.bfs[1] += mode->getBranchingFraction(); - sum += mode->getBranchingFraction(); - } - } - signal->second.modes = bkgModes; - for (int iBf = 0; iBf < (int)signal->second.bfs.size(); ++iBf) - signal->second.bfs[iBf] /= sum; - } - if (final) updated = true; - -} - -//-------------------------------------------------------------------------- - -// Update the Pythia event record with an EvtGen decay tree. - -// The production vertex of each particle (which can also be obtained -// in EvtGen via EvtParticle::get4Pos()) is set by the decay vertex of -// its mother, which in turn is calculated from the mother's -// lifetime. The status code 93 is used to indicate an external decay, -// while the status code 94 is used to indicate an oscillated external -// decay. If the progenitor has a single daughter with the same ID, -// this daughter is used as the progenitor. This is used to prevent -// double oscillations. - -// If the arguments after egPro are no NULL and a particle in the -// decay tree is a signal particle, the decay for this particle is -// removed and the particle is stored as a signal candidate in the -// pySigs and egSigs vectors, to be decayed later. However, if any of -// these arguments is NULL then the entire tree is written. - -void EvtGenDecays::updateEvent(Particle *pyPro, EvtParticle *egPro, - vector *pySigs, vector *egSigs, - vector *bfs, double *wgt) { - - // Set up the mother vector. - if (!pythiaPtr) return; - EvtParticle* egMom = egPro; - if (egPro->getNDaug() == 1 && egPro->getPDGId() == - egPro->getDaug(0)->getPDGId()) egMom = egPro->getDaug(0); - Event &event = pythiaPtr->event; - vector< pair > - moms(1, pair(egMom, pyPro->index())); - - // Loop over the mothers. - while (moms.size() != 0) { - - // Check if particle can decay. - egMom = moms.back().first; - int iMom = moms.back().second; - Particle* pyMom = &event[iMom]; - moms.pop_back(); - if (!checkVertex(pyMom)) continue; - bool osc(checkOsc(egMom)); - - // Set the children of the mother. - pyMom->daughters(event.size(), event.size() + egMom->getNDaug() - 1); - pyMom->statusNeg(); - Vec4 vProd = pyMom->vDec(); - for (int iDtr = 0 ; iDtr < (int)egMom->getNDaug(); ++iDtr) { - EvtParticle *egDtr = egMom->getDaug(iDtr); - int id = egDtr->getPDGId(); - EvtVector4R p = egDtr->getP4Lab(); - int idx = event.append(id, 93, iMom, 0, 0, 0, 0, 0, p.get(1), - p.get(2), p.get(3), p.get(0), egDtr->mass()); - Particle *pyDtr = &event.back(); - pyDtr->vProd(vProd); - pyDtr->tau(egDtr->getLifetime()); - if (pySigs && egSigs && bfs && wgt && checkSignal(pyDtr)) { - pySigs->push_back(pyDtr->index()); - egSigs->push_back(egDtr); - bfs->push_back(signal->second.bfs[0]); - (*wgt) *= 1 - bfs->back(); - egDtr->deleteDaughters(); - } - if (osc) pyDtr->status(94); - if (egDtr->getNDaug() > 0) - moms.push_back(pair(egDtr, idx)); - } - } - -} - -//-------------------------------------------------------------------------- - -// Check if a particle can decay. - -// Modified slightly from ParticleDecays::checkVertex. - -bool EvtGenDecays::checkVertex(Particle *pyPro) { - if (!limitDecay) return true; - if (limitTau0 && pyPro->tau0() > tau0Max) return false; - if (limitTau && pyPro->tau() > tauMax) return false; - if (limitRadius && pow2(pyPro->xDec()) + pow2(pyPro->yDec()) - + pow2(pyPro->zDec()) > pow2(rMax)) return false; - if (limitCylinder && (pow2(pyPro->xDec()) + pow2(pyPro->yDec()) - > pow2(xyMax) || abs(pyPro->zDec()) > zMax) ) return false; - return true; -} - -//-------------------------------------------------------------------------- - -// Check if a particle is signal. - -bool EvtGenDecays::checkSignal(Particle *pyPro) { - signal = signals.find(pyPro->id()); - if (signal != signals.end() && (signal->second.status < 0 || - signal->second.status == pyPro->status())) return true; - else return false; -} - -//-------------------------------------------------------------------------- - -// Check if an EvtGen particle has oscillated. - -// The criteria defined here for oscillation is a single daughter but -// with a different ID from the mother. - -bool EvtGenDecays::checkOsc(EvtParticle *egPro) { - if (egPro && egPro->getNDaug() == 1 && - egPro->getPDGId() != egPro->getDaug(0)->getPDGId()) return true; - else return false; -} - -//========================================================================== - -} // end namespace Pythia8 - -#endif // CMSSW_Pythia8_EvtGen_Patched_H diff --git a/GeneratorInterface/Pythia8Interface/plugins/HMC3Py8EGun.cc b/GeneratorInterface/Pythia8Interface/plugins/HMC3Py8EGun.cc index c2a284cd68f6a..1863b84327e9b 100644 --- a/GeneratorInterface/Pythia8Interface/plugins/HMC3Py8EGun.cc +++ b/GeneratorInterface/Pythia8Interface/plugins/HMC3Py8EGun.cc @@ -8,7 +8,7 @@ // EvtGen plugin // -#include "GeneratorInterface/Pythia8Interface/interface/Pythia8EvtGenPatched.h" +#include "Pythia8Plugins/EvtGen.h" namespace gen { diff --git a/GeneratorInterface/Pythia8Interface/plugins/Pythia8Hadronizer.cc b/GeneratorInterface/Pythia8Interface/plugins/Pythia8Hadronizer.cc index 6494f327cdbf8..4dfad6481224c 100644 --- a/GeneratorInterface/Pythia8Interface/plugins/Pythia8Hadronizer.cc +++ b/GeneratorInterface/Pythia8Interface/plugins/Pythia8Hadronizer.cc @@ -48,7 +48,7 @@ using namespace Pythia8; // EvtGen plugin // -#include "GeneratorInterface/Pythia8Interface/interface/Pythia8EvtGenPatched.h" +#include "Pythia8Plugins/EvtGen.h" #include "FWCore/AbstractServices/interface/RandomNumberGenerator.h" #include "FWCore/Concurrency/interface/SharedResourceNames.h" diff --git a/GeneratorInterface/Pythia8Interface/plugins/Pythia8HepMC3Hadronizer.cc b/GeneratorInterface/Pythia8Interface/plugins/Pythia8HepMC3Hadronizer.cc index 607ff0889aa55..ab5683d735dae 100644 --- a/GeneratorInterface/Pythia8Interface/plugins/Pythia8HepMC3Hadronizer.cc +++ b/GeneratorInterface/Pythia8Interface/plugins/Pythia8HepMC3Hadronizer.cc @@ -47,7 +47,7 @@ using namespace Pythia8; // EvtGen plugin // -#include "GeneratorInterface/Pythia8Interface/interface/Pythia8EvtGenPatched.h" +#include "Pythia8Plugins/EvtGen.h" #include "FWCore/AbstractServices/interface/RandomNumberGenerator.h" #include "FWCore/Concurrency/interface/SharedResourceNames.h" diff --git a/GeneratorInterface/Pythia8Interface/src/Py8GunBase.cc b/GeneratorInterface/Pythia8Interface/src/Py8GunBase.cc index 2957c936e2fee..b560d45205be8 100644 --- a/GeneratorInterface/Pythia8Interface/src/Py8GunBase.cc +++ b/GeneratorInterface/Pythia8Interface/src/Py8GunBase.cc @@ -6,7 +6,7 @@ // EvtGen plugin // -#include "GeneratorInterface/Pythia8Interface/interface/Pythia8EvtGenPatched.h" +#include "Pythia8Plugins/EvtGen.h" using namespace Pythia8; From 4dcf499f1f2c13ad056ee33b8ba5c82257e57c1d Mon Sep 17 00:00:00 2001 From: Ka Wa Ho Date: Fri, 22 May 2026 16:55:08 -0500 Subject: [PATCH 5/8] bump to EvtGen v3.0.0; enable fsr options --- .../plugins/EvtGen/EvtGenInterface.cc | 25 ++++++++++++++++--- 1 file changed, 22 insertions(+), 3 deletions(-) diff --git a/GeneratorInterface/EvtGenInterface/plugins/EvtGen/EvtGenInterface.cc b/GeneratorInterface/EvtGenInterface/plugins/EvtGen/EvtGenInterface.cc index 227ca0e13ff84..3212e9d245f73 100644 --- a/GeneratorInterface/EvtGenInterface/plugins/EvtGen/EvtGenInterface.cc +++ b/GeneratorInterface/EvtGenInterface/plugins/EvtGen/EvtGenInterface.cc @@ -301,7 +301,15 @@ void EvtGenInterface::ensureEvtGenOnThread() { bool usePythia = fPSet->getUntrackedParameter("use_internal_pythia", true); bool useTauola = fPSet->getUntrackedParameter("use_internal_tauola", true); - bool usePhotos = fPSet->getUntrackedParameter("use_internal_photos", true); + + if (fPSet->existsAs("use_internal_fsr", false)) { + edm::LogWarning("EvtGenInterface::ensureEvtGenOnThread") + << "The 'use_internal_fsr' parameter is deprecated and has no effect. " + "Please use 'fsr_model' instead (set to 'none' to disable FSR, or " + "'photos', 'sherpa', or 'vincia' to select the FSR model)."; + } + + std::string fsrModel = fPSet->getParameter("fsr_model"); //Setup evtGen following instructions on http://evtgen.warwick.ac.uk/docs/external/ bool convertPythiaCodes = fPSet->getUntrackedParameter( @@ -322,8 +330,19 @@ void EvtGenInterface::ensureEvtGenOnThread() { // Set up the default external generator list: Photos, Pythia and/or Tauola EvtExternalGenList genList(convertPythiaCodes, pythiaDir, photonType, useEvtGenRandom); EvtAbsRadCorr* radCorrEngine = nullptr; - if (usePhotos) - radCorrEngine = genList.getPhotosModel(); // Get interface to radiative correction engine + if (fsrModel == "none") { + // FSR disabled + } else if (fsrModel == "photos") { + radCorrEngine = genList.getPhotosModel(); + } else if (fsrModel == "sherpa") { + radCorrEngine = genList.getSherpaPhotonsModel(1e-7, 1, 0); + } else if (fsrModel == "vincia") { + radCorrEngine = genList.getVinciaQEDModel(1.0e-7); + } else { + throw cms::Exception("Configuration") + << "EvtGenInterface: unknown fsr_model '" << fsrModel + << "'. Allowed values are 'none', 'photos', 'sherpa', or 'vincia'."; + } std::list extraModels = genList.getListOfModels(); // get interface to Pythia and Tauola std::list myExtraModels; for (unsigned int i = 0; i < extraModels.size(); i++) { From d66589b4fd1853d679b2632d573331ca6bcdb109 Mon Sep 17 00:00:00 2001 From: Ka Wa Ho Date: Tue, 26 May 2026 23:18:52 -0500 Subject: [PATCH 6/8] add fsr models to evtgen tests --- .../plugins/test/Py8_lambadb_evtgen1_Lb2plnuLCSR_cfg.py | 1 + .../EvtGenInterface/test/Py8_bplus_evtgen1_cfg.py | 3 ++- GeneratorInterface/EvtGenInterface/test/Py8_tt_evtgen1_cfg.py | 3 ++- .../EvtGenInterface/test/external_Py8_bplus_evtgen1_cfg.py | 3 ++- 4 files changed, 7 insertions(+), 3 deletions(-) diff --git a/GeneratorInterface/EvtGenInterface/plugins/test/Py8_lambadb_evtgen1_Lb2plnuLCSR_cfg.py b/GeneratorInterface/EvtGenInterface/plugins/test/Py8_lambadb_evtgen1_Lb2plnuLCSR_cfg.py index d9678b66c8d3b..c845c0d3fd33f 100644 --- a/GeneratorInterface/EvtGenInterface/plugins/test/Py8_lambadb_evtgen1_Lb2plnuLCSR_cfg.py +++ b/GeneratorInterface/EvtGenInterface/plugins/test/Py8_lambadb_evtgen1_Lb2plnuLCSR_cfg.py @@ -56,6 +56,7 @@ decay_table = cms.string('GeneratorInterface/EvtGenInterface/data/DECAY_NOLONGLIFE.DEC'), list_forced_decays = cms.vstring('MyLambda_b0','Myanti-Lambda_b0'), particle_property_file = cms.FileInPath('GeneratorInterface/EvtGenInterface/data/evt.pdl'), + fsr_model = cms.string('sherpa'), operates_on_particles = cms.vint32(5122), #only care about signal particle user_decay_embedded = cms.vstring( """ diff --git a/GeneratorInterface/EvtGenInterface/test/Py8_bplus_evtgen1_cfg.py b/GeneratorInterface/EvtGenInterface/test/Py8_bplus_evtgen1_cfg.py index 08dc631309246..289c2d96b9e85 100644 --- a/GeneratorInterface/EvtGenInterface/test/Py8_bplus_evtgen1_cfg.py +++ b/GeneratorInterface/EvtGenInterface/test/Py8_bplus_evtgen1_cfg.py @@ -97,7 +97,8 @@ '0.0168 Mychi_c1 K+ SVS ;', 'Enddecay', 'CDecay MyB-', - 'End') + 'End'), + fsr_model = cms.string('sherpa'), ), parameterSets = cms.vstring('EvtGen1') ), diff --git a/GeneratorInterface/EvtGenInterface/test/Py8_tt_evtgen1_cfg.py b/GeneratorInterface/EvtGenInterface/test/Py8_tt_evtgen1_cfg.py index db0b7f6c948cf..9c2a3e81ae620 100644 --- a/GeneratorInterface/EvtGenInterface/test/Py8_tt_evtgen1_cfg.py +++ b/GeneratorInterface/EvtGenInterface/test/Py8_tt_evtgen1_cfg.py @@ -24,7 +24,8 @@ particle_property_file = cms.FileInPath('GeneratorInterface/EvtGenInterface/data/evt.pdl'), convertPythiaCodes = cms.untracked.bool(False), # this is needed since 1.6 list_forced_decays = cms.vstring(), - operates_on_particles = cms.vint32(0) #will decay all particles coded in interface, it test the whole system + operates_on_particles = cms.vint32(0), #will decay all particles coded in interface, it test the whole system, + fsr_model = cms.string('sherpa') ), parameterSets = cms.vstring('EvtGen1') ), diff --git a/GeneratorInterface/EvtGenInterface/test/external_Py8_bplus_evtgen1_cfg.py b/GeneratorInterface/EvtGenInterface/test/external_Py8_bplus_evtgen1_cfg.py index ceb2a1d1440f9..72d8c8457a2fd 100644 --- a/GeneratorInterface/EvtGenInterface/test/external_Py8_bplus_evtgen1_cfg.py +++ b/GeneratorInterface/EvtGenInterface/test/external_Py8_bplus_evtgen1_cfg.py @@ -97,7 +97,8 @@ '0.0168 Mychi_c1 K+ SVS ;', 'Enddecay', 'CDecay MyB-', - 'End') + 'End'), + fsr_model = cms.string('sherpa') ), parameterSets = cms.vstring('EvtGen1') ), From 150b0c721f6ade0baac5d0f2fc281ef3e6d68824 Mon Sep 17 00:00:00 2001 From: Ka Wa Ho Date: Thu, 28 May 2026 11:58:51 -0500 Subject: [PATCH 7/8] ThreadHandoff.h is moved to FWCore/Utilities --- GeneratorInterface/ExternalDecays/BuildFile.xml | 2 +- .../ExternalDecays/interface/EvtGenThreadOwner.h | 5 ++--- 2 files changed, 3 insertions(+), 4 deletions(-) diff --git a/GeneratorInterface/ExternalDecays/BuildFile.xml b/GeneratorInterface/ExternalDecays/BuildFile.xml index 887efbc72d59e..fbe5ec16d2f99 100644 --- a/GeneratorInterface/ExternalDecays/BuildFile.xml +++ b/GeneratorInterface/ExternalDecays/BuildFile.xml @@ -1,7 +1,7 @@ - + diff --git a/GeneratorInterface/ExternalDecays/interface/EvtGenThreadOwner.h b/GeneratorInterface/ExternalDecays/interface/EvtGenThreadOwner.h index bb0baee995e1d..d3c2573efbedd 100644 --- a/GeneratorInterface/ExternalDecays/interface/EvtGenThreadOwner.h +++ b/GeneratorInterface/ExternalDecays/interface/EvtGenThreadOwner.h @@ -2,9 +2,8 @@ #define gen_EvtGenThreadOwner_h // Pins one dedicated pthread per stream and marshals every EvtGen call onto it. -// Mirrors SimG4Core's OscarMTProducer + omt::ThreadHandoff pattern -#include "SimG4Core/Application/interface/ThreadHandoff.h" +#include "FWCore/Utilities/interface/ThreadHandoff.h" #include "FWCore/ServiceRegistry/interface/ServiceRegistry.h" namespace gen { @@ -26,7 +25,7 @@ namespace gen { } private: - omt::ThreadHandoff m_handoff; + edm::ThreadHandoff m_handoff; }; } From ba5de48be37734de8e11a33d134cddae6fa43813 Mon Sep 17 00:00:00 2001 From: Ka Wa Ho Date: Fri, 29 May 2026 16:21:53 -0500 Subject: [PATCH 8/8] add dummy MPI initialize --- GeneratorInterface/EvtGenInterface/BuildFile.xml | 1 + .../plugins/EvtGen/EvtGenInterface.cc | 14 ++++++++++++-- 2 files changed, 13 insertions(+), 2 deletions(-) diff --git a/GeneratorInterface/EvtGenInterface/BuildFile.xml b/GeneratorInterface/EvtGenInterface/BuildFile.xml index 1c760778eb73e..96c60fae69557 100644 --- a/GeneratorInterface/EvtGenInterface/BuildFile.xml +++ b/GeneratorInterface/EvtGenInterface/BuildFile.xml @@ -4,6 +4,7 @@ + diff --git a/GeneratorInterface/EvtGenInterface/plugins/EvtGen/EvtGenInterface.cc b/GeneratorInterface/EvtGenInterface/plugins/EvtGen/EvtGenInterface.cc index 3212e9d245f73..ed8ff91ee0430 100644 --- a/GeneratorInterface/EvtGenInterface/plugins/EvtGen/EvtGenInterface.cc +++ b/GeneratorInterface/EvtGenInterface/plugins/EvtGen/EvtGenInterface.cc @@ -44,6 +44,8 @@ #include "HepPID/ParticleIDTranslations.hh" +#include "ATOOLS/Org/My_MPI.H" + #include "TString.h" #include #include @@ -54,6 +56,14 @@ thread_local std::unique_ptr gen::EvtGenInterface::m_EvtGen{}; using namespace gen; using namespace edm; +// Sherpa is built with MPI and aborts unless MPI is initialized. +// We only need a dummy init, and it must happen exactly once per process not per thread. +void initDummyMPIOnce() { + [[maybe_unused]] static const bool init = [] { + MPI_Init(0, nullptr); + return true; + }(); +} EvtGenInterface::EvtGenInterface(const ParameterSet& pset) { fPSet = std::make_unique(pset); @@ -287,8 +297,7 @@ void EvtGenInterface::SetDefault_m_PDGs() { } } -// Destructor is a no-op: m_EvtGen is a thread_local std::unique_ptr that cleans up at -// thread exit. +// Destructor is a no-op: m_EvtGen is a thread_local std::unique_ptr that cleans up at thread exit. EvtGenInterface::~EvtGenInterface() {} void EvtGenInterface::ensureEvtGenOnThread() { @@ -335,6 +344,7 @@ void EvtGenInterface::ensureEvtGenOnThread() { } else if (fsrModel == "photos") { radCorrEngine = genList.getPhotosModel(); } else if (fsrModel == "sherpa") { + initDummyMPIOnce(); // Sherpa is built with MPI and complains if not initialized; do this once per process. radCorrEngine = genList.getSherpaPhotonsModel(1e-7, 1, 0); } else if (fsrModel == "vincia") { radCorrEngine = genList.getVinciaQEDModel(1.0e-7);