diff --git a/Framework/Algorithms/src/CalculatePolynomialBackground.cpp b/Framework/Algorithms/src/CalculatePolynomialBackground.cpp
index 472ef3c87cdd7312a1e81eb92319bda955733f17..5d4691b252d205a6d0ca8dd1cfbd4447e04501b1 100644
--- a/Framework/Algorithms/src/CalculatePolynomialBackground.cpp
+++ b/Framework/Algorithms/src/CalculatePolynomialBackground.cpp
@@ -9,18 +9,26 @@
#include "MantidKernel/ArrayOrderedPairsValidator.h"
#include "MantidKernel/ArrayProperty.h"
#include "MantidKernel/BoundedValidator.h"
+#include "MantidKernel/ListValidator.h"
#include <utility>
namespace {
/// String constants for algorithm's properties.
namespace Prop {
+const static std::string COST_FUNCTION = "CostFunction";
const static std::string INPUT_WS = "InputWorkspace";
const static std::string OUTPUT_WS = "OutputWorkspace";
const static std::string POLY_DEGREE = "Degree";
const static std::string XRANGES = "XRanges";
}
+/// String constants for cost function options.
+namespace CostFunc {
+const static std::string UNWEIGHTED_LEAST_SQUARES = "Unweighted least squares";
+const static std::string WEIGHTED_LEAST_SQUARES = "Least squares";
+}
+
/** Filters ranges completely outside the histogram X values.
* @param ranges a vector of start-end pairs to filter
* @param ws a workspace
@@ -155,11 +163,10 @@ std::vector<double> invertRanges(const std::vector<double> &ranges) {
* @param ranges a vector defining the fitting intervals
* @return a vector of final fitted parameters
*/
-std::vector<double> executeFit(Mantid::API::Algorithm &fit,
- const std::string &function,
- Mantid::API::MatrixWorkspace_sptr &ws,
- const size_t wsIndex,
- const std::vector<double> &ranges) {
+std::vector<double>
+executeFit(Mantid::API::Algorithm &fit, const std::string &function,
+ Mantid::API::MatrixWorkspace_sptr &ws, const size_t wsIndex,
+ const std::vector<double> &ranges, const std::string &costFunction) {
const auto fitRanges = histogramRanges(ranges, *ws, wsIndex);
const auto excludedRanges = invertRanges(fitRanges);
fit.setProperty("Function", function);
@@ -168,6 +175,8 @@ std::vector<double> executeFit(Mantid::API::Algorithm &fit,
fit.setProperty("StartX", fitRanges.front());
fit.setProperty("EndX", fitRanges.back());
fit.setProperty("Exclude", excludedRanges);
+ fit.setProperty("Minimizer", "Levenberg-MarquardtMD");
+ fit.setProperty(Prop::COST_FUNCTION, costFunction);
fit.setProperty("CreateOutput", true);
fit.executeAsChildAlg();
Mantid::API::ITableWorkspace_sptr fitResult =
@@ -268,6 +277,12 @@ void CalculatePolynomialBackground::init() {
declareProperty(Kernel::make_unique<Kernel::ArrayProperty<double>>(
Prop::XRANGES, std::vector<double>(), orderedPairs),
"A list of fitting ranges given as pairs of X values.");
+ std::array<std::string, 2> costFuncOpts{
+ {CostFunc::WEIGHTED_LEAST_SQUARES, CostFunc::UNWEIGHTED_LEAST_SQUARES}};
+ declareProperty(
+ Prop::COST_FUNCTION, CostFunc::WEIGHTED_LEAST_SQUARES.c_str(),
+ boost::make_shared<Kernel::ListValidator<std::string>>(costFuncOpts),
+ "The cost function to be passed to the Fit algorithm.");
}
//----------------------------------------------------------------------------------------------
@@ -279,6 +294,7 @@ void CalculatePolynomialBackground::exec() {
API::MatrixWorkspace_sptr outWS{
DataObjects::create<DataObjects::Workspace2D>(*inWS)};
const std::vector<double> inputRanges = getProperty(Prop::XRANGES);
+ const std::string costFunction = getProperty(Prop::COST_FUNCTION);
const auto polyDegree =
static_cast<size_t>(static_cast<int>(getProperty(Prop::POLY_DEGREE)));
const std::vector<double> initialParams(polyDegree + 1, 0.1);
@@ -290,7 +306,8 @@ void CalculatePolynomialBackground::exec() {
PARALLEL_START_INTERUPT_REGION
const bool logging{false};
auto fit = createChildAlgorithm("Fit", 0, 0, logging);
- const auto parameters = executeFit(*fit, fitFunction, inWS, i, inputRanges);
+ const auto parameters =
+ executeFit(*fit, fitFunction, inWS, i, inputRanges, costFunction);
const auto bkgFunction = makeFunctionString(parameters);
evaluateInPlace(bkgFunction, *outWS, i);
progress.report();
diff --git a/Framework/Algorithms/test/CalculatePolynomialBackgroundTest.h b/Framework/Algorithms/test/CalculatePolynomialBackgroundTest.h
index 8106f6d6b51eeacf7d9e516f735049abef51639a..19122f111200a432ca1e11103474c94512bf92a5 100644
--- a/Framework/Algorithms/test/CalculatePolynomialBackgroundTest.h
+++ b/Framework/Algorithms/test/CalculatePolynomialBackgroundTest.h
@@ -98,7 +98,73 @@ public:
const auto &bkgEs = outWS->e(histI);
const auto &bkgXs = outWS->x(histI);
for (size_t binI = 0; binI < nBin; ++binI) {
- TS_ASSERT_DELTA(bkgYs[binI], ys[binI], 1e-12)
+ TS_ASSERT_DELTA(bkgYs[binI], ys[binI], 1e-10)
+ TS_ASSERT_EQUALS(bkgEs[binI], 0)
+ TS_ASSERT_EQUALS(bkgXs[binI], xs[binI])
+ }
+ }
+ }
+
+ void test_costFuctionLeastSquares() {
+ using namespace WorkspaceCreationHelper;
+ const size_t nHist{2};
+ const HistogramData::Counts counts{0, 4, 0, 0};
+ const HistogramData::CountStandardDeviations stdDevs{0, 0.001, 0, 0};
+ const HistogramData::BinEdges edges{0, 1, 2, 3, 4};
+ API::MatrixWorkspace_sptr ws(
+ DataObjects::create<DataObjects::Workspace2D>(
+ nHist, HistogramData::Histogram(edges, counts, stdDevs)).release());
+ auto alg = makeAlgorithm();
+ TS_ASSERT_THROWS_NOTHING(alg->setProperty("InputWorkspace", ws))
+ TS_ASSERT_THROWS_NOTHING(alg->setProperty("OutputWorkspace", "outputWS"))
+ TS_ASSERT_THROWS_NOTHING(alg->setProperty("Degree", 0))
+ TS_ASSERT_THROWS_NOTHING(alg->setProperty("CostFunction", "Least squares"))
+ TS_ASSERT_THROWS_NOTHING(alg->execute())
+ TS_ASSERT(alg->isExecuted())
+ API::MatrixWorkspace_sptr outWS = alg->getProperty("OutputWorkspace");
+ TS_ASSERT(outWS)
+ for (size_t histI = 0; histI < nHist; ++histI) {
+ const auto &xs = ws->x(histI);
+ const auto &bkgYs = outWS->y(histI);
+ const auto &bkgEs = outWS->e(histI);
+ const auto &bkgXs = outWS->x(histI);
+ for (size_t binI = 0; binI < counts.size(); ++binI) {
+ // Number 4 in counts is heavily weighted by the small error.
+ TS_ASSERT_DELTA(bkgYs[binI], 4, 1e-5)
+ TS_ASSERT_EQUALS(bkgEs[binI], 0)
+ TS_ASSERT_EQUALS(bkgXs[binI], xs[binI])
+ }
+ }
+ }
+
+ void test_costFuctionUnweightedLeastSquares() {
+ using namespace WorkspaceCreationHelper;
+ const size_t nHist{2};
+ const HistogramData::Counts counts{0, 4, 0, 0};
+ const HistogramData::BinEdges edges{0, 1, 2, 3, 4};
+ API::MatrixWorkspace_sptr ws(
+ DataObjects::create<DataObjects::Workspace2D>(
+ nHist, HistogramData::Histogram(edges, counts)).release());
+ auto alg = makeAlgorithm();
+ TS_ASSERT_THROWS_NOTHING(alg->setProperty("InputWorkspace", ws))
+ TS_ASSERT_THROWS_NOTHING(alg->setProperty("OutputWorkspace", "outputWS"))
+ TS_ASSERT_THROWS_NOTHING(alg->setProperty("Degree", 0))
+ TS_ASSERT_THROWS_NOTHING(
+ alg->setProperty("CostFunction", "Unweighted least squares"))
+ TS_ASSERT_THROWS_NOTHING(alg->execute())
+ TS_ASSERT(alg->isExecuted())
+ API::MatrixWorkspace_sptr outWS = alg->getProperty("OutputWorkspace");
+ TS_ASSERT(outWS)
+ // Unweighted fitting of zeroth order polynomial is equivalent to the mean.
+ const double result = std::accumulate(counts.cbegin(), counts.cend(), 0.0) /
+ static_cast<double>(counts.size());
+ for (size_t histI = 0; histI < nHist; ++histI) {
+ const auto &xs = ws->x(histI);
+ const auto &bkgYs = outWS->y(histI);
+ const auto &bkgEs = outWS->e(histI);
+ const auto &bkgXs = outWS->x(histI);
+ for (size_t binI = 0; binI < counts.size(); ++binI) {
+ TS_ASSERT_DELTA(bkgYs[binI], result, 1e-5)
TS_ASSERT_EQUALS(bkgEs[binI], 0)
TS_ASSERT_EQUALS(bkgXs[binI], xs[binI])
}
@@ -176,7 +242,7 @@ public:
const auto &bkgXs = outWS->x(0);
const std::vector<double> expected{1.0, 2.0, 3.0, 4.0, 5.0, 6.0};
for (size_t binI = 0; binI < nBin; ++binI) {
- TS_ASSERT_DELTA(bkgYs[binI], expected[binI], 1e-12)
+ TS_ASSERT_DELTA(bkgYs[binI], expected[binI], 1e-10)
TS_ASSERT_EQUALS(bkgEs[binI], 0)
TS_ASSERT_EQUALS(bkgXs[binI], xs[binI])
}
diff --git a/docs/source/algorithms/CalculatePolynomialBackground-v1.rst b/docs/source/algorithms/CalculatePolynomialBackground-v1.rst
index 56d052de6697a9bb1b02acd7872218dd24ba0f13..4887185a85b59e06f4367d95b6008ed308fefab5 100644
--- a/docs/source/algorithms/CalculatePolynomialBackground-v1.rst
+++ b/docs/source/algorithms/CalculatePolynomialBackground-v1.rst
@@ -12,6 +12,8 @@ Description
This algorithm calculates backgrounds for the histograms in a workspace by fitting a polynomial to ranges given by *XRanges*. The fitting is done using :ref:`algm-Fit`. The backgrounds are returned as the *OutputWorkspace* which can be subsequently subtracted from the original workspace. The degree of the polynomial can be given by *Degree*. A zeroth degree polynomial would correspond to fitting a constant value. The *XRanges* property is a list of range start and range end value pairs in the X axis units. Overlapping ranges are merged. If no *XRanges* are given, full histograms are used.
+The value of *CostFunction* is passed to :ref:`algm-Fit` as-is. The default option is 'Least squares' which uses the histogram errors as weights. This might not be desirable, e.g. when there are bins with zero counts and zero errors. An 'Unweighted least squares' option is available to deal with such cases.
+
Usage
-----
diff --git a/docs/source/release/v3.12.0/framework.rst b/docs/source/release/v3.12.0/framework.rst
index 6b9e84934a55dcd01114f7af196ae37a089f3bb1..6453e06fc3f81b8fc841b1fb3bc5e10e075ee190 100644
--- a/docs/source/release/v3.12.0/framework.rst
+++ b/docs/source/release/v3.12.0/framework.rst
@@ -23,6 +23,7 @@ Algorithms
----------
:ref:`NormaliseToMonitor <algm-NormaliseToMonitor>` now supports workspaces with detector scans and workspaces with single-count point data.
+- It is now possible to choose between weighted and unweighted fitting in :ref:`CalculatePolynomialBackground <algm-CalculatePolynomialBackground>`.
- :ref:`CreateWorkspace <algm-CreateWorkspace>` will no longer create a default (and potentially wrong) mapping from spectra to detectors, unless a parent workspace is given. This change ensures that accidental bad mappings that could lead to corrupted data are not created silently anymore. This change does *not* affect the use of this algorithm if: (1) a parent workspace is given, or (2) no instrument is loaded into to workspace at a later point, or (3) an instrument is loaded at a later point but ``LoadInstrument`` is used with ``RewriteSpectraMapping=True``. See also the algorithm documentation for details.
- :ref:`ConjoinWorkspaces <algm-ConjoinWorkspaces>` now supports non-constant bins.
- :ref:`Fit <algm-Fit>` will now respect excluded ranges when ``CostFunction = 'Unweighted least squares'``.