diff --git a/Framework/Algorithms/inc/MantidAlgorithms/SampleCorrections/MCAbsorptionStrategy.h b/Framework/Algorithms/inc/MantidAlgorithms/SampleCorrections/MCAbsorptionStrategy.h
index 4c97318f80d275ad51f3bf528343331566225a65..5d1bc89e2d9ee5f1fe303eb7856842e56d37252a 100644
--- a/Framework/Algorithms/inc/MantidAlgorithms/SampleCorrections/MCAbsorptionStrategy.h
+++ b/Framework/Algorithms/inc/MantidAlgorithms/SampleCorrections/MCAbsorptionStrategy.h
@@ -38,8 +38,8 @@ public:
                        size_t maxScatterPtAttempts);
   std::tuple<double, double> calculate(Kernel::PseudoRandomNumberGenerator &rng,
                                        const Kernel::V3D &finalPos,
-                                       double lambdaBefore,
-                                       double lambdaAfter) const;
+                                       double lambdaBefore, double lambdaAfter,
+                                       std::string &debugString) const;
 
 private:
   const IBeamProfile &m_beamProfile;
diff --git a/Framework/Algorithms/inc/MantidAlgorithms/SampleCorrections/MCInteractionVolume.h b/Framework/Algorithms/inc/MantidAlgorithms/SampleCorrections/MCInteractionVolume.h
index c761b8740c8c6cbc41f2350aaec8f7ae57075b66..a93d3ad1b61e896328c509f2f572f5b5cdaae3d0 100644
--- a/Framework/Algorithms/inc/MantidAlgorithms/SampleCorrections/MCInteractionVolume.h
+++ b/Framework/Algorithms/inc/MantidAlgorithms/SampleCorrections/MCInteractionVolume.h
@@ -46,8 +46,12 @@ public:
                              const Kernel::V3D &startPos,
                              const Kernel::V3D &endPos, double lambdaBefore,
                              double lambdaAfter) const;
+  std::string generateScatterPointStats() const;
+  Kernel::V3D generatePoint(Kernel::PseudoRandomNumberGenerator &rng) const;
 
 private:
+  mutable int m_sampleScatterPoints = 0;
+  mutable std::vector<int> m_envScatterPoints;
   const boost::shared_ptr<Geometry::IObject> m_sample;
   const Geometry::SampleEnvironment *m_env;
   const Geometry::BoundingBox m_activeRegion;
diff --git a/Framework/Algorithms/src/MonteCarloAbsorption.cpp b/Framework/Algorithms/src/MonteCarloAbsorption.cpp
index 777a569ad2a9cf6924f84cf541646bdabd135613..ea990c3c901d05507d30311f9dd08d143fbc98a3 100644
--- a/Framework/Algorithms/src/MonteCarloAbsorption.cpp
+++ b/Framework/Algorithms/src/MonteCarloAbsorption.cpp
@@ -276,8 +276,10 @@ MatrixWorkspace_uptr MonteCarloAbsorption::doSimulation(
       } else {
         // elastic case already initialized
       }
+      std::string debugString;
       std::tie(outY[j], std::ignore) =
-          strategy.calculate(rng, detPos, lambdaIn, lambdaOut);
+          strategy.calculate(rng, detPos, lambdaIn, lambdaOut, debugString);
+      g_log.debug(debugString);
 
       // Ensure we have the last point for the interpolation
       if (lambdaStepSize > 1 && j + lambdaStepSize >= nbins && j + 1 != nbins) {
diff --git a/Framework/Algorithms/src/SampleCorrections/MCAbsorptionStrategy.cpp b/Framework/Algorithms/src/SampleCorrections/MCAbsorptionStrategy.cpp
index 4ead5a72d5c80820a20e1d36c2c1ea0e7cb04bc7..8858f1e9b69138eed602faca4fb61c5c5ca50ad0 100644
--- a/Framework/Algorithms/src/SampleCorrections/MCAbsorptionStrategy.cpp
+++ b/Framework/Algorithms/src/SampleCorrections/MCAbsorptionStrategy.cpp
@@ -43,12 +43,12 @@ MCAbsorptionStrategy::MCAbsorptionStrategy(const IBeamProfile &beamProfile,
  * where it is detected
  * @param lambdaBefore Wavelength, in \f$\\A^-1\f$, before scattering
  * @param lambdaAfter Wavelength, in \f$\\A^-1\f$, after scattering
+ * @param debugString string to return any debug messages
  * @return A tuple of the <correction factor, associated error>.
  */
-std::tuple<double, double>
-MCAbsorptionStrategy::calculate(Kernel::PseudoRandomNumberGenerator &rng,
-                                const Kernel::V3D &finalPos,
-                                double lambdaBefore, double lambdaAfter) const {
+std::tuple<double, double> MCAbsorptionStrategy::calculate(
+    Kernel::PseudoRandomNumberGenerator &rng, const Kernel::V3D &finalPos,
+    double lambdaBefore, double lambdaAfter, std::string &debugString) const {
   const auto scatterBounds = m_scatterVol.getBoundingBox();
   double factor(0.0);
   for (size_t i = 0; i < m_nevents; ++i) {
@@ -74,6 +74,9 @@ MCAbsorptionStrategy::calculate(Kernel::PseudoRandomNumberGenerator &rng,
       }
     } while (true);
   }
+
+  debugString += m_scatterVol.generateScatterPointStats();
+
   using std::make_tuple;
   return make_tuple(factor / static_cast<double>(m_nevents), m_error);
 }
diff --git a/Framework/Algorithms/src/SampleCorrections/MCInteractionVolume.cpp b/Framework/Algorithms/src/SampleCorrections/MCInteractionVolume.cpp
index 030e13b175a097a4e0277e4b76f4eb80381f7cc8..62d8c874eb0a34540b95291935d4aae62eaa9dda 100644
--- a/Framework/Algorithms/src/SampleCorrections/MCInteractionVolume.cpp
+++ b/Framework/Algorithms/src/SampleCorrections/MCInteractionVolume.cpp
@@ -11,6 +11,7 @@
 #include "MantidGeometry/Objects/Track.h"
 #include "MantidKernel/Material.h"
 #include "MantidKernel/PseudoRandomNumberGenerator.h"
+#include <iomanip>
 
 namespace Mantid {
 using Geometry::Track;
@@ -46,6 +47,10 @@ MCInteractionVolume::MCInteractionVolume(
   } catch (std::runtime_error &) {
     // swallow this as no defined environment from getEnvironment
   }
+
+  if (m_env) {
+    m_envScatterPoints.resize(m_env->nelements());
+  }
 }
 
 /**
@@ -56,6 +61,49 @@ const Geometry::BoundingBox &MCInteractionVolume::getBoundingBox() const {
   return m_sample->getBoundingBox();
 }
 
+/**
+ * Generate point randomly across one of the components of the environment
+ * including the sample itself in the selection. The method first selects
+ * a random component and then selects a random point within that component
+ * using Object::generatePointObject
+ * @param rng A reference to a PseudoRandomNumberGenerator where
+ * nextValue should return a flat random number between 0.0 & 1.0
+ * @return The generated point
+ */
+Kernel::V3D MCInteractionVolume::generatePoint(
+    Kernel::PseudoRandomNumberGenerator &rng) const {
+  for (size_t i = 0; i < m_maxScatterAttempts; i++) {
+    Kernel::V3D point;
+    int componentIndex;
+    if (m_env) {
+      componentIndex = rng.nextInt(0, static_cast<int>(m_env->nelements())) - 1;
+    } else {
+      componentIndex = -1;
+    }
+
+    bool pointGenerated;
+    if (componentIndex == -1) {
+      pointGenerated =
+          m_sample->generatePointInObject(rng, m_activeRegion, 1, point);
+    } else {
+      pointGenerated =
+          m_env->getComponent(componentIndex)
+              .generatePointInObject(rng, m_activeRegion, 1, point);
+    }
+    if (pointGenerated) {
+      if (componentIndex == -1) {
+        m_sampleScatterPoints++;
+      } else {
+        m_envScatterPoints[componentIndex]++;
+      }
+      return point;
+    }
+  }
+  throw std::runtime_error("MCInteractionVolume::generatePoint() - Unable to "
+                           "generate point in object after " +
+                           std::to_string(m_maxScatterAttempts) + " attempts");
+}
+
 /**
  * Calculate the attenuation correction factor the volume given a start and
  * end point.
@@ -78,14 +126,8 @@ double MCInteractionVolume::calculateAbsorption(
   // is calculated in reverse, i.e. defining the track from the scatter pt
   // backwards for simplicity with how the Track object works. This avoids
   // having to understand exactly which object the scattering occurred in.
-  V3D scatterPos;
-  if (m_env && (rng.nextValue() > 0.5)) {
-    scatterPos =
-        m_env->generatePoint(rng, m_activeRegion, m_maxScatterAttempts);
-  } else {
-    scatterPos = m_sample->generatePointInObject(rng, m_activeRegion,
-                                                 m_maxScatterAttempts);
-  }
+  V3D scatterPos = generatePoint(rng);
+
   const auto toStart = normalize(startPos - scatterPos);
   Track beforeScatter(scatterPos, toStart);
   int nlinks = m_sample->interceptSurface(beforeScatter);
@@ -120,5 +162,40 @@ double MCInteractionVolume::calculateAbsorption(
          calculateAttenuation(afterScatter, lambdaAfter);
 }
 
+/**
+ * Generate a string summarising which parts of the environment
+ * the simulated scatter points occurred in
+ * @return The generated string
+ */
+std::string MCInteractionVolume::generateScatterPointStats() const {
+  std::stringstream scatterPointSummary;
+  scatterPointSummary << std::fixed;
+  scatterPointSummary << std::setprecision(2);
+
+  scatterPointSummary << "Scatter point counts:" << std::endl;
+
+  int totalScatterPoints = m_sampleScatterPoints;
+  for (std::vector<int>::size_type i = 0; i < m_envScatterPoints.size(); i++) {
+    totalScatterPoints += m_envScatterPoints[i];
+  }
+  scatterPointSummary << "Total scatter points: " << totalScatterPoints
+                      << std::endl;
+
+  double percentage =
+      static_cast<double>(m_sampleScatterPoints) / totalScatterPoints * 100;
+  scatterPointSummary << "Sample: " << m_sampleScatterPoints << " ("
+                      << percentage << "%)" << std::endl;
+
+  for (std::vector<int>::size_type i = 0; i < m_envScatterPoints.size(); i++) {
+    percentage =
+        static_cast<double>(m_envScatterPoints[i]) / totalScatterPoints * 100;
+    scatterPointSummary << "Environment part " << i << " ("
+                        << m_env->getComponent(i).id()
+                        << "): " << m_envScatterPoints[i] << " (" << percentage
+                        << "%)" << std::endl;
+  }
+  return scatterPointSummary.str();
+}
+
 } // namespace Algorithms
 } // namespace Mantid
diff --git a/Framework/Algorithms/test/MCAbsorptionStrategyTest.h b/Framework/Algorithms/test/MCAbsorptionStrategyTest.h
index c23ed83a607d67cf16349cc07d36b3dacb84db96..31fcefa7323a45d1e1976f097f7598fd353fe191 100644
--- a/Framework/Algorithms/test/MCAbsorptionStrategyTest.h
+++ b/Framework/Algorithms/test/MCAbsorptionStrategyTest.h
@@ -57,8 +57,9 @@ public:
     const double lambdaBefore(2.5), lambdaAfter(3.5);
 
     double factor(0.0), error(0.0);
+    std::string debugString;
     std::tie(factor, error) =
-        mcabsorb.calculate(rng, endPos, lambdaBefore, lambdaAfter);
+        mcabsorb.calculate(rng, endPos, lambdaBefore, lambdaAfter, debugString);
     TS_ASSERT_DELTA(0.0043828472, factor, 1e-08);
     TS_ASSERT_DELTA(1.0 / std::sqrt(nevents), error, 1e-08);
   }
@@ -85,8 +86,10 @@ public:
     EXPECT_CALL(rng, nextValue()).WillRepeatedly(Return(0.5));
     const double lambdaBefore(2.5), lambdaAfter(3.5);
     const V3D endPos(0.7, 0.7, 1.4);
-    TS_ASSERT_THROWS(mcabs.calculate(rng, endPos, lambdaBefore, lambdaAfter),
-                     const std::runtime_error &)
+    std::string debugString;
+    TS_ASSERT_THROWS(
+        mcabs.calculate(rng, endPos, lambdaBefore, lambdaAfter, debugString),
+        const std::runtime_error &)
   }
 
 private:
diff --git a/Framework/Algorithms/test/MCInteractionVolumeTest.h b/Framework/Algorithms/test/MCInteractionVolumeTest.h
index fc78ba3f16ba297bcb78b77274ca5a4d8845c622..ccf4ce62e9de9d9208ec41be0aa77020b33fa613 100644
--- a/Framework/Algorithms/test/MCInteractionVolumeTest.h
+++ b/Framework/Algorithms/test/MCInteractionVolumeTest.h
@@ -104,14 +104,15 @@ public:
 
     auto sample = createTestSample(TestSampleType::SamplePlusContainer);
     MockRNG rng;
-    // force scatter in segment can
-    EXPECT_CALL(rng, nextInt(1, 1)).Times(Exactly(1)).WillOnce(Return(1));
+    // force scatter in segment can (1)
+    EXPECT_CALL(rng, nextInt(0, 1))
+        .Times(Exactly(1))
+        .WillOnce(Return(1)); // MCInteractionVolume::generatePoint
     EXPECT_CALL(rng, nextValue())
-        .Times(Exactly(4))
-        .WillOnce(Return(0.75))
-        .WillOnce(Return(0.02))
-        .WillOnce(Return(0.5))
-        .WillOnce(Return(0.5));
+        .Times(Exactly(3))
+        .WillOnce(Return(0.02)) // RandomPoint::bounded r1
+        .WillOnce(Return(0.5))  // r2
+        .WillOnce(Return(0.5)); // r3
 
     MCInteractionVolume interactor(sample,
                                    sample.getEnvironment().boundingBox());
@@ -120,11 +121,13 @@ public:
     TS_ASSERT_DELTA(0.69223681, factorContainer, 1e-8);
     Mock::VerifyAndClearExpectations(&rng);
 
-    // force scatter in sample
+    // force scatter in sample (0)
+    EXPECT_CALL(rng, nextInt(0, 1))
+        .Times(Exactly(1))
+        .WillOnce(Return(0)); // MCInteractionVolume::generatePoint
     EXPECT_CALL(rng, nextValue())
-        .Times(Exactly(4))
-        .WillOnce(Return(0.25))
-        .WillRepeatedly(Return(0.25));
+        .Times(Exactly(3))
+        .WillRepeatedly(Return(0.25)); // RandomPoint::bounded r1, r2, r3
     const double factorSample = interactor.calculateAbsorption(
         rng, startPos, endPos, lambdaBefore, lambdaAfter);
     TS_ASSERT_DELTA(0.73100698, factorSample, 1e-8);
@@ -167,6 +170,97 @@ public:
                                                     lambdaBefore, lambdaAfter),
                      const std::runtime_error &);
   }
+
+  void testGeneratePointConsidersAllComponents() {
+    using namespace ::testing;
+    using namespace Mantid::Kernel;
+    using namespace MonteCarloTesting;
+
+    auto kit = createTestKit();
+    size_t maxAttempts(1);
+
+    using Mantid::API::Sample;
+    Sample sample;
+    sample.setShape(ComponentCreationHelper::createSphere(1));
+    sample.setEnvironment(
+        std::make_unique<Mantid::Geometry::SampleEnvironment>(*kit));
+
+    MCInteractionVolume interactor(
+        sample, kit->getComponent(0).getBoundingBox(), maxAttempts);
+
+    // Generate "random" sequence
+    MockRNG rng;
+    // Selects first component
+    EXPECT_CALL(rng, nextInt(_, _)).Times(Exactly(1)).WillOnce(Return(1));
+    EXPECT_CALL(rng, nextValue())
+        .Times(3)
+        .WillOnce(Return(0.55))
+        .WillOnce(Return(0.65))
+        .WillOnce(Return(0.70));
+    V3D comp1Point;
+    TS_ASSERT_THROWS_NOTHING(comp1Point = interactor.generatePoint(rng));
+    TS_ASSERT(kit->isValid(comp1Point));
+    Mock::VerifyAndClearExpectations(&rng);
+
+    // Selects second component
+    MCInteractionVolume interactor2(
+        sample, kit->getComponent(1).getBoundingBox(), maxAttempts);
+
+    EXPECT_CALL(rng, nextInt(_, _)).Times(Exactly(1)).WillOnce(Return(2));
+    EXPECT_CALL(rng, nextValue())
+        .Times(3)
+        .WillOnce(Return(0.55))
+        .WillOnce(Return(0.65))
+        .WillOnce(Return(0.70));
+    V3D comp2Point;
+    TS_ASSERT_THROWS_NOTHING(comp2Point = interactor2.generatePoint(rng));
+    TS_ASSERT(comp2Point != comp1Point);
+    TS_ASSERT(kit->isValid(comp2Point));
+    Mock::VerifyAndClearExpectations(&rng);
+
+    // Selects third component
+    MCInteractionVolume interactor3(
+        sample, kit->getComponent(2).getBoundingBox(), maxAttempts);
+    EXPECT_CALL(rng, nextInt(_, _)).Times(Exactly(1)).WillOnce(Return(3));
+    EXPECT_CALL(rng, nextValue())
+        .Times(3)
+        .WillOnce(Return(0.55))
+        .WillOnce(Return(0.65))
+        .WillOnce(Return(0.70));
+    V3D comp3Point;
+    TS_ASSERT_THROWS_NOTHING(comp3Point = interactor3.generatePoint(rng));
+    TS_ASSERT(comp3Point != comp2Point);
+    TS_ASSERT(comp3Point != comp1Point);
+    TS_ASSERT(kit->isValid(comp3Point));
+    Mock::VerifyAndClearExpectations(&rng);
+  }
+
+  void testGeneratePointRespectsActiveRegion() {
+    using namespace MonteCarloTesting;
+    using namespace ::testing;
+    using namespace Mantid::Kernel;
+
+    auto kit = createTestKit();
+    size_t maxAttempts(1);
+
+    // Generate "random" sequence
+    MockRNG rng;
+    // Sequence will try to select one of the non-container pieces
+    EXPECT_CALL(rng, nextInt(_, _)).Times(Exactly(1)).WillOnce(Return(2));
+    EXPECT_CALL(rng, nextValue()).Times(3).WillRepeatedly(Return(0.5));
+
+    using Mantid::API::Sample;
+    Sample sample;
+    sample.setShape(ComponentCreationHelper::createSphere(1));
+    sample.setEnvironment(
+        std::make_unique<Mantid::Geometry::SampleEnvironment>(*kit));
+
+    MCInteractionVolume interactor(sample, kit->getContainer().getBoundingBox(),
+                                   maxAttempts);
+    // Restrict region to can
+    TS_ASSERT_THROWS(interactor.generatePoint(rng), const std::runtime_error &);
+    Mock::VerifyAndClearExpectations(&rng);
+  }
 };
 
 #endif /* MANTID_ALGORITHMS_MCINTERACTIONVOLUMETEST_H_ */
diff --git a/Framework/Algorithms/test/MonteCarloAbsorptionTest.h b/Framework/Algorithms/test/MonteCarloAbsorptionTest.h
index cad25ce29122027a88b9d49753c819179f1928b2..cf9b4c634764a3d311e4fa93ccdb42e4d0fc8025 100644
--- a/Framework/Algorithms/test/MonteCarloAbsorptionTest.h
+++ b/Framework/Algorithms/test/MonteCarloAbsorptionTest.h
@@ -117,18 +117,18 @@ public:
     auto outputWS = runAlgorithm(wsProps);
 
     verifyDimensions(wsProps, outputWS);
-    const double delta(1e-05);
+    const double delta(1e-04);
     const size_t middle_index(4);
 
