Measurement of the $\mathrm{t\bar{t}}\mathrm{b\bar{b}}$ production cross section in the all-jet final state in pp collisions at $\sqrt{s} =$ 13 TeV

A measurement of the production cross section of top quark pairs in association with two b jets ($\mathrm{t\bar{t}}\mathrm{b\bar{b}}$) is presented using data collected in proton-proton collisions at $\sqrt{s} =$ 13 TeV by the CMS detector at the LHC corresponding to an integrated luminosity of 35.9 fb$^{-1}$. The cross section is measured in the all-jet decay channel of the top quark pair by selecting events containing at least eight jets, of which at least two are identified as originating from the hadronization of b quarks. A combination of multivariate analysis techniques is used to reduce the large background from multijet events not containing a top quark pair, and to help discriminate between jets originating from top quark decays and other additional jets. The cross section is determined for the total phase space to be 5.5 $\pm$ 0.3 (stat) ${}^{+1.6}_{-1.3}$ (syst) pb and also measured for two fiducial $\mathrm{t\bar{t}}\mathrm{b\bar{b}}$ definitions. The measured cross sections are found to be larger than theoretical predictions by a factor of 1.5-2.4, corresponding to 1-2 standard deviations.


Introduction
At the CERN LHC, top quark pairs are produced with copious amounts of additional jets, including those resulting from the hadronization of b quarks (b jets).Top quark pair production in association with a pair of b jets, ttbb, is challenging to model because of the very different energy scales for the b jets produced in association with the tt system and that of tt system [1], and because of the small but nonnegligible mass of the b quark.Improving the accuracy and the precision of perturbative calculations in quantum chromodynamics (QCD) for this process is crucial, since it represents an important background for numerous searches or other measurements at the LHC.In particular, tt production in association with a Higgs boson (ttH), where the Higgs boson decays to bb, suffers from an irreducible ttbb background [2][3][4][5][6][7].Searches for four top quark production (tttt) are also affected by this background [8][9][10].The two latter processes provide direct access to the top quark Yukawa coupling, a crucial parameter of the standard model [11,12].An improved understanding of the ttbb process would help reduce the uncertainty in such measurements.
Calculations of the production cross section of tt in association with jets have been performed at next-to-leading order (NLO) in QCD and matched with parton showers for up to two additional massless partons in the matrix element [13][14][15].The ttbb cross section at NLO, matched with parton showers, has also been calculated for massless b quarks (five-flavour scheme, 5FS) [16], and has recently become available for massive b quarks (four-flavour scheme, 4FS) [17][18][19].A comparison of the measurements of the ttbb cross section with such calculations provides valuable guidance to improve the different frameworks.The ttbb cross section has been measured previously at √ s = 8 and 13 TeV by the ATLAS and CMS Collaborations, in events containing one or two charged leptons [20][21][22][23][24].This Letter focuses on the all-jet final state of the tt system, where each top quark decays into three jets, leading to a signature of four b jets and four light-quark jets for the ttbb system.This final state is favoured by a large branching fraction and provides a complete reconstruction of top quarks, as opposed to other decay channels of the top quark pairs.Moreover, the main uncertainties affecting the sensitivity in this measurement are different than those affecting final states containing leptons, therefore providing complementary information.However, the all-jet channel also suffers from a large background from multijet production, as well as from the difficulty of identifying jets that originate from decaying top quarks.Multivariate analysis techniques are developed and implemented to mitigate these problems.The ttbb cross section is measured using data collected by the CMS detector in pp collisions at √ s = 13 TeV, corresponding to an integrated luminosity of 35.9 fb −1 [25].

The CMS detector and event simulation
The central feature of the CMS apparatus is a superconducting solenoid of 6 m internal diameter, providing a magnetic field of 3.8 T. A silicon pixel and strip tracker, a lead tungstate crystal electromagnetic calorimeter (ECAL), and a brass and scintillator hadron calorimeter (HCAL), each composed of a barrel and two endcap sections reside within the solenoid field.Forward calorimeters extend the pseudorapidity coverage provided by the barrel and end detectors.Muons are detected in gas-ionization chambers embedded in the steel flux-return yoke outside the solenoid.A more detailed description of the CMS detector, together with a definition of its coordinate system and kinematic variables, can be found in Ref. [26].Samples of tt events are simulated at NLO in QCD using POWHEG (v2) [27][28][29][30].These samples include ttbb events, where the additional b jets are generated by the parton shower.Single top quark production in the t channel or in association with a W boson, and ttH production are simulated at NLO with POWHEG [31][32][33].Production of W or Z bosons in association with jets (V+jets), as well as QCD multijet events, are simulated at leading order (LO) with MADGRAPH5 aMC@NLO (v2.2.2) [14], and the MLM merging scheme [34].The MADGRAPH5 aMC@NLO generator is used at NLO for simulating associated production of top quark pairs with W or Z bosons (ttV).Diboson processes (WW, WZ and ZZ) are simulated at LO using PYTHIA (v8.219) [35].
All simulated events are processed with PYTHIA for modelling of the parton showering, hadronization, and underlying event (UE).The NNPDF 3.0 [36] parton distribution functions (PDFs) are used throughout, at the same perturbative order as used by the event generators.The CUETP8M1 UE tune [37] is used for all processes except for the tt, ttH and single top quark processes.For these, an updated version of the tune is used (CUETP8M2T4), in which an adjusted value of the strong coupling constant is used in the description of initial-state radiation [38].Simulation of the CMS detector response is based on GEANT4 (v9.4) [39].Additional pp interactions in the same or neighbouring bunch crossings (pileup) are simulated with PYTHIA and overlaid with hard-scattering events according to the pileup distribution measured in data.
The various simulated processes are normalized to state-of-the-art predictions for the production cross sections.The tt, V+jets, single top quark, and W + W − samples are normalized to next-to-NLO (NNLO) precision in QCD [40][41][42][43], while remaining processes such as ttV, ttH, and other diboson production are normalized to NLO in QCD [14,44].

Definitions of fiducial phase space
The ttbb production cross section is measured for three different phase space definitions.Two definitions for ttbb events in the fiducial phase space, matching the detector acceptance, are considered: one that is based exclusively on stable generated particles after hadronization (parton-independent), and one that also uses parton-level information after radiation emission (parton-based).The former facilitates comparisons with predictions from event generators, while the latter is closer to the approach taken by searches for ttH production to define the contribution from the ttbb process.The cross section is reported for the total phase space by correcting the parton-based fiducial cross section by the experimental acceptance.
Particle-level jets are defined by clustering stable generated final-state particles, excluding neutrinos, using the anti-k T algorithm [45,46] with a distance parameter of 0.4.These jets are defined unambiguously as b or c jets by rescaling the momenta of generated b and c hadrons to a negligible value, while preserving their direction, and including them in the clustering procedure [47].A jet is labelled b jet if it is matched to at least one b hadron, and labelled c jet if matched with at least one c hadron and no b hadron.
Events in the generated tt sample are divided into exclusive categories according to the flavour of the jets that do not originate from the decay of top quarks, which we refer to as "additional" jets.The b or c jets are considered to originate from a top quark if one of the clustered b or c hadrons features a top quark in its simulation history.Additional jets are required to have a transverse momentum p T > 20 GeV, and absolute pseudorapidity |η| < 2.4.No explicit requirement on the b hadron kinematic variables is used.Events are categorized as ttbb if they contain at least two additional b jets, which defines the total phase space for which the ttbb cross section is measured.Events with a single additional b jet are categorized as ttb (tt2b) if that b jet is matched with exactly one (at least two) b hadron(s).The ttb events correspond to ttbb events where one of the additional b jets fails the above kinematic requirements, while tt2b events arise from collinear gluon splittings.If no b jets are present but at least one additional c jet is present the event is referred to as ttcc; all remaining events are denoted ttjj.
For the parton-based definition of the ttbb fiducial phase space, at least eight jets with p T > 20 GeV and |η| < 2.4 must be present, of which at least six have p T > 30 GeV.At least four of these jets must be b jets, and at least two of those must not originate from top quarks.This last requirement is removed for the parton-independent fiducial definition, in order to be independent of the origin of the b jets, and thus of the simulated parton content.Some ttbb events in the total phase space failing the fiducial requirements may still be reconstructed and selected because of resolution effects, and are referred to as out-of-acceptance.They correspond to 16% of all reconstructed ttbb events.

