Particle-level parton densities and uncertainties for the resolved photon at NLO

Introduction

Resolved, i.e. hadronic, photons are a wonderfully counterintuitive consequence of QCD. Despite being devoid of colour charge, the electromagnetic coupling of photons to quarks enables seeding of a recursive QCD evolution from which hadronic structure emerges in these very unhadronic particles.

The study of resolved-photon structure, which had produced a series of LO and NLO photon-PDF sets [1, 2, 3, 4, 5], essentially paused with the end of HERA operations in 2007, and in the intervening 20 years the focus of the high-energy physics PDF community has instead been on proton PDFs for (and constrained by) LHC data analyses. However, interest is resurgent in anticipation of the Electron–Ion Collider (EIC), scheduled to begin data-taking at Brookhaven in the 2030s [6]. The recent VALO 1.0 fits [7] are an example of this renewed interest: they use historical $F_2^\gamma$ measurements from $e^+e^-$ colliders [8] to produce the first new published photon-PDF sets since 2006 [5]. Fits to $F_2^\gamma$ are the theoretically cleanest and best-established route to the photon’s parton content, since the tagged lepton fixes the probed scale [9, 10].

In this paper we revisit HERA and LEP measurements of jet production in photoproduction and $\gamma\gamma$ collisions. These measurements are numerous and, unlike $F_2^{\gamma\gamma}$, probe the gluon content of the photon directly. As these measurements were typically done at relatively low scales, they are however also subject to large non-perturbative corrections.

Two developments of the past two decades make it worthwhile to revisit these data. First, the much greater precision and technical capability of modern particle-level event generators make it tractable to fit PDFs not to fixed-order partonic calculations but directly to particle-level final states via a shower+hadronisation event generator, so that non-perturbative effects enter the fit itself rather than as external correction factors. Second, PDF error sets have become the standard means of quantifying PDF uncertainties, and the fit presented here provides the first one for the resolved photon from jet data. Using on-the-fly PDF reweighting and a PROFESSOR-style surrogate model [11], the particle-level fit comes at modest computational cost. As the multi-parton-interaction (MPI) model of the generator is untuned for photon-induced collisions, and MPI effects could otherwise be absorbed into the fitted PDFs, we first tune it to the same data; the resulting MPI tunes for HERA and LEP are a further outcome of this study. We describe the technical process and results of these particle-level PDF (PL-PDF) fits, and discuss future prospects for this approach.

The paper is organised as follows: § 2 introduces the particle-level fitting strategy, § 3 the photoproduction and $\gamma\gamma$ processes and the separation of their direct and resolved components, and § 4 the photon-PDF evolution and parametrisation. § 5 describes the event generation, including the MPI model, and § 6 the fitting methodology. We present the results in § 7, discuss prospects in § 8, and conclude in § 9.

Particle-level PDF fits

PDFs are indirectly defined quantities, being not standalone observables but functions related to observables by convolution with a scattering-amplitude function over permitted initial-state configurations. The classic presentation of this is via the factorisation theorem [12],

in which all parton flavours $a,b$ are extracted from initial-state hadrons $A,B$, with all possible momentum fractions $x_a,x_b$ respectively, provided they are compatible with the set of hard-scattering amplitude functions $\mathrm{d}\hat{\sigma}$ mapping to final-state $n$.

In principle, the final state $n$ could be fully exclusive event-topologies as observed at colliders, if there was a full and tractable QFT model through to hadron-level states; in practice, the current limitation is to a partonic final state of very limited multiplicity. The consequence is that fits of PDFs defined under this model, specific to a partonic matrix-element function at fixed loop-order, must be fitted to collider data via theory predictions that are “corrected” from fixed-order parton level to particle level according to some model of radiative corrections and non-perturbative hadronisation, hadron-decay, and multiple-interaction effects. In practice, these correction factors are obtained by use of a particle-level MC generator.

This suggests an alternative picture: rather than restrict to data for which such a non-perturbative correction has been estimated, fit instead to particle-level observables using modern MC-generator tools in which the matrix-element is augmented with radiative corrections via QCD parton-shower algorithms, and generative hadronisation, multiple-interaction and decay models. The loop-order of matrix element supported by this scheme is lower than for the fixed-order approach — currently a maximum of 1- or sometimes 2-loop shower-matched processes as opposed to 2- to 3-loop for traditional PDF fits — but this can be traded against the ability to use a wider range of experimental data, and the fact that a major use of fitted PDFs is as part of the lower-order MC-generation chain. For the resolved photon, the trade-off is particularly favourable: the historical photon-PDF fits, up to and including SAL [5], are at most NLO, which is also the accuracy of the MC@NLO-matched predictions used here, cf. § 5. As differential cross-sections can depend on both the loop- and leg-orders of the scattering process, which typically are respectively lower and higher for ME+PS generators than fixed-order, the current paradigm results in PDF functions that are not obviously optimal for MC-generator predictivity. Dedicated PDF sets for LO event generators, such as MRST LO* [13, 14] and CT09MC [15], address this in part, but are still fitted to fixed-order theory. The PDFs fitted here use a different approach: they are genuine NLO PDFs in the MS scheme with standard DGLAP evolution and can be used in fixed-order calculations. What differs from traditional fits is how non-factorisable corrections are included: instead of ad-hoc external corrections applied to the raw data, we include them on-the-fly using tuned models.

Logistics of particle-level PDF fits

To our knowledge, particle-level PDF fits have not previously been performed, with the practical obstacle being computational tractability. The naïve approach to such a fit would be to run a fit-optimisation algorithm in the space of the PDF parametrisation, with each proposed PDF being used to generate millions of particle-level events and hence compute statistically stable expected observables, before making the next PDF proposal based on their comparison to experimental data. Even with parallel ensembles of optimisers, this approach is intractable due to the multiple hours or days of generation time needed in each serial step. Nor can the interpolated coefficient grids that make fixed-order PDF fits fast [16, 17, 18] be used directly, as the modelling goes beyond fixed order.

However, this limitation has long been worked around for the “tuning” optimisation of generators’ own physics-model parameters, in particular via the PROFESSOR toolkit’s approach [11]. In this framework, an ensemble of efficient algebraic surrogates, or emulators, models each observable-bin value as a function of the model parameters, with established applications at LEP [19] and the LHC [20]. Replacing the generator-parameters of the PROFESSOR methodology with the parameters of a closed-form PDF set requires only minor modification to the established MC-tuning methodology, allowing an MC-coupled PDF set to be “tuned”.

MC tuning typically requires a large up-front (but trivially parallelisable) MC-generation campaign to train the surrogates, as the focus of tuning is on phenomenological components of the models which require explicit re-runs of the generator for each parameter point.¹ For PDF fits, however, a more performant approach is available, as general-purpose event generators support propagation of multiple event-weight streams for PDF variations.² Subject to the usual concerns of statistical support for event-reweighting, the a priori PDF proposals for surrogate-building can be created via single MC-generator runs with a large set of PDF weights, expanded into an ensemble of MC-observable datasets via the RIVET package [23].

Resolved-photon processes

The collinear emission of a quasi-real photon can be computed in the Equivalent Photon Approximation (EPA), for electrons also known as the Weizsäcker–Williams formula [24, 25, 26]. It is based on the observation that photons with small virtualities, $Q^2 < Q_{\max}^2$, can be approximated by real photons, where the virtuality can be integrated over. In terms of the total cross-section this equates to factorising the electron-photon vertex, necessitating a convolution with this flux, additional to the regular PDFs. We use an improved version of the spectrum, following [27], which uses a correction for large $x$, resulting in the following expression:

Here, $x$ denotes the ratio of photon to electron energy, $E_{\gamma}/E_e$, $\alpha_{\mathrm{em}}$ is the electromagnetic coupling constant and $Q_{\max}^2$ and $Q_{\min}^2$ denote the maximal and minimal photon virtuality. The former is effectively a phase-space cut given by the measurement, and the latter can be calculated from kinematic restrictions as

For jet production in photon scattering, we have to take into account both the direct and the resolved part to arrive at a physically meaningful cross-section.

In the direct component, the quasi-real photon couples as a whole, i.e. as a point-like particle, to the hard interaction. In the resolved component, it first fluctuates into a hadronic state, out of which a parton $i$ carrying a fraction $x'$ of the photon momentum enters the hard interaction. The latter is described by the photon PDFs $f_i^{(\gamma)}(x',\mu_F^2)$, which are the subject of this study. The total cross-section is the sum of these components, each of which is written as a convolution of the photon flux(es), the relevant PDFs and the partonic cross-section $\hat{\sigma}$.

Photoproduction at HERA. In electron-proton collisions, a photon emitted by the electron scatters off a parton of the proton, and the cross-section splits into a direct and a resolved component,

with

Here, $x_1$ is the momentum fraction of the photon with respect to the electron, which coincides with the inelasticity $y$ in Eq. (12), $x_1'$ is the momentum fraction of the parton with respect to the photon, and $x_2$ is the momentum fraction of the parton with respect to the proton, described by the proton PDFs $f_j^{(p)}$.

Photon-photon collisions at LEP. In $e^+e^-$ collisions, both leptons emit a quasi-real photon, each of which can either interact directly or be resolved. Correspondingly, the cross-section is the sum of a direct, two single-resolved, and a double-resolved component,

with

where $x_{1,2}$ denote the momentum fractions of the photons with respect to the leptons and $x'_{1,2}$ those of the partons with respect to the photons. While the direct component is independent of the photon PDFs, the single-resolved components are linear and the double-resolved component is quadratic in them.

Beyond leading order, this decomposition is not unique. The direct component contains collinear $\gamma\to q\bar{q}$ splittings, whose singularities are absorbed into the point-like part of the photon PDFs at the factorisation scale, $\mu_F$. Hence, the individual components depend on $\mu_F$ and on the factorisation scheme, and only their sum is a physical observable, in which this dependence cancels up to higher orders. Throughout, we use the MS factorisation scheme, both in the evolution of the photon PDFs and in the NLO calculations.

Separating the resolved and direct components

