-
Notifications
You must be signed in to change notification settings - Fork 122
/
CylinderAbsorption.cpp
139 lines (117 loc) · 5.16 KB
/
CylinderAbsorption.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
#include "MantidAlgorithms/CylinderAbsorption.h"
#include "MantidAPI/MatrixWorkspace.h"
#include "MantidGeometry/Instrument.h"
#include "MantidGeometry/Objects/Object.h"
#include "MantidGeometry/Objects/Track.h"
#include "MantidKernel/BoundedValidator.h"
namespace Mantid {
namespace Algorithms {
// Register the algorithm into the AlgorithmFactory
DECLARE_ALGORITHM(CylinderAbsorption)
using namespace Kernel;
using namespace Geometry;
using namespace API;
CylinderAbsorption::CylinderAbsorption()
: AbsorptionCorrection(), m_cylHeight(0.0), m_cylRadius(0.0),
m_numSlices(0), m_sliceThickness(0), m_numAnnuli(0), m_deltaR(0) {}
void CylinderAbsorption::defineProperties() {
auto mustBePositive = boost::make_shared<BoundedValidator<double>>();
mustBePositive->setLower(0.0);
declareProperty("CylinderSampleHeight", -1.0, mustBePositive,
"The height of the cylindrical sample in centimetres");
declareProperty("CylinderSampleRadius", -1.0, mustBePositive,
"The radius of the cylindrical sample in centimetres");
auto positiveInt = boost::make_shared<BoundedValidator<int>>();
positiveInt->setLower(1);
declareProperty(
"NumberOfSlices", 1, positiveInt,
"The number of slices into which the cylinder is divided for the\n"
"calculation");
declareProperty(
"NumberOfAnnuli", 1, positiveInt,
"The number of annuli into which each slice is divided for the\n"
"calculation");
}
/// Fetch the properties and set the appropriate member variables
void CylinderAbsorption::retrieveProperties() {
m_cylHeight = getProperty("CylinderSampleHeight"); // in cm
m_cylRadius = getProperty("CylinderSampleRadius"); // in cm
m_cylHeight *= 0.01; // now in m
m_cylRadius *= 0.01; // now in m
m_numSlices = getProperty("NumberOfSlices");
m_sliceThickness = m_cylHeight / m_numSlices;
m_numAnnuli = getProperty("NumberOfAnnuli");
m_deltaR = m_cylRadius / m_numAnnuli;
/* The number of volume elements is
* numslices*(1+2+3+.....+numAnnuli)*6
* Since the first annulus is separated in 6 segments, the next one in 12 and
* so on.....
*/
m_numVolumeElements = m_numSlices * m_numAnnuli * (m_numAnnuli + 1) * 3;
if (m_numVolumeElements == 0) {
g_log.error() << "Input properties lead to no defined volume elements.\n";
throw std::runtime_error("No volume elements defined.");
}
m_sampleVolume = m_cylHeight * M_PI * m_cylRadius * m_cylRadius;
if (m_sampleVolume == 0.0) {
g_log.error() << "Defined sample has zero volume.\n";
throw std::runtime_error("Sample with zero volume defined.");
}
}
std::string CylinderAbsorption::sampleXML() {
// Get the sample position, which is typically the origin but we should be
// generic
const V3D samplePos = m_inputWS->getInstrument()->getSample()->getPos();
// Shift so that cylinder is centered at sample position
const double cylinderBase = (-0.5 * m_cylHeight) + samplePos.Y();
std::ostringstream xmlShapeStream;
xmlShapeStream << "<cylinder id=\"detector-shape\"> "
<< "<centre-of-bottom-base x=\"" << samplePos.X() << "\" y=\""
<< cylinderBase << "\" z=\"" << samplePos.Z() << "\" /> "
<< "<axis x=\"0\" y=\"1\" z=\"0\" /> "
<< "<radius val=\"" << m_cylRadius << "\" /> "
<< "<height val=\"" << m_cylHeight << "\" /> "
<< "</cylinder>";
return xmlShapeStream.str();
}
/// Calculate the distances for L1 and element size for each element in the
/// sample
void CylinderAbsorption::initialiseCachedDistances() {
m_L1s.resize(m_numVolumeElements);
m_elementVolumes.resize(m_numVolumeElements);
m_elementPositions.resize(m_numVolumeElements);
int counter = 0;
// loop over slices
for (int i = 0; i < m_numSlices; ++i) {
const double z = (i + 0.5) * m_sliceThickness - 0.5 * m_cylHeight;
// Number of elements in 1st annulus
int Ni = 0;
// loop over annuli
for (int j = 0; j < m_numAnnuli; ++j) {
Ni += 6;
const double R = (j * m_cylRadius / m_numAnnuli) + (m_deltaR / 2.0);
// loop over elements in current annulus
for (int k = 0; k < Ni; ++k) {
const double phi = 2 * M_PI * k / Ni;
// Calculate the current position in the sample in Cartesian
// coordinates.
// Remember that our cylinder has its axis along the y axis
m_elementPositions[counter](R * sin(phi), z, R * cos(phi));
assert(m_sampleObject->isValid(m_elementPositions[counter]));
// Create track for distance in cylinder before scattering point
Track incoming(m_elementPositions[counter], m_beamDirection * -1.0);
m_sampleObject->interceptSurface(incoming);
m_L1s[counter] = incoming.cbegin()->distFromStart;
// Also calculate element volumes here
const double outerR = R + (m_deltaR / 2.0);
const double innerR = outerR - m_deltaR;
const double elementVolume =
M_PI * (outerR * outerR - innerR * innerR) * m_sliceThickness / Ni;
m_elementVolumes[counter] = elementVolume;
counter++;
}
}
}
}
} // namespace Algorithms
} // namespace Mantid