\section{Analysis}
\label{sec:analyse}
After running the code, we need to combine the various {\tt
which\_pieces} into physical results that we will want to use to
create plots. For this purpose, we provide the python tool {\tt
pymule}, though of course other tools can be used as well. Here, we
will only cover the aspects of {\tt pymule} required for the present
analysis as shown in Listing~\ref{lst:pymule}; a full documentation
can be found in {\tt docstrings} used in {\tt pymule} as well as
Appendix~\ref{sec:pymule}. First, we import {\tt pymule}. Next, we
need to point {\tt pymule} to the output directory of {\tt mcmule}
with the {\tt setup} command. In our example this is {\tt
babar-tau-e/out}.
\begin{figure}
\centering
\lstinputlisting[language=python,firstline=10,lastline=51]{figures/mule/babar.py}
\renewcommand{\figurename}{Listing}
\caption{An example code to analyse the results for $\tau\to\nu\bar\nu
e\gamma$ in {\tt pymule}. Note that, in the Fortran code
$G_F=\alpha=1$. In {\tt pymule} they are at their physical
values~\cite{PDG}.}
\label{lst:pymule}
\end{figure}
As a next step, we import the \ac{LO} and \ac{NLO} {\tt
which\_pieces} and combine them using two central {\tt pymule}
commands: {\tt sigma} and {\tt mergefks}. {\tt sigma} takes the {\tt
which\_piece} as an argument and imports matching results, already
merging different random seeds. {\tt mergefks} takes the results of
(multiple) {\tt sigma} invocations, adds results with matching $\xc$
values and combines the result. In the present case, $\sigma_n^{(1)}$
is split into multiple contributions, namely {\tt m2enngV} and {\tt
m2enngC}. This is indicated by the {\tt anyxi} argument.
Users should keep in mind that \mcmule{} ships with a version of {\tt
global\_def} where the couplings $G_F={\tt GF}$ and $\alpha={\tt
alpha}$ are set to $G_F=\alpha=1$. Hence, we use {\tt pymule}'s
function {\tt scaleset} to multiply the result with the correct values
of $G_F$ (in ${\rm MeV}^{-1}$) and $\alpha$ (in the \ac{OS} scheme).
Next, we can use some of {\tt pymule}'s tools (cf.
Listing~\ref{lst:pymule}) to calculate the full \ac{NLO} \ac{BR}s from
the corrections and the \ac{LO} results
\begin{align}
\mathcal{B}|_\text{LO} &= 1.8339(1) \times 10^{-2}\,, &
\mathcal{B}|_\text{NLO} &= 1.6451(1) \times 10^{-2}\,,
\end{align}
which agree with~\cite{Fael:2015gua,Pruna:2017upz}, but
$\mathcal{B}|_\text{NLO}$ is in tension with the value
$\mathcal{B}|_\text{exp} = 1.847(54)\times 10^{-2}$ reported by
{\sc BaBar}~\cite{Lees:2015gea, Oberhof:2015snl}. As discussed
in~\cite{Pruna:2017upz, Ulrich:2017adq} it is very likely that this
tension would be removed if a full \ac{NLO} result was used to take
into account the effects of the stringent experimental cuts to extract
the signal. This issue has been explained in detail
in~\cite{Pruna:2017upz, Ulrich:2017adq, Ulrich:2020phd}.
As a last step, we can use the {\tt matplotlib}-backed {\tt kplot}
command to present the results for the distributions (logarithmic for
$m_{e\gamma}$ and linear for $E_e$). The results are shown in
Figure~\ref{fig:babares}. The upper panel of
Figure~\ref{fig:babares:minv} shows the results for the invariant mass
$m_{e\gamma}$ at \ac{LO} (green) and \ac{NLO} (blue) in the range
$0\le m_{e\gamma} \le 1\,\gev$. Note that this, for the purposes of
the demonstration, does not correspond to the boundaries given in the
run.
The distribution falls sharply for large $m_{e\gamma}$. Consequently,
there are only few events generated in the tail and the statistical
error becomes large. This can be seen clearly in the lower panel,
where the \ac{NLO} $K$ factor is shown. It is defined as
\begin{align}
K^{(1)}
=1+\frac{\D\sigma^{(1)}}{\D\sigma^{(0)}}\,,
\end{align}
and the band represents the statistical error of the Monte Carlo
integration. To obtain a reliable prediction for larger values of
$m_{e\gamma}$, i.e. the tail of the distribution, we would have to
perform tailored runs. To this end, we should introduce a cut
$m_\text{cut}\ll m_\tau$ on $m_{e\gamma}$ to eliminate events with
larger invariant mass. Due to the adaption in the numerical
integration, we then obtain reliable and precise results for values of
$m_{e\gamma} \lesssim m_\text{cut}$.
Figure~\ref{fig:babares:el} shows the electron energy distribution,
again at \ac{LO} (green) and \ac{NLO} (blue). As for $m_{e\gamma}$
the corrections are negative and amount to roughly $10\%$. Since this
plot is linear, they can be clearly seen by comparing \ac{LO} and
\ac{NLO}. In the lower panel once more the $K$ factor is depicted.
Unsurprisingly, at the very end of the distribution, $E_e\sim
900\,{\rm MeV}$, the statistics is out of control.
\begin{figure}
\centering
\subfloat[The invariant mass distribution $m_{e\gamma}$.\label{fig:babares:minv}
]{
\scalebox{0.45}{\input{figures/mule/tau:minv.pgf}}
}
\subfloat[The electron energy distribution
\label{fig:babares:el}
]{
\scalebox{0.45}{\input{figures/mule/tau:energy.pgf}}
}
\caption{Results of the toy run to compute $m_{e\gamma}$ (left) and
$E_e$ (right) for $\tau\to\nu\bar\nu e\gamma$. Upper panels show the
\ac{LO} (green) and \ac{NLO} (blue) results, the lower panels show
the \ac{NLO} $K$ factor.}
\label{fig:babares}
\end{figure}