src/TRestRawSignal.cxx

/*************************************************************************
 * This file is part of the REST software framework.                     *
 *                                                                       *
 * Copyright (C) 2016 GIFNA/TREX (University of Zaragoza)                *
 * For more information see http://gifna.unizar.es/trex                  *
 *                                                                       *
 * REST is free software: you can redistribute it and/or modify          *
 * it under the terms of the GNU General Public License as published by  *
 * the Free Software Foundation, either version 3 of the License, or     *
 * (at your option) any later version.                                   *
 *                                                                       *
 * REST is distributed in the hope that it will be useful,               *
 * but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the          *
 * GNU General Public License for more details.                          *
 *                                                                       *
 * You should have a copy of the GNU General Public License along with   *
 * REST in $REST_PATH/LICENSE.                                           *
 * If not, see http://www.gnu.org/licenses/.                             *
 * For the list of contributors see $REST_PATH/CREDITS.                  *
 *************************************************************************/
 
#include "TRestRawSignal.h"
 
#include <TAxis.h>
#include <TF1.h>
#include <TMath.h>
#include <TRandom3.h>
 
#include <numeric>
 
using namespace std;
 
ClassImp(TRestRawSignal);
 
TRestRawSignal::TRestRawSignal() {
    fGraph = nullptr;
 
    Initialize();
}
 
TRestRawSignal::TRestRawSignal(Int_t nBins) {
    fGraph = nullptr;
 
    Initialize();
 
    fSignalData.resize(nBins, 0);
}
 
TRestRawSignal::~TRestRawSignal() {}
 
void TRestRawSignal::Initialize() {
    fSignalData.clear();
    fPointsOverThreshold.clear();
    fSignalID = -1;
 
    fThresholdIntegral = -1;
 
    fHeadPoints = 0;
    fTailPoints = 0;
 
    fBaseLine = 0;
    fBaseLineSigma = 0;
}
 
void TRestRawSignal::Reset() {
    Int_t nBins = GetNumberOfPoints();
    Initialize();
    fSignalData.resize(nBins, 0);
}
 
void TRestRawSignal::AddPoint(Short_t value) { fSignalData.push_back(value); }
 
void TRestRawSignal::AddPoint(Double_t value) {
    if (value > numeric_limits<Short_t>::max()) {
        value = numeric_limits<Short_t>::max();
    } else if (value < numeric_limits<Short_t>::min()) {
        value = numeric_limits<Short_t>::min();
    }
    AddPoint((Short_t)value);
}
 
Short_t TRestRawSignal::operator[](Int_t n) {
    if (n >= GetNumberOfPoints()) {
        if (fShowWarnings) {
            std::cout << "TRestRawSignal::GetSignalData: outside limits" << endl;
            std::cout << "Warnings at TRestRawSignal have been disabled" << endl;
            fShowWarnings = false;
        }
        return 0xFFFF;
    }
    return fSignalData[n];
}
 
Double_t TRestRawSignal::GetData(Int_t n) const { return (Double_t)fSignalData[n] - fBaseLine; }
 
Double_t TRestRawSignal::GetRawData(Int_t n) const { return (Double_t)fSignalData[n]; }
 
void TRestRawSignal::IncreaseBinBy(Int_t bin, Double_t data) {
    if (bin >= GetNumberOfPoints()) {
        if (fShowWarnings) {
            std::cout << "TRestRawSignal::IncreaseBinBy: outside limits" << endl;
            std::cout << "Warnings at TRestRawSignal have been disabled" << endl;
            fShowWarnings = false;
        }
 
        return;
    }
 
    fSignalData[bin] += data;
}
 
