Skip to content
Snippets Groups Projects
Commit af5e34ab authored by Harriet Brown's avatar Harriet Brown
Browse files

Apply changes from code review

parent ae0988a2
No related branches found
No related tags found
No related merge requests found
......@@ -5,15 +5,13 @@
// & Institut Laue - Langevin
// SPDX - License - Identifier: GPL - 3.0 +
#include "MantidAlgorithms/CalculatePlaczekSelfScattering.h"
#include "MantidAPI/Axis.h"
#include "MantidAPI/Sample.h"
#include "MantidAPI/SpectrumInfo.h"
#include "MantidAPI/WorkspaceFactory.h"
#include "MantidDataObjects/Workspace2D.h"
#include "MantidDataObjects/WorkspaceCreation.h"
#include "MantidGeometry/Instrument/DetectorInfo.h"
#include "MantidKernel/Atom.h"
#include "MantidKernel/Material.h"
#include "MantidKernel/Unit.h"
#include <utility>
......@@ -34,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},
{"concentration", element.multiplicity},
{"stoich", element.multiplicity},
{"bSqrdBar", bSqrdBar}};
atomSpecies[element.atom->symbol] = atomMap;
totalStoich += element.multiplicity;
......@@ -51,12 +49,8 @@ void CalculatePlaczekSelfScattering::init() {
declareProperty(
std::make_unique<API::WorkspaceProperty<API::MatrixWorkspace>>(
"InputWorkspace", "", Kernel::Direction::Input),
"Raw diffraction data workspace for associated correction to be "
"calculated for. Workspace must have instument and sample data.");
declareProperty(
std::make_unique<API::WorkspaceProperty<API::MatrixWorkspace>>(
"IncidentSpecta", "", Kernel::Direction::Input),
"Workspace of fitted incident spectrum with it's first derivative.");
"Workspace of fitted incident spectrum with it's first derivative. "
"Workspace must have instument and sample data.");
declareProperty(
std::make_unique<API::WorkspaceProperty<API::MatrixWorkspace>>(
"OutputWorkspace", "", Kernel::Direction::Output),
......@@ -69,8 +63,8 @@ std::map<std::string, std::string>
CalculatePlaczekSelfScattering::validateInputs() {
std::map<std::string, std::string> issues;
const API::MatrixWorkspace_sptr inWS = getProperty("InputWorkspace");
const API::SpectrumInfo specInfo = inWS->spectrumInfo();
if (specInfo.size() == 0) {
const Geometry::DetectorInfo detInfo = inWS->detectorInfo();
if (detInfo.size() == 0) {
issues["InputWorkspace"] =
"Input workspace does not have detector information";
}
......@@ -87,7 +81,6 @@ CalculatePlaczekSelfScattering::validateInputs() {
*/
void CalculatePlaczekSelfScattering::exec() {
const API::MatrixWorkspace_sptr inWS = getProperty("InputWorkspace");
const API::MatrixWorkspace_sptr incidentWS = getProperty("IncidentSpecta");
API::MatrixWorkspace_sptr outWS = getProperty("OutputWorkspace");
constexpr double factor =
1.0 / 1.66053906660e-27; // atomic mass unit-kilogram relationship
......@@ -102,9 +95,9 @@ void CalculatePlaczekSelfScattering::exec() {
neutronMass / atom.second["mass"];
}
// get incident spectrum and 1st derivative
const MantidVec xLambda = incidentWS->readX(0);
const MantidVec incident = incidentWS->readY(0);
const MantidVec incidentPrime = incidentWS->readY(1);
const MantidVec xLambda = inWS->readX(0);
const MantidVec incident = inWS->readY(0);
const MantidVec incidentPrime = inWS->readY(1);
double dx = (xLambda[1] - xLambda[0]) / 2.0;
std::vector<double> phi1;
for (size_t i = 0; i < xLambda.size() - 1; i++) {
......@@ -131,31 +124,21 @@ void CalculatePlaczekSelfScattering::exec() {
MantidVec xLambdas;
MantidVec placzekCorrection;
size_t nReserve = 0;
const API::SpectrumInfo specInfo = inWS->spectrumInfo();
for (size_t specIndex = 0; specIndex < specInfo.size(); specIndex++) {
if (!specInfo.isMonitor(specIndex)) {
if (!specInfo.l2(specIndex) == 0) {
nReserve += 1;
}
const Geometry::DetectorInfo detInfo = inWS->detectorInfo();
for (size_t detIndex = 0; detIndex < detInfo.size(); detIndex++) {
if (!detInfo.isMonitor(detIndex)) {
nReserve += 1;
}
}
xLambdas.reserve(nReserve);
placzekCorrection.reserve(nReserve);
API::MatrixWorkspace_sptr outputWS =
DataObjects::create<API::HistoWorkspace>(*inWS);
// The algorithm computes the signal values at bin centres so they should
// be treated as a distribution
outputWS->setDistribution(true);
outputWS->setYUnit("");
outputWS->setYUnitLabel("Counts");
int nSpec = 0;
for (size_t specIndex = 0; specIndex < specInfo.size(); specIndex++) {
auto &y = outputWS->mutableY(specIndex);
auto &x = outputWS->mutableX(specIndex);
if (!specInfo.isMonitor(specIndex) && !specInfo.l2(specIndex) == 0) {
const double pathLength = specInfo.l1() + specInfo.l2(specIndex);
const double f = specInfo.l1() / pathLength;
const double sinThetaBy2 = sin(specInfo.twoTheta(specIndex) / 2.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];
......@@ -163,23 +146,24 @@ void CalculatePlaczekSelfScattering::exec() {
const double inelasticPlaczekSelfCorrection =
2.0 * (term1 - term2 + term3) * sinThetaBy2 * sinThetaBy2 *
summationTerm;
x[xIndex] = xLambda[xIndex];
y[xIndex] = inelasticPlaczekSelfCorrection;
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++) {
x[xIndex] = xLambda[xIndex];
y[xIndex] = 0;
}
x.back() = xLambda.back();
nSpec += 1;
}
}
auto incidentUnit = incidentWS->getAxis(0)->unit();
outputWS->getAxis(0)->unit() = incidentUnit;
setProperty("OutputWorkspace", outputWS);
Mantid::API::Algorithm_sptr childAlg =
createChildAlgorithm("CreateWorkspace");
childAlg->setProperty("DataX", xLambdas);
childAlg->setProperty("DataY", placzekCorrection);
childAlg->setProperty("UnitX", "Wavelength");
childAlg->setProperty("NSpec", nSpec);
childAlg->setProperty("ParentWorkspace", inWS);
childAlg->setProperty("Distribution", true);
childAlg->execute();
outWS = childAlg->getProperty("OutputWorkspace");
setProperty("OutputWorkspace", outWS);
}
} // namespace Algorithms
......
......@@ -45,7 +45,7 @@ class FitIncidentSpectrum(PythonAlgorithm):
doc='Output workspace containing the fit and it\'s first derivative.')
self.declareProperty(
name='SpectrumIndex',
name='WorkspaceIndex',
defaultValue=0,
doc='Workspace index of the spectra to be fitted (Defaults to the first index.)')
......@@ -72,7 +72,7 @@ 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._incident_index = self.getProperty('WorkspaceIndex').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
......@@ -96,12 +96,13 @@ class FitIncidentSpectrum(PythonAlgorithm):
x_fit = np.array(rebinned.readX(self._incident_index))
y_fit = np.array(rebinned.readY(self._incident_index))
x_bin_centers = [(x[i+1]+x[i])/2.0 for i in range(x.size-1)]
if len(x_fit) != len(y_fit):
x_fit = x_fit[:-1]
x_fit = [(x_fit[i+1]+x_fit[i])/2.0 for i in range(x_fit.size-1)]
if self._fit_spectrum_with == 'CubicSpline':
# Fit using cubic spline
fit, fit_prime = self.fit_cubic_spline(x_fit, y_fit, x[:-1], s=1e7)
fit, fit_prime = self.fit_cubic_spline(x_fit, y_fit, x_bin_centers, s=1e7)
elif self._fit_spectrum_with == 'CubicSplineViaMantid':
# Fit using cubic spline via Mantid
fit, fit_prime = self.fit_cubic_spline_via_mantid_spline_smoothing(
......@@ -112,7 +113,7 @@ class FitIncidentSpectrum(PythonAlgorithm):
MaxNumberOfBreaks=0)
elif self._fit_spectrum_with == 'GaussConvCubicSpline':
# Fit using Gauss conv cubic spline
fit, fit_prime = self.fit_cubic_spline_with_gauss_conv(x_fit, y_fit, x[:-1], sigma=0.5)
fit, fit_prime = self.fit_cubic_spline_with_gauss_conv(x_fit, y_fit, x_bin_centers, sigma=0.5)
# Create output workspace
unit = self._input_ws.getAxis(0).getUnit().unitID()
......
......@@ -53,7 +53,6 @@ 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"]
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,
......
......@@ -64,22 +64,19 @@ def get_run_details(run_number_string, inst_settings, is_vanadium_run):
def save_unsplined_vanadium(vanadium_ws, output_path):
converted_workspaces = []
if vanadium_ws.id() != "Workspace2D":
for ws_index in range(vanadium_ws.getNumberOfEntries()):
ws = vanadium_ws.getItem(ws_index)
previous_units = ws.getAxis(0).getUnit().unitID()
for ws_index in range(vanadium_ws.getNumberOfEntries()):
ws = vanadium_ws.getItem(ws_index)
previous_units = ws.getAxis(0).getUnit().unitID()
if previous_units != WORKSPACE_UNITS.tof:
ws = mantid.ConvertUnits(InputWorkspace=ws, Target=WORKSPACE_UNITS.tof)
if previous_units != WORKSPACE_UNITS.tof:
ws = mantid.ConvertUnits(InputWorkspace=ws, Target=WORKSPACE_UNITS.tof)
ws = mantid.RenameWorkspace(InputWorkspace=ws, OutputWorkspace="van_bank_{}".format(ws_index + 1))
converted_workspaces.append(ws)
ws = mantid.RenameWorkspace(InputWorkspace=ws, OutputWorkspace="van_bank_{}".format(ws_index + 1))
converted_workspaces.append(ws)
converted_group = mantid.GroupWorkspaces(",".join(ws.name() for ws in converted_workspaces))
mantid.SaveNexus(InputWorkspace=converted_group, Filename=output_path, Append=False)
mantid.DeleteWorkspace(converted_group)
else:
mantid.SaveNexus(InputWorkspace=vanadium_ws, Filename=output_path, Append=False)
converted_group = mantid.GroupWorkspaces(",".join(ws.name() for ws in converted_workspaces))
mantid.SaveNexus(InputWorkspace=converted_group, Filename=output_path, Append=False)
mantid.DeleteWorkspace(converted_group)
def generate_ts_pdf(run_number, focus_file_path, merge_banks=False, q_lims=None):
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment