diff --git a/Detectors/ACTS/CMakeLists.txt b/Detectors/ACTS/CMakeLists.txt new file mode 100644 index 0000000000000..8a6c1adceed70 --- /dev/null +++ b/Detectors/ACTS/CMakeLists.txt @@ -0,0 +1,20 @@ +# Copyright 2019-2026 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. + +# The whole package is optional: ACTS is an optional O2 dependency +# (see dependencies/O2Dependencies.cmake) and standard builds do not ship it. +if(NOT Acts_FOUND) + message(STATUS "ACTS not found, skipping the O2 ACTS interface") + return() +endif() + +add_subdirectory(interface) +add_subdirectory(workflow) diff --git a/Detectors/ACTS/README.md b/Detectors/ACTS/README.md new file mode 100644 index 0000000000000..e8a907aa9e56d --- /dev/null +++ b/Detectors/ACTS/README.md @@ -0,0 +1,123 @@ +# ACTS tracking geometry in O2 + +Provides the shared, immutable `Acts::TrackingGeometry` that ACTS-based +reconstruction needs, built from the TGeo geometry O2 already has live in memory. + +ACTS is an **optional** O2 dependency (`dependencies/O2Dependencies.cmake`). Without +it this package produces no targets and `O2_WITH_ACTS` is not defined, so anything +using it must be guarded. + +## Libraries + +| Target | Contents | +| --- | --- | +| `O2::ACTSInterface` | `TrackingGeometryManager` (the provider), `ITrackingGeometryBuilder` (the detector-specific seam), `SurfaceIndexMap` (sensor ↔ surface lookup), `MagneticFieldAdapter` (`Acts::MagneticFieldProvider` over the O2 field). No DPL dependency. | +| `O2::ACTSWorkflow` | `ActsGeometryService`, the same provider reachable through the DPL `ServiceRegistry`. | +| `O2::ALICE3ACTS` | `Gen3BlueprintBuilder`, the ALICE 3 builder (under `Detectors/Upgrades/ALICE3/ACTS`). | + +## Using it from a DPL task + +The manager fetches nothing itself: it needs `gGeoManager` to be live, which is what +`GRPGeomRequest::Aligned` arranges. Note that a task that only used the material LUT +before will have been passing `GRPGeomRequest::None` and has to be switched. + +```cpp +// spec factory +auto ggRequest = std::make_shared( + false, false, false, true /*GRPMagField*/, false, + o2::base::GRPGeomRequest::Aligned, inputs, true); + +void init(InitContext& ic) +{ + o2::base::GRPGeomHelper::instance().setRequest(mCCDBReq); + + o2::alice3::Gen3BlueprintBuilder::Config cfg; + cfg.geometryConfigFile = ic.options().get("acts-geometry-config"); + o2::acts::TrackingGeometryManager::instance().setBuilder( + std::make_unique(cfg)); +} + +void run(ProcessingContext& pc) +{ + o2::base::GRPGeomHelper::instance().checkUpdates(pc); // gGeoManager now live + const auto& tg = o2::acts::TrackingGeometryManager::instance().get(); // built once +} +``` + +To reach it through the registry instead, put `o2::acts::defaultServicesWithActsGeometry()` +into the `DataProcessorSpec`'s `requiredServices` and use +`pc.services().get()`. Both paths hand out the same object. + +## Mapping O2 clusters onto ACTS surfaces + +`makeSurfaceIndex()` (or `TrackingGeometryManager::getIndex(cache, tol, pathProvider)`, +which caches per detector) maps every sensor of a `DetMatrixCache`-derived geometry helper +onto a sensitive ACTS surface. + +**Pass a `SensorPathProvider` whenever the detector has one.** With it, sensors are +identified by the TGeo node the surface was built from — an identity, so it cannot +mis-assign: + +```cpp +auto* trkGeo = o2::trk::GeometryTGeo::Instance(); +trkGeo->fillMatrixCache(o2::math_utils::bit2Mask(o2::math_utils::TransformType::L2G)); +const auto& index = mgr.getIndex(*trkGeo, 1e-3, + [trkGeo](int chipID) { return std::string(trkGeo->getMatrixPath(chipID).Data()); }); +``` + +Without it the fallback matches on the sensor's local-to-global transform, which is exact +where a transform identifies a sensor but not always: the ALICE 3 vertex-detector petals +are tube segments, so the three layers of a petal share both origin and rotation and +differ only in radius, which the matrix cache does not carry. Those are reported as +ambiguous rather than assigned arbitrarily, and `makeSurfaceIndex()` throws. One surface +per sensor is enforced in both modes. + +## The Gen3 geometry description + +`Detectors/Upgrades/ALICE3/ACTS/config/gen3_geometry_config.json` carries everything +geometry-specific: sensor name globs, region boundaries, per-subsystem clustering +tolerances, passive material structures and the pinned volume IDs. All lengths in mm. + +It is tuned for the layout `run_test.sh` generates — verified against an o2-sim geometry of +that layout: the sensor globs cover every sensitive volume family (and correctly exclude +`FT3Sensor_Inactive_*`), `*EOSCard*` matches the 504 end-of-stave volumes, the passive +cylinder radii and half-lengths match the real support volumes, and the resulting +volume/layer/surface table is identical to the one actsO2's reference geometry produces. + +One known imprecision, inherited and present for both geometries: `TRK_MID_CarbonSupport` +is declared with `halfZ = 1420` while the real `TRK_MID_CARBONSUPPORT` volume has +`dz = 1410 mm`. It only affects the z-extent of that passive layer's material receiver. + +## Material + +Gen3 construction takes no `IMaterialDecorator`, so an ACTS JSON material map is applied +afterwards via `setMaterialMapFile()`. Material maps are keyed on +`Acts::GeometryIdentifier`, so a structural change to a builder invalidates existing +maps; the manager raises rather than silently producing a material-free geometry. + +## Checking a geometry + +`Detectors/Upgrades/ALICE3/ACTS/macros/run_test.sh` does the whole thing: it runs +`o2-sim-serial-run5` with the ALICE 3 layout the ACTS chain is developed against, then +builds the tracking geometry from the result and checks every TRK chip against its ACTS +surface. Run it from a scratch directory — `o2-sim` writes into `$PWD`: + +```bash +mkdir -p /tmp/acts && cd /tmp/acts +bash /Detectors/Upgrades/ALICE3/ACTS/macros/run_test.sh +``` + +The geometry description it needs is shipped in this package +(`Detectors/Upgrades/ALICE3/ACTS/config/gen3_geometry_config.json`, installed to +`$O2_ROOT/share/Detectors/Upgrades/ALICE3/ACTS/config/`) and matches the layout the script +generates. Set `GEN3_CONFIG` to point at a different geometry's config. + +Knobs: `nEvents`, `generator`, `modules`, `GEN3_CONFIG`, `SKIP_ALIGNMENT`. The last one +defaults to 1 because `o2-sim` finishes by fetching alignment from CCDB, which needs a +valid alien token and aborts without one — after the geometry file has already been +written. The ACTS geometry is built from the ideal geometry, so skipping it costs nothing. + +The macro `CheckActsTrackingGeometry.C` can also be run on its own against an existing +`o2sim_geometry.root`; it prints the volume/layer/surface table, which is the authority for +the `(volume, layer)` keys that digitisation and seeding configurations use. It needs ACTS, +Eigen and `$O2_ROOT/include` on `ROOT_INCLUDE_PATH` — `run_test.sh` sets that up. diff --git a/Detectors/ACTS/interface/CMakeLists.txt b/Detectors/ACTS/interface/CMakeLists.txt new file mode 100644 index 0000000000000..e4d167fd73706 --- /dev/null +++ b/Detectors/ACTS/interface/CMakeLists.txt @@ -0,0 +1,26 @@ +# Copyright 2019-2026 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. + +o2_add_library(ACTSInterface + TARGETVARNAME targetName + SOURCES src/TrackingGeometryManager.cxx + src/SurfaceIndexMap.cxx + src/MagneticFieldAdapter.cxx + PUBLIC_LINK_LIBRARIES O2::DetectorsBase + O2::DetectorsCommonDataFormats + O2::Field + O2::Framework + ROOT::Geom + Acts::Core + Acts::PluginRoot + Acts::PluginJson) + +target_compile_definitions(${targetName} PUBLIC O2_WITH_ACTS) diff --git a/Detectors/ACTS/interface/include/ACTSInterface/ITrackingGeometryBuilder.h b/Detectors/ACTS/interface/include/ACTSInterface/ITrackingGeometryBuilder.h new file mode 100644 index 0000000000000..705dfdc18150a --- /dev/null +++ b/Detectors/ACTS/interface/include/ACTSInterface/ITrackingGeometryBuilder.h @@ -0,0 +1,79 @@ +// Copyright 2019-2026 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 ITrackingGeometryBuilder.h +/// \brief Detector-agnostic interface for building an Acts::TrackingGeometry from TGeo +/// \author Paolo Butti +/// + +#ifndef ALICEO2_ACTS_ITRACKINGGEOMETRYBUILDER_H +#define ALICEO2_ACTS_ITRACKINGGEOMETRYBUILDER_H + +#include +#include + +#include "Acts/Geometry/GeometryContext.hpp" +#include "Acts/Geometry/TrackingGeometry.hpp" +#include "Acts/Surfaces/SurfacePlacementBase.hpp" + +#include "ACTSInterface/SurfaceIndexMap.h" + +class TGeoManager; + +namespace o2::acts +{ + +/// Everything a builder produces, kept together because the parts have a +/// lifetime dependency on each other. +struct TrackingGeometryOutput { + /// Non-const so the manager can still decorate it with material before + /// handing it out; it is published as shared_ptr afterwards. + std::shared_ptr geometry; + + /// Detector elements backing the sensitive surfaces. + /// + /// Acts::TrackingGeometry does NOT own these: Acts::Surface keeps a raw, + /// non-owning back-pointer to its placement (Surface::surfacePlacement()). + /// Destroying this store while \a geometry is alive is a dangling-pointer bug, + /// so it must be kept for at least as long as the geometry. + std::vector> elementStore; + + /// Sensor <-> surface lookup. May be empty if the builder cannot provide one. + SurfaceIndexMap index; +}; + +/// Builds an Acts::TrackingGeometry from a TGeo tree. +/// +/// Implementations are detector specific and live with their detector (see +/// o2::alice3::Gen3BlueprintBuilder). They are injected into +/// TrackingGeometryManager at runtime, so this library never depends on any +/// concrete detector. +class ITrackingGeometryBuilder +{ + public: + virtual ~ITrackingGeometryBuilder() = default; + + /// Build from an already-loaded TGeo tree. + /// + /// \param tgeo the geometry manager to read; implementations must not import + /// or otherwise replace it, since it is normally O2's live + /// gGeoManager shared with the rest of the workflow + /// \param gctx nominal context used while placing surfaces + virtual TrackingGeometryOutput build(TGeoManager& tgeo, const Acts::GeometryContext& gctx) = 0; + + /// Short name, used in log messages. + virtual const char* getName() const = 0; +}; + +} // namespace o2::acts + +#endif // ALICEO2_ACTS_ITRACKINGGEOMETRYBUILDER_H diff --git a/Detectors/ACTS/interface/include/ACTSInterface/MagneticFieldAdapter.h b/Detectors/ACTS/interface/include/ACTSInterface/MagneticFieldAdapter.h new file mode 100644 index 0000000000000..881dd2fd2b1e1 --- /dev/null +++ b/Detectors/ACTS/interface/include/ACTSInterface/MagneticFieldAdapter.h @@ -0,0 +1,67 @@ +// Copyright 2019-2026 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 MagneticFieldAdapter.h +/// \brief Acts::MagneticFieldProvider backed by the O2 field +/// \author Paolo Butti +/// + +#ifndef ALICEO2_ACTS_MAGNETICFIELDADAPTER_H +#define ALICEO2_ACTS_MAGNETICFIELDADAPTER_H + +#include "Acts/MagneticField/MagneticFieldProvider.hpp" + +namespace o2::field +{ +class MagneticField; +} + +namespace o2::acts +{ + +/// Exposes the O2 magnetic field to ACTS. +/// +/// Holds no field of its own: it reads the process-global field that +/// o2::base::Propagator::initFieldFromGRP() installs into +/// TGeoGlobalMagField::Instance() from the GRPMagField CCDB object. The +/// underlying object stays valid across field rescalings -- GRPGeomHelper +/// rescales the existing MagneticField in place rather than replacing it -- so +/// an adapter built once keeps returning up-to-date values. +/// +/// Units: O2 works in kGauss and cm, ACTS in its native tesla/mm system, and the +/// conversion is applied here. +class MagneticFieldAdapter : public Acts::MagneticFieldProvider +{ + public: + /// Use the field currently installed in TGeoGlobalMagField. + /// \throw std::runtime_error if no field has been initialised + MagneticFieldAdapter(); + + /// Use an explicitly provided field. + explicit MagneticFieldAdapter(o2::field::MagneticField* field); + + Cache makeCache(const Acts::MagneticFieldContext& mctx) const final; + + Acts::Result getField(const Acts::Vector3& position, Cache& cache) const final; + + const o2::field::MagneticField* getO2Field() const { return mField; } + + private: + struct CacheImpl { + }; + + o2::field::MagneticField* mField = nullptr; +}; + +} // namespace o2::acts + +#endif // ALICEO2_ACTS_MAGNETICFIELDADAPTER_H diff --git a/Detectors/ACTS/interface/include/ACTSInterface/SurfaceIndexMap.h b/Detectors/ACTS/interface/include/ACTSInterface/SurfaceIndexMap.h new file mode 100644 index 0000000000000..26b83c1737bcf --- /dev/null +++ b/Detectors/ACTS/interface/include/ACTSInterface/SurfaceIndexMap.h @@ -0,0 +1,130 @@ +// Copyright 2019-2026 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 SurfaceIndexMap.h +/// \brief Lookup between O2 sensor indices and ACTS surfaces +/// \author Paolo Butti +/// + +#ifndef ALICEO2_ACTS_SURFACEINDEXMAP_H +#define ALICEO2_ACTS_SURFACEINDEXMAP_H + +#include +#include +#include +#include +#include + +#include "Acts/Geometry/GeometryIdentifier.hpp" + +namespace Acts +{ +class Surface; +class TrackingGeometry; +class GeometryContext; +} // namespace Acts + +namespace o2::detectors +{ +class DetMatrixCache; +} + +namespace o2::acts +{ + +/// Lookup between an O2 sensor index (chip ID, as used by the detector's +/// DetMatrixCache-derived geometry helper) and the ACTS surface it sits on. +/// +/// The relation is many-to-one, not one-to-one. A planar sensor becomes its own +/// ACTS surface, but a sub-detector modelled as a whole cylindrical or disc +/// layer -- as the ALICE 3 vertex detector is -- becomes a single surface that +/// backs every chip of that layer. Hence getSurface() takes a sensor and +/// getSensorIDs() gives back a list. +/// +/// ACTS assigns the Acts::GeometryIdentifier inside +/// Acts::Blueprint::construct(), so this map can only be filled once the +/// TrackingGeometry is complete. See makeSurfaceIndex(). +class SurfaceIndexMap +{ + public: + static constexpr int InvalidSensorID = -1; + + /// \throw std::runtime_error if the sensor is already mapped + void add(int sensorID, const Acts::Surface& surface); + + /// \return nullptr if the sensor is not mapped + const Acts::Surface* getSurface(int sensorID) const; + /// \return a default-constructed identifier if the sensor is not mapped + Acts::GeometryIdentifier getGeometryId(int sensorID) const; + + /// All sensors sitting on the given surface, in increasing sensor order. + /// \return an empty vector if the surface backs no sensor + const std::vector& getSensorIDs(Acts::GeometryIdentifier geoId) const; + + bool empty() const { return mBySensor.empty(); } + /// Number of mapped sensors. + std::size_t size() const { return mBySensor.size(); } + /// Number of distinct surfaces the sensors map onto. + std::size_t getNSurfaces() const { return mByGeoId.size(); } + void clear(); + + private: + struct Entry { + const Acts::Surface* surface = nullptr; + Acts::GeometryIdentifier geoId{}; + }; + + std::unordered_map mBySensor; + std::unordered_map> mByGeoId; +}; + +/// Full TGeo node path of one sensor, e.g. o2::trk::GeometryTGeo::getMatrixPath(). +using SensorPathProvider = std::function; + +/// Build the sensor <-> surface index. +/// +/// Two strategies, and \b which \b one \b is \b used \b matters: +/// +/// - **With \p pathProvider (exact, preferred).** Each sensor is resolved to its +/// TGeo node and matched against the node every ACTS surface was built from +/// (ActsPlugins::TGeoDetectorElement::tgeoNode()). This is an identity, not a +/// measurement, so it cannot mis-assign. Repeated placements of one volume +/// share a node and are separated by their accumulated global transform. +/// +/// - **Without it (geometric fallback).** Sensors are matched on their +/// local-to-global transform. Exact where a transform identifies a sensor, but +/// it does not always: an ALICE 3 vertex-detector petal is a tube segment, so +/// all three of its layers share one origin *and* one rotation and differ only +/// in radius, which the matrix cache does not carry. Such sensors are reported +/// as ambiguous rather than assigned arbitrarily. +/// +/// Both read the same TGeo tree that the surfaces were built from, so a successful +/// match is exact rather than approximate; \p toleranceCm only guards floating-point +/// round-trips. +/// +/// \param geometry the constructed tracking geometry +/// \param gctx context used to evaluate the surface transforms +/// \param cache detector geometry helper; requires a filled L2G cache +/// (call fillMatrixCache(bit2Mask(o2::math_utils::TransformType::L2G))) +/// \param toleranceCm matching tolerance, in O2 units (cm) +/// \param pathProvider optional TGeo path lookup, as described above +/// \throw std::runtime_error if a sensor cannot be matched, if a match is +/// ambiguous, or if two sensors would claim the same surface +SurfaceIndexMap makeSurfaceIndex(const Acts::TrackingGeometry& geometry, + const Acts::GeometryContext& gctx, + const o2::detectors::DetMatrixCache& cache, + double toleranceCm = 1e-3, + SensorPathProvider pathProvider = {}); + +} // namespace o2::acts + +#endif // ALICEO2_ACTS_SURFACEINDEXMAP_H diff --git a/Detectors/ACTS/interface/include/ACTSInterface/TrackingGeometryManager.h b/Detectors/ACTS/interface/include/ACTSInterface/TrackingGeometryManager.h new file mode 100644 index 0000000000000..557a3f9ffb6fd --- /dev/null +++ b/Detectors/ACTS/interface/include/ACTSInterface/TrackingGeometryManager.h @@ -0,0 +1,145 @@ +// Copyright 2019-2026 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 TrackingGeometryManager.h +/// \brief Process-wide provider of the ACTS tracking geometry +/// \author Paolo Butti +/// + +#ifndef ALICEO2_ACTS_TRACKINGGEOMETRYMANAGER_H +#define ALICEO2_ACTS_TRACKINGGEOMETRYMANAGER_H + +#include +#include +#include +#include + +#include "Acts/Utilities/Logger.hpp" + +#include "ACTSInterface/ITrackingGeometryBuilder.h" + +namespace o2::acts +{ + +/// Provider of the shared, immutable Acts::TrackingGeometry. +/// +/// Follows the same pattern as o2::base::Propagator::Instance(), +/// o2::base::GRPGeomHelper::instance() and o2::its::GeometryTGeo::Instance(): +/// a process-wide lazily-built singleton fed from the geometry O2 already has +/// live in memory. It fetches nothing by itself -- the caller is responsible for +/// making gGeoManager available first, normally by declaring +/// o2::base::GRPGeomRequest::Aligned and calling +/// o2::base::GRPGeomHelper::instance().checkUpdates(pc). +/// +/// Typical use in a DPL task: +/// \code +/// void init(InitContext&) { +/// o2::base::GRPGeomHelper::instance().setRequest(mCCDBReq); +/// TrackingGeometryManager::instance().setBuilder(std::make_unique(cfg)); +/// } +/// void run(ProcessingContext& pc) { +/// o2::base::GRPGeomHelper::instance().checkUpdates(pc); // gGeoManager now live +/// const auto& tg = TrackingGeometryManager::instance().get(); // built once +/// } +/// \endcode +/// +/// Thread safety: build() is not re-entrant and is expected to run once, before +/// any parallel processing. Once built the geometry is immutable and all +/// accessors are safe to call concurrently. +class TrackingGeometryManager +{ + public: + static TrackingGeometryManager& instance(); + + /// Install the detector-specific builder. Must be called before the first build. + void setBuilder(std::unique_ptr builder); + bool hasBuilder() const { return mBuilder != nullptr; } + + /// Path to an ACTS JSON material map, applied after construction. + /// Empty (the default) means no material decoration. + void setMaterialMapFile(std::string fileName) { mMaterialMapFile = std::move(fileName); } + const std::string& getMaterialMapFile() const { return mMaterialMapFile; } + + void setLogLevel(Acts::Logging::Level level) { mLogLevel = level; } + + /// Build from the live gGeoManager. + /// \throw std::runtime_error if no builder is set or no geometry is loaded + void build(); + + /// Build from a TGeo file, for macros and tests running outside a workflow. + /// Loads the file through o2::base::GeometryManager, so it ends up in + /// gGeoManager exactly as it would inside a workflow. + void buildFromFile(std::string_view geomFileName); + + bool isBuilt() const { return mGeometry != nullptr; } + void clear(); + + /// Builds on first call if not built yet. + std::shared_ptr getShared(); + const Acts::TrackingGeometry& get(); + + /// Sensor <-> surface lookup, as filled by the builder. Empty if the builder + /// did not provide one -- which is the normal case when the geometry spans + /// several detectors, since each has its own sensor numbering. Use the + /// overload below for those. + const SurfaceIndexMap& getIndex(); + + /// Sensor <-> surface lookup for one detector, built from its geometry helper + /// on first request and cached per detector afterwards. + /// + /// Pass \p pathProvider whenever the detector offers one -- it makes the match an + /// identity rather than a geometric inference, see makeSurfaceIndex(): + /// \code + /// auto* trkGeo = o2::trk::GeometryTGeo::Instance(); + /// mgr.getIndex(*trkGeo, 1e-3, + /// [trkGeo](int id) { return std::string(trkGeo->getMatrixPath(id).Data()); }); + /// \endcode + /// + /// \param cache a DetMatrixCache-derived helper with a filled L2G cache, + /// e.g. o2::trk::GeometryTGeo::Instance() + /// \param toleranceCm matching tolerance, see makeSurfaceIndex() + /// \param pathProvider optional TGeo path lookup for exact identification + const SurfaceIndexMap& getIndex(const o2::detectors::DetMatrixCache& cache, + double toleranceCm = 1e-3, + SensorPathProvider pathProvider = {}); + + /// Context the geometry was built with. Carries no alignment payload: O2 bakes + /// alignment into the TGeo matrices before the ACTS geometry is built, so the + /// nominal transforms are already the aligned ones. + const Acts::GeometryContext& getNominalContext() const { return mNominalContext; } + + /// Number of surfaces that received material during decoration. Zero when no + /// material map was configured. + std::size_t getNDecoratedSurfaces() const { return mNDecoratedSurfaces; } + + private: + TrackingGeometryManager(); + + void doBuild(); + void applyMaterialMapImpl(Acts::TrackingGeometry& geometry); + + std::unique_ptr mBuilder; + std::string mMaterialMapFile; + Acts::Logging::Level mLogLevel = Acts::Logging::INFO; + + Acts::GeometryContext mNominalContext; + std::shared_ptr mGeometry; + std::vector> mElementStore; + SurfaceIndexMap mIndex; + /// Per-detector indices, keyed by o2::detectors::DetID. + std::unordered_map mDetIndices; + std::size_t mNDecoratedSurfaces = 0; +}; + +} // namespace o2::acts + +#endif // ALICEO2_ACTS_TRACKINGGEOMETRYMANAGER_H diff --git a/Detectors/ACTS/interface/src/MagneticFieldAdapter.cxx b/Detectors/ACTS/interface/src/MagneticFieldAdapter.cxx new file mode 100644 index 0000000000000..f2e23ba2be554 --- /dev/null +++ b/Detectors/ACTS/interface/src/MagneticFieldAdapter.cxx @@ -0,0 +1,72 @@ +// Copyright 2019-2026 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 MagneticFieldAdapter.cxx +/// \author Paolo Butti +/// + +#include "ACTSInterface/MagneticFieldAdapter.h" + +#include + +#include + +#include "Acts/Definitions/Units.hpp" +#include "Field/MagneticField.h" +#include "Framework/Logger.h" + +namespace o2::acts +{ + +namespace +{ +/// O2 stores the field in kGauss, ACTS in its native unit system. +constexpr double kGaussToActs = 0.1 * Acts::UnitConstants::T; +} // namespace + +MagneticFieldAdapter::MagneticFieldAdapter() + : MagneticFieldAdapter(static_cast(TGeoGlobalMagField::Instance()->GetField())) +{ +} + +MagneticFieldAdapter::MagneticFieldAdapter(o2::field::MagneticField* field) : mField(field) +{ + if (mField == nullptr) { + throw std::runtime_error( + "o2::acts::MagneticFieldAdapter: no O2 magnetic field available. Request GRPMagField in the " + "GRPGeomRequest of the workflow, or initialise the field explicitly with " + "o2::base::Propagator::initFieldFromGRP()."); + } + LOG(info) << "ACTS magnetic field adapter bound to O2 field '" << mField->GetName() << "'"; +} + +Acts::MagneticFieldProvider::Cache MagneticFieldAdapter::makeCache(const Acts::MagneticFieldContext& /*mctx*/) const +{ + // The O2 field parametrisation keeps its own internal caches and offers no + // per-client cache handle, so there is nothing to carry here. + return Acts::MagneticFieldProvider::Cache(std::in_place_type); +} + +Acts::Result MagneticFieldAdapter::getField(const Acts::Vector3& position, + Cache& /*cache*/) const +{ + const double xyz[3] = {position[0] / Acts::UnitConstants::cm, + position[1] / Acts::UnitConstants::cm, + position[2] / Acts::UnitConstants::cm}; + double bxyz[3] = {0., 0., 0.}; + mField->Field(xyz, bxyz); + + return Acts::Result::success( + Acts::Vector3{bxyz[0] * kGaussToActs, bxyz[1] * kGaussToActs, bxyz[2] * kGaussToActs}); +} + +} // namespace o2::acts diff --git a/Detectors/ACTS/interface/src/SurfaceIndexMap.cxx b/Detectors/ACTS/interface/src/SurfaceIndexMap.cxx new file mode 100644 index 0000000000000..f8d74106d6ecd --- /dev/null +++ b/Detectors/ACTS/interface/src/SurfaceIndexMap.cxx @@ -0,0 +1,390 @@ +// Copyright 2019-2026 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 SurfaceIndexMap.cxx +/// \author Paolo Butti +/// + +#include "ACTSInterface/SurfaceIndexMap.h" + +#include +#include +#include +#include +#include + +#include "Acts/Definitions/Units.hpp" +#include "Acts/Geometry/GeometryContext.hpp" +#include "Acts/Geometry/TrackingGeometry.hpp" +#include "Acts/Surfaces/BoundaryTolerance.hpp" +#include "Acts/Surfaces/Surface.hpp" +#include "Acts/Surfaces/SurfacePlacementBase.hpp" +#include "ActsPlugins/Root/TGeoDetectorElement.hpp" + +#include +#include +#include + +#include "DetectorsCommonDataFormats/DetMatrixCache.h" +#include "Framework/Logger.h" + +namespace o2::acts +{ + +namespace +{ +const std::vector gNoSensors{}; +} // namespace + +void SurfaceIndexMap::add(int sensorID, const Acts::Surface& surface) +{ + const auto geoId = surface.geometryId(); + auto [it, inserted] = mBySensor.try_emplace(sensorID, Entry{&surface, geoId}); + if (!inserted) { + throw std::runtime_error("o2::acts::SurfaceIndexMap: sensor " + std::to_string(sensorID) + + " is already mapped"); + } + mByGeoId[geoId].push_back(sensorID); +} + +const Acts::Surface* SurfaceIndexMap::getSurface(int sensorID) const +{ + const auto it = mBySensor.find(sensorID); + return it == mBySensor.end() ? nullptr : it->second.surface; +} + +Acts::GeometryIdentifier SurfaceIndexMap::getGeometryId(int sensorID) const +{ + const auto it = mBySensor.find(sensorID); + return it == mBySensor.end() ? Acts::GeometryIdentifier{} : it->second.geoId; +} + +const std::vector& SurfaceIndexMap::getSensorIDs(Acts::GeometryIdentifier geoId) const +{ + const auto it = mByGeoId.find(geoId); + return it == mByGeoId.end() ? gNoSensors : it->second; +} + +void SurfaceIndexMap::clear() +{ + mBySensor.clear(); + mByGeoId.clear(); +} + +namespace +{ + +/// Uniform grid hash over per-sensor surface centres, so matching stays linear +/// in the number of sensors instead of quadratic. +class CentreLookup +{ + public: + explicit CentreLookup(double cellSize) : mCellSize(cellSize) {} + + void add(const Acts::Surface* surface, const Acts::Vector3& position) + { + mEntries.push_back({surface, position}); + mCells[cellOf(position)].push_back(mEntries.size() - 1); + } + + /// \return the closest surface within \p tolerance, or nullptr + const Acts::Surface* findNearest(const Acts::Vector3& position, double tolerance) const + { + const auto cell = cellOf(position); + const Acts::Surface* best = nullptr; + double bestDist2 = tolerance * tolerance; + for (std::int64_t dx = -1; dx <= 1; ++dx) { + for (std::int64_t dy = -1; dy <= 1; ++dy) { + for (std::int64_t dz = -1; dz <= 1; ++dz) { + const auto it = mCells.find(Cell{cell.x + dx, cell.y + dy, cell.z + dz}); + if (it == mCells.end()) { + continue; + } + for (const auto idx : it->second) { + const double dist2 = (mEntries[idx].position - position).squaredNorm(); + if (dist2 <= bestDist2) { + bestDist2 = dist2; + best = mEntries[idx].surface; + } + } + } + } + } + return best; + } + + std::size_t size() const { return mEntries.size(); } + + private: + struct Entry { + const Acts::Surface* surface = nullptr; + Acts::Vector3 position; + }; + struct Cell { + std::int64_t x = 0, y = 0, z = 0; + bool operator==(const Cell&) const = default; + }; + struct CellHash { + std::size_t operator()(const Cell& c) const noexcept + { + std::size_t h = static_cast(c.x) * 0x9e3779b97f4a7c15ULL; + h ^= static_cast(c.y) + 0x9e3779b97f4a7c15ULL + (h << 6) + (h >> 2); + h ^= static_cast(c.z) + 0x9e3779b97f4a7c15ULL + (h << 6) + (h >> 2); + return h; + } + }; + + Cell cellOf(const Acts::Vector3& p) const + { + return Cell{static_cast(std::floor(p[0] / mCellSize)), + static_cast(std::floor(p[1] / mCellSize)), + static_cast(std::floor(p[2] / mCellSize))}; + } + + double mCellSize = 1.; + std::vector mEntries; + std::unordered_map, CellHash> mCells; +}; + +} // namespace + +namespace +{ + +/// Collect the sensitive surfaces, split by whether their centre identifies them. +struct SensitiveSurfaces { + CentreLookup byCentre; + std::vector layerSurfaces; + std::unordered_map> byNode; + std::size_t nWithoutNode = 0; + + explicit SensitiveSurfaces(double cellSize) : byCentre(cellSize) {} +}; + +SensitiveSurfaces collectSensitive(const Acts::TrackingGeometry& geometry, + const Acts::GeometryContext& gctx, double cellSize) +{ + SensitiveSurfaces out(cellSize); + geometry.apply([&](const Acts::Surface& surface) { + const auto* placement = surface.surfacePlacement(); + if (placement == nullptr || !placement->isSensitive()) { + return; + } + const auto* tgeoElement = dynamic_cast(placement); + if (tgeoElement != nullptr) { + out.byNode[&tgeoElement->tgeoNode()].push_back(&surface); + } else { + ++out.nWithoutNode; + } + const auto type = surface.type(); + if (type == Acts::Surface::SurfaceType::Cylinder || type == Acts::Surface::SurfaceType::Disc) { + out.layerSurfaces.push_back(&surface); + } else { + out.byCentre.add(&surface, surface.center(gctx)); + } + }); + return out; +} + +/// Global translation of an O2 sensor, in ACTS units, plus its rotation. +struct SensorPlacement { + Acts::Vector3 position; + Acts::RotationMatrix3 rotation; +}; + +SensorPlacement sensorPlacementOf(const o2::detectors::DetMatrixCache& cache, int sensorID) +{ + double rot[9] = {0.}; + double tra[3] = {0.}; + cache.getMatrixL2G(sensorID).GetComponents(rot[0], rot[1], rot[2], tra[0], + rot[3], rot[4], rot[5], tra[1], + rot[6], rot[7], rot[8], tra[2]); + SensorPlacement out; + out.position = Acts::Vector3{tra[0] * Acts::UnitConstants::cm, tra[1] * Acts::UnitConstants::cm, + tra[2] * Acts::UnitConstants::cm}; + out.rotation << rot[0], rot[1], rot[2], rot[3], rot[4], rot[5], rot[6], rot[7], rot[8]; + return out; +} + +/// Resolve a sensor through its TGeo node. Exact: the node is an identity, not a +/// measurement. Repeated placements of one volume share a node, so those are +/// separated by the accumulated global translation. +/// \return the surface, or nullptr with \p ambiguous set when several remain +const Acts::Surface* matchByNode(const SensitiveSurfaces& surfaces, const std::string& path, + const Acts::GeometryContext& gctx, const SensorPlacement& placement, + double tolerance, bool& ambiguous) +{ + ambiguous = false; + if (gGeoManager == nullptr) { + throw std::runtime_error("o2::acts::makeSurfaceIndex: no TGeo geometry loaded"); + } + gGeoManager->PushPath(); + const bool ok = gGeoManager->cd(path.c_str()); + const TGeoNode* node = ok ? gGeoManager->GetCurrentNode() : nullptr; + gGeoManager->PopPath(); + if (node == nullptr) { + return nullptr; + } + + const auto it = surfaces.byNode.find(node); + if (it == surfaces.byNode.end()) { + return nullptr; + } + if (it->second.size() == 1) { + return it->second.front(); + } + const Acts::Surface* best = nullptr; + for (const auto* surface : it->second) { + if ((surface->localToGlobalTransform(gctx).translation() - placement.position).norm() <= tolerance) { + if (best != nullptr) { + ambiguous = true; + return nullptr; + } + best = surface; + } + } + return best; +} + +/// Resolve a sensor from its transform alone. Used when no TGeo path is available. +/// \return the surface, or nullptr with \p ambiguous set when several match +const Acts::Surface* matchByPlacement(const SensitiveSurfaces& surfaces, + const Acts::GeometryContext& gctx, + const SensorPlacement& placement, double tolerance, + bool& ambiguous) +{ + ambiguous = false; + + // Rotations are dimensionless; a generous bound on the TGeoHMatrix -> Eigen round trip. + constexpr double kRotationTolerance = 1e-9; + + if (const auto* surface = surfaces.byCentre.findNearest(placement.position, tolerance); + surface != nullptr) { + return surface; + } + + // A sensor modelled as a whole cylinder or disc: its origin is generally not a + // useful key (an ALICE 3 vertex-detector petal is a tube segment, so its origin + // is on the beam axis), so the full transform is compared instead. Every match + // is collected rather than taking the first: the three layers of one petal share + // both origin and rotation and differ only in radius, which the matrix cache does + // not carry, so those are genuinely unresolvable this way and must be reported. + const Acts::Surface* best = nullptr; + for (const auto* surface : surfaces.layerSurfaces) { + const auto& transform = surface->localToGlobalTransform(gctx); + if ((transform.translation() - placement.position).norm() <= tolerance && + (transform.rotation() - placement.rotation).norm() <= kRotationTolerance) { + if (best != nullptr) { + ambiguous = true; + return nullptr; + } + best = surface; + } + } + return best; +} + +} // namespace + +SurfaceIndexMap makeSurfaceIndex(const Acts::TrackingGeometry& geometry, + const Acts::GeometryContext& gctx, + const o2::detectors::DetMatrixCache& cache, + double toleranceCm, + SensorPathProvider pathProvider) +{ + if (!cache.isBuilt()) { + throw std::runtime_error("o2::acts::makeSurfaceIndex: the detector matrix cache is not built"); + } + if (cache.getCacheL2G().getSize() == 0) { + throw std::runtime_error( + "o2::acts::makeSurfaceIndex: the L2G matrix cache is empty. Call fillMatrixCache() with the " + "L2G bit set before building the index."); + } + + // ACTS works in mm, O2 in cm; Acts::UnitConstants::cm is exactly that factor. + const double tolerance = toleranceCm * Acts::UnitConstants::cm; + + const auto surfaces = collectSensitive(geometry, gctx, std::max(tolerance, 1e-6)); + const int nSensors = cache.getSize(); + + LOGP(info, + "Matching {} {} sensors against {} sensitive ACTS surfaces ({} per-sensor, {} whole-layer) " + "by {} (tolerance {} cm)", + nSensors, cache.getDetID().getName(), surfaces.byCentre.size() + surfaces.layerSurfaces.size(), + surfaces.byCentre.size(), surfaces.layerSurfaces.size(), + pathProvider ? "TGeo node" : "placement", toleranceCm); + if (pathProvider && surfaces.nWithoutNode != 0) { + LOGP(warn, "{} sensitive surfaces have no TGeo detector element and cannot be matched by node", + surfaces.nWithoutNode); + } + + SurfaceIndexMap index; + std::size_t nUnmatched = 0; + std::size_t nAmbiguous = 0; + std::unordered_map claimedBy; + std::ostringstream failures; + int nReported = 0; + + for (int sensorID = 0; sensorID < nSensors; ++sensorID) { + const auto placement = sensorPlacementOf(cache, sensorID); + + bool ambiguous = false; + const Acts::Surface* match = + pathProvider + ? matchByNode(surfaces, pathProvider(sensorID), gctx, placement, tolerance, ambiguous) + : matchByPlacement(surfaces, gctx, placement, tolerance, ambiguous); + + if (match == nullptr) { + ambiguous ? ++nAmbiguous : ++nUnmatched; + if (nReported < 10) { + failures << "\n sensor " << sensorID << " at (" + << placement.position[0] / Acts::UnitConstants::cm << ", " + << placement.position[1] / Acts::UnitConstants::cm << ", " + << placement.position[2] / Acts::UnitConstants::cm << ") cm " + << (ambiguous ? "matches several sensitive ACTS surfaces" : "matches no sensitive ACTS surface"); + ++nReported; + } + continue; + } + + // One surface per sensor. Relaxing this hides exactly the failure the + // ambiguity check above is there to catch, so it is enforced. + const auto geoId = match->geometryId(); + const auto [it, inserted] = claimedBy.try_emplace(geoId, sensorID); + if (!inserted) { + throw std::runtime_error( + "o2::acts::makeSurfaceIndex: sensors " + std::to_string(it->second) + " and " + + std::to_string(sensorID) + " both resolve to the ACTS surface with geometry id " + + std::to_string(geoId.value()) + + ". They are indistinguishable by the information available, so the index would be wrong."); + } + index.add(sensorID, *match); + } + + if (nUnmatched != 0 || nAmbiguous != 0) { + std::string hint; + if (nAmbiguous != 0 && !pathProvider) { + hint = + " Pass a SensorPathProvider (e.g. GeometryTGeo::getMatrixPath) so sensors are identified by " + "their TGeo node instead of their placement."; + } + throw std::runtime_error("o2::acts::makeSurfaceIndex: of " + std::to_string(nSensors) + + " sensors, " + std::to_string(nUnmatched) + " matched no surface and " + + std::to_string(nAmbiguous) + " were ambiguous." + hint + failures.str()); + } + + LOGP(info, "Sensor index built: {} sensors on {} surfaces", index.size(), index.getNSurfaces()); + + return index; +} + +} // namespace o2::acts diff --git a/Detectors/ACTS/interface/src/TrackingGeometryManager.cxx b/Detectors/ACTS/interface/src/TrackingGeometryManager.cxx new file mode 100644 index 0000000000000..29a83210952d0 --- /dev/null +++ b/Detectors/ACTS/interface/src/TrackingGeometryManager.cxx @@ -0,0 +1,231 @@ +// Copyright 2019-2026 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 TrackingGeometryManager.cxx +/// \author Paolo Butti +/// + +#include "ACTSInterface/TrackingGeometryManager.h" + +#include +#include +#include + +#include +#include + +#include "Acts/Geometry/GeometryIdentifier.hpp" +#include "Acts/Geometry/TrackingVolume.hpp" +#include "Acts/Surfaces/Surface.hpp" +#include "ActsPlugins/Json/JsonMaterialDecorator.hpp" +#include "ActsPlugins/Json/MaterialMapJsonConverter.hpp" + +#include "CommonUtils/NameConf.h" +#include "DetectorsBase/GeometryManager.h" +#include "DetectorsCommonDataFormats/DetMatrixCache.h" +#include "Framework/Logger.h" + +namespace o2::acts +{ + +namespace +{ +/// Whether the file holds a geometry under one of the names +/// o2::base::GeometryManager::loadGeometry() accepts. +bool hasO2GeometryObject(std::string_view fileName) +{ + TFile file(std::string(fileName).c_str()); + if (file.IsZombie()) { + throw std::runtime_error("o2::acts::TrackingGeometryManager: cannot open '" + std::string(fileName) + "'"); + } + return file.GetKey(std::string(o2::base::NameConf::CCDBOBJECT).c_str()) != nullptr || + file.GetKey(std::string(o2::base::NameConf::GEOMOBJECTNAME_FAIR).c_str()) != nullptr; +} +} // namespace + +TrackingGeometryManager::TrackingGeometryManager() + : mNominalContext(Acts::GeometryContext::dangerouslyDefaultConstruct()) +{ +} + +TrackingGeometryManager& TrackingGeometryManager::instance() +{ + static TrackingGeometryManager inst; + return inst; +} + +void TrackingGeometryManager::setBuilder(std::unique_ptr builder) +{ + if (isBuilt()) { + throw std::runtime_error( + "o2::acts::TrackingGeometryManager: cannot change the builder after the geometry was built. " + "Call clear() first if this is really intended."); + } + mBuilder = std::move(builder); +} + +void TrackingGeometryManager::clear() +{ + // Order matters: the surfaces hold non-owning back-pointers into the element + // store, so the geometry must go first. + mGeometry.reset(); + mElementStore.clear(); + mIndex.clear(); + mDetIndices.clear(); + mNDecoratedSurfaces = 0; +} + +void TrackingGeometryManager::build() +{ + if (isBuilt()) { + LOG(info) << "ACTS tracking geometry is already built, nothing to do"; + return; + } + if (!o2::base::GeometryManager::isGeometryLoaded()) { + throw std::runtime_error( + "o2::acts::TrackingGeometryManager: no TGeo geometry is loaded. Declare " + "o2::base::GRPGeomRequest::Aligned in the workflow's GRPGeomRequest and call " + "GRPGeomHelper::instance().checkUpdates(pc) before asking for the ACTS geometry, or use " + "buildFromFile() outside a workflow."); + } + doBuild(); +} + +void TrackingGeometryManager::buildFromFile(std::string_view geomFileName) +{ + if (isBuilt()) { + LOG(info) << "ACTS tracking geometry is already built, nothing to do"; + return; + } + if (o2::base::GeometryManager::isGeometryLoaded()) { + LOG(warn) << "A TGeo geometry is already loaded, ignoring the requested file and using it"; + build(); + return; + } + + LOGP(info, "Loading TGeo geometry from '{}'", geomFileName); + if (hasO2GeometryObject(geomFileName)) { + o2::base::GeometryManager::loadGeometry(geomFileName); + } else { + // Not an o2-sim / CCDB geometry file. GeometryManager::loadGeometry() only + // accepts those two object names and aborts on anything else, so fall back to + // a plain import; this is the path for geometries exported from GDML, which + // the standalone actsO2 chain uses. + LOGP(warn, "'{}' holds no '{}' or '{}' object: importing it as a plain TGeo file", + geomFileName, o2::base::NameConf::CCDBOBJECT, o2::base::NameConf::GEOMOBJECTNAME_FAIR); + if (TGeoManager::Import(std::string(geomFileName).c_str()) == nullptr) { + throw std::runtime_error("o2::acts::TrackingGeometryManager: cannot import TGeo geometry from '" + + std::string(geomFileName) + "'"); + } + } + build(); +} + +void TrackingGeometryManager::doBuild() +{ + if (mBuilder == nullptr) { + throw std::runtime_error( + "o2::acts::TrackingGeometryManager: no builder installed. Call setBuilder() with the " + "detector-specific ITrackingGeometryBuilder before requesting the geometry."); + } + + LOGP(info, "Building the ACTS tracking geometry with the '{}' builder", mBuilder->getName()); + auto output = mBuilder->build(*gGeoManager, mNominalContext); + if (output.geometry == nullptr) { + throw std::runtime_error(std::string("o2::acts::TrackingGeometryManager: builder '") + + mBuilder->getName() + "' returned no geometry"); + } + + mElementStore = std::move(output.elementStore); + mIndex = std::move(output.index); + + if (!mMaterialMapFile.empty()) { + // Must happen while the geometry is still non-const: Gen3 construction takes + // no IMaterialDecorator, so material can only be injected afterwards. + applyMaterialMapImpl(*output.geometry); + } + + mGeometry = std::const_pointer_cast(output.geometry); + + LOGP(info, "ACTS tracking geometry ready: {} mapped sensors, {} surfaces with material", + mIndex.size(), mNDecoratedSurfaces); +} + +void TrackingGeometryManager::applyMaterialMapImpl(Acts::TrackingGeometry& geometry) +{ + if (!std::filesystem::exists(mMaterialMapFile)) { + throw std::runtime_error("o2::acts::TrackingGeometryManager: material map '" + mMaterialMapFile + + "' does not exist"); + } + + Acts::MaterialMapJsonConverter::Config converterConfig; + converterConfig.context = mNominalContext; + const Acts::JsonMaterialDecorator decorator(converterConfig, mMaterialMapFile, mLogLevel); + + mNDecoratedSurfaces = 0; + geometry.apply([&](Acts::Surface& surface) { + decorator.decorate(surface); + if (surface.surfaceMaterial() != nullptr) { + ++mNDecoratedSurfaces; + } + }); + + if (mNDecoratedSurfaces == 0) { + // Material maps are keyed on Acts::GeometryIdentifier, which is assigned during + // construction. Any structural change to the builder shifts those ids and would + // otherwise silently yield a material-free geometry. + throw std::runtime_error( + "o2::acts::TrackingGeometryManager: material map '" + mMaterialMapFile + + "' matched no surface. It was produced for a different geometry structure and must be " + "regenerated for this one."); + } +} + +std::shared_ptr TrackingGeometryManager::getShared() +{ + if (!isBuilt()) { + build(); + } + return mGeometry; +} + +const Acts::TrackingGeometry& TrackingGeometryManager::get() +{ + return *getShared(); +} + +const SurfaceIndexMap& TrackingGeometryManager::getIndex() +{ + if (!isBuilt()) { + build(); + } + return mIndex; +} + +const SurfaceIndexMap& TrackingGeometryManager::getIndex(const o2::detectors::DetMatrixCache& cache, + double toleranceCm, + SensorPathProvider pathProvider) +{ + if (!isBuilt()) { + build(); + } + const int detID = cache.getDetID(); + const auto it = mDetIndices.find(detID); + if (it != mDetIndices.end()) { + return it->second; + } + auto [inserted, ok] = mDetIndices.emplace( + detID, makeSurfaceIndex(*mGeometry, mNominalContext, cache, toleranceCm, std::move(pathProvider))); + return inserted->second; +} + +} // namespace o2::acts diff --git a/Detectors/ACTS/workflow/CMakeLists.txt b/Detectors/ACTS/workflow/CMakeLists.txt new file mode 100644 index 0000000000000..34c9182cf57f4 --- /dev/null +++ b/Detectors/ACTS/workflow/CMakeLists.txt @@ -0,0 +1,15 @@ +# Copyright 2019-2026 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. + +o2_add_library(ACTSWorkflow + SOURCES src/ActsGeometryService.cxx + PUBLIC_LINK_LIBRARIES O2::ACTSInterface + O2::Framework) diff --git a/Detectors/ACTS/workflow/include/ACTSWorkflow/ActsGeometryService.h b/Detectors/ACTS/workflow/include/ACTSWorkflow/ActsGeometryService.h new file mode 100644 index 0000000000000..020fadf2349c8 --- /dev/null +++ b/Detectors/ACTS/workflow/include/ACTSWorkflow/ActsGeometryService.h @@ -0,0 +1,76 @@ +// Copyright 2019-2026 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 ActsGeometryService.h +/// \brief DPL service exposing the ACTS tracking geometry through the ServiceRegistry +/// \author Paolo Butti +/// + +#ifndef ALICEO2_ACTS_ACTSGEOMETRYSERVICE_H +#define ALICEO2_ACTS_ACTSGEOMETRYSERVICE_H + +#include +#include + +#include "Acts/Geometry/TrackingGeometry.hpp" + +#include "ACTSInterface/SurfaceIndexMap.h" +#include "ACTSInterface/TrackingGeometryManager.h" +#include "Framework/ServiceSpec.h" + +namespace o2::acts +{ + +/// Registry-facing handle on the ACTS tracking geometry. +/// +/// This is a thin shim over TrackingGeometryManager: both access paths hand out +/// the very same geometry object, so a workflow can mix them freely. What the +/// service adds is DPL-owned lifetime (the geometry is released on device exit +/// rather than at static destruction) and retrieval through the registry: +/// +/// \code +/// auto& actsGeo = pc.services().get(); +/// const auto& tg = actsGeo.getGeometry(); +/// \endcode +/// +/// Registered with ServiceKind::DeviceGlobal: one instance per device process, +/// shared by every data processor in it, and safe to read from several threads +/// once built. +class ActsGeometryService +{ + public: + /// Builds on first use. Requires gGeoManager to be live and a builder to have + /// been installed on TrackingGeometryManager. + const Acts::TrackingGeometry& getGeometry() const; + std::shared_ptr getGeometryShared() const; + + const SurfaceIndexMap& getIndex() const; + + const Acts::GeometryContext& getNominalContext() const; + + bool isBuilt() const; + + /// Escape hatch for configuration (builder, material map) from init(). + static TrackingGeometryManager& getManager() { return TrackingGeometryManager::instance(); } +}; + +/// ServiceSpec for ActsGeometryService, to be put into +/// DataProcessorSpec::requiredServices. +o2::framework::ServiceSpec actsGeometryServiceSpec(); + +/// CommonServices::defaultServices() plus the ACTS geometry service. Use as the +/// requiredServices of a DataProcessorSpec that wants registry access. +std::vector defaultServicesWithActsGeometry(); + +} // namespace o2::acts + +#endif // ALICEO2_ACTS_ACTSGEOMETRYSERVICE_H diff --git a/Detectors/ACTS/workflow/src/ActsGeometryService.cxx b/Detectors/ACTS/workflow/src/ActsGeometryService.cxx new file mode 100644 index 0000000000000..7cef408a1bd02 --- /dev/null +++ b/Detectors/ACTS/workflow/src/ActsGeometryService.cxx @@ -0,0 +1,79 @@ +// Copyright 2019-2026 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 ActsGeometryService.cxx +/// \author Paolo Butti +/// + +#include "ACTSWorkflow/ActsGeometryService.h" + +#include "Framework/CommonServices.h" +#include "Framework/Logger.h" +#include "Framework/ServiceHandle.h" +#include "Framework/ServiceRegistryRef.h" +#include "Framework/TypeIdHelpers.h" + +using namespace o2::framework; + +namespace o2::acts +{ + +const Acts::TrackingGeometry& ActsGeometryService::getGeometry() const +{ + return TrackingGeometryManager::instance().get(); +} + +std::shared_ptr ActsGeometryService::getGeometryShared() const +{ + return TrackingGeometryManager::instance().getShared(); +} + +const SurfaceIndexMap& ActsGeometryService::getIndex() const +{ + return TrackingGeometryManager::instance().getIndex(); +} + +const Acts::GeometryContext& ActsGeometryService::getNominalContext() const +{ + return TrackingGeometryManager::instance().getNominalContext(); +} + +bool ActsGeometryService::isBuilt() const +{ + return TrackingGeometryManager::instance().isBuilt(); +} + +ServiceSpec actsGeometryServiceSpec() +{ + return ServiceSpec{ + .name = "acts-geometry", + .init = [](ServiceRegistryRef, DeviceState&, fair::mq::ProgOptions&) -> ServiceHandle { + // Nothing is built here: the TGeo geometry is not available yet at service + // init time. The first getGeometry() call, from run(), triggers the build. + return ServiceHandle{TypeIdHelpers::uniqueId(), new ActsGeometryService(), + ServiceKind::DeviceGlobal, "acts-geometry"}; + }, + .configure = CommonServices::noConfiguration(), + .exit = [](ServiceRegistryRef, void* service) { + TrackingGeometryManager::instance().clear(); + delete reinterpret_cast(service); }, + .kind = ServiceKind::DeviceGlobal}; +} + +std::vector defaultServicesWithActsGeometry() +{ + auto services = CommonServices::defaultServices(); + services.push_back(actsGeometryServiceSpec()); + return services; +} + +} // namespace o2::acts diff --git a/Detectors/CMakeLists.txt b/Detectors/CMakeLists.txt index 289fb3fcd0fa7..a3e3902e6076f 100644 --- a/Detectors/CMakeLists.txt +++ b/Detectors/CMakeLists.txt @@ -65,6 +65,8 @@ endif() o2_data_file(COPY Geometry DESTINATION Detectors) +add_subdirectory(ACTS) + if(ENABLE_UPGRADES) add_subdirectory(Upgrades) else() diff --git a/Detectors/Upgrades/ALICE3/ACTS/CMakeLists.txt b/Detectors/Upgrades/ALICE3/ACTS/CMakeLists.txt new file mode 100644 index 0000000000000..763c9db84b174 --- /dev/null +++ b/Detectors/Upgrades/ALICE3/ACTS/CMakeLists.txt @@ -0,0 +1,24 @@ +# Copyright 2019-2026 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. + +# Concrete ALICE 3 builder for the detector-agnostic O2::ACTSInterface. +# Skipped, like that package, when ACTS is not available. +if(NOT Acts_FOUND) + return() +endif() + + +add_subdirectory(interface) +add_subdirectory(macros) + +# # Geometry description consumed at runtime by Gen3BlueprintBuilder. Installed to +# # $O2_ROOT/share/Detectors/Upgrades/ALICE3/ACTS/config/ +# o2_data_file(COPY config DESTINATION Detectors/Upgrades/ALICE3/ACTS) diff --git a/Detectors/Upgrades/ALICE3/ACTS/interface/CMakeLists.txt b/Detectors/Upgrades/ALICE3/ACTS/interface/CMakeLists.txt new file mode 100644 index 0000000000000..b6b3c8c1d046d --- /dev/null +++ b/Detectors/Upgrades/ALICE3/ACTS/interface/CMakeLists.txt @@ -0,0 +1,20 @@ +# 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. + + +o2_add_library(ALICE3ACTS + SOURCES src/TrkTGeoBlueprintBuilder.cxx + PUBLIC_LINK_LIBRARIES O2::ACTSInterface + nlohmann_json::nlohmann_json) + +o2_target_root_dictionary(ALICE3ACTS + HEADERS include/ALICE3ACTS/TrkTGeoBlueprintBuilder.h + ) diff --git a/Detectors/Upgrades/ALICE3/ACTS/interface/include/ALICE3ACTS/TrkTGeoBlueprintBuilder.h b/Detectors/Upgrades/ALICE3/ACTS/interface/include/ALICE3ACTS/TrkTGeoBlueprintBuilder.h new file mode 100644 index 0000000000000..41abe0a27c88d --- /dev/null +++ b/Detectors/Upgrades/ALICE3/ACTS/interface/include/ALICE3ACTS/TrkTGeoBlueprintBuilder.h @@ -0,0 +1,51 @@ +// Copyright 2019-2026 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 TrkTGeoBlueprintBuilder.h +/// \brief Build an ACTS TRK tracking geometry from the TGeo geometry already loaded by O2 +/// + +#ifndef ALICEO2_ALICE3_TRKTGEOBLUEPRINTBUILDER_H +#define ALICEO2_ALICE3_TRKTGEOBLUEPRINTBUILDER_H + +#include + +#include "Acts/Utilities/Logger.hpp" +#include "ACTSInterface/ITrackingGeometryBuilder.h" + +namespace o2::alice3 +{ + +/// First-step automatic ALICE3 builder: TRK sensitive barrel only, no material. +/// The TGeo volume name is a detector-specific selection rule, not a geometry +/// configuration. Positions, shapes, and layer extents come from the live TGeo. +class TrkTGeoBlueprintBuilder : public o2::acts::ITrackingGeometryBuilder +{ + public: + struct Config { + std::string graphvizFile; + Acts::Logging::Level logLevel = Acts::Logging::INFO; + }; + + TrkTGeoBlueprintBuilder(); + explicit TrkTGeoBlueprintBuilder(Config config); + o2::acts::TrackingGeometryOutput build(TGeoManager& tgeo, + const Acts::GeometryContext& gctx) override; + const char* getName() const override { return "alice3-trk-tgeo"; } + + private: + Config mConfig; +}; + +} // namespace o2::alice3 + +#endif // ALICEO2_ALICE3_TRKTGEOBLUEPRINTBUILDER_H diff --git a/Detectors/Upgrades/ALICE3/ACTS/interface/src/ALICE3ACTSLinkDef.h b/Detectors/Upgrades/ALICE3/ACTS/interface/src/ALICE3ACTSLinkDef.h new file mode 100644 index 0000000000000..c85dd2d378ccb --- /dev/null +++ b/Detectors/Upgrades/ALICE3/ACTS/interface/src/ALICE3ACTSLinkDef.h @@ -0,0 +1,18 @@ +// 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. + +#ifdef __CLING__ + +#pragma link off all globals; +#pragma link off all classes; +#pragma link off all functions; + +#endif diff --git a/Detectors/Upgrades/ALICE3/ACTS/interface/src/TrkTGeoBlueprintBuilder.cxx b/Detectors/Upgrades/ALICE3/ACTS/interface/src/TrkTGeoBlueprintBuilder.cxx new file mode 100644 index 0000000000000..873512e276f06 --- /dev/null +++ b/Detectors/Upgrades/ALICE3/ACTS/interface/src/TrkTGeoBlueprintBuilder.cxx @@ -0,0 +1,241 @@ +// Copyright 2019-2026 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 TrkTGeoBlueprintBuilder.cxx +/// \brief Minimal, material-free TRK-only geometry builder using O2's live TGeo +/// + +#include "ALICE3ACTS/TrkTGeoBlueprintBuilder.h" + +#include "Acts/Definitions/Units.hpp" +#include "Acts/Geometry/Blueprint.hpp" +#include "Acts/Geometry/BlueprintOptions.hpp" +#include "Acts/Geometry/ContainerBlueprintNode.hpp" +#include "Acts/Geometry/CylinderVolumeBounds.hpp" +#include "Acts/Geometry/GeometryIdentifierBlueprintNode.hpp" +#include "Acts/Geometry/LayerBlueprintNode.hpp" +#include "Acts/Geometry/TrackingVolume.hpp" +#include "Acts/Geometry/VolumeAttachmentStrategy.hpp" +#include "Acts/Geometry/VolumeResizeStrategy.hpp" +#include "Acts/Surfaces/CylinderBounds.hpp" +#include "Acts/Surfaces/Surface.hpp" +#include "ActsPlugins/Root/TGeoAxes.hpp" +#include "ActsPlugins/Root/TGeoDetectorElement.hpp" +#include "ActsPlugins/Root/TGeoLayerBuilder.hpp" +#include "ActsPlugins/Root/TGeoParser.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "TGeoManager.h" +#include "TGeoMatrix.h" +#include "TGeoNode.h" +#include "TGeoVolume.h" + +namespace o2::alice3 +{ +using namespace Acts; +using namespace Acts::Experimental; +using namespace Acts::UnitLiterals; + +namespace +{ +constexpr double kTGeoToActs = 10.; // cm -> mm +constexpr double kVertexLayerTolerance = 5. * 1_mm; +constexpr double kBarrelLayerTolerance = 15. * 1_mm; +constexpr double kLayerEnvelopeWidth = 2. * 1_mm; + +using ElementPtr = std::shared_ptr; +using KeyedSurface = std::pair>; + +struct LayerGroup { + double radius = 0.; + std::vector> surfaces; +}; + +std::vector clusterLayers(std::vector elements, double tolerance) +{ + std::vector layers; + if (elements.empty()) { + return layers; + } + + std::ranges::sort(elements, {}, &KeyedSurface::first); + LayerGroup group; + double previous = elements.front().first; + double sum = 0.; + + for (const auto& [radius, surface] : elements) { + if (!group.surfaces.empty() && radius - previous > tolerance) { + group.radius = sum / group.surfaces.size(); + layers.push_back(std::move(group)); + group = LayerGroup{}; + sum = 0.; + } + group.surfaces.push_back(surface); + sum += radius; + previous = radius; + } + group.radius = sum / group.surfaces.size(); + layers.push_back(std::move(group)); + return layers; +} + +} // namespace + +TrkTGeoBlueprintBuilder::TrkTGeoBlueprintBuilder() : TrkTGeoBlueprintBuilder(Config{}) {} + +TrkTGeoBlueprintBuilder::TrkTGeoBlueprintBuilder(Config config) : mConfig(std::move(config)) {} + +o2::acts::TrackingGeometryOutput TrkTGeoBlueprintBuilder::build(TGeoManager& tgeo, + const Acts::GeometryContext& gctx) +{ + ACTS_LOCAL_LOGGER(getDefaultLogger("ALICE3TRKAuto", mConfig.logLevel)); + ACTS_INFO("Building ACTS TRK geometry from the live TGeo '" << tgeo.GetName() << "'"); + + auto* topVolume = tgeo.GetTopVolume(); + if (!topVolume) { + throw std::runtime_error("ALICE3 TRK: no top volume in the loaded TGeo geometry"); + } + + // TGeoParser selects the placed TRK sensitive nodes, including the VD + // PETALCASE surfaces. Avoid parsing every support/passive simulation volume. + ActsPlugins::TGeoParser::Options options; + options.volumeNames = {topVolume->GetName()}; + options.targetNames = {"*TRKSensor*"}; + options.unit = kTGeoToActs; + + ActsPlugins::TGeoParser::State state; + state.volume = topVolume; + state.onBranch = true; + TGeoHMatrix transform = TGeoIdentity(topVolume->GetName()); + ActsPlugins::TGeoParser::select(state, options, transform); + + if (state.selectedNodes.empty()) { + throw std::runtime_error("ALICE3 TRK: no TRK sensitive TGeo nodes found (*TRKSensor*)"); + } + + const auto thinY = ActsPlugins::TGeoAxes::parse("ZXY"); + const auto thinZ = ActsPlugins::TGeoAxes::parse("XYZ"); + std::vector elements; + std::vector vertexSurfaces; + std::vector barrelSurfaces; + elements.reserve(state.selectedNodes.size()); + + for (const auto& snode : state.selectedNodes) { + const std::string_view name = snode.node->GetVolume()->GetName(); + const bool vertex = name.find("PETALCASE") != std::string_view::npos; + const auto axes = vertex ? thinZ : thinY; + + auto element = ActsPlugins::TGeoLayerBuilder::defaultElementFactory( + ActsPlugins::TGeoDetectorElement::Identifier{}, *snode.node, + *snode.transform, axes, kTGeoToActs, nullptr); + auto surface = element->surface().getSharedPtr(); + + if (vertex && surface->type() != Surface::SurfaceType::Cylinder) { + throw std::runtime_error("ALICE3 TRK: expected a cylindrical vertex-detector sensor: " + std::string(name)); + } + // TRK is a barrel detector. A misidentified forward surface must not + // silently enter a cylindrical layer and break ACTS navigation. + const auto center = surface->center(gctx); + const double radius = (surface->type() == Surface::SurfaceType::Cylinder) + ? surface->bounds().values()[CylinderBounds::eR] + : std::hypot(center.x(), center.y()); + if (!vertex) { + const auto normal = surface->normal(gctx, center, Vector3::UnitZ()); + if (std::abs(normal.z()) > 0.7) { + throw std::runtime_error("ALICE3 TRK: a non-barrel surface matched *TRKSensor*: " + std::string(name)); + } + } + + (vertex ? vertexSurfaces : barrelSurfaces).emplace_back(radius, std::move(surface)); + elements.push_back(std::move(element)); + } + + auto layers = clusterLayers(std::move(vertexSurfaces), kVertexLayerTolerance); + auto barrelLayers = clusterLayers(std::move(barrelSurfaces), kBarrelLayerTolerance); + for (auto& layer : barrelLayers) { + layers.push_back(std::move(layer)); + } + std::ranges::sort(layers, {}, &LayerGroup::radius); + + if (layers.empty()) { + throw std::runtime_error("ALICE3 TRK: no barrel layers found in the loaded TGeo"); + } + ACTS_INFO("Found " << elements.size() << " sensors in " << layers.size() << " TRK layers"); + for (std::size_t i = 0; i < layers.size(); ++i) { + ACTS_INFO(" TRK layer " << i << ": R=" << layers[i].radius / 1_mm + << " mm, sensors=" << layers[i].surfaces.size()); + } + + Blueprint::Config config; + // Enlarge the inner world radius to include the beamline (r = 0). + const double innerEnvelope = layers.front().radius + 20. * 1_mm; + config.envelope = ExtentEnvelope{{.z = {20. * 1_mm, 20. * 1_mm}, + .r = {innerEnvelope, 20. * 1_mm}}}; + Blueprint root{config}; + auto& ids = root.withGeometryIdentifier(); + const auto byRadius = [](const TrackingVolume& a, const TrackingVolume& b) { + const auto midRadius = [](const TrackingVolume& volume) { + const auto& bounds = dynamic_cast(volume.volumeBounds()); + return 0.5 * (bounds.get(CylinderVolumeBounds::eMinR) + bounds.get(CylinderVolumeBounds::eMaxR)); + }; + return midRadius(a) < midRadius(b); + }; + // One self-contained TRK region. IDs are deterministic within this builder, + // but they are NOT compatible with the full-detector geometry/material map. + ids.setAllVolumeIdsTo(40).incrementLayerIds(1).sortBy(byRadius); + auto& barrel = ids.addCylinderContainer("TRKBarrel", AxisDirection::AxisR); + barrel.setAttachmentStrategy(VolumeAttachmentStrategy::Gap); + barrel.setResizeStrategies(VolumeResizeStrategy::Gap, VolumeResizeStrategy::Gap); + + const ExtentEnvelope layerEnvelope = ExtentEnvelope{{.z = {kLayerEnvelopeWidth, kLayerEnvelopeWidth}, + .r = {kLayerEnvelopeWidth, kLayerEnvelopeWidth}}}; + for (std::size_t i = 0; i < layers.size(); ++i) { + barrel.addLayer("TRKBarrel_L" + std::to_string(i), [&layers, &layerEnvelope, i](LayerBlueprintNode& node) { + node.setSurfaces(layers[i].surfaces); + node.setLayerType(LayerBlueprintNode::LayerType::Cylinder); + node.setEnvelope(layerEnvelope); + node.setUseCenterOfGravity(false, false, true); + }); + } + + if (!mConfig.graphvizFile.empty()) { + std::ofstream out(mConfig.graphvizFile); + if (!out) { + throw std::runtime_error("ALICE3 TRK: cannot write Graphviz file: " + mConfig.graphvizFile); + } + root.graphviz(out); + } + + BlueprintOptions opts; + auto trackingGeometry = root.construct(opts, gctx, logger()); + + o2::acts::TrackingGeometryOutput output; + output.elementStore.reserve(elements.size()); + for (const auto& element : elements) { + output.elementStore.push_back(element); + } + output.geometry = std::shared_ptr(std::move(trackingGeometry)); + ACTS_INFO("ACTS TRK geometry constructed from the loaded TGeo"); + return output; +} + +} // namespace o2::alice3 diff --git a/Detectors/Upgrades/ALICE3/ACTS/macros/CMakeLists.txt b/Detectors/Upgrades/ALICE3/ACTS/macros/CMakeLists.txt new file mode 100644 index 0000000000000..a15c69d108f68 --- /dev/null +++ b/Detectors/Upgrades/ALICE3/ACTS/macros/CMakeLists.txt @@ -0,0 +1,21 @@ +# Copyright 2019-2026 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. + +# Compile-only: running it needs an o2sim geometry file, which does not exist at +# build time. When ACTS is missing this directory is not added at all, so the +# macro is exempted from the "every macro must be tested" check in +# cmake/O2RootMacroExclusionList.cmake. +o2_add_test_root_macro(CheckActsTrackingGeometry.C + PUBLIC_LINK_LIBRARIES O2::ACTSInterface + O2::ALICE3ACTS + O2::TRKBase + LABELS acts upgrades + COMPILE_ONLY) diff --git a/Detectors/Upgrades/ALICE3/ACTS/macros/CheckActsTrackingGeometry.C b/Detectors/Upgrades/ALICE3/ACTS/macros/CheckActsTrackingGeometry.C new file mode 100644 index 0000000000000..5641781e43f03 --- /dev/null +++ b/Detectors/Upgrades/ALICE3/ACTS/macros/CheckActsTrackingGeometry.C @@ -0,0 +1,135 @@ +// Copyright 2019-2026 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 CheckActsTrackingGeometry.C +/// \brief Build the ACTS tracking geometry from a TGeo file and check it against O2 +/// \author Paolo Butti +/// +/// Standalone check, outside any workflow. It +/// 1. loads the TGeo geometry through o2::base::GeometryManager, +/// 2. builds the ACTS tracking geometry with the ALICE 3 Gen3 builder, +/// 3. prints the volume/layer/surface table -- the authority for the +/// (volume, layer) keys the digitisation and seeding JSON files use, +/// 4. locates every O2 TRK chip on an ACTS surface and, for the chips that own +/// their surface, compares the surface centre with the chip origin. That +/// validates the whole TGeo -> ACTS transform chain including the cm -> mm +/// scaling. +/// +/// Usage: +/// \code +/// root -l -q 'CheckActsTrackingGeometry.C+("o2sim_geometry-aligned.root", \ +/// "gen3_geometry_config.json")' +/// \endcode + +#if !defined(__CLING__) || defined(__ROOTCLING__) + +#include +#include +#include + +#include + +#include "Acts/Definitions/Units.hpp" +#include "Acts/Geometry/TrackingGeometry.hpp" +#include "Acts/Surfaces/Surface.hpp" + +#include "ACTSInterface/SurfaceIndexMap.h" +#include "ACTSInterface/TrackingGeometryManager.h" +#include "ALICE3ACTS/Gen3BlueprintBuilder.h" +#include "DetectorsBase/GeometryManager.h" +#include "MathUtils/Utils.h" +#include "TRKBase/GeometryTGeo.h" + +#endif + +int CheckActsTrackingGeometry(const std::string& geomFile = "o2sim_geometry-aligned.root", + const std::string& gen3ConfigFile = "gen3_geometry_config.json", + const std::string& materialMapFile = "", + double toleranceCm = 1e-3) +{ + auto& manager = o2::acts::TrackingGeometryManager::instance(); + + o2::alice3::Gen3BlueprintBuilder::Config builderConfig; + builderConfig.geometryConfigFile = gen3ConfigFile; + builderConfig.withMaterial = true; + manager.setBuilder(std::make_unique(builderConfig)); + if (!materialMapFile.empty()) { + manager.setMaterialMapFile(materialMapFile); + } + + manager.buildFromFile(geomFile); + const auto& geometry = manager.get(); + + // The TRK geometry helper needs its local-to-global matrices for the check below. + auto* trkGeo = o2::trk::GeometryTGeo::Instance(); + trkGeo->fillMatrixCache(o2::math_utils::bit2Mask(o2::math_utils::TransformType::L2G)); + + // Identify the chips by their TGeo node rather than by position: the vertex + // detector petals are tube segments, so all three layers of a petal share both + // origin and rotation and a placement-only match cannot separate them. + const auto& index = manager.getIndex(*trkGeo, toleranceCm, [trkGeo](int chipID) { + return std::string(trkGeo->getMatrixPath(chipID).Data()); + }); + + std::printf("\n--- ACTS <-> O2 sensor check ---\n"); + std::printf("TRK chips : %d\n", trkGeo->getNumberOfChips()); + std::printf("mapped to ACTS surfaces: %zu (on %zu surfaces)\n", index.size(), index.getNSurfaces()); + if (static_cast(index.size()) != trkGeo->getNumberOfChips()) { + std::printf("FAILED: not every chip is mapped\n"); + return 1; + } + + // Compare the ACTS surface centre with the chip origin, for the chips that own + // their surface. Chips sitting on a whole-layer cylinder or disc (the vertex + // detector petals) are skipped: that surface's centre is on the beam axis, not + // at the chip. makeSurfaceIndex() already enforced the tolerance, so this only + // reports the accuracy actually achieved. + double maxResidualCm = 0.; + int worstChip = -1; + int nPerSensor = 0; + int nOnLayerSurface = 0; + for (int chipID = 0; chipID < trkGeo->getNumberOfChips(); ++chipID) { + const auto* surface = index.getSurface(chipID); + if (surface == nullptr) { + std::printf("FAILED: chip %d has no surface\n", chipID); + return 1; + } + if (surface->type() == Acts::Surface::SurfaceType::Cylinder || + surface->type() == Acts::Surface::SurfaceType::Disc) { + ++nOnLayerSurface; + continue; + } + ++nPerSensor; + double rot[9] = {0.}; + double tra[3] = {0.}; + trkGeo->getMatrixL2G(chipID).GetComponents(rot[0], rot[1], rot[2], tra[0], + rot[3], rot[4], rot[5], tra[1], + rot[6], rot[7], rot[8], tra[2]); + const auto centre = surface->center(manager.getNominalContext()); + const double residual = std::hypot(centre[0] / Acts::UnitConstants::cm - tra[0], + centre[1] / Acts::UnitConstants::cm - tra[1], + centre[2] / Acts::UnitConstants::cm - tra[2]); + if (residual > maxResidualCm) { + maxResidualCm = residual; + worstChip = chipID; + } + } + std::printf("chips on their own surface : %d\n", nPerSensor); + std::printf("chips on a layer surface : %d\n", nOnLayerSurface); + std::printf("max |surface centre - chip origin| : %.3e cm (chip %d)\n", maxResidualCm, worstChip); + + std::printf("surfaces with material : %zu\n", manager.getNDecoratedSurfaces()); + std::printf("geometry version : %s\n", + geometry.geometryVersion() == Acts::TrackingGeometry::GeometryVersion::Gen3 ? "Gen3" : "Gen1"); + std::printf("--- OK ---\n"); + return 0; +} diff --git a/Detectors/Upgrades/ALICE3/ACTS/macros/run_test.sh b/Detectors/Upgrades/ALICE3/ACTS/macros/run_test.sh new file mode 100755 index 0000000000000..7a54162e4dd3b --- /dev/null +++ b/Detectors/Upgrades/ALICE3/ACTS/macros/run_test.sh @@ -0,0 +1,84 @@ +#!/bin/bash +# Copyright 2019-2026 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. + +# Generate an ALICE 3 geometry with o2-sim and build the ACTS tracking geometry +# from it, then check every TRK chip against its ACTS surface. +# +# Run it from a scratch directory - o2-sim writes its output into $PWD: +# mkdir -p /tmp/acts && cd /tmp/acts +# $O2_ROOT/share/macro/... or simply: bash /run_test.sh +# +# Requires an environment with O2 and ACTS (e.g. alienv enter O2/latest). + +set -o errexit +set -o pipefail + +nEvents=${nEvents:-1} +generator=${generator:-pythia8hi} +modules=${modules:-"A3IP TRK FT3 TF3"} + +# Detector layouts. This is the configuration the ACTS chain is developed +# against: a fully cylindrical vertex detector (as in the standalone actsO2 +# geometries), simplified-realistic ML/OT, segmented FT3 staves and segmented +# barrel TOFs with the endcap TOFs switched off. +layout="TRKBase.layoutVD=kIRISFullCyl" +layout+=";TRKBase.layoutMLOT=kSimplifiedRealistic" +layout+=";FT3Base.layoutFT3=kSegmentedStave" +layout+=";IOTOFBase.segmentedInnerTOF=true" +layout+=";IOTOFBase.segmentedOuterTOF=true" +layout+=";IOTOFBase.enableBackwardTOF=false" +layout+=";IOTOFBase.enableForwardTOF=false" + +# o2-sim ends by fetching alignment from CCDB, which needs a valid alien token +# and aborts without one - after the geometry file has already been written. +# Skip it: the ACTS geometry is built from the ideal geometry anyway. Set +# SKIP_ALIGNMENT=0 if you do want the alignment step. +if [ "${SKIP_ALIGNMENT:-1}" = "1" ]; then + layout+=";align-geom.mDetectors=none" +fi + +echo "=== o2-sim: generating the ALICE 3 geometry ===" +o2-sim-serial-run5 -n "${nEvents}" -g "${generator}" -m ${modules} \ + --configKeyValues "${layout}" 2>&1 | tee sim_alice3.log + +macroDir=$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd) + +# The Gen3 blueprint builder needs the geometry description that goes with the +# geometry file. The one shipped with this package matches the layout above; set +# GEN3_CONFIG to use a different geometry's gen3_geometry_config.json. +if [ -n "${GEN3_CONFIG}" ]; then + gen3Config=${GEN3_CONFIG} +elif [ -f "${O2_ROOT}/share/Detectors/Upgrades/ALICE3/ACTS/config/gen3_geometry_config.json" ]; then + gen3Config="${O2_ROOT}/share/Detectors/Upgrades/ALICE3/ACTS/config/gen3_geometry_config.json" +elif [ -f "${macroDir}/../config/gen3_geometry_config.json" ]; then + gen3Config="${macroDir}/../config/gen3_geometry_config.json" # running from a source checkout +else + gen3Config=gen3_geometry_config.json +fi +if [ ! -f "${gen3Config}" ]; then + echo "ERROR: no Gen3 geometry config at '${gen3Config}'." >&2 + echo " Set GEN3_CONFIG to the gen3_geometry_config.json for this geometry." >&2 + exit 1 +fi +echo "Using Gen3 geometry config: ${gen3Config}" + +# The macro includes ACTS and O2 headers that are not on ROOT's search path by +# default. +export ROOT_INCLUDE_PATH="${ACTS_ROOT}/include:${EIGEN3_ROOT}/include/eigen3:${O2_ROOT}/include:${ROOT_INCLUDE_PATH}" + +# Compile the macro here rather than in place: ACLiC drops its .so and dictionary +# next to the macro, which would litter the source tree. +cp "${macroDir}/CheckActsTrackingGeometry.C" . + +echo "=== ACTS: building the tracking geometry and checking the sensor index ===" +root.exe -b -q "CheckActsTrackingGeometry.C+(\"o2sim_geometry.root\",\"${gen3Config}\")" \ + 2>&1 | tee CheckActsTrackingGeometry.log diff --git a/Detectors/Upgrades/ALICE3/CMakeLists.txt b/Detectors/Upgrades/ALICE3/CMakeLists.txt index f587772a1885b..ae9ba918dbdce 100644 --- a/Detectors/Upgrades/ALICE3/CMakeLists.txt +++ b/Detectors/Upgrades/ALICE3/CMakeLists.txt @@ -19,4 +19,5 @@ add_subdirectory(AOD) add_subdirectory(IOTOF) add_subdirectory(RICH) add_subdirectory(MID) +add_subdirectory(ACTS) add_subdirectory(macros) diff --git a/Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/CMakeLists.txt b/Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/CMakeLists.txt index 68afd31835999..f59d0cb6093c4 100644 --- a/Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/CMakeLists.txt +++ b/Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/CMakeLists.txt @@ -10,7 +10,10 @@ # or submit itself to any jurisdiction. if(Acts_FOUND) - set(actsTarget Acts::Core) + set(actsTarget + Acts::Core + O2::ACTSInterface + O2::ALICE3ACTS) endif() o2_add_library(ALICE3GlobalReconstruction @@ -38,5 +41,6 @@ o2_add_library(ALICE3GlobalReconstruction TBB::tbb) if(Acts_FOUND) + target_sources(${targetName} PRIVATE src/TrackerACTS.cxx) target_compile_definitions(${targetName} PUBLIC O2_WITH_ACTS) endif() 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..256f5040c5e40 100644 --- a/Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/include/ALICE3GlobalReconstruction/TrackerACTS.h +++ b/Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/include/ALICE3GlobalReconstruction/TrackerACTS.h @@ -19,28 +19,37 @@ #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 Acts 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; @@ -53,7 +62,6 @@ struct TrackerACTSConfig { float maxSeedsPerMiddleSP = 2; float deltaPhiMax = 0.1f; ///< Maximum phi difference for doublets }; - /// Space point representation for tracking struct SpacePoint { float x{0.f}; @@ -62,12 +70,13 @@ 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); } 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}; @@ -80,7 +89,6 @@ struct SeedACTS { const SpacePoint* top{nullptr}; float quality{0.f}; }; - /// TRK Tracker using ACTS algorithms for seeding and track finding template class TrackerACTS @@ -98,6 +106,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(); @@ -109,7 +121,6 @@ class TrackerACTS /// Set the magnetic field strength void setBz(float bz) { mBz = bz; } - /// Get the magnetic field strength float getBz() const { return mBz; } @@ -125,13 +136,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 @@ -147,15 +159,19 @@ class TrackerACTS // 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, @@ -169,7 +185,6 @@ class TrackerACTS "ACTS seeding", "Track finding"}; }; - template template float TrackerACTS::evaluateTask(Func&& task, std::string_view taskName) @@ -179,11 +194,9 @@ float TrackerACTS::evaluateTask(Func&& task, std::string_view taskName) 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/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 8e841000d483c..0d919ff8bebdf 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)) { @@ -136,7 +150,6 @@ std::vector TrackerDPL::createTrackingParamsFromCon if (paramConfig.contains("StartLayerMask")) { params.StartLayerMask = paramConfig["StartLayerMask"].get(); } - if (paramConfig.contains("NSigmaCut")) { params.NSigmaCut = paramConfig["NSigmaCut"].get(); } @@ -155,7 +168,6 @@ std::vector TrackerDPL::createTrackingParamsFromCon if (paramConfig.contains("MaxChi2NDF")) { params.MaxChi2NDF = paramConfig["MaxChi2NDF"].get(); } - if (paramConfig.contains("UseDiamond")) { params.UseDiamond = paramConfig["UseDiamond"].get(); } @@ -220,7 +232,6 @@ std::vector TrackerDPL::createTrackingParamsFromCon if (paramConfig.contains("DropTFUponFailure")) { params.DropTFUponFailure = paramConfig["DropTFUponFailure"].get(); } - if (paramConfig.contains("LayerZ")) { params.LayerZ = paramConfig["LayerZ"].get>(); } @@ -307,6 +318,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); @@ -418,7 +435,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 diff --git a/Detectors/Upgrades/ALICE3/TRKFT3/TRK/macros/test/CheckClusters.C b/Detectors/Upgrades/ALICE3/TRKFT3/TRK/macros/test/CheckClusters.C index a6adf3c6ba6aa..69d898d9f308a 100644 --- a/Detectors/Upgrades/ALICE3/TRKFT3/TRK/macros/test/CheckClusters.C +++ b/Detectors/Upgrades/ALICE3/TRKFT3/TRK/macros/test/CheckClusters.C @@ -75,17 +75,18 @@ void CheckClusters(const std::string& clusfile = "o2clus_trk.root", // ── Chip response (for hit-segment propagation to charge-collection plane) ── // Fetches the same AlpideSimResponse from CCDB as the digitizer (IT3/Calib/APTSResponse) // and computes Y-intersection planes with the same formulas from Digitizer::init() - auto& ccdbMgr = o2::ccdb::BasicCCDBManager::instance(); - ccdbMgr.setURL(ccdbUrl); - if (ccdbTimestamp > 0) { - ccdbMgr.setTimestamp(ccdbTimestamp); - } - auto* alpResp = ccdbMgr.get("IT3/Calib/APTSResponse"); - if (!alpResp) { - LOGP(fatal, "Cannot retrieve AlpideSimResponse from CCDB at {}", ccdbUrl); - return; - } - const float depthMax = alpResp->getDepthMax(); + // auto& ccdbMgr = o2::ccdb::BasicCCDBManager::instance(); + // ccdbMgr.setURL(ccdbUrl); + // if (ccdbTimestamp > 0) { + // ccdbMgr.setTimestamp(ccdbTimestamp); + // } + // auto* alpResp = ccdbMgr.get("IT3/Calib/APTSResponse"); + // if (!alpResp) { + // LOGP(fatal, "Cannot retrieve AlpideSimResponse from CCDB at {}", ccdbUrl); + // return; + // } + // const float depthMax = alpResp->getDepthMax(); + const float depthMax = 500; // ── Y-plane shifts: why VD and ML/OT need different values ──────────────── // diff --git a/cmake/O2RootMacroExclusionList.cmake b/cmake/O2RootMacroExclusionList.cmake index ab79a802a6e3a..5fea628edd22f 100644 --- a/cmake/O2RootMacroExclusionList.cmake +++ b/cmake/O2RootMacroExclusionList.cmake @@ -79,6 +79,13 @@ if(NOT ENABLE_UPGRADES) list(APPEND O2_ROOT_MACRO_EXCLUSION_LIST ${upgradeMacros}) endif() +if(NOT Acts_FOUND) + # ACTS is an optional dependency: Detectors/Upgrades/ALICE3/ACTS is not added + # to the build without it, so its o2_add_test_root_macro is never called + o2_get_list_of_macros(${CMAKE_SOURCE_DIR}/Detectors/Upgrades/ALICE3/ACTS actsMacros) + list(APPEND O2_ROOT_MACRO_EXCLUSION_LIST ${actsMacros}) +endif() + list(REMOVE_DUPLICATES O2_ROOT_MACRO_EXCLUSION_LIST) # check exclusion list contains only existing macros