\chapter{Data analysis and results}
\label{cha:data_analysis}
This chapter presents the first analysis on subsets of the collected data for
the aluminium 100-\si{\um}-thick target. The analysis use information from
silicon, germanium, and upstream muon detectors. Pulse parameters were
extracted from waveforms by the simplest method of peak sensing (as mentioned
in \cref{sub:offline_analyser}).
Purposes of the analysis include:
\begin{itemize}
  \item testing the analysis chain;
  \item verification of the experimental method, specifically the
    normalisation of number of stopped muons, and particle identification
    using specific energy loss;
  \item extracting a preliminary rate and spectrum of proton emission from 
    aluminium.
\end{itemize}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Number of stopped muons normalisation}
\label{sec:number_of_stopped_muons_normalisation}
The active silicon target runs was used to check for the validity of the
counting of number of stopped muons, where the number can be calculated by two
methods:
\begin{itemize}
  \item counting hits on the active target in coincidence with hits on the
    upstream scintillator counter;
  \item inferred from number of X-rays recorded by the germanium detector.
\end{itemize}
This analysis was done on a subset of the active target runs
\numrange{2119}{2140}, which contains \num{6.43E7} muon events.  
%\num[fixed-exponent=2, scientific-notation = fixed]{6.4293720E7} muon events.  

\subsection{Number of stopped muons from active target counting}
\label{sub:event_selection}
Because of the active target, a stopped muon would cause two coincident hits on
the muon counter and the target. The energy of the muon hit on the active
target is also well-defined as the narrow-momentum-spread beam was used. The
correlation between the energy and timing of all the hits on the active target
is shown in \cref{fig:sir2f_Et_corr}.

\begin{figure}[tbp]
  \centering
    \includegraphics[width=0.85\textwidth]{figs/sir2_sir2f_E_t_corr}
    \includegraphics[width=0.85\textwidth]{figs/sir2_sir2f_amp_1us_slices}
    \caption{Energy - timing correlation of hits on the active target (top),
      and the projections onto the energy axis in 1000-\si{\ns}-long slices
      from \SI{1500}{\ns} (bottom). The prompt peak at roughly \SI{5}{\MeV} in
      the top plot is muon peak. In the delayed energy spectra, the Michel
      electrons dominate at early time, then the beam electrons are more
      clearly seen in longer delay.} 
    \label{fig:sir2f_Et_corr}
\end{figure}

The prompt hits on the active silicon detector are mainly beam particles:
muons and electrons. The most intense spot at time zero 
and about \SI{5}{\MeV} energy corresponds to stopped muons in the thick target.
The band below \SI{1}{\MeV} is due to electrons, either in the beam or from
muon decay in orbits, or emitted during the cascading of muon to the muonic 1S
state. The valley between time zero and 1200~ns shows the minimum distance in
time between two pulses. It is the limitation of the current pulse
parameter extraction method where no pile up or double pulses is accounted for. 

The delayed hits on the active target after 1200~ns are mainly secondary
particles from the stopped muons:
\begin{itemize}
  \item electrons from muon decay in the 1S orbit,
  \item products emitted after nuclear muon capture, including: gamma, neutron,
    heavy charged particles and recoiled nucleus. 
\end{itemize}
It can be seen that there is a faint stripe of muons in the time larger than
1200~ns region, they are scattered muons by other materials without hitting the
muon counter. The electrons in the beam caused the constant band below 1 MeV and 
$t > 5000$ ns (see \cref{fig:sir2f_Et_corr} bottom).
%\begin{figure}[htb]
  %\centering
    %\includegraphics[width=0.85\textwidth]{figs/sir2_sir2f_amp_1us_slices}
  %\caption{Energy deposit on the active target in 1000 ns time slices from the
    %muon hit. The peaks at about 800 keV in large delayed time are from
    %the beam electrons.}
  %\label{fig:sir2_1us_slices}
%\end{figure}

From the energy-timing correlation above, the cuts to select stopped muons are:
\begin{enumerate}
  \item has one hit on muon counter (where a threshold was set to reject MIPs),
    and the first hit on the silicon active target is in coincidence with that
    muon counter hit: 
    \begin{equation}
      \lvert t_{\textrm{target}} - t_{\mu\textrm{ counter}}\rvert \le
      \SI{50}{\ns}\,,
      \label{eqn:sir2_prompt_cut}
    \end{equation}
  \item and the first hit on the target has energy of that of the muons: 
    \begin{equation}
      \SI{3.4}{\MeV} \le E_{\textrm{target}} \le  \SI{5.6}{\MeV}\,.
      \label{eqn:sir2_muE_cut}
    \end{equation}
\end{enumerate}
The two cuts~\eqref{eqn:sir2_prompt_cut} and~\eqref{eqn:sir2_muE_cut} give
a number of stopped muons counted by the active target:
\begin{equation}
  N_{\mu \textrm{ active Si}} = 9.32 \times 10^6 \pm 3.0\times10^3\,.
  \label{eqn:n_stopped_si_count}
\end{equation}

\subsection{Number of stopped muons from the number of X-rays}
\label{sub:number_of_stopped_muons_from_the_number_of_x_rays}
The number of nuclear captures, hence the number of stopped muons in the
active silicon target, can be inferred from the number of emitted
muonic X-rays. The reference energies and intensities of the most prominent
lines of silicon and aluminium are listed in \cref{tab:mucap_pars}.
\begin{table}[htb]
  \begin{center}
    \begin{tabular}{l l l}
      \toprule
      \textbf{Quantity} & \textbf{Aluminium} & \textbf{Silicon}\\
      \midrule 
      Muonic mean lifetime (ns) & $864 \pm 2$ & $758 \pm 2$\\
      Nuclear capture probability (\%) & $60.9 $ & $65.8$\\
      $(2p-1s)$ X-ray energy (keV) & $346.828\pm0.002$ & $400.177\pm0.005$\\
      Intensity (\%) & $79.8\pm0.8$ & $80.3\pm0.8$\\
      \bottomrule
    \end{tabular}
  \end{center}
  \caption{Reference parameters of muon capture in aluminium and silicon taken
  from Suzuki et al.~\cite{SuzukiMeasday.etal.1987} and Measday et
  al.~\cite{MeasdayStocki.etal.2007}.}
  \label{tab:mucap_pars}
\end{table}

The muonic X-rays are emitted during the cascading of the muon to the muonic 1S
state in the time scale of \SI{E-9}{\s}, so the hit caused by the X-rays must
be in coincidence with the muon hit on the active target. Therefore an
additional timing cut is applied for the germanium detector hits: 
\begin{equation}
  \lvert t_{\textrm{Ge}} - t_{\mu\textrm{ counter}} \rvert < \SI{500}{\ns}
  \label{eqn:sir2_ge_cut}
\end{equation}
\begin{figure}[!htb]
  \centering
    \includegraphics[width=0.85\textwidth]{figs/sir2_xray_22}
    \caption{Prompt muonic X-rays spectrum from the active silicon target. The
      $(2p-1s)$ X-ray shows up at \SI{400}{\keV}; higher transitions can also
      be identified.}
    \label{fig:sir2_xray}
\end{figure}

The germanium spectrum after three
cuts~\eqref{eqn:sir2_prompt_cut},~\eqref{eqn:sir2_muE_cut}
and~\eqref{eqn:sir2_ge_cut} is plotted in \cref{fig:sir2_xray}. The $(2p-1s)$
line clearly showed up at \SI{400}{\keV} on a very low background. A peak at
\SI{476}{\keV} is identified as the $(3p-1s)$ transition. Higher transitions
such as $(4p-1s)$, $(5p-1s)$ and $(6p-1s)$ can also be recognised at
\SI{504}{\keV}, \SI{516}{\keV} and \SI{523}{\keV}, respectively.
%The $(2p-1s)$
%line is seen at 399.5~\si{\keV}, 0.7~\si{\keV} off from the reference value of
%400.177~\si{\keV}. This discrepancy is within our detector's resolution,
%and the calculated efficiency is $(4.549 \pm 0.108)\times 10^{-5}$ -- a 0.15\%
%increasing from that of the 400.177~keV line, so no attempt for recalibration
%or correction was made.  

The net area of the $(2p-1s)$ is found to be 2929.7 by fitting a Gaussian
peak on top of a linear background from \SIrange{395}{405}{\keV}. 
Using the same procedure of correcting described in
\cref{sub:germanium_detector}, and taking detector acceptance and X-ray
intensity into account (see \cref{tab:sir2_xray_corr}), the number of muon
stopped is: 
\begin{equation}
  N_{\mu \textrm{ stopped X-ray}} = (9.16 \pm 0.28)\times 10^6\,,
  \label{eqn:n_stopped_xray_count}
\end{equation}
which is consistent with the number of X-rays counted using the active target.

The uncertainty of the number of muons inferred from the X-ray
has equal contributions from statistical uncertainty in peak
area and systematic uncertainty from efficiency calibration. The relative
uncertainty in number of muons is 3\%, good enough for the normalisation in
this measurement.
\begin{table}[btp]
  \begin{center}
    \begin{tabular}{@{}llll@{}}
      \toprule
      \textbf{Measured X-rays}          & \textbf{Value} & \textbf{Absolute error}                                    & \textbf{Relative error}                                   \\ \midrule
      Gross integral                    & 3083           &                                                            &                                                           \\
      Background                        & 101.5          &                                                            &                                                           \\
      Net area $(2p-1s)$                & 2929.7         & 56.4                                                       & 0.02                                                      \\ 
      \vspace{0.03em}\\
      \toprule
      \textbf{Corrections} & \textbf{Value} & \multicolumn{2}{c}{\textbf{Details}}\\
      \midrule
      Random summing                    & 1.06  & \multicolumn{2}{l}{average count rate \SI{491.4}{\Hz},}\\ 
                                        &       & \multicolumn{2}{l}{pulse length \SI{57}{\us}}\\
      TRP reset                         & 1.03 & \multicolumn{2}{l}{\SI{298}{\s} loss during \SI{9327}{\s} run period}\\
      Self-absorption                   & 1.008 & \multicolumn{2}{l}{silicon thickness \SI{750}{\um},}\\ 
                                        &       & \multicolumn{2}{l}{linear attenuation \SI{0.224}{\per\cm}}\\
      True coincidence                  & 1              & \multicolumn{2}{l}{omitted}                                                                                            \\ 
      \vspace{0.03em}\\
      \toprule
      \textbf{Efficiency and intensity} & \textbf{Value} & \textbf{Absolute error}                                    & \textbf{Relative error}                                   \\
      \midrule
      Detector efficiency               & \num{4.40E-4}        & \num{0.10E-4}                                                   & 0.02                                                     \\
      X-ray intensity                   & 0.803          & 0.008                                                      & 0.009                                                     \\ 
      \vspace{0.03em}\\
      \toprule
      \textbf{Results}                  & \textbf{Value} & \textbf{Absolute error}                                    & \textbf{Relative error}                                   \\
      \midrule
      Number of X-rays emitted          & \num{7.36E6}         & \num{0.22E6}                                                     & 0.03                                                      \\
      Number of muons stopped           & \num{9.16E6}         & \num{0.28E6}                                                     & 0.03                                                      \\
      \bottomrule
    \end{tabular}
  \end{center}
  \caption{Corrections, efficiency and intensity used in calculating the number
    of X-rays from the active target.} 
  \label{tab:sir2_xray_corr}
\end{table}

%In order to measure the charged particles after nuclear muon capture, one would
%pick events where the emitted particles are well separated from the
%muon stop time. The energy timing correlation plot suggests a timing window
%starting from at least 1200~ns, therefore another cut is introduced:
%\begin{enumerate}
    %\setcounter{enumi}{2}
  %\item there are at least two hits on the active target, the time
    %difference between the second hit on target (decay or capture product) and
    %the muon counter hit is at least 1300 ns:
    %\begin{equation}
      %t_{\textrm{target 2nd hit}} - t_{\mu\textrm{ counter}} \geq \SI{1300}{\ns}
      %\label{eqn:sir2_2ndhit_cut}
    %\end{equation}
