Search for charged Higgs bosons decaying into top and bottom quarks at s=13\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \sqrt{s}=13 $$\end{document} TeV with the ATLAS detector

A search for charged Higgs bosons heavier than the top quark and decaying via H± → tb is presented. The data analysed corresponds to 36.1 fb−1 of pp collisions at s=13\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \sqrt{s}=13 $$\end{document} TeV and was recorded with the ATLAS detector at the LHC in 2015 and 2016. The production of a charged Higgs boson in association with a top quark and a bottom quark, pp → tbH±, is explored in the mass range from mH± = 200 to 2000 GeV using multi-jet final states with one or two electrons or muons. Events are categorised according to the multiplicity of jets and how likely these are to have originated from hadronisation of a bottom quark. Multivariate techniques are used to discriminate between signal and background events. No significant excess above the background-only hypothesis is observed and exclusion limits are derived for the production cross-section times branching ratio of a charged Higgs boson as a function of its mass, which range from 2.9 pb at mH± = 200 GeV to 0.070 pb at mH± = 2000 GeV. The results are interpreted in two benchmark scenarios of the Minimal Supersymmetric Standard Model.


Introduction
Following the discovery of a Higgs boson, H, with a mass of around 125 GeV and consistent with the Standard Model (SM) [1][2][3] at the Large Hadron Collider (LHC) in 2012 [4] a key question is whether this Higgs boson is the only Higgs boson, or the first observed physical state of an extended Higgs sector. No charged fundamental scalar boson exists in the SM, but many beyond the Standard Model (BSM) scenarios contain an extended Higgs sector with at least one set of charged Higgs bosons, H + and H − , in particular two-Higgs-doublet models (2HDM) [5][6][7][8] and models containing Higgs triplets [9][10][11][12][13].
The production mechanisms and decay modes of a charged Higgs boson1 depend on its mass, m H + . This analysis searches for heavy charged Higgs bosons with m H + > m t + m b , where m t and m b are the masses of the top and bottom quarks, respectively. The dominant production mode is expected to be in association with a top quark and a bottom quark (tbH + ), as illustrated in Figure 1. In the 2HDM, H + production and decay at tree level depend on its mass and two parameters: tanβ and α, which are the ratio of the vacuum expectation values of the two Higgs doublets and the mixing angle between the CP-even Higgs bosons, respectively. The dominant decay mode for heavy charged Higgs bosons is H + → tb in a broad range of models [14,15]. In particular, this is the preferred decay mode in both the decoupling limit scenario and the alignment limit cos(β − α) ≈ 0, where the lightest CP-even neutral Higgs boson of the extended Higgs sector has properties similar to those of the SM Higgs boson [7]. For lower m H + , the dominant decay mode is H + → τν. It is also predicted that this decay mode becomes more relevant as the value of tanβ increases, irrespective of m H + . Therefore, the H + → tb and H + → τν decays naturally complement each other in searches for charged Higgs bosons.
Limits on charged Higgs boson production have been obtained by many experiments, such as the LEP experiments with upper limits on H + production in the mass range 40-100 GeV [16], and CDF and DØ at the Tevatron that set upper limits on the branching ratio B(t → bH + ) for 80 GeV < m H + < 150 GeV [17,18]. The CMS Collaboration has performed direct searches for heavy charged Higgs bosons in 8 TeV protonproton (pp) collisions. By assuming the branching ratio B(H + → tb) = 1, an upper limit of 2.0-0.13 pb was obtained for the production cross-section σ(pp → tbH + ) for 180 GeV < m H + < 600 GeV [19]. The ATLAS Collaboration has searched for similar heavy charged Higgs boson production in the H + → tb decay channel at 8 TeV, setting upper limits on the production cross-section times the H + → tb branching ratio of 6-0.2 pb for 200 GeV < m H + < 600 GeV [20]. Indirect constraints can be obtained from the measurement of flavour-physics observables sensitive to charged Higgs boson exchange. Such observables include the relative branching ratios of B or K meson decays, B meson mixing parameters, the ratio of the Z decay partial widths Γ(Z → bb)/Γ(Z → hadrons), as well as the measurements of b → sγ decays [21,22]. The relative branching ratio R(D ( * ) ) = B(B → D ( * ) τν)/B(B → D ( * ) ν), where denotes e or µ, are especially sensitive to contributions from new physics. Measurements from BaBar [23] exclude H + for all m H + and tanβ values in a Type-II 2HDM. However, more recent measurements from Belle [24][25][26] and LHCb [27] place a weaker constraint on the allowed range of m H + /tanβ values. A global fit combining the most recent flavour-physics results [22] sets a lower limit at 95% confidence level on the charged Higgs boson mass of m H + 600 GeV for tanβ > 1 and m H + 650 GeV for lower tanβ values, assuming a Type-II 2HDM. This paper presents a search for H + production in the H + → tb decay mode using pp collisions at √ s = 13 TeV. Events with one charged lepton ( = e, µ) and jets in the final state ( +jets final state) and events with two charged leptons and jets in the final state ( final state) are considered. Exclusive regions are defined according to the number of jets and those that are tagged as originating from the hadronisation of a b-quark. In order to separate the signal from the SM background, multivariate discriminants are employed in the regions where the signal contributions are expected to be largest. Limits on the H + → tb production cross-section are set by means of a simultaneous fit of binned distributions of multivariate discriminants in the signal-rich regions and inclusive event yields in the signal-depleted regions. The results are interpreted in two benchmark scenarios of the Minimal Supersymmetric Standard Model (MSSM): the m mod− h scenario [28] and the hMSSM [29]. Both scenarios exploit the MSSM in such a way that the light CP-even Higgs boson can be interpreted as the observed Higgs boson with m H = 125 GeV. Limits on the value of tanβ are extracted as a function of the charged Higgs boson mass. Finally, the excluded range of m H + and tanβ values from the H + → tb and H + → τν [30] searches at √ s = 13 TeVare superimposed, providing a summary of the ATLAS sensitivity to H + through the two decay modes.
The paper is organised as follows. Section 2 briefly describes the ATLAS detector. The samples of simulated events used for the analysis are summarised in Section 3. Section 4 presents the reconstruction of objects in ATLAS and the event selection. Section 5 describes the analysis strategy while systematic uncertainties are discussed in Section 6. The statistical analysis of the data is described in Section 7 and the results are presented in Section 8. Finally, a summary is given in Section 9.

