Commit 667965c1 authored by Martyn Gigg's avatar Martyn Gigg
Browse files

Generate points within sample in a uniform manner

parent 8ece4d65
......@@ -2,13 +2,13 @@
#define MANTID_ALGORITHMS_MCINTERACTIONVOLUME_H_
#include "MantidAlgorithms/DllConfig.h"
#include "MantidGeometry/Objects/BoundingBox.h"
namespace Mantid {
namespace API {
class Sample;
}
namespace Geometry {
class BoundingBox;
class Object;
class SampleEnvironment;
}
......@@ -17,6 +17,7 @@ class PseudoRandomNumberGenerator;
class V3D;
}
namespace Algorithms {
class IBeamProfile;
/**
Defines a volume where interactions of Tracks and Objects can take place.
......@@ -54,13 +55,13 @@ public:
const Geometry::BoundingBox &getBoundingBox() const;
double calculateAbsorption(Kernel::PseudoRandomNumberGenerator &rng,
const Kernel::V3D &startPos,
const Kernel::V3D &direc,
const Kernel::V3D &endPos, double lambdaBefore,
double lambdaAfter) const;
private:
const Geometry::Object &m_sample;
const Geometry::SampleEnvironment *m_env;
const Geometry::BoundingBox m_activeRegion;
};
} // namespace Algorithms
......
......@@ -42,18 +42,16 @@ std::tuple<double, double>
MCAbsorptionStrategy::calculate(Kernel::PseudoRandomNumberGenerator &rng,
const Kernel::V3D &finalPos,
double lambdaBefore, double lambdaAfter) const {
// The simulation consists of sampling the beam profile and then computing the
// absorption within the interacting volume defined by the sample (and
// environment. The final correction factor is computed as the simple average
// of m_nevents of this process.
const auto scatterBounds = m_scatterVol.getBoundingBox();
double factor(0.0);
for (size_t i = 0; i < m_nevents; ++i) {
size_t attempts(0);
do {
const auto neutron = m_beamProfile.generatePoint(rng, scatterBounds);
const double wgt = m_scatterVol.calculateAbsorption(
rng, neutron.startPos, neutron.unitDir, finalPos, lambdaBefore,
rng, neutron.startPos, finalPos, lambdaBefore,
lambdaAfter);
if (wgt < 0.0) {
++attempts;
......
......@@ -13,6 +13,10 @@ using Kernel::V3D;
namespace Algorithms {
namespace {
// Maximum number of attempts to generate a scatter point
constexpr size_t MAX_SCATTER_ATTEMPTS = 500;
/**
* Compute the attenuation factor for the given coefficients
* @param rho Number density of the sample in \f$\\A^{-3}\f$
......@@ -27,12 +31,18 @@ double attenuation(double rho, double sigma, double length) {
}
/**
* Construct the volume with only a sample object
* Construct the volume encompassing the sample + any environment kit. The
* beam profile defines a bounding region for the sampling of the scattering
* position.
* @param sample A reference to a sample object that defines a valid shape
* & material
* @param beam A reference to the beam profile to restrict where the position
* of scattering is generated. Used when the beam size is smaller than the
* sample.
*/
MCInteractionVolume::MCInteractionVolume(const API::Sample &sample)
: m_sample(sample.getShape()), m_env(nullptr) {
: m_sample(sample.getShape()), m_env(nullptr),
m_activeRegion(this->getBoundingBox()) {
if (!m_sample.hasValidShape()) {
throw std::invalid_argument(
"MCInteractionVolume() - Sample shape does not have a valid shape.");
......@@ -57,11 +67,11 @@ const Geometry::BoundingBox &MCInteractionVolume::getBoundingBox() const {
}
/**
* Calculate the attenuation correction factor for given track through the
* volume
* @param rng A reference to a PseudoRandomNumberGenerator
* Calculate the attenuation correction factor the volume given a start and
* end point.
* @param rng A reference to a PseudoRandomNumberGenerator producing
* random number between [0,1]
* @param startPos Origin of the initial track
* @param direc Direction of travel of the neutron
* @param endPos Final position of neutron after scattering (assumed to be
* outside of the "volume")
* @param lambdaBefore Wavelength, in \f$\\A^-1\f$, before scattering
......@@ -71,68 +81,60 @@ const Geometry::BoundingBox &MCInteractionVolume::getBoundingBox() const {
*/
double MCInteractionVolume::calculateAbsorption(
Kernel::PseudoRandomNumberGenerator &rng, const Kernel::V3D &startPos,
const Kernel::V3D &direc, const Kernel::V3D &endPos, double lambdaBefore,
double lambdaAfter) const {
// Create track with start position and direction and "fire" it through
// the sample to produce a number of intersections. Choose a random
// intersection and within this section pick a random "depth". This point
// is the scatter point.
// Form a second track originating at the scatter point and ending at endPos
// to give a second set of intersections.
// The total attenuation factor is the product of the attenuation factor
// for each intersection
// Generate scatter point
Track path1(startPos, direc);
int nsegments = m_sample.interceptSurface(path1);
const Kernel::V3D &endPos, double lambdaBefore, double lambdaAfter) const {
// Generate scatter point. If there is an environment present then
// first select whether the scattering occurs on the sample or the
// environment. The attenuation for the path leading to the scatter point
// 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, MAX_SCATTER_ATTEMPTS);
} else {
scatterPos = m_sample.generatePointInObject(rng, m_activeRegion,
MAX_SCATTER_ATTEMPTS);
}
auto toStart = startPos - scatterPos;
toStart.normalize();
Track beforeScatter(scatterPos, toStart);
int nlinks = m_sample.interceptSurface(beforeScatter);
if (m_env) {
nsegments += m_env->interceptSurfaces(path1);
nlinks += m_env->interceptSurfaces(beforeScatter);
}
if (nsegments == 0) {
// This should not happen but numerical precision means that it can
// occasionally occur with tracks that are very close to the surface
if (nlinks == 0) {
return -1.0;
}
int scatterSegmentNo(1);
if (nsegments != 1) {
scatterSegmentNo = rng.nextInt(1, nsegments);
}
double atten(1.0);
V3D scatterPos;
auto segItr(path1.cbegin());
for (int i = 0; i < scatterSegmentNo; ++i, ++segItr) {
double length = segItr->distInsideObject;
if (i == scatterSegmentNo - 1) {
length *= rng.nextValue();
scatterPos = segItr->entryPoint + direc * length;
// Function to calculate total attenuation for a track
auto calculateAttenuation = [](const Track &path, double lambda) {
double factor(1.0);
for (const auto &segment : path) {
const double length = segment.distInsideObject;
const auto &segObj = *(segment.object);
const auto &segMat = segObj.material();
factor *= attenuation(segMat.numberDensity(),
segMat.totalScatterXSection(lambda) +
segMat.absorbXSection(lambda),
length);
}
const auto &segObj = *(segItr->object);
const auto &segMat = segObj.material();
atten *= attenuation(segMat.numberDensity(),
segMat.totalScatterXSection(lambdaBefore) +
segMat.absorbXSection(lambdaBefore),
length);
}
return factor;
};
// Now track to final destination
V3D scatteredDirec = endPos - scatterPos;
scatteredDirec.normalize();
Track path2(scatterPos, scatteredDirec);
m_sample.interceptSurface(path2);
Track afterScatter(scatterPos, scatteredDirec);
m_sample.interceptSurface(afterScatter);
if (m_env) {
m_env->interceptSurfaces(path2);
}
for (const auto &segment : path2) {
double length = segment.distInsideObject;
const auto &segObj = *(segment.object);
const auto &segMat = segObj.material();
atten *= attenuation(segMat.numberDensity(),
segMat.totalScatterXSection(lambdaAfter) +
segMat.absorbXSection(lambdaAfter),
length);
m_env->interceptSurfaces(afterScatter);
}
return atten;
return calculateAttenuation(beforeScatter, lambdaBefore) *
calculateAttenuation(afterScatter, lambdaAfter);
}
} // namespace Algorithms
......
......@@ -39,8 +39,7 @@ public:
using namespace ::testing;
// Testing inputs
const V3D startPos(-2.0, 0.0, 0.0), direc(1.0, 0.0, 0.0),
endPos(0.7, 0.7, 1.4);
const V3D startPos(-2.0, 0.0, 0.0), endPos(0.7, 0.7, 1.4);
const double lambdaBefore(2.5), lambdaAfter(3.5);
MockRNG rng;
EXPECT_CALL(rng, nextInt(1, 1)).Times(Exactly(0));
......@@ -49,7 +48,7 @@ public:
auto sample = createTestSample(TestSampleType::SolidSphere);
MCInteractionVolume interactor(sample);
const double factor = interactor.calculateAbsorption(
rng, startPos, direc, endPos, lambdaBefore, lambdaAfter);
rng, startPos, endPos, lambdaBefore, lambdaAfter);
TS_ASSERT_DELTA(1.06797501e-02, factor, 1e-8);
}
......@@ -59,8 +58,7 @@ public:
using namespace ::testing;
// Testing inputs
const V3D startPos(-2.0, 0.0, 0.0), direc(1.0, 0.0, 0.0),
endPos(2.0, 0.0, 0.0);
const V3D startPos(-2.0, 0.0, 0.0), endPos(2.0, 0.0, 0.0);
const double lambdaBefore(2.5), lambdaAfter(3.5);
auto sample = createTestSample(TestSampleType::Annulus);
......@@ -71,7 +69,7 @@ public:
MCInteractionVolume interactor(sample);
const double factorSeg1 = interactor.calculateAbsorption(
rng, startPos, direc, endPos, lambdaBefore, lambdaAfter);
rng, startPos, endPos, lambdaBefore, lambdaAfter);
TS_ASSERT_DELTA(5.35624555e-02, factorSeg1, 1e-8);
Mock::VerifyAndClearExpectations(&rng);
......@@ -79,7 +77,7 @@ public:
EXPECT_CALL(rng, nextInt(1, 2)).Times(Exactly(1)).WillOnce(Return(2));
EXPECT_CALL(rng, nextValue()).Times(Exactly(1)).WillOnce(Return(0.35));
const double factorSeg2 = interactor.calculateAbsorption(
rng, startPos, direc, endPos, lambdaBefore, lambdaAfter);
rng, startPos, endPos, lambdaBefore, lambdaAfter);
TS_ASSERT_DELTA(7.30835693e-02, factorSeg2, 1e-8);
Mock::VerifyAndClearExpectations(&rng);
}
......@@ -91,8 +89,7 @@ public:
using namespace ::testing;
// Testing inputs
const V3D startPos(-2.0, 0.0, 0.0), direc(1.0, 0.0, 0.0),
endPos(2.0, 0.0, 0.0);
const V3D startPos(-2.0, 0.0, 0.0), endPos(2.0, 0.0, 0.0);
const double lambdaBefore(2.5), lambdaAfter(3.5);
auto sample = createTestSample(TestSampleType::SamplePlusContainer);
......@@ -103,7 +100,7 @@ public:
MCInteractionVolume interactor(sample);
const double factorContainer = interactor.calculateAbsorption(
rng, startPos, direc, endPos, lambdaBefore, lambdaAfter);
rng, startPos, endPos, lambdaBefore, lambdaAfter);
TS_ASSERT_DELTA(6.919239804e-01, factorContainer, 1e-8);
Mock::VerifyAndClearExpectations(&rng);
......@@ -112,7 +109,7 @@ public:
EXPECT_CALL(rng, nextValue()).Times(Exactly(1)).WillOnce(Return(0.35));
const double factorSample = interactor.calculateAbsorption(
rng, startPos, direc, endPos, lambdaBefore, lambdaAfter);
rng, startPos, endPos, lambdaBefore, lambdaAfter);
TS_ASSERT_DELTA(6.9620991317e-01, factorSample, 1e-8);
Mock::VerifyAndClearExpectations(&rng);
}
......@@ -123,15 +120,14 @@ public:
using namespace ::testing;
// Testing inputs
const V3D startPos(-2.0, 0.0, 0.0), direc(0.0, 1.0, 0.0),
endPos(0.7, 0.7, 1.4);
const V3D startPos(-2.0, 0.0, 0.0), endPos(0.7, 0.7, 1.4);
const double lambdaBefore(2.5), lambdaAfter(3.5);
MockRNG rng;
EXPECT_CALL(rng, nextValue()).Times(Exactly(0));
auto sample = createTestSample(TestSampleType::SolidSphere);
MCInteractionVolume interactor(sample);
TS_ASSERT(interactor.calculateAbsorption(rng, startPos, direc, endPos,
TS_ASSERT(interactor.calculateAbsorption(rng, startPos, endPos,
lambdaBefore, lambdaAfter) < 0.0);
}
......
......@@ -8,6 +8,9 @@
#include "MantidGeometry/Instrument/Container.h"
namespace Mantid {
namespace Kernel {
class PseudoRandomNumberGenerator;
}
namespace Geometry {
class Track;
......@@ -52,6 +55,14 @@ public:
inline size_t nelements() const { return m_components.size(); }
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;
......
......@@ -142,9 +142,10 @@ public:
std::vector<Kernel::V3D> const &getCoordSystem() const {
return coord_system;
}
//@}
/// Generate a random point within the box
Kernel::V3D generatePointInside(double r1, double r2, double r3) const;
/** returns the expanded box consisting of all 8 box points,
* shifted into the coordinate system with the observer centre; */
void getFullBox(std::vector<Kernel::V3D> &box,
......
......@@ -15,17 +15,18 @@ namespace Mantid {
// Forward declarations
//----------------------------------------------------------------------
namespace Kernel {
class V3D;
class PseudoRandomNumberGenerator;
class Material;
class V3D;
}
namespace Geometry {
class Rule;
class CacheGeometryHandler;
class CompGrp;
class GeometryHandler;
class Rule;
class Surface;
class Track;
class GeometryHandler;
class CacheGeometryHandler;
class vtkGeometryCacheReader;
class vtkGeometryCacheWriter;
......@@ -153,6 +154,13 @@ public:
// find internal point to object
int getPointInObject(Kernel::V3D &point) const;
/// Select a random point within the object
Kernel::V3D generatePointInObject(Kernel::PseudoRandomNumberGenerator &rng,
const size_t) const;
Kernel::V3D generatePointInObject(Kernel::PseudoRandomNumberGenerator &rng,
const BoundingBox &activeRegion,
const size_t) const;
// Rendering member functions
void draw() const;
// Initialize Drawing
......
......@@ -143,8 +143,8 @@ struct IntersectionPoint {
*/
class MANTID_GEOMETRY_DLL Track {
public:
typedef std::list<Link> LType; ///< Type for the Link storage
typedef std::list<IntersectionPoint> PType; ///< Type for the partial
using LType = std::list<Link>;
using PType = std::list<IntersectionPoint>;
public:
/// Default constructor
......@@ -175,10 +175,24 @@ public:
LType::iterator begin() { return m_links.begin(); }
/// Returns an interator to one-past-the-end of the set of links
LType::iterator end() { return m_links.end(); }
/// Returns an interator to the start of the set of links
/// Returns an interator to the start of the set of links (const version)
LType::const_iterator begin() const { return m_links.begin(); }
/// Returns an interator to one-past-the-end of the set of links (const
/// version)
LType::const_iterator end() const { return m_links.end(); }
/// Returns an interator to the start of the set of links (const version)
LType::const_iterator cbegin() const { return m_links.cbegin(); }
/// Returns an interator to one-past-the-end of the set of links
/// Returns an interator to one-past-the-end of the set of links (const
/// version)
LType::const_iterator cend() const { return m_links.cend(); }
/// Returns a reference to the first link
LType::reference front() { return m_links.front(); }
/// Returns a reference to the last link
LType::reference back() { return m_links.back(); }
/// Returns a reference to the first link (const version)
LType::const_reference front() const { return m_links.front(); }
/// Returns a reference to the last link (const version)
LType::const_reference back() const { return m_links.back(); }
/// Returns the number of links
int count() const { return static_cast<int>(m_links.size()); }
/// Is the link complete?
......
......@@ -5,6 +5,7 @@
#include "MantidGeometry/IObjComponent.h"
#include "MantidGeometry/Objects/Object.h"
#include "MantidGeometry/Objects/Track.h"
#include "MantidKernel/PseudoRandomNumberGenerator.h"
namespace Mantid {
namespace Geometry {
......@@ -38,6 +39,41 @@ 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
......
......@@ -147,7 +147,21 @@ double BoundingBox::angularWidth(const Kernel::V3D &observer) const {
}
return thetaMax;
}
//
/**
* Generate a random point within this box assuming the 3 numbers given
* are random numbers in the range (0,1) & selected from a flat distribution.
* @param r1 Flat random number in range (0,1)
* @param r2 Flat random number in range (0,1)
* @param r3 Flat random number in range (0,1)
* @return A new point within the box such that isPointInside(pt) == true
*/
Kernel::V3D BoundingBox::generatePointInside(double r1, double r2,
double r3) const {
return V3D(xMin() + r1 * (xMax() - xMin()), yMin() + r2 * (yMax() - yMin()),
zMin() + r3 * (zMax() - zMin()));
}
void BoundingBox::getFullBox(std::vector<Kernel::V3D> &box,
const Kernel::V3D &observer) const {
box.resize(8);
......@@ -160,6 +174,7 @@ void BoundingBox::getFullBox(std::vector<Kernel::V3D> &box,
box[6] = Kernel::V3D(xMax(), yMin(), zMax()) - observer;
box[7] = Kernel::V3D(xMax(), yMax(), zMax()) - observer;
}
void BoundingBox::setBoxAlignment(const Kernel::V3D &R0,
const std::vector<Kernel::V3D> &orts) {
this->coord_system.resize(4);
......@@ -169,6 +184,7 @@ void BoundingBox::setBoxAlignment(const Kernel::V3D &R0,
coord_system[3] = orts[2];
is_axis_aligned = false;
}
void BoundingBox::nullify() {
this->m_null = true;
for (int i = 0; i < 3; i++) {
......@@ -176,7 +192,7 @@ void BoundingBox::nullify() {
this->m_maxPoint[i] = -FLT_MAX;
}
}
//
void BoundingBox::realign(std::vector<Kernel::V3D> const *const pCS) {
if (pCS) {
this->coord_system.resize(pCS->size());
......@@ -224,6 +240,7 @@ void BoundingBox::realign(std::vector<Kernel::V3D> const *const pCS) {
this->zMin() = zMin;
this->zMax() = zMax;
}
/**
* Enlarges this bounding box so that it encompasses that given.
* @param other :: The bounding box that should be encompassed
......
......@@ -4,6 +4,7 @@
#include "MantidKernel/Exception.h"
#include "MantidKernel/Material.h"
#include "MantidKernel/MultiThreaded.h"
#include "MantidKernel/PseudoRandomNumberGenerator.h"
#include "MantidKernel/Strings.h"
#include "MantidGeometry/Surfaces/Cone.h"
......@@ -1823,6 +1824,54 @@ int Object::getPointInObject(Kernel::V3D &point) const {
return 0;
}
/**
* Generate a random point within the object. The method simply generates a
* point within the bounding box and tests if this is a valid point within
* the object: if so the point is return otherwise a new point is selected.
* @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
*/
V3D Object::generatePointInObject(Kernel::PseudoRandomNumberGenerator &rng,
const size_t maxAttempts) 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);
}
/**
* Generate a random point within the object that is also bound by the
* activeRegion box.
* @param rng A reference to a PseudoRandomNumberGenerator where
* nextValue should return a flat random number between 0.0 & 1.0
* @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 Object::generatePointInObject(Kernel::PseudoRandomNumberGenerator &rng,
const BoundingBox &activeRegion,
const size_t maxAttempts) const {
size_t attempts(0);
while (attempts < maxAttempts) {
const double r1 = rng.nextValue();
const double r2 = rng.nextValue();
const double r3 = rng.nextValue();
auto pt = activeRegion.generatePointInside(r1, r2, r3);
if (this->isValid(pt))
return pt;
else
++attempts;
};
throw std::runtime_error("Object::generatePointInObject() - Unable to "
"generate point in object after " +
std::to_string(maxAttempts) + " attempts");
}
/**
* Try to find a point that lies within (or on) the object, given a seed point
* @param point :: on entry the seed point, on exit point in object, if found
......
......@@ -209,6 +209,18 @@ public:
TS_ASSERT_EQUALS(box.maxPoint() =