Skip to content

Commit

Permalink
Merge pull request #810 from mantidproject/10909_EnginXFitPeaks_outpu…
Browse files Browse the repository at this point in the history
…t_table_with_difc_zero

EnginXFitPeaks: add new optional property to generate table workspace
  • Loading branch information
Anders-Markvardsen committed Jun 17, 2015
2 parents 829e68c + 8161248 commit 9d4da13
Show file tree
Hide file tree
Showing 2 changed files with 134 additions and 42 deletions.
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
#pylint: disable=no-init,invalid-name
from mantid.kernel import *
from mantid.api import *
import mantid.simpleapi as sapi

import math

Expand All @@ -27,7 +28,15 @@ def PyInit(self):

self.declareProperty(FileProperty(name="ExpectedPeaksFromFile",defaultValue="",
action=FileAction.OptionalLoad,extensions = [".csv"]),
"Load from file a list of dSpacing values to be translated into TOF to find expected peaks.")
"Load from file a list of dSpacing values to be translated into TOF to "
"find expected peaks.")


self.declareProperty('OutputParametersTableName', '', direction=Direction.Input,
doc = 'Name for a table workspace with the fitted values calculated by '
'this algorithm (Difc and Zero parameters) for GSAS. '
'These two parameters are added as two columns in a single row. If not given, '
'the table workspace is not created.')

self.declareProperty("Difc", 0.0, direction = Direction.Output,\
doc = "Fitted Difc value")
Expand All @@ -54,11 +63,91 @@ def PyExec(self):

if foundPeaks.rowCount() < len(expectedPeaksTof):
txt = "Peaks effectively found: " + str(foundPeaks)[1:-1]
raise Exception("Some peaks from the list of expected peaks were not found by the algorithm "
"FindPeaks. " + txt)
raise RuntimeError("Some peaks from the list of expected peaks were not found by the algorithm "
"FindPeaks. " + txt)

difc, zero = self._fitAllPeaks(foundPeaks, expectedPeaksD)

self._produceOutputs(difc, zero)

def _readInExpectedPeaks(self):
""" Reads in expected peaks from the .csv file """
readInArray = []
exPeakArray = []
updateFileName = self.getPropertyValue("ExpectedPeaksFromFile")
if updateFileName != "":
with open(updateFileName) as f:
for line in f:
readInArray.append([float(x) for x in line.split(',')])
for a in readInArray:
for b in a:
exPeakArray.append(b)
if exPeakArray == []:
print "File could not be read. Defaults being used."
expectedPeaksD = sorted(self.getProperty('ExpectedPeaks').value)
else:
print "using file"
expectedPeaksD = sorted(exPeakArray)
else:
expectedPeaksD = sorted(self.getProperty('ExpectedPeaks').value)
return expectedPeaksD

def _getDefaultPeaks(self):
""" Gets default peaks for EnginX algorithm. Values from CeO2 """
defaultPeak = [3.1243, 2.7057, 1.9132, 1.6316, 1.5621, 1.3529, 1.2415,
1.2100, 1.1046, 1.0414, 0.9566, 0.9147, 0.9019, 0.8556,
0.8252, 0.8158, 0.7811]
return defaultPeak

def _produceOutputs(self, difc, zero):
"""
Fills in the output properties as requested. NOTE: this method and the methods that this calls
might/should be merged with similar methods in other EnginX algorithms (like EnginXCalibrate)
and possibly moved into EnginXUtils.py. That depends on how the EnginX algorithms evolve overall.
In principle, when EnginXCalibrate is updated with a similar 'OutputParametersTableName' optional
output property, it should use the 'OutputParametersTableName' of this (EnginXFitPeaks) algorithm.
@param difc :: the difc GSAS parameter as fitted here
@param zero :: the zero GSAS parameter as fitted here
"""
# mandatory outputs
self.setProperty('Difc', difc)
self.setProperty('Zero', zero)

# optional outputs
tblName = self.getPropertyValue("OutputParametersTableName")
if '' != tblName:
self._generateOutputParFitTable(tblName)

def _generateOutputParFitTable(self, name):
"""
Produces a table workspace with the two fitted parameters
@param name :: the name to use for the table workspace that is created here
"""
tbl = sapi.CreateEmptyTableWorkspace(OutputWorkspace=name)
tbl.addColumn('double', 'difc')
tbl.addColumn('double', 'zero')
tbl.addRow([float(self.getPropertyValue('difc')), float(self.getPropertyValue('zero'))])

self.log().information("Output parameters added into a table workspace: %s" % name)

def _fitAllPeaks(self, foundPeaks, expectedPeaksD):
"""
This tries to fit as many peaks as there are in the list of expected peaks passed to the algorithm.
The parameters from the (Gaussian) peaks fitted by FindPeaks elsewhere (before calling this method)
are used as initial guesses.
@param foundPeaks :: list of peaks found by FindPeaks or similar algorithm
@param expectedPeaksD :: list of expected peaks provided as input to this algorithm (in dSpacing units)
@returns difc and zero parameters obtained from fitting a
linear background (in _fitDSpacingToTOF) to the peaks fitted
here individually
"""
fittedPeaks = self._createFittedPeaksTable()

peakType = 'BackToBackExponential'
for i in range(foundPeaks.rowCount()):
row = foundPeaks.row(i)
# Peak parameters estimated by FindPeaks
Expand All @@ -74,7 +163,7 @@ def PyExec(self):
# Approximate peak intensity, assuming Gaussian shape
intensity = height * sigma * math.sqrt(2 * math.pi)

