Calculate veto signals baseline without outliers

inc/TRestRawSignal.h

@@ -41,6 +41,8 @@ class TRestRawSignal {
 
     void CalculateBaseLineSigmaIQR(Int_t startBin, Int_t endBin);
 
+    void CalculateBaseLineSigmaExcludeOutliers(Int_t startBin, Int_t endBin);
+
     std::vector<Float_t> GetSignalSmoothed_ExcludeOutliers(Int_t averagingPoints);
 
    protected:
@@ -198,6 +200,8 @@ class TRestRawSignal {
 
     void CalculateBaseLineMedian(Int_t startBin, Int_t endBin);
 
+    void CalculateBaseLineMedianExcludeOutliers(Int_t startBin, Int_t endBin);
+
     void CalculateBaseLine(Int_t startBin, Int_t endBin, const std::string& option = "");
 
     void GetBaseLineCorrected(TRestRawSignal* smoothedSignal, Int_t averagingPoints);

src/TRestRawPeaksFinderProcess.cxx

@@ -97,7 +97,7 @@ TRestEvent* TRestRawPeaksFinderProcess::ProcessEvent(TRestEvent* inputEvent) {
         } else if (channelType == "veto") {
             // For veto signals the baseline is calculated over the whole range, as we don´t know where the
             // signal will be.
-            signal->CalculateBaseLine(0, 511, "robust");
+            signal->CalculateBaseLine(0, 511, "OUTLIERS");
             // For veto signals the threshold is selected by the user.
             const auto peaks =
                 signal->GetPeaksVeto(signal->GetBaseLine() + fThresholdOverBaseline, fDistance);

src/TRestRawSignal.cxx

@@ -784,6 +784,50 @@ void TRestRawSignal::CalculateBaseLineMedian(Int_t startBin, Int_t endBin) {
     }
 }
 
+///////////////////////////////////////////////
+/// \brief This method is called by CalculateBaseLine with the "OUTLIERS"-option and is used to determine the
+/// value of the baseline as the median of the data points found in the range defined between startBin and
+/// endBin after excluding the outliers. The median is calculated using only the values in the 25-75% removing
+/// big and small outliers.
+///
+void TRestRawSignal::CalculateBaseLineMedianExcludeOutliers(Int_t startBin, Int_t endBin) {
+    if (endBin - startBin <= 0) {
+        fBaseLine = 0.;
+        return;
+    } else if (endBin > static_cast<int>(fSignalData.size())) {
+        cout << "TRestRawSignal::CalculateBaseLine. Error! Baseline range exceeds the rawdata depth!!"
+             << endl;
+        endBin = fSignalData.size();
+    } else {
+        // Extract the data within the interval
+        std::vector<Short_t> data(fSignalData.begin() + startBin, fSignalData.begin() + endBin);
+        std::sort(data.begin(), data.end());
+
+        // Calculate Q1 and Q3 for IQR
+        size_t dataSize = data.size();
+        Short_t Q1 = data[dataSize / 4];
+        Short_t Q3 = data[(3 * dataSize) / 4];
+        Double_t lowerBound = Q1;
+        Double_t upperBound = Q3;
+
+        // Filter out the outliers
+        std::vector<Short_t> filteredData;
+        for (const auto& value : data) {
+            if (value >= lowerBound && value <= upperBound) {
+                filteredData.emplace_back(value);
+            }
+        }
+
+        // Calculate median of filtered data
+        if (filteredData.empty()) {
+            fBaseLine = TMath::Median(data.size(),
+                                      &data[0]);  // Fall back to original median if all values are outliers
+        } else {
+            fBaseLine = TMath::Median(filteredData.size(), &filteredData[0]);
+        }
+    }
+}
+
 ///////////////////////////////////////////////
 /// \brief This method calculates the average and fluctuation of the baseline in the
 /// specified range and writes the values to fBaseLine and fBaseLineSigma respectively.
@@ -792,12 +836,16 @@ void TRestRawSignal::CalculateBaseLineMedian(Int_t startBin, Int_t endBin) {
 ///
 /// \param option By setting this option to "ROBUST", the average is calculated as median,
 /// and the fluctuation as interquartile range (IQR), which are less affected by outliers (e.g. a signal
-/// pulse).
+/// pulse). By setting it to "OUTLIERS" the median and sigma will only take into account the
+/// 25-75% values of the interval.
 ///
 void TRestRawSignal::CalculateBaseLine(Int_t startBin, Int_t endBin, const std::string& option) {
     if (ToUpper(option) == "ROBUST") {
         CalculateBaseLineMedian(startBin, endBin);
         CalculateBaseLineSigmaIQR(startBin, endBin);
+    } else if (ToUpper(option) == "OUTLIERS") {
+        CalculateBaseLineMedianExcludeOutliers(startBin, endBin);
+        CalculateBaseLineSigmaExcludeOutliers(startBin, endBin);
     } else {
         CalculateBaseLineMean(startBin, endBin);
         CalculateBaseLineSigmaSD(startBin, endBin);
@@ -842,6 +890,54 @@ void TRestRawSignal::CalculateBaseLineSigmaIQR(Int_t startBin, Int_t endBin) {
     }
 }
 
+///////////////////////////////////////////////
+/// \brief This method is called by CalculateBaseLine with the "OUTLIERS"-option to
+/// determine the value of the baseline
+/// fluctuation as the standard deviation in the baseline range provided excluding outliers.
+/// Since outliers are strongly suppressed there is no need to use the IQR.
+///
+void TRestRawSignal::CalculateBaseLineSigmaExcludeOutliers(Int_t startBin, Int_t endBin) {
+    if (endBin - startBin <= 0) {
+        fBaseLineSigma = 0.;
+        return;
+    } else if (endBin > static_cast<int>(fSignalData.size())) {
+        cout << "TRestRawSignal::CalculateBaseLineSigma. Error! Range exceeds the rawdata depth!!" << endl;
+        endBin = fSignalData.size();
+    } else {
+        // Extract the data within the interval
+        std::vector<Short_t> data(fSignalData.begin() + startBin, fSignalData.begin() + endBin);
+        std::sort(data.begin(), data.end());
+
+        // Calculate Q1 and Q3 for IQR
+        size_t dataSize = data.size();
+        Short_t Q1 = data[dataSize / 4];
+        Short_t Q3 = data[(3 * dataSize) / 4];
+        Double_t lowerBound = Q1;
+        Double_t upperBound = Q3;
+
+        // Filter out the outliers
+        std::vector<Short_t> filteredData;
+        for (const auto& value : data) {
+            if (value >= lowerBound && value <= upperBound) {
+                filteredData.emplace_back(value);
+            }
+        }
+
+        // Calculate standard deviation of filtered data
+        if (filteredData.empty()) {
+            fBaseLineSigma = 0.;  // If all values are outliers, set sigma to zero
+        } else {
+            double mean =
+                std::accumulate(filteredData.begin(), filteredData.end(), 0.0) / filteredData.size();
+            double variance = 0.0;
+            for (const auto& value : filteredData) {
+                variance += std::pow(value - mean, 2);
+            }
+            fBaseLineSigma = std::sqrt(variance / filteredData.size());  // Standard deviation
+        }
+    }
+}
+
 ///////////////////////////////////////////////
 /// \brief This method adds an offset to the signal data
 ///