Skip to content
Snippets Groups Projects
Commit 43071b8f authored by Peter Parker's avatar Peter Parker
Browse files

Merge remote-tracking branch 'origin/feature/9337_IN16b_reduction'

parents cbf47f05 2260935f
No related branches found
No related tags found
No related merge requests found
"""*WIKI*
A workflow algorithm to perform a data reduction for Indirect ILL instruments. Note that currently only IN16B is supported.
*WIKI*"""
from mantid.simpleapi import *
from mantid.kernel import StringListValidator, Direction
from mantid.api import DataProcessorAlgorithm, PropertyMode, AlgorithmFactory, FileProperty, FileAction, MatrixWorkspaceProperty
from mantid import config, logger, mtd
import numpy as np
import os.path
class IndirectILLReduction(DataProcessorAlgorithm):
def category(self):
return "Workflow\\MIDAS;Indirect;PythonAlgorithms"
def PyInit(self):
#input options
self.declareProperty(FileProperty('Run', '', action=FileAction.Load, extensions=["nxs"]),
doc='File path of run.')
self.declareProperty(name='Analyser', defaultValue='silicon', validator=StringListValidator(['silicon']),
doc='Analyser crystal')
self.declareProperty(name='Reflection', defaultValue='111', validator=StringListValidator(['111']),
doc='Analyser reflection')
self.declareProperty(FileProperty('MapFile', '', action=FileAction.OptionalLoad, extensions=["xml"]),
doc='Filename of the map file to use. If left blank the default will be used.')
self.declareProperty(name='MirrorMode', defaultValue=False, doc='Whether to use mirror mode')
#Output workspace properties
self.declareProperty(MatrixWorkspaceProperty("RawWorkspace", "", direction=Direction.Output),
doc="Name for the output raw workspace created.")
self.declareProperty(MatrixWorkspaceProperty("ReducedWorkspace", "", direction=Direction.Output),
doc="Name for the output reduced workspace created. If mirror mode is used this will be the sum of both"
"the left and right hand workspaces.")
self.declareProperty(MatrixWorkspaceProperty("LeftWorkspace", "", optional=PropertyMode.Optional, direction=Direction.Output),
doc="Name for the left workspace if mirror mode is used.")
self.declareProperty(MatrixWorkspaceProperty("RightWorkspace", "", optional=PropertyMode.Optional, direction=Direction.Output),
doc="Name for the right workspace if mirror mode is used.")
# output options
self.declareProperty(name='Verbose', defaultValue=False, doc='Switch Verbose Off/On')
self.declareProperty(name='Save', defaultValue=False, doc='Switch Save result to nxs file Off/On')
self.declareProperty(name='Plot', defaultValue=False, doc='Whether to plot the output workspace.')
def PyExec(self):
self.log().information('IndirectILLreduction')
run_path = self.getPropertyValue('Run')
self._raw_workspace = self.getPropertyValue('RawWorkspace')
self._red_workspace = self.getPropertyValue('ReducedWorkspace')
self._red_left_workspace = self.getPropertyValue('LeftWorkspace')
self._red_right_workspace = self.getPropertyValue('RightWorkspace')
self._map_file = self.getProperty('MapFile').value
self._use_mirror_mode = self.getProperty('MirrorMode').value
self._verbose = self.getProperty('Verbose').value
self._save = self.getProperty('Save').value
self._plot = self.getProperty('Plot').value
if self._use_mirror_mode:
if self._red_left_workspace == '':
raise ValueError("Mirror Mode requires the LeftWorkspace property to be set to a value")
if self._red_right_workspace == '':
raise ValueError("Mirror Mode requires the RightWorkspace property to be set to a value")
LoadILLIndirect(FileName=run_path, OutputWorkspace=self._raw_workspace)
instrument = mtd[self._raw_workspace].getInstrument()
self._instrument_name = instrument.getName()
self._run_number = mtd[self._raw_workspace].getRunNumber()
self._analyser = self.getPropertyValue('Analyser')
self._reflection = self.getPropertyValue('Reflection')
self._run_name = self._instrument_name + '_' + str(self._run_number)
AddSampleLog(Workspace=self._raw_workspace, LogName="mirror_sense", LogType="String", LogText=str(self._use_mirror_mode))
if self._verbose:
logger.notice('Nxs file : %s' % run_path)
output_workspaces = self._reduction()
if self._save:
workdir = config['defaultsave.directory']
for ws in output_workspaces:
file_path = os.path.join(workdir, ws + '.nxs')
SaveNexusProcessed(InputWorkspace=ws, Filename=file_path)
if self._verbose:
logger.notice('Output file : ' + file_path)
if self._plot:
from IndirectImport import import_mantidplot
mp = import_mantidplot()
graph = mp.newGraph()
for ws in output_workspaces:
mp.plotSpectrum(ws, 0, window=graph)
layer = graph.activeLayer()
layer.setAxisTitle(mp.Layer.Bottom, 'Energy Transfer (meV)')
layer.setAxisTitle(mp.Layer.Left, '')
layer.setTitle('')
self.setPropertyValue('RawWorkspace', self._raw_workspace)
self.setPropertyValue('ReducedWorkspace', self._red_workspace)
if self._use_mirror_mode:
self.setPropertyValue('LeftWorkspace', self._red_left_workspace)
self.setPropertyValue('RightWorkspace', self._red_right_workspace)
def _reduction(self):
"""
Run energy conversion for IN16B
"""
from IndirectCommon import StartTime, EndTime
StartTime('ILLindirect')
if self._verbose:
logger.notice('Input workspace : %s' % self._raw_workspace)
idf_directory = config['instrumentDefinition.directory']
ipf_name = self._instrument_name + '_' + self._analyser + '_' + self._reflection + '_Parameters.xml'
ipf_path = os.path.join(idf_directory, ipf_name)
LoadParameterFile(Workspace=self._raw_workspace, Filename=ipf_path)
AddSampleLog(Workspace=self._raw_workspace, LogName="facility", LogType="String", LogText="ILL")
if self._map_file == '':
# path name for default map file
instrument = mtd[self._raw_workspace].getInstrument()
if instrument.hasParameter('Workflow.GroupingFile'):
grouping_filename = instrument.getStringParameter('Workflow.GroupingFile')[0]
self._map_file = os.path.join(config['groupingFiles.directory'], grouping_filename)
else:
raise ValueError("Failed to find default map file. Contact development team.")
if self._verbose:
logger.notice('Map file : %s' % self._map_file)
grouped_ws = self._run_name + '_group'
GroupDetectors(InputWorkspace=self._raw_workspace, OutputWorkspace=grouped_ws, MapFile=self._map_file,
Behaviour='Average')
monitor_ws = self._run_name + '_mon'
ExtractSingleSpectrum(InputWorkspace=self._raw_workspace, OutputWorkspace=monitor_ws, WorkspaceIndex=0)
if self._use_mirror_mode:
output_workspaces = self._run_mirror_mode(monitor_ws, grouped_ws)
else:
if self._verbose:
logger.notice('Mirror sense is OFF')
self._calculate_energy(monitor_ws, grouped_ws, self._red_workspace)
output_workspaces = [self._red_workspace]
EndTime('ILLindirect')
return output_workspaces
def _run_mirror_mode(self, monitor_ws, grouped_ws):
"""
Runs energy reduction with mirror mode.
@param monitor_ws :: name of the monitor workspace
@param grouped_ws :: name of workspace with the detectors grouped
"""
if self._verbose:
logger.notice('Mirror sense is ON')
x = mtd[grouped_ws].readX(0) # energy array
mid_point = int((len(x) - 1) / 2)
#left half
left_ws = self._run_name + '_left'
left_mon_ws = left_ws + '_left_mon'
CropWorkspace(InputWorkspace=grouped_ws, OutputWorkspace=left_ws, XMax=x[mid_point - 1])
CropWorkspace(InputWorkspace=monitor_ws, OutputWorkspace=left_mon_ws, XMax=x[mid_point - 1])
self._calculate_energy(left_mon_ws, left_ws, self._red_left_workspace)
xl = mtd[self._red_left_workspace].readX(0)
if self._verbose:
logger.notice('Energy range, left : %f to %f' % (xl[0], xl[-1]))
#right half
right_ws = self._run_name + '_right'
right_mon_ws = right_ws + '_mon'
CropWorkspace(InputWorkspace=grouped_ws, OutputWorkspace=right_ws, Xmin=x[mid_point])
CropWorkspace(InputWorkspace=monitor_ws, OutputWorkspace=right_mon_ws, Xmin=x[mid_point])
x = mtd[right_ws].readX(0)
self._calculate_energy(right_mon_ws, right_ws, self._red_right_workspace)
xr = mtd[self._red_right_workspace].readX(0)
if self._verbose:
logger.notice('Energy range, right : %f to %f' % (xr[0], xr[-1]))
xl = mtd[self._red_left_workspace].readX(0)
yl = mtd[self._red_left_workspace].readY(0)
nlmax = np.argmax(np.array(yl))
xlmax = xl[nlmax]
xr = mtd[self._red_right_workspace].readX(0)
yr = mtd[self._red_right_workspace].readY(0)
nrmax = np.argmax(np.array(yr))
xrmax = xr[nrmax]
xshift = xlmax - xrmax
ScaleX(InputWorkspace=self._red_right_workspace, OutputWorkspace=self._red_right_workspace, Factor=xshift, Operation='Add')
RebinToWorkspace(WorkspaceToRebin=self._red_right_workspace, WorkspaceToMatch=self._red_left_workspace, OutputWorkspace=self._red_right_workspace)
#sum both workspaces together
Plus(LHSWorkspace=self._red_left_workspace, RHSWorkspace=self._red_right_workspace, OutputWorkspace=self._red_workspace)
Scale(InputWorkspace=self._red_workspace, OutputWorkspace=self._red_workspace, Factor=0.5, Operation='Multiply')
DeleteWorkspace(monitor_ws)
DeleteWorkspace(grouped_ws)
return [self._red_left_workspace, self._red_right_workspace, self._red_workspace]
def _calculate_energy(self, monitor_ws, grouped_ws, red_ws):
"""
Convert the input run to energy transfer
@param monitor_ws :: name of the monitor workspace to divide by
@param grouped_ws :: name of workspace with the detectors grouped
@param red_ws :: name to call the reduced workspace
"""
x_range = self._monitor_range(monitor_ws)
Scale(InputWorkspace=monitor_ws, OutputWorkspace=monitor_ws, Factor=0.001, Operation='Multiply')
CropWorkspace(InputWorkspace=monitor_ws, OutputWorkspace=monitor_ws, Xmin=x_range[0], XMax=x_range[1])
ScaleX(InputWorkspace=monitor_ws, OutputWorkspace=monitor_ws, Factor=-x_range[0], Operation='Add')
CropWorkspace(InputWorkspace=grouped_ws, OutputWorkspace=grouped_ws, Xmin=x_range[0], XMax=x_range[1])
ScaleX(InputWorkspace=grouped_ws, OutputWorkspace=grouped_ws, Factor=-x_range[0], Operation='Add')
Divide(LHSWorkspace=grouped_ws, RHSWorkspace=monitor_ws, OutputWorkspace=grouped_ws)
formula = self._energy_range(grouped_ws)
ConvertAxisByFormula(InputWorkspace=grouped_ws, OutputWorkspace=red_ws, Axis='X', Formula=formula,
AxisTitle='Energy transfer', AxisUnits='meV')
if self._verbose:
xnew = mtd[red_ws].readX(0) # energy array
logger.notice('Energy range : %f to %f' % (xnew[0], xnew[-1]))
DeleteWorkspace(grouped_ws)
DeleteWorkspace(monitor_ws)
def _monitor_range(self, monitor_ws):
"""
Get sensible values for the min and max cropping range
@param monitor_ws :: name of the monitor workspace
@return tuple containing the min and max x values in the range
"""
x = mtd[monitor_ws].readX(0) # energy array
y = mtd[monitor_ws].readY(0) # energy array
imin = np.argmax(np.array(y[0:20]))
nch = len(y)
im = np.argmax(np.array(y[nch - 21:nch - 1]))
imax = nch - 21 + im
if self._verbose:
logger.notice('Cropping range %f to %f' % (x[imin], x[imax]))
return x[imin], x[imax]
def _energy_range(self, ws):
"""
Calculate the energy range for the workspace
@param ws :: name of the workspace
@return formula for convert axis by formula to convert to energy transfer
"""
x = mtd[ws].readX(0)
npt = len(x)
imid = float(npt / 2 + 1)
gRun = mtd[ws].getRun()
wave = gRun.getLogData('wavelength').value
freq = gRun.getLogData('Doppler.doppler_frequency').value
amp = gRun.getLogData('Doppler.doppler_amplitude').value
if self._verbose:
logger.notice('Wavelength : ' + str(wave))
logger.notice('Doppler frequency : ' + str(freq))
logger.notice('Doppler amplitude : ' + str(amp))
vmax = 1.2992581918414711e-4 * freq * amp * 2.0 / wave # max energy
dele = 2.0 * vmax / npt
formula = '(x-%f)*%f' % (imid, dele)
return formula
AlgorithmFactory.subscribe(IndirectILLReduction) # Register algorithm with Mantid
......@@ -14,6 +14,7 @@ set ( TEST_PY_FILES
FilterLogByTimeTest.py
FindReflectometryLinesTest.py
GetEiT0atSNSTest.py
IndirectILLReductionTest.py
LoadFullprofFileTest.py
LoadLiveDataTest.py
LoadLogPropertyTableTest.py
......
import unittest, os
import mantid
from mantid.simpleapi import *
class IndirectILLReductionTest(unittest.TestCase):
_output_workspaces = []
def setUp(self):
run_name = 'ILLIN16B_034745'
self.kwargs = {}
self.kwargs['Run'] = run_name + '.nxs'
self.kwargs['Analyser'] = 'silicon'
self.kwargs['Reflection'] = '111'
self.kwargs['RawWorkspace'] = run_name + '_' + self.kwargs['Analyser'] + self.kwargs['Reflection'] + '_raw'
self.kwargs['ReducedWorkspace'] = run_name + '_' + self.kwargs['Analyser'] + self.kwargs['Reflection'] + '_red'
self.kwargs['Verbose'] = True
def tearDown(self):
#clean up any files we made
for ws in self._output_workspaces[1:]:
path = os.path.join(config['defaultsave.directory'], ws.name() + '.nxs')
if (os.path.isfile(path)):
try:
os.remove(path)
except IOError, e:
continue
#reset output workspaces list
self._output_workspaces = []
def test_minimal(self):
IndirectILLReduction(**self.kwargs)
red_workspace = mtd[self.kwargs['ReducedWorkspace']]
self.assertTrue(isinstance(red_workspace, mantid.api.MatrixWorkspace), "Should be a matrix workspace")
self.assertEqual("Energy transfer", red_workspace.getAxis(0).getUnit().caption())
def test_mirror_mode(self):
self.kwargs['MirrorMode'] = True
self.kwargs['LeftWorkspace'] = self.kwargs['ReducedWorkspace'] + '_left'
self.kwargs['RightWorkspace'] = self.kwargs['ReducedWorkspace'] + '_right'
IndirectILLReduction(**self.kwargs)
red_workspace = mtd[self.kwargs['ReducedWorkspace']]
left_red_workspace = mtd[self.kwargs['LeftWorkspace']]
right_red_workspace = mtd[self.kwargs['RightWorkspace']]
self.assertTrue(isinstance(red_workspace, mantid.api.MatrixWorkspace), "Should be a matrix workspace")
self.assertEqual("Energy transfer", red_workspace.getAxis(0).getUnit().caption())
self.assertEqual("meV", red_workspace.getAxis(0).getUnit().symbol().ascii())
self.assertTrue(isinstance(left_red_workspace, mantid.api.MatrixWorkspace), "Should be a matrix workspace")
self.assertEqual("Energy transfer", left_red_workspace.getAxis(0).getUnit().caption())
self.assertEqual("meV", left_red_workspace.getAxis(0).getUnit().symbol().ascii())
self.assertTrue(isinstance(right_red_workspace, mantid.api.MatrixWorkspace), "Should be a matrix workspace")
self.assertEqual("Energy transfer", right_red_workspace.getAxis(0).getUnit().caption())
self.assertEqual("meV", right_red_workspace.getAxis(0).getUnit().symbol().ascii())
def test_mirror_mode_default_names(self):
self.kwargs['MirrorMode'] = True
self.assertRaises(RuntimeError, IndirectILLReduction, **self.kwargs)
left = self.kwargs['ReducedWorkspace'] + '_left'
right = self.kwargs['ReducedWorkspace'] + '_right'
self.assertTrue(mtd.doesExist(left))
self.assertTrue(mtd.doesExist(right))
def test_save_results(self):
self.kwargs['Save'] = True
self._output_workspaces = IndirectILLReduction(**self.kwargs)
path = os.path.join(config['defaultsave.directory'], self.kwargs['ReducedWorkspace'] + '.nxs')
self.assertTrue(os.path.isfile(path))
def test_save_mirror_mode_output(self):
self.kwargs['Save'] = True
self.kwargs['MirrorMode'] = True
self.kwargs['LeftWorkspace'] = self.kwargs['ReducedWorkspace'] + '_left'
self.kwargs['RightWorkspace'] = self.kwargs['ReducedWorkspace'] + '_right'
self._output_workspaces = IndirectILLReduction(**self.kwargs)
for ws in self._output_workspaces[1:]:
path = os.path.join(config['defaultsave.directory'], ws.name() + '.nxs')
self.assertTrue(os.path.isfile(path))
def test_no_verbose(self):
self.kwargs['Verbose'] = False
IndirectILLReduction(**self.kwargs)
red_workspace = mtd[self.kwargs['ReducedWorkspace']]
self.assertTrue(isinstance(red_workspace, mantid.api.MatrixWorkspace), "Should be a matrix workspace")
def test_mapping_file_option(self):
# manually get name of grouping file from parameter file
idf = os.path.join(config['instrumentDefinition.directory'], "IN16B_Definition.xml")
ipf = os.path.join(config['instrumentDefinition.directory'], "IN16B_Parameters.xml")
ws = LoadEmptyInstrument(Filename=idf)
LoadParameterFile(ws, Filename=ipf)
instrument = ws.getInstrument()
grouping_filename = instrument.getStringParameter('Workflow.GroupingFile')[0]
DeleteWorkspace(ws)
#supply grouping file as option to algorithm, this tests that a different file can be used
#rather than reading the default directly from the IP File.
self.kwargs['MapFile'] = os.path.join(config['groupingFiles.directory'], grouping_filename)
IndirectILLReduction(**self.kwargs)
red_workspace = mtd[self.kwargs['ReducedWorkspace']]
self.assertEqual(24, red_workspace.getNumberHistograms())
if __name__ == '__main__':
unittest.main()
.. algorithm::
.. summary::
.. alias::
.. properties::
Description
-----------
A workflow algorithm to perform a data reduction for Indirect ILL instruments. Note that currently only IN16B is supported.
Usage
-----
**Example - Running IndirectILLReduction**
.. testcode:: ExIndirectILLReduction
IndirectILLReduction(Run='ILLIN16B_034745.nxs', RawWorkspace='raw_workspace', ReducedWorkspace='reduced_workspace')
print "Reduced workspace has %d spectra" % mtd['reduced_workspace'].getNumberHistograms()
print "Raw workspace has %d spectra" % mtd['raw_workspace'].getNumberHistograms()
Output:
.. testoutput:: ExIndirectILLReduction
Reduced workspace has 24 spectra
Raw workspace has 2057 spectra
**Example - Running IndirectILLReduction in mirror mode**
.. testcode:: ExIndirectILLReductionMirrorMode
IndirectILLReduction(Run='ILLIN16B_034745.nxs', RawWorkspace='raw_workspace', ReducedWorkspace='reduced_workspace', LeftWorkspace='reduced_workspace_left',
RightWorkspace='reduced_workspace_right', MirrorMode=True)
print "Raw workspace has %d spectra" % mtd['raw_workspace'].getNumberHistograms()
print "Reduced workspace has %d spectra" % mtd['reduced_workspace'].getNumberHistograms()
print "Reduced left workspace has %d spectra" % mtd['reduced_workspace_left'].getNumberHistograms()
print "Reduced right workspace has %d spectra" % mtd['reduced_workspace_right'].getNumberHistograms()
Output:
.. testoutput:: ExIndirectILLReductionMirrorMode
Raw workspace has 2057 spectra
Reduced workspace has 24 spectra
Reduced left workspace has 24 spectra
Reduced right workspace has 24 spectra
.. categories::
<?xml version="1.0" encoding="utf-8"?>
<detector-grouping instrument="IN16B">
<group name="psd1"> <ids val="2-129"/> </group>
<group name="psd2"> <ids val="130-257"/> </group>
<group name="psd3"> <ids val="258-385"/> </group>
<group name="psd4"> <ids val="386-513"/> </group>
<group name="psd5"> <ids val="514-641"/> </group>
<group name="psd6"> <ids val="642-769"/> </group>
<group name="psd7"> <ids val="770-897"/> </group>
<group name="psd8"> <ids val="898-1025"/> </group>
<group name="psd9"> <ids val="1026-1153"/> </group>
<group name="psd10"> <ids val="1154-1281"/> </group>
<group name="psd11"> <ids val="1282-1409"/> </group>
<group name="psd12"> <ids val="1410-1537"/> </group>
<group name="psd13"> <ids val="1538-1665"/> </group>
<group name="psd14"> <ids val="1666-1793"/> </group>
<group name="psd15"> <ids val="1794-1921"/> </group>
<group name="psd16"> <ids val="1922-2049"/> </group>
<group name="tube1"><detids val="2049"/></group>
<group name="tube2"><detids val="2050"/> </group>
<group name="tube3"><detids val="2051"/></group>
<group name="tube4"><detids val="2052"/> </group>
<group name="tube5"><detids val="2053"/></group>
<group name="tube6"><detids val="2054"/> </group>
<group name="tube7"><detids val="2055"/></group>
<group name="tube8"><detids val="2056"/> </group>
</detector-grouping>
......@@ -23,6 +23,9 @@
<parameter name="Workflow.beam-height" type="string">
<value val="3.0" />
</parameter>
<parameter name="Workflow.GroupingFile" type="string">
<value val="IN16B_Grouping.xml" />
</parameter>
</component-link>
</parameter-file>
<?xml version="1.0" encoding="UTF-8" ?>
<parameter-file instrument = "IN16B" valid-from = "2014-05-06T00:00:00">
<component-link name = "IN16B">
<parameter name="analysis-type" type="string">
<value val="spectroscopy" />
</parameter>
<parameter name="efixed-val">
<value val="2.082" />
</parameter>
<parameter name="analyser" type="string">
<value val="silicon" />
</parameter>
<parameter name="reflection" type="string">
<value val="111" />
</parameter>
</component-link>
<component-link name="silicon">
<parameter name="Efixed"> <value val="2.082" /> </parameter>
<parameter name="resolution"> <value val="0.0003" /> </parameter>
</component-link>
</parameter-file>
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