1 Introduction

Precise measurements of \(\mathrm {t}\overline{\mathrm {t}}\) production and decay properties [19] provide crucial information for testing the expectations of the standard model (SM) and specifically of calculations in the framework of perturbative quantum chromodynamics (QCD) at high-energy scales. At the energies of the CERN LHC, about half of the \(\mathrm {t}\overline{\mathrm {t}}\) events contain jets with transverse momentum (\(p_{\mathrm {T}}\)) larger than 30\(\,\text {GeV}\) that do not come from the weak decay of the \(\mathrm {t}\overline{\mathrm {t}}\) system [5]. In this paper, these jets will be referred to as “additional jets” and the events as “\(\mathrm {t}\overline{\mathrm {t}}\) +jets”. The additional jets typically arise from initial-state QCD radiation, and their study provides an essential test of the validity and completeness of higher-order QCD calculations describing the processes leading to multijet events.

A correct description of these events is also relevant because \(\mathrm {t}\overline{\mathrm {t}}\) +jets processes constitute important backgrounds in the searches for new physics. These processes also constitute a challenging background in the attempt to observe the production of a Higgs boson in association with a \(\mathrm {t}\overline{\mathrm {t}}\) pair (\(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {H}}\)), where the Higgs boson decays to a bottom (\(\mathrm {b}\)) quark pair (\(\mathrm {b} \overline{\mathrm {b}} \)), because of the much larger cross section compared to the \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {H}}\) signal. Such a process has an irreducible nonresonant background from \(\mathrm {t}\overline{\mathrm {t}}\) pair production in association with a \(\mathrm {b} \overline{\mathrm {b}} \) pair from gluon splitting. Therefore, measurements of \(\mathrm {t}\overline{\mathrm {t}}\) +jets and \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b} \overline{\mathrm {b}}}\) production can give important information about the main background in the search for the \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {H}}\) process and provide a good test of next-to-leading-order (NLO) QCD calculations.

Here, we present a detailed study of the production of \(\mathrm {t}\overline{\mathrm {t}}\) events with additional jets and \(\mathrm {b}\) quark jets in the final state from pp collisions at \(\sqrt{s} = 8\,\text {TeV} \) using the data recorded in 2012 with the CMS detector, corresponding to an integrated luminosity of 19.7 \(\,\text {fb}^\text {-1}\). The \(\mathrm {t}\overline{\mathrm {t}}\) pairs are reconstructed in the dilepton decay channel with two oppositely charged isolated leptons (electrons or muons) and at least two jets. The analysis follows, to a large extent, the strategy used in the measurement of normalized \(\mathrm {t}\overline{\mathrm {t}}\) differential cross sections in the same decay channel described in Ref. [8].

The measurements of the absolute and normalized differential \(\mathrm {t}\overline{\mathrm {t}}\) cross sections are performed as a function of the jet multiplicity for different \(p_{\mathrm {T}}\) thresholds for the jets, in order to probe the momentum dependence of the hard-gluon emission. The results are presented in a visible phase space in which all selected final-state objects are produced within the detector acceptance and are thus measurable experimentally. The study extends the previous measurement at \(\sqrt{s} = 7\,\text {TeV} \) [5], where only normalized differential cross sections were presented.

The absolute and normalized \(\mathrm {t}\overline{\mathrm {t}}\) +jets production cross sections are also measured as a function of the \(p_{\mathrm {T}}\) and pseudorapidity (\(\eta \)) [10] of the leading additional jets, ordered by \(p_{\mathrm {T}}\) . The CMS experiment has previously published a measurement of the inclusive \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b} \overline{\mathrm {b}}}\) production cross section [11]. In the present analysis, the \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b} \overline{\mathrm {b}}}\) and \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b}}\) (referred to as “\(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b} \overline{\mathrm {b}}}\) (\(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b}}\))” in the following) cross sections are measured for the first time differentially as a function of the properties of the additional jets associated with \(\mathrm {b}\) quarks, which will hereafter be called \(\mathrm {b}\) jets. The \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b} \overline{\mathrm {b}}}\) process corresponds to events where two additional \(\mathrm {b}\) jets are generated in the visible phase space, while \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b}}\) represents the same physical process, where only one additional \(\mathrm {b}\) jet is within the acceptance requirements. In cases with at least two additional jets or two \(\mathrm {b}\) jets, the cross section is also measured as a function of the angular distance between the two jets and their dijet invariant mass. The results are reported both in the visible phase space and extrapolated to the full phase space of the \(\mathrm {t}\overline{\mathrm {t}}\) system to facilitate the comparison with theoretical calculations.

Finally, the fraction of events that do not contain additional jets (gap fraction) is determined as a function of the threshold on the leading and subleading additional-jet \(p_{\mathrm {T}}\) , and the scalar sum of all additional-jet \(p_{\mathrm {T}}\) . This was first measured in Refs. [5, 12].

The results are compared at particle level to theoretical predictions obtained with four different event generators: MadGraph [13], mc@nlo [14], powheg [15], and MG5_aMC@NLO [16], interfaced with either pythia [17] or herwig [18], and in the case of powheg with both. Additionally, the measurements as a function of the \(\mathrm {b}\) jet quantities are compared to the predictions from the event generator PowHel [19].

This paper is structured as follows. A brief description of the CMS detector is provided in Sect. 2. Details of the event simulation generators and their theoretical predictions are given in Sect. 3. The event selection and the method used to identify the additional radiation in the event for both \(\mathrm {t}\overline{\mathrm {t}}\) +jets and \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b} \overline{\mathrm {b}}}\) (\(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b}}\)) studies are presented in Sects. 4 and 5. The cross section measurement and the systematic uncertainties are described in Sects. 6 and 7. The results as a function of the jet multiplicity and the kinematic properties of the additional jets and \(\mathrm {b}\) jets are presented in Sects. 810. The definition of the gap fraction and the results are described in Sect. 11. Finally, a summary is given in Sect. 12.

2 The CMS detector

The central feature of the CMS apparatus is a superconducting solenoid of 6 m internal diameter, providing a magnetic field of 3.8 T. Within the solenoid volume are a silicon pixel and strip tracker, a lead tungstate crystal electromagnetic calorimeter (ECAL), and a brass and scintillator hadron calorimeter, each composed of a barrel and two endcap sections. Extensive forward calorimetry complements the coverage provided by the barrel and endcap detectors. Muons are measured in gas-ionization detectors embedded in the steel flux-return yoke outside the solenoid. A more detailed description of the CMS detector, together with a definition of the coordinate system used and the relevant kinematic variables, can be found in Ref. [10].

3 Event simulation and theoretical predictions

Experimental effects coming from event reconstruction, selection criteria, and detector resolution are modelled using Monte Carlo (MC) event generators interfaced with a detailed simulation of the CMS detector response using Geant4 (v. 9.4) [20].

The MadGraph (v. 5.1.5.11) [13] generator calculates the matrix elements at tree level up to a given order in \(\alpha _s\). In particular, the simulated \(\mathrm {t}\overline{\mathrm {t}}\) sample used in this analysis is generated with up to three additional partons. The MadSpin [21] package is used to incorporate spin correlations of the top quark decay products. The value of the top quark mass is chosen to be \(m_{\mathrm {t}} = 172.5\,\text {GeV} \), and the proton structure is described by the CTEQ6L1 [22] set of parton distribution functions (PDF). The generated events are subsequently processed with pythia (v. 6.426) [17] for fragmentation and hadronization, using the MLM prescription for the matching of higher-multiplicity matrix element calculations with parton showers [23]. The pythia parameters for the underlying event, parton shower, and hadronization are set according to the Z2* tune, which is derived from the Z1 tune [24]. The Z1 tune uses the CTEQ5L PDFs, whereas Z2* adopts CTEQ6L.

In addition to the nominal \(\mathrm {t}\overline{\mathrm {t}}\) MadGraph sample, dedicated samples are generated by varying the central value of the renormalization (\(\mu _\mathrm {R}\)) and factorization (\(\mu _\mathrm {F}\)) scales and the matrix element/parton showering matching scale (jet-parton matching scale). These samples are produced to determine the systematic uncertainties in the measurement owing to the theoretical assumptions on the modelling of \(\mathrm {t}\overline{\mathrm {t}}\) events, as well as for comparisons with the measured distributions. The nominal values of \(\mu _\mathrm {R}\) and \(\mu _\mathrm {F}\) are defined by the \(Q^2\) scale in the event: \(\mu _\mathrm {R}^2 =\mu _\mathrm {F}^2 = Q^2 = m_{\mathrm {t}}^2 + \sum {p_{\mathrm {T}} ^2(\text {jet})}\), where the sum runs over all the additional jets in the event not coming from the \(\mathrm {t}\overline{\mathrm {t}}\) decay. The samples with the varied scales use \(\mu _\mathrm {R}^2 =\mu _\mathrm {F}^2 = 4Q^2\) and \(Q^2/4\), respectively. For the nominal MadGraph sample, a jet-parton matching scale of 40\(\,\text {GeV}\) is chosen, while for the varied samples, values of 60 and 30\(\,\text {GeV}\) are employed, respectively. These scales correspond to jet-parton matching thresholds of 20\(\,\text {GeV}\) for the nominal sample, and 40 and 10\(\,\text {GeV}\) for the varied ones.

The powheg (v. 1.0 r1380) and mc@nlo (v. 3.41) generators, along with the CT10 [25] and CTEQ6M [22] PDFs, are used, respectively, for comparisons with the data. The powheg generator simulates calculations of \(\mathrm {t}\overline{\mathrm {t}}\) production to full NLO accuracy, and is matched with two parton shower MC generators: the pythia (v. 6.426) Z2* tune (designated as pythia 6 in the following), and the herwig [18] (v. 6.520) AUET2 tune [26] (referred to as herwig 6 in the following). The parton showering in pythia is based on a transverse-momentum ordering of parton showers, whereas herwig uses angular ordering. The mc@nlo generator implements the hard matrix element to full NLO accuracy, matched with herwig (v. 6.520) for the initial- and final-state parton showers using the default tune. These two generators, powheg and mc@nlo, are formally equivalent up to the NLO accuracy, but they differ in the techniques used to avoid double counting of radiative corrections that may arise from interfacing with the parton showering generators.

The cross section as a function of jet multiplicity and the gap fraction measurements are compared to the NLO predictions of the powheg (v2) [15] and MG5_aMC@NLO [16] generators. The powheg (v2) generator is matched to the pythia (v. 8.205) CUETP8M1 tune [27] (referred to as pythia 8), herwig 6, and pythia 6. In these samples the hdamp parameter of powhegbox, which controls the matrix element and parton shower matching and effectively regulates the high-\(p_{\mathrm {T}}\) radiation, is set to \(m_{\mathrm {t}}= 172.5\,\text {GeV} \). The MG5_aMC@NLO generator simulates \(\mathrm {t}\overline{\mathrm {t}}\) events with up to two additional partons at NLO, and is matched to the pythia 8 parton shower simulation using the FxFx merging prescription [28]. The top quark mass value used in all these simulations is also 172.5\(\,\text {GeV}\) and the PDF set is NNPDF3.0 [29]. In addition, a \(\mathrm {t}\overline{\mathrm {t}}\) MadGraph sample matched to pythia 8 for the parton showering and hadronization is used for comparisons with the data.

The \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b} \overline{\mathrm {b}}}\) production cross sections are also compared with the predictions by the generator PowHel [19] (HELAC-NLO [30] + powhegbox [31]), which implements the full \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b} \overline{\mathrm {b}}}\) process at NLO QCD accuracy, with parton shower matching based on the powheg NLO matching algorithm [15, 32]. The events are further hadronized by means of pythia (v. 6.428), using parameters of the Perugia 2011 C tune [33]. In the generation of the events, the renormalization and factorization scales are fixed to \(\mu _\mathrm {R} = \mu _\mathrm {F} = H_{\mathrm {T}}/4\), where \(H_{\mathrm {T}} \) is the sum of the transverse energies of the final-state partons (\(\mathrm {t}\), \(\overline{\mathrm {t}}\), \(\mathrm {b}\), \(\overline{\mathrm {b}}\)) from the underlying tree-level process, and the CT10 PDFs are used.

The SM background samples are simulated with MadGraph, powheg, or pythia, depending on the process. The MadGraph generator is used to simulate Z/\(\gamma*\) production (referred to as Drell–Yan, DY, in the following), \(\mathrm {t}\overline{\mathrm {t}}\) production in association with an additional boson (referred to as \(\mathrm {t}\overline{\mathrm {t}} \)+\(\mathrm {Z} \), \(\mathrm {t}\overline{\mathrm {t}} \)+\(\mathrm {W}\), and \(\mathrm {t}\overline{\mathrm {t}} \)+\(\mathrm {\gamma }\)), and \(\mathrm {W}\) boson production with additional jets (\(\mathrm {W}\)+jets in the following). Single top quark events (\(\mathrm {t} \mathrm {W}\) channel) are simulated using powheg. Diboson (\(\mathrm {W}\mathrm {W}\), \(\mathrm {W}\mathrm {Z} \), and \(\mathrm {Z} \mathrm {Z} \)) and QCD multijet events are simulated using pythia. For the \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b}}\) and \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b} \overline{\mathrm {b}}}\) measurements, the expected contribution from SM \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {H}}\) processes, simulated with pythia, is also considered, although the final state has not yet been observed.

For comparison with the measured distributions, the events in the simulated samples are normalized to an integrated luminosity of \(19.7{\,\text {fb}^\text {-1}} \) according to their predicted cross sections. These are taken from next-to-next-to-leading-order (NNLO) (\(\mathrm {W}\)+jets [34] and DY [35]), NLO + next-to-next-to-leading logarithmic (NNLL) (single top quark \(\mathrm {t} \mathrm {W}\) channel [36]), NLO (diboson [37], \(\mathrm {t}\overline{\mathrm {t}} \)+\(\mathrm {Z} \) [38], \(\mathrm {t}\overline{\mathrm {t}} \)+\(\mathrm {W}\) [38], and \(\mathrm {t}\overline{\mathrm {t}} \)+\(\mathrm {H} \) [39]), and leading-order (LO) (QCD multijet [17]) calculations. The contribution of QCD multijet events is found to be negligible. The predicted cross section for the \(\mathrm {t}\overline{\mathrm {t}} \)+\(\mathrm {\gamma }\) sample is obtained by scaling the LO cross section obtained with the Whizard event generator [40] by an NLO/LO K-factor correction [41]. The \(\mathrm {t}\overline{\mathrm {t}}\) simulated sample is normalized to the total cross section \(\sigma _{\mathrm {t}\overline{\mathrm {t}}} = 252.9\, \pm \, ^{6.4}_{8.6} \text {(scale)} \pm 11.7 (\mathrm {PDF}+\alpha _s)\,{\mathrm{pb}} \), calculated with the Top++2.0 program to NNLO in perturbative QCD, including soft-gluon resummation to NNLL order [42], and assuming \(m_{\mathrm {t}} = 172.5\,\text {GeV} \). The first uncertainty comes from the independent variation of the factorization and renormalization scales, \(\mu _\mathrm {R}\) and \(\mu _\mathrm {F}\), while the second one is associated with variations in the PDF and \(\alpha _s\), following the PDF4LHC prescription with the MSTW2008 68 % confidence level (CL) NNLO, CT10 NNLO, and NNPDF2.3 5f FFN PDF sets (see Refs. [43, 44] and references therein and Refs. [4547]).

A number of additional pp simulated hadronic interactions (“pileup”) are added to each simulated event to reproduce the multiple interactions in each bunch crossing from the luminosity conditions in the real data taking. Correction factors for detector effects (described in Sects. 4 and 6) are applied, when needed, to improve the description of the data by the simulation.

4 Event reconstruction and selection

The event selection is based on the decay topology of the \(\mathrm {t}\overline{\mathrm {t}}\) events, where each top quark decays into a \(\mathrm {W}\) boson and a \(\mathrm {b}\) quark. Only the cases in which both \(\mathrm {W}\) bosons decayed to a charged lepton and a neutrino are considered. These signatures imply the presence of isolated leptons, missing transverse momentum owing to the neutrinos from \(\mathrm {W}\) boson decays, and highly energetic jets. The heavy-quark content of the jets is identified through \(\mathrm {b}\) tagging techniques. The same requirements are applied to select the events for the different measurements, with the exception of the requirements on the \(\mathrm {b}\) jets, which have been optimized independently for the \(\mathrm {t}\overline{\mathrm {t}}\) +jets and \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b} \overline{\mathrm {b}}}\) (\(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b}}\)) cases. The description of the event reconstruction and selection is detailed in the following.

Events are reconstructed using a particle-flow (PF) algorithm, in which signals from all subdetectors are combined [48, 49]. Charged particles are required to originate from the primary collision vertex [50], defined as the vertex with the highest sum of \(p_{\mathrm {T}} ^2\) of all reconstructed tracks associated with it. Therefore, charged-hadron candidates from pileup events, i.e. originating from additional pp interactions within the same bunch crossing, are removed before jet clustering on an event-by-event basis. Subsequently, the remaining neutral-particle component from pileup events is accounted for through jet energy corrections [51].

Muon candidates are reconstructed from tracks that can be linked between the silicon tracker and the muon system [52]. The muons are required to have \(p_{\mathrm {T}} >20\,\text {GeV} \), be within \(|\eta |<2.4\), and have a relative isolation \(I_{\text {rel}}<0.15\). The parameter \(I_{\text {rel}}\) is defined as the sum of the \(p_{\mathrm {T}}\) of all neutral and charged reconstructed PF candidates, except the muon itself, inside a cone of \(\varDelta R\equiv \sqrt{(\varDelta \eta )^2+(\varDelta \phi )^2} < 0.3\) around the muon direction, divided by the muon \(p_{\mathrm {T}}\), where \(\varDelta \eta \) and \(\varDelta \phi \) are the difference in pseudorapidity and azimuthal angle between the directions of the candidate and the muon, respectively. Electron candidates are identified by combining information from charged-track trajectories and energy deposition measurements in the ECAL [53], and are required to be within \(|\eta |<2.4\), have a transverse energy of at least 20\(\,\text {GeV}\), and fulfill \(I_{\text {rel}} < 0.15\) inside a cone of \(\varDelta R < 0.3\). Electrons from identified photon conversions are rejected. The lepton identification and isolation efficiencies are determined via a tag-and-probe method using \(\mathrm {Z}\) boson events.

Jets are reconstructed by clustering the PF candidates, using the anti-\(k_{\mathrm {T}} \) clustering algorithm [54, 55] with a distance parameter of 0.5. The jet momentum is determined as the vectorial sum of all particle momenta in the jet, and is found in the simulation to be within 5 to 10 % of the true momentum over the entire \(p_{\mathrm {T}}\) range and detector acceptance. Jet energy corrections are derived from the simulation, and are confirmed with in situ measurements with the energy balance of dijet and photon+jet events [56]. The jet energy resolution amounts typically to 15 % at 10\(\,\text {GeV}\) and 8 % at 100\(\,\text {GeV}\). Muons and electrons passing less stringent requirements compared to the ones mentioned above are identified and excluded from the clustering process. Jets are selected in the interval \(|\eta |<2.4\) and with \(p_{\mathrm {T}} >20\,\text {GeV} \). Additionally, the jets identified as part of the decay products of the \(\mathrm {t}\overline{\mathrm {t}}\) system (cf. Sect. 5) must fulfill \(p_{\mathrm {T}} >30\,\text {GeV} \). Jets originating from the hadronization of \(\mathrm {b}\) quarks are identified using a combined secondary vertex algorithm (CSV) [57], which provides a \(\mathrm {b}\) tagging discriminant by combining identified secondary vertices and track-based lifetime information.

The missing transverse energy (\(/\!\!\!\!E_{\mathrm {T}}\)) is defined as the magnitude of the projection on the plane perpendicular to the beams of the negative vector sum of the momenta of all reconstructed particles in an event [58]. To mitigate the effect of contributions from pileup on the \(/\!\!\!\!E_{\mathrm {T}}\) resolution, we use a multivariate correction where the measured momentum is separated into components that originate from the primary and the other collision vertices [59]. This correction improves the \(/\!\!\!\!E_{\mathrm {T}}\) resolution by \({\approx }5~\%\).

Events are triggered by requiring combinations of two leptons (\(\ell \) = e or \(\mu \)), where one fulfills a \(p_{\mathrm {T}}\) threshold of 17\(\,\text {GeV}\) and the other of 8\(\,\text {GeV}\), irrespective of the flavour of the leptons. The dilepton trigger efficiencies are measured using samples selected with triggers that require a minimum \(/\!\!\!\!E_{\mathrm {T}}\) or number of jets in the event, and are only weakly correlated to the dilepton triggers used in the analysis.

Events are selected if there are at least two isolated leptons of opposite charge. Events with a lepton pair invariant mass less than 20\(\,\text {GeV}\) are removed to suppress events from heavy-flavour resonance decays, QCD multijet, and DY production. In the \(\mu \mu \) and ee channels, the dilepton invariant mass is required to be outside a \(\mathrm {Z}\) boson mass window of \(91\pm 15\,\text {GeV} \), and \(/\!\!\!\!E_{\mathrm {T}}\) is required to be larger than 40\(\,\text {GeV}\).

For the \(\mathrm {t}\overline{\mathrm {t}}\) +jets selection, a minimum of two jets is required, of which at least one must be tagged as a \(\mathrm {b}\) jet. A loose CSV discriminator value is chosen such that the efficiency for tagging jets from \(\mathrm {b}\) (\(\mathrm {c}\)) quarks is \({\approx }85~\%\) (40 %), while the probability of tagging jets originating from light quarks (\(\mathrm {u} \), \(\mathrm {d} \), or \(\mathrm {s} \)) or gluons is around 10 %. Efficiency corrections, depending on jet \(p_{\mathrm {T}}\) and \(\eta \), are applied to account for differences in the performance of the \(\mathrm {b}\) tagging algorithm between data and simulation.

For the \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b} \overline{\mathrm {b}}}\) (\(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b}}\)) selection, at least three \(\mathrm {b}\)-tagged jets are required (without further requirements on the minimum number of jets). In this case, a tighter discriminator value [57] is chosen to increase the purity of the sample. The efficiency of this working point is approximately 70 % (20 %) for jets originating from a \(\mathrm {b}\) (\(\mathrm {c}\)) quark, while the misidentification rate for light-quark and gluon jets is around 1 %. The shape of the CSV discriminant distribution in simulation is corrected to better describe the efficiency observed in the data. This correction is derived separately for light-flavour and \(\mathrm {b}\) jets from a tag-and-probe approach using control samples enriched in events with a \(\mathrm {Z}\) boson and exactly two jets, and \(\mathrm {t}\overline{\mathrm {t}}\) events in the \(\mathrm {e}\mu \) channel with no additional jets [60].

5 Identification of additional radiation in the event

To study additional jet activity in the data, the identification of jets arising from the decay of the \(\mathrm {t}\overline{\mathrm {t}}\) system is crucial. In particular, we need to identify correctly the two \(\mathrm {b}\) jets from the top quark decays in events with more than two \(\mathrm {b}\) jets. This is achieved by following two independent but complementary approaches: a kinematic reconstruction [61] and a multivariate analysis, optimized for the two cases under study, \(\mathrm {t}\overline{\mathrm {t}}\) +jets and \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b} \overline{\mathrm {b}}}\) (\(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b}}\)), respectively. The purpose of the kinematic reconstruction is to completely reconstruct the \(\mathrm {t}\overline{\mathrm {t}}\) system based on \(/\!\!\!\!E_{\mathrm {T}}\) and the information on identified jets and leptons, taking into account detector resolution effects. This method is optimized for the case where the \(\mathrm {b}\) jets in the event only arise from the decay of the top quark pair. The multivariate approach is optimized for events with more \(\mathrm {b}\) jets than just those from the \(\mathrm {t}\overline{\mathrm {t}}\) system. This method identifies the two jets that most likely originated from the top quark decays, and the additional \(\mathrm {b}\) jets, but does not perform a full reconstruction of the \(\mathrm {t}\overline{\mathrm {t}}\) system. Both methods are described in the following sections.

5.1 Kinematic reconstruction in \(\mathrm {t}\overline{\mathrm {t}} \)+jets events

