Drell-Yan production with the CCFM-K evolution

We discuss the Drell-Yan dilepton production using the transverse momentum dependent parton distributions evolved with the Catani-Ciafaloni-Fiorani-Marchesini-Kwieci\'nski (CCFM-K) equations in the single loop approximation. Such equations are obtained assuming angular ordering of emitted partons (coherence) for $x\sim 1$ and transverse momentum ordering for $x \ll 1$. This evolution scheme also contains the Collins-Soper-Sterman (CSS) soft gluon resummation. We make a comparison with a broad class of data on transverse momentum spectra of low mass Drell-Yan dileptons.


Introduction
The Drell-Yan (DY) dilepton production [1] is one of the most intensively studied processes in particle physics. The existence of a hard electroweak probe (photon or Z boson), which doesn't interact strongly and decays into a pair of leptons, provides a clear experimental signature of the partonic interactions in the colliding hadrons and significantly simplifies their theoretical description. For these reasons, the DY process is very efficient tool for the investigation of hadronic structure [2], in particular the distributions of partons' transverse momentum. The key concept of the theoretical analysis of DY scattering, called factorization, is a separation between long-range and short-range degrees of freedom. The basic, collinear factorization theorem assumes no transverse momenta of partons in hadrons [3]. Its application to the Drell-Yan production is very well established and commonly used. The present state of art calculation for the DY process includes next-to-next-to-leading order (NNLO) QCD corrections [4][5][6].
There are several kinematical regimes where the fixed-order collinear QCD doesn't provide good description of data. In particular, when transverse momentum of lepton pair is much smaller than it's invariant mass, q T M , the large logarithms log n (M/q T ) occur in all orders of perturbative expansion. These corrections are effectively resumed within the Collins, Soper and Sterman (CSS) approach [7]. As a result, the collinear factorization need to be replaced by a new factorization based on the transverse momentum distributions (TMDs) [8]. It is interesting to ask what is the transition from the small transverse momentum region to the region of q T ∼ M , where fixed order perturbative QCD should apply. As it was shown in [9], when one considers moderate values of M ∼ 5 − 19 GeV, the fixed-order predictions underestimate data in the region of q T ∼ M .
In this paper, we examine in detail the low mass DY productions using the k T -factorization approach combined with the transverse momentum parton distributions defined in the Catani-Ciafaloni Fiorani-Marchesini branching scheme (CCFM) [41][42][43][44]. The original motivation for the CCFM branching was to extend angular ordering (coherence) of soft parton emission from the region of x ∼ 1 into the region of small x 1. In this way, a unified evolution equation for transverse momentum dependent gluon distribution was found with angular ordering in both regions of x in the approximation called all loop. In the fully inclusive case, this equation interpolates between the DGLAP equation at moderate x and the BFKL equation obtained in the small x limit. The Monte Carlo implementation of the all loop CCFM branching scheme was done in [45][46][47]. Further extensions of this implementation and subsequent analyses were presented in [48][49][50][51].
In the region of large or moderate values of x, when the small x coherence can be neglected, the CCFM scheme gives the DGLAP evolution equation with coherence at large x only. This approximation, called single loop, was studied in [44,45] for the gluon distribution. The extension of this scheme in terms of evolution equations for both quark and gluon distributions was proposed by Kwieciński in [52]. This is why we call them the CCFM-K equations. The parton distributions which are obtained by solving these equations depend on transverse momenta and their properties were analysed in [53][54][55]. The first analysis of the weak gauge boson production with the CCFM-K equations was done in [56] while the low mass DY production with photon was addressed in [57].
The main goal of the presented analysis is a comprehensive analysis of all available data on transverse momentum spectra in low mass DY production with the CCFM-K evolved parton distributions and the leading order cross sections computed in k T factorization. We assume the most economical form of the initial conditions for the evolution with only one adjustable parameter. In this way, we concentrate on the most important effects of the CCFM-K evolution which are responsible for good description of data for q T ∼ M . The small q T description is acceptable in most cases, especially for higher masses, although precise comparison should involve more adjustable parameters in the initial conditions like in the CSS approach. This is left for future studies.
The paper is organized as follows. In section 2 we give an overview of the CCFM framework and its version proposed by Kwieciński in which quark and gluon transverse momentum dependent distributions and the CCFM-K evolution equations are introduced. We also discuss the relation of the CCFM-K approach to the CSS formalism (with the full derivation presented in Appendix A). In section 3 we describe the application of the discussed formalisms to the leading order DY cross section with both on-shell and off-shell matrix elements in the CCFM-K case. In section 4 we show numerical results and compare them with the low mass DY data. We summarize in section 5.