%\end{enumerate}

%The three cuts~\eqref{eqn:sir2_prompt_cut},~\eqref{eqn:sir2_muE_cut} and
%~\eqref{eqn:sir2_2ndhit_cut} reduce the sample to the size of \num{9.82E+5}.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Particle identification by specific energy loss}
\label{sec:particle_identification_by_specific_energy_loss}
In this analysis, a subset of runs from \numrange{2808}{2873} with  the
100-\si{\um} aluminium target is used because of following advantages:
\begin{itemize}
  \item it was easier to stop and adjust the muon stopping distribution in
    this thicker target;
  \item a thicker target gives better statistics because of a higher
    muon rate available at a higher momentum and less scattering.
\end{itemize}
Muons with momentum of \SI{30.52}{\MeV\per\hepclight}, 3\%-FWHM spread (scaling factor of
1.09, normalised to \SI{28}{\MeV\per\hepclight}) were used for this target after
a momentum scanning as described in the next subsection.

\subsection{Momentum scan for the 100-\si{\um} aluminium target}
\label{sub:momentum_scan_for_the_100_}
Before deciding to use the momentum scaling factor of 1.09, we have scanned
with momentum scales ranging from 1.04 to 1.12 to maximise the 
observed X-rays rate (and maximising the rate of stopped muons). The X-ray
spectrum at each momentum point was accumulated in about 30 minutes to 
assure a sufficient amount of counts. Details of the scanning runs are listed
in \cref{tab:al100_scan}.
\begin{table}[htb]
  \begin{center}
    \begin{tabular}{cccc}
      \toprule
      \textbf{Momentum (\si{\MeV\per\hepclight})} & \textbf{Scaling factor} & \textbf{Runs}
      & \textbf{Length (s)}\\
      \midrule 
      29.12 & 1.04 & \numrange{2609}{2613} &2299\\
      29.68 & 1.06 & \numrange{2602}{2608} &2563\\
      29.96 & 1.07 & \numrange{2633}{2637} &2030\\
      30.24 & 1.08 & \numrange{2614}{2621} &3232\\
      30.52 & 1.09 & \numrange{2808}{2813} &2120\\
      30.80 & 1.10 & \numrange{2625}{2632} &3234\\
      31.36 & 1.12 & \numrange{2784}{2792} &2841\\
      \bottomrule
    \end{tabular}
  \end{center}
  \caption{Momentum scanning runs for the 100-\si{\um} aluminium target.}
  \label{tab:al100_scan}
\end{table}
The on-site quick analysis suggested the 1.09 scaling factor was the
optimal value so it was chosen for all the runs on this aluminium target. But
the offline analysis later showed that the actual optimal factor was 1.08.
There were two reasons for the discrepancy:
\begin{enumerate}
  \item the X-ray rates were normalised to run length, which is biased
    since there are more muons available at higher momenta;
  \item the $(2p-1s)$ peaks of aluminium at \SI{346.828}{\keV} were not
    fitted properly. The peak is interfered by a background peak at
    \SI{351}{\keV} from $^{214}$Pb, but the X-ray peak area was
    obtained simply by subtracting an automatically estimated background.
\end{enumerate}
In the offline analysis, the X-ray peak and the background peak are fitted by
two Gaussian peaks on top of a linear background. The X-ray peak
area is then normalised to the number of muons hitting the upstream detector
(\cref{fig:al100_xray_fit}).
\begin{figure}[htb]
  \centering
    \includegraphics[width=0.50\textwidth]{figs/al100_xray_fit}
    \includegraphics[width=0.50\textwidth]{figs/al100_xray_musc}
  \caption{Fitting of the $(2p-1s)$ muonic X-ray of aluminium (red) and the
    interfered peak at \SI{351}{\keV} (brown) with a linear background (left).
    The number of muons is integral of the upstream scintillator spectrum
    (right) from \numrange{400}{2000} ADC channels.} 
  \label{fig:al100_xray_fit}
\end{figure}

The ratio between the number of X-rays and the number of muons as a function
of momentum scaling factor is plotted on \cref{fig:al100_scan_rate}. The trend
showed that muons penetrated deeper as the momentum increased, reaching the
optimal value at the scale of 1.08, then decreased as punch through happened
more often from scales of 1.09 and above. The distributions of stopped muons
are illustrated by MC results on the bottom plot in
\cref{fig:al100_scan_rate}. At the 1.09
scale beam, the muons stopped \SI{18}{\um} off-centred to the right silicon
arm, the standard deviation of the depth distribution is \SI{29}{\um}.
\begin{figure}[!htb]
  \centering
    \includegraphics[width=0.77\textwidth]{figs/al100_scan_rate}
    \includegraphics[width=0.77\textwidth]{figs/al100_mu_stop_mc}
  \caption{Number of X-rays per incoming muon as a function of momentum
    scaling factor (top); and muon stopping distributions with scaling factors
    from 1.04 to 1.12 obtained by MC simulation
    (bottom). The depth of muon stopping positions are measured normal to
    the surface of the target facing the muon beam.}
  \label{fig:al100_scan_rate}
\end{figure}

\subsection{Event selection for the passive targets}
\label{sub:event_selection_for_the_passive_targets}
As described in the \cref{sec:analysis_framework}, the hits on all detectors
are re-organised into muon events: central muons; and all hits within
\SI{\pm 10}{\us} from the central muons. The dataset from runs
\numrange{2808}{2873} contains \num{1.17E+9} of such muon events.

\subsubsection{Particle banding identification}
\label{ssub:particle_banding_identification}

Selection of proton (and other heavy charged particles) events starts from
searching for muon event that has at least one hit on thick silicon. If there
is a thin silicon hit within a coincidence window of $\pm 0.5$~\si{\us}\ around
the thick silicon hit, the two hits are considered to belong to one particle.
The thresholds for energy deposited in all silicon channels, except the thin
silicon on the left arm, are set at \SI{100}{\keV} in this analysis. The
threshold on the left $\Delta E$ counter was higher, at roughly
\SI{400}{\keV}, due to higher noise in that channel and it was decided at the
run time to rise its threshold to reduce the triggering rate.
The specific energy loss as a function of total energy of the charged
particles are plotted on \cref{fig:al100_dedx}.
\begin{figure}[p]
  \centering
  \includegraphics[width=0.85\textwidth]{figs/al100_EdE_left}
  \includegraphics[width=0.85\textwidth]{figs/al100_EdE_right}
  \caption{Energy loss in thin silicon detectors as a function of total energy
    recorded by both thin and thick detectors on the left arm (top) and the
    right arm (bottom).}
  \label{fig:al100_dedx}
\end{figure}
With the aid from MC simulation (\cref{fig:pid_sim}), the banding on
\cref{fig:al100_dedx} can be identified as follows: 
\begin{itemize}
  \item the spot at the lower left conner belonged to electron hits; 
  \item the scattered muons formed the small blurry cloud just above the
    electron region;
  \item the most intense band was due to proton hits; 
  \item the less intense, upper band caused by deuteron hits; 
  \item the highest band corresponded to alpha hits; 
  \item the faint stripe above the deuteron band should be triton
    hits, which is consistent with a relatively low probability of emission of
    tritons.  
\end{itemize}
\begin{figure}[htb]
  \centering
    \includegraphics[width=0.85\textwidth]{figs/al100_dedx_overlay}
    \caption{Identifying of charged particles banding: the dots are measured
      points, the histograms are expected bands of protons (red), deuterons
      (green) and tritons (blue), respectively. The MC bands are calculated
      for a pair of 58-\si{\um}-thick and 1535-\si{\um}-thick silicon
      detectors. The error bars on MC bands show the standard deviation of
      $\Delta E$ in E respective bins.
    }
    \label{fig:dummylabel}
\end{figure}

It is not clearly seen in the $\Delta E-E$ plots because of the rather high
thresholds on $\Delta E$, but protons with higher energy would punch through
both silicon detectors. Those events have low $\Delta E$ and $E$, making the
proton bands to go backward to the origin of the $\Delta E-E$ plots. For the
configuration of 58-\si{\um} thin, and 1535-\si{\um} thick detectors, the
effect shows up for protons with energy larger than \SI{16}{\MeV}. The
returning part of the proton band would make the cut described in
the next subsection to include protons with higher energy into lower energy
region. The effect of punch through protons could be eliminate using the veto
plastic scintillators at the back of each silicon arm. But in this initial
analysis, the veto information is not used, therefore the upper limit of
proton energy is set at \SI{8}{\MeV} where there is clear separation between
the protons at lower and higher energies with the same measured total energy
deposition $E$.

\subsubsection{Proton-like probability cut}
\label{ssub:proton_like_probability_cut}
Since protons of interested are at low kinetic energy, their $\Delta E$
distributions do not have long tails as that of the Landau distribution.
For a given $E$, the distribution of $\Delta E$ is more like a Gaussian, and
with slightly deformed high energy tail (see \cref{fig:dE_distribution}).

\begin{figure}[htb]
  \centering
  \includegraphics[width=0.75\textwidth]{figs/dE_distribution}
  \caption{Distributions of $\Delta E$ of protons in a 58-\si{\um}-thick
    silicon detector for given $E$ in various energy ranges.} 
  \label{fig:dE_distribution}
\end{figure}

%In order to select protons, a proton likelihood probability is defined as:
%The band of protons is therefore by cut on likelihood probability
%calculated as:
For a measured event, a proton likelihood probability is defined as:
\begin{equation}
  P_{i} = \dfrac{1}{\sqrt{2\pi}\sigma_{\Delta E}}
  \exp{\left[\dfrac{(\Delta E_{meas.} - \Delta E_i)^2} {2\sigma^2_{\Delta
          E}}\right]}\,,
\end{equation}
where $\Delta E_{\textrm{meas.}}$ and $E_i$ are measured energy deposition in
thin silicon detector and in both detectors, respectively; $\Delta E_i$  and
$\sigma_{\Delta E}$ are the expected value and standard deviation of the energy
loss in the thin detector of protons with energy $E$, calculated by the
MC simulation.  A measured event with higher $P_i$ is more likely to be
a proton event.

The lower threshold of proton-like probability, the more protons will be
selected, but also more contamination from other charged particles would be
classified as protons. The number of protons on the left and right arms at
different cuts on $P_i$ are listed in \cref{tab:nprotons_vs_pcut}. The proton
yields are integrated in the energy range from \SIrange{2.2}{8}{\MeV}. The
lower limit comes from the requirement of having at least one hit on the thick
counter. The upper limit is to avoid the inclusion of punch through particles
as explained in the previous session. 

The cut efficiency depends on actual shape of the proton spectrum, other
charged particles spectra, relative ratio between the yields of different
particle species. The fraction of protons missed out, and the fraction of
contamination from other charged particles at different probability
thresholds, with two different assumptions on spectrum shape: flat distribution
and Sobottka and Wills silicon shape (see \eqref{eqn:EH_pdf}),
are listed in the four last columns of \cref{tab:nprotons_vs_pcut}. The
relative ratio of proton:deuteron:triton:alpha:muon is assumed to be
5:2:1:2:2. The probability threshold is therefore chosen to be \num{1.0E-4} in
in order to have both relatively low missing and contamination fractions
compare to the statistical uncertainty of the measurement. The resulted band of
protons is shown in (\cref{fig:al100_protons}).  

\begin{table}[htb]
  \begin{center}
    \begin{tabular}{c c c c S S S S}
      \toprule
      \textbf{$P_i$} & \textbf{Equiv.} & \multirow{2}{*}{\textbf{Left}} & 
      \multirow{2}{*}{\textbf{Right}} & {\textbf{Miss.}} & \textbf{Contam.}
      & {\textbf{Miss.}} & \textbf{Contam.}\\

      \textbf{threshold} & \textbf{$\sigma$} & & 
      &{\textbf{flat}} & {\textbf{flat}}
      &{\textbf{expo.}} & {\textbf{expo.}}\\
      \midrule 
      \num{4.5E-2} & 2 & 1720 & 2214 &    1.9  \%& 0.03 \%&6.1  \%& 0.06 \%\\
      \num{2.7E-3} & 3 & 1801 & 2340 &    0.7  \%& 0.05 \%&2.8  \%& 0.1 \%\\
      \num{1.0E-4} & 3.89 & 1822 & 2373 & 0.5  \%& 0.1  \%&1.2  \%& 0.3  \%\\
      \num{5.7E-7} & 5 & 1867 & 2421 &    0.4  \%& 0.7  \%&0.7  \%& 0.9  \%\\
      \bottomrule
    \end{tabular}
  \end{center}
  \caption{Proton yields in energy range from \SIrange{2.2}{8}{\MeV} on the two
    silicon arms with different thresholds on proton-like probability $P_i$,
    and the MC calculated missing fractions and contamination levels with two
    different assumptions on spectrum shape: flatly distributed, and
    exponential decay spectrum (see \eqref{eqn:EH_pdf}).}
  \label{tab:nprotons_vs_pcut}
\end{table}

\begin{figure}[htb]
  \centering
    \includegraphics[width=0.47\textwidth]{figs/al100_protons}
    \includegraphics[width=0.47\textwidth]{figs/al100_protons_px_r}
  \caption{Protons (green) selected using the likelihood probability cut of
    \num{1.0E-4} (left). The proton spectrum (right) is obtained by projecting
    the proton band onto the total energy axis.}
  \label{fig:al100_protons}
\end{figure}

\subsubsection{Possible backgrounds}
\label{ssub:possible_backgrounds}

There are several sources of potential backgrounds in this proton measurement:
\begin{enumerate}
  \item Protons emitted after capture of scattered muons in the lead
    shield: the incoming muons could be scattered to other materials
    surrounding the target, emitting protons to the silicon detectors. In
    order to avoid complication of estimating this background, we used lead
    sheets to collimate and shield around the target and detectors. If
    a scattered muon is captured by the lead shielding, the proton from lead
    would be emitted shortly after the muon hit because of the short average
    lifetime of muons in lead (\SI{78.4}{\ns}~\cite{Measday.2001}). In
    comparison, average lifetime of muons in aluminium is
    \SI{864}{\ns}~\cite{Measday.2001}, therefore a simple cut in timing could
    eliminate background of this kind.\\
   The timing of events classified as protons are plotted in
   \cref{fig:al100_proton_timing}. The spectra show no significant fast
   decaying component, which should show up if the background from lead
   shielding were sizeable. A fit of an exponential function on top of a flat
   background gives the average lifetimes of muons as:
   \begin{align}
     \tau_{\textrm{left}} &= \SI{870 \pm 25}{\ns} \,,\\
     \tau_{\textrm{right}} &= \SI{868 \pm 21}{\ns} \,.
   \end{align}
   
  \begin{figure}[htb]
    \centering
    \includegraphics[width=0.85\textwidth]{figs/al100_proton_timing}
    \caption{Timing of protons relative to muon hit. The spectra show the
      characteristic one-component decay shape.}
    \label{fig:al100_proton_timing}
  \end{figure}
  The consistency between fitted lifetimes and the reference value of average
  lifetime of muons in aluminium at \SI{864\pm 2}{\ns} suggests the background
  from the lead shielding is negligible. This smallness can be explained as
  a combination of the two facts: (i) only
  a minority fraction of muons punched through the target and reached the 
  downstream lead shield as illustrated in
  \cref{fig:al100_scan_rate}; and (ii) the probability of emitting protons from
  lead is very low compare to that of aluminium, about 0.4\% per
  capture (see \cref{tab:lifshitzsinger_cal_proton_rate}). 

  \item The protons emitted after scattered muons stopped at the surface of
    the thin silicon detectors: these protons could mimic the signal if they
    appear within \SI{1}{\us} around the time muon hit the upstream counter.
    The $\Delta E$ and $E$ in this case would be sum of energy of a muon and
    energy of the resulted proton. The average energy of scattered muons can be
    seen in \cref{fig:al100_dedx} to be about \SI{1.4}{\MeV}. The measured
    $\Delta E$ and $E$ then would be shifted by \SI{1.4}{\MeV}, makes the
    measured data point move far away from the expected proton band. Therefore
    this kind of background should be small with the current probability cut.

  \item The random background: this kind of background can be 
    examined by the same timing spectrum in
    \cref{fig:al100_proton_timing}. The random events show up at negative time
    difference and large delay time regions and give a negligible contribution
    to the total number of protons. 
\end{enumerate}

It is concluded from above arguments that the backgrounds of this proton
measurement is negligibly small.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Proton emission rate from aluminium}
\label{sec:proton_emission_rate_from_aluminium}
The analysis is done on the same dataset used in
\cref{sec:particle_identification_by_specific_energy_loss}. Firstly, the
number of protons emitted is extracted using specific energy loss. Then
correction for energy loss inside the target is applied. Finally, the number
of protons is normalised to the number of nuclear muon captures.

\subsection{Number of protons emitted}
\label{sub:number_of_protons_emitted}
The numbers of protons in the energy range from \SIrange{2.2}{8.5}{\MeV} after
applying the probability cut are: 
\begin{align}
  N_{\textrm{p meas. left}} = 1822 \pm 42.7 \,,\\
  N_{\textrm{p meas. right}} = 2373 \pm 48.7 \,.
\end{align}
The right arm received significantly more protons than the left arm did, which
is expected as in \cref{sub:momentum_scan_for_the_100_} where it is shown that
muons stopped off-centred to the right arm.

\subsection{Corrections for the number of protons}
\label{sub:corrections_for_the_number_of_protons}
The protons spectra observed by the silicon detectors have been modified by
the energy loss inside the target so correction (or unfolding) is necessary.
The unfolding, essentially, is finding a response function that relates proton's
true energy and measured value. This can be done in MC simulation by generating
protons with a spatial distribution as close as possible to the real
distribution of muons inside the target, then counting the number of protons
reaching the silicon detectors. Such response function conveniently includes
the geometrical acceptance.

For the 100-\si{\um} aluminium target and muons at the momentum scale of 1.09,
the parameters of the initial protons are:
\begin{itemize}
  \item horizontal distribution: Gaussian \SI{26}{\mm} FWHM, centred at the
    centre of the target;
  \item vertical distribution: Gaussian \SI{15}{\mm} FWHM, centred at the
    centre of the target;
  \item depth: Gaussian \SI{69.2}{\um} FWHM, centred at \SI{68.8}{\um}-deep from
    the upstream face of the target;
  \item energy: flatly distributed from \SIrange{1.5}{15}{\MeV}.
\end{itemize}
The calculated response matrices for the two arms are presented in
\cref{fig:al100_resp_matrices}. The different path lengths inside the target
to the two silicon arms causes the difference in the two matrices. The
response matrices are then used as MC truth to train  and test the unfolding
code. The code uses an existing ROOT package called RooUnfold~\cite{Adye.2011}
where the iterative Bayesian unfolding method is implemented.
\begin{figure}[htb]
  \centering
  \includegraphics[width=0.99\textwidth]{./figs/al100_resp}
  \caption{Response functions for the two silicon arms, showing the relation
    between protons energy at birth and as detected by the silicon detector
    arms.}
  \label{fig:al100_resp_matrices}
\end{figure}
%After training, the unfolding code is applied on the measured spectra from the
%left and right arms. The unfolded proton spectra in \cref{fig:al100_unfold}
%reasonably reflect the distribution of initial protons which is off-centred to
%the right arm. The path length to the left arm is longer so less protons at
%energy lower than \SI{5}{\MeV} could reach the detectors. The sharp low-energy
%cut off on the right arm is consistent with the Coulomb barrier for protons,
%which is \SI{4.1}{\MeV} for protons emitted from $^{27}$Mg. 
The unfolded spectra using the two observed spectra at the two arms as input
are shown in \cref{fig:al100_unfold}. The two unfolded spectra generally agree
with each other, except for a few first and last bins. 
In the lower energy region, there is a small probability for such protons to
escape and reach the detectors, therefore the unfolding is generally unstable
and the uncertainties are large.
At the higher end, the jump on the right arm at around \SI{9}{\MeV} can be
explained as the punch-through protons were counted as the proton veto counters
were not used in this analysis. The lower threshold on the thin silicon
detector at the right arm compared with that at the left arm makes this
misidentification worse.
\begin{figure}[!htb]
  \centering
  \includegraphics[width=0.80\textwidth]{figs/al100_unfolded_lr}
  \caption{Unfolded proton spectra from the 100-\si{\um} aluminium target.}
  \label{fig:al100_unfold}
\end{figure}

%Several studies were conducted to assess the performance of the unfolding
%code, including:
%\begin{itemize}
  %\item stability against cut-off energy;
  %\item comparison between the two arms;
  %\item and unfolding of a MC-generated spectrum.
%\end{itemize}
The stability of the unfolding code is tested by varying the lower and upper
cut-off energies of the input spectrum. Plots in \cref{fig:al100_cutoff_study}
show that the shapes of the unfolded spectra are stable after a few first or
last bins.
%The
%lower cut-off energy of the
%output increases as that of the input increases, and the shape is generally
%unchanged after a few bins.
\begin{figure}[!htb]
  \centering
  \includegraphics[width=0.85\textwidth]{figs/al100_cutoff_study}
  \includegraphics[width=0.85\textwidth]{figs/al100_up_cut_off_reco}
  \caption{Unfolded spectra with different lower (top) and upper (bottom)
    cut-off energies.} 
  \label{fig:al100_cutoff_study}
\end{figure}

The proton yields calculated from observed spectra in two arms are compared in
\cref{fig:al100_integral_comparison} where the upper limit of the integrals
is fixed at \SI{8}{\MeV}, and the lower limit is varied in \SI{400}{\keV}
step. The upper limit was chosen to avoid the effects of punched through
protons. The difference is large at cut-off energies less than \SI{4}{\MeV}
due to large uncertainties at the first bins. Above \SI{4}{\MeV}, the two arms
show consistent numbers of protons.
\begin{figure}[htb]
  \centering
  \includegraphics[width=0.85\textwidth]{figs/al100_integral_comparison}
  \caption{Proton yields calculated from two arms. The upper limit of
    integrations is fixed at \SI{8}{\MeV}, the horizontal axis is the lower
    limit of the integrations. The proton yields on the two arm agree well
    with each other from above \SI{4}{\MeV}.}
  \label{fig:al100_integral_comparison}
\end{figure}
The yields of protons from \SIrange{4}{8}{\MeV} are:
\begin{align}
  N_{\textrm{p unfold left}} &= (165.4 \pm 2.7)\times 10^3\\
  N_{\textrm{p unfold right}} &= (173.1 \pm 2.9)\times 10^3
\end{align}
The number of emitted protons is taken as average of the two yields:
\begin{equation}
  N_{\textrm{p unfold}} = (169.3 \pm 1.9) \times 10^3
\end{equation}

\subsection{Number of nuclear captures}
\label{sub:number_of_nuclear_captures}
%\begin{figure}[!htb]
  %\centering
  %\includegraphics[width=0.85\textwidth]{figs/al100_ge_spec}
  %\caption{X-ray spectrum from the aluminium target, the characteristic
  %$(2p-1s)$ line shows up at 346.67~keV}
%\label{fig:al100_ge_spec}
%\end{figure}
 
%The X-ray spectrum on the germanium detector is shown on
%\cref{fig:al100_ge_spec}. 
Fitting the double peaks on top of a linear background
gives the X-ray peak area of $5903.5 \pm 109.2$. With the same
procedure as in the case of the active target, the number stopped muons and
the number of nuclear captures are:
\begin{align}
  N_{\mu \textrm{ stopped}} &= (1.57 \pm 0.05)\times 10^7\,,\\
  N_{\mu \textrm{ nucl. cap.}} &= (9.57\pm 0.31)\times 10^6\,.
\end{align}

