Merge pull request #23 from neutrons/fixed_deat_time

Fixed dead time

reduction/lr_reduction/DeadTimeCorrection.py

@@ -0,0 +1,98 @@
+"""
+    Dead time correction algorithm for single-readout detectors.
+"""
+import time
+import math
+import os
+from mantid.api import *
+from mantid.simpleapi import *
+from mantid.kernel import *
+import numpy as np
+import scipy
+
+def call(InputWorkspace, DeadTime=4.2, TOFStep=100, Paralyzable=False, TOFRange=[0, 0], OutputWorkspace='correction'):
+    """
+        Function to make the algorithm call similar to a normal Mantid call
+    """
+    algo = SingleReadoutDeadTimeCorrection()
+    algo.PyInit()
+    algo.setProperty("InputWorkspace", InputWorkspace)
+    algo.setProperty("DeadTime", DeadTime)
+    algo.setProperty("TOFStep", TOFStep)
+    algo.setProperty("Paralyzable", Paralyzable)
+    algo.setProperty("TOFRange", TOFRange)
+    algo.setProperty("OutputWorkspace", OutputWorkspace)
+    algo.PyExec()
+    return algo.getProperty('OutputWorkspace').value
+
+
+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("DeadTime", 4.2, doc="Dead time in microseconds")
+        self.declareProperty("TOFStep", 100,
+                             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.],
+                                                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
+        dead_time = self.getProperty("DeadTime").value
+        tof_step = self.getProperty("TOFStep").value
+        paralyzing = self.getProperty("Paralyzable").value
+
+        # 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)
+        t_series = np.asarray(_ws_sc.getRun()['proton_charge'].value)
+        non_zero = t_series > 0
+        n_pulses = np.count_nonzero(non_zero)
+        rate = counts_ws.readY(0) / n_pulses
+        tof_bins = counts_ws.readX(0)
+
+         # 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)
+        self.setProperty('OutputWorkspace', counts_ws)
+
+
+AlgorithmFactory.subscribe(SingleReadoutDeadTimeCorrection)

reduction/lr_reduction/event_reduction.py

@@ -7,6 +7,7 @@
 import numpy as np
 
 from . import background
+from . import DeadTimeCorrection
 
 
 def get_wl_range(ws):
@@ -56,13 +57,16 @@ class EventReflectivity(object):
     INSTRUMENT_4B = 1
     DEFAULT_4B_SAMPLE_DET_DISTANCE = 1.83
     DEFAULT_4B_SOURCE_DET_DISTANCE = 15.75
+    DEAD_TIME = 4.2 # Nominally 4.0 microseconds
+    DEAD_TIME_TOF_STEP = 100
 
     def __init__(self, scattering_workspace, direct_workspace,
                  signal_peak, signal_bck, norm_peak, norm_bck,
                  specular_pixel, signal_low_res, norm_low_res,
                  q_min=None, q_step=-0.02, q_max=None,
                  tof_range=None, theta=1.0, instrument=None,
-                 functional_background=False):
+                 functional_background=False, dead_time=False,
+                 paralyzable=False):
         """
             Pixel ranges include the min and max pixels.
 
@@ -80,6 +84,8 @@ def __init__(self, scattering_workspace, direct_workspace,
             :param q_min: value of largest q point
             :param tof_range: TOF range,or None
             :param theta: theta scattering angle in radians
+            :param dead_time: if not zero, dead time correction will be used
+            :param paralyzable: if True, the dead time calculation will use the paralyzable approach
         """
         if instrument in [self.INSTRUMENT_4A, self.INSTRUMENT_4B]:
             self.instrument = instrument
@@ -100,6 +106,8 @@ def __init__(self, scattering_workspace, direct_workspace,
         self._offspec_x_bins = None
         self._offspec_z_bins = None
         self.summing_threshold = None
+        self.dead_time = dead_time
+        self.paralyzable = paralyzable
 
         # Turn on functional background estimation
         self.use_functional_bck = functional_background
@@ -226,6 +234,47 @@ def to_dict(self):
                     dq0=dq0, dq_over_q=dq_over_q, sequence_number=sequence_number,
                     sequence_id=sequence_id)
 
