On the Treatment of Resonances in Next-to-Leading Order Calculations Matched to a Parton Shower

In this work we present a new subtraction method for next-to-leading order calculations that is particularly convenient even when narrow resonances are present. The method is particularly suitable for the implementation of next-to-leading order calculations matched to parton shower generators. It allows at the same time for the inclusion of all finite width effects, including interferences, and for a consistent treatment of resonances in the shower approach, preserving the mass of resonances near their peak. We implement our method, in a fully general and automatic way, within the POWHEG BOX framework, and illustrate it using as a test case the process of $p p \to \mu^+ \nu_\mu j_b j$, that is dominated by $t$-channel single top production.


Contents
1 The problem At present, several methods exist for the computation of next-to-leading order (NLO) corrections in the Standard Model. When strong and/or electromagnetic interactions are present, these calculations must deal properly with collinear and soft divergences, that must cancel when infrared insensitive (IR-safe) observables are computed. The so-called subtraction methods are generally used in order to deal with this problem. In essence, they work as follows. A generic NLO cross section can be written symbolically as where Φ B stands for the Born phase space and Φ R is the real emission phase space. B(Φ B ), V (Φ B ) and R(Φ R ) represent the Born, Virtual and Real cross section respectively. The real emission process corresponds to the Born process in association with an extra parton. 1 For the purpose of this example we assume that we do not have hadrons in the initial state, i.e. we consider processes like Z decays into hadrons. The expression in eq. (1.1), in order to make sense, must be evaluated with some form of regularization for the soft and collinear singularities. Assuming that we are using dimensional regularization, the phase space Φ B and Φ R are evaluated in d = 4 − 2 dimensions. The value of an observable O is then given by We can think of our observable O as the cross section in a given histogram bin of some kinematic distribution. Again, in eq. (1.2) we assume that we have a d-dimensional definition for our observable, with the appropriate 4-dimensional limit. If the observable is IR-safe, soft and collinear singularities will cancel in eq. (1.2), yielding a finite result.

Subtraction method
In the subtraction method, one introduces a parametrization of the real phase space of the where Φ rad has d − 1 dimensions, and parametrizes the emission of the extra parton. The parametrization must have a smooth behaviour in the soft and collinear limit. Thus, in the limit of soft emission, the kinematics of all but the soft parton described by Φ R (Φ B , Φ rad ) must match the Φ B kinematics. In the collinear limit, the kinematics of the system obtained by replacing the two collinear partons in Φ R (Φ B , Φ rad ) with a single parton with the appropriate flavour, having momentum equal to the sum of the momenta of the collinear partons, must match the Φ B kinematics. One also introduces a simplified approximation to the real cross section, R s , that coincides with R in the soft and collinear singular limits. Eq. (1.2) is rewritten as eq. (1.5) is clearly identical to eq. (1.2). It has however the nice property that 1/ divergences in the square bracket of the first term on the right-hand side of the equation (arising in the virtual term and in the dΦ rad integration of R s ) cancel among each other. Furthermore, collinear and soft divergences cancel under the integral sign in the square bracket of the second term. The second term can thus be evaluated in 4 dimensions with numerical methods. The square bracket in the first term can be evaluated analytically once and for all. One usually defines (1.5) and the first term becomes the 4-dimensional expression (1.6) that can be computed numerically. The development of the subtraction method started since the very early QCD computations, already appearing in the bud in the calculation of the Drell-Yan process of ref. [1]. A more systematic use of it was made in ref. [2], in the context of e + e − annihilation into hadrons. In ref. [3] the calculation of ref. [2] was implemented as a parton level generator, such that any given observable could be computed with it without any dedicated analytic work, and was in fact used to compute a number of commonly used IR-safe observables for QCD studies at LEP. Subsequently, the subtraction method implemented in parton level generators was applied also for processes initiated by hadrons [4], and it became common practice to compute the R s term by using the collinear and the soft approximations in d dimensions (see for example [5]).
More recently, fully general formulations of the subtraction method have appeared. The procedure of ref. [6], known as the CS method, uses local subtraction terms for the R s cross section. The formulation given in ref. [7], known as the FKS method, is instead based upon the more traditional phase space parametrizations used in refs. [2] and [4].

The subtraction method and resonances
When resonances are present, in the zero width limit, the cross section factorizes into the product of production and decay terms. In these cases, a standard subtraction method can be applied independently to the production and decay processes. In fact this was done in refs. [8] and [9] for top production and decay. Problems do arise, however, if finite width effects are fully included, so that also interference among radiation produced in production and decays, or among radiation produced in the decay of different resonances, is included.
On the one hand, the presence of a finite width regulates the singularity associated with the resonance peak, so that, strictly speaking, a subtraction method will formally lead to finite and consistent results. On the other hand, taking the zero width limit, a standard subtraction method approach will lead to divergent results. In order to illustrate this problem, we consider the example of t-channel single top production and decay. One Born amplitude for this process is illustrated in fig. 1. The final state is composed by a b and d quark, a muon neutrino and an anti-muon. We assume for the sake of illustration a massless b quark. The system comprising the final state b quark, the neutrino and the anti-muon have an invariant mass close to the top mass, becoming identical to the top mass in the zero width limit. Consider now the real contribution obtained by adding gluon radiation to the final state. As illustrated above, in generic subtraction methods, the subtraction counterterms are obtained by factorizing the real phase space in terms of a Born phase space times a radiation phase space. The subtraction term for the collinear singularity corresponding to the final state gluon being collinear to the final b uses a Born phase space where the collinear bg pair is merged into a single b. The problem with the resonance is better illustrated in the CS subtraction framework, where the kinematics of the subtraction term is built as follows. Calling k ⊕ the 4-momentum of the incoming b, and k b , k g the 4-momenta of the final b and g partons, one defines the momentum of the b quark in the underlying Born configuration as in such a way thatk 2 b = 0. Furthermore, the incoming b quark momentum is redefined as so that the total 4-momentum is conserved. We see that in this way the 4-momentum of the top quark has been altered, near the collinear limit, by an amount Since the CS procedure does not impose that the top 4-momentum is the same in the real and subtraction terms 2 it will turn out that the top virtualities will differ there by an amount of order m 2 bg /E bg . The collinear singularity in the real and subtraction terms will thus match only if where Γ t is the top width. It is easy to see that this is also true in other subtraction methods. For example, in the one used in the POWHEG BOX, the momentum of the collinear counterterm is built by setting the 3-momentum of the b quark parallel to the sum of the 2 It instead imposes that the incoming b momentum minus the momenta of the final b and of the radiated gluon g in the real term equals to the b incoming momentum minus the final b momentum in the subtraction term, and all other momenta remain the same.
3-momenta of the b and g particles in the partonic CM frame. Furthermore, the momentum of the dμν system is boosted along the direction of the merged b quark in order to conserve 3-momentum, and the absolute value of the b quark momentum is chosen in such a way that the final state CM energy is conserved. This procedure is designed to conserve the mass of the final state system, and the mass of the system that recoils with respect to the splitting partons, i.e. theμνd system, while the mass of the top resonance is not conserved. We thus expect that the collinear singularities present in the real and subtraction terms will be exposed in the narrow width limit, spoiling the convergence of the subtraction method. In fact, double logs of the resonance width will arise in different regions of the real cross section, yielding to a failure of convergence in the limit Γ → 0. It is also clear that, in order to overcome this problem, one must devise a subtraction method such that the resonance mass is the same in the real and subtraction terms when approaching the resonance peak even when the resonance is off-shell by an amount greater than its width.

NLO+PS and resonances
If we plan to use an NLO calculation with an interface to a shower generator (NLO+PS from now on), further problems arise due to the resonance treatment.
In the MC@NLO method [10], one should consider the recoil scheme used by the Shower Monte Carlo to build radiation from a decaying resonance and construct the MC counterterms accordingly.
In the POWHEG method [11][12][13], one first computes the inclusive cross section for the production of an event with a given underlying Born configuration. Radiation is then generated according to a Sudakov form factor with the following form: The mapping of the real phase space into a product of an underlying Born times a radiation phase space is the same used in the NLO subtraction procedure. In general, it will not preserve resonance masses, so that in the R/B ratio, unless the condition (1.10) is met, the numerator and the denominator will not be on the resonance peak at the same time.
In case when R is on peak and B is not, this will yield large ratios that badly violate the collinear approximation. A further problem arises when interfacing the NLO+PS calculation to a Shower generator, in order to generate the next-to-hardest radiation. Shower Monte Carlo's should be instructed to preserve the mass of the resonances. Thus, radiation should have a resonance assignment. This is generally not available in processes that include interference among radiation generated by different resonances, or by a resonance and the production process itself.

The method
In order to solve the problems mentioned above, we should separate all contributions to the cross section into terms with definite resonance structure. Each term individually should have resonance peaks only in a single, well defined, resonance cascade chain. The mapping into an underlying Born configuration should be defined for these terms in such a way that the resonance masses are preserved. Thus, when looking for parton pairs that can give rise to a collinear singularity, one should only consider pairs arising from the same resonance decay, or directly from the production process.
In the POWHEG BOX framework, a subtraction method that preserves the resonance masses is already implemented, but it is presently available only for calculations performed in the zero width approximation. In these cases, only one resonance decay chain is possible, and the real emission contributions are separated according to the resonance that originates the radiation. The method is discussed in detail in ref. [14]. In essence, with this method, the subtraction procedure for initial state radiation is the same one used in ref. [12] (the FNO paper from now on). For final state radiation arising from the production process the subtraction procedure is also the same one discussed in Section 5.2 of FNO. This procedure is such that the mapping of the real to the underlying Born configuration does not change the four momentum of the final state. In case of radiation from the decay of a resonance, the subtraction procedure is essentially the same, except that it is applied in the resonance frame, and thus does not alter the resonance four momentum and the momenta of all particles that do not have the resonance as an ancestor.
In the general context when finite width effects are to be included, more than one resonance cascade chain (from now on "resonance history") may be present, and interference between amplitudes with different resonance histories must also be included. We thus need to perform a separation of the cross section into a sum of contributions, each one of them dominating only for a single resonance history. For each of these contributions we should apply the resonance aware subtraction method of ref. [14].
In the following we will describe in great detail the procedure adopted for the separation of the cross section contributions into terms with a definite resonance structure. We will discuss the procedure for the terms that have the Born kinematics (i.e. the Born, Virtual and Collinear Remnant terms), and for real terms. For the latter, subtraction terms having the Born kinematics are also present. We will require that in the collinear and soft limits the separation of the contributions associated with given resonance histories in the real term smoothly matches the corresponding separation in the Born kinematics.

The Born resonance histories
We need to single out contributions from the Born term corresponding to several different resonance histories of the final state. Each resonance history corresponds to a tree graph, where the leaves of the tree are the final state particle, and the intermediate nodes are the resonances. In our case, we include in the tree also the two initial state particles, before the root node.
The root of the tree does not correspond necessarily to any real resonance. For uniformity of treatment, we will however associate to the root a fictitious resonance, and we will refer to it as the "production resonance".
For each given initial and final flavour configurations, we have several possible resonance histories. We will denote with F b the initial and final flavour structure of the Born process, irrespective of the internal nodes of the resonance history. We will instead denote with f b the flavour structure including the resonances decay cascade. We will also refer to it as the resonance history. Summarizing, we will refer to F b as the bare flavour structure of the process, and to f b as the full flavour structure, or simply as the flavour structure.
The Born contributions will be labeled as B F b . Thus, B F b is the square of the amplitude for the production of the final state F b , including all possible resonance histories allowed for the process. We separate the Born contribution in the following way: where T (F ) is the set of all trees having the same bare flavour structure F . The factors Π f b have the property Furthermore, they must be such that Π f b B F b must have resonance peaks compatible with the resonance history of f b . One possible definition for the Π f b is the following. With each resonance i in the resonance history, we associate the factor and define where s i , M i and Γ i are respectively the invariant mass of the decay product system, the mass of the resonance and its width. By Nd(f b ) we denote the set of all nodes of the resonance tree for f b (excluding the root). We then define where we have introduced the notation F b (f b ) to denote the bare flavour structure associated to a given full flavour structure f b . This definition clearly satisfies the property (2.2). Thus B f b exhibits resonance peaks only in correspondence with resonances in its own resonance history. In fact the P f b factors for all alternative resonance histories in the denominator of Π f b cancel the resonance peaks due to alternative resonance histories in Only the peaks compatible with the f b resonance structure, that have a corresponding enhancement factor in the numerator, will remain. It is worth pointing out that our definition of the Π factor is certainly not unique. In particular, there is an alternative possibility that is easily implemented if one has access to the individual sub-amplitudes contributing to the total amplitude characterized by F b : (2.6) The structure of each sub-amplitude represents in this case a resonance history, so that we can create a correspondence i ↔ f b , and define This possibility may prove convenient with current numerical matrix elements programs, where the numerical calculation of the individual amplitude is a necessary step for the computation of the full matrix element. Since this procedure is gauge dependent, care should be taken in the choice of an appropriate gauge.

Implementation of the Born resonance histories in the POWHEG BOX
The internal implementation of the Born flavour structure can be inherited from the present Born level structure in the POWHEG-BOX-V2, starting with the extension of ref. [14] for the inclusion of narrow width resonances. In this implementation, the full flavour structure of a Born term is represented by two arrays, flst born(j,iborn) and flst bornres(j,iborn), where the index iborn labels the particular Born full flavour structure f b . The j index labels the external leg and the internal resonances, with 1 and 2 representing the incoming legs, and the (integer) value of the flst born array represents the corresponding flavour code (that coincides with the PDG code, except for gluons, that are labeled 0). The flst bornres(j,iborn) integer array represents the resonance pointers, so that the whole resonance structure can be reconstructed. For example, for the case of the full flavour structure corresponding to the process gg We see that the resonance pointer list contains zero for particles generated at the production stage (in POWHEG we represent the fictitious production stage resonance as having index 0), while for particles produced in resonance decays the corresponding entry is the position of mother resonance in the list. At variance with the V2 implementation, in the present case we must be prepared to assume that not all Born flavour configurations have the same resonance history and the same number of resonances, so we must admit Born flavour lists of different length. We thus introduce an array flst bornlength(iborn), carrying the length of the flavour list for the Born f b labeled by the iborn index. The entries of this array are set in the user process init processes subroutine.
The POWHEG BOX integration program (the mint integrator [15]) was updated in order to deal with the resonance histories. Since several resonance histories may be present, the mint integrator was also updated to be able to deal with a discrete (summation) variable. It now computes a multidimensional integral in a unit hypercube and the summation over a discrete index. The discrete index is used to label each resonance history. The phase space generator examines the value of this discrete index, identifies the corresponding resonance history, and chooses automatically a phase space parametrization that performs importance sampling over the resonance regions, generating the resonance virtualities with an appropriate Breit-Wigner distribution.

