Search for resonant pair production of Higgs bosons decaying to two bottom quark-antiquark pairs in proton-proton collisions at 8 TeV

A model-independent search for a narrow resonance produced in proton-proton collisions at sqrt(s) = 8 TeV and decaying to a pair of 125 GeV Higgs bosons that in turn each decays into bottom quark-antiquark pairs is performed by the CMS experiment at the LHC. The analyzed data correspond to an integrated luminosity of 17.9 inverse femtobarns. No evidence for a signal is observed. Upper limits at a 95% confidence level on the production cross section for such a resonance, in the mass range from 270 to 1100 GeV, are reported. Using these results, a radion with decay constant of 1 TeV and mass from 300 to 1100 GeV, and a Kaluza-Klein graviton with mass from 380 to 830 GeV are excluded at a 95% confidence level.


Introduction
Following the discovery of a Higgs boson (H) at the CERN LHC [1][2][3], with mass around 125 GeV and properties so far consistent with the standard model (SM) of particle physics, it has become important to search for new resonances that decay into pairs of such Higgs bosons.While non-resonant pair production of the Higgs boson is allowed in the SM, the theoretical production cross section is approximately 10 fb [4] and well beyond the sensitivity of currently acquired data.However, several well-motivated hypotheses of physics beyond the standard model posit narrow-width resonances that decay into pairs of Higgs bosons, and could be produced with large enough cross sections to be probed with existing data.The radion [5] and Kaluza-Klein (KK) gravitons in the Randall-Sundrum (RS1) [6] models of warped extra dimensions are examples of such resonances [7].
This letter reports the results of a model-independent search for the resonant pair production of Higgs bosons.The search for the narrow width resonance, denoted by X, is performed in the 270-1100 GeV mass range.Data from proton-proton collisions at the LHC and recorded by the CMS experiment corresponding to an integrated luminosity of 17.9 ± 0.5 fb −1 at √ s = 8 TeV is used.We perform this search for the case where both Higgs bosons decay into bottom quarkantiquark pairs (bb) [8].The main challenge of this search is to distinguish the signal of four bottom quarks in the final state that hadronize into jets (b jets) from the copious multijet background described by quantum chromodynamics (QCD) in pp collisions.We address this challenge by suitable event selection criteria that include dedicated b-jet identification techniques and a model of the multijet background that is validated in data control regions.Our results may be compared with a search performed by the ATLAS experiment [9] that also probes the physics of resonant Higgs boson pair production, albeit in the channel where one Higgs boson decays to bottom quarks and the other decays to photons.

Detector and Event Reconstruction
A detailed description of the CMS detector, together with a description of the coordinate system used and the relevant kinematic variables, can be found in Ref. [10].The central feature of the CMS apparatus is a superconducting solenoid of 6 m internal diameter that generates an axial magnetic field of 3.8 T. Within the solenoid volume are a silicon pixel and strip tracker, a lead tungstate crystal electromagnetic calorimeter, and a brass and scintillator hadron calorimeter.Muons are detected and their properties measured in gas-ionization detectors embedded in a steel flux-return yoke outside the solenoid.Jets are reconstructed using the anti-k T clustering algorithm [11,12] with a distance parameter of 0.5 applied on the collection of particle candidates reconstructed by the particle-flow (PF) algorithm [13,14].The PF algorithm reconstructs and identifies each individual particle with a combination of information from the various elements of the CMS detector.To mitigate the effect of additional particles that do not originate from the hard interaction in jet reconstruction, we subtract charged hadrons that do not arise from the primary vertex associated with it from the collection of clustered particles.Further, an average neutral energy density from particles not arising from the primary vertex is evaluated and subtracted from the jets [15].Energy corrections for the jets are determined as functions of the jet transverse momentum p T and pseudorapidity |η|.Jet identification criteria [16] to reject detector noise misidentified as jets, and jets not originating from the hard interaction are also applied.
In order to identify (tag) b jets, we rely on the fact that bottom quarks hadronize into b hadrons which have decay lengths of the order of cτ = 450 µm.Thus, their decay products originate from secondary vertices made of tracks that have impact parameters with respect to the primary vertex of a similar scale.The pixel tracker provides an impact parameter resolution of about 15 µm for charged tracks with |η| < 2.4.To maximize the b-tagging performance of the detector, we combine the output discriminants of several b-tagging algorithms described in Ref. [17] with a trained artificial neural network.This we call the combined multivariate (CMVA) algorithm.In particular, we combine the outputs of the combined secondary vertex (CSV) tagger that uses secondary vertices identified by the inclusive vertex finder (IVF) algorithm [18], the jet probability (JP) tagger, and the two soft lepton taggers.We choose a working point for the CMVA algorithm that maximizes the sensitivity of this search.This working point yields a 75% efficiency for tagging jets originating from b hadrons and a mistagging rate of 3% for light-flavor jets with p T > 40 GeV and |η| < 2.4.
The first level of the trigger, consisting of customized processors, collects data for this analysis using information from the calorimeters and requires two jets to exceed p T thresholds of 56 or 64 GeV, depending on luminosity conditions.The second level of the trigger, consisting of software algorithms executed on a farm of commercial processors, uses information from the entire detector to reconstruct PF jets, and requires four PF jets with |η| < 2.4 and p T > 30 GeV, of which two jets must have p T > 80 GeV.Further, to record signal events and reject background QCD multijet events, two jets are required to be tagged by the CSV b-tagging algorithm implemented at the trigger.

Simulated Samples
To model the production of a generic narrow-width spin-0 resonance, we use a Monte Carlo simulation of the RS1 radion produced through gluon fusion.The angular distributions of a spin-2 resonance are distinct from that of a spin-0 resonance, and result in different kinematic distributions.Therefore, we evaluate the signal efficiencies for a narrow-width spin-2 resonance from a separate simulation of the first excitation of the KK graviton produced through gluon fusion in the same extra dimension scenario as the radion.The resonance is forced to decay to a pair of Higgs bosons with parameters reported in Ref. [19], where both Higgs bosons decay to bb.Samples of these signal events, as well as background events from diboson, W+jets, Z+jets and top-quark pair production (tt) processes, are generated using the MADGRAPH 5.1 [20] program interfaced with PYTHIA 6.4 [21] for parton showering and hadronization.QCD multijet event samples are simulated with the PYTHIA 6.4 program.A sample of events where the Higgs boson is produced in association with a Zboson is simulated using the POWHEG event generator [22][23][24] interfaced with the HERWIG++ [25] program for showering and hadronization.We set the PYTHIA 6.4 parameters for the underlying event to the Z2* tune [26].The response of the CMS detector is modeled using GEANT4 [27].
On average, 21 pp interactions occurred per bunch crossing in the data used in this analysis.Additional simulated pp interactions overlapping with the event of interest were added to the simulated samples to reproduce the distribution of the number of primary vertices per event reconstructed in data.

