Mantid
Loading...
Searching...
No Matches
IQTransform.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//----------------------------------------------------------------------
10#include <utility>
11
12#include "MantidAPI/Axis.h"
21#include "MantidHistogramData/Histogram.h"
26#include "MantidKernel/Unit.h"
28
29using namespace Mantid::DataObjects;
30using namespace Mantid::HistogramData;
31
32namespace Mantid::Algorithms {
33
34// Register the algorithm into the AlgorithmFactory
35DECLARE_ALGORITHM(IQTransform)
36
37using namespace Kernel;
38using namespace API;
39
40IQTransform::IQTransform() : API::Algorithm(), m_label(std::make_shared<Units::Label>()) {
41 /* Just for fun, this is implemented as follows....
42 * We fill a map below with the transformation name as the key and
43 * a function pointer to the method that does the transformation as
44 * the value. The 'TransformType' property is filled with the keys
45 * and then we search on that to select the correct function in exec.
46 */
47 m_transforms["Guinier (spheres)"] = &IQTransform::guinierSpheres;
48 m_transforms["Guinier (rods)"] = &IQTransform::guinierRods;
49 m_transforms["Guinier (sheets)"] = &IQTransform::guinierSheets;
51 m_transforms["Debye-Bueche"] = &IQTransform::debyeBueche;
57}
58
60 auto wsValidator = std::make_shared<CompositeValidator>();
61 // Require the input to be in units of Q and to be a distribution
62 // (which the result of a SANS reduction in Mantid will be)
63 wsValidator->add<WorkspaceUnitValidator>("MomentumTransfer");
64 wsValidator->add<RawCountValidator>(false);
65 // Require X data to be increasing from left to right
66 wsValidator->add<IncreasingAxisValidator>();
67
68 declareProperty(std::make_unique<WorkspaceProperty<>>("InputWorkspace", "", Direction::Input, wsValidator),
69 "The input workspace must be a distribution with units of Q");
70 declareProperty(std::make_unique<WorkspaceProperty<>>("OutputWorkspace", "", Direction::Output),
71 "The name of the output workspace");
72
73 // Extract the keys from the transformations map to pass to the property
74 std::set<std::string> plottype;
75 for (TransformMap::const_iterator it = m_transforms.begin(); it != m_transforms.end(); ++it) {
76 plottype.insert(it->first);
77 }
78 declareProperty("TransformType", "", std::make_shared<StringListValidator>(plottype),
79 "The name of the transformation to be performed on the workspace");
80
81 // A background to be subtracted can be a value or a workspace. Both
82 // properties are optional.
83 auto mustBePositive = std::make_shared<BoundedValidator<double>>();
84 mustBePositive->setLower(0.0);
85 declareProperty("BackgroundValue", 0.0, mustBePositive,
86 "A constant value to subtract from the data prior to its transformation");
88 std::make_unique<WorkspaceProperty<>>("BackgroundWorkspace", "", Direction::Input, PropertyMode::Optional),
89 "A workspace to subtract from the input workspace prior to its "
90 "transformation."
91 "Must be compatible with the input (as for the Minus algorithm).");
92
93 declareProperty(std::make_unique<ArrayProperty<double>>("GeneralFunctionConstants"),
94 "A set of 10 constants to be used (only) with the 'General' "
95 "transformation");
96}
97
99 MatrixWorkspace_sptr inputWS = getProperty("InputWorkspace");
100 // Print a warning if the input workspace has more than one spectrum
101 if (inputWS->getNumberHistograms() > 1) {
102 g_log.warning("This algorithm is intended for use on single-spectrum workspaces.\n"
103 "Only the first spectrum will be transformed.");
104 }
105
106 // Do background subtraction from a workspace first because it doesn't like
107 // potential conversion to point data that follows. Requires a temporary
108 // workspace.
110 MatrixWorkspace_sptr backgroundWS = getProperty("BackgroundWorkspace");
111 if (backgroundWS)
112 tmpWS = subtractBackgroundWS(inputWS, backgroundWS);
113 else
114 tmpWS = inputWS;
115
116 // Create the output workspace
117 const size_t length = tmpWS->blocksize();
118 MatrixWorkspace_sptr outputWS = create<MatrixWorkspace>(*inputWS, 1, Points(length));
119 m_label->setLabel("");
120 outputWS->setYUnit("");
121 // Copy the data over. Assume single spectrum input (output will be).
122 outputWS->setPoints(0, tmpWS->points(0));
123 outputWS->setSharedY(0, tmpWS->sharedY(0));
124 outputWS->setSharedE(0, tmpWS->sharedE(0));
125
126 // Subtract a constant background if requested
127 const double background = getProperty("BackgroundValue");
128 if (background > 0.0)
129 outputWS->mutableY(0) -= background;
130
131 // Select the desired transformation function and call it
132 TransformFunc f = m_transforms.find(getProperty("TransformType"))->second;
133 (this->*f)(outputWS);
134
135 // Need the generic label unit on this (unless the unit on the X axis hasn't
136 // changed)
137 if (!m_label->caption().empty())
138 outputWS->getAxis(0)->unit() = m_label;
139 setProperty("OutputWorkspace", outputWS);
140}
141
151 const API::MatrixWorkspace_sptr &background) {
152 g_log.debug() << "Subtracting the workspace " << background->getName() << " from the input workspace.\n";
153 return ws - background;
154}
155
158
165 auto &X = ws->mutableX(0);
166 auto &Y = ws->mutableY(0);
167 auto &E = ws->mutableE(0);
168 std::transform(X.cbegin(), X.cend(), X.begin(), VectorHelper::Squares<double>());
169 std::transform(E.cbegin(), E.cend(), Y.begin(), E.begin(), std::divides<double>());
170 std::transform(Y.cbegin(), Y.cend(), Y.begin(), VectorHelper::LogNoThrow<double>());
171
172 ws->setYUnitLabel("Ln(I)");
173 m_label->setLabel("Q^2");
174}
175
182 auto &X = ws->mutableX(0);
183 auto &Y = ws->mutableY(0);
184 auto &E = ws->mutableE(0);
185
186 std::transform(E.cbegin(), E.cend(), Y.begin(), E.begin(), std::divides<double>());
187 std::transform(Y.cbegin(), Y.cend(), X.begin(), Y.begin(), std::multiplies<double>());
188 std::transform(Y.cbegin(), Y.cend(), Y.begin(), VectorHelper::LogNoThrow<double>());
189 std::transform(X.cbegin(), X.cend(), X.begin(), VectorHelper::Squares<double>());
190
191 ws->setYUnitLabel("Ln(I x Q)");
192 m_label->setLabel("Q^2");
193}
194
201 auto &X = ws->mutableX(0);
202 auto &Y = ws->mutableY(0);
203 auto &E = ws->mutableE(0);
204
205 std::transform(E.cbegin(), E.cend(), Y.begin(), E.begin(), std::divides<double>());
206 std::transform(X.cbegin(), X.cend(), X.begin(), VectorHelper::Squares<double>());
207 std::transform(Y.cbegin(), Y.cend(), X.begin(), Y.begin(), std::multiplies<double>());
208 std::transform(Y.cbegin(), Y.cend(), Y.begin(), VectorHelper::LogNoThrow<double>());
209
210 ws->setYUnitLabel("Ln(I x Q^2)");
211 m_label->setLabel("Q^2");
212}
213
219 auto &X = ws->mutableX(0);
220 auto &Y = ws->mutableY(0);
221 auto &E = ws->mutableE(0);
222 std::transform(X.cbegin(), X.cend(), X.begin(), VectorHelper::Squares<double>());
223 for (size_t i = 0; i < Y.size(); ++i) {
224 if (Y[i] > 0.0) {
225 Y[i] = 1.0 / Y[i];
226 E[i] *= std::pow(Y[i], 2);
227 } else {
228 Y[i] = 0.0;
229 E[i] = 0.0;
230 }
231 }
232
233 ws->setYUnitLabel("1/I");
234 m_label->setLabel("Q^2");
235}
236
242 auto &X = ws->mutableX(0);
243 auto &Y = ws->mutableY(0);
244 auto &E = ws->mutableE(0);
245 std::transform(X.cbegin(), X.cend(), X.begin(), VectorHelper::Squares<double>());
246 for (size_t i = 0; i < Y.size(); ++i) {
247 if (Y[i] > 0.0) {
248 Y[i] = 1.0 / std::sqrt(Y[i]);
249 E[i] *= std::pow(Y[i], 3);
250 } else {
251 Y[i] = 0.0;
252 E[i] = 0.0;
253 }
254 }
255
256 ws->setYUnitLabel("1/sqrt(I)");
257 m_label->setLabel("Q^2");
258}
259
264 auto &X = ws->mutableX(0);
265 auto &Y = ws->mutableY(0);
266 auto &E = ws->mutableE(0);
267 std::transform(Y.cbegin(), Y.cend(), X.begin(), Y.begin(), std::multiplies<double>());
268 std::transform(E.cbegin(), E.cend(), X.begin(), E.begin(), std::multiplies<double>());
269
270 ws->setYUnitLabel("I x Q");
271}
272
277 auto &X = ws->mutableX(0);
278 auto &Y = ws->mutableY(0);
279 auto &E = ws->mutableE(0);
280 MantidVec Q2(X.size());
281 std::transform(X.cbegin(), X.cend(), Q2.begin(), VectorHelper::Squares<double>());
282 std::transform(Y.cbegin(), Y.cend(), Q2.begin(), Y.begin(), std::multiplies<double>());
283 std::transform(E.cbegin(), E.cend(), Q2.begin(), E.begin(), std::multiplies<double>());
284
285 ws->setYUnitLabel("I x Q^2");
286}
287
292 auto &X = ws->mutableX(0);
293 auto &Y = ws->mutableY(0);
294 auto &E = ws->mutableE(0);
295 MantidVec Q4(X.size());
296 std::transform(X.cbegin(), X.cend(), X.cbegin(), Q4.begin(), VectorHelper::TimesSquares<double>());
297 std::transform(Y.cbegin(), Y.cend(), Q4.begin(), Y.begin(), std::multiplies<double>());
298 std::transform(E.cbegin(), E.cend(), Q4.begin(), E.begin(), std::multiplies<double>());
299
300 ws->setYUnitLabel("I x Q^4");
301}
302
309 auto &X = ws->mutableX(0);
310 auto &Y = ws->mutableY(0);
311 auto &E = ws->mutableE(0);
312
313 std::transform(X.cbegin(), X.cend(), X.begin(), VectorHelper::Log<double>());
314 std::transform(E.cbegin(), E.cend(), Y.begin(), E.begin(), std::divides<double>());
315 std::transform(Y.cbegin(), Y.cend(), Y.begin(), VectorHelper::LogNoThrow<double>());
316
317 ws->setYUnitLabel("Ln(I)");
318 m_label->setLabel("Ln(Q)");
319}
320
330 auto &X = ws->mutableX(0);
331 auto &Y = ws->mutableY(0);
332 auto &E = ws->mutableE(0);
333 const std::vector<double> C = getProperty("GeneralFunctionConstants");
334 // Check for the correct number of elements
335 if (C.size() != 10) {
336 std::string mess("The General transformation requires 10 values to be provided.");
337 g_log.error(mess);
338 throw std::invalid_argument(mess);
339 }
340
341 for (size_t i = 0; i < Y.size(); ++i) {
342 double tmpX = std::pow(X[i], C[7]) * std::pow(Y[i], C[8]) * C[9];
343 if (tmpX <= 0.0)
344 throw std::range_error("Attempt to take log of a zero or negative number.");
345 tmpX = std::pow(X[i], C[5]) * std::pow(Y[i], C[6]) * std::log(tmpX);
346 const double tmpY = std::pow(X[i], C[2]) * std::pow(Y[i], C[3]) * C[4];
347 if (tmpY <= 0.0)
348 throw std::range_error("Attempt to take log of a zero or negative number.");
349 const double newY = std::pow(X[i], C[0]) * std::pow(Y[i], C[1]) * std::log(tmpY);
350
351 E[i] *= std::pow(X[i], C[0]) *
352 (C[1] * std::pow(Y[i], C[1] - 1) * std::log(tmpY) +
353 ((std::pow(Y[i], C[1]) * std::pow(X[i], C[2]) * C[4] * C[3] * std::pow(Y[i], C[3] - 1)) / tmpY));
354 X[i] = tmpX;
355 Y[i] = newY;
356 }
357
358 std::stringstream ylabel;
359 ylabel << "Q^" << C[0] << " x I^" << C[1] << " x Ln( Q^" << C[2] << " x I^" << C[3] << " x " << C[4] << ")";
360 ws->setYUnitLabel(ylabel.str());
361 std::stringstream xlabel;
362 xlabel << "Q^" << C[5] << " x I^" << C[6] << " x Ln( Q^" << C[7] << " x I^" << C[8] << " x " << C[9] << ")";
363 m_label->setLabel(xlabel.str());
364}
365
367
368} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
Base class from which all concrete algorithm classes should be derived.
Definition: Algorithm.h:85
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
Kernel::Logger & g_log
Definition: Algorithm.h:451
A validator which checks that the X axis of a workspace is increasing from left to right.
A validator which checks that a workspace contains raw counts in its bins.
A property class for workspaces.
A validator which checks that the unit of the workspace referred to by a WorkspaceProperty is the exp...
API::MatrixWorkspace_sptr subtractBackgroundWS(const API::MatrixWorkspace_sptr &ws, const API::MatrixWorkspace_sptr &background)
Uses the Minus algorithm to subtract the background workspace from the given workspace.
void logLog(const API::MatrixWorkspace_sptr &ws)
Performs a log-log transformation: Ln(I) v Ln(Q)
void kratky(const API::MatrixWorkspace_sptr &ws)
Performs the Kratky transformation: IQ^2 v Q.
void exec() override
Virtual method - must be overridden by concrete algorithm.
Definition: IQTransform.cpp:98
void porod(const API::MatrixWorkspace_sptr &ws)
Performs the Porod transformation: IQ^4 v Q.
void guinierSheets(const API::MatrixWorkspace_sptr &ws)
Performs the Guinier (sheets) transformation: Ln(IQ^2) v Q^2.
void guinierSpheres(const API::MatrixWorkspace_sptr &ws)
Performs the Guinier (spheres) transformation: Ln(I) v Q^2.
void general(const API::MatrixWorkspace_sptr &ws)
Performs a transformation of the form: Q^A x I^B x Ln(Q^C x I^D x E) v Q^F x I^G x Ln(Q^H x I^I x J).
void guinierRods(const API::MatrixWorkspace_sptr &ws)
Performs the Guinier (rods) transformation: Ln(IQ) v Q^2.
void zimm(const API::MatrixWorkspace_sptr &ws)
Performs the Zimm transformation: 1/I v Q^2 The output is set to zero for negative input Y values.
void holtzer(const API::MatrixWorkspace_sptr &ws)
Performs the Holtzer transformation: IQ v Q.
void init() override
Virtual method - must be overridden by concrete algorithm.
Definition: IQTransform.cpp:59
TransformMap m_transforms
A map of transformation name and function pointers.
Definition: IQTransform.h:74
std::shared_ptr< Kernel::Units::Label > m_label
Definition: IQTransform.h:76
void debyeBueche(const API::MatrixWorkspace_sptr &ws)
Performs the Debye-Bueche transformation: 1/sqrt(I) v Q^2 The output is set to zero for negative inpu...
void(IQTransform::*)(const API::MatrixWorkspace_sptr &) TransformFunc
Definition: IQTransform.h:72
Support for a property that holds an array of values.
Definition: ArrayProperty.h:28
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 error(const std::string &msg)
Logs at error level.
Definition: Logger.cpp:77
void warning(const std::string &msg)
Logs at warning level.
Definition: Logger.cpp:86
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
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
Functor giving the product of the squares of the arguments.
Definition: VectorHelper.h:177