CCFM approach 2.1 Branching kinematics
Below we describe kinematics of the CCFM parton branching schemes in two approximationssingle and all loop. We work in the Sudakov base with two light cone vectors S(1, 0, 0, 1) in which the momenta in the branching process shown in Fig. 1 are given by Notice that the momentum fraction are proportional to the plus/minus components, e.g.
The emitted parton momentum can be found from the momentum conservation where and assuming that p 2 i = 0, one can compute the minus component to find The rapidity of the emitted parton is given by where z i = x i /x i−1 < 1. In the massless case y i = − ln tan(θ i /2) where θ i is the emission angle with respect to the z axis defined by the collinear momenta P 1 and P 2 , therefore The CCFM branching scheme [44] is defined with the help of the rescaled transverse momentum of the emitted parton, Thus, the transverse momentum conservation at the vertex i reads From (7) we obtain for the modulus The single loop approximation is defined by the condition Therefore, for given x and √ S, the hard scale determines the maximal emission angle The second theta function in (17) imposes the condition which assures that α s (q) 1 and the CCFM evolution scheme is perturbative. Non-perturbative effects are encoded in the initial condition, f 0 g (x, k T , Q 0 ), imposed at a scale Q 0 . From the third theta function in (17), we find the condition for the real gluon emission, 0 < z < (1 − Q 0 /q), which allows to avoid singularity of P gg (z) at z = 1.
Eq. (17) can be used for the Monte Carlo generation of a parton cascade with angularly ordered emissions which, when resumed, give the gluon distribution f g . However, the mixing between the transverse and longitudinal variables in the theta function θ(Q−zq) prevents writing eq. (17) in the form of the evolution equation. This can be done in the single loop approximation in which the branching scheme leads to evolution equations of the DGLAP type for both quark and gluon distributions with the evolution scale defined by the rescaled momentum Q. We are describing them below

CCFM-Kwieciński evolution equations
In this approximation, the all loop branching conditions in eq. (17) are replaced by [52][53][54][55] Thus, for z → 0, the angular ordering is replaced by the stronger transverse momentum ordering while for z → 1 the angular ordering is still valid. In addition, quark splittings, q → qg,q →qg and g → qq, are taken into account which allow to introduce quark distributions f i=1,...,2N f in addition to the gluon distribution f g . The CCFM-K equations for quark and gluon distributions have the following form [55] where the argument of the parton distributions on the r.h.s. equals The one loop splitting functions are given by where C F = 4/3, C A = 3, T R = 1/2 and N f is the number of active flavours. Notice that after integrating both sides of eqs. (24) over k T , the ordinary DGLAP equations for the collinear quark and gluon distributions, are found. Eqs. (24) can be written with the help of the Fourier transformatioñ which for b = 0 gives the PDFs (27). Since the parton distributions depend on k T = | k T |, we perform the azimuthal angle integration with help of the relation and obtain the parton distributions which depend on b = | b|, Thus, taking the Fourier transform of both sides of eqs. (24), we find the evolution equations which are diagonal in b: These are the CCFM-K equations which we use in our forthcoming analysis. As expected, for b = 0 we obtain the DGLAP evolution equations for the collinear PDFs (27) .

Initial conditions and b-dependence
In order to solve eqs. (31), we need initial conditions specified as functions of x and b at some perturbative scale Q 0 Λ QCD . They have to fulfill the conditions saying that for b = 0 the collinear PDFs are recovered. Thus the simplest possible choice is given in the factorized form where q i and g are the LO collinear quark and gluon distributions and the non-perturbative form factor obeys the condition F (0) = 1. In the forthcoming analysis we will use the gaussian form factor where b 0 is a free parameter. In principle, different form factors can be used for quarks and gluons. However, with the common form factor, we can disentangle the non-perturbative b-dependence given by F (b) and the b-dependence due to the CCFM-K evolution equations (31). Since they are homogeneous, we find their solution f CCFMK i,g (x, b, Q) starting from the initial conditions (32) with F (b) ≡ 1 and then multiplying it by F (b) from (33), In this way, the perturbative and non-perturbative b-dependences are clearly separated. This effect is shown in Fig. 2 for the singlet quark distribution, f S = i f i , plotted as a function of b for fixed x = 10 −2 . The dashed curves are the initial conditions at Q 0 = 1 GeV with the MSTW08 LO PDFs [58] while the solid curves show the distributions evolved to Q = 100 GeV.
On the left plot, the b-dependence of the evolved curve is purely perturbative while on the right plot the curves were multiplied by the form factor F (b). We see that form factor's impact is the strongest for large values of b while for small values the b-dependence of the full solution remains perturbative. After the Fourier transformation to the k T -space, we find broadening of the parton distributions due to evolution, studied in detail in [54]. We also realize that the shape of the non-perturbative form factor, i.e. the parameter b 0 , influences mostly the small k ⊥ behaviour of the parton distributions.

Relation to CSS resummation
The CCFM-K equations contain the Collins, Soper and Sterman (CSS) resummation of soft parton emissions in the limit z → 1. In Appendix A, we present the proof that for the values of  (36) is shown as the dashed lines. The b-range corresponds to (35) the parameter b such that the solution to the CCFM-K equations (31) is given by the CSS formulas [7] f where c ∼ 1 and the parameters A q,g are defined in eq. (86) in Appendix A. Solution (36) is found by picking large logarithms, ln(Qb) and ln(1/Q 0 b), which are present due to the condition for the values of b. It should be noted that eqs. (36) do not contain NLL terms proportional to α 2 s under the integrals (see section 3.3) since the splitting functions in the CCFM-K equations are given in the leading order proportional to α s . In Fig. 3 we show the full numerical solution to the CCFM-K equations for singlet quark and gluon distributions (solid lines) together with the approximate solution (36) (dashed lines) in the range of b which is given by (35).

Drell-Yan cross section with k T -dependent PDFs
The Drell-Yan cross section differential in photon's momentum is given by where (y γ , M, q T ) are photon's rapidity, virtuality and transverse momentum while W µ µ is the trace of the hadronic tensor W µν . With the lowest order matrix element for the process qq → γ * , the hadronic tensor is given by the transverse momentum factorization formula where k 1T , k 2T are quark/antiquark transverse momenta, x 1,2 are their longitudinal momenta and f i , fī are transverse momentum dependent quark/antiquark distributions taken at the scale Q = M .

On-shell matrix element
The trace in (38) is the squared matrix element of the process q(k 1 )q(k 2 ) → γ * in the lowest order. For the on-shell matrix element, we use the quark/antiquark momenta in the collinear approximation what assures gauge invariance of the matrix elements. In such a case and the DY cross section (37) is given by where the quark/antiquark transverse momentum fractions are given by It is easy to check that after integrating (41) over q T , we find the leading order form of the Drell-Yan cross section with collinear PDFs given by eq. (27) dσ DY where Inserting the delta function we find the DY cross section with the b-dependent parton distributions (28) For the parton distributions which depend on b = | b|, the angular integration with the help of relation (29) gives We will use this expression for the comparison with the Drell-Yan data using the parton distributions which are solutions of the CCFM-K equations.

Off-shell matrix elements
In approach with the off-shell matrix, (40) is given by where Γ µ is the Fadin-Sherman photon-quark vertex [59,60] The quark/antiquark momenta k 1,2 in this vertex take into account transverse momenta where k T i = (0, k T i , 0) for i = 1, 2 and the momentum fractions x i are determined from the momentum conservation at the vertex, which gives It is easy to check that the Fadin-Sherman vertex obeys the gauge invariance relation Computing the trace in eq. (48), we obtain where we used momentum conservation (51) to write the last equality. With this result, the DY cross section (37) is given by Inserting the delta function representation (45) and performing the Fourier transformation (28), we obtain the following cross section where ∆ b is the radial part of the two-dimensional Laplacian Note that the first line of (56) equals to the on-shell formula (47) if we use invariant mass M instead of M T in the momentum fractions (52). The term in the second line of (56) occurs due to the off-shellness of the matrix element.
As an illustration of the impact of the off-shell matrix element, in Fig. 4 we show predictions for the DY cross section dσ/dq 2 T against the Fermilab R209 experiment data [36]. The solid line is the prediction obtained from the on-shell formula (47) with the CCFM-K evolved parton distributions. We used initial conditions (32) with the MSTW08 LO PDFs [58] at Q 0 = 1 GeV and the form factor (33) with b 0 = 2.7 GeV −1 . The dash-dotted line is the prediction for the off-shell formula (56) with the same parton distributions. The difference between the two curves for q T ∼ M results from the leading dependence on photon mass M of the corresponding cross sections. For the same values of x 1 and x 2 , the off-shell cross section (55) is M 2 /M 2 T times smaller than in the on-shell case (41). Different momentum fractions which we use, see (42) and (52), make this difference even bigger. In the kinematics of the R209 experiment, the terms in the second line in the off-shell formula (56) are numerically negligible. The dotted line in Figure  4 shows the prediction for the CSS formula (58) in the NLL approximation, presented in detail in the next section.

CSS approach
We also present the results for the DY cross section computed in the CSS approach [7] in which the soft gluon resummation is perfomed for the photon transverse momentum q T M . In this case, the DY cross section is given by where Q = M and x 1,2 = M e ±yγ / √ S. The parton luminosities W iī are given by where f i/ī are effective quark/antiquark distributions with q i/ī and g are being the MS collinear PDFs and The scale µ in (60) is given by µ = c/b * with c = 2e −γ E ≈ 1.12 and Thus, b * interpolates between 0 and b max for b → ∞ such that µ ∈ [c/b max , ∞). Choosing b max = c/Q 0 , where Q 0 is an initial scale for the DGLAP evolution, we ensure that the PDFs are always defined during the integration over b in (58). We use the MSTW08 NLO PDFs [58] and choose Q 0 = 1 GeV. The exponent S in (59) is given by where the coefficients A q and B q are defined by the general perturbative expansion Introducing B = ln(Q 2 b 2 /c 2 ) and L = L(Q) = ln(Q 2 /Λ 2 ), we define the leading logartithmic (LL) approximation by the terms proportional to B(B/L) n while in the next-to-leading logarithmic (NLL) approximation we add terms proportional to (B/L) n . Thus, the NLL approximation is given by the coefficients [61-63] and K = C A ( 67 18 − 1 6 π 2 ) − 10 9 T R N f . By the comparison of the exponent (63) with that in (36), we see that the CCFM-K equations only partially resum the next-to-leading logarithms since they do not lead to the term proportional to A (2) q in (63), which is formally of the NLL accuracy. Such a term can be obtained from the CCFM-K equations with the two-loop splitting functions. Using the two-loop running coupling constant where β 0 = (11C A − 4T R N f )/12π and β 1 = (153 − 19N f )/24π 2 , and performing the integration in (63), one obtains the final NLL form of S [63], The factor W N P = exp{−S N P } in (59) describes the non-perturbative contribution [64][65][66]. In the comparison with the data, we use the form from the BNLY fit [67] to the DY data: where a 1 = 0.21 GeV 2 , a 2 = 0.68 GeV 2 , a 3 = −0.1 GeV 2 and Q 1 = 3.2 GeV. In this paper we do not aim at a comprehensive description of the data using the CSS approach. It only serves for the comparison with the CCFM-K results. For this reason, we do not consider the so-called "Y term" which was proposed in [7] to obtain matching to fixed-order results. For a full analysis within the CSS approach, see for example [9].

Comparison to data
For the comparison with the low mass DY data, we use the CCFM-K approach with both onshell and off-shell matrix elements (see sections 3.1 and 3.2). In this appraoch, we only have one free parameter, b 0 in the non-perturbative form factor (33), which we set to b 0 = 2.7 GeV −1 .
For the initial PDFs we use the MSTW08 LO PDF set [58]. For the comparison, we also use the CSS formulae at the NLL accuracy discussed in the previous section.

DY from fixed target experiments
We start with the data from the fixed target experiments E288 [32], E605 [33], E866 [35] and E722 [34]. The cross section Ed 3 σ/d 3 q measured in these experiments is related to (47) and (58) as follows where ∆M is the bin size of the DY pair mass distribution. In addition, the E605, E866 and E722 experiments also measured the cross section in bins of the Feynman variable where the r.h.s. gives the relation between x F and the DY photon rapidity y γ . The energies of the proton projectile were equal to: 200, 300 and 400 GeV (at E288), and 800 GeV (at E605, E866 and E772). These translate into the center of mass energies √ S = 19.4, 23.8, 27.4 and 38.8 GeV, respectively. The experiments differ by the targets used: E288 used Cu or Pt, E605 used Cu, E866 used H or D and E772 used 2 H. All the cross sections were normalized by the number of nucleons in the target nucleus. In what follows, we neglect nuclear effects of the targets and compare unmodified CCFM-K and CCS approaches with such data.
In Fig. 5 we show the data from the E288 experiment [32] for three values of energies and rapidities. At each panel, the transverse momentum dependence of the DY cross section is shown with fixed mass M and rapidity y (or x F ) 1 . The mass range equals M = 4.5−12.5 GeV. In Fig. 6 we show the data from the E605 [33] and E866 experiments [35] for √ S = 38.8 GeV, x F = 0.1 and the mass range M = 4.7 − 15.5 GeV. In Fig. 7 we show the data from the E866 experiment for √ S = 38.8 GeV, three values of x F and the mass range M = 4.7 − 14.85 GeV. Finally, in Fig. 8 we show the data from the E772 experiment [34] for √ S = 38.8 GeV, 0.1 < x F < 0.3 and the mass range M = 5.5 − 14.5 GeV.
The data in Figs. 5-8 are compared to theoretical curves: CCFM-K on-shell (blue solid curves), CCFM-K off-shell (blue dashed-dotted curves) and CSS (red dashed curves). Comparing CCFM-K to CSS we see that at small q T the CSS resummation predicts higher cross-section than CCFM-K and better agrees with the E288 and E605 data. This comes from the fact that the parameters of the non-perturbative form factor (68) were fitted in [67] to the data while in the CCFM-K approach we fixed just one free parameter b 0 in the form factor (33) The motivation for that was to show the potential of the CCFM-K approach to describe the large q T data without going into details of fitting the parameter b 0 , as in the CSS approach. Thus, at larger q T (2-3 GeV, depending on M and x F ), the CSS curves drop rapidly as we do not match them to the fixed order calculation by adding the "Y term". On the other hand, the CCFM-K curves describe the data reasonably well. Note also that for M ∼ 9 GeV, the data from E288 are significantly above the theoretical predictions which is related to the production of Υ meson, not considered in our calculations.
The E866 and E772 data seems to be systematically above theoretical predictions at small q T , except for a few values of M and x F . As before, CCFM-K provides good description of the data at large q T . In general, one sees better description of data for higher DY masses.
Comparing the CCFM-K on-shell and off-shell approaches one sees that the former approach agrees better with data as providing slower drop with q T . The difference is larger for q T ∼ M , as one should expect, see discussion at the end of section 3.2.

DY in proton-proton collisions
We also consider the data from two experiments measuring the DY production in proton-proton collisions at moderate energies: R209 [36] with √ S = 62 GeV and PHENIX [40] with √ S = 200 GeV. For R209 we apply a change of variables, whereas PHENIX is using the cross section Ed 3 σ/d 3 q given by (69). The theoretical results were compared with the R209 data in Fig. 4. CCFM-K provides very good description of data up to q T ∼ 4 GeV and slightly underestimate cross-section for higher q T while CSS overestimates the cross section at small q T and decreases rapidly at high values. For the PHENIX data shown in Fig 9, CCFM-K gives a better description than CSS, which overestimates the data at moderate values of q T . We note that as for fixed target experiments, the on-shell CCFM-K better describes the data then the off-shell approach.

Summary
Using the CCFM-K parton distributions and the partonic cross-section with on-shell and offshell matrix element, we analyse the transverse momentum spectra of the DY dileptons from all available low mass data. The overall description of these data is quite good, given the simplicity of the non-perturbative Gaussian input (33) with only one free parameter -the Gaussian width. We have chosen somewhat optimised value of this parameter for all experiments to show the potential of the CCFM-K approach in the description of the data for large dilepton transverse momentum q T .
However, for small q T we found less successful description, especially in low mass bins. This calls for an approach with more non-perturbative parameters in initial conditions for the CCFM-K evolution akin to the BNLY fit [67] of the non-perturbative form factor (68) in the CSS approach. This is justified since the CCFM-K evolution includes elements of the CSS resummation. In this sense our paper should be treated as a step towards unified description of the low mass DY data, where the CSS approach matched to the fixed order calculation experiences some troubles [9].
One should also note that the presented analysis is based on the LO matrix elements and the CCFM-K evolution equations in the single loop approximation. For these reasons, we decided to postpone the analysis with more complicated non-perturbative input to future studies with the NLO matrix elements and evolution equations. Finally, it is important to stress that the future analysis should also include the Tevatron and LHC data on the weak bosons production.

