The role of initial gluon emission in double $J/\psi$ production at central rapidities

We consider the process of double $J/\psi$ production in $pp$ collisions at the LHC in the framework of $k_T$-factorization approach. We focus on the gluon fragmentation mechanism which is related to multiple gluon emission in the initial state. The initial state emission is treated according to Catani-Ciafaloni-Fiorani-Marchesini evolution equation, and applies to both single and double parton scattering cases. We show the importance of fragmentation contributions to $J/\psi$ pair production and perform a comparison between theoretical predictions and the latest ATLAS data collected at $\sqrt{s} = 8$ TeV. We find that the effects of multiple gluon emission are essential for both single- and multiparton interaction processes. Finally, we highlight the problem of correct choice of the factorization scale with respect to numerical stability of the calculations and the consistency with non-collinear evolution equations.


Introduction
The production of J/ψ pairs at high energies serves as an important probe testing the quarkonia production mechanisms and their interpretation within nonrelativistic quantum chromodynamics (NRQCD) [1][2][3]. NRQCD provides a rigorous theoretical framework commonly used to describe the production and decay of heavy quark bound states. It implies a factorizable separation between perturbatively calculated short distance cross sections for the production of a heavy quark pair and its subsequent nonperturbative transition into a physical particle. The intermediate QQ state | 2S+1 L (a) J is characterized with its spin S, orbital angular momentum L, total angular momentum J, and color representation a. Its transition to a physical meson proceeds via soft gluon radiation and is described (parametrized) by the long-distance non-perturbative matrix elements (LDMEs), which obey certain hierarchy in powers of the relative heavy quark velocity v [1][2][3]. Combined with next-to-leading order (NLO) short-distance cross sections, NRQCD fits the LHC data on the prompt J/ψ, ψ and χ c transverse momentum distributions (see, for example, [4][5][6][7][8][9][10][11]). A long-standing challenge in explaining the polarization phenomena (the so-called "polarization puzzle") have been solved recently [12] (see also discussions [13][14][15]).
In the last few years, NRQCD has made significant progress in evaluating the prompt double J/ψ production. A complete leading-order (LO) calculation including color singlet (CS) and color octet (CO) terms is done [16]. Relativistic corrections to the J/ψ pair production are carried out [17]. Full NLO contribution to the CS mechanism is known [18], and partial tree-level NLO* contributions to the CS and CO mechanisms are calculated [19]. The latter are found to be essential for both low and large transverse momenta, as compared to the LO results. However, these predictions still suffer from sizeable discrepancies with the latest ATLAS data [20], especially at large transverse momentum p T (J/ψ, J/ψ), large invariant mass m(J/ψ, J/ψ), and large rapidity separation ∆y(J/ψ, J/ψ) in the J/ψ pairs, that motivates looking for additional mechanisms contributing to the double J/ψ events.
In our previous publication [21], we revealed a sizeable combinatorial contribution from multiple gluon radiation originating from the QCD evolution of the initial gluon cascade. This contribution can be efficiently taken into account using the Ciafaloni-Catani-Fiorani-Marchesini (CCFM) evolution equation [22][23][24][25]. Indeed, gluons emitted in a non-collinear evolution cascade have non-zero transverse momenta and give rise to physical J/ψ mesons via color octet fragmentation. The impact of such processes on double J/ψ production at forward rapidities has been investigated [21], and their importance for the associated Z/W ± + J/ψ production at the LHC has been pointed out [26].
It has been shown [21] that the gluon and quark fragmentation into charmonium states could especially play a crucial role in the kinematical region of large invariant masses m(J/ψ, J/ψ) and/or large rapidity separation ∆y(J/ψ, J/ψ) between the J/ψ mesons. Given the fact that this region is covered by the ATLAS experiment [20], we formulate our next goal. In our present study, we are going to investigate the effects of multiple gluon radiation with respect to double J/ψ production at central rapidities and to give a quantitative comparison of our predictions with the available ATLAS data [20]. The effect of multiple gluon radiation in DPS events is to be studied for the first time.
The outline of the paper is the following. In Section 2 we briefly describe the basic steps of our calculations. Section 3 is devoted to a discussion on the different choices of the factorization scale. Section 4 displays our numerical results. Our conclusions are summarised in Section 5. Figure 1: Examples of Feynman diagrams, contributing to the production of: a) J/ψ pair, subprocess (1); b) J/ψ and χ cJ pairs via CS mechanism, subprocesses (2) -(4).

