Skip to content

Commit 1707f30

Browse files
andiwandandiwand
committed
fix
1 parent 3e4ba0a commit 1707f30

1 file changed

Lines changed: 106 additions & 87 deletions

File tree

‎Core/include/Acts/TrackFinding/CombinatorialKalmanFilter.hpp‎

Lines changed: 106 additions & 87 deletions
Original file line numberDiff line numberDiff line change
@@ -26,6 +26,7 @@
2626
#include "Acts/TrackFinding/CombinatorialKalmanFilterExtensions.hpp"
2727
#include "Acts/TrackFitting/BetheHeitlerApprox.hpp"
2828
#include "Acts/TrackFitting/GsfOptions.hpp"
29+
#include "Acts/TrackFitting/detail/GsfComponentMerging.hpp"
2930
#include "Acts/TrackFitting/detail/GsfUtils.hpp"
3031
#include "Acts/Utilities/CalibrationContext.hpp"
3132
#include "Acts/Utilities/Logger.hpp"
@@ -262,8 +263,8 @@ class CombinatorialKalmanFilter {
262263
VolumeConstraintAborter volumeConstraintAborter;
263264

264265
TemporaryStates* temporaryStates{nullptr};
265-
std::vector<BetheHeitlerApprox::Component> betheHeitlerCache;
266-
std::vector<GsfComponent> componentCache;
266+
std::vector<BetheHeitlerApprox::Component>* betheHeitlerCache{nullptr};
267+
std::vector<GsfComponent>* componentCache{nullptr};
267268

268269
/// Actor logger instance
269270
const Logger* actorLogger{nullptr};
@@ -285,7 +286,7 @@ class CombinatorialKalmanFilter {
285286
typename navigator_t>
286287
Result<void> act(propagator_state_t& state, const stepper_t& stepper,
287288
const navigator_t& navigator, result_type& result,
288-
const Logger& /*logger*/) {
289+
const Logger& /*logger*/) const {
289290
ACTS_VERBOSE("CKF Actor called");
290291

291292
assert(result.trackStates && "No MultiTrajectory set");
@@ -442,10 +443,19 @@ class CombinatorialKalmanFilter {
442443
<< currentBranch.tipIndex());
443444

444445
// Reset the stepping state
445-
stepper.initialize(state.stepping, currentState.filtered(),
446-
currentState.filteredCovariance(),
447-
stepper.particleHypothesis(state.stepping),
448-
currentState.referenceSurface());
446+
if constexpr (!IsMultiStepper) {
447+
stepper.initialize(state.stepping, currentState.filtered(),
448+
currentState.filteredCovariance(),
449+
stepper.particleHypothesis(state.stepping),
450+
currentState.referenceSurface());
451+
} else {
452+
const MultiComponentBoundTrackParameters multiBoundParameters(
453+
currentState.referenceSurface().getSharedPtr(),
454+
currentState.filtered(), currentState.filteredCovariance(),
455+
stepper.particleHypothesis(state.stepping));
456+
457+
stepper.initialize(state.stepping, multiBoundParameters);
458+
}
449459

450460
// Reset the navigation state
451461
// Set targetSurface to nullptr for forward filtering
@@ -494,7 +504,7 @@ class CombinatorialKalmanFilter {
494504
typename navigator_t>
495505
Result<void> filter(const Surface& surface, propagator_state_t& state,
496506
const stepper_t& stepper, const navigator_t& navigator,
497-
result_type& result) {
507+
result_type& result) const {
498508
using PM = TrackStatePropMask;
499509

500510
const bool isSensitive = surface.isSensitive();
@@ -849,26 +859,26 @@ class CombinatorialKalmanFilter {
849859
void performMaterialInteraction(propagator_state_t& state,
850860
const stepper_t& stepper,
851861
const Surface& surface,
852-
MaterialUpdateMode updateMode) {
862+
MaterialUpdateMode updateMode) const {
853863
if constexpr (!IsMultiStepper) {
854864
detail::performMaterialInteraction(
855865
state, stepper, surface, updateMode, NoiseUpdateMode::addNoise,
856866
multipleScattering, energyLoss, logger());
857867
} else {
858868
if (ACTS_CHECK_BIT(updateMode, MaterialUpdateMode::PostUpdate)) {
859869
temporaryStates->clear();
860-
betheHeitlerCache.clear();
861-
componentCache.clear();
870+
betheHeitlerCache->clear();
871+
componentCache->clear();
862872
std::size_t nInvalidBetheHeitler = 0;
863873
double maxPathXOverX0 = 0;
864874
double sumPathXOverX0 = 0;
865875

866876
detail::Gsf::convoluteComponents(
867-
surface, state, temporaryStates, *betheHeitlerApprox,
868-
betheHeitlerCache, weightCutoff, componentCache,
877+
state, stepper, *temporaryStates, *betheHeitlerApprox,
878+
*betheHeitlerCache, weightCutoff, *componentCache,
869879
nInvalidBetheHeitler, maxPathXOverX0, sumPathXOverX0, logger());
870880

871-
if (componentCache.empty()) {
881+
if (componentCache->empty()) {
872882
ACTS_WARNING(
873883
"No components left after applying energy loss. "
874884
"Is the weight cutoff "
@@ -880,11 +890,11 @@ class CombinatorialKalmanFilter {
880890
// reduce component number
881891
const auto finalCmpNumber = std::min(
882892
static_cast<std::size_t>(stepper.maxComponents), maxComponents);
883-
extensions.mixtureReducer(componentCache, finalCmpNumber, surface);
893+
extensions.mixtureReducer(*componentCache, finalCmpNumber, surface);
884894

885-
detail::Gsf::removeLowWeightComponents(componentCache, weightCutoff);
895+
detail::Gsf::removeLowWeightComponents(*componentCache, weightCutoff);
886896

887-
detail::Gsf::updateStepper(state, stepper, surface, componentCache,
897+
detail::Gsf::updateStepper(state, stepper, surface, *componentCache,
888898
logger());
889899
}
890900

@@ -899,84 +909,84 @@ class CombinatorialKalmanFilter {
899909
TrackStateProxy& trackState) const {
900910
if constexpr (!IsMultiStepper) {
901911
return extensions.updater(state.geoContext, trackState, *updaterLogger);
902-
}
912+
} else {
913+
using PrtProjector = detail::Gsf::MultiTrajectoryProjector<
914+
detail::Gsf::StatesType::ePredicted, TrackStateContainerBackend>;
915+
using FltProjector = detail::Gsf::MultiTrajectoryProjector<
916+
detail::Gsf::StatesType::eFiltered, TrackStateContainerBackend>;
903917

904-
using PrtProjector = detail::Gsf::MultiTrajectoryProjector<
905-
detail::Gsf::StatesType::ePredicted, TrackStateContainerBackend>;
906-
using FltProjector = detail::Gsf::MultiTrajectoryProjector<
907-
detail::Gsf::StatesType::eFiltered, TrackStateContainerBackend>;
918+
const auto& surface = trackState.referenceSurface();
908919

909-
const auto& surface = trackState.referenceSurface();
920+
temporaryStates->clear();
910921

911-
temporaryStates->clear();
922+
for (auto cmp : stepper.componentIterable(state.stepping)) {
923+
auto& singleState = cmp.singleState(state).stepping;
924+
const auto& singleStepper = cmp.singleStepper(stepper);
912925

913-
for (auto cmp : stepper.componentIterable(state.stepping)) {
914-
auto& singleState = cmp.singleState(state).stepping;
915-
const auto& singleStepper = cmp.singleStepper(stepper);
926+
TrackStatePropMask mask =
927+
TrackStatePropMask::Predicted | TrackStatePropMask::Filtered |
928+
TrackStatePropMask::Jacobian | TrackStatePropMask::Calibrated;
929+
TrackStateProxy trackStateProxy =
930+
temporaryStates->traj.makeTrackState(mask, kTrackIndexInvalid);
916931

917-
TrackStatePropMask mask =
918-
TrackStatePropMask::Predicted | TrackStatePropMask::Filtered |
919-
TrackStatePropMask::Jacobian | TrackStatePropMask::Calibrated;
920-
TrackStateProxy trackStateProxy =
921-
temporaryStates->traj.makeTrackState(mask, kTrackIndexInvalid);
932+
// TODO call calibrator again?
922933

923-
// TODO call calibrator again?
934+
trackStateProxy.setReferenceSurface(surface.getSharedPtr());
935+
// Bind the transported state to the current surface
936+
auto res = singleStepper.boundState(singleState, surface);
937+
if (!res.ok()) {
938+
ACTS_ERROR("Propagate to surface " << surface.geometryId()
939+
<< " failed: " << res.error());
940+
return res.error();
941+
}
942+
const auto& [boundParams, jacobian, pathLength] = *res;
924943

925-
trackStateProxy.setReferenceSurface(surface.getSharedPtr());
926-
// Bind the transported state to the current surface
927-
auto res = singleStepper.boundState(singleState, surface);
928-
if (!res.ok()) {
929-
ACTS_ERROR("Propagate to surface " << surface.geometryId()
930-
<< " failed: " << res.error());
931-
return res.error();
932-
}
933-
const auto& [boundParams, jacobian, pathLength] = *res;
934-
935-
// Fill the track state
936-
trackStateProxy.predicted() = boundParams.parameters();
937-
trackStateProxy.predictedCovariance() = *boundParams.covariance();
938-
trackStateProxy.allocateCalibrated(trackState.calibratedSize());
939-
trackStateProxy.setProjectorSubspaceIndices(
940-
trackState.projectorSubspaceIndices());
941-
trackStateProxy.effectiveCalibrated() =
942-
trackState.effectiveCalibrated();
943-
trackStateProxy.effectiveCalibratedCovariance() =
944-
trackState.effectiveCalibratedCovariance();
945-
946-
const auto updateRes = extensions.updater(
947-
state.geoContext, trackStateProxy, *updaterLogger);
948-
if (!updateRes.ok()) {
949-
return updateRes.error();
950-
}
944+
// Fill the track state
945+
trackStateProxy.predicted() = boundParams.parameters();
946+
trackStateProxy.predictedCovariance() = *boundParams.covariance();
947+
trackStateProxy.allocateCalibrated(trackState.calibratedSize());
948+
trackStateProxy.setProjectorSubspaceIndices(
949+
trackState.projectorSubspaceIndices());
950+
trackStateProxy.effectiveCalibrated() =
951+
trackState.effectiveCalibrated();
952+
trackStateProxy.effectiveCalibratedCovariance() =
953+
trackState.effectiveCalibratedCovariance();
954+
955+
const auto updateRes = extensions.updater(
956+
state.geoContext, trackStateProxy, *updaterLogger);
957+
if (!updateRes.ok()) {
958+
return updateRes.error();
959+
}
951960

952-
temporaryStates->tips.push_back(trackStateProxy.index());
953-
temporaryStates->weights[temporaryStates->tips.back()] = cmp.weight();
954-
}
961+
temporaryStates->tips.push_back(trackStateProxy.index());
962+
temporaryStates->weights[temporaryStates->tips.back()] = cmp.weight();
963+
}
955964

956-
detail::Gsf::computePosteriorWeights(temporaryStates->traj,
957-
temporaryStates->tips,
958-
temporaryStates->weights);
959-
960-
detail::Gsf::normalizeWeights(temporaryStates->tips,
961-
[&](auto idx) -> double& {
962-
return temporaryStates->weights.at(idx);
963-
});
964-
965-
const auto [prtMean, prtCov] = mergeGaussianMixture(
966-
temporaryStates->tips,
967-
PrtProjector{temporaryStates->traj, temporaryStates->weights},
968-
surface, mergeMethod);
969-
trackState.predicted() = prtMean;
970-
trackState.predictedCovariance() = prtCov;
971-
972-
const auto [fltMean, fltCov] = mergeGaussianMixture(
973-
temporaryStates->tips,
974-
FltProjector{temporaryStates->traj, temporaryStates->weights},
975-
surface, mergeMethod);
976-
trackState.filtered() = fltMean;
977-
trackState.filteredCovariance() = fltCov;
965+
detail::Gsf::computePosteriorWeights(temporaryStates->traj,
966+
temporaryStates->tips,
967+
temporaryStates->weights);
968+
969+
detail::Gsf::normalizeWeights(temporaryStates->tips,
970+
[&](auto idx) -> double& {
971+
return temporaryStates->weights.at(idx);
972+
});
973+
974+
const auto [prtMean, prtCov] = detail::Gsf::mergeGaussianMixture(
975+
temporaryStates->tips,
976+
PrtProjector{temporaryStates->traj, temporaryStates->weights},
977+
surface, mergeMethod);
978+
trackState.predicted() = prtMean;
979+
trackState.predictedCovariance() = prtCov;
980+
981+
const auto [fltMean, fltCov] = detail::Gsf::mergeGaussianMixture(
982+
temporaryStates->tips,
983+
FltProjector{temporaryStates->traj, temporaryStates->weights},
984+
surface, mergeMethod);
985+
trackState.filtered() = fltMean;
986+
trackState.filteredCovariance() = fltCov;
978987

979-
return Result<void>::success();
988+
return Result<void>::success();
989+
}
980990
}
981991

982992
template <typename propagator_state_t, typename stepper_t>
@@ -990,8 +1000,13 @@ class CombinatorialKalmanFilter {
9901000
currentState.filteredCovariance(),
9911001
currentState.referenceSurface());
9921002
} else {
993-
detail::Gsf::updateStepper(state, stepper,
994-
currentState.referenceSurface(), logger());
1003+
// Conincidentially the temporaryStates correspond to the currentState
1004+
// but it seems fragile to rely on that. It would be better to store the
1005+
// multi bound parameters on the track state and recover them here. This
1006+
// would also behave better in case of branching.
1007+
// TODO revisit this
1008+
detail::Gsf::updateStepper(state, stepper, *temporaryStates,
1009+
weightCutoff);
9951010
}
9961011
}
9971012
};
@@ -1062,7 +1077,11 @@ class CombinatorialKalmanFilter {
10621077
combKalmanActor.betheHeitlerApprox = tfOptions.betheHeitlerApprox.get();
10631078

10641079
TemporaryStates temporaryStates;
1080+
std::vector<BetheHeitlerApprox::Component> betheHeitlerCache;
1081+
std::vector<GsfComponent> componentCache;
10651082
combKalmanActor.temporaryStates = &temporaryStates;
1083+
combKalmanActor.betheHeitlerCache = &betheHeitlerCache;
1084+
combKalmanActor.componentCache = &componentCache;
10661085

10671086
auto propState =
10681087
m_propagator

0 commit comments

Comments
 (0)