PDF reweighting in the Hessian matrix approach

JHEP12(2014)100 Published for SISSA by Springer Received:May 5, 2014 Revised:November 11, 2014 Accepted:November 21, 2014 Published:December 15, 2014 PDF reweighting in the Hessian matrix approach Hannu Paukkunena,b and Pia Zuritac aDepartment of Physics, University of Jyv¨askyl¨a, P.O. Box 35, FI-40014 University of Jyv¨askyl¨a, Finland bHelsinki Institute of Physics, P.O. Box 64, FI-00014 University of Helsinki, Finland cDepartamento de F´ısica de Part´ıculas and IGFAE, Universidade de Santiago de Compostela, E-15782 Galicia, Spain E-mail: [email protected],[email protected] Abstract: We introduce the Hessian reweighting of parton distribution functions (PDFs). Similarly to the better-known Bayesian methods, its purpose is to address the compatibility of new data and the quantitative modifications they induce within an existing set of PDFs. By construction, the method discussed here applies to the PDF fits that carried out a Hessian error analysis using a non-zero tolerance ∆χ2. The principle is validated by considering a simple, transparent example. We are also able to establish an agreement with the Bayesian technique provided that the tolerance criterion is appropriately accounted for and that a purely exponential Bayesian likelihood is assumed. As a practical example, we discuss the inclusive jet production at the LHC. Keywords: Statistical Methods, Parton Model, QCD ArXiv ePrint: 1402.6623 Open Access,c The Authors. Article funded by SCOAP3.doi:10.1007/JHEP12(2014)100 JHEP12(2014)100 Contents 1 Introduction 1 2 The Hessian method 2 3 The Hessian reweighting 4 4 Bayesian methods 6 5 Simple example 9 6 Equivalence of the Bayesian and Hessian reweighting 14 7 Inclusive jet production at the LHC 15 8 Conclusion 23 1 Introduction A large part of the present high-energy collider physics depends, in one way or another, on the knowledge of the parton distribution functions (PDFs). The use of PDFs leans on a cornerstone theorem of Quantum Chromo Dynamics (QCD), the collinear factorization [1, 2]. Although it can be formally proven only for the most simple cases, it is often assumed to work in general. Ultimately, it is the agreement with the experimental data that decides whether such an assumption is valid. The PDFs are traditionally determined in global analyses [3–6] finding the parametrization that can optimally reproduce a variety of experimental data. This is a complex procedure requiring the ability to efficiently solve the parton evolution equations and to calculate higher-order QCD cross-sections. Nontrivial issues are also the actual way of finding the best fit and quantifying its uncertainties. Although various PDF parametrizations are publicly available for a general user, for a long period of time it was difficult for e.g. an experimental collaboration to understand what would be the implications of their measurements in the context of a global PDF fit. For example, although a given measurement would be known to be most sensitive to, say, the up quark distribution, some other data in a global fit may already provide the up quarks with more stringent constraints, and the real advantage of the measurement would be due to a sub-leading contribution, say, from the strange quarks. The Bayesian reweighting technique, first introduced in [7] and later on elaborated by the Neural Network PDF (NNPDF) collaboration [8,9], provides a way of addressing the consistency and quantitative effects of a new experimental evidence in terms of PDF fits. In – 1 – JHEP12(2014)100 essence, the underlying probability distribution, represented in the NNPDF philosophy [6] by an ensemble (∼1000) of PDF replicas, is updated by assigning each replica a certain weight based on the new data. This method has become an increasingly popular way to estimate the effects of e.g. new LHC measurements [10–16]. The drawback is that it has been proven to work only in conjunction with the NNPDF fits while the majority of the existing PDF fits use a rather different way of quantifying the PDFs and their uncertainties. Along with the best fit found by χ2minimization, they provide a collection (∼50) of Hessian error sets [17] that quantify the neighborhood of the central fit within a certain confidence criterion ∆χ2. An extension of the Bayesian reweighting technique to this particular case was suggested in [18], and has thereafter been used in some occasions [11,19, 20]. However, a recent study [21] revealed clear deviations when comparing the results from reweighting to the ones obtained by a direct fit, indicating that the proposed generalization is not accurate. Here we take a different strategy. Based on the ideas presented in ref. [22], our principal goal is to show how a general user can directly study the consistency and consequences of a new data set within an existing set of PDFs that comes with the Hessian error sets without having to lean on the Bayesian techniques. The method naturally incorporates the confidence criterion ∆χ2defined in the original fit and, by considering a simple numerical example, we argue that it is perfectly compatible with a new fit. Our second objective is to understand how the procedure suggested here relates to the Bayesian reweighting and to the discrepancies found in [21]. We prove that the original Bayesian method, proposed in [7] and advocated recently in [23], is equivalent with the one introduced here once the ∆χ2criterion is properly incorporated. 2 The Hessian method The usual definition of an optimum correspondence between data and a set of PDFs f≡ f(x, Q2) that depends on certain fit parameters {a}, is the minimum of a χ2-function. In its most simple form, we can write it as χ2{a}=X k"Xtheory k[f]−Xdata k δdata k#2 ,(2.1) where Xtheory k[f] are the theory predictions depending on the PDFs. The corresponding experimental values are denoted by Xdata kand their uncertainty by δdata k. Modifications to this definition are necessary if the experimental errors are correlated or some data sets are emphasized in the fit by assigning them an additional weight. In the Hessian approach to quantify the PDF errors [17], the behaviour of χ2around the best fit S0is approximated by a second order polynomial in the space of fit parameters {a} χ2{a} ≈ χ2 0+X ij δaiHijδaj,(2.2) where δaj≡aj−a0 jare the excursions from the best-fit values and χ2 0is the minimum value of χ2. Being symmetric, the Hessian matrix Hij has Neig orthonormal eigenvectors – 2 – JHEP12(2014)100 v(k)and eigenvalues ksatisfying Hijv(k) j=kv(k) i,(2.3) X j v(k) jv(`) j=X j v(j) kv(j) `=δk`.(2.4) Defining a new set of variables as zk≡√kX j v(k) jδaj,(2.5) one easily finds that χ2{a} ≈ χ2 0+X i z2 i.(2.6) That is, the transformation in eq. (2.5) diagonalizes the Hessian matrix. A criterion is needed to specify how much the term Piz2 ican grow while the corresponding PDFs still remain “acceptable”. Those PDF fits that employ the ideal choice ∆χ2= 1 are usually limited a smaller set of data [24,25], while the global fits prefer to take ∆χ2>1 [4, 5] to account for small inconsistencies among the data sets and to compensate for the parametrization bias [28]. It follows [17] that the corresponding uncertainty for a PDF- dependent quantity O=O[f] can be computed as (∆O)2= ∆χ2X k∂O ∂zk2 .(2.7) An essential feature of the Hessian approach is the introduction of the PDF error sets S± k, defined customarily (along with the best fit S0) in the z-space as z(S0) = (0,0,...,0) , z(S± 1) = ±p∆χ2(1,0,...,0) ,(2.8) z(S± 2) = ±p∆χ2(0,1,...,0) , . . . z(S± Neig ) = ±p∆χ2(0,0,...,1) . Using these sets, one can evaluate the derivatives in eq. (2.7) by a linear approximation ∂O ∂zk≈OS+ k−OS− k 2p∆χ2,(2.9) such that (∆O)2=1 4X kOS+ k−OS− k2.(2.10) This formula, or a generalization for asymmetric errors, provides an extremely simple and useful recipe for propagating the PDF-uncertainties to observables. Recently, it has become fashionable [4,5] to abandon the fixed ∆χ2tolerance and define the PDF error sets instead by a “dynamic tolerance” zi(S± k)≡ ±t± kδik,(2.11) which coincides with eq. (2.8) if the tolerance parameters t± kare equal to p∆χ2. – 3 – JHEP12(2014)100 3 The Hessian reweighting Let us now consider a new set of data ~y =y1, y2, . . . , yNdata with covariance matrix C. Our goal here is twofold: to find out whether these new data are consistent within the original set of PDFs and, if so, what would be the effect of incorporating them into the original analysis. In order to answer these questions, we consider a function χ2 new defined as χ2 new ≡χ2 0+ Neig X k z2 k+ Ndata X i,j=1 (yi[f]−yi)C−1 ij (yj[f]−yj),(3.1) where we have simply added the contribution of the new data on top of the “old” χ2in eq. (2.6). Using a similar linear approximation as earlier, we can estimate the theoretical values yi[f] in arbitrary z-space coordinates by yi[f]≈yi[S0] + Neig X k=1 ∂yi[S] ∂zkS=S0 zk≈yi[S0] + Neig X k=1 Dikwk,(3.2) where we have defined Dik ≡yiS+ k−yiS− k 2(3.3) wk≡zk 1 2t+ k+t− k.(3.4) Thus, χ2 new is a continuous, quadratic function of the parameters wk, and its minimum is given simply by ~ wmin =−B−1~ a,(3.5) where the matrix Band vector ~ aare Bkn =X i,j DikC−1 ij Djn +t+ k+t− k 22 δkn ,(3.6) ak=X i,j DikC−1 ij (yj[S0]−yj).(3.7) An important feature of the solution is the “penalty term” P≡ Neig X k=1 t+ k+t− k 2wmin k2t± k→√∆χ2 −−−−−−−→ ∆χ2 Neig X k=1 (wmin k)2,(3.8) which we can use to decide whether the new data set is consistent within the original PDFs. If an overall tolerance ∆χ2was defined in the original fit, having P∆χ2means that the new data could have been incorporated into the original fit without causing a conflict with the other data. On the other hand, if P&∆χ2the new data appears to display significant tension with the considered set of PDFs. In the case that the original fit used a dynamic tolerance a simple interpretation like this is lost, and one can only check whether the new – 4 – JHEP12(2014)100 z-space coordinates remain within the tolerance parameters. However, even a large penalty term does not necessarily mean that the new data would be incompatible with the other data. A situation like this may arise if the new data probe unconstrained components of PDFs whose behaviour was fixed by hand. For example, some recent PDF fits [26] still assume s(x)∝(u(x) + d(x)) for the strange quark distribution and confronted with data sensitive to the strange quarks could lead to this kind of situation. The components of the weight vector ~ wmin also specify the set of PDFs fnew that corresponds to the new global minimum. They can be easily calculated by taking yi= f(x, Q2) in eq. (3.2). That is, fnew ≈fS0+ Neig X k=1 fS+ k−fS− k 2!wmin k.(3.9) The resulting new PDFs are linear combinations of the original ones — they have been “reweighted”. We note that the new PDFs constructed in this way still satisfy the necessary sum rules. For instance, as the original best fit S0and the error sets S± ksatisfy the momentum sum rule Z1 0 dxx X f fS0=Z1 0 dxx X f fS± k= 1,(3.10) then Z1 0 dxx X f fnew =Z1 0 dxx X f fS0+X k wmin k 2 Z1 0 dxx X f fS+ k−Z1 0 dxx X f fS− k  = 1 + X k wmin k 2[1 −1] = 1. Due to the linearity of the parton evolution equations, also fnew satisfies them. Thus, the reweighted distributions comprise a proper set of PDFs which can be consistently utilized in perturbative QCD calculations. One can also construct the new PDF error sets. Indeed, eq. (3.1) can be rewritten as χ2 new =χ2 new~ w=~ wmin +X ij δwiBijδwj,(3.11) where ~ δw=~ w−~ wmin, and the matrix Btakes the role of the Hessian matrix (compare to eq. (2.2)). This can be brought into a diagonal form by an analogue of the transformation in eq. (2.5), and the new error sets defined as δwi(ˆ S± k) = ±ˆv(k) ir1 ˆk ˆ t± k,(3.12) where ˆv(k) iare the eigenvectors and ˆkthe eigenvalues of the matrix B. The original overall tolerance can be set easily by ˆ t± k=p∆χ2. – 5 – JHEP12(2014)100 Non-linear extension of the Hessian reweighting. The linear approximation of eq. (3.2) can be improved by including also quadratic terms of wkas yi[S] = yi[S0] + Neig X k=1 1 2"yiS+ k−yi[S0] t+ k/t− k−yiS− k−yi[S0] t− k/t+ k#wk(3.13) + Neig X k=1 t+ k+t− k 4"yiS+ k−yi[S0] t+ k +yiS− k−yi[S0] t− k#w2 k correcting for the possible non-linear behaviour. Using this in eq. (3.1), χ2 new becomes a quartic function of wkwhose minimum must be found by numerical methods. The corresponding PDFs can be computed by taking yi=f(x, Q2) in eq. (3.13). The matrix Bin eq. (3.11) gets replaced by Bnon.lin kn =X i,j ∂yi[f] ∂wkC−1 ij ∂yj[f] ∂wn(3.14) +X i,j ∂2yi[f] ∂wk∂wnC−1 ij (yj[f]−yj) + t+ k+t− k 22 δkn, where the partial derivatives read ∂yi[f] ∂wk=1 2"yiS+ k−yi[S0] t+ k/t− k−yiS− k−yi[S0] t− k/t+ k#(3.15) +t+ k+t− k 2"yiS+ k−yi[S0] t+ k +yiS− k−yi[S0] t− k#wk ∂2yi[f] ∂wk∂wn=t+ k+t− k 2"yiS+ k−yi[S0] t+ k +yiS− k−yi[S0] t− k#δkn,(3.16) and are understood to be evaluated at the found minimum. 4 Bayesian methods Given a large ensemble of PDFs fk,k= 1 . . . Nrep, such as those of the NNPDF collaboration [6], that represents the underlying probability distribution Pold(f) of the PDFs, one can compute the expectation value hOi and variance δhOi for an observable Oas hOi =1 Nrep Nrep X k=1 O[fk],(4.1) δhOi =v u u t1 Nrep Nrep X k=1 (O[fk]−hOi)2.(4.2) Using the laws of statistics, the initial probability distribution Pold(f) can be updated to include also additional information contained in a new set of data ~y, since, by the Bayes theorem, Pnew(f)∝P(~y|f)Pold(f),(4.3) – 6 – JHEP12(2014)100 where P(~y|f) stands for the conditional probability (the likelihood function) for the new data, given a set of PDFs. It follows that the average value for any observable depending on the PDFs becomes a weighted average hOinew =1 Nrep Nrep X k=1 ωkO[fk],(4.4) δhOinew =v u u t1 Nrep Nrep X k=1 ωk(O[fk]−hOinew)2,(4.5) where the weights ωkturn out to be proportional to the likelihood function P(~y|f). The question of how to choose the likelihood appropriately has been recently revisited in [23] but the conclusive answer, if it exists, remains to be given. Two options, corresponding to different choices of the likelihood, have been discussed in the literature. The one suggested originally by Giele and Keller (GK) [7] follows from taking P(~y|f)dnyas the probability to find the new data to be confined in a differential element dnyaround ~y resulting in ωGK k=exp −χ2 k/2 (1/Nrep)PNrep k=1 exp −χ2 k/2,(4.6) where χ2 k= Ndata X i,j=1 (yi[fk]−yi)C−1 ij (yj[fk]−yj).(4.7) The option advocated by the NNPDF collaboration derives from taking P(~y|f)dχ as the probability for the corresponding χ≡pχ2to be confined in a differential volume dχ around χ, giving instead1 ωchi−squared k=χ2 k(Ndata−1)/2exp −χ2 k/2 (1/Nrep)PNrep k=1 χ2 k(Ndata−1)/2exp −χ2 k/2,(4.8) which has been shown to be consistent with a direct fit in the NNPDF framework [8,9]. It was pointed out in ref. [23] that the former weights contain more information on the new data than the latter ones, as a given data set uniquely determines the value of χ2, while a fixed χ2may correspond to various different data sets. The generic behaviour of these weights with respect to χ2per number of points for Ndata = 10 is shown in figure 1. While the GK weights are always higher for those replicas that give lower χ2, the NNPDF option obviously favors ones with χ2/Ndata ≈1. We note that the latter likelihood may lead to a following situation: if the value of χ2/Ndata (computed with the expectation values of the observables or as an average of the individual χ2 ks) is less than unity before the reweighting, it can happen that the reweighting actually causes χ2for the new data to grow since the replicas with χ2/Ndata ≈1 are favored. However, if new data are directly included into a PDF fit as in eq. (3.1), the value of χ2for the new data should only decrease.2 1We dub these weights ωchi−squared ksince their behaviour is very close to the usual χ2distribution. 2In PDF fits with extremely flexible fit functions, there is a danger that the PDFs that correspond to the minimum χ2unwantedly reproduce also random fluctuations of the data. Here, we assume that this is – 7 – JHEP12(2014)100 Figure 1. Comparison of the likelihoods for the original derivation of the Bayesian re-weighting [7] (red) and the one proposed in [8,9] (blue). In this plot the number of points is n= 10. The large ensemble of PDFs required by the Bayesian approach can be constructed, in analogue to eq. (3.9), by fk≡fS0+ Neig X i fS+ i−fS− i 2!Rik,(4.9) where the coefficients Rik are random numbers drawn from a Gaussian distribution centered at zero and with variance one. An asymmetric version of eq. (4.9) to account for nonlinearities was advocated in ref. [18]. Specifically, it was proposed that the replicas should be generated by fasym k≡fS0+ Neig X ifS± i−fS0|Rik|(4.10) where fS+ ior fS− iis chosen according to the sign of Rik. However, in this case already before the reweighting the expectation values for the observables will not, in general, match those computed directly with the central set of the original fit. To accurately compare with the linear Hessian reweighting, we stick here to the symmetric prescription of eq. (4.9). As pointed out earlier, the replicas built in this way satisfy the PDF sum rules and the parton evolution equations. After computing the weights ωkfor each replica, the reweighted PDFs can be written as fnew =fS0+ Neig X i fS+ i−fS− i 2!  1 Nrep Nrep X k ωkRik ,(4.11) and, similarly to the Hessian case, one can calculate the “penalty” induced in the original fit by P= ∆χ2 Neig X i  1 Nrep Nrep X k ωkRik  2 .(4.12) not the case or that it has been taken into account e.g. by including suitable penalty terms in the original definition of χ2[28,29]. – 8 – JHEP12(2014)100 where χ2[fS0] is the value of χ2computed with the central set S0, and the coefficients Dik and akwere defined in eqs. (3.3) and (3.7). In the limit of infinitely large Nrep, the sum over the replicas above can be replaced by an integral 1 Nrep Nrep X `=1 Nrep→∞ −−−−−→ (2π)−Neig/2Z+∞ −∞ d~ Rexp −1 2~ R2,(6.4) where the additional exponential stems from the probability distribution for the random numbers Rm` being Gaussian. Using this in eq. (6.3) above, we have wGK k=1 N×(2π)−Neig/2e−1 2∆χ2χ2[fS0]Z+∞ −∞ d~ Re−1 2∆χ2~ RTB~ R−1 ∆χ2~ aT~ RRk =−∆χ2 N×(2π)−Neig/2e−1 2∆χ2χ2[fS0]d dakZ+∞ −∞ d~ Re−1 2∆χ2~ RTB~ R−1 ∆χ2~ aT~ R =−1 N×r(∆χ2)Neig det Be−1 2∆χ2χ2[fS0]+ 1 2∆χ2~ aTB−1~ a(B−1~ a)k, where the elements of the matrix Bwere given in eq. (3.6). The corresponding expression for the denominator Nis N=r(∆χ2)Neig det Be−1 2∆χ2χ2[fS0]+ 1 2∆χ2~ aTB−1~ a.(6.5) Upon taking the ratio, the various prefactors cancel, and the coefficients wGK kreduce to wGK k=−(B−1~ a)k,(6.6) which coincides with the eq. (3.5) specifying the coefficients for the linear Hessian reweighting. As shown in ref. [23], the weights wchi−squared kin eq. (4.8) emerge from wGK kby integrating over all possible data sets that give equal χ2. Being “contaminated” by such additional information readily explains their failure in the present context. 7 Inclusive jet production at the LHC In this section, we apply the reweighting methods to the production of inclusive jets in proton+proton collisions at the LHC. In comparison to our example discussed earlier, a new source of non-linearity arises which could potentially decrease the accuracy of the linear reweighting. Namely, the quadratic PDF dependence of the proton+proton cross sections (ˆσdenotes the coefficient functions and jet definitions in general) σpp[S] = f[S]⊗ˆσ⊗f[S],(7.1) may or may not be well approximated by σpp[S0+1 2X kS+ k−S− kwk]≈σpp[S0] + 1 2X kσpp[S+ k]−σpp[S− k]wk,(7.2) – 15 – JHEP12(2014)100 0.96 0.98 1.0 1.02 1.04 0.96 0.98 1.0 1.02 1.04 0 20 40 60 80 100 120 Data point |y|<0.5 0.5<|y|<1.0 1.0|y|<1.5 1.5<|y|<2.0 2.0<|y|<2.5 Linear approximation / True values Replicas with < 1 0.96 0.98 1.0 1.02 1.04 0.96 0.98 1.0 1.02 1.04 0 20 40 60 80 100 120 Data point |y|<0.5 0.5<|y|<1.0 1.0|y|<1.5 1.5<|y|<2.0 2.0<|y|<2.5 Linear approximation / True values Replicas with > 1 Figure 8.Upper panel: ratios between the jet cross sections computed by the linear approximation of eq. (7.2) divided by the “exact” value using a parametrization of the corresponding linear combination of PDFs inside FASTNLO. The results are plotted for 15 random replicas with the weight-vector length ` < 1. The data points have been numbered with increasing transverse momentum. Lower panel: as the upper panel, but with ` > 1. which is the approximation one effectively makes when using eq. (3.2). In some other cases, such as deep-inelastic scattering (linear in PDFs), this would not be an issue whereas for the W-asymmetry, the non-linearities could be even more intricate. Here, we consider the recent CMS √s= 7 TeV jet measurements [30] (Ndata = 133) for which a direct FASTNLO interface [31–33] is available to evaluate the cross sections at NLO accuracy. We use the CTEQ6.6 PDFs [34] whose error sets quantify the uncertainties within 90% confidence level the corresponding tolerance being ∆χ2 CTEQ6.6(90%) = 100. The renormalization scale µr and factorization scale µfwere fixed to the jet transverse momentum as µr=µf=pT/2, and the strong coupling was set to αs(MZ)=0.118 at the Z boson pole. Let us first discuss the adequacy of the approximation in eq. (7.2). To this end, we have prepared some random PDF replicas by eq. (4.9) separating the cases `=PNeig i=1 R2 ik <1 and `=PNeig i=1 R2 ik >1. We compute the jet cross sections “exactly” (by constructing parametrizations corresponding to these PDF replicas and using them in FASTNLO computations) and, on the other hand, by the linear approximation of eq. (7.2). Typical results from such an exercise are shown in figure 8. First, if ` < 1, the linear approximation proves – 16 – JHEP12(2014)100 0.8 1.0 1.2 1.4 0.8 1.0 1.2 1.4 0.8 1.0 1.2 1.4 0.8 1.0 1.2 1.4 0.8 1.0 1.2 1.4 0.8 1.0 1.2 1.4 0.8 1.0 1.2 1.4 0.8 1.0 1.2 1.4 0.8 1.0 1.2 1.4 0.8 1.0 1.2 1.4 100 1000 CTEQ6.6 error (68%) syst. error pT[GeV] Ratio to CTEQ6.6 |y|<0.5 Ratio to CTEQ6.6 0.5<|y|<1.0 Ratio to CTEQ6.6 1.0<|y|<1.5 Ratio to CTEQ6.6 1.5<|y|<2.0 Ratio to CTEQ6.6 1.5<|y|<2.0 0.8 1.0 1.2 1.4 0.8 1.0 1.2 1.4 0.8 1.0 1.2 1.4 0.8 1.0 1.2 1.4 0.8 1.0 1.2 1.4 0.8 1.0 1.2 1.4 0.8 1.0 1.2 1.4 0.8 1.0 1.2 1.4 0.8 1.0 1.2 1.4 0.8 1.0 1.2 1.4 100 1000 CTEQ6.6 error (68%) syst. error pT[GeV] Ratio to CTEQ6.6 |y|<0.5 Ratio to CTEQ6.6 (systematic shifts applied to data) 0.5<|y|<1.0 Ratio to CTEQ6.6 1.0<|y|<1.5 Ratio to CTEQ6.6 1.5<|y|<2.0 Ratio to CTEQ6.6 1.5<|y|<2.0 Figure 9.Left-hand panels: the CMS inclusive jet cross sections for the five rapidity intervals compared to the NLO calculation with CTEQ6.6 PDFs and taking µf=µr=pT/2. The error bars in the data points show the statistical uncertainty, while the total systematic error is indicated by the blue lines. The colored bands show the CTEQ6.6 PDF uncertainty. Right-hand panels: as the left-hand panels, but after applying the systematic shifts. rather accurate, the deviations being normally much less than couple of percents. In the latter case, ` > 1, the linear approximation evidently breaks down. We conclude that if the reweighted PDFs end up sufficiently close to the original ones the linear approximation of eq. (7.2), and thereby the linear Hessian reweighting, should be rather accurate. In the case of Bayesian reweighting, replicas with ` > 1 often occur, and, from the lower panel of figure 8, we can consequently anticipate some differences whether the cross sections are evaluated using the replicas directly inside FASTNLO, or “on the fly” by the linear approximation of eq. (7.2). We account for the correlated systematic errors by constructing a covariance matrix. To be specific, we compute the elements of the covariance matrix Cby Cij =δij (σuncorr i)2+X k βk iβk j,(7.3) – 17 – JHEP12(2014)100 where σuncorr iis the uncorrelated error of data point i, and βk idenotes the absolute shift of this data point corresponding to 1-sigma deviation of the systematic parameter k. In addition to the luminosity, unfolding and jet energy scale uncertainties, we also treat the quoted uncertainty in the multiplicative non-perturbative corrections (underlying event, hadronization) as a correlated systematic error. The uncorrelated errors σuncorr iinclude the statistical and 1% uncorrelated systematic uncertainty added in quadrature. Calculating the χ2using the covariance matrix Cis equivalent to (see e.g. [35–37]) minimizing χ2=X i"ytheory i−ydata i−Pkskβk i σuncorr i#2 +X k s2 k,(7.4) with respect to the systematic parameters sk. This occurs with the parameter values smin k=X j βk j−X i,`,s βk iC−1 i` βs `βs j  ytheory j−ydata j (σuncorr j)2,(7.5) and −Pksmin kβk iis the net systematic shift for the data point ydata i. Figure 9presents a comparison between the CMS data and the NLO predictions, including the 68% PDF error bands (obtained downscaling the original 90% errors by 1/1.645 [11,38]). The χ2value for the central CTEQ6.6 set is χ2/Ndata ≈2.1 accounting for the data correlations as discussed above. The calculation appears to overpredict the experimental cross section by some 5% which, however, gets still easily hidden under the systematic shifts as demonstrated in figure 9. It is a straightforward task to apply the reweighting methods on these jet data. For the needs of the Bayesian techniques we have generated 104PDF replicas using eq. (4.9). As the data uncertainties correspond to one standard deviation, we have rescaled the random numbers in eq. (4.9) by Rik →Rik/1.645 which brings the PDF replicas to the 68% level as well. When computing the GK weights we divide the resulting values of χ2by ∆χ2 CTEQ6.6(68%) = 100/1.6452≈37, to appropriately modify the underlying likelihood. The chi-squared weights are always computed without rescaling the χ2values. We stress that the rescaling of CTEQ6.6 replicas from 90% to 68% confidence level affects only the outcome of chi-squared reweighting as, by construction, the central result of GK reweighting would remain the same by simply using the 90% replicas and rescaling the χ2values by ∆χ2 CTEQ6.6(90%). The Hessian reweighting has been performed directly with the original 90% error sets rescaling the uncertainties by 1/1.645 only at the end. The resulting modifications in the gluon PDFs are presented in the upper panel of figure 10 (as the jet production is predominantly sensitive to the gluons, we find it reasonable to present the results here only in terms of gluon PDFs). From this plot we find that, as expected, the Hessian technique and the Bayesian method with rescaled GK weights are in good agreement also here. The small mismatch at large xoriginates mainly from eq. (7.2) not being exact. As the reweighting penalty P≈21 is less than ∆χ2 CTEQ6.6(68%), we conclude that these data could have been added to the CTEQ6.6 fit without causing a significant disagreement with the original data. In this sense these data are compatible with the – 18 – JHEP12(2014)100 Figure 10.Upper panel: the gluon PDFs at Q2= 10000 GeV2after reweighting using the correlated errors. The red points with error bars correspond to the results using the Bayesian reweighting with rescaled GK weights, and the blue dotted lines mark the corresponding result with chi-squared weights. The result from Hessian reweighting is indicated by the black dashed lines and the colored band is the original CTEQ6.6 uncertainty. All results are normalized to the central set of CTEQ6.6. Lower panel: the distribution of Bayesian weights. – 19 – JHEP12(2014)100 0.8 1.0 1.2 1.4 0.8 1.0 1.2 1.4 0.8 1.0 1.2 1.4 0.8 1.0 1.2 1.4 0.8 1.0 1.2 1.4 0.8 1.0 1.2 1.4 0.8 1.0 1.2 1.4 0.8 1.0 1.2 1.4 0.8 1.0 1.2 1.4 0.8 1.0 1.2 1.4 100 1000 CTEQ6.6 error (68%) Hessian reweighting pT[GeV] Ratio to CTEQ6.6 |y|<0.5 Ratio to CTEQ6.6 0.5<|y|<1.0 Ratio to CTEQ6.6 1.0<|y|<1.5 Ratio to CTEQ6.6 1.5<|y|<2.0 Ratio to CTEQ6.6 1.5<|y|<2.0 0.6 0.8 1.0 1.2 1.4 0.6 0.8 1.0 1.2 1.4 0.6 0.8 1.0 1.2 1.4 0.6 0.8 1.0 1.2 1.4 0.6 0.8 1.0 1.2 1.4 0.6 0.8 1.0 1.2 1.4 0.6 0.8 1.0 1.2 1.4 0.6 0.8 1.0 1.2 1.4 0.6 0.8 1.0 1.2 1.4 0.6 0.8 1.0 1.2 1.4 100 1000 CTEQ6.6 error (68%) chi-squared reweighting pT[GeV] Ratio to CTEQ6.6 |y|<0.5 Ratio to CTEQ6.6 0.5<|y|<1.0 Ratio to CTEQ6.6 1.0<|y|<1.5 Ratio to CTEQ6.6 1.5<|y|<2.0 Ratio to CTEQ6.6 1.5<|y|<2.0 Figure 11.Left-hand panels: the CMS inclusive jet cross sections compared to the predictions after the Hessian reweighting and applying the systematic shifts. All values have been normalized to the central prediction of CTEQ6.6 and only the statistical data errors are shown. Right-hand panels: as the left-hand panels, but with the Bayesian reweighting with chi-squared weights. CTEQ6.6 PDFs despite the largish χ2/Ndata ≈1.75 which could hint, however, that some tension between the Tevatron Run-1 jet data [39,40] (used in CTEQ6.6 fit to constrain large-xgluons) and these new LHC measurements exist.5The line of red points in figure 10 shows the spectrum of rescaled GK weights. For dividing the values of χ2by the tolerance ∆χ2 CTEQ6.6(68%) ≈37, this distribution is rather flat and the final result is affected by replicas with a wide range of χ2explaining also the rather large number of effective replicas, Neff/Nrep ≈0.72. The chi-squared weights give rise to effects which are qualitatively alike but quantitatively much larger. The agreement with the new jet data is admittedly better, χ2/Ndata ≈1.15 (the green points in figure 10 indeed demonstrate that the chi-squared weights assign clearly larger weights to replicas with χ2/Ndata <1.4 than the rescaled GK weights), but at the cost of increasing the original χ2by P≈290 which is way beyond the 5Indeed, the latest CTEQ fit [35] has abandoned the Tevatron Run-1 jet data as they do not completely agree with the Run-2 data. See ref. [5] for a review and ref. [41] for more discussions. – 20 – JHEP12(2014)100 CTEQ6.6 tolerance (both 68% and 90%). The surviving number of replicas Neff is also very low. That is, one could conclude that the jet data considered here were not in agreement with the CTEQ6.6 PDFs. However, as we have demonstrated in the previous sections, the result is misleading since the chi-squared weights do not appropriately account for the other data that were originally included in the global CTEQ6.6 fit. From eq. (5.4) we can estimate the value of ∆χ2that would have brought the results of chi-squared reweighting close to the correct ones: using the original central value χ2/Ndata ≈2.1 and Ndata = 133 we obtain ∆χ2≈1.9. This is much less than ∆χ2 CTEQ6.6(68%) = 37 and concretely explains why the chi-squared weights lead to a result which is so far from the correct one. The new predictions for the cross-section are shown in figure 11. The systematic shifts have been applied to the data and they depend significantly on which method of reweighting was used. The Hessian reweighting has caused a mild downward shift on the cross section which was to be expected given that the original central predictions somewhat overshoot the data (see figure 9). In the case of Bayesian reweighting with chi-squared weights the induced changes with respect to the original CTEQ6.6 predictions are way too large. Reweighting of MSTW2008. In the previous example, we discussed PDF reweighting in the case of a fixed global tolerance ∆χ2assuming the validity of the quadratic approximation in eq. (3.1). It is, however, well known that such quadratic profile is never perfect, and for the fit parameters that are not very well constrained severe deviations from the ideal quadratic behaviour may occur (see e.g. figures 5-6 in [5]). To account for such imperfections and improve the linear approximations in estimating the observables in the space of eigenvectors, the non-linear extension of the Hessian reweighting should provide some improvement. As an example, we consider the MSTW2008 PDFs [5] using the same CMS jet data as earlier. In fact, the effect of incorporating the √s= 7 TeV CMS and ATLAS [42] inclusive jet measurements into the MSTW2008 framework was studied recently via a direct re-fit [21], but due to the large systematic uncertainties of the ATLAS data, the fit was mainly driven by the same CMS measurements we have discussed here. Thus, in what follows, we will contrast our results to what was obtained in that analysis. The MSTW2008 package provides two separate sets of error PDFs, the 68% and 90% confidence-level sets, of which we will use the former one. As the numerical values for the determined tolerances t± kand the corresponding increases T± kof the global χ2are given, we can account also for the non-quadratic behaviour of the original χ2function in the space of eigenvectors. Specifically, in addition to using eq. (3.13) in evaluating the observables, we make the following substitution in eq. (3.1): Neig X k z2 k−→ Neig X k akz2 k+bkz3 k, ak=t− k(T+ k/t+ k)2−t+ k(T− k/t− k)2 t+ k+t− k ,(7.6) bk=(T+ k/t+ k)2−(T− k/t− k)2 t+ k+t− k . – 21 – JHEP12(2014)100 Figure 12. The reweighted gluon distributions at Q2= 10000 GeV2normalized by the central MSTW2008 gluons. The solid black line is the result from the non-linear Hessian reweighting, and the green dashed-dotted line corresponds to the linear Hessian reweighting. Blue dotted line is the prediction of Bayesian reweighting with chi-squared weights, and red dashed line results from using the rescaled GK weights. The shaded band is the original MSTW2008 uncertainty and the purple squares mark the result obtained in ref. [21]. We set αs(M2 Z) = 0.1202, and use µr=µf=pT. As it is not possible (with the available information) to consistently implement the dynamic tolerance after the reweighting, we discuss here only the resulting central values. The gluon PDF after the non-linear Hessian reweighting is shown in figure 12. For comparison, the figure includes also the outcome from the linear Hessian reweighting, Bayesian reweighting with the chi-squared weights, Bayesian reweighting with GK weights rescaled by the average tolerance, ∆χ2 MSTW2008(68%) ≈10, and the result obtained by a refit in ref. [21]. Apart from the Bayesian predictions with chi-squared weights the results of different prescriptions are more or less consistent with each other within 10−4<x<10−1, and also close to the “exact” result.6The differences between the linear and non-linear prescriptions are clearly more important at x < 10−4and x > 10−1, where the original PDF uncertainties become large. The chi-squared weights predict again modifications for the gluon PDFs which are similar to the correct ones but somewhat too pronounced (the estimated penalty stands out from what is obtained with the other methods). The reason can be understood by estimating from eq. (5.4) which value of ∆χ2would have given consistent results: now χ2/Ndata ≈1.34 which translates to ∆χ2≈3.9. While close, this is still less than ∆χ2 MSTW2008(68%) ≈10 and explains why the result is different from the direct re-fit. 6Note that in ref. [21] the open fit parameters were not restricted to those used in generating the 20 eigenvecor sets. Thus, it is natural to expect deviations in comparison to the results of reweighting. – 22 – JHEP12(2014)100 0.9 1.0 1.1 1.2 0.9 1.0 1.1 1.2 0.9 1.0 1.1 1.2 0.9 1.0 1.1 1.2 0.8 0.9 1.0 1.1 1.2 100 1000 MSTW2008 error Hessian reweighting (68% error sets) EPJ C74 (2014) 2934 pT[GeV] Ratio to MSTW2008 |y|<0.5 Ratio to MSTW2008 0.5<|y|<1.0 Ratio to MSTW2008 1.0<|y|<1.5 Ratio to MSTW2008 1.5<|y|<2.0 Ratio to MSTW2008 1.5<|y|<2.0 0.9 1.0 1.1 1.2 0.9 1.0 1.1 1.2 0.9 1.0 1.1 1.2 0.9 1.0 1.1 1.2 0.9 1.0 1.1 1.2 100 1000 MSTW2008 error chi-squared reweighting (68% error sets) EPJ C74 (2014) 2934 pT[GeV] Ratio to MSTW2008 |y|<0.5 Ratio to MSTW2008 0.5<|y|<1.0 Ratio to MSTW2008 1.0<|y|<1.5 Ratio to MSTW2008 1.5<|y|<2.0 Ratio to MSTW2008 1.5<|y|<2.0 Figure 13.Left-hand panels: the jet cross sections after reweighting MSTW2008 PDFs. The results of the non-linear Hessian reweighting are shown by the red lines and the results of ref. [21] are the black dashed lines. All values have been normalized to the central predictions of MSTW2008. Systematic shifts have been applied to the data and only the statistical data errors are shown. Right-hand panels: as the left-hand panels, but with the Bayesian reweighting with chi-squared weights. The predicted jet cross sections after the reweighting are shown in figure 13, normalized to the central MSTW2008 values. The results of Hessian reweighting clearly stay within the original uncertainties and appear to agree with the results of ref. [21]. The predictions of the Bayesian reweighting with chi-squared weights are, however, at the lower limit of the original error bands and disgaree with the true re-fit. 8 Conclusion We have discussed how to test the consistency and estimate the effects that a new set of experimental measurements have within an existing set of Hessian PDFs with non-zero tolerance ∆χ2. To this end, the Hessian reweighting, was introduced as an alternative technique to the prevailing Bayesian methods. While the Hessian reweighting is straight- – 23 – JHEP12(2014)100 forwardly derived considering a new set of data in a global χ2fit, the Bayesian methods are outwardly distinct, based on statistical inference. We compared the different approaches to a direct re-fit through a simple example verifying the adequacy of the Hessian method. In the case of the Bayesian procedure an agreement with a new fit was also established, but only after including the ∆χ2criterion properly into the Bayesian likelihood function which — as we mathematically justified — must be a pure exponential as originally proposed by Giele and Keller. The conditions under which another commonly used (but in this case inadequate) likelihood function gives consistent results was also discussed. The inclusive jet production at the LHC was considered as an additional example. At first, our findings may appear to be in contrast to the works of NNPDF collaboration in which a different functional form for the likelihood is derived. However, as the NNPDF methodology for fitting PDFs is far more involved than the simple χ2minimization considered here, it is possible that a different functional form applies. For example, the NNPDF fits involve extremely flexible fit functions and to avoid fitting random fluctuations the data are divided into “training“ and “control” sets. While the actual χ2minimization is done only for the former subset of data, the latter one is used to decide when the minimization process should be stopped. The outcome of a direct χ2minimization with no stopping criteria would generally be different if the same functional form would be used. However, for the moment it is still unclear what exactly causes the need for using a different Bayesian likelihood. The two types of methods discussed here have their pros and cons: while the Hessian procedure requires evaluating the observables only with the central and error sets (typically around 50 sets in total), in the Bayesian method one needs to deal with a much larger ensemble of PDFs (around 103to find well-converging results). On the other hand, the Hessian reweighting is procedurally a bit more involved requiring e.g. numeric linear algebra, while the Bayesian technique is simpler. The reliability of the both methods depends basically on the accuracy of the quadratic approximation around the minimum χ2made in the original PDF fit and on the adequacy of the linear approximations that one makes. In the case of Hessian reweighting one can improve on these approximations by invoking non-linear corrections as we explained. However, when the PDF fits are updated details like the form of the fit function, data sets included, or value of αs(M2 Z) are often altered. Such modifications cannot be easily accounted for by the reweighting procedures and in this sense the reweighting is always an approximation to a real fit. In any case, it is definitely useful in checking whether the new data appear consistent and which components of the PDFs are bound to undergo a change and how. Acknowledgments We thank J. Rojo and S. Forte for useful discussions on the PDF reweighting and N. Armesto and C. Salgado for their comments regarding the manuscript. We also want to acknowledge the Referee of our paper whose suggestions led to significant improvements especially in section 5 concerning the discussion on the conditions that make the chi-squared reweighting to closely mimic a direct refit. During this work H.P. was supported by the Academy – 24 –