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

The production of a Z boson, decaying to two charged leptons, in association with jets in proton-proton collisions at a centre-of-mass energy of 13 TeV is measured. Data recorded with the CMS detector at the LHC are used that correspond to an integrated luminosity of 2.19 fb$^{-1}$. The cross section is measured as a function of the jet multiplicity and its dependence on the transverse momentum of the 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 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 Z boson and one parton in the final state. The fourth combines the fully differential next-to-next-to-leading order calculation with next-to-next-to-leading logarithm resummation and parton showering.


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 dependencies 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, tt, 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 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-4], 8 [5-7] and 13 [8] TeV, and by the CDF and D0 Collaborations in proton-antiproton collisions at 1.96 TeV [9, 10].
In this paper, we present measurements of the cross section multiplied by the branching fraction for the production of a Z/γ * boson in association with jets and its subsequent decay into a pair of oppositely charged leptons ( + − ) in proton-proton collisions at a centre-of-mass energy of 13 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 fb −1 of integrated luminosity. For convenience, Z/γ * is denoted as Z. In this paper a Z boson is defined as a pair of oppositely charged muons or electrons with invariant mass in the range 91 ± 20 GeV. This range is chosen to have a good balance between the acceptance of the selection and the ratio of Z boson to γ * event yields. It is also consistent with previous measurements [4-6] and eases comparisons.
The cross section is measured as a function of the jet multiplicity (N jets ), transverse momentum (p T ) of the 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 T . Furthermore, the cross section is measured as a function of the scalar sum of the jet transverse momenta (H 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 Z boson for the different jet multiplicities and two Z boson p T regions (p T (Z) < 50 GeV and p T (Z) > 50 GeV).

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 pseudorapidity coverage provided by the barrel and endcap detectors up to |η| = 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 T ≈ 45 GeV from Z → ee decays ranges from 1.7% for nonshowering electrons in the barrel region to 4.5% for showering electrons in the endcaps [11]. When combining information from the entire detector, the jet energy resolution is 15% at 10 GeV, 8% at 100 GeV, and 4% at 1 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 |η| < 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 T < 100 GeV 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 [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 kHz within a time interval of less than 4 µ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 kHz before data storage.

Observables
The cross section is measured for jet multiplicities up to 6 and differentially as a function of the transverse momentum of the 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 Z boson and jet transverse momenta is also studied via the p T balance observable p bal T = | p T (Z) + ∑ jets p T (j i )|, and the so-called jet-Z balance JZB = |∑ jets p T (j i )| − | p T (Z)|, where the sum runs over jets with p T > 30 GeV and |y| < 2.4 [14,15]. The hadronic activity not included in the jets will lead to an imbalance that translates into p bal T and 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 T > 30 GeV will also contribute to the imbalance. The p bal T observable is particularly relevant for the study of the core of the distribution where the unaccounted hadronic activity leads to a shift in the peak to larger values. The JZB variable distinguishes between two configurations, one where transverse momentum due to the unaccounted hadronic activity is in the direction of the 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 JZB distribution, while those in the second configuration populate the negative tail.
The distribution of p bal T is measured for events with minimum jet multiplicities of 1, 2, and 3. To separate low and high jet multiplicity events without p T and y constraints on the jets, the JZB variable is also studied for p T (Z) below and above 50 GeV.
The Z boson transverse momentum p T (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 T (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 be-Finally, the distributions measured for N jets ≥ 1 are compared with the fourth calculation performed at NNLO accuracy for Z + 1 jet using the N-jettiness subtraction scheme (N 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 pp → Z cross section calculated at NNLO in α 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 GeV dilepton mass threshold applied before QED finalstate 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 = (∆η) 2 + (∆φ) 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 pp → + − 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 (µ R ) and factorisation (µ 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 α 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 α S value used in the ME calculation is estimated by varying it by ±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.

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 Z(→ + − ) + jets process is generated with NLO MG5 aMC interfaced with PYTHIA8, using the FxFx merging scheme as described in Section 4. The sample includes the Z → τ + τ − process, which is considered a background. Other processes that can give a final state with two oppositely charged same-flavour leptons and jets are WW, WZ, ZZ, tt pairs, and single top quark production. The tt 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 (WW), MG5 aMC interfaced to PYTHIA8 or PYTHIA8 alone (WZ and ZZ). The background sample corresponding to W + jets production (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 ob-served 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.

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) GeV. After offline reconstruction, the leptons are required to satisfy p T > 20 GeV and |η| < 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 ± 20 GeV. The transition region between the ECAL barrel and endcap (1.444 < |η| < 1.566) is excluded in the reconstruction of electrons and the missing acceptance is corrected to the full |η| < 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 TeV centre-of-mass energy in our analysis. Electrons (muons) are considered isolated if the scalar p T sum of the nearby particles with R < 0.3 (0.4) is less than 15 (25)% of the lepton 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 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τ > 1 cm), neutrinos excluded, clustered with the anti-k T algorithm [44] using a radius parameter of 0.4. The jet fourmomentum 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 η and p 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. Data are further corrected for differences between the simulation and the pileup events and for dijet, Z + jet, and γ + jet events, where the p T -balance is exploited. Since the p T balance in Z + 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, γ, or Z) back-to-back in the transverse plane (|∆φ − π| < 0.34) associated with a soft jet, for various values of ρ = p soft jet T /p ref T , running from 0.1 to 0.3, and is extrapolated to ρ → 0. All jets down to p T = 5 or 10 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 T > 20 GeV and |y| < 2.4 are used in this analysis. The measured cross sections are defined with a 30 GeV threshold on the jet p T .

Backgrounds estimation
The contributions from background processes are estimated using the simulation samples described in Section 5 and are subtracted from the measured distributions. The dominant background, tt, is also measured from data. This tt background contributes mainly due to events with two same-flavour leptons. The production cross sections for e + e − and µ + µ − events from tt are identical to the cross section of e + µ − and e − µ + and can therefore be estimated from the latter. We select events in the tt 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 tt 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 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 in the tt control sample agree with data. The tt contribution to the signal is then obtained from 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 tt production. Multijet and 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 W events is estimated from simulation. Both contributions are found to be negligible. Fig. 2 shows the p bal T distribution separately for electron and muon channels. The tt background does not peak at the same p T balance as the signal, and has a broader spectrum. The JZB distribution is shown in Fig. 3. The tt background is asymmetric, making a larger contribution to the positive side of the distribution because transverse energy is carried away by neutrinos from W boson decays, leading to a reduction in the negative term of the JZB expression. Overall the agreement between data and simulation before the background subtraction is good and differences are within about 10%. ee data  The background distributions are obtained from the simulation, except for the tt 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 Section 5 has been used for the simulation.  Figure 2: Reconstructed data, simulated signal, and background distributions of the transverse momentum balance between the 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 tt 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 Section 5 has been used for the simulation.  : Reconstructed data, simulated signal, and background distributions of the JZB variable for the electron (left) and muon (right) channels. The background distributions are obtained from the simulation, except for the tt 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 Section 5 has been used for the simulation.

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 Z + jets sample simulated with MG5 aMC interfaced with PYTHIA8, using the FxFx merging scheme as described in Section 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.
The particle-level values refer to the stable leptons from the decay of the Z boson and to the jets built from the stable particles (cτ > 1 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 Z boson is taken to be the sum of the momenta of the two highest-p T electrons (or muons). The phase space for the cross section measurement is restricted to events with a Z boson mass between 71 and 111 GeV and both leptons with p T > 20 GeV and |η| < 2.4. Jets are required to have p T > 30 GeV, |y| < 2.4 and a spatial separation from the dressed leptons of R > 0.4.

Systematic uncertainties
The systematic uncertainties are propagated to the measurement by varying the corresponding source by one standard deviation up and down. The uncertainty sources are independent, and the resulting uncertainties are therefore added in quadrature. Tables 3 to 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 tt 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, α 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 mea-  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.

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]. 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 Section 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 ("≤4j LO" for up to four partons at LO and "≤2j 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 Z + 1 jet calculation is denoted as N jetti NNLO in the legends. The measured cross section values along with the uncertainties discussed in Section 9 are given in Tables 3 to 20.  (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 Z(→ + − )+ ≥ 0 jet, m + − > 50 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 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 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 T , below the peak, is NLO MG5 aMC, showing a better agreement than the NNLL τ ' calculation from GENEVA. The shape of the distribution in the region below 10 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 T region, populated with Z+ ≥ 1 jet events, 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 Z+ ≥ 0 jet cross section to its NNLO value.
The jet transverse momenta for the 1 st , 2 nd and 3 rd leading jets can be seen in Figs. 6 and 7 (Tables 6-8). The LO MG5 aMC predicted spectrum differs from the measurement, showing a steeper slope in the low p T region. The same feature was observed in the previous measurements [3,4]. The comparison with NLO MG5 aMC and N jetti NNLO calculation shows that adding NLO terms cures this discrepancy. The GENEVA prediction shows good agreement for the measured p T of the first jet, while it undershoots the data at low p 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-11). The NLO MG5 aMC, GENEVA, and NNLO1j predictions show very good agreement with data.
The total jet activity has been measured via the H 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 T < 400 GeV. For higher jet multiplicities both LO and NLO MG5 aMC are compatible with the measurement, although the contribution in the region H T < 400 GeV is smaller for LO than for NLO MG5 aMC. The description by NLO MG5 aMC is very good starting from 50 GeV. The contribution at lower values of H 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, the comparison is also done with N jetti NNLO and this calculation provides the best description. The uncertainty for N jetti NNLO is larger than in the jet transverse momentum distribution because of the contribution from the additional jets.
The balance in transverse momentum between the jets and the Z boson, p bal T , is shown in Figs. 12 and 13 (Tables 15-17) for inclusive jet multiplicities of 1, 2, and 3. When more jets are  Figure 4: Measured cross section for Z + 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 (α S ). The complete set was computed for one of the predictions. These uncertainties were added together in quadrature (represented by the ⊕ sign in the legend). included, the peak of p bal T 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 T > 30 GeV and |y| > 2.4. An imbalance in the event, i.e. p bal T 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.
The 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 T (Z) below and above 50 GeV. As expected in the high-p T (Z) region, i.e. in the high jet multiplicity sample, the distribution is more symmetric. The NLO MG5 aMC prediction provides a good description of the JZB distribution, while both GENEVA and LO MG5 aMC predictions do not. This applies to both configurations, 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.

Summary
We have measured differential cross sections for the production of a Z boson in association with jets, where the Z boson decays into two charged leptons with p T > 20 GeV and |η| < 2.4. The data sample corresponds to an integrated luminosity of 2.19 fb −1 collected with the CMS detector during the 2015 proton-proton LHC run at a centre-of-mass energy of 13 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 Z boson, jet kinematic variables including jet transverse momentum (p T ), the scalar sum of jet transverse momenta (H 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 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 T using the JZB variable. Jets with p T > 30 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 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 ≈ 10% measurement and 10% theoretical uncertainties, although this part of the calculation is only LO. The multi-     parton 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 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 have been compared with the fixed order calculation. The agreement is good, even for the H 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 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 Z + jets contribution at the LHC in measurements and searches, and its associated uncertainty.