Merge pull request #13920 from mantidproject/13919_update_detector_flood_weighting

DetectorFloodWeighting algorithm additions for SANS

Framework/PythonInterface/plugins/algorithms/WorkflowAlgorithms/DetectorFloodWeighting.py

 from mantid.api import DataProcessorAlgorithm, AlgorithmFactory, MatrixWorkspaceProperty, WorkspaceUnitValidator, \
-    Progress
+PropertyMode, Progress
 
 from mantid.kernel import Direction, FloatArrayProperty, FloatArrayBoundedValidator
 
@@ -22,7 +22,11 @@ class DetectorFloodWeighting(DataProcessorAlgorithm):
 
         self.declareProperty(MatrixWorkspaceProperty('InputWorkspace', '',
                                                      direction=Direction.Input, validator=WorkspaceUnitValidator("Wavelength")),
-                             doc='Flood weighting measurement')
+                                                     doc='Flood weighting measurement')
+        self.declareProperty(MatrixWorkspaceProperty('TransmissionWorkspace', '',
+                                                     direction=Direction.Input, optional=PropertyMode.Optional,
+                                                     validator=WorkspaceUnitValidator("Wavelength")),
+                                                     doc='Flood weighting measurement')
 
         validator = FloatArrayBoundedValidator()
         validator.setLower(0.)
@@ -52,9 +56,6 @@ class DetectorFloodWeighting(DataProcessorAlgorithm):
             issues['Bands'] = 'Even number of Bands boundaries expected'
             return issues # Abort early. Do not continue
 
-        if len(bands) > 2:
-            issues['Bands'] = 'Presently this algorithm only supports one pair of bands'
-
         all_limits=list()
         for i in range(0, len(bands), 2):
             lower = bands[i]
@@ -69,6 +70,13 @@ class DetectorFloodWeighting(DataProcessorAlgorithm):
             all_limits.append(limits)
             if lower >= upper:
                 issues['Bands'] = 'Bands should form lower, upper pairs'
+        input_ws = self.getProperty('InputWorkspace').value
+        trans_ws = self.getProperty('TransmissionWorkspace').value
+        if trans_ws:
+            if not trans_ws.getNumberHistograms() == input_ws.getNumberHistograms():
+                issues['TransmissionWorkspace'] = 'Transmission should have same number of histograms as flood input workspace'
+            if not trans_ws.blocksize() == input_ws.blocksize():
+                issues['TransmissionWorkspace'] = 'Transmission workspace should be rebinned the same as the flood input workspace'
 
         return issues
 
@@ -78,6 +86,32 @@ class DetectorFloodWeighting(DataProcessorAlgorithm):
         divide.setProperty("RHSWorkspace", rhs)
         divide.execute()
         return divide.getProperty("OutputWorkspace").value
+        
+    def _add(self, lhs, rhs):
+        divide = self.createChildAlgorithm("Plus")
+        divide.setProperty("LHSWorkspace", lhs)
+        divide.setProperty("RHSWorkspace", rhs)
+        divide.execute()
+        return divide.getProperty("OutputWorkspace").value
+        
+    def _integrate_bands(self, bands, in_ws):
+        # Formulate bands, integrate and sum
+        accumulated_output = None
+        for i in range(0, len(bands), 2):
+            lower = bands[i]
+            upper = bands[i+1]
+            step = upper - lower
+            rebin = self.createChildAlgorithm("Rebin")
+            rebin.setProperty("Params", [lower, step, upper])
+            rebin.setProperty("InputWorkspace", in_ws) # Always integrating the same input workspace
+            rebin.execute()
+            integrated = rebin.getProperty("OutputWorkspace").value
+            if accumulated_output:
+                accumulated_output = self._add(accumulated_output, integrated)
+            else:
+                # First band
+                accumulated_output = integrated
+        return accumulated_output
 
 
     def PyExec(self):
