TPC: Update Kr cluster finder workflow

* smaller algorithm improvements in Kr cluster finder
* adapt workflow to new moving window in Kr cluster finder to be able to
  process full TFs
* some code cleanup

Author: wiechula

Detectors/TPC/reconstruction/include/TPCReconstruction/KrBoxClusterFinder.h

@@ -87,6 +87,7 @@
 #include <vector>
 #include <array>
 #include <deque>
+#include <gsl/span>
 
 namespace o2
 {
@@ -171,10 +172,12 @@ class KrBoxClusterFinder
     mMaxClusterSizeTime = maxClusterSizeTime;
   }
 
-  void fillADCValue(int cru, int rowInSector, int padInRow, int timeBin, float adcValue);
-  void resetADCMap();
+  void loopOverSector(const gsl::span<const Digit> eventSector, const int sector);
 
-  void loopOverSector(std::vector<o2::tpc::Digit>& eventSector, const int sector);
+  void loopOverSector(const std::vector<Digit>& eventSector, const int sector)
+  {
+    loopOverSector(gsl::span(eventSector.data(), eventSector.size()), sector);
+  }
 
  private:
   // These variables can be varied
@@ -243,10 +246,17 @@ class KrBoxClusterFinder
   /// Time slice four is the interesting one. In there, local maxima are found and clusters are built from it. After it is processed, timeslice number 1 will be dropped and another timeslice will be put at the end of the set.
   std::deque<TimeSliceSector> mSetOfTimeSlices{};
 
-  void createInitialMap(std::vector<o2::tpc::Digit>& eventSector);
-  void popFirstTimeSliceFromMap();
+  /// count the number of ADC values in each slice which satisfy the condition for the minumum Qmax
+  std::deque<int> mNumADCwithMinQmax{};
+
+  void createInitialMap(const gsl::span<const Digit> eventSector);
+  void popFirstTimeSliceFromMap()
+  {
+    mSetOfTimeSlices.pop_front();
+    mNumADCwithMinQmax.pop_front();
+  }
   void fillADCValueInLastSlice(int cru, int rowInSector, int padInRow, float adcValue);
-  void addTimeSlice(std::vector<o2::tpc::Digit>& eventSector, const int timeSlice);
+  void addTimeSlice(const gsl::span<const Digit> eventSector, const int timeSlice);
 
   /// For each ROC, the maximum cluster size has to be chosen
   void setMaxClusterSize(int row);

Detectors/TPC/reconstruction/src/KrBoxClusterFinder.cxx

@@ -38,20 +38,16 @@ void KrBoxClusterFinder::loadGainMapFromFile(const std::string_view calDetFileNa
   LOGP(info, "Loaded gain map object '{}' from file '{}'", calDetFileName, gainMapName);
 }
 
-void KrBoxClusterFinder::createInitialMap(std::vector<o2::tpc::Digit>& eventSector)
+void KrBoxClusterFinder::createInitialMap(const gsl::span<const Digit> eventSector)
 {
   mSetOfTimeSlices.clear();
+  mNumADCwithMinQmax.clear();
 
   for (int iTimeSlice = 0; iTimeSlice <= 2 * mMaxClusterSizeTime; ++iTimeSlice) {
     addTimeSlice(eventSector, iTimeSlice);
   }
 }
 
-void KrBoxClusterFinder::popFirstTimeSliceFromMap()
-{
-  mSetOfTimeSlices.pop_front();
-}
-
 void KrBoxClusterFinder::fillADCValueInLastSlice(int cru, int rowInSector, int padInRow, float adcValue)
 {
   auto& timeSlice = mSetOfTimeSlices.back();
@@ -81,9 +77,10 @@ void KrBoxClusterFinder::fillADCValueInLastSlice(int cru, int rowInSector, int p
   timeSlice[rowInSector][corPad] = adcValue;
 }
 
-void KrBoxClusterFinder::addTimeSlice(std::vector<o2::tpc::Digit>& eventSector, const int timeSlice)
+void KrBoxClusterFinder::addTimeSlice(const gsl::span<const Digit> eventSector, const int timeSlice)
 {
   mSetOfTimeSlices.emplace_back();
+  auto& nADCminQmax = mNumADCwithMinQmax.emplace_back();
 
   for (; mFirstDigit < eventSector.size(); ++mFirstDigit) {
     const auto& digit = eventSector[mFirstDigit];
@@ -97,35 +94,34 @@ void KrBoxClusterFinder::addTimeSlice(std::vector<o2::tpc::Digit>& eventSector,
 
     const int rowInSector = digit.getRow();
     const int padInRow = digit.getPad();
-    float adcValue = digit.getChargeFloat();
+    const float adcValue = digit.getChargeFloat();
+    if (adcValue > mQThresholdMax) {
+      ++nADCminQmax;
+    }
 
     fillADCValueInLastSlice(cru, rowInSector, padInRow, adcValue);
   }
 }
 
-void KrBoxClusterFinder::loopOverSector(std::vector<o2::tpc::Digit>& eventSector, const int sector)
+void KrBoxClusterFinder::loopOverSector(const gsl::span<const Digit> eventSector, const int sector)
 {
   mFirstDigit = 0;
   mSector = sector;
 
   createInitialMap(eventSector);
   for (int iTimeSlice = mMaxClusterSizeTime; iTimeSlice < mMaxTimes - mMaxClusterSizeTime; ++iTimeSlice) {
-    findLocalMaxima(true, iTimeSlice);
+    // only search for a local maximum if the central time slice has at least one ADC above the charge threshold
+    if (mNumADCwithMinQmax[mMaxClusterSizeTime]) {
+      findLocalMaxima(true, iTimeSlice);
+    }
     popFirstTimeSliceFromMap();
     addTimeSlice(eventSector, iTimeSlice + mMaxClusterSizeTime + 1);
-  }
-}
 
-void KrBoxClusterFinder::resetADCMap()
-{
-  // Has to be reimplemented!
-  return;
-}
-
-void KrBoxClusterFinder::fillADCValue(int cru, int rowInSector, int padInRow, int timeBin, float adcValue)
-{
-  // Has to be reimplemented!
-  return;
+    // don't spend unnecessary time looping till mMaxTimes if there is no more data
+    if (mFirstDigit >= eventSector.size()) {
+      break;
+    }
+  }
 }
 
 void KrBoxClusterFinder::init()
@@ -323,6 +319,7 @@ std::vector<std::tuple<int, int, int>> KrBoxClusterFinder::findLocalMaxima(bool
           }
         }
       }
+
       if (!thisIsMax) {
         continue;
       } else {
@@ -338,7 +335,7 @@ std::vector<std::tuple<int, int, int>> KrBoxClusterFinder::findLocalMaxima(bool
       }
     }
   }
-  // }
+
   return localMaximaCoords;
 }
 

Detectors/TPC/workflow/include/TPCWorkflow/KryptonClustererSpec.h

@@ -18,7 +18,6 @@
 #define TPC_KryptonClustererSpec_H_
 
 #include "Framework/DataProcessorSpec.h"
-#include <string_view>
 
 namespace o2
 {
@@ -27,7 +26,7 @@ namespace tpc
 
 /// create a processor spec
 /// read simulated TPC clusters from file and publish
-o2::framework::DataProcessorSpec getKryptonClustererSpec(const std::string inputSpec, int ilane, std::vector<int> const& sectors);
+o2::framework::DataProcessorSpec getKryptonClustererSpec();
 
 } // end namespace tpc
 } // end namespace o2

Detectors/TPC/workflow/src/KryptonClustererSpec.cxx

@@ -10,32 +10,20 @@
 // or submit itself to any jurisdiction.
 
 #include <memory>
-#include <vector>
-#include <string>
-#include <algorithm>
-#include "fmt/format.h"
 
 #include "Framework/Task.h"
-#include "Framework/ControlService.h"
-#include "Framework/ConfigParamRegistry.h"
+#include "Framework/InputRecordWalker.h"
 #include "Framework/Logger.h"
 #include "Framework/DataProcessorSpec.h"
-#include "Framework/WorkflowSpec.h"
-
-#include "DataFormatsTPC/TPCSectorHeader.h"
 #include "Headers/DataHeader.h"
-#include "CCDB/CcdbApi.h"
-#include "DetectorsCalibration/Utils.h"
+#include "DataFormatsTPC/TPCSectorHeader.h"
 
