Merge pull request #13259 from mantidproject/7817_tosca_sample_position_correction

TOSCA sample position correction algorithm

Code/Mantid/Framework/PythonInterface/plugins/algorithms/WorkflowAlgorithms/TOSCABankCorrection.py

+#pylint: disable=no-init
+from mantid.kernel import *
+from mantid.api import *
+from mantid.simpleapi import *
+
+
+class TOSCABankCorrection(DataProcessorAlgorithm):
+
+    _input_ws = None
+    _output_ws = None
+    _search_range = None
+    _peak_tolerance = None
+    _peak_position = None
+    _peak_function = None
+
+
+    def category(self):
+        return 'PythonAlgorithms;Inelastic;CorrectionFunctions'
+
+
+    def summary(self):
+        return 'Corrects TOSCA reductions where the peaks across banks are not in alignment.'
+
+
+    def PyInit(self):
+        self.declareProperty(WorkspaceProperty(name='InputWorkspace', defaultValue='',
+                             direction=Direction.Input),
+                             doc='Input reduced workspace')
+
+        self.declareProperty(FloatArrayProperty(name='SearchRange',
+                                                values=[200, 2000]),
+                             doc='Range over which to find peaks')
+
+        self.declareProperty(name='PeakPosition', defaultValue='',
+                             doc='Specify a particular peak to use')
+
+        self.declareProperty(name='ClosePeakTolerance', defaultValue=20.0,
+                             doc='Tolerance under which peaks are considered to be the same')
+
+        self.declareProperty(name='PeakFunction', defaultValue='Lorentzian',
+                             validator=StringListValidator(['Lorentzian', 'Gaussian']),
+                             doc='Type of peak to search for')
+
+        self.declareProperty(MatrixWorkspaceProperty(name='OutputWorkspace', defaultValue='',
+                             direction=Direction.Output),
+                             doc='Output corrected workspace')
+
+        self.declareProperty(name='TargetPeakCentre', defaultValue=0.0,
+                             direction=Direction.Output,
+                             doc='X position between the centres of the two ' \
+                                 'selected peaks')
+
+        self.declareProperty(name='ScaleFactor1', defaultValue=1.0,
+                             direction=Direction.Output,
+                             doc='Scale factor for the first bank (histogram 0)')
+
+        self.declareProperty(name='ScaleFactor2', defaultValue=1.0,
+                             direction=Direction.Output,
+                             doc='Scale factor for the second bank (histogram 1)')
+
+
+    def validateInputs(self):
+        self._get_properties()
+        issues = dict()
+
+        # Validate search range
+        if len(self._search_range) != 2:
+            issues['SearchRange'] = 'Search range must have two values'
+        elif self._search_range[0] > self._search_range[1]:
+            issues['SearchRange'] = 'Search range must be in format "low,high"'
+
+        # Ensure manual peak position is inside search range
+        if 'SearchRange' not in issues and self._peak_position is not None:
+            if self._peak_position < self._search_range[0] or self._peak_position > self._search_range[1]:
+                issues['PeakPosition'] = 'Peak position must be inside SearchRange'
+
+        return issues
+
+
+    def PyExec(self):
+        self._get_properties()
+
+        # Crop the sample workspace to the search range
+        CropWorkspace(InputWorkspace=self._input_ws,
+                      OutputWorkspace='__search_ws',
+                      XMin=self._search_range[0],
+                      XMax=self._search_range[1])
+
+        peak = self._get_peak('__search_ws')
+        DeleteWorkspace('__search_ws')
+
+        # Ensure there is at least one peak found
+        if peak is None:
+            raise RuntimeError('Could not find any peaks. Try increasing \
+                                width of SearchRange and/or ClosePeakTolerance')
+
+        target_centre = (peak[0] + peak[1]) / 2.0
+        bank_1_scale_factor = target_centre / peak[0]
+        bank_2_scale_factor = target_centre / peak[1]
+
+        logger.information('Bank 1 scale factor: %f' % bank_1_scale_factor)
+        logger.information('Bank 2 scale factor: %f' % bank_2_scale_factor)
+
+        self._apply_correction(bank_1_scale_factor, bank_2_scale_factor)
+
+        self.setPropertyValue('OutputWorkspace', self._output_ws)
+
+        self.setProperty('TargetPeakCentre', target_centre)
+
+        self.setProperty('ScaleFactor1', bank_1_scale_factor)
+        self.setProperty('ScaleFactor2', bank_2_scale_factor)
+
+
+    def _get_properties(self):
+        self._input_ws = self.getPropertyValue('InputWorkspace')
+        self._output_ws = self.getPropertyValue('OutputWorkspace')
+
+        self._search_range = self.getProperty('SearchRange').value
+        self._peak_tolerance = self.getProperty('ClosePeakTolerance').value
+
+        self._peak_function = self.getPropertyValue('PeakFunction')
+
+        try:
+            self._peak_position = float(self.getPropertyValue('PeakPosition'))
+        except ValueError:
+            self._peak_position = None
+
+
+    def _get_peak(self, search_ws):
+        """
+        Finds a matching peak over the two banks.
+
+        @param search_ws Workspace to search
+        @return Peak centres for matching peak over both banks
+        """
+
+        find_peak_args = dict()
+        if self._peak_position is not None:
+            find_peak_args['PeakPositions'] = [self._peak_position]
+
+        # Find the peaks in each bank
+        FindPeaks(InputWorkspace=search_ws,
+                  PeaksList='__bank_1_peaks',
+                  WorkspaceIndex=0,
+                  PeakFunction=self._peak_function,
+                  **find_peak_args)
+
+        FindPeaks(InputWorkspace=search_ws,
+                  PeaksList='__bank_2_peaks',
+                  WorkspaceIndex=1,
+                  PeakFunction=self._peak_function,
+                  **find_peak_args)
+
+        # Sort peaks by height, prefer to match tall peaks
+        SortTableWorkspace(InputWorkspace='__bank_1_peaks',
+                           OutputWorkspace='__bank_1_peaks',
+                           Columns='centre',
+                           Ascending=False)
+
+        bank_1_ws = mtd['__bank_1_peaks']
+        bank_2_ws = mtd['__bank_2_peaks']
+
+        matching_peaks = list()
+
+        # Find the centres of two peaks that are close to each other on both banks
+        for peak_idx in range(0, bank_1_ws.rowCount()):
+            bank_1_centre = bank_1_ws.cell('centre', peak_idx)
+
+            for other_peak_idx in range(0, bank_2_ws.rowCount()):
+                bank_2_centre = bank_2_ws.cell('centre', other_peak_idx)
+
+                if abs(bank_1_centre - bank_2_centre) < self._peak_tolerance:
+                    matching_peaks.append((bank_1_centre, bank_2_centre))
+
+        # Remove temporary workspaces
+        DeleteWorkspace('__bank_1_peaks')
+        DeleteWorkspace('__bank_2_peaks')
+
+        if len(matching_peaks) > 0:
+            selected_peak = matching_peaks[0]
+            logger.debug('Found matching peak: %s' % (str(selected_peak)))
+        else:
+            selected_peak = None
+            logger.warning('No peak found')
+
+        return selected_peak
+
+
+    def _apply_correction(self, bank_1_sf, bank_2_sf):
+        """
+        Applies correction to a copy of the input workspace.
+
+        @param bank_1_sf Bank 1 scale factor
+        @param bank_2_sf Bank 2 scale factor
+        """
+
+        # Create a copy of the workspace with just bank spectra
+        CropWorkspace(InputWorkspace=self._input_ws,
+                      OutputWorkspace=self._output_ws,
+                      StartWorkspaceIndex=0,
+                      EndWorkspaceIndex=1)
+
+        # Scale bank 1 X axis
+        ScaleX(InputWorkspace=self._output_ws,
+               OutputWorkspace=self._output_ws,
+               Operation='Multiply',
+               Factor=bank_1_sf,
+               IndexMin=0,
+               IndexMax=0)
+
+        # Scale bank 2 X axis
+        ScaleX(InputWorkspace=self._output_ws,
+               OutputWorkspace=self._output_ws,
+               Operation='Multiply',
+               Factor=bank_2_sf,
+               IndexMin=1,
+               IndexMax=1)
+
+        # Rebin the corrected workspace to match the input
+        RebinToWorkspace(WorkspaceToRebin=self._output_ws,
+                         WorkspaceToMatch=self._input_ws,
+                         OutputWorkspace=self._output_ws)
+
+        # Recalculate the sum spectra from corrected banks
+        GroupDetectors(InputWorkspace=self._output_ws,
+                       OutputWorkspace='__sum',
+                       SpectraList=[0,1],
+                       Behaviour='Average')
+
+        # Add the new sum spectra to the output workspace
+        AppendSpectra(InputWorkspace1=self._output_ws,
+                      InputWorkspace2='__sum',
+                      OutputWorkspace=self._output_ws)
+
+        DeleteWorkspace('__sum')
+
+
+AlgorithmFactory.subscribe(TOSCABankCorrection)

Code/Mantid/Framework/PythonInterface/test/python/plugins/algorithms/CMakeLists.txt

@@ -74,6 +74,7 @@ set ( TEST_PY_FILES
   UpdatePeakParameterTableValueTest.py
   SANSSubtractTest.py
   TimeSliceTest.py
+  TOSCABankCorrectionTest.py
   TransformToIqtTest.py
   ExportSampleLogsToCSVFileTest.py
   ExportExperimentLogTest.py

Code/Mantid/Framework/PythonInterface/test/python/plugins/algorithms/TOSCABankCorrectionTest.py

+import unittest
+from mantid.simpleapi import *
+from mantid.api import *
+
+
+class TOSCABankCorrectionTest(unittest.TestCase):
+
+    def setUp(self):
+        """
+        Loads sample workspace.
+        """
+
+        self._original = '__TOSCABankCorrectionTest_original'
+
+        Load(Filename='TSC14007_graphite002_red.nxs',
+             OutputWorkspace=self._original)
+
+
+    def tearDown(self):
+        """
+        Removes workspaces.
+        """
+
+        DeleteWorkspace(self._original)
+
+
+    def test_automatic_peak_selection(self):
+        """
+        Tests automatically finding a peak in the default range.
+        """
+
+        corrected_reduction, peak_position, scale_factor_1, scale_factor_2 = \
+          TOSCABankCorrection(InputWorkspace=self._original)
+
+        self.assertAlmostEqual(peak_position, 1077.47222328)
+        self.assertAlmostEqual(scale_factor_1, 1.0059271)
+        self.assertAlmostEqual(scale_factor_2, 0.9941423)
+
+
+    def test_automatic_peak_in_range(self):
+        """
+        Tests automatically finding a peak in a given range.
+        """
+
+        corrected_reduction, peak_position, scale_factor_1, scale_factor_2 = \
+          TOSCABankCorrection(InputWorkspace=self._original,
+                              SearchRange=[200, 800])
+
+        self.assertAlmostEqual(peak_position, 713.20080359)
+        self.assertAlmostEqual(scale_factor_1, 1.006076146)
+        self.assertAlmostEqual(scale_factor_2, 0.993996806)
+
+
+    def test_manual_peak_selection(self):
+        """
+        Tests using a peak provided by the user.
+        """
+
+        corrected_reduction, peak_position, scale_factor_1, scale_factor_2 = \
+          TOSCABankCorrection(InputWorkspace=self._original,
+                              PeakPosition='715')
+
+        self.assertAlmostEqual(peak_position, 713.4430671)
+        self.assertAlmostEqual(scale_factor_1, 1.00611439)
+        self.assertAlmostEqual(scale_factor_2, 0.99395947)
+
+
+    def test_manual_peak_not_found(self):
+        """
+        Tests error handling when a peak cannot be found using a manual peak position.
+        """
+
+        self.assertRaises(RuntimeError,
+                          TOSCABankCorrection,
+                          InputWorkspace=self._original,
+                          OutputWorkspace='__TOSCABankCorrectionTest_output',
+                          PeakPosition='900')
+
+
+    def test_validation_search_range_order(self):
+        """
+        Tests validation to ensure low and high values are entered in correct order.
+        """
+
+        self.assertRaises(RuntimeError,
+                          TOSCABankCorrection,
+                          InputWorkspace=self._original,
+                          OutputWorkspace='__TOSCABankCorrectionTest_output',
+                          SearchRange=[500, 50])
+
+
+    def test_validation_search_range_count(self):
+        """
+        Tests validation to ensure two values exist values are entered in correct order.
+        """
+
+        self.assertRaises(RuntimeError,
+                          TOSCABankCorrection,
+                          InputWorkspace=self._original,
+                          OutputWorkspace='__TOSCABankCorrectionTest_output',
+                          SearchRange=[500])
+
+
+    def test_validation_peak_position_in_search_range(self):
+        """
+        Tests validation to ensure that the PeakPosition falls inside SearchRange.
+        """
+
+        self.assertRaises(RuntimeError,
+                          TOSCABankCorrection,
+                          InputWorkspace=self._original,
+                          OutputWorkspace='__TOSCABankCorrectionTest_output',
+                          SearchRange=[200, 2000],
+                          PeakPosition='2500')
+
+
+if __name__ == '__main__':
+    unittest.main()

Code/Mantid/Testing/Data/DocTest/TSC14007_graphite002_red.nxs.md5

+360b5639a909d42b099693b8d7a99b45

Code/Mantid/Testing/Data/UnitTest/TSC14007_graphite002_red.nxs.md5

+360b5639a909d42b099693b8d7a99b45

Code/Mantid/docs/source/algorithms/IndirectResolution-v1.rst

@@ -23,6 +23,8 @@ Workflow
 Usage
 -----
 
+.. include:: ../usagedata-note.txt
+
 **Example - Running IndirectResolution.**
 
 .. testcode:: ExIndirectResolutionSimple

Code/Mantid/docs/source/algorithms/IndirectTransmissionMonitor-v1.rst

@@ -27,6 +27,8 @@ Workflow
 Usage
 -----
 
+.. include:: ../usagedata-note.txt
+
 **Example - Create mapping file for IRIS**
 
 .. testcode:: exIRISTransmission

Code/Mantid/docs/source/algorithms/TOSCABankCorrection-v1.rst

+.. algorithm::
+
+.. summary::
+
+.. alias::
+
+.. properties::
+
+Description
+-----------
+
+This algorithm attempts to automatically correct TOSCA data in which the
+position of the sample has been moved and has affected the alignment of features
+on the spectra from forward and backscattering detector banks.
+
+The input workspace should be an energy transfer reduction, for the default
+values of SearchRange and ClosePeakTolerance the X axis is assumed to be in
+cm-1, however the X axis is not restricted to this unit.
+
+The algorithm works by finding peaks of a given shape (using the :ref:`FindPeaks
+<algm-FindPeaks>`) on both the forward and backscattering banks, either
+selecting a peak in a given position or selecting the peak with the highest X
+value and attempting to match them to what is believed to be the same feature on
+the other bank.
+
+A scale factor is then calculated for each bank that will align at least the
+selected peak and in doing so will also align the majority of misaligned peaks
+across the two banks.
+
+The sacling factor is calculated as follows:
+
+.. math::
+
+  X_{centre} = \frac{X_{forward peak} + X_{back peak}}{2}
+
+  SF_{forward} = \frac{X_{centre}}{X_{forward peak}}
+
+  SF_{back} = \frac{X_{centre}}{X_{back peak}}
+
+The corrected spectra are then rebinned to the input workspace (using
+:ref:`RebinToWorkspace <algm-RebinToWorkspace>`) to preserve the X range and to
+maintain bin alignment.
+
+The sum spectra (containing both forward and back scattering detectors) is then
+recalculated by averaging the intensities of the two corrected spectra, this
+compensates for the broader peaks seen on the original sum spectra due to the
+misalignment of the peaks.
+
+.. note::
+  This algorithm is only intended to provide an approximation of what the
+  measured spectra would look like if the sample was in the expected sample
+  position.
+
+Usage
+-----
+
+.. include:: ../usagedata-note.txt
+
+**Example - Automatic peak selection.**
+
+.. testcode:: ExTOSCABankCorrectionAutomatic
+
+    original_reduction = Load('TSC14007_graphite002_red.nxs')
+
+    corrected_reduction, peak_position, scale_factor_1, scale_factor_2 = \
+      TOSCABankCorrection(InputWorkspace=original_reduction)
+
+    print 'Target peak centre: %.f' % peak_position
+
+Output:
+
+.. testoutput:: ExTOSCABankCorrectionAutomatic
+
+    Target peak centre: 1077
+
+**Example - Manual peak selection.**
+
+.. testcode:: ExTOSCABankCorrectionManual
+
+    original_reduction = Load('TSC14007_graphite002_red.nxs')
+
+    corrected_reduction, peak_position, scale_factor_1, scale_factor_2 = \
+      TOSCABankCorrection(InputWorkspace=original_reduction,
+                          PeakPosition='715')
+
+    print 'Target peak centre: %.f' % peak_position
+
+Output:
+
+.. testoutput:: ExTOSCABankCorrectionManual
+
+    Target peak centre: 713
+
+.. categories::