QCD parton showers and NLO EW corrections to Drell-Yan

We report on the implementation of an interface between the SANC generator framework for Drell-Yan hard processes, which includes next-to-leading order electroweak (NLO EW) corrections, and the Herwig++ and Pythia8 QCD parton shower Monte Carlos. A special aspect of this implementation is that the initial-state shower evolution in both shower generators has been augmented to handle the case of an incoming photon-in-a-proton, diagrams for which appear at the NLO EW level. The difference between shower algorithms leads to residual differences in the relative corrections of 2-3% in the p_T(mu) distributions at p_T(mu)>~50 GeV (where the NLO EW correction itself is of order 10%).


Introduction
At high energy hadron colliders, studies of Drell-Yan (DY) like processes are of great importance. They are crucial for the understanding of QCD and electroweak (EW) interactions in hadron-hadron collisions. Drell-Yan processes have large cross sections and clean signatures in the detectors. They are used for monitoring of the collider luminosity and calibration of detectors. DY is a reference process for measurements of EW boson properties at hadron colliders. Combination of accurate experimental measurements of these processes with elaborated theoretical predictions allows the extraction of parton density functions (PDFs) in the kinematical regions which were not yet accessed in DIS experiments. DY processes provide an important background to many other processes studied at hadron colliders including searches for Higgs scalar as well as for W ′ and Z ′ bosons in particular. All this gives a strong motivation to have an advanced high precision theoretical description of DY. The experimental precision of DY measurements at the LHC can reach the 1% level. That means that the accuracy of the theoretical predictions needs to be even higher. For this reason it is obvious that QCD and electroweak radiative corrections should be taken into account.
This article presents the results of application of parton shower algorithms to a hard process that was calculated with electroweak radiative corrections in the SANC system [1,2]. The showering procedure was applied to the Drell-Yan processes: pp → (W + ) → l + ν l (γ) + X, pp → (γ, Z) → l + l − (γ) + X, (1.1) where X represents hadrons and l is one of e, µ, τ . The parton shower algorithms implemented in the general-purpose Monte Carlo generators Pythia8 [3] and Herwig++ [4] were used for these processes. It is worth noting that these two programs use essentially different parton shower algorithms: Pythia8 uses an evolution scheme based on transversemomentum ordering [5] and Herwig++ uses the coherent branching algorithm based on angular ordering of emissions in the parton shower [6].
Earlier the combination of the effects due to parton showers (PS) and due to EW radiative corrections for charged current process was considered in [7,8]. In those studies an interface between HORACE [9,10] and fortran event generator Herwig [11] was developed. Here, we present studies in which the SANC generator [12,13] is used for the treatment of complete NLO EW corrections with interfacing it to the Herwig++ v2.4.0 and Pythia8 v.130 generators to apply parton showers.
The paper is organized as follows. In the next sections we describe the chain of simulations and other topics associated with the calculation procedure. In section 2 the relevant features of the SANC MC event generators are described. Section 3 discusses aspects of the parton showers that are added by Pythia8 or Herwig++. Numerical cross checks and results are given in section 4. In section 5 the obtained results and prospects are discussed.

The Drell-Yan processes in SANC
The SANC system [12,13] provides tools for calculating the differential cross sections of the Drell-Yan processes taking into account the complete (real and virtual) O(α) electroweak radiative corrections. Here we give a brief summary of the main properties of this framework.
All calculations are performed within the OMS (on-mass-shell) renormalization scheme [14] in the R ξ gauge which allows an explicit control of the gauge invariance by examining a cancellation of the gauge parameters in the analytical expression of the squared matrix element.
We subdivide the total EW NLO cross section of Drell-Yan process at the partonic level for observables X ( X = (x 1 , ..., x n ) is a generic observable which is a function of the final-state momenta) into four terms: where σ Born is the Born level cross-section, σ virt is a contribution of virtual(loop) corrections, σ soft corresponds to a soft photon emission and σ hard is a contribution of a hard (real) photon emission. The terms with auxiliary parametersω (photon energy which separates phase spaces associated with the soft and hard photon emission) and λ (photon mass which regularizes infrared divergences) cancel out after summation and the differential EW NLO cross-section for infrared-safe observables does not depend on these parameters [15][16][17].
The tree level diagrams for the DY process are shown in figure 1 for neutral and charged currents. Examples of the diagrams corresponding to the electroweak NLO component for neutral and charged currents are shown in figures 2 and 3 respectively.
For real photon emission we separate contributions from initial and final state radiation and their interference in a gauge invariant way. In case of photon emission off the virtual W we introduce the splitting of the W-boson propagators by the following formula: The so-called on-shell singularities which appear in form of logarithms log(ŝ−M 2 W +iǫ) can be regularized by the W -width [18]: Electroweak NLO radiative corrections contain terms proportional to logarithms of the quark masses, log(ŝ/m 2 u,d ). They come from the initial state radiation contributions including hard, soft and virtual photon emission. Such initial state mass singularities are well known, for instance, in the process of e + e − annihilation. However, in the case of hadron collisions these logs have been already effectively taken into account in the parton density functions (PDF). In fact, in the procedure of PDF extraction from experimental data, the  QED radiative corrections to the quark line were not systematically subtracted. Therefore current PDFs effectively include not only the QCD evolution but also the QED one. Moreover, it is known that the leading log behaviours of the QED and QCD DGLAP evolution of the quark density functions are similar (proportional to each other). Consequently one gets an evolution of the PDF with an effective coupling constant where α s is the strong coupling constant, α is the fine structure constant, Q i is the quark charge, and C F is the QCD colour factor. We will use here the MS subtraction scheme, the DIS scheme may be used as well. A solution described in [19] allows to avoid the double counting of the initial quark mass singularities contained in our result for the corrections to the free quark cross section and the ones contained in the corresponding PDF. The latter should also be taken in the same The MS subtraction to the fixed (leading) order in α is given by: where q(x, M 2 ) is the parton density function in the MS scheme computed using the QED DGLAP evolution. The differential hadronic cross section for DY processes [1.1] is given by whereq 1 (x 1 , M 2 ),q 2 (x 2 , M 2 ) are the parton density functions of the incoming quarks modified by the subtraction of the quark mass singularities and σ q 1 q 2 →ℓℓ ′ is the partonic cross section of corresponding hard process. The sum is performed over all possible quark combinations for a given type of process (q 1 q 2 = ud, us, cd, cs for CC and q 1 q 2 = uū, dd, ss, cc, bb for NC). In our calculations we used fixed factorization scales M 2 = M 2 W for CC and M 2 = M 2 Z for NC. The effect of applying different EW schemes in the SANC system is discussed in [13]. In the current study we are using the G µ -scheme [20] since it minimizes EW radiative corrections to the inclusive DY cross section. In this scheme the weak coupling g is related to the Fermi constant and the W boson mass by equation where ∆r represents all radiative corrections to the muon decay amplitude [21]. Since the vertex term between charged particles and photons is proportional to g sin θ W , one can introduce an effective electromagnetic coupling constant which is evaluated from (2.7) in a tree-level approximation by setting ∆r = 0. The total NLO electroweak corrections to the total charged current DY cross section for 14 TeV pp collisions are estimated to be about −2.7% for the G µ -scheme and can reach up to 10% for the differential cross section in certain kinematical regions [12,13].
The EW NLO calculations for the DY processes were performed using semi-analytic calculations with the aid of the FORM symbolic manipulation system [22] and employ LoopTools [23] and SancLib [1] libraries for evaluation of scalar and tensor one-loop integrals. The analytical expressions for different components of the differential EW NLO cross-section for DY processes are realized within standard SANC Fortran modules which are used in our Monte Carlo event generators of unweighted events.

