Merge pull request #114 from rest-for-physics/lobis-fix-noise

Prevent noise from going out of bounds

inc/TRestRawSignal.h

@@ -24,15 +24,16 @@
 #define RestCore_TRestRawSignal
 
 #include <TGraph.h>
-#include <TObject.h>
+#include <TRandom.h>
 #include <TString.h>
 #include <TVector2.h>
 
 #include <iostream>
 #include <string>
+#include <vector>
 
 //! It defines a Short_t array with a physical parameter that evolves in time using a fixed time bin.
-class TRestRawSignal : public TObject {
+class TRestRawSignal {
    private:
     void CalculateThresholdIntegral();
 
@@ -51,6 +52,9 @@ class TRestRawSignal : public TObject {
 
     Bool_t fShowWarnings = true;
 
+    /// Seed used for random number generation
+    UInt_t fSeed = gRandom->GetSeed();
+
    public:
     /// A TGraph pointer used to store the TRestRawSignal drawing
     TGraph* fGraph;  //!
@@ -112,10 +116,7 @@ class TRestRawSignal : public TObject {
     inline TVector2 GetRange() const { return fRange; }
 
     /// Returns false if the baseline and its baseline fluctuation was not initialized.
-    inline Bool_t isBaseLineInitialized() {
-        if (fBaseLineSigma == 0 && fBaseLine == 0) return false;
-        return true;
-    }
+    inline Bool_t isBaseLineInitialized() const { return !(fBaseLineSigma == 0 && fBaseLine == 0); }
 
     Double_t GetData(Int_t n) const;
 
@@ -147,16 +148,16 @@ class TRestRawSignal : public TObject {
 
     void Initialize();
 
-    void AddPoint(Short_t d);
-
-    void AddCharge(Short_t d);
+    void AddPoint(Short_t);
 
-    void AddDeposit(Short_t d);
+    void AddPoint(Double_t);
 
     void IncreaseBinBy(Int_t bin, Double_t data);
 
     void InitializePointsOverThreshold(const TVector2& thrPar, Int_t nPointsOver, Int_t nPointsFlat = 512);
 
+    UInt_t GetSeed() const { return fSeed; }
+
     Double_t GetIntegral();
 
     Double_t GetIntegralInRange(Int_t startBin, Int_t endBin);
@@ -211,12 +212,14 @@ class TRestRawSignal : public TObject {
 
     void Print() const;
 
+    void SetSeed(UInt_t seed) { fSeed = seed; }
+
     TGraph* GetGraph(Int_t color = 1);
 
     TRestRawSignal();
     TRestRawSignal(Int_t nBins);
     ~TRestRawSignal();
 
-    ClassDef(TRestRawSignal, 1);
+    ClassDef(TRestRawSignal, 2);
 };
 #endif

inc/TRestRawSignalEvent.h

@@ -59,6 +59,8 @@ class TRestRawSignalEvent : public TRestEvent {
         return true;
     }
 
+    std::vector<TRestRawSignal> GetSignals() const { return fSignal; }
+
     // Setters
     void AddSignal(TRestRawSignal& s);
 

pipeline/processes/fit/GeneralFit.C

@@ -1,8 +1,19 @@
+
+#include <TRestRawSignal.h>
+#include <TRestRawSignalAddNoiseProcess.h>
+#include <TRestRawSignalEvent.h>
+#include <TRestRawSignalFittingProcess.h>
+#include <TRestRawSignalShapingProcess.h>
+
+using namespace std;
+
 Int_t GeneralFit(Bool_t draw = false) {
     TRestRawSignalEvent* ev = new TRestRawSignalEvent();
 
     TRestRawSignal* sgnl = new TRestRawSignal();
-    for (int n = 0; n < 512; n++) sgnl->AddPoint(0);
+    for (int n = 0; n < 512; n++) {
+        sgnl->AddPoint((Double_t)0);
+    }
     sgnl->IncreaseBinBy(70, 100);
     ev->AddSignal(*sgnl);
 

pipeline/processes/fit/fit.C

@@ -1,3 +1,12 @@
+
+#include <TRestRawSignal.h>
+#include <TRestRawSignalAddNoiseProcess.h>
+#include <TRestRawSignalEvent.h>
+#include <TRestRawSignalFittingProcess.h>
+#include <TRestRawSignalShapingProcess.h>
+
+using namespace std;
+
 Int_t fit(Bool_t draw = false) {
     TRestRawSignalEvent* ev = new TRestRawSignalEvent();
 

src/TRestRawMemoryBufferToSignalProcess.cxx

@@ -304,19 +304,19 @@ TRestEvent* TRestRawMemoryBufferToSignalProcess::ProcessEvent(TRestEvent* inputE
             SemaphoreRed(fSemaphoreId);
 
             for (unsigned int s = 0; s < fShMem_daqInfo->nSignals; s++) {
-                TRestRawSignal sgnl;
-                sgnl.SetSignalID(fShMem_Buffer[s * (maxSamples + 1)]);
+                TRestRawSignal signal;
+                signal.SetSignalID(fShMem_Buffer[s * (maxSamples + 1)]);
 
                 if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug)
-                    cout << "s : " << s << " id : " << sgnl.GetSignalID() << endl;
+                    cout << "s : " << s << " id : " << signal.GetSignalID() << endl;
 
                 for (int n = 0; n < maxSamples; n++) {
-                    sgnl.AddPoint(fShMem_Buffer[s * (maxSamples + 1) + 1 + n]);
+                    signal.AddPoint((Short_t)fShMem_Buffer[s * (maxSamples + 1) + 1 + n]);
                 }
-                fOutputRawSignalEvent->AddSignal(sgnl);
+                fOutputRawSignalEvent->AddSignal(signal);
 
                 if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Extreme) {
-                    sgnl.Print();
+                    signal.Print();
                     GetChar();
                 }
             }

src/TRestRawSignal.cxx

@@ -126,19 +126,19 @@ void TRestRawSignal::Reset() {
 ///////////////////////////////////////////////
 /// \brief Adds a new point to the end of the signal data array
 ///
-void TRestRawSignal::AddPoint(Short_t d) { fSignalData.push_back(d); }
+void TRestRawSignal::AddPoint(Short_t value) { fSignalData.push_back(value); }
 
 ///////////////////////////////////////////////
-/// \brief Adds a new point to the end of the signal data array. Same as
-/// AddPoint.
-///
-void TRestRawSignal::AddCharge(Short_t d) { AddPoint(d); }
-
-///////////////////////////////////////////////
-/// \brief Adds a new point to the end of the signal data array. Same as
-/// AddPoint.
+/// \brief Adds a new point to the end of the signal data array
 ///
-void TRestRawSignal::AddDeposit(Short_t d) { AddPoint(d); }
+void TRestRawSignal::AddPoint(Double_t value) {
+    if (value > numeric_limits<Short_t>::max()) {
+        value = numeric_limits<Short_t>::max();
+    } else if (value < numeric_limits<Short_t>::min()) {
+        value = numeric_limits<Short_t>::min();
+    }
+    AddPoint((Short_t)value);
+}
 
 ///////////////////////////////////////////////
 /// \brief It overloads the operator [] so that we can retrieve a particular
@@ -568,15 +568,19 @@ void TRestRawSignal::GetDifferentialSignal(TRestRawSignal* diffSignal, Int_t sme
     if (smearPoints <= 0) smearPoints = 1;
     diffSignal->Initialize();
 
-    for (int i = 0; i < smearPoints; i++) diffSignal->AddPoint(0);
+    for (int i = 0; i < smearPoints; i++) {
+        diffSignal->AddPoint((Short_t)0);
+    }
 
     for (int i = smearPoints; i < this->GetNumberOfPoints() - smearPoints; i++) {
         Double_t value = 0.5 * (this->GetData(i + smearPoints) - GetData(i - smearPoints)) / smearPoints;
 
         diffSignal->AddPoint((Short_t)value);
     }
 
-    for (int i = GetNumberOfPoints() - smearPoints; i < GetNumberOfPoints(); i++) diffSignal->AddPoint(0);
+    for (int i = GetNumberOfPoints() - smearPoints; i < GetNumberOfPoints(); i++) {
+        diffSignal->AddPoint((Short_t)0);
+    }
 }
 
 ///////////////////////////////////////////////
@@ -588,15 +592,14 @@ void TRestRawSignal::GetDifferentialSignal(TRestRawSignal* diffSignal, Int_t sme
 /// as its standard deviation.
 ///
 void TRestRawSignal::GetWhiteNoiseSignal(TRestRawSignal* noiseSignal, Double_t noiseLevel) {
-    double* dd = new double();
-    uintptr_t seed = (uintptr_t)dd + (uintptr_t)this;
-    delete dd;
-    TRandom3* fRandom = new TRandom3(seed);
+    TRandom3 random(fSeed);
 
     for (int i = 0; i < GetNumberOfPoints(); i++) {
-        noiseSignal->AddPoint(this->GetData(i) + (Short_t)fRandom->Gaus(0, noiseLevel));
+        Double_t value = this->GetData(i) + random.Gaus(0, noiseLevel);
+        // do not cast as short so that there are no problems with overflows
+        // (https://github.com/rest-for-physics/rawlib/issues/113)
+        noiseSignal->AddPoint(value);
     }
-    delete fRandom;
 }
 
 ///////////////////////////////////////////////
@@ -628,7 +631,7 @@ void TRestRawSignal::GetSignalSmoothed(TRestRawSignal* smoothedSignal, Int_t ave
 ///////////////////////////////////////////////
 /// \brief It smoothes the existing signal and returns it in a vector of Float_t values
 ///
-/// \param averagingPoints It defines the number of neightbour consecutive
+/// \param averagingPoints It defines the number of neighbour consecutive
 /// points used to average the signal
 ///
 /// \param option If the option is set to "EXCLUDE OUTLIERS", points that are too far away from the median

src/TRestRawSignalAddNoiseProcess.cxx

@@ -142,7 +142,9 @@ void TRestRawSignalAddNoiseProcess::Initialize() {
 /// corresponding TRestGeant4AnalysisProcess section inside the RML.
 ///
 void TRestRawSignalAddNoiseProcess::LoadConfig(const string& configFilename, const string& name) {
-    if (LoadConfigFromFile(configFilename, name) == -1) LoadDefaultConfig();
+    if (LoadConfigFromFile(configFilename, name) == -1) {
+        LoadDefaultConfig();
+    }
 }
 
 ///////////////////////////////////////////////