From 50142512c0d88ad8e56e6c045911626834cf09dc Mon Sep 17 00:00:00 2001 From: juanan Date: Wed, 15 Feb 2023 18:35:24 +0100 Subject: [PATCH 01/21] New namespace for a generic signal analysis using templates --- .../analysis/inc/TRestSignalAnalysis.h | 267 ++++++++++++++++++ 1 file changed, 267 insertions(+) create mode 100644 source/framework/analysis/inc/TRestSignalAnalysis.h diff --git a/source/framework/analysis/inc/TRestSignalAnalysis.h b/source/framework/analysis/inc/TRestSignalAnalysis.h new file mode 100644 index 000000000..8d0883e26 --- /dev/null +++ b/source/framework/analysis/inc/TRestSignalAnalysis.h @@ -0,0 +1,267 @@ +/************************************************************************* + * 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 +#include +#include +#include + +#ifndef RestCore_TRestSignalAnalysis +#define RestCore_TRestSignalAnalysis + +#include +#include + +/// This namespace define utilities (functions) to calculate different signal parameters +namespace TRestSignalAnalysis { + +/// Generic functions for different calculations + +template +void CalculateBaselineAndSigmaSD(const std::vector& signal, Int_t startBin, Int_t endBin, + Double_t& baseLine, Double_t& baseLineSigma) { + baseLine = 0; + baseLineSigma = 0; + + int nPoints = 0; + + for (auto i = startBin; i < endBin; i++) { + if (i < 0 || i > signal.size()) continue; + baseLine += signal[i]; + baseLineSigma += signal[i] * signal[i]; + nPoints++; + } + + if (nPoints > 0) { + baseLine /= (double)nPoints; + baseLineSigma = TMath::Sqrt(baseLineSigma / (double)nPoints - baseLine * baseLine); + } +} +void CalculateBaselineAndSigmaSD(const std::vector& signal, Int_t startBin, Int_t endBin, + Double_t& baseLine, Double_t& baseLineSigma); +void CalculateBaselineAndSigmaSD(const std::vector& signal, Int_t startBin, Int_t endBin, + Double_t& baseLine, Double_t& baseLineSigma); + +template +void CalculateBaselineAndSigmaIQR(const std::vector& signal, Int_t startBin, Int_t endBin, + Double_t& baseLine, Double_t& baseLineSigma) { + baseLine = 0; + baseLineSigma = 0; + + if (startBin < 0) startBin = 0; + if (endBin >= signal.size()) endBin = signal.size() - 1; + + if (endBin > startBin) return; + + auto first = signal.begin() + startBin; + auto last = signal.begin() + endBin; + std::vector v(first, last); + baseLine = TMath::Median(endBin - startBin, &v[0]); + + std::sort(v.begin(), v.end()); + baseLineSigma = + v[(int)(endBin - startBin) * 0.75] - + v[(int)(endBin - startBin) * 0.25] / + 1.349; // IQR/1.349 equals the standard deviation in case of normally distributed data +} + +void CalculateBaselineAndSigmaIQR(const std::vector& signal, Int_t startBin, Int_t endBin, + Double_t& baseLine, Double_t& baseLineSigma); +void CalculateBaselineAndSigmaIQR(const std::vector& signal, Int_t startBin, Int_t endBin, + Double_t& baseLine, Double_t& baseLineSigma); + +template +Double_t GetAverage(const std::vector& signal, Int_t startBin, Int_t endBin) { + int nPoints = 0; + Double_t avg = 0; + + for (auto i = startBin; i < endBin; i++) { + if (i < 0 || i > signal.size()) continue; + avg += signal[i]; + nPoints++; + } + + if (nPoints > 0) avg /= (double)nPoints; + + return avg; +} + +Double_t GetAverage(const std::vector& signal, Int_t startBin, Int_t endBin); +Double_t GetAverage(const std::vector& signal, Int_t startBin, Int_t endBin); + +/////////////////////////////////////////////// +/// \brief Return smoothing of signal +/// +/// +/// \param neighbours It defines the number of neighbours +/// points used to average the signal +/// +/// \param option If the option is set to "EXCLUDE OUTLIERS", points that are too far away from the median +/// baseline will be ignored to improve the smoothing result +/// +template +std::vector GetSignalSmoothed(const std::vector& signal, int averagingPoints = 3) { + const size_t pulseDepth = signal.size(); + std::vector smoothed(pulseDepth, 0); + + averagingPoints = (averagingPoints / 2) * 2 + 1; // make it odd >= averagingPoints + + Float_t sumAvg = GetAverage(signal, 0, averagingPoints); + + // Points at the beginning, where we can calculate a moving average + for (int i = 0; i <= averagingPoints / 2; i++) smoothed[i] = sumAvg; + + for (int i = averagingPoints / 2 + 1; i < pulseDepth - averagingPoints / 2; i++) { + int index = i - (averagingPoints / 2 + 1); + sumAvg -= signal[index] / averagingPoints; + index = i + (averagingPoints / 2 + 1); + sumAvg += signal[index] / averagingPoints; + smoothed[i] = sumAvg; + } + + // Points at the end, where we can calculate a moving average + for (int i = pulseDepth - averagingPoints / 2; i < pulseDepth; i++) smoothed[i] = sumAvg; + + return smoothed; +} +std::vector GetSignalSmoothed(const std::vector& signal, int averagingPoints = 3); +std::vector GetSignalSmoothed(const std::vector& signal, int averagingPoints = 3); + +template +std::vector GetSignalSmoothed_ExcludeOutliers(const std::vector& signal, int averagingPoints, + Double_t& baseLine, Double_t& baseLineSigma) { + const size_t pulseDepth = signal.size(); + std::vector smoothed(pulseDepth, 0); + + if (baseLine == 0) CalculateBaselineAndSigmaIQR(signal, 5, int(pulseDepth - 5), baseLine, baseLineSigma); + + averagingPoints = (averagingPoints / 2) * 2 + 1; // make it odd >= averagingPoints + + Float_t sumAvg = GetAverage(signal, 0, averagingPoints); + + // Points at the beginning, where we can calculate a moving average + for (int i = 0; i <= averagingPoints / 2; i++) smoothed[i] = sumAvg; + + // Points in the middle + float_t amplitude; + for (int i = averagingPoints / 2 + 1; i < pulseDepth - averagingPoints / 2; i++) { + int index = i - (averagingPoints / 2 + 1); + amplitude = signal[index]; + sumAvg -= (std::abs(amplitude - baseLine) > 3 * baseLineSigma) ? baseLine / averagingPoints + : amplitude / averagingPoints; + index = i + (averagingPoints / 2 + 1); + amplitude = signal[index]; + sumAvg += (std::abs(amplitude - baseLine) > 3 * baseLineSigma) ? baseLine / averagingPoints + : amplitude / averagingPoints; + smoothed[i] = sumAvg; + } + + // Points at the end, where we can calculate a moving average + for (int i = pulseDepth - averagingPoints / 2; i < pulseDepth; i++) smoothed[i] = sumAvg; + + return smoothed; +} + +std::vector GetSignalSmoothed_ExcludeOutliers(const std::vector& signal, + int averagingPoints); +std::vector GetSignalSmoothed_ExcludeOutliers(const std::vector& signal, + int averagingPoints); + +template +std::vector GetPointsOverThreshold(const std::vector& signal, TVector2& range, + const TVector2& thrPar, Int_t nPointsOver, Int_t nPointsFlat, + Double_t baseLine, Double_t baseLineSigma) { + if (range.X() < 0) range.SetX(0); + if (range.Y() <= 0) range.SetY(signal.size()); + + std::vector pointsOverThreshold; + + double pointTh = thrPar.X(); + double signalTh = thrPar.Y(); + + double threshold = pointTh * baseLineSigma; + + for (int i = range.X(); i < range.Y(); i++) { + // Filling a pulse with consecutive points that are over threshold + double data = signal[i] - baseLine; + if (data > threshold) { + int pos = i; + std::vector pulse; + pulse.push_back(data); + 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 < range.Y() && data > threshold) { + if (TMath::Abs(signal[i] - signal[i - 1]) > threshold) { + flatN = 0; + } else { + flatN++; + } + + if (flatN < nPointsFlat) { + pulse.push_back(data); + i++; + } else { + break; + } + data = signal[i] - baseLine; + } + + if (pulse.size() >= (unsigned int)nPointsOver) { + 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 * baseLineSigma) + for (int j = pos; j < i; j++) pointsOverThreshold.push_back(j); + } + } + } + + return pointsOverThreshold; +} + +std::vector GetPointsOverThreshold(const std::vector& signal, TVector2& range, + const TVector2& thrPar, Int_t nPointsOver, Int_t nPointsFlat, + Double_t baseLine, Double_t baseLineSigma); +std::vector GetPointsOverThreshold(const std::vector& signal, TVector2& range, + const TVector2& thrPar, Int_t nPointsOver, Int_t nPointsFlat, + Double_t baseLine, Double_t baseLineSigma); + +template +Int_t GetMaxBin(const std::vector& signal) { + return std::distance(signal.begin(), std::max_element(signal.begin(), signal.end())); +} + +template +Int_t GetMinBin(const std::vector& signal) { + return std::distance(signal.begin(), std::min_element(signal.begin(), signal.end())); +} + +} // namespace TRestSignalAnalysis + +#endif From 1cc427e04b4d9260a65ece4b42e7fb4ab497a352 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Wed, 15 Feb 2023 17:37:35 +0000 Subject: [PATCH 02/21] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- source/framework/analysis/inc/TRestSignalAnalysis.h | 2 ++ 1 file changed, 2 insertions(+) diff --git a/source/framework/analysis/inc/TRestSignalAnalysis.h b/source/framework/analysis/inc/TRestSignalAnalysis.h index 8d0883e26..ecf4504c8 100644 --- a/source/framework/analysis/inc/TRestSignalAnalysis.h +++ b/source/framework/analysis/inc/TRestSignalAnalysis.h @@ -22,6 +22,7 @@ #include #include + #include #include @@ -29,6 +30,7 @@ #define RestCore_TRestSignalAnalysis #include + #include /// This namespace define utilities (functions) to calculate different signal parameters From 498f226353e76dedf585163f4a32e1f660d80364 Mon Sep 17 00:00:00 2001 From: juanan Date: Wed, 15 Feb 2023 20:12:01 +0100 Subject: [PATCH 03/21] Addressing compilation warning --- source/framework/analysis/inc/TRestSignalAnalysis.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/source/framework/analysis/inc/TRestSignalAnalysis.h b/source/framework/analysis/inc/TRestSignalAnalysis.h index ecf4504c8..a34cef6ce 100644 --- a/source/framework/analysis/inc/TRestSignalAnalysis.h +++ b/source/framework/analysis/inc/TRestSignalAnalysis.h @@ -122,7 +122,7 @@ Double_t GetAverage(const std::vector& signal, Int_t startBin, Int_t en /// template std::vector GetSignalSmoothed(const std::vector& signal, int averagingPoints = 3) { - const size_t pulseDepth = signal.size(); + const int pulseDepth = signal.size(); std::vector smoothed(pulseDepth, 0); averagingPoints = (averagingPoints / 2) * 2 + 1; // make it odd >= averagingPoints @@ -151,7 +151,7 @@ std::vector GetSignalSmoothed(const std::vector& signal, int a template std::vector GetSignalSmoothed_ExcludeOutliers(const std::vector& signal, int averagingPoints, Double_t& baseLine, Double_t& baseLineSigma) { - const size_t pulseDepth = signal.size(); + const int pulseDepth = signal.size(); std::vector smoothed(pulseDepth, 0); if (baseLine == 0) CalculateBaselineAndSigmaIQR(signal, 5, int(pulseDepth - 5), baseLine, baseLineSigma); From 64e96092da3623e0146f5d6f13851f06be5d231b Mon Sep 17 00:00:00 2001 From: juanan Date: Wed, 15 Feb 2023 22:45:54 +0100 Subject: [PATCH 04/21] Moving TRestSignalAnalysis implementation to source file --- .../analysis/inc/TRestSignalAnalysis.h | 189 +----------- .../analysis/src/TRestSignalAnalysis.cxx | 286 ++++++++++++++++++ 2 files changed, 289 insertions(+), 186 deletions(-) create mode 100644 source/framework/analysis/src/TRestSignalAnalysis.cxx diff --git a/source/framework/analysis/inc/TRestSignalAnalysis.h b/source/framework/analysis/inc/TRestSignalAnalysis.h index a34cef6ce..a505ff453 100644 --- a/source/framework/analysis/inc/TRestSignalAnalysis.h +++ b/source/framework/analysis/inc/TRestSignalAnalysis.h @@ -40,75 +40,14 @@ namespace TRestSignalAnalysis { template void CalculateBaselineAndSigmaSD(const std::vector& signal, Int_t startBin, Int_t endBin, - Double_t& baseLine, Double_t& baseLineSigma) { - baseLine = 0; - baseLineSigma = 0; - - int nPoints = 0; - - for (auto i = startBin; i < endBin; i++) { - if (i < 0 || i > signal.size()) continue; - baseLine += signal[i]; - baseLineSigma += signal[i] * signal[i]; - nPoints++; - } - - if (nPoints > 0) { - baseLine /= (double)nPoints; - baseLineSigma = TMath::Sqrt(baseLineSigma / (double)nPoints - baseLine * baseLine); - } -} -void CalculateBaselineAndSigmaSD(const std::vector& signal, Int_t startBin, Int_t endBin, - Double_t& baseLine, Double_t& baseLineSigma); -void CalculateBaselineAndSigmaSD(const std::vector& signal, Int_t startBin, Int_t endBin, Double_t& baseLine, Double_t& baseLineSigma); template void CalculateBaselineAndSigmaIQR(const std::vector& signal, Int_t startBin, Int_t endBin, - Double_t& baseLine, Double_t& baseLineSigma) { - baseLine = 0; - baseLineSigma = 0; - - if (startBin < 0) startBin = 0; - if (endBin >= signal.size()) endBin = signal.size() - 1; - - if (endBin > startBin) return; - - auto first = signal.begin() + startBin; - auto last = signal.begin() + endBin; - std::vector v(first, last); - baseLine = TMath::Median(endBin - startBin, &v[0]); - - std::sort(v.begin(), v.end()); - baseLineSigma = - v[(int)(endBin - startBin) * 0.75] - - v[(int)(endBin - startBin) * 0.25] / - 1.349; // IQR/1.349 equals the standard deviation in case of normally distributed data -} - -void CalculateBaselineAndSigmaIQR(const std::vector& signal, Int_t startBin, Int_t endBin, - Double_t& baseLine, Double_t& baseLineSigma); -void CalculateBaselineAndSigmaIQR(const std::vector& signal, Int_t startBin, Int_t endBin, Double_t& baseLine, Double_t& baseLineSigma); template -Double_t GetAverage(const std::vector& signal, Int_t startBin, Int_t endBin) { - int nPoints = 0; - Double_t avg = 0; - - for (auto i = startBin; i < endBin; i++) { - if (i < 0 || i > signal.size()) continue; - avg += signal[i]; - nPoints++; - } - - if (nPoints > 0) avg /= (double)nPoints; - - return avg; -} - -Double_t GetAverage(const std::vector& signal, Int_t startBin, Int_t endBin); -Double_t GetAverage(const std::vector& signal, Int_t startBin, Int_t endBin); +Double_t GetAverage(const std::vector& signal, Int_t startBin, Int_t endBin); /////////////////////////////////////////////// /// \brief Return smoothing of signal @@ -121,136 +60,14 @@ Double_t GetAverage(const std::vector& signal, Int_t startBin, Int_t en /// baseline will be ignored to improve the smoothing result /// template -std::vector GetSignalSmoothed(const std::vector& signal, int averagingPoints = 3) { - const int pulseDepth = signal.size(); - std::vector smoothed(pulseDepth, 0); - - averagingPoints = (averagingPoints / 2) * 2 + 1; // make it odd >= averagingPoints - - Float_t sumAvg = GetAverage(signal, 0, averagingPoints); - - // Points at the beginning, where we can calculate a moving average - for (int i = 0; i <= averagingPoints / 2; i++) smoothed[i] = sumAvg; - - for (int i = averagingPoints / 2 + 1; i < pulseDepth - averagingPoints / 2; i++) { - int index = i - (averagingPoints / 2 + 1); - sumAvg -= signal[index] / averagingPoints; - index = i + (averagingPoints / 2 + 1); - sumAvg += signal[index] / averagingPoints; - smoothed[i] = sumAvg; - } - - // Points at the end, where we can calculate a moving average - for (int i = pulseDepth - averagingPoints / 2; i < pulseDepth; i++) smoothed[i] = sumAvg; - - return smoothed; -} -std::vector GetSignalSmoothed(const std::vector& signal, int averagingPoints = 3); -std::vector GetSignalSmoothed(const std::vector& signal, int averagingPoints = 3); +std::vector GetSignalSmoothed(const std::vector& signal, int averagingPoints = 3); template std::vector GetSignalSmoothed_ExcludeOutliers(const std::vector& signal, int averagingPoints, - Double_t& baseLine, Double_t& baseLineSigma) { - const int pulseDepth = signal.size(); - std::vector smoothed(pulseDepth, 0); - - if (baseLine == 0) CalculateBaselineAndSigmaIQR(signal, 5, int(pulseDepth - 5), baseLine, baseLineSigma); - - averagingPoints = (averagingPoints / 2) * 2 + 1; // make it odd >= averagingPoints - - Float_t sumAvg = GetAverage(signal, 0, averagingPoints); - - // Points at the beginning, where we can calculate a moving average - for (int i = 0; i <= averagingPoints / 2; i++) smoothed[i] = sumAvg; - - // Points in the middle - float_t amplitude; - for (int i = averagingPoints / 2 + 1; i < pulseDepth - averagingPoints / 2; i++) { - int index = i - (averagingPoints / 2 + 1); - amplitude = signal[index]; - sumAvg -= (std::abs(amplitude - baseLine) > 3 * baseLineSigma) ? baseLine / averagingPoints - : amplitude / averagingPoints; - index = i + (averagingPoints / 2 + 1); - amplitude = signal[index]; - sumAvg += (std::abs(amplitude - baseLine) > 3 * baseLineSigma) ? baseLine / averagingPoints - : amplitude / averagingPoints; - smoothed[i] = sumAvg; - } - - // Points at the end, where we can calculate a moving average - for (int i = pulseDepth - averagingPoints / 2; i < pulseDepth; i++) smoothed[i] = sumAvg; - - return smoothed; -} - -std::vector GetSignalSmoothed_ExcludeOutliers(const std::vector& signal, - int averagingPoints); -std::vector GetSignalSmoothed_ExcludeOutliers(const std::vector& signal, - int averagingPoints); + Double_t& baseLine, Double_t& baseLineSigma); template std::vector GetPointsOverThreshold(const std::vector& signal, TVector2& range, - const TVector2& thrPar, Int_t nPointsOver, Int_t nPointsFlat, - Double_t baseLine, Double_t baseLineSigma) { - if (range.X() < 0) range.SetX(0); - if (range.Y() <= 0) range.SetY(signal.size()); - - std::vector pointsOverThreshold; - - double pointTh = thrPar.X(); - double signalTh = thrPar.Y(); - - double threshold = pointTh * baseLineSigma; - - for (int i = range.X(); i < range.Y(); i++) { - // Filling a pulse with consecutive points that are over threshold - double data = signal[i] - baseLine; - if (data > threshold) { - int pos = i; - std::vector pulse; - pulse.push_back(data); - 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 < range.Y() && data > threshold) { - if (TMath::Abs(signal[i] - signal[i - 1]) > threshold) { - flatN = 0; - } else { - flatN++; - } - - if (flatN < nPointsFlat) { - pulse.push_back(data); - i++; - } else { - break; - } - data = signal[i] - baseLine; - } - - if (pulse.size() >= (unsigned int)nPointsOver) { - 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 * baseLineSigma) - for (int j = pos; j < i; j++) pointsOverThreshold.push_back(j); - } - } - } - - return pointsOverThreshold; -} - -std::vector GetPointsOverThreshold(const std::vector& signal, TVector2& range, - const TVector2& thrPar, Int_t nPointsOver, Int_t nPointsFlat, - Double_t baseLine, Double_t baseLineSigma); -std::vector GetPointsOverThreshold(const std::vector& signal, TVector2& range, const TVector2& thrPar, Int_t nPointsOver, Int_t nPointsFlat, Double_t baseLine, Double_t baseLineSigma); diff --git a/source/framework/analysis/src/TRestSignalAnalysis.cxx b/source/framework/analysis/src/TRestSignalAnalysis.cxx new file mode 100644 index 000000000..707e91fc2 --- /dev/null +++ b/source/framework/analysis/src/TRestSignalAnalysis.cxx @@ -0,0 +1,286 @@ +/************************************************************************* + * 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. * + *************************************************************************/ + +////////////////////////////////////////////////////////////////////////// +/// TRestSignalAnalysis defines several functions to calculate different +/// signal parameters. +///-------------------------------------------------------------------------- +/// +/// REST-for-Physics - Software for Rare Event Searches Toolkit +/// +/// History of developments: +/// +/// 2022-December First implementation +/// JuanAn Garcia +/// +/// \class TRestSignalAnalysis +/// \author: JuanAn Garcia e-mail: juanangp@unizar.es +/// +///
+/// + +/////////////////////////////////////////////// +/// \brief Smoothing of the existing signal and returns +/// a vector of Float_t values +/// +/// \param neighbours It defines the number of neighbours +/// points used to average the signal +/// +/// \param option If the option is set to "EXCLUDE OUTLIERS", points that are too far away from the median +/// baseline will be ignored to improve the smoothing result +/// + +#include + +template +void TRestSignalAnalysis::CalculateBaselineAndSigmaSD(const std::vector& signal, Int_t startBin, + Int_t endBin, Double_t& baseLine, + Double_t& baseLineSigma) { + baseLine = 0; + baseLineSigma = 0; + + int nPoints = 0; + + for (int i = startBin; i < endBin; i++) { + if (i < 0 || i > (int)signal.size()) continue; + baseLine += signal[i]; + baseLineSigma += signal[i] * signal[i]; + nPoints++; + } + + if (nPoints > 0) { + baseLine /= (double)nPoints; + baseLineSigma = TMath::Sqrt(baseLineSigma / (double)nPoints - baseLine * baseLine); + } +} +template void TRestSignalAnalysis::CalculateBaselineAndSigmaSD(const std::vector& signal, + Int_t startBin, Int_t endBin, + Double_t& baseLine, + Double_t& baseLineSigma); +template void TRestSignalAnalysis::CalculateBaselineAndSigmaSD(const std::vector& signal, + Int_t startBin, Int_t endBin, + Double_t& baseLine, + Double_t& baseLineSigma); + +template +void TRestSignalAnalysis::CalculateBaselineAndSigmaIQR(const std::vector& signal, Int_t startBin, + Int_t endBin, Double_t& baseLine, + Double_t& baseLineSigma) { + baseLine = 0; + baseLineSigma = 0; + + if (startBin < 0) startBin = 0; + if (endBin >= (int)signal.size()) endBin = signal.size() - 1; + + if (endBin > startBin) return; + + auto first = signal.begin() + startBin; + auto last = signal.begin() + endBin; + std::vector v(first, last); + baseLine = TMath::Median(endBin - startBin, &v[0]); + + std::sort(v.begin(), v.end()); + baseLineSigma = + v[(int)(endBin - startBin) * 0.75] - + v[(int)(endBin - startBin) * 0.25] / + 1.349; // IQR/1.349 equals the standard deviation in case of normally distributed data +} + +template void TRestSignalAnalysis::CalculateBaselineAndSigmaIQR(const std::vector& signal, + Int_t startBin, Int_t endBin, + Double_t& baseLine, + Double_t& baseLineSigma); +template void TRestSignalAnalysis::CalculateBaselineAndSigmaIQR(const std::vector& signal, + Int_t startBin, Int_t endBin, + Double_t& baseLine, + Double_t& baseLineSigma); + +template +Double_t TRestSignalAnalysis::GetAverage(const std::vector& signal, Int_t startBin, Int_t endBin) { + int nPoints = 0; + Double_t avg = 0; + + for (int i = startBin; i < endBin; i++) { + if (i < 0 || i > (int)signal.size()) continue; + avg += signal[i]; + nPoints++; + } + + if (nPoints > 0) avg /= (double)nPoints; + + return avg; +} + +template Double_t TRestSignalAnalysis::GetAverage(const std::vector& signal, Int_t startBin, + Int_t endBin); +template Double_t TRestSignalAnalysis::GetAverage(const std::vector& signal, Int_t startBin, + Int_t endBin); + +/////////////////////////////////////////////// +/// \brief Return smoothing of signal +/// +/// +/// \param neighbours It defines the number of neighbours +/// points used to average the signal +/// +/// \param option If the option is set to "EXCLUDE OUTLIERS", points that are too far away from the median +/// baseline will be ignored to improve the smoothing result +/// +template +std::vector TRestSignalAnalysis::GetSignalSmoothed(const std::vector& signal, + int averagingPoints) { + const int pulseDepth = signal.size(); + std::vector smoothed(pulseDepth, 0); + + averagingPoints = (averagingPoints / 2) * 2 + 1; // make it odd >= averagingPoints + + Float_t sumAvg = TRestSignalAnalysis::GetAverage(signal, 0, averagingPoints); + + // Points at the beginning, where we can calculate a moving average + for (int i = 0; i <= averagingPoints / 2; i++) smoothed[i] = sumAvg; + + for (int i = averagingPoints / 2 + 1; i < pulseDepth - averagingPoints / 2; i++) { + int index = i - (averagingPoints / 2 + 1); + sumAvg -= signal[index] / averagingPoints; + index = i + (averagingPoints / 2 + 1); + sumAvg += signal[index] / averagingPoints; + smoothed[i] = sumAvg; + } + + // Points at the end, where we can calculate a moving average + for (int i = pulseDepth - averagingPoints / 2; i < pulseDepth; i++) smoothed[i] = sumAvg; + + return smoothed; +} +template std::vector TRestSignalAnalysis::GetSignalSmoothed( + const std::vector& signal, int averagingPoints); +template std::vector TRestSignalAnalysis::GetSignalSmoothed( + const std::vector& signal, int averagingPoints); + +template +std::vector TRestSignalAnalysis::GetSignalSmoothed_ExcludeOutliers(const std::vector& signal, + int averagingPoints, + Double_t& baseLine, + Double_t& baseLineSigma) { + const int pulseDepth = signal.size(); + std::vector smoothed(pulseDepth, 0); + + if (baseLine == 0) CalculateBaselineAndSigmaIQR(signal, 5, int(pulseDepth - 5), baseLine, baseLineSigma); + + averagingPoints = (averagingPoints / 2) * 2 + 1; // make it odd >= averagingPoints + + Float_t sumAvg = TRestSignalAnalysis::GetAverage(signal, 0, averagingPoints); + + // Points at the beginning, where we can calculate a moving average + for (int i = 0; i <= averagingPoints / 2; i++) smoothed[i] = sumAvg; + + // Points in the middle + float_t amplitude; + for (int i = averagingPoints / 2 + 1; i < pulseDepth - averagingPoints / 2; i++) { + int index = i - (averagingPoints / 2 + 1); + amplitude = signal[index]; + sumAvg -= (std::abs(amplitude - baseLine) > 3 * baseLineSigma) ? baseLine / averagingPoints + : amplitude / averagingPoints; + index = i + (averagingPoints / 2 + 1); + amplitude = signal[index]; + sumAvg += (std::abs(amplitude - baseLine) > 3 * baseLineSigma) ? baseLine / averagingPoints + : amplitude / averagingPoints; + smoothed[i] = sumAvg; + } + + // Points at the end, where we can calculate a moving average + for (int i = pulseDepth - averagingPoints / 2; i < pulseDepth; i++) smoothed[i] = sumAvg; + + return smoothed; +} + +template std::vector TRestSignalAnalysis::GetSignalSmoothed_ExcludeOutliers( + const std::vector& signal, int averagingPoints, Double_t& baseLine, Double_t& baseLineSigma); +template std::vector TRestSignalAnalysis::GetSignalSmoothed_ExcludeOutliers( + const std::vector& signal, int averagingPoints, Double_t& baseLine, Double_t& baseLineSigma); + +template +std::vector TRestSignalAnalysis::GetPointsOverThreshold(const std::vector& signal, TVector2& range, + const TVector2& thrPar, Int_t nPointsOver, + Int_t nPointsFlat, Double_t baseLine, + Double_t baseLineSigma) { + if (range.X() < 0) range.SetX(0); + if (range.Y() <= 0) range.SetY(signal.size()); + + std::vector pointsOverThreshold; + + double pointTh = thrPar.X(); + double signalTh = thrPar.Y(); + + double threshold = pointTh * baseLineSigma; + + for (int i = range.X(); i < range.Y(); i++) { + // Filling a pulse with consecutive points that are over threshold + double data = signal[i] - baseLine; + if (data > threshold) { + int pos = i; + std::vector pulse; + pulse.push_back(data); + 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 < range.Y() && data > threshold) { + if (TMath::Abs(signal[i] - signal[i - 1]) > threshold) { + flatN = 0; + } else { + flatN++; + } + + if (flatN < nPointsFlat) { + pulse.push_back(data); + i++; + } else { + break; + } + data = signal[i] - baseLine; + } + + if (pulse.size() >= (unsigned int)nPointsOver) { + 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 * baseLineSigma) + for (int j = pos; j < i; j++) pointsOverThreshold.push_back(j); + } + } + } + + return pointsOverThreshold; +} + +template std::vector TRestSignalAnalysis::GetPointsOverThreshold( + const std::vector& signal, TVector2& range, const TVector2& thrPar, Int_t nPointsOver, + Int_t nPointsFlat, Double_t baseLine, Double_t baseLineSigma); +template std::vector TRestSignalAnalysis::GetPointsOverThreshold( + const std::vector& signal, TVector2& range, const TVector2& thrPar, Int_t nPointsOver, + Int_t nPointsFlat, Double_t baseLine, Double_t baseLineSigma); From b3bbf90d54da7a67e4fd3b161f6ea3a6b917460b Mon Sep 17 00:00:00 2001 From: juanan Date: Thu, 16 Feb 2023 10:21:37 +0100 Subject: [PATCH 05/21] Addressing pipeline issue --- source/framework/analysis/src/TRestSignalAnalysis.cxx | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/source/framework/analysis/src/TRestSignalAnalysis.cxx b/source/framework/analysis/src/TRestSignalAnalysis.cxx index 707e91fc2..b0319255c 100644 --- a/source/framework/analysis/src/TRestSignalAnalysis.cxx +++ b/source/framework/analysis/src/TRestSignalAnalysis.cxx @@ -235,11 +235,10 @@ std::vector TRestSignalAnalysis::GetPointsOverThreshold(const std::vector for (int i = range.X(); i < range.Y(); i++) { // Filling a pulse with consecutive points that are over threshold - double data = signal[i] - baseLine; - if (data > threshold) { + if ((signal[i] - baseLine) > threshold) { int pos = i; std::vector pulse; - pulse.push_back(data); + pulse.push_back(signal[i] - baseLine); i++; // If the pulse ends in a flat end above the threshold, the parameter @@ -248,7 +247,7 @@ std::vector TRestSignalAnalysis::GetPointsOverThreshold(const std::vector // decision to cut this anomalous behaviour. And all points over threshold // will be added to the pulse vector. int flatN = 0; - while (i < range.Y() && data > threshold) { + while (i < range.Y() && (signal[i] - baseLine) > threshold) { if (TMath::Abs(signal[i] - signal[i - 1]) > threshold) { flatN = 0; } else { @@ -256,12 +255,11 @@ std::vector TRestSignalAnalysis::GetPointsOverThreshold(const std::vector } if (flatN < nPointsFlat) { - pulse.push_back(data); + pulse.push_back(signal[i] - baseLine); i++; } else { break; } - data = signal[i] - baseLine; } if (pulse.size() >= (unsigned int)nPointsOver) { From d18699e763da504598a13630660799061c143ba4 Mon Sep 17 00:00:00 2001 From: juanan Date: Thu, 16 Feb 2023 18:41:32 +0100 Subject: [PATCH 06/21] Commenting the code for documentation --- .../analysis/inc/TRestSignalAnalysis.h | 14 +--- .../analysis/src/TRestSignalAnalysis.cxx | 82 ++++++++++++++----- 2 files changed, 63 insertions(+), 33 deletions(-) diff --git a/source/framework/analysis/inc/TRestSignalAnalysis.h b/source/framework/analysis/inc/TRestSignalAnalysis.h index a505ff453..0360887d3 100644 --- a/source/framework/analysis/inc/TRestSignalAnalysis.h +++ b/source/framework/analysis/inc/TRestSignalAnalysis.h @@ -33,11 +33,9 @@ #include -/// This namespace define utilities (functions) to calculate different signal parameters +/// This namespace define generic functions to calculate different signal parameters namespace TRestSignalAnalysis { -/// Generic functions for different calculations - template void CalculateBaselineAndSigmaSD(const std::vector& signal, Int_t startBin, Int_t endBin, Double_t& baseLine, Double_t& baseLineSigma); @@ -49,16 +47,6 @@ void CalculateBaselineAndSigmaIQR(const std::vector& signal, Int_t startBin, template Double_t GetAverage(const std::vector& signal, Int_t startBin, Int_t endBin); -/////////////////////////////////////////////// -/// \brief Return smoothing of signal -/// -/// -/// \param neighbours It defines the number of neighbours -/// points used to average the signal -/// -/// \param option If the option is set to "EXCLUDE OUTLIERS", points that are too far away from the median -/// baseline will be ignored to improve the smoothing result -/// template std::vector GetSignalSmoothed(const std::vector& signal, int averagingPoints = 3); diff --git a/source/framework/analysis/src/TRestSignalAnalysis.cxx b/source/framework/analysis/src/TRestSignalAnalysis.cxx index b0319255c..17cac7c7e 100644 --- a/source/framework/analysis/src/TRestSignalAnalysis.cxx +++ b/source/framework/analysis/src/TRestSignalAnalysis.cxx @@ -38,19 +38,14 @@ ///
/// -/////////////////////////////////////////////// -/// \brief Smoothing of the existing signal and returns -/// a vector of Float_t values -/// -/// \param neighbours It defines the number of neighbours -/// points used to average the signal -/// -/// \param option If the option is set to "EXCLUDE OUTLIERS", points that are too far away from the median -/// baseline will be ignored to improve the smoothing result -/// - #include +/////////////////////////////////////////////// +/// \brief This method is used to determine the value +/// of the baseline as average (arithmetic mean) of the +/// data in the range defined between startBin and endBin. +/// The baseline sigma is determined as the standard deviation +/// of the baseline in range provided. template void TRestSignalAnalysis::CalculateBaselineAndSigmaSD(const std::vector& signal, Int_t startBin, Int_t endBin, Double_t& baseLine, @@ -80,7 +75,13 @@ template void TRestSignalAnalysis::CalculateBaselineAndSigmaSD(const st Int_t startBin, Int_t endBin, Double_t& baseLine, Double_t& baseLineSigma); - +/////////////////////////////////////////////// +/// \brief This method is used to determine the value +/// of the baseline as the median of the data in +/// the range defined between startBin and endBin. +/// The baseline sigma is determined as the interquartile +/// range (IQR) in the baseline range provided. The IQR +/// is more robust towards outliers than the standard deviation. template void TRestSignalAnalysis::CalculateBaselineAndSigmaIQR(const std::vector& signal, Int_t startBin, Int_t endBin, Double_t& baseLine, @@ -113,7 +114,10 @@ template void TRestSignalAnalysis::CalculateBaselineAndSigmaIQR(const s Int_t startBin, Int_t endBin, Double_t& baseLine, Double_t& baseLineSigma); - +/////////////////////////////////////////////// +/// \brief This method performs the average of +/// the data points in a given range defined +/// between startBin and endBin template Double_t TRestSignalAnalysis::GetAverage(const std::vector& signal, Int_t startBin, Int_t endBin) { int nPoints = 0; @@ -136,14 +140,11 @@ template Double_t TRestSignalAnalysis::GetAverage(const std::vector std::vector TRestSignalAnalysis::GetSignalSmoothed(const std::vector& signal, @@ -176,6 +177,14 @@ template std::vector TRestSignalAnalysis::GetSignalSmoothed( template std::vector TRestSignalAnalysis::GetSignalSmoothed( const std::vector& signal, int averagingPoints); +/////////////////////////////////////////////// +/// \brief Return smoothing of signal, the points +/// that are too far away from the median baseline +/// will be ignored to improve the smoothing result. +/// +/// \param averagingPoints defines the number of +/// neighbouring points used to average the signal +/// template std::vector TRestSignalAnalysis::GetSignalSmoothed_ExcludeOutliers(const std::vector& signal, int averagingPoints, @@ -218,6 +227,39 @@ template std::vector TRestSignalAnalysis::GetSignalSmoothed_ExcludeOutl template std::vector TRestSignalAnalysis::GetSignalSmoothed_ExcludeOutliers( const std::vector& signal, int averagingPoints, Double_t& baseLine, Double_t& baseLineSigma); +/////////////////////////////////////////////// +/// \brief It returns a vector of the data points +/// that are found over threshold. +/// The parameters provided to this method are +/// used to identify those points. +/// +/// \param thrPar A TVector2 defining two parameters: *pointThreshold* and +/// *signalThreshold*. Both numbers +/// define the number of sigmas over the baseline fluctuation, stored in +/// baseLineSigma. The first parameter, +/// *pointThreshold*, serves to identify if a single point is over threshold by +/// satisfying the condition that +/// is above the baseline by the number of sigmas given in *pointThreshold*. +/// Once a certain number of +/// consecutive points have been identified, the parameter *signalThreshold* +/// will serve to reject the signals +/// (consecutive points over threshold) that their standard deviation is lower +/// that *signalThreshold* times +/// the baseline fluctuation. +/// +/// \param nPointsOver Only data points with at least *nPointsOver* consecutive +/// points will be considered. +/// +/// \param nPointsFlat It will serve to terminate the points over threshold +/// identification in signals where +/// we find an overshoot, being the baseline not returning to zero (or its +/// original value) at the signal tail. +/// +/// \param baseLine value of the signal baseline calculated beforehand +/// +/// \param baseLineSigma value of the baseline fluctuation calculated beforehand +/// + template std::vector TRestSignalAnalysis::GetPointsOverThreshold(const std::vector& signal, TVector2& range, const TVector2& thrPar, Int_t nPointsOver, From 3b0e94dd908dfcbe6002bea4adb164e692a147c8 Mon Sep 17 00:00:00 2001 From: juanan Date: Thu, 16 Feb 2023 18:41:55 +0100 Subject: [PATCH 07/21] Implementing changes after moving pointsOverThreshold to a pair --- .../framework/analysis/inc/TRestSignalAnalysis.h | 9 +++++---- .../analysis/src/TRestSignalAnalysis.cxx | 16 ++++++++-------- 2 files changed, 13 insertions(+), 12 deletions(-) diff --git a/source/framework/analysis/inc/TRestSignalAnalysis.h b/source/framework/analysis/inc/TRestSignalAnalysis.h index 0360887d3..c55d29003 100644 --- a/source/framework/analysis/inc/TRestSignalAnalysis.h +++ b/source/framework/analysis/inc/TRestSignalAnalysis.h @@ -30,7 +30,7 @@ #define RestCore_TRestSignalAnalysis #include - +#include #include /// This namespace define generic functions to calculate different signal parameters @@ -55,9 +55,10 @@ std::vector GetSignalSmoothed_ExcludeOutliers(const std::vector& sig Double_t& baseLine, Double_t& baseLineSigma); template -std::vector GetPointsOverThreshold(const std::vector& signal, TVector2& range, - const TVector2& thrPar, Int_t nPointsOver, Int_t nPointsFlat, - Double_t baseLine, Double_t baseLineSigma); +std::vector > GetPointsOverThreshold(const std::vector& signal, + TVector2& range, const TVector2& thrPar, + Int_t nPointsOver, Int_t nPointsFlat, + Double_t baseLine, Double_t baseLineSigma); template Int_t GetMaxBin(const std::vector& signal) { diff --git a/source/framework/analysis/src/TRestSignalAnalysis.cxx b/source/framework/analysis/src/TRestSignalAnalysis.cxx index 17cac7c7e..a8acc8940 100644 --- a/source/framework/analysis/src/TRestSignalAnalysis.cxx +++ b/source/framework/analysis/src/TRestSignalAnalysis.cxx @@ -261,14 +261,13 @@ template std::vector TRestSignalAnalysis::GetSignalSmoothed_ExcludeOutl /// template -std::vector TRestSignalAnalysis::GetPointsOverThreshold(const std::vector& signal, TVector2& range, - const TVector2& thrPar, Int_t nPointsOver, - Int_t nPointsFlat, Double_t baseLine, - Double_t baseLineSigma) { +std::vector > TRestSignalAnalysis::GetPointsOverThreshold( + const std::vector& signal, TVector2& range, const TVector2& thrPar, Int_t nPointsOver, + Int_t nPointsFlat, Double_t baseLine, Double_t baseLineSigma) { if (range.X() < 0) range.SetX(0); if (range.Y() <= 0) range.SetY(signal.size()); - std::vector pointsOverThreshold; + std::vector > pointsOverThreshold; double pointTh = thrPar.X(); double signalTh = thrPar.Y(); @@ -310,7 +309,8 @@ std::vector TRestSignalAnalysis::GetPointsOverThreshold(const std::vector double stdev = std::sqrt(sq_sum / pulse.size() - mean * mean); if (stdev > signalTh * baseLineSigma) - for (int j = pos; j < i; j++) pointsOverThreshold.push_back(j); + for (unsigned int j = 0; j < pulse.size(); j++) + pointsOverThreshold.push_back(std::make_pair(pos + j, pulse[j])); } } } @@ -318,9 +318,9 @@ std::vector TRestSignalAnalysis::GetPointsOverThreshold(const std::vector return pointsOverThreshold; } -template std::vector TRestSignalAnalysis::GetPointsOverThreshold( +template std::vector > TRestSignalAnalysis::GetPointsOverThreshold( const std::vector& signal, TVector2& range, const TVector2& thrPar, Int_t nPointsOver, Int_t nPointsFlat, Double_t baseLine, Double_t baseLineSigma); -template std::vector TRestSignalAnalysis::GetPointsOverThreshold( +template std::vector > TRestSignalAnalysis::GetPointsOverThreshold( const std::vector& signal, TVector2& range, const TVector2& thrPar, Int_t nPointsOver, Int_t nPointsFlat, Double_t baseLine, Double_t baseLineSigma); From c219a17f62d48bb8cfc12842c385a433ff0b62b3 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 16 Feb 2023 17:42:15 +0000 Subject: [PATCH 08/21] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- source/framework/analysis/inc/TRestSignalAnalysis.h | 1 + 1 file changed, 1 insertion(+) diff --git a/source/framework/analysis/inc/TRestSignalAnalysis.h b/source/framework/analysis/inc/TRestSignalAnalysis.h index c55d29003..a33bb3aaf 100644 --- a/source/framework/analysis/inc/TRestSignalAnalysis.h +++ b/source/framework/analysis/inc/TRestSignalAnalysis.h @@ -30,6 +30,7 @@ #define RestCore_TRestSignalAnalysis #include + #include #include From 847856d4915004d4d235e9b79cad5eabf28b8f18 Mon Sep 17 00:00:00 2001 From: juanan Date: Fri, 31 Mar 2023 10:20:58 +0200 Subject: [PATCH 09/21] Adding further generic methods to TRestSignalAnalysis --- .../analysis/inc/TRestSignalAnalysis.h | 31 +- .../analysis/src/TRestSignalAnalysis.cxx | 355 ++++++++++++++++++ 2 files changed, 381 insertions(+), 5 deletions(-) diff --git a/source/framework/analysis/inc/TRestSignalAnalysis.h b/source/framework/analysis/inc/TRestSignalAnalysis.h index a33bb3aaf..52b70abbb 100644 --- a/source/framework/analysis/inc/TRestSignalAnalysis.h +++ b/source/framework/analysis/inc/TRestSignalAnalysis.h @@ -20,18 +20,17 @@ * For the list of contributors see $REST_PATH/CREDITS. * *************************************************************************/ +#ifndef RestCore_TRestSignalAnalysis +#define RestCore_TRestSignalAnalysis + #include #include #include #include -#ifndef RestCore_TRestSignalAnalysis -#define RestCore_TRestSignalAnalysis - #include - -#include +#include #include /// This namespace define generic functions to calculate different signal parameters @@ -55,6 +54,10 @@ template std::vector GetSignalSmoothed_ExcludeOutliers(const std::vector& signal, int averagingPoints, Double_t& baseLine, Double_t& baseLineSigma); +template +std::vector GetDerivative(const std::vector& signal); + + template std::vector > GetPointsOverThreshold(const std::vector& signal, TVector2& range, const TVector2& thrPar, @@ -71,6 +74,24 @@ Int_t GetMinBin(const std::vector& signal) { return std::distance(signal.begin(), std::min_element(signal.begin(), signal.end())); } +template +Double_t GetIntegral(const std::vector& signal, Int_t startBin, Int_t endBin); + +template +Double_t GetMaxPeakWidth(const std::vector& signal); + +Double_t GetSlopeIntegral(const std::vector > & signal ); +Double_t GetRiseSlope(const std::vector > & signal ); +Double_t GetRiseTime(const std::vector > & signal ); + +std::vector > GetIntWindow(TGraph *signal, double intWindow); +std::array , 3> GetTripleMax(TGraph *signal); +TVector2 GetTripleMaxAverage(TGraph *signal); +Double_t GetTripleMaxIntegral(TGraph *signal); +TVector2 GetMaxGauss(TGraph *signal); +TVector2 GetMaxLandau(TGraph *signal); +TVector2 GetMaxAget(TGraph *signal); + } // namespace TRestSignalAnalysis #endif diff --git a/source/framework/analysis/src/TRestSignalAnalysis.cxx b/source/framework/analysis/src/TRestSignalAnalysis.cxx index a8acc8940..3dc2e1ebd 100644 --- a/source/framework/analysis/src/TRestSignalAnalysis.cxx +++ b/source/framework/analysis/src/TRestSignalAnalysis.cxx @@ -40,6 +40,9 @@ #include +#include +#include + /////////////////////////////////////////////// /// \brief This method is used to determine the value /// of the baseline as average (arithmetic mean) of the @@ -227,6 +230,26 @@ template std::vector TRestSignalAnalysis::GetSignalSmoothed_ExcludeOutl template std::vector TRestSignalAnalysis::GetSignalSmoothed_ExcludeOutliers( const std::vector& signal, int averagingPoints, Double_t& baseLine, Double_t& baseLineSigma); + +/////////////////////////////////////////////// +/// \brief Return derivative of a vector of data +/// points +/// +template +std::vector TRestSignalAnalysis::GetDerivative(const std::vector& signal){ + + std::vector derivative (0, signal.size()-1); + + for(size_t i = 0;i< signal.size()-1; i++){ + derivative[i] = signal[i+1] - signal[i]; + } + + return derivative; +} + +template std::vector TRestSignalAnalysis::GetDerivative(const std::vector& signal); +template std::vector TRestSignalAnalysis::GetDerivative(const std::vector& signal); + /////////////////////////////////////////////// /// \brief It returns a vector of the data points /// that are found over threshold. @@ -324,3 +347,335 @@ template std::vector > TRestSignalAnalysis::GetPoint template std::vector > TRestSignalAnalysis::GetPointsOverThreshold( const std::vector& signal, TVector2& range, const TVector2& thrPar, Int_t nPointsOver, Int_t nPointsFlat, Double_t baseLine, Double_t baseLineSigma); + +/////////////////////////////////////////////// +/// \brief It returns the integral of the signal in the +/// range passed as argument +/// +template +Double_t TRestSignalAnalysis::GetIntegral(const std::vector& signal, Int_t startBin, Int_t endBin) { + Double_t sum = 0; + for (int i = startBin; i < endBin; i++) + if (i > 0 && i < (int)signal.size()) sum += signal[i]; + + return sum; +} +template Double_t TRestSignalAnalysis::GetIntegral(const std::vector& signal, Int_t startBin, Int_t endBin); +template Double_t TRestSignalAnalysis::GetIntegral(const std::vector& signal, Int_t startBin, Int_t endBin); + +/////////////////////////////////////////////// +/// \brief It returns the width of the pulses +/// as the number of bins between the half of the +/// maximum +/// +template +Double_t TRestSignalAnalysis::GetMaxPeakWidth(const std::vector& signal) { + + Int_t maxIndex = GetMaxBin(signal); + Double_t maxValue = signal[maxIndex]; + + const int signalSize = signal.size(); + + Double_t value = maxValue; + Int_t rightIndex = maxIndex; + while (value > maxValue / 2. && rightIndex < signalSize) { + value = signal[rightIndex]; + rightIndex++; + } + Int_t leftIndex = maxIndex; + value = maxValue; + while (value > maxValue / 2. && leftIndex > 0) { + value = signal[rightIndex]; + leftIndex--; + } + + return rightIndex - leftIndex; +} +template Double_t TRestSignalAnalysis::GetMaxPeakWidth(const std::vector& signal); +template Double_t TRestSignalAnalysis::GetMaxPeakWidth(const std::vector& signal); + +/////////////////////////////////////////////// +/// \brief It performs a gaussian fit to the signal +/// which is stored inside a TGraph. It returns +/// a TVector 2 with the maximum and the mean of the +/// gaussian fit +/// +TVector2 TRestSignalAnalysis::GetMaxGauss(TGraph *signal){ + + Int_t maxBin = TMath::LocMax(signal->GetN(), signal->GetY()); + Double_t maxTime = signal->GetPointX(maxBin); + Double_t gaussMax = -1, gaussMean = -1; + Double_t lowerLimit = maxTime - maxTime * 0.1; // us + Double_t upperLimit = maxTime + maxTime * 0.2; // us + + TF1* gaus = new TF1("gaus", "gaus", lowerLimit, upperLimit); + TFitResultPtr fitResult = signal->Fit(gaus, "QNRS"); + + if (fitResult->IsValid()) { + gaussMax = gaus->GetParameter(0); + gaussMean = gaus->GetParameter(1); + } else { + // the fit failed, return -1 to indicate failure + std::cout << std::endl + << "WARNING: bad fit with maximum at time = " << maxTime + << "Failed fit parameters = " << gaus->GetParameter(0) << " || " << gaus->GetParameter(1) + << " || " << gaus->GetParameter(2) << "\n" + << "Assigned fit parameters : energy = " << gaussMax << ", time = " << gaussMean << std::endl; + /* + TCanvas* c2 = new TCanvas("c2", "Signal fit", 200, 10, 1280, 720); + signal->Draw(); + c2->Update(); + getchar(); + delete c2; + */ + } + + delete gaus; + + return TVector2(gaussMean, gaussMax); + +} + +/////////////////////////////////////////////// +/// \brief It performs a landau fit to the signal +/// which is stored inside a TGraph. It returns +/// a TVector 2 with the maximum and the mean of the +/// landau fit +/// +TVector2 TRestSignalAnalysis::GetMaxLandau(TGraph *signal){ + + Int_t maxBin = TMath::LocMax(signal->GetN(), signal->GetY()); + Double_t maxTime = signal->GetPointX(maxBin); + Double_t landauMax = -1, landauMean = -1; + Double_t lowerLimit = maxTime - maxTime * 0.1; // us + Double_t upperLimit = maxTime + maxTime * 0.2; // us + + TF1* landau = new TF1("landau", "landau", lowerLimit, upperLimit); + TFitResultPtr fitResult = signal->Fit(landau, "QNRS"); + + if (fitResult->IsValid()) { + landauMax = landau->GetParameter(0); + landauMean = landau->GetParameter(1); + } else { + // the fit failed, return -1 to indicate failure + std::cout << std::endl + << "WARNING: bad fit with maximum at time = " << maxTime + << "Failed fit parameters = " << landau->GetParameter(0) << " || " << landau->GetParameter(1) + << " || " << landau->GetParameter(2) << "\n" + << "Assigned fit parameters : energy = " << landauMax << ", time = " << landauMean << std::endl; + /* + TCanvas* c2 = new TCanvas("c2", "Signal fit", 200, 10, 1280, 720); + signal->Draw(); + c2->Update(); + getchar(); + delete c2; + */ + } + + delete landau; + + return TVector2(landauMean, landauMax); + +} + +/////////////////////////////////////////////// +/// \brief It performs a fit of the signal to a +/// the aget response function. It returns a +/// TVector 2 with the maximum and the mean of the +/// gaussian fit +/// +TVector2 TRestSignalAnalysis::GetMaxAget(TGraph *signal){ + + Int_t maxBin = TMath::LocMax(signal->GetN(), signal->GetY()); + Double_t maxTime = signal->GetPointX(maxBin); + Double_t agetMax = -1, agetMean = -1; + Double_t lowerLimit = maxTime - maxTime * 0.25; // us + Double_t upperLimit = maxTime + maxTime * 0.35; // us + + // 1.1664 is the x value where the maximum of the base function (i.e. without parameters) + TF1* aget = new TF1("aget", "[&](double *x, double *p){ double arg = (x[0] - par[1] + 1.1664) / par[2]; return par[0] / 0.0440895 * exp(-3 * (arg)) * (arg) * (arg) * (arg)*sin(arg);}", lowerLimit, upperLimit, 3); + TFitResultPtr fitResult = signal->Fit(aget, "QNRS"); + + if (fitResult->IsValid()) { + agetMax = aget->GetParameter(0); + agetMean = aget->GetParameter(1); + } else { + // the fit failed, return -1 to indicate failure + std::cout << std::endl + << "WARNING: bad fit with maximum at time = " << maxTime + << "Failed fit parameters = " << aget->GetParameter(0) << " || " << aget->GetParameter(1) + << " || " << aget->GetParameter(2) << "\n" + << "Assigned fit parameters : energy = " << agetMax << ", time = " << agetMean << std::endl; + /* + TCanvas* c2 = new TCanvas("c2", "Signal fit", 200, 10, 1280, 720); + signal->Draw(); + c2->Update(); + getchar(); + delete c2; + */ + } + + delete aget; + + return TVector2(agetMean, agetMax); + +} + +/////////////////////////////////////////////// +/// \brief It averages the pulse inside the time window +/// passed as argument. It returns a vector of pairs +/// with the integrated window time and energy (charge) +/// +std::vector > TRestSignalAnalysis::GetIntWindow(TGraph *signal, double intWindow){ + + const int nPoints = signal->GetN(); + + std::map > windowMap; + for (int j = 0; j < nPoints; j++) { + int index = signal->GetPointX(j) / intWindow; + auto it = windowMap.find(index); + if (it != windowMap.end()) { + it->second.first++; + it->second.second += signal->GetPointY(j); + } else { + windowMap[index] = std::make_pair(1, signal->GetPointY(j)); + } + } + + std::vector > result; + + for (const auto& [index, pair] : windowMap) { + Double_t hitTime = index * intWindow + intWindow / 2.; + Double_t energy = pair.second / pair.first; + result.push_back(std::make_pair( hitTime, energy ) ); + } + + return result; + +} + +/////////////////////////////////////////////// +/// \brief It calculates the triple max over a TGraph +/// and returns an array of pairs with the time and +/// the energy of the maximum and the neighbouring +/// points(three points in total). +/// +std::array , 3> TRestSignalAnalysis::GetTripleMax(TGraph *signal){ + + Int_t maxBin = TMath::LocMax(signal->GetN(), signal->GetY()); + + std::array , 3> tripleMax; + + for(int i = 0; i < 3; i++){ + int index = maxBin + i - 1; + if(index >0 && index < signal->GetN() ){ + tripleMax[i] = std::make_pair(signal->GetPointX(index), signal->GetPointY(index)); + } else { + tripleMax[i] = std::make_pair(signal->GetPointX(maxBin), signal->GetPointY(maxBin)); + } + } + + return tripleMax; + +} + +/////////////////////////////////////////////// +/// \brief It calculates the triple max average over +/// a TGraph and returns a TVector2 with the average +/// time and energy (charge) +/// +TVector2 TRestSignalAnalysis::GetTripleMaxAverage(TGraph *signal){ + + auto tripleMax = TRestSignalAnalysis::GetTripleMax(signal); + double eAvg = 0; + double hitTimeAvg = 0; + for(const auto & [hitTime, energy] : tripleMax){ + hitTimeAvg += hitTime*energy; + eAvg += energy; + } + + hitTimeAvg /= eAvg; + eAvg /= 3.; + + return TVector2(hitTimeAvg, eAvg); + +} + +/////////////////////////////////////////////// +/// \brief It calculates the triple max integral over +/// a TGraph and returns the addition of the maximum plus +/// the neigbouring bins. +/// +Double_t TRestSignalAnalysis::GetTripleMaxIntegral(TGraph *signal){ + + auto tripleMax = TRestSignalAnalysis::GetTripleMax(signal); + double totEnergy = 0; + for(const auto & [hitTime, energy] : tripleMax){ + totEnergy += energy; + } + + return totEnergy; + +} + +/////////////////////////////////////////////// +/// \brief It returns the integral of the first positive +/// rise (risetime) over a vector of pairs, that should +/// correspond to the points over threshold for a given signal. +/// +Double_t TRestSignalAnalysis::GetSlopeIntegral(const std::vector > & signal ) { + Double_t sum = 0; + /*Double_t pVal = 0; + for (const auto& [index, val] : signal) { + if (val - pVal < 0) break; + sum += val; + pVal = val; + }*/ + auto max = std::max_element(std::begin(signal), std::end(signal), + [](const auto& p1, const auto& p2) { return p1.second < p2.second; }); + + for(auto it = signal.begin(); it != max; ++it) + sum += it->second; + + return sum; +} + +/////////////////////////////////////////////// +/// \brief It returns the slope of the first positive +/// rise (risetime) over a vector of pairs, that should +/// correspond to the points over threshold for a given signal. +/// +Double_t TRestSignalAnalysis::GetRiseSlope(const std::vector > & signal ) { + + if (signal.size() < 2) + return 0; + + auto max = std::max_element(std::begin(signal), std::end(signal), + [](const auto& p1, const auto& p2) { return p1.second < p2.second; }); + + auto maxBin = max->first; + auto startBin = signal.front().first; + Double_t hP = max->second; + Double_t lP = signal.front().second; + + return (hP - lP) / (maxBin - startBin - 1); +} + +/////////////////////////////////////////////// +/// \brief It returns the time of the first positive +/// rise or risetime over a vector of pairs, that should +/// correspond to the points over threshold for a given signal. +/// +Double_t TRestSignalAnalysis::GetRiseTime(const std::vector > & signal ) { + + if (signal.size() < 2) { + return 0; + } + + auto max = std::max_element(std::begin(signal), std::end(signal), + [](const auto& p1, const auto& p2) { return p1.second < p2.second; }); + + auto maxBin = max->first; + auto startBin = signal.front().first; + return maxBin - startBin; +} From bb7ef4f0d89ae4558a201b033bb2b58026ed4f11 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 31 Mar 2023 08:21:25 +0000 Subject: [PATCH 10/21] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- .../analysis/inc/TRestSignalAnalysis.h | 30 +-- .../analysis/src/TRestSignalAnalysis.cxx | 246 +++++++++--------- 2 files changed, 130 insertions(+), 146 deletions(-) diff --git a/source/framework/analysis/inc/TRestSignalAnalysis.h b/source/framework/analysis/inc/TRestSignalAnalysis.h index 52b70abbb..3d94ea651 100644 --- a/source/framework/analysis/inc/TRestSignalAnalysis.h +++ b/source/framework/analysis/inc/TRestSignalAnalysis.h @@ -24,14 +24,13 @@ #define RestCore_TRestSignalAnalysis #include +#include #include +#include #include -#include - -#include -#include #include +#include /// This namespace define generic functions to calculate different signal parameters namespace TRestSignalAnalysis { @@ -57,7 +56,6 @@ std::vector GetSignalSmoothed_ExcludeOutliers(const std::vector& sig template std::vector GetDerivative(const std::vector& signal); - template std::vector > GetPointsOverThreshold(const std::vector& signal, TVector2& range, const TVector2& thrPar, @@ -80,17 +78,17 @@ Double_t GetIntegral(const std::vector& signal, Int_t startBin, Int_t endBin) template Double_t GetMaxPeakWidth(const std::vector& signal); -Double_t GetSlopeIntegral(const std::vector > & signal ); -Double_t GetRiseSlope(const std::vector > & signal ); -Double_t GetRiseTime(const std::vector > & signal ); - -std::vector > GetIntWindow(TGraph *signal, double intWindow); -std::array , 3> GetTripleMax(TGraph *signal); -TVector2 GetTripleMaxAverage(TGraph *signal); -Double_t GetTripleMaxIntegral(TGraph *signal); -TVector2 GetMaxGauss(TGraph *signal); -TVector2 GetMaxLandau(TGraph *signal); -TVector2 GetMaxAget(TGraph *signal); +Double_t GetSlopeIntegral(const std::vector >& signal); +Double_t GetRiseSlope(const std::vector >& signal); +Double_t GetRiseTime(const std::vector >& signal); + +std::vector > GetIntWindow(TGraph* signal, double intWindow); +std::array, 3> GetTripleMax(TGraph* signal); +TVector2 GetTripleMaxAverage(TGraph* signal); +Double_t GetTripleMaxIntegral(TGraph* signal); +TVector2 GetMaxGauss(TGraph* signal); +TVector2 GetMaxLandau(TGraph* signal); +TVector2 GetMaxAget(TGraph* signal); } // namespace TRestSignalAnalysis diff --git a/source/framework/analysis/src/TRestSignalAnalysis.cxx b/source/framework/analysis/src/TRestSignalAnalysis.cxx index 3dc2e1ebd..30fbeb41d 100644 --- a/source/framework/analysis/src/TRestSignalAnalysis.cxx +++ b/source/framework/analysis/src/TRestSignalAnalysis.cxx @@ -38,10 +38,9 @@ ///
/// -#include - #include #include +#include /////////////////////////////////////////////// /// \brief This method is used to determine the value @@ -230,21 +229,19 @@ template std::vector TRestSignalAnalysis::GetSignalSmoothed_ExcludeOutl template std::vector TRestSignalAnalysis::GetSignalSmoothed_ExcludeOutliers( const std::vector& signal, int averagingPoints, Double_t& baseLine, Double_t& baseLineSigma); - /////////////////////////////////////////////// /// \brief Return derivative of a vector of data /// points /// template -std::vector TRestSignalAnalysis::GetDerivative(const std::vector& signal){ +std::vector TRestSignalAnalysis::GetDerivative(const std::vector& signal) { + std::vector derivative(0, signal.size() - 1); - std::vector derivative (0, signal.size()-1); - - for(size_t i = 0;i< signal.size()-1; i++){ - derivative[i] = signal[i+1] - signal[i]; + for (size_t i = 0; i < signal.size() - 1; i++) { + derivative[i] = signal[i + 1] - signal[i]; } - return derivative; + return derivative; } template std::vector TRestSignalAnalysis::GetDerivative(const std::vector& signal); @@ -360,17 +357,18 @@ Double_t TRestSignalAnalysis::GetIntegral(const std::vector& signal, Int_t st return sum; } -template Double_t TRestSignalAnalysis::GetIntegral(const std::vector& signal, Int_t startBin, Int_t endBin); -template Double_t TRestSignalAnalysis::GetIntegral(const std::vector& signal, Int_t startBin, Int_t endBin); +template Double_t TRestSignalAnalysis::GetIntegral(const std::vector& signal, Int_t startBin, + Int_t endBin); +template Double_t TRestSignalAnalysis::GetIntegral(const std::vector& signal, Int_t startBin, + Int_t endBin); /////////////////////////////////////////////// /// \brief It returns the width of the pulses /// as the number of bins between the half of the /// maximum -/// +/// template Double_t TRestSignalAnalysis::GetMaxPeakWidth(const std::vector& signal) { - Int_t maxIndex = GetMaxBin(signal); Double_t maxValue = signal[maxIndex]; @@ -400,27 +398,27 @@ template Double_t TRestSignalAnalysis::GetMaxPeakWidth(const std::vectorGetN(), signal->GetY()); + Double_t maxTime = signal->GetPointX(maxBin); + Double_t gaussMax = -1, gaussMean = -1; + Double_t lowerLimit = maxTime - maxTime * 0.1; // us + Double_t upperLimit = maxTime + maxTime * 0.2; // us - Int_t maxBin = TMath::LocMax(signal->GetN(), signal->GetY()); - Double_t maxTime = signal->GetPointX(maxBin); - Double_t gaussMax = -1, gaussMean = -1; - Double_t lowerLimit = maxTime - maxTime * 0.1; // us - Double_t upperLimit = maxTime + maxTime * 0.2; // us + TF1* gaus = new TF1("gaus", "gaus", lowerLimit, upperLimit); + TFitResultPtr fitResult = signal->Fit(gaus, "QNRS"); - TF1* gaus = new TF1("gaus", "gaus", lowerLimit, upperLimit); - TFitResultPtr fitResult = signal->Fit(gaus, "QNRS"); - - if (fitResult->IsValid()) { + if (fitResult->IsValid()) { gaussMax = gaus->GetParameter(0); gaussMean = gaus->GetParameter(1); } else { // the fit failed, return -1 to indicate failure std::cout << std::endl - << "WARNING: bad fit with maximum at time = " << maxTime - << "Failed fit parameters = " << gaus->GetParameter(0) << " || " << gaus->GetParameter(1) - << " || " << gaus->GetParameter(2) << "\n" - << "Assigned fit parameters : energy = " << gaussMax << ", time = " << gaussMean << std::endl; + << "WARNING: bad fit with maximum at time = " << maxTime + << "Failed fit parameters = " << gaus->GetParameter(0) << " || " << gaus->GetParameter(1) + << " || " << gaus->GetParameter(2) << "\n" + << "Assigned fit parameters : energy = " << gaussMax << ", time = " << gaussMean + << std::endl; /* TCanvas* c2 = new TCanvas("c2", "Signal fit", 200, 10, 1280, 720); signal->Draw(); @@ -433,7 +431,6 @@ TVector2 TRestSignalAnalysis::GetMaxGauss(TGraph *signal){ delete gaus; return TVector2(gaussMean, gaussMax); - } /////////////////////////////////////////////// @@ -442,27 +439,27 @@ TVector2 TRestSignalAnalysis::GetMaxGauss(TGraph *signal){ /// a TVector 2 with the maximum and the mean of the /// landau fit /// -TVector2 TRestSignalAnalysis::GetMaxLandau(TGraph *signal){ +TVector2 TRestSignalAnalysis::GetMaxLandau(TGraph* signal) { + Int_t maxBin = TMath::LocMax(signal->GetN(), signal->GetY()); + Double_t maxTime = signal->GetPointX(maxBin); + Double_t landauMax = -1, landauMean = -1; + Double_t lowerLimit = maxTime - maxTime * 0.1; // us + Double_t upperLimit = maxTime + maxTime * 0.2; // us - Int_t maxBin = TMath::LocMax(signal->GetN(), signal->GetY()); - Double_t maxTime = signal->GetPointX(maxBin); - Double_t landauMax = -1, landauMean = -1; - Double_t lowerLimit = maxTime - maxTime * 0.1; // us - Double_t upperLimit = maxTime + maxTime * 0.2; // us + TF1* landau = new TF1("landau", "landau", lowerLimit, upperLimit); + TFitResultPtr fitResult = signal->Fit(landau, "QNRS"); - TF1* landau = new TF1("landau", "landau", lowerLimit, upperLimit); - TFitResultPtr fitResult = signal->Fit(landau, "QNRS"); - - if (fitResult->IsValid()) { + if (fitResult->IsValid()) { landauMax = landau->GetParameter(0); landauMean = landau->GetParameter(1); } else { // the fit failed, return -1 to indicate failure std::cout << std::endl - << "WARNING: bad fit with maximum at time = " << maxTime - << "Failed fit parameters = " << landau->GetParameter(0) << " || " << landau->GetParameter(1) - << " || " << landau->GetParameter(2) << "\n" - << "Assigned fit parameters : energy = " << landauMax << ", time = " << landauMean << std::endl; + << "WARNING: bad fit with maximum at time = " << maxTime + << "Failed fit parameters = " << landau->GetParameter(0) << " || " + << landau->GetParameter(1) << " || " << landau->GetParameter(2) << "\n" + << "Assigned fit parameters : energy = " << landauMax << ", time = " << landauMean + << std::endl; /* TCanvas* c2 = new TCanvas("c2", "Signal fit", 200, 10, 1280, 720); signal->Draw(); @@ -475,7 +472,6 @@ TVector2 TRestSignalAnalysis::GetMaxLandau(TGraph *signal){ delete landau; return TVector2(landauMean, landauMax); - } /////////////////////////////////////////////// @@ -484,28 +480,31 @@ TVector2 TRestSignalAnalysis::GetMaxLandau(TGraph *signal){ /// TVector 2 with the maximum and the mean of the /// gaussian fit /// -TVector2 TRestSignalAnalysis::GetMaxAget(TGraph *signal){ - - Int_t maxBin = TMath::LocMax(signal->GetN(), signal->GetY()); - Double_t maxTime = signal->GetPointX(maxBin); - Double_t agetMax = -1, agetMean = -1; - Double_t lowerLimit = maxTime - maxTime * 0.25; // us - Double_t upperLimit = maxTime + maxTime * 0.35; // us - - // 1.1664 is the x value where the maximum of the base function (i.e. without parameters) - TF1* aget = new TF1("aget", "[&](double *x, double *p){ double arg = (x[0] - par[1] + 1.1664) / par[2]; return par[0] / 0.0440895 * exp(-3 * (arg)) * (arg) * (arg) * (arg)*sin(arg);}", lowerLimit, upperLimit, 3); - TFitResultPtr fitResult = signal->Fit(aget, "QNRS"); - - if (fitResult->IsValid()) { +TVector2 TRestSignalAnalysis::GetMaxAget(TGraph* signal) { + Int_t maxBin = TMath::LocMax(signal->GetN(), signal->GetY()); + Double_t maxTime = signal->GetPointX(maxBin); + Double_t agetMax = -1, agetMean = -1; + Double_t lowerLimit = maxTime - maxTime * 0.25; // us + Double_t upperLimit = maxTime + maxTime * 0.35; // us + + // 1.1664 is the x value where the maximum of the base function (i.e. without parameters) + TF1* aget = + new TF1("aget", + "[&](double *x, double *p){ double arg = (x[0] - par[1] + 1.1664) / par[2]; return par[0] / " + "0.0440895 * exp(-3 * (arg)) * (arg) * (arg) * (arg)*sin(arg);}", + lowerLimit, upperLimit, 3); + TFitResultPtr fitResult = signal->Fit(aget, "QNRS"); + + if (fitResult->IsValid()) { agetMax = aget->GetParameter(0); agetMean = aget->GetParameter(1); } else { // the fit failed, return -1 to indicate failure std::cout << std::endl - << "WARNING: bad fit with maximum at time = " << maxTime - << "Failed fit parameters = " << aget->GetParameter(0) << " || " << aget->GetParameter(1) - << " || " << aget->GetParameter(2) << "\n" - << "Assigned fit parameters : energy = " << agetMax << ", time = " << agetMean << std::endl; + << "WARNING: bad fit with maximum at time = " << maxTime + << "Failed fit parameters = " << aget->GetParameter(0) << " || " << aget->GetParameter(1) + << " || " << aget->GetParameter(2) << "\n" + << "Assigned fit parameters : energy = " << agetMax << ", time = " << agetMean << std::endl; /* TCanvas* c2 = new TCanvas("c2", "Signal fit", 200, 10, 1280, 720); signal->Draw(); @@ -518,7 +517,6 @@ TVector2 TRestSignalAnalysis::GetMaxAget(TGraph *signal){ delete aget; return TVector2(agetMean, agetMax); - } /////////////////////////////////////////////// @@ -526,32 +524,30 @@ TVector2 TRestSignalAnalysis::GetMaxAget(TGraph *signal){ /// passed as argument. It returns a vector of pairs /// with the integrated window time and energy (charge) /// -std::vector > TRestSignalAnalysis::GetIntWindow(TGraph *signal, double intWindow){ - - const int nPoints = signal->GetN(); - - std::map > windowMap; - for (int j = 0; j < nPoints; j++) { - int index = signal->GetPointX(j) / intWindow; - auto it = windowMap.find(index); - if (it != windowMap.end()) { - it->second.first++; - it->second.second += signal->GetPointY(j); - } else { - windowMap[index] = std::make_pair(1, signal->GetPointY(j)); - } - } - - std::vector > result; +std::vector > TRestSignalAnalysis::GetIntWindow(TGraph* signal, double intWindow) { + const int nPoints = signal->GetN(); + + std::map > windowMap; + for (int j = 0; j < nPoints; j++) { + int index = signal->GetPointX(j) / intWindow; + auto it = windowMap.find(index); + if (it != windowMap.end()) { + it->second.first++; + it->second.second += signal->GetPointY(j); + } else { + windowMap[index] = std::make_pair(1, signal->GetPointY(j)); + } + } - for (const auto& [index, pair] : windowMap) { - Double_t hitTime = index * intWindow + intWindow / 2.; - Double_t energy = pair.second / pair.first; - result.push_back(std::make_pair( hitTime, energy ) ); - } + std::vector > result; - return result; + for (const auto& [index, pair] : windowMap) { + Double_t hitTime = index * intWindow + intWindow / 2.; + Double_t energy = pair.second / pair.first; + result.push_back(std::make_pair(hitTime, energy)); + } + return result; } /////////////////////////////////////////////// @@ -560,23 +556,21 @@ std::vector > TRestSignalAnalysis::GetIntWindow(TGrap /// the energy of the maximum and the neighbouring /// points(three points in total). /// -std::array , 3> TRestSignalAnalysis::GetTripleMax(TGraph *signal){ - - Int_t maxBin = TMath::LocMax(signal->GetN(), signal->GetY()); +std::array, 3> TRestSignalAnalysis::GetTripleMax(TGraph* signal) { + Int_t maxBin = TMath::LocMax(signal->GetN(), signal->GetY()); - std::array , 3> tripleMax; + std::array, 3> tripleMax; - for(int i = 0; i < 3; i++){ + for (int i = 0; i < 3; i++) { int index = maxBin + i - 1; - if(index >0 && index < signal->GetN() ){ - tripleMax[i] = std::make_pair(signal->GetPointX(index), signal->GetPointY(index)); + if (index > 0 && index < signal->GetN()) { + tripleMax[i] = std::make_pair(signal->GetPointX(index), signal->GetPointY(index)); } else { - tripleMax[i] = std::make_pair(signal->GetPointX(maxBin), signal->GetPointY(maxBin)); + tripleMax[i] = std::make_pair(signal->GetPointX(maxBin), signal->GetPointY(maxBin)); } } - return tripleMax; - + return tripleMax; } /////////////////////////////////////////////// @@ -584,21 +578,19 @@ std::array , 3> TRestSignalAnalysis::GetTripleMax(T /// a TGraph and returns a TVector2 with the average /// time and energy (charge) /// -TVector2 TRestSignalAnalysis::GetTripleMaxAverage(TGraph *signal){ - - auto tripleMax = TRestSignalAnalysis::GetTripleMax(signal); - double eAvg = 0; - double hitTimeAvg = 0; - for(const auto & [hitTime, energy] : tripleMax){ - hitTimeAvg += hitTime*energy; - eAvg += energy; - } - - hitTimeAvg /= eAvg; - eAvg /= 3.; +TVector2 TRestSignalAnalysis::GetTripleMaxAverage(TGraph* signal) { + auto tripleMax = TRestSignalAnalysis::GetTripleMax(signal); + double eAvg = 0; + double hitTimeAvg = 0; + for (const auto& [hitTime, energy] : tripleMax) { + hitTimeAvg += hitTime * energy; + eAvg += energy; + } - return TVector2(hitTimeAvg, eAvg); + hitTimeAvg /= eAvg; + eAvg /= 3.; + return TVector2(hitTimeAvg, eAvg); } /////////////////////////////////////////////// @@ -606,24 +598,22 @@ TVector2 TRestSignalAnalysis::GetTripleMaxAverage(TGraph *signal){ /// a TGraph and returns the addition of the maximum plus /// the neigbouring bins. /// -Double_t TRestSignalAnalysis::GetTripleMaxIntegral(TGraph *signal){ - - auto tripleMax = TRestSignalAnalysis::GetTripleMax(signal); - double totEnergy = 0; - for(const auto & [hitTime, energy] : tripleMax){ - totEnergy += energy; - } - - return totEnergy; +Double_t TRestSignalAnalysis::GetTripleMaxIntegral(TGraph* signal) { + auto tripleMax = TRestSignalAnalysis::GetTripleMax(signal); + double totEnergy = 0; + for (const auto& [hitTime, energy] : tripleMax) { + totEnergy += energy; + } + return totEnergy; } /////////////////////////////////////////////// /// \brief It returns the integral of the first positive -/// rise (risetime) over a vector of pairs, that should +/// rise (risetime) over a vector of pairs, that should /// correspond to the points over threshold for a given signal. /// -Double_t TRestSignalAnalysis::GetSlopeIntegral(const std::vector > & signal ) { +Double_t TRestSignalAnalysis::GetSlopeIntegral(const std::vector >& signal) { Double_t sum = 0; /*Double_t pVal = 0; for (const auto& [index, val] : signal) { @@ -634,21 +624,18 @@ Double_t TRestSignalAnalysis::GetSlopeIntegral(const std::vectorsecond; + for (auto it = signal.begin(); it != max; ++it) sum += it->second; return sum; } /////////////////////////////////////////////// /// \brief It returns the slope of the first positive -/// rise (risetime) over a vector of pairs, that should +/// rise (risetime) over a vector of pairs, that should /// correspond to the points over threshold for a given signal. /// -Double_t TRestSignalAnalysis::GetRiseSlope(const std::vector > & signal ) { - - if (signal.size() < 2) - return 0; +Double_t TRestSignalAnalysis::GetRiseSlope(const std::vector >& signal) { + if (signal.size() < 2) return 0; auto max = std::max_element(std::begin(signal), std::end(signal), [](const auto& p1, const auto& p2) { return p1.second < p2.second; }); @@ -663,13 +650,12 @@ Double_t TRestSignalAnalysis::GetRiseSlope(const std::vector > & signal ) { - - if (signal.size() < 2) { - return 0; +Double_t TRestSignalAnalysis::GetRiseTime(const std::vector >& signal) { + if (signal.size() < 2) { + return 0; } auto max = std::max_element(std::begin(signal), std::end(signal), From 6467f33f359621170f64906f77db744759aa9239 Mon Sep 17 00:00:00 2001 From: juanan Date: Fri, 31 Mar 2023 11:40:28 +0200 Subject: [PATCH 11/21] Adding new function to get maxBin and minBin in a given range --- .../framework/analysis/inc/TRestSignalAnalysis.h | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/source/framework/analysis/inc/TRestSignalAnalysis.h b/source/framework/analysis/inc/TRestSignalAnalysis.h index 3d94ea651..2fdbc9f2c 100644 --- a/source/framework/analysis/inc/TRestSignalAnalysis.h +++ b/source/framework/analysis/inc/TRestSignalAnalysis.h @@ -67,11 +67,27 @@ Int_t GetMaxBin(const std::vector& signal) { return std::distance(signal.begin(), std::max_element(signal.begin(), signal.end())); } +template +Int_t GetMaxBin(const std::vector& signal, int startBin, int endBin) { + if (endBin <= 0 || endBin > (int)signal.size()) endBin= signal.size(); + if (startBin < 0) startBin = 0; + + return std::distance(signal.begin(), std::max_element(signal.begin()+startBin, signal.begin()+endBin)); +} + template Int_t GetMinBin(const std::vector& signal) { return std::distance(signal.begin(), std::min_element(signal.begin(), signal.end())); } +template +Int_t GetMinBin(const std::vector& signal, int startBin, int endBin) { + if (endBin <= 0 || endBin > (int)signal.size()) endBin= signal.size(); + if (startBin < 0) startBin = 0; + + return std::distance(signal.begin(), std::min_element(signal.begin()+startBin, signal.begin()+endBin)); +} + template Double_t GetIntegral(const std::vector& signal, Int_t startBin, Int_t endBin); From b545262887bd80441dc0feb91da94233267ac232 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 31 Mar 2023 09:41:06 +0000 Subject: [PATCH 12/21] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- source/framework/analysis/inc/TRestSignalAnalysis.h | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/source/framework/analysis/inc/TRestSignalAnalysis.h b/source/framework/analysis/inc/TRestSignalAnalysis.h index 2fdbc9f2c..c00596441 100644 --- a/source/framework/analysis/inc/TRestSignalAnalysis.h +++ b/source/framework/analysis/inc/TRestSignalAnalysis.h @@ -69,10 +69,11 @@ Int_t GetMaxBin(const std::vector& signal) { template Int_t GetMaxBin(const std::vector& signal, int startBin, int endBin) { - if (endBin <= 0 || endBin > (int)signal.size()) endBin= signal.size(); + if (endBin <= 0 || endBin > (int)signal.size()) endBin = signal.size(); if (startBin < 0) startBin = 0; - return std::distance(signal.begin(), std::max_element(signal.begin()+startBin, signal.begin()+endBin)); + return std::distance(signal.begin(), + std::max_element(signal.begin() + startBin, signal.begin() + endBin)); } template @@ -82,10 +83,11 @@ Int_t GetMinBin(const std::vector& signal) { template Int_t GetMinBin(const std::vector& signal, int startBin, int endBin) { - if (endBin <= 0 || endBin > (int)signal.size()) endBin= signal.size(); + if (endBin <= 0 || endBin > (int)signal.size()) endBin = signal.size(); if (startBin < 0) startBin = 0; - return std::distance(signal.begin(), std::min_element(signal.begin()+startBin, signal.begin()+endBin)); + return std::distance(signal.begin(), + std::min_element(signal.begin() + startBin, signal.begin() + endBin)); } template From 18135b3a251790e1e232ade78f8450aa765cc0b7 Mon Sep 17 00:00:00 2001 From: juanan Date: Fri, 31 Mar 2023 17:30:08 +0200 Subject: [PATCH 13/21] Addressing pipeline failure --- .../analysis/src/TRestSignalAnalysis.cxx | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/source/framework/analysis/src/TRestSignalAnalysis.cxx b/source/framework/analysis/src/TRestSignalAnalysis.cxx index 30fbeb41d..508a66510 100644 --- a/source/framework/analysis/src/TRestSignalAnalysis.cxx +++ b/source/framework/analysis/src/TRestSignalAnalysis.cxx @@ -615,17 +615,19 @@ Double_t TRestSignalAnalysis::GetTripleMaxIntegral(TGraph* signal) { /// Double_t TRestSignalAnalysis::GetSlopeIntegral(const std::vector >& signal) { Double_t sum = 0; - /*Double_t pVal = 0; + Double_t pVal = 0; for (const auto& [index, val] : signal) { - if (val - pVal < 0) break; + if (val - pVal < 0){ + break; + } sum += val; pVal = val; - }*/ - auto max = std::max_element(std::begin(signal), std::end(signal), + } +/* auto max = std::max_element(std::begin(signal), std::end(signal), [](const auto& p1, const auto& p2) { return p1.second < p2.second; }); for (auto it = signal.begin(); it != max; ++it) sum += it->second; - +*/ return sum; } @@ -641,11 +643,12 @@ Double_t TRestSignalAnalysis::GetRiseSlope(const std::vectorfirst; + --maxBin; auto startBin = signal.front().first; Double_t hP = max->second; Double_t lP = signal.front().second; - return (hP - lP) / (maxBin - startBin - 1); + return (hP - lP) / (maxBin - startBin); } /////////////////////////////////////////////// From e30118f21e3af1607231f47edac5cbacf80d13d9 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 31 Mar 2023 15:30:48 +0000 Subject: [PATCH 14/21] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- .../framework/analysis/src/TRestSignalAnalysis.cxx | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/source/framework/analysis/src/TRestSignalAnalysis.cxx b/source/framework/analysis/src/TRestSignalAnalysis.cxx index 508a66510..bc9374bea 100644 --- a/source/framework/analysis/src/TRestSignalAnalysis.cxx +++ b/source/framework/analysis/src/TRestSignalAnalysis.cxx @@ -617,17 +617,17 @@ Double_t TRestSignalAnalysis::GetSlopeIntegral(const std::vectorsecond; -*/ + for (auto it = signal.begin(); it != max; ++it) sum += it->second; + */ return sum; } From 70e5a9a065991aaace2a5871e803dc465bd78909 Mon Sep 17 00:00:00 2001 From: juanan Date: Fri, 31 Mar 2023 23:32:21 +0200 Subject: [PATCH 15/21] Simplifiying points over threshold calculation inside TRestSignalAnalysis, now a vector with the baseline substracted has to be provided --- .../analysis/inc/TRestSignalAnalysis.h | 15 +++------------ .../analysis/src/TRestSignalAnalysis.cxx | 17 ++++++++--------- 2 files changed, 11 insertions(+), 21 deletions(-) diff --git a/source/framework/analysis/inc/TRestSignalAnalysis.h b/source/framework/analysis/inc/TRestSignalAnalysis.h index c00596441..6bf842866 100644 --- a/source/framework/analysis/inc/TRestSignalAnalysis.h +++ b/source/framework/analysis/inc/TRestSignalAnalysis.h @@ -60,15 +60,11 @@ template std::vector > GetPointsOverThreshold(const std::vector& signal, TVector2& range, const TVector2& thrPar, Int_t nPointsOver, Int_t nPointsFlat, - Double_t baseLine, Double_t baseLineSigma); + Double_t baseLineSigma); -template -Int_t GetMaxBin(const std::vector& signal) { - return std::distance(signal.begin(), std::max_element(signal.begin(), signal.end())); -} template -Int_t GetMaxBin(const std::vector& signal, int startBin, int endBin) { +Int_t GetMaxBin(const std::vector& signal, int startBin = 0, int endBin = 0) { if (endBin <= 0 || endBin > (int)signal.size()) endBin = signal.size(); if (startBin < 0) startBin = 0; @@ -77,12 +73,7 @@ Int_t GetMaxBin(const std::vector& signal, int startBin, int endBin) { } template -Int_t GetMinBin(const std::vector& signal) { - return std::distance(signal.begin(), std::min_element(signal.begin(), signal.end())); -} - -template -Int_t GetMinBin(const std::vector& signal, int startBin, int endBin) { +Int_t GetMinBin(const std::vector& signal, int startBin = 0, int endBin = 0) { if (endBin <= 0 || endBin > (int)signal.size()) endBin = signal.size(); if (startBin < 0) startBin = 0; diff --git a/source/framework/analysis/src/TRestSignalAnalysis.cxx b/source/framework/analysis/src/TRestSignalAnalysis.cxx index bc9374bea..91a8d8041 100644 --- a/source/framework/analysis/src/TRestSignalAnalysis.cxx +++ b/source/framework/analysis/src/TRestSignalAnalysis.cxx @@ -283,7 +283,7 @@ template std::vector TRestSignalAnalysis::GetDerivative(const std::vect template std::vector > TRestSignalAnalysis::GetPointsOverThreshold( const std::vector& signal, TVector2& range, const TVector2& thrPar, Int_t nPointsOver, - Int_t nPointsFlat, Double_t baseLine, Double_t baseLineSigma) { + Int_t nPointsFlat, Double_t baseLineSigma) { if (range.X() < 0) range.SetX(0); if (range.Y() <= 0) range.SetY(signal.size()); @@ -296,10 +296,10 @@ std::vector > TRestSignalAnalysis::GetPointsOverThre for (int i = range.X(); i < range.Y(); i++) { // Filling a pulse with consecutive points that are over threshold - if ((signal[i] - baseLine) > threshold) { + if (signal[i] > threshold) { int pos = i; std::vector pulse; - pulse.push_back(signal[i] - baseLine); + pulse.push_back(signal[i]); i++; // If the pulse ends in a flat end above the threshold, the parameter @@ -308,7 +308,7 @@ std::vector > TRestSignalAnalysis::GetPointsOverThre // decision to cut this anomalous behaviour. And all points over threshold // will be added to the pulse vector. int flatN = 0; - while (i < range.Y() && (signal[i] - baseLine) > threshold) { + while (i < range.Y() && signal[i] > threshold) { if (TMath::Abs(signal[i] - signal[i - 1]) > threshold) { flatN = 0; } else { @@ -316,7 +316,7 @@ std::vector > TRestSignalAnalysis::GetPointsOverThre } if (flatN < nPointsFlat) { - pulse.push_back(signal[i] - baseLine); + pulse.push_back(signal[i]); i++; } else { break; @@ -329,8 +329,7 @@ std::vector > TRestSignalAnalysis::GetPointsOverThre double stdev = std::sqrt(sq_sum / pulse.size() - mean * mean); if (stdev > signalTh * baseLineSigma) - for (unsigned int j = 0; j < pulse.size(); j++) - pointsOverThreshold.push_back(std::make_pair(pos + j, pulse[j])); + for (int j = pos; j < i; j++) pointsOverThreshold.push_back(std::make_pair(pos + j, pulse[j])); } } } @@ -340,10 +339,10 @@ std::vector > TRestSignalAnalysis::GetPointsOverThre template std::vector > TRestSignalAnalysis::GetPointsOverThreshold( const std::vector& signal, TVector2& range, const TVector2& thrPar, Int_t nPointsOver, - Int_t nPointsFlat, Double_t baseLine, Double_t baseLineSigma); + Int_t nPointsFlat, Double_t baseLineSigma); template std::vector > TRestSignalAnalysis::GetPointsOverThreshold( const std::vector& signal, TVector2& range, const TVector2& thrPar, Int_t nPointsOver, - Int_t nPointsFlat, Double_t baseLine, Double_t baseLineSigma); + Int_t nPointsFlat, Double_t baseLineSigma); /////////////////////////////////////////////// /// \brief It returns the integral of the signal in the From 0e4c1e7eed26252e4329d3966617d597e2143536 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 31 Mar 2023 21:32:44 +0000 Subject: [PATCH 16/21] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- source/framework/analysis/inc/TRestSignalAnalysis.h | 1 - source/framework/analysis/src/TRestSignalAnalysis.cxx | 3 ++- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/source/framework/analysis/inc/TRestSignalAnalysis.h b/source/framework/analysis/inc/TRestSignalAnalysis.h index 6bf842866..d25135f57 100644 --- a/source/framework/analysis/inc/TRestSignalAnalysis.h +++ b/source/framework/analysis/inc/TRestSignalAnalysis.h @@ -62,7 +62,6 @@ std::vector > GetPointsOverThreshold(const std::vect Int_t nPointsOver, Int_t nPointsFlat, Double_t baseLineSigma); - template Int_t GetMaxBin(const std::vector& signal, int startBin = 0, int endBin = 0) { if (endBin <= 0 || endBin > (int)signal.size()) endBin = signal.size(); diff --git a/source/framework/analysis/src/TRestSignalAnalysis.cxx b/source/framework/analysis/src/TRestSignalAnalysis.cxx index 91a8d8041..f50b48929 100644 --- a/source/framework/analysis/src/TRestSignalAnalysis.cxx +++ b/source/framework/analysis/src/TRestSignalAnalysis.cxx @@ -329,7 +329,8 @@ std::vector > TRestSignalAnalysis::GetPointsOverThre double stdev = std::sqrt(sq_sum / pulse.size() - mean * mean); if (stdev > signalTh * baseLineSigma) - for (int j = pos; j < i; j++) pointsOverThreshold.push_back(std::make_pair(pos + j, pulse[j])); + for (int j = pos; j < i; j++) + pointsOverThreshold.push_back(std::make_pair(pos + j, pulse[j])); } } } From 70550f41c17ba082dc82aba76ffb82397336bbe2 Mon Sep 17 00:00:00 2001 From: juanan Date: Sat, 1 Apr 2023 00:26:07 +0200 Subject: [PATCH 17/21] Addressing bug --- source/framework/analysis/src/TRestSignalAnalysis.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source/framework/analysis/src/TRestSignalAnalysis.cxx b/source/framework/analysis/src/TRestSignalAnalysis.cxx index f50b48929..cd1939dee 100644 --- a/source/framework/analysis/src/TRestSignalAnalysis.cxx +++ b/source/framework/analysis/src/TRestSignalAnalysis.cxx @@ -329,7 +329,7 @@ std::vector > TRestSignalAnalysis::GetPointsOverThre double stdev = std::sqrt(sq_sum / pulse.size() - mean * mean); if (stdev > signalTh * baseLineSigma) - for (int j = pos; j < i; j++) + for (unsigned int j = 0; j < pulse.size(); j++) pointsOverThreshold.push_back(std::make_pair(pos + j, pulse[j])); } } From 894a54b3a80106ef305b12d49f34f797c0399576 Mon Sep 17 00:00:00 2001 From: juanan Date: Mon, 3 Apr 2023 20:30:38 +0200 Subject: [PATCH 18/21] Renaming TRestSignalAnalysis to TRestPulseShapeAnalysis --- ...alAnalysis.h => TRestPulseShapeAnalysis.h} | 8 +- ...alysis.cxx => TRestPulseShapeAnalysis.cxx} | 88 +++++++++---------- 2 files changed, 48 insertions(+), 48 deletions(-) rename source/framework/analysis/inc/{TRestSignalAnalysis.h => TRestPulseShapeAnalysis.h} (96%) rename source/framework/analysis/src/{TRestSignalAnalysis.cxx => TRestPulseShapeAnalysis.cxx} (85%) diff --git a/source/framework/analysis/inc/TRestSignalAnalysis.h b/source/framework/analysis/inc/TRestPulseShapeAnalysis.h similarity index 96% rename from source/framework/analysis/inc/TRestSignalAnalysis.h rename to source/framework/analysis/inc/TRestPulseShapeAnalysis.h index d25135f57..c9dd60be8 100644 --- a/source/framework/analysis/inc/TRestSignalAnalysis.h +++ b/source/framework/analysis/inc/TRestPulseShapeAnalysis.h @@ -20,8 +20,8 @@ * For the list of contributors see $REST_PATH/CREDITS. * *************************************************************************/ -#ifndef RestCore_TRestSignalAnalysis -#define RestCore_TRestSignalAnalysis +#ifndef RestCore_TRestPulseShapeAnalysis +#define RestCore_TRestPulseShapeAnalysis #include #include @@ -33,7 +33,7 @@ #include /// This namespace define generic functions to calculate different signal parameters -namespace TRestSignalAnalysis { +namespace TRestPulseShapeAnalysis { template void CalculateBaselineAndSigmaSD(const std::vector& signal, Int_t startBin, Int_t endBin, @@ -98,6 +98,6 @@ TVector2 GetMaxGauss(TGraph* signal); TVector2 GetMaxLandau(TGraph* signal); TVector2 GetMaxAget(TGraph* signal); -} // namespace TRestSignalAnalysis +} // namespace TRestPulseShapeAnalysis #endif diff --git a/source/framework/analysis/src/TRestSignalAnalysis.cxx b/source/framework/analysis/src/TRestPulseShapeAnalysis.cxx similarity index 85% rename from source/framework/analysis/src/TRestSignalAnalysis.cxx rename to source/framework/analysis/src/TRestPulseShapeAnalysis.cxx index cd1939dee..4ec3b3c3f 100644 --- a/source/framework/analysis/src/TRestSignalAnalysis.cxx +++ b/source/framework/analysis/src/TRestPulseShapeAnalysis.cxx @@ -21,7 +21,7 @@ *************************************************************************/ ////////////////////////////////////////////////////////////////////////// -/// TRestSignalAnalysis defines several functions to calculate different +/// TRestPulseShapeAnalysis defines several functions to calculate different /// signal parameters. ///-------------------------------------------------------------------------- /// @@ -32,7 +32,7 @@ /// 2022-December First implementation /// JuanAn Garcia /// -/// \class TRestSignalAnalysis +/// \class TRestPulseShapeAnalysis /// \author: JuanAn Garcia e-mail: juanangp@unizar.es /// ///
@@ -40,7 +40,7 @@ #include #include -#include +#include /////////////////////////////////////////////// /// \brief This method is used to determine the value @@ -49,7 +49,7 @@ /// The baseline sigma is determined as the standard deviation /// of the baseline in range provided. template -void TRestSignalAnalysis::CalculateBaselineAndSigmaSD(const std::vector& signal, Int_t startBin, +void TRestPulseShapeAnalysis::CalculateBaselineAndSigmaSD(const std::vector& signal, Int_t startBin, Int_t endBin, Double_t& baseLine, Double_t& baseLineSigma) { baseLine = 0; @@ -69,11 +69,11 @@ void TRestSignalAnalysis::CalculateBaselineAndSigmaSD(const std::vector& sign baseLineSigma = TMath::Sqrt(baseLineSigma / (double)nPoints - baseLine * baseLine); } } -template void TRestSignalAnalysis::CalculateBaselineAndSigmaSD(const std::vector& signal, +template void TRestPulseShapeAnalysis::CalculateBaselineAndSigmaSD(const std::vector& signal, Int_t startBin, Int_t endBin, Double_t& baseLine, Double_t& baseLineSigma); -template void TRestSignalAnalysis::CalculateBaselineAndSigmaSD(const std::vector& signal, +template void TRestPulseShapeAnalysis::CalculateBaselineAndSigmaSD(const std::vector& signal, Int_t startBin, Int_t endBin, Double_t& baseLine, Double_t& baseLineSigma); @@ -85,7 +85,7 @@ template void TRestSignalAnalysis::CalculateBaselineAndSigmaSD(const st /// range (IQR) in the baseline range provided. The IQR /// is more robust towards outliers than the standard deviation. template -void TRestSignalAnalysis::CalculateBaselineAndSigmaIQR(const std::vector& signal, Int_t startBin, +void TRestPulseShapeAnalysis::CalculateBaselineAndSigmaIQR(const std::vector& signal, Int_t startBin, Int_t endBin, Double_t& baseLine, Double_t& baseLineSigma) { baseLine = 0; @@ -108,11 +108,11 @@ void TRestSignalAnalysis::CalculateBaselineAndSigmaIQR(const std::vector& sig 1.349; // IQR/1.349 equals the standard deviation in case of normally distributed data } -template void TRestSignalAnalysis::CalculateBaselineAndSigmaIQR(const std::vector& signal, +template void TRestPulseShapeAnalysis::CalculateBaselineAndSigmaIQR(const std::vector& signal, Int_t startBin, Int_t endBin, Double_t& baseLine, Double_t& baseLineSigma); -template void TRestSignalAnalysis::CalculateBaselineAndSigmaIQR(const std::vector& signal, +template void TRestPulseShapeAnalysis::CalculateBaselineAndSigmaIQR(const std::vector& signal, Int_t startBin, Int_t endBin, Double_t& baseLine, Double_t& baseLineSigma); @@ -121,7 +121,7 @@ template void TRestSignalAnalysis::CalculateBaselineAndSigmaIQR(const s /// the data points in a given range defined /// between startBin and endBin template -Double_t TRestSignalAnalysis::GetAverage(const std::vector& signal, Int_t startBin, Int_t endBin) { +Double_t TRestPulseShapeAnalysis::GetAverage(const std::vector& signal, Int_t startBin, Int_t endBin) { int nPoints = 0; Double_t avg = 0; @@ -136,9 +136,9 @@ Double_t TRestSignalAnalysis::GetAverage(const std::vector& signal, Int_t sta return avg; } -template Double_t TRestSignalAnalysis::GetAverage(const std::vector& signal, Int_t startBin, +template Double_t TRestPulseShapeAnalysis::GetAverage(const std::vector& signal, Int_t startBin, Int_t endBin); -template Double_t TRestSignalAnalysis::GetAverage(const std::vector& signal, Int_t startBin, +template Double_t TRestPulseShapeAnalysis::GetAverage(const std::vector& signal, Int_t startBin, Int_t endBin); /////////////////////////////////////////////// @@ -149,14 +149,14 @@ template Double_t TRestSignalAnalysis::GetAverage(const std::vector -std::vector TRestSignalAnalysis::GetSignalSmoothed(const std::vector& signal, +std::vector TRestPulseShapeAnalysis::GetSignalSmoothed(const std::vector& signal, int averagingPoints) { const int pulseDepth = signal.size(); std::vector smoothed(pulseDepth, 0); averagingPoints = (averagingPoints / 2) * 2 + 1; // make it odd >= averagingPoints - Float_t sumAvg = TRestSignalAnalysis::GetAverage(signal, 0, averagingPoints); + Float_t sumAvg = TRestPulseShapeAnalysis::GetAverage(signal, 0, averagingPoints); // Points at the beginning, where we can calculate a moving average for (int i = 0; i <= averagingPoints / 2; i++) smoothed[i] = sumAvg; @@ -174,9 +174,9 @@ std::vector TRestSignalAnalysis::GetSignalSmoothed(const std::vector return smoothed; } -template std::vector TRestSignalAnalysis::GetSignalSmoothed( +template std::vector TRestPulseShapeAnalysis::GetSignalSmoothed( const std::vector& signal, int averagingPoints); -template std::vector TRestSignalAnalysis::GetSignalSmoothed( +template std::vector TRestPulseShapeAnalysis::GetSignalSmoothed( const std::vector& signal, int averagingPoints); /////////////////////////////////////////////// @@ -188,7 +188,7 @@ template std::vector TRestSignalAnalysis::GetSignalSmoothed( /// neighbouring points used to average the signal /// template -std::vector TRestSignalAnalysis::GetSignalSmoothed_ExcludeOutliers(const std::vector& signal, +std::vector TRestPulseShapeAnalysis::GetSignalSmoothed_ExcludeOutliers(const std::vector& signal, int averagingPoints, Double_t& baseLine, Double_t& baseLineSigma) { @@ -199,7 +199,7 @@ std::vector TRestSignalAnalysis::GetSignalSmoothed_ExcludeOutliers(cons averagingPoints = (averagingPoints / 2) * 2 + 1; // make it odd >= averagingPoints - Float_t sumAvg = TRestSignalAnalysis::GetAverage(signal, 0, averagingPoints); + Float_t sumAvg = TRestPulseShapeAnalysis::GetAverage(signal, 0, averagingPoints); // Points at the beginning, where we can calculate a moving average for (int i = 0; i <= averagingPoints / 2; i++) smoothed[i] = sumAvg; @@ -224,9 +224,9 @@ std::vector TRestSignalAnalysis::GetSignalSmoothed_ExcludeOutliers(cons return smoothed; } -template std::vector TRestSignalAnalysis::GetSignalSmoothed_ExcludeOutliers( +template std::vector TRestPulseShapeAnalysis::GetSignalSmoothed_ExcludeOutliers( const std::vector& signal, int averagingPoints, Double_t& baseLine, Double_t& baseLineSigma); -template std::vector TRestSignalAnalysis::GetSignalSmoothed_ExcludeOutliers( +template std::vector TRestPulseShapeAnalysis::GetSignalSmoothed_ExcludeOutliers( const std::vector& signal, int averagingPoints, Double_t& baseLine, Double_t& baseLineSigma); /////////////////////////////////////////////// @@ -234,7 +234,7 @@ template std::vector TRestSignalAnalysis::GetSignalSmoothed_ExcludeOutl /// points /// template -std::vector TRestSignalAnalysis::GetDerivative(const std::vector& signal) { +std::vector TRestPulseShapeAnalysis::GetDerivative(const std::vector& signal) { std::vector derivative(0, signal.size() - 1); for (size_t i = 0; i < signal.size() - 1; i++) { @@ -244,8 +244,8 @@ std::vector TRestSignalAnalysis::GetDerivative(const std::vector& si return derivative; } -template std::vector TRestSignalAnalysis::GetDerivative(const std::vector& signal); -template std::vector TRestSignalAnalysis::GetDerivative(const std::vector& signal); +template std::vector TRestPulseShapeAnalysis::GetDerivative(const std::vector& signal); +template std::vector TRestPulseShapeAnalysis::GetDerivative(const std::vector& signal); /////////////////////////////////////////////// /// \brief It returns a vector of the data points @@ -281,7 +281,7 @@ template std::vector TRestSignalAnalysis::GetDerivative(const std::vect /// template -std::vector > TRestSignalAnalysis::GetPointsOverThreshold( +std::vector > TRestPulseShapeAnalysis::GetPointsOverThreshold( const std::vector& signal, TVector2& range, const TVector2& thrPar, Int_t nPointsOver, Int_t nPointsFlat, Double_t baseLineSigma) { if (range.X() < 0) range.SetX(0); @@ -338,10 +338,10 @@ std::vector > TRestSignalAnalysis::GetPointsOverThre return pointsOverThreshold; } -template std::vector > TRestSignalAnalysis::GetPointsOverThreshold( +template std::vector > TRestPulseShapeAnalysis::GetPointsOverThreshold( const std::vector& signal, TVector2& range, const TVector2& thrPar, Int_t nPointsOver, Int_t nPointsFlat, Double_t baseLineSigma); -template std::vector > TRestSignalAnalysis::GetPointsOverThreshold( +template std::vector > TRestPulseShapeAnalysis::GetPointsOverThreshold( const std::vector& signal, TVector2& range, const TVector2& thrPar, Int_t nPointsOver, Int_t nPointsFlat, Double_t baseLineSigma); @@ -350,16 +350,16 @@ template std::vector > TRestSignalAnalysis::GetPoint /// range passed as argument /// template -Double_t TRestSignalAnalysis::GetIntegral(const std::vector& signal, Int_t startBin, Int_t endBin) { +Double_t TRestPulseShapeAnalysis::GetIntegral(const std::vector& signal, Int_t startBin, Int_t endBin) { Double_t sum = 0; for (int i = startBin; i < endBin; i++) if (i > 0 && i < (int)signal.size()) sum += signal[i]; return sum; } -template Double_t TRestSignalAnalysis::GetIntegral(const std::vector& signal, Int_t startBin, +template Double_t TRestPulseShapeAnalysis::GetIntegral(const std::vector& signal, Int_t startBin, Int_t endBin); -template Double_t TRestSignalAnalysis::GetIntegral(const std::vector& signal, Int_t startBin, +template Double_t TRestPulseShapeAnalysis::GetIntegral(const std::vector& signal, Int_t startBin, Int_t endBin); /////////////////////////////////////////////// @@ -368,7 +368,7 @@ template Double_t TRestSignalAnalysis::GetIntegral(const std::vector& s /// maximum /// template -Double_t TRestSignalAnalysis::GetMaxPeakWidth(const std::vector& signal) { +Double_t TRestPulseShapeAnalysis::GetMaxPeakWidth(const std::vector& signal) { Int_t maxIndex = GetMaxBin(signal); Double_t maxValue = signal[maxIndex]; @@ -389,8 +389,8 @@ Double_t TRestSignalAnalysis::GetMaxPeakWidth(const std::vector& signal) { return rightIndex - leftIndex; } -template Double_t TRestSignalAnalysis::GetMaxPeakWidth(const std::vector& signal); -template Double_t TRestSignalAnalysis::GetMaxPeakWidth(const std::vector& signal); +template Double_t TRestPulseShapeAnalysis::GetMaxPeakWidth(const std::vector& signal); +template Double_t TRestPulseShapeAnalysis::GetMaxPeakWidth(const std::vector& signal); /////////////////////////////////////////////// /// \brief It performs a gaussian fit to the signal @@ -398,7 +398,7 @@ template Double_t TRestSignalAnalysis::GetMaxPeakWidth(const std::vectorGetN(), signal->GetY()); Double_t maxTime = signal->GetPointX(maxBin); Double_t gaussMax = -1, gaussMean = -1; @@ -439,7 +439,7 @@ TVector2 TRestSignalAnalysis::GetMaxGauss(TGraph* signal) { /// a TVector 2 with the maximum and the mean of the /// landau fit /// -TVector2 TRestSignalAnalysis::GetMaxLandau(TGraph* signal) { +TVector2 TRestPulseShapeAnalysis::GetMaxLandau(TGraph* signal) { Int_t maxBin = TMath::LocMax(signal->GetN(), signal->GetY()); Double_t maxTime = signal->GetPointX(maxBin); Double_t landauMax = -1, landauMean = -1; @@ -480,7 +480,7 @@ TVector2 TRestSignalAnalysis::GetMaxLandau(TGraph* signal) { /// TVector 2 with the maximum and the mean of the /// gaussian fit /// -TVector2 TRestSignalAnalysis::GetMaxAget(TGraph* signal) { +TVector2 TRestPulseShapeAnalysis::GetMaxAget(TGraph* signal) { Int_t maxBin = TMath::LocMax(signal->GetN(), signal->GetY()); Double_t maxTime = signal->GetPointX(maxBin); Double_t agetMax = -1, agetMean = -1; @@ -524,7 +524,7 @@ TVector2 TRestSignalAnalysis::GetMaxAget(TGraph* signal) { /// passed as argument. It returns a vector of pairs /// with the integrated window time and energy (charge) /// -std::vector > TRestSignalAnalysis::GetIntWindow(TGraph* signal, double intWindow) { +std::vector > TRestPulseShapeAnalysis::GetIntWindow(TGraph* signal, double intWindow) { const int nPoints = signal->GetN(); std::map > windowMap; @@ -556,7 +556,7 @@ std::vector > TRestSignalAnalysis::GetIntWindow(TGraph /// the energy of the maximum and the neighbouring /// points(three points in total). /// -std::array, 3> TRestSignalAnalysis::GetTripleMax(TGraph* signal) { +std::array, 3> TRestPulseShapeAnalysis::GetTripleMax(TGraph* signal) { Int_t maxBin = TMath::LocMax(signal->GetN(), signal->GetY()); std::array, 3> tripleMax; @@ -578,8 +578,8 @@ std::array, 3> TRestSignalAnalysis::GetTripleMax(T /// a TGraph and returns a TVector2 with the average /// time and energy (charge) /// -TVector2 TRestSignalAnalysis::GetTripleMaxAverage(TGraph* signal) { - auto tripleMax = TRestSignalAnalysis::GetTripleMax(signal); +TVector2 TRestPulseShapeAnalysis::GetTripleMaxAverage(TGraph* signal) { + auto tripleMax = TRestPulseShapeAnalysis::GetTripleMax(signal); double eAvg = 0; double hitTimeAvg = 0; for (const auto& [hitTime, energy] : tripleMax) { @@ -598,8 +598,8 @@ TVector2 TRestSignalAnalysis::GetTripleMaxAverage(TGraph* signal) { /// a TGraph and returns the addition of the maximum plus /// the neigbouring bins. /// -Double_t TRestSignalAnalysis::GetTripleMaxIntegral(TGraph* signal) { - auto tripleMax = TRestSignalAnalysis::GetTripleMax(signal); +Double_t TRestPulseShapeAnalysis::GetTripleMaxIntegral(TGraph* signal) { + auto tripleMax = TRestPulseShapeAnalysis::GetTripleMax(signal); double totEnergy = 0; for (const auto& [hitTime, energy] : tripleMax) { totEnergy += energy; @@ -613,7 +613,7 @@ Double_t TRestSignalAnalysis::GetTripleMaxIntegral(TGraph* signal) { /// rise (risetime) over a vector of pairs, that should /// correspond to the points over threshold for a given signal. /// -Double_t TRestSignalAnalysis::GetSlopeIntegral(const std::vector >& signal) { +Double_t TRestPulseShapeAnalysis::GetSlopeIntegral(const std::vector >& signal) { Double_t sum = 0; Double_t pVal = 0; for (const auto& [index, val] : signal) { @@ -636,7 +636,7 @@ Double_t TRestSignalAnalysis::GetSlopeIntegral(const std::vector >& signal) { +Double_t TRestPulseShapeAnalysis::GetRiseSlope(const std::vector >& signal) { if (signal.size() < 2) return 0; auto max = std::max_element(std::begin(signal), std::end(signal), @@ -656,7 +656,7 @@ Double_t TRestSignalAnalysis::GetRiseSlope(const std::vector >& signal) { +Double_t TRestPulseShapeAnalysis::GetRiseTime(const std::vector >& signal) { if (signal.size() < 2) { return 0; } From 0fda20ce0dc04146c6de583130d7da8f6a2d9534 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 3 Apr 2023 18:30:58 +0000 Subject: [PATCH 19/21] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- .../analysis/src/TRestPulseShapeAnalysis.cxx | 59 +++++++++---------- 1 file changed, 28 insertions(+), 31 deletions(-) diff --git a/source/framework/analysis/src/TRestPulseShapeAnalysis.cxx b/source/framework/analysis/src/TRestPulseShapeAnalysis.cxx index 4ec3b3c3f..f5ad755ee 100644 --- a/source/framework/analysis/src/TRestPulseShapeAnalysis.cxx +++ b/source/framework/analysis/src/TRestPulseShapeAnalysis.cxx @@ -50,8 +50,8 @@ /// of the baseline in range provided. template void TRestPulseShapeAnalysis::CalculateBaselineAndSigmaSD(const std::vector& signal, Int_t startBin, - Int_t endBin, Double_t& baseLine, - Double_t& baseLineSigma) { + Int_t endBin, Double_t& baseLine, + Double_t& baseLineSigma) { baseLine = 0; baseLineSigma = 0; @@ -69,14 +69,12 @@ void TRestPulseShapeAnalysis::CalculateBaselineAndSigmaSD(const std::vector& baseLineSigma = TMath::Sqrt(baseLineSigma / (double)nPoints - baseLine * baseLine); } } -template void TRestPulseShapeAnalysis::CalculateBaselineAndSigmaSD(const std::vector& signal, - Int_t startBin, Int_t endBin, - Double_t& baseLine, - Double_t& baseLineSigma); -template void TRestPulseShapeAnalysis::CalculateBaselineAndSigmaSD(const std::vector& signal, - Int_t startBin, Int_t endBin, - Double_t& baseLine, - Double_t& baseLineSigma); +template void TRestPulseShapeAnalysis::CalculateBaselineAndSigmaSD( + const std::vector& signal, Int_t startBin, Int_t endBin, Double_t& baseLine, + Double_t& baseLineSigma); +template void TRestPulseShapeAnalysis::CalculateBaselineAndSigmaSD( + const std::vector& signal, Int_t startBin, Int_t endBin, Double_t& baseLine, + Double_t& baseLineSigma); /////////////////////////////////////////////// /// \brief This method is used to determine the value /// of the baseline as the median of the data in @@ -86,8 +84,8 @@ template void TRestPulseShapeAnalysis::CalculateBaselineAndSigmaSD(cons /// is more robust towards outliers than the standard deviation. template void TRestPulseShapeAnalysis::CalculateBaselineAndSigmaIQR(const std::vector& signal, Int_t startBin, - Int_t endBin, Double_t& baseLine, - Double_t& baseLineSigma) { + Int_t endBin, Double_t& baseLine, + Double_t& baseLineSigma) { baseLine = 0; baseLineSigma = 0; @@ -108,14 +106,12 @@ void TRestPulseShapeAnalysis::CalculateBaselineAndSigmaIQR(const std::vector& 1.349; // IQR/1.349 equals the standard deviation in case of normally distributed data } -template void TRestPulseShapeAnalysis::CalculateBaselineAndSigmaIQR(const std::vector& signal, - Int_t startBin, Int_t endBin, - Double_t& baseLine, - Double_t& baseLineSigma); -template void TRestPulseShapeAnalysis::CalculateBaselineAndSigmaIQR(const std::vector& signal, - Int_t startBin, Int_t endBin, - Double_t& baseLine, - Double_t& baseLineSigma); +template void TRestPulseShapeAnalysis::CalculateBaselineAndSigmaIQR( + const std::vector& signal, Int_t startBin, Int_t endBin, Double_t& baseLine, + Double_t& baseLineSigma); +template void TRestPulseShapeAnalysis::CalculateBaselineAndSigmaIQR( + const std::vector& signal, Int_t startBin, Int_t endBin, Double_t& baseLine, + Double_t& baseLineSigma); /////////////////////////////////////////////// /// \brief This method performs the average of /// the data points in a given range defined @@ -136,10 +132,10 @@ Double_t TRestPulseShapeAnalysis::GetAverage(const std::vector& signal, Int_t return avg; } -template Double_t TRestPulseShapeAnalysis::GetAverage(const std::vector& signal, Int_t startBin, - Int_t endBin); -template Double_t TRestPulseShapeAnalysis::GetAverage(const std::vector& signal, Int_t startBin, - Int_t endBin); +template Double_t TRestPulseShapeAnalysis::GetAverage(const std::vector& signal, + Int_t startBin, Int_t endBin); +template Double_t TRestPulseShapeAnalysis::GetAverage(const std::vector& signal, + Int_t startBin, Int_t endBin); /////////////////////////////////////////////// /// \brief Return smoothing of signal as the @@ -150,7 +146,7 @@ template Double_t TRestPulseShapeAnalysis::GetAverage(const std::vector /// template std::vector TRestPulseShapeAnalysis::GetSignalSmoothed(const std::vector& signal, - int averagingPoints) { + int averagingPoints) { const int pulseDepth = signal.size(); std::vector smoothed(pulseDepth, 0); @@ -189,9 +185,9 @@ template std::vector TRestPulseShapeAnalysis::GetSignalSmoothed std::vector TRestPulseShapeAnalysis::GetSignalSmoothed_ExcludeOutliers(const std::vector& signal, - int averagingPoints, - Double_t& baseLine, - Double_t& baseLineSigma) { + int averagingPoints, + Double_t& baseLine, + Double_t& baseLineSigma) { const int pulseDepth = signal.size(); std::vector smoothed(pulseDepth, 0); @@ -358,9 +354,9 @@ Double_t TRestPulseShapeAnalysis::GetIntegral(const std::vector& signal, Int_ return sum; } template Double_t TRestPulseShapeAnalysis::GetIntegral(const std::vector& signal, Int_t startBin, - Int_t endBin); + Int_t endBin); template Double_t TRestPulseShapeAnalysis::GetIntegral(const std::vector& signal, Int_t startBin, - Int_t endBin); + Int_t endBin); /////////////////////////////////////////////// /// \brief It returns the width of the pulses @@ -524,7 +520,8 @@ TVector2 TRestPulseShapeAnalysis::GetMaxAget(TGraph* signal) { /// passed as argument. It returns a vector of pairs /// with the integrated window time and energy (charge) /// -std::vector > TRestPulseShapeAnalysis::GetIntWindow(TGraph* signal, double intWindow) { +std::vector > TRestPulseShapeAnalysis::GetIntWindow(TGraph* signal, + double intWindow) { const int nPoints = signal->GetN(); std::map > windowMap; From 5b620cd10bcd38608938096fc186e1ab50d3955b Mon Sep 17 00:00:00 2001 From: juanan Date: Thu, 25 May 2023 13:19:22 +0200 Subject: [PATCH 20/21] =?UTF-8?q?Moving=20GetPointsOverThreshold=20from=20?= =?UTF-8?q?=C2=B4std::pair=C2=B4=20to=20=C2=B4std:?= =?UTF-8?q?:pair=C2=B4?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .../analysis/inc/TRestPulseShapeAnalysis.h | 11 +++++--- .../analysis/src/TRestPulseShapeAnalysis.cxx | 26 +++++++++++++------ 2 files changed, 25 insertions(+), 12 deletions(-) diff --git a/source/framework/analysis/inc/TRestPulseShapeAnalysis.h b/source/framework/analysis/inc/TRestPulseShapeAnalysis.h index c9dd60be8..c53685e6d 100644 --- a/source/framework/analysis/inc/TRestPulseShapeAnalysis.h +++ b/source/framework/analysis/inc/TRestPulseShapeAnalysis.h @@ -57,7 +57,7 @@ template std::vector GetDerivative(const std::vector& signal); template -std::vector > GetPointsOverThreshold(const std::vector& signal, +std::vector > GetPointsOverThreshold(const std::vector& signal, TVector2& range, const TVector2& thrPar, Int_t nPointsOver, Int_t nPointsFlat, Double_t baseLineSigma); @@ -86,9 +86,12 @@ Double_t GetIntegral(const std::vector& signal, Int_t startBin, Int_t endBin) template Double_t GetMaxPeakWidth(const std::vector& signal); -Double_t GetSlopeIntegral(const std::vector >& signal); -Double_t GetRiseSlope(const std::vector >& signal); -Double_t GetRiseTime(const std::vector >& signal); +template +Double_t GetSlopeIntegral(const std::vector >& signal); +template +Double_t GetRiseSlope(const std::vector >& signal); +template +Double_t GetRiseTime(const std::vector >& signal); std::vector > GetIntWindow(TGraph* signal, double intWindow); std::array, 3> GetTripleMax(TGraph* signal); diff --git a/source/framework/analysis/src/TRestPulseShapeAnalysis.cxx b/source/framework/analysis/src/TRestPulseShapeAnalysis.cxx index f5ad755ee..2cb072c3a 100644 --- a/source/framework/analysis/src/TRestPulseShapeAnalysis.cxx +++ b/source/framework/analysis/src/TRestPulseShapeAnalysis.cxx @@ -277,13 +277,13 @@ template std::vector TRestPulseShapeAnalysis::GetDerivative(const std:: /// template -std::vector > TRestPulseShapeAnalysis::GetPointsOverThreshold( +std::vector > TRestPulseShapeAnalysis::GetPointsOverThreshold( const std::vector& signal, TVector2& range, const TVector2& thrPar, Int_t nPointsOver, Int_t nPointsFlat, Double_t baseLineSigma) { if (range.X() < 0) range.SetX(0); if (range.Y() <= 0) range.SetY(signal.size()); - std::vector > pointsOverThreshold; + std::vector > pointsOverThreshold; double pointTh = thrPar.X(); double signalTh = thrPar.Y(); @@ -333,11 +333,10 @@ std::vector > TRestPulseShapeAnalysis::GetPointsOver return pointsOverThreshold; } - -template std::vector > TRestPulseShapeAnalysis::GetPointsOverThreshold( +template std::vector > TRestPulseShapeAnalysis::GetPointsOverThreshold( const std::vector& signal, TVector2& range, const TVector2& thrPar, Int_t nPointsOver, Int_t nPointsFlat, Double_t baseLineSigma); -template std::vector > TRestPulseShapeAnalysis::GetPointsOverThreshold( +template std::vector > TRestPulseShapeAnalysis::GetPointsOverThreshold( const std::vector& signal, TVector2& range, const TVector2& thrPar, Int_t nPointsOver, Int_t nPointsFlat, Double_t baseLineSigma); @@ -610,7 +609,8 @@ Double_t TRestPulseShapeAnalysis::GetTripleMaxIntegral(TGraph* signal) { /// rise (risetime) over a vector of pairs, that should /// correspond to the points over threshold for a given signal. /// -Double_t TRestPulseShapeAnalysis::GetSlopeIntegral(const std::vector >& signal) { +template +Double_t TRestPulseShapeAnalysis::GetSlopeIntegral(const std::vector >& signal) { Double_t sum = 0; Double_t pVal = 0; for (const auto& [index, val] : signal) { @@ -627,13 +627,16 @@ Double_t TRestPulseShapeAnalysis::GetSlopeIntegral(const std::vector >& signal); +template Double_t TRestPulseShapeAnalysis::GetSlopeIntegral(const std::vector >& signal); /////////////////////////////////////////////// /// \brief It returns the slope of the first positive /// rise (risetime) over a vector of pairs, that should /// correspond to the points over threshold for a given signal. /// -Double_t TRestPulseShapeAnalysis::GetRiseSlope(const std::vector >& signal) { +template +Double_t TRestPulseShapeAnalysis::GetRiseSlope(const std::vector >& signal) { if (signal.size() < 2) return 0; auto max = std::max_element(std::begin(signal), std::end(signal), @@ -647,13 +650,17 @@ Double_t TRestPulseShapeAnalysis::GetRiseSlope(const std::vector >& signal); +template Double_t TRestPulseShapeAnalysis::GetRiseSlope(const std::vector >& signal); + /////////////////////////////////////////////// /// \brief It returns the time of the first positive /// rise or risetime over a vector of pairs, that should /// correspond to the points over threshold for a given signal. /// -Double_t TRestPulseShapeAnalysis::GetRiseTime(const std::vector >& signal) { +template +Double_t TRestPulseShapeAnalysis::GetRiseTime(const std::vector >& signal) { if (signal.size() < 2) { return 0; } @@ -665,3 +672,6 @@ Double_t TRestPulseShapeAnalysis::GetRiseTime(const std::vector >& signal); +template Double_t TRestPulseShapeAnalysis::GetRiseTime(const std::vector >& signal); + From 47ed4ff5e4b0a2c19f4cca36001d580294c21c38 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 25 May 2023 11:20:48 +0000 Subject: [PATCH 21/21] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- .../analysis/inc/TRestPulseShapeAnalysis.h | 7 +++---- .../analysis/src/TRestPulseShapeAnalysis.cxx | 17 ++++++++++------- 2 files changed, 13 insertions(+), 11 deletions(-) diff --git a/source/framework/analysis/inc/TRestPulseShapeAnalysis.h b/source/framework/analysis/inc/TRestPulseShapeAnalysis.h index c53685e6d..b2d2b3848 100644 --- a/source/framework/analysis/inc/TRestPulseShapeAnalysis.h +++ b/source/framework/analysis/inc/TRestPulseShapeAnalysis.h @@ -57,10 +57,9 @@ template std::vector GetDerivative(const std::vector& signal); template -std::vector > GetPointsOverThreshold(const std::vector& signal, - TVector2& range, const TVector2& thrPar, - Int_t nPointsOver, Int_t nPointsFlat, - Double_t baseLineSigma); +std::vector > GetPointsOverThreshold(const std::vector& signal, TVector2& range, + const TVector2& thrPar, Int_t nPointsOver, + Int_t nPointsFlat, Double_t baseLineSigma); template Int_t GetMaxBin(const std::vector& signal, int startBin = 0, int endBin = 0) { diff --git a/source/framework/analysis/src/TRestPulseShapeAnalysis.cxx b/source/framework/analysis/src/TRestPulseShapeAnalysis.cxx index 2cb072c3a..c5a0e8745 100644 --- a/source/framework/analysis/src/TRestPulseShapeAnalysis.cxx +++ b/source/framework/analysis/src/TRestPulseShapeAnalysis.cxx @@ -627,8 +627,10 @@ Double_t TRestPulseShapeAnalysis::GetSlopeIntegral(const std::vector >& signal); -template Double_t TRestPulseShapeAnalysis::GetSlopeIntegral(const std::vector >& signal); +template Double_t TRestPulseShapeAnalysis::GetSlopeIntegral( + const std::vector >& signal); +template Double_t TRestPulseShapeAnalysis::GetSlopeIntegral( + const std::vector >& signal); /////////////////////////////////////////////// /// \brief It returns the slope of the first positive @@ -650,9 +652,10 @@ Double_t TRestPulseShapeAnalysis::GetRiseSlope(const std::vector >& signal); -template Double_t TRestPulseShapeAnalysis::GetRiseSlope(const std::vector >& signal); - +template Double_t TRestPulseShapeAnalysis::GetRiseSlope( + const std::vector >& signal); +template Double_t TRestPulseShapeAnalysis::GetRiseSlope( + const std::vector >& signal); /////////////////////////////////////////////// /// \brief It returns the time of the first positive @@ -673,5 +676,5 @@ Double_t TRestPulseShapeAnalysis::GetRiseTime(const std::vector >& signal); -template Double_t TRestPulseShapeAnalysis::GetRiseTime(const std::vector >& signal); - +template Double_t TRestPulseShapeAnalysis::GetRiseTime( + const std::vector >& signal);