\subsection{Proton emission rate}
\label{sub:proton_emission_rate}
The proton emission rate in the range from \SIrange{4}{8}{\MeV} is therefore:
\begin{equation}
  R_{\textrm{p}} = \frac{169.3\times 10^3}{9.57\times 10^6} = 1.7\times
  10^{-2}\,.
  \label{eq:proton_rate_al}
\end{equation}

The total proton emission rate can be estimated by assuming a spectrum shape
with the same parameterisation as in \eqref{eqn:EH_pdf}. The
\eqref{eqn:EH_pdf} function has a power rising edge, and a exponential decay
falling edge. The falling edge has only one decay component and is suitable to
describe the proton spectrum with the equilibrium emission only assumption.
The pre-equilibrium emission contribution should be small for low-$Z$ material,
for aluminium the contribution of this component is 2.2\% of total number of
protons according to Lifshitz and Singer~\cite{LifshitzSinger.1980}.
%%TODO: draw the function and integral
The fitted results
are shown in \cref{fig:al100_parameterisation} and \cref{tab:al100_fit_pars}.
The average spectrum is obtained by taking the average of the two unfolded
spectra from the left and right arms. The fitted parameters are compatible
with each other within their errors.

Using the fitted parameters of the average spectrum, the integration in range
from \SIrange{4}{8}{\MeV} is 51\% of the total number of
protons. The total proton emission rate is therefore estimated to be $3.5\times 10^{-2}$.
\begin{figure}[!p]
  \centering
  \includegraphics[width=0.85\textwidth]{figs/al100_parameterisation}
  \includegraphics[width=0.85\textwidth]{figs/al100_fit_avgspec}
  \includegraphics[width=0.85\textwidth]{figs/al100_fitted_func_integral}
  \caption{Fitting of the unfolded spectra on the left and right arms (top),
    and on the average spectrum (middle). The bottom plot shows the fitted
    function of the average spectrum in the energy range from
    \SIrange{1}{50}{\MeV}. The proton yield in the region from
    \SIrange{4}{8}{\MeV} (shaded) is 51\% of the whole spectral integral.}
  \label{fig:al100_parameterisation}
\end{figure}

\begin{table}[htb]
  \begin{center}
    \begin{tabular}{l S[separate-uncertainty=true] S[separate-uncertainty=true] 
        S[separate-uncertainty=true]}
      \toprule
      \textbf{Parameter} &{\textbf{Left}} & {\textbf{Right}} & {\textbf{Average}}\\
      \midrule 
      $A \times 10^{-5}$ & 2.0 \pm 0.7 & 1.3 \pm 0.1 & 1.5 \pm 0.3\\
      $T_{th}$ (\si{\keV}) & 1301 \pm 490 & 1966 \pm 68 & 1573 \pm 132\\
      $\alpha$             & 3.2 \pm 0.7 & 1.2 \pm 0.1 & 2.0 \pm 1.2\\
      $T_{0}$ (\si{\keV}) & 2469 \pm 203 & 2641 \pm 106 & 2601 \pm 186\\
      \bottomrule
    \end{tabular}
  \end{center}
  \caption{Parameters of the fits on the unfolded spectra, the average spectrum
    is obtained by taking average of the unfolded spectra from left and right
    arms.} 
  \label{tab:al100_fit_pars}
\end{table}


\subsection{Uncertainties of the emission rate}
\label{sub:uncertainties_of_the_emission_rate}
The uncertainties of the emission rate come from:
\begin{itemize}
  \item uncertainties in the number of nuclear captures: these were discussed
    in \cref{sub:number_of_stopped_muons_from_the_number_of_x_rays};
  \item uncertainties in the number of protons:
    \begin{itemize}
      \item statistical uncertainties of the measured spectra which are
        propagated during the unfolding process;
      \item systematic uncertainties due to misidentification: this number is
        small compared to other uncertainties as discussed in
        \cref{sub:event_selection_for_the_passive_targets};
      \item systematic uncertainty from the unfolding
    \end{itemize}
\end{itemize}
The last item is studied by MC method using the parameterisation in
\cref{sub:proton_emission_rate}:
\begin{itemize}
  \item protons with energy distribution obeying the parameterisation are
    generated inside the target. The spatial distribution is the same as that
    of in \cref{sub:corrections_for_the_number_of_protons}. MC truth including
    initial energies and positions are recorded;
  \item the number of protons reaching the silicon detectors are counted,
    the proton yield is set to be the same as the measured yield to make the
    statistical uncertainties comparable;
  \item the unfolding is applied on the observed proton spectra, and the
    results are compared with the MC truth.
\end{itemize}
\begin{figure}[htb]
  \centering
    \includegraphics[width=0.48\textwidth]{figs/al100_MCvsUnfold}
    \includegraphics[width=0.48\textwidth]{figs/al100_unfold_truth_ratio}
  \caption{Comparison between an unfolded spectrum and MC truth. On the left
    hand side, the solid, red line is MC truth, the blue histogram is the
    unfoldede spectrum. The ratio between the two yields is compared in the
    right hand side plot with the upper end of integration is fixed at
    \SI{8}{\MeV}, and a moving lower end of integration. The discrepancy
    is genenerally smaller than 5\% if the lower end energy is smaller than
    \SI{6}{\MeV}.}
  \label{fig:al100_MCvsUnfold}
\end{figure}
\Cref{fig:al100_MCvsUnfold} shows that the yield obtained after unfolding is
in agreement with that from the MC truth. The difference is less than 5\% if
the integration is taken in the range from \SIrange{4}{8}{\MeV}. Therefore
a systematic uncertainty of 5\% is assigned for the unfolding routine.

A summary of uncertainties in the measurement of proton emission rate is
presented in \cref{tab:al100_uncertainties_all}.
\begin{table}[htb]
  \begin{center}
    \begin{tabular}{@{}ll@{}}
      \toprule
      \textbf{Item}& \textbf{Value}          \\ 
      \midrule
      Number of muons                   & 3.2\%          \\
      Statistical from measured spectra & 1.1\%          \\
      Systematic from unfolding         & 5.0\%            \\
      Systematic from PID               & \textless1.0\% \\ 
      \midrule
      Total                             & 6.1\%\\              
      \bottomrule
    \end{tabular}
  \end{center}
  \caption{Uncertainties of the proton emission rate.}
  \label{tab:al100_uncertainties_all}
