Adding time correction in dataset generation

.github/workflows/validation.yml

@@ -149,7 +149,7 @@ jobs:
           python3 -c """import ROOT;import REST;ROOT.TRestRun();ROOT.TRestGeant4Event()"""
 
   metadata:
-    name: "Metadata"
+    name: "Datasets"
     runs-on: ubuntu-latest
     container:
       image: ghcr.io/lobis/root-geant4-garfield:rest-for-physics
@@ -169,28 +169,14 @@ jobs:
         with:
           key: ${{ env.BRANCH_NAME }}-${{ github.sha }}
           path: ${{ env.REST_PATH }}
-      # Not used in gitlab pipeline, checks $CI_SERVER_HOST
-      - name: Load Gas
+      - name: Basic tests
         run: |
           source ${{ env.REST_PATH }}/thisREST.sh
-          cd framework/pipeline/metadata/gas/
-          restRoot -b -q LoadGasFromServerAndValidate.C
-      # Not working, not used in gitlab pipeline
-      #- name: Generate Gas
-      #  run: |
-      #    source ${{ env.REST_PATH }}/thisREST.sh
-      #    source $(root-config --bindir)/thisroot.sh
-      #    cd framework/pipeline/metadata/gas/
-      #    restRoot -b -q GenerateDummyGas.C
-      - name: Basic Readout
-        run: |
-          source ${{ env.REST_PATH }}/thisREST.sh
-          cd framework/projects/basic-examples/
-          git submodule update --init .
-          cd basic-readouts
-          git submodule update --init .
-          restRoot -b -q GenerateReadouts.C'("basic.root")'
-          restRoot -b -q BasicValidation.C'("basic.root", "pixelDecoding")'
+          cd framework/pipeline/datasets/
+          wget https://github.com/rest-for-physics/rest-school/raw/master/data/R11520_00000_RawToTrack_Calibration_30min_jgalan_2.3.15.root
+          wget https://github.com/rest-for-physics/rest-school/raw/master/data/R11521_00000_RawToTrack_Background_23hr_jgalan_2.3.15.root
+          wget https://github.com/rest-for-physics/rest-school/raw/master/data/R11523_00000_RawToTrack_Background_20hr_jgalan_2.3.15.root
+          restRoot -b -q DataSetTests.C
 
   pandax-iii:
     name: "PandaX-III"

pipeline/datasets/DataSetTests.C

@@ -0,0 +1,26 @@
+// We run a couple validations to see if everything in the dataset is ok
+
+Int_t DataSetTests() {
+    TRestDataSet d("dataset.rml");
+    d.GenerateDataSet();
+
+    std::cout << "RunTimeStart: " << d.GetFilterStartTime() << std::endl;
+    std::cout << "RunTimeEnd: " << d.GetFilterEndTime() << std::endl;
+
+    if (d.GetFilterStartTime() != "2021/04/28 00:00") {
+        return 1;
+    }
+
+    if (d.GetTree()->GetEntries() != 3816) {
+        return 2;
+    }
+
+    if (d.GetFileSelection().size() != 1) {
+        return 3;
+    }
+
+    if (d.GetDataFrame().GetColumnNames().size() != 13) {
+        return 4;
+    }
+    return 0;
+}

pipeline/datasets/dataset.rml

@@ -0,0 +1,32 @@
+<datasets>
+
+<TRestDataSet name="Calibrations">
+   <parameter name="filterStartTime" value = "2021/04/28 00:00" />
+   <parameter name="filterEndTime" value = "2021/06/28 23:59" />
+   <parameter name="filePattern" value="*RawToTrack*.root"/>
+   <parameter name="timeCorrection" value = "true" />
+   <parameter name="columnTimeCorrection" value = "rateAna_MeanRate_InHz" />
+   <parameter name="thresholdTimeCorrection" value = "15" />
+   <filter metadata="TRestRun::fRunTag" contains="Calibration" />
+
+   // Will add to the final tree only the specific observables
+   <observables list="sAna_NumberOfGoodSignals,tckAna_xMean,tckAna_yMean,sAna_ThresholdIntegral,sAna_BaseLineMean,sAna_RiseTimeAvg" />
+   // Will add all the observables from the process `rateAna`
+   <processObservables list="rateAna" />
+</TRestDataSet>
+
+<TRestDataSet name="Backgrounds">
+   <parameter name="startTime" value = "2022/04/30 00:00" />
+   <parameter name="endTime" value = "2022/04/31 23:25" />
+   <parameter name="filePattern" value="*RawToTrack*.root"/>
+   <parameter name="timeCorrection" value = "true" />
+   <parameter name="columnTimeCorrection" value = "rateAna_MeanRate_InHz" />
+   <parameter name="thresholdTimeCorrection" value = "15" />
+   <filter metadata="TRestRun::fRunTag" contains="Background" />
+   // Will add to the final tree only the specific observables
+   <observables list="sAna_NumberOfGoodSignals,tckAna_xMean,tckAna_yMean,sAna_ThresholdIntegral,sAna_BaseLineMean,sAna_RiseTimeAvg" />
+   // Will add all the observables from the process `rateAna`
+   <processObservables list="rateAna" />
+</TRestDataSet>
+
+</datasets>

source/framework/core/inc/TRestDataSet.h

@@ -99,6 +99,15 @@ class TRestDataSet : public TRestMetadata {
     /// It keeps track if the generated dataset is a pure dataset or a merged one
     Bool_t fMergedDataset = false;  //<
 
+    /// Keeps track if the time correction algorithm will be applied
+    Bool_t fTimeCorrection = false;  //<
+
+    /// Column of the parameter over which we have to apply the threshold below for time correction
+    std::string fColumnTimeCorrection = "";  //<
+
+    /// Threshold below which to apply time correction
+    Double_t fThresholdTimeCorrection = 0;  //<
+
     /// The list of dataset files imported
     std::vector<std::string> fImportedFiles;  //<
 
@@ -119,6 +128,8 @@ class TRestDataSet : public TRestMetadata {
 
     void InitFromConfigFile() override;
 
+    Double_t GetRunDuration(TRestRun& r);
+
    protected:
     virtual std::vector<std::string> FileSelection();
 
@@ -180,6 +191,7 @@ class TRestDataSet : public TRestMetadata {
     inline auto GetAddedColumns() const { return fColumnNameExpressions; }
     inline auto GetCut() const { return fCut; }
     inline auto IsMergedDataSet() const { return fMergedDataset; }
+    inline auto IsTimeCorrected() const { return fTimeCorrection; }
 
     inline void SetObservablesList(const std::vector<std::string>& obsList) { fObservablesList = obsList; }
     inline void SetFilePattern(const std::string& pattern) { fFilePattern = pattern; }

source/framework/core/src/TRestDataSet.cxx

@@ -469,12 +469,11 @@ std::vector<std::string> TRestDataSet::FileSelection() {
 
             if (properties.strategy == "last") properties.value = value;
         }
-
         if (run.GetStartTimestamp() < fStartTime) fStartTime = run.GetStartTimestamp();
 
         if (run.GetEndTimestamp() > fEndTime) fEndTime = run.GetEndTimestamp();
 
-        fTotalDuration += run.GetEndTimestamp() - run.GetStartTimestamp();
+        fTotalDuration += GetRunDuration(run);
         fFileSelection.push_back(file);
     }
     std::cout << std::endl;
@@ -564,6 +563,12 @@ void TRestDataSet::PrintMetadata() {
     RESTMetadata << " - Accumulated run time (seconds) : " << fTotalDuration << RESTendl;
     RESTMetadata << " - Accumulated run time (hours) : " << fTotalDuration / 3600. << RESTendl;
     RESTMetadata << " - Accumulated run time (days) : " << fTotalDuration / 3600. / 24. << RESTendl;
+    RESTMetadata << "  " << RESTendl;
+    RESTMetadata << " - Time correction activated: " << fTimeCorrection << RESTendl;
+    if (fTimeCorrection) {
+        RESTMetadata << " - Time correction applied to data in column: " << fColumnTimeCorrection << RESTendl;
+        RESTMetadata << " - With a threshold value of: " << fThresholdTimeCorrection << RESTendl;
+    }
 
     RESTMetadata << "  " << RESTendl;
 
@@ -640,6 +645,10 @@ void TRestDataSet::PrintMetadata() {
         for (const auto& fn : fImportedFiles) RESTMetadata << " - " << fn << RESTendl;
     }
 
+    if (fTimeCorrection) {
+        RESTMetadata << "The combined dataset time correction analysis is activated." << RESTendl;
+    }
+
     RESTMetadata << " " << RESTendl;
     if (fMT)
         RESTMetadata << " - Multithreading was enabled" << RESTendl;
@@ -1084,3 +1093,66 @@ void TRestDataSet::Import(std::vector<std::string> fileNames) {
 
     fQuantity.clear();
 }
+
+///////////////////////////////////////////////
+/// \brief This function calculates a corrected time of a given Run,
+/// corresponding to the time where a certain variable has remained above a
+/// given threshold (for example, the time for which the detection rate is above a certain value)
+Double_t TRestDataSet::GetRunDuration(TRestRun& r) {
+    Double_t runTime = 0;
+    Double_t corrected_time = 0;
+    if (fTimeCorrection) {
+        Int_t max_events_below_threshold = 5;
+        Int_t events_below_threshold =
+            0;                      // Number of events in a row (except recovery) below the desired threshold
+        Int_t recovery_events = 0;  // Events above the threshold after a series of bad events
+        Int_t first_bad_event = 0;
+        Int_t last_bad_event = 0;
+        Int_t t1 = 0;
+        Int_t t2 = 0;
+        Double_t val = 0.0;
+        if (r.GetAnalysisTree()->GetObservableType(fColumnTimeCorrection) != "double") {
+            RESTError << fColumnTimeCorrection << "is not a double" << RESTendl;
+            exit(1);
+        }
+        for (int ev = 0; ev < r.GetEntries(); ev++) {
+            r.GetEntry(ev);
+            val = r.GetAnalysisTree()->GetObservableValue<double>(fColumnTimeCorrection);
+
+            if (val < fThresholdTimeCorrection) {
+                recovery_events = 0;
+                if (events_below_threshold == 0) {
+                    first_bad_event = ev;
+                }
+                last_bad_event = ev;
+                events_below_threshold++;
+
+            } else {
+                if (events_below_threshold <= max_events_below_threshold) {
+                    events_below_threshold = 0;
+                } else {
+                    recovery_events++;
+                    if (recovery_events > max_events_below_threshold) {
+                        r.GetEntry(first_bad_event);
+                        t1 = r.GetAnalysisTree()->GetTimeStamp();
+                        r.GetEntry(last_bad_event);
+                        t2 = r.GetAnalysisTree()->GetTimeStamp();
+                        corrected_time += t2 - t1;
+                        events_below_threshold = 0;
+                    }
+                }
+            }
+            if ((ev == r.GetEntries() - max_events_below_threshold) &&
+                (events_below_threshold > max_events_below_threshold)) {
+                last_bad_event = r.GetEntries() - 1;
+                r.GetEntry(first_bad_event);
+                t1 = r.GetAnalysisTree()->GetTimeStamp();
+                r.GetEntry(last_bad_event);
+                t2 = r.GetAnalysisTree()->GetTimeStamp();
+                corrected_time += t2 - t1;
+            }
+        }
+    }
+    runTime = r.GetEndTimestamp() - r.GetStartTimestamp() - corrected_time;
+    return runTime;
+}