---
author:
- |
    [^1]\
    H. Niewodniczański Institute of Nuclear Physics, Polish Academy of Sciences,\
    ul. Radzikowskiego 152, 31-342 Krakow, Poland\
    E-mail:
- |
    A. Kusina\
    Southern Methodist University, Dallas, TX 75275, USA\
    E-mail:
- |
    M. Skrzypek\
    H. Niewodniczański Institute of Nuclear Physics, Polish Academy of Sciences,\
    ul. Radzikowskiego 152, 31-342 Krakow, Poland\
    E-mail:
- |
    M. Slawinska\
    H. Niewodniczański Institute of Nuclear Physics, Polish Academy of Sciences,\
    ul. Radzikowskiego 152, 31-342 Krakow, Poland\
    E-mail:
title: |
    NLO parton shower for LHC physics -\
    hard processes and beyond
---

Introduction
============

The Large Hadron Collider (LHC) at CERN provides rich harvest of experimental data. The proper understanding and interpretation of these data, possibly leading to discovery of new phenomena, requires perfect mastering of the “trivial” effects due to the multiple emissions of soft and collinear gluons and quarks. Perturbative Quantum Chromodynamics (pQCD) [@GWP; @Gross:1974cs; @Georgi:1951sr], supplemented with clever modelling of the low energy nonperturbative effects, is an indispensable tool for disentangling the Standard Model physics component in the data. This work presents part of the global effort of improving quality of the pQCD calculations for LHC experiments.

Most of the results presented here are described in refs. [@Jadach:2011cr] and [@Jadach:2012vs]. Although this work elaborates on the improved method of the pQCD calculation combining NLO-corrected hard process and LO parton shower Monte Carlo (MC), it should be regarded as the first step towards NNLO-corrected hard process combined with the NLO parton shower MC [@IFJPAN-IV-2012-7].

Basic LO parton shower MC
=========================

The multigluon distribution of the single initial state ladder, which is a building block of our parton shower MC, is represented by the integrand of the “exclusive/unintegrated PDF”, which in the LO approximation is the following: $$\label{eq:LOMC}
\begin{split}
& 
D(t,x)
=\int dx_0\; dZ\; 
 \delta_{x=x_0 Z}\;
d_0(\hat{t}_0,x_0)\; G(t, \hat{t}_0- \ln x_0 | Z),
\\&
G(t, t_0 | Z)=
 e^{-S_F} \sum_{n=0}^\infty
\bigg(
\prod_{i=1}^n
 \int d^3{{\cal E}}(\bar{k}_i)\;
 \theta_{\xi_i>\xi_{i-1}}
 \frac{2C_F\alpha_s}{\pi^2} \bar{P}(z_i)
\bigg)
\\&
\qquad\qquad\quad\times
\theta_{t>\xi_n}
\delta_{Z=\prod_{j=1}^n z_j},
\end{split}$$ where evolution kernel is $\bar{P}(z)=\frac{1}{2}(1+z^2)$, evolution time is $\hat{t}_0=\ln(q_0/\Lambda)$ and the “eikonal” phase space integration element is $
d^3{{\cal E}}(k)=\frac{d^3 k}{2k^0}\;\frac{1}{{{\bf{k}}}^2}
            =\pi \frac{d\phi}{2\pi} \frac{d k^+}{k^+} d \xi
$ and $k^\pm = k^0\pm k^3$. We use rapidities $\xi_i =\frac{1}{2}\ln\frac{k^-_i}{k^+_i}\big|_{\rm Rh}$ in the hadron beam rest frame (Rh), and $\eta_i=\frac{1}{2}\ln\frac{k^+_i}{k^-_i}\big|_{\rm RFHP}$ defined in hard process rest frame (RFHP). They are related by $\xi_i=\ln\frac{\sqrt{s}}{m_h}-\eta_i$. Rapidity ordering is now $t=\xi_{\max}>\xi_n>\dots>\xi_i>\xi_{i-1}>\dots>\xi_0=t_0$, where $t_0=\xi_0=\ln(q_0/m_h)-\ln x_0$. The direction of the $z$ axis in the RFHP is pointing out towards the hadron momentum. A lightcone variable of the emitted gluon is defined as $\alpha_i= \frac{2k_i^+}{\sqrt{s}}$ and of the emitter parton (quark) as $x_i=x_0-\sum_{j=0}^{i}\; \alpha_j$ (after $i$ emissions). We also use fractions $z_i=x_i/x_{i-1}$. The Sudakov formfactor $S_F$ comes from the “unitarity” condition[^2] $
\int_0^1 dZ\; G(t, t_0 | Z)=1,
$ which is also instrumental in the Markovian MC implementation used to obtain $D(t,x)$ at any value of $t>t_0$.

The initial distribution $d_0(q_0,x_0)$ related to experiment, to previous steps in the MC ladder, or to PDF in the standard $\overline{MS}$ system is not essential for the following discussion, we only note that the unitarity condition provides baryon number conservation sum rule $\int_0^1 dx\; D(t,x) = \int_0^1 dx_0\; d_0(t_0,x_0)$.

For testing our new method of correcting hard process to the NLO level we use the following simplified MC parton shower, implementing the DY process with two ladders and the hard process: [^3] $$\label{eq:LOMCFBmaster}
\begin{split}
&\sigma_0=
\int d x_{0F} d x_{0B}\;\;
  d_0(\hat{t}_0,x_{0F}) d_0(\hat{t}_0,x_{0B}) 
 \sum_{n_1=0}^\infty\;
 \sum_{n_2=0}^\infty
 \int dx _F\; dx_B\;
\\&~~~~~~~~~\times
e^{-S_{_F}}
\int_{\Xi<\eta_{n_1}}
\bigg(
\prod_{i=1}^{n_1}
 d^3{{\cal E}}(\bar{k}_i)
 \theta_{\eta_i<\eta_{i-1}}
 \frac{2C_F\alpha_s}{\pi^2} \bar{P}(z_{Fi})
\bigg)
\delta_{x_F = x_{0F}\prod_{i=1}^{n_1} z_{Fi}}
\\&~~~~~~~~~\times
e^{-S_{_B}}
\int_{\Xi>\eta_{n_2}}
\bigg(
\prod_{j=1}^{n_2}
 d^3{{\cal E}}(\bar{k}_j)
 \theta_{\eta_j>\eta_{j-1}}
 \frac{2C_F\alpha_s}{\pi^2} \bar{P}(z_{Bj})
\bigg)
\delta_{x_B = x_{0B}\prod_{j=1}^{n_2} z_{Bj}}
\\&~~~~~~~~~\times
 d\tau_2(P-\sum_{j=1}^{n_1+n_2} k_j ;q_1,q_2)\;
\frac{d\sigma_B}{d\Omega}(sx_Fx_B,\hat\theta)\;
W^{NLO}_{MC}.
\end{split}$$ In the LO approximation $W^{NLO}_{MC}=1$. Rapidity $\xi$ is translated into $\eta$ – the center of mass system rapidity, in the forward part (F) of the phase space as $\xi_i=\ln\frac{\sqrt{s}}{m_h}-\eta_i$, $\eta_{0F}>\eta_i>\Xi$, and in the backward (B) part as $\xi_i=-\ln\frac{\sqrt{s}}{m_h}+\eta_i$, $\Xi>\eta_i>\eta_{0B}$. The rapidity boundary between the two hemispheres $\Xi=0$ is used, until a more sophisticated version related to rapidity of the produced $W/Z$ is introduced.

Analytical integration of eq. (\[eq:LOMCFBmaster\]) results in the standard factorization formula ($W^{NLO}_{MC}=1$) $$\label{eq:LOfactDY2LO}
\sigma_0 =
\int_0^1 dx_F\;dx_B\;
D_F(t, x_F)\; D_B(t, x_B)\;
\sigma_B(sx_Fx_B).$$ The distributions $D_F(t, x_F)=(d_0\otimes G_F) (t, x_F)$ and $D_B(t, x_B)=(d_0\otimes G_B) (t, x_B)$ are obtained from separate Markovian LO Monte Carlo runs. The above LO formula is exact, and can be tested with an arbitrary numerical precision.

