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
/*WIKI*
This algorithm calculates the weighted mean from all the spectra in a given
workspace. Monitors and masked spectra are ignored. Also, individual bins with
IEEE values will be excluded from the result. The weighted mean calculated by
the following:
<math>\displaystyle y=\frac{\sum\frac{x_i}{\sigma^{2}_i}}{\sum\frac{1}{\sigma^{2}_i}} </math>
and the variance is calculated by:
<math>\displaystyle \sigma^{2}_y=\frac{1}{\sum\frac{1}{\sigma^{2}_i}} </math>
*WIKI*/
#include "MantidAlgorithms/WeightedMeanOfWorkspace.h"
#include "MantidDataObjects/EventWorkspace.h"
namespace Mantid
{
using namespace API;
using namespace DataObjects;
using namespace Geometry;
using namespace Kernel;
namespace Algorithms
{
// Register the algorithm into the AlgorithmFactory
DECLARE_ALGORITHM(WeightedMeanOfWorkspace)
//----------------------------------------------------------------------------------------------
/** Constructor
*/
WeightedMeanOfWorkspace::WeightedMeanOfWorkspace()
{
}
//----------------------------------------------------------------------------------------------
/** Destructor
*/
WeightedMeanOfWorkspace::~WeightedMeanOfWorkspace()
{
}
//----------------------------------------------------------------------------------------------
/// Algorithm's name for identification. @see Algorithm::name
const std::string WeightedMeanOfWorkspace::name() const { return "WeightedMeanOfWorkspace"; };
/// Algorithm's version for identification. @see Algorithm::version
int WeightedMeanOfWorkspace::version() const { return 1; };
/// Algorithm's category for identification. @see Algorithm::category
const std::string WeightedMeanOfWorkspace::category() const { return "Arithmetic"; }
//----------------------------------------------------------------------------------------------
/// Sets documentation strings for this algorithm
void WeightedMeanOfWorkspace::initDocs()
{
this->setWikiSummary("This algorithm calculates the weighted mean for an entire workspace.");
this->setOptionalMessage("This algorithm calculates the weighted mean for an entire workspace.");
}
//----------------------------------------------------------------------------------------------
/** Initialize the algorithm's properties.
*/
void WeightedMeanOfWorkspace::init()
{
this->declareProperty(new WorkspaceProperty<>("InputWorkspace", "",
Direction::Input), "An input workspace.");
this->declareProperty(new WorkspaceProperty<>("OutputWorkspace", "",
Direction::Output), "An output workspace.");
}
//----------------------------------------------------------------------------------------------
/** Execute the algorithm.
*/
void WeightedMeanOfWorkspace::exec()
{
MatrixWorkspace_sptr inputWS = this->getProperty("InputWorkspace");
g_log.warning() << "Help1" << std::endl;
// Check if it is an event workspace
EventWorkspace_const_sptr eventW = boost::dynamic_pointer_cast<const EventWorkspace>(inputWS);
if (eventW != NULL)
{
throw std::runtime_error("WeightedMeanOfWorkspace cannot handle EventWorkspaces!");
}
g_log.warning() << "Help2" << std::endl;
// Create the output workspace
MatrixWorkspace_sptr singleValued = WorkspaceFactory::Instance().create("WorkspaceSingleValue", 1, 1, 1);
g_log.warning() << "Help3" << std::endl;
// Calculate weighted mean
std::size_t numHists = inputWS->getNumberHistograms();
double averageValue = 0.0;
double weightSum = 0.0;
for (std::size_t i = 0; i < numHists; ++i)
{
g_log.warning() << "Iteration " << i << std::endl;
IDetector_const_sptr det = inputWS->getDetector(i);
if( det->isMonitor() || det->isMasked() )
{
continue;
}
g_log.warning() << "OK Pixel" << std::endl;
MantidVec y = inputWS->dataY(i);
MantidVec e = inputWS->dataE(i);
double weight = 0.0;
g_log.warning() << "A: " << y.size() << std::endl;
for (std::size_t j = 0; j < y.size(); ++j)
{
if (!boost::math::isnan(y[j]) && !boost::math::isinf(y[j]) &&
!boost::math::isnan(e[j]) && !boost::math::isinf(e[j]))
{
weight = 1.0 / (e[j] * e[j]);
g_log.warning() << "B: " << weight << std::endl;
averageValue += (y[j] * weight);
g_log.warning() << "G: " << averageValue << std::endl;
weightSum += weight;
}
}
}
g_log.warning() << "Help4 " << averageValue << ", " << weightSum << std::endl;
singleValued->dataX(0)[0] = 0.0;
singleValued->dataY(0)[0] = averageValue / weightSum;
singleValued->dataE(0)[0] = std::sqrt(weightSum);
g_log.warning() << "Help5" << std::endl;
this->setProperty("OutputWorkspace", singleValued);
g_log.warning() << "Help6" << std::endl;
}
} // namespace Algorithms
} // namespace Mantid