Inclusive J/$\psi$ production at midrapidity in pp collisions at $\sqrt{s}~=~13$ TeV

We report on the inclusive J/$\psi$ production cross section measured at the CERN Large Hadron Collider in proton-proton collisions at a centre-of-mass energy $\sqrt{s}~=~13$ TeV. The J/$\psi$ mesons are reconstructed in the $\rm e^{+} e^{-}$ decay channel and the measurements are performed at midrapidity ($|y|<0.9$) in the transverse-momentum interval $0<p_{\rm T}<40$ GeV/$c$, using a minimum-bias data sample corresponding to an integrated luminosity $L_{\text{int}} = 32.2~\text{nb}^{-1}$ and an Electromagnetic Calorimeter triggered data sample with $L_{\text{int}} = 8.3~\mathrm{pb}^{-1}$. The $p_{\rm T}$-integrated J/$\psi$ production cross section at midrapidity, computed using the minimum-bias data sample, is $\text{d}\sigma/\text{d}y|_{y=0} = 8.97\pm0.24~(\text{stat})\pm0.48~(\text{syst})\pm0.15~(\text{lumi})~\mu\text{b}$. An approximate logarithmic dependence with the collision energy is suggested by these results and available world data, in agreement with model predictions. The integrated and $p_{\rm T}$-differential measurements are compared with measurements in pp collisions at lower energies and with several recent phenomenological calculations based on the non-relativistic QCD and Color Evaporation models.


Introduction
Quarkonium production in hadronic interactions is an excellent case of study for understanding hadronization in quantum chromodynamics (QCD), the theory of strong interactions [1]. In particular, the production of the J/ψ meson, a bound state of a charm and an anti-charm quark and the lightest vector charmonium state, is the subject of many theoretical calculations. The cornerstone of all the theoretical approaches is the factorization theorem, according to which the J/ψ production cross section can be factorized into a short distance part describing the cc production and a long distance part describing the subsequent formation of the bound state. In this way, the cc pair production cross section can be computed perturbatively. The widely used Non-Relativistic QCD (NRQCD) approach [2] describes the transition probabilities of the pre-resonant cc pairs to bound states with a set of long-distance matrix elements (LDME) fitted to experimental data, assumed to be universal. Next-to-leading order (NLO) calculations involving collinear parton densities are able to describe production yields for transverse momentum (p T ) larger than the mass of the bound state [3,4], but have difficulties describing the measured polarization [5,6]. Calculations employing the k T -factorization approach [7] can reach lower p T but have similar difficulties when compared to data [8]. The low-p T range of quarkonium production is modelled also within the Color Glass Condensate effective theory coupled to leading order NRQCD calculations [9], which involves a saturation of the small Bjorken-x gluon densities that dampens the heavy-quark pair production yields. An alternative to the universal LDME approach to hadronization used in the NRQCD framework is provided by the Color Evaporation Model (CEM) [10,11] and its more recent implementation using the k T -factorization approach, the Improved CEM (ICEM) [12]. In the ICEM, the transition probability to a given bound state is proportional to the cc pair production cross section integrated over an invariant-mass range spanning between the mass of the bound state and twice the mass of the lightest charmed meson. Finally, in the Color Singlet Model (CSM) [13][14][15], the preresonant cc pair is produced directly in the color-singlet state with the same quantum numbers as the bound state. Calculations within this model at NLO precision are known to strongly underpredict the measured production cross sections [16]. In this context, a p T -differential measurement of J/ψ production cross section covering a wide p T range, starting from p T =0 and up to high-p T , can discriminate between the different models of quarkonium production.
In this paper, we present the integrated, and the p T and rapidity (y) differential production cross sections of inclusive J/ψ production at midrapidity (|y| < 0.9) in proton-proton (pp) collisions at the center-ofmass energy √ s = 13 TeV. The inclusive J/ψ yields include contributions from directly produced J/ψ, feed-down from prompt decays of higher-mass charmonium states, and non-prompt J/ψ from the decays of beauty hadrons. The p T -differential production cross section of inclusive J/ψ is measured in the 0 < p T < 15 GeV/c interval using a minimum-bias triggered data sample and in the 15 < p T < 40 GeV/c interval using an Electromagnetic Calorimeter triggered data sample. These results complement existing measurements at midrapidity at √ s = 13 TeV performed by the CMS Collaboration [17], which report the prompt J/ψ production cross section for p T > 20 GeV/c. Previous measurements of the J/ψ production cross section in pp collisions performed by the ALICE Collaboration at midrapidity at lower energies were published in Refs. [18][19][20]. The inclusive J/ψ production cross section in pp collisions at √ s = 13 TeV was published by the ALICE Collaboration at forward rapidity in Ref.
[21] and the prompt and non-prompt production cross sections were reported by the LHCb Collaboration in Ref. [22].
In the next sections, the ALICE detector and the data sample are described in Section 2, and the data analysis and the determination of the systematic uncertainties are described in Section 3 and Section 4, respectively. The results are presented and discussed together with recent model calculations in Section 5 and conclusions are drawn in Section 6.
Inclusive J/ψ production at midrapidity in pp collisions at √ s = 13 TeV ALICE collaboration 2 The ALICE detector, data set and event selection A detailed description of the ALICE detector and its performance is provided in Refs. [23,24]. Here we mention only the detector systems used for the reconstruction of the J/ψ mesons decaying in the e + e − channel at midrapidity. Unless otherwise specified, the term electrons will be used throughout the text to refer to both electrons and positrons.
The reconstruction of charged-particle tracks is performed using the Inner Tracking System (ITS) [25] and the Time Projection Chamber (TPC) [26], which are placed inside a solenoidal magnet providing a uniform magnetic field of B = 0.5 T oriented along the beam direction. The ITS is a silicon detector consisting of six cylindrical layers surrounding the beam pipe at radii between 3.9 and 43.0 cm. The two innermost layers consist of silicon pixel detectors (SPD), followed by two layers of silicon drift (SDD) and two layers of silicon strip (SSD) detectors. The TPC is a cylindrical gas drift chamber which extends radially between 85 and 250 cm and longitudinally over 250 cm on each side of the nominal interaction point. Both TPC and ITS have full coverage in azimuth and provide tracking in the pseudorapidity range |η| < 0.9. Additionally, the measurement of the specific energy loss (dE/dx) in the TPC active gas volume is used for electron identification. The Electromagnetic Calorimeter (EMCal) and the Di-jet Calorimeter (DCal) [27][28][29] are employed for triggering and electron identification. The EMCal/DCal is a Shashlik-type lead-scintillator sampling calorimeter located at a radius of 4.5 m from the beam vacuum tube. The EMCal detector covers a pseudorapidity range of |η| < 0.7 over an azimuthal angle of 80 • < ϕ < 187 • , and the DCal covers 0.22 < |η| < 0.7 for 260 • < ϕ < 320 • and |η| < 0.7 for 320 • < ϕ < 327 • . The EMCal and DCal have identical granularity and intrinsic energy resolution, and they form a two-arm electromagnetic calorimeter, which in this paper will be referred to jointly as EMCal.
In addition to these central barrel detectors, the V0 detectors, composed of two scintillator arrays [30] placed along the beam line on either side of the interaction point and covering the pseudorapidity intervals −3.7 < η < −1.7 and 2.8 < η < 5.1, respectively, are used for event triggering. Together with the SPD detector, the V0 is also used to reject background from beam-gas collisions and pileup events.
The measurement of the p T -integrated and p T -differential production cross sections up to p T = 15 GeV/c, the upper limit being determined by the available integrated luminosity, utilizes the minimum-bias (MB) trigger, defined as the coincidence of signals in both V0 scintillator arrays. For the p T interval from 15 up to 40 GeV/c, the EMCal trigger is employed to select events with high-p T electrons. The lower p T limit is chosen such that the trigger efficiency does not vary with p T above this value, thus avoiding systematic uncertainties related to trigger threshold effects. The EMCal trigger is an online trigger which includes a Level 0 (L0) and a Level 1 (L1) component [27]. The calorimeter is segmented into towers and Trigger Region Units (TRUs), the latter being composed of 384 towers each [27]. The L0 trigger is based on the analog charge sum of 4 × 4 groups of adjacent towers evaluated within each TRU, in coincidence with the MB trigger. The L1 trigger decision requires the L0 trigger and, in addition, scans for 4 × 4 groups of adjacent towers across the entire EMCal surface. The EMCal-triggered analysis presented in this paper uses the L1 trigger, and requires that at least one of the charge sums of the 4 × 4 adjacent towers is above 9 GeV.
This analysis includes all the data recorded by the ALICE Collaboration during the LHC Run 2 datataking campaigns of 2016, 2017 and 2018 for pp collisions at √ s = 13 TeV. The maximum interaction rate for the dataset was 260 kHz, with a maximum pileup probability in the same bunch crossing of 0.5×10 −3 . The events selected for analysis were required to have a reconstructed vertex within the interval |z vtx | < 10 cm to ensure a uniform detector acceptance. Beam-gas events and pileup collisions occurring within the readout time of the SPD were rejected offline using timing selections based on the V0 detector information. Pileup collisions occurring within the same LHC bunch crossing were rejected using offline algorithms which identify multiple vertices [24]. The remaining fraction of pileup events Inclusive J/ψ production at midrapidity in pp collisions at √ s = 13 TeV ALICE collaboration surviving the selections is negligible for both the MB and EMCal data samples.
The analyzed MB sample, satisfying all the quality selections, consists of about 2×10 9 events, corresponding to an integrated luminosity L int = 32.2 nb −1 ± 1.6% (syst), and the EMCal-triggered sample consists of approximately 9×10 7 events which corresponds to an integrated luminosity L int = 8.3 pb −1 ± 2.0% (syst). The integrated luminosities are obtained based on the MB trigger cross section (σ MB ), measured in a van der Meer scan [31], separately for each year, as described in Ref. [32]. For each of the used triggers, MB and EMCal, the integrated luminosity is obtained as where N MB is the number of MB-triggered events in the triggered sample, ds trig is the downscaling factor applied to the considered trigger by the ALICE trigger processor and LT trig is the trigger live time, i.e. the fraction of time where the detector cluster 1 assigned to the trigger was available for readout.

