# 1 Introduction

## 1 Introduction

The production of vector bosons $`W^\pm `$, $`Z^0`$ and $`\gamma ^{}`$  in high energy hadronic collisions is one of the most important processes that should be investigated in order to test the Standard Model of electroweak interactions and Quantum Chromodynamics. By measuring the rapidity or the mass distribution of the leptonic decay products one can also investigate the quark and antiquark distributions inside the colliding hadrons. For such processes, generated at the lowest order via a hard scattering like $`q\overline{q}^{}V`$, where $`q`$ and $`q^{}`$ have the same flavour for $`Z/\gamma ^{}`$ production and different flavour for $`W`$ production, higher order corrections due to multiple gluon emission in the initial-state radiation will play a crucial role. Many analyses have been devoted to the phenomenology of vector boson production, particularly to the differential distribution with respect to the transverse momentum $`q_T`$ of the produced vector boson. The approach of resummation of large logarithms of the ratio $`m_V/q_T`$ has been followed in many cases. This was originally proposed by Dokshitzer, Dyakonov and Troyan (DDT)  and then accomplished by Collins, Soper and Sterman (CSS) , who performed the leading logarithmic resummation in the space of the impact parameter $`b`$, which is the Fourier space conjugate to $`q_T`$. Ladinsky and Yuan implemented the CSS results numerically . Resummations of the initial-state multiple emission have been performed in  in the $`b`$-space and more recently in  in both the $`b`$\- and the $`q_T`$-space. In  the resummation is also matched with the exact perturbative first-order result, which is important at high $`q_T`$. In the $`b`$-space approach non-perturbative effects in the region of large values of $`b`$ are taken into account via Gaussian functions in $`b`$, corresponding to a smearing of the transverse momentum distribution , which can also be directly implemented in $`q_T`$-space .

Another possible approach to studying the phenomenology of vector bosons is to use Monte Carlo simulations of the initial-state parton shower. Standard parton showers  are performed in the leading-log approximation, therefore they are reliable only in the soft or collinear region of the phase space, corresponding to low $`q_T`$ values for the produced vector boson. If we wish to study the high $`q_T`$ region of the spectrum it is necessary to provide parton showers with matrix-element corrections. Refs.  implement matrix-element corrections to simulations of vector boson production in the PYTHIA Monte Carlo event generator and compare them with the results obtained at the Tevatron collider by the DØ collaboration.

In this paper we reconsider this problem and apply matrix-element corrections to the initial-state radiation of the HERWIG parton shower, following the general prescription contained in , as we already did for $`e^+e^{}`$ annihilation , Deep Inelastic Scattering  and top quark decays .

It is worth recalling that at present no Monte Carlo program including the full next-to-leading order (NLO) results exists, as it is not known how to set up a full NLO calculation in a probabilistic way. When providing parton showers with matrix-element corrections we still only get the leading-order normalization, because in the initial-state cascade we only include leading logs and not the full one-loop virtual contributions.

In Section 2 we review the basis of the HERWIG parton shower algorithm for the initial-state radiation in hadronic collisions. In Sections 3 and 4 we discuss the hard and soft matrix-element corrections to vector boson production. In Section 5 we plot some relevant phenomenological distributions at the centre-of-mass energy of the Tevatron and of the LHC using the new version of HERWIG. We compare our results with previous versions of HERWIG, resummed calculations and experimental data. Finally, in Section 6 we discuss our results and make some concluding comments.

## 2 The parton shower algorithm

The production of a vector boson $`V`$ in hadronic collisions is given at lowest order by the elementary parton-level process $`q\overline{q}V`$. In the following, we shall assume that the vector boson decays into a lepton pair (Drell–Yan interactions). The first-order tree-level corrections to such a process are given by the processes $`q\overline{q}Vg`$ and $`qgVq`$ ($`\overline{q}gV\overline{q}`$), where the initial-state partons can come from either incoming hadron.

A possible method to implement the initial-state parton shower in a probabilistic way is the Altarelli–Parisi approach, in which the initial energy scale $`Q_0`$ is increased up to the probed value $`Q`$ and all the effect of the emitted partons is integrated out.

On the contrary, standard Monte Carlo programs  explicitly keep track of the accompanying radiation, by implementing the so-called ‘backward evolution’ in which the hard scale is reduced away from the hard vertex, tracing the hard scattering partons back into the original incoming hadrons and explicitly generating the distribution of emitted partons. In the leading infrared approximation, the probability of the emission of an additional parton from a parton $`i`$ is given by the general result for the radiation of a soft/collinear parton:

$$dP=\frac{dq_i^2}{q_i^2}\frac{\alpha _S\left(\frac{1z_i}{z_i}q_i\right)}{2\pi }P_{ab}(z_i)dz_i\frac{\mathrm{\Delta }_{S,a}(q_{i\mathrm{max}}^2,q_c^2)}{\mathrm{\Delta }_{S,a}(q_i^2,q_c^2)}\frac{x_i/z_i}{x_i}\frac{f_b(x_i/z_i,q_i^2)}{f_a(x_i,q_i^2)}.$$
(1)   
The HERWIG parton shower is ordered according to the variable $`q_i^2=E^2\xi _i`$, where $`E`$ is the energy of the parton that split and $`\xi _i=\frac{p_hp_i}{E_hE_i}`$, where $`p_i`$ is the four-momentum of the emitted parton; $`p_h`$ is a lightlike vector with momentum component parallel to the incoming hadron; $`E_h`$ and $`E_i`$ are the energy components of $`p_h`$ and $`p_i`$; and $`z_i`$ is the energy fraction of the outgoing space-like parton (which goes on to participate in the hard process) with respect to the incoming one (i.e. $`z_i=1E_i/E`$). In the approximation of massless partons, we have $`\xi _i=1\mathrm{cos}\theta `$, where $`\theta `$ is the emission angle from the incoming hadron direction. When all emission is soft, the energy of the emitted partons is negligible ($`E_iE`$), therefore ordering according to $`q_i^2`$ corresponds to angular ordering; when the emission is hard, the energy of the radiated parton is similar to that of the splitting parton, so $`q_i^2`$ ordering is equivalent to transverse momentum ordering. In (1) $`f_a(x_i,q_i^2)`$ is the parton distribution function for the partons of type $`a`$ in the initial-state hadron, $`x_i`$ being the parton energy fraction. At each step of the backward evolution a parton of type $`a`$, a quark for example, can evolve back to any other type of parton $`b`$, in this case either a quark of the same flavour or a gluon, having a higher value of $`x_i`$.

The quantity $`\mathrm{\Delta }_{S,a}(q_i^2,q_c^2)`$ is the Sudakov form factor, resulting from the leading-logarithmic resummation and representing the probability that no resolvable radiation is emitted from a parton of type $`a`$ whose upper limit on emission is $`q_i^2`$, with $`q_c^2`$ being, in the case of HERWIG, a cutoff on transverse momentum. This cutoff implies a minimum value of the evolution scale $`q_i^2`$ that can be reached, $`q_i^2>4q_c^2`$, but in practice this is smaller than the smallest scale at which most standard parton distribution function sets are reliable, so an additional cutoff on $`q_i^2`$ has to be applied. The ratio of form factors appearing in Eq. (1) represents the probability that the emission considered is the first, i.e. the one with the highest value of $`q_i^2`$. In terms of Feynman diagrams, the Sudakov form factor sums up all-order virtual and unresolved contributions. $`P_{ab}(z_i)`$ is the Altarelli-Parisi splitting function for a parton of type $`b`$ to evolve to one of type $`a`$ with momentum fraction $`z_i`$. $`\alpha _S`$ is the strong coupling, evaluated at a scale of order the transverse momentum of the emitted parton, which sums large higher order corrections. This, together with the angular ordering condition, makes HERWIG accurate to next-to-leading order at large $`x`$ .

The definition of the variables $`q_i^2`$ and $`z_i`$ is not Lorentz-invariant, but it is frame-dependent. Colour coherence implies that for any pair of colour-connected partons $`i`$ and $`j`$ the maximum values of the $`q`$ variables are related by $`q_{i\mathrm{max}}q_{j\mathrm{max}}=p_ip_j`$. Therefore one is free to choose the frame in which to define the initiating values $`q_{i\mathrm{max}}`$ and $`q_{j\mathrm{max}}`$, with the only prescription being that their product must equal $`p_ip_j`$. The subsequent emissions are then ordered in $`q_i^2`$. For vector boson production, as in most cases, symmetric limits are fixed by HERWIG, i.e. $`q_{i\mathrm{max}}^2=q_{j\mathrm{max}}^2=p_ip_j`$ and the energy of the parton which initiates the cascade is set to $`E=q_{\mathrm{max}}=\sqrt{p_ip_j}`$. Ordering according to $`q_i^2`$ therefore dictates $`\xi _i<z_i^2`$.

After we generate the initial-state shower, the original partons are not on their mass-shell anymore, so their energy and momentum cannot be conserved. Energy-momentum conservation is then achieved by applying a separate boost to each jet along its own direction. As a result of this, the jet momenta are no longer equal to the parton ones, but energy and momentum are globally conserved and the vector boson acquires a transverse momentum from the recoil against the emitted partons. Since the mass shift becomes negligible in the soft and collinear limits, the precise details of this kinematic reshuffling are not fixed a priori, but are free choices of the model.

Once the backward evolution has terminated, a model to reconstruct the original hadron is required. In HERWIG, if the backward evolution has not resulted in a valence quark, additional non-perturbative parton emission is generated to evolve back to a valence quark. Such a valence quark has a Gaussian distribution with respect to the non-perturbative intrinsic transverse momentum in the hadron, with a width that is an adjustable parameter of the model. In the following, when discussing the phenomenological implications of our work, we shall consider both HERWIG’s default value of zero, and an increased value of 1 GeV, bracketing the reasonable range of non-perturbative effects.

The algorithm so far discussed is reliable only in the soft or collinear limits and, since it only describes radiation for $`\xi _i<z_i^2`$, there are regions of the phase space that are completely empty (‘dead zones’). The radiation in such regions, according to the full matrix element, should be suppressed, but not completely absent as happens in HERWIG. We therefore need to improve the HERWIG model by implementing matrix-element corrections. As usual, this method works in two steps: we populate the missing phase space region by generating radiation according to a distribution obtained from the first-order matrix-element calculation (‘hard corrections’); we correct the algorithm in the already-populated region using the matrix-element probability distribution whenever an emission is capable of being the ‘hardest so far’ (‘soft corrections’).

## 3 Hard Corrections

In order to implement the hard and soft matrix-element corrections to simulations of the initial-state radiation in Drell–Yan processes, we firstly have to relate the HERWIG variables $`\xi `$ and $`z`$ to the kinematic ones we use in the matrix-element calculation. For the process $`q(p_1)\overline{q}(p_2)g(p_3)V(q)`$ we parametrize the phase space according to the Mandelstam variables $`\widehat{s}=(p_1+p_2)^2`$, $`\widehat{t}=(p_1p_3)^2`$ and $`\widehat{u}=(p_2p_3)^2`$. Throughout this paper, we neglect the parton masses, so we have $`\widehat{s}+\widehat{t}+\widehat{u}=m_V^2`$.

The phase space limits, in terms of the variables $`\widehat{s}`$ and $`\widehat{t}`$, are:

