Matching NLO QCD computations with PYTHIA using MC@NLO

We present the matching between a next-to-leading order computation in QCD and the PYTHIA parton shower Monte Carlo, according to the MC@NLO formalism. We study the case of initial-state radiation, and consider in particular single vector boson hadroproduction.


Introduction
Perturbation theory offers a systematic way to improve theoretical predictions for any given infrared-safe observable. Depending on whether the expansion parameter in the series is α S , or is α S times a logarithm or a logarithm squared of some function of the kinematics of the hard process, one obtains a fixed-order result or a resummed result, respectively. Cross sections expanded up to a certain order in α S and those resummed are relevant to complementary kinematic regions of the phase-space. It is therefore convenient to combine the features of these two expansions, by defining a matched cross section which is equal to the former or to the latter in the appropriate phase-space region.
Fixed-order results are now generally available at the next-to-leading order (NLO), which corresponds to including in the cross section the coefficients of terms of order α b S (the leading order or Born level) and of order α b+1 S . Thanks to the recent and rapid progress in the automated treatment of both the real and the virtual contributions to NLO computations, it is realistic to assume that phenomenological results at this accuracy will become available in the next few years for all of the reactions of interest to the LHC physics programme. The situation is much less encouraging for fixed-order cross sections at O(α b+2 S ) (or NNLO), where only a handful of results are available, for very small final-state multiplicities.
As far as resummed results are concerned, it is in general understood how to achieve a next-to-leading logarithmic (NLL) accuracy, which is equivalent to including terms proportional to α n S log kn Q and to α n S log kn−1 Q, with k = 1, 2 depending on the nature of Q; some results are also known to next-to-next-to-leading logarithmic accuracy (NNLL). Unfortunately, resummed computations are in general technically complicated, observable-dependent, and error-prone; although for some observable classes a semi-automated algorithm (CAESAR [1]) is available, the overall situation is far less satisfactory than for NLO computations. For this reason, a very appealing alternative is that provided by Parton Shower Monte Carlos (PSMCs). Although formally PSMCs are equivalent to a LL-accurate resummation (which may become NLL in some corners of the phase space, for some PSMCs and subject to certain restrictions), in practice they are known to do much better than that, as comparisons with data and with analytically-resummed results clearly show. Furthermore, PSMCs are fully flexible, give one the possibility of including hadronization models in a consistent manner, and are the workhorses of experimental collaborations, thanks to their capabilities of simulating fully-realistic final-state configurations that can undergo detector simulations.
As is well known, the approximations that form the core of PSMCs severely limit their predictive power in those phase-space regions (corresponding to multi-jet configurations) which are of interest for most of new-physics searches at colliders. These limitations can be alleviated by viewing PSMC predictions as resummed results, to be included with the NLO corrections to the relevant production processes into a matched cross section. The definition of a formalism for matching NLO computations and MC simulations has attracted a considerable amount of attention. There are now several proposals, but only two of them, namely MC@NLO [2] or POWHEG [3], have made it to the stage of actually implementing several hadroproduction processes, and of being routinely used by experimental collaborations.
The MC@NLO formalism requires the computation of the cross section predicted by the PSMC at O(α b+1 S ). Because of the structure of PSMCs, the non-trivial information of this computation is actually process-independent, and is contained in the definitions of the parton branchings; the process-dependent part is entirely factorized in the Born matrix element. Thus, one essentially has to perform one set of computations (since typically initial-and final-state branchings are treated differently by the Monte Carlos) per PSMC, in order to be able to match an NLO computation with a parton shower simulation according to the MC@NLO approach. So far, these computations have been carried out for the case of HERWIG (see refs. [2,4,5]) and, more recently, for HERWIG++ [6]. In this paper, we present the first results relevant to the matching with PYTHIA 6.4. We limit ourselves to considering the case of initial-state branchings, and present some sample results for the Drell-Yan process.
This paper is organized as follows. In sect. 2 we describe the various steps necessary for an MC@NLO matching with PYTHIA, from the definition of the underlying parton-level NLO cross section to the MC@NLO short-distance cross sections used for the generation of hard-subprocess events, to be showered by PYTHIA. We present sample results in sect. 3, and we give our conclusions in sect. 4.

