feat: cosmic track fitter with Acts#4281
Merged
osbornjd merged 14 commits intoJun 6, 2026
Merged
Conversation
…e into cosmics_trk_fitting
Contributor
|
No actionable comments were generated in the recent review. 🎉 ℹ️ Recent review info⚙️ Run configurationConfiguration used: Path: .coderabbit.yaml Review profile: CHILL Plan: Pro Run ID: 📒 Files selected for processing (1)
🚧 Files skipped from review as they are similar to previous changes (1)
📝 WalkthroughWalkthroughRefactors cosmic-track fitting: extracts MaterialSurfaceSelector to a shared header, adds convertSeeds/directNavigator setters and m_directNavigation, precomputes material surfaces optionally, refactors seed kinematics to use signed-radius-sorted cluster positions → circle fit → PCA/momentum helpers, and replaces charge logic with bend-direction counting. ChangesCosmic Track Fitting and Charge Determination
Sequence DiagramsequenceDiagram
participant InitRun
participant Geometry
participant PHActsKalmanFactory
participant PHCosmicsTrkFitter
participant CircleFit
participant PCA
participant Momentum
participant TTree
InitRun->>PHActsKalmanFactory: set logger level, create fitters
InitRun->>Geometry: optionally visit surfaces (m_directNavigation)
Geometry->>PHCosmicsTrkFitter: return m_materialSurfaces
PHCosmicsTrkFitter->>PHCosmicsTrkFitter: loopTracks: build sorted_positions
PHCosmicsTrkFitter->>CircleFit: circleFitByTaubin(sorted_positions)
CircleFit->>PHCosmicsTrkFitter: fit params (R,X0,Y0)
PHCosmicsTrkFitter->>PCA: calculatePCA(seed, sorted_positions)
PCA->>PHCosmicsTrkFitter: pcax/pcay/pcaz
PHCosmicsTrkFitter->>Momentum: calculateMomentum(seed, sorted_positions)
Momentum->>PHCosmicsTrkFitter: px/py/pz
PHCosmicsTrkFitter->>PHCosmicsTrkFitter: getCharge(seed, sorted_positions)
PHCosmicsTrkFitter->>TTree: fill per-seed vectors and m_tree->Fill()
Warning There were issues while running some tools. Please review the errors and either fix the tool's configuration or disable the tool if it's a critical failure. 🔧 Infer (1.2.0)offline/packages/trackreco/PHCosmicsTrkFitter.ccIn file included from offline/packages/trackreco/PHCosmicsTrkFitter.cc:1: ... [truncated 2200 characters] ... from ClangFrontend__CFrontend_errors.protect in file "src/clang/cFrontend_errors.ml", line 48, characters 6-141 Thanks for using CodeRabbit! It's free for OSS, and your support helps us grow. If you like it, consider giving us a shout-out. Comment |
![coderabbitai[bot]](https://avatars.githubusercontent.com/in/347564?s=60&v=4){kind=small}
Contributor
There was a problem hiding this comment.
Actionable comments posted: 6
ℹ️ Review info
⚙️ Run configuration
Configuration used: Path: .coderabbit.yaml
Review profile: CHILL
Plan: Pro
Run ID: a9ba80f1-fffc-499c-97b4-081ccd8fbba1
📒 Files selected for processing (4)
offline/packages/trackbase/ActsTrackFittingAlgorithm.hoffline/packages/trackreco/PHActsTrkFitter.hoffline/packages/trackreco/PHCosmicsTrkFitter.ccoffline/packages/trackreco/PHCosmicsTrkFitter.h
💤 Files with no reviewable changes (1)
- offline/packages/trackreco/PHActsTrkFitter.h
Comment on lines
+30
to
+46
| struct MaterialSurfaceSelector | ||
| { | ||
| std::vector<const Acts::Surface*> surfaces = {}; | ||
|
|
||
| /// @param surface is the test surface | ||
| void operator()(const Acts::Surface* surface) | ||
| { | ||
| if (surface->surfaceMaterial() != nullptr) | ||
| { | ||
| if (std::find(surfaces.begin(), surfaces.end(), surface) == | ||
| surfaces.end()) | ||
| { | ||
| surfaces.push_back(surface); | ||
| } | ||
| } | ||
| } | ||
| }; |
Contributor
There was a problem hiding this comment.
🧩 Analysis chain
🏁 Script executed:
#!/bin/bash
# Confirm <algorithm> is not directly included and std::find is the only dependency on it.
fd -e h -e hpp 'ActsTrackFittingAlgorithm' --exec sh -c '
echo "== {} ==";
grep -nE "`#include` <algorithm>" "{}" || echo "NO direct <algorithm> include";
grep -nE "std::(find|sort|count|copy|transform)\b" "{}"
'Repository: sPHENIX-Collaboration/coresoftware
Length of output: 242
Add #include <algorithm> for std::find in MaterialSurfaceSelector
offline/packages/trackbase/ActsTrackFittingAlgorithm.h uses std::find(...) in MaterialSurfaceSelector::operator() but doesn’t include <algorithm>, so it’s relying on transitive includes.
🛠️ Proposed fix
`#include` <functional>
`#include` <memory>
`#include` <vector>
+#include <algorithm>
Comment on lines
+127
to
+135
| m_fitCfg.dFit = ActsTrackFittingAlgorithm::makeDirectedKalmanFitterFunction( | ||
| m_tGeometry->geometry().tGeometry, | ||
| m_tGeometry->geometry().magField, true, true, 0.0, Acts::FreeToBoundCorrection(), *Acts::getDefaultLogger("DirectedKalman", level)); | ||
|
|
||
| MaterialSurfaceSelector selector; | ||
| if (m_directNavigation) | ||
| { | ||
| m_tGeometry->geometry().tGeometry->visitSurfaces(selector, false); | ||
| m_materialSurfaces = selector.surfaces; |
Contributor
There was a problem hiding this comment.
Wire directNavigator() into the actual fit path.
This initializes m_fitCfg.dFit and caches m_materialSurfaces, but the fitter call later still dispatches (*m_fitCfg.fit)(...) unconditionally. As written, enabling directNavigator() is a no-op, so the new material-aware navigation mode never runs.
Comment on lines
+344
to
+383
| std::vector<Acts::Vector3> pos, sorted_positions; | ||
| // get positions from cluster keys | ||
| // TODO: should implement distortions | ||
| TrackSeedHelper::position_map_t positions; | ||
| for (auto key_iter = tpcseed->begin_cluster_keys(); key_iter != tpcseed->end_cluster_keys(); ++key_iter) | ||
| { | ||
| // get positions from cluster keys | ||
| // TODO: should implement distortions | ||
| TrackSeedHelper::position_map_t positions; | ||
| for( auto key_iter = tpcseed->begin_cluster_keys(); key_iter != tpcseed->end_cluster_keys(); ++key_iter ) | ||
| { | ||
| const auto& key(*key_iter); | ||
| positions.emplace(key, m_tGeometry->getGlobalPosition( key, m_clusterContainer->findCluster(key))); | ||
| } | ||
|
|
||
| TrackSeedHelper::circleFitByTaubin(tpcseed, positions, 0, 58); | ||
| } | ||
|
|
||
| float tpcR = fabs(1. / tpcseed->get_qOverR()); | ||
| float tpcx = tpcseed->get_X0(); | ||
| float tpcy = tpcseed->get_Y0(); | ||
|
|
||
| const auto intersect = | ||
| TrackFitUtils::circle_circle_intersection(m_vertexRadius, | ||
| tpcR, tpcx, tpcy); | ||
| float intx, inty; | ||
|
|
||
| if (std::get<1>(intersect) > std::get<3>(intersect)) | ||
| { | ||
| intx = std::get<0>(intersect); | ||
| inty = std::get<1>(intersect); | ||
| } | ||
| else | ||
| { | ||
| intx = std::get<2>(intersect); | ||
| inty = std::get<3>(intersect); | ||
| } | ||
| std::vector<TrkrDefs::cluskey> keys; | ||
| std::vector<Acts::Vector3> clusPos; | ||
| std::copy(tpcseed->begin_cluster_keys(), tpcseed->end_cluster_keys(), std::back_inserter(keys)); | ||
| TrackFitUtils::getTrackletClusters(m_tGeometry, m_clusterContainer, | ||
| clusPos, keys); | ||
| TrackFitUtils::position_vector_t xypoints, rzpoints; | ||
| for (auto& pos : clusPos) | ||
| { | ||
| float clusr = radius(pos.x(), pos.y()); | ||
| if (pos.y() < 0) | ||
| { | ||
| clusr *= -1; | ||
| } | ||
|
|
||
| // exclude silicon and tpot clusters for now | ||
| if (std::abs(clusr) > 80 || std::abs(clusr) < 30) | ||
| { | ||
| continue; | ||
| } | ||
| xypoints.push_back(std::make_pair(pos.x(), pos.y())); | ||
| rzpoints.push_back(std::make_pair(pos.z(), clusr)); | ||
| const auto& key(*key_iter); | ||
| positions.emplace(key, m_tGeometry->getGlobalPosition(key, m_clusterContainer->findCluster(key))); | ||
| pos.push_back(positions[key]); | ||
| } | ||
| sorted_positions = pos; | ||
|
|
||
| auto rzparams = TrackFitUtils::line_fit(rzpoints); | ||
| float fulllineintz = std::get<1>(rzparams); | ||
| float fulllineslope = std::get<0>(rzparams); | ||
|
|
||
| float slope = tpcseed->get_slope(); | ||
| float intz = m_vertexRadius * slope + tpcseed->get_Z0(); | ||
|
|
||
| Acts::Vector3 inter(intx, inty, intz); | ||
|
|
||
| std::vector<float> tpcparams{tpcR, tpcx, tpcy, tpcseed->get_slope(), | ||
| tpcseed->get_Z0()}; | ||
| auto tangent = TrackFitUtils::get_helix_tangent(tpcparams, | ||
| inter); | ||
|
|
||
| auto tan = tangent.second; | ||
| auto pca = tangent.first; | ||
| std::sort(sorted_positions.begin(), sorted_positions.end(), [](const Acts::Vector3& a, const Acts::Vector3& b) | ||
| { | ||
| float aradius = std::sqrt(a.x()*a.x()+a.y()*a.y()); | ||
| if(a.y() < 0) | ||
| { | ||
| aradius *= -1; | ||
| } | ||
| float bradius = std::sqrt(b.x()*b.x()+b.y()*b.y()); | ||
| if(b.y() < 0) | ||
| { | ||
| bradius *= -1; | ||
| } | ||
| return aradius > bradius; }); | ||
|
|
||
| float p; | ||
| if (m_ConstField) | ||
| { | ||
| p = std::cosh(tpcseed->get_eta()) * fabs(1. / tpcseed->get_qOverR()) * (0.3 / 100) * fieldstrength; | ||
| } | ||
| else | ||
| { | ||
| p = tpcseed->get_p(); | ||
| } | ||
|
|
||
| tan *= p; | ||
| TrackSeedHelper::circleFitByTaubin(tpcseed, positions, 0, 58); | ||
|
|
||
| //! if we got the opposite seed then z will be backwards, so we take the | ||
| //! value of tan.z() multiplied by the sign of the slope determined for | ||
| //! the full cosmic track | ||
| //! same with px/py since a single cosmic produces two seeds that bend | ||
| //! in opposite directions | ||
| float theta = std::atan(fulllineslope); | ||
| /// Normalize to 0<theta<pi | ||
| if (theta < 0) | ||
| { | ||
| theta += M_PI; | ||
| } | ||
| float pz = std::cos(theta) * p; | ||
| if (fulllineslope < 0) | ||
| { | ||
| pz = std::abs(pz); | ||
| } | ||
| else | ||
| { | ||
| pz = std::abs(pz) * -1; | ||
| } | ||
| Acts::Vector3 momentum = Acts::Vector3::Zero(); | ||
| Acts::Vector3 pca = calculatePCA(tpcseed, sorted_positions); | ||
|
|
||
| if (!m_zeroField) | ||
| { | ||
| momentum.x() = charge < 0 ? tan.x() : tan.x() * -1; | ||
| momentum.y() = charge < 0 ? tan.y() : tan.y() * -1; | ||
| } | ||
| else | ||
| { | ||
| auto xyparams = TrackFitUtils::line_fit(xypoints); | ||
| float fulllineslopexy = std::get<0>(xyparams); | ||
| if (fulllineslopexy < 0) | ||
| { | ||
| momentum.x() = fabs(tan.x()); | ||
| } | ||
| else | ||
| { | ||
| momentum.x() = fabs(tan.x()) * -1; | ||
| } | ||
| momentum.y() = fabs(tan.y()) * -1; | ||
| } | ||
| Acts::Vector3 momentum = calculateMomentum(tpcseed, sorted_positions); | ||
|
|
||
| momentum.z() = pz; | ||
| Acts::Vector3 position(pca.x(), pca.y(), | ||
| (m_vertexRadius - fulllineintz) / fulllineslope); | ||
| Acts::Vector3 position = pca * Acts::UnitConstants::cm; | ||
|
|
||
| position *= Acts::UnitConstants::cm; | ||
| if (!is_valid(momentum)) | ||
| { | ||
| continue; | ||
| } | ||
|
|
||
| int charge = getCharge(tpcseed, sorted_positions); |
Contributor
There was a problem hiding this comment.
Guard on TPC cluster count before using sorted_positions.
Line 335 only checks sourceLinks.size(), which counts silicon and TPC measurements together. The downstream circle fit, PCA/momentum helpers, and getCharge use only sorted_positions from the TPC seed, so a seed with fewer than 5 TPC clusters can still pass this gate and then hit undefined behavior on sorted_positions[0..4].
Suggested fix
sorted_positions = pos;
std::sort(sorted_positions.begin(), sorted_positions.end(), [](const Acts::Vector3& a, const Acts::Vector3& b)
{
float aradius = std::sqrt(a.x()*a.x()+a.y()*a.y());
@@
}
return aradius > bradius; });
+
+ if (sorted_positions.size() < 5)
+ {
+ continue;
+ }
TrackSeedHelper::circleFitByTaubin(tpcseed, positions, 0, 58);
Comment on lines
+926
to
+945
| float phi0 = std::atan2(sorted_positions[0].y() - tpcparams[2], sorted_positions[0].x() - tpcparams[1]); | ||
| int posphi = 0; | ||
| int negphi = 0; | ||
| // just take the first 4 outermost clusters as a test to determine the bend angle | ||
| // from the outermost radial cluster | ||
| for (size_t i = 1; i < 5; i++) | ||
| { | ||
| auto key = *clusIter; | ||
| auto cluster = m_clusterContainer->findCluster(key); | ||
| if (!cluster) | ||
| auto cluspos = sorted_positions[i]; | ||
|
|
||
| float phi = std::atan2(cluspos.y() - tpcparams[2], cluspos.x() - tpcparams[1]); | ||
| if (phi > phi0) | ||
| { | ||
| std::cout << "MakeSourceLinks::getCharge: Failed to get cluster with key " << key << " for track seed" << std::endl; | ||
| continue; | ||
| posphi++; | ||
| } | ||
|
|
||
| auto surf = m_tGeometry->maps().getSurface(key, cluster); | ||
| if (!surf) | ||
| else | ||
| { | ||
| continue; | ||
| negphi++; | ||
| } | ||
| } | ||
| int charge = posphi > negphi ? -1 : 1; |
Contributor
There was a problem hiding this comment.
Wrap Δφ before voting on the bend direction.
This compares raw atan2 outputs with phi > phi0. When the outer clusters straddle the -π/π branch cut, a small positive bend becomes a large negative jump and the vote can flip the track charge sign.
Suggested fix
- float phi = std::atan2(cluspos.y() - tpcparams[2], cluspos.x() - tpcparams[1]);
- if (phi > phi0)
+ float phi = std::atan2(cluspos.y() - tpcparams[2], cluspos.x() - tpcparams[1]);
+ float dphi = std::remainder(phi - phi0, 2.f * static_cast<float>(M_PI));
+ if (dphi > 0)
{
posphi++;
}📝 Committable suggestion
‼️ IMPORTANT
Carefully review the code before committing. Ensure that it accurately replaces the highlighted code, contains no missing lines, and has no issues with indentation. Thoroughly test & benchmark the code to ensure it meets the requirements.
Suggested change
| float phi0 = std::atan2(sorted_positions[0].y() - tpcparams[2], sorted_positions[0].x() - tpcparams[1]); | |
| int posphi = 0; | |
| int negphi = 0; | |
| // just take the first 4 outermost clusters as a test to determine the bend angle | |
| // from the outermost radial cluster | |
| for (size_t i = 1; i < 5; i++) | |
| { | |
| auto key = *clusIter; | |
| auto cluster = m_clusterContainer->findCluster(key); | |
| if (!cluster) | |
| auto cluspos = sorted_positions[i]; | |
| float phi = std::atan2(cluspos.y() - tpcparams[2], cluspos.x() - tpcparams[1]); | |
| if (phi > phi0) | |
| { | |
| std::cout << "MakeSourceLinks::getCharge: Failed to get cluster with key " << key << " for track seed" << std::endl; | |
| continue; | |
| posphi++; | |
| } | |
| auto surf = m_tGeometry->maps().getSurface(key, cluster); | |
| if (!surf) | |
| else | |
| { | |
| continue; | |
| negphi++; | |
| } | |
| } | |
| int charge = posphi > negphi ? -1 : 1; | |
| float phi0 = std::atan2(sorted_positions[0].y() - tpcparams[2], sorted_positions[0].x() - tpcparams[1]); | |
| int posphi = 0; | |
| int negphi = 0; | |
| // just take the first 4 outermost clusters as a test to determine the bend angle | |
| // from the outermost radial cluster | |
| for (size_t i = 1; i < 5; i++) | |
| { | |
| auto cluspos = sorted_positions[i]; | |
| float phi = std::atan2(cluspos.y() - tpcparams[2], cluspos.x() - tpcparams[1]); | |
| float dphi = std::remainder(phi - phi0, 2.f * static_cast<float>(M_PI)); | |
| if (dphi > 0) | |
| { | |
| posphi++; | |
| } | |
| else | |
| { | |
| negphi++; | |
| } | |
| } | |
| int charge = posphi > negphi ? -1 : 1; |
Comment on lines
+967
to
+988
| auto arcLength = [&](float x, float y) | ||
| { | ||
| float angle = std::atan2(y - tpcy, x - tpcx); | ||
| return tpcR * angle; | ||
| }; | ||
|
|
||
| float sum_s = 0, sum_z = 0, sum_ss = 0, sum_sz = 0; | ||
| int n = sorted_positions.size(); | ||
| // Compute the arc-length parameter for each cluster, then fit to a line | ||
| // Fit z = a + b*s using simple linear regression | ||
| for (auto& p : sorted_positions) | ||
| { | ||
| float s = arcLength(p.x(), p.y()); | ||
| sum_s += s; | ||
| sum_z += p.z(); | ||
| sum_ss += s * s; | ||
| sum_sz += s * p.z(); | ||
| } | ||
|
|
||
| Acts::Vector3 globalMostOuter(std::numeric_limits<double>::quiet_NaN(), std::numeric_limits<double>::quiet_NaN(), std::numeric_limits<double>::quiet_NaN()); | ||
| Acts::Vector3 globalSecondMostOuter(0, 999999, 0); | ||
| float largestR = 0; | ||
| // loop over global positions | ||
| for (auto& i : global_vec) | ||
| { | ||
| Acts::Vector3 global = i; | ||
| // float r = std::sqrt(square(global.x()) + square(global.y())); | ||
| float r = radius(global.x(), global.y()); | ||
| float denom = n * sum_ss - sum_s * sum_s; | ||
| float b = (n * sum_sz - sum_s * sum_z) / denom; | ||
| float a = (sum_z - b * sum_s) / n; |
Contributor
There was a problem hiding this comment.
Unwrap the circle angle before fitting z(s).
Here s is built from tpcR * atan2(...), which is discontinuous at ±π. Tracks whose clusters cross that cut inject an artificial ~2πR jump into the regression, so the fitted slope and PCA z can be badly wrong even when the XY circle fit is fine.
Comment on lines
+1018
to
+1036
| TrackFitUtils::position_vector_t xypoints, rzpoints; | ||
| for (auto& p : sorted_positions) | ||
| { | ||
| if (i.y() < 0) | ||
| float clusr = radius(p.x(), p.y()); | ||
| if (p.y() < 0) | ||
| { | ||
| continue; | ||
| clusr *= -1; | ||
| } | ||
|
|
||
| float dr = std::sqrt(square(globalMostOuter.x()) + square(globalMostOuter.y())) - std::sqrt(square(i.x()) + square(i.y())); | ||
| //! Place a dr cut to get maximum bend due to TPC clusters having | ||
| //! larger fluctuations | ||
| if (dr < maxdr && dr > 10) | ||
| // exclude silicon and tpot clusters for now | ||
| if (std::abs(clusr) > 80 || std::abs(clusr) < 30) | ||
| { | ||
| maxdr = dr; | ||
| globalSecondMostOuter = i; | ||
| continue; | ||
| } | ||
| xypoints.push_back(std::make_pair(p.x(), p.y())); | ||
| rzpoints.push_back(std::make_pair(p.z(), clusr)); | ||
| } | ||
|
|
||
| //! we have to calculate phi WRT the vertex position outside the detector, | ||
| //! not at (0,0) | ||
| Acts::Vector3 vertex(0, m_vertexRadius, 0); | ||
| globalMostOuter -= vertex; | ||
| globalSecondMostOuter -= vertex; | ||
| auto rzparams = TrackFitUtils::line_fit(rzpoints); |
Contributor
There was a problem hiding this comment.
Validate the post-filter sample size before calling line_fit.
After the 30 < |r| < 80 cut, the new outside-volume cosmics can easily leave fewer than two TPC points. line_fit(rzpoints)—and in zero field, line_fit(xypoints)—then has no well-defined slope, so fulllineslope/pz become garbage even though momentum may still look finite.
Also applies to: 1095-1096
Contributor
Author
|
Here's an example momentum comparison to the truth momentum |
Build & test reportReport for commit 6cb0dd2f91a4d3c5c631aed6c00463b513d37a03:
Automatically generated by sPHENIX Jenkins continuous integration |
Build & test reportReport for commit c4441c0c1294542ae4e5cacf2afc00e3039bfb18:
Automatically generated by sPHENIX Jenkins continuous integration |
Build & test reportReport for commit 3891ae7f5ffad4cfca208eee7f6f4ab9c8829728:
Automatically generated by sPHENIX Jenkins continuous integration |
Build & test reportReport for commit 8201aee906706af3d5bdb3e2027f824bf3cd319e:
Automatically generated by sPHENIX Jenkins continuous integration |
Build & test reportReport for commit 4f61b7334c165724eb5bae2963166cb1e66a1039:
Automatically generated by sPHENIX Jenkins continuous integration |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Overhaul of the cosmic track fitter with Acts.
Fits single muons that originate outside of the nominal tracking volume now in simulation, next to test on data
Types of changes
What kind of change does this PR introduce? (Bug fix, feature, ...)
TODOs (if applicable)
Links to other PRs in macros and calibration repositories (if applicable)
Cosmic Track Fitter Improvements Using Acts Framework
Motivation / context
Cosmic muons that originate outside the nominal tracking volume need different seeding and fitting logic than collision tracks. This PR overhauls the Acts-based cosmic track fitter to improve seed PCA, seed momentum/charge estimation, direct-navigation/material handling, and diagnostics to better reconstruct single-muon cosmic tracks in simulation (data validation planned next).
Key changes
Potential risk areas
Possible future improvements
Note: I used the repository headers/sources to produce this summary; AI can make mistakes or miss subtleties—please review the modified algorithms (charge, PCA, momentum), ROOT tree branch changes, and navigation/material preselection behavior carefully before merging.