%\newcommand{\ttbs}{\char'134}
%\newcommand{\AmS}{{\protect\the\textfont2
%  A\kern-.1667em\lower.5ex\hbox{M}\kern-.125emS}}
% add words to TeX's hyphenation exception list
\hyphenation{author another created financial paper re-commend-ed Post-Script}


%\section{Monte Carlo Generators for $\tau$ Decays}% (Experiment)}


\subsection{Introduction}


%Already soon 
%%%TT
After discovery of the $\tau$ lepton, which is a fundamental lepton,
heavy enough to decay not only into leptons, but also into
dozens of various hadronic final states, it became clear 
that corresponding Monte Carlo (MC) event generators are needed for various 
purposes:

\begin{itemize}
\item
{To calculate detector acceptance, efficiencies 
%$\epsilon$
 and various distributions for signal event selection and comparison
to data. In general the acceptance is small (a few percent) and depends on
the model; in principle, it is a complicated function of invariant
masses, angles, and resolutions. 
Analysis of publications shows that effects of MC signal modelling 
are almost always neglected.}
\item
{To estimate the number of background (BG) events
$N^{\rm BG}_{\rm ev}$ and their distributions;
%its shape;
in addition to background coming from $\tau^+\tau^-$ pairs (so called
cross-feed), there might be BG events from $q\bar{q}$ continuum,
$~\gamma\gamma$ collisions etc.}
\item
To unfold observed distributions 
to get rid of detector effects, 
%again $\epsilon$ is needed; 
important when extracting resonance parameters.
\end{itemize}
 
Various computer packages like, e.g.,  
KORALB~\cite{Jadach:1984iy}, KKMC~\cite{kkcpc:1999}, 
TAUOLA~\cite{Jadach:1990mz,Jezabek:1991qp,Jadach:1993hs} and 
PHOTOS~\cite{photos2:1994} 
were developed to generate events for $\tau$ lepton production
in $e^+e^-$ annihilation and their subsequent decay, taking into
account the possibility of photon emission.
These codes became very important tools 
for experiments at LEP, CLEO, Tevatron and HERA.

Simulation of hadronic decays requires the knowledge of hadronic form 
factors.  Various hadronic final states were
considered 
%by the group of J.~K\"{u}hn 
in the 90's, resulting in 
a large number of specific hadronic currents~\cite{Kuhn:1992nz}. 


%$V-A$ structure of $J_{\rm weak}$ plus phase space production
%in $J_{\rm hadr}$.

However, already experiments at LEP and CLEO show\-ed that with increase
of the collected data sets a more precise description is necessary.
Some attempts were made to improve the parametrisation of various
hadronic currents. One should note the serious efforts of the ALEPH
and CLEO Collaborations, which created their own parametrisations of TAUOLA
hadronic currents already in the late 90's, or 
a  parametrisation of the hadronic current in the $4\pi$ 
decays~\cite{Bondar:2002mw}, based on the experimental information on 
$e^+e^- \to 2\pi^+2\pi^-,~\pi^+\pi^-2\pi^0$ from Novosibirsk~\cite{Akhmetshin:1998df}, which is now implemented in the presently distributed TAUOLA code \cite{Golonka:2003xt}.
% These efforts find their place in the present day TAUOLA distribution \cite{Golonka:2003xt}.



\subsection{Current status of data and MC generators}
In this section we will briefly discuss the most precise recent
experimental data on $\tau$ lepton decays, showing, wherever possible, 
their comparison with the existing MC generators and discussing the 
decay dynamics.
 \begin{center}
\begin{figure}
%\includegraphics[width=0.37\textwidth]{bel21.eps}
\includegraphics[width=0.37\textwidth]{bel21.eps}
\includegraphics[width=0.37\textwidth]{bel22.eps}
\includegraphics[width=0.37\textwidth]{bel23.eps}
\caption{Projections to the missing mass
and missing direction for 
$\tau^- \to  \pi^-\pi^0\nu_{\tau}$ decays at Belle:
(a)--(c) correspond to different ranges of the missing polar angles.
The solid circles represent the data and the histograms the MC
simulation (signal + background). The open histogram shows 
the contribution from $\tau^+\tau^-$ pairs, the vertical
(horizontal) striped area shows that from two-photon leptonic
(hadronic) processes; the wide (narrow) hatched area 
shows that from Bhabha ($\mu^+\mu^-$), and the shaded area that from the 
$q\bar{q}$ continuum.} 
\label{fig:tau21}
\end{figure}
\end{center}
%\vspace*{0.5mm}
%{}
 
\begin{center}
\begin{figure}
%\begin{tabular}{cc}
%\includegraphics[width=0.50\textwidth]{bel24.eps}% &
%\\
%\includegraphics[width=0.45\textwidth]{bel25.eps}
\includegraphics[width=0.40\textwidth]{bel24.eps}% &
\\
\includegraphics[width=0.36\textwidth]{bel25.eps}
%\end{tabular}
\caption{Invariant-mass-squared distribution for
$\tau^- \to  \pi^-\pi^0\nu_{\tau}$ decay at Belle.
(a) Contributions of different background sources. The solid circles 
with error bars represent the data,
and the histogram represents the MC simulation 
(signal + background).
(b) Fully corrected distribution. The solid curve is the result
of a fit to the Gounaris-Sakurai model with the $\rho(770)$,
$\rho(1450)$ and $\rho(1700)$ resonances.}
\label{fig:tau22}
\end{figure}
\end{center}

\subsubsection{$\tau^- \to  \pi^-\pi^0\nu_{\tau}$ at Belle}
Recently results of a study of the
$\tau^- \to  \pi^-\pi^0\nu_{\tau}$ decay by the Belle Collaboration
were published~\cite{Fujikawa:2008ma}.
From less than 10\% of the dataset available the authors selected a
huge statistics of 5.4M events, about two orders of magnitude larger 
than in any previous experiment, determined the branching fraction and 
after the unfolding 
obtained the hadronic mass spectrum, in which for the first time
three $\rho$-like resonances were observed together: 
$\rho(770),~\rho(1450)$ and $\rho(1700)$. Their parameters were also
determined.

%\vspace*{0.5mm}
%{}


The comparison of the obtained missing mass distributions 
with simulations for different polar angle ranges
(Fig.~\ref{fig:tau21}) shows that there exist small discrepancies 
between MC and data. 




Figure~\ref{fig:tau22} shows various background contributions to
the di-pion mass distribution (upper panel) and underlying dynamics
(lower panel), clearly demonstrating a pattern of the three interfering
resonances $\rho(770),~\rho(1450)$ and $\rho(1700)$.



\subsubsection{$\tau^- \to
  \bar{K}^0\pi^-\nu_{\tau},~K^-\pi^0\nu_{\tau}$}

Two high-precision studies of the $\tau$ decay into the $K\pi\nu_{\tau}$ 
final state were recently published.  The BaBar Collaboration
reported  a measurement of the branching fraction of the
$\tau^- \to K^-\pi^0\nu_{\tau}$ decay~\cite{Aubert:2007jh}. They do not study
in detail the $K\pi$ invariant mass distribution, noting only that the
$K^*(892)^-$ resonance is seen prominently above the simulated
background, see Fig.~\ref{fig:tau23}. Near 1.4 GeV$/c^2$ decays 
to higher $K^*$ mesons are
expected, such as the $K^*(1410)^-$ and  $K^*_0(1430)^-$, but their
branching fractions are not yet measured well. These decays are not
included in the BaBar simulation of $\tau$ decays, but seem to be
present in the data around  1.4 GeV$/c^2$. It is also worth noting
that this decay mode is heavily contaminated by cross-feed backgrounds
from other $\tau$ decays. For example, below 0.7 GeV$/c^2$ the
background is dominated by $K^-\pi^0\pi^0\nu_{\tau}$ and  
$K^-K^0\pi^0\nu_{\tau}$ events,
for which the branching fractions are only known with large 
relative uncertainties of $\approx 37\%$ and $\approx 13\%$,
respectively. Non-negligible background may also come from
the $\tau^- \to \pi^-\pi^0\nu_{\tau}$ decay, which has a large
branching fraction and thus should be simulated properly.
    
\begin{center}
\begin{figure}
%\includegraphics[width=0.45\textwidth]{babkpi1.eps}
\includegraphics[width=0.35\textwidth]{babkpi1.eps}
%\vskip -9 mm
\caption{The $K\pi$ invariant mass distribution for the decay 
$\tau^- \to K^-\pi^0\nu_{\tau}$ at BaBar. The dots are the data,
while the histograms are background MC events with selection and 
efficiency corrections: $\tau$ background (dashed line), $q\bar{q}$
(dash-dotted line), $\mu^+\mu^-$ (dotted line).}
\label{fig:tau23}
\end{figure}
\end{center}
%
%
\begin{center}
\begin{figure}
%\begin{tabular}{cc}
\includegraphics[width=0.35\textwidth]{blkpi1.eps} %&
\\
\includegraphics[width=0.35\textwidth]{blkpi2.eps}
%\includegraphics[width=0.45\textwidth]{blkpi1.eps} %&
%\\
%\includegraphics[width=0.45\textwidth]{blkpi2.eps}
%\end{tabular}
%\vskip -9 mm
\caption{The $K\pi$ invariant mass distribution for the decay 
$\tau^- \to K^-\pi^0\nu_{\tau}$ at Belle. Points are experimental 
data, histograms are spectra expected for different models. 
(a) shows the fitted result in the model with the $K^*(892)$ alone.
(b) shows the fitted result in the  $K^*(892)+ K^*_0(800) + K^*(1410)$
model. Also shown are different types of background.} 
\label{fig:tau24}
\end{figure}
\end{center}
%\vspace*{0.5mm}
{}
Another charge combination of the final state particles, i.e.,
$K^0_S\pi^-\nu_{\tau}$, was studied in the Belle
experiment~\cite{:2007rf}. In this case a detailed analysis of the
$K\pi$ invariant mass distribution has been performed. The authors   
also conclude that the decay dynamics differs from pure $K^*(892)$: 
the best fit includes $K^*_0(800)+K^*(892)+K^*(1410)/K^*_0(1430)$,
see Fig.~\ref{fig:tau24}.
\begin{center}
\begin{figure}
\includegraphics[width=0.5\textwidth]{phik.eps}
%\includegraphics[width=0.40\textwidth]{phik.eps}
\caption{The $\phi K$ invariant mass distribution for the decay 
$\tau^- \to \phi K^- \nu_{\tau}$ at Belle. Points with error
bars are the data. The open histogram is the phase-space distributed 
signal MC, and dotted and dot-dashed histograms indicate the signal MC
mediated by a resonance with mass and width of 1650 MeV and 100 MeV,
and 1570 MeV and 150 MeV, respectively.}
\label{fig:tauphik}
\end{figure}
\end{center}

\begin{figure*}
%\begin{figure*}[h!]
\begin{center}
%\begin{figure}
\includegraphics[width=0.75\textwidth]{beleta1.ps}
%\includegraphics[width=0.85\textwidth]{beleta1.ps}
\caption{Invariant mass distributions: 
(a) $\pi\pi^0$ and (b) $\pi\eta\pi^0$ for 
$\tau \to \pi\pi^0\eta\nu_{\tau}$; (c) $\eta K$ for 
$\tau \to K\eta\nu_{\tau}$ and (d) $\pi K^0_S\eta$ for 
$\tau \to \pi K^0_S\eta\nu_{\tau}$ at Belle. The points with error 
bars are the data. The normal and filled histograms indicate the signal and
$\tau^+\tau^-$ background MC distributions, respectively.}
\label{fig:beleta}
%\end{figure}
\end{center}
\end{figure*}

%\clearpage

\subsubsection{$\tau$ decays into three pseudoscalars}

Recently a measurement of the branching fractions of
various particle combinations in the decay to three charged
hadrons (any combination of pions and kaons) was reported by the 
BaBar Collaboration~\cite{Aubert:2007mh}.
A similar study was also performed by the Belle group~\cite{:2008sg}.
However, both groups have not yet analysed the mass spectra in detail.
In the $K^-K^+K^-\nu_{\tau}$ final state BaBar~\cite{Aubert:2007mh} 
and Belle~\cite{Inami:2006vd} reported 
the observation of the decay mode $\phi K^-\nu_{\tau}$, while in the 
$K^-K^+\pi^-\nu_{\tau}$ final state BaBar observed the 
$\phi \pi^- \nu_{\tau}$ decay mode~\cite{Aubert:2007mh}. Belle analysed 
the spectrum of the $\phi K^-$ mass and concluded that it might 
have a complicated dynamics, see Fig.~\ref{fig:tauphik}.



The most detailed previous study of the mass spectra was done by the
CLEO group~\cite{Asner:1999kj}. With the statistics of about 8,000 events 
they conclude that the $3\pi$ mass spectrum is dominated by the
$a_1(1260)$ meson, and confirmed that the decay of the latter 
is not saturated by the $\rho\pi$ intermediate state, having in
addition a significant $f_0(600)\pi^-$ component observed earlier in
$e^+e^-$ annihilation into four charged pions~\cite{Akhmetshin:1998df}.

Recently the Belle Collaboration performed a detailed study of
various decays with the $\eta$ meson in the final
state~\cite{Inami:2008ar}. 
They measured
the branching fractions of the following decay modes:
$\tau^- \to K^-\eta\nu_{\tau}$, $\tau^- \to K^-\pi^0\eta\nu_{\tau}$, 
$\tau^- \to \pi^-\pi^0\eta\nu_{\tau}$, 
$\tau^- \to \pi^-K^0_S\eta\nu_{\tau}$, and 
$\tau^- \to K^{*-}\eta\nu_{\tau}$. They also set upper limits on the
  branching fractions of the  decays into $K^-K^0_S\eta\nu_{\tau}$,
$\pi^-K^0_S\pi^0\eta\nu_{\tau}$, $K^-\eta\eta\nu_{\tau}$,
 $\pi^-\eta\eta\nu_{\tau}$, and non-resonant  $K^-\pi^0\eta\nu_{\tau}$
final states.



Figure~\ref{fig:beleta} shows that there is reasonable agreement for
$\eta\pi^-\pi^0\nu_{\tau}$ (a, b) and a worse one for
$\eta K^- \nu_{\tau}$ (c) and $\eta K^{*-}\nu_{\tau}$ (d). 

\subsubsection{$\tau$ decays to four pions}

There are two possible isospin combinations of this hadronic final state,
$2\pi^-\pi^+\pi^0$ and $\pi^-3\pi^0$. Both have not yet been studied
at $B$ factories, so the best existing results are based on 
ALEPH~\cite{Buskulic:1996qs} and CLEO~\cite{Edwards:1999fj} results.
    
The theoretical description of such decays is based on the CVC relations 
and the available low energy $e^+e^-$ 
data~\cite{Decker:1994af,Czyz:2000wh,Bondar:2002mw,Czyz:2008kw}.

\subsubsection{$\tau^- \to  3h^-2h^+\nu_{\tau}$ at BaBar}

