add signal to sigma post fitting check

Framework/Algorithms/inc/MantidAlgorithms/FitPeaks.h

@@ -216,6 +216,10 @@ class MANTID_ALGORITHMS_DLL FitPeaks final : public API::Algorithm {
   /// calculate peak+background for fitted
   void calculateFittedPeaks(const std::vector<std::shared_ptr<FitPeaksAlgorithm::PeakFitResult>> &fit_results);
 
+  double calculateSignalToSigmaRatio(const size_t &iws, const std::pair<double, double> &peakWindow,
+                                     const API::IPeakFunction_sptr &peakFunction,
+                                     const API::IBackgroundFunction_sptr &backgroundFunction);
+
   /// Get the parameter name for peak height (I or height or etc)
   std::string getPeakHeightParameterName(const API::IPeakFunction_const_sptr &peak_function);
 
@@ -329,6 +333,8 @@ class MANTID_ALGORITHMS_DLL FitPeaks final : public API::Algorithm {
   double m_minSignalToNoiseRatio;
   double m_minPeakTotalCount;
 
+  double m_minSignalToSigmaRatio;
+
   /// flag for high background
   bool m_highBackground;
 

Framework/Algorithms/src/FitPeaks.cpp

@@ -33,6 +33,7 @@
 #include "MantidKernel/IValidator.h"
 #include "MantidKernel/ListValidator.h"
 #include "MantidKernel/StartsWithValidator.h"
+#include "MantidKernel/VectorHelper.h"
 
 #include "boost/algorithm/string.hpp"
 #include "boost/algorithm/string/trim.hpp"
@@ -82,6 +83,7 @@ const std::string OUTPUT_WKSP_PARAM_ERRS("OutputParameterFitErrorsWorkspace");
 const std::string RAW_PARAMS("RawPeakParameters");
 const std::string PEAK_MIN_SIGNAL_TO_NOISE_RATIO("MinimumSignalToNoiseRatio");
 const std::string PEAK_MIN_TOTAL_COUNT("MinimumPeakTotalCount");
+const std::string PEAK_MIN_SIGNAL_TO_SIGMA_RATIO("MinimumSignalToSigmaRatio");
 } // namespace PropertyNames
 } // namespace
 
@@ -431,6 +433,10 @@ void FitPeaks::init() {
                   "Used for validating peaks before fitting. If the total peak window Y-value count "
                   "is under this value, the peak will be excluded from fitting and calibration.");
 
+  declareProperty(PropertyNames::PEAK_MIN_SIGNAL_TO_SIGMA_RATIO, 0.,
+                  "Used for validating peaks before fitting. If the signal-to-sigma ratio is under this value, "
+                  "the peak will be excluded from fitting and calibration.");
+
   const std::string addoutgrp("Analysis");
   setPropertyGroup(PropertyNames::OUTPUT_WKSP_PARAMS, addoutgrp);
   setPropertyGroup(PropertyNames::OUTPUT_WKSP_MODEL, addoutgrp);
@@ -923,6 +929,11 @@ void FitPeaks::processInputPeakTolerance() {
   m_minSignalToNoiseRatio = getProperty(PropertyNames::PEAK_MIN_SIGNAL_TO_NOISE_RATIO);
   if (isEmpty(m_minSignalToNoiseRatio) || m_minSignalToNoiseRatio < 0.)
     m_minSignalToNoiseRatio = 0.;
+
+  // set the signal-to-sigma threshold to zero (default value) if not specified or invalid
+  m_minSignalToSigmaRatio = getProperty(PropertyNames::PEAK_MIN_SIGNAL_TO_SIGMA_RATIO);
+  if (isEmpty(m_minSignalToSigmaRatio) || m_minSignalToSigmaRatio < 0.)
+    m_minSignalToSigmaRatio = 0.;
 }
 
 //----------------------------------------------------------------------------------------------
@@ -1282,6 +1293,14 @@ void FitPeaks::fitSpectrumPeaks(size_t wi, const std::vector<double> &expected_p
                                bkgdfunction, peak_pre_check_result);
       if (peak_pre_check_result->isIndividualPeakRejected())
         fit_result->setBadRecord(peak_index, -1.);
+
+      if (m_minSignalToSigmaRatio > 0) {
+        if (calculateSignalToSigmaRatio(wi, peak_window_i, peakfunction, bkgdfunction); < m_minSignalToSigmaRatio) {
+          fit_result->setBadRecord(peak_index, -1.);
+          cost = DBL_MAX;
+        }
+      }
+
       *pre_check_result += *peak_pre_check_result; // keep track of the rejection count within the spectrum
     }
     pre_check_result->setNumberOfOutOfRangePeaks(number_of_out_of_range_peaks);
@@ -1531,6 +1550,38 @@ void FitPeaks::calculateFittedPeaks(const std::vector<std::shared_ptr<FitPeaksAl
 
   return;
 }
