Filter nexus delayed hits in detsim

invisible_cities/cities/detsim.py

@@ -14,6 +14,7 @@
 """
 
 import os
+import warnings
 import numpy  as np
 import tables as tb
 
@@ -44,6 +45,21 @@
 from .. detsim.detsim_waveforms   import s1_waveforms_creator
 
 
+@check_annotations
+def filter_hits_after_max_time(max_time : float):
+    """
+    Function that filters and warns about delayed hits
+    (hits at times larger that max_time configuration parameter)
+    """
+    def select_hits(x, y, z, energy, time, label, event_number):
+        sel = ((time - min(time)) < max_time)
+        if sel.all(): return x, y, z, energy, time, label
+        else:
+            warnings.warn(f"Delayed hits at event {event_number}")
+            return x[sel], y[sel], z[sel], energy[sel], time[sel], label[sel]
+    return select_hits
+
+
 @check_annotations
 def hits_selector(active_only: bool=True):
     """
@@ -157,6 +173,10 @@ def detsim( *
     lt_sipm  = LT_SiPM(fname=os.path.expandvars(sipm_psf), sipm_database=datasipm)
     el_gap   = lt_sipm.el_gap_width
 
+    filter_delayed_hits = fl.map(filter_hits_after_max_time(buffer_params_["max_time"]),
+                                 args = ('x', 'y', 'z', 'energy', 'time', 'label', 'event_number'),
+                                 out  = ('x', 'y', 'z', 'energy', 'time', 'label'))
+
     select_s1_candidate_hits = fl.map(hits_selector(False),
                                 item = ('x', 'y', 'z', 'energy', 'time', 'label'))
 
@@ -173,6 +193,11 @@ def detsim( *
                                 args = ('x_a', 'y_a', 'z_a', 'time_a', 'energy_a'),
                                 out  = ('x_ph', 'y_ph', 'z_ph', 'times_ph', 'nphotons'))
 
+    count_photons = fl.map(lambda x: np.sum(x) > 0,
+                           args= 'nphotons',
+                           out = 'enough_photons')
+    dark_events   = fl.count_filter(bool, args='enough_photons')
+
     get_buffer_info = buffer_times_and_length_getter(buffer_params_["pmt_width"],
                                                      buffer_params_["sipm_width"],
                                                      el_gap, el_dv,
@@ -182,7 +207,7 @@ def detsim( *
                                          out = ('tmin', 'buffer_length'))
 
     create_pmt_s1_waveforms = fl.map(s1_waveforms_creator(s1_lighttable, ws, buffer_params_["pmt_width"]),
-                                     args = ('x_a', 'y_a', 'z_a', 'time_a', 'energy_a', 'tmin', 'buffer_length'),
+                                     args = ('x', 'y', 'z', 'time', 'energy', 'tmin', 'buffer_length'),
                                      out = 's1_pmt_waveforms')
 
     create_pmt_s2_waveforms = fl.map(s2_waveform_creator(buffer_params_["pmt_width"], lt_pmt, el_dv),
@@ -224,16 +249,21 @@ def detsim( *
                                                        , order_sensors = None)
 
         write_nohits_filter   = fl.sink(event_filter_writer(h5out, "active_hits"), args=("event_number", "passed_active"))
+        write_dark_evt_filter = fl.sink(event_filter_writer(h5out, "dark_events"), args=("event_number", "enough_photons"))
         result = fl.push(source= MC_hits_from_files(files_in, rate),
                          pipe  = fl.pipe( fl.slice(*event_range, close_all=True)
                                         , event_count_in.spy
                                         , print_every(print_mod)
+                                        , filter_delayed_hits
                                         , select_s1_candidate_hits
                                         , select_active_hits
                                         , filter_events_no_active_hits
                                         , fl.branch(write_nohits_filter)
                                         , events_passed_active_hits.filter
                                         , simulate_electrons
+                                        , count_photons
+                                        , fl.branch(write_dark_evt_filter)
+                                        , dark_events.filter
                                         , get_buffer_times_and_length
                                         , create_pmt_waveforms
                                         , create_sipm_waveforms
@@ -242,7 +272,8 @@ def detsim( *
                                         , "event_number"
                                         , evtnum_collect.sink),
                          result = dict(events_in     = event_count_in.future,
-                                       evtnum_list   = evtnum_collect.future))
+                                       evtnum_list   = evtnum_collect.future,
+                                       dark_events   = dark_events   .future))
 
         copy_mc_info(files_in, h5out, result.evtnum_list,
                      detector_db, run_number)

invisible_cities/cities/detsim_test.py

@@ -58,6 +58,28 @@ def test_detsim_filter_active_hits(ICDATADIR, output_tmpdir):
         np.testing.assert_array_equal(filters["passed"], [False, True])
 
 
+def test_detsim_filter_dark_events(ICDATADIR, output_tmpdir):
+    # this file contains delayed hits that are filtered out
+    # leaving a very low energy hit that does not produce electrons
+    # this test shows that without filtering out these events the city crashes
+    PATH_IN  = os.path.join(ICDATADIR    , "nexus_not_enough_energy.h5")
+    PATH_OUT = os.path.join(output_tmpdir, "filtered_dark_events.h5")
+    conf = configure('detsim $ICTDIR/invisible_cities/config/detsim.conf'.split())
+    conf.update(dict(files_in    = PATH_IN,
+                     file_out    = PATH_OUT,
+                     run_number  = 0))
+    conf["buffer_params"]["max_time"] = 3 * units.ms
+
+    result = detsim(**conf)
+
+    assert result.events_in   == 1
+    assert result.evtnum_list == []
+
+    with tb.open_file(PATH_OUT, mode="r") as h5out:
+        filters = h5out.root.Filters.dark_events.read()
+        np.testing.assert_array_equal(filters["passed"], [False])
+
+
 def test_detsim_filter_empty_waveforms(ICDATADIR, output_tmpdir):
     # the first event radius is slighty above NEW active radius of 208.0 mm
     PATH_IN  = os.path.join(ICDATADIR, "nexus_new_kr83m_fast.newformat.sim.emptywfs.h5")

invisible_cities/detsim/detsim_waveforms.py

@@ -1,37 +1,51 @@
 import os
 import numpy as np
 
-from typing import Callable
+from typing          import Callable
+from scipy.integrate import quad
 
 from . light_tables import create_lighttable_function
 
 from . simulate_s1 import compute_scintillation_photons
 from . simulate_s1 import compute_s1_pes_at_pmts
-from . simulate_s1 import generate_s1_times_from_pes
+from . simulate_s1 import s1_times_pdf
 
+from ..core.system_of_units import mus
 
-def histogram_s1_times(s1_times      : list,
-                       buffer_length : float,
-                       bin_width     : float,
-                       start_time    : float=0)->np.ndarray:
+def generate_s1_waveform(pes_at_pmts, hit_times, buffer_length, bin_width, start_time):
     """
-    s1_times are histogramed in a waveform.
-    Parameters
-        :s1_times: list
-            output of generate_s1_times_from_pes, list of size equal to the
-            number of pmts. Each element is an array with the times of the S1
+    Creates the s1-waveform based on the pdf at :s1_times_pdf:
+
+    Parameters:
+        :pes_at_pmts: np.ndarray
+            the pes at each PMT generated by each hit
+        :hit_times: np.ndarray
+            emission time of each hit
         :buffer_length: float
-            waveform buffer lenght in time units
-        :bin_width: float
-            waveform bin width in time units
-        :start_time:
-            waveform time of first bin edge in time units
+            waveform buffer_length, waveform will range from tmin to
+            tmin + buffer_length
+        :bin_width:
+            waveform bin width
+        :start_time: float
+            waveform start time
     Returns:
-        :wfs: np.ndarray
-            waveforms with start_time, buffer_length and bin_width
+        :wfs: list[np.ndarray,..]
+            waveform with the S1 signal
     """
     bins = np.arange(start_time, start_time + buffer_length + bin_width, bin_width)
-    wfs  = np.stack([np.histogram(times, bins=bins)[0] for times in s1_times])
+    wfs_ = np.stack([np.histogram(hit_times, bins=bins, weights=pes)[0] for pes in pes_at_pmts])
+
+    # apply s1-time distribution
+    wfs = np.zeros(wfs_.shape)
+    ## discrete normalization (integrate only up to 2 mus to fasten up)
+    ts = np.arange(0, min(2. * mus, buffer_length) + bin_width, bin_width)
+    p  = np.zeros(int(buffer_length/bin_width)+1)
+    for i, (t0, t1) in enumerate(zip(ts[:-1], ts[1:])): p[i] = quad(s1_times_pdf, t0, t1)[0]
+    ## apply over positive counts
+    for pmt, wf in enumerate(wfs_):
+        idxs = np.argwhere(wf>0).flatten()
+        for idx in idxs:
+            wfs[pmt, idx:] += (wf[idx] * p[:-idx-1])
     return wfs
 
 
@@ -68,8 +82,7 @@ def create_s1_waveforms_from_hits(x, y, z, time, energy, tmin, buffer_length):
         """
         s1_photons     = compute_scintillation_photons(energy, ws)
         s1_pes_at_pmts = compute_s1_pes_at_pmts(x, y, z, s1_photons, s1_lt)
-        s1times = generate_s1_times_from_pes(s1_pes_at_pmts, time)
-        s1_wfs  = histogram_s1_times(s1times, buffer_length, wf_pmt_bin_width, tmin)
-        return s1_wfs
+        s1_wfs = generate_s1_waveform(s1_pes_at_pmts, time, buffer_length, wf_pmt_bin_width, tmin)
+        return np.random.poisson(s1_wfs)
 
     return create_s1_waveforms_from_hits

invisible_cities/detsim/detsim_waveforms_test.py

@@ -1,14 +1,31 @@
 import numpy as np
 
-from . detsim_waveforms import histogram_s1_times
+from . detsim_waveforms import generate_s1_waveform
 
-def test_histogram_s1_times():
-    s1_times = np.array([[1, 2, 3, 4, 5]])
-    buffer_length = 5
-    bin_width     = 1
-    start_time    = 1
+from ..core.system_of_units import ns
 
-    expected = np.array([[1, 1, 1, 1, 1]])
-    h = histogram_s1_times(s1_times, buffer_length, bin_width, start_time)
+def test_generate_s1_waveform_pes():
+    "Test that total pes are conserved in s1 waveform"
+    pes_at_pmts = np.array([[1, 2, 3, 4]])
+    time = np.zeros(pes_at_pmts.shape[-1])
+    buffer_length = 1000 * ns
+    bin_width     = 25  * ns
+    start_time    = 0
 
-    np.testing.assert_allclose(h, expected)
+    wfs = generate_s1_waveform(pes_at_pmts, time, buffer_length, bin_width, start_time)
+    np.testing.assert_allclose(np.sum(pes_at_pmts, axis=1), np.sum(wfs, axis=1), rtol=0.001)
+
+
+def test_generate_s1_waveform_multi():
+    "Test waveform creation for multiple s1s"
+    pes_at_pmts = np.array([[1, 2, 3, 4, 1, 2, 3, 4]])
+    time = np.zeros(pes_at_pmts.shape[-1])
+    time[:-4] = 1000 * ns
+
+    buffer_length = 2000 * ns
+    bin_width     = 25  * ns
+    start_time    = 0
+
+    wfs = generate_s1_waveform(pes_at_pmts, time, buffer_length, bin_width, start_time)
+    np.testing.assert_allclose( wfs[:, :int(buffer_length/bin_width/2)]
+                              , wfs[:,  int(buffer_length/bin_width/2):], rtol=0.001)

invisible_cities/detsim/simulate_s1.py

@@ -66,51 +66,16 @@ def compute_s1_pes_at_pmts(xs      : np.ndarray,
             shape is (number_of_sensors, number_of_hits)
     """
     pes = photons[:, np.newaxis] * lt(xs, ys, zs)
-    pes = np.random.poisson(pes)
     return pes.T
 
 
-def generate_s1_time(size : int=1)->np.ndarray:
+def s1_times_pdf(x):
     """
-    Generate random numbers following two exponential distributions
-    with different probabilities, such that (values in nano-seconds):
-    -with prob 0.1: exp(-x/4.5)
-    -with prob 0.9: exp(-x/100)
+    Implements the s1-times pdf (values in nano-seconds)
+    Reference:
     """
     tau1 = 4.5 * units.ns
     tau2 = 100 * units.ns
     p1 = 0.1
-
-    sel = np.random.uniform(0, 1, size=size)
-    size_1 = np.sum(sel<=p1)
-    size_2 = size - size_1
-    r1 = np.random.exponential(tau1, size_1)
-    r2 = np.random.exponential(tau2, size_2)
-    r  = np.concatenate([r1, r2])
-    np.random.shuffle(r)
-    return r
-
-
-def generate_s1_times_from_pes(s1_pes_at_pmts : np.ndarray,
-                               hit_times      : np.ndarray)->list:
-    """
-    Given the s1_pes_at_pmts, this function returns the times at which the pes
-    are distributed (see generate_S1_time function).
-    It returns a list whose elements are the times at which the photoelectrons in that PMT
-    are generated.
-    Parameters:
-        :s1_pes_at_pmts: np.ndarray
-            the pes at each PMT generated by each hit
-        :hit_times: np.ndarray
-            emission time of each hit
-    Returns:
-        :s1_times: list[np.ndarray,..]
-            Each element are the S1 times for a PMT. If certain sensor
-            do not see any pes, the array is empty.
-    """
-    s1_times = []
-    for pes_at_pmt in s1_pes_at_pmts:
-        repeated_times = np.repeat(hit_times, pes_at_pmt)
-        times = generate_s1_time(np.sum(pes_at_pmt))
-        s1_times.append(times + repeated_times)
-    return s1_times
+    norm = (tau1*p1 + tau2*(1-p1))
+    return (1./norm) * (p1*np.exp(-x/tau1) + (1-p1)*np.exp(-x/tau2))

invisible_cities/detsim/simulate_s1_test.py

@@ -5,8 +5,6 @@
 
 from . simulate_s1 import compute_scintillation_photons
 from . simulate_s1 import compute_s1_pes_at_pmts
-from . simulate_s1 import generate_s1_time
-from . simulate_s1 import generate_s1_times_from_pes
 
 from . light_tables import create_lighttable_function
 
@@ -42,37 +40,4 @@ def test_compute_s1_pes_at_pmts(lighttable_filenames):
     s1_pes_at_pmts = compute_s1_pes_at_pmts(xs, ys, zs, photons, s1_lt)
 
     assert s1_pes_at_pmts.shape == (12, len(xs))
-
-
-def test_generate_s1_time():
-    np.random.seed(1234)
-    times = generate_s1_time(size=2)
-    expected = np.array([ 57.57690962, 153.87864745])
-
-    np.testing.assert_allclose(times, expected)
-
-
-def test_generate_s1_times_from_pes():
-
-    np.random.seed(1234)
-    s1_pes_at_pmts = np.array([[1, 0, 1],
-                               [0, 2, 1]])  # 2 pmts, 3 hits
-    hit_times = np.array([0, 1, 2])
-
-    expected  = [np.array([ 57.57690962, 155.87864745]),
-                 np.array([ 18.29440868,  13.34541421, 170.82805666])]
-    values = generate_s1_times_from_pes(s1_pes_at_pmts, hit_times)
-
-    assert len(values) == len(expected)
-    for exp, val in zip(expected, values):
-        np.testing.assert_allclose(exp, val)
-
-
-@settings(max_examples=500)
-@given(hit_times = arrays(float, 10, elements=floats(min_value=0)))
-def test_generate_s1_times_from_pes_higher_times(hit_times):
-
-    s1_pes_at_pmts = np.ones((1, 10), dtype=int) # just 1 pmt
-    s1_times = generate_s1_times_from_pes(s1_pes_at_pmts, hit_times)[0]
-
-    assert np.all(s1_times>=hit_times)
+