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 "ITSMFTReconstruction/LookUp.h"
#include "SimulationDataFormat/MCCompLabel.h"
#include "CommonConstants/LHCConstants.h"
#include "Rtypes.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 itsmft
{
 
using CompClusCont = std::vector<CompClusterExt>;
using PatternCont = std::vector<unsigned char>;
using ROFRecCont = std::vector<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; }
    auto rowSpan() const { return rowMax - rowMin + 1; }
    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 {
    struct PreCluster {
      int head = 0; // index of precluster head in the pixels
      int index = 0;
    };
    int id = -1;
    Clusterer* parent = nullptr; // parent clusterer
    // buffers for entries in preClusterIndices in 2 columns, to avoid boundary checks, we reserve
    // extra elements in the beginning and the end
    int column1[SegmentationAlpide::NRows + 2];
    int column2[SegmentationAlpide::NRows + 2];
    int* curr = nullptr; // pointer on the 1st row of currently processed columnsX
    int* prev = nullptr; // pointer on the 1st row of previously processed columnsX
    // 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;
    uint16_t currCol = 0xffff;               
    bool noLeftCol = true;                   
    std::array<Label, MaxLabels> labelsBuff; 
    std::vector<PixelData> pixArrBuff;       
    std::vector<PreCluster> preClusters;     
    //
    CompClusCont compClusters;
    PatternCont patterns;
    MCTruth labels;
    std::vector<ThreadStat> stats; // statistics for each thread results, used at merging
    void resetColumn(int* buff) { std::memset(buff, -1, sizeof(int) * SegmentationAlpide::NRows); }
 
    void swapColumnBuffers() { std::swap(prev, curr); }
 
    void expandPreCluster(uint32_t ip, uint16_t row, int preClusIndex)
    {
      auto& firstIndex = preClusters[preClusters[preClusIndex].index].head;
      pixels.emplace_back(firstIndex, ip);
      firstIndex = pixels.size() - 1;
      curr[row] = preClusIndex;
    }
 
    void addNewPrecluster(uint32_t ip, uint16_t row)
    {
      int lastIndex = preClusters.size();
      preClusters.emplace_back(pixels.size(), lastIndex);
      // 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);
      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 ROFRecord& rofPtr);
 
    ClustererThread(Clusterer* par = nullptr, int _id = -1) : parent(par), id(_id), curr(column2 + 1), prev(column1 + 1)
    {
      std::fill(std::begin(column1), std::end(column1), -1);
      std::fill(std::begin(column2), std::end(column2), -1);
    }
  };
  //=========================================================
 
  Clusterer();
  ~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 LookUp& pattIdConverter,
                            VCLUS* compClusPtr, VPAT* patternsPtr, MCTruth* labelsClusPtr, int nlab, bool isHuge = false);
 
  bool isContinuousReadOut() const { return mContinuousReadout; }
  void setContinuousReadOut(bool v) { mContinuousReadout = v; }
 
  bool isDropHugeClusters() const { return mDropHugeClusters; }
  void setDropHugeClusters(bool v) { mDropHugeClusters = 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(int layer = -1) const { return (layer < 0) ? mSquashingDepth : mSquashingLayerDepth[layer]; }
  void setMaxROFDepthToSquash(int v) { mSquashingDepth = v; }
  void addMaxROFDepthToSquash(int v) { mSquashingLayerDepth.push_back(v); }
 
  int getMaxBCSeparationToSquash(int layer = -1) const { return (layer < 0) ? mMaxBCSeparationToSquash : mMaxBCSeparationToSquashLayer[layer]; }
  void setMaxBCSeparationToSquash(int n) { mMaxBCSeparationToSquash = n; }
  void addMaxBCSeparationToSquash(int n) { mMaxBCSeparationToSquashLayer.push_back(n); }
 
  void print(bool showTiming = true) const;
  void clear();
  void reset();
 
  void setNChips(int n)
  {
    mChips.resize(n);
    mChipsOld.resize(n);
  }
 
  void loadDictionary(const std::string& fileName) { mPattIdConverter.loadDictionary(fileName); }
  void setDictionary(const TopologyDictionary* dict) { mPattIdConverter.setDictionary(dict); }
 
  TStopwatch& getTimer() { return mTimer; }           // cannot be const
  TStopwatch& getTimerMerge() { return mTimerMerge; } // cannot be const
 
 private:
  void flushClusters(CompClusCont* compClus, MCTruth* labels);
 
  // clusterization options
  bool mContinuousReadout = true; 
  bool mDropHugeClusters = false; 
 
  int mMaxBCSeparationToMask = 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<int> mSquashingLayerDepth;
  std::vector<int> mMaxBCSeparationToSquashLayer;
 
  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
 
  LookUp mPattIdConverter; 
 
  TStopwatch mTimer;
  TStopwatch mTimerMerge;
};
 
template <typename VCLUS, typename VPAT>
void Clusterer::streamCluster(const std::vector<PixelData>& pixbuf, const std::array<Label, MaxLabels>* lblBuff, const Clusterer::BBox& bbox, const LookUp& pattIdConverter,
                              VCLUS* compClusPtr, VPAT* patternsPtr, MCTruth* labelsClusPtr, int nlab, 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, 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() == 0) ? CompCluster::InvalidPatternID : pattIdConverter.findGroupID(rowSpanW, colSpanW, patt.data());
  uint16_t row = bbox.rowMin, col = bbox.colMin;
  if (pattID == CompCluster::InvalidPatternID || pattIdConverter.isGroup(pattID)) {
    if (pattID != CompCluster::InvalidPatternID) {
      // For grouped topologies, the reference pixel is the COG pixel
      float xCOG = 0., zCOG = 0.;
      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 itsmft
} // namespace o2
#endif /* ALICEO2_ITS_CLUSTERER_H */