From e2af417f67246bc39f5c4513fccc9e1b4095eebc Mon Sep 17 00:00:00 2001
From: Danny Hindson <danny.hindson@stfc.ac.uk>
Date: Wed, 15 Jan 2020 12:23:24 +0000
Subject: [PATCH] Improve sampling of scatter pts in absorption correction
 simulation

The sampling is now performed according to the volume of each part of
the sample\environment that intersects the beam profile

Following changes have been made:

a) the function IObject::generatePointInObject has been modified so that it returns
false rather than raising an exception if it fails to generate a point inside
the object that is also in the active region. This is more efficient when it is
being called with a maxAttempts value of 1 in an attempt to fairly sample the
scatter points among the sample + environment components. This involved a change
in IObject and various child classes (MeshObject, CSGObject, MeshObject2D, Container)

b) the code that calls IObject::generatePointInObject and cycles through the
various parts of the environment\sample has been moved from SampleEnvironment into
MCInteractionVolume so that the sample can be included. There is a new function
MCInteractionVolume::generatePoint that randomly generates a scatter point across the
sample and environment components. Possible the SampleEnvironment class could be retired
entirely and just replaced by a vector of IObject items attached to the sample.

c) change CSGObject::generatePointInObject to stop calling the fallback method when the
maxAttempts parameter equals 1. The fallback method always returns a point if the object's
bounding box is inside the active area which doesn't produce the required sampling across
the env components. This change has modified the random number sequences used in various
tests - including the "sample only" tests

d) added some logging to the simulation to show where the scatter points
occurred. These show that for Pearl around 6% of the scatter points are
in the sample which is less than the 50% assumption previously in the code

e) Several changes to the unit tests (MCInteractionVolumeTest.h, MonteCarloAbsorptionTest.h,
DirectILLSelfShieldingTest.py). The updated sampling means that the absorption corrections
are slightly different than before for cases with a sample + environment.

For DirectILLSelfShieldingTest, an extra parameter has been added to underlying algorithm
(DirectILLSelfShielding) so that this test can continue to use 300 events per point while
the ILLDirectGeometryReductionTest can use 5000 events per point

f) some changes to system tests (ILLDirectGeometryReductionTest, complexShapeAbsorptions)

The calculation used in ILLDirectGeometryReductionTest wasn't converged (changing the seed
gave a ~25% change in the output) so have increased the number of events per point from 300
to 5000. I didn't increase the number of events further because I didn't want to make the
runtime of the system test (esp in debug mode) too large
---
 .../SampleCorrections/MCAbsorptionStrategy.h  |   4 +-
 .../SampleCorrections/MCInteractionVolume.h   |   4 +
 .../Algorithms/src/MonteCarloAbsorption.cpp   |   4 +-
 .../MCAbsorptionStrategy.cpp                  |  11 +-
 .../SampleCorrections/MCInteractionVolume.cpp |  93 ++++++++++++--
 .../test/MCAbsorptionStrategyTest.h           |   9 +-
 .../Algorithms/test/MCInteractionVolumeTest.h | 116 ++++++++++++++++--
 .../test/MonteCarloAbsorptionTest.h           | 109 ++++++++--------
 Framework/Algorithms/test/MonteCarloTesting.h |  18 +++
 .../inc/MantidGeometry/Instrument/Container.h |  17 +--
 .../Instrument/SampleEnvironment.h            |   9 +-
 .../inc/MantidGeometry/Objects/CSGObject.h    |  10 +-
 .../inc/MantidGeometry/Objects/IObject.h      |  15 +--
 .../inc/MantidGeometry/Objects/MeshObject.h   |  10 +-
 .../inc/MantidGeometry/Objects/MeshObject2D.h |  11 +-
 .../src/Instrument/SampleEnvironment.cpp      |  35 ------
 Framework/Geometry/src/Objects/CSGObject.cpp  |  71 ++++++-----
 Framework/Geometry/src/Objects/MeshObject.cpp |  35 +++---
 .../Geometry/src/Objects/MeshObject2D.cpp     |  11 +-
 Framework/Geometry/test/CSGObjectTest.h       |  52 ++++----
 Framework/Geometry/test/MeshObject2DTest.h    |   9 +-
 Framework/Geometry/test/MeshObjectTest.h      |  13 +-
 .../Geometry/test/SampleEnvironmentTest.h     |  73 -----------
 .../DirectILLSelfShielding.py                 |  15 ++-
 .../WorkflowAlgorithms/DirectILL_common.py    |   3 +
 .../ILLDirectGeometryReductionTest.py         |   3 +-
 .../analysis/reference/ILL_IN4_SofQW.nxs.md5  |   2 +-
 .../complexEnvironmentAbsorb.nxs.md5          |   2 +-
 ...complexEnvironmentMultiPartRotated.nxs.md5 |   2 +-
 ...plexEnvironmentMultiPartTranslated.nxs.md5 |   2 +-
 .../complexEnvironmentRotatedAbsorb.nxs.md5   |   2 +-
 ...EnvironmentRotatedTranslatedAbsorb.nxs.md5 |   2 +-
 ...complexEnvironmentTranslatedAbsorb.nxs.md5 |   2 +-
 ...mplexShapeRotatedAbsorbEnvironment.nxs.md5 |   2 +-
 docs/source/release/v4.3.0/framework.rst      |   1 +
 scripts/test/AbsorptionShapesTest.py          |   8 +-
 36 files changed, 459 insertions(+), 326 deletions(-)

diff --git a/Framework/Algorithms/inc/MantidAlgorithms/SampleCorrections/MCAbsorptionStrategy.h b/Framework/Algorithms/inc/MantidAlgorithms/SampleCorrections/MCAbsorptionStrategy.h
index 4c97318f80d..5d1bc89e2d9 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 c761b8740c8..a93d3ad1b61 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 777a569ad2a..ea990c3c901 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 4ead5a72d5c..8858f1e9b69 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 030e13b175a..62d8c874eb0 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 c23ed83a607..31fcefa7323 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 fc78ba3f16b..ccf4ce62e9d 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 cad25ce2912..cf9b4c63476 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 4877acb2b9d..c6ee52715ef 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 ee05a300adf..a18ef6de79a 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 bb47cf66788..5198a15ca40 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 54e8ff63d2c..4eadddaf38b 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 e5c87704b50..c53e199de14 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 7440e2d4a33..55b2c489b47 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 ba80b9e46b8..229f591d159 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 5cd97c88b8d..c88d0b96d65 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 b594a994afb..d6ec6752900 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 be1c15a67ee..150cefb38ca 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 02e14f49bd4..f7ef600ba62 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 4cef5555940..0ffe03d5a58 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 5bf59706f67..9eabaf8328e 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 160da94b438..b18e51afc81 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 0586b9ceb8e..dcbaa942ab3 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 75f82c7026f..1df2865f50b 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 ebc921a8272..a28464a248c 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 15c6d91942e..82c6a6a6e9f 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 8d279153ff2..c27e1dd131d 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 a4fc7f14620..2327a6cdd05 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 f9ea4cd1fff..42c26c0e440 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 9f64039b0fb..37df1c7408a 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 8f1e7c6ea74..a8e70ea1d92 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 7150cd802b1..686f5507457 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 a443c5f6846..7970529cc44 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 116a73c4ad1..065b8e49b47 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 2f8d5dff4a1..f98e767f925 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 67d6370bb0d..9d88ede309a 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)
 
 
 
-- 
GitLab