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

\documentclass[11pt]{article}
%
%----------------------------------------------------------------------------
%
\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{latexsym}

%\newcommand{\st}{\tilde{t}}


%\usepackage{longtable}
% usage : \begin{lontable}{|p{1.5cm}}|c|}  : cree un tableau, la premiere 
%colone fait 1.5 cm et le texte s'arrange dedans tout \footnote{}seul, la second, 
%comme d'abb


%
%  Mes newcommand a moi
%

\newcommand{\beq}[1]{\begin{equation}\label{eq:#1}}
\newcommand{\enq}{\end{equation}}
\newcommand{\refeq}[1]{\ref{eq:#1}}
\newcommand{\und}{\underline}
%\newcommand{\dfrac}{\displaystyle\frac}
\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{.}}



% bon, celle la, elle est pas a moi, je l'ai pique a kossov et faut 
%que je comrenne comment elle marche!
%\newcommand{\eq}[2]{\begin{equation}\label{#1}#2\end{equation}}

\newcommand{\bet}[1]{\begin{table}[htb]\begin{center}\begin{tabular}{#1}}
\newcommand{\ent}[2]{\end{tabular}\caption{#1}\label{tab:#2}\end{center}\end{table}}

\newcommand{\FIG}[4]{\begin{figure}[htb]
                        \begin{center}
                        \begin{sloppypar}
                        \parbox{#4.5cm}{
                        \mbox{
                        \epsfig{file=../picts/#1.eps,width=#4cm}
                                }
                                }
                        \hspace{0.5cm}
                        \parbox{7.cm} {
                                \caption{#2}
                                \label{#3}
                                }
                        \end{sloppypar}
                        \end{center}
                    \end{figure}
                    }

\newcommand{\DBLFIG}[7]{
\begin{figure}[#1]
\begin{center}
  \begin{sloppypar}
    \parbox{8cm}{
       \mbox{\epsfig{file=../picts/#2.eps,width=8cm}}
     }
     \hspace{0.5cm}
     \parbox{8cm}{
        \mbox{\epsfig{file=../picts/#5.eps,width=8cm}}
      }
%\end{sloppypar}
%\begin{sloppypar}
      \parbox{8cm}{
        \caption{#3}
        \label{fig:#4}
      }
      \hspace{0.5cm}
      \parbox{8cm}{
         \caption{#6}
         \label{fig:#7}
      }
  \end{sloppypar}
\end{center}
\end{figure}
}

\newcommand{\DBLFIGDIFF}[6]{
\begin{figure}[pht]
\begin{center}
%  \begin{sloppypar}
    \parbox{6.cm}{
       \mbox{\epsfig{file=../picts/#1.eps,width=5.5cm}}
     }
     \hspace{0.5cm}
     \parbox{9.cm}{
        \mbox{\epsfig{file=../picts/#4.eps,width=8.5cm}}
      }
%  \end{sloppypar}
%  \begin{sloppypar}
      \parbox{6.cm}{
        \begin{flushleft}
          \caption{#2}
          \label{#3}
        \end{flushleft}
      }
      \hspace{0.5cm}
      \parbox{9.cm}{
        \begin{flushleft}
          \caption{#5}
          \label{#6}
        \end{flushleft}
      }
 % \end{sloppypar}
\end{center}
\end{figure}
}

%
%  Piqué sur le web
%
%%%% pour mettre un caption sur un environement non flotant %%%%
\makeatletter
\def\captionof#1#2{{\def\@captype{#1}#2}}
\makeatother



%
%
%----------------------------------------------------------------------------
% Config papier [Package vmargin]
%
\setpapersize{A4}
\setmarginsrb{2cm}{2cm}{2cm}{2cm}{1cm}{1cm}{1cm}{1.5cm}
% \setmarginsrb{2.5cm}{3cm}{2.5cm}{3cm}{1cm}{1cm}{1cm}{1.5cm}
%
% Explication : \setmarginsrb{a}{b}{c}{d}{e}{f}{g}{h}
%       a : marge gauche
%       b : marge haute
%       c : marge droite
%       d : marge basse
%       e : hauteur en-tete
%       f : distance en-tete -- haut du texte 
%       g : hauteur du pied de page
%       h : distance totale bas -- bas pied de page
%
%----------------------------------------------------------------------------
% Quelques definitions utiles
%
%
%
%----------------------------------------------------------------------------
%----------------------------------------------------------------------------
%

\begin{document}
%
% mefiance : le begin{document} fait bien des choses...
%
% \def\figurename{\textrm{Figure}}%
% \def\tablename{\textrm{Table}}%
\def\figurename{\textrm{Fig.}}%
\def\tablename{\textrm{Tab.}}%
\catcode`\@=11
\def\fnum@figure{\figurename~\thefigure\,}
\def\fnum@table{\tablename~\thetable\,}
\catcode`\@=12
%
%
%
%----------------------------------------------------------------------------
%----------------------------------------------------------------------------
%\selectlanguage{francais} 
\thispagestyle{empty}
%
%%
% Titre et resume...
%
\title{The $\theta_{13}$ sensitivity of the SPL-Fréjus project revisited}
\author{Jean Eric Campagne, Antoine Cazes\\
~\\
Laboratoire de l'Accélérateur Lineaire\\
Université Paris-Sud - B\^at. 200 - BP 34\\
91898 Orsay Cedex       }
\date{\today}
%
\maketitle
%
%
%\selectlanguage{francais}
%
\begin{abstract}
  \medskip 

  One of the possible tool for the search of $\theta_{13}$ is a conventional neutrino beam with a high intensity. There is a project of such a beam at CERN using the Super Proton Linac with a 4~MW beam intensity. The presented posted described the simulation of the beam as well as the computation of the neutrino flux at the Fréjus laboratory. Sensitivity is computed for the nominal energy of the SPL ($2.2$~GeV) and the gain in sensitivity is shown in the case of an increase of the SPL energy up to $8$~GeV.
  
\end{abstract}

\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} takes benefit of low energy beam to measure the $\Delta m^2$ 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 experiments in a $\bar{\nu}_e$ disappearance mode ($\sin^22\theta_{13}<0.03$) \cite{Wpaper}, and conventionnal high intensity $\nu_\mu$ beams, called Superbeam, by using the potentiality of $\nu_\mu\rightarrow\nu_e$ oscillations in the appearance mode. Ultimately, Superbeam facilities are foreseen to be extended to produce $\nu_\mu$ beams and $\bar{\nu}_\mu$ beams from muon decays, the so-called Neutrino Factory, in order to study the eventual leptonic CP violation. One of such Neutrino Complex is under study at CERN and details may be found in reference \cite{CERN}. 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 complicated to analyse due to the interplay between the different physics factors as for instance $\theta_{13}$, $\delta_{CP}$, sign$(\Delta m^2_{23})$, sign$\tan(2\theta_{23})$ \cite{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 Fig.~\ref{fig:Sbeam}.  

\begin{figure}[htb]
        \centering
                \includegraphics[scale=0.4]{../picts/cernmeg3_6-02.eps}
        \caption{Sketch of the SPL neutrino Superbeam from CERN to the Fréjus tunnel.}
        \label{fig:Sbeam}
\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 the $\pi^\pm$, $K^{\pm,0}$ and $\mu^\pm$ particles until they decay in the tunnel, the computation of the neutrino flux at the detector site, and finally the $\theta_{13}$ 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{JJG,Mezzetto,DONINI}, 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, 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. 

\begin{table}[htb]
        \centering
                \begin{tabular}{|l|c|}
    \hline
     \multicolumn{2}{|c|}{Hg target}\\
     \hline
%  Materials & Liquid Hg jet \\
      Hg jet speed & $20~m/s$ \\
      density & $13.546$ \\
      Length, radius & $30~cm\mbox{, }7.5~mm$ \\
      \hline                    
                  \end{tabular}
        \caption{Liquid mercury jet parameter}
        \label{tab:targ}
\end{table}

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 Tab.~\ref{tab:targ}). The pencil like simulated proton beam is composed with  $10^6$ mono energetic protons, and 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}.

\bet{|c||c|c|c|c|c|c|c|c|c|c|c|}
\hline
$E_k$ (GeV) & p & n & $\gamma$ & $e^-$ & $e^+$ & $\pi^+$ & $\pi^-$ & $\mu^+$ & $\mu^-$ & $K^+$ & $K^0$ \\
\hline
\hline
2.2 & $1.4$     &  $17$   & $5.0$  &  $0.17$   &  $0.08$    &  $0.24$   &  $0.18$ & $4.~10^{-4}$ & $1.~10^{-4}$ & $7.~10^{-4}$ & $6.~10^{-4}$ \\
\hline
3.5 &$1.8$    &  $17$  & 7.0  &  0.28   &  0.15    &  0.41   &  0.37 & $10.~10^{-4}$ & $3.~10^{-4}$ & $35.~10^{-4}$ & $30.~10^{-4}$ \\
\hline
$4.5$ &  $2.3$  &  $17$ & $7.7$ &  $0.35$  & $0.21$   &  $0.57$ & $0.39$&  $11.~10^{-4}$ & $3.3~10^{-4}$  &  $93.~10^{-4}$ &  $68.~10^{-4}$       \\
\hline
$8$   &  $3.1$  &  $17$ & $11$ & $0.63$  & $0.41$   &  $1.0$  &  $0.85$ &  $30.~10^{-4}$ &  $9.5~10^{-4}$ &  $413.~10^{-4}$ &    $340.~10^{-4}$         \\
\hline
\ent{Average number of the most relevant secondary particles exiting the $30$~cm long, $1.5$~cm diameter mercury target per incident proton (FLUKA). Note that the $K^-$ are at the level of $10^{-5}$ per proton.}{nbPart}

Particle production yields are summarized in Tab.~\ref{tab:nbPart}. The level of the secondary proton and neutron emissions induces important radiation damages and power dissipation problems on the horns which have been addressed in reference \cite{nuFact134}, and which will require specific R\&D effort \cite{nuFact134}. At $2.2$~GeV, kaon yields are very low, but it has a dramatic energy dependence as further studied in Sec.~\ref{sec:kaon}. It is worth to mention that the numbers in Tab.~\ref{tab:nbPart} are not to be taken at face value, 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}.

The cross section uncertainties are the main source of discrepancy between simulator programs. Some comparisons between FLUKA and MARS \cite{MARS} have allready been presented in the same context \cite{nuFact134}. The energy distribution of the pions exiting the target, computed with the two simulator programs (FLUKA in blue dashed line and MARS in red full line), is shown on Fig.~\ref{fig:comp1}. 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, MARS and FLUKA are in better agreement for the energy spectrum computed at the entrance of the decay tunnel (Fig.~\ref{fig:comp2}). Moreover, the horns are designed to focus the high energy part of the spectrum (see Sec.~\ref{sec:horn}), and therefore, the discrepancy  at low energy between MARS and FLUKA does not matter too much for the present application.

\DBLFIG{ht}{compMARSmom_NB}{Momentum distribution of the pion exiting the target simulated by FLUKA (in blue dashed line) and by MARS (in red full line).}{comp1}{compCorneMARS_NB}{Same as figure \ref{fig:comp1} but after pions have been tracked through the magnetic fields of the horns.}{comp2}
%
\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 Tabs.~\ref{tab:BRkp}, \ref{tab:BRk0S} and \ref{tab:BRk0L} in appendix~\ref{sec:kaons}.

The target simulation described in Sec.~\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. Numbers of produced $K^+$ at different proton beam energies are presented on Fig.~\ref{fig:kaons}. Numbers of $K^-$ is also plotted on the same figure (dashed line) but they are almost 40 times less numerous than the $K^+$. For comparison,  the numbers of $\pi^+$ (solid line) and $\pi^-$ (dashed line) produced in the same conditions are presented on Fig.~\ref{fig:pions}. Pion production rate is about two orders of magnitude greater than the kaon production rate.
% ET LES KO???? J.E 25/7/04
\DBLFIG{h}{kaons}{kaon production as a function of the incident proton beam energy for $500~000$ incident protons}{kaons}{piplus}{Same as figure \ref{fig:kaons} for the $\pi^\pm$ production}{pions}

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 experience 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 Tab.~\ref{tab:specif}.

\begin{figure}[htb]
        \centering
                \includegraphics[width=0.75\textwidth]{../picts/corneSchema.eps}
        \caption{Design of the Horn (in blue) and the Reflector (in red). The target is located inside the cylindrical part of the Horn (in grey).}
        \label{fig:plan}
\end{figure}

%\FIG{corne}{}{fig:plan}{9}

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>=240$~MeV$/400$~MeV. This explains the Horn design (Fig.~\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 a cooling system, 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.
\begin{table}[ht]
\begin{center}
\begin{sloppypar}
\parbox{8.0cm}{
\begin{center}
  \begin{tabular}{|l|c|}
  \hline
  \multicolumn{2}{|c|}{inner horn}\\
  \hline
  neck inner radius & $3.7$~cm \\
  neck length & $40$~cm \\
  end cone inner radius & $15.7$~cm \\
  outer radius & $20.3$~cm \\
  total length & $180$~cm \\
  Alu thickness  & $3$~mm \\
  %cooling water thickness & $2~mm$ \\
  %double skin thickness & $2~mm$ \\
  Peak current & $300$~kA \\
  Frequency & $50$~Hz \\
  \hline
  \end{tabular}
\end{center} }
%\hspace{0.5cm}
\parbox{8.0cm}{
\begin{center}
  \begin{tabular}{|l|c|}
  \hline
  \multicolumn{2}{|c|}{reflector}\\
  \hline
  neck inner radius & $20.3$~cm \\
  end cone inner radius & $35.7$~cm \\
  outer radius & $40$~cm \\
  total length & $200$~cm \\
  Peak current & $600$~kA \\
  Frequency & $50$~Hz \\
  \hline
  \end{tabular}
\end{center} }
\end{sloppypar}
\caption{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.}
\label{tab:specif}
\end{center}
\end{table}

%\subsection{Study of the horns designs in function of the energy}

The shape of the horns conductors is a crucial point since it determines the energy spectrum of the neutrino at the detector site. For a $\theta_{13}$ driven oscillation $\nu_\mu \rightarrow \nu_e$, 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. 

Figure \ref{fig:focusing} represents, in the $r-z$ plane, the trajectories of pions emitted at $z=r=0$ with different zenith angles. Particles are deviated\footnote{The magnetic field inside the aluminum conductors themselves is neglected as the material thickness is of the order of 3~mm.} by a toroidal magnetic field for $r>4$~cm:
$$
B_\phi(r) = \frac{\mu_0 I}{2\pi r}
$$
where $I$ is the peak current, which is set to $300$~kA for the Horn ($4~\mathrm{cm}<r<20~\mathrm{cm}$) and to $600$~kA for the Reflector ($20~\mathrm{cm}<r<40~\mathrm{cm}$). The exact profile of the conic part of the horns is determined such that the relevant pions exit parallel to the beam axis \cite{nuFact138}. The Horn cylindrical part is $40$~cm long while the conic part is $80$~cm long with an opening angle of $8.4^\circ$. The Reflector cylindrical part is adjusted to the Horn external cylinder, and the conic part is $80$~cm long with an opening angle of $12.8^\circ$. Another design have been studied with a Reflector having a $60$~cm radius. It shows an improvement by $35\%$ of the fluxes, but no improvement in the $\theta_{13}$ sensitivity, therefore, it is no longer used.

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 with an energy spectrum which of course is affected by a proton beam energy change.    

\begin{figure}[htb]
        \centering
                \includegraphics[width=0.6\textwidth]{../picts/corne600hr.eps}
          \caption{Pion focusing optimization. Lines represent different trajectories of pion having a $600$~MeV momentum emitted at $r=0$, $z=0$ with different angles. The red line joins the points where the trajectories are parallel to the $z$ axis and the dashed line represents the chosen design of the horn, and two possible design of the reflector($20~cm<r-1<40~cm$ and $20~cm<r-1<60~cm$}
                \label{fig:focusing}
\end{figure}
%
\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 used in reference \cite{donega} 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 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/K/\mu$ particle. 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 compare 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 Tab.~\ref{tab:nbPart}).

The probability computation is presented in the appendix. 

\begin{table}[htb]
        \centering
                \begin{tabular}{|c|c|}
                        \hline
                        Beam energy  & number of proton \\
                          (GeV)      & per year ($10^{23}$ p.o.t/y) \\
                        \hline
                        2.2 & 1.1 \\
                        3.5 & 0.7 \\
                        4.5 & 0.56 \\
                        6.5 & 0.4 \\
                        8   & 0.3 \\
                  \hline
                \end{tabular}
        \caption{Number of protons on target for different beam energy at 4~MW constant power.}
        \label{tab:proton}
\end{table}

\subsection{Validation of the algorithms}

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$). A full GEANT simulation of the event (decays included) is performed and presents the advantage to keep all the informations of the neutrino available for further studies. Such method 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 prohibitive CPU time consuming.
\begin{figure}[htb]
        \centering
                \includegraphics[width=0.50\textwidth]{../picts/compGeantDonega_NB.eps}
        \caption{Comparison between the probability method (see Sec.~\ref{sec:algo}) and the full GEANT simulation method for the $\nu_\mu$ and $\bar{\nu}_\mu$ fluxes.}
        \label{fig:compGeantDonega}
\end{figure}

The $\nu_\mu$ and $\bar{\nu}_\mu$ fluxes are displayed on Fig.~\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, but 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 (Sec.~\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. For a $2.2$~GeV beam of $4$~MW, the number of protons is $1.1~10^{16}$~p.o.t/s. 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. In the following, the up to now nominal $2.2$~GeV proton beam is used, if it is not otherwise mentioned.

Three origins are identified in the composition of each neutrino flux (see Tab.~\ref{tab:BRkp}-\ref{tab:BRk0L} in appendix \ref{sec:kaons} for the kaon cases):
\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 as well as the total neutrino fluxes for the four neutrino species are presented on Fig.~\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 neutrino of pion decays, but a tail at high energy (above $500$~MeV) 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 less numerous than the muon neutrino. 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^-$ decay, while the  $\nu_e$ flux is dominated by the $\mu^+$ decays.

On Fig.~\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 Tab.~\ref{tab:nbPart}).

On Figs.~\ref{fig:flux45p} and \ref{fig:flux80p}, one observes the evolution of Fig.~\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 Figs.\ref{fig:flux45m} and \ref{fig:flux80m}. One clearly notices the increase of the kaon induced neutrino contents as the beam energy grows.

\begin{figure}[p]
        \centering
        \includegraphics[scale=0.6]{../picts/flux22p.eps}
        \caption{Neutrino fluxes $100$~km from the decay region and with the collection elements 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. Solid line is the    total flux. Dashed line is the contribution from primary pions and their muons. Dotted line are the neutrino from the charged kaon decay chains, and dotted-dashed line are for the $K^0$ decay chain contribution.}
        \label{fig:flux22p}
\end{figure}

\begin{figure}[p]
        \centering
        \includegraphics[scale=0.6]{../picts/flux22m.eps}
        \caption{Same remarks as for Fig.~\ref{fig:flux22p} but the collection elements are  focusing negative charges.}
        \label{fig:flux22m}
\end{figure}

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

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

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


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

\begin{table}[htb]
        \centering
                \begin{tabular}{|c|c|c|c|c|c|c|c|c|}
                        \hline
                        Settings  & \multicolumn{2}{c|}{$\nu_\mu$} & \multicolumn{2}{c|}{$\nu_e$} 
                                                                &       \multicolumn{2}{c|}{$\bar{\nu}_\mu$} & \multicolumn{2}{c|}{$\bar{\nu}_e$} \\
                                                                \cline{2-9}
                                                                & $+$ & $-$ & $+$ & $-$ & $+$ & $-$ & $+$ & $-$ \\ \hline
(1): $2.2$~GeV  
        & $7.6$          & $4.5\,10^{-1}$ 
        & $3.2\,10^{-2}$ & $2.4\,10^{-3}$
        & $3.1\,10^{-1}$ & $5.8$  
        & $1.1\,10^{-3}$ & $1.6\,10^{-2}$ \\              \hline

(1): $3.5$~GeV  
        & $10.0$          & $8.7\,10^{-1}$ 
        & $4.4\,10^{-2}$  & $6.1\,10^{-3}$
        & $7.0\,10^{-1}$ & $8.5$  
        & $2.8\,10^{-3}$ & $2.2\,10^{-2}$ \\              \hline

(1): $4.5$~GeV  
        & $10.9$          & $1.1$ 
        & $5.0\,10^{-2}$  & $9.6\,10^{-3}$
        & $6.8\,10^{-1}$ & $6.7$  
        & $3.9\,10^{-3}$ & $1.8\,10^{-2}$ \\              \hline

(1): $6.5$~GeV  
        & $10.4$          & $1.4$ 
        & $6.4\,10^{-2}$  & $2.0\,10^{-2}$
        & $9.6\,10^{-1}$ & $7.0$  
        & $9.2\,10^{-3}$ & $2.5\,10^{-2}$ \\              \hline

(1): $8.0$~GeV  
        & $9.7$          & $1.5$ 
        & $6.6\,10^{-2}$  & $2.3\,10^{-2}$
        & $1.1$           & $7.1$  
        & $1.1\,10^{-2}$ & $2.8\,10^{-2}$ \\
                  \hline\hline

(2): $4.5$~GeV  
        & $7.0$          & $1.0$ 
        & $2.7\,10^{-2}$  & $9.1\,10^{-3}$
        & $5.0\,10^{-1}$ & $5.2$  
        & $3.4\,10^{-3}$ & $1.1\,10^{-2}$ \\              \hline

(3): $4.5$~GeV  
        & $12.3$          & $1.2$ 
        & $7.2\,10^{-2}$  & $1.0\,10^{-2}$
        & $7.6\,10^{-1}$ & $7.5$  
        & $4.1\,10^{-3}$ & $2.6\,10^{-2}$ \\              \hline\hline

(4): $4.5$~GeV  
        & $13.2$          & $1.5$ 
        & $6.9\,10^{-2}$  & $1.4\,10^{-2}$
        & $9.0\,10^{-1}$  & $8.1$  
        & $5.6\,10^{-3}$  & $2.7\,10^{-2}$ \\             \hline

                \end{tabular}
        \caption{Integral of the different species fluxes with different settings. Numbers are in $10^{13}/100\mathrm{m}^2/y$ for positive focusing ($+$) and negative focusing ($-$). 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.}
        \label{tab:speciesfluxes}
\end{table}

On Tab.~\ref{tab:speciesfluxes} are reported the integral of the fluxes when one modifies the decay tunnel geometry as for instance the length or the radius, as well as the beam kinetic energy (keeping the beam power constant).  Switching the radius from $1$~m to $1.5$~m (doubling the surface basically) has increased the signal like fluxes ($\nu_\mu$, $\bar{\nu}_\mu$) by a factor $20\%$ but also in paralell the background like fluxes ($\nu_e$, $\bar{\nu}_e$) by a factor $38\%$. Increasing the 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$. 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 Sec.~\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). The indication that a new beam energy base line arround $4.5$~GeV may be proposed, replacing the previous $2.2$~GeV base line design. This proposal will be reenforced in Sec.~\ref{sec:results}.
%
\section{Sensitivity computation ingredients}
%
The sensitivity to $\theta_{13}$ is computed for an $\nu_\mu\rightarrow\nu_e$ appearance experiment. An analysis program described in reference \cite{Mezzetto} has been used for such sensitivity computation. See Tab.~\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^0$ 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 default $2\%$ value and also $5\%$ and $10\%$. 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 Sec.~\ref{sec:horn} for discussion on Horn conductor shape determination).  

\begin{table}[htb]
        \centering
                \begin{tabular}{|l|l|}
                        \hline
                        $\Delta m^2_{12} = 8.2 10^{-5}~\mathrm{eV}^2$ & $\sin^22\theta_{12} = 0.82$ \\
                        $\sigma(\Delta m^2_{12}) = 0.5 10^{-5}~\mathrm{eV}^2$ & $\sigma(\sin^22\theta_{12}) = 9\%$ \\ \hline
                        $\Delta m^2_{23} = 2.5 10^{-3}~\mathrm{eV}^2$ & $\sin^22\theta_{23} = 1.$ \\
                        $\sigma(\Delta m^2_{23}) = 10^{-4}~\mathrm{eV}^2$ & $\sigma(\sin^22\theta_{23}) = 1\%$ \\
                  \hline
                  \hline
                  $L_T = 20$~m & $R_T = 1$~m \\ \hline
                  $M=440$~kT   & $\epsilon_{syst}=2\%$ \\
                  \hline 
                \end{tabular}
        \caption{Default user parameters used to compute the sensitivity curves \cite{Mezzetto}. The quoted errors for the $(12)$ and the $(23)$ parameters 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}.}
        \label{tab:param}
\end{table}

The running time scenario has been fixed to 5 years focusing $\pi^+$ for the $\theta_{13}$ sensitivity studies. This scenario gives comparable $\theta_{13}$ sensitivity than a scenario defined as 2 years focusing $\pi^+$ followed by 8 years focusing $\pi^-$, which is more appropriate for $\delta_{CP}$ sensitivity studies.
%
\section{Results}
\label{sec:results}
The $\theta_{13}$-sensitivity is computed for $\theta_{13} = 0^\circ$ and $\delta_{CP} = 0^\circ$ if not explicitly mentioned. The other default parameters are listed in Tab.~\ref{tab:param}.

Tab.~\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{Mezzetto} 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 prompt $\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.
%
\begin{table}[htb]
        \centering
                \begin{tabular}{|p{4cm}|c|c|c|c|c|}
                        \cline{2-6}
                        \multicolumn{1}{c|}{} & $2.2$~GeV & $3.5$~GeV & $4.5$~GeV & $6.5$~GeV & $8$~GeV \\
                        \hline
                        non oscillated $\nu_\mu$       & $36917.2$  & $60968.6$  & $73202.4$  & $78023.9$ & $76068.2$ \\
                        \hline
                        oscillated $\nu_e$             & $43.0$  & $59.7$  & $63.7$  & $60.6$ & $55.9$ \\
                        \hline
                        beam $\nu_e$, $\bar{\nu}_e$     & $165.9$ & $225.2$ & $244.9$ & $299.8$ & $309.7$ \\
                        \hline
                        other background: $\pi^0$, $\nu_\mu$-elast., $\mu/e$-missId. & $70.3$  & $104.7$ & $126.6$ & $147.6$ & $151.3$ \\
                        \hline
                        Significance                  &  $2.50$ & $2.91$  & $2.93$  & $2.57$ & $2.34$ \\
                        \hline
                \end{tabular}
        \caption{Number of events for 5 years positive focusing scenario with default parameters of Tab.~\ref{tab:param}. The significance parameter is defined by Eq.~\ref{eq:significance}.}
        \label{tab:nbvsE}               
\end{table}

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 Fig.\ref{fig:sensi45} for $4.5$~GeV proton beam kinetic energy. The comparison between the contours at $90\%$~CL but with $2.2$~GeV, $3.5$~GeV, $4.5$~GeV and $8$~GeV is shown on Fig.~\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 Sec.~\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 Tabs.~\ref{tab:thvsE}, \ref{tab:thvseps}. The influance of the systematical level is presented on Fig.~\ref{fig:compEpsSyst}. One should keep in mind that $\epsilon_{syst}$ is considered as an ultimate goal.     
%
\begin{figure}[htb]
        \centering
                \includegraphics{../picts/sensi45.eps}
        \caption{Sensitivity contours obtained with a SPL energy of $4.5$~GeV and default parameters of Tab.~\ref{tab:param}.}
        \label{fig:sensi45}
\end{figure}
%
\begin{figure}[htb]
        \centering
                \includegraphics{../picts/compSensi.eps}
        \caption{Comparison of 90\% sensitivity contours obtained with SPL energies of ($2.2$, $3.5$, $4.5$, $8$)~GeV and default parameters of Tab.~\ref{tab:param}.}
        \label{fig:compSensi}
\end{figure}


\begin{table}[htb]
        \centering
                \begin{tabular}{|l|c|c|c|c|c|}
                        \hline
                        \backslashbox{$L_T$}{$E_k(proton)$} & $2.2$~GeV & $3.5$~GeV & $4.5$~GeV & $6.5$~GeV & $8$~GeV  \\
                        \hline
                        $10$m & $1.10$ & $0.96$ & $1.05$ & $1.10$ & $1.20$ \\
                        \hline
                        $20$m & $1.20$ & $0.96$ & $0.91$ & $1.05$ & $1.15$ \\
                        \hline
                        $40$m & $1.26$ & $1.00$ & $1.00$ & $1.10$ & $1.20$ \\
                        \hline
                \end{tabular}
        \caption{Minimum $\sin^22\theta_{13}\times 10^3$ observable at $90\%$ CL. Other parameters are fixed to default values (Tab.~\ref{tab:param}).}
        \label{tab:thvsE}
\end{table}



\begin{table}[htb]
        \centering
                \begin{tabular}{|l|c|c|c|c|c|}
                        \hline
                        \backslashbox{$\epsilon_{syst}$}{$E_k(proton)$} & $2.2$~GeV & $3.5$~GeV & $4.5$~GeV & $6.5$~GeV & $8$~GeV \\
                        \hline
                        $2\%$ & $1.20$ & $0.96$ & $0.91$ & $1.10$ & $1.20$ \\
                        \hline
                        $5\%$ & $1.51$ & $1.26$ & $1.26$ & $1.51$ & $1.66$ \\
                        \hline
                        $10\%$ & $2.40$ & $2.19$ & $2.29$ & $2.75$ & $3.16$ \\
                        \hline
                \end{tabular}
        \caption{minimum $\sin^22\theta_{13}\times 10^3$ observable at $90\%$ CL. Other parameters are fixed to default values (Tab.~\ref{tab:param}).}
        \label{tab:thvseps}             
\end{table}

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


There are also uncertainties on the minimum $\sin^22\theta_{13}$ value that might be reached in a $\nu_\mu \rightarrow \nu_e$ experiment which are due to the $\mathrm{sign}(\Delta m^2_{23})$ and the $\delta_{CP}$ ambiguities. On Fig.~\ref{fig:compDelta_thVSdm} and Tab.~\ref{tab:sign} are presented the results of these ambiguities making use of the oscillation probability symmetry $P(\Delta m^2_{23},\delta_{CP}) =  P(-\Delta m^2_{23},\pi - \delta_{CP})$. Other ambiguities coming from the sign$(\tan(2\theta_{23}))$ ignorance also take place as studied in reference \cite{DONINI}. From Fig.~9 of this reference, we estimate a 30\% effect on $\sin^2(2\theta_{13})$ sentivity 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}
 
\begin{table}[htb]
        \centering
                \begin{tabular}{|c|c|c|}
                        \hline
                        \backslashbox{sign$(\Delta m^2_{23})$}{$\delta_{CP}$} & $0^\circ$ & $90^\circ$ \\
                        \hline
                        $+$ & $1.20$ & $11.48$ \\
                        \hline
                        $-$ & $1.15$ & $0.44$ \\
                        \hline
                \end{tabular}
        \caption{Minimum $\sin^22\theta_{13}\times 10^3$ observable at $90\%$ CL. Other parameters are fixed to default values (Tab.~\ref{tab:param}).}
        \label{tab:sign}                
\end{table}


The combined $\sin^22\theta_{13}$ and $\delta_{CP}$ sensibility is presented on Fig.~\ref{fig:delta45} for the 5 years positive focusing scenario to appreciate in an other way the above mentioned $\delta_{CP}$ ambiguity. 
\begin{figure}[htb]
        \centering
                \includegraphics{../picts/delta45.eps}
        \caption{Sensitivity contours obtained with a SPL energy of $4.5$~GeV and default parameters of Tab.~\ref{tab:param} and 5 years positive focusing scenario.}
        \label{fig:delta45}
\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{JJG,Mezzetto,DONINI}. The corresponding combined sensitivity contours are presented in Fig.~\ref{fig:delta45FocPM}. With this kind of mixed focusing scenario, $\sin^22\theta_{13} < 10^{-3}$ independently of $\delta_{CP}$.
\begin{figure}[htb]
        \centering
                \includegraphics{../picts/delta45FocPM.eps}
        \caption{Same as Fig.~\ref{fig:delta45} but for 2 years positive focusing and 8 years negative focusing scenario.}
        \label{fig:delta45FocPM}
\end{figure}
The comparison of the two scenarios modifying the SPL beam energy is presented on Figs.~\ref{fig:compDelta5pos}, \ref{fig:compDeltaFocPM}. It is noticeable that for $\delta_{CP}\approx 50^\circ$ ($\Delta m^2_{23}>0$) one gets better sensitivity with a $2.2$~GeV SPL beam than a $4.5$~GeV SPL beam. This is due to 
\begin{figure}[htb]
        \centering
                \includegraphics{../picts/compDelta5pos.eps}
        \caption{Same conditions as Fig.~\ref{fig:delta45} but at different SPL beam energy.}
        \label{fig:compDelta5pos}
\end{figure}
\begin{figure}[htb]
        \centering
                \includegraphics{../picts/compDeltaFocPM.eps}
        \caption{Same conditions as Fig.~\ref{fig:delta45FocPM} but at different SPL beam energy.}
        \label{fig:compDeltaFocPM}
\end{figure}
 
%
\section{Summary and outlook}
A complete chain of simulation has been set up for the SPL super beam project. It includes neutrino production from kaon decay, which is important to test higher SPL energy scenario.
A simulation has been performed with the nominal energy of the SPL and also for an kinetic energy of $3.5$~GeV. The gain in the sensitivity is $30\%$ in this last scenario with a moderate change  in the horn profiles, which leads to a potential sensitivity to $\sin^22\theta_{13} = 2.10^{-3}$ at $90\%$ CL. But, it could be improved since many parameters have not yet been optimized: reflector/horn design, decay tunnel length, $\pi^+$/$\pi^-$ focusing time...
%
\section*{Acknoledgments}
The authors are very grateful to thanks M. Mezzetto to have express 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
\section{Decay probability computations}
This appendix contains the probability formulas and the algorithms used in the flux computation (see Sec.~\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 Fig.~\ref{fig:pionDecay}). This gives:
\beq{probaPi}
\mathcal{P}_\pi = \frac{1}{4\pi}\frac{A}{L^2}\frac{1-\beta^2}{(\beta\cos\alpha-1)^2}
\enq
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}[htb]
        \centering
                \includegraphics{../picts/pionDecay.eps}
        \caption{Pion decay in the tunnel frame. To reach the detector, $\delta = -\alpha$ is needed.}
        \label{fig:pionDecay}
\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}:
\beq{probaMu}
\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]
\enq 
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 Tab.~\ref{tab:Function}. The sign in front of $\Pi_\mu^T$ in Eq.~\refeq{probaMu} is $(-)$ for the $\nu_\mu$ and $(+)$ for the $\bar{\nu}_\mu$, respectively.

\begin{table}[htb]
        \begin{center}
                \begin{tabular}{|c||c|c|}
                \hline
                      &  $f_0(x)$ & $f_1(x)$ \\
                \hline
                \hline
                $\nu_\mu$ & $2x^2(3-2x)$ & $2x^2(1-2x)$ \\
                \hline
                $\nu_e$   & $12x^2(1-x)$ & $12x^2(1-x)$ \\
                \hline
                \end{tabular}
        \end{center}
        \caption{Flux function in the muon rest frame \cite{Gaisser}.}
        \label{tab:Function}
\end{table}

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}:
\beq{pola}
\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 
\enq
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 tables \ref{tab:BRkp}, \ref{tab:BRk0S} and \ref{tab:BRk0L}.

There is a very small amount of kaons produced (Sec.~\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 (Sec.~\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 (Eq.~\refeq{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:
\beq{probaL}
\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}\virg
\enq 
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 Eq.~\refeq{probaPi} is applied to the produced neutrino. In case of a muon, it is treated as explained in Sec.~\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.

\begin{table}
\begin{center}
\begin{sloppypar}
\parbox{5.cm}{
\begin{center}
\mbox{
\begin{tabular}{|c|c|}
\hline
\multicolumn{2}{|c|}{$K^\pm$} \\
\hline
$\mu^\pm\nu_\mu$ & $63,51\%$ \\
\hline
$\pi^\pm\pi^0$ & $21,17\%$ \\                   
\hline
$\pi^\pm\pi^+\pi^-$ & $5,59\%$ \\
\hline
$e^\pm\nu_e\pi^0$ & $4,82\%$ \\                         
\hline  
$\mu^\pm\nu_\mu\pi^0$ & $3,18\%$ \\
\hline
$\pi^\pm\pi^0\pi^0$ & $1,73\%$ \\                       
\hline
\end{tabular}
}
\end{center}
}
\hspace{0.5cm}
\parbox{4cm}{
\begin{center}
\mbox{
\begin{tabular}{|c|c|}
\hline
\multicolumn{2}{|c|}{$K^0$ short} \\
\hline
$\pi^+\pi^-$ & $68,61\%$ \\
\hline
$\pi^0\pi^0$ &  $31.39\%$\\
\hline
\end{tabular}
}
\end{center}
}
\hspace{0.5cm}
\parbox{5.cm}{
\begin{center}
\mbox{
\begin{tabular}{|c|c|}
\hline
\multicolumn{2}{|c|}{$K^0$ long} \\
\hline
$\pi^-e^+\nu_e$ & $19,35\%$ \\
\hline
$\pi^+e^-\bar{\nu}_e$ & $19,35\%$ \\                    
\hline
$\pi^-\mu^+\nu_\mu$ & $13,5\%$ \\
\hline
$\pi^+\mu^-\bar{\nu}_\mu$ & $13,5\%$ \\                         
\hline  
$\pi^0\pi^0\pi^0$ & $21,5\%$ \\
\hline
$\pi^+\pi^-\pi^0$ & $12,38\%$ \\
\hline
\end{tabular}
}
\end{center}
}
%\end{sloppypar}
%    \begin{sloppypar}
      \parbox{5.cm}{
        \caption{Charged kaons decay channels and branching ratios \cite{pdg}}
              \label{tab:BRkp}
      }
      \hspace{0.5cm}
      \parbox{5.cm}{
        \caption{$K^0_S$ decay channels and branching ratios \cite{pdg}}
        \label{tab:BRk0S}
      }
      \hspace{0.5cm}
      \parbox{5.cm}{
        \caption{$K^0_L$ decay channels and branching ratios \cite{pdg}}
        \label{tab:BRk0L}
      }
\end{sloppypar}
\end{center}
\end{table}



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

\bibitem{LSND}
A. Athanassopoulos at al. [LSND Collaboration], Phys. Rev. Lett. {\bf 81} (1998) 1774;\\
A. Aguilar et al. [LSND Collaboration], Phys. Rev. D {\bf 64} (2001) 112007.

\bibitem{MINIBOONE}
I. Stancu et al. [MiniBooNE Collaboration], FERMILAB-TM-2207.

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

\bibitem{SKNU04}
E.~Kearns [SuperKamiokande Collaboration], Talk at Neutrino 2004, 14-19 June 2004 Paris, to be published in Nuclear Physics B (Proceedings Supplement).

\bibitem{K2KNU04}
T.~Nakaya [K2K Collaboration], Talk at Neutrino 2004, 14-19 June 2004 Paris, Nuclear Physics B (Proceedings Supplement). 

\bibitem{MINOS}
M.~Thomson  [MINOS Collaboration], Talk at Neutrino 2004, 14-19 June 2004 Paris, Nuclear Physics B (Proceedings Supplement); \\
P.~Adamson et al., NuMi-L-337.

\bibitem{KAMLAND}
G. Gratta [KamLAND Collaboration], Talk at Neutrino 2004, 14-19 June 2004 Paris, Nuclear Physics B (Proceedings Supplement).

\bibitem{OPERA}
D.~Autiero [OPERA Collaboration], Talk at Neutrino 2004, 14-19 June 2004 Paris, Nuclear Physics B (Proceedings Supplement); \\
M.~Guler et al., Experiment proposal, CERN/SPSC 2000-028, SPSC/P318, LNGS P25/2000 

\bibitem{ICARUS}
M.~Guler et al. [ICARUS Collaboration], Talk at Neutrino 2004, 14-19 June 2004 Paris, Nuclear Physics B (Proceedings Supplement); \\
P.~Aprili et al.  , CERN-SPSC-2002-027.

\bibitem{CNGS}
G.~Acquistapace et al., CERN-98-02.

\bibitem{CHOOZ}
M.~Apollonio et al.  [CHOOZ Collaboration], Phys.\ Lett.\ B {\bf 466} (1999) 415;\\
M.~Apollonio et al.  [CHOOZ Collaboration], Eur. Phys. J. C {\bf 27} (2003) 331.

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


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

\bibitem{SPL}
B. Autin et al, CERN 2000-12.

\bibitem{CERN}
A. Blondel et al., CERN-2004-002.

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

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

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

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

\bibitem{harp}
G. Catanesi et al. [HARP Collaboration], CERN-SPSC 2002/019.

\bibitem{minerva} 
D. Drakoulakos et al. [Minerva Collaboration], Fermilab P-938, arXiv:hep-ex/405002.

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

\bibitem{MMWPSCazes} 
J.E. Campagne, A. Cazes, poster session at the Physics with a Multi MegaWatt Proton Source Workshop held at CERN, 25-27 May 2004. 

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

\bibitem{geant}
Application Software group, Computing and Network Division \& al. --
\newblock GEANT Description and Simulation Tool, CERN Geneva, Switzerland

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

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

\bibitem{donega}
M. Donegà, CERN-NUFACT-Note-78 (2001).

\bibitem{Gaisser}
T.K. Gaisser, Cosmic Rays and Particle Physics, Cambridge University Press (1990) 

\bibitem{picasso}
F. Combley and E. Picasso, Phys. Rept. {\bf 14} (1974).

\bibitem{pdg}
K.~Hagiwara et al., Phys. Rev. D {\bf 66} (2002)  010001.

\bibitem{Mezzetto}
M.~Mezzetto, J. Phys. G {\bf 29} (2003) 1781-1784;\\
arXiv:hep-ex/0302005.

\bibitem{JJG}
J. J. G\'omez Gadenas et al., IFIC/01-31, arXiv:hep-ph/0105297.

\bibitem{DONINI}
A. Donini et al., IFT-UAM/CSIC-04-30, arXiv:hep-ph/0406132.

\bibitem{DONINI-2}
A. Donini, D. Meloni and S. Rigolin, JHEP \textbf{0406} (2004) 011; arXiv:hep-ph/0312072.
 
\bibitem{mosca}
L.~Mosca, Talk in the Physics with a Multi MegaWatt Proton Source Workshop held at CERN, 25-27 May 2004, (http://physicsatmwatt.web.cern.ch/physicsatmwatt/)

\bibitem{NUANCE}
D.~Casper, Nucl. Phys. Proc. Suppl. \textbf{112} (2002) 161-170;\\
arXiv:hep-ph/0208030.

\end{thebibliography}







\end{document}