The kinematic reconstruction method was developed and used for the first time in the analysis from Ref. [8]. In this method the following constraints are imposed: \(/\!\!\!\!E_{\mathrm {T}}\) is assumed to originate solely from the two neutrinos; the \(\mathrm {W}\) boson invariant mass is fixed to \(80.4\,\text {GeV} \) [62]; and the top quark and antiquark masses are fixed to a value of \(172.5\,\text {GeV} \). Each pair of jets and lepton-jet combination fulfilling the selection criteria is considered in the kinematic reconstruction. Effects of detector resolution are accounted for by randomly smearing the measured energies and directions of the reconstructed lepton and \(\mathrm {b}\) jet candidates by their resolutions. These are determined from the simulation of signal events by comparing the reconstructed \(\mathrm {b}\) jets and leptons matched to the generated \(\mathrm {b}\) quarks and leptons from top quark decays. For a given smearing, the solution of the equations for the neutrino momenta yielding the smallest invariant mass of the \(\mathrm {t}\overline{\mathrm {t}}\) system is chosen. For each solution, a weight is calculated based on the expected invariant mass spectrum of the lepton and \(\mathrm {b}\) jet from the top quark decays at the parton level. The weights are summed over 100 randomly smeared reconstruction attempts, and the kinematics of the top quark and antiquark are calculated as a weighted average. Finally, the two jets and lepton-jet combinations that yield the maximum sum of weights are chosen for further analysis. Combinations with two \(\mathrm {b}\)-tagged jets are chosen over those with a single \(\mathrm {b}\)-tagged jet. The efficiency of the kinematic reconstruction, defined as the number of events with a solution divided by the total number of selected \(\mathrm {t}\overline{\mathrm {t}}\) +jets events, is approximately 94 %. The efficiency in simulation is similar to the one in data for all jet multiplicities. Events with no valid solution for the neutrino momenta are excluded from further analysis. In events with additional jets, the algorithm correctly identifies the two jets coming from the \(\mathrm {t}\overline{\mathrm {t}}\) decay in about 70 % of the cases.

After the full event selection is applied, the dominant background in the \(\mathrm {e}\mu \) channel originates from other \(\mathrm {t}\overline{\mathrm {t}}\) decay channels and is estimated using simulation. This contribution corresponds mostly to leptonic \(\tau \) decays, which are considered background in the \(\mathrm {t}\overline{\mathrm {t}}\) +jets measurements. In the \(\mathrm {e}\mathrm {e}\) and \(\mathrm {\mu }\mathrm {\mu }\) channels, the dominant background contribution arises from \(\mathrm {Z}/ \mathrm {\gamma }^*\)+jets production. The normalization of this background contribution is derived from data using the events rejected by the \(\mathrm {Z}\) boson veto, scaled by the ratio of events failing and passing this selection, estimated from simulation [63]. The remaining backgrounds, including the single top quark \(\mathrm {t} \mathrm {W}\) channel, \(\mathrm {W}\)+jets, diboson, and QCD multijet events, are estimated from simulation for all the channels.

In Fig. 1, the multiplicity distributions of the selected jets per event are shown for different jet \(p_{\mathrm {T}}\) thresholds and compared to SM predictions. In this figure and the following ones, the \(\mathrm {t}\overline{\mathrm {t}}\) sample is simulated using MadGraph +pythia 6, where only \(\mathrm {t}\overline{\mathrm {t}}\) events with two leptons (\(\mathrm {e}\) or \(\mu \)) from the \(\mathrm {W}\) boson decay are considered as signal. All other \(\mathrm {t}\overline{\mathrm {t}}\) events, specifically those originating from decays via \(\tau \) leptons, which are the dominant contribution, are considered as background. In the following figures, “Electroweak” corresponds to DY, \(\mathrm {W}\)+jets, and diboson processes, and “\(\mathrm {t}\overline{\mathrm {t}}\) bkg.” includes the \(\mathrm {t}\overline{\mathrm {t}} \)+\(\mathrm {\gamma }/\mathrm {W}/\mathrm {Z} \) events. The data are well described by the simulation, both for the low jet \(p_{\mathrm {T}}\) threshold of 30\(\,\text {GeV}\) and the higher thresholds of 60 and 100\(\,\text {GeV}\). The hatched regions in Figs. 1, 2 and 3 correspond to the uncertainties affecting the shape of the simulated signal and background events (cf. Sect. 6), and are dominated by modelling uncertainties in the former.

Fig. 1
figure 1

Reconstructed jet multiplicity distribution after event selection in data (points) and from signal and background simulation (histograms) for all jets with \(p_{\mathrm {T}}\) of at least 30\(\,\text {GeV}\) (top), 60\(\,\text {GeV}\) (bottom left), and 100\(\,\text {GeV}\) (bottom right). The hatched regions correspond to the uncertainties affecting the shape of the distributions in the simulated signal \(\mathrm {t}\overline{\mathrm {t}}\) events and backgrounds (cf. Sect. 6). The lower plots show the ratio of the data to the MC simulation prediction. Note that in all cases the event selection requires at least two jets with \(p_{\mathrm {T}} > 30\,\text {GeV} \)

Additional jets in the event are defined as those jets within the phase space described in the event selection (cf. Sect. 4) that are not identified by the kinematic reconstruction to be part of the \(\mathrm {t}\overline{\mathrm {t}}\) system. The \(\eta \) and \(p_{\mathrm {T}}\) distributions of the additional jets with the largest and second largest \(p_{\mathrm {T}}\) in the event (referred to as the leading and subleading additional jets in the following) are shown in Fig. 2. Three additional event variables are considered: the scalar sum of the \(p_{\mathrm {T}}\) of all additional jets, \(H_{\mathrm {T}} \), the invariant mass of the leading and subleading additional jets, \(m_{\mathrm {jj}}\), and their angular separation, \(\varDelta R_{\mathrm {jj}} =\sqrt{{(\varDelta \eta )^2+(\varDelta \phi )^2}}\), where \(\varDelta \eta \) and \(\varDelta \phi \) are the pseudorapidity and azimuthal differences between the directions of the two jets. These distributions are shown in Fig. 3. The predictions from the simulation, also shown in the figures, describe the data within the uncertainties.

Fig. 2
figure 2

Distribution of the \(\eta \) (left) and \(p_{\mathrm {T}}\) (right) of the leading (top row) and subleading (bottom row) additional reconstructed jets in data (points) and from signal and background simulation (histograms). The hatched regions correspond to the uncertainties affecting the shape of the simulated distributions in the signal \(\mathrm {t}\overline{\mathrm {t}}\) events and backgrounds (cf. Sect. 6). The lower plots show the ratio of the data to the MC simulation prediction

Fig. 3
figure 3

Distribution of the scalar sum of the \(p_{\mathrm {T}}\) of all additional jets \(H_{\mathrm {T}} \) (top), the invariant mass of the leading and subleading additional jets \(m_{\mathrm {jj}}\) (bottom left), and their angular distance \(\varDelta R_{\mathrm {jj}}\) (bottom right) in data (points) and from signal and background simulation (histograms). The hatched regions correspond to the uncertainties affecting the shape of the distributions in the simulated signal \(\mathrm {t}\overline{\mathrm {t}}\) events and backgrounds (cf. Sect. 6). The lower plots show the ratio of the data to the MC simulation prediction

5.2 Identification of \(\mathrm {t}\overline{\mathrm {t}}\) jets and additional jets in \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b} \overline{\mathrm {b}}}\) events

The multivariate approach uses a boosted decision tree (BDT) to distinguish the \(\mathrm {b}\) jets stemming from the \(\mathrm {t}\overline{\mathrm {t}}\) system from those arising from additional radiation for final states with more than two \(\mathrm {b}\) jets. This method is optimized for \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b} \overline{\mathrm {b}}}\) topologies in the dilepton final state of the \(\mathrm {t}\overline{\mathrm {t}}\) system. The BDT is set up using the TMVA package [64]. To avoid any dependence on the kinematics of the additional jets, and especially on the invariant mass of the two additional jets, the method identifies the jets stemming from the \(\mathrm {t}\overline{\mathrm {t}}\) system by making use of properties of the \(\mathrm {t}\overline{\mathrm {t}}\) system that are expected to be mostly insensitive to the additional radiation. The variables combine information from the two final-state leptons, the jets, and \(/\!\!\!\!E_{\mathrm {T}}\). All possible pairs of reconstructed jets in an event are considered. For each pair, one jet is assigned to the \(\mathrm {b}\) jet and the other to the \(\overline{\mathrm {b}}\) jet. This assignment is needed to define the variables used in the BDT and is based on the measurement of the charge of each jet, which is calculated from the charge and the momenta of the PF constituents used in the jet clustering. The jet in the pair with the largest charge is assigned to the \(\overline{\mathrm {b}}\), while the other jet is assigned to the \(\mathrm {b}\). The efficiency of this jet charge pairing is defined as the fraction of events where the assigned \(\mathrm {b}\) and \(\overline{\mathrm {b}}\) are correctly matched to the corresponding generated b and \(\overline{\mathrm {b}}\) jets, and amounts to 68 %.

A total of twelve variables are included in the BDT. Some examples of the variables used are: the sum and difference of the invariant mass of the \(\mathrm {b} \ell ^+\) and \(\overline{\mathrm {b}} \ell ^-\) systems, \(m^{\mathrm {b} \ell ^+}\pm m^{\overline{\mathrm {b}} \ell ^-}\); the absolute difference in the azimuthal angle between them, \( |\varDelta \phi ^{ \mathrm {b} \ell ^+,\overline{\mathrm {b}} \ell ^- } |\); the \(p_{\mathrm {T}}\) of the \(\mathrm {b} \ell ^+\) and \(\overline{\mathrm {b}} \ell ^-\) systems, \(p_{\mathrm {T}} ^{\mathrm {b} \ell ^+}\) and \(p_{\mathrm {T}} ^{\overline{\mathrm {b}} \ell ^-}\); and the difference between the invariant mass of the two \(\mathrm {b}\) jets and two leptons and the invariant mass of the \(\mathrm {b} \overline{\mathrm {b}} \) pair, \(m^{\mathrm {b} \overline{\mathrm {b}} \ell ^+\ell ^-}-m^{\mathrm {b} \overline{\mathrm {b}} }\). The complete list of variables can be found in Appendix A. The main challenge with this method is the large number of possible jet assignments, given four genuine \(\mathrm {b}\) jets and potential extra jets from additional radiation in each event. The basic methodology is to use the BDT discriminant value of each dijet combination as a measure of the probability that the combination stems from the \(\mathrm{t}\overline{\mathrm{t}}\) system. The jets from the \(\mathrm {t}\overline{\mathrm {t}}\) system are then identified as the pair with the highest BDT discriminant. From the remaining jets, those \(\mathrm {b}\)-tagged jets with the highest \(p_{\mathrm {T}}\) are selected as being the leading additional ones.

The BDT training is performed on a large and statistically independent sample of simulated \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {H}}\) events with the Higgs boson mass varied over the range 110–140\(\,\text {GeV}\). The \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b} \overline{\mathrm {b}}}\) events are not included in the training to avoid the risk of overtraining owing to the limited number of events in the available simulated samples. The simulated \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {H} \,(\mathrm {b} \overline{\mathrm {b}})}\) sample is suited for this purpose since the four \(\mathrm {b}\) jets from the decay of the \(\mathrm {t}\overline{\mathrm {t}}\) system and the Higgs boson have similar kinematic distributions. Since it is significantly harder to identify the jets from the \(\mathrm {t}\overline{\mathrm {t}}\) system in \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {H}}\) events than in \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b} \overline{\mathrm {b}}}\) events, where the additional \(\mathrm {b}\) jets arise from initial- or final-state radiation, a good BDT performance with \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {H}}\) events implies also a good identification in \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b} \overline{\mathrm {b}}}\) events. The distributions of the BDT discriminant in data and simulation are shown in Fig. 4 for all dijet combinations in an event, and for the combination with the highest weight that is assigned to the \(\mathrm {t}\overline{\mathrm {t}}\) system. The subset “Minor bkg.” includes all non-\(\mathrm {t}\overline{\mathrm {t}}\) processes and \(\mathrm {t}\overline{\mathrm {t}} \)+\(\mathrm {Z}/\mathrm {W}/\mathrm {\gamma }\) events. There is good agreement between the data and simulation distributions within the statistical uncertainties.

Fig. 4
figure 4

The BDT discriminant of all dijet combinations in data (points) and from signal and background simulation (histograms) per event (left) and dijet combination with the highest discriminant per event (right) in events with at least four jets and exactly four \(\mathrm {b}\)-tagged jets. The distributions include the correction obtained with the template fit to the \(\mathrm {b}\)-tagged jet multiplicity (cf. Sect. 5.2). The hatched area represents the statistical uncertainty in the simulated samples. “Minor bkg.” includes all non-\(\mathrm {t}\overline{\mathrm {t}}\) processes and \(\mathrm {t}\overline{\mathrm {t}} \)+\(\mathrm {Z}/\mathrm {W}/\mathrm {\gamma }\). The lower plots show the ratio of the data to the MC simulation prediction

The number of simulated events with correct assignments for the additional \(\mathrm {b}\) jets in \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {H}}\) events relative to the total number of events where those jets are selected and matched to the corresponding generator jets, is approximately 34 %. In \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b} \overline{\mathrm {b}}}\) events, this fraction is about 40 %. This efficiency is high enough to allow the measurement of the \(\mathrm {t}\overline{\mathrm {t}}\) cross section as a function of the kinematic variables of the additional \(\mathrm {b}\) jets (the probability of selecting the correct assignments by choosing random combinations of jets is 17 % in events with four jets and 10 % in events with five jets). The relative increase in efficiency with respect to the use of the kinematic reconstruction for \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b} \overline{\mathrm {b}}}\) is about 15 %. Additionally, the BDT approach improves the correlation between the generated and reconstructed variables, especially for the distribution of the invariant mass of the two leading additional \(\mathrm {b}\) jets \(m_{\mathrm {b} \mathrm {b}}\) and their angular separation \(\varDelta R_{\mathrm {b} \mathrm {b}} = \sqrt{{(\varDelta \eta )^2+(\varDelta \phi )^2}}\), where \(\varDelta \eta \) and \(\varDelta \phi \) are the pseudorapidity and azimuthal differences between the directions of the two \(\mathrm {b}\) jets.

The expected fraction of events with additional \(\mathrm {b}\) jets is not properly modelled in the simulation, in agreement with the observation of a previous CMS measurement [11]. This discrepancy between the MadGraph+pythia simulation and data can be seen in the \(\mathrm {b}\) jet multiplicity distribution, as shown in Fig. 5.

Fig. 5
figure 5

The pre-fit distribution of the \(\mathrm {b}\) jet multiplicity in data (points) and from signal and background simulation (histograms) for events fulfilling the lepton selection criteria, having \({\ge } 2\) jets, \({\ge } 1\) \(\mathrm {b}\)-tagged jet (left), and the post-fit distribution (right). The hatched area represents the statistical uncertainty in the simulated samples. “Minor bkg.” includes all non-\(\mathrm {t}\overline{\mathrm {t}}\) processes and \(\mathrm {t}\overline{\mathrm {t}} \)+\(\mathrm {Z}/\mathrm {W}/\mathrm {\gamma }\). The lower plots show the ratio of the data to the MC simulation prediction

To improve the description of the data by the simulation, a template fit to the \(\mathrm {b}\)-tagged jet multiplicity distribution is performed using three different templates obtained from simulation. One template corresponds to the \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b}}\) and \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b} \overline{\mathrm {b}}}\) processes, defined at the generator level as the events where one or two additional \(\mathrm {b}\) jets are generated within the acceptance requirements, \(p_{\mathrm {T}} >20\,\text {GeV} \) and \(|\eta | <2.4\), (referred to as “\(\mathrm {t}\overline{\mathrm {t}}\) +HF”). The \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b} \overline{\mathrm {b}}}\) and \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b}}\) processes are combined into a single template because they only differ by the kinematic properties of the second additional \(\mathrm {b}\) jet. Details about the definition of the \(\mathrm {b}\) jets and the acceptance are given in Sect. 7. The second template includes the background contribution coming from \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {c} \overline{\mathrm {c}}}\) and \(\mathrm {t}\overline{\mathrm {t}}\) +light-jets events (referred to as “\(\mathrm {t}\overline{\mathrm {t}}\)  other”), where \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {c} \overline{\mathrm {c}}}\) events are defined as those that have at least one \(\mathrm {c}\) jet within the acceptance and no additional \(\mathrm {b}\) jets. This contribution is not large enough to be constrained by data, therefore it is combined with the \(\mathrm {t}\overline{\mathrm {t}}\) +light-jets process in a single template. The third template contains the remaining background processes, including \(\mathrm {\mathrm {t} \overline{\mathrm {t}} 2 \mathrm {b}}\), which corresponds to events with two additional \(\mathrm {b}\) hadrons that are close enough in direction to produce a single \(\mathrm {b}\) jet. This process, produced by collinear \(\mathrm {g} \rightarrow \mathrm {b} \overline{\mathrm {b}} \) splitting, is treated separately owing to the large theoretical uncertainty in its cross section and insufficient statistical precision to constrain it with data. The normalizations of the first two templates are free parameters in the fit. The third is fixed to the corresponding cross section described in Sect. 3, except for the cross section for the \(\mathrm {\mathrm {t} \overline{\mathrm {t}} 2 \mathrm {b}}\) process, which is corrected by a factor of \(1.74{}_{-0.74}^{+0.69}\) [65]. The normalization factors obtained for the template fit correspond to \(1.66 \pm 0.43\) (\(\mathrm {t}\overline{\mathrm {t}}\) +HF) and \(1.00 \pm 0.01\) (\(\mathrm {t}\overline{\mathrm {t}}\)  other). Details about the uncertainties in those factors are presented in Sect. 6.1.1. The improved description of the \(\mathrm {b}\) jet multiplicity can be seen in Fig. 5 (right).

Figure 6 (top) shows the \(p_{\mathrm {T}}\) and \(|\eta |\) distributions of the leading additional \(\mathrm {b}\) jet, measured in events with at least three \(\mathrm {b}\)-tagged jets (using the tighter discriminator value described in Sect. 4), after the full selection and including all corrections. The distributions of the \(p_{\mathrm {T}}\) and \(|\eta |\) of the second additional \(\mathrm {b}\) jet in events with exactly four \(\mathrm {b}\)-tagged jets, \(\varDelta R_{\mathrm {b} \mathrm {b}}\), and \(m_{\mathrm {b} \mathrm {b}}\) are also presented. The dominant contribution arises from the \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b} \overline{\mathrm {b}}}\) process. The \(\mathrm {t}\overline{\mathrm {t}}\) decays into \(\tau \) leptons decaying leptonically are included as signal to increase the number of \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b}}\) and \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b} \overline{\mathrm {b}}}\) events both in data and simulation. It has been checked that the distribution of the variables of relevance for this analysis do not differ between the leptons directly produced from \(\mathrm {W}\) boson decays and the leptons from \(\tau \) decays within the statistical uncertainties in the selected \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b}}\) and \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b} \overline{\mathrm {b}}}\) events. In general, the variables presented are well described by the simulation, after correcting for the heavy-flavour content measured in data, although the simulation tends to predict smaller values of \(\varDelta R_{\mathrm {b} \mathrm {b}}\) than the data. After the full selection, the dominant background contribution arises from dilepton \(\mathrm {t}\overline{\mathrm {t}}\) events with additional light-quark, gluon, and \(\mathrm {c}\) jets, corresponding to about 50 and 20 % of the total expected yields for the \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b}}\) and \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b} \overline{\mathrm {b}}}\) cases, respectively. Smaller background contributions come from single top quark production, \(\mathrm {t}\overline{\mathrm {t}}\) in association with \(\mathrm {W}\)or \(\mathrm {Z}\) bosons, and \(\mathrm {t}\overline{\mathrm {t}}\) events in the lepton+jets decay channels. The contribution from \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {H} \,(\mathrm {b} \overline{\mathrm {b}})}\) is also small, amounting to 0.9 and 3 % of the total expected events for the \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b}}\) and \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b} \overline{\mathrm {b}}}\) distributions. The contribution from background sources other than top quark production processes such as DY, diboson, or QCD multijet is negligible.

Fig. 6
figure 6

Distributions of the leading additional \(\mathrm {b}\) jet \(p_{\mathrm {T}}\) (top left) and \(|\eta |\) (top right), subleading additional \(\mathrm {b}\) jet \(p_{\mathrm {T}}\) (middle left) and \(|\eta |\) (middle right), \(\varDelta R_{\mathrm {b} \mathrm {b}}\) (bottom left), and \(m_{\mathrm {b} \mathrm {b}}\) (bottom right) from data (points) and from signal and background simulation (histograms). The hatched area represents the statistical uncertainty in the simulated samples. “Minor bkg.” includes all non-\(\mathrm {t}\overline{\mathrm {t}}\) processes and \(\mathrm {t}\overline{\mathrm {t}} \)+\(\mathrm {Z}/\mathrm {W}/\mathrm {\gamma }\). The lower plots show the ratio of the data to the MC simulation prediction

6 Systematic uncertainties

Different sources of systematic uncertainties are considered arising from detector effects, as well as theoretical uncertainties. Each systematic uncertainty is determined individually in each bin of the measurement by varying the corresponding efficiency, resolution, or model parameter within its uncertainty, in a similar way as in the CMS previous measurement of the \(\mathrm {t}\overline{\mathrm {t}}\) differential cross sections [8]. For each variation, the measured differential cross section is recalculated and the difference with respect to the nominal result is taken as the systematic uncertainty. The overall uncertainty in the measurement is then derived by adding all contributions in quadrature, assuming the sources of systematic uncertainty to be fully uncorrelated.

6.1 Experimental uncertainties

The experimental sources of systematic uncertainty considered are the jet energy scale (JES), jet energy resolution (JER), background normalization, lepton trigger and identification efficiencies, \(\mathrm {b}\) tagging efficiency, integrated luminosity, pileup modelling, and kinematic reconstruction efficiency.

The experimental uncertainty from the JES is determined by varying the energy scale of the reconstructed jets as a function of their \(p_{\mathrm {T}}\) and \(\eta \) by its uncertainty [56]. The uncertainty from the JER is estimated by varying the simulated JER by its \(\eta \)-dependent uncertainty [56].

The uncertainty from the normalization of the backgrounds that are taken from simulation is determined by varying the cross section used to normalize the sample, see Sect. 3, by \({\pm }30~\%\). This variation takes into account the uncertainty in the predicted cross section and all other sources of systematic uncertainty [5, 8, 66]. In the case of the tW background, the variation of \({\pm }30~\%\) covers the theoretical uncertainty in the absolute rate, including uncertainties owing to the PDFs. The contribution from the DY process, as determined from data, is varied in the normalization by \({\pm }30~\%\) [1, 63].

The trigger and lepton identification efficiencies in simulation are corrected by lepton \(p_{\mathrm {T}} \) and \(\eta \) multiplicative data-to-simulation scale factors. The systematic uncertainties are estimated by varying the factors by their uncertainties, which are in the range 1–2 %.

For the \(\mathrm {t}\overline{\mathrm {t}}\) +jets measurements, the \(\mathrm {b}\) tagging efficiency in simulation is also corrected by scale factors depending on the \(p_{\mathrm {T}}\) and \(\eta \) of the jet. The shape uncertainty in the \(\mathrm {b}\) tagging efficiency is then determined by taking the maximum change in the shape of the \(p_{\mathrm {T}}\) and \(|\eta |\) distributions of the \(\mathrm {b}\) jet, obtained by changing the scale factors. This is achieved by dividing the \(\mathrm {b}\) jet distributions in \(p_{\mathrm {T}}\) and \(|\eta |\) into two bins at the median of the respective distributions. The \(\mathrm {b}\) tagging scale factors for \(\mathrm {b}\) jets in the first bin are scaled up by half the uncertainties quoted in Ref. [57], while those in the second bin are scaled down, and vice versa, so that a maximum variation is assumed and the difference between the scale factors in the two bins reflects the full uncertainty. The changes are made separately in the \(p_{\mathrm {T}}\) and \(|\eta |\) distributions, and independently for heavy-flavour (b and c) and light-flavour (s, u, d, and gluon) jets, assuming that they are all uncorrelated. A normalization uncertainty is obtained by varying the scale factors up and down by half the uncertainties. The total uncertainty is obtained by summing in quadrature the independent variations.

The uncertainty in the integrated luminosity is 2.6 % [67]. The effect of the uncertainty in the level of pileup is estimated by varying the inelastic pp cross section in simulation by \({\pm }5~\%\).

The uncertainty coming from the kinematic reconstruction method is determined from the uncertainty in the correction factor applied to account for the small difference in efficiency between the simulation and data, defined as the ratio between the events with a solution and the total number of selected events.

6.1.1 Specific systematic uncertainties associated with the \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b} \overline{\mathrm {b}}}\) (\(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b}}\)) measurements

In the \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b} \overline{\mathrm {b}}}\) (\(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b}}\)) measurements, an additional uncertainty associated with the template fit to the \(\mathrm {b}\)-tagged jet multiplicity distribution is considered. Since the input templates are known to finite precision, both the statistical and systematic uncertainties in the templates are taken into account. The considered systematic uncertainties that affect the shapes of the templates are those of the JES, the CSV discriminant scale factors following the method described in [60], the cross section of the \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {c} \overline{\mathrm {c}}}\) process, which is varied by \({\pm }50~\%\) [60], and the uncertainty in the \(\mathrm {\mathrm {t} \overline{\mathrm {t}} 2 \mathrm {b}}\) cross section. This is taken as the maximum between the largest uncertainty from the measurement described in Ref. [65] and the difference between the corrected cross section and the prediction by the nominal MadGraph simulation used in this analysis. This results in a variation of the cross section of about \({\pm }40~\%\). This uncertainty is included as a systematic uncertainty in the shape of the background template.

6.2 Model uncertainties

The impact of theoretical assumptions on the measurement is determined by repeating the analysis, replacing the standard MadGraph signal simulation by alternative simulation samples. The uncertainty in the modelling of the hard-production process is assessed by varying the common renormalization and factorization scale in the MadGraph signal samples up and down by a factor of two with respect to its nominal value of the Q in the event (cf. Sect. 3). Furthermore, the effect of additional jet production in MadGraph is studied by varying up and down by a factor of two the threshold between jet production at the matrix element level and via parton showering. The uncertainties from ambiguities in modelling colour reconnection (CR) effects are estimated by comparing simulations of an underlying-event (UE) tune including colour reconnection to a tune without it (Perugia 2011 and Perugia 2011 noCR tunes, described in Ref. [33]). The modelling of the UE is evaluated by comparing two different Perugia 11 (P11) pythia tunes, mpiHi and TeV, to the standard P11 tune. The dependency of the measurement on the top quark mass is obtained using dedicated samples in which the mass is varied by \({\pm }1\,\text {GeV} \) with respect to the default value used in the simulation. The uncertainty from parton shower modelling is determined by comparing two samples simulated with powheg and mc@nlo, using either pythia or herwig for the simulation of the parton shower, underlying event, and hadronization. The effect of the uncertainty in the PDFs on the measurement is assessed by reweighting the sample of simulated \(\mathrm {t}\overline{\mathrm {t}}\) signal events according to the 52 CT10 error PDF sets, at the 90 % CL [25].

Since the total uncertainty in the \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b}}\) and \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b} \overline{\mathrm {b}}}\) production cross sections is largely dominated by the statistical uncertainty in the data, a simpler approach than for the \(\mathrm {t}\overline{\mathrm {t}}\) +jets measurements is chosen to conservatively estimate the systematic uncertainties: instead of repeating the measurement, the uncertainty from each source is taken as the difference between the nominal MadGraph +pythia sample and the dedicated simulated sample at generator level. In the case of the uncertainty coming from the renormalization and factorization scales, the uncertainty estimated in the previous inclusive cross section measurement [11] is assigned.

6.3 Summary of the typical systematic uncertainties

Typical values of the systematic uncertainties in the absolute differential cross sections are summarized in Table 1 for illustrative purposes. They are the median values of the distribution of uncertainties over all bins of the measured variables. Details on the impact of the different uncertainties in the results are given in Sects. 811.

In general, for the \(\mathrm {t}\overline{\mathrm {t}}\) +jets case, the dominant systematic uncertainties arise from the uncertainty in the JES, as well as from model uncertainties such as the renormalization, factorization, and jet-parton matching scales and the hadronization uncertainties. For the \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b}}\) and \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b} \overline{\mathrm {b}}}\) cross sections, the total uncertainty, including all systematic uncertainties, is only about 10 % larger than the statistical uncertainty. The experimental uncertainties with an impact on the normalization of the expected number of signal events, such as lepton and trigger efficiencies, have a negligible effect on the final cross section determination, since the normalization of the different processes is effectively constrained by the template fit.

Table 1 Summary of the typical systematic uncertainties in the measurements of the \(\mathrm {t}\overline{\mathrm {t}}\) +jets and \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b} \overline{\mathrm {b}}}\) (\(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b}}\)) absolute differential cross sections and their sources. The median of the distribution of uncertainties over all bins of each measured differential cross section is quoted

7 Differential \(\mathrm {t}\overline{\mathrm {t}}\) cross section

The absolute differential \(\mathrm {t}\overline{\mathrm {t}}\) cross section is defined as:

$$\begin{aligned} \frac{\mathrm {d}\sigma _{\mathrm {t}\overline{\mathrm {t}}}}{\mathrm {d}x_i}=\frac{\sum _j A_{ij}^{-1} (N^j_{\text {data}}-N^j_{\text {bkg}})}{\varDelta _x^i \mathcal {L}}, \end{aligned}$$
(1)

where j represents the bin index of the reconstructed variable x, i is the index of the corresponding generator-level bin, \(N^j_{\text {data}}\) is the number of data events in bin j, \(N^j_{\text {bkg}}\) is the number of estimated background events, \(\mathcal {L}\) is the integrated luminosity, and \(\varDelta _x^i\) is the bin width. Effects from detector efficiency and resolution in each bin i of the measurement are corrected by the use of a regularized inversion of the response matrix (symbolized by \(A_{ij}^{-1}\)) described in this section.

For the measurements of \(\mathrm {t}\overline{\mathrm {t}}\) +jets, the estimated number of background events from processes other than \(\mathrm {t}\overline{\mathrm {t}}\) production (\(N_{{{\rm non} {\rm t}\overline{\mathrm {t}} {\rm bkg}}}\)) is subtracted from the number of events in data (N). The contribution from other \(\mathrm {t}\overline{\mathrm {t}}\) decay modes is taken into account by correcting the difference N\(N_{{{\rm non} \mathrm {t}\overline{\mathrm {t}} {\rm bkg}}} \) by the signal fraction, defined as the ratio of the number of selected \(\mathrm {t}\overline{\mathrm {t}}\) signal events to the total number of selected \(\mathrm {t}\overline{\mathrm {t}}\) events, as determined from simulation. This avoids the dependence on the inclusive \(\mathrm {t}\overline{\mathrm {t}}\) cross section used for normalization. For the \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b}}\) and \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b} \overline{\mathrm {b}}}\) production cross sections, where the different \(\mathrm {t}\overline{\mathrm {t}}\) contributions are fitted to the data, the expected contribution from all background sources is directly subtracted from the number of data events.

The normalized differential cross section is derived by dividing the absolute result, Eq. (1), by the total cross section, obtained by integrating over all bins for each observable. Because of the normalization, the systematic uncertainties that are correlated across all bins of the measurement, e.g. the uncertainty in the integrated luminosity, cancel out.

Effects from the trigger and reconstruction efficiencies and resolutions, leading to migrations of events across bin boundaries and statistical correlations among neighbouring bins, are corrected using a regularized unfolding method [8, 68, 69]. The response matrix \(A_{ij}\) that corrects for migrations and efficiencies is calculated from simulated \(\mathrm {t}\overline{\mathrm {t}}\) events using MadGraph. The generalized inverse of the response matrix is used to obtain the unfolded distribution from the measured distribution by applying a \(\chi ^2\) technique. To avoid nonphysical fluctuations, a smoothing prescription (regularization) is applied. The regularization level is determined individually for each distribution using the averaged global correlation method [70]. To keep the bin-to-bin migrations small, the width of bins in the measurements are chosen according to their purity and stability. The purity is the number of events generated and correctly reconstructed in a certain bin divided by the total number of reconstructed events in the same bin. The stability is the ratio of the number of events generated and reconstructed in a bin to the total number of events generated in that bin. The purity and stability of the bins are typically larger than 40–50 %, which ensures that the bin-to-bin migrations are small enough to perform the measurement. The performance of the unfolding procedure is tested for possible biases from the choice of the input model (the \(\mathrm {t}\overline{\mathrm {t}}\) MadGraph simulation). It has been verified that by reweighting the \(\mathrm {t}\overline{\mathrm {t}}\) simulation the unfolding procedure based on the nominal response matrix reproduces the altered shapes within the statistical uncertainties. In addition, \(\mathrm {t}\overline{\mathrm {t}}\) samples simulated with powheg and mc@nlo are employed to obtain the response matrices used in the unfolding for the determination of systematic uncertainties of the model (Sect. 6.2). Therefore, possible effects from the unfolding procedure are already taken into account in the systematic uncertainties.

The differential cross section is reported at the particle level, where objects are defined as follows. Leptons from \(\mathrm {W}\) boson decays are defined after final-state radiation, and jets are defined at the particle level by applying the anti-\(k_{\mathrm {T}} \) clustering algorithm with a distance parameter of 0.5 [54] to all stable particles, excluding the decay products from \(\mathrm {W}\) boson decays into \(\mathrm {e}\nu \), \(\mu \nu \), and leptonic \(\tau \) final states. A jet is defined as a \(\mathrm {b}\) jet if it has at least one \(\mathrm {b}\) hadron associated with it. To perform the matching between \(\mathrm {b}\) hadrons and jets, the \(\mathrm {b}\) hadron momentum is scaled down to a negligible value and included in the jet clustering (so-called ghost matching [51]). The \(\mathrm {b}\) jets from the \(\mathrm {t}\overline{\mathrm {t}}\) decay are identified by matching the \(\mathrm {b}\) hadrons to the corresponding original \(\mathrm {b}\) quarks. The measurements are presented for two different phase-space regions, defined by the kinematic and geometric attributes of the \(\mathrm {t}\overline{\mathrm {t}}\) decay products and the additional jets. The visible phase space is defined by the following kinematic requirements:

  • Leptons: \(p_{\mathrm {T}} >20\,\text {GeV} \), \(|\eta |<2.4\),

  • \(\mathrm {b}\) jets arising from top quarks: \(p_{\mathrm {T}} >30\,\text {GeV} \), \(|\eta |<2.4\),

  • Additional jets and \(\mathrm {b}\) jets: \(p_{\mathrm {T}} >20\,\text {GeV} \), \(|\eta |<2.4\).

The full phase space is defined by requiring only the additional jets or \(\mathrm {b}\) jets be within the above-mentioned kinematic range, without additional requirements on the decay products of the \(\mathrm {t}\overline{\mathrm {t}}\) system, and including the correction for the corresponding dileptonic branching fraction, calculated using the leptonic branching fraction of the \(\mathrm {W}\) boson [62].

In the following sections, the \(\mathrm {t}\overline{\mathrm {t}}\) differential cross section measured as a function of the jet multiplicity in the visible phase space and the results as a function of the kinematic variables of the additional jets in the event, measured in the visible and the full phase-space regions, are discussed. The absolute cross sections are presented as figures and compared to different predictions. The full results are given in tables in Appendix B, along with the normalized differential cross sections measurements.

8 Differential \(\mathrm {t}\overline{\mathrm {t}}\) cross sections as a function of jet multiplicity

In Fig. 7, the absolute differential \(\mathrm {t}\overline{\mathrm {t}}\) cross section is shown for three different jet \(p_{\mathrm {T}}\) thresholds: \(p_{\mathrm {T}} >30\), 60, and \(100\,\text {GeV} \). The results are presented for a nominal top quark mass of \(172.5\,\text {GeV} \). The lower part of each figure shows the ratio of the predictions from simulation to the data. The light and dark bands in the ratio indicate the statistical and total uncertainties in the data for each bin, which reflect the uncertainties for a ratio of 1.0. All predictions are normalized to the measured cross section in the range shown in the histogram, which is evaluated by integrating over all bins for each observable. The results are summarized in Table 2, together with the normalized cross sections. In general, the MadGraph generator interfaced with pythia 6, and powheg interfaced both with herwig 6 and pythia 6, provide reasonable descriptions of the data. The mc@nlo generator interfaced with herwig 6 does not generate sufficiently large jet multiplicities, especially for the lowest jet \(p_{\mathrm {T}}\) threshold. The sensitivity of MadGraph to scale variations is investigated through the comparison of different renormalization, factorization, and jet-parton matching scales with respect to the nominal MadGraph simulation. Variations in the jet-parton matching threshold do not yield large effects in the cross section, while the shape and normalization are more affected by the variations in the renormalization and factorization scales, which lead to a slightly worse description of the data up to high jet multiplicities, compared to their nominal values.

Fig. 7
figure 7

Absolute differential \(\mathrm {t}\overline{\mathrm {t}}\) cross sections as a function of jet multiplicity for jets with \(p_{\mathrm {T}} >30\,\text {GeV} \) (top row), 60\(\,\text {GeV}\) (middle row), and 100\(\,\text {GeV}\) (bottom row). In the figures on the left, the data are compared with predictions from MadGraph interfaced with pythia 6, mc@nlo interfaced with herwig 6, and powheg with pythia 6 and herwig 6. The figures on the right show the behaviour of the MadGraph generator with varied renormalization, factorization, and jet-parton matching scales. The inner (outer) vertical bars indicate the statistical (total) uncertainties. The lower part of each plot shows the ratio of the predictions to the data

In Fig. 8, the results are compared to the predictions from MadGraph and MG5_aMC@NLO interfaced with pythia 8, and the powheg generator with the hdamp parameter set to \(m_{\mathrm {t}}=172.5\,\text {GeV} \) (labelled powheg (h\(_{\text {damp}}=m_{\mathrm {t}}\)) in the legend), interfaced with pythia 6, pythia 8, and herwig 6. The MadGraph and MG5_aMC@NLO simulations interfaced with pythia 8 predict larger jet multiplicities than measured in the data for all the considered \(p_{\mathrm {T}}\) thresholds. In general, no large deviations between data and the different powheg predictions are observed.

Fig. 8
figure 8

Absolute differential \(\mathrm {t}\overline{\mathrm {t}}\) cross sections as a function of jet multiplicity for jets with \(p_{\mathrm {T}} >30\,\text {GeV} \) (top row), 60\(\,\text {GeV}\) (middle row), and 100\(\,\text {GeV}\) (bottom row). In the figures on the left, the data are compared with predictions from MadGraph interfaced with pythia 6 and pythia 8, and MG5_aMC@NLO interfaced with pythia 8. The figures on the right show the behaviour of the powheg generator without and with hdamp set to \(m_{\mathrm {t}}\), matched with different versions and tunes of pythia and herwig 6. The inner (outer) vertical bars indicate the statistical (total) uncertainties. The lower part of each plot shows the ratio of the predictions to the data

The total systematic uncertainty in the absolute differential cross section ranges between 6 to 30 %, while for the normalized cross section it varies from 2 % up to 20 % for the bins corresponding to the highest number of jets. In both cases, the dominant experimental systematic uncertainty arises from the JES, having a maximum value of 16 % for the absolute cross section bin with at least six jets and \(p_{\mathrm {T}} > 30\,\text {GeV} \). Typical systematic uncertainty values range between 0.5 and 8 %, while the uncertainty in the normalized cross section is 0.5–4 %. Regarding the modelling uncertainties, the most relevant ones are the uncertainty in the renormalization and factorization scales and the parton shower modelling, up to 6 and 10 %, respectively. The uncertainties from the assumed top quark mass used in the simulation and the jet-parton matching threshold amount to 1–2 %. Other modelling uncertainties such as PDF, CR, and UE have slightly smaller impact. These uncertainties cancel to a large extent in the normalized results, with typical contributions below 0.5 %. The total contribution from the integrated luminosity, lepton identification, and trigger efficiency, which only affect the normalization, is 3.5 %. This contribution is below 0.1 % for every bin in the normalized results. The uncertainty from the estimate of the background contribution is around 2 % for the absolute cross sections and typically below 0.5 % for the normalized results.

9 Differential \(\mathrm {t}\overline{\mathrm {t}}\) cross sections as a function of the kinematic variables of the additional jets

The absolute and normalized differential cross sections are measured as a function of the kinematic variables of the additional jets in the visible phase space defined in Sect. 7. The results are compared to predictions from four different generators: powheg interfaced with pythia 6 and herwig 6, mc@nlo +herwig 6, and MadGraph +pythia 6 with varied renormalization, factorization, and jet-parton matching scales. All predictions are normalized to the measured cross section over the range of the observable shown in the histogram in the corresponding figures.

The absolute differential cross sections as a function of the \(p_{\mathrm {T}}\) of the leading and subleading additional jets and \(H_{\mathrm {T}} \), the scalar sum of the \(p_{\mathrm {T}}\) of all additional jets in the event, are shown in Fig. 9. The total uncertainties in the absolute cross sections range from 8–14 % for the leading additional jet \(p_{\mathrm {T}}\) and \(H_{\mathrm {T}} \), and up to 40 % for the subleading additional jet \(p_{\mathrm {T}}\), while the systematic uncertainties in the normalized cross sections for the bins with the larger number of events are about 3–4 %. The dominant sources of systematic uncertainties arise in both cases from model uncertainties, in particular the renormalization and factorization scales, and the parton shower modelling (up to 10 % for the absolute cross sections), and JES (3–6 % for the absolute cross sections). The typical contribution of other uncertainties such as the assumed top quark mass in the simulation, background contribution, etc., amounts to 1–3 % and 0.5–1.5 %, for the absolute and normalized cross sections, respectively.

In general, the simulation predictions describe the behaviour of the data for the leading additional jet momenta and \(H_{\mathrm {T}} \), although some predictions, in particular powheg, favour a harder \(p_{\mathrm {T}}\) spectrum for the leading jet. The mc@nlo +herwig 6 prediction yields the largest discrepancies. The varied MadGraph samples provide similar descriptions of the shape of the data, except for MadGraph with the lower \(\mu _\mathrm {R} = \mu _\mathrm {F}\) scale, which worsens the agreement.

The results as a function of \(|\eta |\) are presented in Fig. 10. The typical total systematic uncertainties in the absolute cross sections vary from 6.5–19 % for the leading additional jet and about 11–20 % for the subleading one. The uncertainty in the normalized cross section ranges from 1.5–9 % and 5–14 %, respectively. The shape of the \(|\eta |\) distribution is well modelled by mc@nlo +herwig 6. The distributions from MadGraph and powheg yield a similar description of the data, being slightly more central than mc@nlo . Variations of the MadGraph parameters have little impact on these distributions.

The differential cross section is also measured as a function of the dijet angular separation \(\varDelta R_{\mathrm {jj}}\) and invariant mass \(m_{\mathrm {jj}}\) for the leading and subleading additional jets (Fig. 11). In general, all simulations provide a reasonable description of the distributions for both variables. All results are reported in Tables 3, 4 and 5 in Appendix B. Representative examples of the migration matrices are presented in Fig. 24 in Appendix C.

Fig. 9
figure 9

Absolute differential \(\mathrm {t}\overline{\mathrm {t}}\) cross section as a function of \(p_{\mathrm {T}}\) of the leading additional jet (top) and the subleading additional jet (middle), and \(H_{\mathrm {T}} \) (bottom) in the visible phase space of the \(\mathrm {t}\overline{\mathrm {t}}\) system and the additional jets. Data are compared to predictions from MadGraph +pythia 6, powheg +pythia 6, powheg +herwig 6, and mc@nlo +herwig 6 (left) and to MadGraph with varied renormalization, factorization, and jet-parton matching scales (right). The inner (outer) vertical bars indicate the statistical (total) uncertainties. The lower part of each plot shows the ratio of the predictions to the data

Fig. 10
figure 10

Absolute differential \(\mathrm {t}\overline{\mathrm {t}}\) cross section as a function of the \(|\eta |\) of the leading additional jet (top) and the subleading additional jet (bottom) in the visible phase space of the \(\mathrm {t}\overline{\mathrm {t}}\) system and the additional jets. Data are compared to predictions from MadGraph +pythia 6, powheg +pythia 6, powheg +herwig 6, and mc@nlo +herwig 6 (left) and to MadGraph with varied renormalization, factorization, and jet-parton matching scales (right). The inner (outer) vertical bars indicate the statistical (total) uncertainties. The lower part of each plot shows the ratio of the predictions to the data

Fig. 11
figure 11

Absolute differential \(\mathrm {t}\overline{\mathrm {t}}\) cross section as a function of \(\varDelta R_{\mathrm {jj}}\) between the leading and subleading additional jets (top) and their invariant mass, \(m_{\mathrm {jj}}\) (bottom). Data are compared to predictions from MadGraph +pythia 6, powheg +pythia 6, powheg +herwig 6, and mc@nlo +herwig 6 (left) and to MadGraph with varied renormalization, factorization, and jet-parton matching scales (right). The inner (outer) vertical bars indicate the statistical (total) uncertainties. The lower part of each plot shows the ratio of the predictions to the data

The absolute and normalized differential cross sections are also measured as a function of the kinematic variables of the additional jets and \(\mathrm {b}\) jets in the event for the full phase space of the \(\mathrm {t}\overline{\mathrm {t}}\) system to facilitate comparison with theoretical calculations. In this case, the phase space is defined only by the kinematic requirements on the additional jets.

Figures 12 and 13 show the absolute cross sections as a function of the \(p_{\mathrm {T}}\) and \(|\eta |\) of the leading and subleading additional jets and \(H_{\mathrm {T}} \), while the results as a function of \(\varDelta R_{\mathrm {jj}}\) and \(m_{\mathrm {jj}}\) are presented in Fig. 14.

The total uncertainties range between 8–12 % for the leading jet \(p_{\mathrm {T}}\) and \(H_{\mathrm {T}} \), 10 % at lower \(p_{\mathrm {T}}\) and 40 % in the tails of distribution of the subleading jet \(p_{\mathrm {T}}\). The uncertainties for \(|\eta |\) are 6–16 % and 10–30 % for the leading and subleading additional jets, respectively. The typical uncertainties in the cross section as a function of \(\varDelta R_{\mathrm {jj}}\) and \(m_{\mathrm {jj}}\) are on the order of 10–20 %. The uncertainties are dominated by the JES, scale uncertainties, and shower modelling.

The numerical values are given in Tables 6, 7 and 8 of Appendix B, together with the normalized results. In the latter, the uncertainties are on average 2–3 times smaller than for the absolute cross sections, owing to the cancellation of uncertainties such as the integrated luminosity, lepton identification, and trigger efficiency, as well as a large fraction of the JES and model uncertainties, as discussed in Sect. 8. The dominant systematic uncertainties are still the model uncertainties, although they are typically smaller than for the absolute cross sections.

The shapes of the distributions measured in the full and visible phase-space regions of the \(\mathrm {t}\overline{\mathrm {t}}\) system are similar, while the absolute differential cross sections are a factor of 2.2 larger than those in the visible phase space of the \(\mathrm {t}\overline{\mathrm {t}}\) system (excluding the factor due to the leptonic branching fraction correction \((4.54 \pm 0.10)~\%\) [62]).

Fig. 12
figure 12

Absolute differential \(\mathrm {t}\overline{\mathrm {t}}\) cross section as a function of \(p_{\mathrm {T}}\) of the leading additional jet (top) and the subleading additional jet (middle) and \(H_{\mathrm {T}} \) (bottom) measured in the full phase space of the \(\mathrm {t}\overline{\mathrm {t}}\) system, corrected for acceptance and branching fractions. Data are compared to predictions from MadGraph +pythia 6, powheg +pythia 6, powheg +herwig 6, and mc@nlo +herwig 6 (left) and to MadGraph with varied renormalization, factorization, and jet-parton matching scales (right). The inner (outer) vertical bars indicate the statistical (total) uncertainties. The lower part of each plot shows the ratio of the predictions to the data

Fig. 13
figure 13

Absolute differential \(\mathrm {t}\overline{\mathrm {t}}\) cross section as a function of the \(|\eta |\) of the leading additional jet (top) and the subleading additional jet (bottom) measured in the full phase space of the \(\mathrm {t}\overline{\mathrm {t}}\) system, corrected for acceptance and branching fractions. Data are compared to predictions from MadGraph +pythia 6, powheg +pythia 6, powheg +herwig 6, and mc@nlo +herwig 6 (left) and to MadGraph with varied renormalization, factorization, and jet-parton matching scales (right). The inner (outer) vertical bars indicate the statistical (total) uncertainties. The lower part of each plot shows the ratio of the predictions to the data

