Mantid
Loading...
Searching...
No Matches
CalculateEfficiency2.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 +
12#include <cmath>
13#include <limits>
14#include <vector>
15
16namespace Mantid::Algorithms {
17
18// Register the class into the algorithm factory
19DECLARE_ALGORITHM(CalculateEfficiency2)
20
21using namespace Kernel;
22using namespace API;
23using namespace Geometry;
24using namespace DataObjects;
25
27namespace PropertyNames {
28const static std::string INPUT_WORKSPACE{"InputWorkspace"};
29const static std::string OUTPUT_WORKSPACE{"OutputWorkspace"};
30const static std::string MIN_THRESHOLD{"MinThreshold"};
31const static std::string MAX_THRESHOLD{"MaxThreshold"};
32const static std::string MERGE_GROUP{"MergeGroup"};
33} // namespace PropertyNames
34
35namespace { // anonymous
36static void applyBadPixelThreshold(MatrixWorkspace &outputWS, double minThreshold, double maxThreshold) {
37
38 // Number of spectra
39 const size_t numberOfSpectra = outputWS.getNumberHistograms();
40 const auto &spectrumInfo = outputWS.spectrumInfo();
41
42 for (size_t i = 0; i < numberOfSpectra; i++) {
43 auto &YOut = outputWS.mutableY(i);
44 auto &EOut = outputWS.mutableE(i);
45
46 // Skip if we have a monitor or if the detector is masked.
47 if (spectrumInfo.isMonitor(i)) {
48 YOut.front() = 1.0;
49 EOut.front() = 0.0;
50 continue;
51 } else if (spectrumInfo.isMasked(i)) {
52 continue;
53 }
54
55 // if the pixel is outside the thresholds let make it EMPTY_DBL
56 // In the documentation is "-inf"
57 const auto y = YOut.front();
58 if (y < minThreshold || y > maxThreshold) {
59 YOut.front() = EMPTY_DBL();
60 EOut.front() = EMPTY_DBL();
61 }
62 }
63}
64} // anonymous namespace
65
71 "The workspace containing the flood data");
74 "The name of the workspace to be created as the output of the algorithm");
75
76 auto positiveDouble = std::make_shared<BoundedValidator<double>>();
77 positiveDouble->setLower(0);
78 declareProperty(PropertyNames::MIN_THRESHOLD, 0.0, positiveDouble, "Minimum threshold for a pixel to be considered");
79 declareProperty(PropertyNames::MAX_THRESHOLD, 2.0, positiveDouble->clone(),
80 "Maximum threshold for a pixel to be considered");
81
83 "Whether to merge entries when WorkspaceGroup is specified as input.");
84}
85
86std::map<std::string, std::string> CalculateEfficiency2::validateInputs() {
87 std::map<std::string, std::string> result;
88
89 // Files from time-of-flight instruments must be integrated in Lambda before
90 // using this algorithm
91
92 auto oneBinMsg = "Input workspace must have only one bin. Consider "
93 "integrating the input over all the bins.";
95 MatrixWorkspace_const_sptr inputWS = std::dynamic_pointer_cast<const MatrixWorkspace>(ws1);
96 if (inputWS == nullptr) {
97 WorkspaceGroup_const_sptr inputGroup = std::dynamic_pointer_cast<const WorkspaceGroup>(ws1);
98 if (inputGroup != nullptr) {
99 for (auto entryNo = 0; entryNo < inputGroup->getNumberOfEntries(); ++entryNo) {
100 auto const entry = std::static_pointer_cast<const MatrixWorkspace>(inputGroup->getItem(entryNo));
101 if (entry->blocksize() > 1) {
102 result[PropertyNames::INPUT_WORKSPACE] = oneBinMsg;
103 break;
104 }
105 }
106 } else {
107 result[PropertyNames::INPUT_WORKSPACE] = "The input property must be either MatrixWorkspace or a "
108 "WorkspaceGroup containing MatrixWorkspaces";
109 }
110 } else if (inputWS->blocksize() > 1) {
111 result[PropertyNames::INPUT_WORKSPACE] = oneBinMsg;
112 }
113
115 result[PropertyNames::OUTPUT_WORKSPACE] = "The output workspace name must be specified.";
116 }
117
118 // get the thresholds once to error check and use in the main function
119 m_minThreshold = getProperty("MinThreshold");
120 m_maxThreshold = getProperty("MaxThreshold");
122 const std::string msg{"MinThreshold must be less than MaxThreshold"};
123 result[PropertyNames::MIN_THRESHOLD] = msg;
124 result[PropertyNames::MAX_THRESHOLD] = msg;
125 }
126
127 return result;
128}
129
135 auto inputWS = std::static_pointer_cast<MatrixWorkspace>(ws1);
136 auto outputWS = calculateEfficiency(inputWS);
138 progress(1.0, "Done!");
139}
140
142 double startProgress, double stepProgress) {
143
144 MatrixWorkspace_sptr outputWS;
145 if (std::dynamic_pointer_cast<EventWorkspace>(inputWorkspace)) {
146 // if event workspace, create the output workspace from the input, while NOT preserving events
147 auto childAlg = createChildAlgorithm("RebinToWorkspace", 0.0, 0.1);
148 childAlg->setProperty("WorkspaceToRebin", inputWorkspace);
149 childAlg->setProperty("WorkspaceToMatch", inputWorkspace);
150 childAlg->setPropertyValue("OutputWorkspace", getPropertyValue(PropertyNames::OUTPUT_WORKSPACE));
151 childAlg->setProperty("PreserveEvents", false);
152 childAlg->executeAsChildAlg();
153 outputWS = childAlg->getProperty("OutputWorkspace");
154 } else {
155 // otherwise just clone the input
156 outputWS = inputWorkspace->clone();
157 }
158
159 // Loop over spectra and sum all the counts to get normalization
160 // Skip monitors and masked detectors
161 // returns tuple with (sum, err, npixels)
162 progress(startProgress + 0.1 * stepProgress, "Computing the counts.");
163 auto counts = sumUnmaskedAndDeadPixels(*outputWS);
164 if (counts.nPixels == 0) {
165 throw std::runtime_error("No pixels being used for calculation");
166 }
167
168 progress(startProgress + 0.3 * stepProgress, "Normalising the detectors.");
169 averageAndNormalizePixels(*outputWS, counts);
170
171 progress(startProgress + 0.5 * stepProgress, "Applying bad pixel threshold.");
172 applyBadPixelThreshold(*outputWS, m_minThreshold, m_maxThreshold);
173
174 // do it again only using the pixels that are within the threshold
175 progress(startProgress + 0.7 * stepProgress, "Computing the counts.");
176 counts = sumUnmaskedAndDeadPixels(*outputWS);
177 if (counts.nPixels == 0) {
178 throw std::runtime_error("All pixels are outside of the threshold values");
179 }
180
181 progress(startProgress + 0.9 * stepProgress, "Normalising the detectors.");
182 averageAndNormalizePixels(*outputWS, counts);
183
184 return outputWS;
185}
186
193 auto results = validateInputs();
194 for (const auto &result : results) {
195 throw std::runtime_error("Issue in " + result.first + " property: " + result.second);
196 }
197}
198
205
206 // if run as a script, validateInputs will not be triggered
207 // for the processGroups, so properties will be validated manually
209
211 WorkspaceGroup_sptr inputWS = std::static_pointer_cast<WorkspaceGroup>(ws1);
212
213 const bool mergeGroups = getProperty(PropertyNames::MERGE_GROUP);
214 if (mergeGroups) {
215 auto mergedWS = mergeGroup(*inputWS);
216 auto outputWS = calculateEfficiency(mergedWS);
218 } else {
219 auto outputGroup = std::make_shared<WorkspaceGroup>();
220 auto const nEntries = inputWS->getNumberOfEntries();
221 auto const stepProgress = 1.0 / nEntries;
222 for (auto entryNo = 0; entryNo < nEntries; ++entryNo) {
223 auto entryWS = std::static_pointer_cast<API::MatrixWorkspace>(inputWS->getItem(entryNo));
224 auto const startProgress = static_cast<double>(entryNo) / nEntries;
225 auto outputWS = calculateEfficiency(entryWS, startProgress, stepProgress);
226 outputGroup->addWorkspace(outputWS);
227 }
228 const std::string groupName = getPropertyValue(PropertyNames::OUTPUT_WORKSPACE);
229 AnalysisDataService::Instance().addOrReplace(groupName, outputGroup);
231 }
232 progress(1.0, "Done!");
233 return true;
234}
235
242 // Number of spectra
243 const size_t numberOfSpectra = workspace.getNumberHistograms();
244 SummedResults results;
245
246 const auto &spectrumInfo = workspace.spectrumInfo();
247 for (size_t i = 0; i < numberOfSpectra; i++) {
248
249 // Retrieve the spectrum into a vector
250 auto &YValues = workspace.y(i);
251 auto &YErrors = workspace.e(i);
252
253 // Skip if we have a monitor, if the detector is masked or if the pixel is
254 // dead
255 if (spectrumInfo.isMonitor(i) || spectrumInfo.isMasked(i) || isEmpty(YValues.front()))
256 continue;
257
258 results.sum += YValues.front();
259 results.error += YErrors.front() * YErrors.front();
260 results.nPixels++;
261 }
262 results.error = std::sqrt(results.error);
263
264 g_log.debug() << "Total of unmasked/dead pixels = " << results.nPixels << " from a total of " << numberOfSpectra
265 << "\n";
266
267 return results;
268}
269
271
272 // Number of spectra
273 const size_t numberOfSpectra = workspace.getNumberHistograms();
274 const auto &spectrumInfo = workspace.spectrumInfo();
275 // Calculate the averages
276 const double averageY = counts.sum / static_cast<double>(counts.nPixels);
277 const double averageE = counts.error / static_cast<double>(counts.nPixels);
278
279 for (size_t i = 0; i < numberOfSpectra; i++) {
280
281 auto &y = workspace.mutableY(i);
282 auto &e = workspace.mutableE(i);
283 const auto yOriginal = y.front();
284
285 // Skip if we have a monitor, the detector is masked, or it has already been
286 // marked as outside of the threashold by being set to EMPTY_DBL
287 if (spectrumInfo.isMasked(i) || spectrumInfo.isMonitor(i) || isEmpty(yOriginal))
288 continue;
289
290 const auto eOriginal = e.front();
291
292 // Normalize counts
293 y.front() = yOriginal / averageY;
294 const double signalToNoiseOrig = eOriginal / yOriginal;
295 const double signalToNoiseAvg = averageE / averageY;
296 e.front() = y.front() * std::sqrt((signalToNoiseOrig * signalToNoiseOrig) + (signalToNoiseAvg * signalToNoiseAvg));
297 }
298
299 g_log.debug() << "Averages :: counts = " << averageY << "; error = " << averageE << "\n";
300}
301
312 auto nEntries = input.getNumberOfEntries();
313 auto mergeRuns = createChildAlgorithm("MergeRuns");
314 mergeRuns->setProperty("InputWorkspaces", input.getName());
315 mergeRuns->executeAsChildAlg();
316 Workspace_sptr mergedWs = mergeRuns->getProperty("OutputWorkspace");
317
318 auto createSingleAlg = createChildAlgorithm("CreateSingleValuedWorkspace");
319 createSingleAlg->setProperty("DataValue", static_cast<double>(nEntries));
320 createSingleAlg->executeAsChildAlg();
321 MatrixWorkspace_sptr normWs = createSingleAlg->getProperty("OutputWorkspace");
322
323 auto divideAlg = createChildAlgorithm("Divide");
324 divideAlg->setProperty("LHSWorkspace", mergedWs);
325 divideAlg->setProperty("RHSWorkspace", normWs);
326 divideAlg->executeAsChildAlg();
327 MatrixWorkspace_sptr mergedNormalisedWs = divideAlg->getProperty("OutputWorkspace");
328
329 auto spectrumInfo = mergedNormalisedWs->spectrumInfo();
330 PARALLEL_FOR_IF(Kernel::threadSafe(*mergedNormalisedWs))
331 for (auto spectrumNo = 0; spectrumNo < static_cast<int>(mergedNormalisedWs->getNumberHistograms()); ++spectrumNo) {
332 if (spectrumInfo.isMasked(spectrumNo)) {
333 auto &detDataY = mergedNormalisedWs->mutableY(spectrumNo);
334 auto &detDataErr = mergedNormalisedWs->mutableE(spectrumNo);
335 auto dataY = 0.0;
336 auto dataE = 0.0;
337 auto nonMaskedEntries = 0;
338 for (auto entryNo = 0; entryNo < nEntries; ++entryNo) {
339 MatrixWorkspace_sptr entry = std::static_pointer_cast<API::MatrixWorkspace>(input.getItem(entryNo));
340 auto spectrumInfoEntry = entry->spectrumInfo();
341 if (!spectrumInfoEntry.isMasked(spectrumNo)) {
342 dataY += entry->readY(spectrumNo)[0];
343 dataE += pow(entry->readE(spectrumNo)[0], 2); // propagate errors
344 nonMaskedEntries++;
345 }
346 }
347 if (nonMaskedEntries != 0) {
348
349 spectrumInfo.setMasked(spectrumNo, false);
350 detDataY.front() = dataY / nonMaskedEntries;
351 detDataErr.front() = sqrt(dataE) / nonMaskedEntries;
352 }
353 }
354 }
355 return mergedNormalisedWs;
356}
357
358} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
IPeaksWorkspace_sptr workspace
Definition: IndexPeaks.cpp:114
#define PARALLEL_FOR_IF(condition)
Empty definitions - to enable set your complier to enable openMP.
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
std::string getPropertyValue(const std::string &name) const override
Get the value of a property as a string.
Definition: Algorithm.cpp:2026
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
static bool isEmpty(const NumT toCheck)
checks that the value was not set by users, uses the value in empty double/int.
const SpectrumInfo & spectrumInfo() const
Return a reference to the SpectrumInfo object.
Base MatrixWorkspace Abstract Class.
virtual std::size_t getNumberHistograms() const =0
Returns the number of histograms in the workspace.
HistogramData::HistogramE & mutableE(const size_t index) &
HistogramData::HistogramY & mutableY(const size_t index) &
Class to hold a set of workspaces.
int getNumberOfEntries() const
Return the number of entries within the group.
Workspace_sptr getItem(const size_t index) const
Return the ith workspace.
A property class for workspaces.
const std::string & getName() const override
Get the workspace name.
Definition: Workspace.cpp:58
void averageAndNormalizePixels(API::MatrixWorkspace &workspace, const SummedResults &results)
void exec() override
Executes the algorithm.
double m_maxThreshold
Maximum efficiency. Pixels with higher efficiency will be masked.
std::map< std::string, std::string > validateInputs() override
Method checking errors on ALL the inputs, before execution.
API::MatrixWorkspace_sptr mergeGroup(API::WorkspaceGroup &)
Merges the input workspace group and tries to remove masked spectra and replace the masked value with...
bool processGroups() override
Process groups and merge entries if MergeGroup property is set.
double m_minThreshold
Minimum efficiency. Pixels with lower efficiency will be masked.
SummedResults sumUnmaskedAndDeadPixels(const API::MatrixWorkspace &workspace)
Sum all detectors, excluding monitors and masked detectors.
void validateGroupInput()
Explicitly calls validateInputs and throws runtime error in case of issues in the input properties.
API::MatrixWorkspace_sptr calculateEfficiency(const API::MatrixWorkspace_sptr &, double startProgress=0.0, double stepProgress=1.0)
void init() override
Initialization method.
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
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
std::shared_ptr< WorkspaceGroup > WorkspaceGroup_sptr
shared pointer to Mantid::API::WorkspaceGroup
std::shared_ptr< Workspace > Workspace_sptr
shared pointer to Mantid::API::Workspace
Definition: Workspace_fwd.h:20
std::shared_ptr< const Workspace > Workspace_const_sptr
shared pointer to Mantid::API::Workspace (const version)
Definition: Workspace_fwd.h:22
std::shared_ptr< const MatrixWorkspace > MatrixWorkspace_const_sptr
shared pointer to the matrix workspace base class (const version)
std::shared_ptr< const WorkspaceGroup > WorkspaceGroup_const_sptr
shared pointer to Mantid::API::WorkspaceGroup, pointer to const version
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
static const std::string OUTPUT_WORKSPACE
static const std::string INPUT_WORKSPACE
static const std::string MIN_THRESHOLD
static const std::string MAX_THRESHOLD
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 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