forked from npshub/mantid
-
Notifications
You must be signed in to change notification settings - Fork 0
/
SpecialWorkspace2D.cpp
432 lines (378 loc) · 13.9 KB
/
SpecialWorkspace2D.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
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
// 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 "MantidDataObjects/SpecialWorkspace2D.h"
#include "MantidAPI/SpectraAxis.h"
#include "MantidAPI/WorkspaceFactory.h"
#include "MantidKernel/IPropertyManager.h"
#include <fstream>
#include <sstream>
using std::set;
using std::size_t;
namespace Mantid::DataObjects {
namespace {
Kernel::Logger g_log("SpecialWorkspace2D");
}
// Register the workspace
DECLARE_WORKSPACE(SpecialWorkspace2D)
//----------------------------------------------------------------------------------------------
/** Constructor, building from an instrument
*
* @param inst :: input instrument that is the base for this workspace
* @param includeMonitors :: If false the monitors are not included
* @return created SpecialWorkspace2D
*/
SpecialWorkspace2D::SpecialWorkspace2D(const Geometry::Instrument_const_sptr &inst, const bool includeMonitors) {
// Init the Workspace2D with one spectrum per detector, in the same order.
this->initialize(inst->getNumberDetectors(!includeMonitors), 1, 1);
// Copy the instrument
this->setInstrument(inst);
// Initialize the spectra-det-map, 1:1 between spectrum number and det ID
this->MatrixWorkspace::rebuildSpectraMapping(includeMonitors);
// Make the mapping, which will be used for speed later.
detID_to_WI.clear();
for (size_t wi = 0; wi < getNumberHistograms(); wi++) {
auto &dets = getSpectrum(wi).getDetectorIDs();
for (auto det : dets) {
detID_to_WI[det] = wi;
}
}
}
//----------------------------------------------------------------------------------------------
/** Constructor, building from a MatrixWorkspace
*
* @param parent :: input workspace that is the base for this workspace
* @return created SpecialWorkspace2D
*/
SpecialWorkspace2D::SpecialWorkspace2D(const API::MatrixWorkspace_const_sptr &parent) {
this->initialize(parent->getNumberHistograms(), 1, 1);
API::WorkspaceFactory::Instance().initializeFromParent(*parent, *this, false);
// Make the mapping, which will be used for speed later.
detID_to_WI.clear();
for (size_t wi = 0; wi < getNumberHistograms(); wi++) {
auto &dets = getSpectrum(wi).getDetectorIDs();
for (auto det : dets) {
detID_to_WI[det] = wi;
}
}
}
//----------------------------------------------------------------------------------------------
/** Sets the size of the workspace and initializes arrays to zero
* @param NVectors :: The number of vectors/histograms/detectors in the
* workspace
* @param XLength :: Must be 1
* @param YLength :: Must be 1
*/
void SpecialWorkspace2D::init(const size_t &NVectors, const size_t &XLength, const size_t &YLength) {
if ((XLength != 1) || (YLength != 1))
throw std::invalid_argument("SpecialWorkspace2D must have 'spectra' of length 1 only.");
// Continue with standard initialization
Workspace2D::init(NVectors, XLength, YLength);
}
void SpecialWorkspace2D::init(const HistogramData::Histogram &histogram) {
if (histogram.xMode() != HistogramData::Histogram::XMode::Points)
throw std::runtime_error("SpecialWorkspace2D can only be initialized with XMode::Points");
if (histogram.x().size() != 1)
throw std::runtime_error("SpecialWorkspace2D can only be initialized with length 1");
Workspace2D::init(histogram);
}
/**
* @return A string containing the workspace description
*/
const std::string SpecialWorkspace2D::toString() const {
std::ostringstream os;
os << "Title: " << getTitle() << "\n";
os << "Histograms: " << getNumberHistograms() << "\n";
return os.str();
}
//----------------------------------------------------------------------------------------------
/** Return the special value (Y) in the workspace at the given detector ID
*
* @param detectorID :: detector ID to look up
* @return the Y value for that detector ID.
* @throw std::invalid_argument if the detector ID was not found
*/
double SpecialWorkspace2D::getValue(const detid_t detectorID) const {
auto it = detID_to_WI.find(detectorID);
if (it == detID_to_WI.end()) {
std::ostringstream os;
os << "SpecialWorkspace2D: " << this->getName() << " Detector ID = " << detectorID
<< " Size(Map) = " << this->detID_to_WI.size() << '\n';
throw std::invalid_argument(os.str());
} else {
return this->dataY(it->second)[0];
}
}
//----------------------------------------------------------------------------------------------
/** Return the special value (Y) in the workspace at the given detector ID,
* but returns a default value instead of throwing if detector is not found.
*
* @param detectorID :: detector ID to look up
* @param defaultValue :: value returned if the ID is not found.
* @return the Y value for that detector ID.
*/
double SpecialWorkspace2D::getValue(const detid_t detectorID, const double defaultValue) const {
auto it = detID_to_WI.find(detectorID);
if (it == detID_to_WI.end())
return defaultValue;
else {
if (it->second < getNumberHistograms()) // don't let it generate an exception
{
return this->dataY(it->second)[0];
} else {
g_log.debug() << "getValue(" << detectorID << "->" << (it->second) << ", " << defaultValue
<< ") index out of range\n";
return defaultValue;
}
}
}
//----------------------------------------------------------------------------------------------
/** Set the special value (Y) in the workspace at the given detector ID
*
* @param detectorID :: detector ID to look up
* @param value :: holder for the Y value
* @param error :: the Y value for that detector ID.
* @throw std::invalid_argument if the detector ID was not found
*/
void SpecialWorkspace2D::setValue(const detid_t detectorID, const double value, const double error) {
auto it = detID_to_WI.find(detectorID);
if (it == detID_to_WI.end()) {
std::stringstream msg;
msg << "SpecialWorkspace2D::setValue(): Input Detector ID = " << detectorID << " Is Invalid";
throw std::invalid_argument(msg.str());
} else {
this->dataY(it->second)[0] = value;
this->dataE(it->second)[0] = error;
}
}
/**
* Set the special value (Y) in the workspace at the given detector IDs. This is
*a
* convenience function that just calls setValue(const detid_t, double, double).
*
* @param detectorIDs :: detector IDs to look up
* @param value :: holder for the Y value
* @param error :: the Y value for the detector IDs.
* @throw std::invalid_argument if any of the detector IDs are not found
*/
void SpecialWorkspace2D::setValue(const set<detid_t> &detectorIDs, const double value, const double error) {
for (auto detectorID : detectorIDs) {
this->setValue(detectorID, value, error);
}
}
//----------------------------------------------------------------------------------------------
/** Return the detector ID at the given workspace index (i.e.,
*spectrum/histogram index)
*
* @param workspaceIndex
* @return
*/
set<detid_t> SpecialWorkspace2D::getDetectorIDs(const std::size_t workspaceIndex) const {
if (size_t(workspaceIndex) > this->getNumberHistograms())
throw std::invalid_argument("SpecialWorkspace2D::getDetectorIDs(): Invalid workspaceIndex given.");
return this->getSpectrum(workspaceIndex).getDetectorIDs();
}
//--------------------------------------------------------------------------------------------
/** Return the result of operator &
* @ parameter
* @ return
*/
void SpecialWorkspace2D::binaryOperation(std::shared_ptr<const SpecialWorkspace2D> &ws,
const unsigned int operatortype) {
// 1. Check compatibility between this and input workspace
if (!this->isCompatible(ws)) {
throw std::invalid_argument("Two SpecialWorkspace2D objects are not compatible!");
}
switch (operatortype) {
case BinaryOperator::AND:
this->binaryAND(ws);
break;
case BinaryOperator::OR:
this->binaryOR(ws);
break;
case BinaryOperator::XOR:
this->binaryXOR(ws);
break;
default:
throw std::invalid_argument("Invalid Operator");
break;
}
}
//--------------------------------------------------------------------------------------------
/** Return the result of operator &
* @ parameter
* @ return
*/
void SpecialWorkspace2D::binaryOperation(const unsigned int operatortype) {
switch (operatortype) {
case BinaryOperator::NOT:
this->binaryNOT();
break;
default:
g_log.error() << "Operator " << operatortype << " Is Not Valid In BinaryOperation(operatortype)\n";
throw std::invalid_argument("Invalid Operator");
break;
}
}
/** AND operator
*
*/
void SpecialWorkspace2D::binaryAND(const std::shared_ptr<const SpecialWorkspace2D> &ws) {
for (size_t i = 0; i < this->getNumberHistograms(); i++) {
double y1 = this->dataY(i)[0];
double y2 = ws->dataY(i)[0];
if (y1 < 1.0E-10 || y2 < 1.0E-10) {
this->dataY(i)[0] = 0.0;
} else {
this->dataY(i)[0] += y2;
}
}
}
/** OR operator
*
*/
void SpecialWorkspace2D::binaryOR(const std::shared_ptr<const SpecialWorkspace2D> &ws) {
for (size_t i = 0; i < this->getNumberHistograms(); i++) {
double y1 = this->dataY(i)[0];
double y2 = ws->dataY(i)[0];
double max = y1;
if (y2 > y1) {
max = y2;
}
this->dataY(i)[0] = max;
/*
if (y1 < 1.0E-10 && y2 < 1.0E-10){
this->dataY(i)[0] = 0.0;
} else {
this->dataY(i)[0] += y2;
}
*/
}
}
/** Excluded Or operator
*
*/
void SpecialWorkspace2D::binaryXOR(const std::shared_ptr<const SpecialWorkspace2D> &ws) {
for (size_t i = 0; i < this->getNumberHistograms(); i++) {
double y1 = this->dataY(i)[0];
double y2 = ws->dataY(i)[0];
if (y1 < 1.0E-10 && y2 < 1.0E-10) {
this->dataY(i)[0] = 0.0;
} else if (y1 > 1.0E-10 && y2 > 1.0E-10) {
this->dataY(i)[0] = 0.0;
} else {
this->dataY(i)[0] = 1.0;
}
}
}
/**
* NOT operator
*/
void SpecialWorkspace2D::binaryNOT() {
for (size_t i = 0; i < this->getNumberHistograms(); i++) {
double y1 = this->dataY(i)[0];
if (y1 < 1.0E-10) {
this->dataY(i)[0] = 1.0;
} else {
this->dataY(i)[0] = 0.0;
}
}
}
//----------------------------------------------------------------------------------------------
/* Check 2 SpecialWorkspace2D are compatible
* @ parameter
* @ return
*/
bool SpecialWorkspace2D::isCompatible(const std::shared_ptr<const SpecialWorkspace2D> &ws) {
// 1. Check number of histogram
size_t numhist1 = this->getNumberHistograms();
size_t numhist2 = ws->getNumberHistograms();
if (numhist1 != numhist2) {
g_log.debug() << "2 Workspaces have different number of histograms: " << numhist1 << " vs. " << numhist2 << '\n';
return false;
}
// 2. Check detector ID
for (size_t ispec = 0; ispec < numhist1; ispec++) {
set<detid_t> ids1 = this->getSpectrum(ispec).getDetectorIDs();
set<detid_t> ids2 = ws->getSpectrum(ispec).getDetectorIDs();
if (ids1.size() != ids2.size()) {
g_log.debug() << "Spectra " << ispec << ": 2 Workspaces have different number of detectors " << ids1.size()
<< " vs. " << ids2.size() << '\n';
return false;
} else if (ids1.empty()) {
g_log.debug() << "Spectra " << ispec << ": 2 Workspaces both have 0 detectors. \n";
return false;
} else if (*ids1.begin() != *ids2.begin()) {
g_log.debug() << "Spectra " << ispec << ": 2 Workspaces have different Detector ID " << *ids1.begin() << " vs. "
<< *ids2.begin() << '\n';
return false;
}
} // false
return true;
}
//----------------------------------------------------------------------------------------------
/** Duplicate SpecialWorkspace2D
*/
void SpecialWorkspace2D::copyFrom(std::shared_ptr<const SpecialWorkspace2D> sourcews) {
// Check
if (this->getNumberHistograms() != sourcews->getNumberHistograms()) {
throw std::invalid_argument("Incompatible number of histograms");
}
// Copy data
for (size_t ispec = 0; ispec < this->getNumberHistograms(); ispec++) {
// 1.1 Check size
const MantidVec &inx = sourcews->readX(ispec);
const MantidVec &iny = sourcews->readY(ispec);
const MantidVec &ine = sourcews->readE(ispec);
MantidVec &outx = this->dataX(ispec);
MantidVec &outy = this->dataY(ispec);
MantidVec &oute = this->dataE(ispec);
if (inx.size() != outx.size() || iny.size() != outy.size() || ine.size() != oute.size()) {
throw std::invalid_argument("X, Y, E size different within spectrum");
}
// 1.2 Copy data
for (size_t i = 0; i < inx.size(); i++) {
outx[i] = inx[i];
}
for (size_t i = 0; i < iny.size(); i++) {
outy[i] = iny[i];
oute[i] = ine[i];
}
}
// Copy detector map
this->detID_to_WI = sourcews->detID_to_WI;
}
} // namespace Mantid::DataObjects
/// @cond TEMPLATE
namespace Mantid::Kernel {
template <>
DLLExport Mantid::DataObjects::SpecialWorkspace2D_sptr
IPropertyManager::getValue<Mantid::DataObjects::SpecialWorkspace2D_sptr>(const std::string &name) const {
auto *prop =
dynamic_cast<PropertyWithValue<Mantid::DataObjects::SpecialWorkspace2D_sptr> *>(getPointerToProperty(name));
if (prop) {
return *prop;
} else {
std::string message =
"Attempt to assign property " + name + " to incorrect type. Expected shared_ptr<SpecialWorkspace2D>.";
throw std::runtime_error(message);
}
}
template <>
DLLExport Mantid::DataObjects::SpecialWorkspace2D_const_sptr
IPropertyManager::getValue<Mantid::DataObjects::SpecialWorkspace2D_const_sptr>(const std::string &name) const {
auto *prop =
dynamic_cast<PropertyWithValue<Mantid::DataObjects::SpecialWorkspace2D_sptr> *>(getPointerToProperty(name));
if (prop) {
return prop->operator()();
} else {
std::string message =
"Attempt to assign property " + name + " to incorrect type. Expected const shared_ptr<SpecialWorkspace2D>.";
throw std::runtime_error(message);
}
}
} // namespace Mantid::Kernel
/// @endcond TEMPLATE