A new study of the $\tau^- \to  3h^-2h^+\nu_{\tau}$ decay ($h=\pi,~K$)
has been performed by the BaBar Collaboration~\cite{Aubert:2005wa}. A large
dataset of over 34,000 events (two orders of magnitude larger than in the
best previous measurement at CLEO~\cite{Gibaut:1994ik}) allows one 
a first search for resonant structures and decay dynamics. 

The invariant mass distribution of the five charged particles
in Fig.~\ref{fig:tau51} shows a clear discrepancy between the data 
and the MC simulation, which uses the phase space distribution
for $\tau^- \to  3\pi^-2\pi^+\nu_{\tau}$.




The mass of the $h^+h^-$ pair combinations in Fig.~\ref{fig:tau52}
(upper panel), 
with a prominent shoulder at 0.77 GeV$/c^2$, suggests a strong
contribution from the $\rho$ meson. Note that there are three allowed
isospin states for this decay, of which two may have a $\rho$ meson. 
The mass of the $2h^+2h^-$ combinations in Fig.~\ref{fig:tau52} (lower
panel) also shows a structure at 1.285 GeV$/c^2$ coming from the 
$\tau^- \to f_1(1285)\pi^-\nu_{\tau}$ decay.         

The first attempt to take into account the dynamics of this decay was
recently performed in Ref.~\cite{Kuhn:2006nw}. 



\subsubsection{$\tau$ decays to six pions}

The six-pion final state was studied by the CLEO 
Collaboration~\cite{Anastassov:2000xu}. Two charge combinations,
 $3\pi^-2\pi^+\pi^0$ and $2\pi^-\pi^+3\pi^0$, 
were observed and it was found that the decays are saturated by
intermediate states with $\eta$ and $\omega$ mesons. Despite the rather
limited statistics (about 260 events altogether), it became clear that
the dynamics of these decays is rather rich.


\subsubsection{Lepton-Flavour Violating Decays}

More than 50 different Lepton-Flavour Violating  (LFV)
decays have been studied by the CLEO, BaBar and Belle
Collaborations.
Publications rarely describe how the simulation of such decays is performed.
Moreover, theoretical papers suggesting LFV in new models usually 
do not provide differential cross sections. In some experimental
papers the authors claim that the production of final state
hadrons with a phase space distribution is assumed.
%production is assumed for final state
%hadrons. 
However, the real meaning of this statement is not very clear 
since LFV assumes New Physics and, therefore, matrix elements are not 
necessarily separated into weak and hadronic parts. 

However, there exist a few theoretical papers considering 
differential cross sections. For example,   
angular correlations for
$\tau^- \to \mu^-\gamma,~\mu^-\mu^+\mu^-$ and $\mu^-e^+e^-$ decays
were studied in Ref.~\cite{Kitano:2000fg}. An attempt to classify different
types of operators entering New Physics Lagrangians for $\tau$ decays 
to three charged leptons was made in~\cite{Dassinger:2007ru}.

\subsection{Status of Monte Carlo event generators for 
$\tau$ production and decays}
%Conclusions}

High-statistics and high-precision experiments, as well as  
searches for rare processes, result in a new challenge: 
Monte Carlo generators based on
an adequate theoretical description of energy and angular distributions.
%based on Monte Carlo generators. 
%In the next Sections
In the following we will 
describe the status of the Monte Carlo programs used by experiments. 
We will review the building blocks used in the simulation 
with the goal in mind to localise the
points requiring most urgent attention.

At present, for the production of $\tau$ pairs, the Monte Carlo programs
KORALB ~\cite{Jadach:1984iy} and
KKMC ~\cite{kkcpc:1999} are the standard
codes to be used. 
For the generation of brems\-strah\-lung in decays, 
the Monte Carlo PHOTOS~\cite{photos2:1994} is used. Finally, $\tau$ decays 
themselves are simulated 
with the program TAUOLA ~\cite{Jadach:1990mz,Jezabek:1991qp,Jadach:1993hs}. The 
EvtGen code was written and maintained for simulation of $B$ meson decays,
see \newline\texttt{www.slac.stanford.edu/\~{}lange/EvtGen/}\ .
It offers a \newline unique opportunity to specify, at run time, a list of the
final state particles\footnote{E.g. $\tau$ lepton decay products 
including neutrinos.},
%EvtGen was written and maintained for simulation of $B$ meson decays,
%see http://www.slac.stanford.edu/~lange/EvtGen/},
without having to change and/or compile the underlying code.  In
a multi-particle final state dominated by
phase space considerations, this generator provides an adequate
description of the final state momenta,
for which the underlying form factor calculation is more involved and
not presently available in a closed form. That is why it is used by 
experiments measuring $\tau$ decays too.
 

%{Spectral functions (shape) are important in various applications: \\
%$m_{\tau},~m_{\nu_{\tau}}$, CVC tests, $(g_{\mu}-2)/2$, QCD parameters
%($\alpha_s$, condensates, $m_s$)}

So far,  our discussion has been based on the comparison of experimental data 
and theory embodied into Monte Carlo 
%%%TT
programs treated as a black box. 
One could see that a typical signature
of any given $\tau$ decay channel is matching rather poorly 
the publicly available Monte Carlo predictions. This 
should be of no surprise as efforts to compare data with
predictions were completed  for the last time 
in late 90's by the ALEPH and CLEO collaborations.
The resulting hadronic currents were afterwards implemented 
in~\cite{Golonka:2003xt}.
% In fact, 
%it was just the archivization  of the  experimental work done 
%by the ALEPH and CLEO collaborations  in late 90's only. 
Since that time no efforts to prepare a complete parametrisation of $\tau$ 
decay simulation for the public use were undertaken seriously.

There is another important message which can be drawn from these comparisons.
Starting from a certain precision level, the study of a given decay mode can 
not be separated from the discussion of others. In the distributions aimed 
at representing the given decay mode, a contribution from the other 
$\tau$ decay modes can be large, up to even 30\%. 

It may be less clear that experiments differ significantly in the way how 
they measure individual decay modes. For instance, ALEPH 
produced $\tau$ samples free of the non-$\tau$ backgrounds,
but, on the other hand, strongly boosted, making the
reconstruction of some angles in the hadronic system more difficult. 
This is important and affects  properties of the decay models 
which will be used for a parametrisation. In particular, when the 
statistics is small, possible fluctuations may affect the picture 
and there are not enough data to complete an estimate of the systematic errors.
In this case, details of the description of the hadronic current,
as the inclusion of intermediate resonances, are not important.
Let us consider, as an example, 
$\tau^- \to K_S^0 \pi^- \pi^0 \nu_{\tau}$. 
The matrix element in the ALEPH parametrisation is saturated by 
 $\rho^-\to \pi^- \pi^0 $ and
$K^{*0} \to K_S^0 \pi^0$, and a 
similar parametrisation is used for $K^{*-}\to K_S^0 \pi^-$. 
In practice, the contribution of the $\rho$ is more significant in 
the ALEPH parametrisation in contrast to the CLEO one where 
the  $K^*$ dominates.
One has to admit that at the time when both collaborations were 
preparing their parametrisations  to be used in 
TAUOLA, the data samples of both experiments were rather
small and the differences were not of much significance. This can, however, 
affect possible estimates of backgrounds for searches of rare decays, 
e.g. of $B$ mesons at LHCb.\footnote{LHCb performed MC studies 
for $B_s^0 \to \mu^+ \mu^-$ and the 
radiative decays $B^0 \to K^* \gamma$ and $B_s^0 \to \phi \gamma$, 
but $\tau$ decays have not yet been taken into account. 
These results are not public and exist only as internal documents
LHCB-ROADMAP1-002 and LHCB-ROADMAP4-001.} 

