1 Introduction

The measurement of the production cross-section of jets, a collimated spray of hadrons, in association with a \(Z\) boson (\(Z \text { + jets}\)), is an important process for testing the predictions of perturbative quantum chromodynamics (pQCD). It provides a benchmark for fixed-order calculations and predictions from Monte Carlo (MC) simulations, which are often used to estimate the \(Z \text { + jets}\) background in the measurements of Standard Model processes, such as Higgs boson production, and in searches for new physics beyond the Standard Model.

Various properties of \(Z \text { + jets}\) production have been measured in proton–antiproton collisions at \(\sqrt{s} =1.96\) TeV at the Tevatron [1,2,3,4]. The differential \(Z \text { + jets}\) cross-section is measured as functions of the \(Z\) boson transverse momentum and the jets’ transverse momenta and rapidities, and as a function of the angular separation between the \(Z\) boson and jets in final states with different jet multiplicities. The experiments at the Large Hadron Collider (LHC) [5] have an increased phase space compared to previous measurements by using proton–proton collision data at \(\sqrt{s} =7\), 8 and 13 TeV [6,7,8,9,10,11,12,13,14,15]. The measurements at the LHC allow state-of-the-art theoretical \(Z \text { + jets}\) predictions to be tested. These have recently been calculated to next-to-next-to-leading-order (NNLO) accuracy in pQCD [16, 17].

This paper studies the double-differential cross-section of inclusive jet production in association with a \(Z\) boson which decays into an electron–positron pair. The cross-section is measured as a function of absolute jet rapidity, \(|y_\text {jet}|\), and jet transverse momentum, \(p_{\text {T}} ^\text {jet}\), using the proton–proton (pp) collision data at \(\sqrt{s}=8\, \text {Te}\text {V}\) collected by the ATLAS experiment. The measured cross-section is unfolded to the particle level.

The cross-section calculated at fixed order for \(Z \text { + jets}\) production in \(pp\) collisions at \(\sqrt{s} =8\) TeV is dominated by quark–gluon scattering. The \(Z \text { + jets}\) cross-section is sensitive to the gluon and sea-quark parton distribution functions (PDFs) in the proton. In the central \(|y_\text {jet}|\) region the \(Z \text { + jets}\) cross-section probes the PDFs in the low x range, where x is the proton momentum fraction, while in the forward \(|y_\text {jet}|\) region it examines the intermediate and high x values. The scale of the probe is set by \(p_{\text {T}} ^\text {jet}\).

The measured cross-section is compared with the next-to-leading-order (NLO) and NNLO \(Z \text { + jets}\) fixed-order calculations, corrected for hadronisation and the underlying event. In addition, the data are compared with the predictions from multi-leg matrix element (ME) MC event generators supplemented with parton shower simulations.

The structure of the paper is as follows. The ATLAS detector is briefly described in Sect. 2. This is followed by a description of the data in Sect. 3 and the simulated samples in Sect. 4. The definition of the object reconstruction, calibration and identification procedures and a summary of the selection criteria are given in Sect. 5. The \(Z \text { + jets}\) backgrounds are discussed in Sect. 6. The correction of the measured spectrum to the particle level is described in Sect. 7. The experimental uncertainties are discussed in Sect. 8. The fixed-order calculations together with parton-to-particle-level corrections are presented in Sect. 9. Finally, the measured cross-section is presented and compared with the theory predictions in Sect. 10. The quantitative comparisons with the fixed-order pQCD predictions are summarised in Sect. 11.

2 The ATLAS detector

The ATLAS experiment [18] at the LHC is a multipurpose particle detector with a forward–backward symmetric cylindrical geometry and nearly \(4\pi \) coverage in solid angle.Footnote 1 It consists of an inner tracking detector surrounded by a thin superconducting solenoid, electromagnetic and hadronic calorimeters, and a muon spectrometer incorporating three large superconducting toroidal magnets.

The inner-detector system (ID) is immersed in a 2 T axial magnetic field and provides charged-particle tracking in the range \(|\eta | < 2.5\). A high-granularity silicon pixel detector covers the \(pp\) interaction region and typically provides three measurements per track. It is followed by a silicon microstrip tracker (SCT), which usually provides four two-dimensional measurement points per track. These silicon detectors are complemented by a transition radiation tracker (TRT), which provides electron identification information.

The calorimeter system covers the pseudorapidity range \(|\eta | < 4.9\). In the region \(|\eta |< 3.2\), electromagnetic calorimetry is provided by barrel and endcap high-granularity lead/liquid-argon (LAr) electromagnetic calorimeters, with an additional thin LAr presampler covering \(|\eta | < 1.8\) to correct for energy loss in material upstream of the calorimeters. Hadronic calorimetry is provided by a steel/scintillator-tile calorimeter, segmented into three barrel structures for \(|\eta | < 1.7\), and two copper/LAr hadronic endcap calorimeters in the range \(1.5<|\eta |< 3.2\). The calorimetry in the forward pseudorapidity region, \(3.1<|\eta |< 4.9\), is provided by the copper-tungsten/LAr calorimeters.

The muon spectrometer surrounds the calorimeters and contains three large air-core toroidal superconducting magnets with eight coils each. The field integral of the toroids ranges between 2.0 and 6.0 T m across most of the detector. The muon spectrometer includes a system of precision tracking chambers and fast detectors for triggering.

A three-level trigger [19] was used to select events for offline analysis. The first-level trigger is implemented in hardware and used a subset of the detector information to reduce the accepted rate to at most 75 kHz. This was followed by two software-based trigger levels that together reduced the average accepted event rate to 400 Hz.

3 Data sample

The data used for this analysis are from proton–proton collisions at \(\sqrt{s}=8\) TeV that were collected by the ATLAS detector in 2012 during stable beam conditions. Events recorded when any of the ATLAS subsystems were defective or non-operational are excluded. Data were selected with a dielectron trigger, which required two reconstructed electron candidates with transverse momenta greater than 12 GeV. Only events with electron energy leakage of less than 1 GeV into the hadronic calorimeter were accepted. The trigger required that reconstructed electron candidates were identified using the ‘loose1’ criteria [20], which are discussed in Sect. 5.1.

The integrated luminosity of the analysis data sample after the trigger selection is \(19.9~\text{ fb }^{-1}\) measured with an uncertainty of \(\pm 1.9\)% [21]. The average number of simultaneous proton–proton interactions per bunch crossing is 20.7.

In addition, a special data sample was selected for a data-driven study of multijet and \(W \text {+ jets}\) backgrounds. For this purpose, the analysis data sample was enlarged by including auxiliary events selected by a logical OR of two single-electron triggers.

The first single-electron trigger required events with at least one reconstructed electron candidate with a transverse momentum greater than 24 GeV and hadronic energy leakage less than 1 GeV. The electron candidate satisfied the ‘medium1’ identification criteria [20], a tightened subset of ‘loose1’. The reconstructed electron track was required to be isolated from other tracks in the event. The isolation requirement rejected an event if the scalar sum of reconstructed track transverse momenta in a cone of size \(\Delta R=0.2\) around the electron track exceeded 10% of the electron track’s transverse momentum.

The second single-electron trigger accepted events with at least one electron candidate with a transverse momentum greater than 60 GeV and identified as ‘medium1’. This trigger reduced inefficiencies in events with high-\(p_{\text {T}}\) electrons that resulted from the isolation requirement used in the first trigger.

Events selected by single-electron triggers include a large number of background events that are normally rejected by the \(Z \text { + jets}\) selection requirements, but these events are used in the data-driven background studies.

4 Monte Carlo simulations

Simulated \(Z \text { + jets}\) signal events were generated using the Sherpa v. 1.4 [22] multi-leg matrix element MC generator. The MEs were calculated at NLO accuracy for the inclusive \(Z\) production process, and additionally with LO accuracy for up to five partons in the final state, using Amegic++ [23]. Sherpa MEs were convolved with the CT10 [24] PDFs. Sherpa parton showers were matched to MEs following the CKKW scheme [25]. The MENLOPS [26] prescription was used to combine different parton multiplicities from matrix elements and parton showers. Sherpa predictions were normalised to the inclusive \(Z\) boson production cross-section calculated at NNLO [27,28,29] and are used for the unfolding to particle level and for the evaluation of systematic uncertainties.

An additional \(Z \text { + jets}\) signal sample with up to five partons in the final state at LO was generated using Alpgen v. 2.14 [30]. The parton showers were generated using Pythia v. 6.426 [31] with the Perugia 2011C [32] set of tuned parameters to model the underlying event’s contribution. The Alpgen MEs were matched to the parton showers following the MLM prescription [30]. The proton structure was described by the CTEQ6L1 [33] PDF. Referred to as the Alpgen+Pythia sample, these predictions were normalised to the NNLO cross-section. This sample is used in the analysis for the unfolding uncertainty evaluation and for comparisons with the measurement.

The five-flavour scheme with massless quarks was used in both the Sherpa and Alpgen+Pythia predictions.

Backgrounds from the \(Z \rightarrow \tau \tau \), diboson (WW, WZ and ZZ), \(t\bar{t}\) and single-top-quark events are estimated using MC simulations. The \(Z \rightarrow \tau \tau \) events were generated using Powheg-Box v. 1.0 [34, 35] interfaced to Pythia v. 8.160 [36] for parton showering using the CT10 PDFs and the AU2 [37] set of tuned parameters. The \(Z \rightarrow \tau \tau \) prediction was normalised to the NNLO cross-section [27,28,29]. The WW, WZ and ZZ events were generated using Herwig v. 6.520.2 [38] with the CTEQ6L1 PDFs and the AUET2 [39] set of tuned parameters. The diboson predictions were normalised to the NLO cross-sections [40, 41]. Samples of single-top-quark events, produced via the s-, t- and \(W t\)-channels, and \(t\bar{t}\) events were generated with Powheg-Box v. 1.0 interfaced to Pythia v. 6.426, which used the CTEQ6L1 PDFs and the Perugia 2011C set of tuned parameters. The prediction for single-top-quark production in s-channel were normalised to the NNLO calculations matched to the next-to-next-to-leading-logarithm (NNLL) calculations (NNLO+NNLL) [42], while predictions in t- and \(W t\)-channel are normalised to the NLO+NNLL calculations [43, 44]. The \(t\bar{t}\) samples were normalised to the NNLO+NNLL calculations [45].

The Photos [46] and Tauola [47] programs were interfaced to the MC generators, excluding Sherpa, to model electromagnetic final-state radiation and \(\tau \)-lepton decays, respectively.

Additional proton–proton interactions, generally called pile-up, were simulated using the Pythia v. 8.160 generator with the MSTW2008 [48] PDFs and the A2 [37] set of tuned parameters. The pile-up events were overlaid onto the events from the hard-scattering physics processes. MC simulated events were reweighted to match the distribution of the average number of interactions per bunch crossing in data.

All MC predictions were obtained from events processed with the ATLAS detector simulation [49] that is based on Geant 4 [50].

5 Object definitions and event selection

The measured objects are the electrons and jets reconstructed in ATLAS. The methods used to reconstruct, identify and calibrate electrons are presented in Sect. 5.1. The reconstruction of jets, their calibration, and background suppression methods are discussed in Sect. 5.2. Finally, all selection requirements are summarised in Sect. 5.3.

5.1 Electron reconstruction and identification

Electron reconstruction in the central region, \(|\eta |<2.5\), starts from energy deposits in calorimeter cells. A sliding-window algorithm scans the central volume of the electromagnetic calorimeter in order to seed three-dimensional clusters. The window has a size of \(3\times 5\) in units of \(0.025\times 0.025\) in \(\eta \)\(\phi \) space. Seeded cells have an energy sum of the constituent calorimeter cells greater than \(2.5\, \text {Ge}\text{V}\). An electron candidate is reconstructed if the cluster is matched to at least one track assigned to the primary vertex, as measured in the inner detector. The energy of a reconstructed electron candidate is given by the energy of a cluster that is enlarged to a size of \(3\times 7\) (\(5\times 5\)) in \(\eta \)\(\phi \) space in the central (endcap) electromagnetic calorimeter in order to take into account the shape of electromagnetic shower energy deposits in different calorimeter regions. The \(\eta \) and \(\phi \) coordinates of a reconstructed electron candidate are taken from the matched track. The details of the electron reconstruction are given in Ref. [51].

A multistep calibration is used to correct the electron energy scale to that of simulated electrons [52]. Cluster energies in data and in MC simulation are corrected for energy loss in the material upstream of the electromagnetic calorimeter, energy lost outside of the cluster volume and energy leakage beyond the electromagnetic calorimeter. The reconstructed electron energy in data is corrected as a function of electron pseudorapidity using a multiplicative scale factor obtained from a comparison of \(Z \rightarrow ee\) mass distributions between data and simulation. In addition, the electron energy in the MC simulation is scaled by a random number taken from a Gaussian distribution with a mean value of one and an \(\eta \)-dependent width, equal to the difference between the electron energy resolution in data and MC simulation, determined in situ using \(Z \rightarrow ee\) events.

A set of cut-based electron identification criteria, which use cluster shape and track properties, is applied to reconstructed electrons to suppress the residual backgrounds from photon conversions, jets misidentified as electrons and semileptonic heavy-hadron decays. There are three types of identification criteria, listed in the order of increasing background-rejection strength but diminishing electron selection efficiency: ‘loose’, ‘medium’ and ‘tight’ [51]. The ‘loose’ criteria identify electrons using a set of thresholds applied to cluster shape properties measured in the first and second LAr calorimeter layers, energy leakage into the hadronic calorimeter, the number of hits in the pixel and SCT detectors, and the angular distance between the cluster position in the first LAr calorimeter layer and the extrapolated track. The ‘medium’ selection tightens ‘loose’ requirements on shower shape variables. In addition, the ‘medium’ selection sets conditions on the energy deposited in the third calorimeter layer, track properties in the TRT detector and the vertex position. The ‘tight’ selection tightens the ‘medium’ identification criteria thresholds, sets conditions on the measured ratio of cluster energy to track momentum and rejects reconstructed electron candidates matched to photon conversions.

Each MC simulated event is reweighted by scale factors that make the trigger, reconstruction and identification efficiencies the same in data and MC simulation. The scale factors are generally close to one and are calculated in bins of electron transverse momenta and pseudorapidity [20, 51].

5.2 Jet reconstruction, pile-up suppression and quality criteria

Jets are reconstructed using the anti-\(k_{t}\) algorithm [53] with a radius parameter \(R=0.4\), as implemented in the FastJet software package [54]. Jet reconstruction uses topologically clustered cells from both the electromagnetic and hadronic calorimeters [55]. The topological clustering algorithm groups cells with statistically significant energy deposits as a method to suppress noise. The energy scale of calorimeter cells is initially established for electromagnetic particles. The local cell weighting (LCW) [56] calibration is applied to topological clusters to correct for the difference between the detector responses to electromagnetic and hadronic particles, energy losses in inactive material and out-of-cluster energy deposits. The LCW corrections are derived using the MC simulation of the detector response to single pions.

The jet energy scale (JES) calibration [57] corrects the energy scale of reconstructed jets to that of simulated particle-level jets. The JES calibration includes origin correction, pile-up correction, MC-based correction of the jet energy and pseudorapidity (MCJES), global sequential calibration (GSC) and residual in situ calibration.

The origin correction forces the four-momentum of the jet to point to the hard-scatter primary vertex rather than to the centre of the detector, while keeping the jet energy constant.

Pile-up contributions to the measured jet energies are accounted for by using a two-step procedure. First, the reconstructed jet energy is corrected for the effect of pile-up by using the average energy density in the event and the area of the jet [58]. Second, a residual correction is applied to remove the remaining dependence of the jet energy on the number of reconstructed primary vertices, \(N_\text {PV}\), and the expected average number of interactions per bunch crossing, \(\langle \mu \rangle \).

The MCJES corrects the reconstructed jet energy to the particle-level jet energy using MC simulation. In addition, a correction is applied to the reconstructed jet pseudorapidity to account for the biases caused by the transition between different calorimeter regions and the differences in calorimeter granularity.

Next, the GSC corrects the jet four-momenta to reduce the response’s dependence on the flavour of the parton that initiates the jet. The GSC is determined using the number of tracks assigned to a jet, the \(p_{\text {T}}\) -weighted transverse distance in the \(\eta \)\(\phi \) space between the jet axis and all tracks assigned to the jet (track width), and the number of muon track segments assigned to the jet.

Finally, the residual in situ correction makes the jet response the same in data and MC simulation as a function of detector pseudorapidity by using dijet events (\(\eta \)-intercalibration), and as a function of jet transverse momentum by using well-calibrated reference objects in \(Z/\gamma \) and multijet events.

Jets originating from pile-up interactions are suppressed using the jet vertex fraction (JVF) [58]. The JVF is calculated for each jet and each primary vertex in the event as a ratio of the scalar sum of \(p_{\text {T}}\) of tracks, matched to a jet and assigned to a given vertex, to the scalar sum of \(p_{\text {T}}\) of all tracks matched to a jet.

Applying jet quality criteria suppresses jets from non-collision backgrounds that arise from proton interactions with the residual gas in the beam pipe, beam interactions with the collimator upstream of the ATLAS detector, cosmic rays overlapping in time with the proton–proton collision and noise in the calorimeter. Jet quality criteria are used to distinguish jets by using the information about the quality of the energy reconstruction in calorimeter cells, the direction of the shower development and the properties of tracks matched to jets. There are four sets of selection criteria that establish jet quality: ‘looser’, ‘loose’, ‘medium’ and ‘tight’. They are listed in the order of increasing suppression of non-collision jet background but decreasing jet selection efficiency. The ‘medium’ jet selection is used in the paper due to its high background rejection rate together with about 100% jet selection efficiency in the \(p_{\text {T}} ^\text {jet}>60\, \text {Ge}\text {V}\) region [57].

5.3 Event selection

Events are required to have a primary vertex with at least three assigned tracks that have a transverse momentum greater than 400 MeV. When several reconstructed primary vertices satisfy this requirement, the hard-scatter vertex is taken to be the one with the highest sum of the squares of the transverse momenta of its assigned tracks.

Each event is required to have exactly two reconstructed electrons, each with transverse momentum greater than \(20 \,\text {Ge}\text {V}\) and an absolute pseudorapidity less than 2.47, excluding the detector transition region, \(1.37< |\eta _{e}|< 1.52\), between barrel and endcap electromagnetic calorimeters. The electrons are required to have opposite charges, be identified using the ‘medium’ [51] criteria and be matched to electron candidates that were selected by the trigger. The ‘medium’ identification ensures the electrons originate from the hard-scatter vertex. The electron-pair invariant mass, \(m_{ee}\), is required to be in the \(66\,\text {Ge}\text {V}< m_{ee}< 116\,\text {Ge}\text {V}\) range.

Jets are required to have a transverse momentum greater than \(25\,\text {Ge}\text {V}\) and an absolute jet rapidity less than 3.4. Jets with \(p_{\text {T}} ^\text {jet}<50\,\text {Ge}\text {V}\), \(|\eta _\text {det}|<2.4\), where \(\eta _\text {det}\) is reconstructed relative to the detector centre, and \(|\text {JVF}|<0.25\) are considered to be from pile-up. Jets originating from pile-up are removed from the measurement. MC simulations poorly describe the effects of high pile-up in the \(p_{\text {T}} ^\text {jet}< 50\,\text {Ge}\text {V}\) and \(|y_\text {jet}|> 2.4\) region, so this region is not included in the measurement. Jets reconstructed within \(\Delta R = 0.4\) of selected electrons are rejected in order to avoid overlap. Jets are required to satisfy the ‘medium’ [57] quality criteria. In addition, jets in regions of the detector that are poorly modelled are rejected in data and MC simulations in order to avoid biasing the measured jet energy [57]. Each jet that meets the selection requirements is used in the measurement.

