CASCADE3 A Monte Carlo event generator based on TMDs

The CASCADE3 Monte Carlo event generator based on Transverse Momentum Dependent (TMD) parton densities is described. Hard processes which are generated in collinear factorization with LO multileg or NLO parton level generators are extended by adding transverse momenta to the initial partons according to TMD densities and applying dedicated TMD parton showers and hadronization. Processes with off-shell kinematics within $k_t$-factorization, either internally implemented or from external packages via LHE files, can be processed for parton showering and hadronization. The initial state parton shower is tied to the TMD parton distribution, with all parameters fixed by the TMD distribution.


Introduction
The simulation of processes for high energy hadron colliders has been improved significantly in the past years by automation of next-to-leading order (NLO) calculations and matching of the hard processes to parton shower Monte Carlo event generators which also include a simulation of hadronization. Among those automated tools are the MADGRAPH5 aMC@NLO [1] generator based on the MC@NLO [2][3][4][5] method or the POWHEG [6,7] generator for the calculation of the hard process. The results from these packages are then combined with either the HERWIG [8] or PYTHIA [9] packages for parton showering and hadronization. Different jet multiplicities can be combined at the matrix element level and then merged with special procedures, like the MLM [10] or CKKW [11] merging for LO processes, the FxFx [12] or MiNLO method [13] for merging at NLO, among others. While the approaches of matching and merging matrix element calculations and parton showers are very successful, two ingredients important for high energy collisions are not (fully) treated: the matrix elements are calculated with collinear dynamics and the inclusion of initial state parton showers results in a net transverse momentum of the hard process; the special treatment of high energy effects (small x) is not included.
The CASCADE Monte Carlo event generator, developed originally for small x processes based on high-energy factorization [14] and the CCFM [15][16][17][18] evolution equation, has been extended to cover the full kinematic range (not only small x) by applying the Parton Branching (PB) method and the corresponding PB Transverse Momentum Dependent (TMD) parton densities [19,20]. The initial state evolution is fully described and determined by the TMD density, as it was in the case of the CCFM gluon density, but now available for all flavor species, including quarks, gluons and photons at small and large x and any scale µ. For a general overview of TMD parton densities, see Ref. [21].
With the advances in determination of PB TMDs [19,20], it is natural to develop a scheme, where the initial parton shower follows as close as possible the TMD parton density and where either collinear (on-shell) or k t -dependent (off-shell) hard process calculations can be included at LO or NLO. In order to be flexible and to use the latest developments in automated matrix element calculations of hard process at higher order in the strong coupling α s , events available in the Les Houches Event (LHE) file format [22], which contains all the information of the hard process including the color structure, can be further processed for parton shower and hadronization in CASCADE3.
In this report we describe the new developments in CASCADE3 for a full PB-TMD parton shower and the matching of TMD parton densities to collinear hard process calculations. We also mention features of the small-x mode of CASCADE3.

The hard process
The cross section for the scattering process of two hadrons A and B can be written in collinear factorization as a convolution of the partonic cross section of partons a and b, a + b → X, and the densities f a(b) (x, µ) of partons a (b) inside the hadrons A (B), where x a (x b ) are the fractions of the longitudinal momenta of hadrons A, B carried by the partons a(b), σ(a + b → X) is the partonic cross section, and µ is the factorization scale of the process. The final state Y contains the partonic final state X and the recoils from the parton evolution and hadron remnants. In CASCADE3 we extend collinear factorization to include transverse momenta in the initial state, either by adding a transverse momentum to an on-shell process or by using offshell processes directly, as described in detail in Sections 2.1 and 2.2. TMD factorization is proven for semi-inclusive deep-inelastic scattering, Drell-Yan production in hadron-hadron collisions and e + e − annihilation [23][24][25][26][27][28][29][30][31][32][33][34][35]. In the high-energy limit (small-x) k T -factorization has been formulated also in hadronic collisions for processes like heavy flavor or heavy boson (including Higgs) production [14,[36][37][38], with so-called unintegrated parton distribution functions (uPDFs), see e.g. Refs. [39][40][41][42][43][44][45][46][47][48][49].

On-shell processes
The hard processes in collinear factorization (with on-shell initial partons, without transverse momenta) can be calculated by standard automated methods like MAD-GRAPH5 aMC@NLO [1] for multileg processes at LO or NLO accuracy. The matrix element processes are calculated with collinear parton densities (PDF), as provided by LHAPDF [50].
We extend the factorization formula given in eq.(1) by replacing the collinear parton densities f (x, µ) by TMD densities A(x, k t , µ) with k t being the transverse momentum of the interacting parton, and integrating over the transverse momenta.
However, when the hard process is to be combined with a TMD parton density, as described later, the integral over k t of the TMD density must agree with the collinear (k tintegrated) density; this feature is guaranteed by construction for the PB-TMDs (also available as integrated PDFs in LHAPDF format).
In a LO partonic calculation the TMD or the parton shower can be included respecting energy momentum conservation, as described below. In an NLO calculation based on the MC@NLO method [2][3][4][5] the contribution from collinear and soft partons is subtracted, as this is added later with the parton shower. For the use with PB TMDs, the HERWIG6 subtraction terms are best suited as the angular ordering conditions coincide with those applied in the PB-method. The PB TMDs play the same role as a parton shower does, in the sense that a finite transverse momentum is created as a result of the parton evolution [51,52].
When transverse momenta of the initial partons from TMDs are to be included to the hard scattering process, which was originally calculated under the assumption of collinear initial partons, care has to be taken that energy and momentum are still conserved. When the initial state partons have transverse momenta, they also acquire virtual masses. The procedure adopted in CASCADE3 is the following: for each initial parton, a transverse momentum is assigned according to the TMD density, and the parton-parton system is boosted to its centerof-mass frame and rotated such that only the longitudinal and energy components are nonzero. The energy and longitudinal component of the initial momenta p a,b are recalculated taking into account the virtual masses Q 2 a = k 2 t,a and Q 2 b = k 2 t,b [53], withŝ = (p a +p b ) 2 with p a (p b ) being the four-momenta of the interacting partons a and b. The partonic system is then rotated and boosted back to the overall center-of-mass system of the colliding particles. By this procedure, the parton-parton mass √ŝ is exactly conserved, while the rapidity of the partonic system is approximately restored, depending on the transverse momenta.
In Fig. 1 a comparison of the Drell-Yan (DY) mass, transverse momentum and rapidity is shown for an NLO calculation of DY production in pp collisions at √ s = 13 TeV in the mass range 30 < m DY < 2000 GeV. The curve labelled NLO(LHE) is the calculation of MAD-GRAPH5 aMC@NLO with the subtraction terms, the curve NLO(LHE+TMD) is the prediction after the transverse momentum is included according to the procedure described above. In the p T spectrum one can clearly see the effect of including transverse momenta from the TMD distribution. The DY mass distribution is not changed, and the rapidity distribution is almost exactly reproduced, only at large rapidities small differences are observed. The transverse momenta k t are generated according to the TMD density A(x, k t , µ), at the original longitudinal momentum fraction x and the hard process scale µ. In a LO calculation, the full range of k t is available, but in an NLO calculation via the MC@NLO method a shower scale defines the boundary between parton shower and real emissions from the matrix element, limiting the transverse momentum k t . Technically the factorization scale µ is calculated within CASCADE3 (see parameter lhescale) as it is not directly accessible from the LHE file, while the shower scale is given by SCALUP. The limitation of the transverse momenta coming from the TMD distribution and TMD shower to be smaller than the shower scale SCALUP guarantees that the overlap with real emissions from the matrix element is minimized according to the subtraction of counterterms in the MC@NLO method.
The advantage of using TMDs for the complete process is that the kinematics are fixed, independent of simulating explicitly the radiation history from the parton shower. For inclusive processes, for example inclusive Drell-Yan processes, the details of the hadronic final state generated by a parton shower do not matter, and only the net effect of the transverse momentum distribution is essential. However, for processes which involve jets, the details of the parton shower become also important. The parton shower, as described below, follows very closely the transverse momentum distribution of the TMD and thus does not change any kinematic distribution after the transverse momentum of the initial partons are included.
All hard processes, which are available in MADGRAPH5 aMC@NLO can be used within CASCADE3. The treatment of multijet merging is described in Section 8.

Off-shell processes
In a region of phase space, where the longitudinal momentum fractions x become very small, the transverse momentum of the partons cannot be neglected and has to be included already at the matrix element level, leading to so-called off-shell processes.
In off-shell processes a natural suppression at large k t [54] (with k t > µ) is obtained, shown explicitly in Fig. 2, where the matrix element for g * g * → QQ, with Q being a heavy quark, is considered. The process is integrated over the final state phase space [55], where dLips is the Lorentz-invariant phase space of the final state, ME is the matrix-element for the process, φ 1,2 is the azimuthal angle between the two initial partons, and a simple scale-independent and k t -independent gluon density xG(x) = (1 − x) 5 is included which suppresses large-x contributions. In Fig. 2 we showσ(k t ) normalized to its on-shell valuẽ σ(0) at √ s = 13000 GeV as a function of the transverse momentum of the incoming gluon k t,2 for different values of x 1 , which are chosen such that the ratio k 2 t,1 /(x 1 s) is kept constant. In Fig. 2 (left) predictions are shown for bottom quarks with mass m = 5 GeV and different k t,1 , in Fig. 2 (right) a comparison is made for different heavy quark masses. Using off-shell matrix elements a suppression at large transverse momenta of the initial partons is obtained, depending on the heavy flavor mass and the transverse momentum. In a collinear approach, with implicit integration over transverse momenta of the initial state partons, the transverse momenta are limited by a theta function at the factorization scale, while off-shell matrix elements give a smooth transition to a high k t tail.
When using off-shell processes, BFKL or CCFM type parton densities should be used to cover the full available phase space in transverse momentum, which can lead to k t 's larger than the transverse momentum of any of the partons of the hard process [56]. Until now, only gluon densities obtained from CCFM [15][16][17][18] or BFKL [57][58][59] are available, thus limiting the advantages of using off-shell matrix elements to gluon induced processes.
Several processes with off-shell matrix elements are implemented in CASCADE3 as listed in Tab. 1, and described in detail in [60]. However, many more processes are accessible via the automated matrix element calculators for off-shell processes, KATIE [61] and PEGA-SUS [62]. The events from the hard process are then read with the CASCADE3 package via LHE files. For processes generated with KATIE or PEGASUS no further corrections need to be performed and the event can be directly passed to the showering procedure, described in the next section.

Initial State Parton Shower based on TMDs
The parton shower, which is described here, follows consistently the parton evolution of the TMDs. By this we mean that the splitting functions P ab , the order and the scale in α s as well as kinematic restrictions are identical to both the parton shower and the evolution of the parton densities (for NLO PB TMD densities, the NLO DGLAP splitting functions [73,74] together with NLO α s is applied, while for the LO TMD densities the corresponding LO splitting functions [75][76][77] and LO α s is used).

From PB TMD evolution to TMD Parton Shower
The PB method describes the TMD parton density as (cf eq.(2.43) in Ref. [19]) with z M < 1 defining resolvable branchings, k (q c ) being the transverse momentum vector of the propagating (emitted) parton, respectively. The transverse momentum of the parton before branching is defined as k t = |k + (1 − z)q| with q = q c /(1 − z) being the rescaled transverse momentum vector of the emitted parton (see Fig. 3, with the notation k t = |k| and q = |q|) and φ being the azimuthal angle between q and k. The argument in α s is in general a function of the evolution scale q. Higher order calculations indicate the transverse momentum of the emitted parton as the preferred scale. The real emission branching probability is denoted by P f (z, q)), z) including α s as described in Ref. [19] (in the following we omit α s in the argument of P (R) ab for easier reading). The Sudakov form factor is given by: Dividing Eq.(5) by ∆ a (µ 2 ) and differentiating with respect to µ 2 gives the differential form of the evolution equation describing the probability for resolving a parton with transverse momentum k and momentum fraction x/z into a parton with momentum fraction x and emitting another parton during a small decrease of µ, The normalized probability is then given by This equation can be integrated between µ 2 i−1 and µ 2 to give the no-branching probability (Sudakov form factor) for the backward evolution ∆ bw , 1 with x = x/z. This Sudakov form factor is very similar to the Sudakov form factor in ordinary parton shower approaches, with the difference that for the PB TMD shower the ratio is applied, which includes a dependence on k t . In Eq.(9) a relation between the Sudakov form factor ∆ a used in the evolution equation and the Sudakov form factor ∆ bw used for the backward evolution of the parton shower is made explicit. A similar relation was also studied in Refs. [78,79]. In Ref [78] the z M limit was identified as a source of systematic uncertainty when using conventional showers with standard collinear pdfs; in the PB approach, the same z M limit is present in the parton evolution as well as in the PB-shower. The PB approach allows a consistent formulation of the parton shower with the PB TMDs, as in both Sudakov form factors ∆ a and ∆ bw the same value of z M is used.
The splitting functions P (R) ab contain the coupling, where the scale f (z, q) in the coupling depends on the ordering condition as discussed later (see Eq. (11)). The advantage of using a PB TMD shower is that as long as the parameters of the parton shower are set through TMD distributions the parton shower uncertainties can be recast as uncertainties of the TMDs, which in turn can be fitted to experimental data in a systematic global manner.

Backward Evolution for initial state TMD Parton Shower
A backward evolution method, as now common in Monte Carlo event generators, is applied for the initial state parton shower, evolving from the large scale of the matrix-element process backwards down to the scale of the incoming hadron. However, in contrast to the conventional parton shower, which generates transverse momenta of the initial state partons during the backward evolution, the transverse momenta of the initial partons of the hard scattering process is fixed by the TMD and the parton shower does not change the kinematics. The transverse momenta during the backward cascade follow the behavior of the TMD. The hard scattering process is obtained as described in section 2. The backward evolution of the initial state parton shower follows very closely the description in [60,80,81], which is based on Ref. [53].
The starting value of the evolution scale µ is calculated from the hard scattering process, as described in Section 2. In case of on-shell matrix elements at NLO, the transverse momentum of the hardest parton in the parton shower evolution is limited by the shower-scale, as described in Section 2.1.
Starting at the hard scale µ = µ i , the parton shower algorithm searches for the next scale µ i−1 at which a resolvable branching occurs (see Fig. 3 left). This scale µ i−1 is selected from the Sudakov form factor ∆ bw as given in Eq.(9) (see also [60]). In the parton shower language, the selection of the next branching comes from solving R = ∆ bw (x, k t , µ i , µ i−1 ) for µ i−1 using uniformly distributed random numbers R for given x and µ i . However, to solve the integrals in Eq.(9) numerically for every branching would be too time consuming, instead the vetoalgorithm [53,82] is applied.
The splitting function P ab as well as the argument f (z, q) in the calculation of α s is chosen exactly as used in the evolution of the parton density. In a parton shower one treats "resolvable" branchings, defined via a cut in z < z M in the splitting function to avoid the singular behavior of the terms 1/(1 − z), and branchings with z > z M are regarded as "nonresolvable" and are treated similarly as virtual corrections: they are included in the Sudakov form factor ∆ bw . The splitting variable z i−1 is obtained from the splitting functions following the standard methods (see Eq.(2.37) in [19]).
The calculation of the transverse momentum k t is sketched in Fig. 3 (right). The transverse momentum q t c can be calculated in case of angular ordering (where the scale q of each branching is associated with the angle of the emission) in terms of the angle Θ of the emitted parton with respect to the beam directions q t,c = (1 − z)E b sin Θ, Once the transverse momentum of the emitted parton q c is known, the transverse momentum of the propagating parton can be calculated from with a uniformly distributed azimuthal angle φ assumed for the vector components of k and q c . The generation of the parton momenta is performed in the center-of-mass frame of the collision (in contrast to conventional parton showers, which are generated in different partonic frames). The whole procedure is iterated until one reaches a scale µ i−1 < q 0 with q 0 being a cutoff parameter, which can be chosen to be the starting evolution scale of the TMD. It is of advantage to continue the parton shower evolution to lower scales q 0 ∼ Λ qcd ∼ 0.3 GeV.
The final transverse momentum of the propagating parton k is the sum of all transverse momenta q c (see Fig. 3 right): with k 0 being the intrinsic transverse momentum. The PB TMD parton shower is selected with PartonEvolution=2 (or ICCF=2).

CCFM parton evolution and parton shower
The CCFM parton evolution and corresponding parton shower follows a similar approach as described in the previous section and in detail also in Refs. [60,80,81,83]. The main difference to the PB-TMD shower are the splitting functions with the non-Sudakov form factor ∆ ns and the allowed phase space for emission. The original CCFM splitting functionP g (z, q, k t ) for branching g → gg is given by 2 where the non-Sudakov form factor ∆ ns is defined as with q t = q 2 t being the magnitude of the transverse vector defined in Eq. (11) and k t the magnitude of the transverse vector in Eq. (12).

The TMD parton densities
In the previous versions of CASCADE the TMD densities were part of the program. With the development of TMDLIB [84,85] there is easy access to all available TMDs, including parton densities for photons (as well as Z, W and H densities, if available). These parton densities can be selected via PartonDensity with a value > 100000. For example the TMDs from the parton branching method [19,20] are selected via PartonDensity=102100 (102200) for PB-NLO-HERAI+II-2018-set1 (set2).
Note that the features of the TMD parton shower are only fully available for the PB-TMD sets and the CCFM shower clearly needs CCFM parton densities (like for instance [86]). PB-TMD parton densities are determined in Ref. [87] from fits to HERA DIS F 2 measurements for Q 2 > 3 GeV 2 , giving very good χ 2 values. In Refs. [88,89] the transverse momentum distribution of Drell-Yan pairs at low and high masses, obtained from PB-TMD densities, are compared with experimental measurements in a wide variety of kinematic regions, from low-energy fixed target experiments to high-energy collider experiments. Good agreement is found between predictions and measurements without the need for tuning of nonperturbative parameters, which illustrates the validity of the approach over a broad kinematic range in energy and mass scales.

