19#include <boost/iterator/counting_iterator.hpp>
47 "Name of the integrated detector vanadium (white beam) workspace.");
50 "A hard mask to apply to the inputworkspace");
52 "A MaskWorkspace containing the masked spectra as zeroes and ones.");
53 auto mustBePosInt = std::make_shared<BoundedValidator<int>>();
54 mustBePosInt->setLower(0);
56 "The index number of the first spectrum to include in the calculation\n"
59 "The index number of the last spectrum to include in the calculation\n"
60 "(default the last histogram)");
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)");
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)");
70 string findDetOutLimGrp(
"Find Detectors Outside Limits");
72 "Spectra whose total number of counts are equal to or below this value\n"
73 "will be marked bad (default 0)");
76 "Spectra whose total number of counts are equal to or above this value\n"
77 "will be marked bad (default off)");
80 string medianDetTestGrp(
"Median Detector Test");
81 auto mustBePositiveDbl = std::make_shared<BoundedValidator<double>>();
82 mustBePositiveDbl->setLower(0);
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 "
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");
94 this->
declareProperty(
"LowThresholdFraction", 0.1,
"Lower acceptable bound as fraction of median value");
96 this->
declareProperty(
"HighThresholdFraction", 1.5,
"Upper acceptable bound as fraction of median value");
98 this->
declareProperty(
"LowOutlier", 0.01,
"Lower bound defining outliers as fraction of median value");
100 this->
declareProperty(
"HighOutlier", 100.,
"Upper bound defining outliers as fraction of median value");
102 this->
declareProperty(
"CorrectForSolidAngle",
false,
"Flag to correct for solid angle efficiency. False by default.");
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");
110 string detEffVarGrp(
"Detector Efficiency Variation");
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 "
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");
122 std::make_unique<EnabledWhenProperty>(
"DetVanCompare",
IS_NOT_DEFAULT));
124 string countsCheck(
"Check Sample Counts");
127 "A sample workspace integrated over the full axis range.");
130 string backgroundCheck(
"Check Sample Background");
133 "A sample workspace integrated over the background region.");
135 this->
declareProperty(
"SampleBkgLowAcceptanceFactor", 0.0, mustBePositiveDbl,
136 "Low threshold for the background check MedianDetectorTest.");
138 this->
declareProperty(
"SampleBkgHighAcceptanceFactor", 5.0, mustBePositiveDbl,
139 "High threshold for the background check MedianDetectorTest.");
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");
148 "Flag to correct for solid angle efficiency for "
149 "background check MedianDetectorTest. False by "
153 string psdBleedMaskGrp(
"Create PSD Bleed Mask");
156 "A sample workspace. This is used in the PSD Bleed calculation.");
159 "The maximum rate allowed for a tube in counts/us/frame.");
161 this->
declareProperty(
"NIgnoredCentralPixels", 80, mustBePosInt,
"The number of pixels about the centre to ignore.");
164 std::make_unique<EnabledWhenProperty>(
"MaxTubeFramerate",
IS_NOT_DEFAULT));
209 md->setProperty(
"Workspace", inputWS);
210 md->setProperty(
"MaskedWorkspace", hardMaskWS);
211 md->executeAsChildAlg();
226 double variation = this->
getProperty(
"DetVanRatioVariation");
231 alg->setProperty(
"WhiteBeamBase", inputWS);
232 alg->setProperty(
"WhiteBeamCompare", input2WS);
233 alg->setProperty(
"StartWorkspaceIndex",
m_minIndex);
234 alg->setProperty(
"EndWorkspaceIndex",
m_maxIndex);
237 alg->setProperty(
"Variation", variation);
238 alg->executeAsChildAlg();
242 int localFails = alg->getProperty(
"NumberOfFailures");
243 numFailed += localFails;
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();
261 int localFails = zeroChk->getProperty(
"NumberOfFailures");
262 numFailed += localFails;
270 double significanceTest = this->
getProperty(
"SampleBkgSignificanceTest");
271 double lowThreshold = this->
getProperty(
"SampleBkgLowAcceptanceFactor");
272 double highThreshold = this->
getProperty(
"SampleBkgHighAcceptanceFactor");
273 bool correctSA = this->
getProperty(
"SampleCorrectForSolidAngle");
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();
291 int localFails = alg->getProperty(
"NumberOfFailures");
292 numFailed += localFails;
298 double maxTubeFrameRate = this->
getProperty(
"MaxTubeFramerate");
299 int numIgnore = this->
getProperty(
"NIgnoredCentralPixels");
304 alg->setProperty(
"InputWorkspace", sampleWS);
305 alg->setProperty(
"MaxTubeFramerate", maxTubeFrameRate);
306 alg->setProperty(
"NIgnoredCentralPixels", numIgnore);
307 alg->executeAsChildAlg();
310 int localFails = alg->getProperty(
"NumberOfFailures");
311 numFailed += localFails;
318 std::vector<int> detList;
320 extract->setProperty(
"InputWorkspace", inputWS);
321 extract->setProperty(
"OutputWorkspace",
"final_mask");
322 extract->setProperty(
"DetectorList", detList);
323 extract->executeAsChildAlg();
324 maskWS = extract->getProperty(
"OutputWorkspace");
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();
354 double lowThreshold = this->
getProperty(
"LowThreshold");
355 double highThreshold = this->
getProperty(
"HighThreshold");
359 fdol->setProperty(
"InputWorkspace", inputWS);
360 fdol->setProperty(
"OutputWorkspace", localMask);
361 fdol->setProperty(
"StartWorkspaceIndex",
m_minIndex);
362 fdol->setProperty(
"EndWorkspaceIndex",
m_maxIndex);
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;
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");
389 mdt->setProperty(
"InputWorkspace", inputWS);
390 mdt->setProperty(
"StartWorkspaceIndex",
m_minIndex);
391 mdt->setProperty(
"EndWorkspaceIndex",
m_maxIndex);
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;
415 :
API::
Algorithm(), m_fracDone(0.0), m_TotalTime(RTTotal), m_parents(0), m_progStepWidth(0.0), m_minIndex(0),
434 const int indexMax,
const double lower,
const double upper,
435 const bool outputWorkspace2D) {
436 g_log.
debug() <<
"Integrating input spectra.\n";
445 childAlg->setProperty(
"InputWorkspace", inputWS);
446 childAlg->setProperty(
"StartWorkspaceIndex", indexMin);
447 childAlg->setProperty(
"EndWorkspaceIndex", indexMax);
450 childAlg->setProperty(
"RangeLower",
lower);
451 childAlg->setProperty(
"RangeUpper",
upper);
452 childAlg->setPropertyValue(
"IncludePartialBins",
"1");
453 childAlg->executeAsChildAlg();
458 if (outputWorkspace2D && std::dynamic_pointer_cast<EventWorkspace>(outputW)) {
459 g_log.
debug() <<
"Converting output Event Workspace into a Workspace2D.\n";
461 childAlg->setProperty(
"InputWorkspace", outputW);
462 childAlg->executeAsChildAlg();
463 finalOutputW = childAlg->getProperty(
"OutputWorkspace");
478 auto maskWS = create<DataObjects::MaskWorkspace>(*inputWS, HistogramData::Points(1));
479 maskWS->setTitle(inputWS->getTitle());
485 return {{boost::counting_iterator<std::size_t>(0),
493 std::multimap<Mantid::Geometry::ComponentID, size_t> mymap;
500 g_log.
warning(
"Workspace has no instrument. LevelsUP is ignored");
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.");
510 const SpectrumInfo &spectrumInfo = countsWS->spectrumInfo();
512 for (
size_t i = 0; i < countsWS->getNumberHistograms(); ++i) {
515 if (anc.size() <
static_cast<size_t>(
m_parents)) {
516 g_log.
warning(
"Too many levels up. Will ignore LevelsUp");
520 mymap.emplace(anc[
m_parents - 1]->getComponentID(), i);
523 std::vector<std::vector<size_t>> speclist;
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) {
528 auto keyRange = mymap.equal_range(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);
534 speclist.emplace_back(std::move(speclistsingle));
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");
558 bool checkForMask =
false;
560 if (instrument !=
nullptr) {
561 checkForMask = ((instrument->getSource() !=
nullptr) && (instrument->getSample() !=
nullptr));
565 for (
const auto &hists : indexmap) {
566 std::vector<double> medianInput;
567 const auto nhists =
static_cast<int>(hists.size());
569 medianInput.reserve(nhists);
572 for (
int i = 0; i < nhists; ++i) {
575 if (checkForMask && spectrumInfo.hasDetectors(hists[i])) {
576 if (spectrumInfo.isMasked(hists[i]) || spectrumInfo.isMonitor(hists[i]))
580 const double yValue = input.
readY(hists[i])[0];
582 throw std::out_of_range(
"Negative number of counts found, could be "
583 "corrupted raw counts or solid angle data");
585 if (!std::isfinite(yValue) || (excludeZeroes && yValue < DBL_EPSILON))
590 PARALLEL_CRITICAL(DetectorDiagnostic_median_d) { medianInput.emplace_back(yValue); }
596 if (medianInput.empty()) {
597 g_log.
information(
"some group has no valid histograms. Will use 0 for median.");
598 medianInput.emplace_back(0.);
602 double median = stats.
median;
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");
608 medianvec.emplace_back(median);
620 g_log.
information() <<
"Workspace already contains a count rate, nothing to do.\n";
631 return childAlg->getProperty(
"Workspace");
#define DECLARE_ALGORITHM(classname)
IPeaksWorkspace_sptr workspace
#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...
#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
Base class from which all concrete algorithm classes should be derived.
void declareProperty(std::unique_ptr< Kernel::Property > p, const std::string &doc="") override
Add a property to the list of managed properties.
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
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.
void interruption_point()
This is called during long-running operations, and check if the algorithm has requested that it be ca...
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....
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.
DetectorDiagnostic()
Default constructor.
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
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.
void warning(const std::string &msg)
Logs at warning level.
void information(const std::string &msg)
Logs at information level.
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
std::shared_ptr< const Mantid::Geometry::IDetector > IDetector_const_sptr
Shared pointer to IDetector (const version)
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.
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.
constexpr int EMPTY_INT() noexcept
Returns what we consider an "empty" integer within a property.
constexpr double EMPTY_DBL() noexcept
Returns what we consider an "empty" double within a property.
@ Input
An input workspace.
@ Output
An output workspace.
Simple struct to store statistics.
double median
Median value.