Top-pair production and decay at NLO matched with parton showers

We present a next-to-leading order (NLO) calculation of $t\bar{t}$ production in hadronic collisions interfaced to shower generators according to the POWHEG method. We start from an NLO result from previous work, obtained in the zero width limit, where radiative corrections to both production and decays are included. The POWHEG interface required an extension of the POWHEG BOX framework, in order to deal with radiation from the decay of resonances. This extension is fully general (i.e. it can be applied in principle to any process considered in the zero width limit), and is here applied for the first time. In order to perform a realistic simulation, we introduce finite width effects using different approximations, that we validated by comparing with published exact NLO results. We have interfaced our POWHEG code to the PYTHIA8 shower Monte Carlo generator. At this stage, we dealt with novel issues related to the treatment of resonances, especially with regard to the initial scale for the shower that needs to be set appropriately. This procedure affects, for example, the fragmentation function of the b quark, that we have studied with particular attention. We believe that the tool presented here improves over previous generators for all aspects that have to do with top decays, and especially for the study of issues related to top mass measurements that involve B hadrons or b jets. The work presented here also constitutes a first step towards a fully consistent matching of NLO calculations involving intermediate resonances decaying into coloured particles, with parton showers.

Top quark pair production is now known at NNLO for total rates [14], and complete differential distributions can be expected to be computed at NNLO in the future too [15,16], at least for on-shell top quarks, and possibly also including factorizable NNLO corrections to the decay [17].
The differential cross section for the on-shell production of tt pairs at NLO has been available for a long time [18], and has been implemented as NLO calculations matched to a parton shower generator (NLO+PS from now on) in the MC@NLO framework [19], and in the POWHEG framework [20]. In both these NLO+PS implementations, the top decay is treated in an approximate way, according to the prescription of ref. [21], and no radiative corrections to the decay processes are included.
In ref. [22][23][24], NLO radiative corrections to production and decay have been computed. The top is treated there as an on-shell particle, with production and decay processes fully factorized. Radiative corrections to W hadronic decays are also included there. Currently, however, the most accurate description of top-quark pair production with fully-exclusive decays is the one obtained at NLO in refs. [25,26], (that will be referred to in the rest of this paper as DDKP), and in refs. [27,28]. These publications contain a computation of NLO corrections to the process pp → e + ν e µ −ν µ bb, hence spin-correlation as well as intermediate top offshellness effects are taken into account exactly. 1 NLO calculations of the complete production process of the top decay products, also including possibly the hadronic decays of the W , can be directly interfaced to parton shower generators like POWHEG, as it was actually done in ref. [29] using the POWHEG BOX framework.
This procedure, however, is bound to fail in the narrow width limit, since (in the current POWHEG BOX) both the POWHEG and the subsequent shower radiation from the resonance decay products may change the resonance virtuality. This problem is particularly evident in cases when a strongly decaying resonance is produced, at the Born level, with nonzero momentum in the partonic CM frame. This is certainly the case for tt production. In order to discuss this problem in a simpler context, we consider the production of a single resonance R, decaying into a coloured parton c plus a massive coloured particle p, in association with a neutral particle n. We also assume colourless incoming particles, like in an e + e − collider, so that we don't need to worry about initial state radiation. The neutral particle n is included in such a way that the resonance has non-zero momentum in the CM frame. The POWHEG formula for this process is where, as usual B(Φ B )dΦ B and R(Φ B , Φ r )dΦ r dΦ B are respectively the Born and real cross section, where the real cross section refers to our basic process with the additional emission of a light parton. V is the virtual contribution. Its soft and collinear divergences cancel against those of the real integral in the expression forB. Again as usual, the real emission phase space is factorized in terms of an underlying Born phase space and radiation variables [30]. The only radiation that we can have now is the one from the decaying resonance. For this example, in the formalism implemented in prior versions of the POWHEG BOX, there is only one singular region, associated with the radiation of a gluon from parton c. In the corresponding underlying Born, the c parton momentum is aligned with the momentum of the c g system. The direction of the momentum of the p n system in the underlying Born configuration is the same as in the real emission configuration, opposite to the c g system (or c, in the underlying Born). The c and p n momenta are rescaled by the same amount, in such a way that the final state energy is the same in the real configuration and in the underlying Born. The rescaling of the p n momentum is achieved with a boost of the p n system along its direction of flight. It is now clear that this procedure does not preserve the virtuality of the resonance. We can easily estimate that the resonance virtuality is changed by an amount of order m 2 /E, where m is the virtuality of the c g system, and E is its energy. Because of this, we expect distortions of the radiation distribution when m 2 /E Γ, where Γ is the resonance width. It is easy to see what can happen by inspecting eq. (1.1). Let us assume that Φ B is on resonance. Then Φ B , Φ r is on resonance only for m 2 /E Γ. The real integral inB is therefore restricted to a vanishing region in the narrow width limit, thus exposing large negative corrections coming from the virtual term. This will yield a large suppression due to theB/B ratio. The Sudakov exponent will yield a substantial suppression only when m is small enough. Let us consider now when Φ B is off resonance. A potentially large contribution inB may be generated when Φ B , Φ r has the resonance on shell, yielding to a largeB/B ratio. This should compensate for the reducedB/B of the previous case. However, in this case the Sudakov exponent will be quite large, yielding a strong suppression. We thus expect considerable distortion of the jet mass spectrum in the region where m 2 /E Γ. We anticipate that there shouldn't be problems if Γ is large enough, but it is not easy to quantify this statement. Future studies comparing the results obtained with the generator introduced in the present work with those of ref [29] will help to clarify this issue.
We remark that problems in the handling of strongly decaying resonances are not specific to POWHEG. In an MC@NLO framework [31], for example, the shower itself can change the resonance virtuality. In order to prevent this one should require that the shower preserves the resonance mass, and MC counterterms would need to be modified accordingly.
In this paper, we present an NLO+PS generator for tt production that includes NLO corrections in production and decay, in the narrow width approximation. In this framework, radiation in production and in decay can be clearly separated, and no ambiguity arises for the assignment of the radiated parton, so that we can avoid the problems mentioned earlier. As we will show, this framework will also help us in understanding and dealing with further problems arising when strongly decaying resonances are involved, and it may eventually lead to a future consistent prescription for a framework that fully accounts for finite width and interference effects. In order to make our generator fully realistic and usable, we include an approximate treatment of finite width and interference effects.
We use the MCFM matrix elements for tt production in the narrow width approximation, introduced in ref. [24] and taken from ref. [32]. The matrix elements are interfaced to the POWHEG BOX V2, that includes a framework for the treatment of decaying resonances, including possible radiation from the decay. The bottom mass is kept finite in our implementation: this allows us to give a realistic description of bottom fragmentation in top decays.
The paper is organized as follows. In section 2 we give a description of our generator. In section 3 we compare our generator at the NLO level with the calculation of DDKP. The purpose of this comparison is to check to what extent our calculation maintains its validity outside of the resonance region, since our implementation of finite width and interference effects is only approximate. Finally, in section 4 we present our phenomenological studies. In section 5 we give our conclusion.

