From 9bd39fe17eeff63965338f6df1a982232c964f03 Mon Sep 17 00:00:00 2001
From: Volker Friese <v.friese@gsi.de>
Date: Tue, 23 Nov 2021 18:22:50 +0000
Subject: [PATCH] Integration of algo::TimeClusterTrigger in the run/task
 framework. Refs #2270.

---
 core/base/CbmDigiManager.cxx      |   4 +-
 macro/reco/reco_digi.C            |  32 ++++--
 reco/tasks/CMakeLists.txt         |   2 +
 reco/tasks/CbmRecoTasksLinkDef.h  |   1 +
 reco/tasks/CbmTaskBuildEvents.cxx |  20 ++--
 reco/tasks/CbmTaskBuildEvents.h   |  11 ++-
 reco/tasks/CbmTaskTriggerDigi.cxx | 155 ++++++++++++++++++++++++++++++
 reco/tasks/CbmTaskTriggerDigi.h   |  95 ++++++++++++++++++
 8 files changed, 298 insertions(+), 22 deletions(-)
 create mode 100644 reco/tasks/CbmTaskTriggerDigi.cxx
 create mode 100644 reco/tasks/CbmTaskTriggerDigi.h

diff --git a/core/base/CbmDigiManager.cxx b/core/base/CbmDigiManager.cxx
index 546b4cc3a0..6871158f3c 100644
--- a/core/base/CbmDigiManager.cxx
+++ b/core/base/CbmDigiManager.cxx
@@ -155,11 +155,11 @@ void CbmDigiManager::SetBranch()
   // --- Add branch object and connect it to the tree
   CbmDigiBranchBase* branch = new CbmDigiBranch<Digi>(branchName.c_str());
   if (branch->ConnectToTree()) {
-    LOG(info) << "DigiManager: Search branch " << branchName << " for class " << className << ": successful";
+    LOG(debug) << "DigiManager: Search branch " << branchName << " for class " << className << ": successful";
     fBranches[systemId] = branch;
   }
   else {
-    LOG(info) << "DigiManager: Search branch " << branchName << " for class " << className << ": failed";
+    LOG(debug) << "DigiManager: Search branch " << branchName << " for class " << className << ": failed";
     delete branch;
   }
 
diff --git a/macro/reco/reco_digi.C b/macro/reco/reco_digi.C
index 88b993476b..8b77b4967d 100644
--- a/macro/reco/reco_digi.C
+++ b/macro/reco/reco_digi.C
@@ -9,9 +9,9 @@
 #include "CbmDefs.h"
 #include "CbmSetup.h"
 #include "CbmTaskBuildEvents.h"
+#include "CbmTaskTriggerDigi.h"
 
 #include <FairFileSource.h>
-#include <FairLogger.h>
 #include <FairMonitor.h>
 #include <FairParAsciiFileIo.h>
 #include <FairParRootFileIo.h>
@@ -20,8 +20,11 @@
 #include <FairSystemInfo.h>
 
 #include <TStopwatch.h>
+
+#include <memory>
 #endif
 
