An NLO+PS generator for $t\bar{t}$ and $Wt$ production and decay including non-resonant and interference effects

We present a Monte Carlo generator that implements significant theoretical improvements in the simulation of top-quark pair production and decay at the LHC. Spin correlations and off-shell effects in top-decay chains are described in terms of exact matrix elements for $p p \to \ell^+\nu_{\scriptscriptstyle\ell}\, l^-\bar{\nu}_{\scriptscriptstyle l} b \bar{b}$ at NLO QCD, where the leptons $\ell$ and $l$ belong to different families, and $b$ quarks are massive. Thus, the contributions from $t\bar{t}$ and $Wt$ single-top production as well as their quantum interference are fully included. Matrix elements are matched to the Pythia8 parton shower using a recently proposed method that allows for a consistent treatment of resonances in the POWHEG framework. These theoretical improvements are especially important for the interpretation of precision measurements of the top-quark mass, for single-top analyses in the $Wt$ channel, and for $t\bar{t}$ and $Wt$ backgrounds in the presence of jet vetoes or cuts that enhance off-shell effects. The new generator is based on a process-independent interface of the OpenLoops amplitude generator with the POWHEGBOX framework.

top-pair production and decay. The vast majority of such calculations rely on the narrowwidth approximation (NWA), where matrix elements for tt production and decay factorise. Various generators based on the NWA approximation [18,20,[23][24][25][26][27][28][29] apply NLO QCD corrections only to tt production and include finite-width effects and spin correlations in an approximate way using the method of Ref. [31]. The best available NWA fixed-order calculations implement NLO QCD corrections to the production and decay parts with exact spin correlations [32][33][34]. The ttb NLO dec 2 generator of Ref. [35] implements the results of Ref. [34] using the POWHEG method [21,22]. Finite width and interference effects are implemented in an approximate way, using LO pp → W + W − bb matrix elements. Thus, in the resonance region it provides NLO corrections to both production and decay, including NLO corrections to W hadronic decays, and implements full spin correlations. In addition, it can be operated both in the five-flavour number scheme (5FNS) and in the four-flavour number scheme (4FNS).
A complete description of tt production and decay beyond the NWA requires the calculation of the full set of Feynman diagrams that contribute to the production of W + W − bb final states, including also leptonic or hadronic W -boson decays. The existing predictions at NLO QCD [36][37][38][39][40][41] deal with the different-flavour dilepton channel, pp → + ν l −ν l bb. Besides an exact NLO treatment of spin correlations and off-shell effects associated with the top-quark and W -boson resonances, such calculations account for non-factorisable NLO effects [42][43][44] and provide an exact NLO description of the top resonance, including quantum corrections to the top propagator. Moreover, in addition to doubly-resonant topologies of tt type, also genuine non-resonant effects stemming from topologies with less than two top or W -propagators are included, as well as quantum interferences between different topologies.
The first NLO calculations of the pp → + ν l −ν l bb process [36][37][38][39] have been performed in the 5FNS, where b quarks are treated as massless particles. In the meanwhile, NLO QCD predictions in the 5FNS are available also for + ν l −ν l bb production in association with one extra jet [45]. Due to the presence of collinear g → bb singularities, the applicability of these calculations in the 5FNS is limited to observables that involve at least two hard b jets. This restriction can be circumvented through NLO calculations 3 in the 4FNS, where b quarks are treated as massive partons [40,41]. In addition to a more reliable description of b-quark kinematics, these calculations give access to the full + ν l −ν l bb phase space, including regions where one or both b quarks become unresolved. This is crucial in order to describe top backgrounds in presence of jet vetoes. Moreover, inclusive + ν l −ν l bb calculations in the 4FNS guarantee a consistent theoretical treatment of single-top W t production at NLO.
In the 5FNS, W t and tt production and decay involve partonic channels of type gb → W + W − b and gg → W + W − bb, respectively. The gg → W + W − bb channel at LO is part of the NLO radiative corrections to the gb → W + W − b one, thus yielding a NLO correction that, being tt mediated, is much larger than the Born term. This led to the proposal of various methods [47][48][49][50] to define single top cross sections not including the resonant tt contribution. However, the separation of tW and tt production breaks gauge invariance and does not allow for a consistent treatment of interference effects. On the other hand, in the 4FNS the pp → + ν l −ν l bb calculations provide a unified NLO description of tt and W t production, with a fully consistent treatment of their quantum interference [41]. Single-top production in the 4FNS is described by topologies with a single top propagator and a collinear g → bb splitting in the initial state. The fact that g → bb splittings are accounted for by the matrix elements guarantees a more precise modelling of the spectator b quark, while the simultaneous presence of W t and tt channels, starting from LO, ensures a perturbatively stable description of both contributions, as well as a NLO accurate prediction for their interference.
A generator based on the POWHEG method and pp → + ν l −ν l bb matrix elements at NLO in the 5FNS has been presented in Ref. [51]. However, the matching of parton showers to matrix elements that involve top-quark resonances poses nontrivial technical and theoretical problems [52] that have not been addressed in Ref. [51] and which cannot be solved within the original formulations of the POWHEG or MC@NLO methods. The problem is twofold. On the one hand, when interfacing a generator to a shower, if we do not specify which groups of final-state particles arise from the decay of the same resonance, the recoil resulting from shower emissions leads to arbitrary shifts of the resonance invariant masses, whose magnitude can largely exceed the top-quark width, resulting in unphysical distortions of the top line shape [52]. On the other hand, in the context of the infraredsubtraction and matching procedures, the standard mappings that connect the Born and real-emission phase spaces affect the top resonances in a way that drastically deteriorates the efficiency of infrared (IR) cancellations and jeopardises the consistency of the matching method [52].
A general NLO+PS matching technique that allows for a consistent treatment of resonances has been introduced, and applied to t-channel single-top production, in Ref. [52]. This approach will be referred to as resonance-aware matching. It is based on the POWHEG 4 method and is implemented in the POWHEG-BOX-RES framework, which represents an extension of the POWHEG-BOX [54]. In this framework each component of the cross section (i.e. Born, virtual and real) is separated into the sum of contributions that are dominated by well-defined resonance histories, such that in the narrow-width limit each parton can be uniquely attributed either to the decay products of a certain resonance or to the production subprocess. Within each contribution the subtraction procedure is organized in such a way that the off-shellness of resonant s-channel propagators is preserved, and resonance information on the final-state particles can be communicated to the shower program that handles further radiation and hadronization. This avoids uncontrolled resonance distortions, ensuring a NLO accurate description of the top line shape. The resonance-aware approach also improves the efficiency of infrared subtraction and phase-space integration in a dramatic way.
In this paper we present a NLO+PS generator, that we dub bb4l in the following, based on NLO matrix elements for pp → + ν l −ν l bb in the 4FNS matched to Pythia8 [55,56] using the resonance-aware POWHEG method. This new generator combines, for the first time, the following physics features: -consistent NLO+PS treatment of top resonances, including quantum corrections to top propagators and off-shell top-decay chains; -exact spin correlations at NLO, interference between NLO radiation from top production and decays, full NLO accuracy in tt production and decays; -unified treatment of tt and W t production with interference at NLO; -improved modelling of b-quark kinematics thanks to b-quark mass effects; -access to phase-space regions with unresolved b quarks and/or jet vetoes.
These improvements are of particular interest for precision top-mass measurements, for W t analyses, and for top backgrounds in the presence of jet vetoes or in the off-shell regime. Technically, the bb4l generator is based on OpenLoops [57] matrix elements. To this end we have developed a general and fully-flexible POWHEG-BOX+OpenLoops interface, which allows one to set up NLO+PS generators for any desired process. The paper is organized as follows. In Section 2 we briefly review the resonance-aware matching method. In Section 3 we discuss new developments in the POWHEG-BOX-RES framework that have been relevant for the present work. In Section 4 we discuss various aspects of the bb4l generator, including scope, usage, interface to Pythia8, and consistency checks. In Section 5 we detail the setup employed for the phenomenological studies presented in the subsequent sections. There we compare the bb4l generator to the previously available POWHEG generators, the hvq and ttb NLO dec ones, and we present technical studies that show the impact of the resonance-aware matching and of other improvements implemented in bb4l. Specifically, in Section 6 we consider observables that are directly sensitive to top-quark resonances and top-decay products, while in Section 7 we investigate the + ν l −ν l bb cross section in the presence of jet vetoes that enhance its single-top content. Our conclusions are presented in Section 8.
The POWHEG-BOX-RES framework together with the bb4l generator can be downloaded at http://powhegbox.mib.infn.it.

Resonance-aware subtraction and matching
In the following we recapitulate the problems that arise in processes where intermediate narrow resonances can radiate as they decay, and summarize the ideas and methodology behind the resonance-aware algorithm of Ref. [52]. We refer the reader to the original publication for the description of the method in full detail.
Commonly used IR subtraction methods for the calculation of NLO corrections [58][59][60] are based upon some procedure of momentum reshuffling for the construction of collinear and infrared counterterms. More specifically, given the kinematics of the real-emission process, and having specified a particular collinear region (i.e. a pair of partons that are becoming collinear), there is a well-defined mapping that constructs a Born-like kinematic configuration (called the "underlying Born" configuration) as a function of the real one. The mapping is such that, in the strict collinear limit, the Born configuration is obtained from the real one by appropriately merging the collinear partons. In the traditional methods, these mappings do not necessarily preserve the virtuality of possible intermediate s-channel resonances. If we consider the collinear region of two partons arising from the decay of the same s-channel resonance, the typical difference in the resonance virtuality between the real kinematics and the underlying-Born one is of order m 2 /E, where m is the mass of the two-parton system, and E is its energy. Because of this, the cancellation between the real contribution and the subtraction term becomes effective only if m 2 /E < Γ, where Γ is the width of the resonance. As long as Γ is above zero, the traditional NLO calculations do eventually converge, thanks to the fact that in the strict collinear limit the cancellation takes place. However, convergence becomes more problematic as the width of the resonance decreases.
The presence of radiation in resonance decays causes even more severe problems in NLO+PS frameworks. In POWHEG, radiation is generated according to the formula The first term in the square bracket corresponds to the probability that no radiation is generated with hardness above an infrared cutoff q cut , and its kinematics corresponds to the Born one. Each α in the sum labels a collinear singular region of the real cross section. The full real matrix element is decomposed into a sum of terms where each R α is singular only in the region labelled by α. The real phase space Φ α (Φ B , Φ rad ) depends upon the singular region α and is given as a function of the Born kinematics Φ B and three radiation variables Φ rad . The inverse of Φ α implements the previously mentioned mapping of the real kinematics into an underlying Born one. Thus, for a given Φ B and Φ rad , each term in the sum inside the square bracket in Eq. (2.1) is associated with a different real phase-space point. For each α, k α T is defined as the hardness of the collinear splitting characterized by the kinematics Φ α (Φ B , Φ rad ). It usually corresponds to the relative transverse momentum of the two collinear partons.
The Sudakov form factor, ∆, is such that the square bracket in Eq. (2.1), after performing the integrals in dΦ rad , becomes exactly equal to one (a property sometimes called unitarity of the real radiation). In general we have In order to achieve NLO accuracy, theB(Φ B ) factor must equal the NLO inclusive cross section at given underlying Born kinematics, where both the second and third term on the right hand side are infrared divergent, but the sum, being an inclusive cross section, is finite. The cancellation of singularities is achieved with the usual subtraction techniques.
We are now in a position to discuss the problems that arise in processes with radiation in decays of resonances. In order to do this, we focus on the W − W + bb production process. As an example of the problem, we consider a real emission contribution where a gluon g is radiated, such that the mass of the W + bg and W −b systems are very close to the top nominal mass. We call α b the singular region corresponding to b and g, and αb the region corresponding to theb and g becoming collinear, respectively. If we consider the case when the b andb partons are relatively close in direction, as g becomes collinear to the b or theb parton, two components will dominate the real cross section, R α b and R αb , in a proportion that is determined by how close the gluon is to the b or to theb partons. If the gluon is not much closer to the b region with respect to theb one, the R αb contribution will be comparable or larger than the R α b one. We now observe that, for the same real kinematic configuration, we have two singular regions and two corresponding underlying-Born configurations. In the α b singular region, the underlying Born is obtained by merging the bg system into a single b, while in the αb region it is thebg system that is merged into a singleb. It is therefore clear that, in the α b merging, the resonance virtualities are nearly preserved in the underlying Born, while in the αb one the resonances will be far off-shell. The R αb /B terms appearing both in Eq. (2.1) and (2.4) will become very large, the top resonances being on-shell in the numerator and off-shell in the denominator. However, in the POWHEG framework, these ratios should be either small (of order α s ) or should approach the Altarelli-Parisi splitting functions for the method to work.
It is thus clear that, if resonances are present, the traditional decomposition into singular regions must be revised. In particular, each α should become associated to a specific resonance structure of the event, such that collinear partons originate from the same resonance. Furthermore, the phase space mapping Φ α (Φ B , Φ rad ) should preserve the virtuality of the intermediate resonances. This is, in brief, what was done in Ref. [52].
The resonance-aware formalism also offers the opportunity to modify and further improve the POWHEG radiation formula. We make, for the moment, the assumption that each decaying resonance has only one singular region, and the radiation not originating from a resonance decay also has only one singular region. This is the case, for example, for the resonance structure of the process gg → (t → W + b)(t → W −b ), since in POWHEG the initial-state-radiation (ISR) regions are combined into a single one. We consider the formula where, by writing Φ α rad , we imply that the radiation variables are now independent for each singular region. By expanding the product, we see that we get a term with no emissions at all, as in Eq. (2.1), plus terms with multiple (up to three) emissions. It can be shown that, as far as the hardest radiation is concerned, formula (2.6) is equivalent to formula (2.1). To this end, one begins by rewriting Eq. (2.6) as a sum of three terms, with appropriate θ functions such that each term represents the case where the hardest radiation comes from one of the three regions. It is easy then to integrate in each term all radiations but the hardest, thus recovering the full Sudakov form factor appearing in the second term in the square bracket of Eq. (2.1).
The bb4l generator can generate radiation using the improved multiple-radiation scheme of formula (2.6) or the conventional single-radiation approach of Eq. (2.1). In events generated with multiple emissions included, the hardest radiation from all sources (i.e. production, t andt decays) may be present. The POWHEG generated event is then completed by a partonic shower Monte Carlo program that attaches further radiation to the event. The interface to the shower must be such that the shower does not generate radiation in production, in t decay and int decay that is harder than the one generated by POWHEG in production, t andt decay, respectively. 5

The POWHEG-BOX-RES framework
In this section we illustrate features that have been added to the POWHEG-BOX-RES package since the publication of Ref. [52], and discuss some issues that were not fully described there.

Automatic generation of resonance histories
In the POWHEG-BOX-RES implementation of Ref. [52], the initial subprocesses and the associated resonance structures were set up by hand. We have now added an algorithm for the automatic generation of all relevant resonance histories for a given process at a specified perturbative order. Thanks to this feature, the user only needs to provide a list of subprocesses, as was the case in the POWHEG-BOX-V2 package. This is a considerable simplification, in view of the fact that, when electroweak processes are considered, the number of resonance histories can increase substantially. Details of this feature are given in Appendix A.1.

Colour assignment
Events that are passed to a shower generator for subsequent showering must include colourflow information in the limit of large number of colours. In the POWHEG-BOX-V2 framework, colours are assigned with a probability proportional to the corresponding component of the colour flow decomposition of the amplitude. The extension of this approach to the POWHEG-BOX-RES framework requires some care due to possible inconsistencies between the colour assignment and the partitioning into resonance histories. This issue and its systematic solution are discussed in detail in Appendix A.2.

POWHEG+OpenLoops interface
All tree and one-loop amplitudes implemented in the bb4l generator are based on the OpenLoops program [57] in combinations with COLLIER [61] or CutTools [62] and OneLOop [63]. In the framework of the present work a new general process-independent interface between the POWHEG-BOX and OpenLoops has been developed. It allows for a straightforward implementation of a multitude of NLO multi-leg processes matched to parton showers including QCD and, in the future, also NLO electroweak corrections [64,65]. Technical details and a brief documentation of this new interface can be found in Appendix A.3.

Description of the generator
The implementation of combined off-shell tt and W t production in the POWHEG-BOX-RES framework presented in this paper is based on all possible Feynman diagrams contributing to the process pp → + ν l −ν l bb+X at NLO accuracy in QCD, i.e. up to order α 3 S α 4 EM . All bottom-mass effects have been fully taken into account and for the consistent treatment of top-, W -, and Z-resonances at NLO we rely on the automated implementation of the complex-mass scheme [66,67] within OpenLoops.

Resonance histories
The Born level resonance structure for pp → + ν l −ν l bb at O(α 2 S α 4 EM ) is actually very simple. Indeed, it is sufficient to consider two kinds of resonance histories. In Fig. 1 we show two corresponding Feynman diagrams for the process pp → µ + ν µ e −ν e bb . Internally, according to the POWHEG-BOX-RES conventions [52], the resonance histories are described by the arrays In flavres, for each particle, we give the position of the resonance from which it originates. For partons associated with the production subprocess flavres is set to zero.  Figure 1. Sample Feynman graphs corresponding to the two resonance histories relevant for pp → µ + ν µ e −ν e bb production. Figure 2. Representative Born diagram for W t production.
The resonance structures that differ only by the external parton flavours are collected into resonance groups, so that, in the present case, we have only two resonance groups. We remark that there is no need of a unique correspondence between resonance structures and possible combinations of resonant propagators in individual Feynman diagrams. What is required is that all resonances present in any given Feynman graph are also present in an associated resonance structure, but not vice versa. For example, in the present implementation of the bb4l generator the consistent treatment of single-top topologies like the one in Fig. 2 is guaranteed through resonance histories of tt type (flav_1,flavres_1), which involve an additionalt →bW − resonance. This does not lead to any problems, since the corresponding subtraction kinematics, which preserves the mass of thebW − system, is perfectly adequate also for single-top topologies.
The POWHEG-BOX-RES code automatically recognizes resonance histories that can be collected into the same resonance group. It also includes a subroutine for the automatic generation of an adequate phase-space sampling for each resonance group. In this context, rather than relying upon standard Breit-Wigner sampling, care is taken that also the offshell regions are adequately populated. This is essential in resonance histories of the kind shown in the right graph of Fig. 1, where the generation of the W virtualities according to their Breit-Wigner shape would well probe the region where an off-shell Z decays into two on-shell W 's, but not the regions where an on-shell Z decays into an on-shell W and an off-shell one. It also guarantees that cases like the diagram in Fig. 2 are properly sampled. The interested reader can find more technical details by inspecting the code itself.

The complex-mass scheme
In our calculation all intermediate massive particles are consistently treated in the complexmass scheme [66,67], where the widths of unstable particles are absorbed into the imaginary part of the corresponding mass parameters, This choice implies a complex-valued weak mixing angle, and guarantees gauge invariance at NLO [67].

The decoupling and MS schemes
When performing a fixed-order calculation with massive quarks, one can define two consistent renormalization schemes that describe the same physics: the usual MS scheme, where all flavours are treated on equal footing, and a mixed scheme [68], that we call decoupling scheme, in which the n lf light flavours are subtracted in the MS scheme, while heavy-flavour loops are subtracted at zero momentum. In this scheme, heavy flavours decouple at low energies.
In the calculation of the + ν l −ν l bb hard scattering cross section we treat the bottom quark as massive and, correspondingly, n lf is equal to four. The renormalization of the virtual contributions is performed in the decoupling scheme with a four-flavour running α S . For consistency, the evolution of parton distribution functions (PDFs) should be performed with four active flavours, so that, in particular, no bottom-quark density is present and no bottom-quark initiated processes have to be considered. However, given that the process at hand is characterised by typical scales far above the b-quark threshold, it is more convenient to convert our results to the MS scheme in such a way that they can be expressed in terms of the MS strong coupling constant, running with five active flavours, and also with fiveflavour PDFs.
The procedure for such a switch of schemes is well known, and was discussed in Ref. [69]. For + ν l −ν l bb production, we need to transform the qq and gg squared Born amplitudes B qq and B gg , computed in the decoupling scheme, in the following way 3) where µ R and µ F are the renormalization and factorization scales, respectively, and m b is the bottom-quark mass. The contribution of the b parton densities, that are present in the five-flavour scheme, should not be included in this context.

