22#include "MantidHistogramData/Interpolate.h"
34using namespace Geometry;
35using HistogramData::interpolateLinearInplace;
36using namespace Kernel;
41constexpr size_t MAX_INTEGRATION_LENGTH{1000};
43inline size_t findMiddle(
const size_t start,
const size_t stop) {
44 auto half =
static_cast<size_t>(floor(.5 * (
static_cast<double>(stop - start))));
51 :
API::
Algorithm(), m_inputWS(), m_sampleObject(nullptr), m_containerObject(nullptr), m_sampleL1s(),
52 m_sample_containerL1s(), m_sampleElementVolumes(), m_sampleElementPositions(), m_numSampleVolumeElements(0),
53 m_sampleVolume(0.0), m_containerL1s(), m_container_sampleL1s(), m_containerElementVolumes(),
54 m_containerElementPositions(), m_numContainerVolumeElements(0), m_containerVolume(0.0),
55 m_ampleLinearCoefTotScatt(0), m_containerLinearCoefTotScatt(0), m_num_lambda(0), m_xStep(0),
56 m_cubeSideSample(0.0), m_cubeSideContainer(0.0) {}
61 auto wsValidator = std::make_shared<CompositeValidator>();
66 "The X values for the input workspace must be in units of wavelength");
68 "Output workspace name");
70 auto positiveInt = std::make_shared<BoundedValidator<int64_t>>();
71 positiveInt->setLower(1);
73 "The number of wavelength points for which the numerical integral is\n"
74 "calculated (default: all points)");
76 auto moreThanZero = std::make_shared<BoundedValidator<double>>();
77 moreThanZero->setLower(0.001);
78 declareProperty(
"ElementSize", 1.0, moreThanZero,
"The size of one side of an integration element cube in mm");
81 "The size of one side of an integration element cube in mm for container."
82 "Default to be the same as ElementSize.");
86 std::map<std::string, std::string> result;
90 const auto &sample = wksp->sample();
91 if (sample.hasEnvironment()) {
92 const auto numComponents = sample.getEnvironment().nelements();
94 if (numComponents == 0) {
95 result[
"InputWorkspace"] =
"Sample does not have a container defined";
98 result[
"InputWorkspace"] =
"Sample does not have a container defined";
116 ass->setDistribution(
true);
118 ass->setYUnitLabel(
"Attenuation factor");
122 assc->setDistribution(
true);
124 assc->setYUnitLabel(
"Attenuation factor");
127 acc->setDistribution(
true);
129 acc->setYUnitLabel(
"Attenuation factor");
133 acsc->setDistribution(
true);
135 acsc->setYUnitLabel(
"Attenuation factor");
139 const auto numHists =
static_cast<int64_t
>(
m_inputWS->getNumberHistograms());
140 const auto specSize =
static_cast<int64_t
>(
m_inputWS->blocksize());
150 std::ostringstream message;
151 message <<
"Numerical integration performed every " <<
m_xStep <<
" wavelength points";
158 throw std::runtime_error(
"Failed to define any initial scattering gauge volume for geometry");
161 const auto &spectrumInfo =
m_inputWS->spectrumInfo();
162 Progress prog(
this, 0.0, 1.0, numHists);
165 for (int64_t i = 0; i < int64_t(numHists); ++i) {
168 ass->setSharedX(i,
m_inputWS->sharedX(i));
169 assc->setSharedX(i,
m_inputWS->sharedX(i));
170 acc->setSharedX(i,
m_inputWS->sharedX(i));
171 acsc->setSharedX(i,
m_inputWS->sharedX(i));
173 if (!spectrumInfo.hasDetectors(i)) {
174 g_log.
information() <<
"Spectrum " << i <<
" does not have a detector defined for it\n";
177 const auto &det = spectrumInfo.detector(i);
192 calculateDistances(det, sample_L2s, sample_container_L2s, container_L2s, container_sample_L2s);
194 const auto wavelengths =
m_inputWS->points(i);
200 auto &assY = ass->mutableY(i);
201 auto &asscY = assc->mutableY(i);
202 auto &accY = acc->mutableY(i);
203 auto &acscY = acsc->mutableY(i);
206 for (int64_t j = 0; j < specSize; j = j +
m_xStep) {
207 double integral = 0.0;
208 double crossIntegral = 0.0;
226 if (
m_xStep > 1 && j +
m_xStep >= specSize && j + 1 != specSize) {
234 auto histnew = ass->histogram(i);
235 interpolateLinearInplace(histnew,
m_xStep);
236 ass->setHistogram(i, histnew);
237 auto histnew2 = assc->histogram(i);
238 interpolateLinearInplace(histnew2,
m_xStep);
239 assc->setHistogram(i, histnew2);
240 auto histnew3 = acc->histogram(i);
241 interpolateLinearInplace(histnew3,
m_xStep);
242 acc->setHistogram(i, histnew3);
243 auto histnew4 = acsc->histogram(i);
244 interpolateLinearInplace(histnew4,
m_xStep);
245 acsc->setHistogram(i, histnew4);
256 const std::string outWSName =
getProperty(
"OutputWorkspace");
258 std::vector<std::string> names;
259 names.emplace_back(outWSName +
"_ass");
261 names.emplace_back(outWSName +
"_assc");
263 names.emplace_back(outWSName +
"_acc");
265 names.emplace_back(outWSName +
"_acsc");
270 group->setProperty(
"InputWorkspaces", names);
271 group->setProperty(
"OutputWorkspace", outWSName);
285 if (
m_inputWS->run().hasProperty(
"GaugeVolume")) {
291 if (raster.l1.size() == 0)
292 throw std::runtime_error(
"Failed to rasterize shape");
301 if (
m_inputWS->run().hasProperty(
"GaugeVolume")) {
307 if (raster.l1.size() == 0)
308 throw std::runtime_error(
"Failed to rasterize shape");
335 g_log.
information(
"Calculating scattering within the gauge volume defined on "
336 "the input workspace");
339 std::shared_ptr<const Geometry::IObject> volume =
349 const auto &sampleObj =
m_inputWS->sample();
376 const std::string mess(
"No shape has been defined for the sample in the input workspace");
378 throw std::invalid_argument(mess);
382 const std::string mess(
"No shape has been defined for the container in the input workspace");
384 throw std::invalid_argument(mess);
390 std::vector<double> &sample_container_L2s,
391 std::vector<double> &container_L2s,
392 std::vector<double> &container_sample_L2s)
const {
394 if (detector.
nDets() > 1) {
398 detector.
getPhi() * 180.0 / M_PI);
443 const double linearCoefAbs,
const double linearCoefTotScatt,
444 const std::vector<double> &elementVolumes,
445 const std::vector<double> &L1s,
const std::vector<double> &L2s,
446 const double linearCoefAbs2,
const double linearCoefTotScatt2,
447 const std::vector<double> &L1s2,
const std::vector<double> &L2s2,
448 const size_t startIndex,
const size_t endIndex)
const {
449 if (endIndex - startIndex > MAX_INTEGRATION_LENGTH) {
450 size_t middle = findMiddle(startIndex, endIndex);
452 doIntegration(integral, crossIntegral, linearCoefAbs, linearCoefTotScatt, elementVolumes, L1s, L2s, linearCoefAbs2,
453 linearCoefTotScatt2, L1s2, L2s2, startIndex, middle);
454 doIntegration(integral, crossIntegral, linearCoefAbs, linearCoefTotScatt, elementVolumes, L1s, L2s, linearCoefAbs2,
455 linearCoefTotScatt2, L1s2, L2s2, middle, endIndex);
459 for (
size_t i = startIndex; i < endIndex; ++i) {
460 double exponent = (linearCoefAbs + linearCoefTotScatt) * (L1s[i] + L2s[i]);
461 integral += (exp(exponent) * (elementVolumes[i]));
462 exponent += (linearCoefAbs2 + linearCoefTotScatt2) * (L1s2[i] + L2s2[i]);
463 crossIntegral += (exp(exponent) * (elementVolumes[i]));
#define DECLARE_ALGORITHM(classname)
#define PARALLEL_START_INTERRUPT_REGION
Begins a block to skip processing is the algorithm has been interupted Note the end of the block if n...
#define PARALLEL_END_INTERRUPT_REGION
Ends a block to skip processing is the algorithm has been interupted Note the start of the block if n...
#define PARALLEL_FOR_IF(condition)
Empty definitions - to enable set your complier to enable openMP.
#define PARALLEL_CHECK_INTERRUPT_REGION
Adds a check after a Parallel region to see if it was interupted.
Base class from which all concrete algorithm classes should be derived.
void declareProperty(std::unique_ptr< Kernel::Property > p, const std::string &doc="") override
Add a property to the list of managed properties.
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
virtual std::shared_ptr< Algorithm > createChildAlgorithm(const std::string &name, const double startProgress=-1., const double endProgress=-1., const bool enableLogging=true, const int &version=-1)
Create a Child Algorithm.
bool isDefault(const std::string &name) const
static bool isEmpty(const NumT toCheck)
checks that the value was not set by users, uses the value in empty double/int.
A validator which checks that a workspace has a valid instrument.
Helper class for reporting progress from algorithms.
This class stores information about the sample used in particular run.
const Geometry::IObject & getShape() const
Return the sample shape.
const Geometry::SampleEnvironment & getEnvironment() const
Get a reference to the sample's environment.
A property class for workspaces.
A validator which checks that the unit of the workspace referred to by a WorkspaceProperty is the exp...
API::MatrixWorkspace_sptr m_inputWS
A pointer to the input workspace.
void exec() override
Execution code.
void constructSample(API::Sample &sample)
Create the sample object using the Geometry classes, or use the existing one.
std::vector< double > m_sampleL1s
Cached sample L1 distances.
double m_containerLinearCoefTotScatt
The total scattering cross-section in 1/m for the container.
PaalmanPingsAbsorptionCorrection()
(Empty) Constructor
void initialiseCachedDistances()
Calculate the distances for L1 (for both self-absorption and absorption by other object) and element ...
std::vector< double > m_containerL1s
Cached container L1 distances.
size_t m_numContainerVolumeElements
The number of container volume elements.
std::vector< double > m_sample_containerL1s
Cached L1 distances through container hitting sample.
void retrieveBaseProperties()
Fetch the properties and set the appropriate member variables.
void calculateDistances(const Geometry::IDetector &detector, std::vector< double > &sample_L2s, std::vector< double > &sample_container_L2s, std::vector< double > &container_L2s, std::vector< double > &container_sample_L2s) const
Calculate the distances traversed by the neutrons within the sample.
std::vector< Kernel::V3D > m_sampleElementPositions
Cached sample element positions.
const Geometry::IObject * m_sampleObject
Local cache of sample object.
Kernel::Material m_material
std::shared_ptr< const Geometry::IObject > constructGaugeVolume()
Create the gague volume for the correction.
int64_t m_num_lambda
The number of points in wavelength, the rest is.
Kernel::Material m_containerMaterial
size_t m_numSampleVolumeElements
The number of sample volume elements.
std::map< std::string, std::string > validateInputs() override
Method checking errors on ALL the inputs, before execution.
const Geometry::IObject * m_containerObject
Local cache of container object.
std::vector< Kernel::V3D > m_containerElementPositions
Cached container element positions.
double m_cubeSideContainer
The length of the side of an element cube in m.
double m_cubeSideSample
The length of the side of an element cube in m.
Kernel::V3D m_beamDirection
The direction of the beam.
void doIntegration(double &integral, double &crossIntegral, const double linearCoefAbs, const double linearCoefTotScatt, const std::vector< double > &elementVolumes, const std::vector< double > &L1s, const std::vector< double > &L2s, const double linearCoefAbs2, const double linearCoefTotScatt2, const std::vector< double > &L1s2, const std::vector< double > &L2s2, const size_t startIndex, const size_t endIndex) const
Carries out the numerical integration over the sample for elastic instruments.
std::vector< double > m_containerElementVolumes
Cached container element volumes.
std::vector< double > m_container_sampleL1s
Cached L1 distances through sample hitting container.
int64_t m_xStep
interpolated linearly
void init() override
Initialisation code.
double m_ampleLinearCoefTotScatt
The total scattering cross-section in 1/m for the sample.
std::vector< double > m_sampleElementVolumes
Cached sample element volumes.
double m_sampleVolume
The total volume of the sample.
double m_containerVolume
The total volume of the container.
virtual Kernel::V3D getPos() const =0
Get the position of the IComponent. Tree structure is traverse through the.
Interface class for detector objects.
virtual double getPhi() const =0
Gives the phi of this detector object in radians.
virtual std::size_t nDets() const =0
Get the number of physical detectors this object represents.
virtual double getTwoTheta(const Kernel::V3D &observer, const Kernel::V3D &axis) const =0
Gives the angle of this detector object with respect to an axis.
virtual int interceptSurface(Geometry::Track &) const =0
virtual IObject * clone() const =0
virtual bool hasValidShape() const =0
const Container & getContainer() const
Class originally intended to be used with the DataHandling 'LoadInstrument' algorithm.
std::shared_ptr< CSGObject > createShape(Poco::XML::Element *pElem)
Creates a geometric object from a DOM-element-node pointing to an element whose child nodes contain t...
Defines a track as a start point and a direction.
void clearIntersectionResults()
Clear the current set of intersection results.
double totalDistInsideObject() const
Returns the sum of all the links distInsideObject in the track.
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void error(const std::string &msg)
Logs at error level.
void information(const std::string &msg)
Logs at information level.
double numberDensityEffective() const
Get the effective number density.
double linearAbsorpCoef(const double lambda=PhysicalConstants::NeutronAtom::ReferenceLambda) const
Returns the linear coefficient of absorption for the material in units of cm^-1 this should match the...
double totalScatterXSection() const
Return the total scattering cross section for a given wavelength in barns.
void report()
Increments the loop counter by 1, then sends the progress notification on behalf of its algorithm.
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
void spherical(const double R, const double theta, const double phi) noexcept
Sets the vector position based on spherical coordinates.
double norm() const noexcept
std::shared_ptr< WorkspaceGroup > WorkspaceGroup_sptr
shared pointer to Mantid::API::WorkspaceGroup
std::shared_ptr< const MatrixWorkspace > MatrixWorkspace_const_sptr
shared pointer to the matrix workspace base class (const version)
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
MANTID_GEOMETRY_DLL Raster calculate(const Kernel::V3D &beamDirection, const IObject &shape, const double cubeSizeInMetre)
std::enable_if< std::is_pointer< Arg >::value, bool >::type threadSafe(Arg workspace)
Thread-safety check Checks the workspace to ensure it is suitable for multithreaded access.
MANTID_KERNEL_DLL V3D normalize(V3D v)
Normalizes a V3D.
constexpr int EMPTY_INT() noexcept
Returns what we consider an "empty" integer within a property.
constexpr double EMPTY_DBL() noexcept
Returns what we consider an "empty" double within a property.
@ Input
An input workspace.
@ Output
An output workspace.