|
| 1 | +// Copyright 2019-2026 CERN and copyright holders of ALICE O2. |
| 2 | +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. |
| 3 | +// All rights not expressly granted are reserved. |
| 4 | +// |
| 5 | +// This software is distributed under the terms of the GNU General Public |
| 6 | +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". |
| 7 | +// |
| 8 | +// In applying this license CERN does not waive the privileges and immunities |
| 9 | +// granted to it by virtue of its status as an Intergovernmental Organization |
| 10 | +// or submit itself to any jurisdiction. |
| 11 | + |
| 12 | +/// \file ClustererACTS.cxx |
| 13 | +/// \brief Implementation of the TRK cluster finder with the ACTS |
| 14 | + |
| 15 | +#include "TRKReconstruction/ClustererACTS.h" |
| 16 | +#include "TRKBase/GeometryTGeo.h" |
| 17 | +#include <Acts/Clusterization/Clusterization.hpp> |
| 18 | + |
| 19 | +#include <algorithm> |
| 20 | +#include <numeric> |
| 21 | + |
| 22 | +using namespace o2::trk; |
| 23 | + |
| 24 | +struct Cell2D { |
| 25 | + Cell2D(int rowv, int colv) : row(rowv), col(colv) {} |
| 26 | + int row, col; |
| 27 | + Acts::Ccl::Label label{Acts::Ccl::NO_LABEL}; |
| 28 | +}; |
| 29 | + |
| 30 | +int getCellRow(const Cell2D& cell) |
| 31 | +{ |
| 32 | + return cell.row; |
| 33 | +} |
| 34 | + |
| 35 | +int getCellColumn(const Cell2D& cell) |
| 36 | +{ |
| 37 | + return cell.col; |
| 38 | +} |
| 39 | + |
| 40 | +bool operator==(const Cell2D& left, const Cell2D& right) |
| 41 | +{ |
| 42 | + return left.row == right.row && left.col == right.col; |
| 43 | +} |
| 44 | + |
| 45 | +bool cellComp(const Cell2D& left, const Cell2D& right) |
| 46 | +{ |
| 47 | + return (left.row == right.row) ? left.col < right.col : left.row < right.row; |
| 48 | +} |
| 49 | + |
| 50 | +struct Cluster2D { |
| 51 | + std::vector<Cell2D> cells; |
| 52 | + std::size_t hash{0}; |
| 53 | +}; |
| 54 | + |
| 55 | +void clusterAddCell(Cluster2D& cl, const Cell2D& cell) |
| 56 | +{ |
| 57 | + cl.cells.push_back(cell); |
| 58 | +} |
| 59 | + |
| 60 | +void hash(Cluster2D& cl) |
| 61 | +{ |
| 62 | + std::ranges::sort(cl.cells, cellComp); |
| 63 | + cl.hash = 0; |
| 64 | + // for (const Cell2D& c : cl.cells) { |
| 65 | + // boost::hash_combine(cl.hash, c.col); |
| 66 | + // } |
| 67 | +} |
| 68 | + |
| 69 | +bool clHashComp(const Cluster2D& left, const Cluster2D& right) |
| 70 | +{ |
| 71 | + return left.hash < right.hash; |
| 72 | +} |
| 73 | + |
| 74 | +template <typename RNG> |
| 75 | +void genclusterw(int x, int y, int x0, int y0, int x1, int y1, |
| 76 | + std::vector<Cell2D>& cells, RNG& rng, double startp = 0.5, |
| 77 | + double decayp = 0.9) |
| 78 | +{ |
| 79 | + std::vector<Cell2D> add; |
| 80 | + |
| 81 | + auto maybe_add = [&](int x_, int y_) { |
| 82 | + Cell2D c(x_, y_); |
| 83 | + // if (std::uniform_real_distribution<double>()(rng) < startp && |
| 84 | + // !rangeContainsValue(cells, c)) { |
| 85 | + // cells.push_back(c); |
| 86 | + // add.push_back(c); |
| 87 | + // } |
| 88 | + }; |
| 89 | + |
| 90 | + // NORTH |
| 91 | + if (y < y1) { |
| 92 | + maybe_add(x, y + 1); |
| 93 | + } |
| 94 | + // NORTHEAST |
| 95 | + if (x < x1 && y < y1) { |
| 96 | + maybe_add(x + 1, y + 1); |
| 97 | + } |
| 98 | + // EAST |
| 99 | + if (x < x1) { |
| 100 | + maybe_add(x + 1, y); |
| 101 | + } |
| 102 | + // SOUTHEAST |
| 103 | + if (x < x1 && y > y0) { |
| 104 | + maybe_add(x + 1, y - 1); |
| 105 | + } |
| 106 | + // SOUTH |
| 107 | + if (y > y0) { |
| 108 | + maybe_add(x, y - 1); |
| 109 | + } |
| 110 | + // SOUTHWEST |
| 111 | + if (x > x0 && y > y0) { |
| 112 | + maybe_add(x - 1, y - 1); |
| 113 | + } |
| 114 | + // WEST |
| 115 | + if (x > x0) { |
| 116 | + maybe_add(x - 1, y); |
| 117 | + } |
| 118 | + // NORTHWEST |
| 119 | + if (x > x0 && y < y1) { |
| 120 | + maybe_add(x - 1, y + 1); |
| 121 | + } |
| 122 | + |
| 123 | + for (Cell2D& c : add) { |
| 124 | + genclusterw(c.row, c.col, x0, y0, x1, y1, cells, rng, startp * decayp, |
| 125 | + decayp); |
| 126 | + } |
| 127 | +} |
| 128 | + |
| 129 | +template <typename RNG> |
| 130 | +Cluster2D gencluster(int x0, int y0, int x1, int y1, RNG& rng, |
| 131 | + double startp = 0.5, double decayp = 0.9) |
| 132 | +{ |
| 133 | + int x0_ = x0 + 1; |
| 134 | + int x1_ = x1 - 1; |
| 135 | + int y0_ = y0 + 1; |
| 136 | + int y1_ = y1 - 1; |
| 137 | + |
| 138 | + int x = std::uniform_int_distribution<std::int32_t>(x0_, x1_)(rng); |
| 139 | + int y = std::uniform_int_distribution<std::int32_t>(y0_, y1_)(rng); |
| 140 | + |
| 141 | + std::vector<Cell2D> cells = {Cell2D(x, y)}; |
| 142 | + genclusterw(x, y, x0_, y0_, x1_, y1_, cells, rng, startp, decayp); |
| 143 | + |
| 144 | + Cluster2D cl; |
| 145 | + cl.cells = std::move(cells); |
| 146 | + |
| 147 | + return cl; |
| 148 | +} |
| 149 | + |
| 150 | +//__________________________________________________ |
| 151 | +void ClustererACTS::process(gsl::span<const Digit> digits, |
| 152 | + gsl::span<const DigROFRecord> digitROFs, |
| 153 | + std::vector<o2::trk::Cluster>& clusters, |
| 154 | + std::vector<unsigned char>& patterns, |
| 155 | + std::vector<o2::trk::ROFRecord>& clusterROFs, |
| 156 | + const ConstDigitTruth* digitLabels, |
| 157 | + ClusterTruth* clusterLabels, |
| 158 | + gsl::span<const DigMC2ROFRecord> digMC2ROFs, |
| 159 | + std::vector<o2::trk::MC2ROFRecord>* clusterMC2ROFs) |
| 160 | +{ |
| 161 | + if (!mThread) { |
| 162 | + mThread = std::make_unique<ClustererThread>(this); |
| 163 | + } |
| 164 | + |
| 165 | + auto* geom = o2::trk::GeometryTGeo::Instance(); |
| 166 | + |
| 167 | + for (size_t iROF = 0; iROF < digitROFs.size(); ++iROF) { |
| 168 | + const auto& inROF = digitROFs[iROF]; |
| 169 | + const auto outFirst = static_cast<int>(clusters.size()); |
| 170 | + const int first = inROF.getFirstEntry(); |
| 171 | + const int nEntries = inROF.getNEntries(); |
| 172 | + |
| 173 | + if (nEntries == 0) { |
| 174 | + clusterROFs.emplace_back(inROF.getBCData(), inROF.getROFrame(), outFirst, 0); |
| 175 | + continue; |
| 176 | + } |
| 177 | + |
| 178 | + // template <typename CellCollection, typename ClusterCollection, |
| 179 | + // std::size_t GridDim = 2, |
| 180 | + // typename Connect = DefaultConnect<typename CellCollection::value_type, GridDim>> |
| 181 | + // requires(GridDim == 1 || GridDim == 2) |
| 182 | + // void createClusters(Acts::Ccl::ClusteringData& data, |
| 183 | + // CellCollection& cells, |
| 184 | + // ClusterCollection& clusters, |
| 185 | + // Connect&& connect = Connect()); |
| 186 | + using Cell = Cell2D; |
| 187 | + using CellCollection = std::vector<Cell>; |
| 188 | + using Cluster = Cluster2D; |
| 189 | + using ClusterCollection = std::vector<Cluster>; |
| 190 | + CellCollection cells; |
| 191 | + Acts::Ccl::ClusteringData data; |
| 192 | + ClusterCollection collection; |
| 193 | + |
| 194 | + Acts::Ccl::createClusters<CellCollection, ClusterCollection, 2>(data, |
| 195 | + cells, |
| 196 | + collection, |
| 197 | + Acts::Ccl::DefaultConnect<Cell, 2>(false)); |
| 198 | + |
| 199 | + // Sort digit indices within this ROF by (chipID, col, row) so we can process |
| 200 | + // chip by chip, column by column -- the same ordering the ALPIDE scanner expects. |
| 201 | + mSortIdx.resize(nEntries); |
| 202 | + std::iota(mSortIdx.begin(), mSortIdx.end(), first); |
| 203 | + std::sort(mSortIdx.begin(), mSortIdx.end(), [&digits](int a, int b) { |
| 204 | + const auto& da = digits[a]; |
| 205 | + const auto& db = digits[b]; |
| 206 | + if (da.getChipIndex() != db.getChipIndex()) { |
| 207 | + return da.getChipIndex() < db.getChipIndex(); |
| 208 | + } |
| 209 | + if (da.getColumn() != db.getColumn()) { |
| 210 | + return da.getColumn() < db.getColumn(); |
| 211 | + } |
| 212 | + return da.getRow() < db.getRow(); |
| 213 | + }); |
| 214 | + |
| 215 | + // Process one chip at a time |
| 216 | + int sliceStart = 0; |
| 217 | + while (sliceStart < nEntries) { |
| 218 | + const int chipFirst = sliceStart; |
| 219 | + const uint16_t chipID = digits[mSortIdx[sliceStart]].getChipIndex(); |
| 220 | + while (sliceStart < nEntries && digits[mSortIdx[sliceStart]].getChipIndex() == chipID) { |
| 221 | + ++sliceStart; |
| 222 | + } |
| 223 | + const int chipN = sliceStart - chipFirst; |
| 224 | + |
| 225 | + mThread->processChip(digits, chipFirst, chipN, &clusters, &patterns, digitLabels, clusterLabels, geom); |
| 226 | + } |
| 227 | + |
| 228 | + clusterROFs.emplace_back(inROF.getBCData(), inROF.getROFrame(), |
| 229 | + outFirst, static_cast<int>(clusters.size()) - outFirst); |
| 230 | + } |
| 231 | + |
| 232 | + if (clusterMC2ROFs && !digMC2ROFs.empty()) { |
| 233 | + clusterMC2ROFs->reserve(clusterMC2ROFs->size() + digMC2ROFs.size()); |
| 234 | + for (const auto& in : digMC2ROFs) { |
| 235 | + clusterMC2ROFs->emplace_back(in.eventRecordID, in.rofRecordID, in.minROF, in.maxROF); |
| 236 | + } |
| 237 | + } |
| 238 | +} |
0 commit comments