-
Notifications
You must be signed in to change notification settings - Fork 122
/
CalculateEfficiency.cpp
370 lines (327 loc) · 14.8 KB
/
CalculateEfficiency.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
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
// 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/CalculateEfficiency.h"
#include "MantidAPI/SpectrumInfo.h"
#include "MantidAPI/WorkspaceFactory.h"
#include "MantidDataObjects/EventList.h"
#include "MantidDataObjects/EventWorkspace.h"
#include "MantidGeometry/ICompAssembly.h"
#include "MantidGeometry/IDTypes.h"
#include "MantidGeometry/IDetector.h"
#include "MantidGeometry/Instrument.h"
#include "MantidGeometry/Instrument/Detector.h"
#include "MantidGeometry/Instrument/RectangularDetector.h"
#include "MantidKernel/ArrayProperty.h"
#include "MantidKernel/BoundedValidator.h"
#include <vector>
namespace Mantid::Algorithms {
// Register the class into the algorithm factory
DECLARE_ALGORITHM(CalculateEfficiency)
using namespace Kernel;
using namespace API;
using namespace Geometry;
using namespace DataObjects;
/** Initialization method.
*
*/
void CalculateEfficiency::init() {
declareProperty(std::make_unique<WorkspaceProperty<>>("InputWorkspace", "", Direction::Input),
"The workspace containing the flood data");
declareProperty(std::make_unique<WorkspaceProperty<>>("OutputWorkspace", "", Direction::Output),
"The name of the workspace to be created as the output of the algorithm");
auto positiveDouble = std::make_shared<BoundedValidator<double>>();
positiveDouble->setLower(0);
declareProperty("MinEfficiency", EMPTY_DBL(), positiveDouble,
"Minimum efficiency for a pixel to be considered (default: no minimum).");
declareProperty("MaxEfficiency", EMPTY_DBL(), positiveDouble,
"Maximum efficiency for a pixel to be considered (default: no maximum).");
declareProperty("MaskedFullComponent", "", "Component Name to fully mask according to the IDF file.");
declareProperty(std::make_unique<ArrayProperty<int>>("MaskedEdges"),
"Number of pixels to mask on the edges: X-low, X-high, Y-low, Y-high");
declareProperty("MaskedComponent", "", "Component Name to mask the edges according to the IDF file.");
}
/** Executes the algorithm
*
*/
void CalculateEfficiency::exec() {
// Minimum efficiency. Pixels with lower efficiency will be masked
double min_eff = getProperty("MinEfficiency");
// Maximum efficiency. Pixels with higher efficiency will be masked
double max_eff = getProperty("MaxEfficiency");
// Get the input workspace
MatrixWorkspace_sptr inputWS = getProperty("InputWorkspace");
MatrixWorkspace_sptr rebinnedWS; // = inputWS;
// BioSANS has 2 detectors and reduces one at the time: one is masked!
// We must use that masked detector
const std::string maskedFullComponent = getPropertyValue("MaskedFullComponent");
if (!maskedFullComponent.empty()) {
g_log.debug() << "CalculateEfficiency: Masking Full Component: " << maskedFullComponent << "\n";
maskComponent(*inputWS, maskedFullComponent);
}
// BioSANS has 2 detectors and the front masks the back!!!!
// We must mask the shaded part to calculate efficency
std::vector<int> maskedEdges = getProperty("MaskedEdges");
if (!maskedEdges.empty() && (maskedEdges[0] > 0 || maskedEdges[1] > 0 || maskedEdges[2] > 0 || maskedEdges[3] > 0)) {
g_log.debug() << "CalculateEfficiency: Masking edges length = " << maskedEdges.size() << ")"
<< " of the component " << maskedFullComponent << "\n";
const std::string maskedComponent = getPropertyValue("MaskedComponent");
maskEdges(inputWS, maskedEdges[0], maskedEdges[1], maskedEdges[2], maskedEdges[3], maskedComponent);
}
// Now create the output workspace
MatrixWorkspace_sptr outputWS; // = getProperty("OutputWorkspace");
// DataObjects::EventWorkspace_const_sptr inputEventWS =
// std::dynamic_pointer_cast<const EventWorkspace>(inputWS);
// Sum up all the wavelength bins
auto childAlg = createChildAlgorithm("Integration", 0.0, 0.2);
childAlg->setProperty<MatrixWorkspace_sptr>("InputWorkspace", inputWS);
childAlg->executeAsChildAlg();
rebinnedWS = childAlg->getProperty("OutputWorkspace");
outputWS = WorkspaceFactory::Instance().create(rebinnedWS);
WorkspaceFactory::Instance().initializeFromParent(*inputWS, *outputWS, false);
for (int i = 0; i < static_cast<int>(rebinnedWS->getNumberHistograms()); i++) {
outputWS->setSharedX(i, rebinnedWS->sharedX(i));
}
setProperty("OutputWorkspace", outputWS);
double sum = 0.0;
double err = 0.0;
int npixels = 0;
// Loop over spectra and sum all the counts to get normalization
// Skip monitors and masked detectors
sumUnmaskedDetectors(rebinnedWS, sum, err, npixels);
// Normalize each detector pixel by the sum we just found to get the
// relative efficiency. If the minimum and maximum efficiencies are
// provided, the pixels falling outside this range will be marked
// as 'masked' in both the input and output workspace.
// We mask detectors in the input workspace so that we can resum the
// counts to find a new normalization factor that takes into account
// the newly masked detectors.
normalizeDetectors(rebinnedWS, outputWS, sum, err, npixels, min_eff, max_eff);
if (!isEmpty(min_eff) || !isEmpty(max_eff)) {
// Recompute the normalization, excluding the pixels that were outside
// the acceptable efficiency range.
sumUnmaskedDetectors(rebinnedWS, sum, err, npixels);
// Now that we have a normalization factor that excludes bad pixels,
// recompute the relative efficiency.
// We pass EMPTY_DBL() to avoid masking pixels that might end up high or low
// after the new normalization.
normalizeDetectors(rebinnedWS, outputWS, sum, err, npixels, EMPTY_DBL(), EMPTY_DBL());
}
}
/*
* Sum up all the unmasked detector pixels.
*
* @param rebinnedWS: workspace where all the wavelength bins have been grouped
*together
* @param sum: sum of all the unmasked detector pixels (counts)
* @param error: error on sum (counts)
* @param nPixels: number of unmasked detector pixels that contributed to sum
*/
void CalculateEfficiency::sumUnmaskedDetectors(const MatrixWorkspace_sptr &rebinnedWS, double &sum, double &error,
int &nPixels) {
// Number of spectra
const auto numberOfSpectra = static_cast<int>(rebinnedWS->getNumberHistograms());
sum = 0.0;
error = 0.0;
nPixels = 0;
const auto &spectrumInfo = rebinnedWS->spectrumInfo();
for (int i = 0; i < numberOfSpectra; i++) {
progress(0.2 + 0.2 * i / numberOfSpectra, "Computing sensitivity");
// Skip if we have a monitor or if the detector is masked.
if (spectrumInfo.isMonitor(i) || spectrumInfo.isMasked(i))
continue;
// Retrieve the spectrum into a vector
auto &YValues = rebinnedWS->y(i);
auto &YErrors = rebinnedWS->e(i);
sum += YValues[0];
error += YErrors[0] * YErrors[0];
nPixels++;
}
g_log.debug() << "sumUnmaskedDetectors: unmasked pixels = " << nPixels << " from a total of " << numberOfSpectra
<< "\n";
error = std::sqrt(error);
}
/*
* Normalize each detector to produce the relative detector efficiency.
* Pixels that fall outside those efficiency limits will be masked in both
* the input and output workspaces.
*
* @param rebinnedWS: input workspace
* @param outputWS: output workspace containing the relative efficiency
* @param sum: sum of all the unmasked detector pixels (counts)
* @param error: error on sum (counts)
* @param nPixels: number of unmasked detector pixels that contributed to sum
*/
void CalculateEfficiency::normalizeDetectors(const MatrixWorkspace_sptr &rebinnedWS,
const MatrixWorkspace_sptr &outputWS, double sum, double error,
int nPixels, double min_eff, double max_eff) {
// Number of spectra
const size_t numberOfSpectra = rebinnedWS->getNumberHistograms();
// Empty vector to store the pixels that outside the acceptable efficiency
// range
std::vector<size_t> dets_to_mask;
const auto &spectrumInfo = rebinnedWS->spectrumInfo();
for (size_t i = 0; i < numberOfSpectra; i++) {
const double currProgress = 0.4 + 0.2 * (static_cast<double>(i) / static_cast<double>(numberOfSpectra));
progress(currProgress, "Computing sensitivity");
// If this spectrum is masked, skip to the next one
if (spectrumInfo.isMasked(i))
continue;
// Retrieve the spectrum into a vector
auto &YIn = rebinnedWS->y(i);
auto &EIn = rebinnedWS->e(i);
auto &YOut = outputWS->mutableY(i);
auto &EOut = outputWS->mutableE(i);
// If this detector is a monitor, skip to the next one
if (spectrumInfo.isMonitor(i)) {
YOut[0] = 1.0;
EOut[0] = 0.0;
continue;
}
// Normalize counts to get relative efficiency
YOut[0] = nPixels / sum * YIn[0];
const double err_sum = YIn[0] / sum * error;
EOut[0] = nPixels / std::abs(sum) * std::sqrt(EIn[0] * EIn[0] + err_sum * err_sum);
// Mask this detector if the signal is outside the acceptable band
if (!isEmpty(min_eff) && YOut[0] < min_eff)
dets_to_mask.emplace_back(i);
if (!isEmpty(max_eff) && YOut[0] > max_eff)
dets_to_mask.emplace_back(i);
}
g_log.debug() << "normalizeDetectors: Masked pixels outside the acceptable "
"efficiency range ["
<< min_eff << "," << max_eff << "] = " << dets_to_mask.size()
<< " from a total of non masked = " << nPixels
<< " (from a total number of spectra in the ws = " << numberOfSpectra << ")\n";
// If we identified pixels to be masked, mask them now
if (!dets_to_mask.empty()) {
// Mask detectors that were found to be outside the acceptable efficiency
// band
try {
auto mask = createChildAlgorithm("MaskDetectors", 0.8, 0.9);
// First we mask detectors in the output workspace
mask->setProperty<MatrixWorkspace_sptr>("Workspace", outputWS);
mask->setProperty<std::vector<size_t>>("WorkspaceIndexList", dets_to_mask);
mask->execute();
mask = createChildAlgorithm("MaskDetectors", 0.9, 1.0);
// Then we mask the same detectors in the input workspace
mask->setProperty<MatrixWorkspace_sptr>("Workspace", rebinnedWS);
mask->setProperty<std::vector<size_t>>("WorkspaceIndexList", dets_to_mask);
mask->execute();
} catch (std::invalid_argument &err) {
std::stringstream e;
e << "Invalid argument to MaskDetectors Child Algorithm: " << err.what();
g_log.error(e.str());
} catch (std::runtime_error &err) {
std::stringstream e;
e << "Unable to successfully run MaskDetectors Child Algorithm: " << err.what();
g_log.error(e.str());
}
}
}
/**
* Fully masks one component named componentName
* @param ws :: workspace with the respective instrument assigned
* @param componentName :: must be a known CompAssembly.
*/
void CalculateEfficiency::maskComponent(MatrixWorkspace &ws, const std::string &componentName) {
auto instrument = ws.getInstrument();
try {
std::shared_ptr<const Geometry::ICompAssembly> component =
std::dynamic_pointer_cast<const Geometry::ICompAssembly>(instrument->getComponentByName(componentName));
if (!component) {
g_log.warning("Component " + componentName + " expected to be a CompAssembly, e.g., a bank. Component " +
componentName + " not masked!");
return;
}
std::vector<detid_t> detectorList;
for (int x = 0; x < component->nelements(); x++) {
std::shared_ptr<Geometry::ICompAssembly> xColumn =
std::dynamic_pointer_cast<Geometry::ICompAssembly>((*component)[x]);
for (int y = 0; y < xColumn->nelements(); y++) {
std::shared_ptr<Geometry::Detector> detector = std::dynamic_pointer_cast<Geometry::Detector>((*xColumn)[y]);
if (detector) {
auto detID = detector->getID();
detectorList.emplace_back(detID);
}
}
}
auto indexList = ws.getIndicesFromDetectorIDs(detectorList);
auto &spectrumInfo = ws.mutableSpectrumInfo();
for (const auto &idx : indexList) {
ws.getSpectrum(idx).clearData();
spectrumInfo.setMasked(idx, true);
}
} catch (std::exception &) {
g_log.warning("Expecting the component " + componentName +
" to be a CompAssembly, e.g., a bank. Component not masked!");
}
}
/**
* Mask edges of a RectangularDetector
* @param ws :: Input workspace
* @param left :: number of columns to mask left
* @param right :: number of columns to mask right
* @param high :: number of rows to mask top
* @param low :: number of rows to mask Bottom
* @param componentName :: Must be a RectangularDetector
*/
void CalculateEfficiency::maskEdges(const MatrixWorkspace_sptr &ws, int left, int right, int high, int low,
const std::string &componentName) {
auto instrument = ws->getInstrument();
std::shared_ptr<Mantid::Geometry::RectangularDetector> component;
try {
component = std::const_pointer_cast<Mantid::Geometry::RectangularDetector>(
std::dynamic_pointer_cast<const Mantid::Geometry::RectangularDetector>(
instrument->getComponentByName(componentName)));
} catch (std::exception &) {
g_log.warning("Expecting the component " + componentName + " to be a RectangularDetector. maskEdges not executed.");
return;
}
if (!component) {
g_log.warning("Component " + componentName + " is not a RectangularDetector. MaskEdges not executed.");
return;
}
std::vector<int> IDs;
int i = 0;
while (i < left * component->idstep()) {
IDs.emplace_back(component->idstart() + i);
i += 1;
}
// right
i = component->maxDetectorID() - right * component->idstep();
while (i < component->maxDetectorID()) {
IDs.emplace_back(i);
i += 1;
}
// low: 0,256,512,768,..,1,257,513
for (int row = 0; row < low; row++) {
i = row + component->idstart();
while (i < component->nelements() * component->idstep() - component->idstep() + low + component->idstart()) {
IDs.emplace_back(i);
i += component->idstep();
}
}
// high # 255, 511, 767..
for (int row = 0; row < high; row++) {
i = component->idstep() + component->idstart() - row - 1;
while (i < component->nelements() * component->idstep() + component->idstart()) {
IDs.emplace_back(i);
i += component->idstep();
}
}
g_log.debug() << "CalculateEfficiency::maskEdges Detector Ids to Mask:" << std::endl;
for (auto id : IDs) {
g_log.debug() << id << " ";
}
g_log.debug() << std::endl;
auto maskAlg = createChildAlgorithm("MaskDetectors");
maskAlg->setChild(true);
maskAlg->setProperty("Workspace", ws);
maskAlg->setProperty("DetectorList", IDs);
maskAlg->execute();
}
} // namespace Mantid::Algorithms