\end{table}

\section{Results of the initial analysis}
\label{sec:results_of_the_initial_analysis}
\subsection{Verification of the experimental method}
\label{sub:verification_of_the_experimental_method}
The experimental method described in \cref{sub:experimental_method} has been
validated:
\begin{itemize}
  \item Number of muon capture normalisation: the number of stopped muons
    calculated from the muonic X-ray spectrum is shown to be consistent with
    that calculated from the active target spectrum. 
  \item Particle identification: the particle identification by specific
    energy loss has been demonstrated. The banding of different particle
    species is clearly visible. The proton extraction method using cut on
    likelihood probability has been established. Since the distribution of
    $\Delta E$ at a given $E$ is not Gaussian, the fraction of protons that do
    not make the cut is 0.5\%, much larger than the threshold at \num{1E-4}.
    The fraction of other charged particles being misidentified as protons is
    smaller than 0.1\%. These uncertainties from particle identification are
    still small in compared with the
    uncertainty of the measurement (2.3\%).
  \item Unfolding of the proton spectrum: the unfolded spectra inferred from
    two measurements at the two silicon arms show good agreement with each
    other, and with the muon stopping distribution obtained in the momentum
    scanning analysis.
\end{itemize}

\subsection{Proton emission rates and spectrum}
\label{sub:proton_emission_rates_and_spectrum}
The proton emission spectrum in \cref{sub:proton_emission_rate} peaks around
\SI{3.7}{\MeV} which is a little below the Coulomb barrier for proton of
\SI{3.9}{\MeV} calculated using \eqref{eqn:classical_coulomb_barrier}. The
spectrum has a decay constant of \SI{2.6}{\MeV} in higher energy region. 

The partial emission rate measured in the energy range from
\SIrange{4}{8}{\MeV} is: 
\begin{equation}
  R_{p \textrm{ meas. }} = (1.7 \pm 0.1)\%. 
  \label{eqn:meas_partial_rate}
\end{equation}

The total emission rate from aluminium assuming the spectrum shape holds for
all energy is:
\begin{equation}
  R_{p \textrm{ total}} = (3.5 \pm 0.2)\%. 
  \label{eqn:meas_total_rate}
\end{equation}

\subsubsection{Comparison to theoretical and other experimental results}
\label{ssub:comparison_to_theoretical_and_other_experimental_results}
There is no existing experimental or theoretical work that could be directly
compared with the obtained proton emission rate.  Indirectly, it is compatible
with the figures calculated by Lifshitz and
Singer~\cite{LifshitzSinger.1978, LifshitzSinger.1980} listed in
\cref{tab:lifshitzsinger_cal_proton_rate}.  It is significantly larger than
the rate of 0.97\% for the $(\mu,\nu p)$ channel, and does not
exceed the inclusion rate for all channels $\Sigma(\mu,\nu p(xn))$ at 4\%,
leaving some room for other modes such as emission of deuterons or tritons.
Certainly, when the full analysis is available, deuterons and tritons emission
rates could be extracted and the combined emission rate could be compared
directly with the inclusive rate.

The result \eqref{eqn:meas_total_rate} is greater than the 
probability of the reaction $(\mu,\nu pn)$ measured by Wyttenbach et
al.~\cite{WyttenbachBaertschi.etal.1978} at 2.8\%. It is expectable because
the contribution from the $(\mu,\nu d)$ channel should be small since it
needs to form a deuteron from a proton and a neutron. 

The rate of 3.5\% was estimated with an assumption that all protons are
emitted in equilibrium. With the exponential constant of \SI{2.6}{\MeV}, the
proton yield in the range from \SIrange{40}{70}{\MeV} is negligibly small
($\sim\num{E-8}$). However, Krane and colleagues reported a significant yield
of 0.1\% in that region~\cite{KraneSharma.etal.1979}. The energetic proton
spectrum shape also has a different exponential constant of \SI{7.5}{\MeV}. One
explanation for these protons is that they are emitted by other mechanisms,
such as capture on two-nucleon cluster suggested by Singer~\cite{Singer.1961}
(see \cref{sub:theoretical_models} and
\cref{tab:lifshitzsinger_cal_proton_rate_1988}). Despite being sizeable, the
yield of high energy protons is still small (3\%) in compared with the result
in \eqref{eqn:meas_total_rate}.

%The $(\mu^-,\nu p):(\mu^-,\nu pn)$ ratio is then roughly 1:1, not 1:6 as in
%\eqref{eqn:wyttenbach_ratio}.

\subsubsection{Comparison to the silicon result}
\label{ssub:comparison_to_the_silicon_result}
The probability of proton emission per nuclear capture of 3.5\% is indeed much
smaller than that of silicon. It is even lower than the rate of the no-neutron
reaction $(\mu,\nu p)$. This can be explained as
the resulted nucleus from muon capture on silicon, $^{28}$Al, is an odd-odd
nucleus and less stable than that from aluminium, $^{27}$Mg. The proton
separation energy for $^{28}$Al is \SI{9.6}{\MeV}, which is significantly
lower than that of $^{27}$Mg at \SI{15.0}{\MeV}~\cite{AudiWapstra.etal.2003}. 

The proton spectrum from aluminium is softer than silicon charged
particles spectrum of Sobottka and Wills~\cite{SobottkaWills.1968} where the
decay constant was \SI{4.6}{\MeV}. Two possible reasons can explain this
difference in shape: 
\begin{enumerate}
  \item The higher proton separation energy of $^{27}$Mg gives less
    phase space for protons at higher energies than that in the case of
    $^{28}$Al if the excitation energies of the two compound nuclei are similar.  
  \item The silicon spectrum includes other heavier
    particles which have higher Coulomb barriers, hence contribute more in the
    higher energy bins, effectively reduces the decay rate.  
\end{enumerate}

