1 Introduction

The proton structure at high momentum-transfer, as encoded in the collinearly factorized parton distribution functions (PDFs), is not only an interesting subject in its own right, but plays a pivotal role in many applications, such as precision electroweak and Higgs physics, searches for new physics, etc. [1]. Likewise, their counterparts for nucleons bound in nuclei, the nuclear PDFs (nPDFs), are essential in e.g. studying the production of hard probes of the Quark Gluon Plasma [2]. In practice, despite the ongoing effort in lattice methods [3], the PDFs are obtained by the well-established means of global analysis using hard-process data. As such, the PDFs have uncertainties which derive from those in the available data and also from the lack of data constraints in certain phase-space regions. It is then often the case that when new data are published or a future experiment is being planned, one would like to study the impact that the measurement could have on the PDFs. A good example of such a case is the recent CMS measurement of dijet pseudorapidity spectra in proton–proton (pp) and proton–lead (pPb) collisions at 5.02 TeV [4], where, on one hand, the measured pp spectra seem to be in a disagreement with next-to-leading order (NLO) perturbative QCD (pQCD) calculations using CT14 [5] and MMHT14 [6] PDFs (see the Supplemental Material of Ref. [4]), while, on the other hand, the nuclear-modification ratio of the pPb and pp spectra appear to have much smaller uncertainties than predictions with various nPDFs. One should therefore study the impact these data could have on both the free-proton PDFs and their nuclear modifications.

As producing a full global fit remains rather involved, even with publicly available tools like the xFitter [7] (built upon the former HERAFitter [8]) coming available, it is in most cases impractical for a general user to try to learn about the constraining power of a data set in this way. For this purpose, approximative methods have been developed, first in the formalism of Bayesian reweighting of Monte Carlo PDF ensembles [9,10,11,12,13] and later in a framework using Hessian error sets [14,15,16]. These methods have their limitations, as the new PDFs rely on all the theoretical assumptions of the original PDF analysis, such as the parametrization form, the value of \(\alpha _{\mathrm{s}}\) and the used heavy-quark scheme. There are also limitations related to how well the methods approximate the true parameter likelihood in the region constrained by the new data. In particular, the applications of Hessian PDF reweighting have resorted to a perfectly quadratic approximation of the \(\chi ^2\) goodness-of-fit function before inclusion of the new data and to linear or up to quadratic terms for responses in the new observables. This applies to the implementation in the xFitter package, where the method is referred to as “Hessian profiling”, as well as to the new software package which has appeared under the name ePump [16]. It is not, however, uncommon that non-quadratic terms in the \(\chi ^2\) function are large (see e.g. Figure 6 of Ref. [17]), and thus it would be beneficial to have a way to take these into account.

The purpose of this article is twofold: first, in Sect. 2, we describe how one can include into Hessian PDF reweighting the first non-quadratic terms in the \(\chi ^2\) function consistently with the original fit, provided that the needed information is available. Second, in Sect. 3, we apply the Hessian PDF reweighting to the aforementioned CMS dijet measurements at 5.02 TeV [4]. We show that the strong disagreement between the pp measurement and next-to-leading order (NLO) calculations using CT14 NLO PDFs [5] can be brought to a much better agreement upon reweighting the CT14 PDFs, but that this requires rather strong modifications for high-x gluons. We demonstrate that such changes in the proton PDFs have also an important impact on predictions for dijet production in pPb. Finally, we then reweight the EPPS16 nPDFs [18] with the nuclear modification ratio of the measured pPb and pp dijet spectra using the non-quadratic approximation developed in Sect. 2 and present a discussion on the importance of these higher-order terms in the reweighting. Preliminary work on this topic can be found in Refs. [19, 20].

2 PDF uncertainties and reweighting in Hessian method

In this section, we first recapitulate the uncertainty determination in the Hessian approach [21], assuming the use of a global tolerance criterion. We then describe how one can perform a reweighting upon such determined error sets, taking into account the first non-quadratic terms in the \(\chi ^2\) function. We end the section with a discussion on the applicability of this method in the case of non-global tolerances.

2.1 Hessian uncertainties with global tolerance criterion

In PDF global analyses, the goodness-of-fit of a parameter vector \(\mathbf {a}\) is dictated by the \(\chi ^2\) function

$$\begin{aligned} \chi ^2(\mathbf {a}) = \sum _{ij}\,(y_i(\mathbf {a}) - y_i^\text {data})\,C^{-1}_{ij}\,(y_j(\mathbf {a}) - y_j^\text {data}), \end{aligned}$$
(1)

where \(y_i(\mathbf {a})\) are theory predictions for the observables included in the analysis, \(y_i^\text {data}\) the corresponding measured values and \(C^{-1}_{ij}\) the elements of the inverse covariance matrix for these data. In the Hessian method for uncertainty estimation, one takes the parameter values \(\mathbf {a}^\text {min}\) which minimize Eq. (1), \(\chi ^2(\mathbf {a}^\text {min}) \equiv \min \chi ^2(\mathbf {a}) \equiv \chi ^2_0\), as the central, best-fit values and studies the behaviour of the \(\chi ^2\) function around this minimum to determine the uncertainty in these parameters.

The leading deviations from the minimum value \(\chi ^2_0\) are given by the quadratic approximation

$$\begin{aligned} \chi ^2 \approx \chi ^2_0 + \sum _{ij}\,(a_i - a_i^\text {min})\,H_{ij}\,(a_j - a_j^\text {min}), \end{aligned}$$
(2)

where \(H_{ij} = \frac{1}{2} \partial ^2 \chi ^2 / \partial a_i \partial a_j|_{\mathbf {a}=\mathbf {a}^\text {min}}\) are the elements of the Hessian matrix. In practice, these elements need to be obtained numerically. Since the Hessian matrix is symmetric, it has a complete set of orthonormal eigenvectors \(\mathbf {v}^{(k)}\) such that

$$\begin{aligned} \sum _j\,H_{ij}\,v^{(k)}_j= & {} \varepsilon _k\,v^{(k)}_i, \end{aligned}$$
(3)
$$\begin{aligned} \sum _i\,v^{(k)}_i\,v^{(\ell )}_i= & {} \delta _{k\ell }, \qquad \sum _k\,v^{(k)}_i\,v^{(k)}_j = \delta _{ij}, \end{aligned}$$
(4)

where \(\varepsilon _k\) are the eigenvalues of the Hessian matrix. With this eigendecomposition we can define new parameters

$$\begin{aligned} z_k = \sum _i \sqrt{\varepsilon _k}\,v^{(k)}_i\,(a_i - a_i^\text {min}) \end{aligned}$$
(5)

such that Eq. (2) becomes

$$\begin{aligned} \chi ^2 \approx \chi ^2_0 + \sum _{k}\,z_k^2. \end{aligned}$$
(6)

Since the new parameters \(z_k\) are uncorrelated in the quadratic approximation, one can use the standard law of error propagation to translate the uncertainties in the parameters \(z_k\) to the uncertainty of any PDF-dependent quantity X as [21]

$$\begin{aligned} \varDelta X = \sqrt{\sum _k \left( \frac{\partial X}{\partial z_k} \varDelta z_k\right) ^2}. \end{aligned}$$
(7)

Given a well justified global tolerance \(\varDelta \chi ^2\) for the allowed growth of \(\chi ^2\) from its minimum, one can determine the allowed parameter variations \(\varDelta z_k\).Footnote 1 If the \(\chi ^2\) function were perfectly quadratic, the uncertainty of the parameter \(z_k\) corresponding to the tolerance \(\varDelta \chi ^2\) would be simply \(\varDelta z_k = \sqrt{\varDelta \chi ^2}\). As this is generally not true, one instead finds \(\delta z^\pm _k\), the positive and negative values of \(z_k\) corresponding to the \(\varDelta \chi ^2\) increase, and assigns \(\varDelta z_k = (\delta z^+_k - \delta z^-_k) / 2\). It is convenient to define error sets \(S^\pm _i\) corresponding to parameter values

$$\begin{aligned} z_k[S^\pm _i] = {\left\{ \begin{array}{ll} \delta z^\pm _i, &{} k = i \\ 0, &{} k \ne i \end{array}\right. } \quad , \end{aligned}$$
(8)

along with the central set \(S_0\), where \(z_k[S_0] = 0\) for all k. Estimating

$$\begin{aligned} \frac{\partial X}{\partial z_k} = \frac{X[S^+_k] - X[S^-_k]}{2\,\varDelta z_k}, \end{aligned}$$
(9)

where \(X[S^\pm _k]\) stands for the quantity X calculated with the parameter set of Eq. (8), yields then a simple form

$$\begin{aligned} \varDelta X = \frac{1}{2}\sqrt{\sum _k \left( X[S^+_k] - X[S^-_k]\right) ^2}. \end{aligned}$$
(10)

As the response in X to the upward and downward parameter shifts can be uneven, one can alternatively specify an upward–downward asymmetric error prescription e.g. with [23]

$$\begin{aligned} \delta X^\pm = \sqrt{\sum _k \left[ \begin{array}{c} \max \\ \min \end{array} \left\{ X[S^+_k]-X[S_0],X[S^-_k]-X[S_0],0 \right\} \right] ^2}. \end{aligned}$$
(11)

2.2 Non-quadratic reweighting

Fig. 1
figure 1

An illustration for the response of \(\chi ^2\) (top) and \(y_i\) (bottom) with respect to a change of parameter \(z_k\) in quadratic–linear (red, long dashed), quadratic–quadratic (blue, short dashed) and cubic–quadratic (black, solid) approximations

Full size image

In the presence of a new data set, the total \(\chi ^2\) can be written as

$$\begin{aligned} \chi ^2_\text {new}(\mathbf {z}) = \chi ^2_\text {old}(\mathbf {z}) + \sum _{ij}\,(y_i(\mathbf {z}) - y_i^\text {data})\,C^{-1}_{ij}\,(y_j(\mathbf {z}) - y_j^\text {data}), \end{aligned}$$
(12)

where \(y_i\) (\(y_i^\text {data}\)) now correspond to the new theoretical (measured) values and \(\chi ^2_\text {old}\) incorporates our knowledge of the original global analysis. Now, as we do not wish to produce a full global analysis with \(\chi ^2_\text {new}\), we need to make suitable approximations. The simplest choice is to use the quadratic approximation in Eq. (6), according to the method introduced in Ref. [15], but if the parameter variations \(\delta z^\pm _k\) and the global tolerance \(\varDelta \chi ^2\) of this fit are known (as is the case with EPPS16 nPDFs, see Table 2 in Ref. [18]), then \(\chi ^2_\text {old}\) can be approximated with a third order polynomial in each of the eigendirections,

$$\begin{aligned} \chi ^2_\text {old} \approx \chi ^2_0 + \sum _{k} (a_k z_k^2 + b_k z_k^3), \end{aligned}$$
(13)

where the coefficients are obtained with

$$\begin{aligned} a_k= & {} \frac{\varDelta \chi ^2}{\delta z^+_k - \delta z^-_k}\left( \frac{\delta z^+_k}{(\delta z^-_k)^2} - \frac{\delta z^-_k}{(\delta z^+_k)^2}\right) , \end{aligned}$$
(14)
$$\begin{aligned} b_k= & {} \frac{\varDelta \chi ^2}{\delta z^+_k - \delta z^-_k}\left( \frac{1}{(\delta z^+_k)^2} - \frac{1}{(\delta z^-_k)^2}\right) . \end{aligned}$$
(15)

This is illustrated in Fig. 1 (upper diagram), where we show an example of a situation where the \(\chi ^2\) grows asymmetrically with respect to \(z_k\). The quadratic approximation fails to acknowledge this fact and a third order polynomial is needed to reproduce the \(\varDelta \chi ^2\) growth at \(\delta z_k^-\) and \(\delta z_k^+\). Similarly, as illustrated in Fig. 1 (lower diagram), the \(y_i\) can be expanded in terms of \(z_k\) as

$$\begin{aligned} y_i(\mathbf {z}) \approx y_i[S_0] + \sum _k (d_{ik} z_k + e_{ik} z_k^2), \end{aligned}$$
(16)

where

$$\begin{aligned} d_{ik}= & {} \frac{1}{\delta z^+_k - \delta z^-_k}\bigg [-\frac{\delta z^-_k}{\delta z^+_k}\left( y_i[S_k^+] - y_i[S_0]\right) \nonumber \\&+ \frac{\delta z^+_k}{\delta z^-_k}\left( y_i[S_k^-] - y_i[S_0]\right) \bigg ], \end{aligned}$$
(17)
$$\begin{aligned} e_{ik}= & {} \frac{1}{\delta z^+_k - \delta z^-_k}\bigg [\frac{1}{\delta z^+_k}\left( y_i[S_k^+] - y_i[S_0]\right) \nonumber \\&- \frac{1}{\delta z^-_k}\left( y_i[S_k^-] - y_i[S_0]\right) \bigg ]. \end{aligned}$$
(18)

One should note that the above approximations do not yield a full Taylor expansion to cubic and quadratic order in \(\chi ^2_\text {old}\) and \(y_i(\mathbf {z})\), respectively, as we have neglected off-diagonal terms proportional to \(z_l z_k^2\) and \(z_l z_k\) for \(l \ne k\). Even so, we will refer to reweighting with these approximations as a cubic–quadratic one.

Changing variables to \(w_k = 2z_k / (\delta z^+_k - \delta z^-_k)\) and defining \(r_k = -\delta z^+_k/\delta z^-_k\), we may alternatively write

$$\begin{aligned} \begin{aligned} \chi ^2_\text {new}(\mathbf {w}) - \chi ^2_0 \approx&\sum _{k} (A_k w_k^2 + B_k w_k^3)\\&+ \sum _{ij}\,(y_i(\mathbf {w}) - y_i^\text {data})\,C^{-1}_{ij}\,(y_j(\mathbf {w}) - y_j^\text {data}), \end{aligned} \end{aligned}$$
(19)

where

$$\begin{aligned} A_k= & {} \frac{\varDelta \chi ^2}{4}\left( \frac{1}{r_k^2} + \frac{1}{r_k} + r_k + r_k^2\right) , \end{aligned}$$
(20)
$$\begin{aligned} B_k= & {} \frac{\varDelta \chi ^2}{8}\left( \frac{1}{r_k^2} + \frac{2}{r_k} - 2\,r_k - r_k^2\right) , \end{aligned}$$
(21)

and

$$\begin{aligned} y_i(\mathbf {w})\approx & {} y_i[S_0] + \sum _k (D_{ik} w_k + E_{ik} w_k^2), \end{aligned}$$
(22)
$$\begin{aligned} D_{ik}= & {} \frac{1}{2}\bigg [\frac{1}{r_k}\left( y_i[S_k^+] - y_i[S_0]\right) \nonumber \\&- r_k\left( y_i[S_k^-] - y_i[S_0]\right) \bigg ], \end{aligned}$$
(23)
$$\begin{aligned} E_{ik}= & {} \frac{1}{4}\bigg [\left( 1 + \frac{1}{r_k}\right) \left( y_i[S_k^+] - y_i[S_0]\right) \nonumber \\&+ (1 + r_k)\left( y_i[S_k^-] - y_i[S_0]\right) \bigg ]. \end{aligned}$$
(24)

Now, it is a simple numerical task to minimize Eq. (19) with respect to \(\mathbf {w}\). We use MINUIT [24] for the practical applications in the following sections. The found minimum should correspond to that of a full global fit, provided that the approximations (19) and (22) are good enough. This is not trivially true, but we should expect the approximations work better the closer we are to the original minimum. Thus it makes sense to define a “penalty term”

$$\begin{aligned} P = \sum _{k} (A_k (w_k^\text {min})^2 + B_k (w_k^\text {min})^3) \approx \chi ^2_\text {old}(\mathbf {w}^\text {min}) - \chi ^2_0, \nonumber \\ \end{aligned}$$
(25)

which essentially counts how much \(\chi ^2_\text {old}\) has grown from its minimum value, \(w_k^\text {min}\) being the values of \(w_k\) at the minimum of \(\chi ^2_\text {new}(\mathbf {w})\). If \(P \ll \varDelta \chi ^2\), the approximations (19) and (22) should work well and the reweighted results can be viewed as a proxy for those of a full global fit. Once P grows close to or above \(\varDelta \chi ^2\), the results of reweighting become more sensitive on the made assumptions and one should be cautious on the interpretations. Moreover, a large P signals a tension between the original fit and the new data, which might be due to incompatibilities of some data sets, but can also be caused by an inflexible PDF parametrization, or other limitations of theory description, such as missing higher-order corrections.

The beauty of the reweighting method lies in the fact that the reweighted result for any quantity can be obtained simply by using Eq. (22). For example, the new, reweighted, PDFs are obtained by replacing \(y_i\) with \(f_i\). One should note that while this expression is quadratic in \(w_k\), the new PDFs retain a linear dependence on the old ones and thus satisfy the PDF sum rules and DGLAP evolution equations.Footnote 2 This applies also to the new error sets, which can be obtained essentially by following the same procedure as in Sect. 2.1, with the exception that the new Hessian matrix \({\hat{H}}_{kl} = \frac{1}{2} \partial ^2 \chi ^2_\text {new} / \partial w_k \partial w_l|_{\mathbf {w}=\mathbf {w}^\text {min}}\) in

$$\begin{aligned} \chi ^2_\text {new}(\mathbf {w}) \approx \chi ^2_\text {new}(\mathbf {w}^\text {min}) + \sum _{kl}\,(w_k - w_k^\text {min})\,{\hat{H}}_{kl}\,(w_l - w_l^\text {min}) \end{aligned}$$
(26)

can be put to an explicit form

$$\begin{aligned} {\hat{H}}_{kl}= & {} (A_k + 3B_k w_k^\text {min})\delta _{kl}\nonumber \\&+ \sum _{ij}(D_{ik} + 2E_{ik} w_k^\text {min})C^{-1}_{ij}(D_{jl} + 2E_{jl} w_l^\text {min})\nonumber \\&+ \sum _{ij}(2E_{ik}\delta _{kl})C^{-1}_{ij}(y_j(\mathbf {w}^\text {min}) - y_j^\text {data}),\nonumber \\ \end{aligned}$$
(27)

by taking second derivatives of Eq. (19). Diagonalizing \({\hat{H}}\) and finding the deviations in the new eigenvector directions corresponding to \(\varDelta \chi ^2\) growth from \(\chi ^2_\text {new}(\mathbf {w}^\text {min})\), one obtains the parameter values for the new error sets, using which the uncertainties of any quantity can again be obtained according to Eq. (22).

The cubic–quadratic approximation considered above is not applicable to all cases, as it requires the knowledge of the \(\delta z^\pm _k\). Lower-order approximations, initially introduced in Ref. [15], can be obtained from the above results by taking appropriate limits. Taking \(r_k \rightarrow 1\) one finds \(A_k = \varDelta \chi ^2\), \(B_k = 0\), thus recovering the quadratic approximation for \(\chi ^2_\text {old}\). In this limit also the definition of the penalty term in Eq. (25) reduces to that of Ref. [15]. As \(y_i\) retains its quadratic parameter dependence in this limit, we call this a quadratic–quadratic approximation. In many cases this is the best option one can resort to, as it only requires access to the PDF error sets and the value of \(\varDelta \chi ^2\). Even simpler, quadratic–linear, approximation can be achieved by taking also \(E_{ik} \rightarrow 0\). This version is very easy to implement, as finding the new central and error sets in this approximation involves only solving a system of linear equations [15].

2.3 Comment on non-global tolerances

The reweighting method can also be extended to non-global tolerances [17], simply by setting

$$\begin{aligned} A_k= & {} \frac{1}{4}\left( (T_k^+)^2\left( \frac{1}{r_k^2} + \frac{1}{r_k}\right) + (T_k^-)^2\left( r_k + r_k^2\right) \right) , \end{aligned}$$
(28)
$$\begin{aligned} B_k= & {} \frac{1}{8}\left( (T_k^+)^2\left( \frac{1}{r_k^2} + \frac{2}{r_k} + 1\right) - (T_k^-)^2\left( 1 + 2\,r_k + r_k^2\right) \right) ,\nonumber \\ \end{aligned}$$
(29)

where \((T_k^\pm )^2 = \chi ^2_\text {old}(\delta z^\pm _k) - \chi ^2_0\) are the tolerances of the individual error sets, determined by requiring acceptable values of \(\chi ^2\) for each individual data set in the original analysis [17]. While the new, reweighted central PDF set can be obtained uniquely in this way, the determination of the new error sets involves additional arbitrariness. As the new eigenvector directions obtained by diagonalizing the Hessian matrix in Eq. (27) are not parallel to the original ones, it is not directly obvious how large tolerances should be allowed in each of these new parameter directions. It was argued in Ref. [16] that if the new eigendirections are not significantly rotated away from the original ones, it would be sufficient to use the original tolerances \((T_k^\pm )^2\) also for obtaining the new error sets. While this can work in some cases, it would be advisable to have a measure on the amount of parameter rotations in the reweighting to test whether the limits of this assumption are met. Another possibility would be to use a global tolerance for the reweighted PDFs, e.g. by taking the average over the \((T_k^\pm )^2\), but this also would lead to changing the error definition from the original one, thus reducing the comparability of the new and old uncertainties. In general, setting the new non-global tolerances reliably would require a complete refit.

3 CMS 5.02 TeV dijets and their impact on PDFs

The CMS dijet data [4] consist of distributions of dijet pseudorapidity

$$\begin{aligned} \eta _{\mathrm{dijet}} = \frac{1}{2}(\eta ^{\mathrm{leading}} + \eta ^\mathrm{subleading}) \end{aligned}$$
(30)

in bins of average transverse momentum of the jet pair

$$\begin{aligned} p_{\mathrm{T}}^{\mathrm{ave}} = \frac{1}{2}(p_{\mathrm{T}}^{\mathrm{leading}} + p_\mathrm{T}^{\mathrm{subleading}}). \end{aligned}$$
(31)

Here, \(\eta ^{\mathrm{(sub)leading}}\) and \(p_{\mathrm{T}}^{\mathrm{(sub)leading}}\) refer to the pseudorapidity and transverse momentum of the jet with (second to) largest transverse momentum of the event. Jets are defined with the anti-\(k_{\mathrm{T}}\) algorithm [25] using a distance parameter \(R = 0.3\). The events used in the analysis are required to have a leading jet with transverse momentum \(p_{\mathrm{T}}^{\mathrm{leading}} > 30 \text {~GeV}\) and a subleading jet with \(p_{\mathrm{T}}^{\mathrm{subleading}} > 20 \text {~GeV}\) and the two jets are required to have an azimuthal angle separation \(\varDelta \phi > 2\pi /3\). In pPb collisions the two jets are required to be in a rapidity interval \(-3< \eta ^{\mathrm{lab}}_{\mathrm{jet}} < 3\) in the laboratory frame. Due to unequal beam energies, \(E_{\mathrm{p}} = 4 \text {~TeV}\) and \(E_{\mathrm{Pb}} = \frac{82}{208}E_{\mathrm{p}}\), the nucleon–nucleon center-of-mass system is boosted in this frame. To attain corresponding coverages in the center-of-mass frames, CMS measured the pp spectra in the interval \(-3.465< \eta ^\mathrm{lab}_{\mathrm{jet}} < 2.535\). Here, as in the CMS publication, the pp data are shifted in pseudorapidity by \(+0.465\), so that the measured dijets cover a pseudorapidity range \(-3< \eta _{\mathrm{dijet}} < 3\) in both pp and pPb.

The CMS data are self-normalized in each bin of \(p_{\mathrm{T}}^\mathrm{ave}\), i.e. given in the form

$$\begin{aligned} \frac{1}{\mathrm {d}\sigma /\mathrm {d}p_\mathrm {T}^\text {ave}}\,\mathrm {d}^2\sigma /\mathrm {d}p_\mathrm {T}^\text {ave}\mathrm {d}\eta _\text {dijet}. \end{aligned}$$
(32)

This is advantageous due to a partial cancellation of correlated experimental (including luminosity-) uncertainties and theoretical hadronization corrections.Footnote 3 Accordingly, we do not apply nonperturbative corrections to our predictions. We work at NLO as the NNLO calculations of Ref. [27] are not publicly available at this moment. Our theory calculations are performed with NLOJet++ [28] using the anti-\(k_{\mathrm{T}}\) algorithm through FastJet package [29]. We fix the factorization and renormalization scales to be the same, \(\mu _\mathrm{F} = \mu _{\mathrm{R}} = \mu \), and use \(\mu = p_{\mathrm{T}}^{\mathrm{ave}}\) as our central scale choice to keep consistency with the CT14 and EPPS16 fits, but study also variations around this central scale choice to approximate the magnitude of missing higher-order uncertainties (MHOUs).Footnote 4 In all figures, PDF uncertainties are presented with the asymmetric prescription of Eq. (11). As the data correlations are not available, we simply add the statistical and systematical uncertainties in quadrature.

3.1 Proton–proton dijet spectra and CT14 reweighting

Fig. 2
figure 2

Upper panels: distributions of dijets in 5.02 TeV proton–proton collisions against \(\eta _{\mathrm{dijet}}\) and normalized to unity in each bin of \(p_{\mathrm{T}}^{\mathrm{ave}}\). The imposed kinematic cuts are discussed in text. Black markers show the data from the CMS measurement [4] with vertical bars showing the statistical and systematical uncertainties added in quadrature. Solid orange lines represent the results from the NLO pQCD calculation using the central set of the CT14 NLO PDFs [5] with \(\mu = p_{\mathrm{T}}^{\mathrm{ave}}\) scale choice, light orange boxes the associated PDF uncertainties from the CT14 NLO error sets. Lower panels: difference to the central CT14 result. Dashed hollow boxes show the dependence of NLO predictions on factor two upward and downward variations of the scale choice. Dotted lines represent the results from the respective LO pQCD calculation. The results with \(\mu = M_{\mathrm{dijet}}\) scale choice and its factor two variations are indicated in green

Full size image

The self-normalized pp dijet spectra measured by CMS are shown in Fig. 2 along with theory calculations using the CT14 NLO PDFs. While the predictions describe well the \(p_\mathrm{T}^{\mathrm{ave}}\) systematics of the data, we see that the predicted pseudorapidity spectra are systematically wider than the measured distributions, with the discrepancy between the data and CT14 central prediction being much larger than the experimental uncertainties, yielding a very poor figure of merit, \(\chi ^2/N_\mathrm{data} = 7.5\). To study the possible source of this discrepancy, we show in Fig. 2 both the uncertainties from CT14 PDFs, as well as factor of two scale variations around the central scale choice \(\mu = p_{\mathrm{T}}^{\mathrm{ave}}\) and results from a leading order (LO) calculation at the central scale.

We see that in most bins, especially towards high \(p_{\mathrm{T}}^\mathrm{ave}\), the discrepancy between the data and theory is larger than the associated scale uncertainty. As the factor of two scale variations often underestimate the true size of higher order corrections (see e.g. Ref. [27]), not much can be learned from this fact alone. However, as the LO-to-NLO corrections shown in Fig. 2 (lower panels) are of the same size as the scale uncertainties, we should not expect the NLO-to-NNLO corrections to be any larger than these. Hence the discrepancy is unlikely to be just due to missing NNLO terms, which in turn points into the direction that the CT14 PDFs need to be modified for a better description of the data. Towards smaller \(p_{\mathrm{T}}^{\mathrm{ave}}\) the scale variation effects become more important, leaving room for improvement with NNLO corrections. Another possible scale choice would be the invariant mass of the dijet, \(\mu = M_{\mathrm{dijet}}\), a choice which was found in Ref. [27] to yield a better perturbative convergence up to leading-color NNLO precision. We have tested this option, shown also in Fig. 2 (lower panels), and report that here at the NLO level it tends to give smaller scale-uncertainty bands especially at low \(p_{\mathrm{T}}^{\mathrm{ave}}\) and that the results do not differ much from the central \(\mu = p_\mathrm{T}^{\mathrm{ave}}\) predictions. This points again towards smallness of the NNLO corrections. With even slightly wider predictions, \(\mu = M_{\mathrm{dijet}}\) gives a worse data description than the \(\mu = p_\mathrm{T}^{\mathrm{ave}}\) scale choice, and thus we work with the latter in what follows.

To see the modifications on the CT14 PDFs the CMS dijet data would indicate, we have performed a reweighting study with these data. As most of the data points lie outside the CT14 uncertainties, we could expect the needed modifications to be rather strong. Nominally, the CT14 uncertainties correspond to a global tolerance \(\varDelta \chi ^2 = 100\), but to enforce a \(90\%\) confidence level agreement individually with each data set used in the analysis, CT14 uses in addition so called “Tier-2 penalties”. Hence, the parameter variations \(\delta z^\pm _k\) in CT14 do not exactly match with \(\pm \sqrt{100}\), but can be somewhat smaller. As no detailed information is available on how large these deviations are, the best we can do is to assume \(\chi ^2\) to be perfectly quadratic and use \(\varDelta \chi ^2 = 100\). For this reason, we perform the CT14 reweighting in the quadratic–quadratic approximation, noting that the reweighted uncertainties might not be directly comparable with the original ones, and that the new central set underestimates the true impact on CT14, as the use of \(\varDelta \chi ^2 = 100\) overestimates the growth of \(\chi ^2_\text {old}\) in varying the PDF parameters.

Fig. 3
figure 3

The impact of reweighting on CT14 NLO PDFs at \(Q^2 = 10^4~\mathrm{GeV}^2\). The original CT14 PDFs are shown in orange, with the solid line representing the central set PDFs, the ratio to which is shown in each panel. The corresponding PDFs obtained with quadratic–quadratic reweighting using \(\varDelta \chi ^2 = 100\) are shown in red and the central set of the reweighting with \(\varDelta \chi ^2 = 10\) is presented with a solid purple line

Full size image
Fig. 4
figure 4

Upper panels: the impact of reweighting on CT14 predictions of pp dijet spectra. The original predictions are shown in orange and the results obtained with quadratic–quadratic reweighting using \(\varDelta \chi ^2 = 100\) are shown in red. In both cases the solid lines corresponding to the central set and the shaded boxes showing the PDF uncertainty. In addition, resulting spectra from reweighting with \(\varDelta \chi ^2 = 10\) are shown as purple lines. Lower panels show again the difference to the original central CT14 results

Full size image
Fig. 5
figure 5

Comparison of the NLO gluon PDFs of the original and reweighted CT14 sets with those from the MMHT14, NNPDF3.1 and 5-flavour ABMP16 analyses. The uncertainty bands of the latter have been scaled with a factor 1.64 to nominally match with the 90% confidence level definition of the CT14 analysis

Full size image

The resulting reweighted PDFs are compared with the original CT14 NLO PDFs in Fig. 3. For all quark flavours, the found modifications are modest compared to the size of PDF uncertainties. Only at very large x we can see a clear downward bend in the central valence-quark PDFs, caused by the fit trying to adapt to the data at large rapidities, where gluon–valence-quark scattering dominates the cross sections. There is a similar, but even more pronounced, large-x depletion for the gluons. In addition, we find an enhancement for gluons at \(x \sim 0.1\), compensating for the excess in data at midrapidity. Such modifications to gluon PDF are not totally unexpected. The MMHT14 gluon PDF [6], which closely resembles that of CT14, acquires rather similar modifications when confronted with the \(7~\mathrm{TeV}\) high-luminosity inclusive jet data [30]. Also, attributed to including \(8~\mathrm{TeV}\) differential top-quark data, the NNPDF3.1 fit has large-x gluons suppressed compared to CT14 and MMHT14 [31]. In addition, a recent reweighting study using multiple top-production data sets found very similar CT14 modifications as we do here [32]. Thus, we have evidence that the CT14 gluon distribution is simply too hard to be able to fully describe jet and top-quark measurements.

Fig. 6
figure 6

Comparison of CMS 7 TeV inclusive jet measurements [34] and NLO predictions obtained using the CT14 NLO PDFs [5] reweighted with the 5.02 TeV dijet data [4]. The optimal systematic shifts in the correlated experimental uncertainties are applied to the data points (similarly as in Ref. [15]) and only statistical uncertainties are shown. Dashed red lines show the ratio of predictions with the original CT14 PDFs to those with the reweighted PDFs

Full size image

Figure 4 shows the reweighted dijet spectra in comparison to data and original CT14 predictions. The reweighting clearly improves compatibility with the data, especially in the midrapidity region, where the data and theory are now in agreement within the associated uncertainties. At \(\eta _{\mathrm{dijet}} \lesssim -1\), the data still deviates from the reweighted results. This is also reflected in the figure of merit, \(\chi ^2/N_{\mathrm{data}} = 2.0\), which is still quite high, but vastly better than before the reweighting. For a comparison, we have calculated the dijet spectra also using the MMHT14 [6], NNPDF3.1 [31] and 5-flavour ABMP16 [33] NLO PDFs. These yield \(\chi ^2/N_\mathrm{data}\) goodness-of-fit values 4.7, 4.0 and 2.7, respectively, showing that less than perfect agreement with the data is not only a problem with CT14. However, the very strong disagreement between data and CT14 before reweighting appears to be a rather extreme case. In Fig. 5 the gluon PDFs of MMHT14, NNPDF3.1 and ABMP16 are compared with the CT14 before and after the reweighting. The reweighting brings the CT14 gluon distribution to a closer agreement with the other PDFs, particularly at small x to the MMHT14 and NNPDF3.1 and, more importantly, at large x to the NNPDF3.1 and ABMP16. Clearly a reduction in high-x gluons compared to CT14 similar to those in the NNPDF3.1 and ABMP16 fits is preferred by the data.

The penalty term for the reweighted CT14 fit is rather high, with \(P/\varDelta \chi ^2 = 1.17\), clearly indicating that we are reaching the limits of the applicability of the reweighting method. This can be interpreted either as a tension between the dijet data and some datasets used in the CT14 analysis, or as an inflexibility of the CT14 fit form in the high-x region which is probed by the dijets at large rapidities, where the data were not well reproduced and where the data would support even stronger suppression in the PDFs. To test if the CT14 parametrization could adapt to the dijet data, we have performed a reweighting also with an artificially low \(\varDelta \chi ^2 = 10\). In a global fit, this would translate to putting an additional tenfold weight on the new data. The results for the new central PDF set are shown as purple lines in Figs. 3 and 4. With stronger low- and high-x suppression and mid-x enhancement for gluons, this fit achieves a much more reasonable goodness-of-fit \(\chi ^2/N_{\mathrm{data}} = 0.9\) for these data. For this, substantial help from valence quarks, which get strong modifications in this case, is also needed. Still, the data at \(\eta _{\mathrm{dijet}} \lesssim -1\) are not perfectly reproduced, which might be a signal of a parametrization issue, as the relative contribution from the original fit to the total \(\chi ^2\) is decreased with the lowered \(\varDelta \chi ^2\). With \(P/\varDelta \chi ^2 = 3.61\), this fit is in a clear tension with the original CT14 analysis. Of course, once the correlations in the dijet data are made available, one should study whether a shift in some of the systematic parameters could improve the fit at \(\eta _{\mathrm{dijet}} \lesssim -1\). It is also conceivable that the residual disagreement is due to the NNLO corrections.

A comprehensive study of possibly conflicting datasets within CT14 is outside the scope of this article, but as a cross check we have tested the compatibility of the reweighted PDFs with the CMS 7 TeV inclusive jet measurements [34] which are included in the CT14 analysis. For these calculations we use the pre-computed fastNLO grids [35], setting the renormalization and factorization scales equal to the transverse momentum \(p_{\mathrm{T}}\) of the individual jet as in the CT14 analysis. Figure 6 shows the data-to-theory ratio for the NLO predictions with the CT14 PDFs reweighted with the dijet data using \(\varDelta \chi ^2 = 100\). Also the ratios of the original CT14 central predictions with the reweighted ones are indicated. The data-to-theory agreement happens to be even slightly better for the reweighted PDFs, with \(\chi ^2/N_{\mathrm{data}} = 1.2\), than for the original set, for which \(\chi ^2/N_{\mathrm{data}} = 1.3\). Thus we find that, in the light of reweighting, the CMS measurements of inclusive jets at 7 TeV and dijets at 5.02 TeV are mutually compatible.

3.2 Significance of proton PDF uncertainties in proton–lead dijet spectra

The pPb dijet spectra, shown in Fig. 7, have a rather similar data-to-theory systematics as we had in the pp case. Here, we use the EPPS16 nuclear modifications along with the CT14 NLO proton PDFs in the predictions, i.e. the PDF of a flavour i in a proton bound in lead at scale \(Q^2\) is obtained with

$$\begin{aligned} f^{\mathrm{p/Pb}}_i(x,Q^2) = R^{\mathrm{Pb}}_i(x,Q^2) f^{\mathrm{p}}_i(x,Q^2), \end{aligned}$$
(33)

where \(R^{\mathrm{Pb}}_i\) is the nuclear modification from the EPPS16 analysis and \(f^{\mathrm{p}}_i\) the corresponding CT14 PDF of the free proton. The total PDF uncertainties in the cross sections are calculated with

$$\begin{aligned} \delta X^\pm _\text {total} = \sqrt{\left( \delta X^\pm _\text {EPPS16}\right) ^2 + \left( \delta X^\pm _\text {CT14}\right) ^2}, \end{aligned}$$
(34)

where \(\delta X^\pm _\text {EPPS16}\) are the upward and downward uncertainties obtained with Eq. (11) using the EPPS16 error sets and keeping the CT14 central set fixed, and \(\delta X^\pm _\text {CT14}\), respectively, the uncertainties from the CT14 error sets keeping the EPPS16 central set fixed.

Again, these predictions give wider distributions than seen in the CMS data, resulting with \(\chi ^2/N_{\mathrm{data}} = 6.9\). While in this case the data points are mostly within the combined nuclear and free-proton PDF uncertainty bands, we can expect that the modifications to the CT14 PDFs, which were found necessary to improve the description of the pp data, play a role also here. Indeed, in Fig. 8 we show results with the PDFs obtained by reweighting CT14 with the pp data, observing a clear improvement in the data to theory agreement. We obtain \(\chi ^2/N_{\mathrm{data}} = 2.8\) for the predictions with CT14 reweighted using \(\varDelta \chi ^2 = 100\) and \(\chi ^2/N_{\mathrm{data}} = 1.6\) when using \(\varDelta \chi ^2 = 10\). These numbers are somewhat higher than what we obtained in the pp case, reflecting the fact that also the EPPS16 nuclear modifications need to be adjusted for optimal description of the data. This can also be seen by comparing the data-to-theory agreement in pPb at \(\eta _{\mathrm{dijet}} \gtrsim 2\) to that in pp: while the CT14 predictions reweighted using \(\varDelta \chi ^2 = 100\) describe well the pp data in these rapidities, the pPb data points lie systematically below the predictions, which hints a preference for deeper nuclear shadowing – the suppression in the gluon PDF, \(R^{\mathrm{Pb}}_g < 1\), at small x – than that in the EPPS16 central set. We will verify this claim in the next section.

Fig. 7
figure 7

As Fig. 2, but now with pPb data and predictions with EPPS16 nuclear modifications imposed on the CT14 NLO proton PDFs and omitting the results with \(\mu = M_\mathrm{dijet}\) for clarity. Light blue boxes show the combined uncertainty from the CT14 and EPPS16 PDFs

Full size image
Fig. 8
figure 8

As Fig. 4, but now with pPb data and with EPPS16 nuclear modifications imposed on the original and reweighted CT14 PDFs. Only uncertainties from the free-proton PDFs are shown

Full size image

An important thing to notice here is that most of the deviations from central theory predictions actually originate from the issues with the free-proton PDFs instead of the nuclear modifications. This large free-proton PDF bias prevents a clean extraction of the PDF nuclear modifications from the pPb spectra. The dijet spectra are certainly not the only pPb observable sensitive to such a free-proton PDF dependence, but the refined proton PDFs found here could also have an effect for example on the predictions for inclusive \(t{\bar{t}}\) production at 8.16 TeV pPb collisions where calculations with CT14+EPPS16 overshoot, but are still compatible with the data [36].

3.3 Nuclear modification ratio and EPPS16 reweighting

Fig. 9
figure 9

The nuclear modification ratio of normalized pPb and pp differential cross sections. Black markers show the data from CMS measurement [4] with vertical bars showing the statistical and systematical uncertainties added in quadrature. Solid orange lines represent the NLO pQCD calculation with \(\mu = p_{\mathrm{T}}^{\mathrm{ave}}\) scale choice using the central set of the CT14 NLO PDFs [5] with EPPS16 [18] nuclear modifications

Full size image
Fig. 10
figure 10

The impact of reweighting on EPPS16 predictions of the nuclear modification ratio of the dijet spectra. The original predictions are shown with solid blue lines and light blue boxes representing the central predictions and the nPDF uncertainties, respectively. The corresponding results after the reweighting are shown with solid black lines and purple boxes

Full size image

Let us now consider the nuclear modification ratio of the normalized dijet spectra discussed above, defined as

$$\begin{aligned} R_{\mathrm{pPb}}^{\mathrm{norm.}} = \frac{\frac{1}{\mathrm {d}\sigma ^\mathrm{pPb}/\mathrm {d}p_\mathrm {T}^{\mathrm{ave}}}\,\mathrm {d}^2\sigma ^\mathrm{pPb}/\mathrm {d}p_\mathrm {T}^{\mathrm{ave}}\mathrm {d}\eta _\mathrm{dijet}}{\frac{1}{\mathrm {d}\sigma ^\mathrm{pp}/\mathrm {d}p_\mathrm {T}^{\mathrm{ave}}}\,\mathrm {d}^2\sigma ^\mathrm{pp}/\mathrm {d}p_\mathrm {T}^{\mathrm{ave}}\mathrm {d}\eta _{\mathrm{dijet}}}. \end{aligned}$$
(35)

As we have seen that the dijet rapidity distributions in pp and pPb have very similar dependence on the free proton PDFs, we can expect this dependence to efficiently cancel in the ratio. This statement is verified in Fig. 9, where we observe the uncertainty band given by CT14 PDFs to be vanishingly small. Also the scale uncertainties, while being larger than the CT14 uncertainties, are small in this observable, implying that MHOUs can be expected to be small as well. This leaves the nuclear modifications as the dominant source of theory uncertainty.

Fig. 11
figure 11

The impact of reweighting the EPPS16 nPDFs with the data on the nuclear modification ratio of the dijet spectra. The original and reweighted EPPS16 nuclear modifications for the lead nucleus are presented at the parametrization scale \(Q^2 = 1.69~\mathrm{GeV}^2\). For better visibility, the s-quark modifications are presented with a different vertical axis scaling

Full size image
Fig. 12
figure 12

The EPPS16 gluon nuclear modifications in Pb at the scales \(Q^2 = 10~\mathrm{GeV}^2\) and \(Q^2 = 10^4~\mathrm{GeV}^2\) before and after reweighting with the dijet data

Full size image

We observe that the CMS data and EPPS16 predictions are in good agreement within the uncertainties. This does not come as a surprise, as part of these data, namely the high-\(p_{\mathrm{T}}^\mathrm{ave}\) part of the pPb cross section [37], were used in the EPPS16 fit. Still, this agreement is not trivial as with the new pp baseline and being a more differential measurement, these \(R_{\mathrm{pPb}}^{\mathrm{norm.}}\) data contain plenty of new information compared to the 7 data points of forward-to-backward ratios included in the EPPS16 analysis. As was anticipated above, the data points at forward rapidities deviate from the central EPPS16 prediction, indicating a preference for a deeper shadowing in the nPDFs.

Compared to the data, the EPPS16 predictions have much larger uncertainties, which promises a good constraining power when fitting to these data. To study the impact these data would have had in the EPPS16 fit, we have performed a reweighting in the cubic–quadratic approximation introduced in Sect. 2.2, using \(\varDelta \chi ^2 = 52\) and taking the values of \(\delta z^\pm _k\) from Table 2 of Ref. [18]. The results for \(R_\mathrm{pPb}^{\mathrm{norm.}}\) are shown in Fig. 10. Most notably, there is a vast reduction in the EPPS16 uncertainties. Also, at forward rapidities the central prediction comes down a bit, as is expected from the low-lying data points in this region. In the backward direction a slight enhancement in the central prediction can be observed, but this is far less prominent than the suppression in the forward bins. In total, we obtain an improvement in the goodness of fit from \(\chi ^2/N_{\mathrm{data}} = 1.7\) to 1.4 with a penalty \(P/\varDelta \chi ^2 = 0.14\).

The corresponding effects on the EPPS16 nuclear modifications in lead at the parametrization scale \(Q^2 = 1.69~\mathrm{GeV}^2\) are presented in Fig. 11. There is a striking impact on gluon modification uncertainties, which are reduced across all x. In the best-constrained mid-x region, the uncertainties are reduced to less than half of their original size. As the uncertainty band lies clearly above unity in this region, we find strong evidence for gluon antishadowing in lead. At small x, the reweighted uncertainty band goes respectively below unity, giving evidence for gluon shadowing. These findings are in accordance with those of Ref. [38], where inclusive heavy-flavour production data from measurements at the LHC were used to study the gluon PDF modifications in nuclei. As expected from inspecting the ratio of the dijet spectra, the new central set seems to support stronger shadowing than in the original EPPS16 central fit.

Even with the increased gluon shadowing, the most forward bins of \(R_{\mathrm{pPb}}^{\mathrm{norm.}}\) are not well reproduced by the reweighted results, which is also the reason why the \(\chi ^2/N_{\mathrm{data}}\) remained somewhat high even after the reweighting. To be consistent with these forward data points, a very deep shadowing for the gluons would be required. Moreover, the probed x region changes very little between the last and second-to-last \(\eta _{\mathrm{dijet}}\) data point, and thus such a steep drop as that suggested by the data is difficult to attain. This is because the DGLAP evolution efficiently smooths out even steep structures in the gluon nuclear modification, as can be seen in Fig. 12 where we show the gluon nuclear modifications evolved to higher scales. We also note that the systematic uncertainty dominates in the last \(\eta _\mathrm{dijet}\) bins, and thus taking into account the data correlations, once available, could improve the fit quality. These findings should, in the future, be contrasted also with the recent ATLAS conditional yield measurement, where an order of 10–\(20\%\) nuclear suppression for dijets was found in the most forward configuration [39].

Also at large x, the reweighted gluon modifications are better constrained than in the original EPPS16 analysis. The new central set has \(R_g^{\mathrm{Pb}}\) closer to unity at x around 0.7. This is partly enforced by momentum sum rule in combination with the stiffness of the EPPS16 fit function and the deepened small-x shadowing. In any case, the uncertainty remains large, and either an enhancement or a suppression for gluons is possible in this region. On this basis, the conclusion made in Ref. [4], that the dijet data would give evidence of large-x gluon suppression, seems premature. This claim was based on comparison of the data with EPS09 [40] and DSSZ [41] nPDFs, where the former, with gluon suppression at large x, agreed well with the data at backward rapidities, but the latter, having the nuclear gluons unmodified, did not. However, going towards backward rapidities, and thus larger x from the Pb side, the contribution of nuclear quarks to the dijet cross section grows rapidly. Hence the difference in predictions with EPS09 and DSSZ in this region has a large contribution from different valence quark modifications. As DSSZ has much smaller large-x suppression for valence quarks than EPS09 (see e.g. Ref. [42]), this also partly explains the difference in the dijet predictions of Ref. [4].

Fig. 13
figure 13

The impact on the average valence and sea quark and gluon modifications under different approximations in the reweighting

Full size image

On these grounds, it might appear surprising that the dijet data are not able to constrain the valence quark modifications at all, as can be seen from the first two panels in Fig. 11. The reason for this is that due to smallness of isospin corrections [43], the backward dijet data mainly probe the average valence modifications,

$$\begin{aligned} R_{u_{\mathrm{V}} + d_{\mathrm{V}}}^{\mathrm{Pb}} = \frac{u^{\mathrm{p/Pb}}_\mathrm{V}+d^{\mathrm{p/Pb}}_{\mathrm{V}}}{u^{\mathrm{p}}_{\mathrm{V}}+d^{\mathrm{p}}_{\mathrm{V}}}, \end{aligned}$$
(36)

shown in Fig. 13. This combination is much better constrained than the individual flavours shown in Fig. 11 and has vastly smaller uncertainties at large x than the gluon modifications. Thus, while large-x valence quarks dominate the dijet cross section at backward rapidities, the uncertainty in the EPPS16 predictions in this region comes dominantly from the less-constrained gluons, and hence it is the gluon modifications which are constrained in the reweighting. Figure 13 shows also the average sea quark modification

$$\begin{aligned} R_{{\overline{u}}+{\overline{d}}+{\overline{s}}}^{\mathrm{Pb}} = \frac{{\overline{u}}^{\mathrm{p/Pb}}+{\overline{d}}^\mathrm{p/Pb}+{\overline{s}}^{\mathrm{p/Pb}}}{{\overline{u}}^\mathrm{p}+{\overline{d}}^{\mathrm{p}}+{\overline{s}}^{\mathrm{p}}}, \end{aligned}$$
(37)

which is the dominant quark combination constrained at forward rapidities. We observe a modest reduction in the small-x uncertainty, much smaller than that for the gluons. At the level of individual flavours, shown in Fig. 11, these constraints affect mostly the s-quark modifications, which were poorly constrained in EPPS16.

3.4 Importance of non-quadratic and non-linear terms in reweighting

We may now ask whether the inclusion of higher-order (non-quadratic and non-linear) components in the reweighting had a sizable effect on our results. Figure 13 shows the impact of the dijet data on the EPPS16 nuclear modifications in all three approximations discussed in Sect. 2.2. While, for simplicity of presentation, we show only the average valence and light-sea-quark modifications in addition to those for gluons, the conclusions below apply to individual flavours as well. We find that the cubic–quadratic and quadratic–quadratic approximations give almost identical results. This is rather easy to understand: The new data are precise enough to dominate the shape of the total \(\chi ^2\) function in the parameter directions that it constrains (mainly those related to gluon degrees of freedom), making the non-quadratic components sub-dominant in the reweighting. Moreover, as the new central set does not divert far from the original, we are working in a region where the quadratic approximation for \(\chi ^2_\text {old}\) is rather good. Under different circumstances this might not be the case and the cubic terms could alter the reweighting results significantly.

Next, we consider the reweighting results in the quadratic–linear approximation. Here, we use the linear approximation for the cross sections, but decide to keep the quadratic dependence in the PDFs for better comparability.Footnote 5 Again, the differences to the results of the cubic–quadratic approximation are rather modest, though for the high-x gluons the quadratic–linear approximation seems to suggest slightly less stringent constraints. The similarity of results in the different approximations can also be seen as a reassuring fact: the results of reweighting do not seem to depend on minute details of our method and we seem to be able to make reliable conclusions based on rather limited information about the original global analysis, at least in this particular case. The obtained results are thus not likely to change if even higher-order contributions are added.

4 Summary and conclusions

In this work, we have presented a non-quadratic extension of the Hessian PDF reweighting introduced in Ref. [15] and applied the method in the context of CMS dijet measurements at 5.02 TeV. This improved method makes use of the knowledge of parameter variations at which the error sets of the original PDFs are defined, to solve for cubic components of the \(\chi ^2\) function before inclusion of new data. Similarly, quadratic components in the responses of observables to parameter variations were taken into account. The additional information needed in this cubic–quadratic approximation prevented us from using it when reweighting the CT14 NLO PDFs with the pp dijet distributions, where we had to resort to a simpler quadratic–quadratic approximation, but we were able to apply it to reweight the EPPS16 nPDFs, for which the needed information is available, with the nuclear modification ratio of the dijet spectra. While no large differences were found in the results of reweighting EPPS16 in the cubic–quadratic or quadratic–quadratic approximation, this observation was limited to one specific case, and under different circumstances the cubic terms could become more important. We thus encourage PDF fitters to publish the details of their analysis to a sufficient accuracy, such that the reweighting including the higher-order terms becomes possible. This can be done by publishing the numerical values of the \(\delta z^\pm _k\) parameters as defined in Sect. 2 in addition to the tolerance \(\varDelta \chi ^2\). Care must be taken in communicating which error set corresponds to each of these values, so that there is no chance of misinterpretation e.g. in what is called a “plus” and what a “minus” direction. A neat way to do this with LHAPDF [44] would be to set in each PDF grid file a custom flag such as “ParamVal” to hold the value \(\delta z^\pm _k\). These parameter values could then be retrieved by using the method info().get_entry("ParamVal") for each of the PDF error sets.

Comparing the measured pp dijet pseudorapidity spectra with theory calculations using the CT14 NLO PDFs revealed a large discrepancy. We showed that at high \(p_{\mathrm{T}}^{\mathrm{ave}}\) this difference is larger than the associated scale uncertainties and exceeds the size of the NLO corrections, thus being unlikely due to missing NNLO terms alone. This suggested the need for modifying the CT14 PDFs to reach a better agreement with the data. In reweighting CT14 with the dijet data, the gluon PDF acquired significant modifications, especially at large x, where a substantial reduction was observed. We discussed also evidence from other studies pointing into the same direction. After reweighting, a much more reasonable \(\chi ^2\) value for the dijet data was found, but this came with a price of a rather high penalty term, i.e. the new central set had diverted quite far from the original minimum. The reason for this apparent discrepancy between CT14 and the dijet data remains elusive. We tested the reweighted PDFs against CMS 7 TeV inclusive jet measurements finding good agreement, and thus no conflict between the considered dijet and inclusive jet data. By performing a reweighting with an artificially low \(\varDelta \chi ^2\), we showed that the CT14 PDFs still had trouble in reproducing the data at \(\eta _{\mathrm{dijet}} \lesssim -1\), signaling a possible parametrization issue, although NNLO corrections and correlated systematics can also play a role here. Solving this issue is beyond the reach of the reweighting tools and should be studied in the context of a global analysis.

Similar discrepancy as seen with the pp spectra is observed also in the case of pPb. We showed that applying the same CT14 modifications as found in the reweighting with pp data substantially improves the data-to-theory agreement also in pPb. As the pPb dijet distributions contain a substantial free-proton PDF dependence, a clean extraction of their nuclear modifications is not possible from these data directly. Taking the ratio of the pPb and pp spectra, however, leads to a very efficient cancellation of not only the free-proton uncertainties but also of the scale uncertainties, thus giving an excellent probe of the nPDFs. We showed that the measured nuclear-modification ratio of dijet spectra is in a good agreement with the NLO predictions using the EPPS16 nPDFs. Some deviation from the EPPS16 central prediction was observed at \(\eta _{\mathrm{dijet}} \gtrsim 2\), supporting a stronger shadowing for gluons than present in the EPPS16 central set. As a whole, these data give compelling evidence of small-x gluon nuclear shadowing and mid-x antishadowing, as was revealed in reweighting the EPPS16 nPDFs. We obtained significant new constraints on the EPPS16 gluon modifications in lead throughout the probed range, reducing the uncertainties even to less than half of their original size.