Adding the possibility to perform smoothing on TRestRawToDetectorSignalProcess

inc/TRestRawToDetectorSignalProcess.h

@@ -73,13 +73,20 @@ class TRestRawToDetectorSignalProcess : public TRestEventProcess {
     /// A parameter to determine if baseline correction has been applied by a previous process
     Bool_t fBaseLineCorrection = false;
 
+    // Perform signal smoothing to the data
+    Bool_t fSignalSmoothing = false;
+
+    // In case signal smoothing is active provide averaging points
+    Int_t fAveragingPoints = 5;
+
    public:
     RESTValue GetInputEvent() const override { return fInputSignalEvent; }
     RESTValue GetOutputEvent() const override { return fOutputSignalEvent; }
 
     TRestEvent* ProcessEvent(TRestEvent* inputEvent) override;
 
-    void ZeroSuppresion(TRestRawSignal* rawSignal, TRestDetectorSignal& sgnl);
+    void ProcessSignal(TRestRawSignal* rawSignal, TRestDetectorSignal& sgnl);
+    void ProcessSignalSmoothed(TRestRawSignal* rawSignal, TRestDetectorSignal& sgnl);
 
     /// It prints out the process parameters stored in the metadata structure
     void PrintMetadata() override {
@@ -90,6 +97,7 @@ class TRestRawToDetectorSignalProcess : public TRestEventProcess {
         RESTMetadata << "Gain : " << fGain << RESTendl;
 
         if (fZeroSuppression) {
+            RESTMetadata << "ZeroSuppression is enabled " << RESTendl;
             RESTMetadata << "Base line range definition : ( " << fBaseLineRange.X() << " , "
                          << fBaseLineRange.Y() << " ) " << RESTendl;
             RESTMetadata << "Integral range : ( " << fIntegralRange.X() << " , " << fIntegralRange.Y()
@@ -99,8 +107,11 @@ class TRestRawToDetectorSignalProcess : public TRestEventProcess {
             RESTMetadata << "Number of points over threshold : " << fNPointsOverThreshold << RESTendl;
         }
 
-        if (fBaseLineCorrection)
-            RESTMetadata << "BaseLine correction is enabled for TRestRawSignalAnalysisProcess" << RESTendl;
+        if (fBaseLineCorrection) RESTMetadata << "BaseLine correction is enabled" << RESTendl;
+        if (fSignalSmoothing) {
+            RESTMetadata << "Signal smoothing is enabled" << RESTendl;
+            RESTMetadata << "Averaging points " << fAveragingPoints << RESTendl;
+        }
 
         EndPrintProcess();
     }
@@ -117,6 +128,6 @@ class TRestRawToDetectorSignalProcess : public TRestEventProcess {
     // Destructor
     ~TRestRawToDetectorSignalProcess();
 
-    ClassDefOverride(TRestRawToDetectorSignalProcess, 2);
+    ClassDefOverride(TRestRawToDetectorSignalProcess, 3);
 };
 #endif

src/TRestRawToDetectorSignalProcess.cxx

@@ -55,6 +55,10 @@
 /// * **signalThreshold**: The number of sigmas a set of consecutive points
 /// identified over threshold must be over the baseline fluctuations to be
 /// finally considered a physical signal.
+/// * **signalSmoothing**: Perform smoothing of the signal before saving as
+/// a TRestDetectorSignalEvent
+/// * **averagingPoints**: In case smoothing is performed set the number of
+/// points to be averaged to perform the smoothing
 ///
 /// List of observables:
 ///
@@ -69,7 +73,8 @@
 /// // A raw signal with 200ns binning will be translated to a
 /// // TRestDetectorSignalEvent. The new signal will start at time=20us, and its
 /// // amplitude will be reduced a factor 50. If zeroSuppression is true it will
-/// // perform
+/// // perform the substraction of the points out of the range of interest. A
+/// // smoothing of the signal will be performed using 3 averaging points.
 ///
 /// <TRestRawToDetectorSignalProcess name="rsTos" title"Raw signal to signal">
 ///     <parameter name="sampling" value="0.2" units="us" />
@@ -81,6 +86,8 @@
 ///     <parameter name="pointThreshold" value="3"/>
 ///     <parameter name="signalThreshold" value="7"/>
 ///     <parameter name="nPointsOverThreshold" value="7"/>
+///     <parameter name="signalSmoothing" value="false"/>
+///     <parameter name="averagingPoints" value="3"/>
 ///     <observable name="NSignalsRejected" value="ON"/>
 ///
 /// </TRestRawToDetectorSignalProcess>
@@ -110,6 +117,8 @@
 ///             Javier Galan
 /// 2022-January: Adding ZeroSuppression method
 ///             JuanAn Garcia
+/// 2023-April: Added new method to process smoothed data
+///             JuanAn Garcia
 ///
 /// \class      TRestRawToDetectorSignalProcess
 /// \author     Javier Gracia
@@ -119,6 +128,8 @@
 /// <hr>
 ///
 #include "TRestRawToDetectorSignalProcess.h"
+
+#include "TRestPulseShapeAnalysis.h"
 using namespace std;
 
 ClassImp(TRestRawToDetectorSignalProcess);
@@ -163,18 +174,17 @@ TRestEvent* TRestRawToDetectorSignalProcess::ProcessEvent(TRestEvent* inputEvent
         TRestRawSignal* rawSignal = fInputSignalEvent->GetSignal(n);
         signal.SetID(rawSignal->GetID());
 
-        if (fZeroSuppression) {
-            ZeroSuppresion(rawSignal, signal);
+        if (fSignalSmoothing) {
+            ProcessSignalSmoothed(rawSignal, signal);
         } else {
-            for (int p = 0; p < rawSignal->GetNumberOfPoints(); p++)
-                if (rawSignal->GetData(p) > fThreshold)
-                    signal.NewPoint(fTriggerStarts + fSampling * p, fGain * rawSignal->GetData(p));
+            ProcessSignal(rawSignal, signal);
         }
 
-        if (signal.GetNumberOfPoints() > 0)
+        if (signal.GetNumberOfPoints() > 0) {
             fOutputSignalEvent->AddSignal(signal);
-        else
+        } else {
             rejectedSignal++;
+        }
     }
 
     SetObservableValue("NSignalsRejected", rejectedSignal);
@@ -184,13 +194,47 @@ TRestEvent* TRestRawToDetectorSignalProcess::ProcessEvent(TRestEvent* inputEvent
     return fOutputSignalEvent;
 }
 
-void TRestRawToDetectorSignalProcess::ZeroSuppresion(TRestRawSignal* rawSignal, TRestDetectorSignal& sgnl) {
-    rawSignal->InitializePointsOverThreshold(TVector2(fPointThreshold, fSignalThreshold),
-                                             fNPointsOverThreshold, 512);
+///////////////////////////////////////////////
+/// \brief This function performs the smoothing
+/// of the TRestRawSignalEvent before transforming
+/// on a TRestDetectorSignalEvent
+///
+void TRestRawToDetectorSignalProcess::ProcessSignalSmoothed(TRestRawSignal* rawSignal,
+                                                            TRestDetectorSignal& sgnl) {
+    auto smoothed = TRestPulseShapeAnalysis::GetSignalSmoothed(rawSignal->GetData(), fAveragingPoints);
+    const double baselineSigma = rawSignal->GetBaseLineSigma();
+    if (fZeroSuppression) {
+        auto pOver = TRestPulseShapeAnalysis::GetPointsOverThreshold(
+            smoothed, fIntegralRange, TVector2(fPointThreshold, fSignalThreshold), fNPointsOverThreshold, 512,
+            baselineSigma);
+        for (const auto& [j, data] : pOver) {
+            sgnl.NewPoint(fTriggerStarts + fSampling * j, fGain * data);
+        }
+    } else {
+        for (size_t p = 0; p < smoothed.size(); p++) {
+            const double data = smoothed[p];
+            if (data > fThreshold) sgnl.NewPoint(fTriggerStarts + fSampling * p, fGain * data);
+        }
+    }
+}
 
-    std::vector<Int_t> pOver = rawSignal->GetPointsOverThreshold();
-    for (unsigned int n = 0; n < pOver.size(); n++) {
-        int j = pOver[n];
-        sgnl.NewPoint(fTriggerStarts + fSampling * j, fGain * rawSignal->GetData(j));
+///////////////////////////////////////////////
+/// \brief This function transform a TRestRawSignalEvent
+/// into a TRestDetectorSignalEvent without
+/// any processing, despite Zero suppression
+///
+void TRestRawToDetectorSignalProcess::ProcessSignal(TRestRawSignal* rawSignal, TRestDetectorSignal& sgnl) {
+    if (fZeroSuppression) {
+        rawSignal->InitializePointsOverThreshold(TVector2(fPointThreshold, fSignalThreshold),
+                                                 fNPointsOverThreshold, 512);
+
+        auto pOver = rawSignal->GetPointsOverThreshold();
+        for (const auto& [j, data] : pOver) {
+            sgnl.NewPoint(fTriggerStarts + fSampling * j, fGain * data);
+        }
+    } else {
+        for (int p = 0; p < rawSignal->GetNumberOfPoints(); p++)
+            if (rawSignal->GetData(p) > fThreshold)
+                sgnl.NewPoint(fTriggerStarts + fSampling * p, fGain * rawSignal->GetData(p));
     }
 }