Figure \[fig:etaW\_LO\_7TeV\] represents a “calibration benchmark” for the overall normalization at the LO level. We show there the properly normalized distribution of the variable $\eta_W^*=\frac{1}{2}\ln(x_F/x_B)$, which in the collinear limit approximates the rapidity of $W$ boson. The distribution in the upper plot of Fig. \[fig:etaW\_LO\_7TeV\], representing eq. (\[eq:LOfactDY2LO\]), is obtained using the general purpose MC program FOAM [@foam:2002]. The collinear PDF $D(t,x)$ there has been obtained from a separate high statistics run ($10^{10}$ events) of a Markovian MC (MMC), creating $D(t,x)$ in a form of the 2-dimensional look-up table[^4]. The other distribution in the upper plot of Fig. \[fig:etaW\_LO\_7TeV\] represents eq. (\[eq:LOMCFBmaster\]) in LO approximation. It comes from the full scale MC generation (with four-momenta conservation). The MC run with $10^8$ events was used. The constrained MC (CMC) technique of ref. [@Jadach:2007qa] is used here because of the narrow Breit-Wigner peak due to a heavy boson propagator[^5]. Two CMC modules and FOAM are combined into one MC generating gluon emissions and the $W$ boson production. FOAM is taking care of the generation of the variables $x_F,x_B,x_{F0},x_{B0}$ and the sharp Breit-Wigner peak in $\hat{s}=s x_F x_B$, then two CMC modules are initialized and generate the gluon four-momenta $\bar{k}^\mu_j$. They are mapped into $k^\mu_j$, following the prescription defined in ref. [@Jadach:2011cr], such that the overall energy-momentum conservation is achieved. Figure \[fig:etaW\_LO\_7TeV\] demonstrates a very good numerical agreement between $d\sigma/ d \eta_W^*$ from our full scale LO parton shower MC of eq. (\[eq:LOMCFBmaster\]) and the simple formula of eq. (\[eq:LOfactDY2LO\]), to within 0.5%, as seen from the ratio of the two results in the lower part of the figure.

Introducing NLO corrections to hard process
===========================================

The NLO corrections to hard process are imposed on top of the LO distributions of eq. (\[eq:LOMCFBmaster\]) using a single “monolithic” weight $W^{NLO}_{MC}$ defined exactly as in ref. [@Jadach:2011cr]: $$\label{eq:NLODYMCwt}
\begin{split}
&W^{NLO}_{MC}=
1+\Delta_{S+V}
+\sum_{j\in F} 
 \frac{{{\tilde{\beta}}}_1(q_1,q_2,\bar{k}_j)}      {\bar{P}(z_{Fj})\;d\sigma_B(\hat{s},\hat\theta)/d\Omega}
+\sum_{j\in B} 
 \frac{{{\tilde{\beta}}}_1(q_1,q_2,\bar{k}_j)}      {\bar{P}(z_{Bj})\;d\sigma_B(\hat{s},\hat\theta)/d\Omega},
\end{split}$$ the NLO soft+virtual correction is $
\Delta_{V+S}
=\frac{C_F \alpha_s}{\pi}\; \left( \frac{2}{3}\pi^2 -\frac{5}{4} \right)
$, and the real correction reads: $$\label{eq:DYbeta1FB}
\begin{split}
&{{\tilde{\beta}}}_1(q_1,q_2,k)=
\Big[
  \frac{(1-\beta)^2}{2}
  \frac{d\sigma_{B}}{d\Omega_q}(\hat{s},\theta_{F})
 +\frac{(1-\alpha)^2}{2}
  \frac{d\sigma_{B}}{d\Omega_q}(\hat{s},\theta_{B})
\Big]
\\&~~~~~~~~~~~~~~~~~~~~
-\theta_{\alpha>\beta}
 \frac{1+(1-\alpha-\beta)^2}{2}
 \frac{d\sigma_{B}}{d\Omega_q}(\hat{s},\hat\theta)
-\theta_{\alpha<\beta}
 \frac{1+(1-\alpha-\beta)^2}{2}
 \frac{d\sigma_{B}}{d\Omega_q}(\hat{s},\hat\theta).
\end{split}$$ The above is the exact ME of the quark-antiquark annihilation into a heavy vector boson with additional single real gluon emission[^6]. The LO component, which is already included in the LO MC, is subtracted here. The variable $\hat{s}=s x_F x_B = (q_1+q_2)^2$ is the effective mass squared of the heavy vector boson. The definition of angle $\hat\theta$ in the LO component is rather arbitrary. We define it in the rest frame of the heavy boson, where $\vec{q}_1+\vec{q}_2=0$, as an angle between the decay lepton momentum $\vec{q}_1$ and the difference of momenta of the incoming quark and antiquark $\hat\theta=\angle(\vec{q}_1,\vec{p}_{0F}-\vec{p}_{0B})$. On the other hand the two angles in the NLO ME are defined quite unambiguously as $\hat\theta_F=\angle(\vec{q}_1,-\vec{p}_{0B})$ and $\hat\theta_B=\angle(\vec{q}_1, \vec{p}_{0F})$. In the above we only need directions of the $\vec{p}_{0F}$ and $\vec{p}_{0B}$ vectors, which are the same as the directions of the hadron beams. The lightcone variables $\alpha_j$ and $\beta_j$ of the emitted gluon are defined in the F and B parts of the phase space as follows[^7]: $$\begin{split}
&\alpha_j=1-z_{Fj},\quad
\beta_j=  \alpha_j\; e^{2(\eta_j-\Xi)},\quad ~~
{\rm for}~~~ j\in F,
\\&
\beta_j=1-z_{Bj},\quad
\alpha_j=\beta_j\;   e^{-2(\eta_j-\Xi)},\quad
{\rm for}~~~ j\in B.
\end{split}$$

Again, the exact phase space integration of eq. (\[eq:LOMCFBmaster\]) including $W^{NLO}_{MC}$ of eq. (\[eq:NLODYMCwt\]) is feasible, and the resulting compact expression for the total cross section is obtained [@Jadach:2011cr]: $$\label{eq:DYanxch}
\begin{split}
\sigma_1 &=
\int_0^1 dx_F\;dx_B\; dz\;
D_F(t, x_F)\;D_B(t, x_B)\;
\sigma_B(szx_Fx_B)
\big\{
\delta_{z=1}(1+\Delta_{S+V})
+C_{2r}(z)
\big\},
\end{split}$$ where $
C_{2r}(z) =\frac{2C_F \alpha_s}{\pi}\; \left[ -\frac{1}{2}(1-z) \right].
$

Numerical test of NLO correction
--------------------------------

Figure \[fig:etaW\_NLOglu\_7TeV\] represents a principal [*proof of concept*]{} of our new methodology for implementing the NLO corrections to the hard process in the parton shower MC. The plotted NLO correction to the $\eta_W^*$ distribution[^8] comes from the parton shower MC with the NLO-corrected hard process according to eqs. (\[eq:LOMCFBmaster\]) and (\[eq:NLODYMCwt\]). Additionally we also plot there result of a simple collinear formula of eq. (\[eq:DYanxch\]), where two collinear PDFs are convoluted with the analytical coefficient function $C_{2r}(z)$ for the hard process. Both results coincide within the statistical error, see their ratio in the lower part of Fig. \[fig:etaW\_NLOglu\_7TeV\].

Technically, the inclusion of the NLO correction in our parton shower MC is rather straightforward, and is obtained by including $W^{NLO}_{MC}$ weight of eq. (\[eq:NLODYMCwt\]). MC is providing both LO and NLO-corrected results in a single run with weighted events. The NLO weight is strongly peaked near $W^{NLO}_{MC}=1$, positive, and without long-range tails. Its distribution is shown in Fig. \[fig:Canv1w\_7TeV\].

In all numerical results we have set $\Delta_{V+S}=0$, as it is completely unimportant for the presented analysis. The initial distributions $d_0(q_0,x_0)$ are defined in ref. [@Jadach:2012vs].

Simplification of the method and comparison with other methodologies
====================================================================

![ The inclusive distribution of gluons on the log Sudakov plane of rapidity $t=\xi$ and $v=\ln(1-z)$ (LHS). Contributions from all gluons weighted with the component weight $W^{NLO}_j$ (RHS). ](mcCanv2_Rho1gluLO.jpg "fig:"){width="75mm"} ![ The inclusive distribution of gluons on the log Sudakov plane of rapidity $t=\xi$ and $v=\ln(1-z)$ (LHS). Contributions from all gluons weighted with the component weight $W^{NLO}_j$ (RHS). ](mcCanv2b_RhoDelNLO.jpg "fig:"){width="75mm"} \[fig:mcRho1glu\]

Our new method for introducing NLO corrections in the hard process, proposed in ref. [@Jadach:2011cr] and tested in ref. [@Jadach:2012vs], is an alternative to the two well established MC@NLO [@Frixione:2002ik] and POWHEG [@Nason:2004rx; @Frixione:2007vw] methodologies. With MC numerical implementation at hand, let us elaborate on the differences with the above two techniques in particular with the POWHEG technique. We shall also see that it is possible to make our method more efficient in terms of CPU time consumption. This improvement is not so critical in the present case of NLO corrected hard process, but may be quite useful in the case of correcting evolution kernels to the NLO in the ladder parts of the MC [@IFJPAN-IV-2012-7].

The most important differences with the POWHEG and MC@NLO techniques are:

-   The summation over all emitted gluons, without deciding which gluon is the one involved in the NLO correction and which ones are merely “LO spectators” in the parton shower.

-   The absence of $(1/(1-z))_+$ distributions in the real part of the NLO corrections (virtual+soft correction is kinematically independent).

To explain more clearly how $W^{NLO}_{MC}$ of eq. (\[eq:NLODYMCwt\]) is distributed over the multigluon phase space, we restrict now to single ladder (hemisphere) with a simplified weight: $$\label{eq:NLODYMCwt_sim}
\begin{split}
&W^{NLO}_{MC}= 1+\sum_{j\in F} W^{NLO}_j,\qquad
W^{NLO}_j=
 \frac{{{\tilde{\beta}}}_1(q_1,q_2,\bar{k}_j)}      {\bar{P}(z_{Fj})\;d\sigma_B(\hat{s},\hat\theta)/d\Omega}.
\end{split}$$ In order to find out the phase space regions specific for NLO corrections we consider inclusive distributions of gluons on the Sudakov logarithmic plane of rapidity $\xi$ and variable $v=\ln(1-z)$. In the left hand side (LHS) of Fig. \[fig:mcRho1glu\] we show gluons inclusive distribution in the LO approximation. The flat plateau there represents IR/collinear singularity[^9] $2C_F\frac{\alpha_S}{\pi} d\xi\frac{dz}{1-z}$ with the drop by factor 1/2 towards $z=0$, due to $\frac{1+z^2}{2}$ factor in the LO kernel. In the right hand side (RHS) of Fig. \[fig:mcRho1glu\] we show contributions from all gluons weighted with the component weight[^10] $-W^{NLO}_j$ of eq. (\[eq:NLODYMCwt\_sim\]). The NLO contribution is concentrated in the area near the hard process rapidity $t=\xi_{\max}$, which has to be true for the genuine NLO contribution [^11]. The completeness of the phase space near this important region ($z=0$, $\xi_{\max}$) is critical for the completeness of the NLO corrections. Both POWHEG and MC@NLO use standard LO MCs which feature an empty “dead zone” in this phase space corner.

Figure \[fig:mcRho1glu\] suggests that the dominant contribution to $\sum_j W^{NLO}_j$ could be from the gluon with the maximum $\ln k_j^T\sim \xi_j+\ln(1-z_j)$, which is closest to the hard process phase space corner. In the MC we may easily relabel generated gluons using new index $K$ such that they are ordered in the variable $\kappa_K=\xi_K+\ln(1-z_K),\; \kappa_{K+1}<\kappa_K$ with $K=1$ being the hardest one.

Figure \[fig:mcCanv1k\] demonstrates a split of the LO inclusive distribution of Fig. \[fig:mcRho1glu\] into the $K=1$ component and the rest $K>1$. The important point is that the $K=1$ component reproduces the original complete distribution over the whole region where the NLO correction is non-negligible! This is exactly the observation on which POWHEG technique is built. According to the POWHEG authors, taking the $K=1$ component is sufficient to reproduce the complete NLO correction (modulo NNLO).

