Mantid
Loading...
Searching...
No Matches
HRPDSlabCanAbsorption.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 +
9#include "MantidAPI/Sample.h"
15
16namespace Mantid::Algorithms {
17
18// Register the algorithm into the AlgorithmFactory
19DECLARE_ALGORITHM(HRPDSlabCanAbsorption)
20
21using namespace Kernel;
22using namespace API;
23using namespace Geometry;
24using namespace Mantid::PhysicalConstants;
25
27 declareProperty(std::make_unique<WorkspaceProperty<>>("InputWorkspace", "", Direction::Input));
28 declareProperty(std::make_unique<WorkspaceProperty<>>("OutputWorkspace", "", Direction::Output));
29
30 auto mustBePositive = std::make_shared<BoundedValidator<double>>();
31 mustBePositive->setLower(0.0);
32 declareProperty("SampleAttenuationXSection", EMPTY_DBL(), mustBePositive,
33 "The ABSORPTION cross-section for the sample material in "
34 "barns if not set with SetSampleMaterial");
35 declareProperty("SampleScatteringXSection", EMPTY_DBL(), mustBePositive,
36 "The scattering cross-section (coherent + incoherent) for "
37 "the sample material in barns if not set with "
38 "SetSampleMaterial");
39 declareProperty("SampleNumberDensity", EMPTY_DBL(), mustBePositive,
40 "The number density of the sample in number of atoms per "
41 "cubic angstrom if not set with SetSampleMaterial");
42 declareProperty("Thickness", 0.2, mustBePositive,
43 "The thickness of the sample in cm. Common values are 0.2, "
44 "0.5, 1.0, 1.5");
45
46 auto positiveInt = std::make_shared<BoundedValidator<int64_t>>();
47 positiveInt->setLower(1);
48 declareProperty("NumberOfWavelengthPoints", int64_t(EMPTY_INT()), positiveInt,
49 "The number of wavelength points for which the numerical integral is\n"
50 "calculated (default: all points)");
51
52 std::vector<std::string> exp_options{"Normal", "FastApprox"};
53 declareProperty("ExpMethod", "Normal", std::make_shared<StringListValidator>(exp_options),
54 "Select the method to use to calculate exponentials, normal or a\n"
55 "fast approximation (default: Normal)");
56
57 auto moreThanZero = std::make_shared<BoundedValidator<double>>();
58 moreThanZero->setLower(0.001);
59 declareProperty("ElementSize", 1.0, moreThanZero, "The size of one side of an integration element cube in mm");
60}
61
63 // First run the numerical absorption correction for the sample
65
66 /* Now do the corrections for the sample holder.
67 Using an analytical formula for a flat plate correction:
68 exp( -mt/cos(theta) )
69 Fullprof says it then divides by another cos(theta), but I
70 can't see where this comes from - it would make the transmission
71 higher at higher angle, which isn't intuitive.
72 In this equation, it's not obvious to me where theta should be
73 calculated from - I'm taking the centre of the plate (in any
74 case this would only make a very small difference)
75 */
76
77 /* The HRPD sample holders consist of a cuboid void for the sample of
78 18x23 mm and some thickness, 0.125mm vanadium foil windows front
79 and back, and an aluminium surround of total width 40mm, meaning
80 that there is 11mm of Al to be traversed en route to the 90 degree
81 bank.
82 */
83
84 const double vanWinThickness = 0.000125; // in m
85 const double vanRho = 0.07192 * 100; // in Angstroms-3 * 100
86 const double vanRefAtten = -5.08 * vanRho / 1.798;
87 const double vanScat = -5.1 * vanRho;
88
89 const size_t numHists = workspace->getNumberHistograms();
90 const size_t specSize = workspace->blocksize();
91
92 const auto &spectrumInfo = workspace->spectrumInfo();
93 Progress progress(this, 0.91, 1.0, numHists);
94 for (size_t i = 0; i < numHists; ++i) {
95
96 if (!spectrumInfo.hasDetectors(i)) {
97 // If a spectrum doesn't have an attached detector go to next one instead
98 continue;
99 }
100
101 // Get detector position
102 V3D detectorPos = spectrumInfo.position(i);
103
104 const int detID = spectrumInfo.detector(i).getID();
105 double angleFactor;
106 // If the low angle or backscattering bank, want angle wrt beamline
107 if (detID < 900000) {
108 double theta = detectorPos.angle(V3D(0.0, 0.0, 1.0));
109 angleFactor = 1.0 / std::abs(cos(theta));
110 }
111 // For 90 degree bank need it wrt X axis
112 else {
113 double theta = detectorPos.angle(V3D(1.0, 0.0, 0.0));
114 angleFactor = 1.0 / std::abs(cos(theta));
115 }
116
117 auto &Y = workspace->mutableY(i);
118 const auto lambdas = workspace->points(i);
119 for (size_t j = 0; j < specSize; ++j) {
120 const double lambda = lambdas[j];
121
122 // Front vanadium window - 0.5-1% effect, increasing with lambda
123 Y[j] *= exp(vanWinThickness * ((vanRefAtten * lambda) + vanScat));
124
125 // 90 degree bank aluminium sample holder
126 if (detID > 900000) {
127 // Below values are for aluminium
128 Y[j] *= exp(-angleFactor * 0.011 * 6.02 * ((0.231 * lambda / 1.798) + 1.503));
129 } else // Another vanadium window
130 {
131 Y[j] *= exp(angleFactor * vanWinThickness * ((vanRefAtten * lambda) + vanScat));
132 }
133 }
134
135 progress.report();
136 }
137
138 setProperty("OutputWorkspace", workspace);
139}
140
142 MatrixWorkspace_sptr m_inputWS = getProperty("InputWorkspace");
143 double sigma_atten = getProperty("SampleAttenuationXSection"); // in barns
144 double sigma_s = getProperty("SampleScatteringXSection"); // in barns
145 double rho = getProperty("SampleNumberDensity"); // in Angstroms-3
146 const Material &sampleMaterial = m_inputWS->sample().getMaterial();
147 if (sampleMaterial.totalScatterXSection() != 0.0) {
148 if (rho == EMPTY_DBL())
149 rho = sampleMaterial.numberDensity();
150 if (sigma_s == EMPTY_DBL())
151 sigma_s = sampleMaterial.totalScatterXSection();
152 if (sigma_atten == EMPTY_DBL())
153 sigma_atten = sampleMaterial.absorbXSection(NeutronAtom::ReferenceLambda);
154 } else // Save input in Sample with wrong atomic number and name
155 {
156 NeutronAtom neutron(0, 0, 0.0, 0.0, sigma_s, 0.0, sigma_s, sigma_atten);
157 auto shape = std::shared_ptr<IObject>(
158 m_inputWS->sample().getShape().cloneWithMaterial(Material("SetInSphericalAbsorption", neutron, rho)));
159 m_inputWS->mutableSample().setShape(shape);
160 }
161
162 // Call FlatPlateAbsorption as a Child Algorithm
163 auto childAlg = createChildAlgorithm("FlatPlateAbsorption", 0.0, 0.9);
164 // Pass through all the properties
165 childAlg->setProperty<MatrixWorkspace_sptr>("InputWorkspace", m_inputWS);
166 childAlg->setProperty<double>("AttenuationXSection", sigma_atten);
167 childAlg->setProperty<double>("ScatteringXSection", sigma_s);
168 childAlg->setProperty<double>("SampleNumberDensity", rho);
169 childAlg->setProperty<int64_t>("NumberOfWavelengthPoints", getProperty("NumberOfWavelengthPoints"));
170 childAlg->setProperty<std::string>("ExpMethod", getProperty("ExpMethod"));
171 childAlg->setProperty<double>("ElementSize", getProperty("ElementSize"));
172 // The height and width of the sample holder are standard for HRPD
173 const double HRPDCanHeight = 2.3;
174 const double HRPDCanWidth = 1.8;
175 childAlg->setProperty("SampleHeight", HRPDCanHeight);
176 childAlg->setProperty("SampleWidth", HRPDCanWidth);
177 const double thickness = getProperty("Thickness");
178 childAlg->setProperty("SampleThickness", thickness);
179 childAlg->executeAsChildAlg();
180 return childAlg->getProperty("OutputWorkspace");
181}
182
183} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
const std::vector< double > * lambda
IPeaksWorkspace_sptr workspace
Definition: IndexPeaks.cpp:114
void declareProperty(std::unique_ptr< Kernel::Property > p, const std::string &doc="") override
Add a property to the list of managed properties.
Definition: Algorithm.cpp:1913
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
Definition: Algorithm.cpp:2076
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.
Definition: Algorithm.cpp:842
void progress(double p, const std::string &msg="", double estimatedTime=0.0, int progressPrecision=0)
Sends ProgressNotification.
Definition: Algorithm.cpp:231
Helper class for reporting progress from algorithms.
Definition: Progress.h:25
A property class for workspaces.
void init() override
Initialisation code.
API::MatrixWorkspace_sptr runFlatPlateAbsorption()
runs the flat plate absorption
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
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
double totalScatterXSection() const
Return the total scattering cross section for a given wavelength in barns.
Definition: Material.cpp:252
Class for 3D vectors.
Definition: V3D.h:34
double angle(const V3D &) const
Angle between this and another vector.
Definition: V3D.cpp:165
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
A namespace containing physical constants that are required by algorithms and unit routines.
Definition: Atom.h:14
constexpr int EMPTY_INT() noexcept
Returns what we consider an "empty" integer within a property.
Definition: EmptyValues.h:25
constexpr double EMPTY_DBL() noexcept
Returns what we consider an "empty" double within a property.
Definition: EmptyValues.h:43
@ 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
static const double ReferenceLambda
The reference wavelength value for absorption cross sections.
Definition: NeutronAtom.h:25