diff --git a/reduction/lr_reduction/event_reduction.py b/reduction/lr_reduction/event_reduction.py
index 354e4c7..ed69866 100644
--- a/reduction/lr_reduction/event_reduction.py
+++ b/reduction/lr_reduction/event_reduction.py
@@ -1,6 +1,7 @@
 """
-    Event based reduction for the Liquids Reflectometer
+Event based reduction for the Liquids Reflectometer
 """
+
 import datetime
 import json
 import os
@@ -17,29 +18,30 @@
 NEUTRON_MASS = 1.675e-27  # kg
 
 # Attenuators: PV name and thickness in cm
-CD_ATTENUATORS = [['BL4B:Actuator:50MRb', 0.0058],
-                  ['BL4B:Actuator:100MRb', 0.0122],
-                  ['BL4B:Actuator:200MRb', 0.0244], # Uncalibrated
-                  ['BL4B:Actuator:400MRb', 0.0488], # Uncalibrated
-                  ]
+CD_ATTENUATORS = [
+    ["BL4B:Actuator:50MRb", 0.0058],
+    ["BL4B:Actuator:100MRb", 0.0122],
+    ["BL4B:Actuator:200MRb", 0.0244],  # Uncalibrated
+    ["BL4B:Actuator:400MRb", 0.0488],  # Uncalibrated
+]
 
 
 def get_wl_range(ws):
     """
-        Determine TOF range from the data
+    Determine TOF range from the data
 
-        :param ws: (Mantid workspace) workspace to work with
+    :param ws: (Mantid workspace) workspace to work with
 
-        :return: (list) [min, max] wavelength range
+    :return: (list) [min, max] wavelength range
     """
     run_object = ws.getRun()
 
-    wl = run_object.getProperty('LambdaRequest').value[0]
-    chopper_speed = run_object.getProperty('SpeedRequest1').value[0]
+    wl = run_object.getProperty("LambdaRequest").value[0]
+    chopper_speed = run_object.getProperty("SpeedRequest1").value[0]
 
     # Cut the edges by using a width of 2.6 A
-    wl_min = (wl - 1.3 * 60.0 / chopper_speed)
-    wl_max = (wl + 1.3 * 60.0 / chopper_speed)
+    wl_min = wl - 1.3 * 60.0 / chopper_speed
+    wl_max = wl + 1.3 * 60.0 / chopper_speed
 
     return [wl_min, wl_max]
 
@@ -61,11 +63,11 @@ def get_q_binning(q_min=0.001, q_max=0.15, q_step=-0.02):
         An array of Q values based on the specified binning.
     """
     if q_step > 0:
-        n_steps = int((q_max-q_min)/q_step)
+        n_steps = int((q_max - q_min) / q_step)
         return q_min + np.asarray([q_step * i for i in range(n_steps)])
     else:
-        _step = 1.0+np.abs(q_step)
-        n_steps = int(np.log(q_max/q_min)/np.log(_step))
+        _step = 1.0 + np.abs(q_step)
+        n_steps = int(np.log(q_max / q_min) / np.log(_step))
         return q_min * np.asarray([_step**i for i in range(n_steps)])
 
 
@@ -94,40 +96,40 @@ def get_attenuation_info(ws):
 
 def read_settings(ws):
     """
-        Read settings file and return values for the given timestamp
+    Read settings file and return values for the given timestamp
 
-        :param ws: Mantid workspace
-        :return: (dict) dictionary with settings
+    :param ws: Mantid workspace
+    :return: (dict) dictionary with settings
     """
     settings = dict()
     package_dir, _ = os.path.split(__file__)
 
-    t = ws.getRun()['start_time'].value.split('T')[0]
+    t = ws.getRun()["start_time"].value.split("T")[0]
     timestamp = datetime.date.fromisoformat(t)
 
-    with open(os.path.join(package_dir, 'settings.json'), 'r') as fd:
+    with open(os.path.join(package_dir, "settings.json"), "r") as fd:
         data = json.load(fd)
         for key in data.keys():
             chosen_value = None
             delta_time = None
             for item in data[key]:
-                valid_from = datetime.date.fromisoformat(item['from'])
+                valid_from = datetime.date.fromisoformat(item["from"])
                 delta = valid_from - timestamp
                 if delta_time is None or (delta.total_seconds() < 0 and delta > delta_time):
                     delta_time = delta
-                    chosen_value = item['value']
+                    chosen_value = item["value"]
             settings[key] = chosen_value
     return settings
 
 
 def process_attenuation(ws, thickness=0):
     """
-        Correct for absorption by assigning weight to each neutron event
+    Correct for absorption by assigning weight to each neutron event
 
-        :param ws: (Mantid workspace) workspace to correct
-        :param thickness: (float) attenuator thickness in cm (default is 0).
+    :param ws: (Mantid workspace) workspace to correct
+    :param thickness: (float) attenuator thickness in cm (default is 0).
 
-        :return: (Mantid workspace) corrected workspace
+    :return: (Mantid workspace) corrected workspace
     """
     settings = read_settings(ws)
     if "source-det-distance" in settings:
@@ -138,19 +140,22 @@ def process_attenuation(ws, thickness=0):
     constant = 1e-4 * NEUTRON_MASS * SDD / PLANCK_CONSTANT
 
     package_dir, _ = os.path.split(__file__)
-    mu_abs = np.loadtxt(os.path.join(package_dir, 'Cd-abs-factors.txt')).T
+    mu_abs = np.loadtxt(os.path.join(package_dir, "Cd-abs-factors.txt")).T
     wl_model = mu_abs[0]
 
     # Turn model into a histogram
     wl_step = wl_model[-1] - wl_model[-2]
     final_wl = wl_model[-1] + wl_step
-    wl_model = np.append(wl_model, [final_wl,])
+    wl_model = np.append(
+        wl_model,
+        [
+            final_wl,
+        ],
+    )
     mu_model = mu_abs[1]
     tof_model = constant * wl_model
-    transmission = 1/np.exp(-mu_model * thickness)
-    transmission_ws = api.CreateWorkspace(OutputWorkspace='transmission',
-                                          DataX=tof_model, DataY=transmission,
-                                          UnitX='TOF', NSpec=1)
+    transmission = 1 / np.exp(-mu_model * thickness)
+    transmission_ws = api.CreateWorkspace(OutputWorkspace="transmission", DataX=tof_model, DataY=transmission, UnitX="TOF", NSpec=1)
 
     ws = api.Multiply(ws, transmission_ws, OutputWorkspace=str(ws))
     return ws
@@ -158,47 +163,50 @@ def process_attenuation(ws, thickness=0):
 
 def get_dead_time_correction(ws, template_data):
     """
-        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.
+    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.
 
-        :param ws: (Mantid worksapce) workspace with raw data to compute correction for
-        :param template_data: (reduction_template_reader.ReductionParameters)
-            reduction parameters
+    :param ws: (Mantid worksapce) workspace with raw data to compute correction for
+    :param template_data: (reduction_template_reader.ReductionParameters)
+        reduction parameters
 
-        :return: (Mantid workspace) workspace with dead time correction to apply
+    :return: (Mantid workspace) workspace with dead time correction to apply
     """
     tof_min = ws.getTofMin()
     tof_max = ws.getTofMax()
 
     run_number = ws.getRun().getProperty("run_number").value
     error_ws = api.LoadErrorEventsNexus("REF_L_%s" % run_number)
-    corr_ws = mantid_algorithm_exec(DeadTimeCorrection.SingleReadoutDeadTimeCorrection,
-                                    InputWorkspace=ws,
-                                    InputErrorEventsWorkspace=error_ws,
-                                    Paralyzable=template_data.paralyzable,
-                                    DeadTime=template_data.dead_time_value,
-                                    TOFStep=template_data.dead_time_tof_step,
-                                    TOFRange=[tof_min, tof_max],
-                                    OutputWorkspace="corr")
+    corr_ws = mantid_algorithm_exec(
+        DeadTimeCorrection.SingleReadoutDeadTimeCorrection,
+        InputWorkspace=ws,
+        InputErrorEventsWorkspace=error_ws,
+        Paralyzable=template_data.paralyzable,
+        DeadTime=template_data.dead_time_value,
+        TOFStep=template_data.dead_time_tof_step,
+        TOFRange=[tof_min, tof_max],
+        OutputWorkspace="corr",
+    )
     corr_ws = api.Rebin(corr_ws, [tof_min, 10, tof_max])
     return corr_ws
 
+
 def apply_dead_time_correction(ws, template_data):
     """
-        Apply dead time correction, and ensure that it is done only once
-        per workspace.
+    Apply dead time correction, and ensure that it is done only once
+    per workspace.
 
-        :param ws: (Mantid workspace) workspace with raw data to compute correction for
-        :param template_data: (reduction_template_reader.ReductionParameters)
-            reduction parameters
+    :param ws: (Mantid workspace) workspace with raw data to compute correction for
+    :param template_data: (reduction_template_reader.ReductionParameters)
+        reduction parameters
 
-        :return: (Mantid workspace) workspace with dead time correction applied
+    :return: (Mantid workspace) workspace with dead time correction applied
     """
-    if 'dead_time_applied' not in ws.getRun():
+    if "dead_time_applied" not in ws.getRun():
         corr_ws = get_dead_time_correction(ws, template_data)
         ws = api.Multiply(ws, corr_ws, OutputWorkspace=str(ws))
-        api.AddSampleLog(Workspace=ws, LogName="dead_time_applied", LogText='1', LogType="Number")
+        api.AddSampleLog(Workspace=ws, LogName="dead_time_applied", LogText="1", LogType="Number")
     return ws
 
 
@@ -244,15 +252,30 @@ class EventReflectivity(object):
     DEFAULT_4B_SAMPLE_DET_DISTANCE = 1.83
     DEFAULT_4B_SOURCE_DET_DISTANCE = 15.75
 
-    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, dead_time=False,
-                 paralyzable=True, dead_time_value=4.2,
-                 dead_time_tof_step=100, use_emission_time=True):
-
+    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,
+        dead_time=False,
+        paralyzable=True,
+        dead_time_value=4.2,
+        dead_time_tof_step=100,
+        use_emission_time=True,
+    ):
         if instrument in [self.INSTRUMENT_4A, self.INSTRUMENT_4B]:
             self.instrument = instrument
         else:
@@ -285,13 +308,13 @@ def __init__(self, scattering_workspace, direct_workspace,
 
         # Process workspaces
         if self.tof_range is not None:
-            self._ws_sc = api.CropWorkspace(InputWorkspace=scattering_workspace,
-                                            XMin=tof_range[0], XMax=tof_range[1],
-                                            OutputWorkspace='_'+str(scattering_workspace))
+            self._ws_sc = api.CropWorkspace(
+                InputWorkspace=scattering_workspace, XMin=tof_range[0], XMax=tof_range[1], OutputWorkspace="_" + str(scattering_workspace)
+            )
             if direct_workspace is not None:
-                self._ws_db = api.CropWorkspace(InputWorkspace=direct_workspace,
-                                                XMin=tof_range[0], XMax=tof_range[1],
-                                                OutputWorkspace='_'+str(direct_workspace))
+                self._ws_db = api.CropWorkspace(
+                    InputWorkspace=direct_workspace, XMin=tof_range[0], XMax=tof_range[1], OutputWorkspace="_" + str(direct_workspace)
+                )
             else:
                 self._ws_db = None
         else:
@@ -303,7 +326,7 @@ def __init__(self, scattering_workspace, direct_workspace,
 
     def extract_meta_data(self):
         """
-            Extract meta data from the data file.
+        Extract meta data from the loaded data file.
         """
         # Get instrument parameters
         settings = read_settings(self._ws_sc)
@@ -324,12 +347,12 @@ def extract_meta_data(self):
         if self.tof_range is None:
             self.wl_range = get_wl_range(self._ws_sc)
         else:
-            self.wl_range = [self.tof_range[0] / self.constant, self.tof_range[1] /  self.constant]
+            self.wl_range = [self.tof_range[0] / self.constant, self.tof_range[1] / self.constant]
 
         # q_min and q_max are the boundaries for the final Q binning
         # We also hold on to the true Q range covered by the measurement
-        self.q_min_meas = 4.0*np.pi/self.wl_range[1] * np.fabs(np.sin(self.theta))
-        self.q_max_meas = 4.0*np.pi/self.wl_range[0] * np.fabs(np.sin(self.theta))
+        self.q_min_meas = 4.0 * np.pi / self.wl_range[1] * np.fabs(np.sin(self.theta))
+        self.q_max_meas = 4.0 * np.pi / self.wl_range[0] * np.fabs(np.sin(self.theta))
 
         if self.q_min is None:
             self.q_min = self.q_min_meas
@@ -341,29 +364,29 @@ def extract_meta_data(self):
 
         # Catch options that can be turned off
         if self.signal_low_res is None:
-            self.signal_low_res = [1, self.n_x-1]
+            self.signal_low_res = [1, self.n_x - 1]
         if self.norm_low_res is None:
-            self.norm_low_res = [1, self.n_x-1]
+            self.norm_low_res = [1, self.n_x - 1]
 
     def extract_meta_data_4A(self):
         """
-            4A-specific meta data
+        4A-specific meta data
         """
         run_object = self._ws_sc.getRun()
-        self.det_distance = run_object['SampleDetDis'].getStatistics().mean
-        source_sample_distance = run_object['ModeratorSamDis'].getStatistics().mean
-        if run_object['SampleDetDis'].units not in ['m', 'meter']:
+        self.det_distance = run_object["SampleDetDis"].getStatistics().mean
+        source_sample_distance = run_object["ModeratorSamDis"].getStatistics().mean
+        if run_object["SampleDetDis"].units not in ["m", "meter"]:
             self.det_distance /= 1000.0
-        if run_object['ModeratorSamDis'].units not in ['m', 'meter']:
+        if run_object["ModeratorSamDis"].units not in ["m", "meter"]:
             source_sample_distance /= 1000.0
         self.source_detector_distance = source_sample_distance + self.det_distance
 
     def extract_meta_data_4B(self):
         """
-            4B-specific meta data
+        4B-specific meta data
 
-            Distance from source to sample was 13.63 meters prior to the source
-            to detector distance being determined with Bragg edges to be 15.75 m.
+        Distance from source to sample was 13.63 meters prior to the source
+        to detector distance being determined with Bragg edges to be 15.75 m.
         """
         settings = read_settings(self._ws_sc)
 
@@ -390,6 +413,11 @@ def extract_meta_data_4B(self):
             self.source_detector_distance = self.DEFAULT_4B_SOURCE_DET_DISTANCE
 
     def __repr__(self):
+        """
+        Generate a string representation of the reduction settings.
+
+        :return: (str) string representation of the reduction settings
+        """
         output = "Reduction settings:\n"
         output += "    sample-det: %s\n" % self.det_distance
         output += "    source-det: %s\n" % self.source_detector_distance
@@ -403,12 +431,14 @@ def __repr__(self):
 
     def to_dict(self):
         """
-            Returns meta-data to be used/stored.
+        Returns meta-data to be used/stored.
+
+        :return: (dict) dictionary with meta-data
         """
         if self._ws_sc.getRun().hasProperty("start_time"):
             start_time = self._ws_sc.getRun().getProperty("start_time").value
         else:
-            start_time = 'live'
+            start_time = "live"
         experiment = self._ws_sc.getRun().getProperty("experiment_identifier").value
         run_number = self._ws_sc.getRun().getProperty("run_number").value
         sequence_number = int(self._ws_sc.getRun().getProperty("sequence_number").value[0])
@@ -421,25 +451,42 @@ def to_dict(self):
             norm_run = 0
 
         dq0 = 0
-        return dict(wl_min=self.wl_range[0], wl_max=self.wl_range[1],
-                    q_min=self.q_min_meas, q_max=self.q_max_meas, theta=self.theta,
-                    start_time=start_time, experiment=experiment, run_number=run_number,
-                    run_title=run_title, norm_run=norm_run, time=time.ctime(),
-                    dq0=dq0, dq_over_q=self.dq_over_q, sequence_number=sequence_number,
-                    sequence_id=sequence_id, q_summing=self.q_summing)
-
-    def specular(self, q_summing=False, tof_weighted=False, bck_in_q=False,
-                 clean=False, normalize=True):
+        return dict(
+            wl_min=self.wl_range[0],
+            wl_max=self.wl_range[1],
+            q_min=self.q_min_meas,
+            q_max=self.q_max_meas,
+            theta=self.theta,
+            start_time=start_time,
+            experiment=experiment,
+            run_number=run_number,
+            run_title=run_title,
+            norm_run=norm_run,
+            time=time.ctime(),
+            dq0=dq0,
+            dq_over_q=self.dq_over_q,
+            sequence_number=sequence_number,
+            sequence_id=sequence_id,
+            q_summing=self.q_summing,
+        )
+
+    def specular(self, q_summing=False, tof_weighted=False, bck_in_q=False, clean=False, normalize=True):
         """
-            Compute specular reflectivity.
+        Compute specular reflectivity.
+
+        For constant-Q binning, it's preferred to use tof_weighted=True.
+
+        :param q_summing: turns on constant-Q binning
+        :param tof_weighted: if True, binning will be done by weighting each event to the DB distribution
+        :param bck_in_q: if True, the background will be estimated in Q space using the constant-Q binning approach
+        :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
 
-            For constant-Q binning, it's preferred to use tof_weighted=True.
+        :return: A tuple containing:
 
-            :param q_summing: turns on constant-Q binning
-            :param tof_weighted: if True, binning will be done by weighting each event to the DB distribution
-            :param bck_in_q: if True, the background will be estimated in Q space using the constant-Q binning approach
-            :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
+            - q_bins: The Q bin boundaries
+            - refl: The reflectivity values
+            - d_refl: The uncertainties in the reflectivity values
         """
         if tof_weighted:
             self.specular_weighted(q_summing=q_summing, bck_in_q=bck_in_q)
@@ -447,7 +494,7 @@ def specular(self, q_summing=False, tof_weighted=False, bck_in_q=False,
             self.specular_unweighted(q_summing=q_summing, normalize=normalize)
 
         # Remove leading zeros
-        r = np.trim_zeros(self.refl, 'f')
+        r = np.trim_zeros(self.refl, "f")
         trim = len(self.refl) - len(r)
         self.refl = self.refl[trim:]
         self.d_refl = self.d_refl[trim:]
@@ -471,15 +518,30 @@ def specular(self, q_summing=False, tof_weighted=False, bck_in_q=False,
 
     def specular_unweighted(self, q_summing=False, normalize=True):
         """
-            Simple specular reflectivity calculation. This is the same approach as the
-            original LR reduction, which sums up pixels without constant-Q binning.
-            The original approach bins in TOF, then rebins the final results after
-            transformation to Q. This approach bins directly to Q.
+        Simple specular reflectivity calculation. This is the same approach as the
+        original LR reduction, which sums up pixels without constant-Q binning.
+        The original approach bins in TOF, then rebins the final results after
+        transformation to Q. This approach bins directly to Q.
+
+        :param q_summing: (bool, optional) If True, sum the data in Q-space (default is False).
+        :param normalize: (bool, optional) If True, normalize the reflectivity by the direct 
+            beam (default is True).
+
+        :return: A tuple containing:
+
+            - q_bins: The Q bin boundaries
+            - refl: The reflectivity values
+            - d_refl: The uncertainties in the reflectivity values
         """
         # Scattering data
-        refl, d_refl = self._reflectivity(self._ws_sc, peak_position=self.specular_pixel,
-                                          peak=self.signal_peak, low_res=self.signal_low_res,
-                                          theta=self.theta, q_summing=q_summing)
+        refl, d_refl = self._reflectivity(
+            self._ws_sc,
+            peak_position=self.specular_pixel,
+            peak=self.signal_peak,
+            low_res=self.signal_low_res,
+            theta=self.theta,
+            q_summing=q_summing,
+        )
 
         # Remove background
         if self.signal_bck is not None:
@@ -493,19 +555,21 @@ def specular_unweighted(self, q_summing=False, normalize=True):
             # we can bin the DB according to the same transform instead of binning and dividing in TOF.
             # This is mathematically equivalent and convenient in terms of abstraction for later
             # use for the constant-Q calculation elsewhere in the code.
-            norm, d_norm = self._reflectivity(self._ws_db, peak_position=0,
-                                              peak=self.norm_peak, low_res=self.norm_low_res,
-                                              theta=self.theta, q_summing=False)
+            norm, d_norm = self._reflectivity(
+                self._ws_db, peak_position=0, peak=self.norm_peak, low_res=self.norm_low_res, theta=self.theta, q_summing=False
+            )
 
             # Direct beam background could be added here. The effect will be negligible.
             if self.norm_bck is not None:
                 norm_bck, d_norm_bck = self.norm_bck_subtraction()
                 norm -= norm_bck
                 d_norm = np.sqrt(d_norm**2 + d_norm_bck**2)
-            db_bins = norm>0
+            db_bins = norm > 0
 
-            refl[db_bins] = refl[db_bins]/norm[db_bins]
-            d_refl[db_bins] = np.sqrt(d_refl[db_bins]**2 / norm[db_bins]**2 + refl[db_bins]**2 * d_norm[db_bins]**2 / norm[db_bins]**4)
+            refl[db_bins] = refl[db_bins] / norm[db_bins]
+            d_refl[db_bins] = np.sqrt(
+                d_refl[db_bins] ** 2 / norm[db_bins] ** 2 + refl[db_bins] ** 2 * d_norm[db_bins] ** 2 / norm[db_bins] ** 4
+            )
 
             # Hold on to normalization to be able to diagnose issues later
             self.norm = norm[db_bins]
@@ -523,25 +587,40 @@ def specular_unweighted(self, q_summing=False, normalize=True):
 
     def specular_weighted(self, q_summing=True, bck_in_q=False):
         """
-            Compute reflectivity by weighting each event by flux.
-            This allows for summing in Q and to estimate the background in either Q
-            or pixels next to the peak.
+        Compute reflectivity by weighting each event by flux.
+        This allows for summing in Q and to estimate the background in either Q
+        or pixels next to the peak.
+
+        :param q_summing: (bool, optional) If True, sum the data in Q-space (default is False).
+        :param bck_in_q: (bool, optional) If True, subtract background along Q lines (default is False).
+
+        :return: A tuple containing:
+
+            - q_bins: The Q bin boundaries
+            - refl: The reflectivity values
+            - d_refl: The uncertainties in the reflectivity values
         """
         # Event weights for normalization
         db_charge = self._ws_db.getRun().getProtonCharge()
         wl_events, wl_weights = self._get_events(self._ws_db, self.norm_peak, self.norm_low_res)
         wl_dist, wl_bins = np.histogram(wl_events, bins=100, weights=wl_weights)
         _bin_width = wl_bins[1:] - wl_bins[:-1]
-        wl_dist = wl_dist/db_charge/_bin_width
-        wl_middle = [(wl_bins[i+1]+wl_bins[i])/2.0 for i in range(len(wl_bins)-1)]
-
-        refl, d_refl = self._reflectivity(self._ws_sc, peak_position=self.specular_pixel,
-                                          peak=self.signal_peak, low_res=self.signal_low_res,
-                                          theta=self.theta, q_summing=q_summing, wl_dist=wl_dist, wl_bins=wl_middle)
+        wl_dist = wl_dist / db_charge / _bin_width
+        wl_middle = [(wl_bins[i + 1] + wl_bins[i]) / 2.0 for i in range(len(wl_bins) - 1)]
+
+        refl, d_refl = self._reflectivity(
+            self._ws_sc,
+            peak_position=self.specular_pixel,
+            peak=self.signal_peak,
+            low_res=self.signal_low_res,
+            theta=self.theta,
+            q_summing=q_summing,
+            wl_dist=wl_dist,
+            wl_bins=wl_middle,
+        )
 
         if self.signal_bck is not None:
-            refl_bck, d_refl_bck = self.bck_subtraction(wl_dist=wl_dist, wl_bins=wl_middle,
-                                                        q_summing=bck_in_q)
+            refl_bck, d_refl_bck = self.bck_subtraction(wl_dist=wl_dist, wl_bins=wl_middle, q_summing=bck_in_q)
             refl -= refl_bck
             d_refl = np.sqrt(d_refl**2 + d_refl_bck**2)
 
@@ -551,28 +630,34 @@ def specular_weighted(self, q_summing=True, bck_in_q=False):
 
     def _roi_integration(self, ws, peak, low_res, q_bins=None, wl_dist=None, wl_bins=None, q_summing=False):
         """
-            Integrate a region of interest and normalize by the number of included pixels.
+        Integrate a region of interest and normalize by the number of included pixels.
 
-            The options are the same as for the reflectivity calculation.
-            If wl_dist and wl_bins are supplied, the events will be weighted by flux.
-            If q_summing is True, the angle of each neutron will be recalculated according to
-            their position on the detector and place in the proper Q bin.
+        The options are the same as for the reflectivity calculation.
+        If wl_dist and wl_bins are supplied, the events will be weighted by flux.
+        If q_summing is True, the angle of each neutron will be recalculated according to
+        their position on the detector and place in the proper Q bin.
         """
         q_bins = self.q_bins if q_bins is None else q_bins
-        refl_bck, d_refl_bck = self._reflectivity(ws, peak_position=0, q_bins=q_bins,
-                                                  peak=peak, low_res=low_res,
-                                                  theta=self.theta, q_summing=q_summing,
-                                                  wl_dist=wl_dist, wl_bins=wl_bins)
-
-        _pixel_area = (peak[1]-peak[0]+1.0)
+        refl_bck, d_refl_bck = self._reflectivity(
+            ws,
+            peak_position=0,
+            q_bins=q_bins,
+            peak=peak,
+            low_res=low_res,
+            theta=self.theta,
+            q_summing=q_summing,
+            wl_dist=wl_dist,
+            wl_bins=wl_bins,
+        )
+
+        _pixel_area = peak[1] - peak[0] + 1.0
         refl_bck /= _pixel_area
         d_refl_bck /= _pixel_area
         return refl_bck, d_refl_bck
 
-    def bck_subtraction(self, normalize_to_single_pixel=False, q_bins=None, wl_dist=None, wl_bins=None,
-                        q_summing=False):
+    def bck_subtraction(self, normalize_to_single_pixel=False, q_bins=None, wl_dist=None, wl_bins=None, q_summing=False):
         """
-            Higher-level call for background subtraction. Hides the ranges needed to define the ROI.
+        Higher-level call for background subtraction. Hides the ranges needed to define the ROI.
         """
         # Sanity check
         if len(self.signal_bck) == 2 and self.use_functional_bck:
@@ -582,63 +667,87 @@ def bck_subtraction(self, normalize_to_single_pixel=False, q_bins=None, wl_dist=
             self.use_functional_bck = False
 
         if self.use_functional_bck:
-            return background.functional_background(self._ws_sc, self, self.signal_peak,
-                                                    self.signal_bck, self.signal_low_res,
-                                                    normalize_to_single_pixel=normalize_to_single_pixel,
-                                                    q_bins=q_bins, wl_dist=wl_dist, wl_bins=wl_bins,
-                                                    q_summing=q_summing)
+            return background.functional_background(
+                self._ws_sc,
+                self,
+                self.signal_peak,
+                self.signal_bck,
+                self.signal_low_res,
+                normalize_to_single_pixel=normalize_to_single_pixel,
+                q_bins=q_bins,
+                wl_dist=wl_dist,
+                wl_bins=wl_bins,
+                q_summing=q_summing,
+            )
         else:
-            return background.side_background(self._ws_sc, self, self.signal_peak, self.signal_bck,
-                                              self.signal_low_res,
-                                              normalize_to_single_pixel=normalize_to_single_pixel,
-                                              q_bins=q_bins, wl_dist=wl_dist, wl_bins=wl_bins,
-                                              q_summing=q_summing)
+            return background.side_background(
+                self._ws_sc,
+                self,
+                self.signal_peak,
+                self.signal_bck,
+                self.signal_low_res,
+                normalize_to_single_pixel=normalize_to_single_pixel,
+                q_bins=q_bins,
+                wl_dist=wl_dist,
+                wl_bins=wl_bins,
+                q_summing=q_summing,
+            )
 
     def norm_bck_subtraction(self):
         """
-            Higher-level call for background subtraction for the normalization run.
+        Higher-level call for background subtraction for the normalization run.
         """
-        return background.side_background(self._ws_db, self, self.norm_peak, self.norm_bck,
-                                          self.norm_low_res, normalize_to_single_pixel=False)
-
-    def slice(self, x_min=0.002, x_max=0.004, x_bins=None, z_bins=None, # noqa A003
-              refl=None, d_refl=None, normalize=False):
+        return background.side_background(
+            self._ws_db, self, self.norm_peak, self.norm_bck, self.norm_low_res, normalize_to_single_pixel=False
+        )
+
+    def slice(
+        self,
+        x_min=0.002,
+        x_max=0.004,
+        x_bins=None,
+        z_bins=None,  # noqa A003
+        refl=None,
+        d_refl=None,
+        normalize=False,
+    ):
         """
-            Retrieve a slice from the off-specular data.
+        Retrieve a slice from the off-specular data.
         """
         x_bins = self._offspec_x_bins if x_bins is None else x_bins
         z_bins = self._offspec_z_bins if z_bins is None else z_bins
         refl = self._offspec_refl if refl is None else refl
         d_refl = self._offspec_d_refl if d_refl is None else d_refl
 
-        i_min = len(x_bins[x_bins<x_min])
-        i_max = len(x_bins[x_bins<x_max])
+        i_min = len(x_bins[x_bins < x_min])
+        i_max = len(x_bins[x_bins < x_max])
 
         _spec = np.sum(refl[i_min:i_max], axis=0)
-        _d_spec = np.sum( (d_refl[i_min:i_max])**2, axis=0)
+        _d_spec = np.sum((d_refl[i_min:i_max]) ** 2, axis=0)
         _d_spec = np.sqrt(_d_spec)
         if normalize:
-            _spec /= (i_max-i_min)
-            _d_spec /= (i_max-i_min)
+            _spec /= i_max - i_min
+            _d_spec /= i_max - i_min
 
         return z_bins, _spec, _d_spec
 
-    def _reflectivity(self, ws, peak_position, peak, low_res, theta, q_bins=None, q_summing=False,
-                      wl_dist=None, wl_bins=None, sum_pixels=True):
+    def _reflectivity(
+        self, ws, peak_position, peak, low_res, theta, q_bins=None, q_summing=False, wl_dist=None, wl_bins=None, sum_pixels=True
+    ):
         """
-            Assumes that the input workspace is normalized by proton charge.
+        Assumes that the input workspace is normalized by proton charge.
         """
         charge = ws.getRun().getProtonCharge()
         _q_bins = self.q_bins if q_bins is None else q_bins
 
-        shape = len(_q_bins)-1 if sum_pixels else ((peak[1]-peak[0]+1), len(_q_bins)-1)
+        shape = len(_q_bins) - 1 if sum_pixels else ((peak[1] - peak[0] + 1), len(_q_bins) - 1)
         refl = np.zeros(shape)
         d_refl_sq = np.zeros(shape)
         counts = np.zeros(shape)
         _pixel_width = self.pixel_width if q_summing else 0.0
 
-        for i in range(low_res[0], int(low_res[1]+1)):
-            for j in range(peak[0], int(peak[1]+1)):
+        for i in range(low_res[0], int(low_res[1] + 1)):
+            for j in range(peak[0], int(peak[1] + 1)):
                 if self.instrument == self.INSTRUMENT_4A:
                     pixel = j * self.n_y + i
                 else:
@@ -660,14 +769,14 @@ def _reflectivity(self, ws, peak_position, peak, low_res, theta, q_bins=None, q_
                 delta_theta_f = np.arctan(x_distance / self.det_distance) / 2.0
 
                 # Sign will depend on reflect up or down
-                ths_value = ws.getRun()['ths'].value[-1]
+                ths_value = ws.getRun()["ths"].value[-1]
                 delta_theta_f *= np.sign(ths_value)
 
-                qz=4.0*np.pi/wl_list * np.sin(theta + delta_theta_f - d_theta)
+                qz = 4.0 * np.pi / wl_list * np.sin(theta + delta_theta_f - d_theta)
                 qz = np.fabs(qz)
 
                 if wl_dist is not None and wl_bins is not None:
-                    wl_weights = 1.0/np.interp(wl_list, wl_bins, wl_dist, np.inf, np.inf)
+                    wl_weights = 1.0 / np.interp(wl_list, wl_bins, wl_dist, np.inf, np.inf)
                     hist_weights = wl_weights * qz / wl_list
                     hist_weights *= event_weights
                     _counts, _ = np.histogram(qz, bins=_q_bins, weights=hist_weights)
@@ -676,14 +785,14 @@ def _reflectivity(self, ws, peak_position, peak, low_res, theta, q_bins=None, q_
                         refl += _counts
                         counts += _norm
                     else:
-                        refl[j-peak[0]] += _counts
-                        counts[j-peak[0]] += _norm
+                        refl[j - peak[0]] += _counts
+                        counts[j - peak[0]] += _norm
                 else:
                     _counts, _ = np.histogram(qz, bins=_q_bins, weights=event_weights)
                     if sum_pixels:
                         refl += _counts
                     else:
-                        refl[j-peak[0]] += _counts
+                        refl[j - peak[0]] += _counts
 
         # The following is for information purposes
         if q_summing:
@@ -692,9 +801,9 @@ def _reflectivity(self, ws, peak_position, peak, low_res, theta, q_bins=None, q_
             delta_theta_f0 = np.arctan(x0 / self.det_distance) / 2.0
             delta_theta_f1 = np.arctan(x1 / self.det_distance) / 2.0
 
-            qz_max = 4.0*np.pi/self.tof_range[1]*self.constant * np.fabs(np.sin(theta + delta_theta_f0))
-            qz_min = 4.0*np.pi/self.tof_range[1]*self.constant * np.fabs(np.sin(theta + delta_theta_f1))
-            mid_point = (qz_max + qz_min)/2.0
+            qz_max = 4.0 * np.pi / self.tof_range[1] * self.constant * np.fabs(np.sin(theta + delta_theta_f0))
+            qz_min = 4.0 * np.pi / self.tof_range[1] * self.constant * np.fabs(np.sin(theta + delta_theta_f1))
+            mid_point = (qz_max + qz_min) / 2.0
             print("Qz range: ", qz_min, mid_point, qz_max)
             self.summing_threshold = mid_point
 
@@ -705,7 +814,7 @@ def _reflectivity(self, ws, peak_position, peak, low_res, theta, q_bins=None, q_
             if not sum_pixels:
                 bin_size = np.tile(bin_size, [counts.shape[0], 1])
 
-            d_refl_sq[non_zero] =  refl[non_zero] / np.sqrt(counts[non_zero]) / charge / bin_size[non_zero]
+            d_refl_sq[non_zero] = refl[non_zero] / np.sqrt(counts[non_zero]) / charge / bin_size[non_zero]
             refl[non_zero] = refl[non_zero] / charge / bin_size[non_zero]
         else:
             d_refl_sq = np.sqrt(np.fabs(refl)) / charge
@@ -715,13 +824,13 @@ def _reflectivity(self, ws, peak_position, peak, low_res, theta, q_bins=None, q_
 
     def _get_events(self, ws, peak, low_res):
         """
-            Return an array of wavelengths for a given workspace.
+        Return an array of wavelengths for a given workspace.
         """
         wl_events = np.asarray([])
         wl_weights = np.asarray([])
 
-        for i in range(low_res[0], int(low_res[1]+1)):
-            for j in range(peak[0], int(peak[1]+1)):
+        for i in range(low_res[0], int(low_res[1] + 1)):
+            for j in range(peak[0], int(peak[1] + 1)):
                 if self.instrument == self.INSTRUMENT_4A:
                     pixel = j * self.n_y + i
                 else:
@@ -737,17 +846,16 @@ def _get_events(self, ws, peak, low_res):
                 wl_weights = np.concatenate((wl_weights, weights))
         return wl_events, wl_weights
 
-    def off_specular(self, x_axis=None, x_min=-0.015, x_max=0.015, x_npts=50,
-                     z_min=None, z_max=None, z_npts=-120, bck_in_q=None):
+    def off_specular(self, x_axis=None, x_min=-0.015, x_max=0.015, x_npts=50, z_min=None, z_max=None, z_npts=-120, bck_in_q=None):
         """
-            Compute off-specular
-            :param x_axis: Axis selection
-            :param x_min: Min value on x-axis
-            :param x_max: Max value on x-axis
-            :param x_npts: Number of points in x (negative will produce a log scale)
-            :param z_min: Min value on z-axis (if none, default Qz will be used)
-            :param z_max: Max value on z-axis (if none, default Qz will be used)
-            :param z_npts: Number of points in z (negative will produce a log scale)
+        Compute off-specular
+        :param x_axis: Axis selection
+        :param x_min: Min value on x-axis
+        :param x_max: Max value on x-axis
+        :param x_npts: Number of points in x (negative will produce a log scale)
+        :param z_min: Min value on z-axis (if none, default Qz will be used)
+        :param z_max: Max value on z-axis (if none, default Qz will be used)
+        :param z_npts: Number of points in z (negative will produce a log scale)
         """
         # Z axis binning
         qz_bins = self.q_bins
@@ -765,23 +873,23 @@ def off_specular(self, x_axis=None, x_min=-0.015, x_max=0.015, x_npts=50,
 
         wl_events, wl_weights = self._get_events(self._ws_db, self.norm_peak, self.norm_low_res)
         wl_dist, wl_bins = np.histogram(wl_events, bins=100, weights=wl_weights)
-        wl_middle = [(wl_bins[i+1]+wl_bins[i])/2.0 for i in range(len(wl_bins)-1)]
+        wl_middle = [(wl_bins[i + 1] + wl_bins[i]) / 2.0 for i in range(len(wl_bins) - 1)]
 
-        _refl, _d_refl = self._off_specular(self._ws_sc, wl_dist, wl_middle, qx_bins, qz_bins,
-                                            self.specular_pixel, self.theta, x_axis=x_axis)
+        _refl, _d_refl = self._off_specular(
+            self._ws_sc, wl_dist, wl_middle, qx_bins, qz_bins, self.specular_pixel, self.theta, x_axis=x_axis
+        )
         db_charge = self._ws_db.getRun().getProtonCharge()
-        _refl *= db_charge * (wl_bins[1]-wl_bins[0])
-        _d_refl *= db_charge * (wl_bins[1]-wl_bins[0])
+        _refl *= db_charge * (wl_bins[1] - wl_bins[0])
+        _d_refl *= db_charge * (wl_bins[1] - wl_bins[0])
 
         # Background
         if self.signal_bck:
             if bck_in_q is None:
                 print("Not implemented")
             else:
-                _, refl_bck, d_refl_bck = self.slice(bck_in_q[0], bck_in_q[1],
-                                                     x_bins=qx_bins, z_bins=qz_bins,
-                                                     refl=_refl, d_refl=_d_refl,
-                                                     normalize=True)
+                _, refl_bck, d_refl_bck = self.slice(
+                    bck_in_q[0], bck_in_q[1], x_bins=qx_bins, z_bins=qz_bins, refl=_refl, d_refl=_d_refl, normalize=True
+                )
                 _refl -= refl_bck
                 _d_refl = np.sqrt(_d_refl**2 + d_refl_bck**2)
 
@@ -794,12 +902,12 @@ def off_specular(self, x_axis=None, x_min=-0.015, x_max=0.015, x_npts=50,
 
     def _off_specular(self, ws, wl_dist, wl_bins, x_bins, z_bins, peak_position, theta, x_axis=None):
         charge = ws.getRun().getProtonCharge()
-        refl = np.zeros([len(x_bins)-1, len(z_bins)-1])
-        counts = np.zeros([len(x_bins)-1, len(z_bins)-1])
+        refl = np.zeros([len(x_bins) - 1, len(z_bins) - 1])
+        counts = np.zeros([len(x_bins) - 1, len(z_bins) - 1])
 
         for j in range(0, self.n_x):
             wl_list = np.asarray([])
-            for i in range(self.signal_low_res[0], int(self.signal_low_res[1]+1)):
+            for i in range(self.signal_low_res[0], int(self.signal_low_res[1] + 1)):
                 if self.instrument == self.INSTRUMENT_4A:
                     pixel = j * self.n_y + i
                 else:
@@ -809,12 +917,12 @@ def _off_specular(self, ws, wl_dist, wl_bins, x_bins, z_bins, peak_position, the
                 wl_list = np.concatenate((wl_events, wl_list))
 
             k = 2.0 * np.pi / wl_list
-            wl_weights = 1.0/np.interp(wl_list, wl_bins, wl_dist, np.inf, np.inf)
+            wl_weights = 1.0 / np.interp(wl_list, wl_bins, wl_dist, np.inf, np.inf)
 
-            x_distance = float(j-peak_position) * self.pixel_width
+            x_distance = float(j - peak_position) * self.pixel_width
             delta_theta_f = np.arctan(x_distance / self.det_distance)
             # Sign will depend on reflect up or down
-            ths_value = ws.getRun()['ths'].value[-1]
+            ths_value = ws.getRun()["ths"].value[-1]
             delta_theta_f *= np.sign(ths_value)
             theta_f = theta + delta_theta_f
 
@@ -827,13 +935,13 @@ def _off_specular(self, ws, wl_dist, wl_bins, x_bins, z_bins, peak_position, the
             _x = qx
             _z = qz
             if x_axis == EventReflectivity.DELTA_KZ_VS_QZ:
-                _x = (ki_z - kf_z)
+                _x = ki_z - kf_z
             elif x_axis == EventReflectivity.KZI_VS_KZF:
                 _x = ki_z
                 _z = kf_z
             elif x_axis == EventReflectivity.THETAF_VS_WL:
                 _x = wl_list
-                _z = np.ones(len(wl_list))*theta_f
+                _z = np.ones(len(wl_list)) * theta_f
 
             histo_weigths = wl_weights * _z / wl_list
             _counts, _, _ = np.histogram2d(_x, _z, bins=[x_bins, z_bins], weights=histo_weigths)
@@ -842,16 +950,16 @@ def _off_specular(self, ws, wl_dist, wl_bins, x_bins, z_bins, peak_position, the
             counts += _counts
 
         bin_size = z_bins[1] - z_bins[0]
-        d_refl_sq =  refl / np.sqrt(counts) / charge / bin_size
+        d_refl_sq = refl / np.sqrt(counts) / charge / bin_size
         refl /= charge * bin_size
 
         return refl, d_refl_sq
 
     def emission_time_correction(self, ws, tofs):
         """
-            Coorect TOF for emission time delay in the moderator
-            :param ws: Mantid workspace
-            :param tofs: list of TOF values
+        Coorect TOF for emission time delay in the moderator
+        :param ws: Mantid workspace
+        :param tofs: list of TOF values
         """
         mt_run = ws.getRun()
         use_emission_delay = False
@@ -867,7 +975,7 @@ def emission_time_correction(self, ws, tofs):
 
     def gravity_correction(self, ws, wl_list):
         """
-            Gravity correction for each event
+        Gravity correction for each event
         """
         # Xi reference would be the position of xi if the si slit were to be positioned
         # at the sample. The distance from the sample to si is then xi_reference - xi.
@@ -878,7 +986,7 @@ def gravity_correction(self, ws, wl_list):
         # Distance between the s1 and the sample
         s1_sample_distance = 1485
         if ws.getInstrument().hasParameter("s1-sample-distance"):
-            s1_sample_distance = ws.getInstrument().getNumberParameter("s1-sample-distance")[0]*1000
+            s1_sample_distance = ws.getInstrument().getNumberParameter("s1-sample-distance")[0] * 1000
 
         xi = 310
         if ws.getInstrument().hasParameter("BL4B:Mot:xi.RBV"):
@@ -888,48 +996,48 @@ def gravity_correction(self, ws, wl_list):
         slit_distance = s1_sample_distance - sample_si_distance
 
         # Angle of the incident beam on a horizontal sample
-        theta_in=-4.0
+        theta_in = -4.0
 
         # Calculation from the ILL paper. This works for inclined beams.
         # Calculated theta is the angle on the sample
 
-        g = 9.8067                        # m/s^2
-        h = 6.6260715e-34                 # Js=kg m^2/s
-        mn = 1.67492749804e-27            # kg
+        g = 9.8067  # m/s^2
+        h = 6.6260715e-34  # Js=kg m^2/s
+        mn = 1.67492749804e-27  # kg
 
-        v = h/(mn*wl_list*1e-10)
-        k = g/(2*v**2)
+        v = h / (mn * wl_list * 1e-10)
+        k = g / (2 * v**2)
 
         # Define the sample position as x=0, y=0. increasing x is towards moderator
-        xs=0
+        xs = 0
 
         # positions of slits
         x1 = sample_si_distance / 1000
         x2 = (sample_si_distance + slit_distance) / 1000
 
-        #height of slits determined by incident theta, y=0 is the sample height
-        y1=x1*np.tan(theta_in*np.pi/180)
-        y2=x2*np.tan(theta_in*np.pi/180)
+        # height of slits determined by incident theta, y=0 is the sample height
+        y1 = x1 * np.tan(theta_in * np.pi / 180)
+        y2 = x2 * np.tan(theta_in * np.pi / 180)
 
         # This is the location of the top of the parabola
-        x0=(y1-y2+k*(x1**2-x2**2))/(2*k*(x1-x2))
+        x0 = (y1 - y2 + k * (x1**2 - x2**2)) / (2 * k * (x1 - x2))
 
         # Angle is arctan(dy/dx) at sample
-        theta_sample = np.arctan(2*k*(x0-xs)) * 180/np.pi
+        theta_sample = np.arctan(2 * k * (x0 - xs)) * 180 / np.pi
 
-        return (theta_sample-theta_in) * np.pi / 180.0
+        return (theta_sample - theta_in) * np.pi / 180.0
 
 
 def compute_resolution(ws, default_dq=0.027, theta=None, q_summing=False):
     """
-        Compute the Q resolution from the meta data.
-        :param theta: scattering angle in radians
-        :param q_summing: if True, the pixel size will be used for the resolution
+    Compute the Q resolution from the meta data.
+    :param theta: scattering angle in radians
+    :param q_summing: if True, the pixel size will be used for the resolution
     """
     settings = read_settings(ws)
 
     if theta is None:
-        theta = abs(ws.getRun().getProperty("ths").value[0]) * np.pi / 180.
+        theta = abs(ws.getRun().getProperty("ths").value[0]) * np.pi / 180.0
 
     if q_summing:
         # Q resolution is reported as FWHM, so here we consider this to be
diff --git a/reduction/lr_reduction/workflow.py b/reduction/lr_reduction/workflow.py
index d51e143..449c8f2 100644
--- a/reduction/lr_reduction/workflow.py
+++ b/reduction/lr_reduction/workflow.py
@@ -1,6 +1,7 @@
 """
-    Autoreduction process for the Liquids Reflectometer
+Autoreduction process for the Liquids Reflectometer
 """
+
 import json
 import os
 
@@ -10,21 +11,21 @@
 from . import event_reduction, output, reduction_template_reader, template
 
 
-def reduce(ws, template_file, output_dir, average_overlap=False,
-           theta_offset: float | None = 0,
-           q_summing=False, bck_in_q=False, is_live=False):
+def reduce(
+    ws, template_file, output_dir, average_overlap=False, theta_offset: float | None = 0, q_summing=False, bck_in_q=False, is_live=False
+):
     """
-        Function called by reduce_REFL.py, which lives in /SNS/REF_L/shared/autoreduce
-        and is called by the automated reduction workflow.
+    Function called by reduce_REFL.py, which lives in /SNS/REF_L/shared/autoreduce
+    and is called by the automated reduction workflow.
 
