Skip to content
Snippets Groups Projects
Commit 9b8e990e authored by Zhou, Wenduo's avatar Zhou, Wenduo
Browse files

Refs #20796. Continued to refactor (unit test passed).

parent 166b4d7d
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Framework/Algorithms/inc/MantidAlgorithms/FindPeakBackground.h
+ 37
− 35
View file @ 9b8e990e
......@@ -3,6 +3,7 @@
#include "MantidAPI/Algorithm.h"
#include "MantidAPI/ITableWorkspace.h"
#include "MantidHistogramData/Histogram.h"
#include "MantidKernel/cow_ptr.h"
namespace Mantid {
......@@ -37,7 +38,6 @@ namespace Algorithms {
File change history is stored at: <https://github.com/mantidproject/mantid>
Code Documentation is available at: <http://doxygen.mantidproject.org>
*/
class DLLExport FindPeakBackground : public API::Algorithm {
public:
/// Algorithm's name for identification overriding a virtual method
......@@ -50,30 +50,9 @@ public:
/// Algorithm's version for identification overriding a virtual method
int version() const override { return 1; }
/// process inputs
void processInputProperties();
/// Algorithm's category for identification overriding a virtual method
const std::string category() const override { return "Utility\\Calculation"; }
/// set histogram data to find background
// void setHistogram(HistogramData &histogram);
/// set sigma constant
void setSigma(const double &sigma);
/// set background order
void setBackgroundOrder(size_t order);
/// set fit window
void setFitWindow(const std::vector<double> &window);
/// find background (main algorithm)
void findPeakBackground();
/// get result
void getBackgroundResult();
private:
std::string m_backgroundType; //< The type of background to fit
......@@ -89,26 +68,30 @@ private:
const bool hasPeak, double &out_bg0, double &out_bg1,
double &out_bg2);
/// process inputs
void processInputProperties();
/// create output workspace
void createOutputWorkspaces();
struct cont_peak {
size_t start;
size_t stop;
double maxY;
};
struct by_len {
bool operator()(cont_peak const &a, cont_peak const &b) {
return a.maxY > b.maxY;
}
};
void findStartStopIndex(size_t &istart, size_t &istop);
// define parameters
/// set histogram data to find background
void setHistogram(const HistogramData::Histogram &histogram);
/// set sigma constant
void setSigma(const double &sigma);
/// set background order
void setBackgroundOrder(size_t order);
/// set fit window
void setFitWindow(const std::vector<double> &window);
/// histogram data to find peak background
// HistogramData::Histogram m_histogram;
boost::shared_ptr<const HistogramData::Histogram> m_histogram;
// HistogramData::Histogram& m_histogram;
/// fit window
std::vector<double> m_vecFitWindows;
/// background order: 0 for flat, 1 for linear, 2 for quadratic
......@@ -117,6 +100,25 @@ private:
double m_sigmaConstant;
/// output workspace (table of result)
API::ITableWorkspace_sptr m_outPeakTableWS;
size_t m_inputWSIndex;
// /// find background (main algorithm)
// void findPeakBackground();
// /// get result
// void getBackgroundResult();
struct cont_peak {
size_t start;
size_t stop;
double maxY;
};
struct by_len {
bool operator()(cont_peak const &a, cont_peak const &b) {
return a.maxY > b.maxY;
}
};
};
} // namespace Algorithms
......
Framework/Algorithms/src/FindPeakBackground.cpp
+ 148
− 183
View file @ 9b8e990e
......@@ -65,117 +65,20 @@ void FindPeakBackground::init() {
}
//----------------------------------------------------------------------------------------------
void FindPeakBackground::processInputProperties() {
// process input workspace and workspace index
MatrixWorkspace_const_sptr inpWS = getProperty("InputWorkspace");
int inpwsindex = getProperty("WorkspaceIndex");
if (isEmpty(inpwsindex)) {
// Default
if (inpWS->getNumberHistograms() == 1) {
inpwsindex = 0;
} else {
throw runtime_error("WorkspaceIndex must be given. ");
}
} else if (inpwsindex < 0 ||
inpwsindex >= static_cast<int>(inpWS->getNumberHistograms())) {
stringstream errss;
errss << "Input workspace " << inpWS->getName() << " has "
<< inpWS->getNumberHistograms() << " spectra. Input workspace index "
<< inpwsindex << " is out of boundary. ";
throw runtime_error(errss.str());
}
setHistogram(inputWS->histogram(wi));
std::vector<double> fitwindow = getProperty("FitWindow");
setFitWindow(fitwindow);
// background
m_backgroundType = getPropertyValue("BackgroundType");
size_t bkgdorder = 0;
if (m_backgroundType.compare("Linear") == 0)
bkgdorder = 1;
else if (m_backgroundType.compare("Quadratic") == 0)
bkgdorder = 2;
setBackgroundOrder(bkgdorder);
// sigma constant
double k = getProperty("SigmaConstant");
setSigma(k);
}
/// set histogram data to find background
void FindPeakBackground::setHistogram(HistogramData &histogram) {
m_histogram = histogram;
}
/// set sigma constant
void FindPeakBackground::setSigma(const double &sigma) {
m_sigmaConstant = sigma;
}
/// set background order
void FindPeakBackground::setBackgroundOrder(size_t order) {
m_backgroundOrder = order;
}
//----------------------------------------------------------------------------------------------
/** find background (main algorithm)
* @brief FindPeakBackground::findPeakBackground
*/
void FindPeakBackground::findPeakBackground() {}
/// get result
void getBackgroundResult();
//----------------------------------------------------------------------------------------------
/** set fit window
* @brief FindPeakBackground::setFitWindow
* @param fitwindow
*/
void FindPeakBackground::setFitWindow(const std::vector<double> &fitwindow) {
// validate input
if (m_vecFitWindows.size() != 2 || m_vecFitWindows[0] >= m_vecFitWindows[1])
throw std::invalid_argument("Fit window has either wrong item number or "
"window value is not in ascending order.");
// m_vecFitWindows.resize(2);
// copy the input to class variable
m_vecFitWindows = fitwindow;
return;
}
//----------------------------------------------------------------------------------------------
/**
* @brief FindPeakBackground::createOutputWorkspaces
*/
void FindPeakBackground::createOutputWorkspaces() {
// Set up output table workspace
m_outPeakTableWS = WorkspaceFactory::Instance().createTable("TableWorkspace");
m_outPeakTableWS->addColumn("int", "wksp_index");
m_outPeakTableWS->addColumn("int", "peak_min_index");
m_outPeakTableWS->addColumn("int", "peak_max_index");
m_outPeakTableWS->addColumn("double", "bkg0");
m_outPeakTableWS->addColumn("double", "bkg1");
m_outPeakTableWS->addColumn("double", "bkg2");
m_outPeakTableWS->addColumn("int", "GoodFit");
m_outPeakTableWS->appendRow();
}
/** Execute body
*/
void FindPeakBackground::exec() {
// Get input and validate
processInputProperties();
void FindPeakBackground::findStartStopIndex(size_t &istart, size_t &istop) {
// Generate output
// auto &inpX = inpWS->x(inpwsindex);
auto inpX = m_histogram.x();
auto inpY = m_histogram.y();
size_t sizex = inpX.size();
size_t sizey = inpY.size();
auto &inpX = m_histogram->x();
auto &inpY = m_histogram->y();
size_t sizex = inpX.size(); // inpWS->x(inpwsindex).size();
size_t sizey = inpY.size(); // inpWS->y(inpwsindex).size();
// initial value of start and stop x index
// determine the fit window with their index in X (or Y)
size_t n = sizey;
size_t l0 = 0;
if (m_vecFitWindows.size() > 1) {
Mantid::Algorithms::FindPeaks fp;
l0 = fp.getIndex(inpX, m_vecFitWindows[0]);
......@@ -183,91 +86,19 @@ void FindPeakBackground::findStartStopIndex(size_t &istart, size_t &istop) {
if (n < sizey)
n++;
}
// m_vecFitWindows won't be used again form this point till end.
istart = l0;
istop = n;
}
//----------------------------------------------------------------------------------------------
/** Execute body
*/
void FindPeakBackground::exec() {
// get input and validate
processInputProperties();
// Get input and validate
// MatrixWorkspace_const_sptr inpWS = getProperty("InputWorkspace");
// int inpwsindex = getProperty("WorkspaceIndex");
// std::vector<double> m_vecFitWindows = getProperty("FitWindow");
// m_backgroundType = getPropertyValue("BackgroundType");
// double k = getProperty("SigmaConstant");
// if (isEmpty(inpwsindex)) {
// // Default
// if (inpWS->getNumberHistograms() == 1) {
// inpwsindex = 0;
// } else {
// throw runtime_error("WorkspaceIndex must be given. ");
// }
// } else if (inpwsindex < 0 ||
// inpwsindex >= static_cast<int>(inpWS->getNumberHistograms())) {
// stringstream errss;
// errss << "Input workspace " << inpWS->getName() << " has "
// << inpWS->getNumberHistograms() << " spectra. Input workspace
// index "
// << inpwsindex << " is out of boundary. ";
// throw runtime_error(errss.str());
// }
// // Generate output
// auto &inpX = inpWS->x(inpwsindex);
// size_t sizex = inpWS->x(inpwsindex).size();
// size_t sizey = inpWS->y(inpwsindex).size();
// size_t n = sizey;
// size_t l0 = 0;
// if (m_vecFitWindows.size() > 1) {
// Mantid::Algorithms::FindPeaks fp;
// l0 = fp.getIndex(inpX, m_vecFitWindows[0]);
// n = fp.getIndex(inpX, m_vecFitWindows[1]);
// if (n < sizey)
// n++;
// }
// Set up output table workspace
createOutputWorkspaces();
size_t istart, istop;
findStartStopIndex(istart, istop);
// // Set up output table workspace
// API::ITableWorkspace_sptr m_outPeakTableWS =
// WorkspaceFactory::Instance().createTable("TableWorkspace");
// m_outPeakTableWS->addColumn("int", "wksp_index");
// m_outPeakTableWS->addColumn("int", "peak_min_index");
// m_outPeakTableWS->addColumn("int", "peak_max_index");
// m_outPeakTableWS->addColumn("double", "bkg0");
// m_outPeakTableWS->addColumn("double", "bkg1");
// m_outPeakTableWS->addColumn("double", "bkg2");
// m_outPeakTableWS->addColumn("int", "GoodFit");
// m_outPeakTableWS->appendRow();
// 3. Get Y values
Progress prog(this, 0.0, 1.0, 1);
// Find background
auto &inpY = inpWS->y(inpwsindex);
double Ymean, Yvariance, Ysigma;
MantidVec maskedY;
auto in = std::min_element(inpY.cbegin(), inpY.cend());
double bkg0 = inpY[in - inpY.begin()];
// WZ:
size_t l0 = istart;
size_t n = istop;
// -----
for (size_t l = l0; l < n; ++l) {
maskedY.push_back(inpY[l] - bkg0);
}
......@@ -291,7 +122,7 @@ void FindPeakBackground::exec() {
const size_t pos = it - maskedY.begin();
maskedY[pos] = 0;
mask[pos] = 1.0;
} while (std::abs(Ymean - Yvariance) > k * Ysigma);
} while (std::abs(Ymean - Yvariance) > m_sigmaConstant * Ysigma);
if (n - l0 > 5) {
// remove single outliers
......@@ -355,12 +186,13 @@ void FindPeakBackground::exec() {
goodfit = 2;
}
estimateBackground(inpX, inpY, l0, n, min_peak, max_peak, (!peaks.empty()),
a0, a1, a2);
// Add a new row
API::TableRow t = m_outPeakTableWS->getRow(0);
t << static_cast<int>(inpwsindex) << static_cast<int>(min_peak)
t << static_cast<int>(m_inputWSIndex) << static_cast<int>(min_peak)
<< static_cast<int>(max_peak) << a0 << a1 << a2 << goodfit;
}
......@@ -538,5 +370,138 @@ double FindPeakBackground::moment4(MantidVec &X, size_t n, double mean) {
sum /= static_cast<double>(n);
return sum;
}
//----------------------------------------------------------------------------------------------
void FindPeakBackground::processInputProperties() {
// process input workspace and workspace index
MatrixWorkspace_const_sptr inpWS = getProperty("InputWorkspace");
int inpwsindex = getProperty("WorkspaceIndex");
if (isEmpty(inpwsindex)) {
// Default
if (inpWS->getNumberHistograms() == 1) {
inpwsindex = 0;
} else {
throw runtime_error("WorkspaceIndex must be given. ");
}
} else if (inpwsindex < 0 ||
inpwsindex >= static_cast<int>(inpWS->getNumberHistograms())) {
stringstream errss;
errss << "Input workspace " << inpWS->getName() << " has "
<< inpWS->getNumberHistograms() << " spectra. Input workspace index "
<< inpwsindex << " is out of boundary. ";
throw runtime_error(errss.str());
}
m_inputWSIndex = static_cast<size_t>(inpwsindex);
setHistogram(inpWS->histogram(inpwsindex));
std::vector<double> fitwindow = getProperty("FitWindow");
setFitWindow(fitwindow);
// background
m_backgroundType = getPropertyValue("BackgroundType");
size_t bkgdorder = 0;
if (m_backgroundType.compare("Linear") == 0)
bkgdorder = 1;
else if (m_backgroundType.compare("Quadratic") == 0)
bkgdorder = 2;
setBackgroundOrder(bkgdorder);
// sigma constant
double k = getProperty("SigmaConstant");
setSigma(k);
}
/// set histogram data to find background
void FindPeakBackground::setHistogram(
const HistogramData::Histogram &histogram) {
m_histogram = boost::make_shared<HistogramData::Histogram>(histogram);
}
/// set sigma constant
void FindPeakBackground::setSigma(const double &sigma) {
m_sigmaConstant = sigma;
}
/// set background order
void FindPeakBackground::setBackgroundOrder(size_t order) {
m_backgroundOrder = order;
}
////----------------------------------------------------------------------------------------------
///** find background (main algorithm)
// * @brief FindPeakBackground::findPeakBackground
// */
// void FindPeakBackground::findPeakBackground() {return;}
///// get result
// void getBackgroundResult();
//----------------------------------------------------------------------------------------------
/** set fit window
* @brief FindPeakBackground::setFitWindow
* @param fitwindow
*/
void FindPeakBackground::setFitWindow(const std::vector<double> &fitwindow) {
// validate input
if (m_vecFitWindows.size() == 0) {
m_vecFitWindows.resize(2);
m_vecFitWindows[0] = m_histogram->y().front();
m_vecFitWindows[1] = m_histogram->y().back();
} else if (m_vecFitWindows.size() != 2 ||
m_vecFitWindows[0] >= m_vecFitWindows[1]) {
throw std::invalid_argument("Fit window has either wrong item number or "
"window value is not in ascending order.");
}
// m_vecFitWindows.resize(2);
// copy the input to class variable
m_vecFitWindows = fitwindow;
return;
}
//----------------------------------------------------------------------------------------------
/**
* @brief FindPeakBackground::createOutputWorkspaces
*/
void FindPeakBackground::createOutputWorkspaces() {
// Set up output table workspace
m_outPeakTableWS = WorkspaceFactory::Instance().createTable("TableWorkspace");
m_outPeakTableWS->addColumn("int", "wksp_index");
m_outPeakTableWS->addColumn("int", "peak_min_index");
m_outPeakTableWS->addColumn("int", "peak_max_index");
m_outPeakTableWS->addColumn("double", "bkg0");
m_outPeakTableWS->addColumn("double", "bkg1");
m_outPeakTableWS->addColumn("double", "bkg2");
m_outPeakTableWS->addColumn("int", "GoodFit");
m_outPeakTableWS->appendRow();
}
void FindPeakBackground::findStartStopIndex(size_t &istart, size_t &istop) {
// Generate output
auto inpX = m_histogram->x();
auto inpY = m_histogram->y();
// size_t sizex = inpX.size();
size_t sizey = inpY.size();
// initial value of start and stop x index
size_t n = sizey;
size_t l0 = 0;
if (m_vecFitWindows.size() > 1) {
Mantid::Algorithms::FindPeaks fp;
l0 = fp.getIndex(inpX, m_vecFitWindows[0]);
n = fp.getIndex(inpX, m_vecFitWindows[1]);
if (n < sizey)
n++;
}
istart = l0;
istop = n;
}
} // namespace Algorithms
} // namespace Mantid
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment