Heavy quark radiation in NLO+PS POWHEG generators

In this paper we deal with radiation from heavy quarks in the context of next-to-leading order calculations matched to parton shower generators. A new algorithm for radiation from massive quarks is presented that has considerable advantages over the one previously employed. We implement the algorithm in the framework of the ${\tt POWHEG-BOX}$, and compare it with the previous one in the case of the ${\tt hvq}$ generator for bottom production in hadronic collisions, and in the case of the ${\tt bb4l}$ generator for top production and decay.


Introduction
The production and detection of bottom quarks play an important rôle in various contexts in LHC physics. Letting aside the very abundant direct production, that is exploited for flavour physics studies, bottom is used to identify top particles and to study their properties. Furthermore, it is the dominant decay mode of the Higgs boson, that can be used to study processes as the associate HV production [1,2] and the large transverse momentum production [3]. In searches for physics beyond the Standard Model, bottom also appears often produced in association with new-physics objects.
Having a mass much larger than the typical hadronic scales, bottom quark production is calculable in perturbative QCD. In cases when the transverse momentum involved in the production is large compared to its mass, as, for example, in high-energy e + e − annihilation, or in production at large transverse momentum in hadronic collisions, bottom can behave as a light parton, and give rise to a hadronic jet. Techniques for dealing with these regimes have been developed in the past [4], and have been applied to the LHC a e-mail: luca.buonocore@na.infn.it b e-mail: paolo.nason@mib.infn.it c e-mail: francesco.tramontano@na.infn.it case [5]. They allow for the computation of the transverse momentum spectrum of promptly produced b quarks at nextto-leading order in QCD, including the resummation of large logarithms of the ratio of the transverse momentum over the bottom mass up to next-to-leading-logarithmic accuracy. These large logarithms can arise both from initial state radiation, when, for instance, an incoming gluon splits into a bb pair, with one of the b undergoing a large-momentumtransfer collision with a parton from the target, and from final state radiation. In the last case, an outgoing gluon can split into a bb pair, or a directly produced b quark can emit a collinear gluon.
The large transverse momentum regime is treated consistently at the leading logarithmic level in parton shower generators. At the most basic level, heavy flavours are treated as light flavours, but with a shower cut-off scale of the order of the heavy quark mass. However, considerable work has been performed to better account for mass effects. In some generators, this is achieved by suitable modifications of the splitting kinematics and splitting kernels [6,7]. The Sherpa dipole shower [8] makes use of the Catani-Seymour dipoles for massive quarks [9]. The Catani-Seymour formalism is also used in the DIRE shower [10]. In ref. [11] a final state dipole-antenna shower for massive fermions is proposed, based upon the corresponding antenna subtraction formalism of ref. [12].
In next-to-leading order (NLO) calculations matched to Shower generators (NLO+PS) for heavy flavour production [13,14], one generally treats the heavy flavour as being very heavy. The heavy quark mass thus acts as a cut-off on collinear singularities, that are thus not resummed. This approach has in fact proven to be quite viable in heavy flavour production even at relatively large momentum transfer [5]. Consider, for example, heavy quark pair production in a POWHEG framework. By neglecting collinear singularities from heavy quarks, the only singular region that we have to consider has to do with initial state radiation involving only light partons. Since the POWHEG procedure guarantees that the matrix elements are given correctly for up to one hard radiation, gluon splitting, flavour excitation and radiation from the heavy flavour are included, so that the logarithmically enhanced terms are correctly reproduced at first order. Higher order leading logarithms, however, are not treated correctly. In particular, there are reasons to give an adequate treatment to final state radiation from a high transverse momentum bottom quark. In fact, this radiation process is intimately related to the physics of the bottom fragmentation function, and may have important effects in processes of considerable interest, like for example in top decay.
In the POWHEG-BOX framework, a facility for the treatment of collinear radiation from a heavy particle was set up in ref. [15], in the framework of electroweak corrections to W production. In that context, the purpose was to deal with mass effects in the final state radiation of photon from the lepton in W decays. The same framework is also appropriate for describing radiation from heavy quarks. In particular, it was adopted in refs. [16], where the ttb_NLO_dec generator was introduced, and in refs. [17] for the b_bbar_4l generator, for dealing with radiation from bottom quarks in top decays.
The purpose of the present paper is twofold: we present a new algorithm for radiation from a heavy quark, that has proven superior to the old one; furthermore we perform a thorough investigation of the behaviour of this component of the POWHEG generator, also by comparing the two methods, both in the framework of bottom quarks generated in top decay, and in inclusive bottom quark pair production. In the last case, such a study was never carried out.
The paper is organized as follows: in section 2 we describe the new algorithm, in section 3 we illustrate our phenomenological studies, and in section 4 we give our conclusions.
2 Description of the new algorithm 2.1 The POWHEG mapping for the massive emitter case Let us assume for definiteness to deal with a scattering process involving n partons in the final state at lowest order in perturbation theory. The generic point in the corresponding Born phase space, referred also as Born configuration, will be denoted with barred momenta The corresponding phase space volume element is given by Fig. 1: Kinematics for a real configuration: k n is the massive emitter, k n+1 is the radiated parton. y denotes the cosine of the angle between the two tri-vectors.
where q is the total incoming 4-momentum. 1 At Next-to-Leading order (NLO), one must also include processes of emission of one more real massless extra parton, resulting in a n + 1-body kinematics which we will denote as The singular regions of the real phase space are separated by means of suitable projection operators; in each of them, the radiated parton phase space is parametrised in terms of the FKS variables [18] (the notations p and p for a generic momentum p denote the tri-impulse and its modulus respectively) as shown in Fig. 1, where we have assumed that the emitter and the FKS partons are respectively the n-th and the n + 1-th parton. The rescaled energy ξ is related to the soft limit (ξ → 0), and the variable y to the collinear one (y → ±1). The definition of the azimuthal angle φ , in the POWHEG framework, departs from the standard FKS definition. It is taken as the polar angle of the splitting around the axis parallel to the momentum of the recoil system, in the rest frame where q = (q 0 , 0).
In what follows, we will construct a one-to-one map from a real configuration with radiation variables (ξ , y, φ ) into a Born one. This leads to a factorisation of the real phase space in term of Born and radiation variables.
The mapping can be reduced to the case of the map from a 3-body phase space into a 2-body one. Inserting into the (n + 1)-body phase space volume element the identities and the phase space is decomposed into a chain of two consecutive processes. With reference to Fig. 1, they are: the decay of a particle with momentum q into the 3-body system formed by the emitter k n , the FKS-parton k n+1 and the "recoil" system, with momentum and invariant mass followed by the decay of the latter into the other n − 1 particles. In formula, we have where We now focus on the 3-body process; under the action of the mapping, the k n and k n+1 partons will be replaced by a single parton with mass m and momentum k n . We define so that We fix the transformation by demanding k n k. Care must be taken to ensure the conservation of energy-momentum also for the resulting Born configuration. This is accomplished by performing a boost Λ in the direction k and defining We determine the velocity parameter β of the boost transformation from the mass-shell condition We get We define the other barred variables as Their mass relations are preserved by the boost transformation and, furthermore, we have which is the energy-momentum conservation for the Born configuration.

Inverse map
We now detail the construction of the inverse map, which is what is actually needed in the applications. Suppose that a Born event has been generated, i.e. the barred variables k i (i = 1, · · · , n) are given. Then, M 2 rec is obtained inverting eq. (13): We want to attach to it a radiation described by the radiation variables ξ , y and φ . For future convenience we introduce the largest allowed value for ξ The energy of the radiated parton is Energy conservation requires that where k 2 rec = k 2 n + k 2 n+1 + 2k n k n+1 y.
We can solve equation (21) for k n in a standard way, by bringing in turn each single square root on one side of the equation and squaring both sides. By doing this we actually find the solutions of all of the following equations for all possible combinations of the signs in front of the square root. The solutions are given by In order for them to exist, the argument of the square root must be positive. This leads to the bound with .
The last equality follows from the fact that We see that ξ (+) is a decreasing function of y 2 . Thus that is larger than the maximum value allowed by energy conservation. Thus, the corresponding k n values should be the solutions of one among equations (23) where some minus signs appear. On the other hand, ξ (−) (y) is an increasing function of y 2 , so that is perfectly acceptable. Furthermore, in the ξ < ξ (−) (y) case the value ξ = 0 is allowed, that lead to the solutions k (±) n = ±k 0 n satisfying eq. (21) with the correct signs of the square roots. Since the k (±) n must always satisfy one of the equations (23), and since they are smooth function of both ξ and y in their allowed range (that includes the ξ = 0 point), we infer by continuity that they satisfy equation (23).
Up to now we have not imposed the positivity of k n . On the other hand, negative k n values still have a physical interpretation, as illustrated in fig. 2. Thus, provided we interpret negative values of k n according to the construction of fig. 2, we have two solutions of equation (21). They are however related, since If we pick just one of them, we have a single-value map from the underlying Born configuration and the radiation variables ξ , y and φ to a real emission configuration. We pick For continuity, k (+) n (ξ , y) vanishes on the boundary line y > 0, ξ = ξ (−) (0) separating the positive and negative regions. The points lying on this curve are degenerate and correspond to the same real configuration with the emitter at rest in the partonic centre-of-mass frame. Apart from them, that constitute a set of zero measure, the map is well defined and bijective. The inverse map is well defined also on the boundary line y > 0, ξ = ξ (−) (0). This means that the corresponding jacobian vanishes on that curve. Then, the inverse map can be safely used both for the integration of the real differential cross section and for the generation of radiation. In fig. 3 we display the ξ , y kinematic region. We remark that the negative k (+) n (ξ , y) region includes neither soft nor collinear singularities, since ξ is large, and since the angular separation of the quark and the radiated gluon is larger than π/2. From now on we will drop the suffix (−) and will use ξ (y) and ξ (0) instead of ξ (−) (y) and ξ (−) (0).
In fig. 4 we show the partition o the kinematic region represented in the more familiar Dalitz plane. Notice that in the massless limit the physical region in the Dalitz plot develops an acute angle in the lower right, corner corresponding to the gluon being anticollinear with the b quark. Thus, the problematic region ξ > ξ (0) is not a singular one.

Full kinematic reconstruction of the real emission
So far, we have got the length of the tri-vectors k n and k n+1 . It is a standard kinematical problem to determine their directions in such a way that their sum k is parallel to k n . We do not enter in further details about it.
The last step is to calculate the β parameter of the boost transformation Λ , eq. (15), and to boost "back" the other barred momenta in the real event The above mapping allows us to write the (n+1)-body phase space element in the factorized form where we have expressed the radiation phase space in terms of the FKS variables with the jacobian function J(ξ , y, φ ) taking into account the change of variables involved in the transformation. In order to extract the jacobian, we have to manipulate and compare the l.h.s and the r.h.s of eq. (33). Recalling eq. (8), we perform the change of variables in the three-body phase space, eq. (9), In polar coordinates, we have and, using as reference direction that of k, where α is the angle between k n+1 and k and φ is the azimuthal angle taking k as the reference direction. Hence On the other hand, following the same arguments that led to eq. (8), we can split the barred Born phase space into a two-body phase space and the phase space of the system recoiling against the emitting parton Since k n = q−Λ k rec , the delta function in eq. (39) constrains the value of k rec to be Then, exploiting the Lorentz invariance of the phase space element, we have where the r.h.s and the l.h.s are related by the boost transformation Λ . In particular, we observe that so that the boost maps the argument of the delta function in the r.h.s into that of the delta function in the l.h.s. Inserting eq.(38) and eq.(39) into eq.(33) and using eq.(41), we get 6 By virtue of the mapping, the vectors k and k n are parallel so that in polar coordinates their angular elements are equal, dΩ = dΩ n . Then, from eq. (43) we have and we are left with the computation of the jacobian of the two-variable-transformation This transformation is implicitly defined by the relations where λ is the kinematical Kallen function: Applying the chain-rule for the derivative, it is straightforward to compute the jacobian. We get The final expression for the full jacobian J is thus Note that the denominator of J vanishes in two regions: when approaching the curve ξ = ξ (0) for y > 0, behaving as ξ (0) − ξ when approaching the curve ξ = ξ (y), as ξ (y) − ξ .
In the first case, the k 3 n term in the numerator vanishes simultaneously as (ξ (0) − ξ ) 3 . It follows that the jacobian vanishes as J ∼ (ξ (0) − ξ ) 2 for ξ → ξ (0) at fixed y > 0. This result is coherent with what has been argued above regarding the degenerate points corresponding to the configuration with the emitter parton at rest in the partonic centre-of-mass frame. In the second region, the jacobian develops an integrable singularity, that can be dealt with by importance sampling techniques in Monte Carlo integration.

Generation of radiation
The POWHEG master formula for the generation of radiation [19,20] is where t min is an infrared cutoff, and the NLO Sudakov form factor is given by In the case of a massless emitter, K ⊥ is a smooth function of the radiation variables, which is required to reduce to the transverse momentum in approaching the soft and collinear limits. For the massive case, in ref. [15] the following definition was proposed y phy denotes the cosine of the physical angle between the emitter and the emitted parton. 2 Eq. (52) has the remarkable property of reducing continuously to the transverse momentum in the massless limit. We assume it as our default scale choice.
According to the standard veto method, we look for a suitable upper bound function U of the integrand in the NLO Sudakov form factor, namely For the sake of simplicity, we have omitted the integration on the azimuthal angle dφ , which results in a constant 2π factor. We model the upper bound function on the asymptotic singular behavior of the real matrix element near the soft and collinear singularities. We recall that the jacobian of the mapping has a divergent behaviour near the curve ξ = ξ (y). The upper bound function should have a behaviour not weaker than the Jacobian near the singular regions, and furthermore, it should be simple enough to allow us to perform an analytical integration in the constrained radiation phase space given by the cut K 2 T > t. It is convenient to perform a change of integration variables from ξ , y to ξ , K 2 T . Indeed, it turns out that K 2 T is a monotonic decreasing function of y at fixed ξ , i.e. ∂ K 2 T /∂ y < 0. The inversion of this mapping is too complex to be performed analytically, 3 but easy to perform numerically. We find that the associated jacobian ∂ K 2 T /∂ y has a behaviour similar to that of the jacobian of the mapping J: so that in the new integration variables the integrand becomes U U must have a simple form, and must have the appropriate behaviour to act as an upper bound for the soft and collinear singularities of the real matrix element.

Upper bound function
The singular behaviour of the real matrix element squared is universal and can be extracted in a straightforward manner by means of the eikonal approximation. In terms of the radiation variables, we get with N a suitable normalization constant. On the other hand, in the soft limit, the jacobian of the mapping behaves as We must also take into account the behaviour in the soft limit of the jacobian term factorized in U: Putting all the three contributions together, we obtain the following expression of the upper bound function U A more complete analysis shows that mapping J is enhanced (although not divergent) at large ξ for y → −1. In order to get a more efficient upper bound, we add the factor 1 1−K 2 T /q 2 to the previous expression. Hence, our final choice for the upper bound function U is

Integral of the upper bound function
In order to integrate the upper bound function analytically, its domain of integration has to be suitably enlarged. This can be done by interpreting the R/B expression as being defined in the larger domain, but as vanishing outside of the physical domain. Since the veto procedure prescribes that a point generated according to the upper bound function should be accepted with a probability proportional to the value of the radiation function divided by the upper bound function, points generated outside the physical domain should always be vetoed according to the above interpretation. From eq. (52), we find the upper bound where β 0 is the velocity of the emitter in the underlying Born configuration (this follows from the fact that we always have β ≤ β 0 ), and we also find We thus take as our domain of integration the region in K T and ξ such that eqs. (63) and (64) are satisfied. We notice that in this way ξ can even become larger than 1. In practice, however, adding also the ξ < 1 or ξ < ξ max limit would render the integration more difficult, so we prefer to deal with it by vetoing. Defining the integral of the upper bound function is then where is the rapidity of the emitter in the underlying Born configuration. Given a number 0 < r < 1, the t value generated by solving the equation

Generation of radiation kinematics
The algorithm for generating the radiation variables proceeds as follows: 1. We set the initial scale t 0 = K 2 max .
2. We generate a uniform random number and get t from eq. (67). If t is below t min , no radiation is generated, and the event is emitted as is. 3. We pick a new uniform random number 0 < r < 1 and we generate a value for ξ as ξ = ξ m (t) exp(y 0 r ).
This is consistent with the distribution of ξ at fixed K 2 T according to eq. (66). 4. If ξ > ξ max , we set t 0 = t, and go back to the step 2. 5. If the veto condition is passed, given t and ξ , we solve numerically for y the implicit equation If a solution does not exist, we set t = t 0 and go back to step 2. 6. Now that ξ and y are available, we generate a random φ , and compute the ratio R = [R/BJ(ξ , y)]/U(ξ , y)], with U given in terms of U in eq. (56), and generate a new random number 0 < r < 1. If r > R we set t 0 = t and go back to the step 2. Otherwise, the event is accepted.

comparison in the bb4l case
We have compared results obtained with the new method presented here, with those obtained with the default POWHEG settings for the bb4l generator of ref. [17]. We found remarkable agreement between the two results for all the distributions that we have examined. Here we show only two of them, to convey the idea of the quality of the agreement. These results were obtained for the 8TeV LHC collider, using the MSTW2008 PDF [21] set for reference only (other sets could be used as well [22,23]). In our simulations we make the B hadrons stable. Jets are reconstructed using the Fastjet [24] implementation of the anti-k T algorithm [25] with R = 0.5. We denote as B (B) the hardest (i.e. largest p T ) b (b) flavoured hadron. The B (B) jet j B ( j¯B) is defined to be the jet that contains the hardest B (B). We discard events where the j B and j¯B coincide. The hardest e + (µ − ) and the hardest ν e (ν µ ) are paired to reconstruct the W + (W − ). The reconstructed top (antitop) quark is identified with the corresponding W + j B (W − jB) pair. We show the invariant mass of the W − b-jet system ( fig. 5) and the B fragmentation function in top decay (Fig. 6), as defined in ref. [17], i.e. the the B energy in the reconstructed top rest frame normalized to the maximum value that it can attain at the given top virtuality.
In the curves, the alt (for "alternative") label stands for our new implementation, while def (for "default") is the current   Fig. 6: B fragmentation function in top quark decay as defined in ref. [17], produced with the bb4l generator for the 8 TeV LHC. The default and alternative implementation of radiation from b quarks are compared. Fig. 7: Example diagrams for the three mechanism that give rise to log-enhanced contributions in heavy flavour production: a) final state radiation from a quark; b) gluon splitting; c) flavour excitation.
POWHEG default. As one can see, the agreement is very good. This also shows that details in the implementation of radiation from the b quark in top decays do not seem to have important impact on physical observables.
We found that the efficiency and the generation rate of the new implementation are comparable with those of the POWHEG default.

b production in hadronic collisions
In this section we study the available POWHEG implementations of radiation from massive quarks for the hvq generator [14], i.e. the default POWHEG implementation and our new one. In spite of the fact that the default formalism has been available for quite some time [15], no such study has been performed so far. We thus discuss it in this work, where we can also compare with our new implementation.
The hvq generator has been available for quite some time as a tool to generate top, bottom and charm pairs in hadronic collisions. It is designed to simulate correctly the production of a heavy flavour pair when the logarithm of the ratio of the transverse momentum of the heavy quark divided by its mass is not too large. This limitation arises because there are three mechanisms, depicted in figure 7, involving radiation from the final state quark, production of a heavy quark-antiquark pair via final state gluon splitting and the splitting of an initial state gluon into a heavy quark-antiquark pair (where one of the two quarks is scattered at large transverse momentum), that can generate large logarithms involving the mass of the heavy quark. In the inclusive cross section for the production of a heavy quark with a given p T , for example, they generate logarithms of p T /m (see ref. [26], eq. (5.1)). The last two mechanisms are commonly referred to as gluon splitting and flavour excitation. In spite of this, the hvq generator has also been used to model relatively large transverse momentum production of heavy flavours, as in ref. [5]. There, the transverse momentum distribution of the heavy flavoured hadron in hvq was compared with the more accurate (but less exclusive) FONLL prediction [4]. It was found to be in rather good agreement. However, the large uncertainties related to the  non-perturbative fragmentation of the heavy quark leads to the suspect that such agreement is at least in part accidental. We will now compare the results obtained with the default hvq generator, that we will label nol (for "no light", meaning that the heavy quark is treated as very heavy), that treats as singular regions only the radiation from massless partons (i.e. initial state radiation); hvq with the inclusion of the radiation from the heavy quark as a singular region will be labeled asl (for "as light", meaning that the heavy quark is treated as a light parton). Furthermore, the default treatment of the heavy quark radiation region will be denoted as def, while the new implementation presented here will be called alt. In figs. 8, 9 and 10 we show a compari- son of def and alt. We can immediately see that we do not find important differences between the two methods, consistently with what was found in the bb4l case. The settings are similar to the bb4l case: we make the B hadrons stable, and define the b (b) jets as the jets containing the hardest b (b) flavoured hadron, with the jets defined as in the bb4l case. However, we do not exclude the case when both hardest bflavoured hadrons are in the same jet. We perform the calculation for the LHC at 8 TeV, using NNPDF30_nlo_as_0118 pdf set [23]. As one can see, the two implementations are in excellent agreement. Observe the jump at 10 GeV in the j B mass. It is due to the case in which the b andb flavoured hadrons are both in the jet cone. From figure 10 we also see that for jet masses above 10 GeV the gluon splitting configuration dominates. We found that the new implementation has a generation efficiency, which is estimated from the numbers of vetoes in FSR generation, three times greater than the default one. This leads to a generation rate of 1316 events per minute, against the 298 events per minute of the POWHEG default, which corresponds to a gain more than a factor of 4.
We now show in the left panels of Figs. 11, 12 and 13 the comparison among the alt and nol. Here we see considerable differences, especially in the large-momentum tail of the B and j B transverse momentum distribution, the alt ones being much harder. The mass of the b jet is also remarkably different. The large difference above 10 GeV hints to the fact that heavy quark pair production via the splitting of a large transverse momentum gluon is treated in a very different way in the two cases, and that this difference may be the cause of the large discrepancy in the transverse momentum distribution of the b hadron.
The difference between the alt and nol cases should not come as a surprise. The generation of radiation is per- where k t is the transverse momentum of the emitted gluon with respect to the beam axis, since the only singular regions that are considered there are the initial-state radiation (ISR) ones. The strong coupling constant and the parton densities are evaluated by default at a scale equal to the transverse mass of the heavy quark at the level of the underlying Born kinematics in theB function, while they are evaluated at a scale k t (or k t ) in the R/B ratios appearing in formula (71). SinceB and B are of order α 2 S , while R is of order α 3 S , this means that in practice two powers of the strong coupling are evaluated at the scale of eq. (71), while one power is evaluated at a scale k t . The mismatch in the scale used inB and in the B appearing in the ratios, combined with the exponential, leads as usual to the correct Sudakov form factor for initial state emission.

Problematic regions
In case the transverse momentum of the gluon is small, the scale assignments and the Sudakov form factor describe the process appropriately. It can happen however, that the real emission kinematics is near the gluon splitting, flavour excitation or quark radiation regimes. In these cases the gluon transverse momentum is not small. Furthermore, the numerator R in the integrand may be enhanced with respect to the denominator, thus yielding a damping of the real cross section that is not justified. Also the scale choices are not appropriate. For example, in the case of production of a high transverse momentum heavy quark pair according to the gluon splitting mechanism, the appropriate scale should correspond to two powers of α S evaluated at the gluon transverse momentum, and one power of α S evaluated at the scale of the order of the invariant mass of the heavy quark pair.
The adoption of the methods illustrated in ref. [15] and in the present work for dealing with radiation from a heavy quark leads to the correct treatment of the radiation from the heavy, quark provided all remaining regions are treated correctly. This is in fact what happens in the case of the bb4l generator, where there is only one enhanced region, but it is not the case for the asl generator, that does not treat in a proper way the two regions of gluon splitting and flavour excitation. Thus, the nol and the asl generators will end up treating the enhanced regions in different (and in both cases incorrect) ways. In fact, while in the nol case the enhanced regions will all be treated as if they were ISR processes, in the asl case they will be split, and treated in part as ISR processes, and in part as radiation from the heavy quarks. In order to test this hypothesis, and in order to explore possible strategies to deal with this problem, we proceed as follows. It is possible in POWHEG to further separate out the real cross section into two terms, such that only one term has singular behaviour, while the remaining term, being finite, can be integrated independently. In the hvq case, this means Eq. (70) is then replaced by We can exploit this mechanism in order to separate out the enhanced regions, in such a way that we can treat them in a more uniform way with our generators. In particular, we separate out the gluon splitting and flavour excitation processes in all cases. In the nol case we also separate out the regions of radiation from the heavy quarks, in such a way that they are treated in a more transparent way. Observe that in performing this separation we rely upon the fact that the three enhanced region are not really singular, since the quark mass cuts off the collinear singularities, and thus the remnant term is actually finite. We define the distance of a real configuration from a given enhanced region as follows where in the first line the distances for ISR and gluon splitting are given, in the second line those for radiation from the heavy quarks, and in the last line the ones for flavour excitation. We then define, for the nol generator For the alt and def generators, we define where the index i labels the three singular regions that POWHEG is handling. In this case, the cross section is damped if the kinematics is near a singular region that is nether ISR nor FSR, i.e. only gluon splitting and flavour excitation kinematics are separated into the (r) component. There is one more issue that needs to be considered when using a damping factor in POWHEG. By default, when evaluating the R (r) component (called "real remnant"), the scale choice is the same as forB, i.e. it is eq. (71) applied to the underlying Born kinematics, that depends upon the considered singular region. This would lead to a different scale choice for the remnants in nol and asl. In order to avoid that, we should set the scale on the basis of the real kinematics. This can be done in POWHEG by setting appropriate flags and by modifying the code that computes the scales for the process. Our scale choice is that has the correct limit to the underlying Born scale both in the ISR and in the FSR case. The result of this procedure is shown in the right panels of Figs. 11, 12 and 13. We notice a remarkable improvement in the agreement, although some important differences do remain. This is not unexpected, since in the two cases radiation from the heavy quark is treated in a very different way. It is interesting to notice that the B and the j B spectra computed with the nol without remnants (which is the default in the standard hvq generator), is in fair agreement with the alt one when the enhanced regions are separated using the remnants. Since the default hvq program gives a description of the transverse momentum distribution of B hadrons that is in fair agreement with the FONLL calculation, we infer that also the alt prediction will display a similar agreement, provided the gluon splitting and flavour excitation region are treated separately as remnants.
The alt (or equivalently the def generator), with the remnant separation discussed above, seems to be at this point the generator that may give the best description of b production data at hadron collider. We should not forget, however, that some flexibility still remains in the treatment of the remnant (in this work we have made a definite scale choice for the remnants in order to have a clearer comparison with the nol generator). We also notice from figs. 11 and 12 that after the remnants are introduced, the B-hadron and b-jet p T spectra become softer. This seems to be in contrast with the discussion at the beginning of sec. 3.2.1. On the other hand, this result may be due to the particular scale choice that we have performed for the real graphs, and that POWHEG applies automatically also to the remnants. This scale turns out to be higher than the typical scale involved in the region discussed at the beginning of sec. 3.2.1. A better approach would be to introduce the possibility of alternative scale choices in the remnants, including the possibility of performing a different scale choice depending upon which enhanced region one is considering. We refrain here from carrying out such study, since we believe that comparison with data on single inclusive b-hadron and b-jet production (see ref. [27][28][29] and references therein) and on correlations of bb pairs [30,31], would be needed in order to make progress in this direction, and this is beyond the scope of the present work. Such study would however be very valuable, and not only for the purpose of testing QCD in bottom production. The production of top pairs is one of the most important background process at the LHC, including also its future high luminosity and eventually high energy developments. An accurate simulation of the tt background would be very valuable for the LHC experimental collaborations, and it is quite clear that a model of b production yielding a good description of the data from low to high transverse momenta should be also well suited to describe top production in all the needed phase space. Furthermore, these studies could lead to an improved description of top production at large transverse momentum, that, as reported currently by CMS and ATLAS (see [32,33] and references therein) seems to be not well described by theoretical models.

Conclusions
In this paper we have presented a method for implementing radiation from a heavy quark in the POWHEG framework. This method is considerably simpler and more transparent than the one presented in ref. [15], and it has a much better numerical performance. The present method overcomes a problem related to the fact that the most natural map from an underlying Born configuration and a set of FKS-like radiation variables for radiation from a massive quark does not have a unique inverse in the whole kinematic region. The POWHEG inverse map has thus two solutions in some region of phase space, and in the present work it is shown how, by giving up the physical connection of the y variable to the gluon emission angle in a very limited region of phase space one can pick one of the two solutions in such a way that we still have a single valued inverse mapping. We have examined the output of the new method in the framework of the generators of ref. [17] and [14]. We found that the new method yields results that are very consistent with the previous one, that is at this moment the POWHEG default. This is reassuring, since it shows that details of the implementation do not impact in a visible way the physics result, and also it shows that the previous implementation, in spite of being quite contrived, is in fact correct.
In this work we have also examined for the first time the impact of the inclusion of the singular region associated with radiation from the heavy quark in the case of the hvq generator. We have shown that, unless one separates the enhanced gluon splitting and flavour excitation contribution from the real contribution that are dealt with by the POWHEG radiation formula, and treats them as remnants, one finds results that are in considerable disagreement with the traditional hvq implementation. On the other hand, it seems that performing this separation is the appropriate thing to do for a consistent modeling of the process. We notice that such modeling would be of great interest also for its potential application to top pair production, that is a very important background to many LHC physics studies.