Clusterer.h

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.
 
#ifndef ALICEO2_ITS_CLUSTERER_H
#define ALICEO2_ITS_CLUSTERER_H
 
#define _PERFORM_TIMING_
 
// uncomment this to not allow diagonal clusters, e.g. like |* |
//                                                          | *|
#define _ALLOW_DIAGONAL_ALPIDE_CLUSTERS_
 
#include <utility>
#include <vector>
#include <cstring>
#include <memory>
#include <gsl/span>
 
#include "ITSMFTBase/SegmentationAlpide.h"
#include "DataFormatsITSMFT/CompCluster.h"
#include "DataFormatsITSMFT/ROFRecord.h"
#include "ITSMFTReconstruction/PixelReader.h"
#include "ITSMFTReconstruction/PixelData.h"
#include "ITS3Reconstruction/LookUp.h"
#include "SimulationDataFormat/MCCompLabel.h"
#include "CommonConstants/LHCConstants.h"
 
#ifdef _PERFORM_TIMING_
#include <TStopwatch.h>
#endif
 
class TTree;
 
namespace o2
{
class MCCompLabel;
namespace dataformats
{
template <typename T>
class ConstMCTruthContainerView;
template <typename T>
class MCTruthContainer;
} // namespace dataformats
 
namespace its3
{
 
using CompClusCont = std::vector<itsmft::CompClusterExt>;
using PatternCont = std::vector<unsigned char>;
using ROFRecCont = std::vector<itsmft::ROFRecord>;
 
// template <class CompClusCont, class PatternCont, class ROFRecCont> // container types (PMR or std::vectors)
 
class Clusterer
{
  using PixelReader = o2::itsmft::PixelReader;
  using PixelData = o2::itsmft::PixelData;
  using ChipPixelData = o2::itsmft::ChipPixelData;
  using CompCluster = o2::itsmft::CompCluster;
  using CompClusterExt = o2::itsmft::CompClusterExt;
  using Label = o2::MCCompLabel;
  using MCTruth = o2::dataformats::MCTruthContainer<o2::MCCompLabel>;
  using ConstMCTruth = o2::dataformats::ConstMCTruthContainerView<o2::MCCompLabel>;
 
 public:
  static constexpr int MaxLabels = 10;
  static constexpr int MaxHugeClusWarn = 5; // max number of warnings for HugeCluster
 
  struct BBox {
    uint16_t chipID = 0xffff;
    uint16_t rowMin = 0xffff;
    uint16_t colMin = 0xffff;
    uint16_t rowMax = 0;
    uint16_t colMax = 0;
    BBox(uint16_t c) : chipID(c) {}
    bool isInside(uint16_t row, uint16_t col) const { return row >= rowMin && row <= rowMax && col >= colMin && col <= colMax; }
    [[nodiscard]] auto rowSpan() const { return rowMax - rowMin + 1; }
    [[nodiscard]] auto colSpan() const { return colMax - colMin + 1; }
    bool isAcceptableSize() const { return colMax - colMin < o2::itsmft::ClusterPattern::MaxColSpan && rowMax - rowMin < o2::itsmft::ClusterPattern::MaxRowSpan; }
    void clear()
    {
      rowMin = colMin = 0xffff;
      rowMax = colMax = 0;
    }
    void adjust(uint16_t row, uint16_t col)
    {
      if (row < rowMin) {
        rowMin = row;
      }
      if (row > rowMax) {
        rowMax = row;
      }
      if (col < colMin) {
        colMin = col;
      }
      if (col > colMax) {
        colMax = col;
      }
    }
  };
 
  //=========================================================
  struct ThreadStat {
    uint16_t firstChip = 0;
    uint16_t nChips = 0;
    uint32_t firstClus = 0;
    uint32_t firstPatt = 0;
    uint32_t nClus = 0;
    uint32_t nPatt = 0;
  };
 
  struct ClustererThread {
    Clusterer* parent = nullptr; // parent clusterer
    int id = -1;
    // buffers for entries in preClusterIndices in 2 columns, to avoid boundary checks, we reserve
    // extra elements in the beginning and the end
    int* column1 = nullptr;
    int* column2 = nullptr;
    int* curr = nullptr; // pointer on the 1st row of currently processed columnsX
    int* prev = nullptr; // pointer on the 1st row of previously processed columnsX
    int size = itsmft::SegmentationAlpide::NRows + 2;
    // pixels[].first is the index of the next pixel of the same precluster in the pixels
    // pixels[].second is the index of the referred pixel in the ChipPixelData (element of mChips)
    std::vector<std::pair<int, uint32_t>> pixels;
    std::vector<int> preClusterHeads; // index of precluster head in the pixels
    std::vector<int> preClusterIndices;
    uint16_t currCol = 0xffff;               
    bool noLeftCol = true;                   
    std::array<Label, MaxLabels> labelsBuff; 
    std::vector<PixelData> pixArrBuff;       
    //
    CompClusCont compClusters;
    PatternCont patterns;
    MCTruth labels;
    std::vector<ThreadStat> stats; // statistics for each thread results, used at merging
    void resetColumn(int* buff) const { std::memset(buff, -1, sizeof(int) * (size - 2)); }
 
    void swapColumnBuffers() { std::swap(prev, curr); }
 
    void expandPreCluster(uint32_t ip, uint16_t row, int preClusIndex)
    {
      auto& firstIndex = preClusterHeads[preClusterIndices[preClusIndex]];
      pixels.emplace_back(firstIndex, ip);
      firstIndex = pixels.size() - 1;
      curr[row] = preClusIndex;
    }
 
    void addNewPrecluster(uint32_t ip, uint16_t row)
    {
      preClusterHeads.push_back(pixels.size());
      // new head does not point yet (-1) on other pixels, store just the entry of the pixel in the ChipPixelData
      pixels.emplace_back(-1, ip);
      int lastIndex = preClusterIndices.size();
      preClusterIndices.push_back(lastIndex);
      curr[row] = lastIndex; // store index of the new precluster in the current column buffer
    }
 
    void fetchMCLabels(int digID, const ConstMCTruth* labelsDig, int& nfilled);
    void initChip(const ChipPixelData* curChipData, uint32_t first);
    void updateChip(const ChipPixelData* curChipData, uint32_t ip);
    void finishChip(ChipPixelData* curChipData, CompClusCont* compClus, PatternCont* patterns,
                    const ConstMCTruth* labelsDig, MCTruth* labelsClus);
    void finishChipSingleHitFast(uint32_t hit, ChipPixelData* curChipData, CompClusCont* compClusPtr,
                                 PatternCont* patternsPtr, const ConstMCTruth* labelsDigPtr, MCTruth* labelsClusPTr);
    void process(uint16_t chip, uint16_t nChips, CompClusCont* compClusPtr, PatternCont* patternsPtr,
                 const ConstMCTruth* labelsDigPtr, MCTruth* labelsClPtr, const itsmft::ROFRecord& rofPtr);
 
    ~ClustererThread()
    {
      delete[] column1;
      delete[] column2;
    }
 
    ClustererThread(Clusterer* par = nullptr, int _id = -1) : parent(par), id(_id), curr(column2 + 1), prev(column1 + 1) {}
    ClustererThread(const ClustererThread&) = delete;
    ClustererThread(ClustererThread&&) = delete;
    ClustererThread& operator=(const ClustererThread&) = delete;
    ClustererThread& operator=(ClustererThread&&) = delete;
  };
  //=========================================================
 
  Clusterer();
  Clusterer(Clusterer&&) = delete;
  Clusterer& operator=(Clusterer&&) = delete;
  ~Clusterer() = default;
  Clusterer(const Clusterer&) = delete;
  Clusterer& operator=(const Clusterer&) = delete;
 
  void process(int nThreads, PixelReader& r, CompClusCont* compClus, PatternCont* patterns, ROFRecCont* vecROFRec, MCTruth* labelsCl = nullptr);
 
  template <typename VCLUS, typename VPAT>
  static void streamCluster(const std::vector<PixelData>& pixbuf, const std::array<Label, MaxLabels>* lblBuff, const BBox& bbox, const its3::LookUp& pattIdConverter,
                            VCLUS* compClusPtr, VPAT* patternsPtr, MCTruth* labelsClusPtr, int nlab, bool isIB, bool isHuge = false);
 
  bool isContinuousReadOut() const { return mContinuousReadout; }
  void setContinuousReadOut(bool v) { mContinuousReadout = v; }
 
  int getMaxBCSeparationToMask() const { return mMaxBCSeparationToMask; }
  void setMaxBCSeparationToMask(int n) { mMaxBCSeparationToMask = n; }
 
  int getMaxRowColDiffToMask() const { return mMaxRowColDiffToMask; }
  void setMaxRowColDiffToMask(int v) { mMaxRowColDiffToMask = v; }
 
  int getMaxROFDepthToSquash() const { return mSquashingDepth; }
  void setMaxROFDepthToSquash(int v) { mSquashingDepth = v; }
 