A Relation of CCFM-K to CSS
Following the method presented in [55], we show the CCFM-K resummation contains the soft gluon resummation of Collins, Soper and Sterman (CSS) [7]. In order to simplify the notation, we apply the Mellin transformf to both sides of eqs. (31) to obtain wheref S = 2N f i=1f i is the singlet quark distribution and the Mellin moments of the splitting functions read Notice that for b = 0 in (73) we obtain the DGLAP evolution equations. This fact motivates the following decomposition of the diagonal splitting functions Thus, for b = 0,P qq2 =P gg2 = 0, andP qq1 andP gg1 become the ordinary Altarelli-Parisi splitting functions. This is why we write (73) in the form where for simplicity we suppressed the arguments. We look for the solutions in the form where Inserting (81) to (80), we find We are interested in the approximate solution when the parameter b obeys the relation (35) written now as Q 0 1/b Q. In such a case, the powers of large logarithms, ln(Qb) and ln(1/Q 0 b), organize the calculations. We will show that the solution in such an approximation is given by the Mellin moments of the CSS formulas (36) where the coefficients Proof. The dominant contribution to the integrals with the Bessel function J 0 (u) comes from the region u < 1, therefore, we use the following approximations where θ is the Heaviside step function and c ∼ 1. The precise value of this parameter is important for numerical studies but for simplicity of this analysis we set c = 1. Thus, the quark exponent S q in (82) withP qq2 given by (77) is given by For qb < 1, the argument of the theta function is negative and S q = 0. Thus, we have From the first theta function q > 1/b Q 0 which sets the lower integration limit to 1/b. In this way, we avoid resummation of large logarithms ln(1/Q 0 b) which are shifted to the functions f i,g in (81) and need to be resummed separately. For 1/b > Q, the variable q > Q and S q = 0. Thus, S q ∼ θ(Qb − 1). Writing P qq in the form we find the result which agrees with the exponent in (84) where in the last line we neglected terms with subleading powers of logs ln(Qb) after the integration. A similar calculation for S g in (82) leads to the exponent in (85) The large logs ln(1/Q 0 b) are resummed using the DGLAP evolution equations. To show this, we write (83) in the integral form The first integral on the r.h.s. is given by where we used (76). Approximating we notice that in order to get the leading logs ln(1/Q 0 b) we have to assume that µb < 1. In such a case, the integration over z in (94) is not constrained but the integration over µ is limited to µ < 1/b Q. In this way Applying the same approximations to second integral in (93), we obtain From (91), in the integration region, µ < 1/b, we have S q (b, µ) = S g (b, µ) = 0. Therefore, we can set the exponent in I 2 equal one and (93) reads These are the DGLAP equation for the moments of quark distributions distribution evolved up to the scale 1/b. The analogous considerations the gluon distributions lead to the gluon counterpart of the DGLAP equations. It can easily be checked that for µb > 1, the theta function (95) imposes the constraint z > 1 − 1/(µb) which gives subleading logarithmic contributions which we neglect.
To make a connection with the collinear PDFs, we note that for the values of b 1/Q 0 which we consider, we obtain b/GeV 1 for Q 0 ∼ 1 GeV. Thus, to good approximation b ≈ 0 in the second argument of the functionsf i,g in (98). Choosing the initial conditions equal to the Mellin moments of the collinear PDFs we obtain the Mellin moments of collinear PDFs at the scale 1/b, due to the DGLAP evolution equations (98), This concludes the proof that (84) and (85) are the approximate solutions to the CCFM-K equations. As a final remark, the parameter c = 1 in (87) leads to the replacement 1/b → c/b in all the formulae above.  Figure 6: Transverse momentum dependence of DY cross-section: data from fixed target experiments E605 [33] (upper panels) and E866 [35] (lower panels) are compared with on-shell CCFM-K (blue solid), off-shell CCFM-K (blue dash-dotted) and CSS (red dashed) approaches.