diff --git a/Framework/Algorithms/inc/MantidAlgorithms/CalculatePolynomialBackground.h b/Framework/Algorithms/inc/MantidAlgorithms/CalculatePolynomialBackground.h
index 93e04fa8f4011a953258365166949ef086e2cdd6..fa1c007f4e81574db5012560ba04fe5a406e7ac7 100644
--- a/Framework/Algorithms/inc/MantidAlgorithms/CalculatePolynomialBackground.h
+++ b/Framework/Algorithms/inc/MantidAlgorithms/CalculatePolynomialBackground.h
@@ -40,8 +40,6 @@ public:
private:
void init() override;
void exec() override;
- std::vector<double> includedRanges(const std::pair<double, double> &totalRange) const;
- std::pair<double, double> totalRange(API::MatrixWorkspace &ws, const size_t wsIndex) const;
};
} // namespace Algorithms
diff --git a/Framework/Algorithms/src/CalculatePolynomialBackground.cpp b/Framework/Algorithms/src/CalculatePolynomialBackground.cpp
index a4e8c07d7f9e276ce7f2cff98f5c273e35035db5..713b468d91dd7ee7762db4fdf72f102f15721d98 100644
--- a/Framework/Algorithms/src/CalculatePolynomialBackground.cpp
+++ b/Framework/Algorithms/src/CalculatePolynomialBackground.cpp
@@ -14,6 +14,7 @@
#include <utility>
namespace {
+/// String constants for algorithm's properties.
namespace Prop {
constexpr char *INPUT_WS = "InputWorkspace";
constexpr char *OUTPUT_WS = "OutputWorkspace";
@@ -21,6 +22,79 @@ constexpr char *POLY_DEGREE = "Degree";
constexpr char *XRANGES = "XRanges";
}
+/** Filters ranges completely outside the histogram X values.
+ * @param ranges a vector of start-end pairs to filter
+ * @param ws a workspace
+ * @param wsIndex a workspace index to specify a histogram in ws
+ * @return a ranges-like vector with filtered pairs removed
+ */
+std::vector<double> filterRangesOutsideX(const std::vector<double> &ranges, const Mantid::API::MatrixWorkspace &ws, const size_t wsIndex) {
+ const auto minX = ws.x(wsIndex).front();
+ const auto maxX = ws.x(wsIndex).back();
+ std::vector<double> filtered;
+ filtered.reserve(ranges.size());
+ for (size_t i = 0; i < ranges.size() / 2; ++i) {
+ const auto first = ranges[2 * i];
+ const auto second = ranges[2 * i + 1];
+ if (!(first > maxX || second < minX)) {
+ filtered.emplace_back(first);
+ filtered.emplace_back(second);
+ }
+ }
+ return filtered;
+}
+
+/** Merges, sorts and limits ranges within totalRange.
+ * @param ranges a vector of start-end pairs to process
+ * @param totalRange a pair of start-end values to limit the output ranges
+ * @return a ranges-like vector of processed ranges
+ */
+std::vector<double> includedRanges(const std::vector<double> &ranges, const std::pair<double, double> &totalRange) {
+ if (ranges.empty()) {
+ return {totalRange.first, totalRange.second};
+ }
+ // Sort the range edges keeping the information whether the edge
+ // 'starts' or 'ends' a range.
+ enum class Edge { start, end };
+ std::vector<std::pair<double, Edge>> edges(ranges.size());
+ for (size_t i = 0; i < ranges.size(); ++i) {
+ edges[i].first = ranges[i];
+ edges[i].second = i % 2 == 0 ? Edge::start : Edge::end;
+ }
+ std::sort(edges.begin(), edges.end(), [](const std::pair<double, Edge> &p1, const std::pair<double, Edge> &p2) {
+ if (p1.first == p2.first)
+ return p1.second == Edge::start;
+ return p1.first < p2.first;
+ });
+ // If an 'end' edge is followed by a 'start', we have a new range. Everything else
+ // can be merged.
+ std::vector<double> mergedRanges;
+ mergedRanges.reserve(ranges.size());
+ auto edgeIt = edges.begin();
+ mergedRanges.emplace_back(std::max(edges.front().first, totalRange.first));
+ while (edgeIt != edges.end()) {
+ auto endEdgeIt = edgeIt + 1;
+ while (endEdgeIt != edges.end()) {
+ const auto val = *endEdgeIt;
+ const auto prevVal = *(endEdgeIt - 1);
+ if (val.second == Edge::start && prevVal.second == Edge::end) {
+ mergedRanges.emplace_back(prevVal.first);
+ mergedRanges.emplace_back(val.first);
+ edgeIt = endEdgeIt;
+ break;
+ }
+ ++endEdgeIt;
+ }
+ ++edgeIt;
+ }
+ mergedRanges.emplace_back(std::min(edges.back().first, totalRange.second));
+ return mergedRanges;
+}
+
+/** Return the gaps between ranges, if any.
+ * @param ranges a vector of start-end pairs to invert
+ * @return a ranges-like vector of gaps between the given ranges.
+ */
std::vector<double> invertRanges(const std::vector<double> &ranges) {
std::vector<double> inversion(ranges.size() - 2);
for (size_t i = 1; i < ranges.size() - 1; ++i) {
@@ -29,6 +103,10 @@ std::vector<double> invertRanges(const std::vector<double> &ranges) {
return inversion;
}
+/** Return a Fit algorithm compatible string representing a polynomial.
+ * @param parameters a vector containing the polynomial coefficients
+ * @return a function string
+ */
std::string makeFunctionString(const std::vector<double> ¶meters) {
const auto degree = parameters.size() - 1;
std::ostringstream s;
@@ -50,6 +128,24 @@ std::string makeFunctionString(const std::vector<double> ¶meters) {
}
return s.str();
}
+
+/** Construct the largest range spanning histogram's X values and ranges.
+ * @param ranges a vector of start-end pairs
+ * @param ws a workspace
+ * @param wsIndex a workspace index identifying a histogram
+ * @return a pair of values spanning a range
+ */
+std::pair<double, double> totalRange(const std::vector<double> &ranges, const Mantid::API::MatrixWorkspace &ws, const size_t wsIndex) {
+ const auto minX = ws.x(wsIndex).front();
+ const auto maxX = ws.x(wsIndex).back();
+ if (ranges.empty()) {
+ return std::pair<double, double>(minX, maxX);
+ }
+ const auto minmaxIt = std::minmax_element(ranges.cbegin(), ranges.cend());
+ const auto minEdge = *minmaxIt.first;
+ const auto maxEdge = *minmaxIt.second;
+ return std::pair<double, double>(std::min(minEdge, minX), std::max(maxEdge, maxX));
+}
}
namespace Mantid {
@@ -103,13 +199,19 @@ void CalculatePolynomialBackground::exec() {
API::MatrixWorkspace_sptr inWS = getProperty(Prop::INPUT_WS);
API::MatrixWorkspace_sptr outWS{DataObjects::create<DataObjects::Workspace2D>(*inWS)};
+ const std::vector<double> ranges = getProperty(Prop::XRANGES);
const auto polyDegree = static_cast<size_t>(static_cast<int>(getProperty(Prop::POLY_DEGREE)));
const std::vector<double> initialParams(polyDegree + 1, 0.1);
const auto fitFunction = makeFunctionString(initialParams);
const auto nHistograms = static_cast<int64_t>(inWS->getNumberHistograms());
const auto nBins = inWS->blocksize();
for (int64_t i = 0; i < nHistograms; ++i) {
- const auto includedR = includedRanges(totalRange(*inWS, i));
+ const auto filteredRanges = filterRangesOutsideX(ranges, *inWS, i);
+ if (filteredRanges.empty()) {
+ throw std::runtime_error("The given XRanges mismatch with the histogram at workspace index " + std::to_string(i));
+ }
+ const auto fullRange = totalRange(filteredRanges, *inWS, i);
+ const auto includedR = includedRanges(filteredRanges, fullRange);
const bool logging{false};
auto fit = createChildAlgorithm("Fit", 0, 0, logging);
fit->setProperty("Function", fitFunction);
@@ -132,6 +234,7 @@ void CalculatePolynomialBackground::exec() {
// We want bkg to directly write to the output workspace.
double *bkgY = const_cast<double *>(outWS->mutableY(i).rawData().data());
bkg->function1D(bkgY, outWS->points(i).rawData().data(), nBins);
+ // Calculate the errors using partial derivatives.
API::BasicJacobian jacobian{nBins, polyDegree + 1};
bkg->functionDeriv1D(&jacobian, outWS->points(i).rawData().data(), nBins);
for (size_t j = 0; j < nBins; ++j) {
@@ -146,60 +249,5 @@ void CalculatePolynomialBackground::exec() {
setProperty(Prop::OUTPUT_WS, outWS);
}
-std::pair<double, double> CalculatePolynomialBackground::totalRange(API::MatrixWorkspace &ws, const size_t wsIndex) const {
- const auto minX = ws.x(wsIndex).front();
- const auto maxX = ws.x(wsIndex).back();
- const std::vector<double> ranges = getProperty(Prop::XRANGES);
- if (ranges.empty()) {
- return std::pair<double, double>(minX, maxX);
- }
- const auto minmaxIt = std::minmax_element(ranges.cbegin(), ranges.cend());
- const auto minEdge = *minmaxIt.first;
- const auto maxEdge = *minmaxIt.second;
- return std::pair<double, double>(std::min(minEdge, minX), std::max(maxEdge, maxX));
-}
-
-std::vector<double> CalculatePolynomialBackground::includedRanges(const std::pair<double, double> &totalRange) const {
- std::vector<double> ranges = getProperty(Prop::XRANGES);
- if (ranges.empty()) {
- return {totalRange.first, totalRange.second};
- }
- // Sort the range edges keeping the information whether the edge
- // 'starts' or 'ends' a range.
- enum class Edge { start, end };
- std::vector<std::pair<double, Edge>> edges(ranges.size());
- for (size_t i = 0; i < ranges.size(); ++i) {
- edges[i].first = ranges[i];
- edges[i].second = i % 2 == 0 ? Edge::start : Edge::end;
- }
- std::sort(edges.begin(), edges.end(), [](const std::pair<double, Edge> &p1, const std::pair<double, Edge> &p2) {
- if (p1.first == p2.first)
- return p1.second == Edge::start;
- return p1.first < p2.first;
- });
- // If an 'end' edge is followed by a 'start', we have a new range. Everything else
- // can be merged.
- std::vector<double> mergedRanges;
- mergedRanges.reserve(ranges.size());
- auto edgeIt = edges.begin();
- mergedRanges.emplace_back(std::max(edges.front().first, totalRange.first));
- while (edgeIt != edges.end()) {
- auto endEdgeIt = edgeIt + 1;
- while (endEdgeIt != edges.end()) {
- const auto val = *endEdgeIt;
- const auto prevVal = *(endEdgeIt - 1);
- if (val.second == Edge::start && prevVal.second == Edge::end) {
- mergedRanges.emplace_back(prevVal.first);
- mergedRanges.emplace_back(val.first);
- edgeIt = endEdgeIt;
- break;
- }
- ++endEdgeIt;
- }
- ++edgeIt;
- }
- mergedRanges.emplace_back(std::min(edges.back().first, totalRange.second));
- return mergedRanges;
-}
} // namespace Algorithms
} // namespace Mantid