Fix bug in CKF and add script for KF timing test

Core/include/Acts/Fitter/KalmanFitter.hpp

@@ -189,6 +189,12 @@ class KalmanFitter {
  public:
   /// Shorthand definition
   using MeasurementSurfaces = std::multimap<const Layer*, const Surface*>;
+  /// The navigator type
+  using KalmanNavigator = typename propagator_t::Navigator;
+
+  /// The navigator has DirectNavigator type or not
+  static constexpr bool isDirectNavigator =
+      std::is_same<KalmanNavigator, DirectNavigator>::value;
 
   /// Default constructor is deleted
   KalmanFitter() = delete;
@@ -220,13 +226,6 @@ class KalmanFitter {
   /// Owned logging instance
   std::shared_ptr<const Logger> m_logger;
 
-  /// The navigator type
-  using KalmanNavigator = typename decltype(m_propagator)::Navigator;
-
-  /// The navigator has DirectNavigator type or not
-  static constexpr bool isDirectNavigator =
-      std::is_same<KalmanNavigator, DirectNavigator>::value;
-
   /// @brief Propagator Actor plugin for the KalmanFilter
   ///
   /// @tparam source_link_t is an type fulfilling the @c SourceLinkConcept
@@ -297,8 +296,9 @@ class KalmanFitter {
         // -> Get the measurement / calibrate
         // -> Create the predicted state
         // -> Perform the kalman update
-        // -> Check outlier behavior (@todo)
-        // -> Fill strack state information & update stepper information
+        // -> Check outlier behavior
+        // -> Fill strack state information & update stepper information if
+        // non-outlier
         if (state.stepping.navDir == forward and not result.smoothed and
             not result.forwardFiltered) {
           ACTS_VERBOSE("Perform forward filter step");
@@ -328,15 +328,15 @@ class KalmanFitter {
           if (backwardFiltering and not result.forwardFiltered) {
             ACTS_VERBOSE("Forward filtering done");
             result.forwardFiltered = true;
-            // Run backward filtering
+            // Start to run backward filtering:
             // Reverse navigation direction and reset navigation and stepping
             // state to last measurement
             ACTS_VERBOSE("Reverse navigation direction.");
             reverse(state, stepper, result);
           } else if (not result.smoothed) {
-            // -> Sort the track states (as now the path length is set)
-            // -> Call the smoothing
-            // -> Set a stop condition when all track states have been handled
+            // --> Search the starting state to run the smoothing
+            // --> Call the smoothing
+            // --> Set a stop condition when all track states have been handled
             ACTS_VERBOSE("Finalize/run smoothing");
             auto res = finalize(state, stepper, result);
             if (!res.ok()) {
@@ -802,34 +802,30 @@ class KalmanFitter {
       // Remember you smoothed the track states
       result.smoothed = true;
 
-      // Get the index of measurement states;
+      // Get the indices of measurement states;
       std::vector<size_t> measurementIndices;
-      auto lastState = result.fittedStates.getTrackState(result.trackTip);
-      if (lastState.typeFlags().test(Acts::TrackStateFlag::MeasurementFlag)) {
-        measurementIndices.push_back(result.trackTip);
-      }
+      measurementIndices.reserve(result.measurementStates);
       // Count track states to be smoothed
       size_t nStates = 0;
       result.fittedStates.applyBackwards(result.trackTip, [&](auto st) {
-        // Smoothing will start from the last measurement state
-        if (measurementIndices.empty()) {
-          // No smoothed parameter for the last few non-measurment states
+        bool isMeasurement =
+            st.typeFlags().test(TrackStateFlag::MeasurementFlag);
+        if (isMeasurement) {
+          measurementIndices.emplace_back(st.index());
+        } else if (measurementIndices.empty()) {
+          // No smoothed parameters if the last measurement state has not been
+          // found yet
           st.data().ismoothed = detail_lt::IndexData::kInvalid;
-        } else {
-          nStates++;
         }
-        size_t iprevious = st.previous();
-        if (iprevious != Acts::detail_lt::IndexData::kInvalid) {
-          auto previousState = result.fittedStates.getTrackState(iprevious);
-          if (previousState.typeFlags().test(
-                  Acts::TrackStateFlag::MeasurementFlag)) {
-            measurementIndices.push_back(iprevious);
-          }
+        // Start count when the last measurement state is found
+        if (not measurementIndices.empty()) {
+          nStates++;
         }
       });
       // Return error if the track has no measurement states (but this should
       // not happen)
       if (measurementIndices.empty()) {
+        ACTS_ERROR("Smoothing for a track without measurements.");
         return KalmanFitterError::SmoothFailed;
       }
       // Screen output for debugging
@@ -982,6 +978,7 @@ class KalmanFitter {
     auto result = m_propagator.template propagate(sParameters, kalmanOptions);
 
     if (!result.ok()) {
+      ACTS_ERROR("Propapation failed: " << result.error());
       return result.error();
     }
 
@@ -997,6 +994,7 @@ class KalmanFitter {
     }
 
     if (!kalmanResult.result.ok()) {
+      ACTS_ERROR("KalmanFilter failed: " << kalmanResult.result.error());
       return kalmanResult.result.error();
     }
 
@@ -1085,6 +1083,7 @@ class KalmanFitter {
     auto result = m_propagator.template propagate(sParameters, kalmanOptions);
 
     if (!result.ok()) {
+      ACTS_ERROR("Propapation failed: " << result.error());
       return result.error();
     }
 
@@ -1100,6 +1099,7 @@ class KalmanFitter {
     }
 
     if (!kalmanResult.result.ok()) {
+      ACTS_ERROR("KalmanFilter failed: " << kalmanResult.result.error());
       return kalmanResult.result.error();
     }
 

Core/include/Acts/TrackFinder/CombinatorialKalmanFilter.hpp

@@ -318,12 +318,19 @@ class CombinatorialKalmanFilter {
       // - Waiting for a current surface
       auto surface = state.navigation.currentSurface;
       if (surface != nullptr and not result.forwardFiltered) {
-        // Check if the surface is in the measurement map
-        // -> Get the measurement / calibrate
-        // -> Create the predicted state
+        // There are three scenarios:
+        // 1) The surface is in the measurement map
         // -> Select source links
-        // -> Perform the kalman update for selected source link
-        // -> Fill strack state information & update stepper information
+        // -> Perform the kalman update for selected non-outlier source links
+        // -> Add track states in multitrajectory. Multiple states mean branch
+        // splitting.
+        // -> Call branch stopper to justify each branch
+        // -> If there is non-outlier state, update stepper information
+        // 2) The surface is not in the measurement map but with material
+        // -> Add a hole or passive material state in multitrajectory
+        // -> Call branch stopper to justify the branch
+        // 3) The surface is neither in the measurement map nor with material
+        // -> Do nothing
         ACTS_VERBOSE("Perform filter step");
         auto res = filter(surface, state, stepper, result);
         if (!res.ok()) {
@@ -399,9 +406,10 @@ class CombinatorialKalmanFilter {
               ACTS_VERBOSE(
                   "Finalize/run smoothing for track with entry index = "
                   << result.trackTips.at(result.iSmoothed));
-              // -> Sort the track states (as now the path length is set)
-              // -> Call the smoothing
-              // -> Set a stop condition when all track states have been handled
+              // --> Search the starting state to run the smoothing
+              // --> Call the smoothing
+              // --> Set a stop condition when all track states have been
+              // handled
               auto res = finalize(state, stepper, result);
               if (!res.ok()) {
                 ACTS_ERROR("Error in finalize: " << res.error());
@@ -541,6 +549,17 @@ class CombinatorialKalmanFilter {
           return sourcelinkSelectionRes.error();
         }
 
+        // Retrieve the previous tip and its state
+        // The states created on this surface will have the common previous tip
+        size_t prevTip = SIZE_MAX;
+        TipState prevTipState;
+        if (not result.activeTips.empty()) {
+          prevTip = result.activeTips.back().first;
+          prevTipState = result.activeTips.back().second;
+          // New state is to be added. Remove the last tip from active tips
+          result.activeTips.erase(result.activeTips.end() - 1);
+        }
+
         // Remember the tip of the neighbor state on this surface
         size_t neighborTip = SIZE_MAX;
         // Loop over the selected source links
@@ -562,7 +581,8 @@ class CombinatorialKalmanFilter {
               sharedTip = index_it->second;
             }
           }
-          // Add a measurement/outlier track state proxy in multi trajectory
+
+          // The mask for adding a state in the multitrajectory
           // No storage allocation for:
           // -> predicted parameter and uncalibrated measurement if
           // already stored
@@ -575,10 +595,10 @@ class CombinatorialKalmanFilter {
               (isOutlier ? ~TrackStatePropMask::Filtered
                          : TrackStatePropMask::All);
 
-          // Add measurement or outlier track state to the multitrajectory
+          // Add measurement/outlier track state to the multitrajectory
           auto addStateRes = addSourcelinkState(
               stateMask, boundState, sourcelinks.at(index), isOutlier, result,
-              state.geoContext, neighborTip, sharedTip);
+              state.geoContext, prevTip, prevTipState, neighborTip, sharedTip);
           if (addStateRes.ok()) {
             const auto& [currentTip, tipState] = addStateRes.value();
             // Remember the track state tip for this stored source link
@@ -622,13 +642,12 @@ class CombinatorialKalmanFilter {
         // first, but could be changed later
         nBranchesOnSurface = 1;
 
-        // Retrieve the tip state and remove the last tip from active tips
+        // Retrieve the previous tip and its state
         size_t prevTip = SIZE_MAX;
         TipState tipState;
         if (not result.activeTips.empty()) {
           prevTip = result.activeTips.back().first;
           tipState = result.activeTips.back().second;
-          result.activeTips.erase(result.activeTips.end() - 1);
         }
 
         // The surface could be either sensitive or passive
@@ -639,11 +658,14 @@ class CombinatorialKalmanFilter {
           // Increment of number of passed sensitive surfaces
           tipState.nSensitiveSurfaces++;
         }
-        // Create state if there is already measurement detected on this branch
-        // For in-sensitive surface, only create state when smoothing is
+        // Add state if there is already measurement detected on this branch
+        // For in-sensitive surface, only add state when smoothing is
         // required
         if (tipState.nMeasurements > 0 and
             (isSensitive or (not isSensitive and smoothing))) {
+          // New state is to be added. Remove the last tip from active tips now
+          result.activeTips.erase(result.activeTips.end() - 1);
+
           // No source links on surface, add either hole or passive material
           // TrackState. No storage allocation for uncalibrated/calibrated
           // measurement and filtered parameter
@@ -723,15 +745,10 @@ class CombinatorialKalmanFilter {
         const TrackStatePropMask::Type& stateMask, const BoundState& boundState,
         const source_link_t& sourcelink, bool isOutlier, result_type& result,
         std::reference_wrapper<const GeometryContext> geoContext,
+        const size_t& prevTip, const TipState& prevTipState,
         size_t neighborTip = SIZE_MAX, size_t sharedTip = SIZE_MAX) const {
-      // Retrieve the tip state and remove the last tip from active tips
-      size_t prevTip = SIZE_MAX;
-      TipState tipState;
-      if (not result.activeTips.empty()) {
-        prevTip = result.activeTips.back().first;
-        tipState = result.activeTips.back().second;
-        result.activeTips.erase(result.activeTips.end() - 1);
-      }
+      // Inherit the tip state from the previous and will be updated later
+      TipState tipState = prevTipState;
 
       // Add a track state
       auto currentTip = result.fittedStates.addTrackState(stateMask, prevTip);
@@ -805,7 +822,7 @@ class CombinatorialKalmanFilter {
         // Increment number of measurements
         tipState.nMeasurements++;
       }
-      return std::make_pair(currentTip, tipState);
+      return std::make_pair(std::move(currentTip), std::move(tipState));
     }
 
     /// @brief CombinatorialKalmanFilter actor operation : add hole track state
@@ -961,7 +978,8 @@ class CombinatorialKalmanFilter {
                           result_type& result) const {
       // The tip of the track being smoothed
       const auto& currentTip = result.trackTips.at(result.iSmoothed);
-      // Get the index of measurement states;
+
+      // Get the indices of measurement states;
       std::vector<size_t> measurementIndices;
       // Count track states to be smoothed
       size_t nStates = 0;

Examples/Algorithms/GeneratorsPythia8/ACTFW/Options/Pythia8Options.cpp

@@ -72,7 +72,7 @@ FW::EventGenerator::Config FW::Options::readPythia8Options(
        Pythia8Generator::makeFunction(hard, lvl)},
       {PoissonMultiplicityGenerator{mu},
        GaussianVertexGenerator{{vtxStdXY, vtxStdXY, vtxStdZ, vtxStdT}},
-       Pythia8Generator::makeFunction(hard, lvl)},
+       Pythia8Generator::makeFunction(pileup, lvl)},
   };
   cfg.shuffle = vm["evg-shuffle"].as<bool>();
 

Examples/Scripts/KF_timing.py

@@ -0,0 +1,55 @@
+import ROOT
+import csv
+import matplotlib.pyplot as plt
+import numpy as np
+
+# Data preparation
+ptDict = {}
+
+# Open the output file
+with open('output.log', mode='r') as csv_file:
+    csv_reader = csv.reader(csv_file, delimiter=',')
+    # read lines and go for it
+    for csv_row in csv_reader : 
+        if len(csv_row) > 1 :                
+            # get the job id
+            jobID = csv_row[0]
+            # etabin, pt value, exec time
+            etabin = float(csv_row[1])
+            ptvalue = float(csv_row[2])
+            exectime = float(csv_row[3])
+
+            # Make sure you have all the keys ready
+            try :
+                pdict =  ptDict[ptvalue]
+            except :
+                ptDict[ptvalue] = {}
+                pdict =  ptDict[ptvalue]
+
+            # Now fill the sub dictionary
+            try :
+                vpdict = pdict[etabin]
+            except:
+                pdict[etabin] = []
+                vpdict = pdict[etabin]
+
+            vpdict += [ exectime ]
+
+# plot the ptDict
+plt.figure(figsize=(7, 5))
+
+ax = plt.subplot(111)
+plt.loglog(ptDict.keys(),[i[0][0] for i in np.array(list(ptDict.values()))],'.-', label='0<$\eta$<0.5')
+plt.loglog(ptDict.keys(),[i[1][0] for i in np.array(list(ptDict.values()))],'.-', label='0.5<$\eta$<1.0')
+plt.loglog(ptDict.keys(),[i[2][0] for i in np.array(list(ptDict.values()))],'.-', label='1.0<$\eta$<1.5')
+plt.loglog(ptDict.keys(),[i[3][0] for i in np.array(list(ptDict.values()))],'.-', label='1.5<$\eta$<2.0')
+plt.loglog(ptDict.keys(),[i[4][0] for i in np.array(list(ptDict.values()))],'.-', label='2.0<$\eta$<2.5')
+ax.set_xlabel('$p_T$ [GeV/c]')
+ax.set_ylabel('time/track [sec]')
+plt.yscale('log')
+ax.set_xlim((0.09,105))
+plt.legend(numpoints=1)
+
+plt.suptitle("KF timing vs. $p_T$")
+
+plt.show()