Let us now go point by point and discuss examples of Monte Carlo programs 
and fitting strategies. We will focus on subjects requiring most 
attention and future work. We will  review the theoretical 
constraints which are useful  in the construction of the models 
used for the data description. 

\begin{center}
\begin{figure}
\includegraphics[width=0.48\textwidth]{bab51.eps}
%\includegraphics[width=0.40\textwidth]{bab51.eps}
\caption{Invariant mass of five charged particles for 
$\tau^- \to  3h^-2h^+\nu_{\tau}$ at BaBar.}
\label{fig:tau51}
\end{figure}

\begin{figure}
%\begin{tabular}{cc}
%\includegraphics[width=0.40\textwidth]{bab52.eps} %&
%\\
%\includegraphics[width=0.40\textwidth]{bab53.eps}
\includegraphics[width=0.48\textwidth]{bab52.eps} %&
\\
\includegraphics[width=0.48\textwidth]{bab53.eps}
%\end{tabular}
\caption{Invariant mass distributions
 for $\tau^- \to  3h^-2h^+\nu_{\tau}$ at BaBar. Points with 
error bars are the data: Upper panel -- $h^+h^-$; the unshaded and shaded 
histograms are the signal and background predicted by MC.  
Lower panel -- $2\pi^+ 2\pi^-$; the solid line is a fit to the data
using a second-order polynomial (dashed line) for the background and 
a Breit-Wigner convoluted with a Gaussian for the peak region.} 
\label{fig:tau52}
\end{figure}
\end{center}


\subsection{Phase space}
Because of the relatively low multiplicity of final state particles, 
it is possible
to separate the description of $\tau$ production and decay into 
segments describing the matrix elements and the phase space.
In the phase space no approximations are used, contrary to the
 matrix elements where all 
approximations and assumptions reside. 
The description of the phase space 
used in TAUOLA is given in detail in~\cite{Jadach:1993hs}. The description 
of the phase space for $\tau$ production is given in~\cite{kkcpc:1999}. Thanks 
to conformal symmetry it is exact for an arbitrary number of 
photons. Using exponentiation, 
see, for example, Yennie-Frautchi-Suura~\cite{yfs:1961}, 
the phase space description 
can be exact and the matrix element can be refined order by order. 
For radiative corrections in the decay PHOTOS can be 
used. Its phase space is described, for example, 
in the journal version of~\cite{Nanava:2006vv} and is exact.
Approximations are made in the matrix element only.
%%%TT
Benchmark comparisons\footnote{The purpose of this type of tests may
vary. If two 
programs differ in their physics assumptions, it may help to control
the physics precision. If the physics assumptions are identical, but
the technical constructions differ, then the comparison checks the
correctness of the implementation of the algorithm. Finally, the
comparison of results from the same program, but installed on
different computers, may check the correctness of the code's
implementation in new software environments. Such comparisons, or just
the data necessary for comparisons, will be referred to as physical,
technical and installation benchmarks, respectively. They are
indispensable for the reliable use of Monte Carlo programs.} 
with other calculations, which are actually 
based on second-order matrix elements and exponentiation, found
%%%TT
excellent agreement \cite{Golonka:2005pn,Golonka:2006tw}. 
%The approximation needs to be present in the way how presamplers for the 
%generation chains of each charged final state particle are matched together.
%They are embedded in the way how the matrix element (or if not available, 
%the interference
%weight) is used for configurations with multiple photon 
%emissions.
%~\cite{Nanava:2006vv}.

\subsection{Spin effects}
The lifetime of the $\tau$ lepton is 
orders of magnitude larger not only than its formation time in high
energy experiments, but also than the time scale of all phenomena
related to higher-order corrections such as bremsstrahlung. 

The separation of $\tau$ production and decay is excellent due to the small width 
of the  $\tau$ lepton. Its propagator can be well approximated by
a delta function for 
phase space and matrix elements.
% thus for the final state
%cross section as well. 
%
The cross section for the process 
$ f \bar f \to \tau^+\tau^-  Y; 
\tau^+ \to X^+ \bar{\nu}_{\tau}; \tau^- \to l^- \nu_l \nu_{\tau}$ reads
 \[ 
{\rm d} \sigma = \sum_{spin }|{\cal M}|^2 
{\rm d}\Omega= \sum_{spin }|{\cal M}|^2 {\rm d}\Omega_{\rm prod} \; 
{\rm d}\Omega_{\tau^+} \; {\rm d}\Omega_{\tau^-}\,,\]
where $Y$ and  $X^{+}$ stand for particles produced together with
 the $\tau^+\tau^-$ and in the $\tau^+$ decay, respectively;
${\rm d}\Omega$,  ${\rm d}\Omega_{\rm prod}$, ${\rm
d}\Omega_{\tau^+}$, 
${\rm d}\Omega_{\tau^-}$ denote
the phase space in the original process, in production and
decay, respectively.

%\item 
This formalism looks simple,
 but because of the over 20 $\tau$ decay channels 
there are more than 400 distinct processes.
% Also pictures of production 
%and decay are mixed.
%\item 

Let us write the spin amplitude separated into the parts for 
$\tau$ pair production and decay:
\[ 
{\cal M}=\sum_{\lambda_1\lambda_2=1}^2{\cal
      M}_{\lambda_1\lambda_2}^{\rm prod} \; 
      {\cal M}_{\lambda_1}^{\tau^+}{\cal M}_{\lambda_2}^{\tau^-}.
\]
%\item 
After integrating out the $\tau$ propagators, the formula for the 
cross section can be rewritten as
 \[ 
{\rm d} \sigma = \Bigl(\sum_{spin }|{\cal M}^{\rm prod}|^2 \Bigr)
 \Bigl(\sum_{spin }|{\cal M}^{\tau^+}|^2 \Bigr)
 \Bigl(\sum_{spin }|{\cal M}^{\tau^-}|^2 \Bigr) 
\]
\[
\hskip 7 mm \times\, wt \; {\rm d}\Omega_{\rm prod} \; {\rm
d}\Omega_{\tau^+} \; {\rm d}\Omega_{\tau^-},
\]
%%%%%\item 
where
\[ 
wt=  \Bigl( \sum_{i,j=0,3} R_{ij} h_+^i h_-^j \Bigr),
\]
\[  
R_{00} =1,~~~ <wt>=1,~~~ 0 \le wt \le 4.
\]  
%{\bf Define $R_{ij}$ and $h_+^i, h_-^j$....} 

$R_{ij}$ can be calculated from  ${\cal M}_{\lambda_1\lambda_2}$, $h_+^i$ and $h_-^j$ from $ {\cal M}^{\tau^+}$ 
and $ {\cal M}^{\tau^-}$, respectively. 
%\item 
Bell inequalities (related to the Einstein-Rosen-Podolsky 
paradox~\cite{Einstein:1935rr}) tell us that in general it is
impossible to rewrite $wt$ in 
 the following factorised form, $wt^{\rm factorized}$:
\[ 
 wt \ne wt^{\rm factorized}= \Bigl( \sum_{i,j=0,3} R^A_{i} h_+^i \Bigr)\Bigl( \sum_{i,j=0,3} 
R^B_{j}  h_-^j \Bigr),
\]
where $R^A_{i}$ and $R^B_{j}$ are four-component objects 
calculated from variables of the process of 
$\tau$ pair production.
In the Monte Carlo construction it is thus impossible to generate a
$\tau^+$ $\tau^-$ pair, where each of the two is in some quantum
state, and later to perform the decays of the $\tau^+$ and the
$\tau^-$ independently. 
%
This holds at all orders of the perturbative expansion. 
$\tau$ production and decay 
%%%TT
are correlated through spin effects, 
which can be represented by the well-behaved factor $wt$ introduced
previously. 
%, $0 \le wt \le 4$ and its average $<wt>=1$.  
The above formulae do not lead to any loss of precision and hold in
presence of radiative corrections as well. Different options for the
formalism, based on these expressions, are used in Monte Carlo
programs and are basically well founded.
This should be confronted with processes where instead of $\tau$
leptons short-lived intermediate states are considered. 
Then, in general, ambiguities appear
and corrections proportional to the ratio of the resonance 
width to its mass (or other energy scales of the process resulting,
for example, from cut-offs) must be included. 
Interfering background diagrams may cause additional 
problems.  For details we refer 
to~\cite{koralb:1985,Jadach:1990mz,kkcpc:1999}.

\subsection{$\tau$ lepton production}
KORALB was 
published~\cite{Jadach:1984iy,Jadach:1985ac} more than 
{\it twenty years ago}. 
It included first-order QED corrections and complete  mass and spin 
effects. It turned out to be very useful, 
and still remains in 
broad use. On the other hand, some of its ingredients
are outdated
and do not match the present day requirements, even for technical tests.
For example the function PIRET(S), which describes the
real part of the photon hadronic vacuum polarisation as measured by the data
collected until the early 80's should be replaced by one of the 
new precise codes (see Section~\ref{sec:4} for details).
% Obviously, this function has to be replaced
%by a more modern one. A possible choice can be the function
%{\tt REPI} of ref.~\cite{Burkhardt:2005se} because it 
%has similar functionality as
%{\tt PIRET} of {\tt KORALB}.

Unfortunately, this replacement
does not solve all
normalisation problems of KORALB. For example, it is well known 
that the one-loop corrections are not sufficient.
The two major improvements which were developed  
during the LEP era are the introduction 
of higher-order QED corrections into Monte Carlo simulation and a better way
to combine loop corrections with the rest of the field theory calculations.
%It was found to be safe  to sum contributions of loop 
%corrections into photon (and $Z$) propagators. Then, 
%terms of all, but incomplete, orders of perturbation expansion
%are taken into account. That is 
%why significant efforts were needed to justify the 
%approach~\cite{Bardin:1999ak}.
For energies up to  10 GeV  (typical of the $B$ factories), 
the KKMC
Monte Carlo~\cite{kkcpc:1999} provides a realisation of the above 
improvements. This program includes higher-order QED
matrix elements with the help of exclusive exponentiation, and 
explicit matrix elements up to the second order. Also in this case 
the function calculating the 
vacuum polarisation must be replaced by a version appropriate for low energy
(see Section~\ref{sec:4}).
 
Once this is completed, and if the 
two-loop photon vacuum polarisation can be neglected, 
KORALB and KKMC can form a base for tests and studies 
of systematic errors for cross section normalisations 
at low energies. Using a strategy 
similar to the one for Bhabha scattering~\cite{Jadach:1991pj},
the results obtained in~\cite{Banerjee:2007is,Jadach:2000ir}
allow to expect a precision  of 0.35--0.45\% using KKMC
at Belle/BaBar energies.
 Certainly, a precision tag similar to that for 
linear colliders can also be achieved for lower energies. 
Work beyond ~\cite{Banerjee:2007is} and explained in that paper
would then be necessary.




\subsection{Separation into leptonic and hadronic current}

 The matrix element used in TAUOLA for semi-leptonic decays, 
 $\tau(P,s)\rightarrow\nu_{\tau}(N)X$,
\begin{equation}
{\cal M}=\frac{G}{\sqrt{2}}\bar{u}(N)\gamma^{\mu}(v+a\gamma_{5})u(P)J_{\mu}
\end{equation}
requires the knowledge of the hadronic current  $J_{\mu}$. The 
expression is easy to manipulate. One obtains:
\begin{eqnarray}
|{\cal M}|^{2}&=& G^{2}\frac{v^{2}+a^{2}}{2}
( \omega + H_{\mu}s^{\mu} ), \nonumber\\
\omega&=&P^{\mu}(\Pi_{\mu}-\gamma_{va}\Pi_{\mu}^{5}), \hskip 5 mm, \nonumber\\
H_{\mu}&=&\frac{1}{M}(M^{2}\delta^{\nu}_{\mu}-P_{\mu}P^{\nu})(\Pi_{\nu}^{5}-
\gamma_{va}\Pi_{\nu}), \nonumber\\
\Pi_{\mu}&=&2[(J^{*}\cdot N)J_{\mu}+(J\cdot N)J_{\mu}^{*}-(J^{*}\cdot J)
N_{\mu}], \nonumber\\
\Pi^{5\mu}&=&2~ {\rm Im} ~\epsilon^{\mu\nu\rho\sigma}
J^{*}_{\nu}J_{\rho}N_{\sigma}, \hskip 5 mm \nonumber\\
\gamma_{va}&=&-\frac{2va}{v^{2}+a^{2}}.
\end{eqnarray}
If the $\tau$ coupling is 
$v+a\gamma_{5}$ and $m_{\nu_{\tau}} \neq 0$ 
is allowed, one has
to add  to $\omega$ and $H_{\mu}$:
\begin{eqnarray}
\hat{\omega}&=&2\frac{v^{2}-a^{2}}{v^{2}+a^{2}}
m_{\nu}M(J^{*} \cdot J), \nonumber \\
\hat{H}^{\mu}&=&-2\frac{v^{2}-a^{2}}
{v^{2}+a^{2}}m_{\nu}~ {\rm Im}~\epsilon^{\mu\nu\rho\sigma}
J_{\nu}^{*}J_{\rho}P_{\sigma}.
\end{eqnarray}
The expressions are useful for Monte Carlo applications and 
are also calculable from first principles. The resulting
expression can be used to the precision level of the order of 0.2--0.3\%. 

In contrast to other parts, the hadronic current $J_{\mu}$ still
can not be calculated reliably from first principles. Some theoretical 
constraints need to be fulfilled, but in general it has to be obtained 
from experimental data. We will return to this point later 
(see Section \ref{hadroniccurrents}). 
 
\subsection{Bremsstrahlung in decays}
The PHOTOS Monte Carlo is widely used for generation
 of radiative corrections 
in cascade decays, starting from the early 
papers~\cite{Barberio:1990ms,Barberio:1994qi}. With time the precision 
of its predictions improved significantly, but the main principle 
remains the same.
Its algorithm is aimed to modify the content of the event record 
filled in with 
complete cascade decays at earlier steps of the generation.  {\tt PHOTOS} 
modifies the content of the event record; it adds additional photons 
to the decay vertices and at the same time modifies the kinematic 
configuration of other decay products.

One could naively expect that this strategy is bound 
to substantial approximations. 
 However, the algorithm is compatible with NLO calculations,
 leads to a complete coverage of the phase space for multi-photon 
final states and provides correct distributions in soft photon
limits. For more details of the program organisation and its phase 
space generation we address the reader to~\cite{Nanava:2006vv}.

The changes introduced over the last few years 
into the PHOTOS Monte Carlo program itself were rather small and
 the work concentrated on its theoretical foundations.
 This wide and complex subject goes beyond the scope of this Review
 and the interested reader can consult~\cite{Was:2008zz}, where  
 some of the topics are discussed.  
