Measurements of differential production cross sections for a Z boson in association with jets in pp collisions at sqrt(s) = 8 TeV

Cross sections for the production of a Z boson in association with jets in proton-proton collisions at a centre-of-mass energy of sqrt(s) = 8 TeV are measured using a data sample collected by the CMS experiment at the LHC corresponding to 19.6 inverse femtobarns. Differential cross sections are presented as functions of up to three observables that describe the jet kinematics and the jet activity. Correlations between the azimuthal directions and the rapidities of the jets and the Z boson are studied in detail. The predictions of a number of multileg generators with leading or next-to-leading order accuracy are compared with the measurements. The comparison shows the importance of including multi-parton contributions in the matrix elements and the improvement in the predictions when next-to-leading order terms are included.


Introduction
The high centre-of-mass energy of proton-proton (pp) collisions at the CERN LHC produces events with large jet transverse momenta (p T ) and high jet multiplicities in association with a Z/γ * boson. For convenience Z/γ * is denoted as Z. The selection of events in which the Z boson decays into two oppositely charged electrons or muons provides a signal sample that is not significantly contaminated by background processes. This decay channel can be reconstructed with high efficiency due to the presence of charged leptons in the final state and is well suited for the validation of calculations within the framework of perturbative quantum chromodynamics (QCD). Furthermore, the production of massive vector bosons with jets is an important background to a number of standard model (SM) processes (single top and tt production, vector boson fusion, WW scattering, Higgs boson production) as well as searches for physics beyond the SM. A good understanding of this background is vital to these searches and measurements. Perturbative QCD calculations of the differential cross sections involve different powers of the strong coupling constant α s and different kinematic scales and are therefore technically challenging. The issue has been addressed over the last 15 years by merging processes with different parton multiplicities before the parton showering, initially at tree level, and more recently with matrix elements calculated at next-to-leading order (NLO) using multileg matrix-element (ME) event generators [1,2].
In this paper we present measurements of the differential cross sections for Z boson production in association with jets at √ s = 8 TeV, in the electron and muon decay channels, using a data sample corresponding to an integrated luminosity of 19.6 fb −1 . Our measurements are compared with calculations obtained from different multileg ME event generators with leading order (LO) MEs (tree level), NLO MEs and a combination of NLO and LO MEs. Measurements of the Z + jets cross section were previously reported by the CDF and D0 Collaborations in proton-antiproton (pp) collisions at a centre-of-mass energy √ s = 1.96 TeV [3,4]. More recent results from proton-proton (pp) collisions at √ s = 7 TeV were published by the ATLAS [5,6] and CMS [7,8] Collaborations.
The cross sections are restricted to the phase space where the lepton transverse momenta are greater than 20 GeV, their absolute pseudorapidities are less than 2.4, and the dilepton mass is in the interval 91 ± 20 GeV. The jets are defined using the infrared and collinear safe anti-k t algorithm applied to all visible particles; the radius parameter is set to 0.5 [9]. The four-momenta of the particles are summed and therefore the jets can be massive. The differential cross sections include only those jets with transverse momentum greater than 30 GeV and further than R = 0.5 from the leptons in the (η, φ)-plane, where φ is the azimuthal angle. In addition, the absolute jet rapidity is required to be smaller than 2.4. The jets are referred to as 1 st , 2 nd , 3 rd , etc. according to their transverse momenta, starting with the highest-p T jet, and denoted as j 1 , j 2 , j 3 , etc. To further investigate the QCD dynamics towards low Bjorken-x values, multidimensional differential cross sections are measured for Z + ≥ 1 jet production in an extended phase space with jet rapidities up to 4.7. The extension of the rapidity coverage from 2.4 to 4.7 is used to tag events from vector boson fusion (e.g., Higgs production). Typically, the Z + jets events constitute a background for such processes, and a good understanding of their production differential cross section including jets in the forward region is important.
For each jet multiplicity (N jets ) a number of measurements are made: the total cross section in the defined phase space, the differential cross sections as functions of the jet transverse momentum scalar sum H T , and the differential cross sections as a function of the individual jet kinematics (transverse momentum p T , and absolute rapidity |y|). For the leading jet a double differential cross section is measured as a function of its absolute rapidity and transverse mo-mentum. Correlations in the jet kinematics are studied with one-dimensional and multidimensional differential cross section measurements, via 1) the distributions in the azimuthal angles between the Z boson and the leading jet and between the two leading jets, and 2) the rapidity distributions of the Z boson and the leading jet. These two rapidities are used as variables of a three-dimensional differential cross section measurement together with the transverse momentum of the jet. The Lorentz boost along the beam axis introduces a large correlation amongst the Z boson and the jet rapidities. The two rapidities are combined to form a variable uncorrelated with the event boost along the beam axis, y diff = 0.5 |y(Z) − y(j i )| and a highly boost-dependent variable, y sum = 0.5 |y(Z) + y(j i )|. The cross section is measured separately as a function of each of these variables. The distribution of y diff is mostly sensitive to the parton scattering, while y sum is expected to be sensitive to the parton distribution functions (PDFs) of the proton.
The Drell-Yan process, where the Z boson can decay into a pair of neutrinos, is a background to searches for new phenomena, such as dark matter, supersymmetry and other theories beyond the SM that predict the presence of invisible particle(s) in the final state. It is particularly important when the Z boson has a large transverse momentum, leading to a large missing transverse energy. The azimuthal angle between the jets is a good handle to suppress backgrounds coming from QCD multijet events, while the H T variable can be used to select events with large jet activity. For such analyses it is important to have a good model of Z + jets production and therefore a good understanding of these observables. This motivates the measurement of the distributions of the azimuthal angles between the jets and between the Z boson and the jets for different thresholds applied to the Z boson transverse momentum, the H T variable, and the jet multiplicity. These angles can be measured with high precision, and thus provide an excellent avenue to test the accuracy of SM predictions [10]. The dijet mass is an essential observable in the selection of Higgs boson events produced by vector boson fusion and it is important to model well both this process and its backgrounds. This observable is measured in Z + jets events for the two leading jets. Section 2 describes the experimental setup and the data samples used for the measurements, while the object reconstruction and the event selection are presented in Section 3. Section 4 is dedicated to the subtraction of the background contribution and the correction of the detector response, and Section 5 to the estimation of the measurement uncertainties. Finally, the results are presented and compared to different theoretical predictions in Section 6 and summarised in Section 7.

The CMS detector, simulation, and data samples
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 covering the range |η| < 2.5 together with a calorimeter covering the range |η| < 3. The latter consists of a lead tungstate crystal electromagnetic calorimeter (ECAL), and a brass and scintillator hadron calorimeter (HCAL). The HCAL is complemented by an outer calorimeter placed outside the solenoid used to measure the tails of hadron showers. The pseudorapidity coverage is extended up to |η| = 5.2 by a forward hadron calorimeter built using radiationhard technology. Gas-ionization detectors exploiting three technologies, drift tubes, cathode strip chambers, and resistive-plate chambers, are embedded in the steel flux-return yoke outside the solenoid and constitute the muon system, used to identify and reconstruct muons over the range |η| < 2.4. A more detailed description of the CMS detector, together with a definition of the coordinate system used and the relevant kinematic variables, can be found in Ref. [11].
Simulated events are used to both subtract the contribution from background processes and to correct for the detector response. The signal and the background (from WW, WZ, ZZ, tt, and single top quark processes) are modelled with the tree-level matrix element event generator MADGRAPH v5.1.3.30 [12] interfaced with PYTHIA 6.4.26 [13]. The PDF CTEQ6L1 [14] and the Z2* PYTHIA 6 tune [15,16] are used. For the ME calculations, α s is set to 0.130 at the Z boson mass scale. The five processes pp → Z + N jets jets, N jets = 0, . . . , 4, are included in the ME calculations. The k t −MLM [17,18] scheme with the merging scale set to 20 GeV is used for the matching of the parton showering (PS) with the ME calculations. The same setup is used to estimate the background from Z + jets → τ + τ − + jets. The signal sample is normalised to the inclusive cross section calculated at next-to-next-to-leading order (NNLO) with FEWZ 2.0 [19] using the CTEQ6M PDF set [14]. Samples of WW, WZ, ZZ events are normalised to the inclusive cross section calculated at NLO using the MCFM 6.6 [20] generator. Finally, an NNLO plus next-to-next-to-leading log (NNLO + NNLL) calculation [21] is used for the normalisation of the tt sample. When comparing the measurements with the predictions from theory, several other event generators are used for the Drell-Yan process. Those, which are not used for the measurement itself, are described in Section 6.
The detector response is simulated with GEANT4 [22]. The events reconstructed by the detector contain several superimposed proton-proton interactions, including one interaction with a high p T track that passes the trigger requirements. The majority of interactions, which do not pass trigger requirements, typically produce low energy (soft) particles because of the larger cross section for these soft events. The effect of this superposition of interactions is denoted as pileup. The samples of simulated events are generated with a distribution of the number of proton-proton interactions per beam bunch crossing close to the one 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 measured distribution of pileup interactions.

