Implementation of electroweak corrections in the POWHEG BOX: single W production

We present a fully consistent implementation of electroweak and strong radiative corrections to single W hadroproduction in the POWHEG BOX framework, treating soft and collinear photon emissions on the same ground as coloured parton emissions. This framework can be easily extended to more complex electroweak processes. We describe how next-to-leading order (NLO) electroweak corrections are combined with the NLO QCD calculation, and show how they are interfaced to QCD and QED shower Monte Carlo. The resulting tool fills a gap in the literature and allows to study comprehensively the interplay of QCD and electroweak effects to W production using a single computational framework. Numerical comparisons with the predictions of the electroweak generator HORACE, as well as with existing results on the combination of electroweak and QCD corrections to W production, are shown for the LHC energies, to validate the reliability and accuracy of the approach


Introduction
The production of a high transverse momentum lepton-neutrino pair in hadronic collisions, a process known as charged current (CC) Drell-Yan (DY), represents one of the most relevant processes at the LHC because of its large cross section and clean signature. This reaction is very well suited for i) a precise determination of two fundamental parameters of the Standard Model (SM), i.e. the W -boson mass M W and width Γ W ; ii) constraining the parton distribution function (PDF) of the proton; iii) background studies in the search for new heavy resonances; iv) a possible determination of the collider luminosity.
Specifically, thanks to the very large statistics, the LHC might measure M W with an accuracy of about 15 MeV or possibly better, even if this measurement appears at present particularly challenging. More generally, the very large statistics accumulated at the LHC implies that the systematic errors, including the theoretical ones, play a dominant role in JHEP04(2012)037 the determination of the measurement error of the total cross sections, as well as of the other observables of experimental interest.
The luminosity and the PDFs are two of the main sources of systematic error for the LHC measurements. As an alternative to ordinary measurements of the PDFs, the process of single W production could provide a handle to constrain the PDFs themselves through the measurement of the W + /W − asymmetry, while the standard luminosity determinations could be cross-checked with the measurement of the total cross section of single vector boson (in particular Z) production.
Last but not least, in the high energy regions of the phase space, where the searches for new physics beyond the SM are focused, the CC DY is the main background to unravel signatures due to the presence of heavy charged bosons predicted by many extensions of the SM. In these high energy tails, the electroweak (EW) loop diagrams containing the exchange of massive gauge bosons give rise to large contributions to the experimental observables, because of the presence of EW Sudakov logarithms ∝ log 2 (ŝ/M 2 V ), log(ŝ/M 2 V ), whereŝ is the squared partonic center of mass (c.m.) energy and M V , V = W, Z the gauge boson mass.
Because of the above reasons, it is essential for the LHC physics programme to have accurate theoretical predictions available for the W production cross section and associated distributions.
In perturbative QCD (pQCD), leading order (LO) calculations have in general a large uncertainties due to the dependence upon the unphysical factorization and renormalization scales. These uncertainties can be reduced including at least pQCD next-to-leading/nextto-next-to-leading order (NLO/NNLO) corrections. The Drell-Yan process was one of the first to be calculated in pQCD at NLO accuracy [1], while NNLO predictions for the integrated cross section became available later in ref. [2]. The NNLO pQCD computation of DY observables in a fully differential form was completed only in recent years by two independent groups [3][4][5].
The size of the NNLO pQCD contributions is naively comparable to the one of the NLO EW radiative corrections, fully computed in refs. [6][7][8][9][10]. The EW corrections are therefore a necessary ingredient in a precise theoretical description of the DY observables. For example, it is known that they induce a shift on M W of the order of 100 MeV [11] and introduce a negative correction to the transverse mass distribution M W T of about 20-30% for M W T above 1 TeV [10]. In particular, realistic phenomenological studies and data analysis require the simulation of all the experimental cuts and of the detector acceptance. Moreover the broad nature of the W physics programme at the LHC relies upon the measurement of a number of observables, which must be precisely predicted to avoid any theoretical bias in data interpretation. Hence, the implementation of the most important theoretical ingredients into a Monte Carlo (MC) event generator is mandatory. The presence of soft and collinear divergences due to the emission of coloured particles by the initial state partons leads any fixed-order calculation to provide accurate predictions only for quantities sufficiently inclusive over QCD radiation, like the total cross section, the W rapidity and transverse mass distribution. However, exclusive quantities, like the W or lepton transverse momentum, are often of great experimental interest, so that Shower Monte Carlo (SMC) programs, like JHEP04(2012)037 HERWIG [12], PYTHIA [13] or SHERPA [14], must be used. Beyond the pure Parton Shower approximation, the generators POWHEG [15,16] and MC@NLO [17] are widely employed at the LHC because they guarantee a better accuracy, being based on a consistent matching of NLO QCD corrections with SMC codes.
Till now, EW and QCD corrections have been typically implemented in event generators separately. Fixed-order QCD programs available for W production are MCFM [18,19] at NLO accuracy, and FEWZ [4,20] and DYNNLO [5] at NNLO. The ResBos code [21] is based on analytical resummation of all transverse momentum logarithms with NNLO accuracy. As already said, MC@NLO [17] and the POWHEG BOX [22] combine NLO QCD corrections with SMC, albeit according to a different methodology. Programs implementing exact NLO EW corrections are DK [6], WGRAD [8], SANC [9] and HORACE [10]. In HORACE NLO EW corrections are matched with a QED Parton Shower [11], while in WINHAC [23] the simulation of multiple photon effects, realized through the YFS framework, is associated with NLO EW contributions through an interface to the SANC module [24]. A QED shower for final state particles can be also described by other codes, using e.g. a universal and widely used tool like PHOTOS [25].
Only during the last few years the interplay of QCD and EW corrections to the CC DY observables has been studied in some detail. In ref. [26] the combined effect of QCD resummation and NLO QED final-state corrections to W boson observables has been addressed, while a more thorough analysis of the EW⊗QCD interplay has been performed in ref. [27], using, among others, the HORACE and MC@NLO programs. First attempts towards a single framework implementation of NLO EW corrections in association with QED and QCD showers are documented in ref. [28] and ref. [29]. The EW corrections to W +jet production (with and without W leptonic decays), which are part of the O(α em α s ) corrections to inclusive W production, were computed in refs. [30][31][32][33]. These studies outlined the importance of the simultaneous control of all the relevant QCD and EW corrections for a sensible description of W boson observables.
In the light of the above motivations, it would be desirable to have at our disposal • a MC event generator for W production • with NLO QCD and NLO EW corrections • interfaced to QCD and QED SMC.
Such a tool is presently unavailable. 1 We have built this program, and the aim of the paper is describing the theoretical and technical steps followed to achieve this goal. We also present first numerical results to cross-check the correct calculation and codification of all the different ingredients. Since we are generally interested in extending the present approach to processes other than W production, we have chosen the POWHEG BOX as a general computer framework for implementing our generator.

