From 1664b48dd3a80bd3db91c54e408de45d958e3736 Mon Sep 17 00:00:00 2001
From: Volker Friese <v.friese@gsi.de>
Date: Sat, 12 Mar 2022 08:19:09 +0100
Subject: [PATCH] Run macro for reconstruction from timeslice: unpacker,
 trigger, event builder. Refs #2275.

---
 algo/evbuild/EventBuilder.cxx     |  16 +---
 macro/reco/reco_ts.C              | 149 ++++++++++++++++++++++++++++++
 reco/tasks/CbmTaskBuildEvents.cxx |  32 +++----
 reco/tasks/CbmTaskTriggerDigi.cxx |   4 +-
 reco/tasks/CbmTaskUnpack.cxx      |  19 ++--
 5 files changed, 174 insertions(+), 46 deletions(-)
 create mode 100644 macro/reco/reco_ts.C

diff --git a/algo/evbuild/EventBuilder.cxx b/algo/evbuild/EventBuilder.cxx
index 04455c4d9b..229c718c11 100644
--- a/algo/evbuild/EventBuilder.cxx
+++ b/algo/evbuild/EventBuilder.cxx
@@ -16,7 +16,6 @@ namespace cbm
     // --- Execution
     std::vector<CbmDigiEvent> EventBuilder::operator()(const CbmDigiTimeslice& ts, const vector<double> triggers) const
     {
-      assert(is_sorted(triggers.begin(), triggers.end()));
       vector<CbmDigiEvent> eventVec(triggers.size());
       std::transform(triggers.begin(), triggers.end(), eventVec.begin(),
                      [&ts, this](const double& trigger) { return BuildEvent(ts, trigger); });
@@ -39,43 +38,30 @@ namespace cbm
         // --- Build the event using trigger window
         switch (system) {
           case ECbmModuleId::kSts: {
-            assert(
-              is_sorted(ts.fData.fSts.fDigis.begin(), ts.fData.fSts.fDigis.end(), IsBefore<CbmStsDigi, CbmStsDigi>));
             event.fData.fSts.fDigis = CopyRange(ts.fData.fSts.fDigis, tMin, tMax);
             break;
           }
           case ECbmModuleId::kRich: {
-            assert(is_sorted(ts.fData.fRich.fDigis.begin(), ts.fData.fRich.fDigis.end(),
-                             IsBefore<CbmRichDigi, CbmRichDigi>));
             event.fData.fRich.fDigis = CopyRange(ts.fData.fRich.fDigis, tMin, tMax);
             break;
           }
           case ECbmModuleId::kMuch: {
-            assert(is_sorted(ts.fData.fMuch.fDigis.begin(), ts.fData.fMuch.fDigis.end(),
-                             IsBefore<CbmMuchDigi, CbmMuchDigi>));
             event.fData.fMuch.fDigis = CopyRange(ts.fData.fMuch.fDigis, tMin, tMax);
             break;
           }
           case ECbmModuleId::kTrd: {
-            assert(
-              is_sorted(ts.fData.fTrd.fDigis.begin(), ts.fData.fTrd.fDigis.end(), IsBefore<CbmTrdDigi, CbmTrdDigi>));
             event.fData.fTrd.fDigis = CopyRange(ts.fData.fTrd.fDigis, tMin, tMax);
             break;
           }
           case ECbmModuleId::kTof: {
-            assert(
-              is_sorted(ts.fData.fTof.fDigis.begin(), ts.fData.fTof.fDigis.end(), IsBefore<CbmTofDigi, CbmTofDigi>));
             event.fData.fTof.fDigis = CopyRange(ts.fData.fTof.fDigis, tMin, tMax);
             break;
           }
           case ECbmModuleId::kPsd: {
-            assert(
-              is_sorted(ts.fData.fPsd.fDigis.begin(), ts.fData.fPsd.fDigis.end(), IsBefore<CbmPsdDigi, CbmPsdDigi>));
             event.fData.fPsd.fDigis = CopyRange(ts.fData.fPsd.fDigis, tMin, tMax);
             break;
           }
-          case ECbmModuleId::kT0: {  //T0 has Tof digis
-            assert(is_sorted(ts.fData.fT0.fDigis.begin(), ts.fData.fT0.fDigis.end(), IsBefore<CbmTofDigi, CbmTofDigi>));
+          case ECbmModuleId::kT0: {
             event.fData.fT0.fDigis = CopyRange(ts.fData.fT0.fDigis, tMin, tMax);
             break;
           }
diff --git a/macro/reco/reco_ts.C b/macro/reco/reco_ts.C
new file mode 100644
index 0000000000..df96f89e59
--- /dev/null
+++ b/macro/reco/reco_ts.C
@@ -0,0 +1,149 @@
+/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Volker Friese [committer] */
+
+
+// --- Includes needed for IDE code analyser
+#if !defined(__CLING__)
+#include "CbmSourceTs.h"
+#include "CbmTaskUnpack.h"
+#include "CbmTsEventHeader.h"
+
+#include <FairRootFileSink.h>
+#include <FairRunOnline.h>
+
+#include <TStopwatch.h>
+#include <TSystem.h>
+
+#include <iostream>
+#include <memory>
+#endif
+
+#include <FairLogger.h>
+
+
+using std::string;
+
+/** @brief Macro for CBM reconstruction from FLES timeslices
+ ** @author Volker Friese <v.friese@gsi.de>
+ ** @since  12 March 2022
+ ** @param tsaFile    Name of input file (w/o extension .tsa)
+ ** @param outFile    Name of output file (w/o extension .digi.root)
+ **
+ ** Reconstruction from timeslice level. Currently included stages:
+ ** - Unpacking (STS only)
+ ** - Event trigger based on STS digis (CbmTaskDigiTrigger)
+ ** - Event building (CbmTaskBuildEvents) (STS only)
+ **
+ ** TODO: To be expanded with the progress of the algo project.
+ **
+ ** If the tsaFile name is not specified, a default file from the reporistory will be used.
+ ** If the outFile name is not specified, the input file name will be used, replacing
+ ** the extension .tsa by .digi.root
+ **/
+
+void reco_ts(TString tsaFile = "", TString outFile = "")
+{
+
+  // ========================================================================
+  //          Adjust this part according to your requirements
+
+  // --- Logger settings ----------------------------------------------------
+  TString logLevel     = "INFO";
+  TString logVerbosity = "LOW";
+  // ------------------------------------------------------------------------
+
+  // -----   Environment   --------------------------------------------------
+  TString myName = "reco_ts";                      // this macro's name for screen output
+  TString srcDir = gSystem->Getenv("VMCWORKDIR");  // top source directory
+  // ------------------------------------------------------------------------
+
+  // In general, the following parts need not be touched
+  // ========================================================================
+
+  tsaFile = "/Users/vfriese/Cbmroot/data/1588_node8_1_0000";
+  outFile = "test";
+
+  // ----- Default file names   ---------------------------------------------
+  if (tsaFile.IsNull()) tsaFile = srcDir + "/input/mcbm_run399_first20Ts";
+  TString inFile = tsaFile + ".tsa";
+  if (outFile.IsNull()) outFile = tsaFile;
+  outFile += ".digi.root";
+  std::cout << std::endl << std::endl;
+  LOG(info) << myName << ": Using input file " << inFile;
+  LOG(info) << myName << ": Output file is   " << outFile;
+  // ------------------------------------------------------------------------
+
+
+  // -----   Timer   --------------------------------------------------------
+  TStopwatch timer;
+  timer.Start();
+  // ------------------------------------------------------------------------
+
+
+  // ----   Run configuration   ---------------------------------------------
+  auto source = new CbmSourceTs(inFile.Data());
+  auto run    = new FairRunOnline(source);
+  auto sink   = new FairRootFileSink(outFile.Data());
+  run->SetSink(sink);
+  auto eventheader = new CbmTsEventHeader();
+  run->SetEventHeader(new CbmTsEventHeader());
+  // ------------------------------------------------------------------------
+
+
+  // -----   Unpacking   ----------------------------------------------------
+  auto unpack = std::make_unique<CbmTaskUnpack>();
+  LOG(info) << myName << ": Added task " << unpack->GetName();
+  run->AddTask(unpack.release());
+  // ------------------------------------------------------------------------
+
+
+  // -----   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);
+  trigger->AddSystem(ECbmModuleId::kSts);
+  LOG(info) << myName << ": Added task " << trigger->GetName();
+  run->AddTask(trigger.release());
+  // ------------------------------------------------------------------------
+
+
+  // -----   Event building   -----------------------------------------------
+  auto evtBuild = std::make_unique<CbmTaskBuildEvents>();
+  evtBuild->SetEventWindow(ECbmModuleId::kSts, -20., 30.);  // event building time window for STS
+  LOG(info) << myName << ": Added task " << evtBuild->GetName();
+  run->AddTask(evtBuild.release());
+  // ------------------------------------------------------------------------
+
+
+  // -----   Logger settings   ----------------------------------------------
+  FairLogger::GetLogger()->SetLogScreenLevel(logLevel.Data());
+  FairLogger::GetLogger()->SetLogVerbosityLevel(logVerbosity.Data());
+  // ------------------------------------------------------------------------
+
+
+  // -----   Run initialisation   -------------------------------------------
+  std::cout << std::endl;
+  LOG(info) << myName << ": Initialise run" << std::endl;
+  run->Init();
+  // ------------------------------------------------------------------------
+
+
+  // -----   Start run   ----------------------------------------------------
+  std::cout << std::endl << std::endl;
+  LOG(info) << myName << ": Starting run" << std::endl;
+  run->Run(-1, 0);
+  // ------------------------------------------------------------------------
+
+
+  // -----   Finish   -------------------------------------------------------
+  timer.Stop();
+  std::cout << std::endl << std::endl;
+  std::cout << "Macro finished successfully." << std::endl;
+  std::cout << "Output file: " << outFile << std::endl;
+  std::cout << "CPU time = " << timer.CpuTime() << " s, real time = " << timer.RealTime() << " s." << std::endl;
+  // ------------------------------------------------------------------------
+
+}  // End of main macro function
diff --git a/reco/tasks/CbmTaskBuildEvents.cxx b/reco/tasks/CbmTaskBuildEvents.cxx
index a3983dd65f..dc18dac91d 100644
--- a/reco/tasks/CbmTaskBuildEvents.cxx
+++ b/reco/tasks/CbmTaskBuildEvents.cxx
@@ -113,8 +113,6 @@ void CbmTaskBuildEvents::Exec(Option_t*)
 
   // --- If the input is already CbmDigiTimeslice (from unpacking), use that directly
   if (fTimeslice) {
-    LOG(info) << "Exec timeslice " << fNumTs << " with " << fTriggers->size() << " triggers"
-              << " from " << *(fTriggers->begin()) << " to " << fTriggers->back();
     timerStep.Start();
     *fEvents = fAlgo(*fTimeslice, *fTriggers);
     timerStep.Stop();
@@ -150,7 +148,7 @@ void CbmTaskBuildEvents::Exec(Option_t*)
   timerTot.Stop();
   fTimeTot += timerTot.RealTime();
   stringstream logOut;
-  logOut << setw(20) << left << GetName() << " [";
+  logOut << setw(15) << left << GetName() << " [";
   logOut << fixed << setw(8) << setprecision(1) << right << timerTot.RealTime() * 1000. << " ms] ";
   logOut << "TS " << fNumTs << ", triggers " << numTriggers << ", events " << numEvents;
 
@@ -182,21 +180,21 @@ void CbmTaskBuildEvents::Finish()
   std::cout << std::endl;
   LOG(info) << "=====================================";
   LOG(info) << GetName() << ": Run summary";
-  LOG(info) << "Timeslices         : " << fNumTs;
-  LOG(info) << "Triggers           : " << fNumTriggers;
-  LOG(info) << "Events             : " << fNumEvents;
-
+  LOG(info) << "Timeslices           : " << fNumTs;
+  LOG(info) << "Triggers             : " << fNumTriggers;
+  LOG(info) << "Events               : " << fNumEvents;
   for (const auto& entry : fEventWindows) {
     auto system = entry.first;
-    LOG(info) << CbmModuleList::GetModuleNameCaps(system) << " digis in timeslice : " << fNumDigisTs[system];
-    LOG(info) << CbmModuleList::GetModuleNameCaps(system) << " digis in events    : " << fNumDigisEv[system] << " = "
-              << fixed << setprecision(2) << 100. * double(fNumDigisEv[system]) / double(fNumDigisTs[system]) << " %";
+    LOG(info) << setw(4) << left << CbmModuleList::GetModuleNameCaps(system)
+              << " digis in TS     : " << fNumDigisTs[system];
+    LOG(info) << setw(4) << left << CbmModuleList::GetModuleNameCaps(system)
+              << " digis in events : " << fNumDigisEv[system] << " = " << fixed << setprecision(2)
+              << 100. * double(fNumDigisEv[system]) / double(fNumDigisTs[system]) << " %";
   }
-
-  LOG(info) << "Time  / TS         : " << fixed << setprecision(2) << 1000. * fTimeTot / double(fNumTs) << " ms";
-  LOG(info) << "Time fill TS       : " << fixed << setprecision(2) << 1000. * fTimeFillTs / double(fNumTs)
+  LOG(info) << "Time  / TS           : " << fixed << setprecision(2) << 1000. * fTimeTot / double(fNumTs) << " ms";
+  LOG(info) << "Time fill TS         : " << fixed << setprecision(2) << 1000. * fTimeFillTs / double(fNumTs)
             << " ms = " << 100. * fTimeFillTs / fTimeTot << " %";
-  LOG(info) << "Time build events  : " << fixed << setprecision(2) << 1000. * fTimeBuildEvt / double(fNumTs)
+  LOG(info) << "Time build events    : " << fixed << setprecision(2) << 1000. * fTimeBuildEvt / double(fNumTs)
             << " ms = " << 100. * fTimeBuildEvt / fTimeTot << " %";
   LOG(info) << "=====================================";
 }
@@ -229,10 +227,6 @@ InitStatus CbmTaskBuildEvents::Init()
   FairRootManager* ioman = FairRootManager::Instance();
   assert(ioman);
 
-  // --- DigiManager instance
-  fDigiMan = CbmDigiManager::Instance();
-  fDigiMan->Init();
-
   std::cout << std::endl;
   LOG(info) << "==================================================";
   LOG(info) << GetName() << ": Initialising...";
@@ -240,7 +234,7 @@ InitStatus CbmTaskBuildEvents::Init()
   // --- Check input data
   // --- DigiTimeslice: Unpacked data from FLES
   fTimeslice = ioman->InitObjectAs<const CbmDigiTimeslice*>("DigiTimeslice.");
-  if (fTimeslice) { LOG(info) << GetName() << ": Found branch DigiTimeslice."; }
+  if (fTimeslice) { LOG(info) << "--- Found branch DigiTimeslice."; }
   // --- DigiManager: Simulated digi data
   else {
     fDigiMan = CbmDigiManager::Instance();
diff --git a/reco/tasks/CbmTaskTriggerDigi.cxx b/reco/tasks/CbmTaskTriggerDigi.cxx
index de60973e9f..82f9cc1517 100644
--- a/reco/tasks/CbmTaskTriggerDigi.cxx
+++ b/reco/tasks/CbmTaskTriggerDigi.cxx
@@ -115,7 +115,7 @@ void CbmTaskTriggerDigi::Exec(Option_t*)
   timerTot.Stop();
   fTimeTot += timerTot.RealTime();
   stringstream logOut;
-  logOut << setw(20) << left << GetName() << " [";
+  logOut << setw(15) << left << GetName() << " [";
   logOut << fixed << setw(8) << setprecision(1) << right << timerTot.RealTime() * 1000. << " ms] ";
   logOut << "TS " << fNumTs << ", digis " << numDigis << ", triggers " << numTriggers;
   LOG(info) << logOut.str();
@@ -229,7 +229,7 @@ InitStatus CbmTaskTriggerDigi::Init()
   // --- Check input data
   // --- DigiTimeslice: Unpacked data from FLES
   fTimeslice = ioman->InitObjectAs<const CbmDigiTimeslice*>("DigiTimeslice.");
-  if (fTimeslice) { LOG(info) << GetName() << ": Found branch DigiTimeslice."; }
+  if (fTimeslice) { LOG(info) << "--- Found branch DigiTimeslice."; }
   // --- DigiManager: Simulated digi data
   else {
     fDigiMan = CbmDigiManager::Instance();
diff --git a/reco/tasks/CbmTaskUnpack.cxx b/reco/tasks/CbmTaskUnpack.cxx
index 727dca5bd5..3250f4d286 100644
--- a/reco/tasks/CbmTaskUnpack.cxx
+++ b/reco/tasks/CbmTaskUnpack.cxx
@@ -107,8 +107,6 @@ void CbmTaskUnpack::Exec(Option_t*)
                     << result.second.fNumErrElinkOutOfRange << " | " << result.second.fNumErrInvalidFirstMessage
                     << " | " << result.second.fNumErrInvalidMsSize << " | " << result.second.fNumErrTimestampOverflow
                     << " | ";
-        //std::move(result.first.begin(), result.first.end(), fTimeslice->fData.fSts.fDigis.end());
-        // TODO: The above usage of std::move does not work (seg. fault). Would need advice.
         const auto it = fTimeslice->fData.fSts.fDigis.end();
         fTimeslice->fData.fSts.fDigis.insert(it, result.first.begin(), result.first.end());
         numDigisInComp += result.first.size();
@@ -127,6 +125,10 @@ void CbmTaskUnpack::Exec(Option_t*)
 
   }  //# component
 
+  // --- Sorting of output digis. Is required by both digi trigger and event builder.
+  std::sort(fTimeslice->fData.fSts.fDigis.begin(), fTimeslice->fData.fSts.fDigis.end(),
+            [](CbmStsDigi digi1, CbmStsDigi digi2) { return digi1.GetTime() < digi2.GetTime(); });
+
 
   // --- Timeslice log
   timer.Stop();
@@ -173,7 +175,7 @@ InitStatus CbmTaskUnpack::Init()
 
   std::cout << std::endl;
   LOG(info) << "==================================================";
-  LOG(info) << GetName() << ": Init";
+  LOG(info) << GetName() << ": Initialising...";
 
   // --- Get source instance
   fSource = dynamic_cast<CbmSourceTs*>(FairRunOnline::Instance()->GetSource());
@@ -181,7 +183,7 @@ InitStatus CbmTaskUnpack::Init()
     LOG(error) << GetName() << ": No valid source class registered!";
     return kFATAL;
   }
-  LOG(info) << GetName() << ": Found CbmSourceTs instance";
+  LOG(info) << "--- Found CbmSourceTs instance";
 
   // --- Get FairRootManager instance
   FairRootManager* ioman = FairRootManager::Instance();
@@ -194,10 +196,7 @@ InitStatus CbmTaskUnpack::Init()
   }
   fTimeslice = new CbmDigiTimeslice();
   ioman->RegisterAny("DigiTimeslice.", fTimeslice, kTRUE);
-  LOG(info) << GetName() << ": Registered branch DigiTimeslice.";
-
-  // --- Initialise STS readout configuration
-  //InitStsConfig();
+  LOG(info) << "--- Registered branch DigiTimeslice.";
 
   // --- Common parameters for all components
   uint32_t numChansPerAsic   = 128;  // R/O channels per ASIC
@@ -222,10 +221,10 @@ InitStatus CbmTaskUnpack::Init()
       par->fElinkParams.push_back(elinkPar);
     }
     fAlgoSts[equip].SetParams(std::move(par));
-    LOG(info) << GetName() << ": configured equipment " << equip << " with " << numElinks << " elinks";
+    LOG(info) << "--- Configured equipment " << equip << " with " << numElinks << " elinks";
   }  //# equipments
 
-  LOG(info) << GetName() << ": configured " << fAlgoSts.size() << " unpacker algorithms for STS.";
+  LOG(info) << "--- Configured " << fAlgoSts.size() << " unpacker algorithms for STS.";
   LOG(debug) << "Readout map:" << fStsConfig.PrintReadoutMap();
   LOG(info) << "==================================================";
   std::cout << std::endl;
-- 
GitLab