Study of Drell-Yan dimuon production in proton-lead collisions at sNN\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}}} $$\end{document} = 8.16 TeV

Differential cross sections for the Drell-Yan process, including Z boson production, using the dimuon decay channel are measured in proton-lead (pPb) collisions at a nucleon-nucleon centre-of-mass energy of 8.16 TeV. A data sample recorded with the CMS detector at the LHC is used, corresponding to an integrated luminosity of 173 nb−1. The differential cross section as a function of the dimuon mass is measured in the range 15–600 GeV, for the first time in proton-nucleus collisions. It is also reported as a function of dimuon rapidity over the mass ranges 15–60 GeV and 60–120 GeV, and ratios for the p-going over the Pb-going beam directions are built. In both mass ranges, the differential cross sections as functions of the dimuon transverse momentum pT and of a geometric variable ϕ* are measured, where ϕ* highly correlates with pT but is determined with higher precision. In the Z mass region, the rapidity dependence of the data indicate a modification of the distribution of partons within a lead nucleus as compared to the proton case. The data are more precise than predictions based upon current models of parton distributions.


Introduction
The annihilation of a quark-antiquark pair into two oppositely charged leptons, through the exchange of a Z boson or a virtual photon (Z/γ * ) in the s-channel, is known as the Drell-Yan (DY) process [1]. The theoretical derivation of the matrix elements is available up to next-to-next-to-leading order in perturbative quantum chromodynamics (QCD) with next-to-leading order (NLO) electroweak (EW) corrections [2][3][4][5]. A precise measurement of this process can add valuable information on its nonperturbative part, including the effect of parton distribution functions (PDFs) [6].
Measurements of EW bosons in proton-nucleus and nucleus-nucleus collisions probe the nuclear modification of the PDFs [7][8][9][10]. The presence of a nuclear environment has been long observed [11] to modify the parton densities in the nucleus, as compared to those in a free nucleon. A first-principle description of such (nonperturbative) nuclear effects remains an open challenge, but they can be modelled using nuclear PDFs (nPDFs) determined with data in the same collinear factorisation approach as for free protons. Global fits of nPDFs [12][13][14][15][16][17][18][19] predict a suppression for small longitudinal momentum fraction x, x 10 −2 (i.e. shadowing [20] region), and an enhancement for intermediate x, 10 −2 x 10 −1 (i.e. antishadowing region).
In this paper, we report the measurement of the differential cross section for µ + µ − production via the DY process, as a function of the following variables: • dimuon mass, m µµ , in the interval 15 < m µµ < 600 GeV; • dimuon transverse momentum, p T , in two dimuon mass intervals  GeV and 60-120 GeV, targeting the continuum at low mass and the Z boson, respectively); • dimuon rapidity in the nucleon-nucleon centre-of-mass (CM) frame, y CM , in the same two mass intervals; and • φ * [36-38] (defined below) in the same two mass intervals.
The dimuon mass and φ * dependencies as well as cross sections in the dimuon mass range  GeV are reported for the first time in proton-nucleus collisions.
The variable φ * , used in numerous Z boson studies, is defined as φ * ≡ tan π − ∆φ 2 sin(θ * η ), (1.1) where ∆φ is the opening angle between the leptons, defined as the difference of their azimuthal angles in the plane transverse to the beam axis, and θ * η is related to the emission angle of the dilepton system with respect to the beam. The variable θ * η is defined in a frame that is Lorentz-boosted along the beam direction such that the two leptons are back-toback in the transverse plane. This angle θ * η is related to the pseudorapidities of the leptons by the relation cos(θ * η ) = tanh(∆η/2), (1.2) where ∆η is the difference in pseudorapidity between the two leptons. By construction, φ * is greater than zero. This quantity strongly correlates with the dimuon p T , while only depending on angular quantities for the leptons. Thus, it is measured with better precision than p T , especially at low p T values. Since φ * ∼ p T /m, where m is the mass of the dilepton system, the range φ * < 1 corresponds to dilepton p T up to about 100 GeV for a dilepton mass close to that of the Z boson. The outline of this paper is as follows. In section 2, the experimental methods are described, from the data and simulation samples used, up to the data analysis description and systematic uncertainties estimation. Results are presented and discussed in section 3, before the summary in section 4.