As a result, 1,486,415 events with two electrons and at least one jet were selected for the analysis.

6 Backgrounds

The majority of irreducible backgrounds in this measurement are studied using MC samples that simulate \(Z {} \rightarrow \tau \tau \), diboson, \(t\bar{t}\) and single top-quark production. The \(Z {} \rightarrow \tau \tau \) process is a background if both \(\tau \)-leptons decay into an electron and neutrino. Diboson production constitutes a background to the \(Z \text { + jets}\) signal if the \(W\) and/or \(Z\) boson decays into electrons. Since the top-quark decays predominantly via \(t \rightarrow Wb\), the \(t\bar{t}\) and single top-quark constitute a background to the \(Z \text { + jets}\) signal when \(W\) bosons decay into an electron or jets are misidentified as electrons.

Multijet production constitutes a background to the \(Z \text { + jets}\) signal when two jets are misidentified as electrons. The \(W \text {+ jets}\) background is due to an electron from \(W\) boson decay and a jet misidentified as electron. A combined background from multijet and \(W \text {+ jets}\) events is studied using a data-driven technique, thus providing a model-independent background estimate.

A background-enriched data sample is used for the combined multijet and \(W \text {+ jets}\) background control region. Its selection requires two reconstructed electrons with at least one electron that satisfies the ‘medium’ identification criteria, but not the ‘tight’ ones. This allows selection of events with at least one jet misidentified as an electron. No identification criteria are applied to the second reconstructed electron, in order to allow for the possibility of \(W \text {+ jets}\) events with a genuine electron from \(W\) boson decay, and multijet events with a jet misidentified as another electron. Both selected electrons are required to have the same charge to suppress the \(Z \text { + jets}\) signal events. The combined multijet and \(W \text {+ jets}\) background template is constructed by subtracting the MC simulated \(Z \text { + jets}\) signal events and \(Z {} \rightarrow \tau \tau \), diboson, \(t\bar{t}\) and single-top-quark background events in the control region from data.

The purity of the template is calculated as the fraction of multijet and \(W \text {+ jets}\) events in the data control region. The purity is about 98% in the tails of the \(m_{ee}\) distribution and is about 80% near the \(m_{ee}\) peak at \(91\,\text {Ge}\text {V}\). The template purity is above \(90\%\) in all \(|y_\text {jet}|\) and \(p_{\text {T}} ^\text {jet}\) bins.

The combined multijet and \(W \text {+ jets}\) background template is normalised to data using the invariant mass distribution of reconstructed electron pairs. A maximum-likelihood fit is used to adjust the normalisation of the combined multijet and \(W \text {+ jets}\) background template relative to the measured \(Z \text { + jets}\) distribution. The normalisations of MC simulated samples are fixed in the fit: the \(Z {} \rightarrow \tau \tau \), diboson, \(t\bar{t}\) and single-top-quark distributions are normalised by their fixed-order cross-sections, whereas the normalisation of the MC simulated \(Z \text { + jets}\) signal events is scaled to data to give the same total number of events near the peak of the \(Z\) mass spectrum in the \(90\,\text {Ge}\text {V}< m_{ee}< 92\,\text {Ge}\text {V}\) range. The combined multijet and \(W \text {+ jets}\) background is fit to data in an extended \(m_{ee}\) region, \(60\,\text {Ge}\text {V}< m_{ee}< 140\,\text {Ge}\text {V}\), excluding the bins under the \(Z\) peak within the \(80\,\text {Ge}\text {V}< m_{ee}< 100\,\text {Ge}\text {V}\) region. The extended \(m_{ee}\) region is used for the normalisation extraction only, as it allows more background events in the tails of the \(Z\) mass spectrum. The normalisation of the multijet and \(W \text {+ jets}\) background template, calculated in the fit, is used to adjust the templates, obtained in the \(|y_\text {jet}|\) and \(p_{\text {T}} ^\text {jet}\) bins, to the \(Z \text { + jets}\) signal region.

The total number of jets in \(Z \text { + jets}\) events are shown as a function of \(|y_\text {jet}|\) and \(p_{\text {T}} ^\text {jet}\) bins in Fig. 1. Data are compared with the sum of signal MC events and all backgrounds.

The Sherpa \(Z \text { + jets}\) simulation, normalised to the NNLO cross-section, is lower than data by about 10% in the \(p_{\text {T}} ^\text {jet}< 200\,\text {Ge}\text {V}\) region. These differences are mostly covered by the variations within electron and jet uncertainties introduced in Sect. 8. In the \(p_{\text {T}} ^\text {jet}>200\,\text {Ge}\text {V}\) region, agreement with data is within the statistical uncertainties.

The Alpgen+Pythia predictions are in agreement with data within 10% for jets with transverse momenta below 100 GeV. However, the level of disagreement increases as a function of the jet transverse momenta, reaching 30% in the \(400\,\text {Ge}\text {V}< p_{\text {T}} ^\text {jet}{} < 1050\,\text {Ge}\text {V}\) region.

The dominant background in the measurement is from \(t\bar{t}\) events. It is 0.3–0.8% in the \(25\,\text {Ge}\text {V}< p_{\text {T}} ^\text {jet}< 50\,\text {Ge}\text {V}\) region and 1–2.5% in the \(50\,\text {Ge}\text {V}< p_{\text {T}} ^\text {jet}< 100\,\text {Ge}\text {V}\) region, with the largest contribution in the central rapidity region. In the \(100\,\text {Ge}\text {V}< p_{\text {T}} ^\text {jet}< 200\,\text {Ge}\text {V}\) region, this background is approximately 3%, while in the \(200\,\text {Ge}\text {V}< p_{\text {T}} ^\text {jet}< 1050\, \text {Ge}\text {V}\) region it is 1.8–8%, increasing for forward rapidity jets.

The combined multijet and \(W \text {+ jets}\) background and the diboson background are similar in size. The contributions of these backgrounds are 0.5–1%.

The \(Z \rightarrow \tau \tau \) and single-top-quark backgrounds are below 0.1%.

Fig. 1
figure 1

The total number of jets in \(Z \text { + jets}\) events as a function of \(|y_\text {jet}|\) in \(p_{\text {T}} ^\text {jet}\) bins for the integrated luminosity of \(19.9~\text{ fb }^{-1}\). Data are presented with markers. The filled areas correspond to the backgrounds stacked. All backgrounds are added to the \(Z \text { + jets}\)Sherpa and Alpgen+Pythia predictions. Lower panels show ratios of MC predictions to data. The grey band shows the sum in quadrature of the electron and jet uncertainties. The statistical uncertainties are shown with vertical error bars. In the lower panels the total data + MC statistical uncertainty is shown

7 Unfolding of detector effects

The experimental measurements are affected by the detector resolution and reconstruction inefficiencies. In order to compare the measured cross-sections with the theoretical \(Z \text { + jets}\) predictions at the particle level, the reconstructed spectrum is corrected for detector effects using the iterative Bayesian unfolding method [59]. The unfolding is performed using the Sherpa \(Z \text { + jets}\) simulation.

The particle-level phase space in the MC simulation is defined using two dressed electrons and at least one jet. For the dressed electron, the four-momenta of any photons within a cone of \(\Delta R=0.1\) around its axis are added to the four-momentum of the electron. Electrons are required to have \(|\eta |<2.47\) and \(p_{\text {T}} >20\,\text {Ge}\text {V}\). The electron pair’s invariant mass is required to be within the range \(66\,\text {Ge}\text {V}<m_{ee}<116\,\text {Ge}\text {V}\).

Jets at the particle level are built by using the anti-\(k_{t}\) jet algorithm with a radius parameter \(R=0.4\) to cluster stable final-state particles with a decay length of \(c\tau >10\) mm, excluding muons and neutrinos. Jets are selected in the \(|y_\text {jet}|<3.4\) and \(p_{\text {T}} ^\text {jet}>25\,\text {Ge}\text {V}\) region. Jets within \(\Delta R = 0.4\) of electrons are rejected.

The closest reconstructed and particle-level jets are considered matched if \(\Delta R\) between their axes satisfies \(\Delta R < 0.4\).

The input for the unfolding is the transfer matrix, which maps reconstructed jets to the particle-level jets in the \(|y_\text {jet}|\)\(p_{\text {T}} ^\text {jet}\) plane, taking into account the bin-to-bin migrations that arise from limited detector resolution. An additional \(p_{\text {T}} ^\text {jet}\) bin, \(17\,\text {Ge}\text {V}<p_{\text {T}} ^\text {jet}<25\,\text {Ge}\text {V}\), is included in the reconstructed and particle-level jet spectra to account for the migrations from the low \(p_{\text {T}} ^\text {jet}\) range. This bin is not reported in the measurement.

Given the significant amount of migration between jet transverse momentum bins, the unfolding is performed in all \(|y_\text {jet}|\) and \(p_{\text {T}} ^\text {jet}\) bins simultaneously. The migration between adjacent \(|y_\text {jet}|\) bins is found to be small.

The transfer matrix is defined for matched jets. Therefore, the reconstructed jet spectrum must be corrected to account for matching efficiencies prior to unfolding. The reconstruction-level matching efficiency is calculated as the fraction of reconstructed jets matching particle-level jets. This efficiency is 80–90% in the \(25\,\text {Ge}\text {V}<p_{\text {T}} ^\text {jet}<100\,\text {Ge}\text {V}\) region and is above 99% in the \(p_{\text {T}} ^\text {jet}>100\,\text {Ge}\text {V}\) region. The particle-level jet matching efficiency is calculated as the fraction of particle-level jets matching reconstructed jets. This efficiency is 45–55% in all bins of the measurement due to the inefficiency of the \(Z\) boson reconstruction.

The backgrounds are subtracted from data prior to unfolding. The unfolded number of jets in data, \(N^{\mathcal {P}}_{i}\), in each bin i of the measurement is obtained as

