Mantid
Loading...
Searching...
No Matches
CreatePSDBleedMask.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 +
8#include "MantidAPI/Run.h"
15
16#include <cfloat>
17#include <iterator>
18#include <list>
19#include <map>
20
21namespace Mantid::Algorithms {
22
23// Register the class
24DECLARE_ALGORITHM(CreatePSDBleedMask)
25
26const std::string CreatePSDBleedMask::category() const { return "Diagnostics"; }
27
31
34
40
41 declareProperty(std::make_unique<WorkspaceProperty<>>("InputWorkspace", "", Direction::Input),
42 "The name of the input workspace.");
43 declareProperty(std::make_unique<WorkspaceProperty<>>("OutputWorkspace", "", Direction::Output),
44 "The name of the output MaskWorkspace which will contain the result "
45 "masks.");
46 auto mustBePosDbl = std::make_shared<BoundedValidator<double>>();
47 mustBePosDbl->setLower(0.0);
48 declareProperty("MaxTubeFramerate", -1.0, mustBePosDbl, "The maximum rate allowed for a tube in counts/us/frame.");
49 auto mustBePosInt = std::make_shared<BoundedValidator<int>>();
50 mustBePosInt->setLower(0);
51 declareProperty("NIgnoredCentralPixels", 80, mustBePosInt, "The number of pixels about the centre to ignore.");
52 declareProperty("NumberOfFailures", 0, std::make_shared<Kernel::NullValidator>(),
53 "An output property containing the number of masked tubes", Direction::Output);
54}
55
59
60 MatrixWorkspace_const_sptr inputWorkspace = getProperty("InputWorkspace");
61 // We require the number of good frames. Check that we have this
62 if (!inputWorkspace->run().hasProperty("goodfrm")) {
63 throw std::invalid_argument("InputWorkspace does not contain the number of \"good frames\".\n"
64 "(The sample log named: goodfrm with value, specifying number of good "
65 "frames)");
66 }
68 dynamic_cast<Kernel::PropertyWithValue<int> *>(inputWorkspace->run().getProperty("goodfrm"));
69 if (!frameProp) {
70 throw std::invalid_argument("InputWorkspace has the number of \"good "
71 "frames\" property (goodfrm log value)"
72 "but this property value is not integer.");
73 }
74
75 int goodFrames = (*frameProp)();
76
77 // Store the other properties
78 double maxFramerate = getProperty("MaxTubeFramerate");
79
80 // Multiply by the frames to save a division for each bin when we loop over
81 // them later
82 double maxRate = maxFramerate * goodFrames;
83 int numIgnoredPixels = getProperty("NIgnoredCentralPixels");
84
85 // This algorithm assumes that the instrument geometry is tube based, i.e. the
86 // parent CompAssembly
87 // of the lowest detector in the tree is a "tube" and that all pixels in a
88 // tube are consecutively ordered
89 // with respect to spectra number
90 const auto numSpectra = static_cast<int>(inputWorkspace->getNumberHistograms());
91 // Keep track of a map of tubes to lists of indices
92 using TubeIndex = std::map<Geometry::ComponentID, std::vector<int>>;
93 TubeIndex tubeMap;
94
95 API::Progress progress(this, 0.0, 1.0, numSpectra);
96
97 const auto &spectrumInfo = inputWorkspace->spectrumInfo();
98
99 // NOTE: This loop is intentionally left unparallelized as the majority of the
100 // work requires a lock around it which actually slows down the loop.
101 // Another benefit of keep it serial is losing the need for a call to 'sort'
102 // when
103 // performing the bleed test as the list of indices will already be in the
104 // correct
105 // order
106 for (int i = 0; i < numSpectra; ++i) {
107 if (!spectrumInfo.hasDetectors(i))
108 continue;
109
110 if (spectrumInfo.isMonitor(i))
111 continue;
112
113 auto &det = spectrumInfo.detector(i);
114 std::shared_ptr<const Geometry::IComponent> parent;
115
116 if (!spectrumInfo.hasUniqueDetector(i)) {
117 const auto &group = dynamic_cast<const Geometry::DetectorGroup &>(det);
118 parent = group.getDetectors().front()->getParent();
119 } else {
120 parent = det.getParent();
121 }
122
123 if (!parent)
124 continue;
125
126 Geometry::ComponentID parentID = parent->getComponentID();
127 // Already have this component
128 if (tubeMap.find(parentID) != tubeMap.end()) {
129 tubeMap[parentID].emplace_back(i);
130 }
131 // New tube
132 else {
133 tubeMap.emplace(parentID, TubeIndex::mapped_type(1, i));
134 }
135
136 progress.report();
137 }
138
139 // Now process the tubes in parallel
140 const auto numTubes = static_cast<int>(tubeMap.size());
141 g_log.information() << "Found " << numTubes << " tubes.\n";
142 int numSpectraMasked(0), numTubesMasked(0);
143 // Create a mask workspace for output
144 MaskWorkspace_sptr outputWorkspace = this->generateEmptyMask(inputWorkspace);
145
146 progress.resetNumSteps(numTubes, 0, 1);
147
148 PARALLEL_FOR_IF(Kernel::threadSafe(*inputWorkspace, *outputWorkspace))
149 for (int i = 0; i < numTubes; ++i) {
151 auto current = std::next(tubeMap.begin(), i);
152 const TubeIndex::mapped_type tubeIndices = current->second;
153 bool mask = performBleedTest(tubeIndices, inputWorkspace, maxRate, numIgnoredPixels);
154 if (mask) {
155 maskTube(tubeIndices, outputWorkspace);
157 numSpectraMasked += static_cast<int>(tubeIndices.size());
159 numTubesMasked += 1;
160 }
161
162 progress.report("Performing Bleed Test");
163
165 }
167
168 g_log.information() << numTubesMasked << " tube(s) failed the bleed tests.";
169 if (numTubesMasked > 0) {
170 g_log.information() << " The " << numSpectraMasked << " spectra have been masked on the output workspace.\n";
171 } else {
172 g_log.information() << '\n';
173 }
174
175 setProperty("NumberOfFailures", numSpectraMasked);
176 setProperty("OutputWorkspace", outputWorkspace);
177}
178
188bool CreatePSDBleedMask::performBleedTest(const std::vector<int> &tubeIndices,
189 const API::MatrixWorkspace_const_sptr &inputWS, double maxRate,
190 int numIgnoredPixels) {
191
192 // Require ordered pixels so that we can define the centre.
193 // This of course assumes that the pixel IDs increase monotonically with the
194 // workspace index
195 // and that the above loop that searched for the tubes was NOT run in parallel
196 const size_t numSpectra(tubeIndices.size());
197 const size_t midIndex(numSpectra / 2);
198 const size_t topEnd(midIndex - numIgnoredPixels / 2);
199 const size_t bottomBegin(midIndex + numIgnoredPixels / 2);
200
203 bool isRawCounts = !(inputWS->isDistribution());
204
205 const auto numBins = static_cast<int>(inputWS->blocksize());
206 std::vector<double> totalRate(numBins, 0.0);
207 size_t top = 0, bot = bottomBegin;
208 for (; top < topEnd; ++top, ++bot) {
209 const int topIndex = tubeIndices[top];
210 const int botIndex = tubeIndices[bot];
211 auto &topY = inputWS->y(topIndex);
212 auto &botY = inputWS->y(botIndex);
213 auto &topX = inputWS->x(topIndex);
214 auto &botX = inputWS->x(botIndex);
215 for (int j = 0; j < numBins; ++j) {
216 double topRate(topY[j]), botRate(botY[j]);
217 if (isRawCounts) {
218 topRate /= (topX[j + 1] - topX[j]);
219 botRate /= (botX[j + 1] - botX[j]);
220 }
221 totalRate[j] += topRate + botRate;
222 // If by now any have hit the allowed maximum then mark this to be masked
223 if (totalRate[j] > maxRate) {
224 return true;
225 }
226 }
227 }
228
229 if (top != topEnd) {
230 g_log.error() << "Error in tube processing, loop variable has an unexpected value.\n";
231 throw std::runtime_error("top != topEnd in CreatePSDBleedMask::performBleedTest()");
232 }
233 if (bot != numSpectra) {
234 g_log.error() << "Error in tube processing, loop variable has an unexpected value.\n";
235 throw std::runtime_error("bot != numSpectra in CreatePSDBleedMask::performBleedTest()");
236 }
237
238 return false;
239}
240
246void CreatePSDBleedMask::maskTube(const std::vector<int> &tubeIndices, const API::MatrixWorkspace_sptr &workspace) {
247 const double deadValue(1.0); // delete the data
248 for (auto tubeIndice : tubeIndices) {
249 workspace->mutableY(tubeIndice)[0] = deadValue;
250 }
251}
252} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
IPeaksWorkspace_sptr workspace
Definition: IndexPeaks.cpp:114
double top
Definition: LineProfile.cpp:78
#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
Helper class for reporting progress from algorithms.
Definition: Progress.h:25
A property class for workspaces.
This algorithm implements a "bleed" diagnostic for PSD detectors (i.e.
void exec() override
Execute the algorithm.
void init() override
Initialize the algorithm properties.
bool performBleedTest(const std::vector< int > &tubeIndices, const API::MatrixWorkspace_const_sptr &inputWS, double maxRate, int numIgnoredPixels)
Process a tube.
CreatePSDBleedMask()
Default constructor.
void maskTube(const std::vector< int > &tubeIndices, const API::MatrixWorkspace_sptr &workspace)
Mask a tube with the given workspace indices.
DataObjects::MaskWorkspace_sptr generateEmptyMask(const API::MatrixWorkspace_const_sptr &inputWS)
Create a masking workspace to return.
Holds a collection of detectors.
Definition: DetectorGroup.h:28
std::vector< IDetector_const_sptr > getDetectors() const
What detectors are contained in the group?
base class for Geometric IComponent
Definition: IComponent.h:51
virtual ComponentID getComponentID() const =0
Returns the ComponentID - a unique identifier of the component.
BoundedValidator is a validator that requires the values to be between upper or lower bounds,...
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void error(const std::string &msg)
Logs at error level.
Definition: Logger.cpp:77
void information(const std::string &msg)
Logs at information level.
Definition: Logger.cpp:105
The concrete, templated class for properties.
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::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
STL namespace.
Describes the direction (within an algorithm) of a Property.
Definition: Property.h:50
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54