From 054013b5c1bf58272b291e15ff3662e7e29b32e8 Mon Sep 17 00:00:00 2001
From: Antti Soininen <soininen@ill.fr>
Date: Tue, 5 Jun 2018 11:27:37 +0200
Subject: [PATCH] Fix floating point inaccuracies in
SparseInstrument::modelHistogram().
MonteCarloAbsorption could fail to interpolate the final wavelengths
thanks to this.
Re #22472
---
.../SampleCorrections/SparseInstrument.h | 5 +--
.../SampleCorrections/SparseInstrument.cpp | 35 ++++++-------------
.../Algorithms/test/SparseInstrumentTest.h | 16 +++++++++
.../release/v3.13.0/direct_inelastic.rst | 1 +
.../release/v3.13.0/indirect_inelastic.rst | 2 ++
5 files changed, 33 insertions(+), 26 deletions(-)
diff --git a/Framework/Algorithms/inc/MantidAlgorithms/SampleCorrections/SparseInstrument.h b/Framework/Algorithms/inc/MantidAlgorithms/SampleCorrections/SparseInstrument.h
index 1ce9cc0de55..f74d7618889 100644
--- a/Framework/Algorithms/inc/MantidAlgorithms/SampleCorrections/SparseInstrument.h
+++ b/Framework/Algorithms/inc/MantidAlgorithms/SampleCorrections/SparseInstrument.h
@@ -55,8 +55,9 @@ geographicalAngles(const Kernel::V3D &p,
const Geometry::ReferenceFrame &refFrame);
DLLExport std::tuple<double, double>
extremeWavelengths(const API::MatrixWorkspace &ws);
-HistogramData::Histogram modelHistogram(const API::MatrixWorkspace &modelWS,
- const size_t wavelengthPoints);
+DLLExport HistogramData::Histogram
+modelHistogram(const API::MatrixWorkspace &modelWS,
+ const size_t wavelengthPoints);
DLLExport API::MatrixWorkspace_uptr
createSparseWS(const API::MatrixWorkspace &modelWS,
const Algorithms::DetectorGridDefinition &grid,
diff --git a/Framework/Algorithms/src/SampleCorrections/SparseInstrument.cpp b/Framework/Algorithms/src/SampleCorrections/SparseInstrument.cpp
index 6c2580cd1e8..3f750c76b35 100644
--- a/Framework/Algorithms/src/SampleCorrections/SparseInstrument.cpp
+++ b/Framework/Algorithms/src/SampleCorrections/SparseInstrument.cpp
@@ -12,6 +12,7 @@
#include "MantidGeometry/Instrument.h"
#include "MantidGeometry/Objects/CSGObject.h"
#include "MantidGeometry/Objects/ShapeFactory.h"
+#include "MantidHistogramData/HistogramIterator.h"
#include "MantidKernel/make_unique.h"
#include <Poco/DOM/AutoPtr.h>
@@ -96,29 +97,12 @@ extremeAngles(const API::MatrixWorkspace &ws) {
std::tuple<double, double> extremeWavelengths(const API::MatrixWorkspace &ws) {
double currentMin = std::numeric_limits<double>::max();
double currentMax = std::numeric_limits<double>::lowest();
- if (ws.histogram(0).xMode() ==
- Mantid::HistogramData::Histogram::XMode::BinEdges) {
- // Deal with point data as in interpolation we need to cover all points,
- // not bin edges.
- for (size_t i = 0; i < ws.getNumberHistograms(); ++i) {
- const auto &xs = ws.x(i);
- const auto x0 = (xs[0] + xs[1]) / 2.0;
- if (x0 < currentMin)
- currentMin = x0;
- const auto x1 = (xs[xs.size() - 2] + xs[xs.size() - 1]) / 2.0;
- if (x1 > currentMax)
- currentMax = x1;
- }
- } else {
- for (size_t i = 0; i < ws.getNumberHistograms(); ++i) {
- const auto &xs = ws.x(i);
- const auto x0 = xs.front();
- if (x0 < currentMin)
- currentMin = x0;
- const auto x1 = xs.back();
- if (x1 > currentMax)
- currentMax = x1;
- }
+ for (size_t i = 0; i < ws.getNumberHistograms(); ++i) {
+ const auto h = ws.histogram(i);
+ const auto first = h.begin();
+ currentMin = std::min(first->center(), currentMin);
+ const auto last = std::prev(h.end());
+ currentMax = std::max(last->center(), currentMax);
}
return std::tie(currentMin, currentMax);
}
@@ -141,9 +125,12 @@ modelHistogram(const API::MatrixWorkspace &modelWS,
if (wavelengthPoints > 1) {
const double step = (maxWavelength - minWavelength) /
static_cast<double>(wavelengthPoints - 1);
- for (size_t i = 0; i < xs.size(); ++i) {
+ for (size_t i = 0; i < xs.size() - 1; ++i) {
xs[i] = minWavelength + step * static_cast<double>(i);
}
+ // Force last point as otherwise it might be slightly off due to
+ // small rounding errors in the calculation above.
+ xs.back() = maxWavelength;
} else {
xs.front() = (minWavelength + maxWavelength) / 2.0;
}
diff --git a/Framework/Algorithms/test/SparseInstrumentTest.h b/Framework/Algorithms/test/SparseInstrumentTest.h
index 3393ec7d9a7..8e15b80deb7 100644
--- a/Framework/Algorithms/test/SparseInstrumentTest.h
+++ b/Framework/Algorithms/test/SparseInstrumentTest.h
@@ -8,6 +8,7 @@
#include "MantidDataObjects/Workspace2D.h"
#include "MantidDataObjects/WorkspaceCreation.h"
#include "MantidHistogramData/Histogram.h"
+#include "MantidHistogramData/LinearGenerator.h"
#include "MantidGeometry/Instrument/ReferenceFrame.h"
#include "MantidGeometry/Instrument.h"
#include "MantidTestHelpers/WorkspaceCreationHelper.h"
@@ -279,6 +280,21 @@ public:
TS_ASSERT_LESS_THAN(lon, grid->longitudeAt(1))
}
+ void test_modelHistogram_coversModelWS() {
+ using namespace Mantid::DataObjects;
+ using namespace Mantid::HistogramData;
+ const BinEdges edges(256, LinearGenerator(-1.33, 0.77));
+ const Counts counts(edges.size() - 1, 0.);
+ auto ws = create<Workspace2D>(2, Histogram(edges, counts));
+ const auto points = ws->points(0);
+ for (size_t nCounts = 2; nCounts < counts.size(); ++nCounts) {
+ const auto histo = modelHistogram(*ws, nCounts);
+ // These have to be equal, don't use DELTA here!
+ TS_ASSERT_EQUALS(histo.x().front(), points.front())
+ TS_ASSERT_EQUALS(histo.x().back(), points.back())
+ }
+ }
+
private:
const Mantid::Geometry::ReferenceFrame m_goofyRefFrame;
const Mantid::Geometry::ReferenceFrame m_standardRefFrame;
diff --git a/docs/source/release/v3.13.0/direct_inelastic.rst b/docs/source/release/v3.13.0/direct_inelastic.rst
index c8f194a9df6..525292d0926 100644
--- a/docs/source/release/v3.13.0/direct_inelastic.rst
+++ b/docs/source/release/v3.13.0/direct_inelastic.rst
@@ -30,6 +30,7 @@ Bug fixes
#########
- Fixed a crash in :ref:`SofQW <algm-SofQW>`, :ref:`SofQWCentre <algm-SofQWCentre>`, :ref:`SofQWNormalisedPolygon <algm-SofQWNormalisedPolygon>` and :ref:`SofQWPolygon <algm-SofQWPolygon>` algorithms when they were supplied with energy or :math:`Q` binning params containing the bin width only.
+- Fixed a failure in the wavelength interpolation of :ref:`MonteCarloAbsorption <algm-MonteCarloAbsorption>` which occurred under certain input property combinations.
Instrument Definitions
----------------------
diff --git a/docs/source/release/v3.13.0/indirect_inelastic.rst b/docs/source/release/v3.13.0/indirect_inelastic.rst
index c1f279e65e7..9d9ab894795 100644
--- a/docs/source/release/v3.13.0/indirect_inelastic.rst
+++ b/docs/source/release/v3.13.0/indirect_inelastic.rst
@@ -45,3 +45,5 @@ Bugfixes
- The MSDFit algorithm now uses the fully specified model given in the interface; previously MSDFit only used the
model specified in the 'Fit Type' drop-down menu.
+- Fixed a failure in the wavelength interpolation of :ref:`MonteCarloAbsorption <algm-MonteCarloAbsorption>` which occurred under certain input property combinations.
+
--
GitLab