-        If average_overlap is used, overlapping points will be averaged, otherwise they
-        will be left in the final data file.
+    If average_overlap is used, overlapping points will be averaged, otherwise they
+    will be left in the final data file.
 
-        :param average_overlap: if True, the overlapping points will be averaged
-        :param q_summing: if True, constant-Q binning will be used
-        :param bck_in_q: if True, and constant-Q binning is used, the background will be estimated
-                         along constant-Q lines rather than along TOF/pixel boundaries.
-        :param theta_offset: Theta offset to apply. If None, the template value will be used.
+    :param average_overlap: if True, the overlapping points will be averaged
+    :param q_summing: if True, constant-Q binning will be used
+    :param bck_in_q: if True, and constant-Q binning is used, the background will be estimated
+                     along constant-Q lines rather than along TOF/pixel boundaries.
+    :param theta_offset: Theta offset to apply. If None, the template value will be used.
     """
     # Get the sequence number
     sequence_number = 1
@@ -42,12 +43,9 @@ def reduce(ws, template_file, output_dir, average_overlap=False,
         template_data.angle_offset = theta_offset
 
     # Call the reduction using the template
-    qz_mid, refl, d_refl, meta_data = template.process_from_template_ws(ws, template_data,
-                                                                        q_summing=q_summing,
-                                                                        tof_weighted=bck_in_q,
-                                                                        clean=q_summing,
-                                                                        bck_in_q=bck_in_q,
-                                                                        info=True)
+    qz_mid, refl, d_refl, meta_data = template.process_from_template_ws(
+        ws, template_data, q_summing=q_summing, tof_weighted=bck_in_q, clean=q_summing, bck_in_q=bck_in_q, info=True
+    )
 
     # Save partial results
     coll = output.RunCollection()
@@ -55,16 +53,15 @@ def reduce(ws, template_file, output_dir, average_overlap=False,
 
     # If this is live data, put it in a separate file to avoid conflict with auto-reduction
     if is_live:
-        reduced_file = os.path.join(output_dir, 'REFL_live_partial.txt')
+        reduced_file = os.path.join(output_dir, "REFL_live_partial.txt")
     else:
-        reduced_file = os.path.join(output_dir, 'REFL_%s_%s_%s_partial.txt' % (meta_data['sequence_id'],
-                                                                               meta_data['sequence_number'],
-                                                                               meta_data['run_number']))
+        reduced_file = os.path.join(
+            output_dir, "REFL_%s_%s_%s_partial.txt" % (meta_data["sequence_id"], meta_data["sequence_number"], meta_data["run_number"])
+        )
     coll.save_ascii(reduced_file, meta_as_json=True)
 
     # Assemble partial results into a single R(q)
-    seq_list, run_list = assemble_results(meta_data['sequence_id'], output_dir,
-                                          average_overlap, is_live=is_live)
+    seq_list, run_list = assemble_results(meta_data["sequence_id"], output_dir, average_overlap, is_live=is_live)
 
     # Save template. This will not happen if the template_file input was
     # template data, which the template processing allows.
@@ -79,7 +76,7 @@ def reduce(ws, template_file, output_dir, average_overlap=False,
 
 def assemble_results(first_run, output_dir, average_overlap=False, is_live=False):
     """
-        Find related runs and assemble them in one R(q) data set
+    Find related runs and assemble them in one R(q) data set
     """
     # Keep track of sequence IDs and run numbers so we can make a new template
     seq_list = []
@@ -88,8 +85,8 @@ def assemble_results(first_run, output_dir, average_overlap=False, is_live=False
 
     file_list = sorted(os.listdir(output_dir))
     for item in file_list:
-        if item.startswith("REFL_%s" % first_run) and item.endswith('partial.txt'):
-            toks = item.split('_')
+        if item.startswith("REFL_%s" % first_run) and item.endswith("partial.txt"):
+            toks = item.split("_")
             if not len(toks) == 5 or int(toks[2]) == 0:
                 continue
             seq_list.append(int(toks[2]))
@@ -97,14 +94,14 @@ def assemble_results(first_run, output_dir, average_overlap=False, is_live=False
 
             # Read the partial data and add to a collection
             _, _, _, _, _meta = output.read_file(os.path.join(output_dir, item))
-            if is_live or not _meta['start_time'] == "live":
+            if is_live or not _meta["start_time"] == "live":
                 coll.add_from_file(os.path.join(output_dir, item))
-        #elif item == "REFL_live_partial.txt":
+        # elif item == "REFL_live_partial.txt":
         #    coll.add_from_file(os.path.join(output_dir, item))
 
-    output_file_name = 'REFL_%s_combined_data_auto.txt' % first_run
+    output_file_name = "REFL_%s_combined_data_auto.txt" % first_run
     if is_live:
-        output_file_name = 'REFL_%s_live_estimate.txt' % first_run
+        output_file_name = "REFL_%s_live_estimate.txt" % first_run
     coll.save_ascii(os.path.join(output_dir, output_file_name))
 
     return seq_list, run_list
@@ -112,8 +109,8 @@ def assemble_results(first_run, output_dir, average_overlap=False, is_live=False
 
 def write_template(seq_list, run_list, template_file, output_dir):
     """
-        Read the appropriate entry in a template file and save an updated
-        copy with the updated run number.
+    Read the appropriate entry in a template file and save an updated
+    copy with the updated run number.
     """
     with open(template_file, "r") as fd:
         xml_str = fd.read()
@@ -122,22 +119,22 @@ def write_template(seq_list, run_list, template_file, output_dir):
         new_data_sets = []
         for i in range(len(seq_list)):
             if len(data_sets) >= seq_list[i]:
-                data_sets[seq_list[i]-1].data_files = [run_list[i]]
-                new_data_sets.append(data_sets[seq_list[i]-1])
+                data_sets[seq_list[i] - 1].data_files = [run_list[i]]
+                new_data_sets.append(data_sets[seq_list[i] - 1])
             else:
                 print("Too few entries [%s] in template for sequence number %s" % (len(data_sets), seq_list[i]))
 
     # Save the template that was used
     xml_str = reduction_template_reader.to_xml(new_data_sets)
-    with open(os.path.join(output_dir, 'REF_L_%s_auto_template.xml' % run_list[0]), 'w') as fd:
+    with open(os.path.join(output_dir, "REF_L_%s_auto_template.xml" % run_list[0]), "w") as fd:
         fd.write(xml_str)
 
 
 def offset_from_first_run(
-        ws,
-        template_file :str,
-        output_dir: str,
-        ):
+    ws,
+    template_file: str,
+    output_dir: str,
+):
     """
     Find a theta offset from the first peak.
     Used when sample is misaligned.
