From be71b8510cbd08ced6f87a213864185a2908d4d4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nicol=C3=B2=20Jacazio?= Date: Sun, 27 Sep 2026 20:49:44 +0200 Subject: [PATCH 1/6] RADIUSMAGNET --- Detectors/Upgrades/ALICE3/Passive/src/Magnet.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Detectors/Upgrades/ALICE3/Passive/src/Magnet.cxx b/Detectors/Upgrades/ALICE3/Passive/src/Magnet.cxx index e9199813d6839..e7013d1a4b225 100644 --- a/Detectors/Upgrades/ALICE3/Passive/src/Magnet.cxx +++ b/Detectors/Upgrades/ALICE3/Passive/src/Magnet.cxx @@ -129,7 +129,7 @@ void Alice3Magnet::ConstructGeometry() mRestMaterialThickness = 6.8f; // cm mOuterWrapInnerRadius = 180.f; // cm mOuterWrapThickness = 3.f; // cm - mZLength = 750.f; // cm Length of the magnet (Z direction) + mZLength = 620.f; // cm Length of the magnet (Z direction) break; case o2::passive::MagnetLayout::MagReducedRadius: mInnerWrapInnerRadius = 125.f; // cm From 19b171df20a3e762b7d06249437ef29f12b31811 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nicol=C3=B2=20Jacazio?= Date: Tue, 17 Mar 2026 08:31:06 +0100 Subject: [PATCH 2/6] dd ACTS tracker --- .../reconstruction/CMakeLists.txt | 10 +- .../TimeFrameMixin.h | 9 +- .../ALICE3GlobalReconstruction/TrackerACTS.h | 46 +- .../ALICE3GlobalReconstruction/TrackerACTS2.h | 236 ++++ .../reconstruction/src/TrackerACTS.cxx | 917 ++++++++++--- .../reconstruction/src/TrackerACTS2.cxx | 1149 +++++++++++++++++ .../workflow/CMakeLists.txt | 4 + .../TrackerSpecImpl.h | 54 +- .../workflow/src/TrackerSpec.cxx | 29 +- .../alice3-global-reconstruction-workflow.cxx | 4 - 10 files changed, 2262 insertions(+), 196 deletions(-) create mode 100644 Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/include/ALICE3GlobalReconstruction/TrackerACTS2.h create mode 100644 Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/src/TrackerACTS2.cxx diff --git a/Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/CMakeLists.txt b/Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/CMakeLists.txt index 68afd31835999..0c665d4ab0411 100644 --- a/Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/CMakeLists.txt +++ b/Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/CMakeLists.txt @@ -10,12 +10,18 @@ # or submit itself to any jurisdiction. if(Acts_FOUND) - set(actsTarget Acts::Core) + set(actsSources src/TrackerACTS.cxx) + set(actsTarget + Acts::Core + O2::ACTSInterface + O2::ALICE3ACTS + ROOT::Tree + ROOT::RIO) endif() o2_add_library(ALICE3GlobalReconstruction TARGETVARNAME targetName - SOURCES src/TimeFrame.cxx + SOURCES src/TimeFrame.cxx ${actsSources} PUBLIC_LINK_LIBRARIES O2::ITStracking O2::GPUCommon diff --git a/Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/include/ALICE3GlobalReconstruction/TimeFrameMixin.h b/Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/include/ALICE3GlobalReconstruction/TimeFrameMixin.h index 4e08b460d999e..4b3cee87a8b70 100644 --- a/Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/include/ALICE3GlobalReconstruction/TimeFrameMixin.h +++ b/Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/include/ALICE3GlobalReconstruction/TimeFrameMixin.h @@ -71,11 +71,12 @@ class TimeFrameMixin : public Base void deriveAndInitTiming(const std::array, nLayers>& layerROFs); const o2::InteractionRecord& getTFAnchorIR() const noexcept { return mTFAnchorIR; } + const auto& getClusterSensorIDs() const noexcept { return mClusterSensorIDs; } protected: void initTimingTables(const std::array& timings); void updateHostROFVertexLookupTable(); - + std::array, nLayers> mClusterSensorIDs; bool mTimingTablesInitialised{false}; o2::InteractionRecord mTFAnchorIR{0, 0}; }; @@ -205,6 +206,7 @@ int TimeFrameMixin::loadROFsFromHitTree(TTree* hitsTree, Geometry this->mUnsortedClusters[iLayer].clear(); this->mTrackingFrameInfo[iLayer].clear(); this->mClusterExternalIndices[iLayer].clear(); + mClusterSensorIDs[iLayer].clear(); this->mClusterSize[iLayer].clear(); } @@ -228,6 +230,7 @@ int TimeFrameMixin::loadROFsFromHitTree(TTree* hitsTree, Geometry this->mUnsortedClusters[iLayer].reserve(clusterCountPerLayer[iLayer]); this->mTrackingFrameInfo[iLayer].reserve(clusterCountPerLayer[iLayer]); this->mClusterExternalIndices[iLayer].reserve(clusterCountPerLayer[iLayer]); + mClusterSensorIDs[iLayer].reserve(clusterCountPerLayer[iLayer]); this->mClusterSize[iLayer].reserve(clusterCountPerLayer[iLayer]); } @@ -294,6 +297,7 @@ int TimeFrameMixin::loadROFsFromHitTree(TTree* hitsTree, Geometry std::array{trkXYZ.y(), trkXYZ.z()}, std::array{resolution[layer] * resolution[layer], 0., resolution[layer] * resolution[layer]}); this->addClusterToLayer(layer, gloXYZ.x(), gloXYZ.y(), gloXYZ.z(), this->mUnsortedClusters[layer].size()); + mClusterSensorIDs[layer].push_back(hit.GetDetectorID()); const int layerHitCounter = hitCounterPerLayer[layer]++; this->addClusterExternalIndexToLayer(layer, layerHitCounter); this->mClusterSize[layer].push_back(1); @@ -344,10 +348,12 @@ int TimeFrameMixin::loadROFrameData(const std::arraymUnsortedClusters[iLayer].clear(); this->mTrackingFrameInfo[iLayer].clear(); this->mClusterExternalIndices[iLayer].clear(); + mClusterSensorIDs[iLayer].clear(); this->mClusterSize[iLayer].clear(); this->mUnsortedClusters[iLayer].reserve(layerClusters[iLayer].size()); this->mTrackingFrameInfo[iLayer].reserve(layerClusters[iLayer].size()); this->mClusterExternalIndices[iLayer].reserve(layerClusters[iLayer].size()); + mClusterSensorIDs[iLayer].reserve(layerClusters[iLayer].size()); this->mClusterSize[iLayer].reserve(layerClusters[iLayer].size()); } @@ -436,6 +442,7 @@ int TimeFrameMixin::loadROFrameData(const std::array{trkXYZ.y(), trkXYZ.z()}, std::array{sigmaY2, 0.f, sigmaZ2}); this->addClusterToLayer(iLayer, gloXYZ.x(), gloXYZ.y(), gloXYZ.z(), this->mUnsortedClusters[iLayer].size()); + mClusterSensorIDs[iLayer].push_back(c.chipID); this->addClusterExternalIndexToLayer(iLayer, clusterId); this->mClusterSize[iLayer].push_back(std::clamp(static_cast(c.size), 0u, 255u)); } diff --git a/Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/include/ALICE3GlobalReconstruction/TrackerACTS.h b/Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/include/ALICE3GlobalReconstruction/TrackerACTS.h index ee69b32a23895..54968b992a2a7 100644 --- a/Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/include/ALICE3GlobalReconstruction/TrackerACTS.h +++ b/Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/include/ALICE3GlobalReconstruction/TrackerACTS.h @@ -19,12 +19,24 @@ #ifndef ALICE3_GLOBALRECONSTRUCTION_INCLUDE_TRACKERACTS_H_ #define ALICE3_GLOBALRECONSTRUCTION_INCLUDE_TRACKERACTS_H_ -#include "Acts/Definitions/Units.hpp" -#include "Framework/Logger.h" - -#include "ITStracking/TimeFrame.h" -#include "TH2F.h" - +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "DataFormatsITS/TrackITS.h" +#include +#include + +namespace Acts +{ +class Surface; +class BoundTrackParameters; +} namespace o2::trk { @@ -62,6 +74,8 @@ struct SpacePoint { int layer{-1}; int clusterId{-1}; int rof{-1}; + int sensorId{-1}; + const Acts::Surface* surface{nullptr}; // Derived quantities float r() const { return std::hypot(x, y); } @@ -98,6 +112,10 @@ class TrackerACTS /// Adopt a TimeFrame for processing void adoptTimeFrame(o2::its::TimeFrame& tf); + void adoptSensorIDs(const std::array, nLayers>& sensorIDs) + { + mSensorIDs = &sensorIDs; + } /// Main tracking entry point: convert clusters to tracks void clustersToTracks(); @@ -125,13 +143,14 @@ class TrackerACTS /// Create seeds (triplets) from space points using ACTS SeedFinder void createSeeds(); - /// Estimate track parameters from a seed using ACTS - bool estimateTrackParams(const SeedACTS& seed, o2::its::TrackITSExt& track) const; - /// Run track finding from seeds void findTracks(); - /// Assign MC labels to tracks + /// Convert fitted ACTS states to a TrackITSExt and publish it in the TimeFrame. + bool acceptTrack(const Acts::BoundTrackParameters& inner, + const Acts::BoundTrackParameters& outer, + const SeedACTS& seed, float chi2); + /// Assign one MC label per accepted track, after all ROFs have been processed. void computeTracksMClabels(); /// Helper: time a task @@ -150,12 +169,17 @@ class TrackerACTS // Seeds std::vector mSeeds; + // Avoid accepting seeds which reuse a cluster already assigned to a track. + std::unordered_set mAcceptedClusters; + // Layer-local indices in mUnsortedClusters, in the same order as getTracks(). + // Keep them separate from the external indices stored in TrackITSExt. + std::vector> mAcceptedTrackClusterIndices; // Tracking state float mBz{0.5f}; ///< Magnetic field in Tesla unsigned int mTimeFrameCounter{0}; double mTotalTime{0.}; - + const std::array, nLayers>* mSensorIDs = nullptr; // Tracking states for logging enum State { SpacePointBuilding = 0, diff --git a/Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/include/ALICE3GlobalReconstruction/TrackerACTS2.h b/Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/include/ALICE3GlobalReconstruction/TrackerACTS2.h new file mode 100644 index 0000000000000..f480f403048d2 --- /dev/null +++ b/Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/include/ALICE3GlobalReconstruction/TrackerACTS2.h @@ -0,0 +1,236 @@ +// 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. + +/// +/// \file TrackerACTS.h +/// \brief TRK tracker using ACTS seeding algorithm +/// \author Nicolò Jacazio, Università del Piemonte Orientale (IT) +/// \since 2026-04-01 +/// + +#ifndef ALICE3_GLOBALRECONSTRUCTION_INCLUDE_TRACKERACTS_H_ +#define ALICE3_GLOBALRECONSTRUCTION_INCLUDE_TRACKERACTS_H_ + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "DataFormatsITS/TrackITS.h" +#include +#include + +class TFile; +class TTree; +namespace Acts +{ +class Surface; +class BoundTrackParameters; +} +namespace o2::trk +{ +/// Configuration for the ACTS-based tracker +struct TrackerACTSConfig { + // Seeding parameters + float minPt = 0.4f * Acts::UnitConstants::GeV; ///< Minimum pT for seeds + float maxImpactParameter = 10.f * Acts::UnitConstants::mm; ///< Maximum impact parameter + float cotThetaMax = std::sinh(4.0f); ///< Maximum cot(theta), corresponds to eta ~4 + // Delta R cuts for doublet/triplet formation + float deltaRMinBottom = 5.f * Acts::UnitConstants::mm; ///< Min deltaR for bottom-middle + float deltaRMaxBottom = 200.f * Acts::UnitConstants::mm; ///< Max deltaR for bottom-middle + float deltaRMinTop = 5.f * Acts::UnitConstants::mm; ///< Min deltaR for middle-top + float deltaRMaxTop = 200.f * Acts::UnitConstants::mm; ///< Max deltaR for middle-top + // Z cuts + float zMin = -3000.f * Acts::UnitConstants::mm; + float zMax = 3000.f * Acts::UnitConstants::mm; + + // Collision region + float collisionRegionMin = -150.f * Acts::UnitConstants::mm; + float collisionRegionMax = 150.f * Acts::UnitConstants::mm; + + // Quality cuts + float maxSeedsPerMiddleSP = 2; + float deltaPhiMax = 0.1f; ///< Maximum phi difference for doublets + // Optional debug output. No ROOT file or histogram is created by default. + bool writeDebugTrees = false; ///< Cluster, seed and ROF summary trees + bool dumpTrackingGeometry = false; ///< Write a JSON volume/surface snapshot on first use + int debugMaxTimeFrames = 1; ///< 1 by default, -1 for all TFs + std::string debugFileName = "TrackerACTS_debug.root"; + std::string geometryDumpFileName = "TrackerACTS_geometry.json"; +}; +/// Space point representation for tracking +struct SpacePoint { + float x{0.f}; + float y{0.f}; + float z{0.f}; + int layer{-1}; + int clusterId{-1}; + int rof{-1}; + int sensorId{-1}; + const Acts::Surface* surface{nullptr}; + + // Derived quantities + float r() const { return std::hypot(x, y); } + float radius() const { return r(); } // required by Acts::CylindricalGridElement concept + float phi() const { return std::atan2(y, x); } + // Variance estimates (can be refined based on cluster properties) + float varianceR{0.01f}; // ~100 um resolution squared + float varianceZ{0.01f}; +}; + +/// Seed (triplet of space points) +struct SeedACTS { + const SpacePoint* bottom{nullptr}; + const SpacePoint* middle{nullptr}; + const SpacePoint* top{nullptr}; + float quality{0.f}; +}; +/// TRK Tracker using ACTS algorithms for seeding and track finding +template +class TrackerACTS +{ + public: + TrackerACTS(); + ~TrackerACTS(); + + /// Adopt a TimeFrame for processing + void adoptTimeFrame(o2::its::TimeFrame& tf); + void adoptSensorIDs(const std::array, nLayers>& sensorIDs) + { + mSensorIDs = &sensorIDs; + } + + /// Main tracking entry point: convert clusters to tracks + void clustersToTracks(); + + /// Configuration + void setConfig(const TrackerACTSConfig& cfg) { mConfig = cfg; } + TrackerACTSConfig& getConfig() { return mConfig; } + const TrackerACTSConfig& getConfig() const { return mConfig; } + /// Configure the magnetic field strength (T). + void setBz(float bz) { mBz = bz; } + float getBz() const { return mBz; } + + /// Print tracking summary + void printSummary() const; + + private: + // Debug resources are allocated only if requested and closed at destruction. + std::unique_ptr mDebugFile; + TTree* mClusterTree = nullptr; // Owned by mDebugFile + TTree* mSeedTree = nullptr; // Owned by mDebugFile + TTree* mROFTree = nullptr; // Owned by mDebugFile + std::unique_ptr mHistSpacePoints; + bool mGeometryDumped = false; + + struct DebugClusterRecord { + int tf = 0, rof = 0, layer = 0, clusterId = -1; + int externalId = -1, sensorId = -1, mappingStatus = 0; + unsigned long long geometryId = 0; + float x = 0.f, y = 0.f, z = 0.f; + float r = 0.f, phi = 0.f; + float varianceR = 0.f, varianceZ = 0.f; + } mDebugCluster; + + struct DebugSeedRecord { + int tf = 0, rof = 0, seedId = 0; + float quality = 0.f; + int layer[3] = {}, clusterId[3] = {}, sensorId[3] = {}; + unsigned long long geometryId[3] = {}; + float x[3] = {}, y[3] = {}, z[3] = {}, r[3] = {}, phi[3] = {}; + } mDebugSeed; + + struct DebugROFRecord { + int tf = 0, rof = 0, nSpacePoints = 0, nSeeds = 0; + int nAcceptedTracks = 0; + } mDebugROF; + + void initDebugOutput(); + void dumpGeometryOnce(); + bool isDebugTimeFrame() const; + void fillDebugSeeds(int rof); + + /// Build space points from clusters in the TimeFrame + void buildSpacePoints(int rof); + + /// Create seeds (triplets) from space points using ACTS SeedFinder + void createSeeds(); + + /// Run track finding from seeds + void findTracks(); + + /// Convert fitted ACTS states to a TrackITSExt and publish it in the TimeFrame. + bool acceptTrack(const Acts::BoundTrackParameters& inner, + const Acts::BoundTrackParameters& outer, + const SeedACTS& seed, float chi2); + /// Assign one MC label per accepted track, after all ROFs have been processed. + void computeTracksMClabels(); + + /// Helper: time a task + template + float evaluateTask(Func&& task, std::string_view taskName); + + // Configuration + TrackerACTSConfig mConfig; + + // TimeFrame data + o2::its::TimeFrame* mTimeFrame = nullptr; + + // Space points built from clusters + std::vector mSpacePoints; + std::vector> mSpacePointsPerLayer; + // Seeds + std::vector mSeeds; + // Avoid accepting seeds which reuse a cluster already assigned to a track. + std::unordered_set mAcceptedClusters; + // Layer-local indices in mUnsortedClusters, in the same order as getTracks(). + // Keep them separate from the external indices stored in TrackITSExt. + std::vector> mAcceptedTrackClusterIndices; + + // Tracking state + float mBz{0.5f}; ///< Magnetic field in Tesla + unsigned int mTimeFrameCounter{0}; + double mTotalTime{0.}; + const std::array, nLayers>* mSensorIDs = nullptr; + // Tracking states for logging + enum State { + SpacePointBuilding = 0, + Seeding, + TrackFinding, + NStates, + }; + State mCurState{SpacePointBuilding}; + static constexpr std::array StateNames{ + "Space point building", + "ACTS seeding", + "Track finding"}; +}; +template +template +float TrackerACTS::evaluateTask(Func&& task, std::string_view taskName) +{ + LOG(debug) << " + Starting " << taskName; + const auto start = std::chrono::high_resolution_clock::now(); + task(); + const auto end = std::chrono::high_resolution_clock::now(); + std::chrono::duration diff{end - start}; + LOG(debug) << " - " << taskName << " completed in: " << std::fixed << std::setprecision(2) << diff.count() << " ms"; + return static_cast(diff.count()); +} + +} // namespace o2::trk +#endif /* ALICE3_GLOBALRECONSTRUCTION_INCLUDE_TRACKERACTS_H_ */ diff --git a/Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/src/TrackerACTS.cxx b/Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/src/TrackerACTS.cxx index e870ee934816f..59d0281846433 100644 --- a/Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/src/TrackerACTS.cxx +++ b/Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/src/TrackerACTS.cxx @@ -18,26 +18,169 @@ #include "ALICE3GlobalReconstruction/TrackerACTS.h" -#include +#include "ACTSInterface/TrackingGeometryManager.h" +#include "MathUtils/Utils.h" +#include "TRKBase/GeometryTGeo.h" + +#include +#include +#include +#include +#include +#include #include -#include -#include -#include -#include -#include -#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include #include -#include + +#include "ACTSInterface/TrackingGeometryManager.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include namespace o2::trk { +namespace +{ +// Convert ACTS surface-bound parameters to the O2 tracking frame +// (x, alpha; y, z, sin(phi-alpha), tan(lambda), q/pT). +// The covariance is transformed by a numerical Jacobian evaluated at the +// fitted state. Position units are converted from ACTS mm to O2 cm. +bool convertTrackParameters(const Acts::BoundTrackParameters& actsTrack, + const Acts::GeometryContext& gctx, + o2::track::TrackParCov& output) +{ + const auto& state = actsTrack.parameters(); + if (!state.allFinite() || !actsTrack.covariance().has_value() || + !actsTrack.covariance()->allFinite()) { + return false; + } + + const auto position = actsTrack.position(gctx); + if (!position.allFinite()) { + return false; + } + const double alpha = std::atan2(position.y(), position.x()); + const double ca = std::cos(alpha); + const double sa = std::sin(alpha); + + auto toO2 = [&](const Acts::BoundVector& p, Eigen::Matrix& dest) -> bool { + const double theta = p[Acts::eBoundTheta]; + const double sinTheta = std::sin(theta); + if (!std::isfinite(sinTheta) || sinTheta <= 1.e-6) { + return false; + } + const auto direction = Acts::makeDirectionFromPhiTheta(p[Acts::eBoundPhi], theta); + const Acts::Vector2 local = p.head<2>(); + const auto xyz = actsTrack.referenceSurface().localToGlobal(gctx, local, direction); + if (!xyz.allFinite()) { + return false; + } + + const double snp = std::sin(p[Acts::eBoundPhi] - alpha); + if (!std::isfinite(snp) || std::abs(snp) >= 0.999999) { + return false; + } + + dest << (-xyz.x() * sa + xyz.y() * ca) / Acts::UnitConstants::cm, + xyz.z() / Acts::UnitConstants::cm, + snp, + std::cos(theta) / sinTheta, + p[Acts::eBoundQOverP] * Acts::UnitConstants::GeV / sinTheta; + return dest.allFinite(); + }; + + Eigen::Matrix pars; + if (!toO2(state, pars)) { + return false; + } + + const double x = (position.x() * ca + position.y() * sa) / Acts::UnitConstants::cm; + if (!std::isfinite(x)) { + return false; + } + + Eigen::Matrix jacobian; + jacobian.setZero(); + for (int i = 0; i < Acts::eBoundSize; ++i) { + if (i == Acts::eBoundTime) { + continue; // O2 TrackParCov has no time coordinate + } + const double step = (i == Acts::eBoundLoc0 || i == Acts::eBoundLoc1) + ? 1.e-3 * Acts::UnitConstants::mm + : (i == Acts::eBoundQOverP) + ? std::max(std::abs(state[i]) * 1.e-4, 1.e-6 / Acts::UnitConstants::GeV) + : 1.e-5; + Acts::BoundVector plus = state; + Acts::BoundVector minus = state; + plus[i] += step; + minus[i] -= step; + Eigen::Matrix fPlus, fMinus; + if (!toO2(plus, fPlus) || !toO2(minus, fMinus)) { + return false; + } + jacobian.col(i) = (fPlus - fMinus) / (2. * step); + } + + const Eigen::Matrix covariance = + jacobian * (*actsTrack.covariance()) * jacobian.transpose(); + if (!covariance.allFinite()) { + return false; + } + + o2::track::TrackParCov::params_t o2Params{}; + o2::track::TrackParCov::covMat_t o2Cov{}; + for (int i = 0; i < 5; ++i) { + if (!(covariance(i, i) > 0.) || !std::isfinite(pars[i])) { + return false; + } + o2Params[i] = static_cast(pars[i]); + for (int j = 0; j <= i; ++j) { + o2Cov[o2::track::CovarMap[i][j]] = + static_cast(0.5 * (covariance(i, j) + covariance(j, i))); + } + } + const int charge = state[Acts::eBoundQOverP] >= 0. ? 1 : -1; + output = o2::track::TrackParCov{static_cast(x), static_cast(alpha), + o2Params, o2Cov, charge, o2::track::PID::Pion}; + return output.isValid(); +} +} // namespace template TrackerACTS::TrackerACTS() { // Initialize space points storage per layer mSpacePointsPerLayer.resize(nLayers); mHistSpacePoints = new TH2F("hSpacePoints", "Space points; x (cm); y (cm)", 200, -100, 100, 200, -100, 100); + + // Initialize ACTS Tracking geometry } template @@ -54,19 +197,48 @@ void TrackerACTS::buildSpacePoints(int rof) layerSPs.clear(); } - // Get clusters from the TimeFrame and convert to space points + if (!mSensorIDs) { + LOG(error) << "Cannot build ACTS space points: Sensor IDs not adopted"; + return; + } + + auto& actsGeometryManager = o2::acts::TrackingGeometryManager::instance(); + auto* trkGeometry = o2::trk::GeometryTGeo::Instance(); + + trkGeometry->fillMatrixCache(o2::math_utils::bit2Mask(o2::math_utils::TransformType::L2G)); + + const auto& surfaceIndex = actsGeometryManager.getIndex(*trkGeometry, 1.e-3, [trkGeometry](int sensorID) { + return std::string(trkGeometry->getMatrixPath(sensorID).Data()); + }); + + int nMapped = 0; + int nMissing = 0; + for (int layer = 0; layer < nLayers; ++layer) { - // For now we take unsorted clusters, as soon as the cluster trackin is in place we can piggy back on it and switch to the clusters - auto clusters = mTimeFrame->getUnsortedClusters()[layer]; - // Resize the clusters to the first 100 clusters for testing - // clusters = clusters.subspan(0, std::min(clusters.size(), 100)); - LOG(debug) << "ACTSTracker: got " << clusters.size() << " clusters"; + const auto clusters = mTimeFrame->getUnsortedClustersOnLayer(rof, layer); + const int firstCluster = mTimeFrame->getSortedStartIndex(rof, layer); + + LOG(debug) << "ACTSTracker: got " << clusters.size() << " clusters on layer " << layer; for (size_t iCluster = 0; iCluster < clusters.size(); ++iCluster) { const auto& cluster = clusters[iCluster]; + const int clusterIndex = firstCluster + static_cast(iCluster); + + if (static_cast(clusterIndex) >= (*mSensorIDs)[layer].size()) { + LOG(error) << "Missing sensor ID for layer " << layer << " cluster " << clusterIndex; + ++nMissing; + continue; + } + + const int sensorID = (*mSensorIDs)[layer][clusterIndex]; + const auto* surface = surfaceIndex.getSurface(sensorID); + if (!surface) { + LOG(error) << "No ACTS surface for sensor " << sensorID << " on layer " << layer; + ++nMissing; + continue; + } SpacePoint sp; - // Check that these are in global coordinates sp.x = cluster.xCoordinate * Acts::UnitConstants::cm; sp.y = cluster.yCoordinate * Acts::UnitConstants::cm; sp.z = cluster.zCoordinate * Acts::UnitConstants::cm; @@ -74,181 +246,606 @@ void TrackerACTS::buildSpacePoints(int rof) if (mHistSpacePoints) { mHistSpacePoints->Fill(sp.x / Acts::UnitConstants::cm, sp.y / Acts::UnitConstants::cm); } + sp.layer = layer; - sp.clusterId = static_cast(iCluster); + sp.clusterId = clusterIndex; sp.rof = rof; - - // Position uncertainties (could be refined based on cluster properties) - sp.varianceR = 0.01f; // ~100 um resolution squared + sp.sensorId = sensorID; + sp.surface = surface; + sp.varianceR = 0.01f; sp.varianceZ = 0.01f; mSpacePoints.push_back(sp); + ++nMapped; } } - // Build per-layer pointers for seeding for (auto& sp : mSpacePoints) { if (sp.layer >= 0 && sp.layer < nLayers) { mSpacePointsPerLayer[sp.layer].push_back(&sp); } } + + LOG(info) << "ACTS surface mapping: " << nMapped << " clusters mapped, " << nMissing << " clusters missing"; } template void TrackerACTS::createSeeds() { + mSeeds.clear(); if (mSpacePoints.empty()) { - LOGF(info, "No space points available for seeding"); + LOG(info) << "No space points available for seeding"; return; } - mSeeds.clear(); - // Backend adaptor that exposes mSpacePoints to Acts::SpacePointContainer - struct SpacePointBackend { - using ValueType = SpacePoint; - explicit SpacePointBackend(const std::vector& sps) : m_sps{&sps} {} - std::size_t size_impl() const { return m_sps->size(); } - float x_impl(std::size_t i) const { return (*m_sps)[i].x; } - float y_impl(std::size_t i) const { return (*m_sps)[i].y; } - float z_impl(std::size_t i) const { return (*m_sps)[i].z; } - float varianceR_impl(std::size_t i) const { return (*m_sps)[i].varianceR; } - float varianceZ_impl(std::size_t i) const { return (*m_sps)[i].varianceZ; } - const SpacePoint& get_impl(std::size_t i) const { return (*m_sps)[i]; } - std::any component_impl(Acts::HashedString /*key*/, std::size_t /*i*/) const - { - LOG(fatal) << "No additional components available for space points"; - throw std::runtime_error("SpacePointBackend: no strip component available"); - } - const std::vector* m_sps; - }; + const float bFieldInZ = mBz * Acts::UnitConstants::T; + + Acts::CylindricalSpacePointGrid::Config gridConfig; + gridConfig.minPt = mConfig.minPt; + gridConfig.rMin = 0.f; + gridConfig.rMax = 100.f * Acts::UnitConstants::cm; + gridConfig.zMin = mConfig.zMin; + gridConfig.zMax = mConfig.zMax; + gridConfig.deltaRMax = std::max(mConfig.deltaRMaxBottom, mConfig.deltaRMaxTop); + gridConfig.cotThetaMax = mConfig.cotThetaMax; + gridConfig.impactMax = mConfig.maxImpactParameter; + gridConfig.bFieldInZ = bFieldInZ; + gridConfig.bottomBinFinder = Acts::GridBinFinder<3ul>{1, std::vector>{}, 0}; + gridConfig.topBinFinder = Acts::GridBinFinder<3ul>{1, std::vector>{}, 0}; + + Acts::CylindricalSpacePointGrid grid{gridConfig}; + for (Acts::SpacePointIndex i = 0; i < mSpacePoints.size(); ++i) { + const auto& sp = mSpacePoints[i]; + grid.insert(i, sp.phi(), sp.z, sp.r()); + } - // Wrap mSpacePoints in an Acts space-point container - SpacePointBackend backend{mSpacePoints}; + // The doublet finder expects each bin to be sorted in increasing radius. + for (size_t i = 0; i < grid.numberOfBins(); ++i) { + auto& bin = grid.at(i); + std::sort(bin.begin(), bin.end(), [this](Acts::SpacePointIndex a, Acts::SpacePointIndex b) { + return mSpacePoints[a].r() < mSpacePoints[b].r(); + }); + } - // Configure the ACTS space point container - Acts::SpacePointContainerConfig spContainerConfig; - Acts::SpacePointContainerOptions spContainerOpts; - spContainerOpts.beamPos = {0.f, 0.f}; - Acts::SpacePointContainer spContainer{spContainerConfig, spContainerOpts, backend}; + // ACTS 47 uses a non-templated, column-based space-point container. + Acts::SpacePointContainer spacePoints{ + Acts::SpacePointColumns::CopiedFromIndex | + Acts::SpacePointColumns::PackedXY | + Acts::SpacePointColumns::PackedZR | + Acts::SpacePointColumns::VarianceR | + Acts::SpacePointColumns::VarianceZ}; + spacePoints.reserve(grid.numberOfSpacePoints()); + + std::vector gridRanges; + gridRanges.reserve(grid.numberOfBins()); + for (size_t i = 0; i < grid.numberOfBins(); ++i) { + const auto first = static_cast(spacePoints.size()); + for (Acts::SpacePointIndex spIndex : grid.at(i)) { + const auto& sp = mSpacePoints[spIndex]; + auto proxy = spacePoints.createSpacePoint(); + proxy.copiedFromIndex() = spIndex; + proxy.xy() = std::array{sp.x, sp.y}; + proxy.zr() = std::array{sp.z, sp.r()}; + proxy.varianceR() = sp.varianceR; + proxy.varianceZ() = sp.varianceZ; + } + const auto last = static_cast(spacePoints.size()); + gridRanges.emplace_back(first, last); + } - // Configure the ACTS seed finder - const unsigned int maxSeeds = static_cast(mConfig.maxSeedsPerMiddleSP); - Acts::SeedFilterConfig filterConfig; - filterConfig.maxSeedsPerSpM = maxSeeds; + Acts::DoubletSeedFinder::Config bottomConfig; + bottomConfig.spacePointsSortedByRadius = true; + bottomConfig.candidateDirection = Acts::Direction::Backward(); + bottomConfig.deltaRMin = mConfig.deltaRMinBottom; + bottomConfig.deltaRMax = mConfig.deltaRMaxBottom; + bottomConfig.collisionRegionMin = mConfig.collisionRegionMin; + bottomConfig.collisionRegionMax = mConfig.collisionRegionMax; + bottomConfig.cotThetaMax = mConfig.cotThetaMax; + bottomConfig.impactMax = mConfig.maxImpactParameter; + bottomConfig.minPt = mConfig.minPt; + auto bottomFinder = Acts::DoubletSeedFinder::create( + Acts::DoubletSeedFinder::DerivedConfig{bottomConfig, bFieldInZ}); + + Acts::DoubletSeedFinder::Config topConfig = bottomConfig; + topConfig.candidateDirection = Acts::Direction::Forward(); + topConfig.deltaRMin = mConfig.deltaRMinTop; + topConfig.deltaRMax = mConfig.deltaRMaxTop; + auto topFinder = Acts::DoubletSeedFinder::create( + Acts::DoubletSeedFinder::DerivedConfig{topConfig, bFieldInZ}); + + Acts::TripletSeedFinder::Config tripletConfig; + tripletConfig.minPt = mConfig.minPt; + tripletConfig.impactMax = mConfig.maxImpactParameter; + tripletConfig.sigmaScattering = 5.f; + tripletConfig.radLengthPerSeed = 0.05f; + auto tripletFinder = Acts::TripletSeedFinder::create( + Acts::TripletSeedFinder::DerivedConfig{tripletConfig, bFieldInZ}); + + Acts::BroadTripletSeedFilter::Config filterConfig; + filterConfig.maxSeedsPerSpM = static_cast(mConfig.maxSeedsPerMiddleSP); + filterConfig.deltaRMin = std::min(mConfig.deltaRMinBottom, mConfig.deltaRMinTop); + Acts::BroadTripletSeedFilter::State filterState; + Acts::BroadTripletSeedFilter::Cache filterCache; + Acts::BroadTripletSeedFilter filter{filterConfig, filterState, filterCache, Acts::getDummyLogger()}; + + Acts::TripletSeeder seeder; + Acts::TripletSeeder::Cache seedingCache; + Acts::SeedContainer seeds; + seeds.assignSpacePointContainer(spacePoints); + + std::vector bottomRanges; + std::vector topRanges; + for (const auto [bottom, middle, top] : grid.binnedGroup()) { + bottomRanges.clear(); + for (const auto bin : bottom) { + bottomRanges.push_back(spacePoints.range(gridRanges.at(bin)).asConst()); + } + topRanges.clear(); + for (const auto bin : top) { + topRanges.push_back(spacePoints.range(gridRanges.at(bin)).asConst()); + } + auto middleRange = spacePoints.range(gridRanges.at(middle)).asConst(); + if (middleRange.empty()) { + continue; + } - // ACTS requires minPt / bFieldInZ >= rMax / 2 (minHelixRadius >= rMax/2). - // Cap rMax so that the constraint is satisfied for the configured minPt and field. - const float bFieldInZ = mBz * Acts::UnitConstants::T; - const float safeRMax = 1.8f * mConfig.minPt / bFieldInZ; // 10% margin below the hard limit - - using SPProxy = typename Acts::SpacePointContainer::SpacePointProxyType; - Acts::SeedFinderConfig finderConfig; - finderConfig.rMin = 0.f; - finderConfig.rMax = 100.f * Acts::UnitConstants::cm; - finderConfig.zMin = mConfig.zMin; - finderConfig.zMax = mConfig.zMax; - finderConfig.deltaRMin = std::min(mConfig.deltaRMinBottom, mConfig.deltaRMinTop); - finderConfig.deltaRMax = std::max(mConfig.deltaRMaxBottom, mConfig.deltaRMaxTop); - finderConfig.deltaRMinBottomSP = mConfig.deltaRMinBottom; - finderConfig.deltaRMaxBottomSP = mConfig.deltaRMaxBottom; - finderConfig.deltaRMinTopSP = mConfig.deltaRMinTop; - finderConfig.deltaRMaxTopSP = mConfig.deltaRMaxTop; - finderConfig.collisionRegionMin = mConfig.collisionRegionMin; - finderConfig.collisionRegionMax = mConfig.collisionRegionMax; - finderConfig.cotThetaMax = mConfig.cotThetaMax; - finderConfig.minPt = mConfig.minPt; - finderConfig.impactMax = mConfig.maxImpactParameter; - finderConfig.maxSeedsPerSpM = maxSeeds; - finderConfig.sigmaScattering = 5.f; - finderConfig.radLengthPerSeed = 0.05f; - finderConfig.seedFilter = std::make_shared>(filterConfig); - finderConfig = finderConfig.calculateDerivedQuantities(); - Acts::SeedFinder> seedFinder{finderConfig, - Acts::getDefaultLogger("Finder", Acts::Logging::Level::VERBOSE)}; - - // Configure and create the cylindrical space-point grid - Acts::CylindricalSpacePointGridConfig gridConfig; - gridConfig.minPt = finderConfig.minPt; - gridConfig.rMin = finderConfig.rMin; - gridConfig.rMax = finderConfig.rMax; - gridConfig.zMin = finderConfig.zMin; - gridConfig.zMax = finderConfig.zMax; - gridConfig.deltaRMax = finderConfig.deltaRMax; - gridConfig.cotThetaMax = finderConfig.cotThetaMax; - gridConfig.impactMax = finderConfig.impactMax; - - Acts::CylindricalSpacePointGridOptions gridOpts; - gridOpts.bFieldInZ = bFieldInZ; - - Acts::SeedFinderOptions finderOpts; - finderOpts.beamPos = spContainerOpts.beamPos; - finderOpts.bFieldInZ = gridOpts.bFieldInZ; - try { - finderOpts = finderOpts.calculateDerivedQuantities(finderConfig); - } catch (const std::exception& e) { - LOG(fatal) << "Error in seed finder configuration: " << e.what(); - return; + seeder.createSeedsFromGroups(seedingCache, *bottomFinder, *topFinder, + *tripletFinder, filter, spacePoints, + bottomRanges, middleRange, topRanges, + {gridConfig.rMin, gridConfig.rMax}, seeds); + } + + // Restore the existing SeedACTS representation using original SP indices. + mSeeds.reserve(seeds.size()); + for (const auto seed : seeds) { + const auto indices = seed.spacePointIndices(); + SeedACTS output; + output.bottom = &mSpacePoints[spacePoints.at(indices[0]).copiedFromIndex()]; + output.middle = &mSpacePoints[spacePoints.at(indices[1]).copiedFromIndex()]; + output.top = &mSpacePoints[spacePoints.at(indices[2]).copiedFromIndex()]; + output.quality = seed.quality(); + mSeeds.push_back(output); } - Acts::CylindricalSpacePointGrid grid = Acts::CylindricalSpacePointGridCreator::createGrid(gridConfig, gridOpts); - try { - Acts::CylindricalSpacePointGridCreator::fillGrid(finderConfig, finderOpts, grid, - spContainer.begin(), spContainer.end()); - } catch (const std::exception& e) { - LOG(fatal) << "Error during grid creation/filling: " << e.what(); + LOGF(info, "Created %zu seeds from %zu space points", mSeeds.size(), mSpacePoints.size()); +} + +template +void TrackerACTS::findTracks() +{ + if (mSeeds.empty()) { return; } - LOG(debug) << "Grid created with " << grid.dimensions(); - // Build the binned group and iterate over triplet combinations - Acts::GridBinFinder<3ul> bottomBinFinder{1, std::vector>{}, 0}; - Acts::GridBinFinder<3ul> topBinFinder{1, std::vector>{}, 0}; - Acts::CylindricalBinnedGroup spGroup{std::move(grid), bottomBinFinder, topBinFinder}; + auto& actsGeometryManager = o2::acts::TrackingGeometryManager::instance(); + const auto& geometryContext = actsGeometryManager.getNominalContext(); + + Acts::MagneticFieldContext magneticFieldContext; + Acts::CalibrationContext calibrationContext; + + // Initial implementation with a constant magnetic field + auto magneticField = std::make_shared( + Acts::Vector3{0., 0., mBz * Acts::UnitConstants::T}); + + using Propagator = Acts::Propagator, Acts::DirectNavigator>; + using Trajectory = Acts::VectorMultiTrajectory; + using TrackContainer = Acts::TrackContainer; + + Propagator propagator{ + Acts::EigenStepper<>{magneticField}, + Acts::DirectNavigator{}}; + + Acts::KalmanFitter fitter{std::move(propagator)}; + + struct HitMeasurement { + const Acts::Surface* surface{nullptr}; + Acts::Vector2 local = Acts::Vector2::Zero(); + Acts::SquareMatrix<2> covariance = Acts::SquareMatrix<2>::Zero(); + }; + + struct HitSourceLink { + const HitMeasurement* hit{nullptr}; + }; + + struct SurfaceAccessor { + const Acts::Surface* getSurface(const Acts::SourceLink& sourceLink) const + { + return sourceLink.get().hit->surface; + } + }; + + struct MeasurementCalibrator { + void calibrate(const Acts::GeometryContext&, + const Acts::CalibrationContext&, + const Acts::SourceLink& sourceLink, + Trajectory::TrackStateProxy trackState) const + { + const auto& hit = *sourceLink.get().hit; + + trackState.setUncalibratedSourceLink(Acts::SourceLink{sourceLink}); + trackState.allocateCalibrated(2); + trackState.calibrated<2>() = hit.local; + trackState.calibratedCovariance<2>() = hit.covariance; + + trackState.setProjectorSubspaceIndices(std::array{static_cast(Acts::eBoundLoc0), + static_cast(Acts::eBoundLoc1)}); + } + }; + + SurfaceAccessor surfaceAccessor; + MeasurementCalibrator calibrator; + Acts::GainMatrixUpdater updater; + Acts::GainMatrixSmoother smoother; + + Acts::KalmanFitterExtensions extensions; + extensions.surfaceAccessor.connect<&SurfaceAccessor::getSurface>(&surfaceAccessor); + extensions.calibrator.connect<&MeasurementCalibrator::calibrate>(&calibrator); + extensions.updater.connect<&Acts::GainMatrixUpdater::operator()>(&updater); + extensions.smoother.connect<&Acts::GainMatrixSmoother::operator()>(&smoother); + + Acts::PropagatorPlainOptions propagatorOptions{geometryContext, magneticFieldContext}; + + Acts::KalmanFitterOptions fitterOptions{geometryContext, + magneticFieldContext, + calibrationContext, + extensions, + propagatorOptions, + nullptr, + false, + false}; + + // Mutable TrackContainer uses reference-holder backends, not rvalue shared_ptrs. + Acts::VectorTrackContainer trackBackend; + Trajectory trajectoryBackend; + TrackContainer tracks{trackBackend, trajectoryBackend}; + + struct FittedCandidate { + Acts::BoundTrackParameters inner; + Acts::BoundTrackParameters outer; + SeedACTS seed; + float chi2; + }; + std::vector candidates; + int nFitted = 0; + int nFailed = 0; + int nInvalid = 0; + + for (const auto& seed : mSeeds) { + const std::array spacePoints{ + seed.bottom, seed.middle, seed.top}; + + if (!spacePoints[0] || !spacePoints[1] || !spacePoints[2]) { + ++nInvalid; + continue; + } + + if (!spacePoints[0]->surface || !spacePoints[1]->surface || + !spacePoints[2]->surface) { + ++nInvalid; + continue; + } - std::vector>> seedsPerGroup; - typename Acts::SeedFinder>::SeedingState seedingState; - seedingState.spacePointMutableData.resize(spContainer.size()); - const Acts::Range1D rMiddleSPRange; - for (auto [bottom, middle, top] : spGroup) { - auto& v = seedsPerGroup.emplace_back(); + // The direct fitter expects at most one measurement per surface. + if (spacePoints[0]->surface == spacePoints[1]->surface || + spacePoints[0]->surface == spacePoints[2]->surface || + spacePoints[1]->surface == spacePoints[2]->surface) { + ++nInvalid; + continue; + } + + Acts::Vector3 pos0{spacePoints[0]->x, spacePoints[0]->y, spacePoints[0]->z}; + Acts::Vector3 pos1{spacePoints[1]->x, spacePoints[1]->y, spacePoints[1]->z}; + Acts::Vector3 pos2{spacePoints[2]->x, spacePoints[2]->y, spacePoints[2]->z}; + + Acts::FreeVector freeParams; try { - seedFinder.createSeedsForGroup(finderOpts, seedingState, spGroup.grid(), v, bottom, middle, top, rMiddleSPRange); + freeParams = Acts::estimateTrackParamsFromSeed( + pos0, 0., pos1, pos2, + Acts::Vector3{0., 0., mBz * Acts::UnitConstants::T}); } catch (const std::exception& e) { - LOG(fatal) << "Error during seed finding for a group: " << e.what(); - return; + LOG(debug) << "ACTS seed parameter estimation failed: " << e.what(); + ++nInvalid; + continue; } - } - LOG(debug) << "Seed finding completed, found " << seedsPerGroup.size() << " groups with seeds"; - // Convert Acts seeds to the internal SeedACTS representation - for (const auto& groupSeeds : seedsPerGroup) { - for (const auto& actsSeed : groupSeeds) { - SeedACTS seed; - seed.bottom = &actsSeed.sp()[0]->externalSpacePoint(); - seed.middle = &actsSeed.sp()[1]->externalSpacePoint(); - seed.top = &actsSeed.sp()[2]->externalSpacePoint(); - seed.quality = actsSeed.seedQuality(); - mSeeds.push_back(seed); + const double qOverP = freeParams[Acts::eFreeQOverP]; + if (!std::isfinite(qOverP) || std::abs(qOverP) < 1.e-12) { + ++nInvalid; + continue; + } + + Acts::Vector3 direction{ + freeParams[Acts::eFreeDir0], + freeParams[Acts::eFreeDir1], + freeParams[Acts::eFreeDir2]}; + + if (!direction.allFinite() || direction.norm() < 1.e-12) { + ++nInvalid; + continue; } + direction.normalize(); + + std::array measurements; + std::vector sourceLinks; + std::vector surfaces; + sourceLinks.reserve(3); + surfaces.reserve(3); + + bool validMeasurements = true; + for (int i = 0; i < 3; ++i) { + const auto* sp = spacePoints[i]; + const auto* surface = sp->surface; + Acts::Vector3 position{sp->x, sp->y, sp->z}; + + auto localResult = surface->globalToLocal( + geometryContext, position, direction, + 1. * Acts::UnitConstants::mm); + + if (!localResult.ok()) { + validMeasurements = false; + break; + } + + auto& measurement = measurements[i]; + measurement.surface = surface; + measurement.local = localResult.value(); + + const double sigma = 100. * Acts::UnitConstants::um; + measurement.covariance.setZero(); + measurement.covariance(0, 0) = sigma * sigma; + measurement.covariance(1, 1) = sigma * sigma; + + sourceLinks.emplace_back(HitSourceLink{&measurement}); + surfaces.push_back(surface); + } + + if (!validMeasurements) { + ++nInvalid; + continue; + } + + Acts::Vector4 position4; + position4 << pos0.x(), pos0.y(), pos0.z(), 0.; + + Acts::BoundMatrix covariance = Acts::BoundMatrix::Zero(); + covariance(Acts::eBoundLoc0, Acts::eBoundLoc0) = + std::pow(1. * Acts::UnitConstants::mm, 2); + covariance(Acts::eBoundLoc1, Acts::eBoundLoc1) = + std::pow(1. * Acts::UnitConstants::mm, 2); + covariance(Acts::eBoundPhi, Acts::eBoundPhi) = 0.1 * 0.1; + covariance(Acts::eBoundTheta, Acts::eBoundTheta) = 0.1 * 0.1; + covariance(Acts::eBoundQOverP, Acts::eBoundQOverP) = + std::pow(0.5 * qOverP, 2); + covariance(Acts::eBoundTime, Acts::eBoundTime) = + std::pow(1. * Acts::UnitConstants::ns, 2); + + auto initialParameters = Acts::BoundTrackParameters::createCurvilinear( + position4, direction, qOverP, covariance, + Acts::ParticleHypothesis::pion()); + + auto result = fitter.fit( + sourceLinks.begin(), sourceLinks.end(), initialParameters, + fitterOptions, surfaces, tracks); + + if (!result.ok()) { + LOG(debug) << "ACTS Kalman fit failed: " << result.error().message(); + ++nFailed; + continue; + } + + const auto track = result.value(); + if (track.nMeasurements() < 3) { + ++nFailed; + continue; + } + + const auto innerState = track.innermostTrackState(); + if (!innerState.has_value()) { + ++nFailed; + continue; + } + const auto outerState = track.outermostTrackState(); + if (!innerState->hasSmoothed() || !outerState.hasSmoothed() || + !std::isfinite(track.chi2()) || track.chi2() < 0.f) { + ++nFailed; + continue; + } + + // Explicitly use the smoothed fit, rather than the original seed estimate. + typename Trajectory::ConstTrackStateProxy innerProxy{*innerState}; + typename Trajectory::ConstTrackStateProxy outerProxy{outerState}; + candidates.push_back({track.createParametersFromState(innerProxy), + track.createParametersFromState(outerProxy), + seed, static_cast(track.chi2())}); + ++nFitted; + LOG(debug) << "ACTS fitted seed: " + << "nMeasurements=" << track.nMeasurements() + << " chi2=" << track.chi2() + << " nDoF=" << track.nDoF(); } - LOGF(info, "Created %zu seeds from %zu space points", mSeeds.size(), mSpacePoints.size()); + // Accept lower-chi2 candidates first: multiple seeds can share clusters. + std::sort(candidates.begin(), candidates.end(), [](const auto& a, const auto& b) { + return a.chi2 < b.chi2; + }); + int nAccepted = 0; + for (const auto& candidate : candidates) { + nAccepted += acceptTrack(candidate.inner, candidate.outer, candidate.seed, candidate.chi2); + } + LOG(info) << "ACTS fitting: " << nFitted << " fitted, " + << nAccepted << " accepted, " << nFailed << " failed, " + << nInvalid << " invalid from " << mSeeds.size() << " seeds"; } template -bool TrackerACTS::estimateTrackParams(const SeedACTS& seed, o2::its::TrackITSExt& track) const +bool TrackerACTS::acceptTrack(const Acts::BoundTrackParameters& inner, + const Acts::BoundTrackParameters& outer, + const SeedACTS& seed, float chi2) { - return true; -} + if (!mTimeFrame || !std::isfinite(chi2) || chi2 < 0.f) { + return false; + } -template -void TrackerACTS::findTracks() -{ + const std::array hits{seed.bottom, seed.middle, seed.top}; + std::array clusterIndices; + clusterIndices.fill(-1); + std::array internalIndices; + internalIndices.fill(-1); + uint32_t pattern = 0; + std::array keys{}; + + for (int i = 0; i < 3; ++i) { + const auto* hit = hits[i]; + if (!hit || hit->layer < 0 || hit->layer >= nLayers || hit->clusterId < 0 || + hit->rof < 0 || hit->rof >= mTimeFrame->getNrof(hit->layer)) { + return false; + } + const uint32_t mask = 1u << hit->layer; + if (pattern & mask) { // the TrackITSExt index array has one slot per layer + return false; + } + pattern |= mask; + keys[i] = (uint64_t(static_cast(hit->layer)) << 32) | + static_cast(hit->clusterId); + if (mAcceptedClusters.contains(keys[i])) { + return false; + } + // Keep the internal index for later MC lookup; publish the external one. + internalIndices[hit->layer] = hit->clusterId; + clusterIndices[hit->layer] = mTimeFrame->getClusterExternalIndex(hit->layer, hit->clusterId); + if (clusterIndices[hit->layer] < 0) { + return false; + } + } + + o2::track::TrackParCov innerO2, outerO2; + const auto& gctx = o2::acts::TrackingGeometryManager::instance().getNominalContext(); + if (!convertTrackParameters(inner, gctx, innerO2) || !convertTrackParameters(outer, gctx, outerO2)) { + return false; + } + + // Intersect the ROF timing intervals exactly as in ITS track acceptance. + const auto& timing = mTimeFrame->getROFOverlapTableView(); + auto nominal = timing.getLayer(hits[0]->layer).getROFTimeBounds(hits[0]->rof); + auto expanded = timing.getLayer(hits[0]->layer).getROFTimeBounds(hits[0]->rof, true); + bool nominalCompatible = true; + for (int i = 1; i < 3; ++i) { + const auto other = timing.getLayer(hits[i]->layer).getROFTimeBounds(hits[i]->rof); + const auto otherExpanded = timing.getLayer(hits[i]->layer).getROFTimeBounds(hits[i]->rof, true); + if (nominalCompatible) { + if (nominal.isCompatible(other)) { + nominal += other; + } else { + nominalCompatible = false; + } + } + if (!expanded.isCompatible(otherExpanded)) { + return false; + } + expanded += otherExpanded; + } + + o2::its::TrackITSExt o2Track{std::move(innerO2), static_cast(hits.size()), chi2, std::move(outerO2), clusterIndices}; + o2Track.setPattern(pattern); + o2Track.getTimeStamp() = (nominalCompatible ? nominal : expanded).makeSymmetrical(); + const float maxHalfROF = timing.getClockLayer().mROFLength * 0.5f; + if (o2Track.getTimeStamp().getTimeStampError() > maxHalfROF) { + o2Track.getTimeStamp().setTimeStampError(maxHalfROF); + } + + // Labels are computed once, after all ROFs, from the internal indices. + mTimeFrame->getTracks().emplace_back(std::move(o2Track)); + mAcceptedTrackClusterIndices.push_back(internalIndices); + for (const auto key : keys) { + mAcceptedClusters.insert(key); + } + return true; } template void TrackerACTS::computeTracksMClabels() { + if (!mTimeFrame || !mTimeFrame->hasMCinformation()) { + return; + } + + auto& tracks = mTimeFrame->getTracks(); + auto& labels = mTimeFrame->getTracksLabel(); + + if (tracks.size() != mAcceptedTrackClusterIndices.size()) { + LOGP(fatal, "ACTS MC labeling: {} tracks but {} internal cluster reference sets", + tracks.size(), mAcceptedTrackClusterIndices.size()); + return; + } + + // There must be exactly one label for each accepted track, in track order. + labels.clear(); + labels.reserve(tracks.size()); + + struct LabelCount { + o2::MCCompLabel label; + int count; + int lastHit; + }; + + for (size_t iTrack = 0; iTrack < tracks.size(); ++iTrack) { + auto& track = tracks[iTrack]; + const auto& indices = mAcceptedTrackClusterIndices[iTrack]; + std::vector counts; + int nHits = 0; + + for (int layer = 0; layer < nLayers; ++layer) { + if (indices[layer] < 0) { + continue; + } + ++nHits; + for (const auto& candidate : mTimeFrame->getClusterLabels(layer, indices[layer])) { + if (!candidate.isValid()) { + continue; + } + + auto it = std::find_if(counts.begin(), counts.end(), [&candidate](const auto& entry) { + return entry.label == candidate; + }); + if (it == counts.end()) { + counts.push_back({candidate, 1, nHits}); + } else if (it->lastHit != nHits) { + // Count a given truth label at most once per attached cluster. + ++it->count; + it->lastHit = nHits; + } + } + } + o2::MCCompLabel winner; + uint32_t pattern = track.getPattern() & 0x0000ffffu; + if (counts.empty()) { + winner.setFakeFlag(); + } else { + const auto best = std::max_element(counts.begin(), counts.end(), [](const auto& a, const auto& b) { + return a.count < b.count; + }); + winner = best->label; + if (best->count != nHits) { + winner.setFakeFlag(); + } + // Preserve the standard TrackITS pattern convention: bits [16, 31] + // indicate clusters compatible with the representative MC label. + for (int layer = 0; layer < nLayers; ++layer) { + if (indices[layer] < 0) { + continue; + } + for (const auto& candidate : mTimeFrame->getClusterLabels(layer, indices[layer])) { + if (candidate.isValid() && candidate == best->label) { + pattern |= (1u << (16 + layer)); + break; + } + } + } + } + track.setPattern(pattern); + labels.emplace_back(winner); + } + LOGP(info, "ACTS MC labeling: {} tracks, {} labels", tracks.size(), labels.size()); } template @@ -260,35 +857,42 @@ void TrackerACTS::clustersToTracks() } double totalTime = 0.; + mAcceptedClusters.clear(); + mAcceptedTrackClusterIndices.clear(); LOG(info) << "==== TRK ACTS Tracking ===="; - LOG(info) << "Processing " << mTimeFrame->getNrof() << " ROFs with B = " << mBz << " T"; - // Process each ROF - for (int iROF = 0; iROF < mTimeFrame->getNrof(); ++iROF) { + auto& actsGeometryManager = o2::acts::TrackingGeometryManager::instance(); + const auto& trackingGeometry = actsGeometryManager.get(); + (void)trackingGeometry; + + auto* trkGeometry = o2::trk::GeometryTGeo::Instance(); + trkGeometry->fillMatrixCache(o2::math_utils::bit2Mask(o2::math_utils::TransformType::L2G)); + const auto& surfaceIndex = actsGeometryManager.getIndex(*trkGeometry, 1.e-3, [trkGeometry](int sensorID) { + return std::string(trkGeometry->getMatrixPath(sensorID).Data()); + }); + + LOG(info) << "ACTS tracking geometry initialized with " << surfaceIndex.size() + << " mapped sensors on " << surfaceIndex.getNSurfaces() << " surfaces"; + + for (int iROF = 0; iROF < mTimeFrame->getNrof(0); ++iROF) { LOG(info) << "Processing ROF " << iROF; - // Build space points + mCurState = SpacePointBuilding; - totalTime += evaluateTask([this, iROF]() { buildSpacePoints(iROF); }, - StateNames[mCurState]); + totalTime += evaluateTask([this, iROF]() { buildSpacePoints(iROF); }, StateNames[mCurState]); - // Run seeding mCurState = Seeding; - totalTime += evaluateTask([this]() { createSeeds(); }, - StateNames[mCurState]); + totalTime += evaluateTask([this]() { createSeeds(); }, StateNames[mCurState]); - // Find tracks mCurState = TrackFinding; - totalTime += evaluateTask([this]() { findTracks(); }, - StateNames[mCurState]); + totalTime += evaluateTask([this]() { findTracks(); }, StateNames[mCurState]); } - // MC labeling + // Run MC labeling once, after every accepted track has been published. if (mTimeFrame->hasMCinformation()) { computeTracksMClabels(); } LOG(info) << "=== TimeFrame " << mTimeFrameCounter << " completed in: " << totalTime << " ms ==="; - ++mTimeFrameCounter; mTotalTime += totalTime; } @@ -297,8 +901,7 @@ template void TrackerACTS::printSummary() const { float avgTF = mTimeFrameCounter > 0 ? static_cast(mTotalTime) / mTimeFrameCounter : 0.f; - LOGP(info, "TrackerACTS summary: Processed {} TFs in TOT={:.2f} ms, AVG/TF={:.2f} ms", - mTimeFrameCounter, mTotalTime, avgTF); + LOGP(info, "TrackerACTS summary: Processed {} TFs in TOT={:.2f} ms, AVG/TF={:.2f} ms", mTimeFrameCounter, mTotalTime, avgTF); } // Explicit template instantiations diff --git a/Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/src/TrackerACTS2.cxx b/Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/src/TrackerACTS2.cxx new file mode 100644 index 0000000000000..17744dff61efb --- /dev/null +++ b/Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/src/TrackerACTS2.cxx @@ -0,0 +1,1149 @@ +// 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. + +/// +/// \file TrackerACTS.cxx +/// \brief TRK tracker using ACTS seeding and track finding +/// \author Nicolò Jacazio, Università del Piemonte Orientale (IT) +/// \since 2026-04-01 +/// + +#include "ALICE3GlobalReconstruction/TrackerACTS.h" + +#include "ACTSInterface/TrackingGeometryManager.h" +#include "MathUtils/Utils.h" +#include "TRKBase/GeometryTGeo.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "ACTSInterface/TrackingGeometryManager.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace o2::trk +{ + +namespace +{ +// Convert ACTS surface-bound parameters to the O2 tracking frame +// (x, alpha; y, z, sin(phi-alpha), tan(lambda), q/pT). +// The covariance is transformed by a numerical Jacobian evaluated at the +// fitted state. Position units are converted from ACTS mm to O2 cm. +bool convertTrackParameters(const Acts::BoundTrackParameters& actsTrack, + const Acts::GeometryContext& gctx, + o2::track::TrackParCov& output) +{ + const auto& state = actsTrack.parameters(); + if (!state.allFinite() || !actsTrack.covariance().has_value() || + !actsTrack.covariance()->allFinite()) { + return false; + } + + const auto position = actsTrack.position(gctx); + if (!position.allFinite()) { + return false; + } + const double alpha = std::atan2(position.y(), position.x()); + const double ca = std::cos(alpha); + const double sa = std::sin(alpha); + + auto toO2 = [&](const Acts::BoundVector& p, Eigen::Matrix& dest) -> bool { + const double theta = p[Acts::eBoundTheta]; + const double sinTheta = std::sin(theta); + if (!std::isfinite(sinTheta) || sinTheta <= 1.e-6) { + return false; + } + const auto direction = Acts::makeDirectionFromPhiTheta(p[Acts::eBoundPhi], theta); + const Acts::Vector2 local = p.head<2>(); + const auto xyz = actsTrack.referenceSurface().localToGlobal(gctx, local, direction); + if (!xyz.allFinite()) { + return false; + } + + const double snp = std::sin(p[Acts::eBoundPhi] - alpha); + if (!std::isfinite(snp) || std::abs(snp) >= 0.999999) { + return false; + } + + dest << (-xyz.x() * sa + xyz.y() * ca) / Acts::UnitConstants::cm, + xyz.z() / Acts::UnitConstants::cm, + snp, + std::cos(theta) / sinTheta, + p[Acts::eBoundQOverP] * Acts::UnitConstants::GeV / sinTheta; + return dest.allFinite(); + }; + + Eigen::Matrix pars; + if (!toO2(state, pars)) { + return false; + } + + const double x = (position.x() * ca + position.y() * sa) / Acts::UnitConstants::cm; + if (!std::isfinite(x)) { + return false; + } + + Eigen::Matrix jacobian; + jacobian.setZero(); + for (int i = 0; i < Acts::eBoundSize; ++i) { + if (i == Acts::eBoundTime) { + continue; // O2 TrackParCov has no time coordinate + } + const double step = (i == Acts::eBoundLoc0 || i == Acts::eBoundLoc1) + ? 1.e-3 * Acts::UnitConstants::mm + : (i == Acts::eBoundQOverP) + ? std::max(std::abs(state[i]) * 1.e-4, 1.e-6 / Acts::UnitConstants::GeV) + : 1.e-5; + Acts::BoundVector plus = state; + Acts::BoundVector minus = state; + plus[i] += step; + minus[i] -= step; + Eigen::Matrix fPlus, fMinus; + if (!toO2(plus, fPlus) || !toO2(minus, fMinus)) { + return false; + } + jacobian.col(i) = (fPlus - fMinus) / (2. * step); + } + + const Eigen::Matrix covariance = + jacobian * (*actsTrack.covariance()) * jacobian.transpose(); + if (!covariance.allFinite()) { + return false; + } + + o2::track::TrackParCov::params_t o2Params{}; + o2::track::TrackParCov::covMat_t o2Cov{}; + for (int i = 0; i < 5; ++i) { + if (!(covariance(i, i) > 0.) || !std::isfinite(pars[i])) { + return false; + } + o2Params[i] = static_cast(pars[i]); + for (int j = 0; j <= i; ++j) { + o2Cov[o2::track::CovarMap[i][j]] = + static_cast(0.5 * (covariance(i, j) + covariance(j, i))); + } + } + const int charge = state[Acts::eBoundQOverP] >= 0. ? 1 : -1; + output = o2::track::TrackParCov{static_cast(x), static_cast(alpha), + o2Params, o2Cov, charge, o2::track::PID::Pion}; + return output.isValid(); +} +} // namespace +template +TrackerACTS::TrackerACTS() +{ + // Initialize space points storage per layer + mSpacePointsPerLayer.resize(nLayers); + // The debug file and histogram are intentionally not allocated here. +} + +template +TrackerACTS::~TrackerACTS() +{ + if (mDebugFile) { + mDebugFile->cd(); + if (mHistSpacePoints) { + mHistSpacePoints->Write("", TObject::kOverwrite); + } + mDebugFile->Write("", TObject::kOverwrite); + mDebugFile->Close(); + } +} + +template +bool TrackerACTS::isDebugTimeFrame() const +{ + return mDebugFile && (mConfig.debugMaxTimeFrames < 0 || + mTimeFrameCounter < static_cast(mConfig.debugMaxTimeFrames)); +} + +template +void TrackerACTS::initDebugOutput() +{ + if (!mConfig.writeDebugTrees || mDebugFile || mConfig.debugMaxTimeFrames == 0) { + return; + } + mDebugFile.reset(TFile::Open(mConfig.debugFileName.c_str(), "RECREATE")); + if (!mDebugFile || mDebugFile->IsZombie()) { + LOG(error) << "Cannot create ACTS debug ROOT file " << mConfig.debugFileName; + mDebugFile.reset(); + return; + } + mDebugFile->cd(); + mClusterTree = new TTree("clusters", "ACTS input clusters and surface mapping; positions in cm; variances in ACTS units"); + mClusterTree->Branch("tf", &mDebugCluster.tf, "tf/I"); + mClusterTree->Branch("rof", &mDebugCluster.rof, "rof/I"); + mClusterTree->Branch("layer", &mDebugCluster.layer, "layer/I"); + mClusterTree->Branch("clusterId", &mDebugCluster.clusterId, "clusterId/I"); + mClusterTree->Branch("externalId", &mDebugCluster.externalId, "externalId/I"); + mClusterTree->Branch("sensorId", &mDebugCluster.sensorId, "sensorId/I"); + mClusterTree->Branch("mappingStatus", &mDebugCluster.mappingStatus, "mappingStatus/I"); + mClusterTree->Branch("geometryId", &mDebugCluster.geometryId, "geometryId/l"); + mClusterTree->Branch("x", &mDebugCluster.x, "x/F"); + mClusterTree->Branch("y", &mDebugCluster.y, "y/F"); + mClusterTree->Branch("z", &mDebugCluster.z, "z/F"); + mClusterTree->Branch("r", &mDebugCluster.r, "r/F"); + mClusterTree->Branch("phi", &mDebugCluster.phi, "phi/F"); + mClusterTree->Branch("varianceR", &mDebugCluster.varianceR, "varianceR/F"); + mClusterTree->Branch("varianceZ", &mDebugCluster.varianceZ, "varianceZ/F"); + + mSeedTree = new TTree("seeds", "ACTS output seed triplets; positions in cm; bottom/middle/top order"); + mSeedTree->Branch("tf", &mDebugSeed.tf, "tf/I"); + mSeedTree->Branch("rof", &mDebugSeed.rof, "rof/I"); + mSeedTree->Branch("seedId", &mDebugSeed.seedId, "seedId/I"); + mSeedTree->Branch("quality", &mDebugSeed.quality, "quality/F"); + mSeedTree->Branch("layer", mDebugSeed.layer, "layer[3]/I"); + mSeedTree->Branch("clusterId", mDebugSeed.clusterId, "clusterId[3]/I"); + mSeedTree->Branch("sensorId", mDebugSeed.sensorId, "sensorId[3]/I"); + mSeedTree->Branch("geometryId", mDebugSeed.geometryId, "geometryId[3]/l"); + mSeedTree->Branch("x", mDebugSeed.x, "x[3]/F"); + mSeedTree->Branch("y", mDebugSeed.y, "y[3]/F"); + mSeedTree->Branch("z", mDebugSeed.z, "z[3]/F"); + mSeedTree->Branch("r", mDebugSeed.r, "r[3]/F"); + mSeedTree->Branch("phi", mDebugSeed.phi, "phi[3]/F"); + + mROFTree = new TTree("rofSummary", "ACTS reconstruction counts per time frame and ROF"); + mROFTree->Branch("tf", &mDebugROF.tf, "tf/I"); + mROFTree->Branch("rof", &mDebugROF.rof, "rof/I"); + mROFTree->Branch("nSpacePoints", &mDebugROF.nSpacePoints, "nSpacePoints/I"); + mROFTree->Branch("nSeeds", &mDebugROF.nSeeds, "nSeeds/I"); + mROFTree->Branch("nAcceptedTracks", &mDebugROF.nAcceptedTracks, "nAcceptedTracks/I"); + + mHistSpacePoints = std::make_unique("hSpacePoints", "Space points; x (cm); y (cm)", + 200, -100, 100, 200, -100, 100); + mHistSpacePoints->SetDirectory(nullptr); // Owned by TrackerACTS, not by TFile. + LOG(info) << "ACTS debug trees enabled: " << mConfig.debugFileName; +} + +template +void TrackerACTS::dumpGeometryOnce() +{ + if (!mConfig.dumpTrackingGeometry || mGeometryDumped) { + return; + } + mGeometryDumped = true; // Even on I/O failure, do not retry on every TF. + auto& manager = o2::acts::TrackingGeometryManager::instance(); + const auto& geometry = manager.get(); + const auto& gctx = manager.getNominalContext(); + auto* trkGeometry = o2::trk::GeometryTGeo::Instance(); + const auto& surfaceIndex = manager.getIndex(*trkGeometry, 1.e-3, [trkGeometry](int sensorID) { + return std::string(trkGeometry->getMatrixPath(sensorID).Data()); + }); + + std::ofstream out(mConfig.geometryDumpFileName); + if (!out) { + LOG(error) << "Cannot open ACTS geometry dump: " << mConfig.geometryDumpFileName; + return; + } + out << "{\n \"description\": \"ACTS volumes and surfaces snapshot, not a complete ACTS geometry serialization\",\n" + " \"lengthUnits\": \"mm (angles in rad)\",\n \"volumes\": [\n"; + bool first = true; + int nVolumes = 0; + geometry.visitVolumes([&](const Acts::TrackingVolume& volume) { + const auto& transform = volume.localToGlobalTransform(gctx); + nlohmann::json item; + item["geometryId"] = volume.geometryId().value(); + item["motherGeometryId"] = volume.motherVolume() ? volume.motherVolume()->geometryId().value() : 0; + item["name"] = volume.volumeName(); + item["boundsType"] = static_cast(volume.volumeBounds().type()); + item["boundsValues"] = volume.volumeBounds().values(); + item["hasMaterial"] = volume.volumeMaterial() != nullptr; + item["translation"] = {transform.translation().x(), transform.translation().y(), transform.translation().z()}; + item["rotation"] = nlohmann::json::array(); + for (int row = 0; row < 3; ++row) { + item["rotation"].push_back({transform.rotation()(row, 0), transform.rotation()(row, 1), transform.rotation()(row, 2)}); + } + if (!first) { + out << ",\n"; + } + out << item.dump(); + first = false; + ++nVolumes; + }); + out << "\n ],\n \"surfaces\": [\n"; + first = true; + int nSurfaces = 0; + geometry.visitSurfaces([&](const Acts::Surface& surface) { + const auto& transform = surface.localToGlobalTransform(gctx); + nlohmann::json item; + item["geometryId"] = surface.geometryId().value(); + item["surfaceType"] = static_cast(surface.type()); + item["boundsType"] = static_cast(surface.bounds().type()); + item["boundsValues"] = surface.bounds().values(); + item["hasMaterial"] = surface.hasMaterial(); + item["translation"] = {transform.translation().x(), transform.translation().y(), transform.translation().z()}; + item["rotation"] = nlohmann::json::array(); + for (int row = 0; row < 3; ++row) { + item["rotation"].push_back({transform.rotation()(row, 0), transform.rotation()(row, 1), transform.rotation()(row, 2)}); + } + item["sensorIds"] = surfaceIndex.getSensorIDs(surface.geometryId()); + if (!first) { + out << ",\n"; + } + out << item.dump(); + first = false; + ++nSurfaces; + }, false); // Also include passive/non-sensitive surfaces where reachable. + out << "\n ]\n}\n"; + if (!out.good()) { + LOG(error) << "Error while writing ACTS geometry dump " << mConfig.geometryDumpFileName; + return; + } + LOG(info) << "ACTS tracking geometry snapshot: " << nVolumes << " volumes and " + << nSurfaces << " surfaces written to " << mConfig.geometryDumpFileName; +} + +template +void TrackerACTS::fillDebugSeeds(int rof) +{ + if (!isDebugTimeFrame()) { + return; + } + for (size_t iSeed = 0; iSeed < mSeeds.size(); ++iSeed) { + const auto& seed = mSeeds[iSeed]; + mDebugSeed = {}; + mDebugSeed.tf = static_cast(mTimeFrameCounter); + mDebugSeed.rof = rof; + mDebugSeed.seedId = static_cast(iSeed); + mDebugSeed.quality = seed.quality; + const std::array points{seed.bottom, seed.middle, seed.top}; + for (int i = 0; i < 3; ++i) { + const auto& sp = *points[i]; + mDebugSeed.layer[i] = sp.layer; + mDebugSeed.clusterId[i] = sp.clusterId; + mDebugSeed.sensorId[i] = sp.sensorId; + mDebugSeed.geometryId[i] = sp.surface ? sp.surface->geometryId().value() : 0; + mDebugSeed.x[i] = sp.x / Acts::UnitConstants::cm; + mDebugSeed.y[i] = sp.y / Acts::UnitConstants::cm; + mDebugSeed.z[i] = sp.z / Acts::UnitConstants::cm; + mDebugSeed.r[i] = sp.r() / Acts::UnitConstants::cm; + mDebugSeed.phi[i] = sp.phi(); + } + mSeedTree->Fill(); + } +} + +template +void TrackerACTS::adoptTimeFrame(o2::its::TimeFrame& tf) +{ + mTimeFrame = &tf; +} + +template +void TrackerACTS::buildSpacePoints(int rof) +{ + mSpacePoints.clear(); + for (auto& layerSPs : mSpacePointsPerLayer) { + layerSPs.clear(); + } + + if (!mSensorIDs) { + LOG(error) << "Cannot build ACTS space points: Sensor IDs not adopted"; + return; + } + + auto& actsGeometryManager = o2::acts::TrackingGeometryManager::instance(); + auto* trkGeometry = o2::trk::GeometryTGeo::Instance(); + + trkGeometry->fillMatrixCache(o2::math_utils::bit2Mask(o2::math_utils::TransformType::L2G)); + + const auto& surfaceIndex = actsGeometryManager.getIndex(*trkGeometry, 1.e-3, [trkGeometry](int sensorID) { + return std::string(trkGeometry->getMatrixPath(sensorID).Data()); + }); + + int nMapped = 0; + int nMissing = 0; + + for (int layer = 0; layer < nLayers; ++layer) { + const auto clusters = mTimeFrame->getUnsortedClustersOnLayer(rof, layer); + const int firstCluster = mTimeFrame->getSortedStartIndex(rof, layer); + + LOG(debug) << "ACTSTracker: got " << clusters.size() << " clusters on layer " << layer; + + for (size_t iCluster = 0; iCluster < clusters.size(); ++iCluster) { + const auto& cluster = clusters[iCluster]; + const int clusterIndex = firstCluster + static_cast(iCluster); + + const bool hasSensorID = clusterIndex >= 0 && static_cast(clusterIndex) < (*mSensorIDs)[layer].size(); + const int sensorID = hasSensorID ? (*mSensorIDs)[layer][clusterIndex] : -1; + const auto* surface = hasSensorID ? surfaceIndex.getSurface(sensorID) : nullptr; + if (isDebugTimeFrame()) { + mDebugCluster = {}; + mDebugCluster.tf = static_cast(mTimeFrameCounter); + mDebugCluster.rof = rof; + mDebugCluster.layer = layer; + mDebugCluster.clusterId = clusterIndex; + mDebugCluster.externalId = mTimeFrame->getClusterExternalIndex(layer, clusterIndex); + mDebugCluster.sensorId = sensorID; + // 0: sensor ID missing; 1: ACTS surface missing; 2: mapped. + mDebugCluster.mappingStatus = !hasSensorID ? 0 : !surface ? 1 : 2; + mDebugCluster.geometryId = surface ? surface->geometryId().value() : 0; + mDebugCluster.x = cluster.xCoordinate; + mDebugCluster.y = cluster.yCoordinate; + mDebugCluster.z = cluster.zCoordinate; + mDebugCluster.r = std::hypot(mDebugCluster.x, mDebugCluster.y); + mDebugCluster.phi = std::atan2(mDebugCluster.y, mDebugCluster.x); + mDebugCluster.varianceR = 0.01f; + mDebugCluster.varianceZ = 0.01f; + mClusterTree->Fill(); + } + if (!hasSensorID) { + LOG(error) << "Missing sensor ID for layer " << layer << " cluster " << clusterIndex; + ++nMissing; + continue; + } + + if (!surface) { + LOG(error) << "No ACTS surface for sensor " << sensorID << " on layer " << layer; + ++nMissing; + continue; + } + + SpacePoint sp; + sp.x = cluster.xCoordinate * Acts::UnitConstants::cm; + sp.y = cluster.yCoordinate * Acts::UnitConstants::cm; + sp.z = cluster.zCoordinate * Acts::UnitConstants::cm; + + if (mHistSpacePoints) { + mHistSpacePoints->Fill(sp.x / Acts::UnitConstants::cm, sp.y / Acts::UnitConstants::cm); + } + + sp.layer = layer; + sp.clusterId = clusterIndex; + sp.rof = rof; + sp.sensorId = sensorID; + sp.surface = surface; + sp.varianceR = 0.01f; + sp.varianceZ = 0.01f; + + mSpacePoints.push_back(sp); + ++nMapped; + } + } + + for (auto& sp : mSpacePoints) { + if (sp.layer >= 0 && sp.layer < nLayers) { + mSpacePointsPerLayer[sp.layer].push_back(&sp); + } + } + + LOG(info) << "ACTS surface mapping: " << nMapped << " clusters mapped, " << nMissing << " clusters missing"; +} + +template +void TrackerACTS::createSeeds() +{ + mSeeds.clear(); + if (mSpacePoints.empty()) { + LOG(info) << "No space points available for seeding"; + return; + } + + const float bFieldInZ = mBz * Acts::UnitConstants::T; + + Acts::CylindricalSpacePointGrid::Config gridConfig; + gridConfig.minPt = mConfig.minPt; + gridConfig.rMin = 0.f; + gridConfig.rMax = 100.f * Acts::UnitConstants::cm; + gridConfig.zMin = mConfig.zMin; + gridConfig.zMax = mConfig.zMax; + gridConfig.deltaRMax = std::max(mConfig.deltaRMaxBottom, mConfig.deltaRMaxTop); + gridConfig.cotThetaMax = mConfig.cotThetaMax; + gridConfig.impactMax = mConfig.maxImpactParameter; + gridConfig.bFieldInZ = bFieldInZ; + gridConfig.bottomBinFinder = Acts::GridBinFinder<3ul>{1, std::vector>{}, 0}; + gridConfig.topBinFinder = Acts::GridBinFinder<3ul>{1, std::vector>{}, 0}; + + Acts::CylindricalSpacePointGrid grid{gridConfig}; + for (Acts::SpacePointIndex i = 0; i < mSpacePoints.size(); ++i) { + const auto& sp = mSpacePoints[i]; + grid.insert(i, sp.phi(), sp.z, sp.r()); + } + + // The doublet finder expects each bin to be sorted in increasing radius. + for (size_t i = 0; i < grid.numberOfBins(); ++i) { + auto& bin = grid.at(i); + std::sort(bin.begin(), bin.end(), [this](Acts::SpacePointIndex a, Acts::SpacePointIndex b) { + return mSpacePoints[a].r() < mSpacePoints[b].r(); + }); + } + + // ACTS 47 uses a non-templated, column-based space-point container. + Acts::SpacePointContainer spacePoints{ + Acts::SpacePointColumns::CopiedFromIndex | + Acts::SpacePointColumns::PackedXY | + Acts::SpacePointColumns::PackedZR | + Acts::SpacePointColumns::VarianceR | + Acts::SpacePointColumns::VarianceZ}; + spacePoints.reserve(grid.numberOfSpacePoints()); + + std::vector gridRanges; + gridRanges.reserve(grid.numberOfBins()); + for (size_t i = 0; i < grid.numberOfBins(); ++i) { + const auto first = static_cast(spacePoints.size()); + for (Acts::SpacePointIndex spIndex : grid.at(i)) { + const auto& sp = mSpacePoints[spIndex]; + auto proxy = spacePoints.createSpacePoint(); + proxy.copiedFromIndex() = spIndex; + proxy.xy() = std::array{sp.x, sp.y}; + proxy.zr() = std::array{sp.z, sp.r()}; + proxy.varianceR() = sp.varianceR; + proxy.varianceZ() = sp.varianceZ; + } + const auto last = static_cast(spacePoints.size()); + gridRanges.emplace_back(first, last); + } + + Acts::DoubletSeedFinder::Config bottomConfig; + bottomConfig.spacePointsSortedByRadius = true; + bottomConfig.candidateDirection = Acts::Direction::Backward(); + bottomConfig.deltaRMin = mConfig.deltaRMinBottom; + bottomConfig.deltaRMax = mConfig.deltaRMaxBottom; + bottomConfig.collisionRegionMin = mConfig.collisionRegionMin; + bottomConfig.collisionRegionMax = mConfig.collisionRegionMax; + bottomConfig.cotThetaMax = mConfig.cotThetaMax; + bottomConfig.impactMax = mConfig.maxImpactParameter; + bottomConfig.minPt = mConfig.minPt; + auto bottomFinder = Acts::DoubletSeedFinder::create( + Acts::DoubletSeedFinder::DerivedConfig{bottomConfig, bFieldInZ}); + + Acts::DoubletSeedFinder::Config topConfig = bottomConfig; + topConfig.candidateDirection = Acts::Direction::Forward(); + topConfig.deltaRMin = mConfig.deltaRMinTop; + topConfig.deltaRMax = mConfig.deltaRMaxTop; + auto topFinder = Acts::DoubletSeedFinder::create( + Acts::DoubletSeedFinder::DerivedConfig{topConfig, bFieldInZ}); + + Acts::TripletSeedFinder::Config tripletConfig; + tripletConfig.minPt = mConfig.minPt; + tripletConfig.impactMax = mConfig.maxImpactParameter; + tripletConfig.sigmaScattering = 5.f; + tripletConfig.radLengthPerSeed = 0.05f; + auto tripletFinder = Acts::TripletSeedFinder::create( + Acts::TripletSeedFinder::DerivedConfig{tripletConfig, bFieldInZ}); + + Acts::BroadTripletSeedFilter::Config filterConfig; + filterConfig.maxSeedsPerSpM = static_cast(mConfig.maxSeedsPerMiddleSP); + filterConfig.deltaRMin = std::min(mConfig.deltaRMinBottom, mConfig.deltaRMinTop); + Acts::BroadTripletSeedFilter::State filterState; + Acts::BroadTripletSeedFilter::Cache filterCache; + Acts::BroadTripletSeedFilter filter{filterConfig, filterState, filterCache, Acts::getDummyLogger()}; + + Acts::TripletSeeder seeder; + Acts::TripletSeeder::Cache seedingCache; + Acts::SeedContainer seeds; + seeds.assignSpacePointContainer(spacePoints); + + std::vector bottomRanges; + std::vector topRanges; + for (const auto [bottom, middle, top] : grid.binnedGroup()) { + bottomRanges.clear(); + for (const auto bin : bottom) { + bottomRanges.push_back(spacePoints.range(gridRanges.at(bin)).asConst()); + } + topRanges.clear(); + for (const auto bin : top) { + topRanges.push_back(spacePoints.range(gridRanges.at(bin)).asConst()); + } + auto middleRange = spacePoints.range(gridRanges.at(middle)).asConst(); + if (middleRange.empty()) { + continue; + } + + seeder.createSeedsFromGroups(seedingCache, *bottomFinder, *topFinder, + *tripletFinder, filter, spacePoints, + bottomRanges, middleRange, topRanges, + {gridConfig.rMin, gridConfig.rMax}, seeds); + } + + // Restore the existing SeedACTS representation using original SP indices. + mSeeds.reserve(seeds.size()); + for (const auto seed : seeds) { + const auto indices = seed.spacePointIndices(); + SeedACTS output; + output.bottom = &mSpacePoints[spacePoints.at(indices[0]).copiedFromIndex()]; + output.middle = &mSpacePoints[spacePoints.at(indices[1]).copiedFromIndex()]; + output.top = &mSpacePoints[spacePoints.at(indices[2]).copiedFromIndex()]; + output.quality = seed.quality(); + mSeeds.push_back(output); + } + fillDebugSeeds(mSpacePoints.front().rof); + + LOGF(info, "Created %zu seeds from %zu space points", mSeeds.size(), mSpacePoints.size()); +} + +template +void TrackerACTS::findTracks() +{ + if (mSeeds.empty()) { + return; + } + + auto& actsGeometryManager = o2::acts::TrackingGeometryManager::instance(); + const auto& geometryContext = actsGeometryManager.getNominalContext(); + + Acts::MagneticFieldContext magneticFieldContext; + Acts::CalibrationContext calibrationContext; + + // Initial implementation with a constant magnetic field + auto magneticField = std::make_shared( + Acts::Vector3{0., 0., mBz * Acts::UnitConstants::T}); + + using Propagator = Acts::Propagator, Acts::DirectNavigator>; + using Trajectory = Acts::VectorMultiTrajectory; + using TrackContainer = Acts::TrackContainer; + + Propagator propagator{ + Acts::EigenStepper<>{magneticField}, + Acts::DirectNavigator{}}; + + Acts::KalmanFitter fitter{std::move(propagator)}; + + struct HitMeasurement { + const Acts::Surface* surface{nullptr}; + Acts::Vector2 local = Acts::Vector2::Zero(); + Acts::SquareMatrix<2> covariance = Acts::SquareMatrix<2>::Zero(); + }; + + struct HitSourceLink { + const HitMeasurement* hit{nullptr}; + }; + + struct SurfaceAccessor { + const Acts::Surface* getSurface(const Acts::SourceLink& sourceLink) const + { + return sourceLink.get().hit->surface; + } + }; + + struct MeasurementCalibrator { + void calibrate(const Acts::GeometryContext&, + const Acts::CalibrationContext&, + const Acts::SourceLink& sourceLink, + Trajectory::TrackStateProxy trackState) const + { + const auto& hit = *sourceLink.get().hit; + + trackState.setUncalibratedSourceLink(Acts::SourceLink{sourceLink}); + trackState.allocateCalibrated(2); + trackState.calibrated<2>() = hit.local; + trackState.calibratedCovariance<2>() = hit.covariance; + + trackState.setProjectorSubspaceIndices(std::array{static_cast(Acts::eBoundLoc0), + static_cast(Acts::eBoundLoc1)}); + } + }; + + SurfaceAccessor surfaceAccessor; + MeasurementCalibrator calibrator; + Acts::GainMatrixUpdater updater; + Acts::GainMatrixSmoother smoother; + + Acts::KalmanFitterExtensions extensions; + extensions.surfaceAccessor.connect<&SurfaceAccessor::getSurface>(&surfaceAccessor); + extensions.calibrator.connect<&MeasurementCalibrator::calibrate>(&calibrator); + extensions.updater.connect<&Acts::GainMatrixUpdater::operator()>(&updater); + extensions.smoother.connect<&Acts::GainMatrixSmoother::operator()>(&smoother); + + Acts::PropagatorPlainOptions propagatorOptions{geometryContext, magneticFieldContext}; + + Acts::KalmanFitterOptions fitterOptions{geometryContext, + magneticFieldContext, + calibrationContext, + extensions, + propagatorOptions, + nullptr, + false, + false}; + + // Mutable TrackContainer uses reference-holder backends, not rvalue shared_ptrs. + Acts::VectorTrackContainer trackBackend; + Trajectory trajectoryBackend; + TrackContainer tracks{trackBackend, trajectoryBackend}; + + struct FittedCandidate { + Acts::BoundTrackParameters inner; + Acts::BoundTrackParameters outer; + SeedACTS seed; + float chi2; + }; + std::vector candidates; + int nFitted = 0; + int nFailed = 0; + int nInvalid = 0; + + for (const auto& seed : mSeeds) { + const std::array spacePoints{ + seed.bottom, seed.middle, seed.top}; + + if (!spacePoints[0] || !spacePoints[1] || !spacePoints[2]) { + ++nInvalid; + continue; + } + + if (!spacePoints[0]->surface || !spacePoints[1]->surface || + !spacePoints[2]->surface) { + ++nInvalid; + continue; + } + + // The direct fitter expects at most one measurement per surface. + if (spacePoints[0]->surface == spacePoints[1]->surface || + spacePoints[0]->surface == spacePoints[2]->surface || + spacePoints[1]->surface == spacePoints[2]->surface) { + ++nInvalid; + continue; + } + + Acts::Vector3 pos0{spacePoints[0]->x, spacePoints[0]->y, spacePoints[0]->z}; + Acts::Vector3 pos1{spacePoints[1]->x, spacePoints[1]->y, spacePoints[1]->z}; + Acts::Vector3 pos2{spacePoints[2]->x, spacePoints[2]->y, spacePoints[2]->z}; + + Acts::FreeVector freeParams; + try { + freeParams = Acts::estimateTrackParamsFromSeed( + pos0, 0., pos1, pos2, + Acts::Vector3{0., 0., mBz * Acts::UnitConstants::T}); + } catch (const std::exception& e) { + LOG(debug) << "ACTS seed parameter estimation failed: " << e.what(); + ++nInvalid; + continue; + } + + const double qOverP = freeParams[Acts::eFreeQOverP]; + if (!std::isfinite(qOverP) || std::abs(qOverP) < 1.e-12) { + ++nInvalid; + continue; + } + + Acts::Vector3 direction{ + freeParams[Acts::eFreeDir0], + freeParams[Acts::eFreeDir1], + freeParams[Acts::eFreeDir2]}; + + if (!direction.allFinite() || direction.norm() < 1.e-12) { + ++nInvalid; + continue; + } + direction.normalize(); + + std::array measurements; + std::vector sourceLinks; + std::vector surfaces; + sourceLinks.reserve(3); + surfaces.reserve(3); + + bool validMeasurements = true; + for (int i = 0; i < 3; ++i) { + const auto* sp = spacePoints[i]; + const auto* surface = sp->surface; + Acts::Vector3 position{sp->x, sp->y, sp->z}; + + auto localResult = surface->globalToLocal( + geometryContext, position, direction, + 1. * Acts::UnitConstants::mm); + + if (!localResult.ok()) { + validMeasurements = false; + break; + } + + auto& measurement = measurements[i]; + measurement.surface = surface; + measurement.local = localResult.value(); + + const double sigma = 100. * Acts::UnitConstants::um; + measurement.covariance.setZero(); + measurement.covariance(0, 0) = sigma * sigma; + measurement.covariance(1, 1) = sigma * sigma; + + sourceLinks.emplace_back(HitSourceLink{&measurement}); + surfaces.push_back(surface); + } + + if (!validMeasurements) { + ++nInvalid; + continue; + } + + Acts::Vector4 position4; + position4 << pos0.x(), pos0.y(), pos0.z(), 0.; + + Acts::BoundMatrix covariance = Acts::BoundMatrix::Zero(); + covariance(Acts::eBoundLoc0, Acts::eBoundLoc0) = + std::pow(1. * Acts::UnitConstants::mm, 2); + covariance(Acts::eBoundLoc1, Acts::eBoundLoc1) = + std::pow(1. * Acts::UnitConstants::mm, 2); + covariance(Acts::eBoundPhi, Acts::eBoundPhi) = 0.1 * 0.1; + covariance(Acts::eBoundTheta, Acts::eBoundTheta) = 0.1 * 0.1; + covariance(Acts::eBoundQOverP, Acts::eBoundQOverP) = + std::pow(0.5 * qOverP, 2); + covariance(Acts::eBoundTime, Acts::eBoundTime) = + std::pow(1. * Acts::UnitConstants::ns, 2); + + auto initialParameters = Acts::BoundTrackParameters::createCurvilinear( + position4, direction, qOverP, covariance, + Acts::ParticleHypothesis::pion()); + + auto result = fitter.fit( + sourceLinks.begin(), sourceLinks.end(), initialParameters, + fitterOptions, surfaces, tracks); + + if (!result.ok()) { + LOG(debug) << "ACTS Kalman fit failed: " << result.error().message(); + ++nFailed; + continue; + } + + const auto track = result.value(); + if (track.nMeasurements() < 3) { + ++nFailed; + continue; + } + + const auto innerState = track.innermostTrackState(); + if (!innerState.has_value()) { + ++nFailed; + continue; + } + const auto outerState = track.outermostTrackState(); + if (!innerState->hasSmoothed() || !outerState.hasSmoothed() || + !std::isfinite(track.chi2()) || track.chi2() < 0.f) { + ++nFailed; + continue; + } + + // Explicitly use the smoothed fit, rather than the original seed estimate. + typename Trajectory::ConstTrackStateProxy innerProxy{*innerState}; + typename Trajectory::ConstTrackStateProxy outerProxy{outerState}; + candidates.push_back({track.createParametersFromState(innerProxy), + track.createParametersFromState(outerProxy), + seed, static_cast(track.chi2())}); + ++nFitted; + LOG(debug) << "ACTS fitted seed: " + << "nMeasurements=" << track.nMeasurements() + << " chi2=" << track.chi2() + << " nDoF=" << track.nDoF(); + } + + // Accept lower-chi2 candidates first: multiple seeds can share clusters. + std::sort(candidates.begin(), candidates.end(), [](const auto& a, const auto& b) { + return a.chi2 < b.chi2; + }); + int nAccepted = 0; + for (const auto& candidate : candidates) { + nAccepted += acceptTrack(candidate.inner, candidate.outer, candidate.seed, candidate.chi2); + } + LOG(info) << "ACTS fitting: " << nFitted << " fitted, " + << nAccepted << " accepted, " << nFailed << " failed, " + << nInvalid << " invalid from " << mSeeds.size() << " seeds"; +} + +template +bool TrackerACTS::acceptTrack(const Acts::BoundTrackParameters& inner, + const Acts::BoundTrackParameters& outer, + const SeedACTS& seed, float chi2) +{ + if (!mTimeFrame || !std::isfinite(chi2) || chi2 < 0.f) { + return false; + } + + const std::array hits{seed.bottom, seed.middle, seed.top}; + std::array clusterIndices; + clusterIndices.fill(-1); + std::array internalIndices; + internalIndices.fill(-1); + uint32_t pattern = 0; + std::array keys{}; + + for (int i = 0; i < 3; ++i) { + const auto* hit = hits[i]; + if (!hit || hit->layer < 0 || hit->layer >= nLayers || hit->clusterId < 0 || + hit->rof < 0 || hit->rof >= mTimeFrame->getNrof(hit->layer)) { + return false; + } + const uint32_t mask = 1u << hit->layer; + if (pattern & mask) { // the TrackITSExt index array has one slot per layer + return false; + } + pattern |= mask; + keys[i] = (uint64_t(static_cast(hit->layer)) << 32) | + static_cast(hit->clusterId); + if (mAcceptedClusters.contains(keys[i])) { + return false; + } + // Keep the internal index for later MC lookup; publish the external one. + internalIndices[hit->layer] = hit->clusterId; + clusterIndices[hit->layer] = mTimeFrame->getClusterExternalIndex(hit->layer, hit->clusterId); + if (clusterIndices[hit->layer] < 0) { + return false; + } + } + + o2::track::TrackParCov innerO2, outerO2; + const auto& gctx = o2::acts::TrackingGeometryManager::instance().getNominalContext(); + if (!convertTrackParameters(inner, gctx, innerO2) || !convertTrackParameters(outer, gctx, outerO2)) { + return false; + } + + // Intersect the ROF timing intervals exactly as in ITS track acceptance. + const auto& timing = mTimeFrame->getROFOverlapTableView(); + auto nominal = timing.getLayer(hits[0]->layer).getROFTimeBounds(hits[0]->rof); + auto expanded = timing.getLayer(hits[0]->layer).getROFTimeBounds(hits[0]->rof, true); + bool nominalCompatible = true; + for (int i = 1; i < 3; ++i) { + const auto other = timing.getLayer(hits[i]->layer).getROFTimeBounds(hits[i]->rof); + const auto otherExpanded = timing.getLayer(hits[i]->layer).getROFTimeBounds(hits[i]->rof, true); + if (nominalCompatible) { + if (nominal.isCompatible(other)) { + nominal += other; + } else { + nominalCompatible = false; + } + } + if (!expanded.isCompatible(otherExpanded)) { + return false; + } + expanded += otherExpanded; + } + + o2::its::TrackITSExt o2Track{std::move(innerO2), static_cast(hits.size()), chi2, std::move(outerO2), clusterIndices}; + o2Track.setPattern(pattern); + o2Track.getTimeStamp() = (nominalCompatible ? nominal : expanded).makeSymmetrical(); + const float maxHalfROF = timing.getClockLayer().mROFLength * 0.5f; + if (o2Track.getTimeStamp().getTimeStampError() > maxHalfROF) { + o2Track.getTimeStamp().setTimeStampError(maxHalfROF); + } + + // Labels are computed once, after all ROFs, from the internal indices. + mTimeFrame->getTracks().emplace_back(std::move(o2Track)); + mAcceptedTrackClusterIndices.push_back(internalIndices); + for (const auto key : keys) { + mAcceptedClusters.insert(key); + } + return true; +} + +template +void TrackerACTS::computeTracksMClabels() +{ + if (!mTimeFrame || !mTimeFrame->hasMCinformation()) { + return; + } + + auto& tracks = mTimeFrame->getTracks(); + auto& labels = mTimeFrame->getTracksLabel(); + + if (tracks.size() != mAcceptedTrackClusterIndices.size()) { + LOGP(fatal, "ACTS MC labeling: {} tracks but {} internal cluster reference sets", + tracks.size(), mAcceptedTrackClusterIndices.size()); + return; + } + + // There must be exactly one label for each accepted track, in track order. + labels.clear(); + labels.reserve(tracks.size()); + + struct LabelCount { + o2::MCCompLabel label; + int count; + int lastHit; + }; + + for (size_t iTrack = 0; iTrack < tracks.size(); ++iTrack) { + auto& track = tracks[iTrack]; + const auto& indices = mAcceptedTrackClusterIndices[iTrack]; + std::vector counts; + int nHits = 0; + + for (int layer = 0; layer < nLayers; ++layer) { + if (indices[layer] < 0) { + continue; + } + ++nHits; + for (const auto& candidate : mTimeFrame->getClusterLabels(layer, indices[layer])) { + if (!candidate.isValid()) { + continue; + } + + auto it = std::find_if(counts.begin(), counts.end(), [&candidate](const auto& entry) { + return entry.label == candidate; + }); + if (it == counts.end()) { + counts.push_back({candidate, 1, nHits}); + } else if (it->lastHit != nHits) { + // Count a given truth label at most once per attached cluster. + ++it->count; + it->lastHit = nHits; + } + } + } + o2::MCCompLabel winner; + uint32_t pattern = track.getPattern() & 0x0000ffffu; + if (counts.empty()) { + winner.setFakeFlag(); + } else { + const auto best = std::max_element(counts.begin(), counts.end(), [](const auto& a, const auto& b) { + return a.count < b.count; + }); + winner = best->label; + if (best->count != nHits) { + winner.setFakeFlag(); + } + // Preserve the standard TrackITS pattern convention: bits [16, 31] + // indicate clusters compatible with the representative MC label. + for (int layer = 0; layer < nLayers; ++layer) { + if (indices[layer] < 0) { + continue; + } + for (const auto& candidate : mTimeFrame->getClusterLabels(layer, indices[layer])) { + if (candidate.isValid() && candidate == best->label) { + pattern |= (1u << (16 + layer)); + break; + } + } + } + } + track.setPattern(pattern); + labels.emplace_back(winner); + } + LOGP(info, "ACTS MC labeling: {} tracks, {} labels", tracks.size(), labels.size()); +} + +template +void TrackerACTS::clustersToTracks() +{ + if (!mTimeFrame) { + LOG(error) << "Cannot run TrackerACTS: No TimeFrame adopted"; + return; + } + + double totalTime = 0.; + mAcceptedClusters.clear(); + mAcceptedTrackClusterIndices.clear(); + LOG(info) << "==== TRK ACTS Tracking ===="; + + // Verify tracking geometry usage + auto& actsGeometryManager = o2::acts::TrackingGeometryManager::instance(); + const auto& trackingGeometry = actsGeometryManager.get(); + + // Dump tracking geometry to json file for debugging + (void)trackingGeometry; + + auto* trkGeometry = o2::trk::GeometryTGeo::Instance(); + trkGeometry->fillMatrixCache(o2::math_utils::bit2Mask(o2::math_utils::TransformType::L2G)); + const auto& surfaceIndex = actsGeometryManager.getIndex(*trkGeometry, 1.e-3, [trkGeometry](int sensorID) { + return std::string(trkGeometry->getMatrixPath(sensorID).Data()); + }); + + LOG(info) << "ACTS tracking geometry initialized with " << surfaceIndex.size() + << " mapped sensors on " << surfaceIndex.getNSurfaces() << " surfaces"; + try { + dumpGeometryOnce(); + } catch (const std::exception& e) { + LOG(error) << "ACTS geometry debug dump failed: " << e.what(); + } + if (mConfig.writeDebugTrees && mTimeFrameCounter == 0) { + initDebugOutput(); + } + + for (int iROF = 0; iROF < mTimeFrame->getNrof(0); ++iROF) { + LOG(info) << "Processing ROF " << iROF; + + mCurState = SpacePointBuilding; + totalTime += evaluateTask([this, iROF]() { buildSpacePoints(iROF); }, StateNames[mCurState]); + + mCurState = Seeding; + totalTime += evaluateTask([this]() { createSeeds(); }, StateNames[mCurState]); + + const auto nTracksBefore = mTimeFrame->getTracks().size(); + mCurState = TrackFinding; + totalTime += evaluateTask([this]() { findTracks(); }, StateNames[mCurState]); + if (isDebugTimeFrame()) { + mDebugROF.tf = static_cast(mTimeFrameCounter); + mDebugROF.rof = iROF; + mDebugROF.nSpacePoints = static_cast(mSpacePoints.size()); + mDebugROF.nSeeds = static_cast(mSeeds.size()); + mDebugROF.nAcceptedTracks = static_cast(mTimeFrame->getTracks().size() - nTracksBefore); + mROFTree->Fill(); + } + } + + // Run MC labeling once, after every accepted track has been published. + if (mTimeFrame->hasMCinformation()) { + computeTracksMClabels(); + } + + LOG(info) << "=== TimeFrame " << mTimeFrameCounter << " completed in: " << totalTime << " ms ==="; + if (isDebugTimeFrame()) { + mDebugFile->cd(); + if (mHistSpacePoints) { + mHistSpacePoints->Write("", TObject::kOverwrite); + } + mDebugFile->Write("", TObject::kOverwrite); + mDebugFile->Flush(); + } + ++mTimeFrameCounter; + mTotalTime += totalTime; +} + +template +void TrackerACTS::printSummary() const +{ + float avgTF = mTimeFrameCounter > 0 ? static_cast(mTotalTime) / mTimeFrameCounter : 0.f; + LOGP(info, "TrackerACTS summary: Processed {} TFs in TOT={:.2f} ms, AVG/TF={:.2f} ms", mTimeFrameCounter, mTotalTime, avgTF); +} + +// Explicit template instantiations +template class TrackerACTS<11>; +} // namespace o2::trk diff --git a/Detectors/Upgrades/ALICE3/GlobalReconstruction/workflow/CMakeLists.txt b/Detectors/Upgrades/ALICE3/GlobalReconstruction/workflow/CMakeLists.txt index 7b5cd0e735802..8526a4b789251 100644 --- a/Detectors/Upgrades/ALICE3/GlobalReconstruction/workflow/CMakeLists.txt +++ b/Detectors/Upgrades/ALICE3/GlobalReconstruction/workflow/CMakeLists.txt @@ -27,6 +27,10 @@ o2_add_library(ALICE3GlobalReconstructionWorkflow O2::CommonUtils nlohmann_json::nlohmann_json) +if(Acts_FOUND) + target_link_libraries(${targetName} PRIVATE O2::ALICE3ACTS) +endif() + if(CUDA_ENABLED OR HIP_ENABLED) target_compile_definitions(${targetName} PUBLIC TRK_HAS_GPU_TRACKING) endif() diff --git a/Detectors/Upgrades/ALICE3/GlobalReconstruction/workflow/include/ALICE3GlobalReconstructionWorkflow/TrackerSpecImpl.h b/Detectors/Upgrades/ALICE3/GlobalReconstruction/workflow/include/ALICE3GlobalReconstructionWorkflow/TrackerSpecImpl.h index 8a2f162019de0..1f661e637d261 100644 --- a/Detectors/Upgrades/ALICE3/GlobalReconstruction/workflow/include/ALICE3GlobalReconstructionWorkflow/TrackerSpecImpl.h +++ b/Detectors/Upgrades/ALICE3/GlobalReconstruction/workflow/include/ALICE3GlobalReconstructionWorkflow/TrackerSpecImpl.h @@ -22,6 +22,9 @@ #include "Field/MagneticField.h" #include "Framework/ControlService.h" #include "ITStracking/Tracker.h" +#ifdef O2_WITH_ACTS +#include "ALICE3GlobalReconstruction/TrackerACTS.h" +#endif #include "SimulationDataFormat/MCCompLabel.h" #include "SimulationDataFormat/MCEventHeader.h" #include "SimulationDataFormat/MCTruthContainer.h" @@ -114,33 +117,50 @@ void TrackerDPL::runTracking(framework::ProcessingContext& pc, TimeFrameT& timeF nRofs = timeFrame.loadROFrameData(layerROFs, layerClusters, layerPatterns, mIsMC ? &layerLabels : nullptr, yPlaneMLOT); timeFrame.addTruthSeedingVertices(); } - const auto trackingLoopStart = std::chrono::steady_clock::now(); - for (size_t iter{0}; iter < mTrackingParams.size(); ++iter) { - LOGP(info, "{}", mTrackingParams[iter].asString()); - trackerTraits.initialiseTimeFrame(iter); - trackerTraits.computeLayerTracklets(iter, -1); - LOGP(info, "Number of tracklets in iteration {}: {}", iter, timeFrame.getNumberOfTracklets()); - trackerTraits.computeLayerCells(iter); - LOGP(info, "Number of cells in iteration {}: {}", iter, timeFrame.getNumberOfCells()); - trackerTraits.findCellsNeighbours(iter); - LOGP(info, "Number of cell neighbours in iteration {}: {}", iter, timeFrame.getNumberOfNeighbours()); - trackerTraits.findRoads(iter); - LOGP(info, "Number of roads in iteration {}: {}", iter, timeFrame.getNumberOfTracks()); +#ifdef O2_WITH_ACTS + if (mUseACTS) { + LOG(info) << "Running ACTS tracking"; + o2::trk::TrackerACTS<11> actsTracker; + actsTracker.adoptTimeFrame(timeFrame); + actsTracker.adoptSensorIDs(timeFrame.getClusterSensorIDs()); + + const auto& config = mHitRecoConfig.empty() ? mClusterRecoConfig : mHitRecoConfig; + // The input magnetic field is in kG; ACTS TrackerACTS expects Tesla. + actsTracker.setBz(config["geometry"]["bz"].get() * 0.1f); + actsTracker.clustersToTracks(); // This will run the full tracking chain: space point building, seeding, and track finding. + actsTracker.printSummary(); + } else +#endif + { + for (size_t iter{0}; iter < mTrackingParams.size(); ++iter) { + LOGP(info, "{}", mTrackingParams[iter].asString()); + trackerTraits.initialiseTimeFrame(iter); + trackerTraits.computeLayerTracklets(iter, -1); + LOGP(info, "Number of tracklets in iteration {}: {}", iter, timeFrame.getNumberOfTracklets()); + trackerTraits.computeLayerCells(iter); + LOGP(info, "Number of cells in iteration {}: {}", iter, timeFrame.getNumberOfCells()); + trackerTraits.findCellsNeighbours(iter); + LOGP(info, "Number of cell neighbours in iteration {}: {}", iter, timeFrame.getNumberOfNeighbours()); + trackerTraits.findRoads(iter); + LOGP(info, "Number of roads in iteration {}: {}", iter, timeFrame.getNumberOfTracks()); + } + if (mIsMC) { + itsTracker.computeTracksMClabels(); + } } const auto trackingLoopElapsedMs = std::chrono::duration_cast(std::chrono::steady_clock::now() - trackingLoopStart).count(); LOGP(info, "Tracking iterations block took {} ms", trackingLoopElapsedMs); - if (mIsMC) { - itsTracker.computeTracksMClabels(); - } - const auto& tracks = timeFrame.getTracks(); const auto& labels = timeFrame.getTracksLabel(); + if (mIsMC && tracks.size() != labels.size()) { + LOGP(fatal, "TRK: {} reconstructed tracks but {} MC labels", tracks.size(), labels.size()); + } std::vector allTracks(tracks.begin(), tracks.end()); std::vector allLabels; - int totalTracks = allTracks.size(); + const int totalTracks = allTracks.size(); int goodTracks = 0; int fakeTracks = 0; diff --git a/Detectors/Upgrades/ALICE3/GlobalReconstruction/workflow/src/TrackerSpec.cxx b/Detectors/Upgrades/ALICE3/GlobalReconstruction/workflow/src/TrackerSpec.cxx index 33dbad1fecb83..85d62f7b95ce2 100644 --- a/Detectors/Upgrades/ALICE3/GlobalReconstruction/workflow/src/TrackerSpec.cxx +++ b/Detectors/Upgrades/ALICE3/GlobalReconstruction/workflow/src/TrackerSpec.cxx @@ -13,6 +13,9 @@ #include #include #include +#include +#include +#include #include #include #include @@ -43,6 +46,8 @@ #ifdef O2_WITH_ACTS #include "ALICE3GlobalReconstruction/TrackerACTS.h" +#include "ALICE3ACTS/TrkTGeoBlueprintBuilder.h" +#include "ACTSInterface/TrackingGeometryManager.h" #endif #include @@ -75,12 +80,16 @@ TrackerDPL::TrackerDPL(std::shared_ptr gr, int trackingThreads) { if (!hitRecoConfigFileName.empty()) { + LOG(debug) << "Loading ALICE3 TRK hit reconstruction config from file: " << hitRecoConfigFileName; std::ifstream configFile(hitRecoConfigFileName); mHitRecoConfig = nlohmann::json::parse(configFile); + LOG(debug) << "Loaded ALICE3 TRK hit reconstruction config: " << mHitRecoConfig.dump(); } if (!clusterRecoConfigFileName.empty()) { + LOG(debug) << "Loading ALICE3 TRK cluster reconstruction config from file: " << clusterRecoConfigFileName; std::ifstream configFile(clusterRecoConfigFileName); mClusterRecoConfig = nlohmann::json::parse(configFile); + LOG(debug) << "Loaded ALICE3 TRK cluster reconstruction config: " << mClusterRecoConfig.dump(); } mIsMC = isMC; mDeviceType = dType; @@ -91,6 +100,10 @@ void TrackerDPL::init(InitContext& ic) { #ifdef O2_WITH_ACTS mUseACTS = ic.options().get("useACTS"); + if (mUseACTS) { + o2::acts::TrackingGeometryManager::instance().setBuilder(std::make_unique()); + LOG(info) << "Registered automatic ACTS TRK geometry builder"; + } #endif } @@ -104,6 +117,7 @@ std::vector TrackerDPL::createTrackingParamsFromCon std::vector trackingParams; auto loadTrackingParamsFromJson = [](std::vector& trackingParams, const nlohmann::json& paramConfigJson) { for (const auto& paramConfig : paramConfigJson) { + LOG(debug) << "Loading ALICE3 TRK tracking parameters from config: " << paramConfig.dump(); o2::its::TrackingParameters params; auto applyPassFlag = [&](const char* name, o2::its::IterationStep step) { if (paramConfig.contains(name)) { @@ -148,7 +162,6 @@ std::vector TrackerDPL::createTrackingParamsFromCon if (paramConfig.contains("StartLayerMask")) { params.StartLayerMask = paramConfig["StartLayerMask"].get(); } - if (paramConfig.contains("NSigmaCut")) { params.NSigmaCut = paramConfig["NSigmaCut"].get(); } @@ -167,7 +180,6 @@ std::vector TrackerDPL::createTrackingParamsFromCon if (paramConfig.contains("MaxChi2NDF")) { params.MaxChi2NDF = paramConfig["MaxChi2NDF"].get(); } - if (paramConfig.contains("UseDiamond")) { params.UseDiamond = paramConfig["UseDiamond"].get(); } @@ -232,7 +244,6 @@ std::vector TrackerDPL::createTrackingParamsFromCon if (paramConfig.contains("DropTFUponFailure")) { params.DropTFUponFailure = paramConfig["DropTFUponFailure"].get(); } - if (paramConfig.contains("LayerZ")) { params.LayerZ = paramConfig["LayerZ"].get>(); } @@ -319,6 +330,12 @@ void TrackerDPL::run(ProcessingContext& pc) mTimer.Start(false); const bool useGPU = mDeviceType != o2::gpu::gpudatatypes::DeviceType::CPU; +#ifdef O2_WITH_ACTS + if (useGPU && mUseACTS) { + LOG(fatal) << "ACTS tracking is only supported with the CPU backend"; + return; + } +#endif if (useGPU) { runGPUTracking(pc); @@ -430,7 +447,11 @@ DataProcessorSpec getTrackerSpec(bool useMC, const std::string& hitRecoConfig, c clusterRecoConfig, dType, trackingThreads)}, - Options{}}; + Options{ +#ifdef O2_WITH_ACTS + {"useACTS", o2::framework::VariantType::Bool, false, {"Use ACTS for tracking"}} +#endif + }}; } } // namespace trk diff --git a/Detectors/Upgrades/ALICE3/GlobalReconstruction/workflow/src/alice3-global-reconstruction-workflow.cxx b/Detectors/Upgrades/ALICE3/GlobalReconstruction/workflow/src/alice3-global-reconstruction-workflow.cxx index b8df418822ed1..98b6d7a9f5fcf 100644 --- a/Detectors/Upgrades/ALICE3/GlobalReconstruction/workflow/src/alice3-global-reconstruction-workflow.cxx +++ b/Detectors/Upgrades/ALICE3/GlobalReconstruction/workflow/src/alice3-global-reconstruction-workflow.cxx @@ -57,10 +57,6 @@ WorkflowSpec defineDataProcessing(ConfigContext const& configcontext) auto disableRootOutput = configcontext.options().get("disable-root-output"); o2::conf::ConfigurableParam::updateFromString(configcontext.options().get("configKeyValues")); - if (hitRecoConfig.empty() && clusterRecoConfig.empty()) { - throw std::invalid_argument("no reconstruction input configured: provide either --tracking-from-hits-config or --tracking-from-clusters-config "); - } - o2::conf::ConfigurableParam::writeINI("o2alice3globalrecoflow_configuration.ini"); return o2::trk::global_reco_workflow::getWorkflow(useMC, hitRecoConfig, clusterRecoConfig, disableRootOutput, gpuDevice, trackingThreads); From 3afa6e9e972de182f7d04314835c8549c16ae625 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nicol=C3=B2=20Jacazio?= Date: Thu, 5 Nov 2026 06:11:53 +0100 Subject: [PATCH 3/6] Add output --- .../ALICE3GlobalReconstruction/TrackerACTS.h | 54 +++- .../reconstruction/src/TrackerACTS.cxx | 250 +++++++++++++++++- 2 files changed, 290 insertions(+), 14 deletions(-) diff --git a/Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/include/ALICE3GlobalReconstruction/TrackerACTS.h b/Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/include/ALICE3GlobalReconstruction/TrackerACTS.h index 54968b992a2a7..a311e4768f2cd 100644 --- a/Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/include/ALICE3GlobalReconstruction/TrackerACTS.h +++ b/Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/include/ALICE3GlobalReconstruction/TrackerACTS.h @@ -23,6 +23,8 @@ #include #include #include +#include +#include #include #include #include @@ -32,6 +34,8 @@ #include #include +class TFile; +class TTree; namespace Acts { class Surface; @@ -64,6 +68,12 @@ struct TrackerACTSConfig { // Quality cuts float maxSeedsPerMiddleSP = 2; float deltaPhiMax = 0.1f; ///< Maximum phi difference for doublets + // Optional debug output. No ROOT file or histogram is created by default. + bool writeDebugTrees = false; ///< Cluster, seed and ROF summary trees + bool dumpTrackingGeometry = false; ///< Write a JSON volume/surface snapshot on first use + int debugMaxTimeFrames = 1; ///< 1 by default, -1 for all TFs + std::string debugFileName = "TrackerACTS_debug.root"; + std::string geometryDumpFileName = "TrackerACTS_geometry.json"; }; /// Space point representation for tracking @@ -101,14 +111,7 @@ class TrackerACTS { public: TrackerACTS(); - ~TrackerACTS() - { - if (mHistSpacePoints) { - mHistSpacePoints->SaveAs("/tmp/mHistSpacePoints.C"); - delete mHistSpacePoints; - mHistSpacePoints = nullptr; - } - } + ~TrackerACTS(); /// Adopt a TimeFrame for processing void adoptTimeFrame(o2::its::TimeFrame& tf); @@ -135,7 +138,40 @@ class TrackerACTS void printSummary() const; private: - TH2F* mHistSpacePoints = nullptr; + // Debug resources are allocated only if requested and closed at destruction. + std::unique_ptr mDebugFile; + TTree* mClusterTree = nullptr; // Owned by mDebugFile + TTree* mSeedTree = nullptr; // Owned by mDebugFile + TTree* mROFTree = nullptr; // Owned by mDebugFile + std::unique_ptr mHistSpacePoints; + bool mGeometryDumped = false; + + struct DebugClusterRecord { + int tf = 0, rof = 0, layer = 0, clusterId = -1; + int externalId = -1, sensorId = -1, mappingStatus = 0; + unsigned long long geometryId = 0; + float x = 0.f, y = 0.f, z = 0.f; + float r = 0.f, phi = 0.f; + float varianceR = 0.f, varianceZ = 0.f; + } mDebugCluster; + + struct DebugSeedRecord { + int tf = 0, rof = 0, seedId = 0; + float quality = 0.f; + int layer[3] = {}, clusterId[3] = {}, sensorId[3] = {}; + unsigned long long geometryId[3] = {}; + float x[3] = {}, y[3] = {}, z[3] = {}, r[3] = {}, phi[3] = {}; + } mDebugSeed; + + struct DebugROFRecord { + int tf = 0, rof = 0, nSpacePoints = 0, nSeeds = 0; + int nAcceptedTracks = 0; + } mDebugROF; + + void initDebugOutput(); + void dumpGeometryOnce(); + bool isDebugTimeFrame() const; + void fillDebugSeeds(int rof); /// Build space points from clusters in the TimeFrame void buildSpacePoints(int rof); diff --git a/Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/src/TrackerACTS.cxx b/Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/src/TrackerACTS.cxx index 59d0281846433..17744dff61efb 100644 --- a/Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/src/TrackerACTS.cxx +++ b/Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/src/TrackerACTS.cxx @@ -22,6 +22,8 @@ #include "MathUtils/Utils.h" #include "TRKBase/GeometryTGeo.h" +#include +#include #include #include #include @@ -52,10 +54,14 @@ #include "ACTSInterface/TrackingGeometryManager.h" +#include +#include +#include #include #include #include #include +#include #include #include #include @@ -178,9 +184,193 @@ TrackerACTS::TrackerACTS() { // Initialize space points storage per layer mSpacePointsPerLayer.resize(nLayers); - mHistSpacePoints = new TH2F("hSpacePoints", "Space points; x (cm); y (cm)", 200, -100, 100, 200, -100, 100); + // The debug file and histogram are intentionally not allocated here. +} + +template +TrackerACTS::~TrackerACTS() +{ + if (mDebugFile) { + mDebugFile->cd(); + if (mHistSpacePoints) { + mHistSpacePoints->Write("", TObject::kOverwrite); + } + mDebugFile->Write("", TObject::kOverwrite); + mDebugFile->Close(); + } +} + +template +bool TrackerACTS::isDebugTimeFrame() const +{ + return mDebugFile && (mConfig.debugMaxTimeFrames < 0 || + mTimeFrameCounter < static_cast(mConfig.debugMaxTimeFrames)); +} + +template +void TrackerACTS::initDebugOutput() +{ + if (!mConfig.writeDebugTrees || mDebugFile || mConfig.debugMaxTimeFrames == 0) { + return; + } + mDebugFile.reset(TFile::Open(mConfig.debugFileName.c_str(), "RECREATE")); + if (!mDebugFile || mDebugFile->IsZombie()) { + LOG(error) << "Cannot create ACTS debug ROOT file " << mConfig.debugFileName; + mDebugFile.reset(); + return; + } + mDebugFile->cd(); + mClusterTree = new TTree("clusters", "ACTS input clusters and surface mapping; positions in cm; variances in ACTS units"); + mClusterTree->Branch("tf", &mDebugCluster.tf, "tf/I"); + mClusterTree->Branch("rof", &mDebugCluster.rof, "rof/I"); + mClusterTree->Branch("layer", &mDebugCluster.layer, "layer/I"); + mClusterTree->Branch("clusterId", &mDebugCluster.clusterId, "clusterId/I"); + mClusterTree->Branch("externalId", &mDebugCluster.externalId, "externalId/I"); + mClusterTree->Branch("sensorId", &mDebugCluster.sensorId, "sensorId/I"); + mClusterTree->Branch("mappingStatus", &mDebugCluster.mappingStatus, "mappingStatus/I"); + mClusterTree->Branch("geometryId", &mDebugCluster.geometryId, "geometryId/l"); + mClusterTree->Branch("x", &mDebugCluster.x, "x/F"); + mClusterTree->Branch("y", &mDebugCluster.y, "y/F"); + mClusterTree->Branch("z", &mDebugCluster.z, "z/F"); + mClusterTree->Branch("r", &mDebugCluster.r, "r/F"); + mClusterTree->Branch("phi", &mDebugCluster.phi, "phi/F"); + mClusterTree->Branch("varianceR", &mDebugCluster.varianceR, "varianceR/F"); + mClusterTree->Branch("varianceZ", &mDebugCluster.varianceZ, "varianceZ/F"); + + mSeedTree = new TTree("seeds", "ACTS output seed triplets; positions in cm; bottom/middle/top order"); + mSeedTree->Branch("tf", &mDebugSeed.tf, "tf/I"); + mSeedTree->Branch("rof", &mDebugSeed.rof, "rof/I"); + mSeedTree->Branch("seedId", &mDebugSeed.seedId, "seedId/I"); + mSeedTree->Branch("quality", &mDebugSeed.quality, "quality/F"); + mSeedTree->Branch("layer", mDebugSeed.layer, "layer[3]/I"); + mSeedTree->Branch("clusterId", mDebugSeed.clusterId, "clusterId[3]/I"); + mSeedTree->Branch("sensorId", mDebugSeed.sensorId, "sensorId[3]/I"); + mSeedTree->Branch("geometryId", mDebugSeed.geometryId, "geometryId[3]/l"); + mSeedTree->Branch("x", mDebugSeed.x, "x[3]/F"); + mSeedTree->Branch("y", mDebugSeed.y, "y[3]/F"); + mSeedTree->Branch("z", mDebugSeed.z, "z[3]/F"); + mSeedTree->Branch("r", mDebugSeed.r, "r[3]/F"); + mSeedTree->Branch("phi", mDebugSeed.phi, "phi[3]/F"); + + mROFTree = new TTree("rofSummary", "ACTS reconstruction counts per time frame and ROF"); + mROFTree->Branch("tf", &mDebugROF.tf, "tf/I"); + mROFTree->Branch("rof", &mDebugROF.rof, "rof/I"); + mROFTree->Branch("nSpacePoints", &mDebugROF.nSpacePoints, "nSpacePoints/I"); + mROFTree->Branch("nSeeds", &mDebugROF.nSeeds, "nSeeds/I"); + mROFTree->Branch("nAcceptedTracks", &mDebugROF.nAcceptedTracks, "nAcceptedTracks/I"); + + mHistSpacePoints = std::make_unique("hSpacePoints", "Space points; x (cm); y (cm)", + 200, -100, 100, 200, -100, 100); + mHistSpacePoints->SetDirectory(nullptr); // Owned by TrackerACTS, not by TFile. + LOG(info) << "ACTS debug trees enabled: " << mConfig.debugFileName; +} + +template +void TrackerACTS::dumpGeometryOnce() +{ + if (!mConfig.dumpTrackingGeometry || mGeometryDumped) { + return; + } + mGeometryDumped = true; // Even on I/O failure, do not retry on every TF. + auto& manager = o2::acts::TrackingGeometryManager::instance(); + const auto& geometry = manager.get(); + const auto& gctx = manager.getNominalContext(); + auto* trkGeometry = o2::trk::GeometryTGeo::Instance(); + const auto& surfaceIndex = manager.getIndex(*trkGeometry, 1.e-3, [trkGeometry](int sensorID) { + return std::string(trkGeometry->getMatrixPath(sensorID).Data()); + }); - // Initialize ACTS Tracking geometry + std::ofstream out(mConfig.geometryDumpFileName); + if (!out) { + LOG(error) << "Cannot open ACTS geometry dump: " << mConfig.geometryDumpFileName; + return; + } + out << "{\n \"description\": \"ACTS volumes and surfaces snapshot, not a complete ACTS geometry serialization\",\n" + " \"lengthUnits\": \"mm (angles in rad)\",\n \"volumes\": [\n"; + bool first = true; + int nVolumes = 0; + geometry.visitVolumes([&](const Acts::TrackingVolume& volume) { + const auto& transform = volume.localToGlobalTransform(gctx); + nlohmann::json item; + item["geometryId"] = volume.geometryId().value(); + item["motherGeometryId"] = volume.motherVolume() ? volume.motherVolume()->geometryId().value() : 0; + item["name"] = volume.volumeName(); + item["boundsType"] = static_cast(volume.volumeBounds().type()); + item["boundsValues"] = volume.volumeBounds().values(); + item["hasMaterial"] = volume.volumeMaterial() != nullptr; + item["translation"] = {transform.translation().x(), transform.translation().y(), transform.translation().z()}; + item["rotation"] = nlohmann::json::array(); + for (int row = 0; row < 3; ++row) { + item["rotation"].push_back({transform.rotation()(row, 0), transform.rotation()(row, 1), transform.rotation()(row, 2)}); + } + if (!first) { + out << ",\n"; + } + out << item.dump(); + first = false; + ++nVolumes; + }); + out << "\n ],\n \"surfaces\": [\n"; + first = true; + int nSurfaces = 0; + geometry.visitSurfaces([&](const Acts::Surface& surface) { + const auto& transform = surface.localToGlobalTransform(gctx); + nlohmann::json item; + item["geometryId"] = surface.geometryId().value(); + item["surfaceType"] = static_cast(surface.type()); + item["boundsType"] = static_cast(surface.bounds().type()); + item["boundsValues"] = surface.bounds().values(); + item["hasMaterial"] = surface.hasMaterial(); + item["translation"] = {transform.translation().x(), transform.translation().y(), transform.translation().z()}; + item["rotation"] = nlohmann::json::array(); + for (int row = 0; row < 3; ++row) { + item["rotation"].push_back({transform.rotation()(row, 0), transform.rotation()(row, 1), transform.rotation()(row, 2)}); + } + item["sensorIds"] = surfaceIndex.getSensorIDs(surface.geometryId()); + if (!first) { + out << ",\n"; + } + out << item.dump(); + first = false; + ++nSurfaces; + }, false); // Also include passive/non-sensitive surfaces where reachable. + out << "\n ]\n}\n"; + if (!out.good()) { + LOG(error) << "Error while writing ACTS geometry dump " << mConfig.geometryDumpFileName; + return; + } + LOG(info) << "ACTS tracking geometry snapshot: " << nVolumes << " volumes and " + << nSurfaces << " surfaces written to " << mConfig.geometryDumpFileName; +} + +template +void TrackerACTS::fillDebugSeeds(int rof) +{ + if (!isDebugTimeFrame()) { + return; + } + for (size_t iSeed = 0; iSeed < mSeeds.size(); ++iSeed) { + const auto& seed = mSeeds[iSeed]; + mDebugSeed = {}; + mDebugSeed.tf = static_cast(mTimeFrameCounter); + mDebugSeed.rof = rof; + mDebugSeed.seedId = static_cast(iSeed); + mDebugSeed.quality = seed.quality; + const std::array points{seed.bottom, seed.middle, seed.top}; + for (int i = 0; i < 3; ++i) { + const auto& sp = *points[i]; + mDebugSeed.layer[i] = sp.layer; + mDebugSeed.clusterId[i] = sp.clusterId; + mDebugSeed.sensorId[i] = sp.sensorId; + mDebugSeed.geometryId[i] = sp.surface ? sp.surface->geometryId().value() : 0; + mDebugSeed.x[i] = sp.x / Acts::UnitConstants::cm; + mDebugSeed.y[i] = sp.y / Acts::UnitConstants::cm; + mDebugSeed.z[i] = sp.z / Acts::UnitConstants::cm; + mDebugSeed.r[i] = sp.r() / Acts::UnitConstants::cm; + mDebugSeed.phi[i] = sp.phi(); + } + mSeedTree->Fill(); + } } template @@ -224,14 +414,35 @@ void TrackerACTS::buildSpacePoints(int rof) const auto& cluster = clusters[iCluster]; const int clusterIndex = firstCluster + static_cast(iCluster); - if (static_cast(clusterIndex) >= (*mSensorIDs)[layer].size()) { + const bool hasSensorID = clusterIndex >= 0 && static_cast(clusterIndex) < (*mSensorIDs)[layer].size(); + const int sensorID = hasSensorID ? (*mSensorIDs)[layer][clusterIndex] : -1; + const auto* surface = hasSensorID ? surfaceIndex.getSurface(sensorID) : nullptr; + if (isDebugTimeFrame()) { + mDebugCluster = {}; + mDebugCluster.tf = static_cast(mTimeFrameCounter); + mDebugCluster.rof = rof; + mDebugCluster.layer = layer; + mDebugCluster.clusterId = clusterIndex; + mDebugCluster.externalId = mTimeFrame->getClusterExternalIndex(layer, clusterIndex); + mDebugCluster.sensorId = sensorID; + // 0: sensor ID missing; 1: ACTS surface missing; 2: mapped. + mDebugCluster.mappingStatus = !hasSensorID ? 0 : !surface ? 1 : 2; + mDebugCluster.geometryId = surface ? surface->geometryId().value() : 0; + mDebugCluster.x = cluster.xCoordinate; + mDebugCluster.y = cluster.yCoordinate; + mDebugCluster.z = cluster.zCoordinate; + mDebugCluster.r = std::hypot(mDebugCluster.x, mDebugCluster.y); + mDebugCluster.phi = std::atan2(mDebugCluster.y, mDebugCluster.x); + mDebugCluster.varianceR = 0.01f; + mDebugCluster.varianceZ = 0.01f; + mClusterTree->Fill(); + } + if (!hasSensorID) { LOG(error) << "Missing sensor ID for layer " << layer << " cluster " << clusterIndex; ++nMissing; continue; } - const int sensorID = (*mSensorIDs)[layer][clusterIndex]; - const auto* surface = surfaceIndex.getSurface(sensorID); if (!surface) { LOG(error) << "No ACTS surface for sensor " << sensorID << " on layer " << layer; ++nMissing; @@ -406,6 +617,7 @@ void TrackerACTS::createSeeds() output.quality = seed.quality(); mSeeds.push_back(output); } + fillDebugSeeds(mSpacePoints.front().rof); LOGF(info, "Created %zu seeds from %zu space points", mSeeds.size(), mSpacePoints.size()); } @@ -861,8 +1073,11 @@ void TrackerACTS::clustersToTracks() mAcceptedTrackClusterIndices.clear(); LOG(info) << "==== TRK ACTS Tracking ===="; + // Verify tracking geometry usage auto& actsGeometryManager = o2::acts::TrackingGeometryManager::instance(); const auto& trackingGeometry = actsGeometryManager.get(); + + // Dump tracking geometry to json file for debugging (void)trackingGeometry; auto* trkGeometry = o2::trk::GeometryTGeo::Instance(); @@ -873,6 +1088,14 @@ void TrackerACTS::clustersToTracks() LOG(info) << "ACTS tracking geometry initialized with " << surfaceIndex.size() << " mapped sensors on " << surfaceIndex.getNSurfaces() << " surfaces"; + try { + dumpGeometryOnce(); + } catch (const std::exception& e) { + LOG(error) << "ACTS geometry debug dump failed: " << e.what(); + } + if (mConfig.writeDebugTrees && mTimeFrameCounter == 0) { + initDebugOutput(); + } for (int iROF = 0; iROF < mTimeFrame->getNrof(0); ++iROF) { LOG(info) << "Processing ROF " << iROF; @@ -883,8 +1106,17 @@ void TrackerACTS::clustersToTracks() mCurState = Seeding; totalTime += evaluateTask([this]() { createSeeds(); }, StateNames[mCurState]); + const auto nTracksBefore = mTimeFrame->getTracks().size(); mCurState = TrackFinding; totalTime += evaluateTask([this]() { findTracks(); }, StateNames[mCurState]); + if (isDebugTimeFrame()) { + mDebugROF.tf = static_cast(mTimeFrameCounter); + mDebugROF.rof = iROF; + mDebugROF.nSpacePoints = static_cast(mSpacePoints.size()); + mDebugROF.nSeeds = static_cast(mSeeds.size()); + mDebugROF.nAcceptedTracks = static_cast(mTimeFrame->getTracks().size() - nTracksBefore); + mROFTree->Fill(); + } } // Run MC labeling once, after every accepted track has been published. @@ -893,6 +1125,14 @@ void TrackerACTS::clustersToTracks() } LOG(info) << "=== TimeFrame " << mTimeFrameCounter << " completed in: " << totalTime << " ms ==="; + if (isDebugTimeFrame()) { + mDebugFile->cd(); + if (mHistSpacePoints) { + mHistSpacePoints->Write("", TObject::kOverwrite); + } + mDebugFile->Write("", TObject::kOverwrite); + mDebugFile->Flush(); + } ++mTimeFrameCounter; mTotalTime += totalTime; } From b03c7448e9d08bec04429f0e73786e9f623fb514 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nicol=C3=B2=20Jacazio?= Date: Thu, 5 Nov 2026 06:33:35 +0100 Subject: [PATCH 4/6] Update --- .../ALICE3GlobalReconstruction/TrackerACTS.h | 79 ++++++---- .../reconstruction/src/TrackerACTS.cxx | 138 ++++++++++++++---- 2 files changed, 165 insertions(+), 52 deletions(-) diff --git a/Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/include/ALICE3GlobalReconstruction/TrackerACTS.h b/Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/include/ALICE3GlobalReconstruction/TrackerACTS.h index a311e4768f2cd..5cc25c53f0f94 100644 --- a/Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/include/ALICE3GlobalReconstruction/TrackerACTS.h +++ b/Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/include/ALICE3GlobalReconstruction/TrackerACTS.h @@ -44,68 +44,90 @@ class BoundTrackParameters; namespace o2::trk { -/// Configuration for the ACTS-based tracker +/// Configuration for ACTS seeding, track finding, and optional diagnostic output. +/// +/// Length and momentum-valued cuts use ACTS units (mm and GeV, respectively). struct TrackerACTSConfig { - // Seeding parameters + /// Minimum transverse momentum accepted by the seeding algorithms. float minPt = 0.4f * Acts::UnitConstants::GeV; ///< Minimum pT for seeds - float maxImpactParameter = 10.f * Acts::UnitConstants::mm; ///< Maximum impact parameter - float cotThetaMax = std::sinh(4.0f); ///< Maximum cot(theta), corresponds to eta ~4 + /// Maximum transverse distance of closest approach used in seeding. + float maxImpactParameter = 10.f * Acts::UnitConstants::mm; + /// Maximum absolute cot(theta); the default corresponds approximately to |eta| = 4. + float cotThetaMax = std::sinh(4.0f); - // Delta R cuts for doublet/triplet formation + /// Minimum and maximum radial separation for bottom-to-middle space-point pairs. float deltaRMinBottom = 5.f * Acts::UnitConstants::mm; ///< Min deltaR for bottom-middle float deltaRMaxBottom = 200.f * Acts::UnitConstants::mm; ///< Max deltaR for bottom-middle + /// Minimum and maximum radial separation for middle-to-top space-point pairs. float deltaRMinTop = 5.f * Acts::UnitConstants::mm; ///< Min deltaR for middle-top float deltaRMaxTop = 200.f * Acts::UnitConstants::mm; ///< Max deltaR for middle-top - // Z cuts + /// Longitudinal range used to select space points for seeding. float zMin = -3000.f * Acts::UnitConstants::mm; float zMax = 3000.f * Acts::UnitConstants::mm; - // Collision region + /// Longitudinal range of the collision region used when constructing seeds. float collisionRegionMin = -150.f * Acts::UnitConstants::mm; float collisionRegionMax = 150.f * Acts::UnitConstants::mm; - // Quality cuts + /// Maximum number of seeds retained for each middle space point. float maxSeedsPerMiddleSP = 2; - float deltaPhiMax = 0.1f; ///< Maximum phi difference for doublets - // Optional debug output. No ROOT file or histogram is created by default. - bool writeDebugTrees = false; ///< Cluster, seed and ROF summary trees - bool dumpTrackingGeometry = false; ///< Write a JSON volume/surface snapshot on first use - int debugMaxTimeFrames = 1; ///< 1 by default, -1 for all TFs + /// Reserved maximum azimuthal angle difference for doublets (radians); currently unused. + float deltaPhiMax = 0.1f; + + /// Enable diagnostic ROOT trees and histograms; disabled by default. + bool writeDebugTrees = false; + /// Write a JSON snapshot of the tracking geometry on the first use. + bool dumpTrackingGeometry = false; + /// Number of time frames with debug output; 0 disables it and -1 means all. + int debugMaxTimeFrames = 1; + /// Destination ROOT file for diagnostic output. std::string debugFileName = "TrackerACTS_debug.root"; + /// Destination JSON file for the tracking-geometry snapshot. std::string geometryDumpFileName = "TrackerACTS_geometry.json"; }; -/// Space point representation for tracking +/// Cluster-derived space point passed to the ACTS seeding algorithms. struct SpacePoint { + /// Cartesian position in the global tracking frame, in mm. float x{0.f}; float y{0.f}; float z{0.f}; + /// Detector layer and cluster indices identifying the source measurement. int layer{-1}; int clusterId{-1}; + /// Readout-frame and sensor indices identifying the source measurement. int rof{-1}; int sensorId{-1}; + /// Geometry surface associated with this measurement; not owned by the space point. const Acts::Surface* surface{nullptr}; - // Derived quantities + /// Transverse radius in mm. float r() const { return std::hypot(x, y); } - float radius() const { return r(); } // required by Acts::CylindricalGridElement concept + /// Radius accessor required by the ACTS cylindrical-grid element concept. + float radius() const { return r(); } + /// Azimuthal angle in radians. float phi() const { return std::atan2(y, x); } - // Variance estimates (can be refined based on cluster properties) - float varianceR{0.01f}; // ~100 um resolution squared + /// Estimated radial and longitudinal position variances, in mm^2. + float varianceR{0.01f}; float varianceZ{0.01f}; }; -/// Seed (triplet of space points) +/// Three-space-point seed ordered from the innermost to outermost measurement. struct SeedACTS { const SpacePoint* bottom{nullptr}; const SpacePoint* middle{nullptr}; const SpacePoint* top{nullptr}; + /// ACTS seed quality value. float quality{0.f}; }; -/// TRK Tracker using ACTS algorithms for seeding and track finding +/// TRK tracker using ACTS algorithms for seeding and track finding. +/// +/// The caller supplies a TimeFrame and, optionally, sensor-ID mappings before +/// calling clustersToTracks(). The adopted data is referenced rather than +/// copied, so it must remain alive while the tracker processes it. template class TrackerACTS { @@ -113,28 +135,33 @@ class TrackerACTS TrackerACTS(); ~TrackerACTS(); - /// Adopt a TimeFrame for processing + /// Set the TimeFrame whose clusters will be reconstructed. + /// The TimeFrame is borrowed and must outlive its use by this tracker. void adoptTimeFrame(o2::its::TimeFrame& tf); + /// Set per-layer sensor-ID mappings used to associate clusters with geometry. + /// The mapping array is borrowed and must remain alive while it is used. void adoptSensorIDs(const std::array, nLayers>& sensorIDs) { mSensorIDs = &sensorIDs; } - /// Main tracking entry point: convert clusters to tracks + /// Build space points, find tracks, and publish accepted tracks in the TimeFrame. void clustersToTracks(); - /// Configuration + /// Replace the current tracker configuration. void setConfig(const TrackerACTSConfig& cfg) { mConfig = cfg; } + /// Access the tracker configuration for modification. TrackerACTSConfig& getConfig() { return mConfig; } + /// Read-only access to the tracker configuration. const TrackerACTSConfig& getConfig() const { return mConfig; } - /// Set the magnetic field strength + /// Set the magnetic field along z, in Tesla. void setBz(float bz) { mBz = bz; } - /// Get the magnetic field strength + /// Get the magnetic field along z, in Tesla. float getBz() const { return mBz; } - /// Print tracking summary + /// Print the accumulated tracking timing and processing summary. void printSummary() const; private: diff --git a/Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/src/TrackerACTS.cxx b/Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/src/TrackerACTS.cxx index 17744dff61efb..2d9751af5511f 100644 --- a/Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/src/TrackerACTS.cxx +++ b/Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/src/TrackerACTS.cxx @@ -15,6 +15,26 @@ /// \author Nicolò Jacazio, Università del Piemonte Orientale (IT) /// \since 2026-04-01 /// +/// \details +/// The tracker processes one readout frame (ROF) at a time. For each ROF it +/// maps clusters to ACTS surfaces, creates ACTS space-point triplet seeds, +/// fits each seed with the ACTS Kalman fitter, and publishes accepted tracks +/// in the O2 TimeFrame. Optional ROOT and JSON output records intermediate +/// information for diagnostics; neither output is required for reconstruction. +/// +/// ACTS geometry is supplied by the process-wide +/// o2::acts::TrackingGeometryManager. Before its first build, the workflow +/// must make the aligned O2 TGeo geometry available and install the +/// detector-specific ACTS geometry builder. The manager builds the geometry +/// lazily and keeps it immutable; surface pointers stored in space points are +/// non-owning references to that manager-owned geometry. +/// +/// Cluster coordinates from the TimeFrame are in cm and are converted to +/// ACTS mm when space points are created. ACTS contexts, local coordinates, +/// field values, and fit parameters must therefore use ACTS units. Conversion +/// back to the O2 track representation is performed only when a track is +/// accepted. +/// #include "ALICE3GlobalReconstruction/TrackerACTS.h" @@ -24,6 +44,7 @@ #include #include +#include #include #include #include @@ -75,10 +96,14 @@ namespace o2::trk namespace { -// Convert ACTS surface-bound parameters to the O2 tracking frame -// (x, alpha; y, z, sin(phi-alpha), tan(lambda), q/pT). -// The covariance is transformed by a numerical Jacobian evaluated at the -// fitted state. Position units are converted from ACTS mm to O2 cm. +/// Convert an ACTS bound state into O2's local track parameterization. +/// +/// O2 stores (x, alpha; y, z, sin(phi-alpha), tan(lambda), q/pT), whereas +/// ACTS stores surface-local coordinates and bound direction/momentum +/// parameters. The covariance is propagated through the same conversion with +/// a central-difference Jacobian. Positions are converted from ACTS mm to O2 +/// cm and q/p is converted to q/pT. ACTS has a time coordinate but +/// TrackParCov does not, so the time row/column is not included in the output. bool convertTrackParameters(const Acts::BoundTrackParameters& actsTrack, const Acts::GeometryContext& gctx, o2::track::TrackParCov& output) @@ -203,6 +228,7 @@ TrackerACTS::~TrackerACTS() template bool TrackerACTS::isDebugTimeFrame() const { + // A negative limit means unlimited debug output; zero disables it. return mDebugFile && (mConfig.debugMaxTimeFrames < 0 || mTimeFrameCounter < static_cast(mConfig.debugMaxTimeFrames)); } @@ -220,6 +246,8 @@ void TrackerACTS::initDebugOutput() return; } mDebugFile->cd(); + // The tree branches point into these member records, which are refilled for + // each entry. The file owns the trees; the tracker owns the histogram. mClusterTree = new TTree("clusters", "ACTS input clusters and surface mapping; positions in cm; variances in ACTS units"); mClusterTree->Branch("tf", &mDebugCluster.tf, "tf/I"); mClusterTree->Branch("rof", &mDebugCluster.rof, "rof/I"); @@ -273,9 +301,15 @@ void TrackerACTS::dumpGeometryOnce() } mGeometryDumped = true; // Even on I/O failure, do not retry on every TF. auto& manager = o2::acts::TrackingGeometryManager::instance(); + // get() builds the shared geometry on first access if necessary. Its + // nominal context is appropriate because O2 alignment was already applied + // to the TGeo matrices from which this ACTS geometry was built. const auto& geometry = manager.get(); const auto& gctx = manager.getNominalContext(); auto* trkGeometry = o2::trk::GeometryTGeo::Instance(); + // Build/reuse the sensor-to-surface map. The matrix cache is populated in + // clustersToTracks(); the path provider makes matching use the TGeo sensor + // path as an identity instead of relying only on transform proximity. const auto& surfaceIndex = manager.getIndex(*trkGeometry, 1.e-3, [trkGeometry](int sensorID) { return std::string(trkGeometry->getMatrixPath(sensorID).Data()); }); @@ -285,19 +319,25 @@ void TrackerACTS::dumpGeometryOnce() LOG(error) << "Cannot open ACTS geometry dump: " << mConfig.geometryDumpFileName; return; } - out << "{\n \"description\": \"ACTS volumes and surfaces snapshot, not a complete ACTS geometry serialization\",\n" + out << "{\n \"description\": \"ACTS volumes and sensitive surfaces snapshot, not a complete ACTS geometry serialization\",\n" " \"lengthUnits\": \"mm (angles in rad)\",\n \"volumes\": [\n"; bool first = true; int nVolumes = 0; - geometry.visitVolumes([&](const Acts::TrackingVolume& volume) { - const auto& transform = volume.localToGlobalTransform(gctx); + // ACTS visitor callbacks receive non-owning pointers. The snapshot records + // global transforms in the nominal context and basic bounds/material + // metadata; it is intended for inspection, not for rebuilding the geometry. + geometry.visitVolumes([&](const Acts::TrackingVolume* volume) { + if (volume == nullptr) { + return; + } + const auto& transform = volume->localToGlobalTransform(gctx); nlohmann::json item; - item["geometryId"] = volume.geometryId().value(); - item["motherGeometryId"] = volume.motherVolume() ? volume.motherVolume()->geometryId().value() : 0; - item["name"] = volume.volumeName(); - item["boundsType"] = static_cast(volume.volumeBounds().type()); - item["boundsValues"] = volume.volumeBounds().values(); - item["hasMaterial"] = volume.volumeMaterial() != nullptr; + item["geometryId"] = volume->geometryId().value(); + item["motherGeometryId"] = volume->motherVolume() ? volume->motherVolume()->geometryId().value() : 0; + item["name"] = volume->volumeName(); + item["boundsType"] = static_cast(volume->volumeBounds().type()); + item["boundsValues"] = volume->volumeBounds().values(); + item["hasMaterial"] = volume->volumeMaterial() != nullptr; item["translation"] = {transform.translation().x(), transform.translation().y(), transform.translation().z()}; item["rotation"] = nlohmann::json::array(); for (int row = 0; row < 3; ++row) { @@ -313,27 +353,33 @@ void TrackerACTS::dumpGeometryOnce() out << "\n ],\n \"surfaces\": [\n"; first = true; int nSurfaces = 0; - geometry.visitSurfaces([&](const Acts::Surface& surface) { - const auto& transform = surface.localToGlobalTransform(gctx); + // The one-argument TrackingGeometry visitor visits sensitive surfaces. + // Associate each geometry identifier with the O2 sensor IDs known to the + // detector-specific index; a surface can have no matching sensor ID. + geometry.visitSurfaces([&](const Acts::Surface* surface) { + if (surface == nullptr) { + return; + } + const auto& transform = surface->localToGlobalTransform(gctx); nlohmann::json item; - item["geometryId"] = surface.geometryId().value(); - item["surfaceType"] = static_cast(surface.type()); - item["boundsType"] = static_cast(surface.bounds().type()); - item["boundsValues"] = surface.bounds().values(); - item["hasMaterial"] = surface.hasMaterial(); + item["geometryId"] = surface->geometryId().value(); + item["surfaceType"] = static_cast(surface->type()); + item["boundsType"] = static_cast(surface->bounds().type()); + item["boundsValues"] = surface->bounds().values(); + item["hasMaterial"] = surface->hasMaterial(); item["translation"] = {transform.translation().x(), transform.translation().y(), transform.translation().z()}; item["rotation"] = nlohmann::json::array(); for (int row = 0; row < 3; ++row) { item["rotation"].push_back({transform.rotation()(row, 0), transform.rotation()(row, 1), transform.rotation()(row, 2)}); } - item["sensorIds"] = surfaceIndex.getSensorIDs(surface.geometryId()); + item["sensorIds"] = surfaceIndex.getSensorIDs(surface->geometryId()); if (!first) { out << ",\n"; } out << item.dump(); first = false; ++nSurfaces; - }, false); // Also include passive/non-sensitive surfaces where reachable. + }); out << "\n ]\n}\n"; if (!out.good()) { LOG(error) << "Error while writing ACTS geometry dump " << mConfig.geometryDumpFileName; @@ -357,6 +403,8 @@ void TrackerACTS::fillDebugSeeds(int rof) mDebugSeed.seedId = static_cast(iSeed); mDebugSeed.quality = seed.quality; const std::array points{seed.bottom, seed.middle, seed.top}; + // Keep the output order consistent with the seeder and seed-triplet + // convention: innermost (bottom), middle, then outermost (top). for (int i = 0; i < 3; ++i) { const auto& sp = *points[i]; mDebugSeed.layer[i] = sp.layer; @@ -376,6 +424,7 @@ void TrackerACTS::fillDebugSeeds(int rof) template void TrackerACTS::adoptTimeFrame(o2::its::TimeFrame& tf) { + // Borrow the TimeFrame; its lifetime must cover all calls that use it. mTimeFrame = &tf; } @@ -395,8 +444,13 @@ void TrackerACTS::buildSpacePoints(int rof) auto& actsGeometryManager = o2::acts::TrackingGeometryManager::instance(); auto* trkGeometry = o2::trk::GeometryTGeo::Instance(); + // Surface matching needs global-to-local transforms from O2's geometry + // helper. Build its L2G cache before requesting the manager's detector index. trkGeometry->fillMatrixCache(o2::math_utils::bit2Mask(o2::math_utils::TransformType::L2G)); + // mSensorIDs is parallel to each layer's cluster ordering. Its sensor ID + // selects the corresponding ACTS surface; missing IDs or surfaces cannot be + // used by seeding/fitting and are skipped below. const auto& surfaceIndex = actsGeometryManager.getIndex(*trkGeometry, 1.e-3, [trkGeometry](int sensorID) { return std::string(trkGeometry->getMatrixPath(sensorID).Data()); }); @@ -450,6 +504,8 @@ void TrackerACTS::buildSpacePoints(int rof) } SpacePoint sp; + // SpacePoint and all subsequent ACTS calculations use mm, while the + // TimeFrame cluster coordinates are stored in cm. sp.x = cluster.xCoordinate * Acts::UnitConstants::cm; sp.y = cluster.yCoordinate * Acts::UnitConstants::cm; sp.z = cluster.zCoordinate * Acts::UnitConstants::cm; @@ -491,6 +547,9 @@ void TrackerACTS::createSeeds() const float bFieldInZ = mBz * Acts::UnitConstants::T; + // ACTS seeding consumes cylindrical (phi, z, r) bins in ACTS units. The + // bins are populated from the owned space-point vector, whose indices are + // retained so ACTS seeds can be translated back to SpacePoint pointers. Acts::CylindricalSpacePointGrid::Config gridConfig; gridConfig.minPt = mConfig.minPt; gridConfig.rMin = 0.f; @@ -630,12 +689,15 @@ void TrackerACTS::findTracks() } auto& actsGeometryManager = o2::acts::TrackingGeometryManager::instance(); + // The same nominal context used to construct/map the geometry must be used + // for surface transforms and propagation. const auto& geometryContext = actsGeometryManager.getNominalContext(); Acts::MagneticFieldContext magneticFieldContext; Acts::CalibrationContext calibrationContext; - // Initial implementation with a constant magnetic field + // Initial implementation with a constant magnetic field; the ACTS geometry + // is used for surface navigation, while this field is uniform in z. auto magneticField = std::make_shared( Acts::Vector3{0., 0., mBz * Acts::UnitConstants::T}); @@ -674,6 +736,8 @@ void TrackerACTS::findTracks() { const auto& hit = *sourceLink.get().hit; + // The hit was converted to this ACTS surface's local coordinates before + // fitting. Calibrate only the two measured local position coordinates. trackState.setUncalibratedSourceLink(Acts::SourceLink{sourceLink}); trackState.allocateCalibrated(2); trackState.calibrated<2>() = hit.local; @@ -789,6 +853,8 @@ void TrackerACTS::findTracks() const auto* surface = sp->surface; Acts::Vector3 position{sp->x, sp->y, sp->z}; + // Convert the cluster's global position to the coordinates expected by + // a measurement on this ACTS surface. auto localResult = surface->globalToLocal( geometryContext, position, direction, 1. * Acts::UnitConstants::mm); @@ -803,6 +869,8 @@ void TrackerACTS::findTracks() measurement.local = localResult.value(); const double sigma = 100. * Acts::UnitConstants::um; + // This initial fitter uses a common 100 um diagonal resolution for all + // three measurements; detector-specific covariances are not yet passed. measurement.covariance.setZero(); measurement.covariance(0, 0) = sigma * sigma; measurement.covariance(1, 1) = sigma * sigma; @@ -906,6 +974,9 @@ bool TrackerACTS::acceptTrack(const Acts::BoundTrackParameters& inner, uint32_t pattern = 0; std::array keys{}; + // TrackITSExt has one cluster slot per detector layer. Retain the internal + // TimeFrame cluster index for MC lookup, but publish its external index in + // the O2 track. for (int i = 0; i < 3; ++i) { const auto* hit = hits[i]; if (!hit || hit->layer < 0 || hit->layer >= nLayers || hit->clusterId < 0 || @@ -931,12 +1002,17 @@ bool TrackerACTS::acceptTrack(const Acts::BoundTrackParameters& inner, } o2::track::TrackParCov innerO2, outerO2; + // Convert the smoothed ACTS parameters at both ends of the fitted + // trajectory. Conversion uses the nominal geometry context so that bound + // local parameters are interpreted on their associated surfaces. const auto& gctx = o2::acts::TrackingGeometryManager::instance().getNominalContext(); if (!convertTrackParameters(inner, gctx, innerO2) || !convertTrackParameters(outer, gctx, outerO2)) { return false; } - // Intersect the ROF timing intervals exactly as in ITS track acceptance. + // Require the three source clusters to have compatible readout timing. Use + // the nominal overlap when possible, otherwise the widened overlap, exactly + // as in ITS track acceptance. const auto& timing = mTimeFrame->getROFOverlapTableView(); auto nominal = timing.getLayer(hits[0]->layer).getROFTimeBounds(hits[0]->rof); auto expanded = timing.getLayer(hits[0]->layer).getROFTimeBounds(hits[0]->rof, true); @@ -1028,6 +1104,8 @@ void TrackerACTS::computeTracksMClabels() } } } + // Choose the truth label occurring in the largest number of attached + // clusters. Mark it fake unless it is represented on every hit. o2::MCCompLabel winner; uint32_t pattern = track.getPattern() & 0x0000ffffu; if (counts.empty()) { @@ -1073,14 +1151,20 @@ void TrackerACTS::clustersToTracks() mAcceptedTrackClusterIndices.clear(); LOG(info) << "==== TRK ACTS Tracking ===="; - // Verify tracking geometry usage + // Trigger (or retrieve) the process-wide immutable ACTS geometry. A builder + // must already have been installed and the aligned O2 geometry must be + // available before this first access. auto& actsGeometryManager = o2::acts::TrackingGeometryManager::instance(); const auto& trackingGeometry = actsGeometryManager.get(); - // Dump tracking geometry to json file for debugging + // Keep the geometry alive and make the explicit initialization step clear; + // detailed lookups are done through the detector-specific surface index. (void)trackingGeometry; auto* trkGeometry = o2::trk::GeometryTGeo::Instance(); + // The index uses this L2G cache and an exact TGeo path provider to map each + // TRK sensor ID to the ACTS surface built for that sensor. It is cached by + // TrackingGeometryManager and reused by the ROF processing below. trkGeometry->fillMatrixCache(o2::math_utils::bit2Mask(o2::math_utils::TransformType::L2G)); const auto& surfaceIndex = actsGeometryManager.getIndex(*trkGeometry, 1.e-3, [trkGeometry](int sensorID) { return std::string(trkGeometry->getMatrixPath(sensorID).Data()); @@ -1097,6 +1181,8 @@ void TrackerACTS::clustersToTracks() initDebugOutput(); } + // Each ROF is processed independently for space-point creation, seeding, + // and fitting. Accepted tracks accumulate in the TimeFrame across ROFs. for (int iROF = 0; iROF < mTimeFrame->getNrof(0); ++iROF) { LOG(info) << "Processing ROF " << iROF; From 369db8b2dd3e3f267553f33a48f1c49e05db9325 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nicol=C3=B2=20Jacazio?= Date: Thu, 5 Nov 2026 06:35:00 +0100 Subject: [PATCH 5/6] U --- .../ALICE3GlobalReconstruction/TrackerACTS.h | 8 +- .../ALICE3GlobalReconstruction/TrackerACTS2.h | 236 ---- .../reconstruction/src/TrackerACTS2.cxx | 1149 ----------------- 3 files changed, 4 insertions(+), 1389 deletions(-) delete mode 100644 Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/include/ALICE3GlobalReconstruction/TrackerACTS2.h delete mode 100644 Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/src/TrackerACTS2.cxx diff --git a/Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/include/ALICE3GlobalReconstruction/TrackerACTS.h b/Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/include/ALICE3GlobalReconstruction/TrackerACTS.h index 5cc25c53f0f94..bfc83ef269536 100644 --- a/Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/include/ALICE3GlobalReconstruction/TrackerACTS.h +++ b/Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/include/ALICE3GlobalReconstruction/TrackerACTS.h @@ -40,7 +40,7 @@ namespace Acts { class Surface; class BoundTrackParameters; -} +} // namespace Acts namespace o2::trk { @@ -49,7 +49,7 @@ namespace o2::trk /// Length and momentum-valued cuts use ACTS units (mm and GeV, respectively). struct TrackerACTSConfig { /// Minimum transverse momentum accepted by the seeding algorithms. - float minPt = 0.4f * Acts::UnitConstants::GeV; ///< Minimum pT for seeds + float minPt = 0.4f * Acts::UnitConstants::GeV; ///< Minimum pT for seeds /// Maximum transverse distance of closest approach used in seeding. float maxImpactParameter = 10.f * Acts::UnitConstants::mm; /// Maximum absolute cot(theta); the default corresponds approximately to |eta| = 4. @@ -59,8 +59,8 @@ struct TrackerACTSConfig { float deltaRMinBottom = 5.f * Acts::UnitConstants::mm; ///< Min deltaR for bottom-middle float deltaRMaxBottom = 200.f * Acts::UnitConstants::mm; ///< Max deltaR for bottom-middle /// Minimum and maximum radial separation for middle-to-top space-point pairs. - float deltaRMinTop = 5.f * Acts::UnitConstants::mm; ///< Min deltaR for middle-top - float deltaRMaxTop = 200.f * Acts::UnitConstants::mm; ///< Max deltaR for middle-top + float deltaRMinTop = 5.f * Acts::UnitConstants::mm; ///< Min deltaR for middle-top + float deltaRMaxTop = 200.f * Acts::UnitConstants::mm; ///< Max deltaR for middle-top /// Longitudinal range used to select space points for seeding. float zMin = -3000.f * Acts::UnitConstants::mm; diff --git a/Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/include/ALICE3GlobalReconstruction/TrackerACTS2.h b/Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/include/ALICE3GlobalReconstruction/TrackerACTS2.h deleted file mode 100644 index f480f403048d2..0000000000000 --- a/Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/include/ALICE3GlobalReconstruction/TrackerACTS2.h +++ /dev/null @@ -1,236 +0,0 @@ -// 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. - -/// -/// \file TrackerACTS.h -/// \brief TRK tracker using ACTS seeding algorithm -/// \author Nicolò Jacazio, Università del Piemonte Orientale (IT) -/// \since 2026-04-01 -/// - -#ifndef ALICE3_GLOBALRECONSTRUCTION_INCLUDE_TRACKERACTS_H_ -#define ALICE3_GLOBALRECONSTRUCTION_INCLUDE_TRACKERACTS_H_ - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include "DataFormatsITS/TrackITS.h" -#include -#include - -class TFile; -class TTree; -namespace Acts -{ -class Surface; -class BoundTrackParameters; -} -namespace o2::trk -{ -/// Configuration for the ACTS-based tracker -struct TrackerACTSConfig { - // Seeding parameters - float minPt = 0.4f * Acts::UnitConstants::GeV; ///< Minimum pT for seeds - float maxImpactParameter = 10.f * Acts::UnitConstants::mm; ///< Maximum impact parameter - float cotThetaMax = std::sinh(4.0f); ///< Maximum cot(theta), corresponds to eta ~4 - // Delta R cuts for doublet/triplet formation - float deltaRMinBottom = 5.f * Acts::UnitConstants::mm; ///< Min deltaR for bottom-middle - float deltaRMaxBottom = 200.f * Acts::UnitConstants::mm; ///< Max deltaR for bottom-middle - float deltaRMinTop = 5.f * Acts::UnitConstants::mm; ///< Min deltaR for middle-top - float deltaRMaxTop = 200.f * Acts::UnitConstants::mm; ///< Max deltaR for middle-top - // Z cuts - float zMin = -3000.f * Acts::UnitConstants::mm; - float zMax = 3000.f * Acts::UnitConstants::mm; - - // Collision region - float collisionRegionMin = -150.f * Acts::UnitConstants::mm; - float collisionRegionMax = 150.f * Acts::UnitConstants::mm; - - // Quality cuts - float maxSeedsPerMiddleSP = 2; - float deltaPhiMax = 0.1f; ///< Maximum phi difference for doublets - // Optional debug output. No ROOT file or histogram is created by default. - bool writeDebugTrees = false; ///< Cluster, seed and ROF summary trees - bool dumpTrackingGeometry = false; ///< Write a JSON volume/surface snapshot on first use - int debugMaxTimeFrames = 1; ///< 1 by default, -1 for all TFs - std::string debugFileName = "TrackerACTS_debug.root"; - std::string geometryDumpFileName = "TrackerACTS_geometry.json"; -}; -/// Space point representation for tracking -struct SpacePoint { - float x{0.f}; - float y{0.f}; - float z{0.f}; - int layer{-1}; - int clusterId{-1}; - int rof{-1}; - int sensorId{-1}; - const Acts::Surface* surface{nullptr}; - - // Derived quantities - float r() const { return std::hypot(x, y); } - float radius() const { return r(); } // required by Acts::CylindricalGridElement concept - float phi() const { return std::atan2(y, x); } - // Variance estimates (can be refined based on cluster properties) - float varianceR{0.01f}; // ~100 um resolution squared - float varianceZ{0.01f}; -}; - -/// Seed (triplet of space points) -struct SeedACTS { - const SpacePoint* bottom{nullptr}; - const SpacePoint* middle{nullptr}; - const SpacePoint* top{nullptr}; - float quality{0.f}; -}; -/// TRK Tracker using ACTS algorithms for seeding and track finding -template -class TrackerACTS -{ - public: - TrackerACTS(); - ~TrackerACTS(); - - /// Adopt a TimeFrame for processing - void adoptTimeFrame(o2::its::TimeFrame& tf); - void adoptSensorIDs(const std::array, nLayers>& sensorIDs) - { - mSensorIDs = &sensorIDs; - } - - /// Main tracking entry point: convert clusters to tracks - void clustersToTracks(); - - /// Configuration - void setConfig(const TrackerACTSConfig& cfg) { mConfig = cfg; } - TrackerACTSConfig& getConfig() { return mConfig; } - const TrackerACTSConfig& getConfig() const { return mConfig; } - /// Configure the magnetic field strength (T). - void setBz(float bz) { mBz = bz; } - float getBz() const { return mBz; } - - /// Print tracking summary - void printSummary() const; - - private: - // Debug resources are allocated only if requested and closed at destruction. - std::unique_ptr mDebugFile; - TTree* mClusterTree = nullptr; // Owned by mDebugFile - TTree* mSeedTree = nullptr; // Owned by mDebugFile - TTree* mROFTree = nullptr; // Owned by mDebugFile - std::unique_ptr mHistSpacePoints; - bool mGeometryDumped = false; - - struct DebugClusterRecord { - int tf = 0, rof = 0, layer = 0, clusterId = -1; - int externalId = -1, sensorId = -1, mappingStatus = 0; - unsigned long long geometryId = 0; - float x = 0.f, y = 0.f, z = 0.f; - float r = 0.f, phi = 0.f; - float varianceR = 0.f, varianceZ = 0.f; - } mDebugCluster; - - struct DebugSeedRecord { - int tf = 0, rof = 0, seedId = 0; - float quality = 0.f; - int layer[3] = {}, clusterId[3] = {}, sensorId[3] = {}; - unsigned long long geometryId[3] = {}; - float x[3] = {}, y[3] = {}, z[3] = {}, r[3] = {}, phi[3] = {}; - } mDebugSeed; - - struct DebugROFRecord { - int tf = 0, rof = 0, nSpacePoints = 0, nSeeds = 0; - int nAcceptedTracks = 0; - } mDebugROF; - - void initDebugOutput(); - void dumpGeometryOnce(); - bool isDebugTimeFrame() const; - void fillDebugSeeds(int rof); - - /// Build space points from clusters in the TimeFrame - void buildSpacePoints(int rof); - - /// Create seeds (triplets) from space points using ACTS SeedFinder - void createSeeds(); - - /// Run track finding from seeds - void findTracks(); - - /// Convert fitted ACTS states to a TrackITSExt and publish it in the TimeFrame. - bool acceptTrack(const Acts::BoundTrackParameters& inner, - const Acts::BoundTrackParameters& outer, - const SeedACTS& seed, float chi2); - /// Assign one MC label per accepted track, after all ROFs have been processed. - void computeTracksMClabels(); - - /// Helper: time a task - template - float evaluateTask(Func&& task, std::string_view taskName); - - // Configuration - TrackerACTSConfig mConfig; - - // TimeFrame data - o2::its::TimeFrame* mTimeFrame = nullptr; - - // Space points built from clusters - std::vector mSpacePoints; - std::vector> mSpacePointsPerLayer; - // Seeds - std::vector mSeeds; - // Avoid accepting seeds which reuse a cluster already assigned to a track. - std::unordered_set mAcceptedClusters; - // Layer-local indices in mUnsortedClusters, in the same order as getTracks(). - // Keep them separate from the external indices stored in TrackITSExt. - std::vector> mAcceptedTrackClusterIndices; - - // Tracking state - float mBz{0.5f}; ///< Magnetic field in Tesla - unsigned int mTimeFrameCounter{0}; - double mTotalTime{0.}; - const std::array, nLayers>* mSensorIDs = nullptr; - // Tracking states for logging - enum State { - SpacePointBuilding = 0, - Seeding, - TrackFinding, - NStates, - }; - State mCurState{SpacePointBuilding}; - static constexpr std::array StateNames{ - "Space point building", - "ACTS seeding", - "Track finding"}; -}; -template -template -float TrackerACTS::evaluateTask(Func&& task, std::string_view taskName) -{ - LOG(debug) << " + Starting " << taskName; - const auto start = std::chrono::high_resolution_clock::now(); - task(); - const auto end = std::chrono::high_resolution_clock::now(); - std::chrono::duration diff{end - start}; - LOG(debug) << " - " << taskName << " completed in: " << std::fixed << std::setprecision(2) << diff.count() << " ms"; - return static_cast(diff.count()); -} - -} // namespace o2::trk -#endif /* ALICE3_GLOBALRECONSTRUCTION_INCLUDE_TRACKERACTS_H_ */ diff --git a/Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/src/TrackerACTS2.cxx b/Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/src/TrackerACTS2.cxx deleted file mode 100644 index 17744dff61efb..0000000000000 --- a/Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/src/TrackerACTS2.cxx +++ /dev/null @@ -1,1149 +0,0 @@ -// 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. - -/// -/// \file TrackerACTS.cxx -/// \brief TRK tracker using ACTS seeding and track finding -/// \author Nicolò Jacazio, Università del Piemonte Orientale (IT) -/// \since 2026-04-01 -/// - -#include "ALICE3GlobalReconstruction/TrackerACTS.h" - -#include "ACTSInterface/TrackingGeometryManager.h" -#include "MathUtils/Utils.h" -#include "TRKBase/GeometryTGeo.h" - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include "ACTSInterface/TrackingGeometryManager.h" - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -namespace o2::trk -{ - -namespace -{ -// Convert ACTS surface-bound parameters to the O2 tracking frame -// (x, alpha; y, z, sin(phi-alpha), tan(lambda), q/pT). -// The covariance is transformed by a numerical Jacobian evaluated at the -// fitted state. Position units are converted from ACTS mm to O2 cm. -bool convertTrackParameters(const Acts::BoundTrackParameters& actsTrack, - const Acts::GeometryContext& gctx, - o2::track::TrackParCov& output) -{ - const auto& state = actsTrack.parameters(); - if (!state.allFinite() || !actsTrack.covariance().has_value() || - !actsTrack.covariance()->allFinite()) { - return false; - } - - const auto position = actsTrack.position(gctx); - if (!position.allFinite()) { - return false; - } - const double alpha = std::atan2(position.y(), position.x()); - const double ca = std::cos(alpha); - const double sa = std::sin(alpha); - - auto toO2 = [&](const Acts::BoundVector& p, Eigen::Matrix& dest) -> bool { - const double theta = p[Acts::eBoundTheta]; - const double sinTheta = std::sin(theta); - if (!std::isfinite(sinTheta) || sinTheta <= 1.e-6) { - return false; - } - const auto direction = Acts::makeDirectionFromPhiTheta(p[Acts::eBoundPhi], theta); - const Acts::Vector2 local = p.head<2>(); - const auto xyz = actsTrack.referenceSurface().localToGlobal(gctx, local, direction); - if (!xyz.allFinite()) { - return false; - } - - const double snp = std::sin(p[Acts::eBoundPhi] - alpha); - if (!std::isfinite(snp) || std::abs(snp) >= 0.999999) { - return false; - } - - dest << (-xyz.x() * sa + xyz.y() * ca) / Acts::UnitConstants::cm, - xyz.z() / Acts::UnitConstants::cm, - snp, - std::cos(theta) / sinTheta, - p[Acts::eBoundQOverP] * Acts::UnitConstants::GeV / sinTheta; - return dest.allFinite(); - }; - - Eigen::Matrix pars; - if (!toO2(state, pars)) { - return false; - } - - const double x = (position.x() * ca + position.y() * sa) / Acts::UnitConstants::cm; - if (!std::isfinite(x)) { - return false; - } - - Eigen::Matrix jacobian; - jacobian.setZero(); - for (int i = 0; i < Acts::eBoundSize; ++i) { - if (i == Acts::eBoundTime) { - continue; // O2 TrackParCov has no time coordinate - } - const double step = (i == Acts::eBoundLoc0 || i == Acts::eBoundLoc1) - ? 1.e-3 * Acts::UnitConstants::mm - : (i == Acts::eBoundQOverP) - ? std::max(std::abs(state[i]) * 1.e-4, 1.e-6 / Acts::UnitConstants::GeV) - : 1.e-5; - Acts::BoundVector plus = state; - Acts::BoundVector minus = state; - plus[i] += step; - minus[i] -= step; - Eigen::Matrix fPlus, fMinus; - if (!toO2(plus, fPlus) || !toO2(minus, fMinus)) { - return false; - } - jacobian.col(i) = (fPlus - fMinus) / (2. * step); - } - - const Eigen::Matrix covariance = - jacobian * (*actsTrack.covariance()) * jacobian.transpose(); - if (!covariance.allFinite()) { - return false; - } - - o2::track::TrackParCov::params_t o2Params{}; - o2::track::TrackParCov::covMat_t o2Cov{}; - for (int i = 0; i < 5; ++i) { - if (!(covariance(i, i) > 0.) || !std::isfinite(pars[i])) { - return false; - } - o2Params[i] = static_cast(pars[i]); - for (int j = 0; j <= i; ++j) { - o2Cov[o2::track::CovarMap[i][j]] = - static_cast(0.5 * (covariance(i, j) + covariance(j, i))); - } - } - const int charge = state[Acts::eBoundQOverP] >= 0. ? 1 : -1; - output = o2::track::TrackParCov{static_cast(x), static_cast(alpha), - o2Params, o2Cov, charge, o2::track::PID::Pion}; - return output.isValid(); -} -} // namespace -template -TrackerACTS::TrackerACTS() -{ - // Initialize space points storage per layer - mSpacePointsPerLayer.resize(nLayers); - // The debug file and histogram are intentionally not allocated here. -} - -template -TrackerACTS::~TrackerACTS() -{ - if (mDebugFile) { - mDebugFile->cd(); - if (mHistSpacePoints) { - mHistSpacePoints->Write("", TObject::kOverwrite); - } - mDebugFile->Write("", TObject::kOverwrite); - mDebugFile->Close(); - } -} - -template -bool TrackerACTS::isDebugTimeFrame() const -{ - return mDebugFile && (mConfig.debugMaxTimeFrames < 0 || - mTimeFrameCounter < static_cast(mConfig.debugMaxTimeFrames)); -} - -template -void TrackerACTS::initDebugOutput() -{ - if (!mConfig.writeDebugTrees || mDebugFile || mConfig.debugMaxTimeFrames == 0) { - return; - } - mDebugFile.reset(TFile::Open(mConfig.debugFileName.c_str(), "RECREATE")); - if (!mDebugFile || mDebugFile->IsZombie()) { - LOG(error) << "Cannot create ACTS debug ROOT file " << mConfig.debugFileName; - mDebugFile.reset(); - return; - } - mDebugFile->cd(); - mClusterTree = new TTree("clusters", "ACTS input clusters and surface mapping; positions in cm; variances in ACTS units"); - mClusterTree->Branch("tf", &mDebugCluster.tf, "tf/I"); - mClusterTree->Branch("rof", &mDebugCluster.rof, "rof/I"); - mClusterTree->Branch("layer", &mDebugCluster.layer, "layer/I"); - mClusterTree->Branch("clusterId", &mDebugCluster.clusterId, "clusterId/I"); - mClusterTree->Branch("externalId", &mDebugCluster.externalId, "externalId/I"); - mClusterTree->Branch("sensorId", &mDebugCluster.sensorId, "sensorId/I"); - mClusterTree->Branch("mappingStatus", &mDebugCluster.mappingStatus, "mappingStatus/I"); - mClusterTree->Branch("geometryId", &mDebugCluster.geometryId, "geometryId/l"); - mClusterTree->Branch("x", &mDebugCluster.x, "x/F"); - mClusterTree->Branch("y", &mDebugCluster.y, "y/F"); - mClusterTree->Branch("z", &mDebugCluster.z, "z/F"); - mClusterTree->Branch("r", &mDebugCluster.r, "r/F"); - mClusterTree->Branch("phi", &mDebugCluster.phi, "phi/F"); - mClusterTree->Branch("varianceR", &mDebugCluster.varianceR, "varianceR/F"); - mClusterTree->Branch("varianceZ", &mDebugCluster.varianceZ, "varianceZ/F"); - - mSeedTree = new TTree("seeds", "ACTS output seed triplets; positions in cm; bottom/middle/top order"); - mSeedTree->Branch("tf", &mDebugSeed.tf, "tf/I"); - mSeedTree->Branch("rof", &mDebugSeed.rof, "rof/I"); - mSeedTree->Branch("seedId", &mDebugSeed.seedId, "seedId/I"); - mSeedTree->Branch("quality", &mDebugSeed.quality, "quality/F"); - mSeedTree->Branch("layer", mDebugSeed.layer, "layer[3]/I"); - mSeedTree->Branch("clusterId", mDebugSeed.clusterId, "clusterId[3]/I"); - mSeedTree->Branch("sensorId", mDebugSeed.sensorId, "sensorId[3]/I"); - mSeedTree->Branch("geometryId", mDebugSeed.geometryId, "geometryId[3]/l"); - mSeedTree->Branch("x", mDebugSeed.x, "x[3]/F"); - mSeedTree->Branch("y", mDebugSeed.y, "y[3]/F"); - mSeedTree->Branch("z", mDebugSeed.z, "z[3]/F"); - mSeedTree->Branch("r", mDebugSeed.r, "r[3]/F"); - mSeedTree->Branch("phi", mDebugSeed.phi, "phi[3]/F"); - - mROFTree = new TTree("rofSummary", "ACTS reconstruction counts per time frame and ROF"); - mROFTree->Branch("tf", &mDebugROF.tf, "tf/I"); - mROFTree->Branch("rof", &mDebugROF.rof, "rof/I"); - mROFTree->Branch("nSpacePoints", &mDebugROF.nSpacePoints, "nSpacePoints/I"); - mROFTree->Branch("nSeeds", &mDebugROF.nSeeds, "nSeeds/I"); - mROFTree->Branch("nAcceptedTracks", &mDebugROF.nAcceptedTracks, "nAcceptedTracks/I"); - - mHistSpacePoints = std::make_unique("hSpacePoints", "Space points; x (cm); y (cm)", - 200, -100, 100, 200, -100, 100); - mHistSpacePoints->SetDirectory(nullptr); // Owned by TrackerACTS, not by TFile. - LOG(info) << "ACTS debug trees enabled: " << mConfig.debugFileName; -} - -template -void TrackerACTS::dumpGeometryOnce() -{ - if (!mConfig.dumpTrackingGeometry || mGeometryDumped) { - return; - } - mGeometryDumped = true; // Even on I/O failure, do not retry on every TF. - auto& manager = o2::acts::TrackingGeometryManager::instance(); - const auto& geometry = manager.get(); - const auto& gctx = manager.getNominalContext(); - auto* trkGeometry = o2::trk::GeometryTGeo::Instance(); - const auto& surfaceIndex = manager.getIndex(*trkGeometry, 1.e-3, [trkGeometry](int sensorID) { - return std::string(trkGeometry->getMatrixPath(sensorID).Data()); - }); - - std::ofstream out(mConfig.geometryDumpFileName); - if (!out) { - LOG(error) << "Cannot open ACTS geometry dump: " << mConfig.geometryDumpFileName; - return; - } - out << "{\n \"description\": \"ACTS volumes and surfaces snapshot, not a complete ACTS geometry serialization\",\n" - " \"lengthUnits\": \"mm (angles in rad)\",\n \"volumes\": [\n"; - bool first = true; - int nVolumes = 0; - geometry.visitVolumes([&](const Acts::TrackingVolume& volume) { - const auto& transform = volume.localToGlobalTransform(gctx); - nlohmann::json item; - item["geometryId"] = volume.geometryId().value(); - item["motherGeometryId"] = volume.motherVolume() ? volume.motherVolume()->geometryId().value() : 0; - item["name"] = volume.volumeName(); - item["boundsType"] = static_cast(volume.volumeBounds().type()); - item["boundsValues"] = volume.volumeBounds().values(); - item["hasMaterial"] = volume.volumeMaterial() != nullptr; - item["translation"] = {transform.translation().x(), transform.translation().y(), transform.translation().z()}; - item["rotation"] = nlohmann::json::array(); - for (int row = 0; row < 3; ++row) { - item["rotation"].push_back({transform.rotation()(row, 0), transform.rotation()(row, 1), transform.rotation()(row, 2)}); - } - if (!first) { - out << ",\n"; - } - out << item.dump(); - first = false; - ++nVolumes; - }); - out << "\n ],\n \"surfaces\": [\n"; - first = true; - int nSurfaces = 0; - geometry.visitSurfaces([&](const Acts::Surface& surface) { - const auto& transform = surface.localToGlobalTransform(gctx); - nlohmann::json item; - item["geometryId"] = surface.geometryId().value(); - item["surfaceType"] = static_cast(surface.type()); - item["boundsType"] = static_cast(surface.bounds().type()); - item["boundsValues"] = surface.bounds().values(); - item["hasMaterial"] = surface.hasMaterial(); - item["translation"] = {transform.translation().x(), transform.translation().y(), transform.translation().z()}; - item["rotation"] = nlohmann::json::array(); - for (int row = 0; row < 3; ++row) { - item["rotation"].push_back({transform.rotation()(row, 0), transform.rotation()(row, 1), transform.rotation()(row, 2)}); - } - item["sensorIds"] = surfaceIndex.getSensorIDs(surface.geometryId()); - if (!first) { - out << ",\n"; - } - out << item.dump(); - first = false; - ++nSurfaces; - }, false); // Also include passive/non-sensitive surfaces where reachable. - out << "\n ]\n}\n"; - if (!out.good()) { - LOG(error) << "Error while writing ACTS geometry dump " << mConfig.geometryDumpFileName; - return; - } - LOG(info) << "ACTS tracking geometry snapshot: " << nVolumes << " volumes and " - << nSurfaces << " surfaces written to " << mConfig.geometryDumpFileName; -} - -template -void TrackerACTS::fillDebugSeeds(int rof) -{ - if (!isDebugTimeFrame()) { - return; - } - for (size_t iSeed = 0; iSeed < mSeeds.size(); ++iSeed) { - const auto& seed = mSeeds[iSeed]; - mDebugSeed = {}; - mDebugSeed.tf = static_cast(mTimeFrameCounter); - mDebugSeed.rof = rof; - mDebugSeed.seedId = static_cast(iSeed); - mDebugSeed.quality = seed.quality; - const std::array points{seed.bottom, seed.middle, seed.top}; - for (int i = 0; i < 3; ++i) { - const auto& sp = *points[i]; - mDebugSeed.layer[i] = sp.layer; - mDebugSeed.clusterId[i] = sp.clusterId; - mDebugSeed.sensorId[i] = sp.sensorId; - mDebugSeed.geometryId[i] = sp.surface ? sp.surface->geometryId().value() : 0; - mDebugSeed.x[i] = sp.x / Acts::UnitConstants::cm; - mDebugSeed.y[i] = sp.y / Acts::UnitConstants::cm; - mDebugSeed.z[i] = sp.z / Acts::UnitConstants::cm; - mDebugSeed.r[i] = sp.r() / Acts::UnitConstants::cm; - mDebugSeed.phi[i] = sp.phi(); - } - mSeedTree->Fill(); - } -} - -template -void TrackerACTS::adoptTimeFrame(o2::its::TimeFrame& tf) -{ - mTimeFrame = &tf; -} - -template -void TrackerACTS::buildSpacePoints(int rof) -{ - mSpacePoints.clear(); - for (auto& layerSPs : mSpacePointsPerLayer) { - layerSPs.clear(); - } - - if (!mSensorIDs) { - LOG(error) << "Cannot build ACTS space points: Sensor IDs not adopted"; - return; - } - - auto& actsGeometryManager = o2::acts::TrackingGeometryManager::instance(); - auto* trkGeometry = o2::trk::GeometryTGeo::Instance(); - - trkGeometry->fillMatrixCache(o2::math_utils::bit2Mask(o2::math_utils::TransformType::L2G)); - - const auto& surfaceIndex = actsGeometryManager.getIndex(*trkGeometry, 1.e-3, [trkGeometry](int sensorID) { - return std::string(trkGeometry->getMatrixPath(sensorID).Data()); - }); - - int nMapped = 0; - int nMissing = 0; - - for (int layer = 0; layer < nLayers; ++layer) { - const auto clusters = mTimeFrame->getUnsortedClustersOnLayer(rof, layer); - const int firstCluster = mTimeFrame->getSortedStartIndex(rof, layer); - - LOG(debug) << "ACTSTracker: got " << clusters.size() << " clusters on layer " << layer; - - for (size_t iCluster = 0; iCluster < clusters.size(); ++iCluster) { - const auto& cluster = clusters[iCluster]; - const int clusterIndex = firstCluster + static_cast(iCluster); - - const bool hasSensorID = clusterIndex >= 0 && static_cast(clusterIndex) < (*mSensorIDs)[layer].size(); - const int sensorID = hasSensorID ? (*mSensorIDs)[layer][clusterIndex] : -1; - const auto* surface = hasSensorID ? surfaceIndex.getSurface(sensorID) : nullptr; - if (isDebugTimeFrame()) { - mDebugCluster = {}; - mDebugCluster.tf = static_cast(mTimeFrameCounter); - mDebugCluster.rof = rof; - mDebugCluster.layer = layer; - mDebugCluster.clusterId = clusterIndex; - mDebugCluster.externalId = mTimeFrame->getClusterExternalIndex(layer, clusterIndex); - mDebugCluster.sensorId = sensorID; - // 0: sensor ID missing; 1: ACTS surface missing; 2: mapped. - mDebugCluster.mappingStatus = !hasSensorID ? 0 : !surface ? 1 : 2; - mDebugCluster.geometryId = surface ? surface->geometryId().value() : 0; - mDebugCluster.x = cluster.xCoordinate; - mDebugCluster.y = cluster.yCoordinate; - mDebugCluster.z = cluster.zCoordinate; - mDebugCluster.r = std::hypot(mDebugCluster.x, mDebugCluster.y); - mDebugCluster.phi = std::atan2(mDebugCluster.y, mDebugCluster.x); - mDebugCluster.varianceR = 0.01f; - mDebugCluster.varianceZ = 0.01f; - mClusterTree->Fill(); - } - if (!hasSensorID) { - LOG(error) << "Missing sensor ID for layer " << layer << " cluster " << clusterIndex; - ++nMissing; - continue; - } - - if (!surface) { - LOG(error) << "No ACTS surface for sensor " << sensorID << " on layer " << layer; - ++nMissing; - continue; - } - - SpacePoint sp; - sp.x = cluster.xCoordinate * Acts::UnitConstants::cm; - sp.y = cluster.yCoordinate * Acts::UnitConstants::cm; - sp.z = cluster.zCoordinate * Acts::UnitConstants::cm; - - if (mHistSpacePoints) { - mHistSpacePoints->Fill(sp.x / Acts::UnitConstants::cm, sp.y / Acts::UnitConstants::cm); - } - - sp.layer = layer; - sp.clusterId = clusterIndex; - sp.rof = rof; - sp.sensorId = sensorID; - sp.surface = surface; - sp.varianceR = 0.01f; - sp.varianceZ = 0.01f; - - mSpacePoints.push_back(sp); - ++nMapped; - } - } - - for (auto& sp : mSpacePoints) { - if (sp.layer >= 0 && sp.layer < nLayers) { - mSpacePointsPerLayer[sp.layer].push_back(&sp); - } - } - - LOG(info) << "ACTS surface mapping: " << nMapped << " clusters mapped, " << nMissing << " clusters missing"; -} - -template -void TrackerACTS::createSeeds() -{ - mSeeds.clear(); - if (mSpacePoints.empty()) { - LOG(info) << "No space points available for seeding"; - return; - } - - const float bFieldInZ = mBz * Acts::UnitConstants::T; - - Acts::CylindricalSpacePointGrid::Config gridConfig; - gridConfig.minPt = mConfig.minPt; - gridConfig.rMin = 0.f; - gridConfig.rMax = 100.f * Acts::UnitConstants::cm; - gridConfig.zMin = mConfig.zMin; - gridConfig.zMax = mConfig.zMax; - gridConfig.deltaRMax = std::max(mConfig.deltaRMaxBottom, mConfig.deltaRMaxTop); - gridConfig.cotThetaMax = mConfig.cotThetaMax; - gridConfig.impactMax = mConfig.maxImpactParameter; - gridConfig.bFieldInZ = bFieldInZ; - gridConfig.bottomBinFinder = Acts::GridBinFinder<3ul>{1, std::vector>{}, 0}; - gridConfig.topBinFinder = Acts::GridBinFinder<3ul>{1, std::vector>{}, 0}; - - Acts::CylindricalSpacePointGrid grid{gridConfig}; - for (Acts::SpacePointIndex i = 0; i < mSpacePoints.size(); ++i) { - const auto& sp = mSpacePoints[i]; - grid.insert(i, sp.phi(), sp.z, sp.r()); - } - - // The doublet finder expects each bin to be sorted in increasing radius. - for (size_t i = 0; i < grid.numberOfBins(); ++i) { - auto& bin = grid.at(i); - std::sort(bin.begin(), bin.end(), [this](Acts::SpacePointIndex a, Acts::SpacePointIndex b) { - return mSpacePoints[a].r() < mSpacePoints[b].r(); - }); - } - - // ACTS 47 uses a non-templated, column-based space-point container. - Acts::SpacePointContainer spacePoints{ - Acts::SpacePointColumns::CopiedFromIndex | - Acts::SpacePointColumns::PackedXY | - Acts::SpacePointColumns::PackedZR | - Acts::SpacePointColumns::VarianceR | - Acts::SpacePointColumns::VarianceZ}; - spacePoints.reserve(grid.numberOfSpacePoints()); - - std::vector gridRanges; - gridRanges.reserve(grid.numberOfBins()); - for (size_t i = 0; i < grid.numberOfBins(); ++i) { - const auto first = static_cast(spacePoints.size()); - for (Acts::SpacePointIndex spIndex : grid.at(i)) { - const auto& sp = mSpacePoints[spIndex]; - auto proxy = spacePoints.createSpacePoint(); - proxy.copiedFromIndex() = spIndex; - proxy.xy() = std::array{sp.x, sp.y}; - proxy.zr() = std::array{sp.z, sp.r()}; - proxy.varianceR() = sp.varianceR; - proxy.varianceZ() = sp.varianceZ; - } - const auto last = static_cast(spacePoints.size()); - gridRanges.emplace_back(first, last); - } - - Acts::DoubletSeedFinder::Config bottomConfig; - bottomConfig.spacePointsSortedByRadius = true; - bottomConfig.candidateDirection = Acts::Direction::Backward(); - bottomConfig.deltaRMin = mConfig.deltaRMinBottom; - bottomConfig.deltaRMax = mConfig.deltaRMaxBottom; - bottomConfig.collisionRegionMin = mConfig.collisionRegionMin; - bottomConfig.collisionRegionMax = mConfig.collisionRegionMax; - bottomConfig.cotThetaMax = mConfig.cotThetaMax; - bottomConfig.impactMax = mConfig.maxImpactParameter; - bottomConfig.minPt = mConfig.minPt; - auto bottomFinder = Acts::DoubletSeedFinder::create( - Acts::DoubletSeedFinder::DerivedConfig{bottomConfig, bFieldInZ}); - - Acts::DoubletSeedFinder::Config topConfig = bottomConfig; - topConfig.candidateDirection = Acts::Direction::Forward(); - topConfig.deltaRMin = mConfig.deltaRMinTop; - topConfig.deltaRMax = mConfig.deltaRMaxTop; - auto topFinder = Acts::DoubletSeedFinder::create( - Acts::DoubletSeedFinder::DerivedConfig{topConfig, bFieldInZ}); - - Acts::TripletSeedFinder::Config tripletConfig; - tripletConfig.minPt = mConfig.minPt; - tripletConfig.impactMax = mConfig.maxImpactParameter; - tripletConfig.sigmaScattering = 5.f; - tripletConfig.radLengthPerSeed = 0.05f; - auto tripletFinder = Acts::TripletSeedFinder::create( - Acts::TripletSeedFinder::DerivedConfig{tripletConfig, bFieldInZ}); - - Acts::BroadTripletSeedFilter::Config filterConfig; - filterConfig.maxSeedsPerSpM = static_cast(mConfig.maxSeedsPerMiddleSP); - filterConfig.deltaRMin = std::min(mConfig.deltaRMinBottom, mConfig.deltaRMinTop); - Acts::BroadTripletSeedFilter::State filterState; - Acts::BroadTripletSeedFilter::Cache filterCache; - Acts::BroadTripletSeedFilter filter{filterConfig, filterState, filterCache, Acts::getDummyLogger()}; - - Acts::TripletSeeder seeder; - Acts::TripletSeeder::Cache seedingCache; - Acts::SeedContainer seeds; - seeds.assignSpacePointContainer(spacePoints); - - std::vector bottomRanges; - std::vector topRanges; - for (const auto [bottom, middle, top] : grid.binnedGroup()) { - bottomRanges.clear(); - for (const auto bin : bottom) { - bottomRanges.push_back(spacePoints.range(gridRanges.at(bin)).asConst()); - } - topRanges.clear(); - for (const auto bin : top) { - topRanges.push_back(spacePoints.range(gridRanges.at(bin)).asConst()); - } - auto middleRange = spacePoints.range(gridRanges.at(middle)).asConst(); - if (middleRange.empty()) { - continue; - } - - seeder.createSeedsFromGroups(seedingCache, *bottomFinder, *topFinder, - *tripletFinder, filter, spacePoints, - bottomRanges, middleRange, topRanges, - {gridConfig.rMin, gridConfig.rMax}, seeds); - } - - // Restore the existing SeedACTS representation using original SP indices. - mSeeds.reserve(seeds.size()); - for (const auto seed : seeds) { - const auto indices = seed.spacePointIndices(); - SeedACTS output; - output.bottom = &mSpacePoints[spacePoints.at(indices[0]).copiedFromIndex()]; - output.middle = &mSpacePoints[spacePoints.at(indices[1]).copiedFromIndex()]; - output.top = &mSpacePoints[spacePoints.at(indices[2]).copiedFromIndex()]; - output.quality = seed.quality(); - mSeeds.push_back(output); - } - fillDebugSeeds(mSpacePoints.front().rof); - - LOGF(info, "Created %zu seeds from %zu space points", mSeeds.size(), mSpacePoints.size()); -} - -template -void TrackerACTS::findTracks() -{ - if (mSeeds.empty()) { - return; - } - - auto& actsGeometryManager = o2::acts::TrackingGeometryManager::instance(); - const auto& geometryContext = actsGeometryManager.getNominalContext(); - - Acts::MagneticFieldContext magneticFieldContext; - Acts::CalibrationContext calibrationContext; - - // Initial implementation with a constant magnetic field - auto magneticField = std::make_shared( - Acts::Vector3{0., 0., mBz * Acts::UnitConstants::T}); - - using Propagator = Acts::Propagator, Acts::DirectNavigator>; - using Trajectory = Acts::VectorMultiTrajectory; - using TrackContainer = Acts::TrackContainer; - - Propagator propagator{ - Acts::EigenStepper<>{magneticField}, - Acts::DirectNavigator{}}; - - Acts::KalmanFitter fitter{std::move(propagator)}; - - struct HitMeasurement { - const Acts::Surface* surface{nullptr}; - Acts::Vector2 local = Acts::Vector2::Zero(); - Acts::SquareMatrix<2> covariance = Acts::SquareMatrix<2>::Zero(); - }; - - struct HitSourceLink { - const HitMeasurement* hit{nullptr}; - }; - - struct SurfaceAccessor { - const Acts::Surface* getSurface(const Acts::SourceLink& sourceLink) const - { - return sourceLink.get().hit->surface; - } - }; - - struct MeasurementCalibrator { - void calibrate(const Acts::GeometryContext&, - const Acts::CalibrationContext&, - const Acts::SourceLink& sourceLink, - Trajectory::TrackStateProxy trackState) const - { - const auto& hit = *sourceLink.get().hit; - - trackState.setUncalibratedSourceLink(Acts::SourceLink{sourceLink}); - trackState.allocateCalibrated(2); - trackState.calibrated<2>() = hit.local; - trackState.calibratedCovariance<2>() = hit.covariance; - - trackState.setProjectorSubspaceIndices(std::array{static_cast(Acts::eBoundLoc0), - static_cast(Acts::eBoundLoc1)}); - } - }; - - SurfaceAccessor surfaceAccessor; - MeasurementCalibrator calibrator; - Acts::GainMatrixUpdater updater; - Acts::GainMatrixSmoother smoother; - - Acts::KalmanFitterExtensions extensions; - extensions.surfaceAccessor.connect<&SurfaceAccessor::getSurface>(&surfaceAccessor); - extensions.calibrator.connect<&MeasurementCalibrator::calibrate>(&calibrator); - extensions.updater.connect<&Acts::GainMatrixUpdater::operator()>(&updater); - extensions.smoother.connect<&Acts::GainMatrixSmoother::operator()>(&smoother); - - Acts::PropagatorPlainOptions propagatorOptions{geometryContext, magneticFieldContext}; - - Acts::KalmanFitterOptions fitterOptions{geometryContext, - magneticFieldContext, - calibrationContext, - extensions, - propagatorOptions, - nullptr, - false, - false}; - - // Mutable TrackContainer uses reference-holder backends, not rvalue shared_ptrs. - Acts::VectorTrackContainer trackBackend; - Trajectory trajectoryBackend; - TrackContainer tracks{trackBackend, trajectoryBackend}; - - struct FittedCandidate { - Acts::BoundTrackParameters inner; - Acts::BoundTrackParameters outer; - SeedACTS seed; - float chi2; - }; - std::vector candidates; - int nFitted = 0; - int nFailed = 0; - int nInvalid = 0; - - for (const auto& seed : mSeeds) { - const std::array spacePoints{ - seed.bottom, seed.middle, seed.top}; - - if (!spacePoints[0] || !spacePoints[1] || !spacePoints[2]) { - ++nInvalid; - continue; - } - - if (!spacePoints[0]->surface || !spacePoints[1]->surface || - !spacePoints[2]->surface) { - ++nInvalid; - continue; - } - - // The direct fitter expects at most one measurement per surface. - if (spacePoints[0]->surface == spacePoints[1]->surface || - spacePoints[0]->surface == spacePoints[2]->surface || - spacePoints[1]->surface == spacePoints[2]->surface) { - ++nInvalid; - continue; - } - - Acts::Vector3 pos0{spacePoints[0]->x, spacePoints[0]->y, spacePoints[0]->z}; - Acts::Vector3 pos1{spacePoints[1]->x, spacePoints[1]->y, spacePoints[1]->z}; - Acts::Vector3 pos2{spacePoints[2]->x, spacePoints[2]->y, spacePoints[2]->z}; - - Acts::FreeVector freeParams; - try { - freeParams = Acts::estimateTrackParamsFromSeed( - pos0, 0., pos1, pos2, - Acts::Vector3{0., 0., mBz * Acts::UnitConstants::T}); - } catch (const std::exception& e) { - LOG(debug) << "ACTS seed parameter estimation failed: " << e.what(); - ++nInvalid; - continue; - } - - const double qOverP = freeParams[Acts::eFreeQOverP]; - if (!std::isfinite(qOverP) || std::abs(qOverP) < 1.e-12) { - ++nInvalid; - continue; - } - - Acts::Vector3 direction{ - freeParams[Acts::eFreeDir0], - freeParams[Acts::eFreeDir1], - freeParams[Acts::eFreeDir2]}; - - if (!direction.allFinite() || direction.norm() < 1.e-12) { - ++nInvalid; - continue; - } - direction.normalize(); - - std::array measurements; - std::vector sourceLinks; - std::vector surfaces; - sourceLinks.reserve(3); - surfaces.reserve(3); - - bool validMeasurements = true; - for (int i = 0; i < 3; ++i) { - const auto* sp = spacePoints[i]; - const auto* surface = sp->surface; - Acts::Vector3 position{sp->x, sp->y, sp->z}; - - auto localResult = surface->globalToLocal( - geometryContext, position, direction, - 1. * Acts::UnitConstants::mm); - - if (!localResult.ok()) { - validMeasurements = false; - break; - } - - auto& measurement = measurements[i]; - measurement.surface = surface; - measurement.local = localResult.value(); - - const double sigma = 100. * Acts::UnitConstants::um; - measurement.covariance.setZero(); - measurement.covariance(0, 0) = sigma * sigma; - measurement.covariance(1, 1) = sigma * sigma; - - sourceLinks.emplace_back(HitSourceLink{&measurement}); - surfaces.push_back(surface); - } - - if (!validMeasurements) { - ++nInvalid; - continue; - } - - Acts::Vector4 position4; - position4 << pos0.x(), pos0.y(), pos0.z(), 0.; - - Acts::BoundMatrix covariance = Acts::BoundMatrix::Zero(); - covariance(Acts::eBoundLoc0, Acts::eBoundLoc0) = - std::pow(1. * Acts::UnitConstants::mm, 2); - covariance(Acts::eBoundLoc1, Acts::eBoundLoc1) = - std::pow(1. * Acts::UnitConstants::mm, 2); - covariance(Acts::eBoundPhi, Acts::eBoundPhi) = 0.1 * 0.1; - covariance(Acts::eBoundTheta, Acts::eBoundTheta) = 0.1 * 0.1; - covariance(Acts::eBoundQOverP, Acts::eBoundQOverP) = - std::pow(0.5 * qOverP, 2); - covariance(Acts::eBoundTime, Acts::eBoundTime) = - std::pow(1. * Acts::UnitConstants::ns, 2); - - auto initialParameters = Acts::BoundTrackParameters::createCurvilinear( - position4, direction, qOverP, covariance, - Acts::ParticleHypothesis::pion()); - - auto result = fitter.fit( - sourceLinks.begin(), sourceLinks.end(), initialParameters, - fitterOptions, surfaces, tracks); - - if (!result.ok()) { - LOG(debug) << "ACTS Kalman fit failed: " << result.error().message(); - ++nFailed; - continue; - } - - const auto track = result.value(); - if (track.nMeasurements() < 3) { - ++nFailed; - continue; - } - - const auto innerState = track.innermostTrackState(); - if (!innerState.has_value()) { - ++nFailed; - continue; - } - const auto outerState = track.outermostTrackState(); - if (!innerState->hasSmoothed() || !outerState.hasSmoothed() || - !std::isfinite(track.chi2()) || track.chi2() < 0.f) { - ++nFailed; - continue; - } - - // Explicitly use the smoothed fit, rather than the original seed estimate. - typename Trajectory::ConstTrackStateProxy innerProxy{*innerState}; - typename Trajectory::ConstTrackStateProxy outerProxy{outerState}; - candidates.push_back({track.createParametersFromState(innerProxy), - track.createParametersFromState(outerProxy), - seed, static_cast(track.chi2())}); - ++nFitted; - LOG(debug) << "ACTS fitted seed: " - << "nMeasurements=" << track.nMeasurements() - << " chi2=" << track.chi2() - << " nDoF=" << track.nDoF(); - } - - // Accept lower-chi2 candidates first: multiple seeds can share clusters. - std::sort(candidates.begin(), candidates.end(), [](const auto& a, const auto& b) { - return a.chi2 < b.chi2; - }); - int nAccepted = 0; - for (const auto& candidate : candidates) { - nAccepted += acceptTrack(candidate.inner, candidate.outer, candidate.seed, candidate.chi2); - } - LOG(info) << "ACTS fitting: " << nFitted << " fitted, " - << nAccepted << " accepted, " << nFailed << " failed, " - << nInvalid << " invalid from " << mSeeds.size() << " seeds"; -} - -template -bool TrackerACTS::acceptTrack(const Acts::BoundTrackParameters& inner, - const Acts::BoundTrackParameters& outer, - const SeedACTS& seed, float chi2) -{ - if (!mTimeFrame || !std::isfinite(chi2) || chi2 < 0.f) { - return false; - } - - const std::array hits{seed.bottom, seed.middle, seed.top}; - std::array clusterIndices; - clusterIndices.fill(-1); - std::array internalIndices; - internalIndices.fill(-1); - uint32_t pattern = 0; - std::array keys{}; - - for (int i = 0; i < 3; ++i) { - const auto* hit = hits[i]; - if (!hit || hit->layer < 0 || hit->layer >= nLayers || hit->clusterId < 0 || - hit->rof < 0 || hit->rof >= mTimeFrame->getNrof(hit->layer)) { - return false; - } - const uint32_t mask = 1u << hit->layer; - if (pattern & mask) { // the TrackITSExt index array has one slot per layer - return false; - } - pattern |= mask; - keys[i] = (uint64_t(static_cast(hit->layer)) << 32) | - static_cast(hit->clusterId); - if (mAcceptedClusters.contains(keys[i])) { - return false; - } - // Keep the internal index for later MC lookup; publish the external one. - internalIndices[hit->layer] = hit->clusterId; - clusterIndices[hit->layer] = mTimeFrame->getClusterExternalIndex(hit->layer, hit->clusterId); - if (clusterIndices[hit->layer] < 0) { - return false; - } - } - - o2::track::TrackParCov innerO2, outerO2; - const auto& gctx = o2::acts::TrackingGeometryManager::instance().getNominalContext(); - if (!convertTrackParameters(inner, gctx, innerO2) || !convertTrackParameters(outer, gctx, outerO2)) { - return false; - } - - // Intersect the ROF timing intervals exactly as in ITS track acceptance. - const auto& timing = mTimeFrame->getROFOverlapTableView(); - auto nominal = timing.getLayer(hits[0]->layer).getROFTimeBounds(hits[0]->rof); - auto expanded = timing.getLayer(hits[0]->layer).getROFTimeBounds(hits[0]->rof, true); - bool nominalCompatible = true; - for (int i = 1; i < 3; ++i) { - const auto other = timing.getLayer(hits[i]->layer).getROFTimeBounds(hits[i]->rof); - const auto otherExpanded = timing.getLayer(hits[i]->layer).getROFTimeBounds(hits[i]->rof, true); - if (nominalCompatible) { - if (nominal.isCompatible(other)) { - nominal += other; - } else { - nominalCompatible = false; - } - } - if (!expanded.isCompatible(otherExpanded)) { - return false; - } - expanded += otherExpanded; - } - - o2::its::TrackITSExt o2Track{std::move(innerO2), static_cast(hits.size()), chi2, std::move(outerO2), clusterIndices}; - o2Track.setPattern(pattern); - o2Track.getTimeStamp() = (nominalCompatible ? nominal : expanded).makeSymmetrical(); - const float maxHalfROF = timing.getClockLayer().mROFLength * 0.5f; - if (o2Track.getTimeStamp().getTimeStampError() > maxHalfROF) { - o2Track.getTimeStamp().setTimeStampError(maxHalfROF); - } - - // Labels are computed once, after all ROFs, from the internal indices. - mTimeFrame->getTracks().emplace_back(std::move(o2Track)); - mAcceptedTrackClusterIndices.push_back(internalIndices); - for (const auto key : keys) { - mAcceptedClusters.insert(key); - } - return true; -} - -template -void TrackerACTS::computeTracksMClabels() -{ - if (!mTimeFrame || !mTimeFrame->hasMCinformation()) { - return; - } - - auto& tracks = mTimeFrame->getTracks(); - auto& labels = mTimeFrame->getTracksLabel(); - - if (tracks.size() != mAcceptedTrackClusterIndices.size()) { - LOGP(fatal, "ACTS MC labeling: {} tracks but {} internal cluster reference sets", - tracks.size(), mAcceptedTrackClusterIndices.size()); - return; - } - - // There must be exactly one label for each accepted track, in track order. - labels.clear(); - labels.reserve(tracks.size()); - - struct LabelCount { - o2::MCCompLabel label; - int count; - int lastHit; - }; - - for (size_t iTrack = 0; iTrack < tracks.size(); ++iTrack) { - auto& track = tracks[iTrack]; - const auto& indices = mAcceptedTrackClusterIndices[iTrack]; - std::vector counts; - int nHits = 0; - - for (int layer = 0; layer < nLayers; ++layer) { - if (indices[layer] < 0) { - continue; - } - ++nHits; - for (const auto& candidate : mTimeFrame->getClusterLabels(layer, indices[layer])) { - if (!candidate.isValid()) { - continue; - } - - auto it = std::find_if(counts.begin(), counts.end(), [&candidate](const auto& entry) { - return entry.label == candidate; - }); - if (it == counts.end()) { - counts.push_back({candidate, 1, nHits}); - } else if (it->lastHit != nHits) { - // Count a given truth label at most once per attached cluster. - ++it->count; - it->lastHit = nHits; - } - } - } - o2::MCCompLabel winner; - uint32_t pattern = track.getPattern() & 0x0000ffffu; - if (counts.empty()) { - winner.setFakeFlag(); - } else { - const auto best = std::max_element(counts.begin(), counts.end(), [](const auto& a, const auto& b) { - return a.count < b.count; - }); - winner = best->label; - if (best->count != nHits) { - winner.setFakeFlag(); - } - // Preserve the standard TrackITS pattern convention: bits [16, 31] - // indicate clusters compatible with the representative MC label. - for (int layer = 0; layer < nLayers; ++layer) { - if (indices[layer] < 0) { - continue; - } - for (const auto& candidate : mTimeFrame->getClusterLabels(layer, indices[layer])) { - if (candidate.isValid() && candidate == best->label) { - pattern |= (1u << (16 + layer)); - break; - } - } - } - } - track.setPattern(pattern); - labels.emplace_back(winner); - } - LOGP(info, "ACTS MC labeling: {} tracks, {} labels", tracks.size(), labels.size()); -} - -template -void TrackerACTS::clustersToTracks() -{ - if (!mTimeFrame) { - LOG(error) << "Cannot run TrackerACTS: No TimeFrame adopted"; - return; - } - - double totalTime = 0.; - mAcceptedClusters.clear(); - mAcceptedTrackClusterIndices.clear(); - LOG(info) << "==== TRK ACTS Tracking ===="; - - // Verify tracking geometry usage - auto& actsGeometryManager = o2::acts::TrackingGeometryManager::instance(); - const auto& trackingGeometry = actsGeometryManager.get(); - - // Dump tracking geometry to json file for debugging - (void)trackingGeometry; - - auto* trkGeometry = o2::trk::GeometryTGeo::Instance(); - trkGeometry->fillMatrixCache(o2::math_utils::bit2Mask(o2::math_utils::TransformType::L2G)); - const auto& surfaceIndex = actsGeometryManager.getIndex(*trkGeometry, 1.e-3, [trkGeometry](int sensorID) { - return std::string(trkGeometry->getMatrixPath(sensorID).Data()); - }); - - LOG(info) << "ACTS tracking geometry initialized with " << surfaceIndex.size() - << " mapped sensors on " << surfaceIndex.getNSurfaces() << " surfaces"; - try { - dumpGeometryOnce(); - } catch (const std::exception& e) { - LOG(error) << "ACTS geometry debug dump failed: " << e.what(); - } - if (mConfig.writeDebugTrees && mTimeFrameCounter == 0) { - initDebugOutput(); - } - - for (int iROF = 0; iROF < mTimeFrame->getNrof(0); ++iROF) { - LOG(info) << "Processing ROF " << iROF; - - mCurState = SpacePointBuilding; - totalTime += evaluateTask([this, iROF]() { buildSpacePoints(iROF); }, StateNames[mCurState]); - - mCurState = Seeding; - totalTime += evaluateTask([this]() { createSeeds(); }, StateNames[mCurState]); - - const auto nTracksBefore = mTimeFrame->getTracks().size(); - mCurState = TrackFinding; - totalTime += evaluateTask([this]() { findTracks(); }, StateNames[mCurState]); - if (isDebugTimeFrame()) { - mDebugROF.tf = static_cast(mTimeFrameCounter); - mDebugROF.rof = iROF; - mDebugROF.nSpacePoints = static_cast(mSpacePoints.size()); - mDebugROF.nSeeds = static_cast(mSeeds.size()); - mDebugROF.nAcceptedTracks = static_cast(mTimeFrame->getTracks().size() - nTracksBefore); - mROFTree->Fill(); - } - } - - // Run MC labeling once, after every accepted track has been published. - if (mTimeFrame->hasMCinformation()) { - computeTracksMClabels(); - } - - LOG(info) << "=== TimeFrame " << mTimeFrameCounter << " completed in: " << totalTime << " ms ==="; - if (isDebugTimeFrame()) { - mDebugFile->cd(); - if (mHistSpacePoints) { - mHistSpacePoints->Write("", TObject::kOverwrite); - } - mDebugFile->Write("", TObject::kOverwrite); - mDebugFile->Flush(); - } - ++mTimeFrameCounter; - mTotalTime += totalTime; -} - -template -void TrackerACTS::printSummary() const -{ - float avgTF = mTimeFrameCounter > 0 ? static_cast(mTotalTime) / mTimeFrameCounter : 0.f; - LOGP(info, "TrackerACTS summary: Processed {} TFs in TOT={:.2f} ms, AVG/TF={:.2f} ms", mTimeFrameCounter, mTotalTime, avgTF); -} - -// Explicit template instantiations -template class TrackerACTS<11>; -} // namespace o2::trk From 2b70e455753f8dbfa0ef7fdb36a59b192151b715 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nicol=C3=B2=20Jacazio?= Date: Thu, 5 Nov 2026 06:37:36 +0100 Subject: [PATCH 6/6] U --- .../GlobalReconstruction/reconstruction/src/TrackerACTS.cxx | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/src/TrackerACTS.cxx b/Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/src/TrackerACTS.cxx index 2d9751af5511f..b7f4f4cf1a11b 100644 --- a/Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/src/TrackerACTS.cxx +++ b/Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/src/TrackerACTS.cxx @@ -288,7 +288,7 @@ void TrackerACTS::initDebugOutput() mROFTree->Branch("nAcceptedTracks", &mDebugROF.nAcceptedTracks, "nAcceptedTracks/I"); mHistSpacePoints = std::make_unique("hSpacePoints", "Space points; x (cm); y (cm)", - 200, -100, 100, 200, -100, 100); + 200, -100, 100, 200, -100, 100); mHistSpacePoints->SetDirectory(nullptr); // Owned by TrackerACTS, not by TFile. LOG(info) << "ACTS debug trees enabled: " << mConfig.debugFileName; } @@ -480,7 +480,8 @@ void TrackerACTS::buildSpacePoints(int rof) mDebugCluster.externalId = mTimeFrame->getClusterExternalIndex(layer, clusterIndex); mDebugCluster.sensorId = sensorID; // 0: sensor ID missing; 1: ACTS surface missing; 2: mapped. - mDebugCluster.mappingStatus = !hasSensorID ? 0 : !surface ? 1 : 2; + mDebugCluster.mappingStatus = !hasSensorID ? 0 : !surface ? 1 + : 2; mDebugCluster.geometryId = surface ? surface->geometryId().value() : 0; mDebugCluster.x = cluster.xCoordinate; mDebugCluster.y = cluster.yCoordinate;