From b79fe8ede9f3fd6fda4da5da47bb6412590ab12b Mon Sep 17 00:00:00 2001
From: Nils Braun <nils.braun@kit.edu>
Date: Thu, 17 Jan 2019 13:29:40 -0800
Subject: [PATCH 01/17] Added a first work-in-progress version of the SKF

---
 .../Acts/Fitter/SequentialKalmanFitter.hpp    | 102 ++++++++++++++++++
 1 file changed, 102 insertions(+)
 create mode 100644 Core/include/Acts/Fitter/SequentialKalmanFitter.hpp

diff --git a/Core/include/Acts/Fitter/SequentialKalmanFitter.hpp b/Core/include/Acts/Fitter/SequentialKalmanFitter.hpp
new file mode 100644
index 000000000..6260763a2
--- /dev/null
+++ b/Core/include/Acts/Fitter/SequentialKalmanFitter.hpp
@@ -0,0 +1,102 @@
+// 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 <boost/variant.hpp>
+#include <memory>
+#include "Acts/EventData/Measurement.hpp"
+#include "Acts/EventData/TrackParameters.hpp"
+#include "Acts/Fitter/detail/VoidKalmanComponents.hpp"
+#include "Acts/Propagator/AbortList.hpp"
+#include "Acts/Propagator/ActionList.hpp"
+#include "Acts/Propagator/Propagator.hpp"
+#include "Acts/Propagator/detail/ConstrainedStep.hpp"
+#include "Acts/Propagator/detail/StandardAborters.hpp"
+#include "Acts/Utilities/Definitions.hpp"
+
+namespace Acts {
+
+    class MapBasedHitLookup {
+    };
+
+    template<typename propagator_t,
+            typename track_state_t,
+            typename updator_t          = VoidKalmanUpdator,
+            typename measurement_selector_t = MapBasedHitLookup,
+            typename calibrator_t       = VoidKalmanComponents>
+    class SequentialKalmanFitter {
+    public:
+        /// Default constructor is deleted
+        SequentialKalmanFitter() = delete;
+
+        /// Constructor from arguments
+        SequentialKalmanFitter(propagator_t pPropagator) : m_propagator(std::move(pPropagator)) {
+        }
+
+        template<typename parameters_t>
+        auto
+        fit(const parameters_t &sParameters,
+            const measurement_selector_t &pHitSelector) 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.hitSelectorPtr = &pHitSelector;
+
+            // 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:
+        propagator_t m_propagator;
+
+        class KalmanActor {
+        public:
+            KalmanActor(updator_t pUpdator = updator_t(), calibrator_t pCalibrator = calibrator_t()) :
+                    m_updator(std::move(pUpdator)), m_calibrator(std::move(pCalibrator)) {
+            }
+
+            struct this_result {
+                // Move the result into the fitted states
+                std::vector<track_state_t> fittedStates = {};
+            };
+
+            using result_type = this_result;
+
+            template<typename propagator_state_t>
+            void
+            operator()(propagator_state_t &state, result_type &result) const {
+                auto surface = state.navigation.currentSurface;
+                if (not surface) {
+                    return;
+                }
+
+                //auto hit = m_measurementSelector(surface, state);
+                std::cout << "Hello!" << std::endl;
+            }
+
+        private:
+            measurement_selector_t m_measurementSelector;
+            updator_t m_updator;
+            calibrator_t m_calibrator;
+        };
+    };
+
+}  // namespace Acts
-- 
GitLab


From c8c30d7d637c935e354907c53ebd325bcf470c0e Mon Sep 17 00:00:00 2001
From: Nils Braun <nils.braun@kit.edu>
Date: Thu, 17 Jan 2019 13:29:58 -0800
Subject: [PATCH 02/17] Started adding a test for the SKF by copying the Kalman
 example

---
 Tests/Core/Fitter/CMakeLists.txt              |   6 +
 .../Fitter/SequentialKalmanFitterTests.cpp    | 392 ++++++++++++++++++
 2 files changed, 398 insertions(+)
 create mode 100644 Tests/Core/Fitter/SequentialKalmanFitterTests.cpp

diff --git a/Tests/Core/Fitter/CMakeLists.txt b/Tests/Core/Fitter/CMakeLists.txt
index 31df73253..059778e39 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/SequentialKalmanFitterTests.cpp b/Tests/Core/Fitter/SequentialKalmanFitterTests.cpp
new file mode 100644
index 000000000..8fcbc08fe
--- /dev/null
+++ b/Tests/Core/Fitter/SequentialKalmanFitterTests.cpp
@@ -0,0 +1,392 @@
+// 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 KalmanFitter 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/Extrapolator/SurfaceCollector.hpp"
+#include "Acts/Fitter/GainMatrixSmoother.hpp"
+#include "Acts/Fitter/GainMatrixUpdator.hpp"
+#include "Acts/Fitter/KalmanFitter.hpp"
+#include "Acts/MagneticField/ConstantBField.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/Tests/CommonHelpers/CubicTrackingGeometry.hpp"
+#include "Acts/Utilities/BinningType.hpp"
+#include "Acts/Utilities/Definitions.hpp"
+#include "Acts/Utilities/GeometryID.hpp"
+#include "Acts/Fitter/SequentialKalmanFitter.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::mt19937 generator(42);
+        std::uniform_real_distribution<double> uniform_dist(0., 1.);
+        std::uniform_int_distribution<int> numberOfMeasurements(1, 5);
+
+        ActsSymMatrixD<1> cov1D;
+        ActsSymMatrixD<2> cov2D;
+
+        bool debugMode = true;
+
+        /// @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::multimap<const Surface *, FittableMeasurement<Identifier>>;
+
+            /// @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()) {
+
+                    //random number to determine how many extra measurements on this surface
+
+                    int nMeasurements = numberOfMeasurements(generator);
+
+                    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);
+
+                        for (int meas = 0; meas < nMeasurements; meas++) {
+
+                            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 (meas != 0) {
+                                        dp *= 5;
+                                    }
+                                    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.insert(std::make_pair(surface, m0));
+                                    } else {
+                                        // push back & move a LOC_1 measurement
+                                        Measurement<Identifier, eLOC_1> m1(
+                                                surface->getSharedPtr(), geoID, cov1D, lPos[eLOC_1] + dp);
+                                        result.insert(std::make_pair(surface, 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.insert(std::make_pair(surface, 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.));
+                }
+            }
+        };
+
+        struct HitSelector {
+
+        };
+
+        ///
+        /// @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();
+
+            // 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;
+
+            // 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);
+
+            //OK, let's make this a vector of start positions instead
+            std::vector<SingleCurvilinearTrackParameters<NeutralPolicy>> ParamsVector;
+            int nTracks = 6;
+            for (int i = 0; i < nTracks; i++) {
+                double xsize = 1;//get dimensions of first surface...
+                double ysize = 1;
+                Vector3D pos1(-3, uniform_dist(generator) * xsize, uniform_dist(generator) * ysize), mom1(
+                        1. * units::_GeV, 0., 0);
+                SingleCurvilinearTrackParameters<NeutralPolicy> start1(nullptr, pos1, mom1);
+                ParamsVector.push_back(start1);
+            }
+
+
+            // 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 meas_map = mResult.template get<MeasurementCreator::result_type>();
+
+            std::cout << "Number of Measurements: " << meas_map.size() << std::endl;
+
+            for(auto it = meas_map.begin(); it != meas_map.end(); it = meas_map.upper_bound(it->first)) {
+                const auto& measurement = *it;
+                std::cout << "New measurement on surface " << measurement.first->geoID().toString() << " and " << meas_map.count(measurement.first) << " in total on this surface" << std::endl;
+                const auto& allMeasurements = meas_map.equal_range(measurement.first);
+                for(auto measIt = allMeasurements.first; measIt != allMeasurements.second; ++measIt) {
+                    boost::apply_visitor([](const auto &concreteMeasurement) {
+                        std::cout << concreteMeasurement.parameters().transpose() << ", ";
+                    }, measIt->second);
+                    std::cout << std::endl;
+                }
+            }
+
+            //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);
+            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;
+
+            auto covPtr = std::make_unique<const ActsSymMatrixD<5>>(cov);
+
+            Vector3D rPos(-3. * units::_m,
+                          10. * units::_um * gauss(generator),
+                          100. * units::_um * gauss(generator));
+            Vector3D rMom(1. * units::_GeV,
+                          0.025 * units::_GeV * gauss(generator),
+                          0.025 * units::_GeV * gauss(generator));
+
+            SingleCurvilinearTrackParameters<ChargedPolicy> rStart(
+                    std::move(covPtr), rPos, rMom, 1.);
+
+            const Surface *rSurface = &rStart.referenceSurface();
+
+            using Updator      = GainMatrixUpdator<BoundParameters>;
+            using KalmanFitter = SequentialKalmanFitter<RecoPropagator, TrackState, Updator, HitSelector>;
+            HitSelector hitSelector;
+
+
+            KalmanFitter fitter(std::move(rPropagator));
+            auto testFit = fitter.fit(rStart, hitSelector);
+
+            /*std::cout<<"do we get anything back? "<<testFit.fittedParameters.get()<<std::endl;
+
+            // Fit the track
+            auto fittedTrack      = kFitter.fit(measurements, rStart, rSurface);
+            auto fittedParameters = fittedTrack.fittedParameters.get();
+
+            // Make sure it is deterministic
+            auto fittedAgainTrack      = kFitter.fit(measurements, rStart, rSurface);
+            auto fittedAgainParameters = fittedAgainTrack.fittedParameters.get();
+
+            BOOST_TEST(fittedParameters.parameters().isApprox(
+                    fittedAgainParameters.parameters()));
+
+            // Change the order of the measurements
+            std::vector<TrackState> shuffledMeasurements = {measurements[3],
+                                                            measurements[2],
+                                                            measurements[1],
+                                                            measurements[4],
+                                                            measurements[5],
+                                                            measurements[0]};
+
+            // Make sure it works for shuffled measurements as well
+            auto fittedShuffledTrack
+                    = kFitter.fit(shuffledMeasurements, rStart, rSurface);
+            auto fittedShuffledParameters = fittedShuffledTrack.fittedParameters.get();
+
+            BOOST_TEST(fittedParameters.parameters().isApprox(
+                    fittedShuffledParameters.parameters()));
+
+            // Remove one measurement and find a hole
+            std::vector<TrackState> measurementsWithHole = {measurements[0],
+                                                            measurements[1],
+                                                            measurements[2],
+                                                            measurements[4],
+                                                            measurements[5]};
+
+            // Make sure it works for shuffled measurements as well
+            auto fittedWithHoleTrack
+                    = kFitter.fit(measurementsWithHole, rStart, rSurface);
+            auto fittedWithHoleParameters = fittedWithHoleTrack.fittedParameters.get();
+
+            // Count one hole
+            BOOST_TEST(fittedWithHoleTrack.missedActiveSurfaces.size() == 1);
+            // And the parameters should be different
+            BOOST_TEST(!fittedParameters.parameters().isApprox(
+                    fittedWithHoleParameters.parameters()));*/
+        }
+
+    }  // namespace Test
+}  // namespace Acts
-- 
GitLab


From fca21d49a01d1ea38f259c220f56802d218ef657 Mon Sep 17 00:00:00 2001
From: Nils Braun <nils.braun@kit.edu>
Date: Thu, 17 Jan 2019 15:37:13 -0800
Subject: [PATCH 03/17] Implemented the hit selection, bot choosing the first
 one every time

---
 .../Acts/Fitter/SequentialKalmanFitter.hpp    | 39 ++++++++---
 .../Fitter/SequentialKalmanFitterTests.cpp    | 69 ++++++++++++++-----
 2 files changed, 79 insertions(+), 29 deletions(-)

diff --git a/Core/include/Acts/Fitter/SequentialKalmanFitter.hpp b/Core/include/Acts/Fitter/SequentialKalmanFitter.hpp
index 6260763a2..f3dbe8e05 100644
--- a/Core/include/Acts/Fitter/SequentialKalmanFitter.hpp
+++ b/Core/include/Acts/Fitter/SequentialKalmanFitter.hpp
@@ -22,13 +22,10 @@
 
 namespace Acts {
 
-    class MapBasedHitLookup {
-    };
-
     template<typename propagator_t,
             typename track_state_t,
+            typename measurement_selector_t,
             typename updator_t          = VoidKalmanUpdator,
-            typename measurement_selector_t = MapBasedHitLookup,
             typename calibrator_t       = VoidKalmanComponents>
     class SequentialKalmanFitter {
     public:
@@ -42,7 +39,7 @@ namespace Acts {
         template<typename parameters_t>
         auto
         fit(const parameters_t &sParameters,
-            const measurement_selector_t &pHitSelector) const {
+            std::shared_ptr<measurement_selector_t> pHitSelector) const {
             // Create the ActionList and AbortList
             using KalmanResult = typename KalmanActor::result_type;
             using Actors       = ActionList<KalmanActor>;
@@ -51,7 +48,7 @@ namespace Acts {
             // Create relevant options for the propagation options
             PropagatorOptions <Actors, Aborters> kalmanOptions;
             auto &kalmanActor = kalmanOptions.actionList.template get<KalmanActor>();
-            //kalmanActor.hitSelectorPtr = &pHitSelector;
+            kalmanActor.m_measurementSelectorPtr = std::move(pHitSelector);
 
             // Run the fitter
             const auto &result
@@ -76,6 +73,7 @@ namespace Acts {
             struct this_result {
                 // Move the result into the fitted states
                 std::vector<track_state_t> fittedStates = {};
+                unsigned int numberOfHoles = 0;
             };
 
             using result_type = this_result;
@@ -83,17 +81,36 @@ namespace Acts {
             template<typename propagator_state_t>
             void
             operator()(propagator_state_t &state, result_type &result) const {
-                auto surface = state.navigation.currentSurface;
+                const Surface *surface = state.navigation.currentSurface;
                 if (not surface) {
                     return;
                 }
 
-                //auto hit = m_measurementSelector(surface, state);
-                std::cout << "Hello!" << std::endl;
+                std::vector<track_state_t> nextStates = (*m_measurementSelectorPtr)(*surface, state);
+                if (nextStates.empty()) {
+                    if (surface->associatedDetectorElement() != nullptr) {
+                        result.numberOfHoles++;
+                    }
+                    return;
+                }
+
+                track_state_t& chosenNextState = nextStates[0];
+
+
+                auto boundStates = state.stepping.boundState(*surface, true);
+                chosenNextState.parameter.predicted = std::get<0>(boundStates);
+                chosenNextState.parameter.jacobian = std::get<1>(boundStates);
+                chosenNextState.parameter.pathLength = std::get<2>(boundStates);
+
+                if (not m_updator(chosenNextState)) {
+                    return;
+                }
+
+                state.stepping.update(*chosenNextState.parameter.filtered);
+                result.fittedStates.push_back(std::move(chosenNextState));
             }
 
-        private:
-            measurement_selector_t m_measurementSelector;
+            std::shared_ptr<measurement_selector_t> m_measurementSelectorPtr;
             updator_t m_updator;
             calibrator_t m_calibrator;
         };
diff --git a/Tests/Core/Fitter/SequentialKalmanFitterTests.cpp b/Tests/Core/Fitter/SequentialKalmanFitterTests.cpp
index 8fcbc08fe..38248c637 100644
--- a/Tests/Core/Fitter/SequentialKalmanFitterTests.cpp
+++ b/Tests/Core/Fitter/SequentialKalmanFitterTests.cpp
@@ -71,7 +71,7 @@ namespace Acts {
             /// The detector resolution
             DetectorResolution detectorResolution;
 
-            using result_type = std::multimap<const Surface *, FittableMeasurement<Identifier>>;
+            using result_type = std::map<const Surface *, std::vector<FittableMeasurement<Identifier>>>;
 
             /// @brief Operater that is callable by an ActionList. The function collects
             /// the surfaces
@@ -110,19 +110,19 @@ namespace Acts {
                                     double sp = lResolution->second[0].second;
                                     cov1D << sp * sp;
                                     double dp = sp * gauss(generator);
-                                    if (meas != 0) {
+                                    /*if (meas != 0) {
                                         dp *= 5;
-                                    }
+                                    }*/
                                     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.insert(std::make_pair(surface, m0));
+                                        result[surface].push_back(m0);
                                     } else {
                                         // push back & move a LOC_1 measurement
                                         Measurement<Identifier, eLOC_1> m1(
                                                 surface->getSharedPtr(), geoID, cov1D, lPos[eLOC_1] + dp);
-                                        result.insert(std::make_pair(surface, m1));
+                                        result[surface].push_back(m1);
                                     }
                                 } else if (lResolution->second.size() == 2) {
                                     // Create the measurment and move it
@@ -138,7 +138,7 @@ namespace Acts {
                                             cov2D,
                                             lPos[eLOC_0] + dx,
                                             lPos[eLOC_1] + dy);
-                                    result.insert(std::make_pair(surface, m01));
+                                    result[surface].push_back(m01);
                                 }
                             }
                         }
@@ -198,7 +198,34 @@ namespace Acts {
         };
 
         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*/) {
+                std::cout << "I am on surface " << surface.geoID().toString() << std::endl;
+                const auto& measurements = m_measurementMap[&surface];
+                if(measurements.empty()) {
+                    return {};
+                }
+                std::cout << "I have found a measurement with ";
+                boost::apply_visitor([](const auto &concreteMeasurement) {
+                    std::cout << concreteMeasurement.parameters().transpose() << ", ";
+                }, measurements[0]);
+                std::cout << std::endl;
+
+                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;
         };
 
         ///
@@ -285,14 +312,14 @@ namespace Acts {
 
             std::cout << "Number of Measurements: " << meas_map.size() << std::endl;
 
-            for(auto it = meas_map.begin(); it != meas_map.end(); it = meas_map.upper_bound(it->first)) {
-                const auto& measurement = *it;
-                std::cout << "New measurement on surface " << measurement.first->geoID().toString() << " and " << meas_map.count(measurement.first) << " in total on this surface" << std::endl;
-                const auto& allMeasurements = meas_map.equal_range(measurement.first);
-                for(auto measIt = allMeasurements.first; measIt != allMeasurements.second; ++measIt) {
+            for(const auto& keyValue : meas_map) {
+                const Surface* surface = keyValue.first;
+                const auto& measurements = keyValue.second;
+                std::cout << "New measurement on surface " << surface->geoID().toString() << " and " << measurements.size() << " in total on this surface" << std::endl;
+                for(const auto& measurement : measurements) {
                     boost::apply_visitor([](const auto &concreteMeasurement) {
                         std::cout << concreteMeasurement.parameters().transpose() << ", ";
-                    }, measIt->second);
+                    }, measurement);
                     std::cout << std::endl;
                 }
             }
@@ -330,15 +357,21 @@ namespace Acts {
             SingleCurvilinearTrackParameters<ChargedPolicy> rStart(
                     std::move(covPtr), rPos, rMom, 1.);
 
-            const Surface *rSurface = &rStart.referenceSurface();
-
             using Updator      = GainMatrixUpdator<BoundParameters>;
-            using KalmanFitter = SequentialKalmanFitter<RecoPropagator, TrackState, Updator, HitSelector>;
-            HitSelector hitSelector;
-
+            using KalmanFitter = SequentialKalmanFitter<RecoPropagator, TrackState, HitSelector, Updator>;
 
             KalmanFitter fitter(std::move(rPropagator));
-            auto testFit = fitter.fit(rStart, hitSelector);
+            auto testFit = fitter.fit(rStart, std::make_shared<HitSelector>(std::move(meas_map)));
+
+            std::cout << "I have found " << testFit.fittedStates.size() << " hits and have skipped " << testFit.numberOfHoles <<  std::endl;
+            for(const auto& state : testFit.fittedStates) {
+                std::cout << "\ton surface " << state.referenceSurface().geoID().toString() << ": ";
+                std::cout << std::endl << *state.parameter.predicted << std::endl << *state.parameter.filtered << ", " << std::endl;
+                boost::apply_visitor([](const auto &concreteMeasurement) {
+                    std::cout << concreteMeasurement.parameters().transpose() << ", ";
+                }, *state.measurement.uncalibrated);
+                std::cout << std::endl;
+            }
 
             /*std::cout<<"do we get anything back? "<<testFit.fittedParameters.get()<<std::endl;
 
-- 
GitLab


From e728df665e0bbe59c5ac74f934eda11059c9da0a Mon Sep 17 00:00:00 2001
From: Nils Braun <nils.braun@kit.edu>
Date: Thu, 17 Jan 2019 16:32:48 -0800
Subject: [PATCH 04/17] Added a filtered chi^2 increment

---
 Core/include/Acts/EventData/TrackState.hpp | 273 +++++++++++----------
 1 file changed, 138 insertions(+), 135 deletions(-)

diff --git a/Core/include/Acts/EventData/TrackState.hpp b/Core/include/Acts/EventData/TrackState.hpp
index d4d861e2e..5b0dbf557 100644
--- a/Core/include/Acts/EventData/TrackState.hpp
+++ b/Core/include/Acts/EventData/TrackState.hpp
@@ -16,7 +16,7 @@
 
 namespace Acts {
 
-class Surface;
+    class Surface;
 
 /// @class TrackState
 ///
@@ -28,138 +28,141 @@ class Surface;
 ///
 /// @note the Surface is only stored as a pointer, i.e. it is
 /// assumed the surface lives longer than the TrackState
-template <typename identifier_t, typename parameters_t>
-class TrackState
-{
-
-public:
-  using Identifier = identifier_t;
-  using Parameters = parameters_t;
-  using Jacobian   = typename Parameters::CovMatrix_t;
-
-  /// Constructor from (uncalibrated) measurement
-  ///
-  /// @tparam measurement_t Type of the measurement
-  /// @param m The measurement object
-  TrackState(FittableMeasurement<identifier_t> m)
-  {
-    m_surface                = MeasurementHelpers::getSurface(m);
-    measurement.uncalibrated = std::move(m);
-  }
-
-  /// Constructor from parameters
-  ///
-  /// @tparam parameters_t Type of the predicted parameters
-  /// @param p The parameters object
-  TrackState(parameters_t p)
-  {
-    m_surface           = &p.referenceSurface();
-    parameter.predicted = std::move(p);
-  }
-
-  /// Virtual destructor
-  virtual ~TrackState() = default;
-
-  /// Copy constructor
-  ///
-  /// @param rhs is the source TrackState
-  TrackState(const TrackState& rhs)
-    : parameter(rhs.parameter)
-    , measurement(rhs.measurement)
-    , m_surface(rhs.m_surface)
-  {
-  }
-
-  /// Copy move constructor
-  ///
-  /// @param rhs is the source TrackState
-  TrackState(TrackState&& rhs)
-    : parameter(std::move(rhs.parameter))
-    , measurement(std::move(rhs.measurement))
-    , m_surface(std::move(rhs.m_surface))
-  {
-  }
-
-  /// Assignment operator
-  ///
-  /// @param rhs is the source TrackState
-  TrackState&
-  operator=(const TrackState& rhs)
-  {
-    parameter   = rhs.parameter;
-    measurement = rhs.measurement;
-    m_surface   = rhs.m_surface;
-    return (*this);
-  }
-
-  /// Assignment move operator
-  ///
-  /// @param rhs is the source TrackState
-  TrackState&
-  operator=(TrackState&& rhs)
-  {
-    parameter   = std::move(rhs.parameter);
-    measurement = std::move(rhs.measurement);
-    m_surface   = std::move(rhs.m_surface);
-    return (*this);
-  }
-
-  /// @brief return method for the surface
-  const Surface&
-  referenceSurface() const
-  {
-    return (*m_surface);
-  }
-
-  /// @brief number of Measured parameters, forwarded
-  /// @note This only returns a value if either of the measurements
-  ///       are set. If not, this returns boost::none
-  ///
-  /// @return number of measured parameters, or boost::none
-  boost::optional<size_t>
-  size()
-  {
-    if (this->measurement.uncalibrated) {
-      return MeasurementHelpers::getSize(*this->measurement.uncalibrated);
-    }
-    if (this->measurement.calibrated) {
-      return MeasurementHelpers::getSize(*this->measurement.calibrated);
-    }
-    return boost::none;
-  }
-
-  /// The parameter part
-  /// This is all the information that concerns the
-  /// the track parameterisation and the jacobian
-  /// It is enough to to run the track smoothing
-  struct
-  {
-    /// The predicted state
-    boost::optional<Parameters> predicted{boost::none};
-    /// The filtered state
-    boost::optional<Parameters> filtered{boost::none};
-    /// The smoothed state
-    boost::optional<Parameters> smoothed{boost::none};
-    /// The transport jacobian matrix
-    boost::optional<Jacobian> jacobian{boost::none};
-    /// The path length along the track - will help sorting
-    double pathLength = 0.;
-  } parameter;
-
-  /// @brief Nested measurement part
-  /// This is the uncalibrated and calibrated measurement
-  /// (in case the latter is different)
-  struct
-  {
-    /// The optional (uncalibrated) measurement
-    boost::optional<FittableMeasurement<identifier_t>> uncalibrated{
-        boost::none};
-    /// The optional calibrabed measurement
-    boost::optional<FittableMeasurement<identifier_t>> calibrated{boost::none};
-  } measurement;
-
-private:
-  /// The surface of this TrackState
-  const Surface* m_surface = nullptr;
-};
+    template <typename identifier_t, typename parameters_t>
+    class TrackState
+    {
+
+    public:
+        using Identifier = identifier_t;
+        using Parameters = parameters_t;
+        using Jacobian   = typename Parameters::CovMatrix_t;
+
+        /// Constructor from (uncalibrated) measurement
+        ///
+        /// @tparam measurement_t Type of the measurement
+        /// @param m The measurement object
+        TrackState(FittableMeasurement<identifier_t> m)
+        {
+          m_surface                = MeasurementHelpers::getSurface(m);
+          measurement.uncalibrated = std::move(m);
+        }
+
+        /// Constructor from parameters
+        ///
+        /// @tparam parameters_t Type of the predicted parameters
+        /// @param p The parameters object
+        TrackState(parameters_t p)
+        {
+          m_surface           = &p.referenceSurface();
+          parameter.predicted = std::move(p);
+        }
+
+        /// Virtual destructor
+        virtual ~TrackState() = default;
+
+        /// Copy constructor
+        ///
+        /// @param rhs is the source TrackState
+        TrackState(const TrackState& rhs)
+                : parameter(rhs.parameter)
+                , measurement(rhs.measurement)
+                , m_surface(rhs.m_surface)
+        {
+        }
+
+        /// Copy move constructor
+        ///
+        /// @param rhs is the source TrackState
+        TrackState(TrackState&& rhs)
+                : parameter(std::move(rhs.parameter))
+                , measurement(std::move(rhs.measurement))
+                , m_surface(std::move(rhs.m_surface))
+        {
+        }
+
+        /// Assignment operator
+        ///
+        /// @param rhs is the source TrackState
+        TrackState&
+        operator=(const TrackState& rhs)
+        {
+          parameter   = rhs.parameter;
+          measurement = rhs.measurement;
+          m_surface   = rhs.m_surface;
+          return (*this);
+        }
+
+        /// Assignment move operator
+        ///
+        /// @param rhs is the source TrackState
+        TrackState&
+        operator=(TrackState&& rhs)
+        {
+          parameter   = std::move(rhs.parameter);
+          measurement = std::move(rhs.measurement);
+          m_surface   = std::move(rhs.m_surface);
+          return (*this);
+        }
+
+        /// @brief return method for the surface
+        const Surface&
+        referenceSurface() const
+        {
+          return (*m_surface);
+        }
+
+        /// @brief number of Measured parameters, forwarded
+        /// @note This only returns a value if either of the measurements
+        ///       are set. If not, this returns boost::none
+        ///
+        /// @return number of measured parameters, or boost::none
+        boost::optional<size_t>
+        size()
+        {
+          if (this->measurement.uncalibrated) {
+            return MeasurementHelpers::getSize(*this->measurement.uncalibrated);
+          }
+          if (this->measurement.calibrated) {
+            return MeasurementHelpers::getSize(*this->measurement.calibrated);
+          }
+          return boost::none;
+        }
+
+        /// The parameter part
+        /// This is all the information that concerns the
+        /// the track parameterisation and the jacobian
+        /// It is enough to to run the track smoothing
+        struct
+        {
+            /// The predicted state
+            boost::optional<Parameters> predicted{boost::none};
+            /// The filtered state
+            boost::optional<Parameters> filtered{boost::none};
+            /// The smoothed state
+            boost::optional<Parameters> smoothed{boost::none};
+            /// The transport jacobian matrix
+            boost::optional<Jacobian> jacobian{boost::none};
+            /// The path length along the track - will help sorting
+            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)
+        struct
+        {
+            /// The optional (uncalibrated) measurement
+            boost::optional<FittableMeasurement<identifier_t>> uncalibrated{
+                    boost::none};
+            /// The optional calibrabed measurement
+            boost::optional<FittableMeasurement<identifier_t>> calibrated{boost::none};
+        } measurement;
+
+    private:
+        /// The surface of this TrackState
+        const Surface* m_surface = nullptr;
+    };
 }
-- 
GitLab


From 3278a7b53d8cfb7cf8e76ff46c399d487663a74f Mon Sep 17 00:00:00 2001
From: Nils Braun <nils.braun@kit.edu>
Date: Thu, 17 Jan 2019 16:33:17 -0800
Subject: [PATCH 05/17] Fill the chi^2 increment

---
 .../include/Acts/Fitter/GainMatrixUpdator.hpp | 211 +++++++++---------
 1 file changed, 106 insertions(+), 105 deletions(-)

diff --git a/Core/include/Acts/Fitter/GainMatrixUpdator.hpp b/Core/include/Acts/Fitter/GainMatrixUpdator.hpp
index bacfc24b4..277eebece 100644
--- a/Core/include/Acts/Fitter/GainMatrixUpdator.hpp
+++ b/Core/include/Acts/Fitter/GainMatrixUpdator.hpp
@@ -25,110 +25,111 @@ namespace Acts {
 ///
 /// This is implemented as a boost vistor pattern for use of the
 /// boost variant container
-template <typename parameters_t, typename calibrator_t = VoidKalmanComponents>
-class GainMatrixUpdator
-{
-
-  using jacobian_t = typename parameters_t::CovMatrix_t;
-
-public:
-  /// Explicit constructor
-  ///
-  /// @param calibrator is the calibration struct/class that converts
-  /// uncalibrated measurements into calibrated ones
-  GainMatrixUpdator(calibrator_t calibrator = calibrator_t())
-    : m_mCalibrator(std::move(calibrator))
-  {
-  }
-
-  /// @brief Public call operator for the boost visitor pattern
-  ///
-  /// @tparam track_state_t Type of the track state for the update
-  ///
-  /// @param trackState the measured track state
-  ///
-  /// @return Bool indicating whether this update was 'successful'
-  /// @note Non-'successful' updates could be holes or outliers,
-  ///       which need to be treated differently in calling code.
-  template <typename track_state_t>
-  bool
-  operator()(track_state_t& trackState) const
-  {
-    using CovMatrix_t = typename parameters_t::CovMatrix_t;
-    using ParVector_t = typename parameters_t::ParVector_t;
-
-    // we should definitely have an uncalibrated measurement here
-    assert(trackState.measurement.uncalibrated);
-    // there should be no calibrated measurement
-    assert(!trackState.measurement.calibrated);
-    // we should have predicted state set
-    assert(trackState.parameter.predicted);
-    // filtring should not have happened yet
-    assert(!trackState.parameter.filtered);
-
-    // read-only prediction handle
-    const parameters_t& predicted = *trackState.parameter.predicted;
-
-    const CovMatrix_t& predicted_covariance = *predicted.covariance();
-
-    ParVector_t filtered_parameters;
-    CovMatrix_t filtered_covariance;
-
-    // we need to remove type-erasure on the measurement type
-    // to access its methods
-    boost::apply_visitor(
-        [&](const auto& uncalibrated) {
-          // type of measurement
-          using meas_t = typename std::remove_const<
-              typename std::remove_reference<decltype(uncalibrated)>::type>::
-              type;
-          // type of projection
-          using projection_t = typename meas_t::Projection_t;
-          // type of gain matrix (transposed projection)
-          using gain_matrix_t = ActsMatrixD<projection_t::ColsAtCompileTime,
-                                            projection_t::RowsAtCompileTime>;
-
-          // Calibrate the measurement
-          meas_t calibrated = m_mCalibrator(uncalibrated, predicted);
-
-          // Take the projector (measurement mapping function)
-          const projection_t& H = calibrated.projector();
-
-          // The Kalman gain matrix
-          gain_matrix_t K = predicted_covariance * H.transpose()
-              * (H * predicted_covariance * H.transpose()
-                 + calibrated.covariance())
-                    .inverse();
-
-          // filtered new parameters after update
-          filtered_parameters
-              = predicted.parameters() + K * calibrated.residual(predicted);
-
-          // updated covariance after filtering
-          filtered_covariance
-              = (CovMatrix_t::Identity() - K * H) * predicted_covariance;
-
-          // plug calibrated measurement back into track state
-          trackState.measurement.calibrated = std::move(calibrated);
-
-        },
-        *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;
-  }
-
-private:
-  /// The measurement calibrator
-  calibrator_t m_mCalibrator;
-};
+    template<typename parameters_t, typename calibrator_t = VoidKalmanComponents>
+    class GainMatrixUpdator {
+
+        using jacobian_t = typename parameters_t::CovMatrix_t;
+
+    public:
+        /// Explicit constructor
+        ///
+        /// @param calibrator is the calibration struct/class that converts
+        /// uncalibrated measurements into calibrated ones
+        GainMatrixUpdator(calibrator_t calibrator = calibrator_t())
+                : m_mCalibrator(std::move(calibrator)) {
+        }
+
+        /// @brief Public call operator for the boost visitor pattern
+        ///
+        /// @tparam track_state_t Type of the track state for the update
+        ///
+        /// @param trackState the measured track state
+        ///
+        /// @return Bool indicating whether this update was 'successful'
+        /// @note Non-'successful' updates could be holes or outliers,
+        ///       which need to be treated differently in calling code.
+        template<typename track_state_t>
+        bool
+        operator()(track_state_t &trackState) const {
+            using CovMatrix_t = typename parameters_t::CovMatrix_t;
+            using ParVector_t = typename parameters_t::ParVector_t;
+
+            // we should definitely have an uncalibrated measurement here
+            assert(trackState.measurement.uncalibrated);
+            // there should be no calibrated measurement
+            assert(!trackState.measurement.calibrated);
+            // we should have predicted state set
+            assert(trackState.parameter.predicted);
+            // filtring should not have happened yet
+            assert(!trackState.parameter.filtered);
+
+            // read-only prediction handle
+            const parameters_t &predicted = *trackState.parameter.predicted;
+
+            const CovMatrix_t &predicted_covariance = *predicted.covariance();
+
+            ParVector_t filtered_parameters;
+            CovMatrix_t filtered_covariance;
+
+            // we need to remove type-erasure on the measurement type
+            // to access its methods
+            boost::apply_visitor(
+                    [&](const auto &uncalibrated) {
+                        // type of measurement
+                        using meas_t = typename std::remove_const<
+                                typename std::remove_reference<decltype(uncalibrated)>::type>::
+                        type;
+                        // type of projection
+                        using projection_t = typename meas_t::Projection_t;
+                        // type of gain matrix (transposed projection)
+                        using gain_matrix_t = ActsMatrixD<projection_t::ColsAtCompileTime,
+                                projection_t::RowsAtCompileTime>;
+
+                        // Calibrate the measurement
+                        meas_t calibrated = m_mCalibrator(uncalibrated, predicted);
+
+                        // Take the projector (measurement mapping function)
+                        const projection_t &H = calibrated.projector();
+
+                        // The Kalman gain matrix
+                        gain_matrix_t K = predicted_covariance * H.transpose()
+                                          * (H * predicted_covariance * H.transpose()
+                                             + calibrated.covariance())
+                                                  .inverse();
+
+                        // filtered new parameters after update
+                        filtered_parameters
+                                = predicted.parameters() + K * calibrated.residual(predicted);
+
+                        // updated covariance after filtering
+                        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);
+
+            // always succeed, no outlier logic yet
+            return true;
+        }
+
+    private:
+        /// The measurement calibrator
+        calibrator_t m_mCalibrator;
+    };
 
 }  // namespace Acts
\ No newline at end of file
-- 
GitLab


From 0f62b613c15c65940292769eecdf95c1c389aa42 Mon Sep 17 00:00:00 2001
From: Nils Braun <nils.braun@kit.edu>
Date: Thu, 17 Jan 2019 16:35:50 -0800
Subject: [PATCH 06/17] Finalized the implementation with a chi^2 ordering

---
 .../Acts/Fitter/SequentialKalmanFitter.hpp    |  33 +++--
 .../Fitter/SequentialKalmanFitterTests.cpp    | 133 ++++++++++--------
 2 files changed, 94 insertions(+), 72 deletions(-)

diff --git a/Core/include/Acts/Fitter/SequentialKalmanFitter.hpp b/Core/include/Acts/Fitter/SequentialKalmanFitter.hpp
index f3dbe8e05..a874b89ea 100644
--- a/Core/include/Acts/Fitter/SequentialKalmanFitter.hpp
+++ b/Core/include/Acts/Fitter/SequentialKalmanFitter.hpp
@@ -94,20 +94,31 @@ namespace Acts {
                     return;
                 }
 
-                track_state_t& chosenNextState = nextStates[0];
-
-
+                std::cout << "Looking for the best among.." << std::endl;
                 auto boundStates = state.stepping.boundState(*surface, true);
-                chosenNextState.parameter.predicted = std::get<0>(boundStates);
-                chosenNextState.parameter.jacobian = std::get<1>(boundStates);
-                chosenNextState.parameter.pathLength = std::get<2>(boundStates);
-
-                if (not m_updator(chosenNextState)) {
-                    return;
+                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);
+                    std::cout << "chi2 = " << nextState.filteredChi2Increment << " ";
+
+                    if(chosenNextState == nullptr) {
+                        chosenNextState = &nextState;
+                        std::cout << "first one" << std::endl;
+                    } else {
+                        if(chosenNextState->filteredChi2Increment > nextState.filteredChi2Increment) {
+                            chosenNextState = &nextState;
+                            std::cout << "is better" << std::endl;
+                        }
+                    }
                 }
 
-                state.stepping.update(*chosenNextState.parameter.filtered);
-                result.fittedStates.push_back(std::move(chosenNextState));
+                assert(chosenNextState != nullptr);
+                state.stepping.update(*(chosenNextState->parameter.filtered));
+                result.fittedStates.push_back(std::move(*chosenNextState));
             }
 
             std::shared_ptr<measurement_selector_t> m_measurementSelectorPtr;
diff --git a/Tests/Core/Fitter/SequentialKalmanFitterTests.cpp b/Tests/Core/Fitter/SequentialKalmanFitterTests.cpp
index 38248c637..87d3647f1 100644
--- a/Tests/Core/Fitter/SequentialKalmanFitterTests.cpp
+++ b/Tests/Core/Fitter/SequentialKalmanFitterTests.cpp
@@ -40,7 +40,7 @@ namespace Acts {
     namespace Test {
 
         // A few initialisations and definitionas
-        using Identifier = GeometryID;
+        using Identifier = unsigned int;
         using Jacobian   = BoundParameters::CovMatrix_t;
 
         using TrackState         = TrackState<Identifier, BoundParameters>;
@@ -84,64 +84,68 @@ namespace Acts {
             operator()(propagator_state_t &state, result_type &result) const {
                 // monitor the current surface
                 auto surface = state.navigation.currentSurface;
-                if (surface and surface->associatedDetectorElement()) {
-
-                    //random number to determine how many extra measurements on this surface
-
-                    int nMeasurements = numberOfMeasurements(generator);
-
-                    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);
-
-                        for (int meas = 0; meas < nMeasurements; meas++) {
-
-                            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 (meas != 0) {
-                                        dp *= 5;
-                                    }*/
-                                    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[surface].push_back(m0);
-                                    } else {
-                                        // push back & move a LOC_1 measurement
-                                        Measurement<Identifier, eLOC_1> m1(
-                                                surface->getSharedPtr(), geoID, cov1D, lPos[eLOC_1] + dp);
-                                        result[surface].push_back(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[surface].push_back(m01);
-                                }
-                            }
+                if (surface == nullptr or not surface->associatedDetectorElement()) {
+                    return;
+                }
+                //random number to determine how many extra measurements on this surface
+                unsigned int nMeasurements = static_cast<unsigned int>(numberOfMeasurements(generator));
+
+                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;
+                }
+
+                for (unsigned int identifier = 0; identifier < nMeasurements; identifier++) {
+                    // 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 (identifier != 0) {
+                            dp = 5*sp + dp*10;
                         }
+                        if (lResolution->second[0].first == eLOC_0) {
+                            // push back & move a LOC_0 measurement
+                            Measurement<Identifier, eLOC_0> m0(
+                                    surface->getSharedPtr(), identifier, cov1D, lPos[eLOC_0] + dp);
+                            result[surface].push_back(m0);
+                        } else {
+                            // push back & move a LOC_1 measurement
+                            Measurement<Identifier, eLOC_1> m1(
+                                    surface->getSharedPtr(), identifier, cov1D, lPos[eLOC_1] + dp);
+                            result[surface].push_back(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);
+                        if (identifier != 0) {
+                            dx = 5*sx * 10*dx;
+                            dy = 5*sy * 10*dy;
+                        }
+                        // push back & move a LOC_0, LOC_1 measurement
+                        Measurement<Identifier, eLOC_0, eLOC_1> m01(
+                                surface->getSharedPtr(),
+                                identifier,
+                                cov2D,
+                                lPos[eLOC_0] + dx,
+                                lPos[eLOC_1] + dy);
+                        result[surface].push_back(m01);
                     }
                 }
             }
@@ -342,8 +346,11 @@ namespace Acts {
 
             // 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);
 
@@ -355,7 +362,7 @@ namespace Acts {
                           0.025 * units::_GeV * gauss(generator));
 
             SingleCurvilinearTrackParameters<ChargedPolicy> rStart(
-                    std::move(covPtr), rPos, rMom, 1.);
+                    std::move(covPtr), mPos, mMom, 1.);
 
             using Updator      = GainMatrixUpdator<BoundParameters>;
             using KalmanFitter = SequentialKalmanFitter<RecoPropagator, TrackState, HitSelector, Updator>;
@@ -363,11 +370,15 @@ namespace Acts {
             KalmanFitter fitter(std::move(rPropagator));
             auto testFit = fitter.fit(rStart, std::make_shared<HitSelector>(std::move(meas_map)));
 
+            BOOST_TEST(testFit.fittedStates.size() == 6);
+            BOOST_TEST(testFit.numberOfHoles == 0);
+
             std::cout << "I have found " << testFit.fittedStates.size() << " hits and have skipped " << testFit.numberOfHoles <<  std::endl;
             for(const auto& state : testFit.fittedStates) {
                 std::cout << "\ton surface " << state.referenceSurface().geoID().toString() << ": ";
                 std::cout << std::endl << *state.parameter.predicted << std::endl << *state.parameter.filtered << ", " << std::endl;
                 boost::apply_visitor([](const auto &concreteMeasurement) {
+                    BOOST_TEST(concreteMeasurement.identifier() == 0);
                     std::cout << concreteMeasurement.parameters().transpose() << ", ";
                 }, *state.measurement.uncalibrated);
                 std::cout << std::endl;
-- 
GitLab


From 70a7e9c04dfc0fab2d653a979f0325b5060da3d7 Mon Sep 17 00:00:00 2001
From: Nils Braun <nils.braun@kit.edu>
Date: Thu, 17 Jan 2019 17:19:53 -0800
Subject: [PATCH 07/17] Clang format changes

---
 Core/include/Acts/EventData/TrackState.hpp    | 276 +++---
 .../include/Acts/Fitter/GainMatrixUpdator.hpp | 217 ++---
 .../Acts/Fitter/SequentialKalmanFitter.hpp    | 215 ++---
 .../Fitter/SequentialKalmanFitterTests.cpp    | 822 +++++++++---------
 4 files changed, 788 insertions(+), 742 deletions(-)

diff --git a/Core/include/Acts/EventData/TrackState.hpp b/Core/include/Acts/EventData/TrackState.hpp
index 5b0dbf557..21fe35e45 100644
--- a/Core/include/Acts/EventData/TrackState.hpp
+++ b/Core/include/Acts/EventData/TrackState.hpp
@@ -16,7 +16,7 @@
 
 namespace Acts {
 
-    class Surface;
+class Surface;
 
 /// @class TrackState
 ///
@@ -28,141 +28,141 @@ namespace Acts {
 ///
 /// @note the Surface is only stored as a pointer, i.e. it is
 /// assumed the surface lives longer than the TrackState
-    template <typename identifier_t, typename parameters_t>
-    class TrackState
-    {
-
-    public:
-        using Identifier = identifier_t;
-        using Parameters = parameters_t;
-        using Jacobian   = typename Parameters::CovMatrix_t;
-
-        /// Constructor from (uncalibrated) measurement
-        ///
-        /// @tparam measurement_t Type of the measurement
-        /// @param m The measurement object
-        TrackState(FittableMeasurement<identifier_t> m)
-        {
-          m_surface                = MeasurementHelpers::getSurface(m);
-          measurement.uncalibrated = std::move(m);
-        }
-
-        /// Constructor from parameters
-        ///
-        /// @tparam parameters_t Type of the predicted parameters
-        /// @param p The parameters object
-        TrackState(parameters_t p)
-        {
-          m_surface           = &p.referenceSurface();
-          parameter.predicted = std::move(p);
-        }
-
-        /// Virtual destructor
-        virtual ~TrackState() = default;
-
-        /// Copy constructor
-        ///
-        /// @param rhs is the source TrackState
-        TrackState(const TrackState& rhs)
-                : parameter(rhs.parameter)
-                , measurement(rhs.measurement)
-                , m_surface(rhs.m_surface)
-        {
-        }
-
-        /// Copy move constructor
-        ///
-        /// @param rhs is the source TrackState
-        TrackState(TrackState&& rhs)
-                : parameter(std::move(rhs.parameter))
-                , measurement(std::move(rhs.measurement))
-                , m_surface(std::move(rhs.m_surface))
-        {
-        }
-
-        /// Assignment operator
-        ///
-        /// @param rhs is the source TrackState
-        TrackState&
-        operator=(const TrackState& rhs)
-        {
-          parameter   = rhs.parameter;
-          measurement = rhs.measurement;
-          m_surface   = rhs.m_surface;
-          return (*this);
-        }
-
-        /// Assignment move operator
-        ///
-        /// @param rhs is the source TrackState
-        TrackState&
-        operator=(TrackState&& rhs)
-        {
-          parameter   = std::move(rhs.parameter);
-          measurement = std::move(rhs.measurement);
-          m_surface   = std::move(rhs.m_surface);
-          return (*this);
-        }
-
-        /// @brief return method for the surface
-        const Surface&
-        referenceSurface() const
-        {
-          return (*m_surface);
-        }
-
-        /// @brief number of Measured parameters, forwarded
-        /// @note This only returns a value if either of the measurements
-        ///       are set. If not, this returns boost::none
-        ///
-        /// @return number of measured parameters, or boost::none
-        boost::optional<size_t>
-        size()
-        {
-          if (this->measurement.uncalibrated) {
-            return MeasurementHelpers::getSize(*this->measurement.uncalibrated);
-          }
-          if (this->measurement.calibrated) {
-            return MeasurementHelpers::getSize(*this->measurement.calibrated);
-          }
-          return boost::none;
-        }
-
-        /// The parameter part
-        /// This is all the information that concerns the
-        /// the track parameterisation and the jacobian
-        /// It is enough to to run the track smoothing
-        struct
-        {
-            /// The predicted state
-            boost::optional<Parameters> predicted{boost::none};
-            /// The filtered state
-            boost::optional<Parameters> filtered{boost::none};
-            /// The smoothed state
-            boost::optional<Parameters> smoothed{boost::none};
-            /// The transport jacobian matrix
-            boost::optional<Jacobian> jacobian{boost::none};
-            /// The path length along the track - will help sorting
-            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)
-        struct
-        {
-            /// The optional (uncalibrated) measurement
-            boost::optional<FittableMeasurement<identifier_t>> uncalibrated{
-                    boost::none};
-            /// The optional calibrabed measurement
-            boost::optional<FittableMeasurement<identifier_t>> calibrated{boost::none};
-        } measurement;
-
-    private:
-        /// The surface of this TrackState
-        const Surface* m_surface = nullptr;
-    };
+template <typename identifier_t, typename parameters_t>
+class TrackState
+{
+
+public:
+  using Identifier = identifier_t;
+  using Parameters = parameters_t;
+  using Jacobian   = typename Parameters::CovMatrix_t;
+
+  /// Constructor from (uncalibrated) measurement
+  ///
+  /// @tparam measurement_t Type of the measurement
+  /// @param m The measurement object
+  TrackState(FittableMeasurement<identifier_t> m)
+  {
+    m_surface                = MeasurementHelpers::getSurface(m);
+    measurement.uncalibrated = std::move(m);
+  }
+
+  /// Constructor from parameters
+  ///
+  /// @tparam parameters_t Type of the predicted parameters
+  /// @param p The parameters object
+  TrackState(parameters_t p)
+  {
+    m_surface           = &p.referenceSurface();
+    parameter.predicted = std::move(p);
+  }
+
+  /// Virtual destructor
+  virtual ~TrackState() = default;
+
+  /// Copy constructor
+  ///
+  /// @param rhs is the source TrackState
+  TrackState(const TrackState& rhs)
+    : parameter(rhs.parameter)
+    , measurement(rhs.measurement)
+    , m_surface(rhs.m_surface)
+  {
+  }
+
+  /// Copy move constructor
+  ///
+  /// @param rhs is the source TrackState
+  TrackState(TrackState&& rhs)
+    : parameter(std::move(rhs.parameter))
+    , measurement(std::move(rhs.measurement))
+    , m_surface(std::move(rhs.m_surface))
+  {
+  }
+
+  /// Assignment operator
+  ///
+  /// @param rhs is the source TrackState
+  TrackState&
+  operator=(const TrackState& rhs)
+  {
+    parameter   = rhs.parameter;
+    measurement = rhs.measurement;
+    m_surface   = rhs.m_surface;
+    return (*this);
+  }
+
+  /// Assignment move operator
+  ///
+  /// @param rhs is the source TrackState
+  TrackState&
+  operator=(TrackState&& rhs)
+  {
+    parameter   = std::move(rhs.parameter);
+    measurement = std::move(rhs.measurement);
+    m_surface   = std::move(rhs.m_surface);
+    return (*this);
+  }
+
+  /// @brief return method for the surface
+  const Surface&
+  referenceSurface() const
+  {
+    return (*m_surface);
+  }
+
+  /// @brief number of Measured parameters, forwarded
+  /// @note This only returns a value if either of the measurements
+  ///       are set. If not, this returns boost::none
+  ///
+  /// @return number of measured parameters, or boost::none
+  boost::optional<size_t>
+  size()
+  {
+    if (this->measurement.uncalibrated) {
+      return MeasurementHelpers::getSize(*this->measurement.uncalibrated);
+    }
+    if (this->measurement.calibrated) {
+      return MeasurementHelpers::getSize(*this->measurement.calibrated);
+    }
+    return boost::none;
+  }
+
+  /// The parameter part
+  /// This is all the information that concerns the
+  /// the track parameterisation and the jacobian
+  /// It is enough to to run the track smoothing
+  struct
+  {
+    /// The predicted state
+    boost::optional<Parameters> predicted{boost::none};
+    /// The filtered state
+    boost::optional<Parameters> filtered{boost::none};
+    /// The smoothed state
+    boost::optional<Parameters> smoothed{boost::none};
+    /// The transport jacobian matrix
+    boost::optional<Jacobian> jacobian{boost::none};
+    /// The path length along the track - will help sorting
+    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)
+  struct
+  {
+    /// The optional (uncalibrated) measurement
+    boost::optional<FittableMeasurement<identifier_t>> uncalibrated{
+        boost::none};
+    /// The optional calibrabed measurement
+    boost::optional<FittableMeasurement<identifier_t>> calibrated{boost::none};
+  } measurement;
+
+private:
+  /// The surface of this TrackState
+  const Surface* m_surface = nullptr;
+};
 }
diff --git a/Core/include/Acts/Fitter/GainMatrixUpdator.hpp b/Core/include/Acts/Fitter/GainMatrixUpdator.hpp
index 277eebece..a82742c24 100644
--- a/Core/include/Acts/Fitter/GainMatrixUpdator.hpp
+++ b/Core/include/Acts/Fitter/GainMatrixUpdator.hpp
@@ -25,111 +25,116 @@ namespace Acts {
 ///
 /// This is implemented as a boost vistor pattern for use of the
 /// boost variant container
-    template<typename parameters_t, typename calibrator_t = VoidKalmanComponents>
-    class GainMatrixUpdator {
-
-        using jacobian_t = typename parameters_t::CovMatrix_t;
-
-    public:
-        /// Explicit constructor
-        ///
-        /// @param calibrator is the calibration struct/class that converts
-        /// uncalibrated measurements into calibrated ones
-        GainMatrixUpdator(calibrator_t calibrator = calibrator_t())
-                : m_mCalibrator(std::move(calibrator)) {
-        }
-
-        /// @brief Public call operator for the boost visitor pattern
-        ///
-        /// @tparam track_state_t Type of the track state for the update
-        ///
-        /// @param trackState the measured track state
-        ///
-        /// @return Bool indicating whether this update was 'successful'
-        /// @note Non-'successful' updates could be holes or outliers,
-        ///       which need to be treated differently in calling code.
-        template<typename track_state_t>
-        bool
-        operator()(track_state_t &trackState) const {
-            using CovMatrix_t = typename parameters_t::CovMatrix_t;
-            using ParVector_t = typename parameters_t::ParVector_t;
-
-            // we should definitely have an uncalibrated measurement here
-            assert(trackState.measurement.uncalibrated);
-            // there should be no calibrated measurement
-            assert(!trackState.measurement.calibrated);
-            // we should have predicted state set
-            assert(trackState.parameter.predicted);
-            // filtring should not have happened yet
-            assert(!trackState.parameter.filtered);
-
-            // read-only prediction handle
-            const parameters_t &predicted = *trackState.parameter.predicted;
-
-            const CovMatrix_t &predicted_covariance = *predicted.covariance();
-
-            ParVector_t filtered_parameters;
-            CovMatrix_t filtered_covariance;
-
-            // we need to remove type-erasure on the measurement type
-            // to access its methods
-            boost::apply_visitor(
-                    [&](const auto &uncalibrated) {
-                        // type of measurement
-                        using meas_t = typename std::remove_const<
-                                typename std::remove_reference<decltype(uncalibrated)>::type>::
-                        type;
-                        // type of projection
-                        using projection_t = typename meas_t::Projection_t;
-                        // type of gain matrix (transposed projection)
-                        using gain_matrix_t = ActsMatrixD<projection_t::ColsAtCompileTime,
-                                projection_t::RowsAtCompileTime>;
-
-                        // Calibrate the measurement
-                        meas_t calibrated = m_mCalibrator(uncalibrated, predicted);
-
-                        // Take the projector (measurement mapping function)
-                        const projection_t &H = calibrated.projector();
-
-                        // The Kalman gain matrix
-                        gain_matrix_t K = predicted_covariance * H.transpose()
-                                          * (H * predicted_covariance * H.transpose()
-                                             + calibrated.covariance())
-                                                  .inverse();
-
-                        // filtered new parameters after update
-                        filtered_parameters
-                                = predicted.parameters() + K * calibrated.residual(predicted);
-
-                        // updated covariance after filtering
-                        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);
-
-            // always succeed, no outlier logic yet
-            return true;
-        }
-
-    private:
-        /// The measurement calibrator
-        calibrator_t m_mCalibrator;
-    };
+template <typename parameters_t, typename calibrator_t = VoidKalmanComponents>
+class GainMatrixUpdator
+{
+
+  using jacobian_t = typename parameters_t::CovMatrix_t;
+
+public:
+  /// Explicit constructor
+  ///
+  /// @param calibrator is the calibration struct/class that converts
+  /// uncalibrated measurements into calibrated ones
+  GainMatrixUpdator(calibrator_t calibrator = calibrator_t())
+    : m_mCalibrator(std::move(calibrator))
+  {
+  }
+
+  /// @brief Public call operator for the boost visitor pattern
+  ///
+  /// @tparam track_state_t Type of the track state for the update
+  ///
+  /// @param trackState the measured track state
+  ///
+  /// @return Bool indicating whether this update was 'successful'
+  /// @note Non-'successful' updates could be holes or outliers,
+  ///       which need to be treated differently in calling code.
+  template <typename track_state_t>
+  bool
+  operator()(track_state_t& trackState) const
+  {
+    using CovMatrix_t = typename parameters_t::CovMatrix_t;
+    using ParVector_t = typename parameters_t::ParVector_t;
+
+    // we should definitely have an uncalibrated measurement here
+    assert(trackState.measurement.uncalibrated);
+    // there should be no calibrated measurement
+    assert(!trackState.measurement.calibrated);
+    // we should have predicted state set
+    assert(trackState.parameter.predicted);
+    // filtring should not have happened yet
+    assert(!trackState.parameter.filtered);
+
+    // read-only prediction handle
+    const parameters_t& predicted = *trackState.parameter.predicted;
+
+    const CovMatrix_t& predicted_covariance = *predicted.covariance();
+
+    ParVector_t filtered_parameters;
+    CovMatrix_t filtered_covariance;
+
+    // we need to remove type-erasure on the measurement type
+    // to access its methods
+    boost::apply_visitor(
+        [&](const auto& uncalibrated) {
+          // type of measurement
+          using meas_t = typename std::remove_const<
+              typename std::remove_reference<decltype(uncalibrated)>::type>::
+              type;
+          // type of projection
+          using projection_t = typename meas_t::Projection_t;
+          // type of gain matrix (transposed projection)
+          using gain_matrix_t = ActsMatrixD<projection_t::ColsAtCompileTime,
+                                            projection_t::RowsAtCompileTime>;
+
+          // Calibrate the measurement
+          meas_t calibrated = m_mCalibrator(uncalibrated, predicted);
+
+          // Take the projector (measurement mapping function)
+          const projection_t& H = calibrated.projector();
+
+          // The Kalman gain matrix
+          gain_matrix_t K = predicted_covariance * H.transpose()
+              * (H * predicted_covariance * H.transpose()
+                 + calibrated.covariance())
+                    .inverse();
+
+          // filtered new parameters after update
+          filtered_parameters
+              = predicted.parameters() + K * calibrated.residual(predicted);
+
+          // updated covariance after filtering
+          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);
+
+    // always succeed, no outlier logic yet
+    return true;
+  }
+
+private:
+  /// The measurement calibrator
+  calibrator_t m_mCalibrator;
+};
 
 }  // namespace Acts
\ No newline at end of file
diff --git a/Core/include/Acts/Fitter/SequentialKalmanFitter.hpp b/Core/include/Acts/Fitter/SequentialKalmanFitter.hpp
index a874b89ea..670b0f06e 100644
--- a/Core/include/Acts/Fitter/SequentialKalmanFitter.hpp
+++ b/Core/include/Acts/Fitter/SequentialKalmanFitter.hpp
@@ -22,109 +22,120 @@
 
 namespace Acts {
 
-    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 from arguments
-        SequentialKalmanFitter(propagator_t pPropagator) : m_propagator(std::move(pPropagator)) {
-        }
+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 from arguments
+  SequentialKalmanFitter(propagator_t pPropagator)
+    : m_propagator(std::move(pPropagator))
+  {
+  }
+
+  template <typename parameters_t>
+  auto
+  fit(const parameters_t&                     sParameters,
+      std::shared_ptr<measurement_selector_t> pHitSelector) 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(pHitSelector);
+
+    // 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:
+  propagator_t m_propagator;
+
+  class KalmanActor
+  {
+  public:
+    KalmanActor(updator_t    pUpdator    = updator_t(),
+                calibrator_t pCalibrator = calibrator_t())
+      : m_updator(std::move(pUpdator)), m_calibrator(std::move(pCalibrator))
+    {
+    }
+
+    struct this_result
+    {
+      // Move the result into the fitted states
+      std::vector<track_state_t> fittedStates  = {};
+      unsigned int               numberOfHoles = 0;
+    };
 
-        template<typename parameters_t>
-        auto
-        fit(const parameters_t &sParameters,
-            std::shared_ptr<measurement_selector_t> pHitSelector) 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(pHitSelector);
-
-            // 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;
+    using result_type = this_result;
+
+    template <typename propagator_state_t>
+    void
+    operator()(propagator_state_t& state, result_type& result) const
+    {
+      const Surface* surface = state.navigation.currentSurface;
+      if (not surface) {
+        return;
+      }
+
+      std::vector<track_state_t> nextStates
+          = (*m_measurementSelectorPtr)(*surface, state);
+      if (nextStates.empty()) {
+        if (surface->associatedDetectorElement() != nullptr) {
+          result.numberOfHoles++;
         }
-
-    private:
-        propagator_t m_propagator;
-
-        class KalmanActor {
-        public:
-            KalmanActor(updator_t pUpdator = updator_t(), calibrator_t pCalibrator = calibrator_t()) :
-                    m_updator(std::move(pUpdator)), m_calibrator(std::move(pCalibrator)) {
-            }
-
-            struct this_result {
-                // Move the result into the fitted states
-                std::vector<track_state_t> fittedStates = {};
-                unsigned int numberOfHoles = 0;
-            };
-
-            using result_type = this_result;
-
-            template<typename propagator_state_t>
-            void
-            operator()(propagator_state_t &state, result_type &result) const {
-                const Surface *surface = state.navigation.currentSurface;
-                if (not surface) {
-                    return;
-                }
-
-                std::vector<track_state_t> nextStates = (*m_measurementSelectorPtr)(*surface, state);
-                if (nextStates.empty()) {
-                    if (surface->associatedDetectorElement() != nullptr) {
-                        result.numberOfHoles++;
-                    }
-                    return;
-                }
-
-                std::cout << "Looking for the best among.." << std::endl;
-                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);
-                    std::cout << "chi2 = " << nextState.filteredChi2Increment << " ";
-
-                    if(chosenNextState == nullptr) {
-                        chosenNextState = &nextState;
-                        std::cout << "first one" << std::endl;
-                    } else {
-                        if(chosenNextState->filteredChi2Increment > nextState.filteredChi2Increment) {
-                            chosenNextState = &nextState;
-                            std::cout << "is better" << std::endl;
-                        }
-                    }
-                }
-
-                assert(chosenNextState != nullptr);
-                state.stepping.update(*(chosenNextState->parameter.filtered));
-                result.fittedStates.push_back(std::move(*chosenNextState));
-            }
-
-            std::shared_ptr<measurement_selector_t> m_measurementSelectorPtr;
-            updator_t m_updator;
-            calibrator_t m_calibrator;
-        };
-    };
+        return;
+      }
+
+      std::cout << "Looking for the best among.." << std::endl;
+      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);
+        std::cout << "chi2 = " << nextState.filteredChi2Increment << " ";
+
+        if (chosenNextState == nullptr) {
+          chosenNextState = &nextState;
+          std::cout << "first one" << std::endl;
+        } else {
+          if (chosenNextState->filteredChi2Increment
+              > nextState.filteredChi2Increment) {
+            chosenNextState = &nextState;
+            std::cout << "is better" << std::endl;
+          }
+        }
+      }
+
+      assert(chosenNextState != nullptr);
+      state.stepping.update(*(chosenNextState->parameter.filtered));
+      result.fittedStates.push_back(std::move(*chosenNextState));
+    }
+
+    std::shared_ptr<measurement_selector_t> m_measurementSelectorPtr;
+    updator_t                               m_updator;
+    calibrator_t                            m_calibrator;
+  };
+};
 
 }  // namespace Acts
diff --git a/Tests/Core/Fitter/SequentialKalmanFitterTests.cpp b/Tests/Core/Fitter/SequentialKalmanFitterTests.cpp
index 87d3647f1..0e95daa34 100644
--- a/Tests/Core/Fitter/SequentialKalmanFitterTests.cpp
+++ b/Tests/Core/Fitter/SequentialKalmanFitterTests.cpp
@@ -23,6 +23,7 @@
 #include "Acts/Fitter/GainMatrixSmoother.hpp"
 #include "Acts/Fitter/GainMatrixUpdator.hpp"
 #include "Acts/Fitter/KalmanFitter.hpp"
+#include "Acts/Fitter/SequentialKalmanFitter.hpp"
 #include "Acts/MagneticField/ConstantBField.hpp"
 #include "Acts/Propagator/EigenStepper.hpp"
 #include "Acts/Propagator/Propagator.hpp"
@@ -34,403 +35,432 @@
 #include "Acts/Utilities/BinningType.hpp"
 #include "Acts/Utilities/Definitions.hpp"
 #include "Acts/Utilities/GeometryID.hpp"
-#include "Acts/Fitter/SequentialKalmanFitter.hpp"
 
 namespace Acts {
-    namespace Test {
-
-        // A few initialisations and definitionas
-        using Identifier = unsigned int;
-        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::mt19937 generator(42);
-        std::uniform_real_distribution<double> uniform_dist(0., 1.);
-        std::uniform_int_distribution<int> numberOfMeasurements(1, 5);
-
-        ActsSymMatrixD<1> cov1D;
-        ActsSymMatrixD<2> cov2D;
-
-        bool debugMode = true;
-
-        /// @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::map<const Surface *, std::vector<FittableMeasurement<Identifier>>>;
-
-            /// @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 == nullptr or not surface->associatedDetectorElement()) {
-                    return;
-                }
-                //random number to determine how many extra measurements on this surface
-                unsigned int nMeasurements = static_cast<unsigned int>(numberOfMeasurements(generator));
-
-                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;
-                }
-
-                for (unsigned int identifier = 0; identifier < nMeasurements; identifier++) {
-                    // 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 (identifier != 0) {
-                            dp = 5*sp + dp*10;
-                        }
-                        if (lResolution->second[0].first == eLOC_0) {
-                            // push back & move a LOC_0 measurement
-                            Measurement<Identifier, eLOC_0> m0(
-                                    surface->getSharedPtr(), identifier, cov1D, lPos[eLOC_0] + dp);
-                            result[surface].push_back(m0);
-                        } else {
-                            // push back & move a LOC_1 measurement
-                            Measurement<Identifier, eLOC_1> m1(
-                                    surface->getSharedPtr(), identifier, cov1D, lPos[eLOC_1] + dp);
-                            result[surface].push_back(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);
-                        if (identifier != 0) {
-                            dx = 5*sx * 10*dx;
-                            dy = 5*sy * 10*dy;
-                        }
-                        // push back & move a LOC_0, LOC_1 measurement
-                        Measurement<Identifier, eLOC_0, eLOC_1> m01(
-                                surface->getSharedPtr(),
-                                identifier,
-                                cov2D,
-                                lPos[eLOC_0] + dx,
-                                lPos[eLOC_1] + dy);
-                        result[surface].push_back(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.));
-                }
-            }
-        };
-
-        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*/) {
-                std::cout << "I am on surface " << surface.geoID().toString() << std::endl;
-                const auto& measurements = m_measurementMap[&surface];
-                if(measurements.empty()) {
-                    return {};
-                }
-                std::cout << "I have found a measurement with ";
-                boost::apply_visitor([](const auto &concreteMeasurement) {
-                    std::cout << concreteMeasurement.parameters().transpose() << ", ";
-                }, measurements[0]);
-                std::cout << std::endl;
-
-                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(kalman_fitter_zero_field) {
-            // Build detector
-            CubicTrackingGeometry cGeometry;
-            auto detector = cGeometry();
-
-            // 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;
-
-            // 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);
-
-            //OK, let's make this a vector of start positions instead
-            std::vector<SingleCurvilinearTrackParameters<NeutralPolicy>> ParamsVector;
-            int nTracks = 6;
-            for (int i = 0; i < nTracks; i++) {
-                double xsize = 1;//get dimensions of first surface...
-                double ysize = 1;
-                Vector3D pos1(-3, uniform_dist(generator) * xsize, uniform_dist(generator) * ysize), mom1(
-                        1. * units::_GeV, 0., 0);
-                SingleCurvilinearTrackParameters<NeutralPolicy> start1(nullptr, pos1, mom1);
-                ParamsVector.push_back(start1);
-            }
-
-
-            // 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 meas_map = mResult.template get<MeasurementCreator::result_type>();
-
-            std::cout << "Number of Measurements: " << meas_map.size() << std::endl;
-
-            for(const auto& keyValue : meas_map) {
-                const Surface* surface = keyValue.first;
-                const auto& measurements = keyValue.second;
-                std::cout << "New measurement on surface " << surface->geoID().toString() << " and " << measurements.size() << " in total on this surface" << std::endl;
-                for(const auto& measurement : measurements) {
-                    boost::apply_visitor([](const auto &concreteMeasurement) {
-                        std::cout << concreteMeasurement.parameters().transpose() << ", ";
-                    }, measurement);
-                    std::cout << std::endl;
-                }
-            }
-
-            //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);
-            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);
-
-            Vector3D rPos(-3. * units::_m,
-                          10. * units::_um * gauss(generator),
-                          100. * units::_um * gauss(generator));
-            Vector3D rMom(1. * units::_GeV,
-                          0.025 * units::_GeV * gauss(generator),
-                          0.025 * units::_GeV * gauss(generator));
-
-            SingleCurvilinearTrackParameters<ChargedPolicy> rStart(
-                    std::move(covPtr), mPos, mMom, 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(meas_map)));
-
-            BOOST_TEST(testFit.fittedStates.size() == 6);
-            BOOST_TEST(testFit.numberOfHoles == 0);
-
-            std::cout << "I have found " << testFit.fittedStates.size() << " hits and have skipped " << testFit.numberOfHoles <<  std::endl;
-            for(const auto& state : testFit.fittedStates) {
-                std::cout << "\ton surface " << state.referenceSurface().geoID().toString() << ": ";
-                std::cout << std::endl << *state.parameter.predicted << std::endl << *state.parameter.filtered << ", " << std::endl;
-                boost::apply_visitor([](const auto &concreteMeasurement) {
-                    BOOST_TEST(concreteMeasurement.identifier() == 0);
-                    std::cout << concreteMeasurement.parameters().transpose() << ", ";
-                }, *state.measurement.uncalibrated);
-                std::cout << std::endl;
-            }
-
-            /*std::cout<<"do we get anything back? "<<testFit.fittedParameters.get()<<std::endl;
-
-            // Fit the track
-            auto fittedTrack      = kFitter.fit(measurements, rStart, rSurface);
-            auto fittedParameters = fittedTrack.fittedParameters.get();
-
-            // Make sure it is deterministic
-            auto fittedAgainTrack      = kFitter.fit(measurements, rStart, rSurface);
-            auto fittedAgainParameters = fittedAgainTrack.fittedParameters.get();
-
-            BOOST_TEST(fittedParameters.parameters().isApprox(
-                    fittedAgainParameters.parameters()));
-
-            // Change the order of the measurements
-            std::vector<TrackState> shuffledMeasurements = {measurements[3],
-                                                            measurements[2],
-                                                            measurements[1],
-                                                            measurements[4],
-                                                            measurements[5],
-                                                            measurements[0]};
-
-            // Make sure it works for shuffled measurements as well
-            auto fittedShuffledTrack
-                    = kFitter.fit(shuffledMeasurements, rStart, rSurface);
-            auto fittedShuffledParameters = fittedShuffledTrack.fittedParameters.get();
-
-            BOOST_TEST(fittedParameters.parameters().isApprox(
-                    fittedShuffledParameters.parameters()));
-
-            // Remove one measurement and find a hole
-            std::vector<TrackState> measurementsWithHole = {measurements[0],
-                                                            measurements[1],
-                                                            measurements[2],
-                                                            measurements[4],
-                                                            measurements[5]};
-
-            // Make sure it works for shuffled measurements as well
-            auto fittedWithHoleTrack
-                    = kFitter.fit(measurementsWithHole, rStart, rSurface);
-            auto fittedWithHoleParameters = fittedWithHoleTrack.fittedParameters.get();
-
-            // Count one hole
-            BOOST_TEST(fittedWithHoleTrack.missedActiveSurfaces.size() == 1);
-            // And the parameters should be different
-            BOOST_TEST(!fittedParameters.parameters().isApprox(
-                    fittedWithHoleParameters.parameters()));*/
+namespace Test {
+
+  // A few initialisations and definitionas
+  using Identifier = unsigned int;
+  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::mt19937                           generator(42);
+  std::uniform_real_distribution<double> uniform_dist(0., 1.);
+  std::uniform_int_distribution<int>     numberOfMeasurements(1, 5);
+
+  ActsSymMatrixD<1> cov1D;
+  ActsSymMatrixD<2> cov2D;
+
+  bool debugMode = true;
+
+  /// @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::map<const Surface*,
+                                 std::vector<FittableMeasurement<Identifier>>>;
+
+    /// @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 == nullptr or not surface->associatedDetectorElement()) {
+        return;
+      }
+      // random number to determine how many extra measurements on this surface
+      unsigned int nMeasurements
+          = static_cast<unsigned int>(numberOfMeasurements(generator));
+
+      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;
+      }
+
+      for (unsigned int identifier = 0; identifier < nMeasurements;
+           identifier++) {
+        // 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 (identifier != 0) {
+            dp = 5 * sp + dp * 10;
+          }
+          if (lResolution->second[0].first == eLOC_0) {
+            // push back & move a LOC_0 measurement
+            Measurement<Identifier, eLOC_0> m0(
+                surface->getSharedPtr(), identifier, cov1D, lPos[eLOC_0] + dp);
+            result[surface].push_back(m0);
+          } else {
+            // push back & move a LOC_1 measurement
+            Measurement<Identifier, eLOC_1> m1(
+                surface->getSharedPtr(), identifier, cov1D, lPos[eLOC_1] + dp);
+            result[surface].push_back(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);
+          if (identifier != 0) {
+            dx = 5 * sx * 10 * dx;
+            dy = 5 * sy * 10 * dy;
+          }
+          // push back & move a LOC_0, LOC_1 measurement
+          Measurement<Identifier, eLOC_0, eLOC_1> m01(surface->getSharedPtr(),
+                                                      identifier,
+                                                      cov2D,
+                                                      lPos[eLOC_0] + dx,
+                                                      lPos[eLOC_1] + dy);
+          result[surface].push_back(m01);
         }
-
-    }  // namespace Test
+      }
+    }
+  };
+
+  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.));
+      }
+    }
+  };
+
+  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*/)
+    {
+      std::cout << "I am on surface " << surface.geoID().toString()
+                << std::endl;
+      const auto& measurements = m_measurementMap[&surface];
+      if (measurements.empty()) {
+        return {};
+      }
+      std::cout << "I have found a measurement with ";
+      boost::apply_visitor(
+          [](const auto& concreteMeasurement) {
+            std::cout << concreteMeasurement.parameters().transpose() << ", ";
+          },
+          measurements[0]);
+      std::cout << std::endl;
+
+      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(kalman_fitter_zero_field)
+  {
+    // Build detector
+    CubicTrackingGeometry cGeometry;
+    auto                  detector = cGeometry();
+
+    // 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;
+
+    // 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);
+
+    // OK, let's make this a vector of start positions instead
+    std::vector<SingleCurvilinearTrackParameters<NeutralPolicy>> ParamsVector;
+    int                                                          nTracks = 6;
+    for (int i = 0; i < nTracks; i++) {
+      double   xsize = 1;  // get dimensions of first surface...
+      double   ysize = 1;
+      Vector3D pos1(
+          -3, uniform_dist(generator) * xsize, uniform_dist(generator) * ysize),
+          mom1(1. * units::_GeV, 0., 0);
+      SingleCurvilinearTrackParameters<NeutralPolicy> start1(
+          nullptr, pos1, mom1);
+      ParamsVector.push_back(start1);
+    }
+
+    // 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 meas_map = mResult.template get<MeasurementCreator::result_type>();
+
+    std::cout << "Number of Measurements: " << meas_map.size() << std::endl;
+
+    for (const auto& keyValue : meas_map) {
+      const Surface* surface      = keyValue.first;
+      const auto&    measurements = keyValue.second;
+      std::cout << "New measurement on surface " << surface->geoID().toString()
+                << " and " << measurements.size() << " in total on this surface"
+                << std::endl;
+      for (const auto& measurement : measurements) {
+        boost::apply_visitor(
+            [](const auto& concreteMeasurement) {
+              std::cout << concreteMeasurement.parameters().transpose() << ", ";
+            },
+            measurement);
+        std::cout << std::endl;
+      }
+    }
+
+    // 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);
+    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);
+
+    Vector3D rPos(-3. * units::_m,
+                  10. * units::_um * gauss(generator),
+                  100. * units::_um * gauss(generator));
+    Vector3D rMom(1. * units::_GeV,
+                  0.025 * units::_GeV * gauss(generator),
+                  0.025 * units::_GeV * gauss(generator));
+
+    SingleCurvilinearTrackParameters<ChargedPolicy> rStart(
+        std::move(covPtr), mPos, mMom, 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(meas_map)));
+
+    BOOST_TEST(testFit.fittedStates.size() == 6);
+    BOOST_TEST(testFit.numberOfHoles == 0);
+
+    std::cout << "I have found " << testFit.fittedStates.size()
+              << " hits and have skipped " << testFit.numberOfHoles
+              << std::endl;
+    for (const auto& state : testFit.fittedStates) {
+      std::cout << "\ton surface "
+                << state.referenceSurface().geoID().toString() << ": ";
+      std::cout << std::endl
+                << *state.parameter.predicted << std::endl
+                << *state.parameter.filtered << ", " << std::endl;
+      boost::apply_visitor(
+          [](const auto& concreteMeasurement) {
+            BOOST_TEST(concreteMeasurement.identifier() == 0);
+            std::cout << concreteMeasurement.parameters().transpose() << ", ";
+          },
+          *state.measurement.uncalibrated);
+      std::cout << std::endl;
+    }
+
+    /*std::cout<<"do we get anything back?
+    "<<testFit.fittedParameters.get()<<std::endl;
+
+    // Fit the track
+    auto fittedTrack      = kFitter.fit(measurements, rStart, rSurface);
+    auto fittedParameters = fittedTrack.fittedParameters.get();
+
+    // Make sure it is deterministic
+    auto fittedAgainTrack      = kFitter.fit(measurements, rStart, rSurface);
+    auto fittedAgainParameters = fittedAgainTrack.fittedParameters.get();
+
+    BOOST_TEST(fittedParameters.parameters().isApprox(
+            fittedAgainParameters.parameters()));
+
+    // Change the order of the measurements
+    std::vector<TrackState> shuffledMeasurements = {measurements[3],
+                                                    measurements[2],
+                                                    measurements[1],
+                                                    measurements[4],
+                                                    measurements[5],
+                                                    measurements[0]};
+
+    // Make sure it works for shuffled measurements as well
+    auto fittedShuffledTrack
+            = kFitter.fit(shuffledMeasurements, rStart, rSurface);
+    auto fittedShuffledParameters = fittedShuffledTrack.fittedParameters.get();
+
+    BOOST_TEST(fittedParameters.parameters().isApprox(
+            fittedShuffledParameters.parameters()));
+
+    // Remove one measurement and find a hole
+    std::vector<TrackState> measurementsWithHole = {measurements[0],
+                                                    measurements[1],
+                                                    measurements[2],
+                                                    measurements[4],
+                                                    measurements[5]};
+
+    // Make sure it works for shuffled measurements as well
+    auto fittedWithHoleTrack
+            = kFitter.fit(measurementsWithHole, rStart, rSurface);
+    auto fittedWithHoleParameters = fittedWithHoleTrack.fittedParameters.get();
+
+    // Count one hole
+    BOOST_TEST(fittedWithHoleTrack.missedActiveSurfaces.size() == 1);
+    // And the parameters should be different
+    BOOST_TEST(!fittedParameters.parameters().isApprox(
+            fittedWithHoleParameters.parameters()));*/
+  }
+
+}  // namespace Test
 }  // namespace Acts
