Search for a heavy resonance decaying into a Z boson and a Z or W boson in 2$\ell$2q final states at $\sqrt{s}=$ 13 TeV

A search has been performed for heavy resonances decaying to ZZ or ZW in 2$\ell$2q final states, with two charged leptons ($\ell=$ e,$\mu$) produced by the decay of a Z boson, and two quarks produced by the decay of a W or Z boson. The analysis is sensitive to resonances with masses in the range from 400 to 4500 GeV. Two categories are defined based on the merged or resolved reconstruction of the hadronically decaying vector boson, optimized for high- and low-mass resonances, respectively. The search is based on data collected during 2016 by the CMS experiment at the LHC in proton-proton collisions with a center-of-mass energy of $\sqrt{s}=$ 13 TeV, corresponding to an integrated luminosity of 35.9 fb$^{-1}$. No excess is observed in the data above the standard model background expectation. Upper limits on the production cross section of heavy, narrow spin-1 and spin-2 resonances are derived as a function of the resonance mass, and exclusion limits on the production of W$'$ bosons and bulk graviton particles are calculated in the framework of the heavy vector triplet model and warped extra dimensions, respectively.


Introduction
The validity of the standard model (SM) of particle physics is corroborated by a wide set of precise experimental results with an impressive level of accuracy. Nonetheless, there are several open points where the SM fails to provide an explanation, either for experimental observations, as in the case of the presence of dark matter in the universe, or for theoretical questions, such as the omission of gravity from the SM, and the hierarchy problem.
Several SM extensions addressing the open questions of the SM predict the presence of new heavy particles with an enhanced branching fraction for decays into pairs of vector bosons. The existence of heavy spin-2 gravitons (G) is predicted in the Randall-Sundrum model with warped extra spatial dimensions (WED) [1][2][3]. In the bulk scenario [4,5], the main free parameters are the mass of the first Kaluza-Klein graviton excitation (the bulk graviton mass), and the ratio κ ≡ κ/M Pl , where κ is a curvature parameter of the WED metric and M Pl ≡ M Pl / √ 8π is the reduced Planck mass. The introduction of a spin-1 triplet of Z and W bosons is described in the heavy vector triplet (HVT) model [6], which generalizes a large number of explicit models in terms of a small set of parameters: c H , controlling the interactions of the triplet with the SM vector and Higgs bosons; c F , which describes the direct interaction with fermions; and g V , which represents the overall strength of the new vector boson triplet interactions.
A variety of searches for heavy resonances decaying to two vector bosons have been carried out in the past. The most recent results from the CERN LHC [7][8][9][10][11], with no evidence of signal, have provided stringent upper limits on signal cross sections in these models. This paper reports on the results of a search for heavy, narrow resonances (collectively indicated as X) decaying into 2 2q final states, with two charged leptons ( = e, µ) produced by the leptonic decay of a Z boson and a pair of quarks produced from the hadronic decay of a vector boson (V = W or Z). In the narrow-width assumption, the width of the heavy resonance is taken to be small in comparison to the experimental resolution. Two complementary search strategies are defined to span the mass range 400 < m X < 4500 GeV, where m X is the mass of the heavy resonance. The first strategy, referred to as the "high-mass analysis", is optimized for the range 850 < m X < 4500 GeV by selecting events where the vector bosons have a large Lorentz boost, resulting in the collimation of their decay products. The high-mass analysis uses dedicated leptonic reconstruction and identification techniques to reconstruct leptons in close proximity to each other in order to retain high signal efficiency, as well as jet substructure techniques to identify the hadronic decay of the W or Z boson into a pair of quarks contained in a single merged reconstructed jet. For lower resonance masses, the quarks produced by the hadronic decay of the V boson may be sufficiently separated to be reconstructed as two single narrow jets (dijet). A second strategy, referred to as the "low-mass analysis", is therefore defined in this regime, exploiting dijet reconstruction in addition to the reconstruction of merged jets to retain signal efficiency in the range 400 < m X < 850 GeV for those events in which no merged V candidate is found. To increase the signal sensitivity, in the low-mass analysis a categorization based on the flavor of the jets is used, to exploit the relatively large decay branching fraction of the Z boson to pairs of b quarks. This paper is organized as follows: in Section 2, a description of the data and simulated samples used in the analysis is provided; Section 3 briefly describes the CMS detector; Section 4 provides a description of the event reconstruction; in Section 5, the event selection is discussed; Section 6 contains the description of the signal and describes the estimation of the SM background; the systematic uncertainties affecting the analysis are presented in Section 7; and the results of the search for heavy spin-1 and spin-2 resonances are presented in Section 8. Finally, results are summarized in Section 9.

