Mantid
Loading...
Searching...
No Matches
ConvertMDHistoToMatrixWorkspace.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 +
7//----------------------------------------------------------------------
8// Includes
9//----------------------------------------------------------------------
17#include "MantidHistogramData/LinearGenerator.h"
18
20#include "MantidKernel/Unit.h"
22
23#include <sstream>
24
25using namespace Mantid::Kernel;
26using namespace Mantid::API;
27
28namespace {
32struct NullDeleter {
33 void operator()(void const * /*unused*/) const { // Do nothing
34 }
35};
36
52size_t findXAxis(const VMD &start, const VMD &end, CoordTransform const *const transform,
53 IMDHistoWorkspace const *const inputWorkspace, Logger &logger, const size_t id,
54 std::string &xAxisLabel) {
55
56 // Find the start and end points in the original workspace
57 VMD originalStart = transform->applyVMD(start);
58 VMD originalEnd = transform->applyVMD(end);
59 VMD diff = originalEnd - originalStart;
60
61 // Now we find the dimension with the biggest change
62 double largest = -1e30;
63
64 size_t dimIndex = id;
65
66 const size_t nOriginalWorkspaces = inputWorkspace->numOriginalWorkspaces();
67 if (nOriginalWorkspaces < 1) {
68 logger.information("No original workspaces. Assume X-axis is Dim0.");
69 return dimIndex;
70 }
71
72 auto originalWS =
73 std::dynamic_pointer_cast<IMDWorkspace>(inputWorkspace->getOriginalWorkspace(nOriginalWorkspaces - 1));
74
75 for (size_t d = 0; d < diff.getNumDims(); d++) {
76 if (fabs(diff[d]) > largest || (originalWS->getDimension(dimIndex)->getIsIntegrated())) {
77 // Skip over any integrated dimensions
78 if (originalWS && !originalWS->getDimension(d)->getIsIntegrated()) {
79 largest = fabs(diff[d]);
80 dimIndex = d;
81 }
82 }
83 }
84 // Use the x-axis label from the original workspace.
85 xAxisLabel = originalWS->getDimension(dimIndex)->getName();
86 return dimIndex;
87}
88} // namespace
89
90namespace Mantid::MDAlgorithms {
91
92// Register the algorithm into the AlgorithmFactory
93DECLARE_ALGORITHM(ConvertMDHistoToMatrixWorkspace)
94
95
97 declareProperty(std::make_unique<WorkspaceProperty<API::IMDHistoWorkspace>>("InputWorkspace", "", Direction::Input),
98 "An input IMDHistoWorkspace.");
99 declareProperty(std::make_unique<WorkspaceProperty<>>("OutputWorkspace", "", Direction::Output),
100 "An output Workspace2D.");
101
102 std::array<std::string, 3> normalizations = {{"NoNormalization", "VolumeNormalization", "NumEventsNormalization"}};
103
104 declareProperty("Normalization", normalizations[0],
106 "Signal normalization method");
107
108 declareProperty(std::make_unique<PropertyWithValue<bool>>("FindXAxis", true, Direction::Input),
109 "If True, tries to automatically determine the dimension to use as the "
110 "output x-axis. Applies to line cut MD workspaces.");
111}
112
115
116 IMDHistoWorkspace_sptr inputWorkspace = getProperty("InputWorkspace");
117 Mantid::Geometry::VecIMDDimension_const_sptr nonIntegDims = inputWorkspace->getNonIntegratedDimensions();
118
119 if (nonIntegDims.size() == 1) {
121 } else if (nonIntegDims.size() == 2) {
123 } else {
124 throw std::invalid_argument("Cannot convert MD workspace with more than 2 dimensions.");
125 }
126}
127
132 IMDHistoWorkspace_sptr inputWorkspace = getProperty("InputWorkspace");
133
134 Mantid::Geometry::VecIMDDimension_const_sptr nonIntegDims = inputWorkspace->getNonIntegratedDimensions();
135
136 std::string alongDim;
137 if (!nonIntegDims.empty())
138 alongDim = nonIntegDims[0]->getDimensionId();
139 else
140 alongDim = inputWorkspace->getDimension(0)->getDimensionId();
141
142 size_t nd = inputWorkspace->getNumDims();
143 Mantid::Kernel::VMD start = VMD(nd);
144 Mantid::Kernel::VMD end = VMD(nd);
145
146 size_t id = 0;
147 for (size_t d = 0; d < nd; d++) {
148 Mantid::Geometry::IMDDimension_const_sptr dim = inputWorkspace->getDimension(d);
149 if (dim->getDimensionId() == alongDim) {
150 // All the way through in the single dimension
151 start[d] = dim->getMinimum();
152 end[d] = dim->getMaximum();
153 id = d; // We take the first non integrated dimension to be the diemnsion
154 // of interest.
155 } else {
156 // Mid point along each dimension
157 start[d] = (dim->getMaximum() + dim->getMinimum()) / 2.0f;
158 end[d] = start[d];
159 }
160 }
161
162 // Unit direction of the line
163 Mantid::Kernel::VMD dir = end - start;
164 dir.normalize();
165
166 std::string normProp = getPropertyValue("Normalization");
167
168 Mantid::API::MDNormalization normalization;
169 if (normProp == "NoNormalization") {
170 normalization = NoNormalization;
171 } else if (normProp == "VolumeNormalization") {
172 normalization = VolumeNormalization;
173 } else if (normProp == "NumEventsNormalization") {
174 normalization = NumEventsNormalization;
175 } else {
176 normalization = NoNormalization;
177 }
178
179 auto line = inputWorkspace->getLineData(start, end, normalization);
180
181 MatrixWorkspace_sptr outputWorkspace =
182 WorkspaceFactory::Instance().create("Workspace2D", 1, line.x.size(), line.y.size());
183 outputWorkspace->mutableY(0) = line.y;
184 outputWorkspace->mutableE(0) = line.e;
185
186 const size_t numberTransformsToOriginal = inputWorkspace->getNumberTransformsToOriginal();
187
188 CoordTransform_const_sptr transform = std::make_shared<NullCoordTransform>(inputWorkspace->getNumDims());
189 if (numberTransformsToOriginal > 0) {
190 const size_t indexToLastTransform = numberTransformsToOriginal - 1;
191 transform = CoordTransform_const_sptr(inputWorkspace->getTransformToOriginal(indexToLastTransform), NullDeleter());
192 }
193
194 assert(line.x.size() == outputWorkspace->x(0).size());
195
196 std::string xAxisLabel = inputWorkspace->getDimension(id)->getName();
197 const bool autoFind = this->getProperty("FindXAxis");
198 if (autoFind) {
199 // We look to the original workspace if possbible to find the dimension of
200 // interest to plot against.
201 id = findXAxis(start, end, transform.get(), inputWorkspace.get(), g_log, id, xAxisLabel);
202 }
203
204 auto &mutableXValues = outputWorkspace->mutableX(0);
205 // VMD inTargetCoord;
206 for (size_t i = 0; i < line.x.size(); ++i) {
207 // Coordinates in the workspace being plotted
208 VMD wsCoord = start + dir * line.x[i];
209
210 VMD inTargetCoord = transform->applyVMD(wsCoord);
211 mutableXValues[i] = inTargetCoord[id];
212 }
213 // outputWorkspace->mutableX(0) = inTargetCoord;
214
215 std::shared_ptr<Kernel::Units::Label> labelX =
216 std::dynamic_pointer_cast<Kernel::Units::Label>(Kernel::UnitFactory::Instance().create("Label"));
217 labelX->setLabel(xAxisLabel);
218 outputWorkspace->getAxis(0)->unit() = labelX;
219
220 outputWorkspace->setYUnitLabel("Signal");
221
222 setProperty("OutputWorkspace", outputWorkspace);
223}
224
229 // get the input workspace
230 IMDHistoWorkspace_sptr inputWorkspace = getProperty("InputWorkspace");
231
232 // find the non-integrated dimensions
233 Mantid::Geometry::VecIMDDimension_const_sptr nonIntegDims = inputWorkspace->getNonIntegratedDimensions();
234
235 auto xDim = nonIntegDims[0];
236 auto yDim = nonIntegDims[1];
237
238 size_t nx = xDim->getNBins();
239 size_t ny = yDim->getNBins();
240
241 size_t xDimIndex = inputWorkspace->getDimensionIndexById(xDim->getDimensionId());
242 size_t xStride = calcStride(*inputWorkspace, xDimIndex);
243
244 size_t yDimIndex = inputWorkspace->getDimensionIndexById(yDim->getDimensionId());
245 size_t yStride = calcStride(*inputWorkspace, yDimIndex);
246
247 // get the normalization of the output
248 std::string normProp = getPropertyValue("Normalization");
249 Mantid::API::MDNormalization normalization;
250 if (normProp == "NoNormalization") {
251 normalization = NoNormalization;
252 } else if (normProp == "VolumeNormalization") {
253 normalization = VolumeNormalization;
254 } else if (normProp == "NumEventsNormalization") {
255 normalization = NumEventsNormalization;
256 } else {
257 normalization = NoNormalization;
258 }
259 auto inverseVolume = static_cast<signal_t>(inputWorkspace->getInverseVolume());
260
261 // create the output workspace
262 MatrixWorkspace_sptr outputWorkspace = WorkspaceFactory::Instance().create("Workspace2D", ny, nx + 1, nx);
263
264 // set the x-values
265 const size_t xValsSize = outputWorkspace->x(0).size();
266 const double dx = xDim->getBinWidth();
267 const double minX = xDim->getMinimum();
268 outputWorkspace->setBinEdges(0, xValsSize, HistogramData::LinearGenerator(minX, dx));
269 // set the y-values and errors
270 for (size_t i = 0; i < ny; ++i) {
271 if (i > 0)
272 outputWorkspace->setSharedX(i, outputWorkspace->sharedX(0));
273
274 size_t yOffset = i * yStride;
275 for (size_t j = 0; j < nx; ++j) {
276 size_t linearIndex = yOffset + j * xStride;
277 signal_t signal = inputWorkspace->getSignalArray()[linearIndex];
278 signal_t error = inputWorkspace->getErrorSquaredArray()[linearIndex];
279 // apply normalization
280 if (normalization != NoNormalization) {
281 if (normalization == VolumeNormalization) {
282 signal *= inverseVolume;
283 error *= inverseVolume;
284 } else // normalization == NumEventsNormalization
285 {
286 signal_t factor = inputWorkspace->getNumEventsArray()[linearIndex];
287 factor = factor != 0.0 ? 1.0 / factor : 1.0;
288 signal *= factor;
289 error *= factor;
290 }
291 }
292 outputWorkspace->mutableY(i)[j] = signal;
293 outputWorkspace->mutableE(i)[j] = sqrt(error);
294 }
295 }
296
297 // set the first axis
298 auto labelX = std::dynamic_pointer_cast<Kernel::Units::Label>(Kernel::UnitFactory::Instance().create("Label"));
299 labelX->setLabel(xDim->getName());
300 outputWorkspace->getAxis(0)->unit() = labelX;
301
302 // set the second axis
303 auto yAxis = std::make_unique<BinEdgeAxis>(ny + 1);
304 for (size_t i = 0; i <= ny; ++i) {
305 yAxis->setValue(i, yDim->getX(i));
306 }
307 auto labelY = std::dynamic_pointer_cast<Kernel::Units::Label>(Kernel::UnitFactory::Instance().create("Label"));
308 labelY->setLabel(yDim->getName());
309 yAxis->unit() = labelY;
310 outputWorkspace->replaceAxis(1, std::move(yAxis));
311
312 // set the "units" for the y values
313 outputWorkspace->setYUnitLabel("Signal");
314
315 // done
316 setProperty("OutputWorkspace", outputWorkspace);
317}
318
325 size_t stride = 1;
326 for (size_t i = 0; i < dim; ++i) {
327 auto dimension = workspace.getDimension(i);
328 stride *= dimension->getNBins();
329 }
330 return stride;
331}
332
333} // namespace Mantid::MDAlgorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
double error
Definition: IndexPeaks.cpp:133
IPeaksWorkspace_sptr workspace
Definition: IndexPeaks.cpp:114
#define fabs(x)
Definition: Matrix.cpp:22
std::string getPropertyValue(const std::string &name) const override
Get the value of a property as a string.
Definition: Algorithm.cpp:2026
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
Definition: Algorithm.cpp:2076
Kernel::Logger & g_log
Definition: Algorithm.h:451
Unique SingleValueParameter Declaration for InputNDimensions.
Mantid::Kernel::VMD applyVMD(const Mantid::Kernel::VMD &inputVector) const
Wrapper for VMD.
Abstract interface to MDHistoWorkspace, for use in exposing to Python.
size_t numOriginalWorkspaces() const
Definition: MDGeometry.cpp:327
std::shared_ptr< Workspace > getOriginalWorkspace(size_t index=0) const
Get the "original" workspace (the workspace that was the source for a binned MDWorkspace).
Definition: MDGeometry.cpp:340
A property class for workspaces.
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
ListValidator is a validator that requires the value of a property to be one of a defined list of pos...
Definition: ListValidator.h:29
The Logger class is in charge of the publishing messages from the framework through various channels.
Definition: Logger.h:52
void information(const std::string &msg)
Logs at information level.
Definition: Logger.cpp:105
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...
TYPE normalize()
Normalize this vector to unity length.
Definition: VMD.cpp:427
size_t getNumDims() const
Definition: VMD.cpp:236
Creates a single spectrum Workspace2D with X,Y, and E copied from an first non-integrated dimension o...
size_t calcStride(const API::IMDHistoWorkspace &workspace, size_t dim) const
Calculate the stride for a dimension.
std::shared_ptr< const CoordTransform > CoordTransform_const_sptr
std::shared_ptr< IMDHistoWorkspace > IMDHistoWorkspace_sptr
shared pointer to Mantid::API::IMDHistoWorkspace
MDNormalization
Enum describing different ways to normalize the signal in a MDWorkspace.
Definition: IMDIterator.h:25
@ VolumeNormalization
Divide the signal by the volume of the box/bin.
Definition: IMDIterator.h:29
@ NumEventsNormalization
Divide the signal by the number of events that contributed to it.
Definition: IMDIterator.h:31
@ NoNormalization
Don't normalize = return raw counts.
Definition: IMDIterator.h:27
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
std::unique_ptr< T > create(const P &parent, const IndexArg &indexArg, const HistArg &histArg)
This is the create() method that all the other create() methods call.
std::vector< IMDDimension_const_sptr > VecIMDDimension_const_sptr
Vector of constant shared pointers to IMDDimensions.
Definition: IMDDimension.h:103
std::shared_ptr< const IMDDimension > IMDDimension_const_sptr
Shared Pointer to const IMDDimension.
Definition: IMDDimension.h:101
std::shared_ptr< IValidator > IValidator_sptr
A shared_ptr to an IValidator.
Definition: IValidator.h:26
VMDBase< VMD_t > VMD
Define the VMD as using the double or float data type.
Definition: VMD.h:86
double signal_t
Typedef for the signal recorded in a MDBox, etc.
Definition: MDTypes.h:36
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54