Pre SN model related changes - doc, unit tests, Pre SN integration test

.github/workflows/test-and-lint.yml

@@ -36,6 +36,19 @@ jobs:
         python -c 'import sntools; sntools.setup()'
         # Compare with sample output
         diff outfile.kin tests/sample_outfile.kin
+    - name: Run pre supernova integration test
+      run: |
+        # Run presnsetup function
+        python -c 'import sntools; sntools.presnsetup()'
+        # compare with sample output
+        diff presnoutfile.kin tests/sample_presnoutfile.kin
+    - name: Run pre supernova binsize test
+      run: |
+        # Run presn binszie test
+        python -c 'import sntools; sntools.presnbinsizetest()'
+        # compare with sample output
+        diff presnbinsizeoutfile.kin tests/sample_presnbinsizeoutfile.kin
+
 
   lint:
     runs-on: 'ubuntu-latest'

README.md

@@ -1,8 +1,8 @@
 # sntools
-sntools is a Monte Carlo event generator for supernova neutrino interactions.
+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 and later extended to support different detectors and detector materials.
+sntools was originally developed for Hyper-Kamiokande 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.
 

doc/documentation.tex

@@ -35,7 +35,7 @@
 \newcommand{\nuxbar}{\ensuremath{\bar{\nu}_x}\xspace}
 
 
-\title{Documentation for sntools 1.5b1\footnote{See \url{https://github.com/SNEWS2/sntools} for the most recent version.}}
+\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}
 
@@ -48,10 +48,10 @@
 
 \section{Introduction}
 
-sntools is a Monte Carlo event generator for supernova neutrino interactions.
+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} and later extended to support different detectors and detector materials.
+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}.
 
@@ -93,7 +93,8 @@ \subsection{Quick Start Guide}
 \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 (in milliseconds). If these aren’t specified, events are generated for the full time range for which fluxes are given in the input file.
+\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.
@@ -104,7 +105,7 @@ \subsection{Quick Start Guide}
 
 
 
-\section{More Details}
+\section{More Details} \label{sec:more-details}
 
 \begin{figure}[htbp]
 	\centering
@@ -120,7 +121,9 @@ \section{More Details}
 
 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 \SI{1}{ms} bin, which is given by
+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}
@@ -323,7 +326,12 @@ \subsection{Input Formats} \label{sec:input-formats}
 \texttt{Tamborra\_2014},
 \texttt{Walk\_2018},
 \texttt{Walk\_2019} and
-\texttt{Zha\_2021}.
+\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}.
 
@@ -396,7 +404,11 @@ \subsubsection{Totani Format}
 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}
 

pyproject.toml

@@ -10,7 +10,7 @@ authors = [
     {name="Jost Migenda"},
     {name="sntools Contributors"},
 ]
-description = "Event generator for supernova burst neutrinos"
+description = "Event generator for supernova burst and pre-supernova neutrinos"
 readme = "README.md"
 license = {file = "LICENSE"}
 
@@ -23,7 +23,7 @@ dependencies = [
     "uproot"
 ]
 
-keywords = ["astrophysics", "supernova", "neutrino", "event generator"]
+keywords = ["astrophysics", "supernova", "pre-supernova", "neutrino", "event generator"]
 classifiers = [
     'License :: OSI Approved :: BSD License',
     'Development Status :: 5 - Production/Stable',

src/sntools/__init__.py

@@ -1,12 +1,12 @@
 """An Event Generator for Supernova Burst Neutrinos
 
-sntools is a Monte Carlo event generator for supernova neutrino interactions.
+sntools is a Monte Carlo event generator for supernova and pre-supernova neutrino interactions.
 Based on detailed time- and energy-dependent neutrino fluxes provided by
 various supernova models, it generates interactions within the detector volume
 and writes them to event files that can be used as an input for a full detector
 simulation.
-sntools was originally developed for Hyper-Kamiokande and later extended to
-support different detectors and detector materials.
+sntools was originally developed for Hyper-Kamiokande using core collapse supernova models and later extended to 
+support different detectors, detector materials and to use pre supernova models.
 
 For usage information, run `sntools -h` or `python sntools/genevts.py -h`.
 For more extensive documentation, to report issues or to contribute code,
@@ -70,6 +70,116 @@ def setup():
         tryprint(u"\u274c", "[ERROR]")
         print("If this persists, please go to https://github.com/SNEWS2/sntools and open a new issue.")
 
+def presnsetup():
+    """
+    Downloads sample flux file from GitHub if necessary and performs pre supernova integration test.
+    """
+    tryprint(u"\u2705", "[SUCCESS]")
+    print("sntools was imported from " + __path__[0])
+    import hashlib
+    import os
+    import sys
+    from . import genevts
+
+    flux_dir = 'fluxes/'
+    flux_file = flux_dir + 'totalLuminosity_15SolarMass.dat'
+    flux_url = 'https://raw.githubusercontent.com/SNEWS2/snewpy-models-presn/master/models/Patton_2017/totalLuminosity_15SolarMass.dat'
+    if os.path.exists(flux_file):
+        tryprint(u"\u2705", "[SUCCESS]")
+        print("Using sample flux file at " + flux_file)
+    else:
+        tryprint(u"\U0001f6e0")
+        print("Downloading sample flux file from " + flux_url)
+        if not os.path.isdir(flux_dir):
+            os.mkdir(flux_dir)
+
+        from urllib.request import urlretrieve
+        try:
+            urlretrieve(flux_url, filename=flux_file)
+            tryprint(u"\u2705", "[SUCCESS]")
+            print("Saved sample flux file to " + flux_file)
+        except IOError:
+            tryprint(u"\u274c", "[ERROR]")
+            print("Cannot download sample flux file.")
+            sys.exit(-1)
+
+    tryprint(u"\U0001f6e0")
+    print("Testing event generation ...")
+    sys.argv += [flux_file, '--format', 'SNEWPY-Patton_2017', '--detector', 'HyperK', '--distance', '0.15', '--transformation', 'AdiabaticMSW_IMO', '--starttime', '-120', '--endtime', '0', '--binsize', '1', '-o', 'presnoutfile.kin', '--randomseed', '205']
+    genevts.main()
+
+    tryprint(u"\U0001f6e0")
+    print("Checking output file ...")
+    with open('presnoutfile.kin', 'r') as f:
+        output_sha = hashlib.sha256(f.read().encode('utf-8')).hexdigest()
+
+    test_sha = "824dffa7ed01016a9044a28b36ca33aacdd81a9db662ca69a08d9ab0eb32bf29"
+    if output_sha == test_sha:
+        tryprint(u"\u2705", "[SUCCESS]")
+        print("Everything seems to work fine. Enjoy using sntools!")
+    else:
+        tryprint(u"\u274c", "[ERROR]")
+        print("Test did not generate the expected events.")
+        tryprint(u"\u274c", "[ERROR]")
+        print("Please ensure you have installed the most recent version of sntools and all dependencies.")
+        tryprint(u"\u274c", "[ERROR]")
+        print("If this persists, please go to https://github.com/SNEWS2/sntools and open a new issue.")
+
+def presnbinsizetest():
+    """
+    Downloads sample flux file from GitHub if necessary and performs pre supernova binsize test.
+    """
+    tryprint(u"\u2705", "[SUCCESS]")
+    print("sntools was imported from " + __path__[0])
+    import hashlib
+    import os
+    import sys
+    from . import genevts
+
+    flux_dir = 'fluxes/'
+    flux_file = flux_dir + 'totalLuminosity_15SolarMass.dat'
+    flux_url = 'https://raw.githubusercontent.com/SNEWS2/snewpy-models-presn/master/models/Patton_2017/totalLuminosity_15SolarMass.dat'
+    if os.path.exists(flux_file):
+        tryprint(u"\u2705", "[SUCCESS]")
+        print("Using sample flux file at " + flux_file)
+    else:
+        tryprint(u"\U0001f6e0")
+        print("Downloading sample flux file from " + flux_url)
+        if not os.path.isdir(flux_dir):
+            os.mkdir(flux_dir)
+
+        from urllib.request import urlretrieve
+        try:
+            urlretrieve(flux_url, filename=flux_file)
+            tryprint(u"\u2705", "[SUCCESS]")
+            print("Saved sample flux file to " + flux_file)
+        except IOError:
+            tryprint(u"\u274c", "[ERROR]")
+            print("Cannot download sample flux file.")
+            sys.exit(-1)
+
+    tryprint(u"\U0001f6e0")
+    print("Testing event generation ...")
+    sys.argv += [flux_file, '--format', 'SNEWPY-Patton_2017', '--detector', 'HyperK', '--distance', '0.15', '--transformation', 'AdiabaticMSW_NMO', '--starttime', '-120', '--endtime', '0', '--binsize', '10', '-o', 'presnbinsizeoutfile.kin', '--randomseed', '100']
+    genevts.main()
+
+    tryprint(u"\U0001f6e0")
+    print("Checking output file ...")
+    with open('presnbinsizeoutfile.kin', 'r') as f:
+        output_sha = hashlib.sha256(f.read().encode('utf-8')).hexdigest()
+
+    test_sha = "7e1e8d14c732d171a9ad9f924bde4ca31f6897ca24aa0c94074e008fd69e76e7"
+    if output_sha == test_sha:
+        tryprint(u"\u2705", "[SUCCESS]")
+        print("Everything seems to work fine. Enjoy using sntools!")
+    else:
+        tryprint(u"\u274c", "[ERROR]")
+        print("Test did not generate the expected events.")
+        tryprint(u"\u274c", "[ERROR]")
+        print("Please ensure you have installed the most recent version of sntools and all dependencies.")
+        tryprint(u"\u274c", "[ERROR]")
+        print("If this persists, please go to https://github.com/SNEWS2/sntools and open a new issue.")
+
 
 def tryprint(default, alternative=''):
     try:

src/sntools/channel.py

@@ -6,7 +6,7 @@
 from scipy import integrate, interpolate
 
 
-def gen_evts(_channel, _flux, mode, binsize, n_targets, seed, verbose):
+def gen_evts(_channel, _flux, format, binsize, n_targets, seed, verbose):
     """Generate events.
 
     * Get event rate by interpolating from time steps in the input data.
@@ -41,16 +41,16 @@ def gen_evts(_channel, _flux, mode, binsize, n_targets, seed, verbose):
                  for t in flux.raw_times]
     event_rate = interpolate.pchip(flux.raw_times, raw_nevts)
 
-    # appropriate bin width is different for each mode
-    if mode == "ccsn":
-        bin_width = 1  # in ms
-    elif mode == "presn":
+    # appropriate bin width is different for ccsn and presn models
+    if format in ("Odrzywolek_2010", "Yoshida_2016", "Patton_2017", "Kato_2017"):
         bin_width = binsize*1000 # converting seconds into ms
+    else:
+        bin_width = 1  # in ms
 
 
     n_bins = int((flux.endtime - flux.starttime) / bin_width)  # number of full-width bins; int() implies floor()
     if verbose:
-        if mode == "presn":
+        if format in ("Odrzywolek_2010", "Yoshida_2016", "Patton_2017", "Kato_2017"):
             print(f"[{tag}] Generating events in {binsize} second bins from {(flux.starttime)/1000} to {(flux.endtime)/1000} seconds ...")
         else:
             print(f"[{tag}] Generating events in {bin_width} ms bins from {flux.starttime} to {flux.endtime} ms ...")

src/sntools/formats/__init__.py

@@ -156,19 +156,20 @@ class SNEWPYCompositeFlux(CompositeFlux):
     """Adapter class to turn a SNEWPY.models.SupernovaModel into an sntools.formats.CompositeFlux"""
 
     @classmethod
-    def from_file(cls, file, mode, format, starttime=None, endtime=None):
+    def from_file(cls, file, format, starttime=None, endtime=None):
         """Create a SNEWPYCompositeFlux from an input file."""
         self = SNEWPYCompositeFlux()
         self._repr = f"SNEWPYCompositeFlux.from_file('{file}', format='{format}', starttime={starttime}, endtime={endtime})"
 
         # snewpy.models.loaders classes treat relative file paths as relative to the snewpy cache directory,
         # so we’ll turn it into an absolute path first.
 
-        # need an if statement here to import the snewpy.models.presn_loaders module if the mode is set to presn.
-        if mode == "ccsn":
-            sn_model = getattr(import_module('snewpy.models.ccsn_loaders'), format)(abspath(file))
-        elif mode == "presn":
+        # need an if statement here to import the snewpy.models.presn_loaders module if the model is presn
+        if format in ("Odrzywolek_2010", "Yoshida_2016", "Patton_2017", "Kato_2017"):
             sn_model = getattr(import_module('snewpy.models.presn_loaders'), format)(abspath(file))
+        else:
+            sn_model = getattr(import_module('snewpy.models.ccsn_loaders'), format)(abspath(file))
+
 
 
         for flv in ('e', 'eb', 'x', 'xb'):