Fig. 14
figure 14

Absolute differential \(\mathrm {t}\overline{\mathrm {t}}\) cross section as a function of \(\varDelta R_{\mathrm {jj}}\) between the leading and subleading additional jets (top) and their invariant mass, \(m_{\mathrm {jj}}\) (bottom) measured in the full phase space of the \(\mathrm {t}\overline{\mathrm {t}}\) system, corrected for acceptance and branching fractions. Data are compared to predictions from MadGraph +pythia 6, powheg +pythia 6, powheg +herwig 6, and mc@nlo +herwig 6 (left) and to MadGraph with varied renormalization, factorization, and jet-parton matching scales (right). The inner (outer) vertical bars indicate the statistical (total) uncertainties. The lower part of each plot shows the ratio of the predictions to the data

10 Differential \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b} \overline{\mathrm {b}}}\) (\(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b}}\)) cross sections as a function of the kinematic variables of the additional \(\mathrm {b}\) jets

Figure 15 shows the absolute \(\mathrm {t}\overline{\mathrm {t}}\) differential cross sections in the visible phase space of the \(\mathrm {t}\overline{\mathrm {t}}\) system and the additional \(\mathrm {b}\) jets as a function of the \(p_{\mathrm {T}}\) and \(|\eta |\) of the leading and subleading additional \(\mathrm {b}\) jets, and \(\varDelta R_{\mathrm {b} \mathrm {b}}\) and \(m_{\mathrm {b} \mathrm {b}}\) of the two \(\mathrm {b}\) jets. The uncertainties in the measured cross sections as a function of the \(\mathrm {b}\) jet kinematic variables are dominated by the statistical uncertainties, with values varying from 20–100 %. The results are quantified in Tables 9 and 10 in Appendix B, together with the normalized results. The corresponding migration matrices between the reconstructed and particle levels for the kinematic properties of the additional \(\mathrm {b}\) jets are presented in Fig. 25 in Appendix C for illustration purposes.

The dominant systematic uncertainties are the \(\mathrm {b}\) tagging efficiency and JES, up to 20 % and 15 %, respectively. Other uncertainties have typical values on the order of or below 5 %. The experimental sources of systematic uncertainties affecting only the normalization, which are constrained in the fit, have a negligible impact. The largest model uncertainty corresponds to that from the renormalization and factorization scales of 8 %. The effect of the assumed top quark mass and the PDF uncertainties have typical values of 1–2 %. On average, the inclusion of all the systematic uncertainties increases the total uncertainties by 10 %.

The measured distributions are compared with the MadGraph +pythia 6 prediction, normalized to the corresponding measured inclusive cross section in the same phase space. The measurements are also compared to the predictions from mc@nlo interfaced with herwig 6 and from powheg with pythia 6 and herwig 6. The normalization factors applied to the MadGraph and powheg predictions are found to be about 1.3 for results related to the leading additional \(\mathrm {b}\) jet. The predictions from both generators underestimate the \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b} \overline{\mathrm {b}}}\) cross sections by a factor 1.8, in agreement with the results from Ref. [11]. The normalization factors applied to mc@nlo are approximately 2 and 4 for the leading and subleading additional \(\mathrm {b}\) jet quantities, respectively, reflecting the observation that the generator does not simulate sufficiently large jet multiplicities. All the predictions have slightly harder \(p_{\mathrm {T}}\) spectra for the leading additional \(\mathrm {b}\) jet than the data, while they describe the behaviour of the \(|\eta |\) and \(m_{\mathrm {b} \mathrm {b}}\) distributions within the current precision. The predictions favour smaller \(\varDelta R_{\mathrm {b} \mathrm {b}}\) values than the measurement, although the differences are in general within two standard deviations of the total uncertainty.

Fig. 15
figure 15

Absolute differential \(\mathrm {t}\overline{\mathrm {t}}\) cross section measured in the visible phase space of the \(\mathrm {t}\overline{\mathrm {t}}\) system and the additional \(\mathrm {b}\) jets, as a function of the leading additional \(\mathrm {b}\) jet \(p_{\mathrm {T}}\) (top left) and \(|\eta |\) (top right), subleading additional \(\mathrm {b}\) jet \(p_{\mathrm {T}}\) (middle left) and \(|\eta |\) (middle right), the angular separation \(\varDelta R_{\mathrm {b} \mathrm {b}}\) between the two leading additional \(\mathrm {b}\) jets (bottom left), and the invariant mass \(m_{\mathrm {b} \mathrm {b}}\) of the two \(\mathrm {b}\) jets (bottom right). Data are compared with predictions from MadGraph interfaced with pythia 6, mc@nlo interfaced with herwig 6, and powheg with pythia 6 and herwig 6, normalized to the measured inclusive cross section. The inner (outer) vertical bars indicate the statistical (total) uncertainties. The lower part of each plot shows the ratio of the predictions to the data

The \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b} \overline{\mathrm {b}}}\) production cross sections are compared to the NLO calculation by PowHel +pythia 6 in Fig. 16. In the figure, the prediction is normalized to the absolute cross section given by the calculation of \(20.8 \pm 0.6 \,\text {(stat)} {}^{+7.9}_{-5.4} \text {(scale)}\,{\mathrm{fb}}\). The prediction describes well the shape of the different distributions, while the predicted absolute \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b} \overline{\mathrm {b}}}\) cross section is about 30 % lower than the measured one, but compatible within the uncertainties.

Fig. 16
figure 16

Absolute differential \(\mathrm {t}\overline{\mathrm {t}}\) cross section measured in the visible phase space of the \(\mathrm {t}\overline{\mathrm {t}}\) system and the additional \(\mathrm {b}\) jets, as a function of the second additional \(\mathrm {b}\) jet \(p_{\mathrm {T}}\) (top left) and \(|\eta |\) (top right), the angular separation \(\varDelta R_{\mathrm {b} \mathrm {b}}\) between the two leading additional \(\mathrm {b}\) jets (bottom left), and the invariant mass \(m_{\mathrm {b} \mathrm {b}}\) of the two \(\mathrm {b}\) jets (bottom right). Data are compared with predictions from PowHel +pythia 6. The inner (outer) vertical bars indicate the statistical (total) uncertainties. The lower part of each plot shows the ratio of the calculation to data

The absolute differential cross sections measured in the visible phase space of the additional \(\mathrm {b}\) jets and the full phase space of the \(\mathrm {t}\overline{\mathrm {t}}\) system are presented in Fig. 17 and given in Tables 11 and 12 of Appendix B. The results are corrected for acceptance and dileptonic branching fractions including \(\tau \) leptonic decays (\(6.43 \pm 0.14\)) % [62]. The results are compared to the same predictions as in Fig. 15, which are scaled to the measured cross section, obtained by integrating all the bins of the corresponding distribution. The normalization factor applied to the simulations is similar to the previous one for the results in the visible phase space of the \(\mathrm {t}\overline{\mathrm {t}}\) system. The description of the data by the simulations is similar as well. The total measured \(\sigma _{\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b} \overline{\mathrm {b}}}}\), as well as the agreement between the data and the simulation, is in agreement with the result obtained in Ref. [11]. In the full phase space, the inclusive \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b} \overline{\mathrm {b}}}\) cross section at NLO given by PowHel +pythia 6 corresponds to \(62 \pm 1 \,\text {(stat)} {}^{+23}_{-17} \text {(scale)}\,{\mathrm{fb}}\) (excluding the dileptonic branching fraction correction). The comparison of the differential \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b} \overline{\mathrm {b}}}\) cross section with the NLO calculation is presented in Fig. 18.

Fig. 17
figure 17

Absolute differential \(\mathrm {t}\overline{\mathrm {t}}\) cross section measured in the full phase space of the \(\mathrm {t}\overline{\mathrm {t}}\) system, corrected for acceptance and branching fractions, and the visible phase space of the additional \(\mathrm {b}\) jets, as a function of the leading additional \(\mathrm {b}\) jet \(p_{\mathrm {T}}\) (top left) and \(|\eta |\) (top right), subleading additional \(\mathrm {b}\) jet \(p_{\mathrm {T}}\) (middle left) and \(|\eta |\) (middle right), the angular separation \(\varDelta R_{\mathrm {b} \mathrm {b}}\) between the leading and subleading additional \(\mathrm {b}\) jets (bottom left), and the invariant mass \(m_{\mathrm {b} \mathrm {b}}\) of the two \(\mathrm {b}\) jets (bottom right). Data are compared with predictions from MadGraph interfaced with pythia 6, mc@nlo interfaced with herwig 6, and powheg intefarced with both pythia 6 and herwig 6, normalized to the measured inclusive cross section. The inner (outer) vertical bars indicate the statistical (total) uncertainties. The lower part of each plot shows the ratio of the predictions to the data

Fig. 18
figure 18

Absolute differential \(\mathrm {t}\overline{\mathrm {t}}\) cross section measured in the full phase space of the \(\mathrm {t}\overline{\mathrm {t}}\) system, corrected for acceptance and branching fractions, and the additional \(\mathrm {b}\) jets, as a function of the second additional \(\mathrm {b}\) jet \(p_{\mathrm {T}}\) (top left) and \(|\eta |\) (top right), the angular separation \(\varDelta R_{\mathrm {b} \mathrm {b}}\) between the leading and subleading additional \(\mathrm {b}\) jets (bottom left), and the invariant mass \(m_{\mathrm {b} \mathrm {b}}\) of the two \(\mathrm {b}\) jets (bottom right). Data are compared with predictions from PowHel +pythia 6. The inner (outer) vertical bars indicate the statistical (total) uncertainties. The lower part of each plot shows the ratio of the calculation to data

Differences between the kinematic properties of the additional jets and b jets are expected owing to the different production mechanisms [71] of both processes. The dominant production mechanism of \(\mathrm {p}\mathrm {p}\rightarrow \mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b} \overline{\mathrm {b}}} \) is gluon-gluon (gg) scattering, while in the case of \(\mathrm {p}\mathrm {p}\rightarrow \mathrm {\mathrm {t} \overline{\mathrm {t}}}\mathrm {jj} \), the quark-gluon (qg) channel is equally relevant. The \(|\eta |\) distributions of the additional b jets seem to be more central than the corresponding distributions of the additional jets, see Figs. 10 and  13. This difference can be attributed mainly to the contribution of the production via the qg channel, which favours the emission of jets at larger \(|\eta | \). The distributions of the differential cross section as a function of \(m_{\mathrm {b} \mathrm {b}}\) peak at smaller invariant masses than those as a function of \(m_{\mathrm {jj}}\), presented in Figs. 11 and 14, because of the larger contribution of the gg channel. Given the large uncertainties in the \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b} \overline{\mathrm {b}}}\) measurements, no statistically significant differences can be observed in the shape of the \(p_{\mathrm {T}}\) distributions of the additional b jets compared to the additional jets, shown in Figs. 9 and 12.

11 Additional jet gap fraction

An alternative way to investigate the jet activity arising from quark and gluon radiation is to determine the fraction of events that do not contain additional jets above a given \(p_{\mathrm {T}}\) threshold [5, 12]. A threshold observable, referred to as the gap fraction, is defined as:

$$\begin{aligned} f(p_{\mathrm {T}} ^j)=\frac{N(p_{\mathrm {T}} ^j)}{N_{\text {total}}}, \end{aligned}$$
(2)

where \(N_{\text {total}}\) is the total number of selected events and \(N(p_{\mathrm {T}} ^j)\) is the number of events that do not contain at least j additional jets (apart from the two jets from the \(\mathrm {t}\overline{\mathrm {t}}\) solution hypothesis) above a \(p_{\mathrm {T}}\) threshold, with j corresponding to one or two jets. The measurements are presented as a function of the \(p_{\mathrm {T}}\) of the leading and subleading additional jets, respectively.

A modified gap fraction can be defined as:

$$\begin{aligned} f(H_{\mathrm {T}})=\frac{N(H_{\mathrm {T}})}{N_{\text {total}}}, \end{aligned}$$
(3)

