Merge pull request #22194 from mantidproject/22191_sofqw3_segfault

SofQW3 segfault

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_*/

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;
+        }
       }
     }
   }