Production of J/psi and Upsilon mesons in pp collisions at sqrt(s) = 8 TeV

The production of J/psi and Upsilon mesons in pp collisions at sqrt(s) = 8 TeV is studied with the LHCb detector. The J/psi and Upsilon mesons are reconstructed in the mu+mu- decay mode and the signal yields are determined with a fit to the mu+mu- invariant mass distributions. The analysis is performed in the rapidity range 2.0<y<4.5 and transverse momentum range 0<p_T<14(15) GeV/c of the J/psi(Upsilon) mesons. The J/psi and Upsilon production cross-sections and the fraction of J/psi mesons from b-hadron decays are measured as a function of the meson p_T and y.


Introduction
Successfully describing heavy quarkonium production is a long-standing problem in QCD. An effective field theory, non-relativistic QCD (NRQCD) [1,2], provides the foundation for much of the current theoretical work. According to NRQCD, the production of heavy quarkonium factorises into two steps: a heavy quark-antiquark pair is first created at short distances and subsequently evolves non-perturbatively into quarkonium at long distances. The NRQCD calculations depend on the colour-singlet (CS) and colour-octet (CO) matrix elements, which account for the probability of a heavy quark-antiquark pair in a particular colour state to evolve into a heavy quarkonium state. The CS model (CSM) [3,4], which provides a leading-order description of quarkonium production, was initially used to describe experimental data. However, it underestimates the observed crosssection for single J/ψ production at high transverse momentum (p T ) at the Tevatron [5]. To resolve this discrepancy, the CO mechanism was introduced [6]. The corresponding matrix elements were determined from the high-p T data, as the CO cross-section decreases more slowly with p T than that predicted by CS. More recent higher-order calculations [7][8][9][10] close the gap between the CS predictions and the experimental data [11], reducing the need for large CO contributions.
In this paper first measurements of quarkonium production at √ s = 8 TeV are reported under the assumption of zero polarisation, an assumption that is discussed in the paper. The differential production cross-sections of prompt J/ψ and Υ mesons, produced at the pp collision point either directly or via feed-down from higher mass charmonium or bottomonium states, are presented in the range of rapidity 2.0 < y < 4.5 and p T < 14 GeV/c (J/ψ ) or p T < 15 GeV/c (Υ ). The fraction of J/ψ mesons from b-hadron decays, abbreviated as "J/ψ from b" in the following, is also measured in the same fiducial region.

The LHCb detector and data set
The LHCb detector [25] is a single-arm forward spectrometer covering the pseudorapidity range 2 < η < 5, designed for the study of particles containing b or c quarks. The detector includes a high precision tracking system consisting of a silicon-strip vertex detector surrounding the pp interaction region, a large-area silicon-strip detector located upstream of a dipole magnet with a bending power of about 4 Tm, and three stations of silicon-strip detectors and straw drift tubes placed downstream. The combined tracking system has a momentum resolution ∆p/p that varies from 0.4% at 5 GeV/c to 0.6% at 100 GeV/c, and an impact parameter resolution of 20 µm for tracks with high p T . Charged hadrons are identified using two ring-imaging Cherenkov detectors [26]. Photon, electron and hadron candidates are identified by a calorimeter system consisting of scintillating-pad and preshower detectors, an electromagnetic calorimeter and a hadronic calorimeter. Muons are identified by a system composed of alternating layers of iron and multiwire proportional chambers, with the exception of the centre of the first station, which uses triple-GEM detectors.
The data sample used in this analysis was collected during the first part of the data taking period at √ s = 8 TeV in April 2012. During this period the average number of interactions per crossing varied. The Υ meson analysis is based on a data sample, corresponding to an integrated luminosity of about 51 pb −1 of pp interactions, collected with an average of 1.3 visible interactions per crossing. The analysis for the more abundant J/ψ mesons is based on data, corresponding to an integrated luminosity of about 18 pb −1 , collected with an average of 1.0 visible interactions per crossing. The trigger [27] consists of a hardware stage, based on information from the calorimeter and muon systems, followed by a software stage, which applies a full event reconstruction. At the hardware stage, events are selected requiring dimuon candidates with a product of their p T larger than 1.68 ( GeV/c) 2 . In the subsequent software trigger, two well reconstructed tracks are required to have hits in the muon system, a p T higher than 500 MeV/c, p higher than 6 GeV/c and to form a common vertex. Only events with a dimuon candidate with an invariant mass m µµ within 120 MeV/c 2 of the known J/ψ meson mass [28] or larger than 4.7 GeV/c 2 are retained for further analysis.

Selection and cross-section determination
The selection is based on the criteria described in Refs. [12,13] and is summarised in Table 1. It starts by combining oppositely-charged particles, identified as muons, with a track p T larger than 700 (1000) MeV/c 2 for the J/ψ (Υ ) meson. Good track quality is ensured by requiring a χ 2 per degree of freedom, χ 2 /ndf, less than 4 for the track fit. Duplicate particles created by the reconstruction are suppressed to the level of 0.5 × 10 −3 using the Kullback-Leibler (KL) distance variable [29]. To ensure good quality vertex reconstruction, the χ 2 probability of the dimuon vertex is required to be larger than 0.5 %. In addition, the primary vertex (PV) associated to the dimuon candidate is required to be within the luminous region, defined as |x PV | < 1 mm , |y PV | < 1 mm and |z PV | < 150 mm.
In the J/ψ analysis additional criteria are applied to the vertex quality. The uncertainty on the pseudo decay time t z , defined in Eq. 2, is required to be less than 0.3 ps, as estimated by the propagation of the uncertainties given by the track reconstruction.
The simulation samples are based on the Pythia 6.4 generator [30] configured with the parameters detailed in Ref. [31]. The EvtGen package [32] is used to generate hadron decays. The interaction of the generated particles with the detector and its response are implemented using the Geant4 toolkit [33] as described in Ref. [34]. Radiative corrections to the decay of the vector meson to dimuons are generated with the Photos package [35].
The differential cross-section for the production of a vector meson V in a bin of (p T , y), where V stands for a J/ψ or Υ meson, decaying into a muon pair, is where N (V → µ + µ − ) is the number of observed V → µ + µ − candidates, tot the total detection efficiency in the given bin, L is the integrated luminosity, B (V → µ + µ − ) is the branching fraction of the V → µ + µ − decay and ∆y = 0.5 and ∆p T = 1 GeV/c are the rapidity and p T bin sizes, respectively. In the case of the J/ψ → µ + µ − decay the branching fraction is well known, B(J/ψ → µ + µ − ) = (5.94 ± 0.06) × 10 −2 [28], and therefore it is chosen to quote an absolute cross-section. On the other hand, the dimuon branching fractions of the Υ mesons are known less precisely [28], and therefore, as in Ref. [13], the product of the cross-section times the dimuon branching fraction is given. The total efficiency tot is the product of the geometric acceptance, the reconstruction and selection efficiency and the trigger efficiency. All efficiency terms are evaluated using simulated samples and validated with data-driven techniques in each (p T , y) bin.
The procedure to measure the integrated luminosity is described in Ref. [36]. For this analysis a van der Meer scan [37] was performed in April 2012, resulting in a measurement of the integrated luminosity of 18.4 ± 0.9 pb −1 for the J/ψ and 50.6 ± 2.5 pb −1 for the Υ samples.

