From 8b3eb50630f0de86d00854f9657203c665182721 Mon Sep 17 00:00:00 2001 From: Mathieu Doucet Date: Wed, 20 Mar 2024 16:55:18 -0400 Subject: [PATCH] add error events to deadtime corr --- reduction/lr_reduction/DeadTimeCorrection.py | 20 +++++++++++++------- reduction/lr_reduction/event_reduction.py | 18 +++++++++++++++--- 2 files changed, 28 insertions(+), 10 deletions(-) diff --git a/reduction/lr_reduction/DeadTimeCorrection.py b/reduction/lr_reduction/DeadTimeCorrection.py index fe10aee..1bbe993 100644 --- a/reduction/lr_reduction/DeadTimeCorrection.py +++ b/reduction/lr_reduction/DeadTimeCorrection.py @@ -1,22 +1,20 @@ """ 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'): +def call(InputWorkspace, InputErrorEventsWorkspace=None, 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("InputErrorEventsWorkspace", InputErrorEventsWorkspace) algo.setProperty("DeadTime", DeadTime) algo.setProperty("TOFStep", TOFStep) algo.setProperty("Paralyzable", Paralyzable) @@ -43,8 +41,10 @@ def summary(self): 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, + 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") @@ -56,6 +56,7 @@ def PyInit(self): 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 @@ -70,13 +71,18 @@ def PyExec(self): # Get the total number of counts on the detector for each TOF bin per pulse counts_ws = SumSpectra(_ws_sc) + + # 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) rate = counts_ws.readY(0) / n_pulses - tof_bins = counts_ws.readX(0) - # Compute the dead time correction for each TOF bin + # 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) diff --git a/reduction/lr_reduction/event_reduction.py b/reduction/lr_reduction/event_reduction.py index b8d6dd7..0771a27 100644 --- a/reduction/lr_reduction/event_reduction.py +++ b/reduction/lr_reduction/event_reduction.py @@ -237,12 +237,17 @@ def to_dict(self): def get_dead_time_correction(self): """ Compute dead time correction to be applied to the reflectivity curve. + The method will also try to load the error events from each of the + data files to ensure that we properly estimate the dead time correction. """ # Scattering workspace tof_min = self._ws_sc.getTofMin() tof_max = self._ws_sc.getTofMax() + run_number = self._ws_sc.getRun().getProperty("run_number").value + error_ws = api.LoadErrorEventsNexus(run_number) corr_ws = DeadTimeCorrection.call(InputWorkspace=self._ws_sc, + InputErrorEventsWorkspace=error_ws, DeadTime=self.DEAD_TIME, TOFStep=self.DEAD_TIME_TOF_STEP, Paralyzable=self.paralyzable, @@ -252,7 +257,10 @@ def get_dead_time_correction(self): wl_bins = corr_ws.readX(0) / self.constant # Direct beam workspace + run_number = self._ws_db.getRun().getProperty("run_number").value + error_ws = api.LoadErrorEventsNexus(run_number) corr_ws = DeadTimeCorrection.call(InputWorkspace=self._ws_db, + InputErrorEventsWorkspace=error_ws, DeadTime=self.DEAD_TIME, TOFStep=self.DEAD_TIME_TOF_STEP, Paralyzable=self.paralyzable, @@ -270,11 +278,15 @@ def get_dead_time_correction(self): # 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) - + + # Cleanup + api.DeleteWorkspace(corr_ws) + api.DeleteWorkspace(error_ws) + return dead_time_corr - + def specular(self, q_summing=False, tof_weighted=False, bck_in_q=False, clean=False, normalize=True): """