Mantid
Loading...
Searching...
No Matches
MedianDetectorTest.cpp
Go to the documentation of this file.
1// Mantid Repository : https://github.com/mantidproject/mantid
2//
3// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
4// NScD Oak Ridge National Laboratory, European Spallation Source,
5// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
6// SPDX - License - Identifier: GPL - 3.0 +
11
12#include <cmath>
13
14namespace Mantid::Algorithms {
15// Register the class into the algorithm factory
16DECLARE_ALGORITHM(MedianDetectorTest)
17
18using namespace Kernel;
19using namespace API;
21using namespace Geometry;
22
25 : DetectorDiagnostic(), m_inputWS(), m_loFrac(0.1), m_hiFrac(1.5), m_minWsIndex(0), m_maxWsIndex(EMPTY_INT()),
26 m_rangeLower(0.0), m_rangeUpper(0.0), m_solidAngle(false) {}
27
28const std::string MedianDetectorTest::category() const { return "Diagnostics"; }
29
32 declareProperty(std::make_unique<WorkspaceProperty<>>("InputWorkspace", "", Direction::Input,
33 std::make_shared<HistogramValidator>()),
34 "Name of the input workspace");
35 declareProperty(std::make_unique<WorkspaceProperty<>>("OutputWorkspace", "", Direction::Output),
36 "A MaskWorkspace where 0 denotes a masked spectra. Any spectra containing"
37 "a zero is also masked on the output");
38
39 auto mustBePositive = std::make_shared<BoundedValidator<double>>();
40 mustBePositive->setLower(0);
41 auto mustBePosInt = std::make_shared<BoundedValidator<int>>();
42 mustBePosInt->setLower(0);
43 declareProperty("LevelsUp", 0, mustBePosInt,
44 "Levels above pixel that will be used to compute the median.\n"
45 "If no level is specified, or 0, the median is over the whole "
46 "instrument.");
47 declareProperty("SignificanceTest", 3.3, mustBePositive,
48 "Error criterion as a multiple of error bar i.e. to fail the "
49 "test, the magnitude of the\n"
50 "difference with respect to the median value must also "
51 "exceed this number of error bars");
52 declareProperty("LowThreshold", 0.1, "Lower acceptable bound as fraction of median value");
53 declareProperty("HighThreshold", 1.5, "Upper acceptable bound as fraction of median value");
54 declareProperty("LowOutlier", 0.01, "Lower bound defining outliers as fraction of median value");
55 declareProperty("HighOutlier", 100., "Upper bound defining outliers as fraction of median value");
56 declareProperty("ExcludeZeroesFromMedian", false,
57 "If false (default) zeroes will be included in "
58 "the median calculation, otherwise they will not be included "
59 "but they will be left unmasked");
60
61 declareProperty("StartWorkspaceIndex", 0, mustBePosInt,
62 "The index number of the first spectrum to include in the calculation\n"
63 "(default 0)");
64 declareProperty("EndWorkspaceIndex", EMPTY_INT(), mustBePosInt,
65 "The index number of the last spectrum to include in the calculation\n"
66 "(default the last histogram)");
67 declareProperty("RangeLower", EMPTY_DBL(),
68 "No bin with a boundary at an x value less than this will be included\n"
69 "in the summation used to decide if a detector is 'bad' (default: the\n"
70 "start of each histogram)");
71 declareProperty("RangeUpper", EMPTY_DBL(),
72 "No bin with a boundary at an x value higher than this value will\n"
73 "be included in the summation used to decide if a detector is 'bad'\n"
74 "(default: the end of each histogram)");
75 declareProperty("CorrectForSolidAngle", false, "Flag to correct for solid angle efficiency. False by default.");
76 declareProperty("NumberOfFailures", 0, Direction::Output);
77}
78
87 // Ensures we have a workspace with a single bin. It will contain any input
88 // masking and will be used to record any
89 // required masking from this algorithm
90 MatrixWorkspace_sptr countsWS =
92
93 // 0. Correct for solid angle, if desired
94 if (m_solidAngle) {
96 if (solidAngle != nullptr) {
97 countsWS = countsWS / solidAngle;
98 }
99 }
100
101 // create a vector of vectors of specIDs that will be used to calculate
102 // medians
103 std::vector<std::vector<size_t>> specmap = makeMap(countsWS);
104 const bool excludeZeroes = getProperty("ExcludeZeroesFromMedian");
105 MaskWorkspace_sptr maskWS;
106
107 // 1. Calculate the median
108 std::vector<double> median = calculateMedian(*countsWS, excludeZeroes, specmap);
109 std::vector<double>::iterator medit;
110 for (medit = median.begin(); medit != median.end(); ++medit) {
111 g_log.debug() << "Median value = " << (*medit) << "\n";
112 }
113 // 2. Mask outliers
114 int numFailed = maskOutliers(median, countsWS, specmap);
115
116 // 3. Recalulate the median
117 median = calculateMedian(*countsWS, excludeZeroes, specmap);
118 for (medit = median.begin(); medit != median.end(); ++medit) {
119 g_log.information() << "Median value with outliers removed = " << (*medit) << "\n";
120 }
121
122 maskWS = this->generateEmptyMask(countsWS);
123 numFailed += doDetectorTests(countsWS, median, specmap, maskWS);
124 g_log.information() << "Median test results:\n"
125 << "\tNumber of failures - " << numFailed << "\n";
126 setProperty("NumberOfFailures", numFailed);
127
128 setProperty("OutputWorkspace", maskWS);
129}
130
137 m_inputWS = getProperty("InputWorkspace");
138 int maxWsIndex = static_cast<int>(m_inputWS->getNumberHistograms()) - 1;
139
140 m_parents = getProperty("LevelsUp");
141 m_minWsIndex = getProperty("StartWorkspaceIndex");
142 if ((m_minWsIndex < 0) || (m_minWsIndex > maxWsIndex)) {
143 g_log.warning("StartSpectrum out of range, changed to 0");
144 m_minWsIndex = 0;
145 }
146 m_maxWsIndex = getProperty("EndWorkspaceIndex");
147 if (m_maxWsIndex == EMPTY_INT())
148 m_maxWsIndex = maxWsIndex;
149 if ((m_maxWsIndex < 0) || (m_maxWsIndex > maxWsIndex)) {
150 g_log.warning("EndSpectrum out of range, changed to max spectrum number");
151 m_maxWsIndex = maxWsIndex;
152 }
153 if ((m_maxWsIndex < m_minWsIndex)) {
154 g_log.warning("EndSpectrum can not be less than the StartSpectrum, changed "
155 "to max spectrum number");
156 m_maxWsIndex = maxWsIndex;
157 }
158
159 m_loFrac = getProperty("LowThreshold");
160 m_hiFrac = getProperty("HighThreshold");
161 if (m_loFrac > m_hiFrac) {
162 throw std::invalid_argument("The threshold for reading high must be "
163 "greater than the low threshold");
164 }
165
166 // Integration Range
167 m_rangeLower = getProperty("RangeLower");
168 m_rangeUpper = getProperty("RangeUpper");
169
170 // Solid angle correction flag
171 m_solidAngle = getProperty("CorrectForSolidAngle");
172}
173
183 g_log.debug("Calculating solid angles");
184 // get percentage completed estimates for now, t0 and when we've finished t1
185 double t0 = m_fracDone, t1 = advanceProgress(RTGetSolidAngle);
186 auto childAlg = createChildAlgorithm("SolidAngle", t0, t1, true);
187 childAlg->setProperty("InputWorkspace", m_inputWS);
188 childAlg->setProperty("StartWorkspaceIndex", firstSpec);
189 childAlg->setProperty("EndWorkspaceIndex", lastSpec);
190 try {
191 // Execute the Child Algorithm, it could throw a runtime_error at this point
192 // which would abort execution
193 childAlg->execute();
194 if (!childAlg->isExecuted()) {
195 throw std::runtime_error("Unexpected problem calculating solid angles");
196 }
197 }
198 // catch all exceptions because the solid angle calculation is optional
199 catch (std::exception &) {
200 g_log.warning("Precision warning: Can't find detector geometry " + name() +
201 " will continue with the solid angles of all spectra set to "
202 "the same value");
204 // The return is an empty workspace pointer, which must be handled by the
205 // calling function
207 // function returns normally
208 return empty;
209 }
210 return childAlg->getProperty("OutputWorkspace");
211}
212
220int MedianDetectorTest::maskOutliers(const std::vector<double> &medianvec, const API::MatrixWorkspace_sptr &countsWS,
221 std::vector<std::vector<size_t>> indexmap) {
222
223 // Fractions of the median
224 const double out_lo = getProperty("LowOutlier");
225 const double out_hi = getProperty("HighOutlier");
226
227 int numFailed(0);
228
229 bool checkForMask = false;
230 Geometry::Instrument_const_sptr instrument = countsWS->getInstrument();
231 if (instrument != nullptr) {
232 checkForMask = ((instrument->getSource() != nullptr) && (instrument->getSample() != nullptr));
233 }
234 auto &spectrumInfo = countsWS->mutableSpectrumInfo();
235
236 for (size_t i = 0; i < indexmap.size(); ++i) {
237 std::vector<size_t> &hists = indexmap[i];
238 double median = medianvec[i];
239
241 for (int j = 0; j < static_cast<int>(hists.size()); ++j) { // NOLINT
242 const double value = countsWS->y(hists[j])[0];
243 if ((value == 0.) && checkForMask) {
244 if (spectrumInfo.hasDetectors(hists[j]) && spectrumInfo.isMasked(hists[j])) {
245 numFailed -= 1; // it was already masked
246 }
247 }
248 if ((value < out_lo * median) && (value > 0.0)) {
249 countsWS->getSpectrum(hists[j]).clearData();
250 PARALLEL_CRITICAL(setMasked) {
251 spectrumInfo.setMasked(hists[j], true);
252 ++numFailed;
253 }
254 } else if (value > out_hi * median) {
255 countsWS->getSpectrum(hists[j]).clearData();
256 PARALLEL_CRITICAL(setMasked) {
257 spectrumInfo.setMasked(hists[j], true);
258 ++numFailed;
259 }
260 }
261 }
263 }
264
265 return numFailed;
266}
267
279int MedianDetectorTest::doDetectorTests(const API::MatrixWorkspace_sptr &countsWS, const std::vector<double> &medianvec,
280 std::vector<std::vector<size_t>> indexmap,
281 const API::MatrixWorkspace_sptr &maskWS) {
282 g_log.debug("Applying the criteria to find failing detectors");
283
284 // A spectra can't fail if the statistics show its value is consistent with
285 // the mean value,
286 // check the error and how many errorbars we are away
287 const double minSigma = getProperty("SignificanceTest");
288
289 // prepare to report progress
290 const int numSpec(m_maxWsIndex - m_minWsIndex);
291 const auto progStep = static_cast<int>(ceil(numSpec / 30.0));
292 int steps(0);
293
294 const double deadValue(1.0);
295 int numFailed(0);
296
297 bool checkForMask = false;
298 Geometry::Instrument_const_sptr instrument = countsWS->getInstrument();
299 if (instrument != nullptr) {
300 checkForMask = ((instrument->getSource() != nullptr) && (instrument->getSample() != nullptr));
301 }
302 const auto &spectrumInfo = countsWS->spectrumInfo();
303
304 PARALLEL_FOR_IF(Kernel::threadSafe(*countsWS, *maskWS))
305 for (int j = 0; j < static_cast<int>(indexmap.size()); ++j) {
306 std::vector<size_t> hists = indexmap.at(j);
307 double median = medianvec.at(j);
308 const size_t nhist = hists.size();
309 g_log.debug() << "new component with " << nhist << " spectra.\n";
310 for (size_t i = 0; i < nhist; ++i) {
311 g_log.debug() << "Counts workspace index=" << i << ", Mask workspace index=" << hists.at(i) << '\n';
313 ++steps;
314 // update the progressbar information
315 if (steps % progStep == 0) {
316 progress(advanceProgress(progStep * static_cast<double>(RTMarkDetects) / numSpec));
317 }
318
319 if (checkForMask && spectrumInfo.hasDetectors(hists.at(i))) {
320 if (spectrumInfo.isMasked(hists.at(i))) {
321 maskWS->mutableY(hists.at(i))[0] = deadValue;
322 continue;
323 }
324 if (spectrumInfo.isMonitor(hists.at(i))) {
325 // Don't include in calculation but don't mask it
326 continue;
327 }
328 }
329
330 const double signal = countsWS->y(hists.at(i))[0];
331
332 // Mask out NaN and infinite
333 if (!std::isfinite(signal)) {
334 maskWS->mutableY(hists.at(i))[0] = deadValue;
336 ++numFailed;
337 continue;
338 }
339
340 const double error = minSigma * countsWS->e(hists.at(i))[0];
341
342 if ((signal < median * m_loFrac && (signal - median < -error)) ||
343 (signal > median * m_hiFrac && (signal - median > error))) {
344 maskWS->mutableY(hists.at(i))[0] = deadValue;
346 ++numFailed;
347 }
348
350 }
352
353 // Log finds
354 g_log.information() << numFailed << " spectra failed the median tests.\n";
355 }
356 return numFailed;
357}
358
359} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
double value
The value of the point.
Definition: FitMW.cpp:51
double error
Definition: IndexPeaks.cpp:133
#define PARALLEL_START_INTERRUPT_REGION
Begins a block to skip processing is the algorithm has been interupted Note the end of the block if n...
Definition: MultiThreaded.h:94
#define PARALLEL_CRITICAL(name)
#define PARALLEL_END_INTERRUPT_REGION
Ends a block to skip processing is the algorithm has been interupted Note the start of the block if n...
#define PARALLEL_ATOMIC
#define PARALLEL_FOR_IF(condition)
Empty definitions - to enable set your complier to enable openMP.
#define PARALLEL_CHECK_INTERRUPT_REGION
Adds a check after a Parallel region to see if it was interupted.
void declareProperty(std::unique_ptr< Kernel::Property > p, const std::string &doc="") override
Add a property to the list of managed properties.
Definition: Algorithm.cpp:1913
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
Definition: Algorithm.cpp:2076
virtual std::shared_ptr< Algorithm > createChildAlgorithm(const std::string &name, const double startProgress=-1., const double endProgress=-1., const bool enableLogging=true, const int &version=-1)
Create a Child Algorithm.
Definition: Algorithm.cpp:842
Kernel::Logger & g_log
Definition: Algorithm.h:451
void progress(double p, const std::string &msg="", double estimatedTime=0.0, int progressPrecision=0)
Sends ProgressNotification.
Definition: Algorithm.cpp:231
A property class for workspaces.
A base class for a detector diagnostic algorithm.
std::vector< double > calculateMedian(const API::MatrixWorkspace &input, bool excludeZeroes, const std::vector< std::vector< size_t > > &indexmap)
Calculate the median of the given workspace.
std::vector< std::vector< size_t > > makeMap(const API::MatrixWorkspace_sptr &countsWS)
method to check which spectra should be grouped when calculating the median
double m_fracDone
An estimate of the percentage of the algorithm runtimes that has been completed.
@ RTGetSolidAngle
An estimate of how much work SolidAngle will do for each spectrum.
@ RTMarkDetects
Time taken to find failing detectors.
API::MatrixWorkspace_sptr integrateSpectra(const API::MatrixWorkspace_sptr &inputWS, const int indexMin, const int indexMax, const double lower, const double upper, const bool outputWorkspace2D=false)
Get the total counts for each spectra.
DataObjects::MaskWorkspace_sptr generateEmptyMask(const API::MatrixWorkspace_const_sptr &inputWS)
Create a masking workspace to return.
double advanceProgress(double toAdd)
Update the fraction complete estimate assuming that the algorithm has completed a task with estimated...
int m_parents
number of parents up, 0 go to instrument
void failProgress(RunTime aborted)
Update the fraction complete estimate assuming that the algorithm aborted a task with estimated RunTi...
double m_rangeLower
Start point for integration.
MedianDetectorTest()
Default constructor initialises all values to zero and runs the base class constructor.
void exec() override
Executes the algorithm that includes calls to SolidAngle and Integration.
double m_rangeUpper
End point for integration.
int m_minWsIndex
The index of the first spectrum to calculate.
const std::string name() const override
Algorithm's name for identification overriding a virtual method.
double m_loFrac
The proportion of the median value below which a detector is considered under-reading.
const std::string category() const override
Algorithm's category for identification.
API::MatrixWorkspace_sptr getSolidAngles(int firstSpec, int lastSpec)
Calculates the sum of solid angles of detectors for each histogram.
int m_maxWsIndex
The index of the last spectrum to calculate.
void init() override
Declare algorithm properties.
void retrieveProperties()
Loads and checks the values passed to the algorithm.
double m_hiFrac
The factor of the median value above which a detector is considered over-reading.
int doDetectorTests(const API::MatrixWorkspace_sptr &countsWS, const std::vector< double > &medianvec, std::vector< std::vector< size_t > > indexmap, const API::MatrixWorkspace_sptr &maskWS)
Do the tests and mask those that fail.
int maskOutliers(const std::vector< double > &medianvec, const API::MatrixWorkspace_sptr &countsWS, std::vector< std::vector< size_t > > indexmap)
Mask the outlier values to get a better median value.
bool m_solidAngle
flag for solid angle correction
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void debug(const std::string &msg)
Logs at debug level.
Definition: Logger.cpp:114
void warning(const std::string &msg)
Logs at warning level.
Definition: Logger.cpp:86
void information(const std::string &msg)
Logs at information level.
Definition: Logger.cpp:105
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
std::shared_ptr< MaskWorkspace > MaskWorkspace_sptr
shared pointer to the MaskWorkspace class
Definition: MaskWorkspace.h:64
std::shared_ptr< const Instrument > Instrument_const_sptr
Shared pointer to an const instrument object.
std::enable_if< std::is_pointer< Arg >::value, bool >::type threadSafe(Arg workspace)
Thread-safety check Checks the workspace to ensure it is suitable for multithreaded access.
Definition: MultiThreaded.h:22
constexpr int EMPTY_INT() noexcept
Returns what we consider an "empty" integer within a property.
Definition: EmptyValues.h:25
constexpr double EMPTY_DBL() noexcept
Returns what we consider an "empty" double within a property.
Definition: EmptyValues.h:43
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54