J/ψ meson signal
As in the previous studies, prompt J/ψ mesons are distinguished from J/ψ from b by means of the pseudo decay time variable defined as  where z J/ψ and z PV are the positions along the beam axis z of the J/ψ decay vertex and of the primary vertex refitted after removing the decay muons of the J/ψ candidate; p z is the measured J/ψ momentum in the beam direction and M J/ψ is the known J/ψ mass [28].
The yields of both prompt J/ψ mesons and J/ψ from b are determined from a twodimensional fit in each (p T , y) bin to the distributions of invariant mass and pseudo decay time of the signal candidates, following the approach described in Ref. [12]. The mass distribution is modelled with a Crystal Ball function [38] for the signal and an exponential function for the combinatorial background.
The signal pseudo decay time distribution is described by a δ-function at t z = 0 for the prompt J/ψ component together with an exponential decay function for the J/ψ from b component. The shape of the tail arising from the association of a J/ψ meson candidate with a wrong primary vertex is derived from the data by combining a J/ψ meson from a given event with the primary vertex of the following event in the sample. The prompt component of the signal function and that from b hadron decays are convolved with a resolution function modelled by the sum of two Gaussian functions. The background distribution is parameterised with an empirical function based on the shape of the t z distribution observed in the J/ψ mass sidebands. It is built as the sum of a δ-function and five exponential components, three for positive t z and two for negative t z , the negative and positive exponential functions with the largest lifetime having their lifetimes τ L fixed to the same value. This function is convolved with a resolution function modelled by the sum of two Gaussian functions. All parameters of the background component are determined independently in each (p T , y) bin from the distribution of the pseudo decay time and are fixed in the final fit. The total fit function is the sum of the products of the mass and t z fit functions for the signal and background. Figure 1 shows the fit projections in mass and t z for one specific bin (3 < p T < 4 GeV/c, 2.5 < y < 3.0) with the fit result superimposed. Summing over all bins, a total signal yield of 2.6 million J/ψ events is obtained.

Υ meson signal
The Υ meson signal yields are determined from a fit to the reconstructed dimuon invariant mass of the selected candidates with 8.5 < m µµ < 11.5 GeV/c 2 . The distribution is described by the sum of three Crystal Ball functions, one for each of the Υ (1S), Υ (2S) and Υ (3S) signals, and an exponential function for the combinatorial background. The parameters α and n of the Crystal Ball function describing the radiative tail are fixed to the values of α = 2 and n = 1 based on simulation studies. The width of the Crystal Ball function describing the Υ (1S) meson is allowed to vary, while the widths of the Υ (2S) and Υ (3S) mesons are constrained to the value of the width of the Υ (1S) signal, scaled by the ratio of the masses of the Υ (2S) and Υ (3S) to the Υ (1S) meson. The peak values of the Υ (1S), Υ (2S) and Υ (3S) mass distributions are allowed to vary in the fit and are consistent with the known values [28]. Figure 2 shows the results of the fit performed over the entire range in p T and y. The obtained signal yields are 43 785 ± 254, 10 976 ± 155 and 5325 ± 122 for the Υ (1S), Υ (2S) and Υ (3S) mesons, respectively, with a mass resolution of the Υ (1S) resonance of 43 MeV/c 2 . The fit is repeated independently for each of the bins in p T and y. When fitting the individual bins, the masses are fixed to the values obtained when fitting the full range, while the mass resolution for the Υ (1S) candidates is parameterised with a linear function of p T and y that reproduces the behaviour observed in data. Bins with too few entries are excluded from the analysis.

Uncorrelated between bins
Production model 1.0 to 6.0 t z fit, for J/ψ from b 1.0 to 12.0

