-
Notifications
You must be signed in to change notification settings - Fork 122
/
DetectorEfficiencyVariation.cpp
263 lines (236 loc) · 10.1 KB
/
DetectorEfficiencyVariation.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
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
#include "MantidAlgorithms/DetectorEfficiencyVariation.h"
#include "MantidAPI/HistogramValidator.h"
#include "MantidAPI/SpectrumInfo.h"
#include "MantidKernel/BoundedValidator.h"
namespace Mantid {
namespace Algorithms {
// Register the class into the algorithm factory
DECLARE_ALGORITHM(DetectorEfficiencyVariation)
const std::string DetectorEfficiencyVariation::category() const {
return "Diagnostics";
}
using namespace Kernel;
using namespace API;
using DataObjects::MaskWorkspace_sptr;
using Geometry::IDetector_const_sptr;
/// Initialize the algorithm
void DetectorEfficiencyVariation::init() {
auto val = boost::make_shared<HistogramValidator>();
declareProperty(make_unique<WorkspaceProperty<MatrixWorkspace>>(
"WhiteBeamBase", "", Direction::Input, val),
"Name of a white beam vanadium workspace");
// The histograms, the detectors in each histogram and their first and last
// bin boundary must match
declareProperty(
make_unique<WorkspaceProperty<MatrixWorkspace>>("WhiteBeamCompare", "",
Direction::Input, val),
"Name of a matching second white beam vanadium run from the same "
"instrument");
declareProperty(make_unique<WorkspaceProperty<MatrixWorkspace>>(
"OutputWorkspace", "", Direction::Output),
"A MaskWorkpace where each spectra that failed the test is "
"masked. Each histogram from the input workspace maps to a "
"histogram in this workspace with one value that indicates "
"if there was a dead detector.");
auto moreThanZero = boost::make_shared<BoundedValidator<double>>();
moreThanZero->setLower(0.0);
declareProperty("Variation", 1.1, moreThanZero,
"Identify histograms whose total number of counts has "
"changed by more than this factor of the median change "
"between the two input workspaces.");
auto mustBePosInt = boost::make_shared<BoundedValidator<int>>();
mustBePosInt->setLower(0);
declareProperty("StartWorkspaceIndex", 0, mustBePosInt,
"The index number of the first spectrum to include in the "
"calculation (default: 0)");
// Mantid::EMPTY_INT() and EMPTY_DBL() are tags that indicate that no
// value has been set and we want to use the default
declareProperty("EndWorkspaceIndex", Mantid::EMPTY_INT(), mustBePosInt,
"The index number of the last spectrum to include in the "
"calculation (default: the last spectrum in the workspace)");
declareProperty(
"RangeLower", Mantid::EMPTY_DBL(),
"No bin with a boundary at an x value less than this will be included "
"in the summation used to decide if a detector is 'bad' (default: the "
"start of each histogram)");
declareProperty(
"RangeUpper", Mantid::EMPTY_DBL(),
"No bin with a boundary at an x value higher than this value will "
"be included in the summation used to decide if a detector is 'bad' "
"(default: the end of each histogram)");
declareProperty("NumberOfFailures", 0, Direction::Output);
}
/** Executes the algorithm that includes calls to SolidAngle and Integration
*
* @throw invalid_argument if there is an incapatible property value and so the
*algorithm can't continue
* @throw runtime_error if a Child Algorithm cannot execute
*/
void DetectorEfficiencyVariation::exec() {
MatrixWorkspace_sptr WB1;
MatrixWorkspace_sptr WB2;
double variation = Mantid::EMPTY_DBL();
int minSpec = 0;
int maxSpec = Mantid::EMPTY_INT();
retrieveProperties(WB1, WB2, variation, minSpec, maxSpec);
const double rangeLower = getProperty("RangeLower");
const double rangeUpper = getProperty("RangeUpper");
MatrixWorkspace_sptr counts1 =
integrateSpectra(WB1, minSpec, maxSpec, rangeLower, rangeUpper);
MatrixWorkspace_sptr counts2 =
integrateSpectra(WB2, minSpec, maxSpec, rangeLower, rangeUpper);
MatrixWorkspace_sptr countRatio;
try {
// Note. This can produce NAN/INFs. Leave for now and sort it out in the
// later tests
countRatio = counts1 / counts2;
} catch (std::invalid_argument &) {
g_log.error() << "The two white beam workspaces size must match.";
throw;
}
double average =
calculateMedian(*countRatio, false, makeInstrumentMap(*countRatio))
.at(0); // Include zeroes
g_log.notice() << name()
<< ": The median of the ratio of the integrated counts is: "
<< average << '\n';
//
int numFailed = doDetectorTests(counts1, counts2, average, variation);
g_log.notice() << "Tests failed " << numFailed << " spectra. "
<< "These have been masked on the OutputWorkspace\n";
// counts1 was overwriten by the last function, now register it
setProperty("NumberOfFailures", numFailed);
}
/** Loads, checks and passes back the values passed to the algorithm
* @param whiteBeam1 :: A white beam vanadium spectrum that will be used to
* check detector efficiency variations
* @param whiteBeam2 :: The other white beam vanadium spectrum from the same
* instrument to use for comparison
* @param variation :: The maximum fractional variation above the median that is
* allowed for god detectors
* @param startWsIndex :: Index number of the first spectrum to use
* @param endWsIndex :: Index number of the last spectrum to use
* @throw invalid_argument if there is an incapatible property value and so the
* algorithm can't continue
*/
void DetectorEfficiencyVariation::retrieveProperties(
API::MatrixWorkspace_sptr &whiteBeam1,
API::MatrixWorkspace_sptr &whiteBeam2, double &variation, int &startWsIndex,
int &endWsIndex) {
whiteBeam1 = getProperty("WhiteBeamBase");
whiteBeam2 = getProperty("WhiteBeamCompare");
if (whiteBeam1->getInstrument()->getName() !=
whiteBeam2->getInstrument()->getName()) {
throw std::invalid_argument("The two input white beam vanadium workspaces "
"must be from the same instrument");
}
int maxWsIndex = static_cast<int>(whiteBeam1->getNumberHistograms()) - 1;
if (maxWsIndex !=
static_cast<int>(whiteBeam2->getNumberHistograms()) -
1) { // we would get a crash later on if this were not true
throw std::invalid_argument("The input white beam vanadium workspaces must "
"be have the same number of histograms");
}
variation = getProperty("Variation");
startWsIndex = getProperty("StartWorkspaceIndex");
if ((startWsIndex < 0) || (startWsIndex > maxWsIndex)) {
g_log.warning("StartWorkspaceIndex out of range, changed to 0");
startWsIndex = 0;
}
endWsIndex = getProperty("EndWorkspaceIndex");
if (endWsIndex == Mantid::EMPTY_INT())
endWsIndex = maxWsIndex;
if ((endWsIndex < 0) || (endWsIndex > maxWsIndex)) {
g_log.warning(
"EndWorkspaceIndex out of range, changed to max Workspace number");
endWsIndex = maxWsIndex;
}
if ((endWsIndex < startWsIndex)) {
g_log.warning(
"EndWorkspaceIndex can not be less than the StartWorkspaceIndex, "
"changed to max Workspace number");
endWsIndex = maxWsIndex;
}
}
/**
* Apply the detector test criterion
* @param counts1 :: A workspace containing the integrated counts of the first
* white beam run
* @param counts2 :: A workspace containing the integrated counts of the first
* white beam run
* @param average :: The computed median
* @param variation :: The allowed variation in terms of number of medians, i.e
* those spectra where
* the ratio of the counts outside this range will fail the tests and will be
* masked on counts1
* @return number of detectors for which tests failed
*/
int DetectorEfficiencyVariation::doDetectorTests(
API::MatrixWorkspace_const_sptr counts1,
API::MatrixWorkspace_const_sptr counts2, const double average,
double variation) {
// DIAG in libISIS did this. A variation of less than 1 doesn't make sense in
// this algorithm
if (variation < 1) {
variation = 1.0 / variation;
}
// criterion for if the the first spectrum is larger than expected
double largest = average * variation;
// criterion for if the the first spectrum is lower than expected
double lowest = average / variation;
const int numSpec = static_cast<int>(counts1->getNumberHistograms());
const int progStep = static_cast<int>(std::ceil(numSpec / 30.0));
// Create a workspace for the output
MaskWorkspace_sptr maskWS = this->generateEmptyMask(counts1);
bool checkForMask = false;
Geometry::Instrument_const_sptr instrument = counts1->getInstrument();
if (instrument != nullptr) {
checkForMask = ((instrument->getSource() != nullptr) &&
(instrument->getSample() != nullptr));
}
const double deadValue(1.0);
int numFailed(0);
const auto &spectrumInfo = counts1->spectrumInfo();
PARALLEL_FOR_IF(Kernel::threadSafe(*counts1, *counts2, *maskWS))
for (int i = 0; i < numSpec; ++i) {
PARALLEL_START_INTERUPT_REGION
// move progress bar
if (i % progStep == 0) {
advanceProgress(progStep * static_cast<double>(RTMarkDetects) / numSpec);
progress(m_fracDone);
interruption_point();
}
if (checkForMask) {
if (spectrumInfo.isMonitor(i))
continue;
if (spectrumInfo.isMasked(i)) {
// Ensure it is masked on the output
maskWS->mutableY(i)[0] = deadValue;
continue;
}
}
const double signal1 = counts1->y(i)[0];
const double signal2 = counts2->y(i)[0];
// Mask out NaN and infinite
if (!std::isfinite(signal1) || !std::isfinite(signal2)) {
maskWS->mutableY(i)[0] = deadValue;
PARALLEL_ATOMIC
++numFailed;
continue;
}
// Check the ratio is within the given range
const double ratio = signal1 / signal2;
if (ratio < lowest || ratio > largest) {
maskWS->mutableY(i)[0] = deadValue;
PARALLEL_ATOMIC
++numFailed;
}
PARALLEL_END_INTERUPT_REGION
}
PARALLEL_CHECK_INTERUPT_REGION
// Register the results with the ADS
setProperty("OutputWorkspace", maskWS);
return numFailed;
}
} // namespace Algorithm
} // namespace Mantid