pythia.
The aforementioned event generators only simulate single proton-nucleon interactions, with the proportion of protons and neutrons found in Pb nuclei. To consider a more realistic distribution of the UE present in pPb collisions, simulated events are embedded into two separate samples of minimum bias (MB) events generated with epos LHC (v3400) [56], one for each pPb boost direction. The epos MC samples provide a good description of the global event properties of the MB pPb data, such as the η distributions of charged hadrons [57] and the transverse energy density [58].
A difference is found between the dimuon p T in powheg MC and that observed in data. To improve the modelling in the simulation, the powheg Z/γ * samples are reweighted event-by-event using an empirical function of the generated boson p T . This weight is applied in Z/γ * MC samples in the derivation of the various corrections described below. However, it is not applied in the figures of this paper, where the original p T spectrum from powheg is used.
The full detector response is simulated for all MC samples, using Geant4 [59], with alignment and calibration conditions tuned to match collision data, and a realistic description of the beam spot. The trigger decisions are also emulated, and the MC events are reconstructed with the standard CMS pp reconstruction algorithms used for the 2016 data.

JHEP05(2021)182
The Z/γ * , W, and tW samples are normalised to their NLO cross sections provided by powheg for pPb collisions, including EPPS16 modifications. The diboson samples are normalised to the cross sections measured by the CMS Collaboration in pp collisions at √ s . The small difference in CM energy with the pPb data is covered by the data-driven correction described in section 2.4 and smaller than the associated systematic uncertainty. The tt background is normalised to the CMS measurement in pPb collisions at √ s NN = 8.16 TeV [63]. All backgrounds receive a correction based on control samples in data, as described in section 2.4.
Simulated events do not feature the same event activity (charged-particle multiplicity or energy density) as the data, mostly because selecting two energetic muons favours higher-activity events (with a larger number of binary nucleon-nucleon collisions), while the epos sample used for embedding simulates MB events. To ensure a proper description of event activity in simulation, the distribution of the energy deposited in both sides of the HF calorimeter is reweighted event-by-event so that it matches that observed in data (selecting Z → µ + µ − events). The corresponding weights have a standard deviation of 0.27 for a mean of 1.

Object reconstruction and event selection
The events used in the analysis are selected with a single-muon trigger, requiring p T > 12 GeV for the muon reconstructed by the HLT. During both online and offline muon reconstruction, the data from the muon detectors are matched and fitted to data from the silicon tracker to form muon candidates. Each muon is required to be within the geometrical acceptance of the detector, |η lab | < 2.4. The leading muon (with highest p T ) is matched to the HLT trigger object and is required to have p T > 15 GeV, in the plateau of the trigger efficiency (around 95%, depending on η lab ). A looser selection of p T > 10 GeV is applied to the other muon.
Muons are selected by applying the standard "tight" selection criteria [46] used, e.g. in refs. [63, 64], with an efficiency of about 98%. Requirements on the impact parameter and the opening angle between the two muons are further imposed to reject cosmic ray muons. Events are selected for further analysis if they contain pairs of oppositely charged muons meeting the above requirements. The χ 2 divided by the number of degrees of freedom (dof) from a fit to the dimuon vertex must be smaller than 20, ensuring that the two muon tracks originate from a common vertex, thus reducing the contribution from heavy-flavour meson decays. In the rare events (about 0.4%) where more than one selected dimuon pair is found, the candidate with the smallest dimuon vertex χ 2 is kept.
To further suppress the background contributions due to muons originating from light and heavy flavour hadron decays, muons are required to be isolated, based on the p T sum of the charged-particle tracks around the muon. Isolation sums are evaluated in a circular region of the (η, φ) plane around the lepton candidate with ∆R < 0.3, where ∆R = (∆η) 2 + (∆φ) 2 . The relative isolation I rel , obtained by dividing this isolation sum by the muon p T , is required to be below 0.2.
In addition to the DY process, lepton pairs can also be produced through photon interactions. Exclusive coherent photon-induced dilepton production is enhanced in pPb -5 -

