
\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}<k_m\ll E_{\rm beam}$).

%\par
The YFS hard photon residuals $\bar\beta_i$ in Eq.~(\ref{eqone}), $i=0,1$,
to $O (\alpha)$ are given exactly 
in Ref.~\cite{Jadach:1995nk} for BHWIDE. Therefore this 
event generator calculates the YFS exponentiated exact
$O (\alpha)$ cross section for $e^+e^-\rightarrow e^+ e^- + n(\gamma)$
with multiple initial, initial-final and final state radiation, 
using a corresponding MC realisation
of Eq.~(\ref{eqone}) in the wide angle regime. The library for $O (\alpha)$ electroweak corrections, relevant
for higher energies, is taken from \cite{Beenakker:1990mb,Bohm:1986fg}.\par
The result (\ref{eqone}) is an exact rearrangement of the loop expansion for the
respective cross section and is independent of the dummy parameter $K_{\rm max}$. To derive this, one may proceed as follows.
Let the amplitude for the emission of $n$ real photons in the
Bhabha process be 
\begin{equation}
{\cal M}^{(n)} 
 = \sum_{\ell}M^{(n)}_{\ell},
\label{subp1}
\end{equation}
where $M^{(n)}_{\ell}$ is the contribution to 
${\cal M}^{(n)}$ from Feynman diagrams with ${\ell}$ virtual loops.
The key result in the YFS theory of Ref.~\cite{Yennie:1961ad,Mahanthappa:1962ex} on virtual corrections
is that we may rewrite Eq.~(\ref{subp1}) as the exact representation
\begin{equation}
{\cal M}^{(n)} = e^{\alpha B}\sum_{j=0}^\infty{\sf m}^{(n)}_j,
\label{yfsrepv}
\end{equation}
where we have defined
\begin{equation}
\alpha B = \int {d^4k\over (k^2-\lambda^2+i\epsilon)}S(k) ,
\label{vbfn}
\end{equation}
with the virtual infrared emission factor given by
\begin{eqnarray}
S(k) && = \frac{-i\alpha}{8\pi^2}\sum_{i'<j}Z_{i'}\theta_{i'}Z_j\theta_j\left(\frac{(2\bar{p}_{i'}\theta_{i'}-k)_\mu}{k^2-2k\bar{p}_{i'}\theta_{i'}+i\epsilon} \right. \nonumber\\
&& \left. +\frac{(2\bar{p}_{j}\theta_{j}+k)_\mu}{k^2+2k\bar{p}_{j}\theta_{j}+i\epsilon}\right)^2.
\label{vbfn1}
\end{eqnarray}
Here, $\lambda$ is an infrared regulator mass, and following Refs.~\cite{Yennie:1961ad,Mahanthappa:1962ex} we identify the sign of the $j$-th external line charge here as $Z_j=Q_j$ and $\theta_j=+(-)$ for outgoing (incoming) 4-momentum $\bar{p}_j$, so that here $\bar{p}_1=p_1,\;\bar{p}_2=q_1,\;\bar{p}_3=p_2,\;\bar{p}_4=q_2$,$\;Z_1=+1,\; \theta_1=-,\; Z_2=-1,\; \theta_2=-,\; Z_3=+1,\;\theta_3=+,\; Z_4=-1,\;\theta_4=+$. The amplitudes $\{{\sf m}^{(n)}_j\}$ are free of all virtual infrared divergences. 

Using the result (\ref{yfsrepv}) for ${\cal M}^{(n)}$, we get the attendant differential cross section by the standard methods as  
\begin{eqnarray}
  {\rm d}\hat\sigma^n && = {e^{2\alpha \mathrm{Re} B}\over {n !}}\int\prod_{l=1}^n
{{\rm d}^3k_l\over (k_l^2+\lambda^2)^{1/2}} \nonumber\\
&& \times \bar\rho^{(n)}(p_1,q_1,p_2,q_2,k_1,\cdots,k_n)\,
{{\rm d}^3p_2{\rm d}^3q_2\over p^0_2 q^0_2} \nonumber\\
&& \times \delta^{(4)}\left(p_1+q_1-p_2-q_2-\sum_{i=1}^nk_i\right) ,
\label{diff1}
\end{eqnarray}
where we have defined
\begin{equation}
\bar\rho^{(n)}(p_1,q_1,p_2,q_2,k_1,\cdots,k_n)=
\sum_{\rm spin} \left|\sum_{j=0}^\infty{\sf m}^{(n)}_j\right|^2 ,
\label{diff2}
\end{equation}
in the incoming $e^+e^-$ c.m. system.
Here we have absorbed the remaining kinematical factors for
the initial state flux and spin averaging into the
normalisation of the amplitudes ${\cal M}^{(n)}$ for pedagogical reasons,
so that the $\bar\rho^{(n)}$ are averaged over initial spins
and summed over final spins. We then use the key result of 
Ref.~\cite{Yennie:1961ad,Mahanthappa:1962ex} on real corrections to write the exact result 
\begin{eqnarray}
%%\begin{split}
&&\bar\rho^{(n)}(p_1,q_1,p_2,q_2,k_1,\cdots,k_n)=\prod_{i=1}^n\tilde{S}(k_i)\bar\beta_0+\cdots  + \nonumber\\
&& \sum_{i=1}^n\tilde{S}(k_i)\bar\beta_{n-1}(k_1,\ldots,k_{i-1},k_{i+1},\ldots,k_n) \nonumber\\
&& + \bar\beta_n(k_1,\ldots,k_n), 
%%\end{split}
\label{realem}
\end{eqnarray}
where the hard photon 
residuals $\bar\beta_j$ are determined recursively~\cite{Yennie:1961ad,Mahanthappa:1962ex}
and are free of all virtual and all real infrared singularities to all orders in $\alpha$. Introducing the result (\ref{realem}) into Eq.~(\ref{diff1}) and summing
over the number of real photons $n$ leads directly to master formula (\ref{eqone}). We see that it allows for exact exclusive treatment of hard photonic effects on
an event-by-event basis.
% 
