From 3068b333b5cbf72674dfba94eef9ed76fa2b02d7 Mon Sep 17 00:00:00 2001
From: Volker Friese <v.friese@gsi.de>
Date: Tue, 11 Jul 2023 10:05:47 +0000
Subject: [PATCH] Simple QA of DigiEvents implemented in algo, with offline
 integration into cbmreco_fairrun.

---
 algo/CMakeLists.txt               |   3 +
 algo/qa/DigiEventQa.cxx           |  56 +++++++++++++
 algo/qa/DigiEventQa.h             | 112 +++++++++++++++++++++++++
 algo/qa/Histo1D.cxx               | 122 ++++++++++++++++++++++++++++
 algo/qa/Histo1D.h                 | 131 ++++++++++++++++++++++++++++++
 reco/tasks/CMakeLists.txt         |   1 +
 reco/tasks/CbmTaskDigiEventQa.cxx | 118 +++++++++++++++------------
 reco/tasks/CbmTaskDigiEventQa.h   |  25 ++++--
 8 files changed, 510 insertions(+), 58 deletions(-)
 create mode 100644 algo/qa/DigiEventQa.cxx
 create mode 100644 algo/qa/DigiEventQa.h
 create mode 100644 algo/qa/Histo1D.cxx
 create mode 100644 algo/qa/Histo1D.h

diff --git a/algo/CMakeLists.txt b/algo/CMakeLists.txt
index f09ec98a1d..ab468ea286 100644
--- a/algo/CMakeLists.txt
+++ b/algo/CMakeLists.txt
@@ -38,6 +38,8 @@ set(SRCS
   detectors/rich/RichReadoutConfig.cxx
   detectors/rich/UnpackRich.cxx
   global/Reco.cxx
+  qa/DigiEventQa.cxx
+  qa/Histo1D.cxx
 )
 
 set(BUILD_INFO_CXX ${CMAKE_CURRENT_BINARY_DIR}/base/BuildInfo.cxx)
@@ -68,6 +70,7 @@ target_include_directories(Algo
          ${CMAKE_CURRENT_SOURCE_DIR}/evselector
          ${CMAKE_CURRENT_SOURCE_DIR}/unpack
          ${CMAKE_CURRENT_SOURCE_DIR}/detectors
+         ${CMAKE_CURRENT_SOURCE_DIR}/qa
  )
 
 target_link_libraries(Algo
diff --git a/algo/qa/DigiEventQa.cxx b/algo/qa/DigiEventQa.cxx
new file mode 100644
index 0000000000..284b1d8cce
--- /dev/null
+++ b/algo/qa/DigiEventQa.cxx
@@ -0,0 +1,56 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Volker Friese [committer] */
+
+#include "DigiEventQa.h"
+
+using std::string;
+using std::vector;
+
+namespace cbm::algo::qa
+{
+
+  // ---   Execution   --------------------------------------------------------
+  DigiEventQaData DigiEventQa::operator()(const vector<CbmDigiEvent>& events, const DigiEventQaConfig& config) const
+  {
+
+    // --- Instantiate return object
+    DigiEventQaData result;
+    for (const auto& entry : config.fData) {
+      ECbmModuleId subsystem = entry.first;
+      string hist_name       = "digi_time_" + ToString(subsystem);
+      auto& detConfig        = config.fData.at(subsystem);
+      result.fDigiTimeHistos.try_emplace(subsystem, detConfig.fNumBins, detConfig.fMinValue, detConfig.fMaxValue,
+                                         hist_name);
+    }
+
+    // --- Event loop. Fill histograms.
+    for (const auto& event : events)
+      for (auto& subsystem : config.fData)
+        QaDigiTimeInEvent(event, subsystem.first, result.fDigiTimeHistos.at(subsystem.first));
+    result.fNumEvents = events.size();
+
+    return result;
+  }
+  // --------------------------------------------------------------------------
+
+
+  // ---  QA: digi time within event   ----------------------------------------
+  void DigiEventQa::QaDigiTimeInEvent(const CbmDigiEvent& event, ECbmModuleId system, Histo1D& histo) const
+  {
+    switch (system) {
+      case ECbmModuleId::kT0: FillDeltaT<CbmBmonDigi>(event.fData.fT0.fDigis, event.fTime, histo); break;
+      case ECbmModuleId::kSts: FillDeltaT<CbmStsDigi>(event.fData.fSts.fDigis, event.fTime, histo); break;
+      case ECbmModuleId::kMuch: FillDeltaT<CbmMuchDigi>(event.fData.fMuch.fDigis, event.fTime, histo); break;
+      case ECbmModuleId::kRich: FillDeltaT<CbmRichDigi>(event.fData.fRich.fDigis, event.fTime, histo); break;
+      case ECbmModuleId::kTrd: FillDeltaT<CbmTrdDigi>(event.fData.fTrd.fDigis, event.fTime, histo); break;
+      case ECbmModuleId::kTrd2d: FillDeltaT<CbmTrdDigi>(event.fData.fTrd2d.fDigis, event.fTime, histo); break;
+      case ECbmModuleId::kTof: FillDeltaT<CbmTofDigi>(event.fData.fTof.fDigis, event.fTime, histo); break;
+      case ECbmModuleId::kPsd: FillDeltaT<CbmPsdDigi>(event.fData.fPsd.fDigis, event.fTime, histo); break;
+      default: throw std::runtime_error("DigiEventQa: Invalid system Id " + ToString(system));
+    }
+  }
+  // --------------------------------------------------------------------------
+
+
+} /* namespace cbm::algo::qa */
diff --git a/algo/qa/DigiEventQa.h b/algo/qa/DigiEventQa.h
new file mode 100644
index 0000000000..1ba0e75e5e
--- /dev/null
+++ b/algo/qa/DigiEventQa.h
@@ -0,0 +1,112 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Volker Friese [committer] */
+
+#ifndef ALGO_QA_DIGIEVENTQA_H
+#define ALGO_QA_DIGIEVENTQA_H 1
+
+#include "CbmDefs.h"
+#include "CbmDigiEvent.h"
+
+#include <map>
+#include <vector>
+
+#include "Histo1D.h"
+
+
+namespace cbm::algo::qa
+{
+
+  /** @struct DigiEventQaData
+   ** @brief QA results for CbmDigiEvent objects
+   ** @author Volker Friese <v.friese@gsi.de>
+   ** @since 16 June 2023
+   **/
+  struct DigiEventQaData {
+    std::map<ECbmModuleId, Histo1D> fDigiTimeHistos = {};
+    size_t fNumEvents                               = 0;
+  };
+
+
+  /** @struct DigiEventQaDetConfig
+   ** @brief Configuration data for the QA of CbmDigiEvents for a given detector
+   ** @author Volker Friese <v.friese@gsi.de>
+   ** @since 16 June 2023
+   **/
+  struct DigiEventQaDetConfig {
+    uint32_t fNumBins;
+    double fMinValue;
+    double fMaxValue;
+    std::string ToString() const
+    {
+      std::stringstream ss;
+      ss << "nbins " << fNumBins << " min " << fMinValue << " max " << fMaxValue << "\n";
+      return ss.str();
+    }
+  };
+
+
+  /** @struct DigiEventQaConfig
+   ** @brief Configuration data for the QA of CbmDigiEvents
+   ** @author Volker Friese <v.friese@gsi.de>
+   ** @since 16 June 2023
+   **/
+  struct DigiEventQaConfig {
+    std::map<ECbmModuleId, DigiEventQaDetConfig> fData;
+    std::string ToString() const
+    {
+      std::stringstream ss;
+      for (const auto& entry : fData)
+        ss << "\n   Subsystem " << ::ToString(entry.first) << "  " << entry.second.ToString();
+      return ss.str();
+    }
+  };
+
+
+  /** @class DigiEventQa
+   ** @brief QA for CbmDigiEvent objects
+   ** @author Volker Friese <v.friese@gsi.de>
+   ** @since 16 June 2023
+   **/
+  class DigiEventQa {
+  public:
+    /** @brief Constructor **/
+    DigiEventQa() = default;
+
+    /** @brief Destructor **/
+    virtual ~DigiEventQa() = default;
+
+    /** @brief Execution
+     ** @param events Vector of DigiEvents to analyse
+     ** @return QA data object
+     **/
+    DigiEventQaData operator()(const std::vector<CbmDigiEvent>& events, const DigiEventQaConfig& config) const;
+
+
+  private:
+    /** @brief Fill histogram with digi time within event
+     ** @param digis  Vector with digi objects
+     ** @param eventTime  Time of event
+     ** @param histo  Histogram to be filled
+     **
+     ** The templated class is required to implement the method double GetTime().
+     **/
+    template<class Digi>
+    void FillDeltaT(const std::vector<Digi>& digis, double eventTime, Histo1D& histo) const
+    {
+      for (const Digi& digi : digis)
+        histo.Add(digi.GetTime() - eventTime);
+    }
+
+    /** @brief Fill histogram with digi time within event
+     ** @param digis  Vector with digi objects
+     ** @param eventTime  Time of event
+     ** @param histo  Histogram to be filled
+     **/
+    void QaDigiTimeInEvent(const CbmDigiEvent& event, ECbmModuleId system, Histo1D& histo) const;
+  };
+
+
+} /* namespace cbm::algo::qa */
+
+#endif /* ALGO_QA_DIGIEVENTQA_H */
diff --git a/algo/qa/Histo1D.cxx b/algo/qa/Histo1D.cxx
new file mode 100644
index 0000000000..9a02940b61
--- /dev/null
+++ b/algo/qa/Histo1D.cxx
@@ -0,0 +1,122 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Volker Friese [committer] */
+
+#include "Histo1D.h"
+
+#include <cassert>
+#include <cstdint>
+#include <sstream>
+#include <stdexcept>
+#include <string>
+#include <vector>
+
+#include <cmath>
+
+
+namespace cbm::algo
+{
+
+  // -----   Standard constructor
+  Histo1D::Histo1D(uint32_t numBins, double minValue, double maxValue, const std::string& name)
+    : fName(name)
+    , fNumBins(numBins)
+    , fMinValue(minValue)
+    , fMaxValue(maxValue)
+  {
+    if (numBins == 0 || maxValue <= minValue)
+      throw std::runtime_error("Histo1D: Invalid specifications " + std::to_string(numBins) + " "
+                               + std::to_string(minValue) + " " + std::to_string(maxValue));
+    fContent.resize(numBins);
+  }
+
+  // -----   Add entry
+  void Histo1D::Add(double value, double weight)
+  {
+    if (value < fMinValue) fUnderflow += weight;
+    else if (!(value < fMaxValue))
+      fOverflow += weight;
+    else {
+      uint32_t bin = uint32_t(double(fNumBins) * (value - fMinValue) / (fMaxValue - fMinValue));
+      assert(bin < fNumBins);
+      fContent[bin] += weight;
+      fNumEntries++;
+    }
+  }
+
+  // -----   Clear content
+  void Histo1D::Clear()
+  {
+    fContent.assign(fNumBins, 0.);
+    fUnderflow  = 0;
+    fOverflow   = 0;
+    fNumEntries = 0;
+  }
+
+  // -----   Content access
+  double Histo1D::Content(uint32_t bin) const
+  {
+    if (bin >= fNumBins) return 0.;
+    else
+      return fContent[bin];
+  }
+
+
+  // -----   First moment of distribution
+  double Histo1D::Mean() const
+  {
+    double sum1    = 0.;
+    double sum2    = 0.;
+    double binsize = (fMaxValue - fMinValue) / double(fNumBins);
+    for (uint32_t bin = 0; bin < fNumBins; bin++) {
+      double x = fMinValue + (bin + 0.5) * binsize;
+      sum1 += x * fContent[bin];
+      sum2 += fContent[bin];
+    }
+    return sum1 / sum2;
+  }
+
+
+  // -----   Operator +=
+  // TODO: Comparison of floating point numbers; probably fishy.
+  Histo1D& Histo1D::operator+=(const Histo1D& other)
+  {
+    if (other.fNumBins == fNumBins && other.fMinValue == fMinValue && other.fMaxValue == fMaxValue) {
+      fUnderflow += other.fUnderflow;
+      fOverflow += other.fOverflow;
+      for (uint32_t bin = 0; bin < fNumBins; bin++)
+        fContent[bin] += other.fContent[bin];
+      fNumEntries += other.fNumEntries;
+    }
+    else
+      throw std::runtime_error("Histo1D: Trying to add incompatible histograms");
+    return *this;
+  }
+
+
+  // -----   Second moment of distribution
+  double Histo1D::Stddev() const
+  {
+    double sum1    = 0.;
+    double sum2    = 0.;
+    double binsize = (fMaxValue - fMinValue) / double(fNumBins);
+    for (uint32_t bin = 0; bin < fNumBins; bin++) {
+      double x = fMinValue + (bin + 0.5) * binsize;
+      sum1 += x * x * fContent[bin];
+      sum2 += fContent[bin];
+    }
+    double mean = Mean();
+    return sqrt(sum1 / sum2 - mean * mean);
+  }
+
+
+  // -----   Properties to string
+  std::string Histo1D::ToString() const
+  {
+    std::stringstream ss;
+    ss << fName << ": bins " << fNumBins << " range " << fMinValue << " to " << fMaxValue << " entries " << fNumEntries
+       << " mean " << Mean() << " stddev " << Stddev() << " out of range " << fUnderflow << " , " << fOverflow;
+    return ss.str();
+  }
+
+} /* namespace cbm::algo */
diff --git a/algo/qa/Histo1D.h b/algo/qa/Histo1D.h
new file mode 100644
index 0000000000..7f6b8f35ac
--- /dev/null
+++ b/algo/qa/Histo1D.h
@@ -0,0 +1,131 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Volker Friese [committer] */
+
+#ifndef ALGO_QA_HISTO1D_H
+#define ALGO_QA_HISTO1D_H 1
+
+#include <cstdint>
+#include <string>
+#include <vector>
+
+
+namespace cbm::algo
+{
+
+  /** @class Histo1D
+   ** @author Volker Friese <v.friese@gsi.de>
+   ** @since 15 June 2023
+   **
+   ** Lightweight one-dimensional histogram class
+   **/
+  class Histo1D {
+  public:
+    /** @brief No default constructor **/
+    Histo1D() = delete;
+
+
+    /** @brief Standard constructor
+     ** @param numBins   Number of bins
+     ** @param minValue  Lower edge of histogram range
+     ** @param maxValue  Upper edge of histogram range
+     ** @param name      Histogram name
+     **/
+    Histo1D(uint32_t numBins, double minValue, double maxValue, const std::string& name = "");
+
+
+    /** @brief Destructor **/
+    virtual ~Histo1D() = default;
+
+
+    /** @brief Add an entry to the histogram
+     ** @param value  Value to be histogrammed
+     ** @param weight Weight for the entry
+     **/
+    void Add(double value, double weight = 1.);
+
+
+    /** @brief Clear histogram contents **/
+    void Clear();
+
+
+    /** @brief Histogram content in a bin
+     ** @param bin Bin index
+     ** @return Histogram content. Zero if bin index is out of scope.
+     **/
+    double Content(uint32_t bin) const;
+
+
+    /** @brief Number of entries **/
+    double NumEntries() const { return fNumEntries; }
+
+
+    /** @brief Upper edge **/
+    double MaxValue() const { return fMaxValue; }
+
+
+    /** @brief First moment of distribution **/
+    double Mean() const;
+
+
+    /** @brief Lower edge **/
+    double MinValue() const { return fMinValue; }
+
+
+    /** @brief Histogram name **/
+    const std::string& Name() const { return fName; }
+
+
+    /** @brief Number of bins **/
+    uint32_t NumBins() const { return fNumBins; }
+
+    /** @brief Add another histogram to an existing one
+     ** @param other Histogram object to be added
+     **/
+    Histo1D& operator+=(const Histo1D& other);
+
+
+    /** @brief Overflow **/
+    double Overflow() const { return fOverflow; }
+
+
+    /** @brief Second moment of distribution **/
+    double Stddev() const;
+
+
+    /** @brief Properties to string **/
+    std::string ToString() const;
+
+
+    /** @brief Underflow **/
+    double UnderFlow() const { return fUnderflow; }
+
+
+  private:
+    /** Properties **/
+    const std::string fName;
+    const uint32_t fNumBins;
+    const double fMinValue;
+    const double fMaxValue;
+
+    /** Content **/
+    std::vector<double> fContent = {};
+    double fUnderflow            = 0;
+    double fOverflow             = 0;
+    size_t fNumEntries           = 0;
+  };
+
+
+  /** @brief Adding two histograms
+   ** @param h1,h2 Histograms to be added
+   ** @return Sum histogram
+   **/
+  inline Histo1D operator+(Histo1D h1, const Histo1D& h2)
+  {
+    h1 += h2;
+    return h1;
+  }
+
+} /* namespace cbm::algo */
+
+#endif /* ALGO_QA_HISTO1D_H */
diff --git a/reco/tasks/CMakeLists.txt b/reco/tasks/CMakeLists.txt
index 6fdef364b6..fe504a8d37 100644
--- a/reco/tasks/CMakeLists.txt
+++ b/reco/tasks/CMakeLists.txt
@@ -3,6 +3,7 @@
 
 set(INCLUDE_DIRECTORIES
   ${CMAKE_CURRENT_SOURCE_DIR}
+  ${CBMROOT_SOURCE_DIR}
   )
 
 set(SRCS
diff --git a/reco/tasks/CbmTaskDigiEventQa.cxx b/reco/tasks/CbmTaskDigiEventQa.cxx
index e7218cd3c0..87fe026458 100644
--- a/reco/tasks/CbmTaskDigiEventQa.cxx
+++ b/reco/tasks/CbmTaskDigiEventQa.cxx
@@ -4,23 +4,28 @@
 
 #include "CbmTaskDigiEventQa.h"
 
-#include "CbmDigiBranchBase.h"
-#include "CbmDigiManager.h"
-#include "CbmDigiTimeslice.h"
-#include "CbmModuleList.h"
 #include "CbmReco.h"  // for CbmRecoConfig
 
 #include <FairRunOnline.h>
 #include <Logger.h>
 
-#include <TH1F.h>
+#include <TH1D.h>
 #include <THttpServer.h>
 #include <TStopwatch.h>
 
 #include <cassert>
 #include <iomanip>
 
+#include "algo/qa/Histo1D.h"
+
 using namespace std;
+using cbm::algo::qa::DigiEventQaConfig;
+using cbm::algo::qa::DigiEventQaData;
+using cbm::algo::qa::DigiEventQaDetConfig;
+
+#define NUM_BINS 100
+#define BORDER 10.
+
 
 // -----   Constructor   -----------------------------------------------------
 CbmTaskDigiEventQa::CbmTaskDigiEventQa() : FairTask("DigiEventQa") {}
@@ -35,11 +40,24 @@ CbmTaskDigiEventQa::~CbmTaskDigiEventQa() {}
 // -----   Configuration   ---------------------------------------------------
 void CbmTaskDigiEventQa::Config(const CbmRecoConfig& config)
 {
-  auto limitIt = config.fEvtbuildWindows.find(ECbmModuleId::kSts);
-  if (limitIt != config.fEvtbuildWindows.end()) {
-    double lower     = limitIt->second.first - 10.;
-    double upper     = limitIt->second.second + 10.;
-    fHistDigiTimeSts = new TH1F("hDigiTimeSts", "STS digi time in event", 100, lower, upper);
+
+  // The histogram ranges are defined by the event building windows. The number of bins
+  // is hard-coded. To be changed on request.
+
+  for (const auto& entry : config.fEvtbuildWindows) {
+    ECbmModuleId system = entry.first;
+
+    // --- Create histogram
+    string name  = "hDigiTime" + ToString(entry.first);
+    string title = ToString(entry.first) + " digi time in event";
+    double lower = entry.second.first - BORDER;   // Lower edge of histogram
+    double upper = entry.second.second + BORDER;  // Upper edge of histogram
+    assert(fDigiTimeHistos.count(system) == 0);
+    fDigiTimeHistos[system] = new TH1D(name.c_str(), title.c_str(), NUM_BINS, lower, upper);
+
+    // --- Set algo configuration
+    assert(fConfig.fData.count(system) == 0);
+    fConfig.fData[system] = {NUM_BINS, lower, upper};
   }
 }
 // ---------------------------------------------------------------------------
@@ -52,55 +70,32 @@ void CbmTaskDigiEventQa::Exec(Option_t*)
   // --- Timer and counters
   TStopwatch timer;
   timer.Start();
-  double sumT          = 0.;
-  double sumTsq        = 0.;
-  double sumNumDigis   = 0.;
-  double sumNumDigisSq = 0.;
-  size_t numEvents     = 0;
-  size_t numDigis      = 0;
-
-  // --- Event loop
-  for (const auto& event : *fEvents) {
-    numEvents++;
-    double numDigisEvt = double(event.fData.fSts.fDigis.size());
-    sumNumDigis += numDigisEvt;
-    sumNumDigisSq += numDigisEvt * numDigisEvt;
-    for (const auto& digi : event.fData.fSts.fDigis) {
-      numDigis++;
-      const double tDigi = digi.GetTime() - event.fTime;
-      if (fHistDigiTimeSts) fHistDigiTimeSts->Fill(tDigi);
-      sumT += tDigi;
-      sumTsq += tDigi * tDigi;
-    }
-  }
 
-  // --- First and second moments of digi times and digi numbers
-  double tMean          = sumT / double(numDigis);
-  double tSqMean        = sumTsq / double(numDigis);
-  double tRms           = sqrt(tSqMean - tMean * tMean);
-  double numDigisMean   = sumNumDigis / double(numEvents);
-  double numDigisSqMean = sumNumDigisSq / double(numEvents);
-  double numDigisRms    = sqrt(numDigisSqMean - numDigisMean * numDigisMean);
+  // --- Algo execution
+  DigiEventQaConfig config;
+  DigiEventQaData result = fAlgo(*fEvents, fConfig);
+
+  // --- Copy QA results (Histo1D) into ROOT output histograms
+  // TODO: Probably not the most efficient implementation. Creates first a ROOT histogram from a CBM one (data copy),
+  // which is then added to the member histogram (another data copy). Should implement a method for direct addition TH1D + Histo1D.
+  for (const auto& entry : result.fDigiTimeHistos) {
+    ECbmModuleId subsystem = entry.first;
+    fDigiTimeHistos[subsystem]->Add(ToTH1D(entry.second));
+  }
 
   // --- Timeslice log
   timer.Stop();
   fExecTime += timer.RealTime();
+  size_t numEvents = result.fNumEvents;
   stringstream logOut;
   logOut << setw(15) << left << GetName() << " [";
   logOut << fixed << setw(8) << setprecision(1) << right << timer.RealTime() * 1000. << " ms] ";
-  logOut << "TS " << fNumTs << ", events " << numEvents << ", digis " << numDigis << ", digis per event ("
-         << numDigisMean << " +- " << numDigisRms << ")"
-         << ",  digi time (" << tMean << " +- " << tRms << ") ns";
+  logOut << "TS " << fNumTs << ", events " << numEvents;
   LOG(info) << logOut.str();
 
   // --- Run statistics
   fNumTs++;
   fNumEvents += numEvents;
-  fNumDigis += numDigis;
-
-
-  // See: macro/run/run_unpack_online.C
-  // See : recp/detectors/sts/unpack/CbmStsUnpackMonitor
 }
 // ----------------------------------------------------------------------------
 
@@ -112,13 +107,15 @@ void CbmTaskDigiEventQa::Finish()
   LOG(info) << GetName() << ": Run summary";
   LOG(info) << "Timeslices : " << fNumTs;
   LOG(info) << "Events     : " << fNumEvents;
-  LOG(info) << "Digis      : " << fNumDigis;
   LOG(info) << "Exec time  : " << fixed << setprecision(2) << 1000. * fExecTime / double(fNumTs) << " ms / TS";
-  if (fHistDigiTimeSts)
-    LOG(info) << "STS digi times : entries " << fHistDigiTimeSts->GetEntries() << ", mean "
-              << fHistDigiTimeSts->GetMean() << ", stddev " << fHistDigiTimeSts->GetStdDev();
+  for (const auto& entry : fDigiTimeHistos) {
+    ECbmModuleId subsystem = entry.first;
+    TH1D* histo            = entry.second;
+    LOG(info) << ToString(subsystem) << " digi times: entries " << histo->GetEntries() << ", mean " << histo->GetMean()
+              << ", stddev " << histo->GetStdDev();
+    histo->Write();
+  }
   LOG(info) << "=====================================";
-  if (fHistDigiTimeSts) fHistDigiTimeSts->Write();
 }
 // ----------------------------------------------------------------------------
 
@@ -145,11 +142,14 @@ InitStatus CbmTaskDigiEventQa::Init()
   THttpServer* server = FairRunOnline::Instance()->GetHttpServer();
   if (server) {
     LOG(info) << "--- Http server present; registering histograms";
-    if (fHistDigiTimeSts) server->Register("DigiEvent", fHistDigiTimeSts);
+    for (const auto& entry : fDigiTimeHistos)
+      server->Register("DigiEvent", entry.second);
   }
   else
     LOG(info) << "--- No Http server present";
 
+  // --- Log configuration
+  LOG(info) << "--- Algo configuration for " << fConfig.fData.size() << " subsystem(s): " << fConfig.ToString();
   LOG(info) << "==================================================";
 
   return kSUCCESS;
@@ -157,4 +157,18 @@ InitStatus CbmTaskDigiEventQa::Init()
 // ----------------------------------------------------------------------------
 
 
+// -----   Convert CBM histogram to ROOT histogram   --------------------------
+TH1D* CbmTaskDigiEventQa::ToTH1D(const cbm::algo::Histo1D& source)
+{
+
+  TH1D* result =
+    new TH1D(source.Name().c_str(), source.Name().c_str(), source.NumBins(), source.MinValue(), source.MaxValue());
+  for (uint32_t bin = 0; bin < source.NumBins(); bin++) {
+    result->SetBinContent(bin, source.Content(bin));
+  }
+  result->SetEntries(source.NumEntries());
+  return result;
+}
+// ----------------------------------------------------------------------------
+
 ClassImp(CbmTaskDigiEventQa)
diff --git a/reco/tasks/CbmTaskDigiEventQa.h b/reco/tasks/CbmTaskDigiEventQa.h
index 32452a9133..1710fb6b85 100644
--- a/reco/tasks/CbmTaskDigiEventQa.h
+++ b/reco/tasks/CbmTaskDigiEventQa.h
@@ -12,19 +12,20 @@
 
 #include <vector>
 
-//#include "EventBuilder.h"
+#include "algo/qa/DigiEventQa.h"
+#include "algo/qa/Histo1D.h"
+
 
-class CbmDigiManager;
 class CbmRecoConfig;
-class TH1F;
+class TH1D;
 
 /** @class CbmTaskDigiEventQa
  ** @brief QA task class for digi events produced by the event builder
  ** @author Volker Friese <v.friese@gsi.de>
  ** @since 15.03.2022
  **
- ** Currently implemented functionality: histogram the STS digi time within each event.
- ** To be expanded for other detector systems and probably more QA figures.
+ ** Currently implemented functionality: histogram the digi time within each event.
+ ** To be expanded for more QA figures.
  ** The histograms are published to the THttpServer.
  **/
 class CbmTaskDigiEventQa : public FairTask {
@@ -67,13 +68,25 @@ private:  // methods
   virtual InitStatus Init();
 
 
+  /** @brief Create a ROOT TH1D from a Histo1D object
+   ** @param Source histogram
+   ** @param ROOT histogram
+   */
+  TH1D* ToTH1D(const cbm::algo::Histo1D& source);
+
+
 private:                                               // members
   const std::vector<CbmDigiEvent>* fEvents = nullptr;  //! Input data (events)
   size_t fNumTs                            = 0;        ///< Number of processed timeslices
   size_t fNumEvents                        = 0;        ///< Number of analysed events
   size_t fNumDigis                         = 0;        ///< Number of analysed digis
   double fExecTime                         = 0.;       ///< Execution time [s]
-  TH1F* fHistDigiTimeSts                   = nullptr;  //!
+
+  cbm::algo::qa::DigiEventQa fAlgo         = {};
+  cbm::algo::qa::DigiEventQaConfig fConfig = {};
+
+  // ---- Histograms with digi times
+  std::map<ECbmModuleId, TH1D*> fDigiTimeHistos = {};
 
 
   ClassDef(CbmTaskDigiEventQa, 1);
-- 
GitLab