void TRestRawSignal::InitializePointsOverThreshold(const TVector2& thrPar, Int_t nPointsOver,
                                                   Int_t nPointsFlat) {
    if (fRange.X() < 0) fRange.SetX(0);
    if (fRange.Y() <= 0) fRange.SetY(GetNumberOfPoints());
 
    fPointsOverThreshold.clear();
 
    double pointTh = thrPar.X();
    double signalTh = thrPar.Y();
 
    double threshold = pointTh * fBaseLineSigma;
 
    for (int i = fRange.X(); i < fRange.Y(); i++) {
        // Filling a pulse with consecutive points that are over threshold
        if (this->GetData(i) > threshold) {
            int pos = i;
            std::vector<double> pulse;
            pulse.push_back(this->GetData(i));
            i++;
 
            // If the pulse ends in a flat end above the threshold, the parameter
            // nPointsFlat will serve to artificially end the pulse.
            // If nPointsFlat is big enough, this parameter will not affect the
            // decision to cut this anomalous behaviour. And all points over threshold
            // will be added to the pulse vector.
            int flatN = 0;
            while (i < fRange.Y() && this->GetData(i) > threshold) {
                if (TMath::Abs(this->GetData(i) - this->GetData(i - 1)) > threshold) {
                    flatN = 0;
                } else {
                    flatN++;
                }
 
                if (flatN < nPointsFlat) {
                    pulse.push_back(this->GetData(i));
                    i++;
                } else {
                    break;
                }
            }
 
            if (pulse.size() >= (unsigned int)nPointsOver) {
                // auto stdev = TMath::StdDev(begin(pulse), end(pulse));
                // calculate stdev
                double mean = std::accumulate(pulse.begin(), pulse.end(), 0.0) / pulse.size();
                double sq_sum = std::inner_product(pulse.begin(), pulse.end(), pulse.begin(), 0.0);
                double stdev = std::sqrt(sq_sum / pulse.size() - mean * mean);
 
                if (stdev > signalTh * fBaseLineSigma)
                    for (int j = pos; j < i; j++) fPointsOverThreshold.push_back(j);
            }
        }
    }
 
    CalculateThresholdIntegral();
}
 
void TRestRawSignal::CalculateThresholdIntegral() {
    if (fRange.X() < 0) fRange.SetX(0);
    if (fRange.Y() <= 0 || fRange.Y() > GetNumberOfPoints()) fRange.SetY(GetNumberOfPoints());
 
    fThresholdIntegral = 0;
 
    for (unsigned int n = 0; n < fPointsOverThreshold.size(); n++) {
        if (fPointsOverThreshold[n] >= fRange.X() && fPointsOverThreshold[n] < fRange.Y()) {
            fThresholdIntegral += GetData(fPointsOverThreshold[n]);
        }
    }
}
 
Double_t TRestRawSignal::GetIntegral() {
    if (fRange.X() < 0) fRange.SetX(0);
    if (fRange.Y() <= 0 || fRange.Y() > GetNumberOfPoints()) fRange.SetY(GetNumberOfPoints());
 
    Double_t sum = 0;
    for (int i = fRange.X(); i < fRange.Y(); i++) sum += GetData(i);
    return sum;
}
 
Double_t TRestRawSignal::GetIntegralInRange(Int_t startBin, Int_t endBin) {
    if (startBin < 0) startBin = 0;
    if (endBin <= 0 || endBin > GetNumberOfPoints()) endBin = GetNumberOfPoints();
 
    Double_t sum = 0;
    for (int i = startBin; i < endBin; i++) sum += GetRawData(i);
    return sum;
}
 
Double_t TRestRawSignal::GetThresholdIntegral() {
    if (fThresholdIntegral == -1)
        if (fShowWarnings) {
            std::cout << "TRestRawSignal::GetThresholdIntegral. "
                         "InitializePointsOverThreshold should be "
                         "called first!"
                      << endl;
            fShowWarnings = false;
        }
    return fThresholdIntegral;
}
 
Double_t TRestRawSignal::GetSlopeIntegral() {
    if (fThresholdIntegral == -1)
        cout << "REST Warning. TRestRawSignal::GetSlopeIntegral. "
                "InitializePointsOverThreshold should be called first."
             << endl;
 
    Double_t sum = 0;
    for (unsigned int n = 0; n < fPointsOverThreshold.size(); n++) {
        if (n + 1 >= fPointsOverThreshold.size()) return sum;
 
        sum += GetData(fPointsOverThreshold[n]);
 
        if (GetData(fPointsOverThreshold[n + 1]) - GetData(fPointsOverThreshold[n]) < 0) break;
    }
    return sum;
}
 
Double_t TRestRawSignal::GetRiseSlope() {
    if (fThresholdIntegral == -1)
        cout << "REST Warning. TRestRawSignal::GetRiseSlope. "
                "InitializePointsOverThreshold should be called first."
             << endl;
 
    if (fPointsOverThreshold.size() < 2) {
        // cout << "REST Warning. TRestRawSignal::GetRiseSlope. Less than 2 points!." << endl;
        return 0;
    }
 
    Int_t maxBin = GetMaxPeakBin() - 1;
 
    Double_t hP = GetData(maxBin);
 
    Double_t lP = GetData(fPointsOverThreshold[0]);
 
    return (hP - lP) / (maxBin - fPointsOverThreshold[0]);
}
 
