Commit 1e9e17c2 authored by Whitfield, Ross's avatar Whitfield, Ross
Browse files

Update HB2C and HB3A workflows using new multiple goniometers

parent 436161e3
......@@ -77,6 +77,10 @@ std::map<std::string, std::string> ConvertHFIRSCDtoMDE::validateInputs() {
inputWS->getExperimentInfo(0)->getInstrument()->getName() !=
"WAND") {
inputWSmsg << "This only works for DEMAND (HB3A) or WAND (HB2C)";
} else if (inputWS->getDimension(2)->getNBins() !=
inputWS->getExperimentInfo(0)->run().getNumGoniometers()) {
inputWSmsg << "goniometers not set correctly, did you run SetGoniometer "
"with Average=False";
} else {
std::string instrument =
inputWS->getExperimentInfo(0)->getInstrument()->getName();
......@@ -84,9 +88,9 @@ std::map<std::string, std::string> ConvertHFIRSCDtoMDE::validateInputs() {
size_t number_of_runs = inputWS->getDimension(2)->getNBins();
std::vector<std::string> logs;
if (instrument == "HB3A")
logs = {"omega", "chi", "phi", "monitor", "time"};
logs = {"monitor", "time"};
else
logs = {"duration", "monitor_count", "s1"};
logs = {"duration", "monitor_count"};
for (auto log : logs) {
if (run.hasProperty(log)) {
if (static_cast<size_t>(run.getLogData(log)->size()) != number_of_runs)
......@@ -175,17 +179,7 @@ void ConvertHFIRSCDtoMDE::exec() {
std::string instrument = expInfo.getInstrument()->getName();
std::vector<double> twotheta, azimuthal;
std::vector<double> s1, omega, chi, phi;
if (instrument == "HB3A") {
auto omegaLog = dynamic_cast<Kernel::TimeSeriesProperty<double> *>(
expInfo.run().getLogData("omega"));
omega = omegaLog->valuesAsVector();
auto chiLog = dynamic_cast<Kernel::TimeSeriesProperty<double> *>(
expInfo.run().getLogData("chi"));
chi = chiLog->valuesAsVector();
auto phiLog = dynamic_cast<Kernel::TimeSeriesProperty<double> *>(
expInfo.run().getLogData("phi"));
phi = phiLog->valuesAsVector();
const auto &di = expInfo.detectorInfo();
for (size_t x = 0; x < 512; x++) {
for (size_t y = 0; y < 512 * 3; y++) {
......@@ -197,9 +191,6 @@ void ConvertHFIRSCDtoMDE::exec() {
}
}
} else { // HB2C
auto s1Log = dynamic_cast<Kernel::TimeSeriesProperty<double> *>(
expInfo.run().getLogData("s1"));
s1 = s1Log->valuesAsVector();
azimuthal =
(*(dynamic_cast<Kernel::PropertyWithValue<std::vector<double>> *>(
expInfo.getLog("azimuthal"))))();
......@@ -245,32 +236,13 @@ void ConvertHFIRSCDtoMDE::exec() {
-sin(twotheta_f) * sin(azimuthal_f) * k,
(1.f - cos(twotheta_f)) * k});
}
const auto run = inputWS->getExperimentInfo(0)->run();
for (size_t n = 0; n < inputWS->getDimension(2)->getNBins(); n++) {
auto gon = run.getGoniometerMatrix(n);
Eigen::Matrix3f goniometer;
if (instrument == "HB3A") {
float omega_radian =
static_cast<float>(omega[n]) * boost::math::float_constants::degree;
float chi_radian =
static_cast<float>(chi[n]) * boost::math::float_constants::degree;
float phi_radian =
static_cast<float>(phi[n]) * boost::math::float_constants::degree;
Eigen::Matrix3f r1;
r1 << cos(omega_radian), 0, -sin(omega_radian), 0, 1, 0,
sin(omega_radian), 0, cos(omega_radian); // omega 0,1,0,-1
Eigen::Matrix3f r2;
r2 << cos(chi_radian), sin(chi_radian), 0, -sin(chi_radian),
cos(chi_radian), 0, 0, 0, 1; // chi 0,0,1,-1
Eigen::Matrix3f r3;
r3 << cos(phi_radian), 0, -sin(phi_radian), 0, 1, 0, sin(phi_radian), 0,
cos(phi_radian); // phi 0,1,0,-1
goniometer = r1 * r2 * r3;
} else { // HB2C
float s1_radian =
static_cast<float>(s1[n]) * boost::math::float_constants::degree;
goniometer << cos(s1_radian), 0, sin(s1_radian), 0, 1, 0, -sin(s1_radian),
0, cos(s1_radian); // s1 0,1,0,1
}
for (int i = 0; i < 3; i++)
for (int j = 0; j < 3; j++)
goniometer(i, j) = static_cast<float>(gon[i][j]);
goniometer = goniometer.inverse().eval();
for (size_t m = 0; m < azimuthal.size(); m++) {
size_t idx = n * azimuthal.size() + m;
......
......@@ -8,6 +8,7 @@
#include <cxxtest/TestSuite.h>
#include "MantidAPI/AlgorithmManager.h"
#include "MantidAPI/AnalysisDataService.h"
#include "MantidAPI/FileFinder.h"
#include "MantidAPI/IMDEventWorkspace.h"
......@@ -48,6 +49,15 @@ public:
.retrieveWS<Mantid::API::IMDHistoWorkspace>(
"ConvertHFIRSCDtoMDETest_data");
auto setGoniometer = AlgorithmManager::Instance().create("SetGoniometer");
setGoniometer->initialize();
setGoniometer->setProperty("Workspace", inputWS);
setGoniometer->setPropertyValue("Axis0", "omega,0,1,0,-1");
setGoniometer->setPropertyValue("Axis1", "chi,0,0,1,-1");
setGoniometer->setPropertyValue("Axis2", "phi,0,1,0,-1");
setGoniometer->setProperty("Average", false);
setGoniometer->execute();
ConvertHFIRSCDtoMDE alg;
// Don't put output in ADS by default
alg.setChild(true);
......
......@@ -96,8 +96,12 @@ class ConvertWANDSCDtoQ(PythonAlgorithm):
# Check that all logs are there and are of correct length
run = inWS.getExperimentInfo(0).run()
# Check number of goniometers
if run.getNumGoniometers() != number_of_runs:
issues["InputWorkspace"] = "goniometers not set correctly, did you run SetGoniometer with Average=False"
if instrument == "HB3A":
for prop in ['omega', 'chi', 'phi', 'monitor', 'time']:
for prop in ['monitor', 'time']:
if run.hasProperty(prop):
p = run.getProperty(prop).value
if np.size(p) != number_of_runs:
......@@ -105,7 +109,7 @@ class ConvertWANDSCDtoQ(PythonAlgorithm):
else:
issues["InputWorkspace"] = "missing log {}".format(prop)
else:
for prop in ['duration', 'monitor_count', 's1']:
for prop in ['duration', 'monitor_count']:
if run.hasProperty(prop):
p = run.getProperty(prop).value
if np.size(p) != number_of_runs:
......@@ -169,14 +173,6 @@ class ConvertWANDSCDtoQ(PythonAlgorithm):
progress = Progress(self, 0.0, 1.0, number_of_runs+4)
# Get rotation array
if instrument == "HB3A":
omega = np.deg2rad(inWS.getExperimentInfo(0).run().getProperty('omega').value)
chi = np.deg2rad(inWS.getExperimentInfo(0).run().getProperty('chi').value)
phi = np.deg2rad(inWS.getExperimentInfo(0).run().getProperty('phi').value)
else:
s1 = np.deg2rad(inWS.getExperimentInfo(0).run().getProperty('s1').value) + np.deg2rad(self.getProperty("S1Offset").value)
normaliseBy = self.getProperty("NormaliseBy").value
if normaliseBy == "Monitor":
if instrument == "HB3A":
......@@ -248,7 +244,7 @@ class ConvertWANDSCDtoQ(PythonAlgorithm):
if inWS.getExperimentInfo(0).getInstrument().getName() == 'HB3A':
polar = polar.reshape(512*3, 512).T.flatten()
if inWS.getExperimentInfo(0).run().hasProperty('twotheta'):
if inWS.getExperimentInfo(0).run().hasProperty('azimuthal'):
azim = np.array(inWS.getExperimentInfo(0).run().getProperty('azimuthal').value)
else:
di = inWS.getExperimentInfo(0).detectorInfo()
......@@ -287,21 +283,7 @@ class ConvertWANDSCDtoQ(PythonAlgorithm):
assert data_array[:,:,0].flags.fnc
for n in range(number_of_runs):
if instrument == "HB3A":
R1 = np.array([[np.cos(omega[n]), 0, -np.sin(omega[n])], # omega 0,1,0,-1
[ 0, 1, 0],
[np.sin(omega[n]), 0, np.cos(omega[n])]])
R2 = np.array([[ np.cos(chi[n]), np.sin(chi[n]), 0], # chi 0,0,1,-1
[-np.sin(chi[n]), np.cos(chi[n]), 0],
[ 0, 0, 1]])
R3 = np.array([[np.cos(phi[n]), 0, -np.sin(phi[n])], # phi 0,1,0,-1
[ 0, 1, 0],
[np.sin(phi[n]), 0, np.cos(phi[n])]])
R = np.dot(np.dot(R1, R2), R3)
else:
R = np.array([[ np.cos(s1[n]), 0, np.sin(s1[n])], # s1 0,1,0,1
[ 0, 1, 0],
[-np.sin(s1[n]), 0, np.cos(s1[n])]])
R = inWS.getExperimentInfo(0).run().getGoniometer(n).getR()
RUBW = np.dot(R,UBW)
q = np.round(np.dot(np.linalg.inv(RUBW),qlab.T)/bin_size-offset).astype(np.int)
q_index = np.ravel_multi_index(q, (dim0_bins+2, dim1_bins+2, dim2_bins+2), mode='clip')
......
......@@ -203,7 +203,7 @@ class HB3AAdjustSampleNorm(PythonAlgorithm):
prog.report()
self.log().information("Processing '{}'".format(in_file))
SetGoniometer(Workspace=scan, Axis0='omega,0,1,0,-1', Axis1='chi,0,0,1,-1', Axis2='phi,0,1,0,-1')
SetGoniometer(Workspace=scan, Axis0='omega,0,1,0,-1', Axis1='chi,0,0,1,-1', Axis2='phi,0,1,0,-1', Average=False)
# If processing multiple files, append the base name to the given output name
if has_multiple:
if load_files:
......
......@@ -99,7 +99,7 @@ class HB3APredictPeaks(PythonAlgorithm):
min_angle = -exp_info.run().getLogData('phi').value.max()
max_angle = -exp_info.run().getLogData('phi').value.min()
# Sometimes you get the 180 degrees off what is expected from the log
phi_log = -exp_info.run().getPropertyAsSingleValueWithTimeAveragedMean('phi')
phi_log = -exp_info.run().getLogData('phi').value[0]
if np.isclose(phi_log + 180, gon[2]):
min_angle += 180
max_angle += 180
......@@ -111,7 +111,7 @@ class HB3APredictPeaks(PythonAlgorithm):
min_angle = -exp_info.run().getLogData('omega').value.max()
max_angle = -exp_info.run().getLogData('omega').value.min()
# Sometimes you get the 180 degrees off what is expected from the log
omega_log = -exp_info.run().getPropertyAsSingleValueWithTimeAveragedMean('omega')
omega_log = -exp_info.run().getLogData('omega').value[0]
if np.isclose(omega_log + 180, gon[0]):
min_angle += 180
max_angle += 180
......
......@@ -173,6 +173,12 @@ class LoadWANDSCD(PythonAlgorithm):
outWS.getExperimentInfo(0).run().addProperty('twotheta', twotheta, True)
outWS.getExperimentInfo(0).run().addProperty('azimuthal', azimuthal, True)
setGoniometer_alg = self.createChildAlgorithm("SetGoniometer", enableLogging=False)
setGoniometer_alg.setProperty("Workspace", outWS)
setGoniometer_alg.setProperty("Axis0", 's1,0,1,0,1')
setGoniometer_alg.setProperty("Average", False)
setGoniometer_alg.execute()
self.setProperty("OutputWorkspace", outWS)
......
......@@ -4,7 +4,8 @@
# NScD Oak Ridge National Laboratory, European Spallation Source,
# Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
# SPDX - License - Identifier: GPL - 3.0 +
from mantid.simpleapi import ConvertWANDSCDtoQ, CreateMDHistoWorkspace, CreateSingleValuedWorkspace, SetUB, mtd
from mantid.simpleapi import ConvertWANDSCDtoQ, CreateMDHistoWorkspace, CreateSingleValuedWorkspace, SetUB, mtd, SetGoniometer
from mantid.kernel import FloatTimeSeriesProperty
import unittest
import numpy as np
......@@ -30,7 +31,10 @@ class ConvertWANDSCDtoQTest(unittest.TestCase):
ConvertWANDSCDtoQTest_data.addExperimentInfo(ConvertWANDSCDtoQTest_dummy)
ConvertWANDSCDtoQTest_data.getExperimentInfo(0).run().addProperty('s1', list(np.arange(0, 50, 0.5)), True)
log = FloatTimeSeriesProperty('s1')
for t, v in zip(range(100), np.arange(0, 50, 0.5)):
log.addValue(t, v)
ConvertWANDSCDtoQTest_data.getExperimentInfo(0).run()['s1'] = log
ConvertWANDSCDtoQTest_data.getExperimentInfo(0).run().addProperty('duration', [60.] * 100, True)
ConvertWANDSCDtoQTest_data.getExperimentInfo(0).run().addProperty('monitor_count', [120000.] * 100, True)
ConvertWANDSCDtoQTest_data.getExperimentInfo(0).run().addProperty('twotheta', list(
......@@ -39,6 +43,7 @@ class ConvertWANDSCDtoQTest(unittest.TestCase):
np.tile(np.linspace(-0.15, 0.15, 32), 240)), True)
SetUB(ConvertWANDSCDtoQTest_data, 5, 5, 7, 90, 90, 120, u=[-1, 0, 1], v=[1, 0, 1])
SetGoniometer(ConvertWANDSCDtoQTest_data, Axis0='s1,0,1,0,1', Average=False)
# Create Normalisation workspace
S = np.ones((32, 240, 1))
......
......@@ -4,7 +4,7 @@
# NScD Oak Ridge National Laboratory, European Spallation Source,
# Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
# SPDX - License - Identifier: GPL - 3.0 +
from mantid.simpleapi import LoadWANDSCD, SetGoniometer
from mantid.simpleapi import LoadWANDSCD
import unittest
......@@ -56,8 +56,10 @@ class LoadWANDTest(unittest.TestCase):
self.assertAlmostEqual(duration[0], 40.05, 5)
self.assertAlmostEqual(duration[1], 40.05, 5)
# test that you can run SetGoniometer, testing the s1 log
SetGoniometer(LoadWANDTest_ws, Axis0='s1,0,1,0,1')
# test that the goniometer has been set correctly
self.assertEqual(run.getNumGoniometers(), 2)
self.assertAlmostEqual(run.getGoniometer(0).getEulerAngles('YZY')[0], -142.6) # s1 from HB2C_7000
self.assertAlmostEqual(run.getGoniometer(1).getEulerAngles('YZY')[0], -142.5) # s1 from HB2C_7001
LoadWANDTest_ws.delete()
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment