-
Notifications
You must be signed in to change notification settings - Fork 122
/
ApplyTransmissionCorrection.cpp
188 lines (160 loc) · 8.08 KB
/
ApplyTransmissionCorrection.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
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
176
177
178
179
180
181
182
183
184
185
186
187
188
// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
// NScD Oak Ridge National Laboratory, European Spallation Source,
// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
// SPDX - License - Identifier: GPL - 3.0 +
#include "MantidAlgorithms/ApplyTransmissionCorrection.h"
#include "MantidAPI/HistogramValidator.h"
#include "MantidAPI/SpectrumInfo.h"
#include "MantidAPI/WorkspaceOpOverloads.h"
#include "MantidAPI/WorkspaceUnitValidator.h"
#include "MantidDataObjects/Workspace2D.h"
#include "MantidDataObjects/WorkspaceCreation.h"
#include "MantidGeometry/IDetector.h"
#include "MantidKernel/BoundedValidator.h"
#include "MantidKernel/CompositeValidator.h"
namespace Mantid {
namespace Algorithms {
// Register the algorithm into the AlgorithmFactory
DECLARE_ALGORITHM(ApplyTransmissionCorrection)
using namespace Kernel;
using namespace API;
using namespace Geometry;
using namespace HistogramData;
using namespace DataObjects;
namespace {
namespace PropertyNames {
const std::string INPUT_WKSP("InputWorkspace");
const std::string OUTPUT_WKSP("OutputWorkspace");
const std::string TRANSMISSION_WKSP("TransmissionWorkspace");
const std::string TRANSMISSION_VALUE("TransmissionValue");
const std::string TRANSMISSION_ERROR("TransmissionError");
} // namespace PropertyNames
} // namespace
void ApplyTransmissionCorrection::init() {
auto wsValidator = std::make_shared<CompositeValidator>();
wsValidator->add<WorkspaceUnitValidator>("Wavelength");
wsValidator->add<HistogramValidator>();
declareProperty(std::make_unique<WorkspaceProperty<>>(PropertyNames::INPUT_WKSP, "", Direction::Input, wsValidator),
"Workspace to apply the transmission correction to");
declareProperty(std::make_unique<WorkspaceProperty<>>(PropertyNames::TRANSMISSION_WKSP, "", Direction::Output,
PropertyMode::Optional),
"Workspace containing the transmission values [optional]");
declareProperty(std::make_unique<WorkspaceProperty<>>(PropertyNames::OUTPUT_WKSP, "", Direction::Output),
"Workspace to store the corrected data in");
// Alternatively, the user can specify a transmission that will ba applied to
// all wavelength bins
declareProperty(PropertyNames::TRANSMISSION_VALUE, EMPTY_DBL(),
"Transmission value to apply to all wavelengths. If specified, "
"TransmissionWorkspace will not be used.");
auto mustBePositive = std::make_shared<BoundedValidator<double>>();
mustBePositive->setLower(0.0);
declareProperty(PropertyNames::TRANSMISSION_ERROR, 0.0, mustBePositive,
"The error on the transmission value (default 0.0)");
declareProperty("ThetaDependent", true, "If true, a theta-dependent transmission correction will be applied.");
}
std::map<std::string, std::string> ApplyTransmissionCorrection::validateInputs() {
std::map<std::string, std::string> result;
// needed to find out what mode we are working in
const double trans_value = getProperty(PropertyNames::TRANSMISSION_VALUE);
if (isEmpty(trans_value)) { // require a transmission workspace
MatrixWorkspace_const_sptr transWS = getProperty(PropertyNames::TRANSMISSION_WKSP);
if (!transWS) {
const std::string msg("Must specify \"TransmissionValue\" or \"TransmissionWorkspace\"");
result[PropertyNames::TRANSMISSION_VALUE] = msg;
result[PropertyNames::TRANSMISSION_WKSP] = msg;
} else {
MatrixWorkspace_const_sptr inputWS = getProperty(PropertyNames::INPUT_WKSP);
if ((transWS->getNumberHistograms() > 1) && (transWS->getNumberHistograms() != inputWS->getNumberHistograms())) {
const std::string msg("Input and transmission workspaces have "
"incompatible number of spectra");
result[PropertyNames::INPUT_WKSP] = msg;
result[PropertyNames::TRANSMISSION_WKSP] = msg;
} else if (transWS->y(0).size() != inputWS->y(0).size()) {
const std::string msg("Input and transmission workspaces have a "
"different number of wavelength bins");
result[PropertyNames::INPUT_WKSP] = msg;
result[PropertyNames::TRANSMISSION_WKSP] = msg;
}
}
}
return result;
}
void ApplyTransmissionCorrection::exec() {
const bool thetaDependent = getProperty("ThetaDependent");
const double trans_value = getProperty(PropertyNames::TRANSMISSION_VALUE);
const double trans_error = getProperty(PropertyNames::TRANSMISSION_ERROR);
MatrixWorkspace_const_sptr transWS = getProperty(PropertyNames::TRANSMISSION_WKSP);
// Check whether we only need to divided the workspace by the transmission.
// theta-dependence modifies the input transmission information
if (thetaDependent) {
MatrixWorkspace_sptr inputWS = getProperty(PropertyNames::INPUT_WKSP);
// initial values of tranmission come from the value
HistogramY TrIn(inputWS->y(0).size(), trans_value);
HistogramE ETrIn(inputWS->y(0).size(), trans_error);
if (isEmpty(trans_value)) {
// override the values specified above
TrIn = transWS->y(0);
ETrIn = transWS->e(0);
}
const auto numHists = static_cast<int>(inputWS->getNumberHistograms());
Progress progress(this, 0.0, 1.0, numHists);
// Create a Workspace2D to hold the correction (applied through
// multiplication)
MatrixWorkspace_sptr corrWS = create<HistoWorkspace>(*inputWS);
const auto &spectrumInfo = inputWS->spectrumInfo();
// Loop through the spectra and apply correction
PARALLEL_FOR_IF(Kernel::threadSafe(*inputWS, *corrWS))
for (int i = 0; i < numHists; i++) {
PARALLEL_START_INTERUPT_REGION
if (!spectrumInfo.hasDetectors(i)) {
g_log.warning() << "Workspace index " << i << " has no detector assigned to it - discarding'\n";
continue;
}
// Copy over the X data
corrWS->setSharedX(i, inputWS->sharedX(i));
// Skip if we have a monitor or if the detector is masked.
if (spectrumInfo.isMonitor(i) || spectrumInfo.isMasked(i))
continue;
// Compute theta-dependent transmission term for each wavelength bin
auto &correctionY = corrWS->mutableY(i);
auto &correctionE = corrWS->mutableE(i);
const double exp_term = 0.5 / cos(spectrumInfo.twoTheta(i)) + 0.5;
for (int j = 0; j < static_cast<int>(inputWS->y(0).size()); j++) {
correctionY[j] = pow(TrIn[j], -1. * exp_term); // 1 / TrIn^exp_term
correctionE[j] = std::fabs(ETrIn[j] * exp_term / pow(TrIn[j], exp_term + 1.0));
}
progress.report("Calculating transmission correction");
PARALLEL_END_INTERUPT_REGION
}
PARALLEL_CHECK_INTERUPT_REGION
// apply the correction
MatrixWorkspace_sptr outputWS = inputWS * corrWS;
;
setProperty(PropertyNames::OUTPUT_WKSP, outputWS);
} else { // !thetaDependent case - divide does the actual work
// this case divide and set output workspace
auto divideAlg = createChildAlgorithm("Divide", 0.0, 1.0);
divideAlg->setProperty("LHSWorkspace", getPropertyValue(PropertyNames::INPUT_WKSP));
if (transWS) { // use the supplied transmission workspace
divideAlg->setProperty("RHSWorkspace", getPropertyValue(PropertyNames::TRANSMISSION_WKSP));
} else {
// set the RHS to be a single value workspace with the requested
// uncertainty
auto createSingleAlg = createChildAlgorithm("CreateSingleValuedWorkspace");
createSingleAlg->setProperty("DataValue", trans_value);
createSingleAlg->setProperty("ErrorValue", trans_error);
createSingleAlg->executeAsChildAlg();
MatrixWorkspace_sptr singleWS = createSingleAlg->getProperty("OutputWorkspace");
divideAlg->setProperty("RHSWorkspace", singleWS);
}
divideAlg->setProperty("OutputWorkspace", getPropertyValue(PropertyNames::OUTPUT_WKSP));
divideAlg->executeAsChildAlg();
// call divide and set output workspace
MatrixWorkspace_sptr outputWS = divideAlg->getProperty("OutputWorkspace");
setProperty(PropertyNames::OUTPUT_WKSP, outputWS);
}
}
} // namespace Algorithms
} // namespace Mantid