MadGraph5_aMC@NLO interfacing (#34)

* Core: uniformised collinear/kt PS generator implementations

* MG5: added an interfacing plugin for external matrix elements computation

* MG5: add capability to define particles aliases

* MG5: making process a collinear one

* MG5: cleanup of the initialisation procedure

* MG5: fixed the particles ID retrieval from process, added new two-photon MG5 processes

* MG5: both kt and collinear example cards for the 3 handled processes

Cards/madgraph_process_aamumu_cfg.py

@@ -0,0 +1,41 @@
+import Config.Core as cepgen
+import Config.ktProcess_cfi as kt
+from Config.generator_cff import generator
+#from Config.Logger_cfi import logger
+#logger.enabledModules += ('MadGraphProcess.eval',)
+
+#--- process definition
+process = kt.process.clone('mg5_aMC',
+    processParameters = cepgen.Parameters(
+        process = 'a a > mu+ mu-',
+        # alternatively, if shared object is already generated
+        #lib = 'libCepGenMadGraphProcess.so',
+        # alternatively, if standalone_cpp directory is already generated
+        #standaloneCppPath = '/tmp/cepgen_mg5_aMC',
+        mode = cepgen.ProcessMode.ElasticElastic,
+    ),
+    inKinematics = cepgen.Parameters(
+        pz = (6500., 6500.),
+        structureFunctions = cepgen.StructureFunctions.LUXlike,
+    ),
+    outKinematics = kt.process.outKinematics.clone(
+        qt = (0., 10.),
+        #eta = (-2.5, 2.5),
+        mx = (1.07, 1000.),
+        pt = (25.,),
+    ),
+)
+
+#--- events generation
+generator.numEvents = 10000
+
+text = cepgen.Module('text',  # histogramming/ASCII output capability
+    #variables = ['nev', 'm(4)', 'tgen'],
+    histVariables={
+        'm(4)': cepgen.Parameters(xrange=(0., 250.), nbins=10, log=True),
+        'pt(4)': cepgen.Parameters(xrange=(0., 25.), nbins=10, log=True),
+        'pt(7):pt(8)': cepgen.Parameters(xrange=(0., 250.), yrange=(0., 250.), log=True)
+    }
+)
+dump = cepgen.Module('dump', printEvery = generator.printEvery)
+output = cepgen.Sequence(text, dump)

Cards/madgraph_process_aamumu_coll_cfg.py

@@ -0,0 +1,34 @@
+import Config.Core as cepgen
+import Config.collinearProcess_cfi as coll
+from Config.generator_cff import generator
+
+#--- process definition
+process = coll.process.clone('mg5_aMC',
+    processParameters = cepgen.Parameters(
+        process = 'a a > mu+ mu-',
+        mode = cepgen.ProcessMode.ElasticElastic,
+    ),
+    inKinematics = cepgen.Parameters(
+        pz = (6500., 6500.),
+        structureFunctions = cepgen.StructureFunctions.LUXlike,
+    ),
+    outKinematics = coll.process.outKinematics.clone(
+        q2 = (0., 10.),
+        #eta = (-2.5, 2.5),
+        mx = (1.07, 1000.),
+        pt = (25.,),
+    ),
+)
+
+#--- events generation
+generator.numEvents = 10000
+
+text = cepgen.Module('text',  # histogramming/ASCII output capability
+    histVariables={
+        'm(4)': cepgen.Parameters(xrange=(0., 250.), nbins=10, log=True),
+        'pt(4)': cepgen.Parameters(xrange=(0., 25.), nbins=10, log=True),
+        'pt(7):pt(8)': cepgen.Parameters(xrange=(0., 250.), yrange=(0., 250.), log=True)
+    }
+)
+dump = cepgen.Module('dump', printEvery = generator.printEvery)
+output = cepgen.Sequence(text, dump)

Cards/madgraph_process_aattbar_cfg.py

@@ -0,0 +1,38 @@
+import Config.Core as cepgen
+import Config.ktProcess_cfi as kt
+from Config.generator_cff import generator
+
+#--- process definition
+process = kt.process.clone('mg5_aMC',
+    processParameters = cepgen.Parameters(
+        process = 'a a > t t~',
+        # alternatively, if shared object is already generated
+        #lib = 'libCepGenMadGraphProcess.so',
+        # alternatively, if standalone_cpp directory is already generated
+        #standaloneCppPath = '/tmp/cepgen_mg5_aMC',
+        mode = cepgen.ProcessMode.ElasticElastic,
+    ),
+    inKinematics = cepgen.Parameters(
+        pz = (6500., 6500.),
+        structureFunctions = cepgen.StructureFunctions.LUXlike,
+    ),
+    outKinematics = kt.process.outKinematics.clone(
+        qt = (0., 10.),
+        #eta = (-2.5, 2.5),
+        mx = (1.07, 1000.),
+        pt = (0.,),
+    ),
+)
+
+#--- events generation
+generator.numEvents = 10000
+
+text = cepgen.Module('text',  # histogramming/ASCII output capability
+    histVariables={
+        'm(4)': cepgen.Parameters(xrange=(0., 250.), nbins=10, log=True),
+        'pt(4)': cepgen.Parameters(xrange=(0., 25.), nbins=10, log=True),
+        'pt(7):pt(8)': cepgen.Parameters(xrange=(0., 250.), yrange=(0., 250.), log=True)
+    }
+)
+dump = cepgen.Module('dump', printEvery = generator.printEvery)
+output = cepgen.Sequence(text, dump)

Cards/madgraph_process_aattbar_coll_cfg.py

@@ -0,0 +1,34 @@
+import Config.Core as cepgen
+import Config.collinearProcess_cfi as coll
+from Config.generator_cff import generator
+
+#--- process definition
+process = coll.process.clone('mg5_aMC',
+    processParameters = cepgen.Parameters(
+        process = 'a a > t t~',
+        mode = cepgen.ProcessMode.ElasticElastic,
+    ),
+    inKinematics = cepgen.Parameters(
+        pz = (6500., 6500.),
+        structureFunctions = cepgen.StructureFunctions.LUXlike,
+    ),
+    outKinematics = coll.process.outKinematics.clone(
+        q2 = (0., 10.),
+        #eta = (-2.5, 2.5),
+        mx = (1.07, 1000.),
+        pt = (0.,),
+    ),
+)
+
+#--- events generation
+generator.numEvents = 10000
+
+text = cepgen.Module('text',  # histogramming/ASCII output capability
+    histVariables={
+        'm(4)': cepgen.Parameters(xrange=(0., 250.), nbins=10, log=True),
+        'pt(4)': cepgen.Parameters(xrange=(0., 25.), nbins=10, log=True),
+        'pt(7):pt(8)': cepgen.Parameters(xrange=(0., 250.), yrange=(0., 250.), log=True)
+    }
+)
+dump = cepgen.Module('dump', printEvery = generator.printEvery)
+output = cepgen.Sequence(text, dump)

Cards/madgraph_process_cfg.py → Cards/madgraph_process_aaww_cfg.py

@@ -5,7 +5,7 @@
 #--- process definition
 process = kt.process.clone('mg5_aMC',
     processParameters = cepgen.Parameters(
-        process = 'a a > mu+ mu-',
+        process = 'a a > w+ w-',
         # alternatively, if shared object is already generated
         #lib = 'libCepGenMadGraphProcess.so',
         # alternatively, if standalone_cpp directory is already generated
@@ -17,8 +17,8 @@
         structureFunctions = cepgen.StructureFunctions.LUXlike,
     ),
     outKinematics = kt.process.outKinematics.clone(
-        #eta = (-2.5, 2.5),
         qt = (0., 10.),
+        #eta = (-2.5, 2.5),
         mx = (1.07, 1000.),
         pt = (0.,),
     ),

Cards/madgraph_process_aaww_coll_cfg.py

@@ -0,0 +1,34 @@
+import Config.Core as cepgen
+import Config.collinearProcess_cfi as coll
+from Config.generator_cff import generator
+
+#--- process definition
+process = coll.process.clone('mg5_aMC',
+    processParameters = cepgen.Parameters(
+        process = 'a a > w+ w-',
+        mode = cepgen.ProcessMode.ElasticElastic,
+    ),
+    inKinematics = cepgen.Parameters(
+        pz = (6500., 6500.),
+        structureFunctions = cepgen.StructureFunctions.LUXlike,
+    ),
+    outKinematics = coll.process.outKinematics.clone(
+        q2 = (0., 10.),
+        #eta = (-2.5, 2.5),
+        mx = (1.07, 1000.),
+        pt = (0.,),
+    ),
+)
+
+#--- events generation
+generator.numEvents = 10000
+
+text = cepgen.Module('text',  # histogramming/ASCII output capability
+    histVariables={
+        'm(4)': cepgen.Parameters(xrange=(0., 250.), nbins=10, log=True),
+        'pt(4)': cepgen.Parameters(xrange=(0., 25.), nbins=10, log=True),
+        'pt(7):pt(8)': cepgen.Parameters(xrange=(0., 250.), yrange=(0., 250.), log=True)
+    }
+)
+dump = cepgen.Module('dump', printEvery = generator.printEvery)
+output = cepgen.Sequence(text, dump)

CepGen/Process/CollinearPhaseSpaceGenerator.cpp

@@ -30,7 +30,9 @@ namespace cepgen {
     CollinearPhaseSpaceGenerator::CollinearPhaseSpaceGenerator(Process* proc) : PhaseSpaceGenerator(proc) {}
 
     void CollinearPhaseSpaceGenerator::initialise() {
-      auto& kin = process().kinematics();
+      const auto& kin = process().kinematics();
+
+      // pick a parton flux parameterisation for each beam
       auto set_flux_properties = [&kin](const Beam& beam, std::unique_ptr<PartonFlux>& flux) {
         auto params = beam.partonFluxParameters();
         const auto params_p_el = CollinearFluxFactory::get().describeParameters(
@@ -57,15 +59,13 @@ namespace cepgen {
         if (!flux)
           throw CG_FATAL("CollinearPhaseSpaceGenerator:init")
               << "Failed to initiate a parton flux object with properties: " << params << ".";
+        if (flux->ktFactorised())
+          throw CG_FATAL("CollinearPhaseSpaceGenerator:init")
+              << "Invalid incoming parton flux: " << flux->name() << ".";
       };
       set_flux_properties(kin.incomingBeams().positive(), pos_flux_);
       set_flux_properties(kin.incomingBeams().negative(), neg_flux_);
 
-      if (pos_flux_->ktFactorised() || neg_flux_->ktFactorised())
-        throw CG_FATAL("CollinearPhaseSpaceGenerator:init")
-            << "Invalid incoming parton fluxes: " << std::vector<std::string>{pos_flux_->name(), neg_flux_->name()}
-            << ".";
-
       // register the incoming partons' virtuality
       const auto log_lim_q2 = kin.cuts().initial.q2.truncate(Limits{1.e-10, 5.}).compute(std::log);
       process().defineVariable(m_t1_, Process::Mapping::exponential, log_lim_q2, "Positive-z parton virtuality");

CepGen/Process/KTPhaseSpaceGenerator.cpp

@@ -34,17 +34,20 @@ namespace cepgen {
       // pick a parton flux parameterisation for each beam
       auto set_flux_properties = [](const Beam& beam, std::unique_ptr<PartonFlux>& flux) {
         auto params = beam.partonFluxParameters();
+        const auto params_p_el = KTFluxFactory::get().describeParameters("BudnevElastic");
+        const auto params_p_inel = KTFluxFactory::get().describeParameters("BudnevInelastic");
+        const auto params_hi_el = KTFluxFactory::get().describeParameters("ElasticHeavyIon");
         if (params.name<std::string>().empty()) {
           if (beam.elastic()) {
             if (HeavyIon::isHI(beam.pdgId()))
-              params = KTFluxFactory::get().describeParameters("ElasticHeavyIon").validate(params);
+              params = params_hi_el.validate(params);
             else
-              params = KTFluxFactory::get().describeParameters("BudnevElastic").validate(params);
+              params = params_p_el.validate(params);
           } else
-            params = KTFluxFactory::get().describeParameters("BudnevInelastic").validate(params);
+            params = params_p_inel.validate(params);
           //TODO: fermions/pions
         }
-        flux.reset(KTFluxFactory::get().build(params).release());
+        flux = std::move(KTFluxFactory::get().build(params));
         if (!flux)
           throw CG_FATAL("KTPhaseSpaceGenerator:init")
               << "Failed to initiate a parton flux object with properties: " << params << ".";

CepGenAddOns/CMakeLists.txt

@@ -13,6 +13,7 @@ add_subdirectory(GnuplotWrapper)
 add_subdirectory(HepMC2Wrapper)
 add_subdirectory(HepMC3Wrapper)
 add_subdirectory(LHAPDFWrapper)
+add_subdirectory(MadGraphWrapper)
 add_subdirectory(MatplotlibWrapper)
 add_subdirectory(PhotosTauolaWrapper)
 add_subdirectory(ProMCWrapper)

CepGenAddOns/MadGraphWrapper/CMakeLists.txt

@@ -0,0 +1,26 @@
+#--- linking with ProMC
+
+set(MADGRAPH_DIRS $ENV{MADGRAPH_DIR} /usr /usr/local)
+find_program(MADGRAPH_BIN NAMES mg5_aMC HINTS ${MADGRAPH_DIRS} PATH_SUFFIXES bin)
+if(MADGRAPH_BIN)
+else()
+  return()
+endif()
+
+file(GLOB tmpl *.tpp)
+
+#----- build the object
+
+set(mg_defs)
+list(APPEND mg_defs "-DMADGRAPH_BIN=\"${MADGRAPH_BIN}\"")
+list(APPEND mg_defs "-DMADGRAPH_PROC_TMPL=\"${tmpl}\"")
+list(APPEND mg_defs "-DCC_CFLAGS=\"${CMAKE_CXX_COMPILER} ${CMAKE_CXX_FLAGS} -I${PROJECT_SOURCE_DIR}\"")
+
+cepgen_build(CepGenMadGraph
+    SOURCES *.cpp
+    EXT_LIBS stdc++fs
+    DEPENDS ${MADGRAPH_BIN}
+    DEFINITIONS ${mg_defs}
+    TESTS test/*.cc
+    INSTALL_COMPONENT madgraph)
+#target_compile_features(CepGenMadGraph PRIVATE cxx_range_for)

CepGenAddOns/MadGraphWrapper/MadGraphDummyProcess.hxx

@@ -0,0 +1,32 @@
+/*
+ *  CepGen: a central exclusive processes event generator
+ *  Copyright (C) 2020-2022  Laurent Forthomme
+ *
+ *  This program is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+// This file allows to link the MadGraph interfacing module without any process
+// generation performed by MG5_aMC.
+// Include it in your source file prior to any linking with libCepGenMadGraph.
+
+#include "CepGenAddOns/MadGraphWrapper/MadGraphProcess.h"
+
+class CPPProcess {};
+namespace cepgen {
+  MadGraphProcess::MadGraphProcess() : incoming_pdgids_{0, 0} {}
+  MadGraphProcess::~MadGraphProcess() {}
+  double MadGraphProcess::eval() { return 0.; }
+  void MadGraphProcess::initialise(const std::string&) {}
+}  // namespace cepgen
+