Search for supersymmetry in final states with photons and missing transverse momentum in proton-proton collisions at 13 TeV

Results are reported for a search for supersymmetry in final states with photons and missing transverse momentum in proton-proton collisions at the LHC. The data sample corresponds to an integrated luminosity of 35.9 fb$^{-1}$ collected at a center-of-mass energy of 13 TeV using the CMS detector. The results are interpreted in the context of models of gauge-mediated supersymmetry breaking. Production cross section limits are set on gluino and squark pair production in this framework. Gluino masses below 1.86 TeV and squark masses below 1.59 TeV are excluded at 95% confidence level.


Introduction
Supersymmetry (SUSY) and the Minimal Supersymmetric Standard Model [1][2][3][4][5][6] are extensions of the standard model (SM) that provide explanations for several outstanding issues with the SM. In particular, SUSY addresses the large quantum corrections to the mass term in the Higgs potential [7] and provides a viable dark matter candidate [8,9]. Models with general gaugemediated (GGM) SUSY breaking [10][11][12][13][14][15][16][17] have the additional benefit of naturally suppressing flavor violations in the SUSY sector. GGM models can have a wide range of features but typically result in final states that include the gravitino ( G) as the lightest supersymmetric particle (LSP). The next-to-lightest supersymmetric particle (NLSP) in these models is often taken to be a neutralino ( χ 0 1 ), which is a mixture of the bino, neutral wino, and neutral higgsinos. The conservation of R parity [18] implies that the gravitino is stable and remains undetected. Therefore, proton-proton (p p) collisions that produce SUSY particles will have an imbalance in the total observed transverse momentum, referred to as missing transverse momentum p miss T and defined as the negative vector sum of the transverse momenta of all visible particles in an event. Its magnitude is referred to as p miss T . If the composition of the neutralino NLSP is primarily bino-like, its main decay will be to a gravitino and a photon (γ), resulting in final states with significant missing transverse momentum and one or more photons. This paper presents a search for GGM SUSY in final states involving two photons and missing transverse momentum. The data sample, corresponding to an integrated luminosity of 35.9 fb −1 of pp collisions at a center-of-mass energy √ s = 13 TeV, was collected with the CMS detector in 2016. The analysis described here achieves a substantial improvement in sensitivity compared to the search performed by the CMS Collaboration on the smaller 2015 data set [19] and is comparable in sensitivity to similar searches from the ATLAS Collaboration [20,21].
Two simplified model frameworks [22][23][24][25][26] are used for the interpretation of the results. The T5gg model assumes gluino ( g) pair production and the T6gg model assumes squark ( q) pair production. The models assume a 100% branching fraction for the gluinos and squarks to decay as shown in Fig. 1. The squarks in the T6gg model can be either first or second generation. We assume a 100% branching fraction for the NLSP neutralino to decay to a nearly massless gravitino and a photon, χ 0 1 → G γ, resulting in characteristic events with large p miss T and two photons. In order for the analysis to be as model independent as possible, we choose not to define the signal region using any hadronic variables such as jet multiplicity or the scalar sum of the transverse momentum of the jets.
Standard model processes such as direct diphoton production or events with jets produced through the strong interaction, referred to as quantum chromodynamics (QCD) multijet events, can result in events with two photons. If the hadronic activity in the event is poorly measured, these processes can mimic the signal topology even though they lack genuine p miss T . For the case of QCD multijet events, there may be real photons in the event, or jets rich in electromagnetic (EM) energy that are misreconstructed as photons. Events with genuine p miss T also contribute to the composition of the candidate sample. These events are mainly from Wγ and W +jet(s) production, where an electron is misidentified as a photon in W → eν decays. A smaller background arises from Zγγ events where the Z boson decays to two neutrinos, Z → νν .