JHEP05(2021)182
collisions compared to pp data, because of the large charge of the lead nucleus. Hadronic collisions are selected by requiring at least one HF calorimeter tower with more than 3 GeV of total energy on either side of the interaction point. In order to further suppress the photon-induced background, characterised by almost back-to-back muons, events are required to contain at least one additional reconstructed track, which completely removes this background. Incoherent photon-induced dimuon production, where the photon is emitted from a parton instead of the whole nucleus and amounting to less than 5% of the total dilepton cross section according to studies in pp collisions at √ s = 13 TeV [29,65], is considered part of the signal and is neither removed nor subtracted.

Background estimation
Various backgrounds are estimated using one of the techniques described below, depending on the nature of the respective background process. Processes involving two isolated muons, such as Z/γ * → τ + τ − , tt, tW, and dibosons, are estimated from simulation and corrected using the "eµ method". Processes with one or more muons in jets, namely W+jets and multijet, are estimated using the "misidentification rate method".
The eµ method takes advantage of the fact that the EW backgrounds, as opposed to the Z/γ * → µ + µ − signal, also contribute to the eµ final state. Events with exactly one electron and one muon of opposite charge are used, where the muon is selected as described previously, matched to the HLT trigger muon and with p T > 15 GeV, while the electron [66] must have p T > 20 GeV and fulfil the same isolation requirement as the muon. The small contribution from heavy-flavour meson decays is estimated from same-sign eµ events. The data-to-simulation ratio with this selection, in each bin of the measured variables, is used to correct the simulated samples in the µ + µ − final state. This ratio is compatible with unity in most bins.
The misidentification rate method estimates the probability for a muon inside a jet and passing the tight selection criteria to pass the isolation requirements. This probability (the misidentification rate) is estimated as a function of p T , separately for |η lab | < 1.2 and |η lab | > 1.2. A sideband in data is selected from opposite-sign dimuon events in which the dimuon vertex χ 2 selection has been inverted. This sample is dominated by contributions from multijet and W+jets production, and the small contribution from EW processes, estimated using simulation, is removed. The misidentification rate is then applied to a control dimuon data sample, passing the dimuon vertex χ 2 selection but in which neither of the two muons passes the isolation requirement, to obtain the multijet contribution in the signal region, where both muons are isolated. The W+jets contribution is estimated with a similar procedure, using events in which exactly one of the two muons passes the isolation requirement. The small contribution from EW processes to these control data samples is estimated using simulation and removed. The multijet contribution in the sample with exactly one isolated muon is also accounted for, using the same technique. The validity of this method is checked in a control sample of same-sign dimuon data, which is also dominated by the multijet and W+jets processes. The same-sign data are found to be compatible with the predictions from the misidentification rate method in most bins, and the residual difference is accounted for as a systematic uncertainty.

JHEP05(2021)182
In figures 1 and 2, data are compared to the prediction from DY simulation and background expectations estimated using the techniques described above. A good overall agreement is found between the data and the expectation, which is dominated by the DY signal. Some hints for the differences will be discussed in terms of potential physics implications in section 3: they include data above expectation for m µµ < 50 GeV, as well as for y CM > 0 when 60 < m µµ < 120 GeV, and trends in dimuon p T and φ * , as mentioned in section 2.2.

Muon momentum scale and resolution corrections
The muon momentum scale and resolution are corrected in both data and simulation following the standard CMS procedure described in ref. [67]. These corrections have been derived using the pp data sample at √ s = 13 TeV recorded in 2016, with the same detector conditions as the pPb data set used in the present analysis.
In addition, the measurement is unfolded to account for finite momentum resolution. No regularisation is found to be needed given the good resolution and modest migrations between the analysis bins, and the maximum likelihood estimate [68] (obtained from the inversion of the response matrix, derived using simulated NLO powheg samples) is used to obtain the unfolded results. The effect of the unfolding is less than 1% in most cases, except for the mass dependence close to the Z boson mass peak, where it can amount to up to 15%.

Acceptance and efficiency
After subtraction of the contributions from different background processes, correction for the muon momentum resolution and scale, and unfolding for the detector resolution, the data need to be corrected for the acceptance and efficiency. The acceptance is defined as the fraction of generated signal events in the full phase space (within the quoted dimuon mass range and −2.87 < y CM < 1.93) passing the kinematic selection defining the so-called fiducial region: leading muon p T > 15 GeV, trailing muon p T > 10 GeV, and |η lab | < 2.4. Results are presented both with and without this acceptance correction, i.e. extrapolated to the full phase space and restricted to the fiducial region, respectively. The efficiency is the fraction of these events passing all other analysis selection criteria, including trigger selection, muon identification and isolation, and dimuon selection.
The efficiency is also checked in data, using Z boson events, with a tag-and-probe technique, as described in ref. [69]. The same procedure and corrections are used as in the measurement of W ± bosons in pPb collisions [64]. The observed differences between the efficiency in data and simulation, estimated separately for the trigger, identification, and isolation, are accounted for as scale factors on a per-muon basis that are applied to the simulated events. These corrections are applied both in the efficiency estimation and in the construction of the background templates described in section 2.4. When both muons in the event have p T > 15 GeV, they can both pass the single-muon trigger used in this data analysis, and the scale factor is computed from the product of inefficiencies. For the muon and central track reconstruction, the data and MC simulation are found to give comparable efficiencies (> 99.9%) and therefore no scale factor is applied for these two components of the efficiency.
Comparison of the data (black points) with the Z/γ * signal and background expectations (filled histograms, where "EW" includes Z/γ * → τ + τ − and diboson), estimated as described in the text, as a function of invariant mass (upper) and rapidity in the centre-of-mass frame for 15 < m µµ < 60 GeV (lower left) and 60 < m µµ < 120 GeV (lower right). Vertical error bars represent statistical uncertainties. The ratios of data over expectations are shown in the lower panels. The boson p T reweighting described in the text is not applied. The shaded regions show the quadratic sum of the systematic uncertainties (including the integrated luminosity, but excluding acceptance and unfolding uncertainties) and the nPDF uncertainties (CT14+EPPS16).    Comparison of the data (black points) with the Z/γ * signal and background expectations (filled histograms, where "EW" includes Z/γ * → τ + τ − and diboson), estimated as described in the text, as a function of p T (upper row) and φ * (lower row), for 15 < m µµ < 60 GeV (left) and 60 < m µµ < 120 GeV (right). The first bins of the p T and φ * distributions start at 0. Vertical error bars represent statistical uncertainties. The ratios of data over expectations are shown in the lower panels. The boson p T reweighting described in the text is not applied. The shaded regions show the quadratic sum of the systematic uncertainties (including the integrated luminosity, but excluding acceptance and unfolding uncertainties) and the nPDF uncertainties (CT14+EPPS16).

Final-state radiation effects
Muons may undergo final-state radiation before being measured in the CMS detector, biasing their momentum and shifting the dimuon mass to lower values. We unfold the measured distributions, after efficiency correction (as well as acceptance, if applicable), to the "pre-FSR" quantities, used for the presentation of our results and defined from a "dressed lepton" definition [28]. Generator-level muon four-momenta are recalculated by adding the four-momenta of all generated photons found inside a cone of radius ∆R = 0.1 around the muon. Again the response matrices for this unfolding procedure, derived using simulated NLO powheg samples, are found to be close to diagonal, thus no regularisation is needed in the unfolding.

