Newer
Older
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
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
#include "MantidAlgorithms/CarpenterSampleCorrection.h"
#include "MantidAPI/AlgorithmManager.h"
#include "MantidAPI/InstrumentValidator.h"
#include "MantidAPI/Sample.h"
#include "MantidAPI/SpectrumInfo.h"
#include "MantidAPI/WorkspaceFactory.h"
#include "MantidAPI/WorkspaceUnitValidator.h"
#include "MantidDataObjects/EventWorkspace.h"
#include "MantidGeometry/Instrument.h"
#include "MantidGeometry/IComponent.h"
#include "MantidKernel/CompositeValidator.h"
#include "MantidKernel/Material.h"
#include "MantidKernel/PhysicalConstants.h"
#include <stdexcept>
namespace Mantid {
namespace Algorithms {
DECLARE_ALGORITHM(CarpenterSampleCorrection) // Register the class
// into the algorithm
// factory
using namespace Kernel;
using namespace API;
using Mantid::DataObjects::EventList;
using Mantid::DataObjects::EventWorkspace;
using Mantid::DataObjects::EventWorkspace_sptr;
using Mantid::DataObjects::WeightedEventNoTime;
using Mantid::HistogramData::HistogramX;
using Mantid::HistogramData::HistogramY;
using Mantid::HistogramData::HistogramE;
using std::vector;
using namespace Mantid::PhysicalConstants;
using namespace Geometry;
// Constants required internally only, so make them static. These are
// Chebyshev expansion coefficients copied directly from Carpenter 1969 Table 1
// clang-format off
const std::string CarpenterSampleCorrection::name() const {
return "CarpenterSampleCorrection";
}
int CarpenterSampleCorrection::version() const { return 1; }
const std::string CarpenterSampleCorrection::category() const {
return "CorrectionFunctions\\AbsorptionCorrections";
}
/**
* Initialize the properties to default values
*/
void CarpenterSampleCorrection::init() {
// The input workspace must have an instrument and units of wavelength
auto wsValidator = boost::make_shared<CompositeValidator>();
wsValidator->add<WorkspaceUnitValidator>("Wavelength");
wsValidator->add<InstrumentValidator>();
declareProperty(make_unique<WorkspaceProperty<API::MatrixWorkspace> >(
"InputWorkspace", "", Direction::Input, wsValidator),
"The name of the input workspace.");
declareProperty(make_unique<WorkspaceProperty<API::MatrixWorkspace> >(
"OutputWorkspace", "", Direction::Output),
"The name of the output workspace.");
declareProperty("AttenuationXSection", 2.8, "Coefficient 1, absorption cross "
"section / 1.81 if not set with "
"SetSampleMaterial");
declareProperty("ScatteringXSection", 5.1, "Coefficient 3, total scattering "
"cross section if not set with "
"SetSampleMaterial");
declareProperty("SampleNumberDensity", 0.0721,
"Coefficient 2, density if not set with SetSampleMaterial");
declareProperty("CylinderSampleRadius", 0.3175, "Sample radius, in cm");
}
/**
* Execute the algorithm
*/
void CarpenterSampleCorrection::exec() {
// common information
MatrixWorkspace_sptr inputWksp = getProperty("InputWorkspace");
double radius = getProperty("CylinderSampleRadius");
double coeff1 = getProperty("AttenuationXSection");
double coeff2 = getProperty("SampleNumberDensity");
double coeff3 = getProperty("ScatteringXSection");
// Calculate the absorption and multiple scattering corrections
WorkspaceGroup_sptr calcOutput = calculateCorrection(
inputWksp, radius, coeff1, coeff2, coeff3, true, true);
Workspace_sptr absPtr = calcOutput->getItem(0);
Workspace_sptr msPtr = calcOutput->getItem(1);
auto absWksp = boost::dynamic_pointer_cast<MatrixWorkspace>(absPtr);
auto msWksp = boost::dynamic_pointer_cast<MatrixWorkspace>(msPtr);
// Clone input -> to apply corrections
MatrixWorkspace_sptr outputWksp = getProperty("OutputWorkspace");
EventWorkspace_sptr inputWkspEvent =
boost::dynamic_pointer_cast<EventWorkspace>(inputWksp);
outputWksp = inputWksp->clone();
// Apply the absorption correction to the sample workspace
outputWksp = divide(outputWksp, absWksp);
// Apply the multiple scattering correction to the sample workspace
auto factorWksp = multiply(inputWksp, msWksp);
outputWksp = minus(outputWksp, factorWksp);
// Output workspace
if (inputWkspEvent) {
auto outputWkspEvent =
boost::dynamic_pointer_cast<EventWorkspace>(outputWksp);
setProperty("OutputWorkspace", outputWkspEvent);
}
setProperty("OutputWorkspace", outputWksp);
}
WorkspaceGroup_sptr CarpenterSampleCorrection::calculateCorrection(
MatrixWorkspace_sptr &inputWksp, double radius, double coeff1,
double coeff2, double coeff3, bool doAbs, bool doMS) {
auto calculate =
this->createChildAlgorithm("CalculateCarpenterSampleCorrection");
calculate->initialize();
calculate->setProperty("InputWorkspace", inputWksp);
calculate->setProperty("CylinderSampleRadius", radius);
calculate->setProperty("AttenuationXSection", coeff1);
calculate->setProperty("SampleNumberDensity", coeff2);
calculate->setProperty("ScatteringXSection", coeff3);
calculate->setProperty("Absorption", doAbs);
calculate->setProperty("MultipleScattering", doMS);
calculate->execute();
WorkspaceGroup_sptr calcOutput =
calculate->getProperty("OutputWorkspaceBaseName");
return calcOutput;
}
MatrixWorkspace_sptr
CarpenterSampleCorrection::divide(const MatrixWorkspace_sptr lhsWS,
const MatrixWorkspace_sptr rhsWS) {
IAlgorithm_sptr divide = this->createChildAlgorithm("Divide");
divide->initialize();
divide->setProperty("LHSWorkspace", lhsWS);
divide->setProperty("RHSWorkspace", rhsWS);
divide->execute();
MatrixWorkspace_sptr outWS = divide->getProperty("OutputWorkspace");
return outWS;
}
MatrixWorkspace_sptr
CarpenterSampleCorrection::multiply(const MatrixWorkspace_sptr lhsWS,
const MatrixWorkspace_sptr rhsWS) {
auto multiply = this->createChildAlgorithm("Multiply");
multiply->initialize();
multiply->setProperty("LHSWorkspace", lhsWS);
multiply->setProperty("RHSWorkspace", rhsWS);
multiply->execute();
MatrixWorkspace_sptr outWS = multiply->getProperty("OutputWorkspace");
return outWS;
}
MatrixWorkspace_sptr
CarpenterSampleCorrection::minus(const MatrixWorkspace_sptr lhsWS,
const MatrixWorkspace_sptr rhsWS) {
auto minus = this->createChildAlgorithm("Minus");
minus->initialize();
minus->setProperty("LHSWorkspace", lhsWS);
minus->setProperty("RHSWorkspace", rhsWS);
minus->execute();
MatrixWorkspace_sptr outWS = minus->getProperty("OutputWorkspace");
return outWS;
}
} // namespace Algorithm
} // namespace Mantid