Mantid
Loading...
Searching...
No Matches
Divide.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
9using namespace Mantid::API;
10using namespace Mantid::Kernel;
11using namespace Mantid::DataObjects;
12
13namespace Mantid::Algorithms {
14// Register the class into the algorithm factory
16
17void Divide::init() {
19 declareProperty("WarnOnZeroDivide", true,
20 "Algorithm usually warns if "
21 "division by 0 occurs. Set this "
22 "value to false if one does not "
23 "want this message appearing ");
24}
25
27 m_warnOnZeroDivide = getProperty("WarnOnZeroDivide");
29}
30
31void Divide::performBinaryOperation(const HistogramData::Histogram &lhs, const HistogramData::Histogram &rhs,
32 HistogramData::HistogramY &YOut, HistogramData::HistogramE &EOut) {
33 const auto bins = static_cast<int>(lhs.e().size());
34
35 for (int j = 0; j < bins; ++j) {
36 // Get references to the input Y's
37 const double leftY = lhs.y()[j];
38 const double rightY = rhs.y()[j];
39
40 // error dividing two uncorrelated numbers, re-arrange so that you don't
41 // get infinity if leftY==0 (when rightY=0 the Y value and the result will
42 // both be infinity)
43 // (Sa/a)2 + (Sb/b)2 = (Sc/c)2
44 // (Sa c/a)2 + (Sb c/b)2 = (Sc)2
45 // = (Sa 1/b)2 + (Sb (a/b2))2
46 // (Sc)2 = (1/b)2( (Sa)2 + (Sb a/b)2 )
47 EOut[j] = sqrt(pow(lhs.e()[j], 2) + pow(leftY * rhs.e()[j] / rightY, 2)) / fabs(rightY);
48
49 // Copy the result last in case one of the input workspaces is also any
50 // output
51 YOut[j] = leftY / rightY;
52 }
53}
54
55void Divide::performBinaryOperation(const HistogramData::Histogram &lhs, const double rhsY, const double rhsE,
56 HistogramData::HistogramY &YOut, HistogramData::HistogramE &EOut) {
57 if (rhsY == 0 && m_warnOnZeroDivide)
58 g_log.warning() << "Division by zero: the RHS is a single-valued vector "
59 "with value zero."
60 << "\n";
61
62 // Do the right-hand part of the error calculation just once
63 const double rhsFactor = pow(rhsE / rhsY, 2);
64 const auto bins = static_cast<int>(lhs.e().size());
65 for (int j = 0; j < bins; ++j) {
66 // Get reference to input Y
67 const double leftY = lhs.y()[j];
68
69 // see comment in the function above for the error formula
70 EOut[j] = sqrt(pow(lhs.e()[j], 2) + pow(leftY, 2) * rhsFactor) / fabs(rhsY);
71 // Copy the result last in case one of the input workspaces is also any
72 // output
73 YOut[j] = leftY / rhsY;
74 }
75}
76
79 if (rhs->YUnit().empty() || !WorkspaceHelpers::matchingBins(*lhs, *rhs, true)) {
80 // Do nothing
81 }
82 // If the Y units match, then the output will be a distribution and will be
83 // dimensionless
84 else if (lhs->YUnit() == rhs->YUnit() && m_rhsBlocksize > 1) {
85 out->setYUnit("");
86 out->setDistribution(true);
87 }
88 // Else we need to set the unit that results from the division
89 else {
90 if (!lhs->YUnit().empty())
91 out->setYUnit(lhs->YUnit() + "/" + rhs->YUnit());
92 else
93 out->setYUnit("1/" + rhs->YUnit());
94 }
95}
96
97// ===================================== EVENT LIST BINARY OPERATIONS
98// ==========================================
106 // We must histogram the rhs event list to divide.
107 MantidVec rhsY, rhsE;
108 rhs.generateHistogram(rhs.readX(), rhsY, rhsE);
109 lhs.divide(rhs.readX(), rhsY, rhsE);
110}
111
121 const MantidVec &rhsE) {
122 // Divide is implemented at the EventList level.
123 lhs.divide(rhsX, rhsY, rhsE);
124}
125
134void Divide::performEventBinaryOperation(DataObjects::EventList &lhs, const double &rhsY, const double &rhsE) {
135 // Multiply is implemented at the EventList level.
136 lhs.divide(rhsY, rhsE);
137}
138
139//---------------------------------------------------------------------------------------------
147 if (m_elhs) {
148 // The lhs workspace is an EventWorkspace. It can be divided while keeping
149 // event-ishness
150 // Output will be EW
152 // Histogram sizes need not match
153 m_matchXSize = false;
154 } else {
155 m_keepEventWorkspace = false;
156 m_matchXSize = true;
157 }
158
159 // Division is not commutative = you can't flip sides.
160 m_flipSides = false;
161 // The RHS operand will be histogrammed first.
163}
164
165//--------------------------------------------------------------------------------------------
179 // --- Check for event workspaces - different than workspaces 2D! ---
180
181 // A SingleValueWorkspace on the right matches anything
182 if (rhs->size() == 1)
183 return "";
184
185 // A SingleValueWorkspace on the left only matches if rhs was single value
186 // too. Why are you using mantid to do simple math?!?
187 if (lhs->size() == 1)
188 return "The left side cannot contain a single value if the right side "
189 "isn't also a single value.";
190
191 // If RHS only has one value (1D vertical), the number of histograms needs to
192 // match.
193 // Each lhs spectrum will be divided by that scalar
194 // Are we allowing the division by different # of spectra, using detector IDs
195 // to match up?
197 (m_rhsBlocksize == 1 && lhs->getNumberHistograms() == rhs->getNumberHistograms())) {
198 return "";
199 }
200
201 if (m_matchXSize) {
202 // Past this point, for a 2D WS operation, we require the X arrays to match.
203 // Note this only checks the first spectrum
204 if (!WorkspaceHelpers::matchingBins(*lhs, *rhs, true)) {
205 return "X arrays must match when dividing 2D workspaces.";
206 }
207 }
208
209 // We don't need to check for matching bins for events. Yay events!
210 const size_t rhsSpec = rhs->getNumberHistograms();
211
212 // If the rhs has a single spectrum, then we can divide. The block size does
213 // NOT need to match,
214 if (rhsSpec == 1)
215 return "";
216
217 // Otherwise, the number of histograms needs to match, but the block size of
218 // each does NOT need to match.
219
220 if (lhs->getNumberHistograms() == rhs->getNumberHistograms()) {
221 return "";
222 } else {
223 return "Number of histograms not identical.";
224 }
225}
226
227} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
const std::vector< double > & rhs
#define fabs(x)
Definition: Matrix.cpp:22
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
size_t m_rhsBlocksize
Cache for RHS workspace's blocksize.
bool m_keepEventWorkspace
Variable set to true if the operation allows the output to stay as an EventWorkspace.
void init() override
Initialisation method.
bool m_AllowDifferentNumberSpectra
The property value.
bool m_useHistogramForRhsEventWorkspace
Are we going to use the histogram representation of the RHS event list when performing the operation?...
bool m_matchXSize
matchXSize set to true if the X sizes of histograms must match.
bool m_flipSides
flipSides set to true if the rhs and lhs operands should be flipped - for commutative binary operatio...
void exec() override
Executes the algorithm.
DataObjects::EventWorkspace_const_sptr m_elhs
Left-hand side EventWorkspace.
Divide performs the division of two input workspaces.
Definition: Divide.h:29
std::string checkSizeCompatibility(const API::MatrixWorkspace_const_sptr lhs, const API::MatrixWorkspace_const_sptr rhs) const override
Performs a simple check to see if the sizes of two workspaces are compatible for a binary operation I...
Definition: Divide.cpp:177
void exec() override
Executes the algorithm.
Definition: Divide.cpp:26
void performEventBinaryOperation(DataObjects::EventList &lhs, const DataObjects::EventList &rhs) override
Carries out the binary operation IN-PLACE on a single EventList, with another EventList as the right-...
Definition: Divide.cpp:105
void setOutputUnits(const API::MatrixWorkspace_const_sptr lhs, const API::MatrixWorkspace_const_sptr rhs, API::MatrixWorkspace_sptr out) override
Should be overridden by operations that need to manipulate the units of the output workspace.
Definition: Divide.cpp:77
void performBinaryOperation(const HistogramData::Histogram &lhs, const HistogramData::Histogram &rhs, HistogramData::HistogramY &YOut, HistogramData::HistogramE &EOut) override
Carries out the binary operation on a single spectrum, with another spectrum as the right-hand operan...
Definition: Divide.cpp:31
void checkRequirements() override
Check what operation will be needed in order to apply the operation to these two types of workspaces.
Definition: Divide.cpp:146
A class for holding :
Definition: EventList.h:56
void divide(const double value, const double error=0.0) override
Divide the weights in this event list by a scalar with an (optional) error.
Definition: EventList.cpp:3633
void warning(const std::string &msg)
Logs at warning level.
Definition: Logger.cpp:86
std::shared_ptr< const MatrixWorkspace > MatrixWorkspace_const_sptr
shared pointer to the matrix workspace base class (const version)
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
static bool matchingBins(const MatrixWorkspace &ws1, const MatrixWorkspace &ws2, const bool firstOnly=false)
Checks whether the bins (X values) of two workspace are the same.