EpIC: novel Monte Carlo generator for exclusive processes

We present the EpIC Monte Carlo event generator for exclusive processes sensitive to generalised parton distributions. EpIC utilises the PARTONS framework, which provides a flexible software architecture and a variety of modelling options for the partonic description of the nucleon. The generator offers a comprehensive set of features, including multi-channel capabilities and radiative corrections. It may be used both in analyses of experimental data, as well as in impact studies, especially for future electron-ion colliders.


Introduction
The objective of this study is to develop a Monte Carlo event generator for exclusive processes involving hadrons remaining coherent during interactions with high-energy leptons. The factorisation theorems developed in the framework of quantum chromodynamics (QCD) allow one to describe such processes in terms of non-perturbative generalised parton distribution functions (GPDs) [1][2][3][4][5] convoluted with perturbatively calculable hard scattering subprocesses.
GPDs are universal, process-independent functions that parametrise the off-forward nucleon matrix elements of quark and gluon bilinear operators with lighta e-mail: elke@bnl.gov b e-mail: varvara.batozskaya@ncbj.gov.pl c e-mail: salvatore.fazio@unical.it d e-mail: k.gates.1@research.gla.ac.uk e e-mail: herve.moutarde@cea.fr f e-mail: daria.sokhan@cea.fr g e-mail: spiesber@uni-mainz.de h e-mail: pawel.sznajder@ncbj.gov.pl i e-mail: ktezgin@bnl.gov like separations. In case there is no momentum transfer to the nucleon, i.e. in the forward limit, certain GPDs become equivalent to PDFs. Additionally, the first Mellin moments of GPDs are related to elastic form factors. In this regard, GPDs may be viewed as a unified concept of elastic form factors studied via elastic scattering processes and one-dimensional parton distribution functions studied via (semi-) inclusive scattering processes. Another key aspect of GPDs is their relation to nucleon tomography. The Fourier transform of GPDs are related to the impact parameter space distributions when there is no collinear, but finite transverse momentum transfer to the nucleon [6][7][8]. This relation enables nucleon tomography, which involves the correlation between impact space distribution functions, parton polarization, and the longitudinal momentum fractions carried by partons. GPDs also exhibit a unique relationship with energy-momentum tensor (EMT) form factors, which encode fundamental properties of the nucleon, such as mass and spin decomposition of the nucleon into its constituent parts (including orbital angular momentum) [2,3,9] as well as its internal structure based on the so-called "mechanical" forces [10][11][12]. For more information on GPDs, we refer to the available reviews on the subject, like Ref. [13].
The extraction of GPDs from data on exclusive processes is not an easy task, mainly due to the difficulty of the inverse problem one must solve. Namely, typically many types of GPDs contribute to a given process, and it is necessary to deconvolute them all from the process amplitude. In order to accomplish this, measurements of many different types of processes are needed in a wide range of kinematic domains. Therefore, measurements of exclusive processes have been performed in several facilities, such as DESY, JLAB, and CERN. In addi-tion to that, GPDs are pillars of scientific programmes for a new generation of machines. This includes both electron-ion colliders, such as EIC [14,15], EIcC [16] and LHeC [17], and fixed target experiments, such as AMBER at CERN [18] and JLAB12 [19].
As part of this global effort, we have developed a Monte Carlo (MC) generator called EpIC, whose logo appears in Fig. 1. The EpIC generator features a novel architecture based on the modularity programming paradigm. In this way, the code structure is kept as simple as possible and the addition of new developments, such as new channels or algorithms to generate random numbers, is made as easy as possible. The architecture of EpIC, including the nomenclature for naming its elements, is based on the PARTONS framework [20]. PARTONS is also used to evaluate the Born cross-section for a given process, which is used to generate MC events after the inclusion of radiative corrections (RCs). EpIC can generate events for a number of exclusive processes. The following are available at the time of this writing: deeply virtual Compton scattering (DVCS), time-like Compton scattering (TCS), and deeply virtual meson production (DVMP) of π 0 and π + mesons. In this article, we choose the case of DVCS whenever an exemplary process is needed. Adapting to other processes is straightforward and typically only requires replacing the acronym, e.g. the equivalence of DVCSGeneratorService for TCS is TCSGeneratorService. EpIC, as in the case of PAR-TONS, follows the so-called Trento convention [21] for the definition of e.g. scattering angles.
The purpose of this article is to provide a reference to be used for a better understanding of EpIC's structure and its functioning. We will mostly focus on the explanation of the architecture, starting with a brief description of the PARTONS framework in Sect. 2. Topics covered in that section are crucial to understand the technical aspects of EpIC. In Sect. 3, we present the architecture of EpIC. In particular, we list all types of modules and show the relation between them. The user interface is introduced in Sect. 4, while the demonstra-tion of the generator's performance is shown in Sect. 5. The summary is given in Sect. 6. In the appendix of the article (Appendix A), we review the basics of radiative corrections and explain how to compute them in the collinear approximation. The purpose of this appendix is to help users familiarise themselves with the parameters of RC used in EpIC.
EpIC is written in the C++ programming language. The code can be accessed at the GitHub platform [22] and is distributed under the GPL 3.0 licence. The project's webpage, Ref.
[23], includes information such as a technical documentation of C++ classes, a guide on how to compile the code, and example code for users.

