Measurement of the mass dependence of the transverse momentum of lepton pairs in Drell–Yan production in proton–proton collisions at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{s} = 13\,\text {Te\hspace{-.08em}V} $$\end{document}s=13TeV

The double differential cross sections of the Drell–Yan lepton pair (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\ell ^+\ell ^-$$\end{document}ℓ+ℓ-, dielectron or dimuon) production are measured as functions of the invariant mass \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m_{\ell \ell }$$\end{document}mℓℓ, transverse momentum \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_{\textrm{T}} (\ell \ell )$$\end{document}pT(ℓℓ), and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varphi ^{*}_{\eta }$$\end{document}φη∗. The \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varphi ^{*}_{\eta }$$\end{document}φη∗ observable, derived from angular measurements of the leptons and highly correlated with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_{\textrm{T}} (\ell \ell )$$\end{document}pT(ℓℓ), is used to probe the low-\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_{\textrm{T}} (\ell \ell )$$\end{document}pT(ℓℓ) region in a complementary way. Dilepton masses up to 1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\,\text {Te\hspace{-.08em}V}$$\end{document}TeV are investigated. Additionally, a measurement is performed requiring at least one jet in the final state. To benefit from partial cancellation of the systematic uncertainty, the ratios of the differential cross sections for various \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m_{\ell \ell }$$\end{document}mℓℓ ranges to those in the Z mass peak interval are presented. The collected data correspond to an integrated luminosity of 36.3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\,\text {fb}^{-1}$$\end{document}fb-1 of proton–proton collisions recorded with the CMS detector at the LHC at a centre-of-mass energy of 13\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\,\text {Te\hspace{-.08em}V}$$\end{document}TeV. Measurements are compared with predictions based on perturbative quantum chromodynamics, including soft-gluon resummation.


Introduction
The Drell-Yan (DY) production of charged-lepton pairs in hadronic collisions [1] provides important insights into the partonic structure of hadrons and the evolution of the parton distribution functions (PDFs). At leading order (LO) in perturbative quantum chromodynamics (pQCD), the DY process is described in terms of an s-channel Z/γ * exchange process convolved with collinear quark and antiquark parton distribution functions of the proton. At LO, the lepton pair transverse momentum p T ( ), corresponding to the exchanged boson transverse momentum, is equal to zero. At higher orders, initial-state QCD radiation gives rise to a sizable p T ( ). Whereas the spectrum for large p T ( ) values is expected to be described through fixed-order calculations in pQCD, at small values (p T < O(m )), where m is the invariant mass of the lepton pair, soft-gluon resummation to all orders is required [2,3]. In addition, the low-p T ( ) region also includes the effects of the intrinsic transverse motion of the partons in the colliding hadrons that has to be extracted from data and parameterized. The resummation functions are universal and obey renormalisation group equations, predicting a simple scale dependence in the leading logarithmic approximation, where the scale is given by m . Therefore, measuring the p T ( ) spectrum in a wide m range tests the validity of the resummation approach and the precision of different predictions. Calculations for inclusive DY production as a function of m and p T ( ) are available up to next-to-next-to-leading order (NNLO) in pQCD [4][5][6][7]. Soft-gluon resummation can be computed analytically, either in transverse momentum dependent parton distributions (TMD) or in parton showers of Monte Carlo (MC) event generators matched with matrix element calculations [8][9][10][11][12][13][14][15].
The p T ( ) resolution is dominated by the uncertainties in the magnitude of the transverse momenta of the leptons, whereas the measurement precision of the lepton angle does not contribute significantly. The kinematic quantity ϕ * η [16][17][18], derived from these lepton angles, is defined by the equation: The variable ∆ϕ is the opening angle between the leptons in the plane transverse to the beam axis. The variable θ * η is the scattering angle of the dileptons with respect to the beam in the longitudinally boosted frame where the leptons are back to back. It is related to the pseudorapidities of the oppositely charged leptons by the relation cos(θ * η ) = tanh [(η − − η + )/2]. The variable ϕ * η , by construction greater than zero, is closely related to the normalized transverse momentum p T ( )/m [16]. Since ϕ * η depends only on angular variables, its resolution is significantly better than that of the transverse momentum, especially at low-p T ( ) values, but its interpretation in terms of initial-state radiation (ISR) is not as direct as that of p T ( ).
The DY process in the presence of one jet is a complementary way to investigate the initial-state QCD radiations. The requirement of a minimal transverse momentum associated with this jet is reflected in the p T ( ) distribution by momentum conservation. When more hadronic activity than a single jet is present in the events, the transverse momentum balance between the leading jet and the lepton pair has a broad distribution. As a consequence, the full p T ( ) spectrum in the presence of jets brings additional information, since at small values it is sensitive to numerous hard QCD radiations. Furthermore, DY production in association with at least one jet also brings up contributions where virtual partons acquire transverse momentum, whose collinear radiations will have a significant angle with respect to the beam, which contributes as a component of the final p T ( ). This paper presents a DY differential cross section measurement in bins of m , over the range of 50 GeV to 1 TeV, as functions of p T ( ) and ϕ * η for inclusive DY production, and in events with at least one jet as a function of p T ( ). The data were collected in 2016 with the CMS detector at the CERN LHC, corresponding to an integrated luminosity of 36.3 fb −1 of protonproton (pp) collisions at a centre-of-mass energy of √ s = 13 TeV. To reduce the uncertainties, the measured cross sections combine measurements of separately extracted cross sections for the electron and the muon channels. The measurements presented in this paper are extensively discussed in Ref. [19].
Complementary measurements of the DY process have been performed recently by the CMS, ATLAS, and LHCb Collaborations at the CERN LHC [20-38] and by the CDF and D0 Collaborations at the Fermilab Tevatron [39][40][41][42][43][44][45]. The cross section measurements presented in this paper extend the mass range below and above the Z boson resonance with respect to the previous CMS measurements of p T ( ) dependence.
The outline of this paper is the following: in Section 2 a brief description of the CMS detector is given. In Section 3 the selection criteria of the measurement are described. The simulation samples used in the measurement are described in Section 4. Section 5 explains the details of the unfolding procedure and the systematic uncertainties are given in Section 6. Theory predictions used for comparison with the measurements are described in Section 7. The results are presented in Section 8 and a summary of the paper is given in Section 9.

The CMS detector
The central feature of the CMS apparatus is a superconducting solenoid of 6 m internal diameter, providing a magnetic field of 3.8 T. Within the solenoid volume are a silicon pixel and strip tracker, a lead tungstate crystal electromagnetic calorimeter (ECAL), and a brass and scintillator hadron calorimeter (HCAL), each composed of a barrel and two endcap sections. Forward calorimeters extend the η coverage provided by the barrel and endcap detectors. Muons are detected in gas-ionization chambers made of detection planes using three technologies: drift tubes, cathode strip chambers, and resistive plate chambers, embedded in the steel flux-return yoke outside the solenoid.
The global event reconstruction (also called particle-flow event reconstruction [46]) reconstructs and identifies each individual particle in an event, with an optimized combination of all subdetector information. In this process, the identification of the particle type (photon, electron, muon, charged or neutral hadron) plays an important role in the determination of the particle direction and energy.
Electrons are identified as a primary charged particle track and potentially many ECAL energy clusters corresponding to this track extrapolation to the ECAL and to possible bremsstrahlung photons emitted along the way. The electron momenta are estimated by combining energy measurements in the ECAL with momentum measurements in the tracker [47]. The momentum resolution for electrons with p T ≈ 45 GeV from Z → ee decays ranges from 1.7 to 4.5%. It is better in the barrel region than in the endcaps, and also depends on the bremsstrahlung energy emitted by the electron as it traverses the material in front of the ECAL.
Muons are identified as tracks in the central tracker consistent with either a track or several hits in the muon system, and associated with calorimeter deposits compatible with the muon hypothesis. The reconstructed muon global track, for muons with 20 < p T < 100 GeV, has a relative transverse momentum resolution of 1.3-2.0% in the barrel and better than 6% in the endcaps. The p T resolution in the barrel is better than 10% for muons with p T up to 1 TeV [48]. The resolution is further improved with corrections derived from the Z mass distribution [49].
Charged hadrons are identified as charged particle tracks not identified as electrons or as muons. Finally, neutral hadrons are identified as HCAL energy clusters not linked to any charged-hadron trajectory, or as a combined ECAL and HCAL energy excess with respect to the expected charged-hadron energy deposit. For each event, hadronic jets are clustered from these reconstructed particles using the infrared-and collinear-safe anti-k T algorithm [50,51] with a distance parameter of 0.4. Jet momentum is determined as the vectorial sum of all particle momenta in the jet, typically within 5 to 10% of the true momentum over the entire p T spectrum and detector acceptance.
The primary vertex (PV) is taken to be the vertex corresponding to the hardest scattering in the event, evaluated using tracking information alone, as described in Section 9.4.1 of Ref. [52].
Events of interest are selected using a two-tiered trigger system. The first level (L1), composed of custom hardware processors, uses information from the calorimeters and muon detectors to select events at a rate of around 100 kHz within a fixed latency of about 4 µs [53]. The second level, known as the high-level trigger, consists of a farm of processors running a version of the full event reconstruction software optimized for fast processing, and reduces the event rate to around 1 kHz before data storage [54].
A more detailed description of the CMS detector is reported in Ref. [55], together with a definition of the coordinate system used and the relevant kinematic variables.

Event selection
The initial event selection requires a dielectron trigger with a p T threshold of 23 and 12 GeV on the two leading electrons in the electron channel. In the muon channel we require a dimuon trigger with p T thresholds of 18 and 7 GeV or a single-muon trigger with a p T threshold of 24 GeV. The final selection is restricted to the region where the triggers are fully efficient: p T > 25 GeV for the leading lepton, p T > 20 GeV for the subleading lepton and |η| < 2.4 for both channels.
An event must contain exactly two isolated leptons of the same flavour (with the isolation criteria as detailed in Ref. [26]). In addition the two leptons must have opposite charges. Events with a third lepton with p T greater than 10 GeV and |η| < 2.4 are vetoed.
Due to the high instantaneous luminosity of the LHC, additional proton-proton interactions occur during the same bunch crossing (pileup) that contribute additional overlapping tracks and energy deposits in the event, and result in an apparent increase of jet momenta. To mitigate this effect, tracks identified as originating from pileup vertices are discarded and an offset correction is applied to correct for the remaining neutral pileup contributions [56]. The two identified leptons can be reconstructed as jets. Those jets are disregarded by requiring a separation, ∆R = √ (∆η) 2 + (∆φ) 2 , between the reconstructed jets and these lepton candidates to be larger than 0.4.
To suppress the contamination of jets coming from pileup, a multivariate discriminant is used. The pileup contamination is also reduced by the choice of the final selection: jets are required to have a minimum transverse momentum of 30 GeV and, to ensure high-quality track information, they are limited to a rapidity range of |y| < 2.4.
To reduce the tt background, events containing one or more b tagged jets are vetoed. The medium discrimination working point of the combined secondary vertex b tagging algorithm [57] is used. The effect on the signal is small and is corrected for in the unfolding procedure.
The effects of finite detector resolution and selection efficiency are corrected by using the unfolding procedure described in Section 5. Scale factors are applied to the simulation used for the unfolding, to correct for differences with respect to the data in the efficiencies of the different selections: trigger, lepton identification, lepton isolation, and b-tagged jet veto. For the trigger, the factor is given as a function of |η| of the two leptons and is applied once per lepton pair. The value of the scale factor is close to one. When dealing with the identification and isolation efficiencies, the scale factor is given per lepton as a function of its p T and |η|, and applied to each of the two selected leptons [26].

