Search for supersymmetry in final states with missing transverse momentum and multiple $b$-jets in proton-proton collisions at $\sqrt{s} = 13$ TeV with the ATLAS detector

A search for supersymmetry involving the pair production of gluinos decaying via third-generation squarks into the lightest neutralino ($\displaystyle\tilde\chi^0_1$) is reported. It uses LHC proton-proton collision data at a centre-of-mass energy $\sqrt{s} = 13$ TeV with an integrated luminosity of 36.1 fb$^{-1}$ collected with the ATLAS detector in 2015 and 2016. The search is performed in events containing large missing transverse momentum and several energetic jets, at least three of which must be identified as originating from $b$-quarks. To increase the sensitivity, the sample is divided into subsamples based on the presence or absence of electrons or muons. No excess is found above the predicted background. For $\displaystyle\tilde\chi^0_1$ masses below approximately 300 GeV, gluino masses of less than 1.97 (1.92) TeV are excluded at 95% confidence level in simplified models involving the pair production of gluinos that decay via top (bottom) squarks. An interpretation of the limits in terms of the branching ratios of the gluinos into third-generation squarks is also provided. These results improve upon the exclusion limits obtained with the 3.2 fb$^{-1}$ of data collected in 2015.


SUSY signal models
Various simplified SUSY models [17,18] are employed to optimise the event selection and/or interpret the results of the search. In terms of experimental signature, they all contain at least four b-jets originating from either gluino or top quark decays, and twoχ 0 1 , which escape the detector unseen, resulting in high E miss T . Gluinos are assumed to be pair-produced and to decay either asg →b 1b org →t 1t (the charge conjugate process is implied throughout this paper). The following top and bottom squark decays are then considered:t 1 → tχ 0 1 ,t 1 → bχ + 1 andb 1 → bχ 0 1 . 2 In all cases, the top or bottom squarks are assumed to be off-shell in order to have simplified models with only two parameters: the gluino andχ 0 1 masses. 3 All other sparticles are decoupled.
Two simplified models are used to optimise the event selection and to interpret the results. In the Gbb (Gtt) model, illustrated in Figure 1(a) (1(b)), each gluino undergoes an effective three-body decayg → bbχ 0 1 (g → ttχ 0 1 ) via off-shell bottom (top) squarks, with a branching ratio of 100%. The Gbb model is the simplest in terms of particle multiplicity, resulting in the minimal common features of four b-jets and twoχ 0 1 . In addition to these particles, the Gtt model produces four W bosons originating from the top quark decays: t → Wb. The presence of these four W bosons motivates the design of signal regions with a higher jet multiplicity than for Gbb models, and in some cases with at least one isolated electron or muon.g This paper includes an interpretation that probes the sensitivity of the search as a function of the gluino branching ratio, in addition to the gluino andχ 0 1 masses. Similar interpretations have been performed by the CMS collaboration [24,27]. For that interpretation a third gluino decay is considered:g → tbχ − 1 (via the off-shell top squark decayt * 1 →bχ − 1 ). Theχ − 1 is then forced to decay asχ ± 1 → W * χ 0 1 → ff χ 0 1 (where f denotes a fermion). To keep the number of model parameters at only two, the mass difference between theχ ± 1 and theχ 0 1 is fixed to 2 GeV. Such a small mass-splitting between theχ ± 1 and theχ 0 1 is typical of models where theχ 0 1 is dominated by the higgsinos, the superpartners of the neutral Higgs boson. Such models are well motivated by naturalness. The products of the decay W * → ff are typically too soft to be detected, except for very large mass differences between the gluino and theχ is not 100%, and where one, two or three top quarks, and thus on-shell W bosons, are possible in the final state, in between the Gbb (no top quarks) and Gtt (four top quarks) decay topologies. The decay topologies that are considered in the variable branching ratio model are illustrated in Figure 2. The model also includes the Gbb and Gtt decay topologies illustrated in 2 The decayb 1 → tχ − 1 is also possible but, followingg →b 1b , it yields the same final state asg →t * 1 t → (bχ The technical implementation of the simulated samples produced from these models is described in Section 4.

ATLAS detector
The ATLAS detector is a multipurpose particle physics detector with a forward-backward symmetric cylindrical geometry and nearly 4π coverage in solid angle. 4 The inner tracking detector (ID) consists 4 ATLAS uses a right-handed coordinate system with its origin at the nominal interaction point in the centre of the detector.
The positive x-axis is defined by the direction from the interaction point to the centre of the LHC ring, with the positive y-axis pointing upwards, while the beam direction defines the z-axis. Cylindrical coordinates (r, φ) are used in the transverse plane, φ being the azimuthal angle around the z-axis. The pseudorapidity η is defined in terms of the polar angle θ by η = − ln tan(θ/2). Rapidity is defined as y = 0.5 ln[(E + p z )/(E − p z )] where E denotes the energy and p z is the component of the momentum along the beam direction.
of silicon pixel and microstrip detectors covering the pseudorapidity region |η| < 2.5, surrounded by a transition radiation tracker, which enhances electron identification in the region |η| < 2.0. Before the start of Run 2, the new innermost pixel layer, the insertable B-layer (IBL) [28], was inserted at a mean sensor radius of 3.3 cm. The ID is surrounded by a thin superconducting solenoid providing an axial 2 T magnetic field and by a fine-granularity lead/liquid-argon (LAr) electromagnetic calorimeter covering |η| < 3.2. A steel/scintillator-tile calorimeter provides coverage for hadronic showers in the central pseudorapidity range (|η| < 1.7). The endcaps (1.5 < |η| < 3.2) of the hadronic calorimeter are made of LAr active layers with either copper or tungsten as the absorber material. The forward region (3.1 < |η| < 4.9) is instrumented with a LAr calorimeter for both the EM and hadronic measurements. A muon spectrometer with an air-core toroidal magnet system surrounds the calorimeters. Three layers of high-precision tracking chambers provide coverage in the range |η| < 2.7, while dedicated fast chambers allow triggering in the region |η| < 2.4. The ATLAS trigger system [29] consists of a hardware-based level-1 trigger followed by a software-based high-level trigger (HLT).

Data and simulated event samples
The data used in this analysis were collected by the ATLAS detector from pp collisions produced by the LHC at a centre-of-mass-energy of 13 TeV and 25 ns proton bunch spacing over the 2015 and 2016 datataking periods. The full dataset corresponds to an integrated luminosity of 36.1 fb −1 after the application of beam, detector and data-quality requirements. The uncertainty in the combined 2015+2016 integrated luminosity is 2.1%. It is derived, following a methodology similar to that detailed in Ref. [30], from a preliminary calibration of the luminosity scale using x-y beam-separation scans performed in August 2015 and May 2016. Events are required to pass an E miss T trigger with thresholds of 70 GeV, 100 GeV and 110 GeV at the HLT level for the 2015, early 2016 and late 2016 datasets, respectively. These triggers are fully efficient for events passing the preselection defined in Section 6, which requires the offline reconstructed E miss T to exceed 200 GeV. There are on average 24 inelastic pp collisions (the interactions other than the hard scatter are referred to as "pile-up") in the dataset.
Samples of Monte Carlo (MC) simulated events are used to model the signal and background processes in this analysis, except multijet processes, which are estimated from data. SUSY signal samples in which each gluino decays into bbχ 0 1 , ttχ 0 1 , or tbχ − 1 were generated with up to two additional partons using MadGraph5_aMC@NLO [31] v2.2.2 at leading order (LO) with the NNPDF 2.3 [32] parton distribution function (PDF) set. These samples were interfaced to Pythia v8.186 [33] for the modelling of the parton showering, hadronisation and underlying event.
The dominant background in the signal regions is the production of tt pairs with additional high transverse momentum (p T ) jets. For the generation of tt and single top quarks in the Wt-and s-channels the Powheg-Box [34] v2 event generator with the CT10 [35] PDF set in the matrix element calculations was used. Electroweak t-channel single-top-quark events were generated using the Powheg-Box v1 event generator. This event generator uses the four-flavour scheme for the next-to-leading order (NLO) matrix elements calculations together with the fixed four-flavour PDF set CT10f4. For all processes involving top quarks, top-quark spin correlations are preserved. In the t-channel, top quarks were decayed using MadSpin [36]. The parton shower, fragmentation, and the underlying event were simulated using Pythia v6.428 [37] with the CTEQ6L1 PDF set [38]. The h damp parameter in Powheg, which controls the p T of the first additional emission beyond the Born level and thus regulates the p T of the recoil emission against the tt system, was set to the mass of the top quark (m top = 172.5 GeV). All events with at least one leptonically decaying W boson are included. Single-top and tt events in which all top quarks decay hadronically do not contain sufficient E miss T to contribute significantly to the background.
Smaller backgrounds in the signal region come from the production of tt pairs in association with W/Z/h bosons and possibly additional jets, and production of tttt, W/Z+jets and WW/WZ/ZZ (diboson) events. Other potential sources of background, such as the production of three top quarks or three gauge bosons, are expected to be negligible. The production of tt pairs in association with electroweak vector bosons W and Z was modelled by samples generated at LO using MadGraph5_aMC@NLO v2.2.2 and showered with Pythia v8.186, while samples to model ttH production were generated using Mad-Graph5_aMC@NLO v2.2.1 and showered with Herwig++ [39] v2.7.1. These samples are described in detail in Ref. [40]. MadGraph5_aMC@NLO was also used to simulate the tttt production and the showering was performed with Pythia v8.186. The W/Z+jets processes were simulated using the Sherpa v2.2.0 [41] event generator, while Sherpa v2.1.1 was used to simulate diboson production processes. Matrix elements for the W/Z+jets and diboson processes were calculated using Comix [42] and Open-Loops [43] and merged with the Sherpa parton shower [44] using the ME+PS@NLO prescription [45]. The Sherpa diboson sample cross-section was scaled down to account for its use of α QED = 1/129 rather than 1/132, corresponding to the use of current Particle Data Group [46] parameters, as input to the G µ scheme [47]. Samples generated using MadGraph5_aMC@NLO v2.2.2 were produced with the NNPDF 2.3 PDF set and W/Z+jets samples were generated with the NNPDF 3.0 PDF set [48], while all other samples used CT10 PDFs.
For all samples, except the ones generated using Sherpa, the EvtGen v1.2.0 program [49] was used to simulate the properties of the bottom-and charm-hadron decays. All Pythia v6.428 samples used the PERUGIA2012 [50] set of tuned parameters (tune) for the underlying event, while Pythia v8.186 and Herwig++ showering were run with the A14 [51] and UEEE5 [52] underlying-event tunes, respectively. In-time and out-of-time pile-up interactions from the same or nearby bunch-crossings were simulated by overlaying additional pp collisions generated by Pythia v8.186 using the A2 tune [53] and the MSTW2008LO parton distribution function set [54] on top of the hard-scattering events. Details of the sample generation and normalisation are summarised in Table 1. Additional samples with different event generators and settings are used to estimate systematic uncertainties in the backgrounds, as described in Section 7.
All simulated event samples were passed through the full ATLAS detector simulation using Geant4 [55], with the exception of signal samples in which at least one gluino decays asg → bbχ 0 1 org → tbχ − 1 , which were passed through a fast simulation that uses a parameterisation for the calorimeter response [56] and Geant4 for the ID and the muon spectrometer. The simulated events are reconstructed with the same algorithm as that used for data.
The signal samples are normalised using the best cross-section calculations at NLO in the strong coupling constant, adding the resummation of soft gluon emission at next-to-leading-logarithm (NLL) accuracy [57][58][59][60][61]. The nominal cross-section and the uncertainty are taken from an envelope of cross-section predictions using different PDF sets and factorisation and renormalisation scales, as described in Ref. [62]. The cross-section of gluino pair-production in these simplified models is 14 ± 3 fb for a gluino mass of 1.5 TeV, falling to 1.0 ± 0.3 fb for 2 TeV mass gluinos. All background processes are normalised using the best available theoretical calculation for their respective cross-sections. The order of this calculation in perturbative QCD (pQCD) for each process is listed in Table 1. For tt, the largest background, this corresponds to a cross-section of 831.8 pb. Finally, contributions from multijet background are estimated from data using a procedure described in Ref. [63], which performs a smearing of the jet response in data events with well-measured E miss T (socalled "seed events"). The response function is derived in Monte Carlo dijet events and is different for b-tagged and non-b-tagged jets.

Event reconstruction
Interaction vertices from the proton-proton collisions are reconstructed from at least two tracks with p T > 0.4 GeV, and are required to be consistent with the beamspot envelope. The primary vertex is identified as the one with the largest sum of squares of the transverse momenta from associated tracks ( |p T,track | 2 ) [71].
Basic selection criteria are applied to define candidates for electrons, muons and jets in the event. An overlap removal procedure is applied to these candidates to prevent double-counting. Further requirements are then made to select the final signal leptons and jets from the remaining candidates. The details of the candidate selections and of the overlap removal procedure are given below.
Candidate jets are reconstructed from three-dimensional topological energy clusters [72] in the calorimeter using the anti-k t jet algorithm [73,74] with a radius parameter of 0.4 (small-R jets). Each topological cluster is calibrated to the electromagnetic scale response prior to jet reconstruction. The reconstructed jets are then calibrated to the particle level by the application of a jet energy scale (JES) derived from √ s = 13 TeV data and simulations [75]. Quality criteria are imposed to reject events that contain at least one jet arising from non-collision sources or detector noise [76]. Further selections are applied to reject jets that originate from pile-up interactions by means of a multivariate algorithm using information about the tracks matched to each jet [77]. Candidate jets are required to have p T > 20 GeV and |η| < 2.8. After resolving overlaps with electrons and muons, selected jets are required to satisfy the stricter requirement of p T > 30 GeV.
A jet is tagged as a b-jet candidate by means of a multivariate algorithm using information about the impact parameters of inner detector tracks matched to the jet, the presence of displaced secondary vertices, and the reconstructed flight paths of b-and c-hadrons inside the jet [78,79]. The b-tagging working point corresponding to an efficiency of 77% to identify b-jets with p T > 20 GeV, as determined from a sample of simulated tt events, is found to be optimal for the statistical significance of this search. The corresponding rejection factors against jets originating from c-quarks, τ-leptons and light quarks and gluons in the same sample at this working point are 6, 22 and 134, respectively.
After resolving the overlap with leptons, the candidate small-R jets are re-clustered [80] into large-R jets using the anti-k t algorithm with a radius parameter of 0.8. The calibration from the input small-R jets propagates directly to the re-clustered jets. These re-clustered jets are then trimmed [80][81][82][83] by removing subjets whose p T falls below 10% of the p T of the original re-clustered jet. The resulting large-R jets are required to have p T > 100 GeV and |η| < 2.0. When it is not explicitly stated otherwise, the term "jets" in this paper refers to small-R jets.
Electron candidates are reconstructed from energy clusters in the electromagnetic calorimeter and inner detector tracks and are required to satisfy a set of "loose" quality criteria [84,85]. They are also required to have |η| < 2.47. Muon candidates are reconstructed from matching tracks in the inner detector and muon spectrometer. They are required to meet "medium" quality criteria, as described in Ref.
Leptons are selected from the candidates that survive the overlap removal procedure if they fulfil a requirement on the scalar sum of p T of additional inner detector tracks in a cone around the lepton track. This isolation requirement is defined to ensure a flat efficiency of around 99% across the whole electron transverse energy and muon transverse momentum ranges. The angular separation between the lepton and the b-jet ensuing from a semileptonic top quark decay narrows as the p T of the top quark increases. This increased collimation is accounted for by setting the radius of the isolation cone to min(0.2, 10 GeV/p lep T ), where p lep T is the lepton p T expressed in GeV. Selected electrons are further required to meet the "tight" quality criteria [84,85]. Electrons (muons) are matched to the primary vertex by requiring the transverse impact parameter d 0 of the associated ID track to satisfy |d 0 |/σ d 0 < 5 (3), where σ d 0 is the measured uncertainty of d 0 , and the longitudinal impact parameter z 0 to satisfy |z 0 sin θ| < 0.5 mm. 5 In addition, events containing one or more muon candidates with |d 0 | (|z 0 |) > 0.2 mm (1 mm) are rejected to suppress cosmic rays.
Overlaps between candidate objects are removed sequentially. Firstly, electron candidates that lie a distance 6 ∆R < 0.01 from muon candidates are removed to suppress contributions from muon bremsstrahlung. Overlaps between electron and jet candidates are resolved next, and finally, overlaps between remaining jets and muon candidates are removed.
Overlap removal between electron and jet candidates aims to resolve two sources of ambiguity: it is designed, firstly, to remove jets that are formed primarily from the showering of a prompt electron and, secondly, to remove electrons that are produced in the decay chains of hadrons. Consequently, any nonb-tagged jet whose axis lies ∆R < 0.2 from an electron is discarded. Electrons with E T < 50 GeV are discarded if they lie ∆R < 0.4 from the axis of any remaining jet and the corresponding jet is kept. For higher-E T electrons, the latter removal is performed using a threshold of ∆R = min(0.4, 0.04+10 GeV/E T ) to increase the acceptance for events with collimated top quark decays.
The procedure to remove overlaps between muon and jet candidates is designed to remove those muons that are likely to have originated from the decay of hadrons and to retain the overlapping jet. Jets and muons may also appear in close proximity when the jet results from high-p T muon bremsstrahlung, and in such cases the jet should be removed and the muon retained. Such jets are characterised by having very few matching inner detector tracks. Therefore, if the angular distance ∆R between a muon and a jet is lower than 0.2, the jet is removed if it is not b-tagged and has fewer than three matching inner detector tracks. Like the electrons, muons with p T below (above) 50 GeV are subsequently discarded if they lie within ∆R = 0.4 (∆R = min(0.4, 0.04 + 10 GeV/p T )) of any remaining jet.
The missing transverse momentum (E miss T ) in the event is defined as the magnitude of the negative vector sum ( p T miss ) of the transverse momenta of all selected and calibrated objects in the event, with an extra term added to account for energy deposits that are not associated with any of these selected objects. This "soft" term is calculated from inner detector tracks matched to the primary vertex to make it more resilient to contamination from pile-up interactions [87, 88] .
Corrections derived from data control samples are applied to simulated events to account for differences between data and simulation in the reconstruction efficiencies, momentum scale and resolution of leptons, in the efficiency and fake rate for identifying b-jets, and in the efficiency for rejecting jets originating from pile-up interactions.

Event selection
The event selection criteria are defined based on kinematic requirements for the objects defined in Section 5. Other discriminating event-based variables, described in Section 6.1, are used to further reject the background. Two sets of preselection criteria targeting the 0-lepton and the 1-lepton channels are presented in Section 6.2. The modelling of the data in these regions is also discussed in that section. The general analysis strategy and the treatment of background sources is presented in Section 6.3. Finally, the event selection for the cut-and-count and multi-bin analyses are discussed in Sections 6.4 and 6.5, respectively.

Discriminating variables
The effective mass variable (m eff ) is defined as: where the first and second sums are over the selected jets (N jet ) and leptons (N lepton ), respectively. It typically has a much higher value in pair-produced gluino events than in background events.
In regions with at least one selected lepton, the transverse mass m T composed of the p T of the leading selected lepton ( ) and E miss T is defined as: It is used to reduce the tt and W+jets background events in which a W boson decays leptonically. Neglecting resolution effects, the m T distribution for these backgrounds has an expected upper bound corresponding to the W boson mass and typically has higher values for Gtt events. Another useful transverse mass variable is m b-jets T,min , the minimum transverse mass formed by E miss T and any of the three highest-p T b-tagged jets in the event: The m b-jets T,min distribution has an expected upper bound corresponding to the top quark mass for tt events with a semileptonic top quark decay, while peaking at higher values for Gbb and Gtt events.
Another powerful variable is the total jet mass variable, defined as: where m J,i is the mass of the large-radius re-clustered jet i in the event. The decay products of a hadronically decaying boosted top quark can be reconstructed in a single large-radius re-clustered jet, resulting in a jet with a high mass. This variable typically has larger values for Gtt events than for background events. This is because Gtt events contain as many as four hadronically decaying top quarks while the background is dominated by tt events with one or two semileptonic top quark decays.
The requirement of a selected lepton, with the additional requirements on jets, E miss T and event variables described above, makes the multijet background negligible for the ≥ 1-lepton signal regions. For the 0-lepton signal regions, the minimum azimuthal angle ∆φ 4j min between p miss T and the p T of the four leading small-R jets in the event, defined as: is required to be greater than 0.4. This requirement supresses the multijet background, which can produce events with large E miss T if containing poorly measured jets or neutrinos emitted close to the axis of a jet. A similar variable, denoted ∆φ j 1 , is also used in the Gbb signal regions targeting small mass differences between the gluino and the neutralino, allowing the identification of events containing a high-p T jet coming from initial-state radiation (ISR) and recoiling against the gluino pair. It is defined as the absolute value of the azimuthal angle separating the p T of the leading jet and p miss T , and is expected to have larger values for the targeted signal than for the background.

Modelling of the data
Preselection criteria in the 0-lepton and 1-lepton channels require E miss T > 200 GeV, in addition to the E miss T trigger requirement, and at least four jets of which at least two must be b-tagged. The 0-lepton (1-lepton) channel requires the event to contain no (at least one) selected lepton.
In this analysis, correction factors need to be extracted to account for shape discrepancies in the m eff spectrum between the data and the expected background for the 1-lepton preselection sample. These factors are defined as the ratio of the number of observed events to the predicted number of background events in a given m eff bin, in a signal-depleted region. This region is defined by applying the 1-lepton preselection criteria and requiring exactly two b-tagged jets and m b-jets T,min < 140 GeV. This kinematic reweighting leads to correction factors ranging from 0.7 to 1.1. They are applied to the background prediction and the full size of the correction is taken as an uncertainty for both the background and signal events. Figures 3 and 4 show the multiplicity of selected jets and b-tagged jets, the distributions of E miss T , m eff , and M Σ J for events passing the 0-lepton or the 1-lepton preselection, respectively. Figure 3 (4) also displays the distribution of m b-jets T,min (m T ) in the 0-lepton (1-lepton) channel. The correction described above is applied in the 1-lepton channel. The uncertainty bands include the statistical and experimental systematic uncertainties, as described in Section 7, but not the theoretical uncertainties in the background modelling.
The data and the predicted background are found to agree reasonably well at the preselection level after the kinematic reweighting described above. A discrepancy between data and prediction is observed for the number of b-tagged jets, but it has a negligible impact on the background estimate after the renormalisation of the simulation in dedicated control regions with the same b-tagged jets requirements as the signal regions, as described in Sections 6.4 and 6.5. Example signal models with enhanced cross-sections are overlaid for comparison.

Analysis strategy and background treatment
In order to enhance the sensitivity to the various signal benchmarks described in Section 2, multiple signal regions (SRs) are defined. The main background in all these regions is the production of a tt pair in association with heavy-and light-flavour jets. A normalisation factor for this background is extracted for each individual SR from a data control region (CR) that has comparable background composition and kinematics. This is ensured by keeping the kinematic requirements similar in the two regions. The CRs and SRs are defined to be mutually exclusive. Signal contributions in the CRs are suppressed by inverting or relaxing some requirements on the kinematic variables (e.g. m T or m b-jets T,min ), leading to a signal contamination in the CRs of 6% at most. The tt normalisation is cross-checked in validation regions (VRs) that share similar background composition, i.e. jet and lepton flavours, with the SR. The signal contamination in the VRs is found to be lower than 30% for benchmark signal mass points above the already excluded mass range. The tt purity is superior to 73% and 53% in the CRs and VRs, respectively.
The non-tt backgrounds mainly consist of single-top, W+jets, Z+jets, tt +W/Z/h, tttt and diboson events. Their normalisation is taken from the simulation normalised using the best available theory prediction. The multijet background is found to be very small or negligible in all regions. It is estimated using a procedure described in Ref. [63], in which the jet response is determined from simulated dijet events. This response function is then used to smear the jet response in low-E miss T events. The jet response is cross-checked with data where the E miss T can be unambiguously attributed to the mismeasurement of one of the jets. Two analysis strategies are followed, and different SR sets are defined for each: • A cut-and-count analysis, using partially overlapping single-bin SRs, optimised to maximise the expected discovery power for benchmark signal models, and allowing for reinterpretation of the results. The SRs are defined to probe the existence of a signal or to assess model-independent upper limits on the number of signal events. T,min for events passing the 0lepton preselection criteria. The statistical and experimental systematic uncertainties (as defined in Section 7) are included in the uncertainty band. The last bin includes overflow events. The lower part of each figure shows the ratio of data to the background prediction. All backgrounds (including tt) are normalised using the best available theoretical calculation described in Section 4. The background category tt + X includes ttW/Z, tth and tttt events. Example signal models with cross-sections enhanced by a factor of 50 are overlaid for comparison. The last bin includes overflow events. The lower part of each figure shows the ratio of data to the background prediction. All backgrounds (including tt) are normalised using the best available theoretical calculation described in Section 4. The background category tt + X includes ttW/Z, tth and tttt events. Example signal models with cross-sections enhanced by a factor of 50 are overlaid for comparison.
• A multi-bin analysis, using a set of non-overlapping SRs and CRs that are combined to strengthen the exclusion limits on the targeted signal benchmarks. This set of regions is used to assess modeldependent interpretations of the various signal models.

Cut-and-count analysis
The SRs are named in the form SR-X-YL-Z, where X indicates the target model, Y indicates the number of leptons and Z labels the type of region targeted. The cut-and-count regions labelled B (for "boosted") are optimised for signals with a large mass difference between the gluino and the neutralino (∆m 1.5 TeV), possibly leading to highly boosted objects in the final state. Conversely, regions C (for "compressed") primarily focus on signals for which the gluino decay products are softer due to the small ∆m (∆m 300 GeV). Regions M (for "moderate") target intermediate values of ∆m. SRs targeting the Gtt model in the 1-and 0-lepton channels are presented in Table 2.
In the 1-lepton channel, these regions differ mainly in their kinematic selections thresholds: m eff , E miss The signal regions of the 0-lepton channel follow a similar strategy to the 1-lepton channel. Background composition studies performed on simulated event samples show that semileptonic tt events, for which the lepton is outside the acceptance or is a hadronically decaying τ-lepton, dominate in the SRs. Thus, CRs to normalise the tt+jets background make use of the 1-lepton channel, requiring the presence of exactly one signal lepton. An inverted selection on m T is applied to remove overlaps with the 1-lepton SRs. The background prediction is validated in a 0-lepton region, inverting the M Σ J selection to remove any overlap with the SRs.
Regions targeting the Gbb model are presented in Table 3. The region definition follows the same pattern as for Gtt-0L regions, in particular for regions B, M and C. For very small values of ∆m, the Gbb signal does not lead to a significant amount of E miss T , except if a hard ISR jet recoils against the gluino pair. Such events are targeted by region VC (for "very compressed") that identifies an ISR-jet candidate as a non-b-tagged high-p T leading jet (j 1 ), with a large azimuthal separation ∆φ j 1 with respect to p T miss .
Similarly, the normalisation factor of the tt background is extracted from a 1-lepton CR, to which an inverted selection on m T is applied to remove the overlaps with Gtt 1-lepton SRs and the corresponding signal contamination. The 0-lepton VRs are constructed in the 0-lepton channel with selections very close to the SR ones. They are mutually exclusive due to an inverted E miss T selection in the VR. Table 2: Definitions of the Gtt SRs, CRs and VRs of the cut-and-count analysis. All kinematic variables are expressed in GeV except ∆φ 4j min , which is in radians. The jet p T requirement is also applied to b-tagged jets.

Gtt 1-lepton
Criteria common to all regions: Region M (Moderate ∆m) Region C (Compressed, moderate ∆m)   Figures 3 and 4 show that a good separation between signal and background can be achieved with various kinematic variables. The distribution of N jet and m eff for different signal benchmarks and ∆m values is used to build a two-dimensional slicing of the phase space in a set of non-overlapping SRs, CRs and VRs that can be statistically combined. The slicing scheme is presented in Figure 5. The SRs are named in the form SR-YL-Z 1 Z 2 , where Y indicates the number of leptons, Z 1 labels the jet multiplicity bin and Z 2 labels the m eff bin. For Z 1 and Z 2 , the letters "H" stands for "high", "I" for "intermediate" and "L" for "low". In the 0-lepton channel, there is also a 0L-ISR region that is a subset of the IL, LL, II and LI regions, and kept mutually exclusive with them as detailed below.

Multi-bin analysis
The low-N jet region probes especially Gbb-like models, for which the number of hard jets is lower than in decay topologies containing top quarks. This category of events is thus only considered in the 0-lepton channel. Gtt events are mostly expected in the high-N jet bin. The intermediate jet multiplicity bin is built to be sensitive to decay topologies with a number of top quarks intermediate between Gbb and Gtt, but also to Gbb (with additional jets originating from radiation) and to Gtt (when some jets fall outside the acceptance). The m eff bins are chosen to provide sensitivity to various kinematic regimes: the low-m eff regions are essentially sensitive to soft signals (low ∆m), while the high-m eff regions are designed to select highly boosted events.
For each N jet -m eff region presented in Figure 5, the selection was optimised over all the other variables to maximise the exclusion power for the Gbb and Gtt models. For each m eff bin, a targeted range of ∆m was used in the optimisation procedure. with the 1-lepton SRs. The other kinematic requirements are kept close to the ones of the SR. One VR is defined for each SR in the corresponding lepton channel. Full independence between the signal and VRs is guaranteed by E miss T and m b-jets T,min requirements. The low-N jet regions are presented in Table 6. Targeting primarily the Gbb model, the transverse momentum of the fourth jet is required to be larger than 90 GeV in all SRs. In the intermediate and low m eff regions, the leading jet is required to be b-tagged or the value of ∆φ j 1 to be lower than 2.9 in order to be mutually exclusive with the 0L-ISR regions. The tt background dominates in all regions, and is normalised in dedicated 1-lepton regions, defined with a low m T requirement, as done for the regions of the cut-and-count analysis. VRs are constructed in the 0-lepton channel, closely reproducing the background composition and kinematics of the SR events.
A dedicated set of regions is designed to target very compressed Gbb scenarios in which a hard ISR jet recoils against the gluino pair. The definition of these regions is presented in Table 6.

High-N jet regions
Criteria common to all regions:

Systematic uncertainties
Figures 6(a) and 6(b) summarise the relative systematic uncertainties in the background estimate for the cut-and-count and multi-bin analyses, respectively. These uncertainties arise from the extrapolation of the tt normalisation obtained in the CRs to the SRs as well as from the yields of the minor backgrounds in the SRs, which are predicted by the simulation. The total systematic uncertainties range from approximately 20% to 80% in the various SRs.
The detector-related systematic uncertainties affect both the background estimate and the signal yield. The largest sources in this analysis relate to the jet energy scale, jet energy resolution (JER) and the b-tagging efficiencies and mistagging rates. The JES uncertainties for the small-R jets are derived from √ s = 13 TeV data and simulations while the JER uncertainties are extrapolated from 8 TeV data using MC simulations [89]. These uncertainties are also propagated to the re-clustered large-R jets, which use them as inputs. The jet mass scale and resolution uncertainties have a negligible impact on the re-clustered jet mass. The impact of the JES uncertainties on the expected background yields is between 4% and 35%, while JER uncertainties affect the background yields by approximately 0-26% in the various regions. Uncertainties in the measured b-tagging efficiencies and mistagging rates are the subleading sources of experimental uncertainty.
The impact of these uncertainties on the expected background yields is 3-24% depending on the considered region. The uncertainties associated with lepton reconstruction and energy measurements have a negligible impact on the final results. All lepton and jet measurement uncertainties are propagated to the calculation of E miss T , and additional uncertainties are included in the scale and resolution of the soft term. The overall impact of the E miss T soft-term uncertainties is also small.
Since the normalisation of the tt background is fit to data in the CRs, uncertainties in the modelling of this background only affect the extrapolation from the CRs to the SRs and VRs. Hadronisation and parton showering model uncertainties are estimated using a sample generated with Powheg and showered by Herwig++ v2.7.1 with the UEEE5 underlying-event tune. Systematic uncertainties in the modelling of initial-and final-state radiation are explored with Powheg samples showered with two alternative settings of Pythia v6.428. The first of these uses the PERUGIA2012radHi tune [50] and has the renormalisation and factorisation scales set to twice the nominal value, resulting in more radiation in the final state. In addition, it has h damp set to 2m top . The second sample, using the PERUGIA2012radLo tune, has h damp = m top and the renormalisation and factorisation scales are set to half of their nominal values, resulting in less radiation in the event. In each case, the uncertainty is taken as the change in the expected yield of tt background with respect to the nominal sample. The uncertainty due to the choice of event generator is estimated by comparing the expected yields obtained using a tt sample generated with MadGraph5_aMC@NLO and one that is generated with Powheg. Both of these samples are showered with Herwig++ v2.7.1. The total theoretical uncertainty in the inclusive tt background is taken as the sum in quadrature of these individual components. An additional uncertainty is assigned to the fraction of tt events produced in association with additional heavy-flavour jets (i.e. tt+ ≥ 1b and tt+ ≥ 1c), a process which suffers from large theoretical uncertainties. Simulation studies show that the heavy-flavour fractions in each set of SR, CR and VR, which have almost identical b-tagged jets requirements, are similar. Therefore, the theoretical uncertainties in this fraction affect these regions in a similar way, and thus largely cancel out in the semi-data-driven tt normalisation based on the observed CR yields. The residual uncertainty in the tt prediction is taken as the difference between the nominal tt prediction and the one obtained after varying the cross-section of tt events with additional heavy-flavour jets by 30%, in accordance with the results of the ATLAS measurement of this cross-section at √ s = 8 TeV [90]. This component typically makes a small contribution (0-8%) to the total impact of the tt modelling uncertainties on the background yields, which ranges between 5% and 76% for the various regions. The statistical uncertainty of the CRs used to extract the tt normalisation factors, which is included in the systematic uncertainties, ranges from 10% to 30% depending on the SR.
Modelling uncertainties affecting the single-top process arise especially from the interference between the tt and Wt processes. This uncertainty is estimated using inclusive WWbb events, generated using MadGraph5_aMC@NLO, which are compared with the sum of tt and Wt processes. Furthermore, as in the tt modelling uncertainties, variations of Pythia v6.428 settings increasing or decreasing the amount of radiation are also used. An additional 5% uncertainty is included in the cross-section of single-top processes [91]. Overall, the modelling uncertainties affecting the single-top process lead to changes of approximately 0-11% in the total yields in the various regions. Uncertainties in the W/Z+jets backgrounds are estimated by varying independently the scales for factorisation, renormalisation and resummation by factors of 0.5 and 2. The scale used for the matching between jets originating from the matrix element and the parton shower is also varied. The resulting uncertainties in the total yield range from approximately 0 to 50% in the various regions. A 50% normalisation uncertainty is assigned to tt + W/Z/h, tttt and diboson backgrounds and are found to have no significant impact on the sensitivity of this analysis. Uncertainties arising from variations of the parton distribution functions were found to affect background yields by less than 2%, and therefore these uncertainties are neglected here. Uncertainties due to the limited number of events in the MC background samples are included if above 5%. They reach approximately 20% in regions targeting large mass-splitting.
The uncertainties in the cross-sections of signal processes are determined from an envelope of different cross-section predictions, as described in Section 4. A systematic uncertainty is also assigned to the kinematic correction described in Section 6. The total size of the correction is used as an uncertainty, and is applied to all simulated event samples for the 1-lepton channel.

Results
The expected SM background is determined separately in each SR with a profile likelihood fit [92] implemented in the HistFitter framework [93], referred to as a background-only fit. The fit uses as a constraint the observed event yield in the associated CR to adjust the tt normalisation, assuming that no signal contributes to this yield, and applies that normalisation factor to the number of tt events predicted by simulation in the SR. The values of the normalisation factors, the expected numbers of background events and the observed data yields in all the CRs are shown in Figures 7(a) and 7(b) for the cut-and-count and multi-bin analyses, respectively.
The inputs to the background-only fit for each SR are the number of events observed in its associated CR and the number of events predicted by simulation in each region for all background processes. The numbers of observed and predicted events in each CR are described by Poisson probability density functions. The systematic uncertainties in the expected values are included in the fit as nuisance parameters. They are constrained by Gaussian distributions with widths corresponding to the sizes of the uncertainties and are treated as correlated, when appropriate, between the various regions. The product of the various probability density functions forms the likelihood, which the fit maximises by adjusting the tt normalisation and the nuisance parameters.
C R _ IS tR _ G tt C R _ L n j_ L m e ff _ G tt C R _ L n j_ Im e ff _ G tt C R _ L n j_ H m e ff _ G tt C R _ In j_ L m e ff _ G tt C R _ In j_ Im e ff _ G tt C R _ H n j_ L m e ff _ G tt C R _ H n j_ Im e ff _ G tt C R _ H n j_ H m e ff _ G tt The background category tt + X includes ttW/Z, ttH and tttt events. All of these regions require at least one signal lepton, for which the multijet background is negligible. All uncertainties describes in Section 7 are included in the uncertainty band. The tt normalisation is obtained from the fit and is displayed in the bottom panel. Figures 8(a) and 8(b) show the results of the background-only fit to the CRs, extrapolated to the VRs for the cut-and-count and multi-bin analyses, respectively. The number of events predicted by the backgroundonly fit is compared to the data in the upper panel. The pull, defined by the difference between the observed number of events (n obs ) and the predicted background yield (n pred ) divided by the total uncertainty (σ tot ), is shown for each region in the lower panel. No evidence of significant background mismodelling is observed in the VRs.
The event yields in the SRs for the cut-and-count and multi-bin analyses are presented in Figure 9, where the pull is shown for each region in the lower panel. No significant excess is found above the predicted background. The maximum deviation is observed in region SR-0L-HH of the multi-bin analysis with a local significance of 2.3 standard deviations. The background is dominated by tt events in all SRs. The subdominant background contributions in the 0-lepton regions are Z(→ νν)+jets and W(→ ν)+jets events, where for W+jets events the lepton is an unidentified electron or muon or a hadronically decaying τ-lepton. In the 1-lepton SRs, the subdominant backgrounds are single-top, ttW and ttZ. Table 7 shows the observed number of events and predicted number of background events from the background-only fit in the Gtt 1-lepton, Gtt 0-lepton and Gbb regions for the cut-and-count analysis. The central value of the fitted background is in general larger than the MC-only prediction. This is in part due to an underestimation of the cross-section of tt+ ≥ 1b and tt+ ≥ 1c processes in the simulation.

Interpretation
Since no significant excess over the expected background from SM processes is observed, the data are used to derive one-sided upper limits at 95% confidence level (CL). Two levels of interpretation are provided in this paper: model-independent exclusion limits and model-dependent exclusion limits set on the Gbb, Gtt and gluino variable branching ratio models.

Model-independent exclusion limits
Model-independent limits on the number of beyond-the-SM (BSM) events for each SR are derived with pseudoexperiments using the CL s prescription [94] and neglecting a possible signal contamination in the CR. Only the single-bin regions from the cut-and-count analysis are used for this purpose, to aid in the reintepretation of these limits. Limits are obtained with a fit in each SR which proceeds in the same way as the fit used to predict the background, except that the number of events observed in the SR is added as an input to the fit. Also, an additional parameter for the non-SM signal strength, constrained to be non-negative, is fit. Upper limits on the visible BSM cross-section (σ 95 vis ) are obtained by dividing the observed upper limits on the number of BSM events with the integrated luminosity. The results are given in Table 8, where the p 0 -values, which represent the probability of the SM background alone to fluctuate to the observed number of events or higher, are also provided.

Model-dependent exclusion limits
The results are used to place exclusion limits on various signal models. The results are obtained using the CL s prescription in the asymptotic approximation [92]. The expected and observed limits were compared to the CL s calculated from pseudoexperiments and found to be compatible. The signal contamination   Figure 7. The upper panel shows the observed number of events and the predicted background yield. All uncertainties defined in Section 7 are included in the uncertainty band. The background category tt + X includes ttW/Z, ttH and tttt events. The lower panel shows the pulls in each VR.   Table 7: Results of the background-only fit extrapolated to the Gtt 1-lepton, Gtt 0-lepton and Gbb SRs in the cutand-count analysis, for the total background prediction and breakdown of the main background sources. The uncertainties shown include all systematic uncertainties. The data in the SRs are not included in the fit. The background category tt + X includes ttW/Z, ttH and tttt events. The row "MC-only background" provides the total background prediction when the tt normalisation is obtained from a theoretical calculation [64].

ATLAS
SR-Gtt-1L  in the CRs and the experimental systematic uncertainties in the signal are taken into account for this calculation. All the regions of the multi-bin analysis are statistically combined to set model-dependent upper limits on the Gbb, Gtt and variable branching ratio models.
The 95% CL observed and expected exclusion limits for the Gtt and Gbb models are shown in the LSP and gluino mass plane in Figures 10(a) and 10(b), respectively. The ±1σ SUS Y theory lines around the observed limits are obtained by changing the SUSY cross-section by one standard deviation (±1σ), as described in Section 4. The yellow band around the expected limit shows the ±1σ uncertainty, including all statistical and systematic uncertainties except the theoretical uncertainties in the SUSY cross-section. Compared to the previous results [19], the gluino mass sensitivities of the current search (assuming massless LSPs) have improved by 300 GeV and 450 GeV for the Gbb and Gtt models, respectively. Gluinos with masses below 1.97 (1.92) TeV are excluded at 95% CL for neutralino masses lower than 300 GeV in the Gtt (Gbb) model. The observed limit for the Gtt model at high gluino mass is weaker than the expected limits due to the mild excesses observed in the signal regions SR-0L-HH and SR-1L-HI of the multi-bin fit analysis. The best exclusion limit on the LSP mass is approximately 1.19 (1.20) TeV, which is reached for a gluino mass of approximately 1.40 (1.68) TeV for Gbb and Gtt models, respectively.
Limits are also set in the signal model described in Section 2 for which the branching ratios of the gluinos to tbχ Similar results are presented in Figure 11(b) assuming a gluino mass of 1.9 TeV and scanning various neutralino masses (1, 600 and 1000 GeV). For neutralino masses between 1 and 600 GeV, most of the branching ratio plane is expected to be excluded at 95% CL. The observed limit is nevertheless worse due to the mild excess observed in the SRs. Thus, for instance, for a massless neutralino hypothesis, only the region with B(g → bbχ   TeV and three neutralino masses (1, 600 and 1000 GeV). In (a), the expected limit for a gluino mass of 1.8 TeV follows the plot axes, meaning that the whole plane is expected to be excluded at 95% CL. The same is true in (b) for a neutralino mass of 600 GeV. The dashed and solid bold lines show the 95% CL expected and observed limits, respectively. The hashing indicates which side of the line is excluded. The upper right half of the plane is forbidden by the requirement that the sum of branching ratios does not exceed 100%.

Conclusion
A search for pair-produced gluinos decaying via bottom or top squarks is presented. LHC proton-proton collision data from the full 2015 and 2016 data-taking periods are analysed, corresponding to an integrated luminosity of 36.1 fb −1 collected at √ s = 13 TeV by the ATLAS detector. The search uses multiple signal regions designed for different scenarios of gluino and LSP masses. The signal regions require several high-p T jets, of which at least three must be b-tagged, large E miss T and either zero or at least one charged lepton. Two strategies are employed: one in which the signal regions are optimised for discovery, and another one in which several non-overlapping signal regions are fitted simultaneously to achieve optimal exclusion limits for benchmark signals. For all signal regions, the background is generally dominated by tt+jets, which is normalised in dedicated control regions. No excess is found above the predicted background in any of the signal regions. Model-independent limits are set on the visible cross-section for new physics processes. Exclusion limits are set on gluino and LSP masses in two simplified models where the gluino decays exclusively asg → bbχ        [26] ATLAS Collaboration, Search for strong production of supersymmetric particles in final states with missing transverse momentum and at least three b-jets at √ s = 8 TeV proton-proton collisions with the ATLAS detector, JHEP 10 (2014) 24 [85] ATLAS Collaboration, Electron efficiency measurements with the ATLAS detector using the 2015 LHC proton-proton collision data, ATLAS-CONF-2016-024, 2016, url: https://cds.cern. ch/record/2157687.