peak = FunctionFactory.createFunction("BackToBackExponential")
peak = FunctionFactory.createFunction(peakType)
peak.setParameter('X0', centre)
peak.setParameter('S', sigma)
peak.setParameter('I', intensity)
Expand Down Expand Up @@ -105,39 +194,23 @@ def PyExec(self):

fittedPeaks.addRow(fittedParams)

(difc, zero) = self._fitDSpacingToTOF(fittedPeaks)

self.setProperty('Difc', difc)
self.setProperty('Zero', zero)

def _readInExpectedPeaks(self):
""" Reads in expected peaks from the .csv file """
readInArray = []
exPeakArray = []
updateFileName = self.getPropertyValue("ExpectedPeaksFromFile")
if updateFileName != "":
with open(updateFileName) as f:
for line in f:
readInArray.append([float(x) for x in line.split(',')])
for a in readInArray:
for b in a:
exPeakArray.append(b)
if exPeakArray == []:
print "File could not be read. Defaults being used."
expectedPeaksD = sorted(self.getProperty('ExpectedPeaks').value)
else:
print "using file"
expectedPeaksD = sorted(exPeakArray)
else:
expectedPeaksD = sorted(self.getProperty('ExpectedPeaks').value)
return expectedPeaksD

def _getDefaultPeaks(self):
""" Gets default peaks for EnginX algorithm. Values from CeO2 """
defaultPeak = [3.1243, 2.7057, 1.9132, 1.6316, 1.5621, 1.3529, 1.2415,
1.2100, 1.1046, 1.0414, 0.9566, 0.9147, 0.9019, 0.8556,
0.8252, 0.8158, 0.7811]
return defaultPeak
# Check if we were able to really fit any peak
if 0 == fittedPeaks.rowCount():
detailTxt = ("Could find " + str(len(foundPeaks)) + " peaks using the algorithm FindPeaks but " +
"then it was not possible to fit any peak starting from these peaks found and using '" +
peakType + "' as peak function.")
raise RuntimeError('Could not fit any peak. Please check the list of expected peaks, as it does not '
'seem to be appropriate for the workspace given. More details: ' +
detailTxt)

# Better than failing to fit the linear function
if 1 == fittedPeaks.rowCount():
raise RuntimeError('Could find only one peak. This is not enough to fit the output parameters '
'difc and zero. Please check the list of expected peaks given and if it is '
'appropriate for the workspace')

difc, zero = self._fitDSpacingToTOF(fittedPeaks)
return (difc, zero)

def _fitDSpacingToTOF(self, fittedPeaksTable):
""" Fits a linear background to the dSpacing <-> TOF relationship and returns fitted difc
Expand Down
29 changes: 24 additions & 5 deletions Code/Mantid/docs/source/algorithms/EnginXFitPeaks-v1.rst
Original file line number Diff line number Diff line change
Expand Up @@ -15,15 +15,25 @@ Description
removed without a notification, should instrument scientists decide to do so.


The pattern is specified by providing a list of dSpacing values where Bragg peaks are expected. The algorithm then fits peaks in those areas using a peak fitting function. The dSpacing values for ExpectedPeaks are then converted by Mantid's convertUnits function to TOF.
The pattern is specified by providing a list of dSpacing values where
Bragg peaks are expected. The algorithm then fits peaks in those areas
using a peak fitting function. The dSpacing values for ExpectedPeaks
are then converted to time-of-flight (TOF) (as in the Mantid
ConvertUnits algorithm).

These values are used as start peak position in fit. It is these adjusted peak TOF value positions that are fitted against ExpectedPeaks dSpacing values according to:
These values are used as start peak position in fit. It is these
adjusted peak TOF value positions that are fitted against
ExpectedPeaks dSpacing values according to:


.. math:: TOF = DifC*d + Zero


ZERO and Difc can then be used within the GSAS program.
ZERO and Difc can then be used within the GSAS program. The
parameters DIFC and ZERO are returned and can be retrieved as output
properties as well. If a name is given in OutputParametersTableName
this algorithm also produces a table workspace with that name,
containing the two fitted (DIFC, ZERO) parameters.

Usage
-----
Expand Down Expand Up @@ -51,19 +61,28 @@ Usage
# Run the algorithm. Defaults are shown below. Files entered must be in .csv format and if both ExpectedPeaks and ExpectedPeaksFromFile are entered, the latter will be used.
# difc, zero = EnginXFitPeaks(InputWorkspace = No default, WorkspaceIndex = None, ExpectedPeaks=[0.6, 1.9], ExpectedPeaksFromFile=None)

difc, zero = EnginXFitPeaks(ws, 0, [0.65, 1.9])

out_tbl_name = 'out_params'
difc, zero = EnginXFitPeaks(ws, 0, [0.65, 1.9], OutputParametersTableName=out_tbl_name)


# Print the results
print "Difc: %.1f" % difc
print "Zero: %.1f" % zero
tbl = mtd[out_tbl_name]
print "The output table has %d row(s)" % tbl.rowCount()
print "Parameters from the table, Difc: %.1f, Zero: %.1f" % (tbl.cell(0,0), tbl.cell(0,1))

Output:

.. testcleanup:: ExTwoPeaks

DeleteWorkspace(out_tbl_name)

.. testoutput:: ExTwoPeaks

Difc: 18400.0
Zero: 46.0
The output table has 1 row(s)
Parameters from the table, Difc: 18400.0, Zero: 46.0

.. categories::

0 comments on commit 9d4da13

Please sign in to comment.