Towards an automated tool to evaluate the impact of the nuclear modification of the gluon density on quarkonium, D and B meson production in proton–nucleus collisions

We propose a simple and model-independent procedure to account for the impact of the nuclear modification of the gluon density as encoded in nuclear collinear PDF sets on two-to-two partonic hard processes in proton–nucleus collisions. This applies to a good approximation to quarkonium, D and B meson production, generically referred to H\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {H}$$\end{document}. Our procedure consists in parametrising the square of the parton scattering amplitude, Agg→HX\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {A}}_{gg \rightarrow {\mathcal {H}} X}$$\end{document} and constraining it from the proton–proton data. Doing so, we have been able to compute the corresponding nuclear modification factors for J/ψ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J/\psi $$\end{document}, Υ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Upsilon $$\end{document} and D0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D^0$$\end{document} as a function of y 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_T$$\end{document} at sNN=5\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{s_\mathrm{NN}}=5$$\end{document} and 8 TeV in the kinematics of the various LHC experiments in a model independent way. It is of course justified since the most important ingredient in such evaluations is the probability of each kinematical configuration. Our computations for D mesons can also be extended to B meson production. To further illustrate the potentiality of the tool, we provide – for the first time – predictions for the nuclear modification factor for ηc\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta _c$$\end{document} production in p\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p$$\end{document}Pb collisions at the LHC.


Introduction
For many years, open and closed heavy-flavour production in hadron-hadron, hadron-nucleus and nucleus-nucleus collisions has been a major subject of investigations, on both experimental and theoretical sides (see [1] for a review in the context of the first LHC results and [2][3][4][5][6] for earlier reviews). In addition of helping us to understand the interface between the perturbative and non-perturbative regimes of QCD in hadron-hadron collisions, these reactions are also sensitive to -and thus probe -the properties of the possible deconfined state of matter (QGP) resulting from nucleus-nucleus (A A) collisions at ultra-relativistic energies. a e-mail: lansberg@in2p3.fr Yet, heavy-flavour production can also be affected by other nuclear effects 1 which are not related to a phase transition; they should in principle be subtracted in a way or another to study the QGP. These are typically believed to be the only ones acting in proton/deuteron-nucleus ( p A) collisions at fixed-target, RHIC and LHC energies. Experimental results from RHIC and the LHC in p A collisions [1] have shown that the yields and the spectra of J/ψ, ϒ, D and B are indeed modified in a magnitude which cannot simply be ignored in QGP studies. Many effects can be at play: break up within the nucleus [7,8] or with comovers for the quarkonia [9][10][11][12], coherent or incoherent energy loss [13][14][15][16][17], colour filtering [18], saturation/small-x/coherence effects [19][20][21][22][23], and the modification of the parton fluxes, as encoded in nuclear Parton Distribution Functions (nPDFs) [24][25][26][27][28].
In what follows, we will focus on the latter effects as a baseline for comparisons with experimental data. Our aim here is not to argue that it is indeed the dominant effect at RHIC and the LHC. Yet, a couple of recent comparisons [1,29,30] have shown that the magnitude of the gluon modification in usual nPDF fits is in reasonnable agreement with quarkonium, D and B meson data in pPb collisions at the LHC. nPDF fits are constantly updated with new data, 2 recently from the LHC, and we have found it useful to propose a simple and model-independent procedure to account for the nPDF impact, in the particular case of gluon-induced 2 → 2 1 In what follows, we will call them "cold nuclear matter effects". 2 We would like to mention that none of the global proton PDF fitting groups has so far used the heavy-flavour data in proton-proton collisions. However, some progress has recently been made in this direction. For example, the cross-section ratios of open-heavy-flavour production data from LHCb have been used to constrain the small-x gluon density in proton PDFs [31]. reactions. 3 Such a procedure, to be encoded in a user friendly forthcoming tool, would then allow anybody to make up one's mind about the typical expected magnitude of the gluon nuclear modifications on a given probe.
In the past, shortcut procedures using simplified kinematics (like the one of Drell-Yan at LO, that is 2 → 1) have widely been used [32][33][34][35][36]. However, it has been shown [37][38][39] that it can yield to systematic differences and, in principle, it cannot account for the P T dependence of the yield. In general, it is just better to rely on a more proper 2 → 2 kinematics, although some higher QCD corrections could involve more than 2 hard particles in the final state at large P T . For this purpose, a probabilistic Glauber Monte Carlo code, Jin [37][38][39], dedicated to the quarkonium case, has been developed to account for the geometry of the nuclear collisions and the impact parameter dependence of the nuclear effects at play along with the nPDF effect with an exact kinematics. However, as for now, the code deals with a limited number of processes (including though b production [40]) and of nPDFs; a simpler tool focusing on a 2 → 2 kinematics as the one we propose here is therefore very complementary. Eventually, both tools could interfaced or merged. Other tools, like APPLgrid [41] and fastNLO [42], with a similar spirit also exist but for different observables.
As will be explained below, the tool which we propose is based on HELAC-Onia [43,44] (but is not restricted to quarkonia) can use any nPDF set included in the library LHAPDF5 [45,46] and LHAPDF6 [47] and does not rely on any model for the hard-probe production, but on pp measurements which are used to tune the partonic-scattering elements, keeping in mind the limitation that it should be tuned for each process and in principle for possible different kinematical region with the chosen proton PDF and factorisation scale.

