Mantid
Loading...
Searching...
No Matches
ConvertAxisByFormula.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 +
8#include "MantidAPI/RefAxis.h"
15
16#include <algorithm>
17#include <memory>
18#include <sstream>
19
20namespace Mantid::Algorithms {
21
22using namespace Kernel;
23using namespace API;
24
25// Register the class into the algorithm factory
26DECLARE_ALGORITHM(ConvertAxisByFormula)
27
28const std::string ConvertAxisByFormula::name() const { return ("ConvertAxisByFormula"); }
29
30int ConvertAxisByFormula::version() const { return (1); }
31
32const std::string ConvertAxisByFormula::category() const { return "Transforms\\Axes"; }
33
38 declareProperty(std::make_unique<WorkspaceProperty<MatrixWorkspace>>("InputWorkspace", "", Direction::Input),
39 "Name of the input workspace");
40 declareProperty(std::make_unique<WorkspaceProperty<MatrixWorkspace>>("OutputWorkspace", "", Direction::Output),
41 "Name of the output workspace");
42
43 std::vector<std::string> axisOptions;
44 axisOptions.emplace_back("X");
45 axisOptions.emplace_back("Y");
46 declareProperty("Axis", "X", std::make_shared<StringListValidator>(axisOptions), "The axis to modify");
47
48 declareProperty("Formula", "",
49 "The formula to use to convert the values, x "
50 "or y may be used to refer to the axis "
51 "values. l1, l2, twotheta and signedtwotheta"
52 "may be used to provide values from the "
53 "instrument geometry.");
54 declareProperty("AxisTitle", "", "The label of he new axis. If not set then the title will not change.");
55 declareProperty("AxisUnits", "", "The units of the new axis. If not set then the unit will not change");
56}
57
62 // get the property values
63 MatrixWorkspace_sptr inputWs = getProperty("InputWorkspace");
64 std::string axis = getProperty("Axis");
65 std::string formula = getProperty("Formula");
66 std::string axisTitle = getProperty("AxisTitle");
67 std::string axisUnits = getProperty("AxisUnits");
68
69 // Just overwrite if the change is in place
70 MatrixWorkspace_sptr outputWs = getProperty("OutputWorkspace");
71 if (outputWs != inputWs) {
72 auto duplicate = createChildAlgorithm("CloneWorkspace", 0.0, 0.6);
73 duplicate->initialize();
74 duplicate->setProperty<Workspace_sptr>("InputWorkspace", std::dynamic_pointer_cast<Workspace>(inputWs));
75 duplicate->execute();
76 Workspace_sptr temp = duplicate->getProperty("OutputWorkspace");
77 outputWs = std::dynamic_pointer_cast<MatrixWorkspace>(temp);
78
79 setProperty("OutputWorkspace", outputWs);
80 }
81
82 // Get the axes - assume X
83 Axis *axisPtr = outputWs->getAxis(0);
84 Axis *otherAxisPtr = outputWs->getAxis(1);
85 if (axis == "Y") {
86 std::swap(axisPtr, otherAxisPtr);
87 }
88
89 if (!axisPtr->isNumeric()) {
90 throw std::invalid_argument("This algorithm only operates on numeric axes");
91 }
92
93 bool isRaggedBins = false;
94 bool isRefAxis = false;
95 auto *refAxisPtr = dynamic_cast<RefAxis *>(axisPtr);
96 if (refAxisPtr != nullptr) {
97 isRaggedBins = !outputWs->isCommonBins();
98 isRefAxis = true;
99 }
100
101 // validate - if formula contains geometry variables then the other axis must
102 // be a spectraAxis
103 auto *spectrumAxisPtr = dynamic_cast<SpectraAxis *>(otherAxisPtr);
104 std::vector<Variable_ptr> variables;
105 variables.reserve(8);
106 // axis value lookups
107 variables.emplace_back(std::make_shared<Variable>("x", false));
108 variables.emplace_back(std::make_shared<Variable>("X", false));
109 variables.emplace_back(std::make_shared<Variable>("y", false));
110 variables.emplace_back(std::make_shared<Variable>("Y", false));
111 // geometry lookups
112 variables.emplace_back(std::make_shared<Variable>("twotheta", true));
113 variables.emplace_back(std::make_shared<Variable>("signedtwotheta", true));
114 variables.emplace_back(std::make_shared<Variable>("l1", true));
115 variables.emplace_back(std::make_shared<Variable>("l2", true));
116
117 bool isGeometryRequired = false;
118 for (auto variablesIter = variables.begin(); variablesIter != variables.end();) {
119 if (formula.find((*variablesIter)->name) != std::string::npos) {
120 if ((*variablesIter)->isGeometric) {
121 if (!spectrumAxisPtr) {
122 throw std::invalid_argument("To use geometry values in the equation, "
123 "the axis not being converted must be a "
124 "Spectra Axis.");
125 }
126 isGeometryRequired = true;
127 }
128 ++variablesIter;
129 } else {
130 // remove those that are not used
131 variablesIter = variables.erase(variablesIter);
132 }
133 }
134
135 // Create muparser
136 mu::Parser p;
137 try {
138 // set parameter lookups for the axis value
139 for (const auto &variable : variables) {
140 p.DefineVar(variable->name, &(variable->value));
141 }
142 // set some constants
143 p.DefineConst("pi", M_PI);
144 p.DefineConst("h", PhysicalConstants::h);
145 p.DefineConst("h_bar", PhysicalConstants::h_bar);
146 p.DefineConst("g", PhysicalConstants::g);
147 p.DefineConst("mN", PhysicalConstants::NeutronMass);
148 p.DefineConst("mNAMU", PhysicalConstants::NeutronMassAMU);
149 p.SetExpr(formula);
150 } catch (mu::Parser::exception_type &e) {
151 std::stringstream ss;
152 ss << "Cannot process the formula"
153 << ". Muparser error message is: " << e.GetMsg();
154 throw std::invalid_argument(ss.str());
155 }
156 if (isRefAxis) {
157 if ((isRaggedBins) || (isGeometryRequired)) {
158 // ragged bins or geometry used - we have to calculate for every spectra
159 size_t numberOfSpectra_i = outputWs->getNumberHistograms();
160 auto &spectrumInfo = outputWs->mutableSpectrumInfo();
161
162 size_t failedDetectorCount = 0;
163 Progress prog(this, 0.6, 1.0, numberOfSpectra_i);
164 for (size_t i = 0; i < numberOfSpectra_i; ++i) {
165 try {
166 MantidVec &vec = outputWs->dataX(i);
167 setGeometryValues(spectrumInfo, i, variables);
168 calculateValues(p, vec, variables);
169 } catch (std::runtime_error &)
170 // two possible exceptions runtime error and NotFoundError
171 // both handled the same way
172 {
173 // could not find the geometry info for this spectra
174 outputWs->getSpectrum(i).clearData();
175 spectrumInfo.setMasked(i, true);
176 failedDetectorCount++;
177 }
178 prog.report();
179 }
180 if (failedDetectorCount != 0) {
181 g_log.information() << "Unable to calculate sample-detector distance for " << failedDetectorCount
182 << " spectra. Masking spectrum.\n";
183 }
184 } else {
185 // common bins - we only have to calculate once
186
187 // Calculate the new (common) X values
188 MantidVec &vec = outputWs->dataX(0);
189 calculateValues(p, vec, variables);
190
191 // copy xVals to every spectra
192 auto numberOfSpectra_i = static_cast<int64_t>(outputWs->getNumberHistograms()); // cast to make openmp happy
193 auto xVals = outputWs->refX(0);
194 Progress prog(this, 0.6, 1.0, numberOfSpectra_i);
196 for (int64_t j = 1; j < numberOfSpectra_i; ++j) {
198 outputWs->setX(j, xVals);
199 prog.report();
201 }
203 }
204 } else {
205 size_t axisLength = axisPtr->length();
206 for (size_t i = 0; i < axisLength; ++i) {
207 setAxisValue(axisPtr->getValue(i), variables);
208 axisPtr->setValue(i, evaluateResult(p));
209 }
210 }
211
212 // If the units conversion has flipped the ascending direction of X, reverse
213 // all the vectors
214 size_t midSpectra = outputWs->getNumberHistograms() / 2;
215 if (!outputWs->dataX(0).empty() && (outputWs->dataX(0).front() > outputWs->dataX(0).back() ||
216 outputWs->dataX(midSpectra).front() > outputWs->dataX(midSpectra).back())) {
217 g_log.information("Reversing data within the workspace to keep the axes in "
218 "increasing order.");
219 this->reverse(outputWs);
220 }
221
222 // Set the Unit of the Axis
223 if ((!axisUnits.empty()) || (!axisTitle.empty())) {
224 try {
225 axisPtr->unit() = UnitFactory::Instance().create(axisUnits);
226 } catch (Exception::NotFoundError &) {
227 if (axisTitle.empty()) {
228 axisTitle = axisPtr->unit()->caption();
229 }
230 if (axisUnits.empty()) {
231 axisUnits = axisPtr->unit()->label();
232 }
233 axisPtr->unit() = std::make_shared<Units::Label>(axisTitle, axisUnits);
234 }
235 }
236}
237
238void ConvertAxisByFormula::setAxisValue(const double value, const std::vector<Variable_ptr> &variables) {
239 for (const auto &variable : variables) {
240 if (!variable->isGeometric) {
241 variable->value = value;
242 }
243 }
244}
245
246void ConvertAxisByFormula::calculateValues(mu::Parser &p, MantidVec &vec, std::vector<Variable_ptr> variables) {
247 MantidVec::iterator iter;
248 for (iter = vec.begin(); iter != vec.end(); ++iter) {
249 setAxisValue(*iter, variables);
250 *iter = evaluateResult(p);
251 }
252}
253
255 const std::vector<Variable_ptr> &variables) {
256 for (const auto &variable : variables) {
257 if (variable->isGeometric) {
258 if (variable->name == "twotheta") {
259 variable->value = specInfo.twoTheta(index);
260 } else if (variable->name == "signedtwotheta") {
261 variable->value = specInfo.signedTwoTheta(index);
262 } else if (variable->name == "l1") {
263 variable->value = specInfo.l1();
264 } else if (variable->name == "l2") {
265 variable->value = specInfo.l2(index);
266 }
267 }
268 }
269}
270
272 double result;
273 try {
274 result = p.Eval();
275 } catch (mu::Parser::exception_type &e) {
276 std::stringstream ss;
277 ss << "Failed while processing axis values"
278 << ". Muparser error message is: " << e.GetMsg();
279 throw std::invalid_argument(ss.str());
280 }
281 return result;
282}
283
284} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
double value
The value of the point.
Definition: FitMW.cpp:51
std::map< DeltaEMode::Type, std::string > index
Definition: DeltaEMode.cpp:19
#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_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
Class to represent the axis of a workspace.
Definition: Axis.h:30
virtual std::size_t length() const =0
Get the length of the axis.
virtual void setValue(const std::size_t &index, const double &value)=0
Sets the value at the specified index.
virtual bool isNumeric() const
Returns true if the axis is numeric.
Definition: Axis.h:52
const std::shared_ptr< Kernel::Unit > & unit() const
The unit for this axis.
Definition: Axis.cpp:28
double getValue(const std::size_t &index, const std::size_t &verticalIndex=0) const
Gets the value at the specified index.
Definition: Axis.cpp:51
Helper class for reporting progress from algorithms.
Definition: Progress.h:25
A class to represent the axis of a 2D (or more) workspace where the value at a given point on the axi...
Definition: RefAxis.h:23
Class to represent the spectra axis of a workspace.
Definition: SpectraAxis.h:31
API::SpectrumInfo is an intermediate step towards a SpectrumInfo that is part of Instrument-2....
Definition: SpectrumInfo.h:53
double signedTwoTheta(const size_t index) const
Returns the signed scattering angle 2 theta in radians (angle w.r.t.
double twoTheta(const size_t index) const
Returns the scattering angle 2 theta in radians (angle w.r.t.
double l2(const size_t index) const
Returns L2 (distance from sample to spectrum).
double l1() const
Returns L1 (distance from source to sample).
A property class for workspaces.
ConvertAxisByFormula : Performs a unit conversion based on a supplied formula.
void setAxisValue(const double value, const std::vector< Variable_ptr > &variables)
const std::string category() const override
function to return a category of the algorithm.
int version() const override
function to return a version of the algorithm, must be overridden in all algorithms
void calculateValues(mu::Parser &p, std::vector< double > &vec, std::vector< Variable_ptr > variables)
void init() override
Initialisation method.
void setGeometryValues(const API::SpectrumInfo &specInfo, const size_t index, const std::vector< Variable_ptr > &variables)
void exec() override
Execution of the algorithm.
void reverse(const API::MatrixWorkspace_sptr &WS)
Reverses the workspace if X values are in descending order.
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 information(const std::string &msg)
Logs at information level.
Definition: Logger.cpp:105
void report()
Increments the loop counter by 1, then sends the progress notification on behalf of its algorithm.
Definition: ProgressBase.h:51
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
std::shared_ptr< Workspace > Workspace_sptr
shared pointer to Mantid::API::Workspace
Definition: Workspace_fwd.h:20
Kernel::Logger g_log("ExperimentInfo")
static logger object
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
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
static constexpr double NeutronMassAMU
Mass of the neutron in AMU.
static constexpr double NeutronMass
Mass of the neutron in kg.
static constexpr double h
Planck constant in J*s.
static constexpr double h_bar
Planck constant in J*s, divided by 2 PI.
static constexpr double g
Standard acceleration due to gravity.
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