Specialized minimal PDFs for optimized LHC calculations.

We present a methodology for the construction of parton distribution functions (PDFs) designed to provide an accurate representation of PDF uncertainties for specific processes or classes of processes with a minimal number of PDF error sets: specialized minimal PDF sets, or SM-PDFs. We construct these SM-PDFs in such a way that sets corresponding to different input processes can be combined without losing information, specifically as regards their correlations, and that they are robust upon smooth variations of the kinematic cuts. The proposed strategy never discards information, so that the SM-PDF sets can be enlarged by the addition of new processes, until the prior PDF set is eventually recovered for a large enough set of processes. We illustrate the method by producing SM-PDFs tailored to Higgs, top-quark pair, and electroweak gauge boson physics, and we determine that, when the PDF4LHC15 combined set is used as the prior, around 11, 4, and 11 Hessian eigenvectors, respectively, are enough to fully describe the corresponding processes.

Introduction

Modern sets of parton distributions (PDFs) [1-6] provide a representation of their associated uncertainties based on either the Hessian [7] or the Monte Carlo (MC) [8] methods, supplementing their central PDF member with additional error members (eigenvectors or MC replicas). The number of PDF members required for an accurate representation of PDF uncertainty can be as large as several hundreds, especially when constructing PDF sets based on the combination of several underlying PDFs fitted to data: for example, the recent PDF4LHC 2015 sets [9] are based on a combined sample of 900 MC PDF replicas.

The usage of such large PDF samples can be computationally unwieldy, and this motivated the development of strategies for reducing the number of PDF members while minimizing accuracy loss. A number of such reduction strategies have been made available recently. Two of these methods provide a Hessian representation of the prior PDF set in terms of a smaller number of eigenvectors: META-PDFs [10] and MCH-PDFs [11]. A third method uses a compression algorithm to reduce the number of replicas of an underlying MC PDF prior: CMC-PDFs [12].

These three methods have been extensively benchmarked in the context of the 2015 PDF4LHC recommendations [9], where it was found that generally a set of about 100 PDFs is required in order to represent PDF uncertainties with percentage accuracy for all PDFs in the complete range of (x, Q) relevant for LHC phenomenology. However, it is well known [13] that, if one is interested only in the description of a specific set of cross sections, the number of PDF error members can be greatly reduced without significant accuracy loss.

In this work we propose a new strategy to achieve this goal. Our methodology, which we denote by Specialized Minimal PDFs (SM-PDFs), is based on the Singular Value Decomposition version of the mc2hessian algorithm, as presented in the appendix of Ref. [11]. Starting from either a Hessian or a Monte Carlo prior set and a list of collider processes, the SM-PDF algorithm leads to a set of eigenvectors optimized for the description of the input processes within some given tolerance.

In comparison to existing methods, such as data set diagonalization [13], our methodology has the advantage that no information is lost in the process of the construction of the specialized set. This is because the specialized set is constructed through a suitable linear transformation, whereby the starting space is separated into a subspace spanned by the optimized SM-PDF set, and its orthogonal subspace. This then implies that any given SM-PDF set can be iteratively expanded in order to maintain a given accuracy for an increasingly large set of processes, and also that SM-PDF sets optimized for different sets of processes can be combined into a single set, either a priori, at the level of PDFs, or a posteriori, at the level of cross sections. This, for example, enables the a posteriori combination of previous independent studies for a signal process and its corresponding backgrounds, with all correlations properly accounted for.

This paper is organized as follows: in Sect. 2 we describe our general strategy and methodology in detail. Then, in Sect. 3, we apply our method to the most important Higgs production channels (ggh, and hV, VBF h) as well as for other standard candles at the LHC, i.e. , Z, and W production. We compute one specific reduced sets for each of them, as well as a single set for all the processes combined. We validate the results by comparing the predictions of these reduced sets to the prior input set. We also show that our method provides an adequate generalization by showing that the predictions are stable when computing similar processes but with different kinematical cuts from those used as input. In Sect. 4 we show how experimental analyses done with different SM-PDFs can be combined together. In Sect. 5 we provide an overview of the deliverables of this work, in particular the code itself which allows one to easily generate reduced sets with personalized configuration and the LHAPDF6 [14] sets of SM-PDFs for the processes described in Sect. 5. Finally, Appendix A presents a graphical illustration of the regions in PDF space which give the dominant contribution to various physical processes, and Appendix B provides some basic instructions for the execution of the SM-PDF code.

Methodology

The SM-PDF methodology is built upon the strategy based on Singular-Value Decomposition (SVD) followed by Principal Component Analysis (PCA) described in the appendix of Ref. [11], in which the MCH method was presented. This SVD+PCA strategy achieves the twofold goal of obtaining a multigaussian representation of a starting (prior) Monte Carlo PDF set, and of allowing for an optimization of this representation for a specific set of input cross sections, which uses the minimal number of eigenvectors required in order to reach a desired accuracy goal. We will now review the SVD+PCA method, and describe how it can be used for the construction of specialized minimal PDF sets, optimized for the description of a specific set of cross sections.

The SVD+PCA method

The main problem we are addressing is the faithful representation of PDF uncertainties, which typically requires a large number of PDF error or Monte Carlo sets. Here we will assume the central value to be the same as in the prior PDF set, from which, if the prior is given as a Monte Carlo, it is typically determined as a mean (though different choices, such as the median, are possible and might be advisable in particular circumstances).

Hence, we are interested in the construction of a multigaussian representation in PDF space: the only information we need is then the corresponding covariance matrix. This is constructed starting with a matrix X which samples over a grid of points the difference between each PDF replica, , and the central set, , namelywhere runs over the independent PDF flavors at the factorization scale , i runs over the points in the x grid where the PDFs are sampled, runs over all grid points, and k runs over the replicas. The sampling is chosen to be fine-grained enough that results will not depend on it.

The desired covariance matrix in PDF space is then constructed asThe key idea which underlies the SVD method is to represent the covariance matrix Eq. (2) over the dimensional linear space spanned by the replicas (assuming ), by viewing its eigenvectors as orthonormal basis vectors in this space, which can thus be represented as linear combinations of replicas. The subsequent PCA optimization then simply consists of picking the subspace spanned by the dominant eigenvectors, i.e., those with largest eigenvalues.

The first step is the SVD of the sampling matrix X, namelywhere U and are orthogonal matrices, with dimensions respectively and , S is a diagonal positive semi-definite matrix, whose elements are the so-called singular values of X, and the initial number of singular values is given by . Note that, because S is diagonal, it can be equivalently viewed as a matrix, since (with ) all its further entries vanish. This point of view was taken in the appendix of [11]. In this case, only the submatrix which actually contributes to the SVD of the matrix V is included. However, for the procedure to be described below, it is more convenient to view V as orthogonal matrix.

The matrix then has the important property thatbut also that it can be expressed asand thus it provides the sought-for representation of the multigaussian covariance matrix in terms of the original PDF replicas: specifically, is the expansion coefficient of the jth eigenvector over the kth replica. We assume henceforth that the singular values are ordered, so that the first diagonal entry of S correspond to the largest value, the second to the second-largest and so forth.

The PCA optimization then consists of only retaining the principal components, i.e. the largest singular values. In this case, U and S are replaced by their sub-matrices, denoted by u and s, respectively, with dimension and , with the number of eigenvectors which have been retained. Due to the ordering, these are the upper left sub-matrices. Because s has only non-vanishing diagonal entries, only the submatrix of V contributes. We call this the principal submatrix P of V:The optimized representation of the original covariance matrix, Eq. (2), is then found by replacing V with its principal submatrix P in Eq. (5). This principal matrix P is thus the output of the SVD+PCA method: it contains the coefficients of the linear combination of the original replicas or error sets which correspond to the principal components, which can be used to compute PDF uncertainties using the Hessian method.

Indeed, given a certain observable (which could be a cross section, the value of a structure function, a bin of a differential distribution, etc.) its PDF uncertainty can be computed in terms of the original Monte Carlo replicas bywhere is the prediction obtained using the kth Monte Carlo PDF replica, is the central prediction, and in the last step we have defined the vector of differenceswith normAssuming linear error propagation and using Eq. (5), the norm of the vector , Eq. (8), can be represented in the eigenvector basis:where the rotated vectorhas the same norm as the original one because of Eq. (4).

Replacing V by the principal matrix P in Eq. (11), i.e., letting j only run up to we getwhere now the vector is both rotated and projectedThe norm of is only approximately equal to that of the starting vector of differences d: . However, it is easy to see that this provides the linear combination of replicas which minimizes the difference in absolute value between the prior and final covariance matrix for the given number of eigenvectors. As the difference decreases monotonically as increases, the value of can be tuned to any desired accuracy goal, with the exact equality Eq. (10) achieved when . Note that, of course, the optimization step can be performed also starting with a symmetric Hessian, rather than Monte Carlo, prior. In such a case, the index k runs over Hessian eigenvectors, Eq. (2) is replaced by , and the rest of the procedure is unchanged.

An interesting feature of this SVD+PCA method is that the matrix V (and thus also the principal matrix P) in Eq. (11) does not depend on the value of the PDF factorization scale Q: the scale dependence is thus entirely given by the DGLAP evolution equation satisfied by the original Monte Carlo replicas. The result of the SVD thus does not depend on the scale at which it is performed. Of course, the subsequent PCA projection may depend on scale if there are level crossings, but this is clearly a minor effect if a large enough number of principal components is retained. Because of this property, the SVD+PCA methodology can be used for the efficient construction [9] of a Hessian representation of combined PDF sets, even when the sets which enter the combination satisfy somewhat different evolution equations, e.g., because of different choices in parameters such as the heavy quark masses, or in the specific solution of the DGLAP equations.

The SM-PDF method

In the SM-PDF method, this same SVD+PCA optimization is performed, but now with the goal of achieving a given accuracy goal not for the full prior PDF set in the complete range of x and , but rather for the aspects of it which are relevant for the determination of a given input set of cross sections, and in such a way that all the information which is not immediately used is stored and can be a posteriori recovered either in part or fully, e.g. if one wishes to add further observables to the input list.

This requires supplementing the SVD+PCA methodology of Ref. [11] with three additional features: a measure of the accuracy goal; a way of singling out the relevant part of the covariance matrix; and a way of keeping the information on the rest of the covariance matrix in such a way that, if needed, the full covariance matrix can be recovered at a later stage.

The main input to the algorithm is the set of observables which we want to reproduce, , with . Theoretical predictions for the cross sections are computed using a prior PDF set, which we assume for definiteness to be given as a Monte Carlo, though the method works with obvious modifications also if the starting PDFs are given in Hessian form. The goal of the SM-PDF methodology is to evaluate the PDF uncertainties , Eq. (7), in terms of a reduced number of Hessian eigenvectors,with the number being as small as possible within a given accuracy. We thus define a measure of the accuracy goal (tolerance) by the conditionin other words, is the maximum relative difference which is allowed between the original and reduced PDF uncertainties, and , respectively, for all the observables .

In order to determine the part of the covariance matrix relevant for the description of the input observables , we define the correlation functionwhere the matrix of PDF differences X(Q) and the grid index have been defined in Eq. (1); is the standard deviation of the PDFs in the prior Monte Carlo representation, given by the usual expression,and , the standard deviation of the ith observable , is given by Eq. (7). The function Eq. (16) measures the correlation between the observables and the lth PDF value (i.e. , with ).

The basic idea of the SM-PDF construction is to apply the SVD to the subset of the covariance matrix which is most correlated to the specific observables that one wishes to reproduce, through a procedure such that information is never discarded, so observables can be added one at a time, or at a later stage. This goal is achieved through an iterative procedure schematically represented in Fig. 1, which we now describe in detail.

Fig. 1

Schematic representation of the SM-PDF strategy

The iteration loop (contained in the dashed box in Fig. 1) is labeled by an iteration index j, such that at each iteration an extra eigenvector is added, thereby increasing the accuracy. If the accuracy goal is achieved for all observables after j iterations, then the final reduced Hessian set contains eigenvectors as error sets. These are obtained as a new principal matrix P, which provides the expansion coefficients of the eigenvectors over the replica basis: namely, is the component of the jth eigenvector in terms of the kth replica. They thus replace the principal matrix of the previous PCA procedure as a final output of the procedure, and they can be used in exactly the same way.

To set off the iterative procedure, we select one of the observables we wish to reproduce from the list, , and compute the correlation coefficient for all grid points and for a suitable choice of the scale Q. We then identify the subset of grid points for which exceeds some threshold value:The threshold value is expressed as a fraction times the maximum value that the correlation coefficient takes over the whole grid, thereby making the criterion independent of the absolute scale of the correlation. The choice of the scale Q and the threshold parameter t should be taken as tunable settings of the procedure, and it will be discussed in Sect. 3 below. For the time being it suffices to say that Q should be of the order of the typical scale of the observable (for example, the average value of the factorization scale).

We then construct a reduced sampling matrix , defined as in Eq. (1), but now only including points in the space which are in the subset . We perform the SVD of the reduced matrixand we only keep the largest principal component, i.e. one single largest eigenvector, which is specified by the coefficients of its expansion over the replica basis, namely, assuming that the singular values are ordered, by the first row of the V matrix. We thus start filling our output principal matrix P by lettingNote that j on the left-hand side labels the eigenvector ( provides expansion coefficients for the jth eigenvector) and on the right-hand side it labels the iteration ( is the first row of the V-matrix at the jth iteration), which we can identify because, as mentioned, at each iteration we will add an eigenvector. The remaining eigenvectors of the principal matrix span the linear subspace orthogonal to P, and we assign them to a residual matrix R:At the first iteration, when there is only one eigenvector, the principal matrix P has just one row, and it coincides with the principal component of V. So far, the procedure is identical to that of the SVD+PCA method, and we can thus use again Eq. (12) to compute uncertainties on observables, check whether the condition Eq. (15) is met, and if it is not, add more eigenvectors. The procedure works in such a way that each time a new eigenvector is selected, exactly the same steps are repeated in the subspace orthogonal to that of the previously selected eigenvectors, thereby ensuring that information is never discarded. This is achieved by a projection method.

Specifically, we project the matrix X and the vector of observable differences on the orthogonal subspace of P, namely, the space orthogonal to that spanned by the eigenvectors which have already been selected (as many as the number of previous iterations). The projections are performed by, respectively, replacing and X by where the first iteration of the residual matrix has been defined in Eq. (21).

After the projection, we proceed as in the first iteration. We first determine again the subset , Eq. (18), of the projected sampling matrix , thereby obtaining a new sampling matrix : this is possible because everything is expressed as a linear combination of replicas anyway. Once the new matrix has been constructed, the procedure is restarted from Eq. (19), leading to a new matrix . The number of columns of the projected matrix (and therefore of ) is , which is the dimension of the subspace of the linear combinations not yet selected by the algorithm (that is, for , and so on). We can now go back to Eq. (20) and proceed as in the previous case, but with the projected matrices: we add another row to the matrix of coefficients to the principal matrix by picking the largest eigenvector of the projected matrix, and determining again the orthogonal subspace.

At the jth iteration, this procedure gives which, respectively, generalize Eqs. (20) and (21) for . The projected row of coefficients Eq. (24) can be used to determine the corresponding unprojected row of coefficients of the principal matrix and of the residual matrix by using the projection R matrix in reverse, i.e., at the jth iteration We thus end up with a principal matrix which has been filled with a further eigenvector, and a new residual matrix and thus a new projection.

In summary, at each iteration we first project onto the residual subspace, Eq. (22), then pick the largest eigenvector in the subspace, Eq. (24), and then re-express the results in the starting space of replicas, Eq. (26), so that P is always the first row of V in each subspace, and Eqs. (13)–(12) remain valid as the P matrix is gradually filled. Determining the correlation and thus after projection ensures that only the correlations with previously unselected linear combinations are kept. The fact that we are always working in the orthogonal subspace implies that the agreement for the observables , which had already been included, can only be improved and not deteriorated by subsequent iterations. It follows that we can always just check the tolerance condition on one observable at a time. The procedure is thus unchanged regardless of whether we are adding a new observable or not. In any case, the subset for Eq. (18) is always determined by only one observable, namely, the one that failed to satisfy the tolerance condition at the previous iteration. The procedure is iterated until the condition is satisfied for all observables in the input list. The number of iterations j until convergence defines the final number of eigenvectors .

The output of the algorithm is the final principal matrix P, which can be used to compute uncertainties on observables using Eqs. (12)–(13). However, for the final result we wish to obtain a set of Hessian eigenvectors. These can be obtained by performing the linear transformation given by P (a rotation and a projection) in the space of PDFs. The X matrix Eq. (1) then becomesso, substituting in Eq. (1), the final eigenvectors are found to be given byThis is the same result as with the SVD+PCA algorithm of Sect. 2.1, but now generally with a smaller number of eigenvectors, namely, those which are necessary to describe the subset of the covariance matrix which is correlated to the input set of observables.

SM-PDF usage and optimization

Upon delivery of the final PDF set, any observable is computed in terms of the resulting Hessian representation Eq. (29). As in the case of the original SVD+PCA methodology, the final result Eq. (29) determines the PDFs for all x and Q. Indeed, Eq. (29) determines the SM-PDF Hessian eigenvectors as linear combinations of replicas, and thus for all values of x and Q for which the original replicas were defined.

Note, however, that in the procedure of Sect. 2.2, in order to test for the tolerance criterion observables have been computed using Eqs. (12)–(13). This is equivalent to using the PDFs Eq. (28) by standard linear error propagation, but it differs from it by nonlinear terms, specifically for hadron collider processes in which observables are quadratic in the PDFs. Even though nonlinear corrections are expected to be small, in principle it could be that the tolerance criterion is no longer satisfied if Eq. (28) is used instead.

We explicitly check for this, and if it is the case for all observables such that the recomputed tolerance criterion is not satisfied, we restart the iteration but now replacing the tolerance with a new value given bywhere is the value of the tolerance that is obtained within the linear approximation, by computing Eq. (15) with Eq. (12). Iterating until the criterion with the new tolerances Eq. (30) is met will be sufficient to ensure that the tolerance criterion is satisfied when using the new PDFs, provided the difference between the linear and exact estimate of is mostly due to the larger eigenvectors that were selected first and remains approximately constant upon addition of smaller eigenvectors in order to correct for this.

In practice, the difference between the linear estimation of the PDF uncertainty and the exact result is generally small, and does not a change the result for target tolerances of 5 % or bigger. This effect can be more important for observables affected by substantial PDF uncertainties, or for processes which depend on a large number of partonic channels (especially when new channels open up at NLO or NNLO). It is, however, not an issue for most practical applications.

Note that this final optimization step may become extremely time consuming if fast grid tools are not available. In view of this, it is possible to disable this check. However, fast interfaces can be obtained for any NLO QCD cross section with arbitrary final-state cuts using the aMCfast interface [15] to Madgraph5_aMC@NLO [16].

The SM-PDF construction can be generally performed at any perturbative order, and specifically starting with an NLO or an NNLO PDF set. The perturbative order enters both in the choice of starting PDF set, and in the computation of the list of observables , specifically used for the determination of the correlation function defined in Eq. (16). Because the NNLO-NLO K factors are usually moderate, for most applications it may be sufficient to compute using NLO theory even when using NNLO PDFs throughout. An obvious exception is the case in which the user is explicitly interested in studying the changes in PDFs when going from NLO to NNLO.

A final issue is whether results depend on the order in which the observables are included, and specifically on the choice of the observable used to start the iteration. Indeed, the eigenvectors selected for a specific iteration depend on the subspace spanned by the previous eigenvectors, and consequently a different ordering will indeed change the particular linear combinations that are selected. However, this does not significantly affect the total number of eigenvectors needed, because the optimal subspace of linear combinations required to describe all observables with a given accuracy remains the same regardless of the order they are presented. We have verified that this is indeed the case, though we observed small fluctuations by one or two units in the final number of eigenvectors due to the discontinuous nature of the tolerance criteria, Eq. (15).

Results and validation

We now present the validation of the SM-PDF algorithm described in the previous section. Using this methodology, we have constructed four specialized minimal PDF sets for different representative cases of direct phenomenological relevance at the LHC:These examples have been chosen since, for each SM-PDF, there is a strong case for the use of optimized PDF sets with a greatly reduced number of eigenvectors. For instance, these SM-PDFs could be of interest for studies of the Higgs Cross-Section Working Group [17] (case 1), the LHC Top Working Group (case 2), and the LHC Electroweak Working Group (case 3), respectively. As an example, the SM-PDFs for W, Z production could be relevant for the determination of the W boson mass [18-20], which is a extremely CPU-time consuming task.

Higgs physics.

Top quark pair production physics.

Electroweak gauge boson production physics.

The combination of all processes included in (1), (2), and (3).

In this section, we will first define the PDF priors and LHC cross sections that have been used to construct the SM-PDF sets listed above, then validate the performance of the algorithm using a variety of figures of merit.

Input PDFs and cross sections

In order to validate the SM-PDF methodology, we have used three different prior PDF sets, all of them in the Monte Carlo representation:These three choices are representative enough for the validation of our methodology; they show that the procedure works regardless of the choice of input PDF set. As already mentioned in Sect. 2.3 the SM-PDF methodology can be applied equally to NLO or NNLO PDFs, and NLO PDFs are chosen here purely for the sake of illustration. Indeed, in Appendix B we provide an example in which NNLO PDFs are used.

The NNPDF3.0 NLO set [6] with replicas,

The MMHT14 NLO set [5] with replicas, obtained from the native Hessian representation using the Watt-Thorne method [21], and

The PDF4LHC 2015 NLO prior set [9], with replicas, built from the combination of 300 replicas from each of the CT14, MMHT14 and NNPDF3.0 NLO sets. This set is denoted by MC900 in the following.

In order to compute the theoretical predictions for all input PDF sets and as many cross sections as possible, we have generated a large number of dedicated APPLgrid grids [22] using the aMCfast [15] interface to MadGraph 5_aMC@NLO [16]. Cross sections and differential distributions have been computed for the LHC Run II kinematics, with a center-of-mass energy of TeV. In particular we have generated fast NLO grids for the following processes:A more detailed description of these processes, including binning and the kinematical cuts applied, is provided in Tables 1, 2, and 3. We also indicate the names of the (publicly available) APPLgrid grids generated for the present validation study. Producing fast NLO grids for additional processes, or with a different binning or set of analysis cuts, is straightforward using the aMC@NLO/aMCfast framework. We adopt the default choice of renormalization and factorization scales in aMC@NLO, namely , withthe scalar sum of the transverse masses of all final-state particles at the matrix-element level.

Table 1

LHC processes and the corresponding differential distributions that have been used as input in the construction of the SM-PDFs dedicated to Higgs physics. In each case we also provide the APPLgrid grid name, the range spanned by each distribution, the number of bins , and the kinematical cuts applied to the final-state particles. For associated production with vector bosons, hW and hZ, we impose basic acceptance cuts on the charged leptons from the weak boson decays. All processes have been computed for the LHC 13 TeV

Process	Input cross sections for SM-PDFs for Higgs physics
	Distribution	Grid name	N_bins	Range	Kin. cuts
gg→h	incl xsec	ggh_13tev	1	–	–
	dσ/dpth	ggh_pt_13tev	10	[0,200] GeV	–
	dσ/dy^h	ggh_y_13tev	10	[-2.5,2.5]	–
VBF hjj	incl xsec	vbfh_13tev	1	–	–
	dσ/dpth	vbfh_pt_13tev	5	[0,200] GeV	–
	dσ/dy^h	vbfh_y_13tev	5	[-2.5,2.5]	–
hW	incl xsec	hw_13tev	1	–	p_T(l)≥10 GeV, |η^l|≤2.5
	dσ/dpth	hw_pt_13tev	10	[0,200] GeV	p_T(l)≥10 GeV, |η^l|≤2.5
	dσ/dy^h	hw_y_13tev	10	[-2.5,2.5]	p_T(l)≥10 GeV, |η^l|≤2.5
hZ	incl xsec	hz_13tev	1	–	p_T(l)≥10 GeV, |η^l|≤2.5
	dσ/dpth	hz_pt_13tev	10	[0,200] GeV	p_T(l)≥10 GeV, |η^l|≤2.5
	dσ/dy^h	hz_y_13tev	10	[-2.5,2.5]	p_T(l)≥10 GeV, |η^l|≤2.5
htt¯	incl xsec	httbar_13tev	1	–	–
	dσ/dpth	httbar_pt_13tev	10	[0,200] GeV	–
	dσ/dy^h	httbar_y_13tev	10	[-2.5,2.5]	–

Table 2

Same as Table 1 for the SM-PDFs dedicated to top-quark pair production physics

Process	Input cross sections for SM-PDFs for tt¯ physics
	Distribution	Grid name	N_bins	Range	Kin. cuts
tt¯	incl xsec	ttbar_13tev	1	–	–
	dσ/dptt¯	ttbar_tbarpt_13tev	10	[40,400] GeV	–
	dσ/dy^t¯	ttbar_tbary_13tev	10	[-2.5,2.5]	–
	dσ/dptt	ttbar_tpt_13tev	10	[40,400] GeV	–
	dσ/dy^t	ttbar_ty_13tev	10	[-2.5,2.5]	–
	dσ/dm^tt¯	ttbar_ttbarinvmass_13tev	10	[300,1000]	–
	dσ/dpttt¯	ttbar_ttbarpt_13tev	10	[20,200]	–
	dσ/dy^tt¯	ttbar_ttbary_13tev	12	[-3,3]	–

Table 3

Same as Table 1 for the SM-PDFs dedicated to electroweak gauge boson production physics. The kinematical cuts are applied to the charged leptons from the weak boson decays

Process	Input cross sections for SM-PDFs for electroweak boson production physics
	Distribution	Grid name	N_bins	Range	Kin. cuts
Z	incl xsec	z_13tev	1	–	p_T(l)≥10 GeV, |η^l|≤2.5
	dσ/dptl^-	z_lmpt_13tev	10	[0,200] GeV	p_T(l)≥10 GeV, |η^l|≤2.5
	dσ/dy^l^-	z_lmy_13tev	10	[-2.5,2.5]	p_T(l)≥10 GeV, |η^l|≤2.5
	dσ/dptl^+	z_lppt_13tev	10	[0,200] GeV	p_T(l)≥10 GeV, |η^l|≤2.5
	dσ/dy^l^-	z_lpy_13tev	10	[-2.5,2.5]	p_T(l)≥10 GeV, |η^l|≤2.5
	dσ/dptz	z_zpt_13tev	10	[0,200] GeV	p_T(l)≥10 GeV, |η^l|≤2.5
	dσ/dy^z	z_zy_13tev	5	[-4,4]	p_T(l)≥10 GeV, |η^l|≤2.5
	dσ/dm^ll	z_lplminvmass_13tev	10	[50,130] GeV	p_T(l)≥10 GeV, |η^l|≤2.5
	dσ/dptll	z_lplmpt_13tev	10	[0,200] GeV	p_T(l)≥10 GeV, |η^l|≤2.5
W	incl xsec	w_13tev	1	–	p_T(l)≥10 GeV, |η^l|≤2.5
	dσ/dϕ	w_cphi_13tev	10	[-1,1]	p_T(l)≥10 GeV, |η^l|≤2.5
	dσ/dEtmiss	w_etmiss_13tev	10	[0,200] GeV	p_T(l)≥10 GeV, |η^l|≤2.5
	dσ/dptl	w_lpt_13tev	10	[0,200] GeV	p_T(l)≥10 GeV, |η^l|≤2.5
	dσ/dy^l	w_ly_13tev	10	[-2.5,2.5]	p_T(l)≥10 GeV, |η^l|≤2.5
	dσ/dm_t	w_mt_13tev	10	[0,200] GeV	p_T(l)≥10 GeV, |η^l|≤2.5
	dσ/dptw	w_wpt_13tev	10	[0,200] GeV	p_T(l)≥10 GeV, |η^l|≤2.5
	dσ/dy^w	w_wy_13tev	10	[-4,4]	p_T(l)≥10 GeV, |η^l|≤2.5

Higgs production Total cross sections and rapidity and differential distributions for gluon fusion, vector-boson fusion, associated production with W and Z bosons and associated production with top-quark pairs. No Higgs decays are included, since we are only interested in the production dynamics.

Top quark pair production Total cross section, , and rapidity distributions of the top and the anti-top quarks, and invariant mass , , and rapidity distributions of the system.

Electroweak gauge boson production For Z production: total cross section, , and rapidity distributions of the two charged leptons and of the Z boson, and and invariant mass distribution of the dilepton pair. For W production: total cross section, , and rapidity distributions of the charged lepton and of the W boson, missing and transverse mass distribution. For the W and Z processes, we apply kinematical cuts to the charged leptons from the weak boson decay to reflect the typical acceptance constraints of the LHC experiments.

Clearly, some of these cross sections contain overlapping information, so our list is partially redundant. For instance, if differential distributions are reproduced, this will also be the case for total inclusive cross sections. Similarly, the rapidity distributions of the W and Z bosons are closely related to the rapidity distributions of the leptons from their decay, so including both distributions will lead to a certain degree of redundancy.

This redundancy can be used to provide a non-trivial check of our methodology. For instance, we have verified that by beginning with the total cross sections, only the most extreme bins of the differential distributions, which contribute less to the cross section, might require extra eigenvectors in order to be reproduced to the desired tolerance. Conversely, if we begin the algorithm using differential distributions as input, no additional eigenvectors are required to describe the corresponding total cross sections.

Choice of settings

The SM-PDF method is fully determined by the choice of kinematic region , Eq. (18), which in turn is fully specified by the correlation function and tolerance . The only tunable parameters are thus the scale Q used for the evaluation of correlations in Eq. (16) and the threshold value t. As the choice of the scale Q, we adopt the mean value of the factorization scale at which the PDFs are evaluated by the corresponding APPLgrid grids, that is, the event-by-event weighted average of the value of used in the calculation of each specific cross section or differential distribution.

The only remaining free parameter is then the threshold t, which specifies according to Eq. (18) which points are included in the reduced matrix : low values of t lead to the inclusion of a wider region in phase space, and conversely. Clearly, if is too wide, the reduction will not be very effective and the ensuing number of eigenvectors will be large. On the other hand, if the region is too small, the number of eigenvectors will be small, but it might be lead to a result which is unstable upon small changes of the input observables.

In order to determine a suitable value of t, we use the full set of cross sections listed in Tables 1, 2, and 3. We will henceforth refer to this specific set of observables (and the associate SM-PDF set) as the “ladder”. In Fig. 2 (left) we plot the number of eigenvectors , which we obtained as a function of the parameter t when the SM-PDF methodology is applied to the MC900 prior set, for a fixed tolerance . We show the results for the Higgs, EW and the “ladder” set of input processes.

Fig. 2

Left final number of eigenvectors obtained applying the SM-PDF algorithm to the MC900 NLO PDF set with 900 Monte Carlo replicas, as a function of the threshold parameter t Eq. (18) for fixed tolerance . We show the results for three choices of input cross sections: Higgs (Table 1), electroweak gauge boson production (Table 3), and “ladder” (all processes in Tables 1, 2, and 3). Right correlation Eq. (16) between all the PDFs and the total cross section for Higgs production in gluon fusion, as a function of x (solid blue lines). The value is shown as a dashed line and the region in which the correlation exceeds the threshold is shown as a shaded band

As expected, decreases as the value of t is raised, since in this case fewer points in the grid are selected. While the specific position of the minimum of the curve depends on the input set of cross sections, we see from Fig. 2 that the curve reaches its minimum around for all processes. Note that, as discussed at the end of Sect. 2.3, the value of can fluctuate, typically by one or two units, depending on the specific ordering of the input processes. We therefore choose : this means that we adopt the smallest value of t (i.e. the widest kinematic region) compatible with having the smallest possible number of eigenvectors.

In Fig. 2 (right) we show the value of the correlation coefficient Eq. (16) between the MC900 prior set and the inclusive cross section for Higgs production in gluon fusion, as a function of x and for the seven independent PDF flavors, evaluated at the average scale Q of the grids. The value of the correlation corresponding to is shown as a dashed red line in the plots; the points for which the correlation coefficient (blue curve) is larger in modulus than the threshold are shown as a shaded region.

We observe that, for this specific cross section, the algorithm in the first iteration will include in the region for Eq. (18) only the gluon PDF for , which corresponds to the region that dominates the total cross section for Higgs production in gluon fusion. In Appendix A we provide additional correlation plots, similar to Fig. 2 (right) but for other Higgs production channels, as well as the correlation plots for subsequent iterations, , of the algorithm, illustrating how the selected regions in the grid vary along the iteration.

Results and validation

We now present the results of applying the SM-PDF procedure to the PDF sets and cross sections described in Sect. 3.1. In Table 4 we show the results for the number of eigenvectors obtained, for each input PDF set, using the three different groups of LHC processes that we consider: Higgs, , and W / Z production. In addition, for the Higgs production processes, we have also studied the results of applying our methodology to each of the Higgs production channels individually, as summarized in Table 5. The algorithm has been applied for two different values of the tolerance , namely 5 and 10 %. We also indicate in the bottom row the results for the “ladder” SM-PDF (i.e. including all the above processes.)

Table 4

Number of eigenvectors obtained by applying the SM-PDF procedure, starting from each of the three input prior PDF sets, to the three families of processes summarized in Tables 1, 2, and 3: Higgs production, production, and W / Z production physics. The final row is based on the inclusion of all the three families of processes, in the same order as they are listed. Results are shown for two different values of the tolerance threshold , 5 and 10 %, respectively

Process	N_eig
	MC900		NNPDF3.0		MMHT14
	T_R=5%	T_R=10%	T_R=5%	T_R=10%	T_R=5%	T_R=10%
h	15	11	13	8	8	7
tt¯	4	4	5	4	3	3
W, Z	14	11	13	8	10	9
ladder	17	14	18	11	10	10

Table 5

Same as Table 4, now for the case where the separate Higgs production channels as used as input to the SM-PDF algorithm

Process	N_eig
	MC900		NNPDF3.0		MMHT14
	T_R=5%	T_R=10%	T_R=5%	T_R=10%	T_R=5%	T_R=10%
gg→h	4	5	4	4	3	3
VBF hjj	7	5	10	5	4	3
hW	6	5	6	4	6	3
hZ	11	7	6	4	8	5
htt¯	3	2	4	4	3	2
Total h	15	11	13	8	8	7

Several comments on Table 4 are in order.

Results are reasonably stable upon a change of tolerance, with differences smaller with the MMHT14 prior, which has smaller underlying number of parameters than NNPDF3.0.

The most dramatic reduction in number of eigenvectors is seen for the production of top pairs, or Higgs in gluon fusion, where only eigenvectors are needed. This can be understood as a consequence of the fact that in both cases the dominant contribution to the cross section arises from the gluon distribution in a narrow region of x.

Total cross sections and differential distributions for all the Higgs production modes can be reproduced, in the case of the MC900 prior, with 11 to 15 eigenvectors (depending on the choice of tolerance ).

The number of eigenvectors required is largest for the Higgs and the W / Z family of processes, as one would expect given that in both cases several PDFs in a wide kinematic range are required.

All the processes that we are including can be described with a SM-PDF set, the “ladder”, which includes about the same number of eigenvectors as needed for the Higgs or for the Drell–Yan and W / Z family of processes. This “ladder” SM-PDF, with only 15 eigenvectors, can be used reliably for a large number of LHC cross sections, including those not included in its construction.

Next, in Fig. 3 we show the total number of eigenvectors , which are required, for a tolerance of , as more and more processes are sequentially included, until the complete list of processes in Tables 1, 2, and 3 has been exhausted. This plot demonstrates the robustness and flexibility of the SM-PDF algorithm, in that it shows how more processes can be added without information loss to a reduced PDF set, thereby allowing for a study of the information brought in by each process. In Fig. 3 results are presented for the three input PDF sets, MC900, NNPDF3.0 and MMHT14. As already seen in Table 4, a smaller number of eigenvectors is required in order to describe the MMHT14 set, which has a smaller underlying number of parameters than the NNPDF3.0 set; the combined MC900 set requires roughly the same number of eigenvectors as NNPDF3.0, which is contained in it. Inspection of Fig. 3 indicates which processes bring in new information in comparison to those already included. For instance, the fact that the number of eigenvectors is unchanged when adding all the observables related to top-quark pair production shows that SM-PDFs based on Higgs processes also describe top production.

Fig. 3

Total number of eigenvectors required by the SM-PDF algorithm to describe a sequentially increasing number of input cross sections and distributions, for a tolerance . Results are presented for the three prior PDF sets, namely MC900, NNPDF3.0, and MMHT14

Fig. 4

Upper plots comparisons of the predictions for the (left) and rapidity (right) differential distributions in Higgs production in gluon fusion between the prior MC900 and the corresponding Higgs SM-PDFs for two different values of the tolerance , 5 and 10 %. Results are shown normalized to the central value of MC900. Lower plots same comparison, now for the SM-PDFs, showing the invariant mass of the pair (left) and the top-quark rapidity (right). See Tables 1 and 2 for the details of the binning and the kinematical cuts in each case

In Figs. 4, 5, 6 we compare various cross sections and differential distributions computed with the MC900 prior PDF set and with the corresponding SM-PDFs for some of the cases discussed above, normalized to the central value of the prior. In the upper plots of Fig. 4, we show the Higgs and y distributions in gluon fusion production, comparing with the Higgs SM-PDF. In the lower plots of Fig. 4, we show the top-quark pair invariant mass and top rapidity distributions, comparing with the SM-PDF. In Fig. 5 we compare various differential distributions in weak gauge boson production with the W, Z SM-PDFs, and finally in Fig. 6 we compare the “ladder” SM-PDFs with various total inclusive cross sections.

Fig. 5

Same as Fig. 4 for representative differential distributions in W and Z production, comparing the MC900 prior with the W, Z SM-PDFs

Fig. 6

Same as Fig. 4 for the “ladder” SM-PDF, now comparing with the total ggH, , Z, and W inclusive cross sections

In these comparisons, results are shown for two values of the tolerance, and . PDF uncertainties are shown as one-sigma confidence intervals; for the MC900 prior, the central 68 % confidence intervals are also shown (inner ticks). In all cases we observe excellent agreement between the prior and the corresponding SM-PDF sets, which provides a further validation of the reliability of the method.

We have also verified that SM-PDFs reproduce well PDF correlations, even though the tolerance criterion Eq. (15) is only imposed on diagonal PDF uncertainties. The PDF-induced correlation between two cross sections computed using a Monte Carlo PDF set is given bywhile for a Hessian set it isIn Fig. 7 we show the difference between the correlations determined using the MC900 prior (from Eq. (32)) and the “ladder” SM-PDF set (from Eq. (33)), with , for all the total inclusive cross sections used as input to the “ladder” SM-PDF set. We find that the deviation in correlation is at the few percent level or better for most cases, and anyway never worse than 20 %.

Fig. 7

Differences in the correlation coefficients between the MC900 prior and the “ladder” SM-PDFs with , computed for all the inclusive cross sections that enter the construction of the latter

An additional validation test can be performed by comparing the predictions for a given SM-PDF outside the kinematic range of the input processes. To illustrate this point, in Fig. 8 we compare the and rapidity distributions in Higgs production via gluon fusion using the Higgs SM-PDF (which uses as input the processes in Table 1) but now with an extended kinematical range: the rapidity distribution now includes , rather than the range used as input, and the distribution covers now GeV as compared to the original input GeV. In both cases, we show both the standard deviation (left) and the full probability distribution obtained with the prior and the two compressed sets with and ; the smoothened probability distributions are obtained using the using the Kernel Density Estimation (KDE) method discussed in Ref. [12]. The good agreement seen in all cases demonstrates the robustness of the SM-PDF method: namely, SM-PDF sets are stable upon variations of kinematic cuts and binning of the input cross sections.

Fig. 8

The and rapidity distributions for Higgs production in gluon fusion, computed with the MC900 prior and with the Higgs SM-PDFs, for two values of the tolerance , this time in a kinematic range that doubles that of the input processes in Table 1 (see text). In the left plot we show the standard deviation in each bin, while in the right plot we show the full probability distributions per bin, reconstructed using the Kernel Density Estimate (KDE) method

While the SM-PDFs are stable upon extrapolation, they will not provide accurate predictions when used for processes dominated by PDFs in an altogether different kinematic range. To illustrate this point, in Fig. 9 we show predictions for inclusive jet distributions obtained using the Higgs and ladder SM-PDF sets, compared to the result obtained using the MC900 prior. Specifically, we show the distributions in the most forward rapidity bin () of the ATLAS 2010 inclusive jet measurement [23]; bins are ordered in increasing . Clearly, the agreement deteriorates at large , where results depend on the large-x quarks and gluon, which are weakly correlated to the processes included in the construction of both the Higgs and the “ladder” SM-PDF sets. This also suggests that good agreement, with a marginally larger number of eigenvectors, could be likely obtained by just widening the range of some of the inputs to the “ladder”, such as, for instance, including the Higgs transverse-momentum distribution up higher values of . In fact, we have explicitly checked [24] that the “ladder” PDF set provides comparable accuracy to the PDF4LHC15 30 eigenvector set when used for the determination of all the hadronic observables included in the NNPDF3.0 PDF determination [6], despite having almost half the number of eigenvectors.

Fig. 9

Same as Fig. 8, but now for the ATLAS inclusive jet distribution in the forward region, and using the “ladder” SM-PDF set

A posteriori combination of SM-PDFs

So far, we considered the construction of a PDF set tailored to a given list of input cross sections. However, one may also encounter the situation in which two SM-PDF sets constructed using different processes as input are already available, and wishes to use them simultaneously, without having to produce a new dedicated SM-PDF set using as input the two processes at the same time. A typical application is a computation in which one of these processes is the signal and the other to a background. The SM-PDF methodology also allows one to deal with this situation: we first discuss how this is done, and then we present an example of an application.

General method

In Sect. 2 we have shown how, starting from a Monte Carlo PDF prior, , Eq. (1), we can construct a specialized minimal Hessian representation, , Eq. (28), in terms of a reasonably small number of eigenvectors. The result of the SM-PDF algorithm can be expressed as a regular Hessian PDF set, with the error parameters given by Eq. (29). Alternatively, one can directly use the final matrix of Hessian coefficients P to express the cross sections computed with each of the replicas of the prior set, Eq. (8), as linear combinations of cross sections computed with the final eigenvector sets, Eq. (13). The two results are equivalent by linear error propagation.

However, we can also read Eq. (13) in reverse: if we definewe can view the set of differences (for each of the observables ) as a Monte Carlo set of cross sections, containing the same information as the reduced SM-PDF set. In other words, the valuesof the observable can be viewed as “pseudo-Monte Carlo” replicas, to be used to compute uncertainties and correlations using the standard Monte Carlo procedure.

If two sets of SM-PDFs corresponding to different processes are available, we can then combine the information contained in them by first turning the predictions obtained from them into replicas using Eq. (35), and then viewing the set of Monte Carlo replica predictions obtained in each case as our best approximation to the Monte Carlo set of predictions for that process obtained with the original PDF replica set. These sets of prediction replicas can then be used in order to compute any quantity which depends on both processes using the standard Monte Carlo methodology, by just making sure that each process is computed using its corresponding replicas.

Validation

We illustrate and validate the methodology presented in Sect. 4.1 with an example. We use as input prior the NNPDF3.0 NLO set with replicas and then generate two SM-PDFs for a fixed choice of the tolerance . The first SM-PDF takes as input the processes from Table 2, while the second is constructed from the W, Z processes of Table 3.

We now use these two SM-PDF sets to calculate the PDF uncertainties on the and the W total inclusive cross sections. This can be done both with the original representation, Eq. (7), or with the new SM-PDF Hessian representation. As shown in Table 4, we find for the SM-PDF and for the W, Z SM-PDF. We obtain the following results for the total cross sections: for the cross section with SM-PDFs and for the W cross section with W, Z SM-PDF Now suppose that we want to compute a quantity which depends both on the and the W cross sections, such as the ratio between the two, . In the computation of the PDF uncertainty on this ratio, it is essential to properly account for the cross-correlations between the two processes. This can be achieved by recasting the results of the two different SM-PDFs into corresponding Monte Carlo sets of predictions through Eq. (35).

Namely, the PDF uncertainty on the cross-section ratio is given bywhere and have been obtained using Eq. (35) with the P matrix that corresponds, respectively, to the and W, Z SM-PDF sets.

Using Eq. (40) we getto be compared to the result obtained from the NNPDF3.0 prior, using the original replicas,which is identical for all practical purposes.

It is important to realize that while Eq. (42) requires the calculation of cross sections, Eq. (41) only requires the knowledge of the cross section differences for the two observables, which is equal to the sum of the number of eigenvectors in the two sets which are being combined; in our case, , with great computational advantage.

As a further cross-check, we have recomputed the same cross section ratio by using the methodology of Sect. 2, namely, by constructing a dedicated SM-PDF set using as input the two families of processes, and W, Z, simultaneously.

This new SM-PDF has now 17 eigenvectors for the case of a tolerance and leads toThis shows that the advantage of constructing a dedicated set in comparison to combining the pre-existing sets is marginal, as the accuracy is the same, and the total number of eigenvectors has only decreased by one unit.

Delivery

Building upon our previous MC2H methodology for the construction of reduced Hessian representations of PDF uncertainties [11], we have presented an algorithm for the construction of a minimal Hessian representation of any given prior PDF set, specialized to reproduce a number of input cross sections. We have shown that the algorithm can be used to construct specialized minimal PDF sets which reproduce with percent accuracy the central values and PDF uncertainties for all input observables in terms of a substantially smaller number of eigenvectors as compared to the prior PDF set. A remarkable advantage of the SM-PDF methodology is that the complete information contained in the original prior set is kept at all stages of the procedure. As a consequence, it is possible to add new processes to any given SM-PDF set with no information loss. Also, it is possible to combine a posteriori SM-PDF sets corresponding to different processes without any new computation.

The SM-PDF code is publicly available from the repository

https://github.com/scarrazza/smpdf/

The code is written in Python using the numerical implementations provided by the NumPy package. Customized interfaces to APPLgrid and LHAPDF6 are also included. The package also includes the APPLgrid grids for all the processes listed in Tables 1, 2, and 3, and additional processes can be easily generated upon request.

The input of the SM-PDF code is the prior PDF set and the list of cross sections to be reproduced. The code settings can be modified by the user by means of a steering card. The cross sections can be provided either by indicating the name of the APPLgrid or by means of a text file (for predictions computed with external codes). An example steering card for the code is presented in Appendix B.

The output of the code is then the corresponding SM-PDF set, directly in the LHAPDF6 format, as well as the corresponding direct and inverse Hessian parameter matrices, P and , respectively as a CSV file. These rotation matrices allow one to easily transform the computed cross sections back and forth from any SM-PDF representation to the prior representation, as well as transforming between different SM-PDF representations, as explained in Sect. 4.

Together with this, a number of additional validation features are included in the SM-PDF package. In particular, comparisons at the level of the input cross sections as those presented in Figs. 2, 4, and 7 can be generated automatically by filling the appropriate options in the YAML configuration file, without the need of writing additional code. The user is encouraged to refer to the documentation for a more extensive description of the different features available. In addition, a web interface to similar to that of APFEL Web on-line PDF plotter [25, 26] is currently under consideration.

Finally, the SM-PDFs constructed in Sect. 3 are also available from the same webpage in the LHAPDF6 format. Users can produce the SM-PDFs that more suitable for specific applications by generating the suitable cross section theory calculations and then running the SM-PDF code. However, users are encouraged to contact the authors for support if assistance is needed. Additional SM-PDFs can be added to this webpage upon request.