Detector, data, and simulated samples
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 covering the pseudorapidity region |η| < 2.5, as well as a lead tungstate crystal In gluino ( g) pair production in the T5gg simplified model (left), the gluino decays to a quarkantiquark pair (q q) and a neutralino ( χ 0 1 ). In squark ( q) pair production in the T6gg simplified model (right), the squark decays to a quark and a neutralino. In both cases, the neutralino subsequently decays to a photon (γ) and a gravitino ( G). electromagnetic calorimeter (ECAL), and a brass and scintillator hadron calorimeter (HCAL), each composed of a barrel and two endcap regions and covering the range |η| < 3.0. Forward calorimeters extend the coverage up to |η| < 5.0. Muons are measured in gas-ionization detectors embedded in the steel flux-return yoke outside the solenoid and cover the range |η| < 2.4. 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. [27].
Events of interest are selected using a two-tiered trigger system [28]. The first level is composed of custom hardware processors and uses information from the calorimeters and muon detectors to select events at a rate of around 100 kHz. The second level, known as the high-level trigger, consists of a farm of processors running a version of the full event reconstruction software optimized for fast processing. This trigger reduces the event rate to around 1 kHz before data storage. This analysis used a diphoton trigger to collect the data. The trigger requires a leading (subleading) photon with transverse momentum p T > 30 (18) GeV, and a combined invariant mass m γγ > 95 GeV. The photons are also required to pass isolation and cluster shape requirements.
Monte Carlo (MC) simulations are used for several purposes in this analysis. Simulations of the signal processes are used to determine signal efficiencies; background process simulation is used for validation of the analysis performance and to model the contribution from Zγγ → νν γγ events. The event generator MADGRAPH5 aMC@NLO 2.3.3 [29] is used to simulate the signal samples at leading order. The background samples are generated at next-toleading order using MADGRAPH5 aMC@NLO 2.4.2. For both signal and background processes, the parton showering, hadronization, SUSY particle decays, multiple-parton interactions, and the underlying event are described by the PYTHIA 8.212 [30] program with the CUETP8M1 [31] generator tune. The signal samples are generated with either two gluinos or two squarks and up to two additional partons in the matrix element calculation. The parton distribution functions (PDFs) are obtained from the NNPDF3.0 [32] set. For the background processes, the detector response is simulated using GEANT4 [33], while the CMS fast simulation [34,35] is used for the signal events. For both signal and background simulated events, additional pp interactions (pileup) are generated with PYTHIA and superimposed on the primary collision process. The simulated events are reweighted to match the pileup distribution observed in data.
The signal events were generated using the T5gg and T6gg simplified models and are characterized by the masses of the particles in the decay chain. For the gluino (squark) mass we simu-late a range of values from 1.4 to 2.5 (1.2 to 2.0) TeV in steps of 50 GeV. These mass ranges were selected to overlap and expand upon the mass ranges excluded by previous searches [19,20]. The neutralino masses range from 10 GeV up to the mass of the gluino or squark. The cross sections are calculated at next-to-leading-order (NLO) accuracy including the resummation of soft gluon emission at next-to-leading-logarithmic (NLL) accuracy [36][37][38][39][40], with all the unconsidered sparticles assumed to be heavy and decoupled. The uncertainties in the cross sections are calculated as described in Ref. [41].

