-
Gigg, Martyn Anthony authoredGigg, Martyn Anthony authored
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
MonteCarloAbsorption.cpp 14.64 KiB
//------------------------------------------------------------------------------
// Includes
//------------------------------------------------------------------------------
#include "MantidAlgorithms/MonteCarloAbsorption.h"
#include "MantidAPI/WorkspaceProperty.h"
#include "MantidAPI/WorkspaceValidators.h"
#include "MantidAPI/SampleEnvironment.h"
#include "MantidKernel/BoundedValidator.h"
#include "MantidGeometry/IDetector.h"
#include "MantidGeometry/Objects/Track.h"
#include "MantidKernel/VectorHelper.h"
#include "MantidKernel/NeutronAtom.h"
// @todo: This needs a factory
#include "MantidKernel/MersenneTwister.h"
/// @cond
namespace
{
/// Number of attempts to choose a random point within the object before it gives up
const int MaxRandPointAttempts(20);
}
/// @endcond
namespace Mantid
{
namespace Algorithms
{
DECLARE_ALGORITHM(MonteCarloAbsorption)
using API::WorkspaceProperty;
using API::WorkspaceUnitValidator;
using API::MatrixWorkspace_sptr;
using API::WorkspaceFactory;
using API::Progress;
using namespace Geometry;
using Kernel::Direction;
//------------------------------------------------------------------------------
// Public methods
//------------------------------------------------------------------------------
/**
* Constructor
*/
MonteCarloAbsorption::MonteCarloAbsorption() :
m_inputWS(), m_sampleShape(NULL), m_container(NULL), m_numberOfPoints(0),
m_xStepSize(0), m_numberOfEvents(1), m_samplePos(), m_sourcePos(),
m_bbox_length(0.0), m_bbox_halflength(0.0), m_bbox_width(0.0),
m_bbox_halfwidth(0.0), m_bbox_height(0.0), m_bbox_halfheight(0.0),
m_randGen(NULL)
{
}
/**
* Destructor
*/
MonteCarloAbsorption::~MonteCarloAbsorption()
{
delete m_randGen;
}
//------------------------------------------------------------------------------
// Private methods
//------------------------------------------------------------------------------
/**
* Initialize the algorithm
*/ void MonteCarloAbsorption::init()
{
declareProperty(new WorkspaceProperty<>("InputWorkspace", "", Direction::Input,
new WorkspaceUnitValidator<> ("Wavelength")),
"The X values for the input workspace must be in units of wavelength");
declareProperty(new WorkspaceProperty<> ("OutputWorkspace", "", Direction::Output),
"Output workspace name");
Kernel::BoundedValidator<int> *positiveInt = new Kernel::BoundedValidator<int> ();
positiveInt->setLower(1);
declareProperty("NumberOfWavelengthPoints", EMPTY_INT(), positiveInt,
"The number of wavelength points for which a simulation is\n"
"performed (default: all points)");
declareProperty("EventsPerPoint", 300, positiveInt->clone(),
"The number of events to simulate per wavelength point used.");
declareProperty("SeedValue", 123456789, positiveInt->clone(),
"A seed for the random number generator");
}
/**
* Execution code
*/
void MonteCarloAbsorption::exec()
{
retrieveInput();
initCaches();
MatrixWorkspace_sptr correctionFactors = WorkspaceFactory::Instance().create(m_inputWS);
correctionFactors->isDistribution(true); // The output of this is a distribution
correctionFactors->setYUnit(""); // Need to explicitly set YUnit to nothing
correctionFactors->setYUnitLabel("Attenuation factor");
const bool isHistogram = m_inputWS->isHistogramData();
const int numHists(m_inputWS->getNumberHistograms());
const int numBins = m_inputWS->blocksize();
// Compute the step size
m_xStepSize = numBins/m_numberOfPoints;
std::ostringstream message;
message << "Simulation performed every " << m_xStepSize << " wavelength points" << std::endl;
g_log.information(message.str());
message.str("");
Progress prog(this,0.0,1.0,numHists);
PARALLEL_FOR2(m_inputWS, correctionFactors)
for( int i = 0; i < numHists; ++i )
{
PARALLEL_START_INTERUPT_REGION
// Copy over the X-values
const MantidVec & xValues = m_inputWS->readX(i);
correctionFactors->dataX(i) = xValues;
// Final detector position
IDetector_const_sptr detector;
try
{
detector = m_inputWS->getDetector(i);
}
catch(Kernel::Exception::NotFoundError&)
{
continue;
}
MantidVec & yValues = correctionFactors->dataY(i);
MantidVec & eValues = correctionFactors->dataE(i);
// Simulation for each requested wavelength point
for( int bin = 0; bin < numBins; bin += m_xStepSize )
{
const double lambda = isHistogram ? (0.5 * (xValues[bin] + xValues[bin + 1]) ) : xValues[bin];
doSimulation(detector.get(), lambda, yValues[bin], eValues[bin]);
// Ensure we have the last point for the interpolation
if ( m_xStepSize > 1 && bin + m_xStepSize >= numBins && bin+1 != numBins)
{
bin = numBins - m_xStepSize - 1;
}
}
// Interpolate through points not simulated
if( m_xStepSize > 1 )
{
Kernel::VectorHelper::linearlyInterpolateY(xValues, yValues, m_xStepSize);
}
prog.report();
PARALLEL_END_INTERUPT_REGION
}
PARALLEL_CHECK_INTERUPT_REGION
// Save the results
setProperty("OutputWorkspace", correctionFactors);
}
/**
* Perform the simulation
* @param detector A pointer to the current detector
* @param lambda The chosen wavelength
* @param attenFactor [Output] The calculated attenuation factor for this wavelength
* @param error [Output] The value of the error on the factor
*/
void MonteCarloAbsorption::doSimulation(const IDetector *const detector, const double lambda,
double & attenFactor, double & error)
{
/**
Currently, assuming square beam profile to pick start position then randomly selecting
a point within the sample using it's bounding box.
This point defines the single scattering point and hence the attenuation path lengths and final
directional vector to the detector
*/
// Absolute detector position
const V3D detectorPos(detector->getPos());
int numDetected(0);
attenFactor = 0.0;
error = 0.0;
while( numDetected < m_numberOfEvents )
{
V3D startPos = sampleBeamProfile();
V3D scatterPoint = selectScatterPoint();
attenFactor += attenuationFactor(startPos, scatterPoint, detectorPos, lambda);
++numDetected;
}
// Attenuation factor is simply the average value
attenFactor /= numDetected;
// Error is 1/sqrt(nevents)
error = 1./sqrt((double)numDetected);
}
/**
* Sample the beam profile for a random start location, assigning a weight
* to the given point
* @returns The starting location
*/
V3D MonteCarloAbsorption::sampleBeamProfile() const
{
return m_sourcePos; }
/**
* Selects a random location within the sample + container environment.
* @returns Selected position as V3D object
*/
V3D MonteCarloAbsorption::selectScatterPoint() const
{
// Generate 3 random numbers, use them to calculate a random location
// within the bounding box of the sample environment + sample
// and then check if this position is within the whole environment.
// If it is return it, if not try again.
V3D scatterPoint;
int nattempts(0);
while( nattempts < MaxRandPointAttempts )
{
scatterPoint = V3D(m_bbox_halflength*(2.0*m_randGen->next() - 1.0),
m_bbox_halfwidth*(2.0*m_randGen->next() - 1.0),
m_bbox_halfheight*(2.0*m_randGen->next() - 1.0) );
++nattempts;
if( m_sampleShape->isValid(scatterPoint) ||
(m_container && m_container->isValid(scatterPoint)) )
{
scatterPoint += m_samplePos;
return scatterPoint;
}
}
// If we got here then the shape is too strange for the bounding box to be of any use.
g_log.error() << "Attempts to generat a random point with the sample/can "
<< "have exceeded the allowed number of tries.\n";
throw std::runtime_error("Attempts to produce random scatter point failed. Check sample shape.");
}
/**
* Return the attenuation factor for the given track
* @param startPos The origin of the track
* @param scatterPoint The point of scatter
* @param finalPos The end point of the track
* @param lambda The wavelength of the neutron
* @returns The attenuation factor for this neutron's track
*/
double
MonteCarloAbsorption::attenuationFactor(const V3D & startPos, const V3D & scatterPoint,
const V3D & finalPos, const double lambda)
{
double factor(1.0);
// Define two tracks, before and after scatter, and trace check their
// intersections with the the environment and sample
Track beforeScatter(scatterPoint, (startPos - scatterPoint));
Track afterScatter(scatterPoint, (finalPos - scatterPoint));
// Theoretically this should never happen as there should always be an intersection
// but do to precision limitations points very close to the surface give
// zero intersection, so just reject
if( m_sampleShape->interceptSurface(beforeScatter) == 0 ||
m_sampleShape->interceptSurface(afterScatter) == 0 )
{
return 0.0;
}
double length = beforeScatter.begin()->distInsideObject;
factor *= attenuation(length, *m_sampleMaterial, lambda);
beforeScatter.clearIntersectionResults();
if( m_container )
{
m_container->interceptSurfaces(beforeScatter);
}
// Attenuation factor is product of factor for each material Track::LType::const_iterator cend = beforeScatter.end();
for(Track::LType::const_iterator citr = beforeScatter.begin();
citr != cend; ++citr)
{
length = citr->distInsideObject;
IObjComponent *objComp = dynamic_cast<IObjComponent*>(citr->componentID);
Material_const_sptr mat = objComp->material();
factor *= attenuation(length, *mat, lambda);
}
length = afterScatter.begin()->distInsideObject;
factor *= attenuation(length, *m_sampleMaterial, lambda);
afterScatter.clearIntersectionResults();
if( m_container )
{
m_container->interceptSurfaces(afterScatter);
}
// Attenuation factor is product of factor for each material
cend = afterScatter.end();
for(Track::LType::const_iterator citr = afterScatter.begin();
citr != cend; ++citr)
{
length = citr->distInsideObject;
IObjComponent *objComp = dynamic_cast<IObjComponent*>(citr->componentID);
Material_const_sptr mat = objComp->material();
factor *= attenuation(length, *mat, lambda);
}
return factor;
}
/**
* Calculate the attenuation for a given length, material and wavelength
* @param length Distance through the material
* @param material A reference to the Material
* @param lambda The wavelength
* @returns The attenuation factor
*/
double
MonteCarloAbsorption::attenuation(const double length, const Geometry::Material& material,
const double lambda) const
{
const double rho = material.numberDensity() * 100.0;
const double sigma_s = material.totalScatterXSection(lambda);
const double sigma_t = sigma_s + material.absorbXSection(lambda);
return exp(-rho*sigma_t*length);
}
/**
* Gather the input values and check validity
* @throws std::invalid_argument If the input is invalid. Currently if there is
* no defined sample shape
*/
void MonteCarloAbsorption::retrieveInput()
{
m_inputWS = getProperty("InputWorkspace");
// // Define a test sample material
// Material *vanadium = new Material("Vanadium", PhysicalConstants::getNeutronAtom(23,0), 0.072);
// m_inputWS->mutableSample().setMaterial(*vanadium);
// // Define test environment
// std::ostringstream xml;
// xml << "<sphere id=\"" << "sp" << "\">"
// << "<centre x=\"" << 0.0 << "\" y=\"" << 0.0 << "\" z=\"" << 0.0 << "\" />"
// << "<radius val=\"" << 0.05 << "\" />"
// << "</sphere>";
// Geometry::ShapeFactory shapeMaker; // Object_sptr p = shapeMaker.createShape(xml.str());
// API::SampleEnvironment * kit = new API::SampleEnvironment("TestEnv");
// kit->add(new ObjComponent("one", p, NULL, boost::shared_ptr<Material>(vanadium)));
// m_inputWS->mutableSample().setEnvironment(kit);
m_sampleShape = &(m_inputWS->sample().getShape());
m_sampleMaterial = &(m_inputWS->sample().getMaterial());
if( !m_sampleShape->hasValidShape() )
{
g_log.debug() << "Invalid shape defined on workspace. TopRule = " << m_sampleShape->topRule()
<< ", No. of surfaces: " << m_sampleShape->getSurfacePtr().size() << "\n";
throw std::invalid_argument("Input workspace has an invalid sample shape.");
}
if( m_sampleMaterial->totalScatterXSection(1.0) == 0.0 )
{
g_log.warning() << "The sample material appears to have zero scattering cross section.\n"
<< "Result will most likely be nonsensical.\n";
}
try
{
m_container = &(m_inputWS->sample().getEnvironment());
}
catch(std::runtime_error&)
{
m_container = NULL;
g_log.information() << "No environment has been defined, continuing with only sample.\n";
}
m_numberOfPoints = getProperty("NumberOfWavelengthPoints");
if( isEmpty(m_numberOfPoints) || m_numberOfPoints > m_inputWS->blocksize() )
{
m_numberOfPoints = m_inputWS->blocksize();
if( !isEmpty(m_numberOfPoints) )
{
g_log.warning() << "The requested number of wavelength points is larger than the spectra size. "
<< "Defaulting to spectra size.\n";
}
}
m_numberOfEvents = getProperty("EventsPerPoint");
}
/**
* Initialise the caches used here including setting up the random
* number generator
*/
void MonteCarloAbsorption::initCaches()
{
if( !m_randGen )
{
m_randGen = new Kernel::MersenneTwister;
int seedValue = getProperty("SeedValue");
m_randGen->setSeed(seedValue);
}
m_samplePos = m_inputWS->getInstrument()->getSample()->getPos();
m_sourcePos = m_inputWS->getInstrument()->getSource()->getPos();
BoundingBox box(m_sampleShape->getBoundingBox());
if( m_container )
{
BoundingBox envBox;
m_container->getBoundingBox(envBox);
box.grow(envBox);
}
//Save the dimensions for quicker calculations later
m_bbox_length = box.xMax() - box.xMin();
m_bbox_halflength = 0.5 * m_bbox_length; m_bbox_width = box.yMax() - box.yMin();
m_bbox_halfwidth = 0.5 * m_bbox_width;
m_bbox_height = box.zMax() - box.zMin();
m_bbox_halfheight = 0.5 * m_bbox_height;
}
}
}