+#include <FairLogger.h>
 
 /** @brief Macro for CBM reconstruction from digi level
  ** @author Volker Friese <v.friese@gsi.de>
@@ -62,7 +65,7 @@ void reco_digi(TString input = "", Int_t nTimeSlices = -1, Int_t firstTimeSlice
   // ------------------------------------------------------------------------
 
   // -----   Environment   --------------------------------------------------
-  TString myName = "run_reco_algo";                // this macro's name for screen output
+  TString myName = "reco_digi";                    // this macro's name for screen log
   TString srcDir = gSystem->Getenv("VMCWORKDIR");  // top source directory
   // ------------------------------------------------------------------------
 
@@ -76,9 +79,9 @@ void reco_digi(TString input = "", Int_t nTimeSlices = -1, Int_t firstTimeSlice
   TString monFile = output + ".moni_reco.root";
   if (paramFile.IsNull()) paramFile = input;
   TString parFile = paramFile + ".par.root";
-  std::cout << "Inputfile " << rawFile << std::endl;
-  std::cout << "Outfile " << outFile << std::endl;
-  std::cout << "Parfile " << parFile << std::endl;
+  LOG(info) << myName << ": Input file is     " << rawFile;
+  LOG(info) << myName << ": Output file is    " << outFile;
+  LOG(info) << myName << ": Parameter file is " << parFile;
 
   // -----   Load the geometry setup   -------------------------------------
   std::cout << std::endl;
@@ -110,11 +113,22 @@ void reco_digi(TString input = "", Int_t nTimeSlices = -1, Int_t firstTimeSlice
   // ------------------------------------------------------------------------
 
 
+  // -----   Digi trigger   -------------------------------------------------
+  auto trigger         = std::make_unique<CbmTaskTriggerDigi>();
+  double triggerWindow = 10.;  // Trigger window size in ns
+  int32_t minNumDigis  = 100;  // Trigger threshold in number of digis
+  double deadTime      = 50.;  // Minimum time between two triggers
+  trigger->SetAlgoParams(triggerWindow, minNumDigis, deadTime);
+  LOG(info) << myName << ": Added task " << trigger->GetName();
+  run->AddTask(trigger.release());
+  // ------------------------------------------------------------------------
+
+
   // -----   Event building   -----------------------------------------------
-  CbmTaskBuildEvents* evtBuild = new CbmTaskBuildEvents();
-  evtBuild->SetTimeWindow(-20., 30.);  // event build window for STS
-  run->AddTask(evtBuild);
-  LOG(info) << myName << ": Addes task " << evtBuild->GetName();
+  auto evtBuild = std::make_unique<CbmTaskBuildEvents>();
+  evtBuild->SetTimeWindow(-20., 30.);  // event building time window for STS
+  LOG(info) << myName << ": Added task " << evtBuild->GetName();
+  run->AddTask(evtBuild.release());
   // ------------------------------------------------------------------------
 
 
diff --git a/reco/tasks/CMakeLists.txt b/reco/tasks/CMakeLists.txt
index 310215ee90..e6e7d6dc31 100644
--- a/reco/tasks/CMakeLists.txt
+++ b/reco/tasks/CMakeLists.txt
@@ -10,6 +10,7 @@ Set(LIBRARY_NAME CbmRecoTasks)
 # -----  Compilation sources   ----------------------------
 set(SRCS
 CbmTaskBuildEvents.cxx
+CbmTaskTriggerDigi.cxx
 )
 # ---------------------------------------------------------
 
@@ -26,6 +27,7 @@ ${CBMROOT_SOURCE_DIR}/core/data/global
 ${CBMROOT_SOURCE_DIR}/core/data/sts
 
 ${CBMROOT_SOURCE_DIR}/algo/evbuild
+${CBMROOT_SOURCE_DIR}/algo/trigger
 
 )
 
diff --git a/reco/tasks/CbmRecoTasksLinkDef.h b/reco/tasks/CbmRecoTasksLinkDef.h
index 8decde767c..45e9ff31f2 100644
--- a/reco/tasks/CbmRecoTasksLinkDef.h
+++ b/reco/tasks/CbmRecoTasksLinkDef.h
@@ -12,5 +12,6 @@
 
 // --- Classes
 #pragma link C++ class CbmTaskBuildEvents + ;
+#pragma link C++ class CbmTaskTriggerDigi + ;
 
 #endif /* __CINT__ */
diff --git a/reco/tasks/CbmTaskBuildEvents.cxx b/reco/tasks/CbmTaskBuildEvents.cxx
index f5d3aa7a7f..b25049c375 100644
--- a/reco/tasks/CbmTaskBuildEvents.cxx
+++ b/reco/tasks/CbmTaskBuildEvents.cxx
@@ -61,17 +61,14 @@ void CbmTaskBuildEvents::Exec(Option_t*)
   timerStep.Stop();
   fTimeFillTs += timerStep.RealTime();
 
-  // --- Construct an artificial trigger list (just preliminary)
-  vector<double> triggerVec {1000., 2000., 3000.};
-
   // --- Call event builder algorithm
   timerStep.Start();
-  *fEvents = fAlgo(ts, triggerVec);
+  *fEvents = fAlgo(ts, *fTriggers);
   timerStep.Stop();
   fTimeBuildEvt += timerStep.RealTime();
 
   // --- Timeslice statistics
-  size_t numTriggers   = triggerVec.size();
+  size_t numTriggers   = fTriggers->size();
   size_t numEvents     = fEvents->size();
   size_t numDigisStsTs = ts.fData.fSts.fDigis.size();
   size_t numDigisStsEv = 0;