ATLAS detector
The ATLAS detector [31] at the LHC is a multipurpose particle detector with a forward-backward symmetric cylindrical geometry and near 4π coverage around the collision point. 2 The ATLAS detector consists of an inner tracking detector (ID) surrounded by a thin superconducting solenoid producing a 2 T axial magnetic field, electromagnetic (EM) and hadronic calorimeters, and an external muon spectrometer (MS) incorporating three large toroid magnet assemblies. The ID contains a high-granularity silicon pixel detector, including an insertable B-layer [32] added in 2014 as a new innermost layer, and a silicon microstrip tracker, providing precision tracking in the pseudorapidity range |η| < 2.5. The silicon detectors are complemented by a transition radiation tracker providing tracking and electron identification information for |η| < 2.0. The EM sampling calorimeter uses lead as the absorber material and liquid argon (LAr) as the active medium, and is divided into barrel (|η| < 1.47) and endcap (1.37 < |η| < 3.20) regions. Hadron calorimetry is also based on the sampling technique, with scintillator tiles or LAr as the active medium, and with steel, copper, or tungsten as the absorber material. The calorimeters cover |η| < 4.9. The MS measures the deflection of muons with |η| < 2.7 using multiple layers of highprecision tracking chambers located in a toroidal field in the central and endcap regions of ATLAS. The field integral of the toroids ranges between 2.0 and 6.0 T m across most of the detector. The MS is also instrumented with separate trigger chambers covering |η| < 2.4. A two-level trigger system, with the first level implemented in custom hardware and followed by a software-based second level, is used to reduce the trigger rate to around 1 kHz for offline storage [33].

Signal and background modelling
The tbH + process was modelled with M G 5_aMC@NLO (MG5_ MC) [34] at next-to-leading order (NLO) in QCD [35] using a four-flavour scheme (4FS) implementation with the NNPDF2.3NLO [36] parton distribution function (PDF).3 Parton showering and hadronisation were modelled by P 8.186 [37] 2 ATLAS uses a right-handed coordinate system with its origin at the nominal interaction point (IP) in the centre of the detector and the z-axis along the beam pipe. The x-axis points from the IP to the centre of the LHC ring, and the y-axis points upwards. 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 θ as η = − ln tan(θ/2). Angular distance is measured in units of ∆R ≡ (∆η) 2 + (∆φ) 2 (pseudorapidity and azimuthal angle). Alternatively, the distance ∆R y ≡ (∆y) 2 + (∆φ) 2 is used, where y = 0.5 ln [(E + p z ) /(E − p z )] is the rapidity of a particle of energy E and momentum component p z along the beam axis. 3 Five-flavour scheme (5FS) PDFs consider b-quarks as a source of incoming partons and the b-quarks are therefore assumed to be massless. In contrast, 4FS PDFs only include lighter quarks and gluons, allowing the b-quark mass to be taken into account properly in the matrix element calculation.
with the A14 [38] set of underlying-event (UE) related parameters tuned to ATLAS data (tune). For the simulation of the tbH + process, the narrow-width approximation was used. This assumption has a negligible impact on the analysis for the models considered in this paper, as the experimental resolution is much larger than the H + natural width. Interference with the SM tt + bb background is neglected.
Altogether 18 H + mass hypotheses are used, with 25 GeV mass steps between an H + mass of 200 GeV and 300 GeV, 50 GeV steps between 300 GeV and 400 GeV, 100 GeV steps between 400 GeV and 1000 GeV and 200 GeV steps from 1000 GeV to 2000 GeV. The step sizes are selected to match the expected resolution of the H + signal. The samples were processed with a fast simulation of the ATLAS detector [39]. Unless otherwise indicated, the cross-section of the signal is set to 1 pb, for easy rescaling to various model predictions. Only the H + decay into tb is considered, and the top quark decays according to the SM predictions.
The nominal sample used to model the tt background was generated using the P -B v2 NLO-in-QCD generator [40][41][42][43], referred to as P in the remainder of this article, with the NNPDF3.0NLO PDF set [44]. The h damp parameter, which controls the transverse momentum p T of the first additional emission beyond the Born configuration, was set to 1.5 times the top quark mass [45]. Parton shower and hadronisation were modelled by P 8.210 [46] with the A14 UE tune. The sample was normalised to the ++2.0 [47] theoretical cross-section of 832 +46 −51 pb, calculated at next-to-next-to-leading order (NNLO) in QCD including resummation of next-to-next-to-leading logarithmic (NNLL) soft gluon terms [48][49][50][51][52]. The generation of the tt sample was performed inclusively, with all possible flavours of additional jets produced. The decay of cand b-hadrons was simulated with the E G v1.2.0 [53] program. The tt + jets background is categorised according to the flavour of additional jets in the event, using the same procedure as described in Ref. [54]. The tt + additional heavy-flavour (HF) jets background is subdivided into the categories tt + ≥ 1b and tt + ≥ 1c, depending on whether the additional HF jets originate from hadrons containing bor c-quarks. Particle jets were reconstructed from stable particles (mean lifetime τ > 3 × 10 −11 seconds) at generator level using the anti-k t algorithm [55] with a radius parameter of 0.4, and were required to have p T > 15 GeV and |η| < 2.5. If at least one particle-level jet in the event is matched (∆R < 0.3) to a b-hadron (not originating from a t-decay) with p T > 5 GeV, the event is categorised as tt + ≥1b. In the remaining events, if at least one jet is matched to a c-hadron (not originating from a W decay) but no b-hadron, the event is categorised as tt + ≥1c. Events with tt + jets that belong to neither the tt + ≥1b nor tt + ≥1c category are called tt + light events.
For the tt + ≥1b process, subcategories are defined in accord with the matching between particle-level jets and the b-hadrons not from t-decay: events where exactly two jets are matched to b-hadrons (tt + bb), events where exactly one jet is matched to a b-hadron (tt + b), events where exactly one jet is matched to two or more b-hadrons (tt + B), and all other events (tt + ≥3b). Events where the additional HF jets can only be matched to b-hadrons from multi-parton interactions and final-state gluon radiation are considered separately and labelled as tt + b (MPI/FSR).
To model the irreducible tt + ≥ 1b background to the highest available precision, the tt + ≥ 1b events from the nominal P +P 8 simulation are reweighted to an NLO prediction of ttbb including parton showering and hadronisation from S 2.1.1 [56,57] with O L [58]. This sample was generated using the 4FS PDF set CT10F4 [59]. The renormalisation scale (µ ) for this sample was set to the µ CMMPS = i=t,t,b,b E 1/4 T,i [57,60], and the factorisation (µ ) and resummation (µ ) scales to H T /2 = 1 2 i=t,t,b,b E T,i . A first type of reweighting is performed in the tt + ≥1b subcategories, using a method similar to the one outlined in Ref. [61]. The reweighting corrects the relative normalisation of the tt + ≥1b subcategories to match the predictions from S , while keeping the overall tt + ≥1b normalisation unchanged. A second type of reweighting is derived and performed on several kinematic variables sequentially. First the p T of the tt system is reweighted, and secondly the p T of the top quarks. The final reweighting is performed depending on the type of tt + ≥1b events. If there is only one additional HF jet, the p T of that jet is used in the final reweighting. If there is more than one additional HF jet, first the ∆R between the HF jets is reweighted and then the p T of the HF dijet system. A closure test is performed on all reweighted kinematic variables, showing a reasonable level of agreement between the reweighted P +P 8 sample and the S sample.
The P -B v1 generator was used to produce the samples of Wt single-top-quark backgrounds, with the CT10 PDF set. Overlaps between the tt and Wt final states were handled using the 'diagram removal' scheme [62]. The t-channel single-top-quark events were generated using the P -B v1 generator with the 4FS for the NLO matrix element calculations and the fixed 4FS PDF set CT10F4. The top quarks were decayed with M S [63], which preserves the spin correlations. The samples were interfaced to P 6.428 [64] with the P 2012 UE tune [65]. The single-top-quark Wt and t-channel samples were normalised to the approximate NNLO (aNNLO) theoretical cross-section [66][67][68].
Samples of W/Z+jets events were generated using S 2.2.1 [56]. Matrix elements were calculated for up to 2 partons at NLO and 4 partons at LO using C [69] and O L and merged with the S parton shower [70] using the ME+PS@NLO prescription [71]. The NNPDF3.0NNLO PDF set was used together with a dedicated parton shower tune developed by the S authors. The W/Z+jets events were normalised to the NNLO cross-sections [72][73][74][75][76].
Samples of ttV (V = W, Z) events were generated at NLO in the matrix elements calculation using MG5_ MC with the NNPDF3.0NLO PDF set interfaced to P 8.210 with the A14 UE tune. The ttH process was modelled using MG5_ MC with NLO matrix elements, NNPDF3.0NLO PDF set and factorisation and renormalisation scales set to µ = µ = m T /2, where m T is defined as the scalar sum of the transverse masses m T = p 2 T + m 2 of all final-state particles. The events were interfaced to P 8.210 with the A14 UE tune. Variations in ttH production due to the extended Higgs sector are not considered in this analysis. Measurements of the ttH production cross-section have been performed by the ATLAS and CMS Collaborations [77,78] at 13 TeV, both in agreement with the SM prediction within 25%.
The minor tH + X backgrounds, consisting of the production of a single top quark in association with a Higgs boson and jets (tH jb), and the production of a single top quark, a W boson and a Higgs boson (WtH), are treated as one background. The tH jb process was simulated with MG5_ MC interfaced to P 8.210 and the CT10 PDF set, and WtH was modelled with MG5_ MC interfaced to Herwig++ [79] using the CTEQ6L1 PDF set [80]. Additional minor SM backgrounds (diboson production, single top s-channel, t Z, tW Z, 4t, ttWW) were also simulated and accounted for, even though they contribute less than 1% in any analysis region.
Except where otherwise stated, all simulated event samples were produced using the full ATLAS detector simulation [81] based on G 4 [82]. Additional pile-up interactions were simulated with P 8.186 using the A2 set of tuned parameters [83] and the MSTW2008LO PDF set [84], and overlaid onto the simulated hard-scatter event. All simulated samples were reweighted such that the average number of interactions per bunch crossing (pile-up) matches that of the data. In the simulation, the top quark mass was set to m t = 172.5 GeV. Decays of band c-hadrons were performed by E G v1.2.0, except in samples simulated by the S event generator.
The samples and their basic generation parameters are summarised in Table 1.

