Measurements of the production cross section of a $Z$ boson in association with jets in pp collisions at $\sqrt{s} = 13$ TeV with the ATLAS detector

Abstract

Measurements of the production cross section of a $Z$ boson in association with jets in proton–proton collisions at $\sqrt{s} = 13$ TeV are presented, using data corresponding to an integrated luminosity of 3.16 fb$^{-1}$ collected by the ATLAS experiment at the CERN Large Hadron Collider in 2015. Inclusive and differential cross sections are measured for events containing a $Z$ boson decaying to electrons or muons and produced in association with up to seven jets with $p_{\text {T}} > 30$ GeV and $|y| <2.5$. Predictions from different Monte Carlo generators based on leading-order and next-to-leading-order matrix elements for up to two additional partons interfaced with parton shower and fixed-order predictions at next-to-leading order and next-to-next-to-leading order are compared with the measured cross sections. Good agreement within the uncertainties is observed for most of the modelled quantities, in particular with the generators which use next-to-leading-order matrix elements and the more recent next-to-next-to-leading-order fixed-order predictions.

1 Introduction

The measurement of the production of a $Z$ boson¹ in association with jets, $Z+\text {jets}$, constitutes a powerful test of perturbative quantum chromodynamics (QCD) [1, 2]. The large production cross section and easily identifiable decays of the $Z$ boson to charged leptonic final states offer clean experimental signatures which can be precisely measured. Such processes also constitute a non-negligible background for studies of the Higgs boson and in searches for new phenomena; typically in these studies, the multiplicity and kinematics of the jets are exploited to achieve a separation of the signal of interest from the Standard Model (SM) $Z+\text {jets}$ process. These quantities are often measured in control regions and subsequently extrapolated to the signal region with the use of Monte Carlo (MC) generators, which are themselves subject to systematic uncertainty and must be tuned and validated using data. Predictions from the most recent generators combine next-to-leading-order (NLO) multi-leg matrix elements with a parton shower (PS) and a hadronisation model. Fixed-order parton-level predictions for $Z+\text {jets}$ production at next-to-next-to-leading order (NNLO) are also available [3,4,5,6].

The $Z+\text {jets}$ production differential cross section was previously measured by the ATLAS [7], CMS [8], and LHCb [9] collaborations at the CERN Large Hadron Collider (LHC) [10] at centre-of-mass energies of $\sqrt{s}=7$ TeV [11,12,13,14,15] and 8 TeV [16,17,18], and by the CDF and D0 collaborations at the Tevatron collider at $\sqrt{s} = 1.96$ TeV [19, 20]. In this paper, proton–proton (pp) collision data corresponding to an integrated luminosity of 3.16 fb$^{-1}$, collected at $\sqrt{s} = 13$ TeV with the ATLAS detector during 2015, are used for measurements of the $Z$-boson production cross section in association with up to seven jets within a fiducial region defined by the detector acceptance. The $Z$ boson is identified using its decays to electron or muon pairs ($Z \rightarrow e^{+}e^{-}$, $Z \rightarrow \mu ^{+}\mu ^{-}$). Cross sections are measured separately for these two channels, and for their combination, as a function of the inclusive and exclusive jet multiplicity $N_\mathrm {jets}$ and the ratio of successive inclusive jet multiplicities $(N_\mathrm {jets}+1)/N_\mathrm {jets}$, the transverse momentum of the leading jet $p_{\text {T}}^{\text {jet}}$ for several jet multiplicities, the jet rapidity $y^{\text {jet}}$, the azimuthal separation between the two leading jets $\Delta \phi _{\text {jj}}$, the two leading jet invariant mass $m_{\text {jj}}$, and the scalar sum $H_{\text {T}}$ of the transverse momenta of all selected leptons and jets.

The paper is organised as follows. Section 2 contains a brief description of the ATLAS detector. The data and simulated samples as well as the $Z+\text {jets}$ predictions used in the analysis are described in Sect. 3. The event selection and its associated uncertainties are presented in Sect. 4, while the methods employed to estimate the backgrounds are shown in Sect. 5. Comparisons between data and Monte Carlo predictions for reconstructed distributions are found in Sect. 6, while the unfolding procedure is described in Sect. 7. Section 8 presents the analysis results, the comparisons to predictions, and a discussion of their interpretation. Conclusions are provided in Sect. 9.

2 The ATLAS detector

The ATLAS experiment at the LHC is a multi-purpose particle detector with a forward-backward symmetric cylindrical geometry and nearly $4\pi$ coverage in solid angle.² It consists of an inner tracking detector, electromagnetic and hadronic calorimeters, and a muon spectrometer. The inner tracker is surrounded by a thin superconducting solenoid magnet and provides precision tracking of charged particles and momentum measurements in the pseudorapidity range $|\eta | < 2.5$. This region is matched to a high-granularity electromagnetic (EM) sampling calorimeter covering the pseudorapidity range $|\eta | < 3.2$, and a coarser granularity calorimeter up to $|\eta | = 4.9$. The hadronic calorimeter system covers the entire pseudorapidity range up to $|\eta | = 4.9$. The muon spectrometer consists of three large superconducting toroids each containing eight coils, a system of trigger chambers, and precision tracking chambers, which provide trigger and tracking capabilities in the range $|\eta | < 2.4$ and $|\eta | < 2.7$, respectively. A two-level trigger system [21] is used to select events. The first-level trigger is implemented in hardware and uses a subset of the detector information. This is followed by the software-based high-level trigger system, which runs offline reconstruction, reducing the event rate to approximately 1 kHz.

3 Data set, simulated event samples, and predictions

3.1 Data set

The data used in this analysis were collected by the ATLAS detector during August to November 2015. During this period, the LHC circulated 6.5 $\text {TeV}$proton beams with a 25 ns bunch spacing. The peak delivered instantaneous luminosity was $L={5 \times 10^{33}}$ cm$^{-2}$ s$^{1}$ and the mean number of pp interactions per bunch crossing (hard scattering and pile-up events) was $\langle \mu \rangle = 13$. The data set used in this measurement corresponds to a total integrated luminosity of 3.16 fb$^{-1}$.

3.2 Simulated event samples

Monte Carlo simulations, normalised to higher-order calculations, are used to estimate most of the contributions from background events, to unfold the data to the particle level, and to compare with the unfolded data distributions. All samples are processed with a Geant4-based simulation [22] of the ATLAS detector [23]. An overview of all signal and background processes considered and of the generators used for the simulation is given in Table 1. Total production cross sections for the samples, their respective uncertainties (mainly coming from parton distribution function (PDF) and factorisation and renormalisation scale variations), and references to higher-order QCD corrections, where available, are also listed in Table 1.

Table 1 Signal and background Monte Carlo samples and the generators used in the simulation. Each sample is normalised to the appropriate production cross section $\sigma$ and multiplied by the relevant branching ratios (${\mathrm {BR}}$) per lepton flavour for this sample, as shown in the third column. For $W$-boson and top-quark production, contributions from higher-order QCD corrections were calculated following the references given in the fifth column for the stated order. Similarly, for $Z$-boson production, higher-order QCD corrections were evaluated in the dilepton invariant mass range $66<m_{\ell \ell }<116$ GeV following the references given in the fifth column for the stated order, and extrapolation scaling factors were applied to match mass ranges used by each simulation as given in the first column. The theory uncertainties as given in the final column correspond to PDF and scale variations. The diboson samples include on-shell and off-shell WW, WZ and ZZ production. Recently, NNLO QCD predictions have been made available for the diboson processes [32, 33]. However, these higher-order corrections have a negligible impact on this analysis

Signal events (i.e. containing a $Z$ boson with associated jets) are simulated using the Sherpa  v2.2.1 [31] generator, denoted by Sherpa 2.2. Matrix elements (ME) are calculated for up to two additional partons at NLO and up to four partons at leading order (LO) using the Comix [34] and OpenLoops [35] matrix element generators. They are merged with the Sherpa parton shower [36] (with a matching scale of 20 GeV) using the ME+PS@NLO prescription [37]. A five-flavour scheme is used for these predictions. The NNPDF30NLO PDF set [38] is used in conjunction with a dedicated set of parton-shower-generator parameters (tune) developed by the Sherpa authors. This sample is used for the nominal unfolding of the data distributions, to compare to the cross-section measurements, and to estimate the systematic uncertainties.

A simulated sample of $Z+\text {jets}$ production is also produced with the MADGRAPH_aMC@NLO (denoted by MG5_aMC@NLO) v2.2.2 generator [39], using matrix elements including up to four partons at leading order and employing the NNPDF30NLO PDF set, interfaced to Pythia  v8.186 [40] to model the parton shower, using the CKKWL merging scheme [41] (with a matching scale of 30 GeV). A five-flavour scheme is used. The A14 [42] parton-shower tune is used together with the NNPDF23LO PDF set [43]. The sample is denoted by MG5_aMC+Py8 CKKWL and is used to provide cross-checks of the systematic uncertainty in the unfolding and to model the small $Z \rightarrow \tau \tau$ background. In addition, a MG5_aMC@NLO sample with matrix elements for up to two jets and with parton showers beyond this, employing the NNPDF30NLO PDF set and interfaced to Pythia  v8.186 to model the parton shower, is generated using the FxFx merging scheme [44] (with a matching scale of 25 GeV [45]) and is denoted by MG5_aMC+Py8 FxFx. This sample also uses a five-flavour scheme and the A14 parton-shower tune with the NNPDF23LO PDF set. Both MG5_aMC@NLO samples are used for comparison with the unfolded cross-section measurements.

The measured cross sections are also compared to predictions from the leading-order matrix element generator Alpgen v2.14 [46] interfaced to Pythia  v6.426 [47] to model the parton shower, denoted by Alpgen+Py6, using the Perugia2011C [48] parton-shower tune and the CTEQ6L1 PDF set [49]. A four-flavour scheme is used. Up to five additional partons are modelled by the matrix elements merged with the MLM prescription [46] (with a matching scale of 20 GeV). The matrix elements for the production of $Z+b\bar{b}$ and $Z+c\bar{c}$ events are explicitly included and a heavy-flavour overlap procedure is used to remove the double counting of heavy quarks from gluon splitting in the parton shower.

The $Z$-boson samples are normalised to the NNLO prediction calculated with the Fewz 3.1 program [24,25,26,27] with CT10nnlo PDFs [50].

Contributions from the top-quark, single-boson, and diboson components of the background (described in Sect. 5) are estimated from the following Monte Carlo samples. Samples of top-quark pair and single top-quark production are generated at NLO with the Powheg-Box generator [51,52,53] [versions v2 (r3026) for top-quark pairs and v1 for single top quarks (r2556 and r2819 for t-and Wt-channels, respectively)] and Pythia  v6.428 (Perugia2012 tune [48]). Samples with enhanced or suppressed additional radiation are generated with the Perugia2012radHi/Lo tunes [48]. An alternative top-quark pair sample is produced using the MG5_aMC@NLO generator interfaced with Herwig++  v2.7.1 [39, 54], using the UE-EE-5 tune [55]. The samples are normalised to the cross section calculated at NNLO+NNLL (next-to-next-to-leading log) with the Top++2.0 program [28].

The $W$-boson backgrounds are modelled using the MG5_aMC+Py8 CKKWL  v2.2.2 generator, interfaced to Pythia  v8.186 and are normalised to the NNLO values given in Table 1. Diboson processes with fully leptonic and semileptonic decays are simulated [56] using the Sherpa  v2.1.1 generator with the CT10nlo PDF set. The matrix elements contain the doubly resonant WW, WZ and ZZ processes, and all other diagrams with four electroweak vertices. They are calculated for one or zero additional partons at NLO and up to three additional partons at LO and merged with the Sherpa parton shower using the ME+PS@NLO prescription.

Events involving semileptonic decays of heavy quarks, hadrons misidentified as leptons, and, in the case of the electron channel, electrons from photon conversions are referred to collectively as “multijet events”. The multijet background was estimated using data-driven techniques, as described in Sect. 5.

Multiple overlaid pp collisions are simulated with the soft QCD processes of Pythia  v.8.186 using the A2 tune [57] and the MSTW2008LO PDF set [58]. All Monte Carlo samples are reweighted so that the pile-up distribution matches that observed in the data.

3.3 Fixed-order predictions

In addition to these Monte Carlo samples, parton-level fixed-order predictions at NLO are calculated by the BlackHat+Sherpa collaboration for the production of $Z$ bosons with up to four partons [59, 60]. The BlackHat+Sherpa predictions use the CT14 PDF set [61] including variations of its eigenvectors at the 68% confidence level, rescaled from 90% confidence level using a factor of 1 / 1.645. The nominal predictions use a factorisation and renormalisation scale of $H_{\text {T}}/2$ with uncertainties derived from the envelope of a common variation of the scales by factors of $0.5, 1/\sqrt{2}, \sqrt{2},$ and 2. The effects of PDF and scale uncertainties range from 1 to 4% and from 0.1 to 10%, respectively, for the cross sections of $Z$-boson production in association with at least one to four partons, and are included in the predictions which are provided by the BlackHat+Sherpa authors for the fiducial phase space of this analysis. Since these predictions are defined before the decay leptons emit photons via final-state radiation (Born-level FSR), corrections to the dressed level (where all photons found within a cone of size $\Delta R = 0.1$ of the lepton from the decay of the $Z$ boson are included) are derived from MG5_aMC+Py8 CKKWL, separately for each kinematic observable used to measure cross sections, with associated systematic uncertainties obtained by comparing to the Alpgen+Py6 generator. This correction is needed in order to match the prediction to the lepton definition used in the measurements. The average size of these corrections is approximately $-2$%. To bring the prediction from parton to particle level, corrections for the non-perturbative effects of hadronisation and the underlying event are also calculated separately for each observable using the Sherpa  v2.2 generator by turning on and off in the simulation both the fragmentation and the interactions between the proton remnants. The net size of the corrections is up to approximately 10% at small values of $p_{\text {T}}^{\text {jet}}$ and vanishes for large values of $p_{\text {T}}^{\text {jet}}$. An uncertainty of approximately 2% for this correction is included in the total systematic uncertainty of the prediction.

Calculations of cross sections at NNLO QCD have recently become available [3,4,5,6]. In this paper, the results are compared to the calculation, denoted by $Z+ \ge 1\text {~jet~N}_{\text {jetti}}$ NNLO [3, 4], which uses a new subtraction technique based on N-jettiness [62] and relies on the theoretical formalism provided in soft-collinear effective theory. The predictions, which are provided by the authors of this calculation for the fiducial phase space of this analysis, use a factorisation and renormalisation scale of $\sqrt{m_{\ell \ell }^2 + \sum \left( p_{\text {T}}^{\text {jet}}\right) ^2}$ (where $m_{\ell \ell }$ is the invariant mass of the dilepton system) and the CT14 PDF set. The QCD renormalisation and factorisation scales were jointly varied by a common factor of two, and are included in the uncertainties. Non-perturbative and FSR corrections and their associated uncertainties as discussed above are also included in the predictions.

4 Event selection

Electron- and muon-candidate events are selected using triggers which require at least one electron or muon with transverse momentum thresholds of $p_{\text {T}} =24\,\text {GeV}$ or 20 $\text {GeV}$, respectively, with some isolation requirements for the muon trigger. To recover possible efficiency losses at high momenta, additional electron and muon triggers which do not make any isolation requirements are included with thresholds of $p_{\text {T}} \ge 60\,\text {GeV}$ and $p_{\text {T}} = 50\,\text {GeV}$, respectively. Candidate events are required to have a primary vertex, defined as the vertex with the highest sum of track $p_{\text {T}} ^2$, with at least two associated tracks with $p_{\text {T}} > 400\,\text {MeV}$.

Electron candidates are required to have $p_{\text {T}} > 25\,\text {GeV}$ and to pass “medium” likelihood-based identification requirements [63, 64] optimised for the 2015 operating conditions, within the fiducial region of $|\eta |<2.47$, excluding candidates in the transition region between the barrel and endcap electromagnetic calorimeters, $1.37<|\eta |<1.52$. Muons are reconstructed for $|\eta |<2.4$ with $p_{\text {T}} > 25\,\text {GeV}$ and must pass “medium” identification requirements [65] also optimised for the 2015 operating conditions. At least one of the lepton candidates is required to match the lepton that triggered the event. The electrons and muons must also satisfy $p_{\text {T}}$-dependent cone-based isolation requirements, using both tracking detector and calorimeter information (described in Refs. [66, 67], respectively). The isolation requirements are tuned so that the lepton isolation efficiency is at least 90% for $p_{\text {T}} >25\,\text {GeV}$, increasing to 99% at 60 GeV. Both the electron and muon tracks are required to be associated with the primary vertex, using constraints on the transverse impact parameter significance $|d_0|/\Delta d_0$, where $d_0$ is the transverse impact parameter and $\Delta d_0$ is its uncertainty, and on the longitudinal impact parameter $z_0$ corrected for the reconstructed position of the primary vertex. The transverse impact parameter significance is required to be less than five for electrons and three for muons, while the absolute value of the corrected $z_0$ multiplied by the sine of the track polar angle is required to be less than 0.5 mm.

Jets of hadrons are reconstructed with the anti-$k_t$ algorithm [68] with radius parameter $R=0.4$ using topological clusters of energy deposited in the calorimeters [69]. Jets are calibrated using a simulation-based calibration scheme, followed by in situ corrections to account for differences between simulation and data [70]. In order to reduce the effects of pile-up contributions, jets with pseudorapidity $|\eta |<2.4$ and $p_{\text {T}} <{60} {\text {GeV}}$ are required to have a significant fraction of their tracks with an origin compatible with the primary vertex, as defined by the jet vertex tagger algorithm [71]. In addition, the expected average energy contribution from pile-up clusters is subtracted according to the $\eta$–$\phi$ catchment area of the jet [72]. Jets used in the analysis are required to have $p_{\text {T}}$ greater than 30 GeV and rapidity $|y|<2.5$.

The overlap between leptons and jets is removed in a two-step process. The first step removes jets closer than $\Delta R = 0.2$ to a selected electron, and jets closer than $\Delta R = 0.2$ to a selected muon, if they are likely to be reconstructed from photons radiated by the muon. In a second step, electrons and muons are discarded if they are located closer than $\Delta R = 0.4$ to a remaining selected jet. This requirement effectively removes events with leptons and jets which are not reliably simulated in the Monte Carlo simulation.

Events containing a $Z$-boson candidate are selected by requiring exactly two leptons of the same flavour but of opposite charge with dilepton invariant mass in the range $71< m_{\ell \ell } < 111$ GeV. The expected and observed numbers of $Z$-boson candidates selected for each inclusive jet multiplicity, for $N_{\mathrm {jets}} \ge 0-7$, are summarised in Table 2, separately for the $Z \rightarrow e^{+}e^{-}$ and the $Z \rightarrow \mu ^{+}\mu ^{-}$ channels. The background evaluation that appears in this table is discussed in Sect. 5. After all requirements, 248,816 and 311,183 $Z+ \ge 1\text {~jet}$ events are selected in the electron and muon channels, respectively.

Table 2 Fraction of signal and background processes in % in the final selection and expected and observed numbers of events for the various inclusive jet multiplicities considered in the electron (top) and muon (bottom) channels

4.1 Correction factors and related systematic uncertainties

Some of the object and event selection efficiencies as well as the energy and momentum calibrations modelled by the simulation must be corrected with simulation-to-data correction factors to better match those observed in the data. These corrections and their corresponding uncertainties fall into the following two categories: dependent and not dependent on lepton flavour.

The corrections and uncertainties specific to each leptonic final state ($Z \rightarrow e^{+}e^{-}$ and $Z \rightarrow \mu ^{+}\mu ^{-}$) are as follows:

• Trigger: The lepton trigger efficiency is estimated in simulation, with a separate data-driven analysis performed to obtain the simulation-to-data trigger correction factors and their corresponding uncertainties [21].

• Lepton reconstruction, identification, and isolation: The lepton selection efficiencies as determined from simulation are also corrected with simulation-to-data correction factors, with corresponding uncertainties [64, 65].

• Energy, momentum scale/resolution: Uncertainties in the lepton calibrations are estimated [65] because they can cause a change of acceptance because of migration of events across the $p_{\text {T}}$ threshold and $m_{\ell \ell }$ boundaries.

The corrections and uncertainties common to the electron and muon final states are as follows:

• Jet energy scale and resolution: Uncertainties in the jet energy-scale calibration and resolution have a significant impact on the measurements, especially for the higher jet multiplicities. The jet energy-scale calibration is based on 13 TeV simulation and on in situ corrections obtained from data [70]. The uncertainties are estimated using a decorrelation scheme, resulting in a set of 19 independent parameters which cover all of the relevant calibration uncertainties. The jet energy scale is the dominant systematic uncertainty for all bins with at least one jet. The jet energy-resolution uncertainty is derived by over-smearing the jet energy in the simulation and using the symmetrised variations as the uncertainty.

• Jet vertex tagger: The modelling of the output variable from the jet vertex tagger is validated using data events where the $Z$ boson recoils against a jet. A percent-level correction is derived and its statistical and systematic uncertainties are used as additional uncertainties in the efficiency to select jets from the primary vertex [71].

• Pile-up: The imperfect modelling of the effects of pile-up leads to acceptance changes at the percent level for different jet multiplicities. To assess this uncertainty, the average number of interactions per bunch crossing $\langle \mu \rangle$ is varied in simulation so that the behaviour of variables sensitive to pile-up matches that observed in data.

• Luminosity: The cross sections have a 2.1% uncertainty from the measurement of the integrated luminosity, which is derived, following a methodology similar to that detailed in Refs. [73, 74], from a calibration of the luminosity using x–y beam-separation scans performed in August 2015.

5 Background estimation

Contributions from the electroweak (single boson and diboson) and top-quark (single top-quark and top-quark pair) components of the background are estimated using the Monte Carlo samples described in Sect. 3 with corresponding uncertainties as listed in Table 1. Contributions from multijet events are evaluated with data-driven techniques as described below. A summary of the composition and relative importance of the backgrounds in the candidate $Z+\text {jets}$ events is given in Table 2. The overall purity of the $Z+\text {jets}$ selections (fraction of signal events in the final selection) ranges from 99% in the inclusive sample to 80–85% in the $\ge \mathrm {7~jets}$ bin.

5.1 Top-quark and electroweak backgrounds

The dominant contribution to the background at high jet multiplicities comes from $t\overline{t}$ production, with the subsequent leptonic decays of the $W$ bosons originating from the top quarks and is evaluated from simulation. An overall uncertainty of 6%, corresponding to the PDF and scale variations on the theoretical predictions of the inclusive cross sections, is assigned (see Table 1). The $t\overline{t}$ background estimate is validated through a cross-section measurement of $t\overline{t}$ production in the dilepton channel at $\sqrt{s} = 13$ TeV [75] as a function of the jet multiplicity, and the modelling of the additional parton radiation in $t\overline{t}$ events by Powheg+Py6 was found to be in good agreement with this measurement. In addition, a systematic uncertainty in the modelling of the shape of the distributions is derived by modifying the parton-shower intensity in the nominal simulation sample and by comparing to the predictions from the alternative generator MG5_aMC@NLO+Herwig++ (both listed in Table 1). The small contribution from single-top-quark events is also estimated using Powheg+Py6 samples and assigned a 6% uncertainty.

Diboson production in leptonic and semileptonic final states with at least two leptons of the same flavour constitutes a co-dominant background for high jet multiplicities (see Table 2). The production of WZ bosons in association with jets at $\sqrt{s}=13$ TeV was found to be well modelled by the Sherpa 2.1 generator [76]. A 6% uncertainty, again corresponding to PDF and scale variations on the predictions, is assessed. Since in Ref. [76] the measurement is limited by the statistical precision for dibosons $+\ge 4$ jets (resulting in $\ge 6$ hadronic jets for semileptonic diboson decays), an additional systematic uncertainty of 50% in the normalisation of the diboson background is added for $Z+ \ge 6\text {~jets}$.

Minor background contributions also arise from single-$W$-boson production decaying to leptonic final states and from single-$Z$-boson production in the $Z \rightarrow \tau ^{+}\tau ^{-}$ process, both estimated with simulation and assigned a 5% uncertainty (as given in Table 1).

5.2 Multijet background

Background-enriched data control regions are used to estimate the multijet contribution in both the electron and muon channels. They are constructed by loosening the lepton identification and isolation requirements. Templates are built from the dilepton invariant mass distribution, a variable that shows discrimination between multijet background and other processes in regions of its kinematic range, but is largely uncorrelated with the variables used to build the multijet control regions. The templates are subsequently normalised to events passing the $Z$-boson signal selection.

