Merge branch 'add_time_converter' into 'master'

Add `TimeConverter` and `has_(jmuon|jshower|aashower|dusjshower)`

See merge request !74

CHANGELOG.rst

@@ -3,6 +3,13 @@ Unreleased changes
 
 Version 0
 ---------
+0.27.0 / 2022-07-20
+~~~~~~~~~~~~~~~~~~~
+* Adds the TimeConverter class to ``src/km3io/tools.py``
+* Update km3net-dataformat requirement to version 0.3.6 or higher
+* Update Black requirement to version 22.3.0 or higher, to prevent ``ImportError: cannot import name '_unicodefun' from 'click'``
+* Remove ``requirements`` folder (all requirements are now configured in ``setup.cfg``)
+
 0.26.1 / 2022-07-06
 ~~~~~~~~~~~~~~~~~~~
 * The warning from OpenMP/Numba is now silenced
@@ -12,7 +19,7 @@ Version 0
 * Added ``km3io.tools.is_nanobeacon()`` to check if the nanobeacon trigger bit is set
 * Added ``km3io.tools.get_w2list_idx()`` to get the w2list index according to the
   simulation program
-
+  
 0.25.2 / 2022-03-27
 ~~~~~~~~~~~~~~~~~~~
 * Fixes the version

requirements/dev.txt

-black==21.6b0
-km3net-testdata>=0.2.26
-ipykernel
-matplotlib
-memory_profiler
-numpydoc==0.9.2
-pillow
-pytest
-pytest-cov
-pytest-flake8
-pytest-pylint
-pytest-watch
-scipy
-sphinx
-sphinx-autoapi
-sphinx-gallery>=0.1.12
-sphinx_rtd_theme
-sphinxcontrib-versioning
-wheel

requirements/install.txt

-docopt
-numba>=0.50
-awkward>=1.0.0rc2
-awkward0
-uproot3>=3.11.1
-uproot>=4.2.2
-setuptools_scm

setup.cfg

@@ -53,8 +53,8 @@ where = src
 [options.extras_require]
 all =
 dev =
-    black==21.6b0
-    km3net-testdata>=0.3.3
+    black>=22.3.0
+    km3net-testdata>=0.3.6
     ipykernel
     matplotlib
     memory_profiler

src/km3io/online.py

@@ -7,9 +7,9 @@ import numpy as np
 
 import numba as nb
 
-TIMESLICE_FRAME_BASKET_CACHE_SIZE = 523 * 1024 ** 2  # [byte]
-SUMMARYSLICE_FRAME_BASKET_CACHE_SIZE = 523 * 1024 ** 2  # [byte]
-BASKET_CACHE_SIZE = 110 * 1024 ** 2
+TIMESLICE_FRAME_BASKET_CACHE_SIZE = 523 * 1024**2  # [byte]
+SUMMARYSLICE_FRAME_BASKET_CACHE_SIZE = 523 * 1024**2  # [byte]
+BASKET_CACHE_SIZE = 110 * 1024**2
 BASKET_CACHE = uproot3.cache.ThreadSafeArrayCache(BASKET_CACHE_SIZE)
 
 # Parameters for PMT rate conversions, since the rates in summary slices are

src/km3io/tools.py

@@ -13,7 +13,7 @@ from km3io.definitions import w2list_genhen as kw2gen
 from km3io.definitions import w2list_gseagen as kw2gsg
 
 # 110 MB based on the size of the largest basket found so far in km3net
-BASKET_CACHE_SIZE = 110 * 1024 ** 2
+BASKET_CACHE_SIZE = 110 * 1024**2
 BASKET_CACHE = uproot3.cache.ThreadSafeArrayCache(BASKET_CACHE_SIZE)
 
 
@@ -404,6 +404,30 @@ def _mask_atleast(arr, atleast):
     return out
 
 
+def has_jmuon(tracks):
+    """Check if given tracks contain JMUON reconstruction."""
+    m = mask(tracks.rec_stages, minmax=(krec.JMUONBEGIN, krec.JMUONEND))
+    return ak.any(m, axis=m.ndim - 1)
+
+
+def has_jshower(tracks):
+    """Check if given tracks contain JSHOWER reconstruction."""
+    m = mask(tracks.rec_stages, minmax=(krec.JSHOWERBEGIN, krec.JSHOWEREND))
+    return ak.any(m, axis=m.ndim - 1)
+
+
+def has_aashower(tracks):
+    """Check if given tracks contain AASHOWER reconstruction."""
+    m = mask(tracks.rec_stages, minmax=(krec.AASHOWERBEGIN, krec.AASHOWEREND))
+    return ak.any(m, axis=m.ndim - 1)
+
+
+def has_dusjshower(tracks):
+    """Check if given tracks contain AASHOWER reconstruction."""
+    m = mask(tracks.rec_stages, minmax=(krec.DUSJSHOWERBEGIN, krec.DUSJSHOWEREND))
+    return ak.any(m, axis=m.ndim - 1)
+
+
 def best_jmuon(tracks):
     """Select the best JMUON track."""
     return best_track(tracks, minmax=(krec.JMUONBEGIN, krec.JMUONEND))