The HERA and LEP measurements used in our fit distinguish the direct and resolved components through the fraction $x_\gamma$ of the photon momentum entering the hard interaction. At leading order, $x_\gamma= 1$ for direct photons and $x_\gamma< 1$ for resolved ones. As $x_\gamma$ is not directly measurable, the HERA experiments reconstruct the proxy $x_\gamma^{\mathrm{obs}}$ from the transverse energies $E_{T,i}^{\mathrm{jet}}$ and pseudorapidities $\eta_i^{\mathrm{jet}}$ of the $n$ leading jets, the inelasticity $y$, and the electron beam-energy $E_e$, as

In the H1 [28, 29] and ZEUS [30, 31, 32] dijet analyses used here, $x_\gamma^{\mathrm{obs}}$ is built from the two jets with the highest transverse energy, i.e. $n = 2$ in Eq. (12); the three- and four-jet analysis ZEUS_2007_I756660 [33] uses $n = 3$ or 4, depending on the observable. The inclusive-jet and charged-particle analyses do not use $x_\gamma^{\mathrm{obs}}$. The LEP dijet analysis uses the analogous variables

where the numerator sums over the two leading jets and the denominator over the full hadronic final state.

Additional radiation outside the jets and hadronisation shift direct events to values below one, so that $x_\gamma^{\mathrm{obs}}$ equals $x_\gamma$ only at leading order. A cut at $x_\gamma^{\mathrm{obs}} \approx0.75$–$0.8$ conventionally separates a direct-enriched region above from a resolved-enriched region below. As the decomposition into direct and resolved components is itself not unique beyond leading order, these regions are enriched in, but never pure samples of, either component. Unlike $x_\gamma$, $x_\gamma^{\mathrm{obs}}$ is built from jets and is well defined at any order and at particle level. Many measurements are presented in separate $x_\gamma^{\mathrm{obs}}$ regions.

However, $x_\gamma^{\mathrm{obs}}$ is known to receive large non-perturbative corrections [34]. In fixed-order fits, these effects enter as external correction factors, typically derived from leading-order MC generators. In our particle-level approach, they are part of the simulation itself, through the hadronisation and MPI models, cf. § 5.

Photon-PDF evolution

Parton distributions in the case of the well-known proton PDFs satisfy the homogeneous evolution equation

where the $*$ denotes the convolution of the purely hadronic splitting function $P_{ij}$ (including all coupling factors) with the PDF $f_j^{(p)}$ in Bjorken-$x$ space:

The case of the photon is more interesting, as the $\gamma\to q\bar{q}$ splitting introduces an additional source term $k_i$, resulting in the inhomogeneous equation

We may retain the QCD splitting functions of $P_{ij}$, and the photon-parton splitting functions $k_\gamma$ are given in $x$-space in Ref. [35]. At leading order, $k_\gamma^{(0)}$ contains only the quark splitting $k_q^{(0)}$. At higher orders, however, $k_\gamma^{(1)}$ contains the splitting to gluons as well [36]. For our NLO PDFs, the evolution includes the QCD splitting functions $P_{ij}$ up to two loops, of order $\alpha_s^2$, and the photon splitting functions up to $k_\gamma^{(1)}$, of order $\alpha_{\mathrm{em}}\alpha_s$.

We must briefly discuss the structure function of the photon here, which is defined as convolutions of the coefficient functions with the respective quark and gluon PDFs of the photon as well as the additional direct-photon coefficient function $C_\gamma$. In the MS scheme,

which enters $F_2^\gamma$ weighted with $e_q^4$ for each quark flavour and leads to significant negative contributions at large $x$, making the decomposition of $F_2^\gamma$ unstable. In the past, an alternative scheme DIS$_\gamma$ [36] has been proposed to avoid this instability, where this direct contribution is absorbed into the quark PDFs.

This change of factorisation scheme, while leaving $F^\gamma$ unchanged, has consequences for the evolution, as the splitting functions $k_\gamma$ are also redefined:

For our evolution, we use a modified version of the APFEL++ DGLAP evolution library [37] that we call APFEL++$^\prime$, which supports the anomalous terms in photon-PDF evolution and is documented in Ref. [38]. This public version of the code is set up to perform the evolution in the DIS$_\gamma$ scheme. We further modify this version by redefining the non-singlet, singlet and gluon splitting functions of Eqs. (18) and (19) to reflect the MS scheme, as the DIS$_\gamma$ scheme is suited only to inclusive DIS and not to the jet production considered here. The pure QCD splitting functions $P_{ij}$ do not receive any corrections from this change of scheme.

Additionally, the evolution code uses the existing APFEL++ machinery to generate heavier quarks dynamically once the scale reaches their thresholds at or above $Q_0$. By default, this results in the 5-flavour scheme; we modify the code to perform the full evolution in the 4-flavour scheme. The strong coupling runs at three loops with four flavours, with $\alpha_s(m_Z) = 0.112$ (for equivalence with 5-flavour running to the standard $\alpha_s(m_Z) = 0.118$ value).

PDF parametrisation

We use an NLO parametrisation for partons in the photon, proposed by Słomiński, Abramowicz and Levy (SAL) [5]. We use the default starting scale of APFEL++$^\prime$, $Q_0 = 1.3$ GeV, as changing this does not significantly affect the fit. As the evolution starts at the charm threshold, also set to 1.3 GeV, we only parametrise the $u$, $d$, $s$ quarks and the gluon. The charm quark is generated radiatively above its threshold. In the 4-flavour scheme used throughout, there are no bottom quarks in the photon. The parton contents are parametrised at $Q_0$ as follows:

• All parametrised quarks are described as a sum of charge-dependent, point-like contributions and hadronic contributions:

where $A_{\rm PL}$ and $B_{\rm PL}$ are free parameters.

• $f_q^{\rm had}$ contains the hadronic term, parametrised for the $u$, $d$ quarks as

where $A^{\rm had}$ and $B^{\rm had}$ are free parameters. As in SAL, the hadronic term is assumed to depend only on the quark mass, so the difference between the $u,d$ quarks comes from the quark charge in the point-like term.

• The heavier strange quark is suppressed by a fixed factor $r_s = 0.3$, identical to SAL,

• The gluon contains only the hadronic term, and is described by

where $A_G^{\rm had}$ and $B_G^{\rm had}$ float freely in the fit.

Following the approach taken by SAL [5], we fix the large-$x$ exponents to $C^{\rm had}=1$, $C_G^{\rm had}=3$, and the strangeness suppression factor to $r_s=0.3$; we revisit the implications of this choice in § 7.2 and § 8.

Event generation

We derive our predictions using a pre-release version of SHERPA-3.1.0 [39], which provides NLO-accurate hadron-level predictions for photoproduction [40, 41]. One-loop amplitudes are provided by OPENLOOPS 2 [42] for all setups, internally using the COLLIER [43] library. The tree-level matrix elements are provided by COMIX [44] and AMEGIC [45]. To regulate infrared divergences, an automated implementation [46, 47] of the Catani–Seymour subtraction [48] is used and the regulated matrix elements are matched to the CSSHOWER parton shower using MC@NLO matching [49, 50]. The photon flux is provided by the EPA interface that is implemented in SHERPA [40] and is described by the Weizsäcker–Williams formula [24, 25, 26, 27, 51], see Eq. (2). Fragmentation of partons into hadrons is modelled with the AHADIC hadronisation model [52], tuned to LEP data [53].

Throughout, we work in a 4-flavour scheme, with the bottom mass set to $m_b = 4.92\,\mathrm{GeV}$. We set the renormalisation and factorisation scales to half the scalar sum of the transverse momenta of the final-state partons, $\mu_F = \mu_R = H_T/2$, and the resummation scale, which sets the starting scale of the parton shower, to its default $\mu_Q = \mu_F$.

For the HERA event generation, we use the PDF4LHC21 set [54] for the proton and take the strong coupling from it. For LEP, we use the same value, $\alpha_s(m_Z) = 0.118$, with three-loop running. This is the five-flavour value corresponding to the four-flavour coupling used in the evolution of the photon PDF, cf. § 4. All other settings, for instance those of the beam remnants and their primordial transverse momentum, are kept at their tuned defaults.

HERA setup

We generate photoproduction events at the HERA beam energies, $E_e = 27.6\,\mathrm{GeV}$ and $E_p = 920\,\mathrm{GeV}$ or $820\,\mathrm{GeV}$, depending on the analysis. For most analyses, the direct and resolved processes are

matched to the CSSHOWER dipole shower with MC@NLO [50]. Here, $j$ is any light parton in the event generation: $u,d,s,c,g$. For the three-jet observables of ZEUS_2007_I756660 [33], the resolved and direct processes are generated with one more jet at NLO accuracy,

and for the four-jet observables of the same analysis, four jets are generated at NLO,

LEP setup

We generate events at the LEP centre-of-mass energies of the analyses, $\sqrt{s} = 198\,\mathrm{GeV}$ for OPAL_2003_I611415 and L3_2004_I661114 and $206\,\mathrm{GeV}$ for OPAL_2008_I754316. For all analyses, we generate direct, single-resolved, and double-resolved events

matched to the CSSHOWER dipole shower with MC@NLO [50]. For the single-resolved process, we only generate the events in one beam orientation, and double the cross-sections. This is sufficient because all LEP observables used here are symmetric under the exchange of the two beams.

Multi-parton interactions

SHERPA simulates multi-parton interactions (MPI) with AMISIC [39], which is based on the model of Sjöstrand and van Zijl [55]. MPI are only simulated if both incoming particles are resolved, i.e. described by PDFs, in which case the resolved photon is treated like a hadron. At HERA, this is the case for the resolved component and at LEP for the double-resolved one, whereas direct and single-resolved events receive no MPI.

The secondary scatters are leading-order QCD $2 \to2$ processes, computed with the PDFs of the incoming photons and protons. They are generated in decreasing order of transverse momentum $p_{\mathrm{T}}$, starting from the scale of the hard process. Their mean number per event is set by the ratio of the integrated $2 \to2$ cross-section above $p_{\mathrm{T},\min}$ to the non-diffractive cross-section $\sigma_{\mathrm{ND}}$. The latter is taken from the parametrisations of Schuler and Sjöstrand [56, 57] and rescaled by a factor $\sigma_{\mathrm{norm}}^{\mathrm{ND}}$. For the photon, it is a sum over $\rho$, $\omega$, $\phi$, and $J/\psi$ states in the vector-meson-dominance picture, weighted with their couplings to the photon. The scatter probability also depends on the impact parameter $b$, through the overlap of the Gaussian matter distributions of the two incoming particles, whose radii are parameters of the model. At each centre-of-mass energy, both radii are multiplied by a common factor such that the $b$-integrated probability for at least one scatter reproduces $\sigma_{\mathrm{ND}}$; where the integrated $2 \to2$ cross-section above $p_{\mathrm{T},\min}$ is smaller than $\sigma_{\mathrm{ND}}$, the radii are used unchanged. In events with a hard process, $b$ is drawn from this overlap, favouring central collisions with more MPI activity. The model does not separate the vector-meson-like and the anomalous, point-like components of the resolved photon: the full photon PDF enters the secondary scatters, and the non-diffractive cross-section and the matter distribution are those of the vector-meson picture.

The divergence of the $2 \to2$ matrix elements at small $p_{\mathrm{T}}$ is regularised by a factor $\left[p_{\mathrm{T}}^2/(p_{\mathrm{T}}^2+p_{\mathrm{T},0}^2)\right]^2$. The strong coupling is evaluated at a scale proportional to $p_{\mathrm{T}}^2+p_{\mathrm{T},0}^2$ and no scatters are generated below the cut-off $p_{\mathrm{T},\min}$. Both scales are parametrised in terms of the centre-of-mass energy $E_{\mathrm{CMS}}$ of the photon–proton or photon–photon system as

with a common energy-scaling exponent $\eta$, a reference energy $E_{\mathrm{CMS}}^{\mathrm{ref}}$, and lower bounds $p_{\mathrm{T},0}^{\mathrm{IR}}$ and $p_{\mathrm{T},\min}^{\mathrm{IR}}$, by default 0.5 GeV and 1 GeV. We set $E_{\mathrm{CMS}}^{\mathrm{ref}}$ to the beam centre-of-mass energy of the respective collider. Since the photon energies follow the EPA spectrum, $E_{\mathrm{CMS}}$ varies from event to event.

Fitting methodology

The fit follows the surrogate-based strategy outlined in § 2; here we describe its implementation.

Surrogate optimisation, or “tuning”

In the surrogate approach, every bin of every histogram is given an independent surrogate function of the continuous parameters being optimised, typically an interpolation between randomly sampled points within the parameter-space hypercube or hypersphere.

The PROFESSOR package uses simple multivariate polynomials, within possibly transformed parameter measures; more advanced emulator techniques such as neural networks or Gaussian processes equally meet the criteria. Armed with an ensemble of surrogate functions, applications such as serial optimisation become tractable, even in real-time for visualisation tasks. In this study, we use PROFESSOR v2.5.9 with third-order polynomials for the PDF fit. The PDF fit has six free parameters, and the MPI tunes have six for HERA and five for LEP, where there is no proton radius.

The prof2-sample program was used to sample parameter points $\{\mathbf{p}\}$ uniformly within the hypercube of PDF parameters documented in § 4.1, with the ranges listed in Table 5, injecting those parameters into steering scripts for a custom PDF grid-writer program linked against the APFEL++$^\prime$ evolution library described in § 4. The grid-writer evaluated the functional forms from § 4.1 to set the PDF distributions at the input scale $Q_0 = 1.3\,\mathrm{GeV}$, and the APFEL++$^\prime$ evolution was used to generate a PDF grid in $(x,Q^2)$ in the LHAPDF 6 data format [58], from $x = 10^{-5}$ to 1, and in $Q^2$ up to $10^4\,\mathrm{GeV}^2$. Parameter points whose evolved PDFs became negative were rejected, so that they entered neither the MC reweighting nor, consequently, the interpolation. The remaining PDF variations were combined into a single LHAPDF set, which was provided to the SHERPA MC generator as if it were a standard, established global PDF with fit-uncertainty variations to be propagated via event-weighting.

Parameter	Central	Mean	Std. dev.	Sampling range
$A_{\mathrm{PL}}$	7.991	5.967	1.505	[0, 8]
$B_{\mathrm{PL}}$	2.378	1.700	0.464	[0, 3.5]
$A^{\mathrm{had}}$	0.5985	0.5383	0.1760	[0, 0.6]
$B^{\mathrm{had}}$	1.4670	1.1274	0.4047	[−2.5, 1.5]
$A^{\mathrm{had}}_{G}$	0.0001	0.6263	0.6523	[0, 2]
$B^{\mathrm{had}}_{G}$	−0.8604	0.4556	0.7518	[−3, 1]

A single set of runs to cover the involved process types was then sufficient to propagate the randomly sampled PDF parametrisations through SHERPA, using its existing on-the-fly reweighting implementation at NLO [59] to predict the differential cross-sections for all input PDFs at once, via the set of photoproduction event-analyses implemented in the RIVET analysis framework. The Python interface to the YODA histogramming package [60] used by RIVET was used to split the resulting nominal and variation observable-histograms into per-variation copies associated with the original sampled parameter variations, and used as input to the standard PROFESSOR interpolation-building and tuning.

The PROFESSOR tuning finds an optimised parameter vector $\mathbf{p}_{\mathrm{tune}}$ that minimises the heuristic goodness-of-fit function

where $b$ indexes the set of summed-over observable bins, $y_b(\mathbf{p})$ is the PROFESSOR surrogate prediction for bin $b$ with parameters $\mathbf{p}$, $d_b$ is the data value of bin $b$, and the denominator contains a quadrature sum of uncertainties for bin $b$. Of these, $\sigma(y_b(\mathbf{p}))$ is a representative uncertainty from the surrogates—in this case the median statistical uncertainty of the MC runs—and $\sigma(d_b)$ and $\epsilon_b d_b$ are respectively the reference-data uncertainty and an ad hoc uncertainty for MC theory; we use a value of 5% for the relative uncertainty $\epsilon_b$ on all bins. This form of the fit-function treats all bins as independent, which is unavoidable as the available data lacks the uncertainty breakdowns and correlations required for a more precise maximum-likelihood fit.

To avoid local minima, we additionally scanned 10 000 points of the parameter space with the PROFESSOR $n$-ball sampler,³ and used the point with the lowest $\chi^2$ as the starting point of the MINUIT minimisation.

As, due to model and correlation-data shortcomings, the fit function is not in practice distributed like a classical $\chi^2$ test-statistic, the replica method is used to estimate fit uncertainties resulting from the observable-measurement uncertainties. To produce the replicas, we start from the reference data in the YODA files and create replica datasets by shifting the data value $d_b$ of each bin $b$ in replica $r$ as

where $\sigma_b$ is the symmetric Gaussian uncertainty of bin $b$ and $\mathcal{N}^{(r)}(0,1)$ is a standard-normal random number. Asymmetric uncertainties are symmetrised for the smearing, but stay asymmetric in the replica, and the uncertainties are not smeared, $\sigma_b^{(r)} = \sigma_b$. Within a replica, the random numbers are generated from a single fixed seed.

The PDF fit was repeated for each of the 100 replica datasets, resulting in 100 replica PDFs. The central PDF is the fit to the unsmeared data, and the replica PDFs can deviate from it in either direction, so that the resulting uncertainties are asymmetric. The uncertainty bands in § 7 show the 68 % interval of the replica PDFs. The replica uncertainties therefore reflect the experimental uncertainties of the individual bins. They do not account for correlations between systematic uncertainties, or for theory variations such as the choice of scales.

MPI tuning

A complicating factor in this analysis is that photoproduction has significant phase-space overlap with the event-generator mechanism of multi-parton interactions (MPI), which model secondary partonic scatterings within hadron-hadron collisions to describe sensitive $ep$ and hadron-collider observables. MPI models are highly phenomenological rather than first-principles, and the setting of MPI parameters has long been a core focus of MC tuning.

As the resolved photon is an effective hadron, and one expected to have quite different properties from the (anti)proton collisions on which the SHERPA default MPI tune was established, poor MPI modelling of resolved $\gamma p$ and $\gamma\gamma$ collisions risks incorrect conclusions about resolved-photon structure: effects that should be described by the MPI model would be “absorbed” into the PDFs. Moreover, as SHERPA-3.1 has an untuned MPI model for HERA and LEP photoproduction, setting a good baseline description for these scenarios with existing PDFs is useful in its own right.

We therefore performed the tuning in two steps. First, we tune the MPI parameters listed in Table 2 with the photon PDF fixed to the historical SaS 1M set [2, 61], the SHERPA default. We then fix the MPI parameters at their tuned values and perform a joint PDF fit to the HERA and LEP data. Both steps use largely the same measurements: the MPI tunes include additional distributions and extended ranges at lower transverse energies, while the PDF fit additionally includes ZEUS_1997_I450085 and the shapes of the four-jet distributions of ZEUS_2007_I756660, given in detail in Appendices A and B.

Parameter	Range	HERA	LEP
$p_{\mathrm{T},0}^{\mathrm{ref}}$ [GeV]	[0.25, 10]	1.77	5.9
$p_{\mathrm{T},\min}^{\mathrm{ref}}$ [GeV]	[0.8, 3]	1.683	1.56
$\eta$	[0, 0.3]	0.2637	0.055
$\sigma_{\mathrm{norm}}^{\mathrm{ND}}$	[0.1, 1.2]	1.069	0.103
$R_{\gamma}$ [fm]	[0.6, 0.9]	0.6743	0.6172
$R_{\mathrm{p}}$ [fm]	[0.6, 1.1]	1.0496	—
$\chi^2/N_{\mathrm{df}}$		1698.2/508 = 3.34	498.2/129 = 3.86

Separate MPI tunes were made for $\gamma p$ collisions at HERA and $\gamma\gamma$ collisions at LEP, both with SaS 1M, and each tune is used for the predictions of its collider in the PDF fit. Due to the different beam energies and particles at the two colliders, these need to be tuned separately. It may be possible to do a combined tune, but this is an active field of research (see Ref. [62]) and for the purposes of this exercise, keeping the complexity to a minimum with simple individual tunes was preferred. With HERA, the additional proton-radius parameter must also be considered; here, we use single Gaussians to model the spatial form-factors of each particle.

For HERA, 200 sets of parameters were used for the MPI variations, and with LEP, reflecting the one less free parameter, we used 150. For both, we performed the interpolation with second-order polynomials.

Observable selection

The event analysis was performed using RIVET 4 and YODA 2, with the analyses listed in tbl:1. We exclude the distributions and bins restricted to the direct-enriched region, as including them changes neither the MPI nor PDF fit. We also exclude the four-jet distributions

caption: Data used in both PDF fitting and MPI tuning. The selected histograms and bins are listed in Appendices A and B.

:::

of ZEUS_2007_I756660 from the MPI tunes, and use them in the PDF fit only as shapes, normalised to unit area, as the NLO predictions lie consistently above the data and would otherwise pull the fit. The $x_{\gamma}^{\mathrm{obs}}$ cuts keep bins that straddle the cut, e.g. $0.7 < x_{\gamma}^{\mathrm{obs}} < 0.85$ in H1_2002_I581409, so the resolved-enriched selections extend into the direct-enriched region. We selected the data by hand, choosing the histograms and bins that are most sensitive to the tuned parameters, mostly from envelope plots, which compare the spread of the predictions for all sampled parameter points with the data. Automated filters, such as the envelope filter of PROFESSOR, serve only as an additional safeguard.

PDF sampling

The six PDF parameters were sampled as described in § 6.1, within the ranges listed in Table 5, using 3000 PDFs at once for all PDF tuning MC generator runs.

Fit results

MPI tunes

Table 2 lists the tuned MPI parameters for HERA and LEP, and Table 3 documents the fit quality of each analysis. The tuned predictions give $\chi^2/N_{\mathrm{df}} = 3.34$ for HERA and 3.86 for LEP. At HERA, all parameters lie inside their ranges, with $\eta$ and $R_{\mathrm{pclose}}$ close to the upper edges. The

Analysis	Bins	$\chi^2$	$\chi^2$/bin
HERA
H1_2002_I581409	67	73.2	1.09
H1_2006_I711847	43	37.7	0.88
ZEUS_2001_I568665	56	216.1	3.86
ZEUS_2007_I753991	62	240.7	3.88
ZEUS_2007_I756660	192	341.4	1.78
ZEUS_2012_I1116258	94	788.9	8.40
LEP
L3_2004_I661114	5	17.8	3.56
OPAL_2003_I611415	113	411.7	3.64
OPAL_2008_I754316	11	68.6	6.24
Total HERA	514	1698.2	3.34
Total LEP	129	498.2	3.86

reference scales $p_{\mathrm{T,0}}^{\mathrm{ref}}$ and $p_{\mathrm{T,min}}^{\mathrm{ref}}$ apply at $E_{\mathrm{CMS}}^{\mathrm{ref}}$, the centre-of-mass energy of the collider, and decrease towards lower photon–proton energies, cf. Eqs. (28) and (29): at $E_{\mathrm{CMS}} = 200\,\mathrm{GeV}$, the tuned values correspond to $p_{\mathrm{T,0}} \approx1.4\,\mathrm{GeV}$ and $p_{\mathrm{T,min}} \approx1.3\,\mathrm{GeV}$.

At LEP, MPI enter only the double-resolved component and the tune prefers a larger regulator than at HERA, $p_{\mathrm{T,0}} \approx4.8\,\mathrm{GeV}$ at $E_{\mathrm{CMS}} = 30\,\mathrm{GeV}$, while $p_{\mathrm{T,min}}$ is similar at $1.3\,\mathrm{GeV}$, and $\sigma_{\mathrm{norm}}^{\mathrm{ND}}$ lies at the lower edge of its range. The photon radius, $R_{\gamma} = 0.62\,\mathrm{fm}$, close to the lower edge of its range, is similar to that of the HERA tune, $0.67\,\mathrm{fm}$. For both colliders, the fits prefer the photon to be smaller than the HERA proton, as expected.

Per analysis, the H1 dijet data are described well, with $\chi^2/N_{\mathrm{bins}} \approx1$. At HERA, the ZEUS inclusive-jet and three-jet data, ZEUS_2012_I1116258 and ZEUS_2007_I756660, account for two thirds of the $\chi^2$, and at LEP, OPAL_2003_I611415 contributes 411.7 of 498.2. At LEP, OPAL_2008_I754316 has the largest $\chi^2$ per bin, 6.24, from 11 bins. ZEUS_2012_I1116258 and OPAL_2008_I754316 are poorly modelled across the full range of distributions. In ZEUS_2012_I1116258, low $E_{\mathrm{T}}$ and most forward $\eta$ bins affect the $\chi^2$ the most, while in OPAL_2008_I754316, intermediate $p_{\mathrm{T}}$ bins between 7.5 GeV and 15 GeV in both $\eta$ regions are the worst described.

Figure 1 compares the tuned predictions with data for a selection of the distributions in the tune, together with predictions without MPI and with the default MPI parameters of SHERPA. At LEP, the tuned predictions are close to those without MPI, consistent with $\sigma_{\mathrm{norm}}^{\mathrm{ND}}$ at the lower edge of its range, whereas the default parameters overshoot the OPAL data at low $x_{\gamma}^{\mathrm{obs}}$ by more than a factor of two. Other than for the ZEUS inclusive-jet photoproduction data [63], the tuning is seen to significantly improve the data description, with the tunes suppressing MPI activity to a level closer to the MPI-off mode.

Vector (PDF)

Vector (PDF)

Vector (PDF)

Vector (PDF)

Vector (PDF)

Vector (PDF)

PDF tune

The PDF tune was performed using the parameter ranges given in Table 5, with the MPI modelled using the SaS 1M PDF and tuned parameters, and the sampled PDFs for the resolved scattering process. The MC statistics were doubled relative to the MPI tunes for analyses where the median $\sigma_{\mathrm{MC}}/\sigma_{\mathrm{data}}$ was greater than 0.50 (detailed in Appendix § 11, Table 10).

Analysis	Bins	Median $\sigma_{\mathrm{MC}}/\sigma_{\mathrm{data}}$	Median $\sigma_{\mathrm{MC}}/\lVert\mathrm{MC}\rVert$
HERA
H1_2002_I581409	67	0.14	3.0%
H1_2006_I711847	43	0.32	4.3%
ZEUS_2001_I568665	56	0.90	8.8%
ZEUS_2007_I753991	62	0.40	4.9%
ZEUS_2007_I756660	192	0.85	11.9%
ZEUS_2012_I1116258	94	0.51	3.9%
LEP
L3_2004_I661114	5	0.22	2.1%
OPAL_2003_I611415	113	0.51	6.1%
OPAL_2008_I754316	11	0.15	1.7%

The final fit metrics are given in Table 4, with significant improvements seen in the ZEUS dijet and inclusive-jet data compared with the MPI tunes with SaS 1M (Table 3): $\chi^2/N_{\mathrm{bins}}$ falls from 3.9 to 0.8 and 1.3 for the dijet data, ZEUS_2001_I568665 and ZEUS_2007_I753991, and from 8.4 to 1.9 for the inclusive-jet data, ZEUS_2012_I1116258. In contrast, the ZEUS three-jet data, ZEUS_2007_I756660, and the L3 inclusive-jet and OPAL dijet data, L3_2004_I661114 and OPAL_2003_I611415, are described less well than in the MPI tunes, and the three-jet and OPAL dijet data together account for three quarters of the total $\chi^2$. As the bin selections of the PDF fit and the MPI tunes differ, these comparisons are indicative only.

Analysis	Bins	$\chi^2$	$\chi^2$/bin
HERA
H1_2002_I581409	52	29.2	0.56
H1_2006_I711847	46	68.3	1.48
ZEUS_1997_I450085	29	17.1	0.59
ZEUS_2001_I568665	55	42.1	0.77
ZEUS_2007_I753991	58	75.6	1.30
ZEUS_2007_I756660 (3-jet)	190	946.5	4.98
ZEUS_2007_I756660 (4-jet)	183	3639.6	19.89
ZEUS_2012_I1116258	96	183.7	1.91
LEP
L3_2004_I661114	9	54.8	6.09
OPAL_2003_I611415	112	507.3	4.53
OPAL_2008_I754316	12	32.9	2.74
Total	659	1957.6	2.97

Figure 2 compares the data for a selection of distributions with the predictions for PL$\gamma$26 and two older sets, GS96 (GSG961) [67] and ACFGP (ACFGPG) [68], all with the tuned MPI parameters. At HERA, PL$\gamma$26 describes the resolved-enriched ZEUS dijet data best: its predictions lie within 5 % of the data for $0.2 < x_{\gamma}^{\mathrm{obs}} < 0.6$, and within 8 % for $x_{\gamma}^{\mathrm{obs}} < 0.75$ with a central leading jet, whereas those with GS96 and ACFGP deviate by up to 20 %. For the ZEUS inclusive-jet data, PL$\gamma$26 and ACFGP agree to within a few per cent, both about 10 % below the data. At LEP, all three sets overshoot the OPAL dijet data with one $x_{\gamma}^{\pm} < 0.75$, PL$\gamma$26 by 20–30 %, ACFGP by 30–40 % and GS96 by 50–60 %. All three also overshoot the L3 inclusive-jet spectrum below $p_{\mathrm{T}} = 10$ GeV, PL$\gamma$26 by up to 40 %, and undershoot it above 15 GeV, independently of the PDF, by a factor of four at the highest $p_{\mathrm{T}}$.

Vector (PDF)

Vector (PDF)

Vector (PDF)

Vector (PDF)

Figure 3 shows the parton distributions of the fitted set, PL$\gamma$26, at $Q^2 = 2$, 5, 20 and 50 GeV$^2$. The up quark is the largest quark distribution, as expected from its larger electric charge, $e_{\mathrm{u}} = \frac{2}{3}$, while the down and strange quarks, with $e_{\mathrm{d,s}} = -\frac{1}{3}$, have similar distributions, except at intermediate $x$, where the down quark is larger as its hadronic component peaks there. The charm quark, generated radiatively above its threshold, grows with $Q^2$, in particular at large $x$. At $Q^2 \ge5$ GeV$^2$, the gluon dominates at small $x$ and falls steeply towards large $x$.

Vector (PDF)

In Figure 4, we compare PL$\gamma$26 with other photon-PDF sets at $Q^2 = 5$ GeV$^2$: the NLO sets SAL [5] in the DIS$_\gamma$ scheme, and ACFGP [68] and VALO 1.0 [7] in the $\overline{\mathrm{MS}}$ scheme. GS96, shown in Figure 2, is not defined at this scale.

Vector (PDF)

Towards $x = 1$, our sets diverge from the SAL PDFs. This is partly an effect of the factorisation scheme: in the DIS$_\gamma$ scheme, the quark distributions fall towards zero at large $x$, whereas in the $\overline{\mathrm{MS}}$ scheme they continue to rise before falling steeply to zero. Part of the instability in this region also comes from the interpolation of the LHAPDF grids, which we mitigate with a high density of grid points at large $x$.

At intermediate $x$, our up- and strange-quark PDFs lie 40–65 % above those of ACFGP and 60–120 % above those of VALO 1.0, and our down quark about a factor two above ACFGP, while it agrees with VALO 1.0 for $x \approx0.6$–0.8 and lies below it towards low $x$. Our gluon is smaller than in both $\overline{\mathrm{MS}}$ sets, by a factor of 1.6–2 relative to VALO 1.0 and about three relative to ACFGP. At low $x$, our gluon diverges significantly from the other sets, as it rises much less steeply towards small $x$ (Figure 4). We note that our quark uncertainties are also greater in the low-$x$ regime, growing from a few per cent at $x = 0.9$ to 25–30 % at $x = 0.1$, simply because we do not have the data constraining this. The low-$x$ gluon band, of about 15%, is narrower, but reflects the restricted gluon parametrisation discussed below rather than a constraint from data. To get a handle on this area, we must look towards the $F_2^\gamma$ measurements made at LEP. The bands are asymmetric, with the central quark PDFs at their upper edges and, above $x \approx0.3$, the central gluon at its lower edge, consistent with $A_{\mathrm{PL}}$, $A^{\mathrm{had}}$ and $B^{\mathrm{had}}$ of the central fit lying at the upper edges of their sampling ranges and $A^{\mathrm{had}}_{\mathrm{G}}$ at the lower edge (Table 5). We note also

that the charm PDF has small (invisible) error bands. As the charm is not parametrised at $Q_{0}$, the central and replica bands rely on the DGLAP evolution itself. Most of the high-$x$ charm is a result of the photon splitting $\gamma\to c\bar{c}$ and to a lesser extent the gluon splitting $g \to c\bar{c}$. Thus, the replica variations in the charm that we would naturally expect to follow from variations in the gluon and light quarks are not large enough to be visible on the plots, but this obviously does not reflect a well-constrained charm quark.

The exponent $B_{G}^{\mathrm{had}}$, which is primarily responsible for the low-$x$ tail, can, at very large negative values, also affect the high-$x$ gluon region that we are sensitive to, but such values make the momentum integral diverge, so we restricted the fit to $B_{G}^{\mathrm{had}}, B_{G}^{\mathrm{had}} > -2.0$ (Table 5). With this restriction, the central fit prefers $A_{G}^{\mathrm{had}}$, the only other free gluon parameter, close to zero, so that the gluon is generated almost entirely by the evolution and $B_{G}^{\mathrm{had}}$ is not constrained. We are content with this because it allows our high-$x$ modelling to describe jet physics at HERA and LEP better than previous sets, and demonstrates that our machinery does work, even with such parametrisation constraints. In future PL$\gamma$26 work, however, this does motivate a more complex parametrisation of the gluon (§ 8).

Differences between VALO 1.0 and our PL$\gamma$26 PDFs may be due to differences in the physics impact of the observable selection. While $F_{2}^{\gamma}$, used in VALO 1.0, fixes the quark content of the photon quite precisely, it has practically no handle on its gluon content, which is therefore largely set by the parametrisation. Conversely, using jet production only, as in the PL$\gamma$26 set, fixes the combined effect of gluons and quarks but it does not provide significant discrimination to separate their relative impact. As a consequence, the relatively larger gluon of VALO 1.0, compared to PL$\gamma$26, is compensated by a larger quark PDF in the latter. This motivates further studies taking into account both jet production and its more differential observables together with the more inclusive $F_{2}^{\gamma}$, as discussed in § 8.

Prospects for further development

As a new method, there are many possible avenues for further investigation and refinement. Firstly, as noted in § 4.1, in this study we have used the parametrisation of SAL [5], in which the large-$x$ exponents are fixed to $C_{\mathrm{had}}=1$ and $C_{G}^{\mathrm{had}}=3$. SAL took these values from counting rules as their data did not reach $x \simeq1$, but the jet data used here have a degree of sensitivity to the high-$x$ region. An extended PDF tune could potentially achieve a better fit by permitting the high-$x$ scaling exponents to vary, freeing other elements of the parametrisation to achieve a better description of lower-$x$ observables. We will pursue this in a follow-up study, in which freeing $C_{G}^{\mathrm{had}}$ will also test whether the differences between our gluon and those of fits to $F_{2}^{\gamma}$ arise from the gluon parametrisations, as discussed in § 7.2. Similarly, the fixed strange-suppression factor $r_{\mathrm{s}}=0.3$ would ideally also be determined from data. The familiar issue of parametrisation bias also raises the possibility of a more parametrisation-free basis, e.g. using the NNPDF approach of relatively unconstrained non-parametric forms⁴ for initial sampling [69], then rendered into a parametrised form for surrogate optimisation by MC-to-Hessian conversion [70].

A further step is a combined fit to the jet data and to the $F_{2}^{\gamma}$ measurements [8] used in fits such as VALO 1.0 [7]. These constrain the quark densities directly, in particular at low $x$, while the jet data fix the combined effect of quarks and gluons. Such a fit would therefore separate the two, and show whether the compensation between quarks and gluon discussed in § 7.2 fully accounts for the differences between PL$\gamma$26 and fits to $F_{2}^{\gamma}$, or whether there is a tension between the two types of data.

Away from the PDF parametrisation, our approach taken to factorise the MPI tuning from the PDF tune was for tractability: fit stability and number of required sample points scale strongly with the number of parameters. As the MPI tunes were performed using the SaS PDF, there is an inconsistency between it and the final use-case. In practice, MC generators are able to use separate PDFs for hard-process and MPI modelling, avoiding the need to create an MPI tune for every prospective resolved-photon PDF, but a single PDF suitable for all aspects of MC generation would be both convenient and a test of model completeness.

The MPI-tune factorisation between HERA and LEP could also be refined, and combined with $pp$ and $p\bar{p}$ underlying-event and minimum-bias data, to attempt a universal description of MPI able to describe observables at $e^{+}e^{-}$, $ep$, and hadron colliders with a single tune. While very interesting, this is a speculative goal and in the current work the pragmatic approach was taken to find less ambitious MPI configurations able to minimise biasing of the resolved-photon PDF extraction.

Another prospect for future refinements in fits of this kind is to incorporate estimations of MC-generator physics uncertainties into the fit. At present, theory uncertainties are aggregated into the blanket $\epsilon=5\%$ regularisation term in Eq. (30). But modern event-generator codes can provide far more detailed breakdowns of sensitivity to modelling systematics such as perturbative scale choices — and increasingly also non-perturbative modelling such as hadronisation and MPI [21, 22] — via event-weighting, a facility not yet used in this fit. Inclusion of these systematic uncertainties would provide a further guard against overfitting to a particular parton shower, or to less-constrained elements of the non-PDF MC modelling, as well as in principle permitting pairing of resolved-photon PDF uncertainties with a particular proton-PDF error set for $ep$ physics. This extension, along with including an informed estimate of surrogate accuracy in the $\chi^{2}$ definition, will require extensions in the Professor toolkit, but is aligned with planned developments for MC tuning.

Conclusions

In this work we have presented a new PDF fit for the resolved photon, $PL\gamma26$, using measurements of jet production in photoproduction at HERA and $\gamma\gamma$ scattering at LEP.

To arrive at the fit we take a novel approach. As these data typically lie at rather low scales, non-factorisable corrections like hadronisation or multi-parton interactions are not negligible. Hence, instead of the usual approach of relying on post-processing of the experimental data to “subtract” event-generator models of these effects and produce data comparable with fixed-order partonic calculations, we perform the PDF fit directly to the particle-level data. This ensures forward application of the non-factorisable modelling specific to each PDF configuration, and is made possible by using the general-purpose event generator SHERPA at NLO+PS accuracy.

As a first step, we tuned the SHERPA multi-parton interactions model separately to HERA and LEP jet data available in the RIVET toolkit (including new analyses added for this purpose), to avoid absorbing such effects into the PDF. The resulting HERA and LEP MPI tunes are themselves of interest for ongoing studies to quantify such effects and optimise generator configurations for photoproduction and DIS measurements, e.g. in ultraperipheral heavy-ion collisions at the LHC [71, 72, 73], at the upcoming Electron–Ion Collider [6, 41, 74], and in studies for FCC-ee [10, 75, 76].

The resolved-photon PDF tune was then performed with the tuned SHERPA generator, on a sensitive subset of HERA and LEP data. Using the 6-parameter PDF ansatz from Ref. [5] at evolution starting-scale $Q_0 = 1.3$ GeV, points within the 6-dimensional parameter space were sampled and evolved using APFEL++′ [38], an appropriately adapted version of the APFEL++ library. This created a set of 3000 random resolved-photon PDFs in the LHAPDF data format; treated as a single PDF set and using the NLO+PS reweighting mechanism in SHERPA, particle-level observable predictions were then efficiently computed for all sampled PDF configurations within a single MC-generator run. The PROFESSOR MC-tuning framework was then used first to construct surrogate models for each observable bin in the 6-parameter PDF space, and secondly to optimise the fit to data with PROFESSOR’s standard least-squares metric. A replica-fit approach with 100 replicas from Gaussian-smeared data was used to estimate the uncertainties. The resulting functions are the first MS resolved-photon PDFs fitted to jet data from both HERA and LEP, and provide an error set for quantitative uncertainty estimation. Compared with fits to $F_2^\gamma$, $PL\gamma26$ has larger up- and strange-quark densities and a smaller gluon. As $F_2^\gamma$ hardly constrains the gluon, while the jet data fix the combined effect of quarks and gluons, we attribute this mainly to the gluon parametrisations of those fits, compensated by larger quarks in $PL\gamma26$; a combined fit to both types of data and a more flexible parametrisation will test this reading.

This generator-based approach to PDF fitting is particularly valuable in settings where data are sparse and the energy scales low enough to be sensitive to non-factorisable effects. While the order of the fit is restricted to NLO, this is competitive for the photon where NLO is currently the standard, as well as by construction resolving any debate about the formal compatibility of PDFs with MC event-generation. These fits provide an important stepping stone to updated phenomenological studies of jet production in ultraperipheral collisions at the LHC, and in photoproduction at the EIC.

This fit uses complementary data compared to the fit performed in Ref. [7], which was performed to $F_2^\gamma$ measurements. To our knowledge, it is the first photon-PDF fit to a wide range of jet measurements from different experiments, combining ten photoproduction and $\gamma\gamma$ analyses from H1 and ZEUS at HERA and from L3 and OPAL at LEP. Earlier fits were constrained mainly by $F_2^\gamma$, adding at most a few jet measurements from a single collider [5, 67]. The resulting set, $PL\gamma26$, and its replicas are provided in the LHAPDF format as PLgamma26, with member 0 the central fit to the unsmeared data and members 1–100 the replicas. Until the MPI are tuned with $PL\gamma26$ itself (§ 8), we recommend using SaS 1M with the tunes of Table 2 for the MPI, as in our fit. This work, while not yet comparable to the scale and depth of global PDF fits and with a number of immediate avenues identified for further investigation, serves both to demonstrate the practical viability of particle-level PDFs via the technical advances in modern MC event generators and to highlight that a combined PDF- and generator-tuning approach may be the most robust and practically useful way to define PDFs for future low-scale physics at EIC and other facilities.

Acknowledgements

Our thanks to Valerio Bertone for advice and assistance with the APFEL++ code, and obtaining a photon-PDF evolution with current tools. We also thank Claire Gwenlan for advice and assistance with the ZEUS 3- and 4-jet analysis, and Ilkka Helenius for providing PYTHIA configurations and the preliminary MS VALO PDF.

Funding acknowledgements for each author follow:

• AB acknowledges funding via the STFC experimental Consolidated Grants programme (grant number ST/K001205/1), and was further supported in this work by the STFC SWIFT-HEP programme, ST/Y005562/1.

• FK acknowledges funding by STFC under grant agreement ST/P006744/1.

• PM has been supported by the Swiss National Science Foundation (SNSF) under contract 240015 and by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme grant agreement 101019620 (ERC Advanced Grant TOPUP).

• VAN acknowledges a PhD studentship funded by the STFC, grant number ST/Y509474/1, and the sub-department of particle physics, University of Oxford.

Observable selection and weights for the MPI tunes

Histograms	Observable	Condition	Bins	Selected range	Weight
H1_2002_I581409 (14 histograms, 67 bins selected)
d01-x01-y01	$M_{jj}$		6/6	[45, ?]	1
d02-x01-y01	$\langle E_{\mathrm{T}}\rangle$		4/4	[20, ?]	1
d03-x01-y01	$E_{\mathrm{T}}^{\max}$		4/4	[25, ?]	1
d04 (2)	$\bar{\eta}$	$25 < E_{\mathrm{T}} < 35\,\mathrm{GeV}$ / $35 < E_{\mathrm{T}} < 80\,\mathrm{GeV}$, $0.1 < y < 0.5$	9/9	[0.6, 2.5]	1
d05 (2)	$\bar{\eta}$	$25 < E_{\mathrm{T}} < 35\,\mathrm{GeV}$ / $35 < E_{\mathrm{T}} < 80\,\mathrm{GeV}$, $0.5 < y < 0.9$	10/10	[0, 2.5]	1
d06 (2)	$x_{\gamma}^{\mathrm{obs}}$		6/8	[0.1, 0.85]	1
d07-x01-y01	$x_{p}^{\mathrm{obs}}$	low $x_{\gamma}^{\mathrm{obs}}$	5/5	[0.05, 0.6]	1
d07-x01-y02	$x_{p}^{\mathrm{obs}}$	high $x_{\gamma}^{\mathrm{obs}}$	5/5	[0.05, 0.6]	0
d08 (2)	$x_{p}^{\mathrm{obs}}$	low $E_{\mathrm{T}}$ / high $E_{\mathrm{T}}$	8/10	[0.1, 0.85]	1
d09-x01-y01	$|\cos\theta^{*}|$	all $M_{jj}$, resolved	7/7	[0, 0.7]	1
d09-x01-y02	$|\cos\theta^{*}|$	all $M_{jj}$, direct	7/7	[0, 0.7]	0
d10-x01-y01	$|\cos\theta^{*}|$	$M_{jj} > 65\,\mathrm{GeV}$, resolved	8/8	[0, 0.85]	1
d10-x01-y02	$|\cos\theta^{*}|$	$M_{jj} > 65\,\mathrm{GeV}$, direct	8/8	[0, 0.85]	0
ZEUS_2007_I753991 (10 histograms, 62 bins selected)
d01-x01-y01	$E_{\mathrm{T}}$	direct	8/8	[17.5, 90.5]	0
d02-x01-y01	$E_{\mathrm{T}}$	resolved	7/7	[17.5, 70.5]	1
d03-x01-y01	$E_{\mathrm{T}}^{\mathrm{jet1}}$	direct	7/7	[20, ?]	0
d04-x01-y01	$E_{\mathrm{T}}^{\mathrm{jet1}}$	resolved	6/6	[20, 70]	1
d05-x01-y01	$\bar{\eta}$	direct	6/6	[$-0.5$, 2.5]	0
d06-x01-y01	$\bar{\eta}$	resolved	6/6	[0, 2.75]	1
d07-x01-y01	$x_{p}^{\mathrm{obs}}$	direct	11/11	[0, 1]	0
d08-x01-y01	$x_{p}^{\mathrm{obs}}$	resolved	11/11	[0, 1]	1
d09-x01-y01	$\Delta\phi$	direct	5/5	[1.83, 3.14]	0
d10-x01-y01	$\Delta\phi$	resolved	7/7	[0, 3.14]	1
d11–d14 (4)	$x_{p}^{\mathrm{obs}}$	direct	14/14	[0, 0.5]	0
d15–d18 (4)	$x_{p}^{\mathrm{obs}}$	resolved	18/18	[0, 1]	1
d19-x01-y01	$x_{\gamma}^{\mathrm{obs}}$		7/8	[0.1, 0.8]	1
ZEUS_1997_I450085 (0 histograms, 0 bins selected)
d01–d04 (4)	$\bar{\eta}$	$E_{\mathrm{T}}^{\mathrm{jet}} > 6\,\mathrm{GeV}$ / $E_{\mathrm{T}}^{\mathrm{jet}} > 8\,\mathrm{GeV}$ / $E_{\mathrm{T}}^{\mathrm{jet}} > 11\,\mathrm{GeV}$ / $E_{\mathrm{T}}^{\mathrm{jet}} > 15\,\mathrm{GeV}$, all	29/29	[$-1.125$, 1.625]	0
d21–d24 (4)	$\bar{\eta}$	$E_{\mathrm{T}}^{\mathrm{jet}} > 6\,\mathrm{GeV}$ / $E_{\mathrm{T}}^{\mathrm{jet}} > 8\,\mathrm{GeV}$ / $E_{\mathrm{T}}^{\mathrm{jet}} > 11\,\mathrm{GeV}$ / $E_{\mathrm{T}}^{\mathrm{jet}} > 15\,\mathrm{GeV}$, direct	25/25	[$-1.125$, 1.625]	0
d25–d28 (4)	$\bar{\eta}$	$E_{\mathrm{T}}^{\mathrm{jet}} > 6\,\mathrm{GeV}$ / $E_{\mathrm{T}}^{\mathrm{jet}} > 8\,\mathrm{GeV}$ / $E_{\mathrm{T}}^{\mathrm{jet}} > 11\,\mathrm{GeV}$ / $E_{\mathrm{T}}^{\mathrm{jet}} > 15\,\mathrm{GeV}$, mid	22/22	[$-0.375$, 1.625]	0

Histograms	Observable	Condition	Bins	Selected range	Weight
ZEUS_2001_I568665 (11 histograms, 56 bins selected)
d01-x01-y01	$\lvert\cos\theta^*\rvert$	resolved	4/8	[0.4, 0.8]	1
d01-x01-y02	$\lvert\cos\theta^*\rvert$	direct	8/8	[0, 0.8]	0
d02–d07 (6)	$E_{\mathrm{T}}^{\mathrm{jet1}}$	direct, $\eta$ bin	38/38	[14, 75]	0
d08–d13 (6)	$E_{\mathrm{T}}^{\mathrm{jet1}}$	resolved, $\eta$ bin	32/32	[14, 65]	1
d14–d19 (6)	$\eta^{\mathrm{jet2}}$	direct, $\eta$ bin	41/41	[$-1$, 2.4]	0
d20–d23 (4)	$x_{\gamma}^{\mathrm{obs}}$	$14 < E_{\mathrm{T}}^{\mathrm{jet1}} < 17\,\mathrm{GeV}$ / $17 < E_{\mathrm{T}}^{\mathrm{jet1}} < 25\,\mathrm{GeV}$ / $25 < E_{\mathrm{T}}^{\mathrm{jet1}} < 35\,\mathrm{GeV}$ / $35 < E_{\mathrm{T}}^{\mathrm{jet1}} < 90\,\mathrm{GeV}$	20/24	[0.1, 0.8]	1
ZEUS_2007_I756660 (19 histograms, 192 bins selected)
d01-x01-y01	$M_{3j}$		11/11	[25, 160]	1
d02-x01-y01	$M_{4j}$		11/11	[25, 160]	0
d03 (2)	$x_{\gamma}^{\mathrm{obs}}$	$25 < M_{3j} < 50\,\mathrm{GeV}$ / $M_{3j} > 50\,\mathrm{GeV}$	13/17	[0, 0.85]	1
d04 (2)	$x_{\gamma}^{\mathrm{obs}}$	low $M_{4j}$ / high $M_{4j}$	17/17	[0, 1]	0
d05 (2)	$y$	low $M_{3j}$ / high $M_{3j}$	18/18	[0.2, 0.85]	1
d06 (2)	$y$	low $M_{4j}$ / high $M_{4j}$	18/18	[0.2, 0.85]	0
d07 (2)	$E_{\mathrm{T}}^{\mathrm{jet1}}$	low $M_{3j}$ / high $M_{3j}$	19/19	[6, 71]	1
d08 (2)	$E_{\mathrm{T}}^{\mathrm{jet2}}$	low $M_{3j}$ / high $M_{3j}$	16/16	[6 60.8]	1
d09 (2)	$E_{\mathrm{T}}^{\mathrm{jet3}}$	low $M_{3j}$ / high $M_{3j}$	12/12	[6, 40]	1
d10 (2)	$E_{\mathrm{T}}^{\mathrm{jet1}}$	low $M_{4j}$ / high $M_{4j}$	16/16	[6 60.8]	0
d11 (2)	$E_{\mathrm{T}}^{\mathrm{jet2}}$	low $M_{4j}$ / high $M_{4j}$	13/13	[6 43.9]	0
d12 (2)	$E_{\mathrm{T}}^{\mathrm{jet3}}$	low $M_{4j}$ / high $M_{4j}$	9/9	[6, 26]	0
d13 (2)	$E_{\mathrm{T}}^{\mathrm{jet4}}$	low $M_{4j}$ / high $M_{4j}$	6/6	[6, 20]	0
d14 (2)	$\eta^{\mathrm{jet1}}$	low $M_{3j}$ / high $M_{3j}$	26/26	[$-1.686$, 2.4]	1
d15 (2)	$\eta^{\mathrm{jet2}}$	low $M_{3j}$ / high $M_{3j}$	27/27	[$-2$, 2.4]	1
d16 (2)	$\eta^{\mathrm{jet3}}$	low $M_{3j}$ / high $M_{3j}$	28/28	[$-2$, 2.4]	1
d17 (2)	$\eta^{\mathrm{jet1}}$	low $M_{4j}$ / high $M_{4j}$	24/24	[$-1.371$, 2.4]	0
d18 (2)	$\eta^{\mathrm{jet2}}$	low $M_{4j}$ / high $M_{4j}$	25/25	[$-1.686$, 2.4]	0
d19 (2)	$\eta^{\mathrm{jet3}}$	low $M_{4j}$ / high $M_{4j}$	25/25	[$-1.686$, 2.4]	0
d20 (2)	$\eta^{\mathrm{jet4}}$	low $M_{4j}$ / high $M_{4j}$	25/25	[$-1.686$, 2.4]	0
d21 (2)	$\cos\psi_3$	low $M_{3j}$ / high $M_{3j}$	22/22	[$-1$, 1]	1
ZEUS_2012_I1116258 (10 histograms, 94 bins selected)
d01-x01-y01	$E_{\mathrm{T}}$	$\eta$ bin	9/9	[17, 95]	1
d02–d03 (2)	$\eta$	$E_{\mathrm{T}} > 17\,\mathrm{GeV}$ / $E_{\mathrm{T}} > 21\,\mathrm{GeV}$	26/26	[$-0.75$, 2.5]	1
d04–d08 (5)	$E_{\mathrm{T}}$	$\eta$ bin	41/41	[17, 95]	1
d09–d10 (2)	$E_{\mathrm{T}}$	anti-$k_{\mathrm{t}}$ / SIScone	18/18	[17, 95]	1
d11–d12 (2)	$\eta$	$E_{\mathrm{T}} > 17\,\mathrm{GeV}$, anti-$k_{\mathrm{t}}$ / SIScone	26/26	[$-0.75$, 2.5]	0

Histograms	Observable	Condition	Bins	Selected range	Weight
H1_2006_I711847 (9 histograms, 43 bins selected)
d01-x01-y01	$\lvert\cos\theta^*\rvert$	resolved	8/8	[0, 0.85]	1
d01-x01-y02	$\lvert\cos\theta^*\rvert$	direct	8/8	[0, 0.85]	0
d02-x01-y01	$\lvert\cos\theta^*\rvert$	$M_{jj}$ cut, resolved	8/8	[0, 0.85]	1
d02-x01-y02	$\lvert\cos\theta^*\rvert$	$M_{jj}$ cut, direct	8/8	[0, 0.85]	0
d03-x01-y01	$x_{\gamma}^{\mathrm{obs}}$		4/4	[0.3, 1]	0
d03-x01-y02	$x_{\gamma}^{\mathrm{obs}}$		4/5	[0.1, 0.85]	1
d04–d06 (3)	$x_{p}^{\mathrm{obs}}$	resolved	8/13	[0.05, 0.32]	1
d07–d09 (3)	$x_{p}^{\mathrm{obs}}$	direct	13/13	[0.05, 0.7]	0
d10–d12 (3)	$E_{\mathrm{T}}^{\max}$	resolved	15/15	[25, 80]	1
d13–d15 (3)	$E_{\mathrm{T}}^{\max}$	direct	16/16	[25, 80]	0

Histograms	Observable	Condition	Bins	Selected range	Weight
L3_2004_I661114 (1 histogram, 9 bins selected)
d01-x01-y01	$p_{\mathrm{T}}^{\mathrm{jet}}$	$\lvert\eta\rvert< 1$	9/9	[3, 50]	1
OPAL_2008_I754316 (2 histograms, 12 bins selected)
d01–d02 (2)	$p_{\mathrm{T}}^{\mathrm{jet}}$	$\lvert\eta\rvert< 1$ / $\lvert\eta\rvert< 1.5$	12/12	[5, 40]	1
OPAL_2003_I611415 (19 histograms, 113 bins selected)
d01–d02 (2)	$\lvert\cos\theta^*\rvert$	$x_{\gamma}^{\pm} > 0.75$ / $x_{\gamma}^{\pm} < 0.75$	16/16	[0, 0.8]	0
d03–d05 (3)	$\bar{E}_{\mathrm{T}}$	all $x_{\gamma}^{\pm}$ / $x_{\gamma}^{+}$ or $x_{\gamma}^{-} < 0.75$ / $x_{\gamma}^{+}$ and $x_{\gamma}^{-} < 0.75$	15/15	[5, 25]	1
d06–d08 (3)	$x_{\gamma}^{\mathrm{obs}}$	$5 < \bar{E}_{\mathrm{T}} < 7\,\mathrm{GeV}$ / $7 < \bar{E}_{\mathrm{T}} < 11\,\mathrm{GeV}$ / $11 < \bar{E}_{\mathrm{T}} < 25\,\mathrm{GeV}$, all	18/24	[0, 0.75]	1
d09–d10 (2)	$x_{\gamma}^{\mathrm{obs}}$	$5 < \bar{E}_{\mathrm{T}} < 7\,\mathrm{GeV}$ / $7 < \bar{E}_{\mathrm{T}} < 11\,\mathrm{GeV}$, $x_{\gamma}^{+}$ or $x_{\gamma}^{-} < 0.75$	12/16	[0, 0.75]	1
d11–d12 (2)	$x_{\gamma}^{\mathrm{obs}}$	$5 < \bar{E}_{\mathrm{T}} < 7\,\mathrm{GeV}$ / $7 < \bar{E}_{\mathrm{T}} < 11\,\mathrm{GeV}$, $x_{\gamma}^{+}$ and $x_{\gamma}^{-} < 0.75$	12/12	[0, 0.75]	1
d13–d15 (3)	$\ln x_{\gamma}^{\mathrm{obs}}$	$5 < \bar{E}_{\mathrm{T}} < 7\,\mathrm{GeV}$, all / $x_{\gamma}^{+}$ or $x_{\gamma}^{-} < 0.75$ / $x_{\gamma}^{+}$ and $x_{\gamma}^{-} < 0.75$	12/15	[−1.65, −0.33]	1
d16–d17 (2)	$\lvert\Delta\eta\rvert$	$x_{\gamma}^{+}$ or $x_{\gamma}^{-} < 0.75$ / $x_{\gamma}^{+}$ and $x_{\gamma}^{-} < 0.75$	14/14	[0, 3.5]	1
d18–d19 (2)	$\min\lvert\eta\rvert$	$x_{\gamma}^{+}$ or $x_{\gamma}^{-} < 0.75$ / $x_{\gamma}^{+}$ and $x_{\gamma}^{-} < 0.75$	14/14	[0, 1.75]	1
d20–d21 (2)	$\max\lvert\eta\rvert$	$x_{\gamma}^{+}$ or $x_{\gamma}^{-} < 0.75$ / $x_{\gamma}^{+}$ and $x_{\gamma}^{-} < 0.75$	16/16	[0, 2]	1

Observable selection and weights for the PDF tune

Histograms	Observable	Condition	Bins	Selected range	Weight
H1_2002_I581409 (11 histograms, 49 bins selected)
d01-x01-y01	$M_{jj}$		4/6	[70, 180]	1
d02-x01-y01	$\langle E_{\mathrm{T}}\rangle$		3/4	[30, 80]	1
d03-x01-y01	$E_{\mathrm{T}}^{\max}$		4/4	[25, 80]	1
d04-x01-y01	$\bar{\eta}$	$25 < E_{\mathrm{T}} < 35\,\mathrm{GeV}$, $0.1 < y < 0.5$	5/5	[0.6, 2.5]	0
d04-x01-y02	$\bar{\eta}$	$35 < E_{\mathrm{T}} < 80\,\mathrm{GeV}$, $0.1 < y < 0.5$	4/4	[0.9, 2.5]	1
d05-x01-y01	$\bar{\eta}$	$25 < E_{\mathrm{T}} < 35\,\mathrm{GeV}$, $0.5 < y < 0.9$	6/6	[0, 2.5]	0
d05-x01-y02	$\bar{\eta}$	$35 < E_{\mathrm{T}} < 80\,\mathrm{GeV}$, $0.5 < y < 0.9$	4/4	[0.6, 2.1]	1
d06 (2)	$x_{\gamma}^{\mathrm{obs}}$		6/8	[0.1, 0.85]	1
d07-x01-y01	$x_{p}^{\mathrm{obs}}$	low $x_{\gamma}^{\mathrm{obs}}$	5/5	[0.05, 0.6]	1
d07-x01-y02	$x_{p}^{\mathrm{obs}}$	high $x_{\gamma}^{\mathrm{obs}}$	5/5	[0.05, 0.6]	0
d08-x01-y01	$x_{\gamma}^{\mathrm{obs}}$	low $E_{\mathrm{T}}$	4/5	[0.1, 0.85]	0
d08-x01-y02	$x_{\gamma}^{\mathrm{obs}}$	high $E_{\mathrm{T}}$	4/5	[0.1, 0.85]	1
d09-x01-y01	$\lvert\cos\theta^{*}\rvert$	all $M_{jj}$, resolved	7/7	[0, 0.7]	1
d09-x01-y02	$\lvert\cos\theta^{*}\rvert$	all $M_{jj}$, direct	7/7	[0, 0.7]	0
d10-x01-y01	$\lvert\cos\theta^{*}\rvert$	$M_{jj} > 65\,\mathrm{GeV}$, resolved	8/8	[0, 0.85]	1
d10-x01-y02	$\lvert\cos\theta^{*}\rvert$	$M_{jj} > 65\,\mathrm{GeV}$, direct	8/8	[0, 0.85]	0
ZEUS_2007_I753991 (10 histograms, 59 bins selected)
d01-x01-y01	$E_{\mathrm{T}}$	direct	8/8	[17.5, 90.5]	0
d02-x01-y01	$E_{\mathrm{T}}$	resolved	5/7	[27.5, 70.5]	1
d03-x01-y01	$E_{\mathrm{T}}^{\mathrm{jet1}}$	direct	7/7	[20, 90]	0
d04-x01-y01	$E_{\mathrm{T}}^{\mathrm{jet1}}$	resolved	5/6	[26, 70]	1
d05-x01-y01	$\bar{\eta}$	direct	6/6	[$-0.5$, 2.5]	0
d06-x01-y01	$\bar{\eta}$	resolved	6/6	[0, 2.75]	1
d07-x01-y01	$x_{p}^{\mathrm{obs}}$	direct	11/11	[0, 1]	0
d08-x01-y01	$x_{p}^{\mathrm{obs}}$	resolved	11/11	[0, 1]	1
d09-x01-y01	$\Delta\phi$	direct	5/5	[1.83, 3.14]	0
d10-x01-y01	$\Delta\phi$	resolved	7/7	[0, 3.14]	1
d11–d14 (4)	$x_{\gamma}^{\mathrm{obs}}$	direct	14/14	[0, 0.5]	0
d15–d18 (4)	$x_{p}^{\mathrm{obs}}$	resolved	18/18	[0, 1]	1
d19-x01-y01	$x_{\gamma}^{\mathrm{obs}}$		7/8	[0.1, 0.8]	1
ZEUS_1997_I450085 (4 histograms, 14 bins selected)
d01–d02 (2)	$\bar{\eta}$	$E_{\mathrm{T}}^{\mathrm{jet}} > 6\,\mathrm{GeV}$ / $E_{\mathrm{T}}^{\mathrm{jet}} > 8\,\mathrm{GeV}$, all $x_{\gamma}^{\mathrm{obs}}$	21/21	[$-1.125$, 1.625]	0
d03–d04 (2)	$\bar{\eta}$	$E_{\mathrm{T}}^{\mathrm{jet}} > 11\,\mathrm{GeV}$ / $E_{\mathrm{T}}^{\mathrm{jet}} > 15\,\mathrm{GeV}$, all $x_{\gamma}^{\mathrm{obs}}$	8/8	[$-0.375$, 1.625]	1
d21–d24 (4)	$\bar{\eta}$	$E_{\mathrm{T}}^{\mathrm{jet}} > 6\,\mathrm{GeV}$ / $E_{\mathrm{T}}^{\mathrm{jet}} > 8\,\mathrm{GeV}$ / $E_{\mathrm{T}}^{\mathrm{jet}} > 11\,\mathrm{GeV}$ / $E_{\mathrm{T}}^{\mathrm{jet}} > 15\,\mathrm{GeV}$, direct	25/25	[$-1.125$, 1.625]	0
d25–d26 (2)	$\bar{\eta}$	$E_{\mathrm{T}}^{\mathrm{jet}} > 6\,\mathrm{GeV}$ / $E_{\mathrm{T}}^{\mathrm{jet}} > 8\,\mathrm{GeV}$, mid $x_{\gamma}^{\mathrm{obs}}$	16/16	[$-0.375$, 1.625]	0
d27–d28 (2)	$\bar{\eta}$	$E_{\mathrm{T}}^{\mathrm{jet}} > 11\,\mathrm{GeV}$ / $E_{\mathrm{T}}^{\mathrm{jet}} > 15\,\mathrm{GeV}$, mid $x_{\gamma}^{\mathrm{obs}}$	6/6	[0.125, 1.625]	1

Histograms	Observable	Condition	Bins	Selected range	Weight
ZEUS_2001_I568665 (9 histograms, 38 bins selected)
d01-x01-y01	$\lvert\cos\theta^{*}\rvert$	resolved	4/8	[0.4, 0.8]	1
d01-x01-y02	$\lvert\cos\theta^{*}\rvert$	direct	8/8	[0, 0.8]	0
d02–d07 (6)	$E_{\mathrm{T}}^{\mathrm{jet1}}$	direct, $\eta$ bin	38/38	[14, 75]	0
d08-x01-y01	$E_{\mathrm{T}}^{\mathrm{jet1}}$	resolved, $\eta$ bin	0/1	none	1
d09–d13 (5)	$E_{\mathrm{T}}^{\mathrm{jet1}}$	resolved, $\eta$ bin	21/31	[21, 65]	1
d14–d19 (6)	$\eta^{\mathrm{jet2}}$	direct, $\eta$ bin	41/41	[$-1$, 2.4]	0
d20-x01-y01	$x_{\gamma}^{\mathrm{obs}}$	$14 < E_{\mathrm{T}}^{\mathrm{jet1}} < 17\,\mathrm{GeV}$	7/8	[0.1, 0.8]	0
d21–d23 (3)	$x_{\gamma}^{\mathrm{obs}}$	$17 < E_{\mathrm{T}}^{\mathrm{jet1}} < 25\,\mathrm{GeV}$ / $25 < E_{\mathrm{T}}^{\mathrm{jet1}} < 35\,\mathrm{GeV}$ / $35 < E_{\mathrm{T}}^{\mathrm{jet1}} < 90\,\mathrm{GeV}$	13/16	[0.1, 0.8]	1
ZEUS_2007_I756660 (40 histograms, 377 bins selected)
d01-x01-y01	$M_{3j}$		11/11	[25, 160]	1
d02-x01-y01	$M_{4j}$		11/11	[25, 160]	$1^{\dagger}$
d03–d04 (4)	$x_{\gamma}^{\mathrm{obs}}$	$25 < M_{3j} < 50\,\mathrm{GeV}$ / $M_{3j} > 50\,\mathrm{GeV}$ / low $M_{4j}$ / high $M_{4j}$	26/34	[0, 0.85]	$1^{\dagger}$
d05–d06 (4)	$y$	low $M_{3j}$ / high $M_{3j}$ / low $M_{4j}$ / high $M_{4j}$	36/36	[0.2, 0.85]	$1^{\dagger}$
d07 (2)	$E_{\mathrm{T}}^{\mathrm{jet1}}$	low $M_{3j}$ / high $M_{3j}$	19/19	[6, 71]	1
d08 (2)	$E_{\mathrm{T}}^{\mathrm{jet2}}$	low $M_{3j}$ / high $M_{3j}$	16/16	[6 60.8]	1
d09 (2)	$E_{\mathrm{T}}^{\mathrm{jet3}}$	low $M_{3j}$ / high $M_{3j}$	12/12	[6, 40]	1
d10 (2)	$E_{\mathrm{T}}^{\mathrm{jet1}}$	low $M_{4j}$ / high $M_{4j}$	16/16	[6 60.8]	$1^{\dagger}$
d11 (2)	$E_{\mathrm{T}}^{\mathrm{jet2}}$	low $M_{4j}$ / high $M_{4j}$	13/13	[6 43.9]	$1^{\dagger}$
d12 (2)	$E_{\mathrm{T}}^{\mathrm{jet3}}$	low $M_{4j}$ / high $M_{4j}$	9/9	[6, 26]	$1^{\dagger}$
d13 (2)	$E_{\mathrm{T}}^{\mathrm{jet4}}$	low $M_{4j}$ / high $M_{4j}$	6/6	[6, 20]	$1^{\dagger}$
d14 (2)	$\eta^{\mathrm{jet1}}$	low $M_{3j}$ / high $M_{3j}$	26/26	[$-1.686$, 2.4]	1
d15 (2)	$\eta^{\mathrm{jet2}}$	low $M_{3j}$ / high $M_{3j}$	27/27	[$-2$, 2.4]	1
d16 (2)	$\eta^{\mathrm{jet3}}$	low $M_{3j}$ / high $M_{3j}$	28/28	[$-2$, 2.4]	1
d17 (2)	$\eta^{\mathrm{jet1}}$	low $M_{4j}$ / high $M_{4j}$	24/24	[$-1.371$, 2.4]	$1^{\dagger}$
d18 (2)	$\eta^{\mathrm{jet2}}$	low $M_{4j}$ / high $M_{4j}$	25/25	[$-1.686$, 2.4]	$1^{\dagger}$
d19 (2)	$\eta^{\mathrm{jet3}}$	low $M_{4j}$ / high $M_{4j}$	25/25	[$-1.686$, 2.4]	$1^{\dagger}$
d20 (2)	$\eta^{\mathrm{jet4}}$	low $M_{4j}$ / high $M_{4j}$	25/25	[$-1.686$, 2.4]	$1^{\dagger}$
d21 (2)	$\cos\psi_{3}$	low $M_{3j}$ / high $M_{3j}$	22/22	[$-1$, 1]	1
ZEUS_2012_I1116258 (9 histograms, 57 bins selected)
d01-x01-y01	$E_{\mathrm{T}}$	$\eta$ bin	6/9	[29, 95]	1
d02-x01-y01	$\eta$	$E_{\mathrm{T}} > 17\,\mathrm{GeV}$	13/13	[$-0.75$, 2.5]	0
d03-x01-y01	$\eta$	$E_{\mathrm{T}} > 21\,\mathrm{GeV}$	13/13	[$-0.75$, 2.5]	1
d04–d08 (5)	$E_{\mathrm{T}}$	$\eta$ bin	26/41	[29, 95]	1
d09–d10 (2)	$E_{\mathrm{T}}$	anti-$k_{\mathrm{t}}$ / SIScone	12/18	[29, 95]	1
d11–d12 (2)	$\eta$	$E_{\mathrm{T}} > 17\,\mathrm{GeV}$, anti-$k_{\mathrm{t}}$ / SIScone	26/26	[$-0.75$, 2.5]	0