In the electron channel, the multijet templates are built for each jet multiplicity from events with two same-charge leptons with no isolation requirement, whose identification criteria are looser than those of the signal selection, which the leptons must not satisfy. In the muon channel, the control region is similarly built from events with two leptons which are selected with looser identification requirements than the signal selection and also fail the nominal isolation requirement. In both cases, dedicated triggers better suited to this purpose are used to populate the templates. The small electroweak and top-quark contamination is subtracted using simulated events.

The normalisation of the multijet template is estimated with a log-likelihood fit to the measured dilepton invariant mass distribution for the inclusive $Z$ selection, using templates for $Z \rightarrow \ell ^{+}\ell ^{-}$ and for the electroweak and top-quark background derived from simulation. The fit is performed in the invariant mass windows of $52< m_{ee} < 148$ GeV and $40< m_{\mu \mu } < 80$ GeV for the electron and muon channels, respectively, in order to benefit from the larger multijet contribution in the mass sidebands. The normalisation of the multijet template is allowed to float freely while the remaining non-multijet templates are constrained to be within 6% of the predicted cross sections for these processes as given in Table 1. The multijet fractions are evaluated separately for each jet multiplicity, except for very high jet multiplicities where the templates are statistically limited, and so these fractions are taken from the estimates of the $\ge \mathrm {5~jets}$ and $\ge \mathrm {4~jets}$ bins in the electron and muon channels, respectively.

The systematic uncertainties on the multijet background are derived by varying the mass range and bin width of the nominal fit, using the lepton transverse impact parameter $d_{0}$ as the fitting variable instead of the invariant mass, using alternative simulation samples for the templates, allowing the normalisations of the non-multijet components to vary independently or within a wider range, and varying the lepton resolution and energy/momentum scales. In addition, given the multiple sources of multijet background in the electron channel, an alternative template is constructed by requiring that the electrons fail to meet an isolation criterion instead of failing to meet the nominal signal selection electron identification criterion.

The resulting estimated multijet fractions in each jet multiplicity bin are given in Table 2. Their corresponding total uncertainties are dominated by their systematic components. These systematic components are approximately 70% of the multijet fraction as estimated in the electron and muon channels.

6 Kinematic distributions

The level of agreement between data and predictions is evaluated from the comparison of kinematic distributions. Figure 1, which presents the dilepton mass for the $Z+ \ge 1\text {~jet}$ topology and the inclusive jet multiplicity, shows how well the Sherpa 2.2 and MG5_aMC+Py8 CKKWL predictions agree with data. The uncertainty bands shown in these distributions include the statistical uncertainties due to the simulation sample sizes, the event-selection uncertainties described in Sect. 4.1 (omitting the common 2.1% luminosity uncertainty), and the background normalisation uncertainties described in Sect. 5.

Fig. 1

Dilepton invariant mass for $Z+ \ge 1\text {~jet}$ (top) and inclusive jet multiplicity (bottom) in the $Z (\rightarrow e^{+}e^{-})+\text {~jets}$ (left) and the $Z (\rightarrow \mu ^{+}\mu ^{-})+\text {~jets}$ (right) channels. All backgrounds and the signal samples are stacked to produce the figures. Systematic uncertainties for the signal and background distributions are combined in the hatched band, and the statistical uncertainty is shown on the data points. The uncertainty in the luminosity and the theory uncertainty in the signal prediction are not included in the uncertainty band

7 Unfolding of detector effects

The cross-section measurements presented in this paper are performed within the fiducial acceptance region defined by the following requirements:

• $p_{\text {T}} ^{\ell }>25$ GeV, $|\eta ^{\ell }|<2.5$

• $p_{\text {T}}^{\text {jet}} > 30$ GeV, $|y^{\text {jet}} |<2.5$

• $\Delta R(\ell ,\text {jet}) > 0.4$

• $71<m_{\ell \ell }<111$ GeV.

The cross sections are defined at particle (“truth”) level, corresponding to dressed electrons and muons from the $Z$ bosons. The particle level also includes jets clustered using the anti-$k_t$ algorithm [68] with radius parameter $R=0.4$ for final-state particles with decay length $c\tau >10$ mm, excluding the dressed $Z$-boson decay products.

The fiducial cross sections are estimated from the reconstructed kinematic observables: jet multiplicity, $p_{\text {T}}^{\text {jet}}$ for different jet multiplicities, $y^{\text {jet}}$, $\Delta \phi _{\text {jj}}$, $m_{\text {jj}}$, and $H_{\text {T}}$, for events that pass the selection described in Sect. 4. The expected background components as described in Sect. 5 are subtracted from the distributions in data. A variable-width binning of these observables is used, such that the purity is at least 50% in each bin and the size of the statistical uncertainty in most of the bins remains below 10%.

Fig. 2

Response matrices corresponding to the exclusive jet multiplicity for $Z+\text {jets}$ events in the electron channel (left) and to the $H_{\text {T}}$ for $Z+ \ge 1\text {~jet}$ events in the muon channel (right). The sum of the entries in each row is normalised to unity. Both matrices are obtained from Sherpa 2.2

An iterative Bayesian unfolding technique [77], as implemented in the RooUnfold package [78], is used to unfold the measurements to the particle level, thereby accounting for detector effects related to inefficiencies, resolution, and systematic biases in the central values of the kinematic variables describing both the leptons and the jets. The iterative unfolding technique updates the initial estimators for the generated (“truth”) distribution in consecutive steps, using Bayes’ theorem in each iteration to derive an unfolding matrix from the initial response matrix (which relates truth and reconstructed distributions of given observables) and the current truth estimator.

The response matrices are constructed using the Sherpa 2.2 $Z (\rightarrow \ell ^{+}\ell ^{-})+\text {~jets}$ samples. Sherpa 2.2 is also used to derive the initial truth estimator. In order to enter the response matrix, events must pass the $Z$-boson selection at generator level and at detector level and contain the number of jets required by the preselection for a given observable at both generator and detector level. Reconstructed jets are required to match the corresponding generator-level jets within a cone of size $\Delta R=0.4$ for all distributions except global quantities such as the jet multiplicity and $H_{\text {T}}$. A given bin (i, j) in the response matrix therefore corresponds to the probability that a true jet object in bin j is reconstructed in bin i of the distribution. Figure 2 illustrates two examples of response matrices. The resulting ratios of detector-level to truth-level event yields are typically 0.65 and 0.8 for the electron and muon channels, respectively.

