Measurement of differential cross sections for ${\text {Z}}$ boson production in association with jets in proton-proton collisions at $\sqrt{s} = 13\,\text {TeV}$

Abstract

The production of a ${\text {Z}}$ boson, decaying to two charged leptons, in association with jets in proton-proton collisions at a centre-of-mass energy of 13$\,\text {TeV}$ is measured. Data recorded with the CMS detector at the LHC are used that correspond to an integrated luminosity of 2.19$\,\text {fb}^\text {-1}$. The cross section is measured as a function of the jet multiplicity and its dependence on the transverse momentum of the ${\text {Z}}$ boson, the jet kinematic variables (transverse momentum and rapidity), the scalar sum of the jet momenta, which quantifies the hadronic activity, and the balance in transverse momentum between the reconstructed jet recoil and the ${\text {Z}}$ boson. The measurements are compared with predictions from four different calculations. The first two merge matrix elements with different parton multiplicities in the final state and parton showering, one of which includes one-loop corrections. The third is a fixed-order calculation with next-to-next-to-leading order accuracy for the process with a ${\text {Z}}$ boson and one parton in the final state. The fourth combines the fully differential next-to-next-to-leading order calculation of the process with no parton in the final state with next-to-next-to-leading logarithm resummation and parton showering.

1 Introduction

Measurements of vector boson production in association with jets provide fundamental tests of quantum chromodynamics (QCD). The high centre-of-mass energy at the CERN LHC allows the production of an electroweak boson along with a large number of jets with large transverse momenta. A precise knowledge of the kinematic distributions in processes with large jet multiplicity is essential to exploit the potential of the LHC experiments. Comparison of the measurements with predictions motivates additional Monte Carlo (MC) generator development and improves our understanding of the prediction uncertainties. Furthermore, the production of a massive vector boson together with jets is an important background to a number of standard model (SM) processes (production of a single top quark, ${{\text {t}}\overline{{\text {t}}}}$, and Higgs boson as well as vector boson fusion and WW scattering) as well as to searches for physics beyond the SM, e.g. supersymmetry. Leptonic decay modes of the vector bosons are often used in the measurement of SM processes and searches for physics beyond the SM since they have a sufficiently high branching fraction and clean signatures that provide a strong rejection of backgrounds. Differential cross sections for the associated production of a ${\text {Z}}$ boson with hadronic jets have been previously measured by the ATLAS, CMS, and LHCb Collaborations in proton-proton collisions at centre-of-mass energies of 7 [1,2,3,4], 8 [5,6,7] and 13 [8] $\,\text {TeV}$, and by the CDF and D0 Collaborations in proton-antiproton collisions at 1.96$\,\text {TeV}$ [9, 10].

In this paper, we present measurements of the cross section multiplied by the branching fraction for the production of a ${\text {Z}}/\gamma ^*$ boson in association with jets and its subsequent decay into a pair of oppositely charged leptons ($\ell ^+\ell ^-$) in proton-proton collisions at a centre-of-mass energy of 13$\,\text {TeV}$. The measurements from the two final states, with an electron–positron pair (electron channel) and with a muon–antimuon pair (muon channel), are combined. The measurements are performed with data from the CMS detector recorded in 2015 at the LHC corresponding to 2.19$\,\text {fb}^\text {-1}$ of integrated luminosity. For convenience, ${\text {Z}}/\gamma ^*$ is denoted as ${\text {Z}}$. In this paper a ${\text {Z}}$ boson is defined as a pair of oppositely charged muons or electrons with invariant mass in the range $91\pm 20\text {GeV}$. This range is chosen to have a good balance between the signal acceptance, the rejection of background processes, and the ratio of ${\text {Z}}$ boson to $\gamma ^*$ event yields. It is also consistent with previous measurements [4,5,6] and eases comparisons.

The cross section is measured as a function of the jet multiplicity ($N_{\text {jets}}$), transverse momentum ($p_{\mathrm {T}}$) of the ${\text {Z}}$ boson, and of the jet transverse momentum and rapidity (y) of the first, second, and third jets, where the jets are ordered by decreasing $p_{\mathrm {T}}$. Furthermore, the cross section is measured as a function of the scalar sum of the jet transverse momenta ($H_{\mathrm {T}}$) for event samples with at least one, two, and three jets. These observables have been studied in previous measurements. In addition, we study the balance in transverse momentum between the reconstructed jet recoil and the ${\text {Z}}$ boson for the different jet multiplicities and two ${\text {Z}}$ boson $p_{\mathrm {T}}$ regions ($p_{\mathrm {T}} ({\text {Z}}) < 50\text {GeV}$ and $p_{\mathrm {T}} ({\text {Z}}) > 50\text {GeV}$).

2 The CMS detector

The central feature of the CMS apparatus is a superconducting solenoid of 6$\text {\,m}$ internal diameter, providing a magnetic field of 3.8$\text {\,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 pseudorapidity coverage provided by the barrel and endcap detectors up to $|\eta |=5$. The electron momentum is estimated by combining the energy measurement in the ECAL with the momentum measurement in the tracker. The momentum resolution for electrons with $p_{\mathrm {T}} \approx 45\text {GeV}$ from ${\text {Z}} \rightarrow \mathrm {e}\mathrm {e}$ decays ranges from 1.7% for nonshowering electrons in the barrel region ($|\eta |< 1.444$) to 4.5% for showering electrons in the endcaps ($1.566< |\eta | < 3$) [11]. When combining information from the entire detector, the jet energy resolution is 15% at 10$\text {GeV}$, 8% at 100$\text {GeV}$, and 4% at 1$\,\text {TeV}$, to be compared to about 40, 12, and 5% obtained when only the ECAL and HCAL calorimeters are used. Muons are measured in the pseudorapidity range $|\eta | < 2.4$, with detection planes made using three technologies: drift tubes, cathode strip chambers, and resistive plate chambers. Matching muons to tracks measured in the silicon tracker results in a relative transverse momentum resolution for muons with $20< p_{\mathrm {T}} < 100\text {GeV}$ of 1.3–2.0% in the barrel and better than 6% in the endcaps. The $p_{\mathrm {T}}$ resolution in the barrel is better than 10% for muons with $p_{\mathrm {T}}$ up to 1$\,\text {TeV}$  [12].

Events of interest are selected using a two-tiered trigger system [13]. 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$\text {\,kHz}$ within a time interval of less than 4$\,\mu \text {s}$. The second level, known as the high-level trigger (HLT), 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$\text {\,kHz}$ before data storage.

3 Observables

The cross section is measured for jet multiplicities up to 6 and differentially as a function of the transverse momentum of the ${\text {Z}}$ boson and as a function of several jet kinematic variables, including the jet transverse momentum, rapidity, and the scalar sum of jet transverse momenta.

Jet kinematic variables are measured for event samples with at least one, two, and three jets. In the following, the jet multiplicity will be referred to as “inclusive” to designate events with at least N jets and as “exclusive” for events with exactly N jets.

The balance between the ${\text {Z}}$ boson and jet transverse momenta is also studied via the $p_{\mathrm {T}}$ balance observable $p_{\mathrm {T}} ^{\text {bal}} = |{\vec p}_{\mathrm {T}} ({\text {Z}}) + \sum _{\text {jets}} {\vec p}_{\mathrm {T}} (\text {j}_i) |$, and the so-called jet-${\text {Z}}$ balance $\text {JZB} = |\sum _{\text {jets}} {\vec p}_{\mathrm {T}} (\text {j}_i) |- |{\vec p}_{\mathrm {T}} ({\text {Z}}) |$, where the sum runs over jets with $p_{\mathrm {T}} > 30\text {GeV}$ and $|y | < 2.4$ [14, 15]. The hadronic activity not included in the jets will lead to an imbalance that translates into $p_{\mathrm {T}} ^{\text {bal}}$ and $\text {JZB}$ values different from zero. It includes the activity in the forward region ($|y | > 2.4$), which is the dominant contribution according to simulation. Gluon radiation in the central region that is not clustered in a jet with $p_{\mathrm {T}} >30\text {GeV}$ will also contribute to the imbalance. Hadronic activity not included in the jets will lead to a shift of the $p_{\mathrm {T}} ^{\text {bal}}$ distribution peak to larger values. The $\text {JZB}$ variable distinguishes between two configurations, one where transverse momentum due to the unaccounted hadronic activity is in the direction of the ${\text {Z}}$ boson and another where it is in the opposite direction. Events in the first configuration that have a large imbalance will populate the positive tail of the $\text {JZB}$ distribution, while those in the second configuration populate the negative tail.

The distribution of $p_{\mathrm {T}} ^{\text {bal}}$ is measured for events with minimum jet multiplicities of 1, 2, and 3. To separate low and high jet multiplicity events without $p_{\mathrm {T}}$ and y constraints on the jets, the $\text {JZB}$ variable is also studied for $p_{\mathrm {T}} ({\text {Z}})$ below and above 50$\text {GeV}$.

The ${\text {Z}}$ boson transverse momentum $p_{\mathrm {T}} ({\text {Z}})$ can be described via fixed-order calculations in perturbative QCD at high values, while at small transverse momentum this requires resummation of multiple soft-gluon emissions to all orders in perturbation theory [16, 17]. The measurement of the distribution of $p_{\mathrm {T}} ({\text {Z}})$ for events with at least one jet, due to the increased phase space for soft gluon radiation, leads to an understanding of the balance in transverse momentum between the jets and the ${\text {Z}}$ boson, and can be used for comparing theoretical predictions that treat multiple soft-gluon emissions in different ways.

4 Phenomenological models and theoretical calculations

The measured ${\text {Z}}+ \text { jets}$ cross section is compared to four different calculations: two merging matrix elements (MEs) with various final-state parton multiplicities together with parton showering; a third with a combination of next-to-next-to-leading order (NNLO) calculation with next-to-next-to-leading logarithmic (NNLL) resummation and with parton showering; and a fourth with fixed-order calculation.

The first two calculations use MadGraph 5_amc@nlo version 2.2.2 (denoted MG5_aMC) [18], which is interfaced with pythia8 (version 8.212) [19]. pythia8 is used to include initial- and final-state parton showers and hadronisation. Its settings are defined by the CUETP8M1 tune [20], in particular the NNPDF 2.3 [21] leading order (LO) parton distribution function (PDF) is used and the strong coupling $\alpha _S (m_{{\text {Z}}})$ is set to 0.130. The first calculation includes MEs computed at LO for the five processes $\mathrm {p}\mathrm {p}\rightarrow {\text {Z}}+ N \text { jets}$, $N=0\ldots 4$ and matched to the parton shower using the $k_{\mathrm {T}}$-MLM [22, 23] scheme with the matching scale set at 19$\text {GeV}$. In the ME calculation, the NNPDF 3.0 LO PDF [24] is used and $\alpha _S (m_{{\text {Z}}})$ is set to 0.130 at the ${\text {Z}}$ boson mass scale. The second calculation includes MEs computed at NLO for the three processes $\mathrm {p}\mathrm {p}\rightarrow {\text {Z}}+ N \text { jets}$, $N=0\ldots 2$ and merged with the parton shower using the FxFx [25] scheme with the merging scale set at 30$\text {GeV}$. The NNPDF 3.0 next-to-leading order (NLO) PDF is used and $\alpha _S (m_{{\text {Z}}})$ is set to 0.118. This second calculation is also employed to derive nonperturbative corrections for the fixed-order prediction discussed in the following.

The third calculation uses the geneva 1.0-RC2 MC program (GE), where an NNLO calculation for Drell–Yan production is combined with higher-order resummation [26, 27]. Logarithms of the 0-jettiness resolution variable, ${\tau }$, also known as beam thrust and defined in Ref. [28], are resummed at NNLL including part of the next-to-NNLL (N$^{3}$LL) corrections. The accuracy refers to the $\tau$ dependence of the cross section and is denoted NNLL’$_\tau$. The PDF set PDF4LHC15 NNLO [29] is used for this calculation and $\alpha _S (m_{{\text {Z}}})$ is set to 0.118. The resulting parton-level events are further combined with parton showering and hadronisation provided by pythia8 using the same tune as for MG5_aMC.

Finally, the distributions measured for $N_{\text {jets}}\ge 1$ are compared with the fourth calculation performed at NNLO accuracy for ${\text {Z}}+1$ jet using the N-jettiness subtraction scheme ($N_{\text {jetti}}$) [30, 31]. The PDF set CT14 [32] is used for this calculation. The nonperturbative correction obtained from MG5_aMC and pythia8 is applied. It is calculated for each bin of the measured distributions from the ratio of the cross section values obtained with and without multiple parton interactions and hadronisation. This correction is less than 7%.

Given the large uncertainty in the LO calculation for the total cross section, the prediction with LO MEs is rescaled to match the $\mathrm {p}\mathrm {p}\rightarrow {\text {Z}}$ cross section calculated at NNLO in $\alpha _S$ and includes NLO quantum electrodynamics (QED) corrections with fewz  [33] (version 3.1b2). The values used to normalise the cross section of the MG5_aMC predictions are given in Table 1. All the numbers correspond to a 50$\text {GeV}$ dilepton mass threshold applied before QED final-state radiation (FSR). With fewz, the cross section is computed in the dimuon channel, using a mass threshold applied after QED FSR, but including the photons around the lepton at a distance $R = \sqrt{(\varDelta \eta )^2+(\varDelta \phi )^2}$ smaller than 0.1. The number given in the table includes a correction computed with the LO sample to account for the difference in the mass definition. This correction is small, $+0.35\%$. When the mass threshold is applied before FSR, the cross section is assumed to be the same for the electron and muon channels.

Table 1 Values of the $\mathrm {p}\mathrm {p}\rightarrow \ell ^+\ell ^-$ total cross section used for the calculation in data-theory comparison plots. The cross section used, the cross section from the MC generator (“native”), and the ratio of the two (k) are provided. The phase space of the sample to which the cross section values correspond is indicated in the second column

Uncertainties in the ME calculation (denoted theo. unc. in the figure legends) are estimated for the NLO MG5_aMC, NNLO, and geneva calculations following the prescriptions recommended by the authors of the respective generators. The uncertainty coming from missing terms in the fixed-order calculation is estimated by varying the normalisation ($\mu _{\mathrm {R}}$) and factorisation ($\mu _{\mathrm {F}}$) scales by factors 0.5 and 2. In the case of the FxFx-merged sample, the envelope of six combinations of the variations is considered, the two combinations where one scale is varied by a factor 0.5 and the other by a factor 2 are excluded. In the case of the NNLO and geneva samples the two scales are varied by the same factor, leading to only two combinations. For geneva, the uncertainty is symmetrised by using the maximum of the up and down uncertainties for both cases. The uncertainty from the resummation is also estimated and added in quadrature. It is estimated using six profile scales [34, 35], as described in Ref. [26]. Uncertainties in PDF and $\alpha _S$ values are also estimated in the case of the FxFx-merged sample. The PDF uncertainty is estimated using the set of 100 replicas of the NNPDF 3.0 NLO PDF, and the uncertainty in the $\alpha _S$ value used in the ME calculation is estimated by varying it by $\pm 0.001$. These two uncertainties are added in quadrature to the ME calculation uncertainties. For geneva and NLO MG5_aMC all these uncertainties are obtained using the reweighting method [26, 36] implemented in these generators.

5 Simulation

MC event generators are used to simulate proton-proton interactions and produce events from signal and background processes. The response of the detector is modeled with Geant4 [37]. The ${\text {Z}}(\rightarrow \ell ^+ \ell ^-) + \text { jets}$ process is generated with NLO MG5_aMC interfaced with pythia8, using the FxFx merging scheme as described in Sect. 4. The sample includes the ${\text {Z}}\rightarrow \mathrm {\tau }^+\mathrm {\tau }^-$ process, which is considered a background. Other processes that can give a final state with two oppositely charged same-flavour leptons and jets are $\mathrm {W}\mathrm {W}$, $\mathrm {W}{\text {Z}}$, ${\text {Z}}{\text {Z}}$, ${\text {t}} \overline{{\text {t}}}$ pairs, and single top quark production. The ${\text {t}} \overline{{\text {t}}}$ and single top quark backgrounds are generated using powheg version 2 [38,39,40,41] interfaced with pythia8. Background samples corresponding to diboson electroweak production (denoted VV in the figure legends) [42] are generated at NLO with powheg interfaced to pythia8 ($\mathrm {W}\mathrm {W}$), MG5_aMC interfaced to pythia8 or pythia8 alone ($\mathrm {W}{\text {Z}}$ and ${\text {Z}}{\text {Z}}$). The background sample corresponding to $\mathrm {W}+ \text { jets}$ production ($\mathrm {W}$) is generated at NLO using MG5_aMC interfaced with pythia8, utilizing the FxFx merging scheme.

The events collected at the LHC contain multiple superimposed proton-proton collisions within a single beam crossing, an effect known as pileup. Samples of simulated pileup are generated with a distribution of proton-proton interactions per beam bunch crossing close to that observed in data. The number of pileup interactions, averaging around 20, varies with the beam conditions. The correct description of pileup is ensured by reweighting the simulated sample to match the number of interactions measured in data.

6 Object reconstruction and event selection

The particle-flow (PF) algorithm [43] is used to reconstruct the events. It combines the information from the various elements of the CMS detector to reconstruct and identify each particle in the event. The reconstructed particles are called PF candidates. If several primary vertices are reconstructed, we use the one with the largest quadratic sum of associated track transverse momenta as the vertex of the hard scattering and the other vertices are assumed to be pileup.

The online trigger selects events with two isolated electrons (muons) with transverse momenta of at least 17 and 12 (17 and 8) $\text {GeV}$. After offline reconstruction, the leptons are required to satisfy $p_{\mathrm {T}} > 20 \text {GeV}$ and $|\eta | < 2.4$. We require that the two electrons (muons) with highest transverse momenta form a pair of oppositely charged leptons with an invariant mass in the range $91\pm 20\text {GeV}$. The transition region between the ECAL barrel and endcap ($1.444< |\eta |< 1.566$) is excluded in the reconstruction of electrons and the missing acceptance is corrected to the full $|\eta | < 2.4$ region. The reconstruction of electrons and muons is described in detail in Refs. [11, 12]. The identification criteria applied for electrons and muons are identical to those described in the Ref. [6] except for the thresholds of the isolation variables, which are optimised for 13$\,\text {TeV}$ centre-of-mass energy in our analysis. Electrons (muons) are considered isolated based on the scalar $p_{\mathrm {T}}$ sum of the nearby PF candidates with a distance $R = \sqrt{(\varDelta \eta )^2+(\varDelta \phi )^2} < 0.3$ (0.4). The scalar $p_{\mathrm {T}}$ sum must be less than 15 (25)% of the electron (muon) transverse momentum. We also correct the simulation for differences from data in the trigger, and the lepton identification, reconstruction and isolation efficiencies. These corrections, which depend on the run conditions, are derived using data taken during the run period, and they typically amount to 1–2% for the reconstruction and identification efficiency and 3–5% for the trigger efficiency.

Jets at the generator level are defined from the stable particles ($c\tau > 1\,\text {cm}$), neutrinos excluded, clustered with the anti-$k_{\mathrm {T}}$ algorithm [44] using a radius parameter of 0.4. The jet four-momentum is obtained according to the E-scheme [45] (vector sum of the four-momenta of the constituents). In the reconstructed data, the algorithm is applied to the PF candidates. The technique of charged-hadron subtraction [43] is used to reduce the pileup contribution by removing charged particles that originate from pileup vertices. The jet four-momentum is corrected for the difference observed in the simulation between jets built from PF candidates and generator-level particles. The jet mass and direction are kept constant for the corrections, which are functions of the jet $\eta$ and $p_{\mathrm {T}}$, as well as the energy density and jet area quantities defined in Ref. [46, 47]. The latter are used in the correction of the energy offset introduced by the pileup interactions. Further jet energy corrections are applied for differences between data and simulation in the pileup in zero-bias events and in the $p_{\mathrm {T}}$ balance in dijet, ${\text {Z}}+\text { jet}$, and $\gamma +\text { jet}$ events. Since the $p_{\mathrm {T}}$ balance in ${\text {Z}}+\text { jet}$ events is one of the observables we are measuring in this paper, it is important to understand how it is used in the jet calibration. The balance is measured for events with two objects (jet, $\gamma$, or ${\text {Z}}$ boson) back-to-back in the transverse plane ($|\varDelta \phi - \pi | < 0.34$) associated with a possible third object, a soft jet. The measurement is made for various values of $\rho =p_{\mathrm {T}} ^{\text {soft jet}}/p_{\mathrm {T}} ^{\text {ref}}$, running from 0.1 to 0.3, and extrapolated to $\rho = 0$. In the case the back-to-back objects are a jet and a boson, $p_{\mathrm {T}} ^{\text {ref}}$ is defined as the transverse momentum of the boson, while in the case of two jets it is defined as the average of their transverse momenta. All jets down to $p_{\mathrm {T}} = 5$ or 10$\text {GeV}$, including jets reconstructed in the forward calorimeter, are considered for the soft jet. The data-simulation adjustment is therefore done for ideal topologies with only two objects, whose transverse momenta must be balanced. The jet calibration procedure is detailed in the Ref. [48]. In this measurement, jets are further required to satisfy the loose identification criteria defined in Ref. [49]. Despite the vertex requirement used in the jet clustering some jets are reconstructed from pileup candidates; these jets are suppressed using the multivariate technique described in Ref. [50]. Jets with $p_{\mathrm {T}} >30\text {GeV}$ and $|y |<2.4$ are used in this analysis.

7 Backgrounds estimation

The contributions from background processes are estimated using the simulation samples described in Sect. 5 and are subtracted from the measured distributions. The dominant background, ${{\text {t}}\overline{{\text {t}}}}$, is also measured from data. This ${{\text {t}}\overline{{\text {t}}}}$ background contributes mainly due to events with two same-flavour leptons. The production cross sections for $\mathrm {e}^+\mathrm {e}^-$ and $\mathrm {\mu ^+}\mathrm {\mu ^-}$ events from ${{\text {t}}\overline{{\text {t}}}}$ are identical to the cross section of $\mathrm {e}^+\mathrm {\mu ^-}$ and $\mathrm {e}^-\mathrm {\mu ^+}$ and can therefore be estimated from the latter. We select events in the ${{\text {t}}\overline{{\text {t}}}}$ control sample using the same criteria as for the measurement, but requiring the two leptons to have different flavours. This requirement rejects the signal and provides a sample enriched in ${\text {t}} \overline{{\text {t}}}$ events. Each of the distributions that we are measuring is derived from this sample and compared with the simulation. This comparison produces a discrepancy for events with at least one jet that we correct by applying a correction factor ${\mathcal {C}}$ to the simulation depending on the event jet multiplicity. These factors, together with their uncertainties, are given in Table 2.

After applying this correction to the simulation, all the distributions considered in this measurement agree with data in the ${{\text {t}}\overline{{\text {t}}}}$ control sample. The agreement is demonstrated with a $\chi ^2$-test. We conclude that a parametrization as a function of the jet multiplicity is sufficient to capture the dependency on the event topology. Remaining sources of uncertainties are the estimate of the lepton reconstruction and selection efficiencies and of the yield of events from processes other than ${{\text {t}}\overline{{\text {t}}}}$ entering in the control region. This yield is estimated from the simulation. Based on the sizes of the statistical uncertainties and background contributions, both these uncertainties are negligible. Therefore, the uncertainty in the correction factor is reduced to the statistical uncertainties in the data and simulation samples.

Table 2 The correction factors (${\mathcal {C}}$) applied to the simulated ${{\text {t}}\overline{{\text {t}}}}$ sample with their uncertainties, which are derived from the statistical uncertainties in the data and simulation samples

Fig. 1

Reconstructed data, simulated signal, and background distributions of the inclusive (left) and exclusive (right) jet multiplicity for the electron (upper) and muon (lower) channels. The background distributions are obtained from the simulation, except for the ${\text {t}} \overline{{\text {t}}}$ contribution which is estimated from the data as explained in the text. The error bars correspond to the statistical uncertainty. In the ratio plots, they include both the uncertainties from data and from simulation. The set of generators described in Sect. 5 has been used for the simulation

Fig. 2

Reconstructed data, simulated signal, and background distributions of the transverse momentum balance between the ${\text {Z}}$ boson and the sum of the jets with at least one jet (left) and three jets (right) for the electron (upper) and muon (lower) channels. The background distributions are obtained from the simulation, except for the ${\text {t}} \overline{{\text {t}}}$ contribution which is estimated from the data as explained in the text. The error bars correspond to the statistical uncertainty. In the ratio plots, they include both the uncertainties from data and from simulation. The set of generators described in Sect. 5 has been used for the simulation

Fig. 3

Reconstructed data, simulated signal, and background distributions of the $\text {JZB}$ variable for the electron (left) and muon (right) channels. The background distributions are obtained from the simulation, except for the ${\text {t}} \overline{{\text {t}}}$ contribution which is estimated from the data as explained in the text. The error bars correspond to the statistical uncertainty. In the ratio plots, they include both the uncertainties from data and from simulation. The set of generators described in Sect. 5 has been used for the simulation

The jet multiplicity distributions in data and simulation are presented in Fig. 1. The background contamination is below 1% for the inclusive cross section, and increases with the number of jets to close to 10% for a jet multiplicity of three and above due to ${\text {t}} \overline{{\text {t}}}$ production. Multijet and $\mathrm {W}$ events could pass the selection if one or two jets are misidentified as leptons. The number of multijet events is estimated from data using a control sample obtained by requiring two same-sign same-flavour lepton candidates, whereas the number of $\mathrm {W}$ events is estimated from simulation. Both contributions are found to be negligible. Figure 2 shows the $p_{\mathrm {T}} ^{\text {bal}}$ distribution separately for electron and muon channels. The ${{\text {t}}\overline{{\text {t}}}}$ background does not peak at the same $p_{\mathrm {T}}$ balance as the signal, and has a broader spectrum. The $\text {JZB}$ distribution is shown in Fig. 3. The ${{\text {t}}\overline{{\text {t}}}}$ background is asymmetric, making a larger contribution to the positive side of the distribution because transverse energy is carried away by neutrinos from $\mathrm {W}$ boson decays, leading to a reduction in the negative term of the $\text {JZB}$ expression. Overall the agreement between data and simulation before the background subtraction is good and differences are within about 10%.

8 Unfolding procedure

The fiducial cross sections are obtained by subtracting the simulated backgrounds from the data distributions and correcting the background-subtracted data distributions back to the particle level using an unfolding procedure, which takes into account detector effects such as detection efficiency and resolution. The unfolding is performed using the D’Agostini iterative method with early stopping [51] implemented in the RooUnfold toolkit [52]. The response matrix describes the migration probability between the particle- and reconstructed-level quantities, including the overall reconstruction efficiencies. It is computed using a ${\text {Z}}+ \text { jets}$ sample simulated with MG5_aMC interfaced with pythia8, using the FxFx merging scheme as described in Sect. 4. The optimal number of iterations is determined separately for each distribution by studying the fluctuations introduced by the unfolding with toy MC experiments generated at each step of the iteration. Final unfolded results have also been checked to be consistent with data-simulation comparisons on detector-folded distributions.

Because of the steep slope at the lower boundary of the jet transverse momentum distributions and in order to improve its accuracy, the unfolding is performed for these distributions using histograms with two additional bins, [20, 24] and $[24, 30]\text {GeV}$, below the nominal $p_{\mathrm {T}}$ threshold. The additional bins are discarded after the unfolding

The particle-level values refer to the stable leptons from the decay of the ${\text {Z}}$ boson and to the jets built from the stable particles ($\hbox {c}\tau >1\,\text {cm}$) other than neutrinos using the same algorithm as for the measurements. The momenta of all the photons whose R distance to the lepton axis is smaller than 0.1 are added to the lepton momentum to account for the effects of the final-state radiation; the leptons are said to be “dressed”. The momentum of the ${\text {Z}}$ boson is taken to be the sum of the momenta of the two highest-$p_{\mathrm {T}}$ electrons (or muons). The phase space for the cross section measurement is restricted to events with a ${\text {Z}}$ boson mass between 71 and 111$\text {GeV}$ and both leptons with $p_{\mathrm {T}} > 20\text {GeV}$ and $|\eta | < 2.4$. Jets are required to have $p_{\mathrm {T}} > 30\text {GeV}$, $|y | < 2.4$ and a spatial separation from the dressed leptons of $R > 0.4$.

9 Systematic uncertainties

The systematic uncertainties are propagated to the measurement by varying the corresponding simulation parameters by one standard deviation up and down when computing the response matrix. The uncertainty sources are independent, and the resulting uncertainties are therefore added in quadrature. Tables 3, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19 and 20 present the uncertainties for each differential cross section.

The dominant uncertainty comes from the jet energy scale (JES). It typically amounts to 5% for a jet multiplicity of one and increases with the number of reconstructed jets. The uncertainty in the jet resolution (JER), which is responsible for the bin-to-bin migrations that is corrected by the unfolding, is estimated and the resulting uncertainty is typically 1%.

The most important uncertainty after the JES arises from the measured efficiency (Eff) of trigger, lepton reconstruction, and lepton identification, which results in a measurement uncertainty of about 2% up to 4% for events with leptons of large transverse momenta. The uncertainty in the measurement of the integrated luminosity (Lumi) is 2.3% [53]. The resulting uncertainty on the measured distributions is 2.3%, although the uncertainty is slightly larger in regions that contain background contributions that are estimated from simulation.

The largest background contribution to the uncertainty (Bkg) comes from the reweighting procedure for the ${{\text {t}}\overline{{\text {t}}}}$ simulation, which is estimated to be less than 1% for jet multiplicity below 4. Theoretical contributions come from the accuracy of the predicted cross sections, and include the uncertainties from PDFs, $\alpha _S$ and the fixed-order calculation. Three other small sources of uncertainty are: (1) the lepton energy scale (LES) and resolution (LER), which are below 0.3% in every bin of the measured distributions; (2) the uncertainty in the pileup model, where the 5% uncertainty in the average number of pileup events results in an uncertainty in the measurement smaller than 1%; and (3) the uncertainty in the input distribution used to build the response matrix used in the unfolding and described as follows.

Because of the finite binning a different distribution will lead to a different response matrix. This uncertainty is estimated by weighting the simulation to agree with the data in each distribution and building a new response matrix. The weighting is done using a finer binning than for the measurement. The difference between the nominal results and the results unfolded using the alternative response matrix is taken as the systematic uncertainty, denoted Unf model. An additional uncertainty comes from the finite size of the simulation sample used to build the response matrix. This source of uncertainty is denoted Unf stat in the table and is included in the systematic uncertainty of the measurement.

Table 3 Cross section in exclusive jet multiplicity for the combination of both decay channels and breakdown of the uncertainties

Table 4 Cross section in inclusive jet multiplicity for the combination of both decay channels and breakdown of the uncertainties

Table 5 Differential cross section in $p_{\mathrm {T}} ({\text {Z}})$ ($N_{\text {jets}}\ge 1$) for the combination of both decay channels and breakdown of the uncertainties

Table 6 Differential cross section in $1^{\text {st}}$ jet $p_{\mathrm {T}}$ ($N_{\text {jets}}\ge 1$) for the combination of both decay channels and breakdown of the uncertainties

Table 7 Differential cross section in $2^{\text {nd}}$ jet $p_{\mathrm {T}}$ ($N_{\text {jets}}\ge 2$) for the combination of both decay channels and breakdown of the uncertainties

Table 8 Differential cross section in $3^{\text {rd}}$ jet $p_{\mathrm {T}}$ ($N_{\text {jets}}\ge 3$) for the combination of both decay channels and breakdown of the uncertainties

Table 9 Differential cross section in $1^{\text {st}}$ jet $\vert y \vert$ ($N_{\text {jets}}\ge 1$) for the combination of both decay channels and breakdown of the uncertainties

Table 10 Differential cross section in $2^{\text {nd}}$ jet $\vert y \vert$ ($N_{\text {jets}}\ge 2$) for the combination of both decay channels and breakdown of the uncertainties

Table 11 Differential cross section in $3^{\text {rd}}$ jet $|y |$ ($N_{\text {jets}}\ge 3$) for the combination of both decay channels and breakdown of the uncertainties

Table 12 Differential cross section in $H_{\mathrm {T}}$ ($N_{\text {jets}}\ge 1$) for the combination of both decay channels and breakdown of the uncertainties

Table 13 Differential cross section in $H_{\mathrm {T}}$ ($N_{\text {jets}}\ge 2$) for the combination of both decay channels and breakdown of the uncertainties

Table 14 Differential cross section in $H_{\mathrm {T}}$ ($N_{\text {jets}}\ge 3$) for the combination of both decay channels and breakdown of the uncertainties

Table 15 Differential cross section in $p_{\mathrm {T}} ^{\text {bal}}$ ($N_{\text {jets}}\ge 1$) for the combination of both decay channels and breakdown of the uncertainties

Table 16 Differential cross section in $p_{\mathrm {T}} ^{\text {bal}}$ ($N_{\text {jets}}\ge 2$) for the combination of both decay channels and breakdown of the uncertainties

Table 17 Differential cross section in $p_{\mathrm {T}} ^{\text {bal}}$ ($N_{\text {jets}}\ge 3$) for the combination of both decay channels and breakdown of the uncertainties

Table 18 Differential cross section in $\text {JZB}$ (full phase space) for the combination of both decay channels and breakdown of the uncertainties

Table 19 Differential cross section in ${\text{JZB}}$ ($p_{\mathrm {T}} ({\text {Z}})<50$ GeV) for the combination of both decay channels and breakdown of the uncertainties

Table 20 Differential cross section in $\text {JZB}$ ($p_{\mathrm {T}} ({\text {Z}})>50$ GeV) for the combination of both decay channels and breakdown of the uncertainties

10 Results

The measurements from the electron and muon channels are found to be consistent and are combined using a weighted average as described in Ref. [6]. For each bin of the measured differential cross sections, the results of each of the two measurements are weighted by the inverse of the squared total uncertainty. The covariance matrix of the combination, the diagonal elements of which are used to extract the measurement uncertainties, is computed assuming full correlation between the two channels for all the sources of uncertainty sources except the statistical uncertainties and those associated with lepton reconstruction and identification, which are taken to be uncorrelated. The integrated cross section is measured for different exclusive and inclusive multiplicities and the results are shown in Tables 3 and 4.

The results for the differential cross sections are shown in Figs. 4 to 15 and are compared to the predictions described in Sect. 4. For the two predictions obtained from MG5_aMC and pythia8 the number of partons included in the ME calculation and the order of the calculation is indicated by distinctive labels (“${\le } 4\hbox {j}$ LO” for up to four partons at LO and “${\le } 2\hbox {j}$ NLO” for up to two partons at NLO). The prediction of geneva is denoted as “GE”. The label “PY8” indicates that pythia8 is used in these calculations for the parton showering and the hadronisation. The NNLO ${\text {Z}}+ 1 \text { jet}$ calculation is denoted as $\hbox {N}_{\text {jetti}}$ NNLO in the legends. The measured cross section values along with the uncertainties discussed in Sect. 9 are given in Tables 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19 and 20.

Figure 4 shows the measured cross section as a function of the exclusive (Table 3) and the inclusive (Table 4) jet multiplicities. Agreement between the measurement and the MG5_aMC prediction is observed. The cross section obtained from LO MG5_aMC tends to be lower than NLO MG5_aMC up to a jet multiplicity of 3. The total cross section for ${\text {Z}}(\rightarrow \ell ^+\ell ^-)+\ge 0 \text { jet}, m_{\ell ^+\ell ^-}>50\text {GeV}$ computed at NNLO and used to normalise the cross section of the LO prediction is similar to the NLO cross section as seen in Table 1. The smaller cross section seen when requiring at least one jet is explained by a steeply falling $p_{\mathrm {T}}$ spectrum of the leading jet in the LO prediction. The geneva prediction describes the measured cross section up to a jet multiplicity of 2, but fails to describe the data for higher jet multiplicities, where one or more jets arise from the parton shower. This effect is not seen in the NLO (LO) MG5_aMC predictions, which give a fair description of the data for multiplicities above three (four).

The measured cross section as a function of the transverse momentum of the ${\text {Z}}$ boson for events with at least one jet is presented in Fig. 5 and Table 5. The best model for describing the measurement at low $p_{\mathrm {T}}$, below the peak, is NLO MG5_aMC, showing a better agreement than the NNLL$_{\tau }$’ calculation from geneva. The shape of the distribution in the region below 10$\text {GeV}$ is better described by geneva than by the other predictions, as shown by the flat ratio plot. This kinematic region is covered by events with extra hadronic activity in addition to the jet required by the event selection. The estimation of the uncertainty in the shape in this region shows that it is dominated by the statistical uncertainty, represented by error bars on the plot since the systematic uncertainties are negligible. In the intermediate region, geneva predicts a steeper rise for the distribution than the other two predictions and than the measurement. The high-$p_{\mathrm {T}}$ region, where geneva and NLO MG5_aMC are expected to have similar accuracy (NLO), is equally well described by the two. The LO predictions undershoot the measurement in this region despite the normalisation of the total ${\text {Z}}+ \ge 0 \text { jet}$ cross section to its NNLO value.

The jet transverse momenta for the $1^{\text {st}}$, $2^{\text {nd}}$ and $3^{\text {rd}}$ leading jets can be seen in Figs. 6 and 7 (Tables 6, 7 and 8). The LO MG5_aMC predicted spectrum differs from the measurement, showing a steeper slope in the low $p_{\mathrm {T}}$ region. The same feature was observed in the previous measurements [3, 4]. The comparison with NLO MG5_aMC and $\hbox {N}_{\text {jetti}}$ NNLO calculation shows that adding NLO terms cures this discrepancy. The geneva prediction shows good agreement for the measured $p_{\mathrm {T}}$ of the first jet, while it undershoots the data at low $p_{\mathrm {T}}$ for the second jet. The jet rapidities for the first three leading jets have also been measured and the distributions are shown in Figs. 8 and 9 (Tables 9, 10 and 11). All the predictions are in agreement with data.

Fig. 4

Measured cross section for ${\text {Z}}+\text { jets}$ as a function of the jet exclusive (left) and inclusive (right) multiplicity. The error bars represent the statistical uncertainty and the grey hatched bands represent the total uncertainty, including the systematic and statistical components. The measurement is compared with different predictions, which are described in the text. The ratio of each prediction to the measurement is shown together with the measurement statistical (black bars) and total (black hatched bands) uncertainties and the prediction (coloured bands) uncertainties. Different uncertainties were considered for the predictions: statistical (stat), ME calculation (theo), and PDF together with the strong coupling constant ($\alpha _S$). The complete set was computed for one of the predictions. These uncertainties were added together in quadrature (represented by the $\oplus$ sign in the legend)

Fig. 5

Measured cross section for ${\text {Z}}+\text { jets}$ as a function of the transverse momentum of the ${\text {Z}}$ boson for events with at least one jet. Other details are as mentioned in the Fig. 4 caption

The total jet activity has been measured via the $H_{\mathrm {T}}$ variable. The differential cross section as a function of this observable is presented in Figs. 10 and 11 (Tables 12–14) for inclusive jet multiplicities of 1, 2, and 3. The LO MG5_aMC calculation predicts fewer events than found in the data for the region $H_{\mathrm {T}} < 400\text {GeV}$. For higher jet multiplicities both LO and NLO MG5_aMC are compatible with the measurement, although the contribution in the region $H_{\mathrm {T}} < 400\text {GeV}$ is smaller for LO than for NLO MG5_aMC. The contribution at lower values of $H_{\mathrm {T}}$ is slightly overestimated, but the discrepancy is compatible with the theoretical and experimental uncertainties. The geneva generator predicts a steeper spectrum than measured. For jet multiplicities of at least one, we also compare with $\hbox {N}_{\text {jetti}}$ NNLO, and the level of agreement is similar to that found with NLO MG5_aMC. The uncertainty for $\hbox {N}_{\text {jetti}}$ NNLO is larger than in the jet transverse momentum distribution because of the contribution from the additional jets.

Fig. 6

Measured cross section for ${\text {Z}}+\text { jets}$ as a function of the transverse momentum of the first jet. Other details are as mentioned in the Fig. 4 caption

Fig. 7

Measured cross section for ${\text {Z}}+\text { jets}$ as a function of the transverse momentum of the second (upper) and third (lower) jet. Other details are as mentioned in the Fig. 4 caption

Fig. 8

Measured cross section for ${\text {Z}}+\text { jets}$ as a function of the absolute rapidity of the first (upper) and second (lower) jet. Other details are as mentioned in the Fig. 4 caption

Fig. 9

Measured cross section for ${\text {Z}}+\text { jets}$ as a function of the absolute rapidity of the third jet. Other details are as mentioned in the Fig. 4 caption

Fig. 10

Measured cross section for ${\text {Z}}+\text { jets}$ as a function of the $H_{\mathrm {T}}$ observable for events with at least one jet. Other details are as mentioned in the Fig. 4 caption

Fig. 11

Measured cross section for ${\text {Z}}+\text { jets}$ as a function of the $H_{\mathrm {T}}$ observable of jets for events with at least two (upper) and three (lower) jets. Other details are as mentioned in the Fig. 4 caption

The balance in transverse momentum between the jets and the ${\text {Z}}$ boson, $p_{\mathrm {T}} ^{\text {bal}}$, is shown in Figs. 12 and 13 (Tables 15–17) for inclusive jet multiplicities of 1, 2, and 3. When more jets are included, the peak of $p_{\mathrm {T}} ^{\text {bal}}$ is shifted to larger values. The measurement is in good agreement with NLO MG5_aMC predictions. The slopes of the distributions for the first two jet multiplicities predicted by LO MG5_aMC do not fully describe the data. This observation indicates that the NLO correction is important for the description of hadronic activity beyond the jet acceptance used in this analysis, $p_{\mathrm {T}} >30\text {GeV}$ and $|y |>2.4$. An imbalance in the event, i.e. $p_{\mathrm {T}} ^{\text {bal}}$ not equal to zero, requires two partons in the final state with one of the two out of the acceptance. Such events are described with NLO accuracy for the NLO MG5_aMC sample and LO accuracy for the two other samples. In the case of the geneva simulation, when at least two jets are required, as in the second plot of Fig. 12, the additional jet must come from parton showering and this leads to an underestimation of the cross section, as in the case of the jet multiplicity distribution. When requiring two jets within the acceptance, the NLO MG5_aMC prediction, which has an effective LO accuracy for this observable, starts to show discrepancies with the measurement. The estimated theoretical uncertainties cover the observed discrepancies.

Fig. 12

Measured cross section for ${\text {Z}}+\text { jets}$ as a function of the transverse momentum balance between the ${\text {Z}}$ boson and the accompanying jets for events with at least one (upper) and two (lower) jets. Other details are as mentioned in the Fig. 4 caption

Fig. 13

Measured cross section for ${\text {Z}}+\text { jets}$ as a function of the transverse momentum balance between the ${\text {Z}}$ boson and the accompanying jets for events with at least three jets. Other details are as mentioned in the Fig. 4 caption

The $\text {JZB}$ distribution is shown in Figs. 14 and 15 (Tables 18–20) for the inclusive one-jet events, in the full phase space, and separately for $p_{\mathrm {T}} ({\text {Z}})$ below and above 50$\text {GeV}$. As expected in the high-$p_{\mathrm {T}} ({\text {Z}})$ region, i.e. in the high jet multiplicity sample, the distribution is more symmetric. The NLO MG5_aMC prediction provides a better description of the $\text {JZB}$ distribution than geneva and LO MG5_aMC. This applies to both configurations, $\text {JZB} <0$ and ${}>0$. This observation indicates that the NLO correction is important for the description of hadronic activity beyond the jet acceptance used in this analysis.

Fig. 14

Measured cross section for ${\text {Z}}+\text { jets}$ as a function of the $\text {JZB}$ variable (see text), with no restriction on $p_{\mathrm {T}} ({\text {Z}})$. Other details are as mentioned in the Fig. 4 caption

Fig. 15

Measured cross section for ${\text {Z}}+\text { jets}$ as a function of the $\text {JZB}$ variable (see text), for $p_{\mathrm {T}} ({\text {Z}})<50\,\hbox {GeV}$ (left) and $p_{\mathrm {T}} ({\text {Z}})>50\,\hbox {GeV}$ (right). Other details are as mentioned in the Fig. 4 caption

11 Summary

We have measured differential cross sections for the production of a ${\text {Z}}$ boson in association with jets, where the ${\text {Z}}$ boson decays into two charged leptons with $p_{\mathrm {T}} > 20\text {GeV}$ and $|\eta |<2.4$. The data sample corresponds to an integrated luminosity of 2.19$\,\text {fb}^\text {-1}$ collected with the CMS detector during the 2015 proton-proton LHC run at a centre-of-mass energy of 13$\,\text {TeV}$.

The cross section has been measured as functions of the exclusive and inclusive jet multiplicities up to 6, of the transverse momentum of the ${\text {Z}}$ boson, jet kinematic variables including jet transverse momentum ($p_{\mathrm {T}}$), the scalar sum of jet transverse momenta ($H_{\mathrm {T}}$), and the jet rapidity (y) for inclusive jet multiplicities of 1, 2, and 3. The balance in transverse momentum between the reconstructed jet recoil and the ${\text {Z}}$ boson has been measured for different jet multiplicities. This balance has also been measured separating events with a recoil smaller and larger than the boson $p_{\mathrm {T}}$ using the $\text {JZB}$ variable. Jets with $p_{\mathrm {T}} >30\text {GeV}$ and $|y | < 2.4$ are used in the definition of the different jet quantities.

The results are compared to the predictions of four different calculations. The first two merge matrix elements with different final-state parton multiplicities. The first is LO for multiplicities up to 4, the second NLO for multiplicities up to 2 and LO for a jet multiplicity of 3, and both are based on MG5_aMC. The third is a combination of NNLO calculation with NNLL resummation, based on geneva. The fourth is a fixed order NNLO calculation of one ${\text {Z}}$ boson and one jet. The first three calculations include parton showering, based on pythia8.

The measurements are in good agreement with the results of the NLO multiparton calculation. Even the measurements for events with more than 2 jets agree within the $\approx 10\%$ measurement and 10% theoretical uncertainties, although this part of the calculation is only LO. The multiparton LO prediction does not agree as well as the NLO multiparton one. It exhibits significant discrepancies with data in jet multiplicity and in both transverse momentum and rapidity distributions of the leading jet.

The transverse momentum balance between the ${\text {Z}}$ boson and the hadronic recoil, which is expected to be sensitive to soft-gluon radiation, has been measured for the first time at the LHC. The multiparton LO prediction fails to describe the measurement, while the multiparton NLO prediction provides a very good description for jet multiplicities computed with NLO accuracy.

Inclusive measurement for events with at least one jet are compared with the NNLO ${\text {Z}}+\ge 1 \text { jet}$ fixed order calculation. The agreement is good, even for the $H_{\mathrm {T}}$ observable, which is sensitive to events of different jet multiplicities.

The NNLO+NNLL predictions provide similar agreement for the measurements of the kinematic variables of the two leading jets, but fail to describe observables sensitive to extra jets. At low transverse momentum of the ${\text {Z}}$ boson, the NLO multiparton calculation provides a better description than the NNLO+NNLL calculation, whereas both calculations provide a similar description at high transverse momentum.

The results suggest using multiparton NLO predictions for the estimation of the ${\text {Z}}+\text { jets}$ contribution at the LHC in measurements and searches, and its associated uncertainty.