NLO parton-level cross section
The starting point for the construction of MC@NLO is that of writing the shortdistance parton-level cross section according to the subtraction formalism of refs. [7,8] (which we shall call FKS subtraction henceforth). The default procedure in FKS is that of treating simultaneously (i.e. in one contribution to the partonic cross section) the two initial-state collinear singularities, due to one given final-state parton being collinear to the initial-state parton coming from the left or from the right. On the other hand, in FKS one can also treat these two singularities independently, by defining two separate contributions to the short-distance cross section, each of which corresponds to one of the initial-state singularities; this procedure is explained in detail in ref. [9]. As we shall discuss in the following, when an NLO computation is matched to PYTHIA according to the MC@NLO formalism, it turns out to be convenient to adopt the latter strategy (at variance with the case of HERWIG, where the simultaneous subtraction suffices).
We shall assign the momenta entering the (real-emission) partonic subprocesses as follows: with V representing a W or a Z boson, and a, b, and c being QCD partons; Born-like processes are simply obtained from eq. (2.1) by removing c(k 2 ) from the r.h.s.. Only soft and initial-state collinear singularities are present in the processes of eq. (2.1). Hence, the S functions required for a separate treatment of the collinear singularities are two and shall be denoted by with the following properties: The expectation value for any observable O will then be written as where (see eqs. (4.14)-(4.16) of ref. [2]) As can be seen in ref. [2], the Σ terms in equation above are defined as the partonic short-distance cross sections, times the luminosity factors. In writing eq. (2.8), we have split (using the S functions) the Born (Σ (b) ab ) and soft-virtual (Σ (sv) ab ) contributions into two terms, and have associated them with the corresponding real-emission contributions (see refs. [5,9]). We have also used the collinear limits of eqs. (2.4) and (2.5), that imply that the S functions associated with the collinear counterterms are trivial. Finally, note that the collinear remainders Σ (c±) ab need not be multiplied by the S functions. Equation (2.8) must be further manipulated in order to use it in MC@NLO; in particular, one has to apply the so-called event projection which allows one to define a unique kinematic configuration associated with all counterevents (for a given realemission configuration) -see sect. A.4 of ref. [2]. As discussed in that paper, although event projection is largely arbitrary, in the context of MC@NLO it turns out to be convenient to derive it from the behaviour of the Monte Carlo one interfaces to. We shall therefore discuss in the next subsection the behaviour of PYTHIA relevant to this issue.

