aMCfast: automation of fast NLO computations for PDF fits

We present the interface between MadGraph5_aMC@NLO, a self-contained program that calculates cross sections up to next-to-leading order accuracy in an automated manner, and APPLgrid, a code that parametrises such cross sections in the form of look-up tables which can be used for the fast computations needed in the context of PDF fits. The main characteristic of this interface, which we dub aMCfast, is its being fully automated as well, which removes the need to extract manually the process-specific information for additional physics processes, as is the case with other matrix element calculators, and renders it straightforward to include any new process in the PDF fits. We demonstrate this by studying several cases which are easily measured at the LHC, have a good constraining power on PDFs, and some of which were previously unavailable in the form of a fast interface.


Introduction
The accurate determination of the Parton Distribution Functions (PDFs) of the proton [1][2][3][4][5] is one of the most important tasks for precision phenomenology at the Large Hadron Collider (LHC). PDFs are a dominant source of theoretical uncertainty in the predictions for Higgs boson production, where the errors that affect them degrade the accuracy of the Higgs characterization in terms of couplings and branching fractions [6]; they induce large uncertainties in the cross sections for processes with very massive new-physics particles [7]; and they substantially affect Standard Model (SM) precision measurements such as those of the mass of the W boson [8] and of the effective lepton mixing angle sin 2 θ l eff [9]. For these reasons, an active program towards better PDFs is being carried out by different groups [10][11][12][13][14], which emphasise the use of new experimental inputs, more accurate theoretical calculations, and improved fitting methodology.
Modern global PDF analyses are based on a variety of hard-scattering data such as deep-inelastic scattering (DIS) structure functions at fixed-target experiments, lepton-proton cross sections from the HERA collider, and inclusive W , Z, and jet production at hadron colliders. Since the beginning of the LHC data taking, these and many other processes for which measurements have become available and can be used to constrain PDFs. For example, LHC data that might provide information on PDFs encompass inclusive electroweak vector boson production [15,16], inclusive jet and dijet production [17][18][19][20][21], direct photon production [22], top quark pair production cross sections [23], W production in association with charm quarks [16,24], low and high-mass Drell-Yan production [25], the W and Z bosons large-p T distributions and their ratios [26], high-mass W production, and single-top production, as well as ratios of cross sections measured at different center-of-mass energies [20,27].
It is thus clear that a wide variety of high-quality measurements sensitive to PDFs are already available; more data will follow in the coming years. Therefore, in order to improve the accuracy of PDF determinations, it is essential to include in PDF fits as many of these measurements as possible, in order to constrain different combinations of PDFs in a wide range of Bjorken x's. The most serious difficulty in doing so is due to the fact that next-to-leading order (NLO) QCD calculations of hadron collider processes with realistic acceptance cuts are much slower than what is needed in the iterative PDF fitting procedure, which requires computing the theoretical predictions a very large number of times. In order to bypass this problem, a popular solution in the past has been that of performing leading order (LO) computations, supplemented by bin-by-bin K factors. Unfortunately, such a solution is not sufficiently accurate to match the precision of present and future LHC data; in particular, it has the very undesirable feature of neglecting the combinations of initial-state partons that do not appear at the LO.
In order to overcome this problem, several solutions have been proposed. The underlying idea common to all these approaches is: interpolating the PDFs in the x, Q 2 -plane with some suitable polynomial basis; precomputing the hadronic cross section by using the basis members as input; and finally reconstructing the original calculation with the numerical convolution of the precomputed cross sections and the actual PDFs. Note that the information about the latter is only required at a finite number of points (x i , Q 2 j ), which are called the interpolating-grid nodes. Therefore, the time-consuming task of the precomputation of the cross section with basis members is performed only once, and the reconstruction of the full result associated with arbitrary PDFs is extremely fast.
The strategy sketched above, which is closely related to the one adopted by x-space PDF evolution codes [28][29][30], is what underpins the two best-known fast interpolators of NLO QCD cross sections, FastNLO [31,32] and APPLgrid [33]. FastNLO is interfaced to the jet-production code NLOjet++ [34], and is thus able to provide fast computations of multijet production at lepton-proton and hadron-hadron colliders 1 . APPLgrid is interfaced to various programs, including NLOjet++ and MCFM [36]. The main drawback of these tools, namely that of being capable of handling a relatively small number of processes, is a consequence of the fact that adding new ones requires an ad-hoc procedure, which is rather time-consuming and error-prone. This is also the principal reason why they are only interfaced to calculations which are of NLO (and, in the one case mentioned above, NNLO) in QCD, but neither to (N)LO results matched to parton showers, nor to NLO results in the electroweak theory.
The goal of this work is that of solving all of these problems in a general manner, which is possible thanks to the fact that NLO(+PS) calculations can now be routinely done by means of automated codes. In fact, among the many features of such codes there are two which are directly relevant to the problems at hand: firstly, given a sufficient amount of CPU power the cross section for any process, however complicated, can be computed; and secondly, the way in which these cross sections are handled and the form in which they are written are completely standardised. This is what renders it possible to construct a generic and automated interface between an automated cross section calculator and a fast interpolator. The main result of this paper is the construction of such an automated interface, that we call aMCfast, and which bridges the automated cross section calculator MadGraph5 aMC@NLO [37] with the fast interpolator APPLgrid. Thus, the chain Mad-Graph5 aMC@NLO -aMCfast -APPLgrid will allow one to include, in a straighforward manner, any present or future LHC measurement in an NLO global PDF analysis. We point out that a strategy analogous to that pursued in the present paper has motivated the recent construction of MCgrid [38], whereby APPLgrid has been interfaced to Rivet [39].
We remind the reader that MadGraph5 aMC@NLO contains all ingredients relevant to the computations of LO and NLO cross sections, with or without matching to parton showers. NLO results not matched to parton showers (called fNLO [37]) are obtained by adopting the FKS method [40,41] for the subtraction of the singularities of the real-emission matrix elements (automated in the module MadFKS [42]), and the OPP integral-reduction procedure [43] for the computation of the one-loop matrix elements (automated in the module MadLoop [44], which makes use of CutTools [45] and of an in-house implementation of the optimisations proposed in ref. [46] (OpenLoops)). Matching with parton showers is achieved by means of the MC@NLO formalism [47]. In the present public version, MadGraph5 aMC@NLO is restricted to computing NLO QCD corrections to SM processes; however, as discussed in ref. [37], all obstacles that enforce such a limitation have been cleared, paving the way to higher-order calculations in the context of arbitrary renormalisable theories in the near future. As far as APPLgrid is concerned, only a few high-level routines have been extended in view of its interface with aMCfast. All of these modifications are of a technical character, except one which is related to the computation of the factorisation and renormalisation scale dependences, as we shall explain in more detail later.
The scope of this paper is that of fNLO QCD computations. Alternatively, in aMCfast the specific nature of the perturbative expansion is used only in a rather trivial way, since it determines a number of interpolating grids and their linear combination that defines the physical cross sections. Therefore, when e.g. electroweak corrections will become publicly available in MadGraph5 aMC@NLO, that feature will be immediately inherited by aMCfast through some straighforward generalisation. Furthermore, what is done here will allow one, with only a few minimal extensions, to construct fast interfaces to NLO+PS predictions. This is expected to have several beneficial effects in the context of PDF fits: a closer connection to the experimentally accessible observables, a wider range of data that 1 Recently, the FastNLO interface has been generalised to processes other than jet production, for instance to the approximate NNLO calculation of differential distributions in top-pair production of ref. [35].
can be used to constrain PDFs, and the possibility of eventually extract PDFs that are specifically tailored to their use with NLO event generators. This paper is organised as follows: in sect. 2 we give a short introduction to the interpolatinggrid techniques employed here, review the computation of cross sections in MadGraph5 aMC@NLO, and discuss how the separation between PDFs, partonic cross sections, and α S dependence can be exploited to construct a fast interface using the APPLgrid routines. The flexibility of aMCfast is illustrated in sect. 3, where we present results for many relevant LHC processes, some of which are obtained for the first time with a fast NLO interface, and whose numerical performance and accuracy are analysed. Finally, in sect. 4 we summarise our findings and briefly discuss the plans for the future developments of aMCfast. An appendix reports some additional information.

Automation of fast NLO computations
In this section we first outline the basics of an interpolation technique based on the expansion of a given function onto a basis of polynomials; we then discuss how short-distance cross sections can be represented, in the most general manner, in terms of interpolating grids, that allow them to be quickly computed with arbitrary PDFs, and factorisation and renormalisation scales. Finally, we show how these formulae can be employed in the construction of the aMCfast bridge that interfaces APPLgrid to MadGraph5 aMC@NLO.

The construction of interpolating grids
The basic idea used by APPLgrid is that of a Lagrange-polynomial expansion. Given a function F (z) of a real variable z, one has the representation: where s is a given integer (the interpolation order), δ is a small number (the grid spacing or binning; pδ for some integer p is called a grid node), the interpolating functions are: and we have denoted by ⌊u⌋ the largest integer smaller than or equal to u: The equality in eq. (1) holds up to functional terms of order I (s+1) i and higher. Such terms vanish identically when the argument of F coincides with a grid node, so that eq. (1) is an identity in that case; this is straightforward to prove, and follows directly from the values that the interpolating functions take when computed with an integer argument: When z is not a grid node (i.e., z = pδ for any integer p), eq. (1) tells one that F (z) is reconstructed by using the values that F takes in the (s + 1) grid nodes which are nearest to z; the number of relevant nodes to the left of z is equal to number of nodes to the right of z, possibly up to one.
For any given function S(z), let us now compute the simple example integral by means of its corresponding Riemann sums (or, equivalently, by Monte Carlo methods). This implies with M points z k ∈ (a, b), and Φ k suitable normalisation factors. By using eq. (1), we obtain where we have defined: which is the integer associated with the leftmost grid node in the set of the (s + 1) nearest neighbours of z. By means of a change of the summation variable i, eq. (7) becomes: with Equation (11) defines the grid values G j . Owing to the Θ function it contains, the sum in eq. (10) features a finite number of non-null contributions (if the range (a, b) is finite). Thus, the meaning of eq. (10) is that the integral J can be computed a posteriori by knowing a finite number of grid values, and the values of the function F at the grid nodes; more importantly, this a-posteriori computation can be done for any function F , because the grid values are independent of F , and can therefore be pre-evaluated and stored for a given function S. This also explains the reason for writing the integrand of eq. (5) as a product of two functions: this is convenient whenever the time spent in evaluating F (the "fast" function) and S (the "slow" function) is small and large respectively. When this is the case, the computation of the grid {G j } may be time-consuming, but it is an operation that has to be carried out only once; on the other hand, the subsequent evaluations of eq. (10) will be quick. We also point out that the derivation above is unchanged in the case where z is not the integration variable, but a function itself of one (or more) integration variable(s) for the problem at hand. This is because the starting point is actually eq. (6), and not eq. (5), and in the former the role of z k as integration variable can be fully ignored.

Generalities on short-distance cross sections
In what follows we use the expressions derived within the FKS subtraction formalism. The same notations as in ref. [48] are adopted; this is particularly convenient in view of the fact that, in that paper, cross sections were represented in terms of PDF-and scale-independent coefficients, and these will be the main ingredients in the definition of the interpolating grids. We point out, however, that the procedure outlined below remains valid regardless of the subtraction method chosen to perform the NLO computations. Lest we clutter the notation with details which are irrelevant here, we work by fixing the partonic process. This implies, in particular, that we do not need to write explicitly the identities of the incoming partons; we shall reinstate them later. The NLO short-distance cross section relevant to a 2 → n + 1 process consists of four terms: called event (α = E), and soft, collinear, and soft-collinear (α = S, C, SC) counterevents, respectively. The quantities dχ Bj and dχ n+1 are the integration measures over the Bjorken x's and the (3n − 1) phase-space variables respectively, while f 1 and f 2 are the PDFs relevant to the colliding partons coming from the left and from the right. The W (α) 's can be parametrised as follows: where the coefficients W are (renormalisation and factorisation) scale-and PDF-independent; the last term on the r.h.s. of eq. (14) is the Born contribution which, as the notation suggests, factorises α b S . In eqs. (13) and (14) we have denoted by x the Bjorken x's, factorisation scale, and renormalisation scale respectively; in general, they may assume different values in the event and various counterevent configurations. Finally, Q is the Ellis-Sexton scale 2 . The integration of the above cross section leads to a set of 4N weighted events: with K (α) n+1;k an (n+1)-body kinematic configuration (possibly degenerate, in which case it is effectively an n-body configuration that can be used to compute Born matrix elements), and the event weights defined as follows: where we understand possible normalization pre-factors (such as 1/N , or adaptive-integration jacobians). Given an observable O, its h th histogram bin defined by O UPP will assume, at the end of the run, the value: In view of eq. (14), it is convenient to recast eq. (17) as follows: where, using eqs. (13) and (16), one has: Here we have defined: namely, the values of the short-distance coefficients and of the scales computed at the kinematic configurations associated with each of the events and counterevents of eq. (15). We now apply the method outlined in sect. 2.1 to eqs. (20)- (23). In order to do so, we recall that the main result of that section is that of computing the number J, whose original expression is that of eq. (6) (regardless of the fact that such an expression was in turn obtained from the integral of eq. (5), for the reason explained at the end of sect. 2.1), by means of the interpolating grid given in eq. (11) and of eq. (10).
One starts by observing that the h th bin value σ O,β with an interpolating grid, as is done for J in eq. (10). In order to do so, one may easily proceed by analogy. First of all, since we are interested in a quick recomputation of the cross section after changing the scales and the PDFs, and given that the latter depend on the Bjorken x's and the factorisation scale, the implication is that the role of the variable z of eqs. (6)-(11) needs to be played by the four variables: In other words, we have the correspondence: The polynomial expansion of eq. (1), whence the interpolating grid is derived, is valid for each of the variables in eq. (28); therefore, the only change w.r.t. the case of the computation of J is the fact that the one-dimensional grid {G j } defined in eq. (11) will be replaced by a four-dimensional grid, corresponding to the variables in eq. (28). For the construction of the latter, it is sufficient to compare directly eq. (6) with eqs. (20)- (23). The role of the "slow" and "fast" functions S and F will be played by the short-distance coefficients W and by the scale-and PDF-dependent terms respectively; that of the factor Φ k by the bin-defining Θ (h,α) k of eq. (18). In other words, we have the following identifications: and It is now sufficient to use eqs. (30)-(34) to rewrite eqs. (10) and (11). We denote by j 1 , j 2 , j 3 , and j 4 the indices that run over the grid nodes relevant to the variables x 1 , x 2 , µ F , and µ R , respectively; δ i and s i , i = 1, . . . 4, are the corresponding grid spacings and interpolating orders. We obtain: with the interpolating grids: Equations (35)- (39) give the most general representation of NLO cross sections in terms of interpolating grids. We shall employ them (in a simplified form) in the next subsection, in the construction of the aMCfast bridge to APPLgrid.

The interface of APPLgrid with MadGraph5 aMC@NLO
APPLgrid is a general-purpose C++ library that provides a suitable number of interpolation and convolution routines that can be used to construct fast interfaces to NLO calculations of lepton-proton and hadron-hadron collider processes. APPLgrid is widely used by various PDF fitting collaborations as well as by ATLAS and CMS in their own PDF studies. Our purpose is that of exploiting APPLgrid for the construction of the interpolating grids of eq. (39), and their subsequent use in the calculation of the histogram bins, eqs. (35)- (38). As those equations show, APPLgrid needs to be given, in an initialisation phase, the observables to be computed and their respective binnings (we stress again that each bin of each observable corresponds to a set of four grids per type of parton luminosity); and, on an event-by-event basis, the values of those observables and the short-distance coefficients W . These tasks are essentially what aMCfast is responsible for: it extracts the relevant information from MadGraph5 aMC@NLO, and feeds them to APPLgrid, in the format required by the latter. The first observation is that, in its current version, APPLgrid supports three-dimensional grids, while those in eq. (39) are four-dimensional. The reason for the latter is that in our derivation we have left the freedom of choosing different functional forms for the factorisation and renormalisation scales. On the other hand, such a flexibility is seldom exploited, and typically one chooses µ F and µ R equal in the whole phase space, up to an overall constant. This is equivalent to setting: with µ a function of the kinematics. When doing this, two of the four variables in eq. (28) are degenerate, and therefore one must simply consider: which reduces the number of grid dimensions from four to three. In this case, it is also convenient to set 3 : so that in eqs. (36) and (37) respectively. The other changes to eqs. (35)-(39) due to the reduction of the grid dimensionality are all trivial: one eliminates the sums over j 4 , formally replaces: and µ F with µ in the next-to-last line on the r.h.s. of eq. (39) (this and eq. (44) are due to the fact that now it is the variable µ which corresponds to one of the dimensions of the grid), and finally eliminates the last line of that equation. The next thing to consider is that APPLgrid does not use the variables in eq. (41) directly, but rather constructs the grids using: which are defined through a change of variables, several of which are available but is usually set to: This is because the grids feature a linear spacing, and given the PDFs and α S dependences on Bjorken x's and scales a linear sampling in terms of the variables of eq. (45) turns out to be more effective than one based the variables of eq. (41). In eq. (46) the parameter κ controls the relative density of grid nodes in the large-w.r.t. the small-x region; in eq. (47), the parameter Λ is typically chosen to be of the order of Λ QCD . Finally, one needs to consider the fact that the formulae presented above are obtained for a given partonic process. In order to compute the actual hadronic cross section, one needs to sum over all possible such processes, so that eq. (13) must be generalised as follows: where r and s are the identities of the incoming partons, and q collectively denotes the identities of all outgoing partons. Note that, for simple-enough cases, the sum over q is trivial, since (r, s) are sufficient to determine unambiguously a partonic process; however, the notation used in eq. (48) is general, and encompasses all possible situations. Since the interpolating grids are additive, the most straightforward solution is that of following the procedure outlined so far, and of creating a grid for each possible partonic process. There is however a better strategy, that helps save disk space and decrease memory footprint, in that it reduces the number of interpolating grids. This is based on the observation that one may find pairs of parton indices (r, s) and (r ′ , s ′ ) such that: This suggests to rewrite eq. (48) in the following way: where and the values of T (l) rs are either zero or one; we have implicitly defined W rsq for (r, s, q) and l such that T (l) rs = 1. We point out that, while there may be more than one way to write the r.h.s. of eq. (50) (in other words, the luminosity factors F (l) may not be uniquely defined), it is always possible to arrive at such a form, thanks to the fact that the r.h.s.'s of eqs. (48) and (50) are strictly identical. In fact, eq. (50) is what is used internally by MadGraph5 aMC@NLO; one of the tasks of aMCfast is that of gathering this piece of information, and of using it to construct the luminosity factors F (l) to be employed at a later stage. In eq. (50) one factors out identical matrix elements; since the computation of these is the most time-consuming operation, this procedure helps increase the overall efficiency, which is larger the larger the number of terms in each luminosity factor F (l) . The fact that, in general, the number of terms in the sum over l in eq. (50) is smaller than that in the sums over (r, s) in eq. (48) is ultimately what allows one to reduce the number of interpolating grids. The counting of the terms in such sums is easy when (r, s) are sufficient to determine uniquely the partonic process. Denoting by n l the largest value assumed by the luminosity index l for a given process: and by N F the number of active light flavours, one has: with the two limiting cases being either that where all of the allowed PDF combinations are associated with the same partonic matrix element, or that where each PDF combination corresponds to a different partonic matrix element. In appendix A we shall give the explicit form of eq. (51) for all of the processes studied. By putting everything together, we finally arrive at the representation of the h th bin value and of its corresponding grids as they are constructed by the MadGraph5 aMC@NLO -aMCfast -APPLgrid chain. We denote by δ y and s y the grid spacing and interpolating order relevant to the variables y 1 and y 2 ; the analogous quantities relevant to the variable τ are denoted by δ τ and s τ respectively. We have: with the grids: × I (sy) In eq. (58) we have inserted the partonic indices l and q in the short-distance quantities W , and we have used eqs. (46) and (47) to introduce: and we have defined (note the factors ξ F and ξ R ): Equations (54)-(58) are the main results of this paper. In the initialisation phase, aMCfast provides APPLgrid with the total number of grids needed (equal to the sum over all observables of the number of bins relevant to each observable, times four), the grid spacings δ y and δ τ , the interpolation orders s y and s τ , and the interpolation ranges in y i and τ . These information are under the user's control; in particular, one must make sure that the latter ranges are sufficiently wide for the process under consideration, and one will want to be careful in the case of a dynamical scale choice for µ. Then, during the course of the run and event-by-event, aMCfast gets Θ i;k , and τ (α) k from MadGraph5 aMC@NLO and feeds 4 them to APPLgrid, whose grid-filling internal routines iteratively construct G (h,β,l) j 1 j 2 j 3 as defined in eq. (58). We conclude this section with two remarks. Firstly, the sums over l and q in the formulae above achieve the sum over subprocesses which is necessary in order to obtain the hadronic cross section. In practice, in MadGraph5 aMC@NLO the number of contributions which are integrated separately and eventually summed to give the physical cross section is often larger than the number of subprocesses, owing to the FKS dynamic partition and to the use of multi-channel techniques (see ref. [37] for more details). When interfacing MadGraph5 aMC@NLO with APPLgrid, each of these contributions generates temporary grids, which are then combined by aMCfast at the end of the run to give the actual grids of eq. (58) that are to be used in fast computations (i.e., in eqs. (54)-(57)). Secondly, thanks to the fact that the information on the cross section is (also) given in terms of the scale-and PDF-independent coefficients W , MadGraph5 aMC@NLO is capable of computing a-posteriori PDF and scale uncertainties independently of its interface to APPLgrid, by means of a reweighting procedure (see ref. [48]). While this feature renders the recomputation of the cross section with alternative PDFs and/or scales much faster than the original calculation, it is still not fast enough to be employed in PDF fits, for which the only viable solution is that of using the interpolating grids discussed here. The reason is related to the sum over events (1 ≤ k ≤ N and α = E, S, C, SC in the previous formulae): while such sum is performed only once in the case of interpolating grids (when the grid is constructed: see eq. (58)), it must be carried out for each new choice of PDFs and scales in the context of reweighting. This difference is crucial, because N must be a large number (especially so in fNLO computations), in order to obtain a good statistical precision 5 .

Scale choices
As one can see from eq. (14), the grids G (h,0,l) j 1 j 2 j 3 and G (h,B,l) j 1 j 2 j 3 are all that is needed when one is interested in computing the cross section that corresponds to setting the factorisation and renormalisation scales equal to their reference value µ, i.e. with ξ F = ξ R = 1 (see eq. (40)). When ξ F = 1 and/or ξ R = 1, then the grids G (h,F,l) j 1 j 2 j 3 and/or G (h,R,l) j 1 j 2 j 3 , respectively, are also necessary (see eqs. (55) and (56)). We point out that these two grids are not constructed when APPLgrid is interfaced to codes other than MadGraph5 aMC@NLO. In that case, APPLgrid is still capable of computing the cross section for non-central scales, since the form of the latter can be determined by using renormalisation group equations (RGEs). By doing so, one arrives at an expression (e.g., eq. (14) of ref. [33]) which is in one-to-one correspondence with the results derived here (in order to see this, one must use the explicit expressions of the W coefficients given in ref. [48]); this is not surprising, since both are ultimately a direct consequence of RGE invariance.
Although fully equivalent from a physics viewpoint, the two-grid and four-grid approaches do not give pointwise-identical results (in other words, their outcomes are strictly equal only in the limit of infinite statistics). The main advantage of using only two grids is that of a smaller memory footprint. Conversely, the two-grid approach has two drawbacks that we consider significant, and are the reason why aMCfast works with four grids. Firstly, the ξ F -dependent term features a convolution (i.e., a one-dimensional integral) of a one-loop Altarelli-Parisi kernel [50] with a PDF, P (0) ⊗ f . This convolution is effectively performed when summing over events in the four-grid approach, in eq. (55). The absence of G (h,F,l) j 1 j 2 j 3 in the two-grid procedure implies that APPLgrid needs to perform this convolution explicitly, and it does so by resorting to an external code, the PDF-evolution package HOPPET [28]. APPLgrid does work without HOPPET being installed, but in this case it is not possible to choose a non-central factorisation scale. Secondly, and related to the previous point: the representation of a cross section in terms of interpolating grids is exact up to terms of higher order in the interpolation-order parameter; this implies, for example, that eq. (55) is an identity up to terms that contain the interpolating functions I (sy+1) and I (sτ +1) -as we shall show later, in all practical cases such missing terms are totally negligible. The crucial point here is that this is an identity regardless of the statistics used in the computation of the cross section (i.e., it is independent of N , up to fluctuations that may affect the grid construction in the case of very small N 's; again, we shall document this later). However, and for the specific case of the ξ F -dependent eq. (55), this is true only because the convolution integral P (0) ⊗ f implicit in there is carried out by the very procedure (i.e., the same number of events, the same seeds) that fills the interpolating grid. Any other way of computing P (0) ⊗ f , and in particular one that uses an external program such as HOPPET, implies that this property does not hold. Therefore, in the two-grid approach the difference between a cross section computed with non-central scales, and its representation in term of grids, is much larger than formally guaranteed by the chosen interpolation orders; this difference tends to zero only in the limit of infinite statistics N → ∞. We did explicitly verify that this is indeed the case.
In summary, we believe that working with four grids in APPLgrid gives a couple of clear advantages over the two-grid procedure: one does not need to install any external convolution package, and the cross section and its grid representation are identical for any scale choices, regardless of the statistics used. The latter feature is beneficial also in view of the fact that results for central and non-central scale choices in MadGraph5 aMC@NLO are fully correlated: hence, the ratios of predictions obtained with different scales are more stable than each of them individually. We point out that both of these advantages are relevant to the factorisation scale dependence. As far as the renormalisation scale is concerned, the situation is simpler, since no convolution integral is involved. Therefore, the optimal approach 6 would be that of using three grids, and construct G (h,R,l) j 1 j 2 j 3 (relevant to the ξ R dependence) on the fly, using G (h,B,l) j 1 j 2 j 3 . However, the gain of doing so w.r.t. our four-grid implementation is extremely marginal, and thus we did not consider it when constructing aMCfast.

Results and validation
As an illustration of the flexibility of aMCfast, in this section we present predictions for a variety of LHC processes, which either are currently or might soon become relevant for PDF determinations. Predictions are given for a 14 TeV LHC, using as inputs the NNPDF2.3 NLO PDF set [51] (associated with α S (M Z ) = 0.1180) in the case of five or four light flavours, and the NNPDF2.1 NLO PDF set [52] (associated with α S (M Z ) = 0.1190) in the case of three light flavours. We have made this choice of input parameters in order to be definite: the pattern of our findings would be unchanged had we used e.g. different PDF sets [2,10,12]. For the same reason, we do not consider PDF uncertainties in our illustrative study, and thus we only employ the central PDF member of the NNPDF2.3 and NNPDF2.1 sets. All such sets are taken from the LHAPDF5 [53] library.

General strategy
The main idea is the following: given a process, an observable O, and a binning for the differential distribution in O, we compute the value of its h th bin (for all bins) σ (h) O , in two different ways: directly, by means of MadGraph5 aMC@NLO, and a posteriori, using the grids constructed with the MadGraph5 aMC@NLO -aMCfast -APPLgrid chain. For the latter chain to be validated, these two results must be identical up to numerical inaccuracies; we shall call the former the reference result, and the latter the reconstructed result. In other words, we regard the l.h.s. of eq. (19) as computed directly with MadGraph5 aMC@NLO, and its r.h.s. as obtained from the interpolating grids, eqs. (54)- (58), so as to establish numerically the accuracy of the equality in eq. (19). As we have already stressed, we expect that inaccuracies are solely due the interpolating procedure, and not to a (possible) lack of statistics in the simulations; importantly, as was discussed in sect. 2.4, in our approach this property must hold for any non-central scale choices as well. In order to document this, all of our results will be presented for three different scale choices 7 : and as obtained from two different runs: one with "low" statistics, and one with "high" statistics. For each choice of scales and type of run, we shall plot the reference result, and the ratio of the reconstructed over the reference result. From that ratio, we shall see that the typical numerical inaccuracies due to the use of the interpolating grids are in the lower 10 −4 , with occasional larger differences in the case of the low-statistics runs, due to a still-insufficient amount of information stored in the grids 8 . The crucial point is the following: even with very limited statistics, the reconstruction of the reference results is extremely good, and this for distributions whose quality is largely insufficient for any phenomenological application. This proves that, for all practical purposes, the level of agreement between reconstructed and reference results is independent of the statistics employed.
In keeping with the strategy at the core of all applications of MadGraph5 aMC@NLO, the only operations to be performed by the user when running with aMCfast are related to the processspecific analysis. In particular, in view of the fast interface to fNLO results, the user's analysis (to be stored, as for any other kind of fNLO user utility, into the directory MYPROC/FixedOrderAnalysis -see sect. 3.2.2 of ref. [37] for more details) will need to contain, for each observable that he/she wants to consider, the instructions for filling the histogram that will contain the reference result (such instructions are identical to those of an fNLO analysis used when no interface with aMCfast is considered), and those for setting up the associated interpolating grids. Several templates of analyses that can be used with aMCfast will be provided in MYPROC/FixedOrderAnalysis (starting from the next version of MadGraph5 aMC@NLO), and the user should be able to easily adapt them to suit the needs of any given calculation. For the whole procedure to work, one will need to install, on top of MadGraph5 aMC@NLO and of the bridge code aMCfast, also APPLgrid version 1.4.56 or higher.
Parameter value Parameter value Table 1: Grid parameters used in the MadGraph5 aMC@NLO -aMCfast -APPLgrid validation runs presented here. In the case of tt production (sect. 3.2), the values of µ max and N τ have been set equal to 10 4 GeV and 50, respectively.
All of the calculations reported below are carried out either with the default physics model of fNLO computations in MadGraph5 aMC@NLO, which corresponds to the SM with N F = 4 light flavours, or with its N F = 5 or N F = 3 variants; the quarks which are not light (the top and, depending on the value of N F , the bottom and the charm) are given a non-zero mass. The parameters used to initialise the interpolating grids are given in table 1: from the first and second lines there, and by using eqs. (46) and (47), one immediately obtains the actual interpolation ranges in the variables y i and τ , and the corresponding grid spacings by taking the ratio of those over (N y − 1) and (N τ − 1), respectively. In the low-and high-statistics runs we employ a number of phase-space points per integration channel of the order of 10 3 and 10 6 respectively. For each process we shall consider a couple of observables to be definite, but we emphasise again that the number and the nature of the observables that one can compute are the user's choice. The kinematical cuts will in general be similar to those used by ATLAS or CMS in some analogous experimental analyses; however, a comparison with experimental data is beyond the scope of this work. In all cases, we shall set µ = H T /2, with H T the scalar sum of the transverse energies of all the final-state particles. We point out that the choice of scales may have some implications on the settings of the corresponding interpolation parameters.
For example, as one can see from table 1, in the case of tt production the values of µ max and N τ have been increased w.r.t. those used for all the other processes, and this in order to attain a comparable numerical precision. On the other hand, we have verified that, by choosing scales which do not depend on the event kinematics, and are of the order of the hardness of the given process at threshold, setting N τ = 10 is amply sufficient. To put things in perspective, however, we emphasise that even by using N τ = 30 in tt production with µ = H T /2, numerical inaccuracies are at most 5·10 −4 and generally smaller than that, i.e. better than any PDF-fit application will ever need. We start by studying top-quark pair production, which provides information on the gluon PDF in the large-x region complementary to those provided by jet-production data. We shall then consider isolated-photon production in association with jets, a process that allows a variety of tests of perturbative QCD and that is also directly sensitive to the gluon PDF. We shall next discuss the production of a lepton pair in association with one extra hard jet. We shall then turn to considering Zbb production, a case-study of how the complete automation of fast NLO QCD calculations allows one to include in global PDF fits arbitrarily complicated processes at no significant extra cost. Finally, we shall discuss in some details another important process at the LHC for PDF extraction, namely the production of a W boson in association with charm quarks, which is directly sensitive to the poorly known strange-quark PDF, and for which data from ATLAS and CMS have recently become available.

Top-quark pair production
Top-quark pair production provides one with unique information on the poorly known large-x gluon PDF [13,23,35]. Unlike the case of the Tevatron, where tt production is mostly driven by qq scattering, at the LHC the gg-initiated process is dominant, contributing to more that 85% of the total inclusive cross section. The recent calculation of the total cross section to NNLO [54], as well as the upcoming NNLO results for differential distributions, make top-quark pair production an essential ingredient for present and future global PDF analyses. A wide range of data from the LHC for this process are available, both for inclusive cross-sections [55][56][57] and for differential distributions [58,59], and much more will be measured in the near future.
As usual, the production of stable tt pairs in MadGraph5 aMC@NLO is simulated by starting with the generation of the process. Since in the default model used by MadGraph5 aMC@NLO there are four light flavours, one needs first to import a proper massless-bottom model. This is done by MG5 aMC> import model loop sm-no b mass MG5 aMC> define p = g u d s c b u~d~s~c~bH ere, the second command instructs MadGraph5 aMC@NLO that the proton contains five light flavours, and does include the bottom quark; by doing so, partonic subprocesses will be generated that feature bottom quarks in the initial state. After this, one may proceed with the generation command; MG5 aMC> generate p p > t t~[QCD] followed by the standard output and launch commands (see sect. 3 of ref. [37] for more details). There are n l = 7 contributing parton luminosities. For each of them, we report in table 2 the corresponding l n rs (r, s) Table 2: Values of l and (r, s) for which T (l) rs = 1, as generated by aMCfast in tt production. For a given l, the number n rs of non-null terms in the sum on the r.h.s. of eq. (51) is also given.
(r, s) contributions, which are the terms on the r.h.s. of eq. (51) with T (l) rs = 1; the number of such terms, denoted by n rs , is also reported in this table. As was already mentioned in sect. 2.3, the assignments of table 2 are automatically determined by aMCfast using the information provided by MadGraph5 aMC@NLO. We point out that each of the lines in table 2 corresponds to the set of four grids of eqs. (54)- (57). It is possible to use them separately, e.g. in order to determine the relative contribution of each parton luminosity to the different bins of the kinematical distributions studied, although this feature has not been used in this paper.
The results of our validation are presented in figs. 1 and 2, where we show the rapidity of the top quark (y(t)), and the tt invariant mass (M (tt)), respectively. Each figure consists of two panels, with that of the left (right) obtained with a low-statistics (high-statistics) run. Each panel has a main frame, where the solid histogram is the MadGraph5 aMC@NLO reference result obtained with the central scales of eq. (63), and the dashed histograms are the MadGraph5 aMC@NLO reference results relevant to the two non-central scale choices of eqs. (64) and (65) -incidentally, we note that the latter two predictions are obtained by using the reweighting technique presented in ref. [48] (i.e., MadGraph5 aMC@NLO has been run only once). The three lower insets show the accuracy obtained with interpolating grids, since they display as dashed histograms the ratio of the reconstructed over the reference result, for each of the three scale choices. These histograms thus represent the validation of the MadGraph5 aMC@NLO -aMCfast -APPLgrid chain, with values equal to one equivalent to a fast computation of the cross section with zero interpolation errors. As one can see from the plots, we obtain in all bins and for all scale choices an accuracy of 3 · 10 −4 at worst in the case of the high-statistics run. In keeping with the general discussion given at the beginning of sect. 3, the reconstruction accuracy of the low-statistics results is also quite good, in spite of the fact that the corresponding differential cross sections display extremely large fluctuations, and are unsuitable for phenomenology.

Photon production in association with one jet
The production of direct photons in hadronic collisions is sensitive to the gluon PDF, owing to the dominance of QCD Compton-scattering-like diagrams. Thanks to tight isolation requirements applied to both theory prediction and data, one is able to significantly reduce or even eliminate completely the cross section of non-direct photons, due to poorly-known fragmentation-fuction contributions in the former, and to the contaminating π 0 → γγ decays in the latter. A recent re-analysis of all available isolated-photon collider data (inclusive in additional QCD radiation) has shown that NLO QCD gives a good description of all the data sets, and in addition it can provide one with some constraints on the gluon PDF in the region of Bjorken-x relevant to the calculation of the gluon-fusion Higgs production cross section [22]. Therefore, there is a solid case to include photon data in the next generation of global PDF analyses. Data for isolated-photon production has been released by ATLAS and CMS [60][61][62]; in particular, the latest ATLAS measurement [60], that makes use of the complete 7 TeV run statistics, extends their kinematical coverage in the photon transverse energy up to a value of E T (γ) ∼ 600 GeV (see also ref. [63] for a study of the sensitivity of this measurement to various sets of PDFs).
On top of inclusive isolated-photon production, that of isolated photons in association with extra jets allows one to extend the sensitivity to the gluon PDF to a wider range of Bjorken-x, and also provides one with some information on the quark PDFs. Some measurements of isolated photon plus jets have been reported by ATLAS and CMS [64][65][66]; while these data are still not precise enough to directly constrain PDFs [67], the increase of the statistics collected in future LHC runs will dramatically change the picture, so that this process will become a useful addition to global PDF fits.
We have thus computed γ+jet production with MadGraph5 aMC@NLO. The computation is performed with five light flavours, so the generation command must be preceded by the same import and define commands which have been used in the case of tt production (see sect. 3.2). On top of those, we must also use here: MG5 aMC> define j = g u d s c b u~d~s~c~bs ince a light jet must have the same parton contents as the proton. The generation command is then: MG5 aMC> generate p p > a j [QCD] We have imposed the photon to be hard and central, p T (γ) ≥ 80 GeV and |η(γ)| ≤ 3. We use the Frixione photon-isolation prescription [68], since in this way one sets identically equal to zero the contribution of the fragmentation-function component in an infrared-safe manner. The parameters of the isolation (which in MadGraph5 aMC@NLO can be easily modified at the run-card level) are set as follows: ǫ γ = 1.0 and n = 1 (see eq. (3.4) of ref. [68]; note that, in an hadronic environment, the energies and angles of that equation have to be formally replaced by transverse energies and  table 3 are not equal, but differ by a trivial overall factor (the charge squared of a up-type vs an down-type quark). The combination of these two contributions could therefore be achieved by simply relaxing the condition that all non-zero T (l) rs elements be equal to one. While this can easily be done in the context of a specific computation, its general implementation in an automated code able to deal with any user-defined process is not straightforward.
The validation plots, presented in figs. 3 and 4, follow the same pattern as those for tt production. The representative observables chosen here are the photon transverse momentum (p T (γ), fig. 3) and the photon pseudorapidity (η(γ), fig. 4). As in sect. 3.2, the agreement between the reference and reconstructed results is at the level of 10 −4 or better for the high-statistics run, and only slightly worse for the low-statistics one.

Dilepton production in association with one jet
We now turn to studying the hadroproduction of a lepton pair in association with one jet; the leptons originate from an intermediate virtual Z or γ boson, and we shall denote the pair with the shorthand notation Z ⋆ in what follows. Inclusive Z-boson production has been used for a long time in PDF fits, since it provides one with a clean handle on the quark flavour separation in the proton, and various quality measurements from the Tevatron [71,72] and the LHC [15,73] are available. On the other hand, by requiring the presence of an additional extra jet one enters in a different domain, since at the Born level one is dominated by qg channels, which thus allow one to probe directly the gluon PDF. In this context, it is convenient to study the large-p T (Z ⋆ ) region: in fact, the fraction of gluon-initiated events increases 9 with p T (Z ⋆ ), and furthermore by doing so one is safely away from the small-p T (Z ⋆ ) region, where fixed-order calculations do not give physically-meaningful predictions. Another advantage associated with this kinematic region is that ratios of cross sections for W over Z production in association with jets can result in clean constraints on quarks and antiquarks PDFs at large-x [26], where they are poorly known (and especially so for the antiquarks).
There are various ATLAS and CMS measurements available for Z+jets production, which typically focus on the comparison between theory and data as a test of perturbative-QCD predictions which are capable of describing a significant amount of hard radiation [74][75][76][77]. However, from the PDF point of view, it is more interesting to perform the analysis as inclusively as possible (bar for the first jet, which one needs in order to be sensitive to the gluon PDF, as mentioned above), and to concentrate on the precise measurement of the Z ⋆ transverse momentum, possibly in bins of different rapidity y(Z ⋆ ). It is particularly convenient to perform such a measurement with a leptonically-decaying Z, which offers a cleaner scenario w.r.t. that of a hadronic decay; there are ongoing analyses in ATLAS and CMS -see e.g. [78].
We have computed dilepton+jet production at fNLO with MadGraph5 aMC@NLO in the fiveflavour scheme, following the generation command 10 : MG5 aMC> generate p p > mu+ mu-j [QCD] which has been preceded by the same import and define commands relevant to sect. 3.3. We have also imposed the following cuts: The list of parton luminosities relevant to this process is the same as that for γ+jet production, and is given in table 3. As in the γ+jet case, jets are reconstructed by using the anti-k t algorithm with R = 0.4, p T (j) ≥ 30 GeV, and |η(j)| ≤ 4.4. For the aMCfast validation we have chosen the transverse momentum of the lepton pair and its invariant mass, p T (µ + µ − ) and M (µ + µ − ), which we display in figs. 5 and 6 respectively. The latter clearly shows the small-M enhancement due to the presence of an intermediate photon. In the former, the feature around p T (µ + µ − ) ∼ 30 GeV is due to the presence of a sharp kinematic threshold only at the Born level (induced by the p T (j) cut), which causes this region to be infrared-sensitive [79]. Although this is an issue from the physics point of view for perturbative calculations, it does not pose any problems for the reconstruction through interpolating grids. We observe again an excellent agreement between the reference and the reconstructed results.

Z production in association with a bb pair
The production of a Z boson in association with a bottom-antibottom quark pair is an interesting example of how the automation of a fast interface to NLO QCD computations achieved by aMCfast can be a valuable tool to easily include arbitrarily complicated process into a global PDF fit. A related application could be that of studying to which degree the recent discrepancies between theory and experiment (see e.g. ref. [80]) for the production of Z bosons in association with bottom quarks might be absorbed by a change in PDFs. The hadroproduction of W bb and (Z/γ * )bb in the MadGraph5 aMC@NLO framework has been previously discussed in ref. [81], without however considering any aspects relevant to PDFs.
The generation of the current process in MadGraph5 aMC@NLO is achieved through the command: MG5 aMC> generate p p > z b b~[QCD] Note that the default model adopted by MadGraph5 aMC@NLO for NLO computations treats the bottom quarks as massive particles, and one works in a four-flavour scheme: hence, no models are explicitly imported, and no multi-particle labels redefined, at variance with was done previously. This is important in order to be able to perform the computation of Zbb production without resorting to any kind of jet-reconstruction algorithm, to probe the b quarks down to zero transverse momentum, and to define a fully-inclusive cross section. The list of partonic luminosities that contribute to this process is given in table 4. We validate the results of the present simulation by considering the rapidity distribution of the Zbb system (shown in fig. 7), and the transverse momentum of the bb pair  (shown in fig. 8). As in the previous cases, we observe no significant deviations of the reconstructed from the reference predictions.

W production in association with charm quarks
The final representative example of an LHC process of interest for PDF fits that we discuss in this work is the hadroproduction of W bosons in association with charm quarks. This process is directly sensitive to the PDF of the strange quark (see e.g. refs. [82,83]), which among the light-quark PDFs is the most poorly determined. In the majority of global PDF fits, the strange PDF is constrained by neutrino data [84][85][86] (in particular by the so-called dimuon process, i.e. charm production in charged-current DIS). The use of neutrino DIS data in PDF fits has several drawbacks, such as the need to model charm fragmentation, and to understand charm-quark mass effects at low scales.
It is therefore quite interesting that W c hadron-collider data provide one with a clean and robust independent check of the strange PDF determined elsewhere. In addition to that, the differences in the production rates of W + and W − bosons in association with charm quarks give information on the strange asymmetry in the proton, that plays an important role in the explanation [87] of the NuTeV anomaly [86]. Data on W production in association with either charm jets or D mesons have been made public by ATLAS [88] and CMS [24], based on the 7 TeV 2011 data sets, and analogous measurements from the 8 TeV run are ongoing. In particular, the CMS W c data has been recently used in a QCD analysis [16], together with HERA data, to show that the strange-quark PDF from collider-only data can be determined with a precision comparable to that of global fits which include neutrino data. Another recent analysis of the compatibility of the LHC W c data with existing fixedtarget DIS and Drell-Yan data has been presented in ref. [89].
Once again, we now compare a couple of reference distributions from MadGraph5 aMC@NLO to those obtained with the a-posteriori aMCfast -APPLgrid convolution. For simplicity, we restrict ourselves to considering W +c production; on the other hand, the process we shall actually deal with is that where an electron and a neutrino, rather than a W , are present in the final state, so as to make sure that all production spin correlations are correctly taken into account. Incidentally, this also shows that one is not restricted to performing computations in the narrow width approximation; in particular, the complex mass scheme [90] is fully supported in MadGraph5 aMC@NLO. Similarly to the case of the b quark in Zbb production discussed in sect. 3.5, in the present case it is best to treat the charm as a massive particle. Since in the default model used by MadGraph5 aMC@NLO the charm is massless, one needs to import a proper massive-charm model before proceeding with the generation of the process. This is done by means of the following commands: MG5 aMC> import model loop sm-c mass MG5 aMC> define p = g u d s u~d~sM G5 aMC> generate p p > e+ ve c~[QCD] In this list, the second command instructs MadGraph5 aMC@NLO that the proton contains three light flavours, and does not feature the charm quark; by doing so, no partonic subprocesses will be generated that feature a charm in the initial state. The following cuts have been imposed: p T (e + ) ≥ 10 GeV , η(e + ) ≤ 2.5 .
The list of partonic luminosities that contribute to this process in given in table 5.
The two representative distributions that we use for validation are the pseudorapidity of the positron, (η(e + ), presented in fig. 9), and its transverse momentum (p T (e + ) , presented in fig. 10). The same conclusions as in all previous cases apply here.

Conclusions and outlook
In this paper we have presented aMCfast, a tool that serves to achieve the complete automation of fast NLO QCD computations for arbitrary processes, by constructing a bridge between Mad-Graph5 aMC@NLO, which performs the core automated cross section calculations, and APPLgrid, which allows one to represent a given observable in term of interpolating grids that can be used for a fast and a-posteriori evaluation, and one in which the possibility is given of changing the scales and PDFs w.r.t. those used in the original computation. Since fast calculations of NLO cross sections are an absolute necessity in the context of global PDF fits, the combination of the automated programs discussed in this paper covers all the present and, more importantly, future needs relevant to PDF analyses.
We have validated our approach by considering a representative sample of production processes at the LHC (some of which not previously available in the form of fast computations), and by  fig. 1, for the e + pseudorapidity, in W +c hadroproduction.
comparing for several observables the bin-by-bin agreement between the original results of Mad-Graph5 aMC@NLO, and those reconstructed by means of aMCfast and APPLgrid. The typical accuracy that we find is of the order of some parts in 10 −4 , which is in fact far more than sufficient for any kind of phenomenological application. Moreover, our approach allows one to consider scale variations without any loss of precision, and without the need for interfacing APPLgrid with an external program to perfom PDF evolution.
A natural outcome of this work is the use of the grids produced using aMCfast in actual global PDF analyses, in view of the possible inclusion into the latter of LHC data, and especially of those which have not been considered so far, for lack of either sufficiently precise measurements, or of theoretical NLO calculations suited to the task. Given that the format of the results of the aMCfast -APPLgrid interface is a trivial extension of that produced by APPLgrid interfaced with other codes, aMCfast should be straighforward to employ by the global-PDF fitters already experienced with the former code.
The outlook for the near future features two important developments. Firstly, there is no conceptual difference between what has been done here, and its analogue in the case of a perturbative expansion in either the electroweak coupling constant α, or simultaneously in α S and α (mixed-coupling expansion). From the technical viewpoint, the only implication of a mixedcoupling expansion is the necessity of introducing additional interpolating grids, which is trivial. This point is important in view of the fact that MadGraph5 aMC@NLO is expected to be able to handle soon both kind of computations; thus, any extension of the MadGraph5 aMC@NLO capabilities will be made immediately available for fast computations, thanks to straightforward modifications to aMCfast. This will pave the way to a systematic inclusion of electroweak NLO corrections into global PDF fits; cases where such an inclusion may prove particularly relevant include high-p T jet production, the high-mass Drell-Yan cross section, as well as the consistent treatment of photon-initiated processes that are necessary for the determination of the photon PDF [91]. Secondly, the scale-and PDFindependent coefficients computed by MadGraph5 aMC@NLO, and used by aMCfast to set up the interpolating grids at fixed order, are quite similar to those relevant to the MC@NLO short-distance cross sections. Therefore, the current structure of aMCfast will need only a minor upgrade to be able to deal with NLO+PS computations. Apart from its obvious phenomenological spinoffs, and the possibility of obtaining PDF fits that include resummation effects well tuned to experimental data, such a project is also interesting from a more theoretical viewpoint, in that it would help spur a thorough investigation of the effects of the PDFs in initial-state showers; there is some evidence that these effects are small, but more systematic studies should be carried out, in view of the fact that they cannot be taken into account by any method based solely on short-distance cross section information. The public version of the aMCfast code will soon be available at: http://amcfast.hepforge.org/ Meanwhile, potential users are kindly requested to contact the authors.