Dead-time correction algorithm and UI elements for settings

docs/developer/documentation.rst

@@ -0,0 +1,13 @@
+=============
+Documentation
+=============
+
+To build the documentation as HTML locally::
+
+   cd docs
+   make html
+
+To test doctest-style code snippets in the documentation using `doctest <https://www.sphinx-doc.org/en/master/usage/extensions/doctest.html>`_::
+
+   cd docs
+   make doctest

docs/developer/environment.rst

@@ -31,6 +31,7 @@ To setup a local development environment, the developers should follow the steps
 The ``environment.yml`` contains all of the dependencies for both the developer and the build servers.
 Update file ``environment.yml`` if dependencies are added to the package.
 
+.. _test-data:
 
 Test Data
 ---------

docs/developer/index.rst

@@ -11,3 +11,4 @@ Development Guide
    integration_test
    ui_test
    release
+   documentation

docs/releasenotes/index.rst

@@ -12,7 +12,7 @@ Notes for major and minor releases. Notes for Patch releases are deferred.
 
 **Of interest to the User**:
 
-- PR #XYZ: one-liner description
+- PR #95: Optional dead-time correction (disabled by default)
 
 **Of interest to the Developer:**
 

docs/user/dead_time_correction.rst

@@ -0,0 +1,111 @@
+.. _dead_time_correction:
+
+SingleReadoutDeadTimeCorrection
+===============================
+
+Dead time is the time after an event that a detector is not able to detect another event.
+For a paralyzable detector, an event that happens during the dead time restarts the dead time. For
+a non-paralyzable detector, the event is simply lost and does not cause additional dead time.
+
+Dead-time correction corrects for detector dead time by weighing the events according to:
+
+.. math:: N = M \frac{1}{(1-\mathrm{rate} \cdot (\frac{t_{\mathrm{dead}}}{t_{\mathrm{bin}}}))}
+
+for non-paralyzable detectors and
+
+.. math:: N = M \frac{\mathrm{Re} (\mathrm{W}(-\mathrm{rate} \cdot (\frac{t_{\mathrm{dead}}}{t_{\mathrm{bin}}})) )}{\frac{\mathrm{rate}}{t_{\mathrm{bin}}}}
+
+for paralyzable detectors, where
+
+| :math:`N` = true count
+| :math:`M` = measured count
+| :math:`t_{\mathrm{dead}}` = dead time
+| :math:`t_{\mathrm{bin}}` = TOF bin width
+| :math:`\mathrm{rate}` = measured count rate
+| :math:`\mathrm{W}` = `Lambert W function <https://en.wikipedia.org/wiki/Lambert_W_function>`_
+
+The class ``SingleReadoutDeadTimeCorrection`` is a Mantid-style algorithm for computing the
+dead-time correction for an event workspace. One can optionally include error events in the
+dead-time computation.
+
+Properties
+----------
+
+.. list-table::
+   :widths: 20 20 20 20 20
+   :header-rows: 1
+
+   * - Name
+     - Direction
+     - Type
+     - Default
+     - Description
+   * - InputWorkspace
+     - Input
+     - EventWorkspace
+     - Mandatory
+     - Input workspace used to compute dead-time correction
+   * - InputErrorEventsWorkspace
+     - Input
+     - EventWorkspace
+     -
+     - Input workspace with error events used to compute dead-time correction
+   * - DeadTime
+     - Input
+     - number
+     - 4.2
+     - Dead time in microseconds
+   * - TOFStep
+     - Input
+     - number
+     - 100.0
+     - TOF bins to compute dead-time correction, in microseconds
+   * - Paralyzable
+     - Input
+     - boolean
+     - False
+     - If True, paralyzable correction will be applied, non-paralyzable otherwise
+   * - TOFRange
+     - Input
+     - dbl list
+     - [0.0, 0.0]
+     - TOF range to use to compute dead-time correction
+   * - OutputWorkspace
+     - Output
+     - MatrixWorkspace
+     - Mandatory
+     - Output workspace containing the dead-time correction factor for each TOF bin
+
+Usage
+-----
+Example using ``SingleReadoutDeadTimeCorrection`` (requires :any:`test-data` files)
+
+.. testcode::
+
+    import mantid.simpleapi as api
+    import os
+    import sys
+    from pathlib import Path
+    from reflectivity_ui.interfaces.data_handling import DeadTimeCorrection
+    from reflectivity_ui.interfaces.data_handling.instrument import mantid_algorithm_exec
+    # Load events
+    path = Path().resolve().parent / "test" / "data" / "reflectivity_ui-data" / "REF_M_42112.nxs.h5"
+    path = path.as_posix()
+    ws = api.LoadEventNexus(Filename=path, OutputWorkspace="raw_events")
+    # Load error events
+    err_ws = api.LoadErrorEventsNexus(path)
+    # Compute dead-time correction
+    tof_min = ws.getTofMin()
+    tof_max = ws.getTofMax()
+    corr_ws = mantid_algorithm_exec(
+        DeadTimeCorrection.SingleReadoutDeadTimeCorrection,
+        InputWorkspace=ws,
+        InputErrorEventsWorkspace=err_ws,
+        Paralyzable=False,
+        DeadTime=4.2,
+        TOFStep=100.0,
+        TOFRange=[tof_min, tof_max],
+        OutputWorkspace="corr",
+    )
+    # Apply dead-time correction
+    ws = api.Multiply(ws, corr_ws, OutputWorkspace=str(ws))

docs/user/index.rst

@@ -4,3 +4,8 @@ User Guide
 ==========
 
 TODO
+
+.. toctree::
+   :maxdepth: 2
+
+   dead_time_correction

reflectivity_ui/interfaces/configuration.py

@@ -37,6 +37,11 @@ class Configuration(object):
     polynomial_stitching = False
     polynomial_stitching_degree = 3
     polynomial_stitching_points = 3
+    # Dead time options
+    apply_deadtime = False
+    paralyzable_deadtime = True
+    deadtime_value = 4.2
+    deadtime_tof_step = 100
 
     def __init__(self, settings=None):
         self.instrument = Instrument()
@@ -227,6 +232,12 @@ def to_q_settings(self, settings):
         settings.setValue("do_final_rebin", self.do_final_rebin)
         settings.setValue("final_rebin_step", self.final_rebin_step)
 
+        # Dead time options
+        settings.setValue("apply_deadtime", self.apply_deadtime)
+        settings.setValue("paralyzable_deadtime", self.paralyzable_deadtime)
+        settings.setValue("deadtime_value", self.deadtime_value)
+        settings.setValue("deadtime_tof_step", self.deadtime_tof_step)
+
         # Off-specular options
         settings.setValue("off_spec_x_axis", self.off_spec_x_axis)
         settings.setValue("off_spec_slice", self.off_spec_slice)
@@ -325,6 +336,12 @@ def _verify_true(parameter, default):
         Configuration.do_final_rebin = _verify_true("do_final_rebin", self.do_final_rebin)
         Configuration.final_rebin_step = float(settings.value("final_rebin_step", self.final_rebin_step))
 
+        # Dead time options
+        Configuration.apply_deadtime = _verify_true("apply_deadtime", self.apply_deadtime)
+        Configuration.paralyzable_deadtime = _verify_true("paralyzable_deadtime", self.paralyzable_deadtime)
+        Configuration.deadtime_value = float(settings.value("deadtime_value", self.deadtime_value))
+        Configuration.deadtime_tof_step = float(settings.value("deadtime_tof_step", self.deadtime_tof_step))
+
         # Off-specular options
         self.off_spec_x_axis = int(settings.value("off_spec_x_axis", self.off_spec_x_axis))
         self.off_spec_slice = _verify_true("off_spec_slice", self.off_spec_slice)
@@ -377,3 +394,7 @@ def setup_default_values(cls):
         cls.polynomial_stitching = False
         cls.polynomial_stitching_degree = 3
         cls.polynomial_stitching_points = 3
+        cls.apply_deadtime = False
+        cls.paralyzable_deadtime = True
+        cls.deadtime_value = 4.2
+        cls.deadtime_tof_step = 100

reflectivity_ui/interfaces/data_handling/DeadTimeCorrection.py

@@ -0,0 +1,131 @@
+"""
+Dead time correction algorithm for single-readout detectors.
+"""
+
+import numpy as np
+import scipy
+from mantid.api import (
+    AlgorithmFactory,
+    IEventWorkspaceProperty,
+    MatrixWorkspaceProperty,
+    PropertyMode,
+    PythonAlgorithm,
+)
+from mantid.kernel import Direction, FloatArrayLengthValidator, FloatArrayProperty
+from mantid.simpleapi import Rebin, SumSpectra, logger
+
+
+class SingleReadoutDeadTimeCorrection(PythonAlgorithm):
+    def category(self):
+        return "Reflectometry\\SNS"
+
+    def name(self):
+        return "SingleReadoutDeadTimeCorrection"
+
+    def version(self):
+        return 1
+
+    def summary(self):
+        return "Single read-out dead time correction calculation"
+
+    def PyInit(self):
+        self.declareProperty(
+            IEventWorkspaceProperty("InputWorkspace", "", Direction.Input),
+            "Input workspace use to compute dead time correction",
+        )
+        self.declareProperty(
+            IEventWorkspaceProperty("InputErrorEventsWorkspace", "", Direction.Input, PropertyMode.Optional),
+            "Input workspace with error events use to compute dead time correction",
+        )
+        self.declareProperty("DeadTime", 4.2, doc="Dead time in microseconds")
+        self.declareProperty(
+            "TOFStep",
+            100.0,
+            doc="TOF bins to compute deadtime correction for, in microseconds",
+        )
+        self.declareProperty(
+            "Paralyzable",
+            False,
+            doc="If true, paralyzable correction will be applied, non-paralyzing otherwise",
+        )
+        self.declareProperty(
+            FloatArrayProperty(
+                "TOFRange",
+                [0.0, 0.0],
+                FloatArrayLengthValidator(2),
+                direction=Direction.Input,
+            ),
+            "TOF range to use",
+        )
+        self.declareProperty(
+            MatrixWorkspaceProperty("OutputWorkspace", "", Direction.Output),
+            "Output workspace",
+        )
+
+    def PyExec(self):
+        # Event data must include error events (all triggers on the detector)
+        ws_event_data = self.getProperty("InputWorkspace").value
+        ws_error_events = self.getProperty("InputErrorEventsWorkspace").value
+        dead_time = self.getProperty("DeadTime").value
+        tof_step = self.getProperty("TOFStep").value
+        paralyzing = self.getProperty("Paralyzable").value
+        output_workspace = self.getPropertyValue("OutputWorkspace")
+
+        # Rebin the data according to the tof_step we want to compute the correction with
+        tof_min, tof_max = self.getProperty("TOFRange").value
+        if tof_min == 0 and tof_max == 0:
+            tof_min = ws_event_data.getTofMin()
+            tof_max = ws_event_data.getTofMax()
+        logger.notice("TOF range: %f %f" % (tof_min, tof_max))
+        _ws_sc = Rebin(
+            InputWorkspace=ws_event_data,
+            Params="%s,%s,%s" % (tof_min, tof_step, tof_max),
+            PreserveEvents=False,
+        )
+
+        # Get the total number of counts on the detector for each TOF bin per pulse
+        counts_ws = SumSpectra(_ws_sc, OutputWorkspace=output_workspace)
+
+        # If we have error events, add them since those are also detector triggers
+        if ws_error_events is not None:
+            _errors = Rebin(
+                InputWorkspace=ws_error_events,
+                Params="%s,%s,%s" % (tof_min, tof_step, tof_max),
+                PreserveEvents=False,
+            )
+            counts_ws += _errors
+
+        t_series = np.asarray(_ws_sc.getRun()["proton_charge"].value)
+        non_zero = t_series > 0
+        n_pulses = np.count_nonzero(non_zero)
+
+        # If we skip pulses, we need to account for them when computing the
+        # instantaneous rate
+        chopper_speed = _ws_sc.getRun()["SpeedRequest1"].value[0]
+        n_pulses = n_pulses * chopper_speed / 60.0
+        rate = counts_ws.readY(0) / n_pulses
+
+        # Compute the dead time correction for each TOF bin
+        if paralyzing:
+            true_rate = -scipy.special.lambertw(-rate * dead_time / tof_step).real / dead_time
+            corr = true_rate / (rate / tof_step)
+            # If we have no events, set the correction to 1 otherwise we will get a nan
+            # from the equation above.
+            corr[rate == 0] = 1
+        else:
+            corr = 1 / (1 - rate * dead_time / tof_step)
+
+        if np.min(corr) < 0:
+            error = "Corrupted dead time correction:\n" + "  Reflected: %s\n" % corr
+            logger.error(error)
+
+        counts_ws.setY(0, corr)
+
+        # We don't compute an error on the dead time correction, so set it to zero
+        counts_ws.setE(0, 0 * corr)
+        counts_ws.setDistribution(True)
+
+        self.setProperty("OutputWorkspace", counts_ws)
+
+
+AlgorithmFactory.subscribe(SingleReadoutDeadTimeCorrection)

reflectivity_ui/interfaces/data_handling/data_set.py

@@ -382,7 +382,7 @@ def load(self, update_parameters=True, progress=None):
             progress(5, "Filtering data...", out_of=100.0)
 
         try:
-            xs_list = self.configuration.instrument.load_data(self.file_path)
+            xs_list = self.configuration.instrument.load_data(self.file_path, self.configuration)
             logging.info("%s loaded: %s xs", self.file_path, len(xs_list))
         except RuntimeError as run_err:
             logging.error("Could not load file(s) {}\n   {}".format(str(self.file_path), run_err))