J/ψ reconstruction
In this work we study the integrated, and the rapidity-and p T -differential inclusive J/ψ production at midrapidity (|y| < 0.9) reconstructing the J/ψ from the e + e − decay channel. The MB sample analysis follows closely the one performed in pp collisions at √ s = 5.02 TeV [20].

Track selection
Electron-track candidates are reconstructed employing the ITS and TPC. They are required to be within the acceptance of the central barrel (|η| < 0.9), and to have a minimum transverse momentum of 1 GeV/c, which suppresses the background with only a moderate J/ψ efficiency loss. The tracks are selected to have at least 2 hits in the ITS, one of which having to be in one of the SPD layers, and share at most one hit with other tracks. A minimum of 70 out of a maximum of 159 clusters are required in the TPC.
In order to reject tracks originating from weak decays and interactions with the detector material, a selection based on the distance-of-closest approach (DCA) to the primary vertex is applied to the tracks. For the MB analysis the tracks are required to have a minimum DCA lower than 0.2 cm in the transverse direction and 0.4 cm along the beam axis. Such tight selection criterion is used in order to improve the signal-to-background ratio and the signal significance. It was checked with Monte Carlo (MC) simulations that these requirements do not lead to efficiency loss for the non-prompt J/ψ relative to the prompt J/ψ. For the EMCal-triggered event analysis, a looser selection on the DCA to the primary vertex at 1 and 3 cm is applied to avoid rejecting non-prompt J/ψ from highly boosted beauty hadron decays.
The electrons are identified using the specific energy loss dE/dx in the TPC gas. Their dE/dx is required to be within a band of [−2, 3] σ relative to the expectation for electrons at the given track momentum, with σ being the dE/dx resolution. The contamination from protons which occurs in the momentum range p < 1.5 GeV/c and from pions for momenta above 2 GeV/c is mitigated by rejecting tracks compatible with the proton or pion hypothesis within 3σ .
For the analysis of EMCal-triggered events, both the TPC and the EMCal are used for electron identification. At least one of the J/ψ decay-electron tracks, initially identified by the TPC, is required to be matched to an EMCal cluster (a group of adjacent towers belonging to the same electromagnetic shower).
In order to ensure a constant trigger efficiency on the selected events, the matched clusters are selected to have a minimal energy of 14 GeV, a value that is significantly higher than the applied online threshold of 9 GeV. Electrons are identified by applying a selection on the energy-to-momentum ratio of the EMCal Secondary electrons from photon conversions, the main background source for both analyses, are rejected using the requirement of a hit in the SPD detector. This requirement rejects most of the electrons from photon conversions occurring beyond the SPD layers. An additional selection based on track pairing, as described in details in Ref.
[20], is applied to further reject conversion electrons, especially those from photons converting in the beam pipe or in the SPD.