Event reconstruction and selection
The particle-flow algorithm [48] aims to reconstruct and identify each particle in an event, with an optimized combination of information from the various elements of the CMS detector.The primary pp interaction vertex is taken to be the reconstructed vertex with the largest sum of the p 2 T of the objects associated to that vertex, where the considered objects are those returned by a jet clustering algorithm [45,46] applied to the tracks assigned to the vertex, and the associated missing transverse momentum, taken as the negative vector sum of the p T of those objects.The energy of photons is obtained from the ECAL measurement.The energy of electrons is determined from a combination of the electron momentum at the primary interaction vertex as determined by the tracker, the energy of the corresponding ECAL cluster, and the energy sum of all bremsstrahlung photons spatially compatible with originating from the electron track.The p T of muons is obtained from the curvature of the corresponding tracks.The energy of charged hadrons is determined from a combination of their momentum measured in the tracker and the matching ECAL and HCAL energy deposits, corrected for zero-suppression effects and for the response function of the calorimeters to hadronic showers.The energy of neutral hadrons is obtained from the corresponding corrected ECAL and HCAL energies.
For each event, hadronic jets are clustered from the reconstructed particles using the anti-k T algorithm with a distance parameter of 0.4.The jet momentum is determined as the vectorial sum of all particle momenta in the jet, and is found from simulation to be within 5 to 10% of the true momentum over the whole p T spectrum and detector acceptance.Pileup interactions can contribute additional tracks and calorimetric energy depositions to the jet momentum.To mitigate this effect, tracks identified to be originating from pileup vertices are discarded and an offset correction is applied to correct for remaining contributions [47].Jet energy corrections are derived from simulation to bring the average measured response of a jet to that of a particle-level jet.In situ measurements of the momentum balance in dijet, photon+jet, Z+jet, and multijet events are used to account for any residual differences in jet energy scale in data and simulation [49].The data used for these measurements are independent of those used for the present Letter.
A combined secondary vertex b tagging algorithm (CSVv2) is used to identify jets originating from the hadronization of b quarks [50], with an efficiency for identifying b jets in simulated tt events of about 65%.The misidentification probability is about 10 and 1% for c and lightflavour jets, respectively, where the latter refers to jets originating from the hadronization of u, d, s quarks or gluons.The distribution of the discriminator score for b and light-flavour jets in the simulation is calibrated to match the distribution measured in control samples of tt events with exactly two leptons (electrons or muons) and two jets, and Z bosons produced in association with jets where the Z bosons decay to pairs of electrons or muons.The calibration is achieved by reweighting events using scale factors that are parameterized by the jet flavour, p T , |η|, and b tagging discriminator score [50].Data are collected using two triggers [51], both requiring at least six jets with |η| < 2.4.The first (second) trigger considers jets with p T > 40 (30) GeV, and requires that the jet scalar p T sum, H T , exceeds 450 (400) GeV and that at least one (two) of the jets is (are) b tagged.The efficiency of these triggers is measured in simulation, as well as in a data control sample collected using independent single-muon triggers.The trigger efficiency in simulation is corrected to match the efficiency observed in the data by reweighting events using scale factors defined as the ratio between the efficiencies in the data and simulation.For events satisfying the preselection criteria detailed below, the trigger efficiency is above 95%.
An offline preselection is applied to data and simulated events, by requiring the presence of at least six jets with p T > 40 GeV and |η| < 2.4, of which at least two are b tagged, and H T > 500 GeV.Additional jets in the events are considered if they satisfy the requirements p T > 30 GeV and |η| < 2.4.Events are vetoed if they contain electrons or muons with p T > 15 GeV and |η| < 2.4 that satisfy highly efficient identification criteria [52,53] and are isolated from hadronic activity.About 20% of the ttbb events in the fiducial phase space pass the offline selection.

Multivariate analysis
The final state considered in this analysis suffers from a large background from multijet production, as well as from the difficulty to identify which jets do not stem from top quark decays.To address these challenges and improve the sensitivity to the ttbb signal, several multivariate analysis tools have been employed.
The multijet background can be discriminated from tt production by observing that the latter is expected to contain four light-quark jets from W boson decays per event, whereas the former is enriched in gluon jets.Gluon and quark jets are separated using a quark-gluon likelihood (QGL) variable, based on jet substructure observables [54,55].Using the individual jet QGL values, the likelihood of an event to contain N q light-quark jets and N g gluon jets is defined as where the sums run over all possible assignments of N q jets to quarks (indices k) and N g jets to gluons (indices m), ζ i is the QGL discriminant of the i th jet, and f q and f g are the probability densities for ζ i under the hypothesis of (u, d, s, or c) quark or gluon origin, respectively.When computing L(N q , N g ), b-tagged jets are not considered.Based on the event likelihoods with N q = 4 and N g = 0, as well as N q = 0 and N g = 4, the QGL ratio (QGLR) is defined as QGLR = L(4, 0)/(L(4, 0) + L(0, 4)).Other values for N q and N g have been tried but led to reduced discrimination between multijet and tt production.We correct the modelling of the QGL in the simulation by reweighting each event based on the quark or gluon origin and the QGL value of all jets in the event, where the weights are measured using data samples enriched in Z+jets and dijet events [55].After applying this correction, a good agreement is found between data and simulation.
To address the large combinatorial ambiguity in identifying the additional jets in the events, we have trained a boosted decision tree (BDT) using the TMVA package [56], henceforth referred to as the "permutation BDT".In events with eight reconstructed jets, there are 28 ways to select six of those as originating from the all-jet decay of a top quark pair, and there are 90 ways to match those six jets to the six partons from the top quark decay chains.Some permutations are indistinguishable and are not considered, i.e. permutations of two jets assigned to a W boson decay are not considered, and neither are the permutations of three jets assigned to a t or t decay.To reduce the large number of permutations, the least favoured ones are rejected using a χ 2 variable quantifying the compatibility of the invariant masses of the different jet pairings with those of the particles they should come from, defined as where m (... ) denotes the invariant mass of the given jets, and σ W = 10.9GeV and σ t = 17.8 GeV are the experimental resolutions in the two-and three-jet invariant masses, respectively.The masses entering the equation are m t = 172.3GeV and m W = 80.2 GeV, measured from the generated tt system after reconstruction.The BDT is trained using simulated tt events after applying the above preselection criteria, requiring the presence of at least seven jets, and reducing the number of permutations by requiring that χ 2 < 33.38, corresponding to a p-value P(χ 2 ) of 10 −6 for a χ 2 distribution with four degrees of freedom.Events for which no permutation satisfies this requirement are rejected.The correct jet-parton assignment is considered as a signal in the training, while all other distinguishable combinations are treated as background.
Input variables used for the BDT include jet b tagging discriminator scores and kinematic quantities, such as invariant masses of pairs and triplets of jets, angular openings between jets, and the transverse momenta of jets.For each permutation, only quantities pertaining to the six jets assumed to originate from the top quarks are used in the training.The permutation yielding the highest BDT score is used for the rest of the analysis.For tt events with eight jets where all six jets from the top quark decays have been selected, the permutation BDT identifies the correct permutation with about 60% efficiency.
As a further handle to reduce the multijet background, we have trained a second BDT to discriminate this background from inclusive tt+jets production.While supervised training of multivariate classifiers relies on samples of simulated events, the poor modelling of multijet production and the insufficient size of the available simulated samples limit the achievable discrimination power.A proposed method to alleviate these shortcomings is a classification without labels (CWoLa) [57].In this weakly supervised approach, the classifier is trained using data, whereby one region in the data is treated as background and another independent region is treated as signal.In the limit of large training sample the resulting classifier converges to the optimal classifier to distinguish between signal and background, provided the two following conditions are fulfilled [57].First, the relative rates of the actual signal and background processes should be different in the two regions.Second, the distributions of the variables entering the CWoLa classifier should be independent of the quantity used to define the two regions, for both the signal and background processes.The CWoLa BDT is trained using a sample of data with exactly seven jets, where two independent regions are defined by requiring that the QGLR is below or above 0.95.The first and second regions are expected to contain about 10 and 20% of tt events, respectively.Variables used for constructing the CWoLa BDT are kinematic quantities similar to those used in the permutation BDT, the output value of the permutation BDT, and the b tagging discriminator scores of the two jets identified by the permutation BDT as the b jets originating from the top quark decays.Only the six jets identified by the permutation BDT as coming from the top quark decays are used to define the CWoLa BDT input variables.The performance of the resulting classifier, measured in the region with at least eight jets, is found to be comparable to that of a supervised classifier trained using simulated samples.

Cross sections
To measure the ttbb cross section we require, in addition to the preselection criteria, the presence of at least eight jets, and P(χ 2 ) > 10 −6 .The distributions in the QGLR and of the CWoLa BDT discriminants for selected events are shown in Fig. 1.The cross section is extracted from a binned maximum likelihood fit to a two-dimensional distribution (referred to as 2DCSV) constructed using the largest and second-largest b tagging discriminator scores among the jets determined to be additional jets by the permutation BDT.In order to increase the signal purity and the precision in the measurement, we define a signal region (SR) by requiring that the CWoLa BDT score be above 0.5, and the QGLR be above 0.8.These thresholds are optimized to obtain the best expected precision in the cross section.About 20% of the ttbb signal that passes the offline preselection is selected into the SR.Both are after preselection, requiring P(χ 2 ) > 10 −6 and at least eight selected jets.All the contributions are based on simulation.The multijet contribution is scaled to match the total yields in data, after the other processes including the ttbb signal have been normalized to their corresponding theoretical cross sections.This choice takes into account only the effect of the shape variation from the multijet background.The small backgrounds include ttV, ttH, single top quark, V+jets, and diboson production.The lower panels show the ratio between the observed data and the predictions.The dashed lines indicate the boundaries between the signal and control regions defined in Section 6. Hatched bands indicate the statistical uncertainty in the predictions without considering the systematic sources, dominated by the uncertainties in the simulated multijet background.Underflow and overflow events were added to the first and last bins, respectively.
The multijet background is also estimated from data.Three independent control regions (CRs), orthogonal to the SR, are defined by inverting the requirements on the CWoLa BDT and the QGLR: the CR1 (BDT > 0.5, QGLR < 0.8), the CR2 (BDT < 0.5, QGLR < 0.8), and the CR3 (BDT < 0.5, QGLR > 0.8).For multijet production, the CWoLa BDT score and the QGLR are nearly independent, so that in each bin i of the 2DCSV distribution the number of multijet events in the SR, N SR i , can be estimated from the number of multijet events in the CRs as This relationship is a consequence of the choice of variables entering the CWoLa BDT, which were required to be independent of the QGLR in order to satisfy the hypotheses of the CWoLa method.In order to properly take into account the small but non-negligible signal contribution in the CRs, the fit to extract the cross section is performed in all four regions, with the multijet rates N CR1 i , N CR2 i , and N CR3 i free to vary in the fit.The assumption of Eq. ( 2) on which this estimation relies is confirmed using the simulation.In addition, we verify that Eq. ( 2) is also satisfied in the data for kinematic distributions, such as the invariant mass of the reconstructed W bosons and top quarks, where for each bin of these distributions the multijet yields are estimated by taking the difference between the observed yields in data and the predicted yields of all simulated processes.Finally, we validate Eq. ( 2) using alternative definitions of the four regions in the plane formed by the QGLR and the CWoLa BDT, excluding the SR as defined above.The outcome of goodness-of-fit tests of the 2DCSV distribution was also positive for each of the alternative region definitions.
The data are fitted using a profiled maximum likelihood technique, where the likelihood is built as a product of independent Poisson likelihoods, defined for each bin i of the 2DCSV distributions in the four event regions using the following expression for the number of events in bin i: where µ is a signal strength parameter, defined by the ratio of observed to expected signal, T k i is the expected yield for process k in bin i, "sig" includes the contributions from ttbb, tt2b, and ttb, and θ is a vector of nuisance parameters affecting the predicted yields of the various processes introduced to model the systematic uncertainties described in the next section.The parameters N i are used to estimate the multijet background from the combined fit of the four regions; they are free parameters in the CRs and are given by Eq. ( 2) in the SR.The likelihood also features constraint terms for each of the nuisance parameters considered in the fit.Different templates are constructed from ttbb events matching the fiducial requirements and from events failing these requirements.For the fiducial ttbb templates, the effect of nuisance parameters corresponding to theoretical uncertainties is normalized such that the ttbb cross section in the fiducial phase space is preserved, i.e.only shape variations within that phase space and their impact on the reconstruction efficiency are taken into account.No such requirement is made for the other templates.The uncertainty in the measured cross section is obtained by profiling the nuisance parameters.As described in the next section, some uncertainties are not profiled and are added in quadrature with the uncertainty obtained from the fit.The fit is repeated for each of the two fiducial phase-space definitions for ttbb events described in Section 3, leading to different in-and out-of-acceptance ttbb templates.The total ttbb cross section is obtained by dividing the cross section for the parton-based fiducial phase space by the acceptance, estimated using POWHEG+PYTHIA to be (29.4± 1.8)%.Uncertainties affecting this acceptance correction are detailed in the next section.

Systematic uncertainties
Several sources of systematic uncertainties affecting the predictions for the signal and background processes entering the analysis are considered.These uncertainties may affect the normalization of the templates entering the fit, or may alter both their shape and their normalization.The migration of events between the four regions is taken into account when relevant.Experimental sources of uncertainties are taken to be fully correlated for all signal and background distributions estimated using the simulation, while only a subset of theoretical uncertainties are correlated among the tt+jets components.
The modelling of the shape of the b tagging discriminator in the simulation represents an important source of systematic uncertainty.Several uncertainties in the calibration of the b tagging discriminator distribution are propagated independently to the shape and normalization of the 2DCSV templates.These are related to the uncertainty in the contamination by light-(heavy-) flavour jets in the control samples used for the measurement of heavy-(light-) jet correction factors, as well as to the statistical uncertainty in these measurements [50].Since no dedicated measurement is performed for c jets, the uncertainty in the shape of the b tagging discriminator distribution for c jets is conservatively taken to be twice the relative uncertainty considered for b jets.In total, six different nuisance parameters are introduced to estimate the uncertainty arising from b tagging.
We evaluate the effect of the uncertainty in the jet energy scale (JES) and jet energy resolution (JER) by shifting the jet four-momenta using correction factors that depend on jet p T and |η| for the JES, and jet |η| for the JER [49].The calibration of the JES is affected by several sources of uncertainty, which are propagated independently to the measurement.The uncertainty in the JES is also propagated to the b tagging calibration, and the resulting effect on the distribution of the b tagging discriminators is taken to be correlated with the effect on the jet momenta.
Uncertainties pertaining to the QGL are estimated conservatively by removing or doubling the scale factors applied to correct the distribution of the QGL in the simulation [55].The uncertainty in the integrated luminosity is evaluated to be 2.5% [25].Uncertainties in the trigger efficiency are estimated by varying the trigger scale factors by their uncertainty, as determined from the efficiency measurements in data and simulation.The uncertainty in the modelling of pileup is estimated by reweighting simulated events to yield different distributions of the expected number of pileup interactions, obtained by varying the total inelastic pp cross section by 4.6% [58].We take into account the limited size of the simulated samples by varying independently the predicted yields in every bin by their statistical uncertainties.
Theoretical uncertainties in the modelling of the tt+jets process enter this analysis both through the efficiency to reconstruct and select ttbb events, and through the contamination from ttcc and ttjj backgrounds.The uncertainties in the renormalization and factorization scales (µ R and µ F , respectively) are estimated by varying both scales independently by a factor of two up or down in the event generation, omitting the two cases where the scales are varied in opposite directions, and taking the envelope of the six resulting variations.Likewise, the uncertainties related to the choice of the scale in the parton shower is evaluated by varying the scale in the initial-state shower by factors of 0.5 and 2, and the scale in the final-state shower by factors of √ 2 and 1/ √ 2. Propagation of the uncertainties associated with the PDFs, as well as with the value of the strong coupling in the PDFs, has been achieved by reweighting generated events using variations of the NNPDF 3.0 set [36].The impact of the choice of the matching scale h damp = 1.58m t between the matrix-element generator and the parton shower in POWHEG is evaluated using simulated samples generated with different choices of h damp = m t and 2.24m t [38].We evaluate the uncertainty related to the UE tune by varying the tune parameters according to their uncertainties.The uncertainty from the modelling of colour reconnection in the final state is evaluated by considering four alternatives to the PYTHIA default, which is based on multiple-parton interactions (MPI) with early resonance decays (ERD) switched off.These alternatives are an MPI-based scheme with ERD switched on, a QCD-inspired scheme [59], and a gluon-move scheme with ERD either off or on [60].All the alternative models were tuned to LHC data [61].It has been verified that the selection efficiency obtained from the nominal tt simulation, in which additional b jets are generated by the parton shower, is in agreement within estimated modelling uncertainties with that obtained using a sample of ttbb events generated at NLO in QCD with massive b quarks (4FS) [19].Since the spectrum of the top quark p T is known to be softer in the data than in the simulation, we evaluate the effect of this mismodelling by reweighting the generated events to match the top quark p T distribution measured in data [62].The latter two uncertainties are not evaluated using profiled nuisance parameters, but by repeating the measurement using varied signal and background predictions.The differences in the measured cross sections are taken as the corresponding uncertainties and are added in quadrature with the uncertainty obtained from the profile likelihood.Uncertainties related to the µ R and µ F scales, the parton shower scale, and the h damp choice are taken to be uncorrelated for the ttbb, ttb, tt2b, ttcc and ttjj templates, while the other modelling uncertainties are taken to be correlated for all tt events.In addition to the aforementioned modelling uncertainties, we assign an uncertainty of 50% to the normalization of the ttcc background to cover the lack of precise measurements of this process.The results are stable when doubling that uncertainty.
Compared to tt+jets and multijet production, the contribution of other background processes such as ttV, ttH, V+jets, diboson, and single top quark production is small.We assign uncertainties to their predicted rates based on the PDF and µ R /µ F scale uncertainties in their theoretical cross sections.
Table 1 summarizes the contributions of the various sources of systematic uncertainties to the total uncertainty in the cross sections measured in the fiducial phase space.The theoretical uncertainty in the acceptance from the various sources listed above is estimated to be 6%, and is added in quadrature with the uncertainty in the parton-based fiducial cross section to yield the systematic uncertainty in the total ttbb cross section.

Results
The result of the maximum likelihood fit described in Section 6 is shown in Fig. 2 for the 2DCSV distributions in the four analysis regions.The contribution from multijet production nearly matches the differences between the yields in data and from the other processes in the CR1, CR2, and CR3 because it is estimated from the data in the four regions according to the method described in the previous section.The measured cross section for the two ttbb definitions in the fiducial phase space, as well as for the total phase space introduced in Section 3, are given in Table 2.The measurement uncertainty is dominated by the systematic effects from the simulation sample sizes, QGL corrections, and µ R and µ F dependences on changes in scale.
Because of the large overlap between the two definitions of the ttbb fiducial phase space, the measured cross sections are numerically equal at the quoted precision.The measurements are compared with NLO predictions from POWHEG for inclusive tt production interfaced with either PYTHIA or HERWIG++ (v2.7.1) [63], using the EE5C UE tune [64] for the latter.Predictions from MADGRAPH5 aMC@NLO at NLO interfaced with PYTHIA for tt production with up to two extra massless partons (5FS) merged using the FxFx scheme [15], and for ttbb production with massive b quarks (4FS), are also compared with the measurements.The predicted cross sections are not rescaled by any NLO to NNLO K-factor, which for inclusive tt production amounts to 1.1-1.15[40].Measured and predicted cross sections are shown in Fig. 3.The predictions underestimate the measured cross section by a factor of 1.5-2.4,corresponding to differences of 1-2 standard deviations.This is consistent with the results from Refs.[20][21][22][23][24].For clarity, the two-dimensional distribution with largest and next-to-largest b tagging discriminant scores for the additional jets have been unrolled to one dimension, and the resulting bins ordered according to increasing values of the ratio between expected signal and background yields in each bin of the SR.The small backgrounds include ttV, ttH, single top quark, V+jets, and diboson production.Hatched bands correspond to uncertainties.The bottom panels show the pull distribution.The pull is defined as the bin by bin difference between data and predicted yields after the fit, divided by the uncertainties accounted for correlations between data and predictions after the fit.Table 1: The considered sources of systematic uncertainties and their respective contributions to the total systematic uncertainty in the measured ttbb cross section for the two defined ttbb fiducial phase spaces.The upper (lower) portion of the table lists uncertainties related to the experimental conditions (theoretical modelling).The numbers are obtained by taking the difference in quadrature of the profile likelihood width when fixing nuisance parameters corresponding to a given source of uncertainty and leaving the others free to vary.

Summary
The first measurement of the ttbb cross section in the all-jet final state was presented, using 35.9 fb −1 of data collected in pp collisions at √ s = 13 TeV.The cross section is first measured in a fiducial region of particle-level phase space by defining two categories of ttbb events, and subsequently this result is corrected to the total phase space.One of the defined fiducial regions corresponds to ignoring parton-level information, while the other uses parton-level information to identify the particle-level jets that do not originate from the decay of top quarks.For both definitions, the cross section is measured to be 1.6 ± 0.1 (stat) +0.5 −0.4 (syst) pb.The cross section in the total phase space is obtained by correcting this measurement for the experimental acceptance on the jets originating from the top quarks, which yields 5.5 ± 0.3 (stat) +1.6 −1.3 (syst) pb.This measurement provides valuable input to studies of the ttH process, where the Higgs boson decays into a pair of b quarks, and for which the normalization and modelling of the ttbb process represent a leading source of systematic uncertainty.Furthermore, these results represent a stringent test of perturbative quantum chromodynamics at the LHC.Predictions from several generators are compared with measurements and found to be smaller than the mea-Table 2: Measured and predicted cross sections for the different definitions of the ttbb phase space considered in this analysis.For measurements, the first uncertainty is statistical, while the second one is from the systematic sources.The uncertainties in the predicted cross sections include the statistical uncertainty, the PDF uncertainties, and the µ R and µ F dependences on changes in scale.The uncertainties in scale for parton showers are not included, and amount to about 15% for POWHEG+PYTHIA.Unless specified otherwise, PYTHIA is used for the modelling the parton shower, hadronization, and the underlying event.sured values by a factor of 1.5-2.4,corresponding to 1-2 standard deviations.This is consistent with previous results for the ttbb cross section and calls for further experimental and theoretical studies of the associated production of top quark pairs and b jets.

Figure 1 :
Figure1: Distributions in the QGLR (left) and the CWoLa BDT discriminants (right).Both are after preselection, requiring P(χ 2 ) > 10 −6 and at least eight selected jets.All the contributions are based on simulation.The multijet contribution is scaled to match the total yields in data, after the other processes including the ttbb signal have been normalized to their corresponding theoretical cross sections.This choice takes into account only the effect of the shape variation from the multijet background.The small backgrounds include ttV, ttH, single top quark, V+jets, and diboson production.The lower panels show the ratio between the observed data and the predictions.The dashed lines indicate the boundaries between the signal and control regions defined in Section 6. Hatched bands indicate the statistical uncertainty in the predictions without considering the systematic sources, dominated by the uncertainties in the simulated multijet background.Underflow and overflow events were added to the first and last bins, respectively.

Figure 2 :
Figure2: Distribution in the 2DCSV in the SR (upper left), CR1 (upper right), CR2 (lower right), and CR3 (lower left) regions.For clarity, the two-dimensional distribution with largest and next-to-largest b tagging discriminant scores for the additional jets have been unrolled to one dimension, and the resulting bins ordered according to increasing values of the ratio between expected signal and background yields in each bin of the SR.The small backgrounds include ttV, ttH, single top quark, V+jets, and diboson production.Hatched bands correspond to uncertainties.The bottom panels show the pull distribution.The pull is defined as the bin by bin difference between data and predicted yields after the fit, divided by the uncertainties accounted for correlations between data and predictions after the fit.

Figure 3 :
Figure3: Comparison of the measured ttbb production cross sections (vertical lines) with predictions from several Monte Carlo generators (squares), for three definitions of our ttbb regions of phase space: fiducial parton-independent (left), fiducial parton-based (middle), total (right).The dark (light) shaded bands show the statistical (total) uncertainties in the measured value.Uncertainty intervals in the theoretical cross sections include the statistical uncertainty as well as the uncertainties in the PDFs and the µ R and µ F scales.