$$\begin{aligned} N^{\mathcal {P}}_{i} = \frac{1}{\mathcal {E}^{\mathcal {P}}_{i}} \sum _{j} U_{ij} \mathcal {E}^{\mathcal {R}}_{j} N^{\mathcal {R}}_{j}, \end{aligned}$$
(1)

where \(N^{\mathcal {R}}_{j}\) is the number of jets reconstructed in bin j after the background subtraction, \(U_{ij}\) is an element of the unfolding matrix, and \(\mathcal {E}^{\mathcal {R}}_{j}\) and \(\mathcal {E}^{\mathcal {P}}_{i}\) are the reconstruction-level and particle-level jet matching efficiencies, respectively.

The transfer matrix and the matching efficiencies are improved using three unfolding iterations to reduce the impact of the particle-level jet spectra mis-modelling on the unfolded data.

8 Experimental uncertainties

8.1 Electron uncertainties

The electron energy scale has associated statistical uncertainties and systematic uncertainties arising from a possible bias in the calibration method, the choice of generator, the presampler energy scale, and imperfect knowledge of the material in front of the EM calorimeter. The total energy-scale uncertainty is calculated as the sum in quadrature of these components. It is varied by \(\pm 1\sigma \) in order to propagate the electron energy scale uncertainty into the measured \(Z \text { + jets}\) cross-sections.

The electron energy resolution has uncertainties associated to the extraction of the resolution difference between data and simulation using \(Z \rightarrow ee\) events, to the knowledge of the sampling term of the calorimeter energy resolution and to the pile-up noise modelling. These uncertainties are evaluated in situ using the \(Z \rightarrow ee\) events, and the total uncertainty is calculated as the sum in quadrature of the different uncertainties. The scale factor for electron energy resolution in MC simulation is varied by \(\pm 1\sigma \) in the total uncertainty in order to propagate this uncertainty into the \(Z \text { + jets}\) cross-section measurements.

The uncertainties in calculations of the electron trigger, reconstruction and identification efficiencies are propagated into the measurements by \(\pm 1\sigma \) variations of the scale factors, used to reweight the MC simulated events, within the total uncertainty of each efficiency [20, 51].

For each systematic variation a new transfer matrix and new matching efficiencies are calculated, and data unfolding is performed. The deviation from the nominal unfolded result is assigned as the systematic uncertainty in the measurements.

8.2 Jet uncertainties

The uncertainty in the jet energy measurement is described by 65 uncertainty components [57]. Of these, 56 JES uncertainty components are related to detector description, physics modelling and sample sizes of the \(Z/\gamma \) and multijet MC samples used for JES in situ measurements. The single-hadron response studies are used to describe the JES uncertainty in high-\(p_{\text {T}}\) jet regions, where the in situ studies have few events. The MC non-closure uncertainty takes into account the differences in the jet response due to different shower models used in MC generators. Four uncertainty components are due to the pile-up corrections of the jet kinematics, and take into account mis-modelling of \(N_\text {PV}\) and \(\langle \mu \rangle \) distributions, the average energy density and the residual \(p_{\text {T}} \) dependence of the \(N_\text {PV}\) and \(\langle \mu \rangle \) terms. Two flavour-based uncertainties take into account the difference between the calorimeter responses to the quark- and gluon-initiated jets. One uncertainty component describes the correction for the energy leakage beyond the calorimeter (‘punch-through’ effect). All JES uncertainties are treated as bin-to-bin correlated and independent of each other.

A reduced set of uncertainties, which combines the uncertainties of the in situ methods into six components with a small loss of correlation information, is used in this measurement. The JES uncertainties are propagated into the measurements in the same way as done for electron uncertainties.

The uncertainty that accounts for the difference in JVF requirement efficiency between data and MC simulation is evaluated by varying the nominal JVF requirement in MC simulation to represent a few percent change in efficiency [60]. The unfolding transfer matrix and the matching efficiencies are re-derived, and the results of varying the JVF requirement are propagated to the unfolded data. The deviations from the nominal results are used as the systematic uncertainty.

Pile-up jets are effectively suppressed by the selection requirements. The jet yields in events with low \(\langle \mu \rangle \) and high \(\langle \mu \rangle \) are compared with the jet yields in events without any requirements on \(\langle \mu \rangle \). These jet yields agree with each other within the statistical uncertainties. The same result is achieved by comparing the jet yields in events that have low or high numbers of reconstructed primary vertices with the jet yields in events from the nominal selection. Consequently, no additional pile-up uncertainty is introduced.

The jet energy resolution (JER) uncertainty accounts for the mis-modelling of the detector jet energy resolution by the MC simulation. To evaluate the JER uncertainty in the measured \(Z \text { + jets}\) cross-sections, the energy of each jet in MC simulation is smeared by a random number taken from a Gaussian distribution with a mean value of one and a width equal to the quadratic difference between the varied resolution and the nominal resolution [61]. The smearing is repeated 100 times and then averaged. The transfer matrix determined from the averaged smearing is used for unfolding. The result is compared with the nominal measurement and the symmetrised difference is used as the JER uncertainty.

The uncertainty that accounts for the mis-modelling of the ‘medium’ jet quality criteria is evaluated using jets, selected with the ‘loose’ and ‘tight’ criteria. The data-to-MC ratios of the reconstructed \(Z \text { + jets}\) distributions, obtained with different jet quality criteria, is compared with the nominal. An uncertainty of 1%, which takes the observed differences into account, is assigned to the measured \(Z \text { + jets}\) cross-section in all bins of \(|y_\text {jet}|\) and \(p_{\text {T}} ^\text {jet}\).

8.3 Background uncertainties

The uncertainties in each background estimation are propagated to the measured \(Z \text { + jets}\) cross-sections.

The data contamination by the \(Z {} \rightarrow \tau \tau \), diboson, \(t\bar{t}\) and single-top-quark backgrounds is estimated using simulated spectra scaled to the corresponding total cross-sections. Each of these background cross-sections has an uncertainty. The normalisation of each background is independently varied up and down by its uncertainty and propagated to the final result. The MC simulation of the dominant \(t\bar{t}\) background describes the shapes of the jet \(p_{\text {T}} ^\text {jet}\) and \(y_\text {jet}\) distributions in data to within a few percent [62], such that possible shape mis-modellings of the jet kinematics in \(t\bar{t}\) events are covered by the uncertainty in the total \(t\bar{t}\) cross-section. The shape mis-modellings in other backgrounds have negligible effect on the final results. Therefore, no dedicated uncertainties due to the background shape mis-modelling are assigned.

The uncertainties in the combined multijet and \(W \text {+ jets}\) background arise from assumptions about the template shape and normalisation. The shape of the template depends on the control region selection and the control region contaminations by the other backgrounds. The template normalisation depends on the \(m_{ee}\) range, used to fit the template to the measured \(Z \text { + jets}\) events, due to different amounts of background contamination in the tails of the \(m_{ee}\) distribution.

To evaluate the template shape uncertainty due to the control region selection, a different set of electron identification criteria is used to derive a modified template. The selection requires two reconstructed electrons with at least one electron that satisfies the ‘loose’ identification criteria, but not the ‘medium’ ones. The difference between the nominal and modified templates is used to create a symmetric template to provide up and down variations of this systematic uncertainty.

To estimate the template shape uncertainty due to the control region contaminations by the other backgrounds, the \(Z {} \rightarrow \tau \tau \), diboson, \(t\bar{t}\) and single-top-quark cross-sections are varied within their uncertainties. The dominant change in the template shape is due to \(t\bar{t}\) cross-section variation, while the contributions from the variation of the other background cross-sections are small. The templates varied within the \(t\bar{t}\) cross-section uncertainties are used to propagate this uncertainty into the measurement.

The uncertainty in the multijet and \(W \text {+ jets}\) background template normalisation to the measured \(Z \text { + jets}\) events is evaluated by fitting the template in the \(66\,\text {Ge}\text {V}<m_{ee}<140\,\text {Ge}\text {V}\) and \(60\,\text {Ge}\text {V}<m_{ee}<116\,\text {Ge}\text {V}\) regions, excluding the bins under the \(Z\) boson peak within the \(80\,\text {Ge}\text {V}<m_{ee}<100\,\text {Ge}\text {V}\) region, and in the \(60\,\text {Ge}\text {V}<m_{ee}<140\,\text {Ge}\text {V}\) region, excluding the bins within the \(70\,\text {Ge}\text {V}<m_{ee}<110\,\text {Ge}\text {V}\). As a result, the normalisation varies up and down depending on the number of background events in both tails of the \(m_{ee}\) distribution. The templates with the largest change in the normalisation are used to propagate this uncertainty into the measurement.

The data unfolding is repeated for each systematic variation of the backgrounds. The differences relative to the nominal \(Z \text { + jets}\) cross-section are used as the systematic uncertainties.

8.4 Unfolding uncertainty

The accuracy of the unfolding procedure depends on the quality of the description of the measured spectrum in the MC simulation used to build the unfolding matrix. Two effects are considered in order to estimate the influence of MC modelling on the unfolding results: the shape of the particle-level spectrum and the parton shower description.

The impact of the particle-level shape mis-modelling on the unfolding is estimated using a data-driven closure test. For this test, the particle-level (\(|y_\text {jet}|\)\(p_{\text {T}} ^\text {jet}\)) distribution in Sherpa is reweighted using the transfer matrix, such that the shape of the matched reconstructed (\(|y_\text {jet}|\)\(p_{\text {T}} ^\text {jet}\)) distribution agrees with the measured spectrum corrected for the matching efficiency. The reweighted reconstructed (\(|y_\text {jet}|\)\(p_{\text {T}} ^\text {jet}\)) distributions are then unfolded using the nominal Sherpa transfer matrix. The results are compared with the reweighted particle-level (\(|y_\text {jet}|\)\(p_{\text {T}} ^\text {jet}\)) spectrum and the relative differences are assigned as the uncertainty.