Event reconstruction and selection
Events with at least two leptons (electrons or muons) are selected. The trigger accepts events with two isolated electrons (muons) with a p T of at least 8 and 17 GeV. After reconstruction these leptons are restricted to a kinematic and geometric acceptance of p T > 20 GeV and |η| < 2.4. We require that the oppositely charged, same-flavor leptons form a pair with an invariant mass within a window of 91 ± 20 GeV. The ECAL barrel-endcap transition region 1.444 < |η| < 1.566 is excluded for the electrons. The acceptance is extended to the full |η| < 2.4 region when correcting for the detector response.
Information from all detectors is combined using the particle-flow (PF) algorithm [23,24] to produce an event consisting of reconstructed particle candidates. The PF candidates are then used to build jets and calculate lepton isolation. The quadratic sum of transverse momenta of the tracks associated to the reconstructed vertices is used to select the primary vertex (PV) of interest. Because pileup involves typically soft particles, the PV with the highest sum is chosen.
The electrons are reconstructed with the algorithm described in Ref. [25]. Identification criteria based on the electromagnetic shower shape and the energy sharing between ECAL and HCAL are applied. The momentum is estimated by combining the energy measurement in the ECAL with the momentum measurement in the tracker. For each electron candidate, an isolation variable, quantifying the energy flow in the vicinity of its trajectory, is built by summing the transverse momenta of the PF candidates within a cone of size ∆R = √ (∆φ) 2 + (∆η) 2 = 0.3, excluding the electron itself and the charged particles not compatible with the primary event vertex. This sum is affected by neutral particles from pileup events, which cannot be rejected with a vertex criterion. An average energy density per unit of ∆R is calculated event by event using the method introduced in Ref. [26] and used to estimate and subtract the neutral particle contribution. The electron is considered isolated if the isolation variable value is less than 15% of the transverse momentum of the electron. The electron candidates are required to be consistent with a particle originating from the PV in the event.
Muon candidates are matched to tracks measured in the tracker, and they are required to satisfy the identification and quality criteria described in Ref. [27] that are based on the number of tracker hits and the response of the muon detectors. The background from cosmic ray muons, which appear as two back-to-back muons, is reduced by criteria on the impact parameter and by requiring that the muon pairs have an acollinearity larger than 2.5 mrad. An isolation variable is defined that is similar to that for electrons, but with a cone size ∆R = 0.4 and a different approach to the subtraction of the contribution from neutral pileup particles. For the muons this contribution is estimated from the sum of the transverse momenta of the charged particles rejected by the vertex requirement, considered as coming from pileup. This sum is multiplied by a factor of 0.5 to take into account the relative fraction of neutral and charged particles. A muon is considered isolated if the isolation variable value is below 20% of its transverse momentum.
The efficiencies for the lepton trigger, reconstruction, and identification are measured with the "tag-and-probe" method [28]. The simulation is corrected using the ratios of the efficiencies obtained in the data sample to those obtained in the simulated sample. These scale factors for lepton reconstruction and identification typically range from 0.95 to 1.05 depending on the lepton transverse momentum and rapidity. The overall efficiency of trigger and event selection is 58% for the electron channel and 88% for the muon channel.
The anti-k T algorithm, with a radius parameter of 0.5, is used to cluster PF candidates to form hadronic jets. The jet momentum is determined as the vectorial sum of all particle momenta. Charged hadrons identified as coming from a pileup event vertex are rejected from the jet clustering. The remaining contribution from pileup events, which comes from neutral hadrons and from charged hadrons whose PV has not been unambiguously identified, is estimated and subtracted event-by-event using a technique based on the jet area method [9,29]. Jet energy corrections are derived from the simulation, and are confirmed with in situ measurements using the energy balance of dijet and photon+jet events [30]. Jets with a transverse momentum less than 30 GeV or overlapping within ∆R = 0.5 with either of the two leptons from the decay of the Z boson are discarded. The single differential cross sections are measured for jet rapidity within |y| < 2.4, which is the region with the best jet resolution and pileup rejection. In this measurement, the jet multiplicity refers to the number of jets fulfilling the jet criteria, within the |y| < 2.4 boundary for the one-dimensional differential cross sections. For the multidimensional differential cross sections reported in Section 6.8 the region for the jet rapidity is extended to |y| < 4.7.

Background subtraction and correction for the detector response
In Fig. 1 the event yield in the electron and muon channels is compared to the simulation. The agreement between simulation and data is excellent up to four jets. Since the Z + jets simulation does not include more than four partons in the ME calculations, we expect a less accurate prediction of the signal for jet multiplicities above four. Background contamination is below 1% for a jet multiplicity of one and increases with the jet multiplicity. The background represents 2% of the event yield for a jet multiplicity of 2 and 20% for a jet multiplicity of 5. The background contribution is estimated from the samples of simulated events described in Section 2 and subtracted bin-by-bin from the data. The simulations are validated using a µ ± e ∓ data control sample, as explained in Section 5. The background contribution from multijet events where the jets are misidentified as leptons is checked with a lepton control sample using two leptons with the same flavor and charge and found to be negligible.
Unfolding the detector response corrects the signal distribution for the migration of events between closely separated bins and across boundaries of the fiducial region. The unfolding procedure also includes a correction for the efficiency of the trigger, and the lepton reconstruction and identification. The unfolding procedure is applied separately to each measured differential cross section. In a first step, the data distribution is corrected to remove the background contribution and the contribution from signal process events outside of the defined phase space. Then, the iterative D'Agostini method [31], as implemented in the statistical analysis toolkit ROOUNFOLD [32], is used to correct for bin-to-bin migration and for efficiency. Using the simulation the method generates a response matrix that relates the probability that an event in bin i of the differential cross section is reconstructed in bin j. These probabilities include the case of bin-i events that do not pass the selection criteria on the reconstructed event or fall outside the distribution boundaries. For the three-dimensional differential cross section, the method is applied within each (y(j 1 ), y(Z)) bin, where the unfolding is performed with respect to the p T (j 1 ) observable. Similarly, the unfolding of the double differential cross sections is performed with respect to the most sensitive variable: p T (j 1 ) for d 2 σ/dp T (j 1 )dy(j 1 ) and y(j 1 ) for d 2 σ/dy(j 1 )dy(Z).

Systematic uncertainties
The response matrices are built from reconstructed and generated quantities using the MAD-GRAPH 5 + PYTHIA 6 Z + jets simulation sample. The generated values refer to the leptons from the decay of the Z boson and to the jets built from the stable particles 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 final-state radiation, and the lepton is said to be "dressed". Although this process does not recover all the final-state radiation; it removes most differences between electrons and muons, and the dilepton mass spectra are identical for the two decay channels after this procedure, as checked with the MADGRAPH 5 + PYTHIA 6 simulation. The Z boson is reconstructed from the dressed lepton momentum vectors.
The phase space for the cross section measurement is restricted to p T > 20 GeV and |η| < 2.4 for both dressed charged leptons, a dilepton mass within 91 ± 20 GeV, and the jet kinematics constrained to p T > 30 GeV and |y| < 2.4. For the extended multidimensional cross sections the phase space is extended to |y| < 4.7. The measured cross section values include the Z branching fraction to a single lepton flavour.
The rejection of jets originating from pileup is more difficult outside the tracker geometric acceptance, since the vertex constraint cannot be used to reject the charged particles coming from pileup. Consequently, despite the jet pileup rejection criterion, a contamination of jets from pileup remains and needs to be subtracted. This region beyond the tracker acceptance, 2.5 < |y| < 4.7, is used only for the Z + ≥ 1 jet multidimensional differential cross section measurements, where only the leading jet is relevant. The fraction of events in which a jet comes from pileup, denoted as f PU , is estimated using a control sample of a Z boson associated with one jet, obtained by requiring one jet with p T above 30 GeV and a veto condition of no other jets with p T > 12 GeV and above 20% of the Z boson transverse momentum. Since a Z boson and a jet coming from two different pp collisions are independent, the distribution of ∆φ(Z, j) is expected to be flat, which is confirmed by the simulation. For the Z boson and the jet from a pp → Z + 1 jet event the distribution is expected to peak at π. The constraint on additional jets enforces the p T balance between the Z boson and the jet, reducing the contribution to low values of ∆φ(Z, j). The simulation shows that this contribution is negligible in the region ∆φ(Z, j) < 1. Therefore, f PU is estimated from the fraction of events in that region. The value of f PU is 30% in the most forward part (3.2 < |y| < 4.7) and lowest p T measurement bin (30 GeV to 40 GeV). The fraction of pileup events estimated from the pileup control sample is used to correct the signal data sample. The same method is applied to the simulation. The ratio of the value of f PU obtained from the simulation to that measured in data decreases monotonically as a function of p T . In the p T bin 30-50 GeV it ranges up to 1.25 (1.35) for |y(j 1 )| between 2.5 and 3.2 (3.2 and 4.7). Beyond p T = 50 GeV the discrepancy is negligible and the results are identical with or without the pileup subtraction.
Since the |y(j)| < 2.4 region contains the bulk of the events and including the |y(j)| > 2.4 region does not improve the precision of the measurement, most of the differential cross section measurements are limited to |y(j)| < 2.4. The subtraction of pileup contributions is not needed when confining measurements to this region.

Systematic uncertainties
The systematic uncertainty in the background-subtracted data distributions is estimated by varying independently each of the contributing factors before the unfolding, and computing the difference induced by the variation in the unfolded distributions. The observed difference between the positive and negative uncertainties is small and the two are therefore averaged.
The different sources of uncertainties are independent and are added in quadrature. The unfolded histogram can be written as a linear combination of the bin contents of the backgroundsubtracted data histogram [31]. This linear combination is used to propagate analytically the statistical uncertainties to the unfolded results and calculate the full covariance matrix for each distribution, separately for each Z boson decay channel. The dominant source of systematic uncertainties is the jet energy correction. The various contributions are listed in Table 1.
The jet energy correction uncertainty (JEC in the table) is calculated by varying this correction by one standard deviation. This uncertainty is p T -and η-dependent and varies from 1.5% up to 5% for |η| < 2.5 and from 7% to 30% for |η| > 2.5. The uncertainty in the measured cross section is between 5.3% and 28% depending on the jet multiplicity. The jet energy resolution (JER) uncertainty is estimated for data and simulation in Ref. [33]. The resulting uncertainty in the measurement is below 1% for all the multiplicities. Other significant background contributions come from tt, diboson, and Z → τ + τ − processes. The related uncertainty (Bkg) is estimated by varying the cross section for each of the background processes (tt, ZZ, WZ, and WW) independently by 10% for tt and 6% for diboson processes. The normalisation variation for the tt events is chosen to cover the maximum observed difference between the simulation and the data in the jet multiplicity, transverse momentum, and rapidity distributions when selecting events with two leptons of oppositely charged, different flavours (µ ± e ∓ ). The uncertainty in diboson cross sections covers theoretical and PDF contributions. The resulting uncertainty in the measurement increases with the jet multiplicity and reaches 4.3%.
Another source of uncertainty is the modelling of the pileup (PU). The number of interactions per bunch crossing in simulated samples is varied by 5%. This covers effects related to the modeling of simulated minimum bias events of 3%, the estimate of the number of interactions per bunch crossing in data based on luminosity measurements of 2.6%, and the experimental uncertainties entering inelastic cross section measurements of 2.9%. The resulting uncertainties range from 0.26% to 5.6% depending on the jet multiplicity. The uncertainty from the pileup subtraction performed in the forward region, |y| > 2.5, is estimated by varying up and down the pileup fraction f PU described in Section 4 by half the difference from the value obtained in the simulation. In the region covered by the tracker and where no correction is applied, it is verified that the jet multiplicity does not depend on the number of vertices reconstructed in the event. This indicates that the jets from pileup events have a negligible impact on the measurement.
The unfolding procedure has an uncertainty due to its dependence on the simulation used to estimate the response matrix (Unf sys) and to the finite size of the simulation sample (Unf stat). The first contribution is estimated using an alternative event generator, SHERPA 1.4 [34], and taking the difference between the two results to represent the uncertainty. The distribution obtained with the alternative generator differs sufficiently from the nominal one to cover the differences with the data. The statistical uncertainty in the response matrix is analytically propagated to the unfolded result [31]. When added in quadrature and depending on the kinematic variable and jet multiplicity, the total unfolding uncertainty varies up to 10%.
The uncertainty in the efficiency of the lepton reconstruction, identification, and isolation is propagated to the measurement by varying the total data-to-simulation scale factor by one standard deviation. It amounts to 2.5% and 2.6% in the dimuon and dielectron channels, respectively.
The uncertainty in the integrated luminosity amounts to 2.6% [35]. Since the background event yield normalisation also depends on the integrated luminosity, the effect of the above uncertainty on the background yield (Lumi) can be larger and amounts to 3.9% in the bins with low signal purity.

Results
The measurements from the electron and muon channels are consistent and are combined using a weighted average. 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 uncertainty sources except for statistical uncertainties and those associated with lepton reconstruction and identification, which are taken to be uncorrelated. The measured differential cross sections are compared to the results obtained from three different calculations as described below.

Theoretical predictions
The measurements are compared to a tree level calculation and two multileg NLO calculations. The first prediction is computed with MADGRAPH 5 [12] interfaced with PYTHIA 6 (denoted as MG5 + PY6 in the figure legends), for parton showering and hadronisation, with the configuration described in Section 2. The total cross section is normalised to the NNLO cross section computed with FEWZ 2.0 [19]. Two multileg NLO predictions including parton showers using the MC@NLO [36] method are used. For these two predictions the total cross section is normalised to the one obtained with the respective event generators. The total cross section values used for the normalisation are summarised in Table 2. the Z boson production processes with 0, 1, and 2 partons at NLO. The FxFx [2] merging scheme is used with a merging scale parameter set to 30 GeV. The NNPDF 3.0 NLO PDF [42] is used for the ME calculations while the NNPDF 2.3 QCD + QED LO [43,44] is used for the backward evolution of the showering. For the ME calculations, α s is set to the current PDG world average [45] rounded to α s (m Z ) = 0.118. For the showering and underlying events the value of the CUETP8M1 tune, α s (m Z ) = 0.130, is used. The larger value is expected to compensate for the missing higher order corrections. NLO accuracy is achieved for pp → Z + N jets jets, N jets = 0, 1, 2 and LO accuracy for N jets = 3. For this prediction, theoretical uncertainties are computed and include the contribution from the fixed-order calculation and from the NNPDF 3.0 PDF set. The two uncertainties are added in quadrature. The fixed-order calculation uncertainties are estimated by varying the renormalisation and the factorisation scales by factors of 1/2 and 2. The envelope of the variations of all factor combinations, excluding the two combinations when one scale is varied by a factor 1/2 and the other by a factor 2, is taken as the uncertainty. The reweighting method [46] provided by the MG5 aMC generator is used to derive the cross sections with the different renormalisation and factorisation scales and with the different PDF replicas used in the PDF uncertainty determination. For the NLO predictions, weighted samples are used (limited to ±1 weights in the case of aMC@NLO), which can lead to larger statistical fluctuations than expected for unweighted samples in some bins of the histograms presented in this section.

Jet multiplicity
The cross sections for jet multiplicities from 0 to 7, and the comparisons with various predictions are presented in this section. Figure 2 shows the cross section for both inclusive and exclusive jet multiplicities and the numbers are compared with the prediction obtained with MG5 aMC + PYTHIA 8 in Table 3 for the exclusive case. The agreement with the predictions is very good for jet multiplicities up to the maximum number of final-state partons included in the ME calculations, namely three for MG5 aMC + PYTHIA 8 and four for both MADGRAPH 5 + PYTHIA 6 and SHERPA 2. The level of precision of the measurement does not allow us to probe the improvement expected from the additional NLO terms. The cross section is reduced by a factor of five for each additional jet.
The predictions already agree well at tree level (MADGRAPH 5 + PYTHIA 6) renormalised to the NNLO inclusive cross section. For N jets = 4, the MG5 aMC + PYTHIA 8 calculation, which does not include this jet multiplicity in the matrix elements, predicts a different cross section from those that do. The predictions that include four jets in the matrix elements are in better agreement with the data, but the difference between the predictions is limited to roughly one standard deviation of the measurement uncertainty. The large uncertainty is due to the sensitivity of the jet p T threshold acceptance to the jet energy scale. The SHERPA 2 prediction for N jets = 5 is closer to the measurement than MADGRAPH 5 + PYTHIA 6, while neither of these in-cludes this multiplicity in the ME calculations. The theoretical uncertainty shown in the figure for the MG5 aMC + PYTHIA 8 prediction uses the standard method described in the previous subsection. In the case of the exclusive jet multiplicity, the presence of large logarithms in the perturbative calculation can lead to an underestimate of this uncertainty, so the Steward and Tackmann prescription (ST) provides a better estimate [47]. The uncertainties calculated with both prescriptions are provided in Table 3. For the calculations considered here, the increase of the ST uncertainty with respect to the standard one is moderate. This is consistent with the observation that the agreement with the measurement and the coverage of the difference by the theoretical uncertainty in Fig. 2 is similar for the inclusive and exclusive jet multiplicities.

Jet transverse momentum
Knowledge of the kinematics of SM events with large jet multiplicity is essential for the LHC experiments since these events are backgrounds to searches for new physics that predict decay chains of heavy coloured particles, such as squarks, gluinos, or heavy top quark partners. The measured differential cross sections as a function of jet p T for the 1 st , 2 nd , 3 rd , 4 th , and 5 th jets are presented in Figs. 3-5. The cross sections fall rapidly with increasing p T . The cross section for the leading jet is measured for p T values between 30 GeV and 1 TeV and decreases by more than five orders of magnitude over this range. The cross section for the fifth jet is measured for p T values between 30 and 100 GeV and decreases even faster, mainly because of the phase space covered.
For the leading jet, the agreement of the MADGRAPH 5 + PYTHIA 6 prediction with the measurement is very good up to ≈150 GeV. Discrepancies are observed from ≈150 to ≈450 GeV. A similar excess in the ratio with the tree-level calculation was observed at √ s = 7 TeV in the CMS measurement [8], using predictions from the same generators, as well as in the ATLAS measurement [5], which used ALPGEN [48] interfaced to HERWIG [49] for the predictions. The calculations that include NLO terms for this jet multiplicity do not show this discrepancy. The theo. ⊕ unc. s α ⊕ PDF ⊕ Figure 2: The cross section for Z (→ ) + jets production measured as a function of the (left) exclusive and (right) inclusive jet multiplicity distributions compared to the predictions calculated with MADGRAPH 5 + PYTHIA 6, SHERPA 2, and MG5 aMC + PYTHIA 8. The lower panels show the ratios of the theoretical predictions to the measurements. Error bars around the experimental points show the statistical uncertainty, while the cross-hatched bands indicate the statistical and systematic uncertainties added in quadrature. The boxes around the MG5 aMC + PYTHIA 8 to measurement ratio represent the uncertainty on the prediction, including statistical, theoretical (from scale variations), and PDF uncertainties. The dark green area represents the statistical and theoretical uncertainties only, while the light green area represents the statistical uncertainty alone.
prediction from SHERPA 2 shows some disagreement with data in the low transverse momentum region. The second jet shows similar behaviour. Both MADGRAPH 5 + PYTHIA 6 and MG5 aMC + PYTHIA 8 are in good agreement with the measurement for the third jet p T spectrum. The shape predicted by the calculations from SHERPA 2 differs from the measurement since the predicted spectrum is harder. For the 4 th jet, the three predictions agree well with the measurements. Calculations from SHERPA 2 and MG5 aMC + PYTHIA 8 predict different spectra. Based on the experimental uncertainties it is difficult to arbitrate between the two, although we expect the one that includes four partons in the matrix elements to be more accurate. The agreement of SHERPA 2 and MADGRAPH 5 + PYTHIA 6 calculations with the measured 5 th jet p T spectrum is similar.
In summary, including many jet multiplicities in the matrix elements provides a good description of the different jet transverse momentum spectra. Including NLO terms improves the agreement with the measured spectra. Nevertheless, some differences are observed between the predictions calculated with SHERPA 2 and MG5 aMC + PYTHIA 8. The two calculations differ in many ways, other than the fixed order: different PDF choices, different jet merging schemes, and different showering models. In Ref. [8] it was shown that the jet p T spectra have little dependence on the PDF choice, therefore the difference between the two generator is likely to be due to the different parton showerings or jet merging schemes.    The differential cross section for Z (→ ) + jets production measured as a function of the 5 th jet p T compared to the predictions calculated with MADGRAPH 5 + PYTHIA 6, SHERPA 2, and MG5 aMC + PYTHIA 8. The lower panels show the ratios of the theoretical predictions to the measurements. Error bars around the experimental points show the statistical uncertainty, while the cross-hatched bands indicate the statistical and systematic uncertainties added in quadrature. The boxes around the MG5 aMC + PYTHIA 8 to measurement ratio represent the uncertainty on the prediction, including statistical, theoretical (from scale variations), and PDF uncertainties. The dark green area represents the statistical and theoretical uncertainties only, while the light green area represents the statistical uncertainty alone.

Jet and Z boson rapidity
The differential cross sections as a function of the absolute rapidity of the first, second, third, fourth, and fifth jets are presented in Figs. 6, 7, and 8, including all events with at least one, two, three, four, and five jets. The differential cross sections in |y| have similar shapes for all jets while they vary by about a factor 2 in the range from 0 to 2.4.
The predictions obtained with SHERPA 2 provide the best overall description regarding the shape of data distributions. The predictions of both MADGRAPH 5 + PYTHIA 6 and MG5 aMC + PYTHIA 8 have a more central distribution than is measured for jets 1 to 4, although this behaviour is less pronounced for the latter. The difference could be attributed to the different showering methods and the different PDF choices for the three predictions. Given the experimental uncertainties, the shape of the spectrum of the 5 th jet rapidity is equally well described by the three calculations.   theo. ⊕ unc. s α ⊕ PDF ⊕ Figure 6: The differential cross section for Z (→ ) + jets production measured as a function of the (left) 1 st and (right) 2 nd jet |y| compared to the predictions calculated with MADGRAPH 5 + PYTHIA 6, SHERPA 2, and MG5 aMC + PYTHIA 8. The lower panels show the ratios of the theoretical predictions to the measurements. Error bars around the experimental points show the statistical uncertainty, while the cross-hatched bands indicate the statistical and systematic uncertainties added in quadrature. The boxes around the MG5 aMC + PYTHIA 8 to measurement ratio represent the uncertainty on the prediction, including statistical, theoretical (from scale variations), and PDF uncertainties. The dark green area represents the statistical and theoretical uncertainties only, while the light green area represents the statistical uncertainty alone.
The Z boson rapidity distribution is presented in Fig. 9 with no requirement on the Z boson transverse momentum. To minimize the uncertainties the measurement is done for the normalized distributions. The relative contributions of matrix elements and parton shower depend on the Z transverse momentum. The measurement is also performed with a lower limit of 150 and theo. ⊕ unc. s α ⊕ PDF ⊕ Figure 8: The differential cross section for Z (→ ) + jets production measured as a function of the 5 th jet |y| compared to the predictions calculated with MADGRAPH 5 + PYTHIA 6, SHERPA 2, and MG5 aMC + PYTHIA 8. The lower panels show the ratios of the theoretical predictions to the measurements. Error bars around the experimental points show the statistical uncertainty, while the cross-hatched bands indicate the statistical and systematic uncertainties added in quadrature. The boxes around the MG5 aMC + PYTHIA 8 to measurement ratio represent the uncertainty on the prediction, including statistical, theoretical (from scale variations), and PDF uncertainties. The dark green area represents the statistical and theoretical uncertainties only, while the light green area represents the statistical uncertainty alone. 300 GeV on the Z boson transverse momentum. Each distribution is normalised to unity. The three calculations are in very good agreement with the measured values. The agreement of the prediction calculated with SHERPA 2 degrades when applying a threshold on the Z boson p T , though it is still consistent with data within the statistical uncertainty.  Figure 9: The normalised differential cross section for Z (→ ) + jets production measured as a function of Z boson rapidity compared to the predictions calculated with MADGRAPH 5 + PYTHIA 6, SHERPA 2, and MG5 aMC + PYTHIA 8. The cross section is measured (left) inclusively with respect to the Z boson p T , (middle) for p T > 150 GeV, and (right) for p T > 300 GeV. The lower panels show the ratios of the theoretical predictions to the measurements. Error bars around the experimental points show the statistical uncertainty, while the cross-hatched bands indicate the statistical and systematic uncertainties added in quadrature. The boxes around the MG5 aMC + PYTHIA 8 to measurement ratio represent the uncertainty on the prediction, including statistical, theoretical (from scale variations), and PDF uncertainties. The dark green area represents the statistical and theoretical uncertainties only, while the light green area represents the statistical uncertainty alone.
The correlations in rapidity between the different objects (Z boson and jets) are shown in Figs. 10 to 14. The normalised cross section is presented as a function of the rapidity difference between the Z boson and the leading jet, y diff (Z, j 1 ) = 0.5|y(Z) − y(j 1 )| in Fig. 10. A large discrepancy is observed between the measured cross section and that predicted by MADGRAPH 5 + PYTHIA 6. Such an effect was previously observed at √ s = 7 TeV [50] and is confirmed here with an increased statistical precision and with an extended range in y diff (Z, j 1 ). The discrepancy is significantly reduced when a threshold is applied to the transverse momentum of the Z boson as shown in the same figure. This observation supports the attribution of the discrepancy to the matching procedure between the ME and PS, as discussed in [50]. By contrast, a quite good agreement is found, independently of any threshold on the Z boson transverse momentum, for the NLO predictions of SHERPA 2, and MG5 aMC + PYTHIA 8. This improvement is expected to come from additional diagrams at NLO with a gluon propagator in the t-channel that populate the forward rapidity regions.
The presence of additional jets in the event should reduce the dependence on the ME/PS matching for the first jet since this jet will have a larger p T on average. Figure 11 shows the normalised cross section for Z production with at least two jets as a function of the rapidity difference between the Z boson and the leading jet, y diff (Z, j 1 ), between the Z boson and the second-leading jet, y diff (Z, j 2 ), and between the Z boson and the system formed by the two leading jets, y diff (Z, dijet). The discrepancies between the measured cross sections and the MADGRAPH 5 + PYTHIA 6 predictions are present in all three cases, but they are less pronounced than in the one-jet case (Fig. 11a compared to Fig. 10a). The NLO predictions from SHERPA 2 and MG5 aMC + PYTHIA 8 reproduce the measured dependencies much better than MADGRAPH 5 + PYTHIA 6 does.   theo. ⊕ unc. theo. ⊕ unc. s α ⊕ PDF ⊕ Figure 10: The normalised differential cross section for Z (→ ) + jets (N jets ≥ 1) production measured as a function of the y diff of the Z boson and the leading jet compared to the predictions calculated with MADGRAPH 5 + PYTHIA 6, SHERPA 2, and MG5 aMC + PYTHIA 8. (left) The cross section is measured inclusively with respect to the Z boson p T and for two different p T (Z) thresholds. The ratio of the prediction to the measurements is shown for (left) p T > 0 GeV, (middle) p T > 150 GeV, and (right) p T > 300 GeV. The lower panels show the ratios of the theoretical predictions to the measurements. Error bars around the experimental points show the statistical uncertainty, while the cross-hatched bands indicate the statistical and systematic uncertainties added in quadrature. The boxes around the MG5 aMC + PYTHIA 8 to measurement ratio represent the uncertainty on the prediction, including statistical, theoretical (from scale variations), and PDF uncertainties. The dark green area represents the statistical and theoretical uncertainties only, while the light green area represents the statistical uncertainty alone.
The rapidity correlation of the two leading jets, independently of the Z boson rapidity, is displayed in Fig. 12, showing the rapidity sum and rapidity difference between the two jets. There is a good agreement between the measured cross section and the three predictions for the rapidity sum dependence. The rapidity difference presents a discrepancy with MADGRAPH 5 + PYTHIA 6 at large values, while the NLO predictions of SHERPA 2, and MG5 aMC + PYTHIA 8 are in good agreement with the data.
The rapidity sum for the system of the Z boson and the leading jet is studied with different thresholds applied to the transverse momentum of the Z boson. Figure 13 shows the normalised cross section as a function of the rapidity sum of the Z boson and the leading jet, y sum (Z, j 1 ) = 0.5|y(Z) + y(j 1 )| for Z boson transverse momentum above 0, 150, and 300 GeV. The observed discrepancy between the measured cross section and that predicted by MAD-GRAPH 5 + PYTHIA 6 is similar to the effect that has been found at 7 TeV [50], and is confirmed here with increased statistical precision. The discrepancy almost vanishes when the transverse momentum of the Z boson is required to be larger than 150 GeV. The NLO predictions of SHERPA 2, and MG5 aMC + PYTHIA 8 are in good agreement with the measured cross section independently of the Z boson transverse momentum. This improvement with respect to MAD-GRAPH 5 + PYTHIA 6 can be attributed to either the different PDF choice, or to the NLO terms. theo. ⊕ unc. s α ⊕ PDF ⊕ Figure 11: The normalised differential cross section for Z (→ ) + jets (N jets ≥ 2) production measured as a function of the y diff of the Z boson and (left) the leading jet, (middle) the secondleading jet, and (right) the system constituted by these two jets. The measurement is compared to the predictions calculated with MADGRAPH 5 + PYTHIA 6, SHERPA 2, and MG5 aMC + PYTHIA 8. The lower panels show the ratios of the theoretical predictions to the measurements. Error bars around the experimental points show the statistical uncertainty, while the cross-hatched bands indicate the statistical and systematic uncertainties added in quadrature. The boxes around the MG5 aMC + PYTHIA 8 to measurement ratio represent the uncertainty on the prediction, including statistical, theoretical (from scale variations), and PDF uncertainties. The dark green area represents the statistical and theoretical uncertainties only, while the light green area represents the statistical uncertainty alone. theo. ⊕ unc. s α ⊕ PDF ⊕ Figure 12: The normalised differential cross section for Z (→ ) + jets (N jets ≥ 2) production measured as a function of the (left) y diff and (right) y sum of the two leading jets. The measurement is compared to the predictions calculated with MADGRAPH 5 + PYTHIA 6, SHERPA 2, and MG5 aMC + PYTHIA 8. The lower panels show the ratios of the theoretical predictions to the measurements. Error bars around the experimental points show the statistical uncertainty, while the cross-hatched bands indicate the statistical and systematic uncertainties added in quadrature. The boxes around the MG5 aMC + PYTHIA 8 to measurement ratio represent the uncertainty on the prediction, including statistical, theoretical (from scale variations), and PDF uncertainties. The dark green area represents the statistical and theoretical uncertainties only, while the light green area represents the statistical uncertainty alone.
For dijet events, Fig. 14 shows cross sections as a function of rapidity sums, for the Z boson and the leading jet, for the Z boson and the second-leading jet, and for the Z boson and the dijet system of the two leading jets. Comparison between the measured cross sections and the MADGRAPH 5 + PYTHIA 6 predictions exhibit a small disagreement for a rapidity sum above 1 for each jet, and the discrepancies increase when the dijet system is considered. Comparison with NLO predictions from SHERPA 2, and from MG5 aMC + PYTHIA 8 shows a very good agreement.
The rapidity correlation study confirms the observations made at √ s = 7 TeV, and shows that the behaviour with respect to the tree-level prediction is similar for the correlation with the second jet and enhanced when considering the dijet system consisting of the two leading jets. The study demonstrates that the two NLO predictions improve the agreement with the measurements, especially for the rapidity difference observables.    theo. ⊕ unc. s α ⊕ PDF ⊕ Figure 13: The normalised differential cross section for Z (→ ) + jets (N jets ≥ 1) production measured as a function of the y sum of the Z boson and the leading jet compared to the predictions calculated with MADGRAPH 5 + PYTHIA 6, SHERPA 2, and MG5 aMC + PYTHIA 8. The cross section is measured inclusively with respect to the Z boson p T and for two different p T (Z) thresholds. The ratio of the prediction to the measurements is shown for (left) p T > 0 GeV, (middle) p T > 150 GeV, and (right) p T > 300 GeV. The lower panels show the ratios of the theoretical predictions to the measurements. Error bars around the experimental points show the statistical uncertainty, while the cross-hatched bands indicate the statistical and systematic uncertainties added in quadrature. The boxes around the MG5 aMC + PYTHIA 8 to measurement ratio represent the uncertainty on the prediction, including statistical, theoretical (from scale variations), and PDF uncertainties. The dark green area represents the statistical and theoretical uncertainties only, while the light green area represents the statistical uncertainty alone.

