Search for high-mass diphoton resonances in proton–proton collisions at 13 TeV and combination with 8 TeV search

we acknowledge the enduring support for the construc-tion and operation of the LHC and the CMS detector provided by the following funding agencies: BMWFW and FWF (Austria); FNRS and FWO (Belgium); CNPq, CAPES, FAPERJ, and FAPESP (Brazil); MES (Bulgaria); CERN; CAS, MOST, and NSFC (China); COLCIEN-CIAS (Colombia); MSES and CSF (Croatia); RPF (Cyprus); SENESCYT (Ecuador); MoER, ERC IUT and ERDF (Estonia); Academy of Fin-land, MEC, and HIP (Finland); CEA and CNRS/IN2P3 (France); BMBF, DFG, and HGF (Germany); GSRT (Greece); OTKA and NIH (Hun-gary); DAE and DST (India); IPM (Iran); SFI (Ireland); INFN (Italy); MSIP and NRF (Republic of Korea); LAS (Lithuania); MOE and UM (Malaysia); BUAP, CINVESTAV, CONACYT, LNS, SEP, and UASLP-FAI (Mexico); MBIE (New Zealand); PAEC (Pakistan); MSHE and NSC (Poland); FCT (Portugal); JINR (Dubna); MON, RosAtom, RAS and RFBR (Russia); MESTD (Serbia); SEIDI and CPAN (Spain); Swiss Funding Agencies (Switzerland); MST (Taipei); ThEPCenter, IPST, STAR and NSTDA (Thailand); TUBITAK and TAEK (Turkey); NASU and SFFR (Ukraine); STFC (United Kingdom); DOE and NSF (USA).


Introduction
The standard model (SM) of particle physics has been highly successful in describing physical phenomena but it is widely considered to be an incomplete theory because of various shortcomings. In particular, the SM suffers from the so-called hierarchy problem [1], which refers to the large difference between the Higgs boson mass of 125 GeV [2] and the highest energy scale up to which the SM must be valid. Many extensions to the SM have been proposed to address the hierarchy problem, including theories with additional space-like dimensions [3] and models with extended Higgs boson sectors [4]. Some of these extensions predict new resonances that decay to a diphoton final state. For example, the Randall-Sundrum (RS) approach [3,5] to extra dimensions postulates massive excitations of spin-2 gravitons that can decay to two photons. A simple extension of the SM Higgs boson sector consists of the addition of a doublet of complex scalar fields. In such models [6], some of these additional scalar resonances can decay to a photon pair [7]. According to the Landau-Yang theorem, the spin of a resonance decaying to two photons can only be zero or an integer larger than one [8,9].

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). The tracking detectors cover the pseudorapidity range |η| < 2.5. The ECAL and HCAL, each composed of a barrel and two endcap sections, cover |η| < 3.0, with the boundary between the barrel and endcaps at around |η| = 1.5. Forward calorimeters extend the coverage to |η| < 5.0. The ECAL consists of 75 848 lead tungstate crystals. The barrel section has a granularity η × φ = 0.0174×0.0174, with φ the azimuthal angle, while the endcap sections have a granularity that coarsens progressively up to η × φ = 0.05×0.05.
Preshower detectors consisting of two planes of silicon sensors interleaved with a total of 3 X 0 of lead are located in front of the endcap sections. Muons are measured within |η| < 2.4 by gasionization detectors embedded in the steel flux-return yoke outside the solenoid. A more detailed description of the CMS detector, together with a definition of the coordinate system and the relevant kinematic variables, can be found in Ref. [30].
In the barrel section of the ECAL, for photons with energies on the scale of tens of GeV, an energy resolution of about 1% is achieved for unconverted photons and for photons that convert "late", i.e., just before entering the ECAL. The remaining barrel photons have an energy resolution of about 1.3% up to |η| = 1.0, rising to about 2.5% for |η| = 1.4. In the endcaps, the corresponding resolution for unconverted and late-converting photons is about 2.5%, while the remaining endcap photons have a resolution between 3% and 4% [31].
The particle-flow algorithm [32,33] reconstructs and identifies each individual particle with an optimized combination of information from the various elements of the CMS detector. Particle candidates are classified as either muons, electrons, photons, τ leptons, charged hadrons, or neutral hadrons.
A two-stage trigger system selects events of interest for the analysis. The level-1 trigger, composed of custom hardware processors, selects events at a maximum readout rate of about 100 kHz using information from the calorimeters and muon detectors. The high-level trigger software algorithms use the full event information to reduce the event rate to less than 1 kHz before data storage.

