Mantid
Loading...
Searching...
No Matches
DetectorDiagnostic.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 +
18
19#include <boost/iterator/counting_iterator.hpp>
20
21#include <algorithm>
22#include <cfloat>
23
24namespace Mantid::Algorithms {
25// Register the class into the algorithm factory
26DECLARE_ALGORITHM(DetectorDiagnostic)
27
31using std::string;
32using namespace Mantid::DataObjects;
33using namespace Mantid::API;
34using namespace Mantid::Kernel;
35
36//--------------------------------------------------------------------------
37// Functions to make this a proper workflow algorithm
38//--------------------------------------------------------------------------
39const std::string DetectorDiagnostic::category() const { return "Diagnostics;Workflow\\Diagnostics"; }
40
41const std::string DetectorDiagnostic::name() const { return "DetectorDiagnostic"; }
42
43int DetectorDiagnostic::version() const { return 1; }
44
46 this->declareProperty(std::make_unique<WorkspaceProperty<>>("InputWorkspace", "", Direction::Input),
47 "Name of the integrated detector vanadium (white beam) workspace.");
48 this->declareProperty(
49 std::make_unique<WorkspaceProperty<>>("HardMaskWorkspace", "", Direction::Input, PropertyMode::Optional),
50 "A hard mask to apply to the inputworkspace");
51 this->declareProperty(std::make_unique<WorkspaceProperty<>>("OutputWorkspace", "", Direction::Output),
52 "A MaskWorkspace containing the masked spectra as zeroes and ones.");
53 auto mustBePosInt = std::make_shared<BoundedValidator<int>>();
54 mustBePosInt->setLower(0);
55 this->declareProperty("StartWorkspaceIndex", 0, mustBePosInt,
56 "The index number of the first spectrum to include in the calculation\n"
57 "(default 0)");
58 this->declareProperty("EndWorkspaceIndex", EMPTY_INT(), mustBePosInt,
59 "The index number of the last spectrum to include in the calculation\n"
60 "(default the last histogram)");
61 this->declareProperty("RangeLower", EMPTY_DBL(),
62 "No bin with a boundary at an x value less than this will be used\n"
63 "in the summation that decides if a detector is 'bad' (default: the\n"
64 "start of each histogram)");
65 this->declareProperty("RangeUpper", EMPTY_DBL(),
66 "No bin with a boundary at an x value higher than this value will\n"
67 "be used in the summation that decides if a detector is 'bad'\n"
68 "(default: the end of each histogram)");
69
70 string findDetOutLimGrp("Find Detectors Outside Limits");
71 this->declareProperty("LowThreshold", 0.0,
72 "Spectra whose total number of counts are equal to or below this value\n"
73 "will be marked bad (default 0)");
74 this->setPropertyGroup("LowThreshold", findDetOutLimGrp);
75 this->declareProperty("HighThreshold", EMPTY_DBL(),
76 "Spectra whose total number of counts are equal to or above this value\n"
77 "will be marked bad (default off)");
78 this->setPropertyGroup("HighThreshold", findDetOutLimGrp);
79
80 string medianDetTestGrp("Median Detector Test");
81 auto mustBePositiveDbl = std::make_shared<BoundedValidator<double>>();
82 mustBePositiveDbl->setLower(0);
83 this->declareProperty("LevelsUp", 0, mustBePosInt,
84 "Levels above pixel that will be used to compute the median.\n"
85 "If no level is specified, or 0, the median is over the whole "
86 "instrument.");
87 this->setPropertyGroup("LevelsUp", medianDetTestGrp);
88 this->declareProperty("SignificanceTest", 0.0, mustBePositiveDbl,
89 "Error criterion as a multiple of error bar i.e. to "
90 "fail the test, the magnitude of the\n"
91 "difference with respect to the median value must also "
92 "exceed this number of error bars");
93 this->setPropertyGroup("SignificanceTest", medianDetTestGrp);
94 this->declareProperty("LowThresholdFraction", 0.1, "Lower acceptable bound as fraction of median value");
95 this->setPropertyGroup("LowThresholdFraction", medianDetTestGrp);
96 this->declareProperty("HighThresholdFraction", 1.5, "Upper acceptable bound as fraction of median value");
97 this->setPropertyGroup("HighThresholdFraction", medianDetTestGrp);
98 this->declareProperty("LowOutlier", 0.01, "Lower bound defining outliers as fraction of median value");
99 this->setPropertyGroup("LowOutlier", medianDetTestGrp);
100 this->declareProperty("HighOutlier", 100., "Upper bound defining outliers as fraction of median value");
101 this->setPropertyGroup("HighOutlier", medianDetTestGrp);
102 this->declareProperty("CorrectForSolidAngle", false, "Flag to correct for solid angle efficiency. False by default.");
103 this->setPropertyGroup("CorrectForSolidAngle", medianDetTestGrp);
104 this->declareProperty("ExcludeZeroesFromMedian", false,
105 "If false (default) zeroes will be included in "
106 "the median calculation, otherwise they will not be "
107 "included but they will be left unmasked");
108 this->setPropertyGroup("ExcludeZeroesFromMedian", medianDetTestGrp);
109
110 string detEffVarGrp("Detector Efficiency Variation");
111 this->declareProperty(
112 std::make_unique<WorkspaceProperty<>>("DetVanCompare", "", Direction::Input, PropertyMode::Optional),
113 "Name of a matching second detector vanadium run from the same\n"
114 "instrument. It must be treated in the same manner as the input detector "
115 "vanadium.");
116 this->setPropertyGroup("DetVanCompare", detEffVarGrp);
117 this->declareProperty("DetVanRatioVariation", 1.1, mustBePositiveDbl,
118 "Identify spectra whose total number of counts has changed by more\n"
119 "than this factor of the median change between the two input workspaces");
120 this->setPropertyGroup("DetVanRatioVariation", detEffVarGrp);
121 this->setPropertySettings("DetVanRatioVariation",
122 std::make_unique<EnabledWhenProperty>("DetVanCompare", IS_NOT_DEFAULT));
123
124 string countsCheck("Check Sample Counts");
125 this->declareProperty(
126 std::make_unique<WorkspaceProperty<>>("SampleTotalCountsWorkspace", "", Direction::Input, PropertyMode::Optional),
127 "A sample workspace integrated over the full axis range.");
128 this->setPropertyGroup("SampleTotalCountsWorkspace", countsCheck);
129
130 string backgroundCheck("Check Sample Background");
131 this->declareProperty(
132 std::make_unique<WorkspaceProperty<>>("SampleBackgroundWorkspace", "", Direction::Input, PropertyMode::Optional),
133 "A sample workspace integrated over the background region.");
134 this->setPropertyGroup("SampleBackgroundWorkspace", backgroundCheck);
135 this->declareProperty("SampleBkgLowAcceptanceFactor", 0.0, mustBePositiveDbl,
136 "Low threshold for the background check MedianDetectorTest.");
137 this->setPropertyGroup("SampleBkgLowAcceptanceFactor", backgroundCheck);
138 this->declareProperty("SampleBkgHighAcceptanceFactor", 5.0, mustBePositiveDbl,
139 "High threshold for the background check MedianDetectorTest.");
140 this->setPropertyGroup("SampleBkgHighAcceptanceFactor", backgroundCheck);
141 this->declareProperty("SampleBkgSignificanceTest", 3.3, mustBePositiveDbl,
142 "Error criterion as a multiple of error bar i.e. to "
143 "fail the test, the magnitude of the\n"
144 "difference with respect to the median value must also "
145 "exceed this number of error bars");
146 this->setPropertyGroup("SampleBkgSignificanceTest", backgroundCheck);
147 this->declareProperty("SampleCorrectForSolidAngle", false,
148 "Flag to correct for solid angle efficiency for "
149 "background check MedianDetectorTest. False by "
150 "default.");
151 this->setPropertyGroup("SampleCorrectForSolidAngle", backgroundCheck);
152
153 string psdBleedMaskGrp("Create PSD Bleed Mask");
154 this->declareProperty(
155 std::make_unique<WorkspaceProperty<>>("SampleWorkspace", "", Direction::Input, PropertyMode::Optional),
156 "A sample workspace. This is used in the PSD Bleed calculation.");
157 this->setPropertyGroup("SampleWorkspace", psdBleedMaskGrp);
158 this->declareProperty("MaxTubeFramerate", 0.0, mustBePositiveDbl,
159 "The maximum rate allowed for a tube in counts/us/frame.");
160 this->setPropertyGroup("MaxTubeFramerate", psdBleedMaskGrp);
161 this->declareProperty("NIgnoredCentralPixels", 80, mustBePosInt, "The number of pixels about the centre to ignore.");
162 this->setPropertyGroup("NIgnoredCentralPixels", psdBleedMaskGrp);
163 this->setPropertySettings("NIgnoredCentralPixels",
164 std::make_unique<EnabledWhenProperty>("MaxTubeFramerate", IS_NOT_DEFAULT));
165
166 this->declareProperty("NumberOfFailures", 0, Direction::Output);
167}
168
170 // get the generic information that everybody uses
171 MatrixWorkspace_sptr inputWS = this->getProperty("InputWorkspace");
172 m_minIndex = this->getProperty("StartWorkspaceIndex");
173 m_maxIndex = this->getProperty("EndWorkspaceIndex");
174 m_rangeLower = this->getProperty("RangeLower");
175 m_rangeUpper = this->getProperty("RangeUpper");
176
177 // integrate the data once to pass to ChildAlgorithms
178 m_fracDone = 0.;
179
180 // Get the other workspaces
181 MatrixWorkspace_sptr input2WS = this->getProperty("DetVanCompare");
182 MatrixWorkspace_sptr totalCountsWS = this->getProperty("SampleTotalCountsWorkspace");
183 MatrixWorkspace_sptr bkgWS = this->getProperty("SampleBackgroundWorkspace");
184 MatrixWorkspace_sptr sampleWS = this->getProperty("SampleWorkspace");
185
186 // calculate the number of tests for progress bar
187 m_progStepWidth = 0;
188 {
189 int numTests(1); // if detector vanadium present, do it!
190 if (input2WS)
191 numTests += 1;
192 if (totalCountsWS)
193 numTests += 1;
194 if (bkgWS)
195 numTests += 1;
196 if (sampleWS)
197 numTests += 1;
198 g_log.information() << "Number of tests requested: " << numTests << '\n';
199 m_progStepWidth = (1. - m_fracDone) / static_cast<double>(numTests);
200 }
201
202 int numFailed(0);
204
205 // Apply hard mask if present
206 MatrixWorkspace_sptr hardMaskWS = this->getProperty("HardMaskWorkspace");
207 if (hardMaskWS) {
208 auto md = createChildAlgorithm("MaskDetectors");
209 md->setProperty("Workspace", inputWS);
210 md->setProperty("MaskedWorkspace", hardMaskWS);
211 md->executeAsChildAlg();
212 }
213
214 // Perform FindDetectorsOutsideLimits and MedianDetectorTest on the
215 // detector vanadium
216 maskWS = this->doDetVanTest(inputWS, numFailed);
217
218 // DetectorEfficiencyVariation (only if two workspaces are specified)
219 if (input2WS) {
220 // apply mask to what we are going to input
221 this->applyMask(input2WS, maskWS);
222
223 maskWS = this->doDetVanTest(input2WS, numFailed);
224
225 // get the relevant inputs
226 double variation = this->getProperty("DetVanRatioVariation");
227
228 // run the ChildAlgorithm
229 auto alg = createChildAlgorithm("DetectorEfficiencyVariation", m_fracDone, m_fracDone + m_progStepWidth);
231 alg->setProperty("WhiteBeamBase", inputWS);
232 alg->setProperty("WhiteBeamCompare", input2WS);
233 alg->setProperty("StartWorkspaceIndex", m_minIndex);
234 alg->setProperty("EndWorkspaceIndex", m_maxIndex);
235 alg->setProperty("RangeLower", m_rangeLower);
236 alg->setProperty("RangeUpper", m_rangeUpper);
237 alg->setProperty("Variation", variation);
238 alg->executeAsChildAlg();
239 MatrixWorkspace_sptr localMaskWS = alg->getProperty("OutputWorkspace");
240 applyMask(inputWS, localMaskWS);
241 applyMask(input2WS, localMaskWS);
242 int localFails = alg->getProperty("NumberOfFailures");
243 numFailed += localFails;
244 }
245
246 // Zero total counts check for sample counts
247 if (totalCountsWS) {
248 // apply mask to what we are going to input
249 applyMask(totalCountsWS, maskWS);
250
251 auto zeroChk = createChildAlgorithm("FindDetectorsOutsideLimits", m_fracDone, m_fracDone + m_progStepWidth);
253 zeroChk->setProperty("InputWorkspace", totalCountsWS);
254 zeroChk->setProperty("StartWorkspaceIndex", m_minIndex);
255 zeroChk->setProperty("EndWorkspaceIndex", m_maxIndex);
256 zeroChk->setProperty("LowThreshold", 1.0e-10);
257 zeroChk->setProperty("HighThreshold", 1.0e100);
258 zeroChk->executeAsChildAlg();
259 MatrixWorkspace_sptr localMaskWS = zeroChk->getProperty("OutputWorkspace");
260 applyMask(inputWS, localMaskWS);
261 int localFails = zeroChk->getProperty("NumberOfFailures");
262 numFailed += localFails;
263 }
264
265 // Background check
266 if (bkgWS) {
267 // apply mask to what we are going to input
268 this->applyMask(bkgWS, maskWS);
269
270 double significanceTest = this->getProperty("SampleBkgSignificanceTest");
271 double lowThreshold = this->getProperty("SampleBkgLowAcceptanceFactor");
272 double highThreshold = this->getProperty("SampleBkgHighAcceptanceFactor");
273 bool correctSA = this->getProperty("SampleCorrectForSolidAngle");
274
275 // run the ChildAlgorithm
276 auto alg = createChildAlgorithm("MedianDetectorTest", m_fracDone, m_fracDone + m_progStepWidth);
278 alg->setProperty("InputWorkspace", bkgWS);
279 alg->setProperty("StartWorkspaceIndex", m_minIndex);
280 alg->setProperty("EndWorkspaceIndex", m_maxIndex);
281 alg->setProperty("SignificanceTest", significanceTest);
282 alg->setProperty("LowThreshold", lowThreshold);
283 alg->setProperty("HighThreshold", highThreshold);
284 alg->setProperty("LowOutlier", 0.0);
285 alg->setProperty("HighOutlier", 1.0e100);
286 alg->setProperty("ExcludeZeroesFromMedian", true);
287 alg->setProperty("CorrectForSolidAngle", correctSA);
288 alg->executeAsChildAlg();
289 MatrixWorkspace_sptr localMaskWS = alg->getProperty("OutputWorkspace");
290 applyMask(inputWS, localMaskWS);
291 int localFails = alg->getProperty("NumberOfFailures");
292 numFailed += localFails;
293 }
294
295 // CreatePSDBleedMask (if selected)
296 if (sampleWS) {
297 // get the relevant inputs
298 double maxTubeFrameRate = this->getProperty("MaxTubeFramerate");
299 int numIgnore = this->getProperty("NIgnoredCentralPixels");
300
301 // run the ChildAlgorithm
302 auto alg = createChildAlgorithm("CreatePSDBleedMask", m_fracDone, m_fracDone + m_progStepWidth);
304 alg->setProperty("InputWorkspace", sampleWS);
305 alg->setProperty("MaxTubeFramerate", maxTubeFrameRate);
306 alg->setProperty("NIgnoredCentralPixels", numIgnore);
307 alg->executeAsChildAlg();
308 MatrixWorkspace_sptr localMaskWS = alg->getProperty("OutputWorkspace");
309 applyMask(inputWS, localMaskWS);
310 int localFails = alg->getProperty("NumberOfFailures");
311 numFailed += localFails;
312 }
313
314 g_log.information() << numFailed << " spectra are being masked\n";
315 setProperty("NumberOfFailures", numFailed);
316
317 // Extract the mask from the vanadium workspace
318 std::vector<int> detList;
319 auto extract = createChildAlgorithm("ExtractMask");
320 extract->setProperty("InputWorkspace", inputWS);
321 extract->setProperty("OutputWorkspace", "final_mask");
322 extract->setProperty("DetectorList", detList);
323 extract->executeAsChildAlg();
324 maskWS = extract->getProperty("OutputWorkspace");
325
326 this->setProperty("OutputWorkspace", maskWS);
327}
328
335 auto maskAlg = createChildAlgorithm("MaskDetectors"); // should set progress bar
336 maskAlg->setProperty("Workspace", inputWS);
337 maskAlg->setProperty("MaskedWorkspace", maskWS);
338 maskAlg->setProperty("StartWorkspaceIndex", m_minIndex);
339 maskAlg->setProperty("EndWorkspaceIndex", m_maxIndex);
340 maskAlg->executeAsChildAlg();
341}
342
350 MatrixWorkspace_sptr localMask;
351
352 // FindDetectorsOutsideLimits
353 // get the relevant inputs
354 double lowThreshold = this->getProperty("LowThreshold");
355 double highThreshold = this->getProperty("HighThreshold");
356 // run the ChildAlgorithm
357 auto fdol = createChildAlgorithm("FindDetectorsOutsideLimits", m_fracDone, m_fracDone + m_progStepWidth);
359 fdol->setProperty("InputWorkspace", inputWS);
360 fdol->setProperty("OutputWorkspace", localMask);
361 fdol->setProperty("StartWorkspaceIndex", m_minIndex);
362 fdol->setProperty("EndWorkspaceIndex", m_maxIndex);
363 fdol->setProperty("RangeLower", m_rangeLower);
364 fdol->setProperty("RangeUpper", m_rangeUpper);
365 fdol->setProperty("LowThreshold", lowThreshold);
366 fdol->setProperty("HighThreshold", highThreshold);
367 fdol->executeAsChildAlg();
368 localMask = fdol->getProperty("OutputWorkspace");
369 int localFails = fdol->getProperty("NumberOfFailures");
370 nFails += localFails;
371
372 // get the relevant inputs for the MedianDetectorTests
373 int parents = this->getProperty("LevelsUp");
374 double significanceTest = this->getProperty("SignificanceTest");
375 double lowThresholdFrac = this->getProperty("LowThresholdFraction");
376 double highThresholdFrac = this->getProperty("HighThresholdFraction");
377 double lowOutlier = this->getProperty("LowOutlier");
378 double highOutlier = this->getProperty("HighOutlier");
379 bool excludeZeroes = this->getProperty("ExcludeZeroesFromMedian");
380 bool correctforSA = this->getProperty("CorrectForSolidAngle");
381
382 // MedianDetectorTest
383 // apply mask to what we are going to input
384 this->applyMask(inputWS, localMask);
385
386 // run the ChildAlgorithm
387 auto mdt = createChildAlgorithm("MedianDetectorTest", m_fracDone, m_fracDone + m_progStepWidth);
389 mdt->setProperty("InputWorkspace", inputWS);
390 mdt->setProperty("StartWorkspaceIndex", m_minIndex);
391 mdt->setProperty("EndWorkspaceIndex", m_maxIndex);
392 mdt->setProperty("RangeLower", m_rangeLower);
393 mdt->setProperty("RangeUpper", m_rangeUpper);
394 mdt->setProperty("LevelsUp", parents);
395 mdt->setProperty("SignificanceTest", significanceTest);
396 mdt->setProperty("LowThreshold", lowThresholdFrac);
397 mdt->setProperty("HighThreshold", highThresholdFrac);
398 mdt->setProperty("LowOutlier", lowOutlier);
399 mdt->setProperty("HighOutlier", highOutlier);
400 mdt->setProperty("ExcludeZeroesFromMedian", excludeZeroes);
401 mdt->setProperty("CorrectForSolidAngle", correctforSA);
402 mdt->executeAsChildAlg();
403 localMask = mdt->getProperty("OutputWorkspace");
404 localFails = mdt->getProperty("NumberOfFailures");
405 nFails += localFails;
406
407 this->applyMask(inputWS, localMask);
408 return localMask;
409}
410
411//--------------------------------------------------------------------------
412// Public member functions
413//--------------------------------------------------------------------------
415 : API::Algorithm(), m_fracDone(0.0), m_TotalTime(RTTotal), m_parents(0), m_progStepWidth(0.0), m_minIndex(0),
416 m_maxIndex(EMPTY_INT()), m_rangeLower(EMPTY_DBL()), m_rangeUpper(EMPTY_DBL()) {}
417
418//--------------------------------------------------------------------------
419// Protected member functions
420//--------------------------------------------------------------------------
421
434 const int indexMax, const double lower, const double upper,
435 const bool outputWorkspace2D) {
436 g_log.debug() << "Integrating input spectra.\n";
437 // If the input spectra only has one bin, assume it has been integrated
438 // already
439 // but we need to pass it to the algorithm so that a copy of the input
440 // workspace is
441 // actually created to use for further calculations
442 // get percentage completed estimates for now, t0 and when we've finished t1
444 auto childAlg = createChildAlgorithm("Integration", t0, t1);
445 childAlg->setProperty("InputWorkspace", inputWS);
446 childAlg->setProperty("StartWorkspaceIndex", indexMin);
447 childAlg->setProperty("EndWorkspaceIndex", indexMax);
448 // pass inputed values straight to this integration trusting the checking done
449 // there
450 childAlg->setProperty("RangeLower", lower);
451 childAlg->setProperty("RangeUpper", upper);
452 childAlg->setPropertyValue("IncludePartialBins", "1");
453 childAlg->executeAsChildAlg();
454
455 // Convert to 2D if desired, and if the input was an EventWorkspace.
456 MatrixWorkspace_sptr outputW = childAlg->getProperty("OutputWorkspace");
457 MatrixWorkspace_sptr finalOutputW = outputW;
458 if (outputWorkspace2D && std::dynamic_pointer_cast<EventWorkspace>(outputW)) {
459 g_log.debug() << "Converting output Event Workspace into a Workspace2D.\n";
460 childAlg = createChildAlgorithm("ConvertToMatrixWorkspace", t0, t1);
461 childAlg->setProperty("InputWorkspace", outputW);
462 childAlg->executeAsChildAlg();
463 finalOutputW = childAlg->getProperty("OutputWorkspace");
464 }
465
466 return finalOutputW;
467}
468
476 // Create a new workspace for the results, copy from the input to ensure that
477 // we copy over the instrument and current masking
478 auto maskWS = create<DataObjects::MaskWorkspace>(*inputWS, HistogramData::Points(1));
479 maskWS->setTitle(inputWS->getTitle());
480
481 return maskWS;
482}
483
484std::vector<std::vector<size_t>> DetectorDiagnostic::makeInstrumentMap(const API::MatrixWorkspace &countsWS) {
485 return {{boost::counting_iterator<std::size_t>(0),
486 boost::counting_iterator<std::size_t>(countsWS.getNumberHistograms())}};
487}
492std::vector<std::vector<size_t>> DetectorDiagnostic::makeMap(const API::MatrixWorkspace_sptr &countsWS) {
493 std::multimap<Mantid::Geometry::ComponentID, size_t> mymap;
494
495 Geometry::Instrument_const_sptr instrument = countsWS->getInstrument();
496 if (m_parents == 0) {
497 return makeInstrumentMap(*countsWS);
498 }
499 if (!instrument) {
500 g_log.warning("Workspace has no instrument. LevelsUP is ignored");
501 return makeInstrumentMap(*countsWS);
502 }
503
504 // check if not grouped. If grouped, it will throw
505 if (countsWS->hasGroupedDetectors()) {
506 throw std::runtime_error("Median detector test: not able to create "
507 "detector to spectra map. Try with LevelUp=0.");
508 }
509
510 const SpectrumInfo &spectrumInfo = countsWS->spectrumInfo();
511
512 for (size_t i = 0; i < countsWS->getNumberHistograms(); ++i) {
513
514 auto anc = spectrumInfo.detector(i).getAncestors();
515 if (anc.size() < static_cast<size_t>(m_parents)) {
516 g_log.warning("Too many levels up. Will ignore LevelsUp");
517 m_parents = 0;
518 return makeInstrumentMap(*countsWS);
519 }
520 mymap.emplace(anc[m_parents - 1]->getComponentID(), i);
521 }
522
523 std::vector<std::vector<size_t>> speclist;
524
525 std::multimap<Mantid::Geometry::ComponentID, size_t>::iterator m_it, s_it;
526 for (m_it = mymap.begin(); m_it != mymap.end(); m_it = s_it) {
527 Mantid::Geometry::ComponentID theKey = m_it->first;
528 auto keyRange = mymap.equal_range(theKey);
529 // Iterate over all map elements with key == theKey
530 std::vector<size_t> speclistsingle;
531 for (s_it = keyRange.first; s_it != keyRange.second; ++s_it) {
532 speclistsingle.emplace_back(s_it->second);
533 }
534 speclist.emplace_back(std::move(speclistsingle));
535 }
536
537 return speclist;
538}
553std::vector<double> DetectorDiagnostic::calculateMedian(const API::MatrixWorkspace &input, bool excludeZeroes,
554 const std::vector<std::vector<size_t>> &indexmap) {
555 std::vector<double> medianvec;
556 g_log.debug("Calculating the median count rate of the spectra");
557
558 bool checkForMask = false;
560 if (instrument != nullptr) {
561 checkForMask = ((instrument->getSource() != nullptr) && (instrument->getSample() != nullptr));
562 }
563 const auto &spectrumInfo = input.spectrumInfo();
564
565 for (const auto &hists : indexmap) {
566 std::vector<double> medianInput;
567 const auto nhists = static_cast<int>(hists.size());
568 // The maximum possible length is that of workspace length
569 medianInput.reserve(nhists);
570
572 for (int i = 0; i < nhists; ++i) { // NOLINT
574
575 if (checkForMask && spectrumInfo.hasDetectors(hists[i])) {
576 if (spectrumInfo.isMasked(hists[i]) || spectrumInfo.isMonitor(hists[i]))
577 continue;
578 }
579
580 const double yValue = input.readY(hists[i])[0];
581 if (yValue < 0.0) {
582 throw std::out_of_range("Negative number of counts found, could be "
583 "corrupted raw counts or solid angle data");
584 }
585 if (!std::isfinite(yValue) || (excludeZeroes && yValue < DBL_EPSILON)) // NaNs/Infs
586 {
587 continue;
588 }
589 // Now we have a good value
590 PARALLEL_CRITICAL(DetectorDiagnostic_median_d) { medianInput.emplace_back(yValue); }
591
593 }
595
596 if (medianInput.empty()) {
597 g_log.information("some group has no valid histograms. Will use 0 for median.");
598 medianInput.emplace_back(0.);
599 }
600
602 double median = stats.median;
603
604 if (median < 0 || median > DBL_MAX / 10.0) {
605 throw std::out_of_range("The calculated value for the median was either "
606 "negative or unreliably large");
607 }
608 medianvec.emplace_back(median);
609 }
610 return medianvec;
611}
612
619 if (workspace->isDistribution()) {
620 g_log.information() << "Workspace already contains a count rate, nothing to do.\n";
621 return workspace;
622 }
623
624 g_log.information("Calculating time averaged count rates");
625 // get percentage completed estimates for now, t0 and when we've finished t1
626 double t0 = m_fracDone, t1 = advanceProgress(RTGetRate);
627 auto childAlg = createChildAlgorithm("ConvertToDistribution", t0, t1);
628 childAlg->setProperty<MatrixWorkspace_sptr>("Workspace", workspace);
629 // Now execute the Child Algorithm but allow any exception to bubble up
630 childAlg->execute();
631 return childAlg->getProperty("Workspace");
632}
633
642 m_fracDone += toAdd / m_TotalTime;
643 // it could go negative as sometimes the percentage is re-estimated backwards,
644 // this is worrying about if a small negative value will cause a problem some
645 // where
646 m_fracDone = std::abs(m_fracDone);
648 return m_fracDone;
649}
650
659 advanceProgress(-aborted);
660 m_TotalTime -= aborted;
661}
662} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
IPeaksWorkspace_sptr workspace
Definition: IndexPeaks.cpp:114
#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_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.
double lower
lower and upper bounds on the multiplier, if known
double upper
Base class from which all concrete algorithm classes should be derived.
Definition: Algorithm.h:85
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 interruption_point()
This is called during long-running operations, and check if the algorithm has requested that it be ca...
Definition: Algorithm.cpp:1687
const SpectrumInfo & spectrumInfo() const
Return a reference to the SpectrumInfo object.
Geometry::Instrument_const_sptr getInstrument() const
Returns the parameterized instrument.
Base MatrixWorkspace Abstract Class.
virtual std::size_t getNumberHistograms() const =0
Returns the number of histograms in the workspace.
const MantidVec & readY(std::size_t const index) const
Deprecated, use y() instead.
API::SpectrumInfo is an intermediate step towards a SpectrumInfo that is part of Instrument-2....
Definition: SpectrumInfo.h:53
const Geometry::IDetector & detector(const size_t index) const
Return a const reference to the detector or detector group of the spectrum with given index.
A property class for workspaces.
const std::string name() const override
Algorithm's name for identification overriding a virtual method.
int version() const override
Algorithm's version for identification overriding a virtual method.
void init() override
Virtual method - must be overridden by concrete algorithm.
double m_rangeUpper
Ending x-axis value for integrations.
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.
int m_TotalTime
An estimate total number of additions or equilivent required to compute a spectrum.
std::vector< std::vector< size_t > > makeMap(const API::MatrixWorkspace_sptr &countsWS)
method to check which spectra should be grouped when calculating the median
std::vector< std::vector< size_t > > makeInstrumentMap(const API::MatrixWorkspace &countsWS)
method to create the map with all spectra
double m_fracDone
An estimate of the percentage of the algorithm runtimes that has been completed.
RunTime
For the progress bar, estimates of how many additions, or equivalent, member functions will do for ea...
@ RTGetTotalCounts
Estimate of the work required from Integrate for each spectrum.
@ RTGetRate
Work required by the ConvertToDistribution algorithm.
int m_maxIndex
Ending workspace index to run tests on.
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.
API::MatrixWorkspace_sptr convertToRate(API::MatrixWorkspace_sptr workspace)
Convert to a distribution.
double m_progStepWidth
The number of tests to be run.
DataObjects::MaskWorkspace_sptr generateEmptyMask(const API::MatrixWorkspace_const_sptr &inputWS)
Create a masking workspace to return.
void applyMask(const API::MatrixWorkspace_sptr &inputWS, const API::MatrixWorkspace_sptr &maskWS)
Apply a given mask.
void exec() override
Virtual method - must be overridden by concrete algorithm.
API::MatrixWorkspace_sptr doDetVanTest(const API::MatrixWorkspace_sptr &inputWS, int &nFails)
Perform checks on detector vanadium.
double m_rangeLower
Starting x-axis value for integrations.
int m_minIndex
Starting workspace index to run tests on.
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...
const std::string category() const override
Algorithm's category for identification.
base class for Geometric IComponent
Definition: IComponent.h:51
virtual std::vector< std::shared_ptr< const IComponent > > getAncestors() const =0
Return an array of all ancestors, the nearest first.
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void setPropertySettings(const std::string &name, std::unique_ptr< IPropertySettings > settings)
void setPropertyGroup(const std::string &name, const std::string &group)
Set the group for a given property.
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< IAlgorithm > IAlgorithm_sptr
shared pointer to Mantid::API::IAlgorithm
std::shared_ptr< const MatrixWorkspace > MatrixWorkspace_const_sptr
shared pointer to the matrix workspace base class (const version)
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 Mantid::Geometry::IDetector > IDetector_const_sptr
Shared pointer to IDetector (const version)
Definition: IDetector.h:102
std::shared_ptr< const Instrument > Instrument_const_sptr
Shared pointer to an const instrument object.
Statistics getStatistics(const std::vector< TYPE > &data, const unsigned int flags=StatOptions::AllStats)
Return a statistics object for the given data set.
Definition: Statistics.cpp:167
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
Simple struct to store statistics.
Definition: Statistics.h:25
double median
Median value.
Definition: Statistics.h:33