Skip to content
Snippets Groups Projects
MSDFit.py 7.89 KiB
Newer Older
#pylint: disable=no-init
Dan Nixon's avatar
Dan Nixon committed
from mantid.simpleapi import *
from mantid.api import *
from mantid.kernel import *


class MSDFit(DataProcessorAlgorithm):
    _output_fit_ws = None
    _spec_range = None
    _x_range = None
    _input_ws = None
    _output_param_ws = None
    _plot = None
    _output_msd_ws = None
Dan Nixon's avatar
Dan Nixon committed

    def category(self):
        return 'Workflow\\MIDAS;PythonAlgorithms'


    def summary(self):
Dan Nixon's avatar
Dan Nixon committed
        return 'Fits log(intensity) vs Q-squared to obtain the mean squared displacement.'
Dan Nixon's avatar
Dan Nixon committed


    def PyInit(self):
        self.declareProperty(MatrixWorkspaceProperty('InputWorkspace', '',direction=Direction.Input),
Dan Nixon's avatar
Dan Nixon committed
                             doc='Sample input workspace')

        self.declareProperty(name='XStart', defaultValue=0.0,
                             doc='Start of fitting range')
        self.declareProperty(name='XEnd', defaultValue=0.0,
                             doc='End of fitting range')

        self.declareProperty(name='SpecMin', defaultValue=0,
                             doc='Start of spectra range to be fit')
        self.declareProperty(name='SpecMax', defaultValue=0,
                             doc='End of spectra range to be fit')

        self.declareProperty(name='Plot', defaultValue=False,
Dan Nixon's avatar
Dan Nixon committed
                             doc='Plots results after fit')
        self.declareProperty(WorkspaceGroupProperty('OutputWorkspace', '',direction=Direction.Output),
Dan Nixon's avatar
Dan Nixon committed
                             doc='Output mean squared displacement')

        self.declareProperty(ITableWorkspaceProperty('ParameterWorkspace', '',
                                                     direction=Direction.Output,
                                                     optional=PropertyMode.Optional),
Dan Nixon's avatar
Dan Nixon committed
                             doc='Output fit parameters table')

        self.declareProperty(WorkspaceGroupProperty('FitWorkspaces', '',
                                                    direction=Direction.Output,
                                                    optional=PropertyMode.Optional),
Dan Nixon's avatar
Dan Nixon committed
                             doc='Output fitted workspaces')


    def validateInputs(self):
        issues = dict()

        workspace = self.getProperty('InputWorkspace').value
        x_data = workspace.readX(0)
Dan Nixon's avatar
Dan Nixon committed

        # Validate X axis fitting range
        x_min = self.getProperty('XStart').value
        x_max = self.getProperty('XEnd').value

        if x_min > x_max:
            msg = 'XStart must be less then XEnd'
            issues['XStart'] = msg
            issues['XEnd'] = msg

        if x_min < x_data[0]:
            issues['XStart'] = 'Must be greater than minimum X value in workspace'

        if x_max > x_data[-1]:
            issues['XEnd'] = 'Must be less than maximum X value in workspace'

Dan Nixon's avatar
Dan Nixon committed
        # Validate specta fitting range
        spec_min = self.getProperty('SpecMin').value
        spec_max = self.getProperty('SpecMax').value

        if spec_min < 0:
            issues['SpecMin'] = 'Minimum spectrum index must be greater than or equal to 0'

        if spec_max > workspace.getNumberHistograms():
            issues['SpecMax'] = 'Maximum spectrum index must be less than number of spectra in workspace'
Dan Nixon's avatar
Dan Nixon committed

        if spec_min > spec_max:
            msg = 'SpecMin must be less then SpecMax'
            issues['SpecMin'] = msg
            issues['SpecMax'] = msg

        return issues


    def PyExec(self):
        self._setup()

        # Fit line to each of the spectra
        function = 'name=LinearBackground, A0=0, A1=0'
        input_params = [self._input_ws + ',i%d' % i for i in xrange(self._spec_range[0],
                                                                    self._spec_range[1] + 1)]
Dan Nixon's avatar
Dan Nixon committed
        input_params = ';'.join(input_params)
        PlotPeakByLogValue(Input=input_params,
                           OutputWorkspace=self._output_msd_ws,
                           Function=function,
                           StartX=self._x_range[0],
                           EndX=self._x_range[1],
                           FitType='Sequential',
                           CreateOutput=True)
Dan Nixon's avatar
Dan Nixon committed

        DeleteWorkspace(self._output_msd_ws + '_NormalisedCovarianceMatrices')
        DeleteWorkspace(self._output_msd_ws + '_Parameters')
        RenameWorkspace(self._output_msd_ws,
                        OutputWorkspace=self._output_param_ws)
