Newer
Older
#include "MantidMDAlgorithms/ConvertWANDSCDtoMDE.h"
#include "MantidAPI/IMDEventWorkspace.h"
#include "MantidAPI/IMDHistoWorkspace.h"
#include "MantidAPI/IMDIterator.h"
#include "MantidAPI/Run.h"
#include "MantidDataObjects/MDBoxBase.h"
#include "MantidDataObjects/MDEventFactory.h"
#include "MantidDataObjects/MDEventInserter.h"
#include "MantidGeometry/MDGeometry/QSample.h"
#include "MantidKernel/ArrayProperty.h"
#include "MantidKernel/PropertyWithValue.h"
#include "MantidKernel/UnitLabelTypes.h"
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
namespace Mantid {
namespace MDAlgorithms {
// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
// NScD Oak Ridge National Laboratory, European Spallation Source
// & Institut Laue - Langevin
// SPDX - License - Identifier: GPL - 3.0 +
using Mantid::API::WorkspaceProperty;
using Mantid::Kernel::Direction;
using namespace Mantid::Geometry;
using namespace Mantid::Kernel;
using namespace Mantid::MDAlgorithms;
using namespace Mantid::API;
using namespace Mantid::DataObjects;
// Register the algorithm into the AlgorithmFactory
DECLARE_ALGORITHM(ConvertWANDSCDtoMDE)
//----------------------------------------------------------------------------------------------
/// Algorithms name for identification. @see Algorithm::name
const std::string ConvertWANDSCDtoMDE::name() const {
return "ConvertWANDSCDtoMDE";
}
/// Algorithm's version for identification. @see Algorithm::version
int ConvertWANDSCDtoMDE::version() const { return 1; }
/// Algorithm's category for identification. @see Algorithm::category
const std::string ConvertWANDSCDtoMDE::category() const {
return "TODO: FILL IN A CATEGORY";
}
/// Algorithm's summary for use in the GUI and help. @see Algorithm::summary
const std::string ConvertWANDSCDtoMDE::summary() const {
return "TODO: FILL IN A SUMMARY";
}
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
std::map<std::string, std::string> ConvertWANDSCDtoMDE::validateInputs() {
std::map<std::string, std::string> result;
std::vector<double> minVals = this->getProperty("MinValues");
std::vector<double> maxVals = this->getProperty("MaxValues");
if (minVals.size() != maxVals.size()) {
std::stringstream msg;
msg << "Rank of MinValues != MaxValues (" << minVals.size()
<< "!=" << maxVals.size() << ")";
result["MinValues"] = msg.str();
result["MaxValues"] = msg.str();
} else {
std::stringstream msg;
size_t rank = minVals.size();
for (size_t i = 0; i < rank; ++i) {
if (minVals[i] >= maxVals[i]) {
if (msg.str().empty())
msg << "max not bigger than min ";
else
msg << ", ";
msg << "at index=" << (i + 1) << " (" << minVals[i]
<< ">=" << maxVals[i] << ")";
}
}
if (!msg.str().empty()) {
result["MinValues"] = msg.str();
result["MaxValues"] = msg.str();
}
}
return result;
}
//----------------------------------------------------------------------------------------------
/** Initialize the algorithm's properties.
*/
void ConvertWANDSCDtoMDE::init() {
declareProperty(
std::make_unique<WorkspaceProperty<API::IMDHistoWorkspace>>(
"InputWorkspace", "", Direction::Input),
"An input workspace.");
declareProperty(
std::make_unique<WorkspaceProperty<API::IMDEventWorkspace>>(
"OutputWorkspace", "", Direction::Output),
"An output workspace.");
declareProperty(std::make_unique<PropertyWithValue<double>>(
"wavelength", Mantid::EMPTY_DBL(), Direction::Input),
"wavelength");
declareProperty(std::make_unique<ArrayProperty<double>>("MinValues"),
"It has to be N comma separated values, where N is the "
"number of dimensions of the target workspace. Values "
"smaller then specified here will not be added to "
"workspace.\n Number N is defined by properties 4,6 and 7 "
"and "
"described on *MD Transformation factory* page. See also "
":ref:`algm-ConvertToMDMinMaxLocal`");
declareProperty(std::make_unique<ArrayProperty<double>>("MaxValues"),
"A list of the same size and the same units as MinValues "
"list. Values higher or equal to the specified by "
"this list will be ignored");
}
//----------------------------------------------------------------------------------------------
/** Execute the algorithm.
*/
void ConvertWANDSCDtoMDE::exec() {
double wavelength = this->getProperty("wavelength");
if (wavelength == Mantid::EMPTY_DBL()) {
throw std::invalid_argument("wavelength not entered!");
}
API::IMDHistoWorkspace_sptr inputWS = this->getProperty("InputWorkspace");
auto &expInfo = *(inputWS->getExperimentInfo(static_cast<uint16_t>(0)));
auto s1 = (*(dynamic_cast<Kernel::PropertyWithValue<std::vector<double>> *>(
expInfo.getLog("s1"))))();
auto azimuthal =
(*(dynamic_cast<Kernel::PropertyWithValue<std::vector<double>> *>(
expInfo.getLog("azimuthal"))))();
auto twotheta =
(*(dynamic_cast<Kernel::PropertyWithValue<std::vector<double>> *>(
expInfo.getLog("twotheta"))))();
Mantid::API::IMDEventWorkspace_sptr outputWS;
outputWS = DataObjects::MDEventFactory::CreateMDWorkspace(3, "MDEvent");
Mantid::Geometry::QSample frame;
std::vector<double> minVals = this->getProperty("MinValues");
std::vector<double> maxVals = this->getProperty("MaxValues");
outputWS->addDimension(boost::make_shared<Geometry::MDHistoDimension>(
"Q_sample_x", "Q_sample_x", frame, static_cast<coord_t>(minVals[0]),
static_cast<coord_t>(maxVals[0]), 1));
outputWS->addDimension(boost::make_shared<Geometry::MDHistoDimension>(
"Q_sample_y", "Q_sample_y", frame, static_cast<coord_t>(minVals[1]),
static_cast<coord_t>(maxVals[1]), 1));
outputWS->addDimension(boost::make_shared<Geometry::MDHistoDimension>(
"Q_sample_z", "Q_sample_z", frame, static_cast<coord_t>(minVals[2]),
static_cast<coord_t>(maxVals[2]), 1));
outputWS->setCoordinateSystem(Mantid::Kernel::QSample);
outputWS->initialize();
BoxController_sptr bc = outputWS->getBoxController();
bc->setSplitThreshold(1000);
bc->setMaxDepth(20);
bc->setSplitInto(5);
outputWS->splitBox();
MDEventWorkspace<MDEvent<3>, 3>::sptr mdws_mdevt_3 =
boost::dynamic_pointer_cast<MDEventWorkspace<MDEvent<3>, 3>>(outputWS);
MDEventInserter<MDEventWorkspace<MDEvent<3>, 3>::sptr> inserter(mdws_mdevt_3);
double k = 2. * M_PI / wavelength;
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
for (size_t n = 0; n < s1.size(); n++) {
Matrix<double> goniometer(3, 3, true);
goniometer[0][0] = cos(s1[n] * M_PI / 180);
goniometer[0][2] = sin(s1[n] * M_PI / 180);
goniometer[2][0] = -sin(s1[n] * M_PI / 180);
goniometer[2][2] = cos(s1[n] * M_PI / 180);
goniometer.Invert();
for (size_t m = 0; m < azimuthal.size(); m++) {
std::vector<double> q_lab(3);
std::vector<Mantid::coord_t> q_sample(3);
q_lab[0] = -sin(twotheta[m]) * cos(azimuthal[m]) * k;
q_lab[1] = -sin(twotheta[m]) * sin(azimuthal[m]) * k;
q_lab[2] = (1 - cos(twotheta[m])) * k;
q_lab = goniometer * q_lab;
q_sample[0] = static_cast<Mantid::coord_t>(q_lab[0]);
q_sample[1] = static_cast<Mantid::coord_t>(q_lab[1]);
q_sample[2] = static_cast<Mantid::coord_t>(q_lab[2]);
size_t idx = n * azimuthal.size() + m;
signal_t signal = inputWS->getSignalAt(idx);
if (signal > 0)
inserter.insertMDEvent(static_cast<float>(signal),
static_cast<float>(signal), 0, 0,
q_sample.data());
}
}
auto *ts = new ThreadSchedulerFIFO();
ThreadPool tp(ts);
outputWS->splitAllIfNeeded(ts);
tp.joinAll();
outputWS->refreshCache();
outputWS->copyExperimentInfos(*inputWS);
setProperty("OutputWorkspace", outputWS);
}
} // namespace MDAlgorithms
} // namespace Mantid