cbmreco: Hit reconstruction in digi-events

algo/base/Options.cxx

@@ -93,6 +93,7 @@ Options::Options(int argc, char** argv)
     ("compress-archive", po::bool_switch(&fCompressArchive)->default_value(false), "Enable compression for output archives")
     ("steps", po::value(&fRecoSteps)->multitoken()->default_value({Step::Unpack, Step::DigiTrigger, Step::LocalReco, Step::Tracking})->value_name("<steps>"),
       "space separated list of reconstruction steps (unpack, digitrigger, localreco, ...)")
+    ("event-reco", po::bool_switch(&fReconstructDigiEvents)->default_value(false), "runs digi event reconstruction (local reco, tracking, trigger)")
     ("systems,s", po::value(&fDetectors)->multitoken()->default_value({Subsystem::STS, Subsystem::TOF, Subsystem::BMON, Subsystem::MUCH, Subsystem::RICH, Subsystem::TRD, Subsystem::TRD2D})->value_name("<detectors>"),
       "space separated list of detectors to process (sts, mvd, ...)")
     ("child-id,c", po::value(&fChildId)->default_value("00")->value_name("<id>"), "online process id on node")

algo/base/Options.h

@@ -68,6 +68,8 @@ namespace cbm::algo
 
     bool Has(QaStep qastep) const;
 
+    bool ReconstructDigiEvents() const { return fReconstructDigiEvents; }
+
    private:                  // members
     std::string fParamsDir;  // TODO: can we make this a std::path?
     std::string fInputLocator;
@@ -97,6 +99,7 @@ namespace cbm::algo
     uint64_t fRunId        = 2391;
     uint64_t fRunStartTime = 0;
     bool fCollectAuxData   = false;
+    bool fReconstructDigiEvents = false;
   };
 
 }  // namespace cbm::algo

algo/global/Reco.cxx

@@ -241,6 +241,15 @@ void Reco::Init(const Options& opts)
     fTrdHitfind  = std::make_unique<trd::Hitfind>(setup, setup2d);
   }
 
+  // Digi event reconstruction:
+  {
+    fbReconstructDigiEvents = Opts().ReconstructDigiEvents();
+    // It makes no sence to reconstruct an event, if there is no STS, TRD or TOF
+    fbReconstructDigiEvents &= Opts().Has(fles::Subsystem::STS);
+    fbReconstructDigiEvents &= Opts().Has(fles::Subsystem::TRD);
+    fbReconstructDigiEvents &= Opts().Has(fles::Subsystem::TOF);
+  }
+
   // Tracking
   if (Opts().Has(Step::Tracking)) {
     if (fQaManager != nullptr && Opts().Has(QaStep::Tracking)) {
@@ -393,27 +402,16 @@ RecoResults Reco::Run(const fles::Timeslice& ts)
       QueueEvbuildMetrics(evbuildMonitor);
     }
 
-    // ***** DEBUG: BEGIN
-    if constexpr (0) {
-      int nEvents = events.size();
-      size_t nBmonHitsOneChannel{0};
-      size_t nBmonHitsTwoChannels{0};
-      for (int iE = 0; iE < nEvents; ++iE) {
-        const auto& event = events[iE];
-        // Calibrate TOF digis:
-        auto [bmonDigis, bmonCalMonitor]         = (*fBmonCalibrator)(event.fBmon);
-        auto [bmonHits, hitmonitor, digiIndices] = (*fBmonHitFinder)(bmonDigis);
-        if (fBmonHitFinderQa != nullptr) {
-          fBmonHitFinderQa->RegisterDigis(&bmonDigis);
-          fBmonHitFinderQa->RegisterHits(&bmonHits);
-          fBmonHitFinderQa->RegisterDigiIndices(&digiIndices);
-          fBmonHitFinderQa->Exec();
+    // --- Reconstruct and select digi events
+    if (Opts().ReconstructDigiEvents()) {
+      size_t nDiscardedEvents{0};
+      for (const auto& event : events) {
+        if (!ReconstructEvent(event)) {
+          ++nDiscardedEvents;
         }
       }
-      L_(info) << "!!!! BMON hits with two channels: " << nBmonHitsTwoChannels << " / "
-               << (nBmonHitsTwoChannels + nBmonHitsOneChannel);
+      L_(info) << "Rate of discarded events " << double(nDiscardedEvents) / events.size();
     }
-    // ***** DEBUG: END
 
     // --- Filter data for output
     if (Opts().HasOutput(RecoData::DigiTimeslice)) {
@@ -491,6 +489,45 @@ void Reco::PrintTimings(xpu::timings& timings)
   }
 }
 
+bool Reco::ReconstructEvent(const DigiEvent& digiEvent)
+{
+  RecoResults recoEvent;
+  //* STS hit reconstruction
+  {
+    auto stsResults = (*fStsHitFinder)(digiEvent.fSts);
+    if (stsResults.hits.NElements() < 2) {  // TODO: Provide a config for cuts (testing mode for now)
+      return false;
+    }
+    recoEvent.stsHits = stsResults.hits;
+  }
+
+  //* TOF hit reconstruction
+  {
+    auto [caldigis, calmonitor]          = (*fTofCalibrator)(digiEvent.fTof);
+    auto [hits, hitmonitor, digiindices] = (*fTofHitFinder)(caldigis);
+    if (hits.NElements() < 1) {  // TODO: Provide a config for cuts (testing mode for now)
+      return false;
+    }
+    recoEvent.tofHits = std::move(hits);
+  }
+
+  //* TRD hit reconstruction
+  {
+    // FIXME: additional copy of digis, figure out how to pass 1d + 2d digis at once to hitfinder
+    const auto& digis1d = digiEvent.fTrd;
+    const auto& digis2d = digiEvent.fTrd2d;
+    PODVector<CbmTrdDigi> allDigis{};
+    allDigis.reserve(digis1d.size() + digis2d.size());
+    std::copy(digis1d.begin(), digis1d.end(), std::back_inserter(allDigis));
+    std::copy(digis2d.begin(), digis2d.end(), std::back_inserter(allDigis));
+    auto trdResults   = (*fTrdHitfind)(allDigis);
+    recoEvent.trdHits = std::move(std::get<0>(trdResults));
+  }
+
+  return true;
+}
+
+
 template<class Unpacker>
 auto Reco::RunUnpacker(const std::unique_ptr<Unpacker>& unpacker, const fles::Timeslice& ts) -> UnpackResult_t<Unpacker>
 {

algo/global/Reco.h

@@ -130,13 +130,16 @@ namespace cbm::algo
 
     void Init(const Options&);
     RecoResults Run(const fles::Timeslice&);
+
+    bool ReconstructEvent(const DigiEvent& event);
     void Finalize();
     void PrintTimings(xpu::timings&);
 
     void QueueProcessingExtraMetrics(const ProcessingExtraMonitor&);
 
    private:
-    bool fInitialized = false;
+    bool fInitialized            = false;
+    bool fbReconstructDigiEvents = false;
     ChainContext fContext;
     xpu::timings fTimesliceTimesAcc;
     std::shared_ptr<HistogramSender> fSender;