Event selection
Photon, electron, muon, charged and neutral hadron candidates are reconstructed with the particle-flow event algorithm [42], which reconstructs particles based on information from all detector subsystems. The energy of photons is directly obtained from the ECAL measurement. The energy of electrons is determined from a combination of the electron momentum at the primary interaction vertex as determined by the tracker, the energy of the corresponding ECAL cluster, and the energy sum of all bremsstrahlung photons spatially compatible with originating from the electron track. The energy of muons is obtained from the curvature of the corresponding track. The energy of charged hadrons is determined from a combination of their momentum 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. Finally, the energy of neutral hadrons is obtained from the corresponding corrected ECAL and HCAL energy.
Photon candidates are required to satisfy a series of identification criteria to ensure a high purity [43]. The shape of the energy deposit in the ECAL must be consistent with that of an EM shower, and the amount of energy present in the corresponding region of the HCAL must not exceed 5% of the ECAL energy, since EM showers are expected to be contained almost entirely within the ECAL. To ensure high trigger efficiency, we require all photons to satisfy p T > 40 GeV. Because the SUSY signal models used in this analysis produce photons primarily in the central region of the detector and because the magnitude of the background increases considerably at high |η|, we consider only photons within the barrel fiducial region of the detector (|η| < 1.44).
To suppress quark and gluon jets that mimic photons, photon candidates are required to be isolated from other reconstructed particles. Separate requirements are made on the scalar p T sums of charged and neutral hadrons and EM objects in a cone of radius ∆R ≡ √ (∆η) 2 + (∆φ) 2 ≡ 0.3 around the photon candidate. Each p T sum is corrected for the effect of pileup, and in each case the momentum of the photon candidate itself is excluded. We further require that the photon candidate has no pixel detector track seed, to distinguish the candidate from an electron.
For the purpose of defining the various control regions used in the analysis, we apply an additional set of selection criteria. A misidentified "fake" photon ( f ) is defined as a photon candidate that satisfies looser requirements on photon isolation and neutral-hadron isolation and fails either the shape requirement for the ECAL clusters or the charged-hadron isolation requirement. In order to ensure that misidentified photons do not differ too much from our photon selection, upper limits are applied to both the charged-hadron isolation and cluster shape requirements. Importantly, because of the large amount of hadronic activity expected in our SUSY signal events, it is possible that real photons from the decay of a neutralino could fail the charged-hadron isolation requirement and therefore fall into the misidentified photon category. In order to avoid this potential signal contamination from SUSY events in the control regions, we additionally require that misidentified photons satisfy R 9 < 0.9, where R 9 is defined as the ratio of the energy deposited in a 3×3 array of ECAL crystals to the total energy in the cluster [43]. Real photons have values of R 9 close to unity, so by requiring R 9 < 0.9 we ensure that real photons from a possible SUSY signal will not enter our control regions.
Because of the similarity of the ECAL response to electrons and photons, Z → ee events are used to measure the photon identification efficiency. The selection of electron candidates is identical to that of photons, with the exception that the candidate is required to be matched to a pixel detector seed consistent with a track, to ensure that the electron selection is orthogonal to that of photons. The photon efficiency is measured via the tag-and-probe method [43]. The ratio of the observed to simulated efficiency is found to be consistent with unity and independent of p T and η. The efficiency of the pixel detector seed veto for photons is measured in Z → µµγ events and is found to agree between data and simulation.
Events are then assigned to one of four mutually exclusive categories depending on the selection of their highest p T EM objects: γγ, ee, f f , and eγ. The two EM objects are required to be separated by ∆R > 0.6. Finally, because of the trigger requirements described in Section 2, the invariant mass of the two EM objects is required to be greater than 105 GeV.
In addition to the requirements already described, any event with a muon satisfying p T > 25 GeV and |η| < 2.4 as well as track quality and isolation requirements is vetoed. Similarly, we veto events with any additional electrons satisfying p T > 25 GeV, |η| < 2.5, and signal shape and isolation requirements.
Events in the candidate γγ sample are divided into the low-p miss T control region (p miss T < 100 GeV) and the high-p miss T signal region (p miss T > 100 GeV). The signal region is further divided into six p miss T bins that were chosen such that there is a sufficient number of events from the f f control sample in each bin.

Estimation of backgrounds
As described in Section 1, there are three primary backgrounds to this analysis. QCD processes such as multijet production can emulate the signal topology if the hadronic activity in the event is mismeasured. A second background arises from electroweak (EWK) processes that have genuine p miss T from the production of neutrinos. There is also a small contribution from Zγγ → γγνν events.
The contribution from the QCD background is estimated from the observed data using the f f control sample. The ratio of the event yield in the candidate γγ sample to that in the f f sample is constructed as a function of p miss T . More f f events are observed at high p miss T relative to the γγ sample. Different functional forms were investigated to model the p miss T dependence, and an exponential function was found to describe the data the best. We fit the γγ to f f ratio in the p miss T < 100 GeV control region. The predicted number of QCD background events (N i QCD ) in bin i of the signal region is then given by the following equation, where N i f f is the number of observed f f events and g i ave is the average value of the fit function g(p miss T ) in that bin: In order to set a systematic uncertainty on the method, we derive a second QCD background prediction by noting that the p miss T distribution of the f f control sample is dependent on the R 9 requirement on the misidentified photons. An alternate f f control sample is built using photon candidates that satisfy all of the requirements for misidentified photons as outlined in Section 3, with the exception that the R 9 requirement is reversed. In the p miss T < 100 GeV control region, we perform an exponential fit to the ratio of the event yield in the high-R 9 f f sample to that of the nominal, low-R 9 f f sample. This function (h(p miss T )) represents the correction required to account for the effect of the R 9 selection on the p miss T distribution. The size of the correction is between 20 and 40% in the p miss T > 100 GeV signal region. Multiplying the number of low-R 9 f f events observed in the signal region by this function gives a proxy high-R 9 f f sample.
For p miss T < 100 GeV, the ratio of the p miss T distribution in the γγ sample to that of the proxy f f sample is fit to a constant C. We multiply this constant value by the proxy f f yield in the signal region to get a second prediction for the QCD background in bin i.
The two background estimation methods give values that are consistent within the uncertainties. All three of the fits used in the two methods are found to represent the data well in the low-p miss T control region. Several studies were performed to verify the procedure, including using a mixed-R 9 f f sample with one misidentified photon satisfying R 9 > 0.9 and one satisfying R 9 < 0.9 to confirm that the exponential fit continues to accurately describe the mixed-R 9 f f to nominal f f ratio in the high-p miss T signal region. As an additional check, a control sample with one photon and one misidentified photon was used as a proxy for the γγ candidate sample in a closure test of the method up to p miss T = 250 GeV. At larger values of p miss T , there is the potential for signal contamination in the γ f control sample.
Another background for this analysis comes from EWK processes with genuine p miss T . This background primarily involves Wγ and W+jets events where the W decays to an electron and a neutrino and the electron is misidentified as photon. This leads to final states with photons and significant p miss T . To obtain an estimate of the EWK background in the signal region, the mass peaks from the Z boson in the ee control sample and the eγ control sample are modeled using an extended likelihood fit for the signal plus background hypothesis. The rate at which electrons are misidentified as photons ( f e→γ ) is calculated using the signal fit integrals N eγ and N ee for each sample. These can be expressed in terms of the number of true Z bosons, N True Z : The factor of 2 in N eγ occurs because either electron in the event could be misidentified as a photon.
Taking the ratio of these two values, we find that the misidentification rate is given by f e→γ = N eγ /(2N ee + N eγ ). The misidentification rate is calculated as a function of several kinematic variables, including the vertex multiplicity and p miss T of the event and the p T of the EM objects. A 30% uncertainty is applied to cover the observed dependencies. The final EWK background prediction is given by scaling the number of events in the eγ control sample by the factor The irreducible Zγγ background is modeled via simulation. A 50% uncertainty is applied to conservatively cover the effects from the statistical uncertainty of the MC sample, the PDF uncertainty in the cross section, NNLO corrections in the simulation, and any other sources of potential mismodeling.

Sources of systematic uncertainty
Systematic uncertainties are calculated for each contribution to the total background prediction. In addition, systematic uncertainties are assigned for the signal efficiency and the integrated luminosity. The value of each uncertainty and the method used to calculate it are described below.
The largest uncertainties in the background prediction come from uncertainties associated with the QCD background estimate. The magnitude of each uncertainty is shown in Table 1 for the six signal bins. The statistical uncertainty from the f f control sample ranges from 7 to 79% in the signal region. The uncertainty obtained from propagating the errors in the fit parameters to the final prediction is between 2 and 5%. Finally, as described in Section 4, a systematic uncertainty in the fitting procedure is calculated by comparing the primary prediction to the cross check prediction derived using the high-R 9 f f sample. The systematic uncertainty is taken as the difference between the two methods or the uncertainty in that difference, whichever is larger, and ranges between 10 and 83% in the signal region. Uncertainties in the EWK background prediction include the statistical uncertainty from the eγ control sample and the 30% uncertainty in the rate at which electrons are misidentified as photons. The statistical uncertainty is less than 9% in each of the six signal bins.
There are also several uncertainties associated with the signal efficiency. The statistical uncertainty from the size of the T5gg or T6gg signal scans ranges from 2 to 44% depending on the mass bin. The PDF uncertainties in the cross sections for signal simulation are between 19 and 35% and are taken from Ref. [41]. Other uncertainties include how well the jet energy scale is known (1 to 30%) and the uncertainty in the photon identification efficiency (2.5%). The uncertainty in the integrated luminosity of the data sample is 2.5% [44].

Results
We determine 95% confidence level (CL) upper limits on gluino pair production and squark pair production cross sections using the modified frequentist CL s method [45,46]. The test statistic is an LHC-style profile likelihood ratio [47], and its distribution is determined using the asymptotic approximation [48]. The likelihood function is constructed from the background and signal p miss T distributions across the six bins described in Section 4. The systematic uncertainties described in Section 5 are included in the test statistic as constrained nuisance parameters. Systematic uncertainties which directly affect the yields of processes are assumed to follow a log-normal probability distribution, while statistical uncertainties from the limited size of the control samples and the signal MC samples are modeled using gamma probability distributions.
The full background prediction and the measured p miss T distribution prior to the fit are shown in Fig. 2. The expected and observed numbers of events for each bin in the signal region are shown in Table 2 for the pre-fit distributions and Table 3 for the post-fit distributions. Notably, in the last bin we observe 12 events and expect 5.4 +1.6 −1.5 background events (pre-fit). The significance of the observed data after the fit across all six bins of the signal region is calculated using the likelihood ratio test for each mass pair value of m versus m q for the T5gg and T6gg models, respectively. The significance does not strongly depend on the SUSY masses, and for all masses in both models, the significance is found to correspond to between 2.35 and 2.45 standard deviations. Several studies were performed to characterize the fit and the excess in the final p miss T bin and to ensure that the statistical treatment of the data is robust. In particular, the pre-and postfit distributions were checked to make sure that the pulls are consistent within the uncertainties. In Fig. 3 we present 95% CL upper limits on the gluino and squark pair production cross sections as a function of the mass pair values for the two models considered in this analysis. From the NLO+NLL predicted signal cross sections and their uncertainties we derive contours representing lower limits in the SUSY mass plane. We also show expected limit contours based on the expected experimental cross section limits and their uncertainties. For values of the neutralino mass between 500 and 1500 GeV, we expect to exclude gluino masses up to 2.02 TeV and squark masses up to 1.74 TeV. This is an improvement of approximately 400 and 300 GeV, respectively, upon the reach of the previous CMS result [19]. We observe exclusions for gluino masses up to 1.86 TeV and squark masses up to 1.59 TeV. The observed exclusions are lower than the expected exclusions because of the observed excess in the data.

Summary
The results of a search for general gauge-mediated supersymmetry breaking in proton-proton collisions with two photons and missing transverse momentum in the final state are reported. The analysis was performed using data corresponding to 35.9 fb −1 of integrated luminosity, recorded with the CMS detector in 2016 at a proton-proton center-of-mass energy of 13 TeV. An excess of events corresponding to 2.4 standard deviations is observed. Limits are determined on the masses of supersymmetric particles in two simplified models using data-driven background estimation methods and NLO+NLL signal cross section calculations.
In both models, the next-to-lightest supersymmetric particle is the neutralino, which decays with a 100% branching fraction to a photon and a gravitino, the lightest supersymmetric particle. The first simplified model assumes gluino pair production, with each gluino decaying to a neutralino and quarks. The second simplified model assumes squark pair production, with each squark decaying to a quark and a neutralino. The expected limits on gluino and squark masses, for the respective models, are 2.02 and 1.74 TeV at 95% confidence level. This is an increase in sensitivity of more than 300 GeV for each model with respect to the analysis performed with 2.3 fb −1 of integrated luminosity collected using the CMS detector in 2015. The observed exclusions are for gluino masses less than 1.86 TeV and squark masses less than 1.59 TeV, where the difference between the expected and observed exclusions is driven by the excess observed in the data. The analysis described in this paper improves the observed limits by 210 GeV for   Figure 3: The 95% confidence level upper limits on the gluino (left) and squark (right) pair production cross sections as a function of gluino or squark and neutralino masses. The contours show the observed and expected exclusions assuming the NLO+NLL cross sections, with their one standard deviation uncertainties. gluino masses and 220 GeV for squark masses with respect to the previous CMS result.

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: