Inclusive search for supersymmetry using razor variables in pp collisions at sqrt(s) = 13 TeV

An inclusive search for supersymmetry using razor variables is performed in events with four or more jets and no more than one lepton. The results are based on a sample of proton-proton collisions corresponding to an integrated luminosity of 2.3 inverse femtobarns collected with the CMS experiment at a center-of-mass energy of sqrt(s) = 13 TeV. No significant excess over the background prediction is observed in data, and 95% confidence level exclusion limits are placed on the masses of new heavy particles in a variety of simplified models. Assuming that pair-produced gluinos decay only via three-body processes involving third-generation quarks plus a neutralino, and that the neutralino is the lightest supersymmetric particle with a mass of 200 GeV, gluino masses below 1.6 TeV are excluded for any branching fractions for the individual gluino decay modes. For some specific decay mode scenarios, gluino masses up to 1.65 TeV are excluded. For decays to first- and second-generation quarks and a neutralino with a mass of 200 GeV, gluinos with masses up to 1.4 TeV are excluded. Pair production of top squarks decaying to a top quark and a neutralino with a mass of 100 GeV is excluded for top squark masses up to 750 GeV.


Introduction
Supersymmetry (SUSY) is a proposed extended spacetime symmetry that introduces a bosonic (fermionic) partner for every fermion (boson) in the standard model (SM) [1][2][3][4][5][6][7][8][9]. Supersymmetric extensions of the SM are particularly compelling because they yield solutions to the gauge hierarchy problem without the need for large fine tuning of fundamental parameters [10][11][12][13][14][15], exhibit gauge coupling unification [16][17][18][19][20][21], and can provide weakly interacting particle candidates for dark matter [22,23]. For SUSY to provide a "natural" solution to the gauge hierarchy problem, the three Higgsinos, two neutral and one charged, are expected to be light, and two top squarks, one bottom squark, and the gluino must have masses below a few TeV, making them potentially accessible at the CERN LHC. Previous searches for SUSY by the CMS [24][25][26][27][28][29][30] and ATLAS [31-37] Collaborations have probed SUSY particle masses near the TeV scale, and the increase in the center-of-mass energy of the LHC from 8 to 13 TeV provides an opportunity to significantly extend the sensitivity to higher SUSY particle masses [38][39][40][41][42][43][44][45][46][47][48][49][50][51]. In R-parity [52] conserving SUSY scenarios, the lightest SUSY particle (LSP) is stable and assumed to be weakly interacting. For many of these models, the experimental signatures at the LHC are characterized by an abundance of jets and a large transverse momentum imbalance, but the exact form of the final state can vary significantly, depending on the values of the unconstrained model parameters. To ensure sensitivity to a broad range of SUSY parameter space, we adopt an inclusive search strategy, categorizing events according to the number of identified leptons and b-tagged jets. The razor kinematic variables M R and R 2 [53,54] are used as search variables and are generically sensitive to pair production of massive particles with subsequent direct or cascading decays to weakly interacting stable particles. Searches for SUSY and other beyond the SM phenomena using razor variables have been performed by both the CMS [53-58] and ATLAS [59,60] Collaborations in the past. We interpret the results of the inclusive search using simplified SUSY scenarios for pair production of gluinos and top squarks. First, we consider models in which the gluino undergoes three-body decay, either to a bottom or top quarkantiquark pair and the lightest neutralino χ 0 1 , assumed to be the lightest SUSY particle; or to a bottom quark (antiquark), a top antiquark (quark), and the lightest chargino χ ± 1 , assumed to be the next-to-lightest SUSY particle (NLSP). The NLSP is assumed to have a mass that is 5 GeV larger than the mass of the LSP, motivated by the fact that in many natural SUSY scenarios the lightest chargino and the two lightest neutralinos are Higgsino-like and quasi-degenerate [61]. The NLSP decays to an LSP and an off-shell W boson, whose decay products mostly have too low momentum to be identifiable. The specific choice of the NLSP-LSP mass splitting does not have a large impact on the results of the interpretation. The full range of branching fractions to the three possible decay modes (bb χ 0 1 , bt χ + 1 or bt χ − 1 , and tt χ 0 1 ) is considered, assuming that these sum to 100%. We also consider a model in which the gluino decays to a first-or second-generation quark-antiquark pair and the LSP. Finally, we consider top squark pair production with the top squark decaying to a top quark and the LSP. Diagrams of these simplified model processes are shown in Fig. 1. This paper is organized as follows. Section 2 presents an overview of the CMS detector. A description of simulated signal and background samples is given in Section 3. Section 4 describes physics object reconstruction and the event selection. Section 5 describes the analysis strategy and razor variables, and the background estimation techniques used in this analysis are described in Section 6. Section 7 covers the systematic uncertainties. Finally, our results and their interpretation are presented in Section 8, followed by a summary in Section 9.  Figure 1: Diagrams displaying the distinct event topologies of gluino (all but last) and top squark (last) pair production considered in this paper. Diagrams corresponding to charge conjugate decay modes are implied. The symbol W * is used to denote a virtual W boson.

