Mantid
Loading...
Searching...
No Matches
FindDetectorsOutsideLimits.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 +
13
14#include <cmath>
15#include <fstream>
16
17namespace Mantid::Algorithms {
18// Register the class into the algorithm factory
19DECLARE_ALGORITHM(FindDetectorsOutsideLimits)
20
21const std::string FindDetectorsOutsideLimits::category() const { return "Diagnostics"; }
22
23using namespace Kernel;
24using namespace API;
25using namespace DataObjects;
27
30 declareProperty(std::make_unique<WorkspaceProperty<MatrixWorkspace>>("InputWorkspace", "", Direction::Input),
31 "Name of the input workspace2D");
32 declareProperty(std::make_unique<WorkspaceProperty<MatrixWorkspace>>("OutputWorkspace", "", Direction::Output),
33 "Each histogram from the input workspace maps to a histogram in this\n"
34 "workspace with one value that indicates if there was a dead detector");
35 declareProperty("HighThreshold", EMPTY_DBL(),
36 "Spectra whose total number of counts are equal to or above this value\n"
37 "will be marked bad (default off)");
38 declareProperty("LowThreshold", 0.0,
39 "Spectra whose total number of counts are equal to or below this value\n"
40 "will be marked bad (default 0)");
41 auto mustBePosInt = std::make_shared<BoundedValidator<int>>();
42 mustBePosInt->setLower(0);
43 declareProperty("StartWorkspaceIndex", 0, mustBePosInt,
44 "The index number of the first spectrum to include in the calculation\n"
45 "(default 0)");
46 declareProperty("EndWorkspaceIndex", EMPTY_INT(), mustBePosInt,
47 "The index number of the last spectrum to include in the calculation\n"
48 "(default the last histogram)");
49 declareProperty("RangeLower", EMPTY_DBL(),
50 "No bin with a boundary at an x value less than this will be used\n"
51 "in the summation that decides if a detector is 'bad' (default: the\n"
52 "start of each histogram)");
53 declareProperty("RangeUpper", EMPTY_DBL(),
54 "No bin with a boundary at an x value higher than this value will\n"
55 "be used in the summation that decides if a detector is 'bad'\n"
56 "(default: the end of each histogram)");
57 declareProperty("NumberOfFailures", 0, Direction::Output);
58}
59
67 double lowThreshold = getProperty("LowThreshold");
68 double highThreshold = getProperty("HighThreshold");
69 bool useHighThreshold = !isEmpty(highThreshold);
70 if (useHighThreshold && highThreshold <= lowThreshold) {
71 throw std::invalid_argument("The high threshold must be higher than the low threshold");
72 }
73
74 MatrixWorkspace_sptr inputWS = getProperty("InputWorkspace");
75 int minIndex = getProperty("StartWorkspaceIndex");
76 int maxIndex = getProperty("EndWorkspaceIndex");
77 const auto inputLength = static_cast<int>(inputWS->getNumberHistograms());
78 if (isEmpty(maxIndex))
79 maxIndex = inputLength - 1;
80 if (maxIndex < minIndex) {
81 std::ostringstream os;
82 os << "Invalid workspace indices. Min=" << minIndex << ",Max=" << maxIndex;
83 throw std::invalid_argument(os.str());
84 }
85
86 if (minIndex > inputLength || maxIndex > inputLength) {
87 std::ostringstream os;
88 os << "Workspace indices must be lower than workspace size. Size=" << inputLength << ", Min=" << minIndex
89 << ",Max=" << maxIndex;
90 throw std::invalid_argument(os.str());
91 }
92
93 const double rangeLower = getProperty("RangeLower");
94 const double rangeUpper = getProperty("RangeUpper");
95 // Get the integrated input workspace; converting to a Workspace2D
96 MatrixWorkspace_sptr countsWS = integrateSpectra(inputWS, minIndex, maxIndex, rangeLower, rangeUpper, true);
97 if (std::dynamic_pointer_cast<EventWorkspace>(countsWS))
98 throw std::runtime_error("Error! Integration output is not a Workspace2D.");
99
100 // Create a new workspace for the results, copy from the input to ensure that
101 // we copy over the instrument and current masking
102 MaskWorkspace_sptr outputWS = this->generateEmptyMask(inputWS);
103
104 const double deadValue(1.0); // delete the data
105
106 const int diagLength = maxIndex - minIndex;
107 const auto progStep = static_cast<int>(std::ceil(diagLength / 100.0));
108
109 bool checkForMask = false;
110 Geometry::Instrument_const_sptr instrument = inputWS->getInstrument();
111 if (instrument != nullptr) {
112 checkForMask = ((instrument->getSource() != nullptr) && (instrument->getSample() != nullptr));
113 }
114
115 int numFailed(0);
116 const auto &spectrumInfo = countsWS->spectrumInfo();
117#pragma omp parallel for if (countsWS->threadSafe() && outputWS->threadSafe()), reduction(+ : numFailed)
118 for (int i = minIndex; i <= maxIndex; ++i) {
119 bool keepData = true;
120 int countsInd = i - minIndex;
121 if (i % progStep == 0) {
122 progress(static_cast<double>(countsInd) / diagLength);
124 }
125
126 if (checkForMask && spectrumInfo.hasDetectors(countsInd)) {
127 if (spectrumInfo.isMonitor(countsInd))
128 continue; // do include or exclude from mask
129 if (spectrumInfo.isMasked(countsInd))
130 keepData = false;
131 }
132
133 const double &yValue = countsWS->y(countsInd)[0];
134 // Mask out NaN and infinite
135 if (!std::isfinite(yValue)) {
136 keepData = false;
137 } else {
138 if (yValue <= lowThreshold) {
139 keepData = false;
140 }
141
142 if (useHighThreshold && (yValue >= highThreshold)) {
143 keepData = false;
144 }
145 }
146
147 if (!keepData) {
148 outputWS->mutableY(i)[0] = deadValue;
149 ++numFailed;
150 }
151 }
152
153 g_log.information() << numFailed << " spectra fell outside the given limits.\n";
154 setProperty("NumberOfFailures", numFailed);
155 // Assign it to the output workspace property
156 setProperty("OutputWorkspace", outputWS);
157}
158} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
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
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
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
static bool isEmpty(const NumT toCheck)
checks that the value was not set by users, uses the value in empty double/int.
A property class for workspaces.
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.
Takes a workspace as input and identifies all spectra where the sum of the counts in all bins is outs...
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
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 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.
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
STL namespace.
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54