Workflow for Krypton Box clusterer

Author: wiechula

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

@@ -99,19 +99,39 @@ class KrBoxClusterFinder
   /// later filled with the recorded charges for each digit
   explicit KrBoxClusterFinder() = default;
 
+  explicit KrBoxClusterFinder(int sector) : mSector(sector){};
+
   /// If a gain map exists, the map can be loaded with this function
   /// The function expects a CalDet file with a gain map (gain entry for each pad).
   void loadGainMapFromFile(const std::string_view calDetFileName, const std::string_view gainMapName = "GainMap");
 
   /// Function used in macro to fill the map with all recorded digits
   void fillAndCorrectMap(std::vector<o2::tpc::Digit>& eventSector, const int sector);
 
+  /// Function to fill single digi
+  void fillADCValue(int cru, int rowInSector, int padInRow, int timeBin, float adcValue);
+
   /// After the map is created, we look for local maxima with this function:
-  std::vector<std::tuple<int, int, int>> findLocalMaxima();
+  std::vector<std::tuple<int, int, int>> findLocalMaxima(bool directFilling = false);
 
   /// After finding the local maxima, we can now build clusters around it.
   /// It works according to the explanation at the beginning of the header file.
-  KrCluster buildCluster(int clusterCenterPad, int clusterCenterRow, int clusterCenterTime);
+  KrCluster buildCluster(int clusterCenterPad, int clusterCenterRow, int clusterCenterTime, bool directFilling = false);
+
+  /// reset the ADC map
+  void resetADCMap();
+
+  /// reset cluster vector
+  void resetClusters() { mClusters.clear(); }
+
+  std::vector<KrCluster>& getClusters() { return mClusters; }
+  const std::vector<KrCluster>& getClusters() const { return mClusters; }
+
+  /// set sector of this instance
+  void setSector(int sector) { mSector = sector; }
+
+  /// ger sector of this instance
+  int getSector() const { return mSector; }
 
  private:
   // These variables can be varied
@@ -134,20 +154,23 @@ class KrBoxClusterFinder
   float mQThreshold = 1.0;        ///< every charge which is added to a cluster must exceed this value or it is discarded
   int mMinNumberOfNeighbours = 1; ///< amount of direct neighbours required for a cluster maximum
 
+  int mSector = -1;                 ///< sector being processed in this instance
   std::unique_ptr<CalPad> mGainMap; ///< Gain map object
 
   /// Maximum Map Dimensions
   /// Here is room for improvements
-  static constexpr size_t MaxPads = 138;  ///< Size of the map in pad-direction
-  static constexpr size_t MaxRows = 152;  ///< Size of the map in row-direction
-  static constexpr size_t MaxTimes = 550; ///< Size of the map in time-direction
+  static constexpr size_t MaxPads = 138;    ///< Size of the map in pad-direction
+  static constexpr size_t MaxRows = 152;    ///< Size of the map in row-direction
+  static constexpr size_t MaxTimes = 20000; ///< Size of the map in time-direction
 
   /// Values to define ROC boundaries
   static constexpr size_t MaxRowsIROC = 63;  ///< Amount of rows in IROC
   static constexpr size_t MaxRowsOROC1 = 34; ///< Amount of rows in OROC1
   static constexpr size_t MaxRowsOROC2 = 30; ///< Amount of rows in OROC2
   static constexpr size_t MaxRowsOROC3 = 25; ///< Amount of rows in OROC3
 
+  std::vector<o2::tpc::KrCluster> mClusters;
+
   /// Need an instance of Mapper to know position of pads
   const Mapper& mMapperInstance = o2::tpc::Mapper::instance();
 

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

@@ -30,10 +30,10 @@ struct KrCluster {
   float maxCharge = 0;            ///< Maximum charge of the cluster (ADC counts)
   float meanPad = 0;              ///< Center of gravity (Pad number)
   float meanRow = 0;              ///< Center of gravity (Row number)
-  float meanTime = 0;             ///< Center of gravity (Time)
   float sigmaPad = 0;             ///< RMS of cluster in pad direction
   float sigmaRow = 0;             ///< RMS of cluster in row direction
-  float sigmaTime = 0;            ///< RMS of cluster in time direction
+  double meanTime = 0;            ///< Center of gravity (Time)
+  double sigmaTime = 0;           ///< RMS of cluster in time direction
 
   /// Used to set all Cluster variables to zero.
   void reset()
@@ -52,7 +52,7 @@ struct KrCluster {
     sigmaTime = 0;
   }
 
-  ClassDefNV(KrCluster, 2);
+  ClassDefNV(KrCluster, 3);
 };
 
 } // namespace tpc

Detectors/TPC/reconstruction/src/KrBoxClusterFinder.cxx

@@ -37,19 +37,26 @@ void KrBoxClusterFinder::loadGainMapFromFile(const std::string_view calDetFileNa
   LOGP(info, "Loaded gain map object '{}' from file '{}'", calDetFileName, gainMapName);
 }
 
-// Fill the map with all digits.
-// You can pass a CalDet file to it so that the cluster finder can already correct for gain inhomogeneities
-// The CalDet File should contain the relative Gain of each Pad.
-void KrBoxClusterFinder::fillAndCorrectMap(std::vector<o2::tpc::Digit>& eventSector, const int sector)
+void KrBoxClusterFinder::resetADCMap()
 {
   // Reset whole map:
   for (int iTime = 0; iTime < MaxTimes; iTime++) {
     for (int iRow = 0; iRow < MaxRows; iRow++) {
-      for (int iPad = 0; iPad < MaxPads; iPad++) {
-        mMapOfAllDigits[iTime][iRow][iPad] = 0;
-      }
+      std::fill(mMapOfAllDigits[iTime][iRow].begin(), mMapOfAllDigits[iTime][iRow].end(), 0.);
+      //for (int iPad = 0; iPad < MaxPads; iPad++) {
+      //mMapOfAllDigits[iTime][iRow][iPad] = 0;
+      //}
     }
   }
+}
+
+// Fill the map with all digits.
+// You can pass a CalDet file to it so that the cluster finder can already correct for gain inhomogeneities
+// The CalDet File should contain the relative Gain of each Pad.
+void KrBoxClusterFinder::fillAndCorrectMap(std::vector<o2::tpc::Digit>& eventSector, const int sector)
+{
+  mSector = sector;
+  resetADCMap();
 
   if (eventSector.size() == 0) {
     // prevents a segementation fault if the envent contains no data
@@ -62,34 +69,49 @@ void KrBoxClusterFinder::fillAndCorrectMap(std::vector<o2::tpc::Digit>& eventSec
     return;
   }
 
-  const auto correctionFactorCalDet = mGainMap.get();
-
   // Fill digits map
   for (const auto& digit : eventSector) {
+    const int cru = digit.getCRU();
     const int time = digit.getTimeStamp();
     const int row = digit.getRow();
     const int pad = digit.getPad();
+    const float adcValue = digit.getChargeFloat();
 
-    float correctionFactor = 1.0;
+    fillADCValue(cru, row, pad, time, adcValue);
+  }
+}
 
-    if (correctionFactorCalDet) {
-      int padNum = mMapperInstance.globalPadNumber(PadPos(row, pad));
-      correctionFactor = correctionFactorCalDet->getValue(sector, padNum);
-    }
-    // Every row starts at pad zero. But the number of pads in a row is not a constant.
-    // If we would just fill the map naively, we would put pads next to each other, which are not neighbours on the pad plane.
-    // Hence, we need to correct for this:
-    const int pads = mMapperInstance.getNumberOfPadsInRowSector(row);
-    const int corPad = pad - (pads / 2) + (MaxPads / 2);
-
-    if (correctionFactor == 0) {
-      LOGP(warning, "Encountered correction factor which is zero.");
-      LOGP(warning, "Digit will be set to 0!");
-      mMapOfAllDigits[time][row][corPad] = 0;
-    } else {
-      mMapOfAllDigits[time][row][corPad] = digit.getChargeFloat() / correctionFactor;
-    }
+void KrBoxClusterFinder::fillADCValue(int cru, int rowInSector, int padInRow, int timeBin, float adcValue)
+{
+  if (timeBin >= MaxTimes) {
+    return;
   }
+
+  // Every row starts at pad zero. But the number of pads in a row is not a constant.
+  // If we would just fill the map naively, we would put pads next to each other, which are not neighbours on the pad plane.
+  // Hence, we need to correct for this:
+  mSector = cru / CRU::CRUperSector;
+  const int pads = mMapperInstance.getNumberOfPadsInRowSector(rowInSector);
+  const int corPad = padInRow - (pads / 2) + (MaxPads / 2);
+
+  const auto correctionFactorCalDet = mGainMap.get();
+  if (!correctionFactorCalDet) {
+    mMapOfAllDigits[timeBin][rowInSector][corPad] = adcValue;
+    return;
+  }
+
+  int padNum = mMapperInstance.globalPadNumber(PadPos(rowInSector, padInRow));
+  float correctionFactor = correctionFactorCalDet->getValue(mSector, padNum);
+
+  if (correctionFactor == 0) {
+    LOGP(warning, "Encountered correction factor which is zero.");
+    LOGP(warning, "Digit will be set to 0!");
+    adcValue = 0;
+  } else {
+    adcValue /= correctionFactor;
+  }
+
+  mMapOfAllDigits[timeBin][rowInSector][corPad] = adcValue;
 }
 
 //#################################################
@@ -114,6 +136,7 @@ void KrBoxClusterFinder::updateTempClusterFinal()
   // Since every padrow is shifted such that neighbouring pads are indeed neighbours, we have to shift once back:
   mTempCluster.meanPad = mTempCluster.meanPad + (corPadsMean / 2.0) - (MaxPads / 2.0);
   mTempCluster.maxChargePad = mTempCluster.maxChargePad + (corPadsMaxCharge / 2.0) - (MaxPads / 2.0);
+  mTempCluster.sector = (decltype(mTempCluster.sector))mSector;
 }
 
 // Function to update the temporal cluster.
@@ -145,7 +168,7 @@ void KrBoxClusterFinder::updateTempCluster(float tempCharge, int tempPad, int te
 
 // This function finds and evaluates all clusters in a 3D mMapOfAllDigits generated by the
 // mMapOfAllDigitsCreator function, this function also updates the cluster tree
-std::vector<std::tuple<int, int, int>> KrBoxClusterFinder::findLocalMaxima()
+std::vector<std::tuple<int, int, int>> KrBoxClusterFinder::findLocalMaxima(bool directFilling)
 {
   std::vector<std::tuple<int, int, int>> localMaximaCoords;
   // loop over whole mMapOfAllDigits the find clusters
@@ -222,7 +245,11 @@ std::vector<std::tuple<int, int, int>> KrBoxClusterFinder::findLocalMaxima()
         if (!thisIsMax) {
           continue;
         } else {
-          localMaximaCoords.emplace_back(std::make_tuple(iPad, iRow, iTime));
+          if (directFilling) {
+            buildCluster(iPad, iRow, iTime, directFilling);
+          } else {
+            localMaximaCoords.emplace_back(std::make_tuple(iPad, iRow, iTime));
+          }
         }
       }
     }
@@ -247,9 +274,8 @@ std::vector<std::tuple<int, int, int>> KrBoxClusterFinder::findLocalMaxima()
 // for loop over whole cluster, to determine if a charge should be added
 // conditions are extrapolation of the 5x5 cluster case to arbitrary
 // cluster sizes in 3 dimensions
