From 09a73befa06241c487bbc26be1d89cdfc7fa48f5 Mon Sep 17 00:00:00 2001
From: Harriet Brown <harriet.brown@stfc.ac.uk>
Date: Fri, 4 Oct 2019 11:58:16 +0100
Subject: [PATCH] Create basic workflow for generate_ts_pdf with merge_banks
---
.../src/CalculatePlaczekSelfScattering.cpp | 41 +++--
scripts/Diffraction/isis_powder/polaris.py | 7 +-
.../polaris_routines/polaris_algs.py | 151 ++++++++++++++++--
.../polaris_routines/polaris_param_mapping.py | 1 +
4 files changed, 171 insertions(+), 29 deletions(-)
diff --git a/Framework/Algorithms/src/CalculatePlaczekSelfScattering.cpp b/Framework/Algorithms/src/CalculatePlaczekSelfScattering.cpp
index 30b0fee11f9..2180760c48f 100644
--- a/Framework/Algorithms/src/CalculatePlaczekSelfScattering.cpp
+++ b/Framework/Algorithms/src/CalculatePlaczekSelfScattering.cpp
@@ -32,7 +32,7 @@ getSampleSpeciesInfo(const API::MatrixWorkspace_const_sptr ws) {
for (const Kernel::Material::Material::FormulaUnit element : formula) {
const std::map<std::string, double> atomMap{
{"mass", element.atom->mass},
- {"stoich", element.multiplicity},
+ {"concentration", element.multiplicity},
{"bSqrdBar", bSqrdBar}};
atomSpecies[element.atom->symbol] = atomMap;
totalStoich += element.multiplicity;
@@ -127,7 +127,9 @@ void CalculatePlaczekSelfScattering::exec() {
const Geometry::DetectorInfo detInfo = inWS->detectorInfo();
for (size_t detIndex = 0; detIndex < detInfo.size(); detIndex++) {
if (!detInfo.isMonitor(detIndex)) {
+ if (!detInfo.l2(detIndex) == 0) {
nReserve += 1;
+ }
}
}
xLambdas.reserve(nReserve);
@@ -136,21 +138,30 @@ void CalculatePlaczekSelfScattering::exec() {
int nSpec = 0;
for (size_t detIndex = 0; detIndex < detInfo.size(); detIndex++) {
if (!detInfo.isMonitor(detIndex)) {
- const double pathLength = detInfo.l1() + detInfo.l2(detIndex);
- const double f = detInfo.l1() / pathLength;
- const double sinThetaBy2 = sin(detInfo.twoTheta(detIndex) / 2.0);
- for (size_t xIndex = 0; xIndex < xLambda.size() - 1; xIndex++) {
- const double term1 = (f - 1.0) * phi1[xIndex];
- const double term2 = f * eps1[xIndex];
- const double term3 = f - 3;
- const double inelasticPlaczekSelfCorrection =
- 2.0 * (term1 - term2 + term3) * sinThetaBy2 * sinThetaBy2 *
- summationTerm;
- xLambdas.push_back(xLambda[xIndex]);
- placzekCorrection.push_back(inelasticPlaczekSelfCorrection);
+ if (!detInfo.l2(detIndex) == 0) {
+ const double pathLength = detInfo.l1() + detInfo.l2(detIndex);
+ const double f = detInfo.l1() / pathLength;
+ const double sinThetaBy2 = sin(detInfo.twoTheta(detIndex) / 2.0);
+ for (size_t xIndex = 0; xIndex < xLambda.size() - 1; xIndex++) {
+ const double term1 = (f - 1.0) * phi1[xIndex];
+ const double term2 = f * eps1[xIndex];
+ const double term3 = f - 3;
+ const double inelasticPlaczekSelfCorrection =
+ 2.0 * (term1 - term2 + term3) * sinThetaBy2 * sinThetaBy2 *
+ summationTerm;
+ xLambdas.push_back(xLambda[xIndex]);
+ placzekCorrection.push_back(inelasticPlaczekSelfCorrection);
+ }
+ xLambdas.push_back(xLambda.back());
+ nSpec += 1;
+ } else {
+ for (size_t xIndex = 0; xIndex < xLambda.size() - 1; xIndex++) {
+ xLambdas.push_back(xLambda[xIndex]);
+ placzekCorrection.push_back(0);
+ }
+ xLambdas.push_back(xLambda.back());
+ nSpec += 1;
}
- xLambdas.push_back(xLambda.back());
- nSpec += 1;
}
}
Mantid::API::Algorithm_sptr childAlg =
diff --git a/scripts/Diffraction/isis_powder/polaris.py b/scripts/Diffraction/isis_powder/polaris.py
index cb078aaf55a..9d82c9948ee 100644
--- a/scripts/Diffraction/isis_powder/polaris.py
+++ b/scripts/Diffraction/isis_powder/polaris.py
@@ -53,13 +53,12 @@ class Polaris(AbstractInst):
# Generate pdf
run_details = self._get_run_details(self._inst_settings.run_number)
focus_file_path = self._generate_out_file_paths(run_details)["nxs_filename"]
- run_details = self._get_run_details(run_number_string=self._inst_settings.run_number)
+ cal_file_name = self._inst_settings.calibration_dir + '/' + self._inst_settings.grouping_file_name
pdf_output = polaris_algs.generate_ts_pdf(run_number=self._inst_settings.run_number,
focus_file_path=focus_file_path,
merge_banks=self._inst_settings.merge_banks,
- # q_lims=self._inst_settings.q_lims,
- # grouping_file_path=run_details.grouping_file_path
- )
+ q_lims=self._inst_settings.q_lims,
+ cal_file_name=cal_file_name)
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 fd6e6062a4f..d8335acdd19 100644
--- a/scripts/Diffraction/isis_powder/polaris_routines/polaris_algs.py
+++ b/scripts/Diffraction/isis_powder/polaris_routines/polaris_algs.py
@@ -5,7 +5,9 @@
# & Institut Laue - Langevin
# SPDX - License - Identifier: GPL - 3.0 +
from __future__ import (absolute_import, division, print_function)
+import numpy as np
+from mantid.api import WorkspaceGroup
import mantid.simpleapi as mantid
from isis_powder.routines import absorb_corrections, common
@@ -81,15 +83,58 @@ 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):
+def generate_ts_pdf(run_number, focus_file_path, merge_banks=False, q_lims=None, cal_file_name=None):
focused_ws = _obtain_focused_run(run_number, focus_file_path)
pdf_output = mantid.ConvertUnits(InputWorkspace=focused_ws.name(), Target="MomentumTransfer")
if merge_banks:
- placzek_self_scattering = mantid.CalculatePlaczekSelfScattering(InputWorkspace=pdf_output)
- pdf_output = mantid.Subtract(LHSWorkspace=pdf_output, RHSWorkspace=placzek_self_scattering)
-
- pdf_output = mantid.MatchSpectra(InputWorkspace=pdf_output, ReferenceSpectrum=1)
+ group_bin_min = None
+ group_bin_max = None
+ group_bin_width = None
+ for i in range(pdf_output.getNumberOfEntries()):
+ x_array = pdf_output.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 = min(group_bin_max, bin_max)
+ group_bin_width = min(group_bin_width, bin_width)
+ fit_spectra = mantid.FitIncidentSpectrum(InputWorkspace=pdf_output.getItem(i),
+ BinningForFit=binning,
+ BinningForCalc=binning,
+ 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)
+ focused_correction = None
+ 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)
+ else:
+ focused_correction = np.add(focused_correction, placzek_self_scattering.dataY(spec_index))
+ mantid.CreateWorkspace(DataX=placzek_self_scattering.dataX(0),
+ DataY=focused_correction,
+ Distribution=True,
+ UnitX='MomentumTransfer',
+ OutputWorkspace='focused_correction_ws')
+ mantid.Rebin(InputWorkspace=pdf_output.getItem(i), OutputWorkspace='pdf_rebined', Params=binning)
+ mantid.Rebin(InputWorkspace='focused_correction_ws', OutputWorkspace='focused_correction_ws', Params=binning)
+ mantid.Subtract(LHSWorkspace='pdf_rebined', RHSWorkspace='focused_correction_ws', OutputWorkspace=pdf_output.getItem(i))
+ binning = [group_bin_min, group_bin_width, group_bin_max]
+ pdf_output = mantid.Rebin(InputWorkspace=pdf_output, Params=binning)
+ pdf_output_combined = mantid.ConjoinSpectra(InputWorkspaces='pdf_output')
+ mantid.MatchSpectra(InputWorkspace=pdf_output_combined, OutputWorkspace=pdf_output_combined, ReferenceSpectrum=1)
if type(q_lims) == str:
q_min = []
q_max = []
@@ -102,16 +147,22 @@ def generate_ts_pdf(run_number, focus_file_path, merge_banks=False, q_lims=None)
q_max.append(value_list[3])
except IOError:
raise RuntimeError("q_lims is not valid")
- elif type(q_lims) == list:
+ elif type(q_lims) == list or type(q_lims) == np.ndarray:
q_min = q_lims[0, :]
q_max = q_lims[1, :]
else:
raise RuntimeError("q_lims is not valid")
- pdf_x_array = pdf_output.readX()
+ pdf_x_array = pdf_output_combined.readX(0)
+ for i in range(q_min.size):
+ 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]
- pdf_output = mantid.CropWorkspaceRagged(InputWorkspace=pdf_output, XMin=q_min, XMax=q_max)
- pdf_output = mantid.Rebin(InputWorkspace=pdf_output, Params=[q_min, bin_width, q_max])
- pdf_output = mantid.SumSpectra(InputWorkspace=pdf_output, WeightedSum=True)
+ pdf_output_combined = mantid.CropWorkspaceRagged(InputWorkspace=pdf_output_combined, XMin=q_min, XMax=q_max)
+ pdf_output_combined = mantid.Rebin(InputWorkspace=pdf_output_combined, Params=[min(q_min), bin_width, max(q_max)])
+ pdf_output_combined = mantid.SumSpectra(InputWorkspace=pdf_output_combined, WeightedSum=True, MultiplyBySpectra=False)
+ pdf_fourier_transform = mantid.PDFFourierTransform(Inputworkspace=pdf_output_combined, InputSofQType="S(Q)",
+ PDFType="G(r)", Filter=True)
+ return pdf_fourier_transform
pdf_output = mantid.PDFFourierTransform(Inputworkspace=pdf_output, InputSofQType="S(Q)", PDFType="G(r)",
Filter=True)
@@ -120,6 +171,86 @@ def generate_ts_pdf(run_number, focus_file_path, merge_banks=False, q_lims=None)
return pdf_output
+def _generate_grouped_ts_pdf(focused_data, q_lims, cal_file_name):
+ group_bin_min = None
+ group_bin_max = None
+ group_bin_width = None
+ for i in range(focused_data.getNumberOfEntries()):
+ x_array = focused_data.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 = min(group_bin_max, bin_max)
+ group_bin_width = min(group_bin_width, bin_width)
+ fit_spectra = mantid.FitIncidentSpectrum(InputWorkspace=focused_data.getItem(i),
+ BinningForFit=binning,
+ BinningForCalc=binning,
+ 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)
+ focused_correction = None
+ 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)
+ else:
+ focused_correction = np.add(focused_correction, placzek_self_scattering.dataY(spec_index))
+ mantid.CreateWorkspace(DataX=placzek_self_scattering.dataX(0),
+ DataY=focused_correction,
+ Distribution=True,
+ UnitX='MomentumTransfer',
+ OutputWorkspace='focused_correction_ws')
+ mantid.Rebin(InputWorkspace=focused_data.getItem(i), OutputWorkspace='focused_rebined', Params=binning)
+ mantid.Rebin(InputWorkspace='focused_correction_ws', OutputWorkspace='focused_correction_ws', Params=binning)
+ mantid.Subtract(LHSWorkspace='focused_rebined', RHSWorkspace='focused_correction_ws',
+ OutputWorkspace=focused_data.getItem(i))
+ binning = [group_bin_min, group_bin_width, group_bin_max]
+ focused_data = mantid.Rebin(InputWorkspace=focused_data, Params=binning)
+ pdf_output_combined = mantid.ConjoinSpectra(InputWorkspaces=focused_data)
+ mantid.MatchSpectra(InputWorkspace=pdf_output_combined, OutputWorkspace=pdf_output_combined, ReferenceSpectrum=1)
+ if type(q_lims) == str:
+ q_min = []
+ q_max = []
+ try:
+ with open(q_lims, 'r') as f:
+ line_list = [line.rstrip('\n') for line in f]
+ for line in line_list[:-1]:
+ value_list = line.split()
+ q_min.append(value_list[2])
+ q_max.append(value_list[3])
+ except IOError:
+ raise RuntimeError("q_lims is not valid")
+ elif type(q_lims) == list or type(q_lims) == np.ndarray:
+ q_min = q_lims[0, :]
+ q_max = q_lims[1, :]
+ else:
+ raise RuntimeError("q_lims is not valid")
+ pdf_x_array = pdf_output_combined.readX(0)
+ for i in range(q_min.size):
+ 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]
+ focused_data_combined = mantid.CropWorkspaceRagged(InputWorkspace=pdf_output_combined, XMin=q_min, XMax=q_max)
+ focused_data_combined = mantid.Rebin(InputWorkspace=pdf_output_combined, Params=[min(q_min), bin_width, max(q_max)])
+ focused_data_combined = mantid.SumSpectra(InputWorkspace=pdf_output_combined, WeightedSum=True,
+ MultiplyBySpectra=False)
+ pdf_output = mantid.PDFFourierTransform(Inputworkspace=pdf_output_combined, InputSofQType="S(Q)",
+ PDFType="G(r)", Filter=True)
+ return pdf_output
+
+
def _obtain_focused_run(run_number, focus_file_path):
"""
Searches for the focused workspace to use (based on user specified run number) in the ADS and then the output
diff --git a/scripts/Diffraction/isis_powder/polaris_routines/polaris_param_mapping.py b/scripts/Diffraction/isis_powder/polaris_routines/polaris_param_mapping.py
index 102354db49e..911791606e0 100644
--- a/scripts/Diffraction/isis_powder/polaris_routines/polaris_param_mapping.py
+++ b/scripts/Diffraction/isis_powder/polaris_routines/polaris_param_mapping.py
@@ -28,6 +28,7 @@ attr_mapping = [
ParamMapEntry(ext_name="mode", int_name="mode", enum_class=POLARIS_CHOPPER_MODES,
optional=True),
ParamMapEntry(ext_name="multiple_scattering", int_name="multiple_scattering", optional=True),
+ ParamMapEntry(ext_name="q_lims", int_name="q_lims"),
ParamMapEntry(ext_name="raw_data_cropping_values", int_name="raw_data_crop_values"),
ParamMapEntry(ext_name="run_number", int_name="run_number"),
ParamMapEntry(ext_name="sample_empty", int_name="sample_empty", optional=True),
--
GitLab