Higgs boson production cross-section measurements and their EFT interpretation in the $4\ell$ decay channel at $\sqrt{s}$ = 13 TeV with the ATLAS detector

Higgs boson properties are studied in the four-lepton decay channel (where lepton = $e$, $\mu$) using 139 fb$^{-1}$ of proton-proton collision data recorded at $\sqrt{s}$ = 13 TeV by the ATLAS experiment at the Large Hadron Collider. The inclusive cross-section times branching ratio for $H\to ZZ^*$ decay is measured to be $1.34 \pm 0.12$ pb for a Higgs boson with absolute rapidity below 2.5, in good agreement with the Standard Model prediction of $1.33 \pm 0.08$ pb. Cross-sections times branching ratio are measured for the main Higgs boson production modes in several exclusive phase-space regions. The measurements are interpreted in terms of coupling modifiers and of the tensor structure of Higgs boson interactions using an effective field theory approach. Exclusion limits are set on the CP-even and CP-odd `beyond the Standard Model' couplings of the Higgs boson to vector bosons, gluons and top quarks.


Introduction
The observation of the Higgs boson by the ATLAS and CMS experiments [1,2] with the Large Hadron Collider (LHC) Run 1 data set at centre-of-mass energies of √ s = 7 TeV and 8 TeV was a major step towards an understanding of the electroweak (EW) symmetry breaking mechanism [3][4][5]. Tests of its spin and CP quantum numbers strongly indicate that the observed particle is of scalar nature and that the dominant coupling structure is CP-even, consistent with the Standard Model (SM) expectation [6][7][8]. The measurements of the Higgs boson production and differential cross-sections, branching ratios, and the derived constraints on couplingstrength modifiers, assuming the SM coupling structure, have also shown no significant deviation from the predictions for the SM Higgs boson with a mass of 125 GeV [9-12]. Furthermore, constraints have been set on various coupling parameters beyond the SM (BSM) that modify the tensor structure of the Higgs boson couplings to SM particles [8,[13][14][15][16][17][18][19][20].
Motivated by a clear Higgs boson signature and a high signal-to-background ratio in the H → Z Z * → 4 decay channel (where = e or μ), the updated measurements of the Higgs boson coupling properties in this channel are presented using the entire Run 2 data set with 139 fb −1 of protonproton ( pp) collision data collected at √ s = 13 TeV by the ATLAS detector between 2015 and 2018. Three types of results are presented in this paper: (i) measurements of the Higgs boson production cross-sections times branching ratio, hereafter referred to as cross-sections, for the main production modes in several exclusive phase-space bins in dedicated fiducial regions; (ii) interpretation of the measurements in terms of constraints on the Higgs boson coupling-strength modifiers within the κ-framework [21]; and (iii) interpretation of the measurements in terms of modifications to the tensor structure of Higgs boson couplings using an effective field theory (EFT) approach.
In addition to a nearly four times higher integrated luminosity, there are several other important differences compared to the previous results in this analysis channel [17]: • an improved lepton isolation to mitigate the impact of additional pp interactions in the same or neighbouring bunch crossings (pile-up), • an improved jet reconstruction using a particle flow algorithm [22], • additional event categories for the classification of Higgs boson candidates, • new discriminants to enhance the sensitivity to distinguish the various production modes of the SM Higgs boson, • the use of data sidebands to constrain the dominant Z Z * background process, • a dedicated control region to constrain the background in the reconstructed event categories probing tt H production, • improved estimates of Z +jets, tt, and W Z backgrounds, and • an EFT interpretation, based on a parameterisation of the cross-sections rather than a direct parameterisation of the reconstructed event yields.

Simplified template cross-sections
In the framework of Simplified Template Cross Sections (STXS) [23][24][25], exclusive regions of phase space are defined for each Higgs boson production mechanism. These phasespace regions, referred to as production bins, are defined to reduce the dependence on theoretical uncertainties that directly fold into the measurements and at the same time maximise the experimental sensitivity to measure the bins, enhance the contribution from possible BSM effects, and allow measurements from different Higgs boson decay modes to be combined. The number of production bins is limited to avoid loss of measurement sensitivity for a given amount of integrated luminosity. The definitions of the production bins used for this measurement are shown in the left panel of Fig. 1 (shaded area). All production bins are defined for Higgs bosons with rapidity |y H | < 2.5 and no requirement is placed on the particlelevel leptons. Two sets of production bins with different granularity are considered, as a trade-off between statistical and theoretical uncertainties.
The first set of production bins (Production Mode Stage) [24] is defined according to the Higgs boson production modes: gluon-gluon fusion (ggF), vector-boson fusion (VBF) and associated production with vector bosons (VH, where V = W or Z ) or top quark pairs (tt H ). Since b-jets from bbH associated production are emitted at small angles relative to the beam axis and usually outside of the detector acceptance, the bbH and ggF Higgs boson production modes have similar signatures and acceptances. Their contributions are considered together with their relative ratio fixed to the SM prediction. In the following, the sum of their contributions is referred to as ggF. Similarly, single top production (tH) is considered together with tt H , with their relative ratio fixed to the SM prediction. In contrast to the Stage-0 production bins described in Ref. [24], the VH events with hadronic decays of the vector boson V are included in the VH production bin rather than in the ggF or VBF bins. In this way, each of the four main Higgs boson production modes can be measured separately.
The second set of production bins (Reduced Stage 1.1) is more exclusive than the first one. Starting from the production bins of a more granular Stage 1.1 set [25], several production bins are merged as the full set of bins cannot be measured separately in the H → Z Z * → 4 channel with the current data sample. The definitions of the bins are based on the multiplicity of particle-level jets, the Higgs boson transverse momentum p H T and the invariant mass m j j of the two jets with the highest transverse momentum. Particle-level jets are built from all stable particles (particles with lifetime cτ >10 mm) including neutrinos, photons, and leptons from hadron decays or those produced in the parton shower. The anti-k t jet reconstruction algorithm [26,27] with a radius parameter R = 0.4 is used. All Higgs boson decay products, as well as the leptons and neutrinos from the decays of the associated V bosons are excluded from the jet building, while the decay products from hadronically decaying associated V bosons, are included. The jets are required to have p T > 30 GeV, with no restrictions on rapidity.
Events from ggF production and gg → Z H production with a hadronically decaying Z boson are split into seven common production bins. Six bins have a Higgs boson transverse momentum below 200 GeV, while the seventh bin with  Events with no jets are split into two bins with p H T below and above 10 GeV. Events with one jet are split into three bins with p H T below 60 GeV, between 60 and 120 GeV, and above 120 GeV. Finally, Higgs boson events with two or more jets are combined into one bin. The bins are respectively denoted by gg2H-0 j-p H T -Low, gg2H-0 j-p H T -High, gg2H-1 j-p H T -Low, gg2H-1 j-p H T -Med, gg2H-1 j-p H T -High and gg2H-2 j. As described in Ref. [25], VBF and VH production with hadronically decaying associated V bosons represent the t-channel and s-channel contributions to the same electroweak qq H production process and are therefore considered together for further splitting. Three bins are defined: one bin, sensitive to BSM contributions (qq2Hqq-BSM), with p H T above 200 GeV and m j j above 350 GeV; one bin (qq2Hqq-VH) with m j j between 60 and 120 GeV to target the VH production mode; and one bin (qq2Hqq-VBF) with the Higgs boson not satisfying these criteria to ensure sensitivity to the VBF process. qq H events in which one or both jets have transverse momenta below the 30 GeV threshold are treated as a part of the qq2Hqq-VBF bin.
The VH process with the associated V boson decaying leptonically is considered separately (VH-Lep). The leptonic decay includes the decays into τ -leptons and neutrino pairs. The tt H production bin remains the same as in the Production Mode Stage.
The middle-right and right panels of Fig. 1 summarise the corresponding categories of reconstructed events in which the cross-section measurements and background estimations are performed. These are described in detail in Sect. 5.

Higgs boson couplings in the κ-framework
To probe physics beyond the SM, the measured production cross-sections are interpreted within a leading-ordermotivated κ-framework [21], in which a set of coupling modifiers κ is introduced to parameterise deviations from the SM predictions of the Higgs boson couplings to SM bosons and fermions. The framework assumes that the data origi-nate from a single CP-even Higgs boson state with a mass of 125 GeV and the tensor coupling structure of the SM for its interactions. Only the coupling strengths are allowed to be modified by the BSM processes. The Higgs boson width is assumed to be small enough such that the narrow-width approximation is valid, allowing the Higgs boson production and decay to be factorised: where σ i is the production cross-section via the initial state i, B and f are the branching ratio and partial decay width for the decay into the final state f , respectively, and H is the total width of the Higgs boson. For a Higgs boson production and decay process via couplings i and f , respectively, coupling-strength modifiers are defined as

Tensor structure of Higgs boson couplings in the effective field theory approach
The κ-framework assumes that the tensor structure of the Higgs boson couplings is the same as in the SM. In order to probe for possible non-SM contributions to the tensor structure of the Higgs boson couplings, the measured simplified template cross-sections are interpreted using an EFT approach. In this approach, which exploits exclusive kinematical regions of the Higgs boson production and decay phase space, the BSM interactions are introduced via additional higher-dimensional operators O (d) i of dimension d, supplementing the SM Lagrangian L SM ,

The parameters C (d)
i specify the strength of new interactions and are known as the Wilson coefficients, and is the scale of new physics. Only dimension-six operators are considered for this paper, since the dimension-five and dimension-seven operators violate lepton and baryon number conservation and the impact of higher-dimensional operators is expected to be suppressed by more powers of the cutoff scale [28]. For energies less than the scale of new physics, only the ratio c i = C (d=6) i / 2 can be constrained by the data. Constraints are set on the Wilson coefficients defined within the Standard Model Effective Field Theory (SMEFT) formalism [29] in the Warsaw basis [30]. The measurements in the H → Z Z * → 4 channel do not provide sensitivity for simultaneous constraints on the full set of these coefficients. To reduce the number of relevant parameters, a minimal flavour-violating scenario is assumed and only operators affecting the Higgs boson cross-section at tree level are considered. Operators affecting only double Higgs boson production and those affecting the Higgs boson couplings to down-type quarks and leptons are neglected due to limited sensitivity. The impact of these operators on the total Higgs boson decay width is also neglected.
The remaining ten operators (see Table 1) comprise five CP-even and five CP-odd ones. The CP-even operators describing interactions between the Higgs boson and gluons and the top-Yukawa interactions are associated with the Wilson coefficients c H G and c u H from Ref. [29], respectively. Similarly, the CP-even Higgs  The constraints on the Wilson coefficients can be derived by comparing the expected with the measured simplified template cross-sections. For that purpose, the corresponding expected signal production cross-sections, the branching ratio and the signal acceptances are parameterised in terms of the Wilson coefficients. The dependence of signal production cross-sections on the EFT parameters can be obtained from its separation into three components: where the first term on the right-hand side is the squared matrix element for the SM, the second term represents the interference between the SM and dimension-six EFT amplitudes and the third term comprises the pure BSM contribution from dimension-six EFT operators alone. Following this expression, the dependence of the Higgs boson cross-section σ p ( c) in a given production bin p on a set of Wilson coefficients c is parameterised relative to the SM prediction σ where the coefficients A p i and B p i j are independent of c and are determined from simulation. A similar procedure is applied to obtain from simulation the EFT parameterisation of the branching ratio B 4 for the H → Z Z * → 4 decay from the partial ( 4 ) and total decay width ( tot ) parame- Table 1 Summary of EFT operators in the SMEFT formalism that are probed in the H → Z Z * → 4 channel. The corresponding tensor structure in terms of the SM fields from Ref. [29] is shown together with the associated Wilson coefficients, the affected production vertices and the impact on the H → Z Z * decay vertex. The Higgs doublet field and its complex conjugate are denoted as H and H , respectively. The left-handed quark doublets of flavour p (the right-handed up-type quarks) are denoted q p (u r ). V μν ( V μν = μνρσ V ρσ ) is the (dual) field strength tensor for a given gauge field V = G, W, B. The bosonic operators with (without) a dual field strength tensor are CP-odd (CP-even where the total decay width is the sum of all partial decay widths f related to the decay mode f . The procedure for the parameterisation of the cross-sections and the branching ratios is described in more detail in Ref. [31]. The criteria employed in the selection of four-lepton candidates introduce an additional dependence of the signal acceptance on the EFT parameters. This is taken into account in the interpretation, as discussed in Sect. 10.

ATLAS detector
The ATLAS detector [32][33][34] at the LHC is a multipurpose particle detector with a forward-backward symmetric cylindrical geometry 1 and a nearly 4π coverage in solid angle. It consists of an inner tracking detector (ID) surrounded by a thin superconducting solenoid, which provides a 2 T axial magnetic field, electromagnetic (EM) and hadron calorimeters, and a muon spectrometer (MS). The inner tracking detector covers the pseudorapidity range |η| < 2.5. It consists of silicon pixel, silicon microstrip, and transition radiation tracking detectors. A lead/liquid-argon (LAr) sampling 1 ATLAS uses a right-handed coordinate system with its origin at the nominal interaction point (IP) in the centre of the detector and the zaxis along the beam pipe. The x-axis points from the IP to the centre of the LHC ring, and the y-axis points upwards. Cylindrical coordinates (r, φ) are used in the transverse plane, φ being the azimuthal angle around the z-axis. The pseudorapidity is defined in terms of the polar angle θ as η = − ln tan(θ/2). Angular distance is measured in units of calorimeter provides electromagnetic energy measurements in the pseudorapidity range |η| < 3.2 with high granularity. A steel/scintillator-tile hadron calorimeter covers the central pseudorapidity range (|η| < 1.7). The endcap and forward regions are instrumented up to |η| = 4.9 with LAr calorimeters for both the EM and hadronic energy measurements. The calorimeters are surrounded by the MS and three large aircore toroidal superconducting magnets with eight coils each. The field integral of the toroid magnets ranges between 2.0 and 6.0 Tm across most of the detector. The MS includes a system of precision tracking chambers and fast detectors for triggering, covering the region |η| < 2.7. Events are selected using a first-level trigger implemented in custom electronics, which reduces the event rate to a maximum of 100 kHz using a subset of detector information. Software algorithms with access to the full detector information are then used in the high-level trigger to yield a recorded event rate of about 1 kHz [35].

Data set and event simulation
The full ATLAS Run 2 data set, consisting of pp collision data at √ s = 13 TeV taken between 2015 and 2018, is used for this analysis. The total integrated luminosity after imposing data quality requirements [36] is 139 fb −1 .
The production of the SM Higgs boson via gluon-gluon fusion, via vector-boson fusion, with an associated vector boson and with a top quark pair was modelled with the Powheg-Box v2 Monte Carlo (MC) event generator [37][38][39]. For ggF, the PDF4LHC next-to-next-to-leading-order (NNLO) set of parton distribution functions (PDF) was used, while for all other production modes, the PDF4LHC next-toleading-order (NLO) set was used [40].
The simulation of ggF Higgs boson production used the Powheg method for merging the NLO Higgs boson + jet cross-section with the parton shower and the multi-scale improved NLO (MINLO) method [41][42][43][44] to simultaneously achieve NLO accuracy for the inclusive Higgs boson production. In a second step, a reweighting procedure (NNLOPS) [45,46], exploiting the Higgs boson rapidity distribution, was applied using the HNNLO program [47,48] to achieve NNLO accuracy in the strong coupling constant α S . The transverse momentum spectrum of the Higgs boson obtained with this sample is compatible with the fixed-order calculation from HNNLO and the resummed calculation at next-tonext-to-leading-logarithm accuracy matched to NNLO fixedorder with Hres2.3 [49,50].
The matrix elements of the VBF, qq → V H, and tt H production mechanisms were calculated up to NLO in QCD. For VH production, the MINLO method was used to merge 0-jet and 1-jet events [41,43,[51][52][53][54]. The gg → Z H contribution was modelled at leading order (LO) in QCD.
The production of a Higgs boson in association with a bottom quark pair (bbH ) was simulated at NLO with Mad-Graph5_aMC@NLO v2.3.3 [55,56], using the CT10 NLO PDF [57]. The production in association with a single top quark (t H+X where X is either jb or W , defined in the following as t H) [58,59] was simulated at NLO with Mad-Graph5_aMC@NLO v2.6.0 using the NNPDF3.0nlo PDF set [60].
For all production mechanisms, the Pythia 8 [61] generator was used for the H → Z Z * → 4 decay with = (e, μ) as well as for parton showering, hadronisation and the underlying event. The contribution of the Z → τ τ decays is shown to have a negligible impact on the final result. The event generator was interfaced to EvtGen v1.2.0 [62] for simulation of the bottom and charm hadron decays. For the ggF, VBF and VH processes, the AZNLO [63] set of tuned parameters was used, while the A14 [64] set was used for tt H , bbH and t H processes. All signal samples were simulated for a Higgs boson mass m H = 125 GeV.
For additional cross-checks, the ggF sample was also generated with MadGraph5_aMC@NLO. This simulation is accurate at NLO QCD accuracy for zero, one and two additional partons merged with the FxFx merging scheme [55,65]. The events were showered using the Pythia 8 generator with the A14 set of tuned parameters.

t H W and t H jb production modes using
MadGraph5_aMC@NLO and the NNPDF23lo PDF. The BSM signal is defined by the flavour symmetric SMEFT-sim_A_U35_MwScheme_UFO_v2.1 model [29,106], which incorporates the SMEFT dimension-six operators in the standard Universal FeynRules Output format created using the FeynRules framework [107,108]. The light quarks (u, d, s and c) and leptons are assumed to be massless in the model. The generated events were showered with Pythia 8, using the CKKW-L matching scheme to match matrix element and parton shower computations with different jet multiplicities [61]. The A14 set of tuned parameters was used. All processes were simulated in the four-flavour scheme, apart from the t H W production, for which the five-flavour scheme was used [55]. Page 7 of 54 957 The Z Z * continuum background from quark-antiquark annihilation was modelled using Sherpa v2.2.2 [109][110][111][112], which provides a matrix element calculation accurate to NLO in α S for 0-jet and 1-jet final states and LO accuracy for 2-jets and 3-jets final states. The merging with the Sherpa parton shower [113] was performed using the ME+PS@NLO prescription [114]. The NLO EW corrections were applied as a function of the invariant mass m Z Z * of the Z Z * system [115,116].
The gluon-induced Z Z * production was modelled by Sherpa v2.2.2 [109][110][111] at LO in QCD for 0-jet and 1-jet final states. The higher-order QCD effects for the gg → Z Z * continuum production cross-section were calculated for massless quark loops [117][118][119] in the heavy topquark approximation [120], including the interference with gg → H * → Z Z processes [121,122]. The gg → Z Z simulation was scaled by a K -factor of 1.7 ± 1.0, which is defined as the ratio of the higher-order to the leading-order cross-section predictions.
Production of Z Z * via vector-boson scattering was simulated with the Sherpa v2.2.2 [112] generator. The LOaccurate matrix elements were matched to a parton shower using the MEPS@LO prescription.
For all Z Z * processes modelled using Sherpa, the NNPDF3.0nnlo PDF set [60] was used, along with a dedicated set of tuned parton-shower parameters.
For additional checks, the qq-initiated Z Z * continuum background was also modelled using Powheg-Box v2 and MadGraph5_aMC@NLO, using the CT10 [57] and the PDF4LHC NLO PDF set, respectively. For the former, the matrix element was generated at NLO accuracy in QCD and effects of singly resonant amplitudes and interference effects due to Z /γ * were included. For the latter, the simulations are accurate to NLO in QCD for zero and one additional parton merged with the FxFx merging scheme. For both, the Pythia 8 generator was used for the modelling of parton showering, hadronisation, and the underlying event. The AZNLO and A14 sets of tuned parameters were used for the simulations performed with Powheg-Box v2 and Mad-Graph5_aMC@NLO generators, respectively.
The WZ background [123] was modelled at NLO accuracy in QCD using Powheg-Box v2 with the CT10 PDF set and was interfaced to Pythia 8, using the AZNLO set of tuned parameters for modelling of parton showering, hadronisation, and the underlying event and to EvtGen v1.2.0 for the simulation of bottom and charm hadron decays. The triboson backgrounds ZZZ, WZZ, and WWZ with four or more prompt leptons (VVV ) were modelled at NLO accuracy for the inclusive process and at LO for up to two additional parton emissions using Sherpa v2.2.2.
The simulation of tt Z events with both top quarks decaying semileptonically and the Z boson decaying leptonically was performed with MadGraph5_aMC@NLO using the NNPDF3.0nlo [60] PDF set interfaced to Pythia 8 using the A14 set of tuned parameters, and the total crosssection was normalised to a prediction computed at NLO in the QCD and EW couplings [98]. For modelling comparisons, Sherpa v2.2.1 was used to simulate tt Z events at LO. The t W Z, tt W W , tt W Z, tt Zγ , tt Z Z, ttt, tttt and t Z background processes were simulated with Mad-Graph5_aMC@NLO interfaced to Pythia 8, using the A14 set of tuned parameters. These processes are collectively referred to as the tXX process.
The modelling of events containing Z bosons with associated jets (Z + jets) was performed using the Sherpa v2.2.1 generator. Matrix elements were calculated for up to two partons at NLO and four partons at LO using Comix [110] and OpenLoops [111], and merged with the Sherpa parton shower [113] using the ME+PS@NLO prescription [114]. The NNPDF3.0nnlo PDF set is used in conjunction with dedicated set of tuned parton-shower parameters.
The tt background was modelled using Powheg-Box v2 with the NNPDF3.0nlo PDF set. This simulation was interfaced to Pythia 8, using the A14 set of tuned parameters, for parton showering, hadronisation, and the underlying event, and to EvtGen v1.2.0 for heavy-flavour hadron decays. Simulated Z + jets and tt background samples were normalised to the data-driven estimates described in Sect. 6.
Generated events were processed through the ATLAS detector simulation [124] within the Geant4 framework [125] and reconstructed in the same way as collision data. Additional pp interactions in the same and nearby bunch crossings were included in the simulation. Pile-up events were generated using Pythia 8 with the A2 set of tuned parameters [126] and the MSTW2008LO PDF set [127]. The simulation samples were weighted to reproduce the distribution of the number of interactions per bunch crossing observed in data.

Event reconstruction
The selection and categorisation of the Higgs boson candidate events rely on the reconstruction and identification of electrons, muons, and jets, closely following the analyses reported in Refs. [17,128].
Proton-proton collision vertices are constructed from reconstructed trajectories of charged particles in the ID with transverse momentum p T > 500 MeV. Events are required to have at least one collision vertex with at least two associated tracks. The vertex with the highest p 2 T of reconstructed tracks is selected as the primary vertex of the hard interaction. The data are subjected to quality requirements to reject Eur. Phys. J. C (2020) 80:957 events in which detector components were not operating correctly. Electron candidates are reconstructed from energy clusters in the electromagnetic calorimeter that are matched to ID tracks [129]. A Gaussian-sum filter algorithm [130] is used to compensate for radiative energy losses in the ID for the track reconstruction, while a dynamical, topological cellbased approach for cluster building is used to improve the energy resolution relative to the previous measurements in Refs. [17,128], in particular for the case of bremsstrahlung photons. Electron identification is based on a likelihood discriminant combining the measured track properties, transition radiation response, electromagnetic shower shapes and the quality of the track-cluster matching. The 'loose' likelihood criteria, applied in combination with track hit requirements, provide an electron reconstruction and identification efficiency of at least 90% for isolated electrons with p T > 30 GeV and 85%-90% below [129]. Electrons are required to have E T > 7 GeV and pseudorapidity |η| < 2.47, with their energy calibrated as described in Ref. [129].
Muon candidate reconstruction [131] within the range |η| < 2.5 is primarily performed by a global fit to fully reconstructed tracks in the ID and the MS, with a 'loose' [131] identification criterion applied. This criterion has an efficiency of at least 98% for isolated muons with p T = 5 GeV and rises to 99.5% at higher p T . At the centre of the detector (|η| < 0.1), which has a reduced MS geometrical coverage, muons are also identified by matching a fully reconstructed ID track to either an MS track segment or a calorimeter energy deposit consistent with a minimum-ionising particle (calorimeter-tagged muons). For these two cases, the muon momentum is measured from the ID track alone. In the forward MS region (2.5 < |η| < 2.7), outside the full ID coverage, MS tracks with hits in the three MS layers are accepted and combined with forward ID tracklets, if they exist (stand-alone muons). Calorimeter-tagged muons are required to have p T > 15 GeV. For all other muon candidates, the transverse momentum is required to be greater than 5 GeV. The muon momentum is calibrated using the procedure described in Ref. [131]. Muons with transverse impact parameter greater than 1 mm are rejected. 2 Additionally, muons and electrons are required to have a longitudinal impact parameter (|z 0 sin θ |) less than 0.5 mm.
Jets are reconstructed using a particle flow algorithm [22] from noise-suppressed positive-energy topological clusters [132] in the calorimeter using the anti-k t algorithm [26,27] with a radius parameter R = 0.4. Energy deposited in the 2 The transverse impact parameter d 0 of a charged-particle track is defined in the transverse plane as the distance from the primary vertex to the track's point of closest approach. The longitudinal impact parameter z 0 is the distance in the z direction between this track point and the primary vertex. calorimeter by charged particles is subtracted and replaced by the momenta of tracks that are matched to those topological clusters. Compared to only using topological clusters, jets reconstructed with the particle flow algorithm with p T > 30 GeV have approximately 10% better transverse momentum resolution. The two different algorithms have similar resolution for p T above 100 GeV. The jet four-momentum is corrected for the calorimeter's non-compensating response, signal losses due to noise threshold effects, energy lost in non-instrumented regions, and contributions from pileup [22,133,134]. Jets are required to have p T > 30 GeV and |η| < 4.5. Jets from pile-up with |η| < 2.5 are suppressed using a jet-vertex-tagger multivariate discriminant [135,136]. Jets with |η| < 2.5 containing b-hadrons are identified using the MV2c10 b-tagging algorithm [137,138], and its 60%, 70%, 77% and 85% efficiency working points are combined into a pseudo-continuous b-tagging weight [139] that is assigned to each jet.
Ambiguities are resolved if electron, muon, or jet candidates overlap in geometry or share the same detector information. If the two calorimeter energy clusters from the two electron candidates overlap, the electron with the higher E T is retained. If a reconstructed electron and muon share the same ID track, the muon is rejected if it is calorimeter-tagged; otherwise the electron is rejected. Reconstructed jets geometrically overlapping in a cone of radial size R = 0.1 (0.2) with a muon (an electron) are also removed.
The missing transverse momentum vector, E miss T , is defined as the negative vector sum of the transverse momenta of all the identified and calibrated leptons, photons and jets and the remaining unclustered energy, where the latter is estimated from lowp T tracks associated with the primary vertex but not assigned to any lepton, photon, hadronically decaying τ -lepton or jet candidate [140,141]. The missing transverse momentum (E miss T ) is defined as the magnitude of E miss T .

Selection of the Higgs boson candidates
A summary of the event selection criteria is given in Table 3. Events were triggered by a combination of single-lepton, dilepton and trilepton triggers with different transverse momentum thresholds. Single-lepton triggers with the lowest thresholds had strict identification and isolation requirements. Both the high-threshold single-lepton triggers and the multilepton triggers had looser selection criteria. Due to an increasing peak luminosity, these thresholds increased slightly during the data-taking periods [142,143]. For singlemuon triggers, the p T threshold ranged from between 20 and 26 GeV, while for single-electron triggers, the p T threshold ranged from 24 to 26 GeV. The global trigger efficiency for signal events passing the final selection is about 98%. In the analysis, at least two same-flavour and oppositecharge lepton pairs (hereafter referred to as lepton pairs) are Table 3 Summary of the criteria applied to the selected Higgs boson candidate in each event. The mass threshold m min is defined in Sect. 4.1

Leptons and jets
Electrons E T > 7 GeV and |η| < 2.47 Muons p T > 5 GeV and |η| < 2.7, calorimeter-tagged: p T > 15 GeV Lepton isolation -The amount of isolation E T after summing the track-based and 40% of the calorimeter-based contribution must be smaller than 16% of the lepton p T Impact parameter -Electrons: To minimise the background contribution from non-prompt muons, at most one calorimeter-tagged or stand-alone muon is allowed per quadruplet. The lepton pair with the invariant mass m 12 (m 34 ) closest (second closest) to the Z boson mass [144] in each quadruplet is referred to as the leading (subleading) lepton pair. Based on the lepton flavour, each quadruplet is classified into one of the following decay final states: 4μ, 2e2μ, 2μ2e and 4e, with the first two leptons always representing the leading lepton pair. In each of these final states, the quadruplet with m 12 closest to the Z boson mass has priority to be considered for the selection of the final Higgs boson candidate. In case additional prompt leptons are present in the event, the priority may change due to the matrix-element based pairing as described later on. All quadruplets are therefore required to pass the following selection criteria.
To ensure that the leading lepton pair from the signal originates from a Z boson decay, the leading lepton pair is required to satisfy 50 GeV < m 12 < 106 GeV. The subleading lepton pair is required to have a mass m min < m 34 < 115 GeV, where m min is 12 GeV for the four-lepton invariant mass m 4 below 140 GeV, rising linearly to 50 GeV at m 4 = 190 GeV and then remaining at 50 GeV for all higher m 4 values. This criterion suppresses the contributions from processes in which an on-shell Z boson is produced in association with a leptonically decaying meson or virtual photon. In the 4e and 4μ final states, the two alternative opposite-charge lepton pairings within a quadruplet are required to have a dilepton mass above 5 GeV to suppress the J/ψ background. All leptons in the quadruplet are required to have an angular separation of R > 0.1.
Each electron (muon) track is required to have a transverse impact parameter significance |d 0 /σ (d 0 )| < 5 (3), to suppress the background from heavy-flavour hadrons. Reducible background from the Z +jets and tt processes is further suppressed by imposing track-based and calorimeter-based isolation criteria on each lepton [131,145]. A scalar p T sum (track isolation) is made from the tracks with p T > 500 MeV which either originate from the primary vertex or have |z 0 sin θ | < 3 mm if not associated with any vertex and lie within a cone of R = 0.3 around the muon or electron. Above a lepton p T of 33 GeV, this cone size falls linearly with p T to a minimum cone size of 0.2 at 50 GeV. Similarly, the scalar E T sum (calorimeter isolation) is calculated from the positive-energy topological clusters that are not associated with a lepton track in a cone of R = 0.2 around the muon or electron. The sum of the track isolation and 40% of the calorimeter isolation is required to be less than 16% of the lepton p T . The calorimeter isolation is corrected for electron shower leakage, pile-up and underlying-event contributions. Both isolations are corrected for track and topological cluster contributions from the remaining three leptons. The pileup dependence of this isolation selection is improved compared with that of the previous measurements [17,128,146] by optimising the criteria used for exclusion of tracks associated with a vertex other than the primary vertex and by the removal of topological clusters associated with tracks. The signal efficiency of the isolation criteria is greater than 80%, improving the efficiency by about 5% compared with the previous analysis for the same background rejection.
The four quadruplet leptons are required to originate from a common vertex point. A requirement corresponding to a signal efficiency of better than 99.5% is imposed on the χ 2 value from the fit of the four lepton tracks to their common vertex.
If there is more than one decay final state per event with the priority quadruplet (m 12 closest to m Z ) satisfying the selection criteria, the quadruplet from the final state with highest selection efficiency, i.e. ordered 4μ, 2e2μ, 2μ2e and 4e, is chosen as the Higgs boson candidate.
In the case of VH or tt H production, there may be additional prompt leptons present in the event, together with the selected quadruplet. Therefore, there is a possibility that one or more of the leptons selected in the quadruplet do not originate from a Higgs boson decay, but rather from the V boson leptonic decay or the top quark semileptonic decay. To improve the lepton pairing in such cases, a matrixelement-based pairing method assuming the SM tensor structure is used for all events containing at least one additional lepton with p T > 12 GeV and satisfying the same identification, isolation and angular separation criteria as the four quadruplet leptons [17,128]. For all possible quadruplet combinations that satisfy the selection, a matrix element for the Higgs boson decay is computed at LO using the MadGraph5_aMC@NLO [55] generator, with the reconstructed lepton momentum vectors as inputs to the calculation. The quadruplet with the largest matrix-element value is selected as the Higgs boson candidate. This method leads to a 50% improvement in correctly identifying the leptons in the quadruplet as those originating from a Higgs boson decay if an extra lepton is identified. The impact of the matrix element on the expected invariant mass distribution is shown in Fig. 2a.
To improve the four-lepton invariant mass reconstruction, the reconstructed final-state radiation (FSR) photons in Z boson decays are accounted for using the same strategy as the previous publications [17,128]. Collinear FSR candidates are defined as candidates with R < 0.15 to the nearest lepton in the quadruplet. Collinear FSR candidates are considered only for muons from the leading lepton pair, while non-collinear FSR candidates are considered for both muons and electrons from leading and subleading Z bosons.
Collinear FSR candidates are selected from reconstructed photon candidates and from electron candidates that share an ID track with the muon. Further criteria are applied to each candidate, based on the following discriminants: the fraction, f 1 , of cluster energy in the front segment of the EM calorimeter divided by the total cluster energy to reduce backgrounds from muon ionisation; the angular distance, R cluster,μ , between the candidate EM cluster and the muon; and the candidate p T , which must be at least 1 GeV. For all selected electron candidates and for photon candidates with p T < 3.5 GeV, a requirement of f 1 > 0.2 and R cluster,μ < 0.08 is imposed. The collinear photon candidates with p T > 3.5 GeV are selected if f 1 > 0.1 and R cluster,μ < 0.15. Non-collinear FSR candidates are selected only from reconstructed isolated photons meeting the 'tight' criteria [129,147] and satisfying p T > 10 GeV and R cluster, > 0.15. Only one FSR candidate is included in the quadruplet, with preference given to collinear FSR and to the candidate with the highest p T . An FSR candidate is added to the lepton pair if the invariant mass of the lepton pair is between 66 and 89 GeV and if the invariant mass of the lepton pair and the photon is below 100 GeV. Approximately 3% of reconstructed Higgs boson candidates have an FSR candidate and its impact on the expected invariant mass distribution is shown in Fig. 2b. The Higgs boson candidates within a mass window of 115 GeV < m 4 < 130 GeV are selected as the signal region. Events failing this requirement but that are within a mass window of 105 GeV < m 4 < 115 GeV or 130 GeV < m 4 < 160 (350) GeV are assigned to the sideband regions used to estimate the leading backgrounds as described in Sect. 6.
The selection efficiencies of the simulated signal in the fiducial region |y H | < 2.5, where y H is the Higgs boson rapidity, are about 33%, 25%, 19% and 16%, in the 4μ, 2e2μ, 2μ2e and 4e final states, respectively.

Event categorisation and production mode discrimination
In order to be sensitive to different production bins in the framework of simplified template cross-sections, the selected Higgs boson candidates in the mass window 115 GeV < m 4 < 130 GeV are classified into several dedicated reconstructed event categories. In addition, the events in the mass sidebands are also categorised for purposes of background estimation described in Sect. 6. In general, more than one production mode contributes to each reconstructed event category, as well as various background processes. For this reason, multivariate discriminants are introduced in most of the mutually exclusive reconstructed event categories to distinguish between these contributions.

Event categorisation
For signal events, the classification is performed in the order shown in the middle-right panel of Fig. 1 (from bottom to top) and as described below. First, those events classified as enriched in the tt H process are split according to the decay mode of the two W bosons from the top quark decays. For semileptonic and dileptonic decays (ttH-Lep-enriched), at least one additional lepton with p T > 12 GeV 3 together with at least two b-tagged jets (with 85% b-tagging efficiency), or at least five jets among which at least one b-tagged jet (with 85% b-tagging efficiency) or at least two jets among which at least one b-tagged jet (with 60% b-tagging efficiency) is required. For the fully hadronic decay (ttH-Hadenriched), there must be either at least five jets among which at least two b-tagged jets (with 85% b-tagging efficiency) or at least four jets among which at least one b-tagged jet (with 60% b-tagging efficiency). Events with additional leptons but not satisfying the jet requirements define the next category enriched in VH production events with leptonic vector-boson decay (VH-Lep-enriched).
The remaining events are classified according to their reconstructed jet multiplicity into events with no jets, exactly one jet or at least two jets. Events with at least two reconstructed jets are divided into two categories: one is a 'BSMlike' category (2 j-BSM-like) and the other (2 j) contains the bulk of events with significant contributions from the VBF and VH production modes in addition to ggF. The 2 j-BSMlike category requires the invariant mass m j j of the two leading jets to be larger than 120 GeV and the four-lepton transverse momentum, p 4 T , to be larger than 200 GeV; the remaining events are placed in the 2 j category.
Events with zero or one jet in the final state are expected to be mostly from the ggF process. Following the particle-level definition of production bins in Sect. 1.1, the 1-jet category is further split into four categories with p 4 T smaller than 60 GeV (1 j-p 4 T -Low), between 60 and 120 GeV (1 j-p 4 T -Med), between 120 and 200 GeV (1 j-p 4 T -High), and larger than 200 GeV (1 j-p 4 T -BSM-like). The largest number of ggF events and the highest ggF purity are expected in the zero-jet category. The zero-jet category is split into three categories with p 4 T smaller than 10 GeV (0 j-p 4 T -Low), between 10 and 100 GeV (0 j-p 4 T -Med) and above 100 GeV (0 j-p 4 T -High). The first two categories follow the production bin splitting, and the last category improves the discrimination between VH (V → ν/νν) and ggF.
As illustrated in Fig. 1, there is a dedicated reconstructed event category for each production bin except for gg2H-2 j, qq2Hqq-VH and qq2Hqq-VBF. These production bins are largely measured from the 2-jet reconstruction category, and to a lesser extent from the 1-jet categories, using multivariate discriminants (see Sect. 5.2). The gg2Hp H T -High production bin is measured simultaneously in all reconstructed event categories with high transverse momentum of the four-lepton system, independent of the reconstructed jet multiplicity.
The rightmost panel of Fig. 1 shows the background event classification. For estimating the tXX process from the mass sideband, a tXX-enriched sideband category (SBt X X-enriched) is defined, which includes events with at least two jets including at least one tagged as a b-jet with 60% efficiency and E miss T > 100 GeV in the m 4 mass range 105-115 GeV or 130-350 GeV. This region is dominated by tt Z (87%) and has small contributions from tt, tttt, t W Z, tt W , tt W W , tt W Z, tt Zγ , tt Z Z and t Z. The tXX process is expected to give the largest contribution in 'tt H -like' categories. The large mass range for this category, larger than for the non-resonant Z Z as discussed next, allows better statistical precision for the estimate of this background.
For the estimation of non-resonant Z Z * production, events not meeting the criteria for the SB-t X X-enriched category and in the m 4 mass range 105-115 GeV or 130-160 GeV are split according to the number of reconstructed jets: exactly zero jets (SB-0 j), exactly one jet (SB-1 j) or at  final-state radiation for candidates with an FSR candidate. For (a), the overflow events are included in the last bin least two jets (SB-2 j). This mass range limits the contribution from the single-resonance process, Z → 4 , and from the on-shell Z Z process. Similarly, events in the same mass range with an extra reconstructed lepton separately form the SB-VH-Lep-enriched category, which is enriched with signal events containing leptons from the associated V leptonic decay or the top quark semileptonic decay. This category is mainly designed to improve the expected sensitivity for VH-Lep by about 5%, having a VH purity of about 19%.
The expected number of signal events is shown in Table 4 for each reconstructed event category separately for each production mode. The ggF and bbH contributions are shown separately to compare their relative contributions, but both belong in the same (ggF) production bin. The highest bbH event yield is expected in the 0 j categories since the jets tend to be more forward than in the tt H process, thus escaping the acceptance of the tt H selection criteria. The sources of uncertainty in these expectations are detailed in Sect. 7. The signal composition in terms of the Reduced Stage-1.1 production bins is shown in Fig. 3.
The separation of the contributions from different production bins, such as the gg2H-2 j, qq2Hqq-VH and qq2Hqq-VBF components contributing in categories with two or more jets, is improved by means of discriminants obtained using multivariate data analysis, as described in the following section.

Multivariate production mode discriminants
To further increase the sensitivity of the cross-section measurements in the production bins (Sect. 1.1), multivariate discriminants using neural networks (NNs) [148] are introduced in many of the reconstructed signal event categories as observables used in the statistical fit, described in Sect. 8.2. The NN architecture and training procedure are defined using Keras with TensorFlow [149,150]. These networks are trained using several discriminating observables, as defined in Table 5, on simulated SM Higgs boson signals with m H = 125 GeV or non-Higgs-boson background. Due to the low number of signal events expected in the 0 j-p 4 T -High, 1 j-p 4 T -BSM-like and ttH-Lep-enriched categories, only the observed yield is used as the discriminant in these categories.
Two types of NNs are used: feed-forward multilayer perceptron (MLP) and recurrent (RNN) [148][149][150][151][152]. Each NN discriminant combines two RNNs, one for the p T -ordered variables related to the four leptons in the quadruplet and one for variables related to jets, and an MLP with additional variables related to the full event. The jet RNN accepts inputs from up to three jets. The outputs of the MLP and the two RNNs are chained into another MLP to complete an NN discriminant, which is trained to approximate the posterior probability for an event to originate from a given process. This is used in each reconstructed event category to discriminate between two or three processes, e.g. ggF, VBF and Z Z background in the 1 j-p 4 T -Low category. The variables used to train the MLP and RNNs for each category along with the processes being separated are summarised in Table 5.
The NN training variables not previously defined are listed as follows. The kinematic discriminant D Z Z * [153], defined as the difference between the logarithms of the squared matrix elements for the signal decay (same as in Sect. 4) and squared matrix elements for the background process, is used Table 4 The Expected Composition to distinguish ggF from the non-resonant Z Z background. Three angles [7] are used to further distinguish these processes: the cosine of the leading Z boson's production angle θ * in the four-lepton rest frame; the cosine of θ 1 defined as the angle between the negatively charged lepton of the leading Z in the leading Z rest frame and the direction of flight of the leading Z in the four-lepton rest frame; and the angle φ Z Z , between the two Z decay planes in the four-lepton rest frame. Page 15 of 54 957 The angular separation of the leading jet from the 4 system, R 4 j , is used to distinguish VBF or tt H from ggF. For categories with two or more jets, kinematic variables that also include the information from the two leading jets are used: the invariant mass, m j j ; the transverse momentum of the 4 and the 2-jet system, p 4 j j T ; and the Zeppenfeld variable, η [154]. The number of reconstructed jets, N jets , the number of b-tagged jets at 70% tagging efficiency, N b-jets,70% , and the scalar sum of the p T of all reconstructed jets, H T , are used to identify the tt H process.
Depending on the category and the number of processes being targeted, the NN has two or three output nodes. The value computed at each node represents the probability, with an integral of one, for the event to originate from the given process. For example, for the 0-jet category, two probabilities are evaluated, NN ggF and NN Z Z . As these two values are a linear transformation of each other, only one output, NN ggF , is used as a discriminant in the fit model. In categories with three targeted processes, only two of the three corresponding output probabilities are independent. In a given category, a selection is applied on one of the three output probabilities to split the events in two subcategories. This output probability is then used as the discriminant for the subcategory of events passing the selection, while for the other subcategory one of the two remaining output probabilities is used. The selection criterion is chosen so as to provide the largest purity of the targeted process for events passing the selection. For example, in the 1-jet category, NN VBF and NN Z Z are used. The subcategory of events with NN Z Z larger than 0.25 uses NN Z Z as the discriminant in the fit model, while NN VBF is used in the remaining subcategory. The subcategory definitions and observables used in all reconstructed event categories are summarised in Table 5.

Background processes with prompt leptons
Non-resonant SM Z Z * production via qq annihilation, gluon-gluon fusion and vector-boson scattering can result in four prompt leptons in the final state and constitutes the largest background for the analysis. While for the previous analyses [17,128], simulation was exclusively used to estimate both the shape and normalisation, in this analysis the normalisation is constrained by a data-driven technique. This allows the systematic uncertainty to be reduced by removing both the theoretical and luminosity uncertainties contributing to the normalisation uncertainty.
As outlined in Sect. 5.1, to estimate the normalisation, sideband categories in the m 4 mass region 105-115 GeV and 130-160 GeV are defined according to the jet multiplicity (SB-0 j, SB-1 j, SB-2 j). The normalisation of the Z Z * background is simultaneously fitted with a common normalisation factor for signal region and sideband categories with the same jet multiplicity. For example, the Z Z * background is scaled by a common factor for 2 j, 2 j-BSM-like and SB-2 j categories. The background shape templates for NN discriminants and the expected fraction of events in relevant reconstructed signal-region event categories are obtained from simulation. As shown in Fig. 4a, good agreement is found between data and simulation for the shape of the NN observable. All expected distributions are shown after the final fit to the data for the Production Mode measurement (see Sect. 8) and are referred to as post-fit distributions in the following. The simulated distributions of the observables p 4 T and m j j employed for the prediction of event fractions in each event category also agree with data, as seen in Figs. 4b, c respectively. The estimation of the Z Z * process in the jet multiplicity bins removes one of the leading theoretical uncertainties [155]. Due to the limited sensitivity and the low expected yield, the normalisation of Z Z * in tt H -like categories is estimated from simulation.
Similarly, backgrounds affecting the tt H -like categories are estimated simultaneously from an enriched sample selected in a dedicated sideband region (SB-t X X-enriched), with the mass cut extended up to 350 GeV to improve the statistical precision of the estimate. The normalisation of the t X X process is simultaneously fitted across the ttH-Lep-enriched, ttH-Had-enriched and SB-t X X-enriched categories. The N jets observable distribution, which is used to predict the event fractions in each category, is shown in Fig. 4d and agrees with data. In all other categories, the sensitivity of the t X X measurement is limited due to a small number of expected t X X events and its normalisation is estimated from simulation.
The contribution from VVV processes is estimated for all categories using the simulated samples presented in Sect. 3.

Background processes with non-prompt leptons
Other processes, such as Z + jets, tt, and W Z, containing at least one jet, photon or lepton from a hadron decay that is misidentified as a prompt lepton, also contribute to the background. These 'reducible' backgrounds are significantly smaller than the non-resonant Z Z * background and are estimated from data using different approaches for the + μμ and + ee final states [17,128].
In the + μμ final states, the normalisation of the Z + jets and tt backgrounds are determined by performing fits to the invariant mass of the leading lepton pair in dedicated independent control regions. The shape of the invariant mass distribution for each region is parameterised using simulated samples. In contrast to the previous analyses [17,128], this fit is performed independently for each reconstructed event  The first two are the primary control regions used to estimate Z + jets and tt, and the latter two improve the estimate by reducing the statistical error of the fitted normalisation.
The background normalisations are obtained separately for the Z + jets and tt background processes using the simultaneous fit in the four control regions. The normalisation n C R i in each control region CR for the background process i is expressed as a fraction, n C R The + ee control-region selection requires the electrons in the subleading lepton pair to have the same charge, and relaxes the identification, impact parameter and isolation requirements on the electron candidate with the lowest transverse energy. This fake electron candidate, denoted by X , can be a light-flavour jet, an electron from photon conversion or an electron from heavy-flavour hadron decay. The heavy-flavour background is determined from simulation. Good agreement is observed between simulation and data in a heavy-flavour enriched control region.
The remaining background is separated into light-flavour and photon conversion background components using the sPlot method [156] which is performed on electron candidates X , separately for each reconstructed category in bins of the jet multiplicity and the transverse momentum of the electron candidate. The size of the two background components is obtained from a fit to the number of hits from the electron candidate X in the innermost ID layer in the + ee data control region, where a hit indicates either a hadron track or an early conversion. A hit in the next-to-innermost pixel layer is used when the electron falls in a region that was either not instrumented with an innermost pixel layer module or where the module was not operating. The templates of the final discriminants for the mentioned fit of the light-flavour and photon conversion background components are obtained from simulated Z + X events with an on-shell Z boson decay candidate accompanied by an electron X selected using the same criteria as in the + ee control region. The simulated Z + X events are also used to obtain the transfer factor for the X candidate for the extrapolation of the light-flavour and photon conversion background contributions from the + ee control region to the signal region, after correcting the simulation to match the data in dedicated control samples of Z + X events. The extrapolation to the signal region is also performed in bins of the electron transverse momentum and the jet multiplicity, separately for each reconstructed event category. A method similar to that for the + μμ final state is used to extract the NN shape, where the fractions of events from light-flavour jets and photon conversions are estimated from simulation and corrected transfer factors are used.

Systematic uncertainties
The systematic uncertainties are categorised into experimental and theoretical uncertainties. The first category includes uncertainties in lepton and jet reconstruction, identification, isolation and trigger efficiencies, energy resolution and scale, and uncertainty in the total integrated luminosity. Uncertainties from the procedure used to derive the data-driven background estimates are also included in this category. The second category includes uncertainties in theoretical modelling of the signal and background processes.
The uncertainties can affect the signal acceptance, selection efficiency and discriminant distributions as well as the background estimates. The dominant sources of uncertainty and their effect are described in the following subsections. The impact of these uncertainties on the measurements is summarised in Table 6.

Experimental uncertainties
The uncertainty in the combined 2015-2018 integrated luminosity is 1.7% [157], obtained using the LUCID-2 detector [158] for the primary luminosity measurements. This uncertainty affects the signal and the normalisation of the simulated background estimates when not constrained by the data sidebands.
The uncertainty in the predicted yields due to pile-up modelling ranges between 1% and 2% and is derived by varying the average number of pile-up events in the simulation to cover the uncertainty in the ratio of the predicted to measured inelastic cross-sections [159].
The electron (muon) reconstruction, isolation and identification efficiencies, and the energy (momentum) scale and resolution are derived from data using large samples of J/ψ → and Z → decays [129,131]. Typical uncertainties in the predicted yields for the relevant decay channels Table 6 The impact of the dominant systematic uncertainties (in percent) on the cross-sections in production bins of the Production Mode Stage and the Reduced Stage 1.1. Similar sources of systematic uncertainties are grouped together: luminosity (Lumi.), electron/muon reconstruction and identification efficiencies and pile-up modelling (e, μ, pile-up), jet energy scale/resolution and b-tagging efficiencies (Jets, flav. tag), uncertainties in reducible background (reducible bkg), theo-retical uncertainties in Z Z * background and tXX background, and theoretical uncertainties in the signal due to parton distribution function (PDF), QCD scale (QCD) and parton showering algorithm (Shower). The uncertainties are rounded to the nearest 0.5%, except for the luminosity uncertainty, which is measured to be 1.7% and increases for the VH signal processes due to the simulation-based normalisation of the Reduced Stage-1.1 production bin cross-sections due to the identification and reconstruction efficiency uncertainties are below 1% for muons and 1%-2% for electrons. The uncertainty in the expected yields due to the muon and electron isolation efficiency is also taken into account, with the typical size being 1%. The uncertainties in the trigger efficiencies have a negligible impact. The uncertainties in the electron and muon energy and momentum scale and resolution are small and also have a negligible impact on the measurements. The uncertainties in the jet energy scale and resolution are in the range 1%-3% [133]. The impact of these uncertainties is more relevant for the VH, VBF and tt H production mode cross-sections (3%-5%) and for all the Reduced Stage-1.1 cross-section measurements, including the ggF process split into the different N jets exclusive production bins (5%-20%).
The uncertainty in the calibration of the b-tagging algorithm, which is derived from dileptonic tt events, amounts to a few percent over most of the jet p T range [138]. This uncertainty is only relevant in the ttH category, with its expected impact being approximately 1% in the tt H cross-section measurement. The uncertainties associated with the E miss T reconstruction have a negligible impact.
A shift in the simulated Higgs boson mass corresponding to the precision of the Higgs boson measurement, m H = 125.09±0.24 GeV [160], is shown to have a negligible impact on the signal acceptance. A small dependency of the NN ggF discriminant shape in the 0 j-p 4 T -Low and 0 j-p 4 T -Med categories on m H is observed for the signal (below 2% in the highest NN score bins) and is included in the signal model. This uncertainty affects the measurement of ggF production, as well as the measurements in other production bins with large ggF contamination.
For the data-driven measurement of the reducible background, three sources of uncertainty are considered: statistical uncertainty, overall systematic uncertainty for each of + μμ and + ee, and a shape systematic uncertainty that varies with the reconstructed event category. Since the yields are estimated by using a statistical fit to a control data region with large statistics, the inclusive background estimate has a relatively small (3%) statistical uncertainty. The systematic uncertainty for + μμ and the heavy-flavour component of + ee is estimated by comparing the lepton identification, isolation and impact parameter significance efficiency between data and simulated events in a separate region, enriched with on-shell Z boson decays accompanied by an electron or a muon. For both the + μμ and + ee estimates, the difference in efficiency is assigned as the uncertainty in the extrapolation of the yield estimate from the control region to the signal region. For the + ee light-flavour component, the efficiency is derived from an enriched control region with a systematic uncertainty estimated by varying the assumed light-and heavy-flavour components. These inclusive uncertainties (6%) are treated as correlated across the reconstructed event categories. Finally, there are additional uncorrelated uncertainties (8%-70%) in the fraction of the reducible background in each event category due to the statistical precision of the simulated samples.

Theoretical uncertainties
The theoretical modelling of the signal and background processes is affected by uncertainties due to missing higherorder corrections, modelling of parton showers and the underlying event, and PDF+α S uncertainties.
The impact of the theory systematic uncertainties on the signal depends on the kind of measurement that is performed. For signal-strength measurements, defined as the measured cross-section divided by the SM prediction, or interpretation of cross-section using the EFT approach, each source of theory uncertainty affects both the acceptance and the predicted SM cross-section. For the cross-section measurements, only effects on the acceptance need to be considered.
The impact of the theory systematic uncertainties on the background depends on the method of estimating the normalisation. If simulation is used, the uncertainties in the acceptance and the predicted SM cross-section are included. If the normalisation is estimated from a data-driven method, only the impact on the relative event fractions between categories is considered.
One of the dominant sources of theoretical uncertainty is the prediction of the ggF process in the different N jets categories. The ggF process gives a large contribution in categories with at least two jets. To estimate the variations due to the impact of higher-order contributions not included in the calculations and migration effects on the N jets ggF crosssections, the approach described in Refs. [24,161] is used, which exploits the latest predictions for the inclusive jet cross-sections. In particular, the uncertainty from the choice of factorisation and renormalisation scales, the choice of resummation scales, and the migrations between the 0-jet and 1-jet phase-space bins or between the 1-jet and ≥ 2-jet bins are considered [24,[162][163][164]. The impact of QCD scale variations on the Higgs boson p T distribution is taken into account as an additional uncertainty. The uncertainty in higher-order corrections to the Higgs boson p T originating from the assumption of infinite top quark mass in the heavyquark loop is also taken into account by comparing the p T distribution predictions to finite-mass calculations. An additional uncertainty in the acceptance of the ggF process in VBF topologies [165] due to missing higher orders in QCD in the calculation is estimated by variations of the renormalisation and factorisation scales using fixed-order calculations with MCFM [166]. An additional uncertainty in the Higgs boson p T distribution, derived by varying the renormalisation, factorisation and NNLOPS scale in the simulation, in the 0-jet topology is considered. This is particularly relevant when measuring the inclusive ggF cross-section using the p 4 T categories for events with no jet activity. To account for higher-order corrections to p H j j T , which is used as an NN input variable, the uncertainty is derived by comparing the predicted distribution obtained using Powheg NNLOPS and MadGraph5_aMC@NLO with the FxFx merging scheme.
For the VBF production mode, the uncertainty due to missing higher orders in QCD is parameterised using the scheme outlined in Ref. [23]. The migration effects due to the selection criteria imposed on the number of jets, transverse momentum of the Higgs boson, transverse momentum of the Higgs boson and the leading dijet system and the invariant mass of the two leading jets, used to define the full Stage 1.1 STXS production bins, are computed by varying the renormalisation and factorisation scales by a factor of two. The uncertainties are cross-checked with fixed-order calculations. Similarly, for the VH production mode with the associated V decaying leptonically, the scale variations are parameterised as migration effects due to the selection criteria imposed on the number of jets and the transverse momentum of the associated boson [167].
For the VH production mode with the associated V decaying hadronically and the tt H production mode, the uncertainty due to missing higher orders in QCD is obtained by varying the renormalisation and factorisation scales by a factor of two. The configuration with the largest impact, as quantified by the relative difference between the varied and the nominal configuration, is chosen to define the uncertainty in each experimental category. These uncertainties are treated as uncorrelated among the different production modes. Due to the limited accuracy of the simulated samples, the uncertainties evaluated using this method for the total cross-sections are larger than those described in Ref. [24].
The uncertainties in the acceptance due to the modelling of parton showers and the underlying event are estimated with AZNLO tune eigenvector variations and by comparing the acceptance using the parton showering algorithm from Pythia 8 with that from Herwig 7 [168] for all signal processes. The uncertainty due to each AZNLO tune variation is taken as correlated among the different production modes while the difference between the parton showering algorithms is treated as an uncorrelated uncertainty. The uncertainties due to higher-order corrections to the Higgs boson decay are modelled using the PROPHECY4F [102,105] and Hto4L [104,169] generators. These corrections are below 2% and have a negligible impact on the results. A 100% uncertainty is assigned to heavy-flavour quark production modelling for the ggF contribution entering in the ttH category. This has a negligible impact on the results.
The impact of the PDF uncertainty is estimated with the thirty eigenvector variations of the PDF4LHC_nlo_30 Hessian PDF set following the PDF4LHC recommendations [40]. The modification of the predictions originating from each eigenvector variation is added as a separate source of uncertainty in the model. The same procedure is applied for the ggF, VBF, VH and tt H processes, enabling correlations to be taken into account in the fit model.
The impacts of the theoretical uncertainties, as described above, on the shape of NN discriminants are also considered. For ggF production, a further cross-check is performed by comparing the NN shapes in the corresponding categories as predicted by Powheg NNLOPS and Mad-Graph5_aMC@NLO with the FxFx merging scheme. All the NN shapes from the two generators agree within the scale variations and, therefore, no additional shape uncertainty is included.
For signal-strength measurements, an additional uncertainty related to the H → Z Z * branching ratio prediction [102,105] is included in the measurement.
Since the normalisation of the Z Z * process in most reconstructed event categories is constrained by performing a simultaneous fit to sideband regions enriched in this contribution together with the signal regions, most of the theoretical uncertainty in the normalisation for this background vanishes. Nevertheless, uncertainties in the shapes of the discriminants for the Z Z * background and in the relative contribution of this background between the sidebands and the signal regions are taken into account. The uncertainties due to missing higher-order effects in QCD are estimated by varying the factorisation and renormalisation QCD scales by a factor of two; the impact of the PDF uncertainty is estimated by using the MC replicas of the NNPDF3.0 PDF set. Uncertainties due to parton shower modelling for the Z Z * process are considered as well. The impact of these uncertainties is below 2% for all production mode cross-sections measured. In addition, a comparison between Sherpa and Powheg is also taken as an additional source of systematic uncertainty. This model uncertainty is treated as uncorrelated among the different sideband-to-signal region extrapolations (in 0-jet, 1-jet and 2-jet categories).
The uncertainty in the gluon-initiated and the vectorboson-initiated Z Z * process is taken into account by changing the relative composition of the quark-initiated, the gluoninitiated and the vector-boson-scattered Z Z * components according to the theoretical uncertainty in the predicted crosssections and the respective K -factors. In addition, the event yield and NN discriminant shapes in each event category are compared with the data in an m 4 sideband around the signal region (105 GeV< m 4 < 115 GeV or 130 GeV< m 4 < 160 GeV). Good agreement between the Sherpa predictions and the data is found.
For the tXX process, uncertainties due to PDF and QCD scale variations are considered in the relative fraction of events present in the tt H -like categories, in the SBt X X-enriched control region and in the NN discriminant shape. Differences between MadGraph5_aMC@NLO and Sherpa are considered as an additional systematic uncertainty. For all other categories where this process is estimated from simulation, the impact of these uncertainties on the SM cross-section and acceptance are also considered.
Uncertainties in the PDF and in missing higher-order corrections in QCD are applied to the VVV background estimate, which is fully taken from MC simulation.
To probe the tensor structure of the Higgs boson coupling in the EFT approach, theoretical uncertainties due to PDF and QCD scale variations are assigned to the signal predictions based on the simulated highest-order SM signal samples. The same uncertainties are assigned to all corresponding BSM signal predictions, since it is shown using the MC signal samples simulated at LO accuracy that the uncertainties change negligibly as a function of the Wilson coefficients.

Observed data
The expected and observed four-lepton invariant mass (postfit) distributions of the selected Higgs boson candidates after the event selection are shown in Fig. 5. The observed and expected (post-fit) distributions of the jet multiplicity, the dijet invariant mass, and the four-lepton transverse momenta in different N jets bins, which are used for the categorisation of reconstructed events, are shown in Fig. 6 for different steps of the event categorisation.
The expected numbers of signal and background events in each reconstructed event category are shown in Table 7 together with the corresponding observed number of events. The expected event yields are in good agreement with the observed ones. The observed and expected (post-fit) distributions of the NN discriminants are shown in Fig. 7 and in Figure 8. In addition, Fig. 8g  is used and in the mass sidebands used to constrain the Z Z * and t X X background, respectively. All distributions are in good agreement with the data.
The statistical interpretation of the results and compatibility with the SM are discussed in the following.

Measurement of simplified template cross-sections
To measure the product σ ·B of the Higgs boson production cross-section and the branching ratio for H → Z Z * decay for the production bins of the Production Mode Stage or the Reduced Stage 1.1, a fit to the discriminant observables introduced in Sect. 5.2 is performed using the likelihood function L( σ , θ) that depends on the Higgs boson production cross-section σ = {σ 1 , σ 2 , . . . , σ N } where σ p is the cross-section in each production bin p and the nuisance parameters θ accounting for the systematic uncertainties. The likelihood function is defined as a product of conditional probabilities over binned distributions of the discriminating observables in each reconstructed signal and sideband event category j, with Poisson distributions P corresponding to the observation of N i, j events in each histogram bin i of the discriminating observable given the expectations for each background process, B i, j ( θ), and for the signal, S i, j ( θ) = L · σ · B · A i, j ( θ), where L is the integrated luminosity and A i, j = {A 1 i, j , A 2 i, j , . . . , A N i, j } is the set of signal acceptances from each production bin. The signal acceptance A p i, j is defined as the fraction of generated signal events in the production bin p that satisfy the event reconstruction and selection criteria in the histogram bin i of the reconstructed event category j. For a given production bin p, the acceptance consists of A p i, j = a p · p i, j , where a p is the particle-level acceptance in the fiducial region defined from requirements listed in Sects. 4 and 5 and p i, j is the reconstruction efficiency of these particle-level events. Constraints on the nuisance parameters corresponding to systematic uncertainties described in Sect. 7 are represented by the functions C m ( θ). The cross-sections are treated as independent parameters for each production bin and correlated among the different reconstructed event categories. The test statistic used to perform the measurements is the ratio of profile likelihoods [170], where σ represents only the cross-section(s) considered as parameter(s) of interest in a given fit. The likelihood in the numerator is the estimator of a conditional fit, i.e. with parameter(s) of interest σ i fixed to a given value, while the remaining cross-sections and nuisance parameters are free-floating parameters in the fit. The values of the nuisance parame-tersˆ θ( σ )) maximise the likelihood on the condition that the parameters of interest are held fixed to a given value. The likelihood in the denominator is the estimator of an unconditional fit in which all σ and θ parameters are free parameters of the fit. The parameter of interest σ in each production bin is alternatively replaced by μ · σ SM ( θ), allowing an interpretation in terms of the signal strength μ relative to the SM prediction σ SM ( θ).
Assuming that the relative signal fractions in each production bin are given by the predictions for the SM Higgs boson, the inclusive H → Z Z * production cross-section for |y H | < 2.5 is measured to be: The measured cross-section and signal strength are in an excellent agreement with the SM prediction, with a p-value of 98.6% for both compatibility tests.
The corresponding likelihood functions are shown in Fig. 9. The dominant systematic uncertainty in the crosssection measurement is the experimental uncertainty in the lepton efficiency and integrated luminosity measurements and theoretical uncertainties related to parton shower mod- elling affecting the acceptance. The signal-strength measurement is also affected by the theoretical uncertainty in the ggF cross-section due to missing higher-order corrections in QCD. The expected SM cross-section, the observed values of σ · B(H → Z Z * ) and their ratio for the inclusive production and in each production bin of the Production Mode Stage and the Reduced Stage 1.1 are shown in Table 8.
The corresponding values are summarised in Fig. 10. In the ratio calculation, uncertainties in the SM expectation are not taken into account. The Production Mode Stage and Reduced Stage-1.1 measurements agree with the predictions for the SM Higgs boson. The p-values of the corresponding compatibility tests are 91% and 77%, respectively.
For the qq2Hqq-VBF bin, most of the sensitivity to the VBF production mode comes from the phase space with m j j > 350 GeV and p H T < 200 GeV. To probe the VBF contribution more directly, the cross-sections in this and in the remaining phase space region of the qq2Hqq-VBF bin are fitted separately to the data, simultaneously with the other Reduced Stage 1.1 bins, using the reconstruction categories described in Sect. 5. The cross-section in the m j j > 350 GeV and p H T < 200 GeV phase space is measured to be 0.060 +0.025 −0.020 pb compared with the predicted cross-section of 0.0335 +0.0007 −0.0011 pb. This measurement has a correlation of 20% with the measurement in the gg2H-2 j bin, while correlations with other bins are up to 50%.
The dominant contribution to the measurement uncertainty in the ggF Production Mode Stage bin originates from the same sources as in the inclusive measurement. For the VBF production bin, the dominant systematic uncertainties are related to parton showering modelling and jet energy scale and resolution uncertainties. The VBF, VH and tt H production bins are also affected by the theoretical uncertainties related to the modelling of the ggF process. For the Reduced Stage-1.1 bins, the dominant cross-section uncertainties are the jet energy scale and resolution, and parton shower uncertainties.

Constraints on the Higgs boson couplings in the κ-framework
The cross-sections measured at the Production Mode Stage are interpreted in the κ-framework described in Sect.
It is assumed that there are no undetected or invisible Higgs boson decays. The observed likelihood contours in the κ V -κ F plane are shown in Fig. 12 (only the quadrant κ F > 0 and κ V > 0 is shown since this channel is not sensitive to the relative sign of the two coupling modifiers). The best-fit value isκ V = 1.02 ± 0.06 and κ F = 0.88 ± 0.16, with the correlation of −0.17. The probability of compatibility with the Standard Model expectation is at the level of 75%.

Constraints on the tensor coupling structure in the EFT approach
To interpret the observed data in the framework of an effective field theory, an EFT signal model is built by parameterising the production cross-sections in each production bin of the Reduced Stage 1.1, as well as the branching ratio and the signal acceptances, as a function of the SMEFT Wilson coefficients introduced in Sect. 1.3. The constraints on the Wilson coefficients are then obtained from the simultaneous fit to the data in all reconstructed signal and sideband event categories. Due to the statistical precision of the data sample, the constraints are always set on one or at most two of the Wilson coefficients at a time, while the values of the remaining coefficients are assumed to be equal to zero.

EFT signal model
The EFT parameterisation of the production cross-sections in each production bin of the Reduced Stage 1.1 is obtained from Eq. (1) using simulated BSM samples introduced in Sect. 3. The contribution from the gg → Z (→ )H process is taken from the SM simulation and assumed to scale with BSM parameters in the same way as the qq → Z (→ )H processes. As in the case of simplified template cross-section measurements, tt H and t H processes are combined into a single tt H production bin. The cut-off scale is set to = 1 TeV. Only LO computation of QCD and SM electroweak processes is provided, with LO effective couplings for the SM Higgs boson to gluon and to photon vertices. An assumption is made that higher-order corrections, applied in a multiplicative way, are the same for both the SM and the BSM LO predictions and therefore no changes in the parameterisation are expected due to higher-order effects [171]. With the current amount of data, the constraints from the VBF, VH and tt H production modes on the relevant Wilson coefficients still allow a rather large range of parameter values in which the quadratic term (the last term in Eq. (1)) cannot be neglected even though its contribution is suppressed by 4 . Such dimension-six quadratic terms are therefore included in the EFT parameterisation. Since the linear terms from dimension-eight operators are suppressed by the same factor, they could in general also give similar non-negligible contributions. Dimension-eight terms are currently not available in the SMEFT model and are thus not taken into account.
The branching ratio for the H →Z Z * →4 decay is parameterised in terms of Wilson coefficients following Eq. (2). The partial and total decay widths are calculated in Mad-Graph5_aMC@NLO. The total decay width is calculated by taking into account the dominant Higgs boson decay modes: γ γ , Z γ , bb, gg, W W and Z Z. Other decay modes are not affected by the probed Wilson coefficients. Their contribution to the total decay width is therefore given by the corresponding SM predictions.
The selection criteria for the four-lepton Higgs boson candidates, in particular the requirements on the minimum invariant mass m 34 of the subleading lepton pair, introduce an additional dependence of the signal acceptance on the BSM coupling parameters. The particle-level signal acceptance A, defined as the fraction of signal events satisfying the Higgs boson candidate selection criteria applied at particle-level, has therefore been simultaneously parameterised in terms of the three

Fig. 12
Likelihood contours at 68% CL (dashed line) and 95% CL (solid line) in the κ V -κ F plane. The best fit to the data (solid cross) and the SM prediction (star) are also indicated even) parameters vanish. The dependence of the acceptance on other EFT coupling parameters is shown to be negligible as these parameters have negligible or no impact on the H → Z Z * decay. The acceptance correction relative to the SM prediction is described by a three-dimensional Lorentzian function with free acceptance parameters α 0 , α 1 , α 2 , β i , δ i , δ (i, j) and δ (i, j,k) , where indices i, j and k run over (HW, HB, HWB) in case of the acceptance correction for the set of CP-even parameters and over (H W , H B, H WB) in case of the CP-odd parameters. A common parameterisation is used for all production bins since the differences between production bins are shown to be negligible. In addition, the reconstructed event categorisation criteria imposed on the selected Higgs boson candidates and the classification in bins of multivariate NN discriminant values do not impact the acceptance parameterisation.
The impact of reconstruction efficiencies on the parameterisation is also negligible, such that Eq. (4) also holds for the ratio A( c)/A SM of reconstruction-level acceptances defined in Sect. 8. The resulting acceptance parameterisation curves are shown in Fig. 13 for the cases in which all but one of the Wilson coefficients are set to zero. For all cases, the acceptance correction is equal to one at the SM point. In the case of the c H W and c H W B Wilson coefficients, the acceptance corrections reach a maximum value slightly larger than one, leading to the shift of the maximum position from the SM point. This shift is compatible with the statistical accuracy of the fit and the impact of linear EFT terms which are not symmetric around the SM point. The final parameterisation of signal yields relative to the SM prediction in each production bin of the Reduced Stage 1.1 is obtained as the product of the corresponding crosssection, branching ratio and acceptance parameterisations. The expected event yields normalised to the SM prediction are shown in Fig. 14 for each of the CP-even Wilson coefficients after setting all other coefficients to zero. Only production bins with the highest sensitivity to a given Wilson coefficient are shown. The impact of the quadratic terms in the EFT parameterisation can clearly be seen as a non-linear dependence on all but the c H G Wilson coefficient. For comparison, the predictions without the acceptance corrections (σ · B), and without both the acceptance and branching ratio corrections (σ ) are also shown. Both the acceptance and the branching ratio parameterisations have a strong impact on the sensitivity to different Wilson coefficients, especially for the c H W , c H B and c H W B parameterisations in gg2H production bins (Fig. 14a, b, c). Since these coefficients do not enter the ggF production vertex, the corresponding sensitivity is entirely driven by their impact on the decay and the acceptance of selected signal events. The acceptance corrections significantly degrade the sensitivity to the c H W coefficient (see Fig. 14a). Additional sensitivity to this coefficient can be gained from the qq2Hqq production bins as shown in Fig. 14d. The Wilson coefficients c H G and c u H , on the other hand, do not affect the acceptance since they are not present in the decay vertex (Fig. 14e, f). The coefficient c H G still has a non-vanishing impact on the branching ratio through its contributions to the total decay width. Similar effects are also seen for the Wilson coefficients of CP-odd operators.

EFT interpretation results
The ratios of the expected signal yield for a chosen EFT parameter value to its SM prediction are shown in Fig. 15 in each production bin of the Reduced Stage 1.1, together with the corresponding measurement.
The EFT parameterisation of signal yields is implemented in the likelihood function of Eq. (3) using the BSMdependent signal-strength parameters μ p ( c) for each given production bin p, This is then fitted to the observed event yields. Default SM predictions at the highest available order are employed for the cross-sections and branching ratios multiplying the signal strengths in the likelihood function. Modifications of background contributions due to EFT effects are not taken into account.
The fit results with only one Wilson coefficient fitted at a time are summarised in Fig. 16 and in Table 9. The results are in good agreement with the SM predictions. The measurements are dominated by the statistical uncertainty. In the case of the CP-odd coupling parameters, each fit gives two degenerate minima since the corresponding EFT parameterisation contains only quadratic terms which are not sensitive to the sign of the fitted parameter. The fit of the CP-even coupling parameter c u H also results in two minima since the corresponding EFT parameterisation curve in the only sensitive tt H production bin crosses the expected SM crosssection value at two different values of the c u H parameter (see Fig. 14f). The same is true also for the observed tt H cross-section. The small degeneracies for other CP-even coupling parameters are removed by the combination of several sensitive production bins.
The strongest constraint, driven mostly by the ggF reconstructed event categories, is obtained on the c H G coefficient related to the CP-even Higgs boson interactions with gluons. The highest sensitivity to this parameter is reached by the measurements in the gg2H-0 j-p H T -Low and gg2H-0 j-p H T -High production bins due to the highest statistical precision. The sensitivity in the gg2Hp H T -High production bin, which is designed to target the BSM physics effects, is limited due to the small number of events observed in the corresponding reconstructed event category. Additional sensitivity in this bin may be provided by the two-loop interactions which are not implemented in the current simulation of the gg H vertex. The constrained range is stringent enough for the linear approximation to hold, i.e. the quadratic terms in the signal parameterisation are small compared with the linear ones (see Fig. 14e). The constraint on the c H G parameter of the related CP-odd operator is worse by about a factor of three since the linear terms from CP-odd operators do not contribute to the total production cross-section. The constraints on the remaining EFT parameters are weaker, such that both the CP-even and CP-odd signals become dominated by the quadratic terms and are therefore comparable in size. The next-strongest constraints are obtained on the c H B , c H W B , c H W , c H B , c H W B and c H W coefficients that mostly affect the H → Z Z * decays. Due to the larger number of events in the 0-jet reconstructed event categories, the corresponding gg2H production bins provide the highest sensitivity to these decays. Additional smaller sensitivity is obtained from the production vertex of the VBF and VH production modes, with the dominant contribution from qq2Hqq-VBF and qq2Hqq-BSM bins. The latter one is designed to enhance the sensitivity to BSM physics. The qq2Hqq production bins improve in particular the sensitivity to the c H W and c H W  Normalised to SM prediction Normalised to SM prediction  Fig. 17 for several combinations of two CP-even EFT parameters and in Fig. 18 for the corresponding CP-odd operators. The best-fit values as well as the deviation from the SM prediction are shown in Table 10. Good agreement with the SM predictions is observed for all such possible combinations.
The anti-correlation between the c H W and c H B coefficients, as well as between c H W and c H B , is driven by their impact on the signal acceptance. The non-ellipsoidal shape is caused by the acceptance correction, which degrades the original branching ratio-driven sensitivity for increasing parameter values, in particular in the case of the c H W (c H W ) coefficient. The sensitivity is, however, partially recouped by the VBF production vertex.
The 'V'-shaped correlation between the c H G and c H B parameters is due to the interplay between the EFT parameterisation in the ggF production vertex and the parameterisation of the branching ratios and acceptances. The ggF production vertex provides the constraint on the c H G parameter alone, independently of c H B . Due to the decay vertex with its acceptance corrections, this constrained range is shifted upward with increasing values of c H B . Close to the SM point, the constrained c H G range remains approximately the same  The best fit to the data (solid cross) and the SM prediction (star) are also indicated. Except for the two fitted Wilson coefficients, all others are set to zero as without the decay constraints. An additional constraint on c H B is provided by the VBF production mode. Around the SM point, the c H B constraints correspond approximately to those from the one-dimensional parameter fit. Additional sensitivity to intermediate values of the c H B parameter is provided by the acceptance corrections, resulting in two additional allowed parameter regions that are disjoint from the region around the SM point. Similar arguments hold also for the CP-odd case with the c H G and c H B parameters. As opposed to the CP-even case, however, the likelihood contours are symmetric around the c H G = 0 axis, since there are no linear terms contributing to the ggF production crosssection.
The correlation between the c H G and c u H (c H G and c u H ) parameters is introduced through the interference term in the tt H vertex. However, the impact of this term on the final result is negligible since the c H G (c H G ) parameter is already constrained to very small values compared with c u H (c u H ). Therefore, the tt H production vertex mainly constrains the c u H and c u H parameters, while the ggF vertex constrains only the other two. The acceptance correction has no impact on these results. The CP-odd parameter range is less constrained than the CP-even one due to the missing linear c u H terms in the cross-section parameterisation.

Conclusion
Higgs boson properties are studied in the four-lepton decay channel using 139 fb −1 of LHC proton-proton collision data at √ s = 13 TeV collected by the ATLAS experiment. The Higgs boson candidate events are categorised into several topologies, providing sensitivity to different production modes in various regions of phase space. Additional multivariate discriminants are used to further improve the sensitivity in reconstructed event categories with a sufficiently The best fit to the data (solid cross) and the SM prediction (star) are also indicated. Except for the two fitted Wilson coefficients, all others are set to zero large number of events. The cross-section times branching ratio for H → Z Z * decay measured in dedicated production bins are in good agreement with the SM predictions. The inclusive cross-section times branching ratio for H → Z Z * decay in the Higgs boson rapidity range of |y H | < 2.5 is measured to be 1.34 ± 0.12 pb compared with the SM prediction of 1.33 ± 0.08 pb. Results are also interpreted within the κ-framework with coupling-strength modifiers κ V and κ F , showing compatibility with the SM. Based on the product of cross-section, branching ratio and acceptance measured in Reduced Stage-1.1 production bins of simplified template cross-sections, constraints are placed on possible CP-even and CP-odd BSM interactions of the Higgs boson to vector bosons, gluons and top quarks within an effective field theory framework in the H → Z Z * decay. The data are found to be consistent with the SM hypothesis.

Data Availability Statement
This manuscript has no associated data or the data will not be deposited. [Authors' comment: All ATLAS scientific output is published in journals, and preliminary results are made available in Conference Notes. All are openly available, without restriction on use by external parties beyond copyright law and the standard conditions agreed by CERN. Data associated with journal publications are also made available: tables and data from plots (e.g. cross section values, likelihood profiles, selection efficiencies, cross section limits, ...) are stored in appropriate repositories such as HEPDATA (http:// hepdata.cedar.ac.uk/). ATLAS also strives to make additional material related to the paper available that allows a reinterpretation of the data in the context of new theoretical models. For example, an extended encapsulation of the analysis is often provided for measurements in the framework of RIVET (http://rivet.hepforge.org/)". This information is taken from the ATLAS Data Access Policy, which is a public document that can be downloaded from http://opendata.cern.ch/record/413 [opendata.cern.ch].] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 .