diff --git a/Core/include/Acts/EventData/TrackState.hpp b/Core/include/Acts/EventData/TrackState.hpp
index d4d861e2eac9cb4abfdf38568e48f385d6c02e06..21fe35e4575c043844611907a0f463f6c54281e3 100644
--- a/Core/include/Acts/EventData/TrackState.hpp
+++ b/Core/include/Acts/EventData/TrackState.hpp
@@ -146,6 +146,9 @@ public:
double pathLength = 0.;
} parameter;
+ /// The increment of the chi^2 after the Kalman update (filtering).
+ double filteredChi2Increment = 0;
+
/// @brief Nested measurement part
/// This is the uncalibrated and calibrated measurement
/// (in case the latter is different)
diff --git a/Core/include/Acts/Fitter/GainMatrixUpdator.hpp b/Core/include/Acts/Fitter/GainMatrixUpdator.hpp
index bacfc24b4b7fc38fc078c2578a5da62c6c437142..a82742c241e047cd405a44fba9befa89f78d05cf 100644
--- a/Core/include/Acts/Fitter/GainMatrixUpdator.hpp
+++ b/Core/include/Acts/Fitter/GainMatrixUpdator.hpp
@@ -108,20 +108,26 @@ public:
filtered_covariance
= (CovMatrix_t::Identity() - K * H) * predicted_covariance;
+ // Create new filtered parameters and covariance
+ parameters_t filtered(
+ std::make_unique<const CovMatrix_t>(filtered_covariance),
+ filtered_parameters,
+ predicted.referenceSurface().getSharedPtr());
+
+ const auto R_filtered = (calibrated.covariance()
+ - H * filtered_covariance * H.transpose());
+ const auto r_filtered = calibrated.residual(filtered);
+ const double chi2Increment
+ = (r_filtered.transpose() * R_filtered * r_filtered).value();
+
// plug calibrated measurement back into track state
trackState.measurement.calibrated = std::move(calibrated);
+ trackState.parameter.filtered = std::move(filtered);
+ trackState.filteredChi2Increment = chi2Increment;
},
*trackState.measurement.uncalibrated);
- // Create new filtered parameters and covariance
- parameters_t filtered(
- std::make_unique<const CovMatrix_t>(std::move(filtered_covariance)),
- filtered_parameters,
- predicted.referenceSurface().getSharedPtr());
-
- trackState.parameter.filtered = std::move(filtered);
-
// always succeed, no outlier logic yet
return true;
}
diff --git a/Core/include/Acts/Fitter/SequentialKalmanFitter.hpp b/Core/include/Acts/Fitter/SequentialKalmanFitter.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..700da210f27dd4405f64fd296b261b4cb262df5b
--- /dev/null
+++ b/Core/include/Acts/Fitter/SequentialKalmanFitter.hpp
@@ -0,0 +1,197 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2016-2018 Acts project team
+//
+// 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/Fitter/detail/VoidKalmanComponents.hpp"
+#include "Acts/Propagator/AbortList.hpp"
+#include "Acts/Propagator/ActionList.hpp"
+#include "Acts/Propagator/Propagator.hpp"
+
+#include <boost/variant.hpp>
+#include <memory>
+
+namespace Acts {
+
+/**
+ * Sequential Kalman Fitter algorithm to traverse through a geometry with
+ * measurements and select a single best combination (= track). For this, a
+ * prepared measurement selector is asked to supply measurements for each
+ * visited surface in the geometry, which are then ranked by chi^2 increment.
+ * The one with the smalles chi^2 increment is taken, the track parameters
+ * updated and the algorithm is continued.
+ *
+ * @tparam propagator_t The type of the propagator to use.
+ * @tparam track_state_t The type of the track state to use.
+ * @tparam measurement_selector_t The type of the measurement selector
+ * algorithm. Should implement an operator returning a vector of track
+ * states to used for a given surface.
+ * @tparam updator_t The type of the updator (used for updating the track
+ * parameters with a measurement)
+ * @tparam calibrator_t The type of the calibrator (used for calibrating the
+ * measurements depending on a track state)
+ */
+template <typename propagator_t,
+ typename track_state_t,
+ typename measurement_selector_t,
+ typename updator_t = VoidKalmanUpdator,
+ typename calibrator_t = VoidKalmanComponents>
+class SequentialKalmanFitter
+{
+public:
+ /// Default constructor is deleted
+ SequentialKalmanFitter() = delete;
+
+ /// Constructor taking a propagator as input
+ explicit SequentialKalmanFitter(propagator_t pPropagator)
+ : m_propagator(std::move(pPropagator))
+ {
+ }
+
+ /**
+ * Main function of the fitter, which takes a seeding start parameters as
+ * input and the prepared measurement selector and outputs a result object
+ * including the found track.
+ * @tparam parameters_t The type of the seed parameters
+ * @param sParameters The seed parameters to start with
+ * @param pMeasurementSelector The prefilled measurement selector, which will
+ * return the measurements for a given surface
+ * @return A result structure with the selected track states.
+ */
+ template <typename parameters_t>
+ auto
+ fit(const parameters_t& sParameters,
+ std::shared_ptr<measurement_selector_t> pMeasurementSelector) const
+ {
+ // Create the ActionList and AbortList
+ using KalmanResult = typename KalmanActor::result_type;
+ using Actors = ActionList<KalmanActor>;
+ using Aborters = AbortList<>;
+
+ // Create relevant options for the propagation options
+ PropagatorOptions<Actors, Aborters> kalmanOptions;
+ auto& kalmanActor = kalmanOptions.actionList.template get<KalmanActor>();
+ kalmanActor.m_measurementSelectorPtr = std::move(pMeasurementSelector);
+
+ // Run the fitter
+ const auto& result
+ = m_propagator.template propagate(sParameters, kalmanOptions);
+
+ // Get the result of the fit
+ auto kalmanResult = result.template get<KalmanResult>();
+
+ // Return the converted Track
+ return kalmanResult;
+ }
+
+private:
+ /// Used propagator instance
+ propagator_t m_propagator;
+
+ /**
+ * Actor used in the SequentialKalmanFitter propagation.
+ * On every call, it will look for compatible measurements on the surface
+ * and select the one with the lowest chi^2 increment.
+ */
+ class KalmanActor
+ {
+ public:
+ /**
+ * Creates a new actor (is called by the propagator automatically).
+ * @param pUpdator The updator instance used during the step.
+ * @param pCalibrator The calibrator instance used during the step.
+ */
+ KalmanActor(updator_t pUpdator = updator_t(),
+ calibrator_t pCalibrator = calibrator_t())
+ : m_updator(std::move(pUpdator)), m_calibrator(std::move(pCalibrator))
+ {
+ }
+
+ /// The type of the result of this actor including the final track and the
+ /// number of holes.
+ struct result_type
+ {
+ /// The resulting track of the algorithm
+ std::vector<track_state_t> resultStateList = {};
+ /// How many surfaces had no measurement on them
+ unsigned int numberOfHoles = 0;
+ };
+
+ /**
+ * Main function of this actor:
+ * * Is only executed if propagated to a surface
+ * * let the measurements selector give us a list of possible next
+ * measurements
+ * * call the update process on them
+ * * rank them by chi^2 and pick the best one.
+ * @tparam propagator_state_t The propagator state type
+ * @param state The state coming from the propagator
+ * @param result The result we use for storing e.g. the resulting track
+ */
+ template <typename propagator_state_t>
+ void
+ operator()(propagator_state_t& state, result_type& result) const
+ {
+ // Do not go on if there is no surface
+ const Surface* surface = state.navigation.currentSurface;
+ if (not surface) {
+ return;
+ }
+
+ // Request a list of possible next states by the measurement selector
+ assert(m_measurementSelectorPtr != nullptr);
+ std::vector<track_state_t> nextStates
+ = (*m_measurementSelectorPtr)(*surface, state);
+ if (nextStates.empty()) {
+ // Only count it as hole if a sensitive surface
+ if (surface->associatedDetectorElement() != nullptr) {
+ result.numberOfHoles++;
+ }
+ return;
+ }
+
+ // Loop through all possible measurement
+ // perform an update and rank by chi^2.
+ auto boundStates = state.stepping.boundState(*surface, true);
+ track_state_t* chosenNextState = nullptr;
+ for (auto& nextState : nextStates) {
+ nextState.parameter.predicted = std::get<0>(boundStates);
+ nextState.parameter.jacobian = std::get<1>(boundStates);
+ nextState.parameter.pathLength = std::get<2>(boundStates);
+
+ m_updator(nextState);
+
+ if (chosenNextState == nullptr) {
+ chosenNextState = &nextState;
+ } else {
+ const double chosenChi2 = chosenNextState->filteredChi2Increment;
+ const double nextChi2 = nextState.filteredChi2Increment;
+ if (chosenChi2 > nextChi2) {
+ chosenNextState = &nextState;
+ }
+ }
+ }
+
+ assert(chosenNextState != nullptr);
+
+ // Update the fitted state with the chosen next state and also store it in
+ // the track list
+ state.stepping.update(*(chosenNextState->parameter.filtered));
+ result.resultStateList.push_back(std::move(*chosenNextState));
+ }
+
+ /// The pointer to the measurement selector (needs to be filled by the user)
+ std::shared_ptr<measurement_selector_t> m_measurementSelectorPtr;
+ /// The updator instance
+ updator_t m_updator;
+ /// The calibrator instance
+ calibrator_t m_calibrator;
+ };
+};
+
+} // namespace Acts
diff --git a/Tests/Core/CommonHelpers/include/Acts/Tests/CommonHelpers/CubicTrackingGeometry.hpp b/Tests/Core/CommonHelpers/include/Acts/Tests/CommonHelpers/CubicTrackingGeometry.hpp
index 594fe8d5ddb6f22e3ce0c0f094afc08af29efe47..195f8df5838d2293c303dcc5797a4a98817c2a3b 100644
--- a/Tests/Core/CommonHelpers/include/Acts/Tests/CommonHelpers/CubicTrackingGeometry.hpp
+++ b/Tests/Core/CommonHelpers/include/Acts/Tests/CommonHelpers/CubicTrackingGeometry.hpp
@@ -5,7 +5,7 @@
// 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 <vector>
#include "Acts/Detector/TrackingGeometry.hpp"
#include "Acts/Detector/TrackingVolume.hpp"
diff --git a/Tests/Core/CommonHelpers/include/Acts/Tests/CommonHelpers/MeasurementCreatorHelper.hpp b/Tests/Core/CommonHelpers/include/Acts/Tests/CommonHelpers/MeasurementCreatorHelper.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..f3962125e385abdce379a48adba560589325dd8e
--- /dev/null
+++ b/Tests/Core/CommonHelpers/include/Acts/Tests/CommonHelpers/MeasurementCreatorHelper.hpp
@@ -0,0 +1,264 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2016-2018 Acts project team
+//
+// 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 <algorithm>
+#include <math.h>
+#include <random>
+#include <vector>
+#include "Acts/Detector/TrackingGeometry.hpp"
+#include "Acts/EventData/Measurement.hpp"
+#include "Acts/EventData/TrackParameters.hpp"
+#include "Acts/EventData/TrackState.hpp"
+#include "Acts/Extrapolator/Navigator.hpp"
+#include "Acts/Fitter/KalmanFitter.hpp"
+#include "Acts/Propagator/EigenStepper.hpp"
+#include "Acts/Propagator/Propagator.hpp"
+#include "Acts/Propagator/StraightLineStepper.hpp"
+#include "Acts/Propagator/detail/DebugOutputActor.hpp"
+#include "Acts/Propagator/detail/StandardAborters.hpp"
+#include "Acts/Surfaces/Surface.hpp"
+#include "Acts/Utilities/BinningType.hpp"
+#include "Acts/Utilities/Definitions.hpp"
+#include "Acts/Utilities/GeometryID.hpp"
+
+namespace Acts {
+namespace Test {
+
+ // A few initialisations and definitionas
+ using Identifier = unsigned int;
+
+ using TrackState = TrackState<Identifier, BoundParameters>;
+ using Resolution = std::pair<ParID_t, double>;
+ using ElementResolution = std::vector<Resolution>;
+ using VolumeResolution = std::map<geo_id_value, ElementResolution>;
+ using DetectorResolution = std::map<geo_id_value, VolumeResolution>;
+
+ using DebugOutput = detail::DebugOutputActor;
+
+ std::normal_distribution<double> gauss(0., 1.);
+ std::uniform_int_distribution<int> numberOfMeasurements(1, 5);
+ std::mt19937 generator(42);
+
+ /// @brief This struct creates FittableMeasurements on the
+ /// detector surfaces, according to the given smearing xxparameters
+ ///
+ struct MeasurementCreator
+ {
+ /// @brief Constructor
+ MeasurementCreator() = default;
+
+ /// The detector resolution
+ DetectorResolution detectorResolution;
+
+ /// Randomize number of hits
+ bool randomizeNumberOfMeasurements = false;
+
+ struct result_type
+ {
+ std::vector<TrackState> trackStates;
+ std::map<const Surface*, std::vector<FittableMeasurement<Identifier>>>
+ measurementMap;
+ };
+
+ /// @brief Operater that is callable by an ActionList. The function collects
+ /// the surfaces
+ ///
+ /// @tparam propagator_state_t Type of the propagator state
+ /// @param [in] state State of the propagator
+ /// @param [out] result Vector of matching surfaces
+ template <typename propagator_state_t>
+ void
+ operator()(propagator_state_t& state, result_type& result) const
+ {
+ // monitor the current surface
+ const Surface* surface = state.navigation.currentSurface;
+ if (surface == nullptr or not surface->associatedDetectorElement()) {
+ return;
+ }
+
+ auto geoID = surface->geoID();
+ geo_id_value volumeID = geoID.value(GeometryID::volume_mask);
+ geo_id_value layerID = geoID.value(GeometryID::layer_mask);
+ // find volume and layer information for this
+ auto vResolution = detectorResolution.find(volumeID);
+ if (vResolution == detectorResolution.end()) {
+ return;
+ }
+ // find layer resolutions
+ auto lResolution = vResolution->second.find(layerID);
+ if (lResolution == vResolution->second.end()) {
+ return;
+ }
+
+ unsigned int nMeasurements = 1;
+
+ if (randomizeNumberOfMeasurements) {
+ nMeasurements
+ = static_cast<unsigned int>(numberOfMeasurements(generator));
+ }
+
+ for (unsigned int identifier = 0; identifier < nMeasurements;
+ identifier++) {
+ addSingleMeasurement(identifier, state, result, surface, lResolution);
+ }
+ }
+
+ template <typename propagator_state_t, typename resolution_t>
+ void
+ addSingleMeasurement(unsigned int identifier,
+ propagator_state_t& state,
+ result_type& result,
+ const Surface* surface,
+ resolution_t lResolution) const
+ {
+ // Apply global to local
+ Acts::Vector2D lPos;
+ surface->globalToLocal(
+ state.stepping.position(), state.stepping.direction(), lPos);
+ if (lResolution->second.size() == 1) {
+ double sp = lResolution->second[0].second;
+ ActsSymMatrixD<1> cov1D;
+ cov1D << sp * sp;
+ double dp = sp * gauss(generator);
+
+ if (identifier != 0) {
+ dp *= 20;
+ }
+
+ if (lResolution->second[0].first == eLOC_0) {
+ // push back & move a LOC_0 measurement
+ Measurement<Identifier, eLOC_0> m(
+ surface->getSharedPtr(), identifier, cov1D, lPos[eLOC_0] + dp);
+ result.measurementMap[surface].push_back(m);
+ result.trackStates.push_back(TrackState(std::move(m)));
+ } else {
+ // push back & move a LOC_1 measurement
+ Measurement<Identifier, eLOC_1> m(
+ surface->getSharedPtr(), identifier, cov1D, lPos[eLOC_1] + dp);
+ result.measurementMap[surface].push_back(m);
+ result.trackStates.push_back(TrackState(std::move(m)));
+ }
+ } else if (lResolution->second.size() == 2) {
+ // Create the measurement and move it
+ double sx = lResolution->second[eLOC_0].second;
+ double sy = lResolution->second[eLOC_1].second;
+ ActsSymMatrixD<2> cov2D;
+ cov2D << sx * sx, 0., 0., sy * sy;
+ double dx = sx * gauss(generator);
+ double dy = sy * gauss(generator);
+
+ if (identifier != 0) {
+ dx *= 20;
+ dy *= 20;
+ }
+
+ // push back & move a LOC_0, LOC_1 measurement
+ Measurement<Identifier, eLOC_0, eLOC_1> m(surface->getSharedPtr(),
+ identifier,
+ cov2D,
+ lPos[eLOC_0] + dx,
+ lPos[eLOC_1] + dy);
+ result.measurementMap[surface].push_back(m);
+ result.trackStates.push_back(TrackState(std::move(m)));
+ }
+ }
+ };
+
+ template <class AStepper>
+ MeasurementCreator::result_type
+ createMeasurements(std::shared_ptr<const TrackingGeometry> detector,
+ AStepper stepper,
+ const Vector3D& startPosition,
+ const Vector3D& startMomentum,
+ DetectorResolution detRes,
+ bool addNoiseMeasurements = false,
+ bool debugMode = false)
+ {
+ // Build navigator for the measurement creation
+ Navigator navigator(detector);
+ navigator.resolvePassive = false;
+ navigator.resolveMaterial = true;
+ navigator.resolveSensitive = true;
+
+ // Define the measurement propagator
+ using MeasurementPropagator = Propagator<AStepper, Navigator>;
+
+ // Build propagator for the measurement creation
+ MeasurementPropagator propagator(stepper, navigator);
+
+ SingleCurvilinearTrackParameters<NeutralPolicy> startParameters(
+ nullptr, startPosition, startMomentum);
+
+ // Create action list for the measurement creation
+ using MeasurementActions = ActionList<MeasurementCreator, DebugOutput>;
+ using MeasurementAborters = AbortList<detail::EndOfWorldReached>;
+
+ // Set options for propagator
+ PropagatorOptions<MeasurementActions, MeasurementAborters> propOps;
+ propOps.debug = debugMode;
+ auto& mCreator = propOps.actionList.get<MeasurementCreator>();
+ mCreator.detectorResolution = detRes;
+ mCreator.randomizeNumberOfMeasurements = addNoiseMeasurements;
+
+ // Launch and collect - the measurements
+ auto result = propagator.propagate(startParameters, propOps);
+ if (debugMode) {
+ const auto debugString
+ = result.template get<DebugOutput::result_type>().debugString;
+ std::cout << ">>>> Measurement creation: " << std::endl;
+ std::cout << debugString;
+ }
+
+ auto measurements = result.template get<MeasurementCreator::result_type>();
+ return measurements;
+ }
+
+ template <class AStepper>
+ MeasurementCreator::result_type
+ createMeasurementsOnCubicDetector(
+ std::shared_ptr<const TrackingGeometry> detector,
+ AStepper stepper,
+ const Vector3D& startPosition,
+ const Vector3D& startMomentum,
+ bool addNoiseMeasurements = false,
+ bool debugMode = false)
+ {
+ auto pixelResX = Resolution(eLOC_0, 25. * units::_um);
+ auto pixelResY = Resolution(eLOC_1, 50. * units::_um);
+ auto stripResX = Resolution(eLOC_0, 100. * units::_um);
+ auto stripResY = Resolution(eLOC_1, 150. * units::_um);
+
+ ElementResolution pixelElementRes = {pixelResX, pixelResY};
+ ElementResolution stripElementResI = {stripResX};
+ ElementResolution stripElementResO = {stripResY};
+
+ VolumeResolution pixelVolumeRes;
+ pixelVolumeRes[2] = pixelElementRes;
+ pixelVolumeRes[4] = pixelElementRes;
+
+ VolumeResolution stripVolumeRes;
+ stripVolumeRes[2] = stripElementResI;
+ stripVolumeRes[4] = stripElementResO;
+ stripVolumeRes[6] = stripElementResI;
+ stripVolumeRes[8] = stripElementResO;
+
+ DetectorResolution detRes;
+ detRes[2] = pixelVolumeRes;
+ detRes[3] = stripVolumeRes;
+
+ return createMeasurements(detector,
+ stepper,
+ startPosition,
+ startMomentum,
+ detRes,
+ addNoiseMeasurements,
+ debugMode);
+ }
+
+} // namespace Test
+} // namespace Acts
diff --git a/Tests/Core/Fitter/CMakeLists.txt b/Tests/Core/Fitter/CMakeLists.txt
index 31df73253cbc42cbe6e1c8845a3ac6a35fb4c40f..059778e3994e88a8ea3c8ba03b4d04f09a1ec05f 100644
--- a/Tests/Core/Fitter/CMakeLists.txt
+++ b/Tests/Core/Fitter/CMakeLists.txt
@@ -4,6 +4,12 @@ target_link_libraries (KalmanFitterTests PRIVATE ActsCore ActsTestsCommonHelpers
add_test (NAME KalmanFitterUnitTests COMMAND KalmanFitterTests)
acts_add_test_to_cdash_project (PROJECT ACore TEST KalmanFitterUnitTests TARGETS KalmanFitterTests)
+add_executable (SequentialKalmanFitterTests SequentialKalmanFitterTests.cpp)
+target_include_directories (SequentialKalmanFitterTests PRIVATE ${Boost_INCLUDE_DIRS})
+target_link_libraries (SequentialKalmanFitterTests PRIVATE ActsCore ActsTestsCommonHelpers)
+add_test (NAME SequentialKalmanFitterUnitTests COMMAND SequentialKalmanFitterTests)
+acts_add_test_to_cdash_project (PROJECT ACore TEST SequentialKalmanFitterUnitTests TARGETS SequentialKalmanFitterTests)
+
add_executable (GainMatrixTests GainMatrixTests.cpp)
target_include_directories (GainMatrixTests PRIVATE ${Boost_INCLUDE_DIRS})
target_link_libraries (GainMatrixTests PRIVATE ActsCore)
diff --git a/Tests/Core/Fitter/KalmanFilterTestUtils.hpp b/Tests/Core/Fitter/KalmanFilterTestUtils.hpp
deleted file mode 100644
index a7ffd9beeb568b9909a7fbf26b5ddc2fe50091df..0000000000000000000000000000000000000000
--- a/Tests/Core/Fitter/KalmanFilterTestUtils.hpp
+++ /dev/null
@@ -1,152 +0,0 @@
-// This file is part of the Acts project.
-//
-// Copyright (C) 2017-2018 Acts project team
-//
-// 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 <cmath>
-#include <fstream>
-#include <iostream>
-#include <memory>
-#include <random>
-#include <vector>
-#include "Acts/Detector/TrackingGeometry.hpp"
-#include "Acts/EventData/Measurement.hpp"
-#include "Acts/Extrapolation/ExtrapolationCell.hpp"
-#include "Acts/Extrapolation/ExtrapolationEngine.hpp"
-#include "Acts/Extrapolation/IExtrapolationEngine.hpp"
-#include "Acts/Extrapolation/MaterialEffectsEngine.hpp"
-#include "Acts/Extrapolation/RungeKuttaEngine.hpp"
-#include "Acts/Extrapolation/StaticEngine.hpp"
-#include "Acts/Extrapolation/StaticNavigationEngine.hpp"
-#include "Acts/Fitter/KalmanFitter.hpp"
-#include "Acts/Fitter/KalmanUpdator.hpp"
-#include "Acts/MagneticField/ConstantBField.hpp"
-#include "Acts/Surfaces/PerigeeSurface.hpp"
-#include "Acts/Utilities/Definitions.hpp"
-#include "Acts/Utilities/Logger.hpp"
-#include "Acts/Utilities/Units.hpp"
-
-// shorthands
-using namespace Acts;
-using FitMeas_t = FittableMeasurement<long int>;
-template <ParID_t... pars>
-using Meas_t = Measurement<long int, pars...>;
-
-/// fit cache
-struct MyCache
-{
- std::unique_ptr<const KF::Step<long int>::JacobianMatrix> jacobian;
- std::unique_ptr<const BoundParameters> parameters;
-
- MyCache() = default;
- MyCache(const MyCache&) = delete;
- MyCache(MyCache&&) = default;
-};
-
-/// extrapolation wrapper
-class MyExtrapolator
-{
-public:
- MyExtrapolator(std::shared_ptr<const IExtrapolationEngine> exEngine
- = nullptr);
-
- MyCache
- operator()(const FitMeas_t& m, const TrackParameters& tp) const;
-
-private:
- std::shared_ptr<const IExtrapolationEngine> m_exEngine;
-};
-
-/// dummy class, returns measurement unchanged
-class NoCalibration
-{
-public:
- std::unique_ptr<const FitMeas_t>
- operator()(const FitMeas_t& m, const BoundParameters&) const;
-};
-
-class CacheGenerator
-{
-public:
- std::unique_ptr<KF::Step<long int>>
- operator()(MyCache m) const;
-};
-
-MyExtrapolator::MyExtrapolator(
- std::shared_ptr<const IExtrapolationEngine> exEngine)
- : m_exEngine(std::move(exEngine)){};
-
-/// wrapper around extrapolate call to exEngine, setting the right flags
-MyCache
-MyExtrapolator::operator()(const FitMeas_t& m, const TrackParameters& tp) const
-{
- auto exCell = std::make_unique<ExtrapolationCell<TrackParameters>>(tp);
- exCell->addConfigurationMode(ExtrapolationMode::CollectJacobians);
- (*exCell).pathLimit = 500;
- const Surface& sf = getSurface(m);
-
- m_exEngine->extrapolate(*exCell, &sf);
- MyCache c;
- auto j = exCell->extrapolationSteps.back().transportJacobian.release();
- c.jacobian.reset(new KF::Step<long int>::JacobianMatrix(*j));
- auto pars
- = static_cast<const BoundParameters*>(exCell->leadParameters->clone());
- c.parameters.reset(pars);
-
- return c;
-};
-
-std::unique_ptr<const FitMeas_t>
-NoCalibration::operator()(const FitMeas_t& m, const BoundParameters&) const
-{
- return std::make_unique<const FitMeas_t>(m);
-};
-
-std::unique_ptr<KF::Step<long int>>
-CacheGenerator::operator()(MyCache m) const
-{
- auto step = std::make_unique<KF::Step<long int>>();
- step->setPredictedState(std::move(m.parameters));
- step->setJacobian(std::move(m.jacobian));
-
- return step;
-};
-
-/// set up extrapolation
-std::shared_ptr<IExtrapolationEngine>
-initExtrapolator(const std::shared_ptr<const TrackingGeometry>& geo)
-{
- auto propConfig = RungeKuttaEngine<>::Config();
- propConfig.fieldService
- = std::make_shared<ConstantBField>(0, 0, 2 * Acts::units::_T);
- auto propEngine = std::make_shared<RungeKuttaEngine<>>(propConfig);
-
- auto matConfig = MaterialEffectsEngine::Config();
- auto materialEngine = std::make_shared<MaterialEffectsEngine>(matConfig);
-
- auto navConfig = StaticNavigationEngine::Config();
- navConfig.propagationEngine = propEngine;
- navConfig.materialEffectsEngine = materialEngine;
- navConfig.trackingGeometry = geo;
- auto navEngine = std::make_shared<StaticNavigationEngine>(navConfig);
-
- auto statConfig = StaticEngine::Config();
- statConfig.propagationEngine = propEngine;
- statConfig.navigationEngine = navEngine;
- statConfig.materialEffectsEngine = materialEngine;
- auto statEngine = std::make_shared<StaticEngine>(statConfig);
-
- auto exEngineConfig = ExtrapolationEngine::Config();
- exEngineConfig.trackingGeometry = geo;
- exEngineConfig.propagationEngine = propEngine;
- exEngineConfig.navigationEngine = navEngine;
- exEngineConfig.extrapolationEngines = {statEngine};
- auto exEngine = std::make_shared<ExtrapolationEngine>(exEngineConfig);
- exEngine->setLogger(getDefaultLogger("ExtrapolationEngine", Logging::INFO));
-
- return exEngine;
-};
\ No newline at end of file
diff --git a/Tests/Core/Fitter/KalmanFitterTests.cpp b/Tests/Core/Fitter/KalmanFitterTests.cpp
index dd5655cfae3a2ae275021cb2f1a0c268dbd3f38d..8dcc2094e9932a971f196a69570294e1e6e84f32 100644
--- a/Tests/Core/Fitter/KalmanFitterTests.cpp
+++ b/Tests/Core/Fitter/KalmanFitterTests.cpp
@@ -30,250 +30,47 @@
#include "Acts/Propagator/detail/StandardAborters.hpp"
#include "Acts/Surfaces/Surface.hpp"
#include "Acts/Tests/CommonHelpers/CubicTrackingGeometry.hpp"
+#include "Acts/Tests/CommonHelpers/MeasurementCreatorHelper.hpp"
#include "Acts/Utilities/BinningType.hpp"
#include "Acts/Utilities/Definitions.hpp"
#include "Acts/Utilities/GeometryID.hpp"
namespace Acts {
namespace Test {
-
- // A few initialisations and definitionas
- using Identifier = GeometryID;
- using Jacobian = BoundParameters::CovMatrix_t;
-
- using TrackState = TrackState<Identifier, BoundParameters>;
- using Resolution = std::pair<ParID_t, double>;
- using ElementResolution = std::vector<Resolution>;
- using VolumeResolution = std::map<geo_id_value, ElementResolution>;
- using DetectorResolution = std::map<geo_id_value, VolumeResolution>;
-
- using DebugOutput = detail::DebugOutputActor;
-
- std::normal_distribution<double> gauss(0., 1.);
- std::default_random_engine generator(42);
-
- ActsSymMatrixD<1> cov1D;
- ActsSymMatrixD<2> cov2D;
-
- bool debugMode = false;
-
- /// @brief This struct creates FittableMeasurements on the
- /// detector surfaces, according to the given smearing xxparameters
- ///
- struct MeasurementCreator
- {
- /// @brief Constructor
- MeasurementCreator() = default;
-
- /// The detector resolution
- DetectorResolution detectorResolution;
-
- using result_type = std::vector<TrackState>;
-
- /// @brief Operater that is callable by an ActionList. The function collects
- /// the surfaces
- ///
- /// @tparam propagator_state_t Type of the propagator state
- /// @param [in] state State of the propagator
- /// @param [out] result Vector of matching surfaces
- template <typename propagator_state_t>
- void
- operator()(propagator_state_t& state, result_type& result) const
- {
- // monitor the current surface
- auto surface = state.navigation.currentSurface;
- if (surface and surface->associatedDetectorElement()) {
- auto geoID = surface->geoID();
- geo_id_value volumeID = geoID.value(GeometryID::volume_mask);
- geo_id_value layerID = geoID.value(GeometryID::layer_mask);
- // find volume and layer information for this
- auto vResolution = detectorResolution.find(volumeID);
- if (vResolution != detectorResolution.end()) {
- // find layer resolutions
- auto lResolution = vResolution->second.find(layerID);
- if (lResolution != vResolution->second.end()) {
- // Apply global to local
- Acts::Vector2D lPos;
- surface->globalToLocal(
- state.stepping.position(), state.stepping.direction(), lPos);
- if (lResolution->second.size() == 1) {
- double sp = lResolution->second[0].second;
- cov1D << sp * sp;
- double dp = sp * gauss(generator);
- if (lResolution->second[0].first == eLOC_0) {
- // push back & move a LOC_0 measurement
- Measurement<Identifier, eLOC_0> m0(
- surface->getSharedPtr(), geoID, cov1D, lPos[eLOC_0] + dp);
- result.push_back(TrackState(std::move(m0)));
- } else {
- // push back & move a LOC_1 measurement
- Measurement<Identifier, eLOC_1> m1(
- surface->getSharedPtr(), geoID, cov1D, lPos[eLOC_1] + dp);
- result.push_back(TrackState(std::move(m1)));
- }
- } else if (lResolution->second.size() == 2) {
- // Create the measurment and move it
- double sx = lResolution->second[eLOC_0].second;
- double sy = lResolution->second[eLOC_1].second;
- cov2D << sx * sx, 0., 0., sy * sy;
- double dx = sx * gauss(generator);
- double dy = sy * gauss(generator);
- // push back & move a LOC_0, LOC_1 measurement
- Measurement<Identifier, eLOC_0, eLOC_1> m01(
- surface->getSharedPtr(),
- geoID,
- cov2D,
- lPos[eLOC_0] + dx,
- lPos[eLOC_1] + dy);
- result.push_back(TrackState(std::move(m01)));
- }
- }
- }
- }
- }
- };
-
- double dX, dY;
- Vector3D pos;
- Surface const* sur;
-
- ///
- /// @brief Simplified material interaction effect by pure gaussian deflection
- ///
- struct MaterialScattering
- {
- /// @brief Constructor
- MaterialScattering() = default;
-
- /// @brief Main action list call operator for the scattering on material
- ///
- /// @todo deal momentum in a gaussian way properly
- ///
- /// @tparam propagator_state_t State of the propagator
- /// @param [in] state State of the propagation
- template <typename propagator_state_t>
- void
- operator()(propagator_state_t& state) const
- {
- // Check if there is a surface with material and a covariance is existing
- if (state.navigation.currentSurface
- && state.navigation.currentSurface->associatedMaterial()
- && state.stepping.cov != ActsSymMatrixD<5>::Zero()) {
- // Sample angles
- std::normal_distribution<double> scatterAngle(
- 0., 0.017); //< \approx 1 degree
- double dPhi = scatterAngle(generator), dTheta = scatterAngle(generator);
-
- // Update the covariance
- state.stepping.cov(ePHI, ePHI) += dPhi * dPhi;
- state.stepping.cov(eTHETA, eTHETA) += dTheta * dTheta;
-
- // Update the angles
- double theta = std::acos(state.stepping.direction().z());
- double phi = std::atan2(state.stepping.direction().y(),
- state.stepping.direction().x());
-
- state.stepping.update(
- state.stepping.position(),
- {std::sin(theta + dTheta) * std::cos(phi + dPhi),
- std::sin(theta + dTheta) * std::sin(phi + dPhi),
- std::cos(theta + dTheta)},
- std::max(state.stepping.p
- - std::abs(gauss(generator)) * units::_MeV,
- 0.));
- }
- }
- };
-
///
/// @brief Unit test for Kalman fitter with measurements along the x-axis
///
BOOST_AUTO_TEST_CASE(kalman_fitter_zero_field)
{
- // Build detector
- CubicTrackingGeometry cGeometry;
- auto detector = cGeometry();
+ CubicTrackingGeometry geom;
+ auto detector = geom();
- // Build navigator for the measurement creatoin
- Navigator mNavigator(detector);
- mNavigator.resolvePassive = false;
- mNavigator.resolveMaterial = true;
- mNavigator.resolveSensitive = true;
-
- // Use straingt line stepper to create the measurements
StraightLineStepper mStepper;
+ Vector3D startPosition(-3. * units::_m, 0., 0.);
+ Vector3D startMomentum(1. * units::_GeV, 0., 0);
- // Define the measurement propagator
- using MeasurementPropagator = Propagator<StraightLineStepper, Navigator>;
-
- // Build propagator for the measurement creation
- MeasurementPropagator mPropagator(mStepper, mNavigator);
- Vector3D mPos(-3. * units::_m, 0., 0.), mMom(1. * units::_GeV, 0., 0);
- SingleCurvilinearTrackParameters<NeutralPolicy> mStart(nullptr, mPos, mMom);
-
- // Create action list for the measurement creation
- using MeasurementActions = ActionList<MeasurementCreator, DebugOutput>;
- using MeasurementAborters = AbortList<detail::EndOfWorldReached>;
-
- auto pixelResX = Resolution(eLOC_0, 25. * units::_um);
- auto pixelResY = Resolution(eLOC_1, 50. * units::_um);
- auto stripResX = Resolution(eLOC_0, 100. * units::_um);
- auto stripResY = Resolution(eLOC_1, 150. * units::_um);
-
- ElementResolution pixelElementRes = {pixelResX, pixelResY};
- ElementResolution stripElementResI = {stripResX};
- ElementResolution stripElementResO = {stripResY};
-
- VolumeResolution pixelVolumeRes;
- pixelVolumeRes[2] = pixelElementRes;
- pixelVolumeRes[4] = pixelElementRes;
-
- VolumeResolution stripVolumeRes;
- stripVolumeRes[2] = stripElementResI;
- stripVolumeRes[4] = stripElementResO;
- stripVolumeRes[6] = stripElementResI;
- stripVolumeRes[8] = stripElementResO;
-
- DetectorResolution detRes;
- detRes[2] = pixelVolumeRes;
- detRes[3] = stripVolumeRes;
-
- // Set options for propagator
- PropagatorOptions<MeasurementActions, MeasurementAborters> mOptions;
- mOptions.debug = debugMode;
- auto& mCreator = mOptions.actionList.get<MeasurementCreator>();
- mCreator.detectorResolution = detRes;
-
- // Launch and collect - the measurements
- auto mResult = mPropagator.propagate(mStart, mOptions);
- if (debugMode) {
- const auto debugString
- = mResult.template get<DebugOutput::result_type>().debugString;
- std::cout << ">>>> Measurement creation: " << std::endl;
- std::cout << debugString;
- }
-
- auto measurements = mResult.template get<MeasurementCreator::result_type>();
+ auto measurements = createMeasurementsOnCubicDetector(
+ detector, mStepper, startPosition, startMomentum)
+ .trackStates;
BOOST_TEST(measurements.size() == 6);
// The KalmanFitter - we use the eigen stepper for covariance transport
- // Build navigator for the measurement creatoin
+
Navigator rNavigator(detector);
rNavigator.resolvePassive = false;
rNavigator.resolveMaterial = true;
rNavigator.resolveSensitive = true;
// Configure propagation with deactivated B-field
- ConstantBField bField(Vector3D(0., 0., 0.));
using RecoStepper = EigenStepper<ConstantBField>;
- RecoStepper rStepper(bField);
+ RecoStepper rStepper;
using RecoPropagator = Propagator<RecoStepper, Navigator>;
RecoPropagator rPropagator(rStepper, rNavigator);
// Set initial parameters for the particle track
ActsSymMatrixD<5> cov;
- cov << 1000. * units::_um, 0., 0., 0., 0., 0., 1000. * units::_um, 0., 0.,
- 0., 0., 0., 0.05, 0., 0., 0., 0., 0., 0.05, 0., 0., 0., 0., 0., 0.01;
+ cov << 1. * units::_um, 0., 0., 0., 0., 0., 1. * units::_um, 0., 0., 0., 0.,
+ 0., 0.001, 0., 0., 0., 0., 0., 0.001, 0., 0., 0., 0., 0., 0.001;
auto covPtr = std::make_unique<const ActsSymMatrixD<5>>(cov);
diff --git a/Tests/Core/Fitter/SequentialKalmanFitterTests.cpp b/Tests/Core/Fitter/SequentialKalmanFitterTests.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..0623d247b6ac2375a5afcffe74d87ea17c5eebc0
--- /dev/null
+++ b/Tests/Core/Fitter/SequentialKalmanFitterTests.cpp
@@ -0,0 +1,120 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2016-2018 Acts project team
+//
+// 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/.
+
+#define BOOST_TEST_MODULE SequentialKalmanFitter Tests
+
+#include <boost/test/included/unit_test.hpp>
+
+#include <algorithm>
+#include <math.h>
+#include <random>
+#include <vector>
+#include "Acts/Detector/TrackingGeometry.hpp"
+#include "Acts/EventData/Measurement.hpp"
+#include "Acts/EventData/TrackParameters.hpp"
+#include "Acts/EventData/TrackState.hpp"
+#include "Acts/Extrapolator/Navigator.hpp"
+#include "Acts/Fitter/GainMatrixSmoother.hpp"
+#include "Acts/Fitter/GainMatrixUpdator.hpp"
+#include "Acts/Fitter/SequentialKalmanFitter.hpp"
+#include "Acts/MagneticField/ConstantBField.hpp"
+#include "Acts/Propagator/EigenStepper.hpp"
+#include "Acts/Propagator/Propagator.hpp"
+#include "Acts/Tests/CommonHelpers/CubicTrackingGeometry.hpp"
+#include "Acts/Tests/CommonHelpers/MeasurementCreatorHelper.hpp"
+
+namespace Acts {
+namespace Test {
+
+ struct HitSelector
+ {
+ public:
+ HitSelector(std::map<const Surface*,
+ std::vector<FittableMeasurement<Identifier>>> map)
+ : m_measurementMap(std::move(map))
+ {
+ }
+
+ template <typename propagator_state_t>
+ std::vector<TrackState>
+ operator()(const Surface& surface, const propagator_state_t& /*state*/)
+ {
+ const auto& measurements = m_measurementMap[&surface];
+ if (measurements.empty()) {
+ return {};
+ }
+
+ std::vector<TrackState> trackStates;
+ for (const auto& measurement : measurements) {
+ trackStates.emplace_back(measurement);
+ }
+
+ return trackStates;
+ }
+
+ private:
+ std::map<const Surface*, std::vector<FittableMeasurement<Identifier>>>
+ m_measurementMap;
+ };
+
+ ///
+ /// @brief Unit test for Kalman fitter with measurements along the x-axis
+ ///
+ BOOST_AUTO_TEST_CASE(sequential_kalman_fitter_zero_field)
+ {
+ CubicTrackingGeometry geom;
+ auto detector = geom();
+
+ StraightLineStepper mStepper;
+ Vector3D startPosition(-3. * units::_m, 0., 0.);
+ Vector3D startMomentum(1. * units::_GeV, 0., 0);
+
+ auto measurementMap
+ = createMeasurementsOnCubicDetector(
+ detector, mStepper, startPosition, startMomentum, true)
+ .measurementMap;
+
+ // The KalmanFitter - we use the eigen stepper for covariance transport
+ // Build navigator for the measurement creatoin
+ Navigator rNavigator(detector);
+ rNavigator.resolvePassive = false;
+ rNavigator.resolveMaterial = true;
+ rNavigator.resolveSensitive = true;
+
+ // Configure propagation with deactivated B-field
+ using RecoStepper = EigenStepper<ConstantBField>;
+ RecoStepper rStepper;
+ using RecoPropagator = Propagator<RecoStepper, Navigator>;
+ RecoPropagator rPropagator(rStepper, rNavigator);
+
+ // Set initial parameters for the particle track
+ ActsSymMatrixD<5> cov;
+ cov << 1. * units::_um, 0., 0., 0., 0., 0., 1. * units::_um, 0., 0., 0., 0.,
+ 0., 0.001, 0., 0., 0., 0., 0., 0.001, 0., 0., 0., 0., 0., 0.001;
+
+ auto covPtr = std::make_unique<const ActsSymMatrixD<5>>(cov);
+
+ SingleCurvilinearTrackParameters<ChargedPolicy> rStart(
+ std::move(covPtr), startPosition, startMomentum, 1.);
+
+ using Updator = GainMatrixUpdator<BoundParameters>;
+ using KalmanFitter = SequentialKalmanFitter<RecoPropagator,
+ TrackState,
+ HitSelector,
+ Updator>;
+
+ KalmanFitter fitter(std::move(rPropagator));
+ auto testFit = fitter.fit(
+ rStart, std::make_shared<HitSelector>(std::move(measurementMap)));
+
+ BOOST_TEST(testFit.resultStateList.size() == 6);
+ BOOST_TEST(testFit.numberOfHoles == 0);
+ }
+
+} // namespace Test
+} // namespace Acts