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
148
149
/*WIKI*
Separates background from signal for spectra of a workspace.
*WIKI*/
#include "MantidAlgorithms/SeparateBackgroundFromSignal.h"
#include "MantidAPI/WorkspaceProperty.h"
#include "MantidAPI/MatrixWorkspace.h"
#include "MantidAPI/WorkspaceFactory.h"
#include "MantidKernel/ArrayProperty.h"
#include "MantidKernel/Statistics.h"
#include "MantidDataObjects/Workspace2D.h"
#include <sstream>
using namespace Mantid;
using namespace Mantid::API;
using namespace Mantid::Kernel;
using namespace Mantid::DataObjects;
using namespace std;
namespace Mantid
{
namespace Algorithms
{
DECLARE_ALGORITHM(SeparateBackgroundFromSignal)
//----------------------------------------------------------------------------------------------
/** Constructor
*/
SeparateBackgroundFromSignal::SeparateBackgroundFromSignal()
{
}
//----------------------------------------------------------------------------------------------
/** Destructor
*/
SeparateBackgroundFromSignal::~SeparateBackgroundFromSignal()
{
}
//----------------------------------------------------------------------------------------------
/** WIKI:
*/
void SeparateBackgroundFromSignal::initDocs()
{
setWikiSummary("Separates background from signal for spectra of a workspace.");
setOptionalMessage("Separates background from signal for spectra of a workspace.");
}
//----------------------------------------------------------------------------------------------
/** Define properties
*/
void SeparateBackgroundFromSignal::init()
{
auto inwsprop = new WorkspaceProperty<MatrixWorkspace>("InputWorkspace", "Anonymous", Direction::Input);
declareProperty(inwsprop, "Name of input MatrixWorkspace to have Z-score calculated.");
auto outwsprop = new WorkspaceProperty<Workspace2D>("OutputWorkspace", "", Direction::Output);
declareProperty(outwsprop, "Name of the output Workspace2D containing the Z-scores.");
declareProperty("WorkspaceIndex", EMPTY_INT(), "Index of the spectrum to have Z-score calculated. "
"Default is to calculate for all spectra.");
}
//----------------------------------------------------------------------------------------------
/** Execute body
*/
void SeparateBackgroundFromSignal::exec()
{
// 1. Get input and validate
MatrixWorkspace_const_sptr inpWS = getProperty("InputWorkspace");
int inpwsindex = getProperty("WorkspaceIndex");
bool separateall = false;
if (inpwsindex == EMPTY_INT())
{
separateall = true;
}
// 2. Generate output
size_t numspec;
if (separateall)
{
numspec = inpWS->getNumberHistograms();
}
else
{
numspec = 1;
}
size_t sizex = inpWS->readX(0).size();
size_t sizey = inpWS->readY(0).size();
Workspace2D_sptr outWS = boost::dynamic_pointer_cast<Workspace2D>(WorkspaceFactory::Instance().create(
"Workspace2D", numspec, sizex, sizey));
// 3. Get Z values
for (size_t i = 0; i < numspec; ++i)
{
// a) figure out wsindex
size_t wsindex;
if (separateall)
{
// Update wsindex to index in input workspace
wsindex = i;
}
else
{
// Use the wsindex as the input
wsindex = static_cast<size_t>(inpwsindex);
if (wsindex >= inpWS->getNumberHistograms())
{
stringstream errmsg;
errmsg << "Input workspace index " << inpwsindex << " is out of input workspace range = "
<< inpWS->getNumberHistograms() << endl;
}
}
// Find background
const MantidVec& inpX = inpWS->readX(wsindex);
const MantidVec& inpY = inpWS->readY(wsindex);
const MantidVec& inpE = inpWS->readE(wsindex);
MantidVec& vecX = outWS->dataX(i);
MantidVec& vecY = outWS->dataY(i);
MantidVec& vecE = outWS->dataE(i);
double Ymean, Yvariance, Ysigma;
MantidVec maskedY = inpWS->readY(wsindex);
size_t n = maskedY.size();
MantidVec mask(n,0.0);
double xn = static_cast<double>(n);
double k = 1.0;
do
{
Statistics stats = getStatistics(maskedY);
Ymean = stats.mean;
Yvariance = stats.standard_deviation * stats.standard_deviation;
Ysigma = std::sqrt((moment(maskedY,n,Ymean,4)-(xn-3.0)/(xn-1.0) * moment(maskedY,n,Ymean,2))/xn);
MantidVec::const_iterator it = std::max_element(maskedY.begin(), maskedY.end());
const size_t pos = it - maskedY.begin();
maskedY[pos] = 0;
mask[pos] = 1.0;
}
while (std::abs(Ymean-Yvariance) > k * Ysigma);
if(n > 5)
{
// remove single outliers
if (mask[1] == mask[2] && mask[2] == mask[3])
mask[0] = mask[1];
if (mask[0] == mask[2] && mask[2] == mask[3])
mask[1] = mask[2];
for (size_t l = 2; l < n-3; ++l)
if (mask[l-1] == mask[l+1] && (mask[l-1] == mask[l-2] || mask[l+1] == mask[l+2]))
{
mask[l] = mask[l+1];
}
if (mask[n-2] == mask[n-3] && mask[n-3] == mask[n-4])
mask[n-1] = mask[n-2];
if (mask[n-1] == mask[n-3] && mask[n-3] == mask[n-4])
mask[n-2] = mask[n-1];
// mask regions not connected to largest region
// for loop can start > 1 for multiple peaks
vector<cont_peak> peaks;
for (size_t l = 1; l < n; ++l)
if (mask[l] != mask[l-1] && mask[l] == 1)
{
peaks.push_back(cont_peak());
peaks[peaks.size()-1].start = l;
}
if (peaks.size() > 0)
{
size_t ipeak = peaks.size()-1;
if (mask[l] != mask[l-1] && mask[l] == 0)
{
peaks[ipeak].stop = l-1;
}
if (inpY[l] > peaks[ipeak].maxY) peaks[ipeak].maxY = inpY[l];
}
if(peaks[peaks.size()-1].stop == 0) peaks[peaks.size()-1].stop = n-1;
std::sort(peaks.begin(), peaks.end(), by_len());
for (size_t l = 1; l < peaks.size(); ++l)
{
for (size_t j = peaks[l].start; j <= peaks[l].stop; ++j)
{
mask[j] = 0;
}
}
// save output of mask * Y
vecX = inpX;
std::transform(mask.begin(), mask.end(), inpY.begin(), vecY.begin(), std::multiplies<double>());
std::transform(mask.begin(), mask.end(), inpE.begin(), vecE.begin(), std::multiplies<double>());
} // ENDFOR
// 4. Set the output
setProperty("OutputWorkspace", outWS);
return;
}
double SeparateBackgroundFromSignal::moment(MantidVec& X, size_t n, double mean, int k)