Measurements of production cross sections of the Higgs boson in the four-lepton final state in proton–proton collisions at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{s} = 13\,\text {TeV} $$\end{document}s=13TeV

Production cross sections of the Higgs boson are measured in the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathrm{H}} \rightarrow {\mathrm{Z}} {\mathrm{Z}} \rightarrow 4\ell $$\end{document}H→ZZ→4ℓ (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\ell ={\mathrm{e}},{{{\upmu }}_{\mathrm{}}^{\mathrm{}}} $$\end{document}ℓ=e,μ) decay channel. A data sample of proton–proton collisions at a center-of-mass energy of 13\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\,\text {TeV}$$\end{document}TeV, collected by the CMS detector at the LHC and corresponding to an integrated luminosity of 137\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\,\text {fb}^{-1}$$\end{document}fb-1 is used. The signal strength modifier \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu $$\end{document}μ, defined as the ratio of the Higgs boson production rate in the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$4\ell $$\end{document}4ℓ channel to the standard model (SM) expectation, is measured to be \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu =0.94 \pm 0.07 \,\text {(stat)} ^{+0.09}_{-0.08} \,\text {(syst)} $$\end{document}μ=0.94±0.07(stat)-0.08+0.09(syst) at a fixed value of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m_{{\mathrm{H}}} = 125.38\,\text {GeV} $$\end{document}mH=125.38GeV. The signal strength modifiers for the individual Higgs boson production modes are also reported. The inclusive fiducial cross section for the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathrm{H}} \rightarrow 4\ell $$\end{document}H→4ℓ process is measured to be \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2.84^{+0.23}_{-0.22} \,\text {(stat)} ^{+0.26}_{-0.21} \,\text {(syst)} \,\text {fb} $$\end{document}2.84-0.22+0.23(stat)-0.21+0.26(syst)fb, which is compatible with the SM prediction of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2.84 \pm 0.15 \,\text {fb} $$\end{document}2.84±0.15fb for the same fiducial region. Differential cross sections as a function of the transverse momentum and rapidity of the Higgs boson, the number of associated jets, and the transverse momentum of the leading associated jet are measured. A new set of cross section measurements in mutually exclusive categories targeted to identify production mechanisms and kinematical features of the events is presented. The results are in agreement with the SM predictions.


