Measurements of Higgs boson production cross sections and couplings in the diphoton decay channel at $\sqrt{s} = $ 13 TeV

Measurements of Higgs boson production cross sections and couplings in events where the Higgs boson decays into a pair of photons are reported. Events are selected from a sample of proton-proton collisions at $\sqrt{s} =$ 13 TeV collected by the CMS detector at the LHC from 2016 to 2018, corresponding to an integrated luminosity of 137 fb$^{-1}$. Analysis categories enriched in Higgs boson events produced via gluon fusion, vector boson fusion, vector boson associated production, and production associated with top quarks are constructed. The total Higgs boson signal strength, relative to the standard model (SM) prediction, is measured to be 1.12 $\pm$ 0.09. Other properties of the Higgs boson are measured, including SM signal strength modifiers, production cross sections, and its couplings to other particles. These include the most precise measurements of gluon fusion and vector boson fusion Higgs boson production in several different kinematic regions, the first measurement of Higgs boson production in association with a top quark pair in five regions of the Higgs boson transverse momentum, and an upper limit on the rate of Higgs boson production in association with a single top quark. All results are found to be in agreement with the SM expectations.


Introduction
Since the discovery of a Higgs boson (H) by the ATLAS and CMS Collaborations in 2012 [1][2][3], an extensive programme of measurements focused on characterising its properties and testing its compatibility with the standard model (SM) of particle physics has been performed. Analysis of data collected during the second run of the CERN LHC at √ s = 13 TeV has already resulted in the observation of Higgs boson production mechanisms and decay modes predicted by the SM [4][5][6][7]. The most precise measurements are obtained by combining results from different Higgs boson decay channels. Such combinations have enabled the total Higgs boson production cross section to be measured with an uncertainty of less than 7% [8,9]. All reported results have so far been consistent with the corresponding SM predictions.
In the SM, the H → γγ decay has a small branching fraction of approximately 0.23% for a Higgs boson mass (m H ) around 125 GeV [10]. However, its clean final-state topology with two well-reconstructed photons provides a narrow invariant mass (m γγ ) peak that can be used to effectively distinguish it from background processes. As a result, H → γγ is one of the most important channels for precision measurements of Higgs boson properties. Furthermore, it is one of the few decay channels that is sensitive to all principal Higgs boson production modes.
The results reported in this paper build upon previous analyses performed by the CMS Collaboration [11,12]. Here, the data collected by the CMS experiment between 2016 and 2018 are analysed together. The resulting statistical power of the combined data set improves the precision on existing measurements and allows new measurements to be made. The structure of this analysis is designed to enable measurements within the simplified template cross section (STXS) framework [10]. Using this structure, various measurements of Higgs boson properties can be performed. These include SM signal strength modifiers, production cross sections, and the Higgs boson's couplings to other particles. Measurements of all these quantities are reported in this paper.
The STXS framework provides a coherent approach with which to perform precision Higgs boson measurements. Its goal is to minimise the theory dependence of Higgs boson measurements, lessening the direct impact of SM predictions on the results, and to provide access to kinematic regions likely to be affected by physics beyond the SM (BSM). At the same time, this approach permits the use of advanced analysis techniques to optimise sensitivity. Reducing theory-dependence is desirable because it makes the measurements both easier to reinterpret and means they are less affected by potential updates to theoretical predictions, making them useful over a longer period of time [13]. The results presented within the STXS framework nonetheless depend on the SM simulation used to model the experimental acceptance of the signal processes, which could be modified in BSM scenarios.
The strategy employed in this analysis is to construct analysis categories enriched in events from as many different kinematic regions as possible, thereby providing sensitivity to the individual regions defined in the STXS framework. This permits measurements to be performed across all the major Higgs boson production modes, including gluon fusion (ggH), vector boson fusion (VBF), vector boson associated production (VH), production associated with a top quark-antiquark pair (ttH), and production in association with a single top quark (tH).
In addition to measurements within the STXS framework, this paper contains several other interpretations of the data. The event categorisation designed to target the individual STXS regions also provides sensitivity to signal strength modifiers, both for inclusive Higgs boson production and for individual production modes, as well as measurements within the κ-framework [14].
The paper is structured as follows. The CMS detector is described in Section 2. An overview of the STXS framework is given in Section 3, together with a summary of the overall strategy of this analysis. In Section 4, details of the data and simulation used to design and perform the analysis are given. The reconstruction of candidate H → γγ events is described in Section 5, before the event categorisation procedure is explained in Section 6. The techniques used to model the signal and background are outlined in Section 7, with the associated systematic uncertainties listed in Section 8. The results are presented in Section 9, with tabulated versions provided in HEPDATA [15]. Finally, the paper is summarised in Section 10.

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 ECAL consists of 75 848 lead tungstate crystals, which provide coverage in pseudorapidity |η| < 1.48 in the barrel region and 1.48 < |η| < 3.00 in the two endcap regions. Preshower detectors consisting of two planes of silicon sensors interleaved with a total of 3 radiation lengths of lead are located in front of each EE detector. Forward calorimeters extend the η coverage provided by the barrel and endcap detectors. Muons are detected in gas-ionisation chambers embedded in the steel flux-return yoke outside the solenoid.
Events of interest are selected using a two-tiered trigger system [16]. The first level, composed of custom hardware processors, uses information from the calorimeters and muon detectors to select events at a rate of around 100 kHz within a fixed time interval of less than 4 µs. The second level, known as the high-level trigger, consists of a farm of processors running a version of the full event reconstruction software optimised for fast processing, and reduces the event rate to around 1 kHz before data storage [17].
The particle-flow (PF) algorithm [18] aims to reconstruct and identify each individual particle (PF candidate) in an event, with an optimised combination of information from the various elements of the CMS detector. The energy of photons is 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 energies.
For each event, hadronic jets are clustered from these reconstructed particles using the infrared and collinear safe anti-k T algorithm [19,20] with a distance parameter of 0.4. Jet momentum is determined as the vectorial sum of all particle momenta in the jet, and is found from simulation to be, on average, within 5 to 10% of the true momentum over the whole transverse momentum (p T ) spectrum and detector acceptance. Additional proton-proton interactions within the same or nearby bunch crossings (pileup) can contribute additional tracks and calorimetric energy depositions to the jet momentum. To mitigate this effect, charged particles identified to be originating from pileup vertices are discarded and an offset correction is applied to correct for remaining contributions. Jet energy corrections are derived from simulation to bring the measured response of jets to that of particle level jets on average. In situ measurements of the momentum balance in dijet, photon + jet, Z + jet, and multijet events are used to account for any residual differences in the jet energy scale between data and simulation [21]. The jet energy resolution amounts typically to 15-20% at 30 GeV, 10% at 100 GeV, and 5% at 1 TeV [21]. Additional selection criteria are applied to each jet to remove jets potentially dominated by anomalous contributions from various subdetector components or reconstruction failures.
The missing transverse momentum vector p miss T is computed as the negative vector p T sum of all the PF candidates in an event, and its magnitude is denoted as p miss T [22]. The p miss T is modified to account for corrections to the energy scale of the reconstructed jets in the event.
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. [23].

Analysis strategy 3.1 The STXS framework
In the STXS framework, kinematic regions based upon Higgs boson production modes are defined. These regions, or bins, exist in varying degrees of granularity, following sequential "stages". At the so-called STXS stage 0, the bins correspond closely to the different Higgs boson production mechanisms. Events where the absolute value of the Higgs boson rapidity, |y H |, is greater than 2.5 are not included in the definition of the bins because they are typically outside of the experimental acceptance. Measurements of stage-0 cross sections in the H → γγ decay channel were presented by the CMS Collaboration in Ref. [12]. Additionally, an analysis probing the coupling between the top quark and Higgs boson in the diphoton decay channel was recently performed by the CMS Collaboration [24]. Several other stage-0 measurements in different decay channels have also been made by both the ATLAS and CMS Collaborations [25][26][27][28][29][30][31]. Each experiment has also presented results combining the various analyses [8,9].
At STXS stage 1, a further splitting of the bins using the events' kinematic properties is performed [32]. This provides additional information for different theoretical interpretations of the measurements, and enhances the sensitivity to possible signatures of BSM physics. Furthermore, increasing the number of independent bins reduces the theory-dependence of the measurements; within each bin SM kinematic properties are assumed, and thus splitting the bins allows these assumptions to be partially lifted.
Measurements at stage 1 of the framework have already been reported by the ATLAS Collaboration [25,26,33]. Following these results, adjustments to the framework and its definitions were made, such that the most recent definition of STXS bins is referred to as STXS stage 1.2. The first measurement of STXS stage-1.2 cross sections was recently performed by the CMS Collaboration [34].
The full set of STXS stage-1.2 bins is described below and an illustration is given in Fig. 1. The ggH region (blue) is split into STXS bins using the Higgs boson transverse momentum (p H T ), the number of jets, and additionally has a VBF-like region with high dijet mass (m jj ). This VBF-like region is split into four STXS bins according to m jj and the transverse momentum of the Higgs boson plus dijet system (p Hjj T ). Events originating from bbH production are grouped with the ggH production mode, as are those from gluon-initiated production in association with a vector boson (ggZH) where the vector boson decays hadronically. The VBF and hadronic VH modes are considered together as electroweak qqH production (orange). Here the STXS bins are defined using the number of jets, p H T , m jj and p Hjj T . The four STXS bins which define the qqH rest region are not explicitly probed in this analysis. The leptonic VH STXS bins (green) are split into three separate regions representing the WH, ZH, and ggZH production modes, which are further divided according to the number of jets and the transverse momentum of the vector boson (p V T ) that decays leptonically. The ttH production mode (pink) is split only by p H T . Finally, the tH STXS bin includes contributions from both the tHq and tHW production modes. All references to STXS bins hereafter imply the STXS stage-1.2 bins. Further details on the exact definitions are contained in Section 6, describing the event categorisation. All the production mechanisms shown in Fig. 1 are measured independently in this analysis.

Analysis categorisation
To perform measurements of Higgs boson properties, analysis categories must first be constructed where the narrow signal peak is distinguishable from the falling background m γγ spectrum. The categorisation procedure uses properties of the reconstructed diphoton system and any additional final-state particles to improve the sensitivity of the analysis. As part of the categorisation, dedicated selection criteria and classifiers are used to select events consistent with the tH, ttH, VH, VBF, and ggH production modes. This both increases the analysis sensitivity and enables measurements of individual production mode cross sections to be performed.
In order to measure cross sections of STXS bins individually, events deemed to be compatible with a given production mode are further divided into analysis categories that differentiate between the various STXS bins. For most production modes, the divisions are made using the detector-level equivalents of the particle-level quantities used to define the STXS bins; an example is using p γγ T to construct analysis categories targeting STXS bins defined by p H T values. Increasing the total number of analysis categories to target individual STXS bins in this way does not degrade the analysis' sensitivity to the individual production mode and total Higgs boson cross sections. For each production mode, the event categorisation is designed to target all of the STXS bins to which some sensitivity can be obtained in the diphoton decay channel with the available data.
Several different machine learning (ML) algorithms are used throughout this analysis for both regression and classification tasks. Examples include regressions that improve the agreement between simulation and data, and classification to improve the discrimination between signal and background processes. The usage of ML techniques for event categorisation is also found to improve the separation between different STXS bins, which further improves the sensitivity of STXS measurements. For the training of boosted decision trees (BDTs), either the XGBOOST [35] or the TMVA package [36] package is used. The TENSORFLOW [37] package is used to train deep neural networks (DNNs).
For the ggH phase space, almost all of the STXS bins can be measured individually, without any bin merging (blue in Fig. 1). The exceptions are the high dijet mass (m jj ) STXS bins, which are difficult to distinguish from VBF events. Furthermore, the sensitivity to STXS bins with particularly high p H T is limited. Analysis categories are constructed using a BDT to assign the most probable STXS bin for each event. The amount of background is reduced using another BDT, referred to as the diphoton BDT. The diphoton BDT is trained to discriminate between all Higgs boson signal events and all other modes of SM diphoton production. Throughout the analysis, events originating from the bbH production mode are grouped together with ggH events.
The VBF production mode and VH production where the vector boson decays hadronically are considered together as (EW) qqH production (orange in Fig. 1). A set of analysis categories enriched in VBF-like events, where a dijet with high m jj is present, is defined. These analysis categories make use of the same diphoton BDT used in the analysis categories targeting ggH to reduce the number of background events. Additionally, a BDT based on the kinematic properties of the characteristic VBF dijet system, known as the dijet BDT, is utilised. The dijet BDT is trained to distinguish between three different classes of events with a VBF-like topology: VBF events, ggH events, and events produced by all other SM processes. This enables VBF events to be effectively separated from both VBF-like ggH events and other SM backgrounds. At least one analysis category is defined to target each VBF-like qqH STXS bin. Additional analysis categories enriched in VH-like events, where the vector boson decays hadronically to give a dijet whose m jj is consistent with a W or Z boson, are defined. These make use of a dedicated VH hadronic BDT to reduce both the number of background events and contamination from ggH events.
Analysis categories targeting VH leptonic production (green in Fig. 1) are divided into three categorisation regions, containing either zero, one, or two reconstructed charged leptons (electrons or muons). Each categorisation region uses a dedicated BDT to reduce the background contamination. It is not possible to measure STXS bins individually with the available data set. Nonetheless, where a sufficient number of events exists, analysis categories are constructed to provide sensitivity to merged groups of STXS bins.
In this analysis, ttH and tH production cross sections are measured independently (ttH STXS bins are purple in Fig. 1, whilst tH is yellow). For this purpose, a dedicated DNN referred to as the top DNN is trained to discriminate between tH and ttH events. An analysis category enriched in tH events is defined that uses the top DNN to reduce the contamination from ttH events, with a BDT used to reject background events from other sources.
The analysis categories targeting ttH production are based on those described in Ref. [24], with separate channels for hadronic and leptonic top quark decays. In each channel, a dedicated BDT is trained to reject background events. Furthermore, the top DNN is used to reduce the amount of contamination from tH events. The analysis categories are divided to provide the sensitivity to the STXS bins, for which four p H T ranges are defined. It is possible for an event to pass the selection criteria for more than one analysis category. To unambiguously assign each event to only one analysis category, a priority sequence is defined. Events that could enter more than one analysis category are assigned to the analysis category with the highest priority. The priority sequence is based on the expected number of signal events, with a higher priority assigned to analysis categories with a lower expected signal yield. This ordering enables the construction of analysis categories containing sufficiently high fractions of the Higgs boson production mechanisms with lower SM cross sections, which is necessary to perform independent measurements of these processes.
Events in data and the corresponding simulation for all three years of data-taking from 2016 to 2018 are grouped together in the final analysis categories. This gives better performance than constructing analysis categories for each year individually, requiring fewer analysis categories in total for a comparable sensitivity. Separating the analysis categories by year would enable differences in the detector conditions -such as the variation in m γγ resolution -to be exploited. However this is found to be less important than the advantage of having a greater number of events with which to train multivariate classifiers and optimise the analysis category definitions. Furthermore, the variations in detector conditions are relatively modest, and in general not substantially greater than variations within a given year of data-taking, which allows all data collected in each of the three years to be analysed together.
Nonetheless, simulated events are generated for each year separately, with the corresponding detector conditions, before they are merged together. This accounts for the variation in the detector itself, in the event reconstruction procedure, and in the LHC beam parameters. Furthermore, corrections to the photon energy scale and other procedures relating to the event reconstruction are also performed for each year individually. Only when performing the final division of selected diphoton events into the analysis categories are the simulated and data events from different years processed together. The full description of all the analysis categories is given in Section 6.
Once the selection criteria for each analysis category are defined, results are obtained by performing a simultaneous fit to the resulting m γγ distributions in all analysis categories. The results of several different measurements with different observables are reported in Section 9. For measurements within the STXS framework, it is not possible to measure each STXS bin individually. Therefore for each fit, a set of observables is defined by merging some STXS bins. In this paper, the results of two scenarios with different parameterisations of the STXS bins are provided. In addition, measurements of SM signal strength modifiers are reported, both for inclusive Higgs boson production and per production mode. Finally, measurements of Higgs boson couplings within the κ-framework are also shown.

Data samples and simulated events
The analysis exploits proton-proton collision data at √ s = 13 TeV, collected in 2016, 2017, and 2018 and corresponding to integrated luminosities of 35.9, 41.5, and 59.4 fb −1 , respectively. The integrated luminosities of the 2016-2018 data-taking periods are individually known with uncertainties in the 2.3-2.5% range [38][39][40], while the total (2016-2018) integrated luminosity has an uncertainty of 1.8%, the improvement in precision reflecting the (uncorrelated) time evolution of some systematic effects. In this section, the data sets and simulated event samples for all three years are described. Any differences between the years are highlighted in the text.
Events are selected using a diphoton high-level trigger with asymmetric photon p T thresholds of 30 (30) and 18 (22) GeV in 2016 (2017 and 2018) data. A calorimetric selection is applied at trigger level, based on the shape of the electromagnetic shower, the isolation of the photon candidate, and the ratio of the hadronic and electromagnetic energy deposits of the shower. The R 9 variable is defined as the energy sum of the 3×3 crystals centred on the most energetic crystal in the candidate electromagnetic cluster divided by the energy of the candidate. The value of R 9 is used to identify photons undergoing a conversion in the material upstream of the ECAL. Unconverted photons typically have narrower transverse shower profiles, resulting in higher values of the R 9 variable, compared to converted photons. The trigger efficiency is measured from Z → ee events using the "tag-and-probe" technique [41]. The efficiency measured in data in bins of p T , R 9 , and η is used to weight the simulated events to replicate the trigger efficiency observed in data.
A Monte Carlo (MC) simulated signal sample for each Higgs boson production mechanism is generated using MADGRAPH5 aMC@NLO (version 2.4.2) at next-to-leading order accuracy [42] in perturbative quantum chromodynamics (QCD). For each production mode, events are generated with m H = 120, 125, and 130 GeV. Events produced via the gluon fusion mechanism are weighted as a function of p H T and the number of jets in the event, to match the prediction from the NNLOPS program [43]. All parton-level samples are interfaced with PYTHIA8 version 8.226 (8.230) [44] for parton showering and hadronization, with the CUETP8M1 [45] (CP5 [46]) tune used for the simulation of 2016 (2017 and 2018) data. Parton distribution functions (PDFs) are taken from the NNPDF 3.0 [47] (3.1 [48]) set, when simulating 2016 (2017 and 2018) data. The production cross sections and branching fractions recommended by the LHC Higgs Working Group [10] are used. The relative fraction of each STXS bin for each inclusive production mode at particle level is taken from simulation and used to compute the SM prediction for the production cross section in each STXS bin. Additional signal samples generated with POWHEG 2.0 [49][50][51][52][53][54] at next-to-leading order accuracy in perturbative QCD are used to train some of the multivariate discriminants described in Section 6.
The dominant source of background events in this analysis is due to SM diphoton production. A smaller component comes from γ +jet or jet+jet events, in which jets are misidentified as photons. In the final fits of the analysis, the background is estimated directly from the diphoton mass distribution in data. Simulated background events from different event generators are only used for the training of multivariate discriminants. The diphoton background is generated with the SHERPA (version 2.2.4) generator [55]. It includes the Born processes with up to 3 additional jets, as well as the box processes at leading order accuracy. The γ +jet and jet+jet backgrounds are simulated at leading order with PYTHIA8, after applying a filter at generator level to enrich the production of jets with a high electromagnetic activity. The filter requires a potential photon signal coming from photons, electrons, or neutral hadrons with p T > 15 GeV. In addition, the filter requires no more than two charged particles (p T > 1.6 GeV and |η| < 2.2) in a cone of radius R = √ (∆η) 2 + (∆ϕ) 2 < 0.2 (where ϕ is the azimuthal angle in radians) around the photon candidate, mimicking the tracker isolation described in Section 5.
A sample of Drell-Yan events is simulated with MADGRAPH5 aMC@NLO, and is used both to derive corrections for simulation and for validation purposes.
The response of the CMS detector is simulated using the GEANT4 package [56]. This includes the simulation of the multiple proton-proton interactions taking place in each bunch crossing. These can occur at the nominal bunch crossing (in-time pileup) or at the crossing of previous and subsequent bunches (out-of-time pileup), and the simulation accounts for both. Simulated out-of-time pileup is limited to a window of [−12, +3] bunch crossings around the nominal, in which the effects on the observables reconstructed in the detector are most relevant. Simulated events are weighted to reproduce the distribution of the number of interaction vertices in data. The average number of interactions per bunch crossing in data in the 2016 (2017 and 2018) data sets is 23 (32).

Photon reconstruction
Efficiently reconstructing photons with an accurate and precise energy determination plays a very important role in the sensitivity of this analysis. This section describes in detail the procedures used to reconstruct the photon energy and the photon preselection criteria.
Photon candidates are reconstructed from energy clusters in the ECAL not linked to any chargedparticle trajectories seeded in the pixel detector. The clusters are built around a "seed" crystal, identified as a local energy maximum above a given threshold. The clusters are grown with a so-called topological clustering, where crystals with at least one side in common with a crystal already in the cluster and with an energy above a given threshold are added to the existing cluster itself. Finally, the clusters are dynamically merged into "superclusters" to ensure good containment of the shower, accounting for geometrical variation along η, and optimising the robustness of the energy resolution against pileup. The energy of the photon is estimated by summing the energy of each crystal in the supercluster, calibrated and corrected for response variations in time [57]. The photon energy is corrected for the imperfect containment of the electromagnetic shower and the energy losses from converted photons. The correction is computed with a multivariate regression technique trained on simulated photons, which estimates simultaneously the energy of the photon and its uncertainty.
After the application of this simulation-derived correction, some differences remain between data and simulation. A sequence of additional corrections are applied to improve the agreement between the two, using Z → ee events where the electrons are reconstructed as photons. First, any residual drift in the energy scale in data over time is corrected for in bins corresponding approximately to the duration of one LHC fill. The second step involves modifying the energy scale in data and the energy resolution in simulation. A set of corrections is derived to align the mean of the dielectron mass spectrum in data with the expected value from simulation, and to smear the resolution in simulation to match that observed in data. These corrections are derived simultaneously in bins of |η| and R 9 . Further details on this procedure are contained in Ref. [58]. Figure 2 shows comparisons between data and simulation after all corrections are applied for two cases where both electrons are reconstructed in the ECAL barrel and endcaps, respectively. In both cases the dielectron invariant mass spectra for the data and simulation are compatible within the uncertainties.
Once the photon energy correction has been applied, photon candidates are preselected before being used to form diphoton candidates. Requirements are placed on the photons' kinematic, shower shape, and isolation variables at values at least as stringent as those applied in the trigger. The preselection criteria are as follows:    • minimum p T of the leading and subleading photons greater than 35 and 25 GeV, respectively; • pseudorapidity of the photons |η| < 2.5 and not in the barrel-endcap transition of 1.44 < |η| < 1.57; • preselection on the R 9 variable and on σ ηη -the lateral extension of the shower, defined as the energy-weighted spread within the 5×5 crystal matrix centred on the crystal with the largest energy deposit in the supercluster -to reject ECAL energy deposits incompatible with a single, isolated electromagnetic shower, such as those coming from neutral mesons; • preselection on the ratio of the energy in the HCAL tower behind the supercluster's seed cluster to the energy in the supercluster (H/E), in order to reject hadronic showers; • electron veto, which rejects the photon candidate if its supercluster in the ECAL is near to the extrapolated path of a track compatible with an electron. Tracks compatible with a reconstructed photon conversion vertex are not considered when applying this veto. • requirement on the photon isolation (I ph ), defined as the p T sum of the particles identified as photons inside a cone of size R = 0.3 around the photon direction; • requirement on the track isolation in a hollow cone (I tk ), the p T sum of all tracks in a cone of size R = 0.3 around the photon candidate direction, excluding tracks in an inner cone of size R = 0.04 to avoid counting tracks arising from photon conversion into electron-positron pairs; • loose requirement on charged-hadron isolation (I ch ), the p T sum of charged hadrons inside a cone of size R = 0.3 around the photon candidate.
The geometrical acceptance requirement is applied to the supercluster position in the ECAL. The requirement on the photon p T is applied after the vertex assignment, which is described in further detail in Section 5.3. The preselection thresholds are shown in Table 1. Additionally, photons are required to satisfy at least one of R 9 > 0.8, I ch /p γ T < 0.3, and I ch < 20 GeV. The preselection efficiency is measured with the tag-and-probe technique using Z → ee events in data, while the efficiency of the electron veto is measured in Z → µµγ events in data.

