Measurement of the Higgs boson inclusive and differential fiducial production cross sections in the diphoton decay channel with pp collisions at $\sqrt{s}$ = 13 TeV

The measurements of the inclusive and differential fiducial cross sections of the Higgs boson decaying to a pair of photons are presented. The analysis is performed using proton-proton collisions data recorded with the CMS detector at the LHC at a centre-of-mass energy of 13 TeV and corresponding to an integrated luminosity of 137 fb$^{-1}$. The inclusive fiducial cross section is measured to be $\sigma_\mathrm{fid}$ = 73.4 $_{-5.3}^{+5.4}$ (stat) ${}_{-2.2}^{+2.4}$ (syst) fb, in agreement with the standard model expectation of 75.4 $\pm$ 4.1 fb. The measurements are also performed in fiducial regions targeting different production modes and as function of several observables describing the diphoton system, the number of additional jets present in the event, and other kinematic observables. Two double differential measurements are performed. No significant deviations from the standard model expectations are observed.

and √ s = 8 TeV [6] and it is designed to provide a measurement of the H → γγ cross section as independent as possible from the theoretical assumptions used for its extraction. To this end, the definition of the fiducial region mimics the detector acceptance and the events are categorized at the reconstruction level using a per-event mass resolution estimator. The large data sets collected so far allow for an approximate doubling of the number of bins for most observables compared with the previous publications and therefore the most precise measurement of the pp → H + X, H → γγ fiducial cross section to date. The measurements presented in this paper can be seen as a complementary and less model dependent addition to the STXS measurements described in Ref. [23].
The fiducial differential cross sections are extracted in observables of the diphoton system, the transverse momentum (p T ) of the jets produced in association with the diphoton system, the dijet system containing the two leading-p T jets, and with respect to several other event-level observables. For the first time, the differential cross section is also measured as a function of the angular observable φ * η , as defined in Eq.
(3), which complements the diphoton p T measurement at low p T values, and τ j C , defined in Eq. (4), which is the jet p T weighted by a function of its rapidity. Furthermore, two double-differential measurements are performed in bins of the diphoton p T and the number of jets, and in bins of the diphoton p T and τ j C . The paper is structured as follows. An overview of the CMS detector is given in Section 2. The data and simulation samples used to perform this analysis are described in Section 3. Section 4 summarizes the reconstruction of H → γγ candidate events and contains a detailed description of a newly developed simulation correction method. The selection criteria for events entering this analysis are shown in Section 5 and their categorization in Section 6. An overview of the observables as a function of which the H → γγ cross section is measured, and the definitions of the fiducial phase spaces are given in Section 7. In Section 8, the statistical procedure used to extract the results from data is presented, while the relevant systematic uncertainties are given in Section 9. The results are reported in Section 10 and the paper is summarized in Section 11. Tabulated results are provided in the HEPData record for this analysis [29].

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, recording the tracks of charged particles up to |η| < 2.5, 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 (EB) and 1.48 < |η| < 3.0 in the two endcap regions (EE). Preshower detectors consisting of two planes of silicon sensors interleaved with a total of three radiation lengths of lead are located in front of each EE detector. Forward calorimeters extend the pseudorapidity coverage up to |η| = 5. Muons are detected in gas-ionization chambers embedded in the steel flux-return yoke outside the solenoid with a coverage up to |η| < 2.4. The excellent resolution of the photon energy provided by the ECAL allows for a relatively large signal significance and therefore a precise Higgs boson cross section measurement.
Events of interest are selected using a two-tiered trigger system [30]. 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 latency of about 4 µs. The second level, known as the high-level trigger (HLT), consists of a farm of processors running a version of the full event reconstruction software optimized for fast processing, and reduces the event rate to around 1 kHz before data storage.
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. [31].

Data samples and simulated events
The analysis uses data corresponding to the integrated luminosities of 36.3 fb −1 , 41.5 fb −1 , and 59.4 fb −1 collected in 2016, 2017, and 2018, respectively. The integrated luminosities of these data-taking periods are individually known with uncertainties in the 1.2-2.5% range [32][33][34], while the total (2016-2018) integrated luminosity has an uncertainty of 1.6%, the improvement in precision reflecting the (uncorrelated) time evolution of some systematic effects. In this section, the data and simulated event samples for all three years are described, highlighting the differences between the years.
Events from pp collisions are selected using a diphoton trigger at the HLT with asymmetric p T thresholds on the leading and subleading photon of 30 (18) and 30 (22) GeV in 2016 (2017-2018) data, respectively. The invariant mass of the photon pair is required to be m γγ > 90 GeV, while each photon is also required to pass loose requirements on the calorimetric isolation and on the shape of the electromagnetic shower of the photon candidates. The trigger efficiency is measured on Z → ee events using the "tag-and-probe" technique [35]. These events are collected with a single-electron trigger. The efficiency measured in data, in bins of p T , the energy fraction in the centre of the electromagnetic cluster R 9 [36], and η, is used to weight the simulated events to replicate the trigger efficiency observed in data. The R 9 variable is defined as the fraction of energy deposited in three by three crystals matrix around the seed E 3×3 with respect the total uncorrected energy E SC,raw in the supercluster: The R 9 variable is useful to identify photons and electrons that started showering in the material upstream of the ECAL.
Simulated signal samples, corresponding to the four production modes with the largest cross section (ggF, VBF, VH, ttH) are generated using MADGRAPH5 aMC@NLO (version 2.6.5) at next-to-leading order (NLO) accuracy [37] in perturbative quantum chromodynamics (QCD). Events produced via the gluon fusion mechanism are weighted as a function of the Higgs boson p T and the number of jets in the event, to match the prediction from the next-to-NLO including parton showering event generator (NNLOPS) program [38][39][40]. The parton-level samples are interfaced with PYTHIA 8 (version 8.240) [41] for parton showering and hadronization, with the CUETP8M1 [42] (CP5 [43]) tune used for the simulation of 2016 (2017-2018) data. The signal cross sections and branching fractions recommended by the LHC Higgs Working Group [15] are used.
The dominant source of background events in this analysis is from QCD diphoton (γγ) production. The majority of the remaining background originates from γ+jet or dijet events, in which jets are misidentified as photons.
The diphoton background (γγ) is generated with the SHERPA (version 2.2.4) generator [44]. It includes the Born processes with up to three additional jets, as well as the box processes at leading order (LO) accuracy. The γ+jet and dijet background are simulated with PYTHIA 8. For the optimization of the event categorization, the γ+jet and γγ background samples are used. The training of the photon identification multivariate analysis (MVA) is performed with the γ+jet sample. As explained in Section 8, no simulated samples for the background are used in the measurement of the fiducial inclusive or differential cross sections.
Samples of Z → ee, Z → µµ, and Z → µµγ simulated events are generated with MAD-GRAPH5 aMC@NLO at NLO precision and used for the derivation of energy scale, resolution, and photon selection correction factors.

