\documentclass{book}
\usepackage{amsmath} % for \text command
\usepackage{fullpage}
\usepackage{tikz}
\usepackage{hyperref}
\usepackage{amsfonts}
\usepackage{algorithmic}
\renewcommand{\algorithmiccomment}[1]{\# #1}
\usepackage{algorithm}
\newcommand*\diff{\mathrm{d}}
\usetikzlibrary{shapes}
\author{Anthony LaTorre}
\date{\today}
\title{Searching for Dark Matter with the Sudbury Neutrino Observatory}
\begin{document}
\maketitle
\chapter{Introduction}
\chapter{Estimating the Event rate in the SNO detector}
The event rate of self destructing dark matter events, $R$, in the SNO detector is given by first integrating over the detector.
\begin{equation}
R = \int_\mathrm{SNO} \mathrm{d}^3r \, R(r) 
\end{equation}
Next, we integrate over the earth where the dark matter annihilates:
\begin{equation}
R = \int_\mathrm{SNO} \mathrm{d}^3r \, \int_{r'} \mathrm{d}^3r' R(r') \mathrm{P}(\text{detect at r} | \text{DM scatters at r'})
\end{equation}
where we have assumed above that the dark matter annihilates immediately after
scattering. The event rate for scattering at a position $r'$ in the earth is:
\begin{equation}
R(r') = \Phi(r') \eta(r') \sigma(r')
\end{equation}
where $\Phi(r')$ is the flux at $r'$, $\eta(r')$ is the number density of
scatterers at $r'$, and $\sigma$ is the cross section for the dark matter to
scatter and annihilate. In general this will be a sum over the elemental
composition of the earth at $r'$, but for notational simplicity we will assume
a single cross section. We will also assume that the cross section is small
enough that the flux is essentially constant over the whole earth so that the
rate may be written as:
\begin{equation}
R(r') = \Phi \eta(r') \sigma(r').
\end{equation}
The rate may then be written as:
\begin{equation}
R = \Phi \int_\mathrm{SNO} \mathrm{d}^3r \, \int_{r'} \mathrm{d}^3r' \, \eta(r') \sigma(r') \mathrm{P}(\text{detect at r} | \text{DM scatters at r'})
\end{equation}
If we assume that the probability of detecting the dark matter is uniform
across the SNO detector we may write it as:
\begin{equation}
R = \Phi \int_\mathrm{SNO} \mathrm{d}^3r \, \int_{r'} \mathrm{d}^3r' \, \eta(r') \sigma(r') \mathrm{P}(\text{detect at SNO} | \text{DM scatters at r'})
\end{equation}
This assumption is pretty well motivated since for most values of the mediator
decay length the probability will be uniform across the detector. The only
value for which it might not be a good approximation is if the mediator decay
length is on the order of the detector radius in which case DM scattering in
the rock of the cavity walls might have a higher event rate at the edge of the
detector. Since the integrand no longer depends on $r$, we may write it as:
\begin{equation}
R = \Phi V_\text{fiducial} \int_{r'} \mathrm{d}^3r' \, \eta(r') \sigma(r') \mathrm{P}(\text{detect at SNO} | \text{DM scatters at r'})
\end{equation}
where $V_\text{fiducial}$ is the fiducial volume of the detector. The
probability that the mediator $V$ is emitted in a direction $\theta$ and
travels a distance $r$ in a spherical coordinate system centered on $r'$ may be
written as:
\begin{equation}
f(r,\theta) = \frac{\sin\theta}{4\pi}\frac{e^{-r/L_V}}{L_V}
\end{equation}
To transform this probability distribution to the coordinate system centered on the SNO detector we first transform it to a cartesian coordinate system:
\begin{align}
f(r,\theta) &= \frac{\sin\theta}{4\pi}\frac{e^{-r/L_V}}{L_V}\frac{1}{r^2sin\theta} \\
&= \frac{1}{4\pi r^2}\frac{e^{-r/L_V}}{L_V}
\end{align}
Then, the distribution is translated to the center of the detector, which
doesn't change the form since the radial coordinate $r$ is the same in both
coordinate systems. Finally, we switch back into spherical coordinates:
\begin{align}
f(r,\theta') &= \frac{1}{4\pi r^2}\frac{e^{-r/L_V}}{L_V}r^2\sin\theta' \\
&= \frac{1}{4\pi}\frac{e^{-r/L_V}}{L_V}\sin\theta'
\end{align}
where $\theta'$ is the polar angle in the SNO coordinate system. We can now write the rate as:
\begin{align}
R &= \Phi V_\text{fiducial} \int_r \mathrm{d}r \, \int_\theta \mathrm{d}\theta \, \int_\phi \mathrm{d}\phi \, \eta(r,\theta,\phi) \sigma(r,\theta,\phi) \frac{1}{4\pi}\frac{e^{-r/L_V}}{L_V}\sin\theta
\end{align}
We now assume that the number density of scatterers $\eta$ and the cross
section $\sigma$ are independent of the position in the earth. This is a good
approximation for certain values of $L_V$ since the integral will be dominated
by a single material. For example, if the mediator decay length $L_V$ is
approximately 1 meter, then the vast majority of the events in the detector
will be caused by dark matter scattering off of water. Similarly if the
mediator decay length is approximately 1 km then the majority of the events in
the detector will be caused by the dark matter scattering off of the norite
rock surrounding the detector. With this approximation, the rate may be
written:
\begin{align}
R &= \Phi V_\text{fiducial} \eta \sigma \int_r \mathrm{d}r \, \int_\theta \mathrm{d}\theta \, \int_\phi \mathrm{d}\phi \, \frac{1}{4\pi}\frac{e^{-r/L_V}}{L_V}\sin\theta \\
&= \Phi V_\text{fiducial} \eta \sigma \int_r \mathrm{d}r \, \frac{1}{2}\frac{e^{-r/L_V}}{L_V} \int_\theta \mathrm{d}\theta \, \sin\theta
\end{align}
The $\theta$ integral goes from $\theta_\text{min}$ to $\pi$:
\begin{align}
R &= \Phi V_\text{fiducial} \eta \sigma \int_r \mathrm{d}r \, \frac{1}{2}\frac{e^{-r/L_V}}{L_V} \int_{\theta_\text{min}}^\pi \mathrm{d}\theta \, \sin\theta
\end{align}
where $\theta_\text{min}$ is equal to:
\begin{equation}
\theta_\text{min} =%
\begin{cases}
0 & \text{if } r < \text{depth} \\
\pi - \arccos\left(\frac{\text{depth}^2 + r^2 - 2R\text{depth}}{2r(R-\text{depth})}\right) & \text{if } \text{depth} < r < 2R-\text{depth} \\
\end{cases}
\end{equation}
where $R$ is the radius of the earth and $\text{depth}$ is the distance from
the surface of the earth to the SNO detector.

\chapter{Cross Section}
In \cite{grossman2017} the differential scattering cross section for dark
matter off a nucleus is calculated as
\begin{equation}
\frac{\diff \sigma_\text{scatter}}{\diff q^2} = \frac{g_V^2 \epsilon^2 e^2}{4\pi v^2 (q^2 + m_V^2)^2} |F_D(q^2)|^2 Z^2 F^2(q),
\end{equation}
where $q$ is the momentum transferred, $g_V$ and $\epsilon$ are coupling
constants (FIXME: is this true?), $v$ is the velocity of the dark matter
particle, $m_V$ is the mass of the mediator, $F_D(q^2)$ is a form factor for
the dark matter to transition from a high angular momentum state to a lower
angular momentum state, $Z$ is the atomic number of the nucleus, and $F^2(q)$
is a nuclear form factor.

In the limit of low momentum transfer, the cross section is approximately
\begin{equation}
\frac{\diff \sigma_\text{scatter}}{\diff q^2} \simeq \frac{g_V^2 \epsilon^2 e^2}{4\pi v^2 m_V^4} |F_D(q^2)|^2 Z^2 F^2(q).
\end{equation}

For existing direct detection dark matter experiments, the relevant cross
section is (FIXME: is this true?)
\begin{equation}
\frac{\diff \sigma_\text{scatter}}{\diff q^2} \simeq \frac{g_V^2 \epsilon^2 e^2}{4\pi v^2 m_V^4} Z^2 F^2(q).
\end{equation}

A standard cross section can be defined as the total cross section in the zero
momentum limit\cite{pepin2016}
\begin{align}
\sigma_0 &= \int_0^{4\mu_T^2 v^2} \frac{\diff \sigma_\text{scatter}}{\diff q^2}\bigg\rvert_{q \rightarrow 0} \diff q^2 \\
&= \frac{\mu_T^2 g_V^2 \epsilon^2 e^2}{\pi m_V^4} Z^2,
\end{align}
where $\mu_T$ is the reduced mass of the WIMP and target nucleus.

Since different experiments use different detector targets, it is also useful
to define a standard cross section, $\sigma_p$ which is independent of the
nuclear target:
\begin{equation}
\sigma_p = \left(\frac{\mu_p}{\mu_T}\frac{1}{Z}\right)^2 \sigma_0.
\end{equation}

The direct detection cross section is then:
\begin{equation}
\frac{\diff \sigma_\text{scatter}}{\diff q^2} \simeq \frac{1}{4 \mu_p^2 v^2} \sigma_p Z^2 F^2(q).
\end{equation}
and the cross section for the dark matter to annihilate is:
\begin{equation}
\frac{\diff \sigma_\text{scatter}}{\diff q^2} \simeq \frac{1}{4 \mu_p^2 v^2} \sigma_p |F_D(q)|^2 Z^2 F^2(q).
\end{equation}

\section{Nuclear Form Factor}
The nuclear form factor, $F(q)$, characterizes the loss of coherence as the de
Broglie wavelength of the WIMP approaches the radius of the
nucleus\cite{caldwell2015}. The most commonly used form factor calculation used
in the direct detection community is that of Helm which is given by:
\begin{equation}
F(q) = 3\frac{j_1(q r_1)}{q r_1} e^{-\frac{(q s)^2}{2}},
\end{equation}
where $j_1$ is the spherical bessel function of the first order, $s$ is a
measure of the nuclear skin thickness, and $r_1$ is a measure of the nuclear
radius\cite{pepin2016}. The values used for these constants were
\begin{align}
s &= 0.9 \text{ fm} \\
a &= 0.52 \text{ fm} \\
c &= 1.23 A^\frac{1}{3} - 0.60 \text{ fm} \\
r_1 &= \sqrt{c^2 + \frac{7}{3}\pi^2 a^2 - 5 s^2}
\end{align}

\begin{figure}
\centering
\begin{tikzpicture}[scale=0.1]
% earth
\draw [thick,domain=120:150] plot[smooth] ({200*cos(\x)},{200*sin(\x)});
\begin{scope}[shift={(-100,100)},rotate=45]
% interaction
\node[star,star points=9,draw] at (20,20){};
% acrylic vessel
\draw [thick,domain={90+asin(75/600)}:{360+90-asin(75/600)}] plot[smooth] ({6*cos(\x)},{6*sin(\x)});
\draw [thick] ({6*cos(90+asin(75/600))},{6*sin(90+asin(75/600))}) -- ({6*cos(90+asin(75/600))},{6*sin(90+asin(75/600))+7.5}) --
      ({6*cos(90+asin(75/600))+2*0.75},{6*sin(90+asin(75/600))+7.5}) --
      ({6*cos(90+asin(75/600))+2*0.75},{6*sin(90+asin(75/600))});
% PSUP
\draw [domain=0:360] plot ({8.89*cos(\x)},{8.89*sin(\x)});
% cavity
\draw (-9.5,-10.5) --
      (-10.6,-10.5+5.6) --
      (-10.6,-10.5+14.93) --
      (-9.5,-10.5+30) --
      (9.5,-10.5+30) --
      (10.6,-10.5+14.93) --
      (10.6,-10.5+5.6) --
      (9.5,-10.5) --
      (-9.5,-10.5);
\draw[->,ultra thick] (-25,0) -- (25,0) node[right]{$x$};
\draw[->,ultra thick] (0,-25) -- (0,25) node[right]{$y$};
\end{scope}
\end{tikzpicture}
\end{figure}

\chapter{Event Reconstruction}
In order to reconstruct the physical parameters associated with an event we
compute a likelihood for that event given a proposed energy, position,
direction, and initial time. The likelihood may be written as:
\begin{equation}
\label{likelihood}
\mathcal{L}(E, \vec{x}, \vec{v}, t_0) = P(\vec{q}, \vec{t} | E, \vec{x}, \vec{v}, t_0)
\end{equation}
where $E$, $\vec{x}$, $\vec{v}$ represent the initial particle's kinetic
energy, position, and direction respectively, $t_0$ represents the initial time
of the event, $\vec{q}$ is the charge seen by each PMT, and $\vec{t}$ is the
time recorded by each PMT.

In general the right hand side of Equation~\ref{likelihood} is not factorable
since for particle tracks which scatter there will be correlations between the
PMT hits. However, to make the problem analytically tractable, we assume that
the probability of each PMT being hit is approximately independent of the
others. With this assumption we can factor the right hand side of the
likelihood as:
\begin{equation}
\mathcal{L}(E, \vec{x}, \vec{v}, t_0) = \prod_i P(\text{not hit} | E, \vec{x}, \vec{v}, t_0) \prod_j P(\text{hit}, q_j, t_j | E, \vec{x}, \vec{v}, t_0)
\end{equation}
where the first product is over all PMTs which weren't hit and the second
product is over all of the hit PMTs.

If we introduce the variable $n$ which represents the number of photoelectrons detected we can write the likelihood as:
\begin{equation}
\mathcal{L}(E, \vec{x}, \vec{v}, t_0) = \prod_i P(n = 0 | E, \vec{x}, \vec{v}, t_0) \prod_j \sum_{n = 1}^{\infty} P(n, q_j, t_j | E, \vec{x}, \vec{v}, t_0)
\end{equation}

We can factor the right hand side of the likelihood as:
\begin{equation}
\mathcal{L}(E, \vec{x}, \vec{v}, t_0) = \prod_i P(n = 0 | E, \vec{x}, \vec{v}, t_0) \prod_j \sum_{n = 1}^{\infty} P(q_j, t_j | n, E, \vec{x}, \vec{v}, t_0) P(n | E, \vec{x}, \vec{v}, t_0)
\end{equation}

If we now assume that the charge and time observed at a given PMT are
independent we can write the likelihood as:
\begin{equation}
\mathcal{L}(E, \vec{x}, \vec{v}, t_0) = \prod_i P(n = 0 | E, \vec{x}, \vec{v}, t_0) \prod_j \sum_{n = 1}^{\infty} P(q_j | n, E, \vec{x}, \vec{v}, t_0) P(t_j | n, E, \vec{x}, \vec{v}, t_0) P(n | E, \vec{x}, \vec{v}, t_0)
\end{equation}

Since there are many photons produced in each event each of which has a small
probability to hit a given PMT, we will assume that the probability of
detecting $n$ photons at a given PMT is poisson distributed, i.e.
\begin{equation}
P(n | E, \vec{x}, \vec{v}, t_0) = e^{-\mu} \frac{\mu^n}{n!}
\end{equation}

We can therefore write the likelihood as:
\begin{equation}
\mathcal{L}(E, \vec{x}, \vec{v}, t_0) = \prod_i e^{-\mu_i} \prod_j \sum_{n = 1}^{\infty} P(q_j | n, E, \vec{x}, \vec{v}, t_0) P(t_j | n, E, \vec{x}, \vec{v}, t_0) e^{-\mu_j} \frac{\mu_j^n}{n!}
\end{equation}
where $\mu_i$ is the expected number of photoelectrons detected at the ith PMT
(given an initial particle's energy, position, and direction).

First, we'll calculate the expected number of photoelectrons for a single non-showering track which undergoes multiple scattering through small angles. In this case, we can calculate the expected number of photoelectrons as:
\begin{equation}
\mu_i = \int_x \diff x \int_\lambda \diff \lambda \frac{\diff^2 N}{\diff x \diff \lambda} P(\text{detected} | E, x, v)
\end{equation}
where $x$ is the position along the track and $\lambda$ is the wavelength of
the light.

If the particle undergoes many small angle Coulomb scatters, the net
angular displacement of the particle after a distance $x$ will be a Gaussian
distribution by the central limit theorem\cite{pdg2017}. The distribution of
the net angular displacement at a distance $x$ along the track is then given
by\footnote{This distribution will be correlated between different points along the track.}:
\begin{equation}
f(\theta,\phi) = \frac{\theta}{2\pi\theta_0^2}e^{-\frac{\theta^2}{2\theta_0^2}}
\end{equation}
where
\begin{equation}
\theta_0 = \frac{13.6 \text{ MeV}}{\beta c p}z\sqrt{\frac{x}{X_0}}\left[1 + 0.038\ln\left(\frac{x z^2}{X_0 \beta^2}\right)\right]
\end{equation}
where $p$, $\beta c$, and $z$ are the momentum, velocity, and charge of the
particle, and $X_0$ is the radiation length of the particle\cite{pdg2017}.

Now, we integrate over the angular displacement of the track around the original velocity:
\begin{align}
\mu_i &= \int_x \diff x \int_\lambda \diff \lambda \frac{\diff^2 N}{\diff x \diff \lambda} \int_\theta \diff \theta \int_\phi \diff \phi P(\text{detected} | \theta, \phi, E, x, v) P(\theta, \phi | E, x, v) \\
\mu_i &= \int_x \diff x \int_\lambda \diff \lambda \frac{\diff^2 N}{\diff x \diff \lambda} \int_\theta \diff \theta \int_\phi \diff \phi P(\text{detected} | \theta, \phi, E, x, v) f(\theta,\phi)
\end{align}
The probability of being detected can be factored into several different compontents:
\begin{align}
\mu_i &= \int_x \diff x \int_\lambda \diff \lambda \frac{\diff^2 N}{\diff x \diff \lambda} \int_\theta \diff \theta \int_\phi \diff \phi P(\text{emitted towards PMT i} | \theta, \phi, E, x, v) f(\theta,\phi) P(\text{not scattered or absorbed} | \lambda, E, x, v) \epsilon(\eta) \mathrm{QE}(\lambda) \\
\label{eq:mui}
\mu_i &= \int_x \diff x \int_\lambda \diff \lambda \frac{\diff^2 N}{\diff x \diff \lambda} P(\text{not scattered or absorbed} | \lambda, E, x, v) \epsilon(\eta) \mathrm{QE}(\lambda) \int_\theta \diff \theta \int_\phi \diff \phi P(\text{emitted towards PMT i} | \theta, \phi, E, x, v) f(\theta,\phi)
\end{align}
where $\eta$ is the angle between the vector connecting the track position $x$
to the PMT position and the normal vector to the PMT, $\epsilon(\eta)$ is the
collection efficiency, and $\mathrm{QE}(\lambda)$ is the quantum efficiency of
the PMT.

The probability that a photon is emitted directly towards a PMT is given by a
delta function (we make the assumption here that the probability is uniform
across the face of the PMT):
\begin{equation}
P(\text{emitted towards PMT i} | \theta, \phi, E, x, v) = \delta\left(\frac{1}{n(\lambda)\beta} - \cos\theta'(\theta,\phi,x)\right) \frac{\Omega(x)}{4\pi}
\end{equation}
where $\theta'$ is the angle between the track and the PMT and $\Omega(x)$ is the solid angle subtended by the PMT.

In a coordinate system with the z axis aligned along the original particle velocity and with the PMT in the x-z plane, the angle $\theta'$ is defined by:
\begin{equation}
\cos\theta' = \sin\theta\cos\phi\sin\theta_1 + \cos\theta\cos\theta_1
\end{equation}
where $\theta_1$ is the angle between the PMT and the original particle velocity.

We can now solve the integral on the right hand side of Equation~\ref{eq:mui} as:
\begin{align}
P(\text{emitted towards PMT i}) &= \int_\theta \diff \theta \int_\phi \diff \phi \delta\left(\frac{1}{n(\lambda)\beta} - \cos\theta'(\theta,\phi,x)\right) \frac{\Omega(x)}{4\pi} \frac{\theta}{2\pi\theta_0^2}e^{-\frac{\theta^2}{2\theta_0^2}} \\
P(\text{emitted towards PMT i}) &= \frac{\Omega(x)}{4\pi} \frac{1}{2\pi\theta_0^2}\int_\theta \diff \theta \int_\phi \diff \phi \delta\left(\frac{1}{n(\lambda)\beta} - \cos\theta'(\theta,\phi,x)\right) \theta e^{-\frac{\theta^2}{2\theta_0^2}} \\
P(\text{emitted towards PMT i}) &= \frac{\Omega(x)}{4\pi} \frac{1}{2\pi\theta_0^2}\int_\theta \diff \theta \int_\phi \diff \phi \delta\left(\frac{1}{n(\lambda)\beta} - \sin\theta\cos\phi\sin\theta_1 - \cos\theta\cos\theta_1\right) \theta e^{-\frac{\theta^2}{2\theta_0^2}}
\end{align}

We now assume $\theta$ is small (which should be valid for small angle scatters), so that we can rewrite the delta function as:
\begin{align}
P(\text{emitted towards PMT i}) &= \frac{\Omega(x)}{4\pi} \frac{1}{2\pi\theta_0^2}\int_\theta \diff \theta \int_\phi \diff \phi \delta\left(\frac{1}{n(\lambda)\beta} - \theta\cos\phi\sin\theta_1 - \cos\theta_1\right) \theta e^{-\frac{\theta^2}{2\theta_0^2}}
\end{align}

We can rewrite the delta function and solve the integral as:
\begin{align}
P(\text{emitted towards PMT i}) &= \frac{\Omega(x)}{4\pi} \frac{1}{2\pi\theta_0^2}\int_\theta \diff \theta \int_\phi \diff \phi \frac{1}{\left|\cos\phi\sin\theta_1\right|}\delta\left(\theta - \frac{\frac{1}{n(\lambda)\beta}-\cos\theta_1}{\cos\phi\sin\theta_1}\right) \theta e^{-\frac{\theta^2}{2\theta_0^2}} \\
&= \frac{\Omega(x)}{4\pi} \frac{1}{2\pi\theta_0^2} \frac{1}{\left|\sin\theta_1\right|} \int_\phi \diff \phi \frac{1}{\left|\cos\phi\right|} \int_\theta \diff \theta \delta\left(\theta - \frac{\frac{1}{n(\lambda)\beta}-\cos\theta_1}{\cos\phi\sin\theta_1}\right) \theta e^{-\frac{\theta^2}{2\theta_0^2}} \\
&= \frac{\Omega(x)}{4\pi} \frac{1}{2\pi\theta_0^2} \frac{1}{\left|\sin\theta_1\right|} \int_\phi \diff \phi \frac{1}{\left|\cos\phi\right|}\frac{\frac{1}{n(\lambda)\beta}-\cos\theta_1}{\cos\phi\sin\theta_1}H\left(\frac{\frac{1}{n(\lambda)\beta}-\cos\theta_1}{\cos\phi\sin\theta_1}\right)e^{-\frac{1}{2\theta_0^2}\left(\frac{\frac{1}{n(\lambda)\beta}-\cos\theta_1}{\cos\phi\sin\theta_1}\right)^2} \\
&= \frac{\Omega(x)}{4\pi} \frac{1}{2\pi\theta_0^2} \frac{1}{\left|\sin\theta_1\right|}\sqrt{2\pi}\theta_0 e^{-\frac{1}{2\theta_0^2}\left(\frac{\frac{1}{n(\lambda)\beta}-\cos\theta_1}{\sin\theta_1}\right)^2} \\
&= \frac{\Omega(x)}{4\pi} \frac{1}{\sqrt{2\pi}\theta_0} \frac{1}{\left|\sin\theta_1\right|} e^{-\frac{1}{2\theta_0^2}\left(\frac{\frac{1}{n(\lambda)\beta}-\cos\theta_1}{\sin\theta_1}\right)^2}.
\end{align}

To simplify this expression we can write
\begin{equation}
P(\text{emitted towards PMT i}) = \frac{\Omega(x)}{4\pi} \frac{1}{\sqrt{2\pi}\theta_0} \frac{1}{\left|\sin\theta_1\right|} e^{-\frac{\Delta^2(\lambda)}{2\theta_0^2}}
\end{equation}
where
\begin{equation}
\Delta(\lambda) = \frac{\frac{1}{n(\lambda)\beta}-\cos\theta_1}{\sin\theta_1}
\end{equation}

Plugging this back into Equation~\ref{eq:mui}
\begin{align}
\label{eq:mui-exact}
\mu_i &= \frac{1}{\sqrt{2\pi}\theta_0} \int_x \diff x \frac{\Omega(x)}{4\pi} \frac{1}{\left|\sin\theta_1\right|} \epsilon(\eta) \int_\lambda \diff \lambda \frac{\diff^2 N}{\diff x \diff \lambda} P(\text{not scattered or absorbed} | \lambda, E, x, v) \mathrm{QE}(\lambda) e^{-\frac{\Delta^2(\lambda)}{2\theta_0^2}}
\end{align}

Ideally we would just evaluate this double integral for each likelihood call,
however the double integral is too computationally expensive to perform for
every likelihood call (FIXME: is this true?). We therefore assume that the
second integral will be dominated by the Bessel function which has a
singularity when it's argument is zero, and rewrite Equation~\ref{eq:mui-exact}
as:
\begin{align}
\mu_i &= 2 \frac{1}{\sqrt{2\pi}\theta_0} \int_x \diff x \Omega(x) \frac{1}{\left|\sin\theta_1\right|} \epsilon(\eta) P(\text{not scattered or absorbed} | \lambda_0, E, x, v) \mathrm{QE}(\lambda_0) e^{-\frac{\Delta^2(\lambda_0)}{4\theta_0^2}} \int_\lambda \diff \lambda \frac{\diff^2 N}{\diff x \diff \lambda} K_0\left(\frac{\Delta^2(\lambda)}{4\theta_0^2}\right)
\end{align}
where $\lambda_0$ is the wavelength at which $\Delta(\lambda) = 0$.

For small values of $\Delta$, the Bessel function may be approximated as:
\begin{equation}
K_0(x) \simeq -\log(x) + \log(2) - \gamma
\end{equation}

We may therefore approximate the expected charge as
\begin{multline}
\label{eq:mui-approx}
\mu_i = 2 \frac{1}{\sqrt{2\pi}\theta_0} \int_x \diff x \Omega(x) \frac{1}{\left|\sin\theta_1\right|} \epsilon(\eta) P(\text{not scattered or absorbed} | \lambda_0, E, x, v) \mathrm{QE}(\lambda_0) e^{-\frac{\Delta^2(\lambda_0)}{4\theta_0^2}} \\
\int_\lambda \diff \lambda \frac{\diff^2 N}{\diff x \diff \lambda} \left(-\log\left(\frac{\Delta^2(\lambda)}{4\theta_0^2}\right) + \log(2) - \gamma\right)
\end{multline}

The number of Cerenkov photons produced per unit length and per unit wavelength
is given by\cite{pdg2017}
\begin{equation}
\frac{\diff^2 N}{\diff x \diff \lambda} = \frac{2\pi\alpha z^2}{\lambda^2}\left(1 - \frac{1}{\beta^2 n^2(\lambda)}\right)
\end{equation}
where $\alpha$ is the fine-structure constant and $z$ is the charge of the
particle in units of the electron charge.

We can therefore write the second integral in Equation~\ref{eq:mui-approx} as
\begin{align}
\int_\lambda \diff \lambda \frac{\diff^2 N}{\diff x \diff \lambda} \left(-\log\left(\frac{\Delta^2(\lambda)}{4\theta_0^2}\right) + \log(2) - \gamma\right) &=
2\pi\alpha z^2 \int_\lambda \diff \lambda \frac{1}{\lambda^2}\left(1 - \frac{1}{\beta^2 n^2(\lambda)}\right) \left(-\log\left(\frac{\Delta^2(\lambda)}{4\theta_0^2}\right) + \log(2) - \gamma\right) \\
\label{eq:lambda}
&\simeq 2\pi\alpha z^2 \left(1 - \frac{1}{\beta^2 n^2(\lambda_0)}\right) \int_\lambda \diff \lambda \frac{1}{\lambda^2}\left(-\log\left(\frac{\Delta^2(\lambda)}{4\theta_0^2}\right) + \log(2) - \gamma\right)
\end{align}

Since the $\Delta$ function only depends on the wavelength through the index
which depends weakly on the wavelength, we can approximate the index of
refraction as:
\begin{equation}
n(\lambda) \simeq a + \frac{b}{\lambda^2}.
\end{equation}

The integral in Equation~\ref{eq:lambda} may then be solved analytically
\begin{multline}
\int_{\lambda_1}^{\lambda_2} \diff \lambda \frac{1}{\lambda^2}\left(-\log\left(\frac{\Delta^2(\lambda)}{4\theta_0^2}\right) + \log(2) - \gamma\right) =
\left[\log(4\theta_0^2) + \log(\sin^2\theta_1) + \log(2) - \gamma\right]\left(\frac{1}{\lambda_1}-\frac{1}{\lambda_2}\right) + \\
\left.\left(-4\sqrt{\frac{a}{b}}\arctan\left(\sqrt{\frac{a}{b}}\lambda\right) +
4\sqrt{\frac{1+a\Delta_0}{b\Delta_0}}\arctan\left(\sqrt{\frac{1+a\Delta_0}{b\Delta_0}}\lambda\right) -
\frac{1}{\lambda}\log\left[\left(\Delta_0+\frac{\lambda^2}{b+a\lambda^2}\right)^2\right]\right)\right|_{\lambda_1}^{\lambda_2}
\end{multline}
where $\lambda_1$ and $\lambda_2$ are chosen to cover the range where the
quantum efficiency is non-zero, typically between 300 nm and 600 nm.

For simplicity we will write this previous expression as $f(x)$
\begin{multline}
f(x) = \left[\log(4\theta_0^2) + \log(\sin^2\theta_1) + \log(2) - \gamma\right]\left(\frac{1}{\lambda_1}-\frac{1}{\lambda_2}\right) + \\
\left.\left(-4\sqrt{\frac{a}{b}}\arctan\left(\sqrt{\frac{a}{b}}\lambda\right) +
4\sqrt{\frac{1+a\Delta_0}{b\Delta_0}}\arctan\left(\sqrt{\frac{1+a\Delta_0}{b\Delta_0}}\lambda\right) -
\frac{1}{\lambda}\log\left[\left(\Delta_0+\frac{\lambda^2}{b+a\lambda^2}\right)^2\right]\right)\right|_{\lambda_1}^{\lambda_2}
\end{multline}

We can now write Equation~\ref{eq:mui-approx} as
\begin{equation}
\mu_i = 2 \frac{1}{\sqrt{2\pi}\theta_0} 2\pi\alpha z^2 \int_x \diff x \Omega(x) \frac{1}{\left|\sin\theta_1\right|} \epsilon(\eta) P(\text{not scattered or absorbed} | \lambda_0, E, x, v) \mathrm{QE}(\lambda_0) e^{-\frac{\Delta^2(\lambda_0)}{4\theta_0^2}} \left(1 - \frac{1}{\beta^2 n^2(\lambda_0)}\right) f(x)
\end{equation}

The probability that a photon travels to the PMT without being scattered or absorbed can be calculated as follows
\begin{align}
P(\text{not scattered or absorbed} | \lambda, x) &=
P(\text{not scattered} | \lambda, x) P(\text{not absorbed} | \lambda, x) \\
&= \int_l^\infty\frac{1}{s(\lambda)}e^{-\frac{x}{s(\lambda)}}\int_l^\infty\frac{1}{a(\lambda)}e^{-\frac{x}{a(\lambda)}} \\
&= e^{-\frac{l}{s(\lambda) + a(\lambda)}}
\end{align}
where $l$ is the distance to the PMT from the position $x$, $s(\lambda)$ is the
scattering length, and $a(\lambda)$ is the absorption length.

We can therefore write the expected charge as
\begin{equation}
\mu_i = 2 \frac{1}{\sqrt{2\pi}\theta_0} 2\pi\alpha z^2 \int_x \diff x \Omega(x) \frac{1}{\left|\sin\theta_1\right|} \epsilon(\eta) e^{-\frac{l(x)}{s(\lambda) + a(\lambda)}} \mathrm{QE}(\lambda_0) e^{-\frac{\Delta^2(\lambda_0)}{4\theta_0^2}} \left(1 - \frac{1}{\beta^2 n^2(\lambda_0)}\right) f(x)
\end{equation}

The last integral is calculated numerically when the likelihood is evaluated.

We now return to the likelihood and calculate the probability of observing a
given time. In principle, this depends on the number of photons hitting a PMT
since the PMT hit will only register the \emph{first} photoelectron which
crosses threshold. However, since this is expected to be a small effect, we
assume that the probability of observing a given time is independent of the
number of photons which hit the PMT, i.e.
\begin{equation}
P(t_j | n, E, \vec{x}, \vec{v}, t_0) \simeq P(t_j | n \geq 1, E, \vec{x}, \vec{v}, t_0)
\end{equation}

We first condition on the \emph{true} time at which the photon hits the PMT
\begin{equation}
P(t_j | n \geq 1, E, \vec{x}, \vec{v}, t_0) = \int_{t_j'} \diff t P(t_j | t_j') P(t_j' | n \geq 1, E, \vec{x}, \vec{v}, t_0)
\end{equation}
where we used the fact that the probability of a measured time only depends on the true PMT hit time.

Now, we integrate over the track
\begin{equation}
P(t_j | n \geq 1, E, \vec{x}, \vec{v}, t_0) = \int_{t_j'} \diff t_j' P(t_j | t_j') \int_x \diff x P(t_j', x | n \geq 1, E, \vec{x}, \vec{v}, t_0)
\end{equation}
where $x$ here stands for the event that a photon emitted at a distance $x$ along the track makes it to the PMT.

We now use Bayes theorem to rewrite the last probability
\begin{align}
P(t_j | n \geq 1, E, \vec{x}, \vec{v}, t_0) &= \int_{t_j'} \diff t_j' P(t_j | t_j') \int_x \diff x P(t_j', x | n \geq 1, E, \vec{x}, \vec{v}, t_0) \\
&= \int_{t_j'} \diff t_j' P(t_j | t_j') \int_x \diff x P(t_j' | x, n \geq 1, E, \vec{x}, \vec{v}, t_0) P(x | n \geq 1, E, \vec{x}, \vec{v}, t_0) \\
\end{align}
The first term in the integral is just a delta function (up to slight differences due to dispersion) since we are assuming direct light
\begin{align}
P(t_j | n \geq 1, E, \vec{x}, \vec{v}, t_0)
&= \int_{t_j'} \diff t_j' P(t_j | t_j') \int_x \diff x \delta\left(\frac{l(x)n(\lambda_0)}{c}-t_j'\right) P(x | n \geq 1, E, \vec{x}, \vec{v}, t_0) \\
\end{align}

We now use Bayes theorem to rewrite the last term
\begin{align}
P(t_j | n \geq 1, E, \vec{x}, \vec{v}, t_0)
&= \int_{t_j'} \diff t_j' P(t_j | t_j') \int_x \diff x \delta\left(\frac{l(x)n(\lambda_0)}{c}-t_j'\right) \frac{P(n \geq 1 | x, E, \vec{x}, \vec{v}, t_0) P(x | E, \vec{x}, \vec{v}, t_0)}{P(n \geq 1 | E, \vec{x}, \vec{v}, t_0)} \\
&= \int_{t_j'} \diff t_j' P(t_j | t_j') \int_x \diff x \delta\left(\frac{l(x)n(\lambda_0)}{c}-t_j'\right) \frac{P(x | E, \vec{x}, \vec{v}, t_0)}{P(n \geq 1 | E, \vec{x}, \vec{v}, t_0)} \\
&= \int_{t_j'} \diff t_j' P(t_j | t_j') \int_x \diff x \delta\left(\frac{l(x)n(\lambda_0)}{c}-t_j'\right) \frac{P(x | E, \vec{x}, \vec{v}, t_0)}{1 - e^{-\mu_j}} \\
&= \int_{t_j'} \diff t_j' P(t_j | t_j') \int_x \diff x \delta\left(\frac{l(x)n(\lambda_0)}{c}-t_j'\right) \frac{\mu_j(x)}{1 - e^{-\mu_j}} \\
&= \frac{1}{1 - e^{-\mu_j}} \int_x \diff x \mu_j(x) \int_{t_j'} \diff t_j' P(t_j | t_j') \delta\left(\frac{l(x)n(\lambda_0)}{c}-t_j'\right) \\
\end{align}

We assume the transit time spread is equal to a gaussian (we ignore the pre and late pulsing)
\begin{align}
P(t_j | n \geq 1, E, \vec{x}, \vec{v}, t_0)
&= \frac{1}{1 - e^{-\mu_j}} \int_x \diff x \mu_j(x) \int_{t_j'} \diff t_j' \frac{1}{\sqrt{2\pi}\sigma_t} e^{-\frac{(t_j-t_j')^2}{2\sigma_t^2}} \delta\left(\frac{l(x)n(\lambda_0)}{c}-t_j'\right) \\
&= \frac{1}{1 - e^{-\mu_j}} \frac{1}{\sqrt{2\pi}\sigma_t} \int_x \diff x \mu_j(x) e^{-\frac{(t_j-t_0(x))^2}{2\sigma_t^2}}
\end{align}
where in the last expression we define
\begin{equation}
t_0(x) \equiv \frac{l(x)n(\lambda_0)}{c}
\end{equation}

\chapter{Backgrounds}
\section{External Muons}
Both cosmic ray muons and muons created from atmospheric neutrinos interacting
in the surrounding rock present a background for this analysis. In both cases,
it is necessary to cut events which start \emph{outside} the PSUP and enter the
detector.

During SNO, these events were cut using the MUON cut which tagged events with
at least 150 hits, 5 or more outward-looking (OWL) PMT hits, and with a time
RMS of less than 90 nanoseconds. This cut would have a negligible sacrifice for
any contained atmospheric or dark matter candidate events, but could
potentially cut events which produce an energetic muon which then exits the
detector. Therefore, I have slightly modified this cut to \emph{also} require
that at least 1 OWL tube is both early and has a high charge relative to the
nearby normal PMTs. We define a early and high charge tube by creating an array
of the ECA calibrated hit times (we can't use PCA calibrated times since the
OWL tubes were never calibrated via PCA) and of the best uncalibrated charge
(FIXME: footnote?) for all normal PMTs within 3 meters of each hit OWL PMT. We
then compute the median charge and time for these normal PMTs. We then compute
how many OWL PMT hits are \emph{both} earlier than the median normal PMT time
and have a higher charge than the surrounding PMTs. If at least 1 OWL PMT hit
satisfies this criteria and all the other criteria from the SNO MUON cut are
satisifed (except the time RMS part) then it's tagged as a muon.

\section{Noise Events}

There are several sources of noise events which refers to events triggered by
sources which do not actually create light in the detector. The two most common
sources are "ringing" after large events and electrical pickup on deck.

These events are tagged by the QvNHIT and ITC cuts which are identical to their
SNO counterparts aside from minor updates\footnote{The ITC cut uses the pt1
time which is the time without the charge walk calibration since otherwise the
cut may fail to tag an event which consists of mostly electronics noise which
has charge too low to apply PCA. The QvNHIT cut does not require good
calibrations for the hits for a similar reason.}.

\section{Neck Events}

Neck events are caused by light produced in or leaking through the glove box on
top of the detector\cite{sonley}. The SNO neck event cut is defined
as\cite{snoman_companion}:

\begin{quotation}
This cuts events containing neck tubes. It requires that either both tubes in
the neck fire, or that one of those tubes fires and it has a high charge and is
early. High charge is defined by a neck tube having a pedestal subtracted
charge greater than 70 or less than -110. Early if defined by the neck tube
having an ECA time 70ns or more before the average ECA time of the PSUP PMTS
with z les than 0. After the cable changes to the neck tubes this time
difference changes to 15ns. 
\end{quotation}

Similarly to the MUON cut, I've used these criteria but added an additional
requirement to avoid tagging high energy upwards going events. The NECK cut I
use also has a requirement that 50\% of the hit PMTs must have a z coordinate
of less than 4.25 meters \emph{or} 50\% of the ECA calibrated QHS charge must
be below z = -4.25 meters.

\section{Flashers}

Flashers are probably the most difficult and common source of instrumental
background for this analysis. A flasher event occurs when there is an
electrical short in the PMT base or dynode stack which causes light to be
emitted from the PMT and hit the opposite side of the
detector\textsuperscript{[citation needed]}. Because this event is caused by
actual light in the detector it is particularly hard to cut while also
maintaining a small signal sacrifice.

The cut algorithm is sufficiently complex that it is easier to describe in
pseudocode. A description of the algorithm is shown in
Algorithm~\ref{flasher_algorithm}.

\begin{algorithm}
\caption{Flasher Cut Algorithm}
\label{flasher_algorithm}
\begin{algorithmic}
  \IF{nhit $< 31$} \RETURN 0 \ENDIF

  \COMMENT{This condition is similar to the SNO QvT cut except we require that 70\% of the normal hit PMTs be 12 meters from the high charge channel and that 70\% of the normal hit PMTs be at least 50 ns after the high charge channel.}

  \IF{highest QLX $>$ second highest QLX $+ 80$}
     \STATE {Collect all hit times from the same slot as the high charge channel and compute the median hit time}
     \IF{At least 4 hits in the slot \AND 70\% of the normal hit PMTs with good calibration are more than 12 meters from the high charge channel \AND 70\% of the normal hit PMTs with good calibration are more than 50 ns after the median hit time in the slot}
       \RETURN 1
     \ENDIF
  \ENDIF
  \FOR{All PC with at least 4 hits}
    \STATE {Collect the QHS, QHL, and QLX charges and the ECA calibrated hit times (EPT) for each PMT in the PC sending charge values below 300 to 4095}
    \STATE {$t \leftarrow \textrm{median}(\textrm{EPT})$}
    \STATE {$\textrm{QHS}_1 \leftarrow \textrm{max}(\textrm{QHS})$}
    \STATE {$\textrm{QHL}_1 \leftarrow \textrm{max}(\textrm{QHL})$}
    \STATE {$\textrm{QLX}_1 \leftarrow \textrm{max}(\textrm{QLX})$}
    \STATE {$\textrm{QHS}_2 \leftarrow \textrm{second highest}(\textrm{QHS})$}
    \STATE {$\textrm{QHL}_2 \leftarrow \textrm{second highest}(\textrm{QHL})$}
    \STATE {$\textrm{QLX}_2 \leftarrow \textrm{second highest}(\textrm{QLX})$}
    \IF{$\textrm{QHS}_1 > \textrm{QHS}_2 + 1000$}
      \IF{70\% of the normal hit PMTs with good calibration are more than 12 meters from the high charge channel \AND 70\% of the normal hit PMTs with good calibration are more than 50 ns after $t$}
        \RETURN 1
      \ENDIF
    \ELSIF{$\textrm{QHL}_1 > \textrm{QHL}_2 + 1000$}
      \IF{70\% of the normal hit PMTs with good calibration are more than 12 meters from the high charge channel \AND 70\% of the normal hit PMTs with good calibration are more than 50 ns after $t$}
        \RETURN 1
      \ENDIF
    \ELSIF{$\textrm{QLX}_1 > \textrm{QLX}_2 + 80$}
      \IF{70\% of the normal hit PMTs with good calibration are more than 12 meters from the high charge channel \AND 70\% of the normal hit PMTs with good calibration are more than 50 ns after $t$}
        \RETURN 1
      \ENDIF
    \ELSE
      \FOR{All normal PMT channels \emph{not} hit in PC}
        \IF{more hits in slot than surrounding 4 meters or median hit time in slot is 10 ns earlier than PMTs within 4 meters}
          \IF{70\% of the normal hit PMTs with good calibration are more than 12 meters from the high charge channel \AND 70\% of the normal hit PMTs with good calibration are more than 50 ns after $t$}
            \RETURN 1
          \ENDIF
        \ENDIF
      \ENDFOR
    \ENDIF
  \ENDFOR
  \RETURN 0
\end{algorithmic}
\end{algorithm}

\section{Breakdowns}

Breakdowns are very similar to flashers except that they produce \emph{much}
more light\footnote{In fact, I think there is a continuous spectrum between
flashers and breakdowns, but the distinction is still helpful since the way to
tag the two are very different}.

Since breakdowns often cause many of the electronics to saturate, it is
\emph{very} difficult to find a single common characteristic on which to cut.
However, the one thing that does seem to be common among almost all breakdowns
is that the crate with the channel that breaks down all has pickup from the
breaking down channel and thus comes much earlier in the event than the rest of
the PMT hits.

Therefore, the breakdown cut tags any event which has at least 1000 PMT hits
and in which the crate with the highest median TAC has at least 256 hits and is
500 TAC counts away from the next highest crate (with at least 20 hits).

Occasionally a breakdown is so big that it causes issues with the TAC
measurement and many of them end up reading outside of the linear TAC region.
Therefore we also tag any events in which less than 70\% of the PMT hits have a
TAC value above 400.

\appendix
\chapter{Poisson Binomial}

Suppose we have a Poisson process whose output is then subject to a binomial
process. For example, we expect $\mu$ background events on average and we can
detect them with probability $p$. What is the probability of detecting $n$
background events?

\begin{align}
p(n) &=                                 \sum_{N=n}^{\infty} P(n|N) P(N) \\
     &=                                 \sum_{N=n}^{\infty} \frac{N!}{n!(N-n)!} p^n (1-p)^{N-n} e^{-\mu}\frac{\mu^N}{N!} \\
     &=                                 \sum_{N=n}^{\infty} \frac{1}{n!(N-n)!} p^n (1-p)^{N-n} e^{-\mu}\mu^N \\
     &= e^{-\mu}   \frac{p^n}{n!}       \sum_{N=n}^{\infty} \frac{1}{(N-n)!} (1-p)^{N-n} \mu^N \\
     &= e^{-\mu}   \frac{(\mu p)^n}{n!} \sum_{N=n}^{\infty} \frac{\left(\mu (1-p)\right)^{N-n}}{(N-n)!}  \\
     &= e^{-\mu p} \frac{(\mu p)^n}{n!}
\end{align}

Therefore the end result is a Poisson distribution with mean $\mu p$.

\begin{thebibliography}{9}
\bibitem{grossman2017}
    Grossman, et al. \textit{Self-Destructing Dark Matter}. \href{https://arxiv.org/abs/1712.00455}{{\tt arXiv:1712.00455}}. Dec 2017.
\bibitem{pepin2016}
    M. Pepin. \textit{Low-Mass Dark Matter Search Results and Radiogenic Backgrounds for the Cryogenic Dark Matter Search}. \url{http://hdl.handle.net/11299/185144}. Dec 2016.
\bibitem{caldwell2015}
    T. Caldwell. \textit{Searching for Dark Matter with Single Phase Liquid Argon}. \url{https://repository.upenn.edu/dissertations/1632}. 2015.
\bibitem{pdg2017}
    C. Patrignani et al. (Particle Data Group), Chin. Phys. C, 40, 100001 (2016) and 2017 update.
\bibitem{sonley}
    T. Sonley. \textit{A Measurement of the Atmospheric Neutrino Flux and Oscillation Parameters at the Sudbury Neutrino Observatory}. Feb 2009. \url{https://www.sno.phy.queensu.ca/sno/papers/Sonley_phd_physics_2009.pdf}. Access Date: Oct 5, 2019.
\bibitem{snoman_companion}
    \textit{SNOMAN Companion}. Last updated: Nov. 8, 2006. \url{http://hep.uchicago.edu/~tlatorre/snoman_companion/}. Access Date: Oct 5, 2019.
\end{thebibliography}
\end{document}