The impact of the differences in the parton shower description between Sherpa and Alpgen+Pythia on the unfolding results is estimated using the following test. The Alpgen+Pythia particle-level (\(|y_\text {jet}|\)\(p_{\text {T}} ^\text {jet}\)) spectrum is reweighted using the Alpgen+Pythia transfer matrix, such that its reconstruction-level distribution agrees with the one in Sherpa. The original reconstructed (\(|y_\text {jet}|\)\(p_{\text {T}} ^\text {jet}\)) distribution in Sherpa is then unfolded using the reweighted Alpgen+Pythia transfer matrix. The results are compared with the original particle-level (\(|y_\text {jet}|\)\(p_{\text {T}} ^\text {jet}\)) spectrum in Sherpa and the differences are assigned as the uncertainty.

Both unfolding uncertainties are symmetrised at the cross-section level.

8.5 Reduction of statistical fluctuations in systematic uncertainties

The systematic uncertainties suffer from fluctuations of a statistical nature.

The statistical components in the electron and jet uncertainties are estimated using toy MC simulations with 100 pseudo-experiments. Each \(Z \text { + jets}\) event in the systematically varied configurations is reweighted by a random number taken from a Poisson distribution with a mean value of one. As a result, 100 replicas of transfer matrix and matching efficiencies are created for a given systematic uncertainty variation, and are used to unfold the data. The replicas of unfolded spectra are then divided by the nominal \(Z \text { + jets}\) distributions to create an ensemble of systematic uncertainty spectra. The statistical component in the systematic uncertainties is calculated as the RMS across all replicas in an ensemble.

The pseudo-experiments are not performed for the JER systematic uncertainty. The statistical errors in the JER systematic uncertainty are calculated, considering the unfolded data in the nominal and JER varied configurations to be independent of each other.

Each component of the unfolding uncertainty is derived using 100 pseudo-experiments to calculate the statistical error.

To reduce the statistical fluctuations, the bins are combined iteratively starting from both the right and left sides of each systematic uncertainty spectrum until their significance satisfies \(\sigma >1.5\). The result with the most bins remaining is used as the systematic uncertainty. A Gaussian kernel is then applied to regain the fine binning and smooth out any additional statistical fluctuations.

If up and down systematic variations within a bin result in uncertainties with the same sign, then the smaller uncertainty is set to zero.

8.6 Statistical uncertainties

Statistical uncertainties are derived using toy MC simulations with 100 pseudo-experiments performed in both data and MC simulation. The data portion of the statistical uncertainty is evaluated by unfolding the replicas of the data using the nominal transfer matrix and matching efficiencies. The MC portion is calculated using the replicas of the transfer matrix and matching efficiencies to unfold the nominal data. To calculate the total statistical uncertainty in the measurement, the \(Z \text { + jets}\) distributions, obtained from pseudo-experiments drawn from the data yields, are unfolded using the transfer matrices and efficiency corrections, calculated using pseudo-experiments in the MC simulation. The covariance matrices between bins of the measurement are computed using the unfolded results. The total statistical uncertainties are calculated using the diagonal elements of the covariance matrices.

8.7 Summary of experimental uncertainties

The \(Z \text { + jets}\) cross-section measurement has 39 systematic uncertainty components. All systematic uncertainties are treated as being uncorrelated with each other and fully correlated among \(|y_\text {jet}|\) and \(p_{\text {T}} ^\text {jet}\) bins.

The systematic uncertainties in the electron energy scale, electron energy resolution, and electron trigger, reconstruction and identification efficiencies are found to be below 1%.

The JES in situ methods uncertainty is 2–5% in most bins of the measurement. The \(\eta \)-intercalibration uncertainty is below 1% in the \(|y_\text {jet}|<1.0\) and \(p_{\text {T}} ^\text {jet}<200\,\text {Ge}\text {V}\) regions, but it increases with \(|y_\text {jet}|\), reaching 6–14% in the most forward rapidity bins. The \(\eta \)-intercalibration uncertainty is below 1.5% for jets with \(p_{\text {T}} ^\text {jet}>200\,\text {Ge}\text {V}\). The flavour-based JES uncertainties are below 3%. The pile-up components of the JES uncertainty are 0.5–1.5%. Other components of the JES uncertainty are below 0.2%.

The JVF uncertainty is below 1%.

The JER is the dominant source of uncertainty in the \(Z \text { + jets}\) cross-section in the \(25\,\text {Ge}\text {V}< p_{\text {T}} ^\text {jet}< 50\,\text {Ge}\text {V}\) region with a 3–10% contribution. In the \(50\,\text {Ge}\text {V}<p_{\text {T}} ^\text {jet}<100\,\text {Ge}\text {V}\) region the JER uncertainty is 1–3%, and below 1% for jets with higher transverse momenta.

The jet quality uncertainty is set constant at 1%, as discussed in Sect. 8.2.

The unfolding uncertainty due to the shape of the particle-level spectrum is 2–5% in the first \(p_{\text {T}} ^\text {jet}\) bin, \(25\,\text {Ge}\text {V}< p_{\text {T}} ^\text {jet}< 50\,\text {Ge}\text {V}\). In the \(50\,\text {Ge}\text {V}< p_{\text {T}} ^\text {jet}< 200\,\text {Ge}\text {V}\) region, this uncertainty is about 1.5% for central jets below \(|y_\text {jet}|=2\), while for forward jets this uncertainty increases to 5%. In the \(p_{\text {T}} ^\text {jet}> 200\,\text {Ge}\text {V}\) region, this uncertainty is below 1.5%. The unfolding uncertainty due to the parton shower description is 0.7% in the \(400\,\text {Ge}\text {V}<p_{\text {T}} ^\text {jet}<1050\,\text {Ge}\text {V}\) region, while for jets with smaller transverse momenta this uncertainty is negligible.

The \(t\bar{t}\) background uncertainty is 0.02–0.6% in all bins of the measurement. The \(Z {} \rightarrow \tau \tau \), diboson and single-top-quark background uncertainties are below 0.05%.

The multijet and \(W \text {+ jets}\) background uncertainty is 0.1–1.2% depending on \(|y_\text {jet}|\) and \(p_{\text {T}} ^\text {jet}\). The uncertainty in the background template normalisation is asymmetric due to different background contributions in the tails of the \(m_{ee}\) distribution in the background normalisation evaluation procedure. This uncertainty is \(_{-0.4}^{+0.1}\)% in the low \(p_{\text {T}} ^\text {jet}\) bins, increasing to \(_{-1.2}^{+0.4}\)% in the high \(p_{\text {T}} ^\text {jet}\) bins. The uncertainty in the multijet and \(W \text {+ jets}\) background control region selection increases from 0.03% to 0.6% as a function of \(p_{\text {T}} ^\text {jet}\). The contribution of the \(t\bar{t}\) cross-section variation to the multijet and \(W \text {+ jets}\) background uncertainty is below 0.1%.

The statistical uncertainties are 0.5–4% in the \(p_{\text {T}} ^\text {jet}<100\,\text {Ge}\text {V}\) region, 2–14% in the \(100\,\text {Ge}\text {V}<p_{\text {T}} ^\text {jet}<300\,\text {Ge}\text {V}\) region, 8–39% in the \(300\,\text {Ge}\text {V}<p_{\text {T}} ^\text {jet}<400\,\text {Ge}\text {V}\) region and 11–18% in the last \(p_{\text {T}} ^\text {jet}\) bin, \(400\,\text {Ge}\text {V}<p_{\text {T}} ^\text {jet}<1050\,\text {Ge}\text {V}\). The smallest statistical uncertainty corresponds to central rapidity regions, while the largest uncertainty corresponds to forward rapidity regions.

The experimental uncertainties are shown in Fig. 2. The largest total systematic uncertainty of 7–12% is in the \(25\,\text {Ge}\text {V}<p_{\text {T}} ^\text {jet}<50\,\text {Ge}\text {V}\) region, where the uncertainty increases from central rapidity jets to the forward rapidity jets, and up to 15% for the forward rapidity jets in the \(100\,\text {Ge}\text {V}<p_{\text {T}} ^\text {jet}<200\,\text {Ge}\text {V}\) region. The total systematic uncertainty decreases with increasing \(p_{\text {T}} ^\text {jet}\). In the \(400\,\text {Ge}\text {V}<p_{\text {T}} ^\text {jet}<1050\,\text {Ge}\text {V}\) region the total systematic uncertainty is 2–5%. The luminosity uncertainty of 1.9% is not shown and not included in the total uncertainty and its components.

Fig. 2
figure 2

Experimental uncertainties in the measured double-differential \(Z \text { + jets}\) production cross-section as a function of \(|y_\text {jet}|\) in \(p_{\text {T}} ^\text {jet}\) bins. The jet energy scale, jet energy resolution, unfolding, ‘other’ and total systematic uncertainties are shown with different colours overlaid. The jet energy scale uncertainty is the sum in quadrature of the jet energy scale uncertainty components. The unfolding uncertainty is the sum in quadrature of two unfolding uncertainties. The ‘other’ systematic uncertainty is the sum in quadrature of the electron uncertainties, background uncertainties, JVF and jet quality uncertainties. The total systematic uncertainty is the sum in quadrature of all systematic uncertainty components except for the luminosity uncertainty of 1.9%. The total statistical uncertainties are shown with vertical error bars

9 Fixed-order predictions and theoretical uncertainties

9.1 Fixed-order calculations

Theoretical \(Z \text { + jets}\) predictions at NLO are calculated using MCFM [40] interfaced to APPLgrid [63] for fast convolution with PDFs. The renormalisation and factorisation scales, \(\mu _\text {R}\) and \(\mu _\text {F}\), are set to

$$\begin{aligned} \mu _\text {R}=\mu _\text {F}=\frac{\sqrt{m_{ee}^2 + p_{\text {T}, Z}^2}+\sum p_\text {T, partons}}{2}, \end{aligned}$$

where \(m_{ee}\) is the electron pair’s invariant mass, \(p_{\text {T}, Z}\) is the transverse momentum of the \(Z\) boson and \(\sum p_\text {T, partons}\) is the sum of the transverse momenta of the partons.

The NLO \(Z \text { + jets}\) predictions are obtained using the CT14 NLO [64], NNPDF3.1 [65], JR14 NLO [66], HERAPDF2.0 [67], MMHT2014 [68], ABMP16 [69] and ATLAS-epWZ16 [70] PDF sets. The PDFs are determined by various groups using the experimental data and are provided with the uncertainties. The PDF uncertainties in the \(Z \text { + jets}\) cross-sections are calculated at the 68% confidence level according to the prescription recommended by the PDF4LHC group [71].

The following variations of factorisation and renormalisation scales are performed to assess the uncertainty due to missing higher-order terms: \(\{\mu _\text {R}/2, \mu _\text {F}\}\), \(\{2\mu _\text {R}, \mu _\text {F}\}\), \(\{\mu _\text {R}, \mu _\text {F}/2\}\), \(\{\mu _\text {R}, 2\mu _\text {F}\}\), \(\{\mu _\text {R}/2, \mu _\text {F}/2\}\), \(\{2\mu _\text {R}, 2\mu _\text {F}\}\). The envelope of the cross-sections calculated with different scales is used as the uncertainty.

The uncertainty due to the strong coupling is estimated using additional PDF sets, calculated with \(\alpha _{\text {S}} (m_Z^2)=0.116\) and \(\alpha _{\text {S}} (m_Z^2)=0.120\). The resulting uncertainty is scaled to the uncertainty of the world average \(\alpha _{\text {S}} (m_Z^2)= 0.118 \pm 0.0012\), as recommended by the PDF4LHC group [71].

The state-of-the-art NNLO \(Z \text { + jets}\) cross-section is calculated by the authors of Ref. [16] using NNLOJET [72]. The NNLO predictions are convolved with the CT14 PDF. The renormalisation and factorisation scales are set similarly to those in NLO calculations.

9.2 Non-perturbative correction

The fixed-order predictions are obtained at the parton level. Bringing fixed-order predictions to the particle level for comparisons with the measured \(Z \text { + jets}\) cross-sections requires a non-perturbative correction (NPC) that accounts for both the hadronisation and underlying-event effects.

The NPCs are studied using several MC generators to account for differences in the modelling of hadronisation and the underlying event. The studies are done using the leading-logarithm parton shower MC generators Pythia v. 8.210 with the A14 [73] underlying-event tune and Herwig++ v. 2.7.1 with the UE-EE5 tune [74], and the multi-leg matrix element MC generators Sherpa v. 1.4.5 with the CT10 PDF, Sherpa v. 2.2.0 with NNPDF 2.3 [75] and MadGraph v. 2.2.3 [76], supplemented with parton showers from Pythia v. 8.210 with the A14 tune.

The NPCs are calculated using the ratios of \(Z \text { + jets}\) cross-sections obtained at the particle level to those at the parton level. The correction derived using Sherpa v. 1.4.5 is the nominal one in this analysis. The envelope of the non-perturbative corrections, calculated with other MC generators, is used as the systematic uncertainty. The NPCs in different MC generators are shown in Fig. 3. The nominal correction for jets with low transverse momenta, \(25\,\text {Ge}\text {V}<p_{\text {T}} ^\text {jet}<50\,\text {Ge}\text {V}\), in the central rapidity regions, \(|y_\text {jet}|<1.5\), is small, but it increases to 5% in the forward rapidity bins. The nominal correction for jets with higher transverse momenta is below 2%. These corrections together with uncertainties are provided in HEPData [77].

9.3 QED radiation correction

The fixed-order \(Z \text { + jets}\) cross-section predictions must be corrected for the QED radiation in order to be compared with data. The correction is determined as the ratio of two \(Z \text { + jets}\) cross-sections, one calculated using dressed electrons after QED final-state radiation (FSR), with all photons clustered within a cone of \(\Delta R = 0.1\) around the electron axis, and the other calculated using Born-level electrons at the lowest order in the electromagnetic coupling \(\alpha _\text {QED}\) prior to QED FSR radiation. The baseline correction is calculated using the Sherpa MC samples, while the correction calculated using Alpgen+Pythia is used to estimate the uncertainty. The uncertainty is calculated as the width of the envelope of corrections obtained with these two MC generators. The results are shown in Fig. 4. The QED correction is largest in the \(25\,\text {Ge}\text {V}<p_{\text {T}} ^\text {jet}<50\,\text {Ge}\text {V}\) region. It is about 5% for jets in the central absolute rapidity regions. In the \(p_{\text {T}} ^\text {jet}>50\,\text {Ge}\text {V}\) regions the QED correction is 1.5–2.5%, decreasing as a function of jet transverse momentum. The QED corrections calculated using Alpgen+Pythia are in good agreement with those from Sherpa. These corrections together with uncertainties are provided in HEPData.

9.4 Summary of theoretical uncertainties

The total theoretical uncertainties are calculated as the sum in quadrature of the effects of the PDF, scale, and \(\alpha _{\text {S}}\) uncertainties, and the uncertainties due to non-perturbative and QED radiation effects.

The uncertainties for the \(Z \text { + jets}\) cross-section calculated at NLO using the CT14 PDF as a function of \(|y_\text {jet}|\) in \(p_{\text {T}} ^\text {jet}\) bins are shown in Fig. 5. The total uncertainties are dominated by the scale and NPC uncertainties in the \(p_{\text {T}} ^\text {jet}<100\,\text {Ge}\text {V}\) region, where they reach \(\pm 15\%\) and \(-10\%\), respectively. In the \(p_{\text {T}} ^\text {jet}>100\,\text {Ge}\text {V}\) region, the scale uncertainty alone dominates, as the NPC uncertainty decreases for high jet transverse momenta. The total uncertainty in this region is 10–20%. Other uncertainties are below 5%.

The NNLO uncertainties are shown in Fig. 6. The scale uncertainty at NNLO is significantly reduced. This uncertainty is below 1% in the \(25\,\text {Ge}\text {V}<p_{\text {T}} ^\text {jet}<50\,\text {Ge}\text {V}\) bin, increasing to 5% in the \(400\,\text {Ge}\text {V}<p_{\text {T}} ^\text {jet}<1050\,\text {Ge}\text {V}\) bin. In the \(p_{\text {T}} ^\text {jet}<200\,\text {Ge}\text {V}\) region, the negative part of the total uncertainty is dominated by the NPC uncertainty and its absolute value reaches 7–15% depending on the jet rapidity. The positive part of the total uncertainty is within 5%, with about equal contributions from PDF, scale and \(\alpha _{\text {S}}\) uncertainties. In the \(p_{\text {T}} ^\text {jet}>200\,\text {Ge}\text {V}\) region, both the negative and positive parts of the total uncertainty are within 6% in most bins.

The uncertainty in the QED correction is below 0.5% and is negligible in the fixed-order theory predictions.

Fig. 3
figure 3

The non-perturbative correction for the \(Z \text { + jets}\) production cross-section as a function of \(|y_\text {jet}|\) in \(p_{\text {T}} ^\text {jet}\) bins. The spread of predictions represents the uncertainty

Fig. 4
figure 4

The correction for QED radiation effects for the \(Z \text { + jets}\) production cross-section as a function of \(|y_\text {jet}|\) in \(p_{\text {T}} ^\text {jet}\) bins. The spread of predictions represents the uncertainty

Fig. 5
figure 5

The uncertainties in NLO pQCD predictions as a function of \(|y_\text {jet}|\) in \(p_{\text {T}} ^\text {jet}\) bins. Total pQCD uncertainty is the sum in quadrature of the PDF, scale and \(\alpha _{\text {S}}\) uncertainties. Total theory uncertainty is the sum in quadrature of the total pQCD uncertainty and the uncertainties from the non-perturbative and QED radiation corrections. The CT14 PDF is used in the calculations

Fig. 6
figure 6

The uncertainties in NNLO pQCD predictions as a function of \(|y_\text {jet}|\) in \(p_{\text {T}} ^\text {jet}\) bins. Total pQCD uncertainty is the sum in quadrature of the PDF, scale and \(\alpha _{\text {S}}\) uncertainties. Total theory uncertainty is the sum in quadrature of the total pQCD uncertainty and the uncertainties from the non-perturbative and QED radiation corrections. The CT14 PDF is used in the calculations

10 Results

The double-differential \(Z \text { + jets}\) cross-section as a function of \(|y_\text {jet}|\) and \(p_{\text {T}} ^\text {jet}\) is calculated as

$$\begin{aligned} \frac{\text {d}^{2} \sigma }{\text {d}p_{\text {T}} ^\text {jet}\text {d}|y_\text {jet}|} = \frac{1}{ \mathcal {L}} \frac{N^{\mathcal {P}}_{i}}{\Delta p_{\text {T}} ^\text {jet}{} \Delta |y_\text {jet}|{}}, \end{aligned}$$

where \(\mathcal {L}\) is the integrated luminosity, \(N^{\mathcal {P}}_{i} \) is the number of jets in data at the particle level as given in Eq. (1), and \(\Delta p_{\text {T}} ^\text {jet}\) and \(\Delta |y_\text {jet}|\) are the widths of the jet transverse momentum and absolute jet rapidity ranges for bin i, respectively. The backgrounds are subtracted before data unfolding is performed to obtain \(N^{\mathcal {P}}_{i} \).

The measured \(Z \text { + jets}\) cross-section covers five orders of magnitude and falls steeply as a function of \(|y_\text {jet}|\) and \(p_{\text {T}} ^\text {jet}\). A summary of measured cross-sections, together with the systematic and statistical uncertainties, is provided in Section A. The measured cross-sections with the full breakdown of all uncertainties are provided in HEPData.

The comparisons with the theoretical predictions are shown in Figs. 7, 8, 9, 10, 11 and 12. The fixed-order theoretical predictions are corrected for the non-perturbative and QED radiation effects. The NLO predictions are lower than the data by approximately 5–10%. However, this difference is covered by the uncertainties. The NNLO calculations compensate for the NLO-to-data differences in most bins of the measurement and show better agreement with the central values of the cross-sections in data. The Sherpa v. 1.4 and Alpgen+Pythia MC-to-data ratios are approximately constant across all \(|y_\text {jet}|\) bins, but a dependence on \(p_{\text {T}} ^\text {jet}\) is observed. The Sherpa v. 1.4 predictions are lower than the data by about 10% in the \(25\,\text {Ge}\text {V}< p_{\text {T}} ^\text {jet}< 200\,\text {Ge}\text {V}\) region, but in the \(p_{\text {T}} ^\text {jet}> 200\,\text {Ge}\text {V}\) region they agree within a few percent. The Alpgen+Pythia predictions agree with data in the \(25\,\text {Ge}\text {V}< p_{\text {T}} ^\text {jet}< 100\,\text {Ge}\text {V}\) region, but exceed the data as a function of \(p_{\text {T}} ^\text {jet}\), the largest difference being about 20% in the highest \(p_{\text {T}} ^\text {jet}\) bin, \(400\,\text {Ge}\text {V}< p_{\text {T}} ^\text {jet}< 1050\,\text {Ge}\text {V}\).

Additionally, data is compared to the Sherpa v. 2.2 prediction. In this prediction, the matrix elements are calculated with NLO accuracy for the inclusive Z production process up to two additional partons in the final state, and with LO accuracy in the final states with up to four partons. Sherpa v. 2.2 MEs are convolved with the NNPDF 3.0 [65] PDFs. The MEs are merged with Sherpa parton shower using the ME+PS@NLO [78] prescription. This prediction shows a good agreement with data in all bins of the measurement.

The ratios between the measured \(Z \text { + jets}\) production cross-sections and the NLO predictions, calculated with various PDF sets, are shown in Figs. 13, 14 and 15. The calculations with MMHT2014 and NNPDF3.1 predict 1–2% larger cross-sections compared to those using the CT14 PDF. The cross-sections calculated with ATLAS-epWZ16 PDF are larger by 2–3%. The ABMP16 and HERAPDF2.0 cross-section predictions in the \(|y_\text {jet}|<2.0\) and \(p_{\text {T}} ^\text {jet}<100\,\text {Ge}\text {V}\) regions are 3–5% larger than those from the CT14 PDF, while in other bins of the measurement their predictions are up to 5% lower than those obtained with the CT14 PDF. The JR14 PDF predictions are 2–5% lower than those from the CT14 PDF in the \(25\,\text {Ge}\text {V}< p_{\text {T}} ^\text {jet}< 200\,\text {Ge}\text {V}\) region and higher by 2% in the \(p_{\text {T}} ^\text {jet}> 200\,\text {Ge}\text {V}\) region. The differences between the cross-sections calculated at NLO accuracy with various PDF sets are covered by the theoretical uncertainties. In the NNLO calculations, the difference between CT14 PDF and NNPDF3.1 predictions is 2–5%, which is comparable to the size of the theoretical uncertainties, as shown in Fig. 16.

Fig. 7
figure 7

The double-differential \(Z \text { + jets}\) production cross-section as a function of \(|y_\text {jet}|\) in the \(25\,\text {Ge}\text {V}< p_{\text {T}} ^\text {jet}< 50\,\text {Ge}\text {V}\) range. The data are compared with the Sherpa v. 1.4, Sherpa v. 2.2 and Alpgen+Pythia parton shower MC generator predictions and with the fixed-order theory predictions. The fixed-order theory predictions are corrected for the non-perturbative and QED radiation effects. The fixed-order calculations are performed using the CT14 PDF. The total statistical uncertainties are shown with error bars. The total uncertainties in the measurement and in the fixed-order theory predictions are represented with shaded bands. The total uncertainty in the measurement is the sum in quadrature of the statistical and systematic uncertainties except for the luminosity uncertainty of 1.9%. The total uncertainty in the fixed-order theory predictions is the sum in quadrature of the effects of the PDF, scale, and \(\alpha _{\text {S}}\) uncertainties, and the uncertainties from the non-perturbative and QED radiation corrections. Lower panels show the ratios of predictions to data

Fig. 8
figure 8

The double-differential \(Z \text { + jets}\) production cross-section as a function of \(|y_\text {jet}|\) in the \(50\,\text {Ge}\text {V}< p_{\text {T}} ^\text {jet}< 100\,\text {Ge}\text {V}\) range. The data are compared with the Sherpa v. 1.4, Sherpa v. 2.2 and Alpgen+Pythia parton shower MC generator predictions and with the fixed-order theory predictions. The fixed-order theory predictions are corrected for the non-perturbative and QED radiation effects. The fixed-order calculations are performed using the CT14 PDF. The total statistical uncertainties are shown with error bars. The total uncertainties in the measurement and in the fixed-order theory predictions are represented with shaded bands. The total uncertainty in the measurement is the sum in quadrature of the statistical and systematic uncertainties except for the luminosity uncertainty of 1.9%. The total uncertainty in the fixed-order theory predictions is the sum in quadrature of the effects of the PDF, scale, and \(\alpha _{\text {S}}\) uncertainties, and the uncertainties from the non-perturbative and QED radiation corrections. Lower panels show the ratios of predictions to data

Fig. 9
figure 9

The double-differential \(Z \text { + jets}\) production cross-section as a function of \(|y_\text {jet}|\) in the \(100\,\text {Ge}\text {V}< p_{\text {T}} ^\text {jet}< 200\,\text {Ge}\text {V}\) range. The data are compared with the Sherpa v. 1.4, Sherpa v. 2.2 and Alpgen+Pythia parton shower MC generator predictions and with the fixed-order theory predictions. The fixed-order theory predictions are corrected for the non-perturbative and QED radiation effects. The fixed-order calculations are performed using the CT14 PDF. The total statistical uncertainties are shown with error bars. The total uncertainties in the measurement and in the fixed-order theory predictions are represented with shaded bands. The total uncertainty in the measurement is the sum in quadrature of the statistical and systematic uncertainties except for the luminosity uncertainty of 1.9%. The total uncertainty in the fixed-order theory predictions is the sum in quadrature of the effects of the PDF, scale, and \(\alpha _{\text {S}}\) uncertainties, and the uncertainties from the non-perturbative and QED radiation corrections. Lower panels show the ratios of predictions to data

Fig. 10
figure 10

The double-differential \(Z \text { + jets}\) production cross-section as a function of \(|y_\text {jet}|\) in the \(200\,\text {Ge}\text {V}< p_{\text {T}} ^\text {jet}< 300\,\text {Ge}\text {V}\) range. The data are compared with the Sherpa v. 1.4, Sherpa v. 2.2 and Alpgen+Pythia parton shower MC generator predictions and with the fixed-order theory predictions. The fixed-order theory predictions are corrected for the non-perturbative and QED radiation effects. The fixed-order calculations are performed using the CT14 PDF. The total statistical uncertainties are shown with error bars. The total uncertainties in the measurement and in the fixed-order theory predictions are represented with shaded bands. The total uncertainty in the measurement is the sum in quadrature of the statistical and systematic uncertainties except for the luminosity uncertainty of 1.9%. The total uncertainty in the fixed-order theory predictions is the sum in quadrature of the effects of the PDF, scale, and \(\alpha _{\text {S}}\) uncertainties, and the uncertainties from the non-perturbative and QED radiation corrections. Lower panels show the ratios of predictions to data

Fig. 11
figure 11

The double-differential \(Z \text { + jets}\) production cross-section as a function of \(|y_\text {jet}|\) in the \(300\,\text {Ge}\text {V}< p_{\text {T}} ^\text {jet}< 400\,\text {Ge}\text {V}\) range. The data are compared with the Sherpa v. 1.4, Sherpa v. 2.2 and Alpgen+Pythia parton shower MC generator predictions and with the fixed-order theory predictions. The fixed-order theory predictions are corrected for the non-perturbative and QED radiation effects. The fixed-order calculations are performed using the CT14 PDF. The total statistical uncertainties are shown with error bars. The total uncertainties in the measurement and in the fixed-order theory predictions are represented with shaded bands. The total uncertainty in the measurement is the sum in quadrature of the statistical and systematic uncertainties except for the luminosity uncertainty of 1.9%. The total uncertainty in the fixed-order theory predictions is the sum in quadrature of the effects of the PDF, scale, and \(\alpha _{\text {S}}\) uncertainties, and the uncertainties from the non-perturbative and QED radiation corrections. Lower panels show the ratios of predictions to data

Fig. 12
figure 12

The double-differential \(Z \text { + jets}\) production cross-section as a function of \(|y_\text {jet}|\) in the \(400\,\text {Ge}\text {V}< p_{\text {T}} ^\text {jet}< 1050\,\text {Ge}\text {V}\) range. The data are compared with the Sherpa v. 1.4, Sherpa v. 2.2 and Alpgen+Pythia parton shower MC generator predictions and with the fixed-order theory predictions. The fixed-order theory predictions are corrected for the non-perturbative and QED radiation effects. The fixed-order calculations are performed using the CT14 PDF. The total statistical uncertainties are shown with error bars. The total uncertainties in the measurement and in the fixed-order theory predictions are represented with shaded bands. The total uncertainty in the measurement is the sum in quadrature of the statistical and systematic uncertainties except for the luminosity uncertainty of 1.9%. The total uncertainty in the fixed-order theory predictions is the sum in quadrature of the effects of the PDF, scale, and \(\alpha _{\text {S}}\) uncertainties, and the uncertainties from the non-perturbative and QED radiation corrections. Lower panels show the ratios of predictions to data

Fig. 13
figure 13

Ratio of the measured \(Z \text { + jets}\) production cross-section and the NLO QCD predictions, obtained using MCFM, corrected for the non-perturbative and QED radiation effects as a function of \(|y_\text {jet}|\) and \(p_{\text {T}} ^\text {jet}\) bins. Theoretical predictions are calculated using various PDF sets. The coloured error bars represent the sum in quadrature of the effects of the PDF, scale, and \(\alpha _{\text {S}}\) uncertainties, and the uncertainties from the non-perturbative and QED radiation corrections. The grey band shows the sum in quadrature of the statistical and systematic uncertainties in the measurement except for the luminosity uncertainty of 1.9%

