CepGen -- A generic central exclusive processes event generator for hadron-hadron collisions

We present an event generator for the simulation of central exclusive processes in hadron-hadron reactions. Among others, it implements the two-photon production of lepton pairs previously introduced in LPAIR. As a proof of principle, we show that the two approaches are numerically consistent. The $k_{\rm T}$-factorized description of this process is also handled, along with the two-photon production of a quark, or a $W^\pm$ gauge boson pair. This toolbox may be used as a common framework for the definition of many other processes following this approach. Additionally, photoproduction and other photon induced processes are also considered, or being implemented.


I. INTRODUCTION
We present a new event simulation tool to facilitate the study of central exclusive processes (CEP) in the scope of high-energy colliding experiments.
Lately, with the experimental observations of photoninduced processes at TeV scale energies with or without the help of forward proton taggers, the need for such a common simulation framework that would allow both the phenomenology and experimental studies of this specific class of events, and easily manageable by the two communities, has been clearly motivated.
In particular, the two-photon production of lepton pairs in proton-proton reactions through a t-channel exchange will be covered here. This choice may be justified by the well demonstrated theoretical interest of this process, namely the precision attached to its purely elastic predictions (a pure tree-level electrodynamic component), and large uncertainties to inelastic contributions. The latter are closely related to the proton's electromagnetic and nuclear structure modeling, thus enabling its precision study in a hadron-hadron collider environment.
Unlike the vast majority of common CEP event generators implementing the incoming photon fluxes through the equivalent photon approximation (EPA) [1], the full unintegrated photon virtuality (both the transverse and longitudinal components) information is used here. This allows a refined treatment of processes in which lowvirtualities are accounting for the biggest fraction of the total cross section, such as the two-photon process quoted above.
We will therefore concentrate on this particular photon-induced process as a proof of principle for the operation and validation of this simulation tool. This paper is organized as follows: in Section II we enumerate and describe some of the processes currently * laurent.forthomme@cern.ch; Now at European Organization for Nuclear Research (CERN), Geneva, Switzerland.
implemented in this code. Then, in Section III we present a review of the technicalities carried in its implementation. In Section IV we list the extensions currently embedded in CepGen to ease the user interaction and allow an increased modularity with its core features. Finally, our conclusions are summarized in Section V.

II. PHYSICS CASES
Beside the two-photon (t-channel) production of lepton (or fermion) pairs already mentioned in this paper, multiple matrix elements are included in the current version of our event generator. For instance, one may quote: • the two-photon production of a W ± gauge boson pair, recently studied experimentally at the LHC [2][3][4], and formulated in [5] ; • the diffractive photoproduction of low-and intermediate mass vector mesons, previously implemented in DIFFVM [6] ; • the inclusive, gluon-induced production of a fermion pair, using unintegrated parton distributions following the KMR procedure [7,8].
A detailed review of the last two processes will appear in a future paper.
If we again concentrate on the two-photon production of a fermion pair, two approaches are implemented within CepGen. For the dilepton final state, one may find a modern retranscription of the LPAIR code [9,10] used intensively in various HERA searches, and introduced in Section II A. Along with the two-photon production of W ± pairs, it can also be formulated through the k Tfactorization technique described in [11,12].
The elastic photon emission through this k T -factorized process may also be formulated for heavy-ion initial beams. In CepGen, this initial state is still being validated while this paper is released.

