text: Matching NLO QCD computations and parton shower simulations – arXiv Vanity
  Matching NLO QCD Computations and Parton Shower Simulations111Work supported in part by the UK Particle Physics and Astronomy Research Council and by the EU Fourth Framework Programme 'Training and Mobility of Researchers', Network 'Quantum Chromodynamics and the Deep Structure of Elementary Particles', contract FMRX-CT98-0194 (DG 12 - MIHT).
  Laboratoire d'Annecy-le-Vieux de Physique de Particules
  Chemin de Bellevue, BP 110, 74941 Annecy-le-Vieux CEDEX, France
  E-mail: On leave of absence from INFN, Sez. di Genova, Italy Stefano.F Bryan R. Webber
  We propose a method for matching the next-to-leading order (NLO) calculation of a given QCD process with a parton shower Monte Carlo (MC) simulation. The method has the following features: fully exclusive events are generated, with hadronization according to the MC model; total exclusive rates are accurate to NLO; NLO results for distributions are recovered upon expansion in αs; hard emissions are treated as in NLO computations while soft/collinear emissions are handled by the MC simulation, with the same logarithmic accuracy as the MC; and matching between the hard- and soft/collinear-emission regions is smooth. A fraction of events with negative weight is generated, but unweighting remains possible with reasonable efficiency. The method is clarified using a simple toy model, and illustrated by application to the hadroproduction of W+W− pairs.
  ††preprint: Cavendish–HEP–02/01
  LAPTH–905/02
  GEF–th–2/2002
  hep-ph/0204244
  The reliable prediction of cross sections and final-state distributions for QCD processes is important not only as a test of QCD but also for the design of collider experiments and new particle searches. All systematic approaches to this problem are based on perturbation theory, usually truncated at next-to-leading order (NLO), or in a few cases next-next-to-leading order (NNLO). Such calculations yield the best available results for sufficiently inclusive observables, in kinematic regions where higher-order corrections are not enhanced. However, in many cases a more exclusive description of final states and/or a wider kinematic coverage are needed.
  For the description of exclusive hadronic final states, perturbative calculations have to be combined with a model for the conversion of partonic final states into hadrons (hadronization). Existing hadronization models are in remarkably good agreement with a wide range of data, after tuning of model parameters [1, 2] . However, these models operate on partonic states with high multiplicity and low relative transverse momenta, which are obtained from a parton shower [3, 4] or dipole cascade [5] approximation to QCD dynamics and not from fixed-order calculations.
  To extend the coverage of perturbative calculations into regions where higher-order contributions are enhanced, one needs to identify and resum the enhanced terms to all orders. The resummed expression can then be combined with the fixed-order result, after subtracting resummed terms of the same order to avoid double-counting. This has been done explicitly for many important observables and processes, but a general procedure has not yet been developed. In principle, the parton shower provides such a procedure, if the problem of double-counting can be overcome. Combining NLO and parton-shower computations also offers the possibility of interfacing to hadronization models, thus overcoming the main deficiencies of existing calculations.
  Methods for combining fixed-order and parton-shower or dipole-cascade computations in QCD have been discussed in refs. [6] - [21] . Many of those papers concern the rather different problem of combining leading-order (tree level) matrix elements for hard multi-jet production with parton showers or dipole cascades for jet fragmentation. The only fully-implemented schemes for matching NLO computations and parton showers that have been used to produce phenomenological results [12] - [15] have been based on the so-called phase-space slicing method for regularising the infrared divergences encountered at NLO. The slicing method is simple but has the disadvantage of introducing a slicing parameter, which should in principle be sent to zero at the end of the calculation but in practice cannot be chosen too small. This can introduce uncontrolled errors and/or large negative event weights. The slicing method has therefore tended to be replaced where possible in NLO calculations by subtraction methods, in which no approximations need to be made and event weights can be controlled by improving the choice of subtraction terms.
  The subtraction method has been considered in refs. [11, 18, 19, 20] . The main concern of those papers was the correct matching of NLO matrix element calculations to the scale and scheme for next-to-leading logarithmic (NLL) parton showers. Up to now, no satisfactory algorithms for NLL showering have been developed. When such algorithms are available, the methods proposed there will make it possible for event generators to generate the NLO evolution of parton distribution functions (PDFs) and hard process cross sections, without the explicit input of NLO PDFs. Existing parton shower generators, however, correspond only to LL evolution, and therefore the matching criteria of refs. [11, 18, 19, 20] operate beyond the level of precision that can be achieved at present.
  In the present paper, we present a method for matching NLO matrix elements and existing parton shower generators based on the subtraction method for NLO calculations. First, in sect. 2, we outline the general objectives and limitations of our approach. The method is then illustrated in sect. 3 in the context of a simplified toy model. The case of QCD is dealt with in sects. 4 and 5. In sect. 4 results are presented for the process of W+W− production in hadronic collisions. In sect. 5 we discuss briefly the application of the method to other, more complex, processes. Concluding comments are given in sect. 6. All the technical details are collected in the appendices, where we also comment on the matching between NLO computations and parton shower approaches in the framework of the slicing method.
  It is undeniable that parton shower Monte Carlos (MC's) provide more general and flexible tools compared to NLO computations. By running an MC we get events, which we believe to be a faithful description of the real events detected by real experiments; this is one of the reasons why the MC's are so popular with experimentalists. On the other hand, NLO computations emphasize the role of (infrared-safe) observables; the presence of negative contributions, and the necessity of considering inclusive quantities in order to get rid of infinities, make it very problematic to even talk about events. Still, NLO computations have their virtues; they can handle hard emissions, and they can estimate total rates with better accuracy than MC's.
  In this paper, we aim to construct an MC that works like a regular MC, but on top of that knows how to treat hard emissions, and can compute rates to NLO accuracy. For brevity, we call this tool an MC@NLO. Clearly, there is a lot of freedom in constructing an MC@NLO. We shall follow a minimal approach. Namely, we want to keep the original MC as unmodified as possible; all major changes have to be carried out on the NLO codes. Amongst other things, this implies that our MC@NLO will not improve the logarithmic accuracy of the original MC, in those regions of the phase space where resummation is needed. As far as the NLO calculations are concerned, we shall not perform any approximations, neither at the level of matrix elements, nor at the level of phase space; thus, we shall adopt the subtraction method in order to cancel the infrared divergences that arise in the intermediate steps of the computation.
  To be more definite, we require MC@NLO to have the following characteristics:
  The output is a set of events, which are fully exclusive.
  Total rates are accurate to NLO.
  NLO results for all observables are recovered upon expansion of MC@NLO results in αs.
  Hard emissions are treated as in NLO computations.
  Soft/collinear emissions are treated as in MC.
  The matching between hard- and soft-emission regions is smooth.
  MC hadronization models are adopted.
  In sect. 4, we shall construct an MC@NLO fulfilling these conditions. However, we find that the technical complications of QCD tend to hide the basic ideas upon which our approach is based. Thus, we prefer to start from an oversimplified (and unrealistic) case, where the structure of the MC@NLO is apparent; this is done in sect. 3. The reason for doing this is that we shall see that the approach that leads to the MC@NLO in a simple case can actually be adopted in the QCD case as well, which will only be technically more involved.
  2.2 Negative weights
  We do not require the MC@NLO to be positive definite; events may have negative weights. However, it should be pointed out that these negative-weight events are of a completely different nature from the negative contributions that appear in NLO computations. In particular, the distributions of the positive- and negative-weight events are separately finite in our MC@NLO (the same is not true in an NLO computation). This implies the possibility of unweighted event generation,222In this paper, we adopt a slightly different notation from that commonly used by experimentalists. Namely, we call unweighted events any set of events whose weights are all identical, up to a sign. In other words, the weights of our unweighted events are of the form ±w, where w is a constant. as is customary in ordinary MC's. The situation is actually completely analogous to that occurring in polarized collisions, where physical observables are obtained as the difference of two positive-definite and finite quantities. Furthermore, we shall see that in our case the number of negative-weight events is reasonably small, and thus the number of MC@NLO events necessary to get a smooth distribution is comparable to that in ordinary MC's.
  2.3 Double counting
  When the MC generates events with real parton emission, it generates kinematical configurations that are also taken into account by the NLO computation: the possibility of having the same kinematical configuration from the MC and from the NLO may lead to double counting. We shall adopt the following quantitative definition of double counting:
  An MC@NLO is affected by double counting if its prediction for any observable, at the first order beyond the Born approximation333By the Born approximation we mean the lowest order at which the hard process contributes to any observable. in the expansion in the coupling constant, is not equal to the NLO prediction.
  According to this definition, double counting may correspond to either an excess or a deficit in the prediction, at any point in phase space. This includes contributions from real emission and virtual corrections. Generally speaking, the MC fills the phase space efficiently only in the soft and collinear regions. In these regions, the leading behaviour of the kinematical distribution of one-particle emission is given by the Born cross section, times a kernel that describes the soft or collinear emission. This holds both for the MC and for the NLO. Thus, in the soft and collinear regions of the one-particle-emission phase space, the MC and NLO results coincide when only the leading terms are considered. Elsewhere, the MC is not reliable and the NLO provides the best estimate. The problem is to merge these two descriptions.
  3 Toy model studies
  In this section we study a toy model that allows a simple discussion of the key features of the NLO computation, of the MC approach, and of the matching between the two. We assume that a system can radiate massless particles (which we call photons), whose energy we denote by x, with 0≤x≤xs≤1, xs being the energy of the system before the radiation. After the radiation, the energy of the system is x′s=xs−x. The system can undergo several further emissions; on the other hand, one photon cannot split further.
  In a perturbative computation, the Born term corresponds to no emissions. The first non-trivial order in perturbation theory gets contribution from those diagrams with one and only one emission, with either virtual or real photons. We write the corresponding contributions to the cross section as follows:
  (dσdx)B = Bδ(x), (1)
  (dσdx)V = a(B2ϵ+V)δ(x), (2)
  (dσdx)R = aR(x)x, (3)
  for the Born, virtual, and real contributions respectively; a is the coupling constant, possibly times a colour factor; B and V are constant with respect to x, and
  limx→0R(x)=B. (4)
  The constant B appears in eqs. (2) and (4) since we expect the residue of the leading singularity of the virtual and real contributions to be given by the Born term, times a suitable kernel. We take this kernel equal to 1, since this simplifies the computations, and it is not restrictive. Finally, ϵ is the parameter entering dimensional regularization in 4−2ϵ dimensions.
  The task of predicting an infrared-safe observable O to NLO accuracy amounts to computing the following quantity
  ⟨O⟩=limϵ→0∫10dxx−2ϵO(x)[(dσdx)B+(dσdx)V+(dσdx)R], (5)
  where O(x) is the observable as a function of x, possibly times a set of Θ functions defining a histogram bin. By "infrared-safe" here we simply mean that the integral in eq. (5) exists. The main technical problem in eq. (5) is due to the presence of the regularising parameter ϵ; in order to have an efficient numerical procedure, it is mandatory to extract the pole in ϵ from the real contribution, in this way cancelling analytically the pole explicitly present in the virtual contribution.
  3.1 NLO: slicing versus subtraction
  Within the context of this toy model, the usual methods of cancelling the pole in ϵ are as follows. In the slicing method, a small parameter δ is introduced into the real contribution in the following way:
  ⟨O⟩R=∫δ0dxx−2ϵO(x)(dσdx)R+∫1δdxx−2ϵO(x)(dσdx)R. (6)
  In the first term on the r.h.s. of this equation we expand O(x) and R(x) in Taylor series around 0, and keep only the first term; the smaller δ, the better the approximation. On the other hand, the second term in eq. (6) does not contain any singularity, and we can just set ϵ=0. We obtain
  ⟨O⟩R = aBO(0)∫δ0dxx−2ϵx+∫1δdxO(x)(dσdx)R+O(δ) (7)
  = a(−12ϵ+logδ)BO(0)+a∫1δdxO(x)R(x)x+O(δ,ϵ). (8)
  Using this result in eq. (5), we get the NLO prediction for O as given in the slicing method:
  ⟨O⟩slice=BO(0)+a[(Blogδ+V)O(0)+∫1δdxO(x)R(x)x]+O(δ). (9)
  The terms O(δ) cannot be computed, and are neglected; in general, one must check that these neglected terms are numerically small, by plotting ⟨O⟩slice versus δ in a suitable range in δ.
  In the subtraction method, no approximation is performed. One writes the real contribution as follows:
  ⟨O⟩R=aBO(0)∫10dxx−2ϵx+a∫10dxO(x)R(x)−BO(0)x1+2ϵ. (10)
  The second term on the r.h.s. does not contain singularities, and we can set ϵ=0:
  ⟨O⟩R=−aB2ϵO(0)+a∫10dxO(x)R(x)−BO(0)x. (11)
  Therefore, the NLO prediction as given in the subtraction method is:
  ⟨O⟩sub=BO(0)+a[VO(0)+∫10dxO(x)R(x)−BO(0)x]. (12)
  We rewrite this in a slightly different form, which will be more suited to discussion of the matching with the MC:
  ⟨O⟩sub=∫10dx[O(x)aR(x)x+O(0)(B+aV−aBx)]. (13)
  3.2 Toy Monte Carlo
  In a treatment based on Monte Carlo methods, the system can undergo an arbitrary number of emissions (branchings), with probability controlled by the Sudakov form factor, defined for our toy model as follows:
  Δ(x1,x2)=exp[−a∫x2x1dzQ(z)z], (14)
  where Q(z) is a monotonic function with the following general properties:
  0≤Q(z)≤1,limz→0Q(z)=1,limz→1Q(z)=0. (15)
  Note that the MC is also well defined without the last condition in eq. (15). We discuss later the specific functional form assumed for Q. If xs is the energy of the system before the first branching occurs, then Δ(x,xs) is the probability that no photon be emitted with energy z such that x≤z≤xs. At fixed x, this probability tends to 1 when a→0; that is, no emission is possible in the zero-coupling limit. On the other hand, when a→∞, we have Δ(x,xs)→0, which means that it is impossible that the system does not emit. The toy MC code we have written implements ordered emissions in the photon energies444This may violate energy-momentum conservation, as is usual in a leading logarithm implementation. and proceeds through the following steps:
  Define x0 as the lower bound on the photon energy after a branching; thus x0 plays the same role as the cutoffs used to define the Sudakov form factor in QCD MC's.
  Define xM to be the maximum energy available to the photon; before the first branching, we have xM=xs, where xs is the energy of the system.
  Pick a random number R1∈[0,1], and solve for z the equation
  exp[−a∫xMzdz′z′]≡(zxM)a=R1⟹z=xMR1/a1. (16)
  Pick a random number R2∈[0,1]. If Q(z)<R2, let xM=z, reject the branching, and return to step 2).
  If z<x0, reject the branching and exit. If z>x0, a photon with energy z has been emitted. We let xM=z, and iterate the procedure going back to 2).
  This procedure ensures that the photon emissions have the distribution prescribed by the Sudakov form factor in eq. (14). As can be seen from that equation, the probability of having soft emission (i.e., with z<x0), and thus of terminating the shower, is Δ(x0,xs), which tends to one when a→0 (weak coupling), and when xs→x0 (no energy left for non-soft branching). In this toy model, the hardest emission is always the first one.
  3.3 Matching NLO and MC
  In order to introduce some useful notation, we start from the LO perturbative result, which does not pose any problems, since it is in fact what is implemented in ordinary MC's. At the LO, the system has not yet lost any energy to perturbative photon emission. Thus, the maximum energy available to the photon in the first branching is 1. For each event, we set xM=1, and proceed as described in sect. 3.2. The total exclusive rate will have to be equal to B, according to eq. (1); thus, if we generate N events, each event has a weight equal to B/N.
  We can formally read this procedure from eq. (13), by setting a=0 and performing the formal substitution
  O(x)⟶IMC(O,xM(x)), (17)
  where IMC(O,xM) stands for interface-to-MC, and indicates symbolically the distribution in the observable O as obtained by running the MC starting from a given xM, and giving each event the weight 1/N; its precise definition is given in app. B.1. For x=0 we need xM=1; we do not need to specify xM(x) for other values of x at present. Interfacing LO with MC, we thus obtain
  (dσdO)MC@LO=BIMC(O,1). (18)
  This equation reminds us of the total rate (B), and the maximum energy available for photons at the first branching (xM=1).
  3.3.1 Naive subtraction
  The first attempt at extending eq. (18) to NLO amounts simply to substituting eq. (17) into eq. (13). We get
  (dσdO)naive=∫10dx[IMC(O,xM(x))aR(x)x+IMC(O,1)(B+aV−aBx)]. (19)
  This equation suggests the following procedure:
  Pick at random 0≤x≤1.
  Generate an MC event with xM(x) as maximum energy available to the photon in the first branching; attach to this event the weight wEV=aR(x)/x.
  Generate another MC event (a "counter-event") with xM=1; attach to this event the weight wCT=B+aV−aB/x.
  Repeat the first three steps N times, and normalize with 1/N.
  Unfortunately, this procedure is bound to fail, since the weights wEV and wCT diverge as x→0. A small cutoff ρ can then be introduced; as x→ρ, we would have wEV≃aB/ρ and wCT≃−aB/ρ. Both these weights are very large in absolute value; this means that the eventual unweighting procedure, performed by the MC, will be highly inefficient. Furthermore, as we shall see in sect. 3.4, eq. (19) has problems of double counting, defined according to the criterion in sect. 2.3.
  3.3.2 Modified subtraction
  To avoid the problems outlined above, we define the modified subtraction method by
  (dσdO)msub = ∫10dx[IMC(O,xM(x))a[R(x)−BQ(x)]x (20)
  +IMC(O,1)(B+aV+aB[Q(x)−1]x)].
  This amounts to saying that, in eq. (19), we subtract and add the quantity
  IMC(O,xM)aBQ(x)x, (21)
  using xM=xM(x) in the first and xM=1 in the second term introduced in this way. We point out that Q is the same function that appears in the definition of the Sudakov form factor, eq. (14). By virtue of the second condition in eq. (15), the weights wEV and wCT attached to events and counter-events are now separately convergent as x→0.
  The two terms involving Q(x), that we inserted in eq. (20), are not identical; therefore, the procedure which leads from eq. (19) to eq. (20) is not a subtraction in the usual sense of an NLO computation. Still, these two terms do not contribute to the observable O at O(a), because they are compensated by analogous terms due to the parton shower BIMC(O,1); this will be shown with explicit examples in sect. 3.4, and in general in app. B.3. We point out that the function aQ(x)/x that appears in eq. (21) does not depend upon the MC cutoff x0. It follows that the dependence upon x0 of the MC@NLO, eq. (20), is identical to that of the ordinary MC, being entirely due to IMC, and is therefore power-like.
  It should be clear from the description of the procedure given above that the dependence upon O in eqs. (19) and (20) is purely formal, serving only to allow us to write such equations in a compact form. We stress that, in the actual implementation of all the MC@NLO's described in this paper, the computation of observables proceeds exactly as in standard MC's: namely, events are generated, showered, and hadronized without reference to any specific observable.
  Equation (20) is our master equation for the toy model. We shall show in what follows that the MC@NLO defined in this way meets the requirements given in sect. 2.1. In sect. 4, we shall use eq. (20) to construct a QCD MC@NLO by analogy with the toy model.
  3.4.1 Exclusive observable
  As an example of an "exclusive" observable in the toy model, we shall consider the quantity
  y=max(x1,…,xn), (22)
  where n denotes the number of emissions in one given event (in the case of no emission, the event is not used to fill the histograms). This observable is the energy of the hardest photon of the event. It is particularly convenient because in the MC approach it coincides with the energy of the first photon emitted in the shower; when MC@NLO is considered, y is the larger of the energy of the photon emitted perturbatively, and that of the first photon in the shower (if any). This fact allows us to obtain analytical results for MC, NLO and MC@NLO.
  The NLO prediction for this observable is (here and in the following, we consider only y>x0 in order to simplify the analytical results)
  (dσdy)NLO=aR(y)y, (23)
  while the MC result, normalized to the Born cross section, is
  (dσdy)MC=aBQ(y)yΔ(y,1), (24)
  which is just the differential cross section at x=y, computed using eq. (14), times the probability that there is no harder emission at x>y.
  In order to study the MC@NLO results, we have to consider the quantities IMC(y,1) and IMC(y,xM(x)) that will be substituted into eq. (19) or (20). In the case of IMC(y,1), the hardest photon comes from the MC and we have
  IMC(y,1)=aQ(y)yΔ(y,1). (25)
  In the case of IMC(y,xM(x)), the hardest photon may come from the NLO contribution or the MC, and so we find
  IMC(y,xM(x)) = Δ(min{x,xM(x)},xM(x))δ(y−x) (26)
  +aQ(y)yΔ(y,xM(x))Θ(y−x)Θ(xM(x)−y).
  As a shorthand notation, in the following we shall denote the contributions to dσ/dy from IMC(y,1) and IMC(y,xM(x)) by dσS/dy and dσH/dy respectively, since they correspond to a Standard MC evolution (i.e., with no prior NLO real emission), and to a Hard MC evolution (whose initial condition includes one NLO real emission). This notation will be extensively used from sect. 4.1 on. In the case of the naive subtraction method, eq. (19), these contributions are separately divergent, but their sum is not. We find
  (dσSdy)naive=aQ(y)yΔ(y,1)[B+aV−aB∫10dxx]. (27)
  The results for dσH/dy depend somewhat on the form of the function xM(x) introduced in eq. (17). We define xe to be the solution of x=xM(x), i.e., xe=xM(xe), and we assume this solution to be unique. We further assume xM(x) to be a monotonically decreasing function; thus, x<xM(x) for x<xe. Then:
  ⋄ for y<xe:
  (dσHdy)naive=aR(y)yΔ(y,xM(y))+a2Q(y)y∫y0dxxR(x)Δ(y,xM(x)); (28)
  ⋄ for y>xe:
  (dσHdy)naive=aR(y)y+a2Q(y)y∫x−1M(y)0dxxR(x)Δ(y,xM(x)). (29)
  Expanding in a, we find
  (dσdy)naive≡(dσSdy)naive+(dσHdy)naive=ay[R(y)+BQ(y)]+O(a2). (30)
  Therefore the naive subtraction method suffers from double counting, according to the definition in sect. 2.3, especially at small y, where Q(y)→1.
  For the modified subtraction method, eq. (20), we find instead that both contributions are finite:
  (dσSdy)msub=aQ(y)yΔ(y,1)[B+aV+aB∫10dxQ(x)−1x], (31)
  (dσHdy)msub = ay[R(y)−BQ(y)]Δ(y,xM(y)) (32)
  + a2Q(y)y∫y0dxx[R(x)−BQ(x)]Δ(y,xM(x));
  (dσHdy)msub=ay[R(y)−BQ(y)]+a2Q(y)y∫x−1M(y)0dxx[R(x)−BQ(x)]Δ(y,xM(x)). (33)
  Expanding in a, we now get
  (dσdy)msub=aR(y)y+O(a2), (34)
  in agreement with the NLO prediction for all values of y. Furthermore, for y→0 we have the expected resummation of leading logarithms,
  (dσdy)msub⟶(aBy+O(a2))Δ(y,1)+PSTfor y→0, (35)
  where PST stands for power-suppressed terms, i.e., terms that are not singular for y→0. We stress that the behaviour of dσ/dy in eq. (35) is dictated by eq. (31), i.e., the contribution of H events is power suppressed for y→0, since (R(y)−BQ(y))/y=O(1) in this limit. We also notice that the functional form in eq. (35) is the same as that in eq. (24). Thus, the MC@NLO resums large logarithms of y in the same way as does the MC. A general derivation of this result for a wide class of observables is given in app. B.4.
  3.4.2 Inclusive observable
  As an "inclusive" observable, we shall consider the fully inclusive distribution of the photon energies; thus, each photon emitted, either at the NLO level or by the MC shower, will contribute to the physical observable, which we shall denote by z. The NLO and MC results (for z>x0) are:
  (dσdz)NLO = aR(z)z, (36)
  (dσdz)MC = aBQ(z)z. (37)
  Eq. (36) coincides with eq. (23), which is a trivial consequence of the fact that only one emission occurs at NLO. On the other hand, the MC result in eq. (37) differs from the exclusive case, eq. (24), since all emissions now contribute, not only the first. As a consequence, eq. (37) has a logarithmically-divergent integral at z→0, whereas this does not happen in eq. (24).
  In order to obtain the predictions of the MC@NLO, we need to find the analogues of eqs. (25) and (26), in the case of the variable z. We obtain
  IMC(z,1) = aQ(z)z, (38)
  IMC(z,xM(x)) = δ(z−x)+aQ(z)zΘ(xM(x)−z). (39)
  The term with the δ-function in eq. (39) accounts for real emission at the NLO level, which always contributes to the inclusive spectrum when it occurs; the other term is the contribution of the MC shower.
  To compute the z distribution for the naive subtraction method, we use eq. (38) and eq. (39) in eq. (19), to obtain
  where σtot is the NLO total cross section
  σtot=B+aV+a∫10dxx[R(x)−B]. (41)
  As before, we may expand in a to find
  (dσdz)naive=az[R(z)+BQ(z)]+O(a2), (42)
  which shows again that the naive subtraction method suffers from double counting.
  For the modified subtraction method we use eq. (38) and eq. (39) in eq. (20) to find
  (dσdz)msub = az[R(z)−BQ(z) (43)
  az+Q(z)(σtot−a∫1x−1M(z)dxR(x)−BQ(x)x)].
  Expanding in a we now find
  (dσdz)msub=aR(z)z+O(a2); (44)
  thus, as in the case of the observable y, the MC@NLO prediction at O(a) coincides with the NLO result for all z. The small-z limit of eq. (43) is
  limz→0z(dσdz)msub=aσtot. (45)
  From eq. (37) we obtain instead
  limz→0z(dσdz)MC=aB. (46)
  Therefore, at z→0 the MC@NLO result is equal to that of the standard MC, times σtot/B, which is the K-factor for the total rate.
  In this section we present toy model results for NLO, MC, and MC@NLO. The following parameters have been used in eqs. (1), (2) and (3):
  a=0.3,B=2,V=1. (47)
  For the real emission distribution in eq. (3) we have taken
  R(x)=B+x(1+x/2+20x2), (48)
  so that the total cross section is
  σtot=B+a(V+9512). (49)
  Note that this implies a large K-factor, σtot/B=2.34. Finally, for the MC cutoff parameter we used
  x0=0.02, (50)
  and for the maximum energy function in eq. (17) we chose the form xM(x)=1−x. In the left panel of fig. 1 we show the NLO result obtained with eq. (13) for the observable y defined in eq. (22) (solid histogram), compared to the result obtained by running the MC (dashed histogram), with Q(x)≡1; the integral of the latter result has been normalized to the NLO rate (i.e., each event is given a weight equal to σtot/N). As a cross check for the numerical implementation against analytical results, we also plot as solid lines the functions given in eqs. (23) and (24) (for the latter, B has been substituted with σtot, in order to be consistent with the normalization of the dashed histogram), which in fact appear to be in perfect agreement with the corresponding histograms. The shapes of the NLO and MC results are completely different. In fact, the function R(x) has been chosen to introduce a clear difference between the NLO and MC results, for a better understanding of the effect of the matching between NLO and MC.
  Figure 1: Left panel: Toy model NLO (solid histogram) and MC results (dashed histogram), overlaid on analytical results of eqs. (23) and (24), with Q(x)≡1. Right panel: MC results, in the case of the dead zone with xdead=0.6. Q is chosen as in eq. (52) (solid), eq. (53) (dashed), and eq. (54) (dotted); the first bin is normalized to unity.
  To study the matching of NLO and MC, we consider an MC with a dead zone. In our toy model, this corresponds to an MC that cannot emit photons with energy x>xdead (in the numerical simulations we set xdead=0.6). We can easily get a toy MC with a dead zone by suitably choosing the function Q in eq. (14). We first introduce a smoothing function of the form
  G(x)=c2(1−x)2βx2α+c2(1−x)2β, (51)
  where α, β and c are free parameters, and then consider three different implementations of the dead zone:
  1.aaQ(x)=Θ(xdead−x); (52)
  2.aaQ(x)=Θ(xdead−x)G(x/xdead),withα=1,β=1,c=1; (53)
  3.aaQ(x)=Θ(xdead−x)G(x/xdead),withα=2,β=1,c=8. (54)
  The results for the pure MC evolution are presented in the right panel of fig. 1; the dependence upon Q is actually fairly big, since we made extreme choices in eqs. (52)–(54). It is thus striking that the corresponding results obtained with the MC@NLO, shown in the left panel of fig. 2, do not display big differences from each other. All three MC@NLO curves are quite close to the pure NLO result (and coincide with it for y>xdead), and differ sizably from it only in the region of small x, where the effect of multiple emission is expected to be important.555Since y is an exclusive variable, and the total rate is that of the NLO result, the differences between the MC@NLO and NLO curves at y>x0 are in fact compensated by the region y<x0 in the total rate. In this region, however, the three MC@NLO curves are close to each other, notwithstanding the differences shown in the right panel of fig. 1. We are here in a region where the leading logarithmic behaviour of the MC is highly dominant, and thus the impact of Q (which is subleading in logarithmic accuracy) is moderate. We see that the curve corresponding to eq. (52) displays a less regular behaviour than the others at y≃xdead; however, the effect is very moderate, being at worst of O(a2).
  Figure 2: Toy model MC@NLO results; a dead zone is present (xdead=0.6) and Q is as in eq. (52) (solid), eq. (53) (dashed), or eq. (54) (dotted). The results are for the exclusive observable y (left panel), and for the fully inclusive observable z (right panel). The normalization is absolute, and the dot-dashed histogram is the pure NLO result.
  The results for the inclusive z distribution are presented in the right panel of fig. 2. Again, we can see that the three MC@NLO curves are remarkably close to each other, in spite of the differences in the implementations of the dead zone. The ratio of the cross section in the leftmost bin of the solid histogram over the NLO result for the same quantity is 2.23. This value is fairly close to the K-factor σtot/B=2.34, that we should find according to eqs. (45) and (46); the small difference is due to fact that this bin is not quite at the limit z→0.
  We can thus conclude that MC@NLO reduces the numerical dependence upon subleading logarithms, which is seen when running the MC alone; this is because those regions where the subleading terms are becoming numerically important get sizable contributions from the hard NLO emissions in the MC@NLO, which are not present in the MC, and which eventually become dominant. On the other hand, in the soft-dominated regions, the shapes of the MC@NLO curves just reproduce those of the MC curves, as we expect.
  4 The QCD case: W+W− production
  We now face the problem of constructing an MC@NLO in QCD. In order to simplify as much as possible all the technical details related to the colour structure of the hard process, we consider here the production of W+W− pairs in hadronic collisions. This process also has a fairly simple structure of soft and collinear singularities. However, we argue that the method we propose is actually valid for any kind of hard production process.
  We shall use HERWIG [1] as our reference parton shower MC. Other choices are of course possible, and our approach is not tailored to work with HERWIG, although some of the technical details, reported in the appendices, do depend on the specific implementation of the hard process and parton showering in that program.
  In order to define our MC@NLO for W pair production, we rely on the results previously presented for the toy model. First we emphasize the main points of the toy MC@NLO. Then we make a parallel between the contributions to the toy model cross section, and the contributions to the W+W− cross section. Using these results, we proceed to construct the QCD MC@NLO, by analogy with the toy MC@NLO.
  4.1 Toy model revisited
  Let us consider again our master equation for the toy model MC@NLO, eq. (20). The MC@NLO approach consists in running the ordinary MC, with initial conditions and weights determined by the NLO calculation. Depending upon the initial conditions, we distinguish two classes of events, which correspond to IMC(O,xM(x)) and IMC(O,1). We denote the class corresponding to IMC(O,xM(x)) by H, which stands for "hard", and reminds us that one emission already occurred, due to NLO, prior to the MC run. The class corresponding to IMC(O,1) is denoted by S, which stands for "standard"; this reminds us that this initial condition is the same as for the ordinary MC runs.
  It is important to realize that for S events the finite part of the virtual contribution, aV, the real counterterm, −aB/x, and the term aBQ(x)/x, serve exclusively to fix the normalization of this class of events, relative to that of the class H. In other words, it is irrelevant that the latter two terms explicitly depend upon x: in the MC@NLO, the spectrum of the observable O resulting from S events is entirely due to MC evolution. In contrast, the NLO hard emissions do contribute to the kinematics of H events, and thus help to determine O.
  In order to obtain in QCD a structure analogous to that outlined above, we therefore have to collect all the terms of the NLO cross section that do not contribute to hard emissions, and use them to compute the normalization of class S events; these events result from MC evolutions whose initial conditions are identical to those of an ordinary MC. On the other hand, the kinematics of NLO hard emissions determine the initial conditions for MC evolutions for H events.
  Finally, we have already seen in the case of the observables y and z (see also app. B.3 for a more formal proof) that at O(a) the dependence upon Q(x) disappears. Roughly speaking, this happens because the contribution due to the Born term B, in class S events, cancels at O(a) with the contribution of the terms ±aBQ(x)/x. We observe that, in the toy model, aBQ(x)/x is the O(a) term in the expansion of the MC result (see eqs. (24) and (37)). We therefore guess that the role of aBQ(x)/x in a QCD MC@NLO is played by the O(αs) term in the expansion of the result of the ordinary MC, when only the shower is taken into account (i.e., the hadronization step is skipped). We shall eventually prove (see app. B.3) that this is in fact the case.
  4.2 NLO cross section
  We now deal explicitly with W+W− production in hadronic collisions, and we start by writing the NLO cross section. In this paper, we use the results of ref. [22] , which we rewrite in a form more suited to our present goal; it should be clear that any W+W− cross section (or at least, one based upon the subtraction method) can be cast in the same form. Owing to the factorization theorem, we have
  dσ(P1,P2)=∑ab∫dx1dx2f(H1)a(x1)f(H2)b(x2)d^σab(x1P1,x2P2), (55)
  for the collision of hadrons H1 and H2 with momenta P1 and P2 respectively. The sum runs over those parton flavours that give partonic subprocesses contributing to W pair production at NLO. We symbolically denote these processes as follows:
  q¯q ⟶ W+W−, (56)
  q¯q ⟶ W+W−g, (57)
  qg ⟶ W+W−q, (58)
  where q=u, d,… and eq. (58) includes the corresponding antiquark processes. Equation (56) contributes to Born and virtual emission terms, whereas eqs. (57) and (58) contribute to real emission terms.
  The information on the hard process in eq. (55) is contained in the subtracted partonic cross sections, d^σab (the hat means that counterterms are implicit in the expression). We write these cross sections as follows (eq. (3.45) of ref. [22] ):
  d^σab=dσ(b)ab+dσ(sv)ab+d^σ(c+)ab+d^σ(c−)ab+d^σ(f)ab, (59)
  dσ(b)ab is the Born contribution. It has 2→2 kinematics, and corresponds to eq. (1) of the toy model.
  dσ(sv)ab is the contribution of the virtual diagrams, and of all those terms, resulting from various subtractions, that have 2→2 kinematics and are of O(αs). This term corresponds to the finite part of eq. (2) of the toy model.
  d^σ(c±)ab are the contributions of the finite remainders resulting from the subtraction of the initial-state collinear divergences (+ and − indicate radiation from partons a and b respectively). These terms do not have an analogue in the toy model, since in the toy model there are no parton densities. They have 2→3 kinematics, with the light outgoing parton always collinear to the beam-line; we refer to these configurations as quasi-2→2 kinematics. We stress that d^σ(c±)ab are subtracted quantities, and that the counterterms have 2→2 kinematics.
  d^σ(f)ab is the contribution of the real emission diagrams, with its divergences suitably subtracted; it has 2→3 kinematics (the counterterms have 2→2 or quasi-2→2 kinematics). The unsubtracted part of this term corresponds to eq. (3) of the toy model.
  We point out that all the terms on the r.h.s. of eq. (59) are separately finite. Also, they do not contain any large logarithms, and so no cancellation of large numbers occurs among them. This property is a feature of the subtraction method, and does not hold in the slicing method.
  Items a)-d) above achieve the identification between terms contributing to the QCD cross section and to the toy model cross section. We now go into more detail, by writing the quantities appearing in eq. (59) more explicitly .
  We denote by dϕ2(s) and dϕ3(s) the phase spaces of the W+W− pair and of the W+W− pair plus a light parton, respectively; they depend upon the kinematical variables of these particles, and upon s=(p1+p2)2, the c.m. energy squared of the partons that initiate the hard scattering, whose momenta we denote by p1=x1P1 and p2=x2P2. The Born term is
  dσ(b)ab(p1,p2)=M(b)ab(p1,p2)dϕ2(s), (60)
  where M(b)ab is the Born matrix element squared, summed over spin and colour degrees of freedom (the average factor for initial-state ones is understood), times the flux factor 1/(2s). The soft-virtual contribution is
  dσ(sv)ab(p1,p2)=[QabM(b)ab(p1,p2)+M(v,reg)ab(p1,p2)]dϕ2(s), (61)
  where Qab is a constant with respect to the kinematics of the partonic process, and M(v,reg)ab is the finite part of the virtual contribution. As with any other finite part resulting from the cancellation of singularities, its definition is to a certain extent arbitrary; however, any redefinition is compensated by suitable changes in the other terms contributing to the cross section. The initial-state collinear terms are:
  d^σ(c±)ab=dσ(c±)ab∣∣ev−dσ(c±)ab∣∣ct, (62)
  dσ(c+)ab(p1,p2)∣∣ev = αs2π[logs(1−x)2μ2P(0)ca(x)−P(1)ca(x)]M(b)cb(xp1,p2)dϕ2(xs)dx, (63)
  dσ(c−)ab(p1,p2)∣∣ev = αs2π[logs(1−x)2μ2P(0)cb(x)−P(1)cb(x)]M(b)ac(p1,xp2)dϕ2(xs)dx.
  In eq. (62), the term dσ(c±)ab|ct (the "ct" stands for counterterm) is constructed in order to cancel the divergences present in the collinear "event" contribution dσ(c±)ab|ev; these divergences are soft (x→1), as can be seen explicitly in eqs. (63) and (LABEL:cmev). The counterterm has a form that is fully specified in a given subtraction scheme (see appendix A). In eqs. (63) and (LABEL:cmev), μ is the factorization mass scale (which also enters in the parton densities), and P(i)ab is the ith term in the ϵ-expansion of the Altarelli-Parisi kernel in 4−2ϵ dimensions, taken for x<1 (i.e., no δ(1−x) terms are included). Finally, we write the real emission term as
  d^σ(f)ab=dσ(f)ab∣∣ev−dσ(f)ab∣∣ct, (65)
  dσ(f)ab(p1,p2)∣∣ev=M(r)ab(p1,p2)dϕ3(s), (66)
  where M(r)ab is the matrix element squared resulting from the real emission diagrams (summed over all spins and colours, and averaged over those relevant to incoming partons), times the flux factor. This quantity has soft and initial-state collinear singularities, and thus dσ(f)ab|ev has divergences, which are cancelled by the counterterms dσ(f)ab|ct in eq. (65). As in the case of dσ(c±)ab|ct, these counterterms are fully specified only in a given subtraction scheme. However, one should keep in mind that the physical cross section is always strictly independent of the choice of the subtraction scheme.
  4.3 Observables at NLO
  As in the case of the toy model, a preliminary step in the construction of the MC@NLO is the prediction of an observable at the NLO. We start by defining an unintegrated, fully-differential cross section
  d^Σab(x1,x2)=f(H1)a(x1)f(H2)b(x2)d^σab(x1P1,x2P2), (67)
  and analogous quantities, for all the contributions to the partonic cross section that appear on the r.h.s. of eq. (59). Any observable O can depend, at the NLO, upon the kinematical variables of the W+, of the W−, and of the light outgoing parton. We denote by O(3) the observable computed in non-soft and non-collinear kinematics; this corresponds to O(x) in the toy model. We also denote by O(2) and O(~2) the observable computed in 2→2 (i.e., soft) and in quasi-2→2 (i.e., collinear non-soft) kinematics respectively; these correspond to O(0) in the toy model. In all cases, we understand the variables of the outgoing particles to be given in the laboratory frame. We then have