upgrade EvtGen interface for v03.00.00

GeneratorInterface/EvtGenInterface/BuildFile.xml

@@ -4,6 +4,7 @@
 <use name="hepmc"/>
 <use name="clhep"/>
 <use name="heppdt"/>
+<use name="sherpa"/>
 <use name="evtgen"/>
 <export>
   <lib name="1"/>

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

GeneratorInterface/EvtGenInterface/plugins/EvtGen/EvtGenInterface.cc

@@ -36,25 +36,38 @@
 #include "EvtGenBase/EvtDecayBase.hh"
 #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"
 
 #include "HepPID/ParticleIDTranslations.hh"
 
+#include "ATOOLS/Org/My_MPI.H"
+
 #include "TString.h"
 #include <string>
 #include <cstdlib>
 #include <cstdio>
 
+thread_local std::unique_ptr<EvtGen> gen::EvtGenInterface::m_EvtGen{};
+
 using namespace gen;
 using namespace edm;
 
-CLHEP::HepRandomEngine* EvtGenInterface::fRandomEngine;
+// 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 = new ParameterSet(pset);
-  the_engine = new myEvtRandomEngine(nullptr);
+  fPSet = std::make_unique<ParameterSet>(pset);
+  the_engine = std::make_unique<myEvtRandomEngine>(nullptr);
 }
 
 void EvtGenInterface::SetDefault_m_PDGs() {
@@ -284,19 +297,28 @@ 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<std::string>("decay_table"));
   edm::FileInPath pdt(fPSet->getParameter<edm::FileInPath>("particle_property_file"));
 
   bool usePythia = fPSet->getUntrackedParameter<bool>("use_internal_pythia", true);
   bool useTauola = fPSet->getUntrackedParameter<bool>("use_internal_tauola", true);
-  bool usePhotos = fPSet->getUntrackedParameter<bool>("use_internal_photos", true);
+
+  if (fPSet->existsAs<bool>("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<std::string>("fsr_model");
 
   //Setup evtGen following instructions on http://evtgen.warwick.ac.uk/docs/external/
   bool convertPythiaCodes = fPSet->getUntrackedParameter<bool>(
@@ -305,8 +327,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);
   }
 
@@ -317,8 +339,20 @@ void EvtGenInterface::init() {
   // 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") {
+    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);
+  } else {
+    throw cms::Exception("Configuration")
+        << "EvtGenInterface: unknown fsr_model '" << fsrModel
+        << "'. Allowed values are 'none', 'photos', 'sherpa', or 'vincia'.";
+  }
   std::list<EvtDecayBase*> extraModels = genList.getListOfModels();  // get interface to Pythia and Tauola
   std::list<EvtDecayBase*> myExtraModels;
   for (unsigned int i = 0; i < extraModels.size(); i++) {
@@ -339,21 +373,12 @@ void EvtGenInterface::init() {
     std::list<EvtDecayBase*>::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<int>("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<EvtGen>(
+      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")) {
@@ -370,9 +395,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++) {
@@ -382,6 +406,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<int>("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")) {
@@ -452,6 +493,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"
@@ -496,8 +541,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 +661,11 @@ void EvtGenInterface::update_particles(HepMC::GenParticle* partHep, HepMC::GenEv
 
 void EvtGenInterface::setRandomEngine(CLHEP::HepRandomEngine* v) {
   the_engine->setRandomEngine(v);
-  fRandomEngine = v;
+  EvtRandom::setRandomEngine(the_engine.get());
 }
 
 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) {

GeneratorInterface/EvtGenInterface/plugins/EvtGen/EvtGenInterface.h

@@ -43,7 +43,7 @@ namespace gen {
     const std::vector<int>& 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);
@@ -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<EvtGen> m_EvtGen;
 
     std::vector<EvtId> forced_id;    // EvtGen Id's of particles  which are to be forced by EvtGen
     std::vector<int> forced_pdgids;  // PDG Id's of particles which are to be forced by EvtGen
@@ -65,10 +67,9 @@ namespace gen {
     std::vector<double> polarize_pol;
     std::map<int, float> polarizations;
     int BmixingOption = 1;
-    edm::ParameterSet* fPSet;
+    std::unique_ptr<edm::ParameterSet> fPSet;
 
-    static CLHEP::HepRandomEngine* fRandomEngine;
-    myEvtRandomEngine* the_engine;
+    std::unique_ptr<myEvtRandomEngine> the_engine;
   };
 }  // namespace gen
 #endif

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"

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(
 """

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')
     ),

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')
     ),

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')
     ),

GeneratorInterface/ExternalDecays/BuildFile.xml

@@ -1,6 +1,7 @@
 <use name="FWCore/Concurrency"/>
 <use name="FWCore/ParameterSet"/>
 <use name="FWCore/Framework"/>
+<use name="FWCore/Utilities"/>
 <use name="heppdt"/>
 <use name="clhep"/>
 <use name="hepmc"/>

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<EvtGenThreadOwner> fThreadOwner;
+    std::unique_ptr<EvtGenInterfaceBase> fEvtGenInterface;
     //std::unique_ptr<TauolaInterfaceBase> fTauolaInterface;
-    //std::unique_ptr<EvtGenInterfaceBase> fEvtGenInterface;
     //std::unique_ptr<PhotosInterfaceBase> fPhotosInterface;
     std::vector<int> fPDGs;
     std::vector<std::string> fSpecialSettings;

GeneratorInterface/ExternalDecays/interface/EvtGenThreadOwner.h

@@ -0,0 +1,33 @@
+#ifndef gen_EvtGenThreadOwner_h
+#define gen_EvtGenThreadOwner_h
+
+// Pins one dedicated pthread per stream and marshals every EvtGen call onto it.
+
+#include "FWCore/Utilities/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 <typename F>
+    void run(F&& f) {
+      auto token = edm::ServiceRegistry::instance().presentToken();
+      m_handoff.runAndWait([token, &f]() {
+        edm::ServiceRegistry::Operate guard{token};
+        f();
+      });
+    }
+
+  private:
+    edm::ThreadHandoff m_handoff;
+  };
+
+} 
+
+#endif