Int_t TRestRawSignal::GetRiseTime() {
    if (fThresholdIntegral == -1)
        cout << "REST Warning. TRestRawSignal::GetRiseTime. "
                "InitializePointsOverThreshold should be called first."
             << endl;
 
    if (fPointsOverThreshold.size() < 2) {
        // cout << "REST Warning. TRestRawSignal::GetRiseTime. Less than 2 points!." << endl;
        return 0;
    }
 
    return GetMaxPeakBin() - fPointsOverThreshold[0];
}
 
Double_t TRestRawSignal::GetTripleMaxIntegral() {
    if (fRange.X() < 0) fRange.SetX(0);
    if (fRange.Y() <= 0 || fRange.Y() > GetNumberOfPoints()) fRange.SetY(GetNumberOfPoints());
 
    if (fThresholdIntegral == -1) {
        cout << "REST Warning. TRestRawSignal::GetTripleMaxIntegral. "
                "InitializePointsOverThreshold should be called first."
             << endl;
        return 0;
    }
 
    if (fPointsOverThreshold.size() < 2) {
        // cout << "REST Warning. TRestRawSignal::GetTripleMaxIntegral. Points over
        // "
        //        "threshold = "
        //     << fPointsOverThreshold.size() << endl;
        return 0;
    }
 
    Int_t cBin = GetMaxPeakBin();
 
    if (cBin + 1 >= GetNumberOfPoints()) return 0;
 
    Double_t a1 = GetData(cBin);
    Double_t a2 = GetData(cBin - 1);
    Double_t a3 = GetData(cBin + 1);
 
    return a1 + a2 + a3;
}
 
Double_t TRestRawSignal::GetAverageInRange(Int_t startBin, Int_t endBin) {
    if (startBin < 0) startBin = 0;
    if (endBin <= 0 || endBin > GetNumberOfPoints()) endBin = GetNumberOfPoints();
 
    Double_t sum = 0;
    for (int i = startBin; i <= endBin; i++) sum += this->GetData(i);
 
    return sum / (endBin - startBin + 1);
}
 
Int_t TRestRawSignal::GetMaxPeakWidth() {
    Int_t mIndex = this->GetMaxPeakBin();
    Double_t maxValue = this->GetData(mIndex);
 
    Double_t value = maxValue;
    Int_t rightIndex = mIndex;
    while (value > maxValue / 2) {
        value = this->GetData(rightIndex);
        rightIndex++;
    }
    Int_t leftIndex = mIndex;
    value = maxValue;
    while (value > maxValue / 2) {
        value = this->GetData(leftIndex);
        leftIndex--;
    }
 
    return rightIndex - leftIndex;
}
 
Double_t TRestRawSignal::GetMaxPeakValue() { return GetData(GetMaxPeakBin()); }
 
Int_t TRestRawSignal::GetMaxPeakBin() {
    Double_t max = numeric_limits<Double_t>::min();
 
    Int_t index = 0;
 
    if (fRange.Y() == 0 || fRange.Y() > GetNumberOfPoints()) fRange.SetY(GetNumberOfPoints());
    if (fRange.X() < 0) fRange.SetX(0);
 
    for (int i = fRange.X(); i < fRange.Y(); i++) {
        if (this->GetData(i) > max) {
            max = GetData(i);
            index = i;
        }
    }
 
    return index;
}
 
Double_t TRestRawSignal::GetMinPeakValue() { return GetData(GetMinPeakBin()); }
 
Int_t TRestRawSignal::GetMinPeakBin() {
    Double_t min = numeric_limits<Double_t>::max();
    Int_t index = 0;
 
    if (fRange.Y() == 0 || fRange.Y() > GetNumberOfPoints()) fRange.SetY(GetNumberOfPoints());
    if (fRange.X() < 0) fRange.SetX(0);
 
    for (int i = fRange.X(); i < fRange.Y(); i++) {
        if (this->GetData(i) < min) {
            min = GetData(i);
            index = i;
        }
    }
 
    return index;
}
 