Systematic uncertainties
Several sources of systematic uncertainties are evaluated. They are estimated in each bin of the measured distributions and added in quadrature. The list of systematic uncertainties is summarised in table 1 and details of the estimation of each source are given below.
Theoretical uncertainties have an impact on the acceptance and efficiency. The renormalisation and factorisation scales have been varied from half to twice their nominal value (set to the dimuon mass), and the envelope of the variations, excluding combinations where both scales are varied in opposite directions, is taken as an uncertainty. In addition, the strong coupling constant value is varied by 0.0015 from its default value, α S (m Z ) = 0.118, as recommended by PDF4LHC [70]. The CT14 and EPPS16 uncertainties are also included, estimated with LHAPDF6 [71] using the PDF4LHC recommendations for Hessian (n)PDF sets [70]. Finally, the full difference between the acceptance and efficiency obtained with and without the Z boson p T reweighting is considered as a systematic uncertainty. The impact of these uncertainties is less than 1% on the efficiency, but up to 10% on the acceptance for low dimuon masses.
We also include uncertainties stemming from the estimation of the efficiencies from data. The statistical component coming from the limited Z boson sample available is treated as a systematic uncertainty in this analysis. We also consider systematic effects associated with the choice of function used to model the p T behaviour of the efficiencies, the dimuon mass fitting procedure to the Z boson peak in the extraction of the efficiencies, a possible data-to-simulation difference in the muon reconstruction efficiency, and the effect of the mismodelling in simulation of the event activity and for additional interactions per bunch crossing. The magnitude of these uncertainties ranges from 1 to 5% at low dimuon mass.
Regarding the estimation of EW backgrounds with the eµ method, the statistical uncertainty in the correction factors is included as a systematic uncertainty, as well as the effect of varying the tt cross section by its uncertainty, 18% [63], the uncertainty in the transfer factor for the heavy-flavour contribution, and the difference between the data and simulation in the eµ distributions. The systematic uncertainty in the multijet and W+jets backgrounds, related to the misidentification rate method, receives several contributions. The statistical uncertainty in the templates derived from data is accounted for, and com--10 -
bined with the full difference between the nominal estimation and an alternative method (based on a different sideband in data, using same-sign dimuon events). The residual nonclosure in the same-sign data sample, as well as its statistical uncertainty, are also both added in quadrature to the other uncertainties related to the misidentification rate method. The total systematic uncertainty in the background estimate, dominated by the residual nonclosure in most bins, ranges from less than 0.5 to 15% (for large dimuon p T ). A different reweighting of the event activity in simulated samples is derived, as a function of the number of offline tracks reconstructed with |η lab | < 2.5 instead of the nominal correction using the total energy deposited in the HF calorimeters, which modifies the efficiency and the background estimation. The observed difference in the measurements, which is less than 1% in most bins, is taken as a systematic uncertainty.
Uncertainties in the muon momentum scale and resolution corrections have been evaluated, based on the 2016 pp data sample at √ s = 13 TeV, from which they are derived. These uncertainties, about 1% or less, arise from the limited data sample size available and variations in the method and its assumptions.
Response matrices used in the muon momentum scale and FSR unfoldings have been re-calculated using the first and second parts of the run alone (accounting for statistical uncertainties in simulation), and using the pyquen generator v1. 5.1 [72] instead of powheg (for a conservative estimation of the model dependence). Differences in the unfolded results, which are up to 2%, are taken into account as a systematic uncertainty.
Finally, the uncertainty in the integrated luminosity measurement is 3.5% [39].
Correlations across bins of these uncertainties have also been evaluated. Theoretical uncertainties are assumed to be fully correlated, with the exception of the nPDF uncer--11 -

JHEP05(2021)182
tainty, whose correlation is calculated using the CTEQ prescription for Hessian sets [73]. Systematic uncertainties in the efficiency scale factors obtained from control samples in data are assumed to be uncorrelated, since they could have different effects in different kinematic regions, while statistical correlations between the scale factors derived in the same region of the detector are accounted for. No correlation is assumed for the uncertainties related to the background estimation. Uncertainties related to the HF energy reweighting, unfolding, and integrated luminosity are treated as fully correlated between the bins and measurements, as well as each of the sources of uncertainty in the muon momentum scale and resolution corrections. The correlation matrices for systematic uncertainties are shown in figures 3 and 4, excluding the fully correlated integrated luminosity uncertainty for clarity. They are derived from the total covariance matrix, obtained from the sum of the covariance matrices for the individual sources, assuming the correlations above. For a given variable, the difference between the matrices in the two mass selections can be explained by the background uncertainty, which is one of the dominant systematic uncertainties for 15 < m µµ < 60 GeV but negligible most of the time for 60 < m µµ < 120 GeV, except at large p T or φ * . Muon efficiency uncertainties, treated as a function of |η lab |, induce a weak anticorrelation visible in systematic uncertainties as a function of rapidity, especially visible in the 60 < m µµ < 120 GeV region where they are the dominant systematic uncertainty.

Results and discussion
Fiducial cross section results, where the fiducial volume is defined from the single-muon p T and η lab selection, are shown in figures 5 and 6, as functions of the dressed lepton kinematic variables (as discussed in section 2.7), together with the expectations from powheg, using the CT14 [51] or CT14+EPPS16 [14] PDF sets. Cross sections in the full phase space, −2.87 < y CM < 1.93, i.e. including the acceptance correction for the single-muon kinematic selections, are presented in figures 7 and 8.
The CT14+EPPS16 predictions suffer from a larger uncertainty than CT14 alone, which is coming from the parametrisation of the nuclear modification of the PDFs. Since the dimuon rapidity is strongly correlated with the longitudinal momentum fraction x Pb of the parton in the lead nucleus, one can identify the shadowing region in the rapidity dependence of the cross section, in the full measured rapidity range for 15 < m µµ < 60 GeV and at positive rapidity for 60 < m µµ < 120 GeV. In the latter mass range, rapidities y CM −1 correspond to the antishadowing region. The inclusion of EPPS16 nuclear PDF modifications tends to provide a better description of the rapidity dependence in data for 60 < m µµ < 120 GeV than the use of the CT14 PDF alone. Uncertainties in the measurement are also smaller than nPDF uncertainties in the Z boson mass region for most analysis bins, showing that these data will impose strong constraints if included in future nPDF fits.
The mass dependence of the cross section sheds further light on the shadowing effects probed at low mass, i.e. at lower x Pb and lower scales than using Z bosons. The cross section measurement extends down to masses close to the Υ meson masses, with potential implications in the understanding of the interplay between nPDF and other effects in quarkonium production in proton-nucleus collisions [74]. The difference between the fiducial cross sections, shown in figures 5 and 6, and the ones corrected to the full phase space, shown in figures 7 and 8, is largest for low masses. The absence of acceptance correction in the former results reduces their model dependence and corresponding theoretical uncertainty, making clearer the trend for a higher cross section in data for low dimuon masses compared to the powheg expectation.
The p T and φ * dependencies of the cross section, especially in the Z boson mass region, both point to a slight mismodelling in powheg, reminiscent of the trend reported previously [35], where the data are softer than powheg predictions. The large sensitivity of these observables to the details of the QCD model, especially nonperturbative effects, is also observed in pp collisions [30] and prevents one from using them to draw unambiguous  Figure 5. Differential fiducial cross section (without the acceptance correction) for the DY process measured in the muon channel, as a function of the dimuon invariant mass (upper) and rapidity in the centre-of-mass frame for 15 < m µµ < 60 GeV (lower left) and 60 < m µµ < 120 GeV (lower right). The error bars on the data represent the quadratic sum of the statistical and systematic uncertainties. Theory predictions from the powheg NLO generator are also shown, using CT14 (blue) or CT14+EPPS16 (red). The boxes show the 68% confidence level (n)PDF uncertainty on these predictions. The ratios of predictions over data are shown in the lower panels, where the data and (n)PDF uncertainties are shown separately, as error bars around one and as coloured boxes, respectively.  Figure 6. Differential fiducial cross sections (without the acceptance correction) for the DY process measured in the muon channel, as functions of p T (upper row) and φ * (lower row), for 15 < m µµ < 60 GeV (left) and 60 < m µµ < 120 GeV (right). The first bin of the p T and φ * measurements starts at 0. The error bars on the data represent the quadratic sum of the statistical and systematic uncertainties. Theory predictions from the powheg NLO generator are also shown, using CT14 (blue) or CT14+EPPS16 (red). The boxes show the 68% confidence level (n)PDF uncertainty on these predictions. The ratios of predictions over data are shown in the lower panels, where the data and (n)PDF uncertainties are shown separately, as error bars around one and as coloured boxes, respectively.  Figure 7. Differential cross section for the DY process measured in the muon channel, as a function of the dimuon invariant mass (upper) and rapidity in the centre-of-mass frame for 15 < m µµ < 60 GeV (lower left) and 60 < m µµ < 120 GeV (lower right). The error bars on the data represent the quadratic sum of the statistical and systematic uncertainties. Theory predictions from the powheg NLO generator are also shown, using CT14 (blue) or CT14+EPPS16 (red). The boxes show the 68% confidence level (n)PDF uncertainty on these predictions. The ratios of predictions over data are shown in the lower panels, where the data and (n)PDF uncertainties are shown separately, as error bars around one and as coloured boxes, respectively.  Figure 8. Differential cross sections for the DY process measured in the muon channel, as functions of p T (upper row) and φ * (lower row), for 15 < m µµ < 60 GeV (left) and 60 < m µµ < 120 GeV (right). The first bin of the p T and φ * measurements starts at 0. The error bars on the data represent the quadratic sum of the statistical and systematic uncertainties. Theory predictions from the powheg NLO generator are also shown, using CT14 (blue) or CT14+EPPS16 (red). The boxes show the 68% confidence level (n)PDF uncertainty on these predictions. The ratios of predictions over data are shown in the lower panels, where the data and (n)PDF uncertainties are shown separately, as error bars around one and as coloured boxes, respectively.  Table 2. χ 2 values between the data and the powheg predictions and associated probability, from the fiducial cross sections, when experimental and theoretical bin-to-bin correlations are taken into account. The integrated luminosity uncertainty is included in the experimental uncertainties.
In tables 2 and 3, the χ 2 values between the data and the predictions are reported, accounting for the bin-to-bin correlations for experimental (systematic uncertainties, shown in figures 3 and 4) and theoretical (from nPDF) uncertainties. The observations discussed above from figures 5 to 8 can be made here more quantitatively and more precisely with fiducial cross sections, thanks to the smaller systematic uncertainty. The inclusion of the EPPS16 modifications to the PDFs of the lead nucleus tends to improve the description for y CM in the Z boson mass region, but conclusions are not clear for other quantities, and could even be opposite in the case of p T and φ * in that region. However, the manifestly imperfect modelling of the cross sections in powheg prevents from drawing strong conclusions about nPDFs using these variables.
Forward-backward ratios (R FB ) are built from the rapidity-dependent cross sections in the two mass regions, defined as the ratio of the y CM > 0 to the y CM < 0 cross sections (p-going to Pb-going). They are shown in figure 9. In both mass regions, the R FB is by construction equal to unity in the absence of nuclear effect (CT14), but decreasing with |y CM | with CT14+EPPS16 and CT14+nCTEQ15WZ [19]. Similar conclusions are drawn as from the rapidity dependence of the cross section, but the construction of these ratios allows for the partial cancellation of theoretical and experimental uncertainties,  Table 3. χ 2 values between the data and the powheg predictions and associated probability, from the full phase space cross sections, when experimental and theoretical bin-to-bin correlations are taken into account. The integrated luminosity uncertainty is included in the experimental uncertainties.
accounting for the correlations described in the previous section. In particular, for 60 < m µµ < 120 GeV and at large |y CM |, an indication of a forward-backward ratio smaller than unity is found, consistent with the expectation from the combination of shadowing and antishadowing effects expected with CT14+EPPS16, as well as with similar results from W bosons [64]. Predictions using CT14+nCTEQ15WZ are found to be in good agreement with the data. The larger amount of shadowing in nCTEQ15 [15], hinted by the recent W boson measurement [64], is not predicted with nCTEQ15WZ. The low mass region is less conclusive, but nPDF uncertainties are smaller in this selection for nCTEQ15WZ than for EPPS16. Finally, experimental uncertainties for 60 < m µµ < 120 GeV are smaller than the nPDF ones, once again showing relevance of these data to the study of nPDF effects.