Simulated samples and backgrounds
For the simulation of the Z/γ * process (including the τ + τ − background), a sample is generated with MADGRAPH5 aMC@NLO [58] version 5.2.2.2 (shortened to MG5 aMC) using the FxFx jet merging scheme [59]. The parton shower, hadronization, and QED final-state radiation (FSR) are calculated with PYTHIA 8.212 [60] using the CUETP8M1 tune [61]. The matrix element calculations include Z/γ * + 0, 1, 2 jets at next-to-leading order (NLO), giving an LO accuracy for Z/γ * + 3 jets. The NLO NNPDF 3.0 [62] is used for the matrix element calculation. In control plots and when comparing to the measurement, this prediction is normalized to the cross section obtained directly from the generator, 1977 pb per lepton channel (for m > 50 GeV).
Other processes that can give a final state with two oppositely charged same-flavour leptons are WW, WZ, ZZ, γγ, tt pairs, and single top quark production. The tt and single top backgrounds are generated at NLO using the POWHEG version 2 [63][64][65][66] interfaced to PYTHIA 8. Background samples corresponding to diboson electroweak production (denoted VV in the figure legends) [67] are generated at NLO with POWHEG interfaced to PYTHIA 8 (WW) or at LO with PYTHIA 8 alone (WZ and ZZ). These samples are generated using NLO NNPDF 3.0 for the matrix element calculation. The γγ background process leading to two charged leptons in the final state, γγ → + − , is simulated using LPAIR [68,69] interfaced with PYTHIA 6 and using the default γ-PDF of Suri-Yennie [70]. This contribution is split into three components, since the interaction at each proton vertex process can be elastic or inelastic.
The total cross sections of WZ and ZZ diboson samples are normalized to the NLO prediction calculated with MCFM v6.6 [71], whereas the cross sections of the WW samples are normalized to the NNLO prediction [72]. The total cross section of the tt production is normalized to the prediction with NNLO accuracy in QCD and next-to-next-to-leading logarithmic (NNLL) accuracy in soft gluon resummation calculated with TOP++ 2.0 [73]. The single top and γγ background distributions are normalized to the cross sections calculated by their respective event generators.
It is possible for hadrons to mimic the signature of an electron in the detector. The main processes that contribute to this background are W + jet production, when the W decays leptonically, and QCD multijet events. Such backgrounds are nonnegligible only in the electron channel.
The contamination of the signal region by events containing hadrons misidentified as electrons is estimated using a control region where two electrons of the same sign are required. This control region mainly contains events with hadrons misidentified as electrons and events originating from the DY process when the charge of one electron is incorrectly attributed. The probability of charge misidentification is obtained as a function of p T and η of each electron in the Z peak region (81 < m < 101 GeV), where the hadron contamination is negligible even in the control region. These probabilities are then used to estimate the charge misidentification rate for other values of m . The difference between the observed number of events in the control region and the estimated charge misidentification rate is assumed to be the contamination from hadron background. We observe that the numbers of misidentified-lepton events in the same-sign electron sample and in the signal (opposite-sign electron) sample are compatible.
The number of events at the reconstructed level is compared with the sum of the contributions from signal and backgrounds. In Fig. 1, the dilepton mass spectrum is shown for both the electron and and the muon channels, whereas Fig. 2 shows the p T ( ) distributions in various m bins for the electron channel only. Globally, the background contamination is lower than 1%. The background becomes around 10% for m outside of the Z boson mass peak and up to 30% in some bins. The simulated samples are processed through a GEANT4 [74] based simulation of the CMS detector, with the same reconstruction algorithms as of data. They also include a pileup profile that is reweighed to match the profile of the data.    Figure 1: Distributions of events passing the selection requirements in the muon (left) and electron channels (right). Each plot also presents in the lower part a ratio of simulation over data. Only statistical uncertainties are shown as error bars on the data points, whereas the ratio presents the statistical uncertainty in the simulation and the data. The plots show the number of events without normalization to the bin width. The different background contributions are discussed in the text.

Measured observables and unfolding procedure
The measurement of the DY cross section is carried out with respect to the p T and ϕ * η of the dilepton pairs produced inclusively, and with respect to p T for pairs produced in association with at least one jet. For the inclusive case, the measurement is divided into five invariant mass bins: 50-76, 76-106, 106-170, 170-350, 350-1000 GeV; the last bin is not included when requiring at least one jet because of the small number of events available. The measurement of the ratio of cross section in mass bins 50-76, 106-170, 170-350, 350-1000 GeV to the cross  Data    Fig. 1. section around the Z mass peak(76-106 GeV) is also performed. The bin widths are chosen to be as small as possible, based on the detector resolution and the number of events.
To correct for the detector resolution and the efficiency of the selection, an unfolding procedure is applied to the measured distributions one dimensionally in each mass bin. To obtain the particle-level distributions from the reconstructed distributions, the unfolding uses a response matrix based on the simulated signal sample. To unfold, the D'Agostini iterative method with early stopping is used as implemented in ROOUNFOLD [75]. The result, converging towards the maximum likelihood estimate, is affected by fluctuations increasing with the number of iterations. The fluctuations are studied using pseudo-experiments for each number of iterations following the method used in Ref. [76]. The procedure is stopped just before the fluctuations become significant with respect to the statistical uncertainty. The number of iterations ranges between 4 and 25.
The particle level refers to stable particles (cτ > 1 cm), other than neutrinos, in the final state. To correct for energy losses due to QED FSR, leptons are "dressed", i.e., all the prompt photons with a distance smaller than ∆R = 0.1 to the lepton axis are added to the lepton momentum. The cross section is extracted in the following phase space: leading and subleading dressed leptons satisfying p T > 25 and 20 GeV and |η| < 2.4. When at least one additional jet is required, it must satisfy p T > 30 GeV, |y| < 2.4, and be spatially separated from the dressed leptons by ∆R > 0.4.
The cross sections are first extracted separately for the electron and muon channels. They are compatible for all studied distributions and the two channels are combined to reduce the statistical uncertainties. The combined differential cross sections are calculated bin-by-bin as the weighted mean values of the differential cross sections of the two channels. The systematic and statistical uncertainties are obtained using the linear combination method described in Ref. [77], considering as fully correlated the uncertainties in the jet energy scale and resolution, the pileup, the background subtraction, b tagging, and the integrated luminosity. Other uncertainties are considered as uncorrelated.