PARTONS framework
PARTONS (PARtonic Tomography Of Nucleon Software) [20] is a software framework for studying the 3D structure of hadrons. It provides essential tools for QCD phenomenology and allows one to compute observables from the models of non-perturbative objects. For instance, it includes various models of GPDs and methods to perform their pQCD evolution. It also includes methods to evaluate amplitudes and cross-sections for a variety of exclusive processes, like DVCS and TCS, described at various orders of pQCD precision. The framework can be used in analyses to predict observables from existing models, see e.g. Refs. [15,16,24], but also to constrain new models from available experimental data [25,26]. Because of its versatility, the framework can also act as a laboratory for studying new concepts for phenomenology, like new ways of modelling, see e.g. Refs. [27][28][29].
What distinguishes PARTONS from other software projects used in the field of particle physics is its architecture utilising the modular programming paradigm. The central role in this architecture, which has been also adopted in EpIC, is played by the modules. A module is a single encapsulated development of a given type, for instance a single GPD model. The only purpose of modules is to provide data for a specific input, like values of GPDs for a given GPD kinematics. The construction of modules utilises both the inheritance and polymorphism mechanisms of C++. Therefore, the development of new modules is remarkably straightforward, as developers need only to provide their code related to the developments they are interested in and adapt it to the predefined classes in the framework. We demonstrate this with the example of GPDGK11 module, which is used to implement the Goloskokov-Kroll (GK) GPD model [30][31][32], from the PARTONS library. The following is an excerpt from the header file: where only the most relevant members of the class are shown, some of which will be discussed further in this section. The purpose of the GPDGK11 module is to evaluate the values of GPDs. This task is predefined in GPDModule, being in terms of inheritance the parent class of GPDGK11. The predefinition means here the declaration of a virtual function, so that each derived classes (GPDGK11, GPDGK16, . . . ) have a function of the same signature among its members, but the actual implementation of this function in those classes is in general different. The evaluation of GPDs takes place, for instance, in the following function: using namespace PARTONS; The function GPDGK11::computeH() returns a container that stores the numerical values of GPDs of type H, including singlet, H (+) , and non-singlet, H (−) , combinations, as defined in Ref. [13]. The input GPD kinematics is accessible in the body of the function via the m_x, m_xi, m_t, m_MuF2 and m_MuR2 variables, which are defined in the GPDModule class. Those variables correspond to x (average longitudinal momentum fraction carried by the active parton), ξ (skewness parameter), t (four-momentum transfer to the hadron target), µ 2 F (factorisation scale squared) and µ 2 R (renormalisation scale squared) GPDs depend on. They are set when one evaluates GPDs via PartonDistribution compute(const GPDKinematic& kinematic, GPDType::Type gpdType), which is also defined in GPDModule. As one may see, in our example quark flavours and GPD types are distinguished by QuarkFlavor and GPDType enums, respectively, while inputs and outputs are encapsulated in the QuarkDistribution, PartonDistribution and GPDKinematic containers.
A single instance of each module is loaded into memory at the beginning of the program execution, namely, when static constant variables are initialised. The addresses of such instances are stored by the registry, which later acts as a phone book, allowing one to access the instances by either their IDs or defined names. We demonstrate this mechanism with the following example: using namespace PARTONS; const unsigned int GPDGK11::classId = BaseObjectRegistry::getInstance()-> registerBaseObject(new GPDGK11("GPDGK11")); → → When the program starts, GPDGK11::classId is initialised with a unique ID assigned by the registry (BaseObjectRegistry class being a singleton). During this process, the registry takes and stores the address of each new instance. All modules inherit from the same basic class called BaseObject, allowing the registry to store addresses of the same type of objects. The constructor of BaseObject requires a unique string of characters that is used as a human-readable name of the module.
The registry stores addresses of single instances of pre-configured modules. Users may get copies of those instances by using the factory (also being a singleton), for instance: The registry and factory are the key architectural ingredients, which allow users to add an unlimited number of new modules in the PARTONS architecture, without having to modify existing modules.
There are many types of modules, but one may configure all of them, i.e. pass additional information to the modules, in a generic way, avoiding, for instance, the static casting. This is possible thanks to the void configure(const ElemUtils::Parameters& parameters) virtual functions, which developers may use to allow transferring additional parameters to modules. For instance: using namespace PARTONS; void configure(const ElemUtils::Parameters& can be used to set myVariable variable to 1 in the following way: //Container for multiple parameters. ElemUtils::Parameters parameters; //Set. parameters.add(ElemUtils::Parameter("myKey", 1)); where "myKey" is the unique name of the parameter. The way of supplying modules with additional parameters presented here plays a crucial role in constructing user interface (UI).
To avoid "manually" repeating some tasks, the services are used. In this way, users can perform complex tasks in a straightforward and robust manner by hiding the complexity of low-level functions. The services link such tasks as parsing input files, running many computations, employing multi-threading computing, filling the database, printing the output to the screen etc. Without services, all those tasks would have to be done explicitly by users, requiring a profound knowledge of the functioning of the project. This way of working would also be time consuming and vulnerable to mistakes.

EpIC's architecture
The generic architecture of the generator, including the flow of data transferred between the modules, is shown in Fig. 2. Each module takes input data, processes it according to predefined tasks, and returns the corresponding output. The actual implementation of some types of modules depends on the simulated process. For instance, RCModule is the parent class for DVCSRCModule and TCSRCModule, which are used for DVCS and TCS processes, respectively. These two types of modules process data according to the same task predefined in RCModule, but they use different inputs. It is either DVCS or TCS kinematics. Because of this difference, RCModule has been developed as a C++ template, allowing DVCSRCModule and TCSRCModule to handle input containers of various types. With those templates, the addition of new processes to the generator is as easy as possible. DVCSRCModule and TCSRCModule are used as base classes for actual modules, where specific physics developments are encoded. In Fig. 3 we present the entire tree of inheritance for modules of RCModule type.
Below, we provide a brief description of each type of module and its purpose. Additionally, we indicate whether the actual implementation of modules of a given type depends on the simulated process. This information is useful in recognising which parts of the generator need to be extended when a new process is added to the generator.

RCModule
Purpose: Simulation of radiative corrections.
Does it depend on the generated process? Yes.
Description: In EpIC events are generated according to the probability distribution given by: Here, dσ 0 (X) is the Born cross-section (given by PAR-TONS), while R(Z) is the radiator function. The set of variables that the Born cross-section depends on is denoted by X, for instance X = {x Bj , t, Q 2 , φ, φ S , E l } for DVCS, where x Bj and Q 2 are the usual DIS variables, φ is the angle between the lepton scattering plane and the production plane, φ S is the angle between the lepton scattering plane and the target spin component perpendicular to the direction of the virtual photon, and  E l is the lepton energy in the fixed target frame [20]. Because of the radiative corrections, the values of those variables may be different from those for a generated event, X . The radiator function, on the other hand, depends on the set of variables Z. The number of variables making this set depends on the approximation one uses. For instance, if one considers the collinear approximation and only the initial and final state radiations of single photons from the incoming and scattered leptons, two variables make the set Z, namely Z = {z 1 , z 3 }, see Appendix A for more details.
Taking the above in mind, modules of RCModule type perform three actions: i ) they define a number of variables that the radiative corrections depend on, i.e. they define the set Z, ii ) they implement a function used for the evaluation of R(Z), but also for the evaluation of "true" kinematics entering the Born crosssection, X = f (X , Z), iii ) they implement a function that stores the four-momenta of radiated photons in the final event records.