The model
To preserve consistency with our previous studies [21,26], here we employ the k Tfactorization approach [27][28][29][30]. This approach is based on the Balitsky-Fadin-Kuraev-Lipatov (BFKL) [30,31] or Catani-Ciafaloni-Fiorani-Marchesini (CCFM) gluon evolution equations and has certain technical advantages in the ease of including higher-order pQCD radiative corrections (namely, the leading-logarithm part of NLO + NNLO + ... terms corresponding to real gluon emissions) in the form of transverse momentum dependent (TMD, or unintegrated) gluon density in a proton. It can be used as a convenient alternative to explicit higher-order pQCD calculations. A detailed description of this approach can be found, for example, in review [32].
First, we consider the O(α 4 s ) off-shell gluon-gluon fusion subprocess (where the initial gluons have nonzero transverse momenta and are then off-shell) which represents the leading order CS contribution: (1) , together with the feeddown contributions from excited states ψ decaying into J/ψ mesons. These subprocesses are formally suppressed by extra powers of α s , or smaller values of P -wave mesonic wave functions, or χ c → J/ψ + γ decay branchings, but have an important kinematic property of filling the region of large invariant masses and large rapidity separation in J/ψ pairs. The typical Feynman diagrams are presented in Fig. 1. The production amplitudes of (1) -(4) contain spin and color projection operators [33][34][35][36] which guarantee the proper quantum numbers of final state charmonia. In accordance with the k T -factorization prescription, initial gluons have non-zero transverse four-momenta k 2 1T = −k 2 1T = 0 and k 2 2T = −k 2 2T = 0, and their polarization vectors have an admixture of longitudinal component. The summation over polarizations is carried out with [27][28][29][30]. This expression converges to the ordinary g µν in the collinear limit k T → 0 after averaging over the azimuthal angle. The gauge invariant expressions for all these amplitudes have been obtained earlier [37] and implemented into the Monte-Carlo event generator pegasus [38]. Now we turn to another class of processes which constitute the key topic of our study, namely, the fragmentation contributions. The fragmentation approach is known to be valid at high transverse momenta p T m ψ . The relevant fragmentation functions (FF) D H a (z, µ 2 ) describing the transition of parton a into charmonium state H through a number of intermediate QQ states can be presented as a series: where n labels the intermediate (CS or CO) state, and O H [n] are the corresponding LDMEs. In this way, the single charmonia production cross section can be written as where p (g * ) , p (c) and p are the momenta of the gluon, charmed quark and outgoing charmonium state H, respectively. We take into consideration the channels g → cc[ 3 S 1 ], g → cc[ 3 P     J ] + c contributing to the P -wave charmonia (χ cJ mesons with J = 1, 2). The χ c0 channel is neglected because of low branching fraction to J/ψ. The present list is more complete in comparison with our previous paper [21] where we refer solely to the g → cc[ 3 S The key point of our consideration is that any hard subprocess (giving or not giving rise to a J/ψ meson) is always accompanied by a number of gluons radiated during Figure 2: Fragmentation contributions to the J/ψ pair production that takes into account multiple gluon emissions for subprocesses: a) g * + g * → g * , b) g * + g * → c +c, c) the non-collinear QCD evolution, and these gluons can fragment into additional J/ψ particles. We calculate the corresponding contributions by collecting all possible parton fragmentation combinations. At high energies, the QCD evolution of gluon cascade can be described by the CCFM equation. This equation smoothly interpolates between the small-x BFKL gluon dynamics and high-x DGLAP one, and, therefore, provides us with a suitable tool for our phenomenological study. The numerical calculations split in two steps. First, we simulate the perturbative production of gluons, quarks, and charm pairs in the corresponding off-shell gluon-gluon fusion subprocesses and then reconstruct the CCFM gluon evolution ladder using the cascade Monte Carlo generator [39]. After that, one can collect J/ψ pairs by looking over all possible combinations of mesons originating from charmed pairs, charmed quarks and gluons (including those formed in the evolution cascade). The combinatorics is rather large, because any fragmenting parton can be paired with any other parton.
For the initial hard subprocesses, we consider g * + g * → g * , g * + g * → c +c, g * + g * → q +q, where q denotes light quarks. The typical diagrams of subprocesses (7) with reconstructed gluon ladders and possible channels of fragmentation are presented in Fig. 2. The CO channel g * + g * → cc[ 3 S 1 ] is excluded from subprocesses (7) in order to avoid double counting with the fragmentation mechanism g * +g * → g * → cc[ 3 S 1 ]. By the way, we have found that these two subprocesses perfectly match one another in the ATLAS kinematic range.
The expressions for fragmentation functions at the initial scale µ 2 0 = m 2 ψ can be found [40]. For the fragmentation g → cc[ 3 P ] + g we use the expression derived very recently [41], with the mass of the emitted gluon (considered as regularization parameter) m g = m c = 1.5 GeV. This expression is positive-definite, smooth and vanishes at the endpoints z = 0 and z = 1. The shapes of fragmentation functions are modified by the final state gluon radiation; these effects can be described in proper way with the DGLAP evolution equation: where P ab are the usual LO DGLAP splitting functions. According to the non-relativistic QCD approximation, we set the charmed quark mass to m c = m ψ /2 and then solve the DGLAP equation (8) numerically with the proper LDME's. The last contribution taken into consideration refers to the double parton scattering (DPS) mechanism. According to the standard factorization formula [42,43], this contribution can be presented in a simple form: where the factor 1/2 prevents double counting between identical particles. The "effective cross section" σ eff is a normalization factor which encodes all "DPS unknowns" into a single parameter and represents the effective transverse overlap of partonic interactions that constitute the DPS process. In general, it can be regarded as free parameter which should be extracted from the data. The decomposition (9) of the DPS cross section into two individual single parton scattering (SPS) factors without correlation and interference between them is acceptable for ATLAS kinematics region. The inclusive cross sections σ(pp → J/ψ + X) involved in (9) are calculated in the k T -factorization approach supplemented by the NRQCD formalism in a standard way (see, for example, [44] and references therein). When calculating the DPS cross section, we also take into account all the accompanying fragmentation contributions (including all possible combinations of radiated partons, see Fig. 3).
In the numerical calculations we use TMD gluon densities in a proton obtained from a numerical solution of CCFM evolution equation, namely, JH'2013 set 1 and JH'2013 set 2 [45]. The input parameters of JH'2013 set 1 gluon distribution have been fitted to the proton structure function F 2 (x, Q 2 ), whereas the input parameters of JH'2013 set 2 gluon Figure 4: The effect of scale variations in TMDs. By changing the factorization scale, we jump between the red and blue lines; in the region µ F k T the difference may be as large as several orders of magnitude.
were fitted to the both structure functions F 2 (x, Q 2 ) and F c 2 (x, Q 2 ). According to [45], we use the two-loop formula for the QCD coupling α s with n f = 4 active quark flavors and Λ QCD = 200 MeV.
The charmonia LDMEs have been defined earlier [44] from a global fit to the LHC data. Using these LDMEs, we reproduce all of the available data on charmonia production at the LHC conditions. A comprehensive information on the fitting procedure can be found in [44]. The masses and the branching fractions of all particles involved into calculations are taken from [46]. The DPS effective cross section is chosen as σ eff = 13.8 mb, which was extracted from a fit to the latest LHCb data on the double J/ψ production (see [21] for more details). This value is very close to the generally accepted value σ eff = 15 mb [47].

Choice of factorization scale
The choice of the factorization scale is a delicate, though very important issue. In the CCFM evolution equation the factorization scale µ F is defined as µ 2 F =ŝ + Q 2 T , wherê s = (k 1 + k 2 ) 2 is the invariant energy of partonic subprocess and Q T is the net transverse momentum of the initial gluon pair. It looks rather natural to use the same definition of µ F throughout all calculations, but this may cause problems in some cases.
A peculiar property of the partonic subprocess (1) is that its cross section drops sharply with increasingŝ, so thatŝ is typically not far from threshold (and, consequently, is small). At the same time, the k T -spectra of the initial gluons are exponentially broad. This means that in an arbitrary pair of gluons, the transverse momentum |k T | of one of the gluons is much larger than that of the other. As a result, the expression µ 2 However, such a definition is self-contradictory. Indeed, the gluon density f g (x, k 2 T , µ 2 F ) must describe the probability distribution for k 2 T at any given µ F . Assume, we set some µ F and generate a random gluon transverse momentum with the probability given by f g (x, k 2 T , µ 2 F ). But then we come to a conflict with setting µ 2 F = k 2 T which requires µ F to be different from what was set originally. In fact, the condition µ 2 F ∝ k 2 T means that the full area of µ 2 F and k 2 T is not accessible, but we are only restricted to a one-dimensional trajectory tracing the functional dependence of µ 2 F on k 2 T . Moreover, this trajectory lies entirely on a steep slope, as the gluon  Double differential distributions in the gluon transverse momenta, dσ/d log 10 |k 1T |d log 10 |k 2T |; toy calculations with differently chosen factorization scales and different toy matrix elements |M| 2 for the partonic subprocess. Left column corresponds to µ 2 F =ŝ + Q 2 T , right column -µ 2 F =ŝ/4. Upper row, |M| 2 ∝ 1; middle row, |M| 2 ∝ 1/ŝ 2 ; lower row, |M| 2 ∝ 1/ŝ 4 . Calculations were performed with JH'2013 set 1.   ŝ, that is, when eitherŝ is small or p T is large. densities steeply change around k 2 T = µ 2 F from high values at k 2 T < µ 2 F to low values at k 2 T > µ 2 F . This makes the calculation very unstable with respect to even small variations in the factorization scale, see Fig.4.
The above properties are further illustrated in Figs. 5 and 6. Shown in Fig. 5 is the double differential distribution in the gluon transverse momenta dσ/d log 10 |k 1T |d log 10 |k 2T | obtained from toy calculations with differently chosen factorization scales and with different toy matrix elements |M| 2 for the partonic subprocess. Here we see that the choice µ 2 F =ŝ/4 favors moderate gluon transverse momenta and the shape of the distribution is insensitive to the properties of the matrix element. On the contrary, the choice µ 2 F =ŝ + Q 2 T makes the k 2 T values large and highly unequal, and the |k T | distribution is sensitive to the properties of |M| 2 .
The effect of numerical instability is shown in Fig. 6. With sharper matrix elements, the leading role in µ F transits fromŝ to Q 2 T thus bringing us to an "unsafe" regime µ 2 F k 2 T . It may be worth noting that the problem is rather general. Whatever the behavior of the matrix elements is, there always exists a kinematic region where p 2 T ŝ. The case of J/ψ pair production was only "lucky" to reveal the problem at smaller energies and smaller p T .
In Fig. 7 we extend our exercises from toy to real matrix elements and inspect the behavior of subprocess (1). As is expected, the region of the biggest numerical sensitivity to the choice of µ F is the region of the smallest invariant masses and, respectively, the region of the smallest rapidity difference ∆y(J/ψ, J/ψ). The divergence between the different predictions is huge. At the same time, the curves obtained with µ 2 F = (ŝ + Q 2 T )/4 are notably close to the experimental points as it is demonstrated at the Section 4. It is necessary to note that these predictions were obtained with the TMD gluon densities JH'2013 set 1 muf−, which use the prefactor of 1/4 in the definition of µ 2 F in the CCFM evolution and in the fitting procedure.
On the one hand, the property that the fullrange gluonic phase space degenerates into a one-dimensional line µ 2 F ∝ k 2 T can be regarded as pathological, though inevitable with increasing p T . It leads to extraordinary sensitivity of the results to the choice of µ F . Since this kind of bad behavior is only seen at certain (rather special) kinematic conditions (subprocesses with very smallŝ), one can try to redefine the scale µ 2 F in this particular case. The definition of µ 2 F used in the TMD fitting procedure is related to the maximum angle between the two quarks formed in the subprocess γ * g * → qq [22][23][24][25] and follows from the angular ordering condition. Processes different in their topology may allow different forms of µ 2 F . On the other hand, one can argue that the definition of µ 2 F must preserve strict consistency between the fitting procedure (from which the gluon densities are obtained) and the actual calculation (to which the gluon densities are applied) and, therefore, have the form µ F ∝ŝ + Q 2 T (probably with some numerical prefactor). The property that the parameter space (µ 2 F versus k 2 T ) is one-dimensional rather than two-dimensional is not dangerous by its own (and is similar to collinear factorization). The choice of µ F is not free, but is part of the proposed TMD parametrization. We simply have to take the gluon densities obtained from the fit (with the given form of µ F ) as they are and directly substitute them into actual calculations.
An interesting parallel can be seen with collinear calculations at the NLO [48]. As one can see from Figs. 6 and 7 in that paper, the integral size of the NLO contribution exceeds the LO result by a significant (huge) factor. The difference between the NLO and LO predictions is most dramatic in the region of small J/ψ invariant masses and, respectively, small rapidity separation. In the k T -factorizaion approach, the NLO contributions are absorbed into the evolution of TMD gluon densities, where the parameter µ F regulates the gluon emission. So, we find it not surprising that the effect of µ F is largest in the same region where the NLO contribution is largest, i.e. at small m(J/ψ, J/ψ) and small ∆y(J/ψ, J/ψ). This has to be compared with our Fig. 7. Were the authors of [48] had calculated the p T (J/ψ, J/ψ) distributions, they would undoubtedly observe the same discrepancy between NLO and LO as we do in Fig. 7 bottom right. We had, in due time, reproduced all of the results [48] with our k T -factorizaion technique.
We would not, however, deny the fact that the numerical instability at µ 2 indicates that we are approaching the applicability limits of the k T -factorization (noncollinear parton evolution). The problem of extraordinary sensitivity of the results to the choice of µ F has probably been first detected in [49], but left without further attention and analysis.