@@ -138,11 +135,20 @@ InitStatus CbmTaskBuildEvents::Init()
   LOG(info) << GetName() << ": Initialising...";
 
 
-  // --- Check input data (STS only for the time being)
+  // --- Check input data (digis, STS only for the time being)
   if (!fDigiMan->IsPresent(ECbmModuleId::kSts)) {
     LOG(fatal) << GetName() << ": No digi branch for STS";
     return kFATAL;
   }
+  LOG(info) << "--- Found branch STS digi";
+
+  // --- Get input data (triggers)
+  fTriggers = ioman->InitObjectAs<std::vector<double> const*>("Trigger");
+  if (!fTriggers) {
+    LOG(fatal) << GetName() << ": No Trigger branch!" << endl;
+    return kFATAL;
+  }
+  LOG(info) << "--- Found branch Trigger";
 
   // --- Register output array (CbmDigiEvent)
   if (ioman->GetObject("DigiEvent")) {
@@ -155,9 +161,11 @@ InitStatus CbmTaskBuildEvents::Init()
     LOG(fatal) << GetName() << ": Output branch could not be created!";
     return kFATAL;
   }
+  LOG(info) << "--- Registered branch DigiEvent";
 
   // --- Configure algorithm
   fAlgo.SetTriggerWindow(ECbmModuleId::kSts, fEvtTimeStsMin, fEvtTimeStsMax);
+  LOG(info) << "--- Use algo EventBuilder with event window [" << fEvtTimeStsMin << ", " << fEvtTimeStsMax << "] ns";
 
 
   LOG(info) << "==================================================";
diff --git a/reco/tasks/CbmTaskBuildEvents.h b/reco/tasks/CbmTaskBuildEvents.h
index 8e7e588de4..41eb11be66 100644
--- a/reco/tasks/CbmTaskBuildEvents.h
+++ b/reco/tasks/CbmTaskBuildEvents.h
@@ -74,11 +74,12 @@ private:  // methods
   virtual InitStatus Init();
 
 
-private:                                         // members
-  CbmDigiManager* fDigiMan = nullptr;            //! Input data
-  std::vector<ECbmModuleId> fSystems {};         //  List of detector systems
-  std::vector<CbmDigiEvent>* fEvents = nullptr;  //! Output data
-  cbm::algo::EventBuilder fAlgo {};              //! Algorithm
+private:                                           // members
+  CbmDigiManager* fDigiMan             = nullptr;  //! Input data (digis)
+  const std::vector<double>* fTriggers = nullptr;  //! Input data (triggers)
+  std::vector<ECbmModuleId> fSystems {};           //  List of detector systems
+  std::vector<CbmDigiEvent>* fEvents = nullptr;    //! Output data (events)
+  cbm::algo::EventBuilder fAlgo {};                //! Algorithm
   double fEvtTimeStsMin = 0.;
   double fEvtTimeStsMax = 0.;
   size_t fNumTs         = 0;  //  Number of processed time slices