The background-subtracted data are corrected for the expected fraction of events with reconstructed objects unmatched to any generator object before entering the iterative unfolding. The number of iterations used for the iterative unfolding of each distribution (two) is chosen by unfolding the Sherpa 2.2 samples reweighted to data and comparing to the generated reweighted distribution. The unfolded event yields are divided by the integrated luminosity of the data sample and the bin width of the distribution in question to provide the final fiducial cross sections. The final result is given by

$$\begin{aligned} \sigma _i = \frac{1}{\epsilon _i L} \sum _j U_{ij} N_j^{\text {data}} \left( 1-f_j^{\text {unmatched}}\right) , \end{aligned}$$

(1)

where L is the integrated luminosity, $\epsilon _i$ is the reconstruction efficiency for truth bin i, $N_j^{\text {data}}$ corresponds to the number of events observed in data in reconstructed bin j and $f_j^{\text {unmatched}}$ is its fraction of unmatched events calculated from simulation, and $U_{ij}$ is the unfolding matrix calculated after two iterations, using the updated prior from the first iteration and the response matrix.

7.1 Systematic uncertainties associated with the unfolding procedure

The limited size of a simulation sample can create biases in the distributions. Systematic uncertainties account for possible residual biases in the unfolding procedure due to, e.g. modelling of the hadronisation in the simulation, migrations into other kinematic distributions not explicitly part of the unfolding, or the finite bin width used in each distribution. The following uncertainties arise from the unfolding procedure.

• The statistical uncertainties of the response matrices derived from Sherpa 2.2 are propagated to the unfolded cross sections with a toy simulation method. A total of 5000 ensembles (pseudo-experiments) of unfolded samples are generated. For each sample, the number of reconstructed events in each bin is generated randomly according to a Gaussian distribution, where the mean is the nominal number of events before unfolding and the width is its corresponding statistical uncertainty. Unfolding is performed for each ensemble. The widths of resulting distributions are taken as a systematic uncertainty of the unfolding.

• The Sherpa 2.2 samples are reweighted at generator level, such that the distribution of the leading jet $p_{\text {T}}$ at detector level matches that observed in the data. The modified Sherpa 2.2 samples are then used to unfold the data again and the variations in the resulting cross sections are used to derive a systematic uncertainty.

• An additional check is performed by unfolding reconstructed MG5_aMC+Py8 CKKWL events using Sherpa 2.2 response matrices. The residual non-closure is accounted for by an additional flat uncertainty of 3% for all distributions.

8 Results

The measured cross sections, presented in Sect. 8.1, are calculated in the electron and muon channels separately and the compatibility of the results of the two channels is evaluated. In order to improve the precision of the measurement, these results are then combined, taking into account the correlations of the systematic uncertainties. The comparisons of the combined results to the predictions are presented in Sect. 8.2.

8.1 Results in the individual channels and the combination

The fiducial cross-section measurements in the $Z (\rightarrow e^{+}e^{-})+\text {~jets}$ and $Z (\rightarrow \mu ^{+}\mu ^{-})+\text {~jets}$ channels as a function of the inclusive jet multiplicities are presented in Table 3. The data statistical uncertainties are propagated through the unfolding by using pseudo-experiments. As mentioned in Sect. 7, the systematic uncertainties are propagated through the unfolding via the migration matrices and via the variation of the subtracted background. Table 4 shows the resulting total relative statistical and systematic uncertainties as well as the systematic components [lepton trigger, lepton selection, jet energy scale and resolution, jet vertex tagging, pile-up, luminosity (all described in Sect. 4.1)], unfolding (described in Sect. 7), and background (described in Sect. 5) as a function of the inclusive jet multiplicity, presented separately for the electron and muon channels. The jet energy scale is the dominant systematic uncertainty for all bins with at least one jet.

Table 3 Measured fiducial cross sections in the electron and muon channels for successive inclusive jet multiplicities. The total statistical and systematic uncertainties are given, along with the uncertainty in the luminosity

Table 4 Relative statistical and systematic uncertainties (in %) in the measured cross sections of $Z+\text {jets}$ production for successive inclusive jet multiplicities in the electron (top) and muon (bottom) channels

Figure 3 shows a comparison of the electron and muon channels for the measured fiducial cross section as a function of the inclusive jet multiplicity and of the leading jet $p_{\text {T}}$ for inclusive $Z+ \ge 1\text {~jet}$ events. This figure demonstrates that the results in the electron and muon channels are compatible and hence can be combined to improve the precision of the measurement. This figure also shows the result of this combination described below.

Fig. 3

Measured fiducial cross section as a function of the inclusive jet multiplicity (left) and the leading jet $p_{\text {T}}$ for inclusive $Z+ \ge 1\text {~jet}$ events (right) in the electron and the muon channels and compared to their combined value. The ratios of the two measurements to the combined results are also shown in the bottom panels. The error bars indicate the statistical uncertainty, and the hatched bands the statistical and the flavour-uncorrelated systematic uncertainties of the combined result, added in quadrature

The results from the electron and muon channels are combined at dressed level for each distribution separately: inclusive and exclusive jet multiplicities, ratio for successive inclusive jet multiplicities, leading jet $p_{\text {T}}$ for $Z+ \ge 1,2,3,4\text {~jet}$ events and jet $p_{\text {T}}$ for exclusive $Z+1\text {~jet}$ events, leading jet rapidity for inclusive $Z+ \ge 1\text {~jet}$ events, $H_{\text {T}}$, $\Delta \phi _{\text {jj}}$, and $m_{\text {jj}}$. A $\chi ^2$ function whose sum runs over all measurement sets (electrons and muons), measurement points, and some of the uncertainty sources, is used for the combination [79, 80] and distinguishes between bin-to-bin correlated and uncorrelated sources of uncertainties, the latter comprising the statistical uncertainty of the data and the statistical unfolding uncertainty. Uncertainties specific to the lepton flavour and to the background are included in the $\chi ^2$ function, while the remaining, flavour-uncorrelated, systematic uncertainties related to jets, pile-up, luminosity, and unfolding are averaged after the combination.

Fig. 4

