17#include "MantidHistogramData/Histogram.h"
32using namespace Kernel;
34using namespace DataObjects;
41 "The input workspace must either have constant width bins or is a "
43 "workspace. It is also assumed that all spectra have the same X bin "
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 "
53 "The X value at which to end the background fit. Mandatory "
54 "for the Linear Fit and Mean modes, ignored by Moving "
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)");
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 "
71 "the InputWorkspace and returned or just returned (default: "
72 "Subtract Background)");
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.",
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 "
88 "The width of the moving average window in bins. Mandatory "
89 "for the Moving Average mode.");
91 std::make_unique<EnabledWhenProperty>(
"Mode",
IS_EQUAL_TO,
"Moving Average"));
99 const size_t numHists = inputWS->getNumberHistograms();
100 const size_t blocksize = inputWS->blocksize();
103 throw std::runtime_error(
"CalculateFlatBackground with only one bin is "
106 const bool skipMonitors =
getProperty(
"SkipMonitors");
107 const bool nullifyNegative =
getProperty(
"NullifyNegativeValues");
108 const std::string modeString =
getProperty(
"Mode");
110 if (modeString ==
"Mean") {
112 }
else if (modeString ==
"Moving Average") {
121 throw std::runtime_error(
"StartX property not set to any value");
124 throw std::runtime_error(
"EndX property not set to any value");
131 throw std::runtime_error(
"AveragingWindowWidth property not set to any value");
134 if (windowWidth <= 0) {
135 throw std::runtime_error(
"AveragingWindowWidth zero or negative");
137 if (blocksize <
static_cast<size_t>(windowWidth)) {
138 throw std::runtime_error(
"AveragingWindowWidth is larger than the number "
139 "of bins in InputWorkspace");
144 std::vector<int> wsInds =
getProperty(
"WorkspaceIndexList");
146 if (wsInds.empty()) {
147 wsInds.resize(numHists);
148 std::iota(wsInds.begin(), wsInds.end(), 0);
152 const bool removeBackground = std::string(
getProperty(
"outputMode")) ==
"Subtract Background";
155 m_progress = std::make_unique<Progress>(
this, 0.0, 0.2, numHists);
160 if (outputWS != inputWS) {
161 outputWS = create<MatrixWorkspace>(*inputWS);
165 double backgroundTotal = 0;
166 size_t calculationCount = 0;
169 for (int64_t i = 0; i < static_cast<int64_t>(numHists); ++i) {
173 auto histogram = inputWS->histogram(i);
174 bool wasCounts =
false;
175 if (histogram.yMode() == HistogramData::Histogram::YMode::Counts) {
177 histogram.convertToFrequencies();
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)) {
187 " Processing background anyway\n");
188 }
else if (spectrumInfo.isMonitor(i)) {
189 skipCalculation =
true;
192 double background = 0;
194 if (!skipCalculation) {
199 LinearFit(histogram, background, variance, startX, endX);
202 Mean(histogram, background, variance, startX, endX);
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 "
214 if (removeBackground) {
215 g_log.
warning() <<
" The spectrum has been left unchanged.\n";
217 g_log.
warning() <<
" The output background has been set to zero.\n";
220 backgroundTotal += background;
222 HistogramData::Histogram outHistogram(histogram);
223 auto &ys = outHistogram.mutableY();
224 auto &es = outHistogram.mutableE();
225 if (removeBackground) {
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)) {
236 err = es[j] > background ? es[j] : background;
243 for (
size_t j = 0; j < ys.size(); ++j) {
244 const double originalVal = histogram.y()[j];
245 if (background < 0) {
248 }
else if (nullifyNegative && (background > originalVal)) {
250 es[j] = es[j] > background ? es[j] : background;
253 es[j] = std::sqrt(variance);
258 outHistogram.convertToCounts();
260 outputWS->setHistogram(i, outHistogram);
266 g_log.
debug() << calculationCount <<
" spectra corrected\n";
267 g_log.
information() <<
"The mean background was " << backgroundTotal /
static_cast<double>(calculationCount) <<
".\n";
285 const std::string failure(
"XMax must be greater than XMin.");
287 throw std::invalid_argument(failure);
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();
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 "
323 ptrdiff_t startInd = std::lower_bound(XS.begin(), XS.end(), startX) - XS.begin() - 1;
324 if (startInd == -1) {
330 const ptrdiff_t endInd = std::lower_bound(XS.begin() + startInd, XS.end(), endX) - XS.begin() - 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");
339 const auto numBins =
static_cast<double>(1 + endInd - startInd);
342 background = std::accumulate(YS.begin() + startInd, YS.begin() + endInd + 1, 0.0) / numBins;
361 const double startX,
const double endX) {
363 WS->setHistogram(0, histogram);
370 childAlg->setProperty<
bool>(
"CreateOutput",
true);
371 childAlg->setProperty<
int>(
"WorkspaceIndex", 0);
372 childAlg->setProperty<
double>(
"StartX", startX);
373 childAlg->setProperty<
double>(
"EndX", endX);
376 childAlg->setProperty<std::string>(
"Minimizer",
"Levenberg-MarquardtMD");
378 childAlg->executeAsChildAlg();
380 std::string outputStatus = childAlg->getProperty(
"OutputStatus");
381 if (outputStatus !=
"success") {
382 g_log.
warning(
"Unable to successfully fit the data: " + outputStatus);
391 output->find(
static_cast<std::string
>(
"A0"), rowA0, 0);
392 output->find(
static_cast<std::string
>(
"A1"), rowA1, 0);
395 const double intercept = output->cell<
double>(rowA0, 1);
396 const double slope = output->cell<
double>(rowA1, 1);
398 const double centre = (startX + endX) / 2.0;
402 background = slope * centre + intercept;
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;
422 for (
size_t i = 0; i < ys.size(); ++i) {
425 for (
size_t j = 0; j < windowWidth; ++j) {
426 size_t index = i + j;
427 if (
index >= ys.size()) {
434 const double average = sum /
static_cast<double>(windowWidth);
435 if (average < currentMin) {
436 currentMin = average;
437 currentVariance = varSqSum /
static_cast<double>(windowWidth * windowWidth);
440 background = currentMin;
441 variance = currentVariance;
#define DECLARE_ALGORITHM(classname)
std::map< DeltaEMode::Type, std::string > index
#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...
#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.
Kernel::Property * getPointerToProperty(const std::string &name) const override
Get a property by name.
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
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.
bool isDefault(const std::string &name) const
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 init() override
Initialisation code.
void checkRange(double &startX, double &endX)
Checks that the range parameters have been set correctly.
void exec() override
Execution code.
Support for a property that holds an array of values.
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.
void error(const std::string &msg)
Logs at error level.
void warning(const std::string &msg)
Logs at warning level.
void information(const std::string &msg)
Logs at information level.
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
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.
constexpr int EMPTY_INT() noexcept
Returns what we consider an "empty" integer within a property.
constexpr double EMPTY_DBL() noexcept
Returns what we consider an "empty" double within a property.
std::string to_string(const wide_integer< Bits, Signed > &n)
@ Input
An input workspace.
@ Output
An output workspace.
Functor to accumulate a sum of squares.