Event reconstruction
Physics objects are reconstructed using the particle-flow algorithm (PF) [45], which aims to identify each individual particle (PF candidate) in an event, with an optimized combination of information from the various elements of the CMS detector.

Vertex identification
Finding the correct primary vertex in the H → γγ final state when the Higgs boson was produced through gluon fusion faces the complication that photons (unless they convert into an electron-positron pair in the material upstream of the ECAL) do not leave any ionization signal in the tracker that can be used for finding their trajectory, and therefore the production vertex.
On the other hand, assigning the right vertex is crucial for measuring the kinematic variables of the diphoton system and keeping a good diphoton mass resolution. It is found that if the reconstructed vertex is within 1 cm in distance along the beam axis of the true vertex, the angular component in the mass resolution is subdominant compared to the ECAL energy resolution.
To obtain the vertex position of an H → γγ event, a dedicated vertex finding algorithm has been developed based on a boosted decision tree (BDT) [13]. The BDT is trained on simulated events produced via gluon fusion, using track information from the objects recoiling against the diphoton system and tracks from electron-positron pairs originating from a converted photon if available. Its performance is validated in data with Z → µµ events, in which the vertices are refitted omitting the muon tracks. The average vertex finding efficiency is 75-80% across the years of data-taking [13].

Photon reconstruction and identification
The reconstruction of a photon starts by clustering the energy deposited in the ECAL. The supercluster formation begins by finding the ECAL particle-flow clusters that locally record the largest energy above a given threshold. Then the clusters spatially compatible with originating from the same prompt photon are merged into a supercluster. Photons are selected as superclusters that are not linked to any charged-particle trajectory associated to a primary interaction vertex. A variable supercluster size in φ recovers the energy of photons initiating an electromagnetic shower in the material upstream of the ECAL. A more detailed description of the photon reconstruction algorithm can be found in Ref. [36].

Photon energy corrections
The response of each ECAL crystal is calibrated individually and adjusted for time-dependent effects. A semi-parametric regression is used to correct for the partial collection of the electromagnetic shower energy in the supercluster and effects from simultaneous soft pp collisions (pileup), and to predict the per-photon energy resolution. This regression is implemented with a BDT and trained on simulated photons using several shower shape variables as input [13].
After this simulation-based correction, a time-dependent correction of the energy response measured in data to the one predicted by simulation is derived from Z → ee events, where the photon reconstruction is applied to the electron candidates. The typical time granularity of this correction is of one LHC fill and is necessary to track the response evolution of the ECAL. In this way any residual energy drift in data is corrected as a function of time. The Z → ee events are used to adjust the energy scale in data and the resolution in simulation: the energy of the electrons in data is corrected to have the peak of their invariant mass distribution matching the one in the Drell-Yan simulation; the energy in simulated electron pairs is instead smeared such that the width of their invariant mass spectrum matches the one observed in data. Further details about this procedure can be found in Ref. [46].

Photon identification
To discriminate prompt photons against those originating from the decay of neutral particles in jets, e.g. from π 0 or η mesons, a classifier implemented with a BDT, called the photon identification MVA, is trained on a γ+jet simulation sample. The BDT input consists of several variables associated with the shape of the supercluster associated with the photon candidate, and isolation sums, built from neutral electromagnetic components and charged hadrons, calculated for different vertex choices [36]. For the calculation of the isolation from charged hadrons, it is necessary to consider different vertex choices, because for a given photon, the tracks (characterizing the charged hadrons) included in the isolation cone around it depend on the considered primary vertex. For photons detected in EE, the energy deposited into the preshower detector is considered by the photon identification MVA, in addition to the variables mentioned above.

Correcting the photon identification MVA imperfect modelling
The photon identification MVA is the primary means of reducing the γ+jet and dijet QCD backgrounds in the analysis. Its imperfect modelling in simulation results in one of the most important sources of systematic uncertainty for this analysis and was the dominant one for the previous CMS differential fiducial H → γγ cross section measurement [28]. The modelling inaccuracy comes from the imperfect modelling of the MVA input variables. Among the reasons for this are the ever changing ECAL conditions due to radiation damage, and the imperfect knowledge of the amount of material traversed by a photon before reaching the ECAL. Together with a campaign of studies to address the root causes of the mismatch [36], an effective correction applied to simulation has been developed for the analyses with diphoton final states [23, 47,48].
A dedicated method has been devised to compute a per-photon correction that morphs the input variable distributions in simulation, and their correlations, to match data. The method, called chained quantile regression (CQR), is used for the cluster shape and isolation variables used as input for the photon identification MVA. It is based on the multitarget regression algorithm described in Ref. [49] and implemented to correct these variables, conditional on the photon kinematic properties and the average spatial energy density of the event.
The method originates from the quantile morphing technique, in which any continuous distribution can be morphed into any other continuous distribution, using the following algorithm: 1. Find the value of the cumulative distribution function (CDF) in simulation F sim Y (y) for value y of the variable Y to be adjusted; 2. Get the point of the cumulative distribution function of the target with the same value F data      Figure 1: The chained approach for the set of input variables for the quantile BDTs. The input variables for the BDT for each variable being corrected are given. The variables to be corrected are used as the target, and the quantile is learned through the learning objective given in Eq. (1). Within one group of variables (y 1 , . . . , y n ), with nonnegligible correlations, an order is set. The quantile BDT for a given variable includes the prior set of variables, within this ordering, as additional inputs. For simulation (right), the additional input variables are corrected before using them as inputs for the quantile BDTs.