A. LPAIR matrix element
The photon emission from initial state interacting protons may be divided in two main sub-components. For the single-, double-dissociative (or semi-exclusive, nonexclusive) scenarios of the two-photon production of leptons pair, the EPA is not sufficient to reproduce the high decorrelation of the two individual leptons observed in the elastic case.
In the scope of an ep collider such as HERA, the LPAIR event generator, with its full two-to-four matrix element definition and its implementation of the proton structure functions modeling, enabled the probe of this additional final state in which the proton would scatter a photon with a sufficient virtuality to excite and fragment in its final state. The strong benefit of LPAIR was its strong numerical treatment of the phase space, ensuring a good stability of its matrix element despite the low-Q 2 range of the incoming state photons accounting for a large fraction of the total cross section. This feature was particularly relevant in the elastic case.
The version implemented in CepGen handles by de-fault the pp initial state. It can however be steered to retrieve the original ep case. Originally interfaced to the JETSET library, the diffractive/dissociated proton(s) could be studied in terms of its/their decay products in LPAIR. Based on this approach, a refreshed implementation of this fragmentation feature has been ported to CepGen, as discussed in Section IV A.
B. kT-factorized matrix elements As previously introduced in [11] for two-photon processes, the k T factorization approach allows to treat separately the hard process definition and the modeling of incoming photon fluxes. Unlike the collinear approximation, photons densities are left unintegrated and their transverse momentum contribution can be accounted for in the central kinematic variables description.
The differential cross section of a generic pp → p ( * ) (γγ → X)p ( * ) (here, p ( * ) conventionally expresses both the final state protonic system after and elastic or dissociative emission of the photon) can therefore be expressed as follows: with F γ el,inel the unintegrated incoming photon fluxes, and dσ the differential cross section for the transverse momentum dependent two-photon process. The earlier may be expressed as a function of the fractional proton momentum ξ carried by the scattered photon, and its squared transverse momentum norm k 2 ⊥ . For instance, using the Budnev prescription quoted in (1), one may define the following modeling for the low-Q 2 range elastic scattering where the proton remains on-shell after photon emission: where Q 2 el , the photon virtuality for elastic emission is obtained from the fractional proton momentum loss and the transverse component photon virtuality: In the formalism of (2), we use: as linear combinations of G E,M (Q 2 ), the electric and magnetic form factors of the proton.
As for the unintegrated inelastic scattering flux, where the initial proton loses a sufficiently high virtuality to be excited into a dissociative system of mass M X , it can be expressed as: with the photon virtuality after its inelastic emission from the beam particle, and x Bj = Q 2 inel /(Q 2 inel + M 2 X − m 2 p ) the Bjorken scaling variable. This inelastic density is hence modeling-dependent through the definition of its F 2,L (x Bj , Q 2 ) proton structure functions, along with their linear combination F 1 (x Bj , Q 2 ). The latter is conventionally defined as: A non-exhaustive list of parametrization implemented in this simulation tool can be found in Section III B.
For completeness, the EPA photon fluxes in the collinear approximation commonly used in other simulation tools for photon-induced processes can be obtained by integrating the Budnev fluxes defined in (2) and (3) over the full transverse virtuality range of interest, i.e.
More generally, it is up to the process developer to add its own modeling of this k T -parameterized flux, including for other intermediate particles exchanges. The earlier photon fluxes definitions are however provided as a core features of CepGen.

III. IMPLEMENTATION
For the sake of simplicity, this event simulation tool may be factorized into several independent building blocks. Its central part is the implementation of the matrix element, as a method computing the scalar event weight for each point computed in the allowed kinematic phase space. Either the full incoming state is needed as an input to this method, as requested by the original LPAIR code, or the already simplified set of kinematic information needed in the k T -factorization scheme.
Any process may be defined by the library developer as a class derived from a pure virtual, generic process object provided by the generator core. It is required to override the three following members: • a tool associating the size of the phase space for each of the modes considered (e.g. elastic, singledissociative, double-dissociative final states for pp processes), and propagated up to the integration algorithm ; • a method computing and returning a scalar event weight for any physical point in the phase space, and null elsewhere ; • an event definition function associating the kinematic variables of all particles in the event for a given point.
This latter populates an internal event object already containing the full incoming state with the central and forward outgoing system for each point. The user may hence retrieve this object at any stage of the computation process and/or apply a more advanced set of physical selections to fit e.g. a given set of experimental constraints.
To ease the integration of new k T -factorized process, an additional derivative object provides the full incoming partons (photons, gluons, . . . ) kinematic variables as an input. Therefore, with the help of a toolbox provided for the definition of inner variables and the Jacobian computation, the processes developer only needs to define a central factorized matrix element and the event structure to allow the process to be fully handled by CepGen.
Furthermore, a FORTRAN interfacing module 1 has been developed to ease the definition of any central k Tfactorized process for a larger fraction of the community.

A. Usage and configuration
Each cross section computation or event generation can be parameterized through a configuration file to be parsed and fed to the internal parameters. Currently, two formats may be parsed by two handler objects: • the standard LPAIR steering cards format, combining a four-letters keyword with a configurable value (either an integer, a floating point value, or a characters string). It is designed to be fully compatible with all cards previously issued and used in the original version of LPAIR [10] ; • a structured Python configuration file, defining all parameters collections and their normal usage.
With the native Python grammar and syntax verification, it allows an easier interaction with the end-user and process developer. Additionally, its nested structure enables an easy import of external configuration parameters while reducing the cards verbosity.
An example of a comparison between the two steering cards formats for the same run requirements may be found in the supporting website (see Appendix).

B. Proton form factors and structure functions
As already introduced, in central exclusive processes, the elastic and semi-elastic scattering of an intermediate particle from the incoming protons may be modeled through its electromagnetic form factors. Our implementation uses the dipole approximation of the Q 2 dependence of the electric and magnetic components introduced in Section II B.
This elastic scattering of photons from the proton is usually characterized by its very low momentum transfer (or virtuality). In particular, a dominant fraction of their momentum is emitted collinear to the initial beam.
With increasing photon virtualities, the probability of protons to leave their bound state to fragment into subcomponents increases as well. As discussed, this behavior is mostly dominated by the inner structure of the proton, characterized by its structure functions.
In the high-energy perturbative limit, these structure functions may be defined according to inner partons distributions in a simple leading-order quark-parton model: with e f the quark electric charge (in units of e), and q (q) the parton density function of the quark (anti-quark) flavor. An interface between CepGen and the LHAPDF [13] library is therefore provided for its vast choice of modelings and its supporting tool for the numerical evaluation of these collinear densities.
Additionally, the following parametrizations of F 2,L structure functions covering a wide range of fractional momentum losses and photon virtualities are implemented in CepGen: • the Suri-Yennie [14] parametrization, extracted from the experimental observations at low-Q 2 of the total γp cross section in several resonances regions, with invariant masses extending from 1.11 to 18.03 GeV ; • the Fiore et al. [15] modeling of the very low Q 2 region from the JLab and SLAC observations of res-onances, already implemented in LPAIR and reimplemented in CepGen for historical reason ; • the Christy-Bosted [16] low-Q 2 (0 < Q 2 < 10 GeV 2 ) and low-diffractive mass (1.1 < M X < 3.2 GeV) modeling, also tuned around the low-mass resonances region ; • the Szczurek-Uleshchenko parametrization [17], valid at low and intermediate Q 2 (up to 100 GeV 2 ), and intermediate-x values. It relies on a shifted value of the factorization scale ensuring the structure functions convergence at Q 2 ≈ 0 ; • the Abramowicz-Levin-Levy-Maor (ALLM) parametrization [18,19] of the intermediate, continuum regime. Updated fits containing additional HERMES data samples (GD07p/GD11p) [20,21] are also provided ; • a F 2/L grid 2 built from the MSTW partonic density functions [22] evaluated at next-to-next-to-leading order.
• a hybrid "LUXlike" set of structure functions with input from the Christy-Bosted resonance, GD11p/ALLM97 continuum, and any perturbative set described earlier. This composite modeling, described in details in [5], is based on the prescription of [23]. A (x Bj , Q 2 ) mapping of the both structures functions for this set is pictured in Fig. 1.
A graphical summary of a good fraction of F 2 parametrizations implemented in CepGen can be found in Fig. 2 for several Q 2 momentum scales, along with a subset of experimental observations used in the global fits. While a fraction of parametrizations is only maintained for historical reasons and ensure backward compatibility with major LPAIR versions, the Suri-Yennie and LUXlike implementations are considered standard for the two-photon production of lepton, gauge boson pairs respectively.
Depending of the modeling, the longitudinal structure function F L can either be extracted from the overall fit, or derived from F 2 and the longitudinal-to-transverse cross sections ratio R(x Bj , Q 2 ) = σ L /σ T through the relation: Again, a collection of data-driven modelings of this ratio is provided natively within the CepGen structure functions library. In this first version of the code, the following are implemented: • the E143 fit [26], valid for small-x Bj values (0.03 < x Bj < 0.1) ; • the R1990 fit [27], extracted in the SLAC energy range (0.6 < Q 2 < 20 GeV 2 , 0.1 < x Bj < 0.9), and known to give its best predictions for Q 2 > 0.3 GeV 2 ; • the Sibirtsev-Blunden parametrization [28], combining measurements from JLab Hall C and SLAC experiments at intermediate Q 2 .
The full set of structure functions (both F 2 and F L , but also the linear combination F 1 ) handled in CepGen can be used in any external user application, being selfcontained in a standalone library. Furthermore, a FOR-TRAN interface is also provided to ease the user interaction with this collection.

