Centrality dependence of inclusive J/$\psi$ production in p-Pb collisions at $\sqrt{s_{\rm NN}}$ = 5.02 TeV

We present a measurement of inclusive J/$\psi$ production in p-Pb collisions at $\sqrt{s_{\rm NN}}$ = 5.02 TeV as a function of the centrality of the collision, as estimated from the energy deposited in the Zero Degree Calorimeters. The measurement is performed with the ALICE detector down to zero transverse momentum, $p_{\rm T}$, in the backward ($-4.46<y_{\rm cms}<-2.96$) and forward ($2.03<y_{\rm cms}<3.53$) rapidity intervals in the dimuon decay channel and in the mid-rapidity region ($-1.37<y_{\rm cms}<0.43$) in the dielectron decay channel. The backward and forward rapidity intervals correspond to the Pb-going and p-going direction, respectively. The $p_{\rm T}$-differential J/$\psi$ production cross section at backward and forward rapidity is measured for several centrality classes, together with the corresponding average $p_{\rm T}$ and $p^2_{\rm T}$ values. The nuclear modification factor, $Q_{\rm pPb}$, is presented as a function of centrality for the three rapidity intervals, and, additionally, at backward and forward rapidity, as a function of $p_{\rm T}$ for several centrality classes. At mid- and forward rapidity, the J/$\psi$ yield is suppressed up to 40% compared to that in pp interactions scaled by the number of binary collisions. The degree of suppression increases towards central p-Pb collisions at forward rapidity, and with decreasing $p_{\rm T}$ of the J/$\psi$. At backward rapidity, the $Q_{\rm pPb}$ is compatible with unity within the total uncertainties, with an increasing trend from peripheral to central p-Pb collisions.


Introduction
Charmonia, bound states of charm and anti-charm quark pairs, are extensively used to study the interplay between the perturbative and the non-perturbative regimes of Quantum ChromoDynamics (QCD) [1]. Charmonium production mechanism can be understood as a hard scattering, describing the charm anticharm quark pair production, followed by the evolution of the pair into a bound state via a non-perturbative process. Models such as colour evaporation (CEM) [2,3], colour singlet (CSM) [4] and non-relativistic QCD (NRQCD) [5] are used to describe the charmonium production in hadronic collisions. None of these models has so far provided a consistent description of the production cross section and polarisation measured in proton-proton (pp) collisions [1,6]. The 1S vector state, the J/ψ meson, is abundantly produced in hadronic collisions at high energy and measurable through its leptonic decays. Its inclusive production contains contributions from direct J/ψ, from decays of higher-mass excited states, ψ(2S) and χ c , as well as from non-prompt J/ψ, from weak decays of beauty hadrons.
In proton-nucleus (p-A) collisions, several effects related to the nuclear medium and commonly denoted as cold nuclear matter effects (CNM) can affect the production of charmonia. The Parton Distribution Functions (PDFs) of nucleons bound in nuclei are modified compared to those of free nucleons [7][8][9]. These functions depend, in particular, on the fraction of the nucleon momentum, Bjorken-x (x Bj ), carried by the probed parton. In the collision energy regime typical of the Large Hadron Collider (LHC), charm quark pairs are produced mainly via the gluon fusion process. The gluon nuclear PDFs (nPDFs) are suppressed at low x Bj (x Bj 0.01), enhanced at intermediate x Bj (0.01 x Bj 0.3) and suppressed again at large x Bj (0.35 x Bj 0.7) compared to those of free nucleons. These three kinematic regions are often referred to as the shadowing, anti-shadowing, and EMC regions, respectively. Alternatively, at low x Bj , the initial colliding nucleus can be described by the Colour Glass Condensate (CGC) effective theory [10,11] as a coherent and dense (saturated) gluonic system. The kinematic distribution of the produced charm quark pairs may be additionally modified by multiple scattering of the incoming gluons and/or the quark pairs with the surrounding nuclear medium [11,12] or by the energy loss via gluon radiation [13,14]. It was also argued that the interference between the gluons radiated before and after the hard scattering can lead to important coherent energy loss effects at large rapidity in the p-going direction [15]. Finally, after their formation, the pre-resonant charm quark pairs or the fully formed resonances may interact with the nucleons when passing through the nucleus (nuclear absorption [16]) or with the other particles produced in the p-Pb collision (comovers [17]). Consequently, they may lose energy or fragment into open charm meson pairs. Due to the short time spent by the charm quark pairs in the nucleus relative to the J/ψ formation time at LHC energies, the effect of nuclear absorption is expected to be small [18,19].
Charmonium production was predicted to be suppressed in a hot medium with a high density of colour charges, the Quark-Gluon Plasma (QGP), as a consequence of the colour screening mechanism [20]. Such a state can be formed in ultra-relativistic nucleus-nucleus collisions. At the LHC, where the charmquark density is large, charmonium may also be created via the (re)combination of charm quarks either during the deconfined phase [21] or at the phase boundary [22], when the system has cooled down and hadronisation takes place. A suppression of J/ψ production in central nucleus-nucleus (A-A) collisions with respect to the one measured in pp collisions scaled by the number of binary nucleon-nucleon collisions was observed at the SPS at √ s NN ∼ 20 GeV [23][24][25], at RHIC at √ s NN = 39, 62.4 and 200 GeV [26][27][28][29] and at the LHC at √ s NN = 2.76 TeV [30][31][32]. However, the J/ψ production measurements at the LHC show a much smaller suppression of the yields integrated over transverse momentum, p T , as compared to the results at lower collision energies. The differential results also indicate a smaller degree of suppression at p T < 3 GeV/c than at higher p T [33,34], at mid-and forward rapidity, in agreement with the expectations from (re)combination models [35,36]. Although, qualitatively, the Pb-Pb measurements by themselves give a strong indication that the (re)combination effect plays a significant role in the J/ψ production at LHC energies, the quantitative understanding of the involved mechanisms The J/ψ production in proton-or deuteron-nucleus collisions was studied at fixed-target (SPS [37,38], HERA [39], Tevatron [40]) and collider experiments (RHIC [41], LHC [42][43][44][45][46]). At the LHC, a suppression of the J/ψ production in p-Pb collisions with respect to binary-scaled pp production has been observed for p T < 5 GeV/c at large rapidity in the p-going direction and at mid-rapidity, while the measurements at high p T as well as at large rapidity in the Pb-going direction are consistent with no suppression. The results are in fair agreement with models based on shadowing or coherent energy loss. While the J/ψ suppression at forward rapidity is overestimated by an early CGC calculation [47], recent calculations [48,49] are in better agreement with the data. The various CNM effects described above should be enhanced at small impact parameters of the collision and thus towards the most central p-Pb collisions. Hence, differential measurements as a function of the p-Pb collision centrality are essential to further constrain the models, in particular their dependence on the impact parameter of the collision.
In this paper, we report on new results in p-Pb collisions at √ s NN = 5.02 TeV for inclusive J/ψ production, measured at backward (−4.46 < y cms < −2.96) and forward (2.03 < y cms < 3.53) center-of-mass rapidity, y cms , in the µ + µ − decay channel, and at mid-rapidity (−1.37 < y cms < 0.43) in the e + e − decay channel. Previous measurements have been carried out as a function of rapidity and p T [42-46]. Here, the measurements are performed as a function of the collision centrality, estimated on the basis of the energy deposited in the Zero Degree Calorimeters (ZDC) [50]. At backward and forward rapidity, the J/ψ cross section is studied as a function of p T for several centrality classes. The corresponding average values p T , and p 2 T , are extracted from the p T -differential cross sections and the p T broadening, defined as ∆ p 2 T = p 2 T pPb − p 2 T pp , is also discussed. The nuclear modification factors are then obtained as a function of centrality for the three rapidity ranges and at backward and forward rapidity, as a function of p T for several classes of centrality.