-- 
GitLab


From c067e73d355a811f431363197aab861670838fc6 Mon Sep 17 00:00:00 2001
From: Nils Braun <nils.braun@kit.edu>
Date: Thu, 17 Jan 2019 17:57:01 -0800
Subject: [PATCH 08/17] Some renaming and adding of comments

---
 .../Acts/Fitter/SequentialKalmanFitter.hpp    | 116 +++++++++++++-----
 .../Fitter/SequentialKalmanFitterTests.cpp    |  70 ++---------
 2 files changed, 91 insertions(+), 95 deletions(-)

diff --git a/Core/include/Acts/Fitter/SequentialKalmanFitter.hpp b/Core/include/Acts/Fitter/SequentialKalmanFitter.hpp
index 670b0f06e..c2be6bc39 100644
--- a/Core/include/Acts/Fitter/SequentialKalmanFitter.hpp
+++ b/Core/include/Acts/Fitter/SequentialKalmanFitter.hpp
@@ -8,20 +8,34 @@
 
 #pragma once
 
-#include <boost/variant.hpp>
-#include <memory>
-#include "Acts/EventData/Measurement.hpp"
-#include "Acts/EventData/TrackParameters.hpp"
 #include "Acts/Fitter/detail/VoidKalmanComponents.hpp"
 #include "Acts/Propagator/AbortList.hpp"
 #include "Acts/Propagator/ActionList.hpp"
 #include "Acts/Propagator/Propagator.hpp"