Signal extraction
The number of reconstructed J/ψ mesons is extracted from the invariant mass (m ee ) distribution of all possible opposite-sign (OS) pairs constructed combining the selected electron tracks within the same event (SE). Besides the J/ψ signal, i.e. pairs of electrons originating from the decay of a common J/ψ mother, the invariant mass distribution contains a background with contributions from combinatorial and correlated sources.
In the MB analysis, the combinatorial background, i.e. pairs of electrons originating from uncorrelated processes, is estimated using the event-mixing technique (ME), in which pairs are built from oppositesign electrons belonging to different events. The mixing is done considering events from the same run (a collection of events taken during a period of time of up to a few hours) with a similar vertex position. The normalized combinatorial background distribution B comb (m ee ) is obtained as where N SE LS , N ME OS and N ME LS are the number of same-event like-sign (LS), mixed-event OS and mixed-event LS pairs, respectively. Here, the mixed-event OS distribution is normalized using the ratio of SE to ME like-sign pairs since these are not expected to contain any significant correlated source. The summation extends over all the mass bins m i between 0 and 5 GeV/c 2 to minimize the statistical uncertainty on the background matching. The correlated background in the mass region relevant for this analysis originates mainly from semi-leptonic decays of heavy-flavor hadrons [33]. In order to extract the number of reconstructed J/ψ, N J/ψ , the combinatorial background-subtracted invariant mass distribution is fitted with a two-component function: one empirical function to describe the correlated background shape, which is a second order polynomial at low pair p T (p ee T < 1 GeV/c) and an exponential at high-p T (p ee T > 1 GeV/c), plus a template shape obtained from MC simulations, described in Section 3.3, for the J/ψ signal.
For the analysis of the EMCal-triggered event sample, due to the relatively large contribution from correlated sources at high p T , the event mixing technique is not used. Instead, a fit of the invariant mass distribution is performed using the MC template for the signal and a third-order polynomial function to describe both the combinatorial and correlated background contributions.
In both analyses, the contribution from ψ(2S) decaying in the dielectron channel is not included in the fit as the expected number of such pairs is ∼1% of the J/ψ raw yield and it is statistically not significant in the analyzed data samples. The number of J/ψ is obtained by counting the number of e + e − pairs in the mass range 2.92 ≤ m ee ≤ 3.16 GeV/c 2 remaining after subtracting the background. The SE-OS dielectron invariant mass distribution for a few of the p T intervals is shown in Fig. 1, together with the estimated signal and background components.