Detectors and data sets
The ALICE apparatus and its performance are described in detail in Ref.
Away from the mid-rapidity region, the J/ψ candidates are reconstructed in the µ + µ − decay channel using the muon spectrometer [51], covering the pseudorapidity range −4 < η lab < −2.5 in the laboratory frame. The muon spectrometer includes a dipole magnet with an integrated field of 3 T·m, five tracking stations comprising two planes of Cathode Pad Chambers each, and two trigger stations consisting of two planes of Resistive Plate Chambers each. A system of absorbers is used for filtering out the hadrons. The front absorber, made of concrete, carbon and steel with a thickness of 4.1 m (10 nuclear interaction lengths, λ int ) is installed between the interaction region and the muon tracking stations. A second absorber, a 1.2 m thick iron wall (7.2 λ int ), is located upstream of the trigger stations and absorbs secondary hadrons escaping from the front absorber and low-momentum muons produced predominantly from π and K decays. Finally, a conical absorber placed around the beam pipe protects the spectrometer from secondary particles produced in interactions of large-η primary particles with the beam pipe.
At mid-rapidity, the J/ψ candidates are measured in the e + e − decay channel with the central barrel detectors in the pseudorapidity range |η lab | < 0.9. The main subsystems used, the Time Projection Chamber (TPC) [53] and the Inner Tracking System (ITS) [54], are placed in a solenoidal magnetic field with a strength of 0.5 T. The TPC, the main tracking and particle identification device, is a gaseous drift detector with a cylindrical geometry extending from 85 to 247 cm in the radial direction and 500 cm longitudinally. The particle identification is performed via the measurement of the specific energy loss, dE/dx, in the gas volume.
The ITS, covering a pseudorapidity range |η lab | < 0.9, consists of 6 layers of silicon detectors placed at radii ranging from 3.9 to 43 cm relative to the beam axis. The two innermost layers are equipped with Centrality dependence of inclusive J/ψ production in p-Pb at √ s NN =5.02 TeV ALICE Collaboration Silicon Pixel Detectors (SPD). The track segments (tracklets) reconstructed from the hits in the two SPD layers are used to reconstruct the interaction vertex position and to reject pile-up events (events with two or more simultaneous interactions per bunch crossing). The position of the interaction vertex is also determined, with better resolution, from the tracks reconstructed in the TPC and the ITS [51].
Two scintillator arrays, V0 [55], placed on both sides of the interaction point (IP) at −3.7 < η lab < −1.7 and 2.8 < η lab < 5.1, are used as trigger detectors and to remove beam-induced background. They are also used for the measurement of luminosity, along with the T0 detector [51], consisting of two quartz Cherenkov counters, placed on each side of the IP covering the ranges −3.3 < η lab < −3.0 and 4.6 < η lab < 4.9. The Zero Degree Calorimeters (ZDC) [56], located along the beam axis at 112.5 m from the IP on both sides, detect protons and neutrons emitted from the nucleus and are used to estimate the centrality of the collision. The neutron calorimeter (ZN) is positioned between the two beam pipes downstream of the first machine dipole that separates the beams. The proton calorimeter (ZP) is installed externally to the outgoing beam pipe. The ZDCs are also used to remove parasitic p-Pb interactions displaced from the nominal position.
The data samples used for the measurements reported in this paper were collected in 2013 in two configurations, obtained by inverting the direction of the p and Pb beams. Due to the asymmetry of the energy per nucleon of the p and Pb beams (E p = 4 TeV and E Pb /208 = 1.58 TeV), the nucleon-nucleon center-of-mass system is shifted with respect to the laboratory system by ∆y = 0.465 in the p-going direction. The two beam configurations allow one to measure the J/ψ production in the backward (−4.46 < y cms < −2.96) and forward (2.03 < y cms < 3.53) centre-of-mass rapidity (y cms ) regions, corresponding to the Pb-going and the p-going directions, respectively. They will be further referred to as Pb-p and p-Pb, for the first and second case. The dielectron analysis was carried out on a data sample corresponding to the p-Pb beam configuration, in the mid-rapidity range −1.37 < y cms < 0.43.
The dielectron analysis was performed on a sample of events satisfying a Minimum Bias (MB) trigger condition and the dimuon analysis used dimuon-triggered events. The MB trigger is defined by a coincidence of the signals from both sides of the V0 detector. The efficiency of the MB trigger in selecting non-single diffractive p-Pb collisions was estimated to be higher than 99% [57], with negligible contamination from diffractive collisions. The dimuon trigger requires, in addition to the MB condition, the detection of two unlike-sign muon candidate tracks in the trigger system of the muon spectrometer. This trigger selects muons with a transverse momentum p T 0.5 GeV/c. This threshold is not sharp in p T and the single-muon trigger efficiency reaches a plateau value of ∼ 96% at p T ∼ 1.5 GeV/c.
In the backward and forward rapidity regions, the measurements are based on a sample of 2.1 × 10 7 and 9.3 × 10 6 dimuon-triggered events, respectively. The MB interaction rate reached a maximum of 200 kHz, corresponding to a maximum pile-up probability of about 3%. The mid-rapidity data sample consists of 1.1 × 10 8 MB-triggered events, collected at a low interaction rate (∼ 10 kHz) and with a fraction of pile-up events lower than 0.6%. Pile-up of collisions from different bunch crossings is negligible considering that the 200 ns bunch-crossing spacing is larger than the integration time of the ZDC and muon trigger detectors, which are used for the track selection for the dimuon analysis. Two independent determinations of the MB trigger cross sections σ MB were carried out in the Pb-p and p-Pb configurations using van der Meer scans [58]. The corresponding cross sections amount to σ Pbp MB = 2.12 ± 0.07 b and σ pPb MB = 2.09 ± 0.07 b, respectively [59]. The integrated luminosity is determined as L = N MB /σ MB where N MB is the number of MB events. The number of MB events corresponding to the dimuontriggered sample is evaluated as N MB = F 2µ/MB · N DIMU , where N DIMU is the number of dimuon-triggered events and F 2µ/MB is the inverse of the probability of having dimuon-triggered events in a MB data sample. The determination of F 2µ/MB is discussed in Sec. 3. The integrated luminosity was also independently measured employing the T0 detector. The two luminosity measurements agree within better than 1% throughout the whole Pb-p and p-Pb data-taking periods [59]. The maximum difference is included as an additional uncertainty for σ MB and thus in the determination of the luminosity uncertainty.
Centrality dependence of inclusive J/ψ production in p-Pb at √ s NN =5.02 TeV ALICE Collaboration The integrated luminosity values used for the results in the backward, forward and mid-rapidity regions are 5.81 ± 0.20 nb −1 , 5.01 ± 0.19 nb −1 and 51.4 ± 1.9 µb −1 , respectively.
The centrality determination in Pb-Pb collisions is usually based on charged-particle multiplicity, estimated using the V0 signal amplitudes [60]. However, in p-Pb collisions, the magnitude of the multiplicity fluctuations at a given impact parameter is comparable to the whole dynamic range of the MB multiplicity distribution. The fluctuations can be related to the various event topologies (e.g. hard collisions with large momentum transfers and/or multiple hard parton-parton interactions, which tend to be associated to high multiplicity events, compared to soft collisions without any high-p T particle), detector acceptance effects (jets fragmenting in or out the experimental coverage), or other effects, as explained in detail in Ref.
[50]. Therefore, a centrality selection based on charged-particle multiplicity may select a sample of p-Pb collisions that contains biases unrelated to the collision geometry. In contrast, a centrality selection based on the energy measured with the ZDC in the Pb-going direction, deposited by nucleons produced in the nuclear de-excitation processes following the collision, or knocked out by wounded nucleons, should not induce such biases. The average number of binary nucleon collisions ( N coll ) or the average nuclear overlap function ( T pPb ) for a given centrality class, defined by a selected range of energy deposited in the Pb-remnant side of ZN, is obtained using the hybrid method described in Ref. [50].
In this method, the N coll determination relies on the assumption that the charged-particle multiplicity measured at mid-rapidity is proportional to the number of participant nucleons ( N part ). The values of N part for a given ZN-energy class, also noted as ZN class in the following, were calculated by scaling the MB value of the number of participant nucleons, N MB part , by the ratio of the average charged-particle multiplicities measured at mid-rapidity for the considered ZN-energy event classes to the corresponding value in MB collisions. The average number of collisions and the average nuclear overlap function were then calculated from N part according to the Glauber model [61], which is generally used to calculate geometrical quantities of nuclear collisions. From here on these values are denoted as N mult coll and T mult pPb to indicate the ansatz used for their derivation. Other assumptions to derive N coll , which are discussed in [50], use the proportionality of N coll to the yield of high-p T charged particles (10 < p T < 20 GeV/c) at mid-rapidity or to the charged-particle multiplicity measured with the V0 detector in the Pb-going direction at forward rapidity. The variations on the N coll values obtained with the three methods do not exceed 6% for any of the centrality classes used for this analysis and are taken into account as a systematic uncertainty uncorrelated over centrality. Uncertainties of 8% and 3.4% on the determination of N MB coll and T MB pPb , respectively, are also included as global systematic uncertainties. These uncertainties are obtained by varying the parameters of the Glauber model. Events without a signal in the ZN detector, which correspond to very peripheral events [50], are assigned to the 80-100% centrality interval. The values of N mult coll and T mult pPb used in this analysis 1 are reported in Tab. 1, together with their uncertainties.