Event simulation
The pythia 8.2 [34] event generator with NNPDF2.3 [35] parton distribution functions (PDFs) is used to produce simulated signal samples of spin-0 and spin-2 resonances decaying to two photons. The samples are generated at leading order (LO), with values of the resonance mass m X in the range 0.5 < m X < 4.5 TeV. Three values of the relative width X /m X are used as benchmarks: 1.4 × 10 −4 , 1.4 × 10 −2 , and 5.6 × 10 −2 , where X is the width of the resonance. These relative widths correspond, respectively, to resonances much narrower than, comparable to, and significantly wider than the detector resolution. In the context of the RS graviton model, for which X /m X = 1.4 k2 [36], the relative widths correspond to the dimensionless coupling parameter k = 0.01, 0.1, and 0.2. The scalar resonances are produced through gluon-gluon fusion, and RS graviton resonances through both gluon-gluon fusion and quark-antiquark annihilation. In the RS model, the first mechanism accounts for approximately 90% of the production cross section.
The SM background mostly arises from the direct production of two photons, the production of γ + jets events in which jet fragments are misidentified as photons, and the production of multijet events with misidentified jet fragments. These backgrounds are simulated with the sherpa 2.1 [37], MadGraph5_aMC@NLO 2.2 [38] (interfaced with pythia 8.2 for parton showering and hadronization), and pythia 8.2 generators, respectively, using the CT10NLO [39], NNPDF3.0 [40], and NNPDF2.3 PDF sets, again respectively. The pythia tune CUETP8M1 [41] is used. For both the signal and background samples, the detector response is simulated using the Geant4 package [42]. The simulated samples incorporate additional pp interactions within the same or a nearby bunch crossing (pileup) and are weighted to reproduce the measured distribution of the number of interactions per bunch crossing. The average number of interactions per bunch crossing is 18, with an RMS of 4.

Event selection and diphoton mass spectrum
The trigger requirements, photon identification criteria, and event selection procedures are described in Ref. [11]. Some details are given below. Energy deposits in the ECAL compatible with the shower shape expected for a photon are clustered together to define a photon candidate. Variations in the crystal transparency during the data collection period are corrected for using a dedicated monitoring system, and the single-channel response is equalized based on collision data [31]. A multivariate regression technique [31] is used to correct the photon energy for the incomplete containment of the shower in the clustered crystals, the shower losses for photons that convert before reaching the calorimeter, and the effects of pileup. The interaction vertex is selected using the algorithm described in Ref. [43], which combines information on the correlation between the diphoton system and the recoiling tracks, the average transverse momentum (p T ) of the recoiling tracks, and, when available, directional information from reconstructed photon conversions. For resonances with a mass above 500 GeV, the fraction of events in which the interaction vertex is correctly assigned is approximately 90%. For each photon candidate, the transverse size of the electromagnetic cluster in the η coordinate must be compatible with that expected for a photon from a hard interaction, and the ratio of the associated energy in the HCAL to the photon energy must be less than 0.05.
Photon candidates are required to have p T > 75 GeV and to be reconstructed within |η| < 2.5. Candidates in the transition region between the barrel and endcap detectors (1. 44 where the reconstruction efficiency is not well described by the simulation, are rejected. Photon candidates associated with electron tracks that are incompatible with conversion tracks are rejected [31]. Photon candidates are required to be isolated. There are two isolation criteria, both of which are imposed: i) the sum of the scalar p T of charged hadron candidates from the interaction vertex that lie within a cone of radius R = ( η) 2 + ( φ) 2 = 0.3 around the photon candidate must be less than 5 GeV, where charged hadron candidates identified as conversion tracks associated with the photon candidate are excluded; ii) the sum of the scalar p T of additional neutral electromagnetic candidates within this same cone must be less than 2.5 GeV, after the contribution of additional interactions in the same bunch crossing has been removed.
The identification and trigger efficiencies are measured as functions of photon p T using data events containing a Z boson decaying to a μ + μ − pair in association with a photon, or to an e + e − pair where the electrons are treated as if they were photons [31]. The efficiency of the photon selection procedure in the kinematic range considered in the analysis is above 90% and 85% for candidates in the barrel and endcap regions, respectively. The ratio between the efficiencies measured in data and simulation is found to be lower than 1 by 3.5% for photons in the barrel region and by 6.5% for photons in the endcap region. No significant p T dependence of the efficiency ratios is observed, and a p T -independent correction is applied to the normalization of the simulated event samples to account for this difference.
The photon candidates in an event are grouped into all possible pairs. At least one photon candidate in the pair must have |η| < 1.44, i.e., be reconstructed in the barrel. Events with both photons in the endcaps are not considered, since their inclusion would increase the signal efficiency by only a few percent, at the cost of introducing a large background. Photon pairs are divided into two categories. The first category, denoted "EBEB", contains pairs for which both candidates lie in the barrel. For the second category, denoted "EBEE", one candidate lies in the barrel and the other in an endcap. The invariant mass m γ γ of the pair must satisfy m γ γ > 230 GeV for EBEB candidates and m γ γ > 330 GeV for EBEE candidates. The fraction of events in which more than one photon pair satisfies the selection criteria is approximately 1%. In these cases, only the pair with the largest scalar sum of photon p T is retained.
The selection efficiency times acceptance for signal events varies between 50% and 70%, depending on the signal hypothesis. Because of the different angular distribution of the decay products, the kinematic acceptance for the RS graviton resonances is lower than that of scalar resonances. For m X < 1 TeV the difference is approximately 20%. The two acceptances are similar for m X > 3 TeV.
The event selection procedure described above is the same as the one documented in [11]. It was finalized on the basis of studies with simulated signal and background event samples prior to inspection of the data in the search region of the diphoton invariant mass distribution, which is defined as m γ γ > 500 GeV.
A total of 6284 (2791) photon pairs are selected in the EBEB (EBEE) category. Of these, 461 (800) pairs have an invariant mass above 500 GeV. According to simulation, the direct production of two photons accounts, respectively, for 90% and 80% of the background events selected in the EBEB and EBEE categories. This prediction is tested in data using the method described in Ref. [44] and good agreement is found between data and simulation.
The diphoton invariant mass distribution of the selected events is shown in Fig. 1, for both the EBEB and EBEE categories. We perform an independent maximum likelihood fit to the data in each category using the function (1) This parametric form is chosen to model the background in the hypothesis tests discussed below. The results of the fits are shown in Fig. 1.

Likelihood fit
A simultaneous fit to the invariant mass spectra of events in the EBEB and EBEE event categories is performed to determine the compatibility of the data with the background-only and the signal+background hypotheses. The test statistic is based on the profile likelihood ratio: where S and B represent the probability density functions for resonant diphoton production and for the SM background, respectively. The parameter μ is the so-called signal strength, while θ represents the nuisance parameters of the model, used to account for systematic uncertainties. The θ notation indicates the best fit value of the parameter θ for any μ value, while θ μ denotes the best fit value of θ for a fixed value μ. To set upper limits on the rate of resonant diphoton production, the modified frequentist method known as CL s [45,46] is used, following the prescription described in Ref. [47]. The compatibility of the observation with the background-only hypothesis is evaluated by computing the background-only p-value. The latter is defined as the probability, in the background-only hypothesis, for q(0) to exceed the value observed in data. This quantity, the "local p-value" p 0 , does not take into account the fact that many signal hypotheses are tested.
Asymptotic formulas [48] are used in the calculations of exclusion limits and local p-values. The accuracy of the asymptotic approximation in the estimation of exclusion limits and significance is studied, using pseudo-experiments, for a subset of the hypothesis tests and is found to be about 10%.
The signal shape in m γ γ is determined from the convolution of the intrinsic shape of the resonance and the CMS detector response to photons. The intrinsic shape is taken from the pythia 8.2 generator. A grid of mass points with 125 GeV spacing, in the range 500-4500 GeV, is used. The resulting shapes are interpolated to intermediate points using a parametric description of the distribution. The detector response is determined using fully simulated signal samples of small intrinsic width, corrected through Gaussian smearing to agree with measurements based on Z → e + e − data. Nine uniformly spaced mass hypotheses in the range 500-4500 GeV are employed. The signal mass resolution, quantified through the ratio of the full width at half maximum of the distribution, divided by 2.35, to the peak position, is approximately 1.0% and 1.5% for the EBEB and EBEE categories, respectively. The signal normalization coefficients are proportional to the product of the kinematic acceptance and the signal efficiency within the acceptance region. These are computed, for each category, in simulated samples and interpolated to intermediate points using quadratic functions of m X and X /m X .
The background shape in m γ γ is described by the parametric function given by Eq. (1). The values of the parameters a and b are determined by the fit to data, with separate values for the EBEB and EBEE categories, and are treated as unconstrained nuisance parameters in the hypothesis tests.
The accuracy of the background parameterization is assessed using simulation and is quantified by studying the difference between the true and predicted numbers of background events in several m γ γ intervals in the search region. The relative widths of the intervals, defined by 2(x 1 − x 2 )/(x 1 + x 2 ) with x 1 and x 2 the lower and upper bin edges, range between 2% and 15%. Pseudoexperiments are drawn from the mass spectrum predicted by the simulation and are fit with the chosen background model. The total number of events in each pseudo-experiment is taken from a Poisson distribution whose mean is set equal to the observation in data. For each interval, the distribution of the pull variable, defined as the difference between the true and predicted numbers of events divided by the estimated statistical uncertainty, is constructed. If the absolute value of the median of this distribution is found to be above 0.5 in an interval, an additional uncertainty is assigned to the background parametrization. A modified pull distribution is then constructed, increasing the statistical uncertainty in the fit by an extra term, denoted the "bias term". The bias term is parametrized as a smooth function of m γ γ , which is tuned in such a manner that the absolute value of the median of the modified pull distribution is less than 0.5 in all intervals. The amplitude of the bias term function is comparable to that of the 1 standard deviation bands in Fig. 1. This additional uncertainty is included in the likelihood function by adding to the background model a component having the same shape as the signal. The normalization coefficient of this component is constrained to have a Gaussian distribution of mean zero, with a width equal to the integral of the bias term function over the full width at half maximum of the tested signal shape. The inclusion of this additional component has the effect of avoiding falsely positive or falsely negative tests that could be induced by a mismodeling of the background shape, and it reduces the sensitivity of the analysis by at most 10%.

Systematic uncertainties
The impact of systematic uncertainties in this analysis is smaller than that of the statistical uncertainties. The parametric background model has no associated systematic uncertainties except for the bias term uncertainty described in the previous section. Since the background shape coefficients a and b [Eq. (1)] are treated as unconstrained nuisance parameters, the associated uncertainties are statistical.
The systematic uncertainties in the signal normalization associated with the integrated luminosity, the selection efficiency, and the PDFs are 6.2%, 6.0%, and 6.0%, respectively. The uncertainty in . The uncertainty associated with the PDFs is evaluated by comparing the overall selection efficiency obtained with the CT10 [39], MSTW08 [50], and NNPDF2.3 [35] PDF sets and taking the largest deviation over all tested signal hypotheses. A 1% uncertainty is associated with the level of knowledge of the energy scale and accounts for the uncertainty in the energy scale at the Z boson peak and its extrapolation to higher masses. A 10% uncertainty is assigned to the knowledge of the photon energy resolution, corresponding to the uncertainty in the estimated additional Gaussian smearing determined at the Z boson peak.

Results for the 2016 data
The observed and expected 95% confidence level (CL) upper limits on the product of the production cross section (σ 13 TeV X ) and branching fraction to two photons (B γ γ ) for scalar and RS graviton resonances are shown in Fig. 2. Using the LO cross sections from pythia 8.2, RS gravitons with masses below 1.75, 3.75, and 4.35 TeV are excluded for k = 0.01, 0.1, and 0.2, respectively, corresponding to X /m X = 1.4 × 10 −4 , 1.4 × 10 −2 , and 5.6 × 10 −2 .
The value of p 0 for different signal hypotheses is shown in Fig. 3. The largest excess is observed for m X ≈ 620 GeV, and has a local significance of approximately 2.4 and 2.7 standard deviations for narrow spin-0 and RS graviton signal hypotheses, respectively. After taking into account the effect of searching for several signal hypotheses, i.e., searching over a range of widths and masses, the significance of the excess is reduced to less than one standard deviation. No excess is observed in the proximity of m X = 750 GeV.

Combination with the 2012 and 2015 data
The results obtained for the 2016 data are combined statistically with those obtained for the data discussed in Ref. Exclusion limits are set on the 13 TeV production cross section for both models, and background-only p-values are computed for the signal hypotheses.
The correlation model between the systematic uncertainties associated with 8 and 13 TeV data is described in Ref. [11]. It assumes all uncertainties to be uncorrelated except for those related to the knowledge of the PDFs, which are taken to be fully correlated, and those related to the knowledge of the photon energy scale, which are taken to have a linear correlation of 0.5. Taking the value of the linear correlation to be 0 or 1 would lead to negligible changes in the results. For the combination of the two 13 TeV data sets, the background shape and the associated bias term uncertainties are assumed to be fully correlated between the corresponding categories of the 2015 (3.8 T) and 2016 data. Independent background normalization coefficients are used for the two data sets. The uncertainty in the signal selection efficiency is taken to be uncorrelated between the 2015 and 2016 data, to account for the large statistical contribution and for the effect on the systematic contribution arising from changes in the data taking conditions, particularly in the instantaneous luminosity. The uncertainty in the knowledge of the integrated luminosity is treated as follows: a 2.3% uncertainty, corresponding to the knowledge of the absolute luminosity scale calibration determined with beam scans, is taken to be fully correlated between the 2015 (3.8 T) and 2016 data, and additional uncertainties of 1.5% and 5.8%, corresponding to the uncertainty in extrapolating the scale calibration to the data collection conditions, are applied, again respectively. Finally, the photon energy scale uncertainties are taken to be fully correlated between the two data sets, being dominated by the extrapolation to high energy. Fig. 4 shows the observed and expected 95% CL upper limits on the 13 TeV production cross section of the different signal hypotheses obtained with the combined analysis of the 13 TeV data recorded in 2015 and 2016. The upper limits on the production of scalar resonances decaying to two photons range from about 10 to 0.2 fb, for resonance masses between 0.5 and 4.5 TeV. Compared to the 2016 data alone, the sensitivity is improved by approximately  Fig. 5. The largest excess is observed for m X ≈ 1.3 TeV and has a local significance of about 2.2 standard deviations, corresponding to less than 1 standard deviation after accounting for the effect of searching for several signal hypotheses. For m X = 750 GeV, the 2.9 standard deviation local significance excess observed in the 2015 data is reduced to 0.8 standard deviations.
The observed and expected 95% CL upper limits on the 13 TeV signal production cross sections obtained through a combined analysis of the 8 TeV data from 2012 and the 13 TeV data from 2015 and 2016 are shown in Fig. 6. Compared to the combined 13 TeV data, the analysis sensitivity improves by about 10% at the low end of the m X range, while the improvement is negligible at the higher end of the range. Thus the lower limits on the mass of RS gravitons obtained by combining the 8 and 13 TeV data coincide with those obtained with the 13 TeV data alone.
The observed p 0 for X /m X = 1.4 × 10 −4 and 5.6 × 10 −2 obtained with the combined 8 and 13 TeV analysis is shown in Fig. 7.
The largest excess, observed for m X ≈ 0.9 TeV, has a local significance of about 2.2 standard deviations, corresponding to less than 1 standard deviation overall. For m X = 750 GeV, the excess with 3.4 standard deviation local significance [11] is reduced to about 1.9 standard deviations.

Summary
A search for the resonant production of high-mass photon pairs has been presented. The analysis is based on a sample of protonproton collisions collected by the CMS experiment in 2016 at √ s =  Century Project Advancement Project (Thailand); and the Welch Foundation, contract C-1845.