-    TS_ASSERT_DELTA(0.6245262704, outputWS->y(0).front(), delta);
-    TS_ASSERT_DELTA(0.2770105008, outputWS->y(0)[middle_index], delta);
-    TS_ASSERT_DELTA(0.1041517761, outputWS->y(0).back(), delta);
-    TS_ASSERT_DELTA(0.6282072570, outputWS->y(2).front(), delta);
-    TS_ASSERT_DELTA(0.2790911215, outputWS->y(2)[middle_index], delta);
-    TS_ASSERT_DELTA(0.1078704876, outputWS->y(2).back(), delta);
-    TS_ASSERT_DELTA(0.6267458002, outputWS->y(4).front(), delta);
-    TS_ASSERT_DELTA(0.2790655428, outputWS->y(4)[middle_index], delta);
-    TS_ASSERT_DELTA(0.1068945921, outputWS->y(4).back(), delta);
+    TS_ASSERT_DELTA(0.6243, outputWS->y(0).front(), delta);
+    TS_ASSERT_DELTA(0.2770, outputWS->y(0)[middle_index], delta);
+    TS_ASSERT_DELTA(0.1042, outputWS->y(0).back(), delta);
+    TS_ASSERT_DELTA(0.6280, outputWS->y(2).front(), delta);
+    TS_ASSERT_DELTA(0.2791, outputWS->y(2)[middle_index], delta);
+    TS_ASSERT_DELTA(0.1079, outputWS->y(2).back(), delta);
+    TS_ASSERT_DELTA(0.6265, outputWS->y(4).front(), delta);
+    TS_ASSERT_DELTA(0.2791, outputWS->y(4)[middle_index], delta);
+    TS_ASSERT_DELTA(0.1069, outputWS->y(4).back(), delta);
   }
 
   void test_Workspace_With_Just_Sample_For_Direct() {
@@ -138,12 +138,12 @@ public:
     auto outputWS = runAlgorithm(wsProps);
 
     verifyDimensions(wsProps, outputWS);
-    const double delta(1e-05);
+    const double delta(1e-04);
     const size_t middle_index(4);
 
-    TS_ASSERT_DELTA(0.5060586383, outputWS->y(0).front(), delta);
-    TS_ASSERT_DELTA(0.3389572854, outputWS->y(0)[middle_index], delta);
-    TS_ASSERT_DELTA(0.2186106403, outputWS->y(0).back(), delta);
+    TS_ASSERT_DELTA(0.5061, outputWS->y(0).front(), delta);
+    TS_ASSERT_DELTA(0.3390, outputWS->y(0)[middle_index], delta);
+    TS_ASSERT_DELTA(0.2186, outputWS->y(0).back(), delta);
   }
 
   void test_Workspace_With_Just_Sample_For_Indirect() {
@@ -153,12 +153,12 @@ public:
     auto outputWS = runAlgorithm(wsProps);
 
     verifyDimensions(wsProps, outputWS);
-    const double delta(1e-05);
+    const double delta(1e-04);
     const size_t middle_index(4);
 
-    TS_ASSERT_DELTA(0.3658064669, outputWS->y(0).front(), delta);
-    TS_ASSERT_DELTA(0.2277021116, outputWS->y(0)[middle_index], delta);
-    TS_ASSERT_DELTA(0.1383621690, outputWS->y(0).back(), delta);
+    TS_ASSERT_DELTA(0.3651, outputWS->y(0).front(), delta);
+    TS_ASSERT_DELTA(0.2277, outputWS->y(0)[middle_index], delta);
+    TS_ASSERT_DELTA(0.1384, outputWS->y(0).back(), delta);
   }
 
   void test_Workspace_With_Sample_And_Container() {
@@ -168,12 +168,12 @@ public:
     auto outputWS = runAlgorithm(wsProps);
 
     verifyDimensions(wsProps, outputWS);
-    const double delta(1e-05);
+    const double delta(1e-04);
     const size_t middle_index(4);
 
-    TS_ASSERT_DELTA(0.6028577409, outputWS->y(0).front(), delta);
-    TS_ASSERT_DELTA(0.2758029301, outputWS->y(0)[middle_index], delta);
-    TS_ASSERT_DELTA(0.1035877536, outputWS->y(0).back(), delta);
+    TS_ASSERT_DELTA(0.5995, outputWS->y(0).front(), delta);
+    TS_ASSERT_DELTA(0.2688, outputWS->y(0)[middle_index], delta);
+    TS_ASSERT_DELTA(0.0981, outputWS->y(0).back(), delta);
   }
 
   void test_Workspace_Beam_Size_Set() {
@@ -183,11 +183,11 @@ public:
     auto outputWS = runAlgorithm(wsProps);
 
     verifyDimensions(wsProps, outputWS);
-    const double delta(1e-05);
+    const double delta(1e-04);
     const size_t middle_index(4);
-    TS_ASSERT_DELTA(0.6279328067, outputWS->y(0).front(), delta);
-    TS_ASSERT_DELTA(0.2770103841, outputWS->y(0)[middle_index], delta);
-    TS_ASSERT_DELTA(0.1041522173, outputWS->y(0).back(), delta);
+    TS_ASSERT_DELTA(0.6243, outputWS->y(0).front(), delta);
+    TS_ASSERT_DELTA(0.2770, outputWS->y(0)[middle_index], delta);
+    TS_ASSERT_DELTA(0.1041, outputWS->y(0).back(), delta);
   }
 
   void test_Linear_Interpolation() {
@@ -199,11 +199,11 @@ public:
     auto outputWS = runAlgorithm(wsProps, nlambda, interpolation);
 
     verifyDimensions(wsProps, outputWS);
-    const double delta(1e-05);
-    TS_ASSERT_DELTA(0.6245262704, outputWS->y(0).front(), delta);
-    TS_ASSERT_DELTA(0.3455351271, outputWS->y(0)[3], delta);
-    TS_ASSERT_DELTA(0.2783583898, outputWS->y(0)[4], delta);
-    TS_ASSERT_DELTA(0.1168965453, outputWS->y(0).back(), delta);
+    const double delta(1e-04);
+    TS_ASSERT_DELTA(0.6243, outputWS->y(0).front(), delta);
+    TS_ASSERT_DELTA(0.3455, outputWS->y(0)[3], delta);
+    TS_ASSERT_DELTA(0.2784, outputWS->y(0)[4], delta);
+    TS_ASSERT_DELTA(0.1169, outputWS->y(0).back(), delta);
   }
 
   void test_CSpline_Interpolation() {
@@ -215,12 +215,12 @@ public:
     auto outputWS = runAlgorithm(wsProps, nlambda, interpolation);
 
     verifyDimensions(wsProps, outputWS);
-    const double delta(1e-05);
-    TS_ASSERT_DELTA(0.6245262704, outputWS->y(0).front(), delta);
+    const double delta(1e-04);
+    TS_ASSERT_DELTA(0.6243, outputWS->y(0).front(), delta);
     // Interpolation gives some negative value due to test setup
-    TS_ASSERT_DELTA(0.3368552576, outputWS->y(0)[3], delta);
-    TS_ASSERT_DELTA(0.2783583898, outputWS->y(0)[4], delta);
-    TS_ASSERT_DELTA(0.1168965453, outputWS->y(0).back(), delta);
+    TS_ASSERT_DELTA(0.3369, outputWS->y(0)[3], delta);
+    TS_ASSERT_DELTA(0.2784, outputWS->y(0)[4], delta);
+    TS_ASSERT_DELTA(0.1169, outputWS->y(0).back(), delta);
   }
 
   //---------------------------------------------------------------------------
@@ -283,17 +283,17 @@ public:
     auto outputWS = runAlgorithm(wsProps, 5, "Linear", true, 3, 3);
 
     verifyDimensions(wsProps, outputWS);
-    const double delta{1e-05};
+    const double delta{1e-04};
     const size_t middle_index{4};
-    TS_ASSERT_DELTA(0.6241295766, outputWS->y(0).front(), delta);
-    TS_ASSERT_DELTA(0.2830812429, outputWS->y(0)[middle_index], delta);
-    TS_ASSERT_DELTA(0.1055135467, outputWS->y(0).back(), delta);
-    TS_ASSERT_DELTA(0.6265926135, outputWS->y(2).front(), delta);
-    TS_ASSERT_DELTA(0.2829539939, outputWS->y(2)[middle_index], delta);
-    TS_ASSERT_DELTA(0.1078604888, outputWS->y(2).back(), delta);
-    TS_ASSERT_DELTA(0.6260645680, outputWS->y(4).front(), delta);
-    TS_ASSERT_DELTA(0.2829361634, outputWS->y(4)[middle_index], delta);
-    TS_ASSERT_DELTA(0.1073283604, outputWS->y(4).back(), delta);
+    TS_ASSERT_DELTA(0.6239, outputWS->y(0).front(), delta);
+    TS_ASSERT_DELTA(0.2831, outputWS->y(0)[middle_index], delta);
+    TS_ASSERT_DELTA(0.1055, outputWS->y(0).back(), delta);
+    TS_ASSERT_DELTA(0.6264, outputWS->y(2).front(), delta);
+    TS_ASSERT_DELTA(0.2830, outputWS->y(2)[middle_index], delta);
+    TS_ASSERT_DELTA(0.1079, outputWS->y(2).back(), delta);
+    TS_ASSERT_DELTA(0.6259, outputWS->y(4).front(), delta);
+    TS_ASSERT_DELTA(0.2829, outputWS->y(4)[middle_index], delta);
+    TS_ASSERT_DELTA(0.1073, outputWS->y(4).back(), delta);
   }
 
   void test_Sparse_Instrument_For_Direct() {
@@ -303,13 +303,12 @@ public:
     auto outputWS = runAlgorithm(wsProps, 5, "Linear", true, 3, 3);
 
     verifyDimensions(wsProps, outputWS);
-    const double delta(1e-05);
+    const double delta(1e-04);
     const size_t middle_index(4);
 
-    TS_ASSERT_DELTA(0.5055988099, outputWS->y(0).front(), delta);
-    TS_ASSERT_DELTA(0.3419146885, outputWS->y(0)[middle_index], delta);
-    const double delta2(1e-08);
-    TS_ASSERT_DELTA(0.2270203007, outputWS->y(0).back(), delta2);
+    TS_ASSERT_DELTA(0.5056, outputWS->y(0).front(), delta);
+    TS_ASSERT_DELTA(0.3419, outputWS->y(0)[middle_index], delta);
+    TS_ASSERT_DELTA(0.2270, outputWS->y(0).back(), delta);
   }
 
   void test_Sparse_Instrument_For_Indirect() {
@@ -319,12 +318,12 @@ public:
     auto outputWS = runAlgorithm(wsProps, 5, "Linear", true, 3, 3);
 
     verifyDimensions(wsProps, outputWS);
-    const double delta(1e-05);
+    const double delta(1e-04);
     const size_t middle_index(4);
 
-    TS_ASSERT_DELTA(0.3651989107, outputWS->y(0).front(), delta);
-    TS_ASSERT_DELTA(0.2308318553, outputWS->y(0)[middle_index], delta);
-    TS_ASSERT_DELTA(0.1385715148, outputWS->y(0).back(), delta);
+    TS_ASSERT_DELTA(0.3646, outputWS->y(0).front(), delta);
+    TS_ASSERT_DELTA(0.2308, outputWS->y(0)[middle_index], delta);
+    TS_ASSERT_DELTA(0.1386, outputWS->y(0).back(), delta);
   }
 
 private:
diff --git a/Framework/Algorithms/test/MonteCarloTesting.h b/Framework/Algorithms/test/MonteCarloTesting.h
index 4877acb2b9d95a614bf59b9f15faf6ad325ebf1d..c6ee52715ef474d6d8aa6bd0c342a6bacbaf7d30 100644
--- a/Framework/Algorithms/test/MonteCarloTesting.h
+++ b/Framework/Algorithms/test/MonteCarloTesting.h
@@ -152,6 +152,24 @@ inline Mantid::API::Sample createTestSample(TestSampleType sampleType) {
   }
   return testSample;
 }