Summary
Differential cross section measurements of the Drell-Yan process in the dimuon channel in proton-lead collisions at √ s NN = 8.16 TeV have been reported, including the transverse momentum (p T ) and rapidity dependencies in the Z boson mass region (60 < m µµ < 120 GeV). In addition, for the first time in collisions including nuclei, the p T and rapidity dependence for smaller masses 15 < m µµ < 60 GeV have been measured. The dependence with φ * (a geometrical variable that highly correlates with dimuon p T but is determined with higher precision) for both 15 < m µµ < 60 GeV and 60 < m µµ < 120 GeV and the mass dependence from 15 to 600 GeV have been presented, also for the first time in proton-nucleus collisions. Finally, forward-backward ratios have been built from the rapidity-dependent cross sections for y CM > 0 to y CM < 0 in both mass regions, highlighting the presence of nuclear effects in the parton distribution functions.
Results for 60 < m µµ < 120 GeV are the most precise to date, featuring smaller uncertainties than the theoretical predictions, and provide novel constraints on the quark and antiquark nuclear parton distribution functions (nPDFs). Measurements in the lower mass . Forward-backward ratios for 15 < m µµ < 60 GeV (left) and 60 < m µµ < 120 GeV (right). The error bars on the data points represent the quadratic sum of the statistical and systematic uncertainties. The theory predictions from the powheg NLO generator are also shown, using CT14 [51] (blue), CT14+EPPS16 [14] (red), or CT14+nCTEQ15WZ [19] (green) PDF sets. The boxes show the 68% confidence level (n)PDF uncertainty in these predictions. The ratios of predictions over data are shown in the lower panels, where the data and the (n)PDF uncertainties are shown separately, as error bars around one and as coloured boxes, respectively. range 15 < m µµ < 60 GeV give access to a new phase space for nPDF studies, extending to lower longitudinal momentum fraction x and lower energy scale Q 2 . The p T -and φ *dependent results are also very sensitive to the details of model details, such as soft quantum chromodynamics phenomena, which they may help to better understand in pPb collisions.

Acknowledgments
We congratulate our colleagues in the CERN accelerator departments for the excellent performance of the LHC and thank the technical and administrative staffs at CERN and at other CMS institutes for their contributions to the success of the CMS effort. In addition, we gratefully acknowledge the computing centres and personnel of the Worldwide LHC Computing Grid and other centres for delivering so effectively the computing infrastructure essential to our analyses. Finally, we acknowledge the enduring support for the construction and operation of the LHC, the CMS detector, and the supporting computing infrastructure provided by the following funding agencies: BMBWF and FWF (Austria);  -26 -