Bool_t TRestRawSignal::IsADCSaturation(int Nflat) {
    if (Nflat <= 0) return false;
    // GetMaxPeakBin() will always find the first max peak bin if multiple bins are in same max value.
    int index = GetMaxPeakBin();
    Short_t value = fSignalData[index];
 
    bool sat = false;
    if (index + Nflat <= (int)fSignalData.size()) {
        for (int i = index; i < index + Nflat; i++) {
            if (fSignalData[i] != value) {
                break;
            }
            if (i == index + Nflat - 1) {
                sat = true;
            }
        }
    }
 
    return sat;
}
 
void TRestRawSignal::GetDifferentialSignal(TRestRawSignal* diffSignal, Int_t smearPoints) {
    if (smearPoints <= 0) smearPoints = 1;
    diffSignal->Initialize();
 
    for (int i = 0; i < smearPoints; i++) {
        diffSignal->AddPoint((Short_t)0);
    }
 
    for (int i = smearPoints; i < this->GetNumberOfPoints() - smearPoints; i++) {
        Double_t value = 0.5 * (this->GetData(i + smearPoints) - GetData(i - smearPoints)) / smearPoints;
 
        diffSignal->AddPoint((Short_t)value);
    }
 
    for (int i = GetNumberOfPoints() - smearPoints; i < GetNumberOfPoints(); i++) {
        diffSignal->AddPoint((Short_t)0);
    }
}
 
void TRestRawSignal::GetWhiteNoiseSignal(TRestRawSignal* noiseSignal, Double_t noiseLevel) {
    TRandom3 random(fSeed);
 
    for (int i = 0; i < GetNumberOfPoints(); i++) {
        Double_t value = this->GetData(i) + random.Gaus(0, noiseLevel);
        // do not cast as short so that there are no problems with overflows
        // (https://github.com/rest-for-physics/rawlib/issues/113)
        noiseSignal->AddPoint(value);
    }
}
 
void TRestRawSignal::GetSignalSmoothed(TRestRawSignal* smoothedSignal, Int_t averagingPoints) {
    smoothedSignal->Initialize();
 
    averagingPoints = (averagingPoints / 2) * 2 + 1;  // make it odd >= averagingPoints
 
    Double_t sumAvg = GetIntegralInRange(0, averagingPoints) / averagingPoints;
 
    for (int i = 0; i <= averagingPoints / 2; i++) smoothedSignal->AddPoint((Short_t)sumAvg);
 
    for (int i = averagingPoints / 2 + 1; i < GetNumberOfPoints() - averagingPoints / 2; i++) {
        sumAvg -= this->GetRawData(i - (averagingPoints / 2 + 1)) / averagingPoints;
        sumAvg += this->GetRawData(i + averagingPoints / 2) / averagingPoints;
        smoothedSignal->AddPoint((Short_t)sumAvg);
    }
 
    for (int i = GetNumberOfPoints() - averagingPoints / 2; i < GetNumberOfPoints(); i++)
        smoothedSignal->AddPoint(sumAvg);
}
 
std::vector<Float_t> TRestRawSignal::GetSignalSmoothed(Int_t averagingPoints, std::string option) {
    std::vector<Float_t> result;
 
    if (option == "") {
        result.resize(GetNumberOfPoints());
 
        averagingPoints = (averagingPoints / 2) * 2 + 1;  // make it odd >= averagingPoints
 
        Float_t sumAvg = (Float_t)GetIntegralInRange(0, averagingPoints) / averagingPoints;
 
        for (int i = 0; i <= averagingPoints / 2; i++) result[i] = sumAvg;
 
        for (int i = averagingPoints / 2 + 1; i < GetNumberOfPoints() - averagingPoints / 2; i++) {
            sumAvg -= this->GetRawData(i - (averagingPoints / 2 + 1)) / averagingPoints;
            sumAvg += this->GetRawData(i + averagingPoints / 2) / averagingPoints;
            result[i] = sumAvg;
        }
 
        for (int i = GetNumberOfPoints() - averagingPoints / 2; i < GetNumberOfPoints(); i++)
            result[i] = sumAvg;
    } else if (ToUpper(option) == "EXCLUDE OUTLIERS") {
        result = GetSignalSmoothed_ExcludeOutliers(averagingPoints);
    } else {
        cout << "TRestRawSignal::GetSignalSmoothed. Error! No such option!" << endl;
    }
    return result;
}
 