C. Matrix element integration
The integration over the full phase space definition is performed using the GSL [29] implementation of the plain, Vegas [30], and MISER [31] Monte-Carlo integration algorithms. The first simply probes randomly the full phase space in a fast but inefficient way, while the latter two rely on an optimized choice of the integration grid.
In Vegas, a first probe allows to extract a grid handling the information on all potential numerical singularities. This operation is critical to ensure the numerical stability of the whole procedure. In the particular case of twophoton induced processes, several peaking distributions are usually expected (for instance, the photon virtuality Q 2 at low values in the case of t-channel exchanges).
In a second iteration, the integrand is integrated over the whole phase space, following a generation density determined by this first preparation stage. The resulting total cross section can hence be estimated within the level of a few percents of precision within only a few iterations. By default in CepGen, the iteration process is stopped once the χ 2 value is compatible with unity.
In the MISER algorithm, the recursive stratified sampling allows to optimize the phase space point generation through a prior estimate of the higher variance regions to be specifically probed. A precise estimate of the overall cross section with a lower multiplicity of points probed can therefore rapidly be obtained.

D. Events generation
The generation of unweighted events to be used as observable states in a simulated process is performed through the hit or miss Monte Carlo technique.
For the Vegas integration algorithm, the gridpreparation stage may be used again to optimize the unweighting of events through importance sampling, and reduce the per-phase space point generation time. This procedure, known as integrand treatment may be steered directly by the user.  [21,24,25]. A tolerance of 1% is set to the Q 2 value to account for slight differences in experimental kinematic ranges.
Since its version 1.0, CepGen introduces an experimental multi-thread implementation of the event generation component. With this approach, the reduced information sharing between each single thread is allowing for them to build unweighted events in a fast optimized way. It allows a noticeable reduction of the per-event generation time of most complex processes currently supported.

E. Validation with LPAIR
As a validation of the CepGen implementation of the LPAIR γγ → + − process introduced in Section II A, the generator level cross section can be evaluated for both generators for multiple definitions of the phase space. In Fig. 3, the scan of this production cross section as a function of the lower transverse momentum cut applied on both the outgoing leptons is pictured in the three possible proton final states for the original LPAIR algorithm (solid lines) and for this implementation (open markers).
The slight drift in the pull distribution observed at higher values of the single lepton transverse momentum cut are dominantly due to the numerical definition of several physics constants in LPAIR, and corrected in this code. However, a good agreement can be seen between the two implementations of this matrix element.
In Fig. 4, both the CepGen and the LPAIR differential cross section distributions for the dilepton invariant mass, transverse momentum, and the normalized azimuthal correlation between the two leptons are shown, along with the ratio to the LPAIR predictions. The error bars displayed for all variables account for the statistical uncertainties for the generation of 10 6 events with Cep-Gen.
Several examples of event generation and cross section computation are provided with the core library to guide the end user towards its integration within the environment. Multiple wrappers provided along with this software suite allow the direct conversion of the internal event and particles structure to a growing multiplicity of event formats commonly used in high-energy physics 3 .

IV. ADDITIONAL FEATURES
Thanks to its modularity, the core features of CepGen can be extended with a set of plugins to ease the end user interaction with external utilities and detector simulation algorithms.

A. Resonances decays and excited proton fragmentation
For this initial release, the PYTHIA 8 [35] library is interfaced to CepGen to allow its branching fractions, decay products, and decay algorithms implementations to decay all unstable particles. An additional interfacing module for Herwig 7 [36,37] is currently being developed to be released for a later version.
The steering of these external libraries can be performed through the Python cards objects definition.
As observed earlier, the large exchanged photon virtualities may also induce the excitation of the scattered proton leading to its dissociation. In CepGen, the inelastic photon emission is parameterized through the splitting of this proton into a quark-diquark system. There, the excited outgoing state is defining the kinematic variables of a valence quark (carrying a fraction x Bj of the outgoing diffractive proton total momentum), thus leaving behind a corresponding color-connected diquark remnant. This 3 A few output formats handled in the current version of CepGen are: HepMC ASCII [32], and Les Houches Event record [33]. An example executable also produces simple tree-structured ROOT files [34] broadly used in the experimental physics community.    latter is hadronized through the Lund fragmentation algorithm implemented in PYTHIA.
Given its overall knowledge of the full event topology, the multi-parton interactions (MPI), or initial-and finalstate radiations frameworks of this latter may be steered to account for larger rapidity gaps suppression in the central detector. It may hence be adjusted to the drop in survival probabilities observed experimentally in many recent searches [2][3][4]38], and lately studied in [39].

B. Taming functions
For an increased modularity, the differential matrix element can furthermore be modified through a set of taming functions steered by the end user. It allows to interact directly with the process kinematic variables and account for additional physics effects in the computation of both the differential and integrated cross sections, while leaving behind details of initial parton or hard matrix element modification.
The differential cross section at a given phase space point x may hence be modified as: with {w i } the collection of analytical taming functions to be applied on the resulting matrix element.
Among the possible effects to simulate, one may quote the soft survival probability corrections resulting from large rapidity gap suppression, and observed experimentally in the high-energy limit of central exclusive productions.
This approach has been studied for instance in [40], where a concept of modified photon parton density functions accounting for large rapidity gap suppression is introduced. The result is an overall modification of the central event kinematic variables.
This taming functions framework hence provides a toolbox to evaluate the effects of this complex procedure to the phase space definition. For instance in the twophoton production of a lepton pair, the rapidity gap survival probability may be treated as a dilepton transverse momentum-dependent correction to follow an exponential reduction, such as: . For the single-dissociative contribution, it corresponds to a survival factor -after its integration over the full phase space considered-of S 2 = 0.76. The effects of this suppression can be studied in terms of other experimental observables, as pictured in Fig. 5.

V. CONCLUSIONS
A new framework for the phenomenological and experimental studies of central exclusive processes in the    scope of hadron-hadron colliders is presented. This tool allows both the communities to develop and parameterize any photon-or color singlet-induced process in terms of simple and highly configurable building blocks.
Some of these blocks, as for instance the large collection of structure functions modelings, may also be included individually in any external user library.
Furthermore, both the C++ and FORTRAN interfaces to processes definition also allows a large fraction of model builders to take part to the development and the study of additional final states, leaving all technicalities of integration, events generation, and external libraries interfacing to be covered by this code.
With an ever increasing number of output formats supported, a direct interfacing to common simulation and reconstruction chains may easily be developed, thus allowing the portage of numerous standalone simulation tools in an embedded framework.