Photon induced contributions
At the O(α) level one can see that there is a non-zero probability to find a quasi-real photon inside one of the colliding protons. This brings up an additional QED contribution to the EW corrections, so called inverse bremsstrahlung. The complete set of O(α) photoninduced contributions for both NC and CC Drell-Yan processes was evaluated in [24]. The charged current results for this component were given by S. Dittmaier and M. Krämer in the proceedings to the Les Houches workshop [25]. The results for the neutral current were presented in [26], using an approach which implies effective resummation of the top quark one-and two-loop corrections in the LO cross section via s W renormalization: and with the function ρ (2) given in [27]. The coupling constant α Gµ is replaced by α Gµs 2 W /s 2 W in this approach. Table 1 presents comparison of SANC results with [26] with corresponding input parameters for the photon-induced contribution without these key differences in calculation schemes taken into account. In SANC the corresponding effects are considered as a part of the first-(and higher) order radiative corrections.
This comparison shows that although the size of the photon-induced contribution can differ by 10% or more between the two approaches, this corresponds to at most a per mille level difference in the total lepton pair cross section. It is therefore smaller than the aimed for accuracy and we do not need to consider this difference in further detail.
The fixed-order diagrams corresponding to the processes are shown in figures 4 and 5. The inverse bremsstrahlung component for the hadronic cross section can be written in a standard way: where f q (x 1 , M 2 ) and f γ (x 2 , M 2 ) are the parton density functions for quark and photon respectively. The quark mass singularity subtraction is performed for this contribution in analogous way to the processes with two quarks in the initial state. The photon induced channels are explicitly included in the SANC event generators. The corresponding correction value defined as δ γq = σ γq /σ 0 , where σ 0 is the tree level process cross section, is below the percent level for the total cross section, but reaches several percents in certain kinematic regions. The corrections for muon-neutrino pair transverse mass and µ + transverse momentum distributions in the charged current process pp → µ + ν for δ γq are shown in figure 6. The corrections for µ + µ − invariant mass and µ + transverse momentum distributions in the neutral current process pp → µ + µ − for δ γq are shown in figure 7. The large corrections for µ + transverse momentum in the charged current process in the region of p T (µ + ) > M W /2 is due to the recoil of a virtual W . The inverse bremsstrahlung contributions can be of resonant and non-resonant type. The latter have the incoming photon coupling to leptons and require a special colour flow interpretation in the code used to apply QCD parton showers to the hard process. As a workaround one can write an event entry for such contributions as a 2 → 3 process without internal structure. The resonant component can be treated in a standard way indicating a Z boson as a virtual propagator.

Parton Showers
In contrast to the fixed-order calculations described above, parton showers rely on an iterative (Markov-chain) branching procedure to reach arbitrary orders in the perturbative expansion. By keeping the total normalization unchanged, the shower explicitly conserves unitarity at each order, generating equal-magnitude but opposite-sign real and virtual corrections. Each branching step is based on universal splitting functions that capture the leading singularities of the full higher-order matrix elements exactly. Subleading terms can usually only be taken into account approximately, and hence different shower models (and "tunings") can give different answers outside the strict soft/collinear limits. Still, in practice, parton showers are reasonably accurate even for finite emission energies and angles, as long as the characteristic scale of each emission is hierarchically smaller than that of the preceding process (strong ordering). As such, they are complementary to the fixed-order truncations discussed above, which are accurate only in the absence of large hierarchies.
Several different shower formulations have been developed. In Herwig++ and Pythia8, which we shall be concerned with here, the shower approximation is cast in terms of evolution equations using DGLAP splitting kernels, which nominally capture only the leadinglogarithmic (LL) behaviour of higher perturbative orders. To further improve the accuracy, parton showers incorporate a number of improvements relative to the naive leading-log picture; 1) they use renormalization-group improved couplings by shifting the argument of α s for shower emissions to α s (p ⊥ ), thereby absorbing the β-function-dependent terms in the one-loop splitting functions into the effective tree-level ones, 2) they approximately incorporate the higher-order interference effect known as coherence by imposing angular ordering either at the level of the evolution variable (Herwig++) or in the construction of the shower phase space (Pythia8), 3) they enforce exact momentum conservation to all orders, albeit in different ways between the two (different "recoil strategies"), and 4) both programs include at least some further corrections due to polarization effects. The resulting approximation is thus significantly better than "pure LL", although it cannot be formally claimed to reach the NLL level.
Prior to the writing of this paper, the initial-state showers in both Herwig++ and Pythia8 included photons as emitted particles, but not as evolving ones. Interfacing of the photon-induced subprocesses to Herwig++ and Pythia8 therefore required a certain modification of these codes. In the two following subsections, we briefly summarize the main properties of these modifications, for each program respectively.