-#include "Acts/Propagator/detail/ConstrainedStep.hpp"
-#include "Acts/Propagator/detail/StandardAborters.hpp"
-#include "Acts/Utilities/Definitions.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,
@@ -33,16 +47,26 @@ public:
   /// Default constructor is deleted
   SequentialKalmanFitter() = delete;
 
-  /// Constructor from arguments
-  SequentialKalmanFitter(propagator_t pPropagator)
+  /// 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> pHitSelector) const
+      std::shared_ptr<measurement_selector_t> pMeasurementSelector) const
   {
     // Create the ActionList and AbortList
     using KalmanResult = typename KalmanActor::result_type;
@@ -52,13 +76,13 @@ public:
     // Create relevant options for the propagation options
     PropagatorOptions<Actors, Aborters> kalmanOptions;
     auto& kalmanActor = kalmanOptions.actionList.template get<KalmanActor>();
-    kalmanActor.m_measurementSelectorPtr = std::move(pHitSelector);
+    kalmanActor.m_measurementSelectorPtr = std::move(pMeasurementSelector);
 
     // Run the fitter
     const auto& result
         = m_propagator.template propagate(sParameters, kalmanOptions);
 
-    /// Get the result of the fit
+    // Get the result of the fit
     auto kalmanResult = result.template get<KalmanResult>();
 
     // Return the converted Track
@@ -66,45 +90,71 @@ public:
   }
 
 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))
     {
     }
 
-    struct this_result
+    /// The type of the result of this actor including the final track and the
+    /// number of holes.
+    struct result_type
     {
-      // Move the result into the fitted states
-      std::vector<track_state_t> fittedStates  = {};
-      unsigned int               numberOfHoles = 0;
+      /// The resulting track of the algorithm
+      std::vector<track_state_t> resultStateList = {};
+      /// How many surfaces had no measurement on them
+      unsigned int numberOfHoles = 0;
     };
 
-    using result_type = this_result;
-
+    /**
+     * 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;
-      }
+      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;
       }
 
-      std::cout << "Looking for the best among.." << std::endl;
+      // 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) {
@@ -113,28 +163,30 @@ private:
         nextState.parameter.pathLength = std::get<2>(boundStates);
 
         m_updator(nextState);
-        std::cout << "chi2 = " << nextState.filteredChi2Increment << " ";
 
         if (chosenNextState == nullptr) {
           chosenNextState = &nextState;
-          std::cout << "first one" << std::endl;
         } else {
-          if (chosenNextState->filteredChi2Increment
-              > nextState.filteredChi2Increment) {
-            chosenNextState = &nextState;
-            std::cout << "is better" << std::endl;
-          }
+          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.fittedStates.push_back(std::move(*chosenNextState));
+      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;
-    updator_t                               m_updator;
-    calibrator_t                            m_calibrator;
+    /// The updator instance
+    updator_t m_updator;
+    /// The calibrator instance
+    calibrator_t m_calibrator;
   };
 };
 
diff --git a/Tests/Core/Fitter/SequentialKalmanFitterTests.cpp b/Tests/Core/Fitter/SequentialKalmanFitterTests.cpp
index 0e95daa34..1d61cbbd1 100644
--- a/Tests/Core/Fitter/SequentialKalmanFitterTests.cpp
+++ b/Tests/Core/Fitter/SequentialKalmanFitterTests.cpp
@@ -99,14 +99,10 @@ namespace Test {
       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;
-      }
+      if (vResolution == detectorResolution.end()) { return; }
       // find layer resolutions
       auto lResolution = vResolution->second.find(layerID);
-      if (lResolution == vResolution->second.end()) {
-        return;
-      }
+      if (lResolution == vResolution->second.end()) { return; }
 
       for (unsigned int identifier = 0; identifier < nMeasurements;
            identifier++) {
@@ -118,9 +114,7 @@ namespace Test {
           double sp = lResolution->second[0].second;
           cov1D << sp * sp;
           double dp = sp * gauss(generator);
-          if (identifier != 0) {
-            dp = 5 * sp + dp * 10;
-          }
+          if (identifier != 0) { dp = 5 * sp + dp * 10; }
           if (lResolution->second[0].first == eLOC_0) {
             // push back & move a LOC_0 measurement
             Measurement<Identifier, eLOC_0> m0(
@@ -223,9 +217,7 @@ namespace Test {
       std::cout << "I am on surface " << surface.geoID().toString()
                 << std::endl;
       const auto& measurements = m_measurementMap[&surface];
-      if (measurements.empty()) {
-        return {};
-      }
+      if (measurements.empty()) { return {}; }
       std::cout << "I have found a measurement with ";
       boost::apply_visitor(
           [](const auto& concreteMeasurement) {
@@ -392,13 +384,13 @@ namespace Test {
     auto         testFit = fitter.fit(
         rStart, std::make_shared<HitSelector>(std::move(meas_map)));
 
-    BOOST_TEST(testFit.fittedStates.size() == 6);
+    BOOST_TEST(testFit.resultStateList.size() == 6);
     BOOST_TEST(testFit.numberOfHoles == 0);
 
-    std::cout << "I have found " << testFit.fittedStates.size()
+    std::cout << "I have found " << testFit.resultStateList.size()
               << " hits and have skipped " << testFit.numberOfHoles
               << std::endl;
-    for (const auto& state : testFit.fittedStates) {
+    for (const auto& state : testFit.resultStateList) {
       std::cout << "\ton surface "
                 << state.referenceSurface().geoID().toString() << ": ";
       std::cout << std::endl
@@ -412,54 +404,6 @@ namespace Test {
           *state.measurement.uncalibrated);
       std::cout << std::endl;
     }
-
-    /*std::cout<<"do we get anything back?
-    "<<testFit.fittedParameters.get()<<std::endl;
-
-    // Fit the track
-    auto fittedTrack      = kFitter.fit(measurements, rStart, rSurface);
-    auto fittedParameters = fittedTrack.fittedParameters.get();
-
-    // Make sure it is deterministic
-    auto fittedAgainTrack      = kFitter.fit(measurements, rStart, rSurface);
-    auto fittedAgainParameters = fittedAgainTrack.fittedParameters.get();
-
-    BOOST_TEST(fittedParameters.parameters().isApprox(
-            fittedAgainParameters.parameters()));
-
-    // Change the order of the measurements
-    std::vector<TrackState> shuffledMeasurements = {measurements[3],
-                                                    measurements[2],
-                                                    measurements[1],
-                                                    measurements[4],
-                                                    measurements[5],
-                                                    measurements[0]};
-
-    // Make sure it works for shuffled measurements as well
-    auto fittedShuffledTrack
-            = kFitter.fit(shuffledMeasurements, rStart, rSurface);
-    auto fittedShuffledParameters = fittedShuffledTrack.fittedParameters.get();
-
-    BOOST_TEST(fittedParameters.parameters().isApprox(
-            fittedShuffledParameters.parameters()));
-
-    // Remove one measurement and find a hole
-    std::vector<TrackState> measurementsWithHole = {measurements[0],
-                                                    measurements[1],
-                                                    measurements[2],
-                                                    measurements[4],
-                                                    measurements[5]};
-
-    // Make sure it works for shuffled measurements as well
-    auto fittedWithHoleTrack
-            = kFitter.fit(measurementsWithHole, rStart, rSurface);
-    auto fittedWithHoleParameters = fittedWithHoleTrack.fittedParameters.get();
-
-    // Count one hole
-    BOOST_TEST(fittedWithHoleTrack.missedActiveSurfaces.size() == 1);
-    // And the parameters should be different
-    BOOST_TEST(!fittedParameters.parameters().isApprox(
-            fittedWithHoleParameters.parameters()));*/
   }
 
 }  // namespace Test
-- 
GitLab


From 4d284b8d957cb6cae804cb287d11f422f1fac0ff Mon Sep 17 00:00:00 2001
From: Nils Braun <nils.braun@kit.edu>
Date: Thu, 17 Jan 2019 18:02:26 -0800
Subject: [PATCH 09/17] Clang format again

---
 .../Acts/Fitter/SequentialKalmanFitter.hpp       |  8 ++++++--
 .../Core/Fitter/SequentialKalmanFitterTests.cpp  | 16 ++++++++++++----
 2 files changed, 18 insertions(+), 6 deletions(-)

diff --git a/Core/include/Acts/Fitter/SequentialKalmanFitter.hpp b/Core/include/Acts/Fitter/SequentialKalmanFitter.hpp
index c2be6bc39..700da210f 100644
--- a/Core/include/Acts/Fitter/SequentialKalmanFitter.hpp
+++ b/Core/include/Acts/Fitter/SequentialKalmanFitter.hpp
@@ -139,7 +139,9 @@ private:
     {
       // Do not go on if there is no surface
       const Surface* surface = state.navigation.currentSurface;
-      if (not surface) { return; }
+      if (not surface) {
+        return;
+      }
 
       // Request a list of possible next states by the measurement selector
       assert(m_measurementSelectorPtr != nullptr);
@@ -169,7 +171,9 @@ private:
         } else {
           const double chosenChi2 = chosenNextState->filteredChi2Increment;
           const double nextChi2   = nextState.filteredChi2Increment;
-          if (chosenChi2 > nextChi2) { chosenNextState = &nextState; }
+          if (chosenChi2 > nextChi2) {
+            chosenNextState = &nextState;
+          }
         }
       }
 
diff --git a/Tests/Core/Fitter/SequentialKalmanFitterTests.cpp b/Tests/Core/Fitter/SequentialKalmanFitterTests.cpp
index 1d61cbbd1..553ec18e8 100644
--- a/Tests/Core/Fitter/SequentialKalmanFitterTests.cpp
+++ b/Tests/Core/Fitter/SequentialKalmanFitterTests.cpp
@@ -99,10 +99,14 @@ namespace Test {
       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; }
+      if (vResolution == detectorResolution.end()) {
+        return;
+      }
       // find layer resolutions
       auto lResolution = vResolution->second.find(layerID);
-      if (lResolution == vResolution->second.end()) { return; }
+      if (lResolution == vResolution->second.end()) {
+        return;
+      }
 
       for (unsigned int identifier = 0; identifier < nMeasurements;
            identifier++) {
@@ -114,7 +118,9 @@ namespace Test {
           double sp = lResolution->second[0].second;
           cov1D << sp * sp;
           double dp = sp * gauss(generator);
-          if (identifier != 0) { dp = 5 * sp + dp * 10; }
+          if (identifier != 0) {
+            dp = 5 * sp + dp * 10;
+          }
           if (lResolution->second[0].first == eLOC_0) {
             // push back & move a LOC_0 measurement
             Measurement<Identifier, eLOC_0> m0(
@@ -217,7 +223,9 @@ namespace Test {
       std::cout << "I am on surface " << surface.geoID().toString()
                 << std::endl;
       const auto& measurements = m_measurementMap[&surface];
-      if (measurements.empty()) { return {}; }
+      if (measurements.empty()) {
+        return {};
+      }
       std::cout << "I have found a measurement with ";
       boost::apply_visitor(
           [](const auto& concreteMeasurement) {
-- 
GitLab


From 5f097432ff58d6fb2941b63014fd29af17da616e Mon Sep 17 00:00:00 2001
From: Nils Braun <nils.braun@kit.edu>
Date: Wed, 23 Jan 2019 15:08:56 -0800
Subject: [PATCH 10/17] Fixed a typo

---
 Tests/Core/Fitter/SequentialKalmanFitterTests.cpp | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/Tests/Core/Fitter/SequentialKalmanFitterTests.cpp b/Tests/Core/Fitter/SequentialKalmanFitterTests.cpp
index 553ec18e8..9e2c34639 100644
--- a/Tests/Core/Fitter/SequentialKalmanFitterTests.cpp
+++ b/Tests/Core/Fitter/SequentialKalmanFitterTests.cpp
@@ -140,8 +140,8 @@ namespace Test {
           double dx = sx * gauss(generator);
           double dy = sy * gauss(generator);
           if (identifier != 0) {
-            dx = 5 * sx * 10 * dx;
-            dy = 5 * sy * 10 * dy;
+            dx = 5 * sx + 10 * dx;
+            dy = 5 * sy + 10 * dy;
           }
           // push back & move a LOC_0, LOC_1 measurement
           Measurement<Identifier, eLOC_0, eLOC_1> m01(surface->getSharedPtr(),
-- 
GitLab


From dcdedff37637abfabccba4e4613ae9398e19ea29 Mon Sep 17 00:00:00 2001
From: Nils Braun <nils.braun@kit.edu>
Date: Wed, 23 Jan 2019 17:01:43 -0800
Subject: [PATCH 11/17] Added a missing pragma once

---
 .../include/Acts/Tests/CommonHelpers/CubicTrackingGeometry.hpp  | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Tests/Core/CommonHelpers/include/Acts/Tests/CommonHelpers/CubicTrackingGeometry.hpp b/Tests/Core/CommonHelpers/include/Acts/Tests/CommonHelpers/CubicTrackingGeometry.hpp
index 594fe8d5d..195f8df58 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"
-- 
GitLab


From c110187480e33382b85eb8ccfeaec1c69d97d240 Mon Sep 17 00:00:00 2001
From: Nils Braun <nils.braun@kit.edu>
Date: Wed, 23 Jan 2019 17:02:02 -0800
Subject: [PATCH 12/17] Started developing a measurement creator helper for
 tests

---
 .../MeasurementCreatorHelper.hpp              | 256 ++++++++++++++++++
 1 file changed, 256 insertions(+)
 create mode 100644 Tests/Core/CommonHelpers/include/Acts/Tests/CommonHelpers/MeasurementCreatorHelper.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 000000000..24aab20cb
--- /dev/null
+++ b/Tests/Core/CommonHelpers/include/Acts/Tests/CommonHelpers/MeasurementCreatorHelper.hpp
@@ -0,0 +1,256 @@
+// 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 = GeometryID;
+
+  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);
+
+  /// @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 == 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; }
+
+      // 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 (lResolution->second[0].first == eLOC_0) {
+          // push back & move a LOC_0 measurement
+          Measurement<Identifier, eLOC_0> m(
+              surface->getSharedPtr(), geoID, cov1D, lPos[eLOC_0] + dp);
+          result.push_back(TrackState(std::move(m)));
+        } else {
+          // push back & move a LOC_1 measurement
+          Measurement<Identifier, eLOC_1> m(
+              surface->getSharedPtr(), geoID, cov1D, lPos[eLOC_1] + dp);
+          result.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);
+        // push back & move a LOC_0, LOC_1 measurement
+        Measurement<Identifier, eLOC_0, eLOC_1> m(surface->getSharedPtr(),
+                                                  geoID,
+                                                  cov2D,
+                                                  lPos[eLOC_0] + dx,
+                                                  lPos[eLOC_1] + dy);
+        result.push_back(TrackState(std::move(m)));
+      }
+    }
+  };
+
+  ///
+  /// @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.));
+      }
+    }
+  };
+
+  template <class ATrackingGeometry, class AStepper>
+  MeasurementCreator::result_type
+  createMeasurements(AStepper           stepper,
+                     const Vector3D&    startPosition,
+                     const Vector3D&    startMomentum,
+                     DetectorResolution detRes,
+                     bool               debugMode = false)
+  {
+    // Build detector
+    ATrackingGeometry cGeometry;
+    auto              detector = cGeometry();
+
+    // 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;
+
+    // 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(AStepper           stepper,
+                                    const Vector3D&    startPosition,
+                                    const Vector3D&    startMomentum,
+                                    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<CubicTrackingGeometry>(stepper, startPosition, startMomentum, detRes, debugMode);
+  }
+
+}  // namespace Test
+}  // namespace Acts
-- 
GitLab


From 44d26c1c2d15dee599ff91236aa55017d5bfa458 Mon Sep 17 00:00:00 2001
From: Nils Braun <nils.braun@kit.edu>
Date: Wed, 23 Jan 2019 17:52:55 -0800
Subject: [PATCH 13/17] Also make it possible to generate noise hits

---
 .../MeasurementCreatorHelper.hpp              | 147 +++++++++---------
 1 file changed, 73 insertions(+), 74 deletions(-)

