27 #include "Acts/TrackFinding/TrackStateCreator.hpp"
28 #include "Acts/Surfaces/PlaneSurface.hpp"
29 #include "Acts/Surfaces/RectangleBounds.hpp"
30 #include "Acts/Utilities/VectorHelpers.hpp"
35 #include "GaudiKernel/PhysicalConstants.h"
80 auto magneticField = std::make_unique<ATLASMagneticFieldWrapper>();
81 std::shared_ptr<const Acts::TrackingGeometry> trackingGeometry =
m_trackingGeometryTool->trackingGeometry();
84 detail::Navigator::Config
cfg{trackingGeometry};
85 cfg.resolvePassive =
true;
86 cfg.resolveMaterial =
true;
87 cfg.resolveSensitive =
true;
95 Acts::TrackSelector::EtaBinnedConfig trackSelectorCfg(std::vector<double>({0, 4}));
96 trackSelectorCfg.cutSets[0].ptMin = 500;
97 trackSelectorCfg.cutSets[0].ptMax = 1000000;
98 trackSelectorCfg.cutSets[0].minMeasurements = 1;
99 trackSelectorCfg.cutSets[0].maxHoles = 4;
100 trackSelectorCfg.cutSets[0].maxOutliers = 4;
101 trackSelectorCfg.cutSets[0].maxHolesAndOutliers = 4;
102 trackSelectorCfg.cutSets[0].maxSharedHits = 4;
103 trackSelectorCfg.cutSets[0].maxChi2 = 10000000.0;
108 std::move(extrapolator),
111 Acts::TrackSelector{trackSelectorCfg}};
113 m_trackFinder = std::make_unique<detail::CKF_config>(std::move(ckfConfig));
115 return StatusCode::SUCCESS;
135 ATH_MSG_DEBUG(
"Size of trackParticles collection " << trackParticles->size());
156 ATH_MSG_DEBUG(
"Retrieved " << uncalibratedMeasurementContainer->size()
172 std::unordered_set<uint32_t> hgtdVolumes {};
173 std::unordered_map<uint32_t, std::unordered_set<uint32_t>> hgtdLayers {};
179 Acts::GeometryIdentifier geoID = surface->geometryId();
183 hgtdVolumes.insert(vol);
184 hgtdLayers[vol].insert(
layer);
188 ATH_MSG_ERROR(
"Failed to retrieve surface for cluster " << cluster->index());
189 return StatusCode::FAILURE;
193 ATH_MSG_DEBUG(
"Found " << hgtdVolumes.size() <<
" different HGTD volumes");
195 ATH_MSG_DEBUG(
"HGTD Volume " << vol <<
" has " << hgtdLayers[vol].
size() <<
" layers:");
205 using DefaultTrackStateCreator = Acts::TrackStateCreator<ActsTrk::detail::UncalibSourceLinkAccessor::Iterator,detail::RecoTrackContainer>;
208 DefaultTrackStateCreator::SourceLinkAccessor slAccessorDelegate;
213 ATH_MSG_DEBUG(
"Create " << uncalibratedMeasurementContainer->size()
214 <<
" source links from measurements in "
219 {uncalibratedMeasurementContainer},
223 Acts::PropagatorPlainOptions plainOptions(geoContext, mfContext);
224 plainOptions.direction = Acts::Direction::Forward();
239 defaultTrackStateCreator.sourceLinkAccessor = slAccessorDelegate;
240 defaultTrackStateCreator.calibrator.template connect<&detail::OnTrackCalibrator<detail::RecoTrackStateContainer>::calibrate>(&calibrator);
243 &DefaultTrackStateCreator::createTrackStates>(&defaultTrackStateCreator);
255 if (!link_to_track.
isValid()) {
256 ATH_MSG_ERROR(
"Invalid ACTS track link for TrackParticle " << trackParticle->index());
257 return StatusCode::FAILURE;
260 std::optional<ActsTrk::TrackContainer::ConstTrackProxy> optional_track = *link_to_track;
261 if (!optional_track.has_value()) {
262 ATH_MSG_ERROR(
"No valid ACTS track associated with TrackParticle " << trackParticle->index());
263 return StatusCode::FAILURE;
266 const ActsTrk::TrackContainer::ConstTrackProxy&
track = optional_track.value();
268 Acts::VectorTrackContainer trackBackend;
269 Acts::VectorMultiTrajectory trackStateBackend;
280 <<
"], skipping extension ------ !!!!!");
283 layerHasExtensionHandle(*trackParticle) = trackData.
hasClusterVec;
284 layerExtensionChi2Handle(*trackParticle) = trackData.
chi2Vec;
285 layerClusterRawTimeHandle(*trackParticle) = trackData.
rawTimeVec;
286 layerClusterTimeHandle(*trackParticle) = trackData.
timeVec;
287 extrapXHandle(*trackParticle) = trackData.
extrapX;
288 extrapYHandle(*trackParticle) = trackData.
extrapY;
289 numHGTDHitsHandle(*trackParticle) = trackData.
numHGTDHits;
294 float trackpT =
track.transverseMomentum();
296 float trackNmeasurements =
track.nMeasurements();
298 ATH_MSG_DEBUG(
"TrackParticle " << trackParticle->index() <<
" has ACTS track with eta: " <<
trackEta <<
", phi = " <<
trackPhi <<
" pT: " << trackpT <<
" and nMeasurements: " << trackNmeasurements);
301 const auto lastMeasurementState = Acts::findLastMeasurementState(
track);
302 if (not lastMeasurementState.ok()) {
303 ATH_MSG_ERROR(
"Problem finding last measurement state for acts track");
304 return StatusCode::FAILURE;
306 const Acts::BoundTrackParameters lastMeasurementStateParameters =
track.createParametersFromState(*lastMeasurementState);
309 const Acts::Surface& refSurface =
track.referenceSurface();
310 const Acts::BoundTrackParameters parametersAtRefSurface(refSurface.getSharedPtr(),
313 track.particleHypothesis());
315 ATH_MSG_DEBUG(
"Initial track parameters for extension - lastMeasurementStateParameters:");
317 ATH_MSG_DEBUG(
" - phi: " << lastMeasurementStateParameters.phi());
318 ATH_MSG_DEBUG(
" - pT: " << std::abs(1./lastMeasurementStateParameters.qOverP() *
std::sin(lastMeasurementStateParameters.theta())));
319 ATH_MSG_DEBUG(
" - theta: " << lastMeasurementStateParameters.theta());
320 ATH_MSG_DEBUG(
" - qOverP: " << lastMeasurementStateParameters.qOverP());
321 ATH_MSG_DEBUG(
" - covariance exists: " << (lastMeasurementStateParameters.covariance().has_value() ?
"yes" :
"no"));
326 ATH_MSG_DEBUG(
"Built " << tracksContainerTemp.size() <<
" tracks from it");
328 for (
const detail::RecoTrackContainer::TrackProxy trackProxy : tracksContainerTemp) {
333 layerHasExtensionHandle(*trackParticle) = trackData.
hasClusterVec;
334 layerExtensionChi2Handle(*trackParticle) = trackData.
chi2Vec;
335 layerClusterRawTimeHandle(*trackParticle) = trackData.
rawTimeVec;
336 layerClusterTimeHandle(*trackParticle) = trackData.
timeVec;
337 extrapXHandle(*trackParticle) = trackData.
extrapX;
338 extrapYHandle(*trackParticle) = trackData.
extrapY;
339 numHGTDHitsHandle(*trackParticle) = trackData.
numHGTDHits;
343 return StatusCode::SUCCESS;
366 if (not surface)
continue;
371 auto localPosition = cluster->localPosition<3>();
372 double clusterTime = localPosition[2];
381 mon_cluster_x, mon_cluster_y, mon_cluster_z,
390 return StatusCode::SUCCESS;
394 const EventContext& ctx,
396 const detail::RecoTrackContainer::TrackProxy& trackProxy)
const {
408 std::size_t nMeasurements = 0;
409 std::size_t nHoles = 0;
410 std::size_t nOutliers = 0;
411 std::size_t nHGTDHits = 0;
413 std::vector<bool> hasHitInLayer = {
false,
false,
false,
false};
414 std::vector<float> chi2PerLayer = {0.0, 0.0, 0.0, 0.0};
415 std::vector<float> timePerLayer = {0.0, 0.0, 0.0, 0.0};
416 std::vector<float> rawTimePerLayer = {0.0, 0.0, 0.0, 0.0};
418 std::vector<int> truthClassPerLayer = {-1, -1, -1, -1};
419 std::vector<bool> isShadowedPerLayer = {
false,
false,
false,
false};
420 std::vector<bool> isMergedPerLayer = {
false,
false,
false,
false};
421 std::vector<bool> primaryExpectedPerLayer = {
false,
false,
false,
false};
427 bool foundExtrapolation =
false;
429 trackProxy.container().trackStateContainer().visitBackwards(
430 trackProxy.tipIndex(),
431 [&](
const auto& state) {
432 auto flags = state.typeFlags();
433 if (flags.test(Acts::TrackStateFlag::HoleFlag)) {
435 }
else if (
flags.test(Acts::TrackStateFlag::OutlierFlag)) {
437 }
else if (
flags.test(Acts::TrackStateFlag::MeasurementFlag)) {
441 if (state.hasReferenceSurface() && flags.test(Acts::TrackStateFlag::MeasurementFlag)) {
442 const auto& surface = state.referenceSurface();
443 Acts::GeometryIdentifier geoID = surface.geometryId();
445 if (isHGTDSurface(geoID)) {
446 std::size_t layerIndex = getHGTDLayerIndex(geoID);
448 const auto& calibrated = state.template calibrated<3>();
449 const auto& predicted = state.predicted();
450 const auto& calibCov = state.template calibratedCovariance<3>();
454 Eigen::Vector2d residual2d;
455 residual2d(0) = calibrated(0) - predicted(Acts::eBoundLoc0);
456 residual2d(1) = calibrated(1) - predicted(Acts::eBoundLoc1);
459 AmgSymMatrix(2) cov_2d{calibCov.template block<2,2>(0,0)};
462 const auto& predictedCov = state.predictedCovariance();
463 AmgSymMatrix(2) predicted_cov_2d{predictedCov.template block<2,2>(0,0)};
466 AmgSymMatrix(2) residual_cov = cov_2d + predicted_cov_2d;
472 if (residual_cov.determinant() != 0)
474 chi2 = residual2d.transpose() * residual_cov.inverse() * residual2d;
481 if (layerIndex < 4) {
483 hasHitInLayer[layerIndex] = true;
484 chi2PerLayer[layerIndex] =chi2/ndf;
488 float rawTime = 0.0f;
489 float calibratedTime = 0.0f;
491 if (state.hasCalibrated()) {
494 const auto& calibrated = state.template calibrated<3>();
495 calibratedTime = static_cast<float>(calibrated(2));
496 ATH_MSG_DEBUG(
"Got time from calibrated<3>: " << calibratedTime);
497 } catch (const std::exception& e) {
498 ATH_MSG_WARNING(
"Failed to extract time from calibrated<3>: " << e.what());
504 const xAOD::HGTDCluster* cluster = getHGTDClusterFromState(state);
507 auto localPos = cluster->localPosition<3>();
508 rawTime = static_cast<float>(localPos[2]);
509 ATH_MSG_DEBUG(
"Got raw time from cluster local position: " << rawTime);
511 ATH_MSG_WARNING(
"Could not get cluster from state");
515 rawTimePerLayer[layerIndex] = rawTime;
519 auto [correctedTime, timeErr] = correctTOF(
524 acts_tracking_geometry,
526 timePerLayer[layerIndex] = correctedTime;
527 ATH_MSG_DEBUG(
"Applied TOF correction: " << calibratedTime <<
" -> " << correctedTime);
530 timePerLayer[layerIndex] = calibratedTime;
531 ATH_MSG_DEBUG(
"No cluster found for TOF correction, using calibrated time: " << calibratedTime);
536 ATH_MSG_DEBUG(
"State does not have calibrated data");
540 if (!foundExtrapolation) {
541 foundExtrapolation = true;
542 if (state.hasPredicted()) {
544 const auto& predicted = state.predicted();
545 Acts::Vector2 localPos(predicted[Acts::eBoundLoc0], predicted[Acts::eBoundLoc1]);
548 Acts::Vector3 globalPos = surface.localToGlobal(
551 Acts::Vector3::Zero());
553 extrapX = globalPos.x();
554 extrapY = globalPos.y();
555 extrapZ = globalPos.z();
557 ATH_MSG_DEBUG(
"Extrapolated position (predicted) at HGTD: x=" << extrapX
558 <<
", y=" << extrapY <<
", z=" << extrapZ);
561 Acts::Vector3 globalPos = surface.center(geoContext);
562 extrapX = globalPos.x();
563 extrapY = globalPos.y();
564 extrapZ = globalPos.z();
566 ATH_MSG_DEBUG(
"Extrapolated position (surface center) at HGTD: x=" << extrapX
567 <<
", y=" << extrapY <<
", z=" << extrapZ);
571 ATH_MSG_DEBUG(
"Found HGTD hit on layer " << layerIndex
572 <<
", chi2=" << chi2PerLayer[layerIndex]
573 <<
", time=" << timePerLayer[layerIndex]);
584 if (nHGTDHits == 0 && !foundExtrapolation) {
587 if (getExtrapolationPosition(ctx, trackProxy, extrapX, extrapY, extrapZ)) {
588 ATH_MSG_DEBUG(
"!!! Didn't find any HGTD hits but calculated extrapolation position: x=" << extrapX
589 <<
", y=" << extrapY <<
", z=" << extrapZ);
590 foundExtrapolation =
true;
595 <<
" nMeasurements=" << nMeasurements
596 <<
" nHGTDHits=" << nHGTDHits
597 <<
" nHoles=" << nHoles
598 <<
" nOutliers=" << nOutliers
599 <<
" extrapolation found: " << (foundExtrapolation ?
"yes" :
"no"));
603 data.hasClusterVec = {hasHitInLayer[0], hasHitInLayer[1], hasHitInLayer[2], hasHitInLayer[3]};
604 data.chi2Vec = {chi2PerLayer[0], chi2PerLayer[1], chi2PerLayer[2], chi2PerLayer[3]};
605 data.timeVec = {timePerLayer[0], timePerLayer[1], timePerLayer[2], timePerLayer[3]};
606 data.rawTimeVec = {rawTimePerLayer[0], rawTimePerLayer[1], rawTimePerLayer[2], rawTimePerLayer[3]};
607 data.extrapX = extrapX;
608 data.extrapY = extrapY;
609 data.extrapZ = extrapZ;
610 data.numHGTDHits = nHGTDHits;
615 std::size_t HGTDTrackExtensionAlg::getHGTDLayerIndex(
const Acts::GeometryIdentifier& geoID)
const {
621 bool isPositiveEndcap = (volume == 25);
622 bool isNegativeEndcap = (volume == 2);
625 if (isPositiveEndcap) {
634 }
else if (isNegativeEndcap) {
648 bool HGTDTrackExtensionAlg::getExtrapolationPosition(
649 const EventContext& ctx,
650 const detail::RecoTrackContainer::TrackProxy&
track,
651 float&
x,
float&
y,
float&
z)
const {
659 bool foundHGTDSurface =
false;
662 Acts::GeometryContext geoContext = m_trackingGeometryTool->getGeometryContext(ctx).context();
663 track.container().trackStateContainer().visitBackwards(
665 [
this, &foundHGTDSurface, &
x, &
y, &
z, geoContext](
const auto& state) {
667 if (foundHGTDSurface || !state.hasReferenceSurface()) {
671 const auto& surface = state.referenceSurface();
672 Acts::GeometryIdentifier geoID = surface.geometryId();
675 if (isHGTDSurface(geoID)) {
677 if (state.hasPredicted()) {
679 const auto& predicted = state.predicted();
683 Acts::Vector2 localPos(predicted[Acts::eBoundLoc0], predicted[Acts::eBoundLoc1]);
686 Acts::Vector3 globalPos = surface.localToGlobal(
689 Acts::Vector3::Zero());
696 ATH_MSG_DEBUG(
"Predicted position on HGTD surface: (" << x <<
", " << y <<
", " << z <<
")");
697 foundHGTDSurface = true;
700 Acts::Vector3 globalPos = surface.center(geoContext);
705 ATH_MSG_DEBUG(
"Fallback to surface center for HGTD: (" << x <<
", " << y <<
", " << z <<
")");
706 foundHGTDSurface = true;
711 return foundHGTDSurface;
715 bool HGTDTrackExtensionAlg::isHGTDSurface(
const Acts::GeometryIdentifier& geoID)
const {
721 if (volume == 2 || volume == 25) {
733 std::pair<float, float> HGTDTrackExtensionAlg::correctTOF(
737 float measuredTimeErr,
738 const Acts::TrackingGeometry* ,
739 const Acts::GeometryContext& geoContext)
const {
743 if (!trackParticle || !cluster) {
745 return {measuredTime, measuredTimeErr};
749 const Acts::Surface* surface =
nullptr;
751 surface = m_surfAcc.get(cluster);
754 return {measuredTime, measuredTimeErr};
758 ATH_MSG_WARNING(
"Could not determine surface for HGTD cluster with id "
760 return {measuredTime, measuredTimeErr};
764 Acts::Vector3 globalHitPos;
769 globalHitPos = surface->localToGlobal(
771 Acts::Vector2(localPos[0], localPos[1]),
776 globalHitPos = surface->center(geoContext);
788 double d0 = trackParticle->
d0();
789 double z0 = trackParticle->
z0();
790 double phi0 = trackParticle->
phi0();
794 ATH_MSG_DEBUG(
"Track origin (perigee): (" << trackOrigin.x() <<
", "
795 << trackOrigin.y() <<
", " << trackOrigin.z() <<
")");
798 float dx = globalHitPos.x() - trackOrigin.x();
799 float dy = globalHitPos.y() - trackOrigin.y();
800 float dz = globalHitPos.z() - trackOrigin.z();
807 float correctedTime = measuredTime - tof;
810 << trackOrigin.y() <<
", " << trackOrigin.z() <<
")");
812 << globalHitPos.y() <<
", " << globalHitPos.z() <<
")");
814 <<
" ns, Corrected time = " << correctedTime);
816 return {correctedTime, measuredTimeErr};
820 if (state.hasUncalibratedSourceLink()) {
821 auto sl = state.getUncalibratedSourceLink().template get<ATLASUncalibSourceLink>();
822 assert( sl !=
nullptr);
831 ATH_MSG_DEBUG(
"Source link contains non-HGTD measurement type: " <<
static_cast<int>(clusterType));
835 if (state.hasReferenceSurface()) {
836 const auto& surface = state.referenceSurface();
837 Acts::GeometryIdentifier geoID = surface.geometryId();
840 if (isHGTDSurface(geoID)) {
841 ATH_MSG_DEBUG(
"This is an HGTD surface with ID: " << geoID.volume() <<
":" << geoID.layer());
844 EventContext ctx = Gaudi::Hive::currentContext();
856 const Acts::GeometryContext& geoContext = m_trackingGeometryTool->getGeometryContext(ctx).context();
857 Acts::Vector3 statePos = surface.center(geoContext);
861 double minDistance = 100.0;
865 const Acts::Surface* clusterSurface = m_surfAcc.get(cluster);
867 if (!clusterSurface)
continue;
870 Acts::GeometryIdentifier clusterGeoID = clusterSurface->geometryId();
871 if (clusterGeoID.volume() == geoID.volume() && clusterGeoID.layer() == geoID.layer()) {
873 Acts::Vector3 clusterPos = clusterSurface->center(geoContext);
876 double dx = clusterPos.x() - statePos.x();
877 double dy = clusterPos.y() - statePos.y();
883 closestCluster = cluster;
889 if (closestCluster) {
890 ATH_MSG_DEBUG(
"Found closest cluster at distance " << minDistance <<
" mm");
891 return closestCluster;
898 ATH_MSG_DEBUG(
"State doesn't have uncalibrated source link");