Measured cross section as a function of the inclusive jet multiplicity (left) and ratio for successive inclusive jet multiplicities (right) for inclusive $Z+\text {jets}$ events. The data are compared to the predictions from BlackHat+Sherpa, Sherpa 2.2, Alpgen+Py6, MG5_aMC+Py8 CKKWL, and MG5_aMC+Py8 FxFx. The error bars correspond to the statistical uncertainty, and the hatched bands to the data statistical and systematic uncertainties (including luminosity) added in quadrature. A constant 5% theoretical uncertainty is used for Sherpa 2.2, Alpgen+Py6, MG5_aMC+Py8 CKKWL, and MG5_aMC+Py8 FxFx. Uncertainties from the parton distribution functions and QCD scale variations are included in the BlackHat+Sherpa predictions, as described in Sect. 3.3

8.2 Comparisons of results to predictions

The cross-section measurement for different inclusive $Z+\text {jets}$ multiplicities and their ratios obtained from the combination are found in Tables 5 and 6. Figure 4 shows the comparison of these results with the NLO QCD fixed-order calculations from BlackHat+Sherpa and with the predictions from Sherpa 2.2, Alpgen+Py6, MG5_aMC+Py8 CKKWL, and MG5_aMC+Py8 FxFx. The plots show the particle-level cross section with the generator predictions normalised to the inclusive NNLO cross sections in the top panel, accompanied by the ratios of the various predictions with respect to the data in the bottom panels. Uncertainties from the parton distributions functions and QCD scale variations are included in the BlackHat+Sherpa predictions, as described in Sect. 3.3. A constant 5% theoretical uncertainty is used for Sherpa 2.2, Alpgen+Py6, MG5_aMC+Py8 CKKWL, and MG5_aMC+Py8 FxFx, as described in Table 1. The inclusive jet multiplicity decreases logarithmically while the ratio is flat in the presence of at least one jet. The predictions are in agreement with the observed cross sections and their ratios, except for Sherpa 2.2, Alpgen+Py6 and MG5_aMC+Py8 FxFx for high jet multiplicity, where a non-negligible fraction of the jets are produced by the parton shower.

Table 5 Measured combined fiducial cross sections for successive inclusive jet multiplicities. The statistical, systematic, and luminosity uncertainties are given

Table 6 Measured combined ratios of the fiducial cross sections for successive inclusive jet multiplicities. The statistical, systematic, and luminosity uncertainties are given

Fig. 5

Measured cross section as a function of the leading jet $p_{\text {T}}$ for inclusive $Z+ \ge 1,2,3,4\text {~jet}$ events (left) and absolute value of the leading jet rapidity for inclusive $Z+ \ge 1\text {~jet}$ events (right). The data are compared to the predictions from $Z+ \ge 1\text {~jet~N}_{\text {jetti}}$ NNLO, BlackHat+Sherpa, Sherpa 2.2, Alpgen+Py6, MG5_aMC+Py8 CKKWL, and MG5_aMC+Py8 FxFx. The error bars correspond to the statistical uncertainty, and the hatched bands to the data statistical and systematic uncertainties (including luminosity) added in quadrature. The details of the prediction uncertainties are given in the text. For clarity, uncertainty bands are not shown for the Monte Carlo predictions in the left-hand plot. Uncertainties from the QCD scale variations for the $Z+ \ge 1\text {~jet~N}_{\text {jetti}}$ NNLO predictions are included, as described in Sect. 3.3

Fig. 6

Measured cross section as a function of jet $p_{\text {T}}$ for exclusive $Z+1\text {~jet}$ events (left) and exclusive jet multiplicity (right). The data are compared to the predictions from BlackHat+Sherpa, Sherpa 2.2, Alpgen+Py6, MG5_aMC+Py8 CKKWL, and MG5_aMC+Py8 FxFx. The error bars correspond to the statistical uncertainty, and the hatched bands to the data statistical and systematic uncertainties (including luminosity) added in quadrature. The details of the prediction uncertainties are given in the text

The jet transverse momentum is a fundamental observable of the $Z+\text {jets}$ process and probes pQCD over a wide range of scales. Moreover, understanding the kinematics of jets in events with vector bosons associated with several jets is essential for the modelling of backgrounds for other SM processes and searches beyond the SM. The leading jet $p_{\text {T}}$ distribution (which is correlated with the $p_{\text {T}}$ of the $Z$ boson) in inclusive $Z+ \ge 1,2,3,4\text {~jet}$ events is shown in Fig. 5 and ranges up to 700 GeV. The LO generator MG5_aMC+Py8 CKKWL models a too-hard jet $p_{\text {T}}$ spectrum. This feature is known from studies of LO generators in pp collisions at lower centre-of-mass energies [11], and can be interpreted as an indication that the dynamic factorisation and renormalisation scale used in the generation is not appropriate for the full jet $p_{\text {T}}$ range. In contrast, the predictions from BlackHat+Sherpa, Sherpa 2.2, and MG5_aMC+Py8 FxFx, which are based on NLO matrix elements, are in agreement with the measured cross section within the systematic uncertainties over the full leading jet $p_{\text {T}}$ range. Alpgen+Py6 also shows good agreement with the measured data. The $Z+ \ge 1\text {~jet~N}_{\text {jetti}}$ NNLO prediction models the spectrum for the $Z+ \ge 1\text {~jet}$ events well. Uncertainties from the QCD scale variations for the $Z+ \ge 1\text {~jet~N}_{\text {jetti}}$ NNLO predictions are included in the uncertainty band, as described in Sect. 3.3. For the leading jet rapidity distribution in inclusive $Z+ \ge 1\text {~jet}$ events, also shown in this figure, all predictions show good agreement with the measured data within the uncertainties.

The exclusive jet $p_{\text {T}}$ distribution probes the validity of $Z+1\text {~jet}$ predictions at increasing QCD scales represented by the jet $p_{\text {T}}$ in the presence of a jet veto at a constant low scale; for a jet $p_{\text {T}}$ range of several hundred GeV, accessible with the current data set, the jet scale is of order ten times larger than the veto scale (30 GeV). Figure 6 demonstrates that all predictions studied are consistent with the data within systematic uncertainties over the full jet $p_{\text {T}}$ range (up to 500 GeV). This figure also shows the measured cross section as a function of the exclusive jet multiplicity, which decreases logarithmically. Similar trends as for the inclusive jet multiplicity (Fig. 4) are observed.

Fig. 7

Measured cross section as a function of $H_{\text {T}}$ for inclusive $Z+ \ge 1\text {~jet}$ events. The data are compared to the predictions from $Z+ \ge 1\text {~jet~N}_{\text {jetti}}$ NNLO, BlackHat+Sherpa, Sherpa 2.2, Alpgen+Py6, MG5_aMC+Py8 CKKWL, and MG5_aMC+Py8 FxFx. The error bars correspond to the statistical uncertainty, and the hatched bands to the data statistical and systematic uncertainties (including luminosity) added in quadrature. The details of the prediction uncertainties are given in the text

Quantities based on inclusive $p_{\text {T}}$ sums of final-state objects, such as $H_{\text {T}}$, the scalar $p_{\text {T}}$ sum of all visible objects in the final state, are often employed in searches for physics beyond the Standard Model, to enrich final states resulting from the decay of heavy particles. The values $H_{\text {T}}$ or $H_{\text {T}}$/2 are also commonly used choices for scales for higher-order perturbative QCD calculations. Large values for this quantity can result either from a small number of very energetic particles or from a large number of less energetic particles. Figure 7 shows the measured cross sections as a function of the $H_{\text {T}}$ distribution (up to 1400 GeV) in inclusive $Z+ \ge 1\text {~jet}$ events. The predictions from Sherpa 2.2, Alpgen+Py6 and MG5_aMC+Py8 FxFx describe well the $H_{\text {T}}$ distribution. The prediction from MG5_aMC+Py8 CKKWL describes well the turn-over in the softer part of the $H_{\text {T}}$ spectrum, but overestimates the contribution at large values of $H_{\text {T}}$, in line with the overestimate of the cross sections for hard jets. The fixed-order $Z+ \ge 1\text {~jet}$ prediction from BlackHat+Sherpa underestimates the cross section for values of $H_{\text {T}} >300$ GeV, as observed in similar measurements at lower centre-of-mass energies [11, 81], due to the missing contributions from events with higher parton multiplicities, which for large values of $H_{\text {T}}$ constitute a substantial portion of the data. Agreement is recovered by adding higher orders in perturbative QCD, as demonstrated by the good description of $H_{\text {T}}$ by $Z+ \ge 1\text {~jet~N}_{\text {jetti}}$ NNLO.

Fig. 8

Measured cross section as a function of $\Delta \phi _{\text {jj}}$ (left) and $m_{\text {jj}}$ (right) for inclusive $Z+ \ge 2\text {~jet}$ events. The data are compared to the predictions from BlackHat+Sherpa, Sherpa 2.2, Alpgen+Py6, MG5_aMC+Py8 CKKWL, and MG5_aMC+Py8 FxFx. The error bars correspond to the statistical uncertainty, and the hatched bands to the data statistical and systematic uncertainties (including luminosity) added in quadrature. The details of the prediction uncertainties are given in the text

Angular relations between the two leading jets and the dijet mass are frequently used to separate either heavier SM particles or beyond-SM physics from the $Z+\text {jets}$ process. Figure 8 shows the differential cross section as a function of azimuthal angular difference between the two leading jets for $Z+ \ge 2\text {~jet}$ events, $\Delta \phi _{\text {jj}}$. The tendency of the two jets to be back-to-back in the transverse plane is well modelled by all predictions. This figure also shows the measured cross sections as a function of the invariant mass $m_{\text {jj}}$ of the two leading jets for $Z+ \ge 2\text {~jet}$ events. The shape of the dijet mass is modelled well by BlackHat+Sherpa, Sherpa 2.2, Alpgen+Py6, and MG5_aMC+Py8 FxFx, whereas MG5_aMC+Py8 CKKWL shows a harder spectrum.

9 Conclusion

Proton–proton collision data at $\sqrt{s}=13$ TeV from the LHC, corresponding to a total integrated luminosity of 3.16 fb$^{-1}$, have been analysed by the ATLAS collaboration to study events with $Z$ bosons decaying to electron or muon pairs, produced in association with one or more jets. The fiducial production cross sections for $Z+ \ge 0$–7 jets have been measured, within the acceptance region defined by $p_{\text {T}} ^{\ell }>25$ GeV, $|\eta ^{\ell }|<2.5$, $71<m_{\ell \ell }<111$ GeV, $p_{\text {T}}^{\text {jet}} > 30$ GeV, $|y^{\text {jet}} |<2.5$, and $\Delta R(\ell ,\text {jet}) > 0.4$, with a precision ranging from 4 to 30%. Ratios of cross sections for successive jet multiplicities and cross-section measurements as a function of different key variables such as the jet multiplicities, jet $p_{\text {T}}$ for exclusive $Z+1$ jet events, leading jet $p_{\text {T}}$ for $Z+ \ge 1,2,3,4\text {~jet}$ events, leading jet rapidity for $Z+ \ge 1\text {~jet}$ events, $H_{\text {T}}$, $\Delta \phi _{\text {jj}}$ and $m_{\text {jj}}$ have also been derived.

The measurements have been compared to fixed-order calculations at NLO from BlackHat+Sherpa and at NNLO from the $Z+ \ge 1\text {~jet~N}_{\text {jetti}}$ NNLO calculation, and to predictions from the generators Sherpa 2.2, Alpgen+Py6, MG5_aMC+Py8 CKKWL, and MG5_aMC+Py8 FxFx. In general, the predictions are in good agreement with the observed cross sections and cross-section ratios within the uncertainties. Distributions which are dominated by a single jet multiplicity are modelled well by fixed-order NLO calculations, even in the presence of a jet veto at a low scale. The ME+PS generator MG5_aMC+Py8 CKKWL, which is based on LO matrix elements, models a too-hard jet spectrum, as observed in $\sqrt{s}=7$ TeV pp collisions. It however models well the inclusive jet multiplicity distribution over the full multiplicity range. The modelling of the jet $p_{\text {T}}$ and related observables is significantly improved by the ME+PS@NLO generators Sherpa 2.2 and MG5_aMC+Py8 FxFx, which use NLO matrix elements for up to two additional partons. The recent $Z+ \ge 1\text {~jet~N}_{\text {jetti}}$ NNLO predictions describe well key distributions such as the leading jet $p_{\text {T}}$ and $H_{\text {T}}$. The results presented in this paper provide essential input for the further optimisation of the Monte Carlo generators of Z+jets production and constitute a powerful test of perturbative QCD for processes with a higher number of partons in the final state.