Hadroproduction of $W^+ W^- b \bar{b}$ at NLO accuracy matched with shower Monte Carlo programs

We present the computation of the differential cross section for the process $pp(\bar{p}) \to (W^+\,W^-\,b\,\bar{b} \to)\;e^+\,\nu_e\,\mu^-\,\bar{\nu}_\mu\,b\, \bar{b}+X$ at NLO QCD accuracy matched to Shower Monte Carlo (SMC) simulations using PowHel, on the basis of the interface between HELAC-NLO and POWHEG-BOX. We include all resonant and non-resonant contributions. This is achieved by fully taking into account the effect of off-shell t-quarks and off-shell W-bosons in the complex mass scheme. We also present a program called DECAYER that can be used to let the t-quarks present in the event files for $pp(\bar{p}) \to {t\,\bar{t}\,X}$ processes decay including both the finite width of the t-quarks and spin correlations. We present predictions for both the Tevatron and the LHC, with emphasis on differences emerging from three different $W^+\,W^-\,b\,\bar{b}$ hadroproduction computations: (i) full implementation of the $p\,p(\bar{p}) \to W^+\,W^-\,b\,\bar{b}$ process, (ii) generating on-shell t-quarks pushed off-shell with a Breit-Wigner finite width and decayed by DECAYER, and (iii) on-shell t-quark production followed by decay in the narrow width approximation, as described by the SMC.


Introduction
Accurate predictions for the production of tt-pairs alone or in association with some hard objects, such as jets, vector and/or scalar bosons are important for many experimental studies at hadron colliders both aiming at better understanding of the Standard Model and searches for new physics. However, the t-quarks and heavy bosons decay quickly and the detectors detect their decay products. The experiments often concentrate on leptons because the leptonic channels offer a much cleaner final state than the hadronic ones. Thus it is important not only to predict cross sections for the production of the heavy quarks and bosons, but also for the spectra of the leptons that emerge in their decays.
There are many different approximations to predict such lepton spectra. The most precise way of treating spin correlations and off-shellness of heavy particles is to use matrix elements that are obtained from all Feynman graphs belonging to the leptonic final state. The corresponding graphs include both resonant and non-resonant ones. This procedure in principle can be implemented in a leading-order (LO) computation automatically using tools like MadGraph4/MadEvent [1] or CompHEP/CalcHEP [2] available for several years now. However, in practice, the complexity of the computation can become forbidding very quickly. For instance, for the simple case of tt-pair production, if the decays are included, the final state already contains six particles due to the t → W + b → + ν b decay chain. Thus LO tools like ALPGEN [3] and HELAC [4] written on the basis of recursion relations, instead of Feynman graphs, and improved, more efficient, generation of Feynman graphs, like in Madgraph5 [5], have been developed to better deal with final states with a large number of particles. These difficulties are magnified if one wishes to perform the same computations at the next-to-leading order (NLO) accuracy. For the process mentioned above, one has to deal with six-point integrals, which makes the computation numerically cumbersome.
The computation can be simplified if performed in the narrow-width approximation (NWA), with only resonant contributions kept. The NWA is a good approximation if the width of the decaying particle is indeed narrow, such as the case of decaying electroweak vector bosons, or t-quarks that decay before hadronization. Furthermore, selection cuts are often chosen to enhance the resonant contributions as compared to the non-resonant ones. Predictions for many important production channels were computed at the nextto-leading order (NLO) accuracy long ago. More recently, the decay of the t-quarks was included in the NWA at NLO accuracy for three processes including spin-correlations [6][7][8][9][10]. A further simplification can be made, called decay chain approximation (DCA), if the resonant graphs are replaced by on-shell production times decays of the heavy particles, neglecting both off-shell effects and spin correlations between production and decay, which may seem a high price for simplicity.
While the NWA is a theoretically well-defined approach, it misses the effects of parton showers and hadronization which are often significant. Another line of research was to make improvements by implementing the production of t-quarks at the NLO accuracy and match those to shower Monte Carlo programs (SMCs) [11][12][13][14]. Within this approach the t-quarks are produced on-shell and the decay is provided by the SMC in the DCA at LO accuracy without any spin-correlation. As spin correlations may affect some lepton spectra significantly, a method was introduced in Ref. [15] to take into account the spin-correlations in the NWA, primarily intended to be applied in matched NLO+SMC computations. In addition to tt-pair production, this method was also implemented to include the effects of spin-correlations in single top [16] and tt+jet production [17], among other processes.
An important goal in this paper is to compare the predictions made using these different approaches within the same NLO+SMC framework. For this purpose we choose the simplest process, tt hadroproduction and the POWHEG method [18,19] as implemented in the PowHel framework [11][12][13][20][21][22] which is based on the HELAC-NLO [23] and POWHEG-BOX [24] programs. The matched NLO+SMC prediction in the DCA and NWA approximations have already been known for some time [25][26][27], we simply implemented those in our framework. The new ingredient of our computation is the implementation of W + W − bb-production within PowHel. This is among the most complex final states for which matched NLO+SMC prediction has been considered so far.
Hadroproduction of W + W − bb is one of the 2 → 4 processes included in the Les Houches 2010 wishlist [28]. It can be considered a background for Higgs boson production in association with a bb-pair and for new physics searches, in cases where cascades of decays of non-SM particles may give rise to leptons and multiple jets in combination with missing energy.
In the Standard Model (SM), most of the contribution to the inclusive cross-section of this process comes from intermediate tt-production, due to the abundance of tt-pairs expected at high-energy hadron colliders and to the fact that t-quarks decay almost exclusively into W b-pairs (with a branching ratio of about 99%). Thus, the easiest way to account for the bulk of this cross-section is to factorize the computation in two parts as in the DCA. Results of computations according to this scheme are available since long time [25] and have been included in Monte Carlo generators like MCFM [29,30]. Effects of spin correlations are also known at the next-to-leading order (NLO) accuracy as well [6].
In order to increase the accuracy of the predictions to the few-percent level, as imposed both by the experimental high-luminosity that can be reached at the LHC and by the need for better understanding and measuring the t-quark properties, a full NLO QCD computation of the p p(p) → W + W − bb process including the effects of off-shell t-quarks, is desirable. Such a computation was presented for the first time in Ref. [31], performed by assigning a width to the t-quarks in the framework of the complex mass scheme. The complex mass scheme was originally introduced at LO accuracy and extended to the NLO accuracy in Refs. [32][33][34]. In this framework, besides doubly-resonant Feynman-graphs even singly-and non-resonant graphs, and also their interference, expected to give a contribution suppressed by powers of O(Γ t /m t ), were accounted for in the cross-section evaluation.
As for W decays, the fully leptonic channel was considered in Ref. [31], by treating the W + → e + ν e and W − → µ −ν µ decays in a spin-correlated NWA. In Refs. [35,36] effects of finite widths of both t-quarks and W bosons were included, also including spin-correlations. The W bosons were treated in a fixed width scheme, since their decay does not receive QCD NLO corrections. One of the main conceptual breakthroughs of these papers was the computation of virtual correction contributions including t-quarks in the loop with finite width, which was achieved by the introduction of complex masses in the one-loop scalar integrals and in the on-shell mass renormalization counterterm. The preservation of gauge invariance was a driving principle to check the reliability of the implementations.
Based on the full computation of the process p p(p) → (W + W − bb →) e + ν e µ −ν µ bb + X at NLO accuracy, as presented in Ref. [35] 1 , in this paper we take a further step by matching the NLO QCD predictions to SMC programs using the POWHEG method. Using an SMC allows for an evolution of the final state to the hadron level. We provide predictions for integrated and differential distributions at this level at both LHC and Tevatron. We explicitly compare the results of the full high-accuracy computation outlined above, with the ones arising in simplified cases, where (a) t-quark production is factorized with respect to their decay, performed either in the DCA, neglecting spin-correlations (as implemented in the SMC), or (b) according to a more advanced approach, where events are produced with on-shell t-quarks, which are pushed off-shell according to a Breit-Wigner distribution in a second step, and then decayed at a mixed LO/NLO accuracy [15], thus including spin-correlations (as implemented in the code Decayer) 2 .
After the implementation of several 2 → 3 processes at NLO+PS accuracy, the first 2 → 4 process, the production of a W + pair in association with two jets, was implemented in the POWHEG-BOX [41,42]. Although seemingly similar, this latter process is much simpler than our case as far as the singularity structure is concerned. The squared matrix element of the leading-order (LO) cross section is finite for W + W + j j production. For the W + W − bb-production, as we consider massless b-quarks, the squared matrix elements are singular even at LO, which necessitates the suppression of such singularities. In this respect our computation is more similar to tt bb-production with massless b-quarks, which was studied at the hadron level by including NLO QCD corrections matched to a SMC program by means of the POWHEG matching procedure [43]. The squared matrix element for this process is also singular already at LO.
The main outcomes of our computations are Les Houches events (LHEs), i.e. partonlevel events in files according to the Les Houches accord [44]. For tt-production such events contain a tt-pair and possible further one parton corresponding to first emission. In the case of tt-production processed by Decayer, as well as for W + W − bb-production the events contain four leptons, a bb-pair (we chose the e + ν e µ −ν µ bb channel) and possible first emission. We perform analysis of these events with some standard selections cuts, but we emphasize that the events are available upon request for almost arbitrary experimental analysis. Our first results on this process on the basis of a preliminary set of events were presented in Ref. [45,46].

