-
Notifications
You must be signed in to change notification settings - Fork 122
/
AverageLogData.cpp
129 lines (111 loc) · 4.97 KB
/
AverageLogData.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
// 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/AverageLogData.h"
#include "MantidAPI/MatrixWorkspace.h"
#include "MantidAPI/Run.h"
#include "MantidKernel/TimeSeriesProperty.h"
using namespace Mantid::Kernel;
using namespace Mantid::API;
namespace Mantid {
namespace Algorithms {
// Register the algorithm into the AlgorithmFactory
DECLARE_ALGORITHM(AverageLogData)
//----------------------------------------------------------------------------------------------
/** Constructor
*/
AverageLogData::AverageLogData() {}
//----------------------------------------------------------------------------------------------
/** Destructor
*/
AverageLogData::~AverageLogData() = default;
//----------------------------------------------------------------------------------------------
/// Algorithm's name for identification. @see Algorithm::name
const std::string AverageLogData::name() const { return "AverageLogData"; }
/// Algorithm's version for identification. @see Algorithm::version
int AverageLogData::version() const { return 1; }
/// Algorithm's category for identification. @see Algorithm::category
const std::string AverageLogData::category() const { return "DataHandling\\Logs"; }
//----------------------------------------------------------------------------------------------
/** Initialize the algorithm's properties.
*/
void AverageLogData::init() {
declareProperty(std::make_unique<WorkspaceProperty<API::MatrixWorkspace>>("InputWorkspace", "", Direction::Input),
"An input workspace that contains a Sample log property, and "
"a proton charge property.");
declareProperty("LogName", "", "Name of the log to be averaged");
declareProperty("FixZero", true,
"If true, the proton charge and the log "
"value time series are assumed to start at "
"the same moment.");
declareProperty("Average", EMPTY_DBL(), "", Direction::Output);
declareProperty("Error", EMPTY_DBL(), "", Direction::Output);
}
//----------------------------------------------------------------------------------------------
/** Execute the algorithm.
*/
void AverageLogData::exec() {
MatrixWorkspace_sptr inputWS = getProperty("InputWorkspace");
std::string logname = this->getProperty("LogName");
if (logname.empty()) {
throw std::runtime_error("Failed to supply a LogName");
}
if (!inputWS->run().hasProperty(logname)) {
throw std::runtime_error("There is no property " + logname + " in the workspace.");
}
auto *slog = dynamic_cast<Kernel::TimeSeriesProperty<double> *>(inputWS->run().getLogData(logname));
if (!slog) {
throw std::runtime_error("Problem reading property " + logname);
}
Kernel::TimeSeriesProperty<double> *pclog =
dynamic_cast<Kernel::TimeSeriesProperty<double> *>(inputWS->run().getLogData("proton_charge"));
if (!pclog) {
throw std::runtime_error("Problem reading the proton charge property");
}
double average(0), error(0), protoncharge(0);
double diffSeconds = static_cast<double>((slog->firstTime() - pclog->firstTime()).total_nanoseconds()) * 1e-9;
if (getProperty("FixZero")) {
diffSeconds = 0.;
}
std::vector<double> stime = slog->timesAsVectorSeconds();
std::vector<double> svalue = slog->valuesAsVector();
std::vector<double> pctime = pclog->timesAsVectorSeconds();
std::vector<double> pcvalue = pclog->valuesAsVector();
stime.emplace_back(EMPTY_DBL());
svalue.emplace_back(0.0);
pctime.emplace_back(EMPTY_DBL() * 1.1); // larger than stime
pcvalue.emplace_back(0.0);
auto isvalue = svalue.begin(), ipctime = pctime.begin(), ipcvalue = pcvalue.begin();
for (auto istime = stime.begin(); istime < (--stime.end()); ++istime) {
// ignore all proton pulses before the lowest time for the log
while ((*ipctime) < (*istime) + diffSeconds) {
++ipctime;
++ipcvalue;
}
// add together proton pulses before the current log time and the next log
// time
while ((*ipctime) < (*(istime + 1)) + diffSeconds) {
protoncharge += (*ipcvalue);
average += (*ipcvalue) * (*isvalue);
error += (*ipcvalue) * (*isvalue) * (*isvalue);
++ipctime;
++ipcvalue;
}
++isvalue;
}
if (protoncharge != 0) {
g_log.warning() << "Proton charge is 0. Average and standard deviations are NANs\n";
}
g_log.debug() << "Sum = " << average << "\nSum squares = " << error << "\nPC = " << protoncharge << '\n';
average /= protoncharge;
error /= protoncharge;
error = std::sqrt(std::fabs(error - average * average));
g_log.information() << "Average value of " << logname << " is " << average << " +/- " << error << '\n';
setProperty("Average", average);
setProperty("Error", error);
}
} // namespace Algorithms
} // namespace Mantid