Fig. 14
figure 14

Ratio of the measured \(Z \text { + jets}\) production cross-section and the NLO QCD predictions, obtained using MCFM, corrected for the non-perturbative and QED radiation effects as a function of \(|y_\text {jet}|\) and \(p_{\text {T}} ^\text {jet}\) bins. Theoretical predictions are calculated using various PDF sets. The coloured error bars represent the sum in quadrature of the effects of the PDF, scale, and \(\alpha _{\text {S}}\) uncertainties, and the uncertainties from the non-perturbative and QED radiation corrections. The grey band shows the sum in quadrature of the statistical and systematic uncertainties in the measurement except for the luminosity uncertainty of 1.9%

Fig. 15
figure 15

Ratio of the measured \(Z \text { + jets}\) production cross-section and the NLO QCD predictions, obtained using MCFM, corrected for the non-perturbative and QED radiation effects as a function of \(|y_\text {jet}|\) and \(p_{\text {T}} ^\text {jet}\) bins. Theoretical predictions are calculated using various PDF sets. The coloured error bars represent the sum in quadrature of the effects of the PDF, scale, and \(\alpha _{\text {S}}\) uncertainties, and the uncertainties from the non-perturbative and QED radiation corrections. The grey band shows the sum in quadrature of the statistical and systematic uncertainties in the measurement except for the luminosity uncertainty of 1.9%

Fig. 16
figure 16

Ratio of the measured \(Z \text { + jets}\) production cross-section and the NNLO QCD predictions, obtained using NNLOJET, corrected for the non-perturbative and QED radiation effects as a function of \(|y_\text {jet}|\) and \(p_{\text {T}} ^\text {jet}\) bins. Theoretical predictions are calculated using various PDF sets. The coloured error bars represent the sum in quadrature of the effects of the PDF, scale, and \(\alpha _{\text {S}}\) uncertainties, and the uncertainties from the non-perturbative and QED radiation corrections. The grey band shows the sum in quadrature of the statistical and systematic uncertainties in the measurement except for the luminosity uncertainty of 1.9%

11 Quantitative data and theory comparison

The fixed-order pQCD predictions at NNLO accuracy, corrected for electroweak and non-perturbative effects, are quantitatively compared with the measured cross-section using a \(\chi ^2\) function that accounts for both the experimental and theoretical uncertainties

$$\begin{aligned}&\chi ^2\left( \mathbf {\beta ^\text {data}},\mathbf {\beta ^\text {theory}}\right) \nonumber \\&\quad = \sum \limits _{i=1}^{N_\text {bins}} \frac{\left( \sigma ^\text {data}_i + \sum _{\mu }\Gamma _{i;\mu }^\text {data}\beta _\mu ^\text {data} - \sigma ^\text {theory}_i - \sum _{\nu }\Gamma _{i;\nu }^\text {theory}\beta _\nu ^\text {theory} \right) ^2}{\Delta _i^2}\nonumber \\&\quad + \sum \limits _{\mu } \left( \beta _\mu ^\text {data}\right) ^2 + \sum \limits _{\nu } \left( \beta _\nu ^\text {theory}\right) ^2, \end{aligned}$$
(2)

where experimental (theoretical) uncertainties are included in the calculation using the nuisance parameter vectors \(\mathbf {\beta ^\text {data}}\left( \mathbf {\beta ^\text {theory}}\right) \) and their influence on the data and the theory predictions is described by the respective \(\Gamma _{l;\rho }\) matrices. The Latin indices run over bins of measurements and the Greek indices render sources of uncertainties. The measured cross-sections and their theory predictions in each bin are represented by \(\sigma ^\text {data}_i\) and \(\sigma ^\text {theory}_i\), respectively. Uncorrelated uncertainties in data are denoted by \(\Delta _i\). The theoretical uncertainties include those arising from renormalisation and factorisation scales variations, PDF uncertainties, uncertainties in calculations of non-perturbative and electroweak effects as well as from the \(\alpha _{\text {S}} (m_Z)\) uncertainty. All experimental and theoretical systematic uncertainties are assumed to be independent of each other, and fully correlated across the bins of the measurement. The negligible correlations of statistical uncertainties are not included in the \(\chi ^2\) tests presented here.

The minimisation of Eq. (2), for the case of symmetric systematic uncertainties, results in a system of linear equations for the shifts of systematic uncertainties, \(\beta _\rho \). The asymmetries in systematic uncertainties are accounted for using an iterative procedure. Here, the influences \(\Gamma _{l;\rho } \) are recalculated as

$$\begin{aligned} \Gamma _{l;\rho } \rightarrow {\Gamma }_{l;\rho } + {{\Omega }}_{l;\rho } \beta _{\rho }, \end{aligned}$$

where \(\Gamma _{l;\rho } = \frac{1}{2} \left( \Gamma _{l;\rho }^+ - \Gamma _{l;\rho }^-\right) \) and \( \Omega _{l;\rho } = \frac{1}{2} \left( \Gamma _{l;\rho }^+ + \Gamma _{l;\rho }^-\right) \), after each iteration using the shifts \(\beta _\rho \) from the previous iteration. The \(\Gamma _{l;\rho }^+\) and \(\Gamma _{l;\rho }^-\) are positive and negative components of systematic uncertainties, respectively. The \(\chi ^2\) values at the minimum provide a measure of the probability of compatibility between the measurements and the predictions.

Table 1 shows a summary of the calculated \(\chi ^2_\text {uncorr}\), the first term in Eq. (2), together with \(\chi ^2_\text {corr}\), the sum of squared shifts of nuisance parameters, for each \(p_{\text {T}} ^\text {jet}\) bin separately. A good agreement between measurements and theory is seen for the fits in individual \(p_{\text {T}} ^\text {jet}\) bins in the \(p_{\text {T}} ^\text {jet}>50\) GeV range, with not so good agreement in the \(25<p_{\text {T}} ^\text {jet}<50\) GeV range. The level of agreement between data and predictions is very similar for different PDF sets.

Table 1 Values of \(\chi ^2_\text {uncorr}\) and \(\chi ^2_\text {corr}\) evaluated for theory predictions corrected for non-perturbative and electroweak effects and measured \(Z \text { + jets}\) cross-sections. The total \(\chi ^2\) is equal to the sum of \(\chi ^2_\text {uncorr}\) and \(\chi ^2_\text {corr}\). The fits are performed individually in each \(p_{\text {T}} ^\text {jet}\) bin. The predictions are calculated using several NNLO QCD PDF sets and one NLO QCD PDF set, CT14nlo

In addition to fits of the predictions to measured cross-sections in the individual \(p_{\text {T}} ^\text {jet}\) bins, all measured data points are fitted simultaneously. Several ranges of \(p_{\text {T}} ^\text {jet}\) are considered. The results of the global fits are presented in Table 2. Very good agreement between measurement and calculation is observed when using the \(p_{\text {T}} ^\text {jet}>50\) GeV bins, while not so good agreement is observed when the \(25<p_{\text {T}} ^\text {jet}<50\) GeV bin is included in the global fit.

The results of the \(\chi ^2\) tests strongly depend on what is assumed about the correlation of systematic uncertainties. In general, the correlations of uncertainties related to detector measurements are carefully studied and well known [57, 79]. In contrast, the assumption of 100% correlations of uncertainties resulting from simple comparisons of two (or more) different MC simulations (two-point systematic uncertainties) are less justified. In order to investigate the impact of these assumptions on the results of \(\chi ^2\) tests performed in this section, the uncertainties that are derived from comparisons of two different MC models, namely uncertainties in the jet flavour composition and jet flavour response, were split into two subcomponents [80, 81]. The first subcomponent is derived by multiplying the original nuisance parameter with a linear function of \(p_{\text {T}} ^\text {jet}\) and jet absolute rapidity and the second subcomponent is constructed such that the sum in quadrature of both subcomponents is equal to the original nuisance parameter. These decorrelations did not result in a large improvement in the \(\chi ^2\) values.

Table 2 Values of \(\chi ^2\) evaluated from the comparison of theory predictions corrected for non-perturbative and electroweak effects with the measured \(Z \text { + jets}\) cross-sections. The fits are performed globally in all bins of the measurement within several \(p_{\text {T}} ^\text {jet}\) ranges. The predictions are calculated using several NNLO QCD PDF sets and one NLO QCD PDF set, CT14nlo

12 Conclusions

The double-differential \(Z \text { + jets}\) cross-section, with the \(Z\) boson decaying into an electron–positron pair, is measured using proton–proton collision data with an integrated luminosity \(19.9~\text{ fb }^{-1}\) collected by the ATLAS experiment at the LHC in 2012 at \(\sqrt{s}=8 \,\text {Te}\text {V}\) centre-of-mass energy. The measurement is performed as a function of the absolute jet rapidity and the jet transverse momentum.

The measured cross-section is corrected for detector effects and the results are provided at the particle level. The measurements are compared with theory predictions, calculated using the multi-leg matrix element MC generators Sherpa and Alpgen+Pythia, supplemented with parton shower simulations. Sherpa v. 1.4 and Alpgen+Pythia describe well the shape of the \(Z \text { + jets}\) distribution as a function of \(|y_\text {jet}|\), but not so well as a function of \(p_{\text {T}} ^\text {jet}\). Sherpa v. 2.2 is in good agreement with data in all bins of the measurement.

The parton-level fixed-order \(Z \text { + jets}\) predictions, corrected for hadronisation, underlying-event and QED radiation effects, agree with the data within the uncertainties.

The uncertainties in the measured cross-sections are about half of the theoretical uncertainties in the NLO calculations in most bins of the measurement and are approximately similar to the uncertainties in the NNLO calculations.

The measured double-differential \(Z \text { + jets}\) cross-section provides a precision input to constrain the parton distribution functions.