Object and event selection
The data used in this analysis were recorded in 2015 and 2016 from √ s = 13 TeV pp collisions with an integrated luminosity of 36.1 fb −1 . Only runs with stable colliding beams and in which all relevant detector components were functional are used. Events are required to have at least one reconstructed vertex with two or more tracks with p T > 0.4 GeV. The vertex with the largest sum of the squared p T of associated tracks is taken as the primary vertex.
Events were recorded using single-lepton triggers, in both the +jets and final states. To maximise the event selection efficiency, multiple triggers were used, with either low p T thresholds and lepton identification and isolation requirements, or with higher p T thresholds but looser identification criteria and no isolation requirements. Slightly different sets of triggers were used for 2015 and 2016 data. For muons, the lowest p T threshold was 20 (26) GeV in 2015 (2016), while for electrons, triggers with a p T threshold of 24 (26) GeV were used. Simulated events were also required to satisfy the trigger criteria.
Electrons are reconstructed from energy clusters in the EM calorimeter associated with tracks reconstructed in the ID [85]. Candidates in the calorimeter transition region 1.37 < |η cluster | < 1.52 are excluded. Electrons are required to satisfy the tight identification criterion described in Ref.
[85], based on showershape and track-matching variables. Muons are reconstructed from track segments in the MS that are matched to tracks in the ID [86]. Tracks are then re-fit using information from both detector systems. The medium identification criterion described in Ref.
[86] is used to select muons. To reduce the contribution of leptons from hadronic decays (non-prompt leptons), both the electrons and muons must satisfy isolation criteria. These criteria include both track and calorimeter information, and have an efficiency of 90% for leptons with a p T of 25 GeV, rising to 99% above 60 GeV, as measured in Z → ee [85] and Z → µµ [86] samples. Finally, the lepton tracks must point to the primary vertex of the event: the longitudinal impact parameter z 0 must satisfy |z 0 sinθ| < 0.5 mm, while the transverse impact parameter significance must satisfy, |d 0 |/σ(|d 0 |) < 5 (3) for electrons (muons).
Jets are reconstructed from three-dimensional topological energy clusters [87] in the calorimeter using the anti-k t jet algorithm [55,88] with a radius parameter of 0.4. Each topological cluster is calibrated to the EM scale response prior to jet reconstruction. The reconstructed jets are then calibrated to the jet energy scale (JES) derived from simulation and in situ corrections based on √ s = 13 TeV data [89]. After energy calibration, jets are required to have p T > 25 GeV and |η| < 2.5. Quality criteria are imposed to identify jets arising from non-collision sources or detector noise, and events containing any such jets are removed [90]. Finally, to reduce the effect of pile-up an additional requirement using information about the tracks and the primary vertex associated to a jet (Jet Vertex Tagger) [91] is applied for jets with p T < 60 GeV and |η| < 2.4.
Jets are identified as containing the decay of a b-hadron (b-tagged) via an algorithm using multivariate techniques to combine information from the impact parameters of displaced tracks with the topological properties of secondary and tertiary decay vertices reconstructed within the jet [92,93]. Jets are b-tagged by directly requiring the output discriminant of the b-tagging algorithm to be above a threshold. A criterion with an efficiency of 70% for b-jets in tt events is used to determine the b-jet multiplicity for all final states and H + masses. For this working point, the c-jet and light-jet rejection factors are 12 and 381, respectively. For m H + ≤ 300 GeV, five exclusive efficiency bins are defined using the same b-tagging discriminant: 0-60%, 60-70%, 70-77%, 77-85% and 85-100%, following the procedure described in Ref. [94]. These step-wise efficiencies are used as input to the kinematic discriminant described in Section 5. When 'a b-tagged jet' is mentioned without any further specification, an efficiency of 70% is implied.
To avoid counting a single detector response as two objects, an overlap removal procedure is used. First, the closest jet within ∆R y = 0.2 of a selected electron is removed. If the nearest jet surviving this selection is within ∆R y = 0.4 of the electron, the electron is discarded, to ensure it is sufficiently separated from nearby jet activity. Muons are removed if they are separated from the nearest jet by ∆R y < 0.4, to reduce the background from muons from HF decays inside jets. However, if this jet has fewer than three associated tracks, the muon is kept and the jet is removed instead; this avoids an inefficiency for high-energy muons undergoing significant energy loss in the calorimeter.
The missing transverse momentum in the event is defined as the negative vector sum of the p T of all the selected electrons, muons and jets described above, with an extra term added to account for energy in the event that is not associated with any of these. This extra term, referred to as the 'soft term' in the following, is calculated from ID tracks matched to the primary vertex to make it resilient to pile-up contamination [95][96][97]. The missing transverse momentum is not used for event selection but is an input to the multivariate discriminants.
Events are required to have at least one electron or muon. The leading lepton must be matched to a lepton with the same flavour reconstructed by the trigger algorithm within ∆R < 0.15, and have a p T > 27 GeV. Additional leptons are required to have p T > 10 GeV, or > 15 GeV for events with two electrons. The latter requirement reduces the bakground due to jets and photons that are misidentified as electrons. Events in the +jets channel and the channel are required to be mutually exclusive. Electrons or muons from τ decays are also included in the analysis.
For the +jets channel, five or more jets, of which at least two jets have to be b-tagged, are required. For the channel, events with two leptons with opposite charge are selected, and at least three jets are required, of which two or more must be b-tagged. In the ee and µµ channels, the dilepton invariant mass must be > 15 GeV and outside the Z boson mass window of 83-99 GeV.