Our approach
As announced, our approach is based on a data-driven modelling of the scattering at the partonic level. Once folded with proton PDFs, they yield pp cross sections and, when folded with one proton PDF and one nuclear PDF, they yield p A cross sections. Such a choice is essentially motivated by the case of inclusive quarkonium hadroproduction. Firstly, it makes the computations faster with a limited loss of generality. Secondly, we have to acknowledge that we do not have at present time a global and consistent theoretical description of inclusive quarkonium production in the whole transverse momentum domain at hadron colliders. Thirdly, most of available models on the market show uncertainties larger than those of the data which they are meant to describe (and which sometimes they do not). Some of these observations also apply to D and B production.
This translates into the following advantages: 1. one can describe single quarkonium, D and B production in pp collisions in a very satisfactory way with only 2-3 tuned parameters for each meson; 2. the uncertainty within our approach is well controlled by the available pp data which, as just said, is much smaller than the theoretical uncertainties of the state-of-the-art calculations; 3. the method is much more efficient to generate events, with significantly reduced Monte Carlo uncertainty, owing to the simplicity of the computation.

pp cross section and partonic amplitude
As for the partonic scattering, we use a functional form for |A(k 1 k 2 → H + k 3 )| 2 initially proposed in [48] and then successfully used in [49][50][51][52] to model single quarkonium production at Tevatron and LHC energies in the context of double-parton scattering (DPS) studies. It reads where k i denote the partons involved in the hard scattering, x 1,2 are the momentum fractions carried by k 1,2 , s is the square of the centre-of-mass energy of the hadron collision, P T (M H ) is the transverse momentum (mass) of the produced particle, H, and θ(x) is the Heaviside step function. |A| 2 is meant to account for the squared amplitude averaged (summed) over the initial (final) helicity/colour factors. It contains 4 parameters λ, κ, P T , n, to be determined from the pp experimental data via a fit after the usual convolution with the PDFs: where f p denotes the proton PDF and 2 is the relativistic two-body phase space measure for the 2 → 2 scattering. The default factorisation scale at which are eveluated the nuclear and proton PDFs is taken as the transverse mass ( M 2 H + P 2 T ) of the particle H in what follows. In what follows, we will only consider processes which are dominated by gluon fusion at LHC energies. All the pro- cedure can readily be generalised to other partonic initial states.
2.2 Accounting for the nuclear PDF impact As announced, we will also only consider the nuclear modification of the PDF among the possible effects acting on quarkonia, D and B mesons. Such a restriction would probably not yield a good description of the quarkonium excited states [2], which we therefore do not discuss. Along the same lines, we will focus on the LHC regime where the nuclear absorption is likely negligible. It may not be so at RHIC and even less at fixed-target energies. Whereas one could think that the proposed procedure can be used to evaluate the sole impact of the nPDF on the excited states or the ground states at lower energies, one may want to be careful that in presence of other significant effects, the impact of the nPDFs may be affected. A clear example is a b-dependent anti-shadowing, which would tend to generate more J/ψ in the centre of the overlap zone, which then may have more chance to be broken up by the nuclear absorption than those produced in the periphery of the overlap zone. Yet, the procedure should give a right order of magnitude of the nPDF impact even if other effects are at play.
As it is customary, the yield of a particle H in p A collisions is obtained from that corresponding to the simple superposition of the equivalent number of pp collisions corrected by a factor encoding the nuclear modification of the parton flux. This is absolutely equivalent to directly using nuclear PDFs (normalised to the nucleus atomic number A) instead of proton PDFs. As aforementioned, our procedure does not currently rely on a Glauber code and we will thus restrict our studies to minimum bias collisions, i.e. integrated on all possible impact parameters b.
As such, the correction factor can be expressed in terms of the ratios R A i of the nuclear PDF (nPDF) in a nucleon belonging to a nucleus A to the PDF in the free nucleon: To illustrate the potentiality of our procedure, we will only use two of the most up-to-date nPDF parametrisations resulting from global analyses with uncertainties. The first is EPS09 [27], which provides the fit uncertainties at both leading order (LO, dubbed EPS09LO) and next-to-leading order (NLO, dubbed EPS09NLO) and is available in the library LHAPDF5 [45]. The nPDF effects are given in terms of R A i (x, Q 2 ) for all the flavours.
A new set, nCTEQ15 [24], has recently been released. It is available in the library LHAPDF6 [47] and provides NLO nuclear PDFs. As such, it is important to use the very same proton PDF as the one used for the fit. We have thus used CT14NLO [53]. In the case of EPS09, which provides ratios, the proton PDF to be used is less critical. In principle, we should have used CTEQ6(L1 or M) by consistency with EPS09, or CT14NLO for a good comparison of the yields with nCTEQ. Since CT14NLO is not available in LHAPDF5 and the code cannot load two PDF libraries at a time, we have preferred to use CT10NLO [54] which anyhow yields very similar gluon PDFs.

Fitting the LHC pp cross sections
At the LHC, we can essentially divide the inclusive (prompt) J/ψ production cross-sections measurements into 2 classes: the slightly forward and low P T (from 0 up to roughly 20 GeV) data of LHCb and ALICE and those from ATLAS and CMS at "high" P T (from 6-8 up to roughly 100 GeV). 4 We have performed 2 times 2 χ 2 fits of d 2 σ/d P T dy of prompt J/ψ production in pp collisions with 2 PDF sets (CT14NLO and CT10LO) using, on the one hand, LHCb data [55,56] and, on the other, ATLAS [57] and CMS [58] data. The fit parameters (λ, κ, P T and n of Eq. 1) are shown in Table 1. A comparison of our fit results with the experimental data is shown in Fig. 1a, d. The procedure is particularly successful,    (e)Υ [59] Inclusive Υ(1S) production at s=7 TeV LHC   Table 2 Results of a fit of d 2 σ/d P T dy of inclusive ϒ(1S) in pp collisions using CT10NLO and CT14NLO, where we have fixed the values of n and P T . The experimental data used in the fit are from ALICE [59], LHCb [61,64], ATLAS [65] and CMS [60]. (The uncertainties from the χ 2 fit below the per cent level are not shown.) CT10NLO 0.687 ± 0.367 0.0864 13.5 (fixed) 2 (fixed) Table 3 Results of a fit of d 2 σ/d P T dy of prompt η c (1S) in pp collisions using CT10NLO and CT14NLO, where we have fixed the values of n and P T . The experimental data used in the fit are from LHCb [62].
(The uncertainties from the χ 2 fit below the per cent level are not shown.) but for a few marginal bins. 5 These will nevertheless do not have a visible impact on the p A observables to be discussed later.
For the ϒ(1S) case, all the experiments have access to low P T data and there is no b feed-down contamination. We have performed 2 fits (with CT14NLO and CT10LO) using data from ALICE [59], LHCb [61,64], ATLAS [65] and CMS [60] altogether. See Table 2 for the fit results and Fig. 1e-h for comparison with the fit spectra.
For the prompt η c case, we have performed 2 fits (with CT14NLO and CT10LO) from the sole LHCb [62] data. See Table 3 for the fit results and Fig. 1h for comparison with the fit spectra.
As for the D 0 , we have performed 2 fits (with CT14NLO and CT10LO) from the LHCb [63] data. See Table 4 for the fit results and Fig. 1i for comparison with the fit spectra.
Just as for the J/ψ fits, the procedure works very well for ϒ(1S), η c and D 0 and gives us confidence that using the corresponding parametrised squared amplitudes will provide us with a reliable mapping of the x 1,2 , y and P T space. 5 Let us in particular note the slight discrepancy with the CMS very high P T data. The very same fits are however consistent with the ATLAS data in the same P T regime, which possibly indicates an underestimation of the systematical experimental uncertainties in that regime.   our results and the measurements of LHCb [66], ALICE [69] and ATLAS [68] 4 Results

