diff --git a/Framework/Algorithms/test/SofQWCutTest.h b/Framework/Algorithms/test/SofQWCutTest.h index ea13b7476613cec13ef681ba983fd72af91300f7..60f86cc05a0ceaa4101f5f2029446d09eac94a50 100644 --- a/Framework/Algorithms/test/SofQWCutTest.h +++ b/Framework/Algorithms/test/SofQWCutTest.h @@ -12,6 +12,7 @@ #include "MantidAPI/WorkspaceGroup.h" #include "MantidKernel/Unit.h" #include "MantidDataHandling/LoadNexusProcessed.h" +#include "MantidDataHandling/CreateSimulationWorkspace.h" #include <boost/make_shared.hpp> @@ -221,6 +222,47 @@ public: TS_ASSERT_DELTA(ws_e->readY(0)[119], 0.045373095, delta); TS_ASSERT_DELTA(ws_e->readE(0)[119], 0.005462714, delta); } + + void test_sofqw3_zerobinwidth() { + // This test sets up a workspace which can yield a bin with zero width + // to check the code FractionalRebinning code handles this correctly + const double delta(1e-08); + Mantid::DataHandling::CreateSimulationWorkspace create_ws; + create_ws.initialize(); + create_ws.setChild(true); + create_ws.setPropertyValue("Instrument", "MARI"); + create_ws.setPropertyValue("BinParams", "-5,0.5,24"); + create_ws.setPropertyValue("OutputWorkspace", "__unused"); + create_ws.execute(); + + MatrixWorkspace_sptr createdWS = create_ws.getProperty("OutputWorkspace"); + auto inWS = boost::dynamic_pointer_cast<MatrixWorkspace>(createdWS); + // Sets one spectrum to zero so final value is not unity. + inWS->setCounts(300, std::vector<double>(58, 0.)); + + Mantid::Algorithms::SofQWNormalisedPolygon sqw; + sqw.initialize(); + sqw.setChild(true); + TS_ASSERT_THROWS_NOTHING(sqw.setProperty("InputWorkspace", inWS)); + TS_ASSERT_THROWS_NOTHING( + sqw.setPropertyValue("OutputWorkspace", "__unused")); + TS_ASSERT_THROWS_NOTHING(sqw.setPropertyValue("QAxisBinning", "0,10,10")); + TS_ASSERT_THROWS_NOTHING(sqw.setPropertyValue("EMode", "Direct")); + TS_ASSERT_THROWS_NOTHING(sqw.setPropertyValue("EFixed", "25")); + TS_ASSERT_THROWS_NOTHING( + sqw.setPropertyValue("EAxisBinning", "-1.5,3,1.5")); + TS_ASSERT_THROWS_NOTHING(sqw.execute()); + TS_ASSERT(sqw.isExecuted()); + + MatrixWorkspace_sptr ws = sqw.getProperty("OutputWorkspace"); + TS_ASSERT_EQUALS(ws->getAxis(0)->length(), 2); + TS_ASSERT_EQUALS(ws->getAxis(0)->unit()->unitID(), "DeltaE"); + TS_ASSERT_EQUALS((*(ws->getAxis(0)))(0), -1.5); + TS_ASSERT_EQUALS((*(ws->getAxis(0)))(1), 1.5); + TS_ASSERT_EQUALS(ws->getAxis(1)->length(), 2); + TS_ASSERT_EQUALS(ws->getAxis(1)->unit()->unitID(), "MomentumTransfer"); + TS_ASSERT_DELTA(ws->readY(0)[0], 0.998910675, delta); + } }; #endif /*SOFQWCUTTEST_H_*/ diff --git a/Framework/DataObjects/src/FractionalRebinning.cpp b/Framework/DataObjects/src/FractionalRebinning.cpp index ff317cc6cce6395dd66c2e835c69e0ea3bd6092b..6e9ee3b1a49dbdad84d1246ae4a98979182efdbc 100644 --- a/Framework/DataObjects/src/FractionalRebinning.cpp +++ b/Framework/DataObjects/src/FractionalRebinning.cpp @@ -152,9 +152,9 @@ double polyArea(T &v1, T &v2, Ts &&... vertices) { * @param yAxis The output data vertical axis * @param inputQ The input quadrilateral * @param y_start The starting y-axis index - * @param y_end The starting y-axis index + * @param y_end The ending y-axis index * @param x_start The starting x-axis index - * @param x_end The starting x-axis index + * @param x_end The ending x-axis index * @param areaInfo Output vector of indices and areas of overlapping bins */ void calcTrapezoidYIntersections( @@ -202,8 +202,8 @@ void calcTrapezoidYIntersections( // Step 1 - construct the left/right bin lims on the lines of the y-grid. const double NaN = std::numeric_limits<double>::quiet_NaN(); const double DBL_EPS = std::numeric_limits<double>::epsilon(); - std::vector<double> leftLim(nx * (ny + 1), NaN); - std::vector<double> rightLim(nx * (ny + 1), NaN); + std::vector<double> leftLim((nx + 1) * (ny + 1), NaN); + std::vector<double> rightLim((nx + 1) * (ny + 1), NaN); auto x0_it = xAxis.begin() + x_start; auto x1_it = xAxis.begin() + x_end + 1; auto y0_it = yAxis.begin() + y_start; @@ -239,7 +239,7 @@ void calcTrapezoidYIntersections( auto right_it = std::upper_bound(x0_it, x1_it, right_val); if (right_it != x1_it) { rightLim[right_it - x0_it - 1 + yjx] = right_val; - } else if (yAxis[yj + y_start] < ur_y) { + } else if (yAxis[yj + y_start] < ur_y && nx > 0) { right_it = x1_it - 1; rightLim[nx - 1 + yjx] = lr_x; } @@ -277,12 +277,15 @@ void calcTrapezoidYIntersections( left_it--; if (right_it == x1_it) right_it--; - size_t li = left_it - x0_it - 1; + size_t li = (left_it > x0_it) ? (left_it - x0_it - 1) : 0; + size_t ri = (right_it > x0_it) ? (right_it - x0_it - 1) : 0; leftLim[li + yjx] = (mTop >= 0) ? val : ll_x; - rightLim[right_it - x0_it - 1 + yjx] = (mTop >= 0) ? lr_x : val; - for (auto x_it = left_it; x_it != right_it; x_it++) { - leftLim[li + 1 + yjx] = *x_it; - rightLim[li++ + yjx] = *x_it; + rightLim[ri + yjx] = (mTop >= 0) ? lr_x : val; + if (left_it < right_it && right_it != x1_it) { + for (auto x_it = left_it; x_it != right_it; x_it++) { + leftLim[li + 1 + yjx] = *x_it; + rightLim[li++ + yjx] = *x_it; + } } } }