Newer
Older
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
/*WIKI*
TODO: Enter a full wiki-markup description of your algorithm here. You can then use the Build/wiki_maker.py script to generate your full wiki page.
*WIKI*/
#include "MantidAlgorithms/AverageLogData.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()
{
}
//----------------------------------------------------------------------------------------------
/// 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";}
//----------------------------------------------------------------------------------------------
/// Sets documentation strings for this algorithm
void AverageLogData::initDocs()
{
std::string summary = "Computes the proton charge averaged value of a given log.";
this->setWikiSummary(summary);
this->setOptionalMessage(summary);
}
//----------------------------------------------------------------------------------------------
/** Initialize the algorithm's properties.
*/
void AverageLogData::init()
{
declareProperty(new 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 checked, 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.");
}
Kernel::TimeSeriesProperty<double> * 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.push_back(EMPTY_DBL());
svalue.push_back(0.0);
pctime.push_back(EMPTY_DBL()*1.1);//larger than stime
pcvalue.push_back(0.0);
std::vector<double>::iterator istime=stime.begin(),isvalue=svalue.begin(), ipctime=pctime.begin(),ipcvalue=pcvalue.begin();
for(;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"<<std::endl;
}
g_log.debug()<<"Sum = "<<average<<std::endl<<"Sum squares = "<<error<<std::endl<<"PC = "<<protoncharge<<std::endl;
average/=protoncharge;
error/=protoncharge;
error=std::sqrt(std::fabs(error-average*average));
g_log.information()<<"Average value of "<<logname<<" is "<<average<<" +/- "<<error<<std::endl;
setProperty("Average",average);
setProperty("Error",error);
}
} // namespace Algorithms
} // namespace Mantid