  int getMaxBCSeparationToSquash() const { return mMaxBCSeparationToSquash; }
  void setMaxBCSeparationToSquash(int n) { mMaxBCSeparationToSquash = n; }
 
  void print() const;
  void clear();
 
  void setDictionary(const its3::TopologyDictionary* dict)
  {
    LOGP(info, "Setting TopologyDictionary with IB size={} & OB size={}", dict->getSize(true), dict->getSize(false));
    mPattIdConverter.setDictionary(dict);
    // dict->print();
  }
 
  TStopwatch& getTimer() { return mTimer; }           // cannot be const
  TStopwatch& getTimerMerge() { return mTimerMerge; } // cannot be const
 
  void setNChips(int n)
  {
    mChips.resize(n);
    mChipsOld.resize(n);
  }
 
 private:
  void flushClusters(CompClusCont* compClus, MCTruth* labels);
 
  // clusterization options
  bool mContinuousReadout = true; 
 
  int mMaxBCSeparationToMask = static_cast<int>(6000. / o2::constants::lhc::LHCBunchSpacingNS + 10);
  int mMaxRowColDiffToMask = 0; 
  int mNHugeClus = 0;           
 
  int mSquashingDepth = 0; 
  int mMaxBCSeparationToSquash = 6000. / o2::constants::lhc::LHCBunchSpacingNS + 10;
 
  std::vector<std::unique_ptr<ClustererThread>> mThreads; // buffers for threads
  std::vector<ChipPixelData> mChips;                      // currently processed ROF's chips data
  std::vector<ChipPixelData> mChipsOld;                   // previously processed ROF's chips data (for masking)
  std::vector<ChipPixelData*> mFiredChipsPtr;             // pointers on the fired chips data in the decoder cache
 
  its3::LookUp mPattIdConverter; 
                                 //
#ifdef _PERFORM_TIMING_
  TStopwatch mTimer;
  TStopwatch mTimerMerge;
#endif
};
 
template <typename VCLUS, typename VPAT>
void Clusterer::streamCluster(const std::vector<PixelData>& pixbuf, const std::array<Label, MaxLabels>* lblBuff, const Clusterer::BBox& bbox, const its3::LookUp& pattIdConverter,
                              VCLUS* compClusPtr, VPAT* patternsPtr, MCTruth* labelsClusPtr, int nlab, bool isIB, bool isHuge)
{
  if (labelsClusPtr && lblBuff) { // MC labels were requested
    auto cnt = compClusPtr->size();
    for (int i = nlab; i--;) {
      labelsClusPtr->addElement(cnt, (*lblBuff)[i]);
    }
  }
  auto colSpanW = bbox.colSpan();
  auto rowSpanW = bbox.rowSpan();
  // add to compact clusters, which must be always filled
  std::array<unsigned char, itsmft::ClusterPattern::MaxPatternBytes> patt{};
  for (const auto& pix : pixbuf) {
    uint32_t ir = pix.getRowDirect() - bbox.rowMin, ic = pix.getCol() - bbox.colMin;
    int nbits = ir * colSpanW + ic;
    patt[nbits >> 3] |= (0x1 << (7 - (nbits % 8)));
  }
  uint16_t pattID = (isHuge || pattIdConverter.size(isIB) == 0) ? CompCluster::InvalidPatternID : pattIdConverter.findGroupID(rowSpanW, colSpanW, isIB, patt.data());
  uint16_t row = bbox.rowMin, col = bbox.colMin;
  LOGP(debug, "PattID: findGroupID({},{},{})={}", row, col, patt[0], pattID);
  if (pattID == CompCluster::InvalidPatternID || pattIdConverter.isGroup(pattID, isIB)) {
    if (pattID != CompCluster::InvalidPatternID) {
      // For groupped topologies, the reference pixel is the COG pixel
      float xCOG = 0., zCOG = 0.;
      itsmft::ClusterPattern::getCOG(rowSpanW, colSpanW, patt.data(), xCOG, zCOG);
      row += round(xCOG);
      col += round(zCOG);
    }
    if (patternsPtr) {
      patternsPtr->emplace_back((unsigned char)rowSpanW);
      patternsPtr->emplace_back((unsigned char)colSpanW);
      int nBytes = rowSpanW * colSpanW / 8;
      if (((rowSpanW * colSpanW) % 8) != 0) {
        nBytes++;
      }
      patternsPtr->insert(patternsPtr->end(), std::begin(patt), std::begin(patt) + nBytes);
    }
  }
  compClusPtr->emplace_back(row, col, pattID, bbox.chipID);
}
 
} // namespace its3
} // namespace o2
#endif /* ALICEO2_ITS_CLUSTERER_H */