@@ -85,47 +119,41 @@ class DetectorFloodWeighting(DataProcessorAlgorithm):
         progress = Progress(self, 0, 1, 4) # Four coarse steps
 
         in_ws = self.getProperty('InputWorkspace').value
+        trans_ws = self.getProperty('TransmissionWorkspace').value
         bands = self.getProperty('Bands').value
 
-        # Formulate bands
-        params = list()
-        for i in range(0, len(bands), 2):
-            lower = bands[i]
-            upper = bands[i+1]
-            step = upper - lower
-            params.append((lower, step, upper))
-        progress.report()
-
-        accumulated_output = None
-        rebin = self.createChildAlgorithm("Rebin")
-        rebin.setProperty("Params", params[0])
-        rebin.setProperty("InputWorkspace", in_ws)
-        rebin.execute()
-        accumulated_output = rebin.getProperty("OutputWorkspace").value
-        progress.report()
-
-        # Determine the max across all spectra
-        y_values = accumulated_output.extractY()
-        max_val = np.amax(y_values)
-
-        # Create a workspace from the single max value
-        create = self.createChildAlgorithm("CreateSingleValuedWorkspace")
-        create.setProperty("DataValue", max_val)
-        create.execute()
-        max_ws = create.getProperty("OutputWorkspace").value
-
-        # Divide each entry by max
-        normalized = self._divide(accumulated_output, max_ws)
+        accumulated_output = self._integrate_bands(bands, in_ws)
+        if trans_ws:
+            accumulated_trans_output = self._integrate_bands(bands, trans_ws)
         progress.report()
 
         # Perform solid angle correction. Calculate solid angle then divide through.
+        normalized=accumulated_output
         if self.getProperty("SolidAngleCorrection").value:
             solidAngle = self.createChildAlgorithm("SolidAngle")
-            solidAngle.setProperty("InputWorkspace", normalized)
+            solidAngle.setProperty("InputWorkspace", accumulated_output)
             solidAngle.execute()
             solid_angle_weighting = solidAngle.getProperty("OutputWorkspace").value
             normalized = self._divide(normalized, solid_angle_weighting)
         progress.report()
+        # Divide through by the transmission workspace provided
+        if trans_ws:
+            normalized = self._divide(normalized, accumulated_trans_output)
+        # Determine the max across all spectra
+        y_values = normalized.extractY()
+        mean_val = np.mean(y_values)
+        # Create a workspace from the single max value
+        create = self.createChildAlgorithm("CreateSingleValuedWorkspace")
+        create.setProperty("DataValue", mean_val)
+        create.execute()
+        mean_ws = create.getProperty("OutputWorkspace").value
+        # Divide each entry by mean
+        normalized = self._divide(normalized, mean_ws)
+        progress.report()
+        # Fix-up ranges
+        for i in range(normalized.getNumberHistograms()):
+            normalized.dataX(i)[0] = bands[0]
+            normalized.dataX(i)[1] = bands[-1]
 
         self.setProperty('OutputWorkspace', normalized)
 

Framework/PythonInterface/test/python/plugins/algorithms/DetectorFloodWeightingTest.py

@@ -68,8 +68,53 @@ class DetectorFloodWeightingTest(unittest.TestCase):
         x_axis = out_ws.readX(0)
         self.assertEquals(x_axis[0], bands[0])
         self.assertEquals(x_axis[-1], bands[-1])
-        print out_ws.readY(0)[0]
         self.assertEquals(out_ws.readY(0)[0], 1.0)
+        
+    def test_execute_multiple_bands_no_solid_angle(self):
+        alg = AlgorithmManager.create("DetectorFloodWeighting")
+        alg.setChild(True)
+        alg.initialize()
+        alg.setProperty("SolidAngleCorrection", False)
+        signal_value = 2
+        in_ws = self._create_ws(units="Wavelength", signal_value=signal_value, data_x=range(0,10,1))
+        alg.setProperty("InputWorkspace", in_ws)
+        bands = [1,2,3,4]
+        alg.setProperty("Bands", bands) # One band
+        alg.setPropertyValue("OutputWorkspace", "dummy")
+        alg.execute()
+
+        out_ws = alg.getProperty("OutputWorkspace").value
+        self.assertEqual(1, out_ws.blocksize())
+        self.assertEqual("Wavelength", out_ws.getAxis(0).getUnit().unitID())
+        self.assertEqual(in_ws.getNumberHistograms(), out_ws.getNumberHistograms(), msg="Number of histograms should be unchanged.")
+        x_axis = out_ws.readX(0)
+        self.assertEquals(x_axis[0], bands[0])
+        self.assertEquals(x_axis[-1], bands[-1])
+        self.assertEquals(out_ws.readY(0)[0], 1.0)
+
+    def test_execute_multiple_bands_no_solid_angle_with_transmission(self):
+        alg = AlgorithmManager.create("DetectorFloodWeighting")
+        alg.setChild(True)
+        alg.initialize()
+        alg.setProperty("SolidAngleCorrection", False)
+        signal_value = 2
+        in_ws = self._create_ws(units="Wavelength", signal_value=signal_value, data_x=range(0,10,1))
+        alg.setProperty("InputWorkspace", in_ws)
+        alg.setProperty("TransmissionWorkspace", in_ws)
+        bands = [1,2,3,4]
+        alg.setProperty("Bands", bands) # One band
+        alg.setPropertyValue("OutputWorkspace", "dummy")
+        alg.execute()
+
+        out_ws = alg.getProperty("OutputWorkspace").value
+        self.assertEqual(1, out_ws.blocksize())
+        self.assertEqual("Wavelength", out_ws.getAxis(0).getUnit().unitID())
+        self.assertEqual(in_ws.getNumberHistograms(), out_ws.getNumberHistograms(), msg="Number of histograms should be unchanged.")
+        x_axis = out_ws.readX(0)
+        self.assertEquals(x_axis[0], bands[0])
+        self.assertEquals(x_axis[-1], bands[-1])
+        self.assertEquals(out_ws.readY(0)[0], 1.0)
+
 
     def test_execute_single_with_solid_angle(self):
         alg = AlgorithmManager.create("DetectorFloodWeighting")
@@ -90,5 +135,6 @@ class DetectorFloodWeightingTest(unittest.TestCase):
         self.assertEqual(in_ws.getNumberHistograms(), out_ws.getNumberHistograms(), msg="Number of histograms should be unchanged.")
 
 
+
 if __name__ == '__main__':
     unittest.main()

docs/source/algorithms/DetectorFloodWeighting-v1.rst

@@ -10,8 +10,7 @@ Description
 -----------
 This algorithm is used to calculate the detector flood weighting workspace use for pixel flood corrections. It was originally developed for the ANSTO Bilby instrument.
 
-This algorithm crops the data over the specified wavelength region, then normalizes each spectrum to the workspace spectrum maxima. The algorithm will then
-perform a solid angle correction on each spectra via :ref:`algm-SolidAngle`.
+This algorithm crops the data over the specified wavelength region, and sums it. The algorithm will then perform a solid angle correction on each spectra via :ref:`algm-SolidAngle` if specified, and divides through by the provided transmission workspace if provided. The result is divided by the mean spectrum value in the previous result.
 
 Usage
 -----
@@ -30,8 +29,6 @@ Usage
    
    print 'Number Histograms',out_ws.getNumberHistograms()
    print 'Min X:', out_ws.readX(0)[0], 'Max X:', out_ws.readX(0)[1]  
-   y_data = out_ws.extractY()
-   print  'Min Y:', np.amin(y_data), 'Max Y:', np.amax(y_data)   
 
 Output:
 
@@ -39,7 +36,6 @@ Output:
 
    Number Histograms 2
    Min X: 0.0 Max X: 10.0
-   Min Y: 0.5 Max Y: 1.0
 
 **Example - With Solid Angle Correction **