From 0ad06d5f8b625c17aa1f3e29b7187e64c96c94ac Mon Sep 17 00:00:00 2001
From: Xiaocong Ai <xiaocong.ai@cern.ch>
Date: Tue, 25 Jun 2024 00:03:24 +0200
Subject: [PATCH] update KalmanFitterTool
---
.../src/FaserActsTrack.h | 24 +
.../src/KalmanFitterTool.cxx | 568 ++++++++++--------
.../src/KalmanFitterTool.h | 75 ++-
.../src/TrackFittingFunction.cxx | 3 +-
4 files changed, 387 insertions(+), 283 deletions(-)
create mode 100644 Tracking/Acts/FaserActsKalmanFilter/src/FaserActsTrack.h
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/FaserActsTrack.h b/Tracking/Acts/FaserActsKalmanFilter/src/FaserActsTrack.h
new file mode 100644
index 000000000..59191a18c
--- /dev/null
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/FaserActsTrack.h
@@ -0,0 +1,24 @@
+// Copyright (C) 2019-2020 CERN for the benefit of the Acts project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#pragma once
+
+#include "Acts/EventData/TrackContainer.hpp"
+#include "Acts/EventData/TrackParameters.hpp"
+#include "Acts/EventData/VectorMultiTrajectory.hpp"
+#include "Acts/EventData/VectorTrackContainer.hpp"
+
+#include <vector>
+
+
+using FaserActsTrackContainer =
+ Acts::TrackContainer<Acts::VectorTrackContainer,
+ Acts::VectorMultiTrajectory, std::shared_ptr>;
+
+using FaserActsConstTrackContainer =
+ Acts::TrackContainer<Acts::ConstVectorTrackContainer,
+ Acts::ConstVectorMultiTrajectory, std::shared_ptr>;
+
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/KalmanFitterTool.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/KalmanFitterTool.cxx
index 5315404c7..c15b29566 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/KalmanFitterTool.cxx
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/KalmanFitterTool.cxx
@@ -1,6 +1,7 @@
#include "KalmanFitterTool.h"
#include "FaserActsGeometry/FASERMagneticFieldWrapper.h"
+#include "FaserActsKalmanFilter/IndexSourceLink.h"
#include "Acts/Propagator/EigenStepper.hpp"
#include "Acts/Propagator/Navigator.hpp"
@@ -21,9 +22,8 @@ StatusCode KalmanFitterTool::initialize() {
ATH_CHECK(m_trackingGeometryTool.retrieve());
ATH_CHECK(m_createTrkTrackTool.retrieve());
ATH_CHECK(detStore()->retrieve(m_idHelper,"FaserSCT_ID"));
- if (m_statesWriter && !m_noDiagnostics) ATH_CHECK(m_trajectoryStatesWriterTool.retrieve());
- if (m_summaryWriter && !m_noDiagnostics) ATH_CHECK(m_trajectorySummaryWriterTool.retrieve());
- m_fit = makeTrackFitterFunction(m_trackingGeometryTool->trackingGeometry());
+//todo if (m_statesWriter && !m_noDiagnostics) ATH_CHECK(m_trajectoryStatesWriterTool.retrieve());
+//todo if (m_summaryWriter && !m_noDiagnostics) ATH_CHECK(m_trajectorySummaryWriterTool.retrieve());
if (m_actsLogging == "VERBOSE" && !m_noDiagnostics) {
m_logger = Acts::getDefaultLogger("KalmanFitter", Acts::Logging::VERBOSE);
} else if (m_actsLogging == "DEBUG" && !m_noDiagnostics) {
@@ -31,6 +31,10 @@ StatusCode KalmanFitterTool::initialize() {
} else {
m_logger = Acts::getDefaultLogger("KalmanFitter", Acts::Logging::INFO);
}
+ auto magneticField = std::make_shared<FASERMagneticFieldWrapper>();
+ double reverseFilteringMomThreshold = 0.1; //@todo: needs to be validated
+ //@todo: the multiple scattering and energy loss are originallly turned off
+ m_fit = makeTrackFitterFunction(m_trackingGeometryTool->trackingGeometry(), magneticField, false, false, reverseFilteringMomThreshold, Acts::FreeToBoundCorrection(false), *m_logger);
return StatusCode::SUCCESS;
}
@@ -48,7 +52,6 @@ KalmanFitterTool::getUnbiasedResidual(const EventContext &ctx, const Acts::Geome
bool /*isMC*/, double origin) const {
std::vector<TSOS4Residual> resi;
resi.clear();
- std::vector<FaserActsRecMultiTrajectory> myTrajectories;
ATH_MSG_DEBUG("start kalmanfilter "<<inputTrack->measurementsOnTrack()->size()<<" "<<origin<<" "<<clusz);
if (!inputTrack->measurementsOnTrack() || inputTrack->measurementsOnTrack()->size() < m_minMeasurements) {
@@ -61,7 +64,6 @@ KalmanFitterTool::getUnbiasedResidual(const EventContext &ctx, const Acts::Geome
return resi;
}
-
auto pSurface = Acts::Surface::makeShared<Acts::PlaneSurface>(
Acts::Vector3 {0, 0, origin}, Acts::Vector3{0, 0, -1});
@@ -76,72 +78,63 @@ KalmanFitterTool::getUnbiasedResidual(const EventContext &ctx, const Acts::Geome
ATH_MSG_DEBUG("charge: " << trackParameters.charge());
ATH_MSG_DEBUG("size of measurements : " << sourceLinks.size());
+
+ auto actsTrackContainer = std::make_shared<Acts::VectorTrackContainer>();
+ auto actsTrackStateContainer = std::make_shared<Acts::VectorMultiTrajectory>();
+ FaserActsTrackContainer tracks(actsTrackContainer, actsTrackStateContainer);
+
+ //@todo: the initialSurface should be targetSurface
+ KalmanFitterTool::GeneralFitterOptions options{
+ gctx, mfContext, calibContext, &(*pSurface),
+ Acts::PropagatorPlainOptions()};
+
+ std::vector<Acts::SourceLink> actsSls;
+ for(const auto& sl: sourceLinks){
+ actsSls.push_back(Acts::SourceLink{sl});
+ }
+
+ ATH_MSG_DEBUG("Invoke fitter");
FaserActsOutlierFinder faserActsOutlierFinder{0};
faserActsOutlierFinder.cluster_z=clusz;
faserActsOutlierFinder.StateChiSquaredPerNumberDoFCut=10000.;
- Acts::KalmanFitterOptions<MeasurementCalibrator, FaserActsOutlierFinder, Acts::VoidReverseFilteringLogic>
- kfOptions(gctx, mfContext, calibContext, MeasurementCalibrator(measurements),
- faserActsOutlierFinder, Acts::VoidReverseFilteringLogic(), Acts::LoggerWrapper{*m_logger},
- //FaserActsOutlierFinder(), Acts::VoidReverseFilteringLogic(), Acts::LoggerWrapper{*m_logger},
- //Acts::VoidOutlierFinder(), Acts::VoidReverseFilteringLogic(), Acts::LoggerWrapper{*m_logger},
- Acts::PropagatorPlainOptions(), &(*pSurface));
- kfOptions.multipleScattering = false;
- kfOptions.energyLoss = false;
-
- ATH_MSG_DEBUG("Invoke fitter");
- auto result = (*m_fit)(sourceLinks, trackParameters, kfOptions);
+ auto result = (*m_fit)(actsSls, trackParameters, options, MeasurementCalibratorAdapter(MeasurementCalibrator(), measurements), faserActsOutlierFinder, tracks);
if (result.ok()) {
- const auto& fitOutput = result.value();
- if (fitOutput.fittedParameters) {
- const auto& params = fitOutput.fittedParameters.value();
- ATH_MSG_DEBUG("ReFitted parameters");
- ATH_MSG_DEBUG(" params: " << params.parameters().transpose());
- ATH_MSG_DEBUG(" position: " << params.position(gctx).transpose());
- ATH_MSG_DEBUG(" momentum: " << params.momentum().transpose());
- ATH_MSG_DEBUG(" charge: " << params.charge());
- using IndexedParams = std::unordered_map<size_t, Acts::BoundTrackParameters>;
- IndexedParams indexedParams;
- indexedParams.emplace(fitOutput.lastMeasurementIndex, std::move(params));
- std::vector<size_t> trackTips;
- trackTips.reserve(1);
- trackTips.emplace_back(fitOutput.lastMeasurementIndex);
- myTrajectories.emplace_back(std::move(fitOutput.fittedStates),
- std::move(trackTips),
- std::move(indexedParams));
- for (const FaserActsRecMultiTrajectory &traj : myTrajectories) {
- const auto& mj = traj.multiTrajectory();
- const auto& trackTips = traj.tips();
- if(trackTips.empty())continue;
- auto& trackTip = trackTips.front();
- mj.visitBackwards(trackTip, [&](const auto &state) {
- auto typeFlags = state.typeFlags();
- if (not typeFlags.test(Acts::TrackStateFlag::OutlierFlag))
- {
- return true;
- }
-
- if (state.hasUncalibrated()&&state.hasSmoothed()) {
- Acts::BoundTrackParameters parameter(
- state.referenceSurface().getSharedPtr(),
- state.smoothed(),
- state.smoothedCovariance());
- // auto covariance = state.smoothedCovariance();
- auto H = state.effectiveProjector();
- auto residual = state.effectiveCalibrated() - H * state.smoothed();
- const Tracker::FaserSCT_Cluster* cluster = state.uncalibrated().hit();
- // const auto& surface = state.referenceSurface();
- Acts::BoundVector meas = state.projector().transpose() * state.calibrated();
- Acts::Vector2 local(meas[Acts::eBoundLoc0], meas[Acts::eBoundLoc1]);
- // const Acts::Vector3 dir = Acts::makeDirectionUnitFromPhiTheta(meas[Acts::eBoundPhi], meas[Acts::eBoundTheta]);
- resi.push_back({local.x(),local.y(),parameter.position(gctx).x(), parameter.position(gctx).y(), parameter.position(gctx).z(), cluster,residual(Acts::eBoundLoc0),parameter});
- ATH_MSG_DEBUG(" residual/global position: " << residual(Acts::eBoundLoc0)<<" "<<parameter.position(gctx).x()<<" "<<parameter.position(gctx).y()<<" "<<parameter.position(gctx).z());
- }
- return true;
- });
-
+ const auto& track = result.value();
+ if (track.hasReferenceSurface()) {
+ //const auto& params = fitOutput.fittedParameters.value();
+ //ATH_MSG_DEBUG("ReFitted parameters");
+ //ATH_MSG_DEBUG(" params: " << params.parameters().transpose());
+ //ATH_MSG_DEBUG(" position: " << params.position(gctx).transpose());
+ //ATH_MSG_DEBUG(" momentum: " << params.momentum().transpose());
+ //ATH_MSG_DEBUG(" charge: " << params.charge());
+ for (const auto& state : track.trackStatesReversed()) {
+ auto flag = state.typeFlags();
+ //@todo: this only check outlier state?
+ if (not flag.test(Acts::TrackStateFlag::OutlierFlag))
+ {
+ continue;
+ }
+
+ if (state.hasUncalibratedSourceLink()&&state.hasSmoothed()) {
+ //todo: make the particle hypothesis configurable
+ Acts::BoundTrackParameters parameter(
+ state.referenceSurface().getSharedPtr(),
+ state.smoothed(),
+ state.smoothedCovariance(), Acts::ParticleHypothesis::muon());
+ // auto covariance = state.smoothedCovariance();
+ auto H = state.effectiveProjector();
+ auto residual = state.effectiveCalibrated() - H * state.smoothed();
+ IndexSourceLink sl = state.getUncalibratedSourceLink().template get<IndexSourceLink>();
+ const Tracker::FaserSCT_Cluster* cluster = sl.hit();
+ // const auto& surface = state.referenceSurface();
+ Acts::BoundVector meas = state.projector().transpose() * state.effectiveCalibrated();
+ Acts::Vector2 local(meas[Acts::eBoundLoc0], meas[Acts::eBoundLoc1]);
+ // const Acts::Vector3 dir = Acts::makeDirectionUnitFromPhiTheta(meas[Acts::eBoundPhi], meas[Acts::eBoundTheta]);
+ resi.push_back({local.x(),local.y(),parameter.position(gctx).x(), parameter.position(gctx).y(), parameter.position(gctx).z(), cluster,residual(Acts::eBoundLoc0),parameter});
+ ATH_MSG_DEBUG(" residual/global position: " << residual(Acts::eBoundLoc0)<<" "<<parameter.position(gctx).x()<<" "<<parameter.position(gctx).y()<<" "<<parameter.position(gctx).z());
+ }
}
-
} else {
ATH_MSG_DEBUG("No fitted parameters for track");
}
@@ -159,7 +152,6 @@ KalmanFitterTool::getUnbiasedResidual(const EventContext &ctx, const Acts::Geome
bool /*isMC*/, double origin, std::vector<const Tracker::FaserSCT_Cluster*>& clusters, const Acts::BoundTrackParameters ini_Param) const {
std::vector<TSOS4Residual> resi;
resi.clear();
- std::vector<FaserActsRecMultiTrajectory> myTrajectories;
if (!inputTrack->measurementsOnTrack() || inputTrack->measurementsOnTrack()->size() < m_minMeasurements) {
ATH_MSG_DEBUG("Input track has no or too little measurements and cannot be fitted");
@@ -171,7 +163,6 @@ KalmanFitterTool::getUnbiasedResidual(const EventContext &ctx, const Acts::Geome
return resi;
}
-
auto pSurface = Acts::Surface::makeShared<Acts::PlaneSurface>(
Acts::Vector3 {0, 0, origin}, Acts::Vector3{0, 0, -1});
@@ -184,69 +175,63 @@ KalmanFitterTool::getUnbiasedResidual(const EventContext &ctx, const Acts::Geome
ATH_MSG_DEBUG("momentum: " << ini_Param.momentum().transpose());
ATH_MSG_DEBUG("charge: " << ini_Param.charge());
+ auto actsTrackContainer = std::make_shared<Acts::VectorTrackContainer>();
+ auto actsTrackStateContainer = std::make_shared<Acts::VectorMultiTrajectory>();
+ FaserActsTrackContainer tracks(actsTrackContainer, actsTrackStateContainer);
+
+ //@todo: the initialSurface should be targetSurface
+ KalmanFitterTool::GeneralFitterOptions options{
+ gctx, mfContext, calibContext, &(*pSurface),
+ Acts::PropagatorPlainOptions()};
+
+ std::vector<Acts::SourceLink> actsSls;
+ for(const auto& sl: sourceLinks){
+ actsSls.push_back(Acts::SourceLink{sl});
+ }
+
+ ATH_MSG_DEBUG("Invoke fitter");
FaserActsOutlierFinder faserActsOutlierFinder{0};
faserActsOutlierFinder.cluster_z=-1000000.;
faserActsOutlierFinder.StateChiSquaredPerNumberDoFCut=10000.;
- Acts::KalmanFitterOptions<MeasurementCalibrator, FaserActsOutlierFinder, Acts::VoidReverseFilteringLogic>
- kfOptions(gctx, mfContext, calibContext, MeasurementCalibrator(measurements),
- faserActsOutlierFinder, Acts::VoidReverseFilteringLogic(), Acts::LoggerWrapper{*m_logger},
- Acts::PropagatorPlainOptions(), &(*pSurface));
- kfOptions.multipleScattering = false;
- kfOptions.energyLoss = false;
-
- ATH_MSG_DEBUG("Invoke fitter");
- auto result = (*m_fit)(sourceLinks, ini_Param, kfOptions);
+
+ auto result = (*m_fit)(actsSls, ini_Param, options, MeasurementCalibratorAdapter(MeasurementCalibrator(), measurements), faserActsOutlierFinder, tracks);
if (result.ok()) {
- const auto& fitOutput = result.value();
- if (fitOutput.fittedParameters) {
- const auto& params = fitOutput.fittedParameters.value();
- ATH_MSG_DEBUG("ReFitted parameters");
- ATH_MSG_DEBUG(" params: " << params.parameters().transpose());
- ATH_MSG_DEBUG(" position: " << params.position(gctx).transpose());
- ATH_MSG_DEBUG(" momentum: " << params.momentum().transpose());
- ATH_MSG_DEBUG(" charge: " << params.charge());
- using IndexedParams = std::unordered_map<size_t, Acts::BoundTrackParameters>;
- IndexedParams indexedParams;
- indexedParams.emplace(fitOutput.lastMeasurementIndex, std::move(params));
- std::vector<size_t> trackTips;
- trackTips.reserve(1);
- trackTips.emplace_back(fitOutput.lastMeasurementIndex);
- myTrajectories.emplace_back(std::move(fitOutput.fittedStates),
- std::move(trackTips),
- std::move(indexedParams));
- for (const FaserActsRecMultiTrajectory &traj : myTrajectories) {
- const auto& mj = traj.multiTrajectory();
- const auto& trackTips = traj.tips();
- if(trackTips.empty())continue;
- auto& trackTip = trackTips.front();
- mj.visitBackwards(trackTip, [&](const auto &state) {
- auto typeFlags = state.typeFlags();
- if (not typeFlags.test(Acts::TrackStateFlag::OutlierFlag))
- {
- return true;
- }
-
- if (state.hasUncalibrated()&&state.hasSmoothed()) {
- Acts::BoundTrackParameters parameter(
- state.referenceSurface().getSharedPtr(),
- state.smoothed(),
- state.smoothedCovariance());
- // auto covariance = state.smoothedCovariance();
- auto H = state.effectiveProjector();
- auto residual = state.effectiveCalibrated() - H * state.smoothed();
- const Tracker::FaserSCT_Cluster* cluster = state.uncalibrated().hit();
- // const auto& surface = state.referenceSurface();
- Acts::BoundVector meas = state.projector().transpose() * state.calibrated();
- Acts::Vector2 local(meas[Acts::eBoundLoc0], meas[Acts::eBoundLoc1]);
- // const Acts::Vector3 dir = Acts::makeDirectionUnitFromPhiTheta(meas[Acts::eBoundPhi], meas[Acts::eBoundTheta]);
- resi.push_back({local.x(),local.y(),parameter.position(gctx).x(), parameter.position(gctx).y(), parameter.position(gctx).z(), cluster,residual(Acts::eBoundLoc0),parameter});
- }
- return true;
- });
-
+ const auto& track = result.value();
+ if (track.hasReferenceSurface()) {
+ //const auto& params = fitOutput.fittedParameters.value();
+ //ATH_MSG_DEBUG("ReFitted parameters");
+ //ATH_MSG_DEBUG(" params: " << params.parameters().transpose());
+ //ATH_MSG_DEBUG(" position: " << params.position(gctx).transpose());
+ //ATH_MSG_DEBUG(" momentum: " << params.momentum().transpose());
+ //ATH_MSG_DEBUG(" charge: " << params.charge());
+ for (const auto& state : track.trackStatesReversed()) {
+ auto flag = state.typeFlags();
+ //@todo: this only check outlier state?
+ if (not flag.test(Acts::TrackStateFlag::OutlierFlag))
+ {
+ continue;
+ }
+
+ if (state.hasUncalibratedSourceLink()&&state.hasSmoothed()) {
+ //todo: make the particle hypothesis configurable
+ Acts::BoundTrackParameters parameter(
+ state.referenceSurface().getSharedPtr(),
+ state.smoothed(),
+ state.smoothedCovariance(), Acts::ParticleHypothesis::muon());
+ // auto covariance = state.smoothedCovariance();
+ auto H = state.effectiveProjector();
+ auto residual = state.effectiveCalibrated() - H * state.smoothed();
+ IndexSourceLink sl = state.getUncalibratedSourceLink().template get<IndexSourceLink>();
+ const Tracker::FaserSCT_Cluster* cluster = sl.hit();
+ // const auto& surface = state.referenceSurface();
+ Acts::BoundVector meas = state.projector().transpose() * state.effectiveCalibrated();
+ Acts::Vector2 local(meas[Acts::eBoundLoc0], meas[Acts::eBoundLoc1]);
+ // const Acts::Vector3 dir = Acts::makeDirectionUnitFromPhiTheta(meas[Acts::eBoundPhi], meas[Acts::eBoundTheta]);
+ resi.push_back({local.x(),local.y(),parameter.position(gctx).x(), parameter.position(gctx).y(), parameter.position(gctx).z(), cluster,residual(Acts::eBoundLoc0),parameter});
+ ATH_MSG_DEBUG(" residual/global position: " << residual(Acts::eBoundLoc0)<<" "<<parameter.position(gctx).x()<<" "<<parameter.position(gctx).y()<<" "<<parameter.position(gctx).z());
+ }
}
-
} else {
ATH_MSG_DEBUG("No fitted parameters for track");
}
@@ -264,7 +249,6 @@ KalmanFitterTool::getUnbiasedResidual(const EventContext &ctx, const Acts::Geome
bool /*isMC*/, double origin) const {
std::vector<TSOS4Residual> resi;
resi.clear();
- std::vector<FaserActsRecMultiTrajectory> myTrajectories;
if (!inputTrack->measurementsOnTrack() || inputTrack->measurementsOnTrack()->size() < m_minMeasurements) {
ATH_MSG_DEBUG("Input track has no or too little measurements and cannot be fitted");
@@ -290,70 +274,63 @@ KalmanFitterTool::getUnbiasedResidual(const EventContext &ctx, const Acts::Geome
ATH_MSG_DEBUG("momentum: " << trackParameters.momentum().transpose());
ATH_MSG_DEBUG("charge: " << trackParameters.charge());
+ auto actsTrackContainer = std::make_shared<Acts::VectorTrackContainer>();
+ auto actsTrackStateContainer = std::make_shared<Acts::VectorMultiTrajectory>();
+ FaserActsTrackContainer tracks(actsTrackContainer, actsTrackStateContainer);
+
+ //@todo: the initialSurface should be targetSurface
+ KalmanFitterTool::GeneralFitterOptions options{
+ gctx, mfContext, calibContext, &(*pSurface),
+ Acts::PropagatorPlainOptions()};
+
+ std::vector<Acts::SourceLink> actsSls;
+ for(const auto& sl: sourceLinks){
+ actsSls.push_back(Acts::SourceLink{sl});
+ }
+ ATH_MSG_DEBUG("Invoke fitter");
FaserActsOutlierFinder faserActsOutlierFinder{0};
faserActsOutlierFinder.cluster_z=-1000000.;
faserActsOutlierFinder.StateChiSquaredPerNumberDoFCut=10000.;
- Acts::KalmanFitterOptions<MeasurementCalibrator, FaserActsOutlierFinder, Acts::VoidReverseFilteringLogic>
- kfOptions(gctx, mfContext, calibContext, MeasurementCalibrator(measurements),
- faserActsOutlierFinder, Acts::VoidReverseFilteringLogic(), Acts::LoggerWrapper{*m_logger},
- Acts::PropagatorPlainOptions(), &(*pSurface));
- kfOptions.multipleScattering = false;
- kfOptions.energyLoss = false;
-
- ATH_MSG_DEBUG("Invoke fitter");
- auto result = (*m_fit)(sourceLinks, trackParameters, kfOptions);
+
+ auto result = (*m_fit)(actsSls, trackParameters, options, MeasurementCalibratorAdapter(MeasurementCalibrator(), measurements), faserActsOutlierFinder, tracks);
if (result.ok()) {
- const auto& fitOutput = result.value();
- if (fitOutput.fittedParameters) {
- const auto& params = fitOutput.fittedParameters.value();
- ATH_MSG_DEBUG("ReFitted parameters");
- ATH_MSG_DEBUG(" params: " << params.parameters().transpose());
- ATH_MSG_DEBUG(" position: " << params.position(gctx).transpose());
- ATH_MSG_DEBUG(" momentum: " << params.momentum().transpose());
- ATH_MSG_DEBUG(" charge: " << params.charge());
- using IndexedParams = std::unordered_map<size_t, Acts::BoundTrackParameters>;
- IndexedParams indexedParams;
- indexedParams.emplace(fitOutput.lastMeasurementIndex, std::move(params));
- std::vector<size_t> trackTips;
- trackTips.reserve(1);
- trackTips.emplace_back(fitOutput.lastMeasurementIndex);
- myTrajectories.emplace_back(std::move(fitOutput.fittedStates),
- std::move(trackTips),
- std::move(indexedParams));
- for (const FaserActsRecMultiTrajectory &traj : myTrajectories) {
- const auto& mj = traj.multiTrajectory();
- const auto& trackTips = traj.tips();
- if(trackTips.empty())continue;
- auto& trackTip = trackTips.front();
- mj.visitBackwards(trackTip, [&](const auto &state) {
- auto typeFlags = state.typeFlags();
- if (not typeFlags.test(Acts::TrackStateFlag::OutlierFlag))
- {
- return true;
- }
-
- if (state.hasUncalibrated()&&state.hasSmoothed()) {
- Acts::BoundTrackParameters parameter(
- state.referenceSurface().getSharedPtr(),
- state.smoothed(),
- state.smoothedCovariance());
- // auto covariance = state.smoothedCovariance();
- auto H = state.effectiveProjector();
- auto residual = state.effectiveCalibrated() - H * state.smoothed();
- const Tracker::FaserSCT_Cluster* cluster = state.uncalibrated().hit();
- // const auto& surface = state.referenceSurface();
- Acts::BoundVector meas = state.projector().transpose() * state.calibrated();
- Acts::Vector2 local(meas[Acts::eBoundLoc0], meas[Acts::eBoundLoc1]);
- // const Acts::Vector3 dir = Acts::makeDirectionUnitFromPhiTheta(meas[Acts::eBoundPhi], meas[Acts::eBoundTheta]);
-resi.push_back({local.x(),local.y(),parameter.position(gctx).x(), parameter.position(gctx).y(), parameter.position(gctx).z(), cluster,residual(Acts::eBoundLoc0),parameter});
- }
- return true;
- });
-
+ const auto& track = result.value();
+ if (track.hasReferenceSurface()) {
+ //const auto& params = track.fittedParameters.value();
+ //ATH_MSG_DEBUG("ReFitted parameters");
+ //ATH_MSG_DEBUG(" params: " << params.parameters().transpose());
+ //ATH_MSG_DEBUG(" position: " << params.position(gctx).transpose());
+ //ATH_MSG_DEBUG(" momentum: " << params.momentum().transpose());
+ //ATH_MSG_DEBUG(" charge: " << params.charge());
+ for (const auto& state : track.trackStatesReversed()) {
+ auto flag = state.typeFlags();
+ //@todo: this only check outlier state?
+ if (not flag.test(Acts::TrackStateFlag::OutlierFlag))
+ {
+ continue;
+ }
+
+ if (state.hasUncalibratedSourceLink()&&state.hasSmoothed()) {
+ //todo: make the particle hypothesis configurable
+ Acts::BoundTrackParameters parameter(
+ state.referenceSurface().getSharedPtr(),
+ state.smoothed(),
+ state.smoothedCovariance(), Acts::ParticleHypothesis::muon());
+ // auto covariance = state.smoothedCovariance();
+ auto H = state.effectiveProjector();
+ auto residual = state.effectiveCalibrated() - H * state.smoothed();
+ IndexSourceLink sl = state.getUncalibratedSourceLink().template get<IndexSourceLink>();
+ const Tracker::FaserSCT_Cluster* cluster = sl.hit();
+ // const auto& surface = state.referenceSurface();
+ Acts::BoundVector meas = state.projector().transpose() * state.effectiveCalibrated();
+ Acts::Vector2 local(meas[Acts::eBoundLoc0], meas[Acts::eBoundLoc1]);
+ // const Acts::Vector3 dir = Acts::makeDirectionUnitFromPhiTheta(meas[Acts::eBoundPhi], meas[Acts::eBoundTheta]);
+ resi.push_back({local.x(),local.y(),parameter.position(gctx).x(), parameter.position(gctx).y(), parameter.position(gctx).z(), cluster,residual(Acts::eBoundLoc0),parameter});
+ ATH_MSG_DEBUG(" residual/global position: " << residual(Acts::eBoundLoc0)<<" "<<parameter.position(gctx).x()<<" "<<parameter.position(gctx).y()<<" "<<parameter.position(gctx).z());
+ }
}
-
} else {
ATH_MSG_DEBUG("No fitted parameters for track");
}
@@ -367,7 +344,6 @@ KalmanFitterTool::fit(const EventContext &ctx, const Acts::GeometryContext &gctx
Trk::Track* inputTrack, const Acts::BoundVector& inputVector,
bool isMC) const {
std::unique_ptr<Trk::Track> newTrack = nullptr;
- std::vector<FaserActsRecMultiTrajectory> myTrajectories;
if (!inputTrack->measurementsOnTrack() || inputTrack->measurementsOnTrack()->size() < m_minMeasurements) {
ATH_MSG_DEBUG("Input track has no or too little measurements and cannot be fitted");
@@ -396,49 +372,58 @@ KalmanFitterTool::fit(const EventContext &ctx, const Acts::GeometryContext &gctx
ATH_MSG_DEBUG("momentum: " << trackParameters.momentum().transpose());
ATH_MSG_DEBUG("charge: " << trackParameters.charge());
+ auto actsTrackContainer = std::make_shared<Acts::VectorTrackContainer>();
+ auto actsTrackStateContainer = std::make_shared<Acts::VectorMultiTrajectory>();
+ FaserActsTrackContainer tracks(actsTrackContainer, actsTrackStateContainer);
+
+ //@todo: the initialSurface should be targetSurface
+ KalmanFitterTool::GeneralFitterOptions options{
+ gctx, mfContext, calibContext, &(*pSurface),
+ Acts::PropagatorPlainOptions()};
+
+ std::vector<Acts::SourceLink> actsSls;
+ for(const auto& sl: sourceLinks){
+ actsSls.push_back(Acts::SourceLink{sl});
+ }
+
+ ATH_MSG_DEBUG("Invoke fitter");
FaserActsOutlierFinder faserActsOutlierFinder{0};
faserActsOutlierFinder.cluster_z=-1000000.;
faserActsOutlierFinder.StateChiSquaredPerNumberDoFCut=10000.;
- Acts::KalmanFitterOptions<MeasurementCalibrator, FaserActsOutlierFinder, Acts::VoidReverseFilteringLogic>
- kfOptions(gctx, mfContext, calibContext, MeasurementCalibrator(measurements),
- faserActsOutlierFinder, Acts::VoidReverseFilteringLogic(), Acts::LoggerWrapper{*m_logger},
- Acts::PropagatorPlainOptions(), &(*pSurface));
- kfOptions.multipleScattering = false;
- kfOptions.energyLoss = false;
-
- ATH_MSG_DEBUG("Invoke fitter");
- auto result = (*m_fit)(sourceLinks, trackParameters, kfOptions);
+
+ auto result = (*m_fit)(actsSls, trackParameters, options, MeasurementCalibratorAdapter(MeasurementCalibrator(), measurements), faserActsOutlierFinder, tracks);
if (result.ok()) {
- const auto& fitOutput = result.value();
- if (fitOutput.fittedParameters) {
- const auto& params = fitOutput.fittedParameters.value();
- ATH_MSG_DEBUG("ReFitted parameters");
- ATH_MSG_DEBUG(" params: " << params.parameters().transpose());
- ATH_MSG_DEBUG(" position: " << params.position(gctx).transpose());
- ATH_MSG_DEBUG(" momentum: " << params.momentum().transpose());
- ATH_MSG_DEBUG(" charge: " << params.charge());
- using IndexedParams = std::unordered_map<size_t, Acts::BoundTrackParameters>;
- IndexedParams indexedParams;
- indexedParams.emplace(fitOutput.lastMeasurementIndex, std::move(params));
- std::vector<size_t> trackTips;
- trackTips.reserve(1);
- trackTips.emplace_back(fitOutput.lastMeasurementIndex);
- myTrajectories.emplace_back(std::move(fitOutput.fittedStates),
- std::move(trackTips),
- std::move(indexedParams));
- } else {
- ATH_MSG_DEBUG("No fitted parameters for track");
- }
- newTrack = m_createTrkTrackTool->createTrack(gctx, myTrajectories.back());
+ const auto& track = result.value();
+ //if (track.hasReferenceSurface()) {
+ //const auto& params = track.fittedParameters.value();
+ //ATH_MSG_DEBUG("ReFitted parameters");
+ //ATH_MSG_DEBUG(" params: " << params.parameters().transpose());
+ //ATH_MSG_DEBUG(" position: " << params.position(gctx).transpose());
+ //ATH_MSG_DEBUG(" momentum: " << params.momentum().transpose());
+ //ATH_MSG_DEBUG(" charge: " << params.charge());
+ //using IndexedParams = std::unordered_map<size_t, Acts::BoundTrackParameters>;
+ //IndexedParams indexedParams;
+ //indexedParams.emplace(fitOutput.lastMeasurementIndex, std::move(params));
+ //std::vector<size_t> trackTips;
+ //trackTips.reserve(1);
+ //trackTips.emplace_back(fitOutput.lastMeasurementIndex);
+ //myTrajectories.emplace_back(std::move(fitOutput.fittedStates),
+ // std::move(trackTips),
+ // std::move(indexedParams));
+ //} else {
+ // ATH_MSG_DEBUG("No fitted parameters for track");
+ //}
+ newTrack = m_createTrkTrackTool->createTrack(gctx, track);
}
- if (m_statesWriter && !m_noDiagnostics) {
- StatusCode statusStatesWriterTool = m_trajectoryStatesWriterTool->write(gctx, myTrajectories, isMC);
- }
- if (m_summaryWriter && !m_noDiagnostics) {
- StatusCode statusSummaryWriterTool = m_trajectorySummaryWriterTool->write(gctx, myTrajectories, isMC);
- }
+//todo
+// if (m_statesWriter && !m_noDiagnostics) {
+// StatusCode statusStatesWriterTool = m_trajectoryStatesWriterTool->write(gctx, myTrajectories, isMC);
+// }
+// if (m_summaryWriter && !m_noDiagnostics) {
+// StatusCode statusSummaryWriterTool = m_trajectorySummaryWriterTool->write(gctx, myTrajectories, isMC);
+// }
return newTrack;
}
@@ -446,25 +431,84 @@ KalmanFitterTool::fit(const EventContext &ctx, const Acts::GeometryContext &gctx
namespace {
-using Updater = Acts::GainMatrixUpdater;
-using Smoother = Acts::GainMatrixSmoother;
using Stepper = Acts::EigenStepper<>;
using Propagator = Acts::Propagator<Stepper, Acts::Navigator>;
-using Fitter = Acts::KalmanFitter<Propagator, Updater, Smoother>;
+using Fitter = Acts::KalmanFitter<Propagator, Acts::VectorMultiTrajectory>;
+
+
+struct SimpleReverseFilteringLogic {
+ double momentumThreshold = 0;
+
+ bool doBackwardFiltering(
+ Acts::VectorMultiTrajectory::ConstTrackStateProxy trackState) const {
+ auto momentum = fabs(1 / trackState.filtered()[Acts::eBoundQOverP]);
+ return (momentum <= momentumThreshold);
+ }
+};
+
struct TrackFitterFunctionImpl
: public KalmanFitterTool::TrackFitterFunction {
Fitter trackFitter;
- TrackFitterFunctionImpl(Fitter &&f) : trackFitter(std::move(f)) {}
+ Acts::GainMatrixUpdater kfUpdater;
+ Acts::GainMatrixSmoother kfSmoother;
+ SimpleReverseFilteringLogic reverseFilteringLogic;
+
+ bool multipleScattering = false;
+ bool energyLoss = false;
+ Acts::FreeToBoundCorrection freeToBoundCorrection;
+
+ IndexSourceLink::SurfaceAccessor slSurfaceAccessor;
+
+ TrackFitterFunctionImpl(Fitter &&f, const Acts::TrackingGeometry& trkGeo) : trackFitter(std::move(f)),slSurfaceAccessor{trkGeo} {}
+
+ template <typename calibrator_t>
+ auto makeKfOptions(const KalmanFitterTool::GeneralFitterOptions& options,
+ const calibrator_t& calibrator, const FaserActsOutlierFinder& outlierFinder) const {
+ Acts::KalmanFitterExtensions<Acts::VectorMultiTrajectory> extensions;
+ extensions.updater.connect<
+ &Acts::GainMatrixUpdater::operator()<Acts::VectorMultiTrajectory>>(
+ &kfUpdater);
+ extensions.smoother.connect<
+ &Acts::GainMatrixSmoother::operator()<Acts::VectorMultiTrajectory>>(
+ &kfSmoother);
+ extensions.reverseFilteringLogic
+ .connect<&SimpleReverseFilteringLogic::doBackwardFiltering>(
+ &reverseFilteringLogic);
+ //@todo: this needs fix!!!
+ //extensions.outlierFinder.connect<&FaserActsOutlierFinder::operator()<Acts::ConstVectorMultiTrajectory>>(
+ // &outlierFinder);
+
+ Acts::KalmanFitterOptions<Acts::VectorMultiTrajectory> kfOptions(
+ options.geoContext, options.magFieldContext, options.calibrationContext,
+ extensions, options.propOptions, &(*options.referenceSurface));
+
+ kfOptions.referenceSurfaceStrategy =
+ Acts::KalmanFitterTargetSurfaceStrategy::first;
+ kfOptions.multipleScattering = multipleScattering;
+ kfOptions.energyLoss = energyLoss;
+ kfOptions.freeToBoundCorrection = freeToBoundCorrection;
+ kfOptions.extensions.calibrator.connect<&calibrator_t::calibrate>(
+ &calibrator);
+ kfOptions.extensions.surfaceAccessor
+ .connect<&IndexSourceLink::SurfaceAccessor::operator()>(
+ &slSurfaceAccessor);
+
+ return kfOptions;
+ }
KalmanFitterTool::TrackFitterResult operator()(
- const std::vector<IndexSourceLink> &sourceLinks,
+ const std::vector<Acts::SourceLink> &sourceLinks,
const KalmanFitterTool::TrackParameters &initialParameters,
- const KalmanFitterTool::TrackFitterOptions &options)
- const override {
- return trackFitter.fit(sourceLinks, initialParameters, options);
- };
+ const KalmanFitterTool::GeneralFitterOptions& options,
+ const MeasurementCalibratorAdapter& calibrator,
+ const FaserActsOutlierFinder& outlierFinder,
+ FaserActsTrackContainer& tracks) const override {
+ const auto kfOptions = makeKfOptions(options, calibrator, outlierFinder);
+ return trackFitter.fit(sourceLinks.begin(), sourceLinks.end(), initialParameters,
+ kfOptions, tracks);
+ }
};
} // namespace
@@ -472,20 +516,38 @@ struct TrackFitterFunctionImpl
std::shared_ptr<KalmanFitterTool::TrackFitterFunction>
KalmanFitterTool::makeTrackFitterFunction(
- std::shared_ptr<const Acts::TrackingGeometry> trackingGeometry) {
- auto magneticField = std::make_shared<FASERMagneticFieldWrapper>();
- auto stepper = Stepper(std::move(magneticField));
- Acts::Navigator::Config cfg{trackingGeometry};
+ std::shared_ptr<const Acts::TrackingGeometry> trackingGeometry,
+ std::shared_ptr<const Acts::MagneticFieldProvider> magneticField,
+ bool multipleScattering, bool energyLoss,
+ double reverseFilteringMomThreshold,
+ Acts::FreeToBoundCorrection freeToBoundCorrection,
+ const Acts::Logger& logger) {
+ // Stepper should be copied into the fitters
+ const Stepper stepper(std::move(magneticField));
+
+ // Standard fitter
+ const auto& geo = *trackingGeometry;
+ Acts::Navigator::Config cfg{std::move(trackingGeometry)};
cfg.resolvePassive = false;
cfg.resolveMaterial = true;
cfg.resolveSensitive = true;
- Acts::Navigator navigator(cfg);
- Propagator propagator(std::move(stepper), std::move(navigator));
- Fitter trackFitter(std::move(propagator));
- return std::make_shared<TrackFitterFunctionImpl>(std::move(trackFitter));
+ Acts::Navigator navigator(cfg, logger.cloneWithSuffix("Navigator"));
+ Propagator propagator(stepper, std::move(navigator),
+ logger.cloneWithSuffix("Propagator"));
+ Fitter trackFitter(std::move(propagator), logger.cloneWithSuffix("Fitter"));
+
+ // build the fitter function. owns the fitter object.
+ auto fitterFunction = std::make_shared<TrackFitterFunctionImpl>(
+ std::move(trackFitter), geo);
+ fitterFunction->multipleScattering = multipleScattering;
+ fitterFunction->energyLoss = energyLoss;
+ fitterFunction->reverseFilteringLogic.momentumThreshold =
+ reverseFilteringMomThreshold;
+ fitterFunction->freeToBoundCorrection = freeToBoundCorrection;
+
+ return fitterFunction;
}
-
Acts::MagneticFieldContext KalmanFitterTool::getMagneticFieldContext(const EventContext& ctx) const {
SG::ReadCondHandle<FaserFieldCacheCondObj> readHandle{m_fieldCondObjInputKey, ctx};
if (!readHandle.isValid()) {
@@ -501,9 +563,9 @@ Acts::MagneticFieldContext KalmanFitterTool::getMagneticFieldContext(const Event
std::tuple<std::vector<IndexSourceLink>, std::vector<Measurement>>
KalmanFitterTool::getMeasurementsFromTrack(Trk::Track *track, std::vector<const Tracker::FaserSCT_Cluster*>& clusters) const {
const int kSize = 1;
- ATH_MSG_DEBUG("clusters in refit "<<clusters.size());
+ ATH_MSG_DEBUG("clusters in refit "<<clusters.size());
std::array<Acts::BoundIndices, kSize> Indices = {Acts::eBoundLoc0};
- using ThisMeasurement = Acts::Measurement<IndexSourceLink, Acts::BoundIndices, kSize>;
+ using ThisMeasurement = Acts::Measurement<Acts::BoundIndices, kSize>;
using IdentifierMap = std::map<Identifier, Acts::GeometryIdentifier>;
std::shared_ptr<IdentifierMap> identifierMap = m_trackingGeometryTool->getIdentifierMap();
@@ -521,7 +583,7 @@ KalmanFitterTool::getMeasurementsFromTrack(Trk::Track *track, std::vector<const
IndexSourceLink sourceLink(geoId, measurements.size(), cluster);
Eigen::Matrix<double, 1, 1> pos {cluster->localPosition().x()};
Eigen::Matrix<double, 1, 1> cov {0.08 * 0.08 / 12};
- ThisMeasurement actsMeas(sourceLink, Indices, pos, cov);
+ ThisMeasurement actsMeas(Acts::SourceLink{std::move(sourceLink)}, Indices, pos, cov);
sourceLinks.push_back(sourceLink);
measurements.emplace_back(std::move(actsMeas));
}
@@ -537,7 +599,7 @@ KalmanFitterTool::getMeasurementsFromTrack(Trk::Track *track, std::vector<const
IndexSourceLink sourceLink(geoId, measurements.size(), cluster);
Eigen::Matrix<double, 1, 1> pos {meas->localParameters()[Trk::locX],};
Eigen::Matrix<double, 1, 1> cov {0.08 * 0.08 / 12};
- ThisMeasurement actsMeas(sourceLink, Indices, pos, cov);
+ ThisMeasurement actsMeas(Acts::SourceLink{std::move(sourceLink)}, Indices, pos, cov);
sourceLinks.push_back(sourceLink);
measurements.emplace_back(std::move(actsMeas));
}
@@ -550,7 +612,7 @@ std::tuple<std::vector<IndexSourceLink>, std::vector<Measurement>>
KalmanFitterTool::getMeasurementsFromTrack(Trk::Track *track, Identifier& /*wafer_id*/) const {
const int kSize = 1;
std::array<Acts::BoundIndices, kSize> Indices = {Acts::eBoundLoc0};
- using ThisMeasurement = Acts::Measurement<IndexSourceLink, Acts::BoundIndices, kSize>;
+ using ThisMeasurement = Acts::Measurement<Acts::BoundIndices, kSize>;
using IdentifierMap = std::map<Identifier, Acts::GeometryIdentifier>;
std::shared_ptr<IdentifierMap> identifierMap = m_trackingGeometryTool->getIdentifierMap();
@@ -568,7 +630,7 @@ KalmanFitterTool::getMeasurementsFromTrack(Trk::Track *track, Identifier& /*wafe
IndexSourceLink sourceLink(geoId, measurements.size(), cluster);
Eigen::Matrix<double, 1, 1> pos {meas->localParameters()[Trk::locX],};
Eigen::Matrix<double, 1, 1> cov {0.08 * 0.08 / 12};
- ThisMeasurement actsMeas(sourceLink, Indices, pos, cov);
+ ThisMeasurement actsMeas(Acts::SourceLink{std::move(sourceLink)}, Indices, pos, cov);
sourceLinks.push_back(sourceLink);
measurements.emplace_back(std::move(actsMeas));
}
@@ -580,7 +642,7 @@ std::tuple<std::vector<IndexSourceLink>, std::vector<Measurement>>
KalmanFitterTool::getMeasurementsFromTrack(Trk::Track *track) const {
const int kSize = 1;
std::array<Acts::BoundIndices, kSize> Indices = {Acts::eBoundLoc0};
- using ThisMeasurement = Acts::Measurement<IndexSourceLink, Acts::BoundIndices, kSize>;
+ using ThisMeasurement = Acts::Measurement<Acts::BoundIndices, kSize>;
using IdentifierMap = std::map<Identifier, Acts::GeometryIdentifier>;
std::shared_ptr<IdentifierMap> identifierMap = m_trackingGeometryTool->getIdentifierMap();
@@ -597,7 +659,7 @@ KalmanFitterTool::getMeasurementsFromTrack(Trk::Track *track) const {
IndexSourceLink sourceLink(geoId, measurements.size(), cluster);
Eigen::Matrix<double, 1, 1> pos {meas->localParameters()[Trk::locX],};
Eigen::Matrix<double, 1, 1> cov {0.08 * 0.08 / 12};
- ThisMeasurement actsMeas(sourceLink, Indices, pos, cov);
+ ThisMeasurement actsMeas(Acts::SourceLink{std::move(sourceLink)}, Indices, pos, cov);
sourceLinks.push_back(sourceLink);
measurements.emplace_back(std::move(actsMeas));
}
@@ -622,16 +684,16 @@ KalmanFitterTool::getParametersFromTrack(const Trk::TrackParameters *inputParame
// 0.};
Acts::BoundVector params;
if (inputVector == Acts::BoundVector::Zero()) {
- params(0.0) = center[Trk::locY];
- params(1.0) = center[Trk::locX];
- params(2.0) = atlasParam[Trk::phi0];
- params(3.0) = atlasParam[Trk::theta];
- params(4.0) = inputParameters->charge() / (inputParameters->momentum().norm() * 1_MeV);
- params(5.0) = 0.;
+ params(0) = center[Trk::locY];
+ params(1) = center[Trk::locX];
+ params(2) = atlasParam[Trk::phi0];
+ params(3) = atlasParam[Trk::theta];
+ params(4) = inputParameters->charge() / (inputParameters->momentum().norm() * 1_MeV);
+ params(5) = 0.;
} else {
params = inputVector;
}
- Acts::BoundSymMatrix cov = Acts::BoundSymMatrix::Identity();
+ Acts::BoundSquareMatrix cov = Acts::BoundSquareMatrix::Identity();
cov.topLeftCorner(5, 5) = *inputParameters->covariance();
for(int i=0; i < cov.rows(); i++){
@@ -644,5 +706,5 @@ KalmanFitterTool::getParametersFromTrack(const Trk::TrackParameters *inputParame
cov(i,i) = m_seedCovarianceScale * cov(i,i);
}
- return Acts::BoundTrackParameters(pSurface, params, inputParameters->charge(), cov);
+ return Acts::BoundTrackParameters(pSurface, params, cov, Acts::ParticleHypothesis::muon());
}
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/KalmanFitterTool.h b/Tracking/Acts/FaserActsKalmanFilter/src/KalmanFitterTool.h
index 4901d324f..b58a31269 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/KalmanFitterTool.h
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/KalmanFitterTool.h
@@ -5,6 +5,7 @@
#include "AthenaBaseComps/AthAlgTool.h"
#include "Acts/EventData/TrackParameters.hpp"
#include "Acts/TrackFitting/KalmanFitter.hpp"
+#include "Acts/EventData/MeasurementHelpers.hpp"
#include "FaserActsGeometryInterfaces/IFaserActsTrackingGeometryTool.h"
#include "FaserActsKalmanFilter/IndexSourceLink.h"
#include "FaserActsKalmanFilter/Measurement.h"
@@ -17,6 +18,8 @@
#include "CreateTrkTrackTool.h"
#include "TrackerIdentifier/FaserSCT_ID.h"
#include "Identifier/Identifier.h"
+#include "FaserActsTrack.h"
+#include "FaserActsGeometry/FASERMagneticFieldWrapper.h"
struct TSOS4Residual{
double fit_local_x;
@@ -31,39 +34,32 @@ struct TSOS4Residual{
};
class FaserSCT_ID;
+
+
//set the cluster to be removed as outlier in order to get the unbiased residual
/// Outlier finder using a Chi2 cut.
struct FaserActsOutlierFinder {
double StateChiSquaredPerNumberDoFCut = 10000.;
double cluster_z = -10000.;
- template <typename track_state_t>
- bool operator()(const track_state_t& state) const {
+
+ template <typename traj_t>
+ bool operator()(typename traj_t::ConstTrackStateProxy state) const {
//remove the whole IFT
if(cluster_z<-10000){
- if(state.uncalibrated().hit()->globalPosition().z()<-100)return true;
+ IndexSourceLink sl = state.getUncalibratedSourceLink().template get<IndexSourceLink>();
+ if(sl.hit()->globalPosition().z()<-100)return true;
+ if(abs(sl.hit()->globalPosition().z()-cluster_z)<3)return true;
}
if (not state.hasCalibrated() or not state.hasPredicted()) {
return false;
}
- return Acts::visit_measurement(
- state.calibrated(), state.calibratedCovariance(),
- state.calibratedSize(),
- [&](const auto calibrated, const auto calibratedCovariance) {
- //remove the cluster
- if(fabs(state.uncalibrated().hit()->globalPosition().z()-cluster_z)<3)return true;
- constexpr size_t kMeasurementSize = decltype(calibrated)::RowsAtCompileTime;
- const auto H =
- state.projector()
- .template topLeftCorner<kMeasurementSize, Acts::BoundIndices::eBoundSize>()
- .eval();
- const auto residual = calibrated - H * state.predicted();
- double chi2 = (residual.transpose() * ((calibratedCovariance + H * state.predictedCovariance() * H.transpose())).inverse() * residual).value();
- return bool(chi2 > StateChiSquaredPerNumberDoFCut * kMeasurementSize);
- });
+
+ return bool(state.chi2() > StateChiSquaredPerNumberDoFCut * state.calibratedSize());
}
};
+
class KalmanFitterTool : virtual public AthAlgTool {
public:
KalmanFitterTool(const std::string &type, const std::string &name, const IInterface *parent);
@@ -71,34 +67,55 @@ public:
virtual StatusCode initialize() override;
virtual StatusCode finalize() override;
+ struct GeneralFitterOptions {
+ std::reference_wrapper<const Acts::GeometryContext> geoContext;
+ std::reference_wrapper<const Acts::MagneticFieldContext> magFieldContext;
+ std::reference_wrapper<const Acts::CalibrationContext> calibrationContext;
+ const Acts::Surface* referenceSurface = nullptr;
+ Acts::PropagatorPlainOptions propOptions;
+ };
+ using TrackFitterResult = Acts::Result<FaserActsTrackContainer::TrackProxy>;
using TrackParameters = Acts::BoundTrackParameters;
- using IndexedParams = std::unordered_map<size_t, TrackParameters>;
- using TrackFitterOptions =
- Acts::KalmanFitterOptions<MeasurementCalibrator, FaserActsOutlierFinder, Acts::VoidReverseFilteringLogic>;
- //Acts::KalmanFitterOptions<MeasurementCalibrator, Acts::VoidOutlierFinder, Acts::VoidReverseFilteringLogic>;
- using TrackFitterResult = Acts::Result<Acts::KalmanFitterResult<IndexSourceLink>>;
+ //using IndexedParams = std::unordered_map<size_t, TrackParameters>;
+
class TrackFitterFunction {
public:
virtual ~TrackFitterFunction() = default;
- virtual TrackFitterResult operator()(const std::vector<IndexSourceLink>&,
- const TrackParameters&,
- const TrackFitterOptions&) const = 0;
+ virtual TrackFitterResult operator()(
+ const std::vector<Acts::SourceLink> &sourceLinks,
+ const TrackParameters &initialParameters,
+ const GeneralFitterOptions & kfOptions,
+ const MeasurementCalibratorAdapter& calibrator,
+ const FaserActsOutlierFinder& outlierFinder,
+ FaserActsTrackContainer& tracks
+ ) const = 0;
};
+
static std::shared_ptr<TrackFitterFunction> makeTrackFitterFunction(
- std::shared_ptr<const Acts::TrackingGeometry> trackingGeometry);
+ std::shared_ptr<const Acts::TrackingGeometry> trackingGeometry,
+ std::shared_ptr<const Acts::MagneticFieldProvider> magneticField,
+ bool multipleScattering, bool energyLoss,
+ double reverseFilteringMomThreshold,
+ Acts::FreeToBoundCorrection freeToBoundCorrection,
+ const Acts::Logger& logger);
+
virtual Acts::MagneticFieldContext getMagneticFieldContext(const EventContext& ctx) const;
+
std::unique_ptr<Trk::Track> fit(const EventContext &ctx, const Acts::GeometryContext &gctx,
Trk::Track *inputTrack,
const Acts::BoundVector& inputVector = Acts::BoundVector::Zero(),
bool isMC=false) const;
+
std::vector<TSOS4Residual> getUnbiasedResidual(const EventContext &ctx, const Acts::GeometryContext &gctx,
Trk::Track *inputTrack,
const Acts::BoundVector& inputVector = Acts::BoundVector::Zero(),
bool isMC=false, double origin=0) const;
+
std::vector<TSOS4Residual> getUnbiasedResidual(const EventContext &ctx, const Acts::GeometryContext &gctx,
Trk::Track *inputTrack, double clusz,
const Acts::BoundVector& inputVector ,
bool isMC, double origin) const;
+
std::vector<TSOS4Residual> getUnbiasedResidual(const EventContext &ctx, const Acts::GeometryContext &gctx,
Trk::Track *inputTrack,
const Acts::BoundVector& inputVector ,
@@ -125,8 +142,8 @@ private:
SG::ReadCondHandleKey<FaserFieldCacheCondObj> m_fieldCondObjInputKey {this, "FaserFieldCacheCondObj", "fieldCondObj", "Name of the Magnetic Field conditions object key"};
ToolHandle<IFaserActsTrackingGeometryTool> m_trackingGeometryTool {this, "TrackingGeometryTool", "FaserActsTrackingGeometryTool"};
- ToolHandle<RootTrajectoryStatesWriterTool> m_trajectoryStatesWriterTool {this, "RootTrajectoryStatesWriterTool", "RootTrajectoryStatesWriterTool"};
- ToolHandle<RootTrajectorySummaryWriterTool> m_trajectorySummaryWriterTool {this, "RootTrajectorySummaryWriterTool", "RootTrajectorySummaryWriterTool"};
+//todo ToolHandle<RootTrajectoryStatesWriterTool> m_trajectoryStatesWriterTool {this, "RootTrajectoryStatesWriterTool", "RootTrajectoryStatesWriterTool"};
+//todo ToolHandle<RootTrajectorySummaryWriterTool> m_trajectorySummaryWriterTool {this, "RootTrajectorySummaryWriterTool", "RootTrajectorySummaryWriterTool"};
ToolHandle<CreateTrkTrackTool> m_createTrkTrackTool {this, "CreateTrkTrackTool", "CreateTrkTrackTool"};
// Acts::KalmanFitterExtensions<Acts::VectorMultiTrajectory> getExtensions();
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/TrackFittingFunction.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/TrackFittingFunction.cxx
index c8d2bbe95..7a1e31dfc 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/TrackFittingFunction.cxx
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/TrackFittingFunction.cxx
@@ -1,4 +1,5 @@
#include "FaserActsKalmanFilterAlg.h"
+#include "FaserActsTrack.h"
//#include "FaserActsGeometry/FASERMagneticFieldWrapper.h"
//#include "FaserActsKalmanFilter/IndexSourceLink.h"
@@ -92,7 +93,7 @@ struct TrackFitterFunctionImpl
const FaserActsKalmanFilterAlg::TrackParameters &initialParameters,
const FaserActsKalmanFilterAlg::GeneralFitterOptions& options,
const MeasurementCalibratorAdapter& calibrator,
- FaserActsKalmanFilterAlg::TrackContainer& tracks) const override {
+ FaserActsTrackContainer& tracks) const override {
const auto kfOptions = makeKfOptions(options, calibrator);
return trackFitter.fit(sourceLinks.begin(), sourceLinks.end(), initialParameters,
kfOptions, tracks);
--
GitLab