Merge pull request #19 from rest-for-physics/lobis-detector-to-raw-calibration

lobis · web-flow · commit 319cc970317b · 2022-09-05T11:18:11.000+02:00
Update `TRestDetectorSignalToRawSignalProcess` to support calibration points and shaping
diff --git a/inc/TRestDetectorSignalToRawSignalProcess.h b/inc/TRestDetectorSignalToRawSignalProcess.h
@@ -24,10 +24,9 @@
 #define RestCore_TRestDetectorSignalToRawSignalProcess
 
 #include <TRestDetectorSignalEvent.h>
+#include <TRestEventProcess.h>
 #include <TRestRawSignalEvent.h>
 
-#include "TRestEventProcess.h"
-
 //! A process to convert a TRestDetectorSignalEvent into a TRestRawSignalEvent
 class TRestDetectorSignalToRawSignalProcess : public TRestEventProcess {
    private:
@@ -52,14 +51,34 @@ class TRestDetectorSignalToRawSignalProcess : public TRestEventProcess {
     TString fTriggerMode = "firstDeposit";
 
     /// The number of time bins the time start is delayed in the resulting output signal.
-    Int_t fTriggerDelay = 100;  // ns
+    Int_t fTriggerDelay = 100;
+
+    /// The starting time for the "fixed" trigger mode (can be offset by the trigger delay)
+    Int_t fTriggerFixedStartTime = 0;
 
-    /// A factor the data values will be multiplied by at the output signal.
-    Double_t fGain = 100.0;
+    /// fCalibrationGain and fCalibrationOffset define the linear calibration.
+    /// output = input * fCalibrationGain + calibrationOffset
+    Double_t fCalibrationGain = 100.0;
+    Double_t fCalibrationOffset = 0.0;  // adc units
 
     /// This parameter is used by integralWindow trigger mode to define the acquisition window.
     Double_t fIntegralThreshold = 1229.0;
 
+    /// two distinct energy values used for calibration
+    TVector2 fCalibrationEnergy = TVector2(0.0, 0.0);
+    /// position in the range corresponding to the energy in 'fCalibrationEnergy'. Values between 0 and 1
+    TVector2 fCalibrationRange = TVector2(0.0, 0.0);
+    /// Usage: fCalibrationEnergy = (0, 100 MeV) and fCalibrationRange = (0.1, 0.9)
+    /// will perform a linear calibration with 0 equal to 0.1 of the range (0.1 * (max - min) + min) and 100
+    /// MeV equal to 0.9 of the range. The range is the one corresponding to a Short_t for rawsignal.
+
+    /// If defined ( > 0 ) we will compute the sin shaping of the signal, this is done in this process to
+    /// avoid artifacts in the signal (e.g. signals not getting cut when they should)
+    Double_t fShapingTime = 0.0;  // us
+
+    Double_t fTimeStart;  //!
+    Double_t fTimeEnd;    //!
+
    public:
     inline Double_t GetSampling() const { return fSampling; }
     inline void SetSampling(Double_t sampling) { fSampling = sampling; }
@@ -73,15 +92,36 @@ class TRestDetectorSignalToRawSignalProcess : public TRestEventProcess {
     inline Int_t GetTriggerDelay() const { return fTriggerDelay; }
     inline void SetTriggerDelay(Int_t triggerDelay) { fTriggerDelay = triggerDelay; }
 
-    inline Double_t GetGain() const { return fGain; }
-    inline void SetGain(Double_t gain) { fGain = gain; }
+    inline Double_t GetGain() const { return fCalibrationGain; }
+    inline void SetGain(Double_t gain) { fCalibrationGain = gain; }
+
+    inline Double_t GetCalibrationOffset() const { return fCalibrationOffset; }
+    inline void SetCalibrationOffset(Double_t offset) { fCalibrationOffset = offset; }
 
     inline Double_t GetIntegralThreshold() const { return fIntegralThreshold; }
     inline void SetIntegralThreshold(Double_t integralThreshold) { fIntegralThreshold = integralThreshold; }
 
+    inline Double_t GetShapingTime() const { return fShapingTime; }
+    inline void SetShapingTime(Double_t shapingTime) { fShapingTime = shapingTime; }
+
+    inline bool IsShapingEnabled() const { return fShapingTime > 0; }
+
+    inline bool IsLinearCalibration() const {
+        // Will return true if two points have been given for calibration
+        return (fCalibrationEnergy.Mod() != 0 && fCalibrationRange.Mod() != 0);
+    }
+
+    inline TVector2 GetCalibrationEnergy() const { return fCalibrationEnergy; }
+    inline void SetCalibrationEnergy(TVector2 calibrationEnergy) { fCalibrationEnergy = calibrationEnergy; }
+
+    inline TVector2 GetCalibrationRange() const { return fCalibrationRange; }
+    inline void SetCalibrationRange(TVector2 calibrationRange) { fCalibrationRange = calibrationRange; }
+
     any GetInputEvent() const override { return fInputSignalEvent; }
     any GetOutputEvent() const override { return fOutputRawSignalEvent; }
 
+    void InitProcess() override;
+
     TRestEvent* ProcessEvent(TRestEvent* inputEvent) override;
 
     void LoadConfig(const std::string& configFilename, const std::string& name = "");
@@ -94,7 +134,19 @@ class TRestDetectorSignalToRawSignalProcess : public TRestEventProcess {
         RESTMetadata << "Points per channel : " << fNPoints << RESTendl;
         RESTMetadata << "Trigger mode : " << fTriggerMode << RESTendl;
         RESTMetadata << "Trigger delay : " << fTriggerDelay << " time units" << RESTendl;
-        RESTMetadata << "ADC gain : " << fGain << RESTendl;
+
+        if (IsLinearCalibration()) {
+            RESTMetadata << "Calibration energy : (" << fCalibrationEnergy.X() << ", "
+                         << fCalibrationEnergy.Y() << ") keV" << RESTendl;
+            RESTMetadata << "Calibration range : (" << fCalibrationRange.X() << ", " << fCalibrationRange.Y()
+                         << ")" << RESTendl;
+        }
+        RESTMetadata << "ADC Gain : " << fCalibrationGain << RESTendl;
+        RESTMetadata << "ADC Offset : " << fCalibrationOffset << RESTendl;
+
+        if (IsShapingEnabled()) {
+            RESTMetadata << "Shaping time : " << fShapingTime << " us" << RESTendl;
+        }
 
         EndPrintProcess();
     }
@@ -112,6 +164,6 @@ class TRestDetectorSignalToRawSignalProcess : public TRestEventProcess {
     // Destructor
     ~TRestDetectorSignalToRawSignalProcess();
 
-    ClassDefOverride(TRestDetectorSignalToRawSignalProcess, 2);
+    ClassDefOverride(TRestDetectorSignalToRawSignalProcess, 3);
 };
 #endif
diff --git a/src/TRestDetectorSignalToRawSignalProcess.cxx b/src/TRestDetectorSignalToRawSignalProcess.cxx
@@ -82,11 +82,11 @@
 /// \htmlonly <style>div.image img[src="trigger.png"]{width:500px;}</style> \endhtmlonly
 ///
 /// The following figure shows the integralThreshold trigger definition for a NLDBD
-/// event. tStart and tEnd define the acquision window, centered on the time
-/// when the signal integral is above the threshold defined. tStart has been
+/// event. fTimeStart and fTimeEnd define the acquisition window, centered on the time
+/// when the signal integral is above the threshold defined. fTimeStart has been
 /// shifted by a triggerDelay = 60 samples * 200ns
 ///
-/// ![An ilustration of the trigger definition](trigger.png)
+/// ![An illustration of the trigger definition](trigger.png)
 ///
 /// * **triggerDelay**: The time start obtained by the trigger mode
 /// definition can be shifted using this parameter. The shift is
@@ -113,7 +113,9 @@
 ///
 /// <hr>
 ///
+
 #include "TRestDetectorSignalToRawSignalProcess.h"
+
 using namespace std;
 
 ClassImp(TRestDetectorSignalToRawSignalProcess);
@@ -181,115 +183,130 @@ void TRestDetectorSignalToRawSignalProcess::Initialize() {
 TRestEvent* TRestDetectorSignalToRawSignalProcess::ProcessEvent(TRestEvent* inputEvent) {
     fInputSignalEvent = (TRestDetectorSignalEvent*)inputEvent;
 
-    if (fInputSignalEvent->GetNumberOfSignals() <= 0) return nullptr;
+    if (fInputSignalEvent->GetNumberOfSignals() <= 0) {
+        return nullptr;
+    }
 
-    if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug) fOutputRawSignalEvent->PrintEvent();
+    if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug) {
+        fOutputRawSignalEvent->PrintEvent();
+    }
 
     fOutputRawSignalEvent->SetID(fInputSignalEvent->GetID());
     fOutputRawSignalEvent->SetSubID(fInputSignalEvent->GetSubID());
     fOutputRawSignalEvent->SetTimeStamp(fInputSignalEvent->GetTimeStamp());
     fOutputRawSignalEvent->SetSubEventTag(fInputSignalEvent->GetSubEventTag());
 
-    // The time event window is defined (tStart, tEnd)
-    Double_t tStart = std::numeric_limits<double>::quiet_NaN();
-    Double_t tEnd = 10000;
     if (fTriggerMode == "firstDeposit") {
-        tStart = fInputSignalEvent->GetMinTime() - fTriggerDelay * fSampling;
-        tEnd = fInputSignalEvent->GetMinTime() + (fNPoints - fTriggerDelay) * fSampling;
+        fTimeStart = fInputSignalEvent->GetMinTime() - fTriggerDelay * fSampling;
+        fTimeEnd = fInputSignalEvent->GetMinTime() + (fNPoints - fTriggerDelay) * fSampling;
     } else if (fTriggerMode == "integralThreshold") {
-        Double_t maxT = fInputSignalEvent->GetMaxTime();
-        Double_t minT = fInputSignalEvent->GetMinTime();
-
-        for (Double_t t = minT - fNPoints * fSampling; t <= maxT + fNPoints * fSampling; t = t + 0.5) {
-            Double_t en = fInputSignalEvent->GetIntegralWithTime(t, t + (fSampling * fNPoints) / 2.);
-
-            if (en > fIntegralThreshold) {
-                tStart = t - fTriggerDelay * fSampling;
-                tEnd = t + (fNPoints - fTriggerDelay) * fSampling;
+        bool thresholdReached = false;
+        for (Double_t t = fInputSignalEvent->GetMinTime() - fNPoints * fSampling;
+             t <= fInputSignalEvent->GetMaxTime() + fNPoints * fSampling; t = t + 0.5) {
+            Double_t energy = fInputSignalEvent->GetIntegralWithTime(t, t + (fSampling * fNPoints) / 2.);
+
+            if (energy > fIntegralThreshold) {
+                fTimeStart = t - fTriggerDelay * fSampling;
+                fTimeEnd = t + (fNPoints - fTriggerDelay) * fSampling;
+                thresholdReached = true;
             }
         }
-    } else {
-        if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Warning) {
-            cout << "REST WARNING. TRestDetectorSignalToRawSignalProcess. fTriggerMode not "
-                    "recognized!"
-                 << endl;
-            cout << "Setting fTriggerMode to firstDeposit" << endl;
-        }
-
-        fTriggerMode = "firstDeposit";
-
-        tStart = fInputSignalEvent->GetMinTime() - fTriggerDelay * fSampling;
-        tEnd = fInputSignalEvent->GetMinTime() + (fNPoints - fTriggerDelay) * fSampling;
-    }
-
-    if (std::isnan(tStart)) {
-        if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Warning) {
-            cout << endl;
-            cout << "REST WARNING. TRestDetectorSignalToRawSignalProcess. tStart was not "
-                    "defined."
-                 << endl;
-            cout << "Setting tStart = 0. tEnd = fNPoints * fSampling" << endl;
-            cout << endl;
+        if (!thresholdReached) {
+            RESTWarning << "Integral threshold for trigger not reached" << RESTendl;
+            fTimeStart = 0;
+            fTimeEnd = fNPoints * fSampling;
         }
-
-        tEnd = fNPoints * fSampling;
     }
 
-    if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug) {
-        cout << "tStart : " << tStart << " us " << endl;
-        cout << "tEnd : " << tEnd << " us " << endl;
-    }
+    RESTDebug << "fTimeStart : " << fTimeStart << " us " << RESTendl;
+    RESTDebug << "fTimeEnd : " << fTimeEnd << " us " << RESTendl;
 
     for (int n = 0; n < fInputSignalEvent->GetNumberOfSignals(); n++) {
-        vector<Double_t> sData(fNPoints);
-        for (int i = 0; i < fNPoints; i++) sData[i] = 0;
-
-        TRestDetectorSignal* sgnl = fInputSignalEvent->GetSignal(n);
-        Int_t signalID = sgnl->GetSignalID();
+        vector<Double_t> data(fNPoints, fCalibrationOffset);
+        TRestDetectorSignal* signal = fInputSignalEvent->GetSignal(n);
+        Int_t signalID = signal->GetSignalID();
 
-        for (int m = 0; m < sgnl->GetNumberOfPoints(); m++) {
-            Double_t t = sgnl->GetTime(m);
-            Double_t d = sgnl->GetData(m);
+        for (int m = 0; m < signal->GetNumberOfPoints(); m++) {
+            Double_t t = signal->GetTime(m);
+            Double_t d = signal->GetData(m);
 
-            if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug && n < 3 && m < 5)
+            if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug && n < 3 && m < 5) {
                 cout << "Signal : " << n << " Sample : " << m << " T : " << t << " Data : " << d << endl;
+            }
 
-            if (t > tStart && t < tEnd) {
+            if (t > fTimeStart && t < fTimeEnd) {
                 // convert physical time (in us) to timeBin
-                Int_t timeBin = (Int_t)round((t - tStart) / fSampling);
+                Int_t timeBin = (Int_t)round((t - fTimeStart) / fSampling);
 
-                if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Warning)
+                if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Warning) {
                     if (timeBin < 0 || timeBin > fNPoints) {
                         cout << "Time bin out of range!!! bin value : " << timeBin << endl;
                         timeBin = 0;
                     }
+                }
 
-                RESTDebug << "Adding data : " << sgnl->GetData(m) << " to Time Bin : " << timeBin << RESTendl;
-                sData[timeBin] += fGain * sgnl->GetData(m);
+                RESTDebug << "Adding data : " << signal->GetData(m) << " to Time Bin : " << timeBin
+                          << RESTendl;
+                data[timeBin] += fCalibrationGain * signal->GetData(m);
             }
         }
 
-        if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Warning)
-            for (int x = 0; x < fNPoints; x++)
-                if (sData[x] < -32768. || sData[x] > 32768.)
-                    cout << "REST Warning : data is outside short range : " << sData[x] << endl;
+        if (IsShapingEnabled()) {
+            const auto shapingFunction = [](Double_t t) -> Double_t {
+                if (t <= 0) {
+                    return 0;
+                }
+                // function is normalized such that its absolute maximum is 1.0
+                // max is at x = 1.1664004483744728
+                return TMath::Exp(-3.0 * t) * TMath::Power(t, 3.0) * TMath::Sin(t) * 22.68112123672292;
+            };
+
+            vector<Double_t> dataAfterShaping(fNPoints, fCalibrationOffset);
+            for (int i = 0; i < fNPoints; i++) {
+                const Double_t value = data[i] - fCalibrationOffset;
+                if (value <= 0) {
+                    // Only positive values are possible, 0 means no signal in this bin
+                    continue;
+                }
+                for (int j = 0; j < fNPoints; j++) {
+                    dataAfterShaping[j] += value * shapingFunction(((j - i) * fSampling) / fShapingTime);
+                }
+            }
+            data = dataAfterShaping;
+        }
 
-        TRestRawSignal rSgnl;
-        rSgnl.SetSignalID(signalID);
+        TRestRawSignal rawSignal;
+        rawSignal.SetSignalID(signalID);
         for (int x = 0; x < fNPoints; x++) {
-            if (sData[x] < -32768. || sData[x] > 32768.) fOutputRawSignalEvent->SetOK(false);
+            double value = round(data[x]);
+            if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Warning) {
+                if (value < numeric_limits<Short_t>::min() || value > numeric_limits<Short_t>::max()) {
+                    RESTWarning << "value (" << value << ") is outside short range ("
+                                << numeric_limits<Short_t>::min() << ", " << numeric_limits<Short_t>::max()
+                                << ")" << RESTendl;
+                }
+            }
 
-            Short_t value = (Short_t)round(sData[x]);
-            rSgnl.AddPoint(value);
+            if (value < numeric_limits<Short_t>::min()) {
+                value = numeric_limits<Short_t>::min();
+                fOutputRawSignalEvent->SetOK(false);
+            } else if (value > numeric_limits<Short_t>::max()) {
+                value = numeric_limits<Short_t>::max();
+                fOutputRawSignalEvent->SetOK(false);
+            }
+            rawSignal.AddPoint((Short_t)value);
         }
 
-        if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug) rSgnl.Print();
+        if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug) {
+            rawSignal.Print();
+        }
         RESTDebug << "Adding signal to raw signal event" << RESTendl;
-        fOutputRawSignalEvent->AddSignal(rSgnl);
+
+        fOutputRawSignalEvent->AddSignal(rawSignal);
     }
 
     RESTDebug << "TRestDetectorSignalToRawSignalProcess. Returning event with N signals "
-          << fOutputRawSignalEvent->GetNumberOfSignals() << RESTendl;
+              << fOutputRawSignalEvent->GetNumberOfSignals() << RESTendl;
 
     return fOutputRawSignalEvent;
 }
@@ -299,10 +316,48 @@ TRestEvent* TRestDetectorSignalToRawSignalProcess::ProcessEvent(TRestEvent* inpu
 /// TRestDetectorSignalToRawSignalProcess metadata section
 ///
 void TRestDetectorSignalToRawSignalProcess::InitFromConfigFile() {
-    fSampling = GetDblParameterWithUnits("sampling");
-    fNPoints = StringToInteger(GetParameter("Npoints", "512"));
-    fTriggerMode = GetParameter("triggerMode", "firstDeposit");
-    fTriggerDelay = StringToInteger(GetParameter("triggerDelay", 100));
-    fGain = StringToDouble(GetParameter("gain", "100"));
-    fIntegralThreshold = StringToDouble(GetParameter("integralThreshold", "1229"));
+    auto nPoints = GetParameter("nPoints");
+    if (nPoints == PARAMETER_NOT_FOUND_STR) {
+        nPoints = GetParameter("Npoints", fNPoints);
+    }
+    fNPoints = StringToInteger(nPoints);
+
+    fTriggerMode = GetParameter("triggerMode", fTriggerMode);
+    const set<string> validTriggerModes = {"firstDeposit", "integralThreshold", "fixed"};
+    if (validTriggerModes.count(fTriggerMode.Data()) == 0) {
+        RESTError << "Trigger mode set to: '" << fTriggerMode
+                  << "' which is not a valid trigger mode. Please use one of the following trigger modes: ";
+        for (const auto& triggerMode : validTriggerModes) {
+            RESTError << triggerMode << " ";
+        }
+        RESTError << RESTendl;
+        exit(1);
+    }
+
+    fSampling = GetDblParameterWithUnits("sampling", fSampling);
+    fTriggerDelay = StringToInteger(GetParameter("triggerDelay", fTriggerDelay));
+    fIntegralThreshold = StringToDouble(GetParameter("integralThreshold", fIntegralThreshold));
+    fTriggerFixedStartTime = GetDblParameterWithUnits("triggerFixedStartTime", fTriggerFixedStartTime);
+
+    fCalibrationGain = StringToDouble(GetParameter("gain", fCalibrationGain));
+    fCalibrationOffset = StringToDouble(GetParameter("offset", fCalibrationOffset));
+    fCalibrationEnergy = Get2DVectorParameterWithUnits("calibrationEnergy", fCalibrationEnergy);
+    fCalibrationRange = Get2DVectorParameterWithUnits("calibrationRange", fCalibrationRange);
+
+    if (IsLinearCalibration()) {
+        const auto range = numeric_limits<Short_t>::max() - numeric_limits<Short_t>::min();
+        fCalibrationGain = range * (fCalibrationRange.Y() - fCalibrationRange.X()) /
+                           (fCalibrationEnergy.Y() - fCalibrationEnergy.X());
+        fCalibrationOffset = range * (fCalibrationRange.X() - fCalibrationGain * fCalibrationEnergy.X()) +
+                             numeric_limits<Short_t>::min();
+    }
+
+    fShapingTime = GetDblParameterWithUnits("shapingTime", fShapingTime);
+}
+
+void TRestDetectorSignalToRawSignalProcess::InitProcess() {
+    if (fTriggerMode == "fixed") {
+        fTimeStart = fTriggerFixedStartTime - fTriggerDelay * fSampling;
+        fTimeEnd = fTimeStart + fNPoints * fSampling;
+    }
 }
