Add test and release note for signal-to-sigma check

mantidproject · Sep 10, 2024 · 46efb92 · 46efb92
1 parent 2f196ea
commit 46efb92
Show file tree

Hide file tree

Showing 4 changed files with 101 additions and 16 deletions.
diff --git a/Framework/Algorithms/inc/MantidAlgorithms/FitPeaks.h b/Framework/Algorithms/inc/MantidAlgorithms/FitPeaks.h
@@ -217,8 +217,7 @@ class MANTID_ALGORITHMS_DLL FitPeaks final : public API::Algorithm {
   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);
+                                     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);

diff --git a/Framework/Algorithms/src/FitPeaks.cpp b/Framework/Algorithms/src/FitPeaks.cpp
@@ -1295,7 +1295,7 @@ void FitPeaks::fitSpectrumPeaks(size_t wi, const std::vector<double> &expected_p
         fit_result->setBadRecord(peak_index, -1.);
 
       if (m_minSignalToSigmaRatio > 0) {
-        if (calculateSignalToSigmaRatio(wi, peak_window_i, peakfunction, bkgdfunction); < m_minSignalToSigmaRatio) {
+        if (calculateSignalToSigmaRatio(wi, peak_window_i, peakfunction) < m_minSignalToSigmaRatio) {
           fit_result->setBadRecord(peak_index, -1.);
           cost = DBL_MAX;
         }
@@ -1552,34 +1552,26 @@ void FitPeaks::calculateFittedPeaks(const std::vector<std::shared_ptr<FitPeaksAl
 }
 
 double FitPeaks::calculateSignalToSigmaRatio(const size_t &iws, const std::pair<double, double> &peakWindow,
-                                             const API::IPeakFunction_sptr &peakFunction,
-                                             const API::IBackgroundFunction_sptr &backgroundFunction) {
+                                             const API::IPeakFunction_sptr &peakFunction) {
   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();
+  peakFunction->function(domain, values);
+  auto peakValues = 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 peakSum = std::accumulate(peakValues.cbegin(), peakValues.cend(), 0.0);
   double sigma = sqrt(std::accumulate(peakErrors.cbegin(), peakErrors.cend(), 0.0, VectorHelper::SumSquares<double>()));
 
-  return (fittedSum - bkgdSum) / ((sigma == 0) ? 1 : sigma);
+  return peakSum / ((sigma == 0) ? 1 : sigma);
 }
 
 namespace {

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", 50.));
+
+    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.