Analysis strategy
After the event selection, the samples in both the and the +jets final states contain mostly tt events. Events passing the event selection are categorised into separate regions according to the number of reconstructed jets and b-tagged jets. The regions where tbH + is enhanced relative to the backgrounds are referred to as signal regions (SRs), whereas the remaining regions are referred to as control regions (CRs).
In the SRs, for each H + mass hypothesis a different discriminating variable based on boosted decision trees (BDTs) is defined. In order to separate the H + signal from the SM background, the binned output of this variable is used together with the total event yields in the CRs in a combined profile likelihood fit. The fit simultaneously determines both the signal and background yields, while constraining the overall background model within the assigned systematic uncertainties. The event yields in the CRs are used to constrain the background normalisation and systematic uncertainties. The profile likelihood fit, including the treatment of backgrounds in the fit, is described in detail in Section 7.
For the +jets final state, two CRs (5j2b and ≥6j2b)4 and four SRs (5j3b, 5j≥4b, ≥6j3b and ≥6j≥4b) are defined, while in the final state, two CRs (3j2b and ≥4j2b) and two SRs (≥4j3b and ≥4j≥4b) are defined for all mass hypotheses. In addition, for the final state, the region with three b-tagged jets and no other jets (3j3b) is considered a SR for m H + < 1 TeV and a CR for m H + ≥ 1 TeV due to the change in expected signal yield for the different H + mass hypotheses.
The background from processes with prompt leptons is estimated using the simulated event samples described in Section 3. For tt production, the number of events with high leading jet p T is overestimated in the simulation, and a reweighting function for the leading jet p T distribution is determined by comparing simulation with data in a +jets CR that requires exactly four jets and at least two b-tagged jets. This function is validated in the dilepton channel and applied to both channels.
The normalisation of the Z+HF jets backgrounds is corrected by a factor of 1.3, extracted from dedicated control regions in data, defined by requiring two opposite-charge same-flavour leptons (e + e − or µ + µ − ) with an invariant mass compatible with the Z boson mass, 83 GeV < m < 99 GeV.
Processes that do not contain prompt electrons or muons from W or Z boson decays (multi-jet events) can still satisfy the selection criteria if they contain non-prompt leptons. The leading sources of non-prompt leptons are from semileptonic hadron decays or jets misidentified as leptons ('fake' leptons). These backgrounds are estimated using data. For the +jets final state a matrix method [98] is employed. An event sample that is enriched in non-prompt leptons is selected by using looser isolation or identification requirements for the lepton. These events are then weighted according to the efficiencies for both the prompt and non-prompt leptons to pass the tighter default selection. These efficiencies are measured using data in dedicated CRs. In the final state, this background is estimated from simulations, and the normalisation is determined by comparing data and simulations in a CR of same-sign dilepton events. The contribution of multi-jet events to the final state is found to be negligible.
The expected event yields of all SM processes and the number of events observed in the data are shown in Figure 2 for the and the +jets final states before performing the fit to data. The expected H + signal yields for m H + = 200 GeV, assuming a cross-section times branching ratio of 1 pb, are also shown.
The training of the BDTs that are used to discriminate signal from background in the SRs is performed with the TMVA toolkit [99]. BDTs are trained separately for each value of the 18 generated H + masses and for each SR against all the backgrounds ( +jets channel) or the tt background ( channel). For the BDT training in the +jets channel, the SRs 5j3b and 5j≥4b are treated as one region, in order to increase the number of simulated events available for training.
The BDT variables include various kinematic quantities with the optimal discrimination against the tt + ≥1b background. For H + masses above 400 GeV the most important variables in the +jets final state Data / Pred.
Non-t Uncertainty Figure 2: Comparison of predicted and observed event yields. Each background process is normalised according to its cross-section and the prediction has not been fitted to the data. The tt + X includes contributions from ttW, tt Z and ttH. A signal with m H + = 200 GeV, normalised to a cross-section times branching ratio for H + → tb of 1 pb, is shown as a dashed line. The lower panel displays the ratio of the data to the total prediction. The hatched bands show uncertainties before the fit to the data, which are dominated by systematic uncertainties as discussed in Section 6. The comparison is shown for all signal and control regions used in the analysis. For the final state: CR 3j2b, CR/SR 3j3b, CR ≥4j2b, SR ≥4j3b, SR ≥4j≥4b. For the +jets final state: CR 5j2b, SR 5j3b, SR 5j≥4b, CR ≥6j2b, SR ≥6j3b, SR ≥6j≥4b.
are the scalar sum of the p T of all jets, H jets T , and the leading jet p T . For a mass at or below 300 GeV, a kinematic discriminant, D, as described below, is used as an input variable for the BDT. The kinematic discriminant, D, and the invariant mass of the pair of jets that are not b-tagged and have the smallest ∆R are the most important variables in the low mass range. The latter variable is not used in the 5j≥4b SR, where it is not well defined.
The kinematic discriminant, D, is a variable reflecting the probability that an event is compatible with the H + → tb and the tt hypotheses, and is defined as D = P H + (x)/(P H + (x) + P tt (x)), where P H + (x) and P tt (x) are probability density functions for x under the signal hypothesis and background (tt) hypothesis, respectively. Here, the event variable x indicates the set of the missing transverse momentum and the four-momenta of reconstructed electrons, muons and jets.
The probability P H + (x) is defined as the product of the probability density functions for each of the reconstructed invariant masses in the event: • the mass of the semileptonically decaying top quark, m b ν , • the mass of the hadronically decaying W boson, m q 1 q 2 , • the difference between the masses of the hadronically decaying top quark and the hadronically decaying W boson m b h q 1 q 2 − m q 1 q 2 , and • the difference between the mass of the charged Higgs boson and the mass of the leptonically or hadronically decaying top quark, depending on whether the top quark from the charged Higgs boson decays leptonically or hadronically.
In this context q 1 or q 2 refer to the quarks from the W boson decay, and ν to the lepton and neutrino from the other W boson decay, b h to the b-quark from the hadronic top quark decay, b to the b-quark from the leptonic top quark decay and b H + to the b-quark directly from the H + decay. The probability P tt (x) is constructed from probability density functions obtained from simulated tt events. For the SRs with five jets, P tt (x) is defined using the same invariant masses as above. The jet that does not originate from a top quark decay is used instead of b H + . For the SRs with at least six jets the power of the discriminant is improved by using the invariant mass of the two highest-p T jets not originating from the hadronisation of The functional form of the probability density functions is obtained from simulation using the reconstructed masses of jets and leptons matched to simulated partons and leptons for H + and tt. The neutrino fourmomentum is derived with the assumption that the missing transverse momentum is solely due to the neutrino; the constraint m 2 W = (p + p ν ) 2 is used to obtain p ν,z . If two real solutions exist, they are sorted according to the absolute value of their p z , i.e., |p z,v 1 | < |p z,v 2 |. In approximately 60% of the cases p z,v 1 is closer than p z,v 2 to the generator-level neutrino p z . Two different probability density functions are constructed, one for each solution, and the probability is defined as a weighted average of the two probability density functions. The weight is taken as the fraction of the corresponding solution being closer to the generated neutrino p z . Also, if no real solution exists, the p x and p y components are scaled by a common factor until the discriminant of the quadratic equation is exactly zero, yielding only one solution.
When evaluating P H + (x) and P tt (x) for the calculation of D, all possible parton-jet assignments are considered since the partonic origin of the jets is not known. In order to suppress the impact from parton-jet assignments that are inconsistent with the correct parton flavours, a weighted average over all parton-jet assignments is used. The value of P H + (x) and P tt (x) for each parton-jet assignment is weighted with a probability based on the b-tagging discriminant value of each jet. The distribution of the step-wise efficiencies of the b-tagging algorithm, as described in Section 4, is used as a probability density function, with the b-jet hypothesis for generated b-quarks and the light-jet hypothesis for other generated partons. Due to the large number of events in which q 1 and q 2 cannot be matched to different jets, the average of two different probability density functions, where either all partons can be matched to jets or only one jet can be matched to q 1 and q 2 , is used. This discriminant gives better background suppression than would be obtained by adding the kinematic input variables directly to the BDT.
In the final state, approximately ten optimal kinematic variables from the analysis objects and their combinations were selected for each SR, independently for the low-mass region (m H + ≤ 600 GeV) and the high-mass region (m H + > 600 GeV). For the high-mass region, the most important variables are the scalar sum of the p T of all jets and leptons, H all T , and the transverse momentum of the jet pair with maximum p T . For the low-mass region, the smallest invariant mass formed by two b-tagged jets and the smallest invariant mass formed by a lepton and a b-tagged jet, are among the most important variables.
All BDT input variables in the +jets and final states are listed in the Appendix. In most regions, the distributions show a reasonable level of agreement between simulation and data within the systematic and statistical uncertainties before the fit to the data (pre-fit). As examples, Figures 3 and 4 show the distribution of the observed and pre-fit expected event yields for H jets T in the +jets channel and H all T in the channel. Figure 5 shows the expected BDT output distributions, normalised to unity, for selected H + signal samples and the background processes in the SRs.     Figure 3: Distributions of the H jets T variable before the fit to the data in the four SRs of the +jets channel: (a) 5j3b, (b) ≥6j3b, (c) 5j≥4b, (d) ≥6j≥4b. Each background process is normalised according to its cross-section and the normalisation of the tt + ≥1b and tt + ≥1c backgrounds corresponds to the prediction from P +P 8 for the fraction of each of these components relative to the total tt prediction. The tt + X includes contributions from ttW, tt Z and ttH. In addition, the expectation for a 200 GeV signal is shown for a cross-section times branching ratio of 1 pb. The lower panels display the ratio of the data to the total prediction. The hatched bands show the pre-fit uncertainties. The level of agreement is improved post-fit due to the adjustment of the normalisation of the tt + ≥1b and tt + ≥1c backgrounds and the other nuisance parameters by the fit.   Figure 4: Distributions of the H all T variable before the fit to the data in the three SRs of the channel: (a) 3j3b, (b) ≥4j3b and (c) ≥4j≥4b. Each background process is normalised according to its cross-section and the normalisation of the tt + ≥1b and tt + ≥1c backgrounds corresponds to the prediction from P +P 8 for the fraction of each of these components relative to the total tt prediction. The tt + X includes contributions from ttW, tt Z and ttH. In addition, the expectation for a 200 GeV signal is shown for a cross-section times branching ratio of 1 pb. The lower panels display the ratio of the data to the total prediction. The hatched bands show the pre-fit uncertainties. The level of agreement is improved post-fit due to the adjustment of the normalisation of the tt + ≥1b and tt + ≥1c backgrounds and the other nuisance parameters by the fit. (d)