The virtual corrections
The virtual contributions have been generated using the new interface of the POWHEG-BOX with the OpenLoops amplitude generator, as described in Appendix A.3. While OpenLoops guarantees a very fast evaluation of one-loop matrix elements, the overall efficiency of the generator can be significantly improved by minimising the number of phase space points that require the calculation of virtual contributions. As detailed in Appendix A.4, this is is achieved by evaluating the virtual and real-emission contributions with independent statistical accuracies optimised according to the respective relative weights. Moreover, when generating events, a reweighting method can be used in order to restrict virtual evaluations to the small fraction of phase space points that survive the unweighting procedure.

Interface to the shower
The generator presented in this work shares many common features with the one of Ref. [35].
In particular, in both generators, Les Houches events include resonance information, and an option for a multiple radiation scheme is implemented, denoted as allrad scheme, according to the corresponding powheg.input flag. As explained in Ref. [35] and reviewed in Section 2, when this scheme is activated, the mechanism of radiation generation is modified. Rather than keeping only the hardest radiation arising from all singular regions, the program stores several "hardest radiations": one that takes place at the production stage, and one for the decay of each resonance that can radiate. All these radiations are assembled into a single Les Houches event. Thus, for example, in events with the t andt resonances, one can have up to three radiated partons: one coming from the initial-state particles, one arising from the b in the t-decay, and one from theb in thet-decay. When generating fully showered events, the hardness 6 of the shower must be limited in a way that depends upon the origin of the radiating parton. If the radiating parton is not son of a resonance, the hardness of the shower arising from it must be limited by the hardness of the Les Houches radiation that arises in production. 7 Radiation arising from partons originating from a resonance must have their hardness limited by the hardness of the parton radiated from the resonance in the Les Houches event. This requires a shower interface that goes beyond the Les Houches approach. In Ref. [35] a suitable procedure has been conceived and implemented in Pythia8 [55,56]. The interested reader can find all details in the Appendix A of Ref. [35]. In essence, the procedure was to examine the showered event, compute the transverse momentum of Pythia8 radiation in top decays, and veto it if higher than the corresponding POWHEG one. Vetoing is performed by rejecting 6 Here and in the following by hardness we mean the relative transverse momentum of two partons arising from a splitting process, either in initial-or in final-state radiation. 7 By radiation in production we mean any radiation that does not arise from a decaying resonance. This can be initial-state radiation, but also radiation from final-state partons, as in the right diagram in Fig. 1 and the one in Fig. 2, where the b's do not arise from a decaying resonance.
the showered event, and generating a new Pythia8 shower, initiated by the same Les Houches event. This procedure was iterated until the showered event passes the veto. In the present work, we have adopted this procedure in order to make a more meaningful comparison with the results of Ref. [35]. However, we have also verified that, by using Pythia8 internal mechanism for vetoing radiation from resonance decay, we get results that are fully compatible with our default approach. 8 This aspect and the comparison among the two methods are shown in Appendix B.2.

