Merge pull request #45 from rest-for-physics/mariajmz_raw_noise

TRestDetectorSignalToRawSignalProcess adding optional noise process
rest-for-physics · May 6, 2024 · a577f65 · a577f65
2 parents da19e6b + 7df364d
commit a577f65
Show file tree

Hide file tree

Showing 2 changed files with 29 additions and 6 deletions.
diff --git a/inc/TRestDetectorSignalToRawSignalProcess.h b/inc/TRestDetectorSignalToRawSignalProcess.h
@@ -83,6 +83,9 @@ class TRestDetectorSignalToRawSignalProcess : public TRestEventProcess {
     /// avoid artifacts in the signal (e.g. signals not getting cut when they should)
     Double_t fShapingTime = 0.0;  // us
 
+    // Noise level
+    Double_t fNoiseLevel = 0.0;
+
    public:
     inline Double_t GetSampling() const { return fSampling; }
 
@@ -118,6 +121,7 @@ class TRestDetectorSignalToRawSignalProcess : public TRestEventProcess {
         Double_t shapingTime = 0.0;
         Double_t calibrationGain = 100;
         Double_t calibrationOffset = 0;
+        Double_t noiseLevel = 0.0;
         TVector2 calibrationEnergy = {0, 0};
         TVector2 calibrationRange = {0, 0};
     };
@@ -137,19 +141,17 @@ class TRestDetectorSignalToRawSignalProcess : public TRestEventProcess {
     /// Returns the name of this process
     const char* GetProcessName() const override { return "signalToRawSignal"; }
 
-    // Constructor
     TRestDetectorSignalToRawSignalProcess();
 
-    TRestDetectorSignalToRawSignalProcess(const char* configFilename);
+    explicit TRestDetectorSignalToRawSignalProcess(const char* configFilename);
 
-    // Destructor
-    ~TRestDetectorSignalToRawSignalProcess();
+    ~TRestDetectorSignalToRawSignalProcess() override;
 
    private:
     std::map<std::string, Parameters> fParametersMap;
     std::set<std::string> fReadoutTypes;
 
-    ClassDefOverride(TRestDetectorSignalToRawSignalProcess, 7);
+    ClassDefOverride(TRestDetectorSignalToRawSignalProcess, 8);
 };
 
 #endif
diff --git a/src/TRestDetectorSignalToRawSignalProcess.cxx b/src/TRestDetectorSignalToRawSignalProcess.cxx
@@ -381,6 +381,7 @@ TRestEvent* TRestDetectorSignalToRawSignalProcess::ProcessEvent(TRestEvent* inpu
             type = "";
         }
 
+        double noiseLevel = fParametersMap.at(type).noiseLevel;
         double sampling = fParametersMap.at(type).sampling;
         double shapingTime = fParametersMap.at(type).shapingTime;
         double calibrationGain = fParametersMap.at(type).calibrationGain;
@@ -420,7 +421,7 @@ TRestEvent* TRestDetectorSignalToRawSignalProcess::ProcessEvent(TRestEvent* inpu
 
             if (t > timeStart && t < timeEnd) {
                 // convert physical time (in us) to timeBin
-                Int_t timeBin = (Int_t)round((t - timeStart) / sampling);
+                auto timeBin = (Int_t)round((t - timeStart) / sampling);
 
                 if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Warning) {
                     if (timeBin < 0 || timeBin > fNPoints) {
@@ -434,6 +435,13 @@ TRestEvent* TRestDetectorSignalToRawSignalProcess::ProcessEvent(TRestEvent* inpu
             }
         }
 
+        // Noise before shaping
+        if (noiseLevel > 0) {
+            for (int i = 0; i < fNPoints; i++) {
+                data[i] += gRandom->Gaus(0, noiseLevel);
+            }
+        }
+
         if (shapingTime > 0) {
             const auto sinShaper = [](Double_t t) -> Double_t {
                 if (t <= 0) {
@@ -466,6 +474,13 @@ TRestEvent* TRestDetectorSignalToRawSignalProcess::ProcessEvent(TRestEvent* inpu
                 }
             }
             data = dataAfterShaping;
+
+            // Noise after shaping
+            if (noiseLevel > 0) {
+                for (int i = 0; i < fNPoints; i++) {
+                    data[i] += gRandom->Gaus(0, noiseLevel);
+                }
+            }
         }
 
         TRestRawSignal rawSignal;
@@ -547,6 +562,7 @@ void TRestDetectorSignalToRawSignalProcess::InitFromConfigFile() {
             Get2DVectorParameterWithUnits("calibrationEnergy" + typeCamelCase, parameters.calibrationEnergy);
         parameters.calibrationRange =
             Get2DVectorParameterWithUnits("calibrationRange" + typeCamelCase, parameters.calibrationRange);
+        parameters.noiseLevel = GetDblParameterWithUnits("noiseLevel" + typeCamelCase, parameters.noiseLevel);
 
         const bool isLinearCalibration =
             (parameters.calibrationEnergy.Mod() != 0 && parameters.calibrationRange.Mod() != 0);
@@ -599,6 +615,7 @@ void TRestDetectorSignalToRawSignalProcess::InitFromConfigFile() {
     // load default parameters (for backward compatibility)
     fSampling = fParametersMap.at(defaultType).sampling;
     fShapingTime = fParametersMap.at(defaultType).shapingTime;
+    fNoiseLevel = fParametersMap.at(defaultType).noiseLevel;
     fCalibrationGain = fParametersMap.at(defaultType).calibrationGain;
     fCalibrationOffset = fParametersMap.at(defaultType).calibrationOffset;
     fCalibrationEnergy = fParametersMap.at(defaultType).calibrationEnergy;
@@ -678,6 +695,10 @@ void TRestDetectorSignalToRawSignalProcess::PrintMetadata() {
         if (shapingTime > 0) {
             RESTMetadata << "Shaping time: " << shapingTime * 1000 << " ns" << RESTendl;
         }
+        const double noiseLevel = fParametersMap.at(readoutType).noiseLevel;
+        if (noiseLevel > 0) {
+            RESTMetadata << "Noise Level: " << noiseLevel << RESTendl;
+        }
 
         if (IsLinearCalibration()) {
             RESTMetadata << "Calibration energies: (" << fParametersMap.at(readoutType).calibrationEnergy.X()