\newcommand{\ta}[1]{#1\hspace{-.42em}/\hspace{-.07em}} % SLASH
\newcommand{\beq}{\begin{equation}}
\newcommand{\eeq}{\end{equation}}
%\newcommand{\bea}{\begin{eqnarray}}
%\newcommand{\eea}{\end{eqnarray}}
\newcommand{\non}{\nonumber}
\newcommand{\nn}{\nonumber}
\newcommand{\eq}[1]{eq. (\ref{#1})}         
\newcommand{\Eq}[1]{Eq. (\ref{#1})}
\newcommand{\e}{\varepsilon} 
%\newcommand{\be}{\beta} 

\def\Li{\hbox{Li}}
\def\DAF{DA\char8NE}


\subsection{History and evolution of radiative return in precision physics}
\label{radret:hist}

The idea to use \emph{Initial State Radiation} to measure
hadronic cross sections from the threshold of a reaction up to the
centre-of-mass (c.m.) energy of colliders with fixed energies $\sqrt{s}$,
to reveal reaction mechanisms and to search for new mesonic states
consists in exploiting the process $e^+e^-\rightarrow hadrons + n
\gamma$, thus reducing the c.m. energy of the colliding
electrons and positrons and consequently the mass squared
$M^2_{\rm had}= s - 2 \sqrt{s} \: \: E_{\gamma}$ of the hadronic
system in the final state by emission of one or more photons. The
method is particularly well suited for modern meson factories
like DA$\mathrm{\Phi}$NE (detector KLOE), running at the
$\mathrm{\phi}$-resonance, BEPC-II (detector BES-III),
commissioned in 2008 and running at the $J/\psi$ and ${\psi}(2S)$-resonances,
PEP-II (detector BaBar) and KEKB (detector
Belle) at the $\Upsilon(4S)$-resonance. Their high
luminosities compensate for the $\alpha / \pi$ suppression
of the photon emission. DA$\mathrm{\Phi}$NE, BEPC-II,
PEP-II and KEKB cover the regions in $M_{\rm had}$ up to
1.02, 3.8 (maximally 4.6) and 10.6 GeV,
respectively (for the latter actually restricted to 4--5 GeV
if hard photons are detected). A big advantage of ISR
is the low point-to-point systematic errors of the
hadronic energy spectra. 
This is because 
the luminosity, the energy of the
electrons and positrons and many other contributions to
the detection efficiencies are determined once for the whole
spectrum. As a consequence, the overall normalisation error is the
same for all energies of the hadronic system. The term
\emph{Radiative Return} alternately used for ISR refers to
the appearance of pronounced resonances (e.g. $\rho,\ \omega,\ \phi,\ 
J/\psi,\ Z$) with energies below the collider energy. Reviews and
updated results can be found in the Proceedings of the
International Workshops in Pisa (2003) \cite{Pisa}, Nara (2004)
\cite{tau04}, Novosibirsk (2006) \cite{Budker}, Pisa (2006)
\cite{tau06}, Frascati (2008) \cite{Frascati08}, and Novosibirsk
(2008) \cite{tau08}.

Calculations of ISR date back to the sixties to seventies of
the $20^{th}$ century. For example, photon emission for muon pair
production in electron-positron collisions has been calculated in
Ref. \cite{Baier:1965bg}, for the $2 \pi $-final state in Refs.\ 
\cite{Baier:1965jz,Pancheri:1969yx}; the resonances $\rho,\ \omega$ and $\phi$ have been
implemented in Ref. \cite{Pancheri:1969yx}, the excitations 
$\psi(3100)$ and $\psi^{\prime}(3700)$ in Ref. \cite{Greco:1975rm}, and
the possibility to determine the pion form factor was discussed in
Ref. \cite{Chen:1974wv}. The application of ISR to the new high
luminosity meson factories, originally aimed at the determination
of the hadronic contribution to the vacuum polarisation, more
specifically the pion form factor, has materialised in the late
nineties. Early calculations of ISR for the colliders
DA$\mathrm{\Phi}$NE, PEP-II and KEKB can be found in
\cite{Spagnolo:1998mt,Khoze:2000fs,Benayoun:1999hm,Arbuzov:1998te}. In Ref. \cite{Arbuzov:1997je}
calculations of radiative corrections for pion and kaon production
below energies of 2 GeV have been reported. An impressive example
of ISR is the \emph{Radiative Return} to the region of the
\emph{Z}-resonance at LEP-2 with collider energies around
200 GeV \cite{Abbiendi:2003dh,Abdallah:2005ph,Achard:2005nb,Schael:2006wu} (see Fig.~\ref{fig:Rsigma.eps}).



\vglue 1.0 cm
\begin{figure}[htb]
%\begin{figure}
%\begin{center}
%\includegraphics[width=8cm,height=7cm]{Rsigma.ps}
\includegraphics[width=8cm,height=7cm]{DELPHI.eps}
%\vspace{0.5cm}
\caption{
The reconstructed distribution of $e^+ e^-   \rightarrow q \bar q$ events as a function of the invariant mass of the quark-antiquark system. The data has been taken for a collider energy range of 182 - 209 GeV. The prominent peak around 90 GeV represents the Z-resonance, populated after emission of photons in the initial state \cite{Abdallah:2005ph}.}
% (With kind permission of The European Physical Journal (EPJ)).} % courtesy PDG Ref. \cite{PDG}}
\label{fig:Rsigma.eps}
%\end{center}
\end{figure}



ISR became a powerful tool for the analysis of experiments
at low and intermediate energies with the development of
EVA-PHOKHARA, a Monte Carlo %event 
generator which is user
friendly, flexible and easy to implement into the software of the
existing detectors
\cite{Binner:1999bt,Czyz:2000wh,Rodrigo:2001jr,Kuhn:2002xg,Rodrigo:2001kf,Czyz:2002np,Czyz:2003ue,Czyz:2004rj,Czyz:2004nq,Czyz:2005as,Czyz:2008kw,Czyz:2004ua,Czyz:2007wi,Czyz:2007ue,Czyz:2006xf,Czyz:2008zz,Grzelinska:2008eb}.

EVA and its successor PHOKHARA  allow to simulate the
process $e^+e^- \to hadrons + \gamma$ for a variety of exclusive
final states.
 As a starting point EVA was constructed
\cite{Binner:1999bt} to simulate leading order ISR and FSR
for the $\pi^+ \pi^-$ channel, and additional soft and collinear ISR
was included on the basis of structure functions taken from
\cite{Caffo:1994dm}. Subsequently EVA was extended to
include the four-pion state \cite{Czyz:2000wh}, albeit without
FSR. Neglecting FSR and radiative corrections, i.e. including one-photon emission from the initial state only, the
cross section for the radiative return can be cast into the
product of a radiator function $H(M^2_{\rm had},s)$ and the cross
section $\sigma(M^2_{\rm had})$ for the reaction $e^+e^- \to hadrons$:

\noindent
$s \:  {{\rm d}\sigma(e^{+} e^{-} \rightarrow hadrons \: \gamma)}/{{\rm d}
M_{\rm had}^2}=\sigma (M_{\rm had}^2) \: H(M_{\rm had}^2,s)$.


However, for a precise evaluation of $\sigma(M^2_{\rm had})$, the
leading logarithmic approximation inherent in EVA is
insufficient. Therefore, in the next step, the exact one-loop
correction to the ISR process was evaluated analytically,
first for large angle photon emission \cite{Rodrigo:2001jr}, then for
arbitrary angles, including collinear configurations \cite{Kuhn:2002xg}. This
was and is one of the key ingredients of the generator called
PHOKHARA \cite{Rodrigo:2001kf,Czyz:2002np}, which also includes
soft and hard real radiation, evaluated using exact matrix
elements formulated within the framework of helicity amplitudes
\cite{Rodrigo:2001kf}. FSR in NLO approximation was
addressed in \cite{Czyz:2003ue} and incorporated in
\cite{Czyz:2004rj,Czyz:2004nq}. The importance of the charge
asymmetry, a consequence of interference between ISR and
FSR amplitudes, for a test of the (model dependent)
description of FSR has been emphasised already in
Ref. \cite{Binner:1999bt} and was further studied in  \cite{Czyz:2004nq}.

%(Frage: gibt es hierzu auch
%veroeffentlichte experimentelle Resultate von KLOE?)

Subsequently the generator was extended to allow for the generation of
many more channels with mesons, like $K^+K^-$, $K^0\bar K^0$,
$\pi^+\pi^-\pi^0$, for an improved description of the $4\pi$ modes
\cite{Czyz:2005as,Czyz:2008kw} and for improvements in the description of
FSR for the $\mu^+\mu^-$ channel
\cite{Czyz:2004rj,Czyz:2004nq}. Also the nucleon channels $p \bar{p}$
and $n \bar{n}$ were implemented \cite{Czyz:2004ua}, and it was
demonstrated that the separation of electric and magnetic proton
form factors is feasible for a wide energy range. In fact, for the
case of $\Lambda \bar{\Lambda}$ and including the 
polarisation-sensitive weak decay of $\Lambda$ into the simulation, it was
shown that even the relative phase between the two independent
form factors could be disentangled \cite{Czyz:2007wi}.

Starting already with \cite{Melnikov:2000gs}, various improvements
were made to include the direct decay $\phi \to \pi^+ \pi^-
\gamma$ as a specific aspect of FSR into the generator, a
contribution of specific importance for data taken on top of the
$\mathrm{\phi}$ resonance. 


This was further pursued in the event generators
FEVA and FASTERD based on EVA-PHOKHARA. FEVA
includes the effects of the direct decay $\phi
\rightarrow \pi^{-} \pi^{+} \gamma$ and the decay via the
$\rho$-resonance $\phi \rightarrow \rho^{\pm} \pi^{\mp}
\rightarrow \pi^{-} \pi^{+} \gamma$
\cite{Dubinsky:2004xv,Pancheri:2006cp,Pancheri:2007xt}. The code
FASTERD takes into account \emph{Final State Radiation} in
the framework of both Resonance Perturbation Theory and sQED,
\emph{Initial State Radiation}, their interference and also the
direct decays $ e^{+} e^{-} \rightarrow \phi \rightarrow (f_{o};
f_{o}+\sigma) \gamma \rightarrow \pi^{+} \pi^{-} \gamma$, $e^{+}
e^{-} \rightarrow \phi \rightarrow \rho^{\pm} \pi^{\mp}
\rightarrow \pi^{+} \pi^{-} \gamma$ and $e^{+} e^{-} \rightarrow
\rho \rightarrow \omega \pi^{o} \rightarrow \pi^{o} \pi^{o}
\gamma$ \cite{Shekhovtsova:2009yn}, with the possibility to include
additional models.


EVA-PHOKHARA was applied for the first time to an
experiment to determine the cross section $e^+e^-\rightarrow
\pi^{+} \pi^{-}$ from the reaction threshold up to the maximum
energy of the collider with the detector KLOE at DA$\mathrm{\Phi}$NE \cite{Cataldi:1999aa,Denig:2001ra,Aloisio:2001xq,Denig:2002ps,Valeriani:2002yk,Venanzoni:2002jb,Muller:2004mb,Aloisio:2003dw,Denig:2003jn,Valeriani:2004mp,Kluge:2004mc,Denig:2005eb,Kluge:2005ac,Aloisio:2005tm,Denig:2006kj,Muller:2006bk,Leone:2006bm,Venanzoni:2007zz,Nguyen:2008rv,Muller:2007zzd,Ambrosino:2007vj,Aloisio:2004bu,:2008en,Kluge:2008fb,Venanzoni:2009aq} (Section
\ref{rr:kloe}). The motivation was
the determination of the $2 \pi$ final state contribution to the hadronic
vacuum polarisation.

The determination of the hadronic contribution to the vacuum
polarisation, which arises from the coupling of virtual photons to
quark-antiquark pairs, $ \gamma ^ \star \rightarrow q \bar {q}
\rightarrow \gamma ^\star $, is possible by measuring the cross
section of electron-positron annihilation into hadrons, $e^+e^-
\rightarrow \gamma^* \rightarrow q \bar {q} \rightarrow hadrons$, and
applying the optical theorem. It is of great importance for the
interpretation of the precision measurement of the anomalous
magnetic moment of the muon $a_\mu$ in Brook\-ha\-ven (E821) \cite{Bennett:2002jb,Bennett:2004pv,Bennett:2006fi,Hertzog:2008zz} and for the determination of the
value of the running QED coupling at the $Z^o$
resonance, $\alpha (m_{Z}^{2})$, which contributes to precision tests of
the \emph{Standard Model} of particle physics, for details see e.g. \emph{Jegerlehner} \cite{Jegerlehner:2008zza}, also \emph{Davie} and
\emph{Marciano} \cite{Davier:2004gb}, or \emph{Teubner et al.}
\cite{Teubner:2008zz,Jegerlehner:2008zz,Stockinger:2008zz}. 
The hadronic contribution to $a_\mu$ below about 2 GeV is
dominated by the $2 \pi$ final state, which contributes about 70\% 
due to the dominance of the $\rho-$resonance. Other major
contributions come from the three- and four-pion final states. These
hadronic final states constitute at present the largest error to
the \emph{Standard Model} values of $a_\mu$ and $\alpha
(m_{Z}^{2})$ and can be determined only experimentally. This is because
calculations within perturbative QCD are unrealistic,
calculations on the lattice are not yet available with the necessary
accuracy, and calculations in the framework of chiral perturbation
theory are restricted to values close to the reaction thresholds.
At energies above about 2 to 2.5 GeV, perturbative QCD
calculations start to become possible and reliable, see e.g. Refs.
\cite{Kuhn:1998ze,Eidelman:1998vc}, and also
\cite{Chetyrkin:1996ia}.

The Novosibirsk groups CMD-2 \cite{Budker,Akhmetshin:2001ig,Akhmetshin:2002vj,Akhmetshin:2003ag,Akhmetshin:2003zn,Aulchenko:2006na,Akhmetshin:2004gw,Akhmetshin:2006wh,Ignatov:2008zz,Akhmetshin:2006bx}
and SND \cite{Achasov:2003ir,Achasov:2005rg,Achasov:2006xc,Achasov:2006vp,Achasov:2007kg,Achasov:2006bv} measured hadronic cross sections below
1.4 GeV by changing the collider energy (\emph{energy scan}, see
the preceding Section \ref{sec:2}). The \emph{Initial State Radiation}
method used by KLOE represents an alternative, independent
and complementary way to determine hadronic cross sections with
different systematic errors. KLOE has determined the cross
section for the reaction $e^{+} e^{-}\rightarrow \pi^{+} \pi^{-}$
in the energy region between 0.63 and 0.958 GeV by measuring the
reaction $e^{+} e^{-}\rightarrow \pi^{+} \pi^{-} \gamma$ and
applying a radiator function based on PHOKHARA. 
For the hadronic contribution to the anomalous magnetic moment of
the muon due to the $2 \pi$ final state it obtained $a_{\mu}^{\pi\pi} = (356.7
\pm 3.1_{{\rm stat}+{\rm syst}}) \cdot 10^{-10}$ \cite{:2008en}. This value
is in good agreement with those from SND \cite{Achasov:2006bv} and
CMD-2 \cite{Akhmetshin:2006bx}, $a_{\mu}^{\pi\pi} = (361.0 \pm
5.1_{{\rm stat}+{\rm syst}}) \cdot 10^{-10}$ and $a_{\mu}^{\pi\pi} = (361.5
\pm 3.4_{stat+syst})\cdot 10^{-10}$, respectively, leading to an evaluation of 
$a_{\mu}$~\cite{Jegerlehner:2008zza,Davier:2004gb,Teubner:2008zz,Jegerlehner:2008zz,Stockinger:2008zz,Davier:2009ag} which differs
%differing
by about three standard deviations from the BNL experiment~\cite{Bennett:2006fi}.
A different evaluation using $\tau$ decays into
two pions results in a reduced discrepancy~\cite{Davier:2004gb,Davier:2009ag}.
The difference between $e^{+}e^{-}$ and $\tau$ based analyses
 is at present not understood.  
But one has to be aware that the evaluation with $\tau$ data needs more theoretical input.

%, however, from the value for
%$a_{\mu}^{\pi\pi}$ obtained by the analysis of $\tau$ -decays into
%2 pions according to $\tau^{-} \rightarrow
%\pi^{-}\pi^{o}\nu_{\tau}$
%\cite{Jegerlehner:2008zza,Davier:2004gb,Teubner:2008zz,Jegerlehner:2008zz,Stockinger:2008zz,Davier:2009ag} which, however, must be
%corrected for isospin symmetry breaking.

Soon after the application of EVA-PHOKHARA to KLOE
\cite{Cataldi:1999aa}, the BaBar collaboration also started the
measurement of hadronic cross sections exploiting ISR
\cite{Solodov:2002xu} and using PHOKHARA (Section
\ref{rr:babar}). In recent years a
plethora of final states has been studied, starting with the
reaction $e^+e^- \rightarrow J/\psi \: \gamma \rightarrow \mu^{+}
\mu^{-}\: \gamma$ \cite{Aubert:2003sv}. While detecting a hard
photon, the upper energy for the hadron cross sections is limited
to roughly 4.5 GeV. Final states with 3, 4, 5, 6 charged and
neutral pions, 2 pions and 2 kaons, 4 kaons, 4 pions and 2 kaons,
with a $\phi$ and an $f_{o}(980)$, $J/\psi$ and 2 pions or 2 kaons,
pions and $\eta$, kaons and $\eta$, but also baryonic final states
with protons and antiprotons, $\Lambda^{o}$ and $\bar
\Lambda^{o}$, $\Lambda^{o}$ and $\bar {\Sigma^{o}}$, $\Sigma^{o}$
and $\bar {\Sigma^{o}}$, $D \bar D$, $D \bar D^*$, and $D^* \bar
D^*$ mesons, etc. have been investigated \cite{Aubert:2004kj,Aubert:2005eg,Aubert:2005cb,Aubert:2006jq,Aubert:2006bu,Aubert:2007ur,Aubert:2007ym,Aubert:2007ef,Aubert:2007uf,Aubert:2006mi,:2008ic,:2009xs,Denig:2008zz}.
In preparation are final states with 2 pions \cite{Davier:2009aa}
and 2 kaons. Particularly important final states are those with 4
pions (including $\omega \pi^o$). They contribute significantly to
the muon anomalous magnetic moment and were poorly known
before the ISR measurements. In many of these channels
additional insights into isospin symmetry breaking are expected from
the comparison between $e^+ e^-$ annihilation and $\tau$ decays.

More recently also Belle joined the ISR programme
with emphasis on final states containing mesons with hidden and
open charm: $J/\psi$ and $\psi(2S)$, 
%$D^*$ and $\bar D^*$,
$D^{(*)}$ and $\bar D^{(*)}$,
$\Lambda_c{^+} \Lambda_c{^-}$ 
%\cite{Abe:2006fj,:2007sj,:2007ea,Pakhlova:2008zza,:2007bt,Pakhlova:2007fq,Pakhlova:2008vn,Pakhlova:2009jv} 
\cite{:2007sj,:2007ea,:2007bt,Abe:2006fj,Pakhlova:2008zza,Pakhlova:2007fq,Pakhlova:2008vn,Pakhlova:2009jv} 
(Section \ref{rr:belle}).



A major surprise in recent years was the opening of a totally new
field of hadron spectroscopy by applying ISR. Several new, 
relatively narrow highly excited states with $J^{PC}= 1^{--}$, the
quantum numbers of the photon, have been discovered (preliminarily
denoted as \emph{X, Y, Z}) at the $B$ factories PEP-II and
KEKB with the detectors BaBar and Belle,
respectively. The first of them was found by BaBar in the
reaction $e^+e^-\rightarrow Y(4260)\: \gamma \rightarrow J/\psi \:
\pi^{+} \pi^{-} \gamma$ \cite{Aubert:2005rm},
 a state around 4260 MeV with
a width of 90 MeV, later confirmed by Belle via ISR
%\cite{Abe:2006hf,:2007sj,Balagura:2008zz} 
\cite{Abe:2006hf,:2007sj} and by CLEO in an direct energy scan
\cite{Coan:2006rv} and a radiative return \cite{He:2006kg}. Another state was detected at 2175 MeV by
BaBar in the reaction $e^+e^-\rightarrow Y(2175)\: \gamma
\rightarrow \phi f_{o}(980)\gamma$ \cite{Aubert:2006bu}. Belle
found new states at 4050, 4360, 4660 MeV in the reactions
$e^+e^-\rightarrow Y\: \gamma \rightarrow J/\psi \: \pi^{+}
\pi^{-}\gamma$ and $e^+e^-\rightarrow Y\: \gamma \rightarrow
\psi(2S)\: \pi^{+} \pi^{-}\gamma$ 
%\cite{Coan:2006rv,:2007ea,:2007sj}
\cite{:2007sj,:2007ea}. The structure of
basically all of these new states (if they will survive) is
unknown so far. Four-quark states, e.g. a $[cs][\bar{c} \bar{s}]$
state for $Y(4260)$, a $[ss][\bar{s}\bar{s}]$ state for $Y(2175)$,
hybrid and molecular structures are discussed, see also
\cite{Kalashnikova:2008zz}.

Detailed analyses allow, in addition, also the identification of
intermediate states, and consequently a study of reaction
mechanisms. For instance, in the case of the final state with 2
charged and 2 neutral pions $(e^+ e^- \rightarrow \pi^{+} \pi^{-}
\pi^{o} \pi^{o} \gamma)$, the dominating intermediate states are
$\omega \pi^{o}$ and $a_{1}(1260) \pi $, while $\rho^{+} \rho^{-}$
and $\rho^{o} f_{o}(980)$ contribute significantly less.

Many more highly excited states with quantum numbers different
from those of the photon have been found in decay chains of the
primarily produced heavy mesons at the $B$ factories
PEP-II and KEKB. These analyses without ISR
have clearly been triggered and encouraged by the unexpected
discovery of highly excited states with $J^{PC}= 1^{--}$ found
with ISR.

Also baryonic final states with protons and antiprotons,
$\Lambda^{o}$ and $\bar \Lambda^{o}$, $\Lambda^{o}$ and $\bar
{\Sigma^{o}}$, $\Sigma^{o}$ and $\bar {\Sigma^{o}}$ have been
investigated using ISR. The effective proton form factor (see Section~\ref{rr:babar}) shows a strong increase down to the $p\bar{p}$ threshold and
nontrivial structures at invariant $p\bar{p}$ masses of 2.25 and
3.0 GeV, so far unexplained \cite{Aubert:2005cb,Maas:2008zza,Salme:2008an,Dmitriev:2008zz,Baldini:2007qg}. Furthermore, it
should be possible to disentangle electric and magnetic form
factors and thus shed light on discrepancies between different
measurements of these quantities in the space-like region
\cite{Arrington:2003df}.

Prospects for the \emph{Radiative Return} at the
Novosibirsk collider VEPP-2000 and BEPC-II are
discussed in Sections \ref{rr:vepp2000} and \ref{rr:besiii}.

\subsection{Radiative return: a theoretical overview}
\label{radret:theo}


\subsubsection{Radiative return at leading order}
\label{sec:LO}

We consider the $e^+ e^-$ annihilation process
\beq
e^+(p_1) + e^-(p_2) \rightarrow {\rm hadrons} + \gamma(k_1)~,
\label{eq:LO}
\eeq
where the real photon is emitted either from the initial 
(Fig.~\ref{fig:born}a) or the final state (Fig.~\ref{fig:born}b).
The former process is denoted initial state radiation (ISR), 
while the latter is called final state radiation (FSR).

The differential rate for the ISR process 
can be cast into the product of a leptonic $L^{\mu\nu}$ 
and a hadronic $H^{\mu\nu}$ tensor and the corresponding 
factorised phase space
\bea
{\rm d}\sigma_{\rm ISR} &=& \frac{1}{2s} L^{\mu \nu}_{\rm ISR} H_{\mu \nu}
\nn \\ && \times {\rm d} \Phi_2(p_1,p_2;Q,k_1) {\rm d} \Phi_n(Q;q_1,\cdot,q_n) 
\frac{{\rm d}Q^2}{2\pi}~,
\eea
where ${\rm d} \Phi_n(Q;q_1,\cdot,q_n)$ denotes the 
hadronic $n$-body phase-space with all the statistical factors 
coming from the hadro\-nic final state included, 
$Q = \sum q_i$ and $s=(p_1+p_2)^2$.

\begin{figure}[h]
\begin{center}
\epsfig{file=diagram_ifs_1.ps,width=8.5cm} 
\caption{Leading order contributions to the reaction
$e^+e^-\to h \, \bar h + \gamma$ from ISR (a) and FSR (b).
Final state particles are pions or muons, or any other
multi-hadron state. The blob represents the hadronic 
form factor.}
\label{fig:born}
\end{center}
\end{figure}


For an arbitrary hadronic final state, the matrix element for 
the diagrams in Fig.~\ref{fig:born}a is given by  
\begin{align}
{\cal A}_{\rm ISR}^{(0)} 
&= M^{(0)}_{\rm ISR} \cdot J^{(0)} =  \non \\
&= - \frac{e^2}{Q^2}  \bar{v}(p_1) \bigg(   
\frac{\ta{\varepsilon}^*(k_1)[\ta{k}_1-\ta{p}_1+m_e]
\gamma^{\mu}}{2 k_1 \cdot p_1} \non \\ 
& \qquad + \frac{\gamma^{\mu}[\ta{p}_2-\ta{k}_1+m_e]\ta{\varepsilon}^*(k_1)}
{2 k_1 \cdot p_2} 
\bigg) u(p_2) \; J_{\mu}^{(0)}~, 
\label{AISR}
\end{align}
where $J_{\mu}$ is the hadronic current.
The superscript $(0)$ indicates that the scattering amplitude 
is evaluated at tree-level. 
Summing over the polarisations of the final real photon,
averaging over the polarisations of the initial $e^+ e^-$ state,
and using current conservation, $Q \cdot J^{(0)} = 0$, 
the leptonic tensor
\begin{equation*}
L_{\rm ISR}^{(0), \mu \nu} = \overline{M_{\rm ISR}^{(0), \, \mu} 
( M_{\rm ISR}^{(0), \, \nu} )^\dagger}
\end{equation*}
can be written in the form
\begin{align}
L_{\rm ISR}^{(0),\, \mu \nu} &= 
\frac{(4 \pi \alpha)^2}{Q^4} \; \bigg[ \left( 
\frac{2 m^2 q^2(1-q^2)^2}{y_1^2 y_2^2}
- \frac{2 q^2+y_1^2+y_2^2}{y_1 y_2} \right) g^{\mu \nu} \non \\ & 
+ \left(\frac{8 m^2}{y_2^2} - \frac{4q^2}{y_1 y_2} \right) 
\frac{p_1^{\mu} p_1^{\nu}}{s} 
+ \left(\frac{8 m^2}{y_1^2} - \frac{4q^2}{y_1 y_2} \right) 
\frac{p_2^{\mu} p_2^{\nu}}{s} \non \\
& - \left( \frac{8 m^2}{y_1 y_2} \right) 
\frac{p_1^{\mu} p_2^{\nu} + p_1^{\nu} p_2^{\mu}}{s} \bigg]~, 
\label{Lmunu0}
\end{align}
with 
\begin{equation}
y_i = \frac{2 k_1 \cdot p_i}{s}~, 
\qquad m^2=\frac{m_e^2}{s}~, \qquad q^2=\frac{Q^2}{s}~.
\label{dimensionless}
\end{equation}
The leptonic tensor is symmetric under the exchange of the electron and 
the positron momenta. Expressing the bilinear products $y_i$ 
by the photon emission angle in the c.m. frame,
\begin{equation*}
y_{1,2} = \frac{1-q^2}{2}(1 \mp \beta \cos \theta)~, 
\qquad \beta = \sqrt{1-4m^2}~,
\end{equation*}
and rewriting the two-body phase space as 
\begin{equation}
{\rm d} \Phi_2(p_1,p_2;Q,k_1) = \frac{1-q^2}{32 \pi^2} {\rm d} \Omega~,
\end{equation}
it is evident that expression (\ref{Lmunu0}) contains several 
singularities: soft singularities for $q^2\rightarrow 1$ and 
collinear singularities for $\cos \theta \rightarrow \pm 1$.
The former are avoided by requiring a minimal photon energy.
The latter are regulated by the electron mass.
For $s \gg m_e^2$ the expression (\ref{Lmunu0}) can  
nevertheless be safely taken in the limit $m_e\rightarrow 0$ if the 
emitted real photon lies far from the collinear region.
In general, however, one encounters spurious singularities in the 
phase space integrations if powers of $m^2=m_e^2/s$ are 
neglected prematurely.

Physics of the hadronic system, whose description is model dependent, 
enters through the hadronic tensor 
\begin{equation}
H_{\mu \nu} = J^{(0)}_{\mu} (J^{(0)}_{\nu})^\dagger~,
\end{equation}
where the hadronic current has to be parametrised through form factors.
For two charged pions in the final state, the current 
\begin{equation}
J^{(0), \, \mu}_{\pi^+\pi^-} = i e F_{2\pi}(Q^2) \; (q_1-q_2)^{\mu}~,
\end{equation}
where $q_1$ and $q_2$ are the momenta of the $\pi^+$ and
$\pi^-$, respectively, is determined by only one function,
the pion form factor $F_{2\pi}$.
The current for the \(\mu^+\mu^-\) final state is 
obviously defined by QED:
\begin{eqnarray}
 J^{(0), \, \mu}_{\mu^+\mu^-} = i e \, \bar u(q_2)\gamma^\mu  v(q_1)~.
\end{eqnarray}



Integrating the hadronic tensor over the hadronic 
phase space, one gets
\begin{equation}
\int H^{\mu \nu} {\rm d}\Phi_n(Q;q_1, \cdot, q_n) 
= \frac{e^2}{6\pi} (Q^{\mu}Q^{\nu} - g^{\mu \nu} Q^2) R(Q^2)~,
\end{equation}
where $R(Q^2) = \sigma(e^+ e^- \rightarrow {\rm hadrons})/
\sigma_0(e^+ e^- \rightarrow \mu^+ \mu^-)$, with 
\beq
\sigma_0(e^+ e^- \rightarrow \mu^+ \mu^-) = \frac{4\pi \, \alpha^2}{3 Q^2}
\eeq
the tree-level muonic cross section in the limit $Q^2\gg 4m_{\mu}^2$.
After the additional integration over the photon angles,  
the differential distribution 
\begin{align}
Q^2 \frac{{\rm d}\sigma_{\rm ISR}}{{\rm d}Q^2} = \frac{4\alpha^3}{3 s} R(Q^2)
\left\{ \frac{s^2+Q^4}{s(s-Q^2)} \left( L - 1 \right) \right\}~,
\label{diff1n}
\end{align}
with $L=\log(s/m_e^2)$ is obtained. 
If instead the photon polar angle is restricted to be 
in the range $\theta_{\rm min}< \theta < \pi-\theta_{\rm min}$,
this differential distribution is given by 
\begin{align}
Q^2 \frac{{\rm d}\sigma_{\rm ISR}}{{\rm d}Q^2} &= \frac{4\alpha^3}{3 s} R(Q^2)
\bigg\{ \frac{s^2+Q^4}{s(s-Q^2)} 
\log \frac{1+\cos \theta_{\rm min}}{1-\cos \theta_{\rm min}} \non \\
& - \frac{s-Q^2}{s} \cos \theta_{\rm min} \bigg\}~.
\label{diff2n}
\end{align}

\begin{figure}[th]
\begin{center}
\epsfig{file=isrfsr_1gev.ps,width=8.5cm}
\caption{Suppression of the FSR contributions
to the cross section by a suitable choice of angular cuts; results from the PHOKHARA generator; no cuts (upper curves) and suitable cuts applied (lower curves). }
\label{fig:isrtofsr}
\end{center}
\end{figure}                                                                   


In the latter case, the electron mass can be taken equal to zero
before integration, since the collinear region is excluded
by the angular cut. The contribution of the two-pion exclusive
channel can be calculated from \Eq{diff1n} and \Eq{diff2n} with 
\beq
R_{\pi^+\pi^-}(Q^2) = \frac14 \left(1-\frac{4m_{\pi}^2}{Q^2}\right)^{3/2} 
|F_{2\pi}(Q^2)|^2~,
\eeq 
and the corresponding muonic contribution with 
\beq
R_{\mu^+\mu^-}(Q^2) = \sqrt{1-\frac{4m_{\mu}^2}{Q^2}} 
\left(1+\frac{2m_\mu^2}{Q^2}\right)~.
\eeq


A potential complication for the measurement of the hadronic cross section 
from the radiative return may arise from the interplay between photons 
from ISR and FSR \cite{Binner:1999bt}.
Their relative strength is strongly dependent
on the photon angle relative to the beam and to the direction of the final 
state particles, the c.m. energy of the reaction and the invariant mass 
of the hadronic system. While ISR is independent of the hadronic final state,
FSR is not. Moreover, it cannot be predicted from first principles and 
thus has to be modelled.


\begin{figure*}[ht]
\begin{center}
\epsfig{file=angularfsr.ps,width=8cm,height=6cm} 
\epsfig{file=mangularfsr.ps,width=8cm,height=6cm} 
\end{center}
\caption{Angular distributions of $\pi^+$ and $\mu^+$ at 
$\sqrt{s}=1.02$~GeV with and without FSR 
for different angular cuts.}
\label{fig:angular}
\end{figure*}

\begin{figure*}[ht]
\begin{center}
\epsfig{file=angular10gev.ps,width=8.5cm,height=6cm} 
\epsfig{file=mangular10gev.ps,width=8.5cm,height=6cm} 
\end{center}
\caption{Angular distributions of $\pi^+$ 
(ISR \(\simeq\) FSR+ISR) and $\mu^+$ at $\sqrt{s}$=10.6~GeV for various $Q^2$ cuts.}
\label{fig:angular10GeV}
\end{figure*}


The amplitude for FSR (Fig.~\ref{fig:born}b) factorises as well as
\beq
{\cal A}_{\rm FSR}^{(0)} = M^{(0)} \cdot J_{\rm FSR}^{(0)}~,
\label{AFSR}
\eeq
where
\beq
M^{(0)}_\mu = \frac{e}{s} \, \bar{v}(p_1) \gamma_{\mu} u(p_2)~.
\eeq
Assuming that pions are point-like, the FSR current for 
two pions in scalar QED (sQED) reads  
\bea
J_{\rm FSR}^{(0), \, \mu} 
&=& - i \, e^2 \, F_{2\pi}(s) \nn \\ &\times&  
\left[ -2 g^{\mu\sigma} 
+ (q_1+k_1-q_2)^{\mu} \, \frac{(2q_1+k_1)^\sigma }{2k_1 \cdot q_1}
\right. \nn \\ && \left.
- (q_1-k_1-q_2)^{\mu} \, \frac{(2q_2+k_1)^\sigma }{2k_1 \cdot q_2}
\right] \epsilon_\sigma^*(k_1)~.  
\eea
Due to momentum conservation, $p_1+p_2=q_1+q_2+k_1$, and 
current conservation, this expression can be simplified further to 
\bea
J_{\rm FSR}^{(0), \, \mu} &=& 2 i \, e^2 \, F_{2\pi}(s) \,
\left[ g^{\mu\sigma} 
+ \frac{q_2^{\mu}\, q_1^\sigma}{k_1 \cdot q_1}
+ \frac{q_1^{\mu}\, q_2^\sigma}{k_1 \cdot q_2}
\right] \epsilon_\sigma^*(k_1)~. \nn \\
\eea
This is the basic model adopted in EVA~\cite{Binner:1999bt} 
and in PHO\-KHARA \cite{Rodrigo:2001jr,Kuhn:2002xg,Rodrigo:2001kf,Czyz:2002np,Czyz:2003ue,Czyz:2004rj,Czyz:2004nq,Czyz:2005as,Czyz:2007wi,Rodrigo:2001cc} 
to simulate FSR off charged pions.
The corresponding FSR current for muons is 
given by QED. 


The fully differential cross section describing photon emission 
at leading order can be split into three pieces
\begin{equation}
{\rm d}\sigma^{(0)} = {\rm d}\sigma_\mathrm{ISR}^{(0)} 
+ {\rm d}\sigma_\mathrm{FSR}^{(0)} + {\rm d}\sigma_\mathrm{INT}^{(0)}~,
\label{LOxsection}
\end{equation}
which originate from the squared ISR and FSR amplitudes and the interference term, respectively. 
The ISR--FSR interference is odd under charge conjugation, 
\beq
{\rm d}\sigma_{\rm INT}^{(0)}(q_1,q_2) = - {\rm d}\sigma_{\rm INT}^{(0)}(q_2,q_1)~, 
\eeq
and its contribution vanishes after angular integration.
It gives rise, however, to a relatively large charge asymmetry and, 
correspondingly, to a forward--backward asymmetry
\begin{equation}
A(\theta) = \frac{N_{h}(\theta)-N_{h}(\pi-\theta)}
{N_{h}(\theta)+N_{h}(\pi-\theta)}~.
\end{equation}
The asymmetry can be used for the calibration of the FSR amplitude,
and fits to the angular distribution $A(\theta)$ can test
details of its model dependence~\cite{Binner:1999bt}. 

The second option to disentangle ISR from FSR exploits
the markedly different angular distribution of the photon 
from the two processes. This observation is completely general
and does not rely on any model like sQED for FSR.
FSR is dominated by photons collinear to the final state 
particles, while ISR is dominated by photons collinear to 
the beam direction. This suggests that we should consider 
only events with photons well separated from the charged 
final state particles and preferentially close to the 
beam~\cite{Binner:1999bt,Rodrigo:2001kf,Czyz:2002np}.
 

\begin{figure}[ht]
\begin{center}
\epsfig{file=diagram_kinematics.ps,width=8.5cm} 
\end{center}
\caption{Typical kinematic configuration of the radiative return at low 
and high energies. }
\label{fig:kinematics}
\end{figure}


This is illustrated in Fig.~\ref{fig:isrtofsr}, 
which has been generated running PHOKHARA at leading order (LO). 
After introducing suitable angular cuts, the contamination of 
events with FSR is easily reduced to less than a few per mill.
The price to pay, however, is a suppression of the threshold 
region too. To have access to that region, 
photons at large angles need to be tagged and a better 
control of FSR is required. In Fig.~\ref{fig:angular}
the angular distribution of $\pi^+$ and $\mu^+$ at 
DA$\mathrm{\Phi}$NE energies, $\sqrt{s}=1.02$~GeV, are shown 
for different angular cuts. The angles are defined with respect 
to the incoming positron. If no angular cut is applied, 
the angular distribution in both cases is highly asymmetric as 
a consequence of the ISR--FSR interference contribution. If cuts 
suitable to suppress FSR, and therefore the ISR--FSR interference,
are applied, the distributions become symmetric. 

Two complementary analyses are therefore possible (for details see 
Section~\ref{rr:kloe}): 
The {\it small photon angle} analysis, where the photon is untagged 
and FSR can be suppressed below some reasonable limit. 
This analysis is suitable 
for intermediate values of the invariant mass of the hadronic 
system. And the {\it large photon angle} analysis, giving access to 
the threshold region, where FSR is more pronounced and 
the charge asymmetry is a useful tool to probe its model
dependence. 

These considerations apply, however, only to low beam energies,
around $1$~GeV. 
At high energies, e.g. at  $B$ facto\-ries, very hard
tagged photons are needed to access the region with low 
hadronic invariant masses, and the hadronic system is mainly 
produced back-to-back to the hard photon. The suppression of FSR 
is naturally accomplished and no special angular cuts are 
needed. This kinematical situation is illustrated 
in Fig.~\ref{fig:kinematics}.  The suppression of FSR contributions 
to $\pi^+\pi^-\gamma$ events is also a consequence of the rapid 
decrease of the form factor above $1$~GeV.
The relative size of FSR is of the order of a few per mill
(see Fig.~\ref{fig:angular10GeV}). 
For $\mu^+\mu^-$ in the final state, the amount of FSR 
depends on the invariant mass of the muons. 
For $\sqrt{Q^2} < 1$~GeV FSR is still tiny, 
but becomes more relevant for larger values of $Q^2$
(see Fig.~\ref{fig:angular10GeV}). 



\subsubsection{Structure functions}
The original and default version of EVA~\cite{Binner:1999bt}, simulating 
the process $e^+ e^- \rightarrow \pi^+ \pi^- \gamma$ at LO, allowed for 
additional initial state radiation of soft and collinear photons by the 
structure function (SF) method~\cite{Caffo:1997yy,Caffo:1994dm}.

In the leading logarithmic approximation (LL), the multiple 
emission of collinear photons off an electron is described 
by the convolution integral 
\beq
\sigma(e^-X\to Y+ n\gamma) = \int_0^1 {\rm d}x \, f_e(x,Q^2) \, 
\sigma(e^- X\to Y)~,
\eeq
where $f_e(x,Q^2)$ is the probability distribution of the electron 
with longitudinal momentum fraction $x$, and $Q$ is the transverse 
momentum of the collinear photons. The function $f_e(x,Q^2)$ fulfils 
the evolution equation 
\bea
&& \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!
\frac{{\rm d}}{{\rm d} \log Q} f_e(x,Q^2) = \int_x^1
\frac{{\rm d} z}{z} \nn \\ && \frac{\alpha}{\pi}
\left( \frac{1+z^2}{(1-z)_+} + \frac32 \delta(1-z)\right)
f_e(\frac{x}{z},Q^2)
\label{evoleq}
\eea
with initial conditions
\beq
\left. f_e(x,Q^2)\right|_{Q^2=m_e^2} = \delta(1-x)~,
\eeq
and the $+$ prescription defined as 
\beq
\int_0^1 {\rm d}x \, \frac{f(x)}{(1-x)_+}
= \int_0^1 {\rm d}x \, \frac{f(x)-f(1)}{(1-x)}~.
\eeq

The analytic solution to \Eq{evoleq} provided in 
Refs.~\cite{Caffo:1997yy,Caffo:1994dm} allows to resum soft 
photons to all orders in perturbation theory, accounting for  
large logarithms of collinear origin, $L=\log(s/m_e^2)$, up to two loops.
The resummed cross section, 
\beq
\sigma_{\rm SF} =
\int_0^1 {\rm d}x_1 \int_0^1  {\rm d}x_2 \, D(x_1) \, D(x_2) 
\, \sigma_{e^+e^-\to {\rm had.}+\gamma}(x_1 x_2 s)~,
\eeq
is thus obtained by convoluting 
the Born cross section of the hard photon emission 
process $e^+e^- \to {\rm hadrons} +\gamma$ 
with the SF distribution~\cite{Caffo:1997yy,Caffo:1994dm} 
\bea
D(x) &=& \left[1+\delta_N  \right]^{1/2} \, 
\frac{\beta_e}{2} (1-x)^{\frac{\beta_e}{2}-1} \nn \\
&& \times\bigg\{ \frac12 (1+x^2) + 
\frac12 \frac{(1-x)^2}{L-1} \nn \\ 
&& + \frac{\beta_e}{8}\left( -\frac{1}{2}(1+3x^2)
\log x - (1-x)^2 \right) \bigg\}~,
\eea
with 
\beq
\beta_e = 2 \, \frac{\alpha}{\pi} \, (L-1)
\eeq
and 
\bea
\delta_N &=& \frac{\alpha}{\pi} \left( \frac{3}{2}L+\frac{\pi^2}{3}-2\right)
\nn \\ &&
+ \beta_e^2 \frac{\pi^2}{8} + \left(\frac{\alpha}{\pi}\right)^2
\left( \frac{11}{8}-\frac{2\pi^2}{3}\right)L^2~. 
\eea

In the SF approach, the additional emission of collinear 
photons reduces the effective c.m. energy of the collision
to $\sqrt{x_1 x_2 s}$.
Momentum conservation is not accomplished because the 
extra radiation is integrated out. 
In order to reduce the kinematic distortion of the events,
a minimal invariant mass of the observed particles, 
hadrons plus the tagged photon, was required in~\cite{Binner:1999bt},
introducing in turn a cut dependence.
Therefore the SF predictions are not accurate enough for 
a high precision measurement of the hadronic cross section 
from radiative return, and a next-to-leading order (NLO) calculation  
is in order. The NLO prediction contains the large logarithms  
$L=\log(s/m_e^2)$ at order $\alpha^3$ and additional sub-leading terms, 
which are not taken into account within the SF method. Furthermore, 
it allows for a better control of the kinematical configurations
because momentum conservation is fulfilled. 
A comparison between SF and NLO predictions can be found 
in~\cite{Rodrigo:2001kf}.

\begin{figure}[ht]
\epsfig{file=diagram_ifs_20.ps,width=8cm} 
\caption{Typical sub-amplitudes describing virtual and real 
corrections to the reaction $e^+e^-\to h \bar h + \gamma (\gamma)$, 
where $h = \pi^-$, $\mu^-$. The superscripts $(0)$ and $(1)$ denote
tree-level and one-loop quantities, respectively. 
ISR and FSR indicate that real photons are emitted from the initial 
or final state. The last two diagrams, with exchange of 
two virtual photons, are non-factorisable. 
Permutations are omitted.}
\label{fig:subamplitudes}
\end{figure}



 
\subsubsection{Radiative return at NLO}
\label{rratnlo}

At NLO, the $e^+ e^-$ annihilation process in~\Eq{eq:LO}
receives contributions from one-loop corrections 
and from the emission of a second real photon 
(see Fig.~\ref{fig:subamplitudes}). 
After renormalisation, the one-loop matrix elements still contain  
infrared divergences. These are cancelled by adding the two-photon contributions to the one-loop corrections. 
There are several well established methods to perform this cancellation. 
The slicing method, where amplitudes are evaluated in dimensional 
regularisation and the two photon contribution is integrated 
analytically in phase space for one of the photon energies
up to an energy cutoff $E_{\gamma}<w\sqrt{s}$ far below $\sqrt{s}$,
was used in \cite{Rodrigo:2001jr,Kuhn:2002xg} 
to calculate the NLO corrections to ISR. 
Here the sum of the virtual and soft contributions is finite, but it depends on the soft photon cutoff. 
The contribution from the emission of the second photon with 
energy $E_{\gamma}>w\sqrt{s}$, which is evaluated numerically, 
completes the calculation and cancels this dependence. 

The size and sign of the NLO corrections do depend on the 
particular choice of the experimental cuts. Hence, only using a 
Monte Carlo event generator one can realistically 
compare theoretical predictions with experiment.
This is the main motivation behind PHOKHARA~\cite{Rodrigo:2001jr,Kuhn:2002xg,Rodrigo:2001kf,Czyz:2002np,Czyz:2003ue,Czyz:2004rj,Czyz:2004nq,Czyz:2005as,Czyz:2007wi,Rodrigo:2001cc}. 

The full set of scattering amplitudes at tree-level and one-loop 
can be constructed from the sub-amplitudes 
depicted in Fig.~\ref{fig:subamplitudes}.
The one-loop amplitude with emission of a single photon is given by 
\bea
 {\cal A}^{(1)}_{1\gamma} &=& 
  {\cal A}^{(1)}_{\rm ISR} + {\cal A}^{(1)}_{\rm FSR} \nn \\ &+&
  M^{(1)}\cdot J^{(0)}_{\rm FSR} 
+ M^{(0)}_{\rm ISR}\cdot J^{(1)}   \nn \\ &+&
  {\cal A}^{2\gamma^*}_{\rm ISR} + {\cal A}^{2\gamma^*}_{\rm FSR}~,
\eea
where
\beq
{\cal A}^{(1)}_{\rm ISR} = M^{(1)}_{\rm ISR}\cdot J^{(0)}~, \qquad
{\cal A}^{(1)}_{\rm FSR} = M^{(0)}\cdot J^{(1)}_{\rm FSR}~,
\eeq
while the amplitude with emission of two real photons reads 
\bea
{\cal A}^{(0)}_{2\gamma} &=&  
  {\cal A}^{(0)}_{\rm 2ISR}
+ {\cal A}^{(0)}_{\rm 2FSR} \nn \\ &+&
  \left( M^{(0)}_{\rm ISR}(k_1)\cdot J^{(0)}_{\rm FSR}(k_2)  
+ (k_1 \leftrightarrow k_2) \right)~,   
\eea
where
\beq
{\cal A}^{(0)}_{\rm 2ISR} = 
M^{(0)}_{\rm 2ISR}\cdot J^{(0)}~, \qquad
{\cal A}^{(0)}_{\rm 2FSR} = 
M^{(0)}\cdot J^{(0)}_{\rm 2FSR}~.
\eeq

PHOKHARA includes the full LO amplitudes and 
the most relevant C-even NLO contributions: 
\beq
{\rm d}\sigma = {\rm d}\sigma^{(0)} + {\rm d}\sigma^{(1)}_{\rm ISR} 
+ {\rm d}\sigma^{(1)}_{\rm IFS}~,
\eeq
where ${\rm d}\sigma^{(0)}$ is the LO differential cross section 
(\Eq{LOxsection}), 
\bea
{\rm d}\sigma^{(1)}_{\rm ISR} &=& \frac{1}{2s} \bigg[
2 {\rm Re} \left\{   
{\cal A}^{(1)}_{\rm ISR}
\left({\cal A}^{(0)}_{\rm ISR}\right)^\dagger \right\} \, 
{\rm d}\Phi_3(p_1,p_2;q_1,q_2,k_1) \nn \\
&& + \left|{\cal A}^{(0)}_{\rm 2ISR}\right|^2 \,
{\rm d}\Phi_4(p_1,p_2;q_1,q_2,k_1,k_2) \bigg]
\eea
is the second order radiative correction to ISR, and 
\bea
{\rm d}\sigma^{(1)}_{\rm IFS} &=& \frac{1}{2s} \bigg[
2 {\rm Re} \bigg\{   
M^{(0)}_{\rm ISR}\cdot J^{(1)} 
\left({\cal A}^{(0)}_{\rm ISR}\right)^\dagger
\nn \\ &&
+ M^{(1)}\cdot J^{(0)}_{\rm FSR} \, 
\left({\cal A}^{(0)}_{\rm FSR}\right)^\dagger
\bigg\} \, {\rm d}\Phi_3(p_1,p_2;q_1,q_2,k_1) \nn \\
&& +
\left( \left|M^{(0)}_{\rm ISR}(k_1)\cdot J^{(0)}_{\rm FSR}(k_2)\right|^2 
+ (k_1 \leftrightarrow k_2) \right) \nn \\ && \times
\, {\rm d}\Phi_4(p_1,p_2;q_1,q_2,k_1,k_2) \bigg]
\label{sigmaIFS}
\eea
is the contribution of events with simultaneous emission of 
one photon from the initial state and another one from the 
final state, together with ISR amplitudes with final state
one-loop vertex corrections, and FSR amplitudes with 
initial state one-loop vertex corrections. We denote these
 corrections as IFS.

Vacuum polarisation corrections are included in the hadronic 
currents multiplicatively: 
\bea
 J^{(i)} &\to& C_{\rm VP}(Q^2) \, J^{(i)}~, \nn \\ 
 J^{(i)}_{\rm FSR}(k_j) &\to& 
C_{\rm VP}((Q+k_j)^2) \, J^{(i)}_{\rm FSR}(k_j)~, \nn \\ 
 J^{(0)}_{\rm 2FSR} &\to& 
C_{\rm VP}(s) \, J^{(0)}_{\rm 2FSR} ~. 
\eea
The virtual photon propagator is by definition included 
in the leptonic sub-amplitudes $M^{(i)}$, $M^{(i)}_{\rm ISR}$ 
and $M^{(0)}_{\rm 2ISR}$:
\bea
&& M^{(i)} \sim \frac{1}{s}~, \nn \\ 
&& M^{(i)}_{\rm ISR}(k_j) \sim \frac{1}{(p_1+p_2+k_j)^2}~, \nn \\ 
&& M^{(0)}_{\rm 2ISR} \sim \frac{1}{Q^2}~. 
\eea

Neither diagrams where two photons are emitted from the final state, nor
final-state vertex corrections with associated real radiation from
the final state are included.
These constitute radiative corrections to FSR and will give 
non-negligible contributions only for those cases where at least one 
photon is collinear with one of the final state particles.
Box diagrams with associated real radiation from the initial- or the 
final-state leptons, as well as pentagon diagrams, are also neglected.
As long as one considers charge symmetric
observables only, their contribution is divergent neither in the soft nor
the collinear limit and is thus of order $\alpha/\pi$ without any enhancement
factor. One should stress that PHOKHARA includes only C-even gauge 
invariant sets of diagrams at NLO. 
The missing contributions are either small or do not
contribute for charge symmetric cuts. Nevertheless their implementation 
is underway.

The calculation of the NLO corrections to ISR,  
${\rm d}\sigma^{(1)}_{\rm ISR}$, is independent of the final 
state.  These corrections are included by default
for all the final state channels implemented in PHOKHARA, 
and can be easily added for any other new channel, 
with the sole substitution of the tree-level final 
state current. The radiative corrections of the IFS process 
depend on the final state. The latest version of PHOKHARA 
(version 6.0~\cite{Czyz:2007wi})
includes these corrections for two charged pions, kaons 
and muons. 



{$\\$}
{\noindent \it Virtual and soft corrections to ISR\\}


The virtual and soft QED corrections to ISR in $e^+ e^-$ annihilation
were originally implemented in PHO\-KHA\-RA through the leptonic tensor. 
For future applications, however, it will be more convenient to 
implement those corrections directly at the amplitude level (in preparation).
In terms of sub-amplitudes, the leptonic tensor is given by  
\bea
 L^{\mu \nu}_{\rm ISR} &=& L^{(0),\mu \nu}_{\rm ISR} 
+ M^{(1),\, \mu}_{\rm ISR} \left(M^{(0),\, \nu}_{\rm ISR}\right)^\dagger 
+ M^{(0),\, \mu}_{\rm ISR} \left(M^{(1),\, \nu}_{\rm ISR}\right)^\dagger \nn \\
&& + \frac{1}{2(2\pi)^{d-1}} \int_0^{w\sqrt{s}} E^{d-3}\, {\rm d}E \, {\rm d}\Omega 
\, M^{(0),\, \mu}_{2\rm ISR} 
\left( M^{(0),\, \nu}_{2\rm ISR}\right)^\dagger~, \nn \\ &&
\eea
where $E$ and $\Omega$ are the energy and the solid angle of the soft photon,
respectively, and $d=4-2\epsilon$ is the number of dimensions in dimensional 
regularisation. 
The leptonic tensor has the general form 
\begin{align}
L^{\mu \nu}_{\rm ISR} &=
\frac{(4 \pi \alpha)^2}{Q^4} \;
\bigg[ a_{00} \; g^{\mu \nu} + a_{11} \; \frac{p_1^{\mu} p_1^{\nu}}{s}
 + a_{22} \; \frac{p_2^{\mu} p_2^{\nu}}{s} \non \\
&+ a_{12} \; \frac{p_1^{\mu} p_2^{\nu} + p_2^{\mu} p_1^{\nu}}{s}
+ i \pi \; a_{-1} \; 
\frac{p_1^{\mu} p_2^{\nu} - p_2^{\mu} p_1^{\nu}}{s} \bigg]~,
\label{generaltensor}
\end{align}
where the scalar coefficients $a_{ij}$ and $a_{-1}$ allow 
the following expansion: 
\begin{equation}
a_{ij} = a_{ij}^{(0)} + \frac{\alpha}{\pi} \; a_{ij}^{(1)}~, \qquad
a_{-1} = \frac{\alpha}{\pi} \; a_{-1}^{(1)}~.
\end{equation}
The imaginary antisymmetric piece, which is proportional to $a_{-1}$,
appears for the first time at second order and is particularly relevant for 
those cases where the hadronic current receives contributions from 
different amplitudes with nontrivial relative phases. This is 
possible, e.g., for final states with three or more mesons, or for 
$p\bar{p}$ production.

The LO coefficients $a_{ij}^{(0)}$ can be read directly from~\Eq{Lmunu0}
\begin{align}
a_{00}^{(0)} &= \frac{2 m^2 q^2(1-q^2)^2}{y_1^2 y_2^2}
- \frac{2 q^2+y_1^2+y_2^2}{y_1 y_2}~, \non \\
a_{11}^{(0)} &= \frac{8 m^2}{y_2^2} - \frac{4q^2}{y_1 y_2}~, \qquad
a_{22}^{(0)} = a_{11}^{(0)} (y_1 \leftrightarrow y_2)~, \non \\
a_{12}^{(0)} &= - \frac{8 m^2}{y_1 y_2} ~. 
\end{align}

The NLO coefficients $a_{ij}^{(1)}$ and $a_{-1}^{(1)}$ are 
obtained by combining the one-loop and the soft contributions. 
It is convenient to split the coefficients $a^{(1)}_{ij}$ 
into a part that contributes at large photon angles and a part 
proportional to $m_e^2$ and $m_e^4$ which is relevant only in the 
collinear regions. These coefficients are denoted by $a^{(1,0)}_{ij}$ 
and $a^{(1,m)}_{ij}$, respectively:
\begin{align}
a^{(1)}_{ij} &= a_{ij}^{(0)} \bigg[ -\log(4w^2) [1+\log(m^2)] \non \\ & 
-\frac{3}{2} \log(\frac{m^2}{q^2}) 
- 2 + \frac{\pi^2}{3} \bigg] 
+ a^{(1,0)}_{ij}+a^{(1,m)}_{ij}~.
\end{align}
The factor proportional to the LO coefficients $a^{(0)}_{ij}$ 
contains the usual soft and collinear logarithms.
The quantity $w$ denotes the dimensionless 
value of the soft photon energy cutoff, $E_{\gamma}<w\sqrt{s}$.
It is enough to present four out of the five coefficients 
because exchanging the positron with the electron momenta
leads to the symmetry relation
\begin{align}
a_{22}^{(1)} = a_{11}^{(1)} (y_1 \leftrightarrow y_2)~.
\label{p2p2}
\end{align}

% gmunu ---------------------------------------------------------------
The large-angle contributions have been calculated 
in Ref.~\cite{Rodrigo:2001jr}.
The coefficient proportional to $g^{\mu \nu}$ reads
\begin{align}
a_{00}^{(1,0)} &= \frac{1}{y_1\, y_2} \bigg[
-\frac{q^2(1-q^2)}{2}  - y_1 y_2
 - \bigg[ q^2 + \frac{2y_1 y_2}{1-q^2} \bigg] \log(q^2) \non \\
&+ \bigg\{ \frac{y_1}{2} \bigg[ 4-y_1-\frac{3(1+q^2)}{1-y_2} \bigg] 
\log(\frac{y_1}{q^2}) \non \\
& - \bigg[ 1 + (1-y_2)^2 + \frac{y_1 q^2}{y_2} \bigg] L(y_1) + 
(y_1 \leftrightarrow y_2)\bigg\}
\bigg]~,
\label{gmunu}
\end{align}
% p1 p1 ---------------------------------------------------------------
where the function $L$ is defined as
\begin{align}
L(y_i) &= \Li_2(-\frac{y_i}{q^2})-\Li_2(1-\frac{1}{q^2}) \non \\ & + 
\log(q^2+y_i) \log(\frac{y_i}{q^2})~, 
\end{align}
with $\Li_2$ the Spence (or dilogarithmic) function defined below
Eq.~(\ref{js}). 
The coefficient in front of the tensor structure 
$p_1^{\mu}p_1^{\nu}$ is given by 
\begin{align}
\label{p1p1}
a_{11}^{(1,0)} &= 
\frac{1}{y_1\, y_2} \bigg[
(1+q^2)^2 \left(\frac{1}{1-y_1}-\frac{1}{1-q^2} \right) -
\frac{4(1-y_2)y_1}{1-q^2} \non \\
&- \frac{2 q^2}{1-q^2} \bigg[(1-y_2)\left(\frac{1}{y_2}+\frac{q^2}{y_1}
+ \frac{2 y_1}{1-q^2} \right) \non \\ & +\frac{2q^2}{1-q^2}
\bigg] \log(q^2) 
- q^2 \bigg[ 1 + \frac{2}{y_2} \bigg] \log(\frac{y_1}{q^2}) \non \\ 
&- q^2 \bigg[ \frac{(2-3y_1)(1-y_2)^2}{y_1(1-y_1)^2} \bigg] 
 \log(\frac{y_2}{q^2}) \non \\
& - 2 q^2 \bigg[ 1+\frac{1}{y_2^2} \bigg] L(y_1)
- 2 q^2 \bigg[ 3+\frac{2q^2}{y_1}+\frac{q^4}{y_1^2} \bigg] L(y_2)
\bigg]~.
\end{align}
% p1 p2 ---------------------------------------------------------------
For the symmetric tensor structure 
($p_1^{\mu}p_2^{\nu}+p_2^{\mu}p_1^{\nu}$) one gets 
\begin{align}
\label{p1p2}
a_{12}^{(1,0)} &= \frac{1}{y_1\, y_2} \bigg[
- \frac{4q^2+(y_1-y_2)^2}{1-q^2} \non \\
&-2 q^2 \bigg[ \frac{q^2}{y_1 y_2}+\frac{1+q^2-2 y_1 y_2}{(1-q^2)^2} \bigg] 
\log(q^2) + \bigg\{ \frac{q^2}{1-y_1} \non \\
& - \frac{2 q^2}{1-y_2} \bigg[1-y_1+\frac{q^2}{y_2}
-\frac{q^2}{2(1-y_2)} \bigg] \log(\frac{y_1}{q^2}) \non \\  
& - 2 q^2 \bigg[ 1+\frac{q^2}{y_2}+\frac{q^2}{y_2^2} \bigg] L(y_1) 
+ (y_1 \leftrightarrow y_2) \bigg\}
\bigg]~. 
\end{align}
% asymmetry ---------------------------------------------------------------
Finally, the antisymmetric coefficient $a_{-1}$ 
accompanying ($p_1^{\mu}p_2^{\nu}-p_2^{\mu}p_1^{\nu}$) reads 
\begin{align}
a_{-1}^{(1,0)} &= \frac{q^2}{y_1\, y_2} 
\bigg[ \frac{2\log(1-y_1)}{y_1} 
+ \frac{1-q^2}{1-y_1} + \frac{q^2}{(1-y_1)^2} \bigg] \non \\ &
 - (y_1\leftrightarrow y_2)~.
\label{asym}
\end{align}


The mass-suppressed coefficients $a_{ij}^{(1,m)}$ 
are given by~\cite{Kuhn:2002xg}
\begin{align}
a_{00}^{(1,m)} &=  
\frac{m^2 q^2}{y_1^2} \bigg[
\log(q^2) \log(\frac{y_1^4}{m^4 q^2}) + 4 \Li_2(1-q^2) \non \\ &
+ \Li_2(1-\frac{y_1}{m^2}) - \frac{\pi^2}{6}  \bigg] 
- \frac{m^2(1-q^2)}{y_1^2} \bigg[ 1 - \log(\frac{y_1}{m^2}) \non \\ &
+ \frac{m^2}{y_1} \bigg( \Li_2(1-\frac{y_1}{m^2})-\frac{\pi^2}{6} \bigg) 
\bigg] + \frac{q^2}{2} n(y_1,\frac{1-3q^2}{q^2}) \non \\ &
+ (y_1 \leftrightarrow y_2)~,
\label{apiupiu}
\end{align}
whereas 
\begin{align}
a_{11}^{(1,m)} &= \frac{q^2}{1-q^2} \bigg\{
\frac{4m^2}{y_1^2} \bigg[1-\log(\frac{y_1}{m^2}) \non \\ &
+ \frac{m^2}{y_1} \bigg( \Li_2(1-\frac{y_1}{m^2})-\frac{\pi^2}{6} \bigg) \bigg]
- n(y_1,1) \non \\ &
+\frac{2m^2 q^2}{y_1(m^2(1-q^2)-y_1)} \bigg[
\frac{1}{q^2} \log(\frac{y_1}{m^2}) + \frac{\log(q^2)}{1-q^2} \non \\ &
+ \bigg(1+\frac{m^2}{m^2(1-q^2)-y_1}\bigg) N(y_1) \bigg] \bigg\} + \non \\ &
+ \frac{1}{1-q^2} \bigg\{
 \frac{4m^2(1-q^2)}{y_2^2} \bigg[ \log(q^2) \log(\frac{y_2^4}{m^4 q^2}) 
\non \\ & + 4 \Li_2(1-q^2) 
+ 2 \bigg( \Li_2(1-\frac{y_2}{m^2})-\frac{\pi^2}{6} \bigg) \bigg] \non \\ & 
+ \frac{4m^2 q^2}{y_2^2} \bigg[1-\log(\frac{y_2}{m^2}) 
+ \bigg(1+\frac{m^2}{y_2}\bigg)
 \bigg( \Li_2(1-\frac{y_2}{m^2}) \non \\ & 
 -\frac{\pi^2}{6} \bigg) \bigg] 
- \frac{1-2q^4}{q^2} n(y_2,\frac{3-8q^2+6q^4}{1-2q^4}) \non \\ &
+\frac{2m^2}{y_2(m^2(1-q^2)-y_2)} \bigg[
\frac{1}{q^2} \log(\frac{y_2}{m^2}) + \frac{\log(q^2)}{1-q^2}  \non \\ &
+ \bigg(3+\frac{m^2}{m^2(1-q^2)-y_2}\bigg) N(y_2) \bigg] \bigg\}~,
\end{align}
and
\begin{align}
a_{12}^{(1,m)} &= \frac{q^2}{1-q^2} \bigg\{
 \frac{4m^2}{y_1^2} \bigg[1-\log(\frac{y_1}{m^2}) \non \\ &
 + \bigg( \frac{1}{2} + \frac{m^2}{y_1} \bigg)
 \bigg( \Li_2(1-\frac{y_1}{m^2})-\frac{\pi^2}{6} \bigg) \bigg] \non \\ &
- \frac{1-q^2}{q^2} n(y_1,\frac{1}{1-q^2})
+ \frac{2m^2}{y_1(m^2(1-q^2)-y_1)} \non \\ & \times \bigg[
\frac{1}{q^2}\log(\frac{y_1}{m^2}) + \frac{\log(q^2)}{1-q^2} \non \\ &
+ \bigg(2+\frac{m^2}{m^2(1-q^2)-y_1}\bigg) N(y_1)  \bigg] \bigg\} 
 + (y_1 \leftrightarrow y_2)~.
\end{align}
The asymmetric coefficient does not get mass corrections,
\begin{align}
a_{-1}^{(1,m)} &= 0~.
\end{align}
The functions $n(y_i,z)$ and $N(y_i)$ are defined through
\begin{align}
n(y_i,z) & = \frac{m^2}{y_i(m^2-y_i)} \bigg[ 1 + z \; \log(\frac{y_i}{m^2}) 
\bigg] \non \\ &
+ \frac{m^2}{(m^2-y_i)^2} \log(\frac{y_i}{m^2})~,
\end{align}
and
\begin{align}
N(y_i) &= \log(q^2) \log(\frac{y_i}{m^2})+\Li_2(1-q^2) \non \\ & +
\Li_2(1-\frac{y_i}{m^2})-\frac{\pi^2}{6}~.
\end{align}
The apparent singularity of the function 
$n(y_i,z)$ inside the phase space limits is 
compensated by the zero in the numerator. 
In the region $y_i$ close to $m^2$ it behaves as 
\begin{align}
n(y_i,z) \big|_{y_i \to m^2}&=
\frac{1}{y_i} \bigg[ 1 + z \; \log(\frac{y_i}{m^2}) \bigg] \non \\ &
- \frac{1}{m^2} \sum_{n=0} \bigg(\frac{1}{n+2}+\frac{z}{n+1}\bigg)
\left( 1-\frac{y_i}{m^2}\right)^n~.
\end{align}
Similarly, the function $N(y_i)$ guarantees that the coefficients 
$a^{(1)}_{ij}$ are finite in the limit $y_i \rightarrow m^2 (1-q^2)$:
\begin{equation}
\frac{m^2 N(y_i)}{m^2(1-q^2)-y_i} \bigg|_{y_i \to m^2(1-q^2)}  = 
-\frac{\log(1-q^2)}{q^2}-\frac{\log(q^2)}{1-q^2}~.
\end{equation}




{$\\$}
{\noindent \it Virtual and soft corrections to IFS\\}



The virtual plus soft photon corrections of 
the initial-state and final-state 
vertex (see \Eq{sigmaIFS}) to FSR and ISR, respectively,  
can be written as~\cite{Berends:1973tz,Berends:1987ab}
\bea
{\rm d}\sigma^{\rm V+S}_{\rm IFS} &=& \frac{\alpha}{\pi} \, 
\left[ \delta^{\rm V+S}(w) \, 
{\rm d}\sigma^{(0)}_{\rm FSR}(s) \right. \nn \\ && + \left. 
\eta^\mathrm{V+S}(s',w) \, 
{\rm d}\sigma^{(0)}_{\mathrm{ISR}}(s') \right]~,
\label{sigvspi}
\eea
where ${\rm d}\sigma^{(0)}_{\mathrm{FSR}}$ and ${\rm d}\sigma^{(0)}_{\mathrm{ISR}}$ 
are the leading order FSR and ISR differential cross sections, 
respectively, $w=E_\gamma^{\rm cut}/\sqrt{s}$ with $E_\gamma^{\rm cut}$
the maximal energy of the soft photon in the $e^+e^-$ 
c.m. rest frame, and $s'$ corresponds to the squared mass of the 
$h \bar h \gamma$ system. The function $\delta^{\mathrm{V+S}}(w)$
is independent of the final state. In the limit $m_e^2\ll s$,
\begin{align}
\delta^{\mathrm{V+S}}(w) = 2 \left[
(L-1)\log{(2w)}+\frac{3}{4}L-1+\frac{\pi^2}{6}\right]~,
\label{delvs}
\end{align}
where $L=\log(s/m_e^2)$.
For two pions in the final state, 
the function $\eta^{\rm V+S}(s',w)$ is given by 
\begin{align}
& \eta^{\rm V+S}(s',w)  = 
- 2 \left[ \frac{1+\beta_\pi^2}
{2\beta_\pi} \log(t_\pi) + 1 \right] \nn \\ & \times
\left[ \log(2w) + 1 + \frac{s'}{s'-s} 
\log \left(\frac{s}{s'} \right) \right]
+ \log \left(\frac{m_\pi^2}{s'}\right) \nn \\ & 
- \frac{1+\beta_\pi^2}{\beta_\pi} \left[  
2 \Li_2 (1-t_\pi) + \log(t_{\pi}) \log (1+t_\pi) 
- \frac{\pi^2}{2} \right] \nn \\ &
- \frac{2+\beta_\pi^2}{\beta_\pi} \log(t_\pi) - 2~, 
\label{etavspi}
\end{align}
where 
\bea
\beta_\pi = \sqrt{1 - \frac{4m_\pi^2}{s'}}~, 
\qquad t_\pi = \frac{1-\beta_\pi}{1+\beta_\pi}~.
\label{betapi}
\eea
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

The function $\eta^\mathrm{V+S}(s',w)$ is equivalent to the familiar
correction factor derived in~\cite{Schwinger:1989ix,Drees:1990te} for the reaction 
$e^+e^- \to \pi^+ \pi^- \gamma$ in the framework of sQED 
(see also \cite{Hoefer:2001mx}) in the limit $s \to s'$:
\beq
\left. \log(2w) + 1 + \frac{s'}{s'-s} 
\log \left(\frac{s}{s'} \right) \right|_{s \to s'} = \log(2w')
\label{changeframe}
\eeq
with $w' = E_\gamma^{\rm cut}/\sqrt{s'}$.
The factor on the right hand side of \Eq{changeframe}
for $s\ne s'$ arises from defining the soft photon cutoff 
in the $e^+e^-$ laboratory frame. 

Correspondingly, the function $\eta^{\rm V+S}(s',w)$
for two muons in the final state reads 
\begin{align}
&\eta^{\rm V+S}(s',w) 
= -2 \left[\frac{1+\beta_{\mu}^2}
{2\beta_{\mu}} \log(t_\mu) + 1 \right] \nn \\ & \times
\left[ \log(2w) + 1 + \frac{s'}{s'-s} 
\log \left(\frac{s}{s'} \right) \right]
+ \log \left(\frac{m_\mu^2}{s'}\right)  \nn \\ &
-\frac{1+\beta_{\mu}^2}{2\beta_{\mu}} \left[
4 \Li_2 ( 1-t_\mu) - 2 \log(t_\mu) \log\left(\frac{1+\beta_{\mu}}{2}\right)
- \pi^2 \right] \nn \\ & 
- \frac{1}{\beta_\mu} 
\left[\frac{3}{3-\beta_\mu^2} + \beta_\mu^2 \right] \log(t_\mu) - 2~, 
\label{etavsmu}
\end{align}
where
\bea
\beta_\mu = \sqrt{1 - \frac{4m_\mu^2}{s'}}~, 
\qquad t_\mu = \frac{1-\beta_{\mu}}{1+\beta_{\mu}}~. 
\label{betamu}
\eea



{$\\$}
{\noindent \it Real corrections \\}


Matrix elements for the emission of two real photons, 
\begin{align}
e^+(p_1) + e^-(p_2) \rightarrow  {\rm hadrons}\, (Q) 
+ \gamma(k_1) + \gamma(k_2)~,
\end{align}
are calculated in PHOKHARA following 
the helicity amplitude method with the conventions 
introduced in~\cite{Jegerlehner:1999wu,Kolodziej:1991pk}.
The Weyl representation for fermions is used where the Dirac matrices 
\begin{align}
\gamma^{\mu}= \left( 
\begin{array}{cc}
0 & \sigma^{\mu}_+ \\
\sigma^{\mu}_- & 0
\end{array} \right)~, \qquad \mu=0,1,2,3~,
\end{align}
are given in terms of the unit $2 \times 2$ matrix $I$ and the Pauli 
matrices $\sigma_i,\ i=1,2,3$, with $\sigma^{\mu}_{\pm}=(I,\pm\sigma_i)$.
The contraction of any four-vector $a^{\mu}$ with the $\gamma^{\mu}$ 
matrices has the form
\begin{align}
\ta{a} = a_{\mu} \gamma^{\mu} = \left( 
\begin{array}{cc}
0   & a^+ \\
a^- & 0
\end{array} \right)~,
\end{align}
where the $2 \times 2$ matrices $a^{\pm}$ are given by 
\begin{align}
a^{\pm} =  a^{\mu} \sigma_{\mu}^{\pm} &= \left( \begin{array}{cc}
a^0 \mp a^3      & \mp (a^1 - i a^2) \\
\mp (a^1 + i a^2) & a^0 \pm a^3
\end{array} \right)~.
\end{align}

The helicity spinors $u$ and $v$ for a particle and an antiparticle 
of four-momentum $p=(E, {\bf p})$ and helicity $\lambda=\pm 1/2$ 
are given by
\begin{align}
u(p,\lambda=\pm 1/2) &= \left( \begin{array}{c}
\sqrt{E \mp |{\bf p}|} \; \chi({\bf p},\pm) \\
\sqrt{E \pm |{\bf p}|} \; \chi({\bf p},\pm)
\end{array} \right) \equiv \left( \begin{array}{c}
u_I \\ u_{II} \end{array} \right)~, \non \\ 
v(p,\lambda=\pm 1/2) &= \left( \begin{array}{c}
\mp \sqrt{E \pm |{\bf p}|}\; \chi({\bf p},\mp) \\
\pm \sqrt{E \mp |{\bf p}|}\; \chi({\bf p},\mp)
\end{array} \right) \equiv \left( \begin{array}{c}
v_I \\ v_{II} \end{array} \right)~.
\end{align}
The helicity eigenstates $\chi({\bf p},\lambda)$ can be expressed
in terms of the polar and azimuthal angles of the momentum
vector ${\bf p}$ as
\begin{align}
\chi({\bf p},+) &= \left( \begin{array}{r}
\cos{(\theta/2)} \\ e^{i\phi} \sin{(\theta/2)}
\end{array} \right)~, \non \\ 
\chi({\bf p},-) &= \left( \begin{array}{r}
-e^{-i\phi}\sin{(\theta/2)} \\ \cos{(\theta/2)}
\end{array} \right)~.
\end{align}
Finally, complex polarisation vectors in the helicity basis are defined
for the real photons: 
\begin{eqnarray}
\varepsilon^{\mu}(k_i,\lambda_i=\pm) &=& \frac{1}{\sqrt{2}} \big(  0,
\mp \cos \theta_i \cos \phi_i + i \sin \phi_i, \non \\  
&\mp& \cos \theta_i \sin \phi_i - i \cos \phi_i,
\pm \sin{\theta_i} \big)\,,
\end{eqnarray}  
with $i=1,2$.





{$\\$}
{\noindent \it Phase space \\}
%\subsubsection{Phase space}
%\label{app:phasespace}

One of the key ingredients of any Monte Carlo simulation is an 
efficient generation of the phase space. 
The generation of the multi-particle phase space in PHOKHARA
is based on the Lorentz-invariant representation 
\begin{align} 
& {\rm d} \Phi_{m+n}(p_1,p_2;k_1,\cdot,k_m,q_1,\cdot,q_n)
= \non \\ & \qquad {\rm d} \Phi_m(p_1,p_2;Q,k_1,\cdot,k_m) {\rm d} \Phi_n(Q;q_1,\cdot,q_n) 
\frac{{\rm d}Q^2}{2\pi}~,
\end{align}
where $p_1$ and $p_2$ are the four-momenta of the initial particles,
$k_1 \ldots k_m$ are the four momenta of the emitted photons and 
$q_1 \ldots q_n$, with $Q = \sum q_i$, label the four-momenta of the 
final state hadrons. 

When two particles of the same mass are produced in the final state, 
$q_i^2=M^2$, their phase space is given by
\begin{equation}
{\rm d} \Phi_2(Q;q_1,q_2) = \frac{\sqrt{1-\frac{4 M^2}{Q^2}}}{32 \pi^2} 
{\rm d} \Omega~,
\end{equation}
where ${\rm d} \Omega$ is the solid angle of one of the final state 
particles at, for instance, the $Q^2$ rest frame.
 
Single photon emission is described by the corresponding 
leptonic part of the phase space, 
\begin{equation}
{\rm d} \Phi_2(p_1,p_2;Q,k_1) = \frac{1-q^2}{32 \pi^2} {\rm d} \Omega_1~,
\end{equation}
with $q^2=Q^2/s$ and $d \Omega_1$ the solid angle of the emitted 
photon at the $e^+ e^-$ rest frame. The polar angle $\theta_1$
is defined with respect to the positron momentum $p_1$.
In order to make the Monte Carlo generation more efficient,
the following substitution is performed:
\begin{align}
\cos \theta_1 = \frac{1}{\beta} \tanh(\beta \; t_1) ~, \quad
t_1 = \frac{1}{2\beta} \log \frac{1+\beta \cos \theta_1}{1-\beta \cos \theta_1}~,
\label{t1}
\end{align}
with $\beta = \sqrt{1-4m_e^2/s}$,
which accounts for the collinear emission peaks  
\begin{equation}
\frac{{\rm d} \cos \theta_1}{1-\beta^2 \cos^2 \theta_1} = {\rm d}t_1~.
\end{equation}
With this the azimuthal angle and the new variable $t_1$ are generated flat. 

Considering the emission of two real photons
in the c.m. of the initial particles, the four-momenta 
of the po\-si\-tron, the electron and the two emitted photons are given by 
\begin{align}
p_1 &= \frac{\sqrt{s}}{2}(1,0,0,\beta)~, \qquad 
p_2 = \frac{\sqrt{s}}{2}(1,0,0,-\beta)~, \non \\
k_1 &= w_1 \sqrt{s} (1,\sin \theta_1 \cos \phi_1,\sin \theta_1 \sin \phi_1,
\cos \theta_1)~, \non \\
k_2 &= w_2 \sqrt{s} (1,\sin \theta_2 \cos \phi_2,\sin \theta_2 \sin \phi_2,
\cos \theta_2)~,
\end{align}
respectively. The polar angles $\theta_1$ and $\theta_2$ are
again defined with respect to the positron momentum $p_1$.
Both photons are generated with energies larger than 
the soft photon cutoff: $w_i>w$ with $i=1,2$. At least 
one of these exceeds the minimal detection energy: 
$w_1 > E_{\gamma}^{\rm min}/\sqrt{s}$ or 
$w_2 > E_{\gamma}^{\rm min}/\sqrt{s}$.
In terms of the solid angles ${\rm d} \Omega_1$ and ${\rm d} \Omega_2$ of the 
two photons and the normalised energy of one of them, e.g. $w_1$,
the leptonic part of the phase space reads 
\begin{align}
{\rm d} \Phi_3(p_1,p_2;Q,k_1,k_2) &= \frac{1}{2!} \;\frac{s}{4(2\pi)^5} 
\non \\ \times &  \frac{w_1 w_2^2}{1-q^2-2w_1}
\; {\rm d}w_1 \; {\rm d}\Omega_1 \; {\rm d}\Omega_2~,
\end{align}
where the limits of the phase space are determined from the constraint 
\begin{equation}
q^2 = 1-2(w_1+w_2)+2w_1 w_2(1-\cos \chi_{12})~,
\label{dpslimits}
\end{equation}
with $\chi_{12}$ being the angle between the two photons
\begin{equation}
\cos \chi_{12} = \sin \theta_1 \sin \theta_2 \cos (\phi_1 - \phi_2) 
+ \cos \theta_1 \cos \theta_2~.
\end{equation}

Again, the matrix element squared contains several peaks, soft and 
collinear, which should be softened by choosing suitable substitutions
in order to achieve an efficient Monte Carlo generator.
The leading behaviour of the matrix element squared is given 
by $1/(y_{11} \; y_{12} \; y_{21} \; y_{22})$, where  
\begin{equation}
y_{ij} = \frac{2 k_i \cdot p_j}{s} = w_i (1 \mp \beta \cos \theta_i)~.
\end{equation}
In combination with the leptonic part of the phase space, we have
\begin{align}
& \frac{{\rm d}\Phi_3(p_1,p_2;Q,k_1,k_2)}{y_{11} \; y_{12} \; y_{21} \; y_{22}}
\sim
\frac{{\rm d}w_1}{w_1(1-q^2-2w_1)} \; \non \\ & \qquad \times 
\frac{{\rm d}\Omega_1}{1-\beta^2 \cos^2{\theta_1}}
\; \frac{{\rm d}\Omega_2}{1-\beta^2 \cos^2{\theta_2}}~.
\end{align}
The collinear peaks are then flattened with the help of \Eq{t1}, with 
one change of variables for each photon polar angle. The remaining soft 
peak, $w_1 \rightarrow w$, is reabsorbed with the following substitution
\begin{align}
w_1 = \frac{1-q^2}{2+e^{-u_1}}~, \quad u_1 = \log \frac{w_1}{1-q^2-2w_1}~,
\end{align}
or
\begin{align}
\frac{{\rm d} w_1}{w_1(1-q^2-2w_1)} = \frac{{\rm d} u_1}{1-q^2}~,
\end{align}
where the new variable $u_1$ is generated flat. 
Multi-channe\-ling is used to absorb simultaneously the soft and 
collinear peaks, and the peaks of the form factors. 




{$\\$}
{\noindent \it NLO cross section and theoretical uncertainty \\}
%\subsubsection{NLO cross-section}
%\label{nlocrosssection}

The LO and NLO predictions for the differential cross section of the 
process $e^+ e^- \rightarrow \pi^+ \pi^- \gamma (\gamma)$ at DA$\mathrm{\Phi}$NE 
energies, $\sqrt{s} = 1.02$~GeV, are presented in 
Fig.~\ref{fig:nlo1gev} as a function of 
the invariant mass of the hadronic system $M_{\pi\pi}$.
We choose the same kinematical cuts as in the small angle analysis 
of KLOE~\cite{:2008en}; pions are restricted to be in the central region, 
$50^o < \theta_{\pi}< 130^o$, with $|p_T|>160$~MeV or $|p_z|>90$~MeV, 
the hard photon is not tagged and the sum of the momenta of the 
two pions, which flows in the opposite direction to the photon's momenta,  
is close to the beam ($\theta_{\pi\pi}<15^o$ or $\theta_{\pi\pi}>165^o$). 
The track mass, which is calculated from the equation 
\bea
&& \!\!\!\!\!\!\!\!\!\!
\left(\sqrt{s}-\sqrt{|\vec{p}_{\pi^+}|^2+M_{\rm trk}^2} 
- \sqrt{|\vec{p}_{\pi^-}|^2+M_{\rm trk}^2} \right)^2 \nn \\ &&
- (\vec{p}_{\pi^+}+\vec{p}_{\pi^-})^2 = 0~,
\label{eq:mtrkdef}
\eea
lies within the limits $130$~MeV$<M_{\rm trk}<220$~MeV 
and $M_{\rm trk}< (250 - 105 \sqrt{1-(M_{\pi\pi}^2/0.85)^2})$~MeV,
with $M_{\pi\pi}$ in GeV,
in order to reject $\mu^+\mu^-$ and $\pi^+\pi^-\pi^0$ events. 
The cut on the track mass, however, does not have any effect for single 
photon emission, as obviously $M_{\rm trk}=m_\pi$ for such 
events. 

The lower plot in Fig.~\ref{fig:nlo1gev} shows the relative 
size, with respect to the LO prediction, of FSR at LO, ISR 
corrections at NLO, and IFS contributions. The NLO ISR radiative 
corrections are almost flat and of the order of $-8$\%, FSR is 
clearly below $1$\%, while IFS corrections are also small 
although they become of the order of a few per cent at high 
values of $M_{\pi\pi}$.

\begin{figure}
\begin{center}
\epsfig{file=pi_track_1GeV.eps,width=9cm} 
\end{center}
\caption{Differential cross section for the process 
$e^+ e^- \rightarrow \pi^+ \pi^- \gamma$ at LO and NLO for 
$\sqrt{s}=1.02$~GeV. The cuts are the same as in the 
small angle analysis of KLOE, including the cut on the 
track mass. The lower plot shows the relative size of FSR 
at LO, ISR at NLO and IFS contributions with respect to 
the full LO prediction.}
\label{fig:nlo1gev}
\end{figure}

To estimate the systematic uncertainty of the NLO prediction, we observe that
leading logarithmic two-loop $O (\alpha^2)$ corrections and the 
associate real emission are not included. 
For samples with untagged photons the process 
$e^+e^- \to e^+e^-\pi^+\pi^-$ might also become a sizable background. 
This process, however, can be simulated with the Monte Carlo event generator
EKHARA~\cite{Czyz:2006dm,Czyz:2005ab}. Its contribution depends on the pion
pair invariant mass, ranges 
from $0.1-0.8$\% for the KLOE event selection,
 and has been taken into account in the 
KLOE analysis~\cite{:2008en}.

From na\"{\i}ve exponentiation one expects that LL corrections 
at next-to-next-to-leading order (NNLO) are of 
the order of $\frac{1}{2} ( \frac{3}{2} (\alpha/\pi) \log(s/m_e^2) )^2 
\approx 0.1$--$0.2\%$ for inclusive observables. 
For less inclusive distributions, a larger error is expected. 
The conservative estimate of the accuracy of PHOKHARA 
from ISR is $0.5$\%. 
This has been confirmed by comparisons with 
KKMC~\cite{Jadach:2006fx,Glosser:2003ux}, where
the biggest observed difference is about
$0.3$\% in the invariant mass regions which are not
close to the nominal energies of the experiments.
Improving the accuracy of PHOKHARA below $0.5$\%,
however, will be required to meet the growing experimental 
requirements in the near future. 


\subsubsection{FSR beyond sQED$*$VMD model}
\label{sQEDVMD}

  The model for FSR from pions described in details in
 Sections \ref{sec:LO} and  \ref{rratnlo} will be called for short
 the sQED$*$VMD model.
  The question arises how well
it
  can reflect the data. 
As shown in \cite{Baier:1965jz}, the first two terms in the
 expansion of the FSR amplitude as a function of $k^0/\sqrt{ Q^2}$
(i.e. the divergence and the constant) are fully given by the pion
 form factor. Thus one could expect that going
  beyond this approximation is necessary only for a hard photon emission. Moreover,
   the pion form factor is extremely big in the $\rho$ resonance
  region, and thus the validity of this approximation is further extended.
  In the kinematical regions where resonance contributions are 
  not contained in  the pion form factor, and also
   near the $\pi^+\pi^-$ threshold, where the emitted photon is hard
  and the pion form factor is relatively small, it is necessary to go
  beyond the sQED$*$VMD model and one needs a more general description
  of the  amplitude  $ M(\gamma ^{\ast
}( Q)\rightarrow \gamma (k)+\pi ^{+}(q_1)+\pi ^{-}(q_2))$.

In the general case the amplitude of the reaction

\noindent
 $\gamma ^{\ast
}( Q)\rightarrow \gamma (k)+\pi ^{+}(q_1)+\pi ^{-}(q_2)$
depends on  three 4-momenta, which can be chosen as $Q$, $k$ and
$l\equiv q_1-q_2$. The second-rank Lorentz
tensor $M^{\mu \nu }( Q,k,l)$ that describes the FSR amplitude 
can be decomposed through ten independent tensors
\cite{Tarrach:1975tu,Drechsel:1996ag}. Taking into account the charge
conjugation symmetry of the S-matrix element

 \noindent
($\langle \gamma
(k),\pi ^{+}(q_1)\pi ^{-}(q_2)|S|\gamma ^{\ast }(Q)\rangle
=$ 

\hskip 3 cm $\langle \gamma (k),\pi ^{-}(q_1)\pi ^{+}(q_2)|S|\gamma ^{\ast
}(Q)\rangle $),

 \noindent
 the photon crossing symmetry ($Q\leftrightarrow
-k$ and $\mu \leftrightarrow \nu $) and the gauge invariance
conditions $Q_{\mu }M^{\mu \nu }(Q,k,l)=0$ and $M_F^{\mu \nu
}(Q,k,l)k_{\nu }=0$, the number of independent tensors
decreases to five. For  a
final  real photon, i.e. $k^{2}=0$ and $k^{\nu }\epsilon _{\nu
}=0$ ($\epsilon _{\nu }$ being the polarisation
vector of the final photon) and  the initial virtual photon produced in
$e^{+}e^{-}$ annihilation ( $Q^{2}\geq
4m_{\pi}^{2}$), the FSR tensor can be rewritten in 
terms of  three gauge invariant tensors \cite{Tarrach:1975tu,Drechsel:1996ag}
\begin{equation}
M^{\mu \nu }(Q,k,l)=\tau _{1}^{\mu \nu }f_{1}+\tau _{2}^{\mu \nu
}f_{2}+\tau _{3}^{\mu \nu }f_{3}, \label{eq:A3}
\end{equation}
where the gauge invariant tensors $\tau_i^{\mu\nu}$ read
\begin{eqnarray}
\tau _{1}^{\mu \nu } &=&k^{\mu }Q^{\nu }-g^{\mu \nu }k\cdot Q,  \\
\tau _{2}^{\mu \nu } &=&k\cdot l(l^{\mu }Q^{\nu }-g^{\mu \nu
}k\cdot
l)+l^{\nu }(k^{\mu }k \cdot l-l^{\mu }k \cdot Q),  \nonumber \\
\tau _{3}^{\mu \nu } &=&Q^{2}(g^{\mu \nu }k\cdot l-k^{\mu }l^{\nu
})+Q^{\mu }(l^{\nu }k\cdot Q-Q^{\nu }k\cdot l).\nonumber  \label{eq:A4}
\end{eqnarray}

It thus follows that the evaluation of the FSR tensor amounts to
the calculation of the scalar functions

 $f_{i}(Q^2,Q\cdot k,k\cdot l) $ \ ($i=1,2,3$).



  As is clear from the above discussion, the extraction of 
  the pion form factor from radiative return experiments is a demanding
  task. The main problem is that in the same experiment one has
  to test the models describing the pion-photon interactions (see 
   Section \ref{expvsth}) and to extract the pion form factor needed
  for the evaluation of the muon anomalous magnetic moment. Fortunately,
  there are event selections, which naturally suppress the FSR contributions,
  independently of their nature. These were already discussed in Section
   \ref{sec:LO} in the context of the sQED$*$VMD model.

  Extensive theoretical studies of the role of the FSR emission
  beyond the sQED$*$VMD model were performed 
  \cite{Czyz:2004nq,Melnikov:2000gs,Pancheri:2006cp,Shekhovtsova:2009yn,Pancheri:2007xt}. They
  are important  mainly for the KLOE measurements at
  DA$\mathrm{\Phi}$NE, 
  as at $B$ factories FSR is naturally suppressed and the accuracy
 needed in its modelling is by far less demanding than that for KLOE purposes. 


For DA$\mathrm{\Phi}$NE, running on or near the $\mathrm{\phi}$ resonance,
the following mechanisms of the
$\pi^+\pi^-$ final state photon emission have to be considered:

 -- bremsstrahlung process
\begin{eqnarray}
e^++e^-&\to&\pi^++\pi^-+\gamma \ ,
  \label{brem}
\end{eqnarray}

\noindent
which is modelled by sQED$*$VMD;

-- direct $\mathrm{\phi}$ decay
\begin{eqnarray}
e^++e^-&\to&\phi\to (f_0;f_0+\sigma)\gamma\to\pi^++\pi^-+\gamma \ ,
\label{phi_direct} 
\end{eqnarray}

\noindent
and 

-- double resonance process
\begin{eqnarray}
e^++e^-&\to&(\phi;\omega')\to \rho\pi\to\pi^++\pi^-+\gamma \ .
\label{phi_vmd}
\end{eqnarray}

\noindent
The resonance chiral theory
(R$\chi$T)~\cite{Ecker:1989yg,Ecker:1988te}
  was used in \cite{Pancheri:2006cp,Pancheri:2007xt}
 to estimate the contributions
beyond sQED$*$VMD. They were implemented at leading order
into the event generator FASTERD \cite{Shekhovtsova:2009yn}.
Having in mind that at present these models still await accurate
 experimental tests, 
 other models \cite{Achasov:1996xz,Achasov:2005hm} were also
 implemented in the event generator FASTERD.
To include both next-to-leading-order radiative corrections
 and the mechanisms discussed 
for FSR, a part of the FASTERD code, based on 
the models \cite{Achasov:1996xz,Achasov:2005hm}, was implemented
by O.~Shekhovtsova in PHOKHARA v6.0 (PHOKHARA v6.1 \cite{OlgaS}) and the studies
presented below are based on this code. The model
used there, even if far from an ideal,
 is the best tested model available in literature.

 We briefly 
describe main features of the models 
 used to describe processes contributing to FSR photon emission
listed above.
 For a more detailed description and the calculation of the
function $f_i$ we refer the reader to 
 \cite{Czyz:2004nq,Melnikov:2000gs,Shekhovtsova:2009yn} (see also references therein).

The sQED$*$VMD part gives contributions to $f_1$  and $f_2$.

The direct $\mathrm{\phi}$ decay is assumed to proceed through the %$f_0$
intermediate scalar meson state:
 $\phi\to (f_0+\sigma)\gamma\to\pi\pi\gamma$.
 Various models  are proposed to
describe the $\phi$-scalar-$\gamma$ vertex: either
it is the direct decay $\phi\to \text(scalar)\gamma$, or the
vertex is generated dynamically through a loop of the charged
kaons.  As shown in~\cite{Melnikov:2000gs}, in the framework of any model, the
direct $\mathrm{\phi}$ decay affects only the form factor $f_1$ of
Eq.~({\ref{eq:A3}}). 

The double resonance contribution consists of the off-shell $\mathrm{\phi}$ meson
decay into ($\rho^\pm \pi^\mp$)  and subsequent decay
$\rho \to\pi\gamma$. In the energy region around $1$ GeV the
tail  of the excited $\omega$ meson  can  also play a role, and 
$\gamma^*\to\omega'\to\rho\pi$ has to be 
considered. The double resonance mechanism affects all three
form factors $f_i$ of Eq.~({\ref{eq:A3}}).

 Assuming isospin symmetry, this part can be deduced from the 
 measurement of the neutral pion pair production. Various models
 \cite{Achasov:1996xz,Achasov:2005hm}
 were confronted with data by KLOE \cite{Ambrosino:2006hb}
  for the neutral mode. The model that was 
 reproducing the data in the best way was adopted to be used for the
 charged pion pair production relying on the isospin symmetry \cite{OlgaS}.

   In \cite{Czyz:2004nq} it was shown that an important tool 
  for testing the various models of FSR is
  the charge asymmetry. At leading order it originates from the fact that 
  the pion pair couples to an even (odd) number of photons if the final state
  photon is emitted  from the final (initial) state. The interference diagrams
  do not give any contribution to the integrated cross section for C--even event
  selections, but produce an asymmetry in the angular distribution.
  The definitions and experimental studies based on the charge asymmetry
  are presented in Section \ref{radret:KLOEasymFSR}.
% The sensitivity of the charge forward backward asymmetry
%\begin{eqnarray}
% {\mathcal A}_{FB}(Q^2) = 
% \frac{N(\theta_{\pi^+}>90^\circ)- N(\theta_{\pi^+}<90^\circ) }
% {N(\theta_{\pi^+}>90^\circ)+ N(\theta_{\pi^+}<90^\circ)}\left(Q^2\right) \ ,
%\label{asymfb}
%\end{eqnarray}
%
%on various model assumptions is shown in Fig.\ref{f2}.
 
%\begin{figure}[ht]
% \vspace{0.5 cm}
%\includegraphics[width=9cm,height=8cm]{v_as_pipl_all.ps}
%\caption{Forward-backward asymmetry for different radiative $\phi$
% decay models compared with the asymmetry calculated within sQED (no $f_0$)
%}
% \vspace{0.5 cm}
%\label{f2}
%\end{figure}

  Few strategies can be adopted to profit in the best way
 from the KLOE data taken on and off peak.
  The 'easiest' part is to look for the event selections where the FSR
 contributions are negligible. This was performed by KLOE \cite{:2008en}
  (see Section \ref{rr:kloe}),
 giving important information on the pion form factor relevant for 
  the prediction of the hadronic contributions to the muon anomalous
   magnetic moment $a_\mu$. Typical contributions of the FSR (1 -- 4\%)
  to the differential cross section (Figs.~\ref{fig:nlo1gev} and \ref{f2})
  allow for excellent control of the accuracy of these corrections.
 One disadvantage of using this event selection is that it does not allow
  to perform measurements near the pion production threshold. 
  %
\begin{figure}[ht]
 \vspace{-0.6 cm}
\includegraphics[width=8cm,height=7cm]{fsr_nlo_sa_mphi_low.eps}
\caption{Relative contribution of the FSR to the differential cross section
 of the reaction $e^+e^-\to \pi^+\pi^-\gamma(\gamma)$ for $\sqrt{s}=m_\phi$
and low invariant masses of pion pairs.
 KLOE small angle event selection \cite{:2008en} was used, 
and for this event selection
 the relative contribution of the FSR is almost identical also
 for the off peak cross section. The effect of a trackmass cut
 (see Section  \ref{rr:kloe}) is  shown. ISRNLO refers to initial state
 corrections at next-to-leading order (NLO). The IFSNLO cross section contains
 the final state emissions at NLO. 
}
 \vspace{-0.5 cm}
\label{f2}
\end{figure}
%

The next step, partly discussed in Section \ref{radret:KLOEasymFSR},
   is to confront the models based on isospin symmetry and the
  neutral channel data with charged pion data taken off-peak, where the
  contributions from models beyond the sQED$*$VMD approximation 
  are relatively small
  (Fig.~\ref{fig3s}).
For the off-peak data \cite{Beltrame:2009zz}
 the region below $Q^2=0.3\ {\rm GeV}^2$ can be
   covered experimentally. However, the small statistics in this region
makes it difficult to perform high-pre\-ci\-sion tests of the models.
 For this analysis an accurate knowledge of the
   pion form factor at the nominal energy of the experiment is
   important, 
  as it defines the sQED$*$VMD predictions and as the FSR corrections 
  (Fig.~\ref{f:fsrla})  are sizeable.

  The last step, which  allows for the most accurate FSR model testing
   and profits from the knowledge of the pion form factor from previous
  analysis, is the on-peak large angle measurement. The large 
  FSR corrections coming from sources  beyond the sQED$*$VMD approximation
 (Figs.~\ref{fig3s} and \ref{f:fsrla})
  make these data \cite{Leone:2007zz} the most valuable source of information
  on these models. In this case, the accumulated data set is much larger than
  the off-peak data set and one is able to cover also the region
  below $Q^2=0.3\ {\rm GeV}^2$.
    
  
\begin{figure}[tbp]
\includegraphics[width=8cm,height=8cm]{fsr_sqed_nlo_la.eps}
\caption{The contributions of FSR beyond the sQED$*$VMD approximation (see 
Eqs.~(\ref{phi_direct}) and~(\ref{phi_vmd}))
 for KLOE large angle event selection \cite{Beltrame:2009zz,Leone:2007zz}
 for $\sqrt{s}=m_\phi$
 and for $\sqrt{s}= 1\ {\rm GeV}$. }
\label{fig3s}
\end{figure}

\begin{figure}[ht]
 \vspace{-0.65 cm}
\includegraphics[width=8cm,height=7cm]{fsr_tot_nlo_la.eps}
\caption{Relative contribution of FSR to the differential cross section
 of the reaction $e^+e^-\to \pi^+\pi^-\gamma(\gamma)$ for $\sqrt{s}=m_\phi$
 and for $\sqrt{s}= 1\ {\rm GeV}$.
 KLOE large angle event selection \cite{Beltrame:2009zz,Leone:2007zz} was used. 
}
 \vspace{-0.6 cm}
\label{f:fsrla}
\end{figure}

%\subsection{Monte Carlo tools}
%\label{radret:tools}


\subsection{Experiment confronting theory } 
\label{expvsth} 

\subsubsection{Study of the process $e^+e^- \to \pi^+\pi^-\gamma$ with FSR with
the CMD-2 detector at VEPP-2M}

%\input{pipigamma}
\input{eidelman_pipig}

 \subsubsection{Study of the process $e^+e^- \to \pi^+\pi^-\gamma$ with FSR with
 KLOE detector}
 \label{radret:KLOEasymFSR}

\input{KLOEasymFSR}
