source: Backup NB/Talks/MEMPHYSetal/SPLFrejusArt/SPLFrejus/SPLFrejus.tex @ 416

\documentclass{iopart}
%
%----------------------------------------------------------------------------
%
\usepackage{slashbox}
\usepackage[T1]{fontenc}
\usepackage[latin1]{inputenc}
\usepackage{graphicx}
\usepackage{epsfig}
%\usepackage{vmargin}
%\usepackage{multicol}
%\usepackage[footnotesize,normal,up]{caption}
%\usepackage{subeqnarray}
%\usepackage[francais,american]{babel}
\usepackage{amssymb}
%\usepackage{amsmath}
\usepackage{setstack}
\usepackage{latexsym}


%\newcommand{\beq}[1]{\begin{equation}\label{eq:#1}}
%\newcommand{\enq}{\end{equation}}
%\newcommand{\refeq}[1]{\ref{eq:#1}}
%\newcommand{\biz}{\begin{itemize}}
%\newcommand{\eiz}{\end{itemize}}
%\newcommand{\beqa}{\begin{eqnarray}}
%\newcommand{\enqa}{\end{eqnarray}}
%\newcommand{\ol}{\overline}
%\newcommand{\virg}{\mbox{,}}
%\newcommand{\pt}{\mbox{.}}



%
%
%----------------------------------------------------------------------------
%----------------------------------------------------------------------------
%

\begin{document}
%
\title{The $\theta_{13}$ and $\delta_{CP}$ sensitivities of the SPL-Fréjus project revisited}
\author{Jean Eric Campagne, Antoine Cazes}
\address{ Laboratoire de l'Accélérateur Linéaire -
Université Paris-Sud - B\^at. 200 - BP 34 -
91898 Orsay Cedex, France}
\eads{\mailto{campagne@lal.in2p3.fr}, \mailto{cazes@lal.in2p3.fr}}
%\date{\today}
%
%
%
%\selectlanguage{francais}
%
\begin{abstract}
An optimization of the CERN SPL beam line has been performed guided by the sensitivity to the $\theta_{13}$ mixing angle and to the $\delta_{CP}$ Dirac CP violating phase performances. A UNO-like 440 ktons water \v{C}erenkov detector located at 130~km from the target in a new foreseen Fréjus laboratory has been used as a generic detector. Concerning $\theta_{13}$ sensitivity, a gain of about $24\%$ may be reached using a $4.5$~GeV proton beam compared to the up to now considered $2.2$~GeV beam energy. This may motivate new machine developments to upgrade the nominal SPL proton beam energy.
\end{abstract}
\pacs{14.60.Pq, 14.60.Lm}

\maketitle


\section{Introduction}
The very near future of the neutrino long baseline experiments is devoted to the study of the oscillation mechanism in the range of $\Delta m^2 = \Delta m^2_{atm} \approx 2.7\,10^{-3}\mathrm{eV}^2$ \cite{SKNU04,K2KNU04} using conventional $\nu_\mu$ beams. The current K2K experiment in Japan \cite{K2KNU04}, and the forthcoming MINOS in the USA \cite{MINOS} take benefit of low energy beam to measure the $\Delta m^2$ parameter using the disappearance mode $\nu_\mu\rightarrow\nu_\mu$, while OPERA/ICARUS experiments \cite{OPERA,ICARUS} using the high energy CNGS beam \cite{CNGS} will be able to detect $\nu_\tau$ appearance. If we do not consider the LSND anomaly \cite{LSND} that will be further studied soon by MiniBooNE experiment \cite{MINIBOONE}, the three flavor family scenario will be confirmed and accommodated by a $3\times 3$ Pontecorvo-Maki-Nakagawa-Sakata (PMNS) mixing matrix \cite{PMNS} with three angles ($\theta_{12}$,$\theta_{13}$,$\theta_{23}$) and one Dirac CP phase $\delta_{CP}$.  

Beyond this medium term plan, two of the next future tasks of neutrino physics are to improve the sensitivity of the last unknown mixing angle parameter, the so-called $\theta_{13}$, and to explore the CP violation mechanism in the leptonic sector. The present upper bound on $\theta_{13}$ is $\sin^22\theta_{13}<0.1$ for large $\Delta m^2$ at $90\%$~CL \cite{CHOOZ}. This sensitivity can be improved using reactor and accelerator experiments. In reactor experiments, one uses $\bar{\nu}_e$ in disappearance mode and may reach $\sin^22\theta_{13}<0.03$ \cite{Wpaper}. In accelerator experiments, one can use $\nu_e$ and $\bar{\nu}_e$ from $\beta$ beams \cite{BETABEAM} in both  disappearance and appearance modes (\textsl{i.e} $\stackrel{\scriptscriptstyle (-)}{\nu}_e\rightarrow \stackrel{\scriptscriptstyle(-)}{\nu}_\mu$), and also $\nu_\mu$ in appearance mode (\textsl{i.e} $\nu_\mu\rightarrow\nu_e$) with conventional beams either with sub-mega watt proton drivers \cite{NOVA,T2K} or with multi-mega watt proton drivers \cite{T2K,BNLHS,CERN}. The later neutrino beam type, called Superbeam, is foreseen to be extended to produce $\nu_\mu$ beam and $\bar{\nu}_\mu$ beam from muon decays, the so-called Neutrino Factory, in order to study the eventual leptonic CP violation. Such Neutrino Complex is under study in Japan, in USA and also in Europe at CERN and details may be found in reference \cite{CERN}. A comparison of the performances of $\beta$ beam and Superbeam may be found for instance in reference \cite{DONINI}.  The reactor experiment result on $\theta_{13}$ is straight forward as compared to Superbeam and Neutrino Factory results that are on one hand reacher but in an other hand more complex to analyse due to the interplay between the different physics factors $\theta_{13}$, $\delta_{CP}$, sign$(\Delta m^2_{23})$, sign$(\tan(2\theta_{23}))$ \cite{DONINI,DOUBLE-CHOOZ}. 

This paper presents results of a new simulation of the SPL (Super Proton Linac) Superbeam that could take place at CERN \cite{SPL}, using for definitiveness a UNO-like 440kT fiducial water \v{C}erenkov detector \cite{UNO} located in a new enlarged underground laboratory under study in the Fréjus tunnel, $130$~km away from CERN \cite{mosca}. The SPL neutrino beam is created by  decays of pions, muons and kaons produced by the interactions of a $4$~MW proton beam impinging a liquid mercury jet \cite{CERN}. Pions, muons and kaons are collected using two concentric electromagnetic lenses (horns), the inner one and the outer one are hereafter called "Horn" and "Reflector" respectively \cite{Meer}. Horns are followed by a decay tunnel where most of the neutrinos are produced. A sketch of the beam line is shown on figure~\ref{fig:Sbeam}.  

\begin{figure}
\centering
\includegraphics[scale=0.2]{../picts/cernmeg3_6-02.eps}
\caption{\label{fig:Sbeam}Sketch of the SPL neutrino Superbeam from CERN to the Fréjus tunnel.} 
\end{figure}

The analysis chain consists of different stages: the simulation of the interactions between the proton beam and the mercury target, the propagation of the resulting secondary particles through the magnetic field and the materials of the horns, the tracking of $\pi^\pm$, $K^{\pm,0}$ and $\mu^\pm$ until they decay, the computation of the neutrino flux at the detector site, and finally the statistical analysis. A part of the simulation chain has already been described in reference \cite{nuFact134,MMWPSCazes}. 


Compared to recent papers on the same subject \cite{DONINI,JJG,Mezzetto}, we have reoptimized the Horn and Reflector shapes \cite{nuFact138}, and introduced the kaon background simulation which allows us to update the SPL beam energy. The organization of this document follows the simulation chain: the interaction between the proton beam and the mercury target is presented in the second section. The kaon production is detailed in the third section. The simulation of the horns is described in the forth section, while the algorithms used to compute the neutrino fluxes are explained in the fifth section. Then, the sensitivities to $\theta_{13}$ and $\delta_{CP}$ are revisited with new studies about the optimization of the proton beam energy, the pion collection and the decay tunnel geometry.

\section{Target simulation}
\label{sec:target}
Since hadronic processes are crucial to describe the interactions of the proton beam on the target, the FLUKA simulator \cite{fluka} has been chosen for this first step of the simulation. 

\Table{\label{tab:targ}Liquid mercury jet parameters.}
    \br
     \centre{2}{Hg target}\\
     \mr
      Hg jet speed & $20$~m/s \\
      density & $13.546$ \\
      Length, radius & $30$~cm, $7.5$~mm \\
      \br
\endTable
The target used in the present study is a mercury liquid jet \cite{CERN} simulated by a cylinder $30$~cm long (representing two hadronic lengths) and $1.5$~cm diameter (see table~\ref{tab:targ}). The pencil like simulated proton beam is composed with  $10^6$ mono-energetic protons. The beam axis is also the symmetry axis of the target and the horns and the decay tunnel. Simulations have been performed for $2.2$~GeV proton kinetic energy, the up to now nominal design \cite{SPL}, as well as for $3.5$~GeV, $4.5$~GeV, $6.5$~GeV and $8$~GeV according to possible new designs \cite{MMWPSGaroby}.


Particle production yields are summarized in table~\ref{tab:nbPart}. The secondary proton and neutron rates induce important radiation damages and power dissipation in the horns which have been addressed in reference \cite{nuFact134}, and which will require specific R\&D effort. At $2.2$~GeV, kaon yields are very low, but it has a dramatic energy dependence as further studied in section~\ref{sec:kaon}. It is worth to mention that the numbers in table~\ref{tab:nbPart} are not to be taken as face values, because the cross sections of pion and kaon productions using proton beam are still under studies as for instance by the HARP experiment \cite{harp} and in the future by the MINERVA experiment \cite{minerva}.

\fulltable{\label{tab:nbPart}Average number of the most relevant secondary particles exiting the $30$~cm long, $1.5$~cm diameter mercury target per incident proton (FLUKA). The $\mu^+/\mu^-$ numbers and the $K^+/K^0$ numbers have been multiplied by $10^4$. Note that the $K^-$ production rate is at the level of $10^{-5}$ per incident proton.}
\br
$E_k$ (GeV) & p & n & $\gamma$ & $e^-$ & $e^+$ & $\pi^+$ & $\pi^-$ & $\mu^+$ & $\mu^-$ & $K^+$ & $K^0$ \\
\mr
2.2 & $1.4$     &  $17$   & $5.0$  &  $0.17$   &  $0.08$    &  $0.24$   &  $0.18$ & $4$ & $1$ & $7$ & $6$ \\

3.5 &$1.8$    &  $17$  & 7.0  &  0.28   &  0.15    &  0.41   &  0.37 & $10$ & $3$ & $35$ & $30$ \\

$4.5$ &  $2.3$  &  $17$ & $7.7$ &  $0.35$  & $0.21$   &  $0.57$ & $0.39$&  $11$ & $3.3$  &  $93$ &  $68$       \\

$8$   &  $3.1$  &  $17$ & $11.0$ & $0.63$  & $0.41$   &  $1.00$  &  $0.85$ &  $30$ &  $9.5$ &  $413$ &    $340$         \\

\br
\endfulltable

The cross section uncertainties are the main source of discrepancy between simulator programs. Some comparisons between FLUKA and MARS \cite{MARS} have already been presented in the same context \cite{nuFact134}. The energy distribution of the pions exiting the target, computed with the two simulator programs FLUKA and MARS, is shown on figure~\ref{fig:compFlukaMars}a. Pions come from $\Delta$ decays which in turn are originated from two different sources. At low energy ($P<200$~MeV/c), $\Delta$ are produced by protons of the target excited by the beam interactions, while the highest energy part of the spectrum is due to transformation of protons of the beam into $\Delta$. The discrepancy is quite large for the low energy part. However, the horns are designed to focus the high energy part of the spectrum (see section~\ref{sec:horn}), and therefore, MARS and FLUKA are in better agreement for the energy spectrum computed at the entrance of the decay tunnel (figure~\ref{fig:compFlukaMars}b). So,  the discrepancy  at low energy between MARS and FLUKA does not matter too much for the present application.

\begin{figure}
\centering
\includegraphics[scale=0.6]{../picts/compMarsFluka.eps}
\caption{\label{fig:compFlukaMars}Pion momentum distribution at the exit the target (a) and at the exit of the horns (b), simulated by FLUKA (\dashed) and by MARS (\full).}
\end{figure}
%
\section{Kaon production}
\label{sec:kaon}
The possibility to increase the SPL energy in order to study the optimization of the physics program has been recently pointed out \cite{MMWPSGaroby}. Then, the kaon production should be clearly addressed because it is a source of $\nu_e$ and $\bar{\nu}_e$ background events. The kaon decay channels and branching ratios are presented in table~\ref{tab:BRKP0SL} in \ref{sec:kaons}.

The target simulation described in section~\ref{sec:target} has been used with $500,000$ p.o.t with kinetic energy uniformly distributed between $2.2$~GeV and $5$~GeV. The momentum of outgoing pions and kaons are recorded when they exit the target. The number of produced $K^{o,\pm}$ at different proton beam energies are presented on figure~\ref{fig:KaonsPions}a. In one hand the $K^o$ production rate is similar to the $K^+$ production rate, but in an other hand the $K^-$ production rate is almost forty times smaller. For comparison, the numbers of $\pi^+$ and $\pi^-$  produced in the same conditions are presented on figure~\ref{fig:KaonsPions}b. Pion production rate is about two orders of magnitude greater than the kaon production rate.

\begin{figure}
\centering
\includegraphics[scale=0.6]{../picts/KaonProd.eps}
\caption{\label{fig:KaonsPions}Kaon production (a) as a function of the incident proton beam energy ($E_p$) for $500~000$ incident protons with (\full) curve for $K^+$, (\dashed) curve for $K^-$ and (\dotted) curve for $K^o$. Pion production (b) in the same conditions with (\full) curve for $\pi^+$ and (\dashed) curve for $\pi^-$.}
\end{figure}

The behavior of the two pion and kaon production rates are quite different. The pion yield grows smoothly with the proton energy while the production of kaons seams to have two origins, which has been confirmed by FLUKA's authors \cite{FLUKAprivate}. For beam energy below approximatively $4$~GeV, the resonance production model is used, and one notices a low production rate with a maximum at about $3.4$~GeV. For beam energy above  $4$~GeV, the dual parton model is used, and the production rate experiences a threshold effect with a rapid rise. The ratio between positive kaon and pion production rates is about $0.5\%$ between $2.2$~GeV and $4$~GeV and grow up to $2.5\%$ at $5$~GeV. One notices that the maching between the two kaon production models may not be optimal.
%
\section{Horn simulation}
\label{sec:horn}
The simulation code of the electromagnetic horns is written using GEANT 3.2.1 \cite{geant} for convenience and since electromagnetic processes are dominant, FLUKA has not been considered as mandatory, but this may be revised in a future work. The geometry of the horns has been inspired by an existing CERN prototype and a Reflector design proposed in reference \cite{SIMONE1}. Depending of the current injection, only positive secondary particles or negative secondary particles are focused. The relevant parameters are detailed in table~\ref{tab:specif}.

The mercury target is localized inside the Horn because of the low energy and the large emittance of the secondary pions produced: $<P_{\pi T}>/<P_\pi> \approx 240$~MeV$/400$~MeV ($2.2$~GeV proton beam energy). This explains the Horn design (figure~\ref{fig:plan}), with a cylindrical part around the target, called the neck, which is larger than the transversal size of the target to simulate the room for target handling, and a conic part designed such that the relevant pions are focused as much as possible to exit the magnetic field parallel to the beam axis.

\Table{\label{tab:specif}Relevant parameters of horns. The shapes of the conductors are not changed by proton beam energy changes, as the focusing has been optimized for a defined pion momentum.}
\br
\centre{2}{inner horn} & \centre{2}{reflector} \\
\mr
  neck inner radius & $3.7$~cm              & neck inner radius & $20.3$~cm \\
  neck length & $40$~cm                     & end cone inner radius & $35.7$~cm \\
  end cone inner radius & $15.7$~cm         & outer radius & $40$~cm \\
  outer radius & $20.3$~cm                  & total length & $200$~cm \\
  total length & $180$~cm                   & Peak current & $600$~kA \\
  Alu thickness  & $3$~mm                   & Frequency & $50$~Hz \\
  Peak current & $300$~kA                   & & \\
  Frequency & $50$~Hz                       & & \\
\br
\endTable

The shape of the horns conductors is a crucial point since it determines the energy spectrum of the neutrino at the detector site. The details of the conductor shape optimization for the present context may be found in reference \cite{nuFact138}. We just recall here some ingredients. For a $\theta_{13}$ driven $\nu_\mu \rightarrow \nu_e$ oscillation, a $\Delta m^2_{23}$ parameter value of $2.5\ 10^{-3}\mathrm{eV}^2$, and a baseline distance of $130$~km, the first oscillation maximum occurs for a neutrino energy of $260$~MeV\footnote{If one includes an energy spectral shape of the neutrino beam with a typical $\sigma_E/E \approx 1/3$, then a 10\% increase of the average neutrino energy is expected. But, the induced change on the $\theta_{13}$ sensitivity is negligeable.}. The   optimization of the physics potential depends at first approximation on the pion neutrino characteristics, which energy is fully determined by the pion 2-body decay and boost. To reach an energy of $260$~MeV, the pion needs a $\beta=0.97$, which in turn induces a pion momentum of $600$~MeV/c. Then, the shape of the conic part of the horns is determined such that these $600$~MeV/c pions exit parallel to the beam axis.

\begin{figure}
\centering
\includegraphics[width=0.75\textwidth]{../picts/corneSchema.eps}
\caption{\label{fig:plan}Design of the Horn and the Reflector conductor shapes implemented in the GEANT simulation. The Hg target is located inside the cylindrical part of the Horn.}
\end{figure}

Before closing this section, it is worth to quote that the Horn/Reflector conductor shapes optimized in the present study to focus a given pion momentum value, is not affected by a proton beam energy change. What is affected is the production rate of the relevant pions. This Horn/Reflector design consideration would be different if one wishes to focus as much as possible all the pions produced for which the mean energy is of course affected by a proton beam energy change.    
%
\section{Particle decay treatment and flux calculation}
%
The decay tunnel representation is a simple cylinder with variable length ($L_T$) and radius ($R_T$) filled with "vacuum". The base line design is a $20$~m long and $1$~m radius cylinder, but simulations have also been conducted with lengths of $10$~m and $40$~m, and radius of $1.5$~m in the spirit of reference \cite{donega}. In the GEANT simulation, to gain in CPU time, only pions, muons and kaons are tracked in the volume of the tunnel, and all particles exiting this volume are discarded. 


Beyond the $1/L^2$ solid angle factor due to the source-detector distance ($L$) which decreases dramatically the fluxes, the neutrino beam focusing is very limited due to the small pion boost factor ($\approx 4$). Therefore, computational algorithms have been used to avoid a too prohibitive CPU time resulting from the simulation of each secondary particle decay. Otherwise, about $10^{15}$ p.o.t would have been necessary to obtain reliable statistics for the estimation of the $\bar{\nu}_e$ flux for instance.

It is worth to mention that the particle decays occurring before the entrance of the decay tunnel are also taken into account and treated in the same manner, which is not the case in reference \cite{donega}. 

\subsection{Algorithm description}
\label{sec:algo}

The decay code has been included in the GEANT code. The basic idea of this algorithm is to compute the neutrino fluxes using the probability of reaching the detector for each neutrino produced by a $\pi$ or a $K$ or a $\mu$ particle (on-axis neutrino beam). This method has already been used in reference \cite{donega} and has been modified and extended to the kaon decay chain for the present study. 

Muon neutrino comes mostly from pion decay. In a first stage, each pion is tracked by GEANT until it decays. Then, the probability for the produced muon neutrino to reach the detector is computed. The flux is obtained applying the probability as a weight for each neutrino. All the pions produced in the simulation are therefore useful to compute the flux, and this allows to reduce the number of events in the simulation to $10^6$ p.o.t. In this computation, the decay region (horns and tunnel) is considered as point like compared to the source-detector distance. 

The same method is applied for neutrino coming from muons and kaons with some modifications  because most of the muons do not decay, and there are very few kaons produced (see table~\ref{tab:nbPart}).
The probability computation is presented in \ref{sec:decayprobcomp}. 

\Table{\label{tab:proton}Number of protons on target for different beam energy at 4~MW constant power.}
\br
                        Beam energy  & number of proton \\
                          (GeV)      & per year ($10^{23}$ p.o.t/y) \\
\mr
                        2.2 & 1.10 \\
                        3.5 & 0.70 \\
                        4.5 & 0.56 \\
                        6.5 & 0.40 \\
                        8.0 & 0.30 \\
\br
\endTable
\subsection{Validation of the algorithm}
The validity of the method presented in the previous section have been tested against a straight forward  algorithm consisting of decaying each pion $N$ times ($N \approx 10^6$) in a full GEANT simulation of the event (decays included). Such method presents the advantage to keep all the information of the neutrino available for further studies. It can be a good approach to compute the muon neutrino flux coming from pion decays. It can also provide the beam profile, but it shows its limits for the muon induced fluxes, especially the $\bar{\nu}_e$ flux. Indeed, this means that each muon is duplicated $N$ times and when a muon decays, it must decay $N$ times again. For $N\approx 10^6$, this is a prohibitive CPU time consuming.
\begin{figure}
\centering
\includegraphics[scale=0.6]{../picts/compGeantDonega.eps}
\caption{\label{fig:compGeantDonega}Comparison between the probability method, (\full) curve, and the full GEANT simulation method, (\dashed) curve, for the $\nu_\mu$ from $\pi^+$ flux (left) and the $\bar{\nu}_\mu$ from $\pi^-$ flux  (right). The horns are set to focus positive particles.}
\end{figure}

The $\nu_\mu$ and $\bar{\nu}_\mu$ fluxes are displayed on figure~\ref{fig:compGeantDonega} for both methods. The two spectra shows a clear agreement, and this makes confidence on the probability method.
%
\subsection{Simulated fluxes}
The fluxes are computed at a distance of $100$~km from the source by convention and can be rescaled at any desired distance. They provide the number of the four neutrino species ($\nu_\mu$, $\bar{\nu}_\mu$, $\nu_e$, $\bar{\nu}_e$) passing through a $100$~m$^2$ fiducial area during $1$~year.

In practice, the fluxes are given as a function of the neutrino energy via histograms composed of $20$~MeV bin width.  These histograms are filled with the energy of each neutrino weighted by the probability to reach the detector (section~\ref{sec:algo}). To obtain the fluxes, the histograms are rescaled to the number of p.o.t per year depending on the beam energy. Table~\ref{tab:proton} reports on the number of p.o.t per year for the different energies studied using the definition of one year being $10^7$~s and keeping the beam power constant (\textsl{i.e} $4$~MW).

Three origins are identified in the composition of each neutrino flux:
\begin{itemize}
        \item neutrinos from pions, which includes neutrinos created by primary pion decays and neutrinos coming from the muons produced by pion decays or muons directly exiting the target. This is the component studied in reference \cite{donega} but with different settings and event generator;
        \item neutrinos emitted during the decay chain of the charged kaons, either by direct production, or produced by the daughter pions and muons;
        \item neutrinos coming from the decay chain of the neutral kaons.
\end{itemize}

The three components of the fluxes for the four neutrino species are presented on figure~\ref{fig:flux22p} for positive particle focusing and a proton beam kinetic energy of $2.2$~GeV. The $\nu_\mu$ flux is dominated by the neutrinos of pion decays, but a tail above $500$~MeV (zoom on the top left part) is created by the $K^+\rightarrow \mu^+\nu_\mu$ channel, which is anyway at least three order of magnitude below the flux maximum. The $\bar{\nu}_\mu$ flux is mostly due to the decays of $\pi^-$ that are not unfocused by the horns, but the higher energy part comes from $\mu^+$ decays. It is noticable that the $\nu_e$ and $\bar{\nu}_e$ fluxes are respectively more than $200$ and more than $7000$ times smaller than the $\nu_\mu$ flux. The $\bar{\nu}_e$ are produced in a large part by the $K^0_L\rightarrow\pi^+e^-\bar{\nu}_e$ decay channel and by $\mu^-$ decays, while the  $\nu_e$ flux is dominated by the $\mu^+$ decays.

On figure~\ref{fig:flux22m}, the horns are set to focus negative particles keeping other parameters identical. By comparison with positive focusing, one can at first approximation translate the results by exchanging particles and anti-particles, except that the $K^+/K^-$ ratio is about 100  in the beam-target interactions (see table~\ref{tab:nbPart}).

On figures~\ref{fig:flux45p} and \ref{fig:flux8p}, one observes the evolution of figure~\ref{fig:flux22p} when the proton beam kinetic energy increases to $4.5$~GeV and $8$~GeV, respectively. Correspondingly, the results for negative particle focusing are presented on figures~\ref{fig:flux45m} and \ref{fig:flux8m}. One clearly notices the increase of the kaon induced neutrino contents as the beam energy grows.

\begin{figure}
\centering
\includegraphics[scale=0.6]{../picts/flux22p.eps}
\caption{       \label{fig:flux22p}Neutrino fluxes $100$~km from the decay region and with the horns focusing the positive particles. The fluxes are computed for a SPL proton beam of $2.2$~GeV (4~MW), a decay tunnel with a length of $20$~m and a radius of $1$~m. The (\full) curve is the contribution from primary pions and the daughter muons, the (\dashed) curve is the contribution from the charged kaon decay chain, and the (\dotted) curve is the contribution from the $K^0$ decay chain.}
\end{figure}

\begin{figure}
\centering
\includegraphics[scale=0.6]{../picts/flux22m.eps}
\caption{\label{fig:flux22m}Same legend as for figure~\ref{fig:flux22p} but the horns are focusing negative particles.}
\end{figure}

\begin{figure}
\centering
\includegraphics[scale=0.6]{../picts/flux45p.eps}
\caption{\label{fig:flux45p}Same legend as for figure~\ref{fig:flux22p} but for proton beam kinetic energy of $4.5$~GeV (4~MW).}
\end{figure}

\begin{figure}
\centering
\includegraphics[scale=0.6]{../picts/flux45m.eps}
\caption{\label{fig:flux45m}Same legend as for figure~\ref{fig:flux22m} but for proton beam kinetic energy of $4.5$~GeV (4~MW).}
\end{figure}

\begin{figure}
\centering
\includegraphics[scale=0.6]{../picts/flux8p.eps}
\caption{\label{fig:flux8p}Same legend as for figure~\ref{fig:flux22p} but for proton beam kinetic energy of $8$~GeV (4~MW).}
\end{figure}

\begin{figure}
\centering
\includegraphics[scale=0.6]{../picts/flux8m.eps}
\caption{\label{fig:flux8m}Same legend as for figure~\ref{fig:flux22m} but for proton beam kinetic energy of $8$~GeV (4~MW).}
\end{figure}

\fulltable{\label{tab:speciesfluxes}Integral of the different species fluxes with different settings. The $\nu_\mu$ and $\bar{\nu}_\mu$ fluxes are expressed in $10^{13}/100\mathrm{m}^2/y$ unit while the $\nu_e$ and $\bar{\nu}_e$ fluxes are expressed in $10^{11}/100\mathrm{m}^2/y$ unit. The positive focusing  and negative focusing  are distinguished by a ($+$) sign and a ($-$) sign, respectively. The settings used corresponds to different values of $L_T$ and $R_T$, the length and radius of the decay tunnel. Setting (1) is the baseline option and means $L_T = 20$~m and $R_T = 1$~m, while setting (2) means $L_T = 10$~m and $R_T = 1$~m and setting (3) means $L_T = 40$~m and $R_T = 1$~m, and finally the setting (4) means $L_T = 20$~m and $R_T = 1.5$~m.}

\br

                        Settings  & \centre{2}{$\nu_\mu$} & \centre{2}{$\nu_e$} 
                                                                &       \centre{2}{$\bar{\nu}_\mu$} & \centre{2}{$\bar{\nu}_e$} \\
                                                                \cline{2-9}
                                                                & $+$ & $-$ & $+$ & $-$ & $+$ & $-$ 
                                                                & $+$ & $-$ \\ 
\mr

(1): $2.2$~GeV  
        & $7.6$         & $0.4$ 
        & $3.2$                                 & $0.2$
        & $0.3$                                 & $5.8$  
        & $0.1$                                 & $1.6$ \\                

(1): $3.5$~GeV  
        & $10.0$          & $0.9$ 
        & $4.4$           & $0.6$
        & $0.7$                                         & $8.5$  
        & $0.3$                                         & $2.2$ \\                

(1): $4.5$~GeV  
        & $10.9$          & $1.1$ 
        & $5.1$           & $1.0$
        & $0.7$                                         & $6.7$  
        & $0.4$                                         & $1.8$ \\                

(1): $6.5$~GeV  
        & $10.4$          & $1.4$ 
        & $6.4$           & $2.0$
        & $1.0$           & $7.1$  
        & $0.9$                                         & $2.5$ \\      

(1): $8.0$~GeV  
        & $9.7$          & $1.5$ 
        & $6.7$           & $2.3$
        & $1.2$           & $7.1$  
        & $1.1$                                         & $2.8$ \\


(2): $4.5$~GeV  
        & $7.1$          & $1.0$ 
        & $2.8$          & $0.9$
        & $0.5$          & $5.2$  
        & $0.3$                                   & $1.1$ \\              

(3): $4.5$~GeV  
        & $12.3$          & $1.2$ 
        & $7.2$           & $1.0$
        & $0.8$           & $7.5$  
        & $0.4$                                         & $2.6$ \\      

(4): $4.5$~GeV  
        & $13.2$          & $1.5$ 
        & $6.9$           & $1.4$
        & $0.9$           & $8.1$  
        & $0.6$                                         & $2.7$ \\
        
\br
\endfulltable

On table~\ref{tab:speciesfluxes} are reported the integral of the fluxes when one modifies the decay tunnel length and radius, as well as the beam kinetic energy.  Switching the radius from $1$~m to $1.5$~m has increased the $\nu_\mu$ and $\bar{\nu}_\mu$ fluxes by $20\%$ but also in paralell the $\nu_e$ and $\bar{\nu}_e$ fluxes by $38\%$. Increasing the decay tunnel length from $10$~m up to $40$~m yields an increase of the signal like flux by $75\%$ and also in paralell an increase of the background like flux by a factor $1.7$. Notice that the $\nu_\mu/\bar{\nu}_\mu$ flux ratio is rather insensitive to the decay tunnel length. The feeling that $L_T = 20$~m and $R_T = 1$~m is a good signal over background compromise is confirmed by sensitivity quantitative studies reported in section~\ref{sec:results}.  

Looking at the evolution of $\nu_\mu$ flux with respect to the beam energy, one notices that a maximum is reached around $4.5$~GeV. This is due to the competition between the cross section rise with respect to the energy and the decrease of the number of p.o.t due to the constant SPL power ($4$~MW).
%
\section{Sensitivity computation ingredients}
%
The sensitivity to $\theta_{13}$ and $\delta_{CP}$ is computed for an $\nu_\mu\rightarrow\nu_e$ appearance experiment. An analysis program described in reference \cite{MEZZETTONUFACT060} has been used for such sensitivity computation. See table~\ref{tab:param} for the default user parameter values used in this paper. We just remind here some key points of the program. 

It is included a full 3-flavors oscillation probability computation with matter effects, but no ambiguities are taken into account. This latest point may be revisited in a future work using reference \cite{DONINI-2}. Concerning the background events, the $\nu_e/\bar{\nu}_e$ from the beam, the $\nu_\mu e^-$ elastic scattering process, the $\pi^o$ production as well as the $\mu/e$ misidentification are taken into account. The cross-sections from the NUANCE program are used \cite{NUANCE}. The systematics error is a user parameter and we have used the $2\%$ value considered as a final goal, but also $5\%$ and $10\%$ \cite{MEZZETTONUFACT060}. The detector considered for definitiveness is similar to the UNO detector, i.e. a $440$~kt fiducial water \v{C}erenkov detector \cite{UNO}. It is located at $L = 130$~km from CERN, in the foreseen new Fréjus laboratory \cite{mosca}. It is worth to mention that if one wants to evaluate the influence of $L$ on the sensitivity, it would mean a re-optimization of the horns for each $L$ envisaged (see section~\ref{sec:horn}).  

\Table{\label{tab:param}Default user parameters used to compute the sensitivity curves \cite{MEZZETTONUFACT060}. The quoted errors in parenthesis for the $(12)$ and the $(23)$ parameters (absolute value for the masses and relative value for the angles) are coming respectively from the up to date combined Solar and KamLAND results \cite{KAMLAND} and from a 200 ktons-years SPL desappearance exposure \cite{JJG}.}
\br
                        $\Delta m^2_{12} = 8.2 (0.5)\times 10^{-5}~\mathrm{eV}^2$ & $\sin^22\theta_{12} = 0.82 (9\%)$ \\
                        $\Delta m^2_{23} = 2.5 (0.1)\times 10^{-3}~\mathrm{eV}^2$ & $\sin^22\theta_{23} = 1.0 (1\%)$ \\
\mr
                  $L_T = 20$~m & $R_T = 1$~m \\ 
                  $M=440$~kT   & $\epsilon_{syst}=2\%$ \\
\br
\endTable

The running time scenario has been fixed to 5 years focusing positive particles for the $\theta_{13}$ sensitivity studies. This scenario gives comparable $\theta_{13}$ sensitivity than a mixed focusing scenario defined as 2 years focusing positive particles and 8 years focusing negative particles which is more appropriate for $\delta_{CP}$ sensitivity studies.
%
\section{Results}
\label{sec:results}
The $\theta_{13}$ and $\delta_{CP}$ sensitivities are computed with $\theta_{13} = 0^\circ$ and $\delta_{CP} = 0^\circ$ if not explicitly mentioned. 

Table~\ref{tab:nbvsE} presents the number of signal and background events for a $5$ years positive focusing experiment, but with different beam energy settings. The significance parameter is defined in reference \cite{MEZZETTONUFACT060} as:
\begin{equation}
\mathcal{S} = \frac{N^{osc}_{\nu_e}}{\sqrt{N^{osc}_{\nu_e} + N^{beam}_{\nu_e} + N^{oth. bkg} + \left(\left(N^{osc}_{\nu_e}+N^{beam}_{\nu_e}\right)\times\epsilon_{syst}\right)^2}}
\label{eq:significance}
\end{equation}
with $N^{osc}_{\nu_e}$ the number of $\nu_e/\bar{\nu}_e$ events due to $\nu_\mu/\bar{\nu}_\mu$ oscillations, $N^{beam}_{\nu_e}$ the number of background events coming from the $\nu_e/\bar{\nu}_e$ contamination of the beam, $N^{oth. bkg}$ the other kinds of background events and $\epsilon_{syst}$ the systematical factor. One can appreciate that $4.5$~GeV beam energy presents better results and may become the new base line energy.
%

\Table{\label{tab:nbvsE}Number of events for 5 years positive focusing scenario with default parameters of table~\ref{tab:param}. : $\pi^0$, $\nu_\mu$-elast., $\mu/e$-missId. The significance parameter is defined by equation~\ref{eq:significance}.}
\br

                         & $2.2$~GeV & $3.5$~GeV & $4.5$~GeV & $6.5$~GeV & $8$~GeV \\
\mr
                        non oscillated $\nu_\mu$       & $36917$  & $60969$  & $73202$  & $78024$ & $76068$ \\ 
                        oscillated $\nu_e$             & $31$  & $44$  & $47$  & $44$ & $41$ \\
                        beam $\nu_e$                   & $165$ & $222$ & $241$ & $287$ & $299$ \\
                        other background & $70$  & $104$ & $126$ & $147$ & $152$ \\
\mr
                        Significance                  &  $1.87$ & $2.18$  & $2.20$  & $1.92$ & $1.74$ \\
\br
\endTable

The contours at $90\%$, $95\%$ and $99\%$~CL of the $\theta_{13}$ sensitivity are presented in the $(\sin^22\theta_{13}, \Delta m^2_{23})$ plane on figure~\ref{fig:sensi45} for $4.5$~GeV proton beam kinetic energy. The comparison between the contours at $90\%$~CL with $2.2$~GeV, $3.5$~GeV, $4.5$~GeV and $8$~GeV beam energies is shown on figure~\ref{fig:compSensi}. One notices a better perfomence reached with a $4.5$~GeV energy beam as a confirmation of significance parameter value. But, in fact there is not much visual difference between a sensitivity obtained with $3.5$~GeV and $4.5$~GeV, even if one should keep in mind that kaon production models are different at these two energies (see section~\ref{sec:kaon}). Quantitative studies of the minimum $\sin^22\theta_{13}$ with respect to the kinetic beam energy $E_k(proton)$, and the decay length $L_T$, and the systematics $\epsilon_{syst}$ are presented in tables~\ref{tab:thvsE} and \ref{tab:thvseps}. The influance of the systematical level is presented on figure~\ref{fig:compEpsSyst}.
%
\begin{figure}
\centering
\includegraphics[scale=0.4]{../picts/sensi45.eps}
\caption{\label{fig:sensi45}Sensitivity contours obtained with a SPL energy of $4.5$~GeV and default parameters of table~\ref{tab:param}.}
\end{figure}
%
\begin{figure}
\centering
\includegraphics[scale=0.4]{../picts/compSensi.eps}
\caption{\label{fig:compSensi}Comparison of 90\% sensitivity contours obtained with SPL energies of ($2.2$, $3.5$, $4.5$, $8$)~GeV and default parameters of table~\ref{tab:param}.}
\end{figure}

\Table{\label{tab:thvsE}Minimum $\sin^22\theta_{13}\times 10^3$ observable at $90\%$ CL computed for diferent decay tunnel length ($L_T$) and kinetic beam energy ($E_k(proton)$). Other parameters are fixed to default values (table~\ref{tab:param}).}
\br
                        & $2.2$~GeV & $3.5$~GeV & $4.5$~GeV & $6.5$~GeV & $8$~GeV  \\
\mr
                        $10$~m & $1.10$ & $0.96$ & $1.05$ & $1.10$ & $1.20$ \\
                        $20$~m & $1.20$ & $0.96$ & $0.91$ & $1.05$ & $1.15$ \\
                        $40$~m & $1.26$ & $1.00$ & $1.00$ & $1.10$ & $1.20$ \\
\br
\endTable

\Table{\label{tab:thvseps}Minimum $\sin^22\theta_{13}\times 10^3$ observable at $90\%$ CL computed for different level of systematics ($\epsilon_{syst}$) and kinetic beam energy ($E_k(proton)$). Other parameters are fixed to default values (table~\ref{tab:param}).}
\br                     
                              & $2.2$~GeV & $3.5$~GeV & $4.5$~GeV & $6.5$~GeV & $8$~GeV \\
\mr                     
                        $2\%$ & $1.20$ & $0.96$ & $0.91$ & $1.10$ & $1.20$ \\
                        $5\%$ & $1.51$ & $1.26$ & $1.26$ & $1.51$ & $1.66$ \\
                        $10\%$ & $2.40$ & $2.19$ & $2.29$ & $2.75$ & $3.16$ \\
\br                     
\endTable

\begin{figure}
\centering
\includegraphics[scale=0.4]{../picts/compEpsSyst.eps}
\caption{\label{fig:compEpsSyst}$90\%$ CL sensitivity contours obtained with a SPL energy of $4.5$~GeV and default parameters of table~\ref{tab:param} but for different $\epsilon_{syst}$ values.}
\end{figure}

There are also variations on the minimum $\sin^22\theta_{13}$ value that may be reached in a $\nu_\mu \rightarrow \nu_e$ experiment which are due to the $\mathrm{sign}(\Delta m^2_{23})$ ambiguity and the $\delta_{CP}$ value. On table~\ref{tab:sign} are presented these kind of variations. Other ambiguities coming from the sign$(\tan(2\theta_{23}))$ ignorance also take place as studied in reference \cite{DONINI}. From figure~9 of this reference, we estimate a 30\% effect on $\sin^2(2\theta_{13})$ sensivity due to these ambiguities.

%\begin{figure}[htb]
%       \centering
%               \includegraphics{../picts/compDelta_thVSdm.eps}
%       \caption{$\theta_{13}$-sensitivity dependence up on the $\delta_CP$ phase obtained for a $4.5$~GeV SPL beam and 5 years running with a positive focusing scenario.}
%       \label{fig:compDelta_thVSdm}
%\end{figure}
 
\Table{\label{tab:sign}Minimum $\sin^22\theta_{13}\times 10^3$ observable at $90\%$ CL computed for different values of sign$(\Delta m^2_{23})$ and $\delta_{CP}$. Other parameters are fixed to default values (table~\ref{tab:param}).}
\br
                            & $-180^\circ$ & $-90^\circ$ & $0^\circ$ & $90^\circ$ & $180^\circ$\\
\mr
                        $+$ & $1.45$ & $0.43$ & $1.20$ & $11.48$ & $1.45$\\
                        $-$ & $1.45$ & $12.02$& $1.15$ & $0.44$  & $1.45$ \\
\br
\endTable

The combined $\sin^22\theta_{13}$ and $\delta_{CP}$ sensitivity is presented on figure~\ref{fig:compDeltaTheta}a for the 5 years positive focusing scenario to appreciate in an other way the above mentioned $\delta_{CP}$ ambiguity. 
%
\begin{figure}
\centering
\includegraphics[scale=0.4]{../picts/compDeltaTheta.eps}
\caption{\label{fig:compDeltaTheta}Sensitivity contours obtained with SPL beam energy of $2.2$~GeV (\dashed), $3.5$~GeV (\chain), $4.5$~GeV (\full) and $8$~GeV (\dotted) at $90\%$ CL. Default parameters of table~\ref{tab:param} are used either with a 5 years positive focusing scenario (a) or a mixed scenario of 2 years positive focusing and 8 years of negative focusing (b).}
\end{figure}
%
To improve the $\delta_{CP}$-independent limit on $\sin^22\theta_{13}$, one may envisaged a combination of 2 years with positive focusing and 8 years negative focusing as in references \cite{DONINI,JJG,Mezzetto}. The corresponding combined sensitivity contours are presented in figure~\ref{fig:compDeltaTheta}b. With this kind of mixed focusing scenario, one gets a limit $\sin^22\theta_{13} < 2.4 \times 10^{-3}$ ($90\%$ CL) independently of $\delta_{CP}$ with a $3.5$~GeV kinetic energy beam. Note that with a $2.2$~GeV beam or a $4.5$~GeV beam, one gets $\sin^22\theta_{13} < 2.7 \times 10^{-3}$ which is just a bit worse, but it is better than the limit of $\sin^22\theta_{13} < 3.2 \times 10^{-3}$ that may be reached with a $8$~GeV beam.
%
\section{Summary and outlook}
A complete chain of simulation has been set up for the SPL-Fréjus project. The neutrino production has been extended to the kaon decay contribution, which is important to test SPL energy scenario above $2.2$~GeV.

The beam line optimization has been performed using the sensitivity to $\sin^22\theta_{13}$. The shape of the focusing system has been updated and we have confirmed that the decay tunnel parameters $20$~m long and $1$~m radius may still be conserved as baseline values. 

With a $4.5$~GeV SPL energy beam, in 5 years of focusing positive particles, one can reach a minimum value of $\sin^22\theta_{13} = 9.1\times 10^{-4}$ ($90\%$ CL, $\delta_{CP} = 0$) and realizes a $24\%$ gain compared to the up to now nominal $2.2$~GeV SPL energy beam. This may be the sign that a $4.5$~GeV SPL beam energy design is to be investigated from the machine development implication point of view. If one includes  a mixed focusing scenario to explore the ($\sin^22\theta_{13}$,$\delta_{CP}$) sensitivity, then their is no much difference between the $\delta_{CP}$-independant $\sin^22\theta_{13}$ minimum value that may be reached with a $2.2$~GeV energy beam and a $4.5$~GeV energy beam.

The authors think that the present study may be extended in two directions. The beam line simulation part may be performed with a single simulator as FLUKA (or GEANT4 \cite{GEANT4} for comparison). The sensitivity analysis may be deeper investigated using the complete set of possible ambiguities as in reference \cite{DONINI}, and the $\theta_{13}$ or $\delta_{CP}$ measurement accuracy with new beam energy scenario may be investigated too.    
%
\ack
The authors would like to thank M.~Mezzetto to have expressed his interest at early stage of this work and to have provided us his sensitivity computation program. Also the authors thank S.~Gilardoni for fruitful discussions.
%
\appendix
\setcounter{section}{1}
\section{Decay probability computations}
\label{sec:decayprobcomp}
This appendix contains the probability formulas and the algorithms used in the flux computation (see section~\ref{sec:algo}).

\subsection{Pion neutrino probability computation}
\label{sec:Ppi}
Pions decay only as $\pi^+\rightarrow \mu^+ + \nu_\mu$ or $\pi^- \rightarrow \mu^- + \bar{\nu}_\mu$ and the neutrinos are emitted isotropically in the pion rest frame, with an energy of about $30$~MeV given by the 2-body decay kinematics. Applying a Lorentz boost knowing the pion momentum and direction, it is possible to compute the probability to reach the detector for the neutrinos. Only neutrino parallel to the beam axis are supposed to pass through the detector fiducial area, and therefore, the neutrino must be emitted by the pion with an angle opposite to the angle between the pion and the beam axis (see figure~\ref{fig:pionDecay}). This gives:
\begin{equation}
\mathcal{P}_\pi = \frac{1}{4\pi}\frac{A}{L^2}\frac{1-\beta^2}{(\beta\cos\alpha-1)^2}
\label{probaPi}
\end{equation}
where $\beta$ is the velocity of the pion in the tunnel frame, $A$ is the fiducial detector surface, $L$ the distance between the neutrino source and the detector, and $\alpha$ the angle between the pion direction and the beam axis in the laboratory frame.

\begin{figure}
\centering
\includegraphics{../picts/pionDecay.eps}
\caption{\label{fig:pionDecay}Pion decay in the tunnel frame. To reach the detector, $\delta = -\alpha$ is needed.}
\end{figure}

\subsection{Muon neutrino probability computation}
\label{sec:Pmu}
Muons decay only as $\mu^+ \rightarrow e^+ + \nu_e + \bar{\nu}_\mu$ or $\mu^- \rightarrow e^- + \bar{\nu}_e + \nu_\mu$, and will produce background events. The mean decay length of the muons is $2$~km, therefore, most of them do not decay in the tunnel. This induces a lake of statistics to estimate the corresponding level of background. This problem has been solved using each muon appearing in the simulation in the following steps:
\begin{enumerate}
        \item the probability for the muon to decay into the tunnel has been computed using a straight line propagation;
        \item the available energy for the neutrino in the tunnel frame has been divided in $20$~MeV energy bins; 
        % (the number of bin depending of the muon energy).
        \item one $\nu_e$ and one $\nu_\mu$ have been simulated in each of the energy bins (step 2). Then, the probability to reach the detector has been computed, and multiplied by the probability computed at step 1.%the neutrino energy is stored in an histogram, weighted with this probability and 
\end{enumerate}
After the probability computation, the non useful muon is discarded by GEANT to gain in CPU time. 

The probability for the muon neutrino and the electron neutrino to be emitted parallel to the beam axis is \cite{donega}:
\begin{equation}
\fl \frac{d\mathcal{P}_\mu}{dE_\nu} = \frac{1}{4\pi}\frac{A}{L^2}\frac{2}{m_\mu}\frac{1}{\gamma_\mu(1+\beta_\mu\cos\theta^*)}\frac{1-\beta_\mu^2}{(\beta_\mu\cos\rho-1)^2}\left[f_0(x)\mp \Pi_\mu^T f_1(x)\cos\theta^*\right]
\label{probaMu}
\end{equation}
where $\beta_\mu$ and $\gamma_\mu$ are the velocity and the Lorentz boost of the muon in the tunnel frame, $\theta^*$ is the angle with respect to the beam axis of the muon in the muon rest frame, $\rho$ is the corresponding angle in the tunnel frame. Like in the pion case, this angle appears because the neutrino must be parallel to the beam axis. $\Pi_\mu^T$ is the muon transverse polarization, the parameter $x$ is defined as $x=2E_\nu^*/m_\mu$ where $E_\nu^*$ is the neutrino energy in the muon rest frame, and the function $f_0(x)$ and $f_1(x)$ coming from the matrix element of the muon decay are given in table~\ref{tab:Function}. The sign in front of $\Pi_\mu^T$ in equation~\ref{probaMu} is $(-)$ for the $\nu_\mu$ and $(+)$ for the $\bar{\nu}_\mu$, respectively.

\Table{\label{tab:Function}Flux function in the muon rest frame \cite{Gaisser}.}
\br
                      &  $f_0(x)$ & $f_1(x)$ \\
\mr
                $\nu_\mu$ & $2x^2(3-2x)$ & $2x^2(1-2x)$ \\
                $\nu_e$   & $12x^2(1-x)$ & $12x^2(1-x)$ \\
\br
\endTable

Muon polarization is computed using the conservation of the transverse component of the velocity four-vector $\gamma(1,\beta)$ between the muon rest frame (where the polarization is computed) and the pion rest frame, where the muon helicity is $-1$, due to the parity non conservation. It yields \cite{picasso}:
\begin{equation}
\Pi_\mu^T = \frac{\gamma_\pi\beta_\pi}{\gamma_\mu\beta_\mu}\sin\theta^*
%\mbox{ and }\Pi_\mu^L = \sqrt{1-\Pi_\mu^{T2}}\virg 
\label{pola}
\end{equation}
where $\gamma_\pi$, $\beta_\pi$, $\gamma_\mu$, and $\beta_\mu$ are the Lorentz boost and velocity of the pion and of the muon in the tunnel frame, and $\theta^*$ the angle with respect to the beam axis of the muon in the pion rest frame.

\subsection{The treatment of the kaons}
\label{sec:kaons}
Contrary to pions and muons, kaons have many decay channels. They are summarized in table~\ref{tab:BRKP0SL}.

\Table{\label{tab:BRKP0SL}Charged and neutral kaon decay channels \cite{pdg}.}
\br
\centre{2}{$K^\pm$}         & \centre{2}{$K^0_L$}     & \centre{2}{$K^0_S$}\\
\mr
$\mu^\pm\nu_\mu$ & $63.51\%$ & $\pi^-e^+\nu_e$ & $19.35\%$        & $\pi^+\pi^-$ & $68.61\%$ \\
$\pi^\pm\pi^0$ & $21.17\%$   & $\pi^+e^-\bar{\nu}_e$ & $19.35\%$  & $\pi^0\pi^0$ &  $31.39\%$ \\                
$\pi^\pm\pi^+\pi^-$ & $5.59\%$ & $\pi^-\mu^+\nu_\mu$ & $13.5\%$ & & \\
$e^\pm\nu_e\pi^0$ & $4.82\%$   & $\pi^+\mu^-\bar{\nu}_\mu$ & $13.5\%$ & & \\                                            
$\mu^\pm\nu_\mu\pi^0$ & $3.18\%$ & $\pi^0\pi^0\pi^0$ & $21.5\%$ & &\\
$\pi^\pm\pi^0\pi^0$ & $1.73\%$ & $\pi^+\pi^-\pi^0$ & $12.38\%$ & & \\                   
\br
\endTable

There is a very small amount of kaons produced (section~\ref{sec:algo}), and this number has been artificially increased in order to obtain statistically satisfactory results.  The multiplicity of decay channels makes impossible the method used for the muon case (\ref{sec:Pmu}). The method chosen for the good compromise between the gain in CPU and the statistical uncertaincy of the results, is to duplicate 300 times each kaon exiting the target. 

All the kaons daughter particles are tracked by GEANT until they decay. Three different types of daughter particles are identified in the kaon decays. The first type corresponds to primary neutrinos, the second type concerns charged pions and muons, and the neutral pions are left for the last type.

In the $K^\pm\rightarrow\mu^\pm \nu_\mu(\bar{\nu}_\mu)$ decay modes, the computation of the probability for a neutrino to reach the detector is the same than the 2-body decay formula used to in the pion decay (equation~\ref{probaPi}), where $\beta$ is now the kaon velocity, and $\alpha$ the angle of the kaon with respect to the beam axis. 

When a neutrino is produced by a kaon 3-body decay, the probability to reach the detector is computed using a pure phase space formula. It yields:
\begin{equation}
\mathcal{P}_K = \frac{1}{4\pi}\frac{A}{L^2}\frac{1}{m_K-m_\pi-m_l}\frac{1}{\gamma_K(1+\beta_K\cos\theta^*)}\frac{1-\beta_K^2}{(\beta_K\cos\delta-1)^2}
\label{probaL}
\end{equation}
where $m_K$ is the kaon mass (charged or neutral), $m_\pi$ is the pion mass ($\pi^0$ mass in $K^\pm$ decays and $\pi^\pm$ mass in $K_L^0$ decays), and $m_l$ is the mass of the lepton associated with the neutrino. The $\beta_K$ and $\gamma_K$ are the velocity and the Lorentz boost of the kaon, $\theta^*$ is the angle between the neutrino direction and the kaon direction, in the kaon rest frame. Finally, $\delta$ is the angle between the kaon direction and the beam axis in the tunnel frame.

When a $\pi^\pm$ is produced in the kaon decay chain, it is tracked by GEANT until it decays, and the probability of equation~\ref{probaPi} is applied to the produced neutrino. In case of a muon, it is treated as explained in \ref{sec:Pmu}. The muon polarization is computed this time using the kaon decay informations. Finally, when a $\pi^0$ is produced, as it cannot create neutrinos, it is simply discarded.

%
\section*{References}
\begin{thebibliography}{99}
%
\bibitem{SKNU04}
Kearns E [SuperKamiokande Collaboration] Talk at Neutrino 2004 to be published in \NP B (Proceedings Supplement)

\bibitem{K2KNU04}
Nakaya T [K2K Collaboration] Talk at Neutrino 2004 to be published in \NP B (Proceedings Supplement)

\bibitem{MINOS}
Thomson T [MINOS Collaboration] Talk at Neutrino 2004 to be published in \NP B (Proceedings Supplement) \\
Adamson P \etal NuMi-L-337

\bibitem{OPERA}
Autiero D [OPERA Collaboration] Talk at Neutrino 2004 to be published in \NP B (Proceedings Supplement)\\
Guler M \etal \textsl{ Experiment proposal} CERN/SPSC 2000-028 SPSC/P318 LNGS P25/2000 

\bibitem{ICARUS}
Guler M \etal [ICARUS Collaboration] Talk at Neutrino 2004 to be published in \NP B (Proceedings Supplement) \\
Aprili P \etal CERN-SPSC-2002-027

\bibitem{CNGS}
Acquistapace G \etal CERN-98-02

\bibitem{LSND}
Athanassopoulos A \etal [LSND Collaboration] 1998 \PRL {\bf 81} 1774\\
Aguilar A \etal [LSND Collaboration] 2001 \PR D {\bf 64} 112007

\bibitem{MINIBOONE}
Stancu I \etal [MiniBooNE Collaboration] FERMILAB-TM-2207

\bibitem{PMNS}
Pontecorvo B  1957 \textsl{Sov. Phys.--JETP} {\bf 6} 429 [1957 \textsl{Zh. Eksp. Teor. Fiz.} \textbf{33} 549]\\
Maki Z, Nakagawa M and Sakata S 1962 \textsl{Prog. Theor. Phys.} \textbf{28} 870\\
Pontecorvo B 1968 \textsl{Sov. Phys.--JETP} \textbf{26}  984 [1967 \textsl{Zh. Eksp. Teor. Fiz.} \textbf{53} 1717]\\
Gribov V N and Pontecorvo B 1969 \PL B \textbf{28} 493

\bibitem{CHOOZ}
Apollonio M \etal  [CHOOZ Collaboration] 1999 \PL B {\bf 466} 415\\
Apollonio M \etal  [CHOOZ Collaboration] 2003 \EJP C {\bf 27}  331

\bibitem{Wpaper}
Anderson K \etal 2004 \textsl{White paper report on using nuclear reactors to search for a value of $\theta_{13}$} http://www.hep.anl.gov/minos/reactor13/reactor13.pdf

\bibitem{BETABEAM}
Zucchelli P 2002 \PL B \textbf{532} 166

\bibitem{NOVA}
Ambats I \etal 2004 NO$\nu$A proposal FNAL P929

\bibitem{T2K}
Itow Y 2001 KEK report 2001-4, ICRR-report-477-2001-7, TRI-PP-01-05\\
see also arXiv:hep-ex/0106019

\bibitem{BNLHS}
Diwan M V 2003 \PR D {\bf 68} 012002

\bibitem{CERN}
Blondel A \etal CERN-2004-002

\bibitem{DONINI}
Donini A \etal 2004 IFT-UAM/CSIC-04-30\\
see also arXiv:hep-ph/0406132.

\bibitem{DOUBLE-CHOOZ}
Dalnoki-Veress F for the Double-CHOOZ project arXiv:hep-ex/0406070

\bibitem{SPL}
Autin B \etal CERN-2000-12

\bibitem{UNO} 
The UNO whitepaper \textsl{Physics Potential and Feasibility of UNO} (preprint SBHEP01-3) (see also http://ale.physics.sunysb.edu/uno/)

\bibitem{mosca}
Mosca L 2004 Talk in the Physics with a Multi MegaWatt Proton Source Workshop at CERN 25-27 May 2004 http://physicsatmwatt.web.cern.ch/physicsatmwatt/

\bibitem{Meer}
van der Meer S 1961 CERN-Report No 61-7

\bibitem{nuFact134} 
Campagne J E and Cazes A 2003 CERN-NUFACT-Note-134\\ http://slac.web.cern.ch/slap/NuFact/NuFact/NFNotes.html

\bibitem{MMWPSCazes} 
Campagne J E and Cazes 2004 Poster session at the Physics with a Multi MegaWatt Proton Source Workshop at CERN 25-27 May 2004 \\
http://physicsatmwatt.web.cern.ch/physicsatmwatt/

\bibitem{JJG}
G\'omez Gadenas J J \etal 2001 IFIC/01-31 \\
see also arXiv:hep-ph/0105297

\bibitem{Mezzetto}
Mezzetto M 2003 \jpg {\bf 29} 1781-1784\\
see also arXiv:hep-ex/0302005


\bibitem{nuFact138} 
Campagne J E 2004 CERN-NUFACT-Note-138\\
 http://slap.web.cern.ch/slap/NuFact/NuFact/NFNotes.html

\bibitem{fluka}
Fasso A \etal, Proceedings of the MonteCarlo 2000 conference, Lisbon, October 26 2000, Kling A \etal - eds., Springer-Verlag Berlin, p159-164 (2001);\\
Fasso A \etal, Proceedings of the MonteCarlo 2000 conference, Lisbon, October 26 2000, Kling A \etal - eds., Springer-Verlag Berlin, p955-960 (2001)

\bibitem{MMWPSGaroby} 
Garoby R 2004 Talk in the Physics with a Multi MegaWatt Proton Source Workshop at CERN 25-27 May 2004 http://physicsatmwatt.web.cern.ch/physicsatmwatt/

\bibitem{harp}
Catanesi C \etal [HARP Collaboration] CERN-SPSC 2002/019

\bibitem{minerva} 
Drakoulakos D \etal [Minerva Collaboration] Fermilab P-938\\
see also arXiv:hep-ex/405002

\bibitem{MARS}
Mokhov N V \etal 1995 The MARS Code System User Guide, Version 13(95), Fermilab-FN-628 \\
Mokhov N V \etal 1998 MARS oed Develepments; LANL Report LA-UR-98-5716, nucl-th/9812038 v2 16 December 1998, http://www-ap.fnal.gov/MARS/

\bibitem{FLUKAprivate}
Battistoni G and Ferrari A \textsl{private communication}

\bibitem{geant}
Application Software group, Computing and Network Division \etal GEANT Description and Simulation Tool, CERN Geneva, Switzerland

\bibitem{SIMONE1}
Gilardoni S 2003 Talk at NuFact03 http://www.cap.bnl.gov/nufact03/

\bibitem{donega}
Blondel A \etal 2001 CERN-NUFACT-Note-78\\
 http://slap.web.cern.ch/slap/NuFact/NuFact/NFNotes.html

\bibitem{DONINI-2}
Donini A, Meloni D and Rigolin S 2004 J. High Energy Phys. \textbf{0406} 011\\
see also arXiv:hep-ph/0312072

\bibitem{NUANCE}
Casper D  2002\textsl{ Nucl. Phys. Proc. Suppl.} \textbf{112} 161-170\\
see also arXiv:hep-ph/0208030

\bibitem{MEZZETTONUFACT060}
Mezzetto M 2000 CERN-NUFACT-Note-60\\
 http://slap.web.cern.ch/slap/NuFact/NuFact/NFNotes.html


\bibitem{KAMLAND}
Gratta G [KamLAND Collaboration] Talk at Neutrino 2004 to be published in \NP B (Proceedings Supplement)

\bibitem{Gaisser}
Gaisser T K 1990 \textsl{Cosmic Rays and Particle Physics} Cambridge University Press 

\bibitem{picasso}
Combley F and Picasso E 1974 \textsl{Phys. Rept.} {\bf 14} 1

\bibitem{pdg}
Hagiwara K \etal 2002 \PR D {\bf 66}  010001

\bibitem{GEANT4}
Agostinelli A \etal [GEANT4 Collaboration] 2003 \NIM A \textbf{506} 250-503
\end{thebibliography}







\end{document}