Mantid
Loading...
Searching...
No Matches
CalculateCarpenterSampleCorrection.cpp
Go to the documentation of this file.
1// Mantid Repository : https://github.com/mantidproject/mantid
2//
3// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
4// NScD Oak Ridge National Laboratory, European Spallation Source,
5// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
6// SPDX - License - Identifier: GPL - 3.0 +
10#include "MantidAPI/Sample.h"
19
20#include <stdexcept>
21
22namespace Mantid::Algorithms {
23DECLARE_ALGORITHM(CalculateCarpenterSampleCorrection) // Register the class
24 // into the algorithm
25 // factory
26using namespace DataObjects;
27using namespace Kernel;
28using namespace API;
31using Mantid::HistogramData::HistogramY;
32using Mantid::HistogramData::Points;
33using std::vector;
34using namespace Mantid::PhysicalConstants;
35using namespace Geometry;
36
37// Constants required internally only, so make them static. These are
38// Chebyshev expansion coefficients copied directly from Carpenter 1969 Table 1
39namespace { // anonymous
40static const double CHEBYSHEV[] = {
41 // l= 0 1 2 3 4 5 // (m,n)
42 0.730284, -0.249987, 0.019448, -0.000006, 0.000249, -0.000004, // (1,1)
43 0.848859, -0.452690, 0.056557, -0.000009, 0.000000, -0.000006, // (1,2)
44 1.133129, -0.749962, 0.118245, -0.000018, -0.001345, -0.000012, // (1,3)
45 1.641112, -1.241639, 0.226247, -0.000045, -0.004821, -0.000030, // (1,4)
46 0.848859, -0.452690, 0.056557, -0.000009, 0.000000, -0.000006, // (2,1)
47 1.000006, -0.821100, 0.166645, -0.012096, 0.000008, -0.000126, // (2,2)
48 1.358113, -1.358076, 0.348199, -0.038817, 0.000022, -0.000021, // (2,3)
49 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, // (2,4)
50 1.133129, -0.749962, 0.118245, -0.000018, -0.001345, -0.000012, // (3,1)
51 1.358113, -1.358076, 0.348199, -0.038817, 0.000022, -0.000021, // (3,2)
52 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, // (3,3)
53 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, // (3,4)
54 1.641112, -1.241639, 0.226247, -0.000045, -0.004821, -0.000030, // (4,1)
55 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, // (4,2)
56 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, // (4,3)
57 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 // (4,4)
58};
59
60static const int Z_size = 36; // Caution, this must be updated if the
61 // algorithm is changed to use a different
62 // size Z array.
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};
67
68static const double LAMBDA_REF = 1.81;
69// Badly named constants, no explanation of the origin of these
70// values. They appear to be used when calculating the multiple
71// scattering correction factor.
72static const double COEFF4 = 1.1967;
73static const double COEFF5 = -0.8667;
74} // namespace
75
76const std::string CalculateCarpenterSampleCorrection::name() const { return "CalculateCarpenterSampleCorrection"; }
77
79
81 return "CorrectionFunctions\\AbsorptionCorrections";
82}
83
88 // The input workspace must have an instrument and units of wavelength
89 auto wsValidator = std::make_shared<CompositeValidator>();
90 wsValidator->add<WorkspaceUnitValidator>("Wavelength");
91 wsValidator->add<InstrumentValidator>();
92
93 declareProperty(
94 std::make_unique<WorkspaceProperty<MatrixWorkspace>>("InputWorkspace", "", Direction::Input, wsValidator),
95 "The name of the input workspace.");
96 declareProperty(
97 std::make_unique<WorkspaceProperty<API::WorkspaceGroup>>("OutputWorkspaceBaseName", "", Direction::Output),
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.",
114}
115
120 // common information
121 MatrixWorkspace_sptr inputWksp = getProperty("InputWorkspace");
122 double radius = getProperty("CylinderSampleRadius");
123 double coeff1 = getProperty("AttenuationXSection");
124 double coeff2 = getProperty("SampleNumberDensity");
125 double coeff3 = getProperty("ScatteringXSection");
126 const bool absOn = getProperty("Absorption");
127 const bool msOn = getProperty("MultipleScattering");
128 const Material &sampleMaterial = inputWksp->sample().getMaterial();
129 if (sampleMaterial.totalScatterXSection() != 0.0) {
130 g_log.information() << "Using material \"" << sampleMaterial.name() << "\" from workspace\n";
131 if (std::abs(coeff1 - 2.8) < std::numeric_limits<double>::epsilon())
132 coeff1 = sampleMaterial.absorbXSection(LAMBDA_REF) / LAMBDA_REF;
133 if ((std::abs(coeff2 - 0.0721) < std::numeric_limits<double>::epsilon()) &&
134 (!isEmpty(sampleMaterial.numberDensity())))
135 coeff2 = sampleMaterial.numberDensity();
136 if (std::abs(coeff3 - 5.1) < std::numeric_limits<double>::epsilon())
137 coeff3 = sampleMaterial.totalScatterXSection();
138 } else // Save input in Sample with wrong atomic number and name
139 {
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);
144 }
145 g_log.debug() << "radius=" << radius << " coeff1=" << coeff1 << " coeff2=" << coeff2 << " coeff3=" << coeff3 << "\n";
146
147 // geometry stuff
148 const auto NUM_HIST = static_cast<int64_t>(inputWksp->getNumberHistograms());
149 Instrument_const_sptr instrument = inputWksp->getInstrument();
150 if (instrument == nullptr)
151 throw std::runtime_error("Failed to find instrument attached to InputWorkspace");
152 IComponent_const_sptr source = instrument->getSource();
153 IComponent_const_sptr sample = instrument->getSample();
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");
158
159 // Initialize progress reporting.
160 Progress prog(this, 0.0, 1.0, NUM_HIST);
161
162 EventWorkspace_sptr inputWkspEvent = std::dynamic_pointer_cast<EventWorkspace>(inputWksp);
163
164 // Create the new correction workspaces
165 MatrixWorkspace_sptr absWksp = createOutputWorkspace(inputWksp, "Attenuation factor");
166 MatrixWorkspace_sptr msWksp = createOutputWorkspace(inputWksp, "Multiple scattering factor");
167
168 // now do the correction
169 const auto &spectrumInfo = inputWksp->spectrumInfo();
170 PARALLEL_FOR_IF(Kernel::threadSafe(*absWksp, *msWksp))
171 for (int64_t index = 0; index < NUM_HIST; ++index) {
173 if (!spectrumInfo.hasDetectors(index))
174 throw std::runtime_error("Failed to find detector");
175 if (spectrumInfo.isMasked(index))
176 continue;
177 const double tth_rad = spectrumInfo.twoTheta(index);
178
179 // absorption
180 if (absOn) {
181 absWksp->setSharedX(index, inputWksp->sharedX(index));
182 const auto lambdas = inputWksp->points(index);
183 auto &y = absWksp->mutableY(index);
184 calculate_abs_correction(tth_rad, radius, coeff1, coeff2, coeff3, lambdas, y);
185 }
186
187 // multiple scattering
188 if (msOn) {
189 msWksp->setSharedX(index, inputWksp->sharedX(index));
190 const auto lambdas = inputWksp->points(index);
191 auto &y = msWksp->mutableY(index);
192 calculate_ms_correction(tth_rad, radius, coeff1, coeff2, coeff3, lambdas, y);
193 }
194
195 prog.report();
197 }
199
200 absWksp->setDistribution(true);
201 absWksp->setYUnit("");
202 absWksp->setYUnitLabel("Attenuation factor");
203
204 msWksp->setDistribution(true);
205 msWksp->setYUnit("");
206 msWksp->setYUnitLabel("Multiple scattering factor");
207
208 // Group and output workspaces we calculated
209 const std::string group_prefix = getPropertyValue("OutputWorkspaceBaseName");
210 auto outputGroup = std::make_shared<API::WorkspaceGroup>();
211 if (absOn) {
212 absWksp = setUncertainties(absWksp);
213 std::string ws_name = group_prefix + std::string("_abs");
214 AnalysisDataService::Instance().addOrReplace(ws_name, absWksp);
215 outputGroup->addWorkspace(absWksp);
216 } else {
217 deleteWorkspace(absWksp);
218 }
219
220 if (msOn) {
221 msWksp = setUncertainties(msWksp);
222 std::string ws_name = group_prefix + std::string("_ms");
223 AnalysisDataService::Instance().addOrReplace(ws_name, msWksp);
224 outputGroup->addWorkspace(msWksp);
225 } else {
226 deleteWorkspace(msWksp);
227 }
228
229 setProperty("OutputWorkspaceBaseName", outputGroup);
230}
231
232namespace { // anonymous namespace
233// Set up the Z table for the specified two theta angle (in degrees).
234vector<double> createZ(const double angle_rad) {
235 vector<double> Z(Z_initial, Z_initial + Z_size);
236
237 const double theta_rad = angle_rad * .5;
238 int l, J;
239 double sum;
240
241 for (int i = 1; i <= 4; i++) {
242 for (int j = 1; j <= 4; j++) {
243 int iplusj = i + j;
244 if (iplusj <= 5) {
245 l = 0;
246 J = 1 + l + 6 * (i - 1) + 6 * 4 * (j - 1);
247 sum = CHEBYSHEV[J - 1];
248
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);
252 }
253 J = 1 + i + 6 * j;
254 Z[J - 1] = sum;
255 }
256 }
257 }
258 return Z;
259}
260
264double AttFac(const double sigir, const double sigsr, const vector<double> &Z) {
265 double facti = 1.0;
266 double att = 0.0;
267
268 for (size_t i = 0; i <= 5; i++) {
269 double facts = 1.0;
270 for (size_t j = 0; j <= 5; j++) {
271 if (i + j <= 5) {
272 size_t J = 1 + i + 6 * j; // TODO J defined in terms of j?
273 att = att + Z[J - 1] * facts * facti;
274 facts = -facts * sigsr / static_cast<double>(j + 1);
275 }
276 }
277 facti = -facti * sigir / static_cast<double>(i + 1);
278 }
279 return att;
280}
281
282double calculate_abs_factor(const double radius, const double Q2, const double sigsct, const vector<double> &Z,
283 const double wavelength) {
284
285 const double sigabs = Q2 * wavelength;
286 const double sigir = (sigabs + sigsct) * radius;
292 const double sigsr = sigir;
293
294 return AttFac(sigir, sigsr, Z);
295}
296
297double calculate_ms_factor(const double radius, const double Q2, const double sigsct, const vector<double> &Z,
298 const double wavelength) {
299
300 const double sigabs = Q2 * wavelength;
301 const double sigir = (sigabs + sigsct) * radius;
307 const double sigsr = sigir;
308
309 const double delta = COEFF4 * sigir + COEFF5 * sigir * sigir;
310 const double deltp = (delta * sigsct) / (sigsct + sigabs);
311
312 double temp = AttFac(sigir, sigsr, Z);
313 return (deltp / temp);
314}
315
316} // namespace
317
333 const double coeff1, const double coeff2,
334 const double coeff3, const Points &wavelength,
335 HistogramY &y_val) {
336
337 const size_t NUM_Y = y_val.size();
338 bool is_histogram = false;
339 if (wavelength.size() == NUM_Y + 1)
340 is_histogram = true;
341 else if (wavelength.size() == NUM_Y)
342 is_histogram = false;
343 else
344 throw std::runtime_error("Data is neither historgram or density");
345
346 // initialize Z array for this angle
347 vector<double> Z = createZ(angle_deg);
348
349 const double Q2 = coeff1 * coeff2;
350 const double sigsct = coeff2 * coeff3;
351
352 for (size_t j = 0; j < NUM_Y; j++) {
353 double wl_val = wavelength[j];
354 if (is_histogram) // average with next value
355 wl_val = .5 * (wl_val + wavelength[j + 1]);
356
357 y_val[j] = calculate_abs_factor(radius, Q2, sigsct, Z, wl_val);
358 }
359}
360
362 const double coeff1, const double coeff2,
363 const double coeff3, const Points &wavelength,
364 HistogramY &y_val) {
365
366 const size_t NUM_Y = y_val.size();
367 bool is_histogram = false;
368 if (wavelength.size() == NUM_Y + 1)
369 is_histogram = true;
370 else if (wavelength.size() == NUM_Y)
371 is_histogram = false;
372 else
373 throw std::runtime_error("Data is neither historgram or density");
374
375 // initialize Z array for this angle
376 vector<double> Z = createZ(angle_deg);
377
378 const double Q2 = coeff1 * coeff2;
379 const double sigsct = coeff2 * coeff3;
380
381 for (size_t j = 0; j < NUM_Y; j++) {
382 double wl_val = wavelength[j];
383 if (is_histogram) // average with next value
384 wl_val = .5 * (wl_val + wavelength[j + 1]);
385
386 y_val[j] = calculate_ms_factor(radius, Q2, sigsct, Z, wl_val);
387 }
388}
389
391 const std::string &ylabel) const {
392 MatrixWorkspace_sptr outputWS = create<HistoWorkspace>(*inputWksp);
393 // The algorithm computes the signal values at bin centres so they should
394 // be treated as a distribution
395 outputWS->setDistribution(true);
396 outputWS->setYUnit("");
397 outputWS->setYUnitLabel(ylabel);
398 return outputWS;
399}
400
402 auto alg = this->createChildAlgorithm("SetUncertainties");
403 alg->initialize();
404 alg->setProperty("InputWorkspace", workspace);
405 alg->execute();
406 return alg->getProperty("OutputWorkspace");
407}
408
410 auto alg = this->createChildAlgorithm("DeleteWorkspace");
411 alg->initialize();
412 alg->setChild(true);
413 alg->setLogging(false);
414 alg->setProperty("Workspace", workspace);
415 alg->execute();
416}
417
418} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
IPeaksWorkspace_sptr workspace
Definition: IndexPeaks.cpp:114
std::map< DeltaEMode::Type, std::string > index
Definition: DeltaEMode.cpp:19
#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...
Definition: MultiThreaded.h:94
#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.
double radius
Definition: Rasterize.cpp:31
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.
Definition: Progress.h:25
A property class for workspaces.
A validator which checks that the unit of the workspace referred to by a WorkspaceProperty is the exp...
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.
Definition: Logger.cpp:114
void information(const std::string &msg)
Logs at information level.
Definition: Logger.cpp:105
A material is defined as being composed of a given element, defined as a PhysicalConstants::NeutronAt...
Definition: Material.h:50
double numberDensity() const
Get the number density.
Definition: Material.cpp:189
double absorbXSection(const double lambda=PhysicalConstants::NeutronAtom::ReferenceLambda) const
Get the absorption cross section at a given wavelength in barns.
Definition: Material.cpp:260
const std::string & name() const
Returns the name of the material.
Definition: Material.cpp:181
double totalScatterXSection() const
Return the total scattering cross section for a given wavelength in barns.
Definition: Material.cpp:252
void report()
Increments the loop counter by 1, then sends the progress notification on behalf of its algorithm.
Definition: ProgressBase.h:51
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.
Definition: IComponent.h:161
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.
Definition: MultiThreaded.h:22
A namespace containing physical constants that are required by algorithms and unit routines.
Definition: Atom.h:14
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54
Structure to store neutronic scattering information for the various elements.
Definition: NeutronAtom.h:22