diff --git a/reco/tasks/CbmTaskTriggerDigi.cxx b/reco/tasks/CbmTaskTriggerDigi.cxx
new file mode 100644
index 0000000000..0e34cf7e13
--- /dev/null
+++ b/reco/tasks/CbmTaskTriggerDigi.cxx
@@ -0,0 +1,155 @@
+/* Copyright (C) 2021 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Volker Friese [committer] */
+
+
+#include "CbmTaskTriggerDigi.h"
+
+#include "CbmDefs.h"
+#include "CbmDigiBranchBase.h"
+#include "CbmDigiManager.h"
+#include "CbmStsDigi.h"
+
+#include <FairLogger.h>
+
+#include "TimeClusterTrigger.h"
+#include <TStopwatch.h>
+
+#include <algorithm>
+#include <cassert>
+#include <iomanip>
+#include <iostream>
+#include <vector>
+
+
+using namespace std;
+
+
+// -----   Constructor   -----------------------------------------------------
+CbmTaskTriggerDigi::CbmTaskTriggerDigi() : FairTask("TriggerDigi") {}
+// ---------------------------------------------------------------------------
+
+
+// -----   Destructor   ------------------------------------------------------
+CbmTaskTriggerDigi::~CbmTaskTriggerDigi()
+{
+  if (fTriggers) delete fTriggers;
+}
+// ---------------------------------------------------------------------------
+
+
+// -----   Execution   -------------------------------------------------------
+void CbmTaskTriggerDigi::Exec(Option_t*)
+{
+
+  // --- Timer and counters
+  TStopwatch timerStep;
+  TStopwatch timerTot;
+  timerTot.Start();
+
+  // --- Get input digi vector
+  CbmDigiBranchBase* stsBranch      = fDigiMan->GetBranch(ECbmModuleId::kSts);
+  const vector<CbmStsDigi>* digiVec = boost::any_cast<const vector<CbmStsDigi>*>(stsBranch->GetBranchContainer());
+  assert(digiVec);
+
+  // --- Extract digi times into to a vector
+  timerStep.Start();
+  std::vector<double> digiTimes(digiVec->size());
+  std::transform(digiVec->begin(), digiVec->end(), digiTimes.begin(),
+                 [](const CbmStsDigi& digi) { return digi.GetTime(); });
+  timerStep.Stop();
+  fTimeExtract += timerStep.RealTime();
+
+  // --- Call the trigger algorithm
+  timerStep.Start();
+  *fTriggers = fAlgo(digiTimes, fTriggerWindow, fMinNumDigis, fDeadTime);
+  timerStep.Stop();
+  fTimeFind += timerStep.RealTime();
+
+  // --- Timeslice statistics
+  size_t numDigis    = digiVec->size();
+  size_t numTriggers = fTriggers->size();
+
+  // --- Timeslice log
+  timerTot.Stop();
+  fTimeTot += timerTot.RealTime();
+  stringstream logOut;
+  logOut << setw(20) << left << GetName() << " [";
+  logOut << fixed << setw(8) << setprecision(1) << right << timerTot.RealTime() * 1000. << " ms] ";
+  logOut << "TS " << fNumTs << ", digis " << numDigis << ", triggers " << numTriggers;
+  LOG(info) << logOut.str();
+
+  // --- Run statistics
+  fNumTs++;
+  fNumDigis += numDigis;
+  fNumTriggers += numTriggers;
+}
+// ----------------------------------------------------------------------------
+
+
+// -----   End-of-timeslice action   ------------------------------------------
+void CbmTaskTriggerDigi::Finish()
+{
+
+  std::cout << std::endl;
+  LOG(info) << "=====================================";
+  LOG(info) << GetName() << ": Run summary";
+  LOG(info) << "Timeslices         : " << fNumTs;
+  LOG(info) << "Digis              : " << fNumDigis;
+  LOG(info) << "Triggers           : " << fNumTriggers;
+  LOG(info) << "Time  / TS         : " << fixed << setprecision(2) << 1000. * fTimeTot / double(fNumTs) << " ms";
+  LOG(info) << "Time extract       : " << fixed << setprecision(2) << 1000. * fTimeExtract / double(fNumTs)
+            << " ms = " << 100. * fTimeExtract / fTimeTot << " %";
+  LOG(info) << "Time find trigger  : " << fixed << setprecision(2) << 1000. * fTimeFind / double(fNumTs)
+            << " ms = " << 100. * fTimeFind / fTimeTot << " %";
+  LOG(info) << "=====================================";
+}
+// ----------------------------------------------------------------------------
+
+
+// -----   Initialisation   ---------------------------------------------------
+InitStatus CbmTaskTriggerDigi::Init()
+{
+
+  // --- Get FairRootManager instance
+  FairRootManager* ioman = FairRootManager::Instance();
+  assert(ioman);
+
+  // --- DigiManager instance
+  fDigiMan = CbmDigiManager::Instance();
+  fDigiMan->Init();
+
+  std::cout << std::endl;
+  LOG(info) << "==================================================";
+  LOG(info) << GetName() << ": Initialising...";
+
+  // --- Check input data (STS only for the time being)
+  if (!fDigiMan->IsPresent(ECbmModuleId::kSts)) {
+    LOG(fatal) << GetName() << ": No digi branch for STS";
+    return kFATAL;
+  }
+  LOG(info) << "--- Found branch STS digi";
+
+  // --- Register output array (Triggers)
+  if (ioman->GetObject("Trigger")) {
+    LOG(fatal) << GetName() << ": Branch Trigger already exists!";
+    return kFATAL;
+  }
+  fTriggers = new vector<double>;
+  ioman->RegisterAny("Trigger", fTriggers, IsOutputBranchPersistent("Trigger"));
+  if (!fTriggers) {
+    LOG(fatal) << GetName() << ": Output branch could not be created!";
+    return kFATAL;
+  }
+  LOG(info) << "--- Registered branch Trigger";
+
+  LOG(info) << "--- Use algo TimeClusterTrigger with TriggerWin " << fTriggerWindow << " ns, threshold " << fMinNumDigis
+            << ", dead time " << fDeadTime << " ns";
+  LOG(info) << "==================================================";
+  std::cout << std::endl;
+
+  return kSUCCESS;
+}
+// ----------------------------------------------------------------------------
+
+ClassImp(CbmTaskTriggerDigi)
diff --git a/reco/tasks/CbmTaskTriggerDigi.h b/reco/tasks/CbmTaskTriggerDigi.h
new file mode 100644
index 0000000000..e5325fa82c
--- /dev/null
+++ b/reco/tasks/CbmTaskTriggerDigi.h
@@ -0,0 +1,95 @@
+/* Copyright (C) 2021 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Volker Friese [committer] */
+
+
+#ifndef CBMTASKTRIGGERDIGI_H
+#define CBMTASKTRIGGERDIGI_H 1
+
+
+#include "CbmDefs.h"
+
+#include <FairTask.h>
+
+#include "TimeClusterTrigger.h"
+
+#include <vector>
+
+class CbmDigiManager;
+
+
+/** @class CbmTaskTriggerDigi
+ ** @brief Task class for minimum-bias event trigger from time-distribution of digi data
+ ** @author Volker Friese <v.friese@gsi.de>
+ ** @since 20.11.2021
+ **
+ ** The tasks calls algo::TimeClusterTrigger with the digi time distribution of the trigger detector.
+ **
+ ** TOFO: The current implementation is for STS only as trigger detector.
+ **/
+class CbmTaskTriggerDigi : public FairTask {
+
+
+public:
+  /** @brief Constructor **/
+  CbmTaskTriggerDigi();
+
+
+  /** @brief Copy constructor (disabled) **/
+  CbmTaskTriggerDigi(const CbmTaskTriggerDigi&) = delete;
+
+
+  /** @brief Destructor **/
+  virtual ~CbmTaskTriggerDigi();
+
+
+  /** @brief Task execution **/
+  virtual void Exec(Option_t* opt);
+
+
+  /** @brief Finish timeslice **/
+  virtual void Finish();
+
+
+  /** @brief Assignment operator (disabled) **/
+  CbmTaskTriggerDigi& operator=(const CbmTaskTriggerDigi&) = delete;
+
+
+  /** @brief Configure the trigger algorithm
+   ** @param triggerWindow    Size of trigger window [ns]
+   ** @param minNumDigis      Trigger threshold on number of digis in trigger windows
+   ** @param deadTime         Minimum time between two triggers [ns]
+   **/
+  void SetAlgoParams(double triggerWindow, int32_t minNumDigis, double deadTime)
+  {
+    fTriggerWindow = triggerWindow;
+    fMinNumDigis   = minNumDigis;
+    fDeadTime      = deadTime;
+  }
+
+
+private:  // methods
+  /** @brief Task initialisation **/
+  virtual InitStatus Init();
+
+
+private:                                     // members
+  CbmDigiManager* fDigiMan = nullptr;        //! Input data
+  std::vector<ECbmModuleId> fSystems {};     //  List of detector systems
+  std::vector<double>* fTriggers = nullptr;  //! Output data
+  cbm::algo::TimeClusterTrigger fAlgo {};    //! Algorithm
+  double fTriggerWindow = 0.;
+  int32_t fMinNumDigis  = 0;
+  double fDeadTime      = 0.;
+  size_t fNumTs         = 0;  //  Number of processed time slices
+  size_t fNumDigis      = 0;  //  Number of digis from trigger detector
+  size_t fNumTriggers   = 0;  //  Number of found triggers
+  double fTimeExtract   = 0.;
+  double fTimeFind      = 0.;
+  double fTimeTot       = 0.;
+
+
+  ClassDef(CbmTaskTriggerDigi, 1);
+};
+
+#endif /* CBMTASKTRIGGERDIGI_H */
-- 
GitLab