26using namespace DataObjects;
27using namespace Kernel;
31using Mantid::HistogramData::HistogramY;
32using Mantid::HistogramData::Points;
35using namespace Geometry;
40static const double CHEBYSHEV[] = {
42 0.730284, -0.249987, 0.019448, -0.000006, 0.000249, -0.000004,
43 0.848859, -0.452690, 0.056557, -0.000009, 0.000000, -0.000006,
44 1.133129, -0.749962, 0.118245, -0.000018, -0.001345, -0.000012,
45 1.641112, -1.241639, 0.226247, -0.000045, -0.004821, -0.000030,
46 0.848859, -0.452690, 0.056557, -0.000009, 0.000000, -0.000006,
47 1.000006, -0.821100, 0.166645, -0.012096, 0.000008, -0.000126,
48 1.358113, -1.358076, 0.348199, -0.038817, 0.000022, -0.000021,
49 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
50 1.133129, -0.749962, 0.118245, -0.000018, -0.001345, -0.000012,
51 1.358113, -1.358076, 0.348199, -0.038817, 0.000022, -0.000021,
52 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
53 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
54 1.641112, -1.241639, 0.226247, -0.000045, -0.004821, -0.000030,
55 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
56 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
57 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
60static const int Z_size = 36;
63static const double Z_initial[] = {
64 1.0, 0.8488263632, 1.0, 1.358122181, 2.0, 3.104279270, 0.8488263632, 0.0, 0.0, 0.0, 0.0, 0.0,
65 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.358122181, 0.0, 0.0, 0.0, 0.0, 0.0,
66 2.0, 0.0, 0.0, 0.0, 0.0, 0.0, 3.104279270, 0.0, 0.0, 0.0, 0.0, 0.0};
68static const double LAMBDA_REF = 1.81;
72static const double COEFF4 = 1.1967;
73static const double COEFF5 = -0.8667;
81 return "CorrectionFunctions\\AbsorptionCorrections";
89 auto wsValidator = std::make_shared<CompositeValidator>();
95 "The name of the input workspace.");
98 "Basename of the output workspace group for corrections."
99 "Absorption suffix = '_abs'. "
100 "Multiple Scattering suffix = '_ms'. ");
101 declareProperty(
"AttenuationXSection", 2.8,
102 "Coefficient 1, absorption cross "
103 "section / 1.81 if not set with "
104 "SetSampleMaterial");
105 declareProperty(
"ScatteringXSection", 5.1,
106 "Coefficient 3, total scattering "
107 "cross section if not set with "
108 "SetSampleMaterial");
109 declareProperty(
"SampleNumberDensity", 0.0721,
"Coefficient 2, density if not set with SetSampleMaterial");
110 declareProperty(
"CylinderSampleRadius", 0.3175,
"Sample radius, in cm");
111 declareProperty(
"Absorption",
true,
"If True then calculates the absorption correction.",
Direction::Input);
112 declareProperty(
"MultipleScattering",
true,
"If True then calculates the multiple scattering correction.",
123 double coeff1 =
getProperty(
"AttenuationXSection");
124 double coeff2 =
getProperty(
"SampleNumberDensity");
127 const bool msOn =
getProperty(
"MultipleScattering");
128 const Material &sampleMaterial = inputWksp->sample().getMaterial();
131 if (std::abs(coeff1 - 2.8) < std::numeric_limits<double>::epsilon())
133 if ((std::abs(coeff2 - 0.0721) < std::numeric_limits<double>::epsilon()) &&
136 if (std::abs(coeff3 - 5.1) < std::numeric_limits<double>::epsilon())
140 NeutronAtom neutron(0, 0, 0.0, 0.0, coeff3, 0.0, coeff3, coeff1);
141 auto shape = std::shared_ptr<IObject>(
142 inputWksp->sample().getShape().cloneWithMaterial(
Material(
"SetInMultipleScattering", neutron, coeff2)));
143 inputWksp->mutableSample().setShape(shape);
145 g_log.
debug() <<
"radius=" <<
radius <<
" coeff1=" << coeff1 <<
" coeff2=" << coeff2 <<
" coeff3=" << coeff3 <<
"\n";
148 const auto NUM_HIST =
static_cast<int64_t
>(inputWksp->getNumberHistograms());
150 if (instrument ==
nullptr)
151 throw std::runtime_error(
"Failed to find instrument attached to InputWorkspace");
154 if (source ==
nullptr)
155 throw std::runtime_error(
"Failed to find source in the instrument for InputWorkspace");
156 if (sample ==
nullptr)
157 throw std::runtime_error(
"Failed to find sample in the instrument for InputWorkspace");
160 Progress prog(
this, 0.0, 1.0, NUM_HIST);
169 const auto &spectrumInfo = inputWksp->spectrumInfo();
173 if (!spectrumInfo.hasDetectors(
index))
174 throw std::runtime_error(
"Failed to find detector");
175 if (spectrumInfo.isMasked(
index))
177 const double tth_rad = spectrumInfo.twoTheta(
index);
181 absWksp->setSharedX(
index, inputWksp->sharedX(
index));
182 const auto lambdas = inputWksp->points(
index);
183 auto &
y = absWksp->mutableY(
index);
189 msWksp->setSharedX(
index, inputWksp->sharedX(
index));
190 const auto lambdas = inputWksp->points(
index);
191 auto &
y = msWksp->mutableY(
index);
200 absWksp->setDistribution(
true);
201 absWksp->setYUnit(
"");
202 absWksp->setYUnitLabel(
"Attenuation factor");
204 msWksp->setDistribution(
true);
205 msWksp->setYUnit(
"");
206 msWksp->setYUnitLabel(
"Multiple scattering factor");
209 const std::string group_prefix =
getPropertyValue(
"OutputWorkspaceBaseName");
210 auto outputGroup = std::make_shared<API::WorkspaceGroup>();
213 std::string ws_name = group_prefix + std::string(
"_abs");
215 outputGroup->addWorkspace(absWksp);
222 std::string ws_name = group_prefix + std::string(
"_ms");
224 outputGroup->addWorkspace(msWksp);
229 setProperty(
"OutputWorkspaceBaseName", outputGroup);
234vector<double> createZ(
const double angle_rad) {
235 vector<double>
Z(Z_initial, Z_initial + Z_size);
237 const double theta_rad = angle_rad * .5;
241 for (
int i = 1; i <= 4; i++) {
242 for (
int j = 1; j <= 4; j++) {
246 J = 1 + l + 6 * (i - 1) + 6 * 4 * (j - 1);
247 sum = CHEBYSHEV[J - 1];
249 for (l = 1; l <= 5; l++) {
250 J = 1 + l + 6 * (i - 1) + 6 * 4 * (j - 1);
251 sum = sum + CHEBYSHEV[J - 1] * cos(l * theta_rad);
264double AttFac(
const double sigir,
const double sigsr,
const vector<double> &
Z) {
268 for (
size_t i = 0; i <= 5; i++) {
270 for (
size_t j = 0; j <= 5; j++) {
272 size_t J = 1 + i + 6 * j;
273 att = att +
Z[J - 1] * facts * facti;
274 facts = -facts * sigsr /
static_cast<double>(j + 1);
277 facti = -facti * sigir /
static_cast<double>(i + 1);
282double calculate_abs_factor(
const double radius,
const double Q2,
const double sigsct,
const vector<double> &
Z,
283 const double wavelength) {
285 const double sigabs = Q2 * wavelength;
286 const double sigir = (sigabs + sigsct) *
radius;
292 const double sigsr = sigir;
294 return AttFac(sigir, sigsr,
Z);
297double calculate_ms_factor(
const double radius,
const double Q2,
const double sigsct,
const vector<double> &
Z,
298 const double wavelength) {
300 const double sigabs = Q2 * wavelength;
301 const double sigir = (sigabs + sigsct) *
radius;
307 const double sigsr = sigir;
309 const double delta = COEFF4 * sigir + COEFF5 * sigir * sigir;
310 const double deltp = (
delta * sigsct) / (sigsct + sigabs);
312 double temp = AttFac(sigir, sigsr,
Z);
313 return (deltp / temp);
333 const double coeff1,
const double coeff2,
334 const double coeff3,
const Points &wavelength,
337 const size_t NUM_Y = y_val.size();
338 bool is_histogram =
false;
339 if (wavelength.size() == NUM_Y + 1)
341 else if (wavelength.size() == NUM_Y)
342 is_histogram =
false;
344 throw std::runtime_error(
"Data is neither historgram or density");
347 vector<double>
Z = createZ(angle_deg);
349 const double Q2 = coeff1 * coeff2;
350 const double sigsct = coeff2 * coeff3;
352 for (
size_t j = 0; j < NUM_Y; j++) {
353 double wl_val = wavelength[j];
355 wl_val = .5 * (wl_val + wavelength[j + 1]);
357 y_val[j] = calculate_abs_factor(
radius, Q2, sigsct,
Z, wl_val);
362 const double coeff1,
const double coeff2,
363 const double coeff3,
const Points &wavelength,
366 const size_t NUM_Y = y_val.size();
367 bool is_histogram =
false;
368 if (wavelength.size() == NUM_Y + 1)
370 else if (wavelength.size() == NUM_Y)
371 is_histogram =
false;
373 throw std::runtime_error(
"Data is neither historgram or density");
376 vector<double>
Z = createZ(angle_deg);
378 const double Q2 = coeff1 * coeff2;
379 const double sigsct = coeff2 * coeff3;
381 for (
size_t j = 0; j < NUM_Y; j++) {
382 double wl_val = wavelength[j];
384 wl_val = .5 * (wl_val + wavelength[j + 1]);
386 y_val[j] = calculate_ms_factor(
radius, Q2, sigsct,
Z, wl_val);
391 const std::string &ylabel)
const {
395 outputWS->setDistribution(
true);
396 outputWS->setYUnit(
"");
397 outputWS->setYUnitLabel(ylabel);
404 alg->setProperty(
"InputWorkspace",
workspace);
406 return alg->getProperty(
"OutputWorkspace");
413 alg->setLogging(
false);
414 alg->setProperty(
"Workspace",
workspace);
#define DECLARE_ALGORITHM(classname)
IPeaksWorkspace_sptr workspace
std::map< DeltaEMode::Type, std::string > index
#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.
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) override
Create a Child Algorithm.
Kernel::IPropertyManager::TypedValue getProperty(const std::string &name) const override
Get the property held by this object.
std::string getPropertyValue(const std::string &name) const override
Get the property held by this object.
A validator which checks that a workspace has a valid instrument.
Helper class for reporting progress from algorithms.
A property class for workspaces.
A validator which checks that the unit of the workspace referred to by a WorkspaceProperty is the exp...
void exec() override
Execute the algorithm.
int version() const override
Algorithm's version for identification overriding a virtual method.
const std::string name() const override
Algorithm's name for identification overriding a virtual method.
void calculate_abs_correction(const double angle_deg, const double radius, const double coeff1, const double coeff2, const double coeff3, const HistogramData::Points &wavelength, HistogramData::HistogramY &y_val)
CalculateCarpenterSampleCorrection correction calculation.
void calculate_ms_correction(const double angle_deg, const double radius, const double coeff1, const double coeff2, const double coeff3, const HistogramData::Points &wavelength, HistogramData::HistogramY &y_val)
const std::string category() const override
Algorithm's category for identification overriding a virtual method.
void deleteWorkspace(const API::MatrixWorkspace_sptr &workspace)
void init() override
Initialize the properties to default values.
API::MatrixWorkspace_sptr createOutputWorkspace(const API::MatrixWorkspace_sptr &inputWS, const std::string &) const
API::MatrixWorkspace_sptr setUncertainties(const API::MatrixWorkspace_sptr &workspace)
This class is intended to fulfill the design specified in <https://github.com/mantidproject/documents...
void debug(const std::string &msg)
Logs at debug level.
void information(const std::string &msg)
Logs at information level.
A material is defined as being composed of a given element, defined as a PhysicalConstants::NeutronAt...
double numberDensity() const
Get the number density.
double absorbXSection(const double lambda=PhysicalConstants::NeutronAtom::ReferenceLambda) const
Get the absorption cross section at a given wavelength in barns.
const std::string & name() const
Returns the name of the material.
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...
Kernel::Logger g_log("ExperimentInfo")
static logger object
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
std::shared_ptr< EventWorkspace > EventWorkspace_sptr
shared pointer to the EventWorkspace class
std::shared_ptr< const IComponent > IComponent_const_sptr
Typdef of a shared pointer to a const IComponent.
std::shared_ptr< const Instrument > Instrument_const_sptr
Shared pointer to an const instrument object.
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.
A namespace containing physical constants that are required by algorithms and unit routines.
@ Input
An input workspace.
@ Output
An output workspace.
Structure to store neutronic scattering information for the various elements.