Differential cross section in jet H T
The hadronic activity of an event can be probed with the scalar sum of the transverse momenta of the jets, H T . Measuring hadronic activity is important in searches for signatures with high jet activity or, by contrast, when wishing to veto such activity, for instance in the central region when looking for vector boson fusion induced processes. In this section we present measurements of the spectra for this variable in Z + jets events. The differential cross sections are shown in Figs. 15-17 for the different inclusive jet multiplicities.     theo. ⊕ unc. s α ⊕ PDF ⊕ Figure 14: The normalised differential cross section for Z (→ ) + jets (N jets ≥ 2) production measured as a function of the y sum of the Z boson and (left) the leading jet, (middle) the secondleading jet, and (right) the system constituted by these two jets. The measurement is compared to the predictions calculated with MADGRAPH 5 + PYTHIA 6, SHERPA 2, and MG5 aMC + PYTHIA 8. The lower panels show the ratios of the theoretical predictions to the measurements. Error bars around the experimental points show the statistical uncertainty, while the cross-hatched bands indicate the statistical and systematic uncertainties added in quadrature. The boxes around the MG5 aMC + PYTHIA 8 to measurement ratio represent the uncertainty on the prediction, including statistical, theoretical (from scale variations), and PDF uncertainties. The dark green area represents the statistical and theoretical uncertainties only, while the light green area represents the statistical uncertainty alone.
The predictions of the generators agree well with the measurements within the experimental uncertainties. For events with three or more jets (Figs. 16), all three simulations predict a distribution that falls more steeply at low values of H T . For the normalised distributions the total uncertainties in the measurements for Z+ ≥ 3 jets reduce for the three first bins to 21%, 10%, and 3.3%, respectively. This indicates that the difference in the shape is significant for MADGRAPH 5 + PYTHIA 6 and MG5 aMC + PYTHIA 8 predictions.  Figure 15: The differential cross section for Z (→ ) + jets production measured as a function of H T for (left) N jets ≥ 1 and (right) N jets ≥ 2 compared to the predictions calculated with MADGRAPH 5 + PYTHIA 6, SHERPA 2, and MG5 aMC + PYTHIA 8. The lower panels show the ratios of the theoretical predictions to the measurements. Error bars around the experimental points show the statistical uncertainty, while the cross-hatched bands indicate the statistical and systematic uncertainties added in quadrature. The boxes around the MG5 aMC + PYTHIA 8 to measurement ratio represent the uncertainty on the prediction, including statistical, theoretical (from scale variations), and PDF uncertainties. The dark green area represents the statistical and theoretical uncertainties only, while the light green area represents the statistical uncertainty alone. Figure 18 shows the differential cross section measurements as a function of the azimuthal angle between the Z boson and the leading jet for three different jet multiplicities. The inclusion of several parton multiplicities in the ME calculations ensures that the Monte Carlo predictions model the data well even at tree level and small differences in azimuthal angles. Differences are observed between tree-level (MADGRAPH 5 + PYTHIA 6) and multileg NLO (SHERPA 2, and MG5 aMC + PYTHIA 8) predictions, the latter being closer to the measurement, but the difference is smaller than one standard deviation in the experimental uncertainties. As the jet multiplicity increases, the ∆φ(Z, j 1 ) distribution flattens out. In an event dominated by the leading jet, the jet is recoiling against the Z boson, resulting in a strong peak at ∆φ(Z, j 1 ) π. theo. ⊕ unc. s α ⊕ PDF ⊕ Figure 16: The differential cross section for Z (→ ) + jets production measured as a function of H T for (left) N jets ≥ 3 and (right) N jets ≥ 4 compared to the predictions calculated with MADGRAPH 5 + PYTHIA 6, SHERPA 2, and MG5 aMC + PYTHIA 8. The lower panels show the ratios of the theoretical predictions to the measurements. Error bars around the experimental points show the statistical uncertainty, while the cross-hatched bands indicate the statistical and systematic uncertainties added in quadrature. The boxes around the MG5 aMC + PYTHIA 8 to measurement ratio represent the uncertainty on the prediction, including statistical, theoretical (from scale variations), and PDF uncertainties. The dark green area represents the statistical and theoretical uncertainties only, while the light green area represents the statistical uncertainty alone.  Figure 17: The differential cross section for Z (→ ) + jets production measured as a function of H T for N jets ≥ 5 compared to the predictions calculated with MADGRAPH 5 + PYTHIA 6, SHERPA 2, and MG5 aMC + PYTHIA 8. The lower panels show the ratios of the theoretical predictions to the measurements. Error bars around the experimental points show the statistical uncertainty, while the cross-hatched bands indicate the statistical and systematic uncertainties added in quadrature. The boxes around the MG5 aMC + PYTHIA 8 to measurement ratio represent the uncertainty on the prediction, including statistical, theoretical (from scale variations), and PDF uncertainties. The dark green area represents the statistical and theoretical uncertainties only, while the light green area represents the statistical uncertainty alone.

[GeV]
As the jet activity increases the Z boson recoils against a combination of several jets and this peak broadens, leading to an overall flattening of the distribution.
For the azimuthal angle between the Z boson and the second-and third-leading jets, as shown in Fig. 19, predictions and measurement agree very well. The differential cross sections are measured for the phase space regions with p T (Z) > 150 GeV and p T (Z) > 300 GeV. The results are shown in Figs. 20-23. The agreement of the predictions with the data is preserved, but the tree-level prediction computed with MADGRAPH 5 + PYTHIA 6 is an overestimate compared to the data at low azimuthal angle for the leading jet. The distributions are more uniform than in the ∆φ(Z, j 1 ) case, but retain a peak close to π. In the ∆φ(Z, j 2 ) case, we also see that the distributions show a larger correlation and a peak emerges at approximately ∆φ(Z, j 2 ) ≈ 2.6. This peak becomes more pronounced as the p T (Z) threshold increases. A similar trend is seen in the ∆φ(Z, j 3 ) distribution: selecting a high Z boson p T increases the fraction of events where the jets recoil against the boson.
Inclusive three-jet production is investigated in regions where both H T and the p T (Z) are large. Good agreement between data and predictions is also present here, as shown in Fig. 24. In this high-p T (Z), high-H T regime, we see a similar behaviour to the other high-p T (Z) selections. The ∆φ (Z, j 2 ) and ∆φ Z, j 3 distributions are also flatter than the corresponding distributions with no H T cut.  Figure 18: The differential cross section as a function of the azimuthal angle between the Z boson and the leading jet for different jet multiplicities, (left) N jets ≥ 1, (middle) N jets ≥ 2, and (right) N jets ≥ 3. The lower panels show the ratios of the theoretical predictions to the measurements. Error bars around the experimental points show the statistical uncertainty, while the cross-hatched bands indicate the statistical and systematic uncertainties added in quadrature. The boxes around the MG5 aMC + PYTHIA 8 to measurement ratio represent the uncertainty on the prediction, including statistical, theoretical (from scale variations), and PDF uncertainties. The dark green area represents the statistical and theoretical uncertainties only, while the light green area represents the statistical uncertainty alone. Figure 25 shows the azimuthal angle between the jets in the three-jet inclusive selections. The bumps seen at ∆φ ∼ 0.5 come from events with the two leading jets close in rapidity, |∆y| 2R, where R is the radius parameter of the jet anti-k t clustering algorithm, R = 0.5. This region is sensitive to the transition from an area of hadronic activity being resolved as one jet to being resolved as two jets. Increasing the p T (Z) threshold value to 150 GeV (Fig. 26) Figure 19: The differential cross section for Z (→ ) + jets production for N jets ≥ 3 as a function of the azimuthal angle between (left) the Z boson and the second leading jet, (right) the Z boson and the third leading jet. The lower panels show the ratios of the theoretical predictions to the measurements. Error bars around the experimental points show the statistical uncertainty, while the cross-hatched bands indicate the statistical and systematic uncertainties added in quadrature. The boxes around the MG5 aMC + PYTHIA 8 to measurement ratio represent the uncertainty on the prediction, including statistical, theoretical (from scale variations), and PDF uncertainties. The dark green area represents the statistical and theoretical uncertainties only, while the light green area represents the statistical uncertainty alone.
a dijet system radiates a Z boson are thus largely suppressed and this is most evident in the ∆φ (j 1 , j 2 ) distribution, where the peak at π is gone. A further increase in the p T (Z) threshold to 300 GeV (Fig. 27) continues this trend. In all cases, the agreement between the measurement and the prediction is still very good.    Figure 21: The differential cross section for Z (→ ) + jets production for N jets ≥ 3 and p T (Z) > 150 GeV as a function of the azimuthal angle (left) between the Z boson and the second leading jet and (right) between Z boson and third-leading jet. The lower panels show the ratios of the theoretical predictions to the measurements. Error bars around the experimental points show the statistical uncertainty, while the cross-hatched bands indicate the statistical and systematic uncertainties added in quadrature. The boxes around the MG5 aMC + PYTHIA 8 to measurement ratio represent the uncertainty on the prediction, including statistical, theoretical (from scale variations), and PDF uncertainties. The dark green area represents the statistical and theoretical uncertainties only, while the light green area represents the statistical uncertainty alone.  The differential cross section for Z (→ ) + jets production for N jets ≥ 3 and p T (Z) > 300 GeV as a function of the azimuthal angle (left) between the Z boson and the second-leading jet and (right) between the Z boson and the third-leading jet. The lower panels show the ratios of the theoretical predictions to the measurements. Error bars around the experimental points show the statistical uncertainty, while the cross-hatched bands indicate the statistical and systematic uncertainties added in quadrature. The boxes around the MG5 aMC + PYTHIA 8 to measurement ratio represent the uncertainty on the prediction, including statistical, theoretical (from scale variations), and PDF uncertainties. The dark green area represents the statistical and theoretical uncertainties only, while the light green area represents the statistical uncertainty alone.  Figure 24: The differential cross section for Z (→ ) + jets production for N jets ≥ 3, p Z T > 150 GeV, and H jet T > 300 GeV as a function of the azimuthal angle between the Z boson and the (left) first-, (middle) second-, and (right) third-leading jet. The lower panels show the ratios of the theoretical predictions to the measurements. Error bars around the experimental points show the statistical uncertainty, while the cross-hatched bands indicate the statistical and systematic uncertainties added in quadrature. The boxes around the MG5 aMC + PYTHIA 8 to measurement ratio represent the uncertainty on the prediction, including statistical, theoretical (from scale variations), and PDF uncertainties. The dark green area represents the statistical and theoretical uncertainties only, while the light green area represents the statistical uncertainty alone.
Overall, the measurements show that Monte Carlo predictions offer a very good description of the azimuthal angles between the jets and the Z boson, achieved when several parton multiplicities are included in the ME calculations and matched with parton showering.  Figure 25: The differential cross section for Z (→ ) + jets production for N jets ≥ 3 as a function of the azimuthal angle between (left) the first-and second-, (middle) the first-and third-, and (right) the second-and third-leading jets. The lower panels show the ratios of the theoretical predictions to the measurements. Error bars around the experimental points show the statistical uncertainty, while the cross-hatched bands indicate the statistical and systematic uncertainties added in quadrature. The boxes around the MG5 aMC + PYTHIA 8 to measurement ratio represent the uncertainty on the prediction, including statistical, theoretical (from scale variations), and PDF uncertainties. The dark green area represents the statistical and theoretical uncertainties only, while the light green area represents the statistical uncertainty alone.

