Search for Higgs boson pair production in the four b quark ﬁnal state in proton-proton collisions at √ s = 13 TeV

A search for pairs of Higgs bosons produced via gluon and vector boson fusion is presented, focusing on the four b quark ﬁnal state. The data sample consists of proton-proton collisions at a center-of-mass energy of 13 TeV, collected with the CMS detector at the LHC, and corresponds to an integrated luminosity of 138 fb − 1 . No deviation from the background-only hypothesis is observed. A 95% conﬁdence level upper limit on the Higgs boson pair production cross section is observed at 3.9 times the standard model prediction for an expected value of 7.8. Constraints are also set on the modiﬁers of the Higgs ﬁeld self-coupling, κ λ , and of the coupling of two Higgs bosons to two vector bosons, κ 2V . The observed (expected) allowed intervals at the 95% conﬁdence level are − 2.3 < κ λ < 9.4 ( − 5.0 < κ λ < 12.0) and − 0.1 < κ 2V < 2.2 ( − 0.4 < κ 2V < 2.5). These are the most stringent observed constraints to date on the HH production cross section and on the κ 2V coupling.


1
The discovery of the Higgs boson (H) by the ATLAS and CMS Collaborations [1][2][3] proves the existence of a fundamental scalar sector of the standard model of particle physics (SM), but the experimental confirmation of the Brout-Englert-Higgs mechanism [4][5][6] requires the determination of the shape of the postulated scalar potential. This shape is governed by a parameter, λ, that drives the strength of the Higgs boson self-couplings and can thus be determined experimentally with a measurement of Higgs boson pair (HH) production.
At the CERN Large Hadron Collider (LHC), at the energy of √ s = 13 TeV, the dominant HH production mode in the SM is through the gluon fusion mechanism (ggF), with a cross section of 31.1 +2.1 −7.2 fb [7][8][9][10][11][12][13][14], followed by the vector boson fusion process (VBF), with a cross section of 1.726 ± 0.036 fb [15] and characterized by the presence of two additional hadronic jets, j, giving a bbbbjj final state. Variations of the Higgs boson self-coupling with respect to the SM prediction are parameterized by the modifier κ λ = λ/λ SM and affect the ggF and VBF production modes. The VBF production mode also depends on the strength of the interaction of pairs of vector bosons V (= W, Z) with a single (VVH) and a pair (VVHH) of Higgs bosons, whose values with respect to the SM prediction are parameterized by the modifiers κ V and κ 2V , respectively. Departures from the relation κ 2V = κ 2 V predicted in the Brout-Englert-Higgs mechanism are possible in models of physics beyond the SM where the Higgs boson is a composite state emerging from the presence of new strong dynamics at the TeV scale [16].
The ATLAS and CMS Collaborations have searched for ggF HH production with a data set corresponding to an integrated luminosity of about 36 fb −1 in a variety of final states [17][18][19][20][21][22][23][24][25][26], whose combinations [27, 28] set an observed (expected) upper limit at the 95% confidence level (CL) on the SM production cross section of 7 (10) and 22 (13) times theoretical prediction, respectively. Updated searches have been performed with an integrated luminosity of about 140 fb −1 [29-31], and the most stringent observed limits are from the ATLAS search in the bb γγ final state and correspond to 4.2 times the SM prediction, with a value of κ λ between −1.5 and 6.7 at the 95% CL. The VBF HH production process has been studied in the bbbb [32] and bb γγ [30] final states by the ATLAS and CMS Collaborations, respectively, and the most stringent observed constraints at the 95% CL correspond to −0.43 < κ 2V < 2.56 from bbbb and 225 times the SM cross section prediction from bb γγ.
This Letter reports on searches for both the ggF and VBF HH production mechanisms in the bbbb decay channel. In the SM, this decay mode is characterized by a combined branching fraction of 0.339 ± 0.008 for m H = 125 GeV [33]. The analysis uses a sample of proton-proton collision (pp) events at √ s = 13 TeV recorded between 2016 and 2018 with the CMS detector, corresponding to an integrated luminosity of 138 fb −1 .
Hadronic jets are clustered from the PF objects using the anti-k T algorithm [46, 47] with a distance parameter of 0.4. Jet energy corrections are derived from simulation studies and corrected with in situ measurements to match the energy scale in data and in simulation [44]. Jets originating from b quarks are identified using as a discriminant the output of a deep neural network algorithm (DEEPJET) [48,49], trained using as input information the properties of the PF con-stituents of the jets and of the secondary vertices associated with them. For the jets in this search, two working points (WPs) of the DEEPJET discriminant are considered: the medium WP, which yields a b jet identification efficiency of 75% with a corresponding misidentification rate of light flavor and gluon (charm) jets of about 1 (10)%, and the tight WP, which corresponds to a b jet identification efficiency of 58% and to a misidentification rate of about 0.1 (2)%.
Signal processes from ggF HH production are simulated at next-to-leading order (NLO) accuracy in quantum chromodynamics (QCD) with POWHEG 2.0 [50-52], and samples for VBF HH production are generated at leading order (LO) accuracy in QCD using MAD-GRAPH5 aMC@NLO 2.6.5 [53] for various combinations of couplings. The distributions are scaled by functions of the couplings defined according to the known dependence of the theoretical cross section [54] and added together to model arbitrary coupling combinations, and the total predictions are normalized to the corresponding next-to-NLO (NNLO) cross section [13] for ggF and by the ratio of the next-to-NNLO [15]  The trigger selection for events collected in 2016 requires the presence of four jets with transverse momentum p T > 45 GeV or of two jets with p T > 30 GeV and two jets with p T > 90 GeV. For 2017 (2018) data, the presence of four jets above the p T thresholds of 40, 45, 60, and 75 GeV is required together with H T > 300 (330) GeV, respectively, where H T denotes the scalar sum of the transverse momentum of the jets reconstructed in the event. As a consequence of the change in jet trigger thresholds, data collected in 2017 and 2018 are analyzed separately from those collected in 2016.
Offline, events are required to contain at least four jets with pseudorapidity |η| < 2.4 (2.5) and p T > 30 (40) GeV for the 2016 (2017-2018) data, respectively. Jets are required to satisfy the tight WP of the PF jet identification algorithm [57,58], corresponding to an efficiency larger than 99%. If their p T is below 50 GeV, the medium WP of the pileup discriminant [57] is also required, for a signal jet efficiency of about 90%. If more than four jets satisfy these criteria, the four objects with the largest DEEPJET output are selected. The p T of these four jets are corrected with a multivariate regression method developed for b jets that improves the determination of the momentum by up to 15% and simultaneously estimates the per jet resolution achieved [59]. After the application of this method, the resolution on the dijet invariant mass for H → bb events reconstructed in this analysis ranges between 11 and 14%. At least three of the selected jets are required to satisfy the medium WP of the DEEPJET discriminant.
Events are rejected if they contain an electron or a muon with p T > 15 and 10 GeV, respectively, and |η| < 2.4, where these objects must satisfy identification discriminants and criteria that include isolation and impact parameter with respect to the primary interaction vertex. This selection suppresses background events containing leptonic top quark decays.
The two Higgs boson candidates are formed by pairing the four jets. There are three possible pairings of jets, and in each the two Higgs boson candidates, denoted as H 1 and H 2 , are defined by the relation p T (H 1 ) > p T (H 2 ). The (H 1 , H 2 ) pairings are ordered according to the increasing value of a distance parameter d = |m H 1 − km H 2 |/ √ 1 + k 2 . The constant k is the ratio of the expected peak positions of the reconstructed Higgs boson masses for events that are correctly paired, k = c 1 /c 2 = (125 GeV)/(120 GeV) = 1.04. Its value differs from 1 because of the residual jet momentum dependence of the multivariate energy regression that more strongly impacts the softer H candidate. If the difference in the distance parameter of the first and second pairing, ∆d, is larger than 30 GeV, corresponding to about two times the resolution on the Higgs boson mass, the pairing with the smallest d is chosen. Conversely, if ∆d ≤ 30 GeV, the experimental resolution limits the capability to identify the correct pairing based on the invariant masses, and a choice is made between the first and second pairing as the one that maximizes the p T of the two Higgs boson candidates in the four-jet center-of-mass reference frame. This procedure results in a correct jet pairing of about 96% of the selected events in a ggF SM HH sample, and amounts to 82-96 (91-98)% for the different couplings studied in ggF (VBF) signal events.
The two non-b jets in the VBF production events are selected with p T > 25 GeV and |η| < 4.7, and they must satisfy the tight WP of the jet identification algorithm and the medium WP of the pileup discriminant if p T < 50 GeV. For the 2017 data, affected by large noise in the endcaps of the electromagnetic calorimeter (ECAL), jets in the region 2.6 < |η| < 3.1 are additionally required to satisfy the tight WP of the pileup discriminant to mitigate the noise effects. The two VBF jet candidates j 1 and j 2 are chosen as the highest p T jet and the second-highest p T jet that has an opposite η sign with respect to the former.
Events that do not contain such a VBF jet pair are assigned to the ggF category. About 26-28% of ggF events contain additional jets that satisfy the above requirements on the VBF jet candidates, and in order to correctly classify them, a boosted decision tree (BDT) discriminant is trained to separate ggF and VBF HH signal events. The discriminant uses p T (H 1 ), p T (H 2 ), p T (j 1 ), p T (j 2 ), the invariant mass and absolute value of pseudorapidity of the jj system, the angular separation ∆R = √ (∆η) 2 + (∆φ) 2 , where φ is the azimuthal angle, between the two H candidates and between each H and VBF jet, the absolute value of the polar angles with respect to the beam direction of the two VBF jets in the centerof-mass frame of the six selected jets and the product of the two Higgs boson centralities , where ∆η = η(j 1 ) − η(j 2 ) and η avg = [η(j 1 ) + η(j 2 )]/2. The discriminant is trained to separate the SM ggF HH signal from the κ 2V = 2 VBF signal, in order to optimize both the sensitivity to the anomalous κ 2V coupling hypotheses and the correct classification of SM ggF signal events. The value κ 2V = 2 is chosen because it is representative of the event kinematics in the presence of anomalous couplings, characterized by the large invariant mass of the jj and HH systems. These signals are associated with a large increase of the total cross section that would make them detectable with the available data set. A threshold on the BDT output is chosen to assign events to either the ggF or VBF category. It results in the correct assignment of about 97% of all ggF HH signal events to the ggF category and of about 60 (80)% of SM (κ 2V = 2) VBF HH events that contain the additional jets to the VBF category.
Events classified as ggF or VBF signal are further divided into subcategories to optimize the sensitivity of the search for anomalous coupling hypotheses. Events in the ggF category are divided into a low-and high-mass category if the reconstructed invariant mass of the HH system, m HH , is below or above 450 GeV, where the boundary is defined according to the kinematic properties of the signal. The latter category efficiently collects SM ggF HH events, while the former increases the acceptance to signals with anomalous κ λ values. Events in the VBF category are instead divided into a "SM-like" and an "anomalous κ 2V -like" category depending on the value of the discriminant trained to separate ggF and VBF. The categorization thresholds were chosen to maximize the expected sensitivity to VBF HH signals, and result in the assignment of about 25-30% of the VBF κ 2V = 2 events to the anomalous κ 2V -like category and 95% of SM VBF events to the SM-like category.
The large multijet background that originates from QCD and tt hadronic processes is estimated from the data using background-dominated regions. Analysis signal (A SR ) and control (A CR ) regions are defined by requiring χ < 25 GeV and 25 ≤ χ < 50 GeV, respectively, where χ is the distance from the expected peak position of the two Higgs boson candidates' invariant masses and is defined as where c 1 and c 2 are as defined for the pairing of the four jets. Both A SR and A CR are divided into a four b jet (4b) and three b jet (3b) region by requiring the b jet candidate with the lowest DEEPJET output to satisfy or fail the medium WP of the discriminant, respectively. There are between 5.5 and 11 times more events in the 3b region than in the 4b region, depending on the topological category and data-taking year considered. The overall efficiency for both ggF and VBF signal events to be selected in the A 4b SR region ranges from 0.3 to 3% depending on the couplings considered and is minimal for the SM VBF production and for the ggF production with κ λ ≈ 5 due to the interference effects in the HH production that result in low momenta of the Higgs bosons. The signal acceptance is mostly limited by the trigger acceptance and the jet b tagging efficiency.
Background events in the A Differences in the distributions of several variables between the 3b and the 4b regions are addressed with the BDT-based reweighting method described in Ref. [60], which uses a dedicated metric to identify the phase space regions with the largest differences in the distributions and compute an event weight to correct for them. This method accurately models multiple variables and their correlations, while minimizing issues related to the statistical uncertainties arising from the limited number of events in the two regions. The BDT is trained in the A 4b CR and A 3b CR regions, and applied to events in A 3b SR to model A 4b SR .
Trainings of this BDT are performed separately for each ggF and VBF category. All trainings use as inputs the following ten variables: p T of the four b jets, m HH , the invariant masses and p T of the H 1 and H 2 candidates, and the absolute value of their pseudorapidity difference (|∆η(H 1 , H 2 )|). In the ggF category, ten additional variables are used: the magnitude of the scalar (∑ p T ) and vector (p T (HH)) sums of the p T of the four b jets, the angular ∆R separations between the two jets that constitute H 1 and H 2 (∆R H 1 (bb), ∆R H 2 (bb)), the minimal ∆R (∆R min ) and the maximal |∆η| (|∆η| max ) between all the possible b jet pairs, the absolute value of the angle with respect to the beamline of one Higgs boson in the four-jet reference frame (|cos θ * |) and of one jet of H 1 in the H 1 candidate reference frame (|cos θ H 1 b |), the sum of the resolution estimators of the three b-tagged jets with the best DEEPJET value (∑ R e ), and the number of these three jets that satisfy the tight DEEPJET WP (N T b ). In the VBF category, four additional variables are used: the absolute value of the φ separation between the two Higgs bosons, the VBF jets invariant mass and absolute value of the η separation, and the output of the production mode BDT discriminant. These variables are chosen as those that best represent the kinematic properties of the events in the 3b and the 4b regions and that provide separation between signal and background and are used in subsequent steps of the analysis.
The training parameters are optimized with a two-step procedure. First, a Kolmogorov-Smirnov distance test is used to ensure that the distributions of the BDT input variables in the target A 4b CR region are compatible with the ones in the reweighted A 3b CR region. Once that is verified, a BDT is trained to separate the A 4b CR and reweighted A 3b CR data, thus testing also the correlations of variables. All the training configurations that are chosen are required to have an area of 0.5 under the receiver operating curve of the discriminant, corresponding to no separation.
The procedure is validated by applying it to a signal-depleted region, defined by shifting the signal and control regions according to the definition of χ, using as values of the center position c 1 = 179 GeV and c 2 = 172 GeV. The center of this validation region is chosen to be along the m H 1 = 1.04 m H 2 line used in the reconstruction of the H candidates to provide an accurate proxy of the analysis region. In analogy to the analysis regions, signal and control validation regions, V SR and V CR , are defined as χ < 25 GeV and 25 ≤ χ < 50 GeV, respectively. After training and applying the reweighting BDT in these regions and computing the normalization transfer factors, the data in V 4b SR were found to be compatible within uncertainties with the predicted background, validating the modeling method. The agreement is quantified with a goodness-of-fit test based on a saturated model [61] performed on the observables used in the analysis. For a fit under the background-only hypothesis, a p-value of 53% is observed, ranging between 12 and 83% for the individual categories.
The impact on the estimated background from the presence of signal events in the A 3b SR region due to jets failing the b tagging requirement is estimated by generating pseudo-data in A 4b SR according to the modeled background plus simulated HH signal, and fitting them under a different background hypothesis that includes the contribution from signal events in A 3b SR weighted as done for background events. This study is repeated for signal yields up to five times larger than the expected sensitivity of this search, and in all cases a signal yield compatible with the true one is observed. We conclude that signal events in A 3b SR do not have any significant impact on the background model and on the results. For the background model, systematic uncertainties are considered for the limited number of events in the A 3b SR . These uncertainties are uncorrelated across the individual bins of the background templates used for the statistical analysis, and correspond to the propagation of a bin-by-bin Poisson uncertainty from the A 3b SR to the A 4b SR region. The uncertainty in the estimation of the transfer factor from the 3b to the 4b region is computed from the statistical uncertainty in A 3b CR and A 4b CR and is 1-2% for the ggF categories, 2-3% for the SM-like VBF category and 18-32% for the anomalous κ 2V -like VBF category. An uncertainty is also considered for the limited number of events in the validation region, in some cases lower than the number of events in the analysis region. It is large (30-33%) for the anomalous κ 2V -like VBF category while it is about 2-3% and below 1% for the other VBF category and the ggF categories, respectively, and represents the inherent limitation on the capability to validate the performance of the background model. For analysis categories where the agreement between the observed and predicted background yields in the validation region differs by more than one standard deviation, an additional uncertainty is included and ranges between 1.5 and 4.7%, depending on the category and year. Finally, the uncertainty in the performance of the reweighting method to interpolate the kinematics from A CR into A SR is estimated by performing alternative trainings in two regions of A CR . The two regions are defined by requiring the product of m ⊥ and m to be either positive or negative, where m ⊥ is defined as the projection of the point in the (m H 1 , m H 2 ) plane onto the axis perpendicular to the one corresponding to m previously defined. The two regions correspond to four quadrants in the (m H 1 , m H 2 ) plane, and allow for tests of the capability of the reweighting method to model A SR , starting from events with kinematic properties that are either similar (m ⊥ m < 0) or that are harder or softer (m ⊥ m > 0) compared to A SR , thus testing the capability of the model to interpolate across different learning domains. The alternative background templates obtained from trainings in these regions represent the uncertainty on the shape of the predicted background distribution. All the uncertainties are independent between the 2016 and 2017-2018 background models. The dominant uncertainties in this search are those associated to the background modeling, and in particular the bin-by-bin and the normalization uncertainties due to the limited number of events in A The effects of the imperfect modeling of the detector response and the inaccurate simulation of signal processes are accounted for as systematic uncertainties. The most important sources of systematic uncertainty are the total integrated luminosity, the jet energy scale and resolution, the efficiency of the trigger and of the b-tagging requirements, the modeling of the pileup distribution, the HH → bbbb branching fraction, and the parameters used for the generators. A specific uncertainty on the parton shower is considered for the VBF production mode. Uncertainties on the theoretically determined HH cross section are considered only when quoting a limit on the HH signal strength (µ), defined as the ratio of the value of the cross section limit relative to the theoretical cross section expectation in the SM (σ HH /σ SM HH ). These uncertainties are negligible in comparison to the background uncertainties. See Appendix A for more details.
A multivariate BDT discriminant is trained with the XGBOOST software [62] in the two ggF subcategories to separate the signal from the weighted A 3b SR background events. The discriminant uses as inputs For each subcategory, a separate training is performed to separate the SM ggF HH signal from the weighted A 3b SR region data. Since the same A 3b SR data are also used to model the background, this data set is divided in two equal-size subsamples. Two trainings are performed on each half and applied to the other half, and the two partial background templates are added together. In this way, the full data set can be used for the modeling, while the BDT discriminant is not evaluated on events used for its training. In the VBF SM-like category, m HH is used as the discriminating variable, while in the anomalous κ 2V -like category, a counting experiment is performed because of the small number of expected background events. The distributions of these variables are shown in Fig. 1. For the VBF anomalous κ 2V -like category in the 2016 (2017-2018) data set, 4 (13) events are observed for a total of 4.0 ± 1.3 (15.0 ± 3.4) background and 1.5 (3.5) VBF κ 2V = 2 signal events expected.
A binned maximum likelihood fit is simultaneously performed in all analysis categories, where the systematic uncertainties previously discussed are introduced as nuisance parameters. No deviation from a background-only hypothesis is observed. Results are used to set 95% CL upper limits on the HH production cross section using the modified frequentist CL s criterion [63,64] with the profile likelihood ratio modified for upper limits [65] as the test statistics, and making use of the asymptotic approximation [66].  the 95% CL, while the value of κ 2V is observed (expected) to be in the range −0.1 < κ 2V < 2.2 (−0.4 < κ 2V < 2.5) at the 95% CL. The total HH production cross section, defined as the sum of the ggF and VBF production modes, is observed (expected) to be smaller than 120 (238) fb, corresponding to 3.9 (7.8) times the SM prediction, when uncertainties on the theoretical production cross section are included. The HH VBF production cross section is observed (expected) to be smaller than 226 (412) times the SM prediction. The deficit in the observed number of events localized around the BDT discriminant values of 0.85-0.9 in the ggF high-mass category in the 2017-2018 data set, which provides the largest sensitivity to the signal, results in the observed limit to be below the expected one. Studies of the background model for the individual BDT input variables and the absence of deficit in the 2016 data in the same high-sensitivity region, and in the other categories suggest that this under fluctuation is of statistical nature.  [216,846]) for the ggF (VBF) production modes, and the observed limit is thus compatible with the expectation within about two standard deviations. The sensitivity is mostly limited by the number of events in the signal and control regions of the analysis. Tabulated results are available in the HEPData record of this analysis [67].
In summary, a search for the production of Higgs boson pairs via gluon and vector boson fusion in the four b quark decay channel has been presented. The data are found to be statistically compatible with the background-only hypothesis, and an observed (expected) upper limit at the 95% confidence level is set to 3.9 (7.8) times the SM prediction for the combined ggF and VBF HH cross section. The value of the Higgs boson self-coupling, normalized to the SM expectation, is observed (expected) to be in the range −2.3 < κ λ < 9.4 (−5.0 < κ λ < 12.0), and the value of the coupling of Higgs boson pairs to vector boson pairs, normalized to the   Figure 2: Observed and expected 95% CL upper limits on the σ ggF+VBF HH cross section as a function of κ λ (left), and on the σ VBF HH cross section as a function of κ 2V (right). The green (yellow) band indicates the regions containing 68% (95%) of the limit values expected under the background-only hypothesis. The red lines denote the theoretical cross section expectation assuming that other couplings are set to the SM prediction. For the cross section limit as a function of κ 2V , the ggF HH production is assumed to correspond to the SM prediction.
SM expectation, to be in the range −0.1 < κ 2V < 2.2 (−0.4 < κ 2V < 2.5). These are the most stringent observed constraints to date on the HH production cross sections and on the κ 2V coupling.
QCD multijet processes are generated with MADGRAPH5 aMC@NLO at LO with MLM merging [70]. The production of pairs of top quarks is simulated at NLO with POWHEG, and normalized to the theoretical cross section at NNLO precision in QCD, including resummation of NNLL soft gluon terms [71]. The production of pairs of Z bosons decaying to hadrons with up to one jet emitted at the matrix element level is simulated with MADGRAPH5 aMC@NLO at NLO with FxFx merging [72] and normalized to the cross section at the same precision [73]. The total 2016-2018 integrated luminosity has an uncertainty of 1.6%, partially correlated for the three years [78][79][80]. The uncertainty on the jet energy scale and resolution is accounted for by modifying the corresponding jet properties within the uncertainties associated to their independent measurement, and affects the shape of the signal distributions and the associated yields by up to 5%, partially correlated for the three years. The trigger efficiency in data is measured in an orthogonal data set enriched in multijet tt events and compared to the same efficiency obtained from the simulation to derive a correction for the simulated signal events that depends on H T and on the p T and DEEPJET discriminant of the selected jets. The uncertainties of these measurements are considered as shape and normalization systematic uncertainties, uncorrelated for the three years, and range between 5 and 40% depending on the coupling and the category, with the largest values for the low-mass ggF category in regions of the phase space with little separation from the background. The impact of this uncertainty on the sensitivity is of about 1%. The b-tagging efficiency and the distribution of the DEEPJET discriminant are measured in tt and QCD multijet events, and independent uncertainties in this measurement for each year are propagated to the simulated events and range between 5 and 12%. For data recorded in 2016 and 2017, a gradual shift in the timing of the inputs of the ECAL Level-1 trigger in the forward endcap region (|η| > 2.4) led to an inefficiency whose impact is evaluated using an unbiased data sample, and an uncertainty of 1-2% on the corresponding correction is applied. The uncertainty associated with the correction of the pileup distribution is considered in the analysis by varying the assumed interaction cross section by 4.6%. Uncertainties arising from the theoretical modeling of the signal and related to the factorization and renormalization scales (2-8%), PDF (1-12%), and parton shower parameters (6-13%) used for the sample generation are considered. For the VBF production mode, two alternative models of signal generation that enable and disable a dipole-ordered parton shower in PYTHIA [81] are compared and their difference in the predicted number of signal events ranges from 2 to 13% and is considered to be a systematic uncertainty. Uncertainties in the HH → bbbb branching fraction and, when quoting a limit relative to the standard model prediction, theoretical uncertainties on the total HH cross section prediction are also considered.  For the former, the output of the BDT discriminant is shown for the low-mass category on the left and for the high-mass category on the right. For the latter, the m HH distribution in the SM-like category is shown. The anomalous κ 2V -like category is not shown because no shape correction is applied for it since the overall number of observed events is used to perform a counting experiment. Data are represented by points with error bars, while the ggF (VBF) signal contribution is shown in blue (red) and not stacked. The background prediction without the shape correction is represented by the shaded blue histograms with the associated systematic uncertainties (gray dashed areas).  For the former, the output of the BDT discriminant is shown for the low-mass category on the left and for the high-mass category on the right. For the latter, the m HH distribution in the SM-like category is shown. The anomalous κ 2V -like category is not shown because no shape correction is applied for it since the overall number of observed events is used to perform a counting experiment. Data are represented by points with error bars, while the ggF (VBF) signal contribution is shown in blue (red) and not stacked. The background prediction without the shape correction is represented by the shaded blue histograms with the associated systematic uncertainties (gray dashed areas). The numbers from 1 to 12 denote the shape benchmarks defined in Ref. [68], corresponding to points in the five-dimensional parameter space with characteristic kinematic properties of the HH system. The leading order simulation is reweighted to match the m HH spectrum at the next-to-leading order as predicted by simulation described in Ref. [69]. The inner (green) band and the outer (yellow) band indicate the regions containing 68 and 95%, respectively, of the distribution of limits expected under the background-only hypothesis.   and observed (right) regions allowed at the 68 and 95% CL as a function of the κ 2V and κ V couplings, assuming the κ t and κ λ couplings to be fixed to 1. The dot denotes the best-fit hypothesis while the star marker represents the SM prediction. and observed (right) regions allowed at the 68 and 95% CL as a function of the κ 2V and κ λ couplings, assuming the κ t and κ V couplings to be fixed to 1. The dot denotes the best-fit hypothesis while the star marker represents the SM prediction.     Figure A.13: Signal strength of the SM ggF and VBF production modes. The solid (dashed) line denotes the contour corresponding to the 68% (95%) confidence levels. The circle represents the best-fit point while the star marker represents the SM expectation. The signal strength is treated as a free parameter that is not bound to be positive, and the best-fit point takes a negative value as a consequence of the deficit in the observed number of the events discussed in the Letter.