NLO QCD+NLO EW corrections to diboson production matched to parton shower

We present the matching of NLO QCD and NLO EW corrections to parton showers for vector-boson pair production at the LHC. We consider leptonic final states, including resonant and non-resonant diagrams, spin correlations and off-shell effects. Our results are obtained interfacing the Recola2-Collier one-loop provider with the POWHEG-BOX-RES framework. We discuss our implementation, we validate it at fixed order, and we show our final results matched to parton shower. A by-product of our work is also a general interface between Recola2-Collier and POWHEG-BOX-RES. This is the first time that EW and QCD corrections to diboson production are consistently matched to parton showers.


Introduction
The pair production of massive vector bosons at the LHC (pp → VV , with VV = {W + W − , ZZ, W ± Z}) is among the most studied Standard Model (SM) processes, both as a signal on its own and as a background to physics beyond the Standard Model (BSM) and Higgs searches. Electroweak boson pair production involving at least a W boson in the final state (W + W − and ZW ± ) is important for collider phenomenology because it is sensitive to the ZWW gauge-boson self interaction, and therefore precision measurements of VV processes provide a test of the electroweak gauge structure. These precision tests are usually carried out by setting bounds on the allowed size of anomalous trilinear gauge couplings (aTGCs) [1], although several other ideas have been proposed to study effects due to BSM physics with VV final states [2][3][4][5][6][7][8][9][10][11][12][13]. Diboson production is also a background for several searches, notably those involving an heavy resonance decaying to a pair of gauge bosons. In particular, pp → W + W − and pp → ZZ are irreducible background for Higgs production, when the Higgs boson decays to gauge bosons. For all the above reasons, it is essential to make accurate predictions for vector boson pair production processes. Among the possible final states, the one where four leptons are present is probably the most interesting one, as it allows precise measurements, due to its clear signature. In the rest of this work, including this section, we will only discuss final states with four leptons, although we will often use the VV intermediate state as a shorthand from POWHEG and the matching to PS, as we will discuss in Sec. 2.
The paper is organized as follows: in Sec. 2 we describe the details of our computation, in Sec. 3 we list the parameters and cuts used throughout this work, and in Sec. 4 we discuss the validation of our implementation. In Sec. 5 we show a selection of the new results, i.e. the matching of NLO EW + NLO QCD corrections to parton shower. We summarize our work in Sec. 6. In the rest of this manuscript we will use the shorthands NLO QCD and NLO EW to denote NLO accuracy in the α S and α perturbative expansion, respectively. We use the notation NLO QCD + NLO EW to denote the additive combination of the hard matrix elements (in the POWHEGB function).

Details of the calculation
In this paper we consider the processes pp → e + ν e µ − ν µ , pp → µ + ν µ e − e + , We stress that the full matrix elements for four fermion production are used and no on-shell or double pole approximation is employed. In the following the three processes will be dubbed as WW , W Z, and ZZ production, and, collectively, as "diboson production". Although we will show results only for W + Z production, our code is fully general and W − Z production can be generated as well.
The calculation of the NLO QCD + NLO EW corrections to diboson production matched to QCD and QED parton shower presented in this paper is performed in the POWHEG-BOX-RES framework [65], which is a framework designed to simulate processes involving intermediate decaying resonances with NLO+PS accuracy. It is a new implementation of the POWHEG method [46,47] that overcomes the limitations of the POWHEG-BOX framework [49]. It has been used in Ref. [66] to simulate the process pp → bb ¯ νν with NLO QCD +PS accuracy, thereby achieving, for the first time, a fully-consistent treatment of tt and Wt production with two leptonic decays, in Ref. [67] to compute the processes pp → HV and pp → HV j production (V = W, Z) with NLO QCD + NLO EW +PS accuracy, and in Ref. [68] to compute the NLO EW +PS corrections to pp → νν j j. In Ref. [69], a simplified version of the POWHEG-BOX-RES algorithm has been implemented also in the W_EW-BMNNP and Z_EW-BMNNPV packages [69][70][71][72] of POWHEG-BOX-V2, in order to simulate neutral and charged Drell-Yan production with NLO QCD + NLO EW +PS accuracy in a fully-consistent manner. 2 A fully independent calculation of NLO QCD + NLO EW +PS corrections to Drell-Yan production was performed also in Ref. [74].
The improvement of the algorithm implemented in POWHEG-BOX-RES with respect to POWHEG-BOX-V2 is twofold. We summarize it briefly in this paragraph, and we refer to Ref. [65] for more details. On the one hand, the calculation of the NLO predictions needed for the event generation uses a modified version of the FKS [75] algorithm for the subtraction of the infrared (IR) singularities, that takes into account the resonant structure of the process under consideration through the concept of "resonance history". Not only does this modification improve the integration stability, but it also allows one to generate the POWHEG hardest radiation preserving the intermediate resonance(s) virtuality everywhere in the POWHEG Sudakov, preventing shape distortions in the matching to PS. On the other hand, POWHEG-BOX-RES can generate up to one "hard" radiation from each resonance (including the hard production process among the resonances): 3 the hardness of each radiation is to be used as a veto scale for PS evolution of the particles belonging to the resonance that emitted the considered radiation. As discussed in Ref. [69], the latter point is crucial when computing the NLO QCD + NLO EW +PS corrections to observables that are very sensitive to final-state QED radiation (FSR QED) but rather insensitive to initial-state QCD corrections (ISR QCD).
In order to implement the NLO QCD and the NLO EW corrections to diboson production in POWHEG-BOX-RES, we had to define the list of all the contributing LO and real partonic subprocesses together with the corresponding resonance histories, and to provide the required Born, virtual, and real matrix elements. We decided to code all the three classes of diboson-production processes (namely, four charged leptons, three charged leptons plus one neutrino, and two charged leptons plus two neutrinos) in the same POWHEG package and let the user select the desired one from the input card.
Concerning the resonance histories, we consider the t-channel ones with two decaying vector bosons (Fig. 1,  left), the s-channel ones involving triple gauge-boson interactions (for WW and W Z production, Fig. 1, center), and the peripheral ones involving the s-channel production of a vector boson that decays into a dilepton pair which radiates a second vector boson (Fig. 1, right). While the first two classes of resonance histories are by far the dominant ones in the typical event selections for diboson production, the third one can be important if a more inclusive analysis is considered (this could be the case, for instance, if the code is used to simulate background contributions to other processes with four final-state leptons). Moreover, we found that the inclusion of the peripheral histories improves the numerical stability of the calculation and strongly reduces the size of the "remnant" cross section. 4 All the needed Born, virtual, and real matrix elements are computed using the RECOLA2-COLLIER one-loop provider. RECOLA2 [77][78][79][80] is a library for the fully automated calculation of tree-level and one-loop matrix elements which relies on the COLLIER [81][82][83][84] library for the reduction of the tensor integrals and the evaluation of the scalar integrals coming from the one-loop diagrams. We use the SM-2.2.2 Recola model file to compute the NLO QCD + NLO EW corrections in the SM, but in principle our code can be easily generalized to use any BSM Recola model file as far as the considered extension of the SM does not involve any modification of the interactions between photons and fermions or among partons (as such modifications might have an impact on the IR subtraction performed by POWHEG-BOX and on the event generation). We developed a completely general interface between POWHEG-BOX-RES and RECOLA2 that can be used for other processes of interest. 5 2 A NLO QCD + NLO EW +PS implementation of the charged Drell-Yan case obtained with the POWHEG-BOX-V2 algorithm was obtained in Ref. [73]. 3 This corresponds to the so called allrad scheme, first introduced in Ref. [76]. 4 We stress the fact that we always employ full matrix elements: the definition of the peripheral resonance histories only affects the way POWHEG-BOX-RES performs the subtraction of the IR singularities and the integration. The concept of the "remnant" cross section in the POWHEG-BOX codes was introduced in Ref. [49]. 5 While the interface to RECOLA2 is general, the current treatment of the NLO EW corrections in POWHEG-BOX-RES is not, as it implies that each virtual process is in one-to-one correspondence with a LO process (so that it can be considered either as a NLO QCD correction or a NLO EW correction to the corresponding LO process), which in general is not the case for complicated processes. See for instance the O(α S α 6 ) and O(α 2 S α 5 ) corrections to pp → νν j j [85].
In order to deal with the presence of unstable particles, RECOLA2 implements the complex-mass scheme (CMS) [86][87][88]. In this scheme, the W and Z boson masses are promoted to complex numbers with the replacement M 2 and all the parameters derived from the gauge-boson masses (like, for instance, the sine of the weak mixing angle) get an imaginary part.
Concerning the calculation of EW corrections, RECOLA2 allows to perform the renormalization of the UV singularities in the SM using three possible input parameter schemes: The results shown in the following are computed in the (G µ , M W , M Z ) scheme, but our code gives the user the possibility to select any of the renormalization schemes mentioned above.
In the calculation of diboson-production cross sections and/or distributions, there are tree-level singularities coming from the presence of s-channel photon propagators that can go on-shell. In order to prevent these singularities, we impose both generation cuts and suitable phase-space suppression factors. For example, if we consider the process pp → e + e − µ + µ − , the generation cuts read: while the suppression factor is: where M e + e − and M µ + µ − are the invariant masses of the underlying-Born electronic and muonic pair, and the actual values of M cut and M supp should be chosen by the user and depend on the cuts applied during the analysis. On top of the suppression factor in eq. (3), we also provide a suppression factor of the form: where H T is the scalar sum of the transverse momenta of the charged leptons, in order to improve the generation efficiency in the typical diboson-event selection. Also for this suppression factor, the actual value of H supp T is chosen by the user. 6 Due to the presence of approximated radiation zeros in W Z production, we use the bornzerodamp option [49,54] of POWHEG-BOX-RES, that we keep activated for all the three processes at hand.
As a final remark, the contribution of the loop induced gg → ZZ and gg → W + W − processes is not included in our calculation. Even though these are O(α 2 S ) effects, their impact is not negligible because of the size of the gluon PDF. These processes can be computed, independently, at LO+PS using tools like GG2ZZ and GG2WW [23,25]. NLO QCD +PS results were presented in Ref. [58]. Photon-induced processes are not included in our calculation. As illustrated in Refs. [30-33, 62, 63], these contributions can be phenomenologically relevant. Dealing with initial-state photons requires extra features in the POWHEG-BOX-RES code, not available while we write. We plan to to include them in a future release of our code.

Input parameters and cuts
The input parameters used in the numerical simulations at √ s = 13 TeV are the following: All fermions are considered as massless, with the exception of the top quark. For this reason, we only provide results for dressed leptons. The on-shell values of the W and Z masses and widths are converted internally to the corresponding pole values with the relations: For W Z and ZZ production, we set the generation cut M cut of eq. (2) at 15 GeV, for each same-flavour oppositecharged lepton pair, and we apply the suppression factors in eqs. (3) and (4) with M supp = 30 GeV and H supp T = 4 GeV. We checked that our results do not depend on these technical parameters in the event selection under consideration.
The UV renormalization for the EW corrections is performed in the on-shell scheme with input parameters (G µ , M W , M Z ) supplemented with the CMS for the treatment of the unstable particles. The MS scheme is used for the renormalization of the NLO QCD corrections. In the following, the factorization and renormalization scales are are the vector bosons that define the signature under consideration), for constant scales, or to the invariant mass of the four-lepton system at the underlying-Born level, when using running scales. The results at NLO+PS accuracy are only computed for running scales.
The Cabibbo-Kobayashi-Maskawa (CKM) matrix is set to the identity matrix. However, in the code, the user can select a non-trivial quark mixing matrix: in this case, the NLO EW corrections are still computed with V CKM = 1 and then multiplied by the actual CKM values coming from the LO part of the amplitude as in Refs. [33,70].
In order to make contact with the results of Ref. [89], the NNPDF23_nlo_as_0118_qed PDF set [90][91][92] is used. However, the user can select any modern PDF set [93][94][95]. The PDF evolution as well as the evolution of the strong coupling constant is provided by the LHAPDF6 library [96].
As in Ref. [89], we always use the same value of α (namely, the one derived from G µ , i.e. α −1 132.357) both for the LO couplings and for the coupling when computing the NLO corrections, both real and virtual. This introduces a small mismatch when POWHEG is interfaced to the QED PS, since the PS uses α 0 for the photonfermion coupling (α −1 0 = 137.03599911). On the one hand, this mismatch is really small and hardly visible on the scale of our plots and, on the other hand, we allow the user to define two different values of α: one to be used for the LO couplings and a second one (corresponding to α 0 ) to be used in the additional power of α in the EW corrections. When this option is selected, POWHEG performs the subtraction of the IR singularities using α 0 , while the virtual and real matrix element are computed by RECOLA2 with a different value of α and then rescaled by a factor α 0 /α. For all diboson-production processes, the b quark is treated as massless, both in the initial and final state, when present. For WW production, we do not include the contribution of initial-state b quarks, in order to remove the real QCD channel gb → W + W − b which is enhanced by the presence of the top-quark resonance, but is usually subtracted in experimental analysis (single-top background).  Table 1: Integrated cross sections for the processes pp → e + ν e µ − ν µ , pp → µ + µ − e − e + , and pp → µ + ν µ e − e + at LO under the event selection in eqs. (7) and (8). The factorization scale is set to Table 2: Integrated cross sections for the processes pp → e + ν e µ − ν µ , pp → µ + µ − e − e + , and pp → µ + ν µ e − e + at LO under the event selection in eqs. (7) and (8). The factorization scale is set to the invariant mass of the four-fermion system.
We provide a dedicated interface for a consistent matching with the PYTHIA8.2 [97,98] PS, that will generate secondary QED and QCD emissions and finally convert partons into hadrons. As we will explain in Sec. 5, a dedicated interface is necessary because we use the allrad scheme in POWHEG.
In this paper we do not consider distributions involving jets, however, we provide a template analysis that can use FASTJET [99,100] to reconstruct them.
In order to make the discussion of the results easier, we use the same basic event selection for all dibosonproduction processes: where and are charged leptons, and ∆R is the separation in rapidity and azimuthal angle. For pp → e + e − µ + µ − and pp → e + e − µ + ν µ we also impose a leptonic mass window around the Z-boson mass: Both muons and electrons are dressed: photons are recombined with charged leptons if their angular distance ∆R( , γ) is less than 0.1.

Cross-checks and validation
In order to validate the implementation of the NLO QCD + NLO EW corrections to diboson-production processes in POWHEG-BOX-RES, we compare the predictions of our code with the ones of the Monte Carlo integrator used in  Table 3: Integrated cross sections for the processes pp → e + ν e µ − ν µ , pp → µ + µ − e − e + , and pp → µ + ν µ e − e + at NLO QCD under the event selection in eqs. (7) and (8) Table 4: Integrated cross sections for the processes pp → e + ν e µ − ν µ , pp → µ + µ − e − e + , and pp → µ + ν µ e − e + at NLO QCD under the event selection in eqs. (7) and (8) Table 5: Integrated cross sections for the processes pp → e + ν e µ − ν µ , pp → µ + µ − e − e + , and pp → µ + ν µ e − e + at NLO EW under the event selection in eqs. (7) and (8) Ref. [89] (MC in the following). Both codes use RECOLA2 for the calculation of the matrix elements, however, the integration and the subtraction of the IR singularities is performed in a completely independent way in the two programs. In particular, POWHEG uses the FKS subtraction (modified to take into account the presence of resonances), while in MC the Catani-Seymour procedure [101,102] is used.
Tables 1-5 collect the results at the integrated cross-section level under the event selection of eqs. (7) and (8) for the processes pp → e + ν e µ − ν µ , pp → µ + ν µ e − e + , and pp → µ + µ − e − e + (dubbed "WW", "WZ", and "ZZ" in the tables). Tables 1 and 2 show the results at LO for fixed and running scales, Tabs. 3 and 4 contain the results at NLO QCD for fixed and running scales, while the predictions at NLO EW accuracy (for fixed scales) are presented in Tab. 5. The NLO QCD corrections are positive, large and they are dominated by real QCD corrections: this is a consequence of the opening of gluon-induced channels (qg → VV q) at NLO QCD .
Since we are focused on the technical comparison of the two programs, we do not perform here scale-variation studies. The effect of scale variation can be read, for instance, from Ref. [89] and turns out to be large, given the size and the nature of the dominant contributions to the NLO QCD corrections.
The NLO EW corrections are negative and moderate at the integrated cross-section level. They are a combination of QED and purely weak effects.
In Tabs. 1-5, the numbers in parenthesis correspond to the statistical integration error on the last digit. As can be seen from Tabs. 1-5, the predictions of the two programs agree within the integration error.
The comparison at the differential distribution level is shown in Figs. 2 and 3 for the NLO QCD (running scales) and the NLO EW (fixed scales) predictions, respectively. For WW and ZZ production we consider the transverse momentum of the positron, p T (e + ), while for W Z production we take the transverse momentum of the e + e − pair, p T (e + e − ), i.e. the reconstructed Z. For ZZ production, we also present the results for the invariant mass of the µ + µ − pair, M(µ + µ − ). In Figs. 2 and 3, the upper panels show the differential distributions at LO (blue) and NLO (red) accuracy, the central panels show the relative NLO corrections (δ = NLO/LO − 1), while, in the lower panels, we plot the ratio of the NLO predictions from POWHEG and MC. As largely discussed in the literature, the NLO QCD corrections to the transverse momentum observables are positive, large, and increase with p T . On the contrary, the NLO QCD corrections to M(µ + µ − ) are flat and correspond to a normalization factor. The NLO EW corrections to the transverse-momentum distributions are negative and show the typical Sudakov behaviour [103][104][105][106][107][108][109] at high p T . The shape of the NLO EW corrections to M(µ + µ − ) is dominated by QED effects, since the radiation of a final-state photon reduces the invariant mass of the lepton pair, shifting the events from the peak of the LO distribution to the region below the Z resonance. As in the case of the cross-section level comparison, from the lower panels of Figs. 2 and 3 we conclude that POWHEG and MC agree within the statistical errors.
For the fixed-order part of the calculation, POWHEG computes the NLO QCD + NLO EW corrections additively. We checked that, when running our code with both QCD and EW corrections, δ QCD+EW from POWHEG is equal to the sum of δ QCD and δ EW computed with MC. We do not show here the plots, since the corresponding information can be read from the combination of Figs. 2 and 3.

Results at NLO+PS accuracy
In this section, we present the results at NLO accuracy matched to PS. For brevity, we only show results for ZZ and WW production, but the code can be used to generate events and perform a similar study for W ± Z, as well. We consider three different levels of accuracy: • NLO QCD + NLO EW + PS QCD,QED : full NLO corrections matched to the full PS with QED and QCD radiation (NLO α+α S + PS α,α S in the plots); • NLO QCD + PS QCD,QED : strong corrections matched to the full PS (NLO α S + PS α,α S in the plots); • NLO QCD + PS QCD : strong corrections matched to a PS without QED radiation (NLO α S + PS α S in the plots).
For the predictions at NLO QCD + NLO EW + PS QCD,QED , according to the allrad scheme, our code generates up to three emissions, namely ISR QCD or QED radiation, and FSR QED radiation from the decay products of each one of the vector bosons. The kinematics of the hard partonic event generated by POWHEG, together with the values of the transverse momenta (with respect to their emitters) of the generated partonic and/or photonic radiation, is then saved in the Les Houches (LH) event file. The transverse momentum of the initial-state radiation, if present, is used by the parton shower algorithm as upper bound for the generation of QED/QCD radiation from the hard production process. The transverse momentum of the photons from the final-state leptons (i.e. from the resonances) is used by the parton shower program as upper bound for further QED radiation. The results presented in this paper have been showered by PYTHIA8. This code allows to veto emissions harder than the ones generated by POWHEG by using dedicated UserHooks. We have also verified that we obtain fully compatible results if we let the PS generate unconstrained emissions and we subsequently check if the transverse momentum of the radiations with respect to the emitting particles is smaller than the POWHEG hardest ones. If this is not the case, we attempt to shower the event again until all constraints are met and the event is accepted.
For the predictions at NLO QCD + PS QCD,QED , the LH events contain at most one initial-state QCD radiation and the transverse momentum of the radiated parton sets the maximum hardness for the QCD PS, while the starting scale for the QED PS is the center of mass energy of the event for ISR, and the virtuality of the resonances for FSR. The predictions at NLO QCD matched to QCD PS are obtained from the same LH events used for the study at NLO QCD + PS QCD,QED accuracy, simply by turning off the QED radiation in PYTHIA.
In Figs. 4 and 5, the upper panels show the differential cross section as a function of the observable under consideration, the central panels contain the ratio of the results at NLO QCD + NLO EW + PS QCD,QED and the ones at NLO QCD + PS QCD,QED , while, in the lower panels, we show the ratio of the predictions at NLO QCD + NLO EW + PS QCD,QED and the ones at NLO QCD + PS QCD .
The calculation at NLO QCD + NLO EW accuracy matched to the full PS (PS QCD,QED ) includes the effect of the     In the lower panels, the ratios are taken with respect to a result where only QCD corrections are included. Therefore, on top of the same effects as in the central panels, these panels also include the effect of all-order photonic corrections (without approximations at O(α), and in PS approximation starting from O(α 2 )).
For the process pp → µ + µ − e − e + (ZZ production, Fig. 4), besides the observables used in Sec. 4 for the validation at NLO (p T (e + ) and M(µ + µ − )), we consider the transverse momentum of the hardest reconstructed Z boson, p T (Z 1 ), and the positron rapidity, y(e + ). For the transverse-momentum distributions (left plots), the predictions at NLO QCD + NLO EW + PS QCD,QED are always lower than the ones at NLO QCD + PS QCD,QED , or at NLO QCD + PS QCD , in particular at high p T , where the EW Sudakov corrections amount to approximately −30% with respect to the LO. The photonic corrections further reduce the predictions. The ratios for the positron rapidity distribution (bottom right plot) are essentially flat and, in the central panel, show an effect of about −3/ − 4% mainly coming from weak corrections that becomes approximately −10% in the lower panel, where the denominator does not include photonic corrections. The ratio of the predictions at NLO QCD + NLO EW + PS QCD,QED and the ones at NLO QCD + PS QCD,QED for the dimuon invariant-mass distribution (top right plot) is rather flat and shows that the effect of the EW corrections beyond QED PS amounts to about −4%. In the lower panel, there is a positive corrections below the Z peak coming from multiple photon radiation (radiative return) very similar to the one observed at fixed order in Fig. 3.
The predictions at NLO+PS accuracy for the process pp → e + ν e µ − ν µ (WW ) are collected in Fig. 5 as a function of the following observables: transverse momentum of the positron, p T (e + ), transverse momentum and invariant mass of the positron-muon system, p T (e + µ − ) and M(e + µ − ), and azimuthal distance between the positron and the muon, ∆φ (e + µ − ). For all the observables under consideration, the inclusion of the NLO EW corrections lowers the predictions with respect to the calculation at NLO QCD + PS QCD,QED (central panels). This effect is more pronounced in the tails of the transverse-momentum and invariant-mass distributions, where the NLO EW corrections are negative and large because of the EW Sudakov logarithms. A comment is in order concerning the p T (e + µ − ) distribution. For this observable, the QCD corrections are positive, large, and increase very steeply starting from about 100 GeV in the absence of a jet veto. 8 As a consequence, even small statistical fluctuations in the differential distributions in the upper panel of the corresponding plot end up in large statistical uncertainties on the NLO α+α S + PS α,α S / NLO α S + PS α,α S and NLO α+α S + PS α,α S / NLO α S + PS α S ratios. From the lower panels of Fig. 5, we conclude that the contribution of multiple photon radiation to the observables under consideration (with the event selection of eq. (7)) is negative with the only exception of the first few bins of the p T (e + ) and M(e + µ − ) distributions.

Conclusions
We computed the NLO QCD + NLO EW corrections to diboson production at hadron colliders matched to a complete parton shower, where QCD and QED radiation is simulated. For diboson production this is the first calculation where the NLO EW corrections have been consistently matched to QED PS. As the considered processes involve the production and the decay of unstable particles, whose decay products can radiate photons, the calcu-lation is based on the POWHEG-BOX-RES framework. The corresponding code is public and all the information for downloading it can be found in the POWHEG-BOX web page. 9 Though we did not perform an extensive phenomenological study that might be the subject of a future publication, we showed the potential of our code, and we pointed out that EW effects, consistently matched to QED parton shower, are relevant for several observables of interest.
The code relies on the RECOLA2 library for the calculation of the matrix elements for four leptons and four leptons plus photon/parton production. In particular, we developed a fully general interface between POWHEG-BOX and RECOLA2 that could be used for other processes. We performed our calculation in the Standard Model. However, given the possibility of RECOLA2 to compute tree-level and one-loop matrix elements in general extensions of the SM, in the future the code could be easily generalized to compute the NLO QCD + NLO EW corrections to diboson production matched to QCD and QED parton shower in the context of models beyond the SM, provided that the one-loop corrections are available in RECOLA2, and that the considered model does not alter the structure of QED and QCD interactions in a non-trivial way. If this were the case, the subtraction of infrared singularities and the Sudakov form factors in the POWHEG-BOX-RES should also be generalized.
As a final remark, the effect of NLO QCD corrections to diboson production is very large and it is dominated by real parton radiation, especially in the regions where one of the two vector bosons is soft with respect to the jet. It is thus important to include QCD corrections beyond NLO accuracy, for instance using a consistent merging of the predictions for VV and VV + jet production (V,V = W, Z) based on the MINLO [60] or MINNLO PS [111] procedures (at NLO and NNLO accuracy, respectively). The code presented in this paper can be taken as the starting point for the inclusion of the EW effects matched to QED PS in the treatment of diboson (+ jet) production in the MINLO or MINNLO PS framework.