-KrCluster KrBoxClusterFinder::buildCluster(int clusterCenterPad, int clusterCenterRow, int clusterCenterTime)
+KrCluster KrBoxClusterFinder::buildCluster(int clusterCenterPad, int clusterCenterRow, int clusterCenterTime, bool directFilling)
 {
-  mTempCluster = KrCluster();
   mTempCluster.reset();
 
   setMaxClusterSize(clusterCenterRow);
@@ -346,6 +372,11 @@ KrCluster KrBoxClusterFinder::buildCluster(int clusterCenterPad, int clusterCent
   // At the end, out mTempCluster should contain all digits that were assigned to the cluster.
   // So before returning it, we update it one last time to calculate the correct means and sigmas.
   updateTempClusterFinal();
+
+  if (directFilling) {
+    mClusters.emplace_back(mTempCluster);
+  }
+
   return mTempCluster;
 }
 
@@ -365,4 +396,4 @@ void KrBoxClusterFinder::setMaxClusterSize(int row)
     mMaxClusterSizePad = mMaxClusterSizePadOROC3;
     mMaxClusterSizeRow = mMaxClusterSizeRowOROC3;
   }
-}
+}

Detectors/TPC/workflow/CMakeLists.txt

@@ -27,6 +27,7 @@ o2_add_library(TPCWorkflow
                        src/ClusterSharingMapSpec.cxx
                        src/FileReaderWorkflow.cxx
                        src/CalDetMergerPublisherSpec.cxx
+                       src/KryptonClustererSpec.cxx
                TARGETVARNAME targetName
                PUBLIC_LINK_LIBRARIES O2::Framework O2::DataFormatsTPC
                                      O2::DPLUtils O2::TPCReconstruction
@@ -58,6 +59,11 @@ o2_add_executable(file-reader
                   SOURCES src/FileReaderWorkflow.cxx
                   PUBLIC_LINK_LIBRARIES O2::TPCWorkflow)
 
+o2_add_executable(krypton-clusterer
+                  COMPONENT_NAME tpc
+                  SOURCES src/tpc-krypton-clusterer.cxx
+                  PUBLIC_LINK_LIBRARIES O2::TPCWorkflow)
+
 o2_add_test(workflow
             COMPONENT_NAME tpc
             LABELS tpc workflow

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

@@ -0,0 +1,34 @@
+// Copyright CERN and copyright holders of ALICE O2. This software is
+// distributed under the terms of the GNU General Public License v3 (GPL
+// Version 3), copied verbatim in the file "COPYING".
+//
+// See http://alice-o2.web.cern.ch/license for full licensing information.
+//
+// In applying this license CERN does not waive the privileges and immunities
+// granted to it by virtue of its status as an Intergovernmental Organization
+// or submit itself to any jurisdiction.
+
+/// @file   RawToDigitsSpec.h
+/// @author Jens Wiechula
+/// @since  2021-02-26
+/// @brief  Processor spec for running TPC Krypton cluster finder
+
+#ifndef TPC_KryptonClustererSpec_H_
+#define TPC_KryptonClustererSpec_H_
+
+#include "Framework/DataProcessorSpec.h"
+#include <string_view>
+
+namespace o2
+{
+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);
+
+} // end namespace tpc
+} // end namespace o2
+
+#endif // TPC_RAWTODIGITSSPEC_H_

Detectors/TPC/workflow/src/CalibProcessingHelper.cxx

@@ -42,6 +42,7 @@ uint64_t calib_processing_helper::processRawData(o2::framework::InputRecord& inp
 
   for (auto const& ref : InputRecordWalker(inputs, filter)) {
     const auto* dh = DataRefUtils::getHeader<o2::header::DataHeader*>(ref);
+    firstOrbit = dh->firstTForbit;
 
     // ---| extract hardware information to do the processing |---
     const auto subSpecification = dh->subSpecification;
@@ -95,7 +96,7 @@ uint64_t calib_processing_helper::processRawData(o2::framework::InputRecord& inp
         }
       }
 
-      firstOrbit = RDHUtils::getHeartBeatOrbit(*rdhPtr);
+      //firstOrbit = RDHUtils::getHeartBeatOrbit(*rdhPtr);
       LOGP(info, "First orbit in present TF: {}", firstOrbit);
       readFirst = true;
     }