Mantid
Loading...
Searching...
No Matches
ConvertHFIRSCDtoMDE.cpp
Go to the documentation of this file.
1// Mantid Repository : https://github.com/mantidproject/mantid
2//
3// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
4// NScD Oak Ridge National Laboratory, European Spallation Source,
5// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
6// SPDX - License - Identifier: GPL - 3.0 +
11#include "MantidAPI/Run.h"
23
24#include "Eigen/Dense"
25#include "boost/math/constants/constants.hpp"
26
27namespace Mantid::MDAlgorithms {
28
31using namespace Mantid::Geometry;
32using namespace Mantid::Kernel;
33using namespace Mantid::MDAlgorithms;
34using namespace Mantid::API;
35using namespace Mantid::DataObjects;
36
37// Register the algorithm into the AlgorithmFactory
39
40//----------------------------------------------------------------------------------------------
41
42
43const std::string ConvertHFIRSCDtoMDE::name() const { return "ConvertHFIRSCDtoMDE"; }
44
46int ConvertHFIRSCDtoMDE::version() const { return 1; }
47
49const std::string ConvertHFIRSCDtoMDE::category() const { return "MDAlgorithms\\Creation"; }
50
52const std::string ConvertHFIRSCDtoMDE::summary() const {
53 return "Convert from the detector vs scan index MDHistoWorkspace into a "
54 "MDEventWorkspace with units in Q_sample.";
55}
56
57std::map<std::string, std::string> ConvertHFIRSCDtoMDE::validateInputs() {
58 std::map<std::string, std::string> result;
59
60 API::IMDHistoWorkspace_sptr inputWS = this->getProperty("InputWorkspace");
61 std::stringstream inputWSmsg;
62 if (inputWS->getNumDims() != 3) {
63 inputWSmsg << "Incorrect number of dimensions";
64 } else if (inputWS->getDimension(0)->getName() != "y" || inputWS->getDimension(1)->getName() != "x" ||
65 inputWS->getDimension(2)->getName() != "scanIndex") {
66 inputWSmsg << "Wrong dimensions";
67 } else if (inputWS->getNumExperimentInfo() == 0) {
68 inputWSmsg << "Missing experiment info";
69 } else if (inputWS->getExperimentInfo(0)->getInstrument()->getName() != "HB3A" &&
70 inputWS->getExperimentInfo(0)->getInstrument()->getName() != "WAND") {
71 inputWSmsg << "This only works for DEMAND (HB3A) or WAND (HB2C)";
72 } else if (inputWS->getDimension(2)->getNBins() != inputWS->getExperimentInfo(0)->run().getNumGoniometers()) {
73 inputWSmsg << "goniometers not set correctly, did you run SetGoniometer "
74 "with Average=False";
75 } else {
76 std::string instrument = inputWS->getExperimentInfo(0)->getInstrument()->getName();
77 const auto run = inputWS->getExperimentInfo(0)->run();
78 size_t number_of_runs = inputWS->getDimension(2)->getNBins();
79 std::vector<std::string> logs;
80 if (instrument == "HB3A")
81 logs = {"monitor", "time"};
82 else
83 logs = {"duration", "monitor_count"};
84 for (auto log : logs) {
85 if (run.hasProperty(log)) {
86 if (static_cast<size_t>(run.getLogData(log)->size()) != number_of_runs)
87 inputWSmsg << "Log " << log << " has incorrect length, ";
88 } else {
89 inputWSmsg << "Missing required log " << log << ", ";
90 }
91 }
92 }
93 if (!inputWSmsg.str().empty())
94 result["InputWorkspace"] = inputWSmsg.str();
95
96 std::vector<double> minVals = this->getProperty("MinValues");
97 std::vector<double> maxVals = this->getProperty("MaxValues");
98
99 if (minVals.size() != 3 || maxVals.size() != 3) {
100 std::stringstream msg;
101 msg << "Must provide 3 values, 1 for every dimension";
102 result["MinValues"] = msg.str();
103 result["MaxValues"] = msg.str();
104 } else {
105 std::stringstream msg;
106
107 size_t rank = minVals.size();
108 for (size_t i = 0; i < rank; ++i) {
109 if (minVals[i] >= maxVals[i]) {
110 if (msg.str().empty())
111 msg << "max not bigger than min ";
112 else
113 msg << ", ";
114 msg << "at index=" << (i + 1) << " (" << minVals[i] << ">=" << maxVals[i] << ")";
115 }
116 }
117
118 if (!msg.str().empty()) {
119 result["MinValues"] = msg.str();
120 result["MaxValues"] = msg.str();
121 }
122 }
123
124 return result;
125}
126
127//----------------------------------------------------------------------------------------------
131
132 declareProperty(std::make_unique<WorkspaceProperty<API::IMDHistoWorkspace>>("InputWorkspace", "", Direction::Input),
133 "An input workspace.");
135 std::make_unique<PropertyWithValue<double>>(
136 "Wavelength", DBL_MAX, std::make_shared<BoundedValidator<double>>(0.0, 100.0, true), Direction::Input),
137 "Wavelength");
138 declareProperty(std::make_unique<ArrayProperty<double>>("MinValues", "-10,-10,-10"),
139 "It has to be 3 comma separated values, one for each dimension in "
140 "q_sample."
141 "Values smaller then specified here will not be added to "
142 "workspace.");
143 declareProperty(std::make_unique<ArrayProperty<double>>("MaxValues", "10,10,10"),
144 "A list of the same size and the same units as MinValues "
145 "list. Values higher or equal to the specified by "
146 "this list will be ignored");
147 // Box controller properties. These are the defaults
148 this->initBoxControllerProps("5" /*SplitInto*/, 1000 /*SplitThreshold*/, 20 /*MaxRecursionDepth*/);
149 declareProperty(std::make_unique<PropertyWithValue<double>>("ObliquityParallaxCoefficient", 1.0, Direction::Input),
150 "Geometrical correction for shift in vertical beam position due to wide beam.");
151 declareProperty(std::make_unique<WorkspaceProperty<API::IMDEventWorkspace>>("OutputWorkspace", "", Direction::Output),
152 "An output workspace.");
153}
154
155//----------------------------------------------------------------------------------------------
159 double wavelength = this->getProperty("Wavelength");
160
161 API::IMDHistoWorkspace_sptr inputWS = this->getProperty("InputWorkspace");
162 auto &expInfo = *(inputWS->getExperimentInfo(static_cast<uint16_t>(0)));
163 std::string instrument = expInfo.getInstrument()->getName();
164
165 std::vector<double> twotheta, azimuthal;
166 if (instrument == "HB3A") {
167 const auto &di = expInfo.detectorInfo();
168 for (size_t x = 0; x < 512; x++) {
169 for (size_t y = 0; y < 512 * 3; y++) {
170 size_t n = x + y * 512;
171 if (!di.isMonitor(n)) {
172 twotheta.push_back(di.twoTheta(n));
173 azimuthal.push_back(di.azimuthal(n));
174 }
175 }
176 }
177 } else { // HB2C
178 azimuthal = (*(dynamic_cast<Kernel::PropertyWithValue<std::vector<double>> *>(expInfo.getLog("azimuthal"))))();
179 twotheta = (*(dynamic_cast<Kernel::PropertyWithValue<std::vector<double>> *>(expInfo.getLog("twotheta"))))();
180 }
181
182 auto outputWS = DataObjects::MDEventFactory::CreateMDWorkspace(3, "MDEvent");
184 std::vector<double> minVals = this->getProperty("MinValues");
185 std::vector<double> maxVals = this->getProperty("MaxValues");
186 outputWS->addDimension(std::make_shared<Geometry::MDHistoDimension>(
187 "Q_sample_x", "Q_sample_x", frame, static_cast<coord_t>(minVals[0]), static_cast<coord_t>(maxVals[0]), 1));
188
189 outputWS->addDimension(std::make_shared<Geometry::MDHistoDimension>(
190 "Q_sample_y", "Q_sample_y", frame, static_cast<coord_t>(minVals[1]), static_cast<coord_t>(maxVals[1]), 1));
191
192 outputWS->addDimension(std::make_shared<Geometry::MDHistoDimension>(
193 "Q_sample_z", "Q_sample_z", frame, static_cast<coord_t>(minVals[2]), static_cast<coord_t>(maxVals[2]), 1));
194 outputWS->setCoordinateSystem(Mantid::Kernel::QSample);
195 outputWS->initialize();
196
197 BoxController_sptr bc = outputWS->getBoxController();
198 this->setBoxController(bc);
199 outputWS->splitBox();
200
201 auto mdws_mdevt_3 = std::dynamic_pointer_cast<MDEventWorkspace<MDEvent<3>, 3>>(outputWS);
202 MDEventInserter<MDEventWorkspace<MDEvent<3>, 3>::sptr> inserter(mdws_mdevt_3);
203
204 double cop = this->getProperty("ObliquityParallaxCoefficient");
205 float coeff = static_cast<float>(cop);
206
207 float k = boost::math::float_constants::two_pi / static_cast<float>(wavelength);
208 // check convention to determine the sign of k
209 std::string convention = Kernel::ConfigService::Instance().getString("Q.convention");
210 if (convention == "Crystallography") {
211 k *= -1.f;
212 }
213 std::vector<Eigen::Vector3f> q_lab_pre;
214 q_lab_pre.reserve(azimuthal.size());
215 for (size_t m = 0; m < azimuthal.size(); ++m) {
216 auto twotheta_f = static_cast<float>(twotheta[m]);
217 auto azimuthal_f = static_cast<float>(azimuthal[m]);
218 q_lab_pre.push_back({-std::sin(twotheta_f) * std::cos(azimuthal_f) * k,
219 -std::sin(twotheta_f) * std::sin(azimuthal_f) * k * coeff, (1.f - std::cos(twotheta_f)) * k});
220 }
221 const auto run = inputWS->getExperimentInfo(0)->run();
222 for (size_t n = 0; n < inputWS->getDimension(2)->getNBins(); n++) {
223 auto gon = run.getGoniometerMatrix(n);
224 Eigen::Matrix3f goniometer;
225 for (int i = 0; i < 3; i++)
226 for (int j = 0; j < 3; j++)
227 goniometer(i, j) = static_cast<float>(gon[i][j]);
228 goniometer = goniometer.inverse().eval();
229 auto goniometerIndex = static_cast<uint16_t>(n);
230 for (size_t m = 0; m < azimuthal.size(); m++) {
231 size_t idx = n * azimuthal.size() + m;
232 coord_t signal = static_cast<coord_t>(inputWS->getSignalAt(idx));
233 if (signal > 0.f && std::isfinite(signal)) {
234 Eigen::Vector3f q_sample = goniometer * q_lab_pre[m];
235 inserter.insertMDEvent(signal, signal, 0, goniometerIndex, 0, q_sample.data());
236 }
237 }
238 }
239
240 auto *ts = new ThreadSchedulerFIFO();
241 ThreadPool tp(ts);
242 outputWS->splitAllIfNeeded(ts);
243 tp.joinAll();
244
245 outputWS->refreshCache();
246 outputWS->copyExperimentInfos(*inputWS);
247 auto &outRun = outputWS->getExperimentInfo(0)->mutableRun();
248 if (outRun.hasProperty("wavelength")) {
249 outRun.removeLogData("wavelength");
250 }
251 outRun.addLogData(new PropertyWithValue<double>("wavelength", wavelength));
252 outRun.getProperty("wavelength")->setUnits("Angstrom");
253
254 auto user_convention = Kernel::ConfigService::Instance().getString("Q.convention");
255 auto ws_convention = outputWS->getConvention();
256 if (user_convention != ws_convention) {
257 auto convention_alg = createChildAlgorithm("ChangeQConvention");
258 convention_alg->setProperty("InputWorkspace", outputWS);
259 convention_alg->executeAsChildAlg();
260 }
261 setProperty("OutputWorkspace", outputWS);
262}
263
264} // namespace Mantid::MDAlgorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
static std::unique_ptr< QThreadPool > tp
void declareProperty(std::unique_ptr< Kernel::Property > p, const std::string &doc="") override
Add a property to the list of managed properties.
Definition: Algorithm.cpp:1913
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
Definition: Algorithm.cpp:2076
virtual std::shared_ptr< Algorithm > createChildAlgorithm(const std::string &name, const double startProgress=-1., const double endProgress=-1., const bool enableLogging=true, const int &version=-1)
Create a Child Algorithm.
Definition: Algorithm.cpp:842
void setBoxController(const Mantid::API::BoxController_sptr &bc, const Mantid::Geometry::Instrument_const_sptr &instrument)
Set the settings in the given box controller.
void initBoxControllerProps(const std::string &SplitInto="5", int SplitThreshold=1000, int MaxRecursionDepth=5)
Initialise the properties.
A property class for workspaces.
static API::IMDEventWorkspace_sptr CreateMDWorkspace(size_t nd, const std::string &eventType="MDLeanEvent", const Mantid::API::MDNormalization &preferredNormalization=Mantid::API::MDNormalization::VolumeNormalization, const Mantid::API::MDNormalization &preferredNormalizationHisto=Mantid::API::MDNormalization::VolumeNormalization)
Create a MDEventWorkspace of the given type.
MDEventInserter : Helper class that provides a generic interface for adding events to an MDWorkspace ...
QSample : Q in the sample frame.
Definition: QSample.h:21
Support for a property that holds an array of values.
Definition: ArrayProperty.h:28
BoundedValidator is a validator that requires the values to be between upper or lower bounds,...
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
The concrete, templated class for properties.
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
A Thread Pool implementation that keeps a certain number of threads running (normally,...
Definition: ThreadPool.h:36
A First-In-First-Out Thread Scheduler.
ConvertHFIRSCDtoMDE : TODO: DESCRIPTION.
int version() const override
Algorithm's version for identification.
const std::string summary() const override
Algorithm's summary for use in the GUI and help.
void exec() override
Execute the algorithm.
const std::string category() const override
Algorithm's category for identification.
void init() override
Initialize the algorithm's properties.
std::map< std::string, std::string > validateInputs() override
Perform validation of ALL the input properties of the algorithm.
std::shared_ptr< IMDHistoWorkspace > IMDHistoWorkspace_sptr
shared pointer to Mantid::API::IMDHistoWorkspace
std::shared_ptr< BoxController > BoxController_sptr
Shared ptr to BoxController.
float coord_t
Typedef for the data type to use for coordinate axes in MD objects such as MDBox, MDEventWorkspace,...
Definition: MDTypes.h:27
STL namespace.
Describes the direction (within an algorithm) of a Property.
Definition: Property.h:50
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54