From 46efb92c8a0a818df26628afd614e3003aad3c29 Mon Sep 17 00:00:00 2001 From: Mohammed Almakki Date: Tue, 10 Sep 2024 09:20:42 +0100 Subject: [PATCH] Add test and release note for signal-to-sigma check --- .../inc/MantidAlgorithms/FitPeaks.h | 3 +- Framework/Algorithms/src/FitPeaks.cpp | 20 ++-- Framework/Algorithms/test/FitPeaksTest.h | 93 +++++++++++++++++++ .../Diffraction/Powder/New_features/37826.rst | 1 + 4 files changed, 101 insertions(+), 16 deletions(-) create mode 100644 docs/source/release/v6.11.0/Diffraction/Powder/New_features/37826.rst diff --git a/Framework/Algorithms/inc/MantidAlgorithms/FitPeaks.h b/Framework/Algorithms/inc/MantidAlgorithms/FitPeaks.h index dd09f237b79c..93d0c527dc7e 100644 --- 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> &fit_results); double calculateSignalToSigmaRatio(const size_t &iws, const std::pair &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 index 1364112fdb56..8b721c7dc55c 100644 --- 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 &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 &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(); - 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 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())); - 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 index 95218ad8dc54..9b6030de86c3 100644 --- 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(1 /*num_specs*/), static_cast(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 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(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 new file mode 100644 index 000000000000..f1785a1f8500 --- /dev/null +++ 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 ` and :ref:`PDCalibration `, where peaks with a signal below the provided threshold will be rejected. \ No newline at end of file