Analysis in the dimuon decay channel
The analysis approach and the selection criteria are similar to those described in detail in Ref. [42]. The primary vertex is reconstructed from the hits in the SPD. No specific requirement is applied on the vertex properties. In pile-up events, the ZN energy of the two (or more) interactions are summed up, increasing the pile-up event contribution in the most central ZN class. This contribution is estimated to be large, of the order of 20-30%, for events belonging to the centrality class 0-2% and this class has therefore been discarded from the analysis. The muon candidate tracks are reconstructed in the muon spectrometer tracking stations using the algorithm described in Ref. [62]. In order to remove particles at the edge of the muon spectrometer acceptance, a fiducial cut on the single-muon pseudorapidity −4 < η lab < −2.5 is applied. An additional selection on the radial coordinate of the track at the exit of the front absorber (17.6 < R abs < 89.5 cm) is required to reject muons crossing the high-density section of the absorber, √ s NN =5.02 TeV ALICE Collaboration ZN class N mult coll T mult pPb 2-10% 11.7 ± 1.2 ± 0.9 0.167 ± 0.012 ± 0.006 10-20% 11.0 ± 0.4 ± 0.9 0.157 ± 0.006 ± 0.005 20-40% 9.6 ± 0.2 ± 0.8 0.136 ± 0.003 ± 0.005 11.4 ± 0.6 ± 0.9 0.164 ± 0.009 ± 0.006 60-100% 3.2 ± 0.2 ± 0.3 0.046 ± 0.002 ± 0.002 Table 1 where energy loss and multiple scattering effects play an important role. Finally, only tracks matching the corresponding tracks reconstructed in the trigger stations are selected.
The J/ψ candidates are obtained by combining pairs of muons of opposite charge that are reconstructed in the rapidity range 2.5 < |y lab | < 4 and with p T < 15 GeV/c. The raw J/ψ yield is estimated for each centrality and p T interval from fits of the dimuon invariant-mass distribution performed with various functions. For the signal component, an extended Crystall Ball function, which includes non-Gaussian tails on either side of the J/ψ peak, as well as a pseudo-Gaussian function with a mass-dependent width [63] are employed. Due to the poor signal-to-background (S/B) ratio in the tail regions, Monte Carlo (MC) simulations are used to constrain the tail parameters in each p T and rapidity interval under study. Since there is no degradation of the tracking resolution due to the large occupancy corresponding to the most central p-Pb collisions [52], the tails are not expected to depend on centrality. The ψ(2S) resonance is also included in the fit function using the strategy described in Ref. [64]. For the background component, two alternative functions are used: a Gaussian with a mass-dependent width and an exponential multiplied by a second-order polynomial. The fits are performed in two different invariant mass ranges, 2 < m µ + µ − < 5 GeV/c 2 and 2.3 < m µ + µ − < 4.7 GeV/c 2 . In the fitting procedure, the mean and width of the J/ψ signal function, the background parameters and the normalisation factors are left free while the tail parameters are fixed to the values estimated from the simulations. The obtained J/ψ mass value agrees with the PDG value [65] within 5 MeV/c 2 . The measured width increases from 59 to 81 MeV/c 2 with increasing p T . It is found to be about 10% larger than in the simulations. The S/B ratio in the 3σ interval around the J/ψ pole increases with increasing p T and towards peripheral events, ranging from 1.1 (1.3) to 5.2 (9) in the Pb-p (p-Pb) configuration. Figure 1 shows examples of fits to the unlike-sign dimuon pair invariant mass distributions for the Pb-p configuration for six centrality classes for p T < 15 GeV/c.
The invariant mass fits are performed using the different combinations of signal and background functions and fitting ranges described above. The number of J/ψ is obtained by integrating the J/ψ signal function over the fitting range. The mean of the distribution of the number of J/ψ obtained from the various fits is used as the central value of the raw yield, while the Root Mean Square (RMS) is used as a systematic uncertainty, which ranges between 0.2% and 3% depending on p T and centrality class. An additional systematic uncertainty of 2% is added to the signal extraction uncertainty. It is estimated from the variation of the raw yield when performing the fit with different tail parameters of the signal function. The raw J/ψ yield varies between about 3000 and 16000 counts for Pb-p and between about 5000 and 17000 counts for p-Pb in the centrality-differential results. In the case of the centrality-and p T -double-differential results, it varies between about 70 and 4000 counts in Pb-p and between about 150 and 4000 counts in p-Pb, where the lower values correspond to the centrality class 80-100% and the highest p T range. √ s NN =5.02 TeV ALICE Collaboration  The J/ψ raw yields are corrected for the detector acceptance and efficiency (A × ε) estimated from simulations of the J/ψ signal. The muon decay products of the J/ψ are propagated through the experimental setup modeled with GEANT 3.21 [66]. The procedure used for track reconstruction is the same in data and simulations. In the latter, the detector conditions and their variation with time during the data-taking period are taken into account. It was checked that the detector occupancy in central collisions does not deteriorate the single muon tracking efficiency and resolution, which justifies that only the J/ψ signal is simulated and not the underlying p-Pb collision. The p T and rapidity distributions of the J/ψ signal in the simulation were tuned to the reconstructed distributions of the p-Pb and Pb-p data using an iterative procedure. The J/ψ production is assumed to be unpolarised, consistent with the observation that no significant J/ψ polarisation has been measured in pp collisions at √ s = 7 TeV [67 -69]. The values of A × ε integrated over p T are 17.1% and 25.4% in the Pb-p and p-Pb configurations, respectively. The lower A × ε for the Pb-p configuration is due to a smaller detector efficiency in the corresponding data-taking period. The A × ε varies as a function of p T from 16% to 33% in Pb-p and from 23% to 48% in p-Pb collisions, where the lowest values correspond to 1 < p T < 2 GeV/c and the largest to 8 < p T < 15 GeV/c. The systematic uncertainty in the choice of the J/ψ kinematic distributions in the simulation is estimated by varying the J/ψ p T and rapidity distributions according to the measured ones over various sub-ranges of y, p T and centrality (see [42] for more details). When integrated over p T , this uncertainty amounts to 1.5% for both the Pb-p and p-Pb configurations, while for the p T -differential studies it does not exceed 1.4%. The uncertainty on the dimuon tracking efficiency amounts to 6% (4%) for Pb-p (p-Pb). It is evaluated using the difference between the single-muon tracking efficiency obtained from simulations and a data-driven approach based on the redundancy of the muon tracking stations [52], assuming that the efficiencies of the two muons are uncorrelated. This uncertainty is correlated over centrality and is taken as constant as a function of p T . The uncertainty on the determination of the dimuon trigger efficiency has three contributions, which are correlated over centrality. The first uncertainty is due to the statistical uncertainty of the trigger detector efficiency which is estimated using data. It is independent of p T and amounts to 2%. The second uncertainty is extracted from the differences observed between data and simulations for the measured trigger response function in the region close to the trigger threshold. This uncertainty varies between 0.5% and 3% and is larger at low p T . The third uncertainty is due to the small Centrality dependence of inclusive J/ψ production in p-Pb at √ s NN =5.02 TeV ALICE Collaboration fraction of opposite-sign pairs which are misidentified as like-sign by the trigger system and increases from 0.5 to 3% with increasing p T . An additional systematic uncertainty results from the choice of the χ 2 cut, which is applied to the matching of tracks reconstructed in the muon tracking and trigger system.
Applied to the number of dimuon pairs, this uncertainty amounts to 1% and is correlated over centrality.
The normalisation factor of dimuon-to MB-triggered events, F 2µ/MB , which is needed to evaluate the integrated luminosity, is determined in a two-step procedure as the product F 2µ/1µ · F 1µ/MB , where F 2µ/1µ (F 1µ/MB ) is the inverse of the probability of having dimuon-triggered events (single-muon-triggered events) in a corresponding data sample of single-muon-triggered events (MB-triggered events). The various quantities are estimated from the number of recorded triggered events in each centrality class. This factor can also be obtained from the centrality-integrated value scaled by N cent MB /N MB and N DIMU /N cent DIMU , the fraction of MB events and the inverse of the fraction of dimuon events in a given centrality class, respectively. The latter method, which is statistically more accurate, is used for the evaluation of F 2µ/MB . The systematic uncertainty is evaluated from the comparison between the two methods and amounts to 1 − 2% depending on the centrality class. The value of F 2µ/MB depends on the centrality class and it smoothly increases from 260 and 3340 in Pb-p and from 660 and 3290 in p-Pb from central to peripheral collisions. The pile-up event contribution to the number of MB events is estimated by using alternatively the interaction vertices reconstructed with the SPD to select pile-up events, or a fast simulation describing the ZN energy distribution. The pile-up event contribution is larger in the 2-10% centrality class where it amounts to 3.5% and 2.7% in the Pb-p and p-Pb beam configurations, respectively. It decreases to less than 2% in all other centrality classes. It has been included in the systematic uncertainty of F 2µ/MB . It was further checked by using the fast simulation that the overall effect of pile-up events, including the shift of events from a given centrality class to a more central one, is covered by the systematic uncertainties quoted for pile-up events.
In order to quantify the nuclear effects in p-Pb collisions, reference measurements in pp collisions at the same energy are needed. Since there are no experimental data available on the J/ψ production cross section in pp collisions at √ s = 5.02 TeV, the procedures described in Ref. [70] for the p T -integrated case and in Ref.
[44] for the p T -differential case are used. These procedures involve an interpolation in energy and an extrapolation in rapidity and are based on existing measurements in pp collisions at different energies. The resulting values of the J/ψ cross section interpolated to √ s = 5.02 TeV are also reported in those references.