Systematic uncertainties
Several sources of uncertainties in the measurement are considered. The integrated luminosity is measured with a precision of 1.2% [78], which results in a relative uncertainty of almost the same value in the measurement. Small variations are caused by the subtraction of the background contributions estimated from the simulation.
The uncertainties coming from the lepton trigger efficiencies are estimated by varying the applied scale factors up and down by one standard deviation. The uncertainties from identification and reconstruction efficiencies are estimated for various sources including QED FSR, resolution, background modeling, and the tag object selection in the tag-and-probe procedure, as well as the statistical component treated separately for each scale factor in p T and η of the lepton [26]. The efficiency uncertainties include a one percent effect in the L1 trigger caused by a timing problem in ECAL endcaps. The lepton energy scale uncertainties are estimated by varying the lepton energy and p T by ±1 standard deviation (reach 0.75% (0.5%) for electrons (muons) depending on η and p T ). Uncertainties coming from the lepton energy resolution are estimated by spreading the lepton energy using the generator-level information.
The uncertainty in the jet energy scale is estimated by varying the jet momenta in data by 2.5% to 5%, depending on the energy and pseudorapidity of the jet. The uncertainty in the jet energy resolution is estimated by varying the smearing factor used to match the simulated jet energy resolution to data by ±1 standard deviation around its central value.
A systematic uncertainty is attributed to the normalisation of the background samples estimated by Monte Carlo event generators. The theoretical uncertainty in the cross section of the dominant tt background is≈6%, using the TOP++ 2.0 program and including scale and PDF variations. The uncertainties in the other background cross sections are smaller. In particular, it has been verified that 6% covers the differences of the γγ → + − samples generated using Suri-Yennie and LuxQED [79,80] photon PDFs. In a conservative way, the uncertainties in all other Monte Carlo based background estimates are also estimated to be 6%. This uncertainty is applied to fully-elastic, semi-and fully-inelastic cases.
The uncertainty in the misidentified electron background estimation using same-sign events is obtained using an uncertainty in the charge misidentification estimation of about 10% per electron at p T (e) = 150 GeV, rising with p T (e). A 20% total uncertainty in the charge misidentification is used and propagated to the estimate of this background.
Alternative pileup profiles are generated by varying the amount of pileup events by 5%, and the difference to the nominal sample is propagated to the final results.
The unfolding model uncertainty is estimated by reweighting the simulated sample to match the data shape for each distribution, and using this as an alternate model for unfolding. The difference with respect to the results obtained with the simulated sample is assigned as the uncertainty. The statistical uncertainty coming from the limited sample size is also included, provided by the ROOUNFOLD package.
The systematic uncertainties are propagated to the measurement through the unfolding procedure by computing new response matrices varying the quantities by one standard deviation up and down. All the experimental uncertainties are symmetrized by taking the average of the deviations from the central value. The uncertainty sources are independent and the resulting uncertainties are added in quadrature.
For the inclusive measurement the main sources of uncertainties are the integrated luminosity measurement, the identification and trigger efficiency corrections of the leptons, and the energy scale of the leptons. For the DY + ≥1 jet case, the major uncertainties come from the jet energy scale and the unfolding model. The estimates of systematic uncertainties for the inclusive differential cross sections in p T ( ) for various m ranges are shown in Fig. 3.
When calculating the cross section ratios, for each p T ( ) or ϕ * η bin, all uncertainties are taken as fully correlated between the numerator and the denominator, except the data and Monte Carlo statistical uncertainties. The total uncertainty corresponds to the quadratic sum of the sources. The estimates of systematic uncertainties for the ratios of the inclusive p T ( ) distributions are shown in Fig. 4.

