-
Notifications
You must be signed in to change notification settings - Fork 122
/
ExtractMaskToTable.cpp
415 lines (365 loc) · 14.5 KB
/
ExtractMaskToTable.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
// NScD Oak Ridge National Laboratory, European Spallation Source,
// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
// SPDX - License - Identifier: GPL - 3.0 +
#include "MantidAlgorithms/ExtractMaskToTable.h"
#include "MantidAPI/TableRow.h"
#include "MantidAPI/WorkspaceProperty.h"
#include "MantidDataObjects/MaskWorkspace.h"
#include "MantidGeometry/Instrument.h"
#include "MantidGeometry/Instrument/DetectorInfo.h"
#include "MantidKernel/ArrayProperty.h"
namespace Mantid {
namespace Algorithms {
using namespace Mantid::API;
using namespace Mantid::DataObjects;
using namespace Mantid::Geometry;
using namespace Mantid::Kernel;
DECLARE_ALGORITHM(ExtractMaskToTable)
//----------------------------------------------------------------------------------------------
/** Declare properties
*/
void ExtractMaskToTable::init() {
declareProperty(
std::make_unique<WorkspaceProperty<MatrixWorkspace>>("InputWorkspace", "",
Direction::Input),
"A workspace whose masking is to be extracted or a MaskWorkspace. ");
declareProperty(
std::make_unique<WorkspaceProperty<TableWorkspace>>(
"MaskTableWorkspace", "", Direction::Input, PropertyMode::Optional),
"A mask table workspace containing 3 columns: "
"XMin, XMax and DetectorIDsList. ");
declareProperty(std::make_unique<WorkspaceProperty<TableWorkspace>>(
"OutputWorkspace", "", Direction::Output),
"A comma separated list or array containing a "
"list of masked detector ID's ");
declareProperty("Xmin", EMPTY_DBL(), "Minimum of X-value.");
declareProperty("Xmax", EMPTY_DBL(), "Maximum of X-value.");
}
//----------------------------------------------------------------------------------------------
/** Main execution body
*/
void ExtractMaskToTable::exec() {
// Get input properties
m_dataWS = getProperty("InputWorkspace");
if (!m_dataWS)
throw std::runtime_error(
"InputWorkspace cannot be cast to a MatrixWorkspace.");
MaskWorkspace_const_sptr maskws =
std::dynamic_pointer_cast<const MaskWorkspace>(m_dataWS);
bool m_inputIsMask = false;
if (maskws) {
g_log.notice() << "InputWorkspace " << m_dataWS->getName()
<< " is a MaskWorkspace.\n";
m_inputIsMask = true;
} else {
m_inputIsMask = false;
}
m_inputTableWS = getProperty("MaskTableWorkspace");
double xmin = getProperty("XMin");
double xmax = getProperty("XMax");
if (xmin == EMPTY_DBL() || xmax == EMPTY_DBL() || xmin >= xmax)
throw std::runtime_error(
"XMin or XMax cannot be empty. XMin must be less than XMax.");
// Create and set up output workspace
auto outws = std::make_shared<TableWorkspace>();
outws->addColumn("double", "XMin");
outws->addColumn("double", "XMax");
outws->addColumn("str", "DetectorIDsList");
setProperty("OutputWorkspace", outws);
// Optionally import the input table workspace
std::vector<detid_t> prevmaskeddetids;
if (m_inputTableWS) {
g_log.notice("Parse input masking table workspace.");
prevmaskeddetids = parseMaskTable(m_inputTableWS);
} else {
g_log.notice("No input workspace to parse.");
}
// Extract mask
std::vector<detid_t> maskeddetids;
if (m_inputIsMask)
maskeddetids = extractMaskFromMaskWorkspace();
else
maskeddetids = extractMaskFromMatrixWorkspace();
g_log.debug() << "[DB] Number of masked detectors = " << maskeddetids.size()
<< ".\n";
// Write out
if (m_inputTableWS) {
g_log.notice()
<< "About to copying input table workspace content to output workspace."
<< ".\n";
copyTableWorkspaceContent(m_inputTableWS, outws);
} else {
g_log.notice() << "There is no input workspace information to copy to "
"output workspace."
<< ".\n";
}
addToTableWorkspace(outws, maskeddetids, xmin, xmax, prevmaskeddetids);
}
//----------------------------------------------------------------------------------------------
/** Parse input TableWorkspace to get a list of detectors IDs of which detector
* are already masked
* @param masktablews :: TableWorkspace containing masking information
* @returns :: vector of detector IDs that are masked
*/
std::vector<detid_t> ExtractMaskToTable::parseMaskTable(
const DataObjects::TableWorkspace_sptr &masktablews) {
// Output vector
std::vector<detid_t> maskeddetectorids;
// Check format of mask table workspace
if (masktablews->columnCount() != 3) {
g_log.error("Mask table workspace must have more than 3 columns. First 3 "
"must be Xmin, Xmax and Spectrum List.");
return maskeddetectorids;
} else {
std::vector<std::string> colnames = masktablews->getColumnNames();
std::vector<std::string> chkcolumans(3);
chkcolumans[0] = "XMin";
chkcolumans[1] = "XMax";
chkcolumans[2] = "DetectorIDsList";
for (int i = 0; i < 3; ++i) {
if (colnames[i] != chkcolumans[i]) {
g_log.error() << "Mask table workspace " << masktablews->getName()
<< "'s " << i << "-th column name is " << colnames[i]
<< ", while it should be " << chkcolumans[i]
<< ". MaskWorkspace is invalid"
<< " and thus not used.\n";
return maskeddetectorids;
}
}
}
// Parse each row
size_t numrows = masktablews->rowCount();
double xmin, xmax;
std::string specliststr;
for (size_t i = 0; i < numrows; ++i) {
TableRow tmprow = masktablews->getRow(i);
tmprow >> xmin >> xmax >> specliststr;
std::vector<detid_t> tmpdetidvec = parseStringToVector(specliststr);
maskeddetectorids.insert(maskeddetectorids.end(), tmpdetidvec.begin(),
tmpdetidvec.end());
}
return maskeddetectorids;
}
//----------------------------------------------------------------------------------------------
/** Parse a string containing list in format (x, xx-yy, x, x, ...) to a vector
* of detid_t
* @param liststr :: string containing list to parse
* @returns :: vector genrated from input string containing the list
*/
std::vector<detid_t>
ExtractMaskToTable::parseStringToVector(const std::string &liststr) {
std::vector<detid_t> detidvec;
// Use ArrayProperty to parse the list
ArrayProperty<int> detlist("i", liststr);
if (!detlist.isValid().empty()) {
std::stringstream errss;
errss << "String '" << liststr
<< "' is unable to be converted to a list of detectors IDs. "
<< "Validation mesage: " << detlist.isValid();
g_log.error(errss.str());
throw std::runtime_error(errss.str());
}
// Convert from ArrayProperty to detectors list
size_t numdetids = detlist.size();
detidvec.reserve(numdetids);
for (size_t i = 0; i < numdetids; ++i) {
int tmpid = detlist.operator()()[i];
detidvec.emplace_back(tmpid);
g_log.debug() << "[DB] Add detector ID: " << tmpid << ".\n";
}
return detidvec;
}
//----------------------------------------------------------------------------------------------
/** Extract mask information from a workspace containing instrument
* @return vector of detector IDs of detectors that are masked
*/
std::vector<detid_t> ExtractMaskToTable::extractMaskFromMatrixWorkspace() {
// Clear input
std::vector<detid_t> maskeddetids;
// Get on hold of instrument
const auto &detectorInfo = m_dataWS->detectorInfo();
if (detectorInfo.size() == 0)
throw std::runtime_error("There is no instrument in input workspace.");
// Extract
const std::vector<detid_t> &detids = detectorInfo.detectorIDs();
for (size_t i = 0; i < detectorInfo.size(); ++i) {
bool masked = detectorInfo.isMasked(i);
if (masked) {
maskeddetids.emplace_back(detids[i]);
}
g_log.debug() << "[DB] Detector No. " << i << ": ID = " << detids[i]
<< ", Masked = " << masked << ".\n";
}
g_log.notice() << "Extract mask: There are " << maskeddetids.size()
<< " detectors that"
" are masked."
<< ".\n";
return maskeddetids;
}
//----------------------------------------------------------------------------------------------
/** Extract masked detectors from a MaskWorkspace
* @return vector of detector IDs of the detectors that are
* masked
*/
std::vector<detid_t> ExtractMaskToTable::extractMaskFromMaskWorkspace() {
// output vector
std::vector<detid_t> maskeddetids;
// Go through all spectra to find masked workspace
MaskWorkspace_const_sptr maskws =
std::dynamic_pointer_cast<const MaskWorkspace>(m_dataWS);
size_t numhist = maskws->getNumberHistograms();
for (size_t i = 0; i < numhist; ++i) {
// Rule out the spectrum without mask
if (maskws->readY(i)[0] < 1.0E-9)
continue;
const auto &detidset = maskws->getSpectrum(i).getDetectorIDs();
std::copy(detidset.cbegin(), detidset.cend(),
std::inserter(maskeddetids, maskeddetids.end()));
}
return maskeddetids;
}
//----------------------------------------------------------------------------------------------
/** Copy table workspace content from one workspace to another
* @param sourceWS :: table workspace from which the content is copied;
* @param targetWS :: table workspace to which the content is copied;
*/
void ExtractMaskToTable::copyTableWorkspaceContent(
const TableWorkspace_sptr &sourceWS, const TableWorkspace_sptr &targetWS) {
// Compare the column names. They must be exactly the same
std::vector<std::string> sourcecolnames = sourceWS->getColumnNames();
std::vector<std::string> targetcolnames = targetWS->getColumnNames();
if (sourcecolnames.size() != targetcolnames.size()) {
std::stringstream errmsg;
errmsg << "Soruce table workspace " << sourceWS->getName()
<< " has different number of columns (" << sourcecolnames.size()
<< ") than target table workspace's (" << targetcolnames.size()
<< ")";
throw std::runtime_error(errmsg.str());
}
for (size_t i = 0; i < sourcecolnames.size(); ++i) {
if (sourcecolnames[i] != targetcolnames[i]) {
std::stringstream errss;
errss << "Source and target have incompatible column name at column " << i
<< ". "
<< "Column name of source is " << sourcecolnames[i] << "; "
<< "Column name of target is " << targetcolnames[i];
throw std::runtime_error(errss.str());
}
}
// Copy over the content
size_t numrows = sourceWS->rowCount();
for (size_t i = 0; i < numrows; ++i) {
double xmin, xmax;
std::string speclist;
TableRow tmprow = sourceWS->getRow(i);
tmprow >> xmin >> xmax >> speclist;
TableRow newrow = targetWS->appendRow();
newrow << xmin << xmax << speclist;
}
}
//----------------------------------------------------------------------------------------------
/** Add a list of spectra (detector IDs) to the output table workspace.
* If a detector is masked in input MaskTableWorkspace, then it will not be
* added to a new row
* @param outws :: table workspace to write
* @param maskeddetids :: vector of detector IDs of which detectors masked
* @param xmin :: minumim x
* @param xmax :: maximum x
* @param prevmaskedids :: vector of previous masked detector IDs
*/
void ExtractMaskToTable::addToTableWorkspace(
const TableWorkspace_sptr &outws, std::vector<detid_t> maskeddetids,
double xmin, double xmax, std::vector<detid_t> prevmaskedids) {
// Sort vector of detectors ID
size_t numdetids = maskeddetids.size();
if (numdetids == 0) {
std::stringstream warnss;
warnss << "Attempting to add an empty vector of masked detectors IDs to "
"output workspace. Operation failed.";
g_log.warning(warnss.str());
return;
} else {
sort(maskeddetids.begin(), maskeddetids.end());
}
// Exclude previously masked detectors IDs from masked detectors IDs
if (!prevmaskedids.empty()) {
sort(prevmaskedids.begin(), prevmaskedids.end());
maskeddetids = subtractVector(maskeddetids, prevmaskedids);
numdetids = maskeddetids.size();
} else {
g_log.debug() << "[DB] There is no previously masked detectors."
<< ".\n";
}
if (numdetids == 0) {
// I don't know what should be done here
throw std::runtime_error("Empty detector ID list");
}
// Convert vector to string
std::stringstream spectralist;
detid_t previd = maskeddetids[0];
detid_t headid = maskeddetids[0];
for (size_t i = 1; i < numdetids; ++i) {
detid_t tmpid = maskeddetids[i];
if (tmpid == previd + 1) {
// Continuous ID
previd = tmpid;
} else if (tmpid > previd + 1) {
// Skipped ID: make a pair
if (previd == headid) {
// Single item
spectralist << " " << headid << ", ";
} else {
// Multiple items
spectralist << " " << headid << "-" << previd << ", ";
}
// New head
headid = tmpid;
previd = tmpid;
} else {
g_log.error() << "Current ID = " << tmpid << ", Previous ID = " << previd
<< ", Head ID = " << headid << ".\n";
throw std::runtime_error("Impossible! Programming logic error!");
}
} // ENDFOR (i)
// Last one
if (previd == headid)
spectralist << " " << headid;
else
spectralist << " " << headid << "-" << previd;
// Add to table workspace
std::string specliststr = spectralist.str();
TableRow newrow = outws->appendRow();
newrow << xmin << xmax << specliststr;
}
//----------------------------------------------------------------------------------------------
/** Remove the detector IDs of one vector that appear in another vector
* @param minuend :: vector with items to be removed from
* @param subtrahend :: vector containing the items to be removed from minuend
*/
std::vector<detid_t>
ExtractMaskToTable::subtractVector(std::vector<detid_t> minuend,
std::vector<detid_t> subtrahend) {
// Define some variables
std::vector<detid_t>::iterator firstsubiter, fiter;
firstsubiter = subtrahend.begin();
// Returned
std::vector<detid_t> diff;
size_t numminend = minuend.size();
for (size_t i = 0; i < numminend; ++i) {
detid_t tmpid = minuend[i];
fiter = lower_bound(firstsubiter, subtrahend.end(), tmpid);
bool exist(false);
if (fiter != subtrahend.end())
exist = *fiter == tmpid;
if (!exist) {
diff.emplace_back(tmpid);
}
firstsubiter = fiter;
}
return diff;
}
} // namespace Algorithms
} // namespace Mantid