NNLOPS accurate Drell-Yan production

We present a next-to-next-to-leading order accurate description of Drell-Yan lepton pair production processes through $\gamma^*/Z$ or $W$ exchange that includes consistently parton shower effects. Results are obtained by upgrading the vector-boson plus one jet NLO calculation in POWHEG with the MiNLO procedure and by applying an appropriate reweighting procedure making use of the DYNNLO program. We compare to existing data and to accurate resummed calculations.


Introduction
During run I at the LHC, at 7 and then 8 TeV centre of mass energy, both the ATLAS and CMS collaboration collected almost 30 fb −1 of data. Because of their very large crosssections and the very small systematic uncertainties, Drell-Yan production through W and Z exchange are standard candles at the LHC. Given the high-statistics reached, kinematic distributions can be studied over many orders of magnitude. With run II at 13-14 TeV even more W and Z events will be available, in particular it will be possible to study distributions over an even larger kinematic range. These distributions provide important input to constrain parton distribution functions. For example the rapidity distribution and the dilepton invariant mass of both neutral and charged Drell-Yan data have been used recently to constrain the photon content of the proton [1].
Higher-order corrections are indispensable for these studies, and Drell-Yan production is to date the theoretically best described process at the LHC. The cross-section is known through next-to-next-to-leading order (NNLO) in QCD including the decay of W and Z bosons to leptons [2][3][4]. Two public codes exist (DYNNLO [4] and FEWZ [5]) that implement QCD NNLO corrections to the hadronic W and Z production. Furthermore electroweak corrections have been the subject of intensive studies [6][7][8][9][10][11]. Electroweak Sudakov effects become more important at large invariant mass [9,[12][13][14], a region that was already interesting at run I and which will be explored even more during the next LHC run. Version 3 of FEWZ [15] implements also NLO EW corrections and the leading photon initiated processes. Recently, a framework for the calculation of the mixed QCD-electroweak O(α s α) corrections to Drell-Yan processes in the resonance region has been developed [16]. The impact of non-factorising (initial-final state) corrections is shown to be very small. Factorisable O(α s α ew ) corrections in the pole approximation have been computed. While QCD corrections are predominantly initial state corrections, and EW contributions are predominantly final-state corrections, because of kinematic effects there are sizable differences between the result of ref. [16] and the naive product of NLO QCD and NLO EW corrections.
While fixed-order predictions provide accurate results for inclusive distributions, there is an obvious advantage in combining state-of-the-art fixed-order perturbative calculations with parton shower generators. POWHEG [17,18] and MC@NLO [19] generators allow to keep NLO accuracy and benefit from the exclusive description from a parton shower (NLO+PS accuracy from now on). Recently, for W and Z production, QCD and EW corrections have been implemented in the POWHEG BOX framework [20][21][22].
In this paper, we take a step forward in further improving the accuracy of the Drell-Yan description by detailing and releasing an NNLO+PS accurate Monte Carlo description of this process. Our implementation relies on the following inputs: • Les Houches events for the Z+one jet or W +one jet process [23] (respectively Zj and Wj from now on), as implemented in POWHEG, upgraded with the improved MiNLO procedure of ref. [24] in such a way that NLO accuracy is guaranteed for inclusive distributions without any jet cut; • NNLO accurate distributions, differential in the Born kinematics of the leptons, as obtained from DYNNLO [4]; • a local reweighting procedure, described in all details in Sec. 2.
Compared to the standard NLOPS implementation of W and Z in POWHEG we find a considerably reduced theoretical uncertainty for inclusive distributions. This is expected since our predictions have full NNLO accuracy. Furthermore the 1-jet region is described at NLO accuracy in our framework, even at small transverse momentum.
Throughout this work we pay particular attention to the issue of assigning a theory uncertainty to our predictions. We describe our procedure in Sec. 3.1. Unlike the standard POWHEG approach (without a separation between singular and finite real contributions), that is known to underestimate the theoretical uncertainty for the W , Z or Higgs boson transverse momentum [25], we believe that our uncertainties are more reliable. Compared to pure NNLO predictions, we find in general a better description of observables sensitive to multiple emissions, such as the boson transverse momentum, φ * , and jet-resolution variables d i (which just vanish at NNLO starting from i = 2). This is both because of the underlying MiNLO procedure, and because of the POWHEG framework.
An important validation of our results comes from comparing with precise calculations of specific observables. For the Z transverse momentum distribution we find a good inclusive in all QCD radiation. After integration over all QCD radiation only the leptonic system is left. This can be characterised by three independent variables, for instance one can choose the invariant mass of the lepton pair, m ll , the rapidity of the boson before decay y V , and the angle of the negatively (or positively) charged lepton with respect to the beam, θ l . Other choices are possible, and later on we will discuss the impact of different choices. In the following we will use Φ B to denote the ensemble of these three Born variables.
We denote schematically the fixed order NNLO cross-section differential over Φ B by dσ NNLO /dΦ B and the cross-section obtained from Vj-MiNLO by dσ MINLO /dΦ B . Since these distributions are identical up to terms of O (α S ) terms, their ratio is equal to one up to O α 2 S terms: where the c i are O (1) coefficients. Since the Vj-MiNLO generator reproduces the inclusive fixed-order result up to and including NLO terms, the NLO accuracy of the cross-section in the presence of one jet (that starts at order α S ) is maintained if the cross-section is reweighed by the factor in eq. (2.1). This is follows from the simple fact that the reweighting factor combined with this cross-section yields spurious terms of order α 3 S and higher. It is also obvious that by reweighting Vj-MiNLO distributions with this ratio, any of the three Φ B distributions acquires NNLO accuracy, and in fact coincides with the NNLO distribution. We will now argue that the Vj-MiNLO generator reweighed with the procedure of eq. (2.1) maintains the original NNLO accuracy of the fixed-order program, used to obtain the dσ NNLO /dΦ B distribution, for all observables. The proof trivially extends the proof given in ref. [31] by replacing the Higgs rapidity with the chosen set of three variables Φ B associated to the Born phase space. We find it useful to recall the way the proof works.
The claim is based on the following theorem, that we will prove in the following.
A parton level V boson production generator that is accurate at O(α 2 S ) for all IR safe observables that vanish with the transverse momenta of all light partons, and that also reaches O(α 2 S ) accuracy for the three Born variables Φ B achieves the same level of precision for all IR safe observables, i.e. it is fully NNLO accurate.
To this end, we consider a generic observable F , including cuts, that is an infrared safe function of the final state kinematics. F could be for instance be a bin of some distribution. Its value will be given by with a sum over final state multiplicities being implicit in the phase space integral. Infrared safety ensures that F has a smooth limit when the transverse momenta of the light partons vanish.
Since the Born kinematics is fully specified by the Born kinetic variables Φ B , such a limit may only depend upon the value of Born kinematic variables Φ B . We generically denote such a limit by F Φ B . The value of F can be considered as the sum of two terms: tends to zero with the transverse momenta of all the light partons, by hypothesis its value is given with O(α 2 S ) accuracy by the parton level generator. On the other hand, which is also exact at the O(α 2 S ) level by hypothesis. Thus, level. The Vj-MiNLO parton level generator (in fact even just Vj) fulfills the first condition of the theorem since it predicts any IR safe observable that vanish when the transverse momentum of the light partons vanish with O(α 2 S ) accuracy. The second hypothesis of the theorem, regarding NNLO accuracy of the Born variables, is simply realised by augmenting the Vj-MiNLO generator by the reweighting procedure described above. We note that MiNLO is crucial to preserve the first property after rescaling. The proof of O(α 2 S ) accuracy for these observables thus corresponds to the general proof of NLO accuracy of the POWHEG procedure, given in refs. [17,32].
Observe, also, that for observables of the type F − F Φ B , adding the full shower development does not alter the O(α 2 S ) accuracy of the algorithm, for the same reasons as in the case of the regular POWHEG method. The only remaining worry one can have concerns the possibility that the inclusive Born variable distributions are modified by the parton shower evolution at the level of O(α 2 S ) terms. However, our algorithm already controls the two hardest emissions with the required α 2 S accuracy. A further emission from the shower is thus bound to lead to corrections of higher order in α S . This concludes then our proof.

Variant schemes
As discussed in detail in ref. [31] the reweighting in eq. (2.1) treats low and high transverse momentum distributions equally, i.e. it spans the virtual correction over the full transverse momentum range considered. On the other hand, if one considers the transverse momentum distribution of the vector boson V or of the leading jet, in the high transverse-momentum region Vj-MiNLO and the NNLO calculation have formally the same (NLO) accuracy. Hence, while it is not wrong to do so, there is no need to "correct" the Vj-MiNLO result in that region. It is therefore natural to introduce a function that determines how the two-loop virtual correction is distributed over the whole transverse momentum region, where β and γ are constant parameters, and m V is the mass of the vector boson. This function has the property that when p T goes to zero it tends to one, while when p T becomes very large, it vanishes.
The function h(p T ) can therefore be used to split the cross-section according to We then reweight the Vj-MiNLO prediction using the following factor which preserves the exact value of the NNLO cross-section (2.9) The proof of NNLO accuracy for this rescaling scheme is completely analogous to the proof given above, so we omit it here.

Phenomenological analysis
In this section we will present some validation plots and comparisons with resummed predictions. Before doing so, we will define the procedure we use to estimate theoretical uncertainties and we will give all details about our practical implementation of the formulae in Sec. 2.

Estimating uncertainties
We detail here the method that we use to estimate the uncertainties in our NNLO event generator. As is standard, the uncertainties in the Vj-MiNLO generator are obtained by varying by a factor 2 up and down independently all renormalisation scales appearing in the MiNLO procedure by K R (simultaneously) and the factorisation scale by K F , keeping 1/2 ≤ K R /K F ≤ 2. This leads to 7 different scale choices given by (3.1) We will consider the variation in our results induced by the above procedure. The seven scale variation combinations have been obtained by using the reweighting feature of the POWHEG BOX.
For the pure NNLO results, the uncertainty band is the envelope of the 7-scale variation obtained by varying the renormalisation and factorization scale by a factor 2 around the central value M V keeping 1/2 ≤ K R /K F ≤ 2. For the NNLOPS results, we have first generated a single of Zj-MiNLO event file with all the weights needed to compute the integrals dσ MINLO A/B /dΦ B , in eq. (2.8), for all 7 scale choices. The differential cross-section dσ NNLO /dΦ was tabulated for each of the three scale variation points corresponding to K R = K F . The analysis is then performed by processing the MiNLO event for given values of (K R , K F ), and multiplying its weight with the factor (3. 2) The central value is obtained by setting (K R , K F ) and (K R , K F ) equal to one, while to obtain the uncertainty band we apply this formula for all the seven (K R , K F ) and three (K R , K F ) choices. This yields 21 variations at NNLOPS level. 3 The reasoning behind varying scales in the NNLO and Zj-MiNLO results independently is that we regard uncertainties in the overall normalisation of distributions, as being independent of the respective uncertainties in the shapes. This is consistent with the recently introduced efficiency method [33], used for estimating errors on cross-sections in the presence of cuts.

Practical implementation
There is some degree of freedom in the way the reweighting procedure described above is carried out in practice. In particular, one is free to choose the three Born variables with respect to which one performs the reweighting, as well as the form of the damping factor h in eq. (2.4). Our choice of the Born variables is driven by the fact that one wants to populate all bins in the three-dimensional histograms sufficiently well. To produce the results presented in the following we used the rapidity of the Z boson, y Z , a variable directly related to the dilepton invariant mass, a mll = arctan((m 2 ll − M 2 Z )/(M Z Γ Z )) and θ * l , where the latter is the angle between the beam and a charged lepton in the frame where the boson has no longitudinal momentum. 4 In the following we will use the values of β = 1 and γ = 2 in eq. (2.4). The choice for β is motivated by the fact that typical resummation scale for Drell-Yan production is set to the boson mass. The second choice is partly driven by the fact that in the case of Higgs production, with this choice the NNLOPS Monte Carlo agrees well with the Higgs transverse momentum and leading jet NNLL+NNLO resumed results. We stress however that while β and γ are arbitrary, the dependence on β and γ is formally O(α 3 S ) (or exactly zero in the case of inclusive quantities). p T in eq. (2.4) denotes the transverse momentum that is used to decide how to distribute the virtual corrections. One could for instance choose the V transverse momentum, the leading jet transverse momentum, or the total transverse momentum of the event. In the following, we will adopt the choice of ref. [31], namely to use the leading jet transverse momentum when clustering events according to the inclusive k T -algorithm with R = 0.7 [34]. This choice ensures that h goes to one when no radiation is present, since in that case the leading jet transverse momentum vanishes. On the contrary, the V transverse momentum can vanish also in the presence of radiation.
Finally we note that for the reweighting we used 25 bins per variable, meaning that our 3-dimensional distributions involve a total of 15625 equal bins. Results presented in the following are based on generating 20 Million events. Even if we perform high-statistics runs, there might be bins that are not well-populated. This can give rise to a rescaling factor that is unphysically large just because too few events ended up in one bin, and this in turn can give rise to spikes in kinematical distributions. To avoid these occurrences, instead of using the local reweighting factor, we use the global K NNLO /K NLO factor to perform the reweighting whenever the local reweighting factor (which is formally 1 + O(α s ) 2 ) exceeds 5. This happens rarely, in about 0.3% of the points. We checked that this procedure has no visible systematic effect on distributions, other than that of removing unphysical spikes. Other workarounds could of course also be adopted.
Before showing validation plots, we list here the settings used for the results obtained in this paper. We used the code DYNNLO [4] to obtain NNLO predictions. 5 Throughout this work we consider the MSTW2008NNLO parton distribution functions [35] and set M Z = 91.1876 GeV, Γ Z = 2.49595 GeV, M W = 80.398 GeV and Γ W = 2.008872 GeV. We choose to use α em = 1/128.94 and sin 2 θ W = 0.22264585. Jets have been constructed using FastJet [36,37]. To compute the h(p T ) factor we use the k T algorithm [34,38] with R = 0.7. To shower partonic events, we have used both Pythia8 [39] (version 8.185) with the "Monash 2013" [40] tune and Pythia6 (version 6.4.28) with the MSTW2008LO variation of the "Perugia" tune ("Perugia P12-M8LO", tune 378). We have generated Z and W events with decays into electrons and positrons, and always switched off QED radiation off leptons and quarks in the showering stage. Moreover, in the following, in order to define the leptons from the boson decays we will always use the Monte Carlo truth, i.e. we disregard complications due to the fact that there might be other leptons in the event.
To obtain the results shown in the following, we have switched on the "doublefsr" option introduced in ref. [41] and used a standard driver for Pythia6 and the Pythia8 driver suggested by the Pythia authors for showering events generated with POWHEG. Although we have also explored the effect of using the alternative prescription first introduced in section 4 of ref. [41] to compute the scale used by the parton shower to veto hard emissions (which is also available as an option in the Pythia8 driver), 6 the plots shown throughout the paper have been obtained keeping the veto scale equal to the default POWHEG prescription, both for Pythia6 and Pythia8.

Validation plots
In this section we consider the case of inclusive Z → e + e − production at 14 TeV, with no cuts on the final state other than requiring the dilepton invariant mass to be in the range 66 GeV ≤ m ll ≤ 116 GeV. In order to validate our results, we will show comparisons 5 Even with high-statistics runs with FEWZ [5], we did not obtain high-quality triple differential distributions, as required here. 6 A factor 2 is missing in eq. (4) of ref. [41], but not in the practical implementation.
to NNLO predictions and to MiNLO-improved Zj-POWHEG results (henceforth denoted as Zj-MiNLO). Since we compare to DYNNLO, for the Zj-MiNLO and NNLOPS results it is useful to consider here pure parton level results, before hadronization (with underlying event and multiparton interaction switched off). Unless otherwise stated, for DYNNLO we have set µ F = µ R = M Z , and the associated uncertainty bands are obtained from a 7-points scale variation envelope. Also unless stated otherwise, we will shower events with Pythia8.
In Fig. 1 we show the Z boson rapidity distribution y Z as predicted at NNLO (green),  with Zj-MiNLO (blue), and at NNLOPS (red). As expected NNLOPS agrees very well with DYNNLO over the whole rapidity range, both for the central value and the uncertainty band, defined as detailed in Sec. 3.1. We also note that as expected the uncertainty band of the NNLOPS result is considerably reduced compared to the one of Zj-MiNLO, which is NLO accurate. In the central region the uncertainty decreases from about (+5:−15) % to about (+2:−2) %. We finally note that because of the positive NNLO corrections, the central value of NNLOPS lies about 5 % above the one of Zj-MiNLO, while no considerably difference in shape is observed in the central rapidity region. Moderate but slightly more pronounced shape differences can be seen at large rapidity.
We proceed by examining the other two distributions that have been used in the reweighting procedure, namely a m ll and θ  invariant mass of the dilepton system m ll which is directly related to a m ll . We notice that the same features observed above hold: NNLOPS agrees well with DYNNLO, it tends to be about 5−10% higher than Zj-MiNLO and the uncertainty band is reduced by about a factor 4. No sizable difference in shape is observed in these two distributions, when comparing Zj-MiNLO and NNLOPS. We now show in Fig. 3 the Z boson transverse momentum in two different ranges. At  At small transverse momenta, while DYNNLO diverges, NNLOPS remains finite because of the Sudakov damping. The difference in shape observed between DYNNLO and NNLOPS at finite p T,Z has to do with the fact that in that region the fixed-order calculation has to compensate for the divergent behaviour at small p T,Z . We also note that the uncertainty band in DYNNLO is far too small when approaching the divergence at p T,Z = 0. The uncertainty band of NNLOPS instead tends to increase at very low transverse momenta, reflecting the fact that one is approaching a non-perturbative region. One can however also note that the uncertainty tends to shrink at about p T,Z = 10 GeV. We have checked that this is not an artifact due to having used a 21-point scale variation as opposed to a 49-point one. We attribute this feature to the fact that the uncertainty band of the fixed-order result shrinks in this region. This is true both for the 3-point and the 7-point scale variation in the fixed order, although in the latter case this effect is slightly less pronounced. We also observe that our Zj-MiNLO result that uses the MiNLO scale prescription does not show this feature. When we upgrade Zj-MiNLO to NNLO accuracy, we necessarily inherit this feature from the NNLO results we are using as input. It is also worth mentioning that this feature has been already observed in several studies where an analytic resummation matched with fixed-order results was performed for this observable [26,42,43].
We end our discussion on neutral Drell-Yan production by looking briefly into the effects of including non-perturbative contributions in the NNLOPS simulation, by turning on hadronization, underlying event and multiparton interaction (MPI). In particular, given the small perturbative uncertainties found with p T,Z , it is interesting to see how much non-perturbative corrections affect the 5-15 GeV region.
In Fig. 4 we show p T,Z (left) and the leading-jet transverse momentum, defined according  to the anti-k T [44] algorithm with R = 0.7 (right), after all non-perturbative stages are included, with Pythia6 (blue) and Pythia8 (red). We observe sizable differences between the two results, in particular for the Z boson transverse momentum at p T,Z < 15 GeV. This is not surprising since this is a region dominated by soft effects, hence the details of the modelling of non-perturbative effects are expected to matter. For the jet-transverse momentum the difference between the two shower models is slightly smaller. This can probably be attributed to the fact that the Sudakov peak is at larger values of the transverse momentum, compared to the Z-boson transverse momentum. Finally, for illustrative purposes, we conclude this section by showing in Fig. 5 predictions for the transverse momentum (left) and pseudorapidity (right) of the charged electron or positron in W production at 7 TeV (combining the W + and W − samples). It is interesting to look at these leptonic observables since they don't coincide with the quantities we are using to perform the NNLO reweighting. Nevertheless, we should recover NNLO accuracy in the regions where the lepton kinematics probes the fully-inclusive phase space: this is precisely what we have found, as illustrated in Fig. 5. Here we include no cuts on the final state other than requiring the transverse W mass m T,W = 2(p T,l p T,miss − p T,l · p T,miss ) to be larger than 40 GeV. All parameters and settings are the same as those used for the Z production case and DYNNLO results have been obtained choosing µ R = µ F = M W as central scale and performing a 3-point scale variation. As expected, in the left panel of Fig. 5 we observe a much more narrow uncertainty on the charged lepton transverse momentum for values of p T,l smaller than M W /2, and a very good agreement with DYNNLO in this region, both in the absolute value of the crosssection as well as in the size of the theoretical uncertainty band. When p T,l is larger than M W /2 all distributions have NLO accuracy only, since in this region the lepton kinematics requires non-vanishing values for p T,W . We observe that in fact all uncertainty bands are larger in this region. Our NNLOPS result reproduces the Zj-MiNLO one well, while there is some difference between DYNNLO and Zj-MiNLO predictions. This is expected since the scales used in the two calculations are effectively different: DYNNLO always uses the mass of the W boson, whereas in MiNLO the transverse momentum of the W boson is used. Therefore when the lepton is just slightly harder than M W /2 we are probing phase-space regions where the bulk of the cross-section typically has 0 p T,W M W , and hence DYNNLO yields smaller cross-sections.
It is also worth mentioning that both the NNLOPS and Zj-MiNLO plots exhibit a smooth behaviour in proximity of the jacobian peak p T,l M W /2, also when thinner bins (not shown) are used. This smooth behaviour is due both to parton-shower effects and to the MiNLO Sudakov form factor. On the contrary, the NNLO prediction displays the typical numerical instability of fixed-order predictions due to the numerical cancellation between real and virtual corrections close to the kinematical boundary. Furthermore, the NNLO prediction has spuriously small uncertainties in this region.
Finally, in the right panel of Fig. 5 we plot the rapidity distribution of the charged lepton. Since in each bin of this distribution we are fully inclusive with respect to QCD radiation, we observe the expected good agreement with the NNLO prediction over the whole range, as well as a quite narrow uncertainty band.

Comparison to analytic resummations
The MiNLO method at the core of the results presented in this paper works by including NLL and (some) NNLL terms in the Sudakov form factors used to improve the validity and accuracy of the underlying NLO computation at hand. Although the formal logarithmic accuracy achieved by MiNLO-improved POWHEG simulations has not been addressed (yet), it is interesting to compare NNLOPS predictions against results obtained with higher-order analytic resummation, for observables where the latter are available. In this subsection we will show results for Z production and focus on three quantities for which NNLL resummation has been performed.
The classical observable to consider in Drell-Yan production to study the effects of softcollinear radiation is the transverse momentum of the dilepton-pair system. This observable has been extensively studied in the past and is now known to NNLO+NNLL level [26,42]. In Fig. 6 we show a comparison between the NNLO+NNLL resummation obtained using DYqT and our NNLOPS result obtained with Pythia8 (left) and Pythia6 (right), switching off all non-perturbative effects (i.e. hadronization and MPI are switched off, primordialk T is set to zero). As usual, the uncertainty band for our results has been obtained as the envelope of a 21-pts scale variation. DYqT uses as resummation scale m ll , and the associated band has been obtained varying µ R and µ F among the usual 7 combinations. On top of this, for the central value of µ R = µ F we varied the resummation scale by a factor 2 up and down. This gives 9 combinations, the envelope of which is used here to define the DYqT uncertainty. The ratio plot shows a pattern quite similar to what was observed in Figs. 4 and 5 of ref. [31], namely differences of up to O(10 − 12%) between analytic resummation and NNLOPS, with a slightly more marked difference in the very first bin. In the region p T,Z ∼ 5 − 15 GeV, the uncertainty bands do not overlap, mainly because of the very narrow uncertainty bands in both predictions, in particular in the NNLO+NNLL result. In the case of Higgs production, instead, uncertainty bands are wider, hence the predictions are more compatible. Changing the β parameter might improve this agreement, although we recall that the NNLOPS prediction does not have NNLL accuracy in this region. By comparing the two NNLOPS results shown in the two panels of fig. 6, we also observe that the spectra obtained with Pythia8 are typically ∼ 5 % harder than those with Pythia6, a feature that was already noticeable in Fig. 4, and which will be present also in other distributions where NNLOPS results are "only" NLO accurate. Few percent differences between different NLO+PS results in these kinematic regions can be due to subleading effects, such as differences in details of the two parton-shower algorithms, as well as the use of different tunes. A comprehensive assessment of these issues goes beyond the purpose of this work, and is therefore left for future work. Another interesting observable to consider is the φ * distribution which is a measure of angular correlations in Drell-Yan lepton pairs [45]. This observable is defined as [43] φ * = tan π − ∆φ 2 sin θ * , where ∆φ is the azimuthal angle between the two leptons and θ * is the scattering angle of the electron with respect to the beam, as computed in the boosted frame where the Z boson is at rest. We note that ATLAS uses a slightly different definition of the angle θ * , and defines it as Since we will compare to ATLAS data in Sec. 4.2, we will use the latter definition throughout this work. In Fig. 7 we compare our NNLOPS simulation with the NLO+NNLL resummation of ref. [43]. 7 From the definition of φ * , it is clear that large values of φ * correspond to events where the Z boson tends to be boosted, while for low values the Z boson is almost at rest. We see that the two predictions agree reasonably well for φ * 0.2, in particular when Pythia6 is used, while the uncertainty bands to not overlap below that point. 8 For high φ * , Pythia8 is slightly harder than Pythia6. Since for large values of φ * the probed phase space regions are dominated by large values of p T,Z , this difference is expected, in view of the discussion at the end of the previous paragraph. We will show a comparison to data for this observable in Sec. 4.1, where we will also comment on the impact of non-perturbative corrections.
Finally we consider the jet-veto efficiency which is defined as T,veto 0 dp T,j 1 dσ(p T,j 1 ) dp T,j 1 . (3.5) Here results have been obtained for the LHC at 8 TeV, using the anti-k t algorithm to construct jets. We restrict ourselves to show NNLOPS results obtained with Pythia8. Analytically resummed results have been obtained using JetVHeto [28], and have NNLL+NNLO accuracy. The JetVHeto results have been obtained using its default setting, i.e. using as a central value for the renormalisation, factorisation and resummation scale M Z /2. 9 As is recommended in ref. [28], for the JetVHeto results we have obtained the uncertainty band as an envelope of eleven curves: we have varied the renormalization and factorization scales independently giving rise to the usual 7-scale choices, additionally, for We observe a very good agreement between the two approaches for R = 0.4 and 0.5, whereas for R = 1 differences are more marked. Few comments are in order here: first, the pattern shown in the plots is consistent with what was already observed in the Higgs case ( fig. 7 of ref. [31]), namely differences up to few percents, and good band overlap, for smaller values of R, and larger differences for R = 1. For very large values of R, the leadingjet momentum will balance against the transverse momentum of the vector boson. Given what we observed for p T,Z , it is therefore no surprise that, when R = 1, we have O(3 − 5)% differences with respect to the resummed result for values of p T,veto ∼ 25 − 30 GeV, as used currently by ATLAS and CMS in Higgs production.

Comparison to data
In this section we compare our predictions with a number of available data from ATLAS, both for Z and for W production.

Z production
We show here a comparison to a number of measurements performed by ATLAS at 7 TeV [46][47][48]. ATLAS applied the following cuts: they consider the leptonic decay of the Z boson to electrons or muons and require an electron (muon) and a positron (anti-muon) with p T > 20 GeV and rapidity |y| < 2.4. The invariant mass of the di-lepton pair should lie in the window 66 GeV < m ll < 116 GeV.
We begin by showing in Fig. 9 a comparison of NNLOPS results (with two versions of  Figure 9. Comparison of NNLOPS prediction obtained with Pythia8 (left) and Pythia6 (right) to data (black) from ref. [46] for the Z boson rapidity distribution at 7 TeV LHC.
Pythia) to data from ref. [46] for the Z boson rapidity distribution. As expected, our result displays a quite narrow uncertainty band, due to having included NNLO corrections. Since this is a fully inclusive observable, the absolute value of the cross-section and the size of uncertainty band will be driven by the NNLO reweighting: hence, Pythia6 and Pythia8 results are almost indistinguishable, as expected. We also observe that we agree with data within the errors for central rapidities. At high rapidity, however, there seems to be a tension between data and our results. This discrepancy between data and pure NNLO was already observed in the original ATLAS paper, although the NNLO results shown in ref. [46] have a slightly larger uncertainty band since they also contain PDF uncertainties. We note that, at the moment, the dominant error is coming from data. We therefore expect the agreement to improve, as more data become available, although systematic errors are non-negligible [46]. In Fig. 10 we now show the same comparison for the Z boson transverse momentum  LHC 7 TeV Figure 11. As in previous figure, but with more luminosity, thinner binnings, and up to larger values of p T,Z . Data are now from taken from ref. [49].
against data from ref. [47]. In the left panel we use Pythia8 to shower events and in the right panel Pythia6. We see that there is a very good agreement between our NNLOPS prediction over the whole p T range. Only in the very first bin we observe a slight tension with data in the case of Pythia8, which can be due to a different modelling of non-perturbative effects. Given the fact that we use two different parton-showers, with different tunes (and even different PDFs for the tunes), an O(20%) difference in the 1 − 3 GeV region can be expected. It is foreseeable that once a tuning of Pythia8 is done in conjunction with POWHEG, such differences will go away. Fig. 11 shows again p T,Z , based now on 4.7 fb −1 of data from ATLAS [49]. Due to the thinner binning, it is now possible to appreciate clearly the differences between Pythia6 and Pythia8: the NNLOPS result obtained with Pythia6 shows a remarkable agreement with data across all the p T,Z range, whereas with Pythia8 we can observe differences of up to 10% for 5 p T,Z 100 GeV. What we observe from these plots is consistent with fig. 10, the difference being that here the improved precision in data allows to conclude that, for the setups and tunes we are using, the best description is obtained when using Pythia6.
Next, in Fig. 12, we consider the comparison to data for the φ * distribution. Although  Fig. 10 for the φ * distribution in Z boson production. Data taken from ref. [48].
both our predictions are consistent with the ATLAS measure, it is clear that events showered with Pythia6 agree better with data, whereas NNLOPS showered with Pythia8 has a slightly different shape, exhibiting a dip in the theory prediction compared to data, at around φ * = 0.06. By comparing with our predictions before the inclusion of non-perturbative effects, we have checked that these effects play a sizable role in the region below φ * < 0.1, therefore it is not completely unexpected that Pythia6 and Pythia8 give slightly different shapes. The fact that non-perturbative effects introduce non-trivial changes on the shape of this distribution was also noted in ref. [27] (for predictions at the Tevatron). We also observe the same pattern shown in ref. [27], namely a moderate increase at φ * ∼ 0.05 and a more pronounced decrease for very low values of φ * when non-perturbative effects are included. Finally, it is worth mentioning that the results we have obtained (especially with Pythia6) clearly show a better agreement with data than what was found in ref. [47], where different tunes were used, both for Pythia6 and Pythia8. It is difficult to draw a solid conclusion, since in ref. [47] POWHEG-Z (as opposed to Zj-MiNLO) was used, and events were also reweighted using ResBos [50]. Nevertheless it seems clear that the best predictions are obtained when higher-order perturbative corrections are included and modern tunes are used.

W production
In this section we compare our predictions to results of refs. [51,52], and in particular we use the combined decay of the W to electrons and muons. The charged lepton is required to have p T > 20 GeV and rapidity |y| < 2.4. The event must have a missing energy p T,miss > 25 GeV and the transverse mass of the W boson defined as m T,W = 2(p T,l p T,miss − p T,l · p T,miss ) must be above 40 GeV. As was the case for the Z studies, ATLAS only provides normalised distributions for the standard observables we will show here. We start by showing in Fig. 13 the transverse momentum of the W boson p T,W . As expected, differences in the parton shower algorithm only play a visible role in the small p T region, where minor differences between Pythia6 and Pythia8 can be observed. For high p T,W , the two predictions are consistent with each other, and agree quite well with data.  Figure 13. Comparison of NNLOPS prediction (red) to data (black) from ref. [51] for the W boson p T,W distribution, using Pythia8 (left) and Pythia6 (right).
We also compared our predictions with the analysis performed by the ATLAS collaboration in ref. [52]. We show results for k T -splitting scales in W +jets events. These observables are defined as the smallest distances found by the (inclusive) k T -algorithm at each step in the clustering sequence. The splitting scale d k is the smallest among all distances found by the algorithm when going from (k + 1) to k objects. Therefore √ d 0 corresponds to the transverse momentum of the leading jet, whereas √ d 1 is the smallest distance among pseudo-jets when clustering from 2 to 1 jet: where are the usual distances used in the k T -algorithm. Among other reasons, these observables are interesting because they can be used as a probe of the details of matching and merging schemes. Due to the underlying Zj-MiNLO simulation, our NNLOPS prediction is NLO accurate for large values of p T,j 1 , and it is at least LL accurate in describing the 1 → 0 jet transition, which is measured in the d 0 distribution. The second jet spectrum and the 2 → 1 jet transition (which is encapsulated in d 1 ) are instead described at LO+LL, due to the underlying POWHEG simulation. Since the definition of d 1 contains d 12 , this observable is a measure of the internal structure of the first jet, and not only of the second jet transverse momentum. In Figs. 14 and 15 we show our NNLOPS predictions against ATLAS data, using as jet radius R = 0.6. We find good agreement, especially when √ d i > 10 GeV. Below this value, we are still compatible with the experimental uncertainty bands, although we are systematically lower than data. Once more, one should consider that the region below 5 − 10 GeV will be affected also by non-perturbative effects. For large values of d i we are instead sensitive to the level of accuracy that we reach in describing hard emissions. In this respect, it is no surprise that we have a better agreement with data than the POWHEG results shown in ref. [52], where d 1 is poorly described since the second emission is only described in the shower approximation. NLO corrections to the W + 1 jet region are included in the NNLOPS simulation, and are very likely the reason why we have a description of d 0 that is better than what was observed in ref. [52].
Finally in Fig. 16 we show the distribution for the ratio d 1 /d 0 , for events with √ d 0 > 20 GeV. Due to the ratio nature of this quantity, a simultaneous over-or underestimation   Figure 16. Comparison to 7 TeV LHC data from ref. [52] for the W boson ratio of the k T splitting scales √ d 0 and √ d 1 using Pythia8 (left) and Pythia6 (right).
in predicting d 1 and d 0 should be partially compensated when plotting d 1 /d 0 . It is therefore no surprise that the agreement with data is better than in Figs. 14 and 15.

W and Z polarization
Recently both ATLAS [53] and CMS [54] have published results on the polarization of the W boson at 7 TeV confirming the Standard Model prediction, that W bosons are mostly left-handed in pp collisions at large transverse momenta [55]. Knowledge about the W boson polarization is important, as it provides a discriminant in searches for new physics. We first very briefly review how to measure the polarization in terms of angular coefficients but refer the reader to the literature for a complete description of the topic [55][56][57][58][59][60][61][62].
Here we will follow the derivation of [55]. We then continue to compare ATLAS data [53] to our NNLOPS results for the W boson polarization and present predictions for the angular coefficients for the Z boson at 8 TeV.
The vector boson cross-section can be expanded in terms of cos θ * and φ * as where the polar angle θ * and the azimuthal angle φ * are defined in some particular rest frame of the dilepton system. Here we will make use of two different frames, the Collins-Soper frame 10 [56] for the Z boson and the helicity frame for the W boson defined as the dilepton rest frame with the z-axis pointing along the direction of flight of the W boson in the lab frame. The cross-section can be differential in any quantity that does not depend on the individual lepton kinematics. The angular coefficients can then be expressed in terms of expectation values defined as by A 0 = 4− < 10 cos 2 θ * >, A 1 = < 5 sin 2θ * cos φ * >, A 2 = < 10 sin 2 θ * cos 2φ * > , A 3 =< 4 sin θ * cos φ * >, A 4 = < 4 cos θ * >, A 5 = < 5 sin 2 θ * sin 2φ * > , A 6 =< 5 sin 2θ * sin φ * >, A 7 = < 4 sin θ * sin φ * > . For both W and Z production it is known that at O(α s ) A 5 = A 6 = A 7 = 0 and that the higher-order corrections are small. We have checked that with 20 million events processed the coefficients do not deviate significantly from zero and we have therefore chosen not to show them here. It is also interesting to notice that as a consequence of the spin-1 structure of the gluon, the Lam-Tung relation A 0 = A 2 11 holds at LO [63]. As we will see, the deviations from this lowest-order result can become quite large, O(20%), in the p T,Z > 10 GeV region. From the angular coefficients we then define the left, f L , right, f R , and longitudinal, f 0 , polarization fractions of the W ± as 10 One defines in the laboratory frame the two beam directions by b+ = (0, 0, 1; 1) and b− = (0, 0, −1; 1).
After boosting to the dilepton centre of mass frame, O , one defines the z-axis as the bisector of b + and − b − such that the z-axis points into the hemisphere of the Z boson direction (in the lab frame). One then defines a q-axis lying in the plane spanned by the b + and b − vectors, orthogonal to the z-axis and pointing in the direction opposite to b + + b − . θ * is now defined with respect to the z-axis and φ * with respect to the q-axis. 11 For qq initiated production the relation is exact to all orders. For qg initiated production it is violated at the NLO level [55]. This is independent of the frame.
It is clear from these equations that the polarization fractions are normalised such that f L + f R + f 0 = 1 and it is therefore sufficient to show results for f L − f R and f 0 . In this way we also separate the A 0 and A 4 dependence. In Table 1 we show our prediction for f L − f R and f 0 for combined W + and W − production at 7 TeV compared to ATLAS data [53] for two different p T regions. We find good agreement between the data and predictions -three of them at the 1σ level and at the 2σ level for f L − f R in the low p T range. It should be noted that the measurements are dominated by systematics uncertainties. 35 GeV < p T,W < 50 GeV p T,W > 50 GeV There is currently no public measurement for the angular coefficients for the Z boson in pp collisions. Previously measurements of A 0 , A 2 , A 3 and A 4 in pp collisions have been published by the CDF collaboration [64] and anticipating that such an analysis might be carried out at the LHC at 8 TeV we present here our predictions for the angular coefficients as defined in eq. (4.5), as a function of the boson transverse momentum. In the left panel of Fig. 17 we plot A 0 (red) and A 2 (blue) along with the difference A 0 − A 2 , while in the right panel we plot A 1 , A 3 and A 4 . For the two coefficients A 0 and A 2 we find that  Figure 17. Predictions for angular coefficients eq. (4.5) in the Collins-Soper frame for Z production at 8 TeV. In the left panel we show A 0 (red) and A 2 (blue) along with the difference A 0 − A 2 (green) and in the right panel we plot A 1 (red), A 3 (blue) and A 4 (green).
they are in very good agreement for low transverse momentum, p T,Z < 10 GeV. For higher transverse momenta deviations start to appear peaking around 60 GeV. For the three other coefficients, we observe that the p T dependence is strongest for larger transverse momenta. We can define polarization fractions for the Z boson analogously to those defined in eq. (4.6) by [65] f L = 1 4 where α = , c L is the coupling to the left-handed lepton and c R is the coupling to the right-handed lepton. As for the W boson these are normalised such that f L +f R +f 0 = 1 and we have in this case that f L − f R = − 1 2 αA 4 . Using c L = sin 2 θ W − 1 2 and c R = sin 2 θ W we get α ≈ 0.032. Comparing with Fig. 17 we see that at low p T the longitudinal polarization of the Z boson is highly suppressed and it is produced almost equally between the left-and right-handed polarization. For high p T the longitudinal polarization starts to dominate.

Conclusions
In this paper we presented an implementation for Z → l + l − and W ± → l ± ν production at NNLO level including parton shower effects. At the core of this method lies the fact that Zj-MiNLO and Wj-MiNLO simulations achieve NLO accuracy also for fully inclusive distributions, i.e. once the jet is integrated out. In the case of Higgs production considered recently in ref. [31], it was enough to rescale the NLOPS results to reproduce the NNLO Higgs rapidity spectrum, thereby achieving NNLOPS accuracy. In the present case instead, to properly take into account the decay of the boson to leptons, NNLOPS accuracy is reached by performing a three-dimensional rescaling of the events generated by Zj-MiNLO and Wj-MiNLO using NNLO distributions as computed, for instance, with DYNNLO. This implies that the calculation is numerically more intensive.
We have validated our procedure by considering observables typically used to study Drell-Yan processes. We have found extremely good agreement with NNLO results for observables fully inclusive over QCD radiation, and all the features expected in a computation matched with parton showers for more exclusive observables. We have also compared our NNLOPS predictions to state-of-the-art analytic resummation for observables sensitive to soft-collinear QCD radiation, which are not described accurately by a fixed-order NNLO calculation. Despite the logarithmic accuracy of our simulation cannot be claimed to reach the same precision as these analytic resummations, we have found reasonably good agreement for the vector boson transverse momentum. Slightly more pronounced discrepancies were instead observed for φ * and for the jet-veto efficiency for large values of the jet radius R.
We have also compared with a number of available experimental results, not only for fully inclusive observables but also for φ * , the transverse momentum of the vector boson, as well as k T splitting scales. The successful outcome of these comparisons is an indication that our computation can be used as a state-of-the art prediction for future studies where a simultaneous inclusion of NNLO corrections and parton shower effects are needed. We have also illustrated how our generator can be used to study the polarization of W and Z bosons produced in Drell-Yan events.
We finally remark that we did not include in this work electroweak corrections. These corrections are known exactly at one loop [8,11] and can be included, at O(α ew ) either additively or multiplicative on top of the QCD corrections included here. Version 3 of the program FEWZ [15] includes directly both EW and QCD corrections, however we could not use it to generate the three dimensional distributions needed here.
The code developed for this paper will soon be made available within the POWHEG BOX-V2 public repository. By parallelising on a cluster the runs needed to obtain the necessary inputs for the "NLO-to-NNLO" reweighter, results can be obtained in a reasonably short timescale. Nevertheless, as mentioned above, running this simulation from scratch requires more computing power than what was needed for the Higgs case. Therefore, together with the code that allows one to compute all the ingredients to achieve NNLOPS accuracy, we find it useful to release also the distributions dσ NNLO /dΦ B for W and Z production for three reference scales µ R = µ F = KM V with K = 0.5, 1, 2.0 (M V denotes the mass of the vector boson) and also provide the dσ MINLO A /dΦ B and dσ MINLO B /dΦ B distributions for 7 reference scales, as discussed in the text. In these cases, NNLOPS accuracy can then be reached by just running POWHEG with MiNLO switched on, and processing the generated events with the "NLO-to-NNLO" reweighter, using as input the aforementioned distributions.