Systematic uncertainties
Systematic uncertainties from various sources affect this search, such as uncertainties in the luminosity measurement, the reconstruction and calibration of physics objects, in particular b-tagged jets, and the modelling of the signal and background processes. Uncertainties can either modify the normalisation of the signal and background processes, change the shape of the final distributions, or both. An overview of the systematic uncertainties is given in Table 2. The impact of the systematic uncertainties is listed in Table 5 in Section 8. The combined uncertainty in the integrated luminosity for the data collected in 2015 and 2016 is 2.1%, and it is applied as a normalisation uncertainty for all processes estimated using simulation. It is derived, following a methodology similar to that detailed in Ref. In the reconstruction of quantities used for the BDT, E miss T is used. The E miss T calculation depends on the reconstruction of leptons and jets, and the E miss T uncertainties are therefore related to the uncertainties associated with these objects, which are propagated to the E miss T uncertainty estimation. Uncertainties due to soft objects (not included in the calculation of the leptons and jets) are also considered [96].
Differences between data and simulation in the b-tagging efficiency for b-jets, c-jets and light jets are taken into account using correction factors. For b-jets, the corrections are derived from tt events with final states containing two leptons, and the corrections are consistent with unity within uncertainties at the level of a few percent over most of the jet p T range. The mis-tag rate for c-jets is also measured in tt events, identifying hadronic decays of W bosons including c-jets. For light jets, the mis-tag rate is measured in multi-jet events using jets containing secondary vertices and tracks with impact parameters consistent with a negative lifetime. Systematic uncertainties affecting the correction factors are derived in the p T and η bins used for extracting the correction factors. They are transformed into uncorrelated components using an eigenvector decomposition, taking into account the bin-to-bin correlations [92, 93, 103]. For m H + > 300 GeV, corrections corresponding to the fixed working point of 70% efficiency are used and a total of 6, 3 and 16 independent uncorrelated eigen-variations are considered as systematic uncertainties for b-, cand light jets, respectively. For m H + ≤ 300 GeV, corrections for the step-wise efficiencies are used to support the kinematic discriminant D and the number of eigen-variations is increased by a factor of five to account for the five b-tagging efficiency bins. In addition, uncertainties due to tagging the hadronic decays of τ-leptons as b-jets are considered. For m H + > 300 GeV, an additional uncertainty is included due to the extrapolation of scale factors for jets with p T > 300 GeV, beyond the kinematic reach of the data calibration samples used [93].
The modelling uncertainty of the H + signal is estimated by varying the renormalisation and factorisation scales up and down by a factor of two. The uncertainty ranges from 7% at low masses to 15% at masses above 1300 GeV for the +jets final state, and from 12% to 16.5% for the final state. Additionally, the PDF uncertainty in the modelling is estimated using the PDF4LHC15_30 PDF set [104], which is based on a combination of the CT14 [105], MMHT14 [106] and NNPDF3.0 [44] PDF sets and contains 30 components obtained using the Hessian reduction method [107][108][109].
The modelling of the tt + jets background is one of the largest sources of uncertainty in the analysis and many different components are considered. The uncertainty in the inclusive tt production cross-section at NNLO+NNLL [47] is 6%, including effects from varying the factorisation and renormalisation scales, the PDF, the QCD coupling constant α s , and the top quark mass. Due to the large difference between the 4FS prediction and the various 5FS predictions for the tt + ≥3b process, an additional 50% normalisation uncertainty is assigned to this background.
The uncertainty due to the choice of NLO generator is derived by comparing the nominal P sample with a sample generated using S 2.2.1 with a 5FS PDF. A P sample with the same settings as in the nominal P +P 8 sample, but using H 7 [79,110] for parton showering, is used to assess the uncertainty due to the choice of parton shower and hadronisation model. Furthermore, the uncertainty due to the modelling of initial-and final-state radiation is evaluated with two different P +P 8 samples in which the radiation is increased or decreased by halving or doubling the renormalisation and factorisation scales in addition to simultaneous changes to the h damp parameter and the A14 tune parameters [111].
For the tt + ≥1b background, an additional uncertainty is assigned by comparing the predictions from P +P 8 and S with 4FS. This takes into account the difference between a 5FS inclusive tt prediction at NLO and a 4FS NLO ttbb prediction. For the tt + ≥1c background, an additional uncertainty is derived by comparing a MG5_ MC sample that is interfaced to Herwig++ [79] with the nominal event sample. In this MG5_ MC event sample, a three-flavour scheme is employed and the ttcc process is generated at the matrix element level [112] using the CT10F3 PDF set, while in the nominal sample the charm jets are primarily produced in the parton shower. All of these uncertainties, with the exception of the inclusive and tt + ≥ 3b cross-sections, are considered to be uncorrelated amongst the tt + ≥ 1b, tt + ≥1c, and tt + light samples. For the modelling of the tt + ≥1b backgrounds, the alternative samples are reweighted to the NLO prediction of ttbb from S before the uncertainty is evaluated.
In addition, uncertainties due to the reweighting to the S NLO prediction of ttbb are considered. For these uncertainties, the tt + ≥1b is reweighted to different S predictions with modified scale parameters, in particular where the renormalisation scale is varied up and down by a factor of two, where the functional form of the resummation scale is changed to µ CMMPS and where a global scale choice µ = µ = µ = µ CMMPS is used. Two alternative PDF sets, MSTW2008NLO [84] and NNPDF2.3NLO [44], are used, and uncertainties in the underlying event and parton shower are estimated from samples with an alternative set of tuned parameters for the underlying event and an alternative shower recoil scheme. Due to the absence of b-jets from multi-parton interactions and final-state gluon radiation in the ttbb prediction from S , a 50% uncertainty is assigned to the tt + b (MPI/FSR) category based on studies of different sets of UE tunes. An uncertainty due to the reweighting of the leading jet p T is determined by comparing a reweighted event sample with an event sample without reweighting. Because the reweighting changes the normalisation for jet p T > 400 GeV by 15%, an additional normalisation uncertainty of 15% is applied in this region. The reweighting factors are derived from the CR with exactly four jets and at least two b-tagged jets and applied to higher jet multiplicity bins. However, the effect of this extrapolation is expected to be small and is covered by the above uncertainties.
An uncertainty of 5% is assigned to the total cross-section for single top-quark production [66][67][68], uncorrelated between Wt and t-channel production. An additional uncertainty due to initial-and finalstate radiation is estimated using samples with factorisation and renormalisation scale variations and appropriate variations of the Perugia 2012 set of tuned parameters. The parton showering and hadronisation modelling uncertainties in the single-top Wt and t-channel production are estimated by comparing with samples where the parton shower generator is Herwig++ instead of P 6.428. The uncertainty in the interference between Wt and tt production at NLO [62] is assessed by comparing the default 'diagram removal' scheme with an alternative 'diagram subtraction' scheme [62,113].
The uncertainty arising from ttV generation is estimated by comparison with samples generated with S . The uncertainty in the ttV production cross-section is about 15%, taken from the NLO predictions [15,[114][115][116], treated as uncorrelated between ttW and tt Z with PDF and QCD scale variations.
The ttH modelling uncertainty is assessed through an uncertainty in the cross-section, uncorrelated between QCD ( +5.8 −9.2 %) and the PDFs (±3.6%) [15,[117][118][119][120][121], and the modelling of the parton shower and hadronisation by comparing P 8 with Herwig++. The minor tH + X backgrounds, tH jb and WtH are treated as one background and its cross-section uncertainty is 6% due to PDF uncertainties and another 10% due to factorisation and renormalisation scale uncertainties [15].
The uncertainties from the data-driven estimation of non-prompt leptons are based on a comparison between data and the non-prompt lepton estimates in CRs. A 50% uncertainty is assigned in the +jets final state. In the final state, where all backgrounds with one or no prompt leptons fall into this category, including W+jets and single top production, an uncertainty of 25% is assigned.
An uncertainty of 40% is assumed for the W+jets cross-section, uncorrelated between jet bins, with an additional 30% for W+HF jets, uncorrelated for two, three and more than three HF jets. These uncertainties are derived from variations of the renormalisation and factorisation scales and matching parameters in S simulations. An uncertainty in Z+jets of 35% is applied, uncorrelated among jet bins in the final state. This uncertainty accounts for both the variation of the scales and matching parameters in S simulations and the data-driven correction factors applied to the Z+HF jets component. In the final state, only the Z+jets component is estimated separately, and the W+jets background is included in the estimation of the background from non-prompt leptons.

Statistical analysis
In order to test for the presence of an H + signal, a binned maximum-likelihood fit to the data is performed simultaneously in all categories, and each mass hypothesis is tested separately. The inputs to the fit include the number of events in the CRs and the binned BDT output in the SRs. Two initially unconstrained fit parameters are used to model the normalisation of the tt + ≥1b and tt + ≥1c backgrounds. The procedures used to quantify the level of agreement with the background-only or background-plus-signal hypothesis and to determine exclusion limits are based on the profile likelihood ratio test and the CL s method [122][123][124]. The parameter of interest is the signal strength, µ, defined as the product of the production cross-section σ(pp → tbH + ) and the branching ratio B(H + → tb).
To estimate the signal strength, a likelihood function, L(µ, θ), is constructed as the product of Poisson probability terms. One Poisson term is included for every CR and every bin of the BDT distribution in the SRs. The expected number of events in the Poisson terms is a function of µ, and a set of nuisance parameters, θ. The nuisance parameters encode effects from the normalisation of backgrounds, including two free normalisation factors for the tt + ≥1b and tt + ≥1c backgrounds, the systematic uncertainties and one parameter per bin to model statistical uncertainties in the simulated samples. All nuisance parameters are constrained with Gaussian or log-normal terms. There are about 170 nuisance parameters considered in the fit, the number varying slightly across the range of mass hypotheses.
To extract the exclusion limit on µ = σ(pp → H + ) × B(H + → tb), the following test statistic is used: The values of the signal strength and nuisance parameters that maximise the likelihood function are represented byμ andθ, respectively. For a given value of µ, the values of the nuisance parameters that maximise the likelihood function are represented byθ(µ). Tables 3 and 4 show the post-fit event yields under the background-plus-signal hypothesis for a signal mass m H + = 200 GeV. A value of σ(pp → tbH + ) × B(H + → tb) = −0.36 pb is obtained from the fit. The corresponding post-fit distributions of the BDT discriminant in the SRs are shown in Figures 6 and 7 for a 200 GeV H + mass hypotheses for the +jets and final state, respectively. Table 3: Event yields of the SM background processes and data in all categories of the +jets final state, after the fit to the data under the background-plus-signal hypothesis (m H + = 200 GeV). The expected event yields for the H + signal masses of 200 GeV and 800 GeV are shown with pre-fit uncertainties and assuming a cross-section times branching ratio of 1 pb. The quoted uncertainties include both the statistical and systematic components. The uncertainties take into account correlations and constraints of the nuisance parameters. 'Other top' includes contributions from Zt as well as sand t-channel production.  A summary of the systematic uncertainties is given in Table 5. Depending on the particular H + mass hypothesis, the total systematic uncertainty is dominated by the uncertainties in the modelling of the tt + ≥1b background, the jet flavour-tagging uncertainties and the uncertainties due to the limited size of simulated event samples.   The 95% confidence level (CL) upper limits on σ(pp → H + ) × B(H + → tb) using the CL s method are presented in Figure 8. The observed (expected) 95% CL upper limits on the pp → tbH + production cross-section times the branching ratio B(H + → tb) range from σ × B = 2.9 (3.0) pb at m H + = 200 GeV to σ × B = 0.070 (0.077) pb at m H + = 2 TeV. The compatibility of the SM hypothesis with the results obtained from the fit to the data is tested. The largest deviation from the SM hypothesis is observed at 300 GeV. Given that a negativeμ is observed under this mass hypothesis, the test statistic t 0 = −2 ln (L(0,θ(0))/L(μ,θ)) is used to quantify the deviation of the fitted result from the SM expectation. A local p 0 value of 1.13% is obtained at 300 GeV, corresponding to the probability to obtain a deviation at least as large as the one observed in data provided that only SM processes are present. In comparison with a previous search for t[b]H + production followed by H + → tb decays [20], more stringent limits on H + masses for particular models and parameter choices can be set. The analysis reach is increased and now also includes H + masses between 600 GeV and 2 TeV. The excluded region of parameter space for the model-dependent interpretation is extended significantly for low tanβ and an additional excluded region is added at high tanβ.

Results
The ATLAS Collaboration has also set limits on the production of H + using the H + → τν decay with the same data [30]. The τν final state can be used to set limits at high tanβ which are more stringent than those from the tb final state, and to probe H + masses below 200 GeV, in both the m mod− h and hMSSM scenarios. Figure 10 shows a superposition of the limits from the two final states, where the limits from the τν final state exclude a larger portion of the parameter space at high tanβ and low H + masses than the tb limits alone.
Non-t Uncertainty (d) Figure 6: Distributions of the BDT output after the fit to the data in the four SRs of the +jets final state: (a) 5j3b, (b) ≥6j3b, (c) 5j≥4b and (d) ≥6j≥4b for the 200 GeV mass hypothesis. Each background process is normalised according to its post-fit cross-section. The tt + X includes contributions from ttW, tt Z and ttH. The total prediction of the BDT distributions includes cases where the signal obtained from the fit is negative. For this particular mass point the fitted signal strength is µ = −0.4 ± 1.5 pb. The pre-fit signal distribution is shown superimposed as a dashed line with arbitrary normalisation. The lower panels display the ratio of the data to the total prediction. The hatched bands show the post-fit uncertainties.    Additionally, taking into consideration the H + → τν decay, even stricter exclusions can be made at high tanβ and low H + masses. In the context of the hMSSM, the H + mass range up to 1100 GeV is excluded at tanβ = 60, and all tanβ values are excluded for m H + below 160 GeV.
(Taiwan), RAL (UK) and BNL (USA), the Tier-2 facilities worldwide and large non-WLCG resource providers. Major contributors of computing resources are listed in Ref. [127]. Table 6: Input variables to the classification BDT in the +jets and channels. The symbols j, b, u, and E miss T represent the four-momenta of jets, b-tagged jets, non-b-tagged jets, the lepton and the missing transverse momentum. All numbered indices refer to ordering in transverse momentum, with 1 as leading. The SRs where the variables are used are indicated for the channel. In the +jets channels the variables are used for all channels. In the +jets channel the discriminant, D, is only used for charged Higgs boson masses m H + ≤ 300 GeV. For the channel a very large set of kinematic variables using combinations of the analysis objects was examined, and approximately ten optimal variables were selected for each SR independently for the low-mass region (m H + ≤ 600 GeV) and the high-mass region (m H + > 600 GeV). Energy difference between the third jet and the subleading lepton Energy of third jet ∆m(j 1 + j 2 , j 1 + j 3 + 2 + E miss T ) Inv. mass difference between j 1 + j 2 and j 1 + j 3 + 2 + E miss T ∆R(j 2 , j 1 + 2 + E miss T ) Angular difference between subleading jet and j 1 + 2 + E miss p T of the pair of lepton and b-tagged jet with largest ∆η m(( , b) ∆φ min ) Inv. mass of the pair of lepton and b-tagged jet with smallest ∆φ ∆E(b 1 , 1 + E miss T ) Energy difference between the leading b-tagged jet and 1 + E miss T ∆m(j 2 + j 3 , j 1 + 1 + 2 ) Inv. mass difference between j 2 + j 3 and j 1 + 1 + 2 ∆m( 1 + j 3 + E miss T , j 1 + j 2 + 2 ) Inv. mass difference between 1 + j 3 + E miss T and j 1 + j 2 + 2 ∆p T (j 1 , j 3 ) p T difference between leading and third jet m min (b-pair) Smallest invariant mass of any b-tagged jet pair m min ( , b) Smallest invariant mass of any pair of lepton and b-tagged jet p T (b 2 + 1 + 2 + E miss T ) p T of b 2 + 1 + 2 + E miss T ∆R( 2 , j 2 + j 3 + 1 + E miss T ) Angular difference between 2 and j 2 + j 3 + 1 + E miss p T difference between leading and third jets ∆m(j 2 + 1 + E miss T , j 1 + j 3 + 1 ) Inv. mass difference between j 2 + 1 + E miss T and j 1 + j 3 + 1 p T (( , b) ∆R min ) p T of the pair of lepton and b-tagged jet with smallest ∆R m(j-pair ∆η min ) Inv. mass of the jet pair with smallest ∆η ∆p T (j 1 , j 2 + E miss T ) p T difference between leading jet and j 2 + E miss Energy difference between 1 + E miss T and j 1 + j 2 E(j 1 ) Energy of the leading jet p T difference between subleading b-tagged jet and b 1 + 2 ∆p T (j 2 , j 3 + 1 + E miss T ) p T difference between subleading jet and j 3 + 1 + E miss T ∆E(j 3 , j 2 + 1 + 2 + E miss T ) Energy difference between third jet and j 2 + 1 + 2 + E miss T ∆m(j 2 + 2 + E miss T , j 1 + 2 + E miss T ) Inv. mass difference between j 2 + 2 + E miss    [27] LHCb Collaboration, Test of lepton flavor universality by the measurement of the B 0 → D * − τ + ν τ branching fraction using three-prong τ decays, Phys.