diff --git a/Cards/madgraph_process_aamumu_cfg.py b/Cards/madgraph_process_aamumu_cfg.py new file mode 100644 index 000000000..471ee3c55 --- /dev/null +++ b/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) diff --git a/Cards/madgraph_process_aamumu_coll_cfg.py b/Cards/madgraph_process_aamumu_coll_cfg.py new file mode 100644 index 000000000..1041c81b2 --- /dev/null +++ b/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) diff --git a/Cards/madgraph_process_aattbar_cfg.py b/Cards/madgraph_process_aattbar_cfg.py new file mode 100644 index 000000000..a7263affa --- /dev/null +++ b/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) diff --git a/Cards/madgraph_process_aattbar_coll_cfg.py b/Cards/madgraph_process_aattbar_coll_cfg.py new file mode 100644 index 000000000..2bd022423 --- /dev/null +++ b/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) diff --git a/Cards/madgraph_process_cfg.py b/Cards/madgraph_process_aaww_cfg.py similarity index 97% rename from Cards/madgraph_process_cfg.py rename to Cards/madgraph_process_aaww_cfg.py index 57583997f..ac087c749 100644 --- a/Cards/madgraph_process_cfg.py +++ b/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.,), ), diff --git a/Cards/madgraph_process_aaww_coll_cfg.py b/Cards/madgraph_process_aaww_coll_cfg.py new file mode 100644 index 000000000..078a830e7 --- /dev/null +++ b/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) diff --git a/CepGen/Process/CollinearPhaseSpaceGenerator.cpp b/CepGen/Process/CollinearPhaseSpaceGenerator.cpp index 8c90d3ce2..9434f832c 100644 --- a/CepGen/Process/CollinearPhaseSpaceGenerator.cpp +++ b/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& 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{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"); diff --git a/CepGen/Process/KTPhaseSpaceGenerator.cpp b/CepGen/Process/KTPhaseSpaceGenerator.cpp index 5a5b0e7d8..079e6af8d 100644 --- a/CepGen/Process/KTPhaseSpaceGenerator.cpp +++ b/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& 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().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 << "."; diff --git a/CepGenAddOns/CMakeLists.txt b/CepGenAddOns/CMakeLists.txt index 1393fbf86..5e0897548 100644 --- a/CepGenAddOns/CMakeLists.txt +++ b/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) diff --git a/CepGenAddOns/MadGraphWrapper/CMakeLists.txt b/CepGenAddOns/MadGraphWrapper/CMakeLists.txt new file mode 100644 index 000000000..157375bd3 --- /dev/null +++ b/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) diff --git a/CepGenAddOns/MadGraphWrapper/MadGraphDummyProcess.hxx b/CepGenAddOns/MadGraphWrapper/MadGraphDummyProcess.hxx new file mode 100644 index 000000000..74b7b7191 --- /dev/null +++ b/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 . + */ + +// 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 + diff --git a/CepGenAddOns/MadGraphWrapper/MadGraphInterface.cpp b/CepGenAddOns/MadGraphWrapper/MadGraphInterface.cpp new file mode 100644 index 000000000..7e77dcb89 --- /dev/null +++ b/CepGenAddOns/MadGraphWrapper/MadGraphInterface.cpp @@ -0,0 +1,303 @@ +/* + * CepGen: a central exclusive processes event generator + * Copyright (C) 2020-2023 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 . + */ + +#ifndef MADGRAPH_BIN +#error "*** MADGRAPH_BIN variable not set! ***" +#endif +#ifndef MADGRAPH_PROC_TMPL +#error "*** MADGRAPH_PROC_TMPL variable not set! ***" +#endif +#ifndef CC_CFLAGS +#error "*** CC_CFLAGS variable not set! ***" +#endif + +#include +#include + +#include "CepGen/Core/Exception.h" +#include "CepGen/Physics/PDG.h" +#include "CepGen/Utils/Caller.h" +#include "CepGen/Utils/String.h" +#include "CepGenAddOns/MadGraphWrapper/MadGraphInterface.h" +#include "CepGenAddOns/MadGraphWrapper/MadGraphProcess.h" + +namespace cepgen { + std::unordered_map MadGraphInterface::mg5_parts_ = { + {"d", (pdgid_t)1}, {"d~", (pdgid_t)1}, {"u", (pdgid_t)2}, {"u~", (pdgid_t)2}, {"s", (pdgid_t)3}, + {"s~", (pdgid_t)3}, {"c", (pdgid_t)4}, {"c~", (pdgid_t)4}, {"b", (pdgid_t)5}, {"b~", (pdgid_t)5}, + {"t", (pdgid_t)6}, {"t~", (pdgid_t)6}, {"e+", (pdgid_t)11}, {"e-", (pdgid_t)11}, {"ve", (pdgid_t)12}, + {"ve~", (pdgid_t)12}, {"mu+", (pdgid_t)13}, {"mu-", (pdgid_t)13}, {"vm", (pdgid_t)14}, {"vm~", (pdgid_t)14}, + {"tau+", (pdgid_t)15}, {"tau-", (pdgid_t)15}, {"vt", (pdgid_t)16}, {"vt~", (pdgid_t)16}, {"g", (pdgid_t)21}, + {"a", (pdgid_t)22}, {"z", (pdgid_t)23}, {"w+", (pdgid_t)24}, {"w-", (pdgid_t)24}, {"h", (pdgid_t)25}, + }; + + MadGraphInterface::MadGraphInterface(const ParametersList& params) + : SteeredObject(params), + proc_(steer("process")), + model_(steer("model")), + tmp_dir_(steerAs("tmpDir")), + card_path_(steerAs("cardPath")), + log_filename_(steer("logFile")), + standalone_cpp_path_(steerAs("standaloneCppPath")), + extra_particles_(steer("extraParticles")) { + if (proc_.empty() && standalone_cpp_path_.empty()) + throw CG_FATAL("MadGraphInterface") << "Neither a 'process' keyword nor a path to a MadGraph process interface " + "already generated ('standaloneCppPath') was set to the parameters!\n" + << params; + std::ofstream log(log_filename_, std::ios::trunc); // clearing the log + parseExtraParticles(); + } + + void MadGraphInterface::parseExtraParticles() { + for (const auto& extra_part : extra_particles_.keysOf()) { + if (mg5_parts_.count(extra_part) > 0) + throw CG_FATAL("MadGraphInterface") + << "Particle with name '" << extra_part << "' is already defined in internal LUT."; + const auto& extra_part_prop = extra_particles_.get(extra_part); + // find the equivalent MadGraph particle to alias + std::string found_mg_equiv; + for (const auto& part : mg5_parts_) + if (part.second == extra_part_prop.pdgid) + found_mg_equiv = part.first; + if (found_mg_equiv.empty()) + throw CG_FATAL("MadGraphInterface") + << "No equivalent for particle with PDG id=" << extra_part_prop.pdgid << " in MadGraph LUT."; + if (found_mg_equiv.at(found_mg_equiv.size() - 1) == '+' || found_mg_equiv.at(found_mg_equiv.size() - 1) == '-') + found_mg_equiv.pop_back(); + if (extra_part_prop.charge != 0) { + extra_part_definitions_ += "\ndefine " + extra_part + "+ = " + found_mg_equiv + "+"; + extra_part_definitions_ += "\ndefine " + extra_part + "- = " + found_mg_equiv + "-"; + mg5_parts_[extra_part + "+"] = mg5_parts_.at(found_mg_equiv + "+"); + mg5_parts_[extra_part + "-"] = mg5_parts_.at(found_mg_equiv + "-"); + } else { + extra_part_definitions_ += "\ndefine " + extra_part + " = " + found_mg_equiv; + mg5_parts_[extra_part] = mg5_parts_.at(found_mg_equiv); + } + //FIXME add extra particles properties (masses, ...) + CG_DEBUG("MadGraphInterface") << "Defined '" << extra_part + << "' as MadGraph alias for CepGen particle with properties: " << extra_part_prop + << "."; + } + CG_LOG << extra_part_definitions_; + } + + std::string MadGraphInterface::run() const { + std::ofstream log(log_filename_, std::ios::app); // appending at the end of the log + + fs::path cpp_path, cg_proc; + if (!standalone_cpp_path_.empty()) { + CG_INFO("MadGraphInterface:run") << "Running on a process already generated by mg5_aMC: " << standalone_cpp_path_; + cpp_path = standalone_cpp_path_; + cg_proc = tmp_dir_ / "cepgen_proc_interface.cpp"; + } else { + CG_INFO("MadGraphInterface:run") << "Running the mg5_aMC process generation."; + prepareCard(); + cpp_path = tmp_dir_; + const auto num_removed_files = fs::remove_all(cpp_path); + CG_DEBUG("MadGraphInterface:run") << "Removed " << utils::s("file", num_removed_files, true) + << " from process directory " << cpp_path << "."; + generateProcess(); + + CG_INFO("MadGraphInterface:run") << "Preparing the mg5_aMC process library."; + log << "\n\n*** mg5_aMC process library compilation ***\n\n"; + cg_proc = prepareMadGraphProcess(); + } + +#ifdef _WIN32 + fs::path lib_path = "CepGenMadGraphProcess.dll"; +#else + fs::path lib_path = "libCepGenMadGraphProcess.so"; +#endif + + generateLibrary(cg_proc, cpp_path, lib_path); + + CG_INFO("MadGraphInterface:run") << "Creating links for all cards in current directory."; + linkCards(); + + return lib_path; + } + + void MadGraphInterface::prepareCard() const { + std::ofstream card(card_path_); + if (!model_.empty()) + card << "set auto_convert_model T\n" + << "import model " << model_ << "\n"; + card << extra_part_definitions_ << "\n" + << "generate " << proc_ << "\n" + << "output standalone_cpp " << tmp_dir_.string() << "\n" + << "exit\n"; + card.close(); + } + + void MadGraphInterface::linkCards() const { + for (const auto& f : fs::directory_iterator(tmp_dir_ / "Cards")) + if (f.path().extension() == ".dat") { + fs::path link_path = f.path().filename(); + if (!fs::exists(link_path)) + fs::create_symlink(f, link_path); + } + } + + std::string MadGraphInterface::prepareMadGraphProcess() const { + //--- open template file + std::ifstream tmpl_file(MADGRAPH_PROC_TMPL); + std::string tmpl = std::string(std::istreambuf_iterator(tmpl_file), std::istreambuf_iterator()); + std::ofstream log(log_filename_, std::ios::app); // appending at the end of the log + log << "\n\n*** mg5_aMC process library compilation ***\n\n"; + + const auto& parts = unpackProcessParticles(proc_); + CG_INFO("MadGraphInterface.prepareMadGraphProcess") + << "Unpacked process particles: " + << "incoming=" << std::vector(parts.first.begin(), parts.first.end()) << ", " + << "outgoing=" << std::vector(parts.second.begin(), parts.second.end()) << "."; + + const auto& in_parts = parts.first; + utils::replace_all(tmpl, "XXX_PART1_XXX", std::to_string(in_parts[0])); + utils::replace_all(tmpl, "XXX_PART2_XXX", std::to_string(in_parts[1])); + + const auto& out_parts = parts.second; + std::ostringstream outparts_str; + std::string sep; + for (const auto& op : out_parts) + outparts_str << sep << std::to_string(op), sep = ", "; + utils::replace_all(tmpl, "XXX_OUT_PART_XXX", outparts_str.str()); + + utils::replace_all(tmpl, "XXX_PROC_NAME_XXX", proc_); + std::string descr = proc_; + if (!model_.empty()) + descr += " (model: " + model_ + ")"; + utils::replace_all(tmpl, "XXX_PROC_DESCRIPTION_XXX", descr); + + std::string src_filename = tmp_dir_ / "cepgen_proc_interface.cpp"; + std::ofstream src_file(src_filename); + src_file << tmpl; + src_file.close(); + return src_filename; + } + + //--------------- static utilities --------------- + + MadGraphInterface::ProcessParticles MadGraphInterface::unpackProcessParticles(const std::string& proc) { + ProcessParticles out; + auto trim_all = [](std::vector coll) -> std::vector { + std::for_each(coll.begin(), coll.end(), [](std::string& it) { it = utils::trim(it); }); + return coll; + }; + //--- dirty fix to specify incoming- and outgoing states + // as extracted from the mg5_aMC process string + const auto prim_proc = utils::split(utils::trim(proc), ',')[0]; + auto parts = trim_all(utils::split(prim_proc, '>')); + if (parts.size() != 2) + throw CG_FATAL("MadGraphInterface:unpackProcessParticles") + << "Unable to unpack particles from process name: \"" << proc << "\" -> " << parts << "!"; + //--- incoming parton-like particles + auto prim_parts = trim_all(utils::split(parts[0], ' ')); + CG_DEBUG("MadGraphInterface:unpackProcessParticles") << "Primary particles: " << prim_parts; + if (prim_parts.size() != 2) + throw CG_FATAL("MadGraphInterface:unpackProcessParticles") + << "Unable to unpack particles from primary particles list: \"" << parts[0] << "\" -> " << prim_parts << "!"; + for (const auto& p : prim_parts) { + if (mg5_parts_.count(p) == 0) + throw CG_FATAL("MadGraphInterface:unpackProcessParticles") + << "Particle with mg5_aMC name '" << p << "' was not recognised!"; + out.first.emplace_back(mg5_parts_.at(p)); + } + //---- outgoing system + auto dec_parts = trim_all(utils::split(parts[1], ' ')); + CG_DEBUG("MadGraphInterface:unpackProcessParticles") << "Outgoing system: " << dec_parts; + for (auto& p : dec_parts) { + if (mg5_parts_.count(p) == 0) + throw CG_FATAL("MadGraphInterface:unpackProcessParticles") + << "Particle with mg5_aMC name '" << p << "' was not recognised!"; + out.second.emplace_back(mg5_parts_.at(p)); + } + return out; + } + + void MadGraphInterface::generateProcess() const { + std::ofstream log(log_filename_, std::ios::app); // appending at the end of the log + log << "\n\n*** mg5_aMC process generation ***\n\n"; + log << utils::Caller::call({MADGRAPH_BIN, "-f", card_path_.string()}); + fs::remove(card_path_); + } + + void MadGraphInterface::generateLibrary(const fs::path& proc_path, + const fs::path& in_path, + const fs::path& out_lib) const { + std::vector src_files; + src_files.emplace_back(proc_path.string()); + + //--- find all processes registered + std::vector processes; + try { + for (const auto& p : fs::directory_iterator(in_path / "SubProcesses")) + if (p.path().filename().string()[0] == 'P') { + processes.emplace_back(p.path()); + for (const auto& f : fs::directory_iterator(p)) + if (f.path().extension() == ".cc") + src_files.emplace_back(f.path()); + } + } catch (const fs::filesystem_error& err) { + throw CG_FATAL("MadGraphInterface:generateLibrary") + << "Failed to retrieve all subprocesses in path " << in_path << "!\n" + << err.what(); + } + + CG_DEBUG("MadGraphInterface:generateLibrary") << "Subprocess list: " << processes << "."; + + if (processes.size() != 1) + throw CG_FATAL("MadGraphInterface:generateLibrary") << "Currently only single-process cases are supported!"; + + //--- find all model source files + for (const auto& f : fs::directory_iterator(in_path / "src")) + if (f.path().extension() == ".cc") + src_files.emplace_back(f.path()); + +#ifdef _WIN32 + throw CG_FATAL("MadGraphInterface:generateLibrary") + << "Library generation not yet implemented for Window$ systems!"; +#else + utils::Caller::call({CC_CFLAGS, + "-fPIC", + "-shared", + "-Wno-unused-variable", + "-Wno-int-in-bool-context", + "-I" + (in_path / "src").string(), + "-I" + processes.at(0), + utils::merge(src_files, " "), + "-o " + out_lib.string()}); +#endif + } + + ParametersDescription MadGraphInterface::description() { + auto desc = ParametersDescription(); + desc.add("process", "").setDescription("MadGraph_aMC process definition"); + desc.add("model", "sm-full").setDescription("MadGraph_aMC model name"); + desc.add("cardPath", "/tmp/cepgen_mg5_input.dat") + .setDescription("Temporary file where to store the input card for MadGraph_aMC"); + desc.add("standaloneCppPath", ""); + desc.add("tmpDir", "/tmp/cepgen_mg5_aMC") + .setDescription("Temporary path where to store the MadGraph_aMC process definition files"); + desc.add("logFile", "/tmp/cepgen_mg5_aMC.log") + .setDescription("Temporary path where to store the log for this run"); + desc.add("extraParticles", ParametersDescription()) + .setDescription("define internal MadGraph alias for a particle name"); + + return desc; + } +} // namespace cepgen diff --git a/CepGenAddOns/MadGraphWrapper/MadGraphInterface.h b/CepGenAddOns/MadGraphWrapper/MadGraphInterface.h new file mode 100644 index 000000000..26fddaa51 --- /dev/null +++ b/CepGenAddOns/MadGraphWrapper/MadGraphInterface.h @@ -0,0 +1,67 @@ +/* + * 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 . + */ + +#ifndef CepGenAddOns_MadGraphWrapper_MadGraphInterface_h +#define CepGenAddOns_MadGraphWrapper_MadGraphInterface_h + +#include +#include + +#include "CepGen/Core/SteeredObject.h" +#include "CepGen/Physics/ParticleProperties.h" +#include "CepGen/Utils/Filesystem.h" + +// forward-declaration of base MadGraph standalone_cpp process +class CPPProcess; + +namespace cepgen { + class MadGraphInterface : public SteeredObject { + public: + explicit MadGraphInterface(const ParametersList&); + + static ParametersDescription description(); + + std::string run() const; + + private: + static constexpr size_t cmd_buffer_size_ = 256; + static std::unordered_map mg5_parts_; + + using ProcessParticles = std::pair, std::vector >; + static ProcessParticles unpackProcessParticles(const std::string&); + + void generateProcess() const; + void generateLibrary(const fs::path&, const fs::path&, const fs::path&) const; + void parseExtraParticles(); + void prepareCard() const; + void linkCards() const; + std::string prepareMadGraphProcess() const; + + const std::string proc_; + const std::string model_; + const fs::path tmp_dir_; + const fs::path card_path_; + const fs::path log_filename_; + const fs::path standalone_cpp_path_; + const ParametersList extra_particles_; + + std::string extra_part_definitions_; + }; +} // namespace cepgen + +#endif diff --git a/CepGenAddOns/MadGraphWrapper/MadGraphProcess.h b/CepGenAddOns/MadGraphWrapper/MadGraphProcess.h new file mode 100644 index 000000000..251bdcaf9 --- /dev/null +++ b/CepGenAddOns/MadGraphWrapper/MadGraphProcess.h @@ -0,0 +1,68 @@ +/* + * 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 . + */ + +#ifndef CepGenAddOns_MadGraphWrapper_MadGraphProcess_h +#define CepGenAddOns_MadGraphWrapper_MadGraphProcess_h + +#include + +#include "CepGen/Physics/Momentum.h" + +// forward-declaration of base MadGraph standalone_cpp process +class CPPProcess; + +namespace cepgen { + /// Wrapper around a generic MadGraph CPPProcess definition + class MadGraphProcess { + public: + MadGraphProcess(); + ~MadGraphProcess(); + + const std::string& name() const { return name_; } + const std::string& description() const { return descr_; } + + void initialise(const std::string&); + inline const std::array& intermediatePartons() const { return incoming_pdgids_; } + inline const std::vector& centralSystem() const { return central_pdgids_; } + double eval(); + + inline MadGraphProcess& setMomentum(size_t i, const Momentum& mom) { + if (i > mom_.size()) + throw CG_FATAL("MadGraphProcess") << "Invalid index for momentum: " << i << "!"; + mom_[i][0] = mom.energy(); + mom_[i][1] = mom.px(); + mom_[i][2] = mom.py(); + mom_[i][3] = mom.pz(); + return *this; + } + const std::vector& momenta(); + + private: + std::unique_ptr proc_; + std::vector mom_; + + const std::string name_; + const std::string descr_; + const std::array incoming_pdgids_; + const std::vector central_pdgids_; + + std::vector momenta_; + }; +} // namespace cepgen + +#endif diff --git a/CepGenAddOns/MadGraphWrapper/MadGraphProcess.tpp b/CepGenAddOns/MadGraphWrapper/MadGraphProcess.tpp new file mode 100644 index 000000000..8ad16c2be --- /dev/null +++ b/CepGenAddOns/MadGraphWrapper/MadGraphProcess.tpp @@ -0,0 +1,88 @@ +/* + * CepGen: a central exclusive processes event generator + * Copyright (C) 2020-2023 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 . + */ + +//============================================================================= +// NOLI SE TANGERE +#include "CPPProcess.h" +#include "CepGen/Core/Exception.h" +#include "CepGenAddOns/MadGraphWrapper/MadGraphProcess.h" + +using namespace cepgen; + +MadGraphProcess::MadGraphProcess() + : proc_(new CPPProcess), + name_("XXX_PROC_NAME_XXX"), + descr_("XXX_PROC_DESCRIPTION_XXX"), + incoming_pdgids_({XXX_PART1_XXX, XXX_PART2_XXX}), + central_pdgids_({XXX_OUT_PART_XXX}) { + CG_INFO("MadGraphProcess") << "Process considered: " << proc_->name() << ". " + << "Incoming particles: " << incoming_pdgids_ << ", outgoing system: " << central_pdgids_ + << "."; +} + +MadGraphProcess::~MadGraphProcess() = default; + +void MadGraphProcess::initialise(const std::string& param_card) { + try { + proc_->initProc(param_card); + } catch (const char* chr) { + throw CG_FATAL("MadGraphProcess:init") + << "Failed to initialise parameters card at \"" << param_card << "\":\n\t" << chr; + } + if (proc_->nprocesses > 1) + throw CG_FATAL("MadGraphProcess:init") << "Multi-processes matrix elements are not (yet) supported!"; + if (proc_->ninitial != 2) + throw CG_FATAL("MadGraphProcess:init") << "Currently only 2->N processes are supported!"; + + CG_DEBUG("MadGraphProcess:init") << "External particles masses (partons + central system): " << proc_->getMasses() + << "."; + + mom_.clear(); + for (size_t i = 0; i < proc_->nexternal; ++i) + mom_.emplace_back(new double[4]{proc_->getMasses().at(i), 0., 0., 0.}); + momenta_.reserve(proc_->nexternal); +} + +double MadGraphProcess::eval() { + proc_->setMomenta(mom_); + proc_->sigmaKin(); + const double* me = proc_->getMatrixElements(); + if (me[0] <= 0) + return 0.; + + CG_DEBUG_LOOP("MadGraphProcess:eval").log([&](auto& log) { + log << "Dump of event kinematics\n\t" + << "Incoming partons 4-momenta: " << std::vector(mom_[0], mom_[0] + 4) << ", " + << std::vector(mom_[1], mom_[1] + 4) << "\n\t" + << "Outgoing particles 4-momenta: "; + std::string sep; + for (size_t i = 0; i < proc_->nexternal - 2; ++i) + log << sep << std::vector(mom_[i + 2], mom_[i + 2] + 4), sep = ", "; + log << "\n\tResulting matrix element: " << me[0] << "."; + }); + return me[0]; +} + +const std::vector& MadGraphProcess::momenta() { + const auto& p4 = proc_->getMomenta(); + // cast it to the member attribute and return it + for (size_t i = 0; i < p4.size(); ++i) + momenta_[i] = Momentum::fromPxPyPzE(p4[i][1], p4[i][2], p4[i][3], p4[i][0]); + return momenta_; +} +//============================================================================= diff --git a/CepGenAddOns/MadGraphWrapper/MadGraphProcessBuilder.cpp b/CepGenAddOns/MadGraphWrapper/MadGraphProcessBuilder.cpp new file mode 100644 index 000000000..9e42ba83d --- /dev/null +++ b/CepGenAddOns/MadGraphWrapper/MadGraphProcessBuilder.cpp @@ -0,0 +1,108 @@ +/* + * CepGen: a central exclusive processes event generator + * Copyright (C) 2020-2023 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 . + */ + +#include "CepGen/Core/Exception.h" +#include "CepGen/Event/Event.h" +#include "CepGen/Generator.h" +#include "CepGen/Modules/ProcessFactory.h" +#include "CepGen/Physics/PDG.h" +#include "CepGen/Process/Process2to4.h" +#include "CepGen/Utils/AbortHandler.h" +#include "CepGenAddOns/MadGraphWrapper/MadGraphInterface.h" +#include "CepGenAddOns/MadGraphWrapper/MadGraphProcess.h" + +using namespace cepgen; + +class MadGraphProcessBuilder : public proc::Process2to4 { +public: + explicit MadGraphProcessBuilder(const ParametersList& params, bool load_library = true) : Process2to4(params, {}) { + if (load_library) + loadMG5Library(); + // once MadGraph process library is loaded into runtime environment, can define its wrapper object + mg5_proc_.reset(new MadGraphProcess); + produced_parts_ = mg5_proc_->centralSystem(); + cs_prop_ = PDG::get()(produced_parts_.at(0)); + CG_DEBUG("MadGraphProcessBuilder") << "MadGraph_aMC process created for:\n\t" + << "* interm. parts.: " << mg5_proc_->intermediatePartons() << "\n\t" + << "* central system: " << produced_parts_ << "."; + } + + proc::ProcessPtr clone() const override { return proc::ProcessPtr(new MadGraphProcessBuilder(parameters(), false)); } + + static ParametersDescription description() { + auto desc = Process2to4::description(); + desc.setDescription("MadGraph_aMC process builder"); + desc.add("lib", "").setDescription("Precompiled library for this process definition"); + desc.add("parametersCard", "param_card.dat").setDescription("Runtime MadGraph parameters card"); + desc += MadGraphInterface::description(); + return desc; + } + + void prepareProcessKinematics() override { + const auto& interm_part = mg5_proc_->intermediatePartons(); + const auto flux_interm_part = std::array{psgen_->positiveFlux().partonPdgId(), + psgen_->negativeFlux().partonPdgId()}; + if (interm_part != flux_interm_part) + throw CG_FATAL("MadGraphProcessBuilder") + << "MadGraph unpacked process incoming state (" << interm_part << ") " + << "is incompatible with user-steered incoming fluxes particles (" << flux_interm_part << ")."; + if (const auto params_card = steer("parametersCard"); !params_card.empty()) { + CG_INFO("MadGraphProcessBuilder") << "Preparing process kinematics for card at \"" << params_card << "\"."; + mg5_proc_->initialise(params_card); + } + } + double computeCentralMatrixElement() const override { + if (!mg5_proc_) + CG_FATAL("MadGraphProcessBuilder:eval") << "Process not properly linked!"; + if (!kinematics().cuts().initial.contain(event()(Particle::Role::Parton1)) || + !kinematics().cuts().initial.contain(event()(Particle::Role::Parton2))) + return 0.; + if (!kinematics().cuts().central.contain(event()(Particle::Role::CentralSystem))) + return 0.; + + CG_DEBUG_LOOP("MadGraphProcessBuilder:eval") + << "Particles content:\n" + << "incoming: " << q1() << " (m=" << q1().mass() << "), " << q2() << " (m=" << q2().mass() << ")\n" + << "outgoing: " << pc(0) << " (m=" << pc(0).mass() << "), " << pc(1) << " (m=" << pc(1).mass() << ")."; + mg5_proc_->setMomentum(0, q1()); // first incoming parton + mg5_proc_->setMomentum(1, q2()); // second incoming parton + mg5_proc_->setMomentum(2, pc(0)); // first outgoing central particle + mg5_proc_->setMomentum(3, pc(1)); // second outgoing central particle + + return mg5_proc_->eval(); + } + +private: + void loadMG5Library() { + utils::AbortHandler(); + try { + const auto& lib_file = steer("lib"); + if (!lib_file.empty()) + loadLibrary(lib_file); + else { + const MadGraphInterface interf(params_); + loadLibrary(interf.run()); + } + } catch (const utils::RunAbortedException&) { + CG_FATAL("MadGraphProcessBuilder") << "MadGraph_aMC process generation aborted."; + } + } + + std::unique_ptr mg5_proc_; +}; +REGISTER_PROCESS("mg5_aMC", MadGraphProcessBuilder); diff --git a/CepGenAddOns/MadGraphWrapper/test/add_particles.cc b/CepGenAddOns/MadGraphWrapper/test/add_particles.cc new file mode 100644 index 000000000..fddf44a93 --- /dev/null +++ b/CepGenAddOns/MadGraphWrapper/test/add_particles.cc @@ -0,0 +1,26 @@ +#include "CepGen/Core/Exception.h" +#include "CepGen/Generator.h" +#include "CepGen/Modules/ProcessFactory.h" +#include "CepGen/Physics/PDG.h" +#include "CepGen/Process/Process.h" +#include "CepGen/Utils/ArgumentsParser.h" +#include "CepGenAddOns/MadGraphWrapper/MadGraphDummyProcess.hxx" +#include "CepGenAddOns/MadGraphWrapper/MadGraphProcess.h" + +int main(int argc, char* argv[]) { + cepgen::ArgumentsParser(argc, argv).parse(); + cepgen::initialise(); + + const cepgen::pdgid_t my_part = 13; + + cepgen::PDG::get().define(cepgen::ParticleProperties(my_part, "la", "laurentino", 0., 42., 0., 3, true)); + CG_LOG << cepgen::PDG::get()(my_part); + + auto mg5 = cepgen::ProcessFactory::get().build( + "mg5_aMC", + cepgen::ParametersList() + .set("extraParticles", cepgen::ParametersList().set("la", cepgen::PDG::get()(my_part))) + .set("process", "a a > la+ la-")); + + return 0; +}