Systematic uncertainties
Previous studies [12,13] have shown that the total efficiency depends on the initial polarisation state of the vector meson. The J/ψ polarisation has been measured at √ s = 7 TeV by the LHCb [39] and ALICE [40] collaborations, in a kinematic range similar to that used in this analysis, and the Υ polarisation has been measured by CMS [41] at large p T and central rapidity. They were both found to be small. Therefore, in this paper results are quoted under the assumption of zero polarisation and no corresponding systematic uncertainty is assigned on the cross-section for this effect. All other systematic uncertainties are summarised in Table 2.
Uncertainties related to the mass model describing the shape of the dimuon mass distribution are estimated by fitting the invariant mass distributions for the J/ψ and Υ mesons with the sum of two Crystal Ball functions. The relative difference in the number of signal events (0.7-2.2%) is taken as a systematic uncertainty. A fraction of events has a lower invariant mass because of the energy lost through bremsstrahlung. Based on simulation studies, about 4% of the signal events are estimated to be outside the analysis mass windows and are not counted as signal. The fitted signal yields are corrected for this effect and an uncertainty of 1% is assigned to the cross-section measurement based on a comparison between the radiative tail observed in data and simulation.
The uncertainty due to the muon identification efficiency is measured on data using a tag-and-probe method. This method reconstructs J/ψ candidates in which one muon is identified by the muon system ("tag") and the other ("probe") is identified selecting a track depositing the energy of minimum-ionising particles in the calorimeters.
The ratio of the muon identification efficiency measured in data to that obtained in the simulation is convolved with the momentum distribution of muons from J/ψ and Υ mesons to obtain an efficiency correction. This is found to be 0.98 ± 0.01; the uncertainty on the correction factor is considered as a systematic uncertainty.
The uncertainty on the reconstruction efficiency of the muon tracks has also been estimated using a data-driven tag-and-probe approach based on partially reconstructed J/ψ decays, and it was found to be 0.9% per muon pair.
Differences between data and simulation in the efficiency of the requirement on the vector meson vertex χ 2 probability lead to a further uncertainty of 1%.
The trigger efficiency is determined using a data-driven method exploiting a sample of events that are still triggered when the signal candidate is removed [27]. The efficiency obtained with this method in each (p T , y) bin is used to check the efficiencies measured in the simulation. The systematic uncertainty associated with the trigger efficiency is the difference between that measured in the data and in the simulation. As a cross-check, the trigger efficiency is also computed using a data sample that has not been required to pass any physics trigger. The results obtained with the two methods are consistent.
The luminosity is determined with an uncertainty of 5%, dominated by differences in the results obtained with a van der Meer scan [37] using the core and off-core parts of the beam.
The dependence of the efficiency calculation on the production model used in the simulation is taken into account by varying the main parameters of the Pythia 6.4 generator related to prompt vector meson production. These parameters define the minimum p T cut-offs for regularising the cross-section. This effect is evaluated in each (p T , y) bin and found to be at most 6%.
Uncertainties related to the t z fitting procedure for the J/ψ mesons are included by changing the parameterisation used to describe the signal and background. A second fitting method based on the sPlot technique [42] is used with the mass as the control variable to unfold the background and to perform an unbinned likelihood fit to the pseudo decay time distribution. The two approaches give consistent results and their difference is taken as an estimate of the systematic uncertainty. These uncertainties are evaluated in each (p T , y) bin and found to be a few percent.

Results on J/ψ meson production
The measured double-differential production cross-sections for prompt J/ψ mesons, under the assumption of zero polarisation, and for J/ψ from b are given in bins of p T and y in Tables 4 and 5, respectively, and are displayed in Fig. 3. The integrated cross-section for prompt J/ψ meson production in the defined fiducial region, summing over all bins of the analysis, is where the first uncertainty is statistical and the second is systematic, computed taking correlations into account. The integrated cross-section for the production of J/ψ from b in the same fiducial region is σ (J/ψ from b, p T < 14 GeV/c, 2.0 < y < 4.5) = 1.28 ± 0.01 ± 0.11 µb.
The total bb production cross-section is computed as  Table 3: Differential production cross-section dσ/dy in nb for prompt J/ψ mesons (assumed unpolarised) and for J/ψ from b, integrated over p T . The first uncertainty is statistical, the second (third) is the part of the systematic uncertainty that is uncorrelated (correlated) between bins. y where the factor α 4π = 5.4 is an extrapolation factor of the cross-section from the measured to the full kinematic region. This factor is obtained using the simulation as described in Sect. 3. The inclusive b→J/ψ X branching fraction is B(b → J/ψ X) = (1.16 ± 0.10)% [28]. The resulting total bb cross-section is σ(pp → bbX) = 298 ± 2 ± 36 µb, where the first uncertainty is statistical and the second is systematic, which includes the uncertainty on B(b → J/ψ X). No systematic uncertainty has been included for the extrapolation factor α 4π estimated from the simulation. For comparison, the value of the extrapolation factor given by NLO calculations is 5.1 [43]. Table 3 and Fig. 4 show the differential production cross-section dσ/dy integrated over p T , for unpolarised prompt J/ψ mesons and J/ψ from b. For both components, the cross-section decreases significantly between the central and forward regions of the acceptance. Table 6 and Fig. 5 give the values of the fraction of J/ψ from b in the different bins of p T and y. The fraction of J/ψ mesons from b-hadron decays increases as a function of p T , and, at constant p T , decreases with increasing y, as seen in the study at √ s = 7 TeV [12].

Results on Υ meson production
The double-differential production cross-sections times the dimuon branching fractions for the Υ mesons in bins of p T and y are given in Tables 7, 8, and 9, with the assumption of no polarisation. The double-differential cross-sections are displayed in Fig. 6. The integrated cross-sections times dimuon branching fractions B iS = B(Υ (iS) → µµ), with i = 1, 2, 3, in the range p T < 15 GeV/c and 2.0 < y < 4.5 are measured to be where the first uncertainty is statistical and the second systematic. The cross-section times dimuon branching fractions for the three Υ states are compared in Fig. 7 as a function of p T and y. These results are used to evaluate the ratios R iS/1S of the Υ (2S) to Υ (1S) and Υ (3S) to Υ (1S) cross-sections times dimuon branching fractions. Most of the uncertainties cancel in the ratio, except those due to the size of the data sample, to the model dependence and to the choice of the fit function. The ratios R iS/1S as a function       of p T and y are given in Tables 10 and 11, respectively, and shown in Fig. 8, with the assumption of no polarisation. For this measurement the p T range has been restricted to p T < 14 GeV/c and the y range to 2.0 < y < 4.0 to ensure enough counts for the three Υ states in all bins. The ratios are constant as a function of y and increase as a function of p T , in agreement with previous observations by LHCb [13] and as reported by ATLAS [19] and CMS [21] at √ s = 7 TeV.

Comparison with theoretical models
The measured differential cross-sections for the production of prompt J/ψ mesons as a function of p T are compared in Fig. 9 to three theoretical models that assume no polarisation. The considered models are -an NRQCD model at next-to-leading order (NLO). The colour-octet matrix elements in this case are determined from a global fit to HERA, Tevatron and LHC data [44,45]; -an NNLO* CSM [9,10]; the notation NNLO* indicates that the calculation at next-to-next leading order is not complete and neglects part of the logarithmic terms; -an NLO CSM [7] with the input parameters related to the choice of scale and charm quark mass given in Ref. [44].
In these comparisons it should be noted that the predictions are for direct J/ψ meson production, whereas the experimental measurements include feed-down from higher charmonium states. In particular, the contribution from J/ψ mesons produced in radiative χ c decays in the considered fiducial range was measured to be at the level of 20% at √ s = 7 TeV [46]. Allowing for this contribution, as was seen in the previous studies, both the NNLO* CSM and the NLO NRQCD models provide reasonable descriptions of the experimental data. In contrast, the CSM at NLO underestimates the cross-section by an order of magnitude.
The results for the production of J/ψ from b can be compared to calculations based on the FONLL formalism [43,47]. This model predicts the b-quark production cross-   section, and includes the fragmentation of the b-quark into b-hadrons and their decay into J/ψ mesons. In Fig. 10 the data for the differential production cross-section as a function of p T and y at √ s = 8 TeV are compared to the FONLL predictions. Good agreement is observed. The prediction for the total cross-section in the fiducial range of this measurement is 1.34 +0.63 −0.49 µb, in good agreement with the result presented here. In Fig. 11 the measurements of the cross-section for J/ψ from b at √ s = 2.76 [14], 7 [12], and 8 TeV are compared to the FONLL prediction. The behaviour as a function of the centre-of-mass energy is in excellent agreement with the prediction. This gives confidence that this model can produce reliable predictions for the b-hadron cross-section at the higher energies expected at the LHC.
In Fig. 12 the cross-sections times dimuon branching fractions for the three Υ meson states are compared to the CSM NLO [7] and NNLO * theoretical predictions [9] as a function of p T . The NNLO* CSM provides a reasonable description of the experimental data, particularly for the Υ (3S) meson, which is expected to be less affected by feeddown. As for the prompt J/ψ meson production, the CSM at NLO underestimates the cross-section by an order of magnitude.

Conclusions
The differential production cross-sections for J/ψ and Υ mesons are measured as a function of p T and y in the forward region, 2.0 < y < 4.5. The analysis is based on a data sample, corresponding to an integrated luminosity of 18 pb −1 and 51 pb −1 for the J/ψ and Υ mesons, respectively, collected in the early part of 2012 at a centre-of-mass energy of √ s = 8 TeV. The production cross-sections of prompt J/ψ mesons and J/ψ from b are individually measured. An estimate of the bb total cross-section is also obtained.
The results are compared with several recent theoretical predictions in the LHCb acceptance. The NNLO* CSM and the NLO NRQCD model (for the J/ψ meson) provide a reasonable description of the experimental data on the production of prompt J/ψ and Υ mesons as a function of p T under the assumption of zero polarisation. This confirms the progress in the theoretical calculations of quarkonium hadroproduction, even if the uncertainties on the predictions are still large. Theoretical predictions based on FONLL calculations are found to describe well the measured cross-section for J/ψ from b and its dependence on centre-of-mass energy.  Table 4: Double-differential cross-section d 2 σ dp T dy in nb/( GeV/c) for prompt J/ψ meson production in bins of of p T and y, with the assumption of no polarisation. The first error is statistical, the second is the component of the systematic uncertainty that is uncorrelated between bins and the third is the correlated component.
p T ( GeV/c) 2.0 < y < 2.5 2.5 < y < 3.0 3.0 < y < 3.  Table 5: Double-differential cross-section d 2 σ dp T dy in nb/( GeV/c) for the production of J/ψ from b in bins of p T and y. The first error is statistical, the second is the component of the systematic uncertainty that is uncorrelated between bins and the third is the correlated component.   Table 7: Double-differential production cross-sections d 2 σ dp T dy × B 1S in pb/( GeV/c) for the Υ (1S) meson in bins of transverse momentum and rapidity, assuming no polarisation. The first error is statistical, the second is the component of the systematic uncertainty that is uncorrelated between bins and the third is the correlated component.
p T ( GeV/c) 2.0 < y < 2.5 2.5 < y < 3.0 3.0 < y < 3.5 Table 8: Double-differential production cross-sections d 2 σ dp T dy × B 2S in pb/( GeV/c) for the Υ (2S) meson in bins of transverse momentum and rapidity, assuming no polarisation. The first error is statistical, the second is the component of the systematic uncertainty that is uncorrelated between bins and the third is the correlated component. Regions where the number of events was not large enough to perform a measurement are indicated with a dash.
p T ( GeV/c) 2.0 < y < 2.5 2.5 < y < 3.0 3.0 < y < 3.  Table 9: Double-differential production cross-sections d 2 σ dp T dy × B 3S in pb/( GeV/c) for the Υ (3S) meson in bins of transverse momentum and rapidity, assuming no polarisation. The first error is statistical, the second is the component of the systematic uncertainty that is uncorrelated between bins and the third is the correlated component. Regions where the number of events was not large enough to perform a measurement are indicated with a dash. Table 10: Ratios of cross-sections Υ (2S) → µ + µ − and Υ (3S) → µ + µ − with respect to Υ (1S) → µ + µ − as a function of p T in the range 2.0 < y < 4.0, assuming no polarisation. The first error is statistical, the second is the component of the systematic uncertainty that is uncorrelated between bins and the third is the correlated component.  Table 11: Ratios of cross-sections Υ (2S) → µ + µ − and Υ (3S) → µ + µ − with respect to Υ (1S) → µ + µ − as a function of y in the range 2.0 < y < 4.0, assuming no polarisation. The first error is statistical, the second is the component of the systematic uncertainty that is uncorrelated between bins and the third is the correlated component.