@@ -153,7 +150,7 @@ def offset_from_first_run(
         sequence_id = ws.getRun().getProperty("sequence_id").value[0]
 
     # Theta value that we are aiming for
-    ths_value = ws.getRun()['ths'].value[-1]
+    ths_value = ws.getRun()["ths"].value[-1]
 
     # Read template so we can load the direct beam run
     template_data = template.read_template(template_file, sequence_number)
@@ -163,31 +160,30 @@ def offset_from_first_run(
     ws_db = mtd_api.LoadEventNexus("REF_L_%s" % template_data.norm_file)
 
     # Look for parameters that might have been determined earlier for this measurement
-    options_file = os.path.join(output_dir, 'REFL_%s_options.json' % sequence_id)
+    options_file = os.path.join(output_dir, "REFL_%s_options.json" % sequence_id)
 
     if sequence_number > 1 and os.path.isfile(options_file):
-        with open(options_file, 'r') as fd:
+        with open(options_file, "r") as fd:
             options = json.load(fd)
-            return options['theta_offset']
+            return options["theta_offset"]
     else:
         # Fit direct beam position
-        x_min=template_data.norm_peak_range[0]
-        x_max=template_data.norm_peak_range[1]
+        x_min = template_data.norm_peak_range[0]
+        x_max = template_data.norm_peak_range[1]
         _, _x, _y = peak_finding.process_data(ws_db, summed=True, tof_step=200)
         peak_center = np.argmax(_y)
-        db_center, db_width, _ = peak_finding.fit_signal_flat_bck(_x, _y, x_min=x_min, x_max=x_max,
-                                                           center=peak_center,
-                                                           sigma=1.)
-        print("    DB center: %g\t Width: %g from [%g %g]" % (db_center, db_width,
-                                                              template_data.norm_peak_range[0],
-                                                              template_data.norm_peak_range[1]))
+        db_center, db_width, _ = peak_finding.fit_signal_flat_bck(_x, _y, x_min=x_min, x_max=x_max, center=peak_center, sigma=1.0)
+        print(
+            "    DB center: %g\t Width: %g from [%g %g]"
+            % (db_center, db_width, template_data.norm_peak_range[0], template_data.norm_peak_range[1])
+        )
 
         # Fit the reflected beam position
-        x_min=template_data.data_peak_range[0]
-        x_max=template_data.data_peak_range[1]
+        x_min = template_data.data_peak_range[0]
+        x_max = template_data.data_peak_range[1]
         _, _x, _y = peak_finding.process_data(ws, summed=True, tof_step=200)
         peak_center = np.argmax(_y[x_min:x_max]) + x_min
-        sc_center, sc_width, _ = peak_finding.fit_signal_flat_bck(_x, _y, x_min=x_min, x_max=x_max, center=peak_center, sigma=3.)
+        sc_center, sc_width, _ = peak_finding.fit_signal_flat_bck(_x, _y, x_min=x_min, x_max=x_max, center=peak_center, sigma=3.0)
         pixel_offset = sc_center - peak_center
         print("    SC center: %g\t Width: %g" % (sc_center, sc_width))
 
@@ -195,36 +191,36 @@ def offset_from_first_run(
         sample_det_distance = settings["sample-det-distance"]
         pixel_width = settings["pixel-width"] / 1000.0
 
-        theta = np.arctan((sc_center-db_center) * pixel_width / sample_det_distance) / 2.0 * 180 / np.pi
+        theta = np.arctan((sc_center - db_center) * pixel_width / sample_det_distance) / 2.0 * 180 / np.pi
         theta_offset = theta - ths_value
 
         # If this is the first angle, keep the value for later
-        options = dict(theta_offset = theta_offset,
-                       pixel_offset = pixel_offset)
-        with open(options_file, 'w') as fp:
+        options = dict(theta_offset=theta_offset, pixel_offset=pixel_offset)
+        with open(options_file, "w") as fp:
             json.dump(options, fp)
 
     return theta_offset
 
 
 def reduce_explorer(ws, ws_db, theta_pv=None, center_pixel=145, db_center_pixel=145, peak_width=10):
-    """
-    """
+    """ """
     from . import peak_finding
 
     if theta_pv is None:
-        if 'BL4B:CS:ExpPl:OperatingMode' in ws.getRun() \
-            and ws.getRun().getProperty('BL4B:CS:ExpPl:OperatingMode').value[0] == 'Free Liquid':
-            theta_pv = 'thi'
+        if (
+            "BL4B:CS:ExpPl:OperatingMode" in ws.getRun()
+            and ws.getRun().getProperty("BL4B:CS:ExpPl:OperatingMode").value[0] == "Free Liquid"
+        ):
+            theta_pv = "thi"
         else:
-            theta_pv = 'ths'
+            theta_pv = "ths"
     print("\nProcessing: %s" % ws.getRunNumber())
 
     # Theta value that we are aiming for
     theta_value = np.fabs(ws.getRun()[theta_pv].value[0])
 
     # Load normalization run
-    tthd_value = ws.getRun()['tthd'].value[0]
+    tthd_value = ws.getRun()["tthd"].value[0]
 
     # Fit direct beam position
     x_min = center_pixel - 25
@@ -235,8 +231,8 @@ def reduce_explorer(ws, ws_db, theta_pv=None, center_pixel=145, db_center_pixel=
     print("    DB center: %g\t Width: %g" % (db_center, db_width))
 
     # Fit the reflected beam position
-    x_min=db_center_pixel-peak_width
-    x_max=db_center_pixel+peak_width
+    x_min = db_center_pixel - peak_width
+    x_max = db_center_pixel + peak_width
     tof, _x, _y = peak_finding.process_data(ws, summed=True, tof_step=200)
     peak_center = np.argmax(_y)
     sc_center, sc_width, _ = peak_finding.fit_signal_flat_bck(_x, _y, x_min=x_min, x_max=x_max, center=peak_center)
@@ -244,41 +240,47 @@ def reduce_explorer(ws, ws_db, theta_pv=None, center_pixel=145, db_center_pixel=
 
     pixel_width = float(ws.getInstrument().getNumberParameter("pixel-width")[0]) / 1000.0
     sample_det_distance = event_reduction.EventReflectivity.DEFAULT_4B_SAMPLE_DET_DISTANCE
-    twotheta = np.arctan((db_center-sc_center)*pixel_width / sample_det_distance) / 2.0 * 180 / np.pi
+    twotheta = np.arctan((db_center - sc_center) * pixel_width / sample_det_distance) / 2.0 * 180 / np.pi
 
     # Store the tthd of the direct beam and account for the fact that it may be
     # different from our reflected beam for this calibration data.
     # This will allow us to be compatible with both fixed and moving detector arm.
-    tthd_db = ws_db.getRun()['tthd'].value[0]
+    tthd_db = ws_db.getRun()["tthd"].value[0]
     twotheta = twotheta + tthd_value - tthd_db
 
     print("    Theta = %g   Two-theta = %g" % (theta_value, twotheta))
 
     # Perform the reduction
     width_mult = 2.5
-    peak = [np.rint(sc_center - width_mult*sc_width).astype(int), np.rint(sc_center + width_mult*sc_width).astype(int)]
-    norm_peak = [np.rint(db_center - width_mult*db_width).astype(int), np.rint(db_center + width_mult*db_width).astype(int)]
-    peak_bck = [peak[0]-3, peak[1]+3]
-    norm_bck = [norm_peak[0]-3, norm_peak[1]+3]
+    peak = [np.rint(sc_center - width_mult * sc_width).astype(int), np.rint(sc_center + width_mult * sc_width).astype(int)]
+    norm_peak = [np.rint(db_center - width_mult * db_width).astype(int), np.rint(db_center + width_mult * db_width).astype(int)]
+    peak_bck = [peak[0] - 3, peak[1] + 3]
+    norm_bck = [norm_peak[0] - 3, norm_peak[1] + 3]
 
     tof_min = ws.getTofMin()
     tof_max = ws.getTofMax()
 
-    theta = theta_value * np.pi / 180.
+    theta = theta_value * np.pi / 180.0
 
-    #TODO: dead time correction should be applied here
-    event_refl = event_reduction.EventReflectivity(ws, ws_db,
-                                                   signal_peak=peak, signal_bck=peak_bck,
-                                                   norm_peak=norm_peak, norm_bck=norm_bck,
-                                                   specular_pixel=sc_center.value,
-                                                   signal_low_res=[65,180], norm_low_res=[65,180],
-                                                   q_min=None, q_max=None,
-                                                   tof_range = [tof_min, tof_max],
-                                                   theta=theta)
+    # TODO: dead time correction should be applied here
+    event_refl = event_reduction.EventReflectivity(
+        ws,
+        ws_db,
+        signal_peak=peak,
+        signal_bck=peak_bck,
+        norm_peak=norm_peak,
+        norm_bck=norm_bck,
+        specular_pixel=sc_center.value,
+        signal_low_res=[65, 180],
+        norm_low_res=[65, 180],
+        q_min=None,
+        q_max=None,
+        tof_range=[tof_min, tof_max],
+        theta=theta,
+    )
 
     # R(Q)
-    qz, refl, d_refl = event_refl.specular(q_summing=False, tof_weighted=False,
-                                           bck_in_q=False, clean=False, normalize=True)
-    qz_mid = (qz[:-1] + qz[1:])/2.0
+    qz, refl, d_refl = event_refl.specular(q_summing=False, tof_weighted=False, bck_in_q=False, clean=False, normalize=True)
+    qz_mid = (qz[:-1] + qz[1:]) / 2.0
 
     return qz_mid, refl, d_refl