Traditional NLO+PS matching
It is possible to run our new generator in a way that is fully equivalent to a standard POWHEG matching algorithm (as implemented in the POWHEG-BOX-V2) ignoring the resonance structure of the processes. This is achieved by including the line nores 1 in the powheg.input file. 9 Such an option is implemented only for the purpose of testing the new formalism with respect to the old one. It turns out that, in the nores 1 mode, the program has much worse convergence properties, most likely because of the less effective cancellation of infrared singularities mentioned in Section 2. We find, for example, that in runs with equal statistics (with about 15 million calls) the absolute error in the nores 1 case is roughly seven times larger than in the nores 0 (default) case. The generation of events also slows down by a similar factor.
We stress again that, in the limit of small widths, the NLO+PS results obtained in the nores 1 mode are bound to become inconsistent, as discussed in Section 2 and, more extensively, in Ref. [52].

Consistency checks
At the level of fixed-order NLO calculations, the traditional machinery of the POWHEG-BOX is well tested and we trust corresponding results to be correct. On the other hand, the NLO subtraction procedure implemented in the POWHEG-BOX-RES code is substantially different and still relatively new. As was done in Ref. [52] for t-channel single-top production, also for the + ν l −ν l bb production presented here, we systematically validated the fixed-order NLO results obtained with the POWHEG-BOX-RES implementation by switching on and off the generation of resonance structures. We found perfect agreement between the two calculations.
Additionally, we performed a detailed comparison against the fixed-order NLO results of Ref. [41] and found agreement at the permil level. Furthermore, via a numerical scan in the limit of the top width going to zero, Γ t → 0, we verified that any α S log (Γ t ) enhanced terms in the soft-gluon limit successfully cancel between real and virtual contributions. This last test was performed for various light-and b-jet exclusive distributions which are subject to sizable non-resonant/off-shell corrections.

Phenomenological setup
In this section we document the input parameters, acceptance cuts and generator settings that have been adopted for the numerical studies presented in Section 6. Moreover we introduce a systematic labelling scheme for the various NLO+PS approximations that are going to be compared.

Input parameters
Masses and widths are assigned the following values The electroweak couplings are derived from the gauge-boson masses and the Fermi constant, where µ W and µ Z are complex masses given by Eq. (4.1). The value of the top-quark width we use is consistently calculated at NLO from all other input parameters by computing the three-body decay widths Γ(t → ff b) into any pair of light fermions f andf and a massive b quark. To this end, we employ a numerical routine of the MCFM implementation of Ref. [34].
As parton distributions we have adopted the five-flavour MSTW2008NLO PDFs [71], as implemented in the Ref. [72], with the corresponding five-flavour strong coupling constant, and for their consistent combination with four-flavour scheme parton-level cross sections the scheme transformation of Section 4.3 was applied. In the evaluation of the matrix elements, only the bottom and the top quarks are massive. All the other quarks are treated as massless. In addition, the Cabibbo-Kobayashi-Maskawa matrix is assumed to be diagonal.
When generating events we adopt the following scale choice: • For resonance histories with a top pair we use where the (anti)top masses and transverse momenta are defined in the underlying Born phase space in terms of final state (off-shell) decay products.
• For resonance histories with an intermediate Z we use In addition, we set the value of the POWHEG-BOX parameter hdamp to the mass of the top quark. This setting yields a transverse-momentum distribution of the top pair that is more sensitive to scale variations and more consistent with data at large transverse momenta. It only affects initial-state radiation. For a detailed description of this parameter, we refer the reader to Ref. [73].

Pythia8 settings
We interface our POWHEG generator to Pythia8.1, 10 as illustrated in Appendix A of Ref. [35], and so we perform the following Pythia8 calls: pythia.readString("SpaceShower:pTmaxMatch = 1"); pythia.readString("TimeShower:pTmaxMatch = 1"); pythia.readString("PartonLevel:MPI = off"); pythia.readString("SpaceShower:QEDshowerByQ = off"); pythia.readString("SpaceShower:QEDshowerByL = off"); pythia.readString("TimeShower:QEDshowerByQ = off"); pythia.readString("TimeShower:QEDshowerByL = off"); The first two calls are required when interfacing Pythia8 to NLO+PS generators. The third call switches off multi-parton interactions and it is only invoked for performance reasons: in fact, the shower of the events is faster when multi-parton interactions are not simulated. The remaining calls switch off the electromagnetic radiation in Pythia8. This makes it easier to reconstruct the W boson momentum, since we do not need to dress the charged lepton, from vector boson decay, with electromagnetic radiation. These settings are appropriate in the present context since we do not make any comparison with data.
Pythia8 provides by default matrix-element corrections [74]. In our case, they are relevant for radiation in the top decays, which are corrected using t → W bg tree level matrix elements. These corrections are also applied in subsequent emissions in order to better model radiation from heavy flavours in general. If not explicitly stated otherwise, we include the following setup calls pythia.readString("TimeShower:MEcorrections = on"); pythia.readString("TimeShower:MEafterFirst = on"); These corrections never modify the Les Houches event weight. They only affect the radiation generated by the shower. Thus, leaving these flags on does not lead to over-counting. If the second flag is off, matrix-element corrections are applied only to the first shower emission. If it is on, they are also applied to subsequent radiation. In fact, even if these corrections cannot fully account for the structure of the matrix elements, they at least better account for mass effects arising in radiation from the off-shell top quarks and from the massive final-state b's.
In our analysis, we keep B hadrons stable, performing the corresponding Pythia8 setup calls. Aside from these, all remaining settings are left to the defaults of Pythia8.1.

Generators and labels
In Section 6 we compare three different generators that implement an increasingly precise treatment of tt production and decay: • the hvq generator of Ref. [20]; • the ttb NLO dec generator of Ref. [35]; • the new bb4l generator, which we consider as our best prediction.
The main physics features of the various generators and the labels that will be used to identify the corresponding predictions are listed in Table 1 In order to quantify the impact of various aspects of the resonance-aware approach, in Section 6 we will compare various settings of the bb4l generator where some resonanceaware improvements are turned on and off or are replaced by certain approximations. Specifically, the following settings will be considered: (a) the resonance-aware formalism is switched on with default settings; (b) the resonance-aware formalism is switched off, which corresponds to using the traditional POWHEG approach; (c) the resonance-aware formalism is switched off, but a resonance assignment is guessed based on the kinematic structure of the events, according to the method described in Appendix B.1; (d) the resonance-aware formalism is switched on, but, instead of applying the multipleradiation scheme of Eq. (2.6), only a single radiation is generated with POWHEG according to Eq.  (e) same as (d), but the resonance information is stripped off in the POWHEG Les Houches event file before passing it to the showering program.
The various bb4l settings and corresponding labels are summarised in Table 2.

Physics objects
In the subsequent sections we study various observables defined in terms of the following physics objects.
(a) We denote as B andB hadron the hardest b-flavoured andb-flavoured hadron in the event.
(b) Final-state hadrons are recombined into jets using the FastJet implementation [75] of the anti-k T jet algorithm [76] with R = 0.5.
(c) We denote as b-jet (j B ) and anti-b-jet (jB) the jet that contains the hardest B and B hadron, respectively. When examining results obtained with the hadronization switched off, jets are b-tagged based on b quarks rather than B hadrons.
(d) Leptons, neutrinos and missing transverse energy are identical to their corresponding objects at matrix-element level, since we switched off QED radiation and hadron decays in Pythia8.
(e) Reconstructed W + and W − bosons are identified with the corresponding off-shell lepton-neutrino pairs in the hard matrix elements. 11 (f) Reconstructed top and anti-top quarks are defined as off-shell W + j B and W − jB pairs, respectively, i.e. b-jets and W -bosons are matched based on charge and b-flavour information at Monte-Carlo truth level. The same approach is used for + j B and l − jB pairs. 11 Similarly as for top resonances, also W resonances are identified with their off-shell decay products according to the resonance history of the event at hand. This information is written in the shower record and propagated through the shower evolution. In this way, possible QED radiation off charged leptons is included into the W -boson momentum at Monte Carlo truth level. However, since electromagnetic radiation from Pythia8 is turned off in our analysis, each W boson coincides with a bare lepton-neutrino system.
Unless stated otherwise, in kinematic distributions we always perform an average over the t andt case (thus also on lepton-antilepton, b-anti-b, etc.). The top-pair observables in Sections 6.2-6.3 are computed by requiring the presence of a b and ab jet with and applying the following leptonic cuts, where l = + , l − and p miss T is obtained from the vector sum of the transverse momentum of the neutrinos in the final state.

Top-pair dominated observables
Here we present numerical predictions for pp → e + ν e µ −ν µ bb + X at √ s =8 TeV. In particular, we study various observables that are sensitive to the shape of top resonances.

Comparison with traditional NLO+PS matching
In the following, we compare nominal bb4l predictions, generated with default settings, with results obtained by switching off the resonance-aware formalism (i.e. setting the flag nores to 1). In this way we get results that are fully equivalent to a POWHEG-BOX-V2 (or "traditional") implementation. For this comparison we do not impose any cuts, i.e. we perform a fully inclusive analysis that involves, besides tt production, also significant contributions from W t single-top production.
Events generated with the traditional implementation do not contain any information whatsoever about their resonance structures. We label the curves obtained by showering these events as res-off. Because the resonance information is not available, the shower generator will not preserve the virtualities of the resonances. In order to further explore the usability of the res-off results, we also consider the possibility of reconstructing the resonance information of the Les Houches event on the basis of its kinematic proximity to one of the possible resonant configurations. Specifically, we perform an educated guess of the resonance structure of the event, assigning it to a tt or to a Z resonance configuration (see Section 4.1), and assigning the radiation either to the initial state or to the outgoing b's. The curves obtained this way are labelled as res-guess and the procedure for reconstructing the resonance information from the event kinematics is detailed in Appendix B.1.
We first consider, in Fig. 3, the invariant mass of the W j B and of the lj B systems. In the res-off case, we observe that the reconstructed mass peak has a wider shape. This is expected, since neither the POWHEG-BOX nor the shower program preserve the virtuality of the top resonances. In the res-guess case the width of the peak is diminished, although not quite at the level of the resonance-aware prediction, labelled as res-default. We also observe a mild shift in the peak in the res-guess case, which improves the agreement with the resdefault result. The distribution in the mass of the lepton-j B system also shows marked  We compare our default resonance-aware predictions (res-default) against the "traditional", i.e. resonance-unaware, implementation (res-off) and a prediction where the event-by-event resonance information is obtained from a guess based on kinematics. In the ratio plot we illustrate relative deviations with respect to res-default.
differences in shape in the region that is most relevant for a top-mass determination, with more pronounced differences in the res-off case.
The above findings suggest that the width of the peak is determined both by the shower generator being aware of the resonances in the Les Houches event, and by the hardest radiation generation being performed in a way that is consistent with the resonance structure. In order to assess the effects that originate solely from resonance-aware matching and showering in a more accurate way, in Fig. 4 we disable the multiple radiation scheme  . Invariant mass of the W j B system obtained with the bb4l generator. We compare our resonance-aware predictions without employing the multiple radiation scheme (res-singlerad) against the "traditional", i.e. resonance-unaware, implementation (res-off) and a prediction where any resonance information is stripped off the Les Houches event file (res-strip). In the ratio plot we illustrate relative deviations with respect to res-singlerad.
of Eq. (2.6) (by setting allrad 0) and compare the resulting resonance-aware predictions (res-singlerad) against the cases where resonance information is removed from the Les Houches event before showering (res-strip) or the case where the resonance-aware system is completely switched off (res-off). We find that the res-strip result lies between the ressinglerad and the res-off ones, somewhat closer to the latter, and the differences between the various predictions are considerable. Therefore, we conclude that the observed widening of the peak in Figures 3-4 can be attributed to both shortcomings of a resonance unaware parton shower matching: the parton shower reshuffling not preserving the resonance masses, and the uncontrolled effects of resonances at the level of the first emission in the traditional POWHEG approach. In Fig. 5 we display the j B mass and profile, defined as This observable corresponds to the cross section weighted by the fraction of the total hadronic transverse momentum of the particles contained in a given cone around the jet axis, with respect to the transverse momentum of the b-jet. Again we observe marked differences among the res-default and the res-off results, and, to a lesser extent, between the res-default and res-guess ones. Both plots suggest that in the res-off case there is less activity around the B hadron, leading to smaller jet masses and to a slightly steeper jet profile. The particularly pronounced shape distortion of the j B mass plot near 10 GeV in the res-guess case can be tentatively attributed to the transition from the region where radiation (generated with the traditional method) does not change the mass of the resonance by an amount comparable or larger than its width, to the region where it does, so that we see the difference between the res-guess and res-default results grow with larger jet masses.
In Fig. 6 we compare the B fragmentation function and the B-hadron transverse momentum computed in the reconstructed top-decay rest frame. The x B variable is defined as the B energy in the reconstructed top rest frame normalized to the maximum value that it can attain at the given top virtuality, while p B T,dec is the transverse momentum of the B relative to the recoiling W in the same frame. We find marked differences also for these distributions. While in the case of the p B T,dec variable we see a reasonable consistency  Figure 6. B fragmentation function and B-hadron transverse momentum in the top decay frame. Absolute predictions and ratios as in Fig. 3.
between the res-guess and res-default results, the agreement deteriorates in the case of the fragmentation function. We conclude that the consistent treatment of resonances implemented in the bb4l generator yields a narrower peak for the reconstructed top distribution with respect to a traditional (resonance-blind) NLO+PS matching approach. Furthermore, a large part of the difference is not related to the lack of resonance information at the level of the shower generator, and thus cannot be reduced by using a more sophisticated interface to the shower based on a resonance-guessing approach of kinematic nature.

