Commit 1a1bd93e authored by Danny Hindson's avatar Danny Hindson
Browse files

Move calculateAttenuation onto Track class

Move the calculateAttenuation method from MCInteractionVolume onto the
Track class so it's available for situations where you're not doing a
Monte Carlo simulation eg non-simulated absorption weighted path length
New X-ray absorption correction algorithm may also benefit from this

Also return tracks from MCInteractionVolume::calculateBeforeAfterTrack in
different way so it's easier to mock the Track class in the MCAbsorptionStrategy
unit tests
Also had to use testing::invoke to return the mocked track class because
SetArgPointee didn't seem to work in a polymorphic way
parent a9cce4b0
......@@ -34,14 +34,12 @@ struct ComponentScatterPoint {
class MANTID_ALGORITHMS_DLL IMCInteractionVolume {
public:
virtual ~IMCInteractionVolume() = default;
virtual bool calculateBeforeAfterTrack(
Kernel::PseudoRandomNumberGenerator &rng, const Kernel::V3D &startPos,
const Kernel::V3D &endPos, Geometry::Track &beforeScatter,
Geometry::Track &afterScatter, MCInteractionStatistics &stats) const = 0;
virtual double calculateAttenuation(const Geometry::Track &beforeScatter,
const Geometry::Track &afterScatter,
double lambdaBefore,
double lambdaAfter) const = 0;
virtual std::tuple<bool, std::shared_ptr<Geometry::Track>,
std::shared_ptr<Geometry::Track>>
calculateBeforeAfterTrack(Kernel::PseudoRandomNumberGenerator &rng,
const Kernel::V3D &startPos,
const Kernel::V3D &endPos,
MCInteractionStatistics &stats) const = 0;
virtual const Geometry::BoundingBox &getBoundingBox() const = 0;
virtual const Geometry::BoundingBox getFullBoundingBox() const = 0;
virtual void setActiveRegion(const Geometry::BoundingBox &region) = 0;
......
......@@ -39,15 +39,12 @@ public:
const Geometry::BoundingBox &getBoundingBox() const override;
const Geometry::BoundingBox getFullBoundingBox() const override;
virtual bool calculateBeforeAfterTrack(
Kernel::PseudoRandomNumberGenerator &rng, const Kernel::V3D &startPos,
const Kernel::V3D &endPos, Geometry::Track &beforeScatter,
Geometry::Track &afterScatter,
MCInteractionStatistics &stats) const override;
virtual double calculateAttenuation(const Geometry::Track &beforeScatter,
const Geometry::Track &afterScatter,
double lambdaBefore,
double lambdaAfter) const override;
virtual std::tuple<bool, std::shared_ptr<Geometry::Track>,
std::shared_ptr<Geometry::Track>>
calculateBeforeAfterTrack(Kernel::PseudoRandomNumberGenerator &rng,
const Kernel::V3D &startPos,
const Kernel::V3D &endPos,
MCInteractionStatistics &stats) const override;
ComponentScatterPoint
generatePoint(Kernel::PseudoRandomNumberGenerator &rng) const;
void setActiveRegion(const Geometry::BoundingBox &region) override;
......
......@@ -152,8 +152,8 @@ void AddAbsorptionWeightedPathLengths::exec() {
Track afterScatter(samplePos, detDir);
sampleShape->interceptSurface(afterScatter);
absFactors[0] = interactionVol.calculateAttenuation(
beforeScatter, afterScatter, lambdas[0], lambdas[0]);
absFactors[0] = beforeScatter.calculateAttenuation(lambdas[0]) *
afterScatter.calculateAttenuation(lambdas[0]);
} else {
MersenneTwister rng(seed + int(i));
......
......@@ -86,17 +86,17 @@ void MCAbsorptionStrategy::calculate(Kernel::PseudoRandomNumberGenerator &rng,
wgtM2(attenuationFactors.size());
for (size_t i = 0; i < m_nevents; ++i) {
Geometry::Track beforeScatter;
Geometry::Track afterScatter;
std::shared_ptr<Geometry::Track> beforeScatter;
std::shared_ptr<Geometry::Track> afterScatter;
for (int j = 0; j < nbins; ++j) {
size_t attempts(0);
do {
bool success = false;
if (m_regenerateTracksForEachLambda || j == 0) {
const auto neutron = m_beamProfile.generatePoint(rng, scatterBounds);
success = m_scatterVol.calculateBeforeAfterTrack(
rng, neutron.startPos, finalPos, beforeScatter, afterScatter,
stats);
std::tie(success, beforeScatter, afterScatter) =
m_scatterVol.calculateBeforeAfterTrack(rng, neutron.startPos,
finalPos, stats);
} else {
success = true;
}
......@@ -112,8 +112,8 @@ void MCAbsorptionStrategy::calculate(Kernel::PseudoRandomNumberGenerator &rng,
} else {
// elastic case already initialized
}
const double wgt = m_scatterVol.calculateAttenuation(
beforeScatter, afterScatter, lambdaIn, lambdaOut);
const double wgt = beforeScatter->calculateAttenuation(lambdaIn) *
afterScatter->calculateAttenuation(lambdaOut);
attenuationFactors[j] += wgt;
// increment standard deviation using Welford algorithm
double delta = wgt - wgtMean[j];
......
......@@ -160,17 +160,17 @@ ComponentScatterPoint MCInteractionVolume::generatePoint(
* @param startPos Origin of the initial track
* @param endPos Final position of neutron after scattering (assumed to be
* outside of the "volume")
* @param beforeScatter Out parameter to return generated before scatter track
* @param afterScatter Out parameter to return generated after scatter track
* @param stats A statistics class to hold the statistics on the generated
* tracks such as the scattering angle and count of scatter points generated in
* each sample or environment part
* @return Whether before/after tracks were successfully generated
* @return A tuple containing a flag to indicate whether before/after tracks
* were successfully generated and (if yes) the before/after tracks
*/
bool MCInteractionVolume::calculateBeforeAfterTrack(
std::tuple<bool, std::shared_ptr<Geometry::Track>,
std::shared_ptr<Geometry::Track>>
MCInteractionVolume::calculateBeforeAfterTrack(
Kernel::PseudoRandomNumberGenerator &rng, const Kernel::V3D &startPos,
const Kernel::V3D &endPos, Geometry::Track &beforeScatter,
Geometry::Track &afterScatter, MCInteractionStatistics &stats) const {
const Kernel::V3D &endPos, MCInteractionStatistics &stats) 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
......@@ -183,56 +183,29 @@ bool MCInteractionVolume::calculateBeforeAfterTrack(
stats.UpdateScatterPointCounts(scatterPos.componentIndex, false);
const auto toStart = normalize(startPos - scatterPos.scatterPoint);
beforeScatter = Track(scatterPos.scatterPoint, toStart);
int nlinks = m_sample->interceptSurface(beforeScatter);
auto beforeScatter =
std::make_shared<Track>(scatterPos.scatterPoint, toStart);
int nlinks = m_sample->interceptSurface(*beforeScatter);
if (m_env) {
nlinks += m_env->interceptSurfaces(beforeScatter);
nlinks += m_env->interceptSurfaces(*beforeScatter);
}
// 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 false;
return {false, nullptr, nullptr};
}
stats.UpdateScatterPointCounts(scatterPos.componentIndex, true);
// Now track to final destination
const V3D scatteredDirec = normalize(endPos - scatterPos.scatterPoint);
afterScatter = Track(scatterPos.scatterPoint, scatteredDirec);
m_sample->interceptSurface(afterScatter);
auto afterScatter =
std::make_shared<Track>(scatterPos.scatterPoint, scatteredDirec);
m_sample->interceptSurface(*afterScatter);
if (m_env) {
m_env->interceptSurfaces(afterScatter);
m_env->interceptSurfaces(*afterScatter);
}
stats.UpdateScatterAngleStats(toStart, scatteredDirec);
return true;
}
/**
* Calculate the attenuation correction factor the volume given a before
* and after track.
* @param beforeScatter Before scatter track
* @param afterScatter After scatter track
* @param lambdaBefore Lambda before scattering
* @param lambdaAfter Lambda after scattering
* @return Absorption factor
*/
double MCInteractionVolume::calculateAttenuation(const Track &beforeScatter,
const Track &afterScatter,
double lambdaBefore,
double lambdaAfter) const {
// Function to calculate total attenuation for a track
auto calculateAttenuationForTrack = [](const Track &path, double lambda) {
double factor(1.0);
for (const auto &segment : path) {
const double length = segment.distInsideObject;
const auto &segObj = *(segment.object);
factor *= segObj.material().attenuation(length, lambda);
}
return factor;
};
return calculateAttenuationForTrack(beforeScatter, lambdaBefore) *
calculateAttenuationForTrack(afterScatter, lambdaAfter);
return {true, beforeScatter, afterScatter};
}
} // namespace Algorithms
......
......@@ -17,6 +17,7 @@
#include "MantidDataObjects/Histogram1D.h"
#include "MantidGeometry/Instrument/ReferenceFrame.h"
#include "MantidGeometry/Objects/BoundingBox.h"
#include "MantidGeometry/Objects/Track.h"
#include "MantidKernel/Logger.h"
#include "MantidKernel/WarningSuppressions.h"
#include "MonteCarloTesting.h"
......@@ -168,22 +169,45 @@ public:
.Times(1)
.WillOnce(ReturnRef(emptyBoundingBox));
EXPECT_CALL(testInteractionVolume,
calculateBeforeAfterTrack(_, _, _, _, _, _))
auto beforeScatter = std::make_shared<MockTrack>();
auto afterScatter = std::make_shared<MockTrack>();
auto ReturnTracksAndTrue = [beforeScatter,
afterScatter](auto &arg0, auto &arg1,
auto &arg2, auto &arg3) {
UNUSED_ARG(arg0);
UNUSED_ARG(arg1);
UNUSED_ARG(arg2);
UNUSED_ARG(arg3);
return std::tuple{true, beforeScatter, afterScatter};
};
auto ReturnNullTracksAndFalse = [beforeScatter,
afterScatter](auto &arg0, auto &arg1,
auto &arg2, auto &arg3) {
UNUSED_ARG(arg0);
UNUSED_ARG(arg1);
UNUSED_ARG(arg2);
UNUSED_ARG(arg3);
return std::tuple{false, nullptr, nullptr};
};
EXPECT_CALL(testInteractionVolume, calculateBeforeAfterTrack(_, _, _, _))
.Times(Exactly(6))
.WillOnce(Return(true))
.WillOnce(Return(false))
.WillOnce(Return(true))
.WillOnce(Return(true))
.WillOnce(Return(true))
.WillOnce(Return(true));
EXPECT_CALL(testInteractionVolume, calculateAttenuation(_, _, _, _))
.WillOnce(Invoke(ReturnTracksAndTrue))
.WillOnce(Invoke(ReturnNullTracksAndFalse))
.WillOnce(Invoke(ReturnTracksAndTrue))
.WillOnce(Invoke(ReturnTracksAndTrue))
.WillOnce(Invoke(ReturnTracksAndTrue))
.WillOnce(Invoke(ReturnTracksAndTrue));
EXPECT_CALL(*beforeScatter, calculateAttenuation(_))
.Times(Exactly(5))
.WillOnce(Return(1.0))
.WillOnce(Return(2.0))
.WillOnce(Return(3.0))
.WillOnce(Return(4.0))
.WillOnce(Return(5.0));
EXPECT_CALL(*afterScatter, calculateAttenuation(_))
.Times(Exactly(5))
.WillRepeatedly(Return(1.0));
testStrategy.calculate(rng, {0., 0., 0.}, {1.0}, 0., attenuationFactors,
attenuationFactorErrors, trackStatistics);
TS_ASSERT_EQUALS(attenuationFactors[0], 3.0);
......@@ -242,23 +266,25 @@ private:
class MockMCInteractionVolume final : public IMCInteractionVolume {
public:
GNU_DIAG_OFF_SUGGEST_OVERRIDE
MOCK_CONST_METHOD6(calculateBeforeAfterTrack,
bool(Mantid::Kernel::PseudoRandomNumberGenerator &rng,
const Mantid::Kernel::V3D &startPos,
const Mantid::Kernel::V3D &endPos,
Mantid::Geometry::Track &beforeScatter,
Mantid::Geometry::Track &afterScatter,
MCInteractionStatistics &stats));
MOCK_CONST_METHOD4(calculateAttenuation,
double(const Mantid::Geometry::Track &beforeScatter,
const Mantid::Geometry::Track &afterScatter,
double lambdaBefore, double lambdaAfter));
MOCK_CONST_METHOD4(
calculateBeforeAfterTrack,
std::tuple<bool, std::shared_ptr<Mantid::Geometry::Track>,
std::shared_ptr<Mantid::Geometry::Track>>(
Mantid::Kernel::PseudoRandomNumberGenerator &rng,
const Mantid::Kernel::V3D &startPos,
const Mantid::Kernel::V3D &endPos, MCInteractionStatistics &stats));
MOCK_CONST_METHOD0(getBoundingBox, Mantid::Geometry::BoundingBox &());
MOCK_CONST_METHOD0(getFullBoundingBox,
const Mantid::Geometry::BoundingBox());
MOCK_METHOD1(setActiveRegion, void(const Mantid::Geometry::BoundingBox &));
GNU_DIAG_ON_SUGGEST_OVERRIDE
};
class MockTrack final : public Mantid::Geometry::Track {
public:
GNU_DIAG_OFF_SUGGEST_OVERRIDE
MOCK_CONST_METHOD1(calculateAttenuation, double(double));
GNU_DIAG_ON_SUGGEST_OVERRIDE
};
Mantid::Kernel::Logger g_log{"MCAbsorptionStrategyTest"};
};
......@@ -60,36 +60,15 @@ public:
auto sample = createTestSample(TestSampleType::SolidSphere);
MCInteractionVolume interactor(sample);
Track beforeScatter;
Track afterScatter;
MCInteractionStatistics trackStatistics(-1, sample);
interactor.calculateBeforeAfterTrack(rng, startPos, endPos, beforeScatter,
afterScatter, trackStatistics);
auto [success, beforeScatter, afterScatter] =
interactor.calculateBeforeAfterTrack(rng, startPos, endPos,
trackStatistics);
TS_ASSERT(success);
TestGeneratedTracks(startPos, endPos, beforeScatter, afterScatter,
sample.getShape());
}
void test_Solid_Sample_Gives_Expected_Absorption() {
auto sample = createTestSample(TestSampleType::SolidSphere);
MCInteractionVolume interactor(sample);
const double lambdaBefore(2.5), lambdaAfter(3.5);
// use tracks generated by previous test
Track beforeScatter({-0.05, -0.05, -0.05},
{-0.999343185, 0.025624184, 0.025624184});
beforeScatter.addLink({-0.05, -0.05, -0.05},
{-0.071481137, -0.049449202, -0.049449202},
0.021495255, sample.getShape());
Track afterScatter({-0.05, -0.05, -0.05},
{0.417472754, 0.417472755, 0.807113993});
afterScatter.addLink({-0.05, -0.05, -0.05},
{0.024407241, 0.024407241, 0.093853999}, 0.021495255,
sample.getShape());
const double factor = interactor.calculateAttenuation(
beforeScatter, afterScatter, lambdaBefore, lambdaAfter);
TS_ASSERT_DELTA(0.0028357258, factor, 1e-8);
}
void test_Sample_With_Hole_Gives_Expected_Tracks() {
using Mantid::Kernel::V3D;
......@@ -104,11 +83,10 @@ public:
.WillRepeatedly(Return(0.25));
MCInteractionVolume interactor(sample);
Track beforeScatter;
Track afterScatter;
MCInteractionStatistics trackStatistics(-1, sample);
interactor.calculateBeforeAfterTrack(rng, startPos, endPos, beforeScatter,
afterScatter, trackStatistics);
auto [success, beforeScatter, afterScatter] =
interactor.calculateBeforeAfterTrack(rng, startPos, endPos,
trackStatistics);
TestGeneratedTracks(startPos, endPos, beforeScatter, afterScatter,
sample.getShape());
TS_ASSERT(Mock::VerifyAndClearExpectations(&rng));
......@@ -117,41 +95,18 @@ public:
EXPECT_CALL(rng, nextValue())
.Times(Exactly(3))
.WillRepeatedly(Return(0.75));
interactor.calculateBeforeAfterTrack(rng, startPos, endPos, beforeScatter,
afterScatter, trackStatistics);
std::tie(success, beforeScatter, afterScatter) =
interactor.calculateBeforeAfterTrack(rng, startPos, endPos,
trackStatistics);
TestGeneratedTracks(startPos, endPos, beforeScatter, afterScatter,
sample.getShape());
}
void test_Sample_With_Hole_Gives_Expected_Absorption() {
auto sample = createTestSample(TestSampleType::Annulus);
MCInteractionVolume interactor(sample);
const double lambdaBefore(2.5), lambdaAfter(3.5);
// use tracks generated by previous test
Track beforeScatter({-0.075, -0.075, -0.0375},
{-0.999052621, 0.038924128, 0.019462064});
beforeScatter.addLink({-0.075, -0.075, -0.0375},
{-0.131142366, -0.072812635, -0.036406318},
0.056195605, sample.getShape());
Track afterScatter({-0.075, -0.075, -0.0375},
{0.999184480, 0.036115102, 0.018057551});
afterScatter.addLink({-0.075, -0.075, -0.0375},
{-0.066490895, -0.074692442, -0.037346221},
0.008516050, sample.getShape());
afterScatter.addLink({0.071709799, -0.069697236, -0.034848618},
{0.133981248, -0.067446460, -0.033723230}, 0.209151816,
sample.getShape());
const double factorSeg1 = interactor.calculateAttenuation(
beforeScatter, afterScatter, lambdaBefore, lambdaAfter);
TS_ASSERT_DELTA(0.030489479, factorSeg1, 1e-8);
}
void test_Sample_And_Environment_Can_Scatter_In_All_Segments() {
using Mantid::Kernel::V3D;
// Testing inputs
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);
MockRNG rng;
......@@ -167,11 +122,10 @@ public:
MCInteractionVolume interactor(sample);
interactor.setActiveRegion(sample.getEnvironment().boundingBox());
Track beforeScatter;
Track afterScatter;
MCInteractionStatistics trackStatistics(-1, sample);
interactor.calculateBeforeAfterTrack(rng, startPos, endPos, beforeScatter,
afterScatter, trackStatistics);
auto [success, beforeScatter, afterScatter] =
interactor.calculateBeforeAfterTrack(rng, startPos, endPos,
trackStatistics);
TestGeneratedTracks(startPos, endPos, beforeScatter, afterScatter,
sample.getEnvironment().getContainer());
TS_ASSERT(Mock::VerifyAndClearExpectations(&rng));
......@@ -183,60 +137,14 @@ public:
EXPECT_CALL(rng, nextValue())
.Times(Exactly(3))
.WillRepeatedly(Return(0.25)); // RandomPoint::bounded r1, r2, r3
interactor.calculateBeforeAfterTrack(rng, startPos, endPos, beforeScatter,
afterScatter, trackStatistics);
std::tie(success, beforeScatter, afterScatter) =
interactor.calculateBeforeAfterTrack(rng, startPos, endPos,
trackStatistics);
TestGeneratedTracks(startPos, endPos, beforeScatter, afterScatter,
sample.getShape());
const double factorSample = interactor.calculateAttenuation(
beforeScatter, afterScatter, lambdaBefore, lambdaAfter);
TS_ASSERT_DELTA(0.73100698, factorSample, 1e-8);
TS_ASSERT(Mock::VerifyAndClearExpectations(&rng));
}
void test_Sample_And_Environment_Gives_Expected_Absorption() {
auto sample = createTestSample(TestSampleType::SamplePlusContainer);
MCInteractionVolume interactor(sample);
const double lambdaBefore(2.5), lambdaAfter(3.5);
// use tracks generated by previous test
// scatter in segment can
Track beforeScatter({-0.0048, 0, 0}, {-1, 0, 0});
beforeScatter.addLink({-0.0048, 0, 0}, {-0.005, 0, 0}, 0.0002,
sample.getEnvironment().getContainer().getShape());
Track afterScatter({-0.0048, 0, 0}, {1, 0, 0});
afterScatter.addLink({-0.0048, 0, 0}, {-0.0046, 0, 0}, 0.0002,
sample.getEnvironment().getContainer().getShape());
afterScatter.addLink({-0.0046, 0, 0}, {0.0046, 0, 0}, 0.0094,
sample.getShape());
afterScatter.addLink({0.0046, 0, 0}, {0.005, 0, 0}, 0.0098,
sample.getEnvironment().getContainer().getShape());
const double factorContainer = interactor.calculateAttenuation(
beforeScatter, afterScatter, lambdaBefore, lambdaAfter);
TS_ASSERT_DELTA(0.69223681, factorContainer, 1e-8);
// scatter in sample
Track beforeScatterSample({-0.0025, -0.0025, -0.0125},
{-0.999979637, 0.001251539, 0.006257695});
beforeScatterSample.addLink({-0.0025, -0.0025, -0.0125},
{-0.003862450, -0.002498295, -0.012491474},
0.001362478, sample.getShape());
beforeScatterSample.addLink(
{-0.003862450, -0.002498295, -0.012491474},
{-0.004331450, -0.002497708, -0.012488539}, 0.001831487,
sample.getEnvironment().getContainer().getShape());
Track afterScatterSample({-0.0025, -0.0025, -0.0125},
{0.999979739, 0.001248414, 0.006242071});
afterScatterSample.addLink({-0.0025, -0.0025, -0.0125},
{0.003866481, -0.002492052, -0.012460259},
0.006366610, sample.getShape());
afterScatterSample.addLink(
{0.003866481, -0.002492052, -0.012460259},
{0.004335042, -0.002491467, -0.012457334}, 0.006835181,
sample.getEnvironment().getContainer().getShape());
const double factorSample = interactor.calculateAttenuation(
beforeScatterSample, afterScatterSample, lambdaBefore, lambdaAfter);
TS_ASSERT_DELTA(0.73100698, factorSample, 1e-8);
}
void test_Construction_With_Env_But_No_Sample_Shape_Does_Not_Throw_Error() {
Mantid::API::Sample sample;
auto kit = createTestKit();
......@@ -269,12 +177,9 @@ public:
rng.setSeed(1);
const size_t maxTries(1);
MCInteractionVolume interactor(sample, maxTries);
Track beforeScatter;
Track afterScatter;
MCInteractionStatistics trackStatistics(-1, sample);
TS_ASSERT_THROWS(interactor.calculateBeforeAfterTrack(
rng, startPos, endPos, beforeScatter, afterScatter,
trackStatistics),
TS_ASSERT_THROWS(interactor.calculateBeforeAfterTrack(rng, startPos, endPos,
trackStatistics),
const std::runtime_error &);
}
......@@ -368,19 +273,19 @@ private:
Mantid::Kernel::Logger g_log{"MCInteractionVolumeTest"};
void TestGeneratedTracks(const V3D startPos, const V3D endPos,
const Track &beforeScatter,
const Track &afterScatter,
const std::shared_ptr<Track> beforeScatter,
const std::shared_ptr<Track> &afterScatter,
const Mantid::Geometry::IObject &shape) {
// check the generated tracks share the same start point (the scatter point)
TS_ASSERT_EQUALS(beforeScatter.startPoint(), afterScatter.startPoint());
TS_ASSERT_EQUALS(shape.isValid(beforeScatter.startPoint()), true);
TS_ASSERT_EQUALS(beforeScatter->startPoint(), afterScatter->startPoint());
TS_ASSERT_EQUALS(shape.isValid(beforeScatter->startPoint()), true);
// check that the generated tracks intersect the supplied startPos and
// endPos
const V3D scatterToEnd = endPos - afterScatter.startPoint();
TS_ASSERT_EQUALS(afterScatter.direction(),
const V3D scatterToEnd = endPos - afterScatter->startPoint();
TS_ASSERT_EQUALS(afterScatter->direction(),
Mantid::Kernel::normalize(scatterToEnd));
const V3D scatterToStart = startPos - beforeScatter.startPoint();
TS_ASSERT_EQUALS(beforeScatter.direction(),
const V3D scatterToStart = startPos - beforeScatter->startPoint();
TS_ASSERT_EQUALS(beforeScatter->direction(),
Mantid::Kernel::normalize(scatterToStart));
}
};
......@@ -186,6 +186,7 @@ public:
Track();
/// Constructor
Track(const Kernel::V3D &startPt, const Kernel::V3D &unitVector);
virtual ~Track() = default;
/// Adds a point of intersection to the track
void addPoint(const TrackDirection direction, const Kernel::V3D &endPoint,
const IObject &obj, const ComponentID compID = nullptr);
......@@ -235,6 +236,8 @@ public:
int surfPointsCount() const { return static_cast<int>(m_surfPoints.size()); }
/// Is the link complete?
int nonComplete() const;
/// Calculate attenuation across all links
virtual double calculateAttenuation(double lambda) const;
private:
Line m_line; ///< Line object containing origin and direction
......
......@@ -6,6 +6,7 @@
// SPDX - License - Identifier: GPL - 3.0 +
#include "MantidGeometry/Objects/Track.h"
#include "MantidGeometry/Surfaces/Surface.h"
#include "MantidKernel/Material.h"
#include "MantidKernel/Matrix.h"
#include "MantidKernel/Tolerance.h"
#include "MantidKernel/V3D.h"
......@@ -286,6 +287,16 @@ std::ostream &operator<<(std::ostream &os, TrackDirection direction) {
return os;
}
double Track::calculateAttenuation(double lambda) const {
double factor(1.0);
for (const auto &segment : m_links) {
const double length = segment.distInsideObject;
const auto &segObj = *(segment.object);
factor *= segObj.material().attenuation(length, lambda);
}
return factor;
}
} // NAMESPACE Geometry
} // NAMESPACE Mantid
......@@ -10,8 +10,10 @@
#include "MantidGeometry/Objects/CSGObject.h"
#include "MantidGeometry/Objects/Track.h"
#include "MantidKernel/Logger.h"
#include "MantidKernel/Material.h"
#include "MantidKernel/System.h"
#include "MantidKernel/V3D.h"
#include "MantidTestHelpers/ComponentCreationHelper.h"
#include <cxxtest/TestSuite.h>
using namespace Mantid;
......@@ -155,4 +157,26 @@ public:
shape); // Entry at -5,-2,0
TS_ASSERT_EQUALS(A.surfPointsCount(), 2);
}
void test_calculateAttenuation() {
const double lambdaBefore(2.5), lambdaAfter(3.5);
auto shape = ComponentCreationHelper::createSphere(0.1);
shape->setMaterial(Kernel::Material(
"Vanadium", Mantid::PhysicalConstants::getNeutronAtom(23), 0.02));
// use tracks generated by
// MCInteractionVolume::test_Solid_Sample_Gives_Expected_Tracks test
Track beforeScatter({-0.05, -0.05, -0.05},
{-0.999343185, 0.025624184, 0.025624184});
beforeScatter.addLink({-0.05, -0.05, -0.05},
{-0.071481137, -0.049449202, -0.049449202},
0.021495255, *shape);