Theory predictions
The measured data are compared with the MG5 aMC + PYTHIA 8 baseline sample described in Section 4. The QCD scale uncertainties are estimated by varying the renormalisation and factorisation scales simultaneously by factors of 2 and 1/2 (omitting the variations in opposite direction and taking the envelope). The strong coupling (α S ) and PDF uncertainties are estimated as the standard deviation of weights from the replicas provided in the NNPDF 3.0 PDF set [62]. Relative   Relative  Figure 4: Estimates of the uncertainties in inclusive differential cross section ratios in p T ( ) for invariant mass ranges with respect to the peak region 76 < m < 106 GeV: 50 < m < 76 GeV (upper left), 106 < m < 170 GeV (upper right), 170 < m < 350 GeV (lower left), and 350 < m < 1000 GeV (lower right). The black line is the quadratic sum of the colored lines.
An event sample at NNLO with a jet merging method is generated with MINNLO PS [81,82]. The coupling α S is evaluated independently at each vertex at a scale that depends on the kinematic configuration. Sudakov form factors are used to interpolate between the scales. The NNLO version of the NNPDF 3.1 PDF set [83] is used along with the PYTHIA version 8 [60] for the parton showers based on the CP5 tune [84] and multiparton interactions (MPI), but including a harder primordial k T of 2.2 GeV obtained from tuning the k T parameter to describe the observed ϕ * η distribution of Ref.
[26]. The results are also compared with a third prediction from the parton branching (PB) TMD method [14,15] obtained from CASCADE 3 [85]. This prediction is of particular interest since the initial-state parton showers are fully determined by TMD and their backward PB evolution, and therefore are free of tuning parameters. The matrix element calculation is performed at NLO for Z + 0 jet using MG5 aMC for the inclusive distributions (labelled MG5 aMC (0 jet at NLO)+ PB (CASCADE)), and for Z + 1 jet for the distributions where one jet is required in the final state (labelled MG5 aMC (1 jet at NLO)+ PB (CASCADE)). Initial-state parton showers, provided by the PB TMD method are matched to the NLO matrix element [86], using the latest TMD PB set: PB-NLO-HERAI+II-2018-set2 [87]. The final parton shower, hadronization, and QED FSR steps are performed with PYTHIA 6 [88]. This approach is equivalent to the inclusion of the next-to-leading logarithmic soft-gluon resummation on top of the fixed-order NLO calculations. The theoretical uncertainties in the cross section are estimated by variation of scales and from TMD uncertainties. This approach is expected to describe the inclusive cross section at low p T ( ) (< 20 GeV) well, and to fail for larger p T ( ), since higher-order matrix element contributions are missing, as already observed for the Z boson mass peak range [26]. Recently, this approach has been developed to include multi-jet merging [89] at LO, which allows a larger p T ( ) region to be described as well.
A fourth prediction is based on an independent approach relying on TMDs obtained from fits to DY and Z boson measurements at different energies [90,91] using an NNLO evolution. The corresponding numerical evaluations are provided by the ARTEMIDE 2.02 code [92]. The resummation corresponds to an N 3 LL approximation. The uncertainty is obtained from scale variations. Due to the approximation of ordering among the scales, the prediction has a limited range of validity for the calculation of: p T ( ) < 0.2 m . Predictions for the ϕ * η cross section dependence as well as the 1 jet case are not provided by ARTEMIDE. The ARTEMIDE sample does not include the QED FSR; a correction is derived from the PYTHIA 8 shower in the MG5 aMC sample. The uncertainty is derived by taking the difference with respect to corrections derived from the POWHEG sample described in Ref. [26]. This uncertainty is smaller than 1% for p T ( ) < 0.2 m Two more predictions are obtained from the GENEVA 1.0-RC3 program [93][94][95] combining higher-order resummation with a DY calculation at NNLO. Originally, the resummation was carried out at NNLL including partially N 3 LL on the 0-jettiness variable τ 0 [96]. More recently it includes the q T resummation at N 3 LL in the Radish formalism [97,98] for the 0 jet case, whereas it keeps the 1-jettiness resummation for the 1 jet case. Two samples are generated, one in the 0-jettiness approach and one in the q T resummation approach. The calculation uses the PDF4LHC15 NNLO [99] PDF set with α S (m Z ) = 0.118, the world average. The events are showered using a specially modified version of PYTHIA 8, which is also used for nonperturbative effects and QED radiation in the initial and final states using a modified tune based on CUETP8M1. The theoretical uncertainties are estimated by variation of scales and from the resummation as described in Ref. [94]. No uncertainty is assigned to the jetiness resummation.

p T ( ) results
The differential cross sections in p T ( ) are shown in Fig. 5 for invariant mass ranges between 50 GeV and 1 TeV. Because of the lack of precision of the muon transverse momentum measurement at high p T , the cross section measurement in the highest mass range is based on the electron channel only. The ratio of the predictions to the data are presented in Figs. 6, 7 and 8. The comparison with different predictions is discussed later in the text. The ratios of the unfolded distributions for invariant masses outside the Z boson peak to the distribution within the Z boson peak (76 < m < 106 GeV) are shown in Fig. 9, and the comparisons to predictions in Figs. 10, 11 and 12.
The measured cross sections are presented in Fig. 13 as a function of p T ( ) for at least one jet, for the same mass ranges except the highest. Ratios of the predictions to the data are presented in Figs. 14 and 15. The ratio of these differential cross sections for various mass ranges with respect to the same distribution in the Z boson peak region are shown in Fig. 16, and the comparisons to predictions in Figs. 17 and 18.
The measurements show that the differential cross sections in p T ( ) are rising from small p T ( ) values up to a maximum between 4 and 6 GeV and then falling towards large p T ( ) (Fig. 5). For these cross sections, the variation of the dilepton invariant mass does not have a visible effect on the peak position (around 5 GeV) or on the rising shape for the values below the peak. However, the increase of m results in a broader distribution for p T ( ) values above the peak. These effects are highlighted by the cross section ratios presented in Fig. 9. It has to be noted that the rising ratio for the lowest m range ( Fig. 9 top left) up to a p T ( ) value of 20 GeV is due to QED radiative effects on the final-state leptons (photon radiations at ∆R( , γ) > 0.1) inducing migrations from the Z mass peak towards lower masses. When a jet with a large transverse momentum is required (Fig. 13), the peak is shifted towards larger p T ( ) values corresponding to the jet selection threshold, here 30 GeV regardless of the m . As in the inclusive case, the distributions become broader for p T ( ) values larger than the peak for increasing m .
A description of these measurements based on QCD requires both multi-gluon resummation and a fixed-order matrix element. The description of the distributions at small p T ( ) values requires an approach taking into account initial-state nonperturbative and perturbative multigluon resummation. The falling behaviour at large p T ( ) is sensitive to hard QCD radiation, which is expected to be well described by matrix element calculations including at least NLO corrections. The size of the QCD radiation is driven by the available kinematic phase space and the value of α S . An increase of m extends the phase space for hard radiations, slightly compensated by the decrease of α S with increasing m . The tail at large p T ( ) is dominated by jet multiplicities above one. For the inclusive cross sections, the resummation effects are concentrated at small p T ( ). The value of the maximum of the distributions is expected to depend weakly on m . In the presence of a hard jet, multiple gluon emissions also affect the perturbative region located in η between the jet and the vector boson. The corresponding cross section measurements therefore provide additional constraints on the resummation treatment in the predictions.
The MG5 aMC + PYTHIA 8 prediction describes the data well globally (Fig. 6), although it predicts a too-small cross section for p T ( ) values below 30 GeV in the inclusive case. This disagreement is more pronounced at higher m and reaches about 20% for masses above 170 GeV. The low-p T ( ) region is sensitive to gluon resummation. In MG5 aMC, the resummation effects are simulated by the parton shower, modelled in PYTHIA 8 depending on parameters tuned on previously published measurements, including DY cross sections in the Z boson mass peak region. It has to be noted that the low p T ( ) spectrum is sensitive to the choice of the tuned parameters [84] and that no related systematic uncertainty is available. The large p T ( ) distributions are well described by MG5 aMC, which relies on NLO matrix elements for 0, 1 and 2 partons in the final state. Nevertheless, MG5 aMC predicts cross sections larger than those observed for the highest p T ( ) values measured in the mass ranges 106 < m < 170 GeV for both the inclusive and 1 jet cases. Since the theoretical uncertainty is dominant in that region, a better agreement might be found using higher-order (e.g., NNLO) multiparton predictions.
The MINNLO PS prediction provides the best global description of the data among the predictions presented in this paper. This approach, based on NNLO matrix element and PYTHIA 8 parton shower and MPI, describes well the large p T ( ) cross sections (Fig. 6) and ratios (Fig. 17), except above 400 GeV, for m around the Z boson peak. The medium and low p T ( ) cross sections are also well described by MINNLO PS which relies on parton showers, a harder primordial k T and Sudakov form factors. The same observation can be made in the one jet case. The inclusion of an NNLO matrix element reduces significantly the scale uncertainties, in particular for the inclusive cross section in for the medium p T ( ) values where the PDF uncertainty becomes significant with respect to other model uncertainties. It has to be noted that no parton shower tune uncertainty is assigned in the case of MINNLO PS as well as in the case of MG5 aMC.
We see that the CASCADE predictions (MG5 aMC + PB(CASCADE)) involving TMDs produce a better description in the low-p T ( ) part than MG5 aMC + PYTHIA 8, which is valid for all m bins. The predicted cross section for medium p T ( ) values is 5 to 10% too low (Fig. 7). What is remarkable is that this prediction is based on TMDs obtained from totally independent data, from a fit to electron-proton deep inelastic scattering measurements performed at HERA. The high p T ( ) part is not described by the Z +0,1 jet matrix element calculations from MG5 aMC with CASCADE due to missing higher fixed-order calculations. The range of p T ( ) values well described extends with increasing m . For the one jet case (Fig. 14), the low-p T ( ) part is mainly dominated by Z + 2 jet events, and the CASCADE predictions are missing the contributions from the double parton scattering. It thus fails to describe the low p T ( ) region. In the low-p T ( ) region of the 1 jet case double parton scattering contributions play a significant role and thus CASCADE without it cannot describe this region. The CASCADE predictions give an overall good description of the ratio measurements (Fig. 17). Recently the predictions have been extended by including multi-jet merging [89] for an improved description of the full p T ( ) spectrum, shown in the Appendix A.
Within its range of validity, p T ( ) < 0.2 m , the ARTEMIDE prediction describes the measurements very well. For all m , the low-p T ( ) distributions predicted by ARTEMIDE, based on TMDs corresponding to an N 3 LL approximation, are in very good agreement with the data, except for the highest masses. Figure 7 shows the prediction with and without QED FSR corrections. This underlines the importance of migrations from the Z boson peak towards lower masses, inducing the peak structure in the p T ( ) ratio distribution of Fig. 9. The remarkable agreement of the ARTEMIDE prediction with the measurement at the Z boson peak is expected since the prediction relies on TMDs fitted on previous DY measurements at the Z boson peak though at lower centre of mass energies. The excellent agreement for higher m confirms the validity of the approach and in particular of the TMD factorization when the mass scale largely dominates over the transverse momentum. No prediction is provided by ARTEMIDE for the 1 jet case nor for the ϕ * η cross section dependence.
Comparisons of the inclusive cross section as a function of p T ( ) with two predictions of GENEVA are presented in Fig. 8 for the inclusive cross sections and in Fig. 15 for the one jet cross sections. The original prediction combining NNLL resummation on the 0-jettiness variable τ 0 (GENEVA-τ) and NNLO corrections does not describe the data well for p T ( ) values below 40 GeV. This too hard p T ( ) spectrum might be related to the choice of α S , as discussed in Ref. [94]. For the high p T ( ) region, which is dominated by the fixed-order effects, the inclusion of NNLO corrections provides a good description of the measured cross section. The more recent GENEVA prediction (GENEVA-q T ), using a q T resummation at N 3 LL, provides a much better description of the measured inclusive cross sections, describing very well the data in the full p T ( ) range except for middle p T ( ) values in the lowest mass bin. Here, as in MINNLO PS case, the inclusion of NNLO corrections provides a significant reduction of the scale uncertainties, leading to very small theory uncertainties in the middle p T ( ) range. The two GENEVA predictions compared with the measured one jet cross sections are similar because both use 1-jettiness in this part of the phase space. This could explain that GENEVA predicts a too hard p T ( ) spectrum, similarly to the 0-jettiness inclusive case.

ϕ * η results
The ϕ * η variable is highly correlated with p T ( ) and it offers a complementary access to the underlying QCD dynamics. Being based only on angle measurements of the final-state charged leptons, the ϕ * η variable can be measured with greater accuracy which allows us to include the muon channel for all m ranges. Figure 19 presents the inclusive differential cross sections in ϕ * η for the same invariant mass ranges as above and comparisons to models. More complete comparisons to model predictions are presented as ratios of the prediction divided by the measurement in Figs. 20 and 21. The results are discussed below. The ratio of these differential cross sections for various m ranges are computed with respect to the same distribution in the Z peak region. They are shown in Fig. 22 and further compared with models in Figs. 23 and 24.
The ϕ * η distributions are monotonic functions, in particular they do not present a peak structure as measured in the p T ( ) distributions. At small values, the ϕ * η distributions contain a plateau whose length decreases with increasing m , and more generally the ϕ * η distributions fall more rapidly with increasing m as clearly shown in Fig. 19. Because the lepton direction is much less affected by QED FSR than the energy, the effect of migrations from the Z boson mass bin towards lower masses is relatively invisible in the ϕ * η shape as highlighted by the ratio distribution in Fig. 22 (upper left).
Since ϕ * η is highly correlated with p T ( ), the comparison of the ϕ * η distributions to theoretical predictions leads to the same basic observations and remarks as related above. The MG5 aMC + PYTHIA 8 prediction describes the measured ϕ * η distributions well globally and predicts a too small cross section in the region sensitive to gluon resummation, i.e., ϕ * η 0.1 on the Z boson mass peak, as shown in Fig. 20. The increase of this disagreement for higher m is also observed, clearly visible in the ratio distributions of Fig. 23.
As for the p T ( ) distributions, the MINNLO PS prediction provides the best global description of the data (Fig. 20). In contrast to the disagreement for p T ( ) above 400 GeV for m around the Z peak that appeared both in the inclusive case (Fig. 6) and in the one jet case (Fig. 17) the large ϕ * η values are well described by MINNLO PS . The inclusion of NNLO corrections reduces scale uncertainties making the PDF uncertainty dominant for medium ϕ * η values in the central m bins. The PDF uncertainty is significantly reduced in the ratio distributions (Fig. 23) leading to remarkable prediction precision of the level of 1.5% in several bins.
The MG5 aMC + PB(CASCADE) prediction describes well the measured shapes for ϕ * η 0.1 in         Details on the presentation of the results are given in Fig. 6 caption.        Figure 15: Comparison of the differential cross sections in p T ( ) to predictions in various invariant mass ranges for the one or more jets case. The measurement is compared with GENEVAτ (left) and GENEVA-q T (right) predictions. Details on the presentation of the results are given in Fig. 8 caption.
all m bins (Fig. 20). This contrasts with the description of the p T ( ) dependence by the same prediction (Fig. 7), owing to the washing out of the details of the p T ( ) distribution in the ϕ * η distribution. The normalisation of the prediction is good for the Z boson mass peak region but underestimates more and more the cross section with increasing m , in a way relatively close to MG5 aMC predictions. The ratio distributions (Fig. 23) also illustrate this, but a compensation effect leads to predictions in agreement over the full ϕ * η range. The measured cross sections as a function of ϕ * η are compared with GENEVA predictions in Fig. 21. Similar to previous discussions of the p T ( ) distributions, GENEVA-q T improves significantly the description of the data with respect to GENEVA-τ. The discrepancy of GENEVA-q T for low p T ( ) values in the two lowest m bins is smoothed here leading to a global agreement everywhere. The cross section ratio distributions of the different m bins over the Z boson mass peak bin, as a function of ϕ * η are shown in Fig. 24. Here both GENEVA predictions provide a good description of the measurements. This indicates that, although the precise shape in ϕ * η is not well reproduced by GENEVA-τ, the scale dependence is well described over the large range covered by the present measurement.
The differential cross section measurements are presented in the HEPData entry [100].  Figure 18: Comparison of the ratios of differential cross sections in p T ( ) for one or more jets in various invariant mass ranges with respect to the peak region 76 < m < 106 GeV. The measured ratio is compared with GENEVA-τ (left) and GENEVA-q T (right) predictions. Details on the presentation of the results are given in Fig. 8 caption.  Figure 21: Comparison of the differential cross sections in ϕ * η ( ) to predictions in various m ranges. The measurement is compared with GENEVA-τ (left) and GENEVA-q T (right) predictions. Details on the presentation of the results are given in Fig. 8 Figure 22: Ratios of differential cross sections in ϕ * η ( ) for invariant mass ranges with respect to the peak region 76 < m < 106 GeV: 50 < m < 76 GeV (upper left), 106 < m < 170 GeV (upper right), 170 < m < 350 GeV (lower left), and 350 < m < 1000 GeV (lower right). The measurement is compared with MG5 aMC (0, 1, and 2 jets at NLO) + PYTHIA 8 (blue dots), MINNLO PS (green diamonds) and MG5 aMC (0 jet at NLO)+ PB (CASCADE) (red triangles). Details on the presentation of the results are given in Fig. 5 Figure 24: Ratios of differential cross sections in ϕ * η ( ) for invariant m with respect to the peak region 76 < m < 106 GeV. Compared to model predictions from GENEVA-τ (left) and GENEVA-q T (right). Details on the presentation of the results are given in Fig. 8 caption.

Summary
Measurements of differential Drell-Yan cross sections in proton-proton collisions at √ s = 13 TeV in the dielectron and dimuon final states are presented, using data collected with the CMS detector, corresponding to an integrated luminosity of 36.3 fb −1 . The measurements are corrected for detector effects and the two leptonic channels are combined. Differential cross sections in the dilepton transverse momentum, p T ( ), and in the lepton angular variable ϕ * η are measured for different values of the dilepton mass, m , between 50 GeV and 1 TeV. To highlight the evolution with the dilepton mass scale, ratios of these distributions for various masses are presented. In addition, dilepton transverse momentum distributions are shown in the presence of at least one jet within the detector acceptance.
The rising behaviour of the Drell-Yan inclusive cross section at small p T ( ) is attributed to soft QCD radiations, whereas the tail at large p T ( ) is only expected to be well described by models relying on higher-order matrix element calculations. Therefore, this variable provides a good sensitivity to initial-state QCD radiations and can be compared with different predictions relying on matrix element calculations at different orders and using different methods to resum the initial-state soft QCD radiations. The measurements show that the peak in the p T ( ) distribution, located around 5 GeV, is not significantly modified by changing the m value in the covered range. However, for higher values of m above the peak, the p T ( ) distributions fall less steeply.
The ϕ * η variable, highly correlated with p T ( ), offers a complementary access to the underlying QCD dynamics. Since it is based only on angle measurements of the final-state charged leptons, it offers, a priori, measurements with greater accuracy. However, these measurements demonstrate that the ϕ * η distributions discriminate between the models less than the p T ( ) distributions, since they wash out the peak structure of the p T ( ) distributions, which reflect the initial-state QCD radiation effects in a more detailed way.
This publication presents comparisons of the measurements to six predictions using different treatments of soft initial-state QCD radiations. Two of them, MG5 aMC + PYTHIA 8 and MINNLO PS , are based on a matrix element calculation merged with parton showers. Two others, ARTEMIDE and CASCADE use transverse momentum dependent parton distributions (TMD). Finally, GENEVA combines a higher-order resummation with a Drell-Yan calculation at next-to-next-to-leading order (NNLO), in two different ways. One carries out the resummation at next-to-next-to-leading logarithm in the 0-jettiness variable τ 0 , the other at N 3 LL in the q T variable.
The comparison of the measurement with the MG5 aMC + PYTHIA 8 Monte Carlo predictions using matrix element calculations including Z + 0, 1, 2 partons at next-to-leading order (NLO) merged with a parton shower, shows generally good agreement, except at p T ( ) values below 10 GeV both for the inclusive and one jet cross sections. This disagreement is enhanced for masses away from the Z mass peak and is more pronounced for the higher dilepton masses, reaching 20% for the highest mass bin.
The MINNLO PS prediction provides the best global description of the data among the predictions presented in this paper, both for the inclusive and the one jet cross sections. This approach, based on NNLO matrix element and PYTHIA 8 parton shower and MPI, describes well the large p T ( ) cross sections and ratios, except for p T ( ) values above 400 GeV for dilepton masses around the Z mass peak. A good description of the medium and low p T ( ) cross sections is obtained using a modified primordial k T parameter of the CP5 parton shower tune. from a fit to electron-proton deep inelastic scattering measurements performed at HERA. These TMDs are merged with NLO matrix element calculations. Low p T ( ) values are globally well described but with too low cross sections at medium p T ( ) values. This discrepancy increases with increasing m in a way similar to the MG5 aMC + PYTHIA 8 predictions. The high part of the p T ( ) distribution is not described by CASCADE due to missing higher fixed-order terms. The model can not describe the low p T ( ) region of the cross section in the presence of one jet due to the missing double parton scattering contributions. The recent inclusion of multi-jet merging allows a larger p T ( ) region to be described as well. ARTEMIDE provides predictions based on TMDs extracted from previous measurements including the Drell-Yan transverse momentum cross section at the LHC at the Z mass peak. By construction, the validity of ARTEMIDE predictions are limited to the range p T ( ) < 0.2 m . In that range, they describe well the present measurements up to the highest dilepton masses.
The GENEVA prediction, combining resummation in the 0-jettiness variable τ 0 (GENEVA-τ) and NNLO matrix element does not describe the measurement well for p T ( ) values below 40 GeV. For the high p T ( ) region the inclusion of NNLO in the matrix element provides a good description of the measured cross section. The recent GENEVA prediction (GENEVA-q T ), using a q T resummation, provides a much better description of the measured inclusive cross sections, describing very well the data in the full p T ( ) range except for middle p T ( ) values in the lowest mass bin. Both GENEVA approaches predict too hard p T ( ) spectra for the one jet cross sections.
The ratio distributions presented in this paper confirm most of the observations based on the comparison between the measurement and the predictions at the cross section level. The observed scale dependence is well described by the different models. Furthermore the partial cancellation of the uncertainties in the cross section ratios allows a higher level of precision to be reached for both the measurement and the predictions.
The present analysis shows the relevance of measuring the Drell-Yan cross section in a wide range in dilepton masses to probe the interplay between the transverse momentum and the mass scales of the process. Important theoretical efforts have been made during the last decade to improve the detailed description of high energy processes involving multiple scales and partonic final states. The understanding of the Drell-Yan process directly benefited from these developments. The present paper shows that they individually describe the measurements well in the regions they were designed for. Nevertheless, no model is able to reproduce all dependencies over the complete covered range. Further progress might come from combining these approaches.     [30] ATLAS Collaboration, "Measurement of the low-mass Drell-Yan differential cross section at √ s = 7 TeV using the ATLAS detector", JHEP 06 (2014) 112, doi:10.1007/JHEP06(2014)112, arXiv:1404.1212.

A Comparisons to other models
In this section, comparisons of the obtained measurement results with predictions from a more recent parton branching (PB) TMD method from CASCADE are presented. The predictions are based on MG5 aMC ME up to three partons at LO in QCD with multi-jet merging [89]. The ratio of the predictions over the data are presented in Fig. A.1. The comparisons to predictions for the ratio of the cross sections for invariant masses outside the Z boson peak to the distribution within the Z boson peak (76 < m < 106 GeV) are shown in Fig. A.2.