From a4ba0c1c6cdcf62638212cdd5489206532372a2c Mon Sep 17 00:00:00 2001 From: JPorron <114648351+JPorron@users.noreply.github.com> Date: Mon, 11 Nov 2024 12:43:33 +0100 Subject: [PATCH] Calculate veto signals baseline without outliers (#146) * Calculate veto signals baseline without outliers * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --------- Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> --- inc/TRestRawSignal.h | 4 ++ src/TRestRawPeaksFinderProcess.cxx | 2 +- src/TRestRawSignal.cxx | 98 +++++++++++++++++++++++++++++- 3 files changed, 102 insertions(+), 2 deletions(-) diff --git a/inc/TRestRawSignal.h b/inc/TRestRawSignal.h index 501cc3da..dc390b0f 100644 --- a/inc/TRestRawSignal.h +++ b/inc/TRestRawSignal.h @@ -41,6 +41,8 @@ class TRestRawSignal { void CalculateBaseLineSigmaIQR(Int_t startBin, Int_t endBin); + void CalculateBaseLineSigmaExcludeOutliers(Int_t startBin, Int_t endBin); + std::vector GetSignalSmoothed_ExcludeOutliers(Int_t averagingPoints); protected: @@ -198,6 +200,8 @@ class TRestRawSignal { void CalculateBaseLineMedian(Int_t startBin, Int_t endBin); + void CalculateBaseLineMedianExcludeOutliers(Int_t startBin, Int_t endBin); + void CalculateBaseLine(Int_t startBin, Int_t endBin, const std::string& option = ""); void GetBaseLineCorrected(TRestRawSignal* smoothedSignal, Int_t averagingPoints); diff --git a/src/TRestRawPeaksFinderProcess.cxx b/src/TRestRawPeaksFinderProcess.cxx index 92c2107a..239e77a4 100644 --- a/src/TRestRawPeaksFinderProcess.cxx +++ b/src/TRestRawPeaksFinderProcess.cxx @@ -97,7 +97,7 @@ TRestEvent* TRestRawPeaksFinderProcess::ProcessEvent(TRestEvent* inputEvent) { } else if (channelType == "veto") { // For veto signals the baseline is calculated over the whole range, as we don´t know where the // signal will be. - signal->CalculateBaseLine(0, 511, "robust"); + signal->CalculateBaseLine(0, 511, "OUTLIERS"); // For veto signals the threshold is selected by the user. const auto peaks = signal->GetPeaksVeto(signal->GetBaseLine() + fThresholdOverBaseline, fDistance); diff --git a/src/TRestRawSignal.cxx b/src/TRestRawSignal.cxx index 72747e3d..ac813f63 100644 --- a/src/TRestRawSignal.cxx +++ b/src/TRestRawSignal.cxx @@ -784,6 +784,50 @@ void TRestRawSignal::CalculateBaseLineMedian(Int_t startBin, Int_t endBin) { } } +/////////////////////////////////////////////// +/// \brief This method is called by CalculateBaseLine with the "OUTLIERS"-option and is used to determine the +/// value of the baseline as the median of the data points found in the range defined between startBin and +/// endBin after excluding the outliers. The median is calculated using only the values in the 25-75% removing +/// big and small outliers. +/// +void TRestRawSignal::CalculateBaseLineMedianExcludeOutliers(Int_t startBin, Int_t endBin) { + if (endBin - startBin <= 0) { + fBaseLine = 0.; + return; + } else if (endBin > static_cast(fSignalData.size())) { + cout << "TRestRawSignal::CalculateBaseLine. Error! Baseline range exceeds the rawdata depth!!" + << endl; + endBin = fSignalData.size(); + } else { + // Extract the data within the interval + std::vector data(fSignalData.begin() + startBin, fSignalData.begin() + endBin); + std::sort(data.begin(), data.end()); + + // Calculate Q1 and Q3 for IQR + size_t dataSize = data.size(); + Short_t Q1 = data[dataSize / 4]; + Short_t Q3 = data[(3 * dataSize) / 4]; + Double_t lowerBound = Q1; + Double_t upperBound = Q3; + + // Filter out the outliers + std::vector filteredData; + for (const auto& value : data) { + if (value >= lowerBound && value <= upperBound) { + filteredData.emplace_back(value); + } + } + + // Calculate median of filtered data + if (filteredData.empty()) { + fBaseLine = TMath::Median(data.size(), + &data[0]); // Fall back to original median if all values are outliers + } else { + fBaseLine = TMath::Median(filteredData.size(), &filteredData[0]); + } + } +} + /////////////////////////////////////////////// /// \brief This method calculates the average and fluctuation of the baseline in the /// specified range and writes the values to fBaseLine and fBaseLineSigma respectively. @@ -792,12 +836,16 @@ void TRestRawSignal::CalculateBaseLineMedian(Int_t startBin, Int_t endBin) { /// /// \param option By setting this option to "ROBUST", the average is calculated as median, /// and the fluctuation as interquartile range (IQR), which are less affected by outliers (e.g. a signal -/// pulse). +/// pulse). By setting it to "OUTLIERS" the median and sigma will only take into account the +/// 25-75% values of the interval. /// void TRestRawSignal::CalculateBaseLine(Int_t startBin, Int_t endBin, const std::string& option) { if (ToUpper(option) == "ROBUST") { CalculateBaseLineMedian(startBin, endBin); CalculateBaseLineSigmaIQR(startBin, endBin); + } else if (ToUpper(option) == "OUTLIERS") { + CalculateBaseLineMedianExcludeOutliers(startBin, endBin); + CalculateBaseLineSigmaExcludeOutliers(startBin, endBin); } else { CalculateBaseLineMean(startBin, endBin); CalculateBaseLineSigmaSD(startBin, endBin); @@ -842,6 +890,54 @@ void TRestRawSignal::CalculateBaseLineSigmaIQR(Int_t startBin, Int_t endBin) { } } +/////////////////////////////////////////////// +/// \brief This method is called by CalculateBaseLine with the "OUTLIERS"-option to +/// determine the value of the baseline +/// fluctuation as the standard deviation in the baseline range provided excluding outliers. +/// Since outliers are strongly suppressed there is no need to use the IQR. +/// +void TRestRawSignal::CalculateBaseLineSigmaExcludeOutliers(Int_t startBin, Int_t endBin) { + if (endBin - startBin <= 0) { + fBaseLineSigma = 0.; + return; + } else if (endBin > static_cast(fSignalData.size())) { + cout << "TRestRawSignal::CalculateBaseLineSigma. Error! Range exceeds the rawdata depth!!" << endl; + endBin = fSignalData.size(); + } else { + // Extract the data within the interval + std::vector data(fSignalData.begin() + startBin, fSignalData.begin() + endBin); + std::sort(data.begin(), data.end()); + + // Calculate Q1 and Q3 for IQR + size_t dataSize = data.size(); + Short_t Q1 = data[dataSize / 4]; + Short_t Q3 = data[(3 * dataSize) / 4]; + Double_t lowerBound = Q1; + Double_t upperBound = Q3; + + // Filter out the outliers + std::vector filteredData; + for (const auto& value : data) { + if (value >= lowerBound && value <= upperBound) { + filteredData.emplace_back(value); + } + } + + // Calculate standard deviation of filtered data + if (filteredData.empty()) { + fBaseLineSigma = 0.; // If all values are outliers, set sigma to zero + } else { + double mean = + std::accumulate(filteredData.begin(), filteredData.end(), 0.0) / filteredData.size(); + double variance = 0.0; + for (const auto& value : filteredData) { + variance += std::pow(value - mean, 2); + } + fBaseLineSigma = std::sqrt(variance / filteredData.size()); // Standard deviation + } + } +} + /////////////////////////////////////////////// /// \brief This method adds an offset to the signal data ///