The CMS detector
The central feature of the CMS detector is a superconducting solenoid of 6 m internal diameter, providing a magnetic field of 3.8 T. Within the solenoid volume are a silicon pixel and a silicon strip tracker, a lead tungstate crystal electromagnetic calorimeter (ECAL), and a brass and scintillator hadron calorimeter (HCAL), each comprising a barrel and two endcap sections. Muons are measured in gas-ionization detectors embedded in the magnet steel flux-return yoke outside the solenoid. Extensive forward calorimetry complements the coverage provided by the barrel and endcap detectors. Jets are reconstructed within the pseudorapidity region |η| < 5 covered by the ECAL and HCAL, where η ≡ − ln[tan(θ/2)] and θ is the polar angle of the trajectory of the particle with respect to the counterclockwise beam direction. Electrons and muons are reconstructed in the region with |η| < 2.5 and 2.4, respectively. Events are selected by a two-level trigger system. The first level is based on a hardware trigger, followed by a software-based high level trigger. A more detailed description of the CMS detector, together with a definition of the coordinate system used and the relevant kinematic variables, can be found in Ref. [62].

Simulated event samples
Simulated Monte Carlo (MC) samples are used for modeling of the SM backgrounds in the search regions and for calculating the selection efficiencies for SUSY signal models. The production of tt+jets, W+jets, Z+jets, γ+jets, and QCD multijet events, as well as production of gluino and top squark pairs, is simulated with the MC generator MADGRAPH v5 [63]. Single top quark events are modeled at next-to-leading order (NLO) with MADGRAPH aMC@NLO v2.2 [64] for the s-channel, and with POWHEG v2 [65,66] for the t-channel and W-associated production. Contributions from ttW, ttZ are also simulated with MADGRAPH aMC@NLO v2.2. Simulated events are interfaced with PYTHIA v8.2 [67] for fragmentation and parton show-

Object reconstruction and selection
Physics objects are defined using the particle-flow (PF) algorithm [77,78]. The PF algorithm reconstructs and identifies each individual particle with an optimized combination of information from the various elements of the CMS detector. All reconstructed PF candidates are clustered into jets using the anti-k T algorithm [79,80] with a distance parameter of 0.4. The jet momentum is determined as the vector sum of all particle momenta in the jet, and jet-energy corrections are derived from simulation and confirmed by in-situ measurements of the energy balance in dijet and photon+jet events. Jets are required to pass loose identification criteria on the jet composition designed to reject spurious signals arising from noise and failures in the event reconstruction [81,82]. For this search, we consider jets with transverse momentum p T > 40 GeV and |η| < 3.0. The missing transverse momentum vector p miss T is defined as the projection on the plane perpendicular to the beams of the negative vector sum of the momenta of all reconstructed PF candidates in an event. Its magnitude is referred to as the missing transverse energy E miss T . Electrons are reconstructed by associating a cluster of energy deposited in the ECAL with a reconstructed track [83], and are required to have p T > 5 GeV and |η| < 2.5. A "tight" selection used to identify prompt electrons is based on requirements on the electromagnetic shower shape, the geometric matching of the track to the calorimeter cluster, the track quality and impact parameter, and isolation. The isolation of electrons and muons is defined as the scalar sum of the transverse momenta of all neutral and charged PF candidates within a cone ∆R = √ (∆η) 2 + (∆φ) 2 along the lepton direction. The variable is corrected for the effects of pileup using an effective area correction [84], and the cone size ∆R shrinks with increasing lepton p T according to The use of the lepton p T dependent isolation cone enhances the efficiency for identifying leptons in events containing a large amount of hadronic energy, such as those with tt production. For tight electrons, the isolation is required to be less than 10% of the electron p T . The selection efficiency for tight electrons increases from 60% for p T around 20 GeV to 70% for p T around by grouping selected leptons and jets in the event into two "megajets", whose four-momenta are defined as the vector sum of the four-momenta of their constituent physics objects [55]. The clustering algorithm selects the grouping that minimizes the sum of the squares of the invariant masses of the two megajets. We define the razor variables M R and M R T as where p j i , p j i T , and p j i z are the momenta of the ith megajet and its transverse and longitudinal components with respect to the beam axis, respectively. The dimensionless variable R is defined as For a typical SUSY decay of a superpartner q decaying into an invisible neutralino χ 0 1 and the standard model partner q, the mass variable M R peaks at a characteristic mass scale [53,54] (m 2 q − m 2 χ 0 1 )/m χ 0 1 . For standard model background processes, the distribution of M R has an exponentially falling shape. The variable R 2 is related to the missing transverse energy and is used to suppress QCD multijet background. The events of interest are triggered either by the presence of a high-p T electron or muon, or through dedicated hadronic triggers requiring the presence of at least two highly energetic jets and with loose thresholds on the razor variables M R and R 2 . The single-electron (single-muon) triggers require at least one isolated electron (muon) with p T > 23 (20) GeV. The isolation requirement is dropped for electrons (muons) with p T > 105 (50) GeV. The efficiencies for the single electron (muon) triggers are above 70% for p T around 25 (20) GeV, and reach a plateau above 97% for p T > 40 GeV. The efficiencies for the single electron trigger were measured in data and simulation and found to be in good agreement, as were the corresponding efficiencies for muons. Corrections for residual difference of trigger efficiency between data and MC simulation are applied to simulated samples. The hadronic razor trigger requires at least two jets with p T > 80 GeV or at least four jets with p T > 40 GeV. The events are also required to pass selections on the razor variables M R > 200 GeV and R 2 > 0.09 and on the product (M R + 300 GeV) × (R 2 + 0.25) > 240 GeV. The efficiency of the hadronic razor trigger for events passing the baseline M R and R 2 selections described below is 97% and is consistent with the prediction from MC simulation. For events in the Electron or Muon Multijet categories, the search region is defined by the selections M R > 400 GeV and R 2 > 0.15. The p T of the electron (muon) is required to be larger than 25 (20) GeV. To suppress backgrounds from the W( ν)+jets and tt processes, we require that the transverse mass M T formed by the lepton momentum and p miss T be larger than 120 GeV. For events in the Multijet category, the search uses a region defined by the selections M R > 500 GeV and R 2 > 0.25 and requires the presence of at least two jets with p T > 80 GeV within |η| < 3.0, for compatibility with the requirements imposed by the hadronic razor triggers. For QCD multijet background events, the E miss T arises mainly from mismeasurement of the energy of one of the leading jets. In such cases, the two razor megajets tend to lie in a back-to-back configuration. Therefore, to suppress the QCD multijet background we require that the azimuthal angle ∆φ R between the two razor megajets be less than 2.8 radians. Finally, events containing signatures consistent with beam-induced background or anomalous noise in the calorimeters are rejected using dedicated filters [90,91].

