Combination of CMS searches for heavy resonances decaying to pairs of bosons or leptons

A statistical combination of searches for heavy resonances decaying to pairs of bosons or leptons is presented. The data correspond to an integrated luminosity of 35.9 fb$^{-1}$ collected during 2016 by the CMS experiment at the LHC in proton-proton collisions at a center-of-mass energy of 13 TeV. The data are found to be consistent with expectations from the standard model background. Exclusion limits are set in the context of models of spin-1 heavy vector triplets and of spin-2 bulk gravitons. For mass-degenerate W' and Z' resonances that predominantly couple to the standard model gauge bosons, the mass exclusion at 95% confidence level of heavy vector bosons is extended to 4.5 TeV as compared to 3.8 TeV determined from the best individual channel. This excluded mass increases to 5.0 TeV if the resonances couple predominantly to fermions.


Introduction
Over the past half century, successive searches for heavy resonances in two-body decays have often led to discoveries of new states. At the CERN LHC, searches for W and Z resonances (collectively referred to as V ) that couple through the electroweak (EW) interaction to standard model (SM) particles have been performed by ATLAS and CMS in final states with two SM bosons [1-18], two leptons [19][20][21][22], and two light-flavor [23,24] or heavy-flavor quarks [25][26][27][28][29]. These new states would couple predominantly to either SM fermions, as in the case of minimal W and Z models [30][31][32], or SM bosons, as in strongly coupled composite Higgs and little Higgs models [33][34][35][36][37][38][39]. In addition, models based on warped extra dimensions provide a candidate for a massive resonance, such as the spin-2 first Kaluza-Klein excitation of the graviton (G) [40][41][42]. In bulk graviton models [43], the SM fermions and gauge bosons are located in the bulk five-dimensional spacetime, and the graviton has a sizable branching fraction to pairs of W, Z, and H bosons.
This Letter describes a statistical combination of CMS searches for heavy resonances that decay to pairs of bosons or leptons [1-10, 19, 20]. The event selection, the simulated samples, the background estimation, and the systematic uncertainties of the individual analyses are unchanged. The results are used to set exclusion limits on models that invoke heavy vector resonances and on a model with a bulk graviton.
The SM bosons produced in the decay of resonances X with masses m X that exceed 1 TeV are expected to have a large Lorentz boost [44]. The decay products of the SM bosons are therefore highly collimated, requiring dedicated techniques for their identification and reconstruction. In the case of hadronic decays, the pair of quarks produce a single large-cone jet with a two-prong structure. In addition, Higgs bosons can also be identified by tagging the b quarks originating from their decays. In models where heavy vector bosons couple predominantly to fermions, the combined contribution from the leptonic decays W → ν and Z → dominates, and these channels [19,20] are included in the combination, providing complementary sensitivity to the searches in diboson channels. The analyses considered in this Letter are based on a data sample of proton-proton (pp) collisions at a center-of-mass energy of 13 TeV collected by the CMS experiment in 2016, and corresponding to an integrated luminosity of 35.9 fb −1 . A similar combination performed on a comparable set of data has been recently published by ATLAS [45].

The CMS detector and event reconstruction
The CMS detector features a silicon pixel and strip tracker, a lead tungstate crystal electromagnetic calorimeter (ECAL), and a brass and scintillator hadron calorimeter (HCAL), each composed of a barrel and two endcap sections. These detectors reside within a superconducting solenoid, which provides a magnetic field of 3.8 T. Forward calorimeters extend the pseudorapidity coverage up to |η| < 5.2. Muons are measured in gas-ionization detectors embedded in the steel flux-return yoke outside the solenoid. A detailed description of the CMS detector, together with a definition of the coordinate system and the kinematic variables, can be found in Ref. [46].
The information from various elements of the CMS detector is used by a particle-flow (PF) algorithm [47] to identify stable particles reconstructed in the detector as electrons, muons, photons, and charged or neutral hadrons. The energy of electrons is determined from a combination of the electron momentum, as determined in the tracker, the energy of the corresponding ECAL cluster, and the energy sum of all bremsstrahlung photons spatially compatible with originat-ing from the electron track. The dielectron mass resolution for Z → ee decays ranges from 1.9% when both electrons are in the ECAL barrel, to 2.9% when both electrons are in the endcaps [48]. The momentum of muons is obtained from the curvature of the corresponding track. The p T resolution in the barrel is better than 7% for muons with p T up to 1 TeV [49]. The τ leptons that decay to hadrons and a neutrino (labeled τ h ) are reconstructed by combining one or three charged PF candidates with up to two neutral pion candidates [50]. The energy of charged hadrons is determined from a combination of their momenta measured in the tracker and the matching ECAL and HCAL energy deposits. Finally, the energy of neutral hadrons is obtained from the corresponding corrected ECAL and HCAL energies.
Jets are reconstructed from PF candidates clustered with the anti-k T algorithm [51], with a distance parameter 0.4 (AK4 jets) or 0.8 (AK8 jets), using the FASTJET 3.0 package [52,53]. The four-momenta of the AK4 and AK8 jets are obtained by clustering candidates passing the charged-hadron subtraction algorithm [54]. The contribution of neutral particles originating from additional pp interactions within the same or neighboring bunch crossings (pileup) is proportional to the jet area, and is estimated using the median area method implemented in FASTJET, and then subtracted from the jet energy. The jet energy resolution, after the application of corrections to the jet energy, is 4% at 1 TeV [55]. The missing transverse momentum vector p miss T is computed as the negative vector sum of the transverse momenta of all the PF candidates in an event, and its magnitude is denoted as p miss T [56]. The p miss T is corrected for adjustments to the energy scale of the reconstructed AK4 jets in the event.
While AK4 jets are used for single quarks, AK8 jets are adopted to reconstruct large momentum SM bosons that decay to quarks. Reconstructing the AK8 jet mass (m j ) and substructure relies on the pileup-per-particle identification algorithm [54,57]. The contributions from soft radiation and additional interactions are removed using the soft-drop algorithm [58,59], with parameters β = 0 and z cut = 0.1. Dedicated mass corrections, obtained from simulation and data in a region enriched with tt events with merged W(qq ) decays, are applied to the jet mass to reduce any residual jet p T dependence [6,60], and to match the jet mass scale and resolution observed in data. The measured soft-drop jet mass resolution is approximately 10% [60]. Exclusive m j intervals labeled m W , m Z , and m H , which range from 65 to 85, 85 to 105, and 105 to 135 GeV, respectively, are defined according to the SM boson masses and the jet mass resolution.
Hadronic decays of W and Z bosons are identified using the ratio between 2-subjettiness and 1-subjettiness [61], τ 21 = τ 2 /τ 1 . The variables m j and τ 21 are calibrated using a top quarkantiquark sample enriched in hadronically decaying W bosons [62]. The decay of a Higgs boson to a pair of b quarks is identified using one of two b tagging algorithms, depending on the background composition. The first consists of a dedicated b tagging discriminator, specifically designed to identify a pair of b quarks clustered in a single jet [63]. The second relies on splitting the AK8 jet into two subjets, then applying the combined secondary vertex algorithm [63] to the subjets. The latter is also applied to AK4 jets to identify isolated b quarks in the event.

Signal modeling
The response of the CMS detector to the production and decay of heavy resonances is evaluated through simulated events, which are reconstructed using the same algorithms as used in collision data. The spin-1 gauge bosons, W and Z , are simulated at leading order (LO) using the MADGRAPH5 aMC@NLO 2.4.2 matrix element generator [64], in the heavy vector triplet (HVT) framework [65,66], which introduces a triplet of heavy vector bosons with sim-ilar mass, of which one is neutral (Z ) and two are electrically charged (W ± ). The coupling strength of the heavy vector bosons to SM bosons and fermions is determined from the combinations g H = g V c H and g f = g 2 c F /g V , respectively. The parameter g V is the strength of the new interaction; c H characterizes the interaction between the HVT bosons, the Higgs boson, and longitudinally polarized SM vector bosons; c F represents the direct interaction between the V bosons and the SM fermions; g is the SM SU(2) L gauge coupling constant. The HVT framework is presented in two scenarios, henceforth referred to as model A and model B, depending on the couplings to the SM particles [66]. In the former, the coupling strengths to the SM bosons and fermions are comparable, and the new particles decay primarily to fermions. In the latter, the couplings to the SM fermions are small, and the branching fraction to the SM bosons is nearly 100%. Events are simulated with different m X hypotheses in the range of 800 to 4500 GeV, assuming a negligible resonance width compared to the experimental resolution (3-20%). The kinematic distributions of the signal do not depend on the choice of benchmark scenario, and the samples are reweighted according only to the cross sections and branching fractions.
The bulk graviton events are also simulated at LO using the same MADGRAPH5 aMC@NLO generator. The cross section and the width of the bulk graviton mainly depend on its mass and the ratio κ ≡ κ/M Pl , where κ is a curvature factor of the model and M Pl is the reduced Planck mass [43,67]. The graviton signals are generated assuming κ = 0.5, which guarantees that the width of the graviton is smaller than the experimental resolution.
The signal events are generated using the NNPDF 3.0 [68] parton distribution functions (PDFs), and are interfaced to PYTHIA 8.205 [69] for parton showering and hadronization, adopting the MLM matching scheme [70]. Pileup interactions are superimposed on the simulated processes, and their frequency distribution is weighted to match the number of interactions per bunch crossing observed in data. Generated events are processed through the CMS detector simulation, based on GEANT4 [71].

Search channels
The search strategies for the analyses contributing to the combination are summarized in this Section. More details are provided in the relevant publications [1-10, 19, 20].

Fully hadronic diboson channels
Diboson resonances have been searched for in several final states, depending on the decay modes of the bosons. The final states targeted were VV → qqqq [1], VH → qqbb [6], and HH → bbbb final states [9, 10]. Each boson is reconstructed as a two-prong, large-cone jet, so that for diboson resonances, the two jets would recoil against each other. The presence of a diboson resonance would be observed in the dijet invariant mass spectrum m jj . The W and Z bosons are identified via the m j and τ 21 variables, and b tagging is used to identify b quarks from Higgs bosons in addition to m j . Although the signal yield is large, because of the large fraction of Higgs boson decays to b quarks, these channels are subject to an overwhelming background from quantum chromodynamics (QCD) multijet production.
In the VV and VH analyses, the background is estimated directly from data, assuming that the invariant mass distribution of the background can be described by a smooth, parameterizable, monotonically decreasing function of m jj . The signal template, based on a Gaussian core, is fitted to the data simultaneously with the background function. The HH analyses also use an additional region in the fit, obtained by inverting the b tagging selection on the H candidates, which constrains the parameters of the background function. ). These final states represent an attractive alternative to all-jet final states, thanks to the large selection efficiencies and natural discrimination against multijet background stemming from the presence in the signal of energetic and isolated leptons or neutrinos.

Semi-leptonic diboson channels
The decay of a Z boson to neutrinos can be identified through its large p miss T , and the resonance mass can be inferred from the transverse mass between the p miss T and the jet originating from the hadronic decay of the other boson. For the W → ν decay, there is a single, isolated lepton associated with a moderate p miss T , and the vector boson can therefore be reconstructed by imposing a constraint from the W boson mass to recover the longitudinal momentum of the neutrino. In Z → decays, two opposite-sign, same-flavor leptons with a combined invariant mass compatible with the Z boson mass are used to determine accurately the Z boson four-momentum. A Higgs boson decaying to τ leptons is identified from dedicated τ h decay reconstruction and isolation techniques [8], and its mass is estimated through the measured momenta of the visible decay products of the τ leptons and the p miss T , with the SVFIT algorithm [72]. The bosons that decay to a pair of quarks are reconstructed as AK8 jets.
In the semi-leptonic analyses, the main V+jets background is estimated from a fit to data in the m j sidebands of the hadronic jet, and extrapolated to the signal region using a transfer function ("α function") obtained from simulation. The top quark pair production is estimated from simulation, but its normalization is rescaled to match the data in control regions obtained by requiring an additional b-tagged AK4 jet in the event [2, 7, 8].
In addition, the WV → νqq analysis introduces a novel signal extraction method based on a two-dimensional (2D) fit to data [3]. The backgrounds are separated into non-resonant and resonant categories depending on the presence or absence of W bosons and top quarks in the jet mass spectrum, and are fitted simultaneously in the space of m j and m WV , accounting for the correlation between the two variables.

