From d9b51bdd8c6d09a4c876ce1aa2a5654b5481672d Mon Sep 17 00:00:00 2001
From: Shreya Roy <shreyaroy@jcbose.ac.in>
Date: Thu, 26 Jan 2023 11:07:07 +0000
Subject: [PATCH] New class for digi event selector and corresponding GTest

---
 algo/CMakeLists.txt                   |   3 +
 algo/detectors/tof/TofConfig.h        |  50 +++++++++++++
 algo/evselector/DigiEventSelector.cxx | 100 ++++++++++++++++++++++++++
 algo/evselector/DigiEventSelector.h   |  45 ++++++++++++
 algo/test/CMakeLists.txt              |  13 ++--
 algo/test/_GTestDigiEventSelector.cxx |  95 ++++++++++++++++++++++++
 reco/tasks/CbmTaskBuildEvents.cxx     |  12 ++++
 reco/tasks/CbmTaskBuildEvents.h       |  44 ++++++++----
 8 files changed, 346 insertions(+), 16 deletions(-)
 create mode 100644 algo/detectors/tof/TofConfig.h
 create mode 100644 algo/evselector/DigiEventSelector.cxx
 create mode 100644 algo/evselector/DigiEventSelector.h
 create mode 100644 algo/test/_GTestDigiEventSelector.cxx

diff --git a/algo/CMakeLists.txt b/algo/CMakeLists.txt
index 3ad6df5dee..5f807e9cbf 100644
--- a/algo/CMakeLists.txt
+++ b/algo/CMakeLists.txt
@@ -5,6 +5,7 @@ add_subdirectory(test)
 set(SRCS
   evbuild/EventBuilder.cxx
   trigger/TimeClusterTrigger.cxx
+  evselector/DigiEventSelector.cxx
   detectors/sts/StsReadoutConfig.cxx
   detectors/sts/UnpackSts.cxx
   detectors/much/MuchReadoutConfig.cxx
@@ -17,8 +18,10 @@ target_include_directories(Algo
   PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}
          ${CMAKE_CURRENT_SOURCE_DIR}/evbuild
          ${CMAKE_CURRENT_SOURCE_DIR}/trigger
+         ${CMAKE_CURRENT_SOURCE_DIR}/evselector
          ${CMAKE_CURRENT_SOURCE_DIR}/detectors/sts
          ${CMAKE_CURRENT_SOURCE_DIR}/detectors/much
+         ${CMAKE_CURRENT_SOURCE_DIR}/detectors/tof
  )
 
 target_link_libraries(Algo PUBLIC OnlineData INTERFACE FairLogger::FairLogger external::fles_ipc)