4.
1 Rapidity and transverse-momentum dependence of the production cross-section in pPb collisions at √ s N N = 5.02 TeV Now that we have described our approach, we can present our results for the cross-section for quarkonium and D 0 production in proton-lead ( pPb) collisions at the LHC. In the following, we show comparisons with all the existing data. Our histograms are calculated under the same cuts as the experimental data. As announced, we have employed three different nPDF EPS09LO, EPS09NLO and nCTEQ15. In all the following plots, the uncertainty bands represent the nuclear PDF uncertainty only. In particular, we have not varied the factorisation scale despite the fact that it can indeed alter our results (see Sect. 4.4 for a short discussion).
The transverse-momentum P T spectra (dσ pPb /d P T ) of promptly produced J/ψ in pPb collisions at √ s N N = 5.02 TeV are shown in Fig. 2. Comparisons are made with the LHCb prompt J/ψ production data [66] in both the forward (1.  Fig. 2a. Figure 2b shows a comparison with the double differential cross sections d 2 σ pPb /d P T dy of J/ψ production of LHCb. Similarly, comparisons with the ALICE data [67] and ATLAS data [68] are given in Fig. 2c, d respectively. We note that ALICE data do not exclude the contribution from b-hadron decays. In general, the agreement with the yields differential in P J/ψ T is satisfactory both at low P T and high P T . In Fig. 3, we have compared the LHCb [66], ALICE [69] and ATLAS data [68] with the J/ψ cross-section differential in y. It is interesting to notice that the results with the three nPDF show different uncertainties. In the forward region (low x 2 ), the result with EPS09NLO has the smallest uncertainty and tend to overshoot the LHCb data [66] (see Fig. 3a). Such a discrepancy does not appear in Fig. 3b. One can also note that the EPS09LO uncertainty can be considered as the combination of both EPS09NLO and nCTEQ15 uncertainties in the forward region. In the backward region, owing to both the significant experimental and nPDF uncertainties, the three nPDFs are compatible with the data. At high P T (in Fig. 3c for the ATLAS data [68]), although the central values of the experimental data are systematically higher than our theoretical bands, they remain compatible within one standard deviation. There could indeed be an overestimation of the nPDF suppression in this region or an offset in the ATLAS data. ity dependence obtained with our procedure with the measurements by LHCb [70], ALICE [71] and ATLAS [72] and e the transversemomentum dependence as measured by ATLAS [72] As for the ϒ(1S), Fig. 4a, b show comparisons with the LHCb data. The agreement is better when the full LHCb range is considered as opposed to that when the LHCb acceptance is restricted to a range where equal positive and negative y can be accessed (Fig. 4b). A good agreement is also obtained with the ALICE data (Fig. 4c) in a similar rapidity domain. In the ATLAS acceptance, all three nPDF magnitudes correctly account for the yield differential in y and P T (Fig. 4d, e). Figures 5 and 6 show comparisons for the D 0 case between our results for the 3 nPDFs and the LHCb and ALICE measurements. The agreement is overall good. The yields tend to lie on the upper half of the uncertainty band. The nPDF uncertainties are however larger than for the quarkonia owing to the smaller value of the factorisation scale. As for the discrepancy in the first P T bin of Fig. 5, one should be careful that our pp parametrisation is not optimal to describe it as well (see Fig. 1i) and tend to undershoot the pp yield.
Finally, Fig. 7 show predictions -the first ever in the literature -for the P T and y differential yield of η c in the LHCb acceptance.

Rapidity and transverse-momentum dependence of
R pPb at √ s N N = 5.02 TeV We now present and discuss our results for the nuclear modification factor R pPb which characterises the yield modification of a given probe, say H, in pPb collisions relative to pp collisions. It is the ratio obtained by normalising the H yield in pPb collisions to the H yield in pp collisions in the same kinematical conditions (y, P T , nucleon-nucleon energy, etc.) times the average number of binary inelastic nucleon-nucleon collisions. When minimum bias collisions are considered, that is when all the possible geometrical configurations are summed over, it simplifies to the ratio of cross sections corrected by the atomic number of the nucleus (A = 208 for Pb): We first discuss the rapidity dependence of R pPb at the LHC with √ s N N = 5.02 TeV for J/ψ production. Our results obtained for the three nPDFs, EPS09LO, EPS09NLO and nCTEQ15 with their associated uncertainties are compared in Fig. 8 to the different experiments. Figure 8a, b show low P T data [66,67,69]. It is expected that the suppression in the forward region is due to the shadowing effect,  while the enhancement in the backward region is due to the anti-shadowing effect. The experimental data are compatible with these expectations. Among the three different nPDFs, the data tend to favour the result obtained with nCTEQ15. It is also interesting to note that the precision of the current data is already better than the nPDF uncertainties, especially in the forward region. This gives some hope that these measurements could ultimately be used to constrain the gluon density in heavy ions, provided that the impact of other nuclear effects could be disentangled. We also note that the shaded boxes on the right of the first two plots refer to the global systematical uncertainty. Such an information is not available for the ATLAS data. A good agreement with the LHCb and ALICE data is obtained; a slight discrepancy with the ATLAS data is observed. It is not clear whether it could be attributed to an offset in the data normalisation. In Fig. 9, we show further comparisons of R pPb vs P J/ψ T between our curves and the ALICE [67] and ATLAS [75] data. Similar to the rapidity distribution, a slight discrepancy is observed in Fig. 9b.
Similar comparisons are shown for ϒ(1S) on Figs. 9c and 10 . The overall agreement is acceptable given the large nPDF and experimental uncertainties. Further comparisons with the D 0 results are presented on Fig. 11. The agreement is also satisfactory and seems to indicate that EPS09 NLO is providing the best predictions. We however postpone further conclusions to the discussion of the R FB results which however do not necessarily confirm this observation. To complete this exhaustive list of comparisons, we present our predictions for R pPb of η c in the LHCb acceptance of its pp analysis on Fig. 12. We are hopeful that it will motivate the first ever experimental studies of η c in pPb collisions at the LHC.

Rapidity and transverse-momentum dependence of
In this section, we discuss the forward-to-backward production ratio R FB which results from the asymmetry of the proton-nucleus collision and is thus also sensitive to the nuclear effects. In addition, it has the advantage to be a ratio in which many of the systematic uncertainties of the data cancel, in particular that from the pp yield or cross section. It is defined as where the "forward" direction was defined as the flight direction of the proton beam. We stress that R FB is identically unity at y c.m.s. = 0. It tends to remain close to one if the nuclear effects cancel between the forward and backward regions, otherwise, it tends to increase more or less quickly for increasing |y c.m.s. |. We further note that in the current implementation of our code the nPDF uncertainties in R FB are generally smaller than in R pPb (or in the cross sections). Indeed, our current code uses the same nPDF eigenset to compute the forward and the backward yields used in a given ratio. This amounts to consider that the uncertainties in the R pPb are correlated. This interpretation (or rather use) of the information given by the nPDF is not unique and we could have considered that the nPDF uncertainties in R pPb in the forward and backward regions are not necessarily correlated (as the widespread use of theory "bands" may suggest). Doing so, the uncertainties in R FB would have been significantly larger. Finally, we recall   ALICE, the magnitude of the asymmetry is well compatible with that of nCTEQ15 and EPS09 LO, at the lower edge of the EPS09 NLO range. As for the ATLAS data with a P T cut, their current uncertainties and the reduced magnitude of the effects (since |y c.m.s. | is smaller) do not allow for any conclusions. Figure 14 shows our results for R FB versus P J/ψ T . A clear trend is seen in the LHCb and ALICE results with a ratio increasing with P T , starting at 0.6. R FB at P T above 10 GeV are compatible with unity, but the larger uncertainties do not exclude values smaller than one. The magnitude of the ratio in the data is compatible with the 3 nPDFs. More advanced studies are needed to go further in the interpretation of the P T dependence using specific eigensets as opposed to bands. We also recall that a given nPDF set can be compatible with R FB and not with R pPb . This can happen due to specific cancellations in the magnitude of the forward and backward nuclear modifications or to a normalisation offset. In particular, we note that there is no tension at all with the ATLAS data for R FB (Figs. 13c, 14c) unlike the case of R pPb (Figs. 8c, 9b). We are inclined to attribute this to a normalisation offset from the pp baseline whose effect disappears in R FB .
We have also computed R FB for ϒ(1S) (Fig. 15) and D 0 (Fig. 16). The same remarks as for the J/ψ case apply. No tension between the data and our computation are found. Just as for the ATLAS J/ψ data, the good agreement with the D 0 LHCb data may indicate that a slight offset in the normalisation affects R pPb as plotted on Fig. 11. Whereas the R pPb values point at a smaller suppression than those encoded in the nPDFs (in particular nCTEQ15), the magnitude of R FB is very well accounted by nCTEQ15 and corresponds to the strongest magnitude encoded in EPS09 NLO. For completeness, we have also computed R FB of prompt η c (1S) (see Fig.  17).

A few words on the (factorisation) scale dependence
As our results have shown, the nPDF uncertainties are significant and, in most cases, larger than those of the data. Moreover, the uncertainties of different nPDF sets do not necessarily overlap. Adding results obtained with the DSSZ set [26] would even probably enlarge the spread of the results. This is a strong motivation to learn to which extent these data could be used in the future to constrain the nPDFs. We are hopeful that our simple procedure can be an useful tool along this agenda.
Yet, as emphasised in [30], one should not forget one additional piece of uncertainty inherent to the use of nPDFs and the collinear factorisation, namely the (factorisation) scale, μ F , at which the parton densities should be evaluated. It is customary to pick up a natural scale related to the process (mass of the produced particle, typical momentum exchange, …) and to vary it by a factor two, or so, about this natural value. It happens that the magnitude of the nuclear effects encoded in the nPDFs do depend on μ F and we are not aware of strong theoretical arguments allowing one to bypass these considerations.
Under our procedure, there is absolutely no difficulty to evaluate the additional uncertainty attached to the unknown value of μ F . 7 Owing to our parametrisation of the pp data, the μ F dependence in the proton PDFs mostly cancels in the fit. If we were to choose a different value of μ F at which we evaluate the proton PDFs and refit the pp data, the denominator of R pPb would nearly be similar.
In Fig. 18, we thus only show the variation of R pPb due to the change of μ F in the evaluation of the nuclear PDFs, or in the ratio R g for EPS09. It is interesting -but not surprising -to note that both nPDF and μ F uncertainties depend 7 In fact, the evaluation of this uncertainty can be automatised as for that of the nPDFs. Comparison between the nPDF and the μ F relative uncertainties on R pPb as a function of rapidity (a) and transverse-momentum (b) for prompt J/ψ production in pPbcollisions at √ s N N = 5.02 TeV on y and P T . For y, it is obvious since it is connected to x 1 and x 2 . The μ F dependence is more marked at low P T which is also not a surprise. Overall, even though it is nearly everywhere smaller that the nPDF uncertainties (of a given set), the scale uncertainty is significant and should be kept in mind when comes the time for more quantitative theory-data comparisons.

Conclusions
We have devised a model-independent procedure to evaluate the impact of the nuclear modification of the gluon densities on hard probes produced in proton-nucleus collisions at colliders energies. It is particularly tailored for two-to-two par-tonic scatterings, relevant for quarkonium and heavy-meson production. The model independence of our procedure lies in the parametrisation of the partonic amplitude squared with parameters fit to pp collision data in similar kinematical conditions as the pPb data to be described; this is in contrast to other tools based on matrix elements calculated at a fixed order in perturbative theory. We have illustrated the capabilities of our approach by computing the cross sections as well as the nuclear modifications factor for J/ψ, ϒ and D 0 production at the LHC. Even though our objective was not to argue that the nPDF effect is the dominant one in this energy range, we have not found out any significant tension between our computations using three common nPDFs (EPS09 LO & NLO and nCTEQ15) and the existing data. To further highlight the potentialities of the approach, we have made predictions for η c production which might be at reach for the LHCb collaboration. We have also made predictions for the 8 TeV pPb run (see the Appendix).
As outlooks for physics studies, our method can easily be transposed to B hadron production. It should also be possible to apply it for non-prompt charmonia provided that the kinematical shift between the b-quark and the charmonium is correctly accounted for. On the side of the tool itself, 8 we plan to improve it such that it could automatically provide the user with the nuclear modification factors starting for measured pp data.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecomm ons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Funded by SCOAP 3 .

Appendix: Predictions for the pPb run at 8 TeV
In this appendix, we show predictions for the pPb run at 8 TeV using the foreseen experimental acceptance of the ALICE, ATLAS, CMS and LHCb experiments. These can serve as baselines to analyse whether nuclear effects beyond those encapsulated in the nPDF are visible at the highest energy possible for proton-nucleus at the LHC (Fig. 19).