From 2f196eacf570faf2fabd6718c18f2c3390da7c4f Mon Sep 17 00:00:00 2001 From: Mohammed Almakki Date: Mon, 19 Aug 2024 14:22:19 +0100 Subject: [PATCH] add signal to sigma post fitting check --- .../inc/MantidAlgorithms/FitPeaks.h | 6 +++ Framework/Algorithms/src/FitPeaks.cpp | 51 +++++++++++++++++++ Framework/Algorithms/src/PDCalibration.cpp | 9 +++- 3 files changed, 65 insertions(+), 1 deletion(-) diff --git a/Framework/Algorithms/inc/MantidAlgorithms/FitPeaks.h b/Framework/Algorithms/inc/MantidAlgorithms/FitPeaks.h index c9fb198fd7bf..dd09f237b79c 100644 --- a/Framework/Algorithms/inc/MantidAlgorithms/FitPeaks.h +++ b/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> &fit_results); + double calculateSignalToSigmaRatio(const size_t &iws, const std::pair &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; diff --git a/Framework/Algorithms/src/FitPeaks.cpp b/Framework/Algorithms/src/FitPeaks.cpp index 44263b85af71..1364112fdb56 100644 --- 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 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 &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 &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(); + + 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 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())); + + return (fittedSum - bkgdSum) / ((sigma == 0) ? 1 : sigma); +} + namespace { bool estimateBackgroundParameters(const Histogram &histogram, const std::pair &peak_window, const API::IBackgroundFunction_sptr &bkgd_function) { diff --git a/Framework/Algorithms/src/PDCalibration.cpp b/Framework/Algorithms/src/PDCalibration.cpp index bafa5308d7fa..071728dfce38 100644 --- 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 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");