Measurements of the production cross-section for a Z boson in association with b-jets in proton-proton collisions at s\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \sqrt{s} $$\end{document} = 13 TeV with the ATLAS detector

This paper presents a measurement of the production cross-section of a Z boson in association with b-jets, in proton-proton collisions at s\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \sqrt{s} $$\end{document} = 13 TeV with the ATLAS experiment at the Large Hadron Collider using data corresponding to an integrated luminosity of 35.6 fb−1. Inclusive and differential cross-sections are measured for events containing a Z boson decaying into electrons or muons and produced in association with at least one or at least two b-jets with transverse momentum pT> 20 GeV and rapidity |y| < 2.5. Predictions from several Monte Carlo generators based on leading-order (LO) or next-to-leading-order (NLO) matrix elements interfaced with a parton-shower simulation and testing different flavour schemes for the choice of initial-state partons are compared with measured cross-sections. The 5-flavour number scheme predictions at NLO accuracy agree better with data than 4-flavour number scheme ones. The 4-flavour number scheme predictions underestimate data in events with at least one b-jet.


Introduction
The measurement of the production rate of a Z boson in association with jets originating from b-quarks1 (Z + b-jets) in proton-proton (pp) collisions provides an important test of perturbative quantum chromodynamics (pQCD). Current predictions for Z + b-jets production are known at next-to-leading-order (NLO) accuracy in pQCD, and they can be derived in either a 4-flavour number scheme (4FNS) or a 5-flavour number scheme (5FNS) [1][2][3][4]. In the 4FNS, b-quarks do not contribute to the parton distribution functions (PDFs) of the proton and, in QCD, they only appear in a massive final state due to gluon splitting (g → bb). In the 5FNS, b-quark density is allowed in the initial state via a b-quark PDF, with the b-quark typically being massless. Therefore, in the 5FNS the Z + b-jets cross-section is sensitive to the b-quark PDF and can be used to constrain it. The ambiguity among the schemes is an intrinsic property of the calculation and is expected to reduce with the inclusion of higher order perturbative corrections.
Furthermore, the measurement of Z + b-jets production provides a benchmark to test predictions from Monte Carlo (MC) simulations. These are commonly used to estimate the background contribution of Z + b-jet events to other topologies, such as the production of a Higgs boson decaying into a b-quark pair in association with a Z boson, or in searches for physics beyond the SM with signatures containing leptons and b-jets in the final state.
The Z + b-jets processes occur more rarely than the production of Z-boson events with inclusive jets (Z+jets) and they are more challenging to measure. The b-jets are identified by exploiting the long lifetime of b-hadrons produced in the quark hadronisation, and a higher level of background affects the measurement. The background is mainly composed of events with a Z boson associated with light-flavour jets or c-jets2, misidentified as b-jets, and events from the dileptonic decay of a tt pair. Inclusive and differential cross-sections of Z + b-jets production have been measured in proton-antiproton collisions at the centre-of-mass energy of √ s = 1.96 TeV by the CDF and D0 experiments [5][6][7][8] and at the Large Hadron Collider (LHC) [9] in √ s = 7 TeV pp collisions by the ATLAS and CMS experiments [10][11][12][13][14][15], as well as in √ s = 8 TeV pp collisions by the CMS experiment [16,17]. The CMS experiment also recently released a measurement of the ratio of Z + b-jets to Z+jets cross-sections and the ratio of Z + c-jets to Z + b-jets cross-sections for events with at least one b-jet or one c-jet in √ s = 13 TeV pp collisions [18].
This paper presents a measurement of the inclusive and differential production cross-sections of a Z boson, decaying into electrons or muons, in association with at least one or at least two b-jets using 35.6 fb −1 of pp collision data collected by the ATLAS experiment at √ s = 13 TeV in 2015 and 2016. For events with at least one b-jet, the differential cross-sections are presented as a function of the transverse momentum3 (p T ) and the absolute value of the rapidity (|y|) of the leading b-jet, the p T and the |y| of the Z boson (Z p T and Z |y|), and as a function of observables correlating the Z boson with the leading b-jet, namely the azimuthal angle between them (∆φ Zb ), the absolute value of their rapidity difference (∆y Zb ), and their angular separation (∆R Zb ). For events with at least two b-jets, the differential cross-sections are presented as a function of the p T of the Z boson and as a function of observables built using the two leading b-jets, namely their p T (p T,bb ), their invariant mass (m bb ), p T,bb divided by their invariant mass (p T,bb /m bb ), the azimuthal angle between them (∆φ bb ), the absolute value of their rapidity difference (∆y bb ), and their angular separation (∆R bb ). The higher √ s leads to a large increase in the measured cross-section in comparison with previous ATLAS publications. This allows more extreme regions of phase space to be explored and new measurements to be performed in the rare two-b-jets configuration (i.e. p T,bb and p T,bb /m bb ). Previous ATLAS measurements were compared with MC predictions based on leading-order matrix elements interfaced with a parton-shower simulation, which showed substantial mismodelling. Recent advances in this field permit this paper to compare the data with the latest MC predictions using next-to-leading-order matrix elements, which are expected to provide a better description of the data.
The experimental apparatus is described in Section 2, and details of the data sample and the MC simulations are provided in Section 3. The object definitions and the event selection at detector level are presented in Section 4. Backgrounds that do not contain a real Z boson are estimated via MC simulations and validated in control regions in data or via data-driven techniques, while backgrounds containing a real Z boson and jets not originating from b-quarks are estimated with a fit to data distributions sensitive to the flavour of the jet (flavour fit); both are described in Section 5. Distributions of the kinematic variables are presented in Section 6. After background subtraction, the data are unfolded to particle level in a fiducial phase space, which is detailed in Section 7. Systematic uncertainties in the unfolded data are discussed in Section 8. The results are presented in Section 9, and conclusions are drawn in Section 10.

The ATLAS detector
The ATLAS detector [19] at the LHC covers nearly the entire solid angle around the collision point. 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 |η| < 2.5. The high-granularity silicon pixel detector covers the vertex region and provides four measurements for most tracks, the first hit normally being in the insertable B-layer [20,21]. It is followed by the silicon microstrip tracker, which provides eight measurements per track. These silicon detectors are complemented by the transition radiation tracker (TRT), which enables radially extended track reconstruction up to |η| = 2.0. The TRT also provides electron identification information based on the fraction of hits (typically 30 in total) with an energy deposit above the transition-radiation threshold.
The calorimeter system covers the pseudorapidity range |η| < 4.9. Within the region |η| < 3.2, electromagnetic calorimetry is provided by barrel and endcap high-granularity lead/liquid-argon (LAr) calorimeters, with an additional thin LAr presampler covering |η| < 1.8 to correct for energy loss in material upstream of the calorimeters. Hadronic calorimetry is provided by the steel/scintillator-tile calorimeter, segmented into three barrel structures within |η| < 1.7, and two copper/LAr hadronic endcap calorimeters. The solid angle coverage is completed with forward copper/LAr and tungsten/LAr calorimeter modules optimised for electromagnetic and hadronic measurements, respectively.
The muon spectrometer (MS) comprises separate trigger and high-precision tracking chambers measuring the deflection of muons in a magnetic field generated by the superconducting air-core toroid magnets. The field integral of the toroid magnets ranges between 2.0 and 6.0 T m across most of the detector. The precision chambers cover the region |η| < 2.7 with three layers of monitored drift tubes, complemented by cathode-strip chambers in the forward region, where the background is highest. The muon trigger system covers the range |η| < 2.4 with resistive-plate chambers in the barrel, and thin-gap chambers in the endcap regions.
Interesting events are accepted by the first-level trigger system implemented in custom hardware, followed by selections made by algorithms implemented in software in the high-level trigger [22]. The first-level trigger accepts events from the 40 MHz bunch crossings at a rate below 100 kHz, which the high-level trigger further reduces in order to record events to disk at about 1 kHz rate.

