Four-lepton production at hadron colliders: aMC@NLO predictions with theoretical uncertainties

We use aMC@NLO to study the production of four charged leptons at the LHC, performing parton showers with both HERWIG and Pythia6. Our underlying matrix element calculation features the full next-to-leading order $O(\alpha_S)$ result and the $O(\alpha_S^2)$ contribution of the $gg$ channel, and it includes all off-shell, spin-correlation, virtual-photon-exchange, and interference effects. We present several key distributions together with the corresponding theoretical uncertainties. These are obtained through a process-independent technique that allows aMC@NLO to compute scale and PDF uncertainties in a fully automated way and at no extra CPU-time cost


Introduction
Light charged leptons constitute a particularly clean trigger for high-energy experiments. It is thus relatively easy to measure, despite their being small, the cross sections of processes that feature several final-state leptons. Prominent among such processes are those mediated by a pair of electroweak vector bosons which, depending on their identities, may give rise to a missing-energy signature as well. Vector boson pair production is interesting in at least two respects. Firstly, it is an irreducible background to Higgs signals, in particular through the W + W − and ZZ channels which are relevant to searches for a standard model Higgs of mass larger than about 140 GeV. While always smaller than the W + W − channel, ZZ decays may provide a cleaner signal due to the possibility of fully reconstructing the decay products of the two Z's. Secondly, di-boson cross sections are quite sensitive to violations of the gauge structure of the Standard Model, and hence are good probes of scenarios where new physics is heavy and not directly accessible at the LHC, yet the couplings in the vector boson sector are affected.
The Higgs being such an elusive particle, searches require a combination of data-driven and theoretical methods in order to overcome irreducible backgrounds. From the theory viewpoint, this essentially implies computations as accurate as possible, and ideally able to realistically reproduce actual experimental events. On the other hand, modifications to tri-linear gauge couplings will tend to manifest themselves at large transverse momenta, and therefore this is the region where theoretical results have to be reliable. The inclusion of next-to-leading order (NLO) QCD corrections into four-lepton cross section predictions is a minimal, and fairly satisfactory, way to fulfill both of the two conditions above.
The aim of this paper is that of studying four-lepton hadroproduction at the NLO accuracy in QCD, by also adding the (finite) contribution due to gg fusion (which is formally of next-to-next-to-leading order, NNLO, yet enhanced by the gluon PDFs), and by including the matching with parton showers according to the MC@NLO formalism [1]; this is done by means of the framework aMC@NLO [2,3] . Since aMC@NLO provides the full automation of the matching procedure, and of the underlying matrix element computations, we can give results for any four-lepton final state. In order to be definite, we have chosen to consider the process: with ℓ , ℓ ′ = e , µ. This choice is simply motivated by the fact that ZZ production as currently implemented in the MC@NLO package [4] does not feature production spin correlation and off-shell effects (which on the other hand are included, although in an approximated way, in the case of W − W + and W ± Z production in the codes of ref. [4]), and virtual-photon contributions with their interference with the Z's. The present aMC@NLO application remedies to these deficiencies by computing in an exact manner all relevant matrix elements (including the contributions of singly-resonant diagrams, which are potentially relevant to analyses in kinematical regions where one of the Z bosons is forced to be off-shell), and by including, as was already mentioned, the gg-initiated contribution as well which, although perturbatively suppressed, can be numerically important at the LHC owing to its dominant parton luminosity. We remind the reader that ZZ production at parton level and NLO accuracy in QCD has been studied for two decades now. The on-shell calculations of refs. [5,6] have been subsequently improved to include leptonic decays [7], singly-resonant diagrams [8], and anomalous couplings [9]; these results have been used recently to match the NLO computation to parton showers according to the POWHEG formalism [10]. The O(α 2 S ) process gg → ZZ was first computed in refs. [11,12,13]. These papers have been later superseded by the inclusion of off-shell effects and Z/γ * interference [14,15]. As far as MC@NLO results are concerned, and on top of the phenomenological results presented here, this paper is a first for three reasons: • The matching with Pythia6 [16] for a kinematically non-trivial process.
• The use of MadLoop [17] for the computation of loop-induced processes (i.e., the contributions of finite one-loop amplitudes squared).
• The computations of scale and PDF uncertainties with a reweighting technique.
In fact, the use of Pythia for the shower phase in MC@NLO has been limited so far to a proof-of-concept case [18], while the squares of loop amplitudes have been considered in MadLoop only for pointwise tests [17], and not in phenomenological applications. Finally, the reweighting procedure which we describe in this paper, and that is used to compute the scale and PDF dependences of the cross section, allows one to determine the corresponding uncertainties without requiring any additional computing time. We point out that the capabilities listed in the three items above are not process dependent and are fully automated.
This paper is organized as follows: in sect. 2 we discuss some of the technicalities relevant to four-lepton production, and the procedure implemented in aMC@NLO for the determination of scale and PDF uncertainties. In sect. 3 we present selected phenomenological results, and in sect. 4 we draw our conclusions. Further details on scale and PDF dependence computations are reported in appendix A.

Technical aspects of the computation
The framework of the aMC@NLO programme used for the phenomenology studies of this paper is unchanged w.r.t. that employed in ref. [3]; in particular, the underlying tree-level computations are performed with MadGraph v4 [19]. We remind the reader that aMC@NLO automates all aspects of an NLO computation and of its matching with partons showers. One-loop amplitudes are evaluated with MadLoop [17], whose core is the OPP integrand reduction method [20] as implemented in CutTools [21]. Real contributions and the corresponding phase-space subtractions, achieved by means of the FKS formalism [22], as well as their combination with the one-loop and Born results and their subsequent integration, are performed by MadFKS [23], which finally takes care of the MC@NLO matching [1] as well.
The novel features of aMC@NLO whose results we present here are the matching with the virtuality-ordered Pythia6 shower, based on ref. [18]; the possibility of determining through reweighting the scale and PDF uncertainties affecting our predictions, which we shall discuss in general in sect. 2.1 (with some further technical details given in appendix A); and the use of MadLoop for the computation of the O(α 2 S ) partonic process The matrix elements relevant to this process are UV-and IR-finite, and therefore the corresponding generated events are treated as an LO sample from the viewpoint of matching with parton showers. The amplitudes can be straightforwardly computed by MadLoop. In the current version, MadLoop assumes that these amplitudes will have to be multiplied by Born ones; we have extented the scope of the code, in order for it to compute amplitudes squared.
The fact that the process in eq. (2.1) is not a genuine virtual correction, i.e. it lacks an underlying Born, implies that event unweighting is essentially a brute-force operation (as opposed to the normal situation where the presence of an underlying Born allows one to pre-determine with good accuracy the peak structure of the one-loop matrix elements, thus significantly increasing the unweighting efficiency). This results in a fairly large number of calls to the matrix elements of eq. (2.1) per unweighted event and, given the computing performances of the current version of MadLoop, renders it very time expensive to obtain a good-sized unweighted-event sample (say, O(1 M)). We have therefore opted for an alternative approach. We have generated unweighted events using the matrix elements of the Born process multiplied by the gg parton luminosity (rather than by the uū one that would be required if computing the Born contribution proper). The kinematical configurations so obtained have been used to compute the matrix elements of eq. (2.1), which thus provide the correct event weights. A similar procedure has been employed in a recent study on Higgs production via heavy-quark loops and has been shown to work well even for more complex final states [24].

Scale and PDF uncertainties
Among the parameters which enter a short-distance cross section, renormalization (µ R ) and factorization (µ F ) scales and PDFs play a special role, since their variations are typically associated with the purely theoretical uncertainty affecting observable predictions (which, in the case of PDFs, is not quite correct -one may say that PDF variations parametrize the uncertainties not directly arising from the process under study). It is therefore fortunate that (the bulk of) the time-consuming matrix element computations can be rendered independent of scales and PDFs, as opposed to what happens in the case of other parameters, e.g. particle masses. This is doable thanks to the fact that short-distance cross sections can be written as linear combinations of scale-and PDF-dependent terms, with coefficients independent of both scales and PDFs; it is thus possible to compute such coefficients once and for all, and to combine them at a later stage with different scales and PDFs at essentially zero cost from the CPU viewpoint. The crucial (and non-trivial) point is that the noteworthy structure mentioned before is a feature not only of the parton-level LO and NLO cross sections, but also of the MC@NLO ones. This implies that from the conceptual point of view the same procedure for determining scale and PDF uncertaintities can be adopted in MC@NLO as in LO-based Monte Carlo simulations; for the former we shall simply need to compute a larger number of coefficients than for the latter. In order to illustrate what is done in aMC@NLO, we start from dealing with LO computations, for which the notation is simpler. Here and in what follows, the expressions for the short-distance cross sections are taken from ref. [23]. The fully-differential cross section is: Here, dχ Bj denotes the integration measure over the Bjorken x's, and J Bj is the (possibly trivial) corresponding jacobian factor. The quantities on the r.h.s. of eq. (2.4) are the scattering amplitude squared M (n,0) (including coupling constants, and colour/spin average and flux factors), a set of cuts J n (B) L which prevent phase-space singularities from appearing at this perturbative order, an average factor N for identical final-state particles, and the n-body phase space dφ n . We shall write the latter as follows: where dχ n is a measure that understands the choice of 3n − 4 independent integration variables, and K n (χ n ) denotes the four-momentum configuration associated with a given choice of the latter. Finally, Φ (n) is the phase-space factor (including jacobian) that arises from the choice of dχ n . One can then rewrite eq. (2.3) as follows: and, consistently with ref. [23], b is the power of α S implicit in M (n,0) . By construction, the quantity w (B,n) defined in eq. (2.7) is independent of scales and PDFs, and is the first example of the coefficients mentioned at the beginning of this section. When integrating eq. (2.3), one obtains the N -event set after shower. In eq. (2.11), MC(K n;i , x 1;i , x 2;i ) denotes the complete final-state configuration obtained by showering the hard subprocess. When making alternative choices for scales and PDFs, say µ ′ R , µ ′ F , and f ′ i , one can simply repeat the procedure outlined above; this is straightforward, but extremely time consuming if the full scale and PDF uncertainties have to be determined. Alternatively, one can compute eqs. (2.10) and (2.11) by performing the replacement where The denominator here is computed using eq. (2.6), i.e. with the original choices of scales and PDFs. As we shall show in the following, an analogous equation will hold for both the NLO and MC@NLO cases. From eq. (2.6) we obtain: , (2.14) which shows explicitly that the computation of R i does not entail any matrix-element calculation. Therefore, after performing the bulk of the calculation, i.e. the determination of the set in eq. (2.8), one can fill a "central" histogram with the weights defined in eq. (2.9), and as many "variation" histograms as one likes with the weights that appear on the r.h.s. of eq. (2.12), by simply recomputing the factors R i for all the different scales and PDFs needed. The procedure outlined above is exact when applied to eq. (2.10), but it implies an approximation in the case of eq. (2.11). This is because the quantity MC(K n;i , x 1;i , x 2;i ) does contain an implicit dependence on f i , and these functions are not replaced by f ′ i when using eq. (2.12) (since such a replacement would imply performing the shower for each new choice of PDFs). However, this approximation is usually a very good one (barring perhaps the corners of the phase space, and in particular the large-rapidity regions). This is empirically well known, and is due to the fact that the Sudakov form factors used in the backward evolution of initial-state partons are sensitive to the ratios of PDFs, whose variation is much smaller than that of their absolute values (the more so within a PDFerror set, which is the typical application of the procedure discussed here). Furthermore, it should be clear that both scale and PDF uncertainties are defined using well-motivated but ultimately arbitrary conventions, and therefore are not quantities that can be computed with arbitrary precision.
We now turn to discussing the case of NLO short-distance cross sections. Equation (2.6) is generalized as follows: where α = E, S, C, and SC correspond to the contributions of the fully-resolved configuration (the event), and of its soft, collinear, and soft-collinear limits (the counterevents) respectively. The quantities W (α) will be written as follows: where Q is the Ellis-Sexton scale; in eq. (2.17), this scale (which is extensively used in the manipulation of the one-loop contribution -see ref. [23] for a discussion in the context of MadFKS) offers a convenient way to parametrize the dependence on the factorization and renormalization scales. As the notation in eqs. (2.16) and (2.17) suggests, the values that such scales assume in the event and counterevents may be different from each other (however, they must tend to the same value when considering the relevant infrared limits). Furthermore, eq. (2.16) takes into account the fact that, when using event projection to write an NLO partonic cross section (see ref. [1]), the Bjorken x's of the event and counterevents need not coincide. The last term on the r.h.s. of eq. (2.17) is the Born contribution, which in MadFKS can be integrated simultaneously with the other terms that enter an NLO cross section; having a soft-type kinematics, it is naturally associated with the soft counterevent. The coefficients W introduced in eq. (2.17) are the analogues of w (B,n) defined in eq. (2.7) -they are scale-and PDF-independent; their explicit forms are given in appendix A. As far as W B is concerned, it is equal to w (B,n) , up to a normalization factor (due to the different integration measures in eq. (2.6) and (2.16)).
The integration of the NLO cross section leads to the set of N weighted events 2 : Equation (2.10) is generalized as follows: As is well known, for any given i the weights defined in eq. (2.19) may diverge, but the sum in eq. (2.20) is finite for any infrared-safe observable. When changing scales and PDFs, we can now adopt the same procedure introduced before for the LO cross sections, and define the rescaling factors: The quantities R (α) i thus defined do not factorize the coefficients W , and hence have a more complicated form than their LO counterpart, eq. (2.14). Nevertheless, precisely as in the case of eq. (2.14), their computations only require the evaluation of PDFs and coupling constants for each new choice of PDFs and scales, the idea being that the coefficients W are computed once and for all when integrating the partonic cross section, and are stored for their later re-use in eq. (2.21).
We have also understood that the formulae presented above apply to the contribution due to one given FKS pair, i.e. they are proportional to a given S function S ij . We remind the reader that while the observable cross section is obtained by summing over all Sfunction contributions, these are fully independent from each other, and can be effectively treated as independent partonic processes.
We finally discuss the case of MC@NLO. The short-distance cross sections can be written as follows [1]: with the Monte Carlo counterterms given by: The index c in eqs. (2.22) and (2.23) runs over all particles which are colour-connected with the particle that branches. The latter is identified with the parent particle of the FKS pair (i, j) which determines the S function S ij implicit in dσ (NLO,α) . This is sufficient to render eqs. ( (2.25) whose kinematic parts are: The final-state kinematic configurations relevant to H and S events are: where we have used the fact that in MadFKS and aMC@NLO the phase-space parametrization is chosen in such a way that the kinematics of the counterevents coincide (up to the irrelevant unresolved partons; see ref. [23]). The weights that appear in eq. (2.25) are: Here, the same remark as made in the case of LO-showered predictions applies, that concerns the PDF dependence implicit in MC(E H;i ) and MC(E S;i ). There is however a further subtlety in the case of MC@NLO. When showering the Born contribution to S events, one cancels the MC counterterm contribution to H events. In the case of initial-state emissions, we can schematically write the contribution to the H events obtained by showering the Born as follows: and therefore the cancellation of the PDFs with argument x 1 does not take place any longer.
Since changing the PDFs in the shower would amount to replacing f 1 with f ′ 1 in eq. (2.36), eq. (2.36) implies a fractional difference w.r.t. the use of f ′ 1 in both the short-distance and shower computations equal to: which is a number generally fairly close to one. In fact, what is discussed here is precisely the situation one faces when running MC@NLO and adopting different PDFs in the shortdistance calculation and shower phases. This is rather commonly done, and no significant differences are observed w.r.t. the results obtained with the strictly correct procedure of using the same PDFs in the two phases. While this approximation may be undesirable when very accurate results are mandatory, it is always acceptable when computing the PDF uncertainty, which is by definition not a precision result. Furthermore, it should be kept in mind that the procedure currently used for determining PDF "error" sets is that of employing parton-level NLO cross sections; hence, such PDFs are not guaranteed to serve their purpose for those observables for which matched predictions differ significantly from the NLO ones (i.e., where resummation effects are important).

Phenomenological results
In this section we present some sample results obtained with aMC@NLO, its LO counterpart (aMC@LO), and for parton-level NLO (O(α S )) and LO (O(α 0 S )) cross sections. We also give showered predictions for the gg O(α 2 S ) contribution. In both HERWIG [25] and Pythia, we have switched QED showers off. We consider the production of the following leptonic final states:  Table 1: Settings of physical parameters used in this work, with dimensionful quantities given in GeV.
at the 7 TeV LHC. The only cut at the generation level is: which is applied to all unlike-sign lepton pairs. This cut, that serves the purpose of avoiding the zero-virtuality divergences of diagrams with resonant photons, is obviously not necessary in the case of eq. (3.1) for e ± µ ∓ pairs, but it is nonetheless applied there in order to perform a comparison on equal footing between the e + e − µ + µ − and e + e − e + e − final states. The parameters used in this paper are reported in table 1, where the values of α S have been chosen as prescribed by the PDF sets we have adopted, MSTW2008(n)lo68cl [26] (with their associated "error" sets). The default renormalization and factorization scales are set equal to the invariant mass of the four leptons, When studying scale dependences we shall vary µ R and µ F independently, i.e. we shall set (µ R , µ F ) = (κ R µ 0 , κ F µ 0 ) and consider the combinations (κ R , κ F ) = (1, 1), (1/2, 1/2), (1/2, 1), (1, 1/2), (1, 2), (2, 1), (2,2). The overall scale uncertainty will be the envelope of all individual results. We begin by presenting in table 2 the fully inclusive cross sections (bar the invariantmass cuts of eq. (3.3)). For the NLO results and the gg contribution we also give the scale (first error) and PDF (second error) uncertainties. The scale dependence at the NLO is fairly small, which is in part due to the fact that the qq parton luminosity is almost µ F -independent in the Bjorken-x range that gives the dominant contribution to the cross section. This implies that at the LO (which, being O(α 0 S ), is independent of µ R ) the scale uncertainty is comparable or smaller, and thus not particularly meaningful: we refrain from quoting it explicitly here. The PDF uncertainty is estimated following the prescription given by the PDF set [26] (asymmetric Hessian); at the LO, we obtain a fractional uncertainty similar to that of the NLO result. On the other hand, the gg channel displays the typical scale dependence of an O(α 2 S ) LO cross section, the figures in the table receiving the dominant contributions from the µ R dependence. The scale uncertainty is largely dominant over the PDF one. Not surprisingly, the ratio between the two leptonic channels considered in this paper is equal to two up to a very small correction (only slightly larger than the integration uncertainty), since differences arise in the offresonance regions (where singly-resonant diagrams are more important than elsewhere)

Cross section (fb)
Process qq/qg channels gg channel  Table 2: Total cross sections for e + e − µ + µ − and e + e − e + e − production at the LHC ( √ S = 7 TeV), for the production channels considered in this paper, and within the cuts given in eq. (3.3). The first and second errors affecting the results are the scale and PDF uncertainties (also given as fractions of the central values). See the text for details. whose contribution to the total integral is small. Given that the gg matrix elements do feature only doubly-resonant diagrams, we have not computed the result relevant to the e + e − e + e − channel. Finally, we point out that the entries of table 2 have been compared with their analogues obtained from MCFM (also for the non-central scale values, which is a check of the reweighting technique discussed in sect. 2.1), and complete agreement has been found.
As was mentioned above, the scale uncertainties reported in table 2 are the envelopes of the results obtained with the seven combinations of (µ R , µ F ) considered here. It is obvious that this is not equivalent to variations in the range µ 0 /2 ≤ µ R , µ F ≤ 2µ 0 if the dependence on either scale is not monotonic. While the study of continuous variations is clearly impossible, by means of reweighting one can probe the relevant scale ranges with more than three values for each scale, at a very modest CPU cost. In fact, each choice of (µ R , µ F ) corresponds to filling an extra set of histograms (each set containing all observables of interest), with the same kinematics as for the central set and with rescaled weights, which is a negligible fraction of the generation-plus-shower computing time.
We now turn to discussing the case of differential distributions, and we start by considering the O(α 0 S ) and O(α S ) contributions. We shall later address the (generally very small) differences between HERWIG and Pythia6, and here we shall use only HERWIG in order to be definite. The aim of figs. 1-3 is that of assessing the perturbative behaviour of our predictions. In the main frame, these figures present the results for the e + e − µ + µ − channel obtained with aMC@NLO (solid black), NLO (dashed red), aMC@LO (solid blue), LO (dashed magenta). The middle insets show the scale (dashed red) and PDF (solid black) uncertainties, both with aMC@NLO, and given as ratios over the central values. Finally, the lower insets present the ratios dσ(e + e − e + e − ) dO dσ(e + e − µ + µ − ) dO , (3.5) with O any of the observables considered, as computed by aMC@NLO. Figure 1 displays two observables constructed with the sum of the four-momenta of the four leptons, the invariant mass (left panel) and the transverse momentum (right panel). These have very different behaviours w.r.t. the extra radiation provided by the parton shower, with the former being (almost) completely insensitive to it, and the latter (almost) maximally sensitive to it. In fact, the predictions for the invariant mass are basically independent of the shower, with NLO (LO) being equal to aMC@NLO (aMC@LO) over the whole range considered. The NLO corrections amount largely to an overall rescaling, with a very minimal tendency to harden the spectrum. The four-lepton p T , on the other hand, is a well known example of an observable whose distribution at the parton-level LO is a delta function (in this case, at p T = 0). Radiation, be it through either showering or hard emission provided by real matrix elements in the NLO computation, fills the phase space with radically different characteristics, aMC@LO being meaningful at small p T and NLO parton level at large p T -aMC@NLO correctly interpolates between the two. The different behaviours under extra radiation of the two observables shown in fig. 1 is reflected in the scale uncertainty: while in the case of the invariant mass the band becomes very marginally wider towards large M (e + e − µ + µ − ) values, the corresponding effect is dramatic in the case of the transverse momentum. This is easy to understand from the purely perturbative point of view, and is due to the fact that, in spite of being O(α S ) for any p T > 0, the transverse momentum in this range is effectively an LO observable (the NLO effects being confined to p T = 0). The matching with shower blurs this picture, and in particular it gives rise to the counterintuitive result where the scale dependence increases, rather than decreasing, when moving towards large p T [18]. Finally, the lower insets of fig. 1 display the ratio defined in eq. (3.5) which, in agreement with the results of table 2, is equal to one half in the whole kinematic ranges considered. The only exception is the small invariant mass region, where off-resonance effects become relevant. A gauge-invariant way to suppress off-resonance effects, and to select doubly-resonant contributions, is that of imposing: on all equal-flavour lepton pairs; we call the cut of eq. (3.6) the Z-id cut. Lepton pairs that pass the Z-id cut are called Z-id matched, and can be roughly seen as coming from the decay of a (generally off-shell) Z boson. While in the case of the e + e − µ + µ − channel there is only one way to choose two same-flavour lepton pairs, there are two different pairings in e + e − e + e − production. In the case both of these pairings result in lepton pairs that fulfill eq. (3.6), we choose that with the smallest pair invariant mass, and assign the Z-id matched pairs according to this choice; in practice, this is a rare event. By imposing the Z-id cuts the M (ℓ + ℓ − ℓ (′)+ ℓ (′)− ) distribution falls steeply below threshold and gets no contributions below 160 GeV.
In fig. 2 we present two transverse momentum distributions, relevant to the positivelycharged leptons (left panel), and to same-charge lepton pairs (right panel); hence, there are two entries in each histogram for any given event. These results are obtained by applying the Z-id cuts, but we have in fact verified that without such cuts we obtain exactly the same patterns. In the case of the p T of the individual lepton, the aMC@NLO (aMC@LO) prediction is fairly close to the NLO (LO) one, but tends to be slightly harder, owing to the extra radiation generated by the shower. This effect is more pronounced at the LO than at the NLO, which is the sign of a behaviour consistent with perturbation theory expectations. In fact, at the LO all hadronic transverse momentum is provided by the shower, while at the NLO this is not the case; therefore, at the NLO the shower will have less necessity to "correct" the prediction obtained at the parton level, a tendency which is naturally embedded in a matching prescription such as aMC@NLO. The scale dependence is quite small over the whole range in p T , but tends to grow larger towards larger p T . This effect has the same origin as that observed in the right panel of fig. 1, but it is much more moderate than there. This is due to the fact that in the present case the whole range in p T is associated with complete NLO corrections. The PDF uncertainty is seen to be similar to or slightly smaller than that due to scale variation; parton densities are well determined in the x range probed here. Finally, there is no difference between the two leptonic channels for this observable; as already mentioned above, this conclusion is independent of whether one applies the Z-id cuts. The p T of the lepton pairs shown in the right panel of fig. 2 follows the same pattern as the one we have just discussed, but the differences between the various predictions are larger in this case. In particular, aMC@LO is closer to NLO than to LO, which is a consequence of the more important role played by extra radiation in this case (as one expects, the present one being a correlation between two particles rather than a single-inclusive observable). Again, the closeness of NLO and aMC@NLO results shows the desired perturbative behaviour. The more significant impact of extra radiation on this variable is reflected in the slightly larger scale dependence at large p T 's w.r.t. what happens for the transverse momentum of the individual leptons discussed before. The two leptonic channels agree well, also when removing the Z-id cuts. Figure 3 shows two observables constructed after applying the Z-id cuts, namely the pseudorapidity of lepton pairs with opposite charge which are also Z-id matched (left panel; this is then the pseudorapidity of would-be Z bosons), and the azimuthal distance between leptons of opposite charge which are not Z-id matched (right panel; thus, these are leptons emerging from different would-be Z bosons). As in the case of fig. 2, there are two entries in each histogram for any given event. These two observables are dominated by small transverse momenta, and therefore it is not suprising that, at both O(α 0 S ) and O(α S ), the predictions are quite independent of whether a shower is generated or not. Slight differences can be seen in the case of the ∆φ distribution, which is indeed known to be more sensitive than pseudorapidity to extra radiation. The small-p T dominance ensures that scale and PDF uncertainties are flat over the whole kinematic ranges, and of the order of those relevant to total cross section. We now discuss the impact of the O(α 2 S ) gg channel on our predictions. The argument for considering such a channel, despite its being of the same perturbative order as all other NNLO contributions which cannot be included, is the dominance of its parton luminosity over those of the qq and qg channels. This dominance grows stronger with decreasing final-state invariant masses, and hence the O(α 2 S ) versus NLO comparison is significantly influenced by the cut in eq. (3.3) -by lowering such a cut, the relative importance of the gg contribution will grow bigger than the 5%-ish reported in table 2. We also discuss in the following the differences that arise when matching our calculation to Pythia6 rather than to HERWIG. We remind the reader that, depending on input parameters, Pythia is rather effective in producing radiation in the whole kinematically-accessible phase space. This is not particularly useful in the context of a matched computation, where hard radiation is provided (in a way fully consistent with perturbation theory) by the underlying realemission matrix elements. Therefore, we have set the maximum virtuality in Pythia equal to the four-lepton invariant mass. For consistency, this setting has been used also when showering the gg-initiated contribution. Figures 4, 5 and 6 present the same observables as figs. 1, 2 and 3 respectively. In the main frame, we show the aMC@NLO predictions plus the gg contribution (including shower), as resulting from HERWIG (solid black) and Pythia (dashed blue) -we shall The HERWIG and Pythia results for the four-lepton invariant mass shown in the left panel of fig. 4 are identical in the case of the gg contribution, and practically so in the case of aMC@NLO. This is because both event generators fix the invariant mass of the final-state system when doing initial-state showers; hence, any difference in the fourlepton invariant mass can only arise from the real-emission contribution to aMC@NLO, where the final-state system does not coincide with the four leptons. The gg channel has a softer distribution than the qq + qg ones, since by moving towards large masses one probes larger values of the Bjorken x's, where quark parton luminosities are relatively more important than at threshold. In the small-invariant-mass region, the gg contribution does not exhibit any peak structure, owing to the lack of singly-resonant diagrams. Ultimately, the only non-negligible effects are around the peak region, where they are larger than 5%. The four-lepton transverse momentum is marginally softer with Pythia6 than with HERWIG, an effect that is more pronounced in the gg contribution (because of a larger amount of radiation there w.r.t. the case of the qq and qg channels, radiation being on average softer in Pythia6 than in HERWIG). On the other hand, this larger amount of radiation implies that the gg channel is significantly harder than aMC@NLO, and is about 10% of the aMC@NLO+gg result at p T = O(100 GeV). The scale uncertainty shows a behaviour typical of an LO computation, and is vastly different from that of fig. 1; it is constant to a good approximation, which is due to a decreasing µ R dependence with increasing p T , compensated by an increasing µ F dependence.
As was already pointed out in the case of fig. 2, the observables presented in fig. 5 are quite insensitive to extra radiation, and especially so for the single-inclusive lepton p T . Therefore, the similarity between the predictions obtained with HERWIG and Pythia has to be expected. To a very good extent, the gg contribution has the same shape as the aMC@NLO one. Things are marginally different in the case of the pair transverse momentum presented in the right panel; however, any difference between HERWIG and Pythia is not significant, being much smaller than the scale uncertainty of the aMC@NLO contribution alone (shown in fig. 2). Likewise, the shape change induced in the aMC@NLO prediction by the gg contribution is very minor.
The small-p T dominance for the observables displayed in fig. 6 implies that again the HERWIG and Pythia results are very similar. As far as the gg contribution is concerned, its shape is almost identical to that of aMC@NLO in the case of ∆φ. On the other hand, the lepton-pair pseudorapidity turns out to be significantly more central than that of the qq + qg channels; the impact on the overall aMC@NLO+gg shape is however modest for the invariant-mass cuts considered here.

Conclusions
In this paper, we have used the aMC@NLO framework to study the hadroproduction of four leptons, by considering, in particular, the process mediated by the exchange of two Z or γ vector bosons. Di-bosons final states are relevant in both Higgs searches and newphysics studies, where precise predictions are an important ingredient to cope with the overwhelming irreducible Standard Model backgrounds. We have included in our calcula-tion all contributions which are expected to be important for accurate phenomenological analyses, such as singly-resonant diagrams, exact spin correlations and off-shell effects, Z/γ * interference, and gg-initiated subprocesses. Our computation being fully automated, any four-lepton final states can be studied in the very same way, including those that feature W/Z interference effects such as e + e − ν eνe .
In any realistic study, precise calculations are as important as the determination of the associated theoretical uncertainties. For this reason, we have also given in this work practical applications of the fully automatic procedures implemented to this aim in aMC@NLO, namely: i) The matching of the hard process to both Pythia6 and HERWIG, paving the way to systematic studies of the influence of different showers on any observable with MC@NLO 3 ; ii) A framework for studying both scale and PDF uncertainties by a simple (and CPUwise costless) reweighting procedure of the unweighted-event samples.
As for item i), this paper provides the first working example of matching with Pythia6 for a kinematically non-trivial process. As for item ii), we stress that our procedure allows one to obtain precise results with very limited statistics. In fact, the central results and all its variations are correlated, since the same events (up to the weights) are used. This is manifest in the smoothness of the scale and PDF fractional uncertainty plots presented in this paper, which may be compared with, for example, those presented in ref. [18], which had been obtained by rerunning the code. As far as our phenomenological results are concerned, they can be summarized as follows. The impact of NLO corrections on total rates is +40%; for differential distributions, they generally correspond to an overall rescaling, but there are cases where they give nontrivial kinematics effects. The uncertainties due to PDFs are of the order of 2%. The scale dependence we obtain by varying µ R and µ F independently is also of the order of 2% for the qq/qg channels, and of the order of 20% in the case of the gg-initiated subprocess. Parton-level NLO and fully showered NLO distributions are generally in good agreement, except in the very few cases in which the former is not expected to give sensible predictions. These findings are rather independent of possible Z-identification cuts applied to enhance doubly-resonant contributions. Generally speaking, tiny differences are observed between HERWIG and Pythia, with the exception of distributions dominated by resummation effects, where they are more pronounced. The contribution of the gg channel is of the order of 5% of the total with a cut of 30 GeV applied to unlike-sign lepton pairs.
We have explicitly verified that the inclusion of the gg → H → ZZ contribution (not considered for producing the results presented here) is straightforward, and does not increase significantly the running time. This opens the way to simulating Higgs production with four-lepton decay channels by including the signal/background interference. Clearly, for this to be phenomenologically sensible, NLO corrections to the gg → H signal should to be included as well, which could be achieved in several ways, with different levels of approximation and automation [27,28,29].
As a last remark, we point out that ready-to-shower four-lepton NLO (i.e., resulting from qq/qg-initiated processes) event files are available at the aMC@NLO web page http://amcatnlo.cern.ch, for all possible (massless) lepton flavour combinations 4 .

Acknowledgments
This research has been supported by the Swiss National Science Foundation (SNF) under contracts 200020-138206 and 200020-129513, by the IAP Program, BELSPO P6/11-P, the IISN convention 4.4511.10, by the Spanish Ministry of education under contract PR2010-0285, and in part by the US National Science Foundation under Grant No. NSF PHY05-51164. F.M. and R.P. thank the financial support of the MEC project FPA2008-02984 (FALCON). R.F. and R.P would like to thank the KITP at UCSB for the kind hospitality offered during the completion of this paper.

A. Definition of W coefficients for NLO cross sections
We collect in this section the expressions of the coefficients that enter eq. (2.17). These are taken from ref. [23], whose equation (n) will be denoted by eq. (I.n) here; the reader can find fuller details in that paper. The master equation for the determination of all the coefficients W (α) is eq. (I.C.14).

• Event
Only the (n + 1)-body cross section contributes to the event. Therefore, from eq. (I.4.29) and eq. (I.6.7) one gets: withΦ (n+1) = Φ (n+1) /ξ i . The quantity defined in eq. (A.5) is the combination of the realemission matrix element, prefactors, and the phase-space factor that is routinely used in the FKS subtraction. We point out that according to our conventions the matrix elements include the coupling constant, hence eq. (A.4) is independent of g S , since we understand that the factorization scale µ R there is computed according to the kinematic configuration K n+1 . The coupling-constant dependence has been factored out in eq. (2.17).
• Collinear counterevent The coefficient W (C) receives contributions from the (n + 1)-body and the degenerate (n + 1)-body cross sections. The latter can be read from eq. (I. 4.41) (to be definite, we assume that the FKS sister is the parton coming from the left). One obtains: The normalization factor F Φ in eqs. (A.6) and (A.7) is due to the fact that the degenerate (n + 1)-body contribution is integrated together with the pure (n + 1)-body one, but is originally defined with a quasi-n-body measure (see eq. (I.4.41)).
• Soft counterevent, n-body contribution, and Born The coefficient W (S) receives contributions from the (n+1)-body and the n-body cross sections. The latter is the sum of the Born (eq. (I.4.4)), collinear (eq. (I.4.5)), soft (eq. (I.4.12)) and finite virtual (eq. (I.4.14)) contributions, plus the last term on the r.h.s. of eq. (I.C.14), which we shall call the RGE term here. The expression for the finite part of the virtual contribution that includes the full scale dependence must be read from eq. (I.C.13). In ref. [23] we have used the Ellis-Sexton scale wherever possible, which now renders it easy to extract the factorization and renormalization scale dependences -more specifically, these dependences are present only in one term in eq. (I.4.5), in one term in eq. (I.C.13), and in the RGE term. This implies: is the analogue of that introduced in eq. (A.6), and its form can be found e.g. in eq. (I.6.14) (note that it includes an S ij factor). In eq. (A.9), by setting µ = Q one eliminates the first term on the r.h.s. of eq. (I.4.6); this is included here in eq. (A.10).
• Soft-collinear counterevent Finally, the results for W (SC) can be obtained by computing the soft limit of W (C) , and by taking into account the fact that the original equations contained plus distributions. One obtains: