From d1162ebaa2e34a3cf6a214d50ccf3cf5e6794fbb Mon Sep 17 00:00:00 2001 From: JPorron Date: Wed, 6 Nov 2024 14:24:40 +0100 Subject: [PATCH 01/11] Use average BaseLine and BaseLineSigma and small changes in peak position --- src/TRestRawPeaksFinderProcess.cxx | 40 ++++++++++++++++++++++++++- src/TRestRawSignal.cxx | 43 +++++++++++++++--------------- 2 files changed, 61 insertions(+), 22 deletions(-) diff --git a/src/TRestRawPeaksFinderProcess.cxx b/src/TRestRawPeaksFinderProcess.cxx index 112fe1e..2bdde60 100644 --- a/src/TRestRawPeaksFinderProcess.cxx +++ b/src/TRestRawPeaksFinderProcess.cxx @@ -30,6 +30,43 @@ TRestEvent* TRestRawPeaksFinderProcess::ProcessEvent(TRestEvent* inputEvent) { vector> eventPeaks; // signalId, time, amplitude + // Calculate average baseline and sigma of all the TPC signals + double BaseLineMean = 0.0; + double BaseLineSigmaMean = 0.0; + int count = 0; + + for (int signalIndex = 0; signalIndex < fInputEvent->GetNumberOfSignals(); signalIndex++) { + const auto signal = fInputEvent->GetSignal(signalIndex); + const UShort_t signalId = signal->GetSignalID(); + + const string channelType = fReadoutMetadata->GetTypeForChannelDaqId(signalId); + const string channelName = fReadoutMetadata->GetNameForChannelDaqId(signalId); + + // Check if channel type is in the list of selected channel types + if (fChannelTypes.find(channelType) == fChannelTypes.end()) { + continue; + } + + // Choose appropriate function based on channel type + if (channelType == "tpc") { + signal->CalculateBaseLine(fBaselineRange.X(), fBaselineRange.Y()); + + double baseline = signal->GetBaseLine(); + double baselinesigma = signal->GetBaseLineSigma(); + + // Accumulate the baseline and sigma values + BaseLineMean += baseline; + BaseLineSigmaMean += baselinesigma; + count++; // Count the signals considered + } + } + + // Calculate the average if there were any matching signals + if (count > 0) { + BaseLineMean /= count; + BaseLineSigmaMean /= count; + } + for (int signalIndex = 0; signalIndex < fInputEvent->GetNumberOfSignals(); signalIndex++) { const auto signal = fInputEvent->GetSignal(signalIndex); const UShort_t signalId = signal->GetSignalID(); @@ -45,8 +82,9 @@ TRestEvent* TRestRawPeaksFinderProcess::ProcessEvent(TRestEvent* inputEvent) { // Choose appropriate function based on channel type if (channelType == "tpc") { signal->CalculateBaseLine(fBaselineRange.X(), fBaselineRange.Y()); + const auto peaks = - signal->GetPeaks(signal->GetBaseLine() + 5 * signal->GetBaseLineSigma(), fDistance); + signal->GetPeaks(BaseLineMean + 10 * BaseLineSigmaMean, fDistance); for (const auto& [time, amplitude] : peaks) { eventPeaks.emplace_back(signalId, time, amplitude); diff --git a/src/TRestRawSignal.cxx b/src/TRestRawSignal.cxx index bf2c4a6..fbb154f 100644 --- a/src/TRestRawSignal.cxx +++ b/src/TRestRawSignal.cxx @@ -977,30 +977,31 @@ vector> TRestRawSignal::GetPeaks(double threshold, UShort } // If it's a peak and it´s above the threshold and further than distance to the previous peak, add - // to peaks + // to peaks the biggest amplitude bin within the next "distance" bins and as amplitude the TripleMaxAverage. + // This is because for flat regions the detected peak is more to the left than the actual one. if (isPeak && smoothedValue > threshold) { if (peaks.empty() || i - peaks.back().first >= distance) { - double fitMinRange = i - 20; - double fitMaxRange = i + 20; - - // Create a Gaussian fit function - TF1 fitFunction("gaussianFit", "gaus", fitMinRange, fitMaxRange); - // Fit the data with the Gaussian function - fitFunction.SetRange(fitMinRange, fitMaxRange); // Initial parameters - - // Create histogram with the values to fit - TH1D histogram("hist", "hist", 40, fitMinRange, fitMaxRange); - for (int k = i - 20; k <= i + 20; ++k) { - histogram.SetBinContent(k - (i - 20) + 1, GetRawData(k)); // Set bin content + + // Initialize variables to find the max amplitude within the next "distance" bins + int maxBin = i; + double maxAmplitude = smoothedValues[i]; + + // Look ahead within the specified distance to find the bin with the maximum amplitude + for (int j = i + 1; j <= i + distance && j < smoothedValues.size(); ++j) { + if (smoothedValues[j] > maxAmplitude) { + maxAmplitude = smoothedValues[j]; + maxBin = j; + } } - histogram.Fit(&fitFunction, "RQ"); - - // Get peak position and amplitude from the fit - double peakPosition = fitFunction.GetParameter(1); - UShort_t formattedPeakPosition = static_cast(peakPosition); - double peakAmplitude = GetRawData(formattedPeakPosition); - - peaks.push_back(std::make_pair(formattedPeakPosition, peakAmplitude)); + + // Calculate the peak amplitude as the average of maxBin and its two neighbors + double amplitude1 = GetRawData(maxBin - 1); + double amplitude2 = GetRawData(maxBin); + double amplitude3 = GetRawData(maxBin + 1); + double peakAmplitude = (amplitude1 + amplitude2 + amplitude3) / 3.0; + + // Store the peak position and amplitude + peaks.push_back(std::make_pair(maxBin, peakAmplitude)); } } } From cff65a935160c584f419f780dbab14b347333038 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Wed, 6 Nov 2024 13:40:03 +0000 Subject: [PATCH 02/11] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/TRestRawPeaksFinderProcess.cxx | 5 ++--- src/TRestRawSignal.cxx | 12 ++++++------ 2 files changed, 8 insertions(+), 9 deletions(-) diff --git a/src/TRestRawPeaksFinderProcess.cxx b/src/TRestRawPeaksFinderProcess.cxx index 2bdde60..da5540d 100644 --- a/src/TRestRawPeaksFinderProcess.cxx +++ b/src/TRestRawPeaksFinderProcess.cxx @@ -50,7 +50,7 @@ TRestEvent* TRestRawPeaksFinderProcess::ProcessEvent(TRestEvent* inputEvent) { // Choose appropriate function based on channel type if (channelType == "tpc") { signal->CalculateBaseLine(fBaselineRange.X(), fBaselineRange.Y()); - + double baseline = signal->GetBaseLine(); double baselinesigma = signal->GetBaseLineSigma(); @@ -83,8 +83,7 @@ TRestEvent* TRestRawPeaksFinderProcess::ProcessEvent(TRestEvent* inputEvent) { if (channelType == "tpc") { signal->CalculateBaseLine(fBaselineRange.X(), fBaselineRange.Y()); - const auto peaks = - signal->GetPeaks(BaseLineMean + 10 * BaseLineSigmaMean, fDistance); + const auto peaks = signal->GetPeaks(BaseLineMean + 10 * BaseLineSigmaMean, fDistance); for (const auto& [time, amplitude] : peaks) { eventPeaks.emplace_back(signalId, time, amplitude); diff --git a/src/TRestRawSignal.cxx b/src/TRestRawSignal.cxx index fbb154f..257a6d1 100644 --- a/src/TRestRawSignal.cxx +++ b/src/TRestRawSignal.cxx @@ -977,15 +977,15 @@ vector> TRestRawSignal::GetPeaks(double threshold, UShort } // If it's a peak and it´s above the threshold and further than distance to the previous peak, add - // to peaks the biggest amplitude bin within the next "distance" bins and as amplitude the TripleMaxAverage. - // This is because for flat regions the detected peak is more to the left than the actual one. + // to peaks the biggest amplitude bin within the next "distance" bins and as amplitude the + // TripleMaxAverage. This is because for flat regions the detected peak is more to the left than + // the actual one. if (isPeak && smoothedValue > threshold) { if (peaks.empty() || i - peaks.back().first >= distance) { - // Initialize variables to find the max amplitude within the next "distance" bins int maxBin = i; double maxAmplitude = smoothedValues[i]; - + // Look ahead within the specified distance to find the bin with the maximum amplitude for (int j = i + 1; j <= i + distance && j < smoothedValues.size(); ++j) { if (smoothedValues[j] > maxAmplitude) { @@ -993,13 +993,13 @@ vector> TRestRawSignal::GetPeaks(double threshold, UShort maxBin = j; } } - + // Calculate the peak amplitude as the average of maxBin and its two neighbors double amplitude1 = GetRawData(maxBin - 1); double amplitude2 = GetRawData(maxBin); double amplitude3 = GetRawData(maxBin + 1); double peakAmplitude = (amplitude1 + amplitude2 + amplitude3) / 3.0; - + // Store the peak position and amplitude peaks.push_back(std::make_pair(maxBin, peakAmplitude)); } From a31c59f2d4ed4b44585df5156a05515f290ffcd6 Mon Sep 17 00:00:00 2001 From: JPorron Date: Thu, 7 Nov 2024 12:18:13 +0100 Subject: [PATCH 03/11] Solve warning that is being treated as error --- src/TRestRawSignal.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/TRestRawSignal.cxx b/src/TRestRawSignal.cxx index 257a6d1..4d482ea 100644 --- a/src/TRestRawSignal.cxx +++ b/src/TRestRawSignal.cxx @@ -987,7 +987,7 @@ vector> TRestRawSignal::GetPeaks(double threshold, UShort double maxAmplitude = smoothedValues[i]; // Look ahead within the specified distance to find the bin with the maximum amplitude - for (int j = i + 1; j <= i + distance && j < smoothedValues.size(); ++j) { + for (std::vector::size_type j = i + 1; j <= i + distance && j < smoothedValues.size(); ++j) { if (smoothedValues[j] > maxAmplitude) { maxAmplitude = smoothedValues[j]; maxBin = j; From 75810dbe9f52999167ace790026c79038d43a895 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 8 Nov 2024 08:30:29 +0000 Subject: [PATCH 04/11] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/TRestRawSignal.cxx | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/TRestRawSignal.cxx b/src/TRestRawSignal.cxx index 4d482ea..f4020dd 100644 --- a/src/TRestRawSignal.cxx +++ b/src/TRestRawSignal.cxx @@ -987,7 +987,8 @@ vector> TRestRawSignal::GetPeaks(double threshold, UShort double maxAmplitude = smoothedValues[i]; // Look ahead within the specified distance to find the bin with the maximum amplitude - for (std::vector::size_type j = i + 1; j <= i + distance && j < smoothedValues.size(); ++j) { + for (std::vector::size_type j = i + 1; + j <= i + distance && j < smoothedValues.size(); ++j) { if (smoothedValues[j] > maxAmplitude) { maxAmplitude = smoothedValues[j]; maxBin = j; From 8825eae8f4d27fd31ec3c26300859120fc114720 Mon Sep 17 00:00:00 2001 From: Luis Antonio Obis Aparicio Date: Fri, 8 Nov 2024 10:28:11 +0100 Subject: [PATCH 05/11] add check just in case --- src/TRestRawPeaksFinderProcess.cxx | 31 ++++++++++++++++++------------ 1 file changed, 19 insertions(+), 12 deletions(-) diff --git a/src/TRestRawPeaksFinderProcess.cxx b/src/TRestRawPeaksFinderProcess.cxx index da5540d..aa3019e 100644 --- a/src/TRestRawPeaksFinderProcess.cxx +++ b/src/TRestRawPeaksFinderProcess.cxx @@ -33,7 +33,7 @@ TRestEvent* TRestRawPeaksFinderProcess::ProcessEvent(TRestEvent* inputEvent) { // Calculate average baseline and sigma of all the TPC signals double BaseLineMean = 0.0; double BaseLineSigmaMean = 0.0; - int count = 0; + unsigned int countTPC = 0; for (int signalIndex = 0; signalIndex < fInputEvent->GetNumberOfSignals(); signalIndex++) { const auto signal = fInputEvent->GetSignal(signalIndex); @@ -50,21 +50,17 @@ TRestEvent* TRestRawPeaksFinderProcess::ProcessEvent(TRestEvent* inputEvent) { // Choose appropriate function based on channel type if (channelType == "tpc") { signal->CalculateBaseLine(fBaselineRange.X(), fBaselineRange.Y()); - - double baseline = signal->GetBaseLine(); - double baselinesigma = signal->GetBaseLineSigma(); - // Accumulate the baseline and sigma values - BaseLineMean += baseline; - BaseLineSigmaMean += baselinesigma; - count++; // Count the signals considered + BaseLineMean += signal->GetBaseLine(); + BaseLineSigmaMean += signal->GetBaseLineSigma(); + countTPC++; // Count the signals considered } } // Calculate the average if there were any matching signals - if (count > 0) { - BaseLineMean /= count; - BaseLineSigmaMean /= count; + if (countTPC > 0) { + BaseLineMean /= countTPC; + BaseLineSigmaMean /= countTPC; } for (int signalIndex = 0; signalIndex < fInputEvent->GetNumberOfSignals(); signalIndex++) { @@ -83,7 +79,18 @@ TRestEvent* TRestRawPeaksFinderProcess::ProcessEvent(TRestEvent* inputEvent) { if (channelType == "tpc") { signal->CalculateBaseLine(fBaselineRange.X(), fBaselineRange.Y()); - const auto peaks = signal->GetPeaks(BaseLineMean + 10 * BaseLineSigmaMean, fDistance); + constexpr double numberOfBaselinesThreshold = 10; + // I think count will never be 0, just in case + const double threshold = + (countTPC > 0) ? BaseLineMean + numberOfBaselinesThreshold * BaseLineSigmaMean + : signal->GetBaseLine() + numberOfBaselinesThreshold * signal->GetBaseLineSigma(); + if (countTPC <= 0) { + cerr << "TRestRawPeaksFinderProcess::ProcessEvent: TPC count is 0 in TPC loop, this should " + "not happen" + << endl; + exit(1); + } + const auto peaks = signal->GetPeaks(threshold, fDistance); for (const auto& [time, amplitude] : peaks) { eventPeaks.emplace_back(signalId, time, amplitude); From 6845cc2731ef7df517c1dc5dfdf454e352cb4ec7 Mon Sep 17 00:00:00 2001 From: Luis Antonio Obis Aparicio Date: Fri, 8 Nov 2024 10:29:37 +0100 Subject: [PATCH 06/11] typo in rml variable name --- src/TRestRawPeaksFinderProcess.cxx | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/TRestRawPeaksFinderProcess.cxx b/src/TRestRawPeaksFinderProcess.cxx index aa3019e..0c4139a 100644 --- a/src/TRestRawPeaksFinderProcess.cxx +++ b/src/TRestRawPeaksFinderProcess.cxx @@ -369,8 +369,8 @@ void TRestRawPeaksFinderProcess::InitFromConfigFile() { fBaselineRange = Get2DVectorParameterWithUnits("baselineRange", fBaselineRange); fDistance = StringToDouble(GetParameter("distance", fDistance)); fWindow = StringToDouble(GetParameter("window", fWindow)); - fRemoveAllVetoes = StringToBool(GetParameter("removeAllVetos", fRemoveAllVetoes)); - fRemovePeaklessVetoes = StringToBool(GetParameter("removePeaklessVetos", fRemovePeaklessVetoes)); + fRemoveAllVetoes = StringToBool(GetParameter("removeAllVetoes", fRemoveAllVetoes)); + fRemovePeaklessVetoes = StringToBool(GetParameter("removePeaklessVetoes", fRemovePeaklessVetoes)); fTimeBinToTimeFactorMultiplier = GetDblParameterWithUnits("sampling", fTimeBinToTimeFactorMultiplier); fTimeBinToTimeFactorOffset = GetDblParameterWithUnits("trigDelay", fTimeBinToTimeFactorOffset); @@ -415,7 +415,7 @@ void TRestRawPeaksFinderProcess::InitFromConfigFile() { if (filterType != "veto" && fRemovePeaklessVetoes) { cerr << "TRestRawPeaksFinderProcess::InitProcess: removing veto signals only makes sense when the " - "process is applied to veto signals. Remove \"removePeaklessVetos\" parameter" + "process is applied to veto signals. Remove \"removePeaklessVetoes\" parameter" << endl; exit(1); } From afea755727697232ec2088cc60a311c5e1298c84 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 8 Nov 2024 09:32:28 +0000 Subject: [PATCH 07/11] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/TRestRawPeaksFinderProcess.cxx | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/TRestRawPeaksFinderProcess.cxx b/src/TRestRawPeaksFinderProcess.cxx index 0c4139a..7036b32 100644 --- a/src/TRestRawPeaksFinderProcess.cxx +++ b/src/TRestRawPeaksFinderProcess.cxx @@ -82,8 +82,9 @@ TRestEvent* TRestRawPeaksFinderProcess::ProcessEvent(TRestEvent* inputEvent) { constexpr double numberOfBaselinesThreshold = 10; // I think count will never be 0, just in case const double threshold = - (countTPC > 0) ? BaseLineMean + numberOfBaselinesThreshold * BaseLineSigmaMean - : signal->GetBaseLine() + numberOfBaselinesThreshold * signal->GetBaseLineSigma(); + (countTPC > 0) + ? BaseLineMean + numberOfBaselinesThreshold * BaseLineSigmaMean + : signal->GetBaseLine() + numberOfBaselinesThreshold * signal->GetBaseLineSigma(); if (countTPC <= 0) { cerr << "TRestRawPeaksFinderProcess::ProcessEvent: TPC count is 0 in TPC loop, this should " "not happen" From 83ea59864c6b10e68261c6543e6e2ec09e5e4233 Mon Sep 17 00:00:00 2001 From: Luis Antonio Obis Aparicio Date: Fri, 8 Nov 2024 10:35:54 +0100 Subject: [PATCH 08/11] emplace_back --- src/TRestRawSignal.cxx | 58 +++++++++++++++++++++++++++++------------- 1 file changed, 40 insertions(+), 18 deletions(-) diff --git a/src/TRestRawSignal.cxx b/src/TRestRawSignal.cxx index f4020dd..72747e3 100644 --- a/src/TRestRawSignal.cxx +++ b/src/TRestRawSignal.cxx @@ -264,12 +264,15 @@ void TRestRawSignal::InitializePointsOverThreshold(const TVector2& thrPar, Int_t if (pulse.size() >= (unsigned int)nPointsOver) { // auto stdev = TMath::StdDev(begin(pulse), end(pulse)); // calculate stdev - double mean = std::accumulate(pulse.begin(), pulse.end(), 0.0) / pulse.size(); + double mean = std::accumulate(pulse.begin(), pulse.end(), 0.0) / double(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); + double stdev = std::sqrt(sq_sum / double(pulse.size()) - mean * mean); - if (stdev > signalTh * fBaseLineSigma) - for (int j = pos; j < i; j++) fPointsOverThreshold.push_back(j); + if (stdev > signalTh * fBaseLineSigma) { + for (int j = pos; j < i; j++) { + fPointsOverThreshold.push_back(j); + } + } } } } @@ -283,14 +286,18 @@ void TRestRawSignal::InitializePointsOverThreshold(const TVector2& thrPar, Int_t /// of fThresholdIntegral. This method is only used internally. /// void TRestRawSignal::CalculateThresholdIntegral() { - if (fRange.X() < 0) fRange.SetX(0); - if (fRange.Y() <= 0 || fRange.Y() > GetNumberOfPoints()) fRange.SetY(GetNumberOfPoints()); + if (fRange.X() < 0) { + fRange.SetX(0); + } + if (fRange.Y() <= 0 || fRange.Y() > GetNumberOfPoints()) { + fRange.SetY(GetNumberOfPoints()); + } fThresholdIntegral = 0; - for (unsigned int n = 0; n < fPointsOverThreshold.size(); n++) { - if (fPointsOverThreshold[n] >= fRange.X() && fPointsOverThreshold[n] < fRange.Y()) { - fThresholdIntegral += GetData(fPointsOverThreshold[n]); + for (int n : fPointsOverThreshold) { + if (n >= fRange.X() && n < fRange.Y()) { + fThresholdIntegral += GetData(n); } } } @@ -301,11 +308,17 @@ void TRestRawSignal::CalculateThresholdIntegral() { /// the integral is calculated in the full range. /// Double_t TRestRawSignal::GetIntegral() { - if (fRange.X() < 0) fRange.SetX(0); - if (fRange.Y() <= 0 || fRange.Y() > GetNumberOfPoints()) fRange.SetY(GetNumberOfPoints()); + if (fRange.X() < 0) { + fRange.SetX(0); + } + if (fRange.Y() <= 0 || fRange.Y() > GetNumberOfPoints()) { + fRange.SetY(GetNumberOfPoints()); + } Double_t sum = 0; - for (int i = fRange.X(); i < fRange.Y(); i++) sum += GetData(i); + for (int i = fRange.X(); i < fRange.Y(); i++) { + sum += GetData(i); + } return sum; } @@ -314,11 +327,17 @@ Double_t TRestRawSignal::GetIntegral() { /// by (startBin,endBin). /// Double_t TRestRawSignal::GetIntegralInRange(Int_t startBin, Int_t endBin) { - if (startBin < 0) startBin = 0; - if (endBin <= 0 || endBin > GetNumberOfPoints()) endBin = GetNumberOfPoints(); + if (startBin < 0) { + startBin = 0; + } + if (endBin <= 0 || endBin > GetNumberOfPoints()) { + endBin = GetNumberOfPoints(); + } Double_t sum = 0; - for (int i = startBin; i < endBin; i++) sum += GetRawData(i); + for (int i = startBin; i < endBin; i++) { + sum += GetRawData(i); + } return sum; } @@ -328,7 +347,7 @@ Double_t TRestRawSignal::GetIntegralInRange(Int_t startBin, Int_t endBin) { /// have been called first. /// Double_t TRestRawSignal::GetThresholdIntegral() { - if (fThresholdIntegral == -1) + if (fThresholdIntegral == -1) { if (fShowWarnings) { std::cout << "TRestRawSignal::GetThresholdIntegral. " "InitializePointsOverThreshold should be " @@ -336,6 +355,7 @@ Double_t TRestRawSignal::GetThresholdIntegral() { << endl; fShowWarnings = false; } + } return fThresholdIntegral; } @@ -917,7 +937,9 @@ vector> TRestRawSignal::GetPeaks(double threshold, UShort 10; // Region to compare for peak/no peak classification. 10 means 5 bins to each side const size_t numPoints = GetNumberOfPoints(); - if (numPoints == 0) return peaks; + if (numPoints == 0) { + return peaks; + } // Pre-calculate smoothed values for all bins using a rolling sum vector smoothedValues(numPoints, 0.0); @@ -1002,7 +1024,7 @@ vector> TRestRawSignal::GetPeaks(double threshold, UShort double peakAmplitude = (amplitude1 + amplitude2 + amplitude3) / 3.0; // Store the peak position and amplitude - peaks.push_back(std::make_pair(maxBin, peakAmplitude)); + peaks.emplace_back(maxBin, peakAmplitude); } } } From 67850507a5f339c794e4a2b848aab7823f453b89 Mon Sep 17 00:00:00 2001 From: JPorron Date: Mon, 11 Nov 2024 10:27:10 +0100 Subject: [PATCH 09/11] Set threshold variable of tpc as input --- inc/TRestRawPeaksFinderProcess.h | 4 +++- src/TRestRawPeaksFinderProcess.cxx | 13 ++++++++++--- 2 files changed, 13 insertions(+), 4 deletions(-) diff --git a/inc/TRestRawPeaksFinderProcess.h b/inc/TRestRawPeaksFinderProcess.h index 85693f9..7e99b12 100644 --- a/inc/TRestRawPeaksFinderProcess.h +++ b/inc/TRestRawPeaksFinderProcess.h @@ -17,6 +17,8 @@ class TRestRawPeaksFinderProcess : public TRestEventProcess { /// \brief threshold over baseline to consider a peak Double_t fThresholdOverBaseline = 2.0; + /// \brief choose times the sigma of the baseline must be overcome to consider a peak + Double_t fSigmaOverBaseline = 10.0; /// \brief range of samples to calculate baseline for peak finding TVector2 fBaselineRange = {0, 10}; /// \brief distance between two peaks to consider them as different (ADC units) @@ -58,7 +60,7 @@ class TRestRawPeaksFinderProcess : public TRestEventProcess { TRestRawPeaksFinderProcess() = default; ~TRestRawPeaksFinderProcess() = default; - ClassDefOverride(TRestRawPeaksFinderProcess, 5); + ClassDefOverride(TRestRawPeaksFinderProcess, 6); }; #endif // REST_TRESTRAWPEAKSFINDERPROCESS_H diff --git a/src/TRestRawPeaksFinderProcess.cxx b/src/TRestRawPeaksFinderProcess.cxx index 7036b32..0c823c4 100644 --- a/src/TRestRawPeaksFinderProcess.cxx +++ b/src/TRestRawPeaksFinderProcess.cxx @@ -79,12 +79,11 @@ TRestEvent* TRestRawPeaksFinderProcess::ProcessEvent(TRestEvent* inputEvent) { if (channelType == "tpc") { signal->CalculateBaseLine(fBaselineRange.X(), fBaselineRange.Y()); - constexpr double numberOfBaselinesThreshold = 10; // I think count will never be 0, just in case const double threshold = (countTPC > 0) - ? BaseLineMean + numberOfBaselinesThreshold * BaseLineSigmaMean - : signal->GetBaseLine() + numberOfBaselinesThreshold * signal->GetBaseLineSigma(); + ? BaseLineMean + fSigmaOverBaseline * BaseLineSigmaMean + : signal->GetBaseLine() + fSigmaOverBaseline * signal->GetBaseLineSigma(); if (countTPC <= 0) { cerr << "TRestRawPeaksFinderProcess::ProcessEvent: TPC count is 0 in TPC loop, this should " "not happen" @@ -367,6 +366,7 @@ void TRestRawPeaksFinderProcess::InitFromConfigFile() { } fThresholdOverBaseline = StringToDouble(GetParameter("thresholdOverBaseline", fThresholdOverBaseline)); + fSigmaOverBaseline = StringToDouble(GetParameter("sigmaOverBaseline", fSigmaOverBaseline)); fBaselineRange = Get2DVectorParameterWithUnits("baselineRange", fBaselineRange); fDistance = StringToDouble(GetParameter("distance", fDistance)); fWindow = StringToDouble(GetParameter("window", fWindow)); @@ -404,6 +404,11 @@ void TRestRawPeaksFinderProcess::InitFromConfigFile() { exit(1); } + if (fSigmaOverBaseline < 0) { + cerr << "TRestRawPeaksFinderProcess::InitProcess: sigma over baseline is < 0" << endl; + exit(1); + } + if (fDistance <= 0) { cerr << "TRestRawPeaksFinderProcess::InitProcess: distance is < 0" << endl; exit(1); @@ -433,6 +438,8 @@ void TRestRawPeaksFinderProcess::PrintMetadata() { RESTMetadata << "Baseline range for tpc signals: " << fBaselineRange.X() << " - " << fBaselineRange.Y() << RESTendl; + RESTMetadata << "Sigma over baseline for tpc signals: " << fSigmaOverBaseline << RESTendl; + RESTMetadata << "Threshold over baseline for veto signals: " << fThresholdOverBaseline << RESTendl; RESTMetadata << "Distance: " << fDistance << RESTendl; From 8a72dde0f34c406c92e34e6f95dc8bfbd852a6d3 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 11 Nov 2024 09:28:39 +0000 Subject: [PATCH 10/11] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/TRestRawPeaksFinderProcess.cxx | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/TRestRawPeaksFinderProcess.cxx b/src/TRestRawPeaksFinderProcess.cxx index 0c823c4..8bb3038 100644 --- a/src/TRestRawPeaksFinderProcess.cxx +++ b/src/TRestRawPeaksFinderProcess.cxx @@ -81,9 +81,8 @@ TRestEvent* TRestRawPeaksFinderProcess::ProcessEvent(TRestEvent* inputEvent) { // I think count will never be 0, just in case const double threshold = - (countTPC > 0) - ? BaseLineMean + fSigmaOverBaseline * BaseLineSigmaMean - : signal->GetBaseLine() + fSigmaOverBaseline * signal->GetBaseLineSigma(); + (countTPC > 0) ? BaseLineMean + fSigmaOverBaseline * BaseLineSigmaMean + : signal->GetBaseLine() + fSigmaOverBaseline * signal->GetBaseLineSigma(); if (countTPC <= 0) { cerr << "TRestRawPeaksFinderProcess::ProcessEvent: TPC count is 0 in TPC loop, this should " "not happen" From d53cec6cfd58092ecea50924b84f76d360b02b56 Mon Sep 17 00:00:00 2001 From: Luis Antonio Obis Aparicio Date: Mon, 11 Nov 2024 10:42:10 +0100 Subject: [PATCH 11/11] bad initialization defaults --- src/TRestRawPeaksFinderProcess.cxx | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/TRestRawPeaksFinderProcess.cxx b/src/TRestRawPeaksFinderProcess.cxx index 8bb3038..92c2107 100644 --- a/src/TRestRawPeaksFinderProcess.cxx +++ b/src/TRestRawPeaksFinderProcess.cxx @@ -364,11 +364,11 @@ void TRestRawPeaksFinderProcess::InitFromConfigFile() { // if no channel type is specified, use all channel types } - fThresholdOverBaseline = StringToDouble(GetParameter("thresholdOverBaseline", fThresholdOverBaseline)); - fSigmaOverBaseline = StringToDouble(GetParameter("sigmaOverBaseline", fSigmaOverBaseline)); + fThresholdOverBaseline = GetDblParameterWithUnits("thresholdOverBaseline", fThresholdOverBaseline); + fSigmaOverBaseline = GetDblParameterWithUnits("sigmaOverBaseline", fSigmaOverBaseline); fBaselineRange = Get2DVectorParameterWithUnits("baselineRange", fBaselineRange); - fDistance = StringToDouble(GetParameter("distance", fDistance)); - fWindow = StringToDouble(GetParameter("window", fWindow)); + fDistance = UShort_t(GetDblParameterWithUnits("distance", fDistance)); + fWindow = UShort_t(GetDblParameterWithUnits("window", fWindow)); fRemoveAllVetoes = StringToBool(GetParameter("removeAllVetoes", fRemoveAllVetoes)); fRemovePeaklessVetoes = StringToBool(GetParameter("removePeaklessVetoes", fRemovePeaklessVetoes));