Dan Nixon's avatar
Dan Nixon committed

        params_table = mtd[self._output_param_ws]

        # MSD value should be positive, but the fit output is negative
        msd = params_table.column('A1')
        for i, value in enumerate(msd):
            params_table.setCell('A1', i, value * -1)

        # Create workspaces for each of the parameters
        parameter_ws_group = []

        # A0 workspace
        ws_name = self._output_msd_ws + '_A0'
        parameter_ws_group.append(ws_name)
        ConvertTableToMatrixWorkspace(self._output_param_ws, OutputWorkspace=ws_name,
                                      ColumnX='axis-1', ColumnY='A0', ColumnE='A0_Err')
        xunit = mtd[ws_name].getAxis(0).setUnit('Label')
        xunit.setLabel('Temperature', 'K')

        # A1 workspace
        ws_name = self._output_msd_ws + '_A1'
        parameter_ws_group.append(ws_name)
        ConvertTableToMatrixWorkspace(self._output_param_ws, OutputWorkspace=ws_name,
                                      ColumnX='axis-1', ColumnY='A1', ColumnE='A1_Err')
        xunit = mtd[ws_name].getAxis(0).setUnit('Label')
        xunit.setLabel('Temperature', 'K')
        SortXAxis(ws_name, OutputWorkspace=ws_name)

        # Group parameter workspaces
Dan Nixon's avatar
Dan Nixon committed
        GroupWorkspaces(InputWorkspaces=','.join(parameter_ws_group),
                        OutputWorkspace=self._output_msd_ws)
Dan Nixon's avatar
Dan Nixon committed

        # Rename fit workspace group
Dan Nixon's avatar
Dan Nixon committed
        original_fit_ws_name = self._output_msd_ws + '_Workspaces'
        if original_fit_ws_name != self._output_fit_ws:
            RenameWorkspace(InputWorkspace=self._output_msd_ws + '_Workspaces',
                            OutputWorkspace=self._output_fit_ws)
Dan Nixon's avatar
Dan Nixon committed

        # Add sample logs to output workspace
Dan Nixon's avatar
Dan Nixon committed
        CopyLogs(InputWorkspace=self._input_ws, OutputWorkspace=self._output_msd_ws)
Dan Nixon's avatar
Dan Nixon committed
        CopyLogs(InputWorkspace=self._output_msd_ws + '_A0', OutputWorkspace=self._output_fit_ws)

        self.setProperty('OutputWorkspace', self._output_msd_ws)
        self.setProperty('ParameterWorkspace', self._output_param_ws)
        self.setProperty('FitWorkspaces', self._output_fit_ws)

Dan Nixon's avatar
Dan Nixon committed
        if self._plot:
            self._plot_result()

Dan Nixon's avatar
Dan Nixon committed

    def _setup(self):
        """
        Gets algorithm properties.
        """
Dan Nixon's avatar
Dan Nixon committed
        self._input_ws = self.getPropertyValue('InputWorkspace')
Dan Nixon's avatar
Dan Nixon committed
        self._output_msd_ws = self.getPropertyValue('OutputWorkspace')
Dan Nixon's avatar
Dan Nixon committed
        self._output_param_ws = self.getPropertyValue('ParameterWorkspace')
Dan Nixon's avatar
Dan Nixon committed
        if self._output_param_ws == '':
            self._output_param_ws = self._output_msd_ws + '_Parameters'

Dan Nixon's avatar
Dan Nixon committed
        self._output_fit_ws = self.getPropertyValue('FitWorkspaces')
Dan Nixon's avatar
Dan Nixon committed
        if self._output_fit_ws == '':
            self._output_fit_ws = self._output_msd_ws + '_Workspaces'
Dan Nixon's avatar
Dan Nixon committed

        self._x_range = [self.getProperty('XStart').value,
                         self.getProperty('XEnd').value]

        self._spec_range = [self.getProperty('SpecMin').value,
                            self.getProperty('SpecMax').value]

Dan Nixon's avatar
Dan Nixon committed
        self._plot = self.getProperty('Plot').value


    def _plot_result(self):
        """
        Handles plotting result workspaces.
        """

        from IndirectImport import import_mantidplot
        mtd_plot = import_mantidplot()

        x_label = ''
        ws_run = mtd[self._input_ws].getRun()
        if 'vert_axis' in ws_run:
            x_label = ws_run.getLogData('vert_axis').value
Dan Nixon's avatar
Dan Nixon committed
        result_ws = mtd[self._output_msd_ws + '_A1']
        if len(result_ws.readX(0)) > 1:
            msd_plot = mtd_plot.plotSpectrum(result_ws, 0, True)
            msd_layer = msd_plot.activeLayer()
            msd_layer.setAxisTitle(mtd_plot.Layer.Bottom, x_label)
            msd_layer.setAxisTitle(mtd_plot.Layer.Left, '<u2>')
Dan Nixon's avatar
Dan Nixon committed


AlgorithmFactory.subscribe(MSDFit)