TopologyClassifier.cxx

Go to the documentation of this file.

// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
// All rights not expressly granted are reserved.
//
// This software is distributed under the terms of the GNU General Public
// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
//
// 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.
 
 
#include "IOTOFReconstruction/TopologyClassifier.h"
 
// Include for bitset
#include <bitset>
 
ClassImp(o2::iotof::TopologyClassifier);
 
using std::array;
 
namespace o2
{
namespace iotof
{
 
o2::iotof::Segmentation* TopologyClassifier::sSegmentation = nullptr;
 
void TopologyClassifier::getTopology(uint16_t bitmask, uint16_t minRow, uint8_t spanRow, uint16_t minCol, uint8_t spanCol, uint32_t& topology)
{
 
  // 1. Guard against spans exceeding 8-bit representation for
  // row, col span and 16-bit bitmasks
  if (spanRow > MaxRowSpan || spanCol > MaxColSpan || bitmask > MaxBitmask) {
    topology = Topologies::kHuge;
    return;
  }
 
  const uint32_t clsTopoKey = packKey(spanRow, spanCol, bitmask);
 
  // Check if the topology is already cached
  auto it = mTopologyCache.find(clsTopoKey);
  if (it != mTopologyCache.end()) {
    topology = it->second.mTopology;
    it->second.mFrequency++;
    return;
  }
 
  // Classify the new topology and cache the result
  accountTopology(bitmask, minRow, spanRow, minCol, spanCol);
  topology = mTopologyCache[clsTopoKey].mTopology;
  LOG(debug) << "Classified new topology: " << static_cast<int>(topology);
}
 
TopologyInfo TopologyClassifier::getTopologyFeatures(uint32_t key)
{
  auto it = mTopologyCache.find(key);
  if (it != mTopologyCache.end()) {
    return it->second;
  } else {
    LOG(debug) << "No cached features found for key: " << key;
    return TopologyInfo(); // Return default-constructed TopologyInfo if not found
  }
}
 
void TopologyClassifier::accountTopology(uint16_t bitmask, uint16_t minRow, uint8_t spanRow, uint16_t minCol, uint8_t spanCol)
{
  LOG(debug) << "Classifying topology for bitmask: " << std::bitset<16>(bitmask) << ", minRow: "
             << static_cast<int>(minRow) << ", spanRow: " << static_cast<int>(spanRow)
             << ", minCol: " << static_cast<int>(minCol) << ", spanCol: " << static_cast<int>(spanCol);
 
  // New cluster topology features
  TopologyInfo newTopo;
  newTopo.mFrequency = 1;
  newTopo.mPattern = bitmask;
  newTopo.mSizeX = spanRow;
  newTopo.mSizeZ = spanCol;
  float xCOG{0.f}, zCOG{0.f}, mXMean{0.f}, mZMean{0.f}, mXSigma2{0.f}, mZSigma2{0.f};
  computeCOG(bitmask, minRow, spanRow, minCol, spanCol, newTopo);
 
  const int maxRow = minRow + spanRow - 1;
  const int maxCol = minCol + spanCol - 1;
 
  const auto hasDigit = [bitmask, minRow, minCol, spanCol](int row, int col) -> bool {
    const int bitIndex = (row - minRow) * spanCol + (col - minCol);
    return (bitmask & (1U << bitIndex)) != 0;
  };
 
  // Basic shapes
  if (spanRow == 1 && spanCol == 1) {
    newTopo.mTopology = Topologies::kSingleDigit;
    mTopologyCache[packKey(spanRow, spanCol, bitmask)] = newTopo;
    return;
  }
  if (spanCol == 1) {
    newTopo.mTopology = Topologies::kLineOnRow;
    mTopologyCache[packKey(spanRow, spanCol, bitmask)] = newTopo;
    return;
  }
  if (spanRow == 1) {
    newTopo.mTopology = Topologies::kLineOnCol;
    mTopologyCache[packKey(spanRow, spanCol, bitmask)] = newTopo;
    return;
  }
 
  // Calculate total active digits in the cluster mask
  int firedDigits = 0;
  for (int r = minRow; r <= maxRow; ++r) {
    for (int c = minCol; c <= maxCol; ++c) {
      if (hasDigit(r, c)) {
        firedDigits++;
      }
    }
  }
 
  // Square and rectangles: all pixels fired
  if (firedDigits == spanRow * spanCol && spanRow == spanCol) {
    newTopo.mTopology = Topologies::kSquare;
    mTopologyCache[packKey(spanRow, spanCol, bitmask)] = newTopo;
    return;
  }
  if (firedDigits == spanRow * spanCol && spanRow != spanCol) {
    newTopo.mTopology = Topologies::kRectangle;
    mTopologyCache[packKey(spanRow, spanCol, bitmask)] = newTopo;
    return;
  }
 
  // Corner occupancy
  const bool hasBottomLeft = hasDigit(minRow, minCol);
  const bool hasBottomRight = hasDigit(minRow, maxCol);
  const bool hasTopLeft = hasDigit(maxRow, minCol);
  const bool hasTopRight = hasDigit(maxRow, maxCol);
 
  // Diagonal and triangles
  if (spanRow == spanCol) {
 
    // Triangles
    const int nCorners = hasTopLeft + hasTopRight + hasBottomLeft + hasBottomRight;
    if (nCorners == 3) {
      const int missing = !hasTopLeft ? 0 : !hasTopRight   ? 1
                                          : !hasBottomLeft ? 2
                                                           : 3;
 
      switch (missing) {
        case 0:
          newTopo.mTopology = Topologies::kLowerTriangleLeft;
          break;
        case 1:
          newTopo.mTopology = Topologies::kLowerTriangleRight;
          break;
        case 2:
          newTopo.mTopology = Topologies::kUpperTriangleLeft;
          break;
        case 3:
          newTopo.mTopology = Topologies::kUpperTriangleRight;
          break;
      }
      mTopologyCache[packKey(spanRow, spanCol, bitmask)] = newTopo;
      return;
    }
 
    if ((firedDigits == spanRow && hasTopLeft && hasBottomRight && !hasTopRight && !hasBottomLeft) ||
        (firedDigits == spanRow && hasTopRight && hasBottomLeft && !hasTopLeft && !hasBottomRight)) {
      newTopo.mTopology = Topologies::kDiagonal;
      mTopologyCache[packKey(spanRow, spanCol, bitmask)] = newTopo;
      return;
    }
  }
 
  // Snake: 3 x 2
  if (spanRow == 3 && spanCol == 2) {
    const bool hasMiddleMin = hasDigit(minRow, minCol + 1);
    const bool hasMiddleMax = hasDigit(minRow, maxCol + 1);
 
    if (hasMiddleMin && hasMiddleMax) {
      if (!hasTopLeft && !hasBottomRight && hasTopRight && hasBottomLeft) {
        newTopo.mTopology = Topologies::kSnake;
        mTopologyCache[packKey(spanRow, spanCol, bitmask)] = newTopo;
        return;
      }
 
      if (hasTopLeft && hasBottomRight && !hasTopRight && !hasBottomLeft) {
        newTopo.mTopology = Topologies::kSnakeRefl;
        mTopologyCache[packKey(spanRow, spanCol, bitmask)] = newTopo;
        return;
      }
    }
  }
 
  // Snake rotated by 90 degrees: 2 x 3
  if (spanRow == 2 && spanCol == 3) {
    const bool hasMiddleLeft = hasDigit(minRow + 1, minCol);
    const bool hasMiddleRight = hasDigit(maxRow + 1, minCol);
 
    if (hasMiddleLeft && hasMiddleRight) {
      if (!hasTopLeft && !hasBottomRight && hasTopRight && hasBottomLeft) {
        newTopo.mTopology = Topologies::kSnakeRot90;
        mTopologyCache[packKey(spanRow, spanCol, bitmask)] = newTopo;
        return;
      }
 
      if (hasTopLeft && hasBottomRight && !hasTopRight && !hasBottomLeft) {
        newTopo.mTopology = Topologies::kSnakeRot90Refl;
        mTopologyCache[packKey(spanRow, spanCol, bitmask)] = newTopo;
        return;
      }
    }
  }
 
  if (newTopo.mTopology == Topologies::kNTopologies) {
    newTopo.mTopology = Topologies::kOther;
    mTopologyCache[packKey(spanRow, spanCol, bitmask)] = newTopo;
    return;
  }
}
 
void TopologyClassifier::computeCOG(uint16_t bitmask, uint16_t minRow, uint8_t spanRow, uint16_t minCol, uint8_t spanCol, TopologyInfo& topoInfo)
{
  int xOffsetCOG = 0;
  int zOffsetCOG = 0;
  int firedPixels = 0;
 
  // Ensure nBits does not exceed the bitmask capacity (16 bits)
  const int nBits = std::min(static_cast<int>(spanRow * spanCol), 16);
 
  for (int iBit = 0; iBit < nBits; ++iBit) {
    // Check if the pixel bit is set
    if (bitmask & (1U << iBit)) {
      int iRow = iBit / spanCol;
      int iCol = iBit % spanCol;
 
      xOffsetCOG += minRow + iRow;
      zOffsetCOG += minCol + iCol;
      ++firedPixels;
    }
  }
 
  topoInfo.mOffsetXToCOG = static_cast<int>((static_cast<float>(xOffsetCOG) / firedPixels) - static_cast<float>(minRow));
  topoInfo.mOffsetZToCOG = static_cast<int>((static_cast<float>(zOffsetCOG) / firedPixels) - static_cast<float>(minCol));
  topoInfo.mNPixels = firedPixels;
 
  const auto& chipSpecs = ChipSpecificsParam::Instance();
  topoInfo.mXMean = -(static_cast<float>(xOffsetCOG) / firedPixels - minRow) * chipSpecs.PitchRow;
  topoInfo.mZMean = (static_cast<float>(zOffsetCOG) / firedPixels - minCol) * chipSpecs.PitchCol;
  topoInfo.mXSigma2 = chipSpecs.PitchRow * chipSpecs.PitchRow / 12. / topoInfo.mSizeX;
  topoInfo.mZSigma2 = chipSpecs.PitchCol * chipSpecs.PitchCol / 12. / topoInfo.mSizeZ;
 
  LOG(debug) << "Computed topology features";
  LOG(debug) << "COG offsets: (" << topoInfo.mOffsetXToCOG << ", " << topoInfo.mOffsetZToCOG << ")";
  LOG(debug) << "Shifts to mean: (" << topoInfo.mXMean << ", " << topoInfo.mZMean << ")";
  LOG(debug) << "Sigmas: (" << topoInfo.mXSigma2 << ", " << topoInfo.mZSigma2 << ")";
  LOG(debug) << "Fired Pixels: " << firedPixels;
}
 
math_utils::Point3D<float> TopologyClassifier::getClusterCoordinates(const Cluster& cluster)
{
  if (!mGeometry) {
    LOG(fatal) << "Geometry not set in TopologyClassifier, cannot execute getClusterCoordinates!";
    return math_utils::Point3D<float>{0.f, 0.f, 0.f};
  }
  if (!sSegmentation) {
    LOG(fatal) << "Segmentation not set in TopologyClassifier, cannot execute getClusterCoordinates!";
    return math_utils::Point3D<float>{0.f, 0.f, 0.f};
  }
  auto refRow = cluster.getRow();
  auto refCol = cluster.getCol();
  float x{0.f};
  float z{0.f};
  int layer = mGeometry->getIOTOFLayer(cluster.getChipID());
  sSegmentation->detectorToLocal(cluster.getRow(), cluster.getCol(), x, z, layer);
 
  uint32_t topoKey = makeKey(cluster.getRowSpan(), cluster.getColSpan(), cluster.getPattern());
  x += this->getTopologyFeatures(topoKey).mXMean;
  z += this->getTopologyFeatures(topoKey).mZMean;
  math_utils::Point3D<float> locCl{x, 0.f, z};
 
  return locCl;
}
 
void TopologyClassifier::saveCacheToFile(const char* filename)
{
  TFile file(filename, "RECREATE");
  // Write directly using TObject::Write syntax with explicit class name handling
  file.WriteObject(&mTopologyCache, "TF3ClusterTopologies");
  file.Close();
}
 
void TopologyClassifier::print()
{
  for (const auto& entry : mTopologyCache) {
    const uint32_t key = entry.first;
    const TopologyInfo& topoInfo = entry.second;
 
    uint8_t spanRow = (key >> 24) & 0xFF;
    uint8_t spanCol = (key >> 16) & 0xFF;
    uint16_t bitmask = key & 0xFFFF;
 
    LOG(info) << "Key: " << key
              << ", SpanRow: " << static_cast<int>(spanRow)
              << ", SpanCol: " << static_cast<int>(spanCol)
              << ", Bitmask: " << std::bitset<16>(bitmask);
    topoInfo.print();
  }
}
 
} // namespace iotof
} // namespace o2