JHEP04(2012)037
The rest of the paper is structured as follows. In section 2 we shortly review the POWHEG BOX framework, and summarize the basic modifications needed for a consistent inclusion of the EW effects. In section 3 the main details of the full NLO EW computation are described, as well as its concrete implementation in the POWHEG BOX. In section 4 we sketch the modifications introduced in the POWHEG BOX to take care of photon radiation. In the same section we also illustrate the interface of NLO corrections with a mixed QCD⊗QED shower. A sample of numerical results, both for NLO EW corrections and for the full combination of QCD and EW contributions, in shown in section 5. Conclusions and perspectives are drawn in section 6.
2 The POWHEG BOX: framework and basic electroweak modifications The POWHEG method (Positive Weight Hardest Emission Generator), first introduced in ref. [15], has been used so far to implement a wide variety of NLO QCD corrections into a NLO + Parton Shower matching framework. 2 It has been described in great detail in ref. [16]. Recently the POWHEG algorithm has been codified into the POWHEG BOX [22], which is a general computer framework for implementing NLO calculations in SMC programs according to the POWHEG method. Using the POWHEG BOX, the implementation of a QCD NLO process requires the following ingredients: • the squared Born matrix elements B; • the colour correlated Born matrix elements B ij ; • the Born phase space; • the squared matrix element R of the real emission process; • the finite part of the virtual correction, V, computed in dimensional regularization and divided by the factor N = (4π) . As specified in the next section, and µ are the usual parameters of the calculation in dimensional regularization, and Q 2 is the squared virtuality scale of the process under consideration.
The remaining parts of the NLO calculation is performed by POWHEG BOX itself, that • identifies all the singular regions (corresponding to soft and collinear emission); • projects the real emission contributions over the singular regions; • implements the subtraction procedure; • computes the soft and collinear remnants; • generates the event with the hardest radiation.

JHEP04(2012)037
In order to include EW corrections together with the NLO QCD ones, several extensions of the POWHEG BOX were needed. In fact, collinear singularities due to the emission of photons must also be consistently included, and one should also foresee that the hardest radiation could be given by a large transverse momentum photon, rather than a coloured parton. 3 In the process we are considering, the only relevant QCD singular region is the one of small transverse momentum gluon radiation, i.e. is the initial state radiation (ISR). We must now also consider the EW singular regions, i.e. the possibility that also a photon may be produced at small transverse momentum by the incoming charged quarks. Furthermore, the final state lepton may emit a collinear photon. Thus, while in W production with only NLO QCD corrections there is only one singular configuration, including EW corrections we end up with three regions: ISR gluon emission, ISR photon emission and final state radiation (FSR) photon emission. The general POWHEG formula when more than one singular region is present is given by eqs. (4.14), (4.16) and (4.17) of ref. [16], that we report here for convenience and 3) The functionB f b (Φ n ) is the NLO inclusive cross section at fixed underlying Born flavour f b and kinematics Φ n . The real contributions are separated into terms labelled by the index α r . Each α r denotes a single flavour structure, and a single singular region. Each term R αr is singular only in the singular region denoted by α r . In our case, for example, the real contribution with flavour structure dū → (W − → eν)g is singular only in the region of small transverse momentum of the gluon, i.e. in the ISR region. Thus it corresponds to a single α r . On the other hand, dū → (W − → eν)γ is singular in the JHEP04(2012)037 region of small transverse momentum of the photon, and in the region where the photon transverse momentum with respect to (w.r.t.) the electron is small, which is a FSR region.
POWHEG separates the real contribution into two terms denoted by a different α r , one term being singular only in the ISR, and the other in the FSR region. The notation α r ∈ {α r |f b } means all the real singular contributions that have f b as underlying Born flavour. The square brackets with subscript α r and superscriptΦ αr n = Φ n mean that everything inside refers to the particular real contribution labelled by α r , and having underlying Born kinematics equal to Φ n .
The Sudakov form factor, eq. (2.3), is a product of individual Sudakov form factors associated with each α r . The POWHEG BOX is already capable of generating radiation using this Sudakov form factor, using the highest bid method (see ref. [16]). The generation of radiation must however be improved to allow for the generation of photon emission.
In summary, the POWHEG BOX was modified in the following points: • the term V in eq. (2.2) is the sum of the soft and virtual contributions. The soft contribution is evaluated by the POWHEG BOX, and now it must also include the terms arising from soft or collinear photons.
• The routines that automatically search for the singular regions in the real contributions have been modified, so that also a photon-quark and a photon-lepton pair of external lines are recognized as a singular region.
• The routines that compute the soft, collinear and soft-collinear limit of the real amplitude had also to be modified with the inclusion of the soft and collinear photon regions.
• The routine that computes the collinear remnants must also deal with photons.
• The computation of the EW virtual corrections are done keeping the mass of the lepton finite. FSR in the original POWHEG BOX considers only strictly massless particles. In fact, massive quarks in the POWHEG BOX are treated as really heavy, and not giving rise to collinear divergences. In the present case, the electron is really too light for this approach to be sensible. Thus, an extension of the POWHEG BOX to allow for collinear radiation from massive particles has been set up. This feature is not specific to EW corrections and can also be adopted for not so heavy quarks.
Notice that we also include the real processes with an incoming gluon, like gd → (W − → eν)u.
We may wonder whether the analogous EW process γd → (W − → eν)u (the so-called photon-induced process) should be included. Although it would be not difficult to include it, we decided to leave it out, basically for two reasons. The first reason is that there are not many PDF sets that include electromagnetic evolution, as further remarked in section 3.3. The second reason is that the photon density is already very small, of the order of α em , and it multiplies a process of the order of α em in this case. One should in fact remember that the enhancement of initial state electromagnetic radiation is a logarithmic factor of log E/λ, where λ is the collinear cutoff.