Data and simulated samples
This analysis uses data collected by the CMS detector during proton-proton (pp) collisions at the LHC at √ s = 13 TeV, corresponding to an integrated luminosity of 35.9 fb −1 . The events were selected online by criteria that require the presence of at least one electron or muon; these criteria are described in Sec. 5.
Simulated signal samples are used in the analysis to optimize the search for the potential production of heavy spin-1 or spin-2 resonances. For this purpose, signal samples are generated according to the HVT and WED scenarios, respectively. For both scenarios, the samples are generated at leading order (LO) in QCD with the MADGRAPH5 aMC@NLO 2.2.2 generator [12]. Two HVT models are considered as benchmarks, "model A" and "model B", with different values of the three defining parameters described earlier: for "model A", g V = 1, c H = −0.556, and c F = −1.316, while for "model B", g V = 3, c H = −0.976, and c F = 1.024.
Different resonance mass hypotheses are considered in the range from 400 to 4500 GeV. The resonance width is predicted to be between 0.4 and 2.3 GeV for a W candidate in HVT model A, and between 14 and 64 GeV for HVT model B, depending on the W mass hypothesis [6]; in the WED model with κ = 0.1, the bulk graviton signal width is predicted to range from 3.6 to 54 GeV [13]. Since the resonance width is small in comparison with the experimental resolution, for simplicity, the width is taken to be 1 MeV in the simulation. In the case of the spin-1 W , the resonance is forced to decay into one Z and one W boson; additionally, the Z boson is then forced to decay to a pair of electrons, muons, or tau leptons, while the W boson is forced to decay into a pair of quarks. The generated spin-2 bulk graviton is instead forced to decay into two Z bosons, one decaying leptonically into any pair of charged leptons, and the other Z boson decaying hadronically into a pair of quarks.
Several SM processes yielding final states with charged leptons and jets are sources of background events for the analysis, and corresponding Monte Carlo (MC) simulated samples have been generated and used in the analysis.
The SM production of a Z boson in association with quarks or gluons in the final state (Z + jets) represents the dominant background process for the analysis, having topological similarities to the signal because of the presence of a pair of charged leptons and jets. However, since the quark-and gluon-induced jets are not associated with the decay of a vector boson, the jet mass spectrum is characterized by a smooth distribution and the distribution of the 2 + jet system invariant mass falls exponentially, in contrast with the peaking distribution expected from the signal in both the jet and 2 + jet mass spectra. The Z + jets MC samples are produced with MADGRAPH5 aMC@NLO at next-to-leading order (NLO), using the FxFx merging scheme [14] between the jets from matrix element calculations and parton showers, and normalized to the next-to-NLO cross section computed using FEWZ v3.1 [15].
Another important source of SM background arises from processes leading to top quark production. Simulated samples describing the production of top quark pairs are generated with MADGRAPH5 aMC@NLO at LO, with the MLM matching scheme [16]. Single top quark production is also considered; sand t-channel single top quark samples are produced in the fourflavor scheme using MADGRAPH5 aMC@NLO and POWHEG v2 [17][18][19][20], respectively, while tW production is simulated at NLO with POWHEG in the five-flavor scheme [21]. Additional top quark background processes, such as the associated production of a Z or W boson with pairproduced top quarks, and the production of tqZ, are also considered in the analysis and produced at NLO with MADGRAPH5 aMC@NLO.
The SM diboson production of VV is an irreducible source of background for the analysis, since the jet mass spectrum will contain a peak from the hadronic decay of W and Z bosons, like the expected jet mass spectrum for the signal; however, this process produces a smoothly falling 2 + jet invariant mass distribution. The SM production of pairs of vector bosons (WW, WZ, and ZZ) is simulated at NLO with MADGRAPH5 aMC@NLO.
For all the simulated samples used in the analysis, the simulation of parton showering and hadronization is described by interfacing the event generators with PYTHIA 8.212 [22] with the CUETP8M1 [23] tune, while the parton distribution functions (PDFs) of the colliding protons are given by the NNPDF 3.0 [24] PDF set. Additional pp interactions occurring in the same or nearby bunch crossings (pileup) are added to the event simulation, with a frequency distribution adjusted to match that observed in data. All samples are processed through a simulation of the CMS detector using GEANT4 [25], and reconstructed using the same algorithms as those for the data collected.

The CMS detector
The central feature of the CMS apparatus is a superconducting solenoid of 6 m internal diameter, providing a magnetic field of 3.8 T. Within the solenoid volume are a silicon pixel and 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. The silicon tracker covers the pseudorapidity range |η| < 2.5, while the ECAL and HCAL cover the range |η| < 3.0. Forward calorimeters extend the coverage provided by the barrel and endcap detectors to |η| < 5.2. Muons are detected in gas-ionization chambers embedded in the steel flux-return yoke outside the solenoid, with detection planes made using three technologies: drift tubes, cathode strip chambers, and resistive-plate chambers.
A more detailed description of the CMS detector, together with a definition of the coordinate system used and the relevant kinematic variables, can be found in Ref. [26].