Analysis in the dielectron decay channel
The analysis method and the selection criteria are similar to those described in detail in Ref. [44]. Events are selected which contain a primary vertex determined from tracks reconstructed in the TPC and the ITS. The vertex position is required to have a distance to the nominal IP smaller than 10 cm along the beamline. Due to the lower interaction rate in the data sample used for the mid-rapidity analysis, the pileup event contribution is negligible and the events belonging to the ZN centrality class 0-2% are included in the analysis. The electron and positron candidate tracks are reconstructed in the pseudorapidity range |η lab | < 0.9. Tracks are required to have at least 70 out of a maximum of 159 clusters in the TPC, a χ 2 normalised to the number of attached TPC clusters smaller than 4 and a distance of closest approach to the primary vertex smaller than 3 cm along the beam axis and 1 cm in the plane transverse to the beam axis. Only tracks with one associated track point in the innermost layer of the ITS and at least a second one in the other layers are selected. This requirement suppresses the secondary electrons from photon conversions in the detector material of the ITS. Tracks are required to be compatible within 3.0σ with the electron hypothesis based on the measured ionisation energy loss of the TPC. In order to reject hadrons, the tracks which are compatible within 3.5σ with the pion or proton hypotheses are excluded. Tracks from identified photon conversions are rejected without any impact on the J/ψ signal efficiency. Finally, tracks are required to have a transverse momentum larger than 1 GeV/c in order to improve the S/B ratio Centrality dependence of inclusive J/ψ production in p-Pb at √ s NN =5.02 TeV ALICE Collaboration The invariant mass distribution of e + e − pairs obtained from this procedure is normalised to the integral of the same-event unlike-sign dielectron pairs in the mass ranges 2.0 − 2.5 GeV/c 2 and 3.2 − 3.7 GeV/c 2 , outside of the signal counting interval. A significant fraction of the J/ψ yield, determined from simulations to be about 30%, falls outside the signal counting window, and is corrected for. This is due to the long tail at low masses caused by the bremsstrahlung of electrons in the detector material and by the radiative soft photon that may be emitted at the J/ψ decay vertex. The systematic uncertainty on the signal extraction procedure, including uncertainties on the mixed-event background scaling and on the invariant-mass shape of the e + e − decay channel, is obtained by varying the mass region used for the scaling of the mixed-event background and by varying the mass window used for counting the signal. Figure 2 illustrates the signal extraction procedure for the centrality classes considered in this analysis. The central value of the J/ψ raw yield is obtained as the average of the raw yields retrieved from the variation of the signal extraction configurations described above and its systematic uncertainty is the RMS of the distribution of the extracted signals. The raw yield has also been evaluated using the like-sign method where the residual background after the like-sign subtraction is estimated by a linear function and the signal is extracted by bin counting. The two methods have been found to provide results compatible within the estimated uncertainties. As a function of centrality, the J/ψ raw yield varies from 73 to 133 counts, with S/B in the interval 2.92 < m e + e − < 3.16 GeV/c 2 ranging from 0.7 to 1.6.
The A × ε correction is estimated with simulations consisting of J/ψ particles set to decay into an e + e − pair added to a p-Pb event generated using HIJING [71]. The produced particles are subsequently propagated through the experimental setup modeled using GEANT 3.21 [66] and the same reconstruction procedure as for the real data is followed. The simulated J/ψ mesons are assumed to be unpolarised. The average value obtained for A × ε is 7.2% with no observed dependence on the collision centrality but with a significant dependence on the J/ψ p T , having a maximum of ∼ 11% at zero p T , a minimum of ∼ 6% at 2 GeV/c, a second maximum of ∼ 10% at 7 GeV/c followed by a slow decrease towards higher momenta. This shape is due to the kinematic selections and the momentum dependence of the particle identification selection efficiency.
The uncertainty on the J/ψ reconstruction efficiency is dominated by the uncertainty on the electron identification. It is estimated from the difference in the TPC specific energy loss distribution of a clean Centrality dependence of inclusive J/ψ production in p-Pb at √ s NN =5.02 TeV ALICE Collaboration sample of electrons from identified photon conversions in data and electrons from the simulation. For dielectron pairs, this uncertainty amounts to 4% and is correlated over centrality. Since the A × ε is not constant as a function of p T , its p T -integrated value depends on the p T shape used to generate the J/ψ mesons in the simulation. To estimate the uncertainty due to the generated p T shape, the p T -differential J/ψ spectrum is varied in the simulations. The p T distribution is parameterised as where C, p 0 and n are parameters constrained by the experimental results in Ref.
[44]. The systematic uncertainty is estimated from the variation of the p T -integrated A × ε when varying the fit parameters within their uncertainties. It amounts to 3% and is correlated among the centrality classes. Furthermore, an uncertainty due to the variation of the spectral shape as a function of centrality is also taken into account, given the fact that statistical uncertainties do not allow for a double-differential measurement.
To evaluate this effect, the relative variation of the fit parameters of the p T distributions measured at forward and backward rapidity for several centrality classes is used. A maximum variation of the A × ε by 1.4% is observed and assigned as an uncorrelated uncertainty over centrality.
The inclusive J/ψ production cross section in pp collisions at √ s = 5.02 TeV, which is needed to quantify the nuclear effects in p-Pb collisions, is obtained using the interpolation procedure described in Ref.
[44]. The method is based on existing measurements in pp collisions at different energies.

Results
The double-differential J/ψ production cross section for a given centrality class is where σ MB is the p-Pb (Pb-p) MB cross section discussed in Sec. 2, BR is the branching ratio of the considered J/ψ dileptonic decay channel, which amounts to (5.96 ± 0.03)% and (5.97 ± 0.03)% for the dimuon and the dielectron decay channels [65], respectively, and Y cent J/ψ→l + l − is the inclusive J/ψ yield per-event. The latter is defined as where N J/ψ→l + l − is the raw number of J/ψ mesons decaying into dileptons for a given centrality class, rapidity and p T range, N MB is the number of MB events for the given centrality class, A × ε is the acceptance times efficiency described in Sec. 3 and 4 and ∆y and ∆p T are the widths of the rapidity and p T intervals, respectively. Table 2 gives a summary of the systematic uncertainties of the J/ψ differential cross section, as well as the correlations of these uncertainties over centrality, collision system and J/ψ p T . The p T -integrated J/ψ cross sections are reported in Tab. 3 for the three rapidity intervals as a function of centrality expressed in percentiles of the non-single diffractive p-Pb cross section. Figure 3 shows the double-differential J/ψ cross sections as a function of p T in the range 0 < p T < 15 GeV/c at backward (left panel) and forward (right panel) rapidity measured for six centrality classes. The vertical error bars represent the statistical uncertainties and the open boxes the systematic uncertainties. The systematic uncertainties correlated over centrality and p T are indicated as a global relative systematic uncertainty. In order to characterise the evolution of the p T -differential cross section shape with centrality, the average values p T and p 2 T , were extracted for each centrality class by performing a fit to the data with the function defined in Eq. 1. The systematic uncertainties on the data points that are correlated over p T are not considered in the fit. The uncertainties on the free parameters obtained Centrality dependence of inclusive J/ψ production in p-Pb at √ s NN =5.  175 ± 4 ± 4 ± 14 220 ± 5 ± 5 ± 14 Table 3: Differential cross sections as a function of centrality. The first quoted uncertainty is statistical while the second and third represent the systematic uncertainties, the latter being fully correlated over centrality.
from the fit of Eq. 1 are propagated to the values of p T and p 2 T . The statistical and systematic uncertainties on p T and p 2 T are obtained by performing the fit using, separately, only the statistical or the uncorrelated systematic uncertainties on the data points, respectively. The range of integration over p T used to compute p T and p 2 T is limited to the measured p T interval 0 < p T < 15 GeV/c. It was verified that extending the integration range to infinity results in an increase of p T and p 2 T values by less than 0.5%. The values of p T and p 2 T obtained for each centrality class are reported in Tab. 4. Both p T and p 2 T values increase with centrality, which indicates a hardening of the p T distributions from peripheral to central collisions in both rapidity intervals.
In order to quantify the nuclear effects on the J/ψ p T spectrum shape, the p T broadening, ∆ p 2 T , defined as is used. Since there are no measurements for pp collisions at √ s = 5.02 TeV, the value of p 2 T pp is evaluated from the p T -differential cross section in pp collisions calculated with the interpolation procedure described in Ref. [44], and using the same p T -integration range as for p-Pb collisions. Figure 4 shows Centrality dependence of inclusive J/ψ production in p-Pb at √ s NN =5.  Table 4: Values of p T and p 2 T of inclusive J/ψ in the range 0 < p T < 15 GeV/c. The first quoted uncertainty is statistical while the second is systematic. The values obtained from the pp cross section interpolated to √ s = 5.02 TeV are also indicated.
∆ p 2 T as a function of the number of binary collisions. Our measurements indicate that ∆ p 2 T increases at backward (forward) rapidity by ∼ 1.1 GeV 2 /c 2 (∼ 1.4 GeV 2 /c 2 ) from peripheral to central p-Pb collisions. At forward rapidity, ∆ p 2 T is larger for all centrality classes and suggests a steeper dependence on centrality compared to backward rapidity values. For the most peripheral collisions, corresponding to N mult coll ∼ 2, the p 2 T value at backward rapidity is compatible with the one in pp collisions, while at forward rapidity it is found to be larger than in pp collisions by 0.7 GeV 2 /c 2 , which corresponds to 1.4 times the total uncertainty on the measured difference. The magnitude of the p T broadening observed by PHENIX [41] in d-Au collisions at √ s NN = 200 GeV in the rapidity ranges −2.2 < y cms < −1.2, |y cms | < 0.35 and 1.2 < y cms < 2.2 is similar to the one measured by ALICE at backward rapidity. At forward rapidity, the ALICE data show a stronger p T broadening and a steeper increase with increasing centrality as compared to PHENIX results.
The calculations from Refs. [72,73] are based on the leading order (LO) CEM production model and include initial and final-state multiple scattering of partons with the nuclear medium (denoted as "Mult. scattering" in Fig. 4). The uncertainties on the theoretical calculations are not available. In this model, the contribution to p T broadening due to final-state multiple scattering is expected to be sensitive to the colour-octet or colour-singlet nature of the pre-resonant cc pair. The calculations are in good agreement with the data at backward and forward rapidity. A second model, which is based on a parameterisation of the prompt J/ψ pp cross section and includes coherent energy loss effects from the incoming and outgoing partons [74] (denoted as "Eloss" in Fig. 4), describes well the centrality dependence of ∆ p 2 T at backward rapidity. The trend predicted by this model at forward rapidity is slightly steeper than the data even considering its uncertainty, evaluated by varying the gluon transport coefficient and the Centrality dependence of inclusive J/ψ production in p-Pb at √ s NN =5.02 TeV ALICE Collaboration parametrisation of the production cross section.
In order to study the modification of the J/ψ production in p-Pb collisions with respect to pp interactions, the J/ψ nuclear modification factor is used. For a given centrality class, rapidity and p T range, it is defined as The T mult pPb values corresponding to the centrality classes used in this analysis are reported in Tab. 1. The J/ψ cross section in pp collisions at √ s = 5.02 TeV, d 2 σ pp /dydp T , is obtained by means of the interpolation procedures outlined in Sec. 3 and 4 and described in Refs. [44,70]. The nuclear modification factor is usually denoted as R pPb but in this analysis the notation Q pPb is used to emphasise the possible bias in the evaluation of T mult pPb , as discussed in Sec. 2 and in Ref.
[50]. The systematic uncertainties on Q pPb are presented in Tab. 2. Figure 5 shows the dependence of the p T -integrated Q pPb on the collision centrality, expressed as N mult coll . The results for the backward, mid-and forward rapidity intervals are displayed in the left, middle and right panel, respectively. At backward rapidity, the Q pPb values show that the measured J/ψ production is compatible, within the total uncertainties, with expectations from binary collision scaling for all centrality classes. When considering only the uncertainties that are not correlated over centrality, an increase from peripheral to central p-Pb collisions is observed in the data. At forward rapidity, the J/ψ yield is suppressed with respect to the binary-scaled pp reference for all the considered centrality classes. The values of Q pPb measured at forward rapidity exhibit a decrease from 0.85 for the 80-100% centrality class down to 0.66 for the 2-10% centrality class. Within the present uncertainties, the mid-rapidity results suggest a similar degree of suppression of the J/ψ yield as at forward rapidity and no conclusion can be drawn on a possible centrality dependence.
Our measurements are compared to several theoretical models including a next-to-leading order (NLO) √ s NN =5.02 TeV ALICE Collaboration CEM calculation [16,75] which contains the EPS09 NLO nPDF parameterisation [7] (denoted as "CEM+EPS09 NLO" in Fig. 5), a model employing the EPS09 LO nPDF with or without effects from the interaction with a comoving medium [17] (denoted as "EPS09 LO+comovers" in Fig. 5), and the coherent energy loss model [74] described above. In the CEM+EPS09 NLO and EPS09 LO+comovers models, assuming the J/ψ production process is gg → J/ψ (2 → 1), the x Bj values of the gluon from the Pb nucleus span a range of about 1 · 10 −2 < x Bj < 5 · 10 −2 at backward rapidity, 4 · 10 −4 < x Bj < 2 · 10 −3 at mid-rapidity, and 2 · 10 −5 < x Bj < 8 · 10 −5 at forward rapidity. The backward rapidity interval therefore corresponds to the x Bj range in the transition between the anti-shadowing and the shadowing region, whereas the mid-and forward rapidity intervals probe a region for which the gluon shadowing is expected to be strong. The CEM+EPS09 NLO model uncertainties are evaluated from the EPS09 uncertainty, which gives the dominant contribution, and from a variation of the values of the charm quark mass, the normalisation and the factorisation scales in the pQCD calculation. The CEM+EPS09 NLO model reproduces well the centrality dependence in each rapidity range. At mid-and forward rapidity, the data are better reproduced when a strong shadowing is considered in the model. In the framework of the EPS09 LO+comovers model, the presence of a comoving medium has only a small effect on J/ψ production at forward rapidity since its density is expected to decrease towards the p-going direction. The effect of comovers is more pronounced at mid-rapidity and especially at backward rapidity and it increases with increasing centrality. The uncertainties on these theoretical calculations are not available. At backward rapidity, the increase of Q pPb towards central collisions observed in the data is better reproduced when the comover effect is not included in the model. Finally, the shape and magnitude of Q pPb is well described by the Eloss model in all rapidity intervals, although the model does not predict an increase with increasing centrality at backward rapidity, as indicated by the data.
It is worth pointing out that the calculations above are done for prompt J/ψ production, while the measurements also include the contribution of J/ψ mesons from b-hadron decays. The Q prompt pPb can be extracted from Q incl pPb using the relation Q prompt where f B is the ratio of non-prompt to prompt J/ψ production cross sections and Q non−prompt pPb is the nuclear modification factor of non-prompt J/ψ mesons. A value of f B of about 0.11 at 2 < y cms < 4.5 and for p T < 14 GeV/c can be calculated from the LHCb measurements in pp collisions at √ s = 7 TeV [76]. The value of f B does not show a strong variation within the quoted rapidity range and with energy, as indicated by the comparison with the results in pp collisions at √ s = 8 TeV [77]. Hence, the value of f B calculated at √ s = 7 TeV in 2 < y cms < 4.5 is used for the following. At mid-rapidity, a value of f B of about 0.17 at |y| < 0.9 and integrated over p T can be extracted from the measurements of ALICE in pp collisions at √ s = 7 TeV [78]. The nuclear modification factor of non-prompt J/ψ was measured to be 0.98 ± 0.06 ± 0.01 (0.83 ± 0.02 ± 0.08) for −4 < y cms < −2.5 (2.5 < y cms < 4) and p T < 14 GeV/c at √ s NN = 5.02 TeV in p-Pb collisions [43]. If the non-prompt J/ψ Q pPb , which has not been measured as Centrality dependence of inclusive J/ψ production in p-Pb at √ s NN =5.02 TeV ALICE Collaboration   [16,74,75]. √ s NN =5.02 TeV ALICE Collaboration a function of the centrality, is conservatively assumed to vary from 0.6 to 1.3 in each centrality interval, then the differences between the inclusive and prompt J/ψ nuclear modification factors cannot exceed 15% in any of the centrality classes and are smaller than the quoted uncertainties. Figure 6 presents the p T -dependence of Q pPb for the 2-10%, 10-20%, 20-40%, 40-60%, 60-80% and 80-100% centrality classes, from the top to the bottom panels, respectively. The left (right) panels show the backward (forward) rapidity results. At backward rapidity and for the most central collisions (2-10% and 10-20% centrality classes), Q pPb is compatible with unity in the full p T interval, and an increase of Q pPb from p T < 1 GeV/c to p T > 1 GeV/c is suggested by the data. For semi-central and peripheral collisions, Q pPb is compatible with unity over the full p T range. At forward rapidity, for all centrality classes other than the most peripheral one, J/ψ production is suppressed compared to binary-scaled pp production at low p T . For these centrality classes, Q pPb increases with increasing p T and is compatible with unity at high p T within the uncertainties. The magnitude of Q pPb slightly increases from central to semi-peripheral collisions over the full p T range. For the most peripheral collisions Q pPb is compatible with unity and, with the current uncertainties, it does not show a significant p T dependence.
The data are compared to the calculations from the CEM+EPS09 NLO [16,75] and the Eloss [74] model. The CEM+EPS09 NLO calculations describe reasonably well the Q pPb results at backward and forward rapidity. The Eloss model reproduces well the p T dependence of Q pPb at backward rapidity for all centrality classes. At forward rapidity, a good agreement of the calculations with the data is observed for peripheral collisions (60-80% and 80-100% centrality classes), while the p T dependence becomes steeper than in data towards more central collisions. It was observed in Ref.
[44] that a better agreement is reached with the p T dependence of the nuclear modification factor in centrality-integrated p-Pb collisions at forward rapidity when shadowing effects are included in the model.