Background modeling
The main background processes in the search regions considered are W( ν)+jets (with = e, µ, τ), Z(νν)+jets, tt, and QCD multijet production. For event categories with zero b-tagged jets, the background is primarily composed of the W( ν)+jets and Z(νν)+jets processes, while for categories with two or more b-tagged jets it is dominated by the tt process. There are also very small contributions from the production of two or three electroweak bosons and from the production of tt in association with a W or Z boson. These contributions are summed and labeled "Other" in Fig. 2-5. We model the background using two independent methods based on control samples in data with entirely independent sets of systematic assumptions. The first method (A) is based on the use of dedicated control regions that isolate a specific background processes in order to control and correct the predictions of the MC simulation. The second method (B) is based on a fit to an assumed functional form for the shape of the observed data distribution in the two-dimensional M R -R 2 plane. These two background predictions are compared and cross-checked against each other in order to significantly enhance the robustness of the background estimate.

Method A: simulation-assisted background prediction from data
The simulation-assisted method defines dedicated control regions that isolate each of the main background processes. Data in these control regions are used to control and correct the accuracy of the MC prediction for each of the background processes. Corrections for the jet energy response and lepton momentum response are applied to the MC, as are corrections for the trigger efficiency and the selection efficiency of electrons, muons, and b-tagged jets. Any disagreement observed in these control regions is then interpreted as an inaccuracy of the MC in predicting the hadronic recoil spectrum and jet multiplicity. Two alternative formulations of the method are typically used in searches for new physics [25,30,31]. In the first formulation, the data control region yields are extrapolated to the search regions via translation factors derived from simulation. In the second formulation, simulation to data correction factors are derived in bins of the razor variables M R and R 2 and are then applied to the simulation prediction of the search region yields. The two formulations are identical and the choice of which formulation is used depends primarily on the convenience of the given data processing sequence. In both cases, the contributions from background processes other than the one under study are subtracted using the MC prediction. We employ the first formulation of the method for the estimate of the QCD background, while the second formulation is used for modeling all other major backgrounds. Details of the control regions used for each of the dominant background processes are described in the subsections below. Finally, the small contribution from rare background processes such as ttZ is modeled using simulation. Systematic uncertainties on the cross sections of these processes are propagated to the final result.

The tt and W( ν)+jets background
The control region to isolate the tt and W( ν)+jets processes is defined by requiring at least one tight electron or muon. To suppress QCD multijet background, the quantities E miss T and M T are both required to be larger than 30 GeV. To minimize contamination from potential SUSY processes and to explicitly separate the control region from the search regions, we require M T < 100 GeV. The tt enhanced control region is defined by requiring that there be at least one b-tagged jet, and the W( ν)+jets enhanced control region is defined by requiring no such b-tagged jets. Other than these b-tagged jet requirements, we place no explicit requirement on the number of jets in the event, in order to benefit from significantly larger control samples. We first derive corrections for the tt background, and then measure correc- tions for the W( ν)+jets process after first applying the corrections already obtained for the tt background in the W( ν)+jets control region. As discussed above, the corrections to the MC prediction are derived in two-dimensional bins of the M R -R 2 plane. We observe that the M R spectrum predicted by the simulation falls off less steeply than the control region data for both the tt and W( ν)+jets processes, as shown in Fig. 2. In Fig. 3, we show the two dimensional M R -R 2 distributions for data and simulation in the W( ν)+jets control region. The statistical uncertainties in the correction factors due to limited event yields in the control region bins are propagated and dominate the total uncertainty of the background prediction. For bins at large M R (near 1000 GeV), the statistical uncertainties range between 15% and 50%. Corrections to the MC simulation are first measured and applied as a function of M R and R 2 , inclusively in the number of selected jets. As our search region requires a higher multiplicity of jets, an additional correction factor is required to accurately model the jet multiplicity. We measure this additional correction factor to be 0.90 ± 0.03 by comparing the data and the MC prediction in the W( ν)+jets and tt control region for events with four or more jets. To control for possible simulation mismodeling that is correlated between the number of jets and the razor variables, we perform additional cross-checks of the M R and R 2 distributions in bins of the number of b-tagged jets in the tt and W( ν)+jets control regions for events with four or more jets. For bins that show statistically significant disagreement, the size of the disagreement is propagated as a systematic uncertainty. The typical range of these additional systematic uncertainties is between 10% and 30%. The tt and W( ν)+jets backgrounds in the zero-lepton Multijet event category are composed of events with at least one lepton in the final state, which is either out of acceptance or fails the veto electron, muon, or τ h selection. Two additional control regions are defined in order to control the accuracy of the modeling of the acceptance and efficiency for selecting electrons or muons, and τ h . We require events in the veto lepton (τ h candidate) control region to have at least one veto electron or muon (τ h candidate) selected. The M T is required to be between 30 and 100 GeV in order to suppress QCD multijet background and contamination from potential new physics processes. At least two jets with p T > 80 GeV and at least four jets with p T > 40 GeV are required, consistent with the search region requirements. Finally, we consider events with M R > 400 GeV and R 2 > 0.25. The distribution of the veto lepton p T for events in the veto lepton and veto τ h control regions are shown in Fig. 4, and demonstrate that the MC models describe well the observed data. The observed discrepancies in any bin are propagated as systematic uncertainties in the prediction of the tt and W( ν)+jets in the Multi-     jet category search region. The tt background in the Electron and Muon Multijet categories is primarily from the dilepton decay mode as the M T requirement highly suppresses the semileptonic decay mode. Corrections to the MC simulation derived from the tt control region primarily arise from semi-leptonic decays. We define an additional control region enhanced in dilepton tt decays to confirm that the MC corrections derived from a region dominated by semi-leptonic decays also apply to dilepton decays. We select events with two tight leptons, both with p T > 30 GeV, E miss T > 40 GeV, and dilepton mass larger than 20 GeV. For events with two leptons of the same flavor, we additionally veto events with a dilepton mass between 76 and 106 GeV in order to suppress background from Z boson decays. At least one b-tagged jet is required to enhance the purity for the tt process. Finally, we mimic the phase space region similar to our search region in the Electron and Muon Multijet categories by treating one lepton as having failed the identification criteria and applying the M T requirement using the other lepton. The correction factors measured in the tt control region are applied to the MC prediction of the dilepton tt cross-check region in bins of M R and R 2 . In Fig. 3 we show the M R -R 2 distribution for the dilepton tt cross-check region in events with four or more jets, and we observe no significant mismodeling by the simulation, indicating that the measured corrections are accurate.

The Z → νν background
Three independent control regions are used to predict the Z(νν)+jets background, relying on the assumption that the hadronic recoil spectrum and the jet multiplicity distribution of the Z(νν)+jets process are similar to those of the W( ν)+jets and γ+jets processes. The primary and most populated control region is the γ+jets control region, defined by selecting events with at least one photon passing loose identification and isolation requirements. The events are triggered using single-photon triggers, and the photon is required to have p T > 50 GeV. The momentum of the photon candidate in the transverse plane is added vectorially to p miss T in order to simulate an invisible particle, as one would have in the case of a Z → νν decay, and the M R and R 2 variables are computed according to this invisible decay scenario. A template fit to the distribution of σ ηη is performed to determine the contribution from misidentified photons to the γ+jets control region and is found to be about 5%, independent of the M R and R 2 . Events from the γ+jets process where the photon is produced within the cone of a jet (labeled as γ+jets fragmentation) are considered to be background and subtracted using the MC predic-   tion. Backgrounds from rarer processes such as Wγ, Zγ, and ttγ are also subtracted similarly. In Fig. 5, we show the M R distribution as well as the two-dimensional M R -R 2 distribution for the γ+jets control region, where we again observe a steeper M R falloff in the data compared to the simulation. Correction factors are derived in bins of M R and R 2 and applied to the MC prediction for the Z → νν background in the search region. The statistical uncertainties for the correction factors range between 10% and 30% and are among the dominant uncertainties for the Z → νν background prediction. Analogous to the procedure for the tt and W( ν)+jets control region, we derive an additional correction factor of 0.87 ± 0.05 to accurately describe the yield in events with four or more jets. Additional cross-checks are performed in bins of the number of b-tagged jets and systematic uncertainties ranging from 4% for events with zero b-tagged jets to 58% for events with three or more b-tagged jets are derived. The second control region, enhanced in the W( ν)+jets, is defined identically to the W( ν)+jets control region described in Section 6.1.1, except that the lepton is treated as invisible by adding its momentum vectorially to p miss T , and the M R and R 2 variables are computed accordingly. Correction factors are computed using events from this control region, and the difference between these correction factors and those computed from the γ+jets control region is propagated as a systematic uncertainty. These uncertainties range between 10% and 40% depending on the M R -R 2 bin. The third control region, enhanced in Z → + − decays, is defined by selecting events with two tight electrons or two tight muons, and requiring that the dilepton mass is between 76 and 106 GeV. Events are required to have no b-tagged jets in order to suppress tt background. The two leptons are treated as invisible by adding their momenta vectorially to p miss T . We apply the correction factors obtained from the γ+jet control region to the Z → + − MC prediction and perform a cross-check against data in this control region. No significant discrepancy between the data and the prediction is observed.