+
+double FitPeaks::calculateSignalToSigmaRatio(const size_t &iws, const std::pair<double, double> &peakWindow,
+                                             const API::IPeakFunction_sptr &peakFunction,
+                                             const API::IBackgroundFunction_sptr &backgroundFunction) {
+  const auto &vecX = m_fittedPeakWS->points(iws);
+  auto startX = std::lower_bound(vecX.begin(), vecX.end(), peakWindow.first);
+  auto stopX = std::lower_bound(vecX.begin(), vecX.end(), peakWindow.second);
+
+  FunctionDomain1DVector domain(startX, stopX);
+  FunctionValues values(domain);
+  CompositeFunction_sptr compFunc = std::make_shared<API::CompositeFunction>();
+
+  compFunc->addFunction(backgroundFunction);
+  compFunc->function(domain, values);
+  auto bkgdValues = values.toVector();
+
+  compFunc->addFunction(peakFunction);
+  compFunc->function(domain, values);
+  auto fittedValues = values.toVector();
+
+  const auto &errors = m_fittedPeakWS->readE(iws);
+  auto startE = std::lower_bound(errors.begin(), errors.end(), peakWindow.first);
+  auto stopE = std::lower_bound(errors.begin(), errors.end(), peakWindow.second);
+  std::vector<double> peakErrors(startE, stopE);
+
+  double fittedSum = std::accumulate(fittedValues.cbegin(), fittedValues.cend(), 0.0);
+  double bkgdSum = std::accumulate(bkgdValues.cbegin(), bkgdValues.cend(), 0.0);
+  double sigma = sqrt(std::accumulate(peakErrors.cbegin(), peakErrors.cend(), 0.0, VectorHelper::SumSquares<double>()));
+
+  return (fittedSum - bkgdSum) / ((sigma == 0) ? 1 : sigma);
+}
+
 namespace {
 bool estimateBackgroundParameters(const Histogram &histogram, const std::pair<size_t, size_t> &peak_window,
                                   const API::IBackgroundFunction_sptr &bkgd_function) {

Framework/Algorithms/src/PDCalibration.cpp

@@ -73,7 +73,7 @@ DECLARE_ALGORITHM(PDCalibration)
 
 namespace { // anonymous
 const auto isNonZero = [](const double value) { return value != 0.; };
-}
+} // namespace
 
 /// private inner class
 class PDCalibration::FittedPeaks {
@@ -315,6 +315,10 @@ void PDCalibration::init() {
                   "Used for validating peaks before fitting. If the total peak window Y-value count "
                   "is under this value, the peak will be excluded from fitting and calibration.");
 
+  declareProperty("MinimumSignalToSigmaRatio", 0.,
+                  "Used for validating peaks before fitting. If the signal-to-sigma ratio is under this value, "
+                  "the peak will be excluded from fitting and calibration.");
+
   // make group for Input properties
   std::string inputGroup("Input Options");
   setPropertyGroup("InputWorkspace", inputGroup);
@@ -336,6 +340,7 @@ void PDCalibration::init() {
   setPropertyGroup("MinimumPeakHeight", fitPeaksGroup);
   setPropertyGroup("MinimumSignalToNoiseRatio", fitPeaksGroup);
   setPropertyGroup("MinimumPeakTotalCount", fitPeaksGroup);
+  setPropertyGroup("MinimumSignalToSigmaRatio", fitPeaksGroup);
   setPropertyGroup("HighBackground", fitPeaksGroup);
   setPropertyGroup("MaxChiSq", fitPeaksGroup);
   setPropertyGroup("ConstrainPeakPositions", fitPeaksGroup);
@@ -489,6 +494,7 @@ void PDCalibration::exec() {
   const double minPeakHeight = getProperty("MinimumPeakHeight");
   const double minPeakTotalCount = getProperty("MinimumPeakTotalCount");
   const double minSignalToNoiseRatio = getProperty("MinimumSignalToNoiseRatio");
+  const double minSignalToSigmaRatio = getProperty("MinimumSignalToSigmaRatio");
   const double maxChiSquared = getProperty("MaxChiSq");
 
   const std::string calParams = getPropertyValue("CalibrationParameters");
@@ -578,6 +584,7 @@ void PDCalibration::exec() {
   algFitPeaks->setProperty("MinimumPeakHeight", minPeakHeight);
   algFitPeaks->setProperty("MinimumPeakTotalCount", minPeakTotalCount);
   algFitPeaks->setProperty("MinimumSignalToNoiseRatio", minSignalToNoiseRatio);
+  algFitPeaks->setProperty("MinimumSignalToSigmaRatio", minSignalToSigmaRatio);
   // some fitting strategy
   algFitPeaks->setProperty("FitFromRight", true);
   const bool highBackground = getProperty("HighBackground");