JHEP04(2012)037
While for an electron beam λ is the mass of the electron, in a hadron beam λ is of the order of a typical hadronic scale, i.e. few hundred MeV, three orders of magnitudes larger than the electron mass. Thus, the one order of magnitude enhancement of ISR in processes with incoming leptons is reduced to a factor of a few for hadron beams. These arguments are supported by the numerical evaluation of photon-induced processes given in refs. [34,35].
In order to be compliant with the general NLO QCD structure of virtual corrections in the POWHEG BOX, it is convenient to calculate also the virtual EW contributions in dimensional regularization of infrared divergences, which is not the scheme under which perturbative EW calculations are usually carried out in the literature. Moreover, such an approach simplifies, at least in principle, the extension of EW corrections to more complicated processes through the adoption of available automated procedures. The calculation of the EW virtual corrections in dimensional regularization is described in section 3.
To summarize, the newly developed tool • ensures normalization with NLO QCD + EW accuracy • combines the complete SM NLO corrections with a mixed QCD⊗QED parton cascade, where the particles present in the shower are coloured particles or photons • consequently, incorporates mixed O(α em α s ) contributions with a better accuracy w.r.t. existing programs. In particular, it can allow to study consistently the interplay between QCD and EW radiation, like e.g. the link between a photon emitted after QCD radiation and viceversa.

NLO electroweak calculation
In the present section we describe the main features of the NLO EW calculation. We first detail the treatment of the virtual contributions, then we discuss finite-width effects, and finally we discuss the subtraction procedure of infrared and collinear divergencies.

Virtual contributions in dimensional regularization
The typical framework used in NLO EW calculations to regularize infrared (IR) and collinear singularities in the loop integrals is the procedure known as mass regularization.
The IR soft singularity is regularized by introducing an infinitesimal photon mass λ. The mass of the lepton acts as physical regulator of the lepton-photon collinear singularities. If quarks take part to the process under consideration, small quark masses are introduced as regulators of the associated collinear singularities. 4 The IR divergences and part of the quark mass singularities cancel in the combination of virtual and real corrections, while the remaining collinear singularities, due to photon radiation off the initial state quarks, are absorbed into the PDFs. Large logarithms of the form log(ŝ/m 2 l ), with m l lepton mass, originated from the emission of collinear photons from leptons, as those coming from the JHEP04(2012)037 decays of the W/Z bosons in the DY process, survive in not fully-inclusive observables. Ultraviolet (UV) divergences are universally treated in dimensional regularization.
In QCD calculations dimensional regularization is used to regularize UV, IR and collinear divergencies. This scheme is implemented in the POWHEG BOX as well. In order to make the treatment of soft and collinear singularities uniform in the QCD and EW sector, we treated the IR and collinear divergencies of EW nature according to the dimensional regularization as well. More precisely, we use a sort of hybrid scheme, adopting dimensional regularization for the singularities associated to the coloured charged particles and the photon, but keeping the mass of the leptons finite, since this mass has a well-defined physical meaning, i.e. it is the true physical regulator of QED mass singularities.
In practice, the squared matrix element (ME) for the CC DY process qq → lν l , including one-loop QCD and EW virtual corrections, can be written as: where M 0 is the ME in the LO approximation and M QCD 1 ≡ δ QCD M 0 is the ME with NLO QCD virtual corrections, calculated in ref. [1] and already present in the available version of the POWHEG BOX. In eq. (3.1) M EW 1 ≡ δ EW M 0 is the novel ingredient, i.e. the NLO EW one-loop ME. Instead of doing ab initio a new complete calculation of the EW virtual corrections, we chose, for simplicity, to rely upon the one-loop structure yielding δ EW as given in ref. [6]. The latter calculation was indeed cross-checked over the years against the predictions of various independent EW codes and found in perfect agreement. According to ref. [6], δ EW is a linear combination of the fundamental 't Hooft-Veltman scalar functions B 0 , C 0 and D 0 [36] and relevant derivatives. Therefore, their expressions in dimensional regularization can be obtained by calculating each divergent scalar function in n = 4 − 2 dimensions. In particular, we calculated the B 0 functions for different energy scale and mass combinations, while for the C 0 and D 0 we referred to ref. [37] and to ref. [38], respectively. In particular we resorted to eqs. (B.2), (B.4), (B.12) and (B.16) from ref. [37] and used eq. (3.78) and eq. (4.19) from ref. [38]. For consistency with the FKS [39] subtraction procedure already present the POWHEG BOX, each scalar function was divided by the overall factor [22] where µ EW is the (arbitrary) mass scale of dimensional regularization and Q EW the squared energy scale of the process. 5 Before matching the virtual corrections with the standard (FKS) subtraction procedure of the POWHEG BOX, we verified, as an intermediate but not trivial check, that all the 1/ 2 and 1/ poles, as well as the arbitrary scale µ EW , cancel out up to numerical precision when adding to the virtual part massive [43] dipole subtraction formulae [44] improved by the inclusion of the singular collinear contributions from the PDFs. Also cancellation of the UV divergences was successfully checked. 5 We checked that the implemented library for the calculation of the scalar form factors agrees with the output of LoopTools [40] for the regular scalar functions and with that of GOLEM [41,42] for the IR ones.

JHEP04(2012)037
We performed further stringent internal cross checks, in order to test also the finite part of the loop corrections. In particular, in order to check the accuracy of our implementation of δ EW in the POWHEG BOX program, we calculated the basic scalar functions in the mass regularization scheme and compared, point by point of the phase space, the corresponding numerical values of the virtual part of δ EW with the ones returned by the HORACE code, which relies upon the same mass regularization scheme. We observed that the two programs perfectly agree, at the level of one part over 10 8 everywhere in the phase space with the only exception of the "low" (below ∼ 10 GeV) invariant mass region. This disagreement can be easily understood because HORACE includes part of the finite mass terms in the virtual correction, while δ EW does not. However, the observed small discrepancy can be considered acceptable, the low invariant mass region giving a negligible contribution to the total cross section and being practically unimportant for physics studies at the LHC. The numerical comparisons between HORACE and the POWHEG BOX at the level of full NLO EW corrections shown in section 5 support this expectation.
Concerning renormalization and the choice of the EW input parameters, it is known that a particularly convenient scheme for the calculation of the EW corrections to the CC DY is provided by the so-called G µ scheme, wherein G µ (the Fermi constant), M W (the W mass) and M Z (the Z mass) play the role of primary quantities. Actually, this choice allows to remove the dependence of the results on the masses of the light quarks entering the self-energy loop diagrams and responsible for the transition from the fine-structure constant α em (0) to the running electromagnetic coupling α em (Q 2 ) at a high-energy scale Q 2 , e.g. Q 2 = M 2 Z . In the G µ scheme the weak coupling constant g is calculated in terms of G µ and sin θ W instead of α em (0) and sin θ W . Consequently, the G µ parametrization of the Born cross section minimizes the EW correction, since the universal corrections induced by the running of α em (as well as by the ρ parameter) are absorbed in the LO amplitude. However, in the code we leave to the user the possibility of choosing between the G µ and α em (0) scheme with the flag scheme in the input file. 6 Obviously, for photon radiation from the charged particles we use the Thomson value α em = α em (0).

Finite-width effects
The relevant regions for DY physics are the regions around the W resonance and the very high energy tails of the distributions above the TeV scale. In particular, around the W peak the treatment of the EW corrections requires particular care because of the importance of this region for precision measurements and the presence of particular logarithms in the electroweak factors. Indeed, δ EW contains logarithms of the form Since these singularities are cured by a Dyson resummation of the W self-energy loop diagrams, what is usually done 6 In the αem(0) scheme, wherein the input parameters are αem(0), MW and MZ , the masses of the light quarks (kept at a finite value only in the gauge invariant subset of fermionic corrections) are chosen in such a way to reproduce the hadronic contribution to the photon vacuum polarization. 7 We consider here for brevity of notation s like the nominal squared c.m. energy entering the partonlevel EW calculation before convolution with the PDFs. It is understood that in the actual hadron-level calculation s is replaced byŝ = x1x2s.

JHEP04(2012)037
in most EW calculations is substituting each log(s − M 2 . This replacement can be properly applied in the one-loop virtual amplitude, because the coefficient of the log(s − M 2 W ) term is gauge invariant. We name this typically adopted procedure [6,10] Changing Logarithm Argument (CLA).
However, as widely discussed in the literature, the proper description of finite-width effects related to a resonance is a delicate point in perturbation theory and, in particular, in NLO EW calculations, that must provide gauge-invariant results.
To resolve this kind of problems at the level of one-loop calculations, in recent years a new, theoretically consistent scheme, known as Complex Mass Scheme (CMS) has been developed [45,46]. With CMS we have at our disposal a universal and easy to implement framework to describe finite-width effects in the phase space regions near and far-off the resonances. Gauge invariance, as well as all the Ward identities, are preserved in the CMS procedure, while unitarity is guaranteed up to O(α 2 em ). A drawback of CMS is that EW spurious corrections at the NNLO level are generated, but they are beyond the required NLO accuracy. In particular, the calculation of the two points scalar functions entering the renormalization constants in CMS implies their evaluation in terms of a complex momentum M 2 W − iΓ W M W . This demands an analytic continuation in the momentum variable to the unphysical Riemann sheet. We have avoided this difficulty as suggested in ref. [45], i.e. by expanding the self-energies appearing in the renormalization constants about real arguments. One-loop accuracy is conserved.
In our calculation we decided to leave the possibility to the user of choosing between CMS and CLA schemes (the option complexmasses in the input file) for the calculation of EW corrections.

Subtraction procedure and remnants
To complete the NLO EW calculation, we calculated the real photon/gluon emission processes qq → lν l γ and qq → lν l g for massive leptons l, to match the phase space generation with massive leptons described in appendix A. We checked that our calculation of the radiative processes agrees numerically with the output of MADGRAPH [47].
The POWHEG BOX implements automatically the FKS subtraction procedure, and its extension to EW processes is straightforward. As already emphasized, the routines that compute the soft and collinear limits of the real amplitude were updated in order to deal also with photons.
The soft-virtual amplitude, described in detail in section 2.4.2 of ref. [16] in the case of massless coloured partons, and in section 4.3 of ref. [22] in the case of massive partons, were extended to include also the emission of a photon by a charged particle. These formulae are easily obtained from those relative to a gluon emitted by a massive or massless quark, that we report here for completeness.
Following the same notation of formula (2.99) in ref. [22], the EW virtual part can be written as the sum of a singular and finite contribution, i.e. 3)

JHEP04(2012)037
where B ij = B for every i, j colour indices, because the colour structure is not relevant in the EW calculation, and V EW f in stands for the finite part. In eq. (3.3) Q EW and I EW depend on the masses, momenta and charges of external particles and are given by (see eq. (A.52) of ref. [22]) In eq. (3.4) q f ⊕ (q f ) denotes the charge of the ⊕ ( ) incoming particle, i is intended as a sum over the final state charged particles, q i is the charge of the i final state particle and β i ≡ | p i |/p i0 . ξ C is an arbitrary parameter which is set equal to one in the code, and µ F is the factorization scale. The factor I EW can be written as a sum of three terms, denoted as H (eq. (2.101) of ref. [16]), J (eq. (A.28) of ref. [22]) and K (eq. (A.40) of ref. [22]) related to a massless-massless, massless-massive and massive-massive charged particle pair, respectively. The H contribution is given explicitly by where the sum over i, j is a sum over all the pairs of charged massless particles. The symbol σ i is defined as in ref. [48]: σ i = +1 for incoming fermions and outgoing anti-fermions, and σ i = −1 for outgoing fermions and incoming anti-fermions. The J term reads where the sum over m, l means a sum over all the massive-massless particles pairs, (k 2 m = 0, k 2 l = 0). I 0 and I are given by eq. (A.23) and eq. (A.24) of ref. [22], respectively. Last, the K factor is given by where the sum runs over all the massive-massive particle pairs. I 0 and I are taken from eq. (A.41) and eq. (A.50) of ref. [22], respectively. For the real emission part, the subtraction formulae are the same as in QCD with the obvious substitutions α s → α em , C F → q 2 /q i q j (see eq. (2.3) of ref. [48]). The collinear remnants are the same as in eq. (2.102) of ref. [16] with the only substitutions necessary for the EW real part. Note that, since the collinear remnants contain a finite part after cancellation of all the singularities present in the NLO calculation and PDFs, this finite part is taken under control correctly only when using a PDF set accounting for both JHEP04(2012)037 QCD and QED radiation off the quarks, like e.g. MRST2004QED [49], in order to be consistent with the full NLO calculation. However, it is known from various studies that the QED contribution to the PDF evolution is actually very small for the typical x values contributing to the DY process at collider energies. As a consequence, PDF sets describing QCD radiation only can be safely used in association with complete QCD/EW collinear remnants, as done in the numerical results shown in section 5. 8 It is worth noting that the implemented subtraction procedure is not strictly related to the CC DY process, but, in principle, it is generally valid for any process involving the interaction of charged particles.

Further issues
In the present section, we briefly describe further modifications of the POWHEG BOX for an efficient and complete inclusion of EW radiation contributions.

Radiation phase space for FSR from massive partons
In principle, we could have treated the leptons as massless, using the same final state radiation setup that is used for final state massless partons in the QCD version of the POWHEG BOX. As mentioned earlier, we have instead preferred to keep as finite the lepton mass in the treatment of final state radiation from the leptons. Therefore, we had to introduce a new final state mapping for FSR, to be used for partons with small mass. This mapping augments the possibilities available in the POWHEG BOX, described in detail in section 5.1 (for ISR) and section 5.2 (for FSR) of ref. [16], and is detailed in appendix A of the present work. There were a few reasons to go along this direction. One reason is that the EW virtual corrections, as previously described, are available with a finite mass of the leptons. A second reason has to do with the fact that the mass of the leptons is the true physical cutoff, and thus by including the mass one has in fact a better result, that may also be used safely for not so light leptons.

Generation of the hardest radiation
As already pointed out in section 2, we may have up to three singular regions associated with the real graphs sharing the same underlying Born configuration. For example, the underlying Born flavour dū → eν has three α r associated with it, one for the dū → eνg real process, with an ISR region, one for dū → eνγ in the ISR region, which is singular only when the transverse momentum of the photon is small, and one for dū → eνγ, with a FSR mass singularity (in the mass of the lepton) when the photon and lepton momenta become parallel. The POWHEG Sudakov form factor for radiation of eq. (2.3) is therefore equal to the product of the Sudakov form factors for each of these regions. The POWHEG BOX handles it by generating a radiation transverse momentum with each of these Sudakov 8 Note also that, because we do not consider in our study PDF sets describing the presence of photons inside the proton, we do not include correspondingly the contribution of the so-called photon-induced processes among the real EW diagrams. However, they could be simply added in the calculation, provided an appropriate PDF set is used.

JHEP04(2012)037
form factors, and then choosing the one with the largest transverse momentum (i.e. using the so-called highest bid procedure). In all cases, a lower radiation p T -cutoff is needed, in order for the generation to terminate in finite time. This is set equal to a typical hadronic scale for gluon or photon radiation from quarks, while it is taken as the mass of the lepton for photon radiation off the leptons.

Numerical results
To validate all the new theoretical and computational ingredients of the POWHEG BOX described above, we performed a number of detailed numerical simulations of W hadroprodution at the LHC energies. More precisely 1. we compared the results of the POWHEG BOX with the predictions of the HORACE generator, in the presence of NLO EW corrections only. Note that HORACE implements completely independent EW form factors and real photon ME, computed in the mass regularization scheme. It can be considered as a benchmark, having been validated against other codes and hadron collider data.
2. We provided full results of the POWHEG BOX for NLO strong and EW corrections interfaced to MC showers, using PYTHIA and PHOTOS for QCD and QED showers, respectively. We compared these predictions for the combined QCD⊗EW corrections with those quoted in ref. [27]. Note that the results of ref. [27] were obtained according to a different methodology, by combining the events generated with HORACE (and showered with the HERWIG Parton Shower) with the QCD predictions of MC@NLO.
We considered the process pp → W + → µ + ν µ + (X), at the c.m. energy √ s = 7 TeV and using the CTEQ6L PDF set with factorization/renormalization scale µ = M W for the comparisons with the HORACE code (µ = M W being the default choice of HORACE) and µ = M lν (the lepton-neutrino invariant mass) for the combination of EW and QCD corrections. We applied the following acceptance cuts: and considered "bare" (i.e. without photon recombination) event selection conditions. The results have been obtained in the G µ scheme, using the following input parameters

Electroweak comparisons with HORACE
The comparisons between the POWHEG BOX and HORACE at the level of pure NLO EW corrections are shown in figure 1 and figure 2, for the energy around the Jacobian peak and far from it, respectively. In the upper panels, the absolute predictions of the two codes for the W transverse mass and lepton transverse momentum distributions are shown, together with the reference Born results. The NLO EW predictions of the POWHEG BOX have been obtained by setting in the code α s at a fictitious, infinitesimal numerical value. The typical some percent reduction of the peak cross section due to EW corrections is clearly seen for both the codes in figure 1, while the significant effect due to the EW Sudakov logarithms can be appreciated in figure 2. The lower panels shows the relative difference in per cent between the two generators. The error bars correspond to 1σ MC errors. It can be seen that the results of the programs agree within the statistical fluctuations, both in the peak region and well above it. This agreement represents a non-trivial check of the correctness of all the improvements and new features realized in the POWHEG BOX, from the calculation of the EW virtual contributions to the improved subtraction procedure and the new FSR mapping of collinear photon singularities. Indeed, it is known that around the W resonance the NLO EW corrections are largely dominated by final state photon radiation, containing logarithms of the type log(ŝ/m 2 l ), which are therefore correctly accounted for by the newly implemented FSR mapping and generalized subtraction scheme. On the other hand, in the high energy tails of the distributions, the NLO EW corrections display large Sudakov logarithms due to the exchange of EW massive gauge bosons in the loops and yielding corrections of tens of per cent. The agreement between the two codes shown in figure 2 demonstrates the correct treatment of the NLO EW virtual contributions in the POWHEG BOX.  We performed further numerical tests not shown here, in particular on the size and shape of the NLO EW corrections to various observables for both the muon and electron final state, to register perfect technical agreement with HORACE.

Combined electroweak and QCD corrections
The full results of the POWHEG BOX for the combined effects of QCD and EW radiation in the resonance region are shown in figure 3 for the W transverse mass, in figure 4 for the lepton p ⊥ and in figure 5 for the lepton rapidity. In figure 6 we show the transverse mass distribution above 1 TeV. The full predictions have been obtained by interfacing the NLO EW and strong corrections with QCD (PYTHIA) and QED (PHOTOS) showers. For the sake of comparison, the pure QCD predictions of the standard POWHEG BOX are also given, together with the pure NLO EW results. These absolute predictions for the various distributions are shown in the upper panels of each plot. The lower panels display the relative difference, in per cent, between the results of the new version of the POWHEG BOX and the standard QCD release, as well as the relative effect due to pure NLO EW corrections. Therefore the comparison between the two lines in each lower panel provides a measure of the combination of QCD and EW corrections and, more precisely, of mixed EW⊗QCD contributions at order α n em α n s , n ≥ 1 in perturbation theory.  From the upper panels, one can clearly see that NLO QCD corrections in association with QCD shower effects are strictly needed for a correct simulation of both the normalization and shape of the distributions. Especially for the lepton p ⊥ , the particularly pronounced smearing of the LO prediction and the rising of a heavy tail above the Jacobian peak, that are the well known effects due to QCD radiation, are clearly visible. However, also EW radiation, when combined with QCD effects, plays a role for precise calculations of the distributions, impacting both on the normalization and shape of the distributions themselves.
The latter observation can be better understood by looking in detail at the lower panels of figures 3-6. Actually, the two lines show, as already emphasized, the relative contribution due to EW radiation in association with QCD effects and the relative contribution due to NLO EW contributions only. For the W transverse mass around the peak (figure 3), the substantial agreement between the two lines indicates that the relative impact of the EW corrections is not significantly altered by the combination with QCD radiation. Stated differently, both corrections are necessary because the mainly positive and dominant QCD effects are partially compensated by the negative EW contributions, but their interplay gives rise to a relative contribution on the same ground as that provided by a pure NLO EW calculation. This means that, not surprisingly, the impact due to mixed EW⊗QCD corrections to the W transverse mass is quite moderate, with the exception of the region M W ⊥ around M W /2, as a consequence of the rather mild dependence of such a distribution on QCD effects. However, the same conclusions can not be drawn for the muon transverse momentum and rapidity distributions ( figure 4 and figure 5) in the peak region and for the transverse mass in the high tail ( figure 6). For the lepton p ⊥ , the two lines in the lower panel show a systematic, statistically significant difference, particularly around the Jacobian peak. In this region the particularly large QCD corrections to the lepton p ⊥ conspire with the per cent level EW contributions to give rise to non-negligible mixed corrections at a few per cent level. Also for the lepton rapidity, as shown in figure 5, the negative EW corrections partially reduce the positive effect of QCD radiation, yielding combined QCD⊗EW contributions of about one per cent. A similar reasoning applies to the transverse mass above 1 TeV shown in figure 6, where the quite large EW corrections, enhanced by Sudakov logarithms, in association with QCD radiation effects produce mixed contributions changing the normalization of some per cent.

JHEP04(2012)037
We compared the results discussed in this section, as well as our predictions for other observables like the W rapidity, with those derived in ref. [27] for the combination of EW and strong corrections to W production. In spite of some different details used in the simulations of ref. [27] (different PDF set, event selection cuts and input parameters) we observed a satisfactory agreement for the size and shape of the effects due to the QCD⊗EW combination. However, whereas in ref. [27] the study required, as still currently done in the simulations performed by the experimental collaborations at the Tevatron and the LHC, a rather time-consuming use of a tandem of generators (HORACE+HERWIG and MC@NLO) and the merging of their output, the results here given were directly obtained by means of a single computational framework.

Conclusions
We realized the first fully consistent implementation of EW radiation effects in the POWHEG BOX. We considered the process of single W hadroproduction as a case study and described how NLO EW corrections have been combined with the already available NLO QCD calculation. We also interfaced NLO EW and QCD corrections with QCD and QED MC showers. The resulting tool allows to study comprehensively the interplay of QCD and EW effects to W production at hadron colliders according to a unified framework, thus facilitating the work of MC simulations in data analysis.
We presented several successful cross-checks of our predictions against benchmark results to emphasize the reliability and accuracy of the new tool.
The approach followed for the implementation of the EW effects mimicked the general, process independent structure of the POWHEG BOX framework. It therefore paves the way for an almost straightforward inclusion in the POWHEG BOX of EW radiation effects in further processes, like single Z hadroprodution and other processes of phenomenological interest.
These perspectives are left to future works.
Note added in proof. During the completion of this work, an independent combination of NLO QCD and EW corrections to W hadroproduction in the POWHEG BOX framework appeared in ref. [51]. This study differs from our approach in several aspects, the most important one being the fact that multiple photon radiation in not dealt with in ref. [51]. It would be interesting to compare the results of the two implementations in the future.

A Final state radiation mapping for massive partons
In this appendix we document the construction of the final state radiation machinery for the radiation of a massless parton off a massive one. We refer throughout this appendix to the notation of ref. [16], where the Φ n refers to the Born phase space, and Φ n+1 refers to the real radiation phase space. The n and n + 1 partons are the radiating and radiated partons, respectively. The real momenta are written k 1 . . . k n+1 and the underlying Born momenta arek 1 . . .k n . We will deal here with the case when the n th parton is massive.
The main ingredient for this construction are the following: • the construction of a factorized phase space dΦ n+1 = dΦ n dΦ rad , where the full (n + 1) particle phase space is expressed in terms of a n particle phase space and a (three dimensional) radiation phase space. This factorization must be unique in both direction, i.e. given a Φ n+1 point we should get a unique (Φ n , Φ rad ) pair and viceversa.
Once the factorized phase space is known, we can compute theB function in POWHEG.
• A definition of a hardness scale for radiation K 2 ⊥ should be given, that coincides with the usual definition in the massless limit. This scale will appear in the Sudakov form factor inside a theta function θ(K 2 ⊥ − p 2 ), where p is the argument of the Sudakov form factor. It should be chosen in such a way that the integral dΦ rad θ(K 2 ⊥ − p 2 ) should be easy to do.
• A form for an upper bounding function of R/B must be given for the generation of radiation, such that the associated Sudakov form factor can be computed analytically. In order for this to be possible, this upper bounding function must be integrable over the radiation phase space with the theta function θ(K 2 ⊥ − p 2 ). Since we are seeking for an upper bound to R/B, we may extend the phase space Φ rad beyond its kinematic limits if that makes the integration easier, assuming later that R/B vanishes when the radiation variables are outside the kinematic bounds.
• The integral of the upper bounding function with θ(K 2 ⊥ − p 2 ) over the (possibly extended) radiation phase space should be performed analytically. Furthermore, once a value of p is generated using the Sudakov form factor constructed with the upper bound of R/B, the generation of Φ rad with the constraint K ⊥ = p should also be set up. Once the full Φ rad point is generated, one can use the veto technique to implement the real Sudakov form factor.

A.1 Factorized phase space
As in the massless case, the underlying Born configuration will be constructed as follows. First of all, we work in the c.m. frame of the k 1 . . . k n+1 system. We define and we will use q also to denote q 2 , when this does not generate confusion.

JHEP04(2012)037
We define the underlying Bornk n to be (spatially) parallel to k n . The underlying Born recoil systemk 1 . . .k n−1 is obtained by boosting the k 1 . . . k n−1 recoil system in a direction parallel to k n , in such a way that its three-momentum balances thek n three-momentum. The modulus of thek n three-momentum is chosen in such a way that the total energy of thek 1 . . .k n system equals the total energy of the k 1 . . . k n+1 system. Defining k rec = n−1 i=1 k i andk rec = n−1 i=1k i , we get immediatelȳ The construction of the radiation phase space is slightly more involved, and it proceeds as follows.
We begin by writing the n + 1 body phase space in factorized form of a three-body phase space involving the n th and (n + 1) th partons, and the recoil system momentum k rec , times the phase space dΦ rec of the recoil system at fixed k rec The three-body phase space part can be written in term of the Dalitz variables: so that the δ function can be integrated in d cos θ, yielding The orientation can be taken relative to any of the three bodies, so we can as well write where Ω is the direction of k rec , and φ is the azimuth of k n or k n+1 relative to k rec . The underlying Born phase space can be factorized into a two-body phase space times the phase space of the system recoiling against the emitting parton  We wish to express the phase space in terms of the underlying Born phase space and the radiation variables, that we take equal to k 0 n+1 , k 0 n and φ. In other words, we must identify eq. (A.3) with eq. (A.8) times dk 0 n+1 dk 0 dφ times a jacobian J that we should infer. This yields Since Φ rec is obtained by boosting Φ rec , the corresponding phase space elements cancel on both sides, so that, canceling other common factors, we get and dΦ n+1 = Jdk 0 n+1 dk 0 n dφdΦ n . (A.12) We thus define dΦ rad = Jdk 0 n+1 dk 0 n dφ , (A. 13) which solves the problem. It is more convenient to introduce variables that have a closer parallel to the FKS viariables. It turns out that k 0 n+1 , k 0 n and k 0 rec live in a convex Dalitz domain shown in figure 7. The boundary of that domain is defined by the configuration regions where the three vectors k n+1 , k n and k rec lie on the same line. The point at k n+1 = 0 certainly belongs to that domain. Notice that when k n+1 = 0 we have k 0 n =k 0 n and k 0 rec =k 0 rec . Since the Dalitz domain is convex, it should be possible to parametrize it as a function of two parameters, k 0 n+1 itself and z, with k 0 n =k 0 n − zk 0 n+1 . (A.14)

JHEP04(2012)037
This corresponds, in figure 7, to parametrize the point X by giving the k 0 n+1 value, and the tangent y of the angle XOA. Thus, the point k 0 n =k 0 n , k n+1 = 0 belongs to the Dalitz region for all values of z. For a given value of z, there is then a maximum value of k n+1 , such that the point is on the boundary of the Dalitz region (the point Y in figure 7). It is characterized by the condition that has to hold for at least one sign combination. Eq. (A.15) corresponds to the boundary of a triangular inequality. It is solved by squaring from which it follows again We now define k n+1 = ξq 2 , z = z 2 − (z 2 − z 1 )(y + 1)/2 , (A. 21) where now z and y will play the role of the FKS variables. Notice that y = 1 corresponds to z = z 1 , which is the lower tangent in the plot. The upper tangent corresponds to k n nearly constant, i.e. to k n recoiling against k n+1 and k rec , that are parallel. The lower tangent corresponds instead to k rec recoiling against k n+1 , k n , that are collinear. Thus, y = 1 corresponds to the mass singularity, which is what we want. Summarizing, the factorized phase space is given by dΦ n+1 = dΦdΦ rad , dΦ rad = 1 (2π) 3 The physical region in z, φ, ξ is delimited as follows: 0 < φ < 2π, z 1 < z < z 2 , with

A.2 Definition of the hardness scale
We need to find an appropriate definition for the hardness of the emission. In the massless case, the transverse momentum is used, but this is not appropriate now, because, in the soft limit, the maximum virtuality available for the emitted gluon is no longer limited by the real transverse momentum. 10 The singularity in the propagator of the massive quark has the structure 1 (k + p) 2 − m 2 , (A. 28) where we have called for ease of notation k n+1 → k, k n → p. The hardness of the emission is usually identified with the maximum virtuality that the emitted gluon can have without perturbing significantly the collinear singularity in eq. (A.28). In fact, the argument of the running strong coupling constant is set to this virtuality. If we assume that k develops a positive virtuality, with k 0 being replaced byk 0 = k 2 + k 2 0 , the variation of the denominator of eq. (A.28) is δ (k + p) 2 − m 2 = (k 0 + p 0 ) 2 − (k 0 + p 0 ) 2 = k 2 + 2p 0 (k 0 − k 0 ) ≈ k 2 + p 0 k 0 k 2 , (A. 29) where we have assumed k 2 k 2 0 . The second term on the r.h.s. of eq. (A.29) is dominant in the soft region, so we require p 0 k 0 k 2 < (k + p) 2 − m 2 , or k 2 < k 0 p 0 2k · p , (A. 30) in order for the virtuality of k not to alter significantly the collinear region. Thus, we can use as hardness definition We thus get A.5 Generation of z and ξ at fixed t We now see how z and ξ can be generated once t has been found. They are distributed according to dξdz Performing first the z integration, we get dξ q 2 t(ξq 2 − t) . (A.66) Thus we first generate ξ uniformly in log(ξq 2 − t). The extremes for ξ are given by ξ min (t) and ξ m (t) = min(ξ max , ξ 1 (t)) ξ = exp log(ξ min (t)q 2 − t) + r log ξ m (t)q 2 − t ξ min (t)q 2 − t + t /q 2 , (A.67) where 0 < r < 1 is a uniform random number. The value of z is then obtained by solving the δ function in eq. (A.65).

JHEP04(2012)037
Open Access. This article is distributed under the terms of the Creative Commons Attribution License which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.