Event Selection
The trigger-level jet p T thresholds confine our search for a narrow-width X → H(bb)H(bb) resonance to masses above 270 GeV.Beyond m X ≈ 800 GeV, the selection efficiency is increasingly limited by the merging of jets from the same Higgs boson, and we curtail this search at 1100 GeV.The kinematic distributions of the decay products vary substantially over this mass range.Therefore, to optimize the search sensitivity, we use different event selection criteria in two main kinematic regions: the low-mass region (LMR) for mass hypotheses from 270 to 450 GeV, and the high-mass region (HMR) for masses from 450 to 1100 GeV.
Event selection begins with the identification of events containing at least four jets in the central region of the detector (|η| < 2.4) that are b-tagged and have p T > 40 GeV.For the LMR, we combine these b jets into pairs to create HH candidates such that |m H − 125 GeV| < 35 GeV for each candidate Higgs boson.The mass resolution on the Higgs boson in the LMR is found to be approximately 9 GeV.Selected HH candidates are required to have at least two jets with p T > 90 GeV.Signal events have large Lorentz factors for the Higgs boson candidates in the HMR, and therefore HH candidates for this case are constructed from four jets such that the ∆R = √ ∆η 2 + ∆φ 2 between the jets associated with an H candidate remain within 1.5, where ∆η and ∆φ are the differences in the pseudorapidities and azimuthal angles of the two jets.For mass hypotheses from 740 to 1100 GeV, we impose an additional requirement of p T > 300 GeV on one of the Higgs boson candidates to better discriminate signal events from background.In case of multiple HH candidates in an event, the combination with the smallest |m Having identified the two Higgs boson candidates in each event, we plot their masses, m H 1 and m H 2 , on a two-dimensional histogram as shown in Fig. 1.H 1 and H 2 are chosen at random from the two reconstructed H candidates.As the final selection criterion, we require events to fall within the signal region (SR) defined as The efficiencies of these selection criteria for spin-0 and spin-2 resonances at representative masses are shown in Table 1.The major loss in efficiency for all mass hypotheses comes from the b-tagging requirement for 4 jets.For the 300 GeV mass hypothesis, this is compounded by  the trigger inefficiency.The distribution of the aforementioned ∆R between jets from a single Higgs boson is narrower for the spin-2 resonance, and thus requiring ∆R < 1.5 results in a higher efficiency for it.

Signal Modeling
For signal events, the aforementioned event selection criteria are expected to produce a sharp peak in the m X distribution over a relatively featureless background from events arising from SM processes.To search for signal events at various mass hypotheses, we fit the m X distribution in data events in the SR to a parametric model for the signal peak on top of parametric models appropriate for components of the SM background.This procedure is performed for the LMR and the HMR separately.
To improve the mass resolution of the signal X → HH resonance, we perform a fit that constrains the invariant masses of the Higgs boson candidates.In the fit, the momenta of the reconstructed b jets are allowed to float within their expected resolutions.Since the uncertainty in the reconstructed mass of the Higgs boson candidate due to the measurement of jet direction is smaller than that due to the measurement of jet energy, this constraint mainly affects the latter.This fit improves the invariant mass resolution of the reconstructed signal resonance by 20-40%, depending on the mass hypothesis.Extensive tests in background-dominated control regions in data show that no artificial structures are introduced in the background mass distributions by this procedure.
We build the parametric model for each signal mass hypothesis by fitting the shape of the m X distribution of simulated events that are accepted by the selection criteria and corrected for differences between data and simulation.A sum of two Gaussian functions, requiring five parameters, is used for the LMR fit to account for tails in the distribution from incorrect combinations of jets.In the HMR, we fit a function with a Gaussian core smoothly extended on both sides to exponential tails, such that the function is continuous both in its value and its first derivative.This requires two parameters for the mean and width of the Gaussian function, and two other parameters for the exponential tails on both sides.An example of a parametric model for the HMR signal obtained through this procedure is shown in Fig. 2. While the model is constructed for the mass hypothesis of 700 GeV, its Gaussian core peaks at 714 GeV and has a width of 21 GeV.This mass shift is found to be linear in m X and occurs due to the aforementioned constraint of jet momenta to m H .