Corrections
The double differential J/ψ production cross section is calculated as where N J/ψ is the number of reconstructed J/ψ in a given interval of rapidity ∆y and transverse momentum ∆p T , BR(J/ψ → e + e − ) is the decay branching ratio into the dielectron channel [34], A × ε is the average acceptance and efficiency factor and L int is the integrated luminosity of the data sample.
The correction for acceptance and efficiency is the product of the kinematical acceptance factor, the reconstruction efficiency, which includes both tracking and particle-identification (PID) efficiency, and the fraction of signal in the signal counting mass window. For the EMCal-triggered events analysis, the efficiency for the EMCal cluster reconstruction is also considered. The efficiency related to the EMCal trigger is estimated using a parameterized simulation of the L1 trigger which, includes decalibration and noise based on measured data, and takes into account the time-dependent detector conditions. With the exception of the PID efficiency, all the corrections are obtained based on a MC simulation of unpolarized J/ψ mesons embedded in inelastic pp collisions simulated using PYTHIA 6.4 [35] with the Perugia 2011 tune [36]. The prompt J/ψ are generated with a flat rapidity distribution and a p T spectrum obtained from a phenomenological interpolation of J/ψ measurements at RHIC, CDF and the LHC at lower energies [37]. For the non-prompt J/ψ, bb pairs are generated using the PYTHIA Perugia 2011 tune. The J/ψ decays are simulated using PHOTOS [38], which includes the radiative component of the J/ψ decay. The generated particles are transported through the ALICE detector setup using the GEANT3 package [39].
The PID efficiency is determined with a data-driven method by using a clean sample of electrons from tagged photon conversion processes, passing the same quality criteria as the electrons selected for the J/ψ reconstruction. The PID selection efficiency for single electrons is propagated to the J/ψ level using a simulation of the J/ψ decay. The acceptance times efficiency correction factor for the MB Inclusive J/ψ production at midrapidity in pp collisions at √ s = 13 TeV ALICE collaboration sample analysis varies with p T between 7.6% and 16% while in the case of the EMCal-triggered sample analysis it increases with p T from 2 to 8%.
Due to the finite size of the p T intervals, there is a mild dependence of the correction factors on the shape of the inclusive J/ψ p T distribution used in the simulation. This is mitigated iteratively by using the corrected J/ψ p T -differential production cross section to reweight the acceptance times efficiency correction factor and obtain an updated corrected cross section. The procedure is stopped when the difference between the input and output corrected p T -differential production cross section drops below 1%, which typically occurs within 1 to 2 iterations, depending on the p T interval. Additionally, to check if the default MC used in the analysis could introduce a bias on the EMCal trigger efficiency, due to the enhancement of J/ψ, another MC simulation, based on a di-jet production generated by PYTHIA8 [40], at √ s = 13 TeV, is used as a cross-check. As a result, the default MC and the di-jet MC lead to a compatible EMCal trigger efficiency.

Systematic uncertainties
There are several sources of systematic uncertainties affecting this analysis, namely the ITS-TPC tracking, the electron identification, the signal extraction procedure, the J/ψ input kinematic distributions used in MC simulations, the determination of the integrated luminosity and the branching ratio of the dielectron decay channel. A summary of these is given in Table 1.
The uncertainty of the ITS-TPC tracking efficiency is one of the dominant sources of systematic uncertainty and has two contributions: one due to the TPC-ITS matching and one related to the track-quality requirements. The former is obtained from the residual difference observed for the ITS-TPC single-track matching between data and MC simulations [41], which is further propagated to J/ψ dielectron pairs. It varies between 2.8% and 5.4%, depending on p T . The uncertainty related to the track-quality requirements is estimated by repeating the analysis with variations of the selection criteria and taking the root mean square (RMS) of the distribution of the results as systematic uncertainty. This uncertainty also depends on p T and is equal to 3.7% for p T < 5 GeV/c and approximately 2% for p T > 5 GeV/c. In Table 1, both contributions are added in quadrature and provided as ranges for the p T -and y-differential results.
As described in Section 3.3, the particle identification efficiency is determined via a data-driven procedure using a sample of identified electrons from tagged photon conversions. For the MB data sample, the uncertainty of this procedure is estimated by repeating the analysis with a looser and a tighter hadron (pion and proton combined) rejection criteria and taking the largest deviation from the results obtained with the standard PID selection divided by √ 12. In addition, the statistical uncertainty of the pure electron sample used for the determination of the efficiency, which becomes non-negligible at high p T , is propagated to the total uncertainty for the PID. The total PID systematic uncertainty is larger than 1% only for p T > 7 GeV/c, reaching 4% in the p T interval 10 < p T < 15 GeV/c. For the EMCal-triggered analysis, the particle identification systematic uncertainty has contributions from the electron identification in both the TPC and the EMCal. The values are estimated by varying the dE/dx range for the electron selection in the TPC, the E/p selection range in the EMCal, and the minimum energy of the matched clusters. The total PID uncertainty, obtained in an analogous way as for the MB sample, increases with p T from 2.8% to 5.4%.
For p T up to 15 GeV/c, the systematic uncertainties associated to the J/ψ signal extraction procedure is dominated by the J/ψ invariant mass signal shape template used to calculate the fraction of reconstructed J/ψ mesons within the signal counting mass window. This uncertainty amounts to 1.9% and is evaluated by repeating the extraction of the corrected yield with different invariant mass intervals used for the signal counting and taking the RMS of the variations as a systematic uncertainty. An additional source of uncertainty, due to the fit of the correlated background, is determined by varying the fit mass range and is typically below 1%. For p T >15 GeV/c, the systematic uncertainty associated to the J/ψ signal Inclusive J/ψ production at midrapidity in pp collisions at √ s = 13 TeV ALICE collaboration Table 1: Summary of contributions to systematic uncertainties of the measured J/ψ production cross section (in percentage).

