KrCluster adjusted for leader pad method and no more clusters found over ROC boundaries. (#4976)

* Added new variables to KrCluster for leader pad method and no more clusters found over ROC boundaries.

* KrCluster Finder can now take gain maps into account.

* Changed some variable types back to unsigned int to save memory.

* Changed the way to specify a CalDet file.

* Included all suggestions regarding KrClusterFinder.

* Search for charge maxima only in exisiting pads.

* A few comments are now more precise.

* Using the default constructor.

Author: HauenderPhilip

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

@@ -53,7 +53,7 @@
 /// the arrays are larger than a sector. For the size in row-direction, a constant is used.
 /// For time and pad direction, a number is set (here is room for improvements).
 ///
-/// When data from a new sector is encountered, the method
+/// ToDo: Find an elegant way to split the huge map into four (IROC, OROC1, OROC2 and OROC3) smaller maps. Unfortunately, this seems to interfere with the rest of the code.
 ///
 /// How to use:
 /// Load tpcdigits.root
@@ -72,6 +72,8 @@
 #include "TPCReconstruction/KrCluster.h"
 #include "TPCBase/Mapper.h"
 
+#include "TPCBase/CalDet.h"
+
 #include <tuple>
 #include <vector>
 #include <array>
@@ -93,7 +95,16 @@ class KrBoxClusterFinder
 {
  public:
   /// Constructor:
-  explicit KrBoxClusterFinder(std::vector<o2::tpc::Digit>& eventSector); ///< Creates a 3D Map
+  /// The constructor allocates a three dimensional array (Pad,Row,Time) which is
+  /// later filled with the recorded charges for each digit
+  explicit KrBoxClusterFinder() = default;
+
+  /// 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);
 
   /// After the map is created, we look for local maxima with this function:
   std::vector<std::tuple<int, int, int>> findLocalMaxima();
@@ -104,26 +115,39 @@ class KrBoxClusterFinder
 
  private:
   // These variables can be varied
+  // They were choses such that the box in each readout chamber is approx. the same size
   int mMaxClusterSizeTime = 3; ///< The "radius" of a cluster in time direction
-  int mMaxClusterSizePad = 3;  ///< radius in pad direction
-  int mMaxClusterSizeRow = 2;  ///< radius in row direction
+  int mMaxClusterSizeRow;      ///< The "radius" of a cluster in row direction
+  int mMaxClusterSizePad;      ///< The "radius" of a cluster in pad direction
+
+  int mMaxClusterSizeRowIROC = 3;  ///< The "radius" of a cluster in row direction in IROC
+  int mMaxClusterSizeRowOROC1 = 2; ///< The "radius" of a cluster in row direction in OROC1
+  int mMaxClusterSizeRowOROC2 = 2; ///< The "radius" of a cluster in row direction in OROC2
+  int mMaxClusterSizeRowOROC3 = 1; ///< The "radius" of a cluster in row direction in OROC3
 
-  // Todo: Differentiate between different ROCS:
-  // int mMaxClusterSizeRowIROC = 3;  // radius in row direction
-  // int mMaxClusterSizeRowOROC1 = 2; // radius in row direction
-  // int mMaxClusterSizeRowOROC2 = 2; // radius in row direction
-  // int mMaxClusterSizeRowOROC3 = 1; // radius in row direction
+  int mMaxClusterSizePadIROC = 5;  ///< The "radius" of a cluster in pad direction in IROC
+  int mMaxClusterSizePadOROC1 = 3; ///< The "radius" of a cluster in pad direction in OROC1
+  int mMaxClusterSizePadOROC2 = 3; ///< The "radius" of a cluster in pad direction in OROC2
+  int mMaxClusterSizePadOROC3 = 3; ///< The "radius" of a cluster in pad direction in OROC3
 
   float mQThresholdMax = 10.0;    ///< the Maximum charge in a cluster must exceed this value or it is discarded
   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
 
+  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
 
+  /// 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
+
   /// Need an instance of Mapper to know position of pads
   const Mapper& mMapperInstance = o2::tpc::Mapper::instance();
 
@@ -132,6 +156,9 @@ class KrBoxClusterFinder
   /// Here the map is defined where all digits are temporarily stored
   std::array<std::array<std::array<float, MaxPads>, MaxRows>, MaxTimes> mMapOfAllDigits{};
 
+  /// For each ROC, the maximum cluster size has to be chosen
+  void setMaxClusterSize(int row);
+
   /// To update the temporary cluster, i.e. all digits are added here
   void updateTempCluster(float tempCharge, int tempPad, int tempRow, int tempTime);
   /// After all digits are assigned to the cluster, the mean and sigmas are calculated here

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

@@ -22,33 +22,37 @@ namespace tpc
 
 struct KrCluster {
  public:
-  float totCharge = 0; ///< Total charge of the cluster (ADC counts)
-  float maxCharge = 0; ///< Maximum charge of the cluster (ADC counts)
-  int size = 0;        ///< Size of the cluster (TPCDigits)
-  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
-  int sector = 0;      ///< Sector number
+  unsigned char size = 0;         ///< Size of the cluster (TPCDigits)
+  unsigned char sector = 0;       ///< Sector number
+  unsigned char maxChargePad = 0; ///< Pad with max. charge in cluster (for leader pad method)
+  unsigned char maxChargeRow = 0; ///< Row with max. charge in cluster (for leader pad method)
+  float totCharge = 0;            ///< Total charge of the cluster (ADC counts)
+  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
 
   /// Used to set all Cluster variables to zero.
   void reset()
   {
+    size = 0;
+    sector = 0;
+    maxChargePad = 0;
+    maxChargeRow = 0;
     totCharge = 0;
     maxCharge = 0;
-    size = 0;
     meanPad = 0;
     meanRow = 0;
     meanTime = 0;
     sigmaPad = 0;
     sigmaRow = 0;
     sigmaTime = 0;
-    sector = 0;
   }
 
-  ClassDefNV(KrCluster, 1);
+  ClassDefNV(KrCluster, 2);
 };
 
 } // namespace tpc

Detectors/TPC/reconstruction/macro/findKrBoxCluster.C

@@ -27,7 +27,7 @@
 #include <vector>
 #endif
 
-void findKrBoxCluster(int lastTimeBin = 1000, int run = -1, int time = -1)
+void findKrBoxCluster(int lastTimeBin = 1000, int run = -1, int time = -1, std::string_view gainMapFile = "")
 {
   // Read the digits:
   TFile* file = new TFile("tpcdigits.root");
@@ -51,22 +51,32 @@ void findKrBoxCluster(int lastTimeBin = 1000, int run = -1, int time = -1)
     tree->SetBranchAddress(Form("TPCDigit_%zu", iSec), &digitizedSignal[iSec]);
   }
 
+  // Create KrBoxClusterFinder object, memory is only allocated once
+  auto clFinder = std::make_unique<o2::tpc::KrBoxClusterFinder>();
+  if (gainMapFile.size()) {
+    clFinder->loadGainMapFromFile(gainMapFile);
+  }
+
   // Now everything can get processed
   // Loop over all events
   for (int iEvent = 0; iEvent < nEntries; ++iEvent) {
     std::cout << iEvent + 1 << "/" << nEntries << std::endl;
     tree->GetEntry(iEvent);
     // Each event consists of sectors (atm only two)
+
     for (int i = 0; i < 36; i++) {
       auto sector = digitizedSignal[i];
       if (sector->size() == 0) {
         continue;
       }
-      // Create ClusterFinder Object on Heap since creation on stack fails
-      // Probably due to too much memory consumption
-      auto clFinder = std::make_unique<o2::tpc::KrBoxClusterFinder>(*sector);
+
+      // Fill map and (if specified) correct with existing gain map
+      clFinder->fillAndCorrectMap(*sector, i);
+
+      // Find all local maxima in sector
       std::vector<std::tuple<int, int, int>> localMaxima = clFinder->findLocalMaxima();
-      // Loop over cluster centers
+
+      // Loop over cluster centers = local maxima
       for (const std::tuple<int, int, int>& coords : localMaxima) {
         int padMax = std::get<0>(coords);
         int rowMax = std::get<1>(coords);
@@ -78,6 +88,7 @@ void findKrBoxCluster(int lastTimeBin = 1000, int run = -1, int time = -1)
         // Build total cluster
         o2::tpc::KrCluster tempCluster = clFinder->buildCluster(padMax, rowMax, timeMax);
         tempCluster.sector = i;
+
         clusters.emplace_back(tempCluster);
       }
     }
@@ -90,3 +101,9 @@ void findKrBoxCluster(int lastTimeBin = 1000, int run = -1, int time = -1)
   fOut->Close();
   return;
 }
+
+int main()
+{
+  findKrBoxCluster();
+  return 0;
+}