Comparison with the ttb NLO dec generator
In this section we compare the bb4l generator against the ttb NLO dec generator of Ref. [35]. The standard tt cuts of Eqs. (5.9)-(5.10) are applied throughout. We examined a large set of distributions, but here we only display the most relevant ones, and those that show the largest discrepancies.
We begin by showing in Fig. 7 the invariant mass distribution of the W j B and lj B systems. We observe remarkable agreement between the bb4l and ttb NLO dec generators, especially in the description of the reconstructed top peak and of the shoulder in the lepton-j B invariant mass. This agreement is quite reassuring. In fact, in the ttb NLO dec generator, the separation of radiation in production and resonance decay is unambiguous, while in bb4l it is based on a probabilistic approach according to a kinematic proximity criterion. Thus, in the light of Fig. 7, the former generator supports the method of separation of resonance histories adopted by the latter. On the other hand, off-shell and non-resonant effects are implemented in the ttb NLO dec generator in LO approximation, by reweighting the on-shell result. Thus the bb4l results support the validity of this approximation in the ttb NLO dec implementation. As an indicative estimate of the potential implications for precision m t determination, we have determined that in a window of ±30 GeV around the peak of the W j B distributions, the average W j B mass computed with the ttb NLO dec generator is roughly 0.1 GeV smaller than the one from bb4l.    The NLO distribution in the mass of the reconstructed top was also examined in Ref. [35] (sec. 3.2, Fig. 3). There, the ttb NLO dec fixed-order NLO result was compared to the fixed-order NLO result of Ref. [38], and the former was found to be enhanced by about 10% in a region of roughly 1 GeV around the peak. This comparison was carried out with massless b quarks, since mass effects were not available in Ref. [38]. We computed the same distribution and carried out the same NLO comparison, using however the bb4l generator instead of the result of Ref. [38] and taking into account b-mass effects. Again, we find the same enhancement in the ttb NLO dec NLO result. However, in the fully showered result we see instead a small suppression of the peak in the ttb NLO dec relative to the bb4l generator, suggesting that the NLO difference tends to be washed out by showering effects.
We examined several distributions involving b-jets (here again we average over the b-andb-jet contributions). We found no appreciable difference for the b-jet transverse momentum, while we did find significant differences in the jet mass and the jet profile, displayed in Fig. 8. Both plots indicate that the bb4l generator yields slightly wider b-jets as compared to the ttb NLO dec one.   In Fig. 9 we plot the B fragmentation function and the p B T,dec observables. We find that the fragmentation function is slightly harder, and the p B T,dec distribution is slightly softer in the bb4l case. Again, this is consistent with the observation of slightly reduced radiation from b's in the bb4l case. We have verified that this feature persists also when hadronization is switched off in Pythia8.
Although the differences in the b-jet structure are quite significant, they are not sufficient to induce an observable shift in the reconstructed mass peak. This could only happen if the difference in the jet profile caused a consistent difference in the jet energy, due to energy loss outside the jet-cone. This does not seem to be the case since the jet profiles become similar in the two generators already for ∆ R < 0.5.

Comparison with the hvq generator
In this section we compare the bb4l generator against the hvq generator of Ref. [22], which is based on on-shell NLO matrix elements for tt production. Again the standard tt cuts of Eqs. (5.9)-(5.10) are applied throughout. The W j B and lj B mass distributions, shown in Fig. 10, show reasonably good agreement between the two generators as far as the shape of the W j B peak and of the lj B shoulder are concerned. However, for large top virtualities, i.e. in the tails of both distributions, sizable differences can be appreciated. As we will see below, such differences originate from the fact that, in this region, the bb4l generator tends to radiate considerably less, which results in narrower b-jets as compared to the hvq generator. We note that the observed deviations with respect to the hvq generator are more drastic than the ones observed in section 6.2 for the ttb NLO dec generator. The m W j B distribution on the left of Fig. 10 additionally suggests a non-negligible shift in the reconstructed top mass between the two generators. In fact, we determined that in a window of ±30 GeV around the peak of the m W j B distributions, the average W j B mass computed with the hvq generator is roughly 0.5 GeV smaller than with the bb4l one.
In Fig. 11 we show distributions in the b-jet mass and profile, as defined in Eq. (6.1). Both plots indicate significantly narrower b-jets in the predictions obtained with the bb4l generator. Similarly, as shown in Fig. 12, the bb4l generator yields a harder B fragmen-  tation function and a softer p B T,dec distribution. The pattern we observe for the structure of b-jets is consistent with the fact that the bb4l generator has a reduced radiation in b-jets with respect to Pythia8. In the hvq generator, radiation from the b's is handled exclusively by Pythia8, while, in the bb4l generator, the hardest radiation from the b is handled by POWHEG. It should be stressed, however, that the B fragmentation function has a considerable sensitivity to the hadronization parameters. It would therefore be desirable to tune these parameters to B production data in e + e − annihilation, within the POWHEG framework, in order to perform a meaningful comparison.
In Fig. 13 we show a summary of the shape of the reconstructed top peak comparing each of the available POWHEG generators for tt production: bb4l, ttb NLO dec and hvq. We notice a fair consistency between the bb4l generator and the ttb NLO dec one, while larger deviations are observed comparing against hvq.  Figure 13. The W j B mass distribution near the top peak for the three generators bb4l (bb4 ), ttb NLO dec (tt ⊗ decay) and hvq (tt). In the ratio plot we illustrate relative deviations with respect to the bb4 prediction.

Jet vetoes and single-top enriched observables
In this section we investigate the behaviour of the bb4l generator in the presence of b-jet and light-jet vetoes. Such kinematic restrictions are widely used in order to reduce top backgrounds in H → W + W − studies and in many other analyses that involve charged leptons and missing energy. Also, jet vetoes play an essential role for experimental studies of W t single-top production [77,78]. In particular, the separation of W t and tt production typically relies upon the requirement that one large transverse-momentum b-jet is missing in the first process. From the theoretical point of view, the separation of W t and tt production is not a clear cut one, since the two processes interfere. As pointed out in the introduction, in the bb4l generator this problem is solved by providing a unified description of tt and W t production and decay, where also interference effects are included at NLO. Thus jet vetoes are expected to enrich the relative single-top content of bb4l samples, resulting in significant differences with respect to other generators that do not include W t contributions and interferences at NLO. The bb4l generator is particularly well-suited for the study of jet vetoes also because it includes b-mass effects, NLO radiation in top-production and -decay subprocesses, as well as resummation of multiple QCD emissions and hadronization effects as implemented in the parton shower. A first picture of the b-jet activity in the three generators, bb4l, ttb NLO dec and hvq (labelled according to Table 1 as bb4 , tt ⊗ decay and tt respectively), is provided by Fig. 14, where we compare NLO+PS distributions in the transverse momentum of the bjet. More precisely, the plotted observable corresponds to the sum of the b-andb-jet spectra and was computed in absence of any acceptance cut. Thus it involves potentially enhanced contributions from single-top topologies, which can lead to significant deviations between the tt prediction 12 and the ones that implement off-shell + ν l −ν l bb matrix elements. At large transverse momentum, the various predictions have rather similar shape, but the tt result features a clear deficit of about 10% with respect the bb4 and tt ⊗ decay ones. This can be attributed to the missing single-top contributions in the hvq generator. At high p T , thanks to the implementation of W t contributions via exact Born matrix elements for pp → + ν l −ν l bb, the tt ⊗ decay prediction is found to be in good agreement with the bb4 one. At small transverse momenta, the relative weight of W t production becomes more important, and the deficit of the tt prediction grows rather quickly, reaching up to 50% for very small transverse momenta. The tt ⊗ decay and bb4 predictions remain in good agreement down to p T,j B 10 GeV, but at smaller transverse momenta the tt ⊗ decay one develops a deficit that grows up to about 25%. This can be attributed, at least in part, to the increased importance of W t channels combined with the fact that these channels are not supplemented by an appropriate NLO correction in the tt ⊗ decay predictions. We also note that the discrepancy at hand can be interpreted as a kinematic shift of a few GeV only, while the enhancement of the resulting correction can be attributed to the pronounced steepness of the absolute p T,j B distribution in the soft region. Its sign is consistent with the fact that radiation arising from + ν l −ν l bb NLO matrix elements is expected to be rather soft in the presence of single-top contributions with initial-state collinear g → bb splittings, while in the tt ⊗ decay generator radiation is always emitted as if all b-quarks would arise from top decays, which results in a harder emission spectrum. The lower frame of Fig. 14 illustrates the relative importance of matching and shower effects in the bb4l generator, comparing against corresponding fixed-order NLO predictions. Again we observe nontrivial shape effects in the soft region. While they are not directly related to the differences observed in the middle frame, such effects highlight the importance of a consistent treatment of radiation and shower effects at small b-jet p T . On the other hand the good agreement between the tt ⊗ decay and bb4 predictions down to 10 GeV suggests that matching and pure shower effects are reasonably well under control in the bulk of the phase space.
Jet-binning and jet-veto effects are studied in Figures 15-16. For this analysis we apply again the lepton selection cuts of Eq. (5.10) and, at variance with the b-jet definition in Section 6, we identify as b-jets those jets containing at least a b-orb-flavoured hadron, irrespectively of its hardness. 13 Events are categorised according to the number of (light or heavy-flavour) jets, n j , and to the number of b-jets, n b , in the rapidity range |η| < 2.5, while we vary the jet transverse-momentum threshold p thr T,jet that defines jets. In Fig. 15, to investigate the effect of a b-jet veto, the integrated cross sections is plotted versus the jet-veto threshold, p thr T,jet . In the left plot the veto acts only on b-jets (n j ≥ n b = 0), while in the right plot a veto against light and b-jets is applied (n j = n b = 0). For p thr T,jet 80 GeV the vetoed cross section is dominated by tt production and quickly converges towards the inclusive result. In this region we observe few-percent level agreement between the tt ⊗ decay and bb4 predictions, while the on-shell tt prediction features a 10% deficit due to the missing single-top topologies. Reducing the jet-veto scale increases this deficit up to −30% in the case of the inclusive n b = 0 cross section. This finding is well consistent with the size of finite-width effects reported in Ref. [41]. In the case of the exclusive zero-jet cross section (n j = n b = 0, shown on the right) the deficit of the tt prediction is even more pronounced and reaches up to −50% at p thr T,jet = 10 GeV. Also the tt ⊗ decay results feature a similar, although less pronounced, deficit as the tt ones in the soft region. This can be attributed to the fact that initial-state radiation in both, the hvq and ttb NLO dec, generators is computed with on-shell tops, and thus overestimates the radiation produced near the single-top kinematic region.
Matching and pure shower effects are illustrated in the lower frames of Fig. 15. Both in the inclusive (n j > n b = 0) and exclusive (n j = n b = 0) case we observe that, down to 20 GeV, NLO+PS predictions feature an increasingly strong enhancement with respect to fixed-order ones. This can be attributed to shower-induced losses of b-jet transverse momentum. In the exclusive case (n j = n b = 0) this enhancement is somewhat milder, which we tentatively attribute to the interplay of parton shower radiation with the additional light-jet veto.
In Fig. 16 we plot the cross section with exactly one b-jet above the threshold p thr T,jet , i.e. we veto additional b-jets above this threshold. Again, inclusive results (n j ≥ n b = 1, shown on the left) are compared with exclusive ones (n j = n b = 1, shown on the right). The one-b-jet bin is typically used in W t single-top analyses. Similarly as for the zero-b-jet case, the difference between the bb4 and tt results points to an increasingly important single-top contribution at small p thr T,jet . Its quantitative impact is consistent with the fixed-order results of Ref. [41], and at p thr T,jet = 30 GeV it amounts to about 10% and 20%, respectively, in the inclusive and exclusive cases. Similarly as for the zero-b-jet case, tt ⊗ decay predictions feature a qualitatively similar but quantitatively less pronounced deficit with respect to the bb4 predictions. Matching and shower effects turn out to be rather mild in the inclusive case, probably due to the fact that the absolute distribution is not particularly steep in the limit of small transverse momentum. In contrast, the exclusive one-jet cross section (n j = n b = 1) is much more sensitive to the jet-veto scale, which leads to sizable matching and shower effects at small p thr T,jet . In summary, jet-vetoed cross sections can involve enhanced single-top contributions that are completely missing in the tt predictions obtained with the hvq generator while they 13 At fixed-order NLO, jet clustering and b-jet tagging are applied at parton level. are significantly underestimated in the tt ⊗ decay predictions of the ttb NLO dec generator, where single-top contributions are implemented via LO reweighting [35]. In practice such a reweighting approach ceases to work in phase-space regions far away from the doubleresonant region.

Conclusions
In this paper we have presented the first Monte Carlo generator that provides a fully consistent NLO+PS simulation of tt production and decay in the different-flavour dilepton channel, including all finite-width and interference effects. This new generator, dubbed bb4l, is based on the full NLO matrix elements for the process pp → + ν l −ν l bb. This guarantees NLO accuracy in tt production and decay, as well as the exact treatment of spin correlations and off-shell effects in top decay. Top resonances are dressed with quantum corrections, and also non-factorisable corrections associated with the interference of radiation in production and decays are taken into account. Bottom-quark masses are consistently included, which is quite important for the accurate modelling of b-quark fragmentation. Moreover, finite b-quark masses permit to avoid collinear singularities from initial-or final-state g → bb splittings. This allows for W + W − bb simulations in the full phase space, including regions with unresolved b quarks, which are indispensable for the simulation of top backgrounds in the presence of jet vetoes. It moreover provides a unified NLO description of tt and single-top W t production, including their quantum interference.
The technical problems that arise from infrared subtractions and NLO+PS matching in the presence of top-quark resonances are addressed by means of the fully general resonance-aware matching method that was proposed in Ref. [52] and implemented in the POWHEG-BOX-RES framework. This framework, besides allowing for a consistent matching to shower Monte Carlo generators, also ameliorates the efficiency of infrared subtraction and phase-space integration in a drastic way, and allows for a factorised treatment of NLO radiation in off-shell top production and decays. This represents a significant improvement (especially for what concerns top decays) with respect to the case where NLO+PS matching is applied to a single QCD emission.
Technically, the bb4l generator was realised by implementing OpenLoops matrix elements in the POWHEG-BOX framework. To this end we have developed a new and fully flexible interface, which allows one to set up POWHEG-BOX+OpenLoops NLO+PS generators for any desired process in a rather straightforward way.
We have carried out a thorough study of the impact of the resonance-aware method. To this end, we have compared our results with those obtained after disabling the resonanceaware formalism in such a way that the bb4l generator becomes fully equivalent to a traditional POWHEG-BOX-V2 implementation. On the one hand we observed that ignoring resonance structures can deteriorate the performance of the generator up to the point of rendering it unusable. On the other hand, we observed considerable distortions in the reconstructed mass of the top resonances with respect to the full resonance-aware result. In essence, the mass distribution becomes wider around the peak and slightly shifted. We were able to track the origin of these effects to two competing causes: the generation of radiation performed by POWHEG, that is considerably modified in the resonance-aware method, and the generation of radiation in the shower stage, where the shower Monte Carlo, being unaware of which groups of particles arise from the same resonance, tends to widen the resonance peaks. We have also shown that it does not seem to be possible to remedy this last problem by reconstructing the resonance structure on the basis of simple kinematic guesses.
Much attention was dedicated to the comparison of the new bb4l generator and the ttb NLO dec generator of Ref. [35]. Both are capable of handling NLO spin correlations and radiation in top decays. However, off-shell effects are only computed at LO in ttb NLO dec by reweighting the NLO cross section using the ratio of the full off-shell Born cross section divided by its zero-width approximation. These two generators are expected to provide similar results in the vicinity of top resonances. In fact, in this region, we find only modest differences between the two. In particular, the top virtuality distribution and distributions involving b jets are in reasonably good agreement. Slightly larger differences are found in distributions involving B hadrons, like for example, the B fragmentation function, in the top-decay frame.
A section of this work was dedicated to a comparison against the hvq generator, which has been heavily used by the LHC experimental collaborations for the generation of tt samples in both Run I and Run II. Close to the mass peak, bb4l and hvq predictions are fairly consistent, but the agreement is quickly spoiled as one moves towards off-shell regions. Furthermore, the ratio of the hvq to the bb4l results exhibits a negative slope across the resonance peak, and we found that the average virtuality of the top resonance in a window of ±30 GeV around the peak differs by about 0.5 GeV for the two generators. This calls for dedicated studies of the implications of resonance-aware matching in the context of precision m t -measurements. More sizable differences have been observed in the structure of the associated b-jets, the bb4l generator leading consistently to narrower jets and a harder fragmentation function for the associated B hadron. The above findings should be interpreted by keeping in mind that within the hvq generator radiation in top decays is solely handled by Pythia8, with matrix-element corrections turned on by default. These matrix-element corrections should improve the overall agreement between hvq and bb4l, and we have verified that disabling them leads to much more pronounced differences between the two generators.
We have included in this work an indicative comparative study of jet-veto effects when using the bb4l, ttb NLO dec and hvq generators. In the presence of jet vetoes, the hvq generator alone is clearly not adequate, since it misses the essential component of associated W t production. Perhaps surprisingly, it turns out that also the ttb NLO dec generator does not perform sufficiently well. Since W t production effects are included in this generator only at the level of a leading-order reweighting, we are led to conclude that the lack of NLO accuracy in the simulation of W t contributions limits the usability of the ttb NLO dec generator in single-top enriched regions. We stress however that the issue of jet-veto effects is complex, and deserves a dedicated future study.
The theoretical improvements implemented in the bb4l generator are relevant for phenomenological studies and experimental analysis that depends on the kinematic details of top-decay products. In particular, this new generator is ideally suited for precision determinations of the top-quark mass, for measurements of W t production, and for analyses where tt and W t production are subject to jet vetoes. The exact treatment of off-shell and non-resonant effects is also important for top backgrounds in Higgs and BSM studies based on kinematic selections with high missing energy or boosted bb pairs. part by the Swiss National Science Foundation (SNF) under contracts BSCGI0-157722 and PP00P2-153027, by the Research Executive Agency of the European Union under the Grant Agreement PITN-GA-2012-316704 (HiggsTools), and by the Kavli Institute for Theoretical Physics through the National Science Foundation's Grant No. NSF PHY11-25915. PN and CO would like to express a special thanks to the Mainz Institute for Theoretical Physics (MITP) for its hospitality and support while part of this work was carried out.

A Technical details
In this appendix we detail the technical improvements to the POWHEG-BOX-RES framework that have been implemented in order to allow for the implementation of pp → + ν l −ν l bb.

A.1 Automatic generation of resonance histories
The algorithm for finding the resonance histories is at present at an experimental level. It has been kept as simple and straightforward as possible in order to allow for future improvements and modifications. The algorithm begins with the lists of particle flavours specified in the user process routine init_processes, where the arrays flst_born and flst_real are filled. At variance with the POWHEG-BOX-V2 version, one also has to specify the length of each flavour list in the arrays. The lengths are stored in the arrays flst_bornlength and flst_reallength. For the process we are considering here (and in most cases) the lengths have all the same values (8 for the Born process and 9 for the real). At this stage, no resonance information is provided for the flavour lists, so the lists of resonance pointers (flst_bornres and flst_realres) remain initialized to zero, and the user does not need to modify them. The powers of the strong and weak coupling constants in the Born amplitudes (res_powst and res_powew) must instead be initialized by the user-process routines. At the moment we do not consider the possibility of having multiple Born-level processes with different orders of the strong and weak coupling constants. This may be required when considering mixed strong and electromagnetic radiation being generated with the POWHEG method, and will require minor modification of the code.
The algorithm proceeds recursively: intermediate particles are added at the end of the flavour list, and the pointers associated with the particles that arise from their splitting are appropriately set.
As an example, we consider the production of a W in association with a quark antiquark pair dū → e −ν e uū, with two powers of the strong coupling constant and two powers of the weak one. The input consists of the following arguments The algorithm proceeds as follows: • The first particle is kept fixed. The second particle is charge reversed, so that the process looks like the decay of the first particle into the remaining ones. At this stage we then have • We look through all (ordered) pairs of particles, excluding the first one, that have flavres equal to zero, and that can be merged into a single particle via a strong or weak interaction vertex. In the example at hand, we would find several cases: the second and last entry (a u and aū) merged into a gluon, a photon or a Z; for the third and fourth entry (an electron and its anti-neutrino) merged into a W − ; the last two entries (a u andū) merged into a gluon, a photon or a Z.
• For each found possible merging, we prepare a new input for the recursive procedure, with a new flavour list including the merged particle and updated values of the resonance pointers and of the power of the couplings. In our example, after the e −ν pair is merged into a W − , the new input for the recursive procedure looks like this Notice that now the flavres third and fourth entries (the e − andν) contain pointers to their mother resonance (the W − ), added in the seventh position. The value of powew has been updated, since one electroweak coupling was used for the W − splitting, and only one is left. Notice also that there are cases where the same particles can merge into a different one, as for the u andū merging case, and all these new inputs are passed to the recursive resonance-searching algorithm.
At this point three conditions are checked: whether no pairs can be further merged, whether no more powers of the coupling constants are available, and whether the last added particle coincides with the first one, meaning that all outgoing particles have been merged into the incoming one. If any of these conditions are not met, the configuration is abandoned.
The list just found represents a tree diagram for the process at hand. As such, there is always a unique path in the tree that joins any two external particles. The path joining the two incoming particle is the t-channel one. It can be found starting from particle 2 and going recursively through its ancestors, until particle 1 is reached.
The list is processed further, by the subroutine clean_resonance_structure, that performs the following operations. It first examines the t-channel structure of the flavour list. If it finds a t-channel fermion line that emits two electroweak bosons (W , Z or H) that can directly couple to each other trough a triple-boson vertex, with any intermediate emission of photons or gluons, it abandons the configuration. This is because another configuration with a richer resonance structure must exist, i.e. the structure where the two electroweak bosons arise from the decay of a single electroweak boson, with splittings involving the trilinear vector coupling, or the Higgs coupling to a vector boson. This richer configuration is well-suited to represent the one where the two electroweak bosons do not originate from another electroweak boson, and thus the latter configurations need not be considered. It then carries out a similar operation on s-channel lines. If we find a fermion line that emits two electroweak bosons, there must be a richer configuration where the two bosons originate from a single electroweak boson, and we thus abandon this configuration. Care is taken to handle the special case when the fermion line becomes a top quark, since the top is treated as a resonance, and the emission of a Higgs from the fermion line, since it must come from a top, in this case.
After the elimination procedure is carried out, all t-channel resonance entries, and all s-channel resonance entries corresponding to a massless particle are deleted from the list. The list is then put in a standard form: the resonances are moved just after the two incoming particles, and the final-state particles follow. The clean_resonance_structure exits. If the examined flavour structure is to be kept, the program calls a subroutine that stores it in a temporary array structure, provided there are no other equivalent configurations already stored. Once all configurations are found and stored, the subroutine pwhg_res_histos_born or pwhg_res_histos_real is called, and the configurations are transferred from the temporary storage to the global arrays flst_born* or flst_real*, that are overwritten with the Born and real flavour structure including resonance-history information.
The procedure that we have illustrated so far should be appropriate for most Standard Model processes. We checked that it works also in the case of single-top production studied in Ref. [52], by replacing the hand-written resonance histories that we used there with those automatically generated with the procedure presented here.

A.2 Colour assignment
In the POWHEG-BOX, colour assignment is mainly performed at the level of the underlying Born process. Given a Born flavour configuration and kinematics, one considers the colour subamplitudes that contribute to the squared Born amplitude, computed in the large colour limit. A color flow is then chosen with a probability proportional to the values of the subamplitudes, for that particular phase-space point. The POWHEG-BOX then generates the QCD radiation for a particular collinear region, through the splitting process. The colour configuration for the generated real-emission amplitude is obtained from the one of the underlying Born by attaching the colour flow that corresponds to that splitting.
Contributions that do not have singular regions (i.e. regular contributions) are instead treated as the Born term itself. In Ref. [79] a corresponding interface was developed such that this colour information could be extracted from MadGraph4. In Section A.3, we describe an analogous implementation in OpenLoops. Figure 17. Different resonance histories (top) for identical processes, and associated colour flows (bottom).
When resonance histories are being considered, a modification of this scheme becomes necessary. In fact, the randomly generated colour assignment may not be compatible with the resonant structure being considered. Consider for example the process qq → qqV V , where V is a vector boson, illustrated in Fig. 17. If the process proceeds via the exchange of a gluon in the s-channel, then the corresponding colour flow is illustrated in the bottom-left diagram in the figure, and initial-and final-state quarks are colour connected. Whereas if the exchanged particle in the s-channel is a colourless vector boson, then the color flow is depicted in the bottom-right diagram, where it is evident that initial-and final-state quarks are not colour connected, but there is a colour connection between the quarks in the initial state, and another colour connection between the quarks in the final state.
In the POWHEG-BOX framework, colour assignment is independent of the resonance structure, and thus one may end up assigning the colour flow in the left of the figure to the resonance history on the right, or vice-versa. In order to remedy to this, we keep gener-ating random colour configurations, and accept the first one that is compatible with the resonance history.

A.3 OpenLoops interface and settings
The OpenLoops program is based on a fast numerical recursion for the generation of tree and one-loop scattering amplitudes [80]. Combined with the OPP reduction method [81] implemented in CutTools [62] and the scalar one-loop library OneLOop [63], or with the tensor integral reduction methods [82][83][84] implemented in COLLIER [61], the employed recursion permits to achieve very high CPU performance and a high degree of numerical stability. The small fraction of numerically unstable one-loop matrix elements is automatically detected and rescued through re-evaluation with CutTools in quadruple precision.
The new POWHEG-BOX+OpenLoops interface is implemented via a Fortran90 module called openloops powheg, which is included in the POWHEG-BOX-RES framework. Internally the POWHEG-BOX+OpenLoops framework automatically compiles, loads and manages all required OpenLoops amplitude libraries. The new interface provides the subroutines openloops born, openloops real, and openloops virtual with interfaces identical to the corresponding POWHEG-BOX routines setborn, setreal, and setvirtual. In particular the openloops born routine returns, besides the squared tree-level amplitude B, the corresponding colour-and spin-correlated tree-level amplitudes B ij and B µν in the format required by the POWHEG-BOX [54]. Additionally, the interface provides the routines openloops init, openloops borncolour and openloops realcolour. The former synchronizes all parameters between OpenLoops and the POWHEG-BOX and should be called at the end of the init processes subroutine of the POWHEG-BOX. The latter two provide the required colour information as outlined in Section A.2, i.e. they return a colour-flow of the squared Born and real matrix elements in the large colour limit, on a probabilistic basis. Since the probability priors are determined from the colour-flow decomposition of the corresponding matrix elements at a given phase-space point, the colour-trace basis employed internally in OpenLoops is converted into a color-flow basis.
Several OpenLoops internal options and switches can be passed directly from the powheg.input file to the code. In particular, OpenLoops offers the possibility to switch between the tensor-integral reduction methods implemented in COLLIER and OPP reduction methods implemented in CutTools. By default COLLIER is used, while inserting the line olpreset 1 in the powheg.input file, reduction via CutTools can be selected. In a similar way, inserting the line olverbose <OpenLoops verbosity level> allows to select the verbosity level of OpenLoops.
While all relevant input parameters are automatically passed by POWHEG-BOX to OpenLoops, further internal OpenLoops parameters can be set directly via the routine (member of the Here, TYPE can either be integer, double or character(*) according to the parameter to be set, as detailed on http://openloops.hepforge.org/parameters.html.

Implementation of new processes
In order to set up a new processes within the POWHEG-BOX+OpenLoops framework, one should run the script <POWHEG-BOX-RES>/COMMON/OpenLoopsStuff/generate process.py with the following arguments ./generate_process.py <library name> -order_ew=<m> -order_qcd=<n> -name=<..> Here, <library name> corresponds to the OpenLoops amplitude library of the desired process and <m> and <n> denote the order of the Born cross section, O(α n S α m EM ), in terms of powers of the strong and weak couplings. This will setup a rudimentary POWHEG-BOX process structure within the directory <POWHEG-BOX-RES>/COMMON/<name>. For example the call ./generate_process.py pplnjjj -order_ew=2 -order_qcd=3 -name=Wjjj will yield the structure for an NLO+PS generator including all required tree and one-loop amplitudes for pp → W (→ ν) + 3 jet production.
A user has only to provide the list of contributing flavour structures of the Born and real subprocesses in the init processes.f file and the number of intermediate resonances in the nlegborn.h file. An automatic generation of these structures is currently being validated and will soon be included. Currently NLO QCD corrections to any SM process are supported by this interface, while NLO electroweak corrections will follow in the future.

Fixed-order NLO calculations
If one is interested in fixed-order NLO results, the most CPU-demanding contributions come from the computation of the real graphs, that also implement the cancellation of the collinear and soft singularities. In the POWHEG-BOX-RES code there are options to separate the virtual contribution from the rest, in such a way that it can be computed with an accuracy that matches the one of the real contribution, thus saving computer time. More specifically, the code can be run twice: in the first run, the user can set the flag novirtual to 1 in the powheg.input file. In this way, no call to the calculation of the virtual corrections is done, and the corresponding distributions do not contain the virtual corrections (plus other soft contributions). The code is then rerun by using the same importance sampling grids used in the first run, with the flag virtonly set to 1 and with lower statistics with respect to the previous run. In this way, the virtual contributions are called fewer times with respect to the Born and real contributions of the first run. Finally, the kinematic distributions obtained in the two steps can then be combined. Details and examples for this procedure are included in the release of the code.

Generation of Monte Carlo events
If one is interested in generating Monte Carlo events, it is more convenient to avoid the computation of the virtual corrections for the large number of events that are vetoed during the generation. This can be done, provided one renounces to generating events with constant weight. In essence, we generate events with settings such that the virtual contribution is not computed, but the cross section and the distributions are sufficiently similar to the exact result. The events are then reweighted with the full cross section including the virtual contribution. With this procedure, the virtual contribution is computed only once for each generated event, instead of the several tens of event that are computed and then vetoed in a standard run. In order to do so, one inserts in the powheg.input file the lines for_reweighting 1 rwl_file '-' <initrwgt> <weight id='xx'> some reweight info </weight> ... </initrwgt> and the program generates events with uniform weight with no virtual corrections. For each <weight line, a new weight is generated and added to the event. These weights are all computed with the inclusion of the virtual corrections. 14 We would like to remark an additional technical issue: the subtraction term in case of a massive fermion emitter (i.e. the subtraction term corresponding to the soft singularity in the b → bg splitting) is modified in such a way that it becomes closer to the P qg (z) Altarelli-Parisi splitting function (that is to say, we give it a weight (1 + z 2 )/(1 − z) rather than the original weight 2/(1 − z)). In fact, we found that if we do not include this modification, the kinematic distributions before reweighting (i.e. those with no soft-virtual contributions) can develop relatively large, negative values near the top-mass peak. After reweighting we do get back the correct results. However, reweighting coefficients can be very large or negative with the original subtraction term, while, with the modified one, we get sensible distributions even before reweighting, and no large reweighting factors.

B.1 Kinematic guess of resonance structures
In this appendix we detail the kinematic procedure for the construction of resonance information from agnostic events, i.e. events where no resonance information is available.
We start at the Les Houches event level and modify each event as follows: the two W 's come from the intermediate Z boson and the gluon is associated to initial-or final-state radiation (from the b's); the W 's and b's come from a tt pair, and the gluon from initial-state radiation; same as before but with the gluon from the top-resonance decay; same as before but with the gluon from the antitop-resonance decay.
If the colour of the b is not consistent with the colour assigned to the gluon, the corresponding weight is set to zero. Then, a resonance history and gluon assignment are chosen among the four configurations considered here with probability proportional to the corresponding w weight.
In order to validate this procedure, we stripped any resonance information from the Les Houches events of a resonance-aware res-singlerad sample, i.e. switching off the multipleradiation scheme. We indicate the corresponding results with the label res-strip. Then, following the procedure outlined above, we added back guessed resonance information. The obtained result, labelled as res-strip-guess, is displayed in Fig. 18 Figure 18. Comparison of the effect of removing resonance information and then adding it back according to a guess, based upon kinematics, in the invariant mass of the W j B system. procedure for the kinematic construction of resonance information reproduces nicely the correct W j B peak.

B.2 Comparisons of shower veto schemes
When generating radiation from resonance decay, the traditional Les Houches generic-user process interface is no longer viable. In fact, the standard [85] contemplates only a single scale (called scalup), and it requires that the shower does not generate any radiation harder than that scale at the production stage. Radiation in resonance decays remains, however, unrestricted, while our generator requires it to be vetoed either by the transverse momentum of the radiation generated by POWHEG in the decay process (allrad 1 case), or by the hardest radiation scale, irrespective of its origin (i.e. either from production or from resonance decay) in the allrad 0 case.
The default method for interfacing the bb4l generator to Pythia8 was taken from Ref. [35], and it is described in appendix A of that paper. In essence, the procedure was to examine the showered event, compute the transverse momentum of Pythia8 radiation in top decays, and veto it if higher than the corresponding POWHEG one.
Pythia8 provides with its own mechanism for vetoing radiation from resonance decay. One should implement a virtual function canSetResonanceScale that returns a true value if Pythia8 is to use this mechanism. Furthermore, one should also implement a function scaleResonance that Pythia8 invokes in each event for each resonance, returning the scale for vetoing radiation in decay. We also implement this mechanism in our generator. It is activated by setting the flag pythiaveto 1 in the powheg.input file.
We show in Fig. 19 the comparison of results obtained with the two veto mechanisms.   Figure 19. Comparison of two veto schemes on the B fragmentation function, on p B T,dec , on the mass of the W j B system and on the mass of j B distributions.
In these plots, as well as in all the others that we have examined, we have found very good agreement between the two veto schemes. We notice that the difference in the ratio of the m j B distribution at small masses (one of the few distributions where we found mild discrepancies) is taking place in a region where the cross section is getting small, and is thus of little relevance.
We conclude that the internal Pythia8 method for vetoing resonance radiation in decay is suitable for use with the bb4l generator, and we can thus recommend its use.

B.3 Impact of the multiple radiation scheme
Difference in tt observables induced by the multiple-radiation scheme of Eq. (2.6) were already discussed at length in Ref. [35] for the ttb NLO dec generator. It was found there that, by switching off the multiple-radiation scheme (allrad 0), radiation from top decays are mostly handled by the shower generator. In fact, the absolute hardest radiation is more often produced by the initial state, in part because of the larger colour charge, and in part due to the wider phase space available. Here we present some comparisons as a brief reminder of the relevant issues. In this section we apply our default cuts defined in Eqs. (5.9)-(5.10). We begin by showing in Fig. 20 Figure 20. Invariant mass of the W j B (left) and of the lj B (right) systems. We compare NLO+PS resonance-aware predictions with (res-default) and without (res-singlerad) employing the multiple radiation scheme. In the ratio plot we illustrate relative deviations with respect to res-default.
W j B and of the lj B systems. There is a good agreement between the two distributions, except for the region of low top virtuality in the left plot. On the other hand, observables that are sensitive to the B and j B properties display larger differences, as can be seen in Fig. 21.
In view of the large differences in the fragmentation function and p B T,dec distribution, we compare in Fig. 22 the same quantities computed using as reference frame the top quark at the level of Monte Carlo truth ("MC truth", usually identified with the last top quark appearing in the shower output list) rather than the reconstructed top. We also add to this comparison the output of the hvq generator. In this last case, we switch off Pythia8 matrix element corrections (MEC), for the purpose of determining whether the use of our generator, even if the allrad feature is switched off, brings out some improvement with respect to a generic shower treatment of top decays. We see from the figures that by using the MC truth for the top reference frame brings the bb4l and res-singlerad results in better agreement, at least as far as the p B T,dec distribution is concerned. The comparison of the bb4l, res-singlerad and hvq results for the p B T,dec distribution is particularly enlightening. If we focus upon radiation in the top decay, in the bb4l case the hardest radiation is always generated by POWHEG. In the res-singlerad case, POWHEG is mostly responsible for radiation with a large value of the p B T,dec observable, since it must be harder than the radiation generated in production. The region of small p B T,dec is thus more often determined by the shower. In the hvq case, radiation in the top decay is handled only     Figure 22. Predictions for the B fragmentation function and transverse momentum distribution of the B hadron in the top decay frame, p B T,dec obtained with the bb4l generator in its default mode employing the multiple radiation scheme (res-default), without employing this scheme (ressinglerad) and corresponding predictions obtained with the hvq generator (tt). In these predictions the top reference frame is determined according to the Monte Carlo truth (MC truth) and the tt predictions are obtained switching off matrix-elements corrections in Pythia8.
by the shower, that has only leading logarithmic accuracy, and thus fails at large values of the p B T,dec observable. This is why we see a large discrepancy between the hvq and the res-singlerad at large values of the p B T,dec observable.