add signal to sigma post fitting check

Author: MohamedAlmaki

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);
@@ -445,7 +451,7 @@ std::map<std::string, std::string> FitPeaks::validateInputs() {
   map<std::string, std::string> issues;
 
   // check that min/max spectra indices make sense - only matters if both are specified
-  if (!(isDefault(PropertyNames::START_WKSP_INDEX) && isDefault(PropertyNames::START_WKSP_INDEX))) {
+  if (!(isDefault(PropertyNames::START_WKSP_INDEX) && isDefault(PropertyNames::STOP_WKSP_INDEX))) {
     const int startIndex = getProperty(PropertyNames::START_WKSP_INDEX);
     const int stopIndex = getProperty(PropertyNames::STOP_WKSP_INDEX);
     if (startIndex > stopIndex) {
@@ -922,6 +928,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.;
 }
 
 //----------------------------------------------------------------------------------------------
@@ -1281,6 +1292,13 @@ 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.);
+
+      auto signalToSigma = calculateSignalToSigmaRatio(wi, peak_window_i, peakfunction, bkgdfunction);
+      if (m_minSignalToSigmaRatio > 0 && signalToSigma < 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);
@@ -1530,6 +1548,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");

buildconfig/CMake/CppCheck_Suppressions.txt.in

@@ -256,7 +256,6 @@ knownConditionTrueFalse:${CMAKE_SOURCE_DIR}/Framework/Algorithms/src/FindPeaks.c
 unreadVariable:${CMAKE_SOURCE_DIR}/Framework/Algorithms/src/InterpolatingRebin.cpp:69
 constVariableReference:${CMAKE_SOURCE_DIR}/Framework/Algorithms/src/GetEi2.cpp:213
 missingOverride:${CMAKE_SOURCE_DIR}/Framework/Algorithms/inc/MantidAlgorithms/GetAllEi.h:32
-duplicateExpression:${CMAKE_SOURCE_DIR}/Framework/Algorithms/src/FitPeaks.cpp:448
 passedByValue:${CMAKE_SOURCE_DIR}/Framework/Algorithms/src/GetTimeSeriesLogInformation.cpp:277
 passedByValue:${CMAKE_SOURCE_DIR}/Framework/Algorithms/src/GetTimeSeriesLogInformation.cpp:501
 constVariableReference:${CMAKE_SOURCE_DIR}/Framework/Algorithms/src/GetTimeSeriesLogInformation.cpp:248