std::vector<Float_t> TRestRawSignal::GetSignalSmoothed_ExcludeOutliers(Int_t averagingPoints) {
    std::vector<Float_t> result(GetNumberOfPoints());
 
    if (fBaseLine == 0) CalculateBaseLine(5, GetNumberOfPoints() - 5, "ROBUST");
 
    averagingPoints = (averagingPoints / 2) * 2 + 1;  // make it odd >= averagingPoints
 
    Float_t sumAvg = (Float_t)GetIntegralInRange(0, averagingPoints) / averagingPoints;
 
    // Points at the beginning, where we can calculate a moving average
    for (int i = 0; i <= averagingPoints / 2; i++) result[i] = sumAvg;
 
    // Points in the middle
    float_t amplitude;
    for (int i = averagingPoints / 2 + 1; i < GetNumberOfPoints() - averagingPoints / 2; i++) {
        amplitude = this->GetRawData(i - (averagingPoints / 2 + 1));
        sumAvg -= (std::abs(amplitude - fBaseLine) > 3 * fBaseLineSigma) ? fBaseLine / averagingPoints
                                                                         : amplitude / averagingPoints;
        amplitude = this->GetRawData(i + averagingPoints / 2);
        sumAvg += (std::abs(amplitude - fBaseLine) > 3 * fBaseLineSigma) ? fBaseLine / averagingPoints
                                                                         : amplitude / averagingPoints;
        result[i] = sumAvg;
    }
 
    // Points at the end, where we can calculate a moving average
    for (int i = GetNumberOfPoints() - averagingPoints / 2; i < GetNumberOfPoints(); i++) result[i] = sumAvg;
    return result;
}
 
void TRestRawSignal::GetBaseLineCorrected(TRestRawSignal* smoothedSignal, Int_t averagingPoints) {
    smoothedSignal->Initialize();
 
    std::vector<Float_t> averagedSignal = GetSignalSmoothed(averagingPoints, "EXCLUDE OUTLIERS");
 
    for (int i = 0; i < GetNumberOfPoints(); i++) {
        smoothedSignal->AddPoint(GetRawData(i) - averagedSignal[i]);
    }
}
 
void TRestRawSignal::CalculateBaseLineMean(Int_t startBin, Int_t endBin) {
    if (endBin - startBin <= 0) {
        fBaseLine = 0.;
    } else if (endBin > (int)fSignalData.size()) {
        cout << "TRestRawSignal::CalculateBaseLine. Error! Baseline range exceeds the rawdata depth!!"
             << endl;
        endBin = fSignalData.size();
    } else {
        Double_t baseLine = 0;
        for (int i = startBin; i < endBin; i++) baseLine += fSignalData[i];
        fBaseLine = baseLine / (endBin - startBin);
    }
}
 
void TRestRawSignal::CalculateBaseLineMedian(Int_t startBin, Int_t endBin) {
    if (endBin - startBin <= 0) {
        fBaseLine = 0.;
    } else if (endBin > (int)fSignalData.size()) {
        cout << "TRestRawSignal::CalculateBaseLine. Error! Baseline range exceeds the rawdata depth!!"
             << endl;
        endBin = fSignalData.size();
    } else {
        vector<Short_t>::const_iterator first = fSignalData.begin() + startBin;
        vector<Short_t>::const_iterator last = fSignalData.begin() + endBin;
        vector<Short_t> v(first, last);
        const Short_t* signalInRange = &v[0];
        fBaseLine = TMath::Median(endBin - startBin, signalInRange);
    }
}
 
void TRestRawSignal::CalculateBaseLine(Int_t startBin, Int_t endBin, const std::string& option) {
    if (ToUpper(option) == "ROBUST") {
        CalculateBaseLineMedian(startBin, endBin);
        CalculateBaseLineSigmaIQR(startBin, endBin);
    } else {
        CalculateBaseLineMean(startBin, endBin);
        CalculateBaseLineSigmaSD(startBin, endBin);
    }
}
 
void TRestRawSignal::CalculateBaseLineSigmaSD(Int_t startBin, Int_t endBin) {
    if (endBin - startBin <= 0) {
        fBaseLineSigma = 0;
    } else {
        Double_t baseLineSigma = 0;
        for (int i = startBin; i < endBin; i++)
            baseLineSigma += (fBaseLine - fSignalData[i]) * (fBaseLine - fSignalData[i]);
        fBaseLineSigma = TMath::Sqrt(baseLineSigma / (endBin - startBin));
    }
}
 