Differential cross section for the dijet invariant mass
The dijet invariant mass is an important variable in the study of the production of a Higgs boson through vector boson fusion, since it can be used to select such events, which contain two jets well-separated in rapidity with a large dijet mass. For this measurement we consider all Z + jets events with at least two jets. The measured cross section as a function of the dijet mass is shown in Fig. 28.
The three predictions considered here agree well with the measurement within the experimental uncertainties, except for a dijet mass below ∼50 GeV, where the predictions made with MADGRAPH 5 + PYTHIA 6 and MG5 aMC + PYTHIA 8 show a deficit with respect to the measurements, while SHERPA 2 has a better agreement with the measurement in this region. In this region there is a relatively small angle between the two jets. The distribution of the difference in the rapidities, which are directly linked to the polar angle θ for massless objects (y = η = − ln[tan(θ/2)]) is well reproduced by the three predictions. The distribution of the angle in the transverse plane between the two jets is also well reproduced by all three calculations.

Multidimensional differential cross sections
The large number of Z+ ≥ 1 jet events allows the measurement of multidimensional cross sections. We focus on three observables, p T (j 1 ), y(Z), and y(j 1 ), that describe the kinematics of the events. Three differential cross sections are measured: d 2 σ/dp T (j 1 )dy(j 1 ), d 2 σ/dy(Z)dy(j 1 ),  Figure 26: The differential cross section for Z (→ ) + jets production for N jets ≥ 3 and p T (Z) > 150 GeV as a function of the azimuthal angle between (left) the first-and second-, (middle) the first-and third-, and (right) the second-and third-leading jets. The lower panels show the ratios of the theoretical predictions to the measurements. Error bars around the experimental points show the statistical uncertainty, while the cross-hatched bands indicate the statistical and systematic uncertainties added in quadrature. The boxes around the MG5 aMC + PYTHIA 8 to measurement ratio represent the uncertainty on the prediction, including statistical, theoretical (from scale variations), and PDF uncertainties. The dark green area represents the statistical and theoretical uncertainties only, while the light green area represents the statistical uncertainty alone. and d 3 σ/dp T (j 1 )dy(j 1 )dy(Z). The symmetry with respect to the transverse plane y = 0 is used to minimise the statistical uncertainties: d 2 σ/dp T (j 1 )dy(j 1 ) is obtained from a two-dimensional histogram of (p T , |y(j 1 )|) and d 2 σ/dy(Z)dy(j 1 ) from a histogram of (|y(j 1 )|, |y(Z)| sign(y(Z)y(j 1 )), where sign(x) = 1 for x ≥ 0 and sign(x) = −1, for x < 0. The three-dimensional differential cross section is calculated similarly.
The d 2 σ/dp T (j 1 )d|y(j 1 )| measurement, shown in Fig. 29, corresponds to the range p T < 550 GeV of the dσ/dp T measurement shown in Fig. 3 and extends the jet absolute rapidity range up to 4.7. The ratios of the theoretical predictions obtained from MADGRAPH 5 + PYTHIA 6, SHERPA 2, and MG5 aMC + PYTHIA 8 to the measurement are presented in Figs. 30-32. The difference in the shapes of the dσ/dp T spectrum between the measurement and the predictions computed with MADGRAPH 5 + PYTHIA 6 increases when moving from the central region, |y(j 1 )| = 0 to the more forward region, |y(j 1 )| = 2.5. The comparison of the SHERPA 2, and MG5 aMC + PYTHIA 8 predictions with the measurement does not show any dependence on the rapidity of the jet for the region |y(j 1 )| < 2.5, within the statistical uncertainty of the prediction, that is larger than for the MADGRAPH 5 + PYTHIA 6 sample. In the region beyond |y(j 1 )| = 2.5 the MADGRAPH 5 + PYTHIA 6 prediction-to-measurement ratio shows the same feature as for |y(j 1 )| < 2.5, despite the large experimental uncertainties due to a larger jet energy scale uncertainty. The SHERPA 2 prediction shows a significant difference with the spectrum of the jet transverse momentum being narrower than in data. The MG5 aMC + PYTHIA 8 shows a similar feature, but less pronounced and covered by the experimental uncertainties.
While the Z boson and jet rapidity distributions are independently well modelled by the simulation, we see in Section 6.4 that it is not the case with the tree-level calculations for the correlations between these two observables. Figures 33-35 Figure 27: The differential cross section for Z (→ ) + jets production for N jets ≥ 3 and p T (Z) > 300 GeV as a function of the azimuthal angle between (left) the first-and second-, (middle) the first-and third-, and (right) the second-and third-leading-jets. The lower panels show the ratios of the theoretical predictions to the measurements. Error bars around the experimental points show the statistical uncertainty, while the cross-hatched bands indicate the statistical and systematic uncertainties added in quadrature. The boxes around the MG5 aMC + PYTHIA 8 to measurement ratio represent the uncertainty on the prediction, including statistical, theoretical (from scale variations), and PDF uncertainties. The dark green area represents the statistical and theoretical uncertainties only, while the light green area represents the statistical uncertainty alone.
with respect to both rapidities. When the Z boson is central, the MADGRAPH 5 + PYTHIA 6 calculation predicts a more central leading jet, while when it is forward, it predicts a more forward leading jet in the same hemisphere (y(Z) y(j 1 ) > 0). These results are consistent with the measurement presented in Section 6.4 which showed that MADGRAPH 5 + PYTHIA 6 predicts a smaller y diff (Fig. 6). The predictions from SHERPA 2, and MG5 aMC + PYTHIA 8 agree well with the measurement when the jet is in the central region |y(j 1 )| < 2.5, while discrepancies start to appear when it is more forward. The tail of the jet rapidity is larger in the prediction, especially when the Z boson and the jets are well separated in rapidity: the discrepancy for y(Z) y(j 1 ) < 0 is larger for higher |y(Z)|. The discrepancies are more pronounced for the prediction obtained with SHERPA 2.
Finally, the measurement of the differential cross section with respect to both jet transverse momentum and rapidity is repeated for two different intervals of the Z boson rapidity as shown in Figs. 37-41. The shape of the ratio of the MADGRAPH 5 + PYTHIA 6 prediction to the measurement of the leading jet transverse momentum spectrum is similar in both intervals, although it shows a more pronounced discrepancy when the boson is in the most forward region. In the jet rapidity region |y(j 1 )| ∈ (1, 2.5) with y(j 1 ) y(Z) > 0, the ratios actually differ between the two Z boson rapidity intervals. However, in view of the measurement uncertainties this discrepancy cannot be considered significant. The behaviour seen previously for d 2 σ/(dy(Z)dy(j 1 )) translates into global shifts of the ratio distributions depending on the y(Z) y(j 1 ) interval. The bottom plots of Figs. 42 and 43 give more insight for the discrepancy with respect to the measurement observed previously for the SHERPA 2, and MG5 aMC + PYTHIA 8 predictions when the jet is in the forward region, |y(j 1 )| ∈ (2.5, 4.7). The observed deficit in the cross section can be attributed to soft jets, since more events with the leading jet below 90 GeV are expected from the prediction. The discrepancy is larger when the Z boson and the leading jet are well The differential cross section for Z (→ ) + jets production as a function of the leading jet transverse momentum and rapidity. The bands around the measurement points represent the total measurement uncertainties. The bands around the prediction points represent the total uncertainty, and its statistical, theoretical, and PDF components for MG5 aMC + PYTHIA 8, and the statistical uncertainty alone for MADGRAPH 5 + PYTHIA 6 and SHERPA 2. separated in rapidity. Indeed, the discrepancy is the smallest for the region |y(Z)| ∈ (1, 2.5) and y(Z) y(j 1 ) > 0, corresponding to the region where the rapidities of the boson and the jet are the closest in the jet rapidity range considered.

Summary
The kinematics of Z + jets events in pp collisions at the centre-of-mass energy of √ s = 8 TeV have been studied and the differential cross sections have been measured as a function of numerous observables. Multidimensional cross section measurements have been performed with respect to up to three variables. The results have been compared with predictions from several multileg generators at different fixed-order accuracies, tree-level and NLO up to 2 partons, and employing different showering algorithms, as implemented in PYTHIA 6, PYTHIA 8, and SHERPA 2.
The comparisons show that it is essential to include a large number of final-state partons in the matrix element calculations in order to correctly describe the kinematics of the leading jets. Besides the individual jet p T , the observable H T , used in searches for physics beyond the SM and defined in this measurement for jets with p T > 30 GeV, is modelled correctly at low values of H T only when a sufficiently large number of partons is included in the matrix element calculations. The discrepancies found for large values of the jet momentum, first observed in  Figure 30: Ratio to the measurement of the differential cross section d 2 σ/dp T (j 1 )dy(j 1 ) obtained with MADGRAPH 5 + PYTHIA 6, with up to four jets at LO. The total experimental uncertainty is shown as a band around 1. Uncertainties in the predictions are shown on the ratio points and include the statistical uncertainty only.  Figure 31: Ratio to the measurement of the differential cross section d 2 σ/dp T (j 1 )dy(j 1 ) obtained with SHERPA 2, with up to two jets at NLO and up to four jets at LO. The total experimental uncertainty is shown as a band around 1. Uncertainties in the predictions are shown on the ratio points and include the statistical uncertainty only.  Figure 32: Ratio to the measurement of the differential cross section d 2 σ/dp T (j 1 )dy(j 1 ) obtained with MG5 aMC + PYTHIA 8, with up to two jets at NLO. The total experimental uncertainty is shown as a band around 1. Uncertainties in the predictions are shown on the ratio points and include the statistical, theoretical, and PDF uncertainties. The dark green area represents the statistical and theoretical uncertainties only, while the light green area represents the statistical uncertainty alone. The differential cross section for Z (→ ) + jets production as a function of the Z boson and leading jet rapidity. The bands around the measurement points represent the total measurement uncertainties. The bands around the prediction points represent the total uncertainty, statistical, theoretical, and PDF components for MG5 aMC + PYTHIA 8, and the statistical uncertainty alone for MADGRAPH 5 + PYTHIA 6 and SHERPA 2.  Figure 36: Ratio to the measurement of the differential cross section d 2 σ/dy(Z)dy(j 1 ) obtained with MG5 aMC + PYTHIA 8, with up to two jets at NLO. The total experimental uncertainty is shown as a band around 1. Uncertainties in the predictions are shown on the ratio points and include the statistical, theoretical, and PDF uncertainties. The dark green area represents the statistical and theoretical uncertainties only, while the light green area represents the statistical uncertainty alone. the √ s = 7 TeV measurements [5,8], are confirmed at √ s = 8 TeV with a larger data set. Such discrepancies are not seen when including the NLO corrections. The differences observed between tree-level predictions and the measurements of the leading jet are larger when the jet is more forward (|y| > 2.5). Discrepancies with LO and NLO predictions have been observed for the dijet mass spectrum at low mass in the region where the angle between the two jet directions is smaller than π/2. Nevertheless, the azimuthal angles between the Z boson and the jet and between the jets are very well reproduced by the predictions including the tree-level one. The excellent agreement remains when restricting the phase space by applying a threshold on the Z boson p T , on H T , or on both. The rapidity distributions of the Z boson and jets are fairly well modelled by the generators, but the correlations between the rapidities, which have been studied by measuring multidimensional differential cross sections and distributions of rapidity differences and sums, are not well reproduced by the multileg tree-level calculation. We have shown that the multileg event generators including NLO terms reproduce the rapidity difference distributions very well. The rapidity sum is also successfully described. For this variable the discrepancy with the tree-level calculation could also be due to a different choice of the parton distribution functions.
In summary, kinematics of Z + jets events have been studied in detail and apart from a few discrepancies, the measurements show a very good agreement with the considered NLO multileg predictions.  37: The differential cross section for Z (→ ) + jets production as a function of the rapidities of the Z boson and leading jet, and of the transverse momentum of the jet for the configuration. The bottom plots correspond to the configuration where the boson and the jet are in different hemispheres (y(Z)y(j 1 ) < 0), while the top plots correspond to both objects in the same hemisphere. The left and right plots show the respective Z boson rapidity ranges, |y(Z)| < 1 and |y(Z)| ∈ (1, 2.5). The bands around the measurement points represent the total measurement uncertainties. The bands around the prediction points represent the total uncertainty, statistical, theoretical, and PDF components for MG5 aMC + PYTHIA 8, and the statistical uncertainty alone for MADGRAPH 5 + PYTHIA 6 and SHERPA 2.  Figure 38: Ratio to the measurement of the differential cross section d 3 σ/dp T (j 1 )dy(j 1 )dy(Z) obtained with MADGRAPH 5 + PYTHIA 6, with up to four jets at LO, for |y(Z)| < 1. Left column corresponds to y(j 1 )y(Z) > 0 and right column to y(j 1 )y(Z) < 0. The total experimental uncertainty is shown as a band around 1. Uncertainties in the predictions are shown on the ratio points and include the statistical uncertainty only.  Figure 39: Ratio to the measurement of the differential cross section d 3 σ/dp T (j 1 )dy(j 1 )dy(Z) obtained with MADGRAPH 5 + PYTHIA 6, with up to four jets at LO, for |y(Z)| ∈ (1, 2.5). Left column corresponds to y(j 1 )y(Z) > 0 and right column to y(j 1 )y(Z) < 0. The total experimental uncertainty is shown as a band around 1. Uncertainties in the predictions are shown on the ratio points and include the statistical uncertainty only.  Figure 40: Ratio to the measurement of the differential cross section d 3 σ/dp T (j 1 )dy(j 1 )dy(Z) obtained with SHERPA 2, with up to two jets at NLO and up to four jets at LO, for |y(Z)| < 1. Left column corresponds to y(j 1 )y(Z) > 0 and right column to y(j 1 )y(Z) < 0. The total experimental uncertainty is shown as a band around 1. Uncertainties in the predictions are shown on the ratio points and include the statistical uncertainty only.  Figure 41: Ratio to the measurement of the differential cross section d 3 σ/dp T (j 1 )dy(j 1 )dy(Z) obtained with SHERPA 2, with up to two jets at NLO and up to four jets at LO, for |y(Z)| ∈ (1, 2.5). Left column corresponds to y(j 1 )y(Z) > 0 and right column to y(j 1 )y(Z) < 0. The total experimental uncertainty is shown as a band around 1. Uncertainties in the predictions are shown on the ratio points and include the statistical uncertainty only.  Figure 42: Ratio to the measurement of the differential cross section d 3 σ/dp T (j 1 )dy(j 1 )dy(Z) obtained with MG5 aMC + PYTHIA 8, with up to two jets at NLO, for |y(Z)| < 1. Left column corresponds to y(j 1 )y(Z) > 0 and right column to y(j 1 )y(Z) < 0. The total experimental uncertainty is shown as a band around 1. Uncertainties in the predictions are shown on the ratio points and include the statistical, theoretical, and PDF uncertainties. The dark green area represents the statistical and theoretical uncertainties only, while the light green area represents the statistical uncertainty alone.  Figure 43: Ratio to the measurement of the differential cross section d 3 σ/dp T (j 1 )dy(j 1 )dy(Z) obtained with MG5 aMC + PYTHIA 8, with up to two jets at NLO, for |y(Z)| ∈ (1, 2.5). Left column corresponds to y(j 1 )y(Z) > 0 and right column to y(j 1 )y(Z) < 0. The total experimental uncertainty is shown as a band around 1. Uncertainties in the predictions are shown on the ratio points and include the statistical, theoretical, and PDF uncertainties. The dark green area represents the statistical and theoretical uncertainties only, while the light green area represents the statistical uncertainty alone.

Acknowledgments
We congratulate our colleagues in the CERN accelerator departments for the excellent performance of the LHC and thank the technical and administrative staffs at CERN and at other CMS institutes for their contributions to the success of the CMS effort. In addition, we gratefully acknowledge the computing centres and personnel of the Worldwide LHC Computing Grid for delivering so effectively the computing infrastructure essential to our analyses. Finally, we acknowledge the enduring support for the construction and operation of the LHC and the CMS detector provided by the following funding agencies: the Austrian Federal Ministry of Science, Research and Economy and the Austrian Science Fund; the Belgian