Source
(p T , y)-integrated y-differential p T -differential (MB) p T -differential (EMCal) 0 < p T < 15 (GeV/c) 15  extraction are evaluated similarly to the MB analysis, however, the components associated to the signal shape and the background fit have similarly large contributions. This is one of the main sources of uncertainty in this p T range. The total uncertainty of the signal extraction varies with p T between about 2 and 7%.
The uncertainty on the EMCal trigger efficiency is studied varying the contribution from random noise (evaluated using energy resolution measurements in data) applied to the 4 × 4 groups of adjacent towers. It was also studied using the di-jet production generated by PYTHIA8 at √ s = 13 TeV mentioned in Sec. 3.3. The latter study showed a difference in the EMCal trigger efficiency of 0.5%, which is assigned as a systematic uncertainty.
Since the J/ψ efficiency has a dependence on p T , the particular J/ψ kinematic distribution used in the MC simulation, from which efficiencies are derived, can have an impact on the average efficiency computed for finite p T intervals. As described in Section 3.3, this is mitigated via an iterative procedure and the remaining uncertainty is related to the precision of the measured J/ψ spectrum. This uncertainty is estimated by fitting the measured p T spectrum with a power law function and allowing the fitted parameters to vary according to the covariance matrix. For each of such a variation, the average efficiency in a given kinematic interval is recomputed and the RMS of the distribution obtained from the variation of efficiency with respect to the central value is taken as a systematic uncertainty. This amounts to 1% for the p T -integrated case and is smaller for all of the considered p T intervals.
The uncertainty of the integrated luminosity arising from the vdM-scan based measurements amounts to 1.6% for both the MB and EMCal-triggered data samples, and is determined as described in Ref. [32]. For the EMCal-triggered data sample only, an additional uncertainty of 1.1% arising from the precision of the trigger downscaling is assigned. The uncertainty on the J/ψ decay branching ratio to the dielectron channel amounts to 0.53% as reported by the Particle Data Group [34].
The uncertainties of the integrated luminosity and branching ratio are treated as global systematic uncertainties. All the other uncertainties are considered as point to point correlated, with the exception of the one due to the background fit which is considered to be fully uncorrelated. The systematic uncertainties are thus dominated by correlated sources. The total correlated, uncorrelated and global uncertainties obtained as the sum in quadrature of the corresponding sources are given in Table 1.
Inclusive J/ψ production at midrapidity in pp collisions at √ s = 13 TeV ALICE collaboration

Results
The inclusive p T -differential J/ψ production cross section in pp collisions at √ s = 13 TeV at midrapidity, obtained from the combined analysis of the MB-triggered (p T < 15 GeV/c) and EMCal-triggered (15 < p T < 40 GeV/c) samples is shown in the upper panels of Fig. 2. Statistical uncertainties are shown as error bars, while the boxes around the data points represent the correlated and uncorrelated systematic uncertainties added in quadrature, excluding the global uncertainty from luminosity determination and branching ratio. The x-axis extent of the boxes illustrates the size of the p T interval with the data points placed in the center. A simple power law function of the type with A, p 0 and n being free parameters, is used to fit the measured distribution. Since the systematic uncertainties are largely correlated, only the statistical ones are used in the fit. Due to the large p T intervals, the fit is performed by considering the mean value of the function in the p T interval rather than its value at the center of the interval. The values obtained for the fitted parameters are A = 2.15 ± 0.18 µb/(GeV/c) 2 , p 0 = 4.09 ± 0.22 GeV/c and n = 3.04 ± 0.09. The fit function, shown in Fig. 2 as a dashed red line, provides a good description of the data points measured with both MB and EMCal data samples, and illustrates the consistency between the two analyses.
In the left panel of Fig. 2, the production cross section measured at midrapidity is compared with the forward rapidity measurement at the same energy [21]. The bottom panel shows the ratio between the forward rapidity data points and the mean cross section at midrapidity obtained by integrating the fit function described above in the p T intervals of the forward-rapidity measurement. The displayed uncertainty boxes include the systematic uncertainty of the forward rapidity measurement and the uncertainty of the function mean, added in quadrature. A monotonic drop of this ratio can be observed towards high p T , indicating a harder J/ψ p T distribution at midrapidity.
The right panel of Fig. 2 shows a comparison with midrapidity measurements performed by the ALICE Collaboration in pp collisions at the lower collision energies of √ s = 7 TeV [18] and 5.02 TeV [20]. In the bottom panel, the ratio between the lower energy measurements and the fitted 13 TeV results is shown. The displayed uncertainty boxes include the systematic uncertainty of the lower energy data points and the uncertainty of the fit function mean, added in quadrature. Although the uncertainties of the measurement at 7 TeV are large, the data indicates an increase of the production cross section with increasing collision energy. In addition, the monotonic drop of the ratio between the 5.02 TeV and 13 TeV measurements indicates a hardening of the p T spectrum with increasing collision energy, as also observed from the energy dependence of the inclusive J/ψ average p T discussed in Ref. [20].
The measured inclusive J/ψ production cross section is compared with several phenomenological calculations of the prompt J/ψ production in the left panel of Fig. 3. In addition, to illustrate the impact of the unaccounted feed-down from beauty decays in the theory predictions, a calculation of the nonprompt J/ψ production by Cacciari et al., using the Fixed-Order Next-to-Leading-Logarithms approach (FONLL) [45], is shown in the same panel. According to the FONLL calculations, the non-prompt contribution to the inclusive J/ψ yield is approximately 10% at low-p T and grows to approximately 50% at p T =40 GeV/c. In the right panel of Fig. 3, the measured inclusive production cross section is compared to predictions for inclusive J/ψ production obtained as the sum of the prompt J/ψ calculations listed above and the beauty feed-down contribution calculated using FONLL. The bottom panel shows the ratio between each theoretical calculation and the fit to the data. The colored bands represent the theoretical uncertainties for each model, centered around the model to data ratio. These uncertainties are typically due to the variation of the renormalization and factorization scales, and of the charm quark mass.
Inclusive J/ψ production at midrapidity in pp collisions at √ s = 13 TeV ALICE collaboration 3.4% ±  The error bars represent statistical uncertainties while the boxes around the data points represent the total systematic uncertainty, excluding the global uncertainty from the luminosity and branching ratio. The lower panels show the ratio between the measurements at different rapidity and energies, and the power law fit, discussed in the text, to the J/ψ production cross section (red dashed line) at midrapidity at √ s = 13 TeV. The boxes in the lower panels include the systematic uncertainty of the data points and the uncertainty of the integral of the fit function in the given p T interval added in quadrature.
Inclusive J/ψ production at midrapidity in pp collisions at √ s = 13 TeV ALICE collaboration  : (Color online) Inclusive J/ψ production cross section compared with calculations for the prompt J/ψ production cross section using ICEM [42], NLO NRQCD [3,4,43], LO NRQCD+CGC [44] and for the nonprompt J/ψ from beauty-hadron feed-down using FONLL [45] (left panel). Inclusive J/ψ production cross section compared with the corresponding calculations obtained as the sum of the prompt J/ψ component shown in the left panel and the non-prompt contribution from FONLL (right panel). The bottom panel shows the ratios between the model calculations and a fit to the data points. The bands illustrate the theoretical uncertainties centered around the ratio between the model calculation and the power-law fit to the data (see text for details).
Inclusive J/ψ production at midrapidity in pp collisions at √ s = 13 TeV ALICE collaboration  Several phenomenological approaches are used for the calculation of the J/ψ yields shown in Fig. 3. The green and blue dashed lines represent NLO NRQCD calculations from Ma et al. [4] and Butenschoen et al. [3], respectively, using collinear gluon parton distribution functions (PDFs). Although the calculation of the short distance terms is very similar, the predictions of these two approaches differ due to the LDME sets which are obtained in separate fits of the Tevatron and HERA data with different lowp T cutoffs. In addition, the calculation from Ref. [3] does not include the feed-down from higher mass charmonium states. The brown solid line represents a calculation obtained with the MC generator PEGA-SUS [43] developed by Lipatov et al. which employs a k T -factorization approach using p T -dependent gluon distribution functions and NRQCD matrix elements combined with LDMEs extracted from an NLO high-transverse momentum analysis [8]. Using the KMR [46] technique to construct the unintegrated gluon PDFs, this calculation can extend down to J/ψ p T = 0. A different model to calculate the low p T J/ψ production cross section, by Ma and Venugopalan [44] (green solid line) is based on a Color-Glass Condensate (CGC) approach coupled to a Leading Order (LO) NRQCD calculation which includes a soft-gluon resummation. The calculations obtained using the ICEM model by Cheung and Vogt [42] within the k T -factorization approach are shown by the violet solid line. In this calculation, LDMEs are not used, however one normalization parameter per charmonium state is used to account for long distance effects [47]. The feed down from the higher mass charmonium states is taken into account in this model.
As shown in the right panel of Fig. 3, all the models provide a reasonable description of the inclusive J/ψ production cross section within the theoretical uncertainties over the entire p T range covered by this measurement. In particular, both the ICEM and NRQCD+CGC calculations show very good agreement with the data in the low-p T range. The NRQCD calculation from Lipatov, which uses the k T -factorization approach, also provides a good description of the data for p T > 2 GeV/c, while it overestimates the measured cross section at lower p T . However, this is a significant progress compared to traditional collinear approaches which tend to diverge towards p T = 0.
The systematic uncertainty includes all the systematic sources mentioned in Table 1 added in quadrature, with the exception of the global ones which are given separately. The fraction of the total cross section covered by the EMCal-triggered event analysis results (p T > 15 GeV/c) is estimated to be less than 0.5% and is not used in this measurement. A comparison of this measurement with previous midrapidity measurements in pp collisions at lower energies from PHENIX [48], The right panel of Fig. 4 shows the rapidity-differential inclusive J/ψ production cross section which includes three data points measured in this analysis around midrapidity and the previously published ALICE measurement at forward rapidity [21]. While no rapidity dependence is observed in the central rapidity range, a steep decrease towards forward rapidity is seen. The two model calculations employed, using ICEM [42] and the NRQCD+CGC [44], combined with non-prompt contributions calculated with FONLL [45], exhibit rather different rapidity dependences. However, both are compatible with the data owing to the large theoretical uncertainties, which are at present much larger than the experimental uncertainties.

Conclusions
The integrated, and the p T -and y-differential inclusive J/ψ production cross sections at midrapidity (|y| < 0.9) in pp collisions at √ s = 13 TeV are presented in the p T range 0 < p T < 40 GeV/c, exceeding the p T range of all the previous measurements reported by the ALICE Collaboration. The measurement up to 15 GeV/c is performed using a minimum-bias triggered data sample and the one at high-p T is performed using an EMCal-triggered data sample, both collected by ALICE during the LHC Run 2.
The data are compared with the ALICE lower collision energy measurements at midrapidity [18,20] and with the ALICE forward-rapidity results [21]. An approximate logarithmic dependence of the integrated J/ψ production cross section with collision energy is suggested by the data, in agreement with the available predictions. The p T -differential cross section measured in this analysis shows a significant hardening with respect to both the forward-rapidity measurement at √ s = 13 TeV and the midrapidity measurement at √ s = 5.02 TeV.
Several calculations within the NRQCD framework [3,4,43,44] and one using the ICEM approach [42] are compared with the measured inclusive J/ψ production cross section. In particular, the ICEM, the NRQCD model based on the CGC approach and a new NRQCD calculation using the k T -factorization approach provide a prediction for the J/ψ production cross section down to p T = 0. All models provide a good description of the measured inclusive J/ψ production cross section, although with large theoretical uncertainties.
Inclusive J/ψ production at midrapidity in pp collisions at √ s = 13 TeV ALICE collaboration The ALICE Collaboration would like to thank all its engineers and technicians for their invaluable contributions to the construction of the experiment and the CERN accelerator teams for the outstanding performance of the LHC complex.  Inclusive J/ψ production at midrapidity in pp collisions at √ s = 13 TeV ALICE collaboration A The ALICE collaboration