Conclusions
The cross sections and nuclear modification factors, Q pPb , of inclusive J/ψ production have been measured with ALICE as a function of rapidity, p T and centrality in p-Pb collisions at √ s NN = 5.02 TeV.
For the most peripheral p-Pb collisions, no modification with respect to pp collisions is observed within the uncertainties of the measurements for both the shape of the p T spectrum of the J/ψ and the Q pPb measurements. On the contrary, the results in central p-Pb collisions suggest sizeable nuclear effects. At both backward (Pb-going direction) and forward (p-going direction) rapidity the ∆ p 2 T measurements show a p T broadening which increases monotonically from peripheral to central p-Pb collisions with larger values at forward rapidity. Our measurements show a stronger p T broadening and a steeper increase with increasing centrality at forward rapidity as compared to PHENIX results in d-Au collisions at √ s NN = 200 GeV [41]. At backward rapidity, a modest increase of J/ψ production compared to a binary-scaled pp reference is suggested by the data in most central collisions. At mid-rapidity, the data indicate that J/ψ production is suppressed compared to binary-scaled pp cross section over the entire centrality range. Within the current uncertainties, the increasing suppression towards central p-Pb collisions suggested by models is compatible with the data. Finally, at forward rapidity, a clear suppression, which increases towards central events, is observed. The p T -and centrality-differential results show that the suppression is stronger at low p T and tends to vanish at high p T . Given the uncertainties of both the measurements and the theoretical calculations, we observe a fair agreement of the models based on coherent energy loss and multiple scattering with the measured p T broadening. Models based on nPDF and coherent energy loss are in fair agreement with the nuclear modification factor measurements. The results presented in this paper provide an important baseline for understanding and constraining the cold nuclear matter effects in p-Pb collisions as well as their centrality dependence. Such an information is essential for a quantitative interpretation of the results obtained in Pb-Pb collisions.