Add PMT dark count simulation

fuse/context.py

@@ -64,6 +64,7 @@
 # Plugins to simulate PMTs and DAQ
 pmt_and_daq_plugins = [
     fuse.pmt_and_daq.PMTAfterPulses,
+    fuse.pmt_and_daq.DarkCounts,
     fuse.pmt_and_daq.PhotonSummary,
     fuse.pmt_and_daq.PulseWindow,
     fuse.pmt_and_daq.PMTResponseAndDAQ,

fuse/dtypes.py

@@ -74,5 +74,5 @@
     (("Photon creates a double photo-electron emission", "dpe"), np.bool_),
     (("Sampled PMT gain for the photon", "photon_gain"), np.int32),
     (("ID of the cluster creating the photon", "cluster_id"), np.int32),
-    (("Type of the photon. S1 (1), S2 (2) or PMT AP (0)", "photon_type"), np.int8),
+    (("Type of the photon. S1 (1), S2 (2), PMT AP (0) or dark count (3)", "photon_type"), np.int8),
 ]

fuse/plugins/pmt_and_daq/__init__.py

@@ -9,3 +9,6 @@
 
 from . import photon_pulses
 from .photon_pulses import *
+
+from . import dark_counts
+from .dark_counts import *

fuse/plugins/pmt_and_daq/dark_counts.py

@@ -0,0 +1,212 @@
+import strax
+import straxen
+import logging
+import numpy as np
+
+from ...plugin import FuseBasePlugin
+from ...dtypes import propagated_photons_fields
+from ...common import pmt_gains, build_photon_propagation_output
+from ...common import (
+    init_spe_scaling_factor_distributions,
+    photon_gain_calculation,
+)
+
+export, __all__ = strax.exporter()
+
+logging.basicConfig(handlers=[logging.StreamHandler()])
+log = logging.getLogger("fuse.pmt_and_daq.dark_counts")
+
+
+@export
+class DarkCounts(FuseBasePlugin):
+    """Plugin to simulate dark counts in a window around the physics
+    interaction."""
+
+    __version__ = "0.0.1"
+
+    depends_on = "microphysics_summary"
+    provides = "dark_count_photons"
+    data_kind = "dark_count_photons"
+
+    save_when = strax.SaveWhen.TARGET
+
+    dtype = propagated_photons_fields + strax.time_fields
+
+    # Config options
+
+    enable_dark_counts = straxen.URLConfig(
+        default=False,
+        type=bool,
+        track=True,
+        help="Decide if you want to to enable dark count simulation",
+    )
+
+    # Get this value from a database. For now lets just set it to 50 Hz per PMT
+    dark_count_rate = straxen.URLConfig(
+        default=50 * 494,
+        type=(int, float),
+        track=True,
+        help="Rate of dark counts in Hz combined for all PMTs",
+    )
+
+    dark_count_left_window = straxen.URLConfig(
+        default=2e6,
+        type=int,
+        track=True,
+        help="Left window of the dark count simulation",
+    )
+
+    dark_count_right_window = straxen.URLConfig(
+        default=2e6,
+        type=int,
+        track=True,
+        help="Right window of the dark count simulation",
+    )
+
+    # Add a default pointing to the database
+    dark_count_probability_per_pmt = straxen.URLConfig(
+        default=None,
+        track=True,
+        help="Probability of dark counts per PMT",
+    )
+
+    # reconsider this one for dark counts....
+    p_double_pe_emision = straxen.URLConfig(
+        default="take://resource://"
+        "SIMULATION_CONFIG_FILE.json?&fmt=json"
+        "&take=p_double_pe_emision",
+        type=(int, float),
+        cache=True,
+        help="Probability of double photo-electron emission",
+    )
+
+    pmt_circuit_load_resistor = straxen.URLConfig(
+        default="take://resource://"
+        "SIMULATION_CONFIG_FILE.json?&fmt=json"
+        "&take=pmt_circuit_load_resistor",
+        type=(int, float),
+        cache=True,
+        help="PMT circuit load resistor [kg m^2/(s^3 A)]",
+    )
+
+    digitizer_bits = straxen.URLConfig(
+        default="take://resource://SIMULATION_CONFIG_FILE.json?&fmt=json&take=digitizer_bits",
+        type=(int, float),
+        cache=True,
+        help="Number of bits of the digitizer boards",
+    )
+
+    digitizer_voltage_range = straxen.URLConfig(
+        default="take://resource://"
+        "SIMULATION_CONFIG_FILE.json?&fmt=json"
+        "&take=digitizer_voltage_range",
+        type=(int, float),
+        cache=True,
+        help="Voltage range of the digitizer boards [V]",
+    )
+
+    gain_model_mc = straxen.URLConfig(
+        default="cmt://to_pe_model?version=ONLINE&run_id=plugin.run_id",
+        infer_type=False,
+        help="PMT gain model",
+    )
+
+    photon_area_distribution = straxen.URLConfig(
+        default="simple_load://resource://simulation_config://"
+        "SIMULATION_CONFIG_FILE.json?"
+        "&key=photon_area_distribution"
+        "&fmt=csv",
+        cache=True,
+        help="Photon area distribution",
+    )
+
+    n_tpc_pmts = straxen.URLConfig(
+        type=int,
+        help="Number of PMTs in the TPC",
+    )
+
+    def setup(self):
+        super().setup()
+
+        self.gains = pmt_gains(
+            self.gain_model_mc,
+            digitizer_voltage_range=self.digitizer_voltage_range,
+            digitizer_bits=self.digitizer_bits,
+            pmt_circuit_load_resistor=self.pmt_circuit_load_resistor,
+        )
+
+        self.pmt_mask = np.array(self.gains) > 0  # Converted from to pe (from cmt by default)
+        self.turned_off_pmts = np.nonzero(np.array(self.gains) == 0)[0]
+
+        self.spe_scaling_factor_distributions = init_spe_scaling_factor_distributions(
+            self.photon_area_distribution
+        )
+
+    def compute(self, interactions_in_roi, start, end):
+        if not self.enable_dark_counts or (len(interactions_in_roi) == 0):
+            return np.zeros(0, dtype=self.dtype)
+
+        single_simulation_window = np.zeros(len(interactions_in_roi), dtype=strax.interval_dtype)
+        single_simulation_window["time"] = interactions_in_roi["time"]
+        single_simulation_window["dt"] = np.ones(len(interactions_in_roi))
+
+        simulation_windows = strax.concat_overlapping_hits(
+            single_simulation_window,
+            extensions=(self.dark_count_left_window, self.dark_count_right_window),
+            pmt_channels=(0, self.n_tpc_pmts),
+            start=start,
+            end=end,
+        )
+
+        # Get the number of dark counts in the simulation window
+        expected_dark_counts_in_simulation_window = (
+            simulation_windows["length"].astype(np.int64) * self.dark_count_rate / 1e9
+        )
+        dark_counts_in_simulation_window = self.rng.poisson(
+            expected_dark_counts_in_simulation_window
+        )
+
+        dark_count_times = get_random_times(
+            simulation_windows["length"], dark_counts_in_simulation_window, self.rng
+        )
+        dark_count_times += np.repeat(simulation_windows["time"], dark_counts_in_simulation_window)
+
+        # distribute the dark counts to the PMTs
+        if self.dark_count_probability_per_pmt is None:
+            log.warning("No dark count probability per PMT provided. Using uniform distribution.")
+            dark_count_channels = self.rng.choice(self.n_tpc_pmts, len(dark_count_times))
+        else:
+            dark_count_channels = self.rng.choice(
+                self.n_tpc_pmts, len(dark_count_times), p=self.dark_count_probability_per_pmt
+            )
+
+        # We for sure need to update the inputs for this step
+        photon_gains, photon_is_dpe = photon_gain_calculation(
+            _photon_channels=dark_count_channels,
+            p_double_pe_emision=self.p_double_pe_emision,
+            gains=self.gains,
+            spe_scaling_factor_distributions=self.spe_scaling_factor_distributions,
+            rng=self.rng,
+        )
+        photon_is_dpe = False
+
+        # now build the output
+        result = build_photon_propagation_output(
+            dtype=self.dtype,
+            _photon_timings=dark_count_times,
+            _photon_channels=dark_count_channels,
+            _photon_gains=photon_gains,
+            _photon_is_dpe=photon_is_dpe,
+            _cluster_id=0,  # rethink this part...
+            photon_type=3,
+        )
+
+        return result
+
+
+# For sure this can be done more efficiently
+def get_random_times(interval_length, number_of_entries, rng):
+    times = []
+    for max_time, n in zip(interval_length, number_of_entries):
+        times.append(rng.uniform(0, max_time, n))
+    return np.concatenate(times)

fuse/plugins/pmt_and_daq/photon_summary.py

@@ -10,7 +10,12 @@ class PhotonSummary(VerticalMergerPlugin):
     """Plugin that concatenates propagated photons for S1, S2 and PMT
     afterpulses."""
 
-    depends_on = ("propagated_s2_photons", "propagated_s1_photons", "pmt_afterpulses")
+    depends_on = (
+        "propagated_s2_photons",
+        "propagated_s1_photons",
+        "pmt_afterpulses",
+        "dark_count_photons",
+    )
 
     provides = "photon_summary"
     data_kind = "propagated_photons"

fuse/plugins/truth_information/peak_truth.py

@@ -26,6 +26,7 @@ class PeakTruth(strax.OverlapWindowPlugin):
         ("s1_photons_in_peak", np.int32),
         ("s2_photons_in_peak", np.int32),
         ("ap_photons_in_peak", np.int32),
+        ("dark_count_photons_in_peak", np.int32),
         ("pi_photons_in_peak", np.int32),
         ("raw_area_truth", np.float32),
         ("observable_energy_truth", np.float32),
@@ -114,6 +115,7 @@ def compute(self, interactions_in_roi, propagated_photons, peaks):
             "s1": 1,
             "s2": 2,
             "ap": 0,
+            "dark_count": 3,
         }
 
         for i in range(n_peaks):

fuse/plugins/truth_information/records_truth.py

@@ -22,6 +22,7 @@ class RecordsTruth(strax.Plugin):
         (("Number of S1 photons in record", "s1_photons_in_record"), np.int32),
         (("Number of S2 photons in record", "s2_photons_in_record"), np.int32),
         (("Number of AP photons in record", "ap_photons_in_record"), np.int32),
+        (("Number of dark counts in record", "dark_count_photons_in_record"), np.int32),
         (("Sum of the photon gains", "raw_area"), np.float32),
     ]
 
@@ -98,6 +99,9 @@ def compute(self, propagated_photons, raw_records):
             result["s1_photons_in_record"][result_mask] = result_buffer["s1_photons_in_record"]
             result["s2_photons_in_record"][result_mask] = result_buffer["s2_photons_in_record"]
             result["ap_photons_in_record"][result_mask] = result_buffer["ap_photons_in_record"]
+            result["dark_count_photons_in_record"][result_mask] = result_buffer[
+                "dark_count_photons_in_record"
+            ]
 
         return result
 
@@ -110,6 +114,7 @@ def fill_result_buffer(list_input, result_buffer):
         result_buffer["s1_photons_in_record"][i] = np.sum(photons["photon_type"] == 1)
         result_buffer["s2_photons_in_record"][i] = np.sum(photons["photon_type"] == 2)
         result_buffer["ap_photons_in_record"][i] = np.sum(photons["photon_type"] == 0)
+        result_buffer["dark_count_photons_in_record"][i] = np.sum(photons["photon_type"] == 3)
 
 
 def split_photons_by_channel(propagated_photons):