+
+inline boost::shared_ptr<Mantid::Geometry::SampleEnvironment> createTestKit() {
+  using namespace Mantid::Geometry;
+  using namespace Mantid::Kernel;
+
+  // at centre
+  ShapeFactory factory;
+  auto can = boost::make_shared<Container>(factory.createShape(
+      ComponentCreationHelper::sphereXML(0.01, V3D(0, 0, 0), "sp-1")));
+  can->setID("8mm");
+  auto kit = boost::make_shared<SampleEnvironment>("TestKit", can);
+  // before sample
+  kit->add(ComponentCreationHelper::createSphere(0.1, V3D(-0.25, 0.0, 0.0)));
+  // after sample
+  kit->add(ComponentCreationHelper::createSphere(0.1, V3D(0.25, 0.0, 0.0)));
+  return kit;
+}
+
 } // namespace MonteCarloTesting
 
 #endif // MANTID_ALGORITHMS_MONTECARLOTESTING_H
diff --git a/Framework/Geometry/inc/MantidGeometry/Instrument/Container.h b/Framework/Geometry/inc/MantidGeometry/Instrument/Container.h
index ee05a300adf9225bbc2fed28da82c8e87b243c6c..a18ef6de79a2968c954791d1aa2bb831268f3354 100644
--- a/Framework/Geometry/inc/MantidGeometry/Instrument/Container.h
+++ b/Framework/Geometry/inc/MantidGeometry/Instrument/Container.h
@@ -79,14 +79,15 @@ public:
   int getPointInObject(Kernel::V3D &point) const override {
     return m_shape->getPointInObject(point);
   }
