Inclusive J/ψ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {J}/\psi $$\end{document} production at midrapidity in pp collisions at s=13\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{s} = 13$$\end{document} TeV

We report on the inclusive J/ψ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {J}/\psi $$\end{document} production cross section measured at the CERN Large Hadron Collider in proton–proton collisions at a center-of-mass energy s=13\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{s}~=~13$$\end{document} TeV. The J/ψ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {J}/\psi $$\end{document} mesons are reconstructed in the e+e-\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {e}^{+}\text {e}^{-}$$\end{document} decay channel and the measurements are performed at midrapidity (|y|<0.9\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$|y|<0.9$$\end{document}) in the transverse-momentum interval 0<pT<40\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0<p_{\mathrm{T}} <40$$\end{document} GeV/c\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c$$\end{document}, using a minimum-bias data sample corresponding to an integrated luminosity Lint=32.2nb-1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L_{\text {int}} = 32.2~\text {nb}^{-1}$$\end{document} and an Electromagnetic Calorimeter triggered data sample with Lint=8.3pb-1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L_{\text {int}} = 8.3~\mathrm {pb}^{-1}$$\end{document}. The pT\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_{\mathrm{T}}$$\end{document}-integrated J/ψ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {J}/\psi $$\end{document} production cross section at midrapidity, computed using the minimum-bias data sample, is dσ/dy|y=0=8.97±0.24(stat)±0.48(syst)±0.15(lumi)μb\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {d}\sigma /\text {d}y|_{y=0} = 8.97\pm 0.24~(\text {stat})\pm 0.48~(\text {syst})\pm 0.15~(\text {lumi})~\mu \text {b}$$\end{document}. 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 pT\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_{\mathrm{T}}$$\end{document}-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 probe-mail: alice-publications@cern.ch abilities 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 lowp 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 invariantmass 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 pre-resonant 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 highp 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-of-mass energy √ s = 13 TeV. The inclusive J/ψ yields include contributions from directly produced J/ψ, feed-down from prompt decays of highermass 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 Sect. 2, and the data analysis and the determination of the systematic uncertainties are described in Sects. 3 and 4, respectively. The results are presented and discussed together with recent model calculations in Sect. 5 and conclusions are drawn in Sect. 6.

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 highp 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 data-taking 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 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 Tdifferential 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 matched track of 0.8 < E/ p < 1.3 and on the dE/dx in the TPC of [−2.25, 3] σ . Due to the additional use of the EMCal for electron identification with respect to the MB based analysis, no explicit hadron rejection was used for the EMCal sample.
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 detail 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 oppositesign (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 opposite-sign 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 highp T ( p ee T > 1 GeV/c), plus a template shape obtained from MC simulations, described in Sect. 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 particleidentification (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 timedependent 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 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 [

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 con-tributions are added in quadrature and provided as ranges for the p T -and y-differential results.
As described in Sect. 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 nonnegligible 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 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) 1  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 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 dijet production generated by PYTHIA8 at √ s = 13 TeV mentioned in Sect. 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 Sect. 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 varia-tion of efficiency with respect to the central value is taken as a systematic uncertainty. This amounts to 1% for the p Tintegrated 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.

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 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 func-tion 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 feeddown from beauty decays in the theory predictions, a calculation of the non-prompt J/ψ production by Cacciari et al., using the Fixed-Order Next-to-Leading-Logarithms approach (FONLL) [45], is shown in the same panel. Accord-  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/ψ com-ponent 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) ing to the FONLL calculations, the non-prompt contribution to the inclusive J/ψ yield is approximately 10% at lowp 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.
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 PEGASUS [43] developed by Lipatov et al. which employs a k T -factorization approach using p Tdependent 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 Tfactorization 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 lowp T range.   [48], STAR [49], CDF [50] and ALICE [18][19][20] collaborations. The forward-rapidity production cross section shown in the right panel is reported by ALICE in Ref. [21] The NRQCD calculation from Lipatov, which uses the k Tfactorization 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 EMCaltriggered 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], STAR [49], CDF [50] and ALICE [18][19][20] is shown in the left panel of Fig. 4. An approximate logarithmic increase of the cross sections with the energy is observed. The collision energy-dependent measurements are also compared with the calculations from the ICEM [42] and the NRQCD + CGC [44] models. Both calculations provide a good description of the energy dependent trend within large theoretical uncertainties, dominated by the lowp T region of the spectrum.
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 highp 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 Tfactorization 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. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 . X. Zhang 7 , Y. Zhang 131 , V. Zherebchevskii 115 , Y. Zhi 11 , N. Zhigareva 95 , D. Zhou 7 , Y. Zhou 92 , J. Zhu 7,110 , Y. Zhu 7 , A. Zichichi 26 , G. Zinovjev 3 , N. Zurlo 59,142