The above statement is checked numerically in Fig. \[fig:mcCanV8k\], where we compare the NLO correction to the $x=\prod_j z_j$ distribution from the complete sum $\sum_j W^{NLO}_j$ and from $ W^{NLO}_{K=1}$. As we see the $K=1$ component saturates the complete sum very well, with the $K=2$ component being negligible in the first approximation.

We can therefore speed up the calculation by means of taking only the $K=1$ contribution. The price will be that the formula of eq. (\[eq:DYanxch\]) will not be exact any more. Our method differs, however, from the POWHEG scheme, where the $K=1$ gluon is generated separately in the first step, and other gluons are generated (by the LO parton shower MC) in the next step. That is easy for LO MC with $k^T$-ordering, while in case of the LO MC with angular-ordering POWHEG requires additional effort of generating the so called vetoed and truncated showers. In our method, there is no need for such vetoed/truncated showers in case of angular ordering.

The reason why POWHEG technique is complicated in case of the angular ordering is illustrated in Fig. \[fig:mcRhoOrdLO\]. We show there the distribution of gluons ordered in rapidity, starting from the gluon with the maximum rapidity, the closest to hard process. The gluon distribution with the highest rapidity $\xi\sim \xi_{\max}$ ($J=1$) has a ridge extending towards the soft region. Notice that, when the IR cut-off $\epsilon\to0$ in $(1-z)<\epsilon$, the width of this ridge also goes to zero. Consequently, the gluon with the highest $\xi$ is unable to reproduce the gluon distribution in the NLO corner, close to hard process. This is why in this case POWHEG requires truncated and vetoed showers, which are not needed in our method.

Summary and outlook
===================

A new method of adding the QCD NLO corrections to the hard process in the initial state Monte Carlo parton shower is tested numerically showing that the basic concept of the new methodology works correctly in the numerical environment of a Monte Carlo parton shower. The differences with the well established methods of MC@NLO and POWHEG are briefly discussed. Also, variants of the new method with better efficiency in terms of CPU time are proposed.

Acknowledgement {#acknowledgement .unnumbered}
===============

This work is partly supported by the Polish National Science Centre grant UMO-2012/04/M/ST2/00240, Foundation for Polish Science grant Homing Plus/2010-2/6, the Research Executive Agency (REA) of the European Union Grant PITN-GA-2010-264564 (LHCPhenoNet), the U.S. Department of Energy under grant DE-FG02-04ER41299 and the Lightner-Sams Foundation.

[99]{}

D. J. Gross and F. Wilczek, [*Phys. Rev. Lett.*]{} [**30**]{} (1973) 1343;\
H. D. Politzer, [*Phys. Rev. Lett.*]{} [**30**]{} (1973) 1346;\
D. J. Gross and F. Wilczek, [*Phys. Rev.*]{} [**D8**]{} (1973) 3633;\
H. D. Politzer, [ *Phys. Rep.*]{} [**14**]{} (1974) 129.

D. J. Gross and F. Wilczek, [*[Asymptotically Free Gauge Theories. 2]{}*]{}, [*Phys. Rev.*]{} [**D9**]{} (1974) 980–993.

H. Georgi and H. D. Politzer, [*[Electroproduction scaling in an asymptotically free theory of strong interactions]{}*]{}, [*Phys. Rev.*]{} [ **D9**]{} (1974) 416–420.

S. Jadach, A. Kusina, W. Placzek, M. Skrzypek, and M. Slawinska, [*[On the inclusion of the QCD NLO corrections in the quark– gluon Monte Carlo shower]{}*]{}, [[1103.5015]{}](http://xxx.lanl.gov/abs/1103.5015).

S. Jadach, M. Jezabek, A. Kusina, W. Placzek, and M. Skrzypek, [*[NLO corrections to hard process in QCD shower – proof of concept]{}*]{}, [[1209.4291]{}](http://xxx.lanl.gov/abs/1209.4291).

S. Jadach, A. Kusina, and Skrzypek, [*[NLO corrections to ladder part of the initial state shower in QCD]{}*]{}, 2012. Report IFJPAN-IV-2012-7, in preparation.

S. Jadach, [*[Foam: A general purpose cellular Monte Carlo event generator]{}*]{}, [*Comput. Phys. Commun.*]{} [**152**]{} (2003) 55–100, \[[[physics/0203033]{}](http://xxx.lanl.gov/abs/physics/0203033)\].

S. Jadach, W. Placzek, M. Skrzypek, and P. Stoklosa, [*[Markovian Monte Carlo program EvolFMC v.2 for solving QCD evolution equations]{}*]{}, [ *Comput. Phys. Commun.*]{} [**181**]{} (2010) 393–412, \[[[0812.3299]{}](http://xxx.lanl.gov/abs/0812.3299)\].

K. Golec-Biernat, S. Jadach, W. Płaczek, and M. Skrzypek, [*[Markovian Monte Carlo solutions of the NLO QCD evolution equations]{}*]{}, [*Acta Phys. Polon.*]{} [**B37**]{} (2006) 1785–1832, \[[[hep-ph/0603031]{}](http://xxx.lanl.gov/abs/hep-ph/0603031)\].

S. Jadach, W. Placzek, M. Skrzypek, P. Stephens, and Z. Was, [*[Constrained MC for QCD evolution with rapidity ordering and minimum kT\*]{}*]{}, [ *Comput.Phys.Commun.*]{} [**180**]{} (2009) 675–698, \[[[hep-ph/0703281]{}](http://xxx.lanl.gov/abs/hep-ph/0703281)\].

T. Sjostrand, [*A model for initial state parton showers*]{}, [*Phys. Lett.*]{} [**B157**]{} (1985) 321.

F. A. Berends and R. Kleiss, [*[Initial State Radiation for e+ e- Annihilation Into Jets]{}*]{}, [*Nucl.Phys.*]{} [**B178**]{} (1981) 141.

S. Alioli, K. Hamilton, and E. Re, [*[Practical improvements and merging of POWHEG simulations for vector boson production]{}*]{}, [*JHEP*]{} [**1109**]{} (2011) 104, \[[[1108.0909]{}](http://xxx.lanl.gov/abs/1108.0909)\].

S. Frixione and B. R. Webber, [*[Matching NLO QCD computations and parton shower simulations]{}*]{}, [*JHEP*]{} [**06**]{} (2002) 029, \[[[hep-ph/0204244]{}](http://xxx.lanl.gov/abs/hep-ph/0204244)\].

P. Nason, [*[A new method for combining NLO QCD with shower Monte Carlo algorithms]{}*]{}, [*JHEP*]{} [**11**]{} (2004) 040, \[[[hep-ph/0409146]{}](http://xxx.lanl.gov/abs/hep-ph/0409146)\].

S. Frixione, P. Nason, and C. Oleari, [*[Matching NLO QCD computations with Parton Shower simulations: the POWHEG method]{}*]{}, [*JHEP*]{} [**0711**]{} (2007) 070, \[[[0709.2092]{}](http://xxx.lanl.gov/abs/0709.2092)\].

[^1]: The partial support of the TH Unit of the CERN PH Division for this author is acknowledged.

[^2]: The usual cutoff $1-z<\epsilon$ regularizing the IR singularity is implicit.

[^3]: Following ref. [@Jadach:2011cr], we adopt $d\tau_2(P;q_1,q_2) =
       \delta^{(4)}(P-q_1-q_2)\frac{d^3q_1}{2q_1^0}\frac{d^3q_2}{2q_2^0}$.

[^4]: This MMC run solves the LO DGLAP equation using the MC method, as in refs. [@Jadach:2008nu; @GolecBiernat:2006xw].

[^5]: A backward evolution algorithm of ref. [@Sjostrand:1985xi] could be also used here.

[^6]: We employ here the compact representation of ref. [@Berends:1980jk], which has also been used in POWHEG [@Alioli:2011nr].

[^7]: See ref. [@Jadach:2011cr] for more explanations.

[^8]: Extra minus sign introduced to facilitate visualization.

[^9]: We use constant $\alpha_S$.

[^10]: We again insert a minus sign in order to facilitate visualization.

[^11]: It also vanishes towards the soft limit $z\to 1$.
