Mantid
Loading...
Searching...
No Matches
CalculateFlatBackground.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 +
17#include "MantidHistogramData/Histogram.h"
22
23#include <algorithm>
24#include <climits>
25#include <numeric>
26
27namespace Mantid::Algorithms {
28
29// Register the algorithm into the AlgorithmFactory
30DECLARE_ALGORITHM(CalculateFlatBackground)
31
32using namespace Kernel;
33using namespace API;
34using namespace DataObjects;
35
38
40 declareProperty(std::make_unique<WorkspaceProperty<>>("InputWorkspace", "", Direction::Input),
41 "The input workspace must either have constant width bins or is a "
42 "distribution\n"
43 "workspace. It is also assumed that all spectra have the same X bin "
44 "boundaries");
45 declareProperty(std::make_unique<WorkspaceProperty<>>("OutputWorkspace", "", Direction::Output),
46 "Name to use for the output workspace.");
48 "The X value at which to start the background fit. Mandatory "
49 "for the Linear Fit and Mean modes, ignored by Moving "
50 "Average.");
51 setPropertySettings("StartX", std::make_unique<EnabledWhenProperty>("Mode", IS_NOT_EQUAL_TO, "Moving Average"));
53 "The X value at which to end the background fit. Mandatory "
54 "for the Linear Fit and Mean modes, ignored by Moving "
55 "Average.");
56 setPropertySettings("EndX", std::make_unique<EnabledWhenProperty>("Mode", IS_NOT_EQUAL_TO, "Moving Average"));
57 declareProperty(std::make_unique<ArrayProperty<int>>("WorkspaceIndexList"),
58 "Indices of the spectra that will have their background removed\n"
59 "default: modify all spectra");
60 std::vector<std::string> modeOptions{"Linear Fit", "Mean", "Moving Average"};
61 declareProperty("Mode", "Linear Fit", std::make_shared<StringListValidator>(modeOptions),
62 "The background count rate is estimated either by taking a "
63 "mean, doing a linear fit, or taking the\n"
64 "minimum of a moving average (default: Linear Fit)");
65 // Property to determine whether we subtract the background or just return the
66 // background.
67 std::vector<std::string> outputOptions{"Subtract Background", "Return Background"};
68 declareProperty("OutputMode", "Subtract Background", std::make_shared<StringListValidator>(outputOptions),
69 "Once the background has been determined it can either be "
70 "subtracted from \n"
71 "the InputWorkspace and returned or just returned (default: "
72 "Subtract Background)");
73 declareProperty("SkipMonitors", false,
74 "By default, the algorithm calculates and removes background "
75 "from monitors in the same way as from normal detectors\n"
76 "If this property is set to true, background is not "
77 "calculated/removed from monitors.",
79 declareProperty("NullifyNegativeValues", true,
80 "When background is subtracted, signals in some time "
81 "channels may become negative.\n"
82 "If this option is true, signal in such bins is nullified "
83 "and the module of the removed signal"
84 "is added to the error. If false, the signal and errors are "
85 "left unchanged",
87 declareProperty("AveragingWindowWidth", Mantid::EMPTY_INT(),
88 "The width of the moving average window in bins. Mandatory "
89 "for the Moving Average mode.");
90 setPropertySettings("AveragingWindowWidth",
91 std::make_unique<EnabledWhenProperty>("Mode", IS_EQUAL_TO, "Moving Average"));
92}
93
95 // Retrieve the input workspace
96 API::MatrixWorkspace_const_sptr inputWS = getProperty("InputWorkspace");
97
98 // Copy over all the data
99 const size_t numHists = inputWS->getNumberHistograms();
100 const size_t blocksize = inputWS->blocksize();
101
102 if (blocksize == 1)
103 throw std::runtime_error("CalculateFlatBackground with only one bin is "
104 "not possible.");
105
106 const bool skipMonitors = getProperty("SkipMonitors");
107 const bool nullifyNegative = getProperty("NullifyNegativeValues");
108 const std::string modeString = getProperty("Mode");
110 if (modeString == "Mean") {
111 mode = Modes::MEAN;
112 } else if (modeString == "Moving Average") {
114 }
115 double startX, endX;
116 int windowWidth = 0;
117 switch (mode) {
119 case Modes::MEAN:
120 if (getPointerToProperty("StartX")->isDefault()) {
121 throw std::runtime_error("StartX property not set to any value");
122 }
123 if (getPointerToProperty("EndX")->isDefault()) {
124 throw std::runtime_error("EndX property not set to any value");
125 }
126 // Get the required X range
127 checkRange(startX, endX);
128 break;
130 if (getPointerToProperty("AveragingWindowWidth")->isDefault()) {
131 throw std::runtime_error("AveragingWindowWidth property not set to any value");
132 }
133 windowWidth = getProperty("AveragingWindowWidth");
134 if (windowWidth <= 0) {
135 throw std::runtime_error("AveragingWindowWidth zero or negative");
136 }
137 if (blocksize < static_cast<size_t>(windowWidth)) {
138 throw std::runtime_error("AveragingWindowWidth is larger than the number "
139 "of bins in InputWorkspace");
140 }
141 break;
142 }
143
144 std::vector<int> wsInds = getProperty("WorkspaceIndexList");
145 // check if the user passed an empty list, if so all of spec will be processed
146 if (wsInds.empty()) {
147 wsInds.resize(numHists);
148 std::iota(wsInds.begin(), wsInds.end(), 0);
149 }
150
151 // Are we removing the background?
152 const bool removeBackground = std::string(getProperty("outputMode")) == "Subtract Background";
153
154 // Initialize the progress reporting object
155 m_progress = std::make_unique<Progress>(this, 0.0, 0.2, numHists);
156
157 MatrixWorkspace_sptr outputWS = getProperty("OutputWorkspace");
158 // If input and output workspaces are not the same, create a new workspace for
159 // the output
160 if (outputWS != inputWS) {
161 outputWS = create<MatrixWorkspace>(*inputWS);
162 }
163
164 // For logging purposes.
165 double backgroundTotal = 0;
166 size_t calculationCount = 0;
167
168 PARALLEL_FOR_IF(Kernel::threadSafe(*inputWS, *outputWS))
169 for (int64_t i = 0; i < static_cast<int64_t>(numHists); ++i) {
171 // Extract a histogram from inputWS, work on it and, finally, put it to
172 // outputWS.
173 auto histogram = inputWS->histogram(i);
174 bool wasCounts = false;
175 if (histogram.yMode() == HistogramData::Histogram::YMode::Counts) {
176 wasCounts = true;
177 histogram.convertToFrequencies();
178 }
179 bool skipCalculation = std::find(wsInds.cbegin(), wsInds.cend(), i) == wsInds.cend();
180 if (!skipCalculation && skipMonitors) {
181 const auto &spectrumInfo = inputWS->spectrumInfo();
182 if (!spectrumInfo.hasDetectors(i)) {
183 // Do nothing.
184 // not every spectra is the monitor or detector, some spectra have no
185 // instrument components attached.
186 g_log.information(" Can not find detector for spectra N: " + std::to_string(i) +
187 " Processing background anyway\n");
188 } else if (spectrumInfo.isMonitor(i)) {
189 skipCalculation = true;
190 }
191 }
192 double background = 0;
193 double variance = 0;
194 if (!skipCalculation) {
195 ++calculationCount;
196 // Now call the function the user selected to calculate the background
197 switch (mode) {
199 LinearFit(histogram, background, variance, startX, endX);
200 break;
201 case Modes::MEAN:
202 Mean(histogram, background, variance, startX, endX);
203 break;
205 MovingAverage(histogram, background, variance, windowWidth);
206 break;
207 }
208 }
209 if (background < 0) {
210 g_log.debug() << "The background for spectra index " << i << "was calculated to be " << background << '\n';
211 g_log.warning() << "Problem with calculating the background number of "
212 "counts spectrum with index "
213 << i << ".";
214 if (removeBackground) {
215 g_log.warning() << " The spectrum has been left unchanged.\n";
216 } else {
217 g_log.warning() << " The output background has been set to zero.\n";
218 }
219 } else {
220 backgroundTotal += background;
221 }
222 HistogramData::Histogram outHistogram(histogram);
223 auto &ys = outHistogram.mutableY();
224 auto &es = outHistogram.mutableE();
225 if (removeBackground) {
226 // When subtracting backgrounds, act only if background is positive.
227 if (background >= 0) {
228 for (size_t j = 0; j < ys.size(); ++j) {
229 double val = ys[j] - background;
230 double err = std::sqrt(es[j] * es[j] + variance);
231 if (nullifyNegative && (val < 0)) {
232 val = 0;
233 // The error estimate must go up in this nonideal situation and the
234 // value of background is a good estimate for it. However, don't
235 // reduce the error if it was already more than that
236 err = es[j] > background ? es[j] : background;
237 }
238 ys[j] = val;
239 es[j] = err;
240 }
241 }
242 } else {
243 for (size_t j = 0; j < ys.size(); ++j) {
244 const double originalVal = histogram.y()[j];
245 if (background < 0) {
246 ys[j] = 0;
247 es[j] = 0;
248 } else if (nullifyNegative && (background > originalVal)) {
249 ys[j] = originalVal;
250 es[j] = es[j] > background ? es[j] : background;
251 } else {
252 ys[j] = background;
253 es[j] = std::sqrt(variance);
254 }
255 }
256 }
257 if (wasCounts) {
258 outHistogram.convertToCounts();
259 }
260 outputWS->setHistogram(i, outHistogram);
261 m_progress->report();
263 }
265
266 g_log.debug() << calculationCount << " spectra corrected\n";
267 g_log.information() << "The mean background was " << backgroundTotal / static_cast<double>(calculationCount) << ".\n";
268 // Assign the output workspace to its property
269 setProperty("OutputWorkspace", outputWS);
270}
271
279void CalculateFlatBackground::checkRange(double &startX, double &endX) {
280 // use the overloaded operator =() to get the X value stored in each property
281 startX = getProperty("StartX");
282 endX = getProperty("EndX");
283
284 if (startX > endX) {
285 const std::string failure("XMax must be greater than XMin.");
286 g_log.error(failure);
287 throw std::invalid_argument(failure);
288 }
289}
290
306void CalculateFlatBackground::Mean(const HistogramData::Histogram &histogram, double &background, double &variance,
307 const double startX, const double endX) const {
308 const auto &XS = histogram.x();
309 const auto &YS = histogram.y();
310 const auto &ES = histogram.e();
311 // the function checkRange should already have checked that startX <= endX,
312 // but we still need to check values weren't out side the ranges
313 if ((endX > XS.back()) || (startX < XS.front())) {
314 throw std::out_of_range("Either the property startX or endX is outside the "
315 "range of X-values present in one of the specified "
316 "spectra");
317 }
318 // Get the index of the first bin contains the X-value, which means this is an
319 // inclusive sum. The minus one is because lower_bound() returns index past
320 // the last index pointing to a lower value. For example if startX has a
321 // higher X value than the first bin boundary but lower than the second
322 // lower_bound returns 1, which is the index of the second bin boundary
323 ptrdiff_t startInd = std::lower_bound(XS.begin(), XS.end(), startX) - XS.begin() - 1;
324 if (startInd == -1) { // happens if startX is the first X-value, e.g. the
325 // first X-value is zero and the user selects zero
326 startInd = 0;
327 }
328
329 // the -1 matches definition of startIn, see the comment above that statement
330 const ptrdiff_t endInd = std::lower_bound(XS.begin() + startInd, XS.end(), endX) - XS.begin() - 1;
331 if (endInd == -1) { //
332 throw std::invalid_argument("EndX was set to the start of one of the "
333 "spectra, it must greater than the first "
334 "X-value in any of the specified spectra");
335 }
336
337 // the +1 is because this is an inclusive sum (includes each bin that contains
338 // each X-value). Hence if startInd == endInd we are still analyzing one bin
339 const auto numBins = static_cast<double>(1 + endInd - startInd);
340 // the +1 here is because the accumulate() stops one before the location of
341 // the last iterator
342 background = std::accumulate(YS.begin() + startInd, YS.begin() + endInd + 1, 0.0) / numBins;
343 // The error on the total number of background counts in the background region
344 // is taken as the sqrt the total number counts. To get the the error on the
345 // counts in each bin just divide this by the number of bins. The variance =
346 // error^2 that is the total variance divide by the number of bins _squared_.
347 variance = std::accumulate(ES.begin() + startInd, ES.begin() + endInd + 1, 0.0, VectorHelper::SumSquares<double>()) /
348 (numBins * numBins);
349}
350
360void CalculateFlatBackground::LinearFit(const HistogramData::Histogram &histogram, double &background, double &variance,
361 const double startX, const double endX) {
362 MatrixWorkspace_sptr WS = create<Workspace2D>(1, histogram);
363 WS->setHistogram(0, histogram);
364 auto childAlg = createChildAlgorithm("Fit");
365
366 IFunction_sptr func = API::FunctionFactory::Instance().createFunction("LinearBackground");
367 childAlg->setProperty<IFunction_sptr>("Function", func);
368
369 childAlg->setProperty<MatrixWorkspace_sptr>("InputWorkspace", WS);
370 childAlg->setProperty<bool>("CreateOutput", true);
371 childAlg->setProperty<int>("WorkspaceIndex", 0);
372 childAlg->setProperty<double>("StartX", startX);
373 childAlg->setProperty<double>("EndX", endX);
374 // Default minimizer doesn't work properly even on the easiest cases,
375 // so Levenberg-MarquardtMD is used instead
376 childAlg->setProperty<std::string>("Minimizer", "Levenberg-MarquardtMD");
377
378 childAlg->executeAsChildAlg();
379
380 std::string outputStatus = childAlg->getProperty("OutputStatus");
381 if (outputStatus != "success") {
382 g_log.warning("Unable to successfully fit the data: " + outputStatus);
383 background = -1;
384 return;
385 }
386
387 Mantid::API::ITableWorkspace_sptr output = childAlg->getProperty("OutputParameters");
388
389 // Find rows with parameters we are after
390 size_t rowA0, rowA1;
391 output->find(static_cast<std::string>("A0"), rowA0, 0);
392 output->find(static_cast<std::string>("A1"), rowA1, 0);
393
394 // Linear function is defined as A0 + A1*x
395 const double intercept = output->cell<double>(rowA0, 1);
396 const double slope = output->cell<double>(rowA1, 1);
397
398 const double centre = (startX + endX) / 2.0;
399
400 // Calculate the value of the flat background by taking the value at the
401 // centre point of the fit
402 background = slope * centre + intercept;
403 // ATM we don't calculate the error here.
404 variance = 0;
405}
406
415void CalculateFlatBackground::MovingAverage(const HistogramData::Histogram &histogram, double &background,
416 double &variance, const size_t windowWidth) const {
417 const auto &ys = histogram.y();
418 const auto &es = histogram.e();
419 double currentMin = std::numeric_limits<double>::max();
420 double currentVariance = 0;
421
422 for (size_t i = 0; i < ys.size(); ++i) {
423 double sum = 0;
424 double varSqSum = 0;
425 for (size_t j = 0; j < windowWidth; ++j) {
426 size_t index = i + j;
427 if (index >= ys.size()) {
428 // Cyclic boundary conditions.
429 index -= ys.size();
430 }
431 sum += ys[index];
432 varSqSum += es[index] * es[index];
433 }
434 const double average = sum / static_cast<double>(windowWidth);
435 if (average < currentMin) {
436 currentMin = average;
437 currentVariance = varSqSum / static_cast<double>(windowWidth * windowWidth);
438 }
439 }
440 background = currentMin;
441 variance = currentVariance;
442}
443
444} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
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
Kernel::Property * getPointerToProperty(const std::string &name) const override
Get a property by name.
Definition: Algorithm.cpp:2033
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
bool isDefault(const std::string &name) const
Definition: Algorithm.cpp:2084
A property class for workspaces.
void MovingAverage(const HistogramData::Histogram &histogram, double &background, double &variance, const size_t windowWidth) const
Utilizes cyclic boundary conditions when calculating the average in the window.
void LinearFit(const HistogramData::Histogram &histogram, double &background, double &variance, const double startX, const double endX)
Uses linear algorithm to do the fitting.
std::unique_ptr< API::Progress > m_progress
Progress reporting.
void Mean(const HistogramData::Histogram &histogram, double &background, double &variance, const double startX, const double endX) const
Gets the mean number of counts in each bin the background region and the variance (error^2) of that n...
void checkRange(double &startX, double &endX)
Checks that the range parameters have been set correctly.
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 setPropertySettings(const std::string &name, std::unique_ptr< IPropertySettings > settings)
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
void information(const std::string &msg)
Logs at information level.
Definition: Logger.cpp:105
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
std::shared_ptr< ITableWorkspace > ITableWorkspace_sptr
shared pointer to Mantid::API::ITableWorkspace
std::shared_ptr< const MatrixWorkspace > MatrixWorkspace_const_sptr
shared pointer to the matrix workspace base class (const version)
std::shared_ptr< IFunction > IFunction_sptr
shared pointer to the function base class
Definition: IFunction.h:732
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
Modes
Enumeration for the different operating modes.
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
constexpr int EMPTY_INT() noexcept
Returns what we consider an "empty" integer within a property.
Definition: EmptyValues.h:25
constexpr double EMPTY_DBL() noexcept
Returns what we consider an "empty" double within a property.
Definition: EmptyValues.h:43
std::string to_string(const wide_integer< Bits, Signed > &n)
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54
Functor to accumulate a sum of squares.
Definition: VectorHelper.h:170