Background Modeling
While the composition of background events in the SR is expected to be dominated by QCD multijet processes, we find through simulation that tt production contributes approximately 22% and 27% in the LMR and HMR, respectively.We also find that Z+jets, ZZ, and ZH processes contribute less than 1% of the background and therefore neglect them in this analysis.The m X distribution of these tt events is found to be somewhat different in shape from that of QCD multijet events, and therefore we treat it as a distinct component of the background and model it with a parametric form.We obtain this parametric form by fitting the shape of the m X distribution of simulated tt events accepted by the event selection criteria to a function with a Gaussian core smoothly extended to an exponential tail on the high side.This function, henceforth referred to as GaussExp, is continuous in its value and its first derivative.It has two parameters for the mean and width of the Gaussian function and one parameter for the decay constant of the exponential tail.This model is normalized to a tt cross section of 234 pb [28], and is allowed to float with a systematic uncertainty of 15% in the final fit to account for theoretical and measurement uncertainties in our kinematically boosted region.
We use the GaussExp parametric model to fit the m X distribution of the QCD multijet component of the background in the SR.With the SR kept blinded, we motivate and validate this choice of parametric model by the fact that it fits well the shape of the m X QCD multijet background distributions in several different regions of the (m H 1 , m H 2 ) plane depicted in Fig. 1 and described below.We do not aim to predict the parameters of the model in the SR from the other regions.These fits are performed for the LMR between 260 and 650 GeV, for the HMR between 400 and 900 GeV, and for the HMR with H p T > 300 GeV between 600 and 1200 GeV.In each case the tt contribution, as expected from simulation, is subtracted.
We define a sideband region (SB) to the SR as 17.5 GeV < √ ∆m 2 H 1 + ∆m 2 H 2 < 35 GeV and ∆m H 1 ∆m H 2 < 0. For events in this region, the m X distribution is expected to be kinematically similar to that for events in the SR-as in each of the sidebands one of the reconstructed Higgs boson masses is slightly higher in value than for events in the SR while the other is slightly lower.As an example, Fig. 3 top left shows the fit performed for events in the SB passing the HMR selection.Another set of events that pass the kinematic requirements of the event selection criteria in the SR region of the (m H 1 , m H 2 ) plane but required to have one of the four jets Here n is the number of degrees of freedom in each fit.The pull, for a given bin, is defined as the number of data events minus the value of the fit model, divided by the uncertainty on the number of data events.not be b-tagged is selected to further test the applicability of the GaussExp model in describing the m X distribution of the QCD multijet background in a different but kinematically similar region.This is called the data Control Region (CR), and the fit for these events, that would have otherwise passed the HMR selection, is also shown in Fig. 3 on the top-right.In both cases, the goodness of the fit, characterized by the χ 2 per degree of freedom, is found to be reasonable.
These two cases already lend significant confidence to the choice of the GaussExp parametric model for the SR.However, we carry out further checks in neighboring validation regions (VR) with a corresponding sideband (VS) that are defined similarly to the SR and SB regions but with m H 1 = m H 2 centered at different values.The good fits for the m X distributions in these regions not only demonstrate the applicability of the GaussExp model to describe these kinematically distinct QCD multijet events, but also that events in the VR are in fact kinematically similar to those in the VS.As examples, Fig. 3 bottom-left and bottom-right plots show the results of these fits for the LMR selection for the VS and the VR, respectively, both centered at m H 1 = m H 2 = 90 GeV.We obtain similar results for the VR centered at m H 1 = m H 2 = 107.5,142.5 and 160 GeV.
While the GaussExp function fits well the m X distribution from QCD multijet events in all these distinct regions and therefore can be expected to be a good approximation of the parametric form of the true parent distribution for events in the SR, other similar parametric models could also be chosen instead of it.Therefore, a systematic uncertainty associated with the choice of this parametric model is evaluated by assuming a 7th order Bernstein polynomial, which also fits the m X distribution well in the SB, to be the true distribution.Pseudo-datasets are generated from this polynomial function and fitted with the GaussExp function as well as other polynomial functions to compute biases in the reconstructed signal strength.This procedure is performed for each mass hypothesis.These biases are found to be of the order of 100 fb for the LMR and 10 fb for the HMR.We account for this bias as a signal-shaped systematic uncertainty on the background model with normalization centered at zero and a Gaussian uncertainty with standard deviation equal to the bias.

Systematic Uncertainties
Sources of systematic uncertainties that affect the selection efficiencies for signal and tt events are listed in Table 2 and described below.We vary the jet energy scale [29] by its uncertainty as a function of jet p T and η and find that this affects signal efficiencies by a relative factor of up to 0.2% and tt efficiencies by up to 0.8%.We evaluate the effect of uncertainty in the jet energy resolution by varying the jet energies according to the measured uncertainty.This is found to affect signal efficiencies by 2-7%, and tt efficiencies by 2-3%.These uncertainties affect not only the normalizations but also the parameters of the signal and tt models, and are taken into account as nuisance parameters in the final fit.
The trigger efficiencies for signal and tt events are evaluated approximately by passing generated events through a trigger simulation.We then correct these efficiencies for differences between simulation and data by computing the difference in a tt-enriched control region obtained using a trigger that requires at least one muon with p T > 24 GeV.Further, event selection criteria requiring at least one muon with p T > 40 GeV and at least four jets in the central region of the detector with p T > 40 GeV are applied.The data-to-simulation correction factor is characterized by the p T and CMVA discriminants of the relevant jets.Uncertainties in this factor impact signal efficiencies by 6-18%, and tt efficiencies by 7-9%.
The b-tagging efficiencies of the CMVA algorithm for signal and tt events are also evaluated approximately through simulation and then corrected by a data-to-simulation comparison.The comparison is performed in the same tt-enriched control region as the calculation of the trigger efficiency.The correction factor for the b-tagging efficiencies is consistent with unity.The uncertainty on this factor is evaluated to be 13%.
Additionally, the yield of signal events for a given production cross section and tt events are both affected by a 2.6% uncertainty in the measurement of the integrated luminosity [30].

Results
The m X distribution that we observe in data within the SR for both the LMR and HMR, along with a binned maximum-likelihood fit with the aforementioned parametric background models are shown in Fig. 4. We compute the observed and expected upper limits on the cross section for pp → X → H(bb)H(bb) at a 95% confidence level (CL) using the modified frequentist CL S method [31,32] by fitting the data with the parametric signal, tt, and QCD multijet models.This is done separately in the disjoint ranges of m X for the LMR and the HMR described in Section 4, and the limits are presented together in Fig. 5.These limits are shown for the spin-0 resonance on the left, and the spin-2 resonance on the right.The green (dark) and yellow (light) bands represent the 1σ and 2σ confidence intervals around the expected limits.The observed upper limits lie within 2σ of the expected upper limits, and thus we conclude that there is no significant deviation from the background-only hypothesis.
The theoretical cross section for the production via gluon fusion of a radion that decays to a pair of Higgs bosons [33] that each in turn decays to a bb pair with a branching fraction of 58% [34], as calculated by MADGRAPH 5.1, is overlaid on the limit for the spin-0 resonance in the plot on the left.In this calculation, the correction factor used to account for next-toleading-order amplitudes of the heavy-quark loop is identical to that used for Higgs boson production through gluon fusion [35].The warped extra dimension scenario for this radion has the product of the curvature, k, and half the circumference of the extra dimension, L, set to 35, a radion decay constant of Λ R = 1 TeV, and no radion-Higgs boson mixing.The theoretical cross section for the radion has an uncertainty of approximately 15% that is not used to compute the experimental limits on spin-0 resonance production shown in Fig. 5. Masses for the radion between 300 and 1100 GeV are excluded at a 95% CL.A similarly calculated theoretical cross section for the KK graviton as the resonance X, in the same warped extra dimension scenario, is overlaid on the limit for the spin-2 resonance in the plot on the right.Masses for such a graviton are excluded at a 95% CL between 380 and 830 GeV. between 400 and 900 GeV of the HMR (top-right), and between 600 and 1200 GeV in the HMR where one of the Higgs boson candidates has p T > 300 GeV (bottom).All distributions are fitted to the background-only hypothesis for illustration, showing the relative contributions of the QCD-multijet (dashed-dotted red) and tt (dashed green) processes.Also for illustration, we overlay the signal models of the spin-0 resonance (dotted blue) corresponding to mass hypotheses and production cross sections of 350 GeV and 653 fb for the LMR, 700 GeV and 17.6 fb for the HMR, and 900 GeV and 8.1 fb for the HMR with H p T > 300 GeV on their respective plots.These cross sections correspond to the observed upper limits, which are computed for signal mass hypotheses from 270 to 450 GeV in the LMR, from 450 to 730 GeV in the HMR, and from 730 to 1100 GeV in the HMR with H p T > 300 GeV.The observed and expected upper limits on the cross section for pp → X → H(bb)H(bb) at a 95% confidence level, where the resonance X has spin-0 (left) and spin-2 (right).The theoretical cross section for the RS1 radion, with Λ R =1 TeV, kL = 35, and no radion-Higgs boson mixing, decaying to four b jets via Higgs bosons is overlaid on the left plot.The theoretical cross section for the first excitation of the KK-graviton for the same parameters is overlaid on the right plot.

Summary
We have presented a model-independent search by the CMS experiment at the LHC for a narrow resonance produced in proton-proton collisions at √ s = 8 TeV and decaying to a pair of 125 GeV Higgs bosons that in turn each decays into bottom quark-antiquark pairs.The analyzed data correspond to an integrated luminosity of 17.9 fb −1 .No evidence for a signal is observed.Upper limits at a 95% CL on the production cross section for such spin-0 and spin-2 resonances, in the mass range from 270 to 1100 GeV, are reported.Using these results, a radion with decay constant of 1 TeV and mass from 300 to 1100 GeV, and a Kaluza-Klein graviton with mass from 380 to 830 GeV are excluded at a 95% confidence level.

Figure 1 :
Figure 1: Illustration of the SR, SB, VR and VS kinematic regions in the (m H 1 , m H 2 ) plane used to motivate and validate the parametric model for the QCD multijet background.The quantities m H 1 and m H 2 are the two reconstructed Higgs boson masses.The distribution in data events after b-tagging and kinematic selections is shown with the SR blinded.

Figure 2 :
Figure2: A maximum likelihood fit to the m X distribution of simulated signal events for the 700 GeV mass hypothesis.The distribution is fitted to a Gaussian core smoothly extended on both sides to exponential tails.Here n is the number of degrees of freedom in the fit.

Figure 3 :
Figure3: The m X distributions in data (after the tt background has been subtracted) in the SB of the HMR (top-left), the CR of the HMR (top-right), the VS of the LMR (bottom-left), and the VR of the LMR (bottom-right).The distributions are fitted to the GaussExp parametric model.The shaded regions correspond to ±1σ variations of this fit.Here n is the number of degrees of freedom in each fit.The pull, for a given bin, is defined as the number of data events minus the value of the fit model, divided by the uncertainty on the number of data events.

Figure 4 :
Figure4: The m X distribution in data in the SR between 260 and 650 GeV of the LMR (top-left), between 400 and 900 GeV of the HMR (top-right), and between 600 and 1200 GeV in the HMR where one of the Higgs boson candidates has p T > 300 GeV (bottom).All distributions are fitted to the background-only hypothesis for illustration, showing the relative contributions of the QCD-multijet (dashed-dotted red) and tt (dashed green) processes.Also for illustration, we overlay the signal models of the spin-0 resonance (dotted blue) corresponding to mass hypotheses and production cross sections of 350 GeV and 653 fb for the LMR, 700 GeV and 17.6 fb for the HMR, and 900 GeV and 8.1 fb for the HMR with H p T > 300 GeV on their respective plots.These cross sections correspond to the observed upper limits, which are computed for signal mass hypotheses from 270 to 450 GeV in the LMR, from 450 to 730 GeV in the HMR, and from 730 to 1100 GeV in the HMR with H p T > 300 GeV.

Figure 5 :
Figure5: The observed and expected upper limits on the cross section for pp → X → H(bb)H(bb) at a 95% confidence level, where the resonance X has spin-0 (left) and spin-2 (right).The theoretical cross section for the RS1 radion, with Λ R =1 TeV, kL = 35, and no radion-Higgs boson mixing, decaying to four b jets via Higgs bosons is overlaid on the left plot.The theoretical cross section for the first excitation of the KK-graviton for the same parameters is overlaid on the right plot.

Table 1 :
Efficiencies of the event selection criteria for generic spin-0 and spin-2 resonances decaying to a pair of Higgs bosons in the four b jet final state at representative masses.

Table 2 :
Relative systematic uncertainties on the selection efficiencies for signal and tt events in the LMR and the HMR.