diff --git a/Core/include/Acts/Fitter/KalmanFitter.hpp b/Core/include/Acts/Fitter/KalmanFitter.hpp index b01fde05390..0fa506367f2 100644 --- a/Core/include/Acts/Fitter/KalmanFitter.hpp +++ b/Core/include/Acts/Fitter/KalmanFitter.hpp @@ -189,6 +189,12 @@ class KalmanFitter { public: /// Shorthand definition using MeasurementSurfaces = std::multimap; + /// The navigator type + using KalmanNavigator = typename propagator_t::Navigator; + + /// The navigator has DirectNavigator type or not + static constexpr bool isDirectNavigator = + std::is_same::value; /// Default constructor is deleted KalmanFitter() = delete; @@ -220,13 +226,6 @@ class KalmanFitter { /// Owned logging instance std::shared_ptr 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::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 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(); } diff --git a/Core/include/Acts/TrackFinder/CombinatorialKalmanFilter.hpp b/Core/include/Acts/TrackFinder/CombinatorialKalmanFilter.hpp index 517a36ea340..f9e70bb48c3 100644 --- a/Core/include/Acts/TrackFinder/CombinatorialKalmanFilter.hpp +++ b/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 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 measurementIndices; // Count track states to be smoothed size_t nStates = 0; diff --git a/Examples/Algorithms/GeneratorsPythia8/ACTFW/Options/Pythia8Options.cpp b/Examples/Algorithms/GeneratorsPythia8/ACTFW/Options/Pythia8Options.cpp index 47ceba8ceaf..c0e2bfd983e 100644 --- a/Examples/Algorithms/GeneratorsPythia8/ACTFW/Options/Pythia8Options.cpp +++ b/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(); diff --git a/Examples/Scripts/KF_timing.py b/Examples/Scripts/KF_timing.py new file mode 100644 index 00000000000..b1b62f950f5 --- /dev/null +++ b/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() diff --git a/Examples/Scripts/KF_timing.sh b/Examples/Scripts/KF_timing.sh new file mode 100644 index 00000000000..1e184b00903 --- /dev/null +++ b/Examples/Scripts/KF_timing.sh @@ -0,0 +1,100 @@ +#!/bin/bash +# +# This script runs the test of KalmanFitter timing vs. pT at different eta +# using particle gun particles +# To run the test:./KF_timing.sh -d -b -t -n +# In default, it will run with Generic detector in constant B field using 1 event and 100 tracks/event + +help() +{ + echo "" + echo "Usage: $0 -d -b -t -n " + echo -e "\t-d The detector type, either 'Generic' or 'DD4hep'. Optional. In default 'Generic'" + echo -e "\t-b The '.txt' or '.root' file for B Field map. Optional. In default using constant BField: (0, 0, 2)" + echo -e "\t-t The number of tracks per event. Optional. In default: 100" + echo -e "\t-n The number of events. Optional. In default: 1" + exit 1 # Exit script after printing help +} + +#Default number of tracks/event and event +detector=Generic +numTracksPerEvent=100 +numEvents=1 +# Get parameters +while getopts "d:b:t:n" opt +do + case "$opt" in + d ) detector="$OPTARG" ;; + b ) bFieldMap="$OPTARG" ;; + t ) numTracksPerEvent="$OPTARG" ;; + n ) numEvents="$OPTARG" ;; + ? ) help ;; # Print help in case unknown parameter + esac +done + +# Check input for B field +bField='--bf-value 0 0 2' +if [ -z "${bFieldMap}" ]; then + echo "bFieldMap is empty. Will use constant B Field (0, 0, 2) then." +elif [ -f "$bFieldMap" ] ; then + echo "Input bField map: $bFieldMap " + bField='--bf-map ${bFieldMap}' +else + echo "Input file for bField map file does not exist!" + exit 1 +fi + +#Print the parameters +echo "Test detector: ${detector}" +echo "B Field: ${bField}" +echo "Number of tracks per event: ${numTracksPerEvent}" +echo "Number of events: ${numEvents}" + +time_stamp=`date +%F%T` +exe_dir=$PWD +run_dir=KalmanFitter_timing_${time_stamp} +mkdir ${run_dir} +cd ${run_dir} + +output_file='output.log' +echo "*****KalmanFitter timing test vs. pT*****" > ${output_file} +echo "Test Detector: ${detector}" >> ${output_file} +echo "BField: ${bField}" >> ${output_file} +echo "Events: ${numEvents}" >> ${output_file} +echo "Tracks_per_event: ${numTracksPerEvent}" >> ${output_file} +echo "****************************************" >> ${output_file} +echo "*" +echo "* job | eta | pt | fit_time_per_event" >> ${output_file} + +jobID=0 + # Loop over the pt bins + for pt in 0.1 0.5 1.0 2.0 3.0 4.0 5.0 8.0 10.0 50.0 100.0 ; do + # Loop over the eta bin number + for etaBin in 0 1 2 3 4; do + etaLow=$(echo "${etaBin}*0.5"|bc) + etaUp=$(echo "${etaBin}*0.5 + 0.5"|bc) + eta=$(echo "${etaBin}*0.5 + 0.25"|bc) + + # Run sim + sim="${exe_dir}/ActsSimFatras${detector} ${bField} -n ${numEvents} --evg-input-type gun --pg-nparticles ${numTracksPerEvent} --pg-pt-range ${pt} ${pt} --pg-eta-range ${etaLow} ${etaUp} --fatras-pmin-gev 0.1 --output-csv=1 --output-dir sim_${detector}_e${numEvents}_t${numTracksPerEvent}_eta${eta}_pt${pt}" + echo "Run sim with '${sim}'" + eval ${sim} + + # Run reco + reco="$exe_dir/ActsRecTruthTracks ${bField} --input-dir sim_${detector}_e${numEvents}_t${numTracksPerEvent}_eta${eta}_pt${pt} --output-dir reco_${detector}_e${numEvents}_t${numTracksPerEvent}_eta${eta}_pt${pt}" + echo "Run reco with '${reco}'" + eval ${reco} + + # Archive with Job ID + mv reco_${detector}_e${numEvents}_t${numTracksPerEvent}_eta${eta}_pt${pt}/timing.tsv timing_${jobID}.tsv + # Extract the fitting time + fit_time_str=`grep "Algorithm:FittingAlgorithm" timing_${jobID}.tsv | awk '{print $3}'` + # Make sure the fit time is fixed-point for calculation with bc + fit_time_per_event=$(echo ${fit_time_str} | awk '{printf("%.10f\n", $1)}') + fit_time_per_track=$(echo "${fit_time_per_event}/${numTracksPerEvent}"|bc -l) + echo "${jobID}, ${etaBin}, ${pt}, ${fit_time_per_track}" >> ${output_file} + + # JobID + let "jobID++" + done + done