Skip to content
Snippets Groups Projects
Commit a6aa461a authored by Doucet, Mathieu's avatar Doucet, Mathieu
Browse files

Re #11555 remove pixel offset

parent ae40bdc0
No related merge requests found
Showing
with 173 additions and 168 deletions
#pylint: disable=no-init,invalid-name
from mantid.api import *
from mantid.simpleapi import *
from mantid.kernel import *
class LRSubtractAverageBackground(PythonAlgorithm):
def category(self):
return "Reflectometry\\SNS"
def name(self):
return "LRSubtractAverageBackground"
def version(self):
return 1
def summary(self):
return "Liquids Reflectometer background subtraction using the average on each side of the peak."
def PyInit(self):
self.declareProperty(WorkspaceProperty("InputWorkspace", "",Direction.Input), "The workspace to check.")
self.declareProperty(IntArrayProperty("PeakRange", [150, 160],
IntArrayLengthValidator(2), direction=Direction.Input),
"Pixel range defining the reflectivity peak")
self.declareProperty(IntArrayProperty("BackgroundRange", [147, 163],
IntArrayLengthValidator(2), direction=Direction.Input),
"Pixel range defining the outer range of the background on each side of the peak")
self.declareProperty(IntArrayProperty("LowResolutionRange", [94, 160],
IntArrayLengthValidator(2), direction=Direction.Input),
"Pixel range defining the low-resolution axis to integrate over")
self.declareProperty("SumPeak", False, doc="If True, the resulting peak will be summed")
self.declareProperty(WorkspaceProperty("OutputWorkspace", "",Direction.Output), "The workspace to check.")
def PyExec(self):
workspace = self.getProperty("InputWorkspace").value
# Signal region
peak_range = self.getProperty("PeakRange").value
peak_min = int(peak_range[0])
peak_max = int(peak_range[1])
# Background outer region
bck_range = self.getProperty("BackgroundRange").value
bck_min = int(bck_range[0])
bck_max = int(bck_range[1])
# Low-resolution range
x_range = self.getProperty("LowResolutionRange").value
x_min = int(x_range[0])
x_max = int(x_range[1])
sum_peak = self.getProperty("SumPeak").value
# Number of pixels in each direction
#TODO: revisit this when we update the IDF
number_of_pixels_x = int(workspace.getInstrument().getNumberParameter("number-of-x-pixels")[0])
number_of_pixels_y = int(workspace.getInstrument().getNumberParameter("number-of-y-pixels")[0])
left_bck = None
if peak_min > bck_min:
left_bck = RefRoi(InputWorkspace=workspace, IntegrateY=False,
NXPixel=number_of_pixels_x,
NYPixel=number_of_pixels_y,
ConvertToQ=False,
XPixelMin=x_min,
XPixelMax=x_max,
YPixelMin=bck_min,
YPixelMax=peak_min - 1,
ErrorWeighting = True,
SumPixels=True, NormalizeSum=True)
right_bck = None
if peak_max < bck_max:
right_bck = RefRoi(InputWorkspace=workspace, IntegrateY=False,
NXPixel=number_of_pixels_x,
NYPixel=number_of_pixels_y,
ConvertToQ=False,
XPixelMin=x_min,
XPixelMax=x_max,
YPixelMin=peak_max + 1,
YPixelMax=bck_max,
ErrorWeighting = True,
SumPixels=True, NormalizeSum=True)
if right_bck is not None and left_bck is not None:
average = (left_bck + right_bck) / 2.0
elif right_bck is not None:
average = right_bck
elif left_bck is not None:
average = left_bck
else:
average = RefRoi(InputWorkspace=workspace, IntegrateY=False,
NXPixel=number_of_pixels_x,
NYPixel=number_of_pixels_y,
ConvertToQ=False,
XPixelMin=x_min,
XPixelMax=x_max,
YPixelMin=bck_min,
YPixelMax=bck_max,
ErrorWeighting = True,
SumPixels=True, NormalizeSum=True)
# Integrate over the low-res direction
workspace = RefRoi(InputWorkspace=workspace, IntegrateY=False,
NXPixel=number_of_pixels_x,
NYPixel=number_of_pixels_y,
XPixelMin=x_min,
XPixelMax=x_max,
ConvertToQ=False,
SumPixels=sum_peak,
OutputWorkspace=str(workspace))
workspace = Minus(LHSWorkspace=workspace, RHSWorkspace=average,
OutputWorkspace=str(workspace))
# Avoid leaving trash behind
average_name = str(average)
if AnalysisDataService.doesExist(str(left_bck)):
AnalysisDataService.remove(str(left_bck))
if AnalysisDataService.doesExist(str(right_bck)):
AnalysisDataService.remove(str(right_bck))
if AnalysisDataService.doesExist(average_name):
AnalysisDataService.remove(average_name)
self.setProperty('OutputWorkspace', workspace)
AlgorithmFactory.subscribe(LRSubtractAverageBackground)
#pylint: disable=no-init,invalid-name
"""
This algorithm is a refactored version of the RefLReduction algorithm.
It was written in an attempt to:
- Not rely on external code but only on algorithms.
- Do work using existing algorithms as opposed to doing everything in arrays.
- Keep the same parameters and work as a drop-in replacement for the old algorithm.
- Reproduce the output of the old algorithm.
"""
import time
import math
import os
......@@ -21,14 +29,17 @@ class LiquidsReflectometryReduction(PythonAlgorithm):
return "Liquids Reflectometer (REFL) reduction"
def PyInit(self):
#TODO: Revisit the choice of names when we are entirely rid of the old code.
self.declareProperty(IntArrayProperty("RunNumbers"), "List of run numbers to process")
self.declareProperty("NormalizationRunNumber", 0, "Run number of the normalization run to use")
self.declareProperty(IntArrayProperty("SignalPeakPixelRange"), "Pixel range defining the data peak")
self.declareProperty(IntArrayProperty("SignalPeakPixelRange", [123, 137],
IntArrayLengthValidator(2), direction=Direction.Input),
"Pixel range defining the data peak")
self.declareProperty("SubtractSignalBackground", True,
doc='If true, the background will be subtracted from the data peak')
self.declareProperty(IntArrayProperty("SignalBackgroundPixelRange", [123, 137],
IntArrayLengthValidator(2), direction=Direction.Input),
"Pixelrange defining the background. Default:(123,137)")
"Pixel range defining the background. Default:(123,137)")
self.declareProperty("NormFlag", True, doc="If true, the data will be normalized")
self.declareProperty(IntArrayProperty("NormPeakPixelRange", [127, 133],
IntArrayLengthValidator(2), direction=Direction.Input),
......@@ -51,7 +62,7 @@ class LiquidsReflectometryReduction(PythonAlgorithm):
self.declareProperty(FloatArrayProperty("TOFRange", [0., 340000.],
FloatArrayLengthValidator(2), direction=Direction.Input),
"TOF range to use")
self.declareProperty("TofRangeFlag", True,
self.declareProperty("TOFRangeFlag", True,
doc="If true, the TOF will be cropped according to the TOF range property")
self.declareProperty("QMin", 0.05, doc="Minimum Q-value")
self.declareProperty("QStep", 0.02, doc="Step size in Q. Enter a negative value to get a log scale")
......@@ -73,11 +84,6 @@ class LiquidsReflectometryReduction(PythonAlgorithm):
"Pixel range to use for calculating the primary fraction correction.")
def PyExec(self):
# The old reflectivity reduction has an offset between the input
# pixel numbers and what it actually uses. Set the offset to zero
# to turn it off.
self.LEGACY_OFFSET = -1
# The old reduction code had a tolerance value for matching the
# slit parameters to get the scaling factors
self.TOLERANCE = 0.020
......@@ -115,7 +121,7 @@ class LiquidsReflectometryReduction(PythonAlgorithm):
SignalRange=signal_range)
# Get the TOF range
crop_TOF = self.getProperty("TofRangeFlag").value
crop_TOF = self.getProperty("TOFRangeFlag").value
tof_step = self.getProperty("TOFSteps").value
if crop_TOF:
TOFrange = self.getProperty("TOFRange").value #microS
......@@ -181,9 +187,11 @@ class LiquidsReflectometryReduction(PythonAlgorithm):
else:
normalized_data = data_cropped
# At this point, the workspace should be considered a distribution of points
normalized_data = ConvertToPointData(InputWorkspace=normalized_data,
OutputWorkspace=str(normalized_data))
normalized_data.setDistribution(True)
# Apply scaling factors
normalized_data = self.apply_scaling_factor(normalized_data)
......@@ -310,93 +318,6 @@ class LiquidsReflectometryReduction(PythonAlgorithm):
angle_offset_deg = self.getProperty("AngleOffset").value
return theta + angle_offset_deg * math.pi / 180.0
def subtract_background(self, workspace, peak_range, background_range,
low_res_range, sum_peak=False, offset=None):
"""
Subtract background in place
@param workspace: Mantid workspace
@param peak_range: range of pixels defining the peak [min, max]
@param background_range: range of pixels defining the background [min, max]
@param low_res_range: low resolution range to integrate over
@param sum_peak: if True, the resulting peak will be summed
"""
if offset is None:
offset = self.LEGACY_OFFSET
peak_min = int(peak_range[0]) + offset
peak_max = int(peak_range[1]) + offset
bck_min = int(background_range[0]) + offset
bck_max = int(background_range[1]) + offset
# Get low-resolution range
x_min = int(low_res_range[0]) + offset
x_max = int(low_res_range[1]) + offset
left_bck = None
if peak_min > bck_min:
left_bck = RefRoi(InputWorkspace=workspace, IntegrateY=False,
NXPixel=self.number_of_pixels_x,
NYPixel=self.number_of_pixels_y,
ConvertToQ=False,
XPixelMin=x_min,
XPixelMax=x_max,
YPixelMin=bck_min,
YPixelMax=peak_min - 1,
ErrorWeighting = True,
SumPixels=True, NormalizeSum=True)
right_bck = None
if peak_max < bck_max:
right_bck = RefRoi(InputWorkspace=workspace, IntegrateY=False,
NXPixel=self.number_of_pixels_x,
NYPixel=self.number_of_pixels_y,
ConvertToQ=False,
XPixelMin=x_min,
XPixelMax=x_max,
YPixelMin=peak_max + 1,
YPixelMax=bck_max,
ErrorWeighting = True,
SumPixels=True, NormalizeSum=True)
if right_bck is not None and left_bck is not None:
average = (left_bck + right_bck) / 2.0
elif right_bck is not None:
average = right_bck
elif left_bck is not None:
average = left_bck
else:
average = RefRoi(InputWorkspace=workspace, IntegrateY=False,
NXPixel=self.number_of_pixels_x,
NYPixel=self.number_of_pixels_y,
ConvertToQ=False,
XPixelMin=x_min,
XPixelMax=x_max,
YPixelMin=bck_min,
YPixelMax=bck_max,
ErrorWeighting = True,
SumPixels=True, NormalizeSum=True)
# Integrate over the low-res direction
workspace = RefRoi(InputWorkspace=workspace, IntegrateY=False,
NXPixel=self.number_of_pixels_x,
NYPixel=self.number_of_pixels_y,
XPixelMin=x_min,
XPixelMax=x_max,
ConvertToQ=False,
SumPixels=sum_peak,
OutputWorkspace=str(workspace))
#TODO Check whether we should multiply by the number of pixels
# in the low-res direction
workspace = Minus(LHSWorkspace=workspace, RHSWorkspace=average,
OutputWorkspace=str(workspace))
# Avoid leaving trash behind
average_name = str(average)
if AnalysisDataService.doesExist(str(left_bck)):
AnalysisDataService.remove(str(left_bck))
if AnalysisDataService.doesExist(str(right_bck)):
AnalysisDataService.remove(str(right_bck))
if AnalysisDataService.doesExist(average_name):
AnalysisDataService.remove(average_name)
return workspace
def process_data(self, workspace, tof_range, crop_low_res, low_res_range,
peak_range, subtract_background, background_range):
"""
......@@ -432,9 +353,11 @@ class LiquidsReflectometryReduction(PythonAlgorithm):
# Subtract background
if subtract_background:
workspace = self.subtract_background(workspace,
peak_range, background_range,
[x_min, x_max])
workspace = LRSubtractAverageBackground(InputWorkspace=workspace,
PeakRange=peak_range,
BackgroundRange=background_range,
LowResolutionRange=[x_min, x_max],
OutputWorkspace=str(workspace))
else:
# If we don't subtract the background, we still have to integrate
# over the low resolution axis
......@@ -453,8 +376,8 @@ class LiquidsReflectometryReduction(PythonAlgorithm):
# Crop to only the selected peak region
cropped = CropWorkspace(InputWorkspace = workspace,
StartWorkspaceIndex=int(peak_range[0]) + self.LEGACY_OFFSET,
EndWorkspaceIndex=int(peak_range[1]) + self.LEGACY_OFFSET,
StartWorkspaceIndex=int(peak_range[0]),
EndWorkspaceIndex=int(peak_range[1]),
OutputWorkspace="%s_cropped" % str(workspace))
# Avoid leaving trash behind
......@@ -595,8 +518,6 @@ class LiquidsReflectometryReduction(PythonAlgorithm):
norm_tof = normalization.dataX(0)
norm_value = normalization.dataY(0)
norm_error = normalization.dataE(0)
#TODO: The following is done on the bin edges.
# Should it not be done for the center of the bin?
for i in range(len(norm_value)):
norm_value[i] = norm_tof[i] * b + a
norm_error[i] = math.sqrt(a_error*a_error + norm_tof[i] * norm_tof[i] * b_error * b_error)
......
......@@ -35,13 +35,6 @@ class LiquidsReflectometryReductionTest(stresstesting.MantidStressTest):
OutputWorkspace='reflectivity_119814')
def validate(self):
# Be more tolerant with the output.
self.tolerance = 0.01
# Skip the first point so we don't have to have a big tolerance
data_y = mtd["reflectivity_119814"].dataY(0)
data_y[1] = 0.631281639115562
self.disableChecking.append('Instrument')
self.disableChecking.append('Sample')
self.disableChecking.append('SpectraMap')
......
......@@ -4,54 +4,6 @@ from mantid import *
from mantid.simpleapi import *
import sys
class LiquidsReflectometryReductionWithBackgroundTest(stresstesting.MantidStressTest):
"""
This test checks that the new liquids reflectometer reduction produces
the same results as the old code. It's more tolerant than the test below.
"""
def runTest(self):
#TODO: The reduction algorithm should not require an absolute path
scaling_factor_file = FileFinder.getFullPath("directBeamDatabaseFall2014_IPTS_11601_2.cfg")
LiquidsReflectometryReduction(RunNumbers=[119816],
NormalizationRunNumber=119692,
SignalPeakPixelRange=[155, 165],
SubtractSignalBackground=True,
SignalBackgroundPixelRange=[146, 165],
NormFlag=True,
NormPeakPixelRange=[154, 162],
NormBackgroundPixelRange=[151, 165],
SubtractNormBackground=True,
LowResDataAxisPixelRangeFlag=True,
LowResDataAxisPixelRange=[99, 158],
LowResNormAxisPixelRangeFlag=True,
LowResNormAxisPixelRange=[118, 137],
TOFRange=[9610, 22425],
IncidentMediumSelected='2InDiamSi',
GeometryCorrectionFlag=False,
QMin=0.005,
QStep=0.01,
AngleOffset=0.009,
AngleOffsetError=0.001,
ScalingFactorFile=scaling_factor_file,
SlitsWidthFlag=True,
CropFirstAndLastPoints=False,
OutputWorkspace='reflectivity_119816')
def validate(self):
# Be more tolerant with the output.
self.tolerance = 0.0002
# Skip the first point so we don't have to have a big tolerance
data_y = mtd["reflectivity_119816"].dataY(0)
data_y[1] = 0.00499601750282373
self.disableChecking.append('Instrument')
self.disableChecking.append('Sample')
self.disableChecking.append('SpectraMap')
self.disableChecking.append('Axes')
return "reflectivity_119816", 'REFL_119816.nxs'
class LiquidsReflectometryReductionWithBackgroundPreciseTest(stresstesting.MantidStressTest):
"""
......@@ -150,7 +102,7 @@ class TOFRangeOFFTest(stresstesting.MantidStressTest):
LowResNormAxisPixelRangeFlag=True,
LowResNormAxisPixelRange=[118, 137],
TOFRange=[9610, 22425],
TofRangeFlag=False,
TOFRangeFlag=False,
IncidentMediumSelected='2InDiamSi',
GeometryCorrectionFlag=False,
QMin=0.005,
......@@ -189,7 +141,7 @@ class NoBackgroundTest(stresstesting.MantidStressTest):
LowResNormAxisPixelRangeFlag=True,
LowResNormAxisPixelRange=[118, 137],
TOFRange=[9610, 22425],
TofRangeFlag=True,
TOFRangeFlag=True,
IncidentMediumSelected='2InDiamSi',
GeometryCorrectionFlag=False,
QMin=0.005,
......@@ -229,7 +181,7 @@ class TOFMismatchTest(stresstesting.MantidStressTest):
LowResNormAxisPixelRangeFlag=True,
LowResNormAxisPixelRange=[118, 137],
TOFRange=[9610, 22425],
TofRangeFlag=True,
TOFRangeFlag=True,
IncidentMediumSelected='2InDiamSi',
GeometryCorrectionFlag=False,
QMin=0.005,
......@@ -268,7 +220,7 @@ class BadDataTOFRangeTest(stresstesting.MantidStressTest):
LowResNormAxisPixelRangeFlag=True,
LowResNormAxisPixelRange=[118, 137],
TOFRange=[29623.0, 42438.0],
TofRangeFlag=True,
TOFRangeFlag=True,
IncidentMediumSelected='2InDiamSi',
GeometryCorrectionFlag=False,
QMin=0.005,
......
791a2b2fe751b8af599a5b97aa2ce3a4
32043ffb3ff3964d1b6c332b670afd7c
f7949eee903277da4a6c1fa957500546
b4185bb9f084032982784e50ededebf7
20bed05c1d84537f743afbb10dfad067
329a3f2503d9cb3bee302efcf68af181
713537859017897277e049c17d7f3fe4
88ffaa8d49e827ad3456599fd7320c7c
2460ea1adf15e55f788c0b7def73dc58
b85f69361aab0d026ed007e4b79a85a8
.. algorithm::
.. summary::
.. alias::
.. properties::
Description
-----------
Used in the Liquids Reflectometer reduction at the SNS, this algorithm
compute and subtracts the background from the reflectivity peak.
.. categories::
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment