Mantid
Loading...
Searching...
No Matches
ConvertAxesToRealSpace.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 +
14#include "MantidKernel/Unit.h"
16
17#include <limits>
18
19namespace Mantid::Algorithms {
20
21using namespace Mantid::Kernel;
22using namespace Mantid::API;
23using namespace Mantid::DataObjects;
24
25// Register the algorithm into the AlgorithmFactory
26DECLARE_ALGORITHM(ConvertAxesToRealSpace)
27
28
29const std::string ConvertAxesToRealSpace::name() const { return "ConvertAxesToRealSpace"; }
30
32int ConvertAxesToRealSpace::version() const { return 1; }
33
35const std::string ConvertAxesToRealSpace::category() const { return "Transforms\\Units;Transforms\\Axes"; }
36
38const std::string ConvertAxesToRealSpace::summary() const {
39 return "Converts the spectrum and TOF axes to real space values, integrating "
40 "the data in the process";
41}
42
46 declareProperty(std::make_unique<WorkspaceProperty<MatrixWorkspace>>("InputWorkspace", "", Direction::Input),
47 "An input workspace.");
48 declareProperty(std::make_unique<WorkspaceProperty<Workspace2D>>("OutputWorkspace", "", Direction::Output),
49 "An output workspace.");
50
51 std::vector<std::string> propOptions;
52 fillUnitMap(propOptions, m_unitMap, "x", "m");
53 fillUnitMap(propOptions, m_unitMap, "y", "m");
54 fillUnitMap(propOptions, m_unitMap, "z", "m");
55 fillUnitMap(propOptions, m_unitMap, "r", "m");
56 fillUnitMap(propOptions, m_unitMap, "theta", "deg");
57 fillUnitMap(propOptions, m_unitMap, "phi", "deg");
58 fillUnitMap(propOptions, m_unitMap, "2theta", "rad");
59 fillUnitMap(propOptions, m_unitMap, "signed2theta", "rad");
60
61 declareProperty("VerticalAxis", "y", std::make_shared<StringListValidator>(propOptions),
62 "What will be the vertical axis ?\n");
63 declareProperty("HorizontalAxis", "2theta", std::make_shared<StringListValidator>(propOptions),
64 "What will be the horizontal axis?\n");
65
66 declareProperty(std::make_unique<Kernel::PropertyWithValue<int>>("NumberVerticalBins", 100),
67 "The number of bins along the vertical axis.");
68 declareProperty(std::make_unique<Kernel::PropertyWithValue<int>>("NumberHorizontalBins", 100),
69 "The number of bins along the horizontal axis.");
70}
71
75 MatrixWorkspace_sptr inputWs = getProperty("InputWorkspace");
76
77 // set up axes data
78 std::vector<AxisData> axisVector = std::vector<AxisData>(2);
79 std::string hAxis = getProperty("HorizontalAxis");
80 axisVector[0].label = hAxis;
81 axisVector[0].bins = getProperty("NumberHorizontalBins");
82 std::string vAxis = getProperty("VerticalAxis");
83 axisVector[1].label = vAxis;
84 axisVector[1].bins = getProperty("NumberVerticalBins");
85 for (int axisIndex = 0; axisIndex < 2; ++axisIndex) {
86 axisVector[axisIndex].max = std::numeric_limits<double>::lowest();
87 axisVector[axisIndex].min = std::numeric_limits<double>::max();
88 }
89
90 // Create the output workspace. Can't re-use the input one because we'll be
91 // re-ordering the spectra.
92 MatrixWorkspace_sptr outputWs =
93 WorkspaceFactory::Instance().create(inputWs, axisVector[1].bins, axisVector[0].bins, axisVector[0].bins);
94
95 // first integrate the data
96 auto alg = createChildAlgorithm("Integration", 0, 0.4);
97 alg->setProperty<MatrixWorkspace_sptr>("InputWorkspace", inputWs);
98 std::string outName = "_" + inputWs->getName() + "_integrated";
99 alg->setProperty("OutputWorkspace", outName);
100 alg->executeAsChildAlg();
101 MatrixWorkspace_sptr summedWs = alg->getProperty("OutputWorkspace");
102
103 auto nHist = static_cast<int>(summedWs->getNumberHistograms());
104 Progress progress(this, 0.4, 1.0, nHist * 4);
105
106 std::vector<SpectraData> dataVector(nHist);
107
108 int failedCount = 0;
109
110 const auto &spectrumInfo = summedWs->spectrumInfo();
111 // for each spectra
112 PARALLEL_FOR_IF(Kernel::threadSafe(*summedWs, *outputWs))
113 for (int i = 0; i < nHist; ++i) {
115 try {
116 V3D pos = spectrumInfo.position(i);
117 double r, theta, phi;
118 pos.getSpherical(r, theta, phi);
119
120 // for each axis
121 for (int axisIndex = 0; axisIndex < 2; ++axisIndex) {
122 double axisValue = std::numeric_limits<double>::min();
123 std::string axisSelection = axisVector[axisIndex].label;
124
125 // get the selected value for this axis
126 if (axisSelection == "x") {
127 axisValue = pos.X();
128 } else if (axisSelection == "y") {
129 axisValue = pos.Y();
130 } else if (axisSelection == "z") {
131 axisValue = pos.Z();
132 } else if (axisSelection == "r") {
133 axisValue = r;
134 } else if (axisSelection == "theta") {
135 axisValue = theta;
136 } else if (axisSelection == "phi") {
137 axisValue = phi;
138 } else if (axisSelection == "2theta") {
139 axisValue = spectrumInfo.twoTheta(i);
140 } else if (axisSelection == "signed2theta") {
141 axisValue = spectrumInfo.signedTwoTheta(i);
142 }
143
144 if (axisIndex == 0) {
145 dataVector[i].horizontalValue = axisValue;
146 } else {
147 dataVector[i].verticalValue = axisValue;
148 }
149
150 // record the max and min values
151 if (axisValue > axisVector[axisIndex].max)
152 axisVector[axisIndex].max = axisValue;
153 if (axisValue < axisVector[axisIndex].min)
154 axisVector[axisIndex].min = axisValue;
155 }
156 } catch (const Exception::NotFoundError &) {
157 g_log.debug() << "Could not find detector for workspace index " << i << '\n';
158 failedCount++;
159 // flag this is the datavector
160 dataVector[i].horizontalValue = std::numeric_limits<double>::min();
161 dataVector[i].verticalValue = std::numeric_limits<double>::min();
162 }
164
165 // take the values from the integrated data
166 dataVector[i].intensity = summedWs->y(i)[0];
167 dataVector[i].error = summedWs->e(i)[0];
168
169 progress.report("Calculating new coords");
170 }
172
173 g_log.warning() << "Could not find detector for " << failedCount << " spectra, see the debug log for more details.\n";
174
175 // set up the axes on the output workspace
176 std::vector<double> x_tmp(axisVector[0].bins);
178 fillAxisValues(x_tmp, axisVector[0], false);
179 HistogramData::Points x(std::move(x_tmp));
180
181 outputWs->getAxis(0)->unit() = UnitFactory::Instance().create("Label");
182 Unit_sptr xUnit = outputWs->getAxis(0)->unit();
183 std::shared_ptr<Units::Label> xlabel = std::dynamic_pointer_cast<Units::Label>(xUnit);
184 xlabel->setLabel(axisVector[0].label, m_unitMap[axisVector[0].label]);
185
186 MantidVec &yRef = y.access();
187 yRef.resize(axisVector[1].bins);
188 fillAxisValues(yRef, axisVector[1], false);
189
190 auto yAxis = std::make_unique<NumericAxis>(yRef);
191 std::shared_ptr<Units::Label> ylabel =
192 std::dynamic_pointer_cast<Units::Label>(UnitFactory::Instance().create("Label"));
193 ylabel->setLabel(axisVector[1].label, m_unitMap[axisVector[1].label]);
194 yAxis->unit() = ylabel;
195 outputWs->replaceAxis(1, std::move(yAxis));
196
197 // work out where to put the data into the output workspace, but don't do it
198 // yet as that needs to be single threaded
200 for (int i = 0; i < nHist; ++i) {
201 // find write index for data point
202 if (dataVector[i].horizontalValue == std::numeric_limits<double>::min()) {
203 dataVector[i].horizontalIndex = -1;
204 dataVector[i].verticalIndex = -1;
205 } else {
206 int xIndex = static_cast<int>(
207 std::distance(x.cbegin(), std::lower_bound(x.cbegin(), x.cend(), dataVector[i].horizontalValue)));
208 if (xIndex > 0)
209 --xIndex;
210 int yIndex = static_cast<int>(
211 std::distance(y->begin(), std::lower_bound(y->begin(), y->end(), dataVector[i].verticalValue)));
212 if (yIndex > 0)
213 --yIndex;
214
215 dataVector[i].horizontalIndex = xIndex;
216 dataVector[i].verticalIndex = yIndex;
217 }
218 progress.report("Calculating Rebinning");
219 }
220
221 // set all the X arrays - share the same vector
222 auto nOutputHist = static_cast<int>(outputWs->getNumberHistograms());
224 for (int i = 0; i < nOutputHist; ++i) {
225 outputWs->setPoints(i, x);
226 }
227
228 // insert the data into the new workspace
229 // single threaded
230 for (int i = 0; i < nHist; ++i) {
231 int xIndex = dataVector[i].horizontalIndex;
232 int yIndex = dataVector[i].verticalIndex;
233
234 // using -1 as a flag for could not find detector
235 if ((xIndex == -1) || (yIndex == -1)) {
236 // do nothing the detector could not be found
237 g_log.warning() << "here " << i << '\n';
238 } else {
239 // update the data
240 auto &yVec = outputWs->mutableY(yIndex);
241 yVec[xIndex] = yVec[xIndex] + dataVector[i].intensity;
242 auto &eVec = outputWs->mutableE(yIndex);
243 eVec[xIndex] = eVec[xIndex] + (dataVector[i].error * dataVector[i].error);
244 }
245
246 progress.report("Assigning to new grid");
247 }
248
249 // loop over the data and sqrt the errors to complete the error calculation
251 for (int i = 0; i < nOutputHist; ++i) {
252 auto &errorVec = outputWs->mutableE(i);
253 std::transform(errorVec.begin(), errorVec.end(), errorVec.begin(), static_cast<double (*)(double)>(sqrt));
254 progress.report("Completing Error Calculation");
255 }
256
257 // Execute the transform and bind to the output.
258 setProperty("OutputWorkspace", outputWs);
259}
260
268void ConvertAxesToRealSpace::fillAxisValues(MantidVec &vector, const AxisData &axisData, bool isHistogram) {
269 int numBins = axisData.bins;
270 double binDelta = (axisData.max - axisData.min) / static_cast<double>(numBins);
271
272 if (isHistogram)
273 numBins++;
274
275 for (int i = 0; i < numBins; ++i) {
276 vector[i] = axisData.min + i * binDelta;
277 }
278}
279
286void ConvertAxesToRealSpace::fillUnitMap(std::vector<std::string> &orderedVector,
287 std::map<std::string, std::string> &unitMap, const std::string &caption,
288 const std::string &unit) {
289 unitMap.emplace(caption, unit);
290 orderedVector.emplace_back(caption);
291}
292
293} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
#define PARALLEL_START_INTERRUPT_REGION
Begins a block to skip processing is the algorithm has been interupted Note the end of the block if n...
Definition: MultiThreaded.h:94
#define PARALLEL_FOR_NO_WSP_CHECK()
#define PARALLEL_END_INTERRUPT_REGION
Ends a block to skip processing is the algorithm has been interupted Note the start of the block if n...
#define PARALLEL_FOR_IF(condition)
Empty definitions - to enable set your complier to enable openMP.
#define PARALLEL_CHECK_INTERRUPT_REGION
Adds a check after a Parallel region to see if it was interupted.
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
Kernel::Logger & g_log
Definition: Algorithm.h:451
void progress(double p, const std::string &msg="", double estimatedTime=0.0, int progressPrecision=0)
Sends ProgressNotification.
Definition: Algorithm.cpp:231
Helper class for reporting progress from algorithms.
Definition: Progress.h:25
A property class for workspaces.
ConvertAxesToRealSpace : Converts the spectrum and TOF axes to real space values, integrating the dat...
std::map< std::string, std::string > m_unitMap
void exec() override
Execute the algorithm.
const std::string category() const override
Algorithm's category for identification.
void fillUnitMap(std::vector< std::string > &orderedVector, std::map< std::string, std::string > &unitMap, const std::string &caption, const std::string &unit)
Fills the unit map and ordered vector with the same data.
int version() const override
Algorithm's version for identification.
void init() override
Initialize the algorithm's properties.
void fillAxisValues(MantidVec &vector, const AxisData &axisData, bool isHistogram)
Fills the values in an axis linearly from min to max for a given number of steps.
const std::string summary() const override
Algorithm's summary for use in the GUI and help.
Exception for when an item is not found in a collection.
Definition: Exception.h:145
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void debug(const std::string &msg)
Logs at debug level.
Definition: Logger.cpp:114
void warning(const std::string &msg)
Logs at warning level.
Definition: Logger.cpp:86
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...
Class for 3D vectors.
Definition: V3D.h:34
constexpr double X() const noexcept
Get x.
Definition: V3D.h:232
constexpr double Y() const noexcept
Get y.
Definition: V3D.h:233
void getSpherical(double &R, double &theta, double &phi) const noexcept
Return the vector's position in spherical coordinates.
Definition: V3D.cpp:117
constexpr double Z() const noexcept
Get z.
Definition: V3D.h:234
Implements a copy on write data template.
Definition: cow_ptr.h:41
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::shared_ptr< Unit > Unit_sptr
Shared pointer to the Unit base class.
Definition: Unit.h:229
std::enable_if< std::is_pointer< Arg >::value, bool >::type threadSafe(Arg workspace)
Thread-safety check Checks the workspace to ensure it is suitable for multithreaded access.
Definition: MultiThreaded.h:22
std::vector< double > MantidVec
typedef for the data storage used in Mantid matrix workspaces
Definition: cow_ptr.h:172
STL namespace.
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54