+    def get_dead_time_correction(self):
+        """
+            Compute dead time correction to be applied to the reflectivity curve.
+        """
+        # Scattering workspace
+        tof_min = self._ws_sc.getTofMin()
+        tof_max = self._ws_sc.getTofMax()
+
+        corr_ws = DeadTimeCorrection.call(InputWorkspace=self._ws_sc,
+                                          DeadTime=self.DEAD_TIME,
+                                          TOFStep=self.DEAD_TIME_TOF_STEP,
+                                          Paralyzable=self.paralyzable,
+                                          TOFRange=[tof_min, tof_max],
+                                          OutputWorkspace="corr")
+        corr_sc = corr_ws.readY(0)
+        wl_bins = corr_ws.readX(0) / self.constant
+
+        # Direct beam workspace
+        corr_ws = DeadTimeCorrection.call(InputWorkspace=self._ws_db,
+                                          DeadTime=self.DEAD_TIME,
+                                          TOFStep=self.DEAD_TIME_TOF_STEP,
+                                          Paralyzable=self.paralyzable,
+                                          TOFRange=[tof_min, tof_max],
+                                          OutputWorkspace="corr")
+        corr_db = corr_ws.readY(0)
+
+        # Flip the correction since we are going from TOF to Q
+        dead_time_per_tof = np.flip(corr_sc / corr_db)
+
+        # Compute Q for each TOF bin
+        d_theta = self.gravity_correction(self._ws_sc, wl_bins)
+        q_values_edges = np.flip(4.0 * np.pi / wl_bins * np.sin(self.theta - d_theta))
+        q_values = (q_values_edges[1:] + q_values_edges[:-1]) / 2
+
+        # Interpolate to estimate the dead time correction at the Q values we measured
+        q_middle = (self.q_bins[1:] + self.q_bins[:-1]) / 2
+
+        dead_time_corr = np.interp(q_middle, q_values, dead_time_per_tof)
+
+        return dead_time_corr
+
     def specular(self, q_summing=False, tof_weighted=False, bck_in_q=False,
                  clean=False, normalize=True):
         """
@@ -239,6 +288,11 @@ def specular(self, q_summing=False, tof_weighted=False, bck_in_q=False,
             :param clean: if True, and Q summing is True, then leading artifact will be removed
             :param normalize: if True, and tof_weighted is False, normalization will be skipped
         """
+        # First, let's compute the dead-time correction if we need it
+        # We do this first because the specular calls below may modify the input workspace
+        if self.dead_time:
+            dead_time_corr = self.get_dead_time_correction()
+
         if tof_weighted:
             self.specular_weighted(q_summing=q_summing, bck_in_q=bck_in_q)
         else:
@@ -252,9 +306,15 @@ def specular(self, q_summing=False, tof_weighted=False, bck_in_q=False,
         self.q_bins = self.q_bins[trim:]
 
         # Dead time correction
-        # dead_time = 4e-6
-        #self.refl = self.refl * t_corr_sc / t_corr_db
-        #self.d_refl = self.d_refl * t_corr_sc / t_corr_db
+        if self.dead_time:
+            i_max = np.argmax(dead_time_corr[trim:])
+            i_min = np.argmin(dead_time_corr[trim:])
+            print("Dead time correction: [%g -> %g] at [%g -> %g]" % (dead_time_corr[trim:][i_min],
+                                                                      dead_time_corr[trim:][i_max],
+                                                                      self.q_bins[i_min],
+                                                                      self.q_bins[i_max]))
+            self.refl *= dead_time_corr[trim:]
+            self.d_refl *= dead_time_corr[trim:]
 
         # Remove leading artifact from the wavelength coverage
         # Remember that q_bins is longer than refl by 1 because

reduction/lr_reduction/reduction_template_reader.py

@@ -65,6 +65,10 @@ def __init__(self):
         self.incident_medium_list = ['air']
         self.incident_medium_index_selected = 0
 
+        # Dead time correction
+        self.dead_time:bool = False
+        self.paralyzable:bool = False
+
     def from_dict(self, data_dict, permissible=True):
         r"""
         Update object's attributes with a dictionary with entries of the type  attribute_name: attribute_value.
@@ -147,6 +151,10 @@ def to_xml(self):
         _xml += "<incident_medium_list>%s</incident_medium_list>\n" % str(self.incident_medium_list[0])
         _xml += "<incident_medium_index_selected>%s</incident_medium_index_selected>\n" % str(self.incident_medium_index_selected)
 
+        # Dead time correction
+        _xml += "<dead_time_correction>%s</dead_time_correction>\n" % str(self.dead_time)
+        _xml += "<dead_time_paralyzable>%s</dead_time_paralyzable>\n" % str(self.paralyzable)
+
         _xml += "</RefLData>\n"
 
         return _xml
@@ -248,6 +256,12 @@ def from_xml_element(self, instrument_dom):
             self.incident_medium_list = ['H2O']
             self.incident_medium_index_selected = 0
 
+        # Dead time correction
+        self.dead_time = getBoolElement(instrument_dom, "dead_time_correction",
+                                        default=self.dead_time)
+        self.paralyzable = getBoolElement(instrument_dom, "dead_time_paralyzable",
+                                          default=self.paralyzable)
+
 
 ###### Utility functions to read XML content ########################
 def getText(nodelist):

reduction/lr_reduction/template.py

@@ -223,6 +223,8 @@ def process_from_template_ws(ws_sc, template_data, q_summing=False,
                                                    q_min=q_min, q_step=q_step, q_max=None,
                                                    tof_range=[tof_min, tof_max],
                                                    theta=np.abs(theta),
+                                                   dead_time=template_data.dead_time,
+                                                   paralyzable=template_data.paralyzable,
                                                    functional_background=functional_background,
                                                    instrument=event_reduction.EventReflectivity.INSTRUMENT_4B)