Numerical results
In this section, we present the results of our calculations and perform a comparison with the recent ATLAS data collected at √ s = 8 TeV [20]. The ATLAS Collaboration has measured the differential cross sections of prompt J/ψ pair production as functions of the transverse momentum p T (J/ψ, J/ψ) and the invariant mass m(J/ψ, J/ψ) of the J/ψ pair, of the rapidity separation ∆y(J/ψ, J/ψ) and of the azimuthal angle difference ∆φ(J/ψ, J/ψ) between the two J/ψ mesons. The following selection criteria were applied: p T (J/ψ) > 8.5 GeV, |y(J/ψ)| < 2.1, |η(µ)| < 2.3 for all muons, p T (µ) > 4 GeV for the muons from the triggered J/ψ meson, and p T (µ) > 2.5 GeV for the muons from the other J/ψ meson. Also presented were the differential cross sections as functions of the subleading J/ψ transverse momentum p T (J/ψ 2 ), of the J/ψ pair transverse momentum, and of the J/ψ pair invariant mass, with the selection criteria p T (J/ψ) > 8.5 GeV, |y(J/ψ)| < 2.1 for each J/ψ, without requirements on the muons in the final state. These measurements were performed in the central |y(J/ψ 2 )| < 1.05 and forward 1.05 < |y(J/ψ 2 )| < 2.1 rapidity intervals. We have implemented the same selection criteria in our calculations.
The calculations are performed with the following setting of hard scales. For the CS contributions (1)-(4) we take µ 2 F = 1 4 (ŝ+Q 2 T ) and µ 2 R = 1 4 (m 2 ψ + 1 2 (p 2 1T +p 2 2T )) with p 1T and p 2T being the transverse momenta of produced particles. For these processes we use the gluon densities JH'2013 set 1 muf− and JH'2013 set 2 muf− obtained with µ 2 F = 1 4 (ŝ+Q 2 T ). So, the numerical prefactor 1/4 in the definition of µ F holds exact correspondence with the fitting procedure [45]; and, in comparison with the choice µ 2 F = (ŝ + Q 2 T ), it softens to some extent the scale dependence of the results.
For the fragmentation contributions and DPS we set for subprocesses with outgoing off-shell gluon or bound charmed pair; and ) for subprocesses with outgoing unbound charmed quarks. The factorization concept implies that the fragmentation scale can only depend on the parameters of the fragmented parton and must not depend on the hard interaction scale. So, we set µ 2 frag = m 2 Q + p 2 T where m Q and p T are the mass and transverse momentum of the fragmented parton.
Our numerical results are shown in Figs. 8 and 9. One can see that the predictions obtained with the JH'2013 set 1 gluon density (used as the default choice) are rather close to the data and can be said compatible with them within theoretical uncertainties (shaded orange bands) for the majority of the measured distributions. As usual, the theoretical uncertainties were estimated by varying the renormalization scales around their default values by a factor of 2.   To highlight the role of DPS mechanism and the combinatorial effects of multiple gluon and/or quark radiation in the initial state, we separately show the relevant contributions. The blue dashed histograms in Figs. 8 and 9 correspond to the estimated DPS terms, whereas the blue solid histograms (labeled as "DPS + Fragm.") represent the sum of DPS and fragmentation contributions. As one can see, the effect of multiple parton radiation is very essential in the region of high invariant mass m(J/ψ, J/ψ) and large rapidity separation ∆y(J/ψ, J/ψ). An accurate account of these contributions is extremely important for ATLAS data. The importance of fragmentation terms have been already demonstrated earlier [21], where, however, only the g * → cc[ 3 S 1 ] transition has been taken into account. Our present calculations include a much larger number of fragmentation channels; and the feeddown contributions from the χ c and ψ decays are also taken into account. Finally, we show the contributions from subprocesses (2) -(4) (blue solid curves), although we find these subprocesses not playing a considerable role for any observable.
To investigate the sensitivity of our results to the choice of TMD gluon density, we have repeated our calculations with JH'2013 set 2 distribution. We find that the latter predictions (red histograms in Figs. 8 and 9) are in good agreement with the ones obtained with the default JH'2013 set 1 gluon density for the majority of observables. There are some discrepancies in the region of high p T (J/ψ, J/ψ) and p T (J/ψ 2 ). They can be promptly attributed to the different behavior of TMD gluon distributions at high transverse momenta k T .
Next, we discuss the role of multiple gluon radiation in the DPS contributions. As it was said above, our calculation scheme implies that the direct CO subprocesses g * + g * → cc[ 3 S (8) 1 ] in the both hard interaction blocks have to be replaced with the g * + g * → g * and/or g * +g * → c+c subprocesses accompanied by gluon evolution ladders (reconstructed according to the CCFM equation) with subsequent fragmentation of all emitted partons into charmonia. This brings additional contributions to the DPS production cross section. To show the role of these terms in more detail, we separately consider three different sources 1 . The first of them contains only [ 3 S J ] combinations, where the symbols g and c denote fragmentation contributions coming from the g * + g * → g * and g * + g * → c +c subprocesses (with multiple gluon radiation taken into account in both cases). The sum of the first and the third contributions represents the DPS cross section in our scheme.
A comparison between the different contributions obtained within the DPS framework is displayed in Fig. 10. It is clearly seen that the yield from modified mechanism (i.e., the third source, orange lines) is approximately three times greater than that from the conventional mechanism (the second source, blue lines). Eventually, it nearly doubles the total estimated DPS contribution. This constitutes in our view a remarkable result. Thus, an accurate treatment of multiple gluon emission in the initial state is indeed very important for evaluating the DPS cross sections. In fact, these additional contributions could also play a role for a number of other processes related to the multiparton interactions. This can, in turn, motivate some revisions of the effective DPS cross section σ eff extracted from the experimental data.   [8] part without 3S ψ + J/ ψ DPS J/ fragm. part ψ + J/ ψ DPS J/ 1 [8] part with 3S ψ + J/ ψ DPS J/