diff --git a/algo/detectors/tof/TofConfig.h b/algo/detectors/tof/TofConfig.h
new file mode 100644
index 0000000000..01d888152f
--- /dev/null
+++ b/algo/detectors/tof/TofConfig.h
@@ -0,0 +1,50 @@
+/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Dominik Smith [committer] */
+
+#ifndef CBM_ALGO_TOFCONFIG
+#define CBM_ALGO_TOFCONFIG 1
+
+#include <cstdint>
+
+namespace cbm::algo
+{
+  /** @class TofConfig
+  ** @author Dominik Smith <d.smith@gsi.de>
+ 
+** @since 2022
+ ** @brief
+ **
+ **/
+  class TofConfig {
+
+  public:
+    static uint32_t GetTofTrackingStation(const uint32_t smType, const uint32_t sm, const uint32_t rpc)
+    {
+      const uint8_t numSmTypes                 = 10;
+      const uint8_t numRpc[numSmTypes]         = {5, 3, 5, 1, 1, 1, 2, 2, 1, 2};
+      const uint8_t numSm[numSmTypes]          = {5, 0, 1, 0, 0, 1, 1, 1, 2, 1};
+      const uint8_t trkStation[numSmTypes][25] = {
+        {0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 3, 3, 3, 3, 3, 0, 0, 0, 0, 0, 3, 3, 3, 3, 3},
+        {0, 0, 0, 0, 0},
+        {2, 2, 2, 2, 2},
+        {0, 0, 0, 0, 0},
+        {0, 0, 0, 0, 0},
+        {0, 0, 0, 0, 0},
+        {1, 1, 0, 0, 0},
+        {0, 0, 0, 0, 0},
+        {0, 0, 0, 0, 0},
+        {1, 1, 2, 2, 0}};
+
+      if (smType < numSmTypes) {
+        if (sm < numSm[smType] && rpc < numRpc[smType]) return trkStation[smType][sm * numRpc[smType] + rpc];
+        else
+          return 0;
+      }
+      else
+        return 0;
+    }
+  };
+}  // namespace cbm::algo
+
+#endif /* CBM_ALGO_TOFCONFIG */
diff --git a/algo/evselector/DigiEventSelector.cxx b/algo/evselector/DigiEventSelector.cxx
new file mode 100644
index 0000000000..813790422b
--- /dev/null
+++ b/algo/evselector/DigiEventSelector.cxx
@@ -0,0 +1,100 @@
+/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Shreya Roy [committer], Pierre-Alain Loizeau, Norbert Herrmann, Volker Friese, Dominik Smith */
+
+#include "DigiEventSelector.h"
+
+#include "CbmStsDigi.h"
+
+#include <Logger.h>
+
+#include "TofConfig.h"
+
+#include <iterator>
+#include <map>
+#include <unordered_set>
+
+
+namespace cbm::algo
+{
+
+  bool DigiEventSelector::operator()(const CbmDigiEvent& event) const
+  {
+    //TO DO: Investigate possible algorithms that use bitset instead of std::unordered_map.
+    //This requires consequtive addresses for modules and stations.
+
+    if (0 < fParams.fStsMinStations) {
+      if (!CheckSts(event.fData.fSts.fDigis)) return false;
+    }
+    if (0 < fParams.fTofMinLayers) {
+      if (!CheckTof(event.fData.fTof.fDigis)) return false;
+    }
+    return true;
+  }
+
+  bool DigiEventSelector::CheckSts(const std::vector<CbmStsDigi>& digis) const
+  {
+    const uint16_t chanPerSide = 1024;
+    /// check for requested number of different STS stations
+    std::unordered_set<uint32_t> setStations;
+    std::unordered_map<int32_t, bool> mModules;
+
+    for (auto& digi : digis) {
+      const int32_t addr = digi.GetAddress();
+      auto itModule      = mModules.find(addr);
+      if (itModule == mModules.end()) {
+        mModules[addr] = digi.GetChannel() / chanPerSide;  // = 0,1 for sides
+      }
+      else {
+        if (digi.GetChannel() / chanPerSide != itModule->second) {
+          /// Found other side => non-zero cluster chance, insert into stations set
+          const uint32_t stationAddr = CbmStsAddress::GetElementId(addr, EStsElementLevel::kStsUnit);
+          auto itStation             = setStations.find(stationAddr);
+          if (itStation == setStations.end()) {
+            setStations.insert(stationAddr);
+            if (setStations.size() == fParams.fStsMinStations) break;
+          }
+        }
+      }
+    }
+    if (setStations.size() < fParams.fStsMinStations) { return false; }
+    else {
+      return true;
+    }
+  }
+
+  bool DigiEventSelector::CheckTof(const std::vector<CbmTofDigi>& digis) const
+  {
+    /// check for requested number of different RPCs
+    std::unordered_set<int32_t> setRpcs;
+    std::unordered_set<int32_t> setTofStation;
+    std::unordered_map<int32_t, bool> mStrips;
+    for (auto& digi : digis) {
+      const int32_t addr      = digi.GetAddress();
+      const int32_t stripAddr = CbmTofAddress::GetStripFullId(addr);
+
+      auto itStrip = mStrips.find(stripAddr);
+      if (itStrip == mStrips.end()) { mStrips[stripAddr] = digi.GetSide(); }
+      else {
+        if (digi.GetSide() != itStrip->second) {
+          /// Found other end => full strip, insert into counter set
+          const int32_t rpcAddr = CbmTofAddress::GetRpcFullId(addr);
+          auto itRpc            = setRpcs.find(rpcAddr);
+          if (itRpc == setRpcs.end()) {
+            const int32_t smId         = CbmTofAddress::GetSmId(addr);
+            const int32_t smType       = CbmTofAddress::GetSmType(addr);
+            const int32_t rpcId        = CbmTofAddress::GetRpcId(addr);
+            const int32_t TofStationId = TofConfig::GetTofTrackingStation(smType, smId, rpcId);
+            setTofStation.insert(TofStationId);
+            if (setTofStation.size() == fParams.fTofMinLayers) break;
+          }
+        }
+      }
+    }
+    if (setTofStation.size() < fParams.fTofMinLayers) { return false; }
+    else {
+      return true;
+    }
+  }
+
+}  // namespace cbm::algo
diff --git a/algo/evselector/DigiEventSelector.h b/algo/evselector/DigiEventSelector.h
new file mode 100644
index 0000000000..f55b9fb48c
--- /dev/null
+++ b/algo/evselector/DigiEventSelector.h
@@ -0,0 +1,45 @@
+/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Shreya Roy [committer], Pierre-Alain Loizeau, Norbert Herrmann, Volker Friese, Dominik Smith */
+
+#ifndef CBM_ALGO_DIGIEVENTSELECTOR_H
+#define CBM_ALGO_DIGIEVENTSELECTOR_H 1
+
+#include "CbmDigiEvent.h"
+
+#include <cstdint>
+
+namespace cbm::algo
+{
+
+  struct DigiEventSelectorParams {
+    uint8_t fStsMinStations = 0;
+    uint8_t fTofMinLayers   = 0;
+  };
+
+  /** @class DigiEventSelector
+  ** @author Dominik Smith <d.smith@gsi.de>
+ ** @author Shreya Roy <s.roy@gsi.de>
+** @author Volker Friese <v.friese@gsi.de>
+ 
+** @since 2022
+ ** @brief
+ **
+ **/
+  class DigiEventSelector {
+
+  public:
+    DigiEventSelector(DigiEventSelectorParams params) { fParams = params; }
+    bool operator()(const CbmDigiEvent& event) const;
+
+    void SetParams(DigiEventSelectorParams params) { fParams = params; }
+
+  private:
+    DigiEventSelectorParams fParams;
+
+    bool CheckSts(const std::vector<CbmStsDigi>& digis) const;
+    bool CheckTof(const std::vector<CbmTofDigi>& digis) const;
+  };
+}  // namespace cbm::algo
+
+#endif /* CBM_ALGO_DIGIEVENTSELECTOR_H */
diff --git a/algo/test/CMakeLists.txt b/algo/test/CMakeLists.txt
index 9a40f30d8d..4cd6184ace 100644
--- a/algo/test/CMakeLists.txt
+++ b/algo/test/CMakeLists.txt
@@ -22,21 +22,26 @@ set(PVT_DEPS
   Gtest
   GtestMain
   )
-  
-  
+
+
 
 Set(TimeClusterTriggerSources
-  _GTestTimeClusterTrigger.cxx 
+  _GTestTimeClusterTrigger.cxx
 )
 
 CreateGTestExeAndAddTest(_GTestTimeClusterTrigger "${INCLUDE_DIRECTORIES}" "${LINK_DIRECTORIES}"
                          "${TimeClusterTriggerSources}" "${PUB_DEPS}" "${PVT_DEPS}" "${INT_DEPS}" "")
 
 Set(EventBuilderSources
-  _GTestEventBuilder.cxx 
+  _GTestEventBuilder.cxx
 )
 
 CreateGTestExeAndAddTest(_GTestEventBuilder "${INCLUDE_DIRECTORIES}" "${LINK_DIRECTORIES}"
                          "${EventBuilderSources}" "${PUB_DEPS}" "${PVT_DEPS}" "${INT_DEPS}" "")
 
+Set(DigiEventSelectorSources
+  _GTestDigiEventSelector.cxx
+)
 
+CreateGTestExeAndAddTest(_GTestDigiEventSelector "${INCLUDE_DIRECTORIES}" "${LINK_DIRECTORIES}"
+                         "${DigiEventSelectorSources}" "${PUB_DEPS}" "${PVT_DEPS}" "${INT_DEPS}" "")
diff --git a/algo/test/_GTestDigiEventSelector.cxx b/algo/test/_GTestDigiEventSelector.cxx
new file mode 100644
index 0000000000..8b528f6a93
--- /dev/null
+++ b/algo/test/_GTestDigiEventSelector.cxx
@@ -0,0 +1,95 @@
+/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Dominik Smith [committer] */
+
+#include "TofConfig.h"
+
+#include <unordered_set>
+
+#include "DigiEventSelector.h"
+#include "gtest/gtest-spi.h"
+#include "gtest/gtest.h"
+
+TEST(_GTestDigiEventSelector, CheckDigiEventSelectorAlgorithmSimple)
+{
+  SCOPED_TRACE("CheckDigiEventSelectorAlgorithSimple");
+
+  //Test STS
+  {
+    const uint maxStsStations = 8;
+    const uint maxStsModules  = 12;
+    const uint maxStsLadders  = 4;
+
+    for (uint numStations = 1; numStations < maxStsStations; numStations++) {
+      //Prepare input
+      CbmDigiEvent eventIn;
+      //Produce digi pairs with valid addresses
+      for (uint station = 0; station < numStations; station++) {
+        for (uint module = 0; module < maxStsModules; module++) {
+          for (uint ladder = 0; ladder < maxStsLadders; ladder++) {
+            for (uint halfladder = 0; halfladder <= 1; halfladder++) {
+              //add digis pairs
+              int32_t address = CbmStsAddress::GetAddress(station, ladder, halfladder, module);
+              eventIn.fData.fSts.fDigis.push_back(CbmStsDigi(address, 0, 0.0, 0.0));
+              eventIn.fData.fSts.fDigis.push_back(CbmStsDigi(address, 1024, 0.0, 0.0));  //other side channel
+
+              //add digis from next station without partner for intentionally failed test
+              int32_t nextAddress = CbmStsAddress::GetAddress(numStations, ladder, 0, module);
+              eventIn.fData.fSts.fDigis.push_back(CbmStsDigi(nextAddress, 1024, 0.0, 0.0));
+            }
+          }
+        }
+      }
+
+      //Prepare event filter
+      cbm::algo::DigiEventSelectorParams filterParams;
+      filterParams.fStsMinStations = numStations;
+      cbm::algo::DigiEventSelector evfilter(filterParams);
+
+      EXPECT_EQ(evfilter(eventIn), true);
+
+      //Test if digi without partner is properly disregarded
+      filterParams.fStsMinStations = numStations + 1;
+      evfilter.SetParams(filterParams);
+
+      EXPECT_EQ(evfilter(eventIn), false);
+    }
+  }
+
+  //Test TOF
+  {
+    //Prepare input
+    CbmDigiEvent eventIn;
+
+    const uint8_t numSmTypes         = 10;
+    const uint8_t numRpc[numSmTypes] = {5, 3, 5, 1, 1, 1, 2, 2, 1, 2};
+    const uint8_t numSm[numSmTypes]  = {5, 0, 1, 0, 0, 1, 1, 1, 2, 1};
+    std::unordered_set<int32_t> setTofStation;
+
+    for (uint smType = 0; smType < numSmTypes; smType++) {
+      for (uint sm = 0; sm < numSm[smType]; sm++) {
+        for (uint rpc = 0; rpc < numRpc[smType]; rpc++) {
+
+          uint32_t addrFront = CbmTofAddress::GetUniqueAddress(sm, rpc, 0, 0, smType, 0);
+          uint32_t addrBack  = CbmTofAddress::GetUniqueAddress(sm, rpc, 0, 1, smType, 0);
+          eventIn.fData.fTof.fDigis.push_back(CbmTofDigi(addrFront, 0.0, 0.0));
+          eventIn.fData.fTof.fDigis.push_back(CbmTofDigi(addrBack, 0.0, 0.0));
+
+          const int32_t TofStationId = cbm::algo::TofConfig::GetTofTrackingStation(smType, sm, rpc);
+          setTofStation.insert(TofStationId);
+
+          cbm::algo::DigiEventSelectorParams filterParams;
+          filterParams.fTofMinLayers = setTofStation.size();
+          cbm::algo::DigiEventSelector evfilter(filterParams);
+
+          EXPECT_EQ(evfilter(eventIn), true);
+
+          filterParams.fTofMinLayers = setTofStation.size() + 1;
+          evfilter.SetParams(filterParams);
+
+          EXPECT_EQ(evfilter(eventIn), false);
+        }
+      }
+    }
+  }
+}
diff --git a/reco/tasks/CbmTaskBuildEvents.cxx b/reco/tasks/CbmTaskBuildEvents.cxx
index 1610d1a569..e752197e02 100644
--- a/reco/tasks/CbmTaskBuildEvents.cxx
+++ b/reco/tasks/CbmTaskBuildEvents.cxx
@@ -135,6 +135,16 @@ void CbmTaskBuildEvents::Exec(Option_t*)
     fTimeBuildEvt += timerStep.RealTime();
   }
 