@@ -566,3 +590,53 @@ def is_nanobeacon(trigger_mask):
       A value or an array of the trigger_mask, either of an event, or a hit.
     """
     return is_bit_set(trigger_mask, ktrg.JTRIGGERNB)
+
+
+class TimeConverter(object):
+    """
+    Auxiliary class to convert Monte Carlo hit times to DAQ/triggered hit times.
+    """
+
+    FRAME_TIME_NS = 1e8  # [ns]
+
+    def __init__(self, event):
+
+        self.__t0 = event.mc_t  # [ns]
+        self.__t1 = self.get_time_of_frame(event.frame_index)  # [ns]
+
+    def get_time_of_frame(self, frame_index):
+        """
+        Get start time of frame in ns since start of run for a given frame index
+
+        Parameters
+        ----------
+        frame_index: int
+          The index of the DAQ frame
+        """
+
+        if frame_index > 0:
+            return (frame_index - 1) * self.FRAME_TIME_NS  # [ns]
+        else:
+            return 0  # [ns]
+
+    def get_DAQ_time(self, t0):
+        """
+        Get DAQ/triggered hit time
+
+        Parameters
+        ----------
+        t0: float or array(float)
+          Simulated time [ns]
+        """
+        return t0 + (self.__t0 - self.__t1)  # [ns]
+
+    def get_MC_time(self, t0):
+        """
+        Get Monte Carlo hit time
+
+        Parameters
+        ----------
+        t0: float or array(float)
+          DAQ/trigger hit time [ns]
+        """
+        return t0 - (self.__t0 - self.__t1)  # [ns]

tests/test_tools.py

@@ -21,6 +21,10 @@ from km3io.tools import (
     best_track,
     get_w2list_param,
     get_multiplicity,
+    has_jmuon,
+    has_jshower,
+    has_aashower,
+    has_dusjshower,
     best_jmuon,
     best_jshower,
     best_aashower,
@@ -32,8 +36,10 @@ from km3io.tools import (
     is_mxshower,
     is_3dmuon,
     is_nanobeacon,
+    TimeConverter,
 )
 
+
 OFFLINE_FILE = OfflineReader(data_path("offline/km3net_offline.root"))
 OFFLINE_USR = OfflineReader(data_path("offline/usr-sample.root"))
 OFFLINE_MC_TRACK_USR = OfflineReader(
@@ -41,6 +47,11 @@ OFFLINE_MC_TRACK_USR = OfflineReader(
         "offline/mcv5.11r2.gsg_muonCChigherE-CC_50-5000GeV.km3_AAv1.jterbr00004695.jchain.aanet.498.root"
     )
 )
+OFFLINE_JMERGEFIT = OfflineReader(
+    data_path(
+        "offline/mcv5.0.gsg_elec-CC_10-100GeV.km3sim.JDK.jte.jmergefit.orca.aanet.909.evtsample.root"
+    )
+)
 GENHEN_OFFLINE_FILE = OfflineReader(
     data_path("offline/mcv5.1.genhen_anumuNC.sirene.jte.jchain.aashower.sample.root")
 )
@@ -286,16 +297,44 @@ class TestBestTrackSelection(unittest.TestCase):
             best_track(self.events.tracks, startend=(1, 4), stages=[1, 3, 5, 4])
 
 
+class TestHasJmuon(unittest.TestCase):
+    def test_has_jmuon(self):
+        assert ak.sum(has_jmuon(OFFLINE_JMERGEFIT.events.tracks)) == len(
+            OFFLINE_JMERGEFIT.events.tracks
+        )
+
+
+class TestHasJshower(unittest.TestCase):
+    def test_has_jshower(self):
+        assert ak.sum(has_jshower(OFFLINE_JMERGEFIT.events.tracks)) == len(
+            OFFLINE_JMERGEFIT.events.tracks
+        )
+
+
+class TestHasAashower(unittest.TestCase):
+    def test_has_aashower(self):
+        # there are no aashower events in this file
+        assert ak.sum(has_aashower(OFFLINE_JMERGEFIT.events.tracks)) == 0
+
+
+class TestHasDusjshower(unittest.TestCase):
+    def test_has_dusjshower(self):
+        # there are no dusj events in this file
+        assert ak.sum(has_dusjshower(OFFLINE_JMERGEFIT.events.tracks)) == 0
+
+
 class TestBestJmuon(unittest.TestCase):
     def test_best_jmuon(self):
         best = best_jmuon(OFFLINE_FILE.events.tracks)
 
         assert len(best) == 10
 
-        assert best.rec_stages[0].tolist() == [1, 3, 5, 4]
-        assert best.rec_stages[1].tolist() == [1, 3, 5, 4]
-        assert best.rec_stages[2].tolist() == [1, 3, 5, 4]
-        assert best.rec_stages[3].tolist() == [1, 3, 5, 4]
+        jmuon_stages = [1, 3, 5, 4]
+
+        assert best.rec_stages[0].tolist() == jmuon_stages
+        assert best.rec_stages[1].tolist() == jmuon_stages
+        assert best.rec_stages[2].tolist() == jmuon_stages
+        assert best.rec_stages[3].tolist() == jmuon_stages
 
         assert best.lik[0] == ak.max(OFFLINE_FILE.events.tracks.lik[0])
         assert best.lik[1] == ak.max(OFFLINE_FILE.events.tracks.lik[1])
@@ -304,19 +343,25 @@ class TestBestJmuon(unittest.TestCase):
 
 class TestBestJshower(unittest.TestCase):
     def test_best_jshower(self):
-        # there are no jshower events in this file
-        best = best_jshower(OFFLINE_FILE.events.tracks)
+        best = best_jshower(OFFLINE_JMERGEFIT.events.tracks)
 
         assert len(best) == 10
 
-        assert best.rec_stages[0] is None
-        assert best.rec_stages[1] is None
-        assert best.rec_stages[2] is None
-        assert best.rec_stages[3] is None
+        jshower_stages = [101, 106, 102, 105, 107, 103]
 
-        assert best.lik[0] is None
-        assert best.lik[1] is None
-        assert best.lik[2] is None
+        assert best.rec_stages[0].tolist() == jshower_stages
+        assert best.rec_stages[1].tolist() == jshower_stages
+        assert best.rec_stages[2].tolist() == jshower_stages
+        assert best.rec_stages[3].tolist() == jshower_stages
+
+        jshower_mask = mask(
+            OFFLINE_JMERGEFIT.events.tracks.rec_stages, sequence=jshower_stages
+        )
+        jshower_tracks = OFFLINE_JMERGEFIT.events.tracks[jshower_mask]
+
+        assert best.lik[0] == ak.max(jshower_tracks.lik[0])
+        assert best.lik[1] == ak.max(jshower_tracks.lik[1])
+        assert best.lik[2] == ak.max(jshower_tracks.lik[2])
 
 
 class TestBestAashower(unittest.TestCase):
@@ -338,7 +383,7 @@ class TestBestAashower(unittest.TestCase):
 
 class TestBestDusjshower(unittest.TestCase):
     def test_best_dusjshower(self):
-        # there are no aashower events in this file
+        # there are no dusj events in this file
         best = best_dusjshower(OFFLINE_FILE.events.tracks)
 
         assert len(best) == 10
@@ -653,3 +698,29 @@ class TestTriggerMaskChecks(unittest.TestCase):
             [False, False, False, False, False, False, False, False, False, False],
             list(is_nanobeacon(GENHEN_OFFLINE_FILE.events.trigger_mask)),
         )
+
+
+class TestTimeConverter(unittest.TestCase):
+    def setUp(self):
+        self.one_event = GENHEN_OFFLINE_FILE.events[5]
+        self.tconverter = TimeConverter(self.one_event)
+
+    def test_get_time_of_frame(self):
+        t_frame = self.tconverter.get_time_of_frame(self.one_event.frame_index)
+
+        assert t_frame < self.one_event.mc_t
+        assert t_frame % self.tconverter.FRAME_TIME_NS == 0
+
+    def test_get_DAQ_time(self):
+        t_DAQ = self.tconverter.get_DAQ_time(self.one_event.mc_hits.t)
+
+        assert all(t_DAQ < self.tconverter.FRAME_TIME_NS)
+
+    def test_get_MC_time(self):
+        t_frame = self.tconverter.get_time_of_frame(self.one_event.frame_index)
+        t_MC = self.tconverter.get_MC_time(self.one_event.hits.t)
+
+        assert all(
+            np.fabs(self.one_event.mc_t + t_MC - t_frame)
+            < self.tconverter.FRAME_TIME_NS
+        )