From 91d5a04d3b760c83ebb0245c696675c691e16078 Mon Sep 17 00:00:00 2001
From: "P.-A. Loizeau" <p.-a.loizeau@gsi.de>
Date: Wed, 13 Apr 2022 16:28:30 +0200
Subject: [PATCH] Add data class + unpacking macro for the SIS18 timing events

---
 core/data/CMakeLists.txt                    |   1 +
 core/data/DataLinkDef.h                     |   3 +
 core/data/raw/AccDataSis18.cxx              | 169 ++++++++++++++
 core/data/raw/AccDataSis18.h                | 115 +++++++++
 macro/beamtime/mcbm2022/UnpackTimingSis18.C | 246 ++++++++++++++++++++
 5 files changed, 534 insertions(+)
 create mode 100644 core/data/raw/AccDataSis18.cxx
 create mode 100644 core/data/raw/AccDataSis18.h
 create mode 100644 macro/beamtime/mcbm2022/UnpackTimingSis18.C

diff --git a/core/data/CMakeLists.txt b/core/data/CMakeLists.txt
index cc3b4c4631..92828e5113 100644
--- a/core/data/CMakeLists.txt
+++ b/core/data/CMakeLists.txt
@@ -119,6 +119,7 @@ set(SRCS
   global/CbmTofTrack.cxx
   global/CbmTrackParam.cxx
 
+  raw/AccDataSis18.cxx
   raw/StsXyterMessage.cxx
   raw/gDpbMessv100.cxx
   raw/CriGet4Mess001.cxx
diff --git a/core/data/DataLinkDef.h b/core/data/DataLinkDef.h
index 223c19ebb0..ecec751d79 100644
--- a/core/data/DataLinkDef.h
+++ b/core/data/DataLinkDef.h
@@ -96,6 +96,9 @@
 #pragma link C++ class CbmVertex + ;
 #pragma link C++ class std::vector < CbmDigiEvent>;
 
+// --- data/raw
+#pragma link C++ class AccTimingEvent;
+#pragma link C++ class AccStatusTs;
 #pragma link C++ class stsxyter::Message;
 #pragma link C++ class gdpbv100::Message;
 #pragma link C++ class gdpbv100::FullMessage;
diff --git a/core/data/raw/AccDataSis18.cxx b/core/data/raw/AccDataSis18.cxx
new file mode 100644
index 0000000000..a5a62af581
--- /dev/null
+++ b/core/data/raw/AccDataSis18.cxx
@@ -0,0 +1,169 @@
+/* Copyright (C) 2022 Facility for Antiproton and Ion Research in Europe, Darmstadt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Pierre-Alain Loizeau [committer] */
+
+#include "AccDataSis18.h"
+
+#include <Logger.h>  // for LOG output
+
+#include <iomanip>
+#include <iostream>
+
+AccTimingEvent::AccTimingEvent(uint64_t ulPlannedUTCIn, uint64_t ulPlannedTAIIn, uint64_t ulRawEventIn,
+                               uint64_t ulRawParamsIn, uint32_t uRawTimingFlagsIn, uint64_t ulExecutedUTCIn,
+                               uint64_t ulExecutedTAIIn)
+  : fulPlannedUTC {ulPlannedUTCIn}
+  , fulPlannedTAI {ulPlannedTAIIn}
+  , fulRawEvent {ulRawEventIn}
+  , fulRawParams {ulRawParamsIn}
+  , fuRawTimingFlags {uRawTimingFlagsIn}
+  , fulExecutedUTC {ulExecutedUTCIn}
+  , fulExecutedTAI {ulExecutedTAIIn}
+{
+}
+
+AccTimingEvent::AccTimingEvent(std::string sLine, bool bVerbose)
+{
+  /*
+   * Format used to generate the timing files
+   * | std::cout << "Planned UTC: " << setw(20) << deadline.getUTC();
+   * | std::cout << " TAI: " << setw(20) <<deadline.getTAI();
+   * | std::cout << " Raw:" << std::hex << std::setfill('0')
+   * |           << " 0x" << std::setw(16) << id
+   * |           << " 0x" << std::setw(16) << param
+   * |           << " 0x" << std::setw(4) << flags
+   * |           << std::dec << std::setfill(' ');
+   * | std::cout << " exec UTC: " << setw(20) <<executed.getUTC();
+   * | std::cout << " TAI: " << setw(20) <<executed.getTAI();
+   * |
+   * | /// Dec  Hex  Name                  Meaning
+   * | ///
+   * | /// Spill limits
+   * | /// 46   2E   EVT_EXTR_START_SLOW   Start of extraction
+   * | /// 51   33   EVT_EXTR_END          End of extraction
+   * | /// 78   4E   EVT_EXTR_STOP_SLOW    End of slow extraction
+   * | ///
+   * | /// Cycle limits
+   * | /// 32   20   EVT_START_CYCLE       First Event in a cycle
+   * | /// 55   37   EVT_END_CYCLE         End of a cycle
+   * | uint32_t uEventNb = ((id >> 36) & 0xfff);
+   * | switch (uEventNb) {
+   * |   case 32:
+   * |     std::cout << " => EVT_START_CYCLE     ";
+   * |     break;
+   * |   case 55:
+   * |     std::cout << " => EVT_END_CYCLE       ";
+   * |     break;
+   * |   case 46:
+   * |     std::cout << " => EVT_EXTR_START_SLOW ";
+   * |     break;
+   * |   case 51:
+   * |     std::cout << " => EVT_EXTR_END        ";
+   * |     break;
+   * |   case 78:
+   * |     std::cout << " => EVT_EXTR_STOP_SLOW  ";
+   * |     break;
+   * | }
+   * | std::cout << tr_formatDate(deadline, pmode);
+   * | std::cout << tr_formatActionFlags(flags, executed - deadline, pmode);
+   * | std::cout << std::endl;
+   */
+
+  std::string sToken = "Planned UTC: ";
+  std::size_t posTok = sLine.find(sToken);
+  sLine              = sLine.substr(posTok + sToken.size());
+
+  sToken        = " TAI: ";
+  posTok        = sLine.find(sToken);
+  fulPlannedUTC = std::stoul(sLine.substr(0, posTok));
+  sLine         = sLine.substr(posTok + sToken.size());
+
+  sToken        = " Raw: 0x";
+  posTok        = sLine.find(sToken);
+  fulPlannedTAI = std::stoul(sLine.substr(0, posTok));
+  sLine         = sLine.substr(posTok + sToken.size());
+
+  sToken      = " 0x";
+  posTok      = sLine.find(sToken);
+  fulRawEvent = std::stoul(sLine.substr(0, posTok), 0, 16);
+  sLine       = sLine.substr(posTok + sToken.size());
+
+  sToken       = " 0x";
+  posTok       = sLine.find(sToken);
+  fulRawParams = std::stoul(sLine.substr(0, posTok), 0, 16);
+  sLine        = sLine.substr(posTok + sToken.size());
+
+  sToken           = " exec UTC: ";
+  posTok           = sLine.find(sToken);
+  fuRawTimingFlags = std::stoul(sLine.substr(0, posTok), 0, 16);
+  sLine            = sLine.substr(posTok + sToken.size());
+
+  sToken         = " TAI: ";
+  posTok         = sLine.find(sToken);
+  fulExecutedUTC = std::stoul(sLine.substr(0, posTok));
+  sLine          = sLine.substr(posTok + sToken.size());
+
+  sToken         = " => ";
+  posTok         = sLine.find(sToken);
+  fulExecutedTAI = std::stoul(sLine.substr(0, posTok));
+  sLine          = sLine.substr(posTok + sToken.size());
+
+  if (bVerbose) Print();
+}
+
+void AccTimingEvent::Print() const
+{
+  /* clang-format off */
+  LOG(info) << "Planned UTC: " << std::setw(20) << fulPlannedUTC
+            << " TAI: " << std::setw(20) << fulPlannedTAI
+            << " Raw:" << std::hex << std::setfill('0')
+            << " 0x" << std::setw(16) << fulRawEvent
+            << " 0x" << std::setw(16) << fulRawParams
+            << " 0x" << std::setw(4) << fuRawTimingFlags
+            << std::dec << std::setfill(' ')
+            << " exec UTC: " << std::setw(20) << fulExecutedUTC
+            << " TAI: " << std::setw(20) << fulExecutedTAI;
+  /* clang-format on */
+}
+//--------------------------------------------------------------------------------------------------------------------//
+
+bool AccStatusTs::IsSpillOnAtTime(uint64_t uTimeUtc)
+{
+  bool bSpillOn = IsSpillOnAtStart();
+
+  std::vector<AccTimingEvent>::iterator it = fvEventsDuringTS.begin();
+  while (it != fvEventsDuringTS.end() && (*it < uTimeUtc)) {
+
+    if (bSpillOn && it->IsExtractionEnd()) {
+      /// We start in the middle of an extraction spill
+      bSpillOn = false;
+    }
+    else if (!bSpillOn && it->IsExtractionStart()) {
+      /// We start in the middle of an extraction spill
+      bSpillOn = true;
+    }
+  }
+  return bSpillOn;
+}
+
+uint32_t AccStatusTs::GetSpillIdxAtTime(uint64_t uTimeUtc)
+{
+  bool bSpillOn      = IsSpillOnAtStart();
+  uint32_t uSpillIdx = fuSpillIndexAtStart;
+
+  std::vector<AccTimingEvent>::iterator it = fvEventsDuringTS.begin();
+  while (it != fvEventsDuringTS.end() && (*it < uTimeUtc)) {
+
+    if (bSpillOn && it->IsExtractionEnd()) {
+      /// We start in the middle of an extraction spill
+      bSpillOn = false;
+    }
+    else if (!bSpillOn && it->IsExtractionStart()) {
+      /// We start in the middle of an extraction spill
+      bSpillOn = true;
+      uSpillIdx++;
+    }
+  }
+  return uSpillIdx;
+}
+//--------------------------------------------------------------------------------------------------------------------//
diff --git a/core/data/raw/AccDataSis18.h b/core/data/raw/AccDataSis18.h
new file mode 100644
index 0000000000..538fdf0840
--- /dev/null
+++ b/core/data/raw/AccDataSis18.h
@@ -0,0 +1,115 @@
+/* Copyright (C) 2022 Facility for Antiproton and Ion Research in Europe, Darmstadt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Pierre-Alain Loizeau [committer] */
+
+#include <Rtypes.h>  // for THashConsistencyHolder, ClassDef
+
+#include <cstdint>
+#include <string>
+#include <vector>
+
+class AccTimingEvent {
+public:
+  AccTimingEvent() = default;
+
+  AccTimingEvent(uint64_t ulPlannedUTCIn, uint64_t ulPlannedTAIIn, uint64_t ulRawEventIn, uint64_t ulRawParamsIn,
+                 uint32_t uRawTimingFlagsIn, uint64_t ulExecutedUTCIn, uint64_t ulExecutedTAIIn);
+
+  AccTimingEvent(std::string sLine, bool bVerbose = false);
+
+  void Print() const;
+
+  uint32_t GetGroupId() const { return ((fulRawEvent >> kOffsetGroupId) & ((1ULL << kNbBitsGroupId) - 1)); }
+  uint32_t GetEventId() const { return ((fulRawEvent >> kOffsetEventId) & ((1ULL << kNbBitsEventId) - 1)); }
+
+  uint64_t GetTime() const { return fulPlannedUTC; }
+
+  inline bool IsCycleStart() const { return (kEventIdStartCycle == GetEventId()); }
+  inline bool IsCycleEnd() const { return (kEventIdEndCycle == GetEventId()); }
+  inline bool IsExtractionStart() const { return (kEventIdExtrStartSlow == GetEventId()); }
+  inline bool IsExtractionEnd() const { return (kEventIdExtrEnd == GetEventId()); }
+
+  inline bool operator<(const AccTimingEvent& rhs) { return fulPlannedUTC < rhs.fulPlannedUTC; }
+  inline bool operator<(const uint64_t& rhs) { return fulPlannedUTC < rhs; }
+
+  /// Needed for finding time position within a vector of AccTimingEvent with lower_bound/upper_bound
+  friend bool operator<(const uint64_t& lhs, const AccTimingEvent& rhs) { return lhs < rhs.fulPlannedUTC; }
+
+private:
+  /// Constants
+  /// --> Bit fields of the Raw event descriptor
+  /// ----> Field size
+  static const uint32_t kNbBitsFormatId   = 4;  // Content of field should always be 1!!
+  static const uint32_t kNbBitsGroupId    = 12;
+  static const uint32_t kNbBitsEventId    = 12;
+  static const uint32_t kNbBitsFlags      = 4;
+  static const uint32_t kNbBitsSequenceId = 12;
+  static const uint32_t kNbBitsBeamProcId = 14;
+  static const uint32_t kNbBitsReserved   = 6;
+  /// ----> Field offset
+  static const uint32_t kOffsetReserved   = 0;
+  static const uint32_t kOffsetBeamProcId = kOffsetReserved + kNbBitsReserved;
+  static const uint32_t kOffsetSequenceId = kOffsetBeamProcId + kNbBitsBeamProcId;
+  static const uint32_t kOffsetFlags      = kOffsetSequenceId + kNbBitsSequenceId;
+  static const uint32_t kOffsetEventId    = kOffsetFlags + kNbBitsFlags;
+  static const uint32_t kOffsetGroupId    = kOffsetEventId + kNbBitsEventId;
+  static const uint32_t kOffsetFormatId   = kOffsetGroupId + kNbBitsGroupId;
+  /// --> Event types
+  ///     Dec  Hex  Name                  Meaning
+  ///
+  ///     Spill limits
+  ///     46   2E   EVT_EXTR_START_SLOW   Start of extraction
+  ///     51   33   EVT_EXTR_END          End of extraction
+  ///     78   4E   EVT_EXTR_STOP_SLOW    End of slow extraction
+  ///
+  ///     Cycle limits
+  ///     32   20   EVT_START_CYCLE       First Event in a cycle
+  ///     55   37   EVT_END_CYCLE         End of a cycle
+  static const uint32_t kEventIdStartCycle    = 32;
+  static const uint32_t kEventIdExtrStartSlow = 46;
+  static const uint32_t kEventIdExtrEnd       = 51;
+  static const uint32_t kEventIdEndCycle      = 55;
+  static const uint32_t kEventIdExtrStopSlow  = 78;
+
+  /// Fields
+  uint64_t fulPlannedUTC    = 0;
+  uint64_t fulPlannedTAI    = 0;
+  uint64_t fulRawEvent      = 0;
+  uint64_t fulRawParams     = 0;
+  uint32_t fuRawTimingFlags = 0;
+  uint64_t fulExecutedUTC   = 0;
+  uint64_t fulExecutedTAI   = 0;
+
+  ClassDef(AccTimingEvent, 1);
+};
+
+//--------------------------------------------------------------------------------------------------------------------//
+
+class AccStatusTs {
+public:
+  AccStatusTs() = default;
+
+  AccStatusTs(uint32_t uSpillIdx, AccTimingEvent lastEvtBefTs)
+    : fuSpillIndexAtStart(uSpillIdx)
+    , fLastEvtBeforeTs(lastEvtBefTs)
+  {
+  }
+
+  void SetLastEvtBefTs(AccTimingEvent lastEvtBefTs) { fLastEvtBeforeTs = lastEvtBefTs; }
+
+  /// True when we start within a spill cycle (not in short interval between cycle end and cycle start))
+  inline bool IsCycleOnAtStart() const { return !(fLastEvtBeforeTs.IsCycleEnd()); }
+  /// True when we start in the middle of an extraction spill
+  inline bool IsSpillOnAtStart() const { return fLastEvtBeforeTs.IsExtractionStart(); }
+
+  bool IsSpillOnAtTime(uint64_t uTimeUtc);
+  uint32_t GetSpillIdxAtTime(uint64_t uTimeUtc);
+
+  /// Members
+  uint32_t fuSpillIndexAtStart                 = 0;
+  AccTimingEvent fLastEvtBeforeTs              = {};
+  std::vector<AccTimingEvent> fvEventsDuringTS = {};
+
+  ClassDef(AccStatusTs, 1);
+};
+//--------------------------------------------------------------------------------------------------------------------//
diff --git a/macro/beamtime/mcbm2022/UnpackTimingSis18.C b/macro/beamtime/mcbm2022/UnpackTimingSis18.C
new file mode 100644
index 0000000000..81d74e4eee
--- /dev/null
+++ b/macro/beamtime/mcbm2022/UnpackTimingSis18.C
@@ -0,0 +1,246 @@
+/* Copyright (C) 2022 Facility for Antiproton and Ion Research in Europe, Darmstadt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Pierre-Alain Loizeau [committer] */
+
+void UnpackTimingSis18(std::string sDigiFileName,
+                       std::string sTimingFilesFolder = "/lustre/cbm/users/ploizeau/mcbm2022", bool bVerbose = false)
+{
+  std::vector<std::string> vTimingFileNames = {"sis18_events_2022_03_24.txt", "sis18_events_2022_03_29.txt",
+                                               "sis18_events_2022_03_30.txt", "sis18_events_2022_03_31.txt",
+                                               "sis18_events_2022_04_01.txt", "sis18_events_2022_04_02.txt",
+                                               "sis18_events_2022_04_03.txt", "sis18_events_2022_04_04.txt"};
+
+  /// Test the decoder
+  /*
+  std::string testStr = "Planned UTC:  1648763972022760000 TAI:  1648764009022760000 Raw: 0x112c0338003005c0 "
+                        "0x0000140000000000 0x0000 exec UTC:  1648763972022760000 TAI:  1648764009022760000 "
+                        "=> EVT_EXTR_END        2022-03-31 22:00:09.022760000";
+  std::cout << testStr << std::endl;
+  AccTimingEvent tester(testStr);
+  */
+
+  std::vector<AccTimingEvent> vAllEvtsBuff = {};
+
+  /// Loop on Timing files
+  for (uint32_t uFile = 0; uFile < vTimingFileNames.size(); ++uFile) {
+    /// Open File
+    std::string sFullPath = sTimingFilesFolder + "/" + vTimingFileNames[uFile];
+    std::ifstream ifsRaw(sFullPath);
+
+    if (ifsRaw.is_open()) {
+      /// Read it line by line
+      std::string sLine;
+      while (std::getline(ifsRaw, sLine)) {
+        /// Convert each line to an event, if not empty
+        if ("\n" != sLine) {
+          if (bVerbose) std::cout << sLine << std::endl;
+          AccTimingEvent newEvent(sLine);
+          vAllEvtsBuff.push_back(newEvent);
+        }
+      }
+      std::cout << "File " << sFullPath << " done" << std::endl;
+    }
+    else {
+      std::cout << "File " << sFullPath << " could not be open!" << std::endl;
+      return;
+    }
+  }
+  std::cout << "Total events in all recorded files: " << vAllEvtsBuff.size() << std::endl;
+
+  /// Open the digi file in update mode
+  TFile* pFile = new TFile(sDigiFileName.data(), "UPDATE");
+  if (kTRUE == pFile->IsZombie()) {
+    std::cout << "Digi file " << sDigiFileName << " could not be open!" << std::endl;
+    return;
+  }
+
+  /// Access the tree in the file
+  TTree* pTree = dynamic_cast<TTree*>(pFile->Get("cbmsim"));
+  if (nullptr == pTree) {
+    std::cout << "Could not find the cbm tree in the digi file" << std::endl;
+    pFile->Close();
+    return;
+  }
+  uint32_t uNbTs = pTree->GetEntries();
+  if (0 == uNbTs) {
+    std::cout << "No entries in the Tree found in the digi file!" << std::endl;
+    pFile->Close();
+    return;
+  }
+
+  /// Access the header of the last TS
+  CbmTsEventHeader* pEventHeader = nullptr;
+  pTree->SetBranchAddress("EventHeader.", &pEventHeader);
+  pTree->GetEntry(uNbTs - 1);
+  if (nullptr == pEventHeader) {
+    std::cout << "Could not find the last TS Event header in the digi file" << std::endl;
+    pFile->Close();
+    return;
+  }
+  /// Get last start time and check that it is later than the first event in the buffer
+  uint64_t ulLastTsStartUtc = pEventHeader->GetTsStartTime();
+  if (ulLastTsStartUtc < vAllEvtsBuff[0].GetTime()) {
+    std::cout << "Last TS start time in this file: " << ulLastTsStartUtc << std::endl;
+    std::cout << "Earlier than the first recorded accelerator Event: " << vAllEvtsBuff[0].GetTime() << std::endl;
+    std::cout << "=> No eccelerator data available for this digi file!" << std::endl;
+    return;
+  }
+
+  /// Access the header of the first TS
+  pTree->GetEntry(0);
+  if (nullptr == pEventHeader) {
+    std::cout << "Could not find the first TS Event header in the digi file" << std::endl;
+    pFile->Close();
+    return;
+  }
+
+  /// Get first start time and search for the closest past event in the buffer
+  uint64_t ulFirstTsStartUtc = pEventHeader->GetTsStartTime();
+  std::cout << "First TS start time: " << ulFirstTsStartUtc << std::endl;
+
+  std::vector<AccTimingEvent>::iterator itEvt =
+    std::upper_bound(vAllEvtsBuff.begin(), vAllEvtsBuff.end(), ulFirstTsStartUtc);
+  if (vAllEvtsBuff.end() == itEvt) {
+    std::cout << "No Events after the start of the first TS => no accelerator data available for this run!"
+              << std::endl;
+    pFile->Close();
+    return;
+  }
+  /*
+  if (itEvt->IsCycleEnd()) {
+    if (itEvt != vAllEvtsBuff.begin()) {
+      itEvt--; // not at end of vector so rewind to previous item
+    }
+    else {
+      std::cout << "No Events before the start of the first TS & next event is the end of a cycle "
+                << "=> unknown state, not handled in current implementation!"
+                << std::endl;
+      pFile->Close();
+      return;
+    }
+  }
+  */
+  std::cout << "Found an event before the first TS start time:" << std::endl;
+  itEvt->Print();
+
+  /// Create a list of spills within the digi file
+  bool bExtractionOn = false;
+  bool bCycleOn      = false;
+  uint32_t uSpillIdx = 0;
+  uint32_t uCycleIdx = 0;
+  if (itEvt->IsExtractionEnd()) {
+    /// We start in the middle of an extraction spill
+    bExtractionOn = true;
+    std::cout << "Spill ongoing when starting at " << ulFirstTsStartUtc << " (" << uSpillIdx << ") " << std::endl;
+  }
+  if (!itEvt->IsCycleStart()) {
+    /// We start within a spill cycle (not in short interval between cycle end and cycle start))
+    bCycleOn = true;
+  }
+  while ((vAllEvtsBuff.end() != itEvt) && itEvt->GetTime() <= ulLastTsStartUtc) {
+    if (itEvt->IsCycleStart()) {
+      //
+      uCycleIdx++;
+      bCycleOn = true;
+    }
+    else if (itEvt->IsCycleEnd()) {
+      //
+      bCycleOn = false;
+    }
+    else if (itEvt->IsExtractionStart()) {
+      //
+      uSpillIdx++;
+      bExtractionOn = true;
+      std::cout << "Spill starting at              " << itEvt->GetTime() << " (" << uSpillIdx << ") " << std::endl;
+    }
+    else if (itEvt->IsExtractionEnd()) {
+      //
+      bExtractionOn = false;
+      std::cout << "Spill end at                   " << itEvt->GetTime() << std::endl;
+    }
+
+    itEvt++;
+  }
+  std::cout << "Number of spill found for this digi file: " << (1 + uSpillIdx) << std::endl;
+
+  /*
+   * TODO: decide what is saved in the digi file
+   * => Should we try looking for the first spill start before the beginning of the first TS?
+   * => Should we save only the start and stop of the spills? or also the index within the file?
+   * => Should we save the spill times in UTC or in ref to first TS time?
+   */
+  /*
+   * FIXME: not compatible with usage of timeslice overlap as relying on "end of TS = start of next TS"
+   * => Need addin the limits of the timeslice core and timeslice overlap in the TsEventHeader class (from TsMetaData)
+   * => Need to add to the RecoUnpack a determination of these parameters from the boundaries of the microslices and
+   *    microslices numbers on component
+   */
+  if (uNbTs < 2) {
+    std::cout << "Only 1 TS in the digi file, cannot determine the TS duration" << std::endl;
+    pFile->Close();
+    return;
+  }
+  pTree->GetEntry(1);
+  uint64_t ulTsDuration = pEventHeader->GetTsStartTime() - ulFirstTsStartUtc;
+
+  /// Create for each TS a list of events, including the last event before its start and all event during its duration
+  AccStatusTs statusTs;
+  TBranch* pBranch = pTree->Branch("AccStatusTs", &statusTs);
+
+  std::vector<AccTimingEvent>::iterator itFirstEventPrevTs = vAllEvtsBuff.begin();
+  std::vector<AccTimingEvent>::iterator itLastEventPrevTs  = vAllEvtsBuff.begin();
+  uSpillIdx                                                = 0;
+  for (uint32_t uTs = 0; uTs < uNbTs; ++uTs) {
+    pTree->GetEntry(uTs);
+
+    /// Find the first event after the start of this TS
+    uint64_t ulTsStartUtc = pEventHeader->GetTsStartTime();
+    std::vector<AccTimingEvent>::iterator itEvtPrev =
+      std::upper_bound(itFirstEventPrevTs, vAllEvtsBuff.end(), ulTsStartUtc);
+
+    /// Loop on events between last TS and the current one to count potentially not recorded spills
+    for (std::vector<AccTimingEvent>::iterator it = itLastEventPrevTs; it != itEvtPrev; ++it) {
+      if (itEvt->IsExtractionStart()) {
+        //
+        uSpillIdx++;
+      }
+    }
+
+    if (vAllEvtsBuff.begin() == itEvtPrev) {
+      std::cout << "No Events before the start of the first TS "
+                << "=> not possible to do a full list of spill status for this digi file!" << std::endl;
+      pFile->Close();
+      return;
+    }
+    itEvtPrev--;
+
+    /// Save iterator to speed up the search for the next event
+    itFirstEventPrevTs = itEvtPrev;
+
+    /// Find the first event after the end of this TS
+    std::vector<AccTimingEvent>::iterator itEvtNext =
+      std::upper_bound(itFirstEventPrevTs, vAllEvtsBuff.end(), ulTsStartUtc + ulTsDuration);
+
+    statusTs.SetLastEvtBefTs(*itEvtPrev);
+    statusTs.fuSpillIndexAtStart             = uSpillIdx;
+    std::vector<AccTimingEvent>::iterator it = itEvtPrev;
+    it++;
+    for (; it != itEvtNext; ++it) {
+      statusTs.fvEventsDuringTS.push_back(*it);
+      if (it->IsExtractionStart()) {
+        //
+        uSpillIdx++;
+        std::cout << "Spill starting at              " << it->GetTime() << " (" << uSpillIdx << ") " << std::endl;
+      }
+    }
+    pBranch->Fill();
+
+    /// Save iterator to allow detection of unrecorded spills in case of missing TS
+    itLastEventPrevTs = itEvtNext;
+  }
+  /// Save only the new version of the tree which includes the new branch
+  pTree->Write("", TObject::kOverwrite);
+
+
+  pFile->Close();
+}
-- 
GitLab