Merge pull request #37826 from mantidproject/35224_add_post_fitting

Add signal to sigma post fitting check
mantidproject · Sep 12, 2024 · c36a1c9 · c36a1c9
2 parents 0a30598 + c115134
commit c36a1c9
Show file tree

Hide file tree

Showing 5 changed files with 150 additions and 1 deletion.
diff --git a/Framework/Algorithms/inc/MantidAlgorithms/FitPeaks.h b/Framework/Algorithms/inc/MantidAlgorithms/FitPeaks.h
@@ -216,6 +216,9 @@ 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);
+
   /// Get the parameter name for peak height (I or height or etc)
   std::string getPeakHeightParameterName(const API::IPeakFunction_const_sptr &peak_function);
 
@@ -329,6 +332,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;
 

diff --git a/Framework/Algorithms/src/FitPeaks.cpp b/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 after 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) < 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,30 @@ 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 auto &vecX = m_inputMatrixWS->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);
+
+  peakFunction->function(domain, values);
+  auto peakValues = values.toVector();
+
+  const auto &errors = m_inputMatrixWS->readE(iws);
+  auto startE = errors.begin() + (startX - vecX.begin());
+  auto stopE = errors.begin() + (stopX - vecX.begin());
+  std::vector<double> peakErrors(startE, stopE);
+
+  double peakSum = std::accumulate(peakValues.cbegin(), peakValues.cend(), 0.0);
+  double sigma = sqrt(std::accumulate(peakErrors.cbegin(), peakErrors.cend(), 0.0, VectorHelper::SumSquares<double>()));
+
+  return peakSum / ((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) {

diff --git a/Framework/Algorithms/src/PDCalibration.cpp b/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 after 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");

diff --git a/Framework/Algorithms/test/FitPeaksTest.h b/Framework/Algorithms/test/FitPeaksTest.h
@@ -1437,6 +1437,99 @@ class FitPeaksTest : public CxxTest::TestSuite {
     API::AnalysisDataService::Instance().remove(fit_window_ws_name);
   }
 
+  //----------------------------------------------------------------------------------------------
+  /** Test that FitPeaks rejects a peak which has a signal-to-sigma ratio below threshold
+   * @brief test_signalToSigmaRatio
+   */
+  void test_signalToSigmaRatio() {
+    g_log.notice() << "TEST SIGNAL TO SIGMA";
+    // create a simple workspace
+    MatrixWorkspace_sptr WS = WorkspaceCreationHelper::create2DWorkspaceWithFullInstrument(
+        static_cast<int>(1 /*num_specs*/), static_cast<int>(400 /*num_data_points*/));
+    WS->getAxis(0)->unit() = Mantid::Kernel::UnitFactory::Instance().create("dSpacing");
+
+    // change the resolution of the binning
+    for (size_t i = 0; i < 1; ++i)
+      WS->mutableX(i) *= 0.05 /*res*/;
+
+    // generate 2 gaussian peaks: a weaker peak at X=5 and a stronger peak at X=10, with some random noise
+    std::random_device rd;
+    std::mt19937 gen(rd());
+    std::normal_distribution<> noise(0, 0.2);
+    const auto &xvals = WS->points(0);
+    std::transform(xvals.cbegin(), xvals.cend(), WS->mutableY(0).begin(), [gen, noise](const double x) mutable {
+      return 20 * exp(-0.5 * pow((x - 10) / 0.1, 2)) + exp(-0.5 * pow((x - 5) / 0.15, 2)) + 1 + noise(gen);
+    });
+
+    // set error values to 2*sqrt(y)
+    const auto &yvals = WS->histogram(0).y();
+    std::transform(yvals.cbegin(), yvals.cend(), WS->mutableE(0).begin(), [](const double y) { return 0.2 * sqrt(y); });
+
+    // set workspace name
+    const std::string data_ws_name("data_s2ns");
+    AnalysisDataService::Instance().addOrReplace(data_ws_name, WS);
+
+    // create peak-center and fit-window workspaces for the 2 peaks generated above (at X=5 and X=10)
+    std::vector<int> peak_index_vec{0, 1};
+    const std::string peak_center_ws_name = genPeakCenterWorkspace(peak_index_vec, "peakcenter_s2nr", 1 /*spectra*/);
+    const std::string fit_window_ws_name =
+        genFitWindowWorkspace(peak_index_vec, "peakwindow_s2nr", 1 /*spectra*/, 1.0 /*fit window halfwidth*/);
+
+    // initialize FitPeaks
+    FitPeaks fitpeaks;
+    fitpeaks.initialize();
+    fitpeaks.setRethrows(true);
+    TS_ASSERT(fitpeaks.isInitialized());
+
+    TS_ASSERT_THROWS_NOTHING(fitpeaks.setProperty("InputWorkspace", data_ws_name));
+    TS_ASSERT_THROWS_NOTHING(fitpeaks.setProperty("StartWorkspaceIndex", 0));
+    TS_ASSERT_THROWS_NOTHING(fitpeaks.setProperty("StopWorkspaceIndex", 1));
+    TS_ASSERT_THROWS_NOTHING(fitpeaks.setProperty("PeakFunction", "Gaussian"));
+    TS_ASSERT_THROWS_NOTHING(fitpeaks.setProperty("PeakCentersWorkspace", peak_center_ws_name));
+    TS_ASSERT_THROWS_NOTHING(fitpeaks.setProperty("FitPeakWindowWorkspace", fit_window_ws_name));
+    TS_ASSERT_THROWS_NOTHING(fitpeaks.setProperty("HighBackground", false));
+    TS_ASSERT_THROWS_NOTHING(fitpeaks.setProperty("MinimumSignalToSigmaRatio", 25.));
+
+    TS_ASSERT_THROWS_NOTHING(fitpeaks.setProperty("OutputWorkspace", "PeakPositionsWS3"));
+    TS_ASSERT_THROWS_NOTHING(fitpeaks.setProperty("OutputPeakParametersWorkspace", "PeakParametersWS3"));
+    TS_ASSERT_THROWS_NOTHING(fitpeaks.setProperty("FittedPeaksWorkspace", "FittedPeaksWS3"));
+    TS_ASSERT_THROWS_NOTHING(fitpeaks.setProperty("MaxFitIterations", 200));
+
+    TS_ASSERT_THROWS_NOTHING(fitpeaks.execute());
+    TS_ASSERT(fitpeaks.isExecuted());
+    if (fitpeaks.isExecuted()) {
+      // check output workspaces
+      TS_ASSERT(API::AnalysisDataService::Instance().doesExist("PeakPositionsWS3"));
+      TS_ASSERT(API::AnalysisDataService::Instance().doesExist("PeakParametersWS3"));
+      TS_ASSERT(API::AnalysisDataService::Instance().doesExist("FittedPeaksWS3"));
+
+      // retrieve fitted parameters
+      API::MatrixWorkspace_sptr peak_params_ws =
+          std::dynamic_pointer_cast<API::MatrixWorkspace>(AnalysisDataService::Instance().retrieve("PeakPositionsWS3"));
+      TS_ASSERT(peak_params_ws);
+      // 1 spectrum
+      TS_ASSERT_EQUALS(peak_params_ws->getNumberHistograms(), 1);
+      // 2 peaks
+      TS_ASSERT_EQUALS(peak_params_ws->histogram(0).x().size(), 2);
+
+      const auto &fitted_positions_0 = peak_params_ws->histogram(0).y();
+      TS_ASSERT_EQUALS(fitted_positions_0.size(), 2); // with 2 peaks to fit
+
+      // When the "signal-to-sigma" check fails, FitPeaks sets the peak position to -4.
+      TS_ASSERT_DELTA(fitted_positions_0[0], -4, 0);
+      TS_ASSERT_DELTA(fitted_positions_0[1], 10.0, 1.E-2);
+
+      // clean algorithm-generated workspaces
+      API::AnalysisDataService::Instance().remove("PeakPositionsWS3");
+      API::AnalysisDataService::Instance().remove("PeakParametersWS3");
+      API::AnalysisDataService::Instance().remove("FittedPeaksWS3");
+    }
+
+    // clean
+    API::AnalysisDataService::Instance().remove(peak_center_ws_name);
+    API::AnalysisDataService::Instance().remove(fit_window_ws_name);
+  }
+
   //--------------------------------------------------------------------------------------------------------------
   /** generate a peak-center workspace compatible with the workspace created by
    * generateTestDataGaussian(), which will have up to 3 spectra up to 2 peaks each

diff --git a/docs/source/release/v6.11.0/Diffraction/Powder/New_features/37826.rst b/docs/source/release/v6.11.0/Diffraction/Powder/New_features/37826.rst
@@ -0,0 +1 @@
+- Add a minimum signal-to-sigma post-fitting check to :ref:`FitPeaks <algm-FitPeaks>` and :ref:`PDCalibration <algm-PDCalibration>`, where peaks with a signal below the provided threshold will be rejected.