Introduction
The discovery of the Higgs boson (H) in 2012 by the ATLAS and CMS collaborations [1][2][3] has been a major step towards the understanding of the electroweak symmetry breaking mechanism [4][5][6][7][8][9]. Further studies by the two experiments [10][11][12][13] have shown that the properties of the new particle are consistent with the standard model (SM) expectations for the H boson.
The H → ZZ → 4 decay channel ( = e, μ) has a large signal-to-background ratio thanks to a low background rate and the complete reconstruction of the final state decay e-mail: cms-publication-committee-chair@cern.ch products, capitalizing on the excellent lepton momentum resolution of the CMS detector. The measurements performed using this decay channel with the LHC Run 1 data set at center-of-mass energies of 7 and 8 TeV, and the Run 2 data set at 13 TeV include the determination of the mass, the spin and the parity of the H boson [14][15][16][17][18][19], its width [20-23], the inclusive and differential fiducial cross sections [18,[24][25][26][27][28], and the tensor structure of the H boson interaction with a pair of neutral gauge bosons in both on-shell and off-shell regions [17,19,21,29,30].
This paper presents the measurement of production cross sections in granular kinematic regions of the H boson in the H → ZZ → 4 decay channel. A data sample of proton-proton (pp) collisions at a center-of-mass energy of √ s = 13 TeV, collected by the CMS detector at the LHC and corresponding to an integrated luminosity of 137 fb −1 is used. The inclusive signal strength modifier, defined as the ratio of the H boson production rate in the 4 channel to the SM expectation, and signal strength modifiers for the individual H boson production modes are measured. The measurements of the inclusive and differential fiducial cross sections are also presented, and their compatibility with the SM predictions is tested. The present analysis is similar to that previously performed by the CMS Collaboration [18], but is based on a larger data sample.
In addition, measurements of the H boson cross sections within the simplified template cross section (STXS) framework [31][32][33] are also presented. The main goals of the STXS framework are to increase the reinterpretability of the precision H boson measurements and to minimize the theory dependence. This is achieved by defining exclusive kinematic regions in the H boson production phase space. 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 beyond the SM (BSM) scenarios. These kinematic regions, referred to as bins, are defined in different stages corresponding to increasing degrees of granularity. This paper presents results in the STXS stage 0 where the bins correspond closely to the different H boson production mechanisms. Previous measurements of cross sections in stage 0 production bins in the H → 4 decay channel were already presented by the CMS Collaboration [18]. In the STXS framework, additional stages are defined by further splitting of the bins enhancing the sensitivity to possible signature of BSM physics at high transverse momentum of the H boson. Measurements of stage 0, stage 1, and stage 1.1 cross sections in the H → 4 decay channel were recently published by the ATLAS Collaboration [27]. The most recent refinement of STXS binning is referred to as STXS stage 1.2. This paper presents a first set of the cross section measurements in the STXS stage 1.2 bins in the H → 4 decay channel.
The paper is organized as follows. A brief introduction of the CMS detector is given in Sect. 2. The data, as well as the simulated signal and background samples, are described in Sect. 3. The event reconstruction and selection, the kinematic discriminants, and the categorization of the H boson candidate events are described in Sects. 4, 5, and 6 , respectively. The background estimation is detailed in Sect. 7 while the signal modeling is described in Sect. 8. The experimental and theoretical systematic uncertainties are described in Sect. 9 and the results are presented in Sect. 10. Concluding remarks are given in Sect. 11.

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. Forward calorimeters extend the pseudorapidity η coverage provided by the barrel and endcap detectors. Muons are detected in gas-ionization chambers embedded in the steel flux-return yoke outside the solenoid.
Events of interest are selected using a two-tiered trigger system. 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 [34]. The second level, known as the high-level trigger, consists of a farm of processors running a version of the full event reconstruction software optimized for fast processing, and reduces the event rate to around 1 kHz before data storage [35].
The candidate vertex with the largest value of summed physics-object squared transverse momentum p 2 T is taken to be the primary pp interaction vertex (PV). The physics objects are the jets, clustered using the jet finding algorithm [36,37] with the tracks assigned to candidate vertices as inputs, and the associated missing transverse momentum, taken as the negative vector sum of the p T of those jets.
The electron momentum is estimated by combining the energy measurement in the ECAL with the momentum measurement in the tracker. The momentum resolution for electrons with p T ≈ 45 GeV from Z → ee decays ranges from 1.7% to 4.5%. It is generally better in the barrel region than in the endcaps, and also depends on the bremsstrahlung energy emitted by the electron as it traverses the material in front of the ECAL [38]. The ECAL consists of 75 848 lead tungstate crystals, which provide coverage of |η| < 1.48 in the barrel region and 1.48 < |η| < 3.00 in the two endcap regions (EE). Preshower detectors consisting of two planes of silicon sensors interleaved with a total of 3X 0 of lead are located in front of each EE detector.
Muons are measured in the pseudorapidity range |η| < 2.4, with detection planes made using three technologies: drift tubes, cathode strip chambers, and resistive plate chambers. The single muon trigger efficiency exceeds 90% over the full η range, and the efficiency to reconstruct and identify muons is greater than 96%. Matching muons to tracks measured in the silicon tracker results in a relative transverse momentum resolution, for muons with p T up to 100 GeV, of 1% in the barrel and 3% in the endcaps. The p T resolution in the barrel is better than 7% for muons with p T up to 1 TeV [39].
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. [40].

Data and simulated samples
This analysis is based on the pp collision data collected by the CMS detector at the LHC in 2016, 2017, and 2018 with integrated luminosities of 35.9, 41.5, and 59.7 fb −1 , respectively [41][42][43]. The collision events are selected by highlevel trigger algorithms that require the presence of leptons passing loose identification and isolation requirements. The main triggers select either a pair of electrons or muons, or an electron and a muon. The minimal transverse momentum of the leading and subleading leptons changed throughout the years to account for the different data-taking conditions and is summarized in Table 1.
To maximize the coverage of the H → 4 phase space, triggers requiring three leptons with relaxed transverse momenta thresholds and no isolation requirement are also used, as are isolated single-electron and single-muon triggers. The overall trigger efficiency for simulated signal events that pass the full event selection (described in Sect. 4) is larger than 99%. The trigger efficiency is derived from data using a sample of 4 events collected by the singlelepton triggers and a method based on the "tag and probe" Table 1 The minimal p T of the leading/subleading leptons for the main di-electron (e/e), di-muon (μ/μ), and electron-muon (e/μ, μ/e) highlevel trigger algorithms used in the H → 4 analysis in 2016, 2017, and 2018 e/e (GeV) μ/μ (GeV) e/μ, μ/e (GeV) technique. One of the four leptons is matched to a candidate reconstructed by the single-lepton trigger and the remaining three leptons in the event are used as probes. The probe leptons are combined in an attempt to reconstruct any of the triggers used in the analysis. The efficiency in data is found to be in agreement with the expectation from the simulation. Monte Carlo (MC) simulation samples for the signals and the relevant background processes are used to evaluate the signal shape, estimate backgrounds, optimize the event selection, and evaluate the acceptance and systematic uncertainties. The SM H boson signals are simulated at next-toleading order (NLO) in perturbative QCD (pQCD) with the powheg 2.0 [44][45][46] generator for the five main production processes: gluon fusion (ggH) [47], vector boson fusion (VBF) [48], associated production with a vector boson (VH, where V is a W or a Z boson) [49], and associated production with a pair of top quarks (tt H) [50]. The ZH production occurs in two ways, qq → ZH and a much smaller contribution from gg → ZH, which is simulated at leading order (LO) using jhugen 7.3.0 [51][52][53][54][55]. In addition to the five main production processes, the contributions due to H boson production in association with a single top quark (tH) and either a quark (tHq) or a W boson (tHW) are simulated at LO using jhugen 7.0.2 and MadGraph5_amc@nlo 2.2.2 [56], respectively. The associated production with a pair of bottom quarks (bbH) is simulated at LO with jhugen 7.0.2. In all cases, the decay of the H boson to four leptons is modeled with jhugen 7.0.2. The theoretical predictions used for the various production and decay modes can be found in Refs.  and are summarized in Ref. [32].
The ZZ background contribution from quark-antiquark annihilation is simulated at NLO pQCD with powheg 2.0 [80], while the gg → ZZ process is generated at LO with mcfm 7.0.1 [81]. The WZ background and the triboson backgrounds ZZZ, WZZ, and WWZ are modeled at NLO using MadGraph5_amc@nlo 2.4.2. The smaller tt Z, tt WW, and tt ZZ background processes are simulated at LO with Mad-Graph5_amc@nlo 2.4.2. The events containing Z bosons with associated jets (Z+jets) are simulated at NLO with MadGraph5_amc@nlo 2.4.2 and the tt background is simulated at NNLO with powheg 2.0. The reducible background determination does not rely on the MC but is based on data, as described in Sect. 7.2.
All signal and background event generators are interfaced with pythia 8.230 [82] using the CUETP8M1 tune [83] for the 2016 data-taking period and the CP5 tune [84] for the 2017 and 2018 data-taking periods, to simulate the multiparton interaction and hadronization effects. The NNPDF3.0 set of parton distribution functions (PDFs) is used [85]. The generated events are processed through a detailed simulation of the CMS detector based on Geant4 [86,87] and are reconstructed with the same algorithms that are used for data. The simulated events include overlapping pp interactions (pileup) and have been reweighted so that the distribution of the number of interactions per LHC bunch crossing in simulation matches that observed in data.

Event reconstruction and selection
The particle-flow (PF) algorithm [88] aims to reconstruct and identify each individual particle (PF candidate) in an event, with an optimized 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 PV 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 the response function of the calorimeters to hadronic showers. Finally, the energy of neutral hadrons is obtained from the corresponding ECAL and HCAL energies.
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 [89]. The p miss T is modified to account for corrections to the energy scale of the reconstructed jets in the event.
Muons with p μ T > 5 GeV are reconstructed within the geometrical acceptance, corresponding to the region |η μ | < 2.4, by combining information from the silicon tracker and the muon system [39]. The matching between the inner and outer tracks proceeds either outside-in, starting from a track in the muon system, or inside-out, starting from a track in the silicon tracker. Inner tracks that match segments in only one or two stations of the muon system are also considered because they may belong to very lowp T muons that do not have sufficient energy to penetrate the entire muon system. The muons are selected among the reconstructed muon track candidates by applying minimal requirements on the track in both the muon system and the inner tracker system, and taking into account the compatibility with small energy deposits in the calorimeters.
To discriminate between prompt muons from Z boson decay and those arising from electroweak (EW) decays of hadrons within jets, an isolation requirement of I μ < 0.35 is imposed, where the relative isolation is defined as In Eq. (1), p charged T is the scalar sum of the transverse momenta of charged hadrons originating from the chosen PV of the event. The quantities p neutral T and p γ T are the scalar sums of the transverse momenta for neutral hadrons and photons, respectively. The isolation sums involved are all restricted to a volume bound by a cone of angular radius ΔR = 0.3 around the muon direction at the PV, where the angular distance between two particles i and j is Since the isolation variable is particularly sensitive to energy deposits from pileup interactions, a p μ,PU T contribution is subtracted, defined as p T , where i runs over the charged hadron PF candidates not originating from the PV, and the factor of 0.5 corrects for the different fraction of charged and neutral particles in the cone [90].
Electrons with p e T > 7 GeV are reconstructed within the geometrical acceptance, corresponding to the pseudorapidity region |η e | < 2.5 [38]. Electrons are identified using a multivariate discriminant which includes observables sensitive to the presence of bremsstrahlung along the electron trajectory, the geometrical and momentum-energy matching between the electron trajectory and the associated cluster in the ECAL, the shape of the electromagnetic shower in the ECAL, and variables that discriminate against electrons originating from photon conversions. Instead of an additional isolation requirement, similar to the one for muons, the electron multivariate discriminant also includes the isolation sums described above ( p charged T , p neutral T , and p γ T ) but computed around the electron direction. The inclusion of isolation sums helps suppressing electrons originating from electroweak decays of hadrons within jets [91] and has a better performance than a simple requirement on the relative isolation observable. The package xgboost [92] is used for the training and optimization of the multivariate discriminant employed for electron identification and isolation. The training is performed with simulated events that are not used at any other stage of the analysis. Events are divided into six regions defined by two transverse momentum ranges (7-10 GeV and > 10 GeV) and three pseudorapidity regions: central barrel (|η e | < 0.8), outer barrel (0.8 < |η e | < 1.479), and endcaps (1.479 < |η e | < 2.5). Separate training is performed for the three different data-taking periods and selection requirements are determined such that the signal efficiency remains the same for all three periods.
The effect of the final-state radiation (FSR) from leptons is recovered as follows. Bremsstrahlung photons already associated to electrons in the reconstruction step are not considered in this procedure. Photons reconstructed by the PF algorithm within |η γ | < 2.4 are considered as FSR candidates if they satisfy the conditions p γ T > 2 GeV and I γ < 1.8, where the photon relative isolation I γ is defined as for the muon in Eq. (1). Every such photon is associated to the closest selected lepton in the event. Photons that do not satisfy the requirements ΔR(γ, )/(p γ T ) 2 < 0.012 GeV −2 and ΔR(γ, ) < 0.5 are discarded. The lowest-ΔR(γ, )/(p γ T ) 2 photon candidate of every lepton, if any, is retained. The photons thus identified are excluded from the isolation computation of the muons selected in the event.
In order to suppress muons from in-flight decays of hadrons and electrons from photon conversions, leptons are rejected if the ratio of their impact parameter in three dimensions, computed with respect to the PV position, to their uncertainty is greater or equal to four.
The momentum scale and resolution of electrons and muons are calibrated in bins of p T and η using the decay products of known dilepton resonances as described in Refs. [38,39].
A "tag and probe" technique [93] based on samples of Z boson events in data and simulation is used to measure the efficiency of the reconstruction and selection for prompt electrons and muons in several bins of p T and η . The difference in the efficiencies measured in simulation and data is used to rescale the yields of selected events in the simulated samples.
For each event, hadronic jets are clustered from the reconstructed particles using the infrared-and collinear-safe antik T algorithm [36,37] with a distance parameter of 0.4. The jet momentum is determined as the vectorial sum of all particle momenta in the jet, and is found from simulation to be within 5-10% of the true momentum over the whole p T spectrum and detector acceptance. Additional pp interactions within the same or nearby bunch crossings can contribute extra tracks and calorimetric energy depositions to the jet. To mitigate this effect, tracks identified as originating from pileup vertices are discarded and an offset correction is applied to correct for the remaining contributions. Jet energy corrections are derived from simulation to match 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 jet energy scale in data and simulation [94]. Jet energies in simulation are smeared to match the resolution in data. The jet energy resolution amounts typically to 16% at 30 GeV, 8% at 100 GeV, and 4% at 1 TeV. Additional selection criteria are applied to remove jets potentially dominated by anomalous contributions from various subdetector components or reconstruction failures. To be considered in the analysis, jets must satisfy the conditions p jet T > 30 GeV and |η jet | < 4.7, and be separated from all selected lepton candidates and any selected FSR photons by ΔR( /γ, jet) > 0.4. Jets are also required to pass the tight identification criteria and the tight working point of pileup jet identification described in Ref. [90].
For event categorization, jets are tagged as b jets using the DeepCSV algorithm [95], which combines information about impact parameter significance, secondary vertex, and jet kinematics. Data to simulation scale factors for the b tagging efficiency are applied as a function of jet p T , η, and flavor.
The event selection is designed to extract signal candidates from events containing at least four well-identified and isolated leptons, each originating from the PV and possibly accompanied by an FSR photon candidate. In what follows, unless otherwise stated, FSR photons are included in invariant mass computations.
First, Z candidates are formed with pairs of leptons of the same flavor and opposite-charge (e + e − , μ + μ − ) that pass the requirement 12 < m + − < 120 GeV. They are then combined into ZZ candidates, wherein we denote as Z 1 the Z candidate with an invariant mass closest to the nominal Z boson mass [96], and as Z 2 the other one. The flavors of the involved leptons define three mutually exclusive subchannels: 4e, 4μ, and 2e2μ.
To be considered for the analysis, ZZ candidates have to pass a set of kinematic requirements that improve the sensitivity to H boson decays. The Z 1 invariant mass must be larger than 40 GeV. All leptons must be separated in angular space by at least ΔR( i , j ) > 0.02. At least two leptons are required to have p T > 10 GeV and at least one is required to have p T > 20 GeV. In the 4μ and 4e subchannels, where an alternative Z a Z b candidate can be built out of the same four leptons, we discard candidates with m Z b < 12 GeV if Z a is closer to the nominal Z boson mass than Z 1 is. This rejects events that contain an on-shell Z and a low-mass dilepton resonance. To further suppress events with leptons originating from hadron decays in jet fragmentation or from the decay of low-mass resonances, all four opposite-charge lepton pairs that can be built with the four leptons (irrespective of flavor) are required to satisfy the condition m + − > 4 GeV, where selected FSR photons are disregarded in the invariant mass computation. Finally, the four-lepton invariant mass m 4 must be larger than 70 GeV, which defines the mass range of interest for the subsequent steps of the analysis.
In events where more than one ZZ candidate passes the above selection, the candidate with the highest value of D kin bkg (defined in Sect. 5) is retained, except if two candidates consist of the same four leptons, in which case the candidate with the Z 1 mass closest to the nominal Z boson mass is retained.

Kinematic discriminants
The full kinematic information from each event using either the H boson decay products and/or the associated particles in the H boson production is extracted by means of matrix element calculations and is used to form several kinematic discriminants. These computations rely on the MELA package [51][52][53]55] and exploit the jhugen matrix elements for the signal and the mcfm matrix elements for the background. Both the H boson decay kinematics and the kinematics of the associated production of H + 1 jet, H + 2 jets, VBF, ZH, and WH are explored. The full event kinematics is described by decay observables Ω H→4 or observables describing the associated production Ω H+jj , which may or may not include the H → 4 decay kinematic information depending on the use case. The definition of these observables can be found in Refs. [51][52][53].
Two types of kinematic discriminants are exploited in the H → 4 analysis. First we construct the three categorization discriminants in order to classify signal events into exclusive categories as defined in Sect. 6.2. Categorization discriminants are designed to increase the purity of the targeted production mechanism in a dedicated event category. In addition, we define another set of three kinematic discriminants that are taken as an observable in the two-dimensional likelihood fits carried out to extract the results, as described in Sect. 10. These kinematic discriminants are designed to separate the targeted H boson production mechanism from its dominant background.
Categorization discriminants are calculated following the prescription in Refs. [18,21,97]. The discriminants sensitive to the VBF signal topology with two associated jets, the VBF signal topology with one associated jet, and the VH (either ZH or WH) signal topology with two associated jets are: where P VBF , P Hjj , P Hj , and P VH are the probabilities for the VBF process, the ggH process in association with two jets (combination of gg/qg/qq parton collisions producing H + 2 jets), the ggH process in association with one jet (H+1 jet), and the VH process, respectively. The quantity dη j P VBF is the integral of the two-jet VBF matrix element probability over the η j values of the unobserved jet, with the constraint that the total transverse momentum of the H + 2 jets system is zero. The discriminant D VH 2jet , used for event categorization, is defined as the maximum value of the two discriminants, D VH 2jet = max(D ZH 2jet , D WH 2jet ). A set of three discriminants used in the likelihood fits is calculated as in Refs. [17,18]. The discriminant sensitive to the gg/qq → 4 process exploits the kinematics of the fourlepton decay system. It is used in most of the event categories described in Sect. 6 to separate signal from background and is defined as: where P gg sig is the probability for the signal and P qq bkg is the probability for the dominant qq → 4 background process, calculated using the LO matrix elements. In the VBF-2jettagged and VH-hadronic-tagged event categories (defined in Sect. 6.2), the background includes the QCD production of ZZ/Zγ * /γ * γ * → 4 in association with two jets, the EW background from the vector boson scattering (VBS), as well as the triboson (VVV) production process. We therefore use dedicated production-dependent discriminants based on the kinematics of the four-lepton decay and information from the associated jets (noted with VBF+ dec or VH+ dec), defined as: where P EW sig is the probability for the VBF and VH signal, P EW bkg is the probability for the VBS and VVV background processes, and P QCD bkg is the probability for ZZ/Zγ * /γ * γ * → 4 QCD production in association with two jets. The quantity c p (m 4 ) for category p is the m 4 -dependent parameter that allows to change the relative normalization of the EW probabilities, separately for the VBF and VH topologies. For each slice of m 4 , the distributions of the signal and back-ground discriminants are plotted, and the c p (m 4 ) value is determined in such a way that the two distributions cross at 0.5. This procedure allows rescaling of the distributions for the linear-scale binning of the templates used in the likelihood fits described in Sect. 10.

