Mantid
Loading...
Searching...
No Matches
FindPeakBackground.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 +
14#include "MantidHistogramData/EstimatePolynomial.h"
18
19#include <sstream>
20
21using namespace Mantid;
22using namespace Mantid::API;
23using namespace Mantid::Kernel;
24using namespace Mantid::DataObjects;
25using namespace std;
26
27namespace Mantid::Algorithms {
28
29DECLARE_ALGORITHM(FindPeakBackground)
30
31//----------------------------------------------------------------------------------------------
34void FindPeakBackground::init() {
35 declareProperty(std::make_unique<WorkspaceProperty<MatrixWorkspace>>("InputWorkspace", "Anonymous", Direction::Input),
36 "Name of input MatrixWorkspace that contains peaks.");
37
38 declareProperty("WorkspaceIndex", EMPTY_INT(),
39 "workspace indices to have peak and background separated. "
40 "No default is taken. ");
41
42 declareProperty("SigmaConstant", 1.0,
43 "Multiplier of standard deviations of the variance for convergence of "
44 "peak elimination. Default is 1.0. ");
45
46 declareProperty(std::make_unique<ArrayProperty<double>>("FitWindow"),
47 "Optional: enter a comma-separated list of the minimum and "
48 "maximum X-positions of window to fit. "
49 "The window is the same for all indices in workspace. The "
50 "length must be exactly two.");
51
52 std::vector<std::string> bkgdtypes{"Flat", "Linear", "Quadratic"};
53 declareProperty("BackgroundType", "Linear", std::make_shared<StringListValidator>(bkgdtypes), "Type of Background.");
54
55 // The found peak in a table
56 declareProperty(std::make_unique<WorkspaceProperty<API::ITableWorkspace>>("OutputWorkspace", "", Direction::Output),
57 "The name of the TableWorkspace in which to store the background found "
58 "for each index. "
59 "Table contains the indices of the beginning and ending of peak "
60 "and the estimated background coefficients for the constant, linear, and "
61 "quadratic terms.");
62}
63
64void FindPeakBackground::findWindowIndex(const HistogramData::Histogram &histogram, size_t &l0, size_t &n) {
65 const auto &inpX = histogram.x();
66 const auto &inpY = histogram.y();
67 const size_t sizey = inpY.size(); // inpWS->y(inpwsindex).size();
68
69 // determine the fit window with their index in X (or Y)
70 n = sizey;
71 l0 = 0;
72 if (m_vecFitWindows.size() > 1) {
74 l0 = fp.getIndex(inpX, m_vecFitWindows[0]);
75 n = fp.getIndex(inpX, m_vecFitWindows[1]);
76 if (n < sizey)
77 n++;
78 }
79}
80
81//----------------------------------------------------------------------------------------------
85 // Get input and validate
87 auto histogram = m_inputWS->histogram(m_inputWSIndex);
88
89 size_t l0, n;
90 findWindowIndex(histogram, l0, n);
91
92 // m_vecFitWindows won't be used again form this point till end.
93
94 // Set up output table workspace
96
97 // 3. Get Y values
98 Progress prog(this, 0.0, 1.0, 1);
99
100 std::vector<size_t> peak_min_max_indexes;
101 std::vector<double> bkgd3;
102 int goodfit = findBackground(histogram, l0, n, peak_min_max_indexes, bkgd3);
103
104 if (goodfit > 0) {
105 size_t min_peak = peak_min_max_indexes[0];
106 size_t max_peak = peak_min_max_indexes[1];
107 double a0 = bkgd3[0];
108 double a1 = bkgd3[1];
109 double a2 = bkgd3[2];
110 API::TableRow t = m_outPeakTableWS->getRow(0);
111 t << static_cast<int>(m_inputWSIndex) << static_cast<int>(min_peak) << static_cast<int>(max_peak) << a0 << a1 << a2
112 << goodfit;
113 }
114
115 prog.report();
116
117 // 4. Set the output
118 setProperty("OutputWorkspace", m_outPeakTableWS);
119}
120
121//----------------------------------------------------------------------------------------------
131int FindPeakBackground::findBackground(const HistogramData::Histogram &histogram, const size_t &l0, const size_t &n,
132 std::vector<size_t> &peak_min_max_indexes, std::vector<double> &bkgd3) {
133 const size_t sizex = histogram.x().size();
134 const auto &inpY = histogram.y();
135 const size_t sizey = inpY.size();
136
137 int goodfit(0);
138
139 // Find background
140 double Ymean, Yvariance, Ysigma;
141 MantidVec maskedY;
142 auto in = std::min_element(inpY.cbegin(), inpY.cend());
143 double bkg0 = inpY[in - inpY.begin()];
144 for (size_t l = l0; l < n; ++l) {
145 maskedY.emplace_back(inpY[l] - bkg0);
146 }
147 MantidVec mask(n - l0, 0.0);
148 auto xn = static_cast<double>(n - l0);
149 if ((0. == xn) || (0. == xn - 1.0))
150 throw std::runtime_error("The number of Y values in the input workspace for the "
151 "workspace index given, minus 'l0' or minus 'l0' minus 1, is 0. This "
152 "will produce a "
153 "divide-by-zero");
154 do {
155 Statistics stats = getStatistics(maskedY);
156 Ymean = stats.mean;
157 Yvariance = stats.standard_deviation * stats.standard_deviation;
158 Ysigma = std::sqrt((moment4(maskedY, static_cast<size_t>(xn), Ymean) - (xn - 3.0) / (xn - 1.0) * Yvariance) / xn);
159 MantidVec::const_iterator it = std::max_element(maskedY.begin(), maskedY.end());
160 const size_t pos = it - maskedY.begin();
161 maskedY[pos] = 0;
162 mask[pos] = 1.0;
163 } while (std::abs(Ymean - Yvariance) > m_sigmaConstant * Ysigma);
164
165 if (n - l0 > 5) {
166 // remove single outliers
167 if (mask[1] == mask[2] && mask[2] == mask[3])
168 mask[0] = mask[1];
169 if (mask[0] == mask[2] && mask[2] == mask[3])
170 mask[1] = mask[2];
171 for (size_t l = 2; l < n - l0 - 3; ++l) {
172 if (mask[l - 1] == mask[l + 1] && (mask[l - 1] == mask[l - 2] || mask[l + 1] == mask[l + 2])) {
173 mask[l] = mask[l + 1];
174 }
175 }
176 if (mask[n - l0 - 2] == mask[n - l0 - 3] && mask[n - l0 - 3] == mask[n - l0 - 4])
177 mask[n - l0 - 1] = mask[n - l0 - 2];
178 if (mask[n - l0 - 1] == mask[n - l0 - 3] && mask[n - l0 - 3] == mask[n - l0 - 4])
179 mask[n - l0 - 2] = mask[n - l0 - 1];
180
181 // mask regions not connected to largest region
182 // for loop can start > 1 for multiple peaks
183 vector<cont_peak> peaks;
184 if (mask[0] == 1) {
185 peaks.emplace_back();
186 peaks.back().start = l0;
187 }
188 for (size_t l = 1; l < n - l0; ++l) {
189 if (mask[l] != mask[l - 1] && mask[l] == 1) {
190 peaks.emplace_back();
191 peaks.back().start = l + l0;
192 } else if (!peaks.empty()) {
193 size_t ipeak = peaks.size() - 1;
194 if (mask[l] != mask[l - 1] && mask[l] == 0) {
195 peaks[ipeak].stop = l + l0;
196 }
197 if (inpY[l + l0] > peaks[ipeak].maxY)
198 peaks[ipeak].maxY = inpY[l + l0];
199 }
200 }
201 size_t min_peak, max_peak;
202 if (!peaks.empty()) {
203 g_log.debug() << "Peaks' size = " << peaks.size() << " -> esitmate background. \n";
204 if (peaks.back().stop == 0)
205 peaks.back().stop = n - 1;
206 std::sort(peaks.begin(), peaks.end(), by_len());
207
208 // save endpoints
209 min_peak = peaks[0].start;
210 // extra point for histogram input - TODO change to use Histogram better
211 max_peak = peaks[0].stop + sizex - sizey;
212 goodfit = 1;
213 } else {
214 // assume the whole thing is background
215 g_log.debug("Peaks' size = 0 -> whole region assumed background");
216 min_peak = n;
217 max_peak = l0;
218
219 goodfit = 2;
220 }
221
222 double a0 = 0., a1 = 0., a2 = 0.;
223 estimateBackground(histogram, l0, n, min_peak, max_peak, (!peaks.empty()), a0, a1, a2);
224
225 // Add a new row
226 peak_min_max_indexes.resize(2);
227 peak_min_max_indexes[0] = min_peak;
228 peak_min_max_indexes[1] = max_peak;
229
230 bkgd3.resize(3);
231 bkgd3[0] = a0;
232 bkgd3[1] = a1;
233 bkgd3[2] = a2;
234 }
235
236 return goodfit;
237}
238
239//----------------------------------------------------------------------------------------------
251void FindPeakBackground::estimateBackground(const HistogramData::Histogram &histogram, const size_t i_min,
252 const size_t i_max, const size_t p_min, const size_t p_max,
253 const bool hasPeak, double &out_bg0, double &out_bg1, double &out_bg2) {
254 double redux_chisq;
255 if (hasPeak) {
256 HistogramData::estimateBackground(m_backgroundOrder, histogram, i_min, i_max, p_min, p_max, out_bg0, out_bg1,
257 out_bg2, redux_chisq);
258 } else {
259 HistogramData::estimatePolynomial(m_backgroundOrder, histogram, i_min, i_max, out_bg0, out_bg1, out_bg2,
260 redux_chisq);
261 }
262
263 g_log.information() << "Estimated background: A0 = " << out_bg0 << ", A1 = " << out_bg1 << ", A2 = " << out_bg2
264 << "\n";
265}
266//----------------------------------------------------------------------------------------------
272double FindPeakBackground::moment4(MantidVec &X, size_t n, double mean) {
273 double sum = 0.0;
274 for (size_t i = 0; i < n; ++i) {
275 sum += (X[i] - mean) * (X[i] - mean) * (X[i] - mean) * (X[i] - mean);
276 }
277 sum /= static_cast<double>(n);
278 return sum;
279}
280
281//----------------------------------------------------------------------------------------------
283 // process input workspace and workspace index
284 m_inputWS = getProperty("InputWorkspace");
285
286 int inpwsindex = getProperty("WorkspaceIndex");
287 if (isEmpty(inpwsindex)) {
288 // Default
289 if (m_inputWS->getNumberHistograms() == 1) {
290 inpwsindex = 0;
291 } else {
292 throw runtime_error("WorkspaceIndex must be given. ");
293 }
294 } else if (inpwsindex < 0 || inpwsindex >= static_cast<int>(m_inputWS->getNumberHistograms())) {
295 stringstream errss;
296 errss << "Input workspace " << m_inputWS->getName() << " has " << m_inputWS->getNumberHistograms()
297 << " spectra. Input workspace index " << inpwsindex << " is out of boundary. ";
298 throw runtime_error(errss.str());
299 }
300 m_inputWSIndex = static_cast<size_t>(inpwsindex);
301
302 std::vector<double> fitwindow = getProperty("FitWindow");
303 setFitWindow(fitwindow);
304
305 // background
306 m_backgroundType = getPropertyValue("BackgroundType");
307 size_t bkgdorder = 0;
308 if (m_backgroundType == "Linear")
309 bkgdorder = 1;
310 else if (m_backgroundType == "Quadratic")
311 bkgdorder = 2;
312 setBackgroundOrder(bkgdorder);
313
314 // sigma constant
315 double k = getProperty("SigmaConstant");
316 setSigma(k);
317}
318
321
324
325//----------------------------------------------------------------------------------------------
330void FindPeakBackground::setFitWindow(const std::vector<double> &fitwindow) {
331 // validate input
332 if ((fitwindow.size() == 2) && fitwindow[0] >= fitwindow[1]) {
333 throw std::invalid_argument("Fit window has either wrong item number or "
334 "window value is not in ascending order.");
335 }
336
337 // m_vecFitWindows.resize(2);
338 // copy the input to class variable
339 m_vecFitWindows = fitwindow;
340}
341
342//----------------------------------------------------------------------------------------------
347 // Set up output table workspace
348 m_outPeakTableWS = std::make_shared<TableWorkspace>();
349 m_outPeakTableWS->addColumn("int", "wksp_index");
350 m_outPeakTableWS->addColumn("int", "peak_min_index");
351 m_outPeakTableWS->addColumn("int", "peak_max_index");
352 m_outPeakTableWS->addColumn("double", "bkg0");
353 m_outPeakTableWS->addColumn("double", "bkg1");
354 m_outPeakTableWS->addColumn("double", "bkg2");
355 m_outPeakTableWS->addColumn("int", "GoodFit");
356
357 m_outPeakTableWS->appendRow();
358}
359
360} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
double sigma
Definition: GetAllEi.cpp:156
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
static bool isEmpty(const NumT toCheck)
checks that the value was not set by users, uses the value in empty double/int.
Helper class for reporting progress from algorithms.
Definition: Progress.h:25
TableRow represents a row in a TableWorkspace.
Definition: TableRow.h:39
A property class for workspaces.
FindPeakBackground : Calculate Zscore for a Matrix Workspace.
void estimateBackground(const HistogramData::Histogram &histogram, const size_t i_min, const size_t i_max, const size_t p_min, const size_t p_max, const bool hasPeak, double &out_bg0, double &out_bg1, double &out_bg2)
Estimate background.
API::MatrixWorkspace_const_sptr m_inputWS
Input workspace.
void exec() override
Implement abstract Algorithm methods.
std::vector< double > m_vecFitWindows
fit window
void setFitWindow(const std::vector< double > &window)
set fit window
void findWindowIndex(const HistogramData::Histogram &histogram, size_t &l0, size_t &n)
find fit window's data point index
double moment4(MantidVec &X, size_t n, double mean)
Calculate 4th moment.
size_t m_backgroundOrder
background order: 0 for flat, 1 for linear, 2 for quadratic
void setBackgroundOrder(size_t order)
set background order
API::ITableWorkspace_sptr m_outPeakTableWS
output workspace (table of result)
int findBackground(const HistogramData::Histogram &histogram, const size_t &l0, const size_t &n, std::vector< size_t > &peak_min_max_indexes, std::vector< double > &bkgd3)
main method to calculate background
void setSigma(const double &sigma)
set sigma constant
void createOutputWorkspaces()
create output workspace
This algorithm searches for peaks in a dataset.
Definition: FindPeaks.h:50
int getIndex(const HistogramData::HistogramX &vecX, double x)
needed by FindPeaksBackground
Definition: FindPeaks.cpp:702
Support for a property that holds an array of values.
Definition: ArrayProperty.h:28
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void debug(const std::string &msg)
Logs at debug level.
Definition: Logger.cpp:114
void information(const std::string &msg)
Logs at information level.
Definition: Logger.cpp:105
void report()
Increments the loop counter by 1, then sends the progress notification on behalf of its algorithm.
Definition: ProgressBase.h:51
Statistics getStatistics(const std::vector< TYPE > &data, const unsigned int flags=StatOptions::AllStats)
Return a statistics object for the given data set.
Definition: Statistics.cpp:167
Helper class which provides the Collimation Length for SANS instruments.
constexpr int EMPTY_INT() noexcept
Returns what we consider an "empty" integer within a property.
Definition: EmptyValues.h:25
std::vector< double > MantidVec
typedef for the data storage used in Mantid matrix workspaces
Definition: cow_ptr.h:172
STL namespace.
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54
Simple struct to store statistics.
Definition: Statistics.h:25
double mean
Mean value.
Definition: Statistics.h:31
double standard_deviation
standard_deviation of the values
Definition: Statistics.h:35