-  Kernel::V3D generatePointInObject(Kernel::PseudoRandomNumberGenerator &rng,
-                                    const size_t i) const override {
-    return m_shape->generatePointInObject(rng, i);
-  }
-  Kernel::V3D generatePointInObject(Kernel::PseudoRandomNumberGenerator &rng,
-                                    const BoundingBox &activeRegion,
-                                    const size_t i) const override {
-    return m_shape->generatePointInObject(rng, activeRegion, i);
+  bool generatePointInObject(Kernel::PseudoRandomNumberGenerator &rng,
+                             const size_t i,
+                             Kernel::V3D &point) const override {
+    return m_shape->generatePointInObject(rng, i, point);
+  }
+  bool generatePointInObject(Kernel::PseudoRandomNumberGenerator &rng,
+                             const BoundingBox &activeRegion, const size_t i,
+                             Kernel::V3D &point) const override {
+    return m_shape->generatePointInObject(rng, activeRegion, i, point);
   }
 
   detail::ShapeInfo::GeometryShape shape() const override {
diff --git a/Framework/Geometry/inc/MantidGeometry/Instrument/SampleEnvironment.h b/Framework/Geometry/inc/MantidGeometry/Instrument/SampleEnvironment.h
index bb47cf66788e3245dad9f86888f5e6c0bab95ef3..5198a15ca407be8d1c30dd18456bf794901b5229 100644
--- a/Framework/Geometry/inc/MantidGeometry/Instrument/SampleEnvironment.h
+++ b/Framework/Geometry/inc/MantidGeometry/Instrument/SampleEnvironment.h
@@ -51,14 +51,7 @@ public:
   const IObject &getComponent(const size_t index) const;
 
   Geometry::BoundingBox boundingBox() const;
-  /// Select a random point within a component
-  Kernel::V3D generatePoint(Kernel::PseudoRandomNumberGenerator &rng,
-                            const size_t maxAttempts) const;
-  /// Select a random point within a component that is also bound by a
-  /// given region
-  Kernel::V3D generatePoint(Kernel::PseudoRandomNumberGenerator &rng,
-                            const BoundingBox &activeRegion,
-                            const size_t maxAttempts) const;
+
   bool isValid(const Kernel::V3D &point) const;
   int interceptSurfaces(Track &track) const;
 
diff --git a/Framework/Geometry/inc/MantidGeometry/Objects/CSGObject.h b/Framework/Geometry/inc/MantidGeometry/Objects/CSGObject.h
index 54e8ff63d2cd80fcd939164dad358a1961d2c20e..4eadddaf38b5ba6bb08438fce5457d8585b7bfd0 100644
--- a/Framework/Geometry/inc/MantidGeometry/Objects/CSGObject.h
+++ b/Framework/Geometry/inc/MantidGeometry/Objects/CSGObject.h
@@ -165,11 +165,11 @@ public:
   int getPointInObject(Kernel::V3D &point) const override;
 
   /// Select a random point within the object
-  Kernel::V3D generatePointInObject(Kernel::PseudoRandomNumberGenerator &rng,
-                                    const size_t) const override;
-  Kernel::V3D generatePointInObject(Kernel::PseudoRandomNumberGenerator &rng,
-                                    const BoundingBox &activeRegion,
-                                    const size_t) const override;
+  bool generatePointInObject(Kernel::PseudoRandomNumberGenerator &rng,
+                             const size_t, Kernel::V3D &point) const override;
+  bool generatePointInObject(Kernel::PseudoRandomNumberGenerator &rng,
+                             const BoundingBox &activeRegion, const size_t,
+                             Kernel::V3D &point) const override;
 
   // Rendering member functions
   void draw() const override;
diff --git a/Framework/Geometry/inc/MantidGeometry/Objects/IObject.h b/Framework/Geometry/inc/MantidGeometry/Objects/IObject.h
index e5c87704b503cc905c5490d3838b6b889286c53b..c53e199de145f3bede302e622cf614934a4dd2a6 100644
--- a/Framework/Geometry/inc/MantidGeometry/Objects/IObject.h
+++ b/Framework/Geometry/inc/MantidGeometry/Objects/IObject.h
@@ -69,13 +69,14 @@ public:
   virtual double volume() const = 0;
 
   virtual int getPointInObject(Kernel::V3D &point) const = 0;
-  virtual Kernel::V3D
-  generatePointInObject(Kernel::PseudoRandomNumberGenerator &rng,
-                        const size_t) const = 0;
-  virtual Kernel::V3D
-  generatePointInObject(Kernel::PseudoRandomNumberGenerator &rng,
-                        const BoundingBox &activeRegion,
-                        const size_t) const = 0;
+
+  virtual bool generatePointInObject(Kernel::PseudoRandomNumberGenerator &rng,
+                                     const size_t,
+                                     Kernel::V3D &point) const = 0;
+  virtual bool generatePointInObject(Kernel::PseudoRandomNumberGenerator &rng,
+                                     const BoundingBox &activeRegion,
+                                     const size_t,
+                                     Kernel::V3D &point) const = 0;
 
   virtual detail::ShapeInfo::GeometryShape shape() const = 0;
   virtual const detail::ShapeInfo &shapeInfo() const = 0;
diff --git a/Framework/Geometry/inc/MantidGeometry/Objects/MeshObject.h b/Framework/Geometry/inc/MantidGeometry/Objects/MeshObject.h
index 7440e2d4a33609a3863e0719074d61e9f3869913..55b2c489b47ed1a0f553eb9ce9da820828d540e6 100644
--- a/Framework/Geometry/inc/MantidGeometry/Objects/MeshObject.h
+++ b/Framework/Geometry/inc/MantidGeometry/Objects/MeshObject.h
@@ -110,11 +110,11 @@ public:
   int getPointInObject(Kernel::V3D &point) const override;
 
   /// Select a random point within the object
-  Kernel::V3D generatePointInObject(Kernel::PseudoRandomNumberGenerator &rng,
-                                    const size_t) const override;
-  Kernel::V3D generatePointInObject(Kernel::PseudoRandomNumberGenerator &rng,
-                                    const BoundingBox &activeRegion,
-                                    const size_t) const override;
+  bool generatePointInObject(Kernel::PseudoRandomNumberGenerator &rng,
+                             const size_t, Kernel::V3D &point) const override;
+  bool generatePointInObject(Kernel::PseudoRandomNumberGenerator &rng,
+                             const BoundingBox &activeRegion, const size_t,
+                             Kernel::V3D &point) const override;
 
   // Rendering member functions
   void draw() const override;
diff --git a/Framework/Geometry/inc/MantidGeometry/Objects/MeshObject2D.h b/Framework/Geometry/inc/MantidGeometry/Objects/MeshObject2D.h
index ba80b9e46b83909e315fd098e8520c0a915aaca8..229f591d1594be9a364073062279dac850ea63c6 100644
--- a/Framework/Geometry/inc/MantidGeometry/Objects/MeshObject2D.h
+++ b/Framework/Geometry/inc/MantidGeometry/Objects/MeshObject2D.h
@@ -64,11 +64,12 @@ public:
   void getBoundingBox(double &xmax, double &ymax, double &zmax, double &xmin,
                       double &ymin, double &zmin) const override;
   int getPointInObject(Kernel::V3D &point) const override;
-  Kernel::V3D generatePointInObject(Kernel::PseudoRandomNumberGenerator &rng,
-                                    const size_t) const override;
-  Kernel::V3D generatePointInObject(Kernel::PseudoRandomNumberGenerator &rng,
-                                    const BoundingBox &activeRegion,
-                                    const size_t) const override;
+
+  bool generatePointInObject(Kernel::PseudoRandomNumberGenerator &rng,
+                             const size_t, Kernel::V3D &point) const override;
+  bool generatePointInObject(Kernel::PseudoRandomNumberGenerator &rng,
+                             const BoundingBox &activeRegion, const size_t,
+                             Kernel::V3D &point) const override;
   detail::ShapeInfo::GeometryShape shape() const override;
   const detail::ShapeInfo &shapeInfo() const override;
   void GetObjectGeom(detail::ShapeInfo::GeometryShape &type,
diff --git a/Framework/Geometry/src/Instrument/SampleEnvironment.cpp b/Framework/Geometry/src/Instrument/SampleEnvironment.cpp
index 5cd97c88b8d8eea9016626962275f24bfa0b71e7..c88d0b96d65386ff4c85b7f594df5ea59628c767 100644
--- a/Framework/Geometry/src/Instrument/SampleEnvironment.cpp
+++ b/Framework/Geometry/src/Instrument/SampleEnvironment.cpp
@@ -53,41 +53,6 @@ Geometry::BoundingBox SampleEnvironment::boundingBox() const {
   return box;
 }
 
-/**
- * Generate a random point within one of the environment's components. The
- * method first selects a random component and then selects a random point
- * within that component using Object::generatePointObject
- * @param rng A reference to a PseudoRandomNumberGenerator where
- * nextValue should return a flat random number between 0.0 & 1.0
- * @param maxAttempts The maximum number of attempts at generating a point
- * @return The generated point
- */
-Kernel::V3D
-SampleEnvironment::generatePoint(Kernel::PseudoRandomNumberGenerator &rng,
-                                 const size_t maxAttempts) const {
-  auto componentIndex = rng.nextInt(1, static_cast<int>(nelements())) - 1;
-  return m_components[componentIndex]->generatePointInObject(rng, maxAttempts);
-}
-
-/**
- * Generate a random point within one of the environment's components. The
- * method first selects a random component and then selects a random point
- * within that component using Object::generatePointObject
- * @param rng A reference to a PseudoRandomNumberGenerator where
- * nextValue should return a flat random number between 0.0 & 1.0
- * @param activeRegion Restrict the generated point to be defined by this box
- * @param maxAttempts The maximum number of attempts at generating a point
- * @return The generated point
- */
-Kernel::V3D
-SampleEnvironment::generatePoint(Kernel::PseudoRandomNumberGenerator &rng,
-                                 const Geometry::BoundingBox &activeRegion,
-                                 const size_t maxAttempts) const {
-  auto componentIndex = rng.nextInt(1, static_cast<int>(nelements())) - 1;
-  return m_components[componentIndex]->generatePointInObject(rng, activeRegion,
-                                                             maxAttempts);
-}
-
 /**
  * Is the point given a valid point within the environment
  * @param point Is the point valid within the environment
diff --git a/Framework/Geometry/src/Objects/CSGObject.cpp b/Framework/Geometry/src/Objects/CSGObject.cpp
index b594a994afbfca7cb6f5a5e86eec10b7b73a214f..d6ec6752900cb8577f73212797757e1e1f926dc7 100644
--- a/Framework/Geometry/src/Objects/CSGObject.cpp
+++ b/Framework/Geometry/src/Objects/CSGObject.cpp
@@ -1982,11 +1982,12 @@ int CSGObject::getPointInObject(Kernel::V3D &point) const {
  * @param rng  A reference to a PseudoRandomNumberGenerator where
  * nextValue should return a flat random number between 0.0 & 1.0
  * @param maxAttempts The maximum number of attempts at generating a point
- * @return The generated point
+ * @param point The generated point
+ * @return whether a point was generated in the object or not
  */
-V3D CSGObject::generatePointInObject(PseudoRandomNumberGenerator &rng,
-                                     const size_t maxAttempts) const {
-  V3D point;
+bool CSGObject::generatePointInObject(PseudoRandomNumberGenerator &rng,
+                                      const size_t maxAttempts,
+                                      V3D &point) const {
   // If the shape fills its bounding box well enough then the most efficient
   // way to get the point is just brute force. We'll try that first with
   // just a few attempts.
@@ -2018,14 +2019,14 @@ V3D CSGObject::generatePointInObject(PseudoRandomNumberGenerator &rng,
       maybePoint = RandomPoint::inGenericShape(
           *this, rng, maxAttempts - bruteForceAttempts);
       if (!maybePoint) {
-        throw std::runtime_error("Unable to generate point in object after " +
-                                 std::to_string(maxAttempts) + " attempts");
+        return false;
+      } else {
+        point = *maybePoint;
+        break;
       }
-      point = *maybePoint;
-      break;
     }
   }
-  return point;
+  return true;
 }
 
 /**
@@ -2036,12 +2037,14 @@ V3D CSGObject::generatePointInObject(PseudoRandomNumberGenerator &rng,
  * @param activeRegion Restrict point generation to this sub-region of the
  * object
  * @param maxAttempts The maximum number of attempts at generating a point
- * @return The newly generated point
- */
-V3D CSGObject::generatePointInObject(Kernel::PseudoRandomNumberGenerator &rng,
-                                     const BoundingBox &activeRegion,
-                                     const size_t maxAttempts) const {
-  boost::optional<V3D> point{boost::none};
+ * @param point The generated point
+ * @return whether a point was generated in the object or not
+ */
+bool CSGObject::generatePointInObject(Kernel::PseudoRandomNumberGenerator &rng,
+                                      const BoundingBox &activeRegion,
+                                      const size_t maxAttempts,
+                                      V3D &point) const {
+  boost::optional<V3D> pointGenerated{boost::none};
   // We'll first try brute force. If the shape fills its bounding box
   // well enough, this should be the fastest method.
   // Increasing the brute force attemps speeds up the shapes that fill
@@ -2050,8 +2053,9 @@ V3D CSGObject::generatePointInObject(Kernel::PseudoRandomNumberGenerator &rng,
   // shape, its dimension and orientation.
   const size_t bruteForceAttempts{
       std::min(static_cast<size_t>(5), maxAttempts)};
-  point = RandomPoint::bounded(*this, rng, activeRegion, bruteForceAttempts);
-  if (!point) {
+  pointGenerated =
+      RandomPoint::bounded(*this, rng, activeRegion, bruteForceAttempts);
+  if (!pointGenerated) {
     detail::ShapeInfo::GeometryShape shape;
     std::vector<Kernel::V3D> shapeVectors;
     double radius;
@@ -2060,32 +2064,37 @@ V3D CSGObject::generatePointInObject(Kernel::PseudoRandomNumberGenerator &rng,
     GetObjectGeom(shape, shapeVectors, innerRadius, radius, height);
     switch (shape) {
     case detail::ShapeInfo::GeometryShape::CUBOID:
-      point = RandomPoint::bounded<RandomPoint::inCuboid>(
-          m_handler->shapeInfo(), rng, activeRegion, maxAttempts);
+      pointGenerated = RandomPoint::bounded<RandomPoint::inCuboid>(
+          m_handler->shapeInfo(), rng, activeRegion,
+          maxAttempts - bruteForceAttempts);
       break;
     case detail::ShapeInfo::GeometryShape::CYLINDER:
-      point = RandomPoint::bounded<RandomPoint::inCylinder>(
-          m_handler->shapeInfo(), rng, activeRegion, maxAttempts);
+      pointGenerated = RandomPoint::bounded<RandomPoint::inCylinder>(
+          m_handler->shapeInfo(), rng, activeRegion,
+          maxAttempts - bruteForceAttempts);
       break;
     case detail::ShapeInfo::GeometryShape::HOLLOWCYLINDER:
-      point = RandomPoint::bounded<RandomPoint::inHollowCylinder>(
-          m_handler->shapeInfo(), rng, activeRegion, maxAttempts);
+      pointGenerated = RandomPoint::bounded<RandomPoint::inHollowCylinder>(
+          m_handler->shapeInfo(), rng, activeRegion,
+          maxAttempts - bruteForceAttempts);
       break;
     case detail::ShapeInfo::GeometryShape::SPHERE:
-      point = RandomPoint::bounded<RandomPoint::inSphere>(
-          m_handler->shapeInfo(), rng, activeRegion, maxAttempts);
+      pointGenerated = RandomPoint::bounded<RandomPoint::inSphere>(
+          m_handler->shapeInfo(), rng, activeRegion,
+          maxAttempts - bruteForceAttempts);
       break;
     default:
-      point = RandomPoint::bounded(*this, rng, activeRegion,
-                                   maxAttempts - bruteForceAttempts);
+      pointGenerated = RandomPoint::bounded(*this, rng, activeRegion,
+                                            maxAttempts - bruteForceAttempts);
       break;
     }
   }
-  if (!point) {
-    throw std::runtime_error("Unable to generate point in object after " +
-                             std::to_string(maxAttempts) + " attempts");
+  if (!pointGenerated) {
+    return false;
+  } else {
+    point = *pointGenerated;
+    return true;
   }
-  return *point;
 }
 
 /**
diff --git a/Framework/Geometry/src/Objects/MeshObject.cpp b/Framework/Geometry/src/Objects/MeshObject.cpp
index be1c15a67ee795c9b5cdd47481d66b77fec4796f..150cefb38ca71cbb7a7aafd0109e661b8f4178ef 100644
--- a/Framework/Geometry/src/Objects/MeshObject.cpp
+++ b/Framework/Geometry/src/Objects/MeshObject.cpp
@@ -383,17 +383,18 @@ int MeshObject::getPointInObject(Kernel::V3D &point) const {
  * @param rng  A reference to a PseudoRandomNumberGenerator where
  * nextValue should return a flat random number between 0.0 & 1.0
  * @param maxAttempts The maximum number of attempts at generating a point
- * @return The generated point
+ * @param point The newly generated point
+ * @return Whether a point was generated or not
  */
-Kernel::V3D
-MeshObject::generatePointInObject(Kernel::PseudoRandomNumberGenerator &rng,
-                                  const size_t maxAttempts) const {
+bool MeshObject::generatePointInObject(Kernel::PseudoRandomNumberGenerator &rng,
+                                       const size_t maxAttempts,
+                                       Kernel::V3D &point) const {
   const auto &bbox = getBoundingBox();
   if (bbox.isNull()) {
     throw std::runtime_error("Object::generatePointInObject() - Invalid "
                              "bounding box. Cannot generate new point.");
   }
-  return generatePointInObject(rng, bbox, maxAttempts);
+  return generatePointInObject(rng, bbox, maxAttempts, point);
 }
 
 /**
@@ -404,21 +405,23 @@ MeshObject::generatePointInObject(Kernel::PseudoRandomNumberGenerator &rng,
  * @param activeRegion Restrict point generation to this sub-region of the
  * object
  * @param maxAttempts The maximum number of attempts at generating a point
- * @return The newly generated point
+ * @param point The newly generated point
+ * @return Whether a point was generated or not
  */
-Kernel::V3D
-MeshObject::generatePointInObject(Kernel::PseudoRandomNumberGenerator &rng,
-                                  const BoundingBox &activeRegion,
-                                  const size_t maxAttempts) const {
+bool MeshObject::generatePointInObject(Kernel::PseudoRandomNumberGenerator &rng,
+                                       const BoundingBox &activeRegion,
+                                       const size_t maxAttempts,
+                                       Kernel::V3D &point) const {
 
-  const auto point =
+  const auto pointGenerated =
       RandomPoint::bounded(*this, rng, activeRegion, maxAttempts);
-  if (!point) {
-    throw std::runtime_error("Object::generatePointInObject() - Unable to "
-                             "generate point in object after " +
-                             std::to_string(maxAttempts) + " attempts");
+
+  if (!pointGenerated) {
+    return false;
+  } else {
+    point = *pointGenerated;
+    return true;
   }
-  return *point;
 }
 
 /**
diff --git a/Framework/Geometry/src/Objects/MeshObject2D.cpp b/Framework/Geometry/src/Objects/MeshObject2D.cpp
index 02e14f49bd4fc44f319834f63fc481944681858e..f7ef600ba62fc5fd1357bbc77e77cd597c58c6dd 100644
--- a/Framework/Geometry/src/Objects/MeshObject2D.cpp
+++ b/Framework/Geometry/src/Objects/MeshObject2D.cpp
@@ -388,18 +388,19 @@ int MeshObject2D::getPointInObject(Kernel::V3D &point) const {
   return this->isValid(point) ? 1 : 0;
 }
 
-Kernel::V3D MeshObject2D::generatePointInObject(
-    Kernel::PseudoRandomNumberGenerator & /*rng*/,
-    const size_t /*unused*/) const {
+bool MeshObject2D::generatePointInObject(
+    Kernel::PseudoRandomNumberGenerator & /*rng*/, const size_t /*unused*/,
+    Kernel::V3D & /*unused*/) const {
   // How this would work for a finite plane is not clear. Points within the
   // plane can of course be generated, but most implementations of this method
   // use the bounding box
   throw std::runtime_error("Not implemented.");
 }
 
-Kernel::V3D MeshObject2D::generatePointInObject(
+bool MeshObject2D::generatePointInObject(
     Kernel::PseudoRandomNumberGenerator & /*rng*/,
-    const BoundingBox & /*activeRegion*/, const size_t /*unused*/) const {
+    const BoundingBox & /*activeRegion*/, const size_t /*unused*/,
+    Kernel::V3D & /*unused*/) const {
 
   // How this would work for a finite plane is not clear. Points within the
   // plane can of course be generated, but most implementations of this method
diff --git a/Framework/Geometry/test/CSGObjectTest.h b/Framework/Geometry/test/CSGObjectTest.h
index 4cef5555940d1cbd401ab9af6ff5a5f0a174396b..0ffe03d5a58e6b6244e0934266c0dad13d6de397 100644
--- a/Framework/Geometry/test/CSGObjectTest.h
+++ b/Framework/Geometry/test/CSGObjectTest.h
@@ -731,8 +731,8 @@ public:
     auto shell = ComponentCreationHelper::createHollowShell(0.5, 1.0);
     constexpr size_t maxAttempts{1};
     V3D point;
-    TS_ASSERT_THROWS_NOTHING(
-        point = shell->generatePointInObject(rng, maxAttempts));
+    TS_ASSERT_EQUALS(shell->generatePointInObject(rng, maxAttempts, point),
+                     true);
 
     constexpr double tolerance{1e-10};
     TS_ASSERT_DELTA(-1. + 2. * 0.55, point.X(), tolerance);
@@ -760,8 +760,8 @@ public:
         ComponentCreationHelper::createCuboid(xLength, yLength, zLength);
     constexpr size_t maxAttempts{0};
     V3D point;
-    TS_ASSERT_THROWS_NOTHING(
-        point = cuboid->generatePointInObject(rng, maxAttempts));
+    TS_ASSERT_EQUALS(cuboid->generatePointInObject(rng, maxAttempts, point),
+                     true);
 
     constexpr double tolerance{1e-10};
     TS_ASSERT_DELTA(xLength - randX * 2. * xLength, point.X(), tolerance);
@@ -794,8 +794,8 @@ public:
         radius, height, bottomCentre, axis, "cyl");
     constexpr size_t maxAttempts{0};
     V3D point;
-    TS_ASSERT_THROWS_NOTHING(
-        point = cylinder->generatePointInObject(rng, maxAttempts));
+    TS_ASSERT_EQUALS(cylinder->generatePointInObject(rng, maxAttempts, point),
+                     true);
     // Global->cylinder local coordinates
     point -= bottomCentre;
     constexpr double tolerance{1e-10};
@@ -833,8 +833,8 @@ public:
         innerRadius, radius, height, bottomCentre, axis, "hol-cyl");
     constexpr size_t maxAttempts{0};
     V3D point;
-    TS_ASSERT_THROWS_NOTHING(
-        point = hollowCylinder->generatePointInObject(rng, maxAttempts));
+    TS_ASSERT_EQUALS(
+        hollowCylinder->generatePointInObject(rng, maxAttempts, point), true);
     // Global->cylinder local coordinates
     point -= bottomCentre;
     constexpr double tolerance{1e-10};
@@ -865,8 +865,8 @@ public:
     auto sphere = ComponentCreationHelper::createSphere(radius);
     constexpr size_t maxAttempts{0};
     V3D point;
-    TS_ASSERT_THROWS_NOTHING(
-        point = sphere->generatePointInObject(rng, maxAttempts));
+    TS_ASSERT_EQUALS(sphere->generatePointInObject(rng, maxAttempts, point),
+                     true);
     // Global->cylinder local coordinates
     constexpr double tolerance{1e-10};
     const double azimuthalAngle{2. * M_PI * randT};
@@ -893,8 +893,9 @@ public:
     // inside hole
     auto shell = ComponentCreationHelper::createHollowShell(0.5, 1.0);
     constexpr size_t maxAttempts{1};
-    TS_ASSERT_THROWS(shell->generatePointInObject(rng, maxAttempts),
-                     const std::runtime_error &);
+    V3D point;
+    TS_ASSERT_EQUALS(shell->generatePointInObject(rng, maxAttempts, point),
+                     false);
   }
 
   void testGeneratePointInsideRespectsActiveRegion() {
@@ -916,8 +917,9 @@ public:
     BoundingBox activeRegion(0.1, 0.1, 0.1, -0.1, -0.1, -0.1);
     constexpr size_t maxAttempts{1};
     V3D point;
-    TS_ASSERT_THROWS_NOTHING(
-        point = ball->generatePointInObject(rng, activeRegion, maxAttempts));
+    TS_ASSERT_EQUALS(
+        ball->generatePointInObject(rng, activeRegion, maxAttempts, point),
+        true);
     // We should get the point generated from the second 'random' triplet.
     constexpr double tolerance{1e-10};
     TS_ASSERT_DELTA(-0.1 + randX * 0.2, point.X(), tolerance)
@@ -1745,44 +1747,52 @@ public:
 
   void test_generatePointInside_Cuboid_With_ActiveRegion() {
     constexpr size_t maxAttempts{500};
+    V3D point;
     for (size_t i{0}; i < m_npoints; ++i) {
-      m_cuboid->generatePointInObject(m_rng, m_activeRegion, maxAttempts);
+      m_cuboid->generatePointInObject(m_rng, m_activeRegion, maxAttempts,
+                                      point);
     }
   }
 
   void test_generatePointInside_Cylinder_With_ActiveRegion() {
     constexpr size_t maxAttempts{500};
+    V3D point;
     for (size_t i{0}; i < m_npoints; ++i) {
-      m_cylinder->generatePointInObject(m_rng, m_activeRegion, maxAttempts);
+      m_cylinder->generatePointInObject(m_rng, m_activeRegion, maxAttempts,
+                                        point);
     }
   }
 
   void test_generatePointInside_Rotated_Cuboid() {
     constexpr size_t maxAttempts{500};
+    V3D point;
     for (size_t i = 0; i < m_npoints; ++i) {
-      m_rotatedCuboid->generatePointInObject(m_rng, maxAttempts);
+      m_rotatedCuboid->generatePointInObject(m_rng, maxAttempts, point);
     }
   }
 
   void test_generatePointInside_Rotated_Cuboid_With_ActiveRegion() {
     constexpr size_t maxAttempts{500};
+    V3D point;
     for (size_t i = 0; i < m_npoints; ++i) {
-      m_rotatedCuboid->generatePointInObject(m_rng, m_activeRegion,
-                                             maxAttempts);
+      m_rotatedCuboid->generatePointInObject(m_rng, m_activeRegion, maxAttempts,
+                                             point);
     }
   }
 
   void test_generatePointInside_Sphere() {
     constexpr size_t maxAttempts{500};
+    V3D point;
     for (size_t i = 0; i < m_npoints; ++i) {
-      m_sphere->generatePointInObject(m_rng, maxAttempts);
+      m_sphere->generatePointInObject(m_rng, maxAttempts, point);
     }
   }
 
   void test_generatePointInside_sphericalShell() {
     constexpr size_t maxAttempts{500};
+    V3D point;
     for (size_t i = 0; i < m_npoints; ++i) {
-      m_sphericalShell->generatePointInObject(m_rng, maxAttempts);
+      m_sphericalShell->generatePointInObject(m_rng, maxAttempts, point);
     }
   }
 
diff --git a/Framework/Geometry/test/MeshObject2DTest.h b/Framework/Geometry/test/MeshObject2DTest.h
index 5bf59706f67c0b49983837acf20564f5d8f9f6e8..9eabaf8328e70d9ffb3558966e9663847d0e758a 100644
--- a/Framework/Geometry/test/MeshObject2DTest.h
+++ b/Framework/Geometry/test/MeshObject2DTest.h
@@ -332,7 +332,8 @@ public:
     // plane
     auto mesh = makeSimpleTriangleMesh();
     testing::NiceMock<MockRNG> generator;
-    TS_ASSERT_THROWS(mesh.generatePointInObject(generator, 10),
+    V3D point;
+    TS_ASSERT_THROWS(mesh.generatePointInObject(generator, 10, point),
                      std::runtime_error &);
   }
   // Characterisation test.
@@ -343,8 +344,10 @@ public:
     auto mesh = makeSimpleTriangleMesh();
     testing::NiceMock<MockRNG> generator;
     Mantid::Geometry::BoundingBox boundingBox;
-    TS_ASSERT_THROWS(mesh.generatePointInObject(generator, boundingBox, 10),
-                     std::runtime_error &);
+    V3D point;
+    TS_ASSERT_THROWS(
+        mesh.generatePointInObject(generator, boundingBox, 10, point),
+        std::runtime_error &);
   }
 
   // Characterisation test.
diff --git a/Framework/Geometry/test/MeshObjectTest.h b/Framework/Geometry/test/MeshObjectTest.h
index 160da94b438fcfbeffbf77880269ae09f0008d4b..b18e51afc8199926e5ebad936dd72696ce84f2f6 100644
--- a/Framework/Geometry/test/MeshObjectTest.h
+++ b/Framework/Geometry/test/MeshObjectTest.h
@@ -848,7 +848,7 @@ public:
     size_t maxAttempts(1);
     V3D point;
     TS_ASSERT_THROWS_NOTHING(
-        point = geom_obj->generatePointInObject(rng, maxAttempts));
+        geom_obj->generatePointInObject(rng, maxAttempts, point));
 
     const double tolerance(1e-10);
     TS_ASSERT_DELTA(0.90, point.X(), tolerance);
@@ -870,8 +870,9 @@ public:
     //  which is outside the octahedron
     auto geom_obj = createOctahedron();
     size_t maxAttempts(1);
-    TS_ASSERT_THROWS(geom_obj->generatePointInObject(rng, maxAttempts),
-                     const std::runtime_error &);
+    V3D point;
+    TS_ASSERT_EQUALS(geom_obj->generatePointInObject(rng, maxAttempts, point),
+                     false);
   }
 
   void testVolumeOfCube() {
@@ -1086,16 +1087,18 @@ public:
   void test_generatePointInside_Convex_Solid() {
     const size_t npoints(6000);
     const size_t maxAttempts(500);
+    V3D point;
     for (size_t i = 0; i < npoints; ++i) {
-      octahedron->generatePointInObject(rng, maxAttempts);
+      octahedron->generatePointInObject(rng, maxAttempts, point);
     }
   }
 
   void test_generatePointInside_NonConvex_Solid() {
     const size_t npoints(6000);
     const size_t maxAttempts(500);
+    V3D point;
     for (size_t i = 0; i < npoints; ++i) {
-      lShape->generatePointInObject(rng, maxAttempts);
+      lShape->generatePointInObject(rng, maxAttempts, point);
     }
   }
 
diff --git a/Framework/Geometry/test/SampleEnvironmentTest.h b/Framework/Geometry/test/SampleEnvironmentTest.h
index 0586b9ceb8e5dab0dc65602eca2b572ee9f4e7a9..dcbaa942ab3605126f7ae03f661c1f4a52e1337a 100644
--- a/Framework/Geometry/test/SampleEnvironmentTest.h
+++ b/Framework/Geometry/test/SampleEnvironmentTest.h
@@ -91,79 +91,6 @@ public:
     TS_ASSERT_DELTA(0.2, widths.Z(), 1e-12);
   }
 
-  void testGeneratePointConsidersAllComponents() {
-    using namespace ::testing;
-    using namespace Mantid::Kernel;
-
-    auto kit = createTestKit();
-    size_t maxAttempts(1);
-
-    // Generate "random" sequence
-    MockRNG rng;
-    // Selects first component
-    EXPECT_CALL(rng, nextInt(_, _)).Times(Exactly(1)).WillOnce(Return(1));
-    EXPECT_CALL(rng, nextValue())
-        .Times(3)
-        .WillOnce(Return(0.55))
-        .WillOnce(Return(0.65))
-        .WillOnce(Return(0.70));
-    V3D comp1Point;
-    TS_ASSERT_THROWS_NOTHING(comp1Point = kit->generatePoint(rng, maxAttempts));
-    TS_ASSERT(kit->isValid(comp1Point));
-    Mock::VerifyAndClearExpectations(&rng);
-
-    // Selects second component
-    EXPECT_CALL(rng, nextInt(_, _)).Times(Exactly(1)).WillOnce(Return(2));
-    EXPECT_CALL(rng, nextValue())
-        .Times(3)
-        .WillOnce(Return(0.55))
-        .WillOnce(Return(0.65))
-        .WillOnce(Return(0.70));
-    V3D comp2Point;
-    TS_ASSERT_THROWS_NOTHING(comp2Point = kit->generatePoint(rng, maxAttempts));
-    TS_ASSERT(comp2Point != comp1Point);
-    TS_ASSERT(kit->isValid(comp2Point));
-    Mock::VerifyAndClearExpectations(&rng);
-
-    // Selects third component
-    EXPECT_CALL(rng, nextInt(_, _)).Times(Exactly(1)).WillOnce(Return(3));
-    EXPECT_CALL(rng, nextValue())
-        .Times(3)
-        .WillOnce(Return(0.55))
-        .WillOnce(Return(0.65))
-        .WillOnce(Return(0.70));
-    V3D comp3Point;
-    TS_ASSERT_THROWS_NOTHING(comp3Point = kit->generatePoint(rng, maxAttempts));
-    TS_ASSERT(comp3Point != comp2Point);
-    TS_ASSERT(comp3Point != comp1Point);
-    TS_ASSERT(kit->isValid(comp3Point));
-    Mock::VerifyAndClearExpectations(&rng);
-  }
-
-  void testGeneratePointRespectsActiveRegion() {
-    using namespace ::testing;
-    using namespace Mantid::Kernel;
-
-    auto kit = createTestKit();
-    size_t maxAttempts(1);
-
-    // Generate "random" sequence
-    MockRNG rng;
-    // Sequence will try to select one of the non-container pieces
-    EXPECT_CALL(rng, nextInt(_, _)).Times(Exactly(1)).WillOnce(Return(2));
-    EXPECT_CALL(rng, nextValue())
-        .Times(3)
-        .WillOnce(Return(0.5))
-        .WillOnce(Return(0.5))
-        .WillOnce(Return(0.5));
-    // Restrict region to can
-    TS_ASSERT_THROWS(kit->generatePoint(rng,
-                                        kit->getContainer().getBoundingBox(),
-                                        maxAttempts),
-                     const std::runtime_error &);
-    Mock::VerifyAndClearExpectations(&rng);
-  }
-
 private:
   boost::shared_ptr<SampleEnvironment> createTestKit() {
     using namespace Mantid::Geometry;
diff --git a/Framework/PythonInterface/plugins/algorithms/WorkflowAlgorithms/DirectILLSelfShielding.py b/Framework/PythonInterface/plugins/algorithms/WorkflowAlgorithms/DirectILLSelfShielding.py
index 75f82c7026fbca519389ca84b44509adf5fb2ce9..1df2865f50bac58e73a74dff8aa5e568c1caea93 100644
--- a/Framework/PythonInterface/plugins/algorithms/WorkflowAlgorithms/DirectILLSelfShielding.py
+++ b/Framework/PythonInterface/plugins/algorithms/WorkflowAlgorithms/DirectILLSelfShielding.py
@@ -62,6 +62,7 @@ class DirectILLSelfShielding(DataProcessorAlgorithm):
     def PyInit(self):
         """Initialize the algorithm's input and output properties."""
         PROPGROUP_SIMULATION_INSTRUMENT = 'Simulation Instrument Settings'
+        positiveInt = IntBoundedValidator(lower=1)
         greaterThanOneInt = IntBoundedValidator(lower=2)
         greaterThanTwoInt = IntBoundedValidator(lower=3)
         inputWorkspaceValidator = CompositeValidator()
@@ -124,6 +125,11 @@ class DirectILLSelfShielding(DataProcessorAlgorithm):
                              validator=greaterThanTwoInt,
                              direction=Direction.Input,
                              doc='Number of wavelength points where the simulation is performed (default: all).')
+        self.declareProperty(name=common.PROP_EVENTS_PER_WAVELENGTH,
+                             defaultValue=300,
+                             validator=positiveInt,
+                             direction=Direction.Input,
+                             doc='The number of neutron events to simulate per wavelength point.')
 
     def validateInputs(self):
         """Check for issues with user input."""
@@ -150,6 +156,7 @@ class DirectILLSelfShielding(DataProcessorAlgorithm):
                                     EMode='Direct',
                                     EnableLogging=self._subalgLogging)
         wavelengthPoints = self.getProperty(common.PROP_NUMBER_OF_SIMULATION_WAVELENGTHS).value
+        eventsPerPoint =self.getProperty(common.PROP_EVENTS_PER_WAVELENGTH).value
         correctionWSName = self._names.withSuffix('correction')
         useFullInstrument = self.getProperty(common.PROP_SIMULATION_INSTRUMENT).value == common.SIMULATION_INSTRUMENT_FULL
         if useFullInstrument:
@@ -158,7 +165,9 @@ class DirectILLSelfShielding(DataProcessorAlgorithm):
                                                 SparseInstrument=False,
                                                 NumberOfWavelengthPoints=wavelengthPoints,
                                                 Interpolation='CSpline',
-                                                EnableLogging=self._subalgLogging)
+                                                EnableLogging=self._subalgLogging,
+                                                EventsPerPoint=eventsPerPoint,
+                                                SeedValue=common.SIMULATION_SEED)
         else:
             rows = self.getProperty(common.PROP_SPARSE_INSTRUMENT_ROWS).value
             columns = self.getProperty(common.PROP_SPARSE_INSTRUMENT_COLUMNS).value
@@ -169,7 +178,9 @@ class DirectILLSelfShielding(DataProcessorAlgorithm):
                                                 NumberOfDetectorColumns=columns,
                                                 NumberOfWavelengthPoints=wavelengthPoints,
                                                 Interpolation='CSpline',
-                                                EnableLogging=self._subalgLogging)
+                                                EnableLogging=self._subalgLogging,
+                                                EventsPerPoint=eventsPerPoint,
+                                                SeedValue=common.SIMULATION_SEED)
         self._cleanup.cleanup(wavelengthWS)
         correctionWS = ConvertUnits(InputWorkspace=correctionWS,
                                     OutputWorkspace=correctionWSName,
diff --git a/Framework/PythonInterface/plugins/algorithms/WorkflowAlgorithms/DirectILL_common.py b/Framework/PythonInterface/plugins/algorithms/WorkflowAlgorithms/DirectILL_common.py
index ebc921a82721c4d952b00e0125de3db6a9bb9ff8..a28464a248c92b7ea6f91067c91b915ef1a53eb0 100644
--- a/Framework/PythonInterface/plugins/algorithms/WorkflowAlgorithms/DirectILL_common.py
+++ b/Framework/PythonInterface/plugins/algorithms/WorkflowAlgorithms/DirectILL_common.py
@@ -91,6 +91,7 @@ PROP_MON_INDEX = 'Monitor'
 PROP_MON_PEAK_SIGMA_MULTIPLIER = 'MonitorPeakWidthInSigmas'
 PROP_NORMALISATION = 'Normalisation'
 PROP_NUMBER_OF_SIMULATION_WAVELENGTHS = 'NumberOfSimulatedWavelengths'
+PROP_EVENTS_PER_WAVELENGTH = 'EventsPerPoint'
 PROP_OUTPUT_DET_EPP_WS = 'OutputEPPWorkspace'
 PROP_OUTPUT_DIAGNOSTICS_REPORT = 'OutputReport'
 PROP_OUTPUT_DIAGNOSTICS_REPORT_WS = 'OutputReportWorkspace'
@@ -123,6 +124,8 @@ PROPGROUP_OPTIONAL_OUTPUT = 'Optional Output'
 SIMULATION_INSTRUMENT_FULL = 'Full Instrument'
 SIMULATION_INSTRUMEN_SPARSE = 'Sparse Instrument'
 
+SIMULATION_SEED = 123456789
+
 SUBALG_LOGGING_OFF = 'Logging OFF'
 SUBALG_LOGGING_ON = 'Logging ON'
 
diff --git a/Testing/SystemTests/tests/analysis/ILLDirectGeometryReductionTest.py b/Testing/SystemTests/tests/analysis/ILLDirectGeometryReductionTest.py
index 15c6d91942ed5c27f6880bb4099144778dfd594d..82c6a6a6e9f0d12bcb99f6218e1b8b0583d9e508 100644
--- a/Testing/SystemTests/tests/analysis/ILLDirectGeometryReductionTest.py
+++ b/Testing/SystemTests/tests/analysis/ILLDirectGeometryReductionTest.py
@@ -58,7 +58,8 @@ class IN4(systemtesting.MantidSystemTest):
         DirectILLSelfShielding(
             InputWorkspace='sample',
             OutputWorkspace='corrections',
-            NumberOfSimulatedWavelengths=10)
+            NumberOfSimulatedWavelengths=10,
+            EventsPerPoint=5000)
         DirectILLApplySelfShielding(
             InputWorkspace='sample',
             OutputWorkspace='corrected',
diff --git a/Testing/SystemTests/tests/analysis/reference/ILL_IN4_SofQW.nxs.md5 b/Testing/SystemTests/tests/analysis/reference/ILL_IN4_SofQW.nxs.md5
index 8d279153ff2d013672c7bbb5c451999d82b4ee75..c27e1dd131d689162bfd987e4978a4cb89e2193b 100644
--- a/Testing/SystemTests/tests/analysis/reference/ILL_IN4_SofQW.nxs.md5
+++ b/Testing/SystemTests/tests/analysis/reference/ILL_IN4_SofQW.nxs.md5
@@ -1 +1 @@
-f074b58b9a4f1afbd81d1ec830e13d53
+23f425fa47b9f5ae9820ffec9b00608f
diff --git a/Testing/SystemTests/tests/analysis/reference/complexEnvironmentAbsorb.nxs.md5 b/Testing/SystemTests/tests/analysis/reference/complexEnvironmentAbsorb.nxs.md5
index a4fc7f14620779e6ee8d4f71c875a828a16f92d6..2327a6cdd059d7b128738c52e426826edf39ece5 100644
--- a/Testing/SystemTests/tests/analysis/reference/complexEnvironmentAbsorb.nxs.md5
+++ b/Testing/SystemTests/tests/analysis/reference/complexEnvironmentAbsorb.nxs.md5
@@ -1 +1 @@
-4de4408058b77acae0714f993f33a09e
+6b853d4629bcccf65617d9b3cdea977d
diff --git a/Testing/SystemTests/tests/analysis/reference/complexEnvironmentMultiPartRotated.nxs.md5 b/Testing/SystemTests/tests/analysis/reference/complexEnvironmentMultiPartRotated.nxs.md5
index f9ea4cd1fffc5a6f57db7b5ff3ab542b563e5fa0..42c26c0e4400dbd8a0f7459168c7c9eff65bd24f 100644
--- a/Testing/SystemTests/tests/analysis/reference/complexEnvironmentMultiPartRotated.nxs.md5
+++ b/Testing/SystemTests/tests/analysis/reference/complexEnvironmentMultiPartRotated.nxs.md5
@@ -1 +1 @@
-fa0d2952a5f6ccc5bf702d51e1e40bb1
+a4e917f964925b19eb09c8e8ea8fecda
diff --git a/Testing/SystemTests/tests/analysis/reference/complexEnvironmentMultiPartTranslated.nxs.md5 b/Testing/SystemTests/tests/analysis/reference/complexEnvironmentMultiPartTranslated.nxs.md5
index 9f64039b0fbf472a4afbc60635bb5910fd7fbf98..37df1c7408a942cdedd3745c977fc5601a5580a2 100644
--- a/Testing/SystemTests/tests/analysis/reference/complexEnvironmentMultiPartTranslated.nxs.md5
+++ b/Testing/SystemTests/tests/analysis/reference/complexEnvironmentMultiPartTranslated.nxs.md5
@@ -1 +1 @@
-b13aac4bcc6d13ad0bd7e3bb9f429d14
+a2a019e073654a7145570dc786386f07
diff --git a/Testing/SystemTests/tests/analysis/reference/complexEnvironmentRotatedAbsorb.nxs.md5 b/Testing/SystemTests/tests/analysis/reference/complexEnvironmentRotatedAbsorb.nxs.md5
index 8f1e7c6ea74dd3b17a46101d24b1de3ecdd05455..a8e70ea1d927558cf53bf246125dc15995a24383 100644
--- a/Testing/SystemTests/tests/analysis/reference/complexEnvironmentRotatedAbsorb.nxs.md5
+++ b/Testing/SystemTests/tests/analysis/reference/complexEnvironmentRotatedAbsorb.nxs.md5
@@ -1 +1 @@
-3576182125013d7aaef52f3255a9af3b
+56ee53abdc27778ca604e266406ac731
diff --git a/Testing/SystemTests/tests/analysis/reference/complexEnvironmentRotatedTranslatedAbsorb.nxs.md5 b/Testing/SystemTests/tests/analysis/reference/complexEnvironmentRotatedTranslatedAbsorb.nxs.md5
index 7150cd802b19261aae6dcf0681b917ed9d13eb14..686f55074576c0d6e427afe88721e120ef827fcf 100644
--- a/Testing/SystemTests/tests/analysis/reference/complexEnvironmentRotatedTranslatedAbsorb.nxs.md5
+++ b/Testing/SystemTests/tests/analysis/reference/complexEnvironmentRotatedTranslatedAbsorb.nxs.md5
@@ -1 +1 @@
-d5a6d85ec779d813c5642581f22e6f43
+db7347e9eecc708c6ebbea4168f106cc
diff --git a/Testing/SystemTests/tests/analysis/reference/complexEnvironmentTranslatedAbsorb.nxs.md5 b/Testing/SystemTests/tests/analysis/reference/complexEnvironmentTranslatedAbsorb.nxs.md5
index a443c5f68460e435e3d4f55e94806d39f223bfc0..7970529cc441cf0c64b72b94e7fe15639a03d384 100644
--- a/Testing/SystemTests/tests/analysis/reference/complexEnvironmentTranslatedAbsorb.nxs.md5
+++ b/Testing/SystemTests/tests/analysis/reference/complexEnvironmentTranslatedAbsorb.nxs.md5
@@ -1 +1 @@
-148984e65e9a5b6e0422c43814129190
+82c237b57482e67b10808fd62f33164a
diff --git a/Testing/SystemTests/tests/analysis/reference/complexShapeRotatedAbsorbEnvironment.nxs.md5 b/Testing/SystemTests/tests/analysis/reference/complexShapeRotatedAbsorbEnvironment.nxs.md5
index 116a73c4ad1222a784f9a6386ddbdafcb22eb64d..065b8e49b47aa45c7b291acd68d3ca42474086b2 100644
--- a/Testing/SystemTests/tests/analysis/reference/complexShapeRotatedAbsorbEnvironment.nxs.md5
+++ b/Testing/SystemTests/tests/analysis/reference/complexShapeRotatedAbsorbEnvironment.nxs.md5
@@ -1 +1 @@
-f1c6070ba64234c00d1ce12f0b218f9a
+9013e0e002bb170bb5af7f726ba402c2
diff --git a/docs/source/release/v4.3.0/framework.rst b/docs/source/release/v4.3.0/framework.rst
index 2f8d5dff4a12d8aab440dc47ece5f5dd01d4253b..f98e767f9251833ecd42466f70328c12d7c504b5 100644
--- a/docs/source/release/v4.3.0/framework.rst
+++ b/docs/source/release/v4.3.0/framework.rst
@@ -32,6 +32,7 @@ Improvements
 - :ref:`SaveAscii <algm-SaveAscii>` can now save table workspaces, and :ref:`LoadAscii <algm-LoadAscii>` can load them again.
 - :ref:`TotScatCalculateSelfScattering <algm-TotScatCalculateSelfScattering>` will calculate a normalized self scattering correction for foccues total scattering data.
 - :ref:`MatchAndMergeWorkspaces <algm-MatchAndMergeWorkspaces>` will merge workspaces in a workspace group withing weighting from a set of limits for each workspace and using `MatchSpectra <algm-MatchSpectra>`.
+- :ref:`MonteCarloAbsorption <algm-MonteCarloAbsorption>` Sampling of scattering points during MC simulation now takes into account relative volume of sample and environment components.
 
 
 Data Objects
diff --git a/scripts/test/AbsorptionShapesTest.py b/scripts/test/AbsorptionShapesTest.py
index 67d6370bb0d217a7422ba944c9925ad8e9052480..9d88ede309ae8c68507e2a9867d7820886221a33 100644
--- a/scripts/test/AbsorptionShapesTest.py
+++ b/scripts/test/AbsorptionShapesTest.py
@@ -85,7 +85,7 @@ class AdsorbtionShapesTest(unittest.TestCase):
         mccor_ws,mc_corr = ash.correct_absorption(test_ws,NumberOfWavelengthPoints=20)
         n_bins = mc_corr.blocksize()
         mccorr_ranges = [n_bins,mc_corr.readY(0)[0],mc_corr.readY(0)[n_bins-1]]
-        np.testing.assert_almost_equal(mccorr_ranges ,[97,0.2657,0.0271],4)
+        np.testing.assert_almost_equal(mccorr_ranges ,[97,0.266,0.026],3)
 
     def test_MARI_axis_cylinder(self):
         """ Test that default axis for MARI is different"""
@@ -140,7 +140,7 @@ class AdsorbtionShapesTest(unittest.TestCase):
         mccor_ws,mc_corr = ash.correct_absorption(test_ws,is_mc=True,NumberOfWavelengthPoints=20)
         n_bins = mc_corr.blocksize()
         mccorr_ranges = [n_bins,mc_corr.readY(0)[0],mc_corr.readY(0)[n_bins-1]]
-        np.testing.assert_almost_equal(mccorr_ranges ,[97,0.5253,0.1296],4)
+        np.testing.assert_almost_equal(mccorr_ranges ,[97,0.52,0.13],2)
 
 
     def test_adsrp_hollow_cylinder(self):
@@ -178,7 +178,7 @@ class AdsorbtionShapesTest(unittest.TestCase):
         mccor_ws,mc_corr = ash.correct_absorption(test_ws,is_mc=True,NumberOfWavelengthPoints=20)
         n_bins = mc_corr.blocksize()
         mccorr_ranges = [n_bins,mc_corr.readY(0)[0],mc_corr.readY(0)[n_bins-1]]
-        np.testing.assert_almost_equal(mccorr_ranges ,[97,0.2657,0.0303],4)
+        np.testing.assert_almost_equal(mccorr_ranges ,[97,0.27,0.03],2)
     #
     def test_string_conversion(self):
         """ check if shape conversion to string representation works"""
@@ -235,7 +235,7 @@ class AdsorbtionShapesTest(unittest.TestCase):
         mccor_ws,mc_corr = ash.correct_absorption(test_ws,is_mc=True,NumberOfWavelengthPoints=20)
         n_bins = mc_corr.blocksize()
         mccorr_ranges = [n_bins,mc_corr.readY(0)[0],mc_corr.readY(0)[n_bins-1]]
-        np.testing.assert_almost_equal(mccorr_ranges ,[97,0.6645,0.5098],4)
+        np.testing.assert_almost_equal(mccorr_ranges ,[97,0.66,0.51],2)