-
Notifications
You must be signed in to change notification settings - Fork 122
/
MaskBins.cpp
154 lines (129 loc) · 5.65 KB
/
MaskBins.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
// 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/MaskBins.h"
#include "MantidAPI/Algorithm.tcc"
#include "MantidAPI/HistogramValidator.h"
#include "MantidKernel/ArrayProperty.h"
#include "MantidKernel/BoundedValidator.h"
#include <limits>
#include <sstream>
using Mantid::HistogramData::BinEdges;
namespace Mantid::Algorithms {
// Register the algorithm into the AlgorithmFactory
DECLARE_ALGORITHM(MaskBins)
using namespace Kernel;
using namespace API;
using DataObjects::EventWorkspace;
using DataObjects::EventWorkspace_const_sptr;
using DataObjects::EventWorkspace_sptr;
void MaskBins::init() {
declareWorkspaceInputProperties<MatrixWorkspace>("InputWorkspace",
"The name of the input workspace. Must contain histogram data.",
std::make_shared<HistogramValidator>());
declareProperty(std::make_unique<WorkspaceProperty<>>("OutputWorkspace", "", Direction::Output),
"The name of the Workspace containing the masked bins.");
// This validator effectively makes these properties mandatory
// Would be nice to have an explicit validator for this, but
// MandatoryValidator is already taken!
auto required = std::make_shared<BoundedValidator<double>>();
required->setUpper(std::numeric_limits<double>::max() * 0.99);
declareProperty("XMin", std::numeric_limits<double>::max(), required, "The value to start masking from.");
declareProperty("XMax", std::numeric_limits<double>::max(), required, "The value to end masking at.");
this->declareProperty(std::make_unique<ArrayProperty<int64_t>>("SpectraList"),
"Deprecated, use InputWorkspaceIndexSet.");
}
/** Execution code.
* @throw std::invalid_argument If XMax is less than XMin
*/
void MaskBins::exec() {
// Check for valid X limits
m_startX = getProperty("XMin");
m_endX = getProperty("XMax");
if (m_startX > m_endX) {
std::stringstream msg;
msg << "XMax (" << m_endX << ") must be greater than XMin (" << m_startX << ")";
g_log.error(msg.str());
throw std::invalid_argument(msg.str());
}
// Copy indices from legacy property
std::vector<int64_t> spectraList = this->getProperty("SpectraList");
if (!spectraList.empty()) {
if (!isDefault("InputWorkspaceIndexSet"))
throw std::runtime_error("Cannot provide both InputWorkspaceIndexSet and "
"SpectraList at the same time.");
setProperty("InputWorkspaceIndexSet", spectraList);
g_log.warning("The 'SpectraList' property is deprecated. Use "
"'InputWorkspaceIndexSet' instead.");
}
MatrixWorkspace_sptr inputWS;
std::tie(inputWS, indexSet) = getWorkspaceAndIndices<MatrixWorkspace>("InputWorkspace");
// Only create the output workspace if it's different to the input one
MatrixWorkspace_sptr outputWS = getProperty("OutputWorkspace");
if (outputWS != inputWS) {
outputWS = inputWS->clone();
setProperty("OutputWorkspace", outputWS);
}
if (std::dynamic_pointer_cast<const EventWorkspace>(inputWS)) {
this->execEvent();
} else {
MantidVec::difference_type startBin(0), endBin(0);
// If the binning is the same throughout, we only need to find the index
// limits once
const bool commonBins = inputWS->isCommonBins();
if (commonBins) {
auto X = inputWS->binEdges(0);
this->findIndices(X, startBin, endBin);
}
Progress progress(this, 0.0, 1.0, indexSet.size());
// Parallel running has problems with a race condition, leading to
// occaisional test failures and crashes
for (const auto wi : indexSet) {
MantidVec::difference_type startBinLoop(startBin), endBinLoop(endBin);
if (!commonBins)
this->findIndices(outputWS->binEdges(wi), startBinLoop, endBinLoop);
// Loop over masking each bin in the range
for (auto j = static_cast<int>(startBinLoop); j < static_cast<int>(endBinLoop); ++j) {
outputWS->maskBin(wi, j);
}
progress.report();
}
}
}
/** Execution code for EventWorkspaces
*/
void MaskBins::execEvent() {
MatrixWorkspace_sptr outputMatrixWS = getProperty("OutputWorkspace");
auto outputWS = std::dynamic_pointer_cast<EventWorkspace>(outputMatrixWS);
Progress progress(this, 0.0, 1.0, outputWS->getNumberHistograms() * 2);
outputWS->sortAll(Mantid::DataObjects::TOF_SORT, &progress);
PARALLEL_FOR_IF(Kernel::threadSafe(*outputWS))
for (int i = 0; i < static_cast<int>(indexSet.size()); // NOLINT
++i) {
PARALLEL_START_INTERRUPT_REGION
outputWS->getSpectrum(indexSet[i]).maskTof(m_startX, m_endX);
progress.report();
PARALLEL_END_INTERRUPT_REGION
}
PARALLEL_CHECK_INTERRUPT_REGION
outputWS->clearMRU();
}
/** Finds the indices of the bins at the limits of the range given.
* @param X :: The X vector to search
* @param startBin :: Returns the bin index including the starting value
* @param endBin :: Returns the bin index after the end value
*/
void MaskBins::findIndices(const BinEdges &X, MantidVec::difference_type &startBin,
MantidVec::difference_type &endBin) {
startBin = std::distance(X.begin(), std::upper_bound(X.cbegin(), X.cend(), m_startX));
if (startBin != 0)
--startBin;
auto last = std::lower_bound(X.cbegin(), X.cend(), m_endX);
if (last == X.cend())
--last;
endBin = std::distance(X.cbegin(), last);
}
} // namespace Mantid::Algorithms