+  // Apply event selector if desired
+  if (fSelector) {
+    timerStep.Start();
+    auto noTrigger = [&](CbmDigiEvent& ev) { return !(*fSelector)(ev); };
+    auto removeIt  = std::remove_if(fEvents->begin(), fEvents->end(), noTrigger);
+    fEvents->erase(removeIt, fEvents->end());
+    timerStep.Stop();
+    fTimeSelectorEvt += timerStep.RealTime();
+  }
+
   // --- Timeslice statistics
   size_t numTriggers = fTriggers->size();
   size_t numEvents   = fEvents->size();
@@ -195,6 +205,8 @@ void CbmTaskBuildEvents::Finish()
             << " ms = " << 100. * fTimeFillTs / fTimeTot << " %";
   LOG(info) << "Time build events    : " << fixed << setprecision(2) << 1000. * fTimeBuildEvt / double(fNumTs)
             << " ms = " << 100. * fTimeBuildEvt / fTimeTot << " %";
+  LOG(info) << "Time selector events   : " << fixed << setprecision(2) << 1000. * fTimeSelectorEvt / double(fNumTs)
+            << " ms = " << 100. * fTimeSelectorEvt / fTimeTot << " %";
   LOG(info) << "=====================================";
 }
 // ----------------------------------------------------------------------------
diff --git a/reco/tasks/CbmTaskBuildEvents.h b/reco/tasks/CbmTaskBuildEvents.h
index 6cb231842e..056d3d1006 100644
--- a/reco/tasks/CbmTaskBuildEvents.h
+++ b/reco/tasks/CbmTaskBuildEvents.h
@@ -12,6 +12,7 @@
 
 #include <vector>
 
+#include "DigiEventSelector.h"
 #include "EventBuilder.h"
 
 class CbmDigiManager;
@@ -64,6 +65,23 @@ public:
     fEventWindows[system] = std::make_pair(tMin, tMax);
   }
 
+  /** @brief Activate event selector which requires a minimum number of fired layers
+   ** @param params Struct with minimum number of layers for different detectors
+   **/
+
+  void SetDigiEventSelector(cbm::algo::DigiEventSelectorParams params)
+  {
+    if (fSelector == nullptr) {
+      // New selector, straightforward
+      fSelector = std::make_unique<cbm::algo::DigiEventSelector>(params);
+    }
+    else {
+      // Re-use existing selector as functor without state
+      fSelector->SetParams(params);
+    }
+  }
+
+
 private:  // methods
   /** @brief Task initialisation **/
   virtual InitStatus Init();
@@ -80,25 +98,27 @@ private:  // methods
    **/
   size_t GetNumDigis(const CbmDigiData& data, ECbmModuleId system);
 
+private:                                                    // members
+  const CbmDigiTimeslice* fTimeslice   = nullptr;           //! Input data (from unpacking)
+  CbmDigiManager* fDigiMan             = nullptr;           //! Input data (from simulation)
+  const std::vector<double>* fTriggers = nullptr;           //! Input data (triggers)
+  std::vector<CbmDigiEvent>* fEvents   = nullptr;           //! Output data (events)
+  std::unique_ptr<cbm::algo::DigiEventSelector> fSelector;  //! Event selector
 
-private:                                           // members
-  const CbmDigiTimeslice* fTimeslice   = nullptr;  //! Input data (from unpacking)
-  CbmDigiManager* fDigiMan             = nullptr;  //! Input data (from simulation)
-  const std::vector<double>* fTriggers = nullptr;  //! Input data (triggers)
-  std::vector<CbmDigiEvent>* fEvents   = nullptr;  //! Output data (events)
-  cbm::algo::EventBuilder fAlgo {};                //! Algorithm
+  cbm::algo::EventBuilder fAlgo {};  //! Algorithm
 
   std::map<ECbmModuleId, std::pair<double, double>> fEventWindows;
 
   // for diagnostics
   std::map<ECbmModuleId, size_t> fNumDigisTs;  //  Number of digis in timeslices
   std::map<ECbmModuleId, size_t> fNumDigisEv;  //  Number of digis in events
-  size_t fNumTs        = 0;                    //  Number of processed time slices
-  size_t fNumTriggers  = 0;                    //  Number of triggers
-  size_t fNumEvents    = 0;                    //  Number of produced events
-  double fTimeFillTs   = 0.;
-  double fTimeBuildEvt = 0.;
-  double fTimeTot      = 0.;
+  size_t fNumTs           = 0;                 //  Number of processed time slices
+  size_t fNumTriggers     = 0;                 //  Number of triggers
+  size_t fNumEvents       = 0;                 //  Number of produced events
+  double fTimeFillTs      = 0.;
+  double fTimeBuildEvt    = 0.;
+  double fTimeSelectorEvt = 0.;
+  double fTimeTot         = 0.;
 
   ClassDef(CbmTaskBuildEvents, 1);
 };
-- 
GitLab