Event projection with PYTHIA
The way in which PYTHIA deals with the simulation of V production is the following [10,11]. First, the hard process is generated. This implies the generation of the two Bjorken x's entering such a process, which we shall denote by ζ 1 and ζ 2 for the parton coming from the left and from the right respectively. The partonic c.m. energy squared is with S the collider energy squared. In the case of single-V production, s 0 = m 2 V , with m V the mass of V (or its virtuality if lepton pair production is considered), but eq. (2.9) is obviously valid regardless of the production process. PYTHIA then begins the showers. Choosing on statistical basis which leg emits "first" (such an emission is the only one that matters as far as MC@NLO is concerned), the shower variables are related to the momenta given in eq. (2.1) as follows: (2.10) (2.11) with i = 1, 2 for the emissions from the parton coming from the left and from the right respectively. The branching generated with eqs. (2.10) and (2.11) is such that the Bjorken x's associated with the initial-state partons of eq. (2.1) (which we shall denote by z 1 and z 2 ) are for an emission from leg 1, and for an emission from leg 2. In the branching, the c.m. energy of the Born-level subprocess is conserved (which, in the present case, is equivalent to the requirement that the virtuality of V be a constant), and thus An easy way to achieve this is that of imposing that the Bjorken x's of the Born-level subprocesses, ζ i , be separately conserved. As shown in ref. [2], the event-projection procedure can be formally determined by constructing two "observables" that are conserved in the branching process. According to what has been discussed above, a possible choice is: It is worth noting that while eq. (2.15) is the same as in HERWIG, eq. (2.16) does not coincide with either of the HERWIG choices considered in ref. [2]; the implication of this fact is that we obtain two different event projections, depending on whether it is leg 1 or 2 that emits. We stress that eqs. (2.15) and (2.16), and the resulting manipulations we are now going to describe, do not depend on the nature of the Born-level final state (a V boson rather than -say -a three-jet configuration), owing to the fact that PYTHIA adopts the so-called s-approach [10] when doing initial-state branchings for all hadroproduction processes. In the context of the sapproach, the invariant mass ( √ s 0 ) of the final-state system at the Born level is kept constant during the branching. This allows one to formally replace V with the set of final-state particles at the Born level, and k 1 with the sum of their four momenta, which does not entail any changes to eqs. (2.9)-(2.14), from which we ultimately derive eqs. (2.15) and (2.16). The procedure described here is thus fully general, and is not restricted to V -boson production. We recall that event projection is a way to re-write the soft and collinear counterterms in a pure-NLO computation, so as all counterterms associated with a given S contribution at the real-emission level (in the present case, S + dΣ (f ) or S − dΣ (f ) ) have the same kinematics. In order to achieve this, it turns out to be convenient to define the counterterm kinematics starting from a fixed real-emission kinematics. This situation is depicted in fig. 2 for soft counterterms (the case of collinear counterterms is essentially identical, the only differences being in the assignments of the momentum fractions of the incoming partons, which we shall specify in what follows). Since the event projections we define here are ultimately motivated by the branching strategy of PYTHIA, the procedure of fig. 2 is by construction the inverse of that depicted in fig. 1. We shall return to this point in sect. 2.4.
In order to actually obtain event projections, we use the master equations (A.36) and (A.37) of ref. [2]. For an emission from leg 1, these equations read in our case 2 )) , (2.18) and their analogues for the collinear case (where one formally replaces s with c+ in the r.h.s. of eqs. (2.17) and (2.18)). Note that, having treated separately the two collinear singularities in the NLO cross section, the case c− need not be considered when studying branching from leg 1). We find: x (s) For an emission from leg 2, we find instead We shall do this in the next subsection.

MC subtraction terms
At the NLO in α S , the cross section resulting from PYTHIA is where the two terms on the r.h.s. account for emissions from leg 1 and 2 respectively. They read: Following what was done in sect. A.5 of ref. [2], we can manipulate the equations above to render them suitable for an integration together with the NLO parton-level cross section. First of all, one notes that the following identifications are valid owing to the construction of event projections: in eq. (2.24), and in eq. (2.25). We can therefore perform a change of integration variables: Here, we have introduced the rescaled virtualitieŝ Furthermore, we have explicitly included the factors Θ ± , which are related to the choice of the maximum virtuality allowed during shower evolution. We have e.g.
with y being the relevant angular variable associated with the FKS parton in the FKS subtraction (here the cosine of the angle between k 2 and p i for emissions from leg i), and 1 − x the (normalized) energy of the FKS parton. In the case of PYTHIA, at variance with HERWIG, x ≡ z.