W + W − bb-production
An important goal in this paper is to check the precision one can achieve with the decay of heavy particles in the decay-chain and narrow-width approximations. In order to make unquestionable statements, we aim at comparing the distributions obtained with approximations to those obtained with an exact computation at the NLO accuracy, when both resonant and non-resonant t-quarks are taken into account in the matrix elements. This amounts to considering all Feynman-graphs to the p p(p) → W + W − bb process that contribute at the NLO accuracy.
The starting point for our computations is the factorization theorem, which gives the cross section for the collision of hadrons A and B as a convolution, of the parton density functions (PDFs) f p/P (x p , µ 2 F ) (p = a, b and P = A, B) with the partonic cross section for the collision of partons a and b, σ ab (p a , p b ; µ 2 F ). Although we set the b-quarks massless when computing the hard-scattering matrix elements, we neglect b-quarks in the initial state, so the summation over a and b runs over u-, d-, c-, s-quarks and the gluon.
The partonic cross section has the perturbative expansion Figure 1. Sample Feynman graphs for the three types of W + W − bb production at Born level: non-resonant, singly-resonant and doubly-resonant ones.
up to NLO accuracy. The LO contribution, is the integral of the spin-and color-averaged squared Born matrix element |M B | 2 over the phase space Φ B of the six final-state particles (b,b, ,ν ¯ , ν ), multiplied with flux normalization, where s = x a x b S, with S = (p A + p B ) 2 being the square of the hadronic center-of-mass energy. We show three sample Feynman-graphs (the decays of the vector bosons being suppressed) of the three types in Fig. 1: (a) non-resonant, (b) singly-resonant and (c) doubly-resonant W b production. There are 38 graphs for the gg-channel and 14 graphs for the qq-channel (with undecayed vector bosons in the final state). The NLO correction is a sum of two separately divergent contributions, the real and virtual corrections, with finite sum for infrared safe observables. The squared matrix element R = |M R | 2 for the real and 2Re(M * B M V ) for the virtual contributions are obtained from 2 → 7 tree graphs and 2 → 6 one-loop graphs, respectively. Sample Feynman-graphs for the three basic types are shown in Figs. 2 and 3, with vector-boson decay suppressed again. There are 795 one-loop graphs in the gg channel and 294 one-loop graphs in the qq-channel. For the real corrections the qg-channel, leading to an extra quark in the final state, opens too.
The finite sum of the two types of radiative corrections can be rewritten as the sum of three finite terms in the following general form: Figure 2. Sample Feynman graphs for the three types of W + W − bb production at one-loop level. Figure 3. Sample Feynman graphs for the three types of W + W − bb+g production at tree level.
where σ NLO{R} is the regularized real correction (real radiation minus subtraction terms) and σ NLO{V } is the regularized virtual corrections (one-loop corrections plus integrated subtraction terms), written in the usual POWHEG terminology [19]. In the same terminology σ LO in Eq.
. The third term in Eq. (2.4) contains the finite remainders from the cancellation of the -poles of the initial-state collinear counterterms, in the POWHEG terminology usually denoted by G ⊕ (Φ B,⊕ ) and G (Φ B, ), where x is the parton-in-parton momentum fraction in the collinear factorization term. There are standard techniques to compute the finite contributions in Eq. (2.4) discussed in the literature in full detail, therefore, we do not deal with this problem here. The necessary formulas for the NLO and PS matching in the two popular subtraction schemes [47,48] are given in Ref. [19]. In the POWHEG-BOX framework used in the present computation, the FKS subtraction scheme [48] is implemented in a process independent way. To obtain the POWHEG cross section for the generation of an event with first (hardest) emission included, one first parametrizes the real-emission phase space in terms of the Born phase space and three more variables that describe the radiation process, leading to where dΦ rad includes the Jacobian of this change of integration variables. Next one defines the NLO-corrected fully differential cross section belonging to the underlying Born configuration 8) and the POWHEG Sudakov form factor The function k ⊥ (Φ R ) is equal to the transverse momentum of the emitted parton relative to the emitting one. Then the POWHEG fully differential cross section (the cross section that can be obtained from the LHEs) is defined as (2.10) The advantage of this formula is that it can be used to generate equal weight events with Born configuration (first term) or including first radiation (second term). Computing the derivative of the Sudakov form factor with respect to p ⊥ , one can prove the unitarity relation [19] Substituting into Eq. (2.10), one finds that the total POWHEG cross section is equal to that at NLO accuracy, In such cases, instead of checking the total cross section one can but check the consistency of the differential distributions obtained from the POWHEG formula with those at NLO accuracy.
The differential distribution of the observable O can be computed as Using Eq. (2.7) and the explicit form of dσ LHE in Eq. (2.10), we can rewrite it as (2.14) We can again use the unitarity relation in Eq. (2.11) to simplify the first integral, and use the expansions of the Sudakov form factor and theB functions in the strong coupling [19], to obtain the relation where Φ R is parametrized as in Eq. (2.7) and we dropped the Θ(k ⊥ (Φ R ) − p min ⊥ ) function, its effect being suppressed by p min ⊥ [19]. Thus any observable has the NLO accuracy up to higher order corrections multiplying the real emission contribution, where the magnitude of the correction is set by the expansion in Eq. (2.15). Therefore, when the NLO K-factor is large, the distributions obtained from dσ LHE may differ from the NLO cross section significantly (see Ref. [50] for more detailed discussion). It is possible to decrease this difference by writing the real correction as a sum, R = R S + R R , where R S is constrained to near singular regions, hence called singular contribution, while R R is called the remnant [50]. Using such a decomposition, the size of the difference between the POWHEG and NLO cross sections is controlled by R S : qq e −ν e µ + ν µ bb → 0 g g e −ν e µ + ν µ bb → 0 qq e −ν e µ + ν µ bb g → 0 g g e −ν e µ + ν µ bb g → 0 Table 1. Born-level (first row) and real emission (second row) subprocesses.

PowHel implementation
The POWHEG-BOX provides a general framework to compute the POWHEG cross section in Eq. (2.10). In this framework, the following ingredients are needed: • The flavor structures of the Born and real radiation emission subprocesses, listed in Table 1.
• The Born-level phase space, that we generate to emphasize the doubly-resonant kinematics (see below).
• We obtain the squared matrix elements for the Born and the real-emission processes and color-correlated Born amplitudes with all incoming momenta using HELAC-NLO [23] (in particular, HELAC-1LOOP [51], based on the OPP method [52] complemented by Feynman-rules for the computation of the QCD R 2 rational terms [53], and Helac-Dipoles [54]). The matrix elements in the physical channels were obtained by crossing. In order to treat the numerical instabilities, we implemented doubledouble precision (dd-precision means 30 digit accuracy) numerics by developing a HELAC-1LOOP@dd version of the HELAC-1LOOP program.
• We project spin-correlated Born amplitudes from the helicity basis to the Lorentz one by using the polarization vectors.
With this input in the POWHEG-BOX we can generate hadronic events. In principle we can choose any SMC program for generating parton showers, decays of heavy particles and hadronization. There is however, a caveat when choosing the PS. We generate events with hardest emission measured by its transverse momentum. If the ordering variable in the PS is the transverse momentum, such as in PYTHIA p ⊥ -ordered, then the final-state system after first radiation is evolved using the SMC and any emission with transverse momentum higher than that of the first emission is impossible (automatically vetoed shower). If the ordering variable in the PS is different from the transverse momentum of the parton splitting (e.g. angular-ordered showers in HERWIG [55], or virtuality-ordered showers in PYTHIA [56]), then the hardest emission is not necessarily the first one. In such cases these SMC codes discard in subsequent splittings shower emissions with larger transverse momentum than that in the real emission correction (vetoed showers). In addition, a truncated shower simulating wide-angle coherent soft emission before the first emission is also needed in principle, but its effect was found small [57]. As there is no implementation of truncated shower in HERWIG using external LHE event files, the effect of the truncated showers is absent from our predictions. To quantify these differences we make predictions in Sect. 5 with both p ⊥ -ordered and virtuality-ordered PYTHIA, as well as with HERWIG.  We take into account the finite width of the t-quarks and W -bosons in the complex mass scheme [32,34] where the Feynman propagator factor is substituted as The presence of these propagators results in Breit-Wigner factors in the squared matrix elements. Thus it is expected that the dominant contributions come from those configurations where the two-particle kinematic invariant assigned to the off-shell particle, e.g. s t = (p b + p W + ) 2 , is around its pole mass, e.g. s t m 2 t . In order to integrate over these resonances efficiently, we built the phase space emphasizing the doubly-resonant channel as shown in Fig. 4.
The POWHEG-BOX uses MINT as its integrator, which performs well in the absence of resonances. For our case with four invariants having Breit-Wigner peaks in the physical region the integrator fails to set up an integration grid, therefore, we performed the optimization for these variables, described as follows. Let us consider a massive particle with momentum p, mass m, and width Γ, then the kinematic invariant can be written as s = p 2 , hence the Breit-Wigner factor is proportional to (2.20) We generate the phase space such that the integration variable is the invariant s itself, and introduce a new integration variable The Jacobian 1/F BW (s) cancels the factor F BW (s) in the squared matrix element. The limits on ρ are related to the limits on s by For the present calculation we used s min = 25 GeV 2 and s max = s. Finally, we map ρ to ξ ∈ [0, 1] using hence the original integration measure over s is expressed through ξ as After performing the change of the integration variable from s to ξ, the integrand is flat in the new variable, and MINT can produce the integration grid. The generation of the matrix elements is straightforward using HELAC-NLO. During integration we had to solve two more problems. The first one is that for vanishing transverse momentum of the b-quarks or vanishing invariant mass of the bb-pair the Born cross section becomes singular. While this can never happen in a LO computation due to the application of selection cuts, it is a problem in the POWHEG method because the selection cuts can only be applied after event generation. The traditional way of treating this problem is the introduction of generation cuts [43,49]. We use p ⊥, b ≥ 2 GeV for both the b-andb-quark and a cut m bb ≥ 1 GeV. With these cuts the LO cross section becomes finite, but the generation of the events is still rather inefficient because most of the events are generated in the region of small p ⊥, b , and are thus lost when the physical selection cuts (usually much higher, in the region of 20 -30 GeV) are applied. In order to make the generation of events more efficient, we introduce suppression factors [58]. As we want to suppress the region of small p ⊥, b , the natural choice for the suppression is where we use i = 3 in our calculation. The second problem is a purely numerical one and is related to the numerical computation of one-loop amplitudes through CUTTOOLS [59], as implemented in HELAC-1LOOP. In order to control numerical instabilities, we developed a dd-precison version of the program, the HELAC-1LOOP@dd code [43], which is a straightforward extension of HELAC-1LOOP to double-double precision using QD [60].

Production of decaying tt-pairs
The simplest way of implementing the decay of heavy particles in a NLO computation matched with SMC is to generate on-shell particles and then let them decay by the SMC, corresponding to decay in the DCA. This procedure has been criticized strongly [10] because it neglects spin correlations. The other extreme is to include the off-shell effects of the heavy particles throughout the computation, which however can easily lead to very cumbersome computations.

Angular correlations in particle production
In Ref. [15] a method was presented to implement decays in a way such that angular correlations between decay products of a heavy particle and the rest of the event are taken into account. It is an improvement to the DCA with including angular correlations at LO accuracy. This method is based on a hit-and-miss technique through the following steps: 1. Generate LHEs with heavy particles in the final state (tt-pair in our case). We denote the collection of the kinematic variables that characterize such an event with Φ u . Compute the squared matrix element |M| 2 (Φ u ) for the event.
2. For each LHE, generate the four-momenta of the decay products uniformly within the decay phase space of the corresponding parent particle. The corresponding kinematics for the event with decay products in the final state.

Generate a uniform random number
If R ≤ r, then discard Φ d and return to step 2.
4. Otherwise, replace the momenta of the heavy particles by those of their decay products as obtained in step 2 (Φ u → Φ d ).
We implemented this method in a program called Decayer, together with a further improvement which takes into account the off-shell effects of the decaying particles, as described in the next subsection.

Reinstating off-shell effects
The program PowHel generates events with the undecayed phase space Φ u (specified explicitly in Eq. (3.7)). We would like to reweight these events with a Breit-Wigner distribution.
In the following we describe the method implemented in the program Decayer to perform the decay of the on-shell heavy quarks, generated by PowHel, with off-shell effects also taken into account. We denote the momenta of the heavy particles in Φ u by {q i } n i=1 . In addition to the n heavy particles, we also allow for m light particles with momenta {p j } m j=1 in the undecayed event. Our aim is to generate a decayed event, with a corresponding phase i obtained from the momenta of the undecayed particles using the procedure described below.
The phase space Φ d for total incoming four-momentum Q µ can be written in a factorized form, (3.1) In Eq. (3.1) the two-particle phase space measures can be written as where dΩ = dφ d cos θ is the surface element of the unit sphere. Here λ is the triangle function which can be used to express the magnitude of the three-momentum for one of the decay products, |k (t) 1 |, in the rest frame of t (hence the upper index (t)) as In this frame k The squared matrix elements contain a Breit-Wigner factor F BW (s i ) for the propagator of each heavy particle i. In the narrow width approximation Γ i /m i → 0, and we can make the replacement lim (s i = t 2 i ), which puts particle i onto its mass-shell. With this replacement we can perform the integrations over all dt 2 i , and obtain where all q i are on shell (q 2 i = m 2 i ). The phase space on the right hand side of Eq. (3.6) is that in DCA. The first factor corresponds to the undecayed phase space Φ u , with ω(p) = p 2 + p 2 , while the second one is We now turn to the reweighting of the undecayed events with momenta q µ i with a Breit-Wigner function. In going from Φ u to Φ d , we reverse the arrow in Eq. (3.6) by inserting and multiplying by the factor Of course, inserting the first factor is formal. We actually generate the virtualities t 2 i according to the Breit-Wigner distribution, following the procedure described in Sect. 2. Using the virtualities t 2 i , we construct the four momenta by keeping the three-momenta of the original event (q i ) and changing the energy components E i of the heavy particles according toẼ This way three-momentum conservation is maintained, but energy conservation is lost, We reinstate energy conservation by rescaling the momentum fractions of the incoming partons x a and x b uniformly,x where √ s = E CM and √s =Ẽ CM . Thus the total partonic center of mass energy is and energy conservation is recovered. The original PowHel weights include the PDFs in Eq. (2.1). The rescaling of the momentum fractions induces a rescaling of the event weights by a factor A practical application of this method to the case of tt-decays is provided in Appendix A.

Checks
The PowHel+SMC framework allows for predictions at different levels of evolution during the collision. First, to check correct implementation of the computations, we compared our PowHel predictions at the NLO accuracy to those computed in Ref. [35] for the LHC at center-of-mass energy √ s = 7 TeV, and found agreement. In making these comparisons, we adopt the PDG [61] values for the parameters, already used by Ref. [35]: However, when making new predictions we changed the values of m t and Γ t as specified in Sect. 5. The renormalization and factorization scales were chosen equal to m t , and we used the CTEQ6m parton distribution functions from LHAPDF. The strong coupling α s was computed with 2-loop running with Λ 5 = 226 MeV. These are the choices adopted in Ref. [35].
Next, to check the reliability of the matching procedure, we compared our PowHel predictions from the LHEs (i.e. distributions obtained from events including no more than the first radiation emission) with the NLO ones. These checks were performed on all observables studied by Ref. [35], by exactly applying the same system of cuts presented in that paper. In general, the overall agreement turned out to be acceptable, even if slightly worse as compared with the typical values obtained in our previous papers in the study of less complicated 2 → 3 hard scattering processes matched to SMC [11][12][13]. The differences can be ascribed to the larger NLO K-factor, which allows for less agreement between the NLO and LHE predictions (see Eq. (2.17)).
The comparisons are shown in Fig. 5, where the transverse momentum (a) and rapidity of the b-quark (b) and the charged lepton (c, d) are plotted together with the σ NLO /σ LHE ratios shown in the lower inset of each figure. By looking at the figures it is clear that the rapidity from the LHEs overestimates the rapidity at the NLO accuracy by 3-5 %, whereas the shape of the distribution is unaffected. The shapes of the p ⊥ -distributions from the LHEs are slightly distorted (i.e. within a few percent) with respect to the NLO ones (taking into account the statistical uncertainties of both computations). These results are within the band of NLO scale dependence and compatible with the expectation about the perturbative accuracy of LHE predictions given in Eq. (2.17) that allows for increasing difference with increasing NLO K-factor [50,62]. For the present case the inclusive NLO K-factor is close to 1.5 [35], which gives the small but noticeable difference between the NLO and LHE predictions.
In Fig. 5 we show two more plots comparing the LHE and NLO predictions. The first is the missing transverse momentum in Fig. 5.e that is relevant for new-physics searches based on missing transverse energy plus jets and leptons. We see similar level of agreement as for the case of the other transverse momenta up to about 150 GeV. In the hard tail the LHEs give significantly smaller cross section than the prediction at NLO. This difference may be due to the higher order corrections exhibited on the right hand side of Eq. (2.17), or simply to the lower statistics in the tail. The NLO K-factor is very large in this region (in the range of 2-4) [63], therefore, this difference is well within the NLO scale-dependence. Nevertheless, in order to have a quantitative understanding of the accuracy of our predictions for the missing transverse momentum above 150 GeV we would need much more events, which is beyond the scope of our computational resources at present. This lack of quantitative understanding of the accuracy is also true for the tail of the p ⊥ -distribution of the hardest bb-dijet system, not shown here. Finally, we show the ∆R-separation of the charged lepton-antilepton pair in Fig. 5.f. This distribution is sensitive to the spin correlations between production and decay of the charged leptons. In this case, we find that the agreement between the distributions from the LHEs and at the NLO accuracy is very good over the whole kinematic range.

Phenomenology
As mentioned before, we consider and compare three different cases: dσ dp 1. one is the computation of the p p(p) → (W + W − bb →) e + ν e µ −ν µ bb+X process at NLO accuracy, including t-quarks with finite width throughout, as well as W bosons. The W bosons decay into lepton pairs through a propagator involving the complex W mass and the lepton current.
2. another is the computation of the p p(p) → tt process at NLO accuracy. The generated t-quarks are on shell and their decays,t → W −b and t → W + b are computed by the SMC in the DCA neglecting spin correlations. The decays of the W bosons are also performed by the SMC.
3. a third corresponds to the computation of the p p(p) → tt process at NLO accuracy. In this case, the heavy quarks, produced on shell, are pushed off shell and decayed by means of the Decayer program, taking into account spin correlation effects, as described in Sect. 3. The decays of the W bosons are also performed by Decayer.
We make predictions at the hadron level. As mentioned earlier, for making predictions we adopted the value of the t-quark mass m t = 173.2 GeV [64]. Correspondingly, we use Γ NLO t = 1.32 GeV for case 1 and Γ LO t = 1.41 GeV for cases 2 and 3 [65]. As for the SMC, to simulate the PS and hadronization, we used PYTHIA 6.4.28 [56]. We used both the untuned version (denoted by PY1), providing a virtuality-ordered PS, and a further version, tuned to the Perugia 2011 set of values [66], one of the recent tunes, providing a p ⊥ -ordered PS, updated on the basis of recent LHC data (denoted by PY2). We also used HERWIG 6.520 [55] (denoted by HW). In order to simplify the analysis, we kept π 0 s, that can be easily reconstructed from their γγ decay products, stable. All other particles and hadrons were allowed to be stable or to decay according to the default implementation of the SMC. In HERWIG (anti-)muons were also set stable, as in PYTHIA by default.
While the b-quark masses were set to zero in generating the hard-scattering amplitudes, these were kept to their default values implemented in PYTHIA in the SMC evolution to Bhadrons. Quark masses in HERWIG were set to the same values as in PYTHIA default. In the configuration of the SMC, for all other mass and width parameters, including light quark masses, we used the default values already implemented. We use the CTEQ6.6m parton distribution functions from LHAPDF with Λ 5 = 226 MeV and strong coupling α s computed with 2-loop running.
For cases 2 and 3 we can potentially consider all W decay channels (leptonic and hadronic ones). In order to be able to compare our predictions to those given in Ref. [35] (without including any SMC effects), we nevertheless forced the decay into the e + ν e µ −ν µ bb channel, adopting in both cases a W branching ratio in e + and µ − corresponding to that implemented in PYTHIA (B(W + → e + ν e ) = B(W − → µ −ν µ ) = 0.108, also enforced in HERWIG).
All hadronic tracks with |η| < 5 were used to build hadronic jets. We used the antik ⊥ algorithm [67] as implemented in FastJet [68], with R = 0.4 and R = 1.2, to study the effect of different jet reconstruction strategies. While in collider experiments b-jets are reconstructed with finite tagging efficiencies from the tracks pointing towards vertices displaced with respect to the primary interaction vertex, in our theoretical simulations bjets were tagged by means of the information included in the MCTRUTH parameter available in the SMC, tracking back their evolution up to their origin as b-quarks, and b-jet tagging efficiencies were neglected.
In the event generation we apply technical cuts of p ⊥, b , p ⊥,b > 2 GeV and m bb > 1 GeV (see Sect. 2.1 for details). Taking into account that physical cuts should always be well above the technical ones, we consider the following set of physical cuts: 1. Each jet is required to have transverse momentum p ⊥, j > 20 GeV and pseudorapidity |η j | < 5, otherwise it is not counted among the jets.
2. Each of the jets satisfying the 1st condition, to be classified as a b-orb-jet, is required to be b-tagged and have |η b | < 3, due to the geometry of the tracking system.
3. We require at least one b-jet and oneb-jet. 4. Each charged lepton is required to have p ⊥, > 20 GeV and |η | < 2.5, otherwise it is not counted among the leptons.

5.
We require at least one charged lepton and one charged anti-lepton, that are isolated from all jets by requiring ∆R( , j) > 0.4 (1.2) in the azimuthal angle-pseudorapidity plane. If there are more leptons that pass cut 4, those are kept without isolation from the jets.
These cuts present some modifications with respect to those in Ref. [35] providing predictions at the NLO accuracy.
In addition to the three cases listed at the beginning of this section, we can also compare predictions after parton shower (PS) and after full SMC (PS + hadronization) simulations, obtained with either PYTHIA, or HERWIG. This makes many options for comparison. We decided to present four groups of plots: In order to make connections among these groups of plots, we always include predictions at the hadron level generated with PY1. We produced at least 30 distributions for each possible case and each combination of comparisons. In order to be able to compare the various effects, we selected six standard plots that we show for each group:  6. distribution of the azimuthal separation of the hardest isolated positron and muon ∆φ e + µ − .
These fall into three categories: (i) distributions of the b-quarks only, (ii) those of the lepton sector, and (iii) those involving a b-quark and a charged lepton.

Predictions at the LHC
In order to study the effect of the SMC, first we compare predictions for W + W − bbproduction at different levels of theoretical description: obtained from the LHEs, after PS and after full SMC. For the PS and SMC we use PY1. In addition to the selection cuts (1-6) we employ a jet veto on the non b-jets, too. The reason for this jet veto is that in the LHEs there can be at most one extra jet besides the b-andb-jets, while after PS and SMC there are usually many more (less energetic ones). Thus the selection cuts affect the latter much more, decreasing the cross sections significantly, which is more a consequence of the selection cuts than the effect of the PS and hadronization. In Table 2 we show the inclusive cross sections after selection cuts. The 10 % decrease of the cross section after PS (with jet veto) is mainly due to the different effect of the selection cuts when applied at different stages of the evolution of the events, and very similar decrease is observed on the distributions below. The additional 2 % decrease after full SMC however, appears in the distributions very differently. In Fig. 6 we show the distributions of p ⊥,b 1 and p ⊥,e + , while in Fig. 7 we show the pseudorapidity distributions of the same objects. We find that the effect of the parton shower is a fairly uniform decrease, it is in the range of 0-20 %, independently of the observable, and in general it is about 10 %, similarly to the decrease of the inclusive cross section. The effect of hadronization however, depends on the observable strongly. For the p ⊥,b 1 -distribution it can reach 50 % for p ⊥ above 150 GeV. Inspection of the curve shows that the large effect is mainly due to the softening of the jets, a shift of the distribution by about 25 GeV. One important reason for this softening is the decay of the unstable hadrons, which often transforms partonic energy to electromagnetic and missing energy. In the case of the positron transverse momentum such decays are absent, and the effect of hadronization is much smaller, at most 5 % apart from the statistical fluctuations. For pseudorapidity distributions this transformation of energy does not influence the direction of the jet significantly, therefore, the effect of hadronization is negligible except for |η b 1 | > 2.    An interesting observable is the invariant mass of the combination of the hardest b-jet and the hardest isolated positron shown in Fig. 8.a. For decaying on-shell t-quarks into W + + b → e + ν e b, neglecting the masses of all final decay products, we have Thus, at lowest order in tt-production, there is a strict kinematic limit for the invariant mass of the b-quark and the positron at m 2 t − m 2 W 153 GeV, that is quite sensitive to m t , which hints that this distribution is useful to measure m t [69]. For off-shell t-quarks (e.g. in a computation at NLO accuracy) this kinematic limit is smeared, nevertheless there is a sharp fall of the cross section in the fixed order predic- tions. 3 In Fig. 8.a we show that the main effect of the PS and also that of the hadronization is to smear this sharp edge observed in the fixed-order computation, as expected. Apart from this region around 150 GeV, the corrections of the SMC are modest. We show the azimuthal separation between the hardest isolated positron and muon in Fig. 8.b. For this observable, the main effect of the PS is almost unaffected by the hadronization effects, and the effect of PS is also modest, varying slightly around 10 %.
We show the effect of PS and SMC for two more variables: the distribution of the missing transverse momentum, / p ⊥ , and that of p ⊥, b 1 b 2 = (p x,b 1 + p x,b 2 ) 2 + (p y,b 1 + p y,b 2 ) 2 of the hardest b-andb-jets in Fig. 9. For / p ⊥ we find notable changes above 100 GeV where the effect of PS is up to 20 % (resulting entirely from the decrease of the cross section due to the different effect of the selection cuts on the LHEs and after PS), and that of the hadronization peaks at / p ⊥ 150 GeV where it reaches 25 %. For p ⊥, b 1 b 2 the decrease due to the PS is again between 0-20 % for the same reason. However, the effect of the hadronization is large, between 30 and 50 % above 150 GeV, the range which is important in boosted-Higgs searches with a large tt background [70]. Similarly to the case of p ⊥ -distribution of the hardest b-jet this large effect is due to the transformation of the hadronic energy into electromagnetic and missing energy during hadronization.
We can draw the following general conclusions for the given set of selection cuts: • The effect of the parton shower is usually a 0-20 % decrease because the jets become softer and fewer events pass the selection cuts.
• The effect of hadronization is small (less than 10 %) in the hard leptonic sector, except for / p ⊥ . These distributions are determined by the hard-scattering process and the decay of the heavy particles, and not influenced by hadronization effects. In fact, even if other leptons can be emitted at lower energy scales, in particular by hadron decays, dσ dp the isolation criterion is quite effective in disentangling just the positron coming from the W + decay (and the muon coming from the W − decay).
• The effect of hadronization can be up to 50 % in the p ⊥ -spectra in the hadronic sector (and negative) due to the transformation of hadronic energy into electromagnetic and missing energy during hadronization.
• The p ⊥ -spectra are softened by PS and further by hadronization, while the angular spectra are modestly influenced by PS and hardly by hadronization, as expected.
Next we compare distributions for the three cases obtained from the LHEs. As the NWA approximation with NLO decays is known to describe most kinematic distributions fairly precisely (the inclusive cross section in NWA is about 0.5 % smaller than the prediction of the full calculation) [36], this comparison gives information mainly about the importance of the NLO decays as compared to the LO decays in the two approximations: DCA and the one obtained with the program Decayer, that includes the treatment of both spin correlations and off-shell effects. In order to compare the importance of the SMC effects to the effect of including the decays, we also show the full SMC predictions for case 1 computed with PY1. We use the selection cuts (1-6).
As seen in Fig. 10 including the decays gives a good description over the whole p ⊥ -range for the distributions of both p ⊥,b 1 and p ⊥,e + , with Decayer performing slightly better at large p ⊥ . The effect of the decays is in general smaller than the effect of the full SMC. In fact, when comparing Figs. 6 and 10 we see that the effect of the hadronization increases significantly if we allow for non-b jets that are also taken into account in the selection cuts (an effect due to different particle content).
Similar, although less drastic effects can be observed in the pseudorapidity distributions (see Fig. 11), where we see a small effect of the approximate treatment of decays, an almost uniform downwards shift amounting to 5-10 %. We attribute these differences to the missing NLO corrections in the decays. dσ dp [fb] Ratio η e + Figure 11. Pseudorapidity distributions of a) the hardest b-jet and b) the hardest isolated positron from the LHEs for the three cases. The lower inset shows the ratio of the predictions with decay to the W + W − bb-prediction.
The decays of the heavy particles are well described by both DCA and Decayer up to 150 GeV in the m b 1 e + -distribution, see Fig. 12.a. For larger values DCA and Decayer fail, leading to an underestimate up to 50 %, which is a result of the missing singly-and non-resonant contributions mostly (about 40 % [? ]), and of the missing NLO corrections in the decays to much less extent (about 10 %). Also the two approximate predictions differ in the range of 150-200 GeV due to the off-shell effects in Decayer, missing from the DCA. The effect of the full SMC is small (within 10 %) below 100 GeV and above 175 GeV, but can become very large in between, especially at 150 GeV, where the LHE cross section has a sharp drop. As one expects, this sharp drop is smeared by the PS and hadronization (see Fig. 8.a).
The ∆φ e + µ − -distribution, shown in Fig. 12.b, is an example where the differences be- tween the three cases are clearly visible. The prediction by Decayer is within 10 % of the full one, with slightly increasing difference towards the separation by 180 • . Nevertheless, the shapes of the LHE predictions are similar, implying that Decayer gives a fairly good approximation to describe spin correlations. The difference in normalization is due to the missing NLO corrections in Decayer. The DCA approximation has a different shape due to the absence of spin correlations. Similar pattern as seen between the DCA and Decayer approximations was already observed in the parton-level calculation of Ref. [71] performed with and without spin correlations. For this distribution the full SMC decreases the cross section by about 30 % almost uniformly.
We conclude that the predictions with decays give the shapes correctly, the NLO corrections in the decays in general cause only a uniform increase up to 10 % except for distributions involving the charged lepton emerging from the decay of the W -bosons (usually the hardest isolated charged lepton). The DCA does not describe the shape of the ∆φ e + µ −distribution due to the lack of spin correlations and both approximations fail in the hard tail of the m b 1 e + -distribution due to the lack of singly-and non-resonant contributions.
We now turn to make predictions at the hadron level that are more interesting from the experimental point of view. For this kind of predictions the selection cuts (1-6) were applied after shower, hadronization, hadron decay and the application of jet algorithms. For the three cases, the integrated cross-sections after cuts are collected in Table 3, for different jet sizes (anti-k ⊥ with R = 0.4 versus anti-k ⊥ with R = 1.2). The cross-section in case 1, the most accurate one, is larger than the cross-section in case 2 and 3 by ∼ 10% due to the NLO accuracy in the decays and, to less extent, to the non-resonant contributions. The cross-section decreases in all cases significantly by using a larger jet radius because we also use a more severe jet-lepton isolation with R = 1.2. By comparing cases 2 and 3 it turns out that the effect of spin-correlations in t-quark and W -boson decays increases the R/case 630 ± 4 573 ± 1 582 ± 1 σ(R = 1.2) (fb) 300 ± 3 253 ± 1 261 ± 1 Table 3. Cross-sections at the hadron level at the LHC after cuts (1-6) for the three cases as a function of the R parameter of the anti-k ⊥ algorithm used to identify jets. The quoted uncertainties are statistical only. The predictions are obtained with the PY1 SMC. dσ dp  cross section by a couple of percent because our selection cuts affect the spin-correlated decays slightly differently.
We also compared almost 50 distributions belonging to the three cases, our standard selection is shown in Figs. 13-15. In the inset in the lower part of each figure we plotted the ratio of the cases 2 and 3 (with decays of heavy particles) to the default one (W + W − bbproduction). We see in Figs. 13 and 14 the general trend that Decayer and DCA give very similar predictions both in shape and normalization for p ⊥ -and η-distributions, while the full W + W − bb-computation followed by the SMC differs, but within 10 % of these. These differences are due to the lack of NLO corrections in the decays, and are present already at the LHE level, as discussed above (see for instance Fig. 10.b). The SMC does not change the main features already seen in the LHEs.
There are some exceptions to this general trend, some shown in Figs. 15 and 16. The m b 1 e + -distribution shows a similar pattern as seen in the LHEs (cf. Fig. 12) but the large effect of the singly-and non-resonant graphs above 150 GeV becomes much reduced after SMC, and also the sharp drop in the cross section at 150 GeV is smeared (already seen in Fig. 8.a).
Pseudorapidity distributions of a) the hardest b-jet and b) the hardest isolated positron after full SMC for the three cases. The lower inset shows the ratio of the predictions with decays of the t-quarks in the DCA and Decayer as compared to the complete W + W − bb computation.
Ratio ∆φ(e + , µ − ) Figure 15. Distributions of a) invariant mass of hardest b-jet and the hardest isolated positron and of b) azimuthal separation between the hardest isolated positron and muon after full SMC for the three cases. The lower inset shows the ratio of the predictions with decays of the t-quarks in the DCA and Decayer as compared to the complete W + W − bb computation. the differences between the three cases were clearly visible in the LHEs. These differences are only slightly altered by the PS, or the full SMC. In particular, the effect of including the spin-correlations leads to an increase of the distribution for small azimuthal separation ∆φ e + µ − , where the distribution of case 3 approaches that of case 1, both including spin correlations. These are both 15-20 % larger than the distribution of case 2, where spin correlations are neglected. At large separations however, after SMC the latter becomes even larger than the predictions from case 1.
We discuss two more, closely related, distributions for these physically most interesting predictions. The first one is the p ⊥ -distribution of the hardest non-b jet in Fig. 16.a. Here we observe a difference up to 50 % between the distributions from W + W − bb-production and the other two cases. The scalar sum of the (W + W − bb) transverse momenta for the events that pass the cuts (1-6) at the hadron level, is plotted in Fig. 16.b. Here we see the Sudakov suppression in the low p ⊥ -region clearly, which is different for case 1 and cases 2, 3 for similar reason as in the p ⊥ -distribution of the non-b jet. Clearly, the differences seen in Fig. 16 have nothing to do with spin correlations (cases 2 and 3 give the same) or nonresonant Feynman graphs present in the W + W − bb-calculation (case 1). The probability of emitting the hardest non-b jet at a given transverse momentum from the initial state is similar in the case of tt and W + W − bb production. However, in the case of tt-pair production the probability of a hard jet from the t-quarks is much larger than from the b-quarks in the W + W − bb final state, the latter being dominated by soft and collinear emissions that generally contribute to the b-jet. As a result the p ⊥,j 1 -distribution is much softer for case 1.
Finally, we check the consistency among three different SMC codes: PY1, PY2 and HW in Figs. [17][18][19]. The lower insets in these plots show the ratios of the predictions with PY2 and HW to that with PY1. Concerning the dependence on the SMC code, in general we find that • in most of the phase space, the three SMC codes give predictions within 10 %; • the shapes predicted with the PY2 and HW codes are very similar, but their normalization differs; • the shapes for p ⊥ -distributions obtained with PY2 and HW are slightly harder than those predicted by the PY1 code. The pseudorapidity and angular distributions differ only in normalization.  [fb] Ratio η e + Figure 18. Same as Fig. 17 as for a) the η b1 -distribution and b) the η e + -distribution.

Predictions at the TeVatron
We also studied the same process at the TeVatron. We do not repeat the whole analysis as some of the conclusions are the same or very similar to those drawn for the LHC, only present examples to illustrate the following general observations: • The effect of PS and hadronization on the shapes of the distributions is very similar as found at the LHC, only their sizes are larger, which is expected due to the smaller colliding energies (examples are shown in Fig. 20). As for the inclusive cross section after cuts, we found the cross section values given in Table 4, obtained with PY1 SMC.
• Similarly to the LHC case, the singly-and non-resonant contributions have a small effect (about 1 %), but this time negative (the NWA predictions being larger) [? ].
Having this in mind and looking at Figs. 21 and 22.b, we find that predictions with decays describe the kinematic distributions better at the TeVatron than at the LHC,  Figure 19.
Same as Fig. 17 as for a) the m b1e + -distribution and b) the distribution of the azimuthal separation between the hardest isolated positron and muon. the effect of the missing NLO corrections in the decays being a fairly uniform decrease of about 4 %. There are some exceptions, one shown in Fig. 22.a. Comparing Fig. 12.a and Fig. 22.a, we see that the resonant contributions and NLO corrections in the decay are similar at the TeVatron as at the LHC in the hard tail of the m b 1 e + -distribution.
• The dependence of the predictions on the three SMC codes, PY1, PY2 and HW is very similar as found at the LHC (for instance, compare the two plots in Fig. 23 to Fig. 17.b and Fig. 18.b).
Thus we focus on making predictions at the hadron level. For this kind of predictions the selection cuts (1-6) were applied after shower, hadronization, hadron decay and the application of jet algorithms. The integrated cross-sections after cuts, in the three cases are collected in Table 5, for different jet sizes (anti-k ⊥ with R = 0.4 versus anti-k ⊥ with      Table 4. Cross-sections from the LHEs, after PS and at the hadron level at the TeVatron after cuts (1-6) (first column) and after an additional jet veto (second column). The quoted uncertainties are statistical only. The PS and SMC predictions are obtained with the PY1 SMC.     Table 5. Cross-sections at the hadron level at the TeVatron after cuts (1-6) for the three cases as a function of the R parameter of the anti-k ⊥ algorithm used to identify jets. The quoted uncertainties are statistical only. The predictions are obtained with the PY1 SMC.    Ratio η e + Figure 25.
Pseudorapidity distributions of a) the hardest b-jet and b) the hardest isolated positron after full SMC for the three cases. The lower inset shows the ratio of the predictions with decays of the t-quarks in DCA and Decayer compared to the complete W + W − bb computation.
Turning to distributions, we show our standard selection in Figs. 24-26, presented similarly to the LHC plots. We see in Figs. 24 and 25 the general trend that Decayer and DCA give very similar predictions both in shape and normalization for p ⊥ -and ηdistributions, while the full W + W − bb-computation followed by the SMC differs, but the difference is smaller than in the case of LHC, usually within 5 %. Thus the features seen in the distributions from the LHEs are kept after full SMC.
The pseudorapidity distributions of the hardest b-jet and of the positron, shown in Fig. 25, exhibit a forward-backward asymmetry computed by all three methods. However, the size of the asymmetry is significantly bigger for W + W − bb-production than that obtained from the DCA or Decayer, although the increase is predicted partially by Decayer.
Looking at the m b 1 e + -distribution in Fig. 26.a, we find that in the hard tail above the kinematic limit at 150 GeV Decayer and DCA approximations give very similar predictions, while that of the full calculation is almost twice as large. Thus the trend observed at the level of the LHE (see Fig. 22.a) survives the SMC, being only slightly attenuated. We present the ∆φ e + µ − -distribution in Fig. 26.b. For this plot the differences are again much smaller than observed for the LHC (cf. with Fig. 15.b).

Conclusions
In this paper we presented the first study of W + W − bb hadroproduction at NLO accuracy matched with SMC. For the matching we used the POWHEG method as implemented in the POWHEG-BOX within the PowHel framework. This framework allows the generation of events including first radiation that can be processed further with the SMC. The events are stored in event files according to the Les Houches accord, and are freely available for experimental analyses. We included all Feynman-graphs in the NLO computation, the doubly-resonant as well as singly-and non-resonant ones. We checked the validity of our computations by comparing cross sections obtained from the LHEs to the NLO predictions available in the literature and found agreement within the expected accuracy. We also studied the effect of the parton shower and hadronization separately. The PS decreases the cross sections by 10-20 % at the LHC (depending on the observable) fairly uniformly. The hadronization results in a further decrease, which however depends on the observable. For p ⊥ -distribution of hadronic objects significant amount of hadronic energy is transformed into electromagnetic and missing energy during hadronization, resulting in large corrections due to a softening shift of the spectrum. This effect is not observed for other observables, such as p ⊥ -distributions of leptons or pseudorapidity distributions. The effect of the hadronization on the inclusive cross section after cuts remains below 5 % both for the LHC and the TeVatron, with somewhat larger effect at the latter machine.
At the end of the SMC the events contain stable hadrons and leptons. Such a final state can also be obtained from LHE files containing events with a tt-pair (including first radiation), and let the SMC decay the heavy quarks. This much simpler computation gives predictions in the decay chain approximation for the original heavy particles of the hard scattering event. With a little more effort we can also perform the decay of the heavy particles including off-shell effects of the t-quarks and spin correlations, implemented in a new program Decayer. An important goal in this paper was to understand the validity of these approximations, with emphasis on the predictions after full SMC.
We found that standard transverse-momentum and pseudorapidity distributions are in general described by the predictions of Decayer and in DCA equally well, although the predictions of the former show the off-shell effects in certain kinematic domains. Both approximations lack the NLO corrections in the decays. The main effect of the latter is to increase the cross section by about 10 % at the LHC and about 5 % at the TeVatron. The obvious remedy to correct for these inaccuracies is to include the NLO corrections in the decays. However, such a project is rather non-trivial in the NLO+PS matching framework and is beyond our present scope.
We observe three types of deviations from the general trend: one that does not require the inclusion of the NLO effects in the decays and two that would not be cured even by decays at NLO accuracy.
The first kind of deviation is related to the different treatment of spin correlations, which is fairly well described by Decayer, but the DCA fails, as expected. The example we discussed is the azimuthal difference between the positron and the muon where the effect can reach up to 20 % at the LHC, while it is smaller at the TeVatron. Any related observable (e.g. ∆R-separation, angular separation between these leptons) shows the same effect.
The second deviation is related to the effect of singly-, and non-resonant contributions in regions of the phase space that are less reachable by the doubly-resonant graphs. Typical example is the m b 1 e + -distribution above m 2 t − m 2 W 153 GeV where there is a sharp fall of the cross section in the fixed order predictions due to kinematic reason.
The third kind of deviation is related to the different probability of emitting a hard jet from a t-quark and from a b-quark (treated massless in our computation). The emission of a parton from a b-quark is dominated by soft and collinear emissions. Such emissions usually become part of the b-jets. As a result the p ⊥ -distribution of the hardest non-b jets is much softer for W + W − bb-production. This effect is very significant both at the LHC and at the TeVatron and can be observed in related distributions that use the transverse momentum of the non-b jets.
Sets of events, both for LHC and for TeVatron, obtained with the same parameters as used for preparing the distributions in Sect. 5 of this paper can be downloaded from http://www.grid.kfki.hu/twiki/bin/view/DbTheory/WwbbProd.

Acknowledgements
events are generated by PowHel and they can be Born-like or after first emission, the final state being a tt-pair or a tt-pair plus one extra parton, respectively. To have a better understanding of this decay approximation it will be illustrated on the process qq → tt g. For this flavor structure the undecayed phase space can be written as dΦ u = dΦ 3 (Q; k t , kt, p g ) , (A.1) while the decayed phase space is dΦ d = dΦ 7 (Q; k b , k e + , k νe , kb, k µ − , kν µ , p g ) .
For a decaying top, two consecutive decays have to be iterated, as the top decays into a b and a W + , followed by the decay of this massive W + into a lepton-antilepton pair (in the leptonic decay channel). Hence the decayed phase space can take the following form The exact formula for the two-particle phase space can be found in Eq. (3.2). To unweight over the decayed phase space we have to calculate F 1 and F 2 of Eqs. (3.10) and (3.15), respectively. While the expression for the latter remains the same as Eq. (3.15), the former is written explicitly as tree-level matrix elements are needed for the undecayed and decayed configurations. The uncorrelated squared matrix element is where M u is the matrix element for tt g production, M t→b W + , Mt →b W − , M W + →¯ ν and M W − → ν are the matrix elements for t, anti-t, W + and W − decays, respectively. The original undecayed phase space is used for computing M u , while the matrix elements concerning the decays use the appropriate momenta from the decayed phase space. The partonic flux factor is included in this matrix element. According to Ref. [15] there exists a U max upper bound, such that where M corr is the matrix element corresponding to bb e + µ − ν eνµ g production, with b e + ν e andb µ −ν µ coming from t-andt-decays, respectively. The upper bound can be calculated analytically, or, like in our case, the decay is started with a pre-defined U max value and, if the event under consideration is proved to have a larger value, the upper bound is replaced by this larger one. All matrix elements needed to perform the decays are taken from HELAC-Phegas.