Processes with incoming photons in Pythia8
For the Pythia8 implementation of incoming photons, we re-use the existing backwardsevolution framework for gluons, with the modification that there is no photon self-coupling and replacing the q → g coupling and colour factor byᾱ = e 2 q α em /2π for q → γ. For future reference, we here summarize the steps specific to the photon backwards-evolution.
Denoting the Pythia8 evolution variable p 2 ⊥evol (see [28]), the evolution equation for a photon-in-a-proton is cast as a standard Sudakov evolution, with subsequent p ⊥evol "trial" scales generated according to an overestimate of the physical evolution probability, obtained by solving the trial evolution equation, where∆ is the trial shower Sudakov, representing a lower bound on the probability that there are no branchings between the two scales p 2 ⊥now and p 2 ⊥next , R is a random number distributed uniformly between zero and one, and the EM coupling used for the trial emission isᾱ =ᾱ(ŝ), withŝ being the CM energy of the two incoming partons (specifically, it is an overestimate of α(p 2 ⊥ ), which will be imposed by veto, below). The trial z integral,Î z , is defined byÎ where x γ is the momentum fraction carried by the incoming photon and the z limits are defined in [28].
Solving eq. (3.1) for p 2 ⊥next , we get Given a trial p ⊥next value obtained from this equation, we may now generate a corresponding trial z value according to where R ′ is a second random number distributed uniformly between zero and one, and the z limits are the same as those used in eq. (3.2). We have now obtained an importance-sampled pair of trial (p ⊥next , z) values. The quark flavour q is chosen with probability proportional to e 2 q f q (x γ /z, p 2 ⊥next ). Since the overestimates used for the importance sampling are not quite identical to the physical distributions we wish to obtain, the second step of the algorithm is to employ the veto algorithm and accept only those trials that lie inside the physical phase space with a probability, where the coupling factor translates the argument of α em in the manner mentioned above, the first PDF factor corrects the factorization scale used in the photon PDF to the new evolution scale, the second PDF factor corrects both the x and Q 2 of the would-be parent quark (or antiquark) to their correct post-branching values, and the last factor (the z factor) corrects for the form of P (z) used for the trial generation (the factor √ z, which may seem to complicate matters unnecessarily, arises since we use a factor 1/ √ z in the trial generation to suppress the high-x bump on valence quark distributions, which could otherwise lead to the trial generator not overestimating the physical splitting probability). If no acceptable branching is found above the global initial-state shower cutoff (cf. the documentation of Pythia8's spacelike showers [28]), the photon is considered as having been extracted directly from the beam remnant. Also note that the maximum expressed by eq. (3.2) could be violated, yielding P > 1 in eq. (3.5), if the photon PDF exhibits any thresholds or other sharp features. Further work would be needed to properly take into account such structures. We note, however, that the code forces the PDF to be bounded from below, so that a vanishing PDF should result in warnings, not crashes.
The (yet higher order) possibility of a fermion backwards-evolving to a photon has not yet been included in this framework. The net effect is therefore only to allow the initialstate shower off an incoming photon to reconstruct back to a quark or antiquark, but not the other way around.

Processes with incoming photons in Herwig++
The physics implementation of processes with incoming photons is very similar in Her-wig++ to that already described for Pythia8, so we do not go into as much detail. However, the practical implementation is somewhat different. Table 2: Event generation conditions referred to in the text as cut1.
When the Herwig++ parton shower is presented with a hard process with an incoming photon, it calls a PreShowerHandler of the IncomingPhotonEvolver class, which is specially written for this purpose. It generates a step of backward evolution from the photon to an incoming quark, in exactly the same way as described above. In particular, the transverse momentum of the q → qγ vertex is required to be smaller than the scale of the hard process, as in Pythia8. However, this backward evolution step is required to be generated with a probability of unity and if no backward step is generated above the infrared cutoff, or if it is generated outside the allowed phase space, then the evolution scale is reset to the hard scale and it loops back to try again. In a very small fraction of events it can happen that, due to mismatch between the hard process (SANC) and parton shower (e.g. in parton mass values and hadron remnant treatment) no backward step is possible. In such cases an EventException is thrown.
Having generated a backward step, the IncomingPhotonEvolver modifies the hard process to include it. That is, it replaces a γ → X event by the corresponding q → qX event, which the rest of Herwig++'s parton shower machinery operates on as normal. The quark line is correctly labeled as colour-disconnected from the rest of the hard process, so the colour coherence inherent to Herwig++'s shower ensures that it only radiates with opening angles smaller than the q → q scattering angle.
The IncomingPhotonEvolver has one parameter that may be of interest to users: minpT, the minimum transverse momentum generated for the q → qγ vertex. All of the plots shown below were generated with the default value of 2.0 GeV.

Cross Checks and Validation
Several cross check simulations were performed in order to verify the new implementation for the processes with incoming photons. The simulations included two steps: i ) hard event generation using the SANC MC generator for charged (CC) and neutral (NC) current cases, and ii ) addition of the parton showers using Herwig++ or Pythia8 generators. At the generation step the event selection shown in table 2 was applied. Events which satisfied these cuts were written in the Les Houches event format (LHEF [29,30]) for further processing.
The generators Pythia8 and Herwig++ in the second step were run with the following non-default configuration. In order to speed up the simulation without significantly influencing final results, multiple interactions and hadronization were turned off in both programs. The QED component of initial and final state radiation was disabled to avoid process cuts pp → W + → l + ν l (γ) + X M inv (µ + ν µ ) > 20 GeV, |η(µ + , ν µ )| < 2.5, p T (µ + , ν µ ) > 20 GeV Table 3: Selection criteria applied after showering procedure referred to in the text as cut2.
2332 (1)   double counting of the radiative corrections, which are calculated in SANC generator in the complete EW NLO approximation. In Herwig++ less strict than default kinematic constraints were set for photon momenta: k T (γ) > 0.0 and |η γ | < 10.
In order to avoid edge effects in the distributions after parton showers were applied the kinematic constraints were strengthened as shown in table 3. The lower limit on invariant mass of the leptons was not changed, which would not lead to edge effects since the transverse momenta constraints for W/Z decay products would indirectly increase the actual threshold for M ℓℓ by a factor of 2.
The physics setup corresponding to the LHC conditions used in this study is specified in [12,13]. The electroweak scheme for the calculations was chosen to be the G µ -scheme. As parton distribution functions the MRST2004QED [31] set was used since it allows to take the photon induced contribution into account. The factorization scale was set to M Z for neutral current case and to M W for the charged current.

Numerical Results
The results presented in this paragraph were obtained for statistics of 7 × 10 7 events for each channel (CC and NC) calculated in both LO and EW NLO approximation. The data produced in the wide selection criteria (cut1) were then subjected to the selection (cut2) with ∼ 50% efficiency. Table 4 shows the effect of this selection on the inclusive cross section and electroweak NLO correction values. Here δ denotes the relative corrections, δ = (σ N LO /σ LO − 1) × 100%. Expressions like "HP SANC + PS Herwig++ " denote the case when the hard process (HP) data were generated with the SANC generator and then processed with Herwig++ to apply parton showers (PS). The first and second rows in the table show the generator-level cross sections calculated with the SANC generator before parton shower algorithms applied in the cut1 and cut2 conditions, respectively. The third and fourth rows show effects of the cut2 selection when parton showers via Herwig++ and Pythia8 were applied.
To compare the parton shower algorithms in Pythia8 and Herwig++ it is convenient to introduce a parameter where dσ Pythia8 /dX represents a differential cross section by an observable X calculated with parton showers applied by Pythia8. Figures 8-13 show distributions for various observables obtained after the cut2 selection. Each figure contains three rows of plots with distribution of the differential cross sections themselves (top row), electroweak K-factor which is defined as usual as K = σ N LO /σ LO (middle row) and the R X value (bottom row). The distributions show that R X values can differ from unity by up to 10% for p T (right columns in figures 8,9,11,12) and by several percent for other observables.
The difference between the shower algorithms is most noticeable in the p T distributions. Nevertheless R X distributions in M inv (µ + µ − ) and M T (µ + ν µ ) are practically flat and differ from unity by only 2-3%.
It should here be emphasized that the prescription presented in this paper only concerns incorporating the first order of EW corrections into a shower framework. In particular, the description of QCD corrections is still handled only with leading-logarithmic precision, and does not include any matching to higher-order QCD matrix elements (see, e.g., [32,33]). Thus, the description of vector boson plus jets can only be expected to be correct for jets with p T ≪ m Z (representing the bulk of the cross section). For harder jets, differences between Herwig++ and Pythia8 reflect the uncertainty associated with QCD corrections beyond LL. Further work would be required to include QCD matrix-element corrections in this region.
The right columns of the plots in figures 8, 11 show the muon transverse momentum in the pp → µ + ν µ + X and pp → µ + µ − + X processes. Although the radiative corrections are washed out by parton showers in the peak region they reach up to 15% for higher p T values. The difference in the parton shower algorithms for EW RC is mildly noticeable at p T > 50 GeV where the K-factors for Pythia8 and Herwig++ diverge with maximum 2%. The small bends in the 20 GeV region are edge effects that appear in the showered events selection and play no role in the physics of the process.
A similar behaviour can be seen in the Z/W transverse momentum distributions in figures 9, 12: the K-factors deviate up to 4%. The (µ + µ − ) invariant mass and (µ + ν µ ) transverse mass plots in figures 10, 13 show no significant effects.

Summary
An interface between the SANC matrix-element generator and the Herwig++ and Pythia8 parton shower Monte Carlo codes has been presented. As part of this work, the new possibility of backwards evolution of photons has been added to both the Herwig++ and Pythia8 initial-state showers. Several numerical crosschecks have been performed, with reasonable results. The addition of parton showering gives a natural smearing effect on the EW K-factor distributions in the Drell-Yan process. The lepton p T distributions are mostly affected in the Z and W peak region.
The remaining difference between Pythia8 and Herwig++ showering algorithms was another focus of this study. The comparative plots included show that the difference in differential cross section can reach up to 10% for certain observables. We expect that this could be further reduced by extending the prescription presented here to include a matching to fixed-order QCD matrix elements for vector boson plus jets.
Since the completion of this work, two implementations of electroweak corrections to W boson production in the POWHEG framework have appeared [34,35], combining both EW and QCD corrections. However these works do not take into account effects of photon-induced processes. We consider that the implementation into a general-purpose electroweak tool like SANC has advantages for precision EW studies, given the importance of EW scheme-dependence in radiative corrections and the need for consistent scheme implementations between different processes. Nevertheless, it is clear that the POWHEG framework is an extremely powerful tool in describing and combining hard QCD and parton shower corrections consistently, and we look forward to making detailed comparisons between the results of these implementations and our own.