Heavy Flavour Production at Tevatron and Parton Shower Effects

We present hadron-level predictions from the Monte Carlo generator Cascade and numerical calculations of charm and beauty production at the Fermilab Tevatron within the framework of the $k_T$-factorization QCD approach. Our consideration is based on the CCFM-evolved unintegrated gluon densities in a proton. The performed analysis covers the total and differential cross sections of open charm and beauty quarks, $B$ and $D$ mesons (or rather muons from their semileptonic decays) and the total and differential cross sections of $b \bar b$ di-jet hadroproduction. We study the theoretical uncertainties of our calculations and investigate the effects coming from parton showers in initial and final states. Our predictions are compared with the recent experimental data taken by the D0 and CDF collaborations. Special attention is put on the specific angular correlations between the final-state particles. We demonstrate that the final state parton shower plays a crucial role in the description of such observables. The decorrelated part of angular separations can be fully described, if the process $gg^*\rightarrow gg$ is included.


Introduction
Charm and beauty production at high energies is subject of intense studies from both theoretical and experimental points of view [1][2][3][4][5][6]. From the theoretical point, these processes provide an opportunity to test the different predictions based on Quantum Chromodynamics (QCD) since the dominant production mechanism at high energies (i.e. small x) is believed to be quark pair production through the gluon-gluon fusion subprocess. At present, the problem of description of the charm and beauty production at the Tevatron within QCD is not fully solved. So, the difference between the D0 and CDF measurements [1][2][3][4] of the b-quark and B-meson production cross sections and the calculations [7] performed in the framework of fixed-order next-to-leading logarithm scheme (FONLL) is about of a factor of 1.7. The central FONLL predictions [8] also lie below the data on the D 0 , D + , D * + and D + s cross sections which has been measured [5] by the CDF collaboration as a functions of their transverse momenta. Recently the calculations in the general-mass variable-flavour-number scheme (GM-VFNS) were performed [9,10] for the transverse momentum distribution of B and D mesons and have been found to be consistent with the data [4,5]. However, there is still no GM-VFNS predictions for other measured quantities (like distributions in rapidity and azimuthal angle difference between the momenta of final mesons).
Heavy flavour production has been considered [11][12][13][14][15][16] also in the framework of the k Tfactorization QCD approach [17]. This approach is based on the Balitsky-Fadin-Kuraev-Lipatov (BFKL) [18] or Ciafaloni-Catani-Fiorani-Marchesini (CCFM) [19] equations for the non-collinear gluon evolution in a proton and gives the possibility to take into account large logarithmic terms (proportional to ln 1/x and, in the case of CCFM, also to ln 1/(1 − x)). A detailed description and discussion of the k T -factorization approach can be found in [20]. A reasonable agreement between the k T -factorization predictions and the Tevatron data on the heavy flavour production has been found in [12][13][14][15][16]. It was demonstrated [11,15,16] that studying the specific angular correlations between the transverse momenta of produced quarks can give an unique information about the non-collinear gluon evolution in a proton since taking into account the non-vanishing initial gluon transverse momentum k T in the k T -factorization approach leads to the violation of back-to-back kinematics event at leading order.
Measurements of bb di-jet total and differential cross sections (as a functions of the leading jet transverse energy E T and the di-jet invariant mass) and azimuthal angle correlations between two b-jets have been performed by the CDF collaboration [21]. These measurements significantly extend the energy range investigated by previous analyses [1][2][3][4][5]. Studying the b-jet production is specially interesting since there are no additional assumptions on the fragmentation of the beauty quark into the B meson. The CDF collaboration has reported the preliminary data [6] on the charm pair production, where the D 0 , D * − pair cross section and the D + , D * − pair cross section as a function of the azimuthal angle between two charmed mesons have been meausred. The angular correlations in B-meson production have been also measured [1][2][3] by the D0 and CDF collaborations.
The main goal of present paper is to give a systematic analysis of all available experimental data [1][2][3][4][5][6]21] on the heavy flavour production at the Tevatron in the framework of the k T -factorization formalism. We produce the relevant numerical calculations in two ways. First, we will perform analytical parton-level calculations (labeled as LZ in the following) similar to that done in [15,16]. In these calculations we will use the CCFM-evolved gluon densities [23] as default sets. The measured cross sections of heavy quark production will be compared with the predictions of Monte Carlo event generator Cascade [24]. Cascade is a full hadron level Monte Carlo event generator for ep, γp, pp 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 1 . In this way we will investigate the influence of parton showers in initial and final states for the description of the data. We will study the possible sources of theoretical uncertainties of our predictions (i.e. uncertainties connected with the gluon evolution scheme, heavy quark mass, hard scale of partonic subprocess and the heavy quark fragmentation functions). To investigate the dependence of our predictions on the non-collinear evolution scheme we will apply the unintegrated gluon densities derived from the usual (DGLAP-evolved) parton distributions (in the Kimber-Martin-Ryskin (KMR) [26] approximation). Our special goal is to study specific kinematic properties of the final heavy quark-antiquark pair which are strongly related to the non-zero initial gluon transverse momentum.
The outline of our paper is following. In Section 2 we recall shortly the basic formulas of the k T -factorization approach with a brief review of calculation steps. In Section 3 we present the numerical results of our calculations and a discussion. Section 4 contains our conclusions.

Theoretical framework
The main formulas have been obtained previously in [14][15][16]. Here we only recall some of them. The cross section of heavy quark hadroproduction at high energies in the k Tfactorization approach is calculated as a convolution of the off-shell (i.e. k T -dependent) partonic cross sectionσ and the unintegrated gluon distributions in a proton. It can be presented in the following form: where A(x, k 2 T , µ 2 ) is the unintegrated gluon distribution in a proton, |M(g * g * → QQ)| 2 is the off-shell (i.e. depending on the initial gluon virtualities k 2 1T and k 2 2T ) matrix element squared and averaged over initial gluon polarizations and colors, and s is the total center-ofmass energy. The produced heavy quark Q and anti-quarkQ have the transverse momenta p 1T and p 2T and the center-of-mass rapidities y 1 and y 2 . The initial off-shell gluons have a fraction x 1 and x 2 of the parent protons longitudinal momenta, non-zero transverse momenta k 1T and k 2T (k 2 1T = −k 2 1T = 0, k 2 2T = −k 2 2T = 0) and azimuthal angles φ 1 and φ 2 . The analytic expression for the |M(g * g * → QQ)| 2 can be found, for example, in [14,17].
The unintegrated gluon distribution in a proton A(x, k 2 T , µ 2 ) in (1) can be obtained from the analytical or numerical solution of the BFKL or CCFM evolution equations. In the numerical calculations we have tested a few different sets. First of them (A0) was obtained in [23] from the CCFM evolution equation. The initial (starting) distribution A 0 (x, k 2 T , µ 2 0 ) has been used in the following form: where all input parameters have been fitted to describe the proton structure function F 2 (x, Q 2 ). An equally good fit can be obtained using different values for the soft cut and a different value for the width of the intrinsic k T distribution (the CCFM set B0). A reasonable description of the F 2 data can be achieved [23] by both these sets.
To evaluate the unintegrated gluon densities in a proton A(x, k 2 T , µ 2 ) we apply also the KMR approach [26]. The KMR approach is a formalism to construct the unintegrated parton (quark and gluon) distributions from the known conventional parton distributions xa(x, µ 2 ), where a = g or a = q. In this scheme, the unintegrated gluon distribution is given by the expression where P ab (z) are the usual unregulated leading order DGLAP splitting functions, q(x, µ 2 ) and g(x, µ 2 ) are the conventional quark and gluon densities 2 and T g (k 2 T , µ 2 ) is the Sudakov form factor. The theta function Θ(∆−z) implies the angular-ordering constraint ∆ = µ/(µ+|k T |) specifically to the last evolution step to regulate the soft gluon singularities.
In Fig. 1 we plot the CCFM set A0 (as a solid lines), the CCFM set B0 (as a dashed lines) and KMR (as a dotted lines) unintegrated gluon densities at probing scale µ 2 = 100 GeV 2 as a function of k 2 T for different values of x. One can see that the shapes of gluon densities under consideration are very different from each other In the following we will study the possible manifestations in the total and differential heavy flavour cross sections.
In our analytical calculations (LZ) the multidimensional integrations in eq.(1) have been performed by the means of Monte Carlo technique, using the routine vegas [30]. The full C++ code is available from the authors on request 3 .

Numerical results
We now are in a position to present our numerical results. First we describe our input and the kinematic conditions. After we fixed the unintegrated gluon distributions, the cross section (1) depends on the renormalization and factorization scales µ R and µ F . In the numerical calculations we set µ 2

