Merge pull request #20857 from mantidproject/findpeakbackground_refactor_20796

Findpeakbackground refactor 20796

Framework/Algorithms/inc/MantidAlgorithms/FindPeakBackground.h

@@ -2,6 +2,8 @@
 #define MANTID_ALGORITHMS_FindPeakBackground_H_
 
 #include "MantidAPI/Algorithm.h"
+#include "MantidAPI/ITableWorkspace.h"
+#include "MantidHistogramData/Histogram.h"
 #include "MantidKernel/cow_ptr.h"
 
 namespace Mantid {
@@ -51,6 +53,25 @@ public:
   /// Algorithm's category for identification overriding a virtual method
   const std::string category() const override { return "Utility\\Calculation"; }
 
+  /// 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 fit window's data point index
+  void findWindowIndex(const HistogramData::Histogram &histogram, size_t &l0,
+                       size_t &n);
+
+  /// main method to calculate background
+  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);
+
 private:
   std::string m_backgroundType; //< The type of background to fit
 
@@ -65,6 +86,29 @@ private:
                           const size_t p_min, const size_t p_max,
                           const bool hasPeak, double &out_bg0, double &out_bg1,
                           double &out_bg2);
+
+  /// process inputs
+  void processInputProperties();
+
+  /// create output workspace
+  void createOutputWorkspaces();
+
+  // Histogram cannot be defined due to lack of default constructor. shared_ptr
+  // will do the copy
+
+  /// fit window
+  std::vector<double> m_vecFitWindows;
+  /// background order: 0 for flat, 1 for linear, 2 for quadratic
+  size_t m_backgroundOrder;
+  /// constant sigma
+  double m_sigmaConstant;
+  /// output workspace (table of result)
+  API::ITableWorkspace_sptr m_outPeakTableWS;
+  /// Input workspace
+  API::MatrixWorkspace_const_sptr m_inputWS;
+  /// workspace index
+  size_t m_inputWSIndex;
+
   struct cont_peak {
     size_t start;
     size_t stop;

Framework/Algorithms/src/FindPeakBackground.cpp

@@ -64,40 +64,15 @@ void FindPeakBackground::init() {
       "quadratic terms.");
 }
 
-//----------------------------------------------------------------------------------------------
-/** Execute body
-  */
-void FindPeakBackground::exec() {
-  // 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;
-
+void FindPeakBackground::findWindowIndex(
+    const HistogramData::Histogram &histogram, size_t &l0, size_t &n) {
+  auto &inpX = histogram.x();
+  auto &inpY = histogram.y();
+  size_t sizey = inpY.size(); // inpWS->y(inpwsindex).size();
+
+  // determine the fit window with their index in X (or Y)
+  n = sizey;
+  l0 = 0;
   if (m_vecFitWindows.size() > 1) {
     Mantid::Algorithms::FindPeaks fp;
     l0 = fp.getIndex(inpX, m_vecFitWindows[0]);
@@ -105,27 +80,70 @@ void FindPeakBackground::exec() {
     if (n < sizey)
       n++;
   }
+}
 
-  // 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");
+//----------------------------------------------------------------------------------------------
+/** Execute body
+  */
+void FindPeakBackground::exec() {
+  // Get input and validate
+  processInputProperties();
+  auto histogram = m_inputWS->histogram(m_inputWSIndex);
 
-  m_outPeakTableWS->appendRow();
+  size_t l0, n;
+  findWindowIndex(histogram, l0, n);
+
+  // m_vecFitWindows won't be used again form this point till end.
+
+  // Set up output table workspace
+  createOutputWorkspaces();
 
   // 3. Get Y values
   Progress prog(this, 0.0, 1.0, 1);
 
-  // Find background
+  std::vector<size_t> peak_min_max_indexes;
+  std::vector<double> bkgd3;
+  int goodfit = findBackground(histogram, l0, n, peak_min_max_indexes, bkgd3);
+
+  if (goodfit > 0) {
+    size_t min_peak = peak_min_max_indexes[0];
+    size_t max_peak = peak_min_max_indexes[1];
+    double a0 = bkgd3[0];
+    double a1 = bkgd3[1];
+    double a2 = bkgd3[2];
+    API::TableRow t = m_outPeakTableWS->getRow(0);
+    t << static_cast<int>(m_inputWSIndex) << static_cast<int>(min_peak)
+      << static_cast<int>(max_peak) << a0 << a1 << a2 << goodfit;
+  }
+
+  prog.report();
+
+  // 4. Set the output
+  setProperty("OutputWorkspace", m_outPeakTableWS);
+}
 
-  auto &inpY = inpWS->y(inpwsindex);
+//----------------------------------------------------------------------------------------------
+/**
+ * @brief FindPeakBackground::findBackground
+ * @param histogram
+ * @param l0
+ * @param n
+ * @param peak_min_max_indexes
+ * @param bkgd3
+ * @return
+ */
+int FindPeakBackground::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) {
+  auto &inpX = histogram.x();
+  auto &inpY = histogram.y();
+  size_t sizex = inpX.size(); // inpWS->x(inpwsindex).size();
+  size_t sizey = inpY.size(); // inpWS->y(inpwsindex).size();
+
+  int goodfit(0);
 
+  // Find background
   double Ymean, Yvariance, Ysigma;
   MantidVec maskedY;
   auto in = std::min_element(inpY.cbegin(), inpY.cend());
@@ -153,7 +171,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
@@ -196,7 +214,6 @@ void FindPeakBackground::exec() {
     }
     size_t min_peak, max_peak;
     double a0 = 0., a1 = 0., a2 = 0.;
-    int goodfit;
     if (!peaks.empty()) {
       g_log.debug() << "Peaks' size = " << peaks.size()
                     << " -> esitmate background. \n";
@@ -217,19 +234,22 @@ 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)
-      << static_cast<int>(max_peak) << a0 << a1 << a2 << goodfit;
+    peak_min_max_indexes.resize(2);
+    peak_min_max_indexes[0] = min_peak;
+    peak_min_max_indexes[1] = max_peak;
+
+    bkgd3.resize(3);
+    bkgd3[0] = a0;
+    bkgd3[1] = a1;
+    bkgd3[2] = a2;
   }
 
-  prog.report();
-
-  // 4. Set the output
-  setProperty("OutputWorkspace", m_outPeakTableWS);
+  return goodfit;
 }
 
 //----------------------------------------------------------------------------------------------
@@ -400,5 +420,92 @@ 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
+  m_inputWS = getProperty("InputWorkspace");
+
+  int inpwsindex = getProperty("WorkspaceIndex");
+  if (isEmpty(inpwsindex)) {
+    // Default
+    if (m_inputWS->getNumberHistograms() == 1) {
+      inpwsindex = 0;
+    } else {
+      throw runtime_error("WorkspaceIndex must be given. ");
+    }
+  } else if (inpwsindex < 0 ||
+             inpwsindex >= static_cast<int>(m_inputWS->getNumberHistograms())) {
+    stringstream errss;
+    errss << "Input workspace " << m_inputWS->getName() << " has "
+          << m_inputWS->getNumberHistograms()
+          << " spectra.  Input workspace index " << inpwsindex
+          << " is out of boundary. ";
+    throw runtime_error(errss.str());
+  }
+  m_inputWSIndex = static_cast<size_t>(inpwsindex);
+
+  std::vector<double> fitwindow = getProperty("FitWindow");
+  setFitWindow(fitwindow);
+
+  // background
+  m_backgroundType = getPropertyValue("BackgroundType");
+  size_t bkgdorder = 0;
+  if (m_backgroundType == "Linear")
+    bkgdorder = 1;
+  else if (m_backgroundType == "Quadratic")
+    bkgdorder = 2;
+  setBackgroundOrder(bkgdorder);
+
+  // sigma constant
+  double k = getProperty("SigmaConstant");
+  setSigma(k);
+}
+
+/// set sigma constant
+void FindPeakBackground::setSigma(const double &sigma) {
+  m_sigmaConstant = sigma;
+}
+
+/// set background order
+void FindPeakBackground::setBackgroundOrder(size_t order) {
+  m_backgroundOrder = order;
+}
+
+//----------------------------------------------------------------------------------------------
+/** set fit window
+ * @brief FindPeakBackground::setFitWindow
+ * @param fitwindow
+ */
+void FindPeakBackground::setFitWindow(const std::vector<double> &fitwindow) {
+  // validate input
+  if ((fitwindow.size() == 2) && fitwindow[0] >= fitwindow[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;
+}
+
+//----------------------------------------------------------------------------------------------
+/**
+ * @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();
+}
+
 } // namespace Algorithms
 } // namespace Mantid