Fully leptonic diboson channels
Searches for diboson resonances decaying to a pair of Z bosons have been performed in fully leptonic final states, with one boson undergoing the decay Z → and the other Z → νν [5]. The presence of the leptons and neutrinos defines a very clean final state with reduced backgrounds, but the small branching fraction makes this channel competitive only for small resonant mass values.

Decays to a pair of fermions
The decay of a heavy resonance to a pair of fermions can be sizable when the couplings to SM fermions are large. If the resonance is electrically charged, as in the case of W ± , the decay to a neutrino and an electron or muon yields a broad excess in the ν transverse mass spectrum [19]. If the new state is neutral, as in the case of Z , a narrow resonance would emerge from the dielectron or dimuon ( ) invariant mass spectra [20]. The analyses of these fermionic decays extend to masses above 5 TeV, and employ selection techniques optimized to identify and measure very energetic electrons and muons in the detector [19,20].

Event selection
The search regions of the analyses entering in the combination are statistically independent because of mutually exclusive selections on the number of leptons and their flavor, number of AK8 jets, and jet mass intervals. Analyses with hadronic final states reject events with isolated leptons or with large p miss T . Overlaps between channels that share the same lepton multiplicity are avoided by selecting different jet mass ranges. The W → ν search does not share any events with the W → VW and W → WH analyses because of the requirements on the angular separation ∆φ( , p miss T ) between the p miss T and the lepton direction. The Z → analysis includes events with a dilepton invariant mass m > 120 GeV, which is incompatible with the 70 < m < 110 GeV selection used in the diboson channels, in which the Z boson is onshell. The two searches for resonant HH bosons that decay to b quarks have common events explicitly removed [10].
In the WV → νqq channel [3] the background is estimated using a 2D fitting technique that scans the full jet mass range, and is therefore not independent of the WH → νbb channel. For this reason, in the W , Z , and V interpretations, where the two signals might be present simultaneously, the "α function" is used to estimate the background instead. This method considers only events in the m j regions of the W and Z bosons, thereby preventing double counting of events in the Higgs boson mass region. The results of the alternative background estimation method are consistent with those obtained in the 2D fit, but the method is about 10% less sensitive. The main selections that define the exclusivity of the analyses are summarized in Table 1. Table 1: Summary of the main selections that guarantee the exclusivity between individual final states. The symbol represents an electron or a muon; τ leptons are considered separately. The AK4 b jets are additional b tagged AK4 jets that do not geometrically overlap with AK8 jets. The symbol "-" implies that no selection is applied.

Ref.
Channel

Systematic uncertainties
The systematic uncertainties originating from the background estimation in the individual channels are considered uncorrelated, because the backgrounds are determined from statis-tically independent control regions. The uncertainties arising from the reconstruction and calibration are instead correlated among the different channels. These include the uncertainties in jet energy and resolution, and in the e, µ, and τ h lepton energy, reconstruction, and identification. The uncertainties in the identification of the SM bosons that decay to quarks are dominant in final states containing at least one such decay, and originate from the m j scale and resolution, the τ 21 selection and its extrapolation at large jet p T , and the b tagging. These uncertainties affect the signal distribution, selection efficiency, and induce event migration effects between search regions. The uncertainty covering the selection of events in the jet mass window of the Higgs boson is also included. The uncertainties in the proton-proton inelastic cross section [73], the integrated luminosity during the 2016 data-taking [74] and the kinematic acceptance of final-state particles, which affect the signal normalization, are also considered as correlated among channels. Theoretical uncertainties in the cross section and in the signal geometric acceptance related to the choice of PDFs used in the event generators [75], and the uncertainties in the factorization and renormalization scales, are evaluated according to the PDF4LHC recommendations [75]. The impact of these uncertainties on the signal cross section can be as large as 78%, depending on the signal mass and the initial state (qq or gg). A summary of the main systematic uncertainties is given in Table 2.

Statistical combination
No significant excess is observed above the SM background expectations. Upper limits are set at 95% confidence level (CL) [76,77] on the cross section of a heavy resonance, which is rescaled by a signal strength modifier parameter µ with uniform prior. Systematic uncertainties are represented by nuisance parameters θ, and affect both signal and background expectations [78]. Systematic uncertainties are considered fully (anti-)correlated when related to a common nuisance parameter, or uncorrelated when different nuisance parameters are used. The prior on these parameters is either flat (represented by the symbol "f" in Table 2) or lognormal distributed (identified with "s", "b", "t" in Table 2). The statistical procedure is based on a likelihood constructed as: where P represents the Poisson probability, and p j (θ j |θ j ) is the frequentist pdf of the nuisance parameter θ j and its default valueθ j associated with the j th uncertainty. The values s c,i (θ) and b c,i (θ) represent the number of signal and background events in the channel c and bin i, respectively. The test statistic is based on the profile likelihood ratioq: and the quantitiesμ,θ µ ,θμ are fixed to their best-fit value, and the range of µ is limited in the 0 ≤μ ≤ µ interval. The uncertainties that affect the signal normalization (PDFs and factorization and renormalization scales, marked with "t" in Table 2) are treated differently depending on how the exclusion is presented. When deriving upper limits on the cross section, these uncertainties are not varied in the fit, but are reported separately as the uncertainty of the theoretical cross sections from the model. When placing limits on the model parameters, these nuisance parameters are fixed at the best-fit values, in the same manner as to the other systematic uncertainties. Table 2: Summary of the main systematic uncertainties. The second column reports whether a systematic uncertainty is considered fully correlated or not across different channels. The third column indicates whether the uncertainty affects the yield, the shape of the distributions, or both, or if it induces migration (migr.) effects across search regions. The fourth column reports the smallest and largest effect of the uncertainty in either the yield or the signal shape parameters. The symbols "s", "b" indicate that the uncertainty affects the signal, the main backgrounds of the analysis, respectively. The treatment of non-dominant backgrounds is often different and not reported here. The symbol "f" indicates that the parameters are not constrained, or associated with large uncertainties as in the case of multi-dimensional fits. The entries labeled with "t" are treated differently depending on the interpretation of the exclusion limit, as discussed in Section 7. Uncertainties marked with "-" are not applicable or are negligible. Bkg. modeling no shape - The upper limit on µ is derived from the 95% CL s criterion, defined as CL s (µ) = p(µ)/(1 − p(0)), such that CL s (µ) = 0.05. The quantities p(µ) and 1 − p(0) represent the probabilities to have a value ofq equal to, or larger than the observed value in the signal or background (µ = 0) hypotheses, respectively, and are derived from analytical functions using the asymptotic approximation [79]. This approximation leads to limits that are up to 30% stronger in regions with a small number of data events, compared to those obtained from the Monte Carlo generation of pseudo-experiments [78].

Results and interpretation
The data are in agreement with the expected background from SM processes. The largest deviation from the expected limit is observed in the V model A at a mass of 1.3 TeV, with a local significance of 2.7 standard deviations, corresponding to a global significance of 1.6 standard deviations. The local significances are obtained in the asymptotic approximation [79], and the global significances are evaluated using the trial factors method [80].
The exclusion limit on the cross section of each diboson channel (WW, WZ, ZZ, WH, ZH, HH) is depicted in Fig. 1 for the combination of all contributing channels, reported in Table 1, according to the spin of the new resonance. The generated signal can be either a spin-1 heavy vector (W or Z , as in the HVT model) or a spin-2 boson (as in the bulk graviton model). In fact, the spin and polarization of the heavy resonance does affect the final state, signal acceptance, and selection efficiencies. The exclusion limits are not presented above 4.5 TeV, since at larger masses the background estimation procedure used in diboson analyses becomes less reliable because of the lack of events in data.    Fig. 2. In this scenario, the contribution of the dilepton channels is negligible because their branching fraction is of the order of a few permil. The contribution from VH decays to VV channels, caused by an underestimation of m j , is also considered. The predictions of the HVT model B are superimposed on the exclusion limits, showing that a W boson of mass below 4.3 TeV, and a Z boson with mass below 3.7 TeV are excluded at 95% CL.
The HVT hypothesis is tested in Fig. 3 [81], which excluded a triplet of heavy resonances with masses up to 2.4 TeV in the same model. The dilepton resonances provide the most stringent results within the HVT model A framework, and are combined with the diboson searches in Fig. 3. A heavy triplet of V resonances is excluded up to a mass of 5.0 TeV.
The exclusion limits on the resonance cross sections shown in Fig. 3 are also interpreted as limits in the [g H , g f ] plane of the HVT parameters. The excluded region of parameter space for narrow resonances obtained from the combination of all the channels is shown in Fig. 4. The dilepton and diboson searches constrain different regions of the parameter space, as the dilepton searches can probe the region where the coupling to the SM bosons approaches zero. In the triplet interpretation, the ratio of the W to Z cross sections is assumed to be determined by the ratio of the partonic luminosities, and to depend only weakly on the model parameters. The fraction of the parameter space where the natural width of the resonances is larger than the  average experimental resolution of 5%, and the narrow-width approximation is thus invalid, is also indicated in Fig. 4.
In the spin-2 bulk graviton model, the WW, ZZ, and HH channels are combined, setting upper limits of up to 1.1 fb on the cross section of a graviton with mass up to 4.5 TeV. In the κ = 0.5 scenario, a graviton with a mass smaller than 850 GeV is excluded at 95% CL, as shown in Fig. 5. Larger κ values increase the production cross sections, but also the graviton natural width, which can be comparable or larger than the experimental resolution. In these cases, the narrow-width approximation is no longer valid.
These results represent the most stringent limits on heavy vector models set by CMS and are comparable at large V mass to the limits obtained the ATLAS in the combination of similar channels [45]. At lower masses, the ATLAS combination excludes smaller cross sections be-cause of the inclusion of final states with three or more leptons. The exclusion limits in the bulk graviton model are not directly comparable because of the large κ = 1.0 parameter adopted by ATLAS, which increases the cross section but also implies a non-negligible natural width.   : Observed and expected 95% CL upper limit on the cross section of the spin-2 bulk graviton as a function of its mass for the statistical combination of the WW, ZZ, and HH channels. The inner green and outer yellow bands represent the ±1 and ±2 standard deviation variations on the expected limit. The solid curve and its shaded area represent the cross section derived with the parameter κ = 0.5 and the associated uncertainty.

Summary
A statistical combination of searches for heavy resonances decaying into pairs of vector bosons, a vector boson and a Higgs boson, two Higgs bosons, or pairs of leptons, has been presented. The results are based on data collected by the CMS experiment at √ s = 13 TeV during 2016 corresponding to an integrated luminosity of 35.9 fb −1 . In models with warped extra dimensions, upper limits of up to 1.1 fb are set at 95% confidence level on the production cross section of the spin-2 bulk graviton. For models with a triplet of narrow spin-1 resonances, heavy vector bosons with masses below 5.0 and 4.5 TeV are excluded at 95% confidence level in models where the W and Z bosons couple predominantly to fermions and bosons, respectively. In the latter, the statistical combination extends the exclusion limit by 700 GeV as compared to the best individual channel.          [27] ATLAS Collaboration, "Search for resonances in the mass distribution of jet pairs with one or two jets identified as b-jets in proton-proton collisions at √ s = 13 TeV with the ATLAS detector", Phys. Rev. D 98 (2018) 032016, doi:10.1103/PhysRevD.98.032016, arXiv:1805.09299.