Conclusion
In the present study, we have addressed a challenging problem of prompt double J/ψ production in pp collisions at the LHC conditions. In addition to the conventional production mechanisms mentioned earlier in the literature, we take into account the effects of multiple gluon radiation in the initial state followed by gluon fragmentation into J/ψ mesons. The evolution of the radiated gluon cascade is described within the k Tfactorization approach, with making use of CCFM equation. We demonstrate the importance of these new contributions both in single and double parton scattering processes.
We paid much attention to inspecting the influence of the factorization scale µ F on the numerical errors and to studying the kinematic conditions under which the computations become unstable. We find that the dangerous region is p 2 T ŝ, that corresponds to µ 2 F k 2 T . The computations can be made stable by adopting either of two alternatives. One can adopt just a different definition of the scale µ F , so that to make it far from the instability region (basically, by taking µ F independent on k T ). As an alternative, one can consider the definition µ F ∝ŝ + Q 2 T as a built-in property of the TMD gluon density. The latter choice preserves the correspondence between the TMD fitting procedure and the actual calculations where the obtained TMD gluon distribution is used.
In the framework of our second hypothesis, we have performed a numerical comparison with the latest measurements reported by the ATLAS collaboration for √ s = 8 TeV. Having the fragmentation mechanism employed in the theory, we greatly reduce the discrepancy between the predictions and the data. This has immediate impact on the value of the effective cross section σ eff that parametrises the DPS contribution.