Mantid
Loading...
Searching...
No Matches
SumEventsByLogValue.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 +
8
9#include "MantidAPI/Axis.h"
10#include "MantidAPI/Column.h"
12#include "MantidAPI/Run.h"
20
21#include <numeric>
22
23namespace Mantid::Algorithms {
24// Register the class into the algorithm factory
25DECLARE_ALGORITHM(SumEventsByLogValue)
26
27using namespace Kernel;
28using namespace API;
29
32 std::make_unique<WorkspaceProperty<DataObjects::EventWorkspace>>("InputWorkspace", "", Direction::Input),
33 "The input EventWorkspace. Must contain 'raw' (unweighted) events");
35 "MonitorWorkspace", "", Direction::Input, PropertyMode::Optional),
36 "A workspace containing the monitor counts relating to the input "
37 "workspace");
38 declareProperty(std::make_unique<WorkspaceProperty<Workspace>>("OutputWorkspace", "", Direction::Output),
39 "The name of the workspace to be created as the output of "
40 "the algorithm. The output workspace will be a "
41 "[[TableWorkspace]] in the case that a log holding integer "
42 "values is given, and a single-spectrum [[Workspace2D]] "
43 "otherwise.");
44
45 declareProperty("LogName", "", std::make_shared<MandatoryValidator<std::string>>(),
46 "The name of the number series log against which the data "
47 "should be summed");
49 std::make_unique<ArrayProperty<double>>("OutputBinning", "", std::make_shared<RebinParamsValidator>(true)),
50 "Binning parameters for the output workspace (see [[Rebin]] for syntax) "
51 "(Optional for logs holding integer values, mandatory otherwise)");
52}
53
54std::map<std::string, std::string> SumEventsByLogValue::validateInputs() {
55 std::map<std::string, std::string> errors;
56
57 // check for null pointers - this is to protect against workspace groups
58 m_inputWorkspace = getProperty("InputWorkspace");
59 if (!m_inputWorkspace) {
60 return errors;
61 }
62
63 // This only works for unweighted events
64 // TODO: Either turn this check into a proper validator or amend the algorithm
65 // to work for weighted events
66 if (m_inputWorkspace->getEventType() != API::TOF) {
67 errors["InputWorkspace"] = "This algorithm only works for unweighted ('raw') events";
68 }
69
70 // Check that the log exists for the given input workspace
71 m_logName = getPropertyValue("LogName");
72 try {
73 auto *log = dynamic_cast<ITimeSeriesProperty *>(m_inputWorkspace->run().getLogData(m_logName));
74 if (log == nullptr) {
75 errors["LogName"] = "'" + m_logName + "' is not a time-series log.";
76 return errors;
77 }
78 if (log->realSize() == 0) {
79 errors["LogName"] = "'" + m_logName + "' is empty.";
80 }
81 } catch (Exception::NotFoundError &) {
82 errors["LogName"] =
83 "The log '" + m_logName + "' does not exist in the workspace '" + m_inputWorkspace->getName() + "'.";
84 return errors;
85 }
86
87 return errors;
88}
89
91 // Get hold of the requested log. Will throw if it doesn't exist (which is
92 // what we want).
93 const Property *const log = m_inputWorkspace->run().getLogData(m_logName);
94
95 // Now we need to know what type of property it is
96 const auto *intLog = dynamic_cast<const TimeSeriesProperty<int> *>(log);
97 const auto *dblLog = dynamic_cast<const TimeSeriesProperty<double> *>(log);
98
99 m_binningParams = getProperty("OutputBinning");
100 // Binning parameters must be provided for floating point logs
101 if (m_binningParams.empty()) {
102 if (intLog != nullptr) {
103 createTableOutput(intLog);
104 } else {
105 throw std::invalid_argument("OutputBinning must be provided for floating-point number logs");
106 }
107 } else // Binning parameters have been given
108 {
109 if (intLog != nullptr) {
110 createBinnedOutput(intLog);
111 } else if (dblLog != nullptr) {
112 createBinnedOutput(dblLog);
113 }
114 // else if ( dynamic_cast<const TimeSeriesProperty<std::string>*>(log) !=
115 // NULL )
116 //{
117 // TODO: Implement (if anyone ever asks for it).
118 //}
119 else {
120 throw std::runtime_error("This algorithm only supports number-series logs");
121 }
122 }
123}
124
129 // This is the version for integer logs when no binning parameters have been
130 // given and has a data point per log value
131 const int minVal = log->minValue();
132 const int maxVal = log->maxValue();
133 const int xLength = maxVal - minVal + 1;
134
135 if (xLength > 10000) {
136 g_log.warning() << "Did you really want to create a " << xLength << " row table? This will take some time!\n";
137 }
138
139 // Accumulate things in a local vector before transferring to the table
140 std::vector<int> Y(xLength);
141 const auto numSpec = static_cast<int>(m_inputWorkspace->getNumberHistograms());
142 Progress prog(this, 0.0, 1.0, numSpec + xLength);
144 for (int spec = 0; spec < numSpec; ++spec) {
146 const IEventList &eventList = m_inputWorkspace->getSpectrum(spec);
147 filterEventList(eventList, minVal, maxVal, log, Y);
148 prog.report();
150 }
152
153 // Create a table workspace to hold the sum.
154 ITableWorkspace_sptr outputWorkspace = WorkspaceFactory::Instance().createTable();
155 auto logValues = outputWorkspace->addColumn("int", m_logName);
156 auto counts = outputWorkspace->addColumn("int", "Counts");
157 auto errors = outputWorkspace->addColumn("double", "Error");
158 outputWorkspace->setRowCount(xLength); // One row per log value across the full range
159 // Set type for benefit of MantidPlot
160 logValues->setPlotType(1); // X
161 counts->setPlotType(2); // Y
162 errors->setPlotType(5); // E
163
164 // Transfer the results to the table
165 for (int i = 0; i < xLength; ++i) {
166 logValues->cell<int>(i) = minVal + i;
167 counts->cell<int>(i) = Y[i];
168 errors->cell<double>(i) = std::sqrt(Y[i]);
169 }
170
171 // Columns for normalisation: monitors (if available), time & proton charge
172 addMonitorCounts(outputWorkspace, log, minVal, maxVal);
173 // Add a column to hold the time duration (in seconds) for which the log had a
174 // certain value
175 auto timeCol = outputWorkspace->addColumn("double", "time");
176 // Add a column to hold the proton charge for which the log had a certain
177 // value
178 auto protonChgCol = outputWorkspace->addColumn("double", "proton_charge");
179 // Get hold of the proton charge log for later
180 const TimeSeriesProperty<double> *protonChargeLog = nullptr;
181 try {
182 protonChargeLog = m_inputWorkspace->run().getTimeSeriesProperty<double>("proton_charge");
183 // Set back to NULL if the log is empty or bad things will happen later
184 if (protonChargeLog->realSize() == 0)
185 protonChargeLog = nullptr;
186 } catch (std::exception &) {
187 // Log and carry on if not found. Column will be left empty.
188 g_log.warning("proton_charge log not found in workspace.");
189 }
190
191 // Get a list of the other time-series logs in the input workspace
192 auto otherLogs = getNumberSeriesLogs();
193 // Add a column for each of these 'other' logs
194 for (auto &otherLog : otherLogs) {
195 auto newColumn = outputWorkspace->addColumn("double", otherLog.first);
196 // For the benefit of MantidPlot, set these columns to be containing X
197 // values
198 newColumn->setPlotType(1);
199 }
200
201 // Now to get the average value of other time-varying logs
202 // Loop through the values of the 'main' log
203 for (int value = minVal; value <= maxVal; ++value) {
204 const int row = value - minVal;
205 // Create a filter giving the times when this log has the current value
206 TimeSplitterType filter;
207 log->makeFilterByValue(filter, value,
208 value); // min & max are the same of course
209
210 // This section ensures that the filter goes to the end of the run
211 if (value == log->lastValue() && protonChargeLog) {
212 TimeInterval timeAfterLastLogValue(log->lastTime(), m_inputWorkspace->getLastPulseTime());
213 log->expandFilterToRange(filter, value, value, timeAfterLastLogValue);
214 }
215
216 // Calculate the time covered by this log value and add it to the table
217 double duration = 0.0;
218 for (auto &time : filter) {
219 duration += time.duration();
220 }
221 timeCol->cell<double>(row) = duration;
222
224 // Sum up the proton charge for this log value
225 if (protonChargeLog)
226 protonChgCol->cell<double>(row) = sumProtonCharge(protonChargeLog, filter);
228
229 for (auto &otherLog : otherLogs) {
230 // Calculate the average value of each 'other' log for the current value
231 // of the main log
232 // Have to (maybe inefficiently) fetch back column by name - move outside
233 // loop if too slow
234 outputWorkspace->getColumn(otherLog.first)->cell<double>(row) = otherLog.second->averageValueInFilter(filter);
235 }
236 prog.report();
237 }
238
239 setProperty("OutputWorkspace", outputWorkspace);
240}
241
253void SumEventsByLogValue::filterEventList(const API::IEventList &eventList, const int minVal, const int maxVal,
254 const Kernel::TimeSeriesProperty<int> *log, std::vector<int> &Y) {
255 if (log->realSize() == 0)
256 return;
257
258 const auto pulseTimes = eventList.getPulseTimes();
259 for (auto pulseTime : pulseTimes) {
260 // Find the value of the log at the time of this event
261 // This algorithm is really concerned with 'slow' logs so we don't care
262 // about
263 // the time of the event within the pulse.
264 // NB: If the pulse time is before the first log entry, we get the first
265 // value.
266 const int logValue = log->getSingleValue(pulseTime);
267
268 if (logValue >= minVal && logValue <= maxVal) {
269 // In this scenario it's easy to know what bin to increment
271 ++Y[logValue - minVal];
272 }
273 }
274}
275
285 const TimeSeriesProperty<int> *log, const int minVal, const int maxVal) {
286 DataObjects::EventWorkspace_const_sptr monitorWorkspace = getProperty("MonitorWorkspace");
287 // If no monitor workspace was given, there's nothing to do
288 if (!monitorWorkspace)
289 return;
290
291 const auto &spectrumInfo = monitorWorkspace->spectrumInfo();
292
293 const int xLength = maxVal - minVal + 1;
294 // Loop over the spectra - there will be one per monitor
295 for (std::size_t spec = 0; spec < monitorWorkspace->getNumberHistograms(); ++spec) {
296 try {
297 // Create a column for this monitor
298 const std::string monitorName = spectrumInfo.detector(spec).getName();
299 auto monitorCounts = outputWorkspace->addColumn("int", monitorName);
300 const IEventList &eventList = monitorWorkspace->getSpectrum(spec);
301 // Accumulate things in a local vector before transferring to the table
302 // workspace
303 std::vector<int> Y(xLength);
304 filterEventList(eventList, minVal, maxVal, log, Y);
305 // Transfer the results to the table
306 for (int i = 0; i < xLength; ++i) {
307 monitorCounts->cell<int>(i) = Y[i];
308 }
309 } catch (Exception::NotFoundError &) {
310 // ADARA-generated nexus files have sometimes been seen to contain
311 // 'phantom' monitors that aren't in the IDF.
312 // This handles that by ignoring those spectra.
313 }
314 }
315}
316
323std::vector<std::pair<std::string, const Kernel::ITimeSeriesProperty *>> SumEventsByLogValue::getNumberSeriesLogs() {
324 std::vector<std::pair<std::string, const Kernel::ITimeSeriesProperty *>> numberSeriesProps;
325 const auto &logs = m_inputWorkspace->run().getLogData();
326 for (auto log : logs) {
327 const std::string logName = log->name();
328 // Don't add the log that's the one being summed against
329 if (logName == m_logName)
330 continue;
331 // Exclude the proton charge log as we have a separate column for the sum of
332 // that
333 if (logName == "proton_charge")
334 continue;
335 // Try to cast to an ITimeSeriesProperty
336 auto tsp = dynamic_cast<const ITimeSeriesProperty *>(log);
337 // Move on to the next one if this is not a TSP
338 if (tsp == nullptr)
339 continue;
340 // Don't keep ones with only one entry
341 // if ( tsp->realSize() < 2 ) continue;
342 // Now make sure it's either an int or double tsp, and if so add log to the
343 // list
344 if (dynamic_cast<TimeSeriesProperty<double> *>(log) || dynamic_cast<TimeSeriesProperty<int> *>(log)) {
345 numberSeriesProps.emplace_back(logName, tsp);
346 }
347 }
348
349 return numberSeriesProps;
350}
351
358 const Kernel::TimeSplitterType &filter) {
359 // Clone the proton charge log and filter the clone on this log value
360 std::unique_ptr<Kernel::TimeSeriesProperty<double>> protonChargeLogClone(protonChargeLog->clone());
361 protonChargeLogClone->filterByTimes(filter);
362 // Seems like the only way to sum this is to yank out the values
363 const std::vector<double> pcValues = protonChargeLogClone->valuesAsVector();
364 return std::accumulate(pcValues.begin(), pcValues.end(), 0.0);
365}
366
372 // If only the number of bins was given, add the min & max values of the log
373 if (m_binningParams.size() == 1) {
374 m_binningParams.insert(m_binningParams.begin(), log->minValue());
375 m_binningParams.emplace_back(log->maxValue() * 1.000001); // Make it a tiny bit larger to cover full range
376 }
377
378 // XValues will be resized in createAxisFromRebinParams()
379 std::vector<double> XValues;
380 const int XLength = VectorHelper::createAxisFromRebinParams(m_binningParams, XValues);
381 assert((int)XValues.size() == XLength);
382
383 // Create the output workspace - the factory will give back a Workspace2D
384 MatrixWorkspace_sptr outputWorkspace = WorkspaceFactory::Instance().create("Workspace2D", 1, XLength, XLength - 1);
385 // Copy the bin boundaries into the output workspace
386 outputWorkspace->mutableX(0) = XValues;
387 outputWorkspace->getAxis(0)->title() = m_logName;
388 outputWorkspace->setYUnit("Counts");
389
390 auto &Y = outputWorkspace->mutableY(0);
391 const auto numSpec = static_cast<int>(m_inputWorkspace->getNumberHistograms());
392 Progress prog(this, 0.0, 1.0, numSpec);
394 for (int spec = 0; spec < numSpec; ++spec) {
396 const IEventList &eventList = m_inputWorkspace->getSpectrum(spec);
397 const auto pulseTimes = eventList.getPulseTimes();
398 for (auto pulseTime : pulseTimes) {
399 // Find the value of the log at the time of this event
400 const double logValue = log->getSingleValue(pulseTime);
401 if (logValue >= XValues.front() && logValue < XValues.back()) {
403 ++Y[VectorHelper::getBinIndex(XValues, logValue)];
404 }
405 }
406
407 prog.report();
409 }
411
412 // The errors are the sqrt of the counts so long as we don't deal with
413 // weighted events.
414 std::transform(Y.cbegin(), Y.cend(), outputWorkspace->mutableE(0).begin(), (double (*)(double))std::sqrt);
415
416 setProperty("OutputWorkspace", outputWorkspace);
417}
418
419} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
double value
The value of the point.
Definition: FitMW.cpp:51
#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_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_ATOMIC
#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
std::string getPropertyValue(const std::string &name) const override
Get the value of a property as a string.
Definition: Algorithm.cpp:2026
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
Definition: Algorithm.cpp:2076
Kernel::Logger & g_log
Definition: Algorithm.h:451
void interruption_point()
This is called during long-running operations, and check if the algorithm has requested that it be ca...
Definition: Algorithm.cpp:1687
IEventList : Interface to Mantid::DataObjects::EventList class, used to expose to PythonAPI.
Definition: IEventList.h:26
virtual std::vector< Mantid::Types::Core::DateAndTime > getPulseTimes() const =0
Return the list of pulse time values.
Helper class for reporting progress from algorithms.
Definition: Progress.h:25
A property class for workspaces.
void init() override
Virtual method - must be overridden by concrete algorithm.
std::vector< double > m_binningParams
The optional binning parameters.
void addMonitorCounts(const API::ITableWorkspace_sptr &outputWorkspace, const Kernel::TimeSeriesProperty< int > *log, const int minVal, const int maxVal)
Looks for monitor event data and, if found, adds columns to the output table corresponding to the mon...
void filterEventList(const API::IEventList &eventList, const int minVal, const int maxVal, const Kernel::TimeSeriesProperty< int > *log, std::vector< int > &Y)
Goes through an event list assigning events to the output vector according to the log value at the ti...
void exec() override
Virtual method - must be overridden by concrete algorithm.
double sumProtonCharge(const Kernel::TimeSeriesProperty< double > *protonChargeLog, const Kernel::TimeSplitterType &filter)
Integrates the proton charge between specified times.
std::map< std::string, std::string > validateInputs() override
Method checking errors on ALL the inputs, before execution.
void createBinnedOutput(const Kernel::TimeSeriesProperty< T > *log)
Create a single-spectrum Workspace2D containing the integrated counts versus log value,...
void createTableOutput(const Kernel::TimeSeriesProperty< int > *log)
Produces the table workspace output for an integer TimeSeriesProperty.
std::string m_logName
The name of the log to sum against.
DataObjects::EventWorkspace_const_sptr m_inputWorkspace
The input workspace.
std::vector< std::pair< std::string, const Kernel::ITimeSeriesProperty * > > getNumberSeriesLogs()
Searches the input workspace for int or double TimeSeriesProperty's other than the one being summed a...
Support for a property that holds an array of values.
Definition: ArrayProperty.h:28
Exception for when an item is not found in a collection.
Definition: Exception.h:145
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
A non-templated interface to a TimeSeriesProperty.
void warning(const std::string &msg)
Logs at warning level.
Definition: Logger.cpp:86
Validator to check that a property is not left empty.
void report()
Increments the loop counter by 1, then sends the progress notification on behalf of its algorithm.
Definition: ProgressBase.h:51
Base class for properties.
Definition: Property.h:94
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
Represents a time interval.
Definition: DateAndTime.h:25
A specialised Property class for holding a series of time-value pairs.
TimeSeriesProperty< TYPE > * clone() const override
"Virtual" copy constructor
TYPE minValue() const
Returns the minimum value found in the series.
TYPE maxValue() const
Returns the maximum value found in the series.
void makeFilterByValue(std::vector< SplittingInterval > &split, double min, double max, double TimeTolerance=0.0, bool centre=false) const override
Fill a TimeSplitterType that will filter the events by matching.
void expandFilterToRange(std::vector< SplittingInterval > &split, double min, double max, const TimeInterval &range) const override
Make sure an existing filter covers the full time range given.
Types::Core::DateAndTime lastTime() const
Returns the last time.
TYPE lastValue() const
Returns the last value.
TYPE getSingleValue(const Types::Core::DateAndTime &t) const
Returns the value at a particular time.
int realSize() const override
Returns the real size of the time series property map:
std::shared_ptr< ITableWorkspace > ITableWorkspace_sptr
shared pointer to Mantid::API::ITableWorkspace
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
std::shared_ptr< const EventWorkspace > EventWorkspace_const_sptr
shared pointer to a const Workspace2D
MANTID_KERNEL_DLL int getBinIndex(const std::vector< double > &bins, const double value)
Return the index into a vector of bin boundaries for a particular X value.
int MANTID_KERNEL_DLL createAxisFromRebinParams(const std::vector< double > &params, std::vector< double > &xnew, const bool resize_xnew=true, const bool full_bins_only=false, const double xMinHint=std::nan(""), const double xMaxHint=std::nan(""), const bool useReverseLogarithmic=false, const double power=-1)
Creates a new output X array given a 'standard' set of rebinning parameters.
std::vector< SplittingInterval > TimeSplitterType
A typedef for splitting events according their pulse time.
Definition: LogManager.h:31
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
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54