Resonance decay in POWHEG BOX V2
The POWHEG BOX V2 framework includes the possibility of implementing radiative corrections to both production and resonance decays. Nested resonances are correctly handled with no limitations. The implementation is fully automatic as far as the NLO calculation and its interface to a shower is concerned, consistent with the general aim of the POWHEG BOX package. It requires external input as far as the Born phase space and the Born, virtual and real matrix elements. The information needed for the description of the flavour structure of the matrix elements is extended in V2 to include resonance information. In the case at hand, for example, the flavour list for the Born subprocess is supplemented with an array that specifies the resonance assignments of the various final state particles: flst_born(1:nlegborn,j) =[ 0, 0, 6, -6, 24,-24,-11, 12, 13,-14, 5, -5] flst_bornres(1:nlegborn,j)=[ 0, 0, 0, 0, 3,4,5,5,6,6,3,4] where the flst_bornres entry for t andt is zero (they come from the production process), for W + and b it is 3 (they come from the decay of the top, which is the 3 rd entry) for the e + and the ν e it is 5 (they come from the W + , which is the 5 th entry), and so on. Similar structures are assigned for the real emission matrix elements. In this case, also the emitted parton has a resonance assignment, specifying whether it originates from the production process, from a top, or from a W . In the present version of the POWHEG BOX V2, it is always assumed that the resonance assignment at the Born level is the same for all subprocesses. On the other hand, the resonance assignment of the radiated parton in real graphs can vary. When computing a real contribution arising from a given singular region, the program can figure out the origin (i.e. from which resonance the radiation arises), and the real emission matrix element returns the appropriate value on that basis. All other features, like the generation of an appropriate real radiation phase space, the soft and collinear limits and the integral of the soft corrections in case of resonance decays, are handled internally by the POWHEG BOX V2. A full description of these features will appear in a forthcoming publication, whose draft is publicly available in the Docs/V2-paper.pdf file in the POWHEG BOX V2 directory.
The kinematics for the mapping of the real phase space to an underlying Born configuration departs from the one described in the POWHEG BOX [33] when radiation from resonance decays are concerned. In this case, the mapping is similar to the standard POWHEG BOX final state mapping, except that the resonance four-momentum is preserved, rather than the momentum of the whole final state. For consistency, the integrated soft counterterms involving the resonance decay are also computed in the resonance frame. If the resonance is coloured, it is also included as the initial state particle in the computation of the corresponding eikonal factors. The soft corrections, arising from the production process, are instead computed treating the resonances as final state particles.

NLO matrix elements
As already mentioned in the introduction, the NLO matrix elements are taken from MCFM [24]. In this calculation, the real contributions due to radiation in production and radiation in decay are separated, and the POWHEG subroutine providing the real matrix element selects the appropriate process, depending upon the origin of the emitting parton.
As illustrated in ref. [24], in the MCFM program no radiative corrections to the total width of the top in the decay t → W b appear. In fact, such a correction enters the virtual contribution in the top decay, but it should also enter the total top width appearing in the top propagator. If the total top width is set to its tree level value (i.e. no strong corrections are included into it), it will exactly cancel the width for the t → W b process. In other words, the total cross section will correspond exactly to the cross section for the production of a tt pair, since the branching ratio for t → W b is one.
The MCFM calculation also computes the width for the W decay. In case of leptonic decays, this width does not receive strong corrections. The MCFM program thus computes its tree level value. If the total W width that appears in the W propagators is an accurate expression including up to fourth order radiative corrections to hadronic W decays, the corresponding leptonic branching fraction will have the same accuracy.
As far as hadronic W decays are concerned, the MCFM program computes also their O(α S ) radiative corrections, that correspond to the usual 1 + α S /π factor. Since radiative corrections to the W hadronic width are known with much better accuracy, we have preferred to modify the standard MCFM behaviour, by subtracting this α S /π term from the virtual corrections. In this way, we get a total W hadronic width equal to its tree-level value, Γ with a = α S (M W )/π, and set the total W width in the MCFM program to Furthermore, for each hadronic decay of the W , we supply the factor .

(2.3)
It is easy to verify that with this prescription we obtain a total cross section that is equal to the tt cross section, and the correct fraction of leptonic over hadronic events for each W , equal to Γ . (2.4)

Finite width and interference effects
We include finite width and interference effects in an approximate way, as described in the following. We generate full Born phase space for off-shell production and decay of a tt pair, and we set up a procedure for mapping this phase space into an on-shell one, that will be used for the calculation of the on-shell matrix elements. We call p t and pt the top and anti-top off-shell four-momenta, and v t = p 2 t , vt = p 2 t their virtualities. The mapping is such that the off-shell top configuration, with top and anti-top virtualities v t , vt, is mapped into an on-shell configuration with top mass equal to (v t + vt)/2. The mapping of the W fourmomentum invariant is performed as follows. The top decay phase space, together with the W Breit Wigner shape, has weight proportional to (neglecting the b mass) where the dimensionful normalization of the argument of the logarithm amounts to an arbitrary integration constant. We thus introduce the following function, corrected for the and compute the new W virtuality v W in such a way that: where m t = (v t + vt)/2 is the value of the top mass that will be adopted in the computation of the matrix elements. With this choice, the range of the W virtuality remains correctly mapped between zero and the top mass adopted, and the W is pushed off resonance only when absolutely necessary. We now must map the phase space with v t , vt, v W + and v W − virtualities of the decaying resonances into one with virtualities m t , m t = m t , v W + and v W − . This mapped phase space is obtained in a straightforward way by a fully general mapping procedure, that is applied recursively, starting with the last decaying resonances (i.e. the W 's) and ending with the full final state before any resonance decay. In this procedure the system comprising the collision final state before any resonance decay, is also treated as if arising from a fictitious resonance. In the following, by direct decay products of a resonance we mean its decay products before any further resonance decay. The mapping procedure is as follows: • Starting from the last decaying resonance in the decay chain, and following recursively the decay chain in the backward direction, for each resonance V , we go to the V rest frame, and rescale all the momenta of the V decay products by the same factor, in such a way that energy is conserved with the newly assigned V virtuality, and with the newly assigned virtualities of the direct decay products of V (in case some decay products are resonances themselves).
• If any of the direct decay products of V is itself a decaying resonance V , the rescaling of the V momentum is accompanied by an appropriate longitudinal boost of all decay products (direct and indirect) of V along the V direction, performed in the V rest frame, such that momentum conservation holds in the decay.
• The last decaying resonance is the fictitious resonance comprising the direct products of the hard collision. Our procedure is applied also in this case, keeping the energy of the system fixed.
The whole POWHEG generation is carried out with the mapped Born phase space. Thus, matrix elements are computed with an on-shell phase space with equal t andt masses. Furthermore, the POWHEG radiation phase space is computed starting from the mapped Born phase space, and it thus has equal masses for the top and the anti-top. So, while we have at hand the original Born phase space and the corresponding mapped one, we end up with a mapped real phase space, and we need to build an "unmapped" phase space out of that. Furthermore, this must be done taking into account the origin of the radiated parton (that is to say, if it was emitted from a resonance or at the production stage), and preserving the virtuality of all decaying resonances. This is done using the following mapping procedure: • A longitudinal boost B L to the longitudinal rest frame of the tt system is applied to the tt system and to all its decay products, followed by transverse boost B T to the tt rest frame (unless the radiation is from the hard process, we have B T = 1).
• The same procedure illustrated earlier for the Born phase space mapping is applied, where now the original resonance masses are m t , m t = m t , v W + and v W − , and we want to end up with v t , vt, v W + and v W − .
• The boost B −1 L B −1 T is applied to the tt system and to all its decay products.
We have also implemented an alternative mapping method, that can be activated in our program by setting the variable altmap 1 in the powheg.input file. In case of radiation in resonance decays, instead of simply rescaling the momenta of the emitter and emitted parton, we consider the emitter-emitted pair as a whole, and rescale its four momentum. The individual emitter-emitted particle momenta are obtained by a suitable longitudinal boost. This alternative scheme maintains a closer relation between the transverse momentum of radiation in the on-shell and off-shell phase spaces. In practice, this feature is not set as default, since we have not found any noticeable effect due to it. Notice that in this case the radiated parton may belong to any decaying resonance, and one must take care to assign it to the appropriate resonance in order to perform the mapping correctly. This is not a problem in POWHEG, since the mapping between the Born and full phase space is always dependent upon the emitting parton, that thus identifies uniquely the resonance from which the radiated parton is emitted.
The off-shell phase space is used for the NLO level analysis, and for the generation of the Les Houches event.
Besides phase space mapping, we need to supply further dynamical information on the distribution of the top and anti-top virtuality. We have considered two options for this: • Weight all matrix elements with the Breit-Wigner weights of the top and anti-top This behaviour is activated by setting the flag mockoffshell equal to 1 in the powheg.input file.
• Default behaviour: we reweight all matrix elements by the ratio of the full off-shell cross section for the underlying Born process at hand, including interference effects, over the on-shell Born matrix elements. Using this procedure we intend to capture, at least at leading order, all off-shell and interference effects that we have left out in our calculation. In order to avoid useless complications, we have computed the exact matrix elements for all decay processes including identical and different fermions in fully leptonic decays, semileptonic decays into first or second generation quarks, and fully hadronic decays into two different quark generations (i.e. we don't include interference effects for decays into the same generation quarks). For any required flavour combination, the program maps it to the the closest one that matches one of the above. The matrix elements were obtained with the MadGraph5 program [34].

Les Houches event generation and interfacing to a shower
Our tt generator can be run with different options controlling the radiation from resonance decays. The simplest case is to exclude radiation from resonance decay, which is achieved in our program by setting the variable nlowhich to 1 in the powheg.input file. In this case, our generator is an improvement over the previous implementations of refs. [19,20] only because it includes spin correlation in production with the full NLO accuracy. The events generated in this way can be passed to any Shower Generator like PYTHIA and HERWIG that is compliant with the standard Les Houches interface for user processes [35]. Shower Monte Carlo's treat radiation from resonance decays using their own shower formalism, and, as far as further radiation from the production process is concerned, radiation hardness is limited by the value of scalup, that is passed to the program via the Les Houches interface. This guarantees that the shower will not generate any radiation harder than the one generated by POWHEG. No limitations are instead imposed on the radiation of decaying resonances. By setting nlowhich to 0, radiation from resonances is also generated. In this case, POWHEG uses the usual highest bid mechanism to decide whether a given event has initial state radiation or radiation from one of the resonances. This implies that the subsequent shower program should also limit the radiation in resonance decays in such a way that no radiation whatsoever, both from the production process and from resonance decays, is generated with a transverse momentum larger than that of the POWHEG generated radiation. The Les Houches Interface does not provide for a standard mechanism to veto radiation in resonance decays. We will explain further on the procedure that we adopted to achieve this goal. A further problem arises in the present case. When using the above procedure, for each event, at most one of the decaying resonances will include an NLO accurate radiation. Furthermore, since initial state radiation turns out to be much more frequent (due to the large center of mass energy available at the production stage), corrections to decays generated by POWHEG itself will not take place so often, with the consequence that the Monte Carlo will be mostly responsible for radiation in resonance decays. On the other hand, a better description of top decay is quite desirable, especially for studies aimed at the determination of the top mass.
In order to deal with the above mentioned problem, a further option is included, that is activated when the variable allrad is set to 1. If this is the case, the POWHEG generator departs from the standard behaviour at the stage of the generation of radiation. Radiation from each allowed singular configuration is generated, that is to say initial state radiation, radiation in the decay of the top and anti-top, and if hadronic decays are required, also radiation from the corresponding decaying W . Normally, only the hardest radiation is kept. In the modified behaviour, all radiations are instead implemented, so that one can find on the Les Houches event structure as many radiated partons as there are resonances decaying into coloured particles, plus possibly one radiation from the initial state. The real radiation phase space configurations are all preserved, and at the end they are merged into a single phase space. The merging can be easily performed by using only two kinds of operations: a boost of the W system in the t (ort) rest frame, in order to account for radiation in top decays, and a longitudinal boost composed with a transverse boost applied to the whole tt system, including its decay products, in order to account for initial state radiation.

Comparison with existing results
Currently, the most accurate description of top-quark pair production with fully-exclusive decays is the one obtained at NLO in refs. [26,27]. These two publications contain a computation of NLO corrections to the process pp → e + ν e µ −ν µ bb, hence spin-correlation effects as well as intermediate top and W offshellness effects are taken into account exactly. Since no narrow width approximation is used, they contain at the amplitude level diagrams which are single-resonant, or non-resonant. In this respect, they are an improvement upon the results of [22,23], which were based on the narrow-width approximation, and an improvement upon the MCFM results that we use here.
The computations in refs. [26,27] have been performed in the 5-flavour scheme, i.e. b quarks are considered massless. More recently, two groups have completed the same computation in the 4-flavour scheme [36,37]. This is phenomenologically important since it allows for the first time a unified description of the tt, W t and "b-quark associated + − νν production" processes (since the b quarks are massive, the result is finite without explicitly requiring two b-jets in the final state, which is obviously needed in a 5-flavour computation).
Since our procedure includes at LO the exact full matrix elements for pp → e + ν e µ −ν µ bb (with massive b quarks), and we also include an approximate treatment of offshellness effects at NLO via the rescaling procedure explained in sec. 2.3, it is interesting to check how well our results compare to an exact NLO computation. To this end, we have performed a comparison (as tuned as possible) with the results presented in DDKP.
The results presented in this section have been obtained for the 8 TeV LHC, using the following input parameters: Since the calculation of DDKP was performed with zero mass for the b quark, the results shown in this section have been obtained setting m b = 0 in our calculation. 2 Moreover, in this section only, we removed all subprocesses with b-quarks in the initial state, since these channels were not included in the DDKP paper. The electroweak coupling α is derived from the Fermi constant in the G µ -scheme: In order to match as closely as possible the procedure adopted in DDKP, we also use a complex-valued expression for the weak mixing angle when evaluating the off-shell MadGraph5 matrix elements, i.e. we define On-shell (MCFM) matrix elements are instead computed with real-valued couplings. The use of a complex-valued weak mixing angle leads to sub-permille differences, which are appreciable at the level of total cross-sections, and it is needed to bring our LO results in agreement to those quoted in DDKP at that level of precision. We use the LO value Γ LO t = 1.4426 GeV at LO. At NLO we use Γ NLO t = 1.3167 GeV in the MadGraph5 matrix elements used for reweighting, whilst Γ LO t is used to evaluate all the core (MCFM) NLO amplitudes. The reweighting procedure described at the end of sec. 2.3 is corrected by an overall factor (Γ NLO t /Γ LO t ) 2 , which is needed to keep the computation consistent with the LO result in the narrow-width limit. Finite width effects for the W boson are fully included in all components of our generator (including the computation of the top width) by default, thus corresponding to the scheme FwW, in the notation of DDKP.
In DDKP two choices for the renormalization and factorization scale are discussed for LHC energies: where E T,i is the transverse energy defined as E 2 T,i = m 2 i + p 2 T,i . As in DDKP, when using the dynamic scale defined in eq. (3.4), we compute the factorization and renormalization scale strictly as a function of the event kinematics. 3 This also guarantees the subtractionscheme independence of the NLO result.
We use the MSTW2008 parton distribution set [40] at LO/NLO order, and for the results presented in this section, we evaluate the strong coupling constant using the corresponding LHAPDF subroutine. In the following we will stick to this set. Other pdf sets, like those of refs. [41,42], can be used as well, but we are not interested in a pdf comparison in the present study.
We cluster final state b quarks and gluons with |η| < 5 into jets using the anti-k T algorithm [43], with R = 0.5, as implemented in the FastJet package [44]. The cross section and distributions shown in the following are defined after these cuts: where we require two b-jets in the final state. At LO, we obtain perfect agreement with DDKP, as expected: Plots in fig. 1 display a comparison between our LO results and DDKP for a selection of observables, obtained using the dynamic scale µ dyn . The observables involving a top or an anti-top shown in this section are "reconstructed", i.e. they are computed by combining the momenta of the W and b-jet associated to the t (t) quark, after the cuts in eq. 3.5. All plots are in good agreement.
At NLO, the procedure we adopted is quite different from the exact approach used in DDKP. As such, we don't expect necessarily a perfect agreement for total cross sections,   In fig. 2 we show the same set of distributions in the NLO case, obtained again with dynamic scales. A noticeable difference between our NLO result and the one of DDKP is seen in the shape of the reconstructed top invariant mass, with our result being about 10% higher very near the peak. This difference does not appreciably alter the peak position and width, and is thus unlikely to cause visible errors in, for example, the determination of the top mass, as can be seen in the bottom-right plot of the figure. Small differences are also observed in the tt mass in the region very close to threshold. The lower cross section there is consistent with our lower cross section for top virtualities below the top mass. We also see small differences in the transverse momentum distribution of the top.
We will show in the following that our description of the peak region turns out to be independent of the reweighting procedure that we adopt, i.e. reweighting using a simple Breit Wigner weight, or including only double resonant graphs in the off-shell matrix elements, would yield the same discrepancy with the DDKP result. This, together with the fact that we reproduce exactly the DDKP Born result, indicates that the differences that we observe at NLO in the peak region must be due to the radiative corrections, that  in the DDKP case also include interference effects. Since the code in DDKP is not public, we cannot investigate this issue further.

Modelling of offshellness and interference effects at fixed-order
We now show a comparison, at the NLO level, of our different approximations to the finite-width and interference effects. The purpose of this subsection is to show the level of consistency of the various approximations we use. We consider: • the results obtained by reweighting the cross sections using a Breit-Wigner function. This is obtained by activating the flag mockoffshell as explained earlier. The corresponding curves are labelled "BW"; • the results obtained by reweighting using the Double Resonant approximation to the LO off-shell matrix elements, labelled as "DR"; • the results obtained by reweighting with the full LO off-shell matrix elements. The corresponding curves are labelled as "full".
In this section we use a finite value for the b-quark mass (m b = 4.75 GeV), which gives Γ LO t = 1.4386 GeV and Γ NLO t = 1.3145 GeV as values for the LO and NLO top width, respectively. Observables labelled with t andt refer to quantities at the "MC truth" level. Notice that in our calculation we can always define an "MC truth" level for the top quarks, that corresponds to the top quarks in the on-shell calculation, mapped with our procedure into the off-shell kinematics.
Before we begin our comparison, we would like to point out that there are phase space regions with a relatively large cross section due to non-resonant contributions: the regions with either of the b quarks collinear to the beam axis, and the region where the b and b quarks are collinear. These regions are kinematically enhanced, due to the processes represented in fig. 3. They can be excluded from the analysis by imposing a lower limit on the b andb transverse momenta (relative to the beam axis), and on the mass of the bb system. Other regions with potentially important contributions from non-resonant graphs also arise if we consider same flavour decays of the W s, and the corresponding singular regions should be limited accordingly. Here we limit ourselves to consider µ + e − decays, thus avoiding further problems. We remark that in our calculation all these non-resonant contributions can only arise from the "full" results.
In fig. 4   with collinear b's, i.e. we require These cuts are applied at the parton level, without any jet definition. As can be seen, after these cuts are applied, differences among the various approximations are much smaller.
Only few percent differences are visible in the tail of the transverse momentum spectrum of the tt system. The effect of the cuts in eq. 3.8 is better displayed in fig. 6, where we compare the invariant mass of the top-quark and tt system (both at the MC truth level). Here we notice that without imposing the cuts on b-quarks, the results from our different approximations are quite different: below the peak, the m t distribution in the "full" case falls off much less rapidly, due to the effect of collinearly-enhanced non-resonant contributions, which affect, to a minor extent, also the tail above the peak. The "DR" and "BW" results are instead quite similar, especially close to the peak, as expected. After including the cuts in eq. 3.8, the three predictions are much closer to each other. Far below the peak, non-resonant contributions are again noticeable, but the "full" distribution now falls off more rapidly, since collinear-enhanced terms are suppressed by the cuts. The results for the invariant mass of the tt system show that there are ambiguities in the final results when approximations like "DR" and "BW" are used, especially below threshold. After cuts, all the three curves are closer, except below threshold, where the "full" result can be considered our best prediction, although it misses genuine NLO effects in this region.
Since in this work we are mainly interested in studying phase space regions close to the peak, in the rest of this section we will show results after the cuts of eq. 3.8.
In fig. 7 we display our results for the top-quark invariant mass, as defined from the MC truth, or as reconstructed from the W + and the b jet. In both cases we see a fair  agreement between the BW and DR results around the peak area, that is worsening as we go above or below the peak (and more so for the m t plot). We notice that the reweighting applied in DR can compensate unwanted features of the mapping procedure. This is not the case for the BW result, where no reweighting is applied. Thus the agreement of the BW result with the DR one near the peak region means that our mapping procedure is quite sound.
The full result displays a larger distance from the DR and the BW result, especially below the resonance region. This is understood as being due to other, non-resonant production mechanisms entering into play, and being favoured by the growth of the parton luminosities for smaller invariant mass of the produced system. Notice that in the m W + j b distribution (right panel of fig. 7), the difference between full and DR/BW schemes is reduced with respect to the corresponding result for "MC truth" tops ( fig. 7, left panel). This is again expected, because QCD final state radiation outside the b jet cone leads to a smaller invariant mass of the reconstructed top. This effect slightly depletes the peak, and raises the region below it, as can be clearly seen in the figure, thus hiding the relative importance of non-resonant production mechanisms.
We notice that the distribution plotted on the right in fig. 7 is fairly similar to the top left plot in fig. 2, where we compare with the DDKP results. As anticipated earlier, we interpret the stability of the distribution for our three levels of approximation as an indication that the differences must have to do with the full radiative corrections (including also interference between production and decay) with respect to the on-shell ones that we apply.
In fig. 8 we display the invariant mass distribution of the lepton-b-quark and lepton-

POWHEG results
In this section we discuss our POWHEG results, both at the level of Les Houches (LH) events and after full parton shower (PS). The observables that we use for our plots are defined as follows: • At the LH level b, t and W stand for the corresponding LH level particles. b(b)jets are obtained clustering all coloured final state particles and then selecting those containing a b(b) quark.
• at NLO+PS level, b is the b-quark after shower but before hadronization (we look for the last b-quark in the shower that has the top as ancestor). We denote as B(B) the hardest (i.e. largest p T ) b(b) flavoured hadron, irrespective of whether it comes from a top or not. We denote as j B (jB) the jet containing the B(B). The t and W are identified at the MC truth level (i.e. the last t or W in the shower record).
We use the same inputs as those used in the last part of the previous section, although from here on, unless specified otherwise, we will impose on all events the same cuts adopted in DDKP (see eq. 3.5), i.e.
in the NLO+PS case, and same cuts with B replaced by b in the LH case. Jets are reconstructed, as before, with the anti-k T algorithm with R = 0.5. The renormalization and factorization scales are set equal to the transverse mass (i.e. m 2 t + p 2 T,t ) of the t or t at the underlying Born level.
In the following, we will first describe the differences between results obtained with or without strictly zero widths for the top and the anti-top. This is to illustrate the need of an off-shell treatment when analyzing observables useful for the determination of the top mass. Then we will compare the different approaches that we have considered for modeling offshellness effects.

Comparison of on-shell vs off-shell results
We first study the differences between our full result, and the on-shell result, obtained in the strict zero width limit. 4 In this section we include radiative corrections both in production and decay. We do however stick to the usual "additive" procedure, that is to say we do not activate the allrad flag discussed in section 2.4. We discuss results at the LH level, and after showering with PYTHIA8.
We begin by showing the lepton and top transverse momentum distributions in fig. 9. Since the cuts in eq. 4.1 remove phase-space regions where non-resonant diagrams can become collinearly enhanced, we expect the predictions for these quantities to be in very good agreement. This is what we observe in the plots, that are shown at the LH level (left column) and after full parton-showering/hadronization but without MPI effects with PYTHIA8 (right column). As expected, one observes a K-factor very close to 1 and constant between the two predictions.
In fig. 10 we show the invariant mass of the lepton-b and lepton-b-jet systems. As expected, width effects are particularly relevant near the end-point, also for the case of the lepton-b-jet mass. The effect is larger at the Les Houches level, where only one radiation can be combined with the b quark to form a jet. After shower, the larger activity of the event reduces the effect, that remains however visible. This is bound to have some impact upon top mass determinations using end-point observables. Notice that, as expected, in the on-shell calculation at the Les Houches level there is a sharp threshold in the bl + invariant mass. Such sharp threshold is not present in the l + j b case, because initial state radiation can combine with the b quark to yield a more energetic jet. Also in the l + B case we observe a tail above threshold, since colour conservation requires that the b quark should hadronize in combination with some other quark not coming from top decay.

Impact of different options on off-shell results
In this section we show a comparison at NLO+PS level between our results obtained using the traditional "additive" POWHEG procedure (i.e. nlowhich 0, denoted as "nlow0" in the following plots), against the results obtained using the new procedure described at the end of sec. 2.4 (i.e. nlowhich 0 and allrad 1, denoted as "allrad" in the figures). For comparison, we also show a curve obtained by setting nlowhich 1 (denoted as "nlow1"), that corresponds to let POWHEG generate only initial state radiation, leaving entirely to the PS program (PYTHIA8 in our case) the task of emitting gluons from the b-quark originating from the top decay. We will also compare results with the setting nlowhich 4 ("nlow4" in the figures), that corresponds to let POWHEG generate only radiation from top decays (i.e., no POWHEG radiation from either production ort decays). We begin by showing in fig. 11 two observables that characterize the kinematics of the [GeV] allrad nlow0 nlow1 nlow4 Figure 11: Properties of the B hadron decay in the top rest frame: the B fragmentation function (left), and the B transverse momentum with respect to the W direction (right), plotted for the four running modes described in the text, and using the cuts in eq. 4.1. Notice that in the "nlow1" case, radiation in decay is fully handled by PYTHIA8 alone.
B hadron coming from the top decay in the top rest frame, looking only at the positively charged top: the B fragmentation function, and the transverse momentum of the B hadron with respect to the direction of the W , for the four running modes described above. We immediately see that the "allrad" and the "nlow4" results are strictly proportional. In fact, in both "allrad" and "nlow4" the hardest radiation in the top decay is generated by POWHEG, while the subsequent ones are generated by the parton shower. On the other hand, since the total cross section in the "nlow4" case is equal to the Born cross section, the two curves differ by the NLO K-factor. Thus, the "allrad" result compares as expected with the "nlow4", and we now consider the remaining two cases with respect to the "allrad" one.
In the "nlow1" case, the radiation in decay is handled by the parton shower alone. We see that for soft radiation, corresponding to large values of x and small p W −axis T,B , the parton shower shape differs appreciably from the "allrad" one, and that for hard radiation (corresponding to small values of x and large p W −axis

T,B
) the shower yields a harder result. In the "nlow0" result, POWHEG generates radiation from the production stage and from t andt decay, and it picks the hardest among them. Initial state radiation is more likely to prevail, because of the wider kinematic range available, and also because of the stronger colour coupling of the initial gluons compared to the final b quark. Radiation in top decays prevails only if relatively hard. This clarifies the behaviour of the "nlow0" result in the figure, that in the hard region approaches the "allrad" one, while in the soft region agrees with the "nlow1", i.e. with the parton shower result. We also notice that the difference between POWHEG and PYTHIA8 in the soft region, where most of the cross section is, forces a normalization difference for the distributions in the hard region, since the total normalization of the curves should be (at least in the absence of cuts) the same.
At this point it is clear that by choosing one of our two most accurate schemes, i.e. the "allrad" or the "nlow0", we decide to rely upon the description of the soft region that is given either by POWHEG+PYTHIA8, or (essentially) by PYTHIA8 alone respectively. Since in the soft region the two frameworks have formally the same accuracy, more studies are needed to clarify the origin of the observed difference, that can be due to the used value of Λ QCD , to subleading terms in leading log resummation, or to other reasons. We also notice that because of the interplay of the soft and non-perturbative regions, the parton shower tune is also bound to have a non-negligible impact. We will not perform such study at this stage. The purpose of the present paper is to provide a tool by which these questions can be addressed.
In figure 12 we compare the mass of W + B and W + j B systems with "nlow1", "nlow0" and "allrad". Both these observables should be peaked near the top mass. However, in the W + B case, some momentum is lost in the radiation of the b quark, and we expect that some events from the top peak will leak to smaller mass values. To a much lesser extent, this should also happen in the W + j B case, since some radiation is lost out of the jet cone. We see that the differences between the results displayed in fig. 12, can be explained in terms of this effect, in the light of the results of fig. 11. In fact, the "nlow1" result is higher than the "allrad one" below the top peak for both observables, since it has a softer fragmentation function (mostly affecting the W B observable), and a harder nlow0/allrad nlow1/allrad Figure 12: Mass of the W -B-hadron and W -B-jet systems obtained with "nlow1", "nlow0" and "allrad", using the cuts in eq. 4.1.
spectrum (yielding more energy loss outside the jet cone, and thus affecting the W + j B observable). The "nlow0" case is similar to "nlow1", except that for regions affected by very hard radiation, it tends to be closer to the "allrad" result.
The same effect can be used to explain the differences observed in the end-point observables m l + B and m l + j B , displayed in fig. 13. nlow0/allrad nlow1/allrad Figure 13: Mass of the l + B and l + j B systems obtained with "nlow1", "nlow0" and "allrad", using the cuts in eq. 4.1.
In fig. 14 we display some observables that are not particularly sensitive to the shape of the top peak. We see that in this case the three approaches yield quite similar results. In the transverse momentum of the tt pair we observe a difference between the "nlow0" and the other two results. This difference is of little significance, being about 2.5% in the second bin, and 5% in the first bin, where however the cross section is rather small. We have reasons to attribute this difference to an imperfect match between the meaning of our hardness scales and the ones adopted by PYTHIA8. We will discuss further this issue in appendix A.

Comparison with the previous POWHEG generator
A POWHEG generator for heavy quark production, capable to describe tt production and decays, but including neither exact LO spin correlations nor radiative corrections to decays, was presented long ago [20], and is now implemented in the POWHEG BOX framework.
Here we compare our new generator with the old one, that we will refer to as the hvq generator. We will limit the comparison to the "nlow1" mode of the new generator, since it is the closest to the hvq one. The differences between the "nlow1" and the other modes have been studied in the previous sections, so no other comparisons are needed. The new "nlow1" generator differs from the hvq one only in the fact that the latter implements spin correlations in an approximate way, according to the method presented in [21]. We will essentially address two questions: • Whether there are observable differences between the new generator and the old generator for quantities that do not depend upon the direction of the top or anti-top decay products (i.e. that are not sensitive to spin correlations).
• Whether there are differences for observables that are sensitive to spin correlations.
As far as the first question is concerned, we have found, among the large set of observables that we have examined, small differences in the transverse momentum distribution of the tt pair and in the bottom fragmentation function, as can be seen in fig. 15 Figure 15: Comparison of our "nlow1" generator and the hvq generator for the tt pair transverse momentum (left) and for the B meson fragmentation function (right). The hvq result is multiplied by a factor of 1.05, to match the normalization of the "nlo" one. No cuts have been imposed to obtain these plots.
cross section in the "nlow1" and hvq case differ, the "nlow1" one being larger by about 5%. This can be attributed to the fact that other production mechanisms can act in the "nlow1" case, yielding a larger cross section below the 2m top threshold for the invariant mass of the tt pair. The transverse momentum difference is visible below 5%, and it affects the transverse momentum of the pair up to 4 GeV. It must be attributed to the different way in which radiation is produced in the two generators. It is clearly a minor effect, not likely to affect significantly realistic analysis. The difference in the fragmentation function appears as a light distortion in shape, for large values of the B energy fraction. It is entirely due to PYTHIA8, since in both the hvq and "nlow1" generators radiation in decay is handled by the shower generator. We tentatively attribute this difference to the fact that the slightly different radiation pattern in production can influence strong radiation in the decay of coloured resonances. We should in fact remember that the bottom quark is colour connected (via the top quark) to particles arising from radiation in production. Notice that these differences in the fragmentation function are much smaller than those found between "allrad" and "nlow1", due to the difference between an NLO exact and a shower approximate treatment of radiation in decay.
In case severe cuts are imposed on the mass of the top decay products, even larger differences are observed between the old and new generator. In fig. 16 we plot the transverse momentum of the tt pair when such cuts are imposed. The relatively important difference is easily traced back to the way in which off-shell effects are simulated in the hvq generator. The event generation starts there with an on-shell tt pair, produced according to the correct on-shell probability. The momenta of the final state particles are then modified with a probability proportional to the tree-level, off-shell matrix elements. This is done, however, without changing the overall initial probability. The mapping of the on-shell to the offshell kinematics is obtained by first assigning new virtualities to the top and anti-top, and then performing a reshuffling procedure, that does not change the pair 3-momentum in the partonic CM frame. Lower top and anti-top virtualities do however imply lower total CM energy, and are therefore favoured, since they imply higher luminosity. A tight cut around the top mass window will thus cut away these events, especially for low invariant mass of the tt pair, and for low value of the transverse momentum, where a reduced top virtuality has more impact on the luminosities. Although it is unlikely that such features may actually spoil realistic analysis, it is also true that they expose limitations of the hvq generator that are not present in the new implementation.
In order to study spin correlation effects, we have examined the azimuthal distance of the two leptons, and the observables cos θ l + , cos θ l − . These are defined as the cosine of the angle between the lepton in the top rest frame, with respect to the flight direction of the top in the centre-of-mass frame of the tt system (see ref. [45]). We have plotted cos θ l + in ten bins of cos θ l − , and we have computed the average of cos θ l + cos θ l − . We have seen no significant differences among the "nlow1" and hvq generators for any of these observables, also when imposing the generic cuts that we have used in the present work. A more detail comparison of the two generators for observables sensitive to spin correlations is left to future work.

Conclusions
In this work we have presented an NLO accurate generator for tt production, implemented in the POWHEG BOX framework, where radiative corrections in decays are also included. This generator is based upon an NLO calculation performed in the narrow width approximation. Finite width effects are however included in an approximate way. This generator improves over a previous one [20] by the inclusion of exact spin correlations in decays, and by the inclusion of NLO corrections to the decay processes.
We have compared the NLO-level output of our generator with the calculation of ref. [26], where all graphs leading to the final state for fully leptonic top decays are included with all their interferences. With respect to this calculation, we reproduce exactly the leading order result, and reproduce fairly the NLO one, provided one does not look too far away from the resonance regions.
We have introduced three different methods for the approximate implementation of off-shell effects. The least accurate one uses a phase space mapping of off-shell configurations into on-shell ones, and mimics the corresponding difference in the matrix element by a simple Breit-Wigner reweighting (here dubbed BW method). The next method (the DR method), reweights the matrix element using the ratio of the exact LO double resonant matrix element divided by the on-shell one, both computed with the kinematics of the underlying Born. Finally our most accurate method is similar to the DR one, but makes use of the full LO cross section for the given final state, including also non-resonant contributions and interference effects. We have found that the DR and BW result are in fairly good agreement, thus giving us confidence that our off-shell-on-shell mapping procedure is sound. The full result is also in good agreement with the DR one, provided that cuts suppressing the non-resonant bb production via gluon splitting are applied.
We have studied in detail the effect of including NLO corrections in the top decay, compared to letting the shower Monte Carlo handle the decay process. By doing this, we have learned that it is better to use an approach where production and decay are factorized, and we have implemented such an approach in our generator. More specifically, since POWHEG normally generates only the hardest radiation, in a standard POWHEG implementation, only the hardest part of radiation in decay would be corrected at the NLO level, while with our factorized prescription, the hardest radiation in decays is always generated by POWHEG.
We found that NLO corrections in top decays mostly affect the B meson fragmentation function, and to a smaller extent observables that measure the hardness of radiation in the top rest frame. We are aware, however, of the fact the B fragmentation properties also depend upon Monte Carlo tuning. Ideally, one should fit the parameters that affect B fragmentation with a POWHEG+Shower setup for e + e − → bb, and then use the same fragmentation parameters in the tt case.
We believe that our generator will be particularly useful to discuss the systematics of top mass determination from end point observables, like the lepton-B-meson and lepton-Bjet invariant mass, and on similar observables that make also use of the measured missing transverse energy.
We have also compared our new generator with the previous POWHEG one [20], in order to see whether the better treatment of off-shell effects and the inclusion of exact spin correlations, leads to other important differences. As far as spin correlations are concerned, we have not seen relevant differences. On the other hand, we have noticed that the previous implementation of off-shell effects may lead to minor differences in the transverse momentum distribution of the tt pair in the small transverse momentum region, where we have reasons to believe that the new implementation is more solid.

A. Interface to PYTHIA8
Several issues have arisen with regard to the interface between POWHEG and PYTHIA8. In the following, we first discuss our relevant choices for the settings, and then the entire procedure that we have used to implement our radiation vetos in the various cases.

A.1 PYTHIA8 and scalup settings
Generally, PYTHIA8 should comply with the requirement that radiation with hardness above the scalup variable [35] should be vetoed. There is, however, an exception that does apply to our case. According to the PYTHIA8 manual, 5 if no light partons (not arising from resonance decays) are present in the final state, the scalup variable is ignored. The PYTHIA8 rational for this behaviour is that, since no radiation is present at the LH level, there is no danger of overcounting by allowing radiation in the full kinematically allowed region. In our case, and especially in the "nlow0" case, radiation in resonance decays and production do compete, and it is therefore necessary to veto ISR even if no other radiation in production is present at the LH level. This can be achieved in PYTHIA8 with the calls pythia.readString("SpaceShower:pTmaxMatch = 1") pythia.readString("TimeShower:pTmaxMatch = 1"), that force PYTHIA8 to blindly honour the scalup requests. All our PYTHIA8 runs have been performed with the above settings.
A.2 Implementation of the radiation vetos in resonance decays.
Although PYTHIA8 does implement a method for vetoing radiation in resonance decays, we have preferred a different approach that allows us to have better control. Our hardness definition in case of radiation from b quarks in t decays, is given by in the top rest frame, that corresponds to the transverse momentum of the gluon E 2 g θ 2 bg in the soft-collinear limit. We found evidence that this definition does not match exactly the one used by PYTHIA8. We thus do the following. We thus look at the last top in the PYTHIA8 event record. This corresponds to the top after any (production stage) radiation, and before its decay. PYTHIA8 provides pointers to the first and last particle in the shower record coming from the top decay. The list of decay products can consist of the following: • Two particles: a W and a b. This happens if only a W and a b where present as decay products of the top in the Les Houches interface. In this case we follow the b sons in the shower history. We have two cases: -The b is pointing to a list containing a b and a gluon. In this case the scale of the shower radiation in the top decay, that we call St, is computed applying formula A.1 to this last bg pair.
-The b is pointing to a list of a single b quark. In this case St is set to a negative number.
• Three particles: W , b and g. This happens if a W , a b and a gluon were present as decay products of the top in the Les Houches interface (i.e. POWHEG has generated the hardest radiation). The following cases may happen: -The b is pointing to a list containing a b and a gluon. In this case the scale of the shower radiation from the b in the top decay, Stb is computed applying formula A.1 to this last bg pair.
-The b is pointing to a list of a single b quark. In this case Stb is set to a negative number.
As far as the gluon is concerned, we also have two cases: -The gluon is pointing to a list containing two light partons 1 and 2. In this case the scale of the shower radiation from the gluon in the top decay Stg is computed applying formula: again corresponding to the square of the transverse momentum of the softest particle in the soft-collinear limit.
-The gluon is pointing to a list containing a single gluon. In this case Stg is set to a negative value.
We apply this procedure to both the top and the anti-top, obtaining two scales, that we call St+ and St-. We use them as follows, in our different running modes: • "allrad": In this case we also compute the hardness of radiation in top and anti-top decay at the Les Houches level, Pst+ and Pst-respectively. We veto the event if either St+ > Pst+ or St-> Pst-.
• "nlow1": we never veto. In this case, the shower builds the resonance radiation according to its own algorithms, with no restrictions. Radiation in production is instead vetoed internally by the shower itself, according to the value of scalup.
• "nlow4": In this case, we veto the event if St+ > scalup. The radiation from the anti-top is handled by the shower. We take care to store the initial scalup value, to be used in vetoing radiation in t decays, and to set it to the mass of the tt pair, in such a way that radiation in production is handled by the shower alone.
When we veto an event, we discard it and re-run the shower with the same Les Houches event. In other words, we keep showering the same Les Houches event until all our conditions are met.