-
Notifications
You must be signed in to change notification settings - Fork 122
/
FilterByLogValue.cpp
235 lines (195 loc) · 8.74 KB
/
FilterByLogValue.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
// 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/FilterByLogValue.h"
#include "MantidAPI/Run.h"
#include "MantidDataObjects/WorkspaceCreation.h"
#include "MantidKernel/BoundedValidator.h"
#include "MantidKernel/ITimeSeriesProperty.h"
#include "MantidKernel/ListValidator.h"
#include "MantidKernel/MandatoryValidator.h"
namespace Mantid {
namespace Algorithms {
// Register the algorithm into the algorithm factory
DECLARE_ALGORITHM(FilterByLogValue)
using namespace Kernel;
using namespace DataObjects;
using namespace API;
using DataObjects::EventList;
using DataObjects::EventWorkspace;
using DataObjects::EventWorkspace_const_sptr;
using DataObjects::EventWorkspace_sptr;
using Types::Core::DateAndTime;
std::string CENTRE("Centre");
std::string LEFT("Left");
void FilterByLogValue::init() {
declareProperty(std::make_unique<WorkspaceProperty<EventWorkspace>>("InputWorkspace", "", Direction::Input),
"An input event workspace");
declareProperty(std::make_unique<WorkspaceProperty<EventWorkspace>>("OutputWorkspace", "", Direction::Output),
"The name to use for the output workspace");
declareProperty("LogName", "", std::make_shared<MandatoryValidator<std::string>>(),
"Name of the sample log to use to filter.\n"
"For example, the pulse charge is recorded in 'ProtonCharge'.");
declareProperty("MinimumValue", Mantid::EMPTY_DBL(), "Minimum log value for which to keep events.");
declareProperty("MaximumValue", Mantid::EMPTY_DBL(), "Maximum log value for which to keep events.");
auto min = std::make_shared<BoundedValidator<double>>();
min->setLower(0.0);
declareProperty("TimeTolerance", 0.0, min,
"Tolerance, in seconds, for the event times to keep. A good "
"value is 1/2 your measurement interval. \n"
"For a single log value at time T, all events between "
"T+-Tolerance are kept.\n"
"If there are several consecutive log values matching the "
"filter, events between T1-Tolerance and T2+Tolerance are "
"kept.");
std::vector<std::string> types(2);
types[0] = CENTRE;
types[1] = LEFT;
declareProperty("LogBoundary", types[0], std::make_shared<StringListValidator>(types),
"How to treat log values as being measured in the centre of "
"the time, or beginning (left) boundary");
declareProperty("PulseFilter", false,
"Optional. Filter out a notch of time for each entry in the "
"sample log named.\n"
"A notch of width 2*TimeTolerance is centered at each log "
"time. The value of the log is NOT used."
"This is used, for example, to filter out veto pulses.");
}
std::map<std::string, std::string> FilterByLogValue::validateInputs() {
std::map<std::string, std::string> errors;
// check for null pointers - this is to protect against workspace groups
EventWorkspace_const_sptr inputWS = this->getProperty("InputWorkspace");
if (!inputWS) {
return errors;
}
// Check that the log exists for the given input workspace
std::string logname = getPropertyValue("LogName");
try {
auto *log = dynamic_cast<ITimeSeriesProperty *>(inputWS->run().getLogData(logname));
if (log == nullptr) {
errors["LogName"] = "'" + logname + "' is not a time-series log.";
return errors;
}
} catch (Exception::NotFoundError &) {
errors["LogName"] = "The log '" + logname + "' does not exist in the workspace '" + inputWS->getName() + "'.";
return errors;
}
const double min = getProperty("MinimumValue");
const double max = getProperty("MaximumValue");
if (!isEmpty(min) && !isEmpty(max) && (max < min)) {
errors["MinimumValue"] = "MinimumValue must not be larger than MaximumValue";
errors["MaximumValue"] = "MinimumValue must not be larger than MaximumValue";
}
return errors;
}
/** Executes the algorithm
*/
void FilterByLogValue::exec() {
// convert the input workspace into the event workspace we already know it is
EventWorkspace_sptr inputWS = this->getProperty("InputWorkspace");
// Get the properties.
double min = getProperty("MinimumValue");
double max = getProperty("MaximumValue");
const double tolerance = getProperty("TimeTolerance");
const std::string logname = getPropertyValue("LogName");
const bool PulseFilter = getProperty("PulseFilter");
// Find the start and stop times of the run, but handle it if they are not
// found.
DateAndTime run_start(0), run_stop("2100-01-01T00:00:00");
bool handle_edge_values = false;
try {
run_start = inputWS->getFirstPulseTime() - tolerance;
run_stop = inputWS->getLastPulseTime() + tolerance;
handle_edge_values = true;
} catch (Exception::NotFoundError &) {
}
// Now make the splitter vector
TimeSplitterType splitter;
// This'll throw an exception if the log doesn't exist. That is good.
auto *log = dynamic_cast<ITimeSeriesProperty *>(inputWS->run().getLogData(logname));
if (log) {
if (PulseFilter) {
// ----- Filter at pulse times only -----
DateAndTime lastTime = run_start;
std::vector<DateAndTime> times = log->timesAsVector();
std::vector<DateAndTime>::iterator it;
for (it = times.begin(); it != times.end(); ++it) {
SplittingInterval interval(lastTime, *it - tolerance, 0);
// Leave a gap +- tolerance
lastTime = (*it + tolerance);
splitter.emplace_back(interval);
}
// And the last one
splitter.emplace_back(lastTime, run_stop, 0);
} else {
// ----- Filter by value ------
// This function creates the splitter vector we will use to filter out
// stuff.
const std::string logBoundary(this->getPropertyValue("LogBoundary"));
log->makeFilterByValue(splitter, min, max, tolerance, (logBoundary == CENTRE));
if (log->realSize() >= 1 && handle_edge_values) {
log->expandFilterToRange(splitter, min, max, TimeInterval(run_start, run_stop));
}
} // (filter by value)
}
g_log.information() << splitter.size() << " entries in the filter.\n";
size_t numberOfSpectra = inputWS->getNumberHistograms();
// Initialise the progress reporting object
Progress prog(this, 0.0, 1.0, numberOfSpectra);
EventWorkspace_sptr outputWS = getProperty("OutputWorkspace");
if (inputWS == outputWS) {
// Filtering in place!
// -------------------------------------------------------------
PARALLEL_FOR_NO_WSP_CHECK()
for (int64_t i = 0; i < int64_t(numberOfSpectra); ++i) {
PARALLEL_START_INTERUPT_REGION
// this is the input event list
EventList &input_el = inputWS->getSpectrum(i);
// Perform the filtering in place.
input_el.filterInPlace(splitter);
prog.report();
PARALLEL_END_INTERUPT_REGION
}
PARALLEL_CHECK_INTERUPT_REGION
// To split/filter the runs, first you make a vector with just the one
// output run
auto newRun = Kernel::make_cow<Run>(inputWS->run());
std::vector<LogManager *> splitRuns = {&newRun.access()};
inputWS->run().splitByTime(splitter, splitRuns);
// Set the output back in the input
inputWS->setSharedRun(newRun);
inputWS->mutableRun().integrateProtonCharge();
// Cast the outputWS to the matrixOutputWS and save it
this->setProperty("OutputWorkspace", inputWS);
} else {
// Make a brand new EventWorkspace for the output
outputWS = create<EventWorkspace>(*inputWS);
// Loop over the histograms (detector spectra)
PARALLEL_FOR_NO_WSP_CHECK()
for (int64_t i = 0; i < int64_t(numberOfSpectra); ++i) {
PARALLEL_START_INTERUPT_REGION
// Get the output event list (should be empty)
std::vector<EventList *> outputs{&outputWS->getSpectrum(i)};
// and this is the input event list
const EventList &input_el = inputWS->getSpectrum(i);
// Perform the filtering (using the splitting function and just one
// output)
input_el.splitByTime(splitter, outputs);
prog.report();
PARALLEL_END_INTERUPT_REGION
}
PARALLEL_CHECK_INTERUPT_REGION
// To split/filter the runs, first you make a vector with just the one
// output run
std::vector<LogManager *> output_runs;
output_runs.emplace_back(&outputWS->mutableRun());
inputWS->run().splitByTime(splitter, output_runs);
// Cast the outputWS to the matrixOutputWS and save it
this->setProperty("OutputWorkspace", outputWS);
}
}
} // namespace Algorithms
} // namespace Mantid