Use average BaseLine and BaseLineSigma and small changes in peak pos

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

src/TRestRawPeaksFinderProcess.cxx

@@ -30,6 +30,39 @@ TRestEvent* TRestRawPeaksFinderProcess::ProcessEvent(TRestEvent* inputEvent) {
 
     vector<tuple<UShort_t, UShort_t, double>> eventPeaks;  // signalId, time, amplitude
 
+    // Calculate average baseline and sigma of all the TPC signals
+    double BaseLineMean = 0.0;
+    double BaseLineSigmaMean = 0.0;
+    unsigned int countTPC = 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());
+            // Accumulate the baseline and sigma values
+            BaseLineMean += signal->GetBaseLine();
+            BaseLineSigmaMean += signal->GetBaseLineSigma();
+            countTPC++;  // Count the signals considered
+        }
+    }
+
+    // Calculate the average if there were any matching signals
+    if (countTPC > 0) {
+        BaseLineMean /= countTPC;
+        BaseLineSigmaMean /= countTPC;
+    }
+
     for (int signalIndex = 0; signalIndex < fInputEvent->GetNumberOfSignals(); signalIndex++) {
         const auto signal = fInputEvent->GetSignal(signalIndex);
         const UShort_t signalId = signal->GetSignalID();
@@ -45,8 +78,18 @@ 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);
+
+            // I think count will never be 0, just in case
+            const double threshold =
+                (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"
+                     << endl;
+                exit(1);
+            }
+            const auto peaks = signal->GetPeaks(threshold, fDistance);
 
             for (const auto& [time, amplitude] : peaks) {
                 eventPeaks.emplace_back(signalId, time, amplitude);
@@ -321,12 +364,13 @@ void TRestRawPeaksFinderProcess::InitFromConfigFile() {
         // if no channel type is specified, use all channel types
     }
 
-    fThresholdOverBaseline = StringToDouble(GetParameter("thresholdOverBaseline", fThresholdOverBaseline));
+    fThresholdOverBaseline = GetDblParameterWithUnits("thresholdOverBaseline", fThresholdOverBaseline);
+    fSigmaOverBaseline = GetDblParameterWithUnits("sigmaOverBaseline", fSigmaOverBaseline);
     fBaselineRange = Get2DVectorParameterWithUnits("baselineRange", fBaselineRange);
-    fDistance = StringToDouble(GetParameter("distance", fDistance));
-    fWindow = StringToDouble(GetParameter("window", fWindow));
-    fRemoveAllVetoes = StringToBool(GetParameter("removeAllVetos", fRemoveAllVetoes));
-    fRemovePeaklessVetoes = StringToBool(GetParameter("removePeaklessVetos", fRemovePeaklessVetoes));
+    fDistance = UShort_t(GetDblParameterWithUnits("distance", fDistance));
+    fWindow = UShort_t(GetDblParameterWithUnits("window", fWindow));
+    fRemoveAllVetoes = StringToBool(GetParameter("removeAllVetoes", fRemoveAllVetoes));
+    fRemovePeaklessVetoes = StringToBool(GetParameter("removePeaklessVetoes", fRemovePeaklessVetoes));
 
     fTimeBinToTimeFactorMultiplier = GetDblParameterWithUnits("sampling", fTimeBinToTimeFactorMultiplier);
     fTimeBinToTimeFactorOffset = GetDblParameterWithUnits("trigDelay", fTimeBinToTimeFactorOffset);
@@ -359,6 +403,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);
@@ -371,7 +420,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);
     }
@@ -388,6 +437,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;

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,14 +347,15 @@ 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 "
                          "called first!"
                       << endl;
             fShowWarnings = false;
         }
+    }
     return fThresholdIntegral;
 }
 
@@ -917,7 +937,9 @@ vector<pair<UShort_t, double>> 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<double> smoothedValues(numPoints, 0.0);
@@ -977,30 +999,32 @@ vector<pair<UShort_t, double>> 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 (std::vector<double>::size_type 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<UShort_t>(peakPosition);
-                    double peakAmplitude = GetRawData(formattedPeakPosition);
+                    // 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;
 
-                    peaks.push_back(std::make_pair(formattedPeakPosition, peakAmplitude));
+                    // Store the peak position and amplitude
+                    peaks.emplace_back(maxBin, peakAmplitude);
                 }
             }
         }