Event reconstruction
The event reconstruction is performed globally using a particle-flow (PF) algorithm [27], which reconstructs and identifies each individual particle with an optimized combination of information from the various elements of the CMS detector.
The reconstructed vertex with the largest value of summed physics-object p 2 T is taken to be the primary pp interaction vertex. The physics objects chosen are those that have been defined using information from the tracking detector. These objects include jets, the associated missing transverse momentum, which was taken as the negative vector sum of the transverse momentum (p T ) of those jets, and charged leptons.
In the silicon tracker, isolated charged particles with p T = 100 GeV and |η| < 1.4 have track resolutions of 2.8% in p T and 10 (30) µm in the transverse (longitudinal) impact parameter [28]. 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, corrected for zero-suppression effects and for the response function of the calorimeters to hadronic showers. The energy of neutral hadrons is obtained from the corresponding corrected ECAL and HCAL energies.
Electrons are required to be within the range |η| < 2.5 covered by the silicon tracker, and are reconstructed from a combination of the deposited energy of the ECAL clusters associated with the track reconstructed from the measurements determined by the inner tracker, and the energy sum of all photons spatially compatible with being bremsstrahlung from the electron track. The identification of electrons is based on selection criteria relying on the direction and momentum of the track in the inner tracker, its compatibility with the primary vertex of the event [27], and on observables sensitive to the shape of energy deposits along the electron trajectory. The momentum resolution for electrons with p T ≈ 45 GeV from Z → ee decays ranges from 1.7% to 4.5% [29]. It is generally better in the barrel region than in the endcaps, and also depends on the amount of bremsstrahlung emitted by the electron as it traverses the material in front of the ECAL.
Muons are reconstructed in the entire CMS muon system acceptance region of |η| < 2.4 by combining in a global fit the information provided by the measurements in the silicon tracker and the muon spectrometer. Candidate muons are selected using criteria based on the degree of compatibility of the inner track, which is reconstructed using the silicon tracker only, and the track reconstructed using the combination of the hits in both the tracker and spectrometer. Further reconstruction requirements include the compatibility of the trajectory with the primary vertex of the event, and the number of hits observed in the tracker and muon systems. The relative p T resolution achieved is 1.3-2.0% for muons with 20 < p T < 100 GeV in the barrel and better than 6% in the endcaps. The p T resolution in the barrel is better than 10% for muons with p T up to 1 TeV [30].
Both electrons and muons are required to be isolated from hadronic activity and other leptons in the event. An isolation variable is defined as the scalar sum of the p T of charged hadrons originating from the primary vertex, plus the scalar sums of the transverse momenta for neutral hadrons and photons, in a cone of ∆R = √ (∆η) 2 + (∆φ) 2 < 0.3 (0.4) around the electron (muon) direction corrected to account for the contribution from neutral candidates originating from pileup, where φ is the azimuthal angle in radians. In the high-mass analysis, a specific muon isolation requirement is implemented to retain signal efficiency for high resonance masses, where the large Z boson boost may result in extremely close pairs of muons. For this reason, muon candidates in the high-mass analysis are retained if they pass an isolation requirement based on the sum of reconstructed p T of all tracks within ∆R < 0.3 from the muon trajectory, ignoring tracks associated with other reconstructed muons.
Hadron jets are clustered from particles reconstructed by the PF algorithm using the infraredand collinear-safe anti-k T algorithm [31,32] with distance parameters of 0.4 (AK4 jets) and 0.8 (AK8 jets). The jet momentum is determined as the vectorial sum of all constituent particle momenta. Contamination from pileup is suppressed using charged hadron subtraction (CHS) which removes from the list of PF candidates all charged particles originating from vertices other than the primary interaction vertex of the event. The residual contribution from neutral and charged particles originating from pileup vertices is removed by means of an event-byevent jet-area-based correction to the jet four-momentum. Identification requirements, based on the estimation of the energy fraction carried by the different types of PF candidates clustered into a jet, along with the multiplicity of the PF candidates, are used to remove jets originating from calorimetric noise. Corrections to the jet energy are derived from the simulation, and are confirmed with in situ measurements with the energy balance of dijet, multijet, photon + jet, and leptonically decaying Z + jet events [33].
A jet grooming technique is used for AK8 jets in this analysis to help identify and discriminate between jets from boosted hadronic V decays, which we refer to as "merged jets", and jets from quarks and gluons. The AK8 jets are groomed by means of the modified mass drop tagger algorithm [34], also known as the soft drop algorithm, with angular exponent β = 0, soft cutoff threshold z cut < 0.1, and characteristic radius R 0 = 0.8 [35]. The soft drop algorithm does not fully reject contributions from the underlying event and pileup. The mass of the AK8 jet (m j ) is therefore defined as the invariant mass associated to the four-momentum of the soft drop jet, after the application of the pileup mitigation corrections provided by the pileup per particle identification (PUPPI) algorithm [36].
Discrimination between AK8 jets originating from vector boson decays and those originating from gluons and quarks is also achieved by the N-subjettiness jet substructure variable [37]. This observable exploits the distribution of the jet constituents found in the proximity of the subjet axes to determine if the jet can be effectively subdivided into a number N of subjets. The generic N-subjettiness variable τ N is defined as the p T -weighted sum of the angular distance of all the k jet constituents from the closest subjet: The normalization factor d 0 is defined as d 0 = ∑ k p T,k R 0 , with R 0 the clustering parameter of the original jet. In this analysis, which aims to select V → qq ( ) decays, the variable that best discriminates V boson jets from those from quarks and gluons is the ratio of the 2-subjettiness to the 1-subjettiness: τ 21 = τ 2 /τ 1 . The τ 21 observable is calculated for the jet before the grooming procedure, and includes the PUPPI algorithm corrections for pileup mitigation.
For the identification of jets originating from the hadronization of bottom quarks, the combined secondary vertex (CSVv2) algorithm [38,39] is used, either directly on the AK4 jets or on the AK8 soft drop subjets with CHS pileup mitigation applied.
Only AK4 and AK8 jets reconstructed centrally in the detector acceptance, within |η| < 2.4, are considered in the analysis.

Event selection
Events are selected online by requiring the reconstruction at trigger [40] level of at least one charged lepton. For the high-mass analysis, p T thresholds of 115 (50) GeV are used for electrons (muons). No isolation requirements are applied at trigger level, to retain efficiency for high-mass signals, where the large boost expected for the leptonically decaying Z boson will cause the two charged leptons to be collimated in the detector. For the low-mass analysis a larger separation between the leptons is expected because of the lower p T of the Z boson, and isolation requirements are included in the trigger selection, allowing the use of lower lepton p T thresholds. The online selection for the low-mass analysis requires at least one electron with p T > 25 GeV and |η| < 2.1 passing tight identification and isolation requirements, or at least one muon with p T > 24 GeV and |η| < 2.4, subject to loose identification and isolation requirements, using the variables described in Ref. [40].
To reconstruct the Z boson candidate, at least two well-identified leptons with opposite charge and the same flavor are required to be present in the event. The leading lepton in the event is required to pass more stringent selection requirements than the online thresholds to avoid inefficiencies induced by the trigger selections. In the high-mass analysis, the leading (subleading) lepton is required to have p T > 135 (35) GeV for electrons, and p T > 55 (20) GeV for muons. Loose isolation and identification requirements are applied to the leptons to retain high signal efficiency. For electrons, we use a set of requirements that have been observed to have an efficiency of about 90% for both low and high mass points. For muons, as the CMS standard requirements [41] only have an efficiency of about 65% for close muons, we instead use a dedicated selection where one of the two muons is allowed to be identified only in the tracker. The isolation variable is calculated removing the contribution of the other muon if it falls within the isolation cone, therefore recovering a signal efficiency of about 90% for high mass resonances. For the low-mass analysis, the leading (subleading) lepton is required to have p T larger than 40 (30) GeV and to fall in the range |η| < 2.1 (2.4).
The selection of the Z boson candidate relies on the invariant mass of the dilepton pair, m . This is required to satisfy 70 < m < 110 GeV, except for the low-mass analysis in the resolved category (discussed below) where the requirement is 76 < m < 106 GeV to enhance the sensitivity to the signal by reducing the nonresonant contribution in the sample with b tagged jets.
Different strategies are used in the low-and high-mass analyses to identify and reconstruct the hadronically decaying V boson, as described below, to cope with the different V boson boost regimes expected for low-and high-mass signal candidates.
In the high-mass analysis a merged jet is required in the event, and its mass m j is used to select the hadronically decaying W or Z. The signal is expected to be almost fully contained in the mass range 65 < m j < 105 GeV, which is thus defined as the signal region (SR). In order to select candidate signal events, where a heavy massive particle decays into a pair of boosted vector bosons, both the dilepton pair and the leading jet selected in an event are required to have p T > 200 GeV; this is motivated by the p T spectrum of the V bosons observed in simulation. Events are divided into categories depending on the flavor of the charged leptons (e or µ) and the value of the jet τ 21 variable. As the signal is expected to have lower values of τ 21 , two different purity categories are defined: events with τ 21 < 0.35 are defined as the high-purity (HP) category, while events with 0.35 < τ 21 < 0.75 fall into a low-purity (LP) category, used to retain some sensitivity to signal although a larger amount of background is expected with respect to the HP category. The τ 21 > 0.75 region is expected to be dominated by the background, and is therefore not used in the high-mass analysis. In total, four exclusive categories (from the two purity and two lepton flavor categories) are defined for the high-mass analysis.
In the low-mass analysis, events are divided into two categories depending on whether the two quarks from the hadronic V decay merge into a single reconstructed jet or can be resolved as two distinct jets. In the merged category, merged jets with p T > 200 GeV and τ 21 < 0.40 are selected. The choice of a looser τ 21 selection with respect to the cutoff applied in the HP category of the high-mass analysis is driven by the higher expected signal efficiency for merged events, which are selected in the low-mass analysis using only one τ 21 category. As in the high-mass analysis, the jet mass is required to be in the range 65 < m j < 105 GeV for the jet to be considered a candidate W or Z boson, which is also defined as the SR for the merged low-mass analysis. The resolved category contains events that do not contain a merged V candidate, but instead two AK4 jets, both with p T > 30 GeV that form a dijet candidate with invariant mass m jj > 30 GeV and p T > 100 GeV. In both the merged and resolved cases, the p T selection is determined by comparing the p T spectrum of simulated signal events with the expected background. Both the merged and resolved categories are further split into two b tag categories. Events in the merged tagged category are required to have at least one subjet satisfying a b tagging requirement corresponding to ≈65% efficiency for b quark identification and ≈1% light-flavor jet mistag rate; events not passing this requirement are placed in the merged untagged category. For the resolved tagged category, events are required to have at least one jet satisfying the same b tagging requirement used in the merged category; a looser b tag selection is instead required for the other jet, with ≈80% efficiency and ≈10% light-flavor jet mistag rate. Events failing these requirements fall in the so-called resolved untagged category. An arbitration procedure is used to select the dijet candidate in case of events containing more than two selected narrow jets: first, if a dijet passing the b tagging requirements is selected in the event, the candidate in the b tag category is chosen; then the dijet candidate closest in mass to the Z boson mass is selected as the candidate V boson. The signal region for the low-mass resolved category accepts events in the dijet mass range 65 < m jj < 110 GeV. Eight categories are defined in the low-mass analysis, based on the lepton flavor, the b-tag category, and the merged or resolved reconstruction of the hadronically decaying V candidate.
The τ 21 and merged jet p T distributions of the V candidate for events selected in the merged category of the low-mass analysis are shown in Fig. 1, where the m j and m jj distributions for events with a V candidate are also shown for the merged and resolved low-mass analysis categories, respectively.
6 Background estimation

High-mass analysis
The main source of background events in the final state of the analysis arises from the production of a leptonically decaying Z boson in association with quark and gluon jets. A second background source relevant for the analysis is SM diboson production, mainly ZZ and ZW, with a leptonically decaying Z boson together with a W or Z boson decaying hadronically. These diboson events are an irreducible background for the analysis, as the mass distribution of the SM V jet peaks in the same region as the signal. Finally, top quark production is considered as a source of background in the analysis, despite having a much smaller contribution with respect to other SM backgrounds in the region probed by this analysis, mostly because of the Z boson invariant mass selection and the large boost required in the event.
All SM background processes are characterized by a smoothly falling distribution of the invariant mass of the dilepton pair and the jet selected (m VZ ), whereas the signal is instead expected to appear as a narrow peak at a value of m VZ close to the actual value of the mass of the resonance m X .
To minimize the dependency on the accuracy of the simulation, the contribution of the dominant background, Z + jets SM production, is estimated using data. Two signal-depleted regions are defined by selecting events with jet mass outside the m j signal mass window defined in Section 5; these are the sideband (SB) regions. A lower sideband (LSB) region is defined for events with 30 < m j < 65 GeV, close to the SR of the analysis, while a higher sideband (HSB) region contains events with 135 < m j < 300 GeV. The region 105 < m j < 135 GeV is not used in the analysis, to exclude events containing the hadronic decays of a SM Higgs boson, which are targeted in other CMS analyses, such as that described in [42].
The Z + jets background m VZ shape and normalization are obtained by extrapolation from fits to data in the SB regions.
The m j distribution for the SM background sources considered in the analysis is modeled by means of analytic functions describing the spectrum of each background in the mass region 30 < m j < 300 GeV. In the LP category, the m j spectrum in Z + jets events is described by a smoothly falling exponential distribution, while a broad structure centered around the mass of the W boson present in the HP category is modeled with an error function convolved with an exponential distribution, which is of particular importance for describing the behavior at large values of m j . The peaking structure of the diboson background, originating from the presence of a jet from a genuine W or Z boson in the event, is described in both the LP and HP categories with a Gaussian distribution. The remaining component of the distribution, consisting of tails extending far from the SR, is modeled in the LP category with an exponential function, similarly to the Z + jets case. In the HP category, the VV events are mostly contained in the SR, and the small fraction of events present in the Higgs boson and LSB regions is described with an additional broad Gaussian contribution. The top quark background (tt, single top quark, tZq, and ttV production) is mostly similar in shape to the Z + jets background; in the LP category, in addition to the exponentially falling component, a Gaussian is included to model the top quark peak appearing in the HSB for m j ≈ 170 GeV.
The expected yield of the Z + jets background in the SR is extracted by a fit of the m j distribution in the SBs taking into account all background contributions. The parameters describing the m j shape and normalization of the subdominant background processes are fixed to those extracted from the simulation. All the parameters used to describe the Z + jets contribution are left free to float in the fit to the data SBs. Alternative functions modeling the m j shape of the main Z + jets background are used to evaluate the impact of the function choice on the signal normalization.
The m j distribution for expected and observed events is shown in Fig. 2. The m j distributions of the events in data, compared to the expected background shape, for the high-mass analysis in the electron (upper) and muon (lower) channels, and for the high-purity (left) and low-purity (right) categories. The expected background shape is extracted from a fit to the data sidebands (Z + jets) or derived from simulation ("top quark" and "VV"). The dashed region around the background sum represents the uncertainty in the Z + jets distribution, while the dashed vertical region ("Higgs") shows the expected SM Higgs boson mass range, excluded from the analysis. The bottom panels show the pull distribution between data and SM background expectation from the fit, where σ data is the Poisson uncertainty in the data.
To describe the shape of the m VZ variable for the Z + jets background in the SR, the following transfer function is defined from simulation: where f

MC,Z+jets SR
(m VZ ) and f

MC,Z+jets SB
(m VZ ) are the probability density functions describing the m VZ spectrum in the SR and SBs, respectively, of the simulated Z + jets sample.
The shape of the Z + jets background in the SR is then extracted from a simultaneous fit to data in the SBs, and to simulation in both the SR and SBs, to correct the functional form obtained from data using the α(m VZ ) ratio. The m VZ shape is described by two-parameter exponential functions for both data and simulation. The final estimate of the background m VZ shape predicted in the SR is then given by the following relation: where N The α(m VZ ) function accounts for differences and correlations in the transfer process from the SB regions to the SR, and is largely unaffected by uncertainties in the overall Z + jets cross section and distribution shapes.
The final m VZ spectra in the SR are shown in Fig. 3, compared to the expected estimated background.
The validity and robustness of the background estimation method is demonstrated by the agreement observed between the shape and normalization for events selected in an intermediate m j mass region (50 < m j < 65 GeV), corresponding to the part of the LSB shown in Fig. 2 above 50 GeV, and the prediction made using the events in the remaining part of the LSB and the full HSB regions.
The description of the signal m VZ shape is extracted from simulated signal samples. Several signal samples generated with resonance mass ranging from 400 to 4500 GeV in the narrow width approximation are modeled independently for each channel with a Crystal Ball (CB) function [43]. The power-law component of the CB function improves the description of the m VZ signal distribution by accounting for the small contribution from lower m VZ tails appearing for high signal masses. The resolution of the reconstructed m VZ can be extracted from the Gaussian core width of the CB function, and is estimated to be 2-3.5% in the electron channel and 3-4% in the muon channel, depending on the mass of the resonance.

Low-mass analysis
For the low-mass analysis, the Z + jets background is characterized using simulated Drell-Yan + jets events. Because of the limited number of simulated events, the m VZ distributions in the btagged categories are susceptible to sizable statistical fluctuations, which affect the quality of the background modeling. It has been observed, however, that within simulation uncertainties, the Z + jets mass shape is the same for events with and without b-tagged jets. Therefore, the Z + jets shape in the b-tagged category is described using the m VZ shape obtained from the simulation without making any b tag requirements.
Sideband regions are defined depending on the mass of the hadronic V boson candidate. The mass ranges 30 < m j < 65 GeV and 135 < m j < 180 GeV are used for the merged category, whereas for the resolved event selection the upper mass threshold is raised to 300 GeV to take advantage of the increased number of events in that region.
In the final fit to the data, the Z + jets background normalization in the SR is constrained by the observed yield in the SBs; this procedure is applied independently to each category. The shape predictions from the NLO Z + jets simulation are taken as a baseline m VZ shape in the SR of every category; additionally, a family of linear correction functions: with individual members of the family defined by the slope parameter s, is considered. Figure 4 shows fits to the SB m VZ distributions where the slope parameter s, allowed to float freely, is constrained by the observed shapes in data. The two-standard-deviation uncertainties in the fitted linear correction functions, which are in the range from 2 × 10 -4 to 6 × 10 -4 GeV −1 , depending on the category, are observed to cover the residual shape differences in the SBs.
In the signal region fit of each category, the SB-constrained slope parameter s is treated as a Z + jets shape systematic effect. In this way the background shape can be corrected to that observed in data. Statistical uncertainties associated with the simulated Z + jets distributions are also taken into account in the fit. The fits in the merged V categories include the peaking region of the background; Fig. 4 shows that the SB data in this particular region are described well by the fit.
Dilepton backgrounds that do not contain a leptonic Z boson decay are estimated from data using eµ events passing the analysis selection. This approach accounts for tt production, WW + jets, Z → ττ + jets, single top quark, and hadrons misidentified as leptons, which we collectively refer to as t + X. The relative yield of ee and µµ events with respect to eµ events has been estimated on a top quark-enriched control sample and shown to be consistent with expectations. Also, the eµ m VZ distribution was compared with the prediction from simulated background events with symmetric lepton flavor, and found to be in agreement. The contribution of this t + X background is 2% and 20% of the total background in the untagged and tagged categories of the resolved analysis, respectively. The merged analysis has a t + X contribution of 0.5% and 1% in the untagged and tagged categories, respectively.
The diboson background (ZZ and ZW, with Z → ) is estimated directly from simulation. The contribution from these events represents 4% and 5% of the total background in the untagged and tagged categories of the resolved analysis, respectively, while in the merged analysis it is about 14% and 16% in the untagged and tagged categories, respectively.
The m VZ distributions for the signal region for the merged and resolved categories are depicted in Fig. 5.

Systematic uncertainties
Several sources of systematic uncertainties influence both the normalization and shape of the backgrounds and signal distributions in the analysis. In the high-mass analysis, where the Z + jets background component is estimated with data, the main systematic uncertainties in the predicted normalization for the Z + jets background arise from the statistical uncertainties in the fit of the m j sidebands in data. Another uncertainty affecting the normalization of the main background is evaluated by taking the difference between the expected Z + jets contribution in the SR obtained by the main function used to describe the m j spectrum, and an alternative function choice. An additional normalization uncertainty is related to the choice of the function used to describe the m j spectrum for the subdominant top quark and VV backgrounds, evaluated from simulation, and propagated to the Z + jets normalization prediction in the SR. Overall, the Z + jets normalization uncertainties contribute from 9 to 15%, depending on the category. The main shape uncertainties in the Z + jets background are extracted from the covariance matrix of the fit to the m VZ data SB spectrum, convolved with the uncertainties provided by the α(m VZ ) ratio, via the simultaneous fit procedure described in Section 6.1. In the low-mass analysis, to account for background shape systematic effects not explicitly evaluated, data and simulation are compared in the sideband region, and the residual shape difference is treated as an additional uncertainty, resulting in the dominant background shape systematic uncertainty of the low-mass analysis.
The top quark and VV background components have a systematic uncertainty in the normalization arising from the degree of knowledge of the respective process production cross sections. The value of the VV production cross section, taken from a recent measurement by the CMS Collaboration [44,45], is assigned an uncertainty of 12%. The top quark background uncertainties are estimated differently in the low-and high-mass analyses: in the low-mass analysis, where a dedicated eµ control region is exploited to measure the t + X background normalization, a 4% uncertainty is estimated by comparing the yield of eµ events with ee + µµ data; in the high-mass analysis, where the top quark production is taken from simulation, a 5% uncertainty in the cross section is used, which is extracted from the recent CMS measurement of top quark pair production in dilepton events [46].
Uncertainties associated with the description in simulation of the trigger efficiencies, as well as the uncertainties in the efficiency for electron and muon reconstruction, identification, and isolation, are extracted from dedicated studies of events with leptonic Z decays, and amount to 1.5-3%, depending on the lepton flavor. The uncertainties in the lepton momentum and energy scales are taken into account, and propagated to the signal shapes and normalization, with a typical impact on the normalization of about 0.5-2%, depending on the lepton flavor.
Uncertainties in the jet energy scale and resolution [47] affect both the normalization and the shape of the background and signal samples. The momenta of the reconstructed jets are varied according to the uncertainties in the jet energy scale, and the selection efficiencies and m VZ signal shapes are reevaluated using these modified samples, resulting in a change of 0.1 to 1.8%, depending on the jet selection. The impact of the jet energy resolution is also propagated, and a smaller impact is observed compared with that due to the uncertainty in the energy scale.
The dominant uncertainty in the signal selection efficiency is the uncertainty in the V boson identification efficiency, corresponding to 11% (23%) for the HP (LP) category in the high-mass analysis, and 6% for the merged category of the low-mass analysis [48]. The V boson identification efficiency, the groomed mass resolution of V jets, and the related systematic uncertainty are measured in data and simulation in an almost pure selection of semileptonic tt events where boosted W bosons produced in the top quark decays are separated from the combinatorial tt background by means of a simultaneous fit to the soft drop mass. The uncertainties in the soft drop mass scale and resolution are propagated to the groomed jet mass, and the impact on the expected selection efficiency of signal and VV background is taken into account. An additional uncertainty affecting the signal normalization is included to account for the extrapolation of the uncertainties extracted from a tt sample at typical jet p T of 200 GeV to higher regimes, estimated from the differences between PYTHIA 8 and HERWIG ++ [49] showering models, yielding an uncertainty from 2.5 to 20% depending on the category. For the high-mass analysis, the uncertainties in the V boson identification efficiency and the extrapolation are treated as anticorrelated between the low-and high-purity categories.
For the low-mass analysis, one of the largest signal selection uncertainties is the uncertainty in the b tagging efficiency for the tagged categories of the analysis. The b tagging efficiencies and their corresponding systematic uncertainties are measured in data using samples enriched in b quark content, and their propagation to the signal region of the low-mass analysis produces an uncertainty of up to 4.3%. The uncertainties in the mistag efficiency are also considered; the uncertainties in the b tagging and mistag efficiencies are treated as anticorrelated between the tagged and untagged categories.
The impact of the uncertainties in the factorization and renormalization scales is propagated both to the normalization and the m VZ shapes for signal, and for the high-mass analysis to top quark and VV backgrounds. The corresponding scales are varied by a factor of 2 to measure the effect, resulting in an uncertainty of 2% for the diboson background normalization and 15% for top quarks. The impact on the signal acceptance is evaluated to be 0.1-3%, depending on the resonance mass and analysis category.
A systematic uncertainty associated with the choice of the set of PDFs used to generate the simulated samples is evaluated by varying the NNPDF 3.0 PDF set within its uncertainties, and its effect is propagated to both the signal and background m VZ shapes and normalization, resulting in a measured uncertainty of approximately 1%.
Additional systematic uncertainties affecting the normalization of backgrounds and signal from the contributions of pileup events and the integrated luminosity [50] are also considered and are reported in Table 1, together with the complete list of uncertainties considered in the analysis. In the high-mass analysis, the typical total uncertainty in the background normalization is in the range 10-60%, depending on the signal mass, and it is 1-5%, depending on the category, in the low-mass analysis. Table 1: Summary of systematic uncertainties, quoted in percent, affecting the normalization of background and signal samples. Where a systematic uncertainty depends on the resonance mass (for signal) or on the category (for background), the smallest and largest values are reported in the table. In the case of a systematic uncertainty applying only to a specific background source, the source is indicated in parentheses. Systematic uncertainties too small to be considered are written as "<0.1", while a dash (-) represents uncertainties not applicable in the specific analysis category.

Results and interpretation
Results are extracted separately for the high-and low-mass analyses from a combined maximum likelihood fit of signal and background to the m VZ distribution, simultaneously in all the categories used in the respective analysis. An unbinned fit is performed in the high-mass analysis, while a binned fit is performed in the low-mass one; this choice is determined by the fact that in the high-mass analysis, the signal and background shapes are described with analytical functions, while in the low-mass analysis, the background shapes are described by binned histograms. The systematic uncertainties discussed in Section 7 are included as nuisance parameters in the maximum likelihood fit, and the background-only hypothesis is tested against the combined background and signal hypothesis [51,52].
The largest excess of events with respect to the background-only hypothesis, with a local significance of 2.5 standard deviations, is observed in the vicinity of m X ≈ 1.2 TeV, and arises predominantly from a localized excess of events in the dimuon HP category of the high-mass analysis.
The limit at 95% confidence level (CL) on the signal cross section for the production of a heavy spin-1 or spin-2 resonance is set using the asymptotic modified frequentist method (CL s ) [51][52][53][54].
The results of the low-and high-mass analyses should agree for the intermediate mass range 800-900 GeV, which is accessible to both strategies with similar expected efficiencies for signal candidates. The results of the analysis are therefore presented based on the low-mass strategy up to resonance masses m X ≤ 850 GeV, and based on the high-mass analysis for m X ≥ 850 GeV. At the intermediate mass point m X = 850 GeV, the results of both strategies are presented, and the expected limits at 95% CL of the low-and high-mass analyses on the signal cross sections are found to be in agreement within 3 and 6% for the W and bulk graviton signal model, respectively.
The observed upper limits on the resonance cross section, multiplied by the branching fraction for the decay into one Z boson and a W or Z boson, σ W B(W → ZW) or σ G B(G → ZZ), are reported as a function of the resonance mass in Fig. 6 assuming a W or G produced in the narrow-width approximation, and the local p-value [55] is shown in Fig. 7.  Figure 6: Observed and expected 95% CL upper limit on σ W B(W → ZW) (left) and σ G B(G → ZZ) (right) as a function of the resonance mass, taking into account all statistical and systematic uncertainties. The electron and muon channels and the various categories used in the analysis are combined together. The green (inner) and yellow (outer) bands represent the 68% and 95% coverage of the expected limit in the background-only hypothesis. The dashed vertical line represents the transition from the low-mass to the high-mass analysis strategy. Theoretical predictions for the signal production cross section are also shown: (left) W produced in the framework of HVT model A with g v = 1 and model B with g v = 3; (right) G produced in the WED bulk graviton model with κ = 0.5.

Summary
A search for a heavy resonance decaying into a Z boson and a Z or a W boson in 2 2q final states has been presented. The data analyzed were collected by the CMS experiment in proton-proton collisions at √ s = 13 TeV during 2016 operations at the LHC, corresponding to an integrated luminosity of 35.9 fb −1 . The final state of interest consists of a Z boson decaying leptonically into an electron or muon pair, and the decay of an additional W or Z boson into a pair of quarks. Two analysis strategies, dedicated to the low-and high-mass regimes (below and above 850 GeV, respectively), have been used to set limits in the range of resonance mass from 400 to 4500 GeV. Depending on the resonance mass, expected upper limits of 3-3000 and 1.5-400 fb have been set on the product of the cross section of a spin-1 W and the ZW branching fraction, and on the product of the cross section of a spin-2 graviton and the ZZ branching fraction, respectively.

Acknowledgments
We congratulate our colleagues in the CERN accelerator departments for the excellent performance of the LHC and thank the technical and administrative staffs at CERN and at other CMS institutes for their contributions to the success of the CMS effort. In addition, we gratefully acknowledge the computing centers and personnel of the Worldwide LHC Computing Grid for delivering so effectively the computing infrastructure essential to our analyses. Finally, we acknowledge the enduring support for the construction and operation of the LHC and the CMS detector provided by the following funding agencies: BMWFW and FWF (Aus