Remove Cherenkov threshold for LS detectors

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

@@ -31,11 +31,7 @@ jobs:
     - name: Run unit tests
       run: python -m unittest discover
     - name: Run integration test
-      run: |
-        # Run setup function
-        python -c 'import sntools; sntools.setup()'
-        # Compare with sample output
-        diff outfile.kin tests/sample_outfile.kin
+      run: python -c 'import sntools; sntools.setup()'
 
   lint:
     runs-on: 'ubuntu-latest'

src/sntools/__init__.py

@@ -58,7 +58,7 @@ def setup():
     with open('outfile.kin', 'r') as f:
         output_sha = hashlib.sha256(f.read().encode('utf-8')).hexdigest()
 
-    test_sha = "0e3444d5a174ee49bd3adad498a2a022d2340736f2e056bfab12df0e5a9b0c54"
+    test_sha = "6a83b302f25e4693bcfa9047876af5a043360b43af0459c26952f78a0127371e"
     if output_sha == test_sha:
         tryprint(u"\u2705", "[SUCCESS]")
         print("Everything seems to work fine. Enjoy using sntools!")
@@ -69,6 +69,7 @@ def setup():
         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.")
+        sys.exit(1)
 
 
 def tryprint(default, alternative=''):

src/sntools/channel.py

@@ -1,33 +1,44 @@
 #!/usr/bin/python
 
+from importlib import import_module
 from math import pi, sin, cos, acos
 import numpy as np
 import random
 from scipy import integrate, interpolate
 
 
-def gen_evts(_channel, _flux, mode, binsize, n_targets, seed, verbose):
+def gen_evts(_channel, flv, _flux, mode, binsize, detector, seed, verbose):
     """Generate events.
 
     * Get event rate by interpolating from time steps in the input data.
     * For each 1ms bin, get number of events from a Poisson distribution.
     * Generate these events from time-dependent energy & direction distribution.
 
     Arguments:
-    _channel -- BaseChannel instance for the current interaction channel
+    _channel -- name of the current interaction channel
+    flv -- incoming neutrino flavor
     _flux -- BaseFlux instance with appropriate flavor and time range (includes weighting due to flux transformation and distance)
-    n_targets -- number of target particles in detector
+    mode -- "ccsn" or "presn"
+    binsize -- bin size in seconds for presupernova mode (ignored for CCSN mode)
+    detector -- Detector instance
     seed -- random number seed to reproducibly generate events
     """
     random.seed(seed)
     np.random.seed(int(seed))
 
+    if "c12e" in detector.material["channel_weights"]:
+        # Detector contains liquid scintillator, so Cherenkov threshold does not apply.
+        # Set threshold slightly above electron mass to avoid ZeroDivisionError in `es` channel.
+        import sntools.interaction_channels
+        sntools.interaction_channels.cherenkov_threshold = 0.512
+
     global channel, cached_flux, flux
     flux = _flux
-    channel = _channel
-    tag = str(channel.__class__).split('.')[-2]
+    channel = import_module("sntools.interaction_channels." + _channel).Channel(flv)
+
+    tag = _channel
     if tag in ('c12nc', 'es'):
-        tag += '-' + str(channel).split("'")[-2]
+        tag += '-' + flv
 
     # ddEventRate(eE, eNu, time) is called hundreds of times for each generated event,
     # often with identical eNu and time values (when integrating over eE).
@@ -37,6 +48,7 @@ def gen_evts(_channel, _flux, mode, binsize, n_targets, seed, verbose):
     # integrate over eE and then eNu to obtain the event rate at time t
     if verbose:
         print(f"[{tag}] Calculating event rate for {flux} ...")
+    n_targets = detector.n_molecules * detector.material["channel_weights"][_channel]
     raw_nevts = [n_targets * integrate.nquad(ddEventRate, [channel.bounds_eE, channel.bounds_eNu], args=[t], opts=[channel._opts, {}])[0]
                  for t in flux.raw_times]
     event_rate = interpolate.pchip(flux.raw_times, raw_nevts)

src/sntools/genevts.py

@@ -4,6 +4,7 @@
 from concurrent.futures import ProcessPoolExecutor, as_completed
 from datetime import datetime
 from importlib import import_module
+import multiprocessing as mp
 import random
 import uproot
 
@@ -36,16 +37,14 @@ def main():
 
     # using the process pool executor for ccsn mode 
     if args.mode == "ccsn":
-        pool = ProcessPoolExecutor(max_workers=args.maxworkers)
+        pool = ProcessPoolExecutor(max_workers=args.maxworkers, mp_context=mp.get_context("spawn"))
         results = []
 
         for channel in sorted(args.channels):
             mod_channel = import_module("sntools.interaction_channels." + channel)
-            n_targets = args.detector.n_molecules * args.detector.material["channel_weights"][channel]
             for flv in mod_channel.possible_flavors:
-                channel_instance = mod_channel.Channel(flv)
                 for flux in flux_at_detector.components[flv]:
-                    results.append(pool.submit(gen_evts, channel_instance, flux, args.mode, args.binsize, n_targets, args.randomseed + random.random(), args.verbose))
+                    results.append(pool.submit(gen_evts, channel, flv, flux, args.mode, args.binsize, args.detector, args.randomseed + random.random(), args.verbose))
 
         events = []
         for result in as_completed(results):
@@ -58,11 +57,9 @@ def main():
 
         for channel in sorted(args.channels):
             mod_channel = import_module("sntools.interaction_channels." + channel)
-            n_targets = args.detector.n_molecules * args.detector.material["channel_weights"][channel]
             for flv in mod_channel.possible_flavors:
-                channel_instance = mod_channel.Channel(flv)
                 for flux in flux_at_detector.components[flv]:
-                    events.extend(gen_evts(_channel=channel_instance, _flux=flux, mode=args.mode, binsize=args.binsize, n_targets=n_targets, seed=args.randomseed + random.random(), verbose=args.verbose))
+                    events.extend(gen_evts(channel, flv, flux, args.mode, args.binsize, args.detector, args.randomseed + random.random(), args.verbose))
 
 
 

src/sntools/interaction_channels/__init__.py

@@ -4,6 +4,7 @@
 
 
 cherenkov_threshold = 0.77  # Cherenkov threshold of electron/positron in water
+# If a detector contains liquid scintillator, we modify this threshold in gen_evts() in channel.py
 
 
 class BaseChannel(ABC):