Data set description
The data used in this measurement were recorded in 2015 and 2016 with the ATLAS detector at the LHC in pp collisions at √ s = 13 TeV. The candidate events were selected by either a single-electron or single-muon trigger that imposed a minimum transverse energy (transverse momentum) threshold for the electron (muon) channel and quality and isolation requirements, which depended on the LHC running conditions. The threshold in 2015 was 24 (20) GeV for the electrons (muons), satisfying loose isolation requirements. Due to the higher instantaneous luminosity in 2016, the threshold was increased to 26 GeV for both the electrons and the muons, and a more restrictive isolation requirement was imposed on both leptons along with more restrictive identification requirements for electrons. Triggers with higher thresholds but with no isolation requirement or with loosened identification criteria were also used to increase the efficiency. Crossings of proton bunches occurred every 25 ns, the collisions achieved a peak instantaneous luminosity of 1.37 × 10 34 cm −2 s −1 , and the mean number of pp interactions per bunch crossing (pile-up) was µ = 24. After applying criteria to ensure good ATLAS detector operation, the total integrated luminosity amounts to 35.6 fb −1 . The uncertainty in the combined 2015-2016 integrated luminosity is 2.1% [23], obtained using the LUCID-2 detector [24] for the primary luminosity measurements.

Simulated event samples for signal and background processes
MC simulations are used to describe signal events, to estimate the contribution of background processes, to unfold the data yield to the particle level, to estimate systematic uncertainties, and to compare predictions with the unfolded data distributions. An overview of all signal and background processes and the generators used for the production of nominal results is given in Table 1 together with the theory uncertainty in the normalisation cross-sections corresponding to QCD scale variations.
Inclusive Z(→ , = e, µ) production in association with both light-and heavy-flavour jets was simulated using the S v2.2. 1 [25] generator. In this set-up, matrix elements at NLO for up to two partons, and matrix elements at LO for up to four partons, were calculated with the Comix [26] and OpenLoops [27,28] libraries. They were matched with the S parton shower [29] using the MEPS@NLO prescription [30][31][32][33]. S uses the 5FNS with massless band c-quarks in the matrix element, but massive quarks in the parton shower. Samples were generated using the NNPDF3.0nnlo PDF set [34], along with the dedicated set of tuned parton-shower parameters developed by the S authors. In Section 9, where several predictions are compared with the unfolded data, these samples are shown with their uncertainties and are referred to as S 5FNS (NLO). Uncertainties from missing higher orders are evaluated [35] using seven variations of the QCD factorisation and renormalisation scales in the matrix elements by factors of 0.5 and 2 and avoiding variations in opposite directions.
Additional Z(→ ) samples were produced with the LO matrix-element generator A v2.14 [36], interfaced with P v6.426 [37] to model parton showers, using the parameter values of the Perugia2011C tune [38] for simulating the underlying event, and the CTEQ6L1 PDF set [39]. Matrix elements were calculated for up to five partons, and merged using the MLM prescription [40] with a matching scale of 15 GeV. A uses the 4FNS with massive band c-quarks in the matrix element and in the parton shower of P . The matrix elements for the production of Z + bb and Z + cc events are explicitly included and a heavy-flavour overlap procedure is used to remove the double counting, between the matrix element and the parton shower, of heavy quarks from gluon splitting. The properties of band c-hadron decays were simulated with EvtGen v1.2.0 [41], as was done in all generated samples where the parton shower was simulated with P . Photos++ v3.52 [42,43] was used to simulate QED final-state radiation (FSR). The A samples are used in the analysis to estimate systematic uncertainties in the unfolding procedure and in backgrounds containing a genuine Z boson. In Section 9 these samples are referred to as A + P 6 4FNS (LO). Samples of Z(→ ττ), W(→ ν), and W(→ τν) events were simulated with S , using the same set-up adopted for the signal samples.
The Z-boson and W-boson samples are normalised to the inclusive next-to-next-to-leading-order (NNLO) cross-section predictions provided by the FEWZ 3.1 program [44][45][46][47] with the CT14 PDF set. The K-factor applied to the Z samples to match the NNLO prediction is 0.975 for S and 1.196 for A .
The production of tt events with at least one W boson decaying leptonically was modelled using the P -B [48][49][50][51] v2 generator at NLO with the NNPDF3.0NLO [34] PDF set. The h damp parameter, which regulates the high-p T emissions against which the tt system recoils, is set to 1.5 m top [52]. The events were interfaced with P v8.230 [53] using the A14 tune [54]. The tt sample is normalised to the theory prediction at NNLO in QCD including the resummation of next-to-next-to-leading logarithmic (NNLL) soft-gluon terms [55][56][57][58][59][60][61]. Four additional tt samples were simulated to evaluate the uncertainty in this process. One sample was produced with M G 5_ MC@NLO [62] and the same parton-shower model of the nominal tt sample in order to estimate the uncertainty due to the modelling of the hard scattering process. A second P -B sample showered with H 7.13 [63,64] was generated to evaluate the uncertainty due to the modelling of the parton shower and hadronization processes. A third sample was produced to simulate higher energy radiation with the factorisation and renormalisation scales changed by a factor of 0.5 while simultaneously increasing the h damp value to 3.0 m top and using the 'Var3c up' variation from the A14 tune. The last sample simulates the lower energy radiation. It was generated with the renormalisation and factorisation scales varied by a factor of 2.0 while keeping the h damp value at 1.5 m top and using the 'Var3c down' variation in the parton shower. The last two samples are also used to estimate the impact of FSR with parton-shower weights that vary the renormalisation scale for QCD emission in the FSR by factors of 0.5 and 2.0.
Single-top-quark events in the Wt-, sand t-channels were generated using the P -B v1 generator interfaced with P v6.4 [37]; the latter simulates parton showers, fragmentation, and the underlying event using the Perugia 2012 tune [38]. The CT10 PDF set was used [65]. The single-top samples for the tand s-channels are normalised to cross-sections from NLO predictions [66,67], while the Wt-channel sample is normalised to cross-sections from approximate NNLO predictions [68].
Diboson processes (WW, W Z, and Z Z) with one of the bosons decaying hadronically and the other leptonically were generated using S v2.
NNLO QCD + NLO EW up to one parton at NLO and up to three additional partons at LO. The samples are normalised to the NLO predictions [69].
Simulated events for qq → V H(→ bb) with V = W or Z plus zero or one jet production at NLO were generated with the P -B v2 + GoSam + MiNLO generator [51,[70][71][72] with the NNPDF3.0NLO PDF set. The contribution from gg → Z H(→ bb) production was simulated using the LO P -B v2 matrix-element generator. The samples of simulated events include all final states where the Higgs boson decays into bb and the vector boson into a leptonic final state. The mass of the Higgs boson is set to 125 GeV and the H → bb branching fraction is set to 58%. The qq → V H(→ bb) cross-section is calculated at NNLO (QCD) and NLO (EW), while the gg → Z H cross-section is calculated at NLO+NLL (QCD). Generated events were processed with the ATLAS detector simulation [76], based on G 4 [77], to simulate the detector response to final-state particles. To account for the effects of pile-up, multiple overlaid pp collisions were simulated with the soft QCD processes of P v8.186 using the A2 tune [78] and the MSTW2008LO PDF set [79]. The distribution of the average number of interactions per bunch crossing in the simulation is weighted to reflect that in the data. Simulated events are processed with the same reconstruction algorithms as for the data.

Theoretical predictions
In addition to the particle-level predictions from the S and A samples described above, unfolded results from data are compared with six other predictions listed in Table 2.
Two particle-level predictions (using specific parton-shower and matching predictions) were produced with the S v2.2.7 generator using NLO matrix elements [80]. The first sample, referred to as S Z 4FNS (NLO), includes Z + bb events generated in the 4FNS at NLO with massive b-quarks. It is interesting to compare this sample, which contains two b-quarks in the matrix elements, with the unfolded data even in the case of distributions with at least one b-jet, to understand if there are regions of the phase space that can be described with such a configuration. The second sample, referred to as S F 4FNS+5FNS (NLO), contains the matrix elements at NLO for up to two partons, and matrix elements at LO for up to three partons. It includes both Z + bb events generated in the 4FNS at NLO with massive b-quarks, and Z+jets events generated in the 5FNS at NLO. They are combined according to the procedure described in Ref. [81]. The combination is achieved by means of a dedicated heavy-flavour overlap removal procedure, the fusing technique, that acts as an additional step after the multijet merging algorithms. This procedure combines the advantages of inclusive 5FNS calculations with the higher precision of 4FNS calculations in regions of phase space where the b-quark mass sets a relevant scale. The two S samples use the NNPDF3.0nnlo PDF set with α S (m Z ) = 0.118 and the corresponding number of active quark flavours. Masses of cand b-quarks are taken into account in the parton shower in all S samples.
Results are also compared with predictions from the LO matrix-element generator M G 5_ MC@NLO v2.2.2 [62] interfaced with P v8.186 [53] with the A14 tune [54] to model the parton shower and underlying event. The matrix element includes up to four partons. Additional jets are produced by the parton shower, which uses the CKKW-L merging procedure [82], with a matching scale of 30 GeV. M G 5_ MC@NLO uses the 5FNS with massless band c-quarks in the matrix element, and massive quarks in the parton shower. The NNPDF3.0nlo PDF set is used with α S (m Z ) = 0.118. This prediction is referred to as MG MC + P 8 5FNS (LO).
Two additional predictions were produced with M G 5_ MC@NLO v2.6.2, using matrix-element calculations with NLO accuracy. The first sample includes Z+jets events generated in the 5FNS with up to one parton at NLO, and massless band c-quarks; the second sample includes Z + bb events generated in the 4FNS at NLO, and massive b-quarks. Both samples were generated using the NNPDF3.0nnlo PDF set with α S = 0.118. They were interfaced to the P v8.186 parton shower using the FxFx merging scheme [83], with a matching scale of 25 GeV. As in the previous case, massive cand b-quarks are produced in the parton shower. The first sample is referred to as MG MC + P 8 5FNS (NLO); the second is referred to as MG MC + P 8 Z 4FNS (NLO).

An additional A
prediction is used to test the sensitivity of the measurements to the parton structure of the proton. The A samples presented in Section 3.2 are reweighted to the NNPDF3.0lo PDF set, using the prescriptions reported in Ref. [84]. These predictions are referred to as A + P 6 (rew. NNPDF3.0lo). The predictions of LO MC generators, such as A + P 6 4FNS (LO) and MG MC + P 8 5FNS (LO), with up to four or five partons in the matrix element, are still an interesting case to study as they allow comparison with the predictions of MC generators at NLO accuracy and with a smaller number of partons in the matrix element. Furthermore, they provide a benchmark in common with past analyses, such as in Ref. [11].

Event selection
Events selected in this analysis are required to have a signature consistent with a Z boson, decaying into two electrons or two muons, in association with at least one or at least two b-jets. Candidate events are required to have a primary vertex (PV), defined as the vertex with the highest sum of track p 2 T with at least two associated tracks measured in the ID (ID tracks), each with p T > 400 MeV.   To select leptons originating from the primary pp interaction, the lepton tracks are required to have a longitudinal impact parameter (z 0 ) satisfying |z 0 sin(θ)| < 0.5 mm relative to the PV. The transverse impact parameter significance of the electron (muon) candidates must satisfy d 0 /σ d 0 < 5 (3). In order to further suppress leptons from non-prompt processes or leptons from hadrons in jets, both the electron and muon candidates are required to satisfy p T -dependent cone-based isolation requirements [86], which use information from ID tracks. The isolation requirements are set so that the scalar sum of the transverse momenta of the tracks in the isolation cone4 around the lepton is less than 6% of the lepton p T .
Jets are reconstructed, using the anti-k t algorithm [87, 88] with radius parameter R = 0.4, from topological clusters of energy deposits in the calorimeter [89]. Jets are calibrated using a simulation-based calibration scheme, followed by in situ corrections to account for differences between simulation and data [90]. Events with jets arising from detector noise or other non-collision sources are discarded [91]. Furthermore, to eliminate jets containing a large energy contribution from pile-up, jets with p T < 60 GeV and |η| < 2.4 are required to have a significant fraction of their tracks with origin compatible with the primary vertex, as defined by a jet vertex tagger discriminant (JVT) [92]. Selected jets must have p T > 20 GeV and rapidity |y| < 2.5.
An overlap removal procedure is applied to electron, muon and jet candidates to prevent double counting. Any jet whose axis lies within ∆R = 0.2 of an electron is removed. If a jet is reconstructed within ∆R = 0.2 of a muon and the jet has fewer than three associated tracks or the muon energy constitutes most of the jet energy, then the jet is removed. Any electron or muon of a given p T reconstructed within ∆R = min(0.4, 0.04 + 10 GeV/p T ) of the axis of any surviving jet is removed. Jets that survive the overlap removal procedure are removed if they are within ∆R = 0.4 of the selected leptons.
The b-jets, defined as the jets containing at least one b-hadron, are identified using a multivariate algorithm, MV2c10 [93,94]. This algorithm uses the impact parameter and reconstructed secondary vertex information of the tracks associated with the jets. Its output lies in the range [−1, +1]. A value close to +1 denotes a higher probability for the jet to be a b-jet. The b-jet candidates are selected if their MV2c10 output is greater than 0.8244. This selection corresponds to an efficiency of 70% for selecting jets containing b-hadrons, and misidentification rates of 0.26% and 8.3%, respectively, for light-flavour (u-, d-, s-quark and gluon) jets and c-jets, as estimated from a sample of simulated tt events. Other working points are defined by different b-tagging discriminant output thresholds; they are used to define control regions and to define the bins used in the flavour fit, as detailed in Section 5.1.
In simulation, reconstructed jets are labelled as b-jets if they lie within ∆R = 0.3 of one or more weakly decaying b-hadrons with p T > 5 GeV. Reconstructed jets not identified as b-jets are considered to be c-jets if they lie within ∆R = 0.3 of any c-hadron with p T > 5 GeV. All other jets are classified as light-jets. Simulated Z+jets events are sequentially categorised depending on the labels of the jets, starting from b-jets, as follows: Z + b when they have exactly one b-jet, Z + bb when they have at least two b-jets, Z + c when they have at least one c-jet, Z + l when they have only light-jets. A similar classification is adopted for simulated W+jets events. In the distributions with at least one b-jet, the sum of Z + b and Z + bb samples is used to define the signal, and the Z+jets background is constituted by the sum of the Z + c and Z + l samples. In the distributions with at least two b-jets, the Z + bb samples alone constitute the signal, while the sum of the Z + b, Z + c, and Z + l samples form the Z+jets background.
The missing transverse momentum (E miss T ), which may correspond to a neutrino escaping interaction with the detector, is defined as the negative vector sum of the transverse momentum of all identified hard physics objects (electrons, muons, jets), as well as an additional track-based soft term defined in Ref. [95].
Events are required to have exactly two leptons5 of the same flavour (ee or µµ) but of opposite charge with their dilepton invariant mass in the range 76 GeV< m <106 GeV. Events with p T < 150 GeV must also have E miss T < 60 GeV. The requirement on the E miss T value reduces by about 55% the background from tt events with dileptonic decay, while the signal is reduced by about 5%. Events passing the above selection and having at least one or at least two jets belong to the region referred to as the pre-tag region. The signal region is a subset of the pre-tag region. Events belonging to the signal region are assigned to two regions: those with at least one b-jet, referred to as the 1-tag region; and those with at least two b-jets, referred to as the 2-tag region, which is a subset of the 1-tag region.
A summary of the object selection and the event selection used in the analysis to define the signal regions and the validation regions for the main backgrounds, which are presented in Section 5, is given in Table 3. 5 At least one of the lepton candidates is required to match the lepton that triggered the event.

Correction factors applied to simulation and corresponding uncertainties
Corrections are applied to simulated samples in order to ensure that the object selection efficiencies and the energy and momentum calibrations agree with data within the uncertainties associated with the corrections.
The electron and muon trigger efficiencies are estimated in data and simulation in order to determine simulation-to-data correction factors and their corresponding uncertainties. The average per-event correction factor is about 0.98 (0.93) for electron (muon) triggers; they are known with an uncertainty below 1% [85,86]. Corrections to efficiencies for lepton reconstruction, identification, isolation and association with the PV in simulated samples are derived from data. Each per-lepton correction factor is close to unity and known with a precision that is better than 1% in the kinematic range considered [85,86].
The energy scale of the electrons and the momentum scale of the muons in simulation are adjusted with correction factors that deviate from unity at the per-mil level and the resolutions are adjusted with correction factors that deviate from unity at the per-cent level in order to match lepton p T and m distributions in data; the corresponding uncertainties are negligible.
The jet energy scale (JES) is calibrated on the basis of the simulation including in situ corrections obtained from data [90]. The JES uncertainties are estimated using a decorrelation scheme comprising a set of 21 independent parameters, the largest of which may reach several per cent in specific corners of the phase space. The jet energy resolution (JER) uncertainty is derived by over-smearing the jet energy in the simulation by about 4% at p T = 20 GeV to about 0.5% at a p T of several hundred GeV. Simulation-to-data corrections and relative uncertainties are also applied to adjust the efficiency of the JVT requirement following the prescriptions of Ref. [96]. The uncertainty in the scale and resolution of E miss T is estimated by propagating the uncertainties in the transverse momenta of reconstructed objects and an uncertainty to account for soft hadronic activity in the event, as described in Ref. [95].
Flavour-tagging efficiencies in simulation are scaled to match those measured in data for jets of all flavours as a function of the different b-tagging discriminant output thresholds, and of the jet p T (and η for light-jets), using weights derived from control samples enriched in jets of each flavour [97]. In the case of b-jets, correction factors and their uncertainties are estimated from data using dileptonic tt events [97]. In the case of c-jets, they are derived using jets from W-boson decays in tt events [98]. In the case of light-flavour jets, correction factors are derived using dijet events [99].
The correction factors for b-jets are close to unity. The uncertainties, described by a set of 28 independent parameters, are as low as 3% for jet p T of about 60 GeV, but reach 10% for jet p T of about 20 GeV and up to 20% beyond 300 GeV. The correction factors for c-jets range from about 1.2 to about 1.6. Uncertainties, described by a set of 28 independent parameters, are about 20%-30% in the bulk of the phase space, but up to 100% for large jet p T and for the b-tagging discriminant output threshold closest to +1. The correction factors for light-jets range from about 2 to about 3, with uncertainties described by a set of 36 independent parameters and ranging from 50% to 100%. An additional uncertainty of 30% is applied to the efficiency of b-tagging for simulated jets originating from pile-up interactions.
A variation in the pile-up reweighting of simulated events (referred to as pile-up uncertainty) is included to account for the uncertainty in the ratio of the predicted and measured inelastic cross-sections in the fiducial volume [100].

Background estimation
The main background in the 1-tag region is constituted by events with a Z boson produced in association with jets, where either a light-jet or a c-jet is misidentified as a b-jet; it is determined using a fit to data as detailed in Section 5.1. Dileptonic tt events dominate in the 2-tag region. Smaller background contributions from the production of dibosons, a Higgs boson, a single top quark, a Z → ττ, or a W → ν are estimated using simulation, as described in Section 3.2. Uncertainties in the normalisation cross-section of these predictions range from 4% to 6% depending on the process, as detailed in Table 1. Background contributions from multijet events are estimated with a data-driven technique and found to be negligible, as described below.
The tt contribution is estimated using simulated events generated with P -B + P normalised to the theoretically predicted cross-section, as discussed in Section 3.2. An uncertainty of about 6% is assigned to the inclusive tt cross-section (see Table 1), following the variation of the renormalisation and factorisation scales by a factor of 2.0, and the variation of the PDFs within their uncertainties. In addition, systematic uncertainties in the modelling of the distributions are derived by comparing the predictions from the nominal tt sample with the ones from the alternative samples described in Section 3.2.
The modelling of tt production in the simulation is validated using a tt-enriched region, which is selected by requiring that events have two leptons of different flavour (eµ); all other selections are the same as in the signal region. As an example, Figure 1 shows the p T,bb and the m bb distributions for events with at least two b-jets. The total background from top quarks is the sum of tt and single-top events, where the latter are about 3% of the tt component in the validation region, and other backgrounds are negligible. Data and simulation agree well within the uncertainties which account for both the yield and shape uncertainties of simulated tt events and the statistical uncertainties of predictions and data. Background contributions from multijet events in the electron and muon channels are estimated using a data-driven technique. Multijet-enriched control regions without b-tag and m requirements are used to derive the expected shape of this background. In the electron channel, the multijet-enriched control region is defined by applying the full signal event selection except for the electron identification and the d 0 /σ d 0 cuts, and inverting the isolation selection for both electron candidates. In the muon channel, the multijet-enriched control region is defined by applying the full signal event selection but requiring both muon candidates to have the same charge. In both channels, contributions from non-multijet sources in the control regions are estimated from simulation and subtracted from the data, with the remaining distributions used as shape templates. A fit of the m distribution to data is then performed within the window of 60 GeV < m < 160 GeV in the one-jet and two-jets pre-tag regions separately and leaving the normalisation of the signal and of the multijet background templates free to float in the fit, while the normalisation of the other processes is fixed in the fit. The multijet background estimate in the pre-tag region is then extrapolated to the two signal regions using normalisation factors equal to the fraction of events in the multijet control region that satisfy the 1-tag and 2-tag requirements. Contributions from non-multijet processes are subtracted before estimating this fraction. Systematic uncertainties are assessed by varying the m range and the binning of the fit, excluding the Z-boson peak from the fit, performing the fit in the tagged regions in place of the pre-tag ones, and by allowing the other processes to be varied independently in the fit. The estimated size of the multijet background is consistent with zero within the statistical uncertainty even after considering all sources of systematic uncertainty. It is therefore neglected in the analysis.

Extraction of the cross-section for Z-boson production in association with light-jets and c-jets
The flavour fit used for the extraction of the yields of Z + light-jets and Z + c-jets backgrounds for the 1-tag and 2-tag selections is a maximum-likelihood fit to data based on flavour-sensitive distributions. The fit is done simultaneously in the electron and muon channels with templates derived from simulation.
In the 1-tag region, the b-tagging discriminant output of the leading b-jet is used as the flavour-sensitive distribution. This observable for events belonging to the signal region is distributed into three intervals that define the bins of the discriminant output distribution. Each bin corresponds to a certain range of b-tagging efficiency. The bins are numbered from 1 to 3, corresponding respectively to efficiencies of 60%-70% (bin 1), 50%-60% (bin 2) and <50% (bin 3) as estimated from simulated tt events. The light-flavour jet (c-jet) misidentification rates for the three bins are respectively 0.195% (5.4%), 0.048% (1.96%), and <0.017% (<0.94%). The signal template is built with simulated Z+ ≥ 1b events. The template shapes of the Z + l and Z + c samples are very similar (as shown in Figure 2), hence those samples are combined to form a single template. All non-Z+jets backgrounds are combined into a single template, determined from the sum of their predicted contributions. The normalisations of the signal and of the Z+jets background are free to float in the fit, while the normalisation of the sum of the non-Z+jets backgrounds is fixed to their estimate.
In the 2-tag region the combination of the three bins of the b-tagging discriminant outputs of the leading and sub-leading b-jets produces a distribution with six bins that is used for the fit to data. The signal template is built with simulated Z + bb events. Templates built with Z + b, Z + c and Z + l simulated events are combined into a single template. Because of the large rejection of light-flavour jets achieved in the 2-tag selection, the simulated Z + l events in this region are not subjected to the b-tagging requirement. Instead they are weighted by a per-event probability that the jets pass the two-b-tags selection. This probability is computed on the basis of the per-jet probabilities, which are assumed to be independent of each other [101].
As for the fit in the 2-tag region, the normalisations of the signal and of the Z+jets background are also free to float, while the normalisation of the other backgrounds is fixed to their estimate. Tables 4 and 5 show the normalisation scale factors in the 1-and 2-tag regions obtained from the fit, together with the post-fit yields for the signal and Z+jet background samples generated with S or A . There is good agreement between the sum of the signal and background post-fit yields of S and A . The differences between S and A in the modelling of the Z+jet backgrounds after the flavour fit are taken into account in the systematic uncertainties as described below. The statistical uncertainty is estimated with pseudo-experiments. The Z+jets backgrounds predicted by S and corrected for the normalisation factor obtained from the fit are used as the nominal estimate in this analysis. Systematic uncertainties due to the object selection efficiencies and calibrations, discussed in Section 4.1, affect the normalisation and the shape of Z+jets backgrounds. They are assessed by repeating the fit with the templates varied according to each of the systematic uncertainties. The fit is also repeated for each of the uncertainties affecting the tt and other backgrounds detailed above. An additional systematic uncertainty (referred to as the flavour fit uncertainty) in the normalisation of the Z+jets backgrounds is estimated by repeating the fit after separating the Z + c from the Z + l template in the 1-tag region, and after separating the Z + b from the Z + c and Z + l templates in the 2-tag region. An uncertainty affecting the shape and rate of the Z+jets background is derived by taking the difference between the post-fit Z+jets background evaluations using S and A samples. Another uncertainty accounts for potential jet-jet correlations that are not covered by the per-event weighting procedure which mitigates the large statistical fluctuations in the 2-tag region for Z + l. A 20% uncertainty is derived by taking the largest difference between the double-tagged event yields obtained with or without the weighting procedure being applied to simulated samples of Z + bb, Z + cc, W + bb, and W + cc.6. These samples suffer less from statistical limitations. The test is done with both the S and A samples.   The post-fit estimate of the S Z+jets background is validated in a region defined by applying the full signal event selection with the exception of b-tagging requirements. Events with at least one b-jet, with the b-tagging discriminant output in the b-jet efficiency range of 70%-77% and light-flavour jet (c-jet) misidentification rates of 0.51% (7.7%), are selected to provide a sample enriched in c-jets and light-flavour jets. As an example, Figure 3 shows the p T of the leading b-jet and the p T of the Z boson in this region. The Z + l and Z + c backgrounds constitute 50% and 28% of the total prediction, respectively. Agreement between data and estimated backgrounds is observed within uncertainties. These include the uncertainties due to the flavour fit and b-tagging efficiency, and the statistical uncertainties of the predictions and data.
The normalisation factors of the signal samples, shown in Tables 4 and 5, are applied in this section to demonstrate the robustness of this procedure, while in the following sections, post-fit normalisation factors are applied only to Z+jets background.

Kinematic distributions
After the signal selection criteria are applied, the measured and expected distributions are compared at the detector level. The Z+jets background is shown for the normalisation factors derived from the flavour fit. Pre-fit distributions are used for the signal samples.     Figure 5 shows the p T of the Z boson and the ∆R bb distributions for events in the 2-tag region. The uncertainty bands include the statistical uncertainties of the simulated sample, the event-selection uncertainties described in Section 4 (omitting the common luminosity uncertainty), and the background uncertainties described in Section 5. Both generators do not describe precisely the data in the full range of the measurement, although the S generator provides the best agreement with data. Entries / GeV  The total numbers of selected events in data and in predictions are presented in Table 6, together with the prediction of each process, expressed as a fraction of the total number of predicted events.  Table 6: The expected size of the signal and backgrounds, expressed as a fraction of the total number of predicted events for inclusive b-jet multiplicities for the signal selection. The signal and Z+jets background predictions are from the S generator, with the Z+jets background estimate obtained after applying the normalisation scale factors obtained from the flavour fit. The total numbers of predicted and observed events are also shown. The uncertainty in the total predicted number of events is statistical only.   The lower panels display the ratio of the predictions for signal plus background to data using either S (red) or A + P 6 (blue) as the signal simulation. The statistical uncertainty of the data is shown as black error bars and the total uncertainty of the prediction as a hatched band. The latter consists of the statistical uncertainty and all systematic uncertainties from the predictions.  Figure 5: Distribution of events passing the signal selection as a function of p T, Z (left) and ∆R bb (right) for events with at least two b-jets. The lower panels display the ratio of the predictions for signal plus background to data using either S (red) or A + P 6 (blue) as the signal simulation. The statistical uncertainty of the data is shown as black error bars and the total uncertainty of the prediction as the hatched band. The latter consists of the statistical uncertainty and all systematic uncertainties from the predictions.

Correction to particle level
The signal event yields are determined by subtracting the estimated background contributions from the data. The resulting distributions are corrected for detector-level effects to the fiducial phase space at particle level defined in Table 7. The procedure, based on simulated samples, corrects for Z-boson, jet, and b-jet selection efficiencies, resolution effects, and small differences between the fiducial and detector-level phase spaces. The pre-fit distributions of the S signal samples are used to perform the unfolding procedure. The signal samples for the simulation of Z events with at least one or at least two b-jets are defined in Section 4. Particle-level objects are selected with requirements close to the corresponding requirements for reconstructed signal candidate objects, in order to limit the dependence of the measurement on theoretical predictions. In this definition, the lepton kinematic variables are computed using final-state leptons from the Z-boson decay. Photons radiated by the boson decay products within a cone of size ∆R = 0.1 around the direction of a final-state lepton are added to the lepton, and the sum is referred to as the 'dressed' lepton. Particle-level jets are identified by applying the anti-k t algorithm with R = 0.4 to all final-state particles with a lifetime longer than 30 ps, excluding the dressed Z-boson decay products. A jet is identified as b-tagged if it lies within ∆R = 0.3 of one or more weakly decaying b-hadrons with p T > 5 GeV. If a b-hadron matches more than one jet, only the closest jet in ∆R is labelled as a b-jet.
The correction of differential distributions is implemented using an iterative Bayesian method of unfolding [102] with two iterations. Simulated events are used to generate a response matrix for each distribution to account for bin-to-bin migration effects between the detector-level and particle-level distributions. The matrix is filled with the events that pass both the detector-level and particle-level selections. The particlelevel prediction is used as the initial prior to determine the first estimate of the unfolded data distribution. For the second iteration, the new estimate of unfolded data is obtained using the background-subtracted data and an unfolding matrix, which is derived on the basis of the Bayes' theorem from the response matrix and the current prior. The background-subtracted data are corrected for the expected fraction of events that pass the detector-level selection, but not the particle-level one, before entering the iterative unfolding. For each bin of each differential distribution, the unfolded event yields are divided by the integrated luminosity of the data sample and by the bin width, to obtain the cross-section measurement. The differential cross-section measurement of a given observable in the i-th bin is given by: where L is the integrated luminosity, i is the reconstruction efficiency in i-th bin, N bsD j is the number of background-subtracted data events in the j-th bin, f j is the factor that corrects for unmatched events in the j-th bin, and U i j is the element (i, j) of the unfolding matrix calculated after two iterations, using the updated prior from the first iteration and the response matrix.
The measurement of the inclusive cross-section for Z-boson events with at least one or at least two b-jets is obtained by applying a particle-level correction to the number of events in data with at least one or at least two b-jets, after background subtraction. The correction, which is applied as a divisor of the background-subtracted data, is derived from the ratio of the total number of reconstructed events in the detector-level phase space to the number of particle-level events in the fiducial phase space. It is 0.399 ± 0.001 for Z-boson events with at least one b-jet and 0.258 ± 0.002 for Z-boson events with at least two b-jets, using S signal samples and quoting the statistical error.
Since the electron and muon decay channels are combined to increase the precision of the signal fits to data, the corrections and response matrices are made using electron and muon signal samples to obtain combined particle-level yields. To validate this procedure, the analysis is performed for each of the two lepton channels separately. The results obtained from the individual channels are compatible within 1.4σ and 1.6σ with the inclusive cross-section of Z-boson events with at least one b-jet and at least two b-jets, respectively. This comparison uses only the sum in quadrature of the statistical and uncorrelated systematic uncertainties. The differential cross-section measurements in the two channels also agree over the full range of each distribution. Table 8 summarises the systematic uncertainties of the inclusive Z + b-jets cross-sections in the one-and two-b-tag regions. Figure 6 shows as an example the breakdown of the systematic uncertainties in the cross-section as a function of Z-boson p T for events with at least one b-jet and as a function of ∆R bb for events with at least two b-jets.

Uncertainties in the cross-section measurements
The systematic uncertainties in the cross-sections associated with the detector-level uncertainty sources described in Section 4.1 are derived for each observable by propagating systematic shifts from each source through both the response matrices (unfolding factor) and the subtracted background contributions into the unfolded data for the differential (inclusive) cross-section measurements. The dominant source of uncertainty is the modelling of the b-tagging efficiency. Its impact on the inclusive cross-section ranges from 7.0% for Z-boson events with at least one b-jet to 14% for Z-boson events with at least two b-jets. Its effect on differential cross-section measurements ranges from 5% to 10% for Z-boson events with at least one b-jet and from 10% to 15% for Z-boson events with at least two b-jets. The impact of the mistag rate of cand light-jets is smaller; it is 2.4% for Z-boson events with at least one b-jet and 1% for Z-boson events with at least two b-jets.
The uncertainty from each background source is determined by applying shifts to the subtracted background contributions and to the nominal response matrices or unfolding factors. The sources of uncertainty considered for Z + l and Z + c (and Z + 1b in the Z + ≥ 2b-jets measurement), tt and single-top, diboson and other minor backgrounds are described in Section 5. The dominant uncertainty in the background to events with at least one b-jet originates from Z+jets events. This uncertainty contributes 4.5% to the uncertainty in the inclusive cross-section. An uncertainty of 3.7% derives from the difference between the modelling in A and S , while 2.6% is due to the flavour fit uncertainty. The impact of this uncertainty on the differential cross-sections ranges from a few per cent up to 25% in the extreme corners Source of uncertainty Total [%] 10 16 Table 8: Relative systematic uncertainties in the measured production cross-sections of Z(→ ) + ≥ 1 b-jet and Z(→ ) + ≥ 2 b-jets events. The "Jet" term includes the JES, JER and JVT uncertainties. The "Lepton" term includes the lepton trigger, efficiency, scale and resolution uncertainties. The "Z + c Z + l backgrounds" term also includes the Z + 1b background in the Z + ≥ 2 b-jets measurement.  Table 8 are shown in different colours.
of the phase space. For a Z-boson p T value of about 500 GeV, the difference between the modelling in A and S contributes 18% to this uncertainty, and the flavour fit uncertainty is 12%.
In contrast, the uncertainty in the estimation of background from tt events is the dominant source of uncertainty in the background to Z-boson events with at least two b-jets. It contributes 3.8% to the inclusive cross-section and ranges from 1% to 9% in the differential cross-sections.
The uncertainty due to modelling of the Z + b-jets signal samples in the events with at least one and at least two b-jets are also accounted for. This is evaluated for each observable by reweighting the generator-level distribution in the S samples to provide a better description of the data at detector level. The modified S samples are then used to emulate data and are unfolded with the nominal simulated sample. An additional source accounts for the possible mismodelling of an observable that is not one of the unfolded observables (i.e. a hidden variable). This uncertainty is evaluated by reweighting, in the S samples, the generator-level distribution of the leading lepton's p T , which is one of the observables showing the largest mismodelling, to provide a better description of the data at detector level. The modified S samples are used to unfold the data. The effect of the hidden variable's mismodelling is negligible for all considered variables and all bins. A third uncertainty source accounts for the different hadronisation and parton-shower models used for the signal simulation. This uncertainty is evaluated by unfolding the A signal samples, which emulate the background-subtracted data, with the S signal samples. The generator-level distributions from the A samples are first reweighted to agree with S in order to remove effects related to shape differences. The difference between the generator-level distribution and the unfolded A reweighted distribution is taken as the uncertainty. For the inclusive cross-section, the modelling uncertainty is estimated by replacing the unfolding factor computed with S with the one computed with A . The dependence on the size of the simulated sample is derived using pseudo-experiments, and the spread of the results is taken as an uncertainty. The statistical term is typically less than a few per cent. It reaches 5% in the last bin of the ∆R bb distribution and 15% only in the last bin of the ∆y bb distribution.
The total unfolding uncertainty in the inclusive cross-sections is at the level of 4% in each of the two signal regions. In the differential distributions it is less than 5% in the 1-tag region and at a level of 5%-10% in the 2-tag region, except in some bins of the angular variables and in the tail of the p T and m bb distributions, where it reaches 20%.

Results
The inclusive and differential cross-section measurements for Z + ≥ 1 b-jet and Z + ≥ 2 b-jets are shown in Figures 7-15. The statistical uncertainty of the data is propagated through the unfolding by using 1000 pseudo-experiments, repeating the flavour fit for each of them. The statistical uncertainty in the inclusive cross-sections of Z + ≥ 1 b-jet and Z + ≥ 2 b-jets is 0.3% and 0.8% respectively. As mentioned in

Inclusive cross-sections
The measured inclusive cross-sections for Z + ≥ 1 b-jet and Z + ≥ 2 b-jets, shown in Figure 7, are 10.90 ± 0.03(stat.) ± 1.08(syst.) ± 0.25(lumi.) pb and 1.32 ± 0.01(stat.) ± 0.21(syst.) ± 0.04(lumi.) pb, respectively. The 4FNS MC predictions are systematically lower than data in the inclusive one-b-jet case, both for MC generators with LO matrix elements, as implemented in A + P 6 4FNS (LO), and for Z bb predictions at NLO, as implemented in S Z 4FNS (NLO) and MG MC + P 8 Z 4FNS (NLO). The 4FNS predictions agree well with data in the inclusive two-b-jet case. Even though the LO A + P 6 4FNS (LO) underestimates the data, the predictions and data agree within two standard deviations (2σ) of the experimental uncertainty. Use of the NNPDF3.0lo PDF set in A predictions gives better agreement with data because of a higher acceptance in the fiducial region. The 5FNS simulations, in general, adequately predict the inclusive cross-sections for both Z + ≥ 1 b-jet and Z + ≥ 2 b-jets. Overall, this is consistent with the results presented in the ATLAS measurement at √ s = 7 TeV [11].

Differential cross-sections for Z + ≥ 1 b-jet
The differential cross-section measurements for the Z + ≥ 1 b-jet process are shown in Figures 8-11. Each distribution is presented and discussed in detail in this section.
The distributions of the transverse momentum of the Z boson and of the jets probe pQCD over a wide range of scales and provide important input to the background prediction for other SM processes, including Higgs boson production and searches beyond the SM. The differential cross-section as a function of the Z-boson p T for events with at least one b-jet is shown in Figure 8 (left). In the low p T region, up to 100 GeV, where soft radiative effects play a role, all the predicted shapes except that of MG MC + P 8 Z 4FNS (NLO) exhibit trends different from those in the data. Overall, the predictions from S 5FNS (NLO) and S F 4FNS+5FNS (NLO) show the best agreement with data. Predictions from MG MC + P 8 5FNS (LO) and MG MC + P 8 5FNS (NLO) are within the experimental uncertainty band for most of the bins. The harder Z-boson p T in A predictions than in data has already been reported by ATLAS for data collected at √ s = 7 TeV [11]. Figure 8 (right) shows the leading b-jet p T . MG MC + P 8 5FNS (LO) provides a satisfactory description within the uncertainty of the data, while MG MC + P 8 5FNS (NLO) underestimates the data in the high p T region. This region is populated by additional hard radiation, which in MG MC + P 8 5FNS (NLO) is simulated only via parton shower. S 5FNS (NLO) exhibits the best agreement with data. The contrasting behaviour of S F 4FNS+5FNS (NLO), which underestimates the data at high p T , may be interesting to investigate further in the future. The NLO 4FNS predictions of Z bb, as implemented in S and MG MC, show a softer leading b-jet p T , while the inclusive LO 4FNS prediction, as implemented in A , describes the shape of the data quite well despite the large underestimation of the normalisation already discussed for Figure 7.
The distributions of the Z-boson rapidity, the leading b-jet rapidity, and their separation, ∆y Zb , are directly sensitive to the b-quark PDFs and to higher-order diagram contributions, and they may show differences for different flavour schemes. The differential cross-sections as a function of the Z-boson rapidity and of the leading b-jet rapidity for events with at least one b-jet are shown in Figure 9. All MC predictions provide a satisfactory description of the shape of the data. Some modulation relative to data is observed in the leading b-jet |y| distribution, in some cases beyond the experimental uncertainty. Figure 10 (right) shows the differential cross-section as a function of ∆y Zb . S 5FNS (NLO) and S F 4FNS+5FNS (NLO) describe the data quite well, while all other predictions exhibit a slightly smaller rapidity separation than data, even if within the uncertainty of the data. Use of a different PDF set as in A predictions leads to a change in the distribution, but the differences are small compared with the experimental uncertainties.
The distribution of ∆φ Zb is sensitive to the presence of additional radiation in the event. In fixed order calculations of the Z + 1b process, the LO matrix element provides contributions only for ∆φ Zb = π, while the NLO matrix element is the first order which populates the region of ∆φ Zb < π. In MC simulations the region below π is populated via parton shower and via merging of parton shower with multi-parton matrix elements. Therefore the region of small azimuthal separation between the Z boson and the leading b-jet is the most sensitive to additional QCD radiation and soft corrections. It is also sensitive to the presence of boosted particles decaying into a Z boson and b-quarks. The differential cross-section as a function of ∆φ Zb for events with at least one b-jet is shown in Figure 10 (left). The S 5FNS (NLO) generator provides the best agreement with data. S F 4FNS+5FNS (NLO) is still consistent with data within the experimental uncertainty in most of the bins, but a small difference between the two simulations is observed for small values. This result is highly correlated with the difference observed in the leading b-jet p T distribution. It confirms that the current performance of S F 4FNS+5FNS (NLO) in the regime of high-p T jets with a Z boson emitted collinearly is slightly worse than the S 5FNS (NLO) configuration. All MG MC simulations predict too many large azimuthal separations, with a consequent deficit at small angles. Also, in this case the modelling in MG MC + P 8 5FNS (NLO) is slightly worse than in MG MC + P 8 5FNS (LO). The differential cross-section as a function of ∆R Zb , as shown in Figure 11, contains the convolution of effects discussed for the ∆y Zb and ∆φ Zb distributions.

Differential cross-sections for Z + ≥ 2 b-jets
Events with a Z boson produced in association with two b-jets constitute an important background to other SM and beyond-SM processes. Furthermore, they probe the mechanism of a gluon splitting into heavy quarks. The differential cross-section measurements for Z + ≥ 2 b-jet are shown in Figures 12-15. Each distribution is presented and discussed in detail in this section.
The distributions of angular separation between the two leading b-jets allow characterisation of the hard radiation at large angles and the soft radiation for collinear emissions. The differential cross-sections as a function of ∆φ bb and of ∆y bb are shown in Figure 12. Most of the predictions provide satisfactory descriptions of the data within the large experimental uncertainties. Disagreement between data and MG MC + P 8 Z 4FNS (NLO) is observed at low values of ∆φ bb . Mismodelling of ∆y bb is observed for A . This observable has some sensitivity to PDFs, but that is below the experimental uncertainties. The ∆R bb observable is sensitive to the various production mechanisms of the Z bb final state. The region at low ∆R bb is dominated by the production of two b-jets from gluon splitting. Probing this region requires two b-jets in the final state, so it is not sensitive to very small angles of the splitting. The interplay of the modelling of ∆φ bb and ∆y bb in A + P 6 4 FNS (LO) influences the prediction of the ∆R bb distribution shown in Figure 13 (left). All S predictions describe the shape of this observable quite well, featuring a substantial improvement at low ∆R bb relative to the LO version reported by ATLAS using data at √ s = 7 TeV. Overall, this is consistent with the results presented in the ATLAS measurement of gluon-splitting properties at √ s = 13 TeV [11]. MG MC + P 8 Z 4FNS (NLO) presents a large mismodelling at low ∆R bb , which is the part of the phase space dominated by gluon splitting.
The invariant mass of the two leading b-jets is an important observable in the measurement of associated Z H production with Higgs boson decays into bb, and in searches for physics beyond the SM in the same final state. The differential cross-section as a function of m bb for events with at least two b-jets is shown in Figure 13 (right). All S predictions provide a quite good model of the shape of this observable's distribution up to about 300 GeV, while the other predictions show various discrepancies in this region. This is particularly evident for MG MC + P 8 Z 4FNS (NLO), and it is consistent with the mismodelling observed at low ∆R bb , the region dominated by gluon splitting. In the high mass range all predictions underestimate the data, resulting in a sizeable mismodelling. Hence the use of these predictions for the background estimate in searches for physics beyond the SM in this final state could be problematic.
The differential cross-sections as a function of the Z-boson p T and of the p T of the di-b-jet system (p T,bb ) for events with at least two b-jets are shown in Figure 14. Most of the predictions agree with data within the large experimental uncertainties, which are about 25% in most of the bins, and large statistical uncertainties of the predictions, which for some MC samples reach 25% in the highest bins. A shows a harder Z-boson p T spectrum than data, as was observed in the distribution of events with at least one b-jet. The Z bb simulation at NLO with 4FNS, as implemented in MG MC + P 8 Z 4FNS (NLO) and S Z 4FNS (NLO), shows better agreement with data with respect to the p T distributions for events with at least one b-jet, but significant disagreement is still observed. Finally, the ratio of the p T of the di-b-jet system to its invariant mass (p T,bb /m bb ) is sensitive to gluon splitting: a small value indicates a hard splitting and a large value is a consequence of soft splitting. The differential cross-section as a function of p T,bb /m bb is shown in Figure 15

Conclusion
This paper presents a measurement of the cross-sections for Z-boson production in association with one or more b-jets in pp collisions at √ s = 13 TeV. The analysed data correspond to an integrated luminosity of 35.6 fb −1 recorded by the ATLAS detector at the LHC.
The cross-sections are measured using the electron and muon decay modes of the Z boson in a fiducial phase space. In addition to the inclusive cross-sections, differential cross-sections of several kinematic observables are measured, extending the range of jet transverse momenta to higher values than reported in previous ATLAS publications, which used data at lower centre-of-mass energies.
The measurements are compared with predictions from a variety of Monte Carlo generators. In general, 5-flavour number scheme (5FNS) calculations at NLO accuracy predict the inclusive cross-sections well, while inclusive 4-flavour number scheme (4FNS) LO calculations largely underestimate the data. Predictions of Z bb at NLO accuracy agree with data only in the two-b-jets case, and underestimate the data in the case of events with at least one b-jet. Overall, S 5FNS (NLO), a 5FNS generator with matrix elements at NLO for up to two partons and matrix elements at LO for up to four partons, describes the various differential distributions within the experimental uncertainties. A significant discrepancy, common to all generators, is found for large values of m bb . The S F 4FNS+5FNS (NLO) simulation, which combines 4FNS with 5FNS at NLO accuracy using a novel technique, agrees with S 5FNS (NLO), showing that in general at the scales tested by this measurement the effects of this merging are minor. A disagreement of about 20 30% is observed for large values of the leading b-jet transverse momentum, and for small angular separations between the Z boson and the leading b-jet.
The 5FNS simulation with matrix elements for up to four partons at LO, as implemented in MG MC + P 8 (LO), describes the data within the experimental uncertainties in most cases. In some cases this simulation is even better than predictions from MG MC + P 8 5FNS (NLO), which has matrix elements with only one parton at NLO. This indicates the importance of simulations with several partons in the matrix element for a fair description of the data. The pure Z bb simulation at NLO in the 4FNS, as generated by S and MG MC, shows significant deviations from the data even in the two-b-jets configuration, and this is more pronounced in MG MC.
This measurement provides essential input for the improvement of theoretical predictions and Monte Carlo generators of Z-boson production in association with b-jets, allowing a better quantitative understanding of perturbative QCD.    [6] D0 Collaboration, Measurement of the ratio of differential cross sections σ(pp → Z+bjet)/σ(pp → Z + jet) in pp collisions at          [103] ATLAS Collaboration, ATLAS Computing Acknowledgements, ATL-GEN-PUB-2016-002, : https://cds.cern.ch/record/2202407.