\documentstyle[12pt,epsfig]{article}

\setlength{\textwidth}{6.5in}
\setlength{\textheight}{8.5in}
%\setlength{\footskip}{0.35in}
\setlength{\topmargin}{0.0in}
\setlength{\oddsidemargin}{0in}
\setlength{\parindent}{5ex}
\setlength{\parskip}{.12in}
\newcommand{\detw}{\mbox{$\frac{\Delta E}{2}$}}
\newcommand{\hs}{\hspace*{5ex}}


%\def\decayright#1{\kern#1em\raise1.1ex\hbox{$|$}\kern-.5em\rightarrow}

\begin{document}

\begin{center}
{\Large
Summary of a Simulation of the ANITA Detector
}
~\\
Amy Connolly\\
\today\\

~\\
~\\
{\bf Abstract}\\

\parbox[t]
{5.5in}
{
\hs We outline our code written to simulate the ANITA
detection system.
}
\end{center}

\tableofcontents


\section{Introduction}

The code outlined here is meant to be a simulation of the
ANITA and ANITA-lite detection systems.  It should provide
the collaboration with a tool to study the impact of various
effects on the experiment's sensitivity.  It will also 
provide the energy-dependent sensitivity that will be used
for setting limits on models for UHE neutrinos from
ANITA-lite data and ultimately from the full ANITA flights.
It also serves as a consistency check with the Monte Carlo
simulation being developed by our collaborators at the
University of Hawaii.

\section{Overall Strategy}
\label{sec:strategy}
\hs The simulation models interactions from neutrinos
of a specified energy.
The Askaryan pulse is parameterized by the current
theoretical model.
We require an interaction to occur within
the Antarctic ice volume within the balloon's horizon.  
We generate two different types of events.  The first is where the
ray seen by the balloon is direct; the ray is emitted from the
interaction upward.  The second is where the ray is emitted downward
and is reflected from the ice-rock interface before being reaching
the surface.  For half of the neutrinos, we force the ray to be the
first type and for the other half, the second type.
For a given type of ray, we find the unique path along which
an RF signal would travel from the interaction to the balloon, snelled
through ice layers 
at the ice-air surface.  Next, we pick a direction for the neutrino
path, only considering directions such that the Cerenkov cone is
close enough to the unique ray from interaction to balloon that
the signal is still detectable under a best-case scenario.
Depth-dependent attenuation lengths and indices of refraction
in the ice are based on recent South Pole measurements.  Surface
slopeyness is taken into account under a simple model.
Frequency-dependent 
antenna response is based on the manufacturer's specifications.
The trigger configuration is based on the most
current design.
The restrictions in neutrino phase space 
that come about in our selection of neutrino interaction positions
and neutrino directions 
are corrected for in the
final calculation (see Section~\ref{sec:volstr}).


\subsection{Coordinate System}


\hs Where Cartesian coordinates are used, we define $x$,$y$ and $z$ to 
be zero at the center of the earth.  The $+z$ axis 
runs from the center of the earth
to the south pole (so, our coordinate system is defined so that
the earth is ``upside-down'').  The $x$ and $y$ axes lie in the
0$^{\circ}$ latitude plane, with $+x$ pointing to
the 90$^{\circ}$ E longitudinal line, and $+y$ points to
0$^{\circ}$ longitude.  This coordinate system is right-handed.

The $\phi$ coordinate is defined as usual, zero along the $+x$
axis, and increasing moving counter-clockwise from the view 
looking down on the
$+z$ axis.  The $\theta$ coordinate is measured relative to the
$+z$ axis.

Bedrock and ice surface elevations are measured radially from
the center of the earth, and are quoted relative to sea level.
Sea level $r(\theta)$ is latitude-dependent and is quoted as a
distance from the center of the earth.  
We use a geoid shape given by~\cite{geoid}:
\begin{equation}
r(\theta)=\frac{r_{min}\cdot r_{max}}{\sqrt{r_{min}^2-(r_{min}^2-r_{max}^2) \cos^2{\theta}}}  
\end{equation}
where $r_{min}=6356.752$ km and $r_{max}=6378.137$ km.

To accommodate the global crust model 2.0 \cite{crust2}
for thicknesses of sediment layers through the crust as well as ice
thicknesses, which is on a $2^{\circ} \times 2^{\circ}$ grid, we bin
the earth's  surface into 180 bins in longitude and
90 bins in latitude.  Note that the resulting bins are not uniform
in area.
 The bins in latitude are
numbered such that the 0$^{th}$ bin is closest to the $+z$ axis (at the
south pole), and the 89$^{th}$ bin is closest to the north pole ($-z$).


\subsection{CSEDI Crust 2.0}
\label{sec:crust2}
\hs Crust 2.0 is the latest model of the Earth's interior near the 
surface published as the result of
an initiative called Cooperative Studies of the Earth's Deep
Interior (CSEDI).  It is based on seismological data, and
the model gives thicknesses and densities of each of seven
layers in $2^{\circ} \times 2^{\circ}$ bins:  ice, water, soft sediments, hard sediments, upper crust,
middle crust, and lower crust.  The ice thicknesses in Crust 2.0 are 
claimed to be within 250 m of the true ice thickness. Sediment thicknesses 
in each cell are to within 1.0 km of the true sediment thickness 
and crustal thickness are within 5 km of the true crustal thicknesses. 




%\hs  The user may input a list of
%bedrock and ice surface altitudes along with their coordinates
%in longitude and latitude, and the gaps are filled in with a 
%smoothing process.  In other words, for the bins where no data
%is given,  the altitudes for the bedrock and ice surface are
%interpolated.  The data that we are using now for bedrock and 
%ice surface elevations comes from: \\
%\centerline{http://classroomantarctica.aad.gov.au/materials/Elevations.pdf }
%which gives rockbed and ice surface elevations only along a single
%line across the continent. 

%The smoothing procedure is as follows.  The longitude and
%latitude coordinates for the inputted data points are translated
%to their appropriate $\phi$ and $\cos{\theta}$ bins and the
%bedrock and ice surface altitudes are set for these bins according
%to the data.  Also, the bins along the edge of the circle that
%defines our coordinate system (latitude=-60$^{\circ}$ and all bins
%in $\phi$) are set to be at sea level (altitudes set to zero
%relative to the earth radius).  All of the bins that have been
%set up to this point are called \emph{reference bins}.

%For a given bin that is not a reference bin, we need to
%set its bedrock and ice surface altitude.  First, we find
%the closest reference bin that is non-trivial.  A trivial reference
%bin is one set to sea level on the outer edge of the circle.
%Next, we find the closest reference bin (this time allowing trivial
%reference bins) that is on the other side of the bin of interest.
%So we want to be able to draw a straight line whose endpoints are
%reference points (at least one non-trivial) and that includes the
%bin of interest somewhere along the line.  Then, the bedrock altitude
%for the bin of interest is set to a value between those of the
%two reference points, determined by the bin's relative distance 
%from each of the reference points.  The same is done for the ice
%surface altitude.

Figures~\ref{fig:upperlayers} and~\ref{fig:lowerlayers} show the 
elevations of the seven crustal layers included in Crust 2.0.
\begin{figure}
\begin{center}
\epsfig{figure=upperlayers.eps,height=3in}
\caption{Altitude of the upper four layers given in Crust 2.0 along the
the $75^{\circ}$ S latitude line.  The horizontal axis is degrees in
longitude.
\label{fig:upperlayers}}
\end{center}
\end{figure}
\begin{figure}
\begin{center}
\epsfig{figure=lowerlayers.eps,height=3in}
\caption{Altitude of the lower three layers given in Crust 2.0 along the
the $75^{\circ}$ S latitude line.
\label{fig:lowerlayers}}
\end{center}
\end{figure}

We compute
the total Antarctic ice volume by summing the product of ice
thickness and surface area for each bin within the Antarctic continent.
For the area of each bin, we use:
\begin{equation}
\int_{\phi_1}^{\phi_2}\int_{\theta_1}^{\theta_2}\sin{\theta}~d\theta~d\phi=\left(\phi_1-\phi_2\right)\times\left(\cos{\theta_1}-\cos{\theta_2}\right)
\end{equation}
Where the limits of the integrals define the edges of the bin in
latitude and longitude.
We find $2.7 \times 10^{16}$ m$^2$ of Antarctic ice in this model.
Compare this to the $3.01098 \times 10^{16}$ m$^3$ volume of ice
in Antarctica reported by the US Geological Survey~\cite{usgs};
they are different by 10\%.
Unfortunately, Crust 2.0 calls any ice that sits above water (ice
shelves) just water.  Therefore, the Ross Ice Shelf needs to be
added into the model by hand.  Ice shelves make up 2.4\% of the
ice volume in Antarctica, so it won't account for the entire
$10\%$ difference. 

%\begin{center}
%\begin{tabular}{|l|l|l|l|}
%\hline\hline
%term      &   symbol &  with arguments & units for neutrinos\\
%\hline\hline
%neutrinos &
%  $N_{\nu}$    & 
%  $N_{\nu}$ &
%  (unitless) \\
%\hline
%neutrino intensity & 
%  $I$ &
%  $I(\theta,\phi)$  &
%  [length]$^{-2}$ [sr]$^{-1}$ [time]$^{-1}$\\
%\hline
%neutrino brightness, {\it or} & 
%  $I_{E}$, &
%  $I_{E}(E, \theta, \phi)$ &
%  [length]$^{-2}$ [sr]$^{-1}$ [time]$^{-1}$ [energy]$^{-1}$\\
%neutrino specific intensity  &
%  &
%  &
%  \\
%\hline
%net flux & 
%  $F_{E}$ &
%  $F_{E} (E)$ &
%  [length]$^{-2}$ [time]$^{-1}$ [energy]$^{-1}$\\
%\hline 
%total integrated flux &
%  $F$ &
%  $F$ &
%  [length]$^{-2}$ [time]$^{-1}$\\
%\hline\hline
%\end{tabular}
%\end{center}


\section{The Askaryan Signal}

\subsection{Magnitude of the Pulse}

We use the parameterization outlined in~\cite{jaime} for the peak of the
Askaryan signal:
\begin{equation}
\label{eq:vmmhz}
\frac{\cal E^{(\mathrm{@ 1m})} }{\mathrm{V/m/MHz}}=2.53\times 10^{-7}\cdot \frac{\sin{\theta_v}}{\sin{\theta_c}}\cdot  \frac{E_{em}}{\mathrm{TeV}} \cdot \frac{\nu}{\nu_0} \cdot \frac{1}{1.+\left( \frac{\nu}{\nu_0} \right)^{1.44}}
\end{equation}
where $\nu$ is the frequency, $\nu_0=1.15$ GHz, $E_{em}$ is the shower energy, $\theta_v$ is the viewing angle and $\theta_c$ is the Cerenkov angle. 

\subsection{Electromagnetic and Hadronic Components of the Shower}
 We assume flavor democracy; that by the time the neutrinos get
here, the flavors are fully mixed.  We also keep track of each
flavor individually so the sensitivity to each flavor may be
quoted separately.

We pick an inelasticity $y$ according to~\cite{ghandi},
which we have approximated by a double-exponential.
If we call the electromagnetic component $f_e$ and the 
hadronic component $f_h$, then the dependence of $f_e$ and
$f_h$ on $y$ depends on the whether it is a charged current
interaction that occurred or a neutral current interaction.
We choose
$70.64\%$ of the events at random to be charged-current
and the remainder to be neutral current.
If the incident neutrino was a $\nu_e$ and it is a charged current
event, then $f_e=1-y$ and $f_h=y$.
If the incident neutrino was $\nu_\tau$ or $\nu_\mu$, or
if the $\nu_e$ interacts through a neutral current, then
the electromagnetic component is for now treated as negligible, and
$f_h=y$.  A future version of the code will include $\mu/\tau$
bremsstrahlung and photo-nuclear interactions.


\subsection{Width of Cerenkov Cone}
\label{sec:cerenkov_width}

The user can select one of two different parameterizations of the
width of the Cerenkov cone in inputs.txt.  Here we refer to the
two parameterizations as the ``old'' parameterization and 
``new'' parameterization, for lack of a better terminology.
The user should select 
0 (1) for the old (new) parameterizations in the appropriate line in