Inclusive charm and beauty production
We begin the discussion by studying the role of non-zero gluon transverse momentum k T in the off-shell matrix elements involved in (1). In Fig. 2 we plot the differential cross section for cc and bb pair production as a function of transverse momentum p Q T , rapidity y Q and the azimuthal angle difference between the transverse momenta of produced quarks ∆φ QQ at √ s = 1960 GeV. The solid histograms correspond to the results obtained according to the master formula in eq.(1). The dotted histograms are obtained by using the same formula but without virtualities of the incoming gluons in partonic amplitude and with the additional requirement k 2 1,2 T < µ 2 R . There are no cuts applied on the phase space of the produced quarks. As it was expected, in the back-to-back region ∆φ QQ ∼ π both results coincide with each other. However, we find that a sizeable effect appears at low ∆φ QQ . Therefore the non-zero gluon transverse momentum in the hard matrix element is important for the description of the data at low and mediate ∆φ QQ . This effect is more significant for charm production due to smaller x. We have checked that the predictions between the LZ and Cascade calculations agree well at parton level. Now we turn to the transverse momentum distributions of charm and beauty production. In the case of inclusive b-quark production, the transverse momentum distribution has been measured [1] by the D0 collaboration and has been presented in form of integrated cross section (as a function of b-quark minimal transverse momentum p b T min ) at the total pp energy √ s = 1800 GeV for |y b | < 1 (note that there is no cut on the rapidity ofb). Our predictions are shown in Fig. 3 and compared to the data. We find good agreement between the LZ and Cascade predictions. We find a good description of the data when using the CCFM-evolved (A0) gluon distribution. The results obtained by using the KMR gluon density are rather close to the NLO pQCD ones [31] (not shown) but lie below the data. The difference between the CCFM and KMR predictions comes from the different behaviour of these gluon densities (see Fig. 1) which is due to absence of small-x resummation in the KMR distributions.
The CDF collaboration has measured [4,5] the transverse momentum distributions of B + and several D mesons (namely, D 0 , D + , D * + and D + s ) with |y| < 1 (where y is the meson rapidities in the center-of-mass frame) at √ s = 1960 GeV. Our predictions are shown in Figs. 4 -6 in comparison with the data [4,5]. The fragmentation of the charm and beauty quarks into a B and D mesons is described with the Peterson fragmentation function [32] with ǫ b = 0.006 and ǫ c = 0.06. According to [33,34], the following branching fractions are used: The observed difference between the LZ and Cascade predictions is due to the missing parton shower effects in the LZ calculations. We address this point in more detail in Section 3.3.
Since the predicted meson transverse momentum distributions depend on the quark-tohadron fragmentation, we have repeated our calculations for B + and D 0 mesons with the shifted values of the Peterson shape parameter ǫ, namely ǫ b = 0.003 and ǫ c = 0.03. These values are also often used in the NLO pQCD calculations. Additionally, we have applied the non-perturbative fragmentation functions which have been proposed in [7,8,35] and which have been used in the FONLL calculations. The input parameters were determined [8,35] by a fit to LEP data. The results of our calculations are shown in Fig. 7. We find that the predicted cross sections (in the considered p T region) are larger for smaller values of parameter ǫ or if the fragmentation function from [7,8,35] is used. However, the typical dependence of numerical predictions on the fragmentation scheme is much smaller than the dependence on the unintegrated gluon density.
The D0 experiment has measured muons originating from the semileptonic decays of b-quarks at √ s = 1800 GeV [1] for 4 < p µ T < 25 GeV, |η µ | < 0.8 (for both muons) and 6 < m µµ < 35 GeV, where η µ is the muon pseudo-rapidity and m µµ is the invariant mass of the produced muon pair. In [2] the measurements are extended to the forward muon rapidity region, namely 2.4 < |y µ | < 3.2. To produce muons from b-quarks in the LZ calculations, we first convert b-quarks into B-mesons (using the Peterson fragmentation function with default value ǫ b = 0.006) and then simulate their semileptonic decay according to the standard electroweak theory. The branching of b → µ as well as the cascade decay b → c → µ are taken into account with the branching fraction taken from [34]. The predictions of the LZ and Cascade calculations are shown in Figs. 8 and 9. We find that our predictions with both CCFM-evolved unintegrated gluon densities describe the experimental data for both the transverse momentum and rapidity distributions of muons reasonably well.
The calculated total cross sections of the b-quarks, B and D mesons and their decay muons compared to the CDF experimental data [4,5] are listed in Table 1. In Table 2 the systematic uncertainties of our calculations are summarized. To estimate the uncertainty coming from the renormalization scale µ R , we used the CCFM set A0+ and A0− instead of the default density function A0 in Cascade. These two sets represent a variation of the scale used in α s in the off-shell matrix element. The A0+ stands for a variation of 2 µ R , while set A0− reflects µ R /2. In all heavy flavour analyses studied here, we observe a deviation of roughly +10% for set A0+. The uncertainty coming from set A0− is generally smaller, but still positive. The dependence on the heavy flavour masses is investigated as well using We turn now to the investigation of the angular correlations between the produced particles. As it was mentioned above, such correlations have been not studied yet in the framework of GM-VFNS scheme (which has been applied for transverse momentum distributions of the charm and beauty mesons only [9,10]). Experimental data on the azimuthal correlations in charm and beauty production come from both the CDF and D0 collaborations. In the case of b-quark production, CDF data [3] refer to the B,B azimuthal angle distribution measured at |y| < 1, p T (B) > 14 GeV, p T (B) > 7.5 GeV and √ s = 1800 GeV. The D0 data [1] refer to the muon-muon correlation measured in the region of 4 < p µ T < 25 GeV, |η µ | < 0.8 and 6 < m µµ < 35 GeV at the same energy. In the case of charm production, the CDF collaboration has presented the measurement [6] of D 0 , D * − and D + , D * − pair cross section as a function of angle between the two charm mesons for the kinematic range |y| < 1, 5.5 < p T (D 0 ) < 20 GeV, 5.5 < p T (D * − ) < 20 GeV and 7 < p T (D + ) < 20 GeV. Our predictions are shown in Figs. 10 and 11 in comparison with the data [1,3,6]. We observe that the predicted shapes of azimuthal angle distributions are very different for different unintegrated gluon density functions. This is in a contrast to the cross sections as a function of transverse momenta or rapidities where all unintegrated gluon densities gave a similar behaviour. We can conclude that the cross section as a function of ∆φ is very sensitive to the details of the non-collinear gluon evolution in a proton. As it was pointed out in [11,15,16], such observ-  [4,5] was obtained in the central rapidity region (|y| < 1) at 6 < p T < 25 GeV and 5.5 < p T < 20 GeV in the case of B and D mesons, respectively.
Source   [4,5] was obtained in the central rapidity region (|y| < 1) at 6 < p T < 25 GeV and 5.5 < p T < 20 GeV in the case of B and D mesons, respectively.
ables can serve as an important and crucial test discriminating the different approaches of the small-x physics. The CCFM-evolved gluon densities overshoot the data at ∆φ ∼ π and tends to underestimate them at ∆φ ∼ 0. We observe that the peak at ∆φ → 0 is not at all reproduced by the CCFM unintegrated gluon densities; we address this point in more detail in Section 3.3 where we discuss the the process gg * → gg. However the KMR gluon density has obviously a very different k T distribution (see Fig. 1) and therefore provides a better description.

bb di-jet production
Recently, the experimental study of the bb di-jet production in pp collisions at √ s = 1960 GeV has been presented [21,22]. The total cross section has been measured in the kinematical region defined by |y 1 | < 1.2, |y 2 | < 1.2, E 1T > 35 GeV and E 2T > 32 GeV. The differential cross sections as a functions of the leading jet transverse energy and of the di-jet invariant mass have been measured with the azimuthal angle correlation between the two jets. The b-jets are reconstructed with the JETCLU cone algorithm with radius R > 0.4. Our predictions are shown in Figs. 12 -14. We observe that at large E T and at large dijet invariant masses M the predictions fall below the measurement. Note, however, that in this kinematical region the quark induced subprocesses (such as qq → bb) becomes important but are not taken into account here. At small and moderate E T and M, where gluon induced subprocesses play the leading role, the overall agreement of our predictions and the data is quite good. The measured ∆φ distribution has a significantly different tail at small ∆φ compared to the predictions; this will be addressed in more detail in the next section.

Parton shower effects and additional gluon processes
As shown above, the full set of the data on the heavy flavour production is in general reasonably well described. However the measured cross sections as a function of ∆φ show significant differences. These distributions are directly sensitive to additional parton radiation treated by parton showers or by other additional processes. In the following we investigate the influence of the details of the parton shower and the process gg * → gg with subsequent branching of one of the outgoing (off-shell) g → bb on the ∆φ distribution. For these studies we use the CCFM set A0.
In Cascade, the initial state parton shower is angular ordered and based on the CCFM evolution equations. The final state parton shower is also angular ordered, but based on the DGLAP evolution equations. To investigate the dependence on the parton shower, we calculated the cross sections without parton shower, only initial state, only final state and both initial and final state parton shower, as shown in Fig. 15. Here, b dijet production is selected as an example. We observe a very small contribution of the initial state parton shower, since in k T -factorisation the initial state parton shower does not influence the k T of the gluons (since it is determined from the unintegrated gluon density). However, a significant contribution comes from the final state parton shower. The prediction with full parton shower underestimates the back-to-back region of the azimuthal separation of the two b jets. This suggests that the final state parton showering generates too much gluon radiation, which causes less correlated jets as shown in Fig. 15 (a). To investigate this further we changed the final state parton shower scale Q 2 max from Q 2 max = 4m 2 T to Q 2 max = m 2 . As shown in Figure 15 (b), this leads to a higher correlation of the b jets and a very good description of the back-to-back region of the jets is achieved.
At lower values of the azimuthal separation of the dijet system, higher order processes are expected to contribute significantly. Therefore, we repeated the b dijet analysis for the process gg * → gg. In the matrix element calculation, which was performed in [35], one gluon in the initial state is on-shell and the other one is off-shell, as shown in Fig. 16. The heavy quark pair is then produced via parton showers in the final state as g → bb. As shown in Fig. 17 (a), this process contributes significantly to the tail of the ∆Φ jj distribution very well, without any additional adjustment of parameters. Part of the contribution from gg * → gg with g → bb is already included in the matrix element g * g * → bb, so simply adding both contributions would result in double counting [37]. In [22] the influence of multi-parton interactions in the collinear frame is investigated, and it was found, that the contribution even in the small ∆Φ jj region is small, because a large E T of the jets is required.
In conclusion, the azimuthal cross section measurements can be very well described by adjusting the scale for the final state parton shower to Q 2 max = m 2 and by including the process gg * → gg with subsequent branching g → bb, which in a collinear calculation would contribute only at NNLO.