The QCD Multijet background
The QCD multijet processes contribute about 10% of the total background in the zero-lepton Multijet event category for bins with zero or one b-tagged jets. Such events enter the search regions in the tails of the E miss T distribution when the energy of one of the jets in the event is significantly under-or over-measured. In most such situations, the p miss T points either toward or away from the leading jets and therefore the two megajets tend to be in a back-to-back configuration. The search region is defined by requiring that the azimuthal angle between the two megajets ∆φ R be less than 2.8, which was found to be an optimal selection based on studies of QCD multijet and signal simulated samples. We define the control region for the QCD background process to be events with ∆φ R > 2.8, keeping all other selection requirements identical to those for the search region. The purity of the QCD multijet process in the control region is more than 70%. After subtracting the non-QCD background, we project the observed data yield in the control region to the search region using the translation factor ζ: where the numerator and denominator are the number of events passing and failing the selection on |∆φ R | < 2.8, respectively. We find that the translation factor calculated from the MC simulation decreases as a function of M R and is, to a large degree, constant as a function of R 2 . Using data events in the low R 2 region (0.15 to 0.25), dominated by QCD multijet background, we measure the translation factor ζ as a function of M R to cross-check the values obtained from the simulation. The M R dependence of ζ is modeled as the sum of a power law and a constant. This functional shape is fitted to the values of ζ calculated from the MC. A systematic uncertainty of 87% is propagated, covering both the spread around the fitted model as a function of M R and R 2 in simulation, and the difference between the values measured in simulation and data. The function used for ζ and the values measured in data and simulation are shown in Fig. 6. We perform two additional cross-checks on the accuracy of the MC prediction for ζ in control regions dominated by processes similar to the QCD multijet background with no invisible neutrinos in the final state. The first cross-check is performed on a dimuon control region enhanced in Z → µ + µ − decays, and the second cross-check is performed on a dijet control region enhanced in QCD dijet events. In both cases, the events at large R 2 result from cases similar to our search region where the energy of a leading jet is severely mismeasured. We compare the values of ζ measured in these data control regions to the values predicted by the simulation and observe an agreement at the 20% level, well within the systematic uncertainty of 87% assigned to the QCD background estimate.

Method B: fit-based background prediction
The second background prediction method is based on a fit to the data with an assumed functional form for the shape of the background distribution in the M R -R 2 plane. Based on past studies [54, 56], the shape of the background in the M R and R 2 variables is found to be well described by the following functional form: where M 0 R , R 2 0 , b, and n are free parameters. In the original study [54], this function with n fixed to 1 was used to model the data in each category. The function choice was motivated by the observation that for n = 1, the function projects to an exponential both on R 2 and M R , and b is proportional to the exponential rate parameter in each one-dimensional projection. The generalized function in Eq. (6) was found to be in better agreement with the SM backgrounds over a larger range of R 2 and M R [56]. The two parameters b and n determine the tail of the distribution in the two-dimensional plane, while the M 0 R (R 2 0 ) parameter affects the tail of the one-dimensional projection on R 2 (M R ). Background estimation is performed using an extended, binned, maximum likelihood fit to the M R and R 2 distribution in one of two ways: • A fit to the data in the sideband regions in M R and R 2 , defined more precisely below, as a model-independent way to look for excesses or discrepancies. The fit is performed using only the data in the sideband, and the functional form is extrapolated to the full M R and R 2 plane.
• A fit to the data in the full search region in M R and R 2 under background-only and signal-plus-background hypotheses, following a modified frequentist approach (LHC CL s ) [92][93][94][95][96] to interpret the data in the context of particular SUSY simplified models.
The sideband region is defined to be 100 GeV in width in M R and 0.05 in R 2 . Explicitly, for the Multijet event category, it comprises the region 500 GeV < M R < 600 GeV and R 2 > 0.3, plus the region M R > 500 GeV and 0.25 < R 2 < 0.3. For the Muon and Electron Multijet event categories, it comprises the region 400 GeV < M R < 500 GeV and R 2 > 0.2, plus the region M R > 400 GeV and 0.15 < R 2 < 0.2. For each event category, we fit the two-dimensional distribution of M R and R 2 in the sideband region using the above functional form, separately for events with zero, one, two, and three or more b-tagged jets. The normalization in each event category and each b-tagged jet bin is independently varied in the fit. Due to the lack of data events in the category with three or more b-tagged jets, we constrain the shape in this category to be related to the shape for events with two b-tagged jets as follows: where f 2b SM (M R , R 2 ) and f ≥3b SM (M R , R 2 ) are the probability density functions for events with two and with three or more b-tagged jets, respectively; M offset R is the lowest M R value in a particular event category; and m M R is a floating parameter constrained by a Gaussian distribution centered at the value measured using the simulation and with a 100% uncertainty. The above form for the shape of the background events with three or more b-tagged jets is verified in simulation. Numerous tests are performed to establish the robustness of the fit model in adequately describing the underlying distributions. To demonstrate that the background model gives an accurate description of the background distributions, we construct a representative data set using MC samples, and perform the background fit using the form given by Eq. (6). Goodness of fit is evaluated by comparing the background prediction from the fit with the prediction from the simulation. This procedure is performed separately for each of the search categories and we find that the fit function yields an accurate representation of the background predicted by the simulation. We also observe that the fit model is insensitive to variations of the background composition predicted by the simulation in each event category by altering relative contributions of the dominant backgrounds, performing a new fit with the alternative background composition, and comparing the new fit results to the nominal fit result. The contributions of the main tt, W( ν)+jets, and Z(νν) backgrounds are varied by 30%, and the rare backgrounds from QCD multijet and ttZ processes are varied by 100%. For the Muon and Electron Multijet event categories, we also vary the contributions from the dileptonic and semi-leptonic decays of the tt background separately by 30%. In each of these tests, we observe that the chosen functional form can adequately describe the shapes of the M R and R 2 distributions as predicted by the modified MC simulation. Additional pseudo-experiment studies are performed comparing the background prediction from the sideband fit and the full region fit to evaluate the average deviation between the two fit predictions. We observe that the sideband fit and the full region fit predictions in the signal-sensitive region differ by up to 15% and we propagate an additional systematic uncertainty to the sideband fit background prediction to cover this average difference. To illustrate method B, we present the data and fit-based background predictions in Fig. 7, for events in the 2 b-tag and ≥3 b-tag Multijet categories. The number of events observed in data is compared to the prediction from the sideband fit in the M R and R 2 bins. To quantify the agreement between the background model and the observation, we generate alternative sets of background shape parameters from the covariance matrix calculated by the fit. An ensemble of pseudo-experiment data sets is created, generating random (M R , R 2 ) pairs distributed according to each of these alternative shapes. For each M R -R 2 bin, the distribution of the predicted yields from the ensemble of pseudo-experiments is compared to the observed yield in data. The agreement between the predicted and the observed yields is described as a two-sided p-value and translated into the corresponding number of standard deviations for a normal distribution. Positive (negative) significance indicates the observed yield is larger (smaller) than the predicted one. We find that the pattern of differences between data and background predictions in the different bins considered is consistent with statistical fluctuations. To demonstrate that the model-independent sideband fit procedure used in the analysis would be sensitive to the presence of a signal, we perform a signal injection test. We sample a signal-plus-background pseudo-data set and perform a background-only fit in the sideband. We show one illustrative example of such a test in Fig. 8        corresponding to gluino pair production, in which each gluino decays to a neutralino and a bb pair with m g = 1.4 TeV and m χ 0 1 = 100 GeV. The deviations with respect to the fit predictions are shown for the 2 b-tag and ≥3 b-tag Multijet categories. We observe characteristic patterns of excesses in two adjacent groups of bins neighboring in M R .

