INTEGRATING NOVEL GPU-ACCELERATED MONTE CARLO EVENT GENERATORS INTO THE HIGH-ENERGY PHYSICS TOOLCHAIN

INTEGRATING NOVEL GPU-ACCELERATED MONTE CARLO EVENT GENERATORS INTO THE HIGH ENERGY PHYSICS TOOLCHAIN October 2025 AUTHOR(S): Carla Judith López Zurita ETH Zürich SUPERVISOR(S): Enrico Bothmann Zenny Wettersten Stefan Roiser CERN openlab report - Oct 2025 PROJECT SPECIFICATION Monte-Carlo event generators are used in High-Energy Physics (HEP) to simulate collider events. They are a cornerstone in the physics programmes of CERN experiments such as ATLAS and CMS. Given the very high experimental requirements for precision and the strong trend of increasingly relying on hardware accelerators such as GPUs for high-performance computing, the Physics Software Engineering section of CERN’s IT department develops a new generation of novel parallelized event generators that are accelerated by GPU devices and CPU vector instructions. This project will focus on the portable parton-level event generator framework Pepper. The student will design, implement, and test the C++ and / or Python API for the Pepper framework to provide its deeper integration into the HEP toolchain. This will facilitate the use of accelerated event generation in the production of large-scale simulated event samples for the LHC and enable the use of the framework as a building block for cutting-edge theoretical calculations at next-to-leading and next-to-next-to-leading order in perturbation theory of quantum field theory. In addition, such an interface would facilitate using Pepper for generating training data for modern machine learning applications, such as studying Monte Carlo sampling based on deep neural networks or training high-fidelity surrogate models for the evaluation of computationally expensive scattering matrix elements. INTEGRATING NOVEL GPU-ACCELERATED MC EVENT GENERATORS INTO THE HEP TOOLCHAIN 1 CERN openlab report - Oct 2025 ABSTRACT Pyper is a Python interface to the portable standalone tree-level matrix element generator Pepper, a C++ framework for computing scattering amplitudes and cross sections in highenergy physics (HEP). Pyper enables Pepper’s matrix element calculations and cross section evaluations to be performed directly from Python, relying on Kokkos for abstraction on CPU and GPU backends. It also removes the previous bottleneck of I/O operations through intermediary files when developing new Python applications with Pepper. Through batch-based calculations, Pyper integrates naturally into machine learning workflows, where it can serve as a learning function within training loops, for example, in the development of more efficient Monte Carlo phase-space sampling strategies. The implementation preserves Pepper’s functionalities such as optimization, diagnostics, and timing features, while exposing a user-friendly API that makes the system more accessible and extensible, enabling Pepper’s deeper integration into the HEP toolchain. This report documents the design of the interface, the underlying architecture, and the testing and validation strategies used to ensure correctness and reliability. INTEGRATING NOVEL GPU-ACCELERATED MC EVENT GENERATORS INTO THE HEP TOOLCHAIN 2 CERN openlab report - Oct 2025 TABLE OF CONTENTS 1 Introduction 4 1.1 Monte Carlo event generators . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.2 The Pepper framework ............................... 5 1.3 Currentlimitations.................................. 6 2 Python bindings for Pepper 6 2.1 Implementation.................................... 6 2.1.1 Installation .................................. 7 2.1.2 Initialization and error handling . . . . . . . . . . . . . . . . . . . . . . . 7 2.1.3 Dataarrays .................................. 7 2.1.4 Functiondescriptions............................. 8 2.1.5 Outputsandtimings ............................. 9 2.2 TestingandValidation ................................ 9 3 Conclusion 10 4 References 12 A Appendix 14 A.1 Pepper::Pyper..................................... 14 INTEGRATING NOVEL GPU-ACCELERATED MC EVENT GENERATORS INTO THE HEP TOOLCHAIN 3 CERN openlab report - Oct 2025 1 Introduction Monte Carlo event generators (MCEGs) play a central role in high-energy physics (HEP), where they are used to simulate collider events. They form a crucial component of the physics programmes at CERN, in particular Large Hadron Collider (LHC) experiments, such as ATLAS and CMS, in the study of hadron–hadron collisions [6,18]. For both ATLAS and CMS, a significant portion of the computing time devoted to MCEGs is due to the routine production of large event samples for hard subprocesses such as pp →V+jets and pp →t¯ t+jets [1]. In recent years, the planned move towards the High Luminosity LHC phase as well as the current excellent understanding of the ATLAS and CMS detectors calls for an increased precision of theoretical predictions, such as those provided by Monte Carlo (MC) event simulations. Precision for MCEG is comprised of both the individual event precision (i.e. the order of the calculation) and the MC error, determined by the number of events that are generated for a given observable in a given region of phase space. MC event simulation campaigns for LHC physics are now limited by the poor computational performance of high-precision event generators [4]. Given the trend of growing reliance on hardware accelerators such as GPUs for high-performance computing, an effective approach to these requirements consists of shifting the computational bottleneck onto high-throughput devices. 1.1 Monte Carlo event generators An event generator is composed of the following key ingredients: •the phase space Φn, which specifies the set of quantum numbers (momenta, flavours, . . . ) of the outgoing particles, •the differential cross section dσn/dΦn, describing the transition probability to a given scattering final state, •and the value of an observable O, such as the number of jets, rapidity, or thrust, among many others. Theoretical predictions for an observable Oare expressed as ⟨O⟩ =ZdΦn dσn dΦn O(Φn).(1) This equation can be compared to the computation of an expectation value in statistics ⟨g⟩=Zdxf(x)g(x).(2) However, analytic evaluation of these integrals is generally not feasible due to the substantial manual work required for complex processes or even lack of closed-form solutions. This difficulty goes hand in hand with the fact that each observable must be treated separately, since O appears directly within the integral. Consequently, evaluation is typically performed using statistical sampling methods such as MC integration, which estimates integrals by randomly sampling points from the integration domain and averaging the function values. MC methods are particularly well suited for event generators because the phase space can easily extend to O(20) dimensions, yet MC maintains a convergence rate of O(N−1/2)that is independent of dimensionality [8]. The equation then becomes ⟨O⟩ ≈ V N N X i=1 dσn dΦnΦ(i) nOΦ(i) n,(3) INTEGRATING NOVEL GPU-ACCELERATED MC EVENT GENERATORS INTO THE HEP TOOLCHAIN 4 CERN openlab report - Oct 2025 where Nis the number of sampling points, and Vis the phase-space volume, V=RdΦn. In practice, the generated events Φ(i) nare stored together with their corresponding weights wi=dσn/dΦn(Φ(i) n), allowing the evaluation of the observable O(Φ(i) n)for different analysis without regenerating the events. This also makes the complexity of the problem independent of the observable. Typically, MCEGs use the Vegas algorithm [13] to efficiently sample the phase space. The key idea is to generate variables preferentially in regions of the integration domain that contribute most to the result, a technique known as stratified sampling. The algorithm achieves this by iteratively evaluating the function and adapting the integration grid, thus improving the accuracy of the sampling. In a typical setup, an additional accept–reject step is performed with probability P=wi/wmax, where wmax is the maximum event weight, i.e., the total cross section (the integral in Eq. 3 for O= 1). This procedure, known as unweighting or rejection sampling, produces a reduced event sample in which all events carry unit weight [15]. This results in a sample that uses a minimal number of events while still following the same physical distribution. Thus, considerably less storage is needed, and any downstream processing of the sample that scales with the number of events is computationally cheaper. This is in particular true for detector simulations, which process each event and simulate the corresponding electronic readout of the detectors of a real experiment, thus allowing for direct comparisons between simulated and real data. Detector simulations are typically computationally even more expensive compared to event generation, such that a minimal event sample is highly desirable. Therefore, in practice, large event sample production is always done with unweighting enabled. 1.2 The Pepper framework The portable standalone tree-level matrix element generator Pepper, an acronym for Portable Engine for the Production of Parton-level Event Records, is a novel framework for event generators written in C++ [3]. Pepper can be compiled for single-threaded CPU execution, for multithreaded CPU execution, and for execution on a GPU device. This is achieved with a single codebase using the Kokkos C++ Performance Portability framework [17], and MPI for distributed memory communication. In Pepper, each computing thread generates one event, which allows parallelization over independent MC events. To maximize performance on the GPU, the same partonic process group and helicity configuration are set for groups of threads/events within a given event batch and then reshuffled to remove correlation. The computation can be broken down into a series of steps that are all performed in parallel [3]: 1. Generate random numbers. 2. Generate external momenta with an optional phase-space bias. 3. Apply phase-space cuts (i.e. set the weight of an event to zero if its external momenta do not pass the cuts). 4. Evaluate phase-space sampling weight. 5. Evaluate dynamical unphysical renormalization and factorization scales. 6. Evaluate the running coupling, the sum over initial states and the corresponding partonic fluxes. 7. Sample (or sum) helicities, evaluate helicity sampling weight and calculate external polarization vectors. INTEGRATING NOVEL GPU-ACCELERATED MC EVENT GENERATORS INTO THE HEP TOOLCHAIN 5 CERN openlab report - Oct 2025 8. Evaluate amplitudes recursively and sum the squared amplitudes over color configurations. 9. Unweight events against the weight maximum (set event weight to zero if the event is rejected). 10. Optionally, project onto a leading color configuration. 11. Copy non-zero events from the device to the host. 1.3 Current limitations To obtain statistically significant results, a large number of events must be generated. However, with the current approach, event generators often suffer from very low unweighting efficiencies, which represents a major bottleneck. Unweighting efficiency typically sits at the level of 1%, but can be much lower still at higher (but still phenomenologically relevant) jet multiplicities. Recently, more modern approaches based on machine learning (ML) have been explored, for instance, deep neural networks for MC sampling, or the use of surrogate models for the evaluation of computationally expensive scattering matrix elements [5]. Such models are typically developed and trained in Python, due to both its rapid prototyping capabilities and mature frameworks such as PyTorch [16] and TensorFlow [14] for ML. Without a public interface, a user working in Python and using Pepper would need to write intermediary files to disk, perform ML calculations, write the results back to disk, and then have them read again by Pepper. This implies two full write–read cycles, making memory I/O a substantial performance bottleneck, which complicates the workflow and introduces severe inefficiencies in the overall process. An example of this phenomenon can be found in [5]. The introduction of a dedicated Python/C++ API will overcome this limitation, enabling accelerated event generation for the production of large-scale simulated samples at the LHC. Moreover, it will allow Pepper to serve as a flexible building block for state-of-the-art theoretical predictions at next-to-leading and next-to-next-to-leading order in perturbative quantum field theory. 2 Python bindings for Pepper Pybind11 is a lightweight header-only library that enables interoperability between C++ and Python by exposing C++ types in Python and vice versa [12]. Pyper, the Python bindings for Pepper, provides a Python API to the Pepper framework using Pybind11. Its purpose is to make Pepper directly accessible from Python while at the same time allowing phase-space sampling to be optional rather than tied to the default Vegas algorithm. After initialization, Pyper offers a set of functions that replicate the behavior of Pepper, enabling the computation of physical observables. Whereas Pepper originally generated uniformly distributed random numbers internally and used Vegas to perform phase-space sampling, the framework can now accept externally provided random numbers, including those drawn from custom distributions. 2.1 Implementation The objective of this section is to present the main components of Pyper and the essential concepts users should understand to use it effectively. We also describe the key considerations that guided the implementation of the interface and outline potential issues or scenarios that may result in failure. INTEGRATING NOVEL GPU-ACCELERATED MC EVENT GENERATORS INTO THE HEP TOOLCHAIN 6 CERN openlab report - Oct 2025 2.1.1 Installation Pepper is built and compiled using CMake [9]. It was straightforward to add Pyper as a new target within the existing Pepper CMake setup. However, a few considerations had to be addressed. The code must be compiled with the position-independent code compiler flag enabled, as this is required for generating shared libraries that can be dynamically loaded by Python. It was also necessary to ensure that the relevant Python installation is correctly identified, which is particularly important on systems with multiple Python versions or virtual environments. During the CMake setup, Pyper queries the current environment via a Python system call and generates distinct shared object files based on the detected version. The Pyper shared object is built using the same installation prefix as the Pepper executable when the Pybind11 library is found. Once the targets are built, the user must update their PYTHONPATH environment variable appropriately so that the library can be located and used. This can be achieved using the code snippet in Listing 1. export PYPERPATH=$(python -c "from distutils import sysconfig as sc; \ print(sc.get_python_lib(prefix='', plat_specific=True))") export PYTHONPATH="$PYPERPATH:$PYTHONPATH" Listing 1: Setting the Pyper installation path and appending it to the Python search path. 2.1.2 Initialization and error handling Special care should be taken when initializing the Pyper object, since most errors occur at this stage. In a standalone script, if an error triggers an abort, the process simply ends and prints the error message. In the context of ML training, however, this behavior can represent a serious issue, since training is often carried out in Jupyter notebooks. If an abort is triggered here, the entire Jupyter kernel crashes without printing an error message, leaving all other cells unusable and discarding any stored variables. A prominent example is the initialization of MPI and Kokkos, which is managed by the Pyper layer. Only one Pyper object can be created in a Python process, since MPI and Kokkos are initialized by the first instance and can only be initialized once per process. To mitigate this, the Pyper layer provides additional error handling that detects whether MPI is already initialized, surfacing these issues as Python errors in a Jupyter cell rather than crashing the notebook. Additional error handling has been added where possible. However, not all exceptions were caught, especially those originating within Pepper. As mentioned before, most of these issues arise during initialization and are usually caused by invalid user input. 2.1.3 Data arrays Pepper relies on Kokkos to manage arrays through its Views class and memory space specification, which provides a unified abstraction over CPU and GPU usage. Pyper was designed to ensure that the arrays allocated in Python retain their original device placement and can be used to perform calculations in C++. In practice, input arrays (random numbers) and output arrays (squared matrix elements or cross section) are created in Python and passed by reference through the process_batch function. Some copying of the input data is unavoidable, since the random numbers corresponding to phase space points are only a subset of a larger array used INTEGRATING NOVEL GPU-ACCELERATED MC EVENT GENERATORS INTO THE HEP TOOLCHAIN 7 CERN openlab report - Oct 2025 by Pepper. However, for outputs, no copying is necessary: the pointer of the output array in Pepper is redirected to the output array pointer provided by Pyper. 2.1.4 Function descriptions The list of functions that are included in the Python bindings are explicitly listed in a file similar to Listing 2. Because the bindings are explicit, it is straightforward for a developer to identify which functions belong to the Pyper library and to further extend functionality. Adding a new function to the library requires writing only a single new line in the bindings file, provided that the function is declared and implemented within the corresponding header and cpp files. The full set of functions and explicit file locations can be found in Appendix A.1. PYBIND11_MODULE(_core, m) { m.doc() ="Python bindings for python_interface.h"; py::class_<Pyper>(m, "Pyper") .def(py::init([](std::vector<std::string>args){ // constructor return createSettings(args);})) .def("batch_size",&Pyper::batch_size) // class variables .def("n_batches",&Pyper::n_batches) ... .def("optimize",&Pyper::optimize) // functions .def("process_batch",&Pyper::process_batch) .def("finalize",&Pyper::finalize); } Listing 2: Python bindings for Pepper. During the interface design, a main concern was to reuse existing library functions as much as possible, with the goal of maintaining the code’s readability and ease of maintenance. This approach led to some refactoring within Pepper, resulting in more modular and flexible functions and classes that now support a wider range of use cases. A notable example was the transformation of the original generate_batches function into three separate functions by extracting the initialization and finalization routines. This change allows generate_batch to be called as often as needed, making it directly usable for the Pyper use case and avoiding unnecessary code duplication. We provide a more detailed description of the most important functions: •constructor: The constructor of the Python object calls the function createSettings to parse the arguments from Python to C++ and to instantiate a new object of the class. Pyper’s constructor requires a list of arguments, where the first one must be the string pepper. If a runcard is specified, it should be included as −−pepper−runcard−path= <runcard_path>. To simplify usage and prevent crashes, Pyper provides a helper function, parse_cml_args. A file under the sudirectory examples in the Pepper repository contains an example that demonstrates the use of the library, which can be used running the command python main.py <runcard_path> arguments.... •optimize: This optional step follows initialization and performs tasks related to event sampling optimization. These include adjusting the Vegas bins (if enabled), setting the helicity-configuration weights, and determining the channel-selection weights for the various subprocesses. Once the optimization is complete, the Vegas grids and selection INTEGRATING NOVEL GPU-ACCELERATED MC EVENT GENERATORS INTO THE HEP TOOLCHAIN 8