Conclusions
We have studied charm and beauty production at Tevatron energies within the framework of the k T -factorization. Our calculations are based on the CCFM-evolved unintegrated gluon densities in a proton. The analysis covers the total and differential cross sections of open charm and beauty quarks, B and D mesons (or rather muons from their semileptonic decays). The cross sections of bb di-jet production have been studied also. Special attention was put on the specific angular correlations between the final-state particles. Using full hadron-level Monte Carlo generator Cascade, we investigated the effects coming from the parton showers in initial and final states. Different sources of theoretical uncertainties have been specially studied.
We obtain good agreement of our calculations for all observables and the recent experimental data taken by the D0 and CDF collaborations. We have demonstrated, that the parton shower plays a significant role in the description of the cross section. We obtain a good description of the correlations once the higher order process gg * → gg with subsequent g → bb splitting is included.

Acknowledgements
We thank S.P. Baranov for his encouraging interest and useful discussions. We are very grateful o S. Vallecorsa for many discussions on the CDF results. We are also grateful to S. Baranov [7,8,35]. The first column shows the LZ numerical results while the second one depicts the Cascade predictions. Here we use CCFM set A0 gluon density for illustration. The experimental data are from CDF [4,5].    The left histogram shows the LZ numerical results while the right plot depicts the Cascade predictions. The kinematical cuts applied are described in the text. Notation of all histograms is the same as in Fig. 3. The experimental data are from CDF [21].
[rad] jj Φ ∆ . In all cases the CCFM set A0 is used and the data are taken from CDF [21] and CDF [6]. g* g ME g* g ME + Figure 16: The process gg * → gg is shown schematically. The matrix elements are calculated in [36].  Figure 17: The processes g * g * → QQ (dashed line) and gg * → gg for b dijet production in (a) and D 0 , D * − meson production in (b). In all cases the CCFM set A0 is used and the data are taken from CDF [21] and CDF [6].