Event categorization
In order to improve the sensitivity to the H boson production mechanisms, the selected events are classified into mutually exclusive categories based on the features of the reconstructed objects associated with the H → 4 candidates. Event categorization is organized in two steps with increasing granularity of the categories. First step is primarily designed to separate the ggH, VBF, VH, and tt H processes. There is little sensitivity to bbH or tH, even though these production modes are considered explicitly in the analysis. The reconstructed event categories from the first step are further subdivided (as discussed in Sect. 6.2) in order to study each production mechanism in more detail. This subdivision is carried out by matching the recommended binning of the framework of STXS described in the following section.

STXS production bins
The STXS framework has been adopted by the LHC experiments as a common framework for studies of the H boson. It has been developed to define fine-grained measurements of the H boson production modes in various kinematic regions, and to reduce the theoretical uncertainties that are folded into the measurements. It also allows for the use of advanced categorization techniques and provides a common scheme for combining measurements in different decay channels or by different experiments. The regions of phase space defined by this framework are referred to as production bins and are determined by using generator-level information for H bosons with rapidity |y H | < 2.5. Generator-level jets are defined as anti-k T jets with a distance parameter of 0.4 and a p T threshold of 30 GeV; no requirement is placed on the generator-level leptons.
The STXS framework has been designed to complement the Run I measurements of the production signal strength modifiers and fiducial differential cross sections of the H boson by combining their advantages. The sensitivity to theoretical uncertainties in the signal strength modifier results is suppressed by excluding dominant theoretical uncertainties causing production bin migration effects from the STXS measurements. They are included only when comparing the results with the theoretical predictions. In contrast to the fiducial differential cross section measurements, in the STXS framework measurements are optimized for sensitivity by means of event categories and matrix element discriminants.
To account for the evolving experimental sensitivity, different stages of production bins with increasing granularity are developed.
The stage 0 production bins are called ggH, qqH, VH-lep, and ttH and are designed to closely match the main H boson production mechanisms. The qqH bin includes the EW production of the H boson in association with two quarks from either VBF or VH events with hadronic decays of the vector boson V. The VH-lep production bin includes VH events with leptonic decays of the vector boson V. The low rate bbH and tH production processes are considered together with the ggH and ttH production bins, respectively. In this analysis, a modified version of the stage 0 production bins is also studied, where instead of VH-lep and qqH bins we define the WH, ZH, and VBF bins that map the H boson production mechanisms without the splitting of the VH events in leptonic and hadronic decays.
Stage 1 of the STXS framework was designed by further splitting the bins from the stage 0, one of the main motives being the enhanced sensitivity to possible signatures of BSM physics. This is achieved by dividing stage 0 bins with additional requirements on generator-level quantities that include the transverse momentum of the H boson ( p H T ), the number of associated jets (N j ), the dijet invariant mass (m jj ), the transverse momentum of the H boson and the leading jet ( p Hj T ), and the transverse momentum of the H boson and the two leading jets ( p Hjj T ). These bins were designed in order to maximize sensitivity to new physics while also taking into account the current experimental sensitivity limited mostly by the amount of collected data. The most recent set of bins defined in the STXS framework is referred to as stage 1.2. This paper presents a first set of cross section measurements in the H → 4 channel for the stage 1.2 of the STXS framework. However, several stage 1.2 production bins are merged as the full set of bins cannot be measured with the current data sample. The merging scheme results in 19 production bins; it is illustrated in Fig. 1 and discussed in more detail below.
The ggH production process is split into events with p H T < 200 GeV and p H T > 200 GeV. The events with p H T > 200 GeV are placed into one single production bin called ggH/p T > 200. The events with p H T < 200 GeV are split in events with zero, one, and two or more jets. The events with zero or one jets are split into the following production bins according to the H boson p T : ggH-0j/p T [0, 10], ggH-0j/p T [10, 200], ggH-1j/p T [0, 60], ggH-1j/ p T [60,120], and ggH-1j/p T [120,200]. The events with two or more jets are split according to the dijet invariant mass as follows. The events with m jj < 350 GeV are split into three production bins according to the H boson p T : ggH-2j/p T [0, 60], ggH-2j/p T [60,120], and ggH-2j/p T [120,200]. The events with m jj > 350 GeV are all placed into one production bin ggH-2j/m jj > 350, which merges four bins originally suggested in stage 1.2 of the STXS framework.
The merging scheme of the electroweak qqH production bins is as follows. The events with zero jets, one jet, or with two or more jets with m j j < 60 GeV or 120 < m jj < 350 GeV correspond to production bins with insufficient statistics; they are all merged into one bin called qqH-rest. The events with two or more jets and 60 < m jj < 120 GeV are placed in the qqH-2j/m jj [60,120]  The three production processes qq → WH, gg → ZH, and qq → ZH are combined to build VH-lep reduced stage 1.2 production bins. Several proposed production bins are merged into just two bins according to the p T of the H boson: In stage 1.2 of the STXS framework the ttH stage 0 production bin is split in five different bins according to the p T of the H boson. Because of the very low expected yields all these bins are merged into a single bin that includes the tH production process as well.
Finally, in stage 1.2 the bbH production process, which has small cross section, is classified in the ggH-0j/p T [10, 200] production bin.
The first measurement of STXS stage 1.2 cross sections was recently performed by the CMS Collaboration [98].

Reconstructed event categories
In order to be sensitive to different production bins, the ZZ candidates that pass the event selection described in Sect. 4 are classified into several dedicated reconstructed event categories. The category definitions are based on the multiplicity of jets, b-tagged jets, and additional leptons in the event. Additional leptons are not involved in the ZZ candidate selection but, if present, should satisfy the identification, vertex compatibility, and isolation requirements. Requirements on the categorization discriminants described in Sect. 5, the invariant mass of the two leading jets, and the transverse momentum of the ZZ candidate are also exploited.
The event categorization is carried out in two steps. In the first step, the ZZ candidates are split into seven initial categories to target the main H boson production mechanisms corresponding to the stage 0 production bins. The first step of the categorization closely follows the analysis strategy from Fig. 1 Binning of the gluon fusion production process, the electroweak production process (combines VBF and VH with hadronic V decay), the VH production process with leptonic V decay (combining WH, ZH, and gluon fusion ZH production), and the tt H production process in the merged stage 1.2 of the STXS framework used in the H → ZZ → 4 analysis the previous publication [18]. To ensure exclusive categories, an event is considered for the subsequent category only if it does not satisfy the requirements of the previous one.
In the first categorization step, the following criteria are applied: -The VBF-2jet-tagged category requires exactly 4 leptons.
In addition there must be either 2 or 3 jets of which at most 1 is b-tagged, or at least 4 jets and no b-tagged jets. Finally, D VBF 2jet > 0.5 is required. -The VH-hadronic-tagged category requires exactly 4 leptons. In addition there must be 2 or 3 jets with no btagging requirements, or at least 4 jets and no b-tagged jets. Finally, D VH 2jet > 0.5 is required. -The VH-leptonic-tagged category requires no more than 3 jets and no b-tagged jets in the event, and exactly 1 additional lepton or 1 additional pair of opposite sign, same flavor leptons. This category also includes events with no jets and at least 1 additional lepton.
-The tt H-hadronic-tagged category requires at least 4 jets, of which at least 1 is b-tagged, and no additional leptons. -Thett H-leptonic-tagged category requires at least 1 additional lepton in the event. -The VBF-1jet-tagged category requires exactly 4 leptons, exactly 1 jet and D VBF 1jet > 0.7. -The untagged category consists of the remaining events.
Reconstructed events are further subdivided in the second step of the categorization that is designed to closely match the merged stage 1.2 production bins described in the previous section. In the second categorization step, the untagged, VBF-2jet-tagged, VH-hadronic-tagged, and VHleptonic-tagged categories are further split exploiting additional variables like the invariant mass of the two leading jets and the transverse momentum of the ZZ candidate. A total number of twenty-two reconstructed event categories is defined and details of the categorization are presented in Table 2. The irreducible background to the H boson signal in the 4 channel, which comes from the production of ZZ via qq annihilation or gluon fusion, is estimated using simulation. The fully differential cross section for the qq → ZZ process is computed at NNLO [99], and the NNLO/NLO K factor as a function of m ZZ is applied to the powheg sample. This K factor varies from 1.0 to 1.2 and is 1.1 at m ZZ = 125 GeV. Additional NLO electroweak corrections that depend on the initial state quark flavor and kinematics are also applied in the region m ZZ > 2m Z following the prescription in Ref. [100]. The production of ZZ via gluon fusion contributes at NNLO in pQCD. It has been shown that the soft collinear approximation is able to describe the cross section for this process and the interference term at NNLO [101]. Further calculations also show that the K factors are very similar at NLO for signal and background [102] and at NNLO for signal and interference terms [103]. Therefore, the same K factor is used for signal and background [104]. The NNLO K factor for the signal is obtained as a function of m ZZ using the hnnlo v2 program [105][106][107] by calculating the NNLO and LO gg → H → 2 2 cross sections for the H boson decay width of 4.07 MeV and taking their ratios. The NNLO/LO K factor for gg → ZZ varies from ≈2.0 to 2.6 and is 2.27 at m ZZ = 125 GeV; a systematic uncertainty of 10% is assigned to it when applied to the background process.
The triboson background processes ZZZ, WZZ, and WWZ, as well as tt Z, tt WW, and tt ZZ are also considered. These rare backgrounds are all estimated from simulation and are jointly referred to as the EW backgrounds.
Simulated samples are used to obtain shapes of the fourlepton invariant mass that are later used to build the likelihood function. For each irreducible background contribution, events are divided in three final states (4μ, 4e, and 2e2μ) and 22 event sub-categories defined in Sect. 6.1. To extract the shape of the m 4 distribution, expected yields are fitted to empirical functional forms built from a third order Bernstein polynomial. In sub-categories with not enough statistics to perform a fit, the shape is extracted from the inclusive distribution in the corresponding final state.

Reducible backgrounds
Additional backgrounds to the H boson signal in the 4 channel arise from processes in which decays of heavy-flavor hadrons, in-flight decays of light mesons within jets, or (for electrons) charged hadrons overlapping with π 0 decays are misidentified as leptons. The main processes leading to these backgrounds are Z+jets, tt +jets, Zγ +jets, WW+jets, and WZ+jets production. We denote these reducible backgrounds as "Z+X" since they are dominated by the Z+jets process. The contribution from the reducible background is estimated with two independent methods, each with dedicated control regions in data. The control regions are defined by the presence of both a lepton pair satisfying all the requirements of a Z 1 candidate and two additional opposite sign (OS) or same sign (SS) leptons; the two additional leptons satisfy identification requirements looser than those used in the analysis. These four leptons are then required to pass the analysis ZZ candidate selection. The event yield in the signal region is obtained by weighting the control region events by the lepton misidentification probability f e ( f μ ), defined as the fraction of non-signal electrons (muons) that are identified by the analysis selection criteria. A detailed description of both methods can be found in Ref. [18].
The lepton misidentification rates f e and f μ are measured as a function of p T and η by means of a sample that includes a Z 1 candidate consisting of a pair of leptons, both passing the selection requirements used in the analysis, and exactly one additional lepton passing the relaxed selection. Furthermore, the p miss T is required to be less than 25 GeV in order to suppress contamination from WZ and tt processes.
For the OS method, the mass of the Z 1 candidate is required to satisfy the condition |Z 1 − m Z | < 7 GeV in order to reduce the contribution of (asymmetric) photon conversions, which is estimated separately. In the SS method, the contribution of photon conversions to the misidentification rate is estimated with dedicated samples.
The predicted yields of the reducible background from the two methods are in agreement within their uncertainties for each final state (4μ, 4e, and 2e2μ). The final yield used in the analysis is a weighted average of the two independent estimates. To extract the shape of the m 4 distribution for the reducible background a maximum-likelihood fit is performed in each of the 22 event sub-categories defined in Sect. 6.1. For each sub-category, the expected "Z+X" yields from the OS and SS methods are binned as a function of m 4 and fitted to empirical functional forms built from Landau distributions [108]. In sub-categories with not enough statistics to perform a fit, the shape is extracted from the inclusive distribution in the corresponding final state.
The dominant systematic uncertainty on the reducible background estimation arises from the difference in the composition of the sample from which the misidentification rate is computed and that of the control regions. It is determined from the MC simulation and is found to be around 30%, depending on the final state. Additional sources of systematic uncertainty arise from the limited number of events in the control regions as well as in the region where the misidentification rates are computed.