Photon identification
Photons in events passing the preselection criteria are further required to satisfy a photon identification criterion based on a BDT trained to separate genuine ("prompt") photons from jets mimicking a photon signature. This ID BDT is trained on a simulated sample of γ +jet events, where prompt photons are used as the signal, while jets are used as the background. Input variables to the ID BDT include shower shape variables, isolation variables, the photon energy and η, and global event variables sensitive to pileup, such as the median energy density per unit area ρ [12].
Simulated inputs for the photon ID BDT, both shower shape and isolation variables, are corrected to agree with data using a chained quantile regression (CQR) method [59]. This method was developed to improve the agreement in the photon ID BDT output between data and simulation, thus reducing the size of the associated systematic uncertainty relative to previous analyses. Corrections are derived using an unbiased set of electrons from Z → ee events selected with a tag-and-probe method. The CQR comprises a set of BDTs that predict the cumulative distribution function (CDF) of a given input variable. Its prediction is conditional upon three electron kinematic variables (p T , |η|, φ) and ρ. The CDFs extracted in this way from data and simulated events are then used to derive a correction factor to be applied to any given simulated electron. These correction factors morph the CDF of the simulated shower shape onto the one observed in data.
The CQR method accounts for correlations among the shower shape variables and adjusts the correlation in the simulation to match the one observed in data. To achieve this, an ordered chain of the shower shape variables is constructed. The CDF of the first shower shape variable is predicted solely from the electron kinematic variables and event ρ values, while the corrected values of the previously processed shower shape variables are also added as inputs for subsequent predictions. The order of the different shower shape variables in the chain is optimised to minimise the final discrepancy of the ID BDT score between data and simulation.
The isolation variables are not included in the chain since their correlation with the shower shape variables is negligible. Furthermore, there is a p T threshold on the particle candidates included in the computation of the isolation variables. This causes these variables to follow a disjoint distribution, with a gap present between the peak at zero and a tail at positive values. The CDF of the isolation variables are therefore constant over the range of values between zero and the start of the tail, which prevents the use of the same technique used for the shower shape variables. The CQR method is thus extended with additional BDTs that are used to match, again based on the electron kinematic variables and the event ρ value, the relative population of the peak and tail between data and simulation. The tails of the isolation variable distributions themselves are then morphed using the same technique for the shower shape variables.
A systematic uncertainty associated with the corrections is also included in the analysis. This is estimated by rederiving the corrections with equally sized subsets of the Z → ee events used for training. Its magnitude corresponds to the standard deviation of the event-by-event differences in the corrected ID BDT output score obtained with the two training subsets. This uncertainty reflects the limited capacity of the network arising from the finite size of the training set. The size of the resulting experimental uncertainty is smaller than that required to cover discrepancies between data and simulation in previous versions of this analysis.
The distribution of the photon ID BDT for the lowest scoring photon for signal events and the different background components is shown in Fig. 3, together with a comparison of data and simulation using Z → ee events where the electrons are reconstructed as photons. These Z → ee events are chosen because of the similarity in the detector signature and reconstruction procedures for electrons and photons. Here, the electrons being reconstructed as photons means that the track information is not used, and the energy is determined using the algorithm and corrections corresponding to photons rather than electrons. The photon ID BDT distribution is also checked with photons using Z → µµγ events, where data and simulation are found to agree within uncertainties.
As an additional preselection criterion, photons are required to have a photon identification BDT score of at least −0.9. Both photons pass this additional requirement in more than 99% of simulated signal events. The efficiency of the requirement in simulation is corrected to match that in data using Z → ee events, and a corresponding systematic uncertainty is introduced.

Diphoton vertex identification
The determination of the primary vertex from which the two photons originate has a direct impact on the m γγ resolution. If the position along the beam axis (z) of the interaction producing the diphoton is known to better than around 1 cm, the m γγ resolution is dominated by the photon energy resolution.
The RMS of the distribution in z of the reconstructed vertices in data in 2016-2018 varies in the range 3.4-3.6 cm. The corresponding distribution in each year's simulation is reweighted to match that in data.
The diphoton vertex assignment is performed using a BDT (the vertex identification BDT) whose inputs are observables related to tracks recoiling against the diphoton system [12]. It is trained on simulated ggH events and identifies a single vertex in each event.
The performance of the vertex identification BDT is validated using Z → µ + µ − events. The vertices are refitted with the muon tracks omitted from the fit, to mimic a diphoton system. Figure 4 (left plot) shows the efficiency of correctly assigning the vertex, as a function of the dimuon p T . The data and simulation agree to within approximately 2% across the entire p T range. Nonetheless, the simulation is subsequently corrected to match the efficiencies measured in data, whilst preserving the total number of events. A systematic uncertainty is introduced with a magnitude equal to the size of this correction.  Figure 3: The left plot shows the distribution of the photon identification BDT score of the lowest scoring photon in diphoton pairs with 100 < m γγ < 180 GeV, for data events passing the preselection (black points), and for simulated background events (red band). Histograms are also shown for different components of the simulated background. The blue histogram corresponds to simulated Higgs boson signal events. The right plot shows the same distribution for Z → ee events in data and simulation, where the electrons are reconstructed as photons. The statistical and systematic uncertainty in simulation is also shown (pink band). Photons with an identification BDT score in the grey shaded region (below −0.9) are not considered in the analysis. The full data set collected in 2016-2018 and the corresponding simulation are shown.
The efficiency of assigning the diphoton vertex to be within 1 cm of the true vertex in simulated H → γγ events is approximately 79%. The events with an incorrectly-assigned vertex are primarily ggH events with zero additional jets, and the associated systematic uncertainty affects ggH events only.
A second vertex-related multivariate discriminant, the vertex probability BDT, estimates the probability that the vertex, chosen by the vertex identification BDT, is within 1 cm of the vertex from which the diphoton originated. The vertex probability BDT is trained on simulated H → γγ events using input variables relating to the vertices in the event, their vertex identification BDT scores, the number of photons with associated conversion tracks, and the p T of the diphoton system. Agreement is observed between the average vertex probability and the vertex efficiency in simulation, as shown in Fig. 4 (right plot).

Additional objects
Objects in the event other than the two photons are reconstructed as described in Section 2. Charged hadrons originating from interaction vertices other than the one chosen by the vertex identification BDT are removed from the analysis. In addition, all jets are required to have p T > 25 GeV, be within |η| < 4.7, and be separated from both photons by ∆R(jet, γ) > 0.4. Depending on the analysis category, more stringent constraints on the jet p T and |η| may be imposed; this is described in the text where relevant. In addition, some analysis categories require that jets also pass an identification criterion designed to reduce the number of selected jets originating from pileup collisions [60]. Jets from the hadronization of bottom quarks are  tagged using a DNN that takes secondary vertices and PF candidates as inputs [61].
Electrons and muons are used in the analysis categories targeting ttH and leptonic VH production. Electrons are required to have p T > 10 GeV and be within |η| < 2.4, excluding the barrel-endcap transition region. Muons must have p T > 5 GeV and fall within |η| < 2.4. In addition, isolation and identification requirements are imposed on both [62,63].

Event categorisation
The event selection in all analysis categories requires the two leading preselected photon candidates to have p γ1 T > m γγ /3 and p γ2 T > m γγ /4, respectively, with an invariant mass in the range 100 < m γγ < 180 GeV. The requirements on the scaled photon p T prevent distortions at the lower end of the m γγ spectrum. As described in Section 3, events are divided into analysis categories to provide sensitivity to different production mechanisms and STXS bins. Each analysis category is designed to select as many events as possible from a given STXS bin, or set of bins, referred to here as the target bin or bins. The requirements for each analysis category should also select as few events from other, non-targeted STXS bins as possible, to enable simultaneous measurements of different cross sections. Finally, the selection should also reject as many background events as possible, to maximise the measurements' eventual sensitivity. This section describes the several different categorisation schemes used for different event topologies, and the relevant STXS bins for each.
The STXS bins themselves are defined using particle-level quantities. In all targeted bins, |y H | is required to be less than 2.5. Jets are clustered using the anti-k T algorithm [19] with a distance parameter of 0.4. All stable particles, except for those arising from the decay of the Higgs boson or the leptonic decay of an associated vector boson, are included in the clustering. Jets are also required to have p T > 30 GeV. The definition of leptons includes electrons, muons, and tau leptons. Further details of the objects used to define the STXS bins can be found in Ref. [10].
In many of the categorisation schemes, ML algorithms are used to classify signal events or discriminate between signal and background processes. The output scores of the algorithms can then form part of the selection criteria used to define analysis categories. Where these ML techniques are used to classify events, two types of validation are performed. Firstly, in the typical case where simulated signal and background events are used to train the algorithm, a comparison of the simulated background to the corresponding data is performed. Good agreement between the two gives confidence that the background processes are accurately modelled and therefore that the ML algorithm performs well in its classification task. Since the background model used in the final maximum likelihood fit is derived directly from data, poor agreement in background-like regions cannot induce any biases, but only result in sub-optimal performance of the classifier. The second form of validation involves finding a signal-like region in which to compare the classifier output scores in simulation and data. Here the aim of the comparison is to instil confidence that simulated Higgs boson signal events, which do enter the final measurement, are sufficiently well-modelled. Therefore simulation and data should be expected to agree within statistical and systematic uncertainties in these cases. Furthermore, for all of the classifiers, the input variables are chosen such that m γγ cannot be inferred. This prevents distortion of the m γγ spectrum when applying selection thresholds on the output scores.
A summary of all the analysis categories, together with the STXS bin or bins each analysis category targets, is given in Section 6.6.

Event categories for ggH production
The definitions of the ggH STXS bins are given in Table 2, corresponding to the blue entries in Fig. 1. The bins are defined using p H T , the number of jets, and m jj . Those bins with p H T > 200 GeV are referred to as "BSM" bins because they have a cross section that is predicted to be low in the SM, but which could be enhanced by the presence of additional BSM particles. Events originating from ggZH production in which the Z boson decays hadronically are included in the definition of ggH. Analysis categories are defined to target each ggH STXS bin independently, except for those in the VBF-like phase space. Events from the VBF-like bins are categorised separately, as described in Section 6.2.
The ggH categorisation procedure can be summarised as follows. First, events are classified using the so-called ggH BDT. The ggH BDT predicts the probability that a diphoton event belongs to a given ggH STXS class. Each class corresponds either to an individual STXS bin or to a set of multiple STXS bins. The first eight classes considered by the ggH BDT are individual STXS bins. These comprise the zero, one, and two jet bins with p  Table 2: Definition of the ggH STXS bins. The product of the cross section and branching fraction (B), evaluated at √ s = 13 TeV and m H = 125 GeV, is given for each bin in the last column. The fraction of the total production mode cross section from each STXS bin is also shown. Events originating from ggZH production, in which the Z decays hadronically, are grouped together with the corresponding ggH STXS bin in the STXS measurements and are shown as a separate column in the table. The bbH production mode, whose σ SM B = This assignment is performed using the event's reconstructed p γγ T value. Finally, the analysis' sensitivity is maximised by further dividing the analysis categories using the diphoton BDT, which is trained to discriminate between signal and background processes and described in further detail below.
The ggH BDT is trained using simulated ggH events only. Input features to the ggH BDT are properties of the photons and quantities related to the kinematic properties of up to three p T > 20 GeV jets. The photon features used are the photon kinematic variables, ID BDT scores, m γγ resolution estimates, and the vertex probability estimate. The p γγ T value is also included as an input. As previously mentioned, the set of variables is chosen such that m γγ cannot be inferred from the inputs; for this reason, the p T /m γγ values of each photon, rather than p T , are used. The variables related to jets include the kinematic variables and pileup ID scores of the three leading jets in the event.
The resulting STXS bin assignment performs better than simply using the reconstructed p γγ T and number of jets. The fraction of selected ggH events in simulation that are assigned to the correct STXS bin increases from 77 to 82% when using the ggH BDT rather than the reconstructed p γγ T and number of jets. This improvement can be explained by the fact that the ggH BDT is able to exploit the correlations between the photon and jet kinematic properties. In this way, the well-measured photon quantities can be used to infer information about the less well-measured jets. As a result, the contamination of analysis categories due to migration across jet bins is reduced; the migration across p γγ T boundaries is much smaller and essentially unchanged by the ggH BDT. The ggH BDT therefore slightly improves the analysis sensitivity, most noticeably in the zero-and one-jet bins. Furthermore, the correlations between cross section parameters in the final fits are reduced.
To validate the modelling of the ggH BDT and its input variables, the agreement in the STXS class prediction between data and simulation in Z → ee events, with electrons reconstructed as photons, is checked. Figure 5 shows the number of events predicted to belong to each event class. The uncertainties in the photon ID BDT, the photon energy resolution, and the jet energy scale and resolution are included. There is good agreement between data and simulation in this signal-like region.
The diphoton BDT is used, after events are classified by the ggH BDT, to reduce the background from SM diphoton production, thereby maximising the analysis sensitivity. The diphoton BDT is trained with all Higgs boson signal events against SM diphoton production as background. A high score is assigned to events with photons showing signal-like kinematic properties, good m γγ resolution, and high photon identification BDT score. The input variables to the classifier are the photon kinematic variables, ID BDT scores, m γγ resolution estimates and the vertex probability estimate. Figure 6 shows the output score of the diphoton BDT for signal and background events, together with corresponding data from the m γγ sidebands, meaning 100 < m γγ < 120 GeV or 130 < m γγ < 180 GeV. A validation of the diphoton BDT obtained in Z → ee events, where the electrons are reconstructed as photons, is also shown in Fig. 6. Here the data and simulation agree within the statistical and systematic uncertainties.
After being classified by the ggH BDT, events are divided into analysis categories using the diphoton BDT, with the boundaries chosen to maximise the expected sensitivity. The resulting analysis categories are referred to as "tags". For ggH production, there is at least one tag targeting each individual STXS bin, except for the VBF-like bins. The tag names are given in decreasing order of the expected ratio of signal-to-background events (S/B). For example, the Events / class CMS 137 fb -1 (13 TeV) Predicted STXS class  Table 3. The yields shown in this and subsequent tables correspond to those in the final analysis categories, meaning that events selected by analysis categories with higher priority are not considered.

Event categories for VBF production
In the STXS framework, the qqH production mode includes both VBF events and VH events where the vector boson decays hadronically. Within qqH production, there are five STXS bins that correspond to typical VBF-like events, with a single bin for VH-like events. The precise definitions of the qqH STXS bins are given in Table 4. These correspond to the orange entries in Fig. 1.
Events with a dijet system characteristic of the VBF production mode have a dedicated categorisation scheme in this analysis, described in this section. Those events where the dijet is instead consistent with the decay of a vector boson are categorised separately, as described in Section 6.3. No analysis categories are constructed to target the zero or one jet qqH STXS bins, nor those with m jj < 60 GeV or 120 < m jj < 350 GeV.
Following the STXS binning scheme, the particle-level definition of the VBF-like dijet system Table 3: The expected number of signal events for m H = 125 GeV in analysis categories targeting ggH production, excluding those targeting the VBF-like phase space, shown for an integrated luminosity of 137 fb −1 . The fraction of the total number of events arising from each production mode in each analysis category is provided, as is the fraction of events originating from the targeted STXS bin or bins. Entries with values less than 0.05% are not shown. Here qqH includes contributions from both VBF and hadronic VH production, whilst "Top" includes ttH and tH together. The σ eff , defined as the smallest interval containing 68.3% of the m γγ distribution, is listed for each analysis category. The final column shows the expected ratio of signal to signal-plus-background, S/(S+B), where S and B are the numbers of expected signal and background events in a ±1σ eff window centred on m H .  requires two jets with p T > 30 GeV, and whose m jj > 350 GeV. These bins are defined analogously for EW qqH production as well as from ggH production. When constructing the corresponding analysis categories at reconstruction level, a dijet preselection is applied that requires two jets within |η| < 4.7, with p T > 40(30) GeV for the leading (subleading) jet, in addition to m jj > 350 GeV. Jets are also required to pass a threshold on a pileup identification score.
The so-called dijet BDT is trained to estimate the probability that an event passing the VBF preselection originated from VBF, ggH, or non-Higgs boson SM diphoton production. The inputs to the dijet BDT include various jet kinematic and angular variables, as well as the p T /m γγ of each photon and angular variables involving both jets and photons. These inputs for VBF, ggH, and non-Higgs boson SM production of two prompt photons are taken from simulation. However, the modelling of backgrounds, where at least one of the two photons is a misreconstructed jet, is poor, predominantly due to the fact that very few simulated events pass the selection criteria. In this analysis, an approach is adopted whereby the simulated background events with nonprompt photons are replaced with data from a dedicated control sample.
The control sample is defined using the sideband of the photon ID BDT distribution, by requiring at least one photon ID BDT score to be below −0.5. The events in this control sample can potentially have both a different normalisation and different kinematic properties from those in the signal region. To correct this, the events are reweighted in bins of the p T and |η| of the photon passing the ID BDT requirement. The required weights are derived from simulation, by estimating both the fraction of background events that contain nonprompt photons and the ratio of the expected number of events in the signal region to the control sample. The product of these two factors is applied as a weight to each data event in the control sample, and these reweighted events are subsequently used to train the dijet BDT.
The resulting distributions of the dijet BDT input variables are compared to the m γγ sideband data and are found to be in reasonable agreement. Furthermore, the increase in the number of events available for the training of the dijet BDT leads to an improvement in its discrimination power.
The two independent output probabilities of the dijet BDT, taken to be the VBF probability and the ggH probability, are validated in Z → ee + jets events with the electrons reconstructed as photons. The dijet preselection criteria required to enter the VBF-like analysis categories are also applied to the Z → ee + jets events. The VBF probability distribution in simulation and data sidebands is shown in the left plot of Fig. 7, while the right plot demonstrates good agreement between data and simulation in Z → ee + jets events. A similar level of agreement is observed in the ggH probability distribution.
Due to the use of the data control sample with photon ID BDT score below −0.5 in the dijet BDT training, an additional requirement that the two photons have a photon ID BDT score of larger than −0.2 is placed on events entering the VBF-like analysis categories. The final analysis categories are constructed following the structure of the STXS binning scheme. Events can be assigned to analysis categories targeting one of five VBF-like STXS bins, as shown in Events are further divided into analysis categories using both the dijet BDT output probabilities and the diphoton BDT score. For each of the five STXS bins, a set of analysis categories is constructed with events originating from VBF production considered as signal. An optimisation is performed defining lower bounds on the dijet VBF probability and diphoton BDT score, with an upper bound on the dijet ggH probability. Two analysis categories are constructed to target each STXS bin, the expected composition of which is given in Table 5.
An additional set of analysis categories is defined covering the four STXS bins with p H T < 200 GeV, but considering ggH events as signal instead of VBF. Two analysis categories targeting the set of four STXS bins together are constructed. Here lower bounds are set on the dijet ggH probability and diphoton BDT score, with an upper bound placed on the dijet VBF probability. The expected composition of these is also given in Table 5.

Event categories for hadronic VH production
In the EW qqH STXS binning scheme, there is a bin representing hadronic VH production, defined at the particle level by 60 < m jj < 120 GeV. Analysis categories targeting this bin are constructed in a similar way to those targeting VBF-like dijet events. The principal difference is in the selection of the two jets. The hadronic VH preselection requires two jets within |η| < 2.4 and with p T > 30 GeV, and satisfying a pileup jet identification criterion. In addition, the reconstructed m jj is required to be consistent with a decay of a vector boson, 60 < m jj < 120 GeV.
A BDT referred to as the VH hadronic BDT is trained with VH hadronic events as signal, against ggH and non-Higgs boson SM diphoton production together as background. The training events of VH, ggH, and SM production of two prompt photons are taken from simulation. The remaining background containing nonprompt photons is derived from a control sample in the same way as that employed for the dijet BDT training. The control sample is defined by requiring that at least one photon has a photon ID BDT score of less than −0.5, but otherwise passes the VH hadronic preselection. The resulting events are weighted to reproduce the expected number of background events and used in the BDT training of the VH hadronic BDT.
The input variables for the VH hadronic BDT are similar to those for the dijet BDT. Variables that aid in identifying events consistent with the vector boson decay are added, including the cosine of the difference of two angles: that of the diphoton system in the diphoton-dijet centreof-mass frame, and that of the diphoton-dijet system in the lab frame.
The final two analysis categories use the output scores of both the VH hadronic BDT and the diphoton BDT to increase sensitivity, in addition to requiring a photon ID BDT score of greater than −0.2.
The output score of the VH hadronic BDT in simulation and data sidebands is shown in the left plot of Fig. 8. The VH hadronic BDT is also validated in Z → ee + jets events with electrons reconstructed as photons, after the VH hadronic preselection is applied. The two distributions in simulation and data are shown in the right plot Fig. 8 and exhibit good agreement.
The expected signal and background yields in each VBF and hadronic VH analysis category are shown in Table 5.

Event categories for leptonic VH production
The analysis categories described here target events in which the Higgs boson is produced in association with a W or Z vector boson that subsequently decays leptonically. Depending on the particular leptonic decay mode of the vector boson, the possible final states include zero, one, or two charged leptons. The full definitions of each VH leptonic STXS bin are given in Table 6, corresponding to the green entries in Fig. 1. The bins are defined using p V T and the number of jets in the event.
For each of the three channels, a dedicated BDT classifier is used to discriminate between the VH signal and background events. Each of these three BDTs are trained on simulated signal and background events. The exception is the zero-lepton final state, for which some simulated backgrounds are replaced by events derived from data, as described below. The simulated SM background processes include photons plus jets, Drell-Yan, diboson production, and top quark pair production. The production modes of the Higgs boson other than VH are also treated as backgrounds. Where there are a sufficient number of expected signal events, the categorisation regions are further split into analysis categories sensitive to merged groups of STXS bins.
The categorisation region with two same-flavour reconstructed leptons in the final state focuses on the Z( )H production mode. Additional selection criteria are imposed to select two leptons consistent with the decay of a Z boson, including a requirement that the dilepton mass (m ) is between 60 and 120 GeV.
The so-called ZH leptonic BDT is used to discriminate the Z( )H signal events from backgrounds including both other Higgs boson production modes and non-Higgs-boson SM pro- cesses. Its input variables are kinematic properties of the photons, leptons, and jets present in the event, including angular variables describing the separation between the photons and leptons. In addition, jet identification variables such as the b tag score are used as inputs, which helps to discriminate against backgrounds containing top quarks.
The distributions of the ZH leptonic BDT score for simulated signal and background events, along with the same for the data sidebands, are shown in Fig. 9. With the available data set, this categorisation region is not sensitive to the corresponding individual STXS bins. For this reason, further splitting of the analysis categories is not performed. The sensitivity to inclusive leptonic ZH production is maximised by defining two analysis categories using the BDT score.
To gain sensitivity to the W( ν)H production mode, events with one reconstructed lepton are selected. Additional selection criteria are applied on the photon ID BDT to further reject background events containing nonprompt photons, and on the invariant mass of the reconstructed lepton with each photon to reduce the contamination of Drell-Yan events with an electron misidentified as a photon.
With this selection, the WH leptonic BDT is trained with simulated W( ν)H signal events against other Higgs boson modes and non-Higgs-boson SM backgrounds. The input features of the WH leptonic BDT are similar to those used in the ZH leptonic BDT, including photon, lepton, and jet kinematic variables. In addition, the transverse mass of the leading lepton and p miss T are used. The distributions of the WH leptonic BDT score for the signal and background  In each case, the signal and background simulation are shown as histograms with the data as black points. Events are taken from the m γγ sidebands, satisfying either 100 < m γγ < 120 GeV or 130 < m γγ < 180 GeV. The statistical uncertainty in the data points is denoted as vertical bars and that on the background simulation by the pink band. The simulated signal and background distributions are normalised to the luminosity of the data. To increase its visibility, the signal is scaled by a factor of 500 for the VH MET BDT, with a factor of 50 applied for both ZH leptonic and WH leptonic BDTs. The regions shaded grey are not considered in the analysis. The full data set collected in 2016-2018 and the corresponding simulation are shown. Table 5: The expected number of signal events for m H = 125 GeV in analysis categories targeting VBF-like phase space and VH production in which the vector boson decays hadronically, shown for an integrated luminosity of 137 fb −1 . The fraction of the total number of events arising from each production mode in each analysis category is provided, as is the fraction of events originating from the targeted STXS bin or bins. Entries with values less than 0.05% are not shown. Here ggH includes contributions from the ggZ(qq )H and bbH production modes, whilst "Top" represents both ttH and tH production together. The σ eff , defined as the smallest interval containing 68.3% of the m γγ distribution, is listed for each analysis category. The final column shows the expected ratio of signal to signal-plus-background, S/(S+B), where S and B are the numbers of expected signal and background events in a ±1σ eff window centred on m H . simulation samples and data sidebands is shown Fig. 9.
This single-lepton final state is sensitive to a reduced set of STXS bins. Three sets of analysis categories are defined, with p γγ T thresholds at 75 and 150 GeV. The p γγ T variable is used because it provides the most accurate estimate of the particle level p V T used to define the STXS bins; the presence of a neutrino in the final state means that the vector boson itself cannot be fully reconstructed. The sensitivity to each set of STXS bins is optimised by deriving analysis categories based on the WH leptonic BDT score. Two analysis categories are constructed with p γγ T < 75 GeV and 75 < p γγ T < 150 GeV, whilst one analysis category is defined with p γγ T > 150 GeV. The analysis categories targeting VH production where there are no reconstructed leptons in the event are referred to as the VH MET tags. These analysis categories receive contributions from both the Z(νν)H and W( ν)H production modes. In addition to vetoing events with leptons, p miss T > 50 GeV is required and the azimuthal angle between the diphoton system and p miss T must be greater than two radians.
With this selection the VH MET BDT is trained to discriminate between signal and background processes. The input features of the VH MET BDT rely on the same diphoton variables as in the ZH and WH leptonic BDTs, together with p miss T and jet variables. One of the dominant backgrounds in this final state consists of γ +jets events where one of the jets is misidentified as a photon. The simulation does not model this process well and the number of such events Table 6: Definition of the VH leptonic STXS bins. The product of the cross section and branching fraction (B), evaluated at √ s = 13 TeV and m H = 125 GeV, is given for each bin in the last column. The fraction of the total production mode cross section from each STXS bin is also shown. Unless stated otherwise, the STXS bins are defined for |y H | < 2.5. Events with |y H | > 2.5 are mostly outside of the experimental acceptance and therefore have a negligible contribution to all analysis categories. Only leptonic decays of the W and Z bosons are included in these definitions.

STXS bin
Definition units of p available is limited. Hence the γ +jets background component is modelled from a sample of data events where one of the photon candidates fails to satisfy the photon ID BDT requirement of −0.9. To enable this, a control sample is constructed by inverting the requirement on the photon ID BDT score. These events otherwise fulfil the full set of selection requirements for the VH MET BDT channel. A new value of the photon ID BDT score is generated for each event. This is achieved by assigning a random value drawn from the photon ID BDT distribution of simulated γ +jets events which pass the full set of selection criteria. The events are then appropriately weighted and used in the VH MET BDT training instead of the corresponding simulated samples. This is the same method first developed for the analysis described in Ref. [24], but differs to the method used in the training of the VBF and VH hadronic BDTs. The resulting increased number of events on which to train, as well as the improved modelling of the input variable distributions, improves the performance of the VH MET BDT.
The distributions of the VH MET BDT output score for the signal and background simulation samples together with the same for the data sidebands are shown Fig. 9.
The final expected signal and background yields for each ZH leptonic, WH leptonic, and VH MET analysis category are shown in Table 7.

Event categories for top quark associated production
The coupling between the Higgs boson and the top quark affects H → γγ cross sections both via ggH production, entering in the gluon loop, and via decay in the diphoton decay loop. In addition, the coupling can be accessed directly by measuring the rate of H → γγ events when the Higgs boson is produced in association with one or more top quarks. The obser- Table 7: The expected number of signal events for m H = 125 GeV in analysis categories targeting Higgs boson production in association with a leptonically decaying W or Z boson, shown for an integrated luminosity of 137 fb −1 . The fraction of the total number of events arising from each production mode in each analysis category is provided, as is the fraction of events originating from the targeted STXS bin or bins. Entries with values less than 0.05% are not shown.
Here ggH includes contributions from the ggZ(qq )H and bbH production modes, qqH incorporates both VBF and VH production with hadronic vector boson decays, and "Top" represents both ttH and tH production together. The σ eff , defined as the smallest interval containing 68.3% of the m γγ distribution, is listed for each analysis category. The final column shows the expected ratio of signal to signal-plus-background, S/(S+B), where S and B are the numbers of expected signal and background events in a ±1σ eff window centred on m H . vation of ttH production in the diphoton decay channel was recently reported by CMS and ATLAS [24,64]. There, multivariate discriminants are trained separately for hadronic and leptonic decays of the top quarks to construct analysis categories enriched in ttH events. In this analysis, the same techniques for the event categorisation are used. Additional analysis categories are constructed to provide sensitivity to individual STXS bins, the definitions of which are given in Table 8. These correspond to the purple entries in Fig. 1 for ttH, and the single yellow entry for tH.
Production of the Higgs boson in association with a single top quark is also measured in this analysis. A dedicated analysis category enriched in tHq events where the top decays leptonically is constructed. The tHq leptonic and ttH leptonic final states are very similar; an effort is therefore made to distinguish between the two.
A DNN referred to as the top DNN is trained with ttH as signal and tHq as background. It is used both by the tHq leptonic tag to reduce ttH contamination, and by the ttH leptonic analysis categories to reduce the contamination from tHq. The tHq leptonic tag is considered first in the tag priority sequence because of its lower expected signal yield. Each of the three categorisation regions (tHq leptonic, ttH leptonic, and ttH hadronic) then uses a dedicated discriminant referred to as BDT-bkg. The purpose of the BDT-bkg is to reduce backgrounds from non-Higgs-boson SM diphoton production and split events further by expected S/B into the final analysis categories.
For an event to be considered for the tHq leptonic analysis category, it must have at least one lepton, at least one b-tagged jet, and at least one additional jet. The top DNN and the tHq leptonic BDT-bkg are trained with these selection criteria applied. The top DNN takes both kinematic information from individual objects characteristic of top decays and global event information as inputs. The objects considered are the six leading jets and two leading leptons in p T . The four-momenta, along with the b tagging score and lepton identification scores, are included for each object. The global event features include the p miss T , number of jets, and photon kinematic variables and identification scores.
The tHq leptonic BDT-bkg uses similar input variables to distinguish tHq events from non-Higgs boson SM backgrounds, both of which are taken from simulation to perform the training. Kinematic variables and b tag scores for the three leading jets and b-tagged jets in p T are considered, as well as photon kinematic variables, and angular variables relating the jet and photon directions.
The distributions of the output scores for both the top DNN and the tHq BDT-bkg are shown in Fig. 10. In both cases, the agreement between data and simulation in this background-like region is imperfect. However, this does not affect the results of this analysis because the final background model is derived directly from data.
The final analysis category is defined by placing a requirement on both the score of the top DNN and the tHq leptonic BDT-bkg. Due to the low expected tHq signal yield, only one analysis category is constructed.
The analysis categories targeting ttH production are divided into two channels, representing either fully hadronic or leptonic decays of the tt system. The hadronic channel is defined by having zero isolated leptons, whilst the leptonic channel requires one or more isolated leptons, meaning it includes events where one or both top quarks decay leptonically. In the hadronic channel, three or more jets must be present, of which at least one is tagged as originating from a bottom quark. The leptonic channel requires the presence of one or more jets, and also includes a loose requirement on the top DNN to reject tHq events. This loose preselection for both channels maximises the available number of events for the training of the BDT-bkg in each channel and the top DNN in the leptonic channel.
For each channel, the BDT-bkg is trained on simulated signal and background events. The exception is that in the hadronic channel, γ+jets events are modelled from data. This provides both an improved description of the input features and a greater number of events on which to train the BDT-bkg. The procedure used to derive these events is identical to that described in Section 6.4. The inputs to the ttH BDT-bkg discriminants in each channel are kinematic properties of the jets, leptons, photons, and diphoton pair. It is not possible to infer the diphoton mass from the inputs. In addition to these features, the outputs of dedicated DNNs designed to reject specific backgrounds and the output of a dedicated "top quark tagger BDT" are used [65].
The additional DNNs are trained with ttH signal events against one source of background only. There are three such DNNs in total: one for each of the γγ+ jets and tt +γγ backgrounds in the hadronic channel, and one for the tt +γγ background in the leptonic channel. These backgrounds are chosen because both they are well-modelled in simulation and because it is possible to generate a high number of simulated events on which to train. Furthermore, tt +γγ events in particular are the principal background in the analysis categories most sensitive to ttH production. With these sufficiently large training samples, the background-specific DNNs are able to exploit features such as the full four-momentum vectors of physics objects. Adding these features directly to the inputs of the BDT-bkg do not improve its performance; the DNNs are required as an intermediate step to utilise this information effectively.
The top quark tagger BDT is designed to distinguish events with top quarks decaying into three jets from events that do not contain top quarks. It is trained on jet triplets from simulation of tt events, with inputs related to the kinematics, b tag scores, and jet shape information. The signal is jet triplets matched at generator-level to a top quark, and background is taken as random jet triplets [65].
The output distributions of the BDT-bkg for both the hadronic and leptonic channels are shown in Fig. 11. To validate the modelling of the BDT-bkg in each channel, a ttZ, Z → ee control region is used. The ttZ events have similar kinematical properties to ttH events, and are there- Table 9: The expected number of signal events for m H = 125 GeV in analysis categories targeting Higgs boson production in association with top quark, shown for an integrated luminosity of 137 fb −1 . The fraction of the total number of events arising from each production mode in each analysis category is provided, as is the fraction of events originating from the targeted STXS bin or bins. Entries with values less than 0.05% are not shown. Here ggH includes contributions from the ggZ(qq )H and bbH production modes, whilst qqH incorporates both VBF and hadronic VH production. The σ eff , defined as the smallest interval containing 68.3% of the m γγ distribution, is listed for each analysis category. The final column shows the expected ratio of signal to signal-plus-background, S/(S+B), where S and B are the numbers of expected signal and background events in a ±1σ eff window centred on m H . fore suitable for testing the agreement between data and simulation in the BDT-bkg score distributions. Additional requirements on the dielectron kinematics, number of jets, and number of b-tagged jets are imposed to increase the ttZ purity. The resulting comparisons between data and simulation are shown in Fig. 11 for the hadronic and leptonic channels.
Finally, events are split using the reconstructed p γγ T value to targeting specific STXS bins. Analysis categories are then defined through requirements placed on the BDT-bkg output, with the boundaries chosen to maximise the expected sensitivity to each bin.
The expected signal and background yields in each analysis category targeting top quark associated Higgs boson production are shown in Table 9.
BDT-bkg   syst. unc. ⊕ Stat. Figure 11: Distributions of BDT-bkg output used in the analysis categories targeting ttH production, for the leptonic (left) and the hadronic (right) channels. The upper two plots show events taken from the m γγ sidebands, satisfying either 100 < m γγ < 120 GeV or 130 < m γγ < 180 GeV. The lower two contain events from the ttZ control regions, described in the text. The grey region contains BDT-bkg scores below the lowest threshold for the ttH analysis categories. Total background uncertainties (statistical ⊕ systematic) are represented by the black (pink) shaded bands.

Summary of the event categorisation
The full set of analysis categories targeting the ggH, VBF, hadronic and leptonic VH, ttH, and tHq production mechanisms are summarised in Table 10. The different categorisation regions are shown in descending order of tag priority, starting with the tHq leptonic tag. If an event passes the selection criteria for more than one analysis category, it is assigned to the tag with the highest priority. Each STXS bin, or merged group of bins, and the number of analysis categories targeting it are shown.

Statistical procedure
The statistical procedure used in this analysis is identical to that described in Ref. [66], as developed by the ATLAS and CMS Collaborations. Simultaneous binned maximum likelihood fits are performed to the m γγ distributions of all analysis categories, in the range 100 < m γγ < 180 GeV. A likelihood function is defined for each analysis category using analytic models to describe the m γγ distributions of signal and background events, with nuisance parameters to account for the experimental and theoretical systematic uncertainties. The analytic signal model is derived from simulation, with a model constructed for each particle level STXS bin in each reconstructed analysis category. Both the shape and normalisation of the model are parametrised as functions of m H .
The background model is determined directly from the observed m γγ distribution in data. The analytic model for each analysis category can take one of a range of different functional forms, all of which represent a smoothly falling spectrum.
The best fit values and confidence intervals for the parameters of interest are estimated using a profile likelihood test statistic . (1) The likelihood functions in the numerator and denominator of Eq. (1) are constructed using the product over the likelihood functions defined for each analysis category. The quantitiesˆ α andˆ θ describe the unconditional maximum likelihood estimates for the parameters of interest and the nuisance parameters, respectively, whereasˆ θ α corresponds to the conditional maximum likelihood estimate for fixed values of the parameters of interest, α. In this analysis, the parameters of interest can be signal strengths, cross sections or coupling modifiers, depending on the fit being performed. In all fits, m H is fixed to its most precisely measured value of 125.38 GeV [58]. This choice is made to ensure that all measurements are reported with respect to the theoretical predictions consistent with the best available knowledge of m H . Further discussion of the implications of this choice and the difference with respect to profiling m H is given in Section 9.1.
The best fit parameter values,ˆ α, are identified as those that maximise the likelihood. For onedimensional measurements, such as the signal strength and STXS fits, the 68 and 95% confidence intervals are defined by the union of intervals for which q( α) < 0.99 and < 3.84, respectively. In the case where there are multiple parameters of interest in the fit, the intervals are determined treating the other parameters as nuisance parameters. For two-dimensional measurements, such as those performed to coupling modifiers in the κ-framework, the 68 and 95% confidence regions are defined by the set of parameter values that satisfy q( α) < 2.30 and < 5.99, respectively. To compute the SM expected results, the observed data is replaced by an Asimov data set generated with all parameter values set to the SM expectation [67].
The methods used to construct the signal and background models are described in detail in the remainder of this section.

Signal model
The signal shape for the m γγ distribution in each analysis category and for a nominal m H is constructed from the simulation of each production process.
Since the distribution of m γγ depends on whether the vertex associated with the candidate diphoton was correctly identified within 1 cm, the correct vertex and wrong vertex scenarios are considered separately when constructing the signal model. In a given analysis category, a separate function is constructed for events originating from each STXS bin in each vertex scenario, by fitting the m γγ distribution using a sum of at most five Gaussian functions. This choice provides sufficient flexibility in the fit whilst maintaining computational efficiency. The number of Gaussian functions is determined using an F -test [68], avoiding overfitting statistical fluctuations due to the limited size of the simulated samples.
The final fit function for each analysis category is obtained by summing the individual functions for all STXS bins in both vertex scenarios. Figure 12 shows signal models for each year individually, and for the sum of the three years together. The σ eff is defined as half of the smallest interval containing 68.3% of the m γγ distribution.

Background model
The model used to describe the background is extracted from data using the discrete profiling method [11,69]. This technique estimates the systematic uncertainty associated with choosing a particular analytic function to fit the background m γγ distribution. The choice of the background function is treated as a discrete nuisance parameter in the likelihood fit to the data.
A large set of candidate function families is considered, including exponential functions, Bernstein polynomials, Laurent series, and power law functions [69]. For each family of functions, an F -test [68] is performed to determine the maximum order of parameters to be used, while the minimum order is determined by placing a requirement on the goodness-of-fit to the data.
When fitting these functions to the m γγ distribution, the value of twice the negative logarithm of the likelihood (−2∆ ln L) is minimised. A penalty is added to −2∆ ln L to take into account the number of floating parameters in each candidate function. When making a measurement of a given parameter of interest, the discrete profiling method minimises the overall −2∆ ln L considering all allowed functions for each analysis category. Checks are performed to ensure that describing the background m γγ distribution in this way introduces negligible bias to the final results.

Systematic uncertainties
In this analysis, the systematic uncertainty associated with the background estimation from data is handled using the discrete profiling method, as described above. There are many systematic uncertainties that affect the signal model; these are handled in one of two ways. Uncertainties that modify the shape of the m γγ distribution are incorporated into the signal model as nuisance parameters, where the mean and width of each Gaussian function can be affected. These uncertainties are typically experimental uncertainties relating to the energy of the individual photons. Conversely if the shape of the m γγ distribution is unaffected, the uncertainty is treated as a log-normal variation in the event yield. These uncertainties include theoretical sources and experimental uncertainties such as those affecting the BDTs used to categorise events. The magnitude of each uncertainty's impact is determined individually for each STXS bin in each analysis category.

Theoretical uncertainties
Theoretical uncertainties affect both the overall cross section prediction for a given STXS bin, and the distributions of kinematic variables used in the event selection and categorisation. When measurements of cross sections are performed, the uncertainties in the overall cross sections are omitted, and instead are considered as uncertainties in the SM predictions. The uncertainties related to the event kinematic properties, which affect the efficiency and acceptance of the analysis, are still taken into account. Uncertainties affecting the overall cross section normalisations and those affecting the event kinematic properties are included when measuring signal strengths and coupling modifiers. When deriving the effect on the kinematic distributions, the impact on the STXS bin cross section normalisation is factored out to avoid double counting.
The sources of theoretical uncertainty considered in this analysis are as follows: • Renormalisation and factorisation scale uncertainties: the uncertainty arising from variations of the renormalisation and factorisation scales used when computing the expected SM cross section and event kinematic properties. These account for the missing higher-order terms in perturbative calculations. The recommendations provided in Ref. [10] are followed. The uncertainty in the overall normalisation is estimated by: varying the renormalisation by a factor of two; varying the factorisation scale by a factor of two; varying both in the same direction simultaneously. Depending on the production process, the size of the uncertainty in the overall normalisation varies from around 0.5% for VBF production to 15% for tHq production.
To estimate the uncertainty in the event kinematic properties, the distribution of events falling into each analysis category is recalculated when varying both the renormalisation and factorisation scales by a factor of two in the same direction simultaneously. The overall cross section for a given STXS bin is kept constant. These uncertainties, representing migrations between analysis categories, are decorrelated for different production modes and different regions of the Higgs boson phase space, resulting in 22 independent nuisance parameters. The migration uncertainties are in general around 1%. • Uncertainties in the ggH STXS fractions: for ggH production, additional sources are included that account for the uncertainty in the modelling of the p H T distributions, the number of jets in the event, and the ggH contamination in the VBF categories. A number of sources are introduced to reflect the migration of events around the p H T bin boundaries, at 10 GeV for zero-jet events and 60 and 120 GeV for events with at-least one jet, such that their magnitude depends on the number of jets and the p H T . An additional source covers the uncertainty in p H T in the Lorentz-boosted region arising from the assumption of infinite top quark mass in the ggH loop. This is determined by comparing the p H T distribution to the prediction from finite-mass calculations. Two further sources account for the migration between the zero, one, and two or more jet bins. The uncertainty in the ggH production of events with a VBF-like dijet system is covered by two sources corresponding to the prediction in the two-and three-jet-like bins. In addition, two nuisance parameters are introduced to account for migrations across the m jj bin boundaries, at 350 and 700 GeV. The total magnitudes of these uncertainties vary from around 5 to 30%, with events that have one or more jets and high values of p H T typically having the greatest associated uncertainty. These uncertainties affect the overall cross section normalisations, and so are attributed to the SM prediction in the cross section measurements. • Uncertainties in the qqH STXS fractions: similarly for qqH production, additional sources are introduced to account for the uncertainty in the modelling of the p H T , m jj , and p Hjj T distributions, and the number of jets in the event. A total of six sources are defined to reflect migrations of events across m jj boundaries at 60, 120, 350, 700, 1000, and 1500 GeV. Two additional nuisance parameters account for migrations across the p H T = 200 GeV and p Hjj T = 25 GeV bin boundaries. Finally, a single source is defined to account for a migration between the zero and one, and the two or more jet bins. In each case, the uncertainty is computed by varying the renormalisation and factorisation scales and recalculating the fractional breakdown of qqH STXS stage-1.2 cross sections. The total magnitude varies between bins but is at most 8%. Again, these are considered as uncertainties in the SM predictions when performing cross section measurements. • Uncertainties in the tt H STXS fractions: for ttH production, four nuisance parameters are used to account for the uncertainty in the p H T distributions. Each nuisance parameter represents migration across one of the boundaries at the p H T values of 60, 120, 200, and 300 GeV that define the ttH STXS bins. The magnitudes of these uncertainties are derived by varying the renormalisation and factorisation scales, and have values of up to 9%. When performing cross section measurements, the sources are treated as uncertainties on the SM prediction.
• Uncertainties in the VH leptonic STXS fractions: for VH leptonic production, additional sources are introduced to account for the uncertainty in the modelling of the p V T distributions, and the number of jets in the event. Four independent sources are defined to reflect the migrations of events across the p V T boundaries at 75, 150, and 250 GeV, in addition to the migration between the zero and greater than one-jet bins for events with p V T of 150-250 GeV. These sources are defined separately for the WH leptonic, ZH leptonic, and ggZH leptonic production modes, leading to 12 independent nuisance parameters. In each case, the uncertainty is computed by varying the renormalisation and factorisation scales and recalculating the fractional breakdown of VH leptonic STXS stage-1.2 cross sections. The total magnitude varies between bins but is at most 5% for the dominant WH and ZH leptonic production modes. Again, these are considered as uncertainties in the SM predictions when performing cross section measurements. • Uncertainty in the ggH contamination of the top quark associated categories: the theoretical predictions for ggH are less reliable in a regime where the Higgs boson is produced in association with a large number of jets. Three different contributions are considered: the uncertainty from the parton shower modelling, estimated by taking the observed difference in the jet multiplicity between MADGRAPH5 aMC@NLO predictions and data in tt + jets events [70], the uncertainty in the gluon splitting modelling, estimated by scaling the fraction of events from ggH with real b quark jets in simulation by the measured difference between data and simulation of σ(ttbb )/σ(ttjj) [71] and the uncertainty due to the limited size of the simulated samples.
The combined impact of these uncertainties in the top quark associated signal strength is about 2%. • Parton distribution function uncertainties: these account for the uncertainty due to imperfect knowledge of the composition of the proton, which affects the partons that are most likely to initiate high energy events. The overall normalisation uncertainties for each Higgs boson production process also include the uncertainty in the value of the strong coupling constant α S , and are taken from Ref. [10]. Uncertainties in the event kinematic properties are calculated following the PDF4LHC 100 prescription [47,[72][73][74] using the MC2HESSIAN procedure [75,76]. As with the renormalisation and factorisation scale uncertainties, the normalisation for a given STXS bin is kept constant when calculating the migrations between analysis categories. The overall normalisation uncertainties are 1-5%, with the migrations significantly smaller, usually less than 1%. • Uncertainty in the strong coupling constant: the uncertainty in the value of α S is included in the treatment of the PDF uncertainties, following the PDF4LHC prescription. The impact on the overall normalisation is largest for ggH production, with a value of 2.6%. An additional source is included to account for changes in the event kinematic properties due to the uncertainty in α S . This is calculated using a similar procedure to the renormalisation and factorisation scale migration uncertainties, but instead varying the value of α S , and corresponds to uncertainties that are in general less than 1%. • Uncertainty in the H → γγ branching fraction: the probability of the Higgs boson decaying to two photons is required to calculate the SM expected cross section, but this branching fraction is not known exactly. The uncertainty is currently estimated to be around 3% [10]. This uncertainty is included in the signal strength and coupling modifier measurements, and is considered an uncertainty in the SM predictions for cross section measurements. • Underlying-event and parton shower uncertainties: these uncertainties are obtained using dedicated simulated samples. The parton shower uncertainties originating from the modelling of the hadronization are evaluated by varying the renormalisation scale for QCD emissions in initial-state and final-state radiation by a factor of 2 and 0.5. The uncertainties in the modelling of the underlying-event are evaluated by varying the PYTHIA8 tune from that used in the nominal simulation samples, introduced in Section 4. Both these uncertainties are treated as variations in the relative contributions from each STXS bin for a given production mode, and therefore affect the STXS bin cross section normalisation. The impact is in general around 5%, but can be as large as 30% for bins corresponding to high p H T and high jet multiplicity. As described in Section 9.2, it is necessary to merge certain STXS bins when measuring cross sections to avoid large uncertainties or very high correlations between parameters. If two bins are measured individually, the theoretical uncertainty representing event migrations between the two bins are not included since both cross sections are being fitted. The act of merging bins across a boundary means the measurement is sensitive to the relative fraction of the two bins, and an uncertainty must be included to model this. As a result, the uncertainty sources accounting for migrations across the merged boundaries are included in the relevant cross section measurements.

Experimental uncertainties
The uncertainties that affect the shape of the signal m γγ distribution are listed below. These include uncertainties that account for the difference between photon showers and the electron showers used to derive the energy scale corrections. The combined effect of all signal model shape uncertainties in the measurement of the inclusive Higgs boson signal strength modifier is found to be about 2%.
• Photon energy scale and resolution: the uncertainties associated with the corrections applied to the photon energy scale in data and the resolution in simulation are evaluated using Z → ee events. The estimate is computed by varying the regression training scheme, the distribution of R 9 , and the electron selection criteria. For the majority of photons the resulting uncertainty in the energy scale is 0.05-0.15%, although for those with very high p T the effect can be 0.5-3.0%. • Nonlinearity of the photon energy scale: a further source of uncertainty covers possible remaining differences in the linearity of the photon energy scale between data and simulation. The uncertainty is estimated using Lorentz-boosted Z → ee events. In this analysis, an uncertainty of 0.2% on the photon energy scale is assigned, which accounts for the nonlinearity across the full range of photon p T values. • Shower shape corrections: an uncertainty in the shower shape corrections accounts for the imperfect modelling of shower shapes in simulation. The impact is estimated by comparing the energy scale before and after the corrections to shower shape variables, as described in Section 5, are applied. The magnitude of the uncertainty in the energy scale ranges from 0.01-0.15%, depending on the photon |η| and R 9 values. • Longitudinal nonuniformity of light collection: an uncertainty is associated with the modelling of the light collection as a function of emission depth within a given ECAL crystal. The calculation of this uncertainty is described in detail in Ref. [58]. The uncertainty is 0.16-0.25% for photons with R 9 > 0.96, whilst the magnitude for low R 9 photons is below 0.07%. • Modelling of material in front of the ECAL: the amount of material through which objects pass before reaching the ECAL affects the behaviour of the electromagnetic showers, and may not be perfectly modelled in simulation. Dedicated samples with variations in the amount of upstream material are used to estimate the impact on the photon energy scale. The magnitude of the resulting uncertainty ranges from 0.02-0.05% for the most central photons, increasing to as much as 0.24% for those in the endcap. • Vertex assignment: the largest contribution to the uncertainty in the fraction of events where the chosen vertex is smaller than 1 cm from the true vertex comes from the modelling of the underlying-event. In addition, the uncertainty in the ratio of data and simulation obtained using Z → µ + µ − events is incorporated. A nuisance parameter is included in the signal model that allows the fraction of events in each vertex scenario to vary by ±2%.
The uncertainties that only modify the event yield have an effect of around 4% on the inclusive Higgs boson signal strength modifier measurement. They include the set of sources described below.
• Integrated luminosity: uncertainties of 2.5, 2.3, and 2.5% are determined by the CMS luminosity monitoring for the 2016, 2017, and 2018 data sets [38][39][40], respectively, whilst the uncertainty on the total integrated luminosity of the three years together is 1.8%. The uncertainties for each data set are partially correlated to account for common sources in the luminosity measurement schemes. • Photon identification BDT score: the uncertainty arising from the photon identification BDT score is estimated by varying the set of events used to train the quantile regression corrections. It is seen to cover the residual discrepancies between data and simulation. The uncertainty in the signal yields is estimated by propagating this uncertainty through the full category selection procedure. The impact in the most sensitive analysis categories is around 3%. • Jet energy scale and smearing corrections: The energy scale of jets is measured using the p T balance of jets with Z bosons and photons in Z → ee, Z → µ + µ − , and γ +jets events, as well as the p T balance between jets in dijet and multijet events [77]. The uncertainty in the jet energy scale is a few percent and depends on p T and η. The impact of jet energy scale uncertainties in event yields is evaluated by varying the jet energy corrections within their uncertainties and propagating the effect to the final result. Correlations between years are introduced for the different jet energy scale uncertainty sources, ranging between 0 and 100%. The impact on the category yields is largest for those targeting VBF, hadronic VH and top quark associated production and can be as high as 22% for the scale uncertainties, but is less than around 8% for the resolution. • Per-photon energy resolution estimate: the uncertainty in the per-photon resolution is parametrised as a rescaling of the resolution by ±5% about its nominal value. This is designed to cover all differences between data and simulation in the distribution, which is an output of the energy regression. The maximum yield variation in an analysis category is around 5%, however for most categories the impact is below the percent level. • Trigger efficiency: the efficiency of the trigger selection is measured with Z → ee events using the tag-and-probe technique. The size of its uncertainty is less than 1%.
An additional uncertainty is introduced to account for a gradual shift in the timing of the inputs of the ECAL first level trigger in the region at |η| > 2.0, which caused a specific trigger inefficiency during 2016 and 2017 data taking [78]. Both photons and to a greater extent jets can be affected by this inefficiency. The resulting uncertainty is largest for the categories targeting VBF production, with a maximum impact on the yield of 1.4%. • Photon preselection: the uncertainty in the preselection efficiency is computed as the ratio between the efficiency measured in data and in simulation. Its magnitude is less than 1%. • Missing transverse momentum: this uncertainty is computed by shifting the reconstructed p T of the particle candidates entering the p miss T computation, within the momentum scale and resolution uncertainties appropriate to each type of reconstructed object, as described in Ref. [77]. In this analysis, the impact on the category yields is never larger than 5%, even for analysis categories that explicitly use the p miss T in their definition. • Pileup jet identification: the uncertainty in the pileup jet classification output score is estimated by comparing the score of jets in events with a Z boson and one balanced jet in data and simulation. The magnitude is of the order 1%. • Lepton isolation and identification: this uncertainty affecting electrons and muons is computed by varying the ratio of the efficiency in simulation to the efficiency in data and using the tag-and-probe technique in Z → ee events. The resulting impact on the categories selecting leptons is up to around 1%. • b jet tagging: uncertainties in the b tagging efficiency are evaluated by comparing data and simulated distributions for the b tagging discriminator. The uncertainties include the statistical component in the estimate of the fraction of heavy-and light-flavour jets in data and simulation. Its magnitude is around 3% for the analysis categories targeting top quark associated production, which make use of the b tagging discriminant.
Most of the experimental uncertainties are left uncorrelated among the different years. The exceptions are the partial correlations introduced for the integrated luminosity and jet energy correction uncertainties.

Results
The expected signal composition of the analysis categories in terms of a set of merged STXS bins is shown in Fig. 13. In the plot, the analysis categories targeting a common STXS region are summed, such that the signal compositions of the individual analysis categories are weighted according to the ratio of the numbers of signal to signal-plus-background events (S/S+B). The fractional contribution of the total signal yield in a given analysis category group arising from each process is shown.
The best fit signal-plus-background model is shown with data for the sum of all analysis categories in Fig. 14. Again each analysis category is weighted by (S/S+B), such that the absolute signal yield is kept constant.    Figure 13: The composition of the analysis categories in terms of a merged set of STXS bins is shown. The granularity of the STXS bin merging corresponds to the finest granularity used for the cross section measurements in this analysis. Analysis categories targeting a common STXS region are summed, where the signal compositions of the individual categories are weighted in the sum by the expected ratio of signal to signal-plus-background events. The colour scale corresponds to the fractional yield in each analysis category group (rows) accounted for by each STXS process (columns). Each row therefore sums to 100%. Entries with values less than 0.5% are not shown. Simulated events for each year in the period 2016-2018 are combined with appropriate weights corresponding to their relative integrated luminosity in data. The column labelled as "qqH rest" includes contributions from the qqH 0J, qqH 1J, qqH m jj < 60 GeV and qqH 120 < m jj < 350 GeV STXS bins.

Signal strength modifiers
A common signal strength modifier, µ, is defined as the ratio of the observed product of the Higgs boson cross section and diphoton branching fraction to the SM expectation. It is measured to be µ = 1.12 +0.09 −0.09 = 1.12 +0.06 −0.06 (theo) +0.03 −0.03 (syst) +0.07 −0.06 (stat). The uncertainty is decomposed into theoretical systematic, experimental systematic, and statistical components. The statistical component includes the uncertainty in the background modelling. The compatibility of this fit with respect to the SM prediction, expressed as a p-value, is approximately 17%.
In this fit, and in all subsequent fits, m H is fixed to its most precisely measured value of 125.38 GeV [58]. The precise determination of m H and the systematic uncertainties that enter its measurement are beyond the scope of this analysis. Nonetheless, the dependence of the measured signal strengths on m H is checked. Profiling m H without constraint, rather than fixing it to 125.38 GeV, has a small impact on the measured results; the best fit signal strength values change by 0.7-1.8%. In each case, the change is less than 10% of the measured uncertainty.
Signal strength modifiers for each Higgs boson production mode are also measured. Unlike the subsequent STXS fits described in Section 9.2, the VH hadronic and VH leptonic processes are grouped to scale according to µ VH , whereas the VBF production mode scales with µ VBF . The parameter µ top scales the ttH, tHq and tHW production modes equally and µ ggH scales both ggH and bbH production.
The resulting signal and background m γγ distributions after the fit using this parameter scheme are shown in Fig. 15. Analysis categories are divided into four groups, corresponding to those targeting the ggH, VBF, VH, and top quark production modes. In each group, the individual analysis categories are summed after weighting by S/(S+B).
The values of the production mode signal strength modifiers and their uncertainties are displayed in Fig. 16. The precision of these measurements is significantly improved from previous analyses performed by the CMS Collaboration in the H → γγ decay channel. In particular, the measurement of the µ VH signal strength modifier has improved substantially from that shown in Ref. [12], beyond what would be expected from the increase in the size of the data set alone. The p-value of the production mode signal strength modifier fit with respect to the SM prediction is approximately 50%.
The main sources of systematic uncertainty affecting the signal strength modifier in each production mode are summarised in Fig. 17. The dominant contributions to the measurement uncertainty in the µ ggH , µ VH and µ top signal strength modifiers originate from the corresponding renormalisation and factorisation scale uncertainties, whereas the underlying event and parton shower uncertainties are the dominant sources of uncertainty in the µ VBF measurement. The largest experimental uncertainties originate from the integrated luminosity, the photon identification, and the photon energy measurement for the µ ggH and µ VH signal strength modifiers. The uncertainties in the jet energy scale and resolution have a larger impact on µ VBF and µ top , where µ top has an additional large contribution from the uncertainty in the b tagging.

Simplified template cross sections
This section details the fits performed to extract cross sections within the STXS framework and their respective 68% confidence level (CL) intervals. The theoretical uncertainties in the normalisation of the signal parameters are not included in the cross section measurements. In each fit, ggZH events in which the Z boson decays hadronically are grouped with the corresponding ggH STXS bin. All bbH events are treated as part of the ggH 0J high p H T bin. The hadronic VH processes are grouped with VBF production to form the qqH parameters. Parameters which are not measured are constrained to their SM prediction, within theoretical uncertainties. These are the zero jet, one jet, m jj < 60 GeV, and 120 < m jj < 350 GeV bins in the qqH binning scheme.
Two different parameterisations are considered, with varying levels of granularity defined by the merging of certain STXS bins. It is necessary to merge bins to avoid either very large uncertainties in some parameters or very high correlations between parameters. Merging fewer bins keeps the model-dependence of the results as low as possible, as no additional assumptions are made about the relative contributions of different STXS bins. The results with reduced modeldependence however have larger uncertainties in the measured cross section parameters.
In this paper, the results of two different fits to the cross sections of partially merged STXS bins are reported. The first is referred to as the "maximal" merging scenario, where in general STXS bins are merged until their expected uncertainty is less than 150% of the SM prediction. The second "minimal" merging fit instead merges as few bins as possible whilst ensuring that parameters do not become too anti-correlated, meaning values of less than around 90%.
The maximal merging scheme defines 17 parameters of interest. The VBF-like regions (≥2-jets,    Figure 16: Observed results of the fit to signal strength modifiers of the four principal production modes. The contributions to the total uncertainty in each parameter from the theoretical systematic, experimental systematic, and statistical components are shown. The colour scheme is chosen to match the diagram presented in Fig. 1. The compatibility of this fit with respect to the SM prediction, expressed as a p-value, is approximately 50%. Also shown in black is the result of the fit to the inclusive signal strength modifier, which has a p-value of 17%. m jj > 350 GeV) in the ggH and qqH schemes are merged to define the ggH VBF-like and qqH VBF-like parameters, respectively. The four bins with p H T > 200 GeV in the ggH scheme are merged into a single bin, labelled as ggH BSM. Additionally, the WH leptonic, ZH leptonic and ttH bins are all fully merged into single parameters. The ZH leptonic parameter groups both ZH and ggZH production.
The minimal merging scheme defines a more granular fit with 27 parameters of interest. The qqH VBF-like region is fully split into the four STXS bins defined by the boundaries at m jj = 700 GeV and p Hjj T = 25 GeV. To avoid large correlations between the fitted parameters, the four ggH VBF-like bins are merged with the corresponding bins in the qqH scheme. Additional splittings are included in the ggH scheme at p H T = 300 and 450 GeV, and the WH leptonic scheme at p V T = 75 and 150 GeV. Furthermore, the ttH region is split into five parameters according to the boundaries at p H T = 60, 120, 200, and 300 GeV. Table 11 summarises the maximal and minimal merging schemes by listing the STXS bins that contribute to each parameter of interest. The STXS bins that are constrained to their respective SM predictions in both fits are also listed.
The best fit cross sections and 68% CL intervals are shown for the two merging schemes in Figs. 18 and 20. The corresponding numerical values are given in Tables 12 and 13. For both the maximal and minimal fits, the statistical component of the uncertainty dominates for all measured cross sections. Overall, the results from both merging scenario fits are in agreement with SM predictions; the p-values with respect to the SM predictions are 31 and 70% for the maximal and minimal merging scenarios, respectively.
In the maximal merging scenario, ggH production with p H T > 200 GeV, which is particularly sensitive to BSM physics entering the ggH loop, is measured to a precision of less than 50%, relative to the SM prediction. The cross section is found to be consistent with the SM, with a measured value of 0.9 +0.4 −0.3 relative to the SM prediction. In addition, the product of the tH production cross section times H → γγ branching fraction is measured to be 1.3 +0. 8 −0.7 fb, corresponding to an excess of 6.3 +3.4 −3.7 times the SM expectation. Using the CL s procedure [79], a rate of tH production of 14 (8) times the SM expectation is observed (expected) to be excluded at the 95% CL.
The minimal merging scenario fit represents the current most granular cross section measurement performed in a single Higgs boson decay channel, showing reasonable sensitivity to many different regions of Higgs boson production phase space. In particular, the results contain the first measurements of ttH production in bins of p H T . The size of the uncertainty in each of the four ttH bins with p H T < 300 GeV is less than 100% of the SM prediction. Additionally, ggH production with p H T > 200 GeV is measured in three separate regions. The three corresponding cross sections are all measured to be within one standard deviation of the respective SM expectations.
Correlations between the fitted parameters are presented in Figs. 19 and 21. The correlations for the ggH parameters are observed to be small between adjacent p H T bins and larger between adjacent number of jet bins. This results from the fact that p γγ T is a well-measured quantity, whereas reconstructing the number of jets in an event is more difficult. Nevertheless, the application of the ggH BDT in the event categorisation helps to minimise these correlations. The largest correlations in the maximal merging scheme exist between the ggH VBF-like and qqH VBF-like cross sections and the ttH and tH cross sections, with values of −0.76 and −0.59, respectively. These result from the sizeable contamination of ggH VBF-like events in the qqH analysis categories, and the contamination of ttH events in the tHq leptonic category. The act of splitting ttH production into five separate parameters in the minimal merging scenario introduces larger correlations into the measurement.  Here the tH cross section ratio has a different scale, due to its high best fit value and uncertainty. The cross sections are constrained to be non-negative, as indicated by the hashed pattern below zero. The parameters whose best fit values are at zero are known to have 68% CL intervals which slightly under-cover; this is checked to be a small effect using pseudo-experiments. The colour scheme is chosen to match the diagram presented in Fig. 1. The compatibility of this fit with respect to the SM prediction, expressed as a p-value, is approximately 31%.

Coupling modifiers
The κ-framework defines coupling modifiers to directly parametrise deviations from the SM expectation in the couplings of the Higgs boson to other particles [14]. Two different likelihood scans, each with two dimensions, are performed. Full details of each parameterisation are given in Ref. [8].
In the first fit, the resolved κ model is used. Here the scaling factors of loops present in Higgs boson production and decay are resolved into their SM components, in terms of the other κ parameters. The most important of these are in ggH production and H → γγ decay, but others, such as the loop in ggZH production, are also resolved. The results of a two-dimensional scan in κ V and κ F , scaling the Higgs boson coupling to vector bosons and to fermions, respectively, are shown in the upper plot of Fig. 22. The H → γγ decay rate contains an interference term proportional to κ V κ F . This means that the rate of ggH and ttH production (∝ κ 2 F ), relative     Here the tH cross section ratio has a different scale, due to its high best fit value and uncertainty. The cross sections are constrained to be non-negative, as indicated by the hashed pattern below zero. The parameters whose best fit values are at zero are known to have 68% CL intervals which slightly under-cover; this is checked to be a small effect using pseudo-experiments. The colour scheme is chosen to match the diagram presented in Fig. 1. The orange lines dashed with blue for the VBF-like parameters represent contributions from both the ggH and the qqH STXS bins. The compatibility of this fit with respect to the SM prediction, expressed as a p-value, is approximately 70%.   to the rate of VBF and VH production (∝ κ 2 V ), can be used to gain sensitivity to the relative sign of the tt-H and VV-H couplings. In addition, the tHq and tHW production modes also include a term proportional to κ V κ F . This analysis explicitly targets tHq production via the tHq leptonic analysis category, the inclusion of which helps to further reduce the degeneracy between positive and negative κ F values. The region with negative values of κ F is observed (expected) to be excluded with a significance of 0.5 (2.4) standard deviations. The reduction in the observed significance with respect to the expected is due to the observed excess in the tH production mode cross section.
A second fit is performed using the unresolved κ model, where the ggH and H → γγ loops are given their own effective scaling factors denoted κ g and κ γ , respectively. The κ g and κ γ parameters are particularly sensitive to additional BSM states, that contribute towards the rate of Higgs boson production and decay via loop processes. The observed result of a two-dimensional scan in these two parameters is shown in the lower plot of Fig. 22. In the scan, the other κ parameters in the unresolved model are fixed to unity. The best fit point is consistent with the SM expectation at approximately the 68% CL.

Summary
Measurements of Higgs boson properties with the Higgs boson decaying into a pair of photons are reported. Events with two photons are selected from a sample of proton-proton collisions at a centre-of-mass energy √ s = 13 TeV collected with the CMS detector at the LHC from 2016 to 2018, corresponding to an integrated luminosity of 137 fb −1 . Analysis categories enriched in events produced via gluon fusion, vector boson fusion, vector boson associated production, production associated with two top quarks, and production associated with one top quark are constructed.
A range of production and coupling properties of the Higgs boson are measured. The total Higgs boson signal strength, relative to the standard model (SM) prediction, is measured to be 1.12 ± 0.09. A simultaneous measurement of the signal strengths of the four principal Higgs boson production mechanisms is performed and found to be compatible with the SM prediction with a p-value of 50%. Two different measurements are performed within the simplified template cross section framework, in which 17 and 27 independent kinematic regions are measured simultaneously, with corresponding p-values with respect to the SM of 31 and 70%, respectively. Many of these kinematic regions are measured for the first time, including a simultaneous measurement of Higgs boson production in association with two top quarks in five different regions of the Higgs boson transverse momentum p H T . Furthermore, several additional measurements are the most precise made in a single channel to date. These include cross sections of vector boson fusion in different kinematic regions, gluon fusion in association with jets, and the region of gluon fusion production with p H T > 200 GeV, which is particularly sensitive to physics beyond the SM. The gluon fusion cross section with p H T > 200 GeV is found to be consistent with the SM, with a measured value of 0.9 +0.4 −0.3 relative to the SM prediction. An upper limit on the rate of Higgs boson production in association with a single top quark is also presented. The observed (expected) limit at 95% confidence level is found to be 14 (8) times the SM prediction. All other results, such as measurements of the Higgs boson's couplings to vector bosons and to fermions, are also in agreement with the SM expectations.  Figure 22: Observed two-dimensional likelihood scans performed in the κ-framework: κ V -vsκ F in the resolved κ model (upper) and κ γ -vs-κ g in the unresolved κ model (lower) . The 68 and 95% CL regions are given by the solid and dashed contours, respectively. The best fit and SM points are shown by the black cross and red diamond, respectively. The colour scale indicates the value of the test statistic.