Final state parton showers
The final state parton shower uses the parton shower routine PYSHOW of PYTHIA. Leptons in the final state, coming for example from Drell-Yan decays, can radiate photons, which are also treated in the final state parton shower. Here the method from PYADSH of PYTHIA is applied, with the scale for the QED shower being fixed at the virtuality of the decaying particle (for example the mass of the Z-boson).

Hadronization
The hadronization (fragmentation of the partons in colorless systems) is done exclusively by PYTHIA. Hadronization (fragmentation) is switched off by Hadronization = 0 (or NFRA = 0 for the older steering cards). All parameters of the hadronization model can be changed via the steering cards.

Uncertainties
Uncertainties of QCD calculations mainly arise from missing higher order corrections, which are estimated by varying the factorization and renormalization scales up and down by typically a factor of 2. The scale variations are performed when calculating the matrix elements and are stored as additional weights in the LHE file, which are then passed directly via CAS-CADE3 to the HEPMC [90] output file for further processing.
The uncertainties coming from the PDFs can also be calculated as additional weight factors during the matrix element calculation. However, when using TMDs, additional uncertainties arise from the transverse momentum distribution of the TMD. The PB-TMDs come with uncertainties from the experimental uncertainties as well as from model uncertainties, as discussed in Ref. [87]. These uncertainties can be treated and applied as additional weight factors with the parameter Uncertainty_TMD=1.

Multi-jet merging
Showered multijet LO matrix element calculations can be merged using the prescription discussed in Ref. [91]. The merging performance is controlled by the three parameters Rclus, Etclus, Etaclmax. Final-state partons with pseudorapidity η <Etaclmax present in the event record after the shower step but before hadronization are passed to the merging machinery if Imerge = 1. Partons are clustered using the kt-jet algorithm with a cone radius Rclus and matched to the PB evolved matrix element partons if the distance between the parton and the jet is R < 1.5×Rclus. The hardness of the reconstructed jets is controlled by its minimum transverse energy Etclus (merging scale).
The number of light flavor partons is defined by the NqmaxMerge parameter. Heavy flavor partons and their corresponding radiation are not passed to the merging algorithm. All jet multiplicities are treated in exclusive mode except for the highest multiplicity MaxJetsMerge which is treated in inclusive mode.

Program description
In CASCADE3 all variables are declared as Double Precision. With CASCADE3 the source of PYTHIA 6.428 is included to avoid difficulties in linking.

Random Numbers
CASCADE3 uses the RANLUX random number generator, with luxory level LUX = 4. The random number seed can be set via the environment variable CASEED, the default value is CASEED=12345.

Event Output
When HEPMC is included, generated events are written out in HEPMC [90] format for further processing. The environment variable HEPMCOUT is used to specify the file name, by default this variable is set to HEPMCOUT=/dev/null.
The HEPMC events can be further processed, for example with Rivet [92].

Input parameters
The input parameters are steered via steering files. The new format of steering is discussed in Section 9.3.1 and should be used when reading LHE files, while the other format, which is appropriate for the internal off-shell processes, is discussed in Section 9.3.2.

Method of solution:
Since measurements involve complex cuts and multi-particle final states, the ideal tool for any theoretical description of the data is a Monte Carlo event generator which generates initial state parton showers according to Transverse Momentum Dependent (TMD) parton densities, in a backward evolution, which follows the evolution equation as used for the determination of the TMD.
Restrictions on the complexity of the problem: Any LHE file (with on-shell or off-shell) initial state partons can be processed.
Other Program used: PYTHIA (version > 6.4) for final state parton shower and hadronization, BASES/SPRING 5.1 for integration (both supplied with the program package), TMDLIB as a library for TMD parton densities.
Download of the program: http://www.desy.de/˜jung/cascade Unusual features of the program: None