Add ACTS tracker

UpdateClsACTS

U

Author: njacazio

Detectors/Upgrades/ALICE3/GlobalReconstruction/reconstruction/src/TrackerACTS.cxx

@@ -28,6 +28,8 @@
 #include <Acts/Seeding/detail/CylindricalSpacePointGrid.hpp>
 #include <Acts/Utilities/GridBinFinder.hpp>
 #include <Acts/Utilities/RangeXD.hpp>
+#include <Acts/Definitions/TrackParametrization.hpp>
+#include <Acts/Seeding/EstimateTrackParamsFromSeed.hpp>
 
 namespace o2::trk
 {
@@ -238,17 +240,143 @@ void TrackerACTS<nLayers>::createSeeds()
 template <int nLayers>
 bool TrackerACTS<nLayers>::estimateTrackParams(const SeedACTS& seed, o2::its::TrackITSExt& track) const
 {
+  const SpacePoint* sp0 = seed.bottom;
+  const SpacePoint* sp1 = seed.middle;
+  const SpacePoint* sp2 = seed.top;
+
+  // Use ACTS parameter estimation
+  Acts::Vector3 pos0{sp0->x, sp0->y, sp0->z};
+  Acts::Vector3 pos1{sp1->x, sp1->y, sp1->z};
+  Acts::Vector3 pos2{sp2->x, sp2->y, sp2->z};
+
+  // Magnetic field vector (along z-axis)
+  Acts::Vector3 bField{0., 0., mBz * Acts::UnitConstants::T};
+
+  // Use the ACTS function with time parameter (t0 = 0)
+  LOG(info) << "Calling ACTS estimateTrackParamsFromSeed with mag field " << mBz << " T";
+  LOG(info) << "Seed space points: (" << pos0.transpose() << "), (" << pos1.transpose() << "), (" << pos2.transpose() << ")";
+
+  Acts::FreeVector params;
+  try {
+    params = Acts::estimateTrackParamsFromSeed(pos0, 0.0, pos1, pos2, bField);
+  } catch (const std::exception& e) {
+    LOG(fatal) << "ACTS parameter estimation failed: " << e.what();
+    return false;
+  }
+  LOG(info) << "ACTS parameter estimation successful: x=" << params[Acts::eFreePos0] << " y=" << params[Acts::eFreePos1]
+            << " z=" << params[Acts::eFreePos2] << " q/p=" << params[Acts::eFreeQOverP];
+
+  // Extract parameters from ACTS format
+  const auto& p = params;
+  // ACTS FreeVector: x, y, z, t, dir_x, dir_y, dir_z, q/|p|
+  // Direction components are normalized (unit vector)
+  float x = p[Acts::eFreePos0];
+  float y = p[Acts::eFreePos1];
+  float z = p[Acts::eFreePos2];
+  float qOverP = p[Acts::eFreeQOverP];
+  float pMag = 1.0f / std::abs(qOverP); // |p| = 1 / |q/p| (for unit charge)
+  float px = p[Acts::eFreeDir0] * pMag;
+  float py = p[Acts::eFreeDir1] * pMag;
+  float pz = p[Acts::eFreeDir2] * pMag;
+
+  // Calculate track parameters in O2 format
+  float pt = std::hypot(px, py);
+  float phi = std::atan2(py, px);
+  float theta = std::atan2(pt, pz);
+  float eta = -std::log(std::tan(theta / 2.0f));
+
+  // Set cluster indices from seed
+  track.setExternalClusterIndex(sp0->layer, sp0->clusterId);
+  track.setExternalClusterIndex(sp1->layer, sp1->clusterId);
+  track.setExternalClusterIndex(sp2->layer, sp2->clusterId);
+
+  // Set track parameters
+  track.setX(x);
+  track.setAlpha(phi);
+  track.setY(y * std::cos(phi) - x * std::sin(phi));
+  track.setZ(z);
+  track.setSnp(std::sin(phi - track.getAlpha()));
+  track.setTgl(std::tan(o2::constants::math::PIHalf - theta));
+
+  // q/pT = q/|p| * |p|/pT = qOverP * |p| / pT = qOverP / (pT / |p|) = qOverP / sin(theta)
+  // Or simply: charge / pT where charge = sign(qOverP)
+  float charge = (qOverP > 0) ? 1.0f : -1.0f;
+  track.setQ2Pt(charge / pt);
+
+  LOG(info) << "Estimated track parameters";
   return true;
 }
 
 template <int nLayers>
 void TrackerACTS<nLayers>::findTracks()
 {
+  return; // For now we only create seeds, track finding and fitting will be implemented in the next iterations
+  int nTracks = 0;
+
+  for (const auto& seed : mSeeds) {
+    o2::its::TrackITSExt track;
+
+    LOG(info) << "Estimating track parameters for seed with quality (pT) = " << seed.quality;
+    if (!estimateTrackParams(seed, track)) {
+      continue;
+    }
+
+    // Add track to TimeFrame
+    const int rof = seed.middle->rof;
+    if (mTimeFrame && rof >= 0 && rof < mTimeFrame->getNrof(0)) {
+      LOG(info) << "Adding track to ROF " << rof;
+      auto& tracks = mTimeFrame->getTracks();
+      // tracks.emplace_back(track);
+      ++nTracks;
+    }
+  }
+  LOG(info) << "Created " << nTracks << " tracks from " << mSeeds.size() << " seeds";
 }
 
 template <int nLayers>
 void TrackerACTS<nLayers>::computeTracksMClabels()
 {
+  return; // For now we skip MC labeling, will be implemented in the next iterations once we have track candidates to label
+  if (!mTimeFrame || !mTimeFrame->hasMCinformation()) {
+    return;
+  }
+
+  // MC labeling using majority voting on cluster labels
+  for (int iROF = 0; iROF < mTimeFrame->getNrof(0); ++iROF) {
+    for (auto& track : mTimeFrame->getTracks()) {
+      std::vector<std::pair<MCCompLabel, size_t>> labelCounts;
+
+      for (int iCluster = 0; iCluster < o2::its::TrackITSExt::MaxClusters; ++iCluster) {
+        const int clusterIdx = track.getClusterIndex(iCluster);
+        if (clusterIdx == o2::its::constants::UnusedIndex) {
+          continue;
+        }
+
+        auto clusterLabels = mTimeFrame->getClusterLabels(iCluster, clusterIdx);
+        for (const auto& label : clusterLabels) {
+          auto it = std::find_if(labelCounts.begin(), labelCounts.end(),
+                                 [&label](const auto& p) { return p.first == label; });
+          if (it != labelCounts.end()) {
+            ++(it->second);
+          } else {
+            labelCounts.emplace_back(label, 1);
+          }
+        }
+      }
+
+      if (!labelCounts.empty()) {
+        // Find label with most occurrences
+        auto maxIt = std::max_element(labelCounts.begin(), labelCounts.end(),
+                                      [](const auto& a, const auto& b) { return a.second < b.second; });
+
+        MCCompLabel trackLabel = maxIt->first;
+        if (maxIt->second < static_cast<size_t>(track.getNumberOfClusters())) {
+          trackLabel.setFakeFlag();
+        }
+        mTimeFrame->getTracksLabel().emplace_back(trackLabel);
+      }
+    }
+  }
 }
 
 template <int nLayers>

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<o2::itsmft::AlpideSimResponse>("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<o2::itsmft::AlpideSimResponse>("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 ────────────────
   //

Detectors/Upgrades/ALICE3/TRKFT3/TRK/simulation/src/TRKLayer.cxx

@@ -91,7 +91,8 @@ void TRKCylindricalLayer::createLayer(TGeoVolume* motherVolume)
 TRKSegmentedLayer::TRKSegmentedLayer(int layerNumber, std::string layerName, float rInn, float tiltAngle, int numberOfStaves, int numberOfModules, float thickOrX2X0, MatBudgetParamMode mode)
   : TRKCylindricalLayer(layerNumber, layerName, rInn, numberOfModules * sModuleLength, thickOrX2X0, mode), mTiltAngle(tiltAngle), mNumberOfStaves(numberOfStaves), mNumberOfModules(numberOfModules)
 {
-  assert(numberOfStaves % 2 == 0 && "Error: numberOfStaves must be even!");
+  LOG(info) << "Creating segmented layer: id: " << mLayerNumber << " rInner: " << mInnerRadius << " rOuter: " << mOuterRadius << " zLength: " << mLength
+         << " x2X0: " << mX2X0 << " number of modules: " << mNumberOfModules;
 }
 
 TGeoVolume* TRKSegmentedLayer::createSensor()

Detectors/Upgrades/ALICE3/TRKFT3/common/reconstruction/include/TRKReconstruction/Clusterer.h

@@ -176,6 +176,9 @@ class Clusterer
   static o2::math_utils::Point3D<float> getClusterLocalCoordinates(const ClusterType& cluster, const uint8_t* patt,
                                                                    float yPlaneMLOT = 0.f) noexcept;
 
+  static o2::math_utils::Point3D<float> getClusterLocalCoordinates(const Cluster& cluster, const uint8_t* patt,
+                                                                    float yPlaneMLOT = 0.f) noexcept;
+
  protected:
   int mNHugeClus = 0;
   std::unique_ptr<ClustererThread> mThread;