$`m_V^2<`$ $`\widehat{s}`$ $`<s,`$ (2)   
$`m_V^2\widehat{s}<`$ $`\widehat{t}`$ $`<\mathrm{\hspace{0.33em}\hspace{0.33em}0},`$ (3)   
where $`s`$ is the squared energy in the centre-of-mass frame. Note that the point $`\widehat{s}=m_V^2`$ corresponds to the soft singularity, and the lines $`\widehat{t}=0`$ and $`\widehat{t}=m_V^2\widehat{s}`$ to collinear emission.

In order to relate $`\widehat{s}`$ and $`\widehat{t}`$ to $`\xi `$ and $`z`$ we use the property that the mass $`m`$ and the transverse momentum $`p_t`$ of the $`q`$$`g`$ ($`\overline{q}`$$`g`$) jets are conserved in the showering frame. In doing this, we observe that, in the approximation of massless partons, the energy of the annihilating $`q\overline{q}`$ pair which produces the vector boson $`V`$ is equal to $`E^{}=\sqrt{p_qp_{\overline{q}}}=\sqrt{m_V^2/2}`$.

In terms of the showering variables, we obtain:

$`m^2`$ $`=`$ $`{\displaystyle \frac{1z}{z^2}}\xi m_V^2,`$ (4)   
$`p_t^2`$ $`=`$ $`{\displaystyle \frac{(1z)^2}{2z^2}}\xi (2\xi )m_V^2,`$ (5)   
and in terms of the matrix-element variables:

$`m^2`$ $`=`$ $`\widehat{t},`$ (6)   
$`p_t^2`$ $`=`$ $`{\displaystyle \frac{\widehat{u}\widehat{t}}{\widehat{s}}}.`$ (7)   
Combining them we get the following equations:

$`z`$ $`=`$ $`{\displaystyle \frac{m_V^2}{\widehat{t}}}+\sqrt{\left({\displaystyle \frac{m_V^2}{\widehat{t}}}\right)^2{\displaystyle \frac{2m_V^2}{\widehat{s}\widehat{t}}}(m_V^2\widehat{t})},`$ (8)   
$`\xi `$ $`=`$ $`2{\displaystyle \frac{\frac{m_V^2}{\widehat{t}}\frac{m_V^2\widehat{t}}{\widehat{s}}+m_V^2\sqrt{\frac{1}{\widehat{t}^2}\frac{2}{\widehat{s}\widehat{t}}+\frac{2}{m_V^2\widehat{s}}}}{1\frac{m_V^2}{\widehat{t}}m_V^2\sqrt{\frac{1}{\widehat{t}^2}\frac{2}{\widehat{s}\widehat{t}}+\frac{2}{m_V^2\widehat{s}}}}}.`$ (9)   
The region where HERWIG does not allow gluon radiation can be derived by solving the equation $`\xi >z^2`$:

$`\widehat{s}_{\mathrm{min}}<`$ $`\widehat{s}`$ $`<s`$ (10)   
$`\widehat{t}_{\mathrm{min}}<`$ $`\widehat{t}`$ $`<\widehat{t}_{\mathrm{max}},`$ (11)   
where $`\widehat{t}_{\mathrm{max}}`$ can be obtained by solving the equation

$$\widehat{t}^2+3m_V^2\widehat{t}+2m_V^4\left(1\frac{m_V^2}{\widehat{s}}\right)=0,$$
(12)   
or:

$$\widehat{t}_{\mathrm{max}}=\frac{m_V^2}{2}\left(3\sqrt{1+8m_V^2/\widehat{s}}\right).$$
(13)   
It is straightforward to write $`\widehat{t}_{\mathrm{min}}`$ as

$$\widehat{t}_{\mathrm{min}}=m_V^2\widehat{s}\widehat{t}_{\mathrm{max}},$$
(14)   
while $`\widehat{s}_{\mathrm{min}}`$ can be determined by the condition $`\widehat{t}_{\mathrm{min}}(\widehat{s})<\widehat{t}_{\mathrm{max}}(\widehat{s})`$:

$$\widehat{s}_{\mathrm{min}}=\frac{m_V^2}{2}\left(7\sqrt{17}\right).$$
(15)   
In Fig. 1 we plot the total phase space and HERWIG’s limits for a vector boson mass of $`m_V=80`$ GeV and a centre-of-mass energy of $`\sqrt{s}=200`$ GeV, in terms of the normalized variables $`\widehat{s}/m_V^2`$ and $`\widehat{t}/m_V^2`$.

We can see that, as in cases  and  and differently from , the soft and the collinear singularities are well inside the HERWIG phase space: we also have an overlapping region, corresponding to a kinematic configuration in $`\widehat{s}`$ and $`\widehat{t}`$ where radiation can come from either parton. Note that the only dependence on the external physical parameters is the position of the edge at large $`\widehat{s}`$, i.e. $`\widehat{s}_{\mathrm{max}}=s`$, while the values of $`\widehat{s}_{\mathrm{min}}/m_V^2`$ and the limits in $`\widehat{t}/m_V^2`$ are independent of the centre-of-mass energy and of other kinematic conditions like the vector boson rapidity.

Once we have the total and HERWIG phase space limits, in order to implement matrix-element corrections, we have to apply the exact differential cross section in the dead zone. In  a general prescription is given to allow first-order corrections to quark scattering and annihilation processes once a generator of the lowest order process is available. For the Drell–Yan case, assuming that the virtuality and the rapidity of the produced vector boson are fixed by the Born process $`q\overline{q}V`$, the first-order differential cross section $`d\sigma `$ is proportional to the parton-level lowest order $`\sigma _0`$ according to the relation:

$$d^2\sigma =\sigma _0\frac{f_{q/1}(\chi _1)f_{\overline{q}/2}(\chi _2)}{f_{q/1}(\eta _1)f_{\overline{q}/2}(\eta _2)}\frac{C_F\alpha _S}{2\pi }\frac{d\widehat{s}d\widehat{t}}{\widehat{s}^2\widehat{t}\widehat{u}}\left[(m_V^2\widehat{u})^2+(m_V^2\widehat{t})^2\right].$$
(16)   
In the above equation $`f_{q/1}(\chi _1)`$ and $`f_{\overline{q}/2}(\chi _2)`$ are the parton distribution functions of the scattering partons inside the incoming hadrons 1 and 2 for energy fractions $`\chi _1`$ and $`\chi _2`$ in the process $`q\overline{q}Vg`$, while $`f_{q/1}(\eta _1)`$ and $`f_{\overline{q}/2}(\eta _2)`$ refer to the Born process and cancel off the factors that are already in $`\sigma _0`$. The assumption that the rapidity and mass of the vector boson are the same as in the process $`q\overline{q}V`$ allows us to recover the Born result in the limit of an extremely soft gluon radiation.

As stated here, Eq. (16) is a trivial rewriting of the first-order differential cross section, but the main point of  is that if the azimuth of the emitted gluon is generated in the right way, Eq. (16) correctly describes the full process including the vector boson decay. Thus, to implement our matrix-element corrections, we do not need to know anything about the final state of the vector boson – its properties are correctly inherited from the Born process.

In a similar way we deal with the Compton process $`q(p_1)g(p_3)q(p_2)V(q)`$. We define the Mandelstam variables $`\widehat{s}=(p_1+p_3)^2`$, $`\widehat{t}=(p_3p_2)^2`$ and $`\widehat{u}=(p_1p_2)^2`$, and we find the same expressions for the variables $`\xi `$ and $`z`$ and for the phase space limits. We obtain for the differential cross section :

$$d^2\sigma =\sigma _0\frac{f_{q/1}(\chi _1)f_{g/2}(\chi _2)}{f_{q/1}(\eta _1)f_{\overline{q}/2}(\eta _2)}\frac{T_R\alpha _S}{2\pi }\frac{d\widehat{s}d\widehat{t}}{\widehat{s}^3\widehat{t}}\left[(m_V^2\widehat{t})^2+(m_V^2\widehat{s})^2\right].$$
(17)   
Extending this formula to processes where we have an incoming antiquark or where the gluon belongs to hadron 1 is straightforward. We then generate events according to the above distributions in the dead zone using standard techniques.

When applying the hard corrections, in principle one should also implement the form factor, but, since we are quite far from the soft and collinear singularities, it is actually not important and we shall neglect it in the following. This is justified by the fact that the total fraction of events that receive an emission from the hard correction is small. For example for $`W`$ production it is 3.9% at the Tevatron and 9.2% at the LHC. Also, the fact that such fractions are quite small allows us to neglect multiple emissions in the dead zone and makes the use of the exact first-order result reliable.

Among these events, the fraction of $`q\overline{q}^{}Wg`$ processes is 53.5% at the Tevatron and 24.5% at the LHC. The reason for the differences between the two machines can be understood in terms of the parton distribution functions. The gluon density inside the protons is higher when the colliding energy is increased because $`x`$ is decreased; moreover at the LHC we have $`pp`$ interactions instead of $`p\overline{p}`$, therefore a $`q\overline{q}^{}`$ annihilation requires an antiquark $`\overline{q}^{}`$ to be taken from the ‘sea’.

The equivalent numbers for $`Z`$ production are essentially identical.

## 4 Soft Corrections

According to , we should also correct the emission in the region that is already populated by HERWIG using the exact first-order calculation for every emission that is the hardest so far. This can be performed by multiplying the parton shower distribution by a factor that is equal to the ratio of HERWIG’s differential distribution to the matrix-element one. The only non-trivial part of this is in calculating the Jacobian factor $`J(\widehat{s},\widehat{t};z,\xi )`$ of the transformation $`(z,\xi )(\widehat{s},\widehat{t})`$. HERWIG’s cross section is then given by

$$\frac{d^2\sigma }{d\widehat{s}d\widehat{t}}=\frac{d^2\sigma }{dzd\xi }J(\widehat{s},\widehat{t};z,\xi ),$$
(18)   
where $`d^2\sigma /dzd\xi `$ is given by the elementary emission probability given in Eq. (1). The Jacobian factor $`J`$ can be simply calculated from the relations given earlier:

$$J(\widehat{s},\widehat{t};z,\xi )=\frac{\widehat{t}(m_V^2\widehat{t})}{\widehat{s}^2}\frac{z^5}{m_V^4\xi (1z)^2\left(z+\xi (1z)\right)}.$$
(19)   
At this point we are able to make some comparisons with the approach that is followed in  where matrix-element corrections are added to the PYTHIA simulation of vector boson production. The parton shower probability distribution is applied over the whole phase space (in its older versions PYTHIA had a cutoff on the virtuality $`k^2`$ of the hard scattering parton that was constrained to be $`k^2<m_W^2`$ in order to avoid double counting) and the exact matrix-element correction is applied only to the branching that is closest to the hard vertex. Unlike , we have complementary phase space regions where we apply either the parton shower distribution (1) or the exact matrix-element ones (16,17), while in the parton shower region ($`\xi <z^2`$) we use the exact amplitude to generate the hardest emission so far instead of just the first emission. Correcting only the first emission can lead to problems due to the implementation of the Sudakov form factor whenever a subsequent harder emission occurs, as we would get the unphysical result that the probability of the hard emission would depend on the infrared cutoff that appears in the expression of the form factor. See  for more details on this point.

## 5 Results

Having implemented the hard and soft matrix-element corrections to the initial-state radiation in Drell–Yan interactions, we wish to investigate the impact they have on relevant phenomenological observables that can be measured at the Tevatron and in future at the LHC. In the following, we shall mostly concentrate on $`W`$ production although $`\gamma ^{}`$ and $`Z`$ events are treated in exactly the same way.

