diff --git a/Framework/Algorithms/inc/MantidAlgorithms/FindPeakBackground.h b/Framework/Algorithms/inc/MantidAlgorithms/FindPeakBackground.h
index e318a2f00947adfe79d8821b5efda0c71559335e..62ce8bda1e0f32d9677e690ff1b173b7a76bcb88 100644
--- a/Framework/Algorithms/inc/MantidAlgorithms/FindPeakBackground.h
+++ b/Framework/Algorithms/inc/MantidAlgorithms/FindPeakBackground.h
@@ -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
diff --git a/Framework/Algorithms/src/FindPeakBackground.cpp b/Framework/Algorithms/src/FindPeakBackground.cpp
index 324ab4492871de6ee7ca31750141147c9768766b..ec60398fc8b0cfe996924c09fe254520b4c830bf 100644
--- a/Framework/Algorithms/src/FindPeakBackground.cpp
+++ b/Framework/Algorithms/src/FindPeakBackground.cpp
@@ -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