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..bfc83ef269536 100644 --- a/Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/include/ALICE3GlobalReconstruction/TrackerACTS.h +++ b/Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/include/ALICE3GlobalReconstruction/TrackerACTS.h @@ -19,105 +19,186 @@ #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 +#include +#include "DataFormatsITS/TrackITS.h" +#include +#include + +class TFile; +class TTree; +namespace Acts +{ +class Surface; +class BoundTrackParameters; +} // namespace Acts 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 - 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 + /// Minimum transverse momentum accepted by the seeding algorithms. + 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. + float cotThetaMax = std::sinh(4.0f); + + /// 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 - float deltaRMinTop = 5.f * Acts::UnitConstants::mm; ///< Min deltaR for middle-top - float deltaRMaxTop = 200.f * Acts::UnitConstants::mm; ///< Max deltaR for middle-top + /// 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 + /// 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 { public: TrackerACTS(); - ~TrackerACTS() - { - if (mHistSpacePoints) { - mHistSpacePoints->SaveAs("/tmp/mHistSpacePoints.C"); - delete mHistSpacePoints; - mHistSpacePoints = nullptr; - } - } + ~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: - 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); @@ -125,13 +206,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 +232,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/src/TrackerACTS.cxx b/Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/src/TrackerACTS.cxx index e870ee934816f..b7f4f4cf1a11b 100644 --- a/Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/src/TrackerACTS.cxx +++ b/Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/src/TrackerACTS.cxx @@ -15,34 +15,416 @@ /// \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" -#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 #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 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) +{ + 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); + // 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 +{ + // A negative limit means unlimited debug output; zero disables it. + 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(); + // 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"); + 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(); + // 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()); + }); + + std::ofstream out(mConfig.geometryDumpFileName); + if (!out) { + LOG(error) << "Cannot open ACTS geometry dump: " << mConfig.geometryDumpFileName; + return; + } + 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; + // 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["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; + // 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["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; + }); + 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}; + // 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; + 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) { + // Borrow the TimeFrame; its lifetime must cover all calls that use it. mTimeFrame = &tf; } @@ -54,19 +436,77 @@ 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(); + + // 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()); + }); + + 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); + + 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; - // Check that these are in global coordinates + // 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; @@ -74,181 +514,629 @@ 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 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; + 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 std::vector* m_sps; - }; + const auto last = static_cast(spacePoints.size()); + gridRanges.emplace_back(first, last); + } - // Wrap mSpacePoints in an Acts space-point container - SpacePointBackend backend{mSpacePoints}; + 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}); - // 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::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}); - // Configure the ACTS seed finder - const unsigned int maxSeeds = static_cast(mConfig.maxSeedsPerMiddleSP); - Acts::SeedFilterConfig filterConfig; - filterConfig.maxSeedsPerSpM = maxSeeds; + 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 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; + 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); } - 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(); + // 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; } - 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(); + // 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; 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}); + + 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; + + // 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; + 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}; - 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(); + 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 { - 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; + } + + 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}; + + // 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); + + if (!localResult.ok()) { + validMeasurements = false; + break; + } + + auto& measurement = measurements[i]; + measurement.surface = surface; + 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; + + 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; } - } - 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 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{}; + + // 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 || + 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; + // 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; + } + + // 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); + 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; + } + } + } + // 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()) { + 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 +1148,78 @@ 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) { + // 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(); + + // 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()); + }); + + 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(); + } + + // 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; - // 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 + const auto nTracksBefore = mTimeFrame->getTracks().size(); mCurState = TrackFinding; - totalTime += evaluateTask([this]() { findTracks(); }, - StateNames[mCurState]); + 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(); + } } - // 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 ==="; - + if (isDebugTimeFrame()) { + mDebugFile->cd(); + if (mHistSpacePoints) { + mHistSpacePoints->Write("", TObject::kOverwrite); + } + mDebugFile->Write("", TObject::kOverwrite); + mDebugFile->Flush(); + } ++mTimeFrameCounter; mTotalTime += totalTime; } @@ -297,8 +1228,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/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); 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