ProcessModule
Purpose: Evaluation of Born cross-section.
Does it depend on the generated process? Yes.
Description: The Born cross-section is evaluated in EpIC using the modules of ProcessModule type, which are defined in the PARTONS library. In the most general case, where the Born cross-section is directly calculated from GPDs, ProcessModule requires modules of other types for proper functioning: ConvolCoeffFunctionModule for the evaluation of the process amplitudes, GPDModule for the evaluation of GPDs at a reference scale, GPDEvolutionModule for the evolution of GPDs, XiConverterModule for the evaluation of GPD skewness variable, ξ, from the process kinematics, and ScalesModule for the evaluation of renormalisation and factorisation scales from the process kinematics. However, this general case can hardly be used in the Monte Carlo generation because of the nested integrations evaluating the Born cross-section. Instead, one may start the evaluation, for instance, from the level of amplitudes parametrised in lookup tables. In this case, GPDModule and GPDEvolutionModule are not used as there is no need for a convolution of GPDs with hard scattering coefficient functions. More details about the ProcessModule modules can be found in Ref. [20].  Fig. 3 Inheritance tree of modules used for the implementation of radiative corrections. Each box represents a single C++ class. Two first classes are the part of PARTONS library, the rest are the part of EpIC project. Classes acting as modules to be used be users are indicated by boldface fonts.

EventGeneratorModule
Purpose: Generation of kinematic configurations according to the probability distribution given by the product of the Born cross-section and the radiator function.
Does it depend on the generated process? No.
Description: Modules of this type play a central role in the generation of events. Three actions are implemented in those modules: i ) Initialisation, where the algorithm implemented in the module probes the probability distribution, for instance to "memorise" it using an interpolation method, or to capture its key features, like the location of places in which the probability distribution tends to change rapidly. For some complicated probability distribution functions, the initialisation process can be time-consuming, but the results of the process can be stored in a file for use in subsequent jobs.
ii ) Generation of kinematic configurations according to the probability distribution, based on the data collected during the initialisation stage. iii ) Evaluation of the total cross-section, which is needed to normalise the event distributions to a given integrated luminosity.

KinematicModule
Purpose: Evaluation of four-momenta for a given kinematic configuration.
Does it depend on the generated process? Yes.
Description: Modules of this type are responsible for the evaluation of four-momenta for a given kinematic configuration. In addition, they implement a predefined function that determines if a given kinematic configuration is physical, i.e. if it does not break any kinematic limit. This function is utilised by GeneratorModule to probe only valid kinematics.

WriterModule
Purpose: Composing and saving event records.
Does it depend on the generated process? No.
Description: The purpose of modules of this type is to create and store event records in output files. Different formats may be implemented by different modules. Additional information, like the integrated cross-section, can also be stored in output files, e.g. in their headers.

User interface
Both compilation and linking of the project make use of the CMake tool [33]. In the current version (1.0.0) EpIC requires the following external libraries: PAR-TONS [20] (providing elements of the architecture and used for the evaluation of Born cross-sections), ROOT [34] (used in one of the EventGeneratorModule modules), HepMC3 [35] (used in one of the WriterModule modules), GSL [36] (used for the generation of random numbers). For more detailed and always up-to-date information we refer to the online documentation [23]. The executable of the project, epic, must be invoked with two arguments, like this: where SEED is the random seed (unsigned integer) to be used in the initialisation of modules that deal with random numbers, and SCENARIO PATH is the relative or absolute path to the scenario containing all options used in the generation.
The EpIC scenarios are written in XML markup language which, thanks to the usage of tags based on standard words, is an easy-to-read format for humans while being straightforward for machines to process. EpIC assumes input data to be provided in units of GeV (and its powers, whenever applicable), and radians for angles. The general structure of a single scenario is the following:

Demonstration of EpIC
Here, we demonstrate EpIC in a detailed manner by using specific components of the PARTONS framework. The demonstration is based on one million MC events generated for the DVCS sub-process, i.e. exclusive leptoproduction of a single photon without the Bethe-Heitler contribution, at 10 GeV positive helicity electron and 100 GeV unpolarized proton energies. The following modules are used during the generation process: DVCSCFFCMILOU3DTables for the parameterisation of CFFs obtained from the GK GPD model [30][31][32] and LO coefficients functions, DVCSProcessBMJ12 for the evaluation of DVCS crosssection based on the set of expressions published in Ref. [38], DVCSRCNull for the simulation without radiative corrections, EventGeneratorFOAM for the generation of DVCS kinematics based on cross-section with the FOAM algorithm [39], DVCSKinematicDefault for the default evaluation of four-momenta from DVCS kinematics, and WriterHepMC3 to save four-momenta in HepMC3 format.
The following cuts are used in the generation process: 0.0001 ≤ x B ≤ 0.6, 0.01 ≤ y ≤ 0.95 (here, y is the inelasticity parameter), 1 ≤ Q 2 ≤ 100 GeV 2 , 0 ≤ |t| ≤ 1 GeV 2 , 0 ≤ φ ≤ 2π, and 0 ≤ φ S ≤ 2π. On the basis of these cuts, the distribution of events is shown in Fig. 4. The following FOAM library parameters are used in producing those events [39]: nCells = 3000, nSamples = 600, nBins = 600. For the FOAM library, the initialisation lasted about 40 minutes, whereas the generation time after initialisation took around 0.0052 seconds per event at BNL computer farms. The initialisation time can be shortened by changing the param-eters of the FOAM library, which may, however, lead to a longer time needed by EpIC to generate a single event. We remind that the result of the initialisation can be stored in a file, allowing to avoid this step in subsequent MC runs.
To ensure our results are consistent, we also plot the theory expectation values, as solid lines, on top of each histogram. The number of events associated with the theory curves is calculated as follows: at each bin, the cross-section for the given bin width is calculated, and then the result is divided by the total cross-section of the phase space. One can then multiply this ratio by the total number of events to find the number of events in that specific bin. To obtain those curves in Fig. 4, first we assign the number of events to the middle point of each bin, then interpolate all these results. As a result, we observe that EpIC generates events in accordance with the theory values.

Summary
In this work, we described in detail the overall structure and essential features of the novel MC generator EpIC. This generator is based on the PARTONS framework in the model selection. It is designed to generate events for exclusive processes, such as DVCS, TCS, DVMP and more broadly to all exclusive processes implemented in the PARTONS framework. As these processes have the potential to reveal the 3D structure of the nucleon through GPDs, they have been studied extensively at facilities such as DESY, JLab, and CERN, and they will remain crucial for future EIC programmes.
The EpIC MC generator stands out with its flexible architecture that utilises a modular programming paradigm. Consequently, the code is easier to understand, as each element of the architecture is designed to perform a specific task. The distribution of events generated by EpIC is highly precise since it utilises an accurate representation of the cross-section of the underlying process in the FOAM library. Besides generating accurate events, EpIC also allows users to generate events with radiative corrections as well as to run multi-channel analyses. On top of this, EpIC is released through GitHub as an open-source project under a GPL licence. In particular it means that users can trace code modifications and new contributors can easily join the development team. To the best of our knowledge, a tool such as EpIC is unique in the GPD community by its very nature. The comprehensive set of features and the variety of models available in EpIC make it ideal for the systematic study of the 3D structure of the nucleon.
The generator was designed in a way that facilitates easy expansion. The encapsulation of elements is one of the fundamental characteristics of the generator when it comes to maintaining the project for a long time and with a changing team of developers. This allows us to easily adapt the architecture to the latest developments. Finally, to make this project as useful to the particle physics community as possible, we promote open-source standards and offer easy access to its development. This perspective is well suited to the experimental timeline of the 3D hadron structure community. good accuracy in the collinear approximation [40,41]. In this approximation the transverse component of the momenta of emitted photons are neglected. Therefore, the radiated photons are restricted to move along the direction of the source leptons.
The collinear approximation allows one to introduce a single parameter, z 1 or z 3 , for the initial or final state radiation, respectively. The parameter describes the energy of the emitted photon: where E e (E e ) and E γ (E γ ) are the energies of the incoming (outgoing) lepton and the initial (final) state emitted photon, respectively. The DIS cross section can be then expressed by [41]: where d 2σ Born is the differential Born cross-section evaluated for "true" variables, which are related to the observed ones by [41]: and where due to kinematic limits one has: with L = ln(Q 2 /m 2 l ). Here, α is the fine-structure constant, and m l is the lepton mass. The parameter acts as a cutoff avoiding the generation of a very soft photons characterised by energies E γ ≤ E e and E γ ≤ E e . It is crucial to choose the cutoff carefully so that E e and E e are much smaller than the energy resolution of the experiment. Due to the cutoff Eq. (A.5) consists of two parts: the term proportional to the Dirac delta function collectively describes the contribution of very soft photons and virtual corrections, while the term proportional to the step function describes the effects of real radiation of photons with sizeable energies.
The generalisation of leading order radiative corrections for the DVCS case can be achieved as follows: