\subsection{Multiple photon effects and matching with NLO corrections
\label{MP}}
\subsubsection{Universal methods for leading logarithmic corrections
\label{HO}}
From inspection of Eqs.~(\ref{bsv}) and (\ref{bsvgg}) for the SV NLO QED
corrections to the cross section of the Bhabha scattering and $e^+ e^- \to \gamma\gamma$
process, it can be seen that large logarithms $L=\ln(s/m_e^2)$, due to
collinear photon emission, are present. Similar large logarithmic terms arise after
integration of the hard photon contributions from the kinematical domains of photon
emission at small angles with respect to charged particles. For the energy range of
meson factories
the logarithm is large numerically, i.e. $L \sim 15$ at the $\mathrm{\phi}$ factories and
$L \sim 20$ at the $B$ factories, and the corresponding terms give
the bulk of the total radiative correction. These contributions represent also the dominant
part of the NNLO effects discussed in Section \ref{NNLO}. Therefore, to achieve the required theoretical accuracy, the logarithmically enhanced contributions due to emission of soft and collinear photons must be taken into account at all orders in perturbation theory.
The methods for the calculation of higher-order (HO) QED corrections on the
basis of the generators employed nowadays at flavour factories were already
widely and successfully used in the 90s at LEP/SLC for electroweak tests of the SM.
They were adopted for the calculation of both the small-angle
Bhabha scattering cross section (necessary for the high-precision luminosity
measurement) and $Z$-boson observables. Hence, the theory accounting for
the control of HO QED corrections at meson factories can be considered
particularly robust, having passed the very stringent tests of the LEP/SLC era.
The most popular and standard methods to keep multiple photon effects under control are the QED Structure Function (SF)
approach \cite{Kuraev:1985hb,Altarelli:1986kq,Nicrosini:1986sm,Nicrosini:1987sw}
and Yennie-Frautschi-Suura (YFS) exponentiation \cite{Yennie:1961ad}.
The former is used in all the versions of the generator BabaYaga
\cite{CarloniCalame:2000pz,CarloniCalame:2003yt,Balossini:2006wc} and
MCGPJ \cite{Arbuzov:2005pt}
(albeit according to different realisations), while the latter is the theoretical recipe adopted
in BHWIDE \cite{Jadach:1995nk}. Actually, analytical QED SFs $D(x,Q^2)$, valid in the strictly collinear approximation,
are implemented in MCGPJ, whereas BabaYaga is based on
a MC Parton Shower (PS) algorithm to reconstruct $D(x,Q^2)$ numerically. \\
\noindent {\em{The Structure Function approach \label{SF}}} \\
\vspace*{-2mm}
Let us consider the annihilation process $e^- e^+ \to X$, where $X$ is some given final
state and $\sigma_0 (s)$ its LO cross section.
Initial-state (IS) QED radiative corrections can be described according to the following
picture. Before arriving at the annihilation point, the incoming electron (positron) of
four-momentum $p_{-(+)}$ radiates real and virtual photons. These
photons, due to the dynamical features of QED, are mainly radiated along the direction of
motion of the radiating particles, and their effect is mainly to reduce
the original four-momentum of the incoming electron (positron) to $ x_{1(2)} p_{-(+)}$.
After this pre-emission, the hard scattering process $e^- (x_1 p_-) e^+ (x_2 p_+) \to X$
takes place, at a reduced squared c.m. energy ${\hat s} = x_1 x_2 s$. The
resulting cross section, corrected for IS QED radiation, can be represented in the form~\cite{Kuraev:1985hb,Altarelli:1986kq,Nicrosini:1986sm}
\begin{equation}
\sigma (s) = \int_0^1 {\rm d} x_1 {\rm d} x_2 D(x_1, s) D(x_2, s) \sigma_0 (x_1 x_2 s)
\Theta({\rm cuts}),
\label{eq:masterd}
\end{equation}
where $D(x,s)$ is the electron SF, representing the probability that an
incoming electron (positron) radiates a collinear photon, retaining a fraction $x$ of its
original momentum at the energy scale $Q^2 = s$, and $\Theta({\rm cuts})$
stands for a rejection algorithm taking care of experimental cuts. When considering photonic radiation
only the non-singlet part of the SF is of interest. If the running of the QED
coupling constant is neglected, the non-singlet part of the SF is the solution
of the following Renormalisation Group (RG) equation, analogous to the
Dokshitzer-Gribov-Li\-pa\-tov-Altarelli-Parisi (DGLAP) equation of QCD
\cite{Gribov:1972rt,Altarelli:1977zs,Dokshitzer:1977sg}:
\begin{equation}
s {{\partial} \over {\partial s}} D(x,s) = {{\alpha} \over {2 \pi}} \int_x^1 {{{\rm d}z} \over {z}}
P_+(z) D \left( {{x} \over {z}} , s \right) ,
\label{eq:apeqdiff}
\end{equation}
where $P_+(z) $ is the regularised Altarelli-Parisi (AP) splitting function for the process ${\rm electron} \to {\rm electron} + {\rm photon}$, given by
\begin{eqnarray}
&& P_+(z) = P(z) - \delta (1-z) \int_0^1 {\rm d}x P(x) , \nonumber \\
&& P(z) = {{1 + z^2} \over {1-z}} .
\label{eq:apvert}
\end{eqnarray}
Equation (\ref{eq:apeqdiff}) can be also transformed into an integral equation, subject to the
boundary condition $D(x, m_e^2) = \delta(1-x)$:
\begin{equation}
D(x,s) = \delta (1-x) + {{\alpha} \over {2 \pi}} \int_{m_e^2}^{s} {{{\rm d} Q^2} \over {Q^2}}
\int_x^1 {{{\rm d}z} \over {z}} P_+(z) D \left( {{x} \over {z}} , Q^2 \right).
\label{eq:apeq}
\end{equation}
Equation~(\ref{eq:apeq}) can be solved exactly by means of numerical methods, such as the inverse
Mellin transform method. However, this derivation of $D(x,s)$ turns out be problematic
in view of phenomenological applications. Therefore, approximate (but very accurate)
analytical representations of the solution of the evolution equation are of major interest
for practical purposes. This type of solution was the one typically adopted in the context
of LEP/SLC phenomenology. A first analytical solution can be obtained in the
soft photon approximation, i.e. in the limit $x \simeq 1$. This solution, also known as
Gribov-Lipatov (GL) approximation, exponentiates the
large logarithmic contributions of infrared and collinear origin at all perturbative orders, but
it does not take into account hard-photon (collinear) effects. This drawback can be
overcome by solving the evolution equation iteratively. At the $n$-th step of the iteration,
one obtains the $O (\alpha^n)$ contribution to the SF for any value
of $x$. By combining the GL solution with the iterative one, in which
the soft-photon part has been eliminated in order to avoid
double counting, one can build a hybrid solution of the evolution equation. It exploits
all the positive features of the two kinds of solutions and is not affected by the limitations
intrinsic to each of them. Two classes of hybrid solutions, namely
the additive and factorised ones, are known in the literature, and both were adopted
for applications to LEP/SLC precision physics. A typical additive solution, where
the GL approximation $D_{GL} (x,s)$ is supplemented by finite-order terms present in the
iterative solution, is given by \cite{Cacciari:1992pz}
\begin{eqnarray}
&& D_{A} (x,s) = \sum_{i=0}^3 d_A^{(i)} (x,s) , \nonumber \\
&& d_A^{(0)} (x,s) = {{\exp \left[ {{{1} \over {2}} \beta \left( {{3} \over {4}} - \gamma_E \right)} \right] }
\over { \Gamma \left( 1 + {{1} \over {2}} \beta \right) }} {{1} \over {2}} \beta (1 - x)^{{{1} \over {2}} \beta - 1}, \nonumber \\
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%G \, {\eta \over 2} (1 - x)^{{\eta \over 2 } - 1} , \nonumber \\
&& d_A^{(1)} (x,s) = - {1 \over 4 } \beta (1+x) , \nonumber \\
&& d_A^{(2)} (x,s) = {1 \over 32 } \beta^2 \left[ (1+x) \left( -4 \ln (1-x) + 3 \ln x
\right) \right. \nonumber \\
&& \qquad \qquad \left. - 4 {{\ln x } \over {1-x}} - \, 5 - x \right] , \nonumber \\
&& d_A^{(3)} (x,s) = {1 \over 384 } \beta^3 \left\{ (1+x) \left[ 18 \zeta (2) - 6 \hbox{\rm Li}_2 (x)
\right. \right. \nonumber \\
&& \qquad \qquad \left. -12 \ln^2 (1-x) \right] + {1 \over {1-x}} \left[ - {3 \over 2} (1 + 8 x + 3 x^2) \ln x
\right. \nonumber \\
&& \qquad \qquad + {1 \over 2} (1 + 7 x^2 ) \ln^2 x -12 (1+ x^2) \ln x \ln (1-x) \nonumber \\
&& \qquad \qquad \left. \left.
- 6 (x+5) (1-x) \ln (1-x) \right. \right. \nonumber \\
&& \qquad \qquad \left. \left. - {1 \over 4} (39 - 24 x - 15 x^2 ) \right] \right\} ,
\label{eq:strucfuna}
\end{eqnarray}
where $\Gamma$ is the Euler gamma-function, $\gamma_E
\approx 0.5772$ the Euler-Mascheroni constant,
$\zeta $ the Riemann $\zeta$-function and $\beta$ is the large collinear factor
\begin{eqnarray}
\beta = {{2 \alpha} \over {\pi}} \left[ \ln \left( {{s} \over {m_e^2}} \right)
- 1 \right] .
\end{eqnarray}
Explicit examples of factorised solutions, which are obtained by multiplying the GL solution
by finite-order terms in such a way that, order by order, the iterative contributions are exactly recovered,
can be found in \cite{Skrzypek:1992vk}. For the calculation of HO corrections with a per mill accuracy
analytical SFs in additive and factorised form containing up to $O (\alpha^3)$
finite-order terms are sufficient and in excellent agreement. They also agree
with an accuracy much better than 0.1 with the exact numerical solution of the QED evolution equation.
Explicit solutions %for the QED LLA case in the collinear approximation have been found
up to the fifth order in $\alpha$ were calculated in~\cite{Przybycien:1992qe,Arbuzov:1999cq}.
%%%%%%%%% Andrej
The RG method described above was applied in~\cite{Arbuzov:1997pj}
for the treatment of LL QED radiative corrections to various processes of
interest for physics at meson factories. Such a formulation
was later implemented in the generator MCGPJ.
For example, according to~\cite{Arbuzov:1997pj}, the Bhabha scattering cross section,
accounting for LL terms in all orders, $O (\alpha^nL^n),\ n=1,2,\ldots$,
of perturbation theory, is given by
\begin{eqnarray} \label{master}
&& \dd \sigma^{\mathrm{Bhabha}}_{\mathrm{LLA}} = \!\!\!\!
\sum\limits_{a,b,c,d=e^\pm,\gamma}^{} \int^{1}_{\bar{z}_1}\dd z_1
\int^{1}_{\bar{z}_2} \dd z_2
D^{\mathrm{str}}_{ae^-} (z_1) D^{\mathrm{str}}_{be^+} (z_2)
\nonumber \\
&& \quad \times \dd \sigma_{0}^{ab\to cd}(z_1,z_2)
\int^{1}_{\bar{y}_1}\! \frac{\dd y_1}{Y_1}D^{\mathrm{frg}}_{e^-c} (\frac{y_1}{Y_1})
\int^{1}_{\bar{y}_2}\!\! \frac{\dd y_2}{Y_2}D^{\mathrm{frg}}_{e^+d} (\frac{y_2}{Y_2})
\nonumber \\
&& \quad + \, O \left(\alpha^2L,\alpha\frac{m_e^2}{s}\right)\,.
\label{eq:sfmcgpj}
\end{eqnarray}
Here $\dd\sigma_{0}^{ab\to cd}(z_1,z_2)$ is the differential LO cross
section of the process $ab\to cd$, with energy fractions of the incoming particles
being scaled by factors $z_1$ and $z_2$ with respect to the initial
electron and positron, respectively. In the notation of~\cite{Arbuzov:1997pj},
the electron SF $D^{\mathrm{str}}_{ab}(z)$ is distinguished from the
electron fragmentation function $D^{\mathrm{frg}}_{ab}(z)$ to point out the
role played by IS radiation (described by $D^{\mathrm{str}}_{ab}(z)$)
with respect to the one due to final-state radiation (described by
$D^{\mathrm{frg}}_{ab}(z)$). However, because of their probabilistic
meaning, the electron structure and fragmentation functions coincide.
In Eq.~(\ref{eq:sfmcgpj}) the quantities $Y_{1,2}$ are the energy fractions of particles
$c$ and $d$ with respect to the beam energy. Explicit expressions for
$Y_{1,2}=Y_{1,2}(z_1,z_2,\cos\theta)$ and other details on the kinematics can
be found in~\cite{Arbuzov:1997pj}.
The lower limits of the integrals, $\bar{z}_{1,2}$ and
$\bar{y}_{1,2}$, should be defined according to the experimental conditions of
particle detection and kinematical constraints. For the case of
the $e^+e^-\to \gamma\gamma$ process
one has to change the master formula~(\ref{master}) by picking up the two-photon final state.
Formally this can be done by just choosing the proper fragmentation functions,
$D^{\mathrm{frg}}_{\gamma c}$ and $D^{\mathrm{frg}}_{\gamma d}$.
The photonic part of the non-singlet electron structure (fragmentation)
function in $O (\alpha^nL^n)$ considered in~\cite{Arbuzov:1997pj} reads
\begin{eqnarray} \nonumber
&& D_{ee}^{NS,\gamma}(z)=\delta(1-z)+\sum_{i=1}^{n}\left(\frac{\alpha}{2\pi}(L-1)\right)^i
\frac{1}{i!}\left[P^{(0)}_{ee}(z)\right]^{\otimes i} \!\!\!\! ,
\\ \nonumber
&& D_{\gamma e}(z)= \frac{\alpha}{2\pi}(L-1)P_{\gamma e}(z) + O (\alpha^2L^2),
\\ \nonumber
&& D_{e \gamma}(z)= \frac{\alpha}{2\pi}L P_{e\gamma}(z) + O (\alpha^2L^2),
\\ \nonumber
&& P^{(0)}_{ee}(z)=\left[\frac{1+z^2}{1-z}\right]_{+}
\\ \nonumber
&&\quad = \lim\limits_{\Delta \to 0}\left\{\delta(1-z)(2\ln\Delta
+ \frac{3}{2}) + \Theta(1-z-\Delta)\frac{1+z^2}{1-z}\right\},
\\
&& \left[P^{(0)}_{ee}(z)\right]^{\otimes i} = \int\limits_{z}^{1}\frac{\dd t}{t}
P^{(i-1)}_{ee}(t)P^{(0)}_{ee}
\left(\frac{z}{t}\right),
\\ \nonumber
&& P_{\gamma e}(z) = z^2 + (1-z)^2, \quad
P_{e\gamma}(z) = \frac{1+(1-z)^2}{z}\, .
\end{eqnarray}
Starting from the second order in $\alpha$ there appear also non-singlet and singlet
$e^+e^-$ pair contributions to the structure function:
\begin{eqnarray}
&& D_{ee}^{NS,e^+e^-}(z) = \frac{1}{3}\left(\frac{\alpha}{2\pi}L\right)^2 P^{(1)}_{ee}(z)
+ O (\alpha^3L^3),
\nonumber \\ \nonumber
&& D_{ee}^{S,e^+e^-}(z) = \frac{1}{2!}\left(\frac{\alpha}{2\pi}L\right)^2R(z)
+ O (\alpha^3L^3),
\\ \nonumber
&& R(z) = P_{e \gamma}\otimes P_{\gamma e}(z)
= \frac{1-z}{3z}(4+7z+4z^2) \nonumber \\
&& \qquad \qquad \qquad \qquad \quad + 2(1+z)\ln z.
\end{eqnarray}
Note that radiation of a real pair, {\it i.e.} appearance of additional
electrons and positrons in the final state, require the application of nontrivial
conditions of experimental particle registration. Unambiguously, that can be done
only within a MC event generator based on
four-particle matrix elements,
as already discussed in Section \ref{NNLO}.
In the same way as in QCD, the LL cross sections depend on the choice of
the factorisation scale $Q^2$ in the argument of the large logarithm $L=\ln(Q^2/m_e^2)$,
which is not fixed a priori by the theory. However, the scale should be taken of the order
of the characteristic energy transfer in the
process under consideration. Typical choices are $Q^2=s$, $Q^2=-t$ and $Q^2=s t/u$.
The first one is good for annihilation channels like $e^+e^-\to \mu^+\mu^-$, the
second one is optimal for small-angle Bhabha scattering where the $t$-channel exchange
dominates, see~\cite{Arbuzov:2006mu}.
The last choice allows to exponentiate the leading contribution due to initial-final state
interference \cite{Greco:1990on} and is particularly suited for large-angle Bhabha scattering
in QED. The option $Q^2=s t/u$ is adopted in all the versions of the generator
BabaYaga. Reduction of the scale dependence can be achieved by taking into account
next-to-leading corrections in $O (\alpha^nL^{n-1})$, next-to-next-to-leading
ones in $O (\alpha^nL^{n-2})$ {\it etc.} \\
\noindent {\em{The Parton Shower algorithm \label{PS}}} \\
\vspace*{-2mm}
The PS algorithm is a method for providing a MC iterative
solution of the evolution equation and, at the same time, for generating the four-momenta of the
electron and photon at a given step of the iteration. It was developed
within the context of QCD %(see for instance refs.~\cite{o80} and~\cite{mw84}) and
and later applied in QED too.
%(see refs.~\cite{km89}, \cite{fsm93}, \cite{m95} and references therein).
%This approach is implemented in all the versions of the event generator BabaYaga
%to account for multiple photon emission effects.
In order to implement the algorithm, it
is first necessary to assume the existence of an upper limit for the energy fraction
$x$ in such a way that the AP splitting function is regularised by writing
\begin{equation}
P_+(z) = \theta(x_+ - z) P(z) - \delta (1-z) \int_0^{x_+} {\rm d}x P(x) .
\label{eq:apvertmod}
\end{equation}
Of course, in the limit $x_+ \to 1$, Eq.~(\ref{eq:apvertmod}) recovers the usual definition of
the AP splitting function given in Eq.~(\ref{eq:apvert}). By inserting the modified AP vertex into
Eq.~(\ref{eq:apeqdiff}), one obtains
\begin{eqnarray}
s {{\partial} \over {\partial s}} D(x,s) = && {{\alpha} \over {2 \pi}} \int_x^{x_+}
{{{\rm d}z} \over {z}} P(z) D \left( {{x} \over {z}} , s \right) \nonumber\\
&& - {{\alpha} \over {2 \pi}} D(x,s) \int_x^{x_+} {\rm d}z P(z) .
\label{eq:apeqmod}
\end{eqnarray}
Separating the variables and introducing the Sudakov form factor
\begin{equation}
\Pi (s_1, s_2) = \exp \left[
- {{\alpha} \over {2 \pi}} \int_{s_2}^{s_1} {{{\rm d} s'} \over {s'}} \int_0^{x_+} {\rm d}z P(z)
\right] ,
\label{eq:sudakov}
\end{equation}
which is the probability that the electron evolves from virtuality $-s_2$ to $-s_1$ without
emitting photons of energy fraction larger than $1-x_+ \equiv \epsilon$ ($\epsilon \ll 1$),
Eq.~(\ref{eq:apeqmod}) can be recast into the integral form
\begin{eqnarray}
D(x,s) && = \Pi (s,m_e^2) D(x, m_e^2) \nonumber\\
&& \quad \! \! \! \! \! \! + {{\alpha} \over {2 \pi}} \int_{m_e^2}^{s}
{{{\rm d}s'} \over {s'}} \Pi (s,s') \int_{x}^{x_+} {{{\rm d}z} \over {z}} P(z)
D \left( {{x} \over {z}} , s' \right) . \nonumber\\\
\label{eq:apeqmodi}
\end{eqnarray}
The formal iterative solution of Eq.~(\ref{eq:apeqmodi}) can be represented by the infinite series
\begin{eqnarray}
&& D(x,s) = \sum_{n=0}^{\infty} \prod_{i=1}^{n} \left\{
\int_{m_e^2}^{s_{i-1}} {{{\rm d} s_i} \over {s_i}} \Pi (s_{i-1}, s_i) \right. \nonumber \\
&& \! \! \! \left. \times {{\alpha} \over {2 \pi}}
\int_{x/(z_1 \cdots z_{i-1})}^{x_+} {{{\rm d} z_i} \over {z_i}} P(z_i) \right\}
\Pi (s_n, m_e^2) D \left( {{x} \over {z_1 \cdots z_n}} , m_e^2 \right) . \nonumber \\
\label{eq:dmc}
\end{eqnarray}
The particular form of Eq.~(\ref{eq:dmc}) allows to exploit a MC method for building the
solution iteratively. The steps of the algorithm are as follows:
\begin{itemize}
\item[1 --]{set $Q^2 = m_e^2$, and fix $x=1$ according to the boundary condition $D(x,m_e^2) =
\delta(1-x)$; }
\item[2 --]{generate a random number $\xi$ in the interval $[0,1]$; }
\item[3 --]{if $\xi < \Pi (s,Q^2)$} stop the evolution; otherwise
\item[4 --]{compute $Q'^2$ as solution of the equation $\xi = \Pi (Q'^2, Q^2)$};
\item[5 --]{generate a random number $z$ according to the probability density $P(z)$ in the
interval $[0,x_+]$; }
\item[6 --]{substitute $x \to xz$ and $Q^2 \to Q'^2$}; go to 2.
\end{itemize}
The $x$ distribution of the electron SF as obtained by means of the PS algorithm and
a numerical solution (based on the inverse Mellin transform method) of the QED evolution
equation is shown in Fig.~\ref{figmp:1}. Perfect agreement is seen. Once $D(x,s)$ has been reconstructed
by the algorithm, the master formula of
Eq.~(\ref{eq:masterd}) can be used for the calculation of LL corrections
to the cross section of interest. This cross section must be independent of the soft-hard
photon separator $\epsilon$ in the limit of small values for $\epsilon$. This can be clearly seen
in Fig.~\ref{figmp:2}, where the QED corrected Bhabha cross section as a function of the fictitious
parameter $\varepsilon$ is shown for DA$\mathrm{\Phi}$NE energies with the
cuts of Eq.~(\ref{eq:cutsmf}), but for an angular acceptance
$\theta_{\pm}$ of $55^\circ \div 125^\circ$. The cross section reaches a plateau for $\epsilon$ smaller
than $10^{-4}$.
\begin{figure}
\begin{center}
\resizebox{0.45\textwidth}{!}{%
\includegraphics{fig1_ps.eps}
}
\caption{Comparison for the $x$ distribution of the electron SF %at %$\sqrt{s} = 190$ GeV,
as obtained by means of a numerical solution of the QED evolution equation
(solid line) and the PS algorithm (histogram). From \cite{CarloniCalame:2000pz}.}
%\label{fig3aa-lumi}
\label{figmp:1}
\end{center}
\end{figure}
\begin{figure}
\begin{center}
\resizebox{0.5\textwidth}{!}{%
\includegraphics{eps-ind.eps}
}
\caption{QED corrected Bhabha cross section at DA$\mathrm{\Phi}$NE as a function of the infrared regulator
$\varepsilon$ of the PS approach, according to
the setup of Eq.~(\ref{eq:cutsmf}). The error bars correspond to $1\sigma$ MC errors. From
\cite{Balossini:2006wc}.}
\label{figmp:2}
\end{center}
\end{figure}
The main advantage of the PS algorithm with respect to the %collinear
analytical solutions
of the electron evolution is the possibility of going beyond the
strictly collinear approximation and generating transverse momentum $p_\perp$
of electrons and photons at each branching. In fact, the kinematics
of the branching process $e(p) \to e'(p') + \gamma(q)$ can be written as
\begin{eqnarray}
&& p=(E, \vec{0}, p_z)\,, \nonumber\\
&& p'=(zE, \vec{p}_\perp, p'_z)\,, \nonumber\\
&& q=((1-z)E, - \vec{p}_\perp, q_z)\,.
\label{eq:kinealt}
\end{eqnarray}
Once the variables $p^2$, ${p'}^2$ and $z$ are generated by the PS algorithm,
the on-shell condition $q^2=0$, together with the longitudinal momentum
con\-ser\-va\-tion, al\-lows to obtain an expression for the $p_{\perp}$
variable:
%\eqnarray
%&&p_z=E-\frac{k^2}{2E} \nonumber \\
%&&p'_z=zE-\frac{(1-z)k^2+{k'}^2}{2E} \nonumber \\
%&&q_z=(1-z)E-\frac{zk^2-{k'}^2}{2E} \nonumber \\
\begin{equation}
p^{2}_{\perp}=(1-z)(zp^2-p^{\prime 2}), \ \
\label{eq:kinep}
\end{equation}
%\endeqnarray
valid at first order in $p^2 / E^2 \ll 1$, $p^2_\perp / E^2 \ll 1$.
%%%TT
However, due to the approximations inherent to Eq.~(\ref{eq:kinep}), this PS approach can lead to an incorrect behaviour of the reconstruction of the exclusive photon kinematics.
%%%%%%%%%%%%%%%%%%%%%
First of all, since
within the PS algorithm the generation of $p^{\prime 2}$ and $z$ are
independent, it can happen that in some branchings the $p^2_{\perp}$ as given
by Eq.~(\ref{eq:kinep}) is negative. In order to avoid this problem, the
introduction of any kinematical cut on the $p^2$ or $z$ generation (or the
regeneration of the whole event) would prevent the correct reconstruction of
the SF $x$ distribution, which is important for a precise cross section
calculation. Furthermore, in the PS scheme, each fermion produces its photon
cascade independently of the other ones, missing the effects due to the
interference of radiation coming from different charged particles. As far as
inclusive cross sections
%
%%%%%%%%%%% REVISED %%%%%%%%%%%%%%
(i.e. cross sections with no cuts imposed on the generated photons)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
are concerned, these effects are largely integrated out.
%and we do not expect them to be large,
However, as shown in \cite{CarloniCalame:2001ny}, they become important when more exclusive variables distributions are considered.
The first problem can be overcome by
%
%%%%%%%% REVISED %%%%%%%%%%%%%
%leaving away from a pure PS picture.
choosing the generated $p_\perp$ of the photons different from Eq.~(\ref{eq:kinep}).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
For example, one can choose to extract the photon
$\cos\vartheta_\gamma$ according to the universal leading poles $1/p\cdot k$
present in the matrix element for photon emission. %
%%%%%%%% REVISED
%as done for example in \cite{psohl}.
%%%%%%%%%%%%%%%%%%%
%
Namely, one can generate $\cos\vartheta_\gamma$ as
\begin{equation}
\cos\vartheta_\gamma \propto \frac{1}{1-\beta\cos\vartheta_\gamma} \ \ ,
\label{cthg_LL}
\end{equation}
where $\beta$ is the speed of the emitting particle. In this way, photon
energy and angle are generated independently, different from Eq.~(\ref{eq:kinep}).
The nice feature of this prescription is that
$p^2_{\perp}=E^2_\gamma\sin^2\vartheta_\gamma$ is always well defined, and the
$x$ distribution reproduces exactly the SF, because
%%%TT
no further kinematical cuts have to be imposed to avoid unphysical events. At this stage, the
PS is used only to generate the
energies and multiplicity of the photons.
%, we simply
%do not ``believe'' in the reliability of eq. (\ref{eq:kinep}) for $p_\perp$
%reconstruction.
The problem of including the radiation interference is still unsolved, because
the variables of photons emitted by a fermion are still uncorrelated with those of the
other charged particles. The issue of including photon interference can be successfully worked out
looking at the YFS formula \cite{Yennie:1961ad}:
\begin{equation}
{\rm d}\sigma_n\approx {\rm d}\sigma_0\frac{e^{2n}}{n!}\prod_{l=1}^n\frac{{\rm d}^{\it 3}
\mathbf{k}_l}{(2\pi)^32k^0_l}\sum_{i,j=1}^N\eta_i\eta_j\frac{-p_i\cdot p_j}
{(p_i\cdot k_l)(p_j\cdot k_l)} \ \ .
\label{yfse}
\end{equation}
It %is a very general formula which
gives the differential cross section
${\rm d}\sigma_n$ for the emission of $n$ photons, whose momenta are
$k_1,\cdots,k_n$, from a
kernel process described by ${\rm d}\sigma_0$ and involving $N$ fermions, whose
momenta are $p_1,\cdots,p_N$. In Eq.~(\ref{yfse}) $\eta_i$ is a charge
factor, which is $+1$ for incoming $e^-$ or outgoing $e^+$ and $-1$ for
incoming $e^+$ or outgoing $e^-$. Note that Eq.~(\ref{yfse}) is valid in the
soft limit ($k_i\to 0$). The important point is that it also accounts for
coherence effects. From the YFS formula it is straightforward to read out the
angular spectrum of the $l^{th}$ photon:
\begin{equation}
\cos\vartheta_l\propto-\sum_{i,j=1}^N\eta_i\eta_j
\frac{1-\beta_i\beta_j\cos\vartheta_{ij}}{(1-\beta_i\cos\vartheta_{il})
(1-\beta_j\cos\vartheta_{jl})} \ \ .
\label{cthYFS}
\end{equation}
It is worth noticing that in the LL prescription the same quantity can be written as
\begin{equation}
\cos\vartheta_l\propto\sum_{i=1}^N\frac{1}{1-\beta_i\cos\vartheta_{il}},
\label{cthLL}
\end{equation}
whose terms are of course contained in Eq.~(\ref{cthYFS}).
In order to consider also coherence effects in the angular distribution of the
photons, one can generate $\cos\vartheta_\gamma$ according to Eq.~(\ref{cthYFS}), rather
than to Eq.~(\ref{cthLL}). This recipe \cite{CarloniCalame:2001ny}
is adopted in BabaYaga v3.5 and BabaYaga@NLO. \\
\noindent {\em{Yennie-Frautschi-Suura exponentiation \label{YFS}}} \\
\vspace*{-2mm}
%Alternative method: Yennie-Frautschi-Suura exponentiation (as in BHWIDE).
The YFS exponentiation procedure,
implemented in the code BHWIDE, is %essentially
a technique for summing up all the infrared (IR) singularities present in any %given
process
accompanied by photonic radiation~\cite{Yennie:1961ad}. It is inherently exclusive, i.e. all the
summations of the IR singular
contributions are done before any phase-space integration over the virtual or real photon
four-momenta are performed. %~\cite{jw90}.
The method was mainly developed by S. Jadach, B.F.L. Ward and collaborators to
realise precision MC tools.
In the following, the general ideas underlying the procedure are summarised.
%(see for instance~\cite{wascern} for a more detailed description).
Let us consider the scattering process $e^+ (p_1) e^- (p_2) \to f_1 (q_1) \cdots f_n (q_n)$,
where $f_1 (q_1) \cdots f_n (q_n)$ represents a given arbitrary final state, and let
${\cal M}_0$ be its tree-level matrix element.
By using standard Feynman-diagram techniques, it is possible
to show that the same process, when accompanied by $l$
additional real photons radiated by the
IS particles, and under the assumption that the $l$ additional photons are soft,
i.e. their energy is much smaller that any energy scale involved in the process, can be
described by the factorised matrix element built up by the LO one, ${\cal M}_0$,
times the product of $l$ eikonal currents, namely
\begin{equation}
{\cal M} \simeq {\cal M}_0 \prod_{i=1}^{l} \left[ e \left(
{{\varepsilon_i (k_i) \cdot p_2} \over {k_i \cdot p_2}} -
{{\varepsilon_i (k_i) \cdot p_1} \over {k_i \cdot p_1}}
\right) \right] ,
\label{eq:lrphot}
\end{equation}
where $e$ is the electron charge, $k_i$ are the momenta of the photons and
$\varepsilon_i (k_i)$ their polarisation vectors.
%(see Fig.~\ref{fig:realg} for a representation of the
%Feynman diagrams relative to the initial-state ${\cal O} (\alpha)$ bremsstrahlung correction).
Taking the square of the matrix element
in Eq.~(\ref{eq:lrphot}) and multiplying by the proper flux factor and the Lorentz-invariant
phase space volume,
%phase space volume, neglecting the four-momenta of the radiated photons in the
%overall energy-momentum conservation, summing over the final state photons polarisations
%and combining properly all the factors,
the cross section for the process
$e^+ (p_1) e^- (p_2) \to f_1 (q_1) \cdots f_n (q_n) + \, l \, {\rm real} \, {\rm photons}$ can be written as
\begin{eqnarray}
{\rm d} \sigma^{(l)}_{r} = && {\rm d} \sigma_0 {{1} \over {l!}} \prod_{i=1}^{l} \left[
k_i {\rm d} k_i {\rm d} \cos \vartheta_i {\rm d} \varphi_i {{1} \over {2 (2 \pi)^3}} \right. \nonumber \\
&& \left. \times\,\sum_{\varepsilon_i} e^2 \left(
{{\varepsilon_i (k_i) \cdot p_2} \over {k_i \cdot p_2}} -
{{\varepsilon_i (k_i) \cdot p_1} \over {k_i \cdot p_1}}
\right)^2 \right] .
\label{eq:lreal}
\end{eqnarray}
By summing over the number of final-state photons, one obtains the cross section for the
original process accompanied by an arbitrary number of real photons, namely
\begin{eqnarray}
{\rm d} \sigma^{(\infty)}_{r} =&& \sum_{l=0}^{\infty} {\rm d} \sigma^{(l)}_{r} \nonumber \\
=&& {\rm d} \sigma_0 \exp \left[
k {\rm d} k {\rm d} \cos \vartheta {\rm d} \varphi {{1} \over {2 (2 \pi)^3}} \right. \nonumber\\
&& \left. \times\,\sum_{\varepsilon} e^2 \left(
{{\varepsilon (k) \cdot p_2} \over {k \cdot p_2}} -
{{\varepsilon (k) \cdot p_1} \over {k \cdot p_1}}
\right)^2 \right] .
\label{eq:expr}
\end{eqnarray}
Equation~(\ref{eq:expr}), being limited to real radiation only, is IR divergent once the phase
space integrations are performed down to zero photon energy. This problem, as is well known,
finds its solution in the matching between real and virtual photonic radiation.
Equation~(\ref{eq:expr}) already shows the key feature of exclusive exponentiation, i.e. summing
up all the perturbative contributions before performing any phase space integration.
In order to get meaningful radiative corrections it is necessary to consider, besides IS real photon corrections, also IS virtual photon
corrections, i.e. the corrections due to additional internal photon lines connecting the
IS electron and positron. %(see Fig.~\ref{fig:virtualg}).
For a vertex-type
amplitude, the result can be written as
\begin{eqnarray}
{\cal M}_{V_1} = - i {{e^2} \over {(2 \pi)^4}} \int {\rm d}^4 k && \!\! \!\! \!\! {{1} \over {k^2 + i \varepsilon}}
{\bar v} (p_1) \gamma^\mu
{{- ( {\rlap / {p_1}} + {\rlap / {k}} ) + m } \over {2 p_1 \cdot k + k^2 + i
\varepsilon}} \nonumber \\
&\times & \!\! \Gamma
{{ ( {\rlap / {p_2}} + {\rlap / {k}} ) + m } \over {2 p_2 \cdot k + k^2 + i \varepsilon}}
\gamma_\mu u(p_2) ,
\label{eq:virt}
\end{eqnarray}
where $\Gamma$ stands for the Dirac structure of the LO process, in such a
way that ${\cal M}_0 = {\bar v} (p_1) \Gamma u(p_2)$. The soft-photon part of the amplitude
can be extracted by taking $k^\mu \simeq 0 $ in all the numerators. In this approximation, the
amplitude of Eq.~(\ref{eq:virt}) becomes
\begin{eqnarray}
&& {\cal M}_{V_1} = {\cal M}_0 \times V , \nonumber \\
&& V = {{2 i \alpha} \over {(2 \pi)^3}} \int {\rm d}^4 k
{{4 p_1 \cdot p_2} \over {(2 p_1 \cdot k + k^2 + i \varepsilon)
(2 p_2 \cdot k + k^2 + i \varepsilon)}} \nonumber\\
&& \qquad \qquad \qquad \quad \times {{1} \over {k^2 + i \varepsilon}} .
\label{eq:onevirt}
\end{eqnarray}
%Some comments are in order here. First,
It can be seen that, as in the real case, the IR virtual correction
factorises off the LO matrix element so that it is universal, i.e. independent
of the details of the process under consideration, and divergent in the IR portion of the
phase space.
The correction given by $n$ soft virtual photons can be seen to factorise with an additional factor $1/n!$, namely
\begin{equation}
{\cal M}_{V_n} = {\cal M}_0 \times {{1} \over {n!}} V^n ,
\label{eq:nvirt}
\end{equation}
so that by summing over all the additional soft virtual photons one obtains
\begin{equation}
{\cal M}_V = {\cal M}_0 \times \exp [V] .
\label{eq:virtexp}
\end{equation}
As already noticed both the real and virtual factors are IR divergent. In order to obtain
meaningful expressions one has to adopt some regularisation procedure. One
possibility is to give the photon a (small) mass $\lambda$ and to modify
Eqs.~(\ref{eq:lreal}) and (\ref{eq:onevirt}) accordingly. Once all the expressions are
properly regularised, one can write down a YFS master formula that takes into account
real and virtual photonic corrections to the LO process. In virtue of the
factorisation properties discussed above, the master formula can be obtained from
Eq.~(\ref{eq:expr}) with the substitution ${\rm d} \sigma_0 \to {\rm d} \sigma_0 \vert \exp (V) \vert^2$, i.e.
\begin{eqnarray}
{\rm d} \sigma && = {\rm d} \sigma_0 \vert \exp (V) \vert^2 \exp \left[
k {\rm d} k {\rm d} \cos \vartheta {\rm d} \varphi {{1} \over {2 (2 \pi)^3}} \right. \nonumber\\
&& \left. \times\,\sum_{\varepsilon} e^2 \left(
{{\varepsilon (k) \cdot p_2} \over {k \cdot p_2}} -
{{\varepsilon (k) \cdot p_1} \over {k \cdot p_1}}
\right)^2 \right] .
\label{eq:yfsmaster}
\end{eqnarray}
As a last step it is possible to analytically perform the IR cancellation between virtual
and very soft real photons. Actually, since very soft real photons do not affect the
kinematics of the process, the real photon exponent can be split into a contribution coming
from photons with energy less than a cutoff $k_{min}$ plus a contribution from photons
with energy above it. The first contribution can be integrated over all its phase
space and can then be combined with the virtual exponent. After this step it is possible to remove the regularising photon mass by taking the limit
$\lambda \to 0$, so that Eq.~(\ref{eq:yfsmaster}) becomes
\begin{eqnarray}
{\rm d} \sigma && = {\rm d} \sigma_0 \exp (Y) \exp \left[
k {\rm d} k {\rm d} \Theta (k - k_{min})\cos \vartheta {\rm d} \varphi {{1} \over {2 (2 \pi)^3}} \right. \nonumber\\
&& \left. \times\,\sum_{\varepsilon} e^2 \left(
{{\varepsilon (k) \cdot p_2} \over {k \cdot p_2}} -
{{\varepsilon (k) \cdot p_1} \over {k \cdot p_1}}
\right)^2 \right] ,
\label{eq:yfsfinal}
\end{eqnarray}
where $Y$ is given by
\begin{eqnarray}
Y && = 2 V + \int k {\rm d} k {\rm d} \Theta ( k_{min} - k )
\cos \vartheta {\rm d} \varphi {{1} \over {2 (2 \pi)^3}} \nonumber\\
&& \times\,\sum_{\varepsilon} e^2 \left(
{{\varepsilon (k) \cdot p_2} \over {k \cdot p_2}} -
{{\varepsilon (k) \cdot p_1} \over {k \cdot p_1}}
\right)^2 .
\end{eqnarray}
The explicit form of $Y$ can be derived by performing all the details of the calculation, and
reads
\begin{eqnarray}
Y &=& \beta \ln {{k_{min}} \over {E}} + \delta_{YFS} , \nonumber \\
\delta_{YFS} &=& {{1} \over {4}} \beta + {{\alpha} \over {\pi}} \left(
{{\pi^2} \over {3}} - {{1} \over {2}} \right) . %\nonumber \\
%&&\beta = {{2 \alpha} \over {\pi}} \left[ \ln \left( {{s} \over {m_e^2}} \right) - 1 \right] .
\end{eqnarray}
\subsubsection{Matching NLO and higher-order corrections
\label{NLO-HO}}
As will be shown numerically in Section \ref{NUMERICS}, NLO corrections must be combined with multiple
photon emission effects to achieve a theoretical accuracy at the per mill level.
%%%TT
This combination, technically known as {\em matching}, is a fundamental ingredient of
the most precise generators used for luminosity monitoring, i.e. BabaYaga@NLO,
BHWIDE and MCGPJ. Although the matching is implemented according to
different theoretical details, some general aspects are common to all the recipes
and must be emphasised:
\begin{enumerate}
\item It is possible to match NLO and HO corrections consistently, avoiding double
counting of LL contributions at order $\alpha$ and preserving the advantages of
resummation of soft and collinear effects beyond $O (\alpha)$.
\item The convolution of NLO corrections with HO terms allows
to include the dominant part of NNLO corrections, given by
infrared-enhanced $\alpha^2 L$ sub-leading contributions. This was argued
and demonstrated analytically and numerically in \cite{Montagna:1996gw} through comparison with the available
$O (\alpha^2)$ corrections to $s$-channel processes and $t$-channel
Bhabha scattering. Such an aspect of the matching procedure is crucial to
settle the theoretical accuracy of the generators by means of explicit comparisons
with the exact NNLO perturbative corrections discussed in Section \ref{NNLO}, and
will be addressed in Section \ref{TH}.
\item BabaYaga@NLO and BHWIDE implement a fully factorised
matching recipe, while MCGPJ includes some terms in additive form,
as will be visible in
the formulae reported below.
%This can give rise to some (minor) differences, when performing tuned comparisons
%between the programs predictions.
\end{enumerate}
\vspace{-0.2cm}
In the following we summarise the basic features of the matching
procedure as implemented in the codes MCGPJ,
BabaYaga@NLO and BHWIDE.
The matching approach realised in the MC event generator MCGPJ was
developed in~\cite{Arbuzov:2005pt}. In particular, Bha\-bha scattering with complete
$O (\alpha)$ and HO LL photonic corrections can written as
\begin{eqnarray}
&& \frac{\dd\sigma^{e^+e^-\to e^+e^-(\gamma)}}{\dd\Omega_-}
= \int\limits_{\bar{z}_1}^{1}\dd z_1\;\int\limits_{\bar{z}_2}^{1}\dd z_2\;
D_{ee}^{NS,\gamma}(z_1) D_{ee}^{NS,\gamma}(z_2)
\nonumber \\ && \quad \times
\frac{\dd\hat{\sigma}_{0}^{\mathrm{Bhabha}}(z_1,z_2)}{\dd\Omega_-}
\left(1+\frac{\alpha}{\pi}K_{SV}\right)\Theta({\mathrm{cuts}})
\nonumber \\ \nonumber && \quad
\times \int\limits^{Y_1}_{y_{\mathrm{th}}}\frac{\dd y_1}{Y_1}\;
\int\limits^{Y_2}_{y_{\mathrm{th}}}\frac{\dd y_2}{Y_2}\;
D_{ee}^{NS,\gamma}(\frac{y_1}{Y_1})D_{ee}^{NS,\gamma}(\frac{y_2}{Y_2})
\\ \nonumber && \quad
+ \frac{\alpha}{\pi}\int\limits_{\Delta}^{1}\frac{\dd x}{x}
\Biggl\{\biggl[\left(1-x+\frac{x^2}{2}\right)\ln\frac{\theta_0^2(1-x)^2}{4}
+ \frac{x^2}{2}\biggr]
\\ \nonumber && \quad \times
2 \frac{\dd\sigma_0^{\mathrm{Bhabha}}}{\dd \Omega_-}
+ \biggl[\left(1-x+\frac{x^2}{2}\right)\ln\frac{\theta_0^2}{4}
+ \frac{x^2}{2}\biggr]
\\ \nonumber && \quad \times
\Biggl[ \frac{\dd\hat{\sigma}_0^{\mathrm{Bhabha}}(1-x,1)}{\dd \Omega_-}
+ \frac{\dd\hat{\sigma}_0^{\mathrm{Bhabha}}(1,1-x)}{\dd \Omega_-}
\Biggr]\Biggr\}\Theta({\mathrm{cuts}})
\\ && \quad
- \frac{\alpha^2}{4s}
\left(\frac{3+c^2}{1-c}\right)^2
\frac{8\alpha}{\pi}\ln(\mbox{ctg}\frac{\theta}{2})
\ln\frac{\Delta\eps}{\eps}
\nonumber \\ && \quad
+ \frac{\alpha^3}{2\pi^2s}\!\!\!\!\!
\int\limits_{\stackrel{k^0>\Delta\eps}{\theta_i >\theta_0}}
\!\! \frac{WT}{4}\Theta({\mathrm{cuts}}) \frac{\dd \Gamma_{e\bar{e}\gamma}}{\dd \Omega_-}\, .
\label{bhabha_nlo}
\end{eqnarray}
Here the step functions $\Theta({\mathrm{cuts}})$ stand for the particular
cuts applied.
The auxiliary parameter $\theta_0$ defines cones around the directions of
the motion of the charged particles
in which the emission of hard photons is approximated by the factorised form
by convolution of collinear radiation factors~\cite{Arbuzov:2007wu}
with the Born cross section.
The dependence on the parameters $\Delta$ and $\theta_0$ cancels out in the sum with the last term
of Eq.~(\ref{bhabha_nlo}), where the photon energy and emission angles with respect to all charged
particles are limited from below $(k^0 > \Delta\eps, \theta_i > \theta_0)$.
Taking into account vacuum polarisation, the Born level Bhabha cross section with reduced
energies of the incoming electron and positron can be cast in the form
\begin{eqnarray}
&& \frac{\dd\hat{\sigma}_{0}^{\mathrm{Bhabha}}(z_1,z_2)}{\dd\Omega_-} =
\frac{4\alpha^2}{sa^2}\biggl\{\frac{1}{|1-\Pi(\hat{t})|^2}\;
\frac{a^2+z_2^2(1+c)^2}{2z_1^2(1-c)^2}
\nonumber \\ && \quad
+ \frac{1}{|1-\Pi(\hat{s})|^2}\;
\frac{z_1^2(1-c)^2+z_2^2(1+c)^2}{2a^2}
\nonumber \\ && \quad
- \mbox{Re}\;\frac{1}{(1-\Pi(\hat{t}))(1-\Pi(\hat{s}))^{*}}\;
\frac{z_2^2(1+c)^2}{az_1(1-c)} \biggr\}\dd\Omega_-\, ,
\nonumber \\ &&
\hat{s} = z_1z_2s, \quad \hat{t} = -\frac{sz_1^2z_2(1-c)}{z_1+z_2-(z_1-z_2)c}\,,
\end{eqnarray}
where $\Pi(Q^2)$ is the photon self-energy correction. Note that in the cross section
above the cosine of the scattering angle, $c$, is
given for the original c.m. reference frame of the colliding beams.
For the two-photon production channel, a similar representation is used in MCGPJ:
\begin{eqnarray}
&& \dd\sigma^{e^+e^-\to\gamma\gamma(\gamma)}
= \int\limits_{\bar{z}_1}^{1}\dd z_1 D_{ee}^{NS,\gamma}(z_1)
\int\limits_{\bar{z}_2}^{1}\dd z_2 D_{ee}^{NS,\gamma}(z_2)
\nonumber \\ && \quad \times
\dd\hat{\sigma}_0^{\gamma\gamma}(z_1,z_2)
\left(1+\frac{\alpha}{\pi}K_{SV}^{\gamma\gamma}\right)
+\ \frac{\alpha}{\pi}\int\limits_{\Delta}^{1}\frac{\dd x}{x}
\nonumber \\ && \quad \times
\left[\left(1-x+\frac{x^2}{2}\right)\ln\frac{\theta^2_0}{4}
+\frac{x^2}{2}\right]
\biggl[\dd\hat{\sigma}_0(1-x,1)
\nonumber \\ && \quad
+ \dd\hat{\sigma}_0(1,1-x)\biggr]
+\ \frac{1}{3}\,\frac{4\alpha^3}{\pi^2s^2} \!\!\!\!\int\limits_{\stackrel{z_i\geq\Delta}
{\pi-\theta_0\geq\theta_i\geq\theta_0}}\!\!
\dd \Gamma_{3\gamma}
\nonumber \\ && \quad \times
\biggl[\frac{z_3^2(1+c_3^2)}{z_1^2z_2^2(1-c_1^2)(1-c_2^2)}
+ \mbox{two cyclic permutations}\biggr],
\nonumber \\
&& z_i=\frac{q_i^0}{\eps},\quad c_i=\cos\theta_i,\quad
\theta_i=\widehat{\vec{p}_-\vec{q}}_i \, ,
\label{eq:agg}
\end{eqnarray}
where
the cross section with reduced energies has the form
\begin{eqnarray} \nonumber
\frac{\dd\hat{\sigma}_0^{\gamma\gamma}(z_1,z_2)}{\dd\Omega_1}
= \frac{2\alpha^2}{s}\frac{z_1^2(1-c_1)^2+z_2^2(1+c_1)^2}
{(1-c_1^2)(z_1+z_2+(z_2-z_1)c_1)^2} \, ,
\end{eqnarray}
and the factor $1/3$ in the last term of Eq.~(\ref{eq:agg}) takes into
account the identity of the final-state photons.
The sum of the last two terms does not
depend on $\Delta$ and $\theta_0$.
%Note that the annihilation process
%is a pure QED one, hadronic contributions as well as weak interaction
%effects are far beyond the required accuracy.
Concerning BabaYaga@NLO, the matching starts from the observation that
Eq.~(\ref{eq:masterd}) for the QED corrected all-order cross section
can be rewritten in terms of the PS ingredients as
\begin{equation}
{\rm d}\sigma^{\infty}_{LL}=%\frac1{\mbox{flux}}
{\Pi}(Q^2,\varepsilon)~
\sum_{n=0}^\infty \frac{1}{n!}~|{\ourcal M}_{n,LL}|^2~{\rm d}\Phi_n\,.
\label{generalLL}
\end{equation}
By construction, the expansion of Eq.~(\ref{generalLL}) at \oa\ does not coincide with the exact \oa\ result. In fact
\bea
{\rm d}\sigma^{\alpha}_{LL}&=&
\left[
%%%TT: \ln not \log
1-\frac{\alpha}{2\pi}~I_+~\ln\frac{Q^2}{m^2}
\right] |{\ourcal M}_0|^2 {\rm d}\Phi_0+
|{\ourcal M}_{1,LL} |^2 {\rm d}\Phi_1 \nonumber\\
&\equiv&
\left[
1+C_{\alpha,LL}
\right] |{\ourcal M}_0|^2 {\rm d}\Phi_0
+
|{\ourcal M}_{1,LL} |^2 {\rm d}\Phi_1\,,
\label{LL1}
\eea
where $I_+ \equiv \int_0^{1-\epsilon} P(z) {\rm d}z$, whereas the exact NLO cross section can always be cast in the form
\be
{\rm d}\sigma^{\alpha}
=
\left[
1+C_{\alpha}
\right] |{\ourcal M}_0|^2 {\rm d}\Phi_0
+
|{\ourcal M}_{1} |^2 {\rm d}\Phi_1\,.
\label{exact1}
\ee
The coefficients $C_{\alpha}$ contain the complete \oa\ virtual
and soft-bremsstrahlung corrections in units of
the squared Born amplitude,
and $|{\ourcal M}_1|^2$ is the exact squared matrix element with the
emission of one hard photon.
We remark that $C_{\alpha,LL}$ has the same logarithmic structure as
$C_\alpha$ and that $|{\ourcal M}_{1,LL}|^2$ has the same singular
behaviour as $|{\ourcal M}_1|^2$.
In order to match the LL and NLO calculations,
the following correction factors, which are by construction
infrared safe and free of collinear logarithms, are introduced:
\be
F_{SV}~=~
1+\left(C_\alpha-C_{\alpha,LL}\right),~~~
F_H~=~
1+\frac{|{\ourcal M}_1|^2-|{\ourcal M}_{1,LL}|^2}{|{\ourcal M}_{1,LL}|^2}\,.
\label{FSVH}
\ee
With them the exact \oa\ cross section can be expressed, up
to terms of ${\ourcal O}(\alpha^2)$,
in terms of its LL approximation as
\be
{\rm d}\sigma^\alpha~=~
F_{SV} (1+C_{\alpha,LL} ) |{\ourcal M}_0|^2 {\rm d}\Phi_0
~+~
F_H |{\ourcal M}_{1,LL}|^2 {\rm d}\Phi_1 .
\label{matchedalpha}
\ee
Driven by Eq.~\myref{matchedalpha}, Eq.~\myref{generalLL} can be improved
by writing the resummed matched cross section as
\begin{eqnarray}
{\rm d}\sigma^{\infty}_{\rm matched} = &&
F_{SV}~\Pi(Q^2,\varepsilon) \nonumber\\
&& \times\,\sum_{n=0}^\infty \frac{1}{n!}~
\left( \prod_{i=0}^n F_{H,i}\right)~
|{\ourcal M}_{n,LL}|^2~
{\rm d}\Phi_n .
\label{matchedinfty}
\end{eqnarray}
The correction factors $F_{H,i}$ follow from the definition \myref{FSVH}
for each photon emission.
The \oa\ expansion of Eq.~\myref{matchedinfty} now coincides with
the exact NLO cross section of Eq.~\myref{exact1}, and
all HO LL contributions
are the same as in Eq.~\myref{generalLL}. This formulation is implemented in BabaYaga@NLO for both Bhabha scattering and
photon pair production, using, of course, the appropriate SV and hard
bremsstrah\-lung formulae. This matching formulation has also been
applied to the study of Drell-Yan-like processes,
by combining the complete $O (\alpha)$ electroweak corrections
with QED shower evolution in the
generator HORACE
\cite{CarloniCalame:2003ux,CarloniCalame:2005vc,CarloniCalame:2006zq,CarloniCalame:2007cd}.
As far as BHWIDE is concerned, this MC event generator realises
the process
\begin{equation}
e^+(p_1) + e^-(q_1)\; \longrightarrow\; e^+(p_2) + e^-(q_2) \;
+ \gamma_1(k_1) + \ldots + \gamma_n(k_n)
\label{process}
\end{equation}
via the YFS exponentiated cross section formula
\begin{eqnarray}
{\rm d}\sigma= && e^{2\alpha\mathrm{Re}B+2\alpha
\tilde B}\,\sum_{n=0}^\infty{1\over n!}\int\prod_{j=1}^n{{\rm d}^3k_j\over k_j^0
}\int{{\rm d}^4y\over(2\pi)^4} \nonumber\\
&& \times\,e^{iy(p_1+q_1-p_2-q_2-\sum_jk_j)+D}
%\nonumber\cr
\bar\beta_n(k_1,\dots,k_n)\,{{\rm d}^3p_2{\rm d}^3q_2\over p_2^0q_2^0}, \nonumber\\
\label{eqone}
\end{eqnarray}
where the
real infrared function $\tilde B$ and the virtual infrared function $B$ are
given in~\cite{Jadach:1995nk}.
%refs.~\cite{yfs:1961,bhlumi1:1989,sjbw:1988,bhlumi2:1992,bennie:1987}
%(for definiteness,
%we record them in the form we shall use presently, as this representation
%is not readily available in the current literature and is essential
%for practical applications of the type we pursue here),
%%%TT
Here we note the usual connections
\begin{eqnarray}
2\alpha\tilde B && = \int^{k\le K_{max}}{{\rm d}^3k\over k_0}\, \tilde S(k)\,,
\nonumber\\
D && =\int {\rm d}^3k\,{\tilde S(k)\over k^0}
\left(e^{-iy\cdot k}-\theta(K_{\rm max}-k)\right)
\label{eqtwo}
\end{eqnarray}
for the standard YFS infrared real emission factor
\begin{equation}
\tilde S(k)= {\alpha\over4\pi^2}\left[Q_fQ_{f'}
%Q_{\overset{\tiny(-)}{f'}}
\left({p_1\over p_1\cdot k}-{q_1
\over q_1\cdot k}\right)^2+\ldots\right]\,,
\label{eqthree}
\end{equation}
and where $Q_f$ is the electric charge of $f$ in units of the positron charge.
In Eq.~(\ref{eqthree}) the ``$\ldots$''
represent the remaining terms in $\tilde S(k)$, obtained from the given one
by respective of $Q_f$, $p_1$,
%$Q_{\overset{\tiny(-)}{f'}}$
$Q_{f'}$, $q_1$ with corresponding values for the
other pairs of the external charged legs according to the YFS
prescription of Ref.~\cite{Yennie:1961ad,Mahanthappa:1962ex} (wherein due attention is taken to obtain
the correct relative sign of each of the terms in $\tilde S(k)$ according to
this latter prescription). The explicit representation is given by
\begin{eqnarray}
&& 2\alpha\mathrm{Re} B(p_1,q_1,p_2,q_2) + 2\alpha\tilde{B}(p_1,q_1,p_2,q_2;k_m)
= \nonumber\\
&& R_1(p_1,q_1;k_m) + R_1(p_2,q_2;k_m)
%\nonumber\\
+ R_2(p_1,p_2;k_m) + \nonumber\\
&& R_2(q_1,q_2;k_m) - R_2(p_1,q_2;k_m) - R_2(q_1,p_2;k_m)\,,
\label{YFSff}
\end{eqnarray}
with
\begin{equation}
R_1(p,q;k_m) = R_2(p,q;k_m) + \left(\frac{\alpha}{\pi}\right)
\frac{\pi^2}{2}
\label{R1}
\end{equation}
and
\begin{eqnarray}
R_2(p,q;k_m) && = \frac{\alpha}{\pi} \left\{ \bigg(\ln\frac{2pq}{m_e^2} -1 \bigg)
\ln\frac{k_m^2}{p^0q^0} + \frac{1}{2}\ln\frac{2pq}{m_e^2} \right. \nonumber \\
&& \left. -\frac{1}{2}\ln^2\frac{p^0}{q^0}
-\frac{1}{4}\ln^2\frac{(\Delta+\delta)^2}{4p^0q^0}
-\frac{1}{4}\ln^2\frac{(\Delta-\delta)^2}{4p^0q^0} \right. \nonumber \\
&& \left. -\mathrm{Re}\, {\rm Li}_2\left(\frac{\Delta+\omega}{\Delta+\delta}\right)
-\mathrm{Re}\, {\rm Li}_2\left(\frac{\Delta+\omega}{\Delta-\delta}\right) \right. \nonumber \\
&& \left. -\mathrm{Re}\, {\rm Li}_2\left(\frac{\Delta-\omega}{\Delta+\delta}\right) %\right. \nonumber \\
-\mathrm{Re}\, {\rm Li}_2\left(\frac{\Delta-\omega}{\Delta-\delta}\right) \right. \nonumber\\
&& \left. +\frac{\pi^2}{3} -1 \right\},
\label{R2}
\end{eqnarray}
%\end{align}
where $\Delta=\sqrt{2pq+(p^0-q^0)^2}$, $\omega=p^0+q^0$, $\delta=p^0-q^0$,
and $k_m$ is a soft photon cut-off in the c.m. system
($E_{\gamma}^{\rm soft}