Adding time correction in dataset generation

.github/workflows/validation.yml

@@ -182,6 +182,25 @@ jobs:
       #    source $(root-config --bindir)/thisroot.sh
       #    cd framework/pipeline/metadata/gas/
       #    restRoot -b -q GenerateDummyGas.C
+      - name: Generate Readout
+        run: |
+          source ${{ env.REST_PATH }}/thisREST.sh
+          cd framework/pipeline/metadata/readout/
+          restManager --c generateReadout.rml --o readout.root
+          restRoot -b -q PrintReadout.C'("readout.root")' > /dev/null
+      # We need to introduce basic validation here
+      #    - diff validation.txt print.txt
+
+      - name: DataSet Tests
+        run: |
+          source ${{ env.REST_PATH }}/thisREST.sh
+          cd framework/pipeline/datasets/
+          wget https://github.com/rest-for-physics/rest-school/blob/master/data/R11520_00000_RawToTrack_Calibration_30min_jgalan_2.3.15.root
+          wget https://github.com/rest-for-physics/rest-school/blob/master/data/R11521_00000_RawToTrack_Background_23hr_jgalan_2.3.15.root
+          wget https://github.com/rest-for-physics/rest-school/blob/master/data/R11523_00000_RawToTrack_Background_20hr_jgalan_2.3.15.root
+          restRoot -b -q DataSetTests.C'
+
+
       - name: Basic Readout
         run: |
           source ${{ env.REST_PATH }}/thisREST.sh

source/framework/core/inc/TRestDataSet.h

@@ -98,6 +98,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;  //<
 
@@ -109,6 +118,8 @@ class TRestDataSet : public TRestMetadata {
 
     void InitFromConfigFile() override;
 
+    Double_t GetRunDuration(TRestRun& r);
+
    protected:
     virtual std::vector<std::string> FileSelection();
 
@@ -159,6 +170,7 @@ class TRestDataSet : public TRestMetadata {
     inline auto GetQuantity() const { return fQuantity; }
     inline auto GetCut() const { return fCut; }
     inline auto IsMergedDataSet() const { return fMergedDataset; }
+    inline auto IsTimeCorrected() const { return fTimeCorrection; }
 
     inline void SetFilePattern(const std::string& pattern) { fFilePattern = pattern; }
 
@@ -180,6 +192,6 @@ class TRestDataSet : public TRestMetadata {
     TRestDataSet(const char* cfgFileName, const std::string& name = "");
     ~TRestDataSet();
 
-    ClassDefOverride(TRestDataSet, 3);
+    ClassDefOverride(TRestDataSet, 4);
 };
 #endif

source/framework/core/src/TRestDataSet.cxx

@@ -425,12 +425,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);
     }
     RESTInfo << RESTendl;
@@ -519,6 +518,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;
 
@@ -588,6 +593,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;
 }
 
@@ -917,3 +926,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()->GetObservableValue<double>("rateAna_SecondsFromStart");
+                        r.GetEntry(last_bad_event);
+                        t2 = r.GetAnalysisTree()->GetObservableValue<double>("rateAna_SecondsFromStart");
+                        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()->GetObservableValue<double>("rateAna_SecondsFromStart");
+                r.GetEntry(last_bad_event);
+                t2 = r.GetAnalysisTree()->GetObservableValue<double>("rateAna_SecondsFromStart");
+                corrected_time += t2 - t1;
+            }
+        }
+    }
+    runTime = r.GetEndTimestamp() - r.GetStartTimestamp() - corrected_time;
+    return runTime;
+}