text:
\section{The NuMI{} Beam}
\label{sec:beam}
\begin{figure*}
  \centering
  \includegraphics[width=\textwidth,keepaspectratio,clip]{fig/numi-beam-figure.eps}
  \caption{ \label{fig:numi} Plan and elevation views of the NuMI{} beam facility.  A proton beam is directed onto a target, where the secondary pions and kaons are focused into an evacuated decay volume via magnetic horns.  Ionization chambers at the end of the beam-line measure the secondary hadron beam and tertiary muon beam.}
\end{figure*}

The NuMI{} Neutrino Beam~\cite{numitdr,Kopp:2005bt} generates neutrinos mainly from the decays of pions and kaons produced in the NuMI{} target, with a smaller contribution from muon decays. A schematic diagram of the NuMI{} beam-line is shown in Fig.~\ref{fig:numi}.  Protons of \unit[120]{\ensuremath{\mathrm{GeV}/\mathit{c}}{}} momentum are extracted from the Main Injector (MI) accelerator in a \unit[10]{$\mu$s} spill, bent downward by \unit[58]{mrad} to point at the Far Detector, and impinged upon the NuMI{} hadron production target.  Positioned downstream of the target, two toroidal magnets called ``horns'' sign-select and focus the secondary mesons from the target, as shown in Fig.~\ref{fig:targ-horns}.  The mesons are directed into a \unit[675]{m}{} long evacuated volume, where they may decay to muons and neutrinos.  At the end of the decay volume, a beam absorber stops the remnant hadrons. The absorber is followed by approximately \unit[240]{m} of un-excavated rock which stops the tertiary muons, leaving only neutrinos. Figure~\ref{fig:numi-spectra} shows \ensuremath{\nu_{\mu}} charged-current{} spectra from three configurations of the NuMI{} beam.  The following section provides a brief overview of the beam-line and instrumentation.

\subsection{Beam Hardware and Performance}
\subsubsection{Primary Proton Beam}

The Main Injector accepts batches of protons from the \unit[8]{\ensuremath{\mathrm{GeV}/\mathit{c}}{}} Booster accelerator, which are then accelerated to \unit[120]{\ensuremath{\mathrm{GeV}/\mathit{c}}{}}. In most MI cycles seven batches were accelerated. Protons were removed from the MI ring in two increments, with a pair of batches being sent to the Antiproton Accumulator for the Tevatron program, while the remaining five batches were directed into the NuMI{} primary-proton line and transported \unit[350]{m} to the NuMI{} target.  The NuMI{} extractions typically contained a total of $\unit[2.1\times10^{13}]{protons}$ with a cycle time of $2.2-\unit[2.4]{s}$. Extractions of up to six batches and \unit[$2.5\times10^{13}$]{protons} were achieved during the first year of operations.  The proton beam centroid at the target was stable to within $\unit[\pm 0.1]{mm}$ and the area of the beam-spot varied within the range $\unit[3.3-4.5]{mm^{2}}$.  Further details of the primary proton beam delivery system and its performance are given in Ref.~\cite{bishai2006}.

\subsubsection{Target and Horns}

\begin{figure}
  \centering
  \includegraphics[width=\columnwidth]{fig/horn-fig.eps}
  \caption{\label{fig:targ-horns} In the NuMI{} neutrino beam,  horns 1 and 2 are separated by \unit[10]{m}.  A collimating baffle upstream of the target protects the horns from direct exposure to misdirected proton beam pulses.  The target and baffle system can be moved further upstream of the horns to produce higher energy neutrino beams~\cite{kostin2002}.  The vertical scale is four times that of the horizontal (beam axis) scale.}
\end{figure}

\begin{figure*}
  \centering
\includegraphics[width=\columnwidth]{fig/three_beams_nearCC.eps}
  \caption{\label{fig:numi-spectra} Calculated rate of \ensuremath{\nu_{\mu}} charged-current{} interactions in the MINOS Near Detector, located \unit[1040]{m}{} from the NuMI{} target. The figures were made by combining the flux calculation discussed in \ref{sec:beam_sim} with the cross-sections discussed in Sec.~\ref{sec:mc}. Three spectra are shown, corresponding to the LE10/185kA, LE100/200kA, and LE250/200kA configurations of Tab.~\ref{tab:beam-config}.}
\end{figure*}

The production target is a rectangular graphite rod, segmented longitudinally into 47 segments, or fins. The target dimensions are \unit[6.4]{mm} in width, \unit[15]{mm} in height and \unit[940]{mm} in length. The beam spot size at the target is \unit[1.2-1.5]{mm}~\cite{bishai2006}.  A collimating baffle upstream of the target provides protection for the target, its cooling lines, as well as downstream beam components.  The baffle is a \unit[1.5]{m} long graphite rod with an \unit[11]{mm} diameter inner bore.

The particles produced in the target are focused by two magnetic horns~\cite{Abramov2002}, shown schematically in Fig.~\ref{fig:targ-horns}.  The \unit[200]{kA} pulsed current produces a maximum \unit[30]{kG} toroidal field which selects particles produced at the target by charge-sign and momentum.  Measurements of the horns show the expected $1/r$ fall-off in the magnetic field to within a percent and the field at distances from the beam axis smaller than the inner conductors is observed to be less than \unit[20]{G}.  The absolute value of the current flowing through the horns was calibrated to within $\pm 0.5\%$ and was observed to vary less than 0.2\% over the course of the data collection period.  The horns have parabolic-shaped inner conductors~\cite{Abramov2002}, for which the focal length for point-to-parallel focusing of particles from the target is proportional to momentum~\cite{budker1961,Baratov1977d}. The alignment of the target and horn system relative to the beam axis was checked using the proton beam itself~\cite{zwaska2006}.

The relative longitudinal positions of the two horns and the target optimizes the momentum focus for pions and therefore the typical neutrino energy.  To fine-tune the beam energy, the target is mounted on a rail-drive system with \unit[2.5]{m} of longitudinal travel, permitting remote change of the beam energy without directly accessing the horns and target~\cite{kostin2002}.  In its furthest downstream location, the target is cantilevered approximately \unit[65]{cm} into the first parabolic horn.  Characteristics of the predicted neutrino spectra from several target position and horn current settings are listed in Tab.~\ref{tab:beam-config} and \ensuremath{\nu_{\mu}} charged-current{} energy distributions are shown in Fig.~\ref{fig:numi-spectra}.  Data from the LE10/185kA configuration were used in the oscillation search. Data from the other beam configurations were used to constrain the beam Monte Carlo (see Sec.~\ref{sec:skzp}).

\begin{table*}
\begin{ruledtabular}
\begin{tabular}{lcccd}
Beam		&	Target          &       Horn  & Peak  & \mbox{Exposure}  \\
Configuration	&	Position (cm) 	&Current (kA) & $\ensuremath{E_{\nu}}{}\pm \mathrm{r.m.s.}$ (GeV)	& \unit[10^{18}]{POT}\\ \hline\hline
LE10/0kA	&	10		&	0	& $7.4 \pm 4.1$\footnote{The \unit[0]{kA} ``horn-off'' beam is unfocused and has a broad energy distribution} &  2.69 \\
LE10/170kA	&	10		&	170	& $3.1\pm 1.1$ &  1.34 \\
LE10/185kA	&	10		&	185	& $3.3\pm 1.1$ & 127. \\
LE10/200kA	&	10		&	200	& $3.5\pm 1.1$ & 1.26 \\
LE100/200kA	&	100		&	200	& $5.6\pm 1.5$ & 1.11 \\
LE250/200kA	&	250		&	200	& $8.6\pm 2.7$ & 1.55 \\

\end{tabular}
\end{ruledtabular}
    \caption{Beam configurations and data sets used in this publication. The target position refers to the distance the target was displaced upstream of its default position inside the first focusing horn.  The peak (i.e., most probable) neutrino energy \ensuremath{E_{\nu}}{} is determined after multiplying the muon-neutrino flux predicted by the beam Monte Carlo simulation by charged-current cross-section{}. The r.m.s.\ refers to the root mean square of the peak of the neutrino energy distribution. \label{tab:beam-config}}
\end{table*}

\subsubsection{Decay Volume and Absorber}

Particles are focused by the horns into a \unit[675]{m}{} long, \unit[2]{m} diameter steel pipe, evacuated to \unit[0.5]{Torr} to reduce meson absorption and scattering.  This length is approximately the decay length of a \unit[10]{GeV} pion. The entrance to the decay pipe is sealed by a two-piece aluminum-steel window. The central ($\mbox{radius} < \unit[50]{cm}$) portion of the window is made of \unit[1]{mm} thick aluminum and is strengthened by an outer ($\mbox{radius} > \unit[50]{cm}$) section made of \unit[1.8]{cm} thick steel. The design reduces scattering in the window while maintaining vacuum integrity. The decay volume is surrounded by \unit[2.5-3.5]{m} of concrete shielding.  At the end of the decay volume is a beam absorber consisting of a water-cooled aluminum and steel core followed by steel and concrete blocks.

\subsubsection{Instrumentation}

The primary proton beam position is monitored along the transport line by 24 capacitive beam position monitors (BPMs), the beam intensity by two toroidal beam current transformers and 44 loss monitors, and the beam position and spot size by ten retractable segmented foil secondary emission monitors (SEMs)~\cite{kopp2006a}.  During normal operations only the last SEM upstream of the target is inserted in the beam.  The absolute toroid uncertainty  was determined to be $\pm1.0\%$  by precision current pulses and monitored for drift throughout the run by comparison of the toroids in the NuMI{} line with monitors in the Main Injector.  The beam position and size at the target are measured to within \unit[$\pm50$]{$\mu$m}~\cite{bishai2006}.

In the target hall, the current flowing in the horns was monitored spill-to-spill, as was the temperature of the current-delivering strip-line.  Both were seen to change by approximately $0.2\%$ due to thermal variations in the target hall.  The temperature of the upstream collimating baffle was continuously monitored, as this was observed to be a good measure of the proton beam halo scraping on the baffle. During NuMI{} operations it was found that approximately 0.3\% of the beam was obstructed by the baffle.

\begin{figure*}[h]
  \centering
\includegraphics[width=0.8\columnwidth]{fig/mu_mon_1}
  \caption{\label{fig:mu_mon} Horizontal and vertical beam centroid measured by the first muon station. The observed variations in the centroid correspond to \unit[20]{$\mu$rad} variations in the beam angle. Such deviations are expected to have less than 0.2\% effect on the \ensuremath{\nu_{\mu}}{} energy spectrum at the Far Detector . }
\end{figure*}

Ionization chambers are used to monitor the secondary and tertiary particle beams~\cite{kopp2006b}.  An array is located immediately upstream of the beam absorber to monitor the remnant hadrons at the end of the decay pipe. There are also three muon monitoring stations, one downstream of the absorber, one after \unit[12]{m} of rock, and a third after an additional \unit[18]{m} of rock.  The muon stations confirmed the relative alignment of the neutrino beam direction to approximately \unit[20]{$\mu$rad} throughout the run as shown in Fig.~\ref{fig:mu_mon}. This is expected to limit spectral variations at the Far Detector to less than 0.2\%.  In addition, the charge (per spill per proton on target) deposited in the muon monitors varied by only $\pm 2\%$, indicating a similar level of stability in the neutrino flux.

The relative position of a point on the surface above the Far Detector and the NuMI{} target were measured by GPS survey with an accuracy of \unit[1]{cm} in the horizontal plane, and \unit[6]{cm} in the vertical~\cite{Bocean:2006mh,numitdr}. Translating the surface coordinates to the underground Far Detector hall introduces an additional uncertainty of \unit[70]{cm} in the vertical. The positions of beam-line components with respect to the target were surveyed and are known with an accuracy of \unit[0.5]{mm} (95\% C.L.).   These data indicate that the beam angle with respect to the ideal Near to Far Detector vector is less than \unit[15]{$\mu$rad} (68\%~C.L.). This introduces a negligible uncertainty into the calculation of the neutrino flux at the Far Detector.

\subsection{Simulation of the Neutrino Beam}
\label{sec:beam_sim}
\begin{figure}[h]
  \centering
  \subfloat[\label{fig:near_spec}]{\includegraphics[width=0.55\columnwidth]{fig/LE10_CCNear_log_v3.eps}}\\
  \subfloat[\label{fig:near_far_overlay}]{\includegraphics[width=0.55\columnwidth]{fig/far_near_cc_over_le010z185i.eps}}
  \caption{
  \protect\subref{fig:near_spec} Expected energy distribution of \ensuremath{\nu_{\mu}}{} charged-current interactions in the Near Detector. The contributions from $\pi$ and $K$ decays are shown separately. \protect\subref{fig:near_far_overlay} Comparison of the neutrino energy spectrum at the Near and the Far Detectors.  The two are not identical, due to solid angle differences between the two detectors (see Fig.~\ref{fig:near-far-diagram}).  \label{fig:near-far}}

\end{figure}

\begin{figure*}[h]
  \centering
  \includegraphics[width=1.0\textwidth]{fig/near-far-angles.eps}
  \caption{
The differences in the Near Detector and Far Detector spectra are due to the different angular acceptance of the two detectors and the proximity of the Near Detector to the beam-line \label{fig:near-far-diagram}}
\end{figure*}

The neutrino beam-line is modeled in three stages: (1) a simulation of the hadrons produced by \unit[120]{\ensuremath{\mathrm{GeV}/\mathit{c}}{}} protons incident on the NuMI{} target, (2) the propagation of those hadrons and their progeny through the magnetic focusing elements, along the \unit[675]{m}{} decay pipe and into the primary beam absorber allowing for decay of unstable particles and (3) the calculation of the probability that any produced neutrino traverses the Near and Far detectors.

The production of secondary mesons in the NuMI{} target was calculated using the {\tt FLUKA05}{}~\cite{fluka4} Monte Carlo.  Particles exiting the target are recorded and later propagated in a {\tt GEANT3}{}~\cite{g3man} simulation of the NuMI{} beam-line. The simulation describes the magnetic focusing horns, surrounding shielding and decay pipe. The magnetic field inside the horns was modeled as a perfect $1/r$ toroidal field.  The magnetic field in the inner conductor of the horn was calculated assuming that the current was uniformly distributed throughout the \unit[3]{mm} thickness of the aluminum conductor.  This assumption was motivated by the large skin depth $\delta\approx\unit[7.7]{mm}$ for the approximately $\unit[1]{ms}$ current pulse.  The \unit[10-30]{G} magnetic field in the inner apertures of the horns was neglected, as was the approximately $\unit[1]{G}$ field in the decay pipe due to residual magnetization of the iron vessel. The {\tt GEANT-FLUKA}{} code is used to describe hadronic interactions in the beam-line and the associated production of secondary particles as well as the full particle decay chains. Decays in which a neutrino is produced are saved and later used as input for neutrino event simulation in the Near and Far Detectors.

The neutrino event simulation uses each decay recorded by the NuMI{} simulation with a probability and neutrino energy determined by the decay kinematics and the (randomly chosen) trajectory through the Near or Far Detector. For two-body $\pi/K\rightarrow\mu\nu$ decays of relativistic mesons, the neutrino energy is given by

\begin{equation}
E_\nu\appro
\frac{(1-\frac{m_\mu^2}{M^2}) E}{1+\gamma^2\tan^2\theta_\nu} \quad .
\label{eq:enu-eqn}
\end{equation}

\noindent where $m_\mu$ and $M$ are the muon and parent hadron masses, $E$ the parent hadron energy, $\gamma=E/M$ is the parent's Lorentz boost, and $\theta_\nu$ is the angle (in the lab) between the neutrino and parent hadron directions.  The neutrino is forced to pass through either the Near or Far MINOS Detector, with probability for the particular meson decay given by