void TRestRawSignal::CalculateBaseLineSigmaIQR(Int_t startBin, Int_t endBin) {
    if (endBin - startBin <= 0) {
        fBaseLineSigma = 0;
    } else {
        vector<Short_t>::const_iterator first = fSignalData.begin() + startBin;
        vector<Short_t>::const_iterator last = fSignalData.begin() + endBin;
        vector<Short_t> v(first, last);
        std::sort(v.begin(), v.end());
        Short_t Q1 = v[(int)(endBin - startBin) * 0.25];
        Short_t Q3 = v[(int)(endBin - startBin) * 0.75];
        Double_t IQR = Q3 - Q1;
        fBaseLineSigma =
            IQR / 1.349;  // IQR/1.349 equals the standard deviation in case of normally distributed data
    }
}
 
void TRestRawSignal::AddOffset(Short_t offset) {
    if (fBaseLine != 0 || fBaseLineSigma != 0) fBaseLineSigma += (Double_t)offset;
    for (int i = 0; i < GetNumberOfPoints(); i++) fSignalData[i] = fSignalData[i] + offset;
}
 
void TRestRawSignal::Scale(Double_t value) {
    for (int i = 0; i < GetNumberOfPoints(); i++) {
        Double_t scaledValue = value * fSignalData[i];
        fSignalData[i] = (Short_t)scaledValue;
    }
}
 
void TRestRawSignal::SignalAddition(const TRestRawSignal& signal) {
    if (this->GetNumberOfPoints() != signal.GetNumberOfPoints()) {
        cout << "ERROR : TRestRawSignal::SignalAddition." << endl;
        cout << "I cannot add two signals with different number of points" << endl;
        return;
    }
 
    for (int i = 0; i < GetNumberOfPoints(); i++) {
        fSignalData[i] += signal.GetData(i);
    }
}
 
void TRestRawSignal::WriteSignalToTextFile(const TString& filename) {
    // We should check it is writable
    FILE* file = fopen(filename.Data(), "w");
    for (int i = 0; i < GetNumberOfPoints(); i++) fprintf(file, "%d\t%lf\n", i, GetData(i));
    fclose(file);
}
 
void TRestRawSignal::Print() const {
    cout << "---------------------" << endl;
    cout << "Signal id : " << this->GetSignalID() << endl;
    cout << "Baseline : " << fBaseLine << endl;
    cout << "Baseline sigma : " << fBaseLineSigma << endl;
    cout << "+++++++++++++++++++++" << endl;
    for (int i = 0; i < GetNumberOfPoints(); i++)
        cout << "Bin : " << i << " amplitude : " << GetRawData(i) << endl;
    cout << "---------------------" << endl;
}
 
TGraph* TRestRawSignal::GetGraph(Int_t color) {
    delete fGraph;
    fGraph = new TGraph();
 
    fGraph->SetLineWidth(2);
    fGraph->SetLineColor(color % 8 + 1);
    fGraph->SetMarkerStyle(7);
 
    for (int i = 0; i < GetNumberOfPoints(); i++) {
        fGraph->SetPoint(i, i, GetData(i));
    }
 
    fGraph->GetXaxis()->SetLimits(0, GetNumberOfPoints() - 1);
 
    /*
     * To draw x axis in multiples of 2
    for (int i = 0; i < values.size(); i++) {
        if (i % 32 != 0 && i != values.size() - 1){
            continue;
        }
        fGraph->GetXaxis()->SetBinLabel(fGraph->GetXaxis()->FindBin(i), std::to_string(i).c_str());
    }
     */
 
    return fGraph;
}
 
vector<pair<UShort_t, double>> TRestRawSignal::GetPeaks(double threshold, UShort_t distance) const {
    vector<pair<UShort_t, double>> peaks;
 
    for (UShort_t i = 0; i < GetNumberOfPoints(); i++) {
        const double point = GetRawData(i);
        if (i > 0 && i < GetNumberOfPoints() - 1) {
            double prevPoint = GetRawData(i - 1);
            double nextPoint = GetRawData(i + 1);
 
            if (point > threshold && point >= prevPoint && point >= nextPoint) {
                // Check if the peak is spaced far enough from the previous peak
                if (peaks.empty() || i - peaks.back().first >= distance) {
                    peaks.push_back(std::make_pair(i, point));
                }
            }
        }
    }
    return peaks;
}