Mantid
Loading...
Searching...
No Matches
DetectorEfficiencyVariation.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
12namespace Mantid::Algorithms {
13
14// Register the class into the algorithm factory
15DECLARE_ALGORITHM(DetectorEfficiencyVariation)
16
17const std::string DetectorEfficiencyVariation::category() const { return "Diagnostics"; }
18
19using namespace Kernel;
20using namespace API;
23
26 auto val = std::make_shared<HistogramValidator>();
27 declareProperty(std::make_unique<WorkspaceProperty<MatrixWorkspace>>("WhiteBeamBase", "", Direction::Input, val),
28 "Name of a white beam vanadium workspace");
29 // The histograms, the detectors in each histogram and their first and last
30 // bin boundary must match
31 declareProperty(std::make_unique<WorkspaceProperty<MatrixWorkspace>>("WhiteBeamCompare", "", Direction::Input, val),
32 "Name of a matching second white beam vanadium run from the same "
33 "instrument");
34 declareProperty(std::make_unique<WorkspaceProperty<MatrixWorkspace>>("OutputWorkspace", "", Direction::Output),
35 "A MaskWorkpace where each spectra that failed the test is "
36 "masked. Each histogram from the input workspace maps to a "
37 "histogram in this workspace with one value that indicates "
38 "if there was a dead detector.");
39 auto moreThanZero = std::make_shared<BoundedValidator<double>>();
40 moreThanZero->setLower(0.0);
41 declareProperty("Variation", 1.1, moreThanZero,
42 "Identify histograms whose total number of counts has "
43 "changed by more than this factor of the median change "
44 "between the two input workspaces.");
45 auto mustBePosInt = std::make_shared<BoundedValidator<int>>();
46 mustBePosInt->setLower(0);
47 declareProperty("StartWorkspaceIndex", 0, mustBePosInt,
48 "The index number of the first spectrum to include in the "
49 "calculation (default: 0)");
50
51 // Mantid::EMPTY_INT() and EMPTY_DBL() are tags that indicate that no
52 // value has been set and we want to use the default
53 declareProperty("EndWorkspaceIndex", Mantid::EMPTY_INT(), mustBePosInt,
54 "The index number of the last spectrum to include in the "
55 "calculation (default: the last spectrum in the workspace)");
56 declareProperty("RangeLower", Mantid::EMPTY_DBL(),
57 "No bin with a boundary at an x value less than this will be included "
58 "in the summation used to decide if a detector is 'bad' (default: the "
59 "start of each histogram)");
60 declareProperty("RangeUpper", Mantid::EMPTY_DBL(),
61 "No bin with a boundary at an x value higher than this value will "
62 "be included in the summation used to decide if a detector is 'bad' "
63 "(default: the end of each histogram)");
64 declareProperty("NumberOfFailures", 0, Direction::Output);
65}
66
76 double variation = Mantid::EMPTY_DBL();
77 int minSpec = 0;
78 int maxSpec = Mantid::EMPTY_INT();
79 retrieveProperties(WB1, WB2, variation, minSpec, maxSpec);
80
81 const double rangeLower = getProperty("RangeLower");
82 const double rangeUpper = getProperty("RangeUpper");
83
84 MatrixWorkspace_sptr counts1 = integrateSpectra(WB1, minSpec, maxSpec, rangeLower, rangeUpper);
85 MatrixWorkspace_sptr counts2 = integrateSpectra(WB2, minSpec, maxSpec, rangeLower, rangeUpper);
86 MatrixWorkspace_sptr countRatio;
87 try {
88 // Note. This can produce NAN/INFs. Leave for now and sort it out in the
89 // later tests
90 countRatio = counts1 / counts2;
91 } catch (std::invalid_argument &) {
92 g_log.error() << "The two white beam workspaces size must match.";
93 throw;
94 }
95 double average = calculateMedian(*countRatio, false, makeInstrumentMap(*countRatio)).at(0); // Include zeroes
96 g_log.notice() << name() << ": The median of the ratio of the integrated counts is: " << average << '\n';
97 //
98 int numFailed = doDetectorTests(counts1, counts2, average, variation);
99
100 g_log.notice() << "Tests failed " << numFailed << " spectra. "
101 << "These have been masked on the OutputWorkspace\n";
102
103 // counts1 was overwriten by the last function, now register it
104 setProperty("NumberOfFailures", numFailed);
105}
106
120 API::MatrixWorkspace_sptr &whiteBeam2, double &variation,
121 int &startWsIndex, int &endWsIndex) {
122 whiteBeam1 = getProperty("WhiteBeamBase");
123 whiteBeam2 = getProperty("WhiteBeamCompare");
124 if (whiteBeam1->getInstrument()->getName() != whiteBeam2->getInstrument()->getName()) {
125 throw std::invalid_argument("The two input white beam vanadium workspaces "
126 "must be from the same instrument");
127 }
128 int maxWsIndex = static_cast<int>(whiteBeam1->getNumberHistograms()) - 1;
129 if (maxWsIndex !=
130 static_cast<int>(whiteBeam2->getNumberHistograms()) - 1) { // we would get a crash later on if this were not true
131 throw std::invalid_argument("The input white beam vanadium workspaces must "
132 "be have the same number of histograms");
133 }
134
135 variation = getProperty("Variation");
136
137 startWsIndex = getProperty("StartWorkspaceIndex");
138 if ((startWsIndex < 0) || (startWsIndex > maxWsIndex)) {
139 g_log.warning("StartWorkspaceIndex out of range, changed to 0");
140 startWsIndex = 0;
141 }
142 endWsIndex = getProperty("EndWorkspaceIndex");
143 if (endWsIndex == Mantid::EMPTY_INT())
144 endWsIndex = maxWsIndex;
145 if ((endWsIndex < 0) || (endWsIndex > maxWsIndex)) {
146 g_log.warning("EndWorkspaceIndex out of range, changed to max Workspace number");
147 endWsIndex = maxWsIndex;
148 }
149 if ((endWsIndex < startWsIndex)) {
150 g_log.warning("EndWorkspaceIndex can not be less than the StartWorkspaceIndex, "
151 "changed to max Workspace number");
152 endWsIndex = maxWsIndex;
153 }
154}
155
170 const API::MatrixWorkspace_const_sptr &counts2, const double average,
171 double variation) {
172 // DIAG in libISIS did this. A variation of less than 1 doesn't make sense in
173 // this algorithm
174 if (variation < 1) {
175 variation = 1.0 / variation;
176 }
177 // criterion for if the the first spectrum is larger than expected
178 double largest = average * variation;
179 // criterion for if the the first spectrum is lower than expected
180 double lowest = average / variation;
181
182 const auto numSpec = static_cast<int>(counts1->getNumberHistograms());
183 const auto progStep = static_cast<int>(std::ceil(numSpec / 30.0));
184
185 // Create a workspace for the output
186 MaskWorkspace_sptr maskWS = this->generateEmptyMask(counts1);
187
188 bool checkForMask = false;
189 Geometry::Instrument_const_sptr instrument = counts1->getInstrument();
190 if (instrument != nullptr) {
191 checkForMask = ((instrument->getSource() != nullptr) && (instrument->getSample() != nullptr));
192 }
193
194 const double deadValue(1.0);
195 int numFailed(0);
196 const auto &spectrumInfo = counts1->spectrumInfo();
197 PARALLEL_FOR_IF(Kernel::threadSafe(*counts1, *counts2, *maskWS))
198 for (int i = 0; i < numSpec; ++i) {
200 // move progress bar
201 if (i % progStep == 0) {
202 advanceProgress(progStep * static_cast<double>(RTMarkDetects) / numSpec);
205 }
206
207 if (checkForMask) {
208 if (spectrumInfo.isMonitor(i))
209 continue;
210 if (spectrumInfo.isMasked(i)) {
211 // Ensure it is masked on the output
212 maskWS->mutableY(i)[0] = deadValue;
213 continue;
214 }
215 }
216
217 const double signal1 = counts1->y(i)[0];
218 const double signal2 = counts2->y(i)[0];
219
220 // Mask out NaN and infinite
221 if (!std::isfinite(signal1) || !std::isfinite(signal2)) {
222 maskWS->mutableY(i)[0] = deadValue;
224 ++numFailed;
225 continue;
226 }
227
228 // Check the ratio is within the given range
229 const double ratio = signal1 / signal2;
230 if (ratio < lowest || ratio > largest) {
231 maskWS->mutableY(i)[0] = deadValue;
233 ++numFailed;
234 }
235
237 }
239
240 // Register the results with the ADS
241 setProperty("OutputWorkspace", maskWS);
242
243 return numFailed;
244}
245
246} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
#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_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
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
A property class for workspaces.
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 > > 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.
@ 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...
const std::string name() const override
Algorithm's name for identification overriding a virtual method.
void retrieveProperties(API::MatrixWorkspace_sptr &whiteBeam1, API::MatrixWorkspace_sptr &whiteBeam2, double &variation, int &minSpec, int &maxSpec)
Loads and checks the values passed to the algorithm.
int doDetectorTests(const API::MatrixWorkspace_const_sptr &counts1, const API::MatrixWorkspace_const_sptr &counts2, const double average, double variation)
Apply the detector test criterion.
void exec() override
Executes the algorithm that includes calls to SolidAngle and Integration.
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void notice(const std::string &msg)
Logs at notice level.
Definition: Logger.cpp:95
void error(const std::string &msg)
Logs at error level.
Definition: Logger.cpp:77
void warning(const std::string &msg)
Logs at warning level.
Definition: Logger.cpp:86
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.
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
STL namespace.
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54