Mantid
Loading...
Searching...
No Matches
MaskPeaksWorkspace.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 "MantidAPI/TableRow.h"
18
19#include <boost/math/special_functions/round.hpp>
20
21namespace Mantid::Crystal {
22
23// Register the class into the algorithm factory
24DECLARE_ALGORITHM(MaskPeaksWorkspace)
25
26using namespace Kernel;
27using namespace API;
28using namespace DataObjects;
29using namespace Geometry;
30using std::string;
31
33MaskPeaksWorkspace::MaskPeaksWorkspace() : m_xMin{0}, m_xMax{0}, m_yMin{0}, m_yMax{0}, m_tofMin{0}, m_tofMax{0} {}
34
39
40 declareProperty(std::make_unique<WorkspaceProperty<MatrixWorkspace>>("InputWorkspace", "", Direction::Input,
41 std::make_shared<InstrumentValidator>()),
42 "A workspace containing one or more rectangular area "
43 "detectors. Each spectrum needs to correspond to only one "
44 "pixelID (e.g. no grouping or previous calls to "
45 "SumNeighbours).");
46 declareProperty(std::make_unique<WorkspaceProperty<PeaksWorkspace>>("InPeaksWorkspace", "", Direction::Input),
47 "The name of the workspace that will be created. Can replace "
48 "the input workspace.");
49 declareProperty("XMin", -2, "Minimum of X (col) Range to mask peak relative to peak's center");
50 declareProperty("XMax", 2, "Maximum of X (col) Range to mask peak relative to peak's center");
51 declareProperty("YMin", -2, "Minimum of Y (row) Range to mask peak relative to peak's center");
52 declareProperty("YMax", 2, "Maximum of Y (row) Range to mask peak relative to peak's center");
53 declareProperty("TOFMin", EMPTY_DBL(),
54 "Optional(all TOF if not specified): "
55 "Minimum TOF relative to peak's "
56 "center TOF.");
57 declareProperty("TOFMax", EMPTY_DBL(),
58 "Optional(all TOF if not specified): "
59 "Maximum TOF relative to peak's "
60 "center TOF.");
61}
62
70
71 PeaksWorkspace_const_sptr peaksW = getProperty("InPeaksWorkspace");
72
73 // To get the workspace index from the detector ID
74 const detid2index_map pixel_to_wi = m_inputW->getDetectorIDToWorkspaceIndexMap();
75 // Get some stuff from the input workspace
76 Geometry::Instrument_const_sptr inst = m_inputW->getInstrument();
77
78 // Init a table workspace
79 DataObjects::TableWorkspace_sptr tablews = std::make_shared<DataObjects::TableWorkspace>();
80 tablews->addColumn("double", "XMin");
81 tablews->addColumn("double", "XMax");
82 tablews->addColumn("str", "SpectraList");
83
84 // Loop over peaks
85 const std::vector<Peak> &peaks = peaksW->getPeaks();
86 PARALLEL_FOR_IF(Kernel::threadSafe(*m_inputW, *peaksW, *tablews))
87 for (int i = 0; i < static_cast<int>(peaks.size()); i++) { // NOLINT
89 const Peak &peak = peaks[i];
90 // get the peak location on the detector
91 double col = peak.getCol();
92 double row = peak.getRow();
93 int xPeak = boost::math::iround(col) - 1;
94 int yPeak = boost::math::iround(row) - 1;
95 g_log.debug() << "Generating information for peak at x=" << xPeak << " y=" << yPeak << "\n";
96
97 // the detector component for the peak will have all pixels that we mask
98 const string &bankName = peak.getBankName();
99 if (bankName == "None")
100 continue;
101 Geometry::IComponent_const_sptr comp = inst->getComponentByName(bankName);
102 if (!comp) {
103 g_log.debug() << "Component " + bankName + " does not exist in instrument\n";
104 continue;
105 }
106
107 // determine the range in time-of-flight
108 double x0;
109 double xf;
110 bool tofRangeSet(false);
111 size_t wi = this->getWkspIndex(pixel_to_wi, comp, xPeak, yPeak);
112 if (wi != static_cast<size_t>(EMPTY_INT())) { // scope limit the workspace index
113 this->getTofRange(x0, xf, peak.getTOF(), m_inputW->x(wi));
114 tofRangeSet = true;
115 }
116
117 // determine the spectrum numbers to mask
118 std::set<size_t> spectra;
119 for (int ix = m_xMin; ix <= m_xMax; ix++) {
120 for (int iy = m_yMin; iy <= m_yMax; iy++) {
121 // Find the pixel ID at that XY position on the rectangular detector
122 size_t wj = this->getWkspIndex(pixel_to_wi, comp, xPeak + ix, yPeak + iy);
123 if (wj == static_cast<size_t>(EMPTY_INT()))
124 continue;
125 spectra.insert(wj);
126 if (!tofRangeSet) { // scope limit the workspace index
127 this->getTofRange(x0, xf, peak.getTOF(), m_inputW->x(wj));
128 tofRangeSet = true;
129 }
130 }
131 }
132
133 // sanity check the results
134 if (!tofRangeSet) {
135 g_log.warning() << "Failed to set time-of-flight range for peak (x=" << xPeak << ", y=" << yPeak
136 << ", tof=" << peak.getTOF() << ")\n";
137 } else if (spectra.empty()) {
138 g_log.warning() << "Failed to find spectra for peak (x=" << xPeak << ", y=" << yPeak << ", tof=" << peak.getTOF()
139 << ")\n";
140 continue;
141 } else
142 PARALLEL_CRITICAL(tablews) {
143 // append to the table workspace
144 API::TableRow newrow = tablews->appendRow();
145 newrow << x0 << xf << Kernel::Strings::toString(spectra);
146 }
148 } // end loop over peaks
150
151 // Mask bins
152 auto maskbinstb = createChildAlgorithm("MaskBinsFromTable", 0.5, 1.0, true);
153 maskbinstb->setProperty("InputWorkspace", m_inputW);
154 maskbinstb->setPropertyValue("OutputWorkspace", m_inputW->getName());
155 maskbinstb->setProperty("MaskingInformation", tablews);
156 maskbinstb->execute();
157}
158
160 m_inputW = getProperty("InputWorkspace");
161
162 m_xMin = getProperty("XMin");
163 m_xMax = getProperty("XMax");
164 if (m_xMin >= m_xMax)
165 throw std::runtime_error("Must specify Xmin<Xmax");
166
167 m_yMin = getProperty("YMin");
168 m_yMax = getProperty("YMax");
169 if (m_yMin >= m_yMax)
170 throw std::runtime_error("Must specify Ymin<Ymax");
171
172 // Get the value of TOF range to mask
173 m_tofMin = getProperty("TOFMin");
174 m_tofMax = getProperty("TOFMax");
175 if ((!isEmpty(m_tofMin)) && (!isEmpty(m_tofMax))) {
176 if (m_tofMin >= m_tofMax)
177 throw std::runtime_error("Must specify TOFMin < TOFMax");
178 } else if ((!isEmpty(m_tofMin)) || (!isEmpty(m_tofMax))) // check if only one is empty
179 {
180 throw std::runtime_error("Must specify both TOFMin and TOFMax or neither");
181 }
182}
183
185 const int x, const int y) {
186 Geometry::RectangularDetector_const_sptr det = std::dynamic_pointer_cast<const Geometry::RectangularDetector>(comp);
187 Geometry::Instrument_const_sptr Iptr = m_inputW->getInstrument();
188
189 if (det) {
190 if (x >= det->xpixels() || x < 0 || y >= det->ypixels() || y < 0) {
191 // throw std::runtime_error("Failed to find workspace index for x=" + std::to_string(x) + " y=" +
192 // std::to_string(y) + "(max x=" + std::to_string(det->xpixels()) +
193 // ", max y=" + std::to_string(det->ypixels()) + ")"); // Useful for debugging
194 return EMPTY_INT();
195 }
196
197 int pixelID = det->getAtXY(x, y)->getID();
198
199 // Find the corresponding workspace index, if any
200 auto wiEntry = pixel_to_wi.find(pixelID);
201 if (wiEntry == pixel_to_wi.end()) {
202 std::stringstream msg;
203 msg << "Failed to find workspace index for x=" << x << " y=" << y;
204 throw std::runtime_error(msg.str());
205 }
206 return wiEntry->second;
207 } else {
208 std::vector<Geometry::IComponent_const_sptr> children;
209 std::shared_ptr<const Geometry::ICompAssembly> asmb =
210 std::dynamic_pointer_cast<const Geometry::ICompAssembly>(comp);
211 asmb->getChildren(children, false);
212
213 std::shared_ptr<const Geometry::ICompAssembly> asmb2 =
214 std::dynamic_pointer_cast<const Geometry::ICompAssembly>(children[0]);
215 std::vector<Geometry::IComponent_const_sptr> grandchildren;
216 asmb2->getChildren(grandchildren, false);
217
218 auto NROWS = static_cast<int>(grandchildren.size());
219 auto NCOLS = static_cast<int>(children.size());
220
221 std::ostringstream msg;
222 if (Iptr->getName().compare("CORELLI") == 0) {
223 msg << "Instrument is CORELLI\n";
224 // CORELLI has one extra layer than WISH
225 std::shared_ptr<const Geometry::ICompAssembly> asmb3 =
226 std::dynamic_pointer_cast<const Geometry::ICompAssembly>(grandchildren[0]);
227 std::vector<Geometry::IComponent_const_sptr> greatgrandchildren;
228 asmb3->getChildren(greatgrandchildren, false);
229 // update for CORELLI
230 NCOLS = static_cast<int>(grandchildren.size());
231 NROWS = static_cast<int>(greatgrandchildren.size());
232 } else {
233 msg << "Instrument is WISH\n";
234 }
235
236 // Wish pixels and tubes start at 1 not 0
237 if (x - 1 >= NCOLS || x - 1 < 0 || y - 1 >= NROWS || y - 1 < 0) {
238 // useful for future dev in plan
239 // msg << "--(x,y) = (" << x << "," << y << ")\n"
240 // << "--NCOLS = " << NCOLS << "\n"
241 // << "--NROWS = " << NROWS << "\n";
242 // g_log.warning() << msg.str();
243 return EMPTY_INT();
244 }
245
246 std::string bankName = comp->getName();
247 auto it = pixel_to_wi.find(findPixelID(bankName, x, y));
248 if (it == pixel_to_wi.end())
249 return EMPTY_INT();
250 return (it->second);
251 }
252}
253
261void MaskPeaksWorkspace::getTofRange(double &tofMin, double &tofMax, const double tofPeak,
262 const HistogramData::HistogramX &tof) {
263 tofMin = tof.front();
264 tofMax = tof.back() - 1;
265 if (!isEmpty(m_tofMin)) {
266 tofMin = tofPeak + m_tofMin;
267 }
268 if (!isEmpty(m_tofMax)) {
269 tofMax = tofPeak + m_tofMax;
270 }
271}
272
281int MaskPeaksWorkspace::findPixelID(const std::string &bankName, int col, int row) {
282 Geometry::Instrument_const_sptr Iptr = m_inputW->getInstrument();
283 std::shared_ptr<const IComponent> parent = Iptr->getComponentByName(bankName);
284
285 if (parent->type() == "RectangularDetector") {
286 std::shared_ptr<const RectangularDetector> RDet = std::dynamic_pointer_cast<const RectangularDetector>(parent);
287
288 std::shared_ptr<Detector> pixel = RDet->getAtXY(col, row);
289 return pixel->getID();
290 } else if (Iptr->getName().compare("CORELLI") == 0) {
291 // Checking for CORELLI
292 // pixel full name example
293 // /CORELLI/A row/bank10/sixteenpack/tube10/pixel23
294 // ^ the extra layer that makes CORELLI different from WISH
295 std::ostringstream pixelString;
296 pixelString << parent->getFullName() // /CORELLI/A row/bank10
297 << "/sixteenpack" // /sixteenpack
298 << "/tube" << col // /tube10
299 << "/pixel" << row; // /pixel23
300 std::shared_ptr<const Geometry::IComponent> component = Iptr->getComponentByName(pixelString.str());
301 std::shared_ptr<const Detector> pixel = std::dynamic_pointer_cast<const Detector>(component);
302 //
303 return pixel->getID();
304 } else {
305 std::string bankName0 = bankName;
306 // Only works for WISH
307 bankName0.erase(0, 4);
308 std::ostringstream pixelString;
309 pixelString << Iptr->getName() << "/" << bankName0 << "/" << bankName << "/tube" << std::setw(3)
310 << std::setfill('0') << col << "/pixel" << std::setw(4) << std::setfill('0') << row;
311 std::shared_ptr<const Geometry::IComponent> component = Iptr->getComponentByName(pixelString.str());
312 std::shared_ptr<const Detector> pixel = std::dynamic_pointer_cast<const Detector>(component);
313 return pixel->getID();
314 }
315}
316} // namespace Mantid::Crystal
#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_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.
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
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
static bool isEmpty(const NumT toCheck)
checks that the value was not set by users, uses the value in empty double/int.
TableRow represents a row in a TableWorkspace.
Definition: TableRow.h:39
A property class for workspaces.
double m_tofMin
The start of the box around the peak in tof.
int findPixelID(const std::string &bankName, int col, int row)
std::size_t getWkspIndex(const detid2index_map &pixel_to_wi, const Geometry::IComponent_const_sptr &comp, const int x, const int y)
void init() override
Initialisation method.
int m_xMin
The start of the X range for fitting.
int m_yMax
The end of the Y range for fitting.
void exec() override
Executes the algorithm.
double m_tofMax
The end of the box around the peak in tof.
void getTofRange(double &tofMin, double &tofMax, const double tofPeak, const HistogramData::HistogramX &tof)
API::MatrixWorkspace_sptr m_inputW
A pointer to the input workspace.
void retrieveProperties()
Read in all the input parameters.
int m_xMax
The end of the X range for fitting.
int m_yMin
The start of the Y range for fitting.
Structure describing a single-crystal peak.
Definition: Peak.h:34
int getCol() const override
For RectangularDetectors only, returns the column (x) of the pixel of the detector or -1 if not found...
Definition: Peak.cpp:341
void debug(const std::string &msg)
Logs at debug level.
Definition: Logger.cpp:114
void warning(const std::string &msg)
Logs at warning level.
Definition: Logger.cpp:86
std::shared_ptr< const PeaksWorkspace > PeaksWorkspace_const_sptr
Typedef for a shared pointer to a const peaks workspace.
std::shared_ptr< TableWorkspace > TableWorkspace_sptr
shared pointer to Mantid::DataObjects::TableWorkspace
std::shared_ptr< const IComponent > IComponent_const_sptr
Typdef of a shared pointer to a const IComponent.
Definition: IComponent.h:161
std::shared_ptr< const RectangularDetector > RectangularDetector_const_sptr
std::shared_ptr< const Instrument > Instrument_const_sptr
Shared pointer to an const instrument object.
std::string toString(const T &value)
Convert a number to a string.
Definition: Strings.cpp:703
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
std::unordered_map< detid_t, size_t > detid2index_map
Map with key = detector ID, value = workspace index.
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