\begin{equation}
\frac{dP}{d\Omega_\nu}\approx\frac{1}{4\pi}\frac{4\gamma^2(1+\tan^2\theta_\nu)^{3/2}}{(1+\gamma^2\tan^2\theta_\nu)^2} \quad ,
\label{eq:angle}
\end{equation}

\noindent In deriving both expressions, $\beta=\sqrt{1-(1/\gamma^2)}\approx1$ is assumed. We have also accounted for the effect of $\mu$ polarization (e.g., induced in $\pi\rightarrow\mu\nu$) on the neutrino flux from $\mu$ decays.  Three-body kaon decays are also included, but contribute $<0.1\%$ to the \ensuremath{\nu_{\mu}}{} event rate.

\subsection{Expected Neutrino Energy Spectra}
\label{sec:expected_spectra}

Figure~\ref{fig:numi-spectra} shows the energy spectra of \ensuremath{\nu_{\mu}} charged-current{} interactions in the MINOS Near Detector in three of the beam configurations of Tab.~\ref{tab:beam-config}. In the LE10/185kA beam configuration, the neutrino in 87\% of \ensuremath{\nu_{\mu}} charged-current{} interactions was produced by $\pi^+\rightarrow\mu^+\nu_\mu$, with $K^+$ decays contributing the additional 13\%. $K^0$ and $\mu$ decays contribute $<0.1\%$ to \ensuremath{\nu_{\mu}}{} event rate, though they do contribute significantly to the \ensuremath{\nu_{e}}{} rate. The contributions of $\pi$ and $K$ to the neutrino energy spectrum is shown in Fig.~\ref{fig:near-far}.  The predicted spectrum for the LE10/185kA beam is composed of 92.9\% \ensuremath{\nu_{\mu}}{}, 5.8\% \ensuremath{\overline{\nu}_{\mu}}{} and 1.3\% \ensuremath{\nu_{e}}{}+\ensuremath{\overline{\nu}_{e}}{}.

The uncertainties on the neutrino flux due to focusing effects are expected to be largest at the edges of apertures of the horns.  The uncertainties are shown in Fig~\ref{fig:miserr_n}-\ref{fig:miserr_fon} and described in more detail in Sec~\ref{sec:flux_uncertainties}.  Uncertainties at higher neutrino energies are smaller because higher energy neutrinos arise from high energy $\pi$ and $K$ mesons which pass through the field-free apertures of the horns without focusing.

The neutrino energy spectra at the MINOS Near and Far Detectors are not identical, due to the different solid angles subtended by the two detectors and due to the energy-dependence of the meson decay locations.  As indicated in Fig.~\ref{fig:near-far-diagram}, higher-energy pions decay further down the length of the decay volume, and hence a variety of decay angles $\theta_\nu$ result in neutrinos which strike the Near Detector.  The Far Detector, by contrast, can only be struck for very small-angle decays.  By Eq.~\ref{eq:enu-eqn}, the different allowed decay angles, as well as the close proximity of fast pion decays to the Near Detector, results in different energy spectra at the two detectors. These differences are conveniently characterized in terms of the Far/Near spectral ratio shown in Fig.~\ref{fig:fn_before_after}.  Thus, while the Near Detector in principle measures the neutrino energy spectrum in the absence of oscillations, corrections must be applied to this spectrum to derive a predicted spectrum at the Far Detector, as is discussed in Sec.~\ref{sec:extrap}.

\subsection{Uncertainties In The Neutrino Flux}
\label{sec:flux_uncertainties}

\begin{figure*}[h]
  \centering
 \subfloat[\label{fig:miserr_n}]{\includegraphics[width=0.45\textwidth]{fig/err_near_le10.eps} }\\
 \subfloat[\label{fig:miserr_fon}]{\includegraphics[width=0.45\textwidth]{fig/err_fn_le10.eps} }
  \caption{
The effect that beam focusing uncertainties have on \protect\subref{fig:miserr_n} the Near Detector \ensuremath{\nu_{\mu}} charged-current{} spectrum and \protect\subref{fig:miserr_fon} on the ratio (Far/Near) of \ensuremath{\nu_{\mu}} charged-current{} spectra at the two detectors.  The lines show the fractional change in the number of events \protect\subref{fig:miserr_n} or event ratio \protect\subref{fig:miserr_fon} expected in each energy bin due to a one standard deviation shift in various beam parameters. These include the horn and target alignment, proton beam halo scraping on the collimator baffle, knowledge of the horn current, and the modeling of the current distribution in the horn conductor.  The total error is the sum in quadrature of the individual uncertainties. \label{fig:miserr}}
\end{figure*}

The dominant contribution to the uncertainty in the flux at both detector sites is caused by uncertainty in the yield of hadrons off the target as a function of \ensuremath{p_{z}}{} and \ensuremath{p_{T}}{}, the components of the hadron momentum along and transverse to the beam-line. The magnitude of the uncertainty on the energy spectra is difficult to estimate given the lack of hadron production data.  As described in Sec.~\ref{sec:skzp}, the data in the Near Detector can be used to constrain the beam flux calculation, and in particular the yield of hadrons off the target. The constrained calculation improves agreement with the Near Detector data and reduces uncertainties in the prediction of the Far Detector flux. The procedure used to apply this constraint is described in Sec.~\ref{sec:skzp}. Unless otherwise noted the Monte Carlo simulation will use the flux calculation constrained by the Near Detector data.

Systematic uncertainties on the predicted neutrino flux from beam focusing effects are readily calculable, and are constrained using data from the instrumentation in the primary, secondary, and tertiary beams. Figure~\ref{fig:miserr_n} shows the expected uncertainty on the neutrino energy spectrum at the Near Detector. Focusing uncertainties produce spectral distortions in both detectors that occur predominantly in the portions of the neutrino energy spectrum which correspond to parent hadrons which cross the edges of the horn apertures (see Fig.~\ref{fig:targ-horns}).  The largest beam focusing uncertainties include our knowledge of the absolute current flowing through the horns, the alignment of the horns and target to the rest of the beam-line~\cite{zwaska2006}, the shielding geometry, the uncertainty in the fraction of the proton beam scraping on the upstream collimating baffle (which acts as a target further upstream of the horns), and the modeling of the current distribution  in the horns' inner conductors (e.g., skin-depth effect). The expected distortion in the ratio of the Far and Near Detector spectra are smaller (see Fig.~\ref{fig:miserr_fon}) than the individual Near or Far Detector uncertainties (Fig.~\ref{fig:miserr_n}).  Hence, use of the Near Detector to measure the spectrum results in a prediction for the spectrum at the Far Detector with smaller uncertainty.

\subsection{Constraining the Beam Model}
\label{beammodeling}

\begin{figure*}
  \begin{centering}
    \includegraphics[width=0.35\textwidth]{fig/ptxf_le010z170i.eps}
    \includegraphics[width=0.35\textwidth]{fig/ptxf_le010.eps}
\\
    \includegraphics[width=0.35\textwidth]{fig/ptxf_le010z200i.eps}
    \includegraphics[width=0.35\textwidth]{fig/ptxf_le100.eps}
\\
    \includegraphics[width=0.35\textwidth]{fig/ptxf_le250.eps}
    \includegraphics[width=0.35\textwidth]{fig/ptxf_le010z000i.eps}
      \caption{The yield, \ensuremath{\ensuremath{\mathrm{d}}^{2}N/\ensuremath{\mathrm{d}} p_{z}\ensuremath{\mathrm{d}} p_{T}}{}, of $\pi^+$ off the target that contribute to the \ensuremath{\nu_{\mu}}{} charged-current event rate in each of the beam configurations of Tab.~\ref{tab:beam-config}. Here \ensuremath{p_{T}}{} and \ensuremath{p_{z}}{} refer to the transverse and longitudinal momentum of pions leaving the NuMI{} target before entering the focusing horns.  The box areas are proportional to the probability that a pion has the given (\ensuremath{\pz{},\pt{}}{}) and results in a \ensuremath{\nu_{\mu}}{} charged-current interaction in the Far Detector. As is evident, each beam configuration samples a different region of (\ensuremath{\pz{},\pt{}}{}).
      \label{fig:ptxf_allbeams}}
  \end{centering}
\end{figure*}

\begin{figure}
  \centering
  \includegraphics[width=\columnwidth]{fig/ptpz_weights_color.eps}
  \caption{Pion production scale factors, as a function of \ensuremath{\pz{},\pt{}}{}, derived by fitting Eq.~\ref{eq:skzp} to the Near Detector data according to the procedure described in the text. The scale factors are relative to the {\tt FLUKA05}{} prediction.}
  \label{fig:pzpt_weights}
\end{figure}

\begin{figure}
  \begin{centering}
    \includegraphics[width=\columnwidth]{fig/farovernear.eps}
    \caption{The predicted ratio (Far/Near) of \ensuremath{\nu_{\mu}}{} flux at the two detectors before and after the tuning procedure described in the text. The height of the box around each point denotes the uncertainty on the ratio.  For the points labeled ``Fluka 2005'' the uncertainties were calculated from the {\it a priori} uncertainties shown in Fig.~\ref{fig:miserr_fon} as well as the range of predictions from several hadron production codes~\cite{fluka4,mars1,malensek,bmpt}. For the points labeled ``Tuned MC'' the uncertainties are computed from the post-fit uncertainties on the functions $A',B'$ and $C'$ of Eq.~\ref{eq:skzp} as well as those on the additional parameters in Tab.~\ref{tab:fit_pars}. Unless specified otherwise, the tuned Monte Carlo is used in all further calculations. \label{fig:fn_before_after}}
  \end{centering}
\end{figure}

\begin{table}[t]
\begin{center}
\begin{ruledtabular}
\begin{tabular}{lcc}
parameter & initial uncertainty & best fit \\
\hline\hline
Horn 1 misalignment & \unit[1]{mm} & $\unit[0.1\pm 0.4]{mm}$  \\
POT Normalization & 2\% & $0.0\pm 0.8\% $\\
Horn Current Distribution $k$ & $\unit[0.17]{mm^{-1}}$ & $\unit[0.07\pm 0.05]{mm^{-1}}$\\
Shower Energy Offset & \unit[100]{MeV} & $\unit[24.7\pm 100.0]{MeV}$\\
Baffle scraping & 0.25\% & $ 0.1 \pm 0.1 \% $\\
Horn current calibration & 1\% & $0.9 \pm 0.3 \% $\\
Neutrino Energy Mis-calibration & 10\% & $-4.9 \pm 10.0\%$\\
Neutral-Current Contamination & 30\% & $-14\pm 7\%$\\
\end{tabular}
\end{ruledtabular}
\caption{Detector and beam modeling effects included when fitting the Monte Carlo to the Near Detector data. The effect of a $\unit[1]{\sigma}$ variation in each of the beam focusing parameters is shown in Fig.~\ref{fig:miserr_n}. The horn current distribution as a function of depth $z$ in the conductor is modeled as $I(z)\propto \exp(-k z)$ and we tabulate $k$ here.
\label{tab:fit_pars}}
\end{center}
\end{table}

 As noted in Sec~\ref{sec:flux_uncertainties} uncertainties in the neutrino flux calculation arise from insufficient knowledge of hadron production off the NuMI{} target as well as from several beam focusing effects described in Sec.~\ref{sec:flux_uncertainties}.  We constrain our flux calculation using \ensuremath{\nu_{\mu}}{} charged-current energy spectra measured in the Near Detector. The uncertainty in hadron production can be constrained because the position of the target and the horn current determine the region of pion (\ensuremath{\pz{},\pt{}}{}) which contributes to \ensuremath{\nu_{\mu}}{} at the Near Detector. This is shown in Fig.~\ref{fig:ptxf_allbeams} for each of the beam configurations in which we collected data (see Tab.~\ref{tab:beam-config}).  Our method works by representing the underlying production yield, \ensuremath{\ensuremath{\mathrm{d}}^{2}N/\ensuremath{\mathrm{d}} p_{z}\ensuremath{\mathrm{d}} p_{T}}{}, as a parametric function which we use to tune the Monte Carlo in a \ensuremath{\chi^{2}}{} fit to the Near Detector data.  We add terms to the \ensuremath{\chi^{2}}{} which describe the influence that beam focusing and detector modeling uncertainties have on the \ensuremath{\nu_{\mu}}{} spectrum (see Fig.~\ref{fig:miserr_n}).  This technique is similar to those used in previous experiments~\cite{ahrens1986,allen1993,astier2003,mcfarland1998,zeller2001} but the multiple beam configurations allow us to selectively enhance different regions of \ensuremath{p_{z}}{} and \ensuremath{p_{T}}{}.

The {\tt FLUKA05}{} prediction of $\pi^{+}$ yields off the NuMI{} target is well described by the function

\begin{equation}
f(\ensuremath{\pz{},\pt{}}{}) = \frac{\ensuremath{\mathrm{d}}^2N}{\ensuremath{\mathrm{d}} \ensuremath{p_{z}}{} \ensuremath{\mathrm{d}} p_T} = \left[A(\ensuremath{p_{z}}{})+ B(\ensuremath{p_{z}}{}) p_T \right]\exp\left(-C(\ensuremath{p_{z}}{})p_T^{3/2}\right)
\label{eq:skzp}
\end{equation}

\noindent where  $A$, $B$ and $C$ are functions of \ensuremath{p_{z}}{}. Equation~\ref{eq:skzp} and $A$, $B$ and $C$ are similar to the functions advocated by~\cite{bmpt} but we have modified them to better describe the thick target yields. The \ensuremath{\nu_{\mu}}{} spectra are then fit by warping $A,B,C$ as linear functions of \ensuremath{p_{z}}{}. The fit outputs the warped functions $A',B',C'$ and we calculate hadron yield scale factors relative to {\tt FLUKA05}{} as
\begin{equation}
 w(\ensuremath{\pz{},\pt{}}{})=\frac{f(\ensuremath{\pz{},\pt{}}{};A',B',C')}{f(\ensuremath{\pz{},\pt{}}{};A,B,C)}\quad .
\end{equation}
These scale factors are shown in Fig.~\ref{fig:pzpt_weights}.

In our fit we include a term to constrain the mean transverse momentum of pions, $\langle p_T\rangle$, to the {\tt FLUKA05}{} value of \unit[364]{\ensuremath{\mathrm{MeV}/\mathit{c}}{}}, with an uncertainty of \unit[15]{\ensuremath{\mathrm{MeV}/\mathit{c}}{}}, obtained from the variation of independent hadron production calculations~\cite{fluka4,mars1,malensek,wang,ckp}. The $K^{+}$ yield does not contribute significantly to the \ensuremath{\nu_{\mu}}{} flux below \unit[30]{GeV} and was permitted to vary by a scale factor. We augment the fit to include terms which represent distortions of the \ensuremath{\nu_{\mu}}{} charged-current spectra caused by uncertainties in the horn current, horn position, distribution of current in the horn conductors, and proton beam scraping on the baffle. The effect of these uncertainties are shown in Fig.~\ref{fig:miserr_n}. We also include a 2\% uncertainty in the normalization due to proton on target counting.  To avoid forcing the fit to over-correct for detector acceptance uncertainties we include a 10\% error on the neutrino energy, a \unit[100]{MeV} uncertainty in the shower energy and a 30\% uncertainty in the normalization of the neutral-current contamination.

The results of the fit are shown in Fig.~\ref{fig:6beamfit_all_enu}, Fig.~\ref{fig:pzpt_weights} and Tab.~\ref{tab:fit_pars}.  The agreement between the measured and calculated neutrino energy spectra in Fig.~\ref{fig:6beamfit_all_enu} is significantly improved with residual discrepancies reduced to 5\% or less in all beam configurations. The \ensuremath{\chi^{2}}{} per degree of freedom improves from $3087/360$ to $569/345$.  Figure~\ref{fig:pzpt_weights} shows the scale factors applied to the hadron production yield \ensuremath{\ensuremath{\mathrm{d}}^{2}N/\ensuremath{\mathrm{d}} p_{z}\ensuremath{\mathrm{d}} p_{T}}{}. The most significant adjustment is made to the region $\unit[\ensuremath{p_{T}}{}\lesssim 150]{\ensuremath{\mathrm{MeV}/\mathit{c}}{}}$ and $\unit[\ensuremath{p_{z}}{}\gtrsim 15]{\ensuremath{\mathrm{GeV}/\mathit{c}}{}}$ to increase the flux of unfocused hadrons which dominate the flux for neutrino energies larger than $\unit[8]{GeV}$ in the LE10 beam configurations. Only modest adjustments are made in the region focused in the LE10 configurations ($\unit[100\lesssim \ensuremath{p_{T}}{}\lesssim 500]{\ensuremath{\mathrm{MeV}/\mathit{c}}{}}$, $\unit[\ensuremath{p_{z}}{}\lesssim 10]{\ensuremath{\mathrm{GeV}/\mathit{c}}{}}$ ). Hadrons produced with high \ensuremath{p_{T}}{} and \ensuremath{p_{z}}{} do not contribute to neutrinos in the MINOS{} detectors and the scale factors are not constrained.  The best fit parameters, shown in Tab.~\ref{tab:fit_pars}, are all within $1\sigma$ of their initial values. Figure~\ref{fig:fn_before_after} shows the predicted ratio Far/Near of \ensuremath{\nu_{\mu}}{} charged-current energy spectra before and after the fitting procedure.  The fit improves our understanding of the relationship between the Far and Near spectra and we use the results in all further Monte Carlo calculations. Uncertainties in the ratio are shown as boxes which are reduced after the fit as a result of the constraints the neutrino data places on the functions $A,B,C$ of Eq.~\ref{eq:skzp} and the parameters of Tab.~\ref{tab:fit_pars}.  We have checked that alternate parametrizations of \ensuremath{\ensuremath{\mathrm{d}}^{2}N/\ensuremath{\mathrm{d}} p_{z}\ensuremath{\mathrm{d}} p_{T}}{} yield similar Far/Near ratios. This is important since the ratio characterizes the way in which the neutrino flux changes between the two detectors and indeed provides a way in which the Near Detector data can be used to predict the Far.

\subsection{Detector Calibration}
\label{sec:calib}

\begin{figure}
\begin{center}
\includegraphics[width=\columnwidth]{fig/blessed_drift.eps}
\caption{\label{fig:drift} Variations in the median signal per plane deposited by through-going cosmic-ray muons observed during the data-taking period covered by this paper. The time dependence is largely due to variations in the environmental conditions in the Near and Far Detector halls and aging of the scintillator.  The zero point on the ordinate is arbitrary.}
\end{center}
\end{figure}

\begin{figure}
\begin{center}
\includegraphics[width=\columnwidth]{fig/muon_calibration.eps}
\caption{\label{fig:muon_dedx}  Mean signal per scintillator strip induced by beam muons. Corrections are made with cosmic-ray muons and account for non-uniformities in channel gain $D(t)$, attenuation in optical fibers $A(i,x)$, and strip light output $U(i)$. }
\end{center}
\end{figure}

The principal tools for calibrating the detector response are an LED based light-injection (LI) system, a test-bench scan of the scintillator modules with a radioactive source and cosmic-ray muons.  The detector is calibrated in a multi-stage procedure that converts the raw photomultiplier signal $Q_{raw}\left(i,t,x\right)$ measured by channel $i$ at time $t$ for an energy deposition at position $x$ into a fully corrected signal $Q_{cor}$. Each calibration stage produces a ``calibration constant''.  The fully corrected $Q_{cor}$ is defined as the product of $Q_{raw}\left(i,t,x\right)$ and the calibration constant from each stage:

\begin{equation*}
Q_{cor} = Q_{raw}\left(i,t,x\right) \times D\left(t\right)  \times U\left(i\right) \times A\left(i,x\right)  \times E
\end{equation*}

\noindent where $D,U,A$ and $E$ refer to:

\begin{description}

\item[Drift correction $D\left(t\right)$:] The channel gains and their variation over time are measured with the LED based light-injection system, demonstrating that short-term ($<\unit[24]{h}$) gain variations are small and occurred mostly due to environmental changes in the detector halls. Light-injection data were eventually superseded by measurements of the mean signal per plane induced by through-going cosmic-ray muons since the muon data were also able to correct for variations in the scintillator light-output. The detector response varies by $<2\%$ over the data-taking period as shown in Fig.~\ref{fig:drift}. The decreasing response in the Fear Detector is likely due to aging of the scintillator.

\item[Uniformity correction $U\left(i\right)$:] Through-going cosmic-ray muons are used to account for differences in light output between individual strips as well as attenuation in the optical fibers and connections.  Event-by-event corrections are applied to account for the muon track angle and the expected inefficiency due to statistical fluctuations at low light-levels. The calibration reduces the uncorrected 30\% strip-to-strip response variations to approximately 8\%.
\item[Attenuation correction $A\left(i,x\right)$:] A radioactive source was used to map out the response of each scintillator module at many positions along each strip. This was done on a test-bench setup prior to installation of the modules in the Near and Far Detectors. The data were then fit to an empirical model of optical attenuation in the wavelength-shifting and clear optical fibers. The resulting parametrizations are used to correct the signals from cosmic-ray muons during the calculation of the uniformity calibration constants  $U\left(i\right)$ and also to correct the reconstructed shower energy for attenuation based on the reconstructed shower position. The relative size of the correction varies by about 30\% over the \unit[8]{m} length of a Far Detector scintillator strip when signals from both ends are added, and by about 50\% over the \unit[3]{m} length of a Near Detector strip.

\item[Signal scale calibration $E$:] The overall scale of the signals is established by the detector's response to stopping muons. This provides the standard which fixes the absolute calibration of the Near and Far Detectors, allowing their signals to be compared to one another.  The calibration is done by tabulating the detector's response to muon crossings using only the portion of each track in which the muon energy is between \unit[0.5]{GeV} and \unit[1.1]{GeV}, deduced from the distance to the track endpoint. This energy window avoids the rapid variation in \ensuremath{\D E /\D x}{} near the track's end. Prior to this calibration the signals have already been corrected for $D,A,U$ as described above. For each muon a correction is applied to account for the muon's path length in each scintillator plane. The mean response is then calculated for each individual strip and a single constant representing the entire detector is derived from the median over all strips.

\end{description}

The effect of the uniformity and attenuation calibration stages is shown in Fig.~\ref{fig:muon_dedx}.  After calibration, the fully corrected signal is expressed in muon equivalent units (MEU) as defined by the signal scale ($E$) calibration procedure. One MEU corresponds to approximately 3.7 photoelectrons at the center of the Far Detector (7.5 photoelectrons when signals from both strip-ends are summed) and 5.4 photoelectrons at the beam center in the Near Detector. Based on muon stopping-power tables~\cite{Groom:2001kq} we find that one MEU corresponds to $\unit[2.00 \pm 0.02]{MeV}$ of muon energy loss in scintillator~\cite{Adamson:2006xv,jeff_thesis}. This calibration was independently derived at the Near and Far Detectors with a relative uncertainty of $\pm 2\%$.

\label{sec:shwr_calib}Data from a dedicated, test-beam calibration detector~\cite{Adamson:2006xv} are used to benchmark the hadronic and electromagnetic shower simulation done by {\tt GEANT}{}/{\tt GCALOR}{}~\cite{mike_thesis,tricia_thesis}. The measurements demonstrate that the simulated response agrees with the data at the level of 1--5\%, depending on the energy and particle type. This level of agreement is inclusive of energy dependent offsets in the beam momentum of the test beam-lines and uncertainties in the simulation of energy loss upstream of the calibration detector. In light of the much larger uncertainty associated with final state interactions, no attempt was made to resolve these discrepancies when determining the conversion from MEU to GeV for neutrino induced showers. Instead a conservative uncertainty of 5.6\% is assigned to the absolute shower energy scale to account for the data/MC discrepancies from the calibration detector as well as the precision of the stopping muon calibration performed there. An energy dependent MEU to GeV conversion is then extracted from the simulation such that the reconstructed shower energy estimates the energy transferred to the hadronic system.

\section{Please Note}

\section{Introduction}
\input{intro}

\section{The Neutrino Beam}
\input{beam_description}
\input{beam_mc}

\section{The Neutrino Detectors}

\input{det.tex}

\subsection{Event Reconstruction}
\input{reco.tex}

\input{mc.tex}
\input{calib.tex}

\section{CC Event Selection}
\input{evsel.tex}

\section{Constraints on the Neutrino Flux From Near Detector Data}
\input{beam_tuning}

\section{Prediction of the Far Detector Spectrum}

\label{sec:extrap}
\input{extrap_methods.tex}
\input{extrap_syst.tex}

\section{Oscillation Analysis}
\input{oscillationfit.tex}

\section{Summary}

\input{summary.tex}

\begin{acknowledgments}
\input{ack.tex}
\end{acknowledgments}

\ifthenelse{\isundefined{\beginnumbering}}{}{\endnumbering}

\subsection{Near Detector Distributions}

The basic reconstructed quantities, including vertex distributions and track angular distributions are observed to be in good agreement between data and Monte Carlo (as shown in fig~\ref{fig:nd_lowlevel} for events that satisfy the pre-selection cuts).
The distributions of the three PDF distributions that are inputs to the CC/NC PID parameter are shown in figure~\ref{fig:nd_pid} for data and Monte Carlo. In general the shapes of these distributions are well-modelled by the MC. The resulting PID distributions are also shown in figure~\ref{fig:nd_pid} - the shapes of these are also in good agreement - especially in the region $PID>-0.1$, which is used for the flux extrapolation to the Far Detector. The residual shape difference between data and MC below $PID=-0.1$ is used to estimate the uncertainty on the NC contamination in the CC-like sample, as will be described later.

\begin{figure}[h]
\centering
\includegraphics[width=\columnwidth]{./fig/NDComp_4plots_bw.eps}
\caption{Distributions of the three input Probability Density Functions to the CC/NC Particle Identification (PID) parameter for Near Detector data and Monte Carlo. The bottom-right plot shows the PID distribution that derives from these three PDFs.
  } \label{fig:nd_pid}
\end{figure}

\subsection{Far Detector Distributions}

We also performed the same checks on the selected sample of Far Detector Data events. A visual scan of the events was carried out in order to ensure that there was no evidence of background contamination in the selected event sample and that the events were well-reconstructed. Figure~\ref{fig:fd_vtx} shows that the reconstructed vertex distributions of the Far Detector events that pass the pre-selection cuts are consistent with the Monte Carlo expectation (the Monte Carlo distributions are normalised to the same number of events as the data to account for the possible effects of oscillations).

\begin{figure}[h]
\centering

\includegraphics[width=\columnwidth]{./fig/prd_fd_vertices.eps}

\caption{Distributions of track vertices for Far Detector data and Monte Carlo after pre-selection cuts. The Monte Carlo distributions are normalised to the same number of events as the data} \label{fig:fd_vtx}
\end{figure}

\vspace{0.1in}
\par

Figure~\ref{fig:fd_pid} shows the distributions of the three input variables that are used in the charged-current event selection, and the resulting PID parameter itself. Since these variables are sensitive to neutrino oscillations, the Monte Carlo prediction is shown with and without oscillations (where the mixing parameters are chosen to be the best fit values obtained from the fit to the Far Detector reconstructed energy spectrum, described below). The data is well-described by the oscillated Monte Carlo prediction in all of these distributions.

\begin{figure}[h]
\centering

\includegraphics[width=\columnwidth]{./fig/prd_fd_pid.eps}

\caption{Distributions of the three input Probability Density Functions to the CC/NC Particle Identification (PID) parameter for Far Detector data and Monte Carlo. The bottom-right plot shows the PID distribution that derives from these three PDFs. The Monte Carlo distributions are shown both without oscillations (solid histograms) and weighted by the best-fit neutrino mixing parameters from the oscillation analysis described below (dashed histograms).
  } \label{fig:fd_pid}
\end{figure}

In conclusion, we have observed a satisfactory level of agreement between basic reconstructed distributions for data and Monte Carlo in both Near and Far detectors - which indicate that no significant reconstruction or intensity-related biases exist in the data. We have observed differences in the shape and overall rate of the reconstructed energy spectra in several beam configurations, which is consistent with known systematic uncertainties in beam flux and cross-section predictions. Since we do not beleve that these differences are due to reconstruction biases, the measured spectum can be used ``as-is'' to predict the Far Detector flux. However, we can use the information contained in the measured Near Detector spectra for the various beam configurations to provide constraints on our beam model, and therefore improve the level of agreement between data and Monte Carlo - the procedure for which will now be described.

\subsection{Detector Technology}

The MINOS{} detectors are magnetized steel/scintillator tracking/sampling calorimeters designed to measure muon-neutrinos produced by the NuMI{} beam.  The active medium comprises \unit[4.1]{cm}-wide, \unit[1.0]{cm}-thick plastic scintillator strips arranged side and encased within aluminum sheets to form light-tight modules of 20 or 28 strips. Modules are combined to form scintillator planes which are then mounted on a \unit[2.54]{cm}-thick steel absorber plate. The detectors are composed of a series of these steel/scintillator planes hung vertically at a \unit[5.94]{cm} pitch. The scintillator strips in successive planes are rotated by $90^{\circ}$ to measure the three dimensional event topology.   The strips were produced in an extrusion process during which a \unit[2.0]{mm}-wide, \unit[2.0]{mm}-deep groove was driven along each strip. The entire strip, apart from the groove, was co-extruded with a \unit[0.25]{mm}-thick $\mathrm{TiO_{2}}$ doped polystyrene reflective layer.  Scintillation light is collected by a \unit[1.2]{mm}-diameter wavelength-shifting fiber that was glued into the groove with optical epoxy. The fiber transports the light to an optical connector located on the edge of the scintillator layer where it is transferred to a clear polystyrene fiber and routed to multi-anode photomultipliers, as shown in Fig.~\ref{fig:optical_readout}. Light emitting diodes are used to produce controlled pulses of light during the data-taking to track gain changes in and characterize the response of the photomultipliers and electronics~\cite{Adamson:2002ze}. The data acquisition (DAQ)~\cite{Belias:2004bj} and timing system synchronize and continuously read out the front-end electronics. Software triggering in the DAQ provides flexible event selection and data processing. GPS timestamps allow data from the two detectors to be synchronized with the beam pulses. The two detectors have different front-end electronics due to the disparate rates of neutrino interactions and cosmic-ray crossings at the two sites.

\subsection{The Far Detector}

The Far Detector has a mass of \unit[5400]{metric tons} and consists of 486 steel plates arranged in two modules separated by a \unit[1.1]{m} gap. The plates are \unit[2.54]{cm} thick, \unit[8]{m} wide octagons. A \unit[15.2]{kA} magnet coil is routed through a hole in the center of each plane and induces an average field of \unc{\unit[1.27]{T}} in the steel. Each scintillator plane covers the \unit[8]{m} octagon and consists of 192 scintillator strips read out from both ends by Hamamatsu 16-anode photomultipliers. The signals from eight strip-ends, separated by approximately $\unit[1]{m}$ within a single plane, are optically summed onto a single photomultiplier channel. The optical summing pattern is different for the two sides of the detector, allowing for resolution of the eight-fold ambiguity associated with each hit.

The front-end electronics~\cite{Oliver:2004ek} digitize signals from the photomultipliers with a 14 bit ADC (\unit[2]{fC} precision) when the common dynode signal exceeds approximately $\unit[0.25]{photoelectrons}$ and time-stamps them with a \unit[1.5625]{ns} least significant bit.  To reduce the electronics dead-time, hits are only digitized if more than one of the photomultipliers serving a contiguous group of 20 or 24 planes on one side of the detector are above threshold in coincidence. Data selected for further processing are pedestal suppressed and sent to a trigger farm where software triggering and data processing are performed and GPS timestamps are applied.

The DAQ has data buffering that is large enough to allow it to wait for the GPS time-stamp of the spill to arrive from the Near Detector over the Internet. This is used to form a bias-free beam trigger by recording all hits in the detector within a $\unit[100]{\mu s}$ window around the time of the spill. To avoid splitting events, the window size is extended for each spill to ensure that the entire window is bounded on both sides by an activity free period of at least \unit[156]{ns}. Finally, all hits in the $\unit[30]{\mu s}$ prior to each trigger are added to the event to provide a mapping of channels that were engaged in digitization at the time of the trigger. Fake spill triggers are also generated to monitor backgrounds. When spill information is not available at the Far Detector site, and for all out-of-spill triggering, candidate events are formed from time sequential blocks of hits bounded on each side by at least \unit[156]{ns} in which no detector activity occurred; trigger algorithms based on the spatial and energy clustering of hits in the detector are then applied to select events of potential interest. The integrated trigger rate in the Far Detector is typically \unit[4]{Hz} and is dominated by cosmic-rays and single photoelectron noise.

\subsection{The Near Detector}

The \unit[980]{metric ton} Near Detector consists of 282 steel/scintillator planes arranged in a single magnetized module.  Each steel plate has a ``squashed-octagon'' shape \unit[6.2]{m} wide and \unit[3.8]{m} high with a $\unit[30\times 30]{cm^{2}}$ hole offset \unit[56]{cm} from the horizontal center of the detector to accommodate a magnet coil.  The \unit[40]{kA} carried by the coil induces a \unc{\unit[1.17]{T}} field at the neutrino beam center, located \unit[1.49]{m} to the left of the coil. The coil and detector geometry was designed to provide a magnetic field in the region around the beam center that is similar to the field in the Far Detector. The Near Detector has two different types of scintillator planes: partially instrumented and fully instrumented. The planes are smaller than those in the Far Detector and are read out on only one side using Hamamatsu 64 anode photomultipliers. Fully instrumented Near Detector planes have 96 strips and cover a $\unit[13.2]{m^{2}}$ area. These planes are attached to one in every five steel plates along the entire length of the detector. Partially instrumented planes are attached to four out of five plates in the upstream portion of the detector. These planes are comprised of 64 strips covering a $\unit[6.0]{m^{2}}$ area.

The Near Detector is organized in two sections. The upstream 121 steel plates form the calorimeter, which is used to define the interaction vertex, find the upstream portion of muon tracks, and measure the energy of the neutrino induced hadronic shower. In the calorimeter each plate is instrumented with scintillator and the signals from each strip are read out independently. The downstream 161 plates are used as a muon spectrometer and only one in every five is instrumented. Furthermore, in the spectrometer, the signals from four strips are summed onto one electronics channel. The four-fold ambiguity is resolved in the event reconstruction program by extrapolating the muon track found in the upstream portion of the detector.

At the typical beam intensity of $\unit[2.2\times 10^{13}]{POT/spill}$, an average of 16 (44) neutrino interactions~\footnote{This includes interactions occurring throughout the detector, many of which cannot be fully reconstructed but which nevertheless induce some activity in the detector.}  occurred in the Near Detector during each $\unit[10]{\mu{}s}$ spill in the LE10/185kA (LE250/200kA) beam configuration. About half of these events occur in the calorimeter region and may be fully reconstructed. The Near Detector readout electronics~\cite{Cundiff:2006yz} is designed to measure neutrino interactions throughout the spill, without dead-time, and with a timing resolution that allows efficient separation of time-adjacent events.  Analog to digital conversion has a floating precision with a minimum least significant bit of \unit[1.4]{fC} and occurs in contiguous \unit[18.8]{ns} time-intervals corresponding to the \unit[53]{MHz} RF of the Main Injector. The electronics has two primary operating modes, ``spill-gate'' and ``dynode'', which are switched in real time.  In spill-gate mode the output from every photomultiplier pixel is digitized continuously in a $\unit[13]{\mu{}s}$ period starting about $\unit[1.5]{\mu{}s}$ before the arrival of neutrinos at the Near Detector. In dynode mode, used for out-of-spill acquisition of cosmic-rays, continuous digitization for a period of \unit[150]{ns} is initiated independently for each photomultiplier when the dynode signal exceeds a programmable threshold. The integrated trigger rate is typically \unit[30]{Hz}.

\subsection{Event Classification}
\label{sec:evtsel}

In the MINOS experiment, neutrino oscillations are expected to cause a deficit of \ensuremath{\nu_{\mu}} charged-current{} events at energies $E_{\nu}\lesssim\unit[5]{GeV}$ for $\ensuremath{\Delta m^2} \approx \unit[2.5\times 10^{-3}]{\ensuremath{\mathrm{eV}^{2}/\mathit{c}^{4}}{}}$ and our baseline of \unit[735]{km}. Neutral-current events which have been misidentified as \ensuremath{\nu_{\mu}} charged-current{} tend to populate this energy range and could obscure the oscillation signal. We remove neutral-current events from the oscillation sample with a technique based on the event topology. The technique uses three probability density functions (PDFs), which are constructed for the following variables: (a) the event length, expressed in units of the number of detector planes, (b) the fraction of the total event signal in the reconstructed track, and (c) the average signal per plane induced by the reconstructed track. These quantities are related to the muon-range, the event inelasticity, and the average energy loss \ensuremath{\D E /\D x}{} of the muon track and are distributed differently for \ensuremath{\nu_{\mu}} charged-current{} and neutral-current events as shown in Fig.~\ref{fig:pdf}. The probability that a particular event is consistent with the \ensuremath{\nu_{\mu}} charged-current{} or neutral-current PDFs is given by the product of the three individual probabilities

\begin{equation} \label{eq:P}
P_{CC,NC}=\prod_{i=1}^{3}f_{i}(x_{i})_{CC,NC} \quad,
\end{equation}

\noindent where the $f_{i}(x_{i})$ are the individual PDFs for \ensuremath{\nu_{\mu}} charged-current{} and neutral-current events respectively. The PDF distributions are well modeled by the simulation in the regions in which the neutral-current and \ensuremath{\nu_{\mu}} charged-current{} samples overlap and for neutrino energies above approximately \unit[10]{GeV} where the Monte Carlo simulation is better constrained by data from previous experiments and the neutrino flux has a weaker dependence on the energy.

\begin{figure}
\centering
  \subfloat[\label{fig:pdf1}]{\includegraphics[height=0.25\textheight,clip]{fig/pdf1.eps}}\\
  \subfloat[\label{fig:pdf2}]{\includegraphics[height=0.25\textheight,clip]{fig/pdf2.eps}}\\
  \subfloat[\label{fig:pdf3}]{\includegraphics[height=0.25\textheight,clip]{fig/pdf3.eps}}
  \caption{Distributions of the variables used to the define the event separation parameter $S$ in Eq.~\ref{eq:P} and Eq.~\ref{eq:S} for events satisfying the criteria of Sec.~\ref{sec:reco}. Near Detector data collected in the LE10/185kA beam configuration is shown.  \label{fig:pdf} }
\end{figure}

An event selection parameter $S$ is derived from $P_{CC}$ and $P_{NC}$ according to

\begin{equation} \label{eq:S}
S=-\left(\sqrt{-\ln(P_{CC})}-\sqrt{-\ln(P_{NC})}\right) \quad.
\end{equation}

\noindent Events that are more likely to originate from the \ensuremath{\nu_{\mu}} charged-current{} interactions give positive values, and those that are more likely to be neutral-current give negative values. Figure~\ref{fig:pid} shows the distribution of $S$ for data collected in the Near Detector in the LE10/185kA, LE100/200kA and LE250/200kA beam configurations. The behavior of the data is well-reproduced by the Monte Carlo simulation especially in the region $S>-0.1$ in which \ensuremath{\nu_{\mu}} charged-current{} events dominate.

\begin{figure}
\centering
  \subfloat[\label{fig:pid_le}]{\includegraphics[height=0.25\textheight]{fig/nd_pid_le10.eps}}\\
  \subfloat[\label{fig:pid_me}]{\includegraphics[height=0.25\textheight]{fig/nd_pid_le100.eps}}\\
  \subfloat[\label{fig:pid_he}]{\includegraphics[height=0.25\textheight]{fig/nd_pid_le250.eps}}
  \caption{ The event separation parameter $S$, as defined in Eq.~\ref{eq:S}, plotted for Near Detector data and simulation, and for three beam configurations in Tab.~\ref{tab:beam-config}.  \label{fig:pid} }
\end{figure}

\begin{figure}
\centering
  \subfloat[\label{fig:eff_nd}]{\includegraphics[width=0.7\textwidth]{fig/pideffic_nd.eps}}
\\
  \subfloat[\label{fig:eff_fd}]{\includegraphics[width=0.7\textwidth]{fig/pideffic_fd.eps}}
  \caption{Efficiency for identifying \ensuremath{\nu_{\mu}} charged-current{} events using the event separation parameter defined in Eq.~\ref{eq:S}, along with the purity of the resulting sample. The efficiency is defined with respect to the events satisfying the criteria of Sec.~\ref{sec:reco}. \label{fig:pideffic} \label{fig:pid_eff} }
\end{figure}

 We isolate an enriched sample of \ensuremath{\nu_{\mu}} charged-current{} events for use in the oscillation analysis by requiring $S>-0.2$  in the Far Detector and $S>-0.1$ in the Near Detector. These values were chosen to optimize the statistical sensitivity to the oscillation parameters $\ensuremath{\Delta m^2}$ and $\ensuremath{\sin^{2}2{\theta}}$ in the Far Detector, and to provide a selected event sample in the Near Detector that has the same purity, defined as the fraction of selected events that are true \ensuremath{\nu_{\mu}} charged-current{} interactions, as the Far Detector.  Figure~\ref{fig:pid_eff} shows the selection efficiency, defined as the fraction of true \ensuremath{\nu_{\mu}} charged-current{} events that have $S>-0.1$ in the Near Detector ($S>-0.2$ in the Far Detector), and the resulting purity of the sample in the two detectors. The efficiency is calculated with respect to the events passing the criteria of Sec.~\ref{sec:reco} and varies slowly as a function of reconstructed energy above \unit[0.5]{GeV}. The contamination from misidentified neutral-current interactions is greater at low reconstructed energy, hence the purity of the selected sample drops below \unit[1]{GeV}. The Monte Carlo simulation indicates that the average selection efficiency and purity (integrated over the unoscillated Far Detector LE010/185kA spectrum) are $90.0\%$ and $98.2\%$ respectively.

\begin{table*}[]
\begin{center}

\begin{tabular}{lc}
\hline
\hline
Selection Criterion & Far Detector Events \\
\hline

Reco. track in fiducial vol. & 427\\

Data quality cuts & 408\\

Event Timing cut & 404\\

Beam Quality cuts & 390\\

Track direction cut & 384 \\

Track quality cuts & 365 \\

Negative track charge & 306 \\

Reco. energy $<$ 30 GeV & 275\\

$S>-0.2$ & 215\\
\hline
\hline

\end{tabular}
\caption{Effect of selection criteria on the Far Detector neutrino event sample. The first eight criteria are described in Sec.~\ref{sec:reco}. In accordance with our analysis strategy, the total number of events passing all cuts was not revealed until all of the data quality checks were complete and our analysis method was fully defined.\label{tab:fdsel}}
\end{center}
\end{table*}

\begin{figure}[h]
\centering
\subfloat[\label{fig:fd_pid1}]{\includegraphics[width=0.45\textwidth]{./fig/prd_fd_pid_1.eps}}
\hfill
\subfloat[\label{fig:fd_pid2}]{\includegraphics[width=0.45\textwidth]{./fig/prd_fd_pid_2.eps}} \\

\subfloat[\label{fig:fd_pid3}]{\includegraphics[width=0.45\textwidth]{./fig/prd_fd_pid_3.eps}}
\hfill
\subfloat[\label{fig:fd_pid4}]{ \includegraphics[width=0.45\textwidth]{./fig/prd_fd_pid_4.eps}}

\caption{Distributions of the variables that are used in \ensuremath{\nu_{\mu}} charged-current{} and neutral-current event classification, for Far Detector data and Monte Carlo simulation. The bottom-right plot shows the distribution of the event separation parameter $S$ that is derived from these three PDFs. Oscillations are accounted for according to the mixing parameters extracted from our data (see Sec.~\ref{sec:osc}). These distributions are shown here to illustrate the level of agreement between data and simulation for $S$ and its input variables. They were not examined until the analysis procedure had been fully defined.\label{fig:fd_pid} }
\end{figure}

The PDFs used to discriminate between \ensuremath{\nu_{\mu}} charged-current{} and neutral-current interactions depend on the \ensuremath{\nu_{\mu}}{} energy spectrum and are therefore sensitive to neutrino oscillations. As such, the discrimination technique was developed and validated by comparing the Monte Carlo simulation to data collected in the Near Detector. A fraction of the  Far Detector dataset was available for data quality checks. Examination of oscillation-sensitive distributions in the full Far Detector data set was only performed after these checks had been made and we had fully defined our analysis procedure. Table~\ref{tab:fdsel} shows the effect of our selection on those Far Detector events that are recorded in coincidence with NuMI{} beam spills. After requiring $S>-0.2$ we retain 215 of the 275 events passing the criteria of Sec.~\ref{sec:reco}. The variables used in Eq.~\ref{eq:P} and Eq.~\ref{eq:S} are shown in Fig.~\ref{fig:fd_pid} for the 275 event sample. Two versions of the simulation are shown. The first does not account for neutrino oscillations while in the second we apply the effect of oscillations according to the mixing parameters extracted from our data using the procedure described in Sec.~\ref{sec:osc}. The measured deficit of events with respect to the unoscillated simulation is consistent with the partial disappearance of \ensuremath{\nu_{\mu}} charged-current{} events at energies less than about \unit[10]{GeV} and is well described when we account for oscillations.

\subsection{Backgrounds}

\subsubsection{Neutral-Current Contamination}
\label{sec:nc}
The neutral-current background in the Far Detector is estimated from the Monte Carlo simulation as the number of neutral-current events that pass the event selection cut. These events comprise only $<2\%$ of the selected unoscillated \ensuremath{\nu_{\mu}}{} CC sample but tend to congregate in the lowest reconstructed energy bins, as shown in Fig.~\ref{fig:pideffic}, and so constitute a significant background to the \ensuremath{\nu_{\mu}} charged-current{} oscillation signal. In this section, we describe how we have used Near Detector data to directly provide an estimate of the systematic error on the number of neutral-current events that are selected in the final Far Detector event sample.

To handle uncertainties that arise from the modeling of hadronic showers in the Monte Carlo, showers from cleanly identified \ensuremath{\nu_{\mu}} charged-current{}  events in the Near Detector with the muon tracks removed were used to model neutral current events.  The ratio between the track-removed shower $S$ distributions for data/MC provides a $S$-dependent Monte Carlo scaling factor for the shape of the neutral-current event selection parameter distribution.

 To determine the uncertainty on the overall neutral-current signal normalization, the event selection parameter distribution $S$ in the data was fitted to determine the amplitude of the \ensuremath{\nu_{\mu}} charged-current{} and neutral-current signals.  In this fit, the shape of the charged-current component was obtained from the Monte Carlo simulation, while the neutral-current shape was scaled according to the data/MC ratio observed in the track-removed \ensuremath{\nu_{\mu}} charged-current{} event sample.  The fits to the data for the signal amplitude and the data/MC comparison of the track-removed \ensuremath{\nu_{\mu}} charged-current{} events were performed in six reconstructed energy slices.

A comparison of the event separation distribution $S$ for the simulated neutral-current and selected muon removed \ensuremath{\nu_{\mu}} charged-current{} events yield reasonable agreement, as shown in Fig.~\ref{fig:mrpid}. However, the muon removed \ensuremath{\nu_{\mu}} charged-current{} events from Near Detector data, which are represented by the black points in Fig.~\ref{fig:mrpid}, peak at lower values of $S$ than the corresponding Monte Carlo distribution. This discrepancy implies that the data contain shorter events and have higher average pulse height per plane than the Monte Carlo events.

\begin{figure}
\centering
\includegraphics[width=\columnwidth]{fig/mrpid.eps}
  \caption{Distribution of the event separation parameter $S$ for true neutral-current Monte Carlo events (shaded histogram), for selected $\nu_\mu$  Monte Carlo events with the reconstructed muon track removed (open circles) and for selected Near Detector data events with the muon track removed (filled circles). \label{fig:mrpid} }
\end{figure}

The ratio of the data to simulated $S$ distributions obtained from the muon-removed events is used to provide a reweighting of the shape of the $S$ distribution for neutral-current Monte Carlo events. This procedure is performed in six bins of reconstructed energy and in all cases, the reweighted distributions predict fewer neutral-current events for $S>-0.1$ than the {\it a priori} simulation, by an amount that ranges from 20-50\%, as shown in Table~\ref{table:nc_summary}.

\begin{table}[]

\centering
\begin{tabular}{lc}
\hline
\hline
Energy Range\hspace{0.3in} &  Data/MC ratio \\
\hline

0-1 GeV & 0.642 \\
1-2 GeV & 0.503 \\
2-4 GeV & 0.804 \\
4-6 GeV & 0.583 \\
6-12 GeV & 0.434 \\
12-30 GeV & 0.579 \\
\hline
\hline
\end{tabular}

\caption{The ratio of the number of neutral-current events with $S>-0.1$ in the simulation and the best fit to the data using the neutral-current shape derived from the track-removed \ensuremath{\nu_{\mu}} charged-current{} events. The data indicates that the neutral-current background is over-predicted by the simulation. \label{table:nc_summary} }
\end{table}

The origin of this difference is most likely due to an unknown combination of shower modeling and/or neutral-current cross-section uncertainties. We have thus assigned 50\%, the largest difference observed in the fit above $S=-0.1$, as the error on the neutral-current background.

\subsubsection{Cosmic Ray Background}

We have estimated the background rate from cosmic-ray muons in the Far Detector selected sample using two independent methods. First, we examined the rate of selected events in the spill trigger window that are not within the expected $\unit[10]{\mu s}$ wide window around the time of true beam spill (see Fig~\ref{fig:fdtiming}). Secondly, we applied the \ensuremath{\nu_{\mu}} charged-current{} selection criteria described above to a sample of $2.3\times 10^{6}$ spill triggers taken in anti-coincidence with the beam spill. Both of these methods yielded an upper limit (90\% C.L) of 0.5 background events in a $\unit[10]{\mu s}$ wide window around the time of true beam spill for an exposure of $\unit[1.27\times10^{20}]{POT}$. The background due to neutrino-induced \ensuremath{\nu_{\mu}} charged-current{} interactions in the rock upstream of the Far Detector is estimated from the Monte Carlo simulation to be \unit[0.38]{events}. The rate of cosmic-ray muons is about 30 times larger in the Near Detector than it is in the Far Detector. The cosmic-ray background is negligible, however, due to the much greater ($\times 10^{6}$) neutrino flux at the Near Detector.

\subsection{Extrapolation Techniques}
\label{sec:extrap_methods}

We refer to the process of predicting the Far Detector spectrum as ``extrapolation'' because of the crucial role played by the Near Detector data, and the expected similarity~\footnote{Modulo the ratio of solid angles $\Omega_{FD}/\Omega_{ND} = d^{2}_{ND}/d^{2}_{FD}\approx(\unit[1]{km}/\unit[735]{km})^{2}$ where $d_{ND},d_{FD}$ refer to the distance from the target to the detectors.} between the neutrino fluxes at the two sites (see Fig.~\ref{fig:near-far}, Fig.~\ref{fig:fn_before_after} and the discussion in Sec.~\ref{sec:expected_spectra}). The extrapolation process may then, in one approach, be viewed as making $\sim 30\%$\footnote{These adjustments are relative to the solid angle correction $\Omega_{FD}/\Omega_{ND}$.} adjustments based on the energy spectrum measured in the Near Detector in order to predict the Far Detector spectrum. These adjustments are made, in part, by using the Monte Carlo simulation but are relative between the two detectors and are thus less sensitive to uncertainties in the absolute flux, cross-sections and detector acceptance. We refer to this approach as ``direct'' extrapolation because the Near Detector data are used as measured without further constraining the simulation. In a second approach Near to Far Detector extrapolation may be viewed as using the Near Detector data to constrain the Monte Carlo calculation of the neutrino flux, cross-sections and detector acceptance. The improved simulation is then used to calculate the energy spectrum expected at the Far Detector. This approach is referred to as ``indirect'' extrapolation in contrast to the ``direct'' approach. In practice the two approaches are complementary and result in very similar predictions of the Far Detector spectrum. The present section describes the specific techniques that were developed to implement each approach. We have pursued a total of four techniques in order to better understand the extrapolation process and study the robustness of our oscillation measurement.

\subsubsection{Direct Extrapolation Methods}

In the absence of oscillations the neutrino fluxes at the Far and Near Detectors are similar, but not identical, and the differences can be most transparently characterized in terms of the Far to Near flux ratio shown in Fig.~\ref{fig:fn_before_after}. The energy dependence in the ratio comes from two main sources~\cite{Szleper:2001nj,Kopp:2006ky}. First, the solid angle acceptance of the Near Detector varies as a function of the distance from the neutrino production point while the acceptance of the Far Detector is essentially constant. The two detectors subtend different solid angles and by Eq.~\ref{eq:angle} observe different neutrino angle distributions which, according to Eq.~\ref{eq:enu-eqn}, result in different neutrino energy distributions. Apertures in the beam-line (for example, the decay pipe) and under or over focusing in the horns enhance this effect by attenuating pions in certain energy ranges and locations along the beam-line. According to Eq.~\ref{eq:enu-eqn}, this introduces an energy dependent suppression of the neutrino flux which, because of the disparate angular acceptance, is not the same in the two detectors. Second, neutrinos produced by decays at finite radii from the beam axis will cover different angular ranges when intersecting the fiducial volumes of the Near and Far Detectors. This effect also introduces different energy distributions in the two detectors. These effects cause the peaks and dips in Fig~\ref{fig:fn_before_after}. Thus, the major sources of {\it relative} Near to Far differences in the neutrino flux (and hence, energy spectrum) are largely due to beam-line geometry, focusing, angular acceptance and decay kinematics. This suggests that the Monte Carlo simulation may be used, with relatively small uncertainties, to derive a transfer function that extrapolates the neutrino energy spectrum measured in the Near Detector to the Far Detector.

The Far to Near flux ratio of Fig.~\ref{fig:fn_before_after} is itself nearly a transfer function except that it is expressed in true neutrino energy and does not account for detector energy resolution, small disparities in efficiency between the two detectors, and their different fiducial masses.  A suitable replacement, however, can be evaluated by applying the \ensuremath{\nu_{\mu}} charged-current{} selection criteria of Sec.~\ref{sec:evtsel} to fully simulated events in both detectors to derive neutrino event rates, $\{n_{i}\}$ and $\{f_{i}\}$ in bins $i$ of reconstructed neutrino energy. The Near Detector data $\{N_{i}\}$ are then used to predict the Far spectrum:
\begin{equation}
F^{predicted}_{i}=N_{i}\times\frac{f_{i}}{n_{i}} \qquad \quad \mbox{``The F/N Method''}\quad .
\end{equation}
\noindent This technique is referred to as the ``F/N method'' but is equivalent to scaling each bin in the simulated Far Detector reconstructed neutrino energy spectrum by the ratio of the number of observed to expected events in the corresponding Near Detector reconstructed neutrino energy bin. The effect of neutrino oscillations may be accounted for by modifying $f_{i}$ according to Eq.~\ref{eq:oscprob}.

\begin{figure}
\centering
\subfloat[\label{fig:matrix_near}]{\includegraphics[width=0.45\textwidth]{fig/PrettyMatrixPlot-bw-near.eps}}
\subfloat[\label{fig:matrix_far}]{\includegraphics[width=0.45\textwidth]{fig/PrettyMatrixPlot-bw-far.eps}}
\caption{\label{fig:energy_spread} The relationship between the energy of neutrinos observed in the Near Detector with those observed in the Far Detector. Decays producing neutrinos with a given energy in the Near Detector would produce a range of energies in the Far Detector, yielding the energy smearing seen here. This relationship may be represented algebraically by treating the joint distribution in Fig.~\ref{fig:the_matrix} as a matrix, and the distributions \protect\subref{fig:matrix_near} and \protect\subref{fig:matrix_far} as column vectors as in Eq.~\ref{eq:matrix}. }
\end{figure}

\begin{figure}
\centering
\includegraphics[width=\columnwidth]{./fig/messier_matrix_plot.eps}
\caption{\label{fig:the_matrix} The joint distribution of neutrino energies observed in the Near and Far Detectors. The contents of each cell represent the mean number of \ensuremath{\nu_{\mu}}{} events expected in the Far Detector for one event in the Near Detector.  This distribution may be treated as a matrix, as in Eq.~\ref{eq:matrix}, to relate the energy spectra measured in the Near Detector to those in the Far Detector. }
\end{figure}

Neutrinos having a given energy in the Near Detector come from decays which would, collectively, yield neutrinos covering a range of energies in the Far Detector as indicated in Fig.~\ref{fig:near-far}, Eq.~\ref{eq:enu-eqn} and Eq.~\ref{eq:angle}. This effect is shown at several energies in Fig.~\ref{fig:energy_spread} and is again a consequence of the differences in angular acceptance of the two detectors, beam-line apertures, decay kinematics and off-axis decays. The fact that a single energy in the Near Detector corresponds to a range in the Far Detector suggests that the neutrino energy spectra may be related by a two-dimensional matrix rather than a one dimensional ratio.  This ``beam matrix'' $\{B_{ij}\}$ is shown in Fig.~\ref{fig:the_matrix}. Each cell represents the number of neutrinos expected in energy bin $i$ at the Far Detector for one neutrino in bin $j$ in the Near Detector. Neutrinos with energies in the range \unit[1-6]{GeV} are well focused and the contents of the cell along the diagonal is approximately proportional to the ratio of solid angles $\Omega_{FD}/\Omega_{ND} = d^{2}_{ND}/d^{2}_{FD}\approx(\unit[1]{km}/\unit[735]{km})^{2}\approx 1.8\times10^{-6}$. The matrix is constructed from the beam simulation using the known geometric acceptance of the two detectors. The energy dependence of the \ensuremath{\nu_{\mu}} charged-current{} cross-section{} is included in the calculation but is most relevant for the small off-diagonal elements. As with the F/N ratio, the Far Detector spectrum predicted by the matrix is relatively insensitive to uncertainties in the hadron production calculation, neutrino cross-sections and detector effects.

The matrix can only be employed after accounting for detector acceptance and inefficiencies. The Near Detector measurement is first corrected to remove non-  \ensuremath{\nu_{\mu}} charged-current{} contributions. A second correction, derived from a sample of fully simulated events in the Near Detector,  is then applied to deconvolve the detector efficiency and energy resolution. Both corrections are derived from the simulation described in Sec.~\ref{sec:mc}. The resulting distribution estimates the true \ensuremath{\nu_{\mu}} charged-current{} energy spectrum at the Near Detector. This distribution, organized in energy bins, is treated as an $m$-dimensional column vector $\{N_{i}\}$ and multiplied by the $m\times m$ dimensional matrix $\{B_{ij}\}$ to estimate the true \ensuremath{\nu_{\mu}} charged-current{} energy at the Far Detector
\begin{equation}
\label{eq:matrix}
F_{i}=\sum_{j=0}^{m} B_{ij} N_{j} \qquad \qquad \mbox{``The Beam Matrix Method''} \quad .
\end{equation}
The effects of neutrino oscillations may be applied to the derived Far Detector spectrum after which detector energy resolution, efficiency, and non-\ensuremath{\nu_{\mu}} charged-current{} backgrounds are added back in. The resulting spectrum may be directly compared to the Far Detector data.

The direct extrapolation methods described here reduce, but do not entirely eliminate, the effect that uncertainties in hadron production, neutrino cross-sections and detector acceptance have on the prediction of the neutrino energy spectrum at the Far Detector. We will discuss these uncertainties in greater detail in Sec~\ref{sec:extrap_syst}.

\subsubsection{Indirect Extrapolation Methods}

In these methods, a fit is performed to observed Near Detector distributions using a Monte Carlo dataset  with parametrized uncertainties due to beam modeling, neutrino cross-sections and neutrino energy scales. The result of this fit is a tuned Monte Carlo calculation which can then be applied to the Far Detector in order to obtain the predicted neutrino spectrum. These fits have the advantage that they can, in principle, separate systematic effects according to their effect on the fitted distributions with the disadvantage that the fitting functions may not adequately describe the data or (due to degeneracies in the fitted parameters) correctly describe the underlying physics.

Two indirect extrapolation methods have been developed: ``NDFit'' and ``2DFit''. In both techniques a fit is performed to the measured Near Detector energy spectra in order to improve agreement between the observed distributions and those predicted from Monte Carlo simulation. Both methods use data collected in the six beam configurations of Tab.~\ref{tab:beam-config} with approximately equal numbers of events in each configuration. In total, approximately $2.5\times 10^{5}$ Near Detector data events are used. The fitting techniques are permitted to re-extract the neutrino energy scale and normalization parameters originally determined when tuning the flux calculation (see Tab.~\ref{tab:fit_pars}). The two techniques employ very similar fitting procedures. In the case of the NDFit method, the fit is performed to one-dimensional reconstructed neutrino energy histograms, whereas in the case of the 2DFit, the reconstructed neutrino energy distributions are sub-divided in bins of inelasticity, $y = 1-E_{\mu}/E_{\nu}$.

The NDFit method varies the parameters $\boldsymbol{\alpha}=\{\alpha_{j}\}$ to minimize

\begin{equation}\label{eq:ndfit}
\chi^{2}=\sum_{k}^{E_{\nu}\mathrm{-bins}} \frac{\left(N_{k}-n_{k}\left(\boldsymbol{\alpha}\right)\right)^{2}}{\sigma_{k}^{2}+{S_{k}}^{2}} +\sum_{j=1}^{\mathrm{syst}}\left(\frac{\Delta\alpha_{j}}{\sigma^{\alpha}_{j} }\right)^{2} \qquad \qquad \mbox{``The NDFit Method''} \quad ,
\end{equation}

where $N_{k}$ and $n_{k}$ are the numbers of data and simulated events in bin $k$ of the reconstructed neutrino energy distribution. There are 38 energy bins for each of the six beam configurations (see Tab.~\ref{tab:beam-config} and Fig.~\ref{fig:6beamfit_all_enu}). The $\{\Delta\alpha_{j}\}$ are the deviations from nominal of the systematic parameters $\boldsymbol{\alpha}$, with associated uncertainties $\{\sigma^{\alpha}_{j}\}$. These are described below. The second term in Eq.~\ref{eq:ndfit} constrains the fit by increasing $\chi^{2}$ as the parameters are varied away from their nominal values. The error on the number of predicted events is the sum in quadrature of the statistical $\{\sigma_{k}\}$ and systematic $\{S_{k}\}$ uncertainties in the neutrino flux, where the latter are derived from the calculation described in Sec.~\ref{sec:skzp}.

 The five systematic parameters $\boldsymbol{\alpha}=\{\alpha_{1}\ldots\alpha_{5}\}$ describe uncertainties in the shape and normalization of the $\nu_\mu$ charged-current cross-section and the muon and hadron energy scales of the detector. These parameters are: $\alpha_1$, charged-current axial vector masses for quasi-elastic and resonance production processes, varied coherently; $\alpha_{2}$, cross-section scale factor for the non-resonant $1\pi$ and $2\pi$ production cross-section at invariant masses $W<\unit[1.7]{\ensuremath{\mathrm{GeV}/\mathit{c}^{2}}{}}$; $\alpha_3$, the absolute hadronic energy scale; $\alpha_4$, the absolute muon energy scale; $\alpha_5$, the overall normalization. The magnitude of these uncertainties is shown in Tab.~\ref{tab:ndfitres}.  The hadronic energy scale uncertainty is dominated by final state interactions (see Sec.~\ref{sec:mc} and \ref{sec:shwr_calib}). The muon energy scale was estimated from a comparison of range and curvature measurements for stopping muon tracks observed in the Near Detector. The normalization uncertainty accounts for the remaining uncertainties left after the improved flux calculation of Sec.~\ref{sec:skzp}.

The result of the NDfit is shown in Fig.~\ref{fig:ndfit}. Here the ratio of observed over simulated reconstructed neutrino energy spectra are shown before and after the fit for three of the six beam configurations. The large excursions in the ratio formed with the nominal simulation are substantially reduced by the improved hadron production calculation of Sec.~\ref{sec:skzp}. Further improvements are observed after minimizing Eq.~\ref{eq:ndfit}, especially in the 0-\unit[5]{GeV} region of the LE10/185kA configuration, and around the focusing peak of the LE100/200kA and LE250/200kA spectra. Initial excursions in the data/MC ratio of $\sim30\%$ have been reduced to $\sim 10\%$ or less in all six spectra. Remaining deviations are generally within the uncertainties of the hadron production calculation.

\begin{figure*}
\centering
\subfloat{\includegraphics[width=0.33\textwidth,clip]{fig/ndfit_summary_le.eps}}
\subfloat{\includegraphics[width=0.33\textwidth,clip]{fig/ndfit_summary_me.eps}}
\subfloat{\includegraphics[width=0.33\textwidth,clip]{fig/ndfit_summary_he.eps}}

\caption{\label{fig:ndfit} Ratio of Near Detector data to Monte Carlo simulation as a function of reconstructed neutrino energy for three beam spectra before and after the NDfit. The dotted gray histograms show the ratio relative to the nominal simulation, the dashed black lines show the ratio after applying the improved hadron production calculation (Sec.~\ref{sec:skzp}), and the black points with error bars show the ratio after the best-fit NDfit parameters are also applied to the simulation. The errors on the data points represent statistical uncertainties and the shaded regions represent the residual systematic uncertainty from the beam-re-weighting procedure.}
\end{figure*}

Table~\ref{tab:ndfitres} shows the best-fit values of the five systematic parameters used in the NDfit. All five parameters lie within one standard deviation of their nominal values. It is important to note that the results do not constitute a direct measurement of these quantities, as strong correlations exist between the fit variables. The value of $\chi^2$ at the best-fit point is 186.0 for 228 degrees of freedom, of which the penalty terms contribute four units of $\chi^{2}$ for five parameters.

\begin{table}[h]
\begin{center}
\begin{ruledtabular}
\begin{tabular}{clcc}
\multicolumn{2}{c}{Systematic parameter}  & Uncertainty &  Best-fit \\
\hline
$\quad \alpha_1$ & CC+RES axial-vector mass & $5\%$ & $+5\%$ \\

$\quad \alpha_2$ & non-resonant 1 and 2$\pi$ production & $20\%$ & $+12\%$\\

$\quad \alpha_3$ & Hadronic energy scale & $10\%$ & $-4\%$ \\

$\quad \alpha_4$ & Muon energy scale & $2\%$ & $-3\%$ \\

$\quad \alpha_5$ & Normalization  & $5\%$ & $-2.5\%$

\end{tabular}
\end{ruledtabular}

\caption{Systematic parameter values returned from the NDFit (Eq.~\ref{eq:ndfit}).\label{tab:ndfitres}}
\end{center}
\end{table}

The 2DFit method treats the Near Detector data as a function of both the reconstructed neutrino energy $\ensuremath{E_{\nu}}{}$ and inelasticity $y$ in order to allow degeneracies between muon and hadron energy scale parameters to be broken. In the fit, the joint distribution of reconstructed neutrino energy and $y$ is divided into a two dimensional grid with 40 energy bins spaced between \unit[0]{GeV} and \unit[30]{GeV} and eight reconstructed $y$ bins between zero and one. Data from each of the six beam configurations is tabulated separately, yieldings a total of 1920 analysis bins which are employed simultaneously in the fit. The fit attempts to minimize

\begin{equation} \label{eq:2dfit}
\chi^{2}=\sum_{k}^{E_{\nu}\mathrm{-bins}}\sum_{l}^{y\mathrm{-bins}} \frac{\left(N_{k,l}-n_{k,l}\left(\boldsymbol{\beta}\right)\right)^{2}}{\sigma_{k,l}^{2}} +\sum_{j=1}^{\mathrm{syst}}\left(\frac{\Delta\beta_{j}}{\sigma^{\beta}_{j} }\right)^{2} \qquad \qquad \mbox{``The 2DFit Method''} \quad ,
\end{equation}

\noindent where the variables $N_{k,l}$ and $n_{k,l}$ refer to the contents of the \ensuremath{E_{\nu}}{},$y$ bins for the selected \ensuremath{\nu_{\mu}} charged-current{} data and Monte Carlo samples respectively. Like the $\alpha_{j}$ in Eq.~\ref{eq:ndfit}, the $\beta_{j}$ are parameters which account for the systematic uncertainties in the Monte Carlo simulation and are varied in the fitting procedure to improve agreement with the data.

The $n_{k,l}$ are implicit functions of the parameters $\boldsymbol{\beta}$ which account for the following systematic uncertainties: $\beta_1$, the normalization of the neutral-current background; $\beta_2$, an energy dependent variation in the event rate dictated by the uncertainty in the tuned hadron production calculation; $\beta_3$, the absolute hadronic energy scale; $\beta_4$ the absolute muon energy scale; In addition, the true neutrino energy versus $y$ distribution for Monte Carlo events is divided into a two dimensional grid composed of eight neutrino energy bins (0-2, 2-4, 4-6, 6-8, 8-10, 10-15, 15-20, $>\unit[20]{GeV}$) and five $y$ bins (0.0-0.2, 0.2-0.4, 0.4-0.6, 0.6-1.0). This yields an additional 40 fit parameters, $\beta_5\ldots\beta_{44}$, which allow the normalization of the Monte Carlo to vary independently in each bin with an uncertainty of 20\%. This uncertainty is comparable to the {\em a priori} uncertainties in the neutrino flux and charged-current cross-sections. Table~\ref{tab:2dfitres} shows the magnitude of the systematic parameters. In total, 45 fit parameters are varied to minimize Eq.~\ref{eq:2dfit}.

\begin{table}[h]
\begin{center}
\begin{ruledtabular}
\begin{tabular}{clcc}
\multicolumn{2}{c}{Systematic parameter}  & Uncertainty &  Best-fit \\
\hline
$\quad \beta_1$ & Neutral-current background & $50\%$ & $-37\%$ \\

$\quad \beta_2$ & Flux uncertainty  & energy dependent  & $0.88\sigma$ \\

$\quad \beta_3$ & Hadronic energy scale & $10\%$ & $+2.5\%$\\

$\quad \beta_4$ & Muon energy scale & $2\%$ & $-4.8\%$ \\

$\quad \beta_5\ldots\beta_{44}$ & $(E_{\nu},y)$ normalization factors  & $20\%$ & 0.64--1.61
\end{tabular}
\end{ruledtabular}

\caption{Systematic parameter values returned from the 2DFit (Eq.~\ref{eq:2dfit}). The uncertainty on the flux parametrized by $\beta_2$ is energy dependent, approximately $\pm 9\%$ ($\pm 6\%$) for $E_{\nu} \lesssim \unit[6]{GeV}$ ($E_{\nu} \gtrsim \unit[6]{GeV}$). The fit parameter is expressed in units $\sigma$ of the uncertainty. The parameters $\beta_5\ldots\beta_{44}$ take values between 0.64 and 1.61 with a mean of 1.02 and an r.m.s. of 0.21, approximately equal to the {\it a priori} uncertainty. \label{tab:2dfitres}}
\end{center}
\end{table}

Figure~\ref{fig:2dfit_nd_le10} shows the neutrino energy distribution from the LE10~\unit[185]{kA} configuration split into four ranges of $y$. The value of $\chi^{2}$ at the best-fit point is 2606.3 for 1919 degrees of freedom, compared to a $\chi^{2}$ of 5432.1 before the fit. The best-fit values of the 40 $(E_{\nu},y)$ normalization factors range from 0.64 to 1.61, with a mean of 1.02 and an r.m.s. of 0.21 (approximately equal to the {\it a priori} uncertainty). The best-fit values for the other parameters are -37\% for neutral-current normalization and -4.8\% and +2.5\% for muon and shower energy scales, respectively. The overall normalization was reduced according to a $0.88\sigma$ shift in the tuned hadron production calculation.

\begin{figure}[h]
\begin{center}
  \includegraphics[width=\columnwidth]{fig/2dfit_4plots.eps}
  \caption{Near Detector reconstructed neutrino energy spectra collected in the LE10/185kA beam configuration and for different reconstructed $y$ regions. Open circles show the data. The Monte Carlo simulation is shown before and after the tuning procedure performed by the 2DFit technique (Eq.~\ref{eq:2dfit}). To simplify the presentation we have grouped the eight reconstructed $y$ bins used in the fit into four bins here.}
   \label{fig:2dfit_nd_le10}
\end{center}
\end{figure}

\subsection{Predicted Far Detector spectrum}

Figure~\ref{fig:4spectra} shows a comparison of the predicted Far Detector unoscillated visible energy spectra obtained from the four extrapolation methods, for an exposure of $1.27\times10^{20}$ protons on target. The ratios of the three other spectra relative to the Beam Matrix prediction are shown in Fig.~\ref{fig:4ratios}. The agreement between the predicted spectra is better than $\pm 5\%$ in the region 1-\unit[15]{GeV}, where the oscillation signal is expected to lie. The spread on the predicted spectra is significantly smaller than the statistical error when the data are binned in \unit[1]{GeV} wide energy bins. Table~\ref{tab:fdprednum} shows the number of Far Detector events predicted between 0-30 GeV for each extrapolation method. The systematic errors on the predicted rate, which are described below, are dominated by the $4\%$ relative normalization error.

\begin{figure}[h]
\begin{center}
\subfloat[\label{fig:4spectra}]{\includegraphics[width=0.7\columnwidth,clip]{fig/4methods_spectra_1.eps}}\\
\subfloat[\label{fig:4ratios}]{\includegraphics[width=0.7\columnwidth,clip]{fig/4methods_spectra_2.eps}}
\caption{\label{fig:spectra_4methods} Neutrino energy spectra at the Far Detector in the absence of neutrino oscillations as predicted by the four extrapolation methods.}
\end{center}
\end{figure}

\begin{table*}[]
\begin{center}
\begin{ruledtabular}
\begin{tabular}{cc}
Extrapolation Method & Predicted Number of Far Detector Events \\
\hline

Beam Matrix & $336.0\pm 14.4$ \\

Far/Near Ratio & $330.0\pm 14.2$ \\

NDFit & $332.7\pm 14.3$ \\

2DFit & $323.9\pm 13.9$ \\

\end{tabular}
\end{ruledtabular}

\caption{Predicted numbers of Far Detector events for the four extrapolation methods for an exposure $\unit[1.27\times10^{20}]{POT}$. The systematic errors on the predictions, which are dominated by the 4\% relative normalization uncertainty, are also shown.\label{tab:fdprednum}}
\end{center}
\end{table*}

\subsection{Sensitivity to systematic errors for the different extrapolation methods}

\label{sec:extrap_syst}

In this section we describe the various sources of systematic uncertainty and their impact on the neutrino spectrum predicted at the Far Detector.  Each source was considered separately and Monte Carlo datasets were generated with $\pm1\sigma$ shifts of the systematic parameter applied along with oscillations according to $\ensuremath{\Delta m^2}=\unit[2.72\times10^{-3}]{\ensuremath{\mathrm{eV}^{2}/\mathit{c}^{4}}{}}$ and $\ensuremath{\sin^{2}2{\theta}}=1$, chosen to agree with our best fit to the data.  These ``mock datasets'' were produced for both Near and Far Detectors and were analyzed in the same way as the actual data. In particular the tuning procedure of Sec.~\ref{sec:skzp} was run on each dataset and the extrapolation methods were able to use the Near Detector mock data to predict the spectra at the Far Detector. These predictions were then used to fit the Far Detector mock data for neutrino oscillations, allowing us to examine the way in which systematic uncertainties affect our measurements of \ensuremath{\Delta m^2}{} and \ensuremath{\sin^{2}2{\theta}}{}. The largest sources of uncertainty were identified with this procedure and then accounted for when fitting the Far Detector data for neutrino oscillations.

\subsubsection{Systematic error sources}

Figure~\ref{fig:systnceshw} illustrates the effect that four representative sources of uncertainty have on the Far Detector neutrino energy spectrum, and the way in which the extrapolation methods are able to account for them with Near Detector measurements. Black points show the ratio of reconstructed energy spectra at the Far Detector, where the denominator corresponds to the nominal simulation and the numerator was generated with a particular systematic uncertainty source shifted with respect to the default value. Deviations from unity indicate the way in which the systematic uncertainties modify the neutrino energy spectrum. The size of the deviations does not necessarily indicate the potential effect on \ensuremath{\Delta m^2}{} and \ensuremath{\sin^{2}2{\theta}}{} because the extrapolation methods use similarly shifted Near Detector Monte Carlo samples as ``mock data''. Their prediction of the ratio, shown as lines, would follow the black points if the systematic shift could be completely corrected with Near Detector data. In that case one would extract the same oscillation parameters as were put into the study and the particular source of uncertainty would not affect the results. We describe the sources of uncertainty below.

{\em Charged-current cross-sections:} We assign a $\pm10\%$ uncertainty in the axial vector masses for quasi-elastic and resonance production processes. Fig.~\ref{fig:syst_ccma} shows the effect on the FD spectrum when a $+10\%$ shift in those quantities is applied in the simulation. The predicted spectra from the four extrapolation methods successfully reproduce the shift.  This is expected, since cross-section changes are common to both Near and Far detector events, the Near to Far extrapolation methods should provide a significant cancellation of such uncertainties. We have also conducted a similar study to verify that our analyses are able to correct for the $\pm20\%$ uncertainty in the magnitude of the non-resonant $1\pi$ and $2\pi$ production cross-sections for invariant masses $W<\unit[1.7]{\ensuremath{\mathrm{GeV}/\mathit{c}^{2}}{}}$. See Sec.~\ref{sec:mc} for further details.

\begin{figure}[h]
\centering
\subfloat[\label{fig:syst_ccma}]{\includegraphics[width=0.45\textwidth]{./fig/systematics_test_ccma.eps}}
\hfill
\subfloat[\label{fig:syst_eshw}]{ \includegraphics[width=0.45\textwidth]{./fig/systematics_test_showerenergy.eps}}
\\
\subfloat[\label{fig:syst_skzp}]{\includegraphics[width=0.45\textwidth]{./fig/systematics_test_skzp.eps}}
\hfill
\subfloat[\label{fig:syst_nc}]{\includegraphics[width=0.45\textwidth]{./fig/systematics_test_ncbackground.eps}}

 \caption{Effect of systematic uncertainties on the predicted Far Detector energy spectrum. The following systematic shifts are shown: \protect\subref{fig:syst_ccma} $+10\%$ increase in \ensuremath{\nu_{\mu}}{} charged-current quasi-elastic and resonance production cross-sections;  \protect\subref{fig:syst_eshw} $-10\%$ shift in reconstructed shower energy scale; \protect\subref{fig:syst_skzp} $+50\%$ uncertainty on the beam tuning correction; \protect\subref{fig:syst_nc} $+50\%$ increase in neutral-current background rate. Further details are described in the text. \label{fig:systnceshw}}
\end{figure}

{\em Final-state interactions:} As described in Sec.~\ref{sec:mc} we estimate that uncertainties in final-state interactions cause a  $\pm10\%$ uncertainty on the shower energy scale. Figure~\ref{fig:syst_eshw} shows the effect that a $-10\%$ shift in the energy scale has on the Far Detector neutrino spectrum. The distortion of the spectrum below \unit[10]{GeV} is significant. Because the same shift is present in both Near and Far Detectors, the predicted spectra from the four methods provide an improved description of the distorted spectrum and thus reduce the effect of this systematic error. However, the neutrino energy and \ensuremath{\Delta m^2}{} are coupled in Eq.~\ref{eq:oscprob} and a change in the energy scale may be incorrectly interpreted as a shift in \ensuremath{\Delta m^2}{}. It is therefore necessary to explicitly account for energy scale uncertainties when extracting the oscillation parameters from our data.

{\em Beam-related uncertainties:} The beam tuning procedure described in Sec.~\ref{sec:skzp} applies corrections to the (\ensuremath{\pz{},\pt{}}{}) distribution in the Monte Carlo that result in changes in the visible energy spectrum of approximately $10\%$ around the peak in the LE10/185kA beam configuration and about $30\%$ in the high energy tail. We have conservatively assigned a $50\%$ error on the magnitude of this correction as an estimate of the uncertainty on the beam tuning. The effect of this uncertainty on the FD spectrum is shown in Fig.~\ref{fig:syst_skzp}. There is a significant shift in the high-energy tail of the spectrum ($\ensuremath{E_{\nu}}{} > \unit[5]{GeV}$) where the beam tuning correction is largest. The plot shows that the direct extrapolation methods are in general more robust than the fitting methods to this type of uncertainty, as they are much less dependent on the requirement that the tuned Monte Carlo closely match the data in order to provide an accurate prediction of the FD spectrum.

{\em Normalization:} We estimate an uncertainty of $\pm4\%$ in the relative normalization of the energy spectra measured in the Near and Far Detectors. The estimate is composed from a $2\%$ uncertainty in the fiducial mass of both detectors, a $3\%$ uncertainty in the relative Near and Far Detector reconstruction efficiencies, estimated from a visual scan of Near and Far Detector data and Monte Carlo events, and a $1\%$ uncertainty in the live time.

{\em Neutral-current background:} As described in Sec.~\ref{sec:nc}, we estimate an uncertainty of $\pm 50\%$ on the rate of neutral-current events misidentified as \ensuremath{\nu_{\mu}}{} charged-current. Figure~\ref{fig:syst_nc} shows that the effect of this uncertainty is largest in the lowest reconstructed energy bins ($E_\nu <\unit[2]{GeV}$) and is only partially corrected by the F/N, Beam Matrix and NDFit extrapolation methods. The 2DFit method includes the neutral-current rate as a fit parameter and recovers the correct prediction in this special case. In the general case however, where several systematic effects will be present, no method will be able to completely correct for this source of uncertainty.

{\em Shower energy scale calibration:} Based on a comparison of test beam measurements in the Calibration Detector and the Monte Carlo simulation (see Sec.~\ref{sec:shwr_calib}) we estimate a $6 \%$ uncertainty in the response to single pions and protons. The uncertainty in the response to electrons and muons is negligible. As described in Sec.~\ref{sec:calib}, stopping muons were used to cross-calibrate the two detectors with a relative uncertainty of $\pm 2\%$.

{\em Muon energy scale:} A $2\%$ muon energy scale uncertainty is assumed based on studies of fully contained muon tracks reconstructed in the Near Detector. The difference between the momentum obtained from the track range and the momentum obtained from a fit to the curvature of the track due to the magnetic field of the detector was examined for both real and simulated events, as a function of the range momentum. The deviation of this quantity between data and simulation was approximately 2\%, which was taken as an estimate of the uncertainty on the magnetic field calibration in the Near Detector, and hence the error on the relative muon energy scale between Near and Far Detectors. Our estimate utilizes the muon range/energy relation which has an uncertainty of approximately 2\% based on a material assay of the detectors and a comparison of our muon stopping power tables with those in~\cite{Groom:2001kq}.

As discussed in Sec.~\ref{sec:mc}, our magnetic field was recalibrated after we performed this analysis. This predominantly affects the reconstruction of muons with a momentum larger than $\sim \unit[7]{\ensuremath{\mathrm{GeV}/\mathit{c}}{}}$ which arise from neutrinos outside the energy region in which we observe oscillations.

{\em CC selection efficiency:} We varied the requirement on the event selection parameter $S$ (see Eq.~\ref{eq:S}) for mock data sets by $\pm 0.02$ while holding the cut applied in the nominal Monte Carlo dataset constant. This changes both the number of true charged-current and neutral current events that are classified as \ensuremath{\nu_{\mu}}{} charged-current. The magnitude of this shift was obtained from a comparison of the data and Monte Carlo $S$ distributions in the Near Detector.

We used mock datasets including the systematic uncertainties listed above to examine the accuracy with which the different methods are able to predict the Far Detector spectrum and study the \ensuremath{\sin^{2}2{\theta}}{} and \ensuremath{\Delta m^2}{} extracted from our oscillation fits. As described below, this procedure was used to explore the capabilities of the different extrapolation techniques and identify the most important sources of systematic uncertainty.

\subsubsection{Choice of primary analysis method}

 The indirect extrapolation methods (NDFit, 2DFit) attempt to adjust the Monte Carlo prediction of the neutrino energy spectrum in the Near Detector to improve agreement with the data and then use the improved Monte Carlo calculation to predict the Far Detector spectrum. These procedures are not able to arrive at a perfect description of the data. In particular, discrepancies of $\pm5\%$ are present when data collected in the LE10/185kA beam configuration are compared with the Monte Carlo prediction after the best-fit NDFit and 2DFit systematic parameters have been applied. We also notice that somewhat larger discrepancies are still present in the LE100/200kA beam configuration. If the major source of these distortions is caused by uncertainties in the neutrino flux or neutrino cross-sections{} then the direct (Beam Matrix and the F/N) extrapolation methods simply translate the measurements to the Far Detector and the discrepancies do not affect the oscillation measurement. As a consequence we have found that the predictions made by the direct methods are less sensitive to the absolute agreement between data and the Monte Carlo simulation in the Near Detector than are the those made by the indirect methods.

In light of these issues, we have decided to use one of the direct extrapolation methods to obtain our primary oscillation result. We chose the Beam Matrix method as our primary extrapolation technique because it had smaller systematic errors than the F/N ratio method when all sources of uncertainty were considered.  The results that we will present in Sec.~\ref{sec:osc} therefore use the Beam Matrix method to predict the Far Detector spectrum and perform the oscillation fit. As a cross-check of these results, we present the best fit oscillation parameters and allowed regions obtained from the other three extrapolation methods as well.

\subsubsection{Systematic Uncertainties in the Oscillation Fit}

The three largest contributions to the systematic error for the Beam Matrix method are a) the uncertainty in the relative normalization of the energy spectra measured in the two detectors, b) uncertainties in the absolute hadronic energy scale, and c) uncertainties in the neutral current background rate. The systematic shifts calculated for the other sources of uncertainty are small in comparison. The magnitude of these shifts is summarized in Tab.~\ref{tab:oscsyst}. As expected from the above discussion, because uncertainties due to beam modeling and cross-sections are common to the two detectors they largely cancel out in the extrapolation. The three largest uncertainties are included as systematic nuisance parameters in the Far Detector oscillation fit. By fitting for these systematic parameters simultaneously using both Near and Far Detector data, the effect of these uncertainties is substantially reduced, due to significant cancellations of these errors between the two detectors.

\subsection{Calculation of systematic errors on the predicted FD spectrum}

A procedure to identify the most important sources of systematic uncertainty on the measurement of $\ensuremath{\Delta m^2}$ and $\ensuremath{\sin^{2}2{\theta}}$ was carried out using fake datasets generated from Monte Carlo. These were produced for both Near and Far detectors and were analysed in the same way as the real data. Each source of systematic error was considered separately and fake datasets were generated with $\pm1\sigma$ shifts of the systematic parameter applied and, in the case of the simulated FD data, oscillations with $\ensuremath{\Delta m^2}=2.72\times10^{-3}$~eV$^{2}$ and $\ensuremath{\sin^{2}2{\theta}}=1$.

\vspace{0.1in}
\par

The systematic error sources considered were the following:

\begin{itemize}

\item Neutral Current Background - estimated uncertainty $\pm50\%$ from a comparison of the shape and normalisation of the observed Near detector data PID distribution to the expectation from Monte Carlo, as described above.

\item Intranuclear Re-scattering - estimated uncertainty of $\pm10\%$ on the absolute shower energy scale, which is based on the change in the visible shower energy spectrum from the Monte Carlo, when the limiting cases of `full pion absorption' and no intranuclear re-scattering are considered.

\item Absolute Normalization - estimated uncertainty of $\pm4\%$, based on uncertainties in fiducial volume uncertainties ($\sim2\%$), relative Near/Far reconstruction efficiencies (estimated from a visual scan of ND and FD data and Monte Carlo - $\sim 2-3\%$)  and POT counting ($\sim 1\%$)

\item Beam Tuning - The beam tuning procedure described in the previous section applies corrections to the $\ensuremath{x_{F}},\ensuremath{p_{T}}$ distributions in the Monte Carlo that result in changes in the $LE$ visible energy spectrum of $\sim10\%$ around the peak, and $\sim30\%$ in the high energy tail. We have conservatively assigned a $50\%$ error on the magnitude of this correction as an estimate of the uncertainty on the beam tuning procedure.

\item Charged-current cross-sections - uncertainties of $\pm10\%$ in the axial vector masses for quasi-elastic and resonance production processes are assigned. Uncertainties of $\pm20\%$ are applied to the cross-sections of low-multiplicity DIS events in the Resonance-DIS overlap region. See the previous section on modelling neutrino interactions for details of these error estimates.

\item Shower energy scale calibration - A $\pm6\%$ uncertainty on the absolute shower energy scale is assumed, based on a comparison of test beam measurements in the Calibration detector and the Monte Carlo simulation. Also, a $\pm2\%$ relative uncertainty between Near and Far detectors is assumed, based on uncertainties in calibrating the MINOS detetors using stopping muons.

\item Muon energy scale - assumed $2\%$ relative muon energy scale uncertainty based on studies of fully contained muon tracks recontructed in the Near Detector. The difference between the momentum obtained from the track range and the momentum obtained from a fit to the curvature of the track due to the magnetic field of the detector was examined for both data and MC events, as a function of the range momentum. The deviation of this quantity between data and Monte Carlo was observed to be less than 2\%, which was taken as an estimate of the uncertainty on the magnetic field calibration in the Near Detector, and hence the error on the relative muon energy scale between Near and Far.

\item CC selection efficiency - applied a shift of $\pm0.02$ in the PID cut that is applied to the fake data samples. This changes both the number of true charged-current and neutral current events that are classified as CC-like. The magnitude of this shift was obtained from a comparison of the data and Monte Carlo PID distributions in the Near Detector.

\end{itemize}

\subsubsection{Sensitivity to systematic errors for the different extrapolation methods}

The main advantage of the direct extrapolation methods (Beam Matrix and F/N ratio) is that the two largest sources of error on the predicted spectrum, due to
cross section uncertainties and beam modelling (mainly hadron
production) related uncertainties cancel out when the ND data is used directly to predict the FD spectrum, without the need for prior tuning of the Monte Carlo to better match the measured ND spectrum.
This is clearly illustrated in figure~\ref{fig2} which shows the extent to which beam and cross-section related uncertainties are reduced using the Beam matrix extrapolation method. Fake datasets were prepared with $\pm1\sigma$ shifts in the beam model (left-hand plot) and cross-section uncertainties (right-hand plot), using the Monte Carlo procedure described above. The black lines in the plot show the ratio of the visible energy spectra of these shifted datasets relative to the nominal Monte Carlo in the Near Detector. These shifted spectra are then used to predict the unoscillated spectra in the Far Detector, and the red histograms show the resulting ratio of the predicted FD spectra to the expected FD spectrum (from the FD fake dataset that is generated with the same systematic shift that is applied to the ND fake dataset). The plot shows that the large initial differences between fake data and Monte Carlo are substantially reduced using this extrapolation technique.

Figure~\ref{fig:systndfit} shows the same comparison as Figure~\ref{fig2}, but for the NDFit method. The right-hand panel of this plot shows that the NDFit extrapolation is quite successful in reducing the effect of cross-section errors on the extrapolated spectrum. However, the left-hand plot, which is for beam-related uncertainties and does not include a prior beam-tuning step in the extrapolation, shows less complete cancellation of systematic errors than the direct extrapolation methods, and illustrates that Monte Carlo tuning is a much more important component of the fit methods.

With additional information that will enable us to more cleanly factorise beam and cross-section related uncertainties, the fit methods can potentially provide more powerful constraints on the systematic errors. However, given that our current largest uncertainties on the predicted spectrum are due to beam and cross-section errors, we have decided to use one of the direct extrapolation methods to obtain our primary oscillation result. We chose the Beam Matrix method as our primary extrapolation technique, since it was shown to be more robust to systematic errors that the F/N ratio method, using the fake data MC technique that we outlined above.

\begin{figure}[h]
\centering
\begin{tabular}{cc}

\includegraphics[width=0.45\columnwidth]{./fig/farnear_ratio_beam.eps}

& \includegraphics[width=0.45\columnwidth]{./fig/farnear_ratio_xsec.eps} \\

\end{tabular}
 \caption{Monte Carlo evaluation of the sensitivity of the Beam Matrix method to beam and cross-section related  uncertainties. See accompanying text for details  \label{fig2} }
\end{figure}

\begin{figure}[h]
\centering
\begin{tabular}{cc}

\includegraphics[width=0.45\columnwidth]{./fig/NDFD_SKZP2Par0.eps}

& \includegraphics[width=0.45\columnwidth]{./fig/NDFD_CCMA.eps} \\

\end{tabular}
 \caption{As Figure~\ref{fig2} but for the NDFit extrapolation method. \ecom{The current plot on the left is not exactly ideal, but can serve placeholder figure until a better alternative can be found to illustrate that the fit methods require some form of beam tuning to reduce beam-related systematic errors. Also need to make the line styles for this plot and figure ~\ref{fig2} consistent}} \label{fig:systndfit}
\end{figure}

\subsubsection{Leading Systematic Uncertainties}

As a result of this study, the three largest contributions to the systematic error for the Beam Matrix method were found to be a) the 4\% uncertainty in the overall normalisation; b)  uncertainties in the absolute hadronic energy scale due to intranuclear rescattering, and c) uncertainties in the neutral current background rate. The systematic shifts calculated for the other sources of uncertainty were found to be small in comparison. Table~\ref{tab:oscsyst} shows the magnitudes of the systematic shifts obtained from this procedure. As expected from the above discussion, uncertainties due to beam modelling and cross-sections (entries d) and e) in the table) - which are common to the two detectors - largely cancel out in the extrapolation. The three largest uncertainties (a, b, and c) are included as systematic nuisance parameters in the Far Detector oscillation fit, which is described in the following section. By fitting for these systematic parameters simultaneously using both Near and Far detector data, the effect of these uncertainties are substantially reduced, due to significant cancellations of these errors between Near and Far.

\begin{table}[h]
\begin{center}
      \begin{tabular}{llcc}
	\hline
	\hline
	Uncertainty &&\ensuremath{|\Delta m^{2}|} & \ensuremath{\sin^{2}2{\theta}} \\
	  && $(10^{-3}\,\mathrm{eV}^2)$ & \\
	\hline
	(a) Normalization & $(\pm 4\%)$          & 0.050 & 0.005\\
	(b) Abs. hadronic E scale & $(\pm 11\%)$ & 0.057 & 0.048\\
	(c) NC contamination &$(\pm 50\%)$       & 0.090 & 0.050\\
	\hline
	(d) Beam uncertainties &    & 0.015 & $<$0.005\\
	(e) Cross sections  &       & 0.011 & 0.005\\
	All other systematics &                  & 0.041 & 0.013\\
	\hline
Statistical Error &   & 0.35 & 0.13\\
        \hline
	\hline
      \end{tabular}
\caption{Systematic shifts on the measurement
	of \ensuremath{|\Delta m^{2}|}{} and \ensuremath{\sin^{2}2{\theta}}{} for various sources of systematic error, calculated using Monte Carlo generated with the best-fit oscillation parameters obtained from the FD data.
	The values quoted are the average shifts for $\pm$ 1 standard deviation variations in each of the systematic parameters.
        The last row of the table shows the expected statistical uncertainty on the measurement of \ensuremath{\Delta m^2} and \ensuremath{\sin^{2}2{\theta}} for an exposure of $1.27\times 10^{20}$~POT.
	The shifts on \ensuremath{|\Delta m^{2}|} and \ensuremath{\sin^{2}2{\theta}} are treated as uncorrelated, and correlations between the various systematic effects are not taken
	into account.\label{tab:oscsyst}}
\end{center}
\end{table}

\subsection{Modeling Neutrino Interactions}
\label{sec:mc}

\begin{figure}
\centering
\includegraphics[width=\columnwidth]{fig/sig_prl}
\caption{\label{fig:ng_xsec}{\tt NEUGEN-v3}{} calculation of  muon-neutrino charged-current cross-sections{} per nucleon on an isoscalar target.  The cross-section{} per GeV is shown as a function of the neutrino energy for inclusive scattering, quasi-elastic scattering and single pion production.  The calculation is compared with experimental data tabulated by~\cite{durham}.  The shaded band corresponds to the cross-section{} uncertainties described in the text.}
\end{figure}

Neutrino interactions are modeled by the {\tt NEUGEN-v3}{}~\cite{Gallagher:2002sf} program. {\tt NEUGEN}{} simulates both quasi-elastic and inelastic neutrino scattering. The latter includes a Rein-Sehgal~\cite{Rein:1980wg} based treatment of neutrino induced resonance production, charged- and neutral-current coherent pion production and a modified leading order deep-inelastic scattering (DIS) model~\cite{Bodek:2004ea} extended to improve the treatment in the transition region between DIS and resonant production. KNO scaling~\cite{Koba:1972ng} is used to calculate the final state multiplicity in the DIS regime. Hadrons produced in the neutrino scattering are allowed to interact while exiting the target nucleus (``final state interactions''). The final state interaction calculation incorporates pion elastic and inelastic scattering, single charge exchange and absorption~\cite{Ransome:2005vb}.  The calculation is benchmarked by a comparison of final states in $\nu + d$ and $\nu + \mathrm{Ne}$ interactions as measured in the BEBC and ANL--\unit[12]{ft} bubble chambers~\cite{Merenyi:1992gf}.

Figure~\ref{fig:ng_xsec} shows the \ensuremath{\nu_{\mu}} charged-current{} cross-section{} as a function of neutrino energy in the laboratory frame. Based on a comparison of the model predictions to independent data, some of which~\cite{durham} is shown in Fig.~\ref{fig:ng_xsec}, we assign a systematic uncertainty of 3\% on the normalization of the DIS ($W>\unit[1.7]{\ensuremath{\mathrm{GeV}/\mathit{c}^{2}}{}}$) cross-section{}, and a 10\% uncertainty in the normalization of the single-pion and quasi-elastic cross-sections{}. We estimate a 20\% uncertainty in the relative contribution of non-resonant states to the 1$\pi$ and 2$\pi$ production cross-sections{} for $W<\unit[1.7]{\ensuremath{\mathrm{GeV}/\mathit{c}^{2}}{}}$. This uncertainty was determined from the parameter uncertainties and variations observed in fits to both inclusive and exclusive channel data, and in fits to data in different invariant mass regions. Final state interactions are expected to have a significant effect on the visible energy of the hadronic final state~\cite{Kordosky:2006gt}.  In particular there are significant uncertainties in the rate of pion absorption, the mechanism for transferring the pion's energy to a nucleon cluster, and the amount of energy eventually visible to the detector. We account for these uncertainties by studying the shift in the reconstructed shower energy when we turn the simulation of final-state interactions off, and when we modify the simulation so that all of an absorbed pion's energy is lost. We find that the predicted response to hadronic showers changes by approximately 10\%~\cite{Kordosky:2006gt} in these two extreme cases and use this as a conservative estimate of the uncertainty on the absolute hadronic energy scale.

 The MINOS{} detector simulation is based on {\tt GEANT3}{}~\cite{g3man} and is used to generate raw energy depositions ({\tt GEANT}{} hits) which serve as the input to our detector response model. The simulation randomly samples neutrinos from the flux predicted by the beam simulation (Sec.~\ref{sec:beam_sim}) and traces them through the Near and Far Detector halls. Events are generated inside the detectors as well as in the surrounding support structure and rock. The simulation includes a detailed geometric model of the detector which describes the material crossed by neutrinos and neutrino induced tracks to within \unc{1\%} plane-to-plane and \unc{0.3\%} averaged over the detector. The position of individual scintillator strips was determined with a precision of approximately \unit[1]{mm} using cosmic-ray tracks. The magnetic field is modeled via finite element analysis driven by bench measurements of the steel B--H curve. Since performing this analysis we have recalibrated our field and found that it increased by 12.3\% and 9.2\% averaged over the fiducial volumes in the Near and Far Detectors, respectively. This recalibrated field shifts the momentum scale for muons exiting the Near Detector by approximately 6.2\% (4.6\% in the Far Detector) but does not significantly affect the scale for muons which stop in the detector.

The detector response simulation has been tuned by directly incorporating the channel-by-channel calibration constants determined by the detector calibration procedure (see Sec.~\ref{sec:calib} below). Additional data, measured in pre-installation tests, were used to constrain details of the scintillator and photomultiplier response models. The simulation has been benchmarked against through-going and stopping cosmic rays in the Near and Far Detectors, beam neutrino events in the Near Detector, as well as a series of test-beam measurements~\cite{mike_thesis,tricia_thesis,jeff_thesis,anatael_thesis} collected with a scaled down version of the MINOS{} detectors in the CERN PS East Hall. The test-beam measurements are used to fix the energy scale and validate the simulation of electromagnetic and hadronic showers. The hadronic shower code {\tt GCALOR}{}~\cite{Zeitnitz:1994bs} was found to be in relatively good agreement with the data (see Sec.~\ref{sec:shwr_calib}) and is used in our simulations.

\subsubsection{Initial Selection of Muon Neutrino Induced Events}

The total reconstructed energy of each event is estimated by summing the energy of the most energetic track with the energy of any shower present at the upstream end of the track. We select \ensuremath{\nu_{\mu}} charged-current{} interactions by requiring that the event has at least one well reconstructed track with a starting-point, interpreted as the neutrino interaction point (vertex), in the fiducial volume. In the Near Detector, the fiducial volume is a cylinder of radius \unit[1]{m} from the beam center~\footnote{The radius is calculated from the beam center in the first detector plane. The cylinders follow the z axis of the detectors which differs from the beam axis by \unit[58]{mrad}.} and length \unit[4]{m} beginning \unit[1]{m} downstream of the front face of the detector. The Far Detector fiducial volume is a cylinder of radius \unit[3.7]{m} from the detector center. Vertices are required to be at least \unit[50]{cm} from the front and rear planes of the two detector modules and, to assure a track that is long enough to be analyzed, greater than \unit[2]{m} upstream of the last plane of module two. Tracks are required to have a negative charge to suppress \ensuremath{\overline{\nu}_{\mu}}{} and only events with neutrino energy less than \unit[30]{GeV} are used so as to preferentially select neutrinos from  $\pi$, rather than $K$, decays.  A final set of criteria are applied to remove rare periods in which the magnetic field coil was not energized or in which the high voltage was not on in some portion of the detector.

The resulting sample in the Near Detector is 89\% pure \ensuremath{\nu_{\mu}} charged-current{} (91\% in the Far Detector, assuming no neutrino oscillations) with the dominant background coming from neutral-current events in which a (usually short) track was reconstructed. Section~\ref{sec:evtsel} describes the way in which event-topological quantities were used to improve the sample purity. In the following subsections we describe some characteristics of the dataset before purification.

\subsubsection{Near Detector Data}
\label{sec:nd_dist}

We have reconstructed more than $2.6\times10^{6}$ neutrino events in the Near Detector in several different beam configurations (see Tab.~\ref{tab:beam-config}). These data have been used to verify and improve the quality of our neutrino interaction simulation, detector calibration, event reconstruction, neutral current rejection procedure (described in Sec.~\ref{sec:evtsel}), and beam flux calculation (Sec.~\ref{sec:skzp}). Since we are ultimately interested in predicting the Far Detector energy spectrum, it is not necessary that the Monte Carlo simulation reproduce the Near Detector data exactly. Instead we must have confidence that any discrepancies have causes which are common to both detectors such that the Near Detector data may be used to improve the Far Detector prediction. For example, uncertainties in the neutrino flux, neutrino cross-sections{} and energy resolutions (though not the absolute energy scale) may be mitigated by using Near Detector measurements.

\begin{figure}
\centering
\subfloat[]{\includegraphics[width=0.4\textwidth,clip]{fig/prd_nd_lowlev_1.eps}}
\\
\subfloat[]{\includegraphics[width=0.4\textwidth,clip]{fig/prd_nd_lowlev_2.eps}}
\\
\subfloat[]{\includegraphics[width=0.4\textwidth,clip]{fig/prd_nd_lowlev_3.eps}}
 \caption{Distributions of basic reconstructed quantities in the Near Detector for events reconstructed and selected using the procedure described in the text.  The distributions are normalized to the same numbers of events. The left-hand and center plot show the distribution of reconstructed track vertices for $X$ (perpendicular to detector axis) and $Z$ (parallel to detector axis) detector co-ordinates, and the right-hand plot shows the distribution of the reconstructed track direction relative to the incident neutrino beam. \label{fig:nd_lowlevel}}
\end{figure}

The Monte Carlo simulation reproduces many reconstructed quantities in the Near Detector, including the vertex distributions and track angular distribution shown in Fig.~\ref{fig:nd_lowlevel}. Figure~\ref{fig:events_per_ppp} shows the mean number of reconstructed events as a function of beam intensity ranging between $5\times 10^{11}$ and \unit[$2.7\times 10^{13}$]{POT/spill}. The linearity of the curve indicates that the Near Detector is able to measure individual neutrino interactions with little dependence on the interaction rate. The behavior of the curve near the origin suggests a negligible background from non-beam related events. The reconstructed neutrino energy spectrum measured during several months is overlaid in Fig.~\ref{fig:enu_by_month}. The average spill intensity in June 2005 is \unit[$1.65\times 10^{13}$]{POT/spill} and rises to a maximum of \unit[$2.4\times 10^{13}$]{POT/spill} in November 2005. The distributions are consistent within statistical errors and show no rate dependence. We therefore conclude that the Near Detector event reconstruction is stable over this range of spill intensities, and does not introduce any observable biases in the reconstructed neutrino energy spectrum. Consequently, we do not assign any intensity based uncertainty in the detector efficiency.

\begin{figure}[h]
\centering
\includegraphics[width=\columnwidth]{./fig/nevent_vs_ppp.eps}
\caption{\label{fig:events_per_ppp} Mean number of reconstructed events per Near Detector spill as a function of spill intensity.}
\end{figure}

\begin{figure*}[h]
\centering
\includegraphics[width=\textwidth]{./fig/enu_by_month_narrow.eps}
\caption{\label{fig:enu_by_month} Reconstructed neutrino energy distributions by calendar month. The distributions are normalized to protons on target and only data obtained in the LE10/185kA configuration is included. Markers representing individual months are offset to clarify the presentation and the line represents the average over all months. }
\end{figure*}

\subsubsection{Far Detector Data}
\label{sec:fd_dist}

Several criteria were applied to assure the integrity of the Far Detector data. These explicitly exclude data taken during magnet coil and high voltage power failures and periods where the GPS timing system was not operational. The POT-weighted live-time of the Far Detector was 98.9\% during the period May 20, 2005 -- Feb 25, 2006 in which the oscillation dataset was accumulated. Two additional conditions are applied to Far Detector data in order to reduce contamination from cosmic ray events to a negligible level. First, the direction cosine of the primary reconstructed track with respect to the beam direction must be greater than 0.6 in order to reject the high angle tracks typical of cosmic rays. Furthermore, the event time must lie within a $\unit[50]{\mu s}$ window around the predicted time of the beam spill at the Far Detector site.

Figure~\ref{fig:fdtiming} shows the time (relative to the spill time) of the 384 selected neutrino events that satisfy the criteria described above. The width of this distribution is consistent with the Main Injector spill length, and there is no evidence of background contamination within the $\unit[50]{\mu s}$ window.  Figure~\ref{fig:fd_pot} shows the cumulative distributions of the number of total protons on target and the number of Far Detector neutrino events as a function of time. The two distributions follow each other closely. Figure~\ref{fig:nu_per_pot} shows the number of neutrinos per $10^{17}$ protons on target and exhibits no significant time dependence.

\begin{figure}
\begin{center}
  \includegraphics[width=\columnwidth]{fig/fd_timing.eps}
  \caption{Distribution of the event time of the 384 Far Detector neutrino events relative to the time of the nearest beam spill~\cite{tofpaper}.}
   \label{fig:fdtiming}
\end{center}
\end{figure}

\begin{figure}[h]
\centering
\subfloat[\label{fig:fd_pot}]{\includegraphics[width=0.7\textwidth]{./fig/fdpot.eps} }
\\
\subfloat[\label{fig:nu_per_pot}]{\includegraphics[width=0.7\textwidth]{./fig/neus_per_pot_prd.eps}}
 \caption{\label{fig:ccpot} \protect\subref{fig:fd_pot} Cumulative distributions of the total number of protons on target and the number of Far Detector neutrino events observed as a function of time. \protect\subref{fig:nu_per_pot} Number of observed Far Detector neutrino interactions per $10^{17}$ protons on target as a function of time.}
\end{figure}

A visual scan of the events was carried out in order to ensure that there was no background contamination in the selected event sample and that the events were well reconstructed. Figure~\ref{fig:fd_vtx} shows the reconstructed vertex distributions of the Far Detector events.  The Monte Carlo distributions are normalized to the same number of events as the data to account for the possible effects of oscillations and are in agreement with the data. Additional comparisons between the data and simulation are shown in the sections that follow.

\begin{figure}
\centering
\subfloat[]{\includegraphics[width=0.4\textwidth,clip]{fig/prd_fd_vertices_1.eps}}
\\
\subfloat[]{\includegraphics[width=0.4\textwidth,clip]{fig/prd_fd_vertices_2.eps}}
\\
\subfloat[]{\includegraphics[width=0.4\textwidth,clip]{fig/prd_fd_vertices_3.eps}}
\caption{Distributions of track vertices for Far Detector data and Monte Carlo simulation. The Monte Carlo distributions are normalized to the same number of events as the data. \label{fig:fd_vtx}}
\end{figure}
\clearpage

redpajama_set_name: RedPajamaArXiv