diff --git a/Framework/Algorithms/src/CalculatePlaczekSelfScattering.cpp b/Framework/Algorithms/src/CalculatePlaczekSelfScattering.cpp index 2180760c48fae408c19f66fb98c6b39269a5d02c..a27463048337d86337fe5bd9d95aaf3711b17451 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 f0d253168239af4e21c369e1a9d3c789ead2d937..94590e018be46e4d3d9fb49b999de9c66fdf7e84 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 484544b37a256b992ab5668b75db073984a93e0b..5acba2f709a9722159549480cbcce6b8e4d2a7bb 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 0cea1734feb3c1b1f86cccc633a9da007a008381..9d557356daed70143c260dc2913ac6713a180069 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