MC@NLO short-distance cross sections
At this point, we have all the ingredients needed to generate hard-subprocess events according to the MC@NLO formalism, which will be subsequently showered by PYTHIA. Recalling that we denote by Σ (α) the contributions Σ (α) to the NLO partonlevel cross sections after event projection, event generation will be performed following the procedure outlined in sect. 4.5 of ref. [2] (see also sect. 3 of ref. [5]), using the following integrals: We remind the reader that all integrals defined in these equations are separately finite, and that σ tot = I is an exact equation (with σ tot the fully-inclusive NLO rate). Although eqs. (2.35) and (2.36) solve the problem of the generation of hardsubprocess events, a further simplification is possible. We start by considering the integrals I (±) H , that are used to generate H events. These two integrals have to be computed separately. However, one observes that, for a given choice of (z 1 , z 2 , φ 2 ), the kinematics associated with these two contributions are actually identical, and the integrals can therefore be computed together. This is equivalent to using the following integral for the generation of H events: The first term in the integrand has been simplified thanks to eq. (2.3 ab are evaluated in two different kinematics configurations, that eventually result after the first branching in the same real-emission kinematics (z 1 , z 2 , φ 2 ). This situation is precisely the one depicted in fig. 2.
A simplification analogous to that of eq. (2.38) is also possible in the case of S events, although the argument is slightly more involved. If one fixes (z 1 , z 2 , φ 2 ) (i.e., the real-emission kinematics), the kinematic configurations associated with I As we have already stressed, although there is an ample freedom in choosing P

(±)
H→S , it is best to adopt a form motivated by what the PSMC does when branching. In practice, it is therefore convenient first to obtain P (±) S→H from the PSMC, and then (by inverting them) P (±) H→S . In the case of PYTHIA, P (±) S→H amount to performing a boost of the Born-level four momenta in the transverse direction, to balance the transverse momentum of the parton produced in the branching (c(k 2 ) in eq. (2.1)), followed by a boost in the longitudinal direction, to e.g. the lab frame. The information on the latter is equivalent to the assignments of the momentum fractions of the incoming particles. The prescription for the transverse boost given above is obviously trivial in the case of V production (since k 1T = −k 2T ), but can be applied to final states with arbitrary multiplicity. Thus, what done here is also valid for final states more complicated than single V .
The two procedures related to eq. (2.39) and (2.40) are equivalent, since the integrals are performed over the whole phase space (possibly subject to kinematic restrictions, which are however identical in the two cases). Equation (2.40) corresponds to the situation depicted in fig. 1, except for the seemingly different assignments of the momentum fractions of the incoming particles. These differences are however immaterial, since such fractions are simply integration variables which can be manipulated and renamed at will. We shall now proceed to perform such manipulations, starting from a change of variables in eq. (2.36) Furthermore, consistently with eq. (2.40), the integrands have to be computed with the suitable kinematic configurations, 2 + or 2 − . We can take these two operations into account by the formal replacements It is easy to realize that, in the case of the Born, soft-virtual, and soft counterterms contributions, we have by construction since for these terms the jacobian appearing in eq. (2.42) exactly compensates the one in the definition of the Σ (α) ab terms, see eq. (4.18) of [2]. In general, this is not true for contributions with a purely collinear structure (such as the collinear counterterms to the real-emission matrix elements). Equations (2.41), (2.42), and (2.43) allow us to rewrite  for the generation of S events. One observes that in the integrand of eq. (2.46) there will be a term where we have made use of eq. (2.3). The same kind of simplification will occur in the sum of the two soft counterterms. As far as the collinear counterterms are concerned, we stress that the S ± terms that multiply them are equal to one (by construction of the S functions -see eqs. (2.4) and (2.5)). We arrive therefore at the following form: Equations (2.38) and (2.48) are our final expressions, used for the generation of H and S events respectively. The net result of the various manipulations carried out in this section is that it is still possible to treat simultaneously the two initial-state collinear singularities, as was the case when matching with HERWIG. We point out that the fact that eqs. (2.38) and (2.48) do not depend upon the S functions is not a property of V -boson production, but applies as well to all processes whose Born-level final-state particles are all colour singlets 1 . In fact, for this independence of S functions to occur, we only need eq. (2.3) to hold, which is true in the cases just mentioned (see ref. [9] for further details on the construction of S functions). Clearly, in order to make use of eq.  fig. 2). On the other hand, the short-distance cross sections for S events we started from, eq. (2.36), are associated with two different kinematic configurations, denoted by 1 ± in eq. (2.39) (i.e. the lower part of fig. 2). Thanks to the procedure described above, we have manipulated eq. (2.36) precisely to be able to associate the two ± contributions with the same kinematic configuration, denoted by 1 in eq. (2.40) (i.e. the upper part of fig. 1). Owing to the properties of event projections that we have discussed above, the whole procedure can obviously be carried out for any final-state multiplicity, with 2 and 1 formally replaced by the real-emission and Born-level configurations respectively. Finally, we stress that the manipulations performed here will also be valid in the case of strongly-interacting particles in the final state. The only difference is that, in such a case, eq. (2.3) will not hold any longer, and therefore the quantity S + + S − will appear in eqs. (2.38) and (2.48) as a factor multiplying the real-emission, soft counterterms, soft-virtual, and Born contributions (on the other hand, the purely collinear terms will be unchanged, owing to the fact that eqs. (2.4) and (2.5) hold regardless of the nature of the final state). Clearly, the contributions due to branchings of Born-level final-state strongly-interacting particles will be given by short-distance cross sections analogous to those of eqs. (2.38) and (2.48) (see e.g. sect. 3 of ref. [5]). The computation of such contributions is beyond the scope of the present work, and we postpone it to a future publication.
We conclude this section by stressing that, due to the form of the shower variables used by PYTHIA, the soft limit of the sum of the MC subtraction terms coincides with that of the real-emission matrix elements, which was not the case for HERWIG. This implies that here we can set G ≡ 1, where G is the function introduced in sect. A.5 of ref.
[2] -we refer the reader to that paper for a discussion on this issue.

Results
In this section, we present sample results relevant to W + production in pp collisions at √ S = 14 TeV. Our aim is not that of performing a phenomenological study, but rather that of presenting a few control plots that show that the matching of the NLO results for single vector boson production with PYTHIA according to the MC@NLO formalism works as we expect. This is non trivial, given the differences between HERWIG and PYTHIA, which in turn result in different short-distance MC@NLO cross sections, as discussed previously. We set m W = 80.4 GeV, Γ W = 2.14 GeV, and use CTEQ6.6 [12] PDFs. Our default scale choices are µ F = µ R = m T , where m T is the transverse mass of the W . When reconstructing jets, we adopt the k T -jetfinding algorithm of ref. [13], with Y cut = (10 GeV) 2 .
In each of the plots we present in figs. 3-5, we display three histograms. The solid (black) histograms are the results of MC@NLO with PYTHIA (which we shall call MC@NLO/PY) i.e. what has been computed in this paper. The dotted (red) histograms are the results of MC@NLO with HERWIG (which we shall call MC@NLO/HW). Finally, the dashed (blue) histograms are the results obtained with PYTHIA stan-dalone. When using PYTHIA for showering hard events, we set MSTP(81)=0 and MSTP(91)=0, which corresponds to switching off multiple interactions and primordial k T respectively. We switch off matrix element corrections by setting MSTP(68)=0, which also forces the maximum virtuality in the shower to be equal to the vector boson mass, when standalone generation is performed. On the other hand, hard subprocess events generated by MC@NLO/PY are given to PYTHIA in the standard Les Houches format [14]; we set PARP(67)=1 and SCALUP=m W to have the same maximum scale in the shower as in the standalone generation. In order to facilitate the visual comparisons between MC@NLO/PY and PYTHIA, the results of the latter are rescaled (by different factors, depending on the observables considered). As has already been discussed at length in the literature, we expect MC@NLO to be identical in shape to the PSMC results in the regions of the phase space dominated by those large logarithms the PSMC is able to resum; this is the motivation for comparing MC@NLO/PY with PYTHIA standalone. Also, we expect MC@NLO to coincide, in shape and normalization, with the NLO results in the regions where hard emissions are dominant. Given that it has already been shown in the past that this is the case for MC@NLO/HW, here we compare MC@NLO/PY directly with MC@NLO/HW, rather than with the pure NLO results, since this will also give us the opportunity to observe the different behaviours of the two underlying PSMCs in the soft/collinear regions. As far as negative weights are concerned, for W production at the LHC their fraction with MC@NLO/PY is about 0.6%, while in the case of MC@NLO/HW is about 8%. Since the phase-space parametrizations used in the two codes are identical, this large difference is essentially due to the different choices of shower variables made by the two PSMCs.
In fig. 3 we consider the transverse momentum distributions of the W + boson (left pane), and of the hardest jet of the event (right pane). These distributions are interesting since in the large-p T region they are dominated by hard emissions, whereas in the small-p T one there are logarithmically-enhanced terms that need be resummed (the pure NLO results diverge there). The comparisons among the three results for the two distributions follow the same pattern, as we expect. At small p T 's, the shapes of MC@NLO/PY and of PYTHIA are identical. To show this clearly, the latter results have been rescaled so as their first bins coincide with those resulting from MC@NLO/PY. On the other hand, the MC@NLO/PY and MC@NLO/HW results are quite different in this p T region, owing to the different treatment of soft and collinear emissions by the two PSMCs. It is known that PYTHIA tends to give softer spectra than HERWIG, which is therefore what we expect to and do find with MC@NLO/PY and MC@NLO/HW. At large p T 's, MC@NLO/PY coincides with MC@NLO/HW -both in shape and in absolute normalization, which is again what we expect.
In fig. 3 we also show the scale dependence of MC@NLO/PY, determined according to the following procedure. The renormalization-scale variations are de-  fig. 3 as the upper and lower bounds of the shaded areas, for the two observables considered there. In spite of the lack of statistics in the high-p T tails 2 the trend is clear: the fractional scale dependence grows from about ±5% at low p T 's to about ±10% at large p T 's. This has to be compared with the pure-NLO result (not shown here) that features a decrease in the scale dependence, from ±15% at low p T 's to ±10% at large p T 's; this behaviour results from a dependence on µ F almost constant w.r.t. p T , and a dependence on µ R decreasing with p T (due to the running of α S ). The scale dependences of MC@NLO/PY and of the pure-NLO predictions are therefore consistent at large p T 's. This is what we expect, since there the MC@NLO/PY result is basically coincident with the NLO one, the effect of the shower being negligible on these inclusive variables. The size of the scale dependence is also compatible with the fact that, in fixed-order perturbation theory, O(α S ) is actually the first order that contributes to p T > 0. This observation therefore applies also to low p T 's in the case of the pure-NLO predictions, but it does not in the case of MC@NLO/PY. In fact, MC@NLO/PY fills the low-and intermediate-p T regions mostly through the showering of S events. The shower, however, determines only the kinematics of the final-state configuration (i.e., the p T of the W and of the hardest jet here), but the weights are given by the short-distance cross section of eq. (2.48), and therefore receive both O(α 0 S ) and O(α S ) contributions. This implies that in the low-p T region one expects MC@NLO/PY to have a scale dependence smaller than that of the treelevel O(α S ) term alone, which is what we observe. We conclude this discussion by noting that other definitions can be given of the uncertainties associated with mass scales (e.g., variations may be summed linearly rather than in quadrature), without this changing the pattern found here. We also point out that the scale dependences of the rapidity observables which we shall discuss below are rather featureless (i.e., scale variations are independent of rapidities), and therefore they will not be shown in what follows. In fig. 4 we consider the rapidity distributions of the W + boson (left pane), and of the hardest jet of the event (right pane). In the absence of any cuts in p T , the boson rapidity is an inclusive variable unaffected by large logarithms, and both NLO computations and MC simulations should predict it relatively well. We do indeed see an overall consistency among MC@NLO/PY, MC@NLO/HW, and PYTHIA standalone. Differences among the three predictions are larger in the case of the rapidity of the hardest jet, since this is a less inclusive variable w.r.t. the W rapidity, and is also more sensitive to hadronization corrections. The PYTHIA standalone results have been rescaled by the NLO K-factor, σ NLO /σ LO , before any cuts are applied and jets are reconstructed. Finally, in fig. 5 we consider the difference between the rapidities of the W + boson and of the hardest jet, for two different cuts on the transverse momentum of the hardest jet. The PYTHIA standalone results have been rescaled in the same way as in fig. 4. This observable and its analogues (with y W replaced by y S , S being the system emerging from the hard process at the Born level -e.g. the tt pair in top-pair production) has attracted some attention in the past (see e.g. ref. [15] for a recent discussion), owing to the fact that MC@NLO/HW has a different behaviour around y S − y j ≃ 0 w.r.t. the underlying NLO computation -the former being flatter than the latter, or having a dip, depending on the nature of the system S. The "flatness" or the presence of a dip is a feature of MC simulations, and more specifically is a consequence of the choices made for initial conditions, as was done here for PYTHIA in eq. (2.34). It has to be stressed that this is true for both HERWIG and PYTHIA, as the dashed histogram on the right pane of fig. 5 shows, and therefore has nothing to do with the presence of dead zones 3 as such in HERWIG. In PYTHIA, the choice of initial conditions is much more flexible than in HERWIG, and one can make the dip of fig. 5 disappear by choosing a large-enough scale as the maximum virtuality allowed for the shower. In doing so, however, one may extend the collinear approximation, that is the core of both HERWIG and PYTHIA, outside its proper range of validity, and we consider the choice made in eq. (2.34) a sensible one. H events in MC@NLO, or matrix element corrections in HERWIG and PYTHIA, will give a substantial contribution to the region y S − y j ≃ 0 (see the solid and dotted histograms in fig. 5), where the predictions will become closer (w.r.t. those of MC simulations without matrix element corrections) to the pure-NLO results. Differences will in general remain, that can be formally ascribed to effects beyond NLO.
In order to further this argument, we consider an extension of the choice made in eq. (2.34). Namely, we parametrize the maximum virtuality allowed in the shower in terms of a number f introduced as follows: The plots shown so far have therefore been obtained with f = 1. In fig. 6 we consider again the difference in rapidity of fig. 5, with f = 1/2 (dotted red histograms) and with f = 2 (dashed blue histograms), together with our default choice f = 1 (solid black histograms); we point out that for the values of f considered here the fraction of negative weights is basically a constant. Two additional p T cuts are also considered on top of those of fig. 5. The results obtained with PYTHIA feature a very large sensitivity to the choice of f , while those obtained with MC@NLO/PY are fairly stable. One may be tempted to use the results of fig. 6 in order to "tune" the parameter f , and obtain a PYTHIA prediction in decent agreement with that of MC@NLO/PY. However, this is an a posteriori procedure, which can be justified phenomenologically, but which implies that any possibility is given up of a sensible estimate of the theoretical uncertainties affecting the observable considered here; also, such a procedure is in general observable-dependent. Indeed, the left pane of fig. 6 shows that, if no information on 2 → 2 matrix elements is included, this rapidity difference basically cannot be predicted by Monte Carlos at moderate and large p T 's. This is not surprising, given that one is attempting to use the collinear approximation outside its range of validity: the correct result can be recovered (which is equivalent to choosing the argument of a logarithm so as its numerical value coincides with a given constant), but only in a heuristic way. The situation improves if the correct information on matrix elements is used, as in MC@NLO/PY. This is analogous to what happens when one varies the renormalization and factorization scales, at the LO and the NLO levels. As in the case of scale variations, extreme values for f will lead to problems; however, f = O(1) appears to be a safe choice, allowing one to realistically estimate theoretical uncertainties. We conclude by mentioning that other observables (such as those shown in fig. 3 and 4) display the same pattern of dependence on f as the observable in fig. 6. We shall discuss this issue further [16] by considering also the case of Higgs production, where the effects are more pronounced.

Conclusions
We have presented the construction of the matching between an NLO QCD computation and the virtuality-ordered PYTHIA Monte Carlo, according to the MC@NLO formalism. We have limited ourselves to considering only the case of initial-state radiation, and applied the formalism to the study of W hadroproduction. Owing to the different structures of HERWIG and PYTHIA, the short-distance MC@NLO cross sections used to generate hard events are different in the two cases. However, the FKS subtraction, which is the method used in MC@NLO for dealing with infrared singularities, needs no modifications, and is able to treat both MCs. Likewise, no changes are needed in the MC@NLO formalism, and the matching with PYTHIA is achieved by performing a process-independent calculation. For the process considered in this paper, the fraction of negative weights in MC@NLO/PY is much smaller than that in MC@NLO/HW. From the physical viewpoint, the pattern of the comparison between MC@NLO, pure-NLO, and PSMC results is the same for MC@NLO/PY as for MC@NLO/HW -MC@NLO shows the same behaviour as the NLO or MC where either one is most reliable, with a smooth transition between the hard and soft-collinear emission regions.
Although the process studied in this paper is particularly simple, the formulae given here will be basically sufficient for performing the matching in the case of more complicated reactions (and limited to initial-state branchings only). Therefore, the results presented here are the first step towards the construction of MC@NLO/PY for generic processes, where also final-state radiation is present, and towards the extension of this formalism to the p T -ordered version of PYTHIA.