### 5.1 Vector boson transverse momentum

A particularly significant phenomenological quantity is the transverse momentum of the $`W`$ to the beam axis, which has been object of many theoretical and experimental analyses. In the soft/collinear limit, the transverse momentum of the $`W`$ is constrained to be $`q_T<m_W`$, since in the hard process $`q\overline{q}^{}W`$ the $`W`$ is produced with no transverse momentum and it can acquire some $`q_T`$ only as a result of the initial-state parton showering. When the emission is generated according to the exact matrix element of processes like $`q\overline{q}^{}Wg`$, $`qgq^{}W`$ or $`\overline{q}g\overline{q}^{}W`$, the $`W`$ produced in the hard process is allowed to have a non-zero $`q_T`$ and events with $`q_T>m_W`$ are expected.

In Fig. 2 we compare the differential cross sections with respect to the $`W`$ $`q_T`$ for $`p\overline{p}`$ collisions at the Tevatron energy<sup>*</sup><sup>*</sup>*The next Tevatron run will be at the slightly higher energy of 2 TeV. For the sake of comparison with existing data we use 1.8 TeV, but the results would not be qualitatively different for 2 TeV., $`\sqrt{s}=1.8`$ TeV, obtained using HERWIG 5.9, the latest public version, with 6.1, the new version in progress where we include for the first time matrix-element corrections to the initial-state parton shower in vector boson production.

We set an intrinsic transverse momentum equal to zero and use the MRS (R1) parton distribution functions . We can see from the plots that the impact of matrix-element corrections is negligible at low values of $`q_T`$ but it is quite relevant at high $`q_T`$, where we have many more events with respect to the 5.9 version. Above some value of $`q_T`$ HERWIG 5.9 does not generate events anymore, while the 6.1 version still gives a non-zero differential cross section thanks to the events generated via the exact hard matrix element. As in $`e^+e^{}`$ annihilation and DIS, it is actually the hard corrections that have a marked impact on our distributions, while the effect of the soft ones is quite negligible.

It is also interesting to plot the $`q_T`$ spectrum obtained running the ‘$`W`$ \+ jets’ process of HERWIG forcing the produced $`W`$ to decay leptonically. This generates the hard process $`q\overline{q}^{}Wg`$ (or the equivalent ones with an initial-state gluon) for all events. As this matrix element diverges when the transverse momentum of the $`W`$ approaches zero, HERWIG applies a user-defined cut on the $`q_T`$ generated in the hard process, which we set to 10 GeV. In our plot, this does not appear as a sharp cutoff since the $`W`$ gets some recoil momentum due to the initial-state parton shower which can increase or decrease its transverse momentum. If, on the contrary, we had plotted the $`q_T`$ of the $`W`$ generated in the hard $`22`$ process, we would have got a sharp peak at $`q_T=`$ 10 GeV. The agreement we find between the simulations of Drell–Yan processes provided with matrix-element corrections and the ‘$`W`$ \+ jet’ events for large $`q_T`$ reassures us that the implementation of the hard corrections is reliable.

We also see that the HERWIG distributions for Drell–Yan processes show a sharp peak in the first bin, which includes the value $`q_T=0`$: it corresponds to a fraction of events with no initial-state radiation and so $`W`$ bosons produced with zero transverse momentum. With any fixed infrared cutoff value, one expects a non-zero (though exponentially suppressed) fraction of events to give no resolvable radiation. These would normally be smeared out by non-perturbative effects like the intrinsic transverse momentum, as we shall see in later plots, but since we set the width of its distribution to zero by default, all such events appear in the lowest $`q_T`$ bin. This is actually a technical deficiency of the Monte Carlo simulation and not a detectable physical effect.

The DØ collaboration recently published data on the transverse momentum distribution of $`W`$ bosons at the Tevatron. In Fig. 3 we compare the HERWIG predictions with it.

In order to contribute to the investigation of possible effects of non-perturbative physics, we also run HERWIG setting an intrinsic partonic transverse momentum equal to 1 GeV, which we consider to be the maximum reasonable value.

We see that the data has a significantly broader distribution at small $`q_T`$ than HERWIG without intrinsic transverse momentum and that increasing the r.m.s. $`p_t`$ to 1 GeV is nowhere near enough to account for this. Furthermore the description of the data in the intermediate $`q_T`$ range, 30–70 GeV, is also rather poor. The intrinsic $`p_t`$ does not affect the predictions significantly for $`q_T`$ values above about 10 GeV, so there is no obvious way to improve the fit for intermediate values.

However, these effects are actually because the DØ data are uncorrected for detector effects. We have run HERWIG through DØ’s fast simulation program, CMS, and show results in Fig. 4. We see that HERWIG now describes the data rather well. The detector smearing is so strong at low $`q_T`$ that the additional smearing produced by the intrinsic transverse momentum becomes irrelevant.

At this point it is worthwhile commenting on the results shown in Refs. . Both compare generator-level results with the DØ data. Ref. ’s actually look rather similar to our generator-level results, so it is likely that after applying detector corrections they will describe the data as well as HERWIG. Ref.  found good agreement with the DØ data, but only after increasing the intrinsic transverse momentum to 4 GeV. It seems likely that this accounts for the smearing at low $`q_T`$ and would not be necessary after including detector smearing. However, in the intermediate $`q_T`$ range the results of  are significantly lower than HERWIG. Since we find a detector correction of around a factor of two in this region, it would be very interesting to see the results of  at detector level to see whether they are still able to fit the data.

In Fig. 5 we show the $`q_T`$ distributions for $`pp`$ collisions at the energy of the LHC, $`\sqrt{s}=14`$ TeV, and find that the impact of the corrections is even bigger once the energy is increased.

Unlike in the Tevatron transverse momentum distributions, we do not have the previously-mentioned sharp peak at $`q_T=0`$: this is because at the LHC we have $`pp`$ interactions and the protons do not have valence antiquarks, while in order to produce a $`W`$ we do need a $`q\overline{q}^{}`$ hard scattering. As a result, the backward evolution has to produce at least one splitting, which always gives the $`W`$ itself some transverse momentum. From the window in the top-right corner, we also see that at very low $`q_T`$ the uncorrected version, 5.9, has a few percent more events than 6.1, particularly in the case of the LHC. As we said in the introduction, although we have matched to the tree-level NLO matrix elements, we still get the LO normalization, therefore the total cross sections obtained from versions 5.9 and 6.1 are the same. Since at the energy of the LHC we are generating a higher fraction of events at large $`q_T`$ via the exact matrix-element distribution, it is reasonable that this enhancement is partially compensated by a slight suppression in the low $`q_T`$ region.

In the region of low $`q_T`$, it is worthwhile comparing the HERWIG distributions with some resummed calculations that are available in the literature. All these calculations are based on the approach suggested in  where the differential cross section with respect to the vector boson $`q_T`$ is expressed as the resummation of logarithms $`l=\mathrm{log}(m_V^2/q_T^2)`$ to all orders in $`\alpha _S`$. Two conflicting nomenclatures are used in the literature to denote which logarithms are summed: in the differential cross section, at each order in $`\alpha _S`$ the largest term is $`1/q_T^2\alpha _S^nl^{2n1}`$, which are sometimes known as the leading logarithms, $`1/q_T^2\alpha _S^nl^{2n2}`$ being known as the next-to-leading logarithms, and so on. According to this classification, the results in  and  are NNLL and NNNLL respectively.

However, these logarithms ‘exponentiate’, allowing the differential cross section to be written in terms of the exponential (the ‘form factor’) of a series in $`\alpha _S`$ whose largest term is $`\alpha _S^nl^{n+1}`$, which are also sometimes known as the leading logarithms, $`\alpha _S^nl^n`$ being known as the next-to-leading logarithms, and so on. In  all NLL terms according to this nomenclature are summed in the form factor, which is evaluated either in the impact parameter $`b`$-space or in the $`q_T`$-space. The non-perturbative contribution is taken into account in the $`b`$-space formalism following the general ideas in  and setting Gaussian functions in the impact parameter $`b`$ to quantify these effects.

In Fig. 6 we compare the HERWIG 6.1 differential cross sectionAs the resummed calculations deal with a fixed value of $`m_W`$, Fig. 6 is obtained running HERWIG with a vanishingly small $`W`$ width, so $`m_W80.4`$ GeV, its default value. $`W`$ width effects are nevertheless fully included in HERWIG and in the other plots we show. At low $`q_T`$ this assumption does not change the results dramatically, at high $`q_T`$ the effect of the $`W`$ width is important as it allows values of $`q_T`$ larger that the default $`W`$ mass even in the parton shower approximation, otherwise they could come only via the exact matrix-element generated events. with the ones obtained from these approaches  all normalized to the corresponding total cross section. HERWIG clearly lies well within the range of the resummed approaches, except at very small $`q_T`$ where they become unreliable. To be more precise, the agreement is better between HERWIG and the two resummations in the $`q_T`$-space, but even the $`b`$-space result is not too far from the Monte Carlo distribution.

If we now wish to compare the HERWIG simulation with matrix-element corrections to the resummed calculations even for larger values of $`q_T`$ we need to match the latter with the exact $`𝒪(\alpha _S)`$ results to make them reliable there. This has already been done in the literature within the approach of , while in  the analysis is limited to the low $`q_T`$ regime.

Many prescriptions exist concerning how to perform such a matching. Ours is to simply add the exact matrix-element cross sections for the parton level processes $`q\overline{q}^{}Wg`$ and $`q(\overline{q})gWq^{}(\overline{q}^{})`$, already calculated in (16) and (17), to the resummed expressions and, in order to avoid double counting, subtract off the terms they have in common. It is straightforward to show that these are simply those terms in the exact $`𝒪(\alpha _S)`$ result that do not vanish as $`q_T0`$. This prescription works fine if the resulting distribution is continuous at the point $`q_T=m_W`$, which means that the resummation and the low $`q_T`$ $`𝒪(\alpha _S)`$ result exactly compensate each other and only the exact ‘hard’ matrix-element contribution survives.

As discussed in , it is not trivial to implement such a matching: the authors in fact do not succeed in obtaining a continuous distribution at the crucial matching point $`q_T=m_W`$, but rather a step of size $`\alpha _S^2`$ was found. This comes about because the derivative of the Sudakov exponent is not required to go smoothly to zero at that point.

We independently implement the matching for all the resummed calculations with which we wish to compare the HERWIG results and we find that the matching works well only for the $`q_T`$-space resummation performed in . For the $`b`$-space and the approaches in  we do indeed find a step at $`q_T=m_W`$. In fact even for the $`q_T`$-space method of  we find a ‘kink’ at $`q_T=m_W`$, i.e. although the curve is continuous, its derivative changes discontinuously there, albeit by an amount that is too small to notice on the figure.

In Fig. 7 we compare the HERWIG 6.1 distributions with  after the matching over the whole $`q_T`$ spectrum.

We see that for the resummed distribution that is well matched and continuous for $`q_T=m_W`$ the agreement with HERWIG is pretty good everywhere. We do have slight discrepancies for medium values of $`q_T`$, but they are well within the range that could be expected from the differences between the approaches followed by the Monte Carlo program and the calculation which keeps all the next-to-leading logarithms in the form factor.

While the plots shown so far refer to $`W`$ production, it is also worth comparing with some recent preliminary CDF data on $`Z`$ production . In Fig. 8, we have the CDF distribution with respect to the transverse momentum of the $`\gamma ^{}/Z`$ boson produced at the Tevatron and decaying into an $`e^+e^{}`$ pair with invariant mass in the range 66 GeV $`<m_Z<`$ 116 GeV.

We compare the data with HERWIG before and after matrix-element corrections; we also normalize the HERWIG distribution to the experimental value of the cross section, 245.3 pb. The result is that we obtain good agreement with the experimental data only thanks to the application of the hard and soft corrections, otherwise the predictions would have been badly wrong for values $`q_T>50`$ GeV. There is perhaps some evidence that HERWIG does not produce enough smearing at low $`q_T`$, even with an intrinsic $`p_t`$ of 1 GeV, with HERWIG peaking at about 2 GeV and the data peaking at about 3 GeV, but the overall fit is nevertheless acceptable.

We have however found that better agreement can be obtained with an intrinsic $`p_t`$ of 2 GeV.

### 5.2 Jet distributions

We now look at the impact the matrix-element corrections have on the jet activity at the Tevatron and at the LHC. An interesting object to analyse is the hardest jet in transverse energy (the so-called ‘first jet’). In Fig. 9 we plot the differential spectrum for the transverse energy of the first jet for $`\sqrt{s}=1.8`$ TeV and $`\sqrt{s}=14`$ TeV, using HERWIG 5.9 and 6.1 and running the inclusive version of the $`k_T`$ algorithm  for a radius $`R=0.5`$ at the Tevatron and $`R=1`$ for the LHC The correspondence between the radius $`R`$ in the $`k_T`$ algorithm and $`R_{\mathrm{cone}}`$ in an iterative cone algorithm is $`R_{\mathrm{cone}}0.75\times R`$ . The Tevatron experimentalists run an iterative cone algorithm with radius $`R_{\mathrm{cone}}=0.4`$, so we choose $`R=0.5`$ for the radius parameter when we consider jet events at $`\sqrt{s}=1.8`$ TeV. For the LHC we stick to the recommended value of $`R=1`$..

The result is that the improvement introduced does have a significant effect since the number of events in which the first jet has high $`E_T`$ is significantly increased. The 5.9 and 6.1 distributions are similar for small values of $`E_T`$, but for increasing $`E_T`$ the effect of the corrections introduced in HERWIG gets more and more relevant and at very high $`E_T`$ only events generated via the exact hard amplitude survive. The impact is really enormous in the case of the LHC as can be seen from Fig. 9b.

In Fig. 10 we plot the inclusive number of jets $`n_{\mathrm{jets}}`$ that pass a transverse energy cut $`E_T>10`$ GeV.

We see that implementing the matrix-element corrections significantly shifts the distribution towards larger $`n_{\mathrm{jets}}`$. If we look at events with three or four jets having $`E_T>10`$ GeV, we see that their number is increased considerably both at the Tevatron and at the LHC. We have roughly an enhancement of a factor of 2 for three-jet events at both the energies, while for events with four high transverse energy jets, at the LHC we still get an enhancement of 2, while at the Tevatron the difference is almost a factor of 4.

### 5.3 Rapidity distributions

Fig. 11 shows the distribution of the rapidity $`y`$ of the dilepton pair at $`\sqrt{s}=1.8`$ TeV and $`\sqrt{s}=14`$ TeV.

As we sum over $`W^+`$ and $`W^{}`$ they are symmetric in $`\pm y`$. We see that the matrix-element corrections do not significantly affect the rapidity distributions.

Fig. 12 shows the comparison between HERWIG and CDF  for the rapidity of the produced $`e^+e^{}`$ pair; the agreement is again good and the contribution of the matrix-element corrections is insignificant.

## 6 Conclusions

We have analysed Drell–Yan processes in hadron collisions in the Monte Carlo parton shower approach. This is accurate in the soft/collinear approximation, but leaves empty regions in phase space. We implemented matrix-element corrections generating radiation according to the first-order amplitude in the dead zone and for every hardest so far emission in the already populated region in the HERWIG parton shower. We compared our results with the previous version HERWIG 5.9, with experimental Tevatron data from the DØ and CDF collaborations and existing resummed calculations of the spectrum of the transverse momentum $`q_T`$ of the vector boson. We found that the implemented corrections have a marked impact on the phenomenological distributions for high values of $`q_T`$ and the new version of HERWIG fits the DØ data for $`W`$ production well over the whole $`q_T`$ spectrum after we correct the HERWIG results to detector level. At large $`q_T`$ it is crucial to provide the Monte Carlo algorithm with matrix-element corrections in order to succeed in obtaining such an agreement.

We also compared the HERWIG results after matrix-element corrections to some existing calculations based on a resummation of the initial-state radiation in the $`q_T`$-space and in the $`b`$-space. We found that in the range of low $`q_T`$, where actually the effect of matrix-element corrections is not so relevant and the resummed calculations are quite reliable, the parton shower distribution is in reasonable agreement with all of them, with discrepancies due to the methods followed by these different approaches. We also matched the resummed results to the exact $`𝒪(\alpha _S)`$ result, in such a way to allow them to be trustworthy at all $`q_T`$ values and found that the matching works well for a resummation performed in the $`q_T`$-space keeping all the next-to-leading logarithms in the Sudakov exponent. In this case, we also obtained good agreement with the HERWIG 6.1 $`q_T`$ distribution. The other approaches considered showed a discontinuity at the point $`q_T=m_W`$ once we match them to the exact first-order perturbative result.

We also studied $`W+`$ jet events at the Tevatron and at the LHC and found a significant effect of the new improvement of HERWIG as a larger number of jets of large transverse energy passes the typical experimental cuts.

We compared the new version of HERWIG with the experimental data of the CDF collaboration on the transverse momentum and rapidity of $`Z`$ bosons. We found good agreement after implementing the corrections. As a result, we feel confident that the simulation of vector boson production is now reliable. Using the new version of HERWIG 6.1 to fit the experimental data will therefore provide us with better tests of the Standard Model and of QCD for the following Run II at the Tevatron and, ultimately, at the LHC.

For the sake of completeness we have however to say that our analysis has been performed forcing the vector boson to decay into a lepton pair, as most of the experimental studies do. For hadronic channels (i.e. $`Wq\overline{q}^{}`$), also the decay products are allowed to emit gluons and to give rise to a parton shower that is still described in the leading soft/collinear approximation by HERWIG. The implementation of matrix-element corrections to hadronic $`W`$ decays is straightforward as they are very similar to the corrections to the process $`Zq\overline{q}`$ that are discussed in  and is in progress. It is also worth remarking that the method applied in this paper to implement matrix-element corrections to the initial-state radiation in $`W`$ and $`Z`$ production can be extended to a wide range of processes that are relevant for the phenomenology of hadron colliders. Among these, the inclusion of matrix-element corrections to simulations of heavy quark production and particularly of top production in $`p\overline{p}`$ or $`pp`$ interactions is expected to have a marked impact on the top mass reconstruction and many observables that are relevant for heavy quark phenomenology. This work is also in progress.

## Acknowledgements

We acknowledge Stefano Frixione, Michelangelo Mangano, Stefano Moretti, Willis Sakumoto and Bryan Webber for discussions of these and related topics. We are indebted to Giovanni Ridolfi who provided us with the code to obtain the plots in Fig. 6. We are also grateful to the DØ Collaboration for making their detector simulation, CMS, available to us, and especially to Cecilia Gerber for the considerable effort it took to make it run outside the usual DØ environment.

## References

1. G. Altarelli, R.K. Ellis and G. Martinelli, Nucl. Phys. B143 (1978) 521; Nucl. Phys. B146 (1978) 544 (erratum); Nucl. Phys. B157 (1979) 461.
2. Yu.L. Dokshitzer, D.I. Dyakonov and S.I. Troyan, Phys. Rep. 58 (1980) 269.
3. J. Collins and D. Soper, Nucl. Phys. B193 (1981) 381; Erratum Nucl. Phys. B213 (1983) 454; Nucl. Phys. B197 (1982) 446;   
J. Collins, D. Soper and G. Sterman, Nucl. Phys. B250 (1985) 199.
4. G.A. Ladinsky and C.P. Yuan, Phys. Rev. D50 (1994) 4239.
5. C.T.H. Davies, W.J. Stirling and B.R. Webber, Nucl. Phys. B256 (1985) 413.
6. P.B. Arnold and R. Kauffman, Nucl. Phys. B349 (1991) 381.
7. R.K. Ellis, D.A. Ross and S. Veseli, Nucl. Phys. B503 (1997) 309.
8. R.K. Ellis and S. Veseli, Nucl. Phys. B511 (1998) 649.
9. S. Frixione, P. Nason and G. Ridolfi, Nucl. Phys. B542 (1999) 311.
10. A. Kulesza and W.J. Stirling, DTP-99-02, hep-ph/9902234.
11. T. Sjöstrand, Comp. Phys. Comm. 46 (1987) 367.
12. G. Marchesini et al., Comput. Phys. Commun. 67 (1992) 465.
13. G. Miu and T. Sjöstrand, Phys. Lett. B449 (1999) 313.
14. S. Mrenna, UCD-99-13, hep-ph/9902471.
15. M.H. Seymour, Comput. Phys. Commun. 90 (1995) 95.
16. M.H. Seymour, Z. Phys. C56 (1992) 161.
17. M.H. Seymour, Matrix Element Corrections to Parton Shower Simulation of Deep Inelastic Scattering, contributed to 27th International Conference on High Energy Physics (ICHEP), Glasgow, 1994, Lund preprint LU-TP-94-12, unpublished.
18. G. Corcella and M.H. Seymour, Phys. Lett. B442 (1998) 417.
19. T. Sjöstrand, Phys. Lett. B157 (1985) 231.
20. G. Marchesini and B.R. Webber, Nucl. Phys. B310 (1988) 461.
21. S. Catani, G. Marchesini and B.R. Webber, Nucl. Phys. B349 (1991) 635.
22. M.H. Seymour, Nucl. Phys. B436 (1995) 443.
23. A.D. Martin, R.G. Roberts and W.J. Stirling, Phys. Lett. B387 (1996) 419.
24. DØ Collaboration, B. Abbott et al., Phys. Rev. Lett. 80 (1998) 5498.
25. DØ Collaboration, B. Abbott et al., Phys. Rev. Lett. 80 (1998) 3000;   
DØ Collaboration, B. Abbott et al., Phys. Rev. D58 (1998) 092003;   
I. Adam, Ph.D. thesis, Columbia University, 1997, Nevis Report #294,   
http://www-d0.fnal.gov/results/publications\_talks/thesis/adam/ian\_thesis\_all.ps;   
E. Flattum, Ph.D. thesis, Michigan State University, 1996,   
http://www-d0.fnal.gov/results/publications\_talks/thesis/flattum/eric\_thesis.ps.
26. CDF Collaboration, T. Affolder et al., Fermilab-Pub-99/220-E;   
CDF Collaboration, A. Bodek et al., Fermilab-Conf-99/160-E.
27. S. Catani, Yu.L. Dokshitzer, M.H. Seymour and B.R. Webber, Nucl. Phys. B406 (1993) 187   
M.H. Seymour, Nucl. Phys. B421 (1994) 545.
28. S.D. Ellis and D.E. Soper, Phys. Rev. D48 (1993) 3160.