From 69e7cae0046d4b9636cb0931a4cc39a8f5acde8f Mon Sep 17 00:00:00 2001
From: Laurent Forthomme <laurent.forthomme@cern.ch>
Date: Mon, 25 Apr 2022 10:32:43 +0200
Subject: [PATCH] Interface to the CepGen event generator

Along with the "generator filter" implementation linking everything together, this contains
generator fragments for a LPAIR and a PPtoFF single-dissociative generation.
The PPtoFF fragment also demonstrates the capability to use Pythia 6 for the fragmentation
of the dissociative beam remnants.

Requires CepGen v1.2.0 at least.
---
 ...n_LPAIR_GamGamToLL_singlediss_13TeV_cfi.py | 25 +++++
 ..._pptoff_GamGamToLL_singlediss_13TeV_cfi.py | 28 ++++++
 .../CepGenInterface/BuildFile.xml             | 10 ++
 .../interface/CepGenEventGenerator.h          | 43 +++++++++
 .../interface/CepGenGeneratorFilter.h         | 24 +++++
 .../interface/CepGenParametersConverter.h     | 41 +++++++++
 .../CepGenInterface/plugins/BuildFile.xml     |  7 ++
 .../plugins/CepGenGeneratorFilter.cc          |  8 ++
 .../python/cepgenDefaultParameters_cff.py     | 23 +++++
 .../src/CepGenEventGenerator.cc               | 92 +++++++++++++++++++
 .../src/CepGenGeneratorFilter.cc              |  4 +
 11 files changed, 305 insertions(+)
 create mode 100644 Configuration/Generator/python/CepGen_LPAIR_GamGamToLL_singlediss_13TeV_cfi.py
 create mode 100644 Configuration/Generator/python/CepGen_pptoff_GamGamToLL_singlediss_13TeV_cfi.py
 create mode 100644 GeneratorInterface/CepGenInterface/BuildFile.xml
 create mode 100644 GeneratorInterface/CepGenInterface/interface/CepGenEventGenerator.h
 create mode 100644 GeneratorInterface/CepGenInterface/interface/CepGenGeneratorFilter.h
 create mode 100644 GeneratorInterface/CepGenInterface/interface/CepGenParametersConverter.h
 create mode 100644 GeneratorInterface/CepGenInterface/plugins/BuildFile.xml
 create mode 100644 GeneratorInterface/CepGenInterface/plugins/CepGenGeneratorFilter.cc
 create mode 100644 GeneratorInterface/CepGenInterface/python/cepgenDefaultParameters_cff.py
 create mode 100644 GeneratorInterface/CepGenInterface/src/CepGenEventGenerator.cc
 create mode 100644 GeneratorInterface/CepGenInterface/src/CepGenGeneratorFilter.cc

diff --git a/Configuration/Generator/python/CepGen_LPAIR_GamGamToLL_singlediss_13TeV_cfi.py b/Configuration/Generator/python/CepGen_LPAIR_GamGamToLL_singlediss_13TeV_cfi.py
new file mode 100644
index 0000000000000..f089fb29c5c3e
--- /dev/null
+++ b/Configuration/Generator/python/CepGen_LPAIR_GamGamToLL_singlediss_13TeV_cfi.py
@@ -0,0 +1,25 @@
+import FWCore.ParameterSet.Config as cms
+
+from GeneratorInterface.CepGenInterface.cepgenDefaultParameters_cff import *
+
+generator = cms.EDFilter("CepGenGeneratorFilter",
+    process = cms.PSet(
+        name = cms.string('lpair'),
+        kinematics = cms.PSet(
+            mode = cms.uint32(2),  # 1 = elastic, 2-3 = SD, 4 = DD
+            structureFunctions = cms.PSet(
+                name = cms.int32(301)  # see https://cepgen.hepforge.org/raw-modules
+            ),
+            cmEnergy = cms.double(13000.0),
+            eta = cms.vdouble(-2.5, 2.5),
+            pt = cms.vdouble(25., 9999.)
+        ),
+        processParameters = cms.PSet(
+            pair = cms.uint32(13)
+        )
+    ),
+    modifierModules = cepgenPythia6BeamFragmenter,
+    outputModules = cepgenOutputModules,
+    maxEventsToPrint = cms.untracked.int32(0),
+    verbosity = cms.untracked.int32(0)
+)
diff --git a/Configuration/Generator/python/CepGen_pptoff_GamGamToLL_singlediss_13TeV_cfi.py b/Configuration/Generator/python/CepGen_pptoff_GamGamToLL_singlediss_13TeV_cfi.py
new file mode 100644
index 0000000000000..d9d7c5ed1f51e
--- /dev/null
+++ b/Configuration/Generator/python/CepGen_pptoff_GamGamToLL_singlediss_13TeV_cfi.py
@@ -0,0 +1,28 @@
+import FWCore.ParameterSet.Config as cms
+
+from GeneratorInterface.CepGenInterface.cepgenDefaultParameters_cff import *
+
+generator = cms.EDFilter("CepGenGeneratorFilter",
+    process = cms.PSet(
+        name = cms.string('pptoff'),
+        kinematics = cms.PSet(
+            mode = cms.uint32(2),  # 1 = elastic, 2-3 = SD, 4 = DD
+            structureFunctions = cms.PSet(
+                name = cms.int32(301)  # see https://cepgen.hepforge.org/raw-modules
+            ),
+            cmEnergy = cms.double(13000.0),
+            qt = cms.vdouble(0., 50.),
+            rapidity = cms.vdouble(-6., 6.),
+            eta = cms.vdouble(-2.5, 2.5),
+            pt = cms.vdouble(25., 9999.),
+            ptdiff = cms.vdouble(0., 500.)
+        ),
+        processParameters = cms.PSet(
+            ktFactorised = cms.bool(True),
+            pair = cms.uint32(13)
+        )
+    ),
+    outputModules = cepgenOutputModules,
+    maxEventsToPrint = cms.untracked.int32(0),
+    verbosity = cms.untracked.int32(0)
+)
diff --git a/GeneratorInterface/CepGenInterface/BuildFile.xml b/GeneratorInterface/CepGenInterface/BuildFile.xml
new file mode 100644
index 0000000000000..940486bc97ba8
--- /dev/null
+++ b/GeneratorInterface/CepGenInterface/BuildFile.xml
@@ -0,0 +1,10 @@
+<use name="CepGen"/>
+<lib name="CepGenHepMC2"/>
+<use name="FWCore/Framework"/>
+<use name="GeneratorInterface/Core"/>
+<use name="clhep"/>
+<use name="hepmc"/>
+<flags CXXFLAGS="-Wno-error=missing-braces -Wno-missing-braces -Wno-non-virtual-dtor"/>
+<export>
+  <lib name="1"/>
+</export>
diff --git a/GeneratorInterface/CepGenInterface/interface/CepGenEventGenerator.h b/GeneratorInterface/CepGenInterface/interface/CepGenEventGenerator.h
new file mode 100644
index 0000000000000..6ef187abc9602
--- /dev/null
+++ b/GeneratorInterface/CepGenInterface/interface/CepGenEventGenerator.h
@@ -0,0 +1,43 @@
+// CepGen-CMSSW interfacing module
+//   2022-2024, Laurent Forthomme
+
+#ifndef GeneratorInterface_CepGenInterface_CepGenEventGenerator_h
+#define GeneratorInterface_CepGenInterface_CepGenEventGenerator_h
+
+#include "FWCore/Framework/interface/ConsumesCollector.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "GeneratorInterface/Core/interface/BaseHadronizer.h"
+
+#include <CepGen/Generator.h>
+
+namespace gen {
+  class CepGenEventGenerator : public BaseHadronizer {
+  public:
+    explicit CepGenEventGenerator(const edm::ParameterSet&, edm::ConsumesCollector&&);
+    virtual ~CepGenEventGenerator();
+
+    bool readSettings(int) { return true; }
+    bool declareStableParticles(const std::vector<int>&) { return true; }
+    bool declareSpecialSettings(const std::vector<std::string>&) { return true; }
+
+    bool initializeForInternalPartons();
+    bool generatePartonsAndHadronize();
+    bool decay() { return true; }  // NOT used - let's call it "design imperfection"
+    bool residualDecay() { return true; }
+
+    void finalizeEvent() {}
+    void statistics() {}
+
+    const char* classname() const { return "CepGenEventGenerator"; }
+    const std::vector<std::string>& doSharedResources() const override { return shared_resources_; }
+
+  private:
+    cepgen::Generator* gen_{nullptr};
+    const cepgen::ParametersList proc_params_;
+    std::vector<std::pair<std::string, cepgen::ParametersList> > modif_modules_, output_modules_;
+    const std::vector<std::string> shared_resources_;
+    HepMC::GenCrossSection xsec_;
+  };
+}  // namespace gen
+
+#endif
diff --git a/GeneratorInterface/CepGenInterface/interface/CepGenGeneratorFilter.h b/GeneratorInterface/CepGenInterface/interface/CepGenGeneratorFilter.h
new file mode 100644
index 0000000000000..a922d9d73ff34
--- /dev/null
+++ b/GeneratorInterface/CepGenInterface/interface/CepGenGeneratorFilter.h
@@ -0,0 +1,24 @@
+// CepGen-CMSSW interfacing module
+//   2022-2024, Laurent Forthomme
+
+#ifndef GeneratorInterface_CepGenInterface_CepGenGeneratorFilter_h
+#define GeneratorInterface_CepGenInterface_CepGenGeneratorFilter_h
+
+#include "GeneratorInterface/Core/interface/GeneratorFilter.h"
+#include "GeneratorInterface/CepGenInterface/interface/CepGenEventGenerator.h"
+#include "GeneratorInterface/ExternalDecays/interface/ExternalDecayDriver.h"
+
+namespace edm {
+  template <>
+  inline GeneratorFilter<gen::CepGenEventGenerator, gen::ExternalDecayDriver>::GeneratorFilter(
+      const ParameterSet& iConfig)
+      : hadronizer_(iConfig, consumesCollector()) {
+    init(iConfig);
+  }
+}  // namespace edm
+
+namespace gen {
+  using CepGenGeneratorFilter = edm::GeneratorFilter<CepGenEventGenerator, gen::ExternalDecayDriver>;
+}
+
+#endif
diff --git a/GeneratorInterface/CepGenInterface/interface/CepGenParametersConverter.h b/GeneratorInterface/CepGenInterface/interface/CepGenParametersConverter.h
new file mode 100644
index 0000000000000..240f442b59c4e
--- /dev/null
+++ b/GeneratorInterface/CepGenInterface/interface/CepGenParametersConverter.h
@@ -0,0 +1,41 @@
+// CepGen-CMSSW interfacing module
+//   2022-2024, Laurent Forthomme
+
+#ifndef GeneratorInterface_CepGenInterface_CepGenParametersConverter_h
+#define GeneratorInterface_CepGenInterface_CepGenParametersConverter_h
+
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+
+#include <CepGen/Core/ParametersList.h>
+
+namespace cepgen {
+  ParametersList fromParameterSet(const edm::ParameterSet& iConfig) {
+    ParametersList params;
+    for (const auto& param : iConfig.getParameterNames()) {
+      const auto cepgen_param = param == "name" ? MODULE_NAME : param;
+      if (iConfig.existsAs<bool>(param))
+        params.set(cepgen_param, iConfig.getParameter<bool>(param));
+      if (iConfig.existsAs<int>(param))
+        params.set(cepgen_param, iConfig.getParameter<int>(param));
+      if (iConfig.existsAs<unsigned>(param))
+        params.set<unsigned long long>(cepgen_param, iConfig.getParameter<unsigned>(param));
+      if (iConfig.existsAs<double>(param))
+        params.set(cepgen_param, iConfig.getParameter<double>(param));
+      if (iConfig.existsAs<std::string>(param))
+        params.set(cepgen_param, iConfig.getParameter<std::string>(param));
+      if (iConfig.existsAs<std::vector<double> >(param)) {
+        const auto& vec = iConfig.getParameter<std::vector<double> >(param);
+        if (vec.size() == 2)
+          params.set<Limits>(cepgen_param, Limits{vec.at(0), vec.at(1)});
+        params.set(cepgen_param, iConfig.getParameter<std::vector<double> >(param));
+      }
+      if (iConfig.existsAs<std::vector<std::string> >(param))
+        params.set(cepgen_param, iConfig.getParameter<std::vector<std::string> >(param));
+      if (iConfig.existsAs<edm::ParameterSet>(param))
+        params.set(cepgen_param, fromParameterSet(iConfig.getParameter<edm::ParameterSet>(param)));
+    }
+    return params;
+  }
+}  // namespace cepgen
+
+#endif
diff --git a/GeneratorInterface/CepGenInterface/plugins/BuildFile.xml b/GeneratorInterface/CepGenInterface/plugins/BuildFile.xml
new file mode 100644
index 0000000000000..2d00a611c2704
--- /dev/null
+++ b/GeneratorInterface/CepGenInterface/plugins/BuildFile.xml
@@ -0,0 +1,7 @@
+<use name="FWCore/Framework"/>
+<use name="GeneratorInterface/Core"/>
+<use name="GeneratorInterface/CepGenInterface"/>
+<use name="GeneratorInterface/ExternalDecays"/>
+<library file="CepGenGeneratorFilter.cc" name="GeneratorInterfaceCepGenGenerator">
+  <flags EDM_PLUGIN="1"/>
+</library>
diff --git a/GeneratorInterface/CepGenInterface/plugins/CepGenGeneratorFilter.cc b/GeneratorInterface/CepGenInterface/plugins/CepGenGeneratorFilter.cc
new file mode 100644
index 0000000000000..ae3fdc0d8789a
--- /dev/null
+++ b/GeneratorInterface/CepGenInterface/plugins/CepGenGeneratorFilter.cc
@@ -0,0 +1,8 @@
+// CepGen-CMSSW interfacing module
+//   2022-2024, Laurent Forthomme
+
+#include "FWCore/Framework/interface/MakerMacros.h"
+#include "GeneratorInterface/CepGenInterface/interface/CepGenGeneratorFilter.h"
+
+using gen::CepGenGeneratorFilter;
+DEFINE_FWK_MODULE(CepGenGeneratorFilter);
diff --git a/GeneratorInterface/CepGenInterface/python/cepgenDefaultParameters_cff.py b/GeneratorInterface/CepGenInterface/python/cepgenDefaultParameters_cff.py
new file mode 100644
index 0000000000000..22401fac2f296
--- /dev/null
+++ b/GeneratorInterface/CepGenInterface/python/cepgenDefaultParameters_cff.py
@@ -0,0 +1,23 @@
+import FWCore.ParameterSet.Config as cms
+
+cepgenOutputModules = cms.untracked.PSet(
+    # place here all CepGen modules steering for output
+    # e.g. for a printout of every 1000th event:
+    #dump = cms.PSet(printEvery = cms.uint32(1000))
+)
+
+cepgenPythia6BeamFragmenter = cms.untracked.PSet(
+    pythia6 = cms.PSet(
+        maxTrials = cms.int32(10),
+        seed = cms.int32(42),
+        preConfiguration = cms.vstring(
+            'MSTU(21)=1     ! Check on possible errors during program execution',
+            'MSTU(25)=0     ! No warnings are written',
+            'MSTJ(22)=2     ! Decay those unstable particles',
+            'PARJ(71)=10 .  ! for which ctau  10 mm',
+            'MSTP(33)=0     ! no K factors in hard cross sections',
+            'MSTP(2)=1      ! which order running alphaS',
+            'MSTP(81)=0     ! multiple parton interactions 1 is Pythia default'
+        )
+    )
+)
diff --git a/GeneratorInterface/CepGenInterface/src/CepGenEventGenerator.cc b/GeneratorInterface/CepGenInterface/src/CepGenEventGenerator.cc
new file mode 100644
index 0000000000000..4507f5d7ec415
--- /dev/null
+++ b/GeneratorInterface/CepGenInterface/src/CepGenEventGenerator.cc
@@ -0,0 +1,92 @@
+// CepGen-CMSSW interfacing module
+//   2022-2024, Laurent Forthomme
+
+#include "FWCore/MessageLogger/interface/MessageLogger.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "FWCore/Utilities/interface/StreamID.h"
+
+#include "GeneratorInterface/CepGenInterface/interface/CepGenEventGenerator.h"
+#include "GeneratorInterface/CepGenInterface/interface/CepGenParametersConverter.h"
+
+#include <CepGen/Core/Exception.h>
+#include <CepGen/Core/RunParameters.h>
+#include <CepGen/Event/Event.h>
+#include <CepGen/EventFilter/EventExporter.h>
+#include <CepGen/EventFilter/EventModifier.h>
+#include <CepGen/Generator.h>
+#include <CepGen/Modules/EventExporterFactory.h>
+#include <CepGen/Modules/EventModifierFactory.h>
+#include <CepGen/Modules/ProcessFactory.h>
+#include <CepGen/Process/Process.h>
+#include <CepGenAddOns/HepMC2Wrapper/HepMC2EventInterface.h>
+
+using namespace gen;
+
+CepGenEventGenerator::CepGenEventGenerator(const edm::ParameterSet& iConfig, edm::ConsumesCollector&& iC)
+    : BaseHadronizer(iConfig),
+      proc_params_(cepgen::fromParameterSet(iConfig.getParameter<edm::ParameterSet>("process"))) {
+  // specify the overall module verbosity
+  cepgen::utils::Logger::get().setLevel(
+      static_cast<cepgen::utils::Logger::Level>(iConfig.getUntrackedParameter<int>("verbosity", 0)));
+
+  // build the process
+  edm::LogInfo("CepGenEventGenerator") << "Process to be generated: " << proc_params_ << ".";
+
+  const auto modif_mods = cepgen::fromParameterSet(
+      iConfig.getUntrackedParameter<edm::ParameterSet>("modifierModules", edm::ParameterSet{}));
+  edm::LogInfo("CepGenEventGenerator") << "Event modifier modules: " << modif_mods << ".";
+  for (const auto& mod : modif_mods.keys())
+    modif_modules_.emplace_back(std::make_pair(mod, modif_mods.get<cepgen::ParametersList>(mod)));
+
+  const auto output_mods =
+      cepgen::fromParameterSet(iConfig.getUntrackedParameter<edm::ParameterSet>("outputModules", edm::ParameterSet{}));
+  edm::LogInfo("CepGenEventGenerator") << "Output modules: " << output_mods << ".";
+  for (const auto& mod : output_mods.keys())
+    output_modules_.emplace_back(std::make_pair(mod, output_mods.get<cepgen::ParametersList>(mod)));
+}
+
+CepGenEventGenerator::~CepGenEventGenerator() { edm::LogInfo("CepGenEventGenerator") << "Destructor called."; }
+
+bool CepGenEventGenerator::initializeForInternalPartons() {
+  gen_ = new cepgen::Generator(true /* "safe" mode: start without plugins */);
+
+  auto pproc = proc_params_;
+  {  // little treatment to allow for standard CepGen configurations to be copy-pasted in place
+    pproc += proc_params_.get<cepgen::ParametersList>("processParameters");
+    pproc.erase("processParameters");
+    auto& pkin = pproc.operator[]<cepgen::ParametersList>("kinematics");
+    pkin += pproc.get<cepgen::ParametersList>("inKinematics");
+    pproc.erase("inKinematics");
+    pkin += pproc.get<cepgen::ParametersList>("outKinematics");
+    pproc.erase("outKinematics");
+    if (pproc.has<unsigned long long>("mode"))
+      pkin.set<int>("mode", pproc.get<unsigned long long>("mode"));
+  }
+
+  gen_->runParameters().setProcess(cepgen::ProcessFactory::get().build(pproc));
+  if (!gen_->runParameters().hasProcess())
+    throw cms::Exception("CepGenEventGenerator") << "Failed to retrieve a process from the configuration";
+  for (const auto& mod : modif_modules_) {
+    auto modifier = cepgen::EventModifierFactory::get().build(mod.first, mod.second);
+    for (const auto& cfg : mod.second.get<std::vector<std::string> >("preConfiguration"))
+      modifier->readString(cfg);
+    for (const auto& cfg : mod.second.get<std::vector<std::string> >("processConfiguration"))
+      modifier->readString(cfg);
+    gen_->runParameters().addModifier(std::move(modifier));
+  }
+  for (const auto& mod : output_modules_)
+    gen_->runParameters().addEventExporter(cepgen::EventExporterFactory::get().build(mod.first, mod.second));
+
+  edm::LogInfo("CepGenEventGenerator") << "Run parameters:\n" << gen_->runParameters();
+  const auto xsec = gen_->computeXsection();
+  xsec_.set_cross_section(xsec, xsec.uncertainty());
+  runInfo().setInternalXSec(GenRunInfoProduct::XSec(xsec, xsec.uncertainty()));
+  return true;
+}
+
+bool CepGenEventGenerator::generatePartonsAndHadronize() {
+  event().reset(new HepMC::CepGenEvent(gen_->next()));
+  event()->set_cross_section(xsec_);
+  event()->weights().push_back(1.);
+  return true;
+}
diff --git a/GeneratorInterface/CepGenInterface/src/CepGenGeneratorFilter.cc b/GeneratorInterface/CepGenInterface/src/CepGenGeneratorFilter.cc
new file mode 100644
index 0000000000000..d0ca0a8ebc136
--- /dev/null
+++ b/GeneratorInterface/CepGenInterface/src/CepGenGeneratorFilter.cc
@@ -0,0 +1,4 @@
+// CepGen-CMSSW interfacing module
+//   2022-2024, Laurent Forthomme
+
+#include "GeneratorInterface/CepGenInterface/interface/CepGenGeneratorFilter.h"