Signal modeling
In order to generate an accurate signal model, the p T spectrum of the H boson, p H T , was tuned in the powheg simulation of the dominant gluon fusion production mode to better match the predictions from full phase space calculations implemented in the HRes generator [107,109,110].
In order to take advantage of the accuracy of the nnlops [111] simulation available for the ggH process, an event reweighting procedure is used. Events originating from the ggH process are subdivided into classes with 0, 1, 2, and ≥3 jets; the jets with p T > 30 GeV are clustered from all stable particles using the anti-k T algorithm with a distance parameter of 0.4, excluding the decay products of the H boson or associated vector bosons. The weights are obtained as the ratios of the p H T distributions from the nnlops and the powheg generators for each event class; the sum of the weights in each sample is normalized to the inclusive cross section.
The signal shape is parametrized by means of a doublesided Crystal Ball function [112] around m H ≈ 125 GeV. A Landau function is added in the total probability density function for the non-resonant part of the signal for the case of WH, ZH and tt H production modes. The signal shape is parametrized as a function of m H by performing a simultaneous fit of several mass points for all production modes in the 105 to 140 GeV mass range. Each parameter of the doublesided Crystal Ball function has a linear dependence on m H , for a total of 12 free parameters. An examples of the fit is shown in Fig. 2.

Systematic uncertainties
The systematic uncertainties are divided into experimental and theoretical. The main experimental uncertainties originate from the imperfect knowledge of the detector; the dominant sources are the uncertainties in the luminosity measurement, the lepton reconstruction and selection efficiency, the lepton and jet energy scale and resolution, the b tagging efficiency, and the reducible background estimate. The theoretical uncertainties account for the uncertainties in the modeling of the signal and background processes.
Both types of uncertainties can affect the signal selection, cause migrations between the event categories, and affect the signal or background shapes used in the fit. All the uncertainties affecting this analysis are modeled as nuisance parameters (NPs) that are profiled in the maximum likelihood fit described in Sect. 10.
In the combination of the three data-taking periods, all theoretical uncertainties are treated as correlated across these periods. The experimental uncertainties related to reconstruction and selection efficiency, the lepton energy scale Fig. 2 The shape of the parametric signal model for each year of simulated data, and for the sum of all years together. The black points represent weighted simulation events of the ggH production mechanism for m H = 125 GeV and the blue line the corresponding model. Also shown is the σ CB value (half the width of the narrowest interval containing 68% of the invariant mass distribution) in the gray shaded area. The contribution of the signal model from each year of data-taking is illustrated with the dotted lines. The models are shown for the 4e (upper) and 4μ (lower) final states in the untagged event category and resolution, and the b-tagging efficiency are also considered correlated across data-taking periods. Luminosity uncertainty is treated as partially correlated. All other experimental uncertainties are treated as uncorrelated. Correlated sources of uncertainty are assigned the same NP and uncorrelated sources have a dedicated NP in the likelihood fit described in Sect. 10.
The dominant sources of uncertainties and their effect on the analysis are discussed in detail in the following subsections. The impact of a NP on a parameter of interest (POI) is defined as the shift induced on POI when NP is varied by a ±1 standard deviation from its post-fit value, with all other parameters profiled as usual. The relative impact of the dominant systematic uncertainties on some of the measurements discussed in Sect. 10 is illustrated in Fig. 3. Here selection efficiency includes all the steps from trigger to impact parameter significance and finally identification and isolation requirements. The uncertainty ranges from 1   Lepton momentum scale and resolution uncertainties are estimated from dedicated studies on the Z → + − mass distribution in data and simulation. Events are classified according to the p T and η of one of the two leptons, determined randomly, and integrated over the other. The dilepton mass distributions are then fit by a Breit-Wigner parameterization convolved with the double-sided Crystal Ball function described in Sect. 8. The scale uncertainty is found to be 0.04% in the 4μ channel and 0.3% in the 4e channel, while the resolution uncertainty is 20% for both channels. In both cases full correlation between the leptons in the event is assumed. Both scale and resolution uncertainties alter the signal shape by allowing the corresponding parameters of the double-sided Crystal Ball function to vary. The impact is found to be non-negligible only in the case of fiducial cross section measurements.
The effects of the jet energy corrections are studied in a similar manner. While jet energy scale and smearing do not alter signal selection efficiency, they cause event migrations between the categories. They can also alter the shape of the discriminants, but the effect on the shape is negligible. The uncertainty in the jet energy scale ranges from 1% in the high jet p T range and increases up to 5% in the low jet p T range. The uncertainty in jet energy resolution ranges from 1 to 2%. A detailed description of the determination of the jet energy scale and smearing uncertainties can be found in [113]. The effect on the analysis is studied in detail by propagating the uncertainties and estimating the effect on event migration in each of the 22 sub-categories. Their impact on the inclusive measurements is found to be negligible. However, the impact is significant in measurements of the VBF and VH production modes and differential cross section measurements as a function of jet kinematics, where it is one of the leading sources of uncertainty. The uncertainty in the b-tagging efficiency is found to be 1% in the high jet p T range and increases up to 3% in the low jet p T range. The impact from the category migration is found to be negligible in all categories.
Finally, experimental uncertainties in the reducible background estimation, described in Sect. 7.2, originating from the background composition and misidentification rate uncertainties vary between 30 and 45% depending on the final state and category. However, the impact of this uncertainty on the measurements is found to be negligible.
Other sources of experimental uncertainties are also studied but their impact is negligible compared to the sources described above.

Theoretical uncertainties
Theoretical uncertainties that affect both the signal and background estimation include those related to the renormaliza-tion and factorization scales, and the choice of the PDF set. The uncertainty from the renormalization and factorization scales is determined by varying these scales between 0.5 and 2 times their nominal value, while keeping their ratio between 0.5 and 2. The uncertainty due to the PDF set is determined following the PDF4LHC recommendations by taking the root mean square of the variation of the results when using different replicas of the default NNPDF set [114,115]. The uncertainties just described have an effect both on the signal and background yields, as well as on the migration of events between the categories. An additional 10% uncertainty in the K factor used for the gg → ZZ prediction is applied as described in Sect. 7.1. A systematic uncertainty of 2% [32] in the branching fraction of H → 4 only affects the signal yield.
Theoretical uncertainties that affect the predictions of the STXS production bins are described in Ref. [32]. From here  ), and the ZZ backgrounds and rare electroweak backgrounds are normalized to the SM expectation, the Z+X background to the estimation from data. The vertical dashed lines denote the working points used in the event categorization. The SM H boson signal is separated into two components: the production mode which is targeted by the specific discriminant, and other production modes, where the gluon fusion process dominates. on we will refer to these uncertainties as the theoretical uncertainty scheme (THU).
The THU for the ggH process includes 10 NPs, which account for uncertainties in the cross section prediction for exclusive jet bins (including the migration between the 0 and 1-jet, as well as between the 1 and ≥2-jet bins), the 2 jet and ≥3 jet VBF phase space, migrations around the p H T bin boundaries at 10, 60, and 120 GeV, and the uncertainty in the p H T distribution due to missing higher order finite top quark mass corrections.
In the THU uncertainties for VBF and VH production, additional sources are introduced to account for the uncertainty in the modeling of the p H T , m jj and p Hjj T distributions, as well as that of the number of jets in the event. A total of 6 NPs account for the migrations of events across the m jj boundaries at 60, 120, 350, 700, 1000, and 1500 GeV. Two additional NPs account for migrations across the p H T = 200 GeV and p Hjj T = 25 GeV bin boundaries. Finally, a single source is introduced to account for migrations between the zero and one jet, as well as the the two or more jet bins. In each case, the uncertainty is computed by varying the renormalization and factorization scales and recalculating the fractional breakdown of the qqH STXS stage 1.2 cross sections.
A set of THU uncertainties is considered as NPs in the likelihood fit when signal strength modifiers, rather than STXS, are measured. In the STXS framework, THU uncertainties only enter at the interpretation step and are thus applied only to the SM cross section predictions.
Additional theoretical effects that only cause migration of signal and background events between categories originate from the modeling of the hadronization and the underlying event. The underlying event modeling uncertainty is determined by varying initial-and final-state radiation scales between 0.25 and 4 times their nominal value. The effects of the modeling of hadronization are determined by simulating additional events with the variation of the nominal pythia tune described in Sect. 3.

Results
The reconstructed four-lepton invariant mass distribution is shown in Fig. 4 for the 4e, 4μ and 2e2μ events together, and is compared with the expectations for signal and background processes. The error bars on the data points correspond to the intervals at 68% confidence level (CL) [116]. The observed distribution agrees with the expectation within the statistical uncertainties over the whole spectrum.
The reconstructed four-lepton invariant mass distribution is shown in Fig. 5 for the three 4 final states and is compared with the expectations from signal and background processes.  . Points with error bars represent the data and stacked histograms represent expected distributions of the signal and background processes. The yields of the different H boson production mechanisms with m H = 125 GeV, denoted as H(125), and those of the ZZ and rare electroweak backgrounds are normalized to the SM expectations, while the Z+X background yield is normalized to the estimate from the data. In the center and lower figures the SM H boson signal is separated into two components: the production mode which is targeted by the specific discriminant, and other production modes, where the gluon fusion process dominates The number of candidates observed in the data and the expected yields for 137 fb −1 , for the backgrounds and H boson signal after the full event selection, are given in Table 3 for each of the 22 reconstructed event categories (described in Sect. 6.2) for the 105 < m 4 < 140 GeV mass window around the Higgs boson peak. Fig. 6 shows the number of expected and observed events for each of the categories.
The reconstructed invariant masses of the Z 1 and Z 2 dilepton systems are shown in Fig. 7 for 118 < m 4 < 130 GeV, together with their 2D distribution in the 105 < m 4 < 140 GeV mass region. The distribution of the discriminants used for event categorization along with the corresponding working point values are shown in Fig. 8.
The results presented in Sects. 10.1 and 10.2 are extracted with a two-dimensional likelihood fit that relies on two variables, the four-lepton invariant mass m 4 and the matrix ele-  . 12 Result of the 2D likelihood scan for the μ f ≡ μ ggH, t t H,bb H,tH and μ V ≡ μ VBF,VH signal strength modifiers. The solid and dashed contours show the 68 and 95% CL regions, respectively. The cross indicates the best fit value, and the diamond represents the expected value for the SM Higgs boson ment kinematic discriminant D. The fiducial cross section measurements are extracted with a one-dimensional likelihood fit that relies only on the four-lepton invariant mass. The fit procedure and results are presented in Sect. 10.3. The fit is performed in the 105 < m 4 < 140 GeV mass region. The parameters of interest (POIs) are estimated with their corresponding confidence intervals using a profile likelihood ratio test statistic [117,118], in which the experimental and theoretical uncertainties are incorporated via NPs. The choice of the POIs depends on the specific measurement under consideration, while the remaining parameters are treated as NPs.
All the POIs considered in the analysis are forced to be greater than or equal to zero; this reflects the fact that the signal yield is substantially larger than the background yield in the mass range studied. Negative POIs would imply negative signal strength modifiers and a negative probability density function (pdf). We define a two-dimensional pdf as the product of two one-dimensional pdfs: The first term, P(m 4 ), is the unbinned analytical shape described in Sect. 8 for signals and Sect. 7 for backgrounds. The second term, P(D|m 4 ), is a binned template of D that is conditional to m 4 . This is achieved by creating a twodimensional template of m 4 vs. D and normalizing it to 1 for each bin of m 4 . In almost all sub-categories we use a decay-only kinematic discriminant (D = D kin bkg ) to separate the H boson signal from the background as defined in Eq. (3). Conversely, in the subcategories of the VBF-2jet-tagged, the D = D VBF+dec bkg discriminant (defined in Eq. (4)) is used, which is sensitive to the VBF production mechanism. Similarly, in two sub-categories of the VH-hadronic-tagged category, the D = D VH+dec bkg discriminant (defined in Eq. (5)) is used.
The ggH, VBF, WH, ZH and tt H samples are used to build different signal templates for each of the nineteen STXS production bins described in Sect. 6.1. Irreducible background templates are built starting from qq → ZZ and gg → ZZ samples. Finally, reducible background templates are built using data driven methods described in Sect. 7.2. Following the described procedure, P(D|m 4 ) templates are obtained for the twenty-two event categories and the three final states (4μ, 4e, 2e2μ).
The unbinned likelihood function, L( μ), is defined as the product over N observed events: where μ i is the signal strength modifier for the production bin with the four-lepton invariant mass is shown in Fig. 9 for the mass interval 105 < m 4 < 140 GeV. Their distributions for the mass interval 118 < m 4 < 130 GeV are shown in Fig. 10.

Signal strength modifier
A simultaneous fit to all categories is performed to extract the signal strength modifier, defined as the ratio of the observed H boson yield in the H → 4 decay channel to the standard model expectation.
The combined measurement of the inclusive signal strength modifier is measured to be μ = 0.94 +0.12 −0.11 or μ = 0.94 ± 0.07 (stat) +0.07 −0.06 (theo) +0.06 −0.05 (exp) at a fixed mass value m H = 125.38 GeV, which is the current most precise measurement of the H boson mass published by the CMS Collaboration [119]. In all subsequent fits, m H is fixed to this value. The dominant experimental sources of systematic uncertainty are the uncertainties in the lepton identification efficiencies and luminosity measurement, while the dominant theoretical source is the uncertainty in the total gluon fusion cross sec-tion. The contributions to the total uncertainty from experimental and theoretical sources are found to be similar in magnitude. The signal strength modifiers are further studied in terms of the five main SM Higgs boson production mechanisms, namely ggH, VBF, ZH, WH, and tt H. The contributions of the bb H and tH production modes are also taken into account. The relative normalizations of the bbH and the gluon fusion contributions are kept fixed in the fit, and so are the tH and tt H ones. The results are shown in Fig. 11 for the observed and expected profile likelihood scans of the inclusive signal strength modifier and those for the signal strength modifiers of the five main SM Higgs boson production mechanisms. The corresponding numerical values, including the decomposition of the uncertainties into statistical and systematic components, as well as the expected uncertainties, are given in Table 4.
The dependence of the measured signal strengths on the profiling of m H is checked and found to have a small impact both on the inclusive results and those in terms of the five main H boson production mechanisms, well within the measurement uncertainties. The best fit signal value changes at most by 4% and the profiled value of the mass is found to be m H = 125.09 +0.15 −0.14 (stat) GeV. It is important to note here that the precise determination of m H and the systematic uncer-  Table 5 Best fit values and ±1 standard deviation uncertainties for the measured cross sections (σ B) obs , the SM predictions (σ B) SM , and their ratio for the stage 0 STXS production bins at m H = 125.38 GeV for H → ZZ decay tainties that enter its measurement are beyond the scope of this analysis. Two signal strength modifiers, μ f ≡ μ ggH, tt H,bb H,tH and μ V ≡ μ VBF,VH , are introduced for the fermion and vectorboson induced contributions to the expected SM cross sec-tion. A two-parameter fit is performed simultaneously to the events reconstructed in all categories, leading to μ f = 0.96 +0.14 −0.12 and μ V = 0.82 +0.36 −0.31 . The expected values for m H = 125. 38 GeV are μ f = 1.00 +0.15 −0.13 and μ V = 1.00 +0.39 −0.33 . The 68 and 95% CL contours in the (μ f , μ V ) plane are shown in Fig. 12 and the SM predictions lie within the 68% CL regions of this measurement.

Simplified template cross section
The results for the H boson product of cross section times branching fraction for H → ZZ decay, (σ B) obs , and comparisons with the SM expectation, (σ B) SM , for the stages of production bins defined in Sect. 6.1, are shown in Fig. 13 for the stage 0 and in Fig. 14 for the merged stage 1.2. The corresponding numerical values are given in Tables 5 and 6 .
As discussed, the set of THU uncertainties described in Sect. 9.2 is not considered for the STXS measurements: THU uncertainties are model dependent and should be only considered in the interpretation of the results. Therefore, the THU uncertainties are included in the SM predictions of the cross section. The correlation matrices are shown in Fig. 15. The dominant experimental sources of systematic uncertainty are the same as for the signal strength modifiers measurement, while the dominant theoretical source is the uncertainty in the category migration for the ggH process.

Fiducial cross section
In this section the cross section measurement for the process pp → H → 4 within a fiducial volume that closely matches the reconstruction level selection is presented. In particular, the integrated fiducial cross section is measured as well as differential cross sections as a function of the transverse momentum of the H boson ( p H T ), its rapidity (|y H |), the number of associated jets (N j ), and the transverse momentum of the leading jet ( p j T ). These measurements are largely independent of the assumptions on the relative fractions and kinematic distributions of the individual production modes. The definition of the fiducial volume is based on generator-level quantities and is identical to that in Ref.
[18]. In order to reduce the experimental uncertainties, only jets with p j T > 30 GeV and |η j | < 2.5 are considered for the differential cross sections as a function of jet observables. An increase in model dependence compared to Ref.
[25] is observed when using the ZZ candidate selection at reconstruction level where the candidate with the best D kin bkg discriminant value is chosen. Therefore, the fiducial cross section measurement is performed using the event selection algorithm in Ref. [25]. Specifically, the Z 1 candidate is chosen to be the one with m(Z 1 ) closest to the nominal Z boson mass, and in cases where multiple Z 2 candidates satisfy all criteria, the pair of leptons with the largest sum of the transverse momenta magnitudes is chosen. The full fiducial volume definition is detailed in Table 7 and the acceptance for various SM production modes is given in Table 8.
A maximum likelihood fit of the signal and background parameterizations to the observed 4 mass distribution, N obs (m 4 ), is performed to extract the integrated fiducial cross section for the process pp → H → 4 (σ fid ). The fit is carried out inclusively (i.e., without any event categorization) and does not use the D kin bkg observable in order to minimize the model dependence. The fit is performed simultaneously in all final states and assumes a H boson mass m H = 125.38 GeV, while the branching fractions of the H boson to different final states (4e, 4μ, 2e2μ) are free parameters in the fit. The systematic uncertainties described in Sect. 9 are included in the form of NPs and the results are obtained using an asymptotic approach [118] with a test statistic based on the profile likelihood ratio [117]. This procedure accounts for the unfolding of detector effects from the observed distributions and is the same as in Refs. [25,120].   The number of expected events in each final state f and in each bin i of a given observable is expressed as a function of m 4 as: The shape of the resonant signal contribution, P res (m 4 ), is described by a double-sided Crystal Ball function as discussed in Sect. 8, and the normalization is used to extract the fiducial cross section. The non-resonant signal function, P nonres (m 4 ), is determined by the WH, ZH, and tt H contributions where one of the leptons from the H boson decay is lost or not selected. It is modeled by a Landau distribution with shape parameters constrained in the fit to be within a range determined from simulation. This contribution is referred to as the "combinatorial signal" and is treated as a background in this measurement.
The quantity f i, j represents the detector response matrix that maps the number of expected events in bin j of a given observable at the fiducial level to the number of expected events in bin i at the reconstruction level. This response matrix is determined using simulated signal samples and includes corrections for residual differences between data and simulation. In the case of the integrated fiducial cross section measurement, the response matrices become single numbers, which are listed in Table 8 for different SM production mechanism.
An additional resonant contribution arises from events which are accepted but do not originate from the fiducial phase space. These events are due to detector effects that cause differences between the quantities used for the fiducial phase space definition and the corresponding quantities at the reconstruction level. This contribution is treated as background and is referred to as the "non-fiducial signal" contribution. A simulated sample is used to verify that the shape of the distribution for these events is identical to that of the fiducial signal, and its normalization is fixed to the corresponding fraction of the fiducial signal. The value of   Fig. 16. The corresponding numerical values, including the decomposition of the uncertainties into statistical and systematic components, and the corresponding expected uncertainties, are given in Table 9.
The measured differential cross sections as a function of the H boson transverse momentum and rapidity are shown in Fig. 17. The corresponding numerical values are given in Tables 10 and 11 . Finally, the measured differential cross sections as a function of the number of associated jets and the transverse momentum of the leading jet are shown in Fig. 18. The corresponding numerical values are given in Tables 12  and 13.
For all the fiducial measurements the dominant systematic uncertainties are those on the lepton identification efficiencies and luminosity measurement, while the theoretical uncertainties are smaller. In order to assess the model dependence of the measurement, the unfolding procedure is repeated using different response matrices created by varying the relative fraction of each SM production mode within its experimental constraints. The uncertainty is negligible with respect to the experimental systematic uncertainties.

Summary
Several measurements of the Higgs boson production in the four-lepton final state at √ s = 13 TeV have been presented, using data samples corresponding to an integrated luminosity of 137 fb −1 . Thanks to a large signal-to-background ratio and the complete reconstruction of the final state decay products, this channel enables a detailed study of the Higgs boson production properties. The measured signal strength modifier is μ = 0.94 ± 0.07 (stat) +0.07 −0.06 (theo) +0.06 −0.05 (exp) and the integrated fiducial cross section is measured to be σ fid = 2.84 +0.23 −0.22 (stat) +0.26 −0.21 (syst) fb with a standard model prediction of 2.84±0.15 fb for the same fiducial region.. The signal strength modifiers for the main Higgs boson production modes are also reported. A new set of measurements, designed to quantify the different Higgs boson production processes in specific kinematical regions of phase space, have also been presented. The differential cross sections as a function of the transverse momentum and rapidity of the Higgs boson, the number of associated jets, and the transverse momentum of the leading associated jet are determined. All results are consistent, within their uncertainties, with the expectations for the standard model Higgs boson.