where \(N(H_{\mathrm {T}})\) is the number of events in which the sum of the scalar \(p_{\mathrm {T}}\) of the additional jets \((H_{\mathrm {T}})\) is less than a certain threshold. In both cases, detector effects are unfolded using the MadGraph simulation to obtain the results at the particle level. The additional jets at the generator level are defined as all jets within the kinematic acceptance, excluding the two \(\mathrm {b}\) jets originating from the \(\mathrm {b}\) quarks from top quark decay (see Sect. 7). For each value of the \(p_{\mathrm {T}}\) and \(H_{\mathrm {T}} \) thresholds the gap fraction at the generator level is evaluated, along with the equivalent distributions after the detector simulation and analysis requirements. Given the high purity of the selected events, above 70 % for any bin for the leading additional jet \(p_{\mathrm {T}}\) and \(H_{\mathrm {T}} \), and above 85 % for any bin for the subleading additional jets, a correction for detector effects is applied by following a simpler approach than the unfolding method used for other measurements presented here. The data are corrected to the particle level by applying the ratio of the generated distributions at particle level to the simulated ones at the reconstruction level, using the nominal MadGraph simulation.

The measured gap fraction distributions are compared to predictions from MadGraph interfaced with pythia 6, powheg 6 interfaced with pythia 6 and herwig 6, mc@nlo interfaced with herwig 6, and to the MadGraph predictions with varied renormalization, factorization, and jet-parton matching scales. Figure 19 displays the gap fraction distribution as a function of the \(p_{\mathrm {T}}\) of the leading and subleading additional jets, and \(H_{\mathrm {T}} \). The lower part of the figures shows the ratio of the predictions to the data. The light band indicates the total uncertainty in the data in each bin. The threshold, defined at the value where the data point is shown, is varied from 25\(\,\text {GeV}\) (lower value compared to previous measurements [5]) to 190\(\,\text {GeV}\). In general, MadGraph interfaced with pythia 6 agrees with the data distributions of the three variables, while powheg interfaced with pythia 6 and herwig 6 also provide a good description of the data, though they tend to predict a lower gap fraction than the measured ones. The mc@nlo generator interfaced with herwig 6 describes the data well as a function of the leading additional jet \(p_{\mathrm {T}}\). However, it predicts higher values of the gap fraction as a function of the subleading jet \(p_{\mathrm {T}}\) and \(H_{\mathrm {T}} \). Modifying the renormalization and factorization scales in MadGraph worsens the agreement with data, while variations of the jet-parton matching threshold provide similar predictions as the nominal MadGraph simulation, in agreement with the results shown before.

Fig. 19
figure 19

Measured gap fraction as a function of the leading additional jet \(p_{\mathrm {T}}\) (top row), subleading additional jet \(p_{\mathrm {T}}\) (middle row), and of \(H_{\mathrm {T}} \) (bottom row). Data are compared to predictions from MadGraph, powheg interfaced with pythia and herwig, and mc@nlo interfaced with herwig (left), and to MadGraph with varied renormalization, factorization, and jet-parton matching scales (right). For each bin the threshold is defined at the value where the data point is placed. The vertical bars on the data points indicate the statistical uncertainty. The shaded band corresponds to the statistical and the total systematic uncertainty added in quadrature. The lower part of each plot shows the ratio of the predictions to the data

The results are also compared in Fig. 20 with the recently available simulations, described in Sect. 3, matched to different versions of the parton showering models. The MadGraph and MG5_aMC@NLO generators interfaced with pythia 8 predict up to 10 % lower values of the gap fraction for all the variables, which reflects the fact that those simulations generate larger jet multiplicities, as discussed in Sect. 8. Within the uncertainties, the predictions of the powheg +pythia 8 simulation agree well with data, while the powheg generator (with \(\textsc {hdamp} = m_{\mathrm {t}}\)) interfaced with pythia 6 and herwig 6 tends to overestimate and underestimate the measured values, respectively.

Fig. 20
figure 20

Measured gap fraction as a function of the leading additional jet \(p_{\mathrm {T}}\) (top row), subleading additional jet \(p_{\mathrm {T}}\) (middle row), and of \(H_{\mathrm {T}} \) (bottom row). Data are compared to predictions from MadGraph, interfaced with pythia 6 and pythia 8, and MG5_aMC@NLO interfaced with herwig 6 (left), and to powheg interfaced with different versions of pythia and herwig 6 (right). For each bin the threshold is defined at the value where the data point is placed. The vertical bars on the data points indicate the statistical uncertainty. The shaded band corresponds to the statistical and the total systematic uncertainty added in quadrature. The lower part of each plot shows the ratio of the predictions to the data

The gap fraction is also measured in different \(|\eta |\) regions of the additional jets, with the results presented in Figs. 21, 22 and 23 as a function of the leading additional jet \(p_{\mathrm {T}}\), subleading additional jet \(p_{\mathrm {T}}\), and \(H_{\mathrm {T}} \), respectively. In general, the gap fraction values predicted by the simulations describe the data better in the higher \(|\eta |\) ranges. The values given by MadGraph and powheg interfaced with pythia 6 are slightly below the measured ones in the central region for the leading \(p_{\mathrm {T}}\) jet and \(H_{\mathrm {T}} \), while mc@nlo +herwig 6 yields higher values of the gap fraction. In the case of the subleading jet \(p_{\mathrm {T}}\), all predictions agree with the data within the uncertainties, except for mc@nlo +herwig 6 in the more central regions. Variations of the jet-parton matching threshold do not have a noticeable impact on the gap fraction, while MadGraph with the varied renormalization and factorization scales provides a poorer description of the data.

Fig. 21
figure 21

Measured gap fraction as a function of the leading additional jet \(p_{\mathrm {T}}\) in different \(\eta \) regions. Data are compared to predictions from MadGraph, powheg interfaced with pythia 6 and herwig 6, and mc@nlo interfaced with herwig 6 (left) and to MadGraph with varied renormalization, factorization, and jet-parton matching scales (right). For each bin the threshold is defined at the value where the data point is placed. The vertical bars on the data points indicate the statistical uncertainty. The shaded band corresponds to the statistical uncertainty and the total systematic uncertainty added in quadrature. The lower part of each plot shows the ratio of the predictions to the data

Fig. 22
figure 22

Measured gap fraction as a function of the subleading additional jet \(p_{\mathrm {T}}\) in different \(|\eta |\) regions. Data are compared to predictions from MadGraph, powheg interfaced with pythia 6 and herwig 6, and mc@nlo interfaced with herwig 6 (left) and to MadGraph with varied with varied renormalization, factorization, and jet-parton matching scales (right). For each bin the threshold is defined at the value where the data point is placed. The vertical bars on the data points indicate the statistical uncertainty. The shaded band corresponds to the statistical uncertainty and the total systematic uncertainty added in quadrature. The lower part of each plot shows the ratio of the predictions to the data

Fig. 23
figure 23

Measured gap fraction as a function of \(H_{\mathrm {T}} \) in different \(\eta \) regions. Results in data are compared to the nominal MadGraph signal sample, powheg and mc@nlo (left) and to the samples with varied renormalization, factorization, and jet-parton matching scales (right). For each bin the threshold is defined at the value where the data point is placed. The vertical bars on the data points indicate the statistical uncertainty. The shaded band corresponds to the statistical uncertainty and the total systematic uncertainty added in quadrature. The lower part of each plot shows the ratio of the predictions to the data

The total systematic uncertainty in the gap fraction distributions is about 5 % for low values of the threshold (\(p_{\mathrm {T}}\) or \(H_{\mathrm {T}} \)) and decreases to \({<}0.5~\%\) for the highest values. The measurement of the gap fraction as a function of \(H_{\mathrm {T}} \) has larger uncertainties because of the impact of the lower-momentum jets that have a significantly larger uncertainty, as discussed in Sect. 9. The uncertainty in JES is the dominant source of systematic uncertainty, corresponding to approximately 4 % for the smallest \(p_{\mathrm {T}}\) and \(H_{\mathrm {T}} \) values. Other sources with a smaller impact on the total uncertainty are the \(\mathrm {b}\) tagging efficiency, JER, pileup, and the simulated sample used to correct the data to the particle level.

12 Summary

Measurements of the absolute and normalized differential top quark pair production cross sections have been presented using pp collisions at a centre-of-mass energy of 8\(\,\text {TeV}\), corresponding to an integrated luminosity of 19.7\(\,\text {fb}^\text {-1}\), in the dilepton decay channel as a function of the number of jets in the event, for three different jet \(p_{\mathrm {T}}\) thresholds, and as a function of the kinematic variables of the leading and subleading additional jets. The results have been compared to the predictions from MadGraph interfaced with pythia 6, powheg interfaced with both pythia 6 and herwig 6, mc@nlo interfaced with herwig 6, and MadGraph samples with varied renormalization, factorization, and jet-parton matching scales. In general, all these generators are found to give a reasonable description of the data.

The MadGraph and powheg generators interfaced with pythia 6 describe the data well for all measured jet multiplicities; while mc@nlo interfaced with herwig 6 generates lower multiplicities than observed for the lower-\(p_{\mathrm {T}}\) thresholds. The prediction from MadGraph with varied renormalization and factorization scales does not provide an improved description of the data compared to the nominal simulation.

These results are also compared to the predictions from powheg with the hdamp parameter set to the top quark mass interfaced with pythia 6, pythia 8, and herwig 6, which provide a reasonable description of the data within the uncertainties, and the predictions from MadGraph and MG5_aMC@NLO interfaced with pythia 8, which generate higher jet multiplicities for all the \(p_{\mathrm {T}}\) thresholds.

The measured kinematic variables of the leading and subleading additional jets are consistent with the various predictions. The simulations also describe well the data distributions of the leading additional jet \(p_{\mathrm {T}}\) and \(H_{\mathrm {T}} \), although they tend to predict higher \(p_{\mathrm {T}}\) values and more central values in \(\eta \). MadGraph with varied parameters yields similar predictions, except for varying the renormalization and factorization scales, which tends to give higher \(H_{\mathrm {T}} \) values. The mc@nlo generator predicts lower yields than observed for the subleading additional jet \(p_{\mathrm {T}}\).

The uncertainties in the measured \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b} \overline{\mathrm {b}}}\) (\(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b}}\)) absolute and normalized differential cross sections as a function of the \(\mathrm {b}\) jet kinematic variables are dominated by the statistical uncertainties. In general, the predictions describe well the shape of the measured cross sections as a function of the variables studied, except for \(\varDelta R_{\mathrm {b} \mathrm {b}}\), where they favour smaller values than the measurement. The predictions underestimate the total \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b} \overline{\mathrm {b}}}\) cross section by approximately a factor of 2, in agreement with previous measurements [11]. The calculation by PowHel [19] describes well the shape of the distributions, while the predicted absolute cross section is about 30 % lower, but compatible with the measurements within the uncertainties.

The gap fraction has been measured as a function of the \(p_{\mathrm {T}}\) of the leading and subleading additional jets and \(H_{\mathrm {T}} \) of the additional jets in different \(\eta \) ranges. For a given threshold value, the gap fraction as a function of \(H_{\mathrm {T}} \) is lower than the gap fraction as a function of the \(p_{\mathrm {T}}\) of the leading additional jet, showing that the measurement is probing multiple quark and gluon emission. Within the uncertainties, all predictions describe the gap fraction well as a function of the momentum of the first additional jet, while mc@nlo interfaced with herwig fails to describe the gap fraction as a function of the subleading additional jet \(p_{\mathrm {T}}\) and \(H_{\mathrm {T}} \). In general, MadGraph with decreased renormalization and factorization scales more poorly describes the observed gap fraction, while varying the jet-parton matching threshold provides a similar description of the data. The MadGraph and MG5_aMC@NLO generators interfaced with pythia 8 predict lower values than measured. The powheg simulation with \(\textsc {hdamp} = m_{\mathrm {t}}\) interfaced with pythia 8 is consistent with the data, while the simulation interfaced with herwig 6 and pythia 6 tends to worsen the comparison with the measurement.

In general, the different measurements presented are in agreement with the SM predictions as formulated by the various event generators, within their uncertainties. The correct description of \(\mathrm {t}\overline{\mathrm {t}}\) +jets production is important since it constitutes a major background in searches for new particles in several supersymmetric models and in \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {H}}\) processes, where the Higgs boson decays into \(\mathrm {b} \overline{\mathrm {b}} \). The \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b} \overline{\mathrm {b}}}\) (\(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {b}}\)) differential cross sections, measured here for the first time, also provide important information about the main irreducible background in the search for \(\mathrm {\mathrm {t} \overline{\mathrm {t}} \mathrm {H} \,(\mathrm {b} \overline{\mathrm {b}})}\).