diff --git a/Tests/Core/CommonHelpers/include/Acts/Tests/CommonHelpers/MeasurementCreatorHelper.hpp b/Tests/Core/CommonHelpers/include/Acts/Tests/CommonHelpers/MeasurementCreatorHelper.hpp
index 24aab20cb..25d7bc98a 100644
--- a/Tests/Core/CommonHelpers/include/Acts/Tests/CommonHelpers/MeasurementCreatorHelper.hpp
+++ b/Tests/Core/CommonHelpers/include/Acts/Tests/CommonHelpers/MeasurementCreatorHelper.hpp
@@ -30,7 +30,7 @@ namespace Acts {
 namespace Test {
 
   // A few initialisations and definitionas
-  using Identifier = GeometryID;
+  using Identifier = unsigned int;
 
   using TrackState         = TrackState<Identifier, BoundParameters>;
   using Resolution         = std::pair<ParID_t, double>;
@@ -40,8 +40,9 @@ namespace Test {
 
   using DebugOutput = detail::DebugOutputActor;
 
-  std::normal_distribution<double> gauss(0., 1.);
-  std::default_random_engine       generator(42);
+  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
@@ -54,7 +55,15 @@ namespace Test {
     /// The detector resolution
     DetectorResolution detectorResolution;
 
-    using result_type = std::vector<TrackState>;
+    /// 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
@@ -67,7 +76,7 @@ namespace Test {
     operator()(propagator_state_t& state, result_type& result) const
     {
       // monitor the current surface
-      auto surface = state.navigation.currentSurface;
+      const Surface* surface = state.navigation.currentSurface;
       if (surface == nullptr or not surface->associatedDetectorElement()) {
         return;
       }
@@ -82,6 +91,29 @@ namespace Test {
       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
+    {
+      auto geoID = surface->geoID();
+
       // Apply global to local
       Acts::Vector2D lPos;
       surface->globalToLocal(
@@ -91,16 +123,21 @@ namespace Test {
         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(), geoID, cov1D, lPos[eLOC_0] + dp);
-          result.push_back(TrackState(std::move(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(), geoID, cov1D, lPos[eLOC_1] + dp);
-          result.push_back(TrackState(std::move(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
@@ -110,77 +147,34 @@ namespace Test {
         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(),
-                                                  geoID,
+                                                  identifier,
                                                   cov2D,
                                                   lPos[eLOC_0] + dx,
                                                   lPos[eLOC_1] + dy);
-        result.push_back(TrackState(std::move(m)));
+        result.measurementMap[surface].push_back(m);
+        result.trackStates.push_back(TrackState(std::move(m)));
       }
     }
   };
 
-  ///
-  /// @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.));
-      }
-    }
-  };
-
-  template <class ATrackingGeometry, class AStepper>
+  template <class AStepper>
   MeasurementCreator::result_type
-  createMeasurements(AStepper           stepper,
-                     const Vector3D&    startPosition,
-                     const Vector3D&    startMomentum,
-                     DetectorResolution detRes,
-                     bool               debugMode = false)
+  createMeasurements(std::shared_ptr<const TrackingGeometry> detector,
+                     AStepper                                stepper,
+                     const Vector3D&                         startPosition,
+                     const Vector3D&                         startMomentum,
+                     DetectorResolution                      detRes,
+                     bool addNoiseMeasurements = false,
+                     bool debugMode            = false)
   {
-    // Build detector
-    ATrackingGeometry cGeometry;
-    auto              detector = cGeometry();
-
     // Build navigator for the measurement creation
     Navigator navigator(detector);
     navigator.resolvePassive   = false;
@@ -205,6 +199,7 @@ namespace Test {
     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);
@@ -221,10 +216,13 @@ namespace Test {
 
   template <class AStepper>
   MeasurementCreator::result_type
-  createMeasurementsOnCubicDetector(AStepper           stepper,
-                                    const Vector3D&    startPosition,
-                                    const Vector3D&    startMomentum,
-                                    bool               debugMode = false)
+  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);
@@ -249,7 +247,8 @@ namespace Test {
     detRes[2] = pixelVolumeRes;
     detRes[3] = stripVolumeRes;
 
-    return createMeasurements<CubicTrackingGeometry>(stepper, startPosition, startMomentum, detRes, debugMode);
+    return createMeasurements(
+        detector, stepper, startPosition, startMomentum, detRes, addNoiseMeasurements, debugMode);
   }
 
 }  // namespace Test
-- 
GitLab


From b6f543583bb2062cf0bdb72f43d9a35c728766ca Mon Sep 17 00:00:00 2001
From: Nils Braun <nils.braun@kit.edu>
Date: Thu, 24 Jan 2019 11:11:10 -0800
Subject: [PATCH 14/17] Use the new measurement creator helper for the Kalman
 and sequential kalman test

---
 Tests/Core/Fitter/KalmanFitterTests.cpp       | 226 +-----------
 .../Fitter/SequentialKalmanFitterTests.cpp    | 330 +-----------------
 2 files changed, 26 insertions(+), 530 deletions(-)

diff --git a/Tests/Core/Fitter/KalmanFitterTests.cpp b/Tests/Core/Fitter/KalmanFitterTests.cpp
index dd5655cfa..89fd7887f 100644
--- a/Tests/Core/Fitter/KalmanFitterTests.cpp
+++ b/Tests/Core/Fitter/KalmanFitterTests.cpp
@@ -30,250 +30,46 @@
 #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
index 9e2c34639..3ec34e822 100644
--- a/Tests/Core/Fitter/SequentialKalmanFitterTests.cpp
+++ b/Tests/Core/Fitter/SequentialKalmanFitterTests.cpp
@@ -6,7 +6,7 @@
 // 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 KalmanFitter Tests
+#define BOOST_TEST_MODULE SequentialKalmanFitter Tests
 
 #include <boost/test/included/unit_test.hpp>
 
@@ -19,194 +19,18 @@
 #include "Acts/EventData/TrackParameters.hpp"
 #include "Acts/EventData/TrackState.hpp"
 #include "Acts/Extrapolator/Navigator.hpp"
-#include "Acts/Extrapolator/SurfaceCollector.hpp"
 #include "Acts/Fitter/GainMatrixSmoother.hpp"
 #include "Acts/Fitter/GainMatrixUpdator.hpp"
-#include "Acts/Fitter/KalmanFitter.hpp"
 #include "Acts/Fitter/SequentialKalmanFitter.hpp"
 #include "Acts/MagneticField/ConstantBField.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/Tests/CommonHelpers/CubicTrackingGeometry.hpp"
-#include "Acts/Utilities/BinningType.hpp"
-#include "Acts/Utilities/Definitions.hpp"
-#include "Acts/Utilities/GeometryID.hpp"
+#include "Acts/Tests/CommonHelpers/MeasurementCreatorHelper.hpp"
 
 namespace Acts {
 namespace Test {
 
-  // A few initialisations and definitionas
-  using Identifier = unsigned int;
-  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::mt19937                           generator(42);
-  std::uniform_real_distribution<double> uniform_dist(0., 1.);
-  std::uniform_int_distribution<int>     numberOfMeasurements(1, 5);
-
-  ActsSymMatrixD<1> cov1D;
-  ActsSymMatrixD<2> cov2D;
-
-  bool debugMode = true;
-
-  /// @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::map<const Surface*,
-                                 std::vector<FittableMeasurement<Identifier>>>;
-
-    /// @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 == nullptr or not surface->associatedDetectorElement()) {
-        return;
-      }
-      // random number to determine how many extra measurements on this surface
-      unsigned int nMeasurements
-          = static_cast<unsigned int>(numberOfMeasurements(generator));
-
-      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;
-      }
-
-      for (unsigned int identifier = 0; identifier < nMeasurements;
-           identifier++) {
-        // 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 (identifier != 0) {
-            dp = 5 * sp + dp * 10;
-          }
-          if (lResolution->second[0].first == eLOC_0) {
-            // push back & move a LOC_0 measurement
-            Measurement<Identifier, eLOC_0> m0(
-                surface->getSharedPtr(), identifier, cov1D, lPos[eLOC_0] + dp);
-            result[surface].push_back(m0);
-          } else {
-            // push back & move a LOC_1 measurement
-            Measurement<Identifier, eLOC_1> m1(
-                surface->getSharedPtr(), identifier, cov1D, lPos[eLOC_1] + dp);
-            result[surface].push_back(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);
-          if (identifier != 0) {
-            dx = 5 * sx + 10 * dx;
-            dy = 5 * sy + 10 * dy;
-          }
-          // push back & move a LOC_0, LOC_1 measurement
-          Measurement<Identifier, eLOC_0, eLOC_1> m01(surface->getSharedPtr(),
-                                                      identifier,
-                                                      cov2D,
-                                                      lPos[eLOC_0] + dx,
-                                                      lPos[eLOC_1] + dy);
-          result[surface].push_back(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.));
-      }
-    }
-  };
-
   struct HitSelector
   {
   public:
@@ -220,19 +44,8 @@ namespace Test {
     std::vector<TrackState>
     operator()(const Surface& surface, const propagator_state_t& /*state*/)
     {
-      std::cout << "I am on surface " << surface.geoID().toString()
-                << std::endl;
       const auto& measurements = m_measurementMap[&surface];
-      if (measurements.empty()) {
-        return {};
-      }
-      std::cout << "I have found a measurement with ";
-      boost::apply_visitor(
-          [](const auto& concreteMeasurement) {
-            std::cout << concreteMeasurement.parameters().transpose() << ", ";
-          },
-          measurements[0]);
-      std::cout << std::endl;
+      if (measurements.empty()) { return {}; }
 
       std::vector<TrackState> trackStates;
       for (const auto& measurement : measurements) {
@@ -250,106 +63,19 @@ namespace Test {
   ///
   /// @brief Unit test for Kalman fitter with measurements along the x-axis
   ///
-  BOOST_AUTO_TEST_CASE(kalman_fitter_zero_field)
+  BOOST_AUTO_TEST_CASE(sequential_kalman_fitter_zero_field)
   {
-    // Build detector
-    CubicTrackingGeometry cGeometry;
-    auto                  detector = cGeometry();
-
-    // Build navigator for the measurement creatoin
-    Navigator mNavigator(detector);
-    mNavigator.resolvePassive   = false;
-    mNavigator.resolveMaterial  = true;
-    mNavigator.resolveSensitive = true;
+    CubicTrackingGeometry geom;
+    auto                  detector = geom();
 
-    // 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);
-
-    // OK, let's make this a vector of start positions instead
-    std::vector<SingleCurvilinearTrackParameters<NeutralPolicy>> ParamsVector;
-    int                                                          nTracks = 6;
-    for (int i = 0; i < nTracks; i++) {
-      double   xsize = 1;  // get dimensions of first surface...
-      double   ysize = 1;
-      Vector3D pos1(
-          -3, uniform_dist(generator) * xsize, uniform_dist(generator) * ysize),
-          mom1(1. * units::_GeV, 0., 0);
-      SingleCurvilinearTrackParameters<NeutralPolicy> start1(
-          nullptr, pos1, mom1);
-      ParamsVector.push_back(start1);
-    }
-
-    // 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 meas_map = mResult.template get<MeasurementCreator::result_type>();
-
-    std::cout << "Number of Measurements: " << meas_map.size() << std::endl;
-
-    for (const auto& keyValue : meas_map) {
-      const Surface* surface      = keyValue.first;
-      const auto&    measurements = keyValue.second;
-      std::cout << "New measurement on surface " << surface->geoID().toString()
-                << " and " << measurements.size() << " in total on this surface"
-                << std::endl;
-      for (const auto& measurement : measurements) {
-        boost::apply_visitor(
-            [](const auto& concreteMeasurement) {
-              std::cout << concreteMeasurement.parameters().transpose() << ", ";
-            },
-            measurement);
-        std::cout << std::endl;
-      }
-    }
-
-    // BOOST_TEST(measurements.size() == 6);
+    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
@@ -359,9 +85,8 @@ namespace Test {
     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);
 
@@ -372,15 +97,8 @@ namespace Test {
 
     auto covPtr = std::make_unique<const ActsSymMatrixD<5>>(cov);
 
-    Vector3D rPos(-3. * units::_m,
-                  10. * units::_um * gauss(generator),
-                  100. * units::_um * gauss(generator));
-    Vector3D rMom(1. * units::_GeV,
-                  0.025 * units::_GeV * gauss(generator),
-                  0.025 * units::_GeV * gauss(generator));
-
     SingleCurvilinearTrackParameters<ChargedPolicy> rStart(
-        std::move(covPtr), mPos, mMom, 1.);
+        std::move(covPtr), startPosition, startMomentum, 1.);
 
     using Updator      = GainMatrixUpdator<BoundParameters>;
     using KalmanFitter = SequentialKalmanFitter<RecoPropagator,
@@ -390,28 +108,10 @@ namespace Test {
 
     KalmanFitter fitter(std::move(rPropagator));
     auto         testFit = fitter.fit(
-        rStart, std::make_shared<HitSelector>(std::move(meas_map)));
+        rStart, std::make_shared<HitSelector>(std::move(measurementMap)));
 
     BOOST_TEST(testFit.resultStateList.size() == 6);
     BOOST_TEST(testFit.numberOfHoles == 0);
-
-    std::cout << "I have found " << testFit.resultStateList.size()
-              << " hits and have skipped " << testFit.numberOfHoles
-              << std::endl;
-    for (const auto& state : testFit.resultStateList) {
-      std::cout << "\ton surface "
-                << state.referenceSurface().geoID().toString() << ": ";
-      std::cout << std::endl
-                << *state.parameter.predicted << std::endl
-                << *state.parameter.filtered << ", " << std::endl;
-      boost::apply_visitor(
-          [](const auto& concreteMeasurement) {
-            BOOST_TEST(concreteMeasurement.identifier() == 0);
-            std::cout << concreteMeasurement.parameters().transpose() << ", ";
-          },
-          *state.measurement.uncalibrated);
-      std::cout << std::endl;
-    }
   }
 
 }  // namespace Test
-- 
GitLab


From 46819426139f743eeddadc42a491fb6cd5d636b9 Mon Sep 17 00:00:00 2001
From: Nils Braun <nils.braun@kit.edu>
Date: Thu, 24 Jan 2019 11:20:45 -0800
Subject: [PATCH 15/17] clang-format

---
 .../MeasurementCreatorHelper.hpp              | 21 ++++++++++++++-----
 Tests/Core/Fitter/KalmanFitterTests.cpp       |  3 ++-
 .../Fitter/SequentialKalmanFitterTests.cpp    |  4 +++-
 3 files changed, 21 insertions(+), 7 deletions(-)

diff --git a/Tests/Core/CommonHelpers/include/Acts/Tests/CommonHelpers/MeasurementCreatorHelper.hpp b/Tests/Core/CommonHelpers/include/Acts/Tests/CommonHelpers/MeasurementCreatorHelper.hpp
index 25d7bc98a..99aa285fb 100644
--- a/Tests/Core/CommonHelpers/include/Acts/Tests/CommonHelpers/MeasurementCreatorHelper.hpp
+++ b/Tests/Core/CommonHelpers/include/Acts/Tests/CommonHelpers/MeasurementCreatorHelper.hpp
@@ -86,10 +86,14 @@ namespace Test {
       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; }
+      if (vResolution == detectorResolution.end()) {
+        return;
+      }
       // find layer resolutions
       auto lResolution = vResolution->second.find(layerID);
-      if (lResolution == vResolution->second.end()) { return; }
+      if (lResolution == vResolution->second.end()) {
+        return;
+      }
 
       unsigned int nMeasurements = 1;
 
@@ -124,7 +128,9 @@ namespace Test {
         cov1D << sp * sp;
         double dp = sp * gauss(generator);
 
-        if (identifier != 0) { dp *= 20; }
+        if (identifier != 0) {
+          dp *= 20;
+        }
 
         if (lResolution->second[0].first == eLOC_0) {
           // push back & move a LOC_0 measurement
@@ -247,8 +253,13 @@ namespace Test {
     detRes[2] = pixelVolumeRes;
     detRes[3] = stripVolumeRes;
 
-    return createMeasurements(
-        detector, stepper, startPosition, startMomentum, detRes, addNoiseMeasurements, debugMode);
+    return createMeasurements(detector,
+                              stepper,
+                              startPosition,
+                              startMomentum,
+                              detRes,
+                              addNoiseMeasurements,
+                              debugMode);
   }
 
 }  // namespace Test
diff --git a/Tests/Core/Fitter/KalmanFitterTests.cpp b/Tests/Core/Fitter/KalmanFitterTests.cpp
index 89fd7887f..8dcc2094e 100644
--- a/Tests/Core/Fitter/KalmanFitterTests.cpp
+++ b/Tests/Core/Fitter/KalmanFitterTests.cpp
@@ -50,7 +50,8 @@ namespace Test {
     Vector3D            startMomentum(1. * units::_GeV, 0., 0);
 
     auto measurements = createMeasurementsOnCubicDetector(
-        detector, mStepper, startPosition, startMomentum).trackStates;
+                            detector, mStepper, startPosition, startMomentum)
+                            .trackStates;
     BOOST_TEST(measurements.size() == 6);
 
     // The KalmanFitter - we use the eigen stepper for covariance transport
diff --git a/Tests/Core/Fitter/SequentialKalmanFitterTests.cpp b/Tests/Core/Fitter/SequentialKalmanFitterTests.cpp
index 3ec34e822..0623d247b 100644
--- a/Tests/Core/Fitter/SequentialKalmanFitterTests.cpp
+++ b/Tests/Core/Fitter/SequentialKalmanFitterTests.cpp
@@ -45,7 +45,9 @@ namespace Test {
     operator()(const Surface& surface, const propagator_state_t& /*state*/)
     {
       const auto& measurements = m_measurementMap[&surface];
-      if (measurements.empty()) { return {}; }
+      if (measurements.empty()) {
+        return {};
+      }
 
       std::vector<TrackState> trackStates;
       for (const auto& measurement : measurements) {
-- 
GitLab


From f807af414dee210f7ba498ee68d17b6857235ca8 Mon Sep 17 00:00:00 2001
From: Nils Braun <nils.braun@kit.edu>
Date: Thu, 24 Jan 2019 11:23:52 -0800
Subject: [PATCH 16/17] Removed unused file

---
 Tests/Core/Fitter/KalmanFilterTestUtils.hpp | 152 --------------------
 1 file changed, 152 deletions(-)
 delete mode 100644 Tests/Core/Fitter/KalmanFilterTestUtils.hpp

diff --git a/Tests/Core/Fitter/KalmanFilterTestUtils.hpp b/Tests/Core/Fitter/KalmanFilterTestUtils.hpp
deleted file mode 100644
index a7ffd9bee..000000000
--- 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
-- 
GitLab


From 84794f89d524b91797b96fdd165ed5b2f50b5d08 Mon Sep 17 00:00:00 2001
From: Nils Braun <nils.braun@kit.edu>
Date: Thu, 24 Jan 2019 12:30:53 -0800
Subject: [PATCH 17/17] Removed unused line

---
 .../Acts/Tests/CommonHelpers/MeasurementCreatorHelper.hpp       | 2 --
 1 file changed, 2 deletions(-)

diff --git a/Tests/Core/CommonHelpers/include/Acts/Tests/CommonHelpers/MeasurementCreatorHelper.hpp b/Tests/Core/CommonHelpers/include/Acts/Tests/CommonHelpers/MeasurementCreatorHelper.hpp
index 99aa285fb..f3962125e 100644
--- a/Tests/Core/CommonHelpers/include/Acts/Tests/CommonHelpers/MeasurementCreatorHelper.hpp
+++ b/Tests/Core/CommonHelpers/include/Acts/Tests/CommonHelpers/MeasurementCreatorHelper.hpp
@@ -116,8 +116,6 @@ namespace Test {
                          const Surface*      surface,
                          resolution_t        lResolution) const
     {
-      auto geoID = surface->geoID();
-
       // Apply global to local
       Acts::Vector2D lPos;
       surface->globalToLocal(
-- 
GitLab