The real resonance histories and singular regions
In the case of real graphs, we have more resonance histories, because we have one more final state particle that can belong to resonances. In analogy with the Born case, we introduce F r and f r as before, labeling the bare and full flavour structure for a real graph. We will now introduce a label α r , that labels a singular region 3 compatible (in a sense that we will specify in the following) with a given f r α r ∈ Sr(f r ). (2.8) Also here we will use the notation f r (α r ) and F r (α r ) to denote the full and bare flavour structure associated with a given singular region. We only consider singular regions that are compatible with the given resonance history in the following sense: the particles that become collinear should be siblings, i.e. should arise directly from the decay of the same resonance or from the root (if they are directly produced in the hard reaction).
We now perform the separation of the real cross section with a given bare flavour structure into singular region contributions: where F r (α r ) stands for the bare flavour structure associated with α r . We require that the real weights P fr(αr) are compatible with the Born weights, in the sense that, in the soft or collinear limit, the P fr(αr) must approach smoothly a P f b factor of the corresponding underlying Born. This is certainly the case if they are defined as in the Born case. We notice that in the standard POWHEG scheme, the real contribution to a given region is enhanced if the collinear pair has a smaller transverse momentum than all other possible collinear pairings. In the present scheme, the relative transverse momentum of the pair is no longer the only element that decides about the partition of singular regions. As an example, consider three final state partons i, j and k. The cross section is parted among the i, j and i, k singular regions, depending upon how small are the relative transverse momenta in the two cases, and how far from the resonance peaks are the resonances containing respectively the i, j and i, k partons.
The d −1 factors used in the POWHEG BOX have the form where b is a positive constant parameter. Eq. (2.10) is used for the collinear region characterized by parton i, with energy E i and angle θ i (relative to the beam) in the partonic CM, becoming collinear to either incoming hadrons. Eq. (2.11) is again for initial state collinear regions, but distinguishes among the two collinear directions. Eq. (2.12) is used for the region characterised by final state partons i and j becoming collinear. They are commonly evaluated in the partonic rest frame. In the present case, however, in case of final state singularities associated with the decay products of a resonance, it seems more natural to compute them in the resonance rest frame. They thus become dependent upon the full flavour structure f r of the real contribution. It is however important for the following developments that the d ij factors do not depend upon the resonance structure in the collinear limit. This is in fact the case with our definition, since where lim ij denotes the limit for particles i and j becoming collinear and z is the energy fraction. The last expression is obviously Lorentz invariant in the collinear limit. Thanks to this property, it will turn out that the sum of all R αr associated with the same underlying Born full flavour structure factorizes in the collinear limit 4 This follows from the fact that all (and only) the α r associated with particles i, j becoming collinear dominate in this limit, and, being all equal, they simplify out in the numerator and denominator of eq. (2.9). We emphasize, however, that the d ij terms are not frame independent in the soft limit. This is quite clear from eq. (2.12), that in the E i → 0 limit becomes that is clearly frame dependent. As in the Born case, the scheme discussed here is not the only alternative for the partition of the singular regions and of the resonance structure. Using weights equal to the square of individual sub-amplitudes is still a valid alternative, as long as one computes the amplitudes in a physical gauge, in such a way that squared amplitudes also retain the full collinear singularity structure. In this case one does not need to introduce the d ij factors, since the squared amplitudes already have the appropriate singular behaviour in the collinear limit. In order to further pursue this alternative, issues related to the lack of gauge invariance of the individual amplitudes squared should be addressed. In the present work we did not investigate this alternative any further, since we prefer to assume that in general the individual amplitude for the process may not be available.

Example: electroweak uū → udūd
We illustrate the separation of the resonance structures in the process uū → udūd, considering only electroweak interactions. In order to simplify the discussion, we will (wrongfully!) assume that only the diagrams illustrated in fig. 2 contribute to it. We remark that this process is chosen only for illustration purposes. We are aware of the fact that it has no physical relevance and that we are omitting other relevant resonance histories. There is only one F b , corresponding to the bare flavour structure uū → udūd. We have two f b , represented in fig. 3, corresponding respectively to uū The P factors for the two configurations are Notice that we have assigned the values 1 and 2 to the f b index of the two flavour configurations depicted in the figure. Particles are labeled by an integer, starting from the lower incoming line, and going through all other particles clockwise. Summarizing, we have two (full) flavour structures for the given bare flavour structure uū → udūd. The corresponding Born contributions will be given by Notice that B is the full Born contribution, given by the square of the sum of the graphs in fig. 2. However, B 1 will be dominated by the square of the first graph, and B 2 by the second.
The number of real graphs is already quite large, and we do not show the corresponding figures. They are obtained by adding one final state gluon to the Born flavour configuration, and by replacing one of the initial lines with a gluon, adding a corresponding quark of opposite flavour to the final state. Here we focus upon the bare flavour configuration uū → udūdg. The corresponding full flavour configuration trees are depicted in fig. 4.
We will now label the gluon as 7, and keep the same labels used in the Born case for all other particles. The P factors are now The singular regions α r are displayed in tab. 1. Notice that the final state radiation d factors carry in the subscript the position of the two partons that become collinear, and, after a comma, an index specifying the resonance history. We are in fact assuming that the d factors are computed in the frame of the resonance that owns the two collinear partons. Notice also that the standard (non resonance aware) POWHEG implementation would have found 5 regions, one for the initial state radiation, and 4 for final state radiation, corresponding to a gluon being emitted by each final state parton. It is interesting to see how the singular part of the cross section is shared among the various resonance histories. We consider as an example the gluon emission from particle 3, carrying the d −1 37 singularity. In the standard POWHEG formulation this region corresponds to a single α r . On the other hand, in our resonance aware extension, that singularity is shared by the α r number 2 and 7. The one of the two that is more enhanced by resonant propagators (i.e. by its P factor) will dominate over the other. We have Notice that near the 3,7 collinear singularity, we have where the last equality follows from our requirement that the d factors are Lorentz invariant in the collinear limit. We thus see that in the collinear limit the collinear contribution is distributed among the 2 and 5 resonance histories, favouring the one that is nearer the resonance peaks. The underlying Born corresponding to the real flavour configurations f r ∈ {1, 2, 3} is the Born flavour configuration f b = 1. In the limit of vanishing momentum of the additional radiation in leg 7, i.e. when the momenta on the legs 1-6 in the real diagrams can be mapped to a given set of momenta of its underlying Born diagram, all the P r factors of the real flavour configuration reduce to the P b factors of the corresponding underlying Born contributions. In our case: This implies that in the same limit: (2.21) so that, in the soft limit, for example , (2.22) and a similar relation holds for all other α r contributions. As shown in eq. (2.14) we cannot drop the resonance history dependence in the d ij factors. Thus, unlike the case of collinear singularities, in the soft limit a full factorization of the P and d −1 factors does not hold in general. We will see in the following sections that this fact leads to a minor complication in the evaluation of the soft contribution.

Soft-collinear contributions
In the V2 implementation of narrow resonance decays, the soft collinear contributions (to be added to the virtual one) are computed assuming that no interference terms arise between the different resonances (or between a resonance and the direct production). In the finite width case we are considering now, this restriction has to be removed, because such interference terms do arise. Furthermore, the soft-collinear contributions depend upon the adopted subtraction procedure, and we are now departing from the default one used in the POWHEG BOX. We thus need to discuss in detail and compute the soft-collinear contributions in the present framework.
As specified previously, each singular region α r is associated with a single full flavour structure f r (α r ). The treatment of initial state singularities remain the same as in the standard case, since no resonance decays are involved in ISR (initial state radiation). We thus focus upon FSR (final state radiation). Thus, from now on, the singular region α r corresponds to final state particles i and j becoming collinear. Furthermore, the soft singularity is associated with particle i becoming soft. This is the case when i is a gluon and j is a quark. If also j is a gluon, in POWHEG, the R αr contribution is multiplied by a factor of the form where h is typically a power. This factor damps the soft singularity when particle j becomes soft. Since the cross section is symmetric in the exchange of the two gluons, this procedure leads to the correct result.
Given the singular region, POWHEG selects a phase space mapping from the Born phase space Φ B and the three radiation variables ξ i , y ij and φ to the full real emission phase space. In the α r region, partons i and j will arise from the same resonance k res . The phase space mapping will thus be chosen in such a way that only the momenta of the decay product of the resonance will be affected. For example, in case of partons i and j corresponding to a gluon and a b quark arising from top decay, the phase space mapping will build the radiation phase space starting from the momenta of the top decay products (i.e. the b and the W + ), maintaining fixed all remaining momenta together with the top four-momentum. The mapping procedure will correspond to the prescription described in the FNO paper (ref. [12]), applied to the top decay product in the rest frame of the top. It is the same procedure that is applied in ref. [14] for the case of tt production and decay in the factorized approach.
Following the notation of the FNO paper, we thus write the phase space as where d = 4 − 2 is the dimensionality of spacetime. Furthermore, we introduce the parametrization where ξ is defined as , computed in the rest frame of resonance k res . In the soft limit, the phase space becomes where, following the notation of FNO, the underlying Born kinematics is expressed in terms of the barred variables, and dΦ B is the underlying Born phase space. The α r contribution to the cross section can be written as We now introduce the expansion i.e. is defined as a distribution with a vanishing integral between 0 and 1. We get with (3.8) where the meaning of the second equation is simply to replace ξ −1−2 with ξ −1−2 + in the real cross section integral, since the corresponding δ(ξ) contribution has been subtracted out.

Soft terms
We now discuss explicitly the computation of the soft term I s,αr . In the standard treatment, by summing over all singular regions one recovers the full R, that can be approximated in the soft limit by the eikonal formula. We cannot follow this procedure now, since it requires that the soft limit is taken in the same frame for all α r , which is not our case. In order to deal with this complication, we employ the following trick. We use the identity to rewrite I s,αr as where withR αr we denote the Taylor expansion of R αr in the soft limit for k i : We notice thatR αr is obviously independent from the frame used to define ξ. It is in fact obtained from R αr by linearizing it in the k i momentum. In formula (3.11) the frame dependence of the soft contribution is all contained in the exponential, the rest of the expression being fully Lorentz invariant. In order to perform the integral we proceed as follows. We rewrite (3.11) as where m is an arbitrary timelike momentum. For definiteness, we choose m = q, the total four-momentum of the final state particles. The I s,αr term in eq. (3.13) is infrared finite. The I (2) s,αr term can now be integrated in any frame we like. We then just pick a common frame for all α r that have the same underlying Born bare flavour structure F b , and sum over all of them. We get (3.14) In the sum ofR αr , all dependencies upon the partition of the resonance regions and upon the d −1 factors cancel out, yielding the fullR for a soft gluon emission with an underlying Born flavour configuration equal to F b . In fact, the bare flavour structure of the α r such that F b (α r ) = F b consist of the same flavour assignment F b plus one (soft) gluon. Thus, from eq. (2.9), it follows that in the sum in eq. (3.14) all resonance history and d −1 factors simplify. We thus have Performing the integration in the rest frame of m, we can use again the replacement that leads to the standard calculation of the soft contributions, regardless of their resonance assignment. The corresponding result is reported in Appendix A.1 of the POWHEG BOX paper (ref. [13]). At this stage, we can split the Born terms again in terms of their full flavour structures, and thus compute each contribution using the appropriate (importance sampled according to the resonance structure) phase space.
The treatment of the I s,αr term of eq. (3.13) requires some care, since although the soft singularity is no longer there, it has still a collinear singularity corresponding to the α r region. We evaluate the I (1) s,αr integral in the CM rest frame. Other choices are possible, but there is no reason to make more complex choices, since the resonance virtualities inR αr do not depend upon the soft momentum k i , and thus no particular importance sampling is needed in the k i integration. We thus write where y = cos θ ij , (3.18) and j is the emitter associated with the α r region. By expanding we can write where The I (1) s+,αr term has to to be computed numerically. It has no analogue in the previous POWHEG BOX implementation.
We now work through the I sδ,αr . We have where C j is the Casimir invariant associated with the particle that underwent the splitting for the region α r . Observe that in deriving the identity (3.23) we have assumed that in the collinear limit the d −1 factors associated with a given pair of final state particles all coincide, irrespective of the resonance structure, as we have remarked earlier. Using we get where we have written, in the collinear limit andk j is the momentum of the emitter in the soft limit, i.e. at the underlying Born level. Performing the ξ integration (from 0 to ∞), we get where we have introduced the common normalization factor (3.28) The normalization factor in (3.28) should be the same one adopted in the computation of the virtual term, as defined in the FNO and POWHEG BOX papers. If this is the case, the 1/ singularities in the virtual and soft virtual contributions cancel exactly, and one needs only to retain the finite terms.

Collinear terms
We now turn to the collinear integral We are considering the region where i, j become collinear, and where, as discussed previously, k j does not lead to a soft singularity. This is the case for gq pairs i should correspond to g. In the gg case, the d −1 ij factors are supplemented with an energy damping factor where, as in ref. [12], h is typically defined to be a simple power law. Since there is no soft singularity in k j , in the following we can assume that we have a lower cutoff on the k j energy, that can be smoothly removed at the end of the calculation, so that in our manipulation we can always assume that k j is not vanishingly small. We now write the I +,αr term as where we remind that ξ is evaluated in the frame of the resonance to which i and j belong. We introduce for Φ i the phase space (always defined in the rest frame of the resonance) where the angular integration is done with respect to the j direction, i.e. y = 1 − cos θ ij . We separate out the collinear divergent term by using eq. (3.19), that yields I +,αr = I +δ,αr + I ++,αr , and Eq. (3.33) should be interpreted as following: take the dΦ n+1 R αr expression, work it out in y and ξ variables, and replace 1/(1 − y) and 1/ξ by the corresponding + distributions. We now write the Φ j integral in terms of angular and radial variables, and introduce a variable k 0 = k 0 Because of the δ(1 − y) factor we have now that k j and k i are proportional, and thus Ω j represents their common direction. Defining and defining performing the dk 0 j integration using the energy δ function we get k 0 j = zk 0 , and trading k 0 i for z we get We notice that the expression on the second line corresponds to the phase space of the parton into which i and j have merged, i.e. thek j phase space. The whole expression thus becomes where we have replaced k 0 withk 0 j , that is the energy of the underlying Born emitter. We have and where with P αr (z) we mean We thus get Let us begin by evaluating the integrals where in the first two steps we have used twice the z → 1 − z symmetry of the integral. For the q → qg case we have Finally, for the g → qq case: (3.50) We now define, as usual and find where by j(f b ) we mean the flavour of the j th parton in the f b flavour structure. We get where we have used for ξ max the covariant expression ξ max = 2k j · k res k 2 res . (3.54) We now expand it as We now combine this term with the I sδ,αr integral: where f b stands for f b (α r ), and j is the emitter for the region α r . Notice also that nowk 0 j represents the energy of the emitter in our common reference frame, while earlier, with the same symbol we denoted its energy in the resonance frame. Notice also that ξ max is now frame dependent. We will denote as I A the finite part of I A .

Summary
We now summarize the real and soft-collinear terms that need to be included in the calculation: A. Real integral: B. Collinear terms: where j is the emitter and k res is the momentum of the resonance that contains the emitter for the region α r . Withk 0 j we denote the energy of the emitter in our common reference frame, that is the CM frame of the final state. On the other hand, ξ max is computed in the resonance frame.
C. Soft terms remain the same as in the standard treatment, and are reported in Appendix A.1 of the POWHEG BOX paper (ref. [13]).
D. Soft mismatch: to be summed over all the α r with the emitter belonging to a resonance.
The integration phase space is defined in the partonic CM frame, with the third axis pointing along the direction of the emitter j.
E. Collinear terms related to initial state radiation remain the same as in the standard treatment.
In tab. 2, we list the terms that needed to be newly implemented in the new, resonance aware version of the POWHEG BOX that we are presenting here.

Soft log Γ terms
In the procedure that we have illustrated, collinear singular regions arise only among partons produced in the decay of the same resonance. This property arises because, in the separation of the singular regions, we restrict ourselves to singular structures that are compatible with the resonance history. While this feature guarantees a smooth cancellation of the collinear logarithms in the subtraction procedure, we cannot expect a corresponding cancellation of all soft, non collinear logarithms. There are in fact two sources of soft radiation with a lower or upper cut off of the order of the resonance virtualities: • Soft radiation arising from the interference of soft emissions from coloured partons belonging to different resonances. These terms have an upper cut-off of the order of the resonance width.
• Soft emission involving amplitudes with radiation arising from the resonances internal lines. These terms have a lower cut off of the order of the resonance width.
We thus expect that in our procedure log Γ terms will arise in the integration of the real cross section. The virtual corrections will also have corresponding log Γ terms, that cancel the real ones when summed together.
In this section we discuss the structure of these soft terms. As we will see, it is possible, in principle, to remove them from the integration of the real cross section, and include them in the soft term, in such a way that their cancellation takes place in the softvirtual contribution. However, we have not attempted to implement this in the POWHEG BOX. In view of the relatively large size of the resonance widths in the typical processes that we consider, it is unlikely that they may cause problems in practical NLO calculations. Furthermore, as far as NLO+PS implementations are concerned, these terms are in fact properly treated in our resonance-aware POWHEG framework, and do not require any further action.
We now discuss the structure of the soft logarithms in the presence of narrow resonances. For simplicity, we assume that all the resonances have comparable widths of order Γ 0 . We consider two regions: • Region a: is characterized by soft emissions with energy ω larger than Γ 0 . In this region Γ 0 plays the role of an infrared cutoff. The dominant region of integration has log ω uniformly distributed between log Γ 0 (lower cut-off) and the log of some hard scale in the process (high cut-off), typically of the order of the mass of the resonances.
• Region b: is characterized by soft emissions with energy ω less than Γ 0 . The lower limit in this region is regulated in the usual technical ways (like dimensional regularization). Its upper cut-off is Γ 0 .
In region a, since Γ 0 acts as an infrared cutoff, the emissions from resonances internal lines near their mass shell should also be considered as soft. In fig. 5 we illustrate the insertion of a soft emission in an internal resonance line. The product of the resonance propagators will be given by Under the assumption that ω = k 0 > Γ, the two denominators cannot be near their mass shell at the same time. When the first term in the square bracket is near its mass shell, the process corresponds to the resonance radiating during production. In fact, in this case the p momentum is far from the mass shell by a scale of order k 0 , while the p − k momentum is near the mass shell by a scale of order Γ. In coordinate space, this means that the line carrying momentum p has a length of order 1/k 0 , much shorter than the length of order 1/Γ of the p − k line. Conversely, if the second term is on-shell, radiation is taking place during decay. When squaring the amplitude, interference between these two terms is suppressed, since the two propagators cannot be on-shell at the same time, and the integration is effectively cut off by at a scale of order Γ, leaving no phase space for soft logarithms to build up. For the same reason, interference from emissions arising at production with emissions from resonance decay, as well as from emissions arising from the decay of different resonances, do not yield soft logarithms, since they also lead to propagators off the resonance peaks in the interfering amplitudes. Reasoning in terms of radiation and decay times, by assuming ω > Γ 0 we are assuming that radiation time is shorter than the resonances lifetimes. Thus, soft radiation in production cannot interfere with radiation in decays, since they happen at different times, and for the same reason radiation from different resonances cannot interfere. So, as far as soft singularities are concerned, the process can be though of as the product of independent production and decay processes, each one of them with resonances appearing only as initial or final state particles, but not as internal lines. For all these independent components, soft emissions is given by the usual eikonal formula applied only to initial and final state particles, that in this case can also be unstable resonances.
The structure of the soft singularity in region b is determined by the initial and final state particles after the decay of all resonances. The resonances are considered as off-shell particles, as far as soft emissions are concerned, and interference terms from emissions arising from different resonances are not suppressed by small Breit-Wigner weights, since the emission energy is below the resonance widths. In terms of time, this is the case when the time for soft radiation is longer than the resonance widths, so that only particles that live longer than the resonances can contribute.
The form of the soft subtraction term in the b region is the same one that we have adopted in the present method. However, that in our present treatment we are considering unrestricted emissions, while the b region is defined to involve soft energies below Γ 0 . When considering a given underlying Born resonance history, we will thus have the following cases: • In the emissions from pairs of coloured massless partons belonging to the same resonance, the terms in regions a and b will combine, yielding an unrestricted soft energy integral, and no log Γ terms.
• In the emission from pairs of coloured massless partons belonging to different resonances, only the terms in the b region will be present. These contributions will be cut-off at energies above Γ 0 , since for larger energies they will push one of the two resonances out of its mass shell, thus damping the cross section. They will thus lead to log Γ contributions to the cross section.
• In any emission from an internal resonance leg, the emission energy will have a lower infrared cutoff Γ 0 , and will yield other log Γ terms.
It is conceivable that our method may be modified, by adding further soft subtraction terms to the real cross section and corresponding integrated soft terms to the soft-virtual cross section, in such a way that the log Γ terms cancel within the soft-virtual contribution. This procedure may make the NLO calculation more convergent in the zero width limit. However, it would have no effect in the generation of radiation according to the POWHEG method. In POWHEG, the cancellation of the log Γ terms takes place numerically in the calculation of theB function, between the real and the soft-virtual integral. In the generation of radiation, the cross section is unitarized by construction, so that no further log Γ terms arise in inclusive quantities. We thus did not attempt to implement such an improvement in the present work.

Code organization
The implementation of the subtraction scheme described in the present paper has required an extensive rewriting of several parts of the POWHEG BOX framework. While we postpone writing a full documentation for the new code, we will describe in the present section the structures that are used to describe the various components of the cross section in terms of flavour and resonance histories.
The flavour structures used to implement our subtraction scheme are organized as follows. The process specific code provides the flavour structure in terms of arrays carrying the flavour of the particles involved in the process, including intermediate resonances. These have the purpose of grouping together the Born full flavour configurations that have similar resonance structure so that they can be integrated together. Thus, the value of flst bornresgroup(iborn) (an integer from 1 to flst nbornresgroup) labels the resonance structure group of the born full flavour structure iborn. This is needed because the POWHEG BOX groups together flavour structures with similar resonance histories when performing the integration, since these configuration can be integrated with the same importance sampling.
In the user process arrays, each flavour structure can appear only once, following the traditional approach of FNO. In the case of resonances, this leads to a non-trivial subtlety that we describe now. Consider the flavour structure for the subprocess qq → e + e − e + e − . According to the traditional approach of the POWHEG BOX there is only one flavour structure associated with this process, i.e. only a single ordering of the final state electron and positron will appear. When resonance structures are considered, we realize that we have two ways of pairing the electrons-anti-electrons to build an intermediate Z boson. These two pairings are fully equivalent up to permutations of the final state particles, and thus only one resonance assignment will appear for them. We should not forget, however, that the contribution with a given resonance assignment carries a resonance projection factor. Assigning the ordering 1 to 8 to the qq → ZZ → e + e − e + e − , assuming that the Z in position 3 decays into the e + e − pair in position 5-6, and assuming that we have only these two resonance histories, we will have a factor of the form P (5, 6; 7, 8) P (5, 6; 7, 8) + P (5, 8; 7, 6) , where, for example, It is clear now that we should supply a factor of 2 for this graph, since by assigning the resonances we break part of the exchange symmetry for final state identical particles. Thus, the user process should provide the symmetry factor appropriate for the given final state irrespective of resonance assignment. The POWHEG machinery will take care of supplying the appropriate factors arising from the resonance assignment specification. Thus, the user process list is built from the list of flavour structures for all distinct final states (where distinct means that final state differing by permutations are not allowed). For each final state, the list will be expanded by assigning all possible resonance histories. But, again, full flavour structure differing only by a permutation of the resonance histories will not be allowed. The POWHEG BOX checks explicitly that no full flavour structures equivalent up to a permutation will appear. Notice that the factor of two will lead to a total resonance factor of 2P (5, 6; 7, 8) P (5, 6; 7, 8) + P (5, 8; 7, 6) , that is not 1. However, by symmetry, an analysis of generated events that does not distinguish among identical final state particles will lead to the same results as if we included both weights For example, if we have several gluons and a quark in the final state, there will be regions associated with each gluon being collinear to the quark, and the program will find as many regions of this type as there are gluon. It will recognize that all these regions are equivalent, and it will emit a single region with a multiplicity factor equal to the number of equivalent regions.
If there are real contributions that do not have any singular region, they are collected into the "regular" arrays An array flst regularmult(1:flst nregular) is provided also in this case.
The task of finding out the singular and regular contributions is carried out by the subroutine genflavreglist, in the file find regions.f. In the traditional POWHEG BOX implementation, at this stage, for each alr a list of competing singular regions was found. These were needed since each alr contribution is obtained by multiplying the corresponding real graph by the ratio of the d −1 αr factor divided by the sum of all the d −1 αr associated with the other competing singular regions. In the current implementation, we also need to list together the competing resonance histories that lead to the same final state, in order to compute the corresponding resonance projection factor. So, for either the alr, the regular or the Born terms, we have an array of pointers representing the full flavour structure of the competing resonance histories for the Born graphs. We stress that we cannot use flst real and flst born in place of flst allrh and flst allbornrh, because in the latter also configurations that differ by a permutation of the intermediate resonances are included. The rh arrays described above are filled by the subroutine fill res histories, in the fill res histories.f file, that also computes the multiplicity factor associated with alternative resonance histories that differ only by a permutation of the resonances from the contribution being considered. These arrays are required for the computation of the weights needed to project out a given resonance history contribution from the real and Born amplitudes. Besides these, in the case of the alr contributions, we also need to know the singular regions associated with a given resonance history. This information is contained in the arrays flst_allrhnumreg(maxreshists), flst_allrhreg(1:2,maxregions,maxreshists).
For each entry of the realrh arrays, flst allrhnumreg gives the number of singular regions, and flst allrhreg gives the indices of the two particles becoming collinear in the corresponding singular region.

The example of single top, t-channel production
We consider the Born level process bq → be + ν e q . In the following we will label for conciseness as q and q all light quarks or antiquarks (excluding the b) with the appropriate flavour structure that can appear in the process, and we imply also the presence of the corresponding processes with exchanged initial state particles (i.e. qb in this case). This process is dominated by the single top production process bq → (t → b(W + → e + ν e ))q , such that the top quark is not produced at rest in the partonic centre-of-mass. Therefore this process is relevant for testing our formalism. In fact, the standard momentum mapping leading to the underlying Born configuration, in case of collinear radiation from the b quark arising from top decay, would conserve the incoming partons 4-momentum by adjusting the 3-momenta of the b and the W + q systems with appropriate boosts. This procedure would preserve the mass of the W + q system, but not the mass of the top.
At the Born level, it is enough to consider a single resonance history, namely bq → (t → b(W + → e + ν e ))q . Alternatively, one may consider two different resonance histories: The second one is actually not needed, since treating the bW + system as a resonance (i.e. preserving its mass in the underling Born mapping) rather than preserving the mass of the W + q system, as the POWHEG BOX would do for the resonance history of eq. (5.2), does not lead to any inaccuracy. We did however include this resonance assignment as an option, and used it to test that our setup works also in the case when more than one resonance history is present at the Born level. We are considering the following real processes: bq → e + ν e q g and bg → e + ν e qq . We do not consider processes of the form qg → be + ν e q b , that include also s-channel contributions. This is adequate for the purposes of the present paper, where we would like to present and validate a method, rather then provide a realistic simulation of single top production.
We will now list the resonance histories for the real contributions corresponding to the choice of a single resonance history at the Born level:

4)
bg → (t → bg(W + → e + ν e ))qq , Notice that the last two processes are really regular ones, since for them no collinear singularity can arise by pairing particles belonging to the same resonance.
The Born (together with the colour correlated Born) and real matrix elements for the process were easily generated using the MadGraph4-POWHEG interface described in ref. [16]. The virtual contribution was extracted by hand from code generated using MadGraph5_aMC@NLO [17].

Test at the NLO level
We first tested our method (that we refer to as POWHEG-BOX-RES) by comparing its NLO level results with a (traditional) POWHEG-BOX-V2 implementation of the same process. More specifically, we implemented the bq → be + ν +eq process in the POWHEG-BOX-V2 framework, without the inclusion of any resonance information, and using exactly the same matrix elements used with the new method. The comparison was carried out by using exactly the same choice of scales and parton density functions.
We observe that the process (5.6) has initial state collinear singularities, due to a contribution arising from an initial state g → bb splitting followed by the subprocess bb → (Z/γ → (W + → e + ν e )(W − → qq )). This process is not subtracted in our procedure, since we do not have a corresponding underlying Born contribution. In order to avoid the associated collinear divergence, and only in the framework of the NLO tests that we are discussing in this section, we supplied our cross section with a damping factor of the form that damps low transverse momentum b emissions. It is applied to all terms of the cross section, as a function of the corresponding b quark kinematics. We stress that this damping factor is not used when doing phenomenological calculations. It is needed here only to guarantee that the NLO cross sections computed with the traditional POWHEG BOX implementation is finite and agrees formally with the one computed with the new method. In order to perform the NLO test we did not need the virtual contributions, since they are the same in the two approaches. In the "traditional" implementation the Born phase space was generated with importance sampling on the dominant decay chain bq → (t → b(W + → e + ν e ))q . It was found that the new implementation yielded better convergence, speeding up the calculation by about a factor of 2 or more. Very high statistics runs of both implementations were performed, and full numerical agreement was found for both the total cross section and the differential distributions.

Results and comparisons at the full shower level
As mentioned earlier, the process we are considering is singular when final state b quarks have very small transverse momenta. Thus, event generation requires a generation cut or a Born suppression factor [18]. We adopt the latter method in the present work, using a suppression factor of the form As a result, less and less events are generated as the b transverse momentum becomes small, but the event weight is increased correspondingly, thus yielding a potentially divergent cross section for observables that do not suppress the contribution of small transverse momentum b quarks. The shower was performed using Pythia8 version 81.85 [19][20][21]. In all cases, Pythia8 was run with its default initialization, supplemented with the following calls: pythia.readString("SpaceShower:pTmaxMatch = 1"); pythia.readString("TimeShower:pTmaxMatch = 1"); pythia.readString("SpaceShower:QEDshowerByQ = off"); // From quarks. pythia.readString("SpaceShower:QEDshowerByL = off"); // From Leptons. pythia.readString("TimeShower:QEDshowerByQ = off"); // From quarks. pythia.readString("TimeShower:QEDshowerByL = off"); // From Leptons.
In order to test our generator, we generated four samples of one million of events each, and compared the relative output. The samples are obtained as follows: • NORES Sample. This is obtained using the traditional POWHEG-BOX-V2 implementation for the process. The events are fed to Pythia8, with the setting listed earlier.
Pythia8 is required to veto radiations at scales harder than the value of the Les Houches variable scalup, set equal to the transverse momentum of the POWHEG generated radiation for each Les Houches event.
• RES-HR Sample. This sample is obtained using the POWHEG-BOX-RES implementation of the process. The Les Houches events include the hardest radiation generated by POWHEG (the HR in RES-HR stands for "hardest radiation"). Since POWHEG is generating the hardest radiation, besides vetoing radiation in production with the usual scalup mechanism, we also forbid any Pythia8 radiation from top decays harder than scalup. We do this by explicitly examining the showered events. If a radiation generated by Pythia8 in top decay has a transverse momentum greater than scalup, the program discards it, and runs Pythia8 again on the same Les Houches partonic event. This procedure is repeated indefinitely, thus explicitly vetoing any event with radiation harder than scalup.
• RES-AR Sample. This sample is obtained using the POWHEG-BOX-RES implementation of the process. However, rather than keeping only the hardest radiation, we kept both the hardest radiation in top decay and the hardest radiation in production (the AR in RES-AR stands for "all radiation") . These radiations are composed into a single event, using the usual POWHEG mapping mechanism. In this case, besides the normal scalup veto for radiation in production, we must forbid Pythia8 radiation in top decays if harder than the POWHEG generated one in top decay. There are thus two different veto scales, one for production (i.e. the scalup value) and one for decay. There are no provisions in the Les Houches Interface for User Processes to store the radiation scale for decaying resonances. The program thus computes explicitly the transverse momentum of radiation in top decay at the Les Houches level. If the showered event contains shower generated radiation in top decay harder than this computed scale, the shower is discarded, and a new shower is generated on the same Les Houches event, repeating the procedure indefinitely until the event passes the required condition.
The RES-AR implementation is fully analogous to the allrad procedure illustrated in ref. [14]. The method (and the software) for vetoing Pythia8 radiation is also borrowed from that reference.
• ST-tch Sample. This sample is generated using the ST-tch POWHEG generator of ref. [22]. Radiation in decay is not included in this generator, and thus we let Pythia8 shower the event according to its default setup, vetoing events with radiation in production harder than scalup, and with no veto in top decays. In order to match more closely what we include in RES-HR, we deleted in the ST-tch the real processes initiated by a light (i.e. not b) quark and a gluon.

Phenomenological analysis
We have considered the LHC 8 TeV configuration for our phenomenological runs. We have used throughout the MSTW2008 set [23] at NLO order. Other PDF sets, like those of refs. [24,25], can be used as well, but we are not interested in a PDF comparison in the present study. We only consider the b µ + ν µ final state (i.e. not the conjugate one). We set the top mass to m t = 172.5 GeV. For this value of the mass and PDF choice the computed top width, including NLO strong corrections, is 1.3306 GeV. We use the same NLO value of the width also in the ST-tch generator. In this generator it only affects the top line-shape, since the cross section is determined by the top cross section multiplied by a user supplied branching fraction. On the other hand, in our generator the width must be computed with the same Standard Model parameters that are used in the matrix elements, since the cross section will be proportional to a partial width (depending upon the couplings that are used in the matrix elements) divided by the total width that appears in the denominator of the top propagator.
Since we will compare generators that do not include top resonance information, our analysis will be performed (unless explicitly stated otherwise) without using "Monte Carlo truth" information as far as the top particle is concerned, thus relying solely upon a particle level reconstruction of the top kinematics. We thus define the following objects: • The lepton. This is the hardest µ + in the event. • The neutrino. This is the hardest ν µ in the event.
• The W + . This is the system formed by the lepton and the neutrino.
• The b hadron. This is the hardest hadron with a b quark content (not ab!).
• The b-jet. Jets are reconstructed using the anti-k t algorithm [26], as implemented in the FastJet package [27], with R = 0.5. The b-jet is the jet that contains the b hadron.
• The top quark. This is defined as the system comprising the W and the b-jet. Only b-jets with a p t of at least 25 GeV and |η| < 4.5 are accepted for this purpose. In case such a b-jet is not found, no top is reconstructed.

RES-AR and ST-tch comparison
We begin by comparing the RES-AR and the ST-tch results. For this purpose (and only in this case) we have deleted from the RES-AR generator the real amplitudes with the final non-b light partons in a colour singlet. This excludes in particular contributions with tW − associated production, leading to a more meaningful comparison and to a better agreement on the total cross sections, since these contributions are not included in the ST-tch generator. In fig. 6 we plot the transverse momentum and the rapidity distributions of the reconstructed top. In these plots no requirement is made on the mass of the reconstructed top. As one can see, reasonable agreement is found with these distributions. In fig. 7 we show the mass peak, both for the reconstructed top and for the top particle in the Monte Carlo record (more specifically, we pick the last top in the Monte Carlo event record). It is apparent that the line-shape of the reconstructed top are not in complete agreement. Assuming that this is due to differences in the structure of the b-jet, we plot in fig. 8 the b-jet mass and profile, for b-jets with transverse momentum above 15 GeV and |η| < 4.5. The jet profile is defined as follows where N is chosen in such a way that where 0.5 is the ∆ R value that defines the b-jet. Thus, for ∆ R < 0.5, P b -jet (∆ R ) is the fractional distribution of the transverse momentum in the jet. In fig. 9 we show the transverse momentum of the b-jet. We see that these plots show consistently that the b-jet is harder and more massive in the RES-AR case than in the ST-tch one. In particular, the jet profile plot shows that there are more partons sharing the jet momentum in the region with ∆ R near 0.1 in the RES-AR case. In the ST-tch case more momentum is concentrated at very small ∆ R (presumably due to the b meson), and there are more partons at larger values of ∆ R , also outside of the jet cone. We interpret these fact as consistent with the reconstructed top mass peak being slightly shifted to the right in the RES-AR case.
In order to quantify the shift in the mass extraction that one would get using one or the other Monte Carlo, we define an observable M trec , equal to the average value of the reconstructed top mass in a window of ±15 GeV around m t , in order to mimic the typical experimental resolution on the reconstructed top mass. We get M trec = 170.54(2) GeV for the RES-AR, and M trec = 169.59(1) GeV for the ST-tch generator. Thus, extracting the top mass with the ST-tch generator we would get a value 1 GeV larger than if we used the RES-AR one.
As a further comment on our findings, we remind the reader that, as far as the reconstructed top line-shape is concerned, the RES-AR and ST-tch generators differ mainly in the way that radiation from the b quark is treated. In the RES-AR generator the hardest radiation from the b quark is always handled by POWHEG, with Pythia8 handling the remaining radiation. In the ST-tch generator, on the other hand, POWHEG generates no radiation from the decaying top. Thus, all radiation from the b quark is handled by Pythia8. We must therefore ascribe the differences that we find to the different treatment of radiation in POWHEG and Pythia8. This issue was also discussed in ref. [14]. In view of its impact on the top mass determination, this topic deserves a more detailed phenomenological study, that goes beyond the scope of the present work.

RES-AR and RES-HR comparison
We now compare the RES-AR and RES-HR generators. As mentioned earlier, the two generators differ in the way that the Les Houches record is formed after the stage of generation of radiation. In the former, the hardest radiation is kept for both production and top decay independently. So, events with up to two more partons with respect to the Born kinematics are stored in the Les Houches record, and are passed to Pythia8 for showering. RES-HR generator the b-jet hardest radiation is in part controlled by POWHEG and in part by Pythia8. We remark that also in the RES-HR case POWHEG should correct the hardest radiation from the b quark to yield an NLO accurate result, at least for sufficiently hard radiation.
It is again interesting to quantify the shift in the mass extraction that one would get using one or the other Monte Carlo. Computing, as before, our M trec observable, we get M trec = 170.06(3) GeV for the RES-HR generator, and M trec = 170.55(2) GeV for the RES-AR. The small difference in the RES-AR result with respect to the one given in the previous subsection was due to the fact that there a set of real contributions was left out, as explained earlier.

NORES and RES-HR comparison
We now compare the NORES and RES-HR generators. The purpose of this comparison is to see if and how a generator that is not aware of resonance structures can exhibit visible distortions. In this case, MC resonance truth information is not available for the NORES generator. It is possible however to try to guess the resonance information from the structure of the event, as we will detail later. We begin by showing results with the NORES contribution without performing any resonance assignment. In fig. 14 we plot the transverse momentum and the rapidity distributions of the reconstructed top. In fig. 15 we show the mass peak for the reconstructed top. In fig. 16 we plot the b-jet mass and profile. In fig. 17 we show the b-jet transverse momentum.
We see marked distortions in the mass peak, in the b-jet mass and profile, and in the transverse momentum distribution. In this case we get M trec = 170.54(2) GeV for the NORES and M trec = 170.06(3) GeV for the RES-HR generator. We remark that in this case no MC-truth was available for the top quark mass in the NORES case, and thus the corresponding plot is missing. Finally, we performed a guess resonance assignment on the NORES output record, in the following way: • The b µ + ν µ system is assigned to the top in all events.
• If the radiated parton is a gluon, and the g b system has the colour of a quark, then we compute the transverse momentum of the gluon relative to the beam axis, k T,isr , relative to the final state light quark, k T,fsr , and relative to the b quark in the g b µ + ν µ frame, k T,b . Furthermore we compute the quantities (5.13) The gluon is assigned or not assigned to the top resonance with probabilities propor- (5.14) • If the radiated parton is not a gluon, it is not assigned to the top.
We label as NORES-i (i for "improved") the corresponding generator, and show its comparison with the RES-HR output in figs. 18 through 21. We see that now the differences are much less pronounced, and do not seem to affect significantly the determination of the top mass. In fact, the average value of our M trec observable is now M trec = 170.07(2) GeV for the NORES-i and M trec = 170.06(3) GeV for the RES-HR generator. Mild distortions are observed for the remaining distributions. On the other hand, we see that the MC-truth top line-shape seems to exhibit unphysical features, presumably due to the way that resonance assignment was performed.

Conclusions
In this work we have presented a formalism for dealing with intermediate resonances in NLO+PS generators built within the POWHEG framework. We have formulated a subtraction method such that no double logarithms of the resonance's width arise separately in the integrated real and in the soft-virtual term of the NLO calculation. Single logarithms of the widths do however arise in the soft-virtual term (in fact in the virtual contribution) and in the integrated real cross section, and cancel only when adding them up. Thus, in the framework of a POWHEG generator, these single log terms cancel in theB function, so that both double and single logarithms of the resonance widths are absent there. In POWHEG, the generation of radiation is unitarized by construction. Therefore, all soft divergences (including also those that are cut-off by the resonance width) are properly regulated there by Sudakov form factors and/or by finite width effects.
Our formalism is fully general, and has been implemented in a general way in a mod-ified version of the POWHEG-BOX-V2. The framework is such that in order to implement a specific process, one must supply the Born, including spin and colour correlated amplitudes, the virtual and the real amplitudes. The framework takes care of everything else: it finds the possible resonance histories and singular regions, it builds a Born phase space consistent with the resonance histories, and it constructs the subtraction terms, the softvirtual contributions and the mismatch terms described in sec. 3. It then performs the various stages of the POWHEG event generation.
In the present work we have considered, as a test case, the process pp → µ + ν µ j b j, that is dominated by single-top t-channel production, in the 5-flavours scheme. Within this framework, we have examined the output of our generator, with particular attention to observables that can have an impact on the top mass measurements.
We have compared two variants of our generators. In the first one we use the traditional POWHEG method, dubbed RES-HR in this paper, retaining only the hardest radiation, feeding the corresponding partonic event to a shower Monte Carlo, and vetoing any showergenerated radiation harder than the POWHEG one. In the second one, dubbed RES-AR, the hardest radiations in production and in resonance decays are both kept and combined in the final partonic event. In this case, the veto scales on the hardness of the Shower radiation in production and in decays are different, being set to the corresponding scales of radiation in POWHEG. We also compare our generator to the previous single-top, t-channel generator, the ST-tch process of ref. [22] in the POWHEG-BOX-V2. We can briefly summarize our findings as follows. We find differences in the reconstructed top, mostly due to the structure of the b-jets, and we ascribe these differences to the fact that the hardest radiation in the b-jet is fully determined by the Shower Monte Carlo in the ST-tch generator, it is in part determined by POWHEG and in part by the Shower Monte Carlo in the RES-HR generator, and it is fully determined by POWHEG in the RES-AR generator.
We have also considered the output of a generator using the same, full matrix elements for the pp → µ + ν µ j b j process that we have used in our new generator, implemented however in the traditional POWHEG-BOX-V2 framework, that is to say without resonance aware formalism. In this context we have again considered two alternative options. In the first one we pass the events to the Shower Monte Carlo without any resonance information. In the second one we reconstruct a most probable resonance history of the event based upon kinematics. The aim of such test is to search for distortions in the generated radiation due to the lack of proper treatment of resonance decays. We have found that if no resonance information is passed to the shower important differences are in fact observed. On the other hand, if we make an educated guess of the event resonance history, and pass it to the shower, smaller differences are present at the reconstructed level, although the top line-shape at the MC-truth level exhibits unphysical features.
We remark that one comparison is still missing in the present work: the comparison to a generator using the on-shell approximation (i.e. a single top generator analogous to the ttb_NLO_dec tt generator of ref. [14]), in order to check if off-shell and non-resonant contributions at the radiation level are relevant for an accurate simulation of the production of a top resonance.
The relevant code for the present work is available in the POWHEG BOX repository at the url http://powhegbox.mib.infn.it/POWHEG-BOX-RES-beta. It has to be considered as very preliminary at the present stage, since only one relatively simple process has been implemented with it so far.