The CCFM Monte Carlo generator CASCADE 2.2.0

CASCADE is a full hadron level Monte Carlo event generator for ep, \gamma p and p\bar{p} and pp processes, which uses the CCFM evolution equation for the initial state cascade in a backward evolution approach supplemented with off - shell matrix elements for the hard scattering. A detailed program description is given, with emphasis on parameters the user wants to change and variables which completely specify the generated events.


The CCFM evolution equation
The formulation of the CCFM [1][2][3][4] parton evolution for the implementation into a full hadron level Monte Carlo program is described in detail in [5,6]. Here only the main results are summarized and discussed. The pattern of QCD initial state radiation in a small-x   Fig. 1 together with labels for the kinematics. According to the CCFM evolution equation, the emission of partons during the initial cascade is only allowed in an angular ordered region of phase space. In terms of Sudakov variables Υ and Ξ the quark pair momentum is written as: where p (1) and p (2) are the four-vectors of incoming particles (electron-proton, proton-antiproton or proton-proton), respectively and Q t is the transverse momentum of the quark pair in the center of mass frame of p (1) and p (2) (cms). The variable Ξ is related to the rapidity Y in the center of mass (CMS) frame via Using E = E q + Eq and p z = p q z + pq z gives E + p z = Υ √ s, E − p z = ΥΞ √ s with E = √ s/2 and s = (p (1) + p (2) ) 2 being the squared center of mass energy. Therefore Ξ can be used to define the maximum allowed angle in the evolution. The momenta p i of the gluons emitted during the initial state cascade are given by (here treated massless): with υ i = (1 − z i )x i−1 and x i = z i x i−1 . The variables x i and υ i are the momentum fractions of the exchanged and emitted gluons, while z i is the momentum fraction in the branching (i − 1) → i and p ti is the transverse momentum of the emitted gluon i. Again the rapidities y i are given by y i = −0.5 log ξ i in the CMS frame. The angular ordered region is then specified by ( Fig. 1a and the lower part of the cascade in Fig. 1b, for the upper part the variables have to be changed accordingly): which becomes: where the rescaled transverse momentum q i of the emitted gluon is defined by: The scaleq (related to the maximum angle) can be written as: withŝ = (p q + pq) 2 and the relation ofq to a particular choice of the factorisation scale µ f in the collinear approach becomes obvious. The CCFM evolution equation can be written in a differential form [4], which is best suited for a backward evolution approach adopted in the Monte Carlo generator CASCADE [5,6]:q where A(x, k t ,q) is the unintegrated gluon density, depending on x, k t and the evolution variableq. The splitting variable is z = x/x ′ and k t ′ = (1 − z)/z q + k t , where the vector q is at an azimuthal angle φ. The Sudakov form factor ∆ s is given by: withᾱ s = C A αs π = 3αs π . For inclusive quantities at leading-logarithmic order the Sudakov form factor cancels against the 1/(1 − z) collinear singularity of the splitting function.
The original splitting functionP g (z i , q i , k ti ) for branching i is given by (neglecting finite terms as they are not obtained in CCFM at the leading infrared accuracy (cf p. 72 in [3]): where the non-Sudakov form factor ∆ ns is defined as: The implementation of the full splitting function including non singular terms can lead to inconsistencies. Replacing naively only 1 1−z → 1 1−z − 2 + z(1 − z) in the CCFM splitting function can lead to negative branching probabilities.
In [7] it was suggested to use: where B is a parameter to be chosen arbitrarily between 0 and 1, we take B = 0.5. As a consequence of the replacement, the Sudakov form factor will change, but also the non-Sudakov form factor needs to be replaced by:

Backward evolution: CCFM and CASCADE
The idea of a backward evolution [8,9] is to first generate the hard scattering process with the initial parton momenta distributed according to the parton distribution functions. This involves in general only a fixed number of degrees of freedom, and the hard scattering process can be generated quite efficiently. The initial state cascade is generated by going backwards from the hard scattering process towards the beam particles. According to the CCFM equation the probability of finding a gluon in the proton depends on three variables, the momentum fraction x, the transverse momentum squared k 2 t of the exchanged gluons and the scaleq = x n √ sΞ, which is related to the maximum angle Ξ allowed for any emission. To solve eq.(8) the unintegrated parton distribution A(x, k t ,q) has to be determined beforehand.
Given A(x, k t ,q), the generation of a full hadronic event is separated into three steps, as implemented in the hadron level Monte Carlo program CASCADE: • The hard scattering process is generated, with k 1 (k 2 ) being the momenta of the incoming partons to the subprocess k 1 + k 2 → X with X being the final state. The definition ofσ follows [10]. The available processes are shown in Tab 2. The momenta of the incoming partons are given in Sudakov representation: where the last expression comes from the high energy approximation (x 1, 2 ≪ 1), which then gives −k 2 ≃ k 2 t .
• The initial state cascade is generated according to CCFM in a backward evolution approach.
• The hadronisation is performed using the Lund string fragmentation implemented in PYTHIA /JETSET [11].
The parton virtuality enters the hard scattering process and also influences the kinematics of the produced particles (Z 0 , W, Higgs and quarks) and therefore the maximum angle allowed for any further emission in the initial state cascade. This virtuality is only known after the whole cascade has been generated, since it depends on the history of the parton evolution (asx in eq.(15) may not be neglected for exact kinematics). In the evolution equations itself it does not enter, since there only the longitudinal energy fractions z and the transverse momenta are involved. This problem can only approximately be overcome by using for the virtuality which is correct in the case of no further parton emission in the initial state.
The Monte Carlo program CASCADE can be used to generate unweighted full hadron level events, including initial state parton evolution according to the CCFM equation and the off-shell matrix elements for the hard scattering process. It is applicable for pp, pp, photoproduction as well as for deep inelastic scattering. A discussion of the phenomenological applications of CASCADE can be found in [12].
The typical time needed to generate one event is similar to the time needed by standard Monte Carlo event generators like PYTHIA [11]. The CCFM unintegrated parton density xA(x, k t ,q) can be obtained from a forward evolution procedure as implemented in SMALLX [16,17] by a fit to the measured structure function F 2 as described i.e. in [5,6]. From the initial parton distribution, which includes a Gaussian intrinsic k t distribution, a set of values x and k t is obtained by evolving up to a given scale logq using a forward evolution procedure. Technically the parton density is stored on a grid in log x, log k t and logq and a linear interpolation is used to obtain the parton density for values in between the grid points. The data file (i.e. ccfm-xxxx.dat) containing the grid points is read in at the beginning of the program.

The unintegrated parton density
Several sets (J2003 set 1 -3 [15] ( IGLU=1001-1003), set A [13] IGLU=1010-1013 and set B [13] IGLU=1020-1023) of unintegrated gluon densities are available with the input parameters fitted to describe the structure function F 2 (x, Q 2 ) in the range x < 5 · 10 −3 and Q 2 > 4.5 GeV 2 as measured at H1 [18,19] and ZEUS [20,21]. Set JS [6] (IGLU=1) is fitted only to F 2 (x, Q 2 ) of Ref. [18]. The collinear cutoff k cut t = Q 0 which regulates the region of z → 1, is applied both to the real emissions as well as inside the Sudakov form factor. In JS, J2003 set 1 -3 and set A we have k cut t = Q 0 = 1.3 GeV. Similarly, fits can be obtained using different values for the soft cut k cut t = 0.25 GeV, which are available in set B. The set C (IGLU=1101) [14] uses the full splitting function, as described in eq.(12), with a value for Λ (4) QCD = 0.13 GeV. This set was obtained by fitting simultaneously the inclusive F 2 (x, Q 2 ) and jet measurements in DIS resulting in a changed intrinsic k t distribution. A CCFM parametrisation of the valence quark distribution is available (using as starting distribution CTEQ 5 [22] and evolved with a splitting function P qq [3] including angular ordering of the emitted gluon). The available CCFM uPDF sets with the parameters of the starting distributions are listed in Tab. 2.1.
Initial state parton showers can be only generated for the CCFM unintegrated gluon density (with IGLU=1 and IGLU=1001-1023,1101). For all other sets only the cross section can be calculated without explicit inclusion of initial state parton showers, since the angular variable, essential for angular ordering in the initial state cascade, is not available in the uPDFs. However, the transverse momenta of the incoming partons are properly treated. Only the KMR set (IGLU=6) provides a prescription for the emission of at most one additional gluon.

Hard processes in CASCADE
Different sets of hard processes applicable for lepto (photo) -and hadroproduction have been calculated and are implemented in CASCADE. The available processes are listed in Tab. 2.

Lepto(photo)production
CASCADE can be used to simulate leptoproduction events over the whole Q 2 range. By fixing the light quark masses to m q = 0.25 GeV and α s for small µ, the hard scattering matrix element remains finite over the full phase space. The total cross section is simulated by selecting IPRO=10 and NFLAV=4(5). With IPRO=10 light quarks (u, d, s) are selected and with NFLAV>3 the program automatically includes heavy flavour production via the process IPRO=11 and IHFLA=4 up to IHFLA=NFLAV. The flag IRE1 indicates, whether beam 1 has a internal structure: IRE1=1 is used to generate resolved photon events.
Heavy flavour production can be generated separately via IPRO=11. The value of IHFLA determines the heavy flavour to be generated.
The matrix element for γg * → J/ψ(Υ)g calculated in [31][32][33] is available for quasi-real γ's via the process IPRO=2. The flavour of the Onium is selected via IHFLA, i.e. IHFLA=4 for Lepto(photo)production process J/ψ and IHFLA=5 for Υ. The matrixelement including J/ψ(Υ) polarisation and subsequent leptonic decay can be selected with IPSIPOL=1. CASCADE can be used to simulate real photoproduction events by using KBE1=22. The same options as for leptoproduction are available. Resolved photon events can be generated with IRE1=1.

Hadroproduction
The hadroproduction processes available are listed in Tab 2. The flavour code for beam 1 (2) can be chosen as KBE1=2212 for proton or KBE1=-2212 for anti -proton, for beam 2 KBE2 is changed accordingly.
CASCADE can be used to simulate heavy quark production in pp or pp collisions (g * g * → QQ IPRO=11 for heavy flavour production, and IHFLA=4(5) for charm (bottom) quarks), but also for light quarks with IPRO=10. The matrix element for g * g * → J/ψg calculated in [31][32][33] is available via the process IPRO=2. The matrixelement including J/ψ polarisation and subsequent leptonic decay can be selected with IPSIPOL=1. The process g * g * → χ is available with IPRO=3 including all three χ states with appropriate spin and angular momentum. The flavour of the Onium is selected via IHFLA, i.e. IHFLA=4 for J/ψ(χ c ) and IHFLA=5 for Υ(χ b ).
The process g * g * → h 0 with the matrix element calculated in [34] is available via IPRO=102, the Higgs mass can be selected via PMAS (25).
The process g * g * → Zqq, calculated in [35,36], is available via IPRO=503 for light quarks and IPRO=504 for the heavy quarks. The flavor of the heavy quark is selected via IHFLA=4 (5,6) for charm(bottom, top). The process g * g * → W q i q j , calculated in [35,36], is available via IPRO=513 for light quarks and IPRO=514 for heavy quarks, where the flavour of the heaviest quark is defined by IHFLA=3, (4,5). Which of the quarks are produced depends on the charge of the W which is randomly selected.
The processes g * g → gg and g * q → gq [38] are included via process IPRO=10. The individual processes can be selected for g * g * → qq via IRPA=1, g * g → gg via IRPB=1 and g * q → gq via IRPC=1. Note that here one of the partons is treated on-shell. For the quarks the unintegrated quark distribution (for valence quarks) is used.

α s and the choice of scales
The strong coupling α s is calculated via the PYTHIA [11] subroutine PYALPS. Maximal and minimal number of flavours used in α s are set by MSTU(113) and MSTU(114), Λ QCD = PARU(112) with respect to the number of flavours given in MSTU(112) and stored in the PYTHIA common block COMMON/PYDAT1/. In the initial state cascade according to CCFM, the transverse momenta of the t-channel gluons are allowed to perform a random walk for small z values and k t can become very small. In the 1/z part of the splitting function we use µ = k t as the scale in α s (µ) and in the 1/(1 − z) part µ = p t is used. In addition we require µ > Q 0 , resulting in α s (µ > Q 0 ) < 0.6.
The scale µ, which is used in α s in the hard scattering matrix element, can be changed with the parameter IQ2, the default choice is µ 2 = p 2 t . The renormalisation scale dependence of the final cross section can be estimated by changing the scale used in α s in the off-shell matrix element. Since here we are using the LO α s matrix elements, any scale variation will change the cross section. In order to obtain a reasonable result, the uPDF was fitted to describe F 2 by varying the scale µ r . The set A0-,B0correspond to a scale µ r = 0.5p t whereas set A0+,B0+ correspond to a scale µ r = 2p t .
In order to investigate the uncertainties coming from the specific choice of the evolution scale, another definition is applied, relating the factorisation scale only to the quark (or antiquark): µ f = pt 1−z with p t being the transverse momentum of the quark (anti-quark) and z =k t yxgs . The set A1,B1 correspond to a scale µ f = pt 1−z . In the PDF set C Λ QCD was fixed to Λ (4) QCD = 0.13 GeV.

Quark masses
The quark mass for light quarks (u, d, s) is fixed to m q = 0.25 GeV. This, together with the treatment of α s at small scales µ, gives also a reasonable total cross-section for photoproduction at HERA energies. The masses for heavy quarks are given by the JETSET / PYTHIA defaults (m c = 1.5 GeV, m b = 4.8 GeV) and can be changed according to the PYTHIA prescription.

Initial and final state parton showers
Initial state parton showers are generated in a backward evolution approach described in detail in [5,6]. The initial state parton shower consists only of gluon branchings and is generated in an angular ordered region in the laboratory frame. The gluons emitted in the initial state can undergo further timelike branchings. The maximum timelike mass m max is calculated using the angular constraint. With this mass, the parton which can further undergo a timelike branching is boosted to its rest frame with m max but keeping the original energy. The timelike branching is performed with the PYTHIA routine PYSHOW. After successful timelike branching, the proper mass is associated to the parton and the kinematics are calculated appropriately. Gluon radiation from the valence quarks is also included.
All parameters (like the scale µ in α s , the collinear cut-off k cut t = Q 0 ) for the initial state cascade are fixed from the determination of the unintegrated gluon density. The transverse momenta of the partons which enter the hard scattering matrix element are already generated in the beginning and are not changed during the whole initial and final state parton showering.

Remnant treatment
In CASCADE version 1 the proton remnant was built in subroutine CAREMN, which is a slightly modified version of the PYTHIA/LEPTO subroutine PYREMN. No intrinsic transverse momentum, in addition to the transverse momentum from the initial state cascade, was included.
From version 2.0 on the proton remnant can be generated directly via PYTHIA by selecting ILHA=10 (which is now the default). The structure of the event record is then identical to the one obtained from a standard PYTHIA run.

Hadronization
In CASCADE version 1 the hadronization was done exclusively by PYTHIA. From version 2 onwards, events can be written into a file (via switch ILHA=1) according the LHA Accord [39], which can be read by any hadronization program (like PYTHIA or HERWIG ), generating the remnants and performing the hadronization. With ILHA=10 the hadronization is performed within PYTHIA. The old CASCADE format is obtained with ILHA=0. Please note, that top decays can only be simulated properly within the PYTHIA fragmentation and are therefore available only with ILHA=10.

Description of the program components
In CASCADE all variables are declared as Double Precision. The Lund string model is used for hadronization as implemented in PYTHIA [11]. The final state QCD radiation is performed via PYSHOW from PYTHIA . The treatment of the proton remnant follows very closely the ones in LEPTO [40] for the leptoproduction case and the one in PYTHIA for the proton -proton case. However slight modifications were needed to adapt to the cascade treatment here.
The unintegrated gluon density is stored on data files (ccfm.dat,kms.dat,kmr.dat), and is read in at the beginning of the program.
The program has to be compiled and linked together with PYTHIA 6, to ensure that the double precision code of JETSET is loaded.

Random number generator
Since the variables are declared as double precision, also a double precision random number generator has to be used to avoid any bias. The function DCASRN gives a single random number, the function DCASRNV returns an array of length LEN of random numbers. The default random number generator is RANLUX (called in DCASRN and DCASRNV) The source code of RANLUX (extracted from CERNLIB) is included in the distribution. The user can change this to any preferred Double Precision random number generator.

Integration and event generation
The integration of the total cross section and the generation of unweighted events is performed with the help of BASES/SPRING [41], which is included in the distribution package.

Subroutines and functions
The source code of CASCADE and this manual can be found under: http://www.desy.de/˜jung/cascade/ CAMAIN main program. CASINI to initialises the program. CASCADE to perform integration of the cross section. This routine has to be called before event generation can start. CAEND to print the cross section and the number of events. CAUNIGLU(KF,X,KT,P,XPQ) to extract the unintegrated gluon density xA(x, k t ,q) for a proton with KF=2212, as a function of x =X, k 2 t =KT andq =P. The gluon density is returned in XPQ(0), where XPQ is an array with XPQ(-6:6). EVENT to perform the event generation.

FXN1
to call routines for selected processes: XSEC1. CUTG (IPRO) to cut on p t for 2 → 2 process in integration and event generation. MEOFFSH matrix element for γ * g * → qq and g * g * → qq including quark masses. q can be light or heavy quarks. MEHIGGS matrix element for γ * g * → h 0 .

DOT(A,B)
four-vector dot product of A and B.

DOT1(I,J)
four-vector dot product of vectors I and J in PYJETS common. PHASE to generate momenta of final partons in a 2 → 2 subprocess according to phase space. P SEMIH to generate kinematics and the event record for ep, γp and pp processes. CAREMN(IPU1,IPU2) to generate the beam remnants. Copied from LEPTO 6.1 [40] and updated for the use in CASCADE. CASPLI(KF,KPA,KFSP,KFCH) to give the spectator KFSP and KFCH partons when a parton KPA is removed from particle KF. Copied from LEPTO 6.1 [40] and updated for the use CASCADE.

CAPS
to generate color flow for all processes and prepare for initial and final state parton showers. CASCPS(IPU1,IPU2) to generate initial state radiation. COLORFLOW to generate color configuration for g * g → gg and g * q → qg processes.