Previous tests  of two-body decays of the $Z$  
into a pair of charged leptons~\cite{Golonka:2006tw}
and a pseudoscalar $B$  into a pair of scalars~\cite{Nanava:2006vv} 
were recently supplemented~\cite{Photos_tests} with the study 
of $W^\pm \to l^\pm \nu \gamma$. 
The study of the process $\gamma^* \to \pi^+\pi^-$ is on-going~\cite{Xu}.
In all of these cases a universal kernel of PHOTOS was replaced 
with the one matching the exact first-order matrix element. 
In this way terms for the NLO/NLL 
level are implemented. The algorithm covers the full multi-photon 
phase space and it is exact in the infrared region of the phase space. 
 One should
not forget that PHOTOS generates weight-one events.


The results of all tests of PHOTOS with an NLO kernel are
at a sub-per mill level. No differences with 
benchmarks were found, even for samples of $10^9$ events. 
When simpler physics assumptions were used,
differences  between total rates at sub-per mill level 
 were observed  
 or they were matching 
a precision of the programs used for tests. 

This is very encouraging and points to the possible extension of the
approach beyond (scalar) QED, and in particular to QCD and/or models
with phenomenological Lagrangians for interactions of photons with
ha\-drons. For this work to be completed, spin amplitudes have 
to be further studied~\cite{vanHameren:2008dy}.
 
The  refinements discussed above
affect the practical side of simulations for $\tau$ physics only indirectly. 
Changes  in the kernels necessary for NLO  may remain as 
options  for tests only. They are available from the PHOTOS 
web page~\cite{Photos_tests}, but are not recommended for wider use. 
The corrections are small, and distributions visualising their size are available.
On the other hand, their use could be perilous, as it requires control of the
decaying particle spin state. 
%%%TT
It is known (see, e.g., ~\cite{Was:2004bk}) that this is not easy
%@Article{Was:2004bk,
%     author    = "Was, Z.",
%     title     = "{Observables with tau leptons at LHC and LC: Structure of
%                  event records  and Monte Carlo algorithms}",
%     journal   = "Nucl. Instrum. Meth.",
%     volume    = "A534",
%     year      = "2004",
%     pages     = "260-264",
%     eprint    = "hep-ph/0402129",
%     archivePrefix = "arXiv",
%     doi       = "10.1016/j.nima.2004.07.097",
%     SLACcitation  = "%%CITATION = HEP-PH/0402129;%%"
%}
because of technical reasons. 


We will show later that radiative corrections do not provide a limitation 
in the quest for  improved precision of matching theoretical models 
to experimental data until issues discussed in subsection \ref{Thechallenges} 
are solved.

\subsection{Hadronic currents}
\label{hadroniccurrents}
So far all discussed contributions to the predictions were found 
to be controlled to the precision level of 0.5\% with respect to the decay 
rate under study.\footnote{This $0.5\%$ uncertainty is for QED 
radiative effects.  One should bear in mind other mechanisms involving 
the production of
photons, like, for example, the decay channel 
$\omega\to \pi \gamma$, which occurs with a probability of 
$(8.28 \pm 0.28)$\% and  does not belong to 
the category of radiative corrections. }



This is not the case for the hadronic current, which is the main source 
of our difficulties.
It can not be obtained from perturbative QCD as the energy scales involved are too
small. On the other hand, for the low energy limits the scale is too
large. 
%%%TT
Despite these difficulties one can obtain a theoretically clear object
if enough effort is devoted. This may lead to a better understanding of
the boundaries of the perturbative domain of QCD as well.

The unquestionable property which hadronic currents must fulfil 
is Lorentz invariance. For example, if the final state consists of
three scalars with momenta $p_1$, $p_2$, $p_3$, respectively, it must
take the form 
{\small{
\begin{eqnarray}
\label{fiveF}
&J^\mu & =N \bigl\{T^\mu_\nu \bigl[ c_1 (p_2-p_3)^\nu F_1  + c_2 (p_3-p_1)^\nu
 F_2  \nonumber \\
 &+& c_3  (p_1-p_2)^\nu F_3 \bigr] + c_4  q^\mu F_4  -{ i  c_5\epsilon^\mu_{.\ \nu\rho\sigma} \over 4 \pi^2 f_\pi^2}     
 p_1^\nu p_2^\rho p_3^\sigma F_5      \bigr\}\,,
\end{eqnarray}
}}
%\begin{equation}
%\label{fiveF}
%\begin{array}{l}
%J^\mu =N \bigl\{T^\mu_\nu \bigl[ c_1 (p_2-p_3)^\nu F_1  + c_2 (p_3-p_1)^\nu
% F_2  \nonumber \\
%+ c_3  (p_1-p_2)^\nu F_3 \bigr] + c_4  q^\mu F_4  -{ i  c_5\epsilon^\mu_{.\ \nu\rho\sigma} \over 4 \pi^2 f_\pi^2}     
% p_1^\nu p_2^\rho p_3^\sigma F_5      \bigr\},
%\end{array}\end{equation}
where  $T_{\mu\nu} = g_{\mu\nu} - q_\mu q_\nu/q^2$ is the transverse
projector and $q=p_1+p_2+p_3$.
The functions $F_i$ depend on three variables that can be chosen 
as $q^2=(p_1+p_2+p_3)^2$
and two of the following three, $s_1=(p_2+p_3)^2$, $s_2=(p_1+p_3)^2$, 
$s_3=(p_1+p_2)^2$.
 This form is obtained  from Lorentz invariance only. 

Among the first four  hadronic structure functions 
($F_1$, $F_2$, $F_3$, $F_4$),   only three  are independent.
%, as it can be deducted from the chiral symmetry. 
We leave the structure function $F_4$ in the basis because, 
neglecting the pseudoscalar resonance production mechanism, 
the contribution due to $F_4$ is negligible ($\sim
m_{\pi}^2/q^2$)~\cite{GomezDumm:2003ku} and (depending on the decay
channel) one of $F_1$, $F_2$ and $F_3$ drops out, exactly as it is in
TAUOLA since long. 

In each case, 
the number of independent functions  is  four (rather than  five) 
and not larger 
than the dimension of our space-time.
That is why  
the projection operators can be defined, for two- and three-scalar 
final states.
Work in that direction has already been done in 
Ref.~\cite{Kuhn:1992nz} and then 
implemented in tests of TAUOLA too.
Thanks to such a method, hadronic currents can be obtained from data 
without any 
need of phenomenological assumptions. Since long such methods  were useful for 
data analysis, but only in part. Experimental samples were simply too small.

At present, for high statistics and precision the method may  be revisited. 
That is why it is of great interest to verify whether 
detector deficiencies will invalidate the method or if adjustments due
to incomplete phase space coverage are necessary. We will return to that 
question later. In the mean time let us return to other theoretical 
considerations which  constrain the form of hadronic currents, 
but not always  to the precision of today's data.


\subsection{The resonance chiral approximation and its result for the currents}
Once the allowed Lorentz structures are determined and a proper 
minimal set of them is chosen, one should impose the QCD symmetries 
valid at low energies. The chiral symmetry of massless QCD allows
 to develop an 
effective field theory description valid for momenta much smaller 
than the $\rho$ mass, $\chi PT$~\cite{Gasser:1983yg,Gasser:1984gg}.\\ 
Although $\chi PT$ cannot provide predictions valid over the full $\tau$ 
decay phase space, it constrains the form and the normalisation of the 
form factors in such limits.\\
The model, proposed in~\cite{Kuhn:1990ad} for $\tau$ decaying to pions
  and used also 
 for extensions to other 
decay channels, employs weighted products of Breit-Wigner functions to 
take into account resonance exchange. The form factors
used there have the right chiral limit at LO. However, as it 
was demonstrated in~\cite{GomezDumm:2003ku},
they  do not reproduce the  NLO chiral limit.
 \\