SIM
3. Assign the new value y C as the corrected value instead of y in the distribution to be corrected.
This can be written as: To be able to perform these steps, the CDFs (F data Y , F sim Y ) of the distribution to be corrected and the target have to be known.
For the purpose of correcting the set of input variables of the photon identification MVA, quantile morphing is extended to be conditional on the photon kinematic variables (p T , η, φ) and the average spatial event energy density ρ [13]. The conditional shape of the CDFs in simulation and data are approximated using 21 quantile regressors (implemented with BDTs, called quantile BDTs from here), each predicting the quantile for a different value of the CDF.
Apart from the kinematics of the photon and the event energy density, the photon identification MVA input variables can be sorted into three groups with nonnegligible correlations among the variables within them: the variables describing the shape of the electromagnetic shower, the isolation sums of charged candidates reconstructed by the PF algorithm, and the isolation sum of electromagnetic neutral objects. The quantile BDTs are trained using a set of input variables that depends on the variable to be corrected and the group it is assigned to. Within each group, a sequence of variables is defined and all variables appearing before the one to be corrected are added to the input variables used for the BDTs on top of the photon p T , η, and φ, and the event ρ. The set of input variables for the quantile BDTs for the variables within one of the groups introduced above for simulation and data is illustrated in Fig. 1. This chained approach allows the correlations within the groups of variables to be corrected at the same time as their conditional shape. The order in which the variables are inserted in the chain has been optimized on the correction performance. The training objective of the quantile BDTs is defined as [50]: where τ corresponds to the CDF value for which the quantile q Y is computed.
The training is performed independently in data and simulation using the SCIKIT-LEARN package [51] on electrons (reconstructed as photons) from simulated Z → ee events and data collected using a single-electron trigger. The tag-and-probe method is then applied to data, while for simulation the reconstructed electrons are matched to their true generated counterparts before the tag-and-probe method is applied. This gives the unbiased sample of probe electrons needed to derive the conditional CDF shapes.
The largest group of variables are the electromagnetic cluster shapes, for which the CQR approach introduced above is sufficient. The other groups of variables to be corrected include the photon and charged isolation sums which are discontinuous, preventing the immediate use of the CQR algorithm. The discontinuity is caused by detector thresholds applied to the energy of the constituents of the isolation sums, resulting in a fraction of events with exactly zero isolation energy. The isolation distributions present a peak at zero followed by a continuous tail. The number of events in the peak transforms into a plateau in the cumulative distribution function. Since the population of events with zero isolation sum is not necessarily the same in data and simulation, the codomains of the CDFs are not covering the same range of values, which is a necessary condition to apply the quantile morphing algorithm.
A dedicated procedure has been developed to overcome this complication and still apply the quantile morphing to the tail parts of the distributions. The mitigation strategy differs between the two types of isolation sum: the photon isolation is corrected independently from all other variables, while the two charged isolation sums, which are correlated, are corrected simultaneously.
For the photon isolation sum, the first step is equalizing the number of events in the first bin of the isolation distribution between data and simulation using a stochastic procedure to move events from the peak (zero isolation) to the tail of the distribution and vice-versa. The events to be moved are randomly picked with a probability calculated from the relative population between the peak and the tail. The probability to have an isolation sum of zero is estimated with a BDT trained with the binary cross-entropy loss function conditionally on the photon kinematic variables and ρ using the XGBOOST package [52]. For all isolation variables, a newly assigned value for an event moved to the tail part (isolation greater than zero) is obtained by sampling from the tail part of the distribution conditional on the photon kinematic variables and ρ. A set of example distributions for the photon isolation sum, one being uncorrected, one after equalizing the number of events with zero isolation, and the third one after applying the equalizing step and the CQR technique to the continuous tail, together with an illustration of the equalizing step are shown in Fig. 2.
The charged isolation computed with respect to the vertex giving it its largest value, called worst charged isolation, is by construction larger than or equal to the charged isolation calculated with respect to the chosen vertex (hereafter called charged isolation). This leads to three possible categories: both charged isolation sums being exactly zero, the charged isolation being zero and the worst charged isolation being larger than zero, and both of them being larger than zero. A classifier implemented with a BDT, using (p T , η, φ, ρ) as input, is trained to distinguish between them. The events to be moved between the three categories are randomly picked with a probability calculated from the output of this classifier such that their populations in  Figure 2: The upper pane shows the distribution of the photon isolation sum (I ph ) in 2018 data (black dots) and simulation (coloured histograms). The green histogram shows the uncorrected distribution, the orange one the distribution after equalizing the number of events with zero isolation in simulation with data, and the purple one the distribution after applying the equalizing step and the CQR technique to its tail part. The arrows show the two ways events can be shifted. From peak to tail (green) and from tail to peak (yellow) with their respective probabilities p(peak to tail) and p(tail to peak). The lower pane shows the ratio of the three simulation distributions to the one from data. The distributions shown in this figure are taken from a set of events distinct from the ones used for the derivation of the correction method.
simulation match the ones observed in data.
Finally, the continuous distributions of the two charged-isolation sums greater than zero are corrected using a CQR: the morphed charged-isolation is fed as an additional input to the quantile BDTs for the worst charged isolation, to correctly model their correlation.
After correcting the photon identification MVA input variables in simulation with the method just described, its output distribution matches the corresponding one in data within 1% for the ECAL barrel and within 5% for the ECAL endcap, as shown in Fig. 3. Consequently, the impact of the associated systematic uncertainty on the inclusive fiducial cross section measurement is significantly reduced from 5%, for the previous analysis [28], to 1.5%, when adding the per-year contributions to the uncertainty in quadrature.

Other objects
Electrons are identified from a charged-particle track and the ECAL energy clusters consistent with this track's extrapolation to the ECAL. Additionally, possible bremsstrahlung photons emitted along the way through the tracker material are taken into account. 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.  Figure 3: Distribution of the output of the photon identification MVA for the probe candidate in a Z → ee tag-and-probe sample for data and the MADGRAPH5 aMC@NLO simulation. The electrons have been reconstructed as photons and a selection to reduce the number of misidentified electrons in data is applied. The simulated events have been reweighted with respect to p T , η, φ, and ρ to match data in order to remove effects from mismodelled kinematic variables.
Electrons that are detected in the barrel (|η| < 1.4442) or endcap (|η| > 1.566) part of the ECAL and the corresponding distributions are shown on the left or right, respectively. The blue band shows the systematic uncertainty assigned to the data simulation mismatch of the output of the photon identification MVA. The orange histogram and points in the upper and lower plots, respectively, show the photon identification MVA distribution and its ratio to data evaluated using the uncorrected version of its input variables in simulation.
Muons are identified from tracks in the central tracker consistent with either a track or several hits in the muon system, and associated with calorimeter deposits compatible with the muon hypothesis. Their momentum 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 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.
Particle flow candidates are clustered to jets using the infrared and collinear-safe anti-k T clustering algorithm [53,54] with a distance parameter of 0.4. In this process, the contribution from each calorimeter tower is assigned a momentum, the absolute value and the direction of which are given by the energy measured in the tower, and the coordinates of the tower. The raw jet energy is obtained from the sum of the PF candidate energies, and the raw jet momentum by the vectorial sum of the candidate momenta, which results in a nonzero jet mass. Additional proton-proton interactions within the same or nearby bunch crossings can contribute additional tracks and calorimetric energy depositions, increasing the apparent jet momentum. To mitigate this effect, tracks identified to be originating from pileup vertices are discarded and an offset correction is applied to correct for remaining contributions. The raw jet energies are then corrected to establish a uniform relative response in η and a calibrated absolute response in p T . Jets originating from the hadronization of b quarks are identified using a tight selection on the output of the DEEPJET b tagger [55], corresponding to an efficiency of 56% for 2016 and 62% for 2017 and 2018.
The missing transverse momentum vector p miss T is computed as the negative vector sum of the transverse momenta of all the PF candidates in an event, and its magnitude is denoted as p miss T [56]. The p miss T is modified to account for corrections to the energy scale of the reconstructed jets in the event.

Event selection
The selection criteria for photons forming the diphoton candidates that enter the analysis are designed to be tighter than the trigger thresholds based on p T , isolation, and cluster shape variables [36]. The details of this photon preselection are given in Ref. [23].
Each candidate in the photon pair must have a supercluster with |η| < 2.5, where the region 1.4442 < |η| < 1.566 is excluded. The transverse momentum relative to the invariant mass of the diphoton pair (p T /m γγ ) of the leading-p T (subleading-p T ) photon has to be greater than 1/3 (1/4). Both photons must satisfy a minimal score of the photon identification MVA. If an event contains more than one photon pair, the one with the highest p T is selected.
Two selections of jets enter the analysis, one with |η| < 2.5, the other with |η| < 4.7. In both cases the jets are required to have p T > 30 GeV. The first selection is applied in differential measurements with respect to observables calculated when exactly one additional jet is present in the event, or with event-level observables, such as the number of jets and the number of b-tagged jets. The looser η criterion is used for observables taking into account variables involving two or more jets. Events with the tighter η requirement benefit from a better jet energy resolution because of tracker information being available for jets with |η| < 2.5. To ensure a mutually exclusive selection of photon and jet candidates, each jet is required to have a ∆R > 0.4, which is defined as the quadratic sum of the differences in azimuthal angle and pseudorapidity between a photon candidate and the jet: A jet reconstructed within |η| < 2.4, and satisfying the identification criteria of b jets described in Section 4.3, is considered a b jet at reconstruction level. At the particle level, a jet is considered as originating from a b quark if, among its constituent hadrons, the decay products of at least one b hadron are present.
Both electrons and muons are required to have p T > 20 GeV and |η| < 2.4. They have to be separated by ∆R > 0.35 from the photon candidates. As with the photons, electrons are not reconstructed in the ECAL region 1.4442 < |η| < 1.566. To reject Z+γ → e + e − γ events with one electron misidentified as a photon, the invariant mass of any electron-photon pair cannot be within 5 GeV of the Z boson mass. The loose identification criteria defined in Ref. [36] are applied to all electrons. Muons are selected when fulfilling a tight selection criterion [57] based on the quality of the track and the relative isolation, which is calculated as the sum of the transverse energy of charged and neutral hadrons, and photons, in a cone with radius parameter ∆R = 0.4 around the muon divided by the muon p T .

Event categorization
The photon pairs that enter this analysis are categorized according to the decorrelated mass resolution estimator, described below, and the photon identification MVA. The same boundary for the photon identification MVA selection is applied to the p T -leading and p T -subleading photons in events falling into any resolution category. These categorization criteria are designed to ensure minimal model dependence, while yielding the best possible analysis sensitivity.

Mass resolution estimator
The mass resolution estimator σ m /m for a photon pair is defined as the quadratic sum of the relative energy resolution of each photon as estimated by the photon energy regression BDT. The per-photon energy resolution in simulation is smeared to match the one observed in data using the Z → ee decay as a reference. The relative ECAL energy resolution (σ E /E) improves with larger photon energies, for which the impact of noise and stochastic effects on the resolution decreases. This effect leads to a correlation between the mass resolution estimator (σ m /m) and the invariant mass of the diphoton system. Such a correlation, combined with a selection on σ m /m, can distort the shape of the background by depleting the low-mass region for highresolution categories, creating a rising edge in the background shape that spoils the assumption of a monotonically falling background shape, which is necessary to model it using the discrete profiling method presented in Section 8.2. To avoid this effect, the mass resolution estimator is decorrelated with respect to the diphoton invariant mass. This procedure is based on applying the quantile morphing algorithm described in Section 4.2.3 to the mass dependent σ m /m spectrum and begins by dividing the σ m /m distribution into mass bins. The σ m /m distribution in each mass bin is then morphed to match the shape of the target, which is the σ m /m distribution at a reference mass value of m γγ = 125 GeV. In the following, the decorrelated relative mass resolution estimator is referred to as σ D m . More details on this procedure can be found in Ref. [6].

Decorrelated relative mass resolution categorization
The boundaries of the σ D m categories are optimized simultaneously with the photon identification MVA selection, using the expected signal significance as a figure of merit. The optimization is performed using simulated samples of signal events at m H = 125 GeV, and diphoton and γ+jet events to model the background, simulating different detector conditions for each year of data taking. For all years, it is found that the optimal number of σ D m categories is three. Photon pairs above the highest σ D m boundary are discarded. The σ D m boundaries roughly group the photon pairs by the ECAL regions where each of the two photons were detected. The overall efficiency of the photon identification criterion ranges from 75 to 80% across the three years of data taking. The effect of the event categorization is to increase the signal significance by between 13 and 18% depending on the particular reconstruction level bin. A detailed breakdown, together with the efficiencies of the σ D m categorization, is given in Table 1. Examples of the signal shape in different categories are shown in Fig. 4.

Observables and fiducial phase space definitions
The measurements of the inclusive and differential production cross sections are performed in fiducial phase spaces to reduce their model dependence. A set of kinematic and identification selections is applied at the level of stable particles obtained from the event generator and parton shower simulation (particle-level) to the single photons and to the diphoton system, as summarized in Table 2. The particle-level information does include QCD+EWK final state radiation effects but not the simulation of the interaction of the stable particles with the detector. The I γ gen is defined as the sum of the energy of all stable hadrons produced in a cone of radius ∆R = 0.3 around the photon. This requirement is introduced to mimic the photon identification MVA selection criterion at the reconstruction level. In addition to the baseline fiducial phase space, other fiducial phase spaces are defined specifically to target certain observables of interest. There are several observables that are calculated with variables of the leading-p T , or the leading-and subleading-p T additional jets. All additional jets considered in this analysis are required to satisfy p j T > 30 GeV. If the measured observable requires one additional jet to Table 1: Efficiencies of the photon identification MVA and σ D m categories for events taken from the signal sample for all three years of data taking. The second row shows the efficiency of the photon identification MVA selection in the three σ D m categories and for the full sample (Overall). The third row shows the efficiencies of the selections for the three σ D m categories without the photon identification MVA selection applied. The forth row reports the effective width (σ eff ) of the Higgs boson signal in each category. The four dominant Higgs boson production modes considered for this analysis are included in the sample and m H = 125 GeV is used. Only events satisfying the fiducial selection are included.   be evaluated, the 1-jet phase space is used. It requires at least one particle-level jet being in the event and the leading-p T jet to be within |η j | < 2.5. An observable that is calculated involving two additional jets is measured in the 2-jet phase space, for which the |η| selection is relaxed to |η j | < 4.7 for all jets. A region within the 2-jet phase space is designed to predominantly accept events originating from the vector boson fusion (VBF) Higgs boson production process and is called VBF-enriched. The additional requirements for events to enter this phase space region are: n jets ≥ 2, m jj > 200 GeV, and ∆η jj > 3.5, where m jj is the invariant mass of the dijet system containing the two leading-p T jets and ∆η jj is the difference in pseudorapidity between these two jets. The selection for the VBF-enriched region is designed with loose cut-based kinematic criteria to ensure the model dependence of the measurement is kept as small as possible.
All differential cross section observables, along with the phase space regions in which they are measured, are given in Table 3. The bin boundaries have been chosen such that the statistical uncertainty is similar for the cross section measured in each bin, with the goal of a per-bin statistical uncertainty of around 30%. Where applicable, the bin boundaries have been unified with the ones for the ATLAS H → γγ differential cross section measurement [25]. The same bin boundaries are used for the corresponding reconstruction-level quantities. The measurements are performed binned in the transverse momentum and absolute rapidity of the diphoton system, which are sensitive to the Higgs boson production mechanism, the modelling of the QCD radiation, and the parton distribution functions (PDFs) of the proton. Two additional observables are defined for the diphoton system. One is the cosine of the polar angle in the Collins-Soper reference frame of the diphoton system cos θ * [58], which is sensitive to the spin and CP properties of the diphoton resonance. The other is |φ * η |, designed to probe the low p T region while minimizing the impact of experimental uncertainties. It is defined as [59,60] where the acoplanarity angle is related to the angle between the photons in the transverse plane, ∆φ, as φ acop = π − ∆φ, and θ * η is the scattering angle of the two photons with respect to the proton beam in the reference frame in which the two particles are back to back in the (r, θ) plane. Observables involving the two p T -leading jets, such as their transverse momentum or absolute rapidity, are sensitive to the Higgs boson production in QCD. In the system of the two photons originating from the Higgs boson and the leading jet, the difference in azimuthal angle and rapidity between the diphoton system and the jet are taken as observables. These specifically probe the properties of the VBF production mechanism. The τ j C observable is also considered. It is defined as [61]: and calculated for the (up to) six largest p T jets in the event, taking the maximum among them. It represents the transverse momentum of the jet weighted by a function that depends on the rapidity of the jet and smoothly suppresses the contribution from forward jets. In the denominator of Eq. (4), the difference between the rapidity of the jet y j and the rapidity of the Higgs boson y H is taken. This means that τ j C is calculated in the frame where y H = 0. Contrary to requiring an unweighted p j T veto on jets, applying a kinematic selection using or binning events in τ j C for jets does not introduce extra logarithms (or minimizes their contribution) in the resummation region (low-p T ) of gluon fusion Higgs boson production cross section calculations including additional jets. Therefore, measuring the Higgs boson production cross section with respect to τ j C represents a test of QCD resummation [61]. In the γγj 1 j 2 system, the difference in φ and η between the photon pair and the two jets, and between the two jets themselves, are the quantities of interest, where the η difference between the photon pair and the jets is taken as the η difference between the photon pair and the average η of the two jets, defined as |η j 1 j 2 − η γγ | [62]. This observable serves as an additional probe of the properties of the VBF production mode. The invariant mass of the two leading-p T jets is also considered.
In the VBF-enriched phase space the cross section is measured with respect to the transverse momentum of the diphoton system and of the subleading jet, as well as the difference in φ between the dijet and the diphoton and the two leading jets themselves. Two double-differential measurements of the H → γγ cross section are performed: one with respect to the transverse momentum of the diphoton system and the number of jets, and the other one with respect to the diphoton transverse momentum and τ j C . Measurements are also performed with respect to the number of jets, the number of leptons, the number of b-tagged jets, as well as the missing transverse momentum in the H → γγ event.
Finally, the H → γγ cross section is measured in dedicated regions of the fiducial phase space designed to loosely target specific production modes. Their selection criteria require the presence of one extra lepton, and p miss T > 100 GeV (WH) or p miss T < 100 GeV (complementary to WH). For the third region, targeting the ttH production mode, events need to have at least one extra lepton and at least one b-tagged jet.

Statistical analysis
The events are categorized in bins of the observable of interest and in categories of the relative mass resolution estimator σ D m , as described in Section 6. To unfold the bins defined at reconstruction-level to the particle-level bins, the signal simulation sample is further split into bins in the particle-level observable of interest.
The signal production cross section is extracted through a simultaneous extended maximum likelihood fit of parametric signal and background probability density functions to the diphoton invariant mass spectrum of events falling into each reconstruction-level bin further divided into σ D m categories. The shape of the signal contribution is determined separately for each of the two-dimensional bins originating from dividing the events in reconstruction-level bins and particle-level bins, as described in Section 8.1. This makes the signal model sensitive to the fraction of events from each particle-level bin being reconstructed in a reconstruction-level bin and category, thus performing the maximum likelihood unfolding of the detector response. The background shape is derived through the discrete profiling method [63] in each reconstructionlevel bin further split in categories of σ D m , as described in Section 8.2. The imperfect alignment of the fiducial phase space at the particle level and the reconstruction selection means that there are events that do not enter the fiducial phase space, but are selected at reconstruction level and vice versa. This effect has to be taken into account by the analysis since the theoretical predictions for the cross sections that are to be measured are calculated at the particle level. This is done by deriving a dedicated probability density function (pdf) for signal events originating outside the fiducial phase space at particle level, and taking them into account in the maximum likelihood fit as an additional signal-like component.
For observables that are measured in restricted regions of the phase space, the cross section of the events that enter the baseline fiducial phase space, but not the restricted region, is left floating in the fit and treated as an additional bin. Letting the cross section in this additional bin float and taking migrations into and out of it into account allows for a tighter constraint of the cross sections within the restricted phase space region. The efficiency multiplied with the acceptance per particle-level bin, reconstruction-level bin, and resolution category combination is shown in Fig. 5 for the year 2018 and the observables p T and n jets .
The likelihood in the i-th σ D m category and j-th reconstruction-level observable bin can be written as: The σ D m categories also run over the three years of data taking, with three σ D m categories per year, as presented in Section 6. In Eq. (5), L i stands for the total analyzed integrated luminosity for the year with which the respective category is associated. The number of mass bins of the m γγ distribution is denoted with n m γ γ , and n b stands for the number of reconstruction-level bins for a given variable. The parameters of interest are encoded in the vector ∆ σ fid = (∆σ fid 1 , . . . , ∆σ fid n b ), which contains the per particle-level bin fiducial cross sections multiplied by the branching fraction of the Higgs boson decaying into two photons. The particle-level cross sections are directly obtained from the maximum likelihood fit and are not constrained to be positive, as   this would invalidate the method to extract the confidence intervals described below.
The events originating from a particle-level bin k are reconstructed in the reconstruction-level bin j, and category i. The magnitude of the migration is encoded in the detector response matrix K ij k . The pdfs in m γγ modelling the signal are denoted with S ij k for the particle-level bin k and reconstruction-level bin j and category i. The corresponding background pdf for reconstruction-level bin j and category i is written as B ij . The sum of the signal and background pdfs is used to fit the number of observed events n lij ev in m γγ bin l of reconstructionlevel bin j and category i. The numbers of signal and background events are denoted as n ij sig and n ij bkg , respectively. For signal events that enter the reconstruction-level analysis selection, while originating from outside the fiducial phase space, the label OOA (outside of acceptance) is used, where n ij OOA is the number of events that fulfil this condition and S ij OOA the pdf that describes their shape in m γγ . The number of events, n ij OOA , is fixed to the value predicted by the MADGRAPH5 aMC@NLO generator including the NNLOPS reweighting for the gluon fusion production mode. Finally, the nuisance parameters associated with the signal and background model are encoded as θ S and θ B , respectively.
The complete extended likelihood can be written as: with n cat being the number of σ D m categories across all years. Here, pdf stands for the probability density function of the nuisance parameters, while Pois denotes the Poisson probability distribution. From Eq. (6) the fiducial cross section ∆ σ fid can be extracted with a maximum likelihood fit. The mass of the Higgs boson m H is treated as a Gaussian constraint with a central value and width set according to the most precise CMS measurement of m H = 125.38 ± 0.14 GeV [46]. This means that effectively the mass values are independent per observable, but all particlelevel bins of one observable share the same value.
The negative log-likelihood ratio q(∆ σ) is used as the test statistic to obtain the uncertainties and correlation matrices [64]: where the number of background events and both sets of nuisance parameters are written as θ = (n bkg , θ S , θ B ); ∆ˆ σ andˆ θ denote the best fit estimates of the cross sections and nuisance parameters, whileˆ θ ∆ σ denotes the best fit estimates of the nuisance parameters at a fixed value of ∆ σ.

Signal model
The signal model is derived in a parametric form for each observable separately. The simulated m γγ signal shape at reconstruction level, including efficiency corrections estimated with the help of control samples, is fit with a sum of up to four Gaussian pdfs in each bin of the particle-level observable, reconstruction-level observable bins, and reconstruction categories. The shapes for right and wrong vertices are fit individually and added together to give the full signal pdf per bin. A signal shape is derived for each of the three mass hypotheses with m H ∈ {120, 125, 130} GeV. For mass values between these points, the signal pdf parameters are linearly interpolated. The sample for deriving the signal model is obtained by combining the four dominating SM Higgs boson production modes simulated with the MAD-GRAPH5 aMC@NLO generator, reweighted to the NNLOPS prediction for gluon fusion. Examples of signal pdfs used for the fiducial cross section measurement for different σ D m categories are shown in Fig. 4.

Background model
The determination of the background model uses a strategy called the discrete profiling method [63]. Since the exact form of the background is not known, a number of different functions to fit the smoothly falling shape of the m γγ background are tried, and their choice is treated as a discrete nuisance parameter in the maximum likelihood fit used to extract the differential cross section. The background fit is done in each reconstruction-level observable and reconstruction category bin over the range of 100 GeV < m γγ < 180 GeV. The families of functions considered for B m γγ | θ B are exponentials, power-law functions, Laurent polynomials, and polynomials in the Bernstein basis. A penalty term equal to the number of degrees of freedom is added to the double negative log-likelihood, hence allowing the fit to determine the number of degrees of freedom for all considered families.

Systematic uncertainties
In the following, the sources of systematic uncertainty that are considered in the analysis are described. Since the discrete profiling method is used to derive both the shape and the normalization of the background, the systematic uncertainty related to the background description is included in the final results as a component of the statistical uncertainty. Uncertainties that influence the shape of the signal model are implemented as nuisance parameters built into the signal model. Systematic variations that leave the shape of the m γγ distribution unaltered are taken into account as log-normal nuisances affecting only the predicted signal rate. If a systematic variation influences the event selection and thus changes the event yields in several categories, it is considered as a category migration systematic uncertainty, i.e. as a correlated log-normal uncertainty in the category yields. The systematic uncertainty arising from the imperfect modelling of the pileup distribution in the simulation is incorporated in to the uncertainty of the affected quantities (chiefly the photon MVA identification and vertex identification). Uncertainties of a theoretical nature are fully correlated among all years of data taking. The overall normalization of the particle-level bins is kept constant while evaluating the effect of the theoretical systematic uncertainties. This means that their effect is only considered as causing event migrations between reconstruction level bins and categories within the particular year of data taking. All yield and simple log-normal systematic variations of an experimental nature are treated as independent across years. The shape systematic uncertainties related to the scale and smearing corrections of the photon energy are fully correlated among years.
The theoretical systematic uncertainties include: • PDF uncertainty: this accounts for the imperfect knowledge of the parton distribution functions of the proton. These PDF uncertainties are derived through the relative yield variation per particle-level bin after reweighting the events with several alternative sets of weights. These weights are derived according to the prescriptions in PDF4LHC [65] and NNPDF3.1 [66] for 2017 and 2018, and NNPDF3.0 [67] with the MC2HESSIAN [68] procedure for 2016. This systematic uncertainty is treated as fully correlated between 2017 and 2018 and uncorrelated between 2016 and the other years of data taking. The largest effect on the cross section from reweighting across all sets of PDFs is at the per-mille level; • α S uncertainty: this uncertainty is related to the variation of the value of the strong force coupling constant used for the evaluation of the parton distribution functions, following the prescription in Ref. [65]. The impact of the nuisance parameter associated with this uncertainty in the cross section is at the per-mille level, similar to the PDF uncertainty; • Renormalization and factorization scales uncertainties: the migration uncertainty between analysis bins caused by the scale variations are estimated by moving both scales by a factor of 2, excluding the (2, 1/2) and (1/2, 2) combinations, as they are unphysical. The associated nuisance parameters can have an impact of up to 0.5% on the measured cross section.
The experimental uncertainties affecting the signal shape include: • Photon energy scale and resolution: this is the uncertainty related to the scale of the photon energy in data and its smearing correction in simulation [46]. The variations are evaluated using Z → ee events by varying the R 9 distribution and the selection criteria. The statistical uncertainty from the Z → ee sample used to derive the scale and smearing corrections is also taken into account. The uncertainty amounts to 0.05-0.15% for most of the photons, while they can go up to 3% in some categories. • Photon energy scale nonlinearity: the difference between the linearity of the energy scale in data and simulation is covered by this variation. It is estimated using Lorentzboosted Z → ee events and it amounts to 0.2%. • Cluster shape uncertainties: the energy regression used to correct the photon energy, as well as the data-to-simulation corrections of the photon identification MVA input variables, are all evaluated on electrons from Z boson decays. Therefore an uncertainty in the shower shape variables that accounts for the difference between electrons and photons has to be taken into account in an analysis that considers photons as signal. The resulting uncertainty in the photon energy scale is 0.01-0.15%, depending on the |η| and R 9 values of the photon [69]. • Nonuniformity of light collection: this uncertainty is related to the modelling of the light collection depending on the emission depth in the ECAL crystals. It amounts to 0.16-0.25% for photons with R 9 > 0.96. For low-R 9 photons, the uncertainty is below 0.07% [46]. • Modelling of the material upstream of ECAL: the showering behaviour of photons is influenced by the material upstream of the ECAL. This systematic variation is associated with the uncertainty in the tracker material model. Its impact on the photon energy scale is evaluated by varying the amount of material in the simulation by 5%. The uncertainty amounts to 0.02-0.05% for central photons, increasing to 0.25% for photons in the endcap [13]. • Vertex assignment: this variation covers the uncertainty in the fraction of events that have a reconstructed vertex that lies within 1 cm of the true vertex in data and simulation. The signal model includes a nuisance parameter that allows the number of events in which the vertex is at a distance of more than 1 cm from the true vertex to vary by ±2%. The corresponding nuisance parameter can lead to variations of the cross section by a few parts per mille.
The systematic uncertainties that only affect the event yield are: • Integrated luminosity: the uncertainty in the measured value of the integrated luminosity and its correlation across the years is applied following the CMS recommendations in Refs. [32][33][34]. Its magnitude is 1.2% for 2016, 2.3% for 2017, and 2.5% for 2018, while the overall uncertainty for the 2016-2018 period is 1.6%. • Per-photon energy resolution: the impact on the per-category event yields from the per-photon energy resolution mismodelling in simulation is evaluated by shifting the σ E /E value for each photon up and down by 2% in all events and repeating the full analysis. The magnitude of the per-event shift covers the data-to-simulation mismatch in the invariant mass resolution for Z → ee events. This uncertainty has one of the largest systematic effects on the measurement, which can be up to 2% in certain cross sections. • Jet energy scale and smearing corrections: the uncertainty in the jet energy scale is at the percent level and depends on the p T and η of the jet. It is propagated to the migrations between jet bins by varying the corrections within their uncertainties. The systematic uncertainties related to jet energy scale and smearing are partly correlated across years [70]. For observables involving jets, this uncertainty can impact the cross section measured as a function of those by up to 20%. • Trigger efficiency: the trigger efficiency is measured with Z → ee events using the tag-and-probe method. Its uncertainty affects bin migrations and has an effect of less than 1%. An additional source of uncertainty is included to cover the gradual shift of the inputs to the ECAL L1 trigger timing for values of |η| > 2 in 2016 and 2017. This affected mostly jets but also photons. The measured cross sections are barely affected by this systematic uncertainty, its effect being below the permille level. • Photon identification MVA score: the input variables to the photon identification MVA are corrected in simulation to match data as described in Section 4.2.3. This leads to a better agreement of the photon identification MVA score between simulation and data. The residual mismatch is covered by evaluating the uncertainty of the correction method and shifting the value of the photon identification BDT accordingly, as shown in Fig. 3. The uncertainty of the correction method is evaluated by varying the size of the data set used for its derivation. Finally, the uncertainty in the category yields is evaluated by propagating the variation of the photon identification MVA through the event selection. The nuisance parameters associated with this systematic variation have an effect of up to 1% on the measured cross section. Summing these effects in quadrature yields an impact of 1.5% on the inclusive fiducial cross section measurement. • Photon preselection: the uncertainty in the efficiency of the photon preselection is evaluated by varying the signal and background shapes that are used to compute the preselection efficiency. This is done in data and simulation using Z → ee tagand-probe events. The resulting variation is propagated to the corresponding scale factor. Its magnitude is below 1%. • Missing transverse momentum: the uncertainty in the missing transverse momentum is computed by shifting the scale and resolution of the particles included in the computation of the missing transverse momentum. This leads to yield variations for bins in observables involving missing transverse momentum. This systematic uncertainty affects only the measurement of the H → γγ cross section differential in p miss T . Its impact on the cross section is negligible.
• Pileup jet identification: the variation in the pileup jet identification score leads to bin migration in analysis bins related to jet variables. The resulting effect on the results presented here is negligible. • Lepton identification: the electron and muon identification efficiencies are varied within their uncertainties leading to bin migrations in lepton observables. These are propagated through the full analysis chain to calculate the related impact on the measured cross sections. This uncertainty only affects the measurement of the H → γγ cross section as a function of the number of leptons produced alongside the diphoton system. The magnitude of its effect is roughly 5%. • b tagging efficiency: the b jet tagging efficiency for the working point used in this analysis is varied within its uncertainty. This leads to bin migrations for observables involving b jets and in the phase spaces with a requirement on the number of b jets. This only affects the cross section measurement with respect to the number of b tagged jets. The effect of the corresponding nuisance parameter on the cross section is between 5 and 20%.
A summary of the dominant systematic impacts on the inclusive fiducial H → γγ cross section measurement are given in Table 4.

Results
In this section, the results of the analysis are presented. The cross section measured for the fiducial region defined in Table 2 is: The nominal theoretical predictions for all results presented in this section have been extracted from MADGRAPH5 aMC@NLO (version 2.6.5), reweighted to match the NNLOPS prediction for gluon fusion and interfaced with PYTHIA 8 (version 8.240) using the CP5 tune. Figure 6 shows the likelihood scan for the fiducial cross section and the nominal theoretical prediction from MADGRAPH5 aMC@NLO reweighted to match the NNLOPS prediction for gluon fusion. Figure 7 shows the m γγ histogram with the hypotheses from the signal and background models at the best fit point. The events from all categories are added, weighted by the The uncertainty in the prediction of the fiducial inclusive H → γγ cross section, shown as the hatched area in Fig. 6 and of the cross section in the dedicated phase space regions, shown as the hatched areas in Fig. 8, are the combination of several components. The dominating one of these is the uncertainty in the Higgs boson production cross section, with the second largest contribution coming from the uncertainty in the H → γγ branching fraction, both taken from [15]. The remaining contributions come from the variation of the fiducial acceptance due to varying the set of PDF replicas [66], the value of α S by 0.002 around its nominal values of 0.118 and finally the variation of renormalization and factorization scales by a factor of 2, while excluding the (2, 1/2) and (1/2, 2) variations. The interference between the H → γγ signal and the continuous diphoton background [71] has not been taken into account for the predictions presented here. The measured and predicted cross sections for selected special regions of the fiducial phase space are shown in Fig. 8. Two examples of correlation matrices, calculated using Eq. (7), are shown in Fig. 9. The correlation matrices for all other observables are available in the HEPData record for this analysis [29].
The results for the differential fiducial cross section, with respect to all observables listed in Table 3, are shown in Figs. 10-15. The results of the double-differential cross section measurement with respect to p T and n jets , and p T and τ  several additional predictions. The HX component, denoting the sum of the Higgs boson production cross sections from the VBF, VH, and ttH production processes, is taken from the MADGRAPH5 aMC@NLO simulation and is common to the different SM predictions shown. The predicted cross section for the gluon-fusion production mode is taken from the MAD-GRAPH5 aMC@NLO simulation, with and without NNLOPS reweighting, and the POWHEG event generator [72][73][74][75], and added to the HX component to obtain the total per-bin predictions shown. The uncertainty in the theoretical predictions only takes into account the variation in the predicted differential cross section shape coming from varying the set of PDF replicas, the renormalization and factorization scales, and α S . The uncertainty in the total Higgs boson production cross section and branching fraction is not taken into account for the results on the differential cross sections.
Overall, the differential cross section results agree within uncertainties with the nominal SM prediction. For each observable, a p-value is calculated using the test statistic given in Eq. (7) evaluated at the SM point, where the H → γγ cross section is set to the nominal SM value in all particle-level bins, extracted using the MADGRAPH5 aMC@NLO simulation with the NNLOPS reweighting for ggF. The p-value is then computed on the χ 2 pdf with the number of degrees of freedom set to the number of particle-level bins for the respective observable. The cross section measured for the underflow bin is not reported on the results figure if not specified otherwise, but always taken into account for the calculation of the p-value. The observed pvalues for the SM point are from 0.004 to 0.96, with the cross section measured as a function of p γγ T having a p-value of 0.24, and as a function of n jets a p-value of 0.69. For the fiducial cross section a p-value of 0.73 is observed for the SM point. The lowest p-value of 0.004 is seen for the difference between the η coordinate of the leading and subleading jet |∆η j 1 j 2 |, shown in Fig. 14. The per-bin uncertainties for observables of the diphoton system range from 10% to 40%. For observables that involve the leading-p T jet, uncertainties reach around 100%.
For the two-jet phase space for observables being calculated with the two leading-p T jets, uncertainties are at the same level, because of the larger bin sizes. The cross section measurements as a function of p γγ T , p j 1 T , |∆φ γγ,j 1 j 2 |, and |∆φ j 1 ,j 2 | are also performed in the VBF-enriched phase space, where the per-bin uncertainties can reach 150%. Overall, an under-fluctuation is observed for events that match the criteria for this phase space (see Fig. 8 Figure 9: The correlation matrices for the cross sections σ fid per particle-level bin for p γγ T (upper), and n jets (lower), as given in Table 3, extracted from the simultaneous maximum likelihood fit for the cross sections. For the blue lines no NNLOPS reweighting is applied and the green lines take the prediction for the ggH production mode from POWHEG. The hatched areas show the uncertainties in theoretical predictions on both the ggF and HX components. Only effects coming from varying the set of PDF replicas, the α S value, and the renormalization and factorization scales that impact the shape are taken into account here, the total cross section is kept constant at the value from Ref. [15]. The given p-values are calculated for the nominal SM prediction and the bottom panes show the ratio to the same prediction. If the last particle-level bin expands to infinity is explicitly marked on the plot together with the normalization of this bin.

Summary
The measurement of the fiducial inclusive Higgs boson (H) production cross section with the H → γγ decay mode has been presented. The fiducial phase space is defined by the ratio of the transverse momentum (p T ) of the leading (subleading) photon to diphoton invariant mass satisfying p T /m γγ > 1/3 (1/4), their pseudorapidity being within |η| < 2.5, and both photons being isolated. The production cross section for the Higgs boson decaying into two photons in the aforementioned phase space is measured to σ fid = 73.4 +6.1 −5.9 fb, in agreement with the theoretical prediction from the standard model (SM) of 75.4 ± 4.1 fb.
Furthermore, the pp → H + X, H → γγ cross section in the fiducial phase space has been measured as a function of observables of the diphoton system, as well as several others involving properties of the leading-p T and subleading-p T jets. Observables corresponding to the number of jets, leptons, and b-tagged jets are included as well. For the first time using the CMS detector, the cross section has been measured as a function of the rapidity weighted jet-p T (τ j C ), using up to six additional jets in the event. The cross section as a function of a measure for the deviation from "back-to-backness" |φ * η | for the diphoton system has been measured for the first time using the H → γγ channel. Two double-differential cross section measurements have been performed: one in bins of p T and the number of jets, the other in bins of p T and τ j C . A selected set of differential measurements has been performed in a dedicated phase space enriched with events compatible with vector boson fusion Higgs boson production. Finally, the production cross section has been measured in three fiducial phase spaces loosely targeting the vector boson and tt associated production modes.
Overall, the performed differential fiducial cross section measurements of the Higgs boson production in proton-proton collisions are found to be in agreement with the SM prediction within the uncertainties. [4] ATLAS and CMS Collaborations, "Measurements of the Higgs boson production and decay rates and constraints on its couplings from a combined ATLAS and CMS analysis of the LHC pp collision data at √ s = 7 and 8 TeV", JHEP 08 (2016) 045, doi:10.1007/JHEP08(2016)045, arXiv:1606.02266.
[21] ATLAS Collaboration, "Measurement of the properties of Higgs boson production at √ s = 13 TeV in the H → γγ channel using 139 fb −1 of pp collision data with the ATLAS experiment", 2022. arXiv:2207.00348. Submitted to JHEP.
[25] ATLAS Collaboration, "Measurements of the Higgs boson inclusive and differential fiducial cross-sections in the diphoton decay channel with pp collisions at √ s = 13 TeV with the ATLAS detector", 2022. arXiv:2202.00487. Submitted to JHEP.