Comparison of two methods
The background predictions obtained from methods A and B are systematically compared in all of the search region categories. For method B, the model-independent fit to the sideband is used for this comparison. In Fig. 9, we show the comparison of the two background predictions for two example event categories. The predictions from the two methods agree within the uncertainties of each method. The uncertainty from the fit-based method tends to be slightly larger at high M R and R 2 due to the additional uncertainty in the exact shape of the tail of the distribution, as the n and b parameters are not strongly constrained by the sideband data. The two background predictions use methods based on data that make very different systematic assumptions. Method A assumes that corrections to the simulation prediction measured in control regions apply also to the signal regions, while method B assumes that the shape of the background distribution in M R and R 2 is well described by a particular exponentially falling functional form. The agreement observed between predictions obtained using these two very different methods significantly enhances the confidence of the background modeling, and also validates the respective assumptions.

Systematic uncertainties
Various systematic uncertainties are considered in the evaluation of the signal and background predictions. Different types of systematic uncertainties are considered for the two different background models. For method A, the largest uncertainties arise from the precision with which the MC corrections are measured. The dominant uncertainties in the correction factors result from statistical uncertainties due to the limited size of the control region event sample. We also propagate systematic uncertainties in the theoretical cross-section for the small residual backgrounds present in the control regions, and they contribute 2-5% to the correction factor uncertainty. Additional systematic uncertainties are computed from the procedure that tests that the accuracy of the MC corrections as a function of (M R , R 2 ), and the number of b-tagged jets in events with four or more jets. The total uncertainty from this procedure ranges from 10% for the most populated bins to 50% and 100% for the least populated bins. For the Z → νν process, we also propagate the difference in the correction factors measured in the three alternative control regions as a systematic uncertainty, intended to estimate the possible differences in the simulation mismodeling of the hadronic recoil for the γ+jets process and the Z(νν)+jets process. These systematic uncertainties range from 10 to 40%. For the QCD background prediction the statistical uncertainty due to limited event counts in the ∆φ R > 2.8 control regions and the systematic uncertainty of 87% in the translation factor ζ are propagated. For method B, the systematic uncertainties in the background are propagated as part of the maximum likelihood fit procedure. For each event category, the background shape in M R and R 2 is described by four independent parameters: two that control the exponential fall off and two that control the behavior of the nonexponential tail. Systematic uncertainties in the background are propagated through the profiling of these unconstrained shape parameters. For more populated bins, such as the 0 b-tag and 1 b-tag bins in the Multijet category, the systematic uncertainties range from about 30% at low M R and R 2 to about 70% at high M R and R 2 . For sparsely populated bins such as the 3-or-more b-tag bin in the Muon Multijet or Electron Multijet categories, the systematic uncertainties range from about 60% at low M R and R 2 to more than 200% at high M R and R 2 .    Systematic uncertainties due to instrumental and theoretical effects are propagated as shape uncertainties in the signal predictions for methods A and B, and on the background predictions for method A. The background prediction from method B is not affected by these uncertainties as the shape and normalization are measured from data. Uncertainties in the trigger and lepton selection efficiency, and the integrated luminosity [97] primarily affect the total normalization. Uncertainties in the b-tagging efficiency affect the relative yields between different b-tag categories. The uncertainties from missing higher-order corrections and the uncertainties in the jet energy and lepton momentum scale affect the shapes of the M R and R 2 distributions. For the signal predictions, we also propagate systematic uncertainties due to possible inaccuracies of the fast simulation in modeling the lepton selection and b tagging efficiencies. These uncertainties were evaluated by comparing the tt and signal GEANT based MC samples with those that used fast simulation. Finally, we propagate an uncertainty in the modeling of initial-state radiation for signal predictions, that ranges from 15% for signal events with recoil between 400 and 600 GeV to 30% for events with recoil above 600 GeV. The systematic uncertainties and their typical impact on the background and signal predictions are summarized in Table 1.

Results and interpretations
We present results of the search using method A as it provides slightly better sensitivity. The two-dimensional M R -R 2 distributions for the search regions in the Multijet, Electron Multijet, and Muon Multijet categories observed in data are shown in Figures 10-15, along with the background prediction from method A. We observe no statistically significant discrepancies and interpret the null search result using method A by determining the 95% confidence level (CL) upper limits on the production cross sections of the SUSY models presented in Section 1 using a global likelihood determined by combining the likelihoods of the different search boxes and sidebands. Following the LHC CL s procedure [96], we use the profile likelihood ratio test statistic and the asymptotic formula to evaluate the 95% CL observed and expected limits on the SUSY production cross section σ. Systematic uncertainties are taken into account by incor-porating nuisance parameters θ, representing different sources of systematic uncertainty, into the likelihood function L(σ, θ). For each signal model the simulated SUSY events are used to estimate the effect of possible signal contamination in the analysis control regions, and the method A background prediction is corrected accordingly. To determine a confidence interval for σ, we construct the profile likelihood ratio test statistic −2 ln[L(σ,θ σ )/L(σ,θ)] as a function of σ, whereθ σ refers to the conditional maximum likelihood estimators of θ assuming a given value σ, andσ andθ correspond to the global maximum of the likelihood. Then for example, a 68% confidence interval for σ can be taken as the region for which the test statistic is less than 1. By allowing each nuisance parameter to vary, the test statistic curve is wider, reflecting the systematic uncertainty arising from each source, and resulting in a larger confidence interval for σ. First, we consider the scenario of gluino pair production decaying to third-generation quarks. Gluino decays to the third-generation are enhanced if the masses of the third-generation squarks are significantly lighter than those of the first two generations, a scenario that is strongly motivated in natural SUSY models [61,[98][99][100]. Prompted by this, we consider the three decay modes: • g → bb χ 0 ; • g → tt χ 0 ; • g → bt χ + 1 → btW * + χ 0 1 or charge conjugate, where W * denotes a virtual W boson. Due to a technical limitation inherent in the event generator, we consider these three decay modes for |m g − m χ 0 1 | ≥ 225 GeV. For |m g − m χ 0 1 | < 225 GeV, we only consider the g → bb χ 0 decay mode. The three-body gluino decays considered here capture all of the possible final states within this natural SUSY context including those of twobody gluino decays with intermediate top or bottom squarks. Past studies have shown that LHC searches exhibit a similar sensitivity to three-body and two-body gluino decays with a only a weak dependence on the intermediate squark mass [40]. We perform a scan over all possible branching fractions to these three decay modes and compute limits on the production cross section under each such scenario. The production cross section limits for a few characteristic branching fraction scan points are shown on the left of Fig. 16 as a function of the gluino and neutralino masses. We find a range of excluded regions for different branching fraction assumptions and generally observe the strongest limits for the g → bb χ 0 1 decay mode over the full two-dimensional mass plane and the weakest limits for the g → tt χ 0 1 decay mode. For scenarios that include the intermediate decay χ ± 1 → W * ± χ 0 1 and small values of m χ 0 1 the sensitivity is reduced because the LSP carries very little momentum in both the NLSP rest frame and the laboratory frame, resulting in small values of E miss T and R 2 . By considering the limits obtained for all scanned branching fractions, we calculate the exclusion limits valid for any assumption on the branching fractions, presented on the right of Fig. 16. For an LSP with mass of a few hundred GeV, we exclude pair production of gluinos decaying to third-generation quarks for mass below about 1.6 TeV. This result represents a unique attempt to obtain a branching fraction independent limit on gluino pair production at the LHC for the scenario in which gluino decays are dominated by three-body decays to third-generation quarks and a neutralino LSP. In Figure 17, we present additional interpretations for simplified model scenarios of interest. On the left, we show the production cross section limits on gluino pair production where the gluino decays to two light-flavored quarks and the LSP, and on the right we show the production cross section limits on top squark pair production where the top squark decays to a top quark and the LSP. For a very light LSP, we exclude top squark production with mass below 750 GeV.      Figure 11: The M R -R 2 distribution observed in data is shown along with the background prediction obtained from method A for the Multijet event category in the 2 b-tag (upper) and ≥3 b-tag (lower) bins. A detailed explanation of the panels is given in the caption of Fig. 10.     Figure 13: The M R -R 2 distribution observed in data is shown along with the background prediction obtained from method A for the Muon Multijet event category in the 2 b-tag (upper) and ≥3 b-tag (lower) bins. A detailed explanation of the panels is given in the caption of Fig. 10.   Figure 14: The M R -R 2 distribution observed in data is shown along with the background prediction obtained from method A for the Electron Multijet event category in the 0 b-tag (upper) and 1 b-tag (lower) bins. A detailed explanation of the panels is given in the caption of Fig. 10.   Figure 15: The M R -R 2 distribution observed in data is shown along with the background prediction obtained from method A for the Electron Multijet event category in the 2 b-tag (upper) and ≥3 b-tag (lower) bins. A detailed explanation of the panels is given in the caption of Fig. 10.
Observed Expected  Figure 16: (Left) the expected and observed 95% confidence level (CL) upper limits on the production cross section for gluino pair production decaying to third-generation quarks under various assumptions of the branching fractions. The two gray dashed diagonal lines correspond to |m g − m χ 0 1 | = 25 GeV, which is where the scan ends for the g → bb χ 0 1 decay mode, and |m g − m χ 0 1 | = 225 GeV, which is where the scan ends for the remaining modes due to a technical limitation inherent in the event generator. For |m g − m χ 0 1 | < 225 GeV, we only consider the g → bb χ 0 1 decay mode. (Right) the analogous upper limits on the gluino pair production cross section valid for any values of the gluino decay branching fractions.  , and as a result the precise determination of the cross section upper limit is uncertain because of the finite granularity of the available MC samples in this region of the (m t , m χ 0 1 ) plane.

Summary
We have presented an inclusive search for supersymmetry in events with no more than one lepton, a large multiplicity of energetic jets, and missing transverse energy. The search is sensitive to a broad range of SUSY scenarios including pair production of gluinos and top squarks. The event categorization in the number of leptons and the number of b-tagged jets enhances the search sensitivity for a variety of different SUSY signal scenarios. Two alternative background estimation methods are presented, both based on transfer factors between data control regions and the search regions, but having very different systematic assumptions: one relying on the simulation and associated corrections derived in the control regions, and the other relying on the accuracy of an assumed functional form for the shape of background distribution in the M R and R 2 variables. The two predictions agree within their uncertainties, thereby demonstrating the robustness of the background modeling. No significant deviations from the predicted standard model background are observed in any of the search regions, and this result is interpreted in the context of simplified models of gluino or top squark pair production. For decays to a top quark and an LSP with a mass of 100 GeV, we exclude top squarks with masses below 750 GeV. Considering separately the decays to bottom quarks and the LSP or first-and second-generation quarks and the LSP, gluino masses up to 1.65 TeV or 1.4 TeV are excluded, respectively. Furthermore, this search goes beyond the existing simplified model paradigm by interpreting results in a broader context inspired by natural SUSY, with multiple gluino decay modes considered simultaneously. By scanning over all possible branching fractions for three-body gluino decays to third generation quarks, exclusion limits are derived on gluino pair production that are valid for any values of the gluino decay branching fractions. For a chargino NLSP nearly degenerate in mass with the LSP and LSP masses in the range between 200 and 600 GeV, we exclude gluinos with mass below 1.55 to 1.6 TeV, regardless of their decays. This result is a more generic constraint on gluino production than previously reported at the LHC.

Acknowledgments
We congratulate our colleagues in the CERN accelerator departments for the excellent performance of the LHC and thank the technical and administrative staffs at CERN and at other CMS institutes for their contributions to the success of the CMS effort. In addition, we gratefully acknowledge the computing centers and personnel of the Worldwide LHC Computing Grid for delivering so effectively the computing infrastructure essential to our analyses. Finally, we acknowledge the enduring support for the construction and operation of the LHC and the CMS detector provided by the following funding agencies: the Austrian Federal Ministry of Science, Research and Economy and the Austrian Science Fund; the Belgian    [50] ATLAS Collaboration, "Search for bottom squark pair production in proton-proton collisions at √ s = 13 TeV with the ATLAS detector", (2016). arXiv:1606.08772. Submitted to EPJC.
[51] ATLAS Collaboration, "Search for supersymmetry in a final state containing two photons and missing transverse momentum in √ s = 13 TeV pp collisions at the LHC using the ATLAS detector", (2016). arXiv:1606.09150. Submitted to EPJC.
[52] G. R. Farrar and P. Fayet, "Phenomenology of the production, decay, and detection of new hadronic states associated with supersymmetry", Phys.

A Results of method B fit-based background prediction
In Section 6.2, we detail the fit-based background prediction methodology and present the model-independent SUSY search results in the 2 b-tag and ≥3 b-tag bins of the Multijet category in Fig. 7. In Figs. [18][19][20][21][22] in this Appendix, we present the results of the search for SUSY signal events in the remaining categories, namely the 0 b-tag and 1 b-tag bins of the Multijet, the Muon Multijet, and Electron Multijet categories. No statistically significant deviations from the expected background predictions are observed in these categories in data.