The step towards incorporating the right low-energy limit up to  NLO
 and the contributions from meson resonances which reflect the experimental
  data was
  done within Resonance Chiral Theory 
($R\chi T$~\cite{Ecker:1988te,Ecker:1989yg}). 
 The current state-of-the-art for the 
hadronic form factors ($F_i$) appearing in the $\tau$ decays is 
described in~\cite{Dumm:2009kj,Dumm:2009va}. Apart from the correct low energy
properties, 
it includes the right falloff~\cite{Brodsky:1973kr,Lepage:1980fj} 
at high energies.\\ 
%$R\chi T$ also takes advantage of the $N_C\to\infty$ limit of $QCD$ \cite{Nc}.
The energy-dependent imaginary parts
in the propagators 
of the vector and the axial-vector mesons, 
$~1/(m^2-q^2- i m\Gamma(q^2))$, were
  calculated in ~\cite{GomezDumm:2000fz}  at one-loop,
exploiting the optical theorem that relates the appropriate hadronic 
matrix elements of $\tau$ decays and the cuts with on-shell mesons 
in the (axial-) vector-(axial-) vector correlators.\\
This formalism has been shown to successfully describe the invariant
mass spectra of
experimental data in $\tau$ decays for the following hadronic 
systems: $\pi\pi$ \cite{Guerrero:1997ku,Pich:2001pj,Pich:2002ne}, 
$\pi K$ \cite{Jamin:2006tk,Jamin:2008qg}, 
$3\pi$ \cite{GomezDumm:2003ku,Dumm:2009kj,Dumm:2009va,Roig:2008xt} and 
$KK\pi$\cite{Dumm:2009kj,Roig:2008xt}. 
Other channels will be worked out along the same lines.
\\
It has already been checked that the $R\chi T$ results 
provide also a good description of the three-meson processes 
$\Gamma(\tau \to 3\pi \nu_\tau)$~\cite{Barate:1998uf} and 
$\sigma(e^+e^-\to KK\pi)_{I=1}$~\cite{Aubert:2007ym}.\\
Both the spin-one resonance widths and the form factors of 
the decays $\tau^- \to (\pi\pi,\ \pi K,\ 3\pi,\ KK\pi)^- \nu_\tau$ 
computed within $R\chi T$ are being implemented in TAUOLA only now.

Starting from a certain precision level, the predictions, 
like  the ones presented above,
may turn out to be not sufficiently precise.
Nonetheless, even in such a case they can provide some essential 
constraints on the form of the functions $F_i$. 
Further refinements will  require  large and combined efforts of 
experimental and theoretical physicists. We will elaborate on possible 
technical solutions later in the review.
Such attempts turned out to be difficult in the past
and a long time was needed 
for parametrisations given 
in~\cite{Golonka:2003xt} to become
public. Even now they are semi-official and are not based on the
final ALEPH and/or CLEO data.


\subsection{Isospin symmetry of the hadronic currents}

If one neglects quark masses, QCD is invariant under a 
transformation replacing quark flavours.
As a consequence, hadronic currents describing vector $\tau$ decays
($2\pi,\ 4\pi,\ \eta\pi\pi,\ \ldots$) and low energy
$e^+e^-$ annihilation into corresponding iso\-vector final states are
related and can be obtained from one another~\cite{Tsai:1971vv,Thacker:1971hy}.
This property, often referred to as conservation of the vector current
(CVC) in $\tau$ decays, results in the possibility to predict invariant
mass distributions of the hadronic system, as well as the corresponding
branching fractions in $\tau$ decays using $e^+e^-$ data.
 A systematic check of these predictions showed
that at the (5--10)\% level they work rather well~\cite{Eidelman:1990pb}.

In principle, the corrections due to mass  and
charge differences between $u$ and $d$ quarks are not expected 
to provide significant and impossible to control 
effects~\cite{Cirigliano:2001er,Cirigliano:2002pv}. However,
the high-precision data of the CLEO \cite{Anderson:1999ui},
 ALEPH \cite{Schael:2005am}, 
 OPAL \cite{Ackerstaff:1998yj},  Belle \cite{Fujikawa:2008ma},
 CMD-2 \cite{Akhmetshin:2003zn,Aulchenko:2006na,Akhmetshin:2006wh,Akhmetshin:2006bx}, SND \cite{Achasov:2006vp} and KLOE \cite{:2008en}
 collaborations in the $2\pi$ channel challenged this statement, and
 as it was shown 
 in~\cite{Davier:2002dy,Jegerlehner:2003qp,Jegerlehner:2003rx,Ghozzi:2003yn,Davier:2003pw,Jegerlehner:2008zza,Jegerlehner:2009ry} that 
 the
spectral functions for $\tau^- \to \pi^-\pi^0\nu_{\tau}$ significantly differ
 from those obtained using $e^+e^- \to \pi^+\pi^-$ data.
 Some evidence for a similar
discrepancy is also observed in the  $\tau^- \to 2\pi^-\pi^+\pi^0\nu_{\tau}$
decay  \cite{Davier:2005xq,Druzhinin:2007cs,Czyz:2008kw}.
 This effect remains unexplained. The magnitude of the isospin-breaking
corrections has been updated recently, making the discrepancy in the
$2\pi$ channel smaller~\cite{Davier:2009ag}.

These CVC based relations were originally used in the TAUOLA
form factors parametrisation,
but they were often modified to improve fits to the data.
Let us point here to an example where experimental $e^+e^-$ data were
used for the model of the $\tau \to 4\pi \nu_{\tau}$ decay 
channels~\cite{Bondar:2002mw}.
In this case, only a measurement of the distribution
  in the total invariant mass
of the hadronic system was available. 
This is not enough to fix the distribution over the multidimensional phase space.
For other dimensions one had to rely on
phenomenological models or other experiments.
In the future, this may not be necessary, but will always 
remain as a method of benchmarks construction.



%\subsection{Here the problem resides}
\subsection{The challenges}
\label{Thechallenges}

As we have argued before, refined techniques for fits,
 involving simultaneous fits to many $\tau$ decay channels,  are necessary
 to improve the phenomenological description 
of $\tau$ decays. Complex backgrounds (where each channel contributes 
to signatures of other decay modes as well),
different sensitivities of experiments for measurements of some 
angular distributions within the same hadronic system, and sometimes even 
an incomplete reconstruction of final states, 
are the main cause of this necessity.
Moreover, theoretical models based on the Lagrangian approach
simultaneously describe more than one  $\tau$ decay channel with the same
set of parameters, and only simultaneous fits allow to establish
 their experimental constraints in a consistent way.
Significant efforts are thus necessary and
  close collaboration between phenomenologists and 
experimental physicists is indispensable.
 As a result,
techniques of automated calculations of hadronic currents may become 
 necessary~\cite{Korchin}. 

\subsection{Technical solutions for fits}

For the final states of up to three scalars, 
the use of projection operators~\cite{Kuhn:1992nz} is popular since 
long~\cite{Davier:2005xq}. 
%%%TT
It enables, at least in principle, to obtain form factors used in
hadronic currents directly from the data, for one scalar function
defined in Eq.~(\ref{fiveF}) at a time. Only recently experimental samples became sufficiently large.
However, to exploit this method one may have to 
improve it first by systematically including the effects of a limited detector 
acceptance.  
Implementation of the projection operators into packages like
MC-TESTER~\cite{Davidson:2008ma} may be useful. Efforts
in that direction are being pursued now\footnote{ This may help to
embed the method in the modern software for fits, see, e.g.,
\cite{Brun:2008zza}.}~\cite{OlgaArtur}.

On the theoretical side one may need to choose  
predictions from many  models, before a sufficiently good
agreement with data will be achieved. 
Some automated methods of calculations may then become 
useful~\cite{Automated}. This is especially important 
for hadronic multiplicities larger than three,  when 
projector operators have never been defined.

Certain automation of the methods is
thus advisable. To discriminate from the broad spectrum of choices, 
new methods of data analysis may become useful~\cite{Hocker:2007ht}. 
Such methods may require simulating samples
of events where several options for the matrix element calculation 
are used simultaneously.\footnote{Attempts to code such methods into TAUOLA, 
combined with programs for $\tau$ pair production and experimental 
detector environment, were recently
performed~\cite{VladimirTomasz}, but they were applied so far 
as prototypes only, see Fig. 1  of Ref.~\cite{Was:2009iy}.}
 


The neutrino coming from $\tau$ decays escapes detection 
and as a result the $\tau$ rest frame 
can not be reconstructed. Nevertheless, as was 
shown in Ref.~\cite{Kuhn:1992nz}, angular distributions can be used 
for the construction of projection operators, which allow the extraction
 of the hadronic structure functions from the data.
 This is possible as they
depend on $s_1$, $s_2$ and $q^2$ only.

A dedicated module for the MC-TESTER~\cite{Davidson:2008ma}, 
implementing the moments of different angular functions defined in 
Eqs.~(39)--(47) of Ref.~\cite{Kuhn:1992nz}, is under development.
The moments are proportional to combinations of the type 
$\alpha |\mathrm{F_{i}}|^{2} + \beta |\mathrm{F_{j}}|^{2} + 
\gamma \mathrm{Re(F_{i} F_{j}^{*}})$,
where the coefficients $\alpha$, $\beta$ and $\gamma$ are functions 
of hadron four-momentum components in the hadronic rest frame.
Preliminary results obtained with large statistics of five million 
$\tau \rightarrow a_{1} \nu_{\tau} \rightarrow 3\pi \nu_{\tau}$ decays, 
and assuming vanishing $\mathrm F_{3}$ and $\mathrm F_{5}$ form factors,
show that it is possible to extract 
$\mathrm |F_{1}|^{2}$, $\mathrm |F_{2}|^{2}$ and 
$\mathrm |F_{1}\cdot F_{2}^{*}|^{2}$ as functions of 
$s_{1}$, $s_{2}$ and $Q^{2}$. 
This extraction requires solving a set of equations. Since the
solution is sensitive to the precision of the estimation of 
the moments entering the equation, large data samples of the order of 
$O (10^{6}-10^{7})$ are necessary. The calculation of the moments 
also requires the knowledge of the initial $\sqrt{s}$ of the 
$\tau$ pair, which makes the analysis sensitive to 
initial state radiation (ISR) effects. The same studies show 
that the analysis is easier if one, instead of extracting
 the form factors $\mathrm |F_{i}|^{2}$, compares the 
moments obtained from the experimental data with theoretical
predictions. Such a comparison does not require repetition of the Monte Carlo
simulation of $\tau$ decays with different form factors, and only 
the calculation of combinations
of  $|\mathrm{F_{i}}|^{2}$ and $\mathrm{Re(F_{i} F_{j}^{*}})$ is necessary. 
This is much simpler than 
comparing the kinematic distributions obtained from data
with distributions coming from Monte Carlo simulations
  with various theoretical models. 
Further complications, for example, due to the presence 
of an initial state bremsstrahlung or an incomplete acceptance of decay 
phase space, were not yet taken into account. 

\subsection{Prospects}

 
Definitely the improvements of $\tau$ decay simulation packages and 
fit strategies
are of interest for phenomenology of low energy.
As a consequence, their input for such domains 
like phenomenology of the muon $g-2$ or $\alpha_{\rm QED}$,
$\alpha_{\rm QCD}$ 
and their use in constraints of new physics would improve. 
%This topic is covered in references \cite{Simon??,Hagiwara} 
%and/or Section \ref{section} of present review. 

In this section, let us argue if  possible benefits for LHC 
phenomenology may arise from a better 
understanding of $\tau$ decay channels in measurements as well. In 
the papers~\cite{Bower:2002zx,Desch:2003rw} it was shown that spin effects 
can indeed be useful to measure 
properties of the Higgs boson such as parity. Moreover, 
such methods were verified to work well when detector effects 
as proposed for a future linear collider 
were taken into account.
Good control of the decay properties is helpful. For example,
in Ref.~\cite{Privitera:1993pr}  it was shown that for
the $\tau \to \ a_1 \nu_{\tau} \to 3\pi\nu_{\tau}$ decay the
sensitivity to the $\tau$ polarisation increases about four times when
all angular variables are used compared with the usual ${\rm
d}\Gamma/{\rm d}q^2$, see also~\cite{Kuhn:1995nn}. 

Even though $\tau$ decays provide some of the  most prominent signatures
for the LHC physics program, see, e.g., Ref.~\cite{RichterWas:2009wx},
for some time it was  expected that methods exploiting detailed
properties of $\tau$ cascade decays are not practical
for LHC studies. Thanks to
efforts on reconstruction of $\pi^0$ and  $\rho$ 
invariant mass peaks, this opinion evolves. 
Such work was done for studies of the CMS ECAL
detector inter-calibration~\cite{Bayatian:2006zz}, and 
in a relatively narrow $p_T$ range (5--10 GeV) some potentially 
encouraging results were obtained.
Some work in context of searches for new particles 
started recently~\cite{Nattermann:2009gh}. There, improved knowledge of
distinct $\tau$ decay modes may become important at a certain point.

One can conclude that the situation is  similar to that at the start
 of LEP, and
 some control of all $\tau$ decay channels is  important.
Nonetheless,  only if detector  studies of $\pi^0$ and  $\rho$ reconstruction
will provide positive results, the gate to improve the 
sensitivity of $\tau$ spin measurements with most of its decay modes, 
as at LEP~\cite{:2008zzm,Aad:2009wy,Heister:2001uh}, will be open. 
At this moment, however, it is difficult to judge about the importance
of such improvements in the description of $\tau$ decays for LHC perspectives.
The experience of the first years of LHC must be consolidated first.
In any case such an activity is important for the physics of future 
Linear Colliders.  



\subsection{Summary}

We have shown  that the most urgent challenge in the quest
for a better understanding 
of $\tau$ decays is the development of efficient techniques for fitting 
multidimensional distributions, which take into account
 realistic detector conditions.
  This includes cross 
contamination of different $\tau$ decay modes, their respective signatures 
 and detector acceptance effects, which have to be simultaneously 
taken into account
 when fitting experimental data. 
 Moreover, at the current experimental precision, theoretical concepts have 
to be reexamined. In contrast to the
past, the precision of predictions based on chiral
Lagrangians and/or isospin symmetry 
can not be expected to always match the precision of the data. 
The use of model-independent data analyses
should be encouraged whenever possible in 
realistic conditions.

Good understanding of $\tau$ decays is crucial for understanding 
the low energy regime of strong interactions and
the matching between the non-perturbative and the perturbative domains. 
Further work on better simulations of $\tau$ decays at the LHC is needed 
 to improve its potential
to study processes of new physics, especially in the Higgs sector.
In addition, an accurate simulation of $\tau$ decays is important   
for the control 
of backgrounds for very rare decays.
 For the project to be successful, this should lead  
 to the encapsulation of our  knowledge on $\tau$ decays
  in form of 
a Monte Carlo library to be used by  low-energy as well as  
high-energy applications.


