diff --git a/Framework/PythonInterface/plugins/algorithms/LRDirectBeamSort.py b/Framework/PythonInterface/plugins/algorithms/LRDirectBeamSort.py
new file mode 100644
index 0000000000000000000000000000000000000000..8060917f7338aecd3f137b4ce4282fb94e25513e
--- /dev/null
+++ b/Framework/PythonInterface/plugins/algorithms/LRDirectBeamSort.py
@@ -0,0 +1,247 @@
+#pylint: disable=no-init,invalid-name
+import math
+import time
+import numpy as np
+import mantid
+from mantid.api import *
+from mantid.simpleapi import *
+from mantid.kernel import *
+import logging
+
+class CompareTwoNXSDataForSFcalculator(object):
+ """
+ will return -1, 0 or 1 according to the position of the nexusToPosition in relation to the
+ nexusToCompareWith based on the following criteria
+ #1: number of attenuators (ascending order)
+ #2: lambda requested (descending order)
+ #3: S2W (ascending order)
+ #4: S2H (descending order)
+ #5 if everything up to this point is identical, return 0
+ """
+ nexusToCompareWithRun = None
+ nexusToPositionRun = None
+ resultComparison = 0
+
+ def __init__(self, nxsdataToCompareWith, nxsdataToPosition):
+ self.nexusToCompareWithRun = nxsdataToCompareWith.getRun()
+ self.nexusToPositionRun = nxsdataToPosition.getRun()
+
+ compare1 = self.compareParameter('LambdaRequest', 'descending')
+ if compare1 != 0:
+ self.resultComparison = compare1
+ return
+
+ compare2 = self.compareParameter('vATT', 'ascending')
+ if compare2 != 0:
+ self.resultComparison = compare2
+ return
+
+ pcharge1 = self.nexusToCompareWithRun.getProperty('gd_prtn_chrg').value/nxsdataToCompareWith.getNEvents()
+ pcharge2 = self.nexusToPositionRun.getProperty('gd_prtn_chrg').value/nxsdataToPosition.getNEvents()
+
+ self.resultComparison = -1 if pcharge1 < pcharge2 else 1
+
+ def compareParameter(self, param, order):
+ """
+ Compare parameters for the two runs
+ """
+ _nexusToCompareWithRun = self.nexusToCompareWithRun
+ _nexusToPositionRun = self.nexusToPositionRun
+
+ _paramNexusToCompareWith = float(_nexusToCompareWithRun.getProperty(param).value[0])
+ _paramNexusToPosition = float(_nexusToPositionRun.getProperty(param).value[0])
+
+ if order == 'ascending':
+ resultLessThan = -1
+ resultMoreThan = 1
+ else:
+ resultLessThan = 1
+ resultMoreThan = -1
+
+ if (_paramNexusToPosition < _paramNexusToCompareWith):
+ return resultLessThan
+ elif (_paramNexusToPosition > _paramNexusToCompareWith):
+ return resultMoreThan
+ else:
+ return 0
+
+ def result(self):
+ return self.resultComparison
+
+def sorter_function(r1, r2):
+ """
+ Sorter function used by with the 'sorted' call to sort the direct beams.
+ """
+ return CompareTwoNXSDataForSFcalculator(r2, r1).result()
+
+
+class LRDirectBeamSort(PythonAlgorithm):
+
+ def category(self):
+ return "Reflectometry\\SNS"
+
+ def name(self):
+ return "LRDirectBeamSort"
+
+ def version(self):
+ return 1
+
+ def summary(self):
+ return "Sort a set of direct beams for the purpose of calculating scaling factors."
+
+ def PyInit(self):
+ self.declareProperty(IntArrayProperty("RunList", [], direction=Direction.Input),
+ "List of run numbers (integers) to be sorted - takes precedence over WorkspaceList")
+ self.declareProperty(StringArrayProperty("WorkspaceList", [], direction=Direction.Input),
+ "List of workspace names to be sorted")
+ self.declareProperty("UseLowResCut", False, direction=Direction.Input,
+ doc="If True, an x-direction cut will be determined and used")
+ self.declareProperty("ComputeScalingFactors", True, direction=Direction.Input,
+ doc="If True, the scaling factors will be computed")
+ self.declareProperty("TOFSteps", 200.0, doc="TOF bin width")
+ self.declareProperty(FileProperty("ScalingFactorFile","",
+ action=FileAction.OptionalSave,
+ extensions=['cfg']),
+ "Scaling factor file to be created")
+
+ self.declareProperty(IntArrayProperty("OrderedRunList", [], direction=Direction.Output),
+ "Ordered list of run numbers")
+ self.declareProperty(StringArrayProperty("OrderedNameList", [], direction=Direction.Output),
+ "Ordered list of workspace names corresponding to the run list")
+
+ def PyExec(self):
+ compute = self.getProperty("ComputeScalingFactors").value
+ lr_data = []
+ run_list = self.getProperty("RunList").value
+ if len(run_list) > 0:
+ for run in run_list:
+ workspace = LoadEventNexus(Filename="REF_L_%s" % run, OutputWorkspace="__data_file_%s" % run, MetaDataOnly=not compute)
+ lr_data.append(workspace)
+ else:
+ ws_list = self.getProperty("WorkspaceList").value
+ for ws in ws_list:
+ lr_data.append(mtd[ws])
+
+ lr_data_sorted = sorted(lr_data, cmp=sorter_function)
+
+ # Set the output properties
+ run_numbers = [r.getRunNumber() for r in lr_data_sorted]
+ ws_names = [str(r) for r in lr_data_sorted]
+ self.setProperty("OrderedRunList", run_numbers)
+ self.setProperty("OrderedNameList", ws_names)
+
+ # Compute the scaling factors if requested
+ if compute:
+ sf_file = self.getProperty("ScalingFactorFile").value
+ if len(sf_file)==0:
+ logger.error("Scaling factors were requested but no output file was set")
+ else:
+ self._compute_scaling_factors(lr_data_sorted)
+
+ def _compute_scaling_factors(self, lr_data_sorted):
+ """
+ If we need to compute the scaling factors, group the runs by their wavelength request
+ @param lr_data_sorted: ordered list of workspaces
+ """
+ #TODO: change LRScalingFactor so that it reads the number of attenuators from the file.
+ group_list = []
+ current_group = []
+ group_wl = None
+ for r in lr_data_sorted:
+ wl = r.getRun().getProperty('LambdaRequest').value[0]
+
+ if not group_wl == wl:
+ # New group
+ group_wl = wl
+ if len(current_group)>0:
+ group_list.append(current_group)
+ current_group = []
+
+ current_group.append(r)
+
+ # Add in the last group
+ group_list.append(current_group)
+
+ tof_steps = self.getProperty("TOFSteps").value
+ scaling_file = self.getProperty("ScalingFactorFile").value
+ use_low_res_cut = self.getProperty("UseLowResCut").value
+ summary = ""
+ for g in group_list:
+ if len(g) == 0:
+ continue
+ runs = [r.getRunNumber() for r in g]
+ logger.notice(str(runs))
+
+ direct_beam_runs = []
+ attenuators = []
+ peak_ranges = []
+ x_ranges = []
+ bck_ranges = []
+
+ for run in g:
+ # Create peak workspace
+ number_of_pixels_x = int(run.getInstrument().getNumberParameter("number-of-x-pixels")[0])
+ number_of_pixels_y = int(run.getInstrument().getNumberParameter("number-of-y-pixels")[0])
+
+ # Direct beam signal
+ workspace = RefRoi(InputWorkspace=run, ConvertToQ=False,
+ NXPixel=number_of_pixels_x,
+ NYPixel=number_of_pixels_y,
+ IntegrateY=False,
+ OutputWorkspace="__ref_peak")
+ workspace = Transpose(InputWorkspace=workspace)
+ peak, _, _ = LRPeakSelection(InputWorkspace=workspace)
+
+ # Low resolution cut
+ if use_low_res_cut:
+ workspace = RefRoi(InputWorkspace=run, ConvertToQ=False,
+ NXPixel=number_of_pixels_x,
+ NYPixel=number_of_pixels_y,
+ IntegrateY=True,
+ OutputWorkspace="__ref_peak")
+ workspace = Transpose(InputWorkspace=workspace)
+ _, low_res, _ = LRPeakSelection(InputWorkspace=workspace)
+ else:
+ low_res = [0, number_of_pixels_x]
+
+ # TOF range
+ att = run.getRun().getProperty('vATT').value[0]-1
+ attenuators.append(int(att))
+
+ direct_beam_runs.append(run.getRunNumber())
+ peak_ranges.append(int(peak[0]))
+ peak_ranges.append(int(peak[1]))
+ x_ranges.append(int(low_res[0]))
+ x_ranges.append(int(low_res[1]))
+ bck_ranges.append(int(peak[0])-3)
+ bck_ranges.append(int(peak[1])+3)
+
+ summary += "%10s %s %5s,%5s %5s,%5s\n" % (run.getRunNumber(), att, peak[0], peak[1], low_res[0], low_res[1])
+
+ # Determine TOF range from first file
+ sample = g[0].getInstrument().getSample()
+ source = g[0].getInstrument().getSource()
+ source_sample_distance = sample.getDistance(source)
+ detector = g[0].getDetector(0)
+ sample_detector_distance = detector.getPos().getZ()
+ source_detector_distance = source_sample_distance + sample_detector_distance
+ h = 6.626e-34 # m^2 kg s^-1
+ m = 1.675e-27 # kg
+ wl = g[0].getRun().getProperty('LambdaRequest').value[0]
+ tof_min = source_detector_distance / h * m * (wl + 0.5 - 1.7) * 1e-4
+ tof_max = source_detector_distance / h * m * (wl + 0.5 + 1.7) * 1e-4
+ tof_range = [tof_min, tof_max]
+
+ summary += "TOF: %s\n\n" % tof_range
+
+ # Compute the scaling factors
+ LRScalingFactors(DirectBeamRuns=direct_beam_runs,
+ Attenuators = attenuators,
+ TOFRange=tof_range, TOFSteps=tof_steps,
+ SignalPeakPixelRange=peak_ranges,
+ SignalBackgroundPixelRange=bck_ranges,
+ LowResolutionPixelRange=x_ranges,
+ ScalingFactorFile=scaling_file)
+ logger.notice(summary)
+
+AlgorithmFactory.subscribe(LRDirectBeamSort)
diff --git a/Framework/PythonInterface/plugins/algorithms/LRPeakSelection.py b/Framework/PythonInterface/plugins/algorithms/LRPeakSelection.py
new file mode 100644
index 0000000000000000000000000000000000000000..21cfd23cdd6d100e54dfc1ac3fc20e7f74bd4582
--- /dev/null
+++ b/Framework/PythonInterface/plugins/algorithms/LRPeakSelection.py
@@ -0,0 +1,252 @@
+#pylint: disable=no-init,invalid-name
+import math
+import time
+import numpy as np
+import mantid
+from mantid.api import *
+from mantid.simpleapi import *
+from mantid.kernel import *
+
+class PeakFinderDerivation(object):
+ """
+ Determine various types of peak for reflectivity.
+ Those include specular peaks and low-resolution direction signal range.
+ """
+ xdata_firstderi = []
+ ydata_firstderi = []
+ five_highest_ydata = []
+ five_highest_xdata = []
+ sum_peak_counts = -1
+ sum_peak_counts_time_pixel = -1
+ peak_pixel = -1
+ deri_min = 1
+ deri_max = -1
+ deri_min_pixel_value = -1
+ deri_max_pixel_value = -1
+ mean_counts_firstderi = -1
+ std_deviation_counts_firstderi = -1
+ peak_max_final_value = -1
+ peak_min_final_value = -1
+
+ def __init__(self, workspace, back_offset=4):
+ self.back_offset = back_offset
+ self.ydata = workspace.dataY(0)
+ self.xdata = np.arange(len(self.ydata))
+
+ self.compute()
+
+ def compute(self):
+ """
+ Perform the computation
+ """
+ self.initArrays()
+ self.calculate_five_highest_points()
+ self.calculate_peak_pixel()
+ self.calculate_first_derivative()
+ self.calculate_min_max_derivative_pixels()
+ self.calculate_avg_and_derivative()
+ self.calculate_min_max_signal_pixels()
+ self.low_resolution_range()
+
+ def initArrays(self):
+ """
+ Initialize internal data members
+ """
+ self.xdata_firstderi = []
+ self.ydata_firstderi = []
+ self.peak = [-1, -1]
+ self.low_res = [-1, -1]
+ self.five_highest_ydata = []
+ self.five_highest_xdata = []
+ self.sum_five_highest_ydata = -1
+ self.peak_pixel = -1
+ self.deri_min = 1
+ self.deri_max = -1
+ self.deri_min_pixel_value = -1
+ self.deri_max_pixel_value = -1
+ self.mean_counts_firstderi = -1
+ self.std_deviation_counts_firstderi = -1
+ self.peak_max_final_value = -1
+ self.peak_min_final_value = -1
+
+ def calculate_five_highest_points(self):
+ _xdata = self.xdata
+ _ydata = self.ydata
+
+ _sort_ydata = np.sort(_ydata)
+ _decreasing_sort_ydata = _sort_ydata[::-1]
+ self.five_highest_ydata = _decreasing_sort_ydata[0:5]
+
+ _sort_index = np.argsort(_ydata)
+ _decreasing_sort_index = _sort_index[::-1]
+ _5decreasing_sort_index = _decreasing_sort_index[0:5]
+ self.five_highest_xdata = _xdata[_5decreasing_sort_index]
+
+ def calculate_peak_pixel(self):
+ self.sum_peak_counts = sum(self.five_highest_ydata)
+ _sum_peak_counts_time_pixel = -1
+ for index,yvalue in enumerate(self.five_highest_ydata):
+ _sum_peak_counts_time_pixel += yvalue * self.five_highest_xdata[index]
+ self.sum_peak_counts_time_pixel = _sum_peak_counts_time_pixel
+ self.peak_pixel = round(self.sum_peak_counts_time_pixel / self.sum_peak_counts)
+
+ def calculate_first_derivative(self):
+ xdata = self.xdata
+ ydata = self.ydata
+
+ _xdata_firstderi = []
+ _ydata_firstderi = []
+ for i in range(len(xdata)-1):
+ _left_x = xdata[i]
+ _right_x = xdata[i+1]
+ _xdata_firstderi.append(np.mean([_left_x, _right_x]))
+
+ _left_y = ydata[i]
+ _right_y = ydata[i+1]
+ _ydata_firstderi.append((_right_y - _left_y)/(_right_x - _left_x))
+
+ self.xdata_firstderi = _xdata_firstderi
+ self.ydata_firstderi = _ydata_firstderi
+
+ def calculate_min_max_derivative_pixels(self):
+ _pixel = self.xdata_firstderi
+ _counts_firstderi = self.ydata_firstderi
+
+ _sort_counts_firstderi = np.sort(_counts_firstderi)
+ self.deri_min = _sort_counts_firstderi[0]
+ self.deri_max = _sort_counts_firstderi[-1]
+
+ _sort_index = np.argsort(_counts_firstderi)
+ self.deri_min_pixel_value = int(min([_pixel[_sort_index[0]], _pixel[_sort_index[-1]]]))
+ self.deri_max_pixel_value = int(max([_pixel[_sort_index[0]], _pixel[_sort_index[-1]]]))
+
+ def calculate_avg_and_derivative(self):
+ _counts_firstderi = np.array(self.ydata_firstderi)
+ self.mean_counts_firstderi = np.mean(_counts_firstderi)
+ _mean_counts_firstderi = np.mean(_counts_firstderi * _counts_firstderi)
+ self.std_deviation_counts_firstderi = math.sqrt(_mean_counts_firstderi)
+
+ def calculate_min_max_signal_pixels(self):
+ """
+ Determine specular peak region
+ """
+ _counts = self.ydata_firstderi
+ _pixel = self.xdata_firstderi
+
+ _deri_min_pixel_value = self.deri_min_pixel_value
+ _deri_max_pixel_value = self.deri_max_pixel_value
+
+ _std_deviation_counts_firstderi = self.std_deviation_counts_firstderi
+
+ px_offset = 0
+ while abs(_counts[int(_deri_min_pixel_value - px_offset)]) > _std_deviation_counts_firstderi:
+ px_offset += 1
+ _peak_min_final_value = _pixel[int(_deri_min_pixel_value - px_offset)]
+
+ px_offset = 0
+ while abs(_counts[int(round(_deri_max_pixel_value + px_offset))]) > _std_deviation_counts_firstderi:
+ px_offset += 1
+ _peak_max_final_value = _pixel[int(round(_deri_max_pixel_value + px_offset))]
+
+ self.peak = [int(_peak_min_final_value), int(np.ceil(_peak_max_final_value))]
+
+ def low_resolution_range(self):
+ """
+ Determine the x range of the signal
+ """
+ y_integrated = []
+ total = 0.0
+ for y_value in self.ydata:
+ total += y_value
+ y_integrated.append(total)
+
+ for i in range(len(y_integrated)):
+ y_integrated[i] /= total
+
+ # Derivative of the flipped integrated distribution
+ deriv = []
+ offset = 1
+ for i in range(offset, len(y_integrated)):
+ value = (self.xdata[i]-self.xdata[i-offset]) / (y_integrated[i]-y_integrated[i-offset])
+ deriv.append(value)
+
+ # Find lower edge of the main peak
+ center = int(len(deriv)/2.0)
+ middle_value = deriv[center]
+ i_min = 0
+ for i in range(center, 0, -1):
+ if deriv[i]/middle_value>3:
+ i_min = i
+ break
+ # Find upper edge of the main peak
+ i_max = len(deriv)
+ for i in range(center, i_max):
+ if deriv[i]/middle_value>3:
+ i_max = i
+ break
+
+ self.low_res = [int(self.xdata[i_min])-self.back_offset, int(self.xdata[i_max])+self.back_offset]
+ return self.low_res
+
+class LRPeakSelection(PythonAlgorithm):
+
+ def category(self):
+ return "Reflectometry\\SNS"
+
+ def name(self):
+ return "LRPeakSelection"
+
+ def version(self):
+ return 1
+
+ def summary(self):
+ return "Find reflectivity peak and return its pixel range."
+
+ def PyInit(self):
+ self.declareProperty(WorkspaceProperty("InputWorkspace", "",Direction.Input), "Workspace to select peak from")
+ self.declareProperty(IntArrayProperty("PeakRange", [0,0], direction=Direction.Output))
+ self.declareProperty(IntArrayProperty("LowResRange", [0,0], direction=Direction.Output))
+ self.declareProperty(IntArrayProperty("PrimaryRange", [0,0], direction=Direction.Output))
+ self.declareProperty("ComputePrimaryRange", False, doc="If True, the primary fraction range will be determined")
+
+ def PyExec(self):
+ workspace = self.getProperty("InputWorkspace").value
+
+ # Main peak finding algorithm
+ pf = PeakFinderDerivation(workspace)
+ self.setProperty("PeakRange", pf.peak)
+ logger.information("Peak: [%s, %s]" % (pf.peak[0], pf.peak[1]))
+ self.setProperty("LowResRange", pf.low_res)
+ logger.information("Low Res: [%s, %s]" % (pf.low_res[0], pf.low_res[1]))
+
+ # Primary fraction range
+ compute_primary = self.getProperty("ComputePrimaryRange").value
+ if compute_primary:
+ primary_range = self.clocking_range(workspace)
+ self.setProperty("PrimaryRange", primary_range)
+ logger.information("Primary: [%s, %s]" % (primary_range[0], primary_range[1]))
+
+ def clocking_range(self, workspace):
+ """
+ Determine the primary fraction range
+ @param workspace: workspace to determine the range from
+ """
+ # Get the full range
+ pf = PeakFinderDerivation(workspace, back_offset=0)
+ [left_max, right_min] = pf.low_res
+
+ # Process left-end data
+ pf.ydata = workspace.dataY(0)[0: left_max]
+ pf.xdata = np.arange(len(pf.ydata))
+ pf.compute()
+ left_clocking = pf.calculate_low_res_range()[0]
+
+ pf.ydata = workspace.dataY(0)[right_min: -1]
+ pf.xdata = np.arange(len(pf.ydata))
+ pf.compute()
+ right_clocking = pf.calculate_low_res_range()[1]
+
+ return [left_clocking, right_clocking]
+
+
+AlgorithmFactory.subscribe(LRPeakSelection)
diff --git a/docs/source/algorithms/LRDirectBeamSort-v1.rst b/docs/source/algorithms/LRDirectBeamSort-v1.rst
new file mode 100644
index 0000000000000000000000000000000000000000..1fa9d3863959e77d94a321add929fd29f35dbb81
--- /dev/null
+++ b/docs/source/algorithms/LRDirectBeamSort-v1.rst
@@ -0,0 +1,17 @@
+.. algorithm::
+
+.. summary::
+
+.. alias::
+
+.. properties::
+
+Description
+-----------
+
+Sort a set of direct beams for the purpose of calculating scaling factors.
+By setting ComputeScalingFactors to True, the scaling factors will be computed.
+
+.. categories::
+
+.. sourcelink::
diff --git a/docs/source/algorithms/LRPeakSelection-v1.rst b/docs/source/algorithms/LRPeakSelection-v1.rst
new file mode 100644
index 0000000000000000000000000000000000000000..c295738b6f502dd6fa58f24de7417a44b428dbdf
--- /dev/null
+++ b/docs/source/algorithms/LRPeakSelection-v1.rst
@@ -0,0 +1,17 @@
+.. algorithm::
+
+.. summary::
+
+.. alias::
+
+.. properties::
+
+Description
+-----------
+
+Determine various types of peak for reflectivity.
+Those include specular peaks and low-resolution direction signal range.
+
+.. categories::
+
+.. sourcelink::
diff --git a/docs/source/algorithms/LRReflectivityOutput-v1.rst b/docs/source/algorithms/LRReflectivityOutput-v1.rst
new file mode 100644
index 0000000000000000000000000000000000000000..dd826d0041a663301f254b4b339fe993016c6f98
--- /dev/null
+++ b/docs/source/algorithms/LRReflectivityOutput-v1.rst
@@ -0,0 +1,17 @@
+.. algorithm::
+
+.. summary::
+
+.. alias::
+
+.. properties::
+
+Description
+-----------
+
+Produces a single reflectivity curve from multiple reflectivity ranges. The algorithm can scale the specular plateau to unity.
+It will also write out the Q resolution.
+
+.. categories::
+
+.. sourcelink::