documentation.tex

\documentclass[11pt, oneside]{article}
\usepackage{geometry}
\geometry{a4paper}
\usepackage{graphicx}
\usepackage{amssymb}

\usepackage{amsmath}
\usepackage{hyperref}
\usepackage[utf8]{inputenc}
\usepackage{multirow}
\usepackage{relsize}
\usepackage{siunitx}
\usepackage{xspace}


\renewcommand{\d}{\mathrm{d}}
\newcommand{\tdiff}[3][\empty]{\ifx\empty#1
		\frac{\d\,#2}{\d #3}
	\else
		\frac{\d^{#1}\,#2}{\d #3^{#1}}
	\fi} % (#1-th) total differential (in equations)
\newcommand{\tdiffx}[3][\empty]{\ifx\empty#1
		\d #2 / \d #3
	\else
		\d^{#1} #2 / \d #3^{#1}
	\fi} % (#1-th) total differential (in text)
\newcommand{\mean}[1]{\ensuremath{\langle #1 \rangle}}
\newcommand{\nue}{\ensuremath{\nu_e}\xspace}
\newcommand{\numu}{\ensuremath{\nu_\mu}\xspace}
\newcommand{\nutau}{\ensuremath{\nu_\tau}\xspace}
\newcommand{\nux}{\ensuremath{\nu_x}\xspace}
\newcommand{\nuebar}{\ensuremath{\bar{\nu}_e}\xspace}
\newcommand{\numubar}{\ensuremath{\bar{\nu}_\mu}\xspace}
\newcommand{\nutaubar}{\ensuremath{\bar{\nu}_\tau}\xspace}
\newcommand{\nuxbar}{\ensuremath{\bar{\nu}_x}\xspace}


\title{Documentation for sntools 1.5\footnote{See \url{https://github.com/SNEWS2/sntools} for the most recent version.}}
\author{Jost Migenda\footnote{King’s College London, \url{jost.migenda@kcl.ac.uk}}}
\date{\today}

\begin{document}
\maketitle
\setcounter{tocdepth}{2}
\tableofcontents
\clearpage


\section{Introduction}

sntools is a Monte Carlo event generator for supernova and pre-supernova neutrino interactions.

It is a command line application that uses detailed time- and energy-dependent neutrino fluxes (provided by various supernova models) to generate interactions within the detector volume and write them to event files that can be used as an input for a full detector simulation.
sntools was originally developed for Hyper-Kamiokande~\cite{Migenda2019} using core-collapse supernova models and later extended to support different detectors, detector materials and to use pre-supernova models.

Additionally, sntools can be used as a Python library that implements neutrino cross sections. See section~\ref{sec:interaction-channels}.

\subsection{Installation Instructions}
% sntools should support *at least* the python/scipy/numpy versions required by NEP 29:
% https://numpy.org/neps/nep-0029-deprecation_policy.html

First, make sure you have Python 3.9 or higher installed on your computer.
(If you don’t, you can get it from \href{https://www.python.org}{python.org} or as part of the \href{https://www.anaconda.com/products/individual}{Anaconda} distribution.)

In a terminal, run \texttt{pip install sntools} to install the latest version of sntools and all dependencies.
Then, run \texttt{python -c 'import sntools; sntools.setup()'} to check whether sntools is working and to download the sample input file used below.


\subsection{Quick Start Guide}

To generate your first supernova neutrino events with sntools, open a terminal window and type in the following command:

\texttt{sntools fluxes/intp2001.data --format nakazato --endtime 50}

This generates events for the first \SI{50}{ms} of a supernova based on the neutrino flux given in the input file \texttt{fluxes/intp2001.data}\footnote{This sample file, which is for a \SI{20}{M_\odot} progenitor from~\cite{Nakazato2013}, is included as part of sntools. See \url{http://asphwww.ph.noda.tus.ac.jp/snn/} for more information.} (which is in the \texttt{nakazato} format) and writes them to a file named \texttt{outfile.kin} in the current directory.

A more realistic usage that demonstrates more of sntools’ capabilities looks like this:

\texttt{sntools fluxes/intp2001.data --format nakazato --channel es --detector SuperK --distance 20 --verbose --output intp2001es.kin}

This uses the neutrino flux given in the same input file to generate neutrino-electron elastic scattering events in Super-Kamiokande for a supernova at \SI{20}{kpc} distance.
It also produces more verbose output, which lets you see that sntools generates events separately for \nue, \nuebar, \nux and \nuxbar (which have different fluxes and cross sections), before merging them into an output file named \texttt{intp2001es.kin}.


The input file is a required argument. The following optional arguments are available:

\begin{description}
\item[\texttt{--format <value>}] Format of input file(s). See section~\ref{sec:input-formats}.
\item[\texttt{--output <value>}] Name of the output file.
\item[\texttt{--mcformat <value>}] Format of output file. Can be \texttt{NUANCE} (used e.\,g. by Hyper-Kamiokande), \texttt{RATPAC} (used e.\,g. by SNO+ and WATCHMAN) or \texttt{JUNO\_ROOT} (used by the JUNO internal software).\footnote{For a description of these formats see \url{http://neutrino.phy.duke.edu/nuance-format/} and \url{https://rat.readthedocs.io/en/latest/generators.html\#external}. For more details on the JUNO SNEvent internal structure, contact \url{marta.colomer@ulb.be}.}

\item[\texttt{--detector <value>}] Detector configuration. See section~\ref{sec:detector-configurations}.
\item[\texttt{--channel <value>}] Interaction channel to generate events for. See section~\ref{sec:interaction-channels}.
\item[\texttt{--transformation <value>}] Transformation between initial flux (inside the supernova) and detected flux on Earth. See section~\ref{sec:transformation}.
\item[\texttt{--distance <value>}] Distance to supernova (in kpc).
\item[\texttt{--starttime <value>} and \texttt{--endtime <value>}] Generate events in a certain time range. These times are given in milliseconds (for core-collapse supernova) or minutes (for pre-supernova). When using a pre-supernova model, all times before the core collapse of the star are given as negative. If these arguments aren’t specified, events are generated for the full time range for which fluxes are given in the input file.
\item[\texttt{--binsize <value>}] For use with presn models, specify the width of time bins used during event generation (in seconds). See section~\ref{sec:pre-supernova-time-binning}.
\item[\texttt{--randomseed <value>}] Set an integer as a seed for the random number generator, to generate events reproducibly.
\item[\texttt{--maxworkers <value>}] Maximum number of parallel processes. sntools will generate events for up to this many channels in parallel, which may improve performance on multi-core CPUs.
\item[\texttt{--verbose}] Print more detailed output.
\item[\texttt{--version}] Print the current version number and exit.
\end{description}

The command \texttt{sntools -h} displays a more detailed overview over all arguments, including their default values.



\section{More Details} \label{sec:more-details}

\begin{figure}[htbp]
	\centering
	\includegraphics[scale=0.52]{sntools-overview.pdf}
	\caption[Overview over the structure of sntools]{Overview over the structure of sntools. See text for a detailed description.}
	\label{fig-sim-sntools-overview}
\end{figure}

An overview over the structure of sntools is given in figure~\ref{fig-sim-sntools-overview}.
% sntools is written in Python and uses the scipy and numpy libraries~\cite{Virtanen2020,Walt2011}, which implement many numerical calculations in Fortran for performance reasons.
% Additionally, critical parts of sntools were tuned for increased performance.
%On a current desktop computer, sntools needs \SI[parse-numbers=false]{\mathcal{O}(10)}{min} to generate events for a supernova at the fiducial distance of \SI{10}{kpc} in Hyper-Kamiokande.

The main user interface is provided by the file \texttt{genevts.py}, which parses command line arguments and then calls code in \texttt{channel.py} to generate events.
Events are generated separately for each combination of interaction channel and input species.\footnote{sntools uses the species \nue, \nuebar, \nux and \nuxbar, with \nux (\nuxbar) representing any one of \numu and \nutau (\numubar and \nutaubar). Since the energy of supernova neutrinos is too low to produce $\mu^\pm$ or $\tau^\pm$ in a detector, neutrinos with those two flavours interact in the same ways.}
The code first calculates the total number of events expected in each specified time bin. For core-collapse supernova models, the bin size is fixed at \SI{1}{ms}. For pre-supernova model use, the binsize is an optional argument with a default value of \SI{1}{s}.
The number of events per ms is calculated first and then multiplied by the bin width.
The number of events per ms is given by
\begin{equation}
N(t) = \mathop{\mathlarger{\mathlarger{\iint}}} \tdiff{\Phi (t, E_\nu)}{E_\nu} \tdiff{\sigma (E_\nu, E_e)}{E_e}\, \d E_e\,\d E_\nu,
\end{equation}
where $\Phi (t, E_\nu)$ is the neutrino flux and $\sigma (E_\nu, E_e)$ is the cross section of the current interaction channel.
It then picks the actual number of events to generate within that time bin from a Poisson distribution with expectation value $N(t)$.
Finally, it generates events by rejection sampling from the energy spectrum of neutrino interactions at that time and the distribution of outgoing particle directions.

The event generation code relies on a plug-in architecture to support various different input formats and interaction channels.
Input format plug-ins provide functions that read in the data from an input file and return the number luminosity as a function of time and energy.
Interaction channel plug-ins specify properties of the interaction channel, like the neutrino species that undergo this interaction, and provide functions to calculate quantities like the differential cross section $\tdiffx{\sigma (E_\nu, E_e)}{E_e}$ or the kinematically allowed energy range.
This modular design makes sntools easily extensible, with roughly 100 lines of code required to add a new input format or interaction channel.

Finally, \texttt{genevts.py} collects the events generated in all interaction channels and writes them to an output file, which can be used as input for a full detector simulation.

The rest of this chapter describes the various components of the code and the physics behind them in more detail.


\subsection{Detector Materials}
sntools currently supports three detector materials: water, liquid scintillator and water-based liquid scintillator.
Water is assumed to consist of an oxygen-16 nucleus and two hydrogen nuclei (free protons).
Generic liquid scintillator is assumed to consist of $n$ carbon-12 nuclei and $2n$ hydrogen nuclei (free protons).
In addition, two specific liquid scintillators\footnote{both linear alkylbenzene (LAB), a mixture of C$_6$H$_5$C$_n$H$_{2n+1}$ molecules with $9<n<14$} are implemented, which match the materials used in SNO+~\cite{anderson2021development} and JUNO.
Water-based liquid scintillator is assumed to be a mixture of $x\,\%$ generic liquid scintillator and $(100-x)\,\%$ water.
%Each material supports the appropriate detection channels presented in section~\ref{sec:interaction-channels}.

The detector material doesn’t need to be specified explicitly; it is determined by the detector configuration.


\subsection{Detector Configurations}\label{sec:detector-configurations}
Select a detector configuration by using the \texttt{--detector <value>} command line option. The following detector configurations are currently implemented:

\begin{description}
\item[\texttt{HyperK}] Inner detector of Hyper-Kamiokande~\cite{HKDR2018}. Cylinder with a height of \SI{65.8}{m} and a diameter of \SI{64.8}{m}, filled with water.\footnote{In September 2019, the dimensions were changed slightly compared to the 2018 Design Report. The old configuration is available as \texttt{HyperKDR}.}
\item[\texttt{SuperK}] Inner detector of Super-Kamiokande~\cite{Fukuda2003}. Cylinder with a height of \SI{36.2}{m} and a diameter of \SI{33.7}{m}, filled with water.
\item[\texttt{WATCHMAN}] WATCHMAN detector~\cite{Askins2015}. Cylinder with a height and diameter of \SI{12.8}{m}, filled with water.
\item[\texttt{WATCHMAN-LS}] Same as \texttt{WATCHMAN}, but filled with liquid scintillator.
\item[\texttt{WATCHMAN-WbLS}] Same as \texttt{WATCHMAN}, but filled with water-based liquid scintillator (3\,\% LS, 97\,\% water).
\item[\texttt{THEIA25}] THEIA-25 geometry~\cite{Askins2020}. Box with dimensions $\SI{20}{m}\times\SI{18}{m}\times\SI{70}{m}$, filled with \SI{25}{kt} of water-based liquid scintillator (10\,\% LS, 90\,\% water).
\item[\texttt{THEIA100}] THEIA-100 geometry~\cite{Askins2020}. Cylindrical detector filled with \SI{100}{kt} of water-based liquid scintillator (10\,\% LS, 90\,\% water).
\item[\texttt{SNOplusAV}] Inner acrylic vessel of SNO+ \cite{albanese2021sno+}. Sphere with a radius of \SI{6}{m}, filled with linear alkylbenzene (LAB) liquid scintillator.
\item[\texttt{SNOplusEW}] External water of SNO+ \cite{albanese2021sno+}. Two concentric spheres of radii \SI{9}{m} and \SI{6.05}{m}. The outer sphere is filled with water, the inner sphere is empty.
\item[\texttt{JUNO}] Central detector of JUNO \cite{JUNO:2021vlw}. Sphere with a radius of \SI{17.7}{m}, filled with \SI{20}{kton} of LAB liquid scintillator.
\end{description}


\subsection{Interaction Channels} \label{sec:interaction-channels}
sntools supports multiple different interaction channels described in this section.
By default, it will generate events across all channels that are available in the selected detector, but it can be restricted to a single channel by using the \texttt{--channel <value>} command line argument, where \texttt{<value>} can be one of \texttt{ibd}, \texttt{es}, \texttt{ps}, \texttt{o16e}, \texttt{o16eb}, \texttt{c12e}, \texttt{c12eb} and \texttt{c12nc}.

For all supported channels, sntools includes code to calculate the differential cross sections, outgoing particle energy as a function of neutrino energy and other quantities.
These could be used as a library from other Python code, e.\,g. as follows:\\
\texttt{>> from sntools.interaction\_channels import ibd}\\
\texttt{>> c = ibd.Channel('eb')}\\
\texttt{>> c.bounds\_eE(eNu=20)  \# min/max energy of outgoing positron in MeV}\\
\texttt{[17.941220926704954, 18.70577898514169]}

Type \texttt{help(ibd)} on the Python command line for documentation or see \href{https://github.com/SNOwGLoBES/snowglobes/pull/12/commits/f1cf96d8c326b99ac35474bb6a36ef4d1fde809d}{this SNOwGLoBES pull request} for a full usage example.


\subsubsection{\texttt{ibd}: Inverse Beta Decay}
In both water and liquid scintillator, inverse beta decay (IBD; $\nuebar + p \rightarrow n + e^+$) is the dominant interaction channel for supernova neutrinos due to its relatively high cross section and low energy threshold of $E_\nu^\text{thr} \approx \SI{1.8}{MeV}$, as well as the large number of free protons in the detector.
The observed energy of IBD events is closely related to the neutrino energy, making this an excellent channel to reconstruct the \nuebar spectrum.

In sntools, IBD is implemented using the full tree-level cross section calculated in reference~\cite{Strumia2003} and including radiative corrections based on the approximation from reference~\cite{Kurylov2003}.\footnote{That calculation uses the limit $m_e \rightarrow 0$. This approximation is accurate to better than 0.1\,\% above $E_e = \SI{1}{MeV}$.}


\subsubsection{\texttt{es}: Neutrino-Electron Scattering}
In both water and liquid scintillator, elastic neutrino-electron scattering  ($\nu + e^- \rightarrow \nu + e^-$) is a subdominant interaction channel due to its low cross section, which is only partially compensated by the large number of electrons in the detector.
Elastic scattering events make up only a few per cent of all events but their angular distribution is strongly peaked into a forward direction, pointing away from the supernova.
In detectors that can reconstruct the direction of scattered electrons, this channel can therefore be used to determine the direction of a supernova~\cite{Abe2016,HKDR2018}.

Elastic scattering is sensitive to all neutrino flavours.
However, the cross section of \nue and \nuebar, which can interact through both neutral and charged currents, is higher than that of \nux and \nuxbar, which can interact only through neutral currents.

In sntools, elastic scattering is implemented using the tree-level cross sections from standard electroweak theory calculated by ’t~Hooft~\cite{t-Hooft1971} and including one-loop electroweak and QCD corrections as well as QED radiative corrections as calculated in reference~\cite{Bahcall1995}.

\subsubsection{\texttt{ps}: Neutrino-Proton Scattering}
In liquid scintillator, neutrino-proton scattering ($\nu + p \rightarrow \nu + p$) is an additional available subdominant interaction channel. The recoiling proton energy spectrum is soft meaning the deposited energy in the detector is relatively low ($<5$~MeV). These protons are invisible in water as they're below the Cherenkov threshold. Furthermore, the slow heavily ionizing protons lose energy very quickly, quenching the effective signal. Despite this, over a realistic detection threshold ($\mathcal{O}(100)$s~keV), the neutrino-proton scattering yield can be the second largest available to liquid scintillator detectors below IBD. 

Neutrino-proton scattering is a neutral-current interaction and is available to all (anti-)neutrino flavours. The proton recoil spectrum can also provide spectral information about the incoming neutrino, which could solve a long-standing problem of how to separately measure the total energy and temperature of $\nu_{\mu}$, $\nu_{\tau}$, $\bar{\nu}_{\mu}$, and $\bar{\nu}_{\tau}$ \cite{beacom2002detection}.

In sntools, the implementation of neutrino-proton scattering is based a prediction by Beacom, Farr, and Vogel \cite{beacom2002detection}. Here the cross-section is calculated directly from the Standard Model, which has been extensively verified by GeV-scale experiments \cite{ahrens1987measurement}. 

\subsubsection{\texttt{o16e} and \texttt{o16eb}: Charged-Current Interactions on $^{16}$O}
In water, charged-current interactions of \nue and \nuebar on $^{16}$O nuclei,
\begin{align}
\nue + ^{16}\text{O} &\rightarrow e^- + X \\
\nuebar + ^{16}\text{O} &\rightarrow e^+ + X,
\end{align}
are a subdominant interaction channel.\footnote{\label{fn:cc-interactions}Charged-current interactions of other neutrino species do not occur, since the energy of supernova neutrinos is too small to produce muons or $\tau$ leptons.}
Due to the high energy threshold of both interactions of approximately \SI{15}{MeV} and \SI{11}{MeV}, respectively, as well as the steep energy dependence of the cross sections, the number of events in each channel is a very sensitive probe of the high-energy tail of the supernova neutrino flux.
It may vary by more than two orders of magnitude depending on the supernova models and oscillation scenario.

In sntools, the cross section for both interaction channels is based on a new shell model calculation~\cite{Suzuki2018} which selected 42 different nuclear states and calculated their respective partial cross sections.
To simplify the evaluation of the cross section, the implementation in sntools uses the four-group fit presented in reference~\cite{Nakazato2018}.
This fit matches the cross sections calculated from the full set of nuclear states to within a few per cent at neutrino energies of up to \SI{100}{MeV}.
For a typical supernova neutrino flux, the difference in the resulting event spectra when using the four groups instead of all 42 nuclear states is very small.


\subsubsection{\texttt{c12e} and \texttt{c12eb}: Charged-Current Interactions on $^{12}$C}
In liquid scintillator, charged-current interactions of \nue and \nuebar on $^{12}$C nuclei,
\begin{align}
\nue + ^{12}\text{C} &\rightarrow e^- + X \\
\nuebar + ^{12}\text{C} &\rightarrow e^+ + X,
\end{align}
are a subdominant interaction channel.\footnote{Footnote \ref{fn:cc-interactions} applies here as well.}
Due to the high energy threshold of both interactions of approximately \SI{17}{MeV} and \SI{14}{MeV}, respectively, as well as the steep energy dependence of the cross sections, the number of events in each channel is a very sensitive probe of the high-energy tail of the supernova neutrino flux.

In sntools, the cross section for both interaction channels is based on the calculation in reference~\cite{Kolbe1999}.
To interpolate between the tabulated cross sections, sntools uses the fit
\begin{equation}
\sigma (E_\nu) = \sigma_0 \cdot \left[ a_1 (E_\nu - E_\text{thr}) + a_2 (E_\nu - E_\text{thr})^2 + a_3 (E_\nu - E_\text{thr})^3 \right],
\end{equation}
where $E_\nu$ is the neutrino energy, $E_\text{thr}$ is the energy threshold of the interaction and $\sigma_0$ and $a_i$ are fit parameters.
The energy of the outgoing lepton is approximated as $E_e = E_\nu - E_\text{thr}$.


\subsubsection{\texttt{c12nc}: Neutral-Current Interactions on $^{12}$C}
In liquid scintillator, neutral-current interactions on $^{12}$C nuclei,
\begin{align}
\nu + ^{12}\text{C} \rightarrow \nu' + &^{12}\text{C}^*\\
	&^{12}\text{C}^* \rightarrow ^{12}\text{C} + \gamma
\end{align}
are a subdominant interaction channel that neutrinos and antineutrinos of all flavours contribute to equally.
In sntools, the cross section for this interaction channel is based on the theoretical calculation in~\cite{Donnelly1979}
\begin{equation}
\sigma (E_\nu) = \SI{1.08e-38}{cm^2} \left( \frac{E_\nu - E_\text{thr}}{m_N} \right)^2 \cdot \beta \kappa,
\end{equation}
where $E_\nu$ is the neutrino energy, $E_\text{thr} \approx \SI{15}{MeV}$ is the energy threshold of the interaction, $m_N \approx \SI{939}{MeV}$ is the nucleon mass and $\beta \kappa = 1.11$~\cite{Armbruster1998}.
The energy of the outgoing $\gamma$ is $E_\gamma = E_\text{thr}$.


\subsection{Treatment of Neutrino Flavour Conversion} \label{sec:transformation}

sntools implements different transformations between the original neutrino flux generated in the supernova and the flux detected on Earth, which can be selected by using the \texttt{--transformation <value>} command line argument.

In the simplest case, \texttt{NoTransformation}, no flavor transformation takes place. Apart from the geometrical factor $\frac{1}{4 \pi d^2}$ (which depends only on the distance $d$ of the supernova\footnote{This geometrical factor is always taken into account by sntools. For simplicity, I will omit it in the following.}), the flux of a neutrino species $\nu_i$ observed by a detector on Earth, $\Phi_{\nu_i}$, is identical to the flux originating within the supernova, $\Phi_{\nu_i}^0$, which is given by the input file(s).

\subsubsection{Adiabatic MSW Effect}
The two transformations \texttt{AdiabaticMSW\_NMO} and \texttt{AdiabaticMSW\_IMO} assume that neutrinos traverse a smoothly varying density profile%
\footnote{After the accretion phase, the revived shock front travels outwards and passes through the layer within the star where the adiabatic flavour conversion takes place.
This causes a sudden change in the matter and electron density and can severely impact the flavour conversion processes~\cite{Schirato2002}.
Since this occurs after the shock wave is revived, it mainly affects the late-time part of the supernova neutrino signal (e.\,g. for the Totani model, this effect is expected to become relevant more than \SI{1}{s} after core-bounce~\cite{Fogli2005}). Many simulations only cover the first few \SI{100}{ms} and would thus likely not exhibit this effect; making the assumption of a smoothly varying density profile reasonable.}
while exiting the star and undergo adiabatic flavour conversion via the MSW effect.


The resulting detected fluxes are linear combinations of the initial fluxes, which, for normal mass ordering (NMO), are given by~\cite{Dighe2000}
\begin{align}
\begin{split}
\Phi_{\nue} &= \sin^2 \theta_{13} \cdot \Phi_{\nue}^0 + \cos^2 \theta_{13} \cdot \Phi_{\nux}^0\\
\Phi_{\nuebar} &= \cos^2 \theta_{12} \cos^2 \theta_{13} \cdot \Phi^0_{\nuebar} + (1 - \cos^2 \theta_{12} \cos^2 \theta_{13}) \cdot \Phi^0_{\nuxbar} \\
2 \Phi_{\nux} &= \cos^2 \theta_{13} \cdot \Phi^0_{\nue} + (1 + \sin^2 \theta_{13}) \cdot \Phi^0_{\nux} \\
2 \Phi_{\nuxbar} &= (1 - \cos^2 \theta_{12} \cos^2 \theta_{13}) \cdot \Phi^0_{\nuebar} + (1 + \cos^2 \theta_{12} \cos^2 \theta_{13}) \cdot \Phi^0_{\nuxbar},
\end{split}
\end{align}
while for inverted mass ordering (IMO), they are
\begin{align}
\begin{split}
\Phi_{\nue} &= \sin^2 \theta_{12} \cos^2 \theta_{13} \cdot \Phi_{\nue}^0 + (1 - \sin^2 \theta_{12} \cos^2 \theta_{13}) \cdot \Phi_{\nux}^0\\
\Phi_{\nuebar} &= \sin^2 \theta_{13} \cdot \Phi_{\nuebar}^0 + \cos^2 \theta_{13} \cdot \Phi_{\nuxbar}^0\\
2 \Phi_{\nux} &= (1 - \sin^2 \theta_{12} \cos^2 \theta_{13}) \cdot \Phi_{\nue}^0 + (1 + \sin^2 \theta_{12} \cos^2 \theta_{13}) \cdot \Phi_{\nux}^0\\
2 \Phi_{\nuxbar} &= \cos^2 \theta_{13} \cdot \Phi_{\nuebar}^0 + (1 + \sin^2 \theta_{13}) \cdot \Phi^0_{\nuxbar}.
\end{split}
\end{align}

In both cases, the factor of 2 in the last two equations accounts for the fact that \nux (\nuxbar) includes the fluxes of \numu and \nutau (\numubar and \nutaubar).
These equations assume purely adiabatic transition (corresponding to $P_H = 0$ in~\cite{Dighe2000,Fogli2005}) as explained below.
Values for $\theta_{12}$ and $\theta_{13}$ are taken from the Particle Data Group~\cite{PDG2020}.

In cases where the detected flux is a mixture of original fluxes of different species, sntools generates events for each original species separately with the appropriate weighting factor applied.
For example, when generating inverse beta decay events using the \texttt{AdiabaticMSW\_NMO} transformation, sntools will generate events first using the input flux $\Phi_{\nuebar} = \cos^2 \theta_{12} \cos^2 \theta_{13} \cdot \Phi^0_{\nuebar}$ and then using the input flux $\Phi_{\nuebar} = (1 - \cos^2 \theta_{12} \cos^2 \theta_{13}) \cdot \Phi^0_{\nuxbar}$, before finally combining both sets of events into one output file.

\subsubsection{Additional Transformations}
sntools also supports the following flavour transformations, which are implemented by SNEWPY~\cite{Baxter2021}:
\texttt{CompleteExchange},
\texttt{NonAdiabaticMSWH},
\texttt{TwoFlavorDecoherence\_NMO},
\texttt{TwoFlavorDecoherence\_IMO},
\texttt{ThreeFlavorDecoherence}.
See reference~\cite{SNEWS:2021ezc} for a full description.
To use these transformations in sntools, prefix them with \texttt{SNEWPY-}. For example, to apply three flavor decoherence use sntools with the command line argument \texttt{--transformation SNEWPY-ThreeFlavorDecoherence}.

Some additional flavour transformations, like \texttt{NeutrinoDecay}, are implemented in SNEWPY but require some more work to be supported by sntools. If you are interested in these, \href{https://github.com/SNEWS2/sntools/issues/28}{please join the discussion on GitHub}!


\subsection{Input Formats} \label{sec:input-formats}
sntools supports multiple different input formats for the neutrino fluxes from a simulation.
To select one, use the \texttt{--format <value>} command line argument, where \texttt{<value>} can be one of \texttt{gamma}, \texttt{warren2020}, \texttt{nakazato}, \texttt{princeton} or \texttt{totani}.
This section will briefly describe each format and the processing steps necessary to calculate the spectral number luminosity.

In addition, sntools supports the following input formats, which are implemented by SNEWPY~\cite{Baxter2021,SNEWS:2021ezc}:
\texttt{Bollig\_2016},
\texttt{Fornax\_2021},
\texttt{Fornax\_2022},
\texttt{Kuroda\_2020},
\texttt{Mori\_2023},
\texttt{Nakazato\_2013},
\texttt{OConnor\_2015},
\texttt{Sukhbold\_2015},
\texttt{Tamborra\_2014},
\texttt{Walk\_2018},
\texttt{Walk\_2019} and
\texttt{Zha\_2021} for core-collapse supernova neutrino event generation,
and 
\texttt{Odrzywolek\_2010},
\texttt{Patton\_2017},
\texttt{Yoshida\_2016} and 
\texttt{Kato\_2017} for pre-supernova neutrino event generation.
SNEWPY also lets you download neutrino flux files from each of the corresponding simulations. See SNEWPY documentation for instructions.
To use these flux files in sntools, prefix them with \texttt{SNEWPY-}. For example, to use an \texttt{OConnor\_2015} input file, use sntools with the command line argument \texttt{--format SNEWPY-OConnor\_2015}.


All formats contain separate information on the three species \nue, \nuebar and \nux.
The fluxes of \nux and \nuxbar are typically assumed to be equal and most simulations don’t provide separate fluxes for \nux and \nuxbar.
For simplicity, references to each species are omitted in the following.

\subsubsection{Gamma Format}\label{sec:format-gamma}
Files in this format contain, for each time step $t_n$, the luminosity $L$, mean energy $\mean{E_\nu}$ and mean squared energy $\mean{E_\nu^2}$ of neutrinos.
A sample file is provided under \texttt{fluxes/sample-gamma.txt}.
%Each line is in the format:
%
%\texttt{$t$, $L_{\nue}$, $\mean{E_{\nue}}$, $\mean{E_{\nue}^2}$, $L_{\nuebar}$, $\mean{E_{\nuebar}}$, $\mean{E_{\nuebar}^2}$, $L_{\nux}$, $\mean{E_{\nux}}$, $\mean{E_{\nux}^2}$}

To reconstruct the neutrino spectrum from this, we assume that it is described by a normalized Gamma distribution~\cite{Keil2003,Tamborra2012} given by
\begin{equation}
f (E_\nu) = \frac{E_\nu^\alpha}{\Gamma (\alpha + 1)} \left( \frac{\alpha + 1}{A} \right)^{\alpha + 1} \exp \left[ - \frac{(\alpha + 1) E_\nu}{A} \right].
\end{equation}
In this formula, $A$ is an energy scale, while $\alpha$ determines the shape of the distribution: $\alpha = 2$ corresponds to a Maxwell-Boltzmann distribution, while $\alpha > 2$ corresponds to a “pinched” spectrum, which is more typical for neutrino spectra from supernovae.

The first two energy moments of the distribution are
\begin{align}
\mean{E_\nu}	&= \int_0^\infty \d E_\nu\, E_\nu f(E_\nu) = A\\
\mean{E_\nu^2}&= \int_0^\infty \d E_\nu\, E_\nu^2 f(E_\nu) = \frac{\alpha + 2}{\alpha + 1} A^2,
\end{align}
and therefore,
\begin{equation}
\alpha = \frac{\mean{E_\nu^2} - 2 \mean{E_\nu}^2}{\mean{E_\nu}^2 - \mean{E_\nu^2}}.
\end{equation}
Thus, the shape of the spectral number luminosity is uniquely determined by the mean energy \mean{E_\nu} and the mean squared energy \mean{E_\nu^2}, while the normalization is provided by $L / \mean{E_\nu}$.
To determine the spectral number luminosity at arbitrary times, each of the three parameters is interpolated separately before calculating the spectral number luminosity using the interpolated values.

\subsubsection{Warren2020 Format}
Files in this format contain, for each time step $t_n$, the luminosity $L$, mean energy $\mean{E_\nu}$ and RMS energy $\sqrt{\mean{E_\nu^2}}$ of neutrinos in the HDF5 format.
After reading data from the input file, it is processed as described in section~\ref{sec:format-gamma}.

\subsubsection{Nakazato Format}
Files in this format\footnote{See \url{http://asphwww.ph.noda.tus.ac.jp/snn/guide.pdf} for more details.} contain, for 20 energy bins $E_k$ during each time step $t_n$, the quantities $\Delta N_k (t_n) / \Delta E_k$ and $\Delta L_k (t_n) / \Delta E_k$, which reflect the number luminosity and luminosity at those energies, respectively.
For each energy bin, sntools calculates the mean energy within that bin, which is given by
\begin{equation}
\mean{E_k} = \frac{\quad  \frac{\Delta L_k (t_n)}{\Delta E_k}  \quad}{  \frac{\Delta N_k (t_n)}{\Delta E_k}  },
\end{equation}
and set the differential neutrino number flux at that energy to $\Delta N_k (t_n) / \Delta E_k$.
Finally, a linear interpolation in time and cubic spline interpolation in energy are used to determine the spectral number luminosity at an arbitrary time and energy. % as recommended by Nakazato\footnote{private communication, April 7, 2018}.

\subsubsection{Princeton Format}
Files in this format contain, for each time step $t_n$, the spectral luminosity $\tdiffx{L}{E}$ for 20 logarithmically spaced energy bins $E_k$.
I divide this by the central energy $\sqrt{E_k E_{k+1}}$ of the respective bin to get the spectral number luminosity at that energy.
Finally, a linear interpolation in time and cubic interpolation in energy are used to determine the spectral number luminosity at an arbitrary time and energy.

This follows the procedure described in reference~\cite{Seadrow2018}.
It is similar to that used for the Nakazato format described above, though with a different definition of the bin energy.

\subsubsection{Totani Format}
Files in this format contain, for each time step $t_n$, the total number of neutrinos emitted until that time, which makes it possible to calculate the number $N_n$ of neutrinos emitted since the previous time step.

For 20 energy bins $E_k$ per time step, a quantity $X_k$ is provided, which is proportional to the number of neutrinos emitted during that time step and in that energy bin.
I divide this by the width of each energy bin to get
\begin{equation}
X_k^\text{spec} = \frac{X_k}{E_{k+1} - E_k},
\end{equation}
which is proportional to the spectral number emission during that time step.

Integrating $X_k^\text{spec}$ over all energy bins and dividing it by that integral gives the spectral number emission during that time step normalized to 1, $X_k^\text{norm}$.
The spectral number luminosity at time $t_n$ and energy $E_k$ is given by
\begin{equation}
\tdiff{\text{NL} (t_n, E_k)}{E} = \frac{N_n}{t_n - t_{n-1}} \cdot X_k^\text{norm}.
\end{equation}
Finally, a linear interpolation in time and log cubic spline interpolation in energy are used to determine the spectral number luminosity at an arbitrary time and energy.%
\footnote{This approach closely follows the one used in code provided by Totani. There is excellent agreement of the calculated fluxes between Totani’s code and sntools, with differences of at most a few per mille due to slight differences in the numerical interpolation algorithms used.}

\subsection{pre-supernova Time Binning}\label{sec:pre-supernova-time-binning}
As mentioned in section~\ref{sec:more-details}, when using pre-supernova models, the size of the time bins used for calculations within SNTools can be specified through the optional argument \texttt{--binsize <value>}.
During the development of support for pre-supernova models, a study was performed to optimse the time bin size for Hyper-Kamiokande. The study concluded that a bin size of \SI{1}{s} was the best choice,  
limiting any artificial step features in the distribution of events arising from the nature of the time binning. As a result, \SI{1}{s} is chosen as the default value of this argument.
It is advisable to use the default time bin size when working with pre-supernova models unless a detailed time bin optimsation study has been performed for the detector geometry you are using.  

\section{Getting Help and Contributing}

To report problems or ask for help, \href{https://github.com/SNEWS2/sntools/issues}{open an issue on GitHub} or email the lead developer (see title page of this document).

sntools is designed to be versatile and support many different detectors and input formats.
\textbf{Your contributions are welcome!}
To help extend sntools, follow these steps:

\begin{itemize}
	\item Clone the GitHub repository to your local computer.
	\item Edit the code---see the subsections below for help with the most common ways to extend sntools.
	\item Use \texttt{pip install .} to install the modified code; then try it out.
	\item To test your changes, run \texttt{python -m unittest discover} from the top-level directory of the repository. (Tests will also run automatically whenever you submit a pull request on GitHub.)
	\item Finally, please \href{https://docs.github.com/en/github/collaborating-with-issues-and-pull-requests/about-pull-requests}{submit a pull request} with your contributions!
\end{itemize}

sntools is released under a BSD license (see \texttt{LICENSE.md}) and with a Contributor Code of Conduct (see \texttt{CODE\_OF\_CONDUCT.md}).
By participating in this project you agree to abide by its terms.


\subsection{Add a New Detector Configuration}\label{sec:new-detector-config}

In \texttt{detectors.py}:
\begin{itemize}
\item Add the name of your detector configuration to the list \texttt{supported\_detectors}.
\item Add an \texttt{elif} clause for your detector in the \texttt{\_\_init\_\_} function of the \texttt{Detector} class and set the detector dimensions and material there.
\end{itemize}

If your new detector configuration uses a detector material that is not yet implemented, see section~\ref{sec:new-detector-material}.

\subsection{Add a New Detector Material}\label{sec:new-detector-material}

In \texttt{detectors.py}, create a new dictionary for your detector material which contains three entries:
\begin{itemize}
\item \texttt{molecular\_weight}: Molecular weight of the material in \si{g/mol}.
\item \texttt{density}: Density of the material in \si{g/cm^3}.
\item \texttt{channel\_weight}: Dictionary with supported interaction channels as keys and their respective weighting factors as values. The weighting factor is the number of potential targets per molecule. For example, water has 2 protons and 10 electrons per molecule, so the weighting factor is 2 for \texttt{ibd} and 10 for \texttt{es}.
\end{itemize}
% In more advanced cases, you can also write a function, which generates and returns such a dictionary. -> See implementation of WbLS, which lets the user give the fraction of LS as an argument.

If necessary, add new interaction channels which may occur in this detector material (see section~\ref{sec:new-interaction-channel}).
Finally, add a new detector configuration which uses this detector material (see section~\ref{sec:new-detector-config}).

\subsection{Add a New Interaction Channel}\label{sec:new-interaction-channel}

\begin{itemize}
\item In the \texttt{interaction\_channels/} folder, create a new file for your interaction channel. It must define two things:
	\begin{itemize}
		\item A list of neutrino flavors that can interact in this channel (e.g. if all flavors interact in your channel: \texttt{possible\_flavors = ["e", "eb", "x", "xb"]})
		\item A class named \texttt{Channel}, which inherits from the \texttt{BaseChannel} class defined in \texttt{interaction\_channels/\_\_init\_\_.py}. Overwrite all functions and properties that are marked \texttt{@abstractmethod} in the \texttt{BaseChannel} class to implement the relevant physics.
		(See docstrings for explanations.)
	\end{itemize}
\item In \texttt{genevts.py}, add the channel in the \texttt{parse\_command\_line\_options()} function.
\item In \texttt{detectors.py}, add the channel to all relevant detector materials.
\end{itemize}


\subsection{Add a New Input Format}\label{sec:new-input-format}
First, check if it is straightforward to transform the input files into one of the supported formats described in section~\ref{sec:input-formats}.
If not:
\begin{itemize}
\item In the \texttt{formats/} folder, create a new file for your interaction channel. It must contain a class named \texttt{Flux}, which inherits from the \texttt{BaseFlux} class defined in \texttt{formats/\_\_init\_\_.py}. Overwrite all functions that are marked \texttt{@abstractmethod} in the \texttt{BaseFlux} class.
(See docstrings for explanations.)
\item In \texttt{genevts.py}, add the format in the \texttt{parse\_command\_line\_options()} function.
\end{itemize}

%\appendix
%\include{Appendix} % add copy of license here?

\clearpage
\footnotesize % use smaller font size so bibliography doesn’t span as many pages
\bibliographystyle{unsrt}
\bibliography{documentation}

\end{document}