From cb36dffebc7d70a7d58f2ffc9a87b9b131fce70f Mon Sep 17 00:00:00 2001 From: xai Date: Thu, 23 Apr 2020 17:14:44 -0700 Subject: [PATCH 1/6] Clean-up the KalmanFitter --- Core/include/Acts/Fitter/KalmanFitter.hpp | 55 +++++++++---------- .../TrackFinder/CombinatorialKalmanFilter.hpp | 15 +++-- 2 files changed, 36 insertions(+), 34 deletions(-) diff --git a/Core/include/Acts/Fitter/KalmanFitter.hpp b/Core/include/Acts/Fitter/KalmanFitter.hpp index b01fde05390..bb844fa9752 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 decltype(m_propagator)::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 @@ -328,15 +327,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 +801,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 parameter if the last measurment 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 +977,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 +993,7 @@ class KalmanFitter { } if (!kalmanResult.result.ok()) { + ACTS_ERROR("KalmanFilter failed: " << kalmanResult.result.error()); return kalmanResult.result.error(); } @@ -1085,6 +1082,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 +1098,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..5478fe0b6f3 100644 --- a/Core/include/Acts/TrackFinder/CombinatorialKalmanFilter.hpp +++ b/Core/include/Acts/TrackFinder/CombinatorialKalmanFilter.hpp @@ -399,9 +399,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()); @@ -562,7 +563,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,7 +577,7 @@ 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); @@ -961,7 +963,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; From d7cddf69f917f8446e1678692490c5b2dad59c0b Mon Sep 17 00:00:00 2001 From: xai Date: Thu, 23 Apr 2020 23:43:13 -0700 Subject: [PATCH 2/6] Fix the pythia8 option --- .../GeneratorsPythia8/ACTFW/Options/Pythia8Options.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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(); From 84a8f40136e8ec23c0baf9b8ede6f145d00cfd94 Mon Sep 17 00:00:00 2001 From: xai Date: Thu, 23 Apr 2020 23:43:34 -0700 Subject: [PATCH 3/6] Add the script for KF timing test --- Examples/Scripts/KF_timing.sh | 119 ++++++++++++++++++++++++++++++++++ 1 file changed, 119 insertions(+) create mode 100644 Examples/Scripts/KF_timing.sh diff --git a/Examples/Scripts/KF_timing.sh b/Examples/Scripts/KF_timing.sh new file mode 100644 index 00000000000..28f8a1b178c --- /dev/null +++ b/Examples/Scripts/KF_timing.sh @@ -0,0 +1,119 @@ +#!/bin/bash +# +# This script runs the KalmanFitter timing test with different smoothing method at different |eta| or pT bins +# using particle gun particles +# +helpFunction() +{ + echo "" + echo "Usage: $0 -v -d -b -t -n " + echo -e "\t-v The binning variable, either 'eta' or 'pt'. Required" + echo -e "\t-d The detector type, either 'Generic' or 'DD4hep'. Required" + 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 +numTracksPerEvent=100 +numEvents=1 + +# Get parameters +while getopts "v:d:b:t:n" opt +do + case "$opt" in + v ) binVal="$OPTARG" ;; + d ) detector="$OPTARG" ;; + b ) bFieldMap="$OPTARG" ;; + t ) numTracksPerEvent="$OPTARG" ;; + n ) numEvents="$OPTARG" ;; + ? ) helpFunction ;; # Print helpFunction in case unknown parameter + esac +done + +#Check input for binning variable +if [ -z "${binVal}" ]; then + echo "binVal is empty"; + helpFunction +elif [[ ${binVal} -ne "eta" ]] && [[ ${binVal} -ne "pt" ]]; then + echo "binVal is not valid. It should be either 'eta' or 'pt'!" + exit 1 +fi + +# Check input for detector +if [ -z "${detector}" ]; then + echo "detector is empty"; + helpFunction +fi + +# 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 "Binning Variable: ${binVal}" +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_directory=$PWD +run_directory=KalmanFitter_timing_${time_stamp} +mkdir ${run_directory} +cd ${run_directory} + +output_file='output.log' +echo "*****KalmanFitter timing test vs. ${binVal}*****" > ${output_file} +echo "* Test Detector: ${detector}" >> ${output_file} +echo "* B Field: ${bField}" >> ${output_file} +echo "* Events: ${numEvents}" >> ${output_file} +echo "* Tracks/event: ${numTracksPerEvent}" >> ${output_file} +echo "****************************************" >> ${output_file} +echo "*" +echo "* job | ${bFieldMap} " >> ${output_file} + +jobID=0 +if [[ ${binVal} -eq "eta" ]];then + # Loop over the eta bins + for eta in 0 0.5 1 1.5 2 2.5 ; do + # Run simulation + simulation="${exe_directory}/ActsSimFatras${detector} ${bField} -n ${numEvents} --evg-input-type gun --pg-nparticles ${numTracksPerEvent} --pg-pt-range 0.1 100 --pg-eta-range ${eta} ${eta} --fatras-pmin-gev 0.1 --output-csv=1 --output-dir sim_${detector}_e${numEvents}_t${numTracksPerEvent}_eta${eta}" + echo ${simulation} + eval ${simulation} + # Run reconstruction + reconstruction="$exe_directory/ActsRecTruthTracks ${bField} --input-dir sim_${detector}_e${numEvents}_t${numTracksPerEvent}_eta${eta} --output-dir reco_${detector}_e${numEvents}_t${numTracksPerEvent}_eta${eta}" + eval ${reconstruction} + + echo "${jobID}, ${eta}" >> ${output_file} + # Archive with Job ID + mv reco_${detector}_e${numEvents}_t${numTracksPerEvent}_eta${eta}/timing.tsv timing_${jobID}.tsv + # JobID + let "jobID++" + done +else + # Loop over the pt bins + for pt in 0.5 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 50.0 100.0 ; do + # Run simulation + simulation="${exe_directory}/ActsSimFatras${detector} ${bField} -n ${numEvents} --evg-input-type gun --pg-nparticles ${numTracksPerEvent} --pg-pt-range ${pt} ${pt} --pg-eta-range -2.5 2.5 --fatras-pmin-gev 0.1 --output-csv=1 --output-dir sim_${detector}_e${numEvents}_t${numTracksPerEvent}_pt${pt}" + eval ${simulation} + # Run reconstruction + reconstruction="$exe_directory/ActsRecTruthTracks ${bField} --input-dir sim_${detector}_e${numEvents}_t${numTracksPerEvent}_pt${pt} --output-dir reco_${detector}_e${numEvents}_t${numTracksPerEvent}_pt${pt}" + eval ${reconstruction} + + echo "${jobID}, ${pt}" >> ${output_file} + # Archive with Job ID + mv reco_${detector}_e${numEvents}_t${numTracksPerEvent}_pt${pt}/timing.tsv timing_${jobID}.tsv + # JobID + let "jobID++" + done +fi From dbb0c9e6888625b597ba29f62c59f27acb64b812 Mon Sep 17 00:00:00 2001 From: xai Date: Fri, 24 Apr 2020 13:35:24 -0700 Subject: [PATCH 4/6] Add the python script for plotting KF timing test result --- Core/include/Acts/Fitter/KalmanFitter.hpp | 2 +- Examples/Scripts/KF_timing.py | 55 +++++++++++ Examples/Scripts/KF_timing.sh | 115 +++++++++------------- 3 files changed, 104 insertions(+), 68 deletions(-) create mode 100644 Examples/Scripts/KF_timing.py diff --git a/Core/include/Acts/Fitter/KalmanFitter.hpp b/Core/include/Acts/Fitter/KalmanFitter.hpp index bb844fa9752..94d963f89d4 100644 --- a/Core/include/Acts/Fitter/KalmanFitter.hpp +++ b/Core/include/Acts/Fitter/KalmanFitter.hpp @@ -190,7 +190,7 @@ class KalmanFitter { /// Shorthand definition using MeasurementSurfaces = std::multimap; /// The navigator type - using KalmanNavigator = typename decltype(m_propagator)::Navigator; + using KalmanNavigator = typename propagator_t::Navigator; /// The navigator has DirectNavigator type or not static constexpr bool isDirectNavigator = 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 index 28f8a1b178c..1e184b00903 100644 --- a/Examples/Scripts/KF_timing.sh +++ b/Examples/Scripts/KF_timing.sh @@ -1,14 +1,15 @@ #!/bin/bash # -# This script runs the KalmanFitter timing test with different smoothing method at different |eta| or pT bins +# This script runs the test of KalmanFitter timing vs. pT at different eta # using particle gun particles -# -helpFunction() +# 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 -v -d -b -t -n " - echo -e "\t-v The binning variable, either 'eta' or 'pt'. Required" - echo -e "\t-d The detector type, either 'Generic' or 'DD4hep'. Required" + 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" @@ -16,37 +17,21 @@ helpFunction() } #Default number of tracks/event and event +detector=Generic numTracksPerEvent=100 numEvents=1 - # Get parameters -while getopts "v:d:b:t:n" opt +while getopts "d:b:t:n" opt do case "$opt" in - v ) binVal="$OPTARG" ;; d ) detector="$OPTARG" ;; b ) bFieldMap="$OPTARG" ;; t ) numTracksPerEvent="$OPTARG" ;; n ) numEvents="$OPTARG" ;; - ? ) helpFunction ;; # Print helpFunction in case unknown parameter + ? ) help ;; # Print help in case unknown parameter esac done -#Check input for binning variable -if [ -z "${binVal}" ]; then - echo "binVal is empty"; - helpFunction -elif [[ ${binVal} -ne "eta" ]] && [[ ${binVal} -ne "pt" ]]; then - echo "binVal is not valid. It should be either 'eta' or 'pt'!" - exit 1 -fi - -# Check input for detector -if [ -z "${detector}" ]; then - echo "detector is empty"; - helpFunction -fi - # Check input for B field bField='--bf-value 0 0 2' if [ -z "${bFieldMap}" ]; then @@ -60,60 +45,56 @@ else fi #Print the parameters -echo "Binning Variable: ${binVal}" 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_directory=$PWD -run_directory=KalmanFitter_timing_${time_stamp} -mkdir ${run_directory} -cd ${run_directory} +exe_dir=$PWD +run_dir=KalmanFitter_timing_${time_stamp} +mkdir ${run_dir} +cd ${run_dir} output_file='output.log' -echo "*****KalmanFitter timing test vs. ${binVal}*****" > ${output_file} -echo "* Test Detector: ${detector}" >> ${output_file} -echo "* B Field: ${bField}" >> ${output_file} -echo "* Events: ${numEvents}" >> ${output_file} -echo "* Tracks/event: ${numTracksPerEvent}" >> ${output_file} +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 | ${bFieldMap} " >> ${output_file} +echo "* job | eta | pt | fit_time_per_event" >> ${output_file} jobID=0 -if [[ ${binVal} -eq "eta" ]];then - # Loop over the eta bins - for eta in 0 0.5 1 1.5 2 2.5 ; do - # Run simulation - simulation="${exe_directory}/ActsSimFatras${detector} ${bField} -n ${numEvents} --evg-input-type gun --pg-nparticles ${numTracksPerEvent} --pg-pt-range 0.1 100 --pg-eta-range ${eta} ${eta} --fatras-pmin-gev 0.1 --output-csv=1 --output-dir sim_${detector}_e${numEvents}_t${numTracksPerEvent}_eta${eta}" - echo ${simulation} - eval ${simulation} - # Run reconstruction - reconstruction="$exe_directory/ActsRecTruthTracks ${bField} --input-dir sim_${detector}_e${numEvents}_t${numTracksPerEvent}_eta${eta} --output-dir reco_${detector}_e${numEvents}_t${numTracksPerEvent}_eta${eta}" - eval ${reconstruction} - - echo "${jobID}, ${eta}" >> ${output_file} - # Archive with Job ID - mv reco_${detector}_e${numEvents}_t${numTracksPerEvent}_eta${eta}/timing.tsv timing_${jobID}.tsv - # JobID - let "jobID++" - done -else # Loop over the pt bins - for pt in 0.5 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 50.0 100.0 ; do - # Run simulation - simulation="${exe_directory}/ActsSimFatras${detector} ${bField} -n ${numEvents} --evg-input-type gun --pg-nparticles ${numTracksPerEvent} --pg-pt-range ${pt} ${pt} --pg-eta-range -2.5 2.5 --fatras-pmin-gev 0.1 --output-csv=1 --output-dir sim_${detector}_e${numEvents}_t${numTracksPerEvent}_pt${pt}" - eval ${simulation} - # Run reconstruction - reconstruction="$exe_directory/ActsRecTruthTracks ${bField} --input-dir sim_${detector}_e${numEvents}_t${numTracksPerEvent}_pt${pt} --output-dir reco_${detector}_e${numEvents}_t${numTracksPerEvent}_pt${pt}" - eval ${reconstruction} + 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} - echo "${jobID}, ${pt}" >> ${output_file} - # Archive with Job ID - mv reco_${detector}_e${numEvents}_t${numTracksPerEvent}_pt${pt}/timing.tsv timing_${jobID}.tsv - # JobID - let "jobID++" + # JobID + let "jobID++" + done done -fi From 70e9d10a170eeaabc2fa890609e0a736fe3054db Mon Sep 17 00:00:00 2001 From: xai Date: Sat, 25 Apr 2020 14:26:57 -0700 Subject: [PATCH 5/6] Fix a bug in CKF --- Core/include/Acts/Fitter/KalmanFitter.hpp | 5 +- .../TrackFinder/CombinatorialKalmanFilter.hpp | 53 ++++++++++++------- 2 files changed, 37 insertions(+), 21 deletions(-) diff --git a/Core/include/Acts/Fitter/KalmanFitter.hpp b/Core/include/Acts/Fitter/KalmanFitter.hpp index 94d963f89d4..0fc7991ca2a 100644 --- a/Core/include/Acts/Fitter/KalmanFitter.hpp +++ b/Core/include/Acts/Fitter/KalmanFitter.hpp @@ -296,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"); diff --git a/Core/include/Acts/TrackFinder/CombinatorialKalmanFilter.hpp b/Core/include/Acts/TrackFinder/CombinatorialKalmanFilter.hpp index 5478fe0b6f3..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()) { @@ -542,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 @@ -580,7 +598,7 @@ class CombinatorialKalmanFilter { // 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 @@ -624,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 @@ -641,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 @@ -725,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); @@ -807,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 From 3038d5f4838f4f6bc6d143b29df6f3fddd66f760 Mon Sep 17 00:00:00 2001 From: xai Date: Mon, 27 Apr 2020 11:48:38 -0700 Subject: [PATCH 6/6] Fix the comment --- Core/include/Acts/Fitter/KalmanFitter.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Core/include/Acts/Fitter/KalmanFitter.hpp b/Core/include/Acts/Fitter/KalmanFitter.hpp index 0fc7991ca2a..0fa506367f2 100644 --- a/Core/include/Acts/Fitter/KalmanFitter.hpp +++ b/Core/include/Acts/Fitter/KalmanFitter.hpp @@ -813,7 +813,7 @@ class KalmanFitter { if (isMeasurement) { measurementIndices.emplace_back(st.index()); } else if (measurementIndices.empty()) { - // No smoothed parameter if the last measurment state has not been + // No smoothed parameters if the last measurement state has not been // found yet st.data().ismoothed = detail_lt::IndexData::kInvalid; }