-#include "TPCBase/RDHUtils.h"
-#include "TPCBase/Mapper.h"
 #include "TPCReconstruction/KrBoxClusterFinder.h"
 #include "TPCReconstruction/KrCluster.h"
-#include "TPCReconstruction/RawReaderCRU.h"
-#include "TPCWorkflow/CalibProcessingHelper.h"
 #include "TPCWorkflow/KryptonClustererSpec.h"
 
 using namespace o2::framework;
+using namespace o2::header;
 using SubSpecificationType = o2::framework::DataAllocator::SubSpecificationType;
 
 namespace o2
@@ -46,132 +34,69 @@ namespace tpc
 class KrBoxClusterFinderDevice : public o2::framework::Task
 {
  public:
-  KrBoxClusterFinderDevice(int lane, const std::vector<int>& sectors) : mClusters(36), mLane{lane}, mSectors(sectors), mClusterFinder{std::make_unique<KrBoxClusterFinder>()} {}
+  KrBoxClusterFinderDevice() : mClusterFinder{std::make_unique<KrBoxClusterFinder>()} {}
 
   void init(o2::framework::InitContext& ic) final
   {
-    // set up ADC value filling
-    mRawReader.createReader("");
-
-    mRawReader.setLinkZSCallback([this](int cru, int rowInSector, int padInRow, int timeBin, float adcValue) -> bool {
-      const int sector = cru / 10;
-      if ((mLastSector > -1) && (sector != mLastSector)) {
-        LOGP(debug, "analysing sector {} ({})", mLastSector, sector);
-        mClusterFinder->findLocalMaxima(true);
-        LOGP(info, "found {} clusters in sector {}", mClusterFinder->getClusters().size(), mLastSector);
-        std::swap(mClusters[mLastSector], mClusterFinder->getClusters());
-        mClusterFinder->resetADCMap();
-        mClusterFinder->resetClusters();
-      }
-
-      mClusterFinder->fillADCValue(cru, rowInSector, padInRow, timeBin, adcValue);
-
-      mLastSector = sector;
-      return true;
-    });
-
-    mMaxEvents = static_cast<uint32_t>(ic.options().get<int>("max-events"));
-    mForceQuit = ic.options().get<bool>("force-quit");
+    mClusterFinder->init();
   }
 
   void run(o2::framework::ProcessingContext& pc) final
   {
-    // in case the maximum number of events was reached don't do further processing
-    if (mReadyToQuit) {
-      return;
-    }
+    for (auto const& inputRef : InputRecordWalker(pc.inputs())) {
+      auto const* sectorHeader = DataRefUtils::getHeader<TPCSectorHeader*>(inputRef);
+      if (sectorHeader == nullptr) {
+        LOGP(error, "sector header missing on header stack for input on ", inputRef.spec->binding);
+        continue;
+      }
 
-    std::for_each(mClusters.begin(), mClusters.end(), [](auto& cl) { cl.clear(); });
+      const int sector = sectorHeader->sector();
+      auto inDigits = pc.inputs().get<gsl::span<o2::tpc::Digit>>(inputRef);
 
-    auto& reader = mRawReader.getReaders()[0];
-    mActiveSectors = calib_processing_helper::processRawData(pc.inputs(), reader, false, mSectors);
+      mClusterFinder->loopOverSector(inDigits, sector);
 
-    // analyse the final sector
-    if (mLastSector > -1) {
-      LOGP(debug, "analysing sector {}", mLastSector);
-      mClusterFinder->findLocalMaxima(true);
-      if (mClusterFinder->getClusters().size()) {
-        LOGP(info, "found {} clusters in sector {}", mClusterFinder->getClusters().size(), mLastSector);
-        std::swap(mClusters[mLastSector], mClusterFinder->getClusters());
-      }
-    }
-    mClusterFinder->resetADCMap();
-    mClusterFinder->resetClusters();
+      snapshotClusters(pc.outputs(), mClusterFinder->getClusters(), sector);
 
-    ++mProcessedTFs;
-    LOGP(info, "Number of processed time frames: {} ({})", mProcessedTFs, mMaxEvents);
-
-    snapshotClusters(pc.outputs());
-
-    // TODO: is this still needed?
-    if (mMaxEvents && (mProcessedTFs >= mMaxEvents)) {
-      LOGP(info, "Maximm number of time frames reached ({}), no more processing will be done", mMaxEvents);
-      mReadyToQuit = true;
-      if (mForceQuit) {
-        pc.services().get<ControlService>().endOfStream();
-        pc.services().get<ControlService>().readyToQuit(QuitRequest::All);
-      } else {
-        pc.services().get<ControlService>().readyToQuit(QuitRequest::Me);
-      }
-    }
-  }
+      LOGP(info, "processed sector {} with {} digits and {} reconstructed clusters", sector, inDigits.size(), mClusterFinder->getClusters().size());
 
-  void endOfStream(o2::framework::EndOfStreamContext& ec) final
-  {
-    LOGP(info, "endOfStream");
-    if (mActiveSectors) {
-      snapshotClusters(ec.outputs());
+      mClusterFinder->resetClusters();
     }
-    ec.services().get<ControlService>().readyToQuit(QuitRequest::Me);
+
+    ++mProcessedTFs;
+    LOGP(info, "Number of processed time frames: {}", mProcessedTFs);
   }
 
  private:
-  std::vector<std::vector<KrCluster>> mClusters;
   std::unique_ptr<KrBoxClusterFinder> mClusterFinder;
-  int mLastSector{-1};
-  rawreader::RawReaderCRUManager mRawReader;
-  int mLane{0};                ///< lane number of processor
-  std::vector<int> mSectors{}; ///< sectors to process in this instance
-  uint32_t mMaxEvents{0};
   uint32_t mProcessedTFs{0};
-  bool mReadyToQuit{false};
-  bool mCalibDumped{false};
-  bool mForceQuit{false};
-  uint64_t mActiveSectors{0}; ///< bit mask of active sectors
 
   //____________________________________________________________________________
-  void snapshotClusters(DataAllocator& output)
+  void snapshotClusters(DataAllocator& output, const std::vector<o2::tpc::KrCluster>& clusters, int sector)
   {
-    for (const int sector : mSectors) {
-      o2::tpc::TPCSectorHeader header{sector};
-      header.activeSectors = mActiveSectors;
-      // digit for now are transported per sector, not per lane
-      output.snapshot(Output{"TPC", "KRCLUSTERS", static_cast<SubSpecificationType>(sector), Lifetime::Timeframe, header},
-                      mClusters[sector]);
-    }
-    mActiveSectors = 0;
+    o2::tpc::TPCSectorHeader header{sector};
+    header.activeSectors = (0x1 << sector);
+    output.snapshot(Output{gDataOriginTPC, "KRCLUSTERS", static_cast<SubSpecificationType>(sector), Lifetime::Timeframe, header}, clusters);
   }
 };
 
-o2::framework::DataProcessorSpec getKryptonClustererSpec(const std::string inputSpec, int ilane, std::vector<int> const& sectors)
+o2::framework::DataProcessorSpec getKryptonClustererSpec()
 {
   using device = o2::tpc::KrBoxClusterFinderDevice;
 
+  std::vector<InputSpec> inputs{
+    InputSpec{"digits", gDataOriginTPC, "DIGITS", 0, Lifetime::Timeframe},
+  };
+
   std::vector<OutputSpec> outputs;
-  for (auto isector : sectors) {
-    outputs.emplace_back("TPC", "KRCLUSTERS", static_cast<SubSpecificationType>(isector), Lifetime::Timeframe);
-  }
+  outputs.emplace_back(gDataOriginTPC, "KRCLUSTERS", 0, Lifetime::Timeframe);
 
   return DataProcessorSpec{
-    fmt::format("tpc-krypton-clusterer-{}", ilane),
-    select(inputSpec.data()),
+    "tpc-krypton-clusterer",
+    inputs,
     outputs,
-    AlgorithmSpec{adaptFromTask<device>(ilane, sectors)},
-    Options{
-      {"max-events", VariantType::Int, 0, {"maximum number of events to process"}},
-      {"force-quit", VariantType::Bool, false, {"force quit after max-events have been reached"}},
-    } // end Options
-  };  // end DataProcessorSpec
+    AlgorithmSpec{adaptFromTask<device>()},
+    Options{} // end Options
+  };          // end DataProcessorSpec
 }
 } // namespace tpc
 } // namespace o2