From 1993fb1a7ce796001f8e7bad8eaf006ec016cb89 Mon Sep 17 00:00:00 2001
From: Harriet Brown <harriet.brown@stfc.ac.uk>
Date: Tue, 8 Oct 2019 13:14:18 +0100
Subject: [PATCH] Fix FitIncidentSpectrum fitting wrong data by loading raw.
Also fixes FitIncidentSpectrum.py and CalculatePlaczekSelfScattering so
that output workspace has same x units
---
.../src/CalculatePlaczekSelfScattering.cpp | 5 +-
.../plugins/algorithms/FitIncidentSpectrum.py | 25 +++---
scripts/Diffraction/isis_powder/polaris.py | 3 +-
.../polaris_routines/polaris_algs.py | 89 ++++++++-----------
4 files changed, 58 insertions(+), 64 deletions(-)
diff --git a/Framework/Algorithms/src/CalculatePlaczekSelfScattering.cpp b/Framework/Algorithms/src/CalculatePlaczekSelfScattering.cpp
index 2180760c48f..a2746304833 100644
--- a/Framework/Algorithms/src/CalculatePlaczekSelfScattering.cpp
+++ b/Framework/Algorithms/src/CalculatePlaczekSelfScattering.cpp
@@ -5,6 +5,7 @@
// & Institut Laue - Langevin
// SPDX - License - Identifier: GPL - 3.0 +
#include "MantidAlgorithms/CalculatePlaczekSelfScattering.h"
+#include "MantidAPI/Axis.h"
#include "MantidAPI/Sample.h"
#include "MantidAPI/WorkspaceFactory.h"
#include "MantidDataObjects/Workspace2D.h"
@@ -12,6 +13,7 @@
#include "MantidGeometry/Instrument/DetectorInfo.h"
#include "MantidKernel/Atom.h"
#include "MantidKernel/Material.h"
+#include "MantidKernel/Unit.h"
#include <utility>
@@ -164,11 +166,12 @@ void CalculatePlaczekSelfScattering::exec() {
}
}
}
+ std::string unit = inWS->getAxis(0)->unit()->unitID();
Mantid::API::Algorithm_sptr childAlg =
createChildAlgorithm("CreateWorkspace");
childAlg->setProperty("DataX", xLambdas);
childAlg->setProperty("DataY", placzekCorrection);
- childAlg->setProperty("UnitX", "Wavelength");
+ childAlg->setProperty("UnitX", unit);
childAlg->setProperty("NSpec", nSpec);
childAlg->setProperty("ParentWorkspace", inWS);
childAlg->setProperty("Distribution", true);
diff --git a/Framework/PythonInterface/plugins/algorithms/FitIncidentSpectrum.py b/Framework/PythonInterface/plugins/algorithms/FitIncidentSpectrum.py
index f0d25316823..94590e018be 100644
--- a/Framework/PythonInterface/plugins/algorithms/FitIncidentSpectrum.py
+++ b/Framework/PythonInterface/plugins/algorithms/FitIncidentSpectrum.py
@@ -44,6 +44,11 @@ class FitIncidentSpectrum(PythonAlgorithm):
direction=Direction.Output),
doc='Output workspace containing the fit and it\'s first derivative.')
+ self.declareProperty(
+ name='SpectrumIndex',
+ defaultValue=0,
+ doc='Workspace index of the spectra to be fitted (Defaults to the first index.)')
+
self.declareProperty(FloatArrayProperty(name="BinningForCalc",
validator=RebinParamsValidator(AllowEmpty=True),
direction=Direction.Input),
@@ -67,29 +72,28 @@ class FitIncidentSpectrum(PythonAlgorithm):
def _setup(self):
self._input_ws = self.getProperty('InputWorkspace').value
self._output_ws = self.getProperty('OutputWorkspace').valueAsStr
+ self._incident_index = self.getProperty('SpectrumIndex').value
+ self._binning_for_calc = self.getProperty('BinningForCalc').value
self._binning_for_fit = self.getProperty('BinningForFit').value
self._fit_spectrum_with = self.getProperty('FitSpectrumWith').value
- self._binning_for_calc = self.getProperty('BinningForCalc').value
- if not self._binning_for_calc.all():
+ if self._binning_for_calc.size == 0:
x = self._input_ws.readX(0)
- self._binning_for_calc = [str(i) for i in [min(x), x[1] - x[0], max(x) + x[1] - x[0]]]
+ self._binning_for_calc = [i for i in [min(x), x[1] - x[0], max(x) + x[1] - x[0]]]
def PyExec(self):
self._setup()
-
x = np.arange(self._binning_for_calc[0], self._binning_for_calc[2], self._binning_for_calc[1])
- incident_index = 0
if self._binning_for_fit.size == 0:
- x_fit = np.array(self._input_ws.readX(incident_index))
- y_fit = np.array(self._input_ws.readY(incident_index))
+ x_fit = np.array(self._input_ws.readX(self._incident_index))
+ y_fit = np.array(self._input_ws.readY(self._incident_index))
else:
rebinned = Rebin(
self._input_ws,
Params=self._binning_for_fit,
PreserveEvents=True,
StoreInADS=False)
- x_fit = np.array(rebinned.readX(incident_index))
- y_fit = np.array(rebinned.readY(incident_index))
+ x_fit = np.array(rebinned.readX(self._incident_index))
+ y_fit = np.array(rebinned.readY(self._incident_index))
if len(x_fit) != len(y_fit):
x_fit = x_fit[:-1]
@@ -110,10 +114,11 @@ class FitIncidentSpectrum(PythonAlgorithm):
fit, fit_prime = self.fit_cubic_spline_with_gauss_conv(x_fit, y_fit, x, sigma=0.5)
# Create output workspace
+ unit = self._input_ws.getAxis(0).getUnit().unitID()
output_workspace = CreateWorkspace(
DataX=x,
DataY=np.append(fit, fit_prime),
- UnitX='Wavelength',
+ UnitX=unit,
NSpec=2,
Distribution=False,
ParentWorkspace=self._input_ws,
diff --git a/scripts/Diffraction/isis_powder/polaris.py b/scripts/Diffraction/isis_powder/polaris.py
index 484544b37a2..5acba2f709a 100644
--- a/scripts/Diffraction/isis_powder/polaris.py
+++ b/scripts/Diffraction/isis_powder/polaris.py
@@ -58,7 +58,8 @@ class Polaris(AbstractInst):
focus_file_path=focus_file_path,
merge_banks=self._inst_settings.merge_banks,
q_lims=q_lims,
- cal_file_name=cal_file_name)
+ cal_file_name=cal_file_name,
+ sample_details=self._sample_details)
return pdf_output
def set_sample_details(self, **kwargs):
diff --git a/scripts/Diffraction/isis_powder/polaris_routines/polaris_algs.py b/scripts/Diffraction/isis_powder/polaris_routines/polaris_algs.py
index 0cea1734feb..9d557356dae 100644
--- a/scripts/Diffraction/isis_powder/polaris_routines/polaris_algs.py
+++ b/scripts/Diffraction/isis_powder/polaris_routines/polaris_algs.py
@@ -82,13 +82,13 @@ def save_unsplined_vanadium(vanadium_ws, output_path):
mantid.SaveNexus(InputWorkspace=vanadium_ws, Filename=output_path, Append=False)
-def generate_ts_pdf(run_number, focus_file_path, merge_banks=False, q_lims=None, cal_file_name=None):
+def generate_ts_pdf(run_number, focus_file_path, merge_banks=False, q_lims=None, cal_file_name=None, sample_details=None):
focused_ws = _obtain_focused_run(run_number, focus_file_path)
- focused_ws = mantid.ConvertUnits(InputWorkspace=focused_ws.name(), Target="MomentumTransfer")
if merge_banks:
- pdf_output = _generate_grouped_ts_pdf(focused_ws, q_lims, cal_file_name)
+ pdf_output = _generate_grouped_ts_pdf(run_number, focused_ws, q_lims, cal_file_name, sample_details)
else:
+ focused_ws = mantid.ConvertUnits(InputWorkspace=focused_ws.name(), Target="MomentumTransfer")
pdf_output = mantid.PDFFourierTransform(Inputworkspace=focused_ws, InputSofQType="S(Q)", PDFType="G(r)",
Filter=True)
pdf_output = mantid.RebinToWorkspace(WorkspaceToRebin=pdf_output, WorkspaceToMatch=pdf_output[4],
@@ -123,58 +123,42 @@ def _obtain_focused_run(run_number, focus_file_path):
return focused_ws
-def _generate_grouped_ts_pdf(focused_ws, q_lims, cal_file_name):
- group_bin_min = None
- group_bin_max = None
- group_bin_width = None
+def _generate_grouped_ts_pdf(run_number, focused_ws, q_lims, cal_file_name, sample_details):
+ raw_ws = mantid.LoadRaw(Filename='POL'+str(run_number))
+ mantid.SetSample(InputWorkspace=raw_ws,
+ Geometry=common.generate_sample_geometry(sample_details),
+ Material=common.generate_sample_material(sample_details))
+ fit_spectra = mantid.FitIncidentSpectrum(InputWorkspace=raw_ws,
+ SpectrumIndex=11,
+ FitSpectrumWith="GaussConvCubicSpline")
+ placzek_self_scattering = mantid.CalculatePlaczekSelfScattering(InputWorkspace=fit_spectra)
+ cal_workspace = mantid.LoadCalFile(InputWorkspace=placzek_self_scattering,
+ CalFileName=cal_file_name,
+ Workspacename='cal_workspace',
+ MakeOffsetsWorkspace=False,
+ MakeMaskWorkspace=False)
for i in range(focused_ws.getNumberOfEntries()):
- x_array = focused_ws.getItem(i).readX(0)
- bin_min = x_array[0]
- bin_max = x_array[-1]
- bin_width = (x_array[-1] - x_array[0]) / x_array.size
- binning = [bin_min, bin_width, bin_max]
- if not group_bin_min:
- group_bin_min = bin_min
- group_bin_max = bin_max
- group_bin_width = bin_width
- else:
- group_bin_min = min(group_bin_min, bin_min)
- group_bin_max = max(group_bin_max, bin_max)
- group_bin_width = min(group_bin_width, bin_width)
- fit_spectra = mantid.FitIncidentSpectrum(InputWorkspace=focused_ws.getItem(i),
- BinningForFit=binning,
- BinningForCalc=binning,
- FitSpectrumWith="GaussConvCubicSpline")
- # find the flattest part of the fit spectra for background subtraction
- lambda_prime = fit_spectra.dataY(1)
- n_bg = int(lambda_prime.size*0.05)
- idx = np.argpartition(np.absolute(lambda_prime), n_bg)
- background = np.mean(fit_spectra.dataY(0)[idx[:n_bg]])
- placzek_self_scattering = mantid.CalculatePlaczekSelfScattering(InputWorkspace=fit_spectra)
- cal_workspace = mantid.LoadCalFile(InputWorkspace=placzek_self_scattering,
- CalFileName=cal_file_name,
- Workspacename='cal_workspace',
- MakeOffsetsWorkspace=False,
- MakeMaskWorkspace=False)
focused_correction = None
+ axis = focused_ws.getItem(i).getAxis(0)
+ unit = axis.getUnit().unitID()
for spec_index in range(placzek_self_scattering.getNumberHistograms()):
if cal_workspace.dataY(spec_index)[0] == i + 1:
if focused_correction is None:
- focused_correction = placzek_self_scattering.dataY(spec_index) + background
- else:
- focused_correction = np.add(focused_correction, placzek_self_scattering.dataY(spec_index))
+ focused_correction = placzek_self_scattering.dataY(spec_index)
+ #else:
+ #focused_correction = np.add(focused_correction, placzek_self_scattering.dataY(spec_index))
focused_correction_ws = mantid.CreateWorkspace(DataX=placzek_self_scattering.dataX(0),
DataY=focused_correction,
Distribution=True,
- UnitX='MomentumTransfer')
- mantid.Rebin(InputWorkspace=focused_ws.getItem(i), OutputWorkspace=focused_ws.getItem(i), Params=binning)
- focused_correction_ws = mantid.Rebin(InputWorkspace=focused_correction_ws, Params=binning)
+ UnitX=unit)
+ mantid.RebinToWorkspace(WorkspaceToRebin=focused_correction_ws,
+ WorkspaceToMatch=focused_ws.getItem(i),
+ OutputWorkspace=focused_correction_ws)
mantid.Subtract(LHSWorkspace=focused_ws.getItem(i),
RHSWorkspace=focused_correction_ws,
OutputWorkspace=focused_ws.getItem(i))
- binning = [group_bin_min, group_bin_width, group_bin_max]
- focused_data = mantid.Rebin(InputWorkspace=focused_ws, Params=binning)
- focused_data_combined = mantid.ConjoinSpectra(InputWorkspaces=focused_data)
+ focused_data_combined = mantid.ConjoinSpectra(InputWorkspaces=focused_ws)
+ focused_data_combined = mantid.ConvertUnits(InputWorkspace=focused_data_combined, Target="MomentumTransfer")
mantid.MatchSpectra(InputWorkspace=focused_data_combined,
OutputWorkspace=focused_data_combined,
ReferenceSpectrum=5)
@@ -195,27 +179,28 @@ def _generate_grouped_ts_pdf(focused_ws, q_lims, cal_file_name):
q_max = q_lims[1, :]
else:
raise RuntimeError("q_lims is not valid")
- pdf_x_array = focused_data_combined.readX(0)
+ bin_width = np.inf
for i in range(q_min.size):
+ pdf_x_array = focused_data_combined.readX(i)
q_min[i] = pdf_x_array[np.amin(np.where(pdf_x_array >= q_min[i]))]
q_max[i] = pdf_x_array[np.amax(np.where(pdf_x_array <= q_max[i]))]
- bin_width = pdf_x_array[1] - pdf_x_array[0]
+ bin_width = min(pdf_x_array[1] - pdf_x_array[0], bin_width)
focused_data_combined = mantid.CropWorkspaceRagged(InputWorkspace=focused_data_combined, XMin=q_min, XMax=q_max)
focused_data_combined = mantid.Rebin(InputWorkspace=focused_data_combined,
Params=[min(q_min), bin_width, max(q_max)])
focused_data_combined = mantid.SumSpectra(InputWorkspace=focused_data_combined,
WeightedSum=True,
- MultiplyBySpectra=False)
+ MultiplyBySpectra=False)
pdf_output = mantid.PDFFourierTransform(Inputworkspace=focused_data_combined,
InputSofQType="S(Q)",
PDFType="G(r)",
Filter=True)
- common.remove_intermediate_workspace(cal_workspace)
- common.remove_intermediate_workspace(fit_spectra)
- common.remove_intermediate_workspace(focused_correction_ws)
- common.remove_intermediate_workspace(focused_data)
- common.remove_intermediate_workspace(placzek_self_scattering)
+ #common.remove_intermediate_workspace(cal_workspace)
+ #common.remove_intermediate_workspace(fit_spectra)
+ #common.remove_intermediate_workspace(focused_correction_ws)
+ #common.remove_intermediate_workspace(placzek_self_scattering)
+ #common.remove_intermediate_workspace(raw_ws)
return pdf_output
--
GitLab