diff --git a/.github/workflows/builds.yml b/.github/workflows/builds.yml index 32072613347..b42793a85b8 100644 --- a/.github/workflows/builds.yml +++ b/.github/workflows/builds.yml @@ -22,115 +22,6 @@ env: CCACHE_KEY_SUFFIX: r1 jobs: - lcg: - runs-on: ubuntu-latest - container: ghcr.io/acts-project/${{ matrix.image }}:v41 - strategy: - matrix: - image: - - centos7-lcg100-gcc10 - - centos7-lcg101-gcc11 - - centos8-lcg100-gcc10 - - centos8-lcg101-gcc11 - env: - SETUP: source /opt/lcg_view/setup.sh - INSTALL_DIR: ${{ github.workspace }}/install - ACTS_LOG_FAILURE_THRESHOLD: WARNING - steps: - - uses: actions/checkout@v3 - - - name: Restore ccache - uses: actions/cache/restore@v3 - id: ccache-restore - with: - path: ${{ github.workspace }}/ccache - key: ${{ runner.os }}-ccache-${{ matrix.image }}_${{ env.CCACHE_KEY_SUFFIX }}_${{ github.sha }} - restore-keys: | - ${{ runner.os }}-ccache-${{ matrix.image }}_${{ env.CCACHE_KEY_SUFFIX }}_ - - - name: Configure - # setting CMAKE_CXX_STANDARD=17 is a workaround for a bug in the - # dd4hep CMake configuration that gets triggered on recent CMake - # versions - run: > - ${SETUP} && - ccache -z && - cmake -B build -S . - -GNinja - -DCMAKE_CXX_COMPILER_LAUNCHER=ccache - -DCMAKE_BUILD_TYPE=Release - -DCMAKE_CXX_FLAGS=-Werror - -DCMAKE_CXX_STANDARD=17 - -DCMAKE_INSTALL_PREFIX="${INSTALL_DIR}" - -DACTS_LOG_FAILURE_THRESHOLD=WARNING - -DACTS_BUILD_EXAMPLES_PYTHON_BINDINGS=ON - -DACTS_FORCE_ASSERTIONS=ON - -DACTS_BUILD_EXAMPLES=ON - -DACTS_BUILD_PLUGIN_DD4HEP=OFF - -DACTS_BUILD_PLUGIN_TGEO=ON - -DACTS_BUILD_PLUGIN_IDENTIFICATION=ON - -DACTS_BUILD_PLUGIN_JSON=ON - -DACTS_BUILD_FATRAS=ON - -DACTS_BUILD_PLUGIN_LEGACY=ON - -DACTS_BUILD_PLUGIN_AUTODIFF=OFF - -DACTS_BUILD_BENCHMARKS=ON - -DACTS_BUILD_UNITTESTS=ON - -DACTS_BUILD_INTEGRATIONTESTS=ON - -DACTS_BUILD_EXAMPLES_DD4HEP=OFF - -DACTS_BUILD_PLUGIN_EDM4HEP=OFF - -DACTS_BUILD_EXAMPLES_GEANT4=ON - -DACTS_BUILD_EXAMPLES_HEPMC3=ON - -DACTS_BUILD_EXAMPLES_PYTHIA8=ON - -DACTS_BUILD_FATRAS_GEANT4=ON - -DACTS_BUILD_FATRAS=ON - -DACTS_BUILD_ALIGNMENT=ON - -DACTS_BUILD_ANALYSIS_APPS=ON - -DACTS_BUILD_PLUGIN_ACTSVG=ON - - - name: Build - run: ${SETUP} && cmake --build build - - - name: ccache stats - run: ${SETUP} && ccache -s - - - name: Save ccache - uses: actions/cache/save@v3 - if: always() - with: - path: ${{ github.workspace }}/ccache - key: ${{ steps.ccache-restore.outputs.cache-primary-key }} - - - name: Unit tests - run: ${SETUP} && cmake --build build --target test - - - name: Integration tests - run: ${SETUP} && cmake --build build --target integrationtests - - - name: Install - run: ${SETUP} && cmake --build build --target install - - - uses: actions/upload-artifact@v3 - with: - name: acts-${{ matrix.image }} - path: ${{ env.INSTALL_DIR }} - - - name: Downstream configure - run: > - ${SETUP} && - cmake -B build-downstream -S Tests/DownstreamProject - -GNinja - -DCMAKE_BUILD_TYPE=Release - -DCMAKE_CXX_FLAGS=-Werror - -DCMAKE_CXX_STANDARD=17 - -DCMAKE_PREFIX_PATH="${INSTALL_DIR}" - -DDD4HEP=OFF - - - name: Downstream build - run: ${SETUP} && cmake --build build-downstream - - - name: Downstream run - run: ${SETUP} && ./build-downstream/bin/ShowActsVersion - linux_ubuntu: runs-on: ubuntu-latest container: ghcr.io/acts-project/ubuntu2204:v41 diff --git a/CI/physmon/phys_perf_mon.sh b/CI/physmon/phys_perf_mon.sh index 8b4e023be7c..fc6cc6f3d3f 100755 --- a/CI/physmon/phys_perf_mon.sh +++ b/CI/physmon/phys_perf_mon.sh @@ -1,6 +1,7 @@ #!/bin/bash set -e +set -x mode=${1:all} @@ -296,6 +297,22 @@ if [[ "$mode" == "all" || "$mode" == "fullchains" ]]; then --config CI/physmon/vertexing_config.yml ec=$(($ec | $?)) + Examples/Scripts/generic_plotter.py \ + $outdir/tracksummary_ckf_ttbar.root \ + tracksummary \ + $outdir/tracksummary_ckf_ttbar_hist.root \ + --config CI/physmon/tracksummary_ckf_config.yml + ec=$(($ec | $?)) + + # remove ntuple file because it's large + rm $outdir/tracksummary_ckf_ttbar.root + + run_histcmp \ + $outdir/tracksummary_ckf_ttbar_hist.root \ + $refdir/tracksummary_ckf_ttbar_hist.root \ + "Track Summary CKF ttbar" \ + tracksummary_ckf_ttbar + # remove ntuple file because it's large rm $outdir/performance_amvf_ttbar.root diff --git a/CI/physmon/reference/tracksummary_ckf_ttbar_hist.root b/CI/physmon/reference/tracksummary_ckf_ttbar_hist.root new file mode 100644 index 00000000000..d9592b163a9 Binary files /dev/null and b/CI/physmon/reference/tracksummary_ckf_ttbar_hist.root differ diff --git a/Core/include/Acts/EventData/ParticleHypothesis.hpp b/Core/include/Acts/EventData/ParticleHypothesis.hpp index 994f5179e85..3fb33b27664 100644 --- a/Core/include/Acts/EventData/ParticleHypothesis.hpp +++ b/Core/include/Acts/EventData/ParticleHypothesis.hpp @@ -110,6 +110,9 @@ class NonNeutralChargedParticleHypothesis pion().mass(), absQ); } + static NonNeutralChargedParticleHypothesis chargedGeantino() { + return chargedGeantino(Acts::UnitConstants::e); + } static NonNeutralChargedParticleHypothesis chargedGeantino(float absQ) { return NonNeutralChargedParticleHypothesis(PdgParticle::eInvalid, 0, absQ); } @@ -154,6 +157,9 @@ class ParticleHypothesis : public GenericParticleHypothesis { static ParticleHypothesis geantino() { return NeutralParticleHypothesis::geantino(); } + static ParticleHypothesis chargedGeantino() { + return chargedGeantino(Acts::UnitConstants::e); + } static ParticleHypothesis chargedGeantino(float absQ) { return ParticleHypothesis(PdgParticle::eInvalid, 0, absQ); } diff --git a/Core/include/Acts/Vertexing/AdaptiveGridDensityVertexFinder.hpp b/Core/include/Acts/Vertexing/AdaptiveGridDensityVertexFinder.hpp index 4e0085517f1..4896c7c6757 100644 --- a/Core/include/Acts/Vertexing/AdaptiveGridDensityVertexFinder.hpp +++ b/Core/include/Acts/Vertexing/AdaptiveGridDensityVertexFinder.hpp @@ -77,7 +77,7 @@ class AdaptiveGridDensityVertexFinder { /// /// Only needed if cacheGridStateForTrackRemoval == true struct State { - // Map from z from the z bin values to the corresponding track density + // Map from the z bin values to the corresponding track density DensityMap mainDensityMap; // Map from input track to corresponding track density map diff --git a/Core/include/Acts/Vertexing/AdaptiveGridDensityVertexFinder.ipp b/Core/include/Acts/Vertexing/AdaptiveGridDensityVertexFinder.ipp index bc0737d0ad7..f33cf2947a7 100644 --- a/Core/include/Acts/Vertexing/AdaptiveGridDensityVertexFinder.ipp +++ b/Core/include/Acts/Vertexing/AdaptiveGridDensityVertexFinder.ipp @@ -60,22 +60,22 @@ auto Acts::AdaptiveGridDensityVertexFinder::find( if (not state.mainDensityMap.empty()) { if (not m_cfg.estimateSeedWidth) { // Get z value of highest density bin - auto maxZres = m_cfg.gridDensity.getMaxZPosition(state.mainDensityMap); + auto maxZTRes = m_cfg.gridDensity.getMaxZTPosition(state.mainDensityMap); - if (!maxZres.ok()) { - return maxZres.error(); + if (!maxZTRes.ok()) { + return maxZTRes.error(); } - z = *maxZres; + z = (*maxZTRes).first; } else { // Get z value of highest density bin and width - auto maxZres = - m_cfg.gridDensity.getMaxZPositionAndWidth(state.mainDensityMap); + auto maxZTResAndWidth = + m_cfg.gridDensity.getMaxZTPositionAndWidth(state.mainDensityMap); - if (!maxZres.ok()) { - return maxZres.error(); + if (!maxZTResAndWidth.ok()) { + return maxZTResAndWidth.error(); } - z = (*maxZres).first; - width = (*maxZres).second; + z = (*maxZTResAndWidth).first.first; + width = (*maxZTResAndWidth).second; } } diff --git a/Core/include/Acts/Vertexing/AdaptiveGridTrackDensity.hpp b/Core/include/Acts/Vertexing/AdaptiveGridTrackDensity.hpp index c6befda103f..31c53af0d64 100644 --- a/Core/include/Acts/Vertexing/AdaptiveGridTrackDensity.hpp +++ b/Core/include/Acts/Vertexing/AdaptiveGridTrackDensity.hpp @@ -13,41 +13,77 @@ #include +#include + namespace Acts { /// @class AdaptiveGridTrackDensity -/// @brief Implements a 1-dim density grid to be filled with -/// track Gaussian distributions. Each single track is modelled -/// as a 2-dim Gaussian distribution grid in the d0-z0 plane, -/// but only the overlap with the z-axis (i.e. a 1-dim density -/// vector) needs to be calculated. +/// @brief Implements a 1D (no time seeding) / 2D (time seeding) +/// grid that is filled with track densities. +/// Each track is modelled by a 2D / 3D Gaussian distribution in the +/// d0-z0 / d0-z0-t0 plane, which is evaluated at d0=0. Therefore, +/// each track effectively lives in 1D / 2D. /// The position of the highest track density (of either a single /// bin or the sum of a certain region) can be determined. /// Single tracks can be cached and removed from the overall density. -/// Unlike the GaussianGridTrackDensity, the overall density vector -/// grows adaptively with the tracks densities being added to the grid. +/// Unlike in the GaussianGridTrackDensity, the overall density map +/// grows adaptively when tracks densities are added to the grid. /// -/// @tparam trkGridSize The 2-dim grid size of a single track, i.e. -/// a single track is modelled as a (trkGridSize x trkGridSize) grid -/// in the d0-z0 plane. Note: trkGridSize has to be an odd value. -template +/// @tparam spatialTrkGridSize Number of bins per track in z direction +/// @tparam temporalTrkGridSize Number of bins per track in t direction +/// @note In total, a track is represented by a grid of size +/// spatialTrkGridSize * temporalTrkGridSize +template class AdaptiveGridTrackDensity { - // Assert odd trkGridSize - static_assert(trkGridSize % 2); + // Assert odd spatial and temporal track grid size + static_assert(spatialTrkGridSize % 2); + static_assert(temporalTrkGridSize % 2); public: - using DensityMap = std::unordered_map; + // The first (second) integer indicates the bin's z (t) position + using Bin = std::pair; + // Mapping between bins and track densities + using DensityMap = std::unordered_map>; + // Coordinates in the z-t plane; the t value will be set to 0 if time + // vertex seeding is disabled + using ZTPosition = std::pair; + // z-t position of a maximum and its width + using ZTPositionAndWidth = std::pair; /// The configuration struct struct Config { - /// @param binSize_ The binSize in mm - Config(float binSize_ = 0.1) : binSize(binSize_) {} - - // Z size of one single bin in grid - float binSize; // mm + /// @param spatialBinExtent_ The spatial extent of a bin in mm + Config(float spatialBinExtent_) : spatialBinExtent(spatialBinExtent_) { + if constexpr (temporalTrkGridSize > 1) { + throw std::invalid_argument( + "temporalBinExtent must be provided if temporalTrkGridSize > 1 " + "(i.e., if time vertex seeding is enabled)."); + } + } + + /// @param spatialBinExtent_ The spatial extent of a bin in mm + /// @param temporalBinExtent_ The temporal extent of a bin in mm + /// @note The speed of light is set to 1, hence the unit. + Config(float spatialBinExtent_, float temporalBinExtent_) + : spatialBinExtent(spatialBinExtent_), + temporalBinExtent(temporalBinExtent_) { + if constexpr (temporalTrkGridSize == 1) { + throw std::invalid_argument( + "temporalBinExtent must not be provided if temporalTrkGridSize == " + "1 (i.e., if time vertex seeding is disabled)."); + } + } + + // Spatial extent of a bin in d0 and z0 direction, should always be set to a + // positive value + float spatialBinExtent = 0.; // mm + + // Temporal extent of a bin, should be set to 0 if time vertex seeding is + // disabled (i.e., if temporalTrkGridSize = 1) + float temporalBinExtent = 0.; // mm // Do NOT use just the z-bin with the highest - // track density, but instead check the (up to) + // track density, but instead check (up to) // first three density maxima (only those that have // a maximum relative deviation of 'relativeDensityDev' // from the main maximum) and take the z-bin of the @@ -64,100 +100,120 @@ class AdaptiveGridTrackDensity { /// @brief Calculates the bin center from the bin number /// @param bin Bin number + /// @param binExtent Bin extent /// @return Bin center - float getBinCenter(int bin) const; + static float getBinCenter(int bin, float binExtent); - /// @brief Calculates the bin number corresponding to a d or z value - /// @param value d or z value + /// @brief Calculates the bin number corresponding to a d, z, or time value + /// @param value d, z, or time value + /// @param binExtent Bin extent /// @return Bin number - int getBin(float value) const; + static int getBin(float value, float binExtent); /// @brief Finds the maximum density of a DensityMap - /// @param densityMap Map between z bins and corresponding density value - /// @return Iterator of the map entry with the highest density value + /// @param densityMap Map between bins and corresponding density + /// values + /// @return Iterator of the map entry with the highest density DensityMap::const_iterator highestDensityEntry( const DensityMap& densityMap) const; - /// @brief Returns the z position of maximum (surrounding) track density + /// @brief Returns the z and t coordinate of maximum (surrounding) + /// track density + /// @note if time vertex seeding is not enabled, the t coordinate + /// will be set to 0. /// - /// @param densityMap Map from z bins to corresponding track density + /// @param densityMap Map between bins and corresponding density + /// values /// - /// @return The z position of maximum track density - Result getMaxZPosition(DensityMap& densityMap) const; + /// @return The z and t coordinates of maximum track density + Result getMaxZTPosition(DensityMap& densityMap) const; - /// @brief Returns the z position of maximum track density and - /// the estimated width + /// @brief Returns the z-t position of maximum track density + /// and the estimated z-width of the maximum /// - /// @param densityMap Map from z bins to corresponding track density + /// @param densityMap Map between bins and corresponding density + /// values /// - /// @return The z position of maximum track density and width - Result> getMaxZPositionAndWidth( + /// @return The z-t position of the maximum track density and + /// its width + Result getMaxZTPositionAndWidth( DensityMap& densityMap) const; /// @brief Adds a single track to the overall grid density /// - /// @param trk The track to be added. - /// @param mainDensityMap Map from z bins to corresponding track density. + /// @param trk The track to be added + /// @param mainDensityMap Map between bins and corresponding density /// /// @return The density map of the track that was added DensityMap addTrack(const BoundTrackParameters& trk, DensityMap& mainDensityMap) const; - /// @brief Removes a track from the overall grid density + /// @brief Removes a track from the overall grid density. /// - /// @param trackDensityMap Map from z bins to corresponding track density. - /// @note The track density comes from a single track. - /// @param mainDensityMap Map from z bins to corresponding track density. - /// @note The track density comes an arbitrary number of tracks. + /// @param trackDensityMap Map between bins and corresponding density + /// @note The track density comes from a single track + /// @param mainDensityMap Map between bins and corresponding density + /// @note The track density comes from an arbitrary number of tracks void subtractTrack(const DensityMap& trackDensityMap, DensityMap& mainDensityMap) const; private: - /// @brief Function that creates a track density map, i.e., a map of z bins - /// to corresponding density values coming from a single track. - /// - /// @param d0 Transverse impact parameter - /// @param z0 Longitudinal impact parameter - /// @param centralZBin Central z bin of the track (where its density is the highest) - /// @param cov 2x2 impact parameter covariance matrix - DensityMap createTrackGrid(float d0, float z0, int centralZBin, - const Acts::SquareMatrix2& cov) const; - - /// @brief Function that estimates the seed width based on the full width - /// at half maximum (FWHM) of the maximum density peak + /// @brief Function that creates a track density map, i.e., a map from bins + /// to the corresponding density values for a single track. + /// + /// @param impactParams vector containing d0, z0, and t0 of the track + /// @param centralBin Central z and t bin of the track (where its + /// density is the highest) + /// @param cov 3x3 impact parameter covariance matrix + DensityMap createTrackGrid(const Acts::Vector3& impactParams, + const Bin& centralBin, + const Acts::SquareMatrix3& cov) const; + + /// @brief Function that estimates the seed width in z direction based + /// on the full width at half maximum (FWHM) of the maximum density peak /// @note This only works if the maximum is sufficiently isolated since /// overlapping neighboring peaks might lead to an overestimation of the /// seed width. /// - /// @param densityMap Map from z bins to corresponding track density - /// @param maxZ z-position of the maximum density value + /// @param densityMap Map from bins to corresponding track density + /// @param maxZT z-t position of the maximum density value /// /// @return The width Result estimateSeedWidth(const DensityMap& densityMap, - float maxZ) const; + const ZTPosition& maxZT) const; - /// @brief Helper to retrieve values according to a 2-dim normal distribution - /// @note This function is defined in coordinate system centered around d0 and z0 - float normal2D(float d, float z, const SquareMatrix2& cov) const; - - /// @brief Checks the (up to) first three density maxima (only those that have - /// a maximum relative deviation of 'relativeDensityDev' from the main - /// maximum) and take the z-bin of the maximum with the highest surrounding - /// density + /// @brief Helper to retrieve values of an nDim-dimensional normal + /// distribution + /// @note The constant prefactor (2 * pi)^(- nDim / 2) is discarded + /// + /// @param args Coordinates where the Gaussian should be evaluated + /// @note args must be in a coordinate system with origin at the mean + /// values of the Gaussian + /// @param cov Covariance matrix + /// + /// @return Multivariate Gaussian evaluated at args + template + static float multivariateGaussian(const Acts::ActsVector& args, + const Acts::ActsSquareMatrix& cov); + + /// @brief Checks (up to) first three density maxima that have a + /// maximum relative deviation of 'relativeDensityDev' from the + /// global maximum. Returns the bin of the maximum that has the + /// highest surrounding density in z direction. /// - /// @param densityMap Map between z bins and corresponding density value + /// @param densityMap Map between bins and corresponding density values /// - /// @return The z-bin position - int highestDensitySumBin(DensityMap& densityMap) const; + /// @return The bin corresponding to the highest surrounding density + Bin highestDensitySumBin(DensityMap& densityMap) const; - /// @brief Calculates the density sum of a z-bin and its two neighboring bins - /// as needed for 'highestDensitySumBin' + /// @brief Calculates the density sum of a bin and its two neighboring bins + /// in z direction /// - /// @param densityMap Map between z bins and corresponding density value - /// @param zBin The center z-bin whose neighbors we want to sum up + /// @param densityMap Map between bins and corresponding density values + /// @param bin Bin whose neighbors in z we want to sum up /// - /// @return The sum - float getDensitySum(const DensityMap& densityMap, int zBin) const; + /// @return The density sum + float getDensitySum(const DensityMap& densityMap, const Bin& bin) const; Config m_cfg; }; diff --git a/Core/include/Acts/Vertexing/AdaptiveGridTrackDensity.ipp b/Core/include/Acts/Vertexing/AdaptiveGridTrackDensity.ipp index 8f527132213..8f106eb5ac5 100644 --- a/Core/include/Acts/Vertexing/AdaptiveGridTrackDensity.ipp +++ b/Core/include/Acts/Vertexing/AdaptiveGridTrackDensity.ipp @@ -10,20 +10,25 @@ #include -template -float Acts::AdaptiveGridTrackDensity::getBinCenter(int bin) const { - return bin * m_cfg.binSize; +template +float Acts::AdaptiveGridTrackDensity< + spatialTrkGridSize, temporalTrkGridSize>::getBinCenter(int bin, + float binExtent) { + return bin * binExtent; } -template -int Acts::AdaptiveGridTrackDensity::getBin(float value) const { - return static_cast(std::floor(value / m_cfg.binSize - 0.5) + 1); +template +int Acts::AdaptiveGridTrackDensity< + spatialTrkGridSize, temporalTrkGridSize>::getBin(float value, + float binExtent) { + return static_cast(std::floor(value / binExtent - 0.5) + 1); } -template -typename Acts::AdaptiveGridTrackDensity::DensityMap::const_iterator -Acts::AdaptiveGridTrackDensity::highestDensityEntry( - const DensityMap& densityMap) const { +template +typename Acts::AdaptiveGridTrackDensity< + spatialTrkGridSize, temporalTrkGridSize>::DensityMap::const_iterator +Acts::AdaptiveGridTrackDensity:: + highestDensityEntry(const DensityMap& densityMap) const { auto maxEntry = std::max_element( std::begin(densityMap), std::end(densityMap), [](const auto& densityEntry1, const auto& densityEntry2) { @@ -32,118 +37,173 @@ Acts::AdaptiveGridTrackDensity::highestDensityEntry( return maxEntry; } -template -Acts::Result -Acts::AdaptiveGridTrackDensity::getMaxZPosition( - DensityMap& densityMap) const { +template +Acts::Result::ZTPosition> +Acts::AdaptiveGridTrackDensity:: + getMaxZTPosition(DensityMap& densityMap) const { if (densityMap.empty()) { return VertexingError::EmptyInput; } - int zBin = -1; + Bin bin; if (!m_cfg.useHighestSumZPosition) { - zBin = highestDensityEntry(densityMap)->first; + bin = highestDensityEntry(densityMap)->first; } else { // Get z position with highest density sum // of surrounding bins - zBin = highestDensitySumBin(densityMap); + bin = highestDensitySumBin(densityMap); } // Derive corresponding z value - return getBinCenter(zBin); + float maxZ = getBinCenter(bin.first, m_cfg.spatialBinExtent); + + ZTPosition maxValues = std::make_pair(maxZ, 0.); + + // Get t value of the maximum if we do time vertex seeding + if constexpr (temporalTrkGridSize > 1) { + float maxT = getBinCenter(bin.second, m_cfg.temporalBinExtent); + maxValues.second = maxT; + } + + return maxValues; } -template -Acts::Result> -Acts::AdaptiveGridTrackDensity::getMaxZPositionAndWidth( - DensityMap& densityMap) const { - // Get z maximum value - auto maxZRes = getMaxZPosition(densityMap); - if (not maxZRes.ok()) { - return maxZRes.error(); +template +Acts::Result::ZTPositionAndWidth> +Acts::AdaptiveGridTrackDensity:: + getMaxZTPositionAndWidth(DensityMap& densityMap) const { + // Get z value where the density is the highest + auto maxZTRes = getMaxZTPosition(densityMap); + if (not maxZTRes.ok()) { + return maxZTRes.error(); } - float maxZ = *maxZRes; + ZTPosition maxZT = *maxZTRes; // Get seed width estimate - auto widthRes = estimateSeedWidth(densityMap, maxZ); + auto widthRes = estimateSeedWidth(densityMap, maxZT); if (not widthRes.ok()) { return widthRes.error(); } float width = *widthRes; - std::pair returnPair{maxZ, width}; - return returnPair; + ZTPositionAndWidth maxZTAndWidth{maxZT, width}; + return maxZTAndWidth; } -template -typename Acts::AdaptiveGridTrackDensity::DensityMap -Acts::AdaptiveGridTrackDensity::addTrack( - const Acts::BoundTrackParameters& trk, DensityMap& mainDensityMap) const { - SquareMatrix2 cov = trk.spatialImpactParameterCovariance().value(); - float d0 = trk.parameters()[0]; - float z0 = trk.parameters()[1]; +template +typename Acts::AdaptiveGridTrackDensity::DensityMap +Acts::AdaptiveGridTrackDensity:: + addTrack(const Acts::BoundTrackParameters& trk, + DensityMap& mainDensityMap) const { + ActsVector<3> impactParams = trk.impactParameters(); + ActsSquareMatrix<3> cov = trk.impactParameterCovariance().value(); // Calculate bin in d direction - int centralDBin = getBin(d0); + int centralDBin = getBin(impactParams(0), m_cfg.spatialBinExtent); // Check if current track affects grid density - if (std::abs(centralDBin) > (trkGridSize - 1) / 2.) { + if (std::abs(centralDBin) > (spatialTrkGridSize - 1) / 2.) { DensityMap emptyTrackDensityMap; return emptyTrackDensityMap; } // Calculate bin in z direction - int centralZBin = getBin(z0); + int centralZBin = getBin(impactParams(1), m_cfg.spatialBinExtent); + + // If we don't do time vertex seeding, the time index is set to 0 + Bin centralBin = std::make_pair(centralZBin, 0.); + + // Calculate bin in t direction if we do time vertex seeding + if constexpr (temporalTrkGridSize > 1) { + int centralTBin = getBin(impactParams(2), m_cfg.temporalBinExtent); + centralBin.second = centralTBin; + } - DensityMap trackDensityMap = createTrackGrid(d0, z0, centralZBin, cov); + DensityMap trackDensityMap = createTrackGrid(impactParams, centralBin, cov); for (const auto& densityEntry : trackDensityMap) { - int zBin = densityEntry.first; + Bin bin = densityEntry.first; float trackDensity = densityEntry.second; // Check if z bin is already part of the main grid - if (mainDensityMap.count(zBin) == 1) { - mainDensityMap.at(zBin) += trackDensity; + if (mainDensityMap.count(bin) == 1) { + mainDensityMap.at(bin) += trackDensity; } else { - mainDensityMap[zBin] = trackDensity; + mainDensityMap[bin] = trackDensity; } } return trackDensityMap; } -template -void Acts::AdaptiveGridTrackDensity::subtractTrack( - const DensityMap& trackDensityMap, DensityMap& mainDensityMap) const { +template +void Acts::AdaptiveGridTrackDensity:: + subtractTrack(const DensityMap& trackDensityMap, + DensityMap& mainDensityMap) const { for (auto it = trackDensityMap.begin(); it != trackDensityMap.end(); it++) { mainDensityMap.at(it->first) -= it->second; } } -template -typename Acts::AdaptiveGridTrackDensity::DensityMap -Acts::AdaptiveGridTrackDensity::createTrackGrid( - float d0, float z0, int centralZBin, const Acts::SquareMatrix2& cov) const { +template +typename Acts::AdaptiveGridTrackDensity::DensityMap +Acts::AdaptiveGridTrackDensity:: + createTrackGrid(const Acts::Vector3& impactParams, const Bin& centralBin, + const Acts::SquareMatrix3& cov) const { DensityMap trackDensityMap; - int halfTrkGridSize = (trkGridSize - 1) / 2; - int firstZBin = centralZBin - halfTrkGridSize; - // Loop over columns - for (int j = 0; j < trkGridSize; j++) { - int zBin = firstZBin + j; - float z = getBinCenter(zBin); - trackDensityMap[zBin] = normal2D(-d0, z - z0, cov); + int halfSpatialTrkGridSize = (spatialTrkGridSize - 1) / 2; + int firstZBin = centralBin.first - halfSpatialTrkGridSize; + + // If we don't do time vertex seeding, firstTBin will be 0. + int halfTemporalTrkGridSize = (temporalTrkGridSize - 1) / 2; + int firstTBin = centralBin.second - halfTemporalTrkGridSize; + + // Loop over bins + for (int i = 0; i < temporalTrkGridSize; i++) { + int tBin = firstTBin + i; + // If we don't do vertex time seeding, we set the time to 0 since it will be + // discarded in the for loop below anyways + float t = 0; + if constexpr (temporalTrkGridSize > 1) { + t = getBinCenter(tBin, m_cfg.temporalBinExtent); + } + for (int j = 0; j < spatialTrkGridSize; j++) { + int zBin = firstZBin + j; + float z = getBinCenter(zBin, m_cfg.spatialBinExtent); + // Bin coordinates in the d-z-t plane + Acts::Vector3 binCoords(0., z, t); + // Transformation to coordinate system with origin at the track center + binCoords -= impactParams; + Bin bin = std::make_pair(zBin, tBin); + if constexpr (temporalTrkGridSize == 1) { + trackDensityMap[bin] = multivariateGaussian<2>( + binCoords.head<2>(), cov.topLeftCorner<2, 2>()); + } else { + trackDensityMap[bin] = multivariateGaussian<3>(binCoords, cov); + } + } } return trackDensityMap; } -template -Acts::Result -Acts::AdaptiveGridTrackDensity::estimateSeedWidth( - const DensityMap& densityMap, float maxZ) const { +template +Acts::Result Acts::AdaptiveGridTrackDensity< + spatialTrkGridSize, + temporalTrkGridSize>::estimateSeedWidth(const DensityMap& densityMap, + const ZTPosition& maxZT) const { if (densityMap.empty()) { return VertexingError::EmptyInput; } - // Get z bin of max density and max density value - int zMaxBin = getBin(maxZ); - const float maxValue = densityMap.at(zMaxBin); + // Get z bin of max density + int zMaxBin = getBin(maxZT.first, m_cfg.spatialBinExtent); + int tMaxBin = 0; + // Fill the time bin with a non-zero value if we do time vertex seeding + if constexpr (temporalTrkGridSize > 1) { + tMaxBin = getBin(maxZT.second, m_cfg.temporalBinExtent); + } + const float maxValue = densityMap.at(std::make_pair(zMaxBin, tMaxBin)); int rhmBin = zMaxBin; float gridValue = maxValue; @@ -152,21 +212,21 @@ Acts::AdaptiveGridTrackDensity::estimateSeedWidth( bool binFilled = true; while (gridValue > maxValue / 2) { // Check if we are still operating on continuous z values - if (densityMap.count(rhmBin + 1) == 0) { + if (densityMap.count(std::make_pair(rhmBin + 1, tMaxBin)) == 0) { binFilled = false; break; } rhmBin += 1; - gridValue = densityMap.at(rhmBin); + gridValue = densityMap.at(std::make_pair(rhmBin, tMaxBin)); } // Use linear approximation to find better z value for FWHM between bins float rightDensity = 0; if (binFilled) { - rightDensity = densityMap.at(rhmBin); + rightDensity = densityMap.at(std::make_pair(rhmBin, tMaxBin)); } - float leftDensity = densityMap.at(rhmBin - 1); - float deltaZ1 = m_cfg.binSize * (maxValue / 2 - leftDensity) / + float leftDensity = densityMap.at(std::make_pair(rhmBin - 1, tMaxBin)); + float deltaZ1 = m_cfg.spatialBinExtent * (maxValue / 2 - leftDensity) / (rightDensity - leftDensity); int lhmBin = zMaxBin; @@ -174,25 +234,25 @@ Acts::AdaptiveGridTrackDensity::estimateSeedWidth( binFilled = true; while (gridValue > maxValue / 2) { // Check if we are still operating on continuous z values - if (densityMap.count(lhmBin - 1) == 0) { + if (densityMap.count(std::make_pair(lhmBin - 1, tMaxBin)) == 0) { binFilled = false; break; } lhmBin -= 1; - gridValue = densityMap.at(lhmBin); + gridValue = densityMap.at(std::make_pair(lhmBin, tMaxBin)); } // Use linear approximation to find better z value for FWHM between bins - rightDensity = densityMap.at(lhmBin + 1); + rightDensity = densityMap.at(std::make_pair(lhmBin + 1, tMaxBin)); if (binFilled) { - leftDensity = densityMap.at(lhmBin); + leftDensity = densityMap.at(std::make_pair(lhmBin, tMaxBin)); } else { leftDensity = 0; } - float deltaZ2 = m_cfg.binSize * (rightDensity - maxValue / 2) / + float deltaZ2 = m_cfg.spatialBinExtent * (rightDensity - maxValue / 2) / (rightDensity - leftDensity); - float fwhm = (rhmBin - lhmBin) * m_cfg.binSize - deltaZ1 - deltaZ2; + float fwhm = (rhmBin - lhmBin) * m_cfg.spatialBinExtent - deltaZ1 - deltaZ2; // FWHM = 2.355 * sigma float width = fwhm / 2.355f; @@ -200,76 +260,89 @@ Acts::AdaptiveGridTrackDensity::estimateSeedWidth( return std::isnormal(width) ? width : 0.0f; } -template -float Acts::AdaptiveGridTrackDensity::normal2D( - float d, float z, const Acts::SquareMatrix2& cov) const { - float det = cov.determinant(); - float coef = 1 / std::sqrt(det); - float expo = - -1 / (2 * det) * - (cov(1, 1) * d * d - (cov(0, 1) + cov(1, 0)) * d * z + cov(0, 0) * z * z); - return coef * safeExp(expo); +template +template +float Acts::AdaptiveGridTrackDensity:: + multivariateGaussian(const Acts::ActsVector& args, + const Acts::ActsSquareMatrix& cov) { + float expo = -0.5 * args.transpose().dot(cov.inverse() * args); + float gaussianDensity = safeExp(expo) / std::sqrt(cov.determinant()); + return gaussianDensity; } -template -int Acts::AdaptiveGridTrackDensity::highestDensitySumBin( - DensityMap& densityMap) const { - // Checks the first (up to) 3 density maxima, if they are close, checks which - // one has the highest surrounding density sum (the two neighboring bins) - +template +typename Acts::AdaptiveGridTrackDensity::Bin +Acts::AdaptiveGridTrackDensity:: + highestDensitySumBin(DensityMap& densityMap) const { // The global maximum auto firstMax = highestDensityEntry(densityMap); - int zBinFirstMax = firstMax->first; + Bin binFirstMax = firstMax->first; float valueFirstMax = firstMax->second; - float firstSum = getDensitySum(densityMap, zBinFirstMax); + float firstSum = getDensitySum(densityMap, binFirstMax); + // Smaller maxima must have a density of at least: + // valueFirstMax - densityDeviation float densityDeviation = valueFirstMax * m_cfg.maxRelativeDensityDev; // Get the second highest maximum - densityMap.at(zBinFirstMax) = 0; + densityMap.at(binFirstMax) = 0; auto secondMax = highestDensityEntry(densityMap); - int zBinSecondMax = secondMax->first; + Bin binSecondMax = secondMax->first; float valueSecondMax = secondMax->second; float secondSum = 0; if (valueFirstMax - valueSecondMax < densityDeviation) { - secondSum = getDensitySum(densityMap, zBinSecondMax); + secondSum = getDensitySum(densityMap, binSecondMax); + } else { + // If the second maximum is not sufficiently large the third maximum won't + // be either + densityMap.at(binFirstMax) = valueFirstMax; + return binFirstMax; } // Get the third highest maximum - densityMap.at(zBinSecondMax) = 0; + densityMap.at(binSecondMax) = 0; auto thirdMax = highestDensityEntry(densityMap); - int zBinThirdMax = thirdMax->first; + Bin binThirdMax = thirdMax->first; float valueThirdMax = thirdMax->second; float thirdSum = 0; if (valueFirstMax - valueThirdMax < densityDeviation) { - thirdSum = getDensitySum(densityMap, zBinThirdMax); + thirdSum = getDensitySum(densityMap, binThirdMax); } // Revert back to original values - densityMap.at(zBinFirstMax) = valueFirstMax; - densityMap.at(zBinSecondMax) = valueSecondMax; + densityMap.at(binFirstMax) = valueFirstMax; + densityMap.at(binSecondMax) = valueSecondMax; - // Return the z-bin position of the highest density sum + // Return the z bin position of the highest density sum if (secondSum > firstSum && secondSum > thirdSum) { - return zBinSecondMax; + return binSecondMax; } if (thirdSum > secondSum && thirdSum > firstSum) { - return zBinThirdMax; + return binThirdMax; } - return zBinFirstMax; + return binFirstMax; } -template -float Acts::AdaptiveGridTrackDensity::getDensitySum( - const DensityMap& densityMap, int zBin) const { - float sum = densityMap.at(zBin); +template +float Acts::AdaptiveGridTrackDensity:: + getDensitySum(const DensityMap& densityMap, const Bin& bin) const { + // Add density from the bin. + float sum = densityMap.at(bin); // Check if neighboring bins are part of the densityMap and add them (if they - // are not part of the map, we assume them to be 0) Note that each key in a - // map is unique; the .count() function therefore returns either 0 or 1 - if (densityMap.count(zBin - 1) == 1) { - sum += densityMap.at(zBin - 1); + // are not part of the map, we assume them to be 0). + // Note that each key in a map is unique; the .count() function therefore + // returns either 0 or 1. + Bin binShifted = bin; + // Add density from the neighboring bin in -z direction. + binShifted.first -= 1; + if (densityMap.count(binShifted) == 1) { + sum += densityMap.at(binShifted); } - if (densityMap.count(zBin + 1) == 1) { - sum += densityMap.at(zBin + 1); + + // Add density from the neighboring bin in +z direction. + binShifted.first += 2; + if (densityMap.count(binShifted) == 1) { + sum += densityMap.at(binShifted); } return sum; } diff --git a/Core/include/Acts/Vertexing/ImpactPointEstimator.hpp b/Core/include/Acts/Vertexing/ImpactPointEstimator.hpp index f12a3054e91..30d385859ad 100644 --- a/Core/include/Acts/Vertexing/ImpactPointEstimator.hpp +++ b/Core/include/Acts/Vertexing/ImpactPointEstimator.hpp @@ -40,7 +40,7 @@ struct ImpactParametersAndSigma { /// /// @brief Estimator for impact point calculations /// A description of the underlying mathematics can be found here: -/// https://github.com/acts-project/acts/pull/2414 +/// https://github.com/acts-project/acts/pull/2506 /// TODO: Upload reference at a better place template > @@ -81,10 +81,6 @@ class ImpactPointEstimator { int maxIterations = 20; /// Desired precision of deltaPhi in Newton method double precision = 1.e-10; - /// Minimum q/p value - double minQoP = 1e-15; - /// Maximum curvature value - double maxRho = 1e+15; }; /// @brief Constructor @@ -147,6 +143,8 @@ class ImpactPointEstimator { /// corresponding 3D PCA. Returns also the momentum direction at the 3D PCA. /// The template parameter nDim determines whether we calculate the 3D /// distance (nDim = 3) or the 4D distance (nDim = 4) to the 3D PCA. + /// @note For straight tracks we use an analytical solution; for helical + /// tracks we use the Newton method. /// /// @tparam nDim Number of dimensions used to compute compatibility /// @note If nDim = 3 we only consider spatial dimensions; if nDim = 4, we diff --git a/Core/include/Acts/Vertexing/ImpactPointEstimator.ipp b/Core/include/Acts/Vertexing/ImpactPointEstimator.ipp index 6e023f384d9..a5a698e0295 100644 --- a/Core/include/Acts/Vertexing/ImpactPointEstimator.ipp +++ b/Core/include/Acts/Vertexing/ImpactPointEstimator.ipp @@ -241,13 +241,65 @@ Acts::ImpactPointEstimator:: // Reference point R Vector3 refPoint = trkParams.referenceSurface().center(gctx); - // Perigee parameters (parameters of 2D PCA) + // Extract charge-related particle parameters + double absoluteCharge = trkParams.particleHypothesis().absoluteCharge(); + double qOvP = trkParams.parameters()[BoundIndices::eBoundQOverP]; + + // Z-component of the B field at the reference position. + // Note that we assume a constant B field here! + auto fieldRes = m_cfg.bField->getField(refPoint, state.fieldCache); + if (!fieldRes.ok()) { + return fieldRes.error(); + } + double bZ = (*fieldRes)[eZ]; + + // The particle moves on a straight trajectory if its charge is 0 or if there + // is no B field. In that case, the 3D PCA can be calculated analytically, see + // Sec 3.2 of the reference. + if (absoluteCharge == 0. or bZ == 0.) { + // Momentum direction (constant for straight tracks) + Vector3 momDirStraightTrack = trkParams.direction(); + + // Current position on the track + Vector3 positionOnTrack = trkParams.position(gctx); + + // Distance between positionOnTrack and the 3D PCA + ActsScalar distanceToPca = + (vtxPos.template head<3>() - positionOnTrack).dot(momDirStraightTrack); + + // 3D PCA + ActsVector pcaStraightTrack; + pcaStraightTrack.template head<3>() = + positionOnTrack + distanceToPca * momDirStraightTrack; + if constexpr (nDim == 4) { + // Track time at positionOnTrack + double timeOnTrack = trkParams.parameters()[BoundIndices::eBoundTime]; + + ActsScalar m0 = trkParams.particleHypothesis().mass(); + ActsScalar p = std::abs(absoluteCharge / qOvP); + + // Speed in units of c + ActsScalar beta = p / std::hypot(p, m0); + + pcaStraightTrack[3] = timeOnTrack + distanceToPca / beta; + } + + // Vector pointing from the vertex position to the 3D PCA + ActsVector deltaRStraightTrack = pcaStraightTrack - vtxPos; + + return std::make_pair(deltaRStraightTrack, momDirStraightTrack); + } + + // Charged particles in a constant B field follow a helical trajectory. In + // that case, we calculate the 3D PCA using the Newton method, see Sec 4.2 in + // the reference. + + // Spatial Perigee parameters (i.e., spatial parameters of 2D PCA) double d0 = trkParams.parameters()[BoundIndices::eBoundLoc0]; double z0 = trkParams.parameters()[BoundIndices::eBoundLoc1]; + // Momentum angles at 2D PCA double phiP = trkParams.parameters()[BoundIndices::eBoundPhi]; double theta = trkParams.parameters()[BoundIndices::eBoundTheta]; - double qOvP = trkParams.parameters()[BoundIndices::eBoundQOverP]; - double tP = trkParams.parameters()[BoundIndices::eBoundTime]; // Functions of the polar angle theta for later use double sinTheta = std::sin(theta); double cotTheta = 1. / std::tan(theta); @@ -256,22 +308,8 @@ Acts::ImpactPointEstimator:: // Note that phi corresponds to phiV in the reference. double phi = phiP; - // B-field z-component at the reference position. - // Note that we assume a constant B field here! - auto fieldRes = m_cfg.bField->getField(refPoint, state.fieldCache); - if (!fieldRes.ok()) { - return fieldRes.error(); - } - double bZ = (*fieldRes)[eZ]; - // Signed radius of the helix on which the particle moves - double rho = 0; - // Curvature is infinite w/o b field - if (bZ == 0. || std::abs(qOvP) < m_cfg.minQoP) { - rho = m_cfg.maxRho; - } else { - rho = sinTheta * (1. / qOvP) / bZ; - } + double rho = sinTheta * (1. / qOvP) / bZ; // Position of the helix center. // We can set the z-position to a convenient value since it is not fixed by @@ -308,9 +346,11 @@ Acts::ImpactPointEstimator:: helixCenter + rho * Vector3(-sinPhi, cosPhi, -cotTheta * phi); if constexpr (nDim == 4) { + // Time at the 2D PCA P + double tP = trkParams.parameters()[BoundIndices::eBoundTime]; + ActsScalar m0 = trkParams.particleHypothesis().mass(); - ActsScalar p = - std::abs(trkParams.particleHypothesis().absoluteCharge() / qOvP); + ActsScalar p = std::abs(absoluteCharge / qOvP); // Speed in units of c ActsScalar beta = p / std::hypot(p, m0); diff --git a/Core/src/Geometry/CylinderVolumeBuilder.cpp b/Core/src/Geometry/CylinderVolumeBuilder.cpp index ceb72c7affb..a4520d3b828 100644 --- a/Core/src/Geometry/CylinderVolumeBuilder.cpp +++ b/Core/src/Geometry/CylinderVolumeBuilder.cpp @@ -161,6 +161,15 @@ Acts::CylinderVolumeBuilder::trackingVolume( wConfig.cVolumeConfig = analyzeContent(gctx, centralLayers, centralVolumes); wConfig.pVolumeConfig = analyzeContent(gctx, positiveLayers, {}); // TODO + bool hasLayers = wConfig.nVolumeConfig.present || + wConfig.cVolumeConfig.present || + wConfig.pVolumeConfig.present; + + if (!hasLayers) { + ACTS_INFO("No layers present, returning nullptr"); + return nullptr; + } + std::string layerConfiguration = "|"; if (wConfig.nVolumeConfig) { // negative layers are present diff --git a/Core/src/Geometry/TrackingGeometryBuilder.cpp b/Core/src/Geometry/TrackingGeometryBuilder.cpp index 911909ad1b6..90a10b0ba04 100644 --- a/Core/src/Geometry/TrackingGeometryBuilder.cpp +++ b/Core/src/Geometry/TrackingGeometryBuilder.cpp @@ -47,7 +47,14 @@ Acts::TrackingGeometryBuilder::trackingGeometry( for (auto& volumeBuilder : m_cfg.trackingVolumeBuilders) { // assign a new highest volume (and potentially wrap around the given // highest volume so far) - highestVolume = volumeBuilder(gctx, highestVolume, nullptr); + auto volume = volumeBuilder(gctx, highestVolume, nullptr); + if (!volume) { + ACTS_INFO( + "Received nullptr volume from builder, keeping previous highest " + "volume"); + } else { + highestVolume = std::move(volume); + } } // create the TrackingGeometry & decorate it with the material diff --git a/Examples/Python/src/EventData.cpp b/Examples/Python/src/EventData.cpp index ffbcaba41be..455b6c541cc 100644 --- a/Examples/Python/src/EventData.cpp +++ b/Examples/Python/src/EventData.cpp @@ -36,9 +36,20 @@ void addEventData(Context& ctx) { [](py::object /* self */) { return Acts::ParticleHypothesis::pion(); }) - .def_property_readonly_static("electron", [](py::object /* self */) { - return Acts::ParticleHypothesis::electron(); - }); + .def_property_readonly_static( + "electron", + [](py::object /* self */) { + return Acts::ParticleHypothesis::electron(); + }) + .def_property_readonly_static( + "geantino", + [](py::object /* self */) { + return Acts::ParticleHypothesis::geantino(); + }) + .def_property_readonly_static( + "chargedGeantino", [](py::object /* self */) { + return Acts::ParticleHypothesis::chargedGeantino(); + }); } } // namespace Acts::Python diff --git a/Examples/Python/src/Propagation.cpp b/Examples/Python/src/Propagation.cpp index 15f04f74d6e..9876481d553 100644 --- a/Examples/Python/src/Propagation.cpp +++ b/Examples/Python/src/Propagation.cpp @@ -116,9 +116,9 @@ void addPropagation(Context& ctx) { propagatorImpl, randomNumberSvc, mode, sterileLogger, debugOutput, energyLoss, multipleScattering, recordMaterialInteractions, ntests, d0Sigma, z0Sigma, phiSigma, thetaSigma, qpSigma, tSigma, phiRange, - etaRange, ptRange, ptLoopers, maxStepSize, propagationStepCollection, - propagationMaterialCollection, covarianceTransport, covariances, - correlations); + etaRange, ptRange, particleHypothesis, ptLoopers, maxStepSize, + propagationStepCollection, propagationMaterialCollection, + covarianceTransport, covariances, correlations); py::class_>( diff --git a/Tests/UnitTests/Core/Vertexing/AdaptiveGridTrackDensityTests.cpp b/Tests/UnitTests/Core/Vertexing/AdaptiveGridTrackDensityTests.cpp index 2b33c5d789c..87f2e141070 100644 --- a/Tests/UnitTests/Core/Vertexing/AdaptiveGridTrackDensityTests.cpp +++ b/Tests/UnitTests/Core/Vertexing/AdaptiveGridTrackDensityTests.cpp @@ -37,21 +37,30 @@ namespace Test { using Covariance = BoundSquareMatrix; +Covariance makeRandomCovariance(int seed = 31415) { + std::srand(seed); + Covariance randMat((Covariance::Random() + 1.5 * Covariance::Identity()) * + 0.05); + + // symmetric covariance matrix + Covariance covMat = 0.5 * (randMat + randMat.transpose()); + + return covMat; +} + BOOST_AUTO_TEST_CASE(compare_to_analytical_solution_for_single_track) { using Vector2 = Eigen::Matrix; using Matrix2 = Eigen::Matrix; // Using a large track grid so we can choose a small bin size - const int trkGridSize = 4001; + const int spatialTrkGridSize = 4001; // Arbitrary (but small) bin size - const float binSize = 3.1e-4; + const float binExtent = 3.1e-4; // Arbitrary impact parameters const float d0 = 0.4; const float z0 = -0.2; Vector2 impactParameters{d0, z0}; - Covariance covMat(Covariance::Identity() * 0.05); - covMat(0, 1) = -0.02; - covMat(1, 0) = -0.02; + Covariance covMat = makeRandomCovariance(); Matrix2 subCovMat = covMat.block<2, 2>(0, 0).cast(); BoundVector paramVec; paramVec << d0, z0, 0, 0, 0, 0; @@ -63,11 +72,11 @@ BOOST_AUTO_TEST_CASE(compare_to_analytical_solution_for_single_track) { BoundTrackParameters params1(perigeeSurface, paramVec, covMat, ParticleHypothesis::pion()); - AdaptiveGridTrackDensity::Config cfg(binSize); - AdaptiveGridTrackDensity grid(cfg); + AdaptiveGridTrackDensity::Config cfg(binExtent); + AdaptiveGridTrackDensity grid(cfg); // Empty map - AdaptiveGridTrackDensity::DensityMap mainDensityMap; + AdaptiveGridTrackDensity::DensityMap mainDensityMap; // Add track auto trackDensityMap = grid.addTrack(params1, mainDensityMap); @@ -85,10 +94,10 @@ BOOST_AUTO_TEST_CASE(compare_to_analytical_solution_for_single_track) { for (auto const& it : mainDensityMap) { // Extract variables for better readability - int zBin = it.first; + int zBin = it.first.first; float density = it.second; // Argument for 2D gaussian - Vector2 dzVec{0., grid.getBinCenter(zBin)}; + Vector2 dzVec{0., grid.getBinCenter(zBin, binExtent)}; // Compute correct density... float correctDensity = gaussian2D(dzVec, impactParameters, subCovMat); // ... and check if our result is equivalent @@ -106,11 +115,11 @@ BOOST_AUTO_TEST_CASE(compare_to_analytical_solution_for_single_track) { subCovMat.determinant() / subCovMat(0, 0)); // Estimate maximum z position and seed width - auto res = grid.getMaxZPositionAndWidth(mainDensityMap); + auto res = grid.getMaxZTPositionAndWidth(mainDensityMap); BOOST_CHECK(res.ok()); // Extract variables for better readability... - float maxZ = res.value().first; + float maxZ = res.value().first.first; float fwhm = res.value().second * 2.355f; // ... and check if they are correct (note: the optimization is not as exact @@ -120,15 +129,116 @@ BOOST_AUTO_TEST_CASE(compare_to_analytical_solution_for_single_track) { CHECK_CLOSE_REL(fwhm, correctFWHM, relTolOptimization); } -BOOST_AUTO_TEST_CASE(check_seed_width_estimation) { +BOOST_AUTO_TEST_CASE( + compare_to_analytical_solution_for_single_track_with_time) { + // Number of bins in z- and t-direction + const int spatialTrkGridSize = 401; + const int temporalTrkGridSize = 401; + // Bin extents + const float spatialBinExtent = 3.1e-3; + const float temporalBinExtent = 3.1e-3; + // Arbitrary impact parameters + const float d0 = -0.1; + const float z0 = -0.2; + const float t0 = 0.1; + Vector3 impactParameters{d0, z0, t0}; + + // symmetric covariance matrix + Covariance covMat = makeRandomCovariance(); + + BoundVector paramVec; + paramVec << d0, z0, 0, 0, 0, t0; + // Create perigee surface + std::shared_ptr perigeeSurface = + Surface::makeShared(Vector3(0., 0., 0.)); + + BoundTrackParameters params(perigeeSurface, paramVec, covMat, + ParticleHypothesis::pion()); + + ActsSquareMatrix<3> ipCov = params.impactParameterCovariance().value(); + + AdaptiveGridTrackDensity::Config cfg( + spatialBinExtent, temporalBinExtent); + AdaptiveGridTrackDensity grid(cfg); + + // Empty map + AdaptiveGridTrackDensity::DensityMap + mainDensityMap; + + // Add track + auto trackDensityMap = grid.addTrack(params, mainDensityMap); + + float relTol = 1e-5; + float small = 1e-5; + + auto gaussian3D = [&](const Vector3& args, const Vector3& mus, + const SquareMatrix3& sigmas) { + Vector3 diffs = args - mus; + float coef = 1 / std::sqrt(sigmas.determinant()); + float expo = -0.5 * diffs.transpose().dot(sigmas.inverse() * diffs); + return coef * std::exp(expo); + }; + + for (auto const& it : mainDensityMap) { + // Extract variables for better readability + float z = grid.getBinCenter(it.first.first, spatialBinExtent); + float t = grid.getBinCenter(it.first.second, temporalBinExtent); + float density = it.second; + // Argument for 3D gaussian + Vector3 dztVec{0., z, t}; + + // Compute correct density... + float correctDensity = gaussian3D(dztVec, impactParameters, ipCov); + + // ... and check if our result is equivalent + CHECK_CLOSE_OR_SMALL(density, correctDensity, relTol, small); + } + + // The analytical calculations of the following can be found here: + // https://github.com/acts-project/acts/pull/2460. + // TODO: upload reference at a better place. + // Analytical maximum of the Gaussian + ActsSquareMatrix<3> ipWeights = ipCov.inverse(); + ActsScalar denom = + ipWeights(1, 1) * ipWeights(2, 2) - ipWeights(1, 2) * ipWeights(1, 2); + + ActsScalar zNom = + ipWeights(0, 1) * ipWeights(2, 2) - ipWeights(0, 2) * ipWeights(1, 2); + ActsScalar correctMaxZ = zNom / denom * d0 + z0; + + ActsScalar tNom = + ipWeights(0, 2) * ipWeights(1, 1) - ipWeights(0, 1) * ipWeights(1, 2); + ActsScalar correctMaxT = tNom / denom * d0 + t0; + + // Analytical FWHM of the Gaussian + ActsScalar correctFWHM = 2. * std::sqrt(2 * std::log(2.) / ipWeights(1, 1)); + + // Estimate maximum z position and seed width + auto res = grid.getMaxZTPositionAndWidth(mainDensityMap); + BOOST_CHECK(res.ok()); + + // Extract variables for better readability... + float maxZ = res.value().first.first; + float maxT = res.value().first.second; + float fwhm = res.value().second * 2.355f; + + // ... and check if they are correct (note: the optimization is not as exact + // as the density values). + float relTolOptimization = 1e-1; + CHECK_CLOSE_REL(maxZ, correctMaxZ, relTolOptimization); + CHECK_CLOSE_REL(maxT, correctMaxT, relTolOptimization); + CHECK_CLOSE_REL(fwhm, correctFWHM, relTolOptimization); +} + +BOOST_AUTO_TEST_CASE(seed_width_estimation) { // Dummy track grid size (not needed for this unit test) - const int trkGridSize = 1; - float binSize = 2.; - AdaptiveGridTrackDensity::Config cfg(binSize); - AdaptiveGridTrackDensity grid(cfg); + const int spatialTrkGridSize = 1; + float binExtent = 2.; + AdaptiveGridTrackDensity::Config cfg(binExtent); + AdaptiveGridTrackDensity grid(cfg); // Empty map - AdaptiveGridTrackDensity::DensityMap mainDensityMap; + AdaptiveGridTrackDensity::DensityMap mainDensityMap; // z-position of the maximum density float correctMaxZ = -2.; @@ -138,16 +248,16 @@ BOOST_AUTO_TEST_CASE(check_seed_width_estimation) { // of 20 mm. The linear approximation we use during the seed width estimation // should be exact in this case. for (int i = -6; i <= 4; i++) { - mainDensityMap[i] = - 1.0 - 0.1 * std::abs(correctMaxZ - grid.getBinCenter(i)); + mainDensityMap[std::make_pair(i, 0)] = + 1.0 - 0.1 * std::abs(correctMaxZ - grid.getBinCenter(i, binExtent)); } // Get maximum z position and corresponding seed width - auto res = grid.getMaxZPositionAndWidth(mainDensityMap); + auto res = grid.getMaxZTPositionAndWidth(mainDensityMap); BOOST_CHECK(res.ok()); // Check if we found the correct maximum - float maxZ = res.value().first; + float maxZ = res.value().first.first; BOOST_CHECK_EQUAL(maxZ, correctMaxZ); // Calculate full width at half maximum from seed width and check if it's @@ -156,13 +266,13 @@ BOOST_AUTO_TEST_CASE(check_seed_width_estimation) { BOOST_CHECK_EQUAL(fwhm, 10.); } -BOOST_AUTO_TEST_CASE(adaptive_gaussian_grid_density_track_adding_test) { - const int trkGridSize = 15; +BOOST_AUTO_TEST_CASE(track_adding) { + const int spatialTrkGridSize = 15; - double binSize = 0.1; // mm + double binExtent = 0.1; // mm - AdaptiveGridTrackDensity::Config cfg(binSize); - AdaptiveGridTrackDensity grid(cfg); + AdaptiveGridTrackDensity::Config cfg(binExtent); + AdaptiveGridTrackDensity grid(cfg); // Create some test tracks in such a way that some tracks // e.g. overlap and that certain tracks need to be inserted @@ -192,7 +302,7 @@ BOOST_AUTO_TEST_CASE(adaptive_gaussian_grid_density_track_adding_test) { ParticleHypothesis::pion()); // Empty map - AdaptiveGridTrackDensity::DensityMap mainDensityMap; + AdaptiveGridTrackDensity::DensityMap mainDensityMap; // Track is too far away from z axis and was not added auto trackDensityMap = grid.addTrack(params0, mainDensityMap); @@ -200,38 +310,49 @@ BOOST_AUTO_TEST_CASE(adaptive_gaussian_grid_density_track_adding_test) { // Track should have been entirely added to both grids trackDensityMap = grid.addTrack(params1, mainDensityMap); - BOOST_CHECK_EQUAL(mainDensityMap.size(), trkGridSize); + BOOST_CHECK_EQUAL(mainDensityMap.size(), spatialTrkGridSize); // Track should have been entirely added to both grids trackDensityMap = grid.addTrack(params2, mainDensityMap); - BOOST_CHECK_EQUAL(mainDensityMap.size(), 2 * trkGridSize); + BOOST_CHECK_EQUAL(mainDensityMap.size(), 2 * spatialTrkGridSize); // Track 3 has overlap of 2 bins with track 1 trackDensityMap = grid.addTrack(params3, mainDensityMap); - BOOST_CHECK_EQUAL(mainDensityMap.size(), 3 * trkGridSize - 2); + BOOST_CHECK_EQUAL(mainDensityMap.size(), 3 * spatialTrkGridSize - 2); // Add first track again, should *not* introduce new z entries trackDensityMap = grid.addTrack(params1, mainDensityMap); - BOOST_CHECK_EQUAL(mainDensityMap.size(), 3 * trkGridSize - 2); + BOOST_CHECK_EQUAL(mainDensityMap.size(), 3 * spatialTrkGridSize - 2); } -BOOST_AUTO_TEST_CASE(adaptive_gaussian_grid_density_max_z_and_width_test) { - const int trkGridSize = 29; +BOOST_AUTO_TEST_CASE(max_z_t_and_width) { + const int spatialTrkGridSize = 29; + const int temporalTrkGridSize = 29; + + // spatial and temporal bin extent + double binExtent = 0.05; - double binSize = 0.05; // mm + // 1D grid of z values + AdaptiveGridTrackDensity::Config cfg1D(binExtent); + AdaptiveGridTrackDensity grid1D(cfg1D); - AdaptiveGridTrackDensity::Config cfg(binSize); - AdaptiveGridTrackDensity grid(cfg); + // 2D grid of z and t values + AdaptiveGridTrackDensity::Config + cfg2D(binExtent, binExtent); + AdaptiveGridTrackDensity grid2D( + cfg2D); // Create some test tracks Covariance covMat(Covariance::Identity() * 0.005); float z0Trk1 = 0.25; + float t0Trk1 = 0.05; float z0Trk2 = -10.95; + float t0Trk2 = 0.1; BoundVector paramVec1; - paramVec1 << 0.02, z0Trk1, 0, 0, 0, 0; + paramVec1 << 0.02, z0Trk1, 0, 0, 0, t0Trk1; BoundVector paramVec2; - paramVec2 << 0.01, z0Trk2, 0, 0, 0, 0; + paramVec2 << 0.01, z0Trk2, 0, 0, 0, t0Trk2; // Create perigee surface std::shared_ptr perigeeSurface = @@ -242,37 +363,65 @@ BOOST_AUTO_TEST_CASE(adaptive_gaussian_grid_density_max_z_and_width_test) { BoundTrackParameters params2(perigeeSurface, paramVec2, covMat, ParticleHypothesis::pion()); - // Empty map - AdaptiveGridTrackDensity::DensityMap mainDensityMap; - - auto trackDensityMap = grid.addTrack(params1, mainDensityMap); - auto res1 = grid.getMaxZPosition(mainDensityMap); - BOOST_CHECK(res1.ok()); - // Maximum should be at z0Trk1 position - BOOST_CHECK_EQUAL(*res1, z0Trk1); - - trackDensityMap = grid.addTrack(params2, mainDensityMap); - auto res2 = grid.getMaxZPosition(mainDensityMap); - BOOST_CHECK(res2.ok()); - // Trk 2 is closer to z-axis and should yield higher density values - // New maximum is therefore at z0Trk2 - BOOST_CHECK_EQUAL(*res2, z0Trk2); - - auto resWidth1 = grid.getMaxZPositionAndWidth(mainDensityMap); - BOOST_CHECK(resWidth1.ok()); - BOOST_CHECK_EQUAL((*resWidth1).first, z0Trk2); - BOOST_CHECK((*resWidth1).second > 0); + // Empty maps + AdaptiveGridTrackDensity::DensityMap mainDensityMap1D; + AdaptiveGridTrackDensity::DensityMap + mainDensityMap2D; + + // Add first track to spatial grid + auto trackDensityMap = grid1D.addTrack(params1, mainDensityMap1D); + auto firstRes = grid1D.getMaxZTPosition(mainDensityMap1D); + BOOST_CHECK(firstRes.ok()); + // Maximum should be at z0Trk1 position ... + BOOST_CHECK_EQUAL((*firstRes).first, z0Trk1); + // ... and the corresponding time should be set to 0 + BOOST_CHECK_EQUAL((*firstRes).second, 0.); + + // Add first track to 2D grid + auto trackDensityMap2D = grid2D.addTrack(params1, mainDensityMap2D); + auto firstRes2D = grid2D.getMaxZTPosition(mainDensityMap2D); + BOOST_CHECK(firstRes2D.ok()); + // Maximum should be at z0Trk1 position ... + BOOST_CHECK_EQUAL((*firstRes2D).first, z0Trk1); + // ... and the corresponding time should be at t0Trk1 + BOOST_CHECK_EQUAL((*firstRes2D).second, t0Trk1); + + // Add second track to spatial grid + trackDensityMap = grid1D.addTrack(params2, mainDensityMap1D); + // Calculate maximum and the corresponding width + auto secondRes = grid1D.getMaxZTPositionAndWidth(mainDensityMap1D); + BOOST_CHECK(secondRes.ok()); + // Trk 2 is closer to z-axis and should thus yield higher density values. + // Therefore, the new maximum is at z0Trk2 ... + BOOST_CHECK_EQUAL((*secondRes).first.first, z0Trk2); + // ... the corresponding time should be set to 0... + BOOST_CHECK_EQUAL((*secondRes).first.second, 0.); + // ... and it should have a positive width + BOOST_CHECK((*secondRes).second > 0); + + // Add second track to 2D grid + trackDensityMap = grid2D.addTrack(params2, mainDensityMap2D); + // Calculate maximum and the corresponding width + auto secondRes2D = grid2D.getMaxZTPositionAndWidth(mainDensityMap2D); + BOOST_CHECK(secondRes2D.ok()); + // Trk 2 is closer to z-axis and should thus yield higher density values. + // Therefore, the new maximum is at z0Trk2 ... + BOOST_CHECK_EQUAL((*secondRes2D).first.first, z0Trk2); + // ... the corresponding time should be at t0Trk2 ... + BOOST_CHECK_EQUAL((*secondRes2D).first.second, t0Trk2); + // ... and it should have approximately the same width in z direction + CHECK_CLOSE_OR_SMALL((*secondRes2D).second, (*secondRes).second, 1e-5, 1e-5); } -BOOST_AUTO_TEST_CASE(adaptive_gaussian_grid_density_highest_density_sum_test) { - const int trkGridSize = 29; +BOOST_AUTO_TEST_CASE(highest_density_sum) { + const int spatialTrkGridSize = 29; - double binSize = 0.05; // mm + double binExtent = 0.05; // mm - AdaptiveGridTrackDensity::Config cfg(binSize); + AdaptiveGridTrackDensity::Config cfg(binExtent); cfg.useHighestSumZPosition = true; - AdaptiveGridTrackDensity grid(cfg); + AdaptiveGridTrackDensity grid(cfg); // Create some test tracks Covariance covMat(Covariance::Identity() * 0.005); @@ -294,56 +443,69 @@ BOOST_AUTO_TEST_CASE(adaptive_gaussian_grid_density_highest_density_sum_test) { ParticleHypothesis::pion()); // Empty map - AdaptiveGridTrackDensity::DensityMap mainDensityMap; + AdaptiveGridTrackDensity::DensityMap mainDensityMap; // Fill grid with track densities auto trackDensityMap = grid.addTrack(params1, mainDensityMap); - auto res1 = grid.getMaxZPosition(mainDensityMap); + auto res1 = grid.getMaxZTPosition(mainDensityMap); BOOST_CHECK(res1.ok()); // Maximum should be at z0Trk1 position - BOOST_CHECK_EQUAL(*res1, z0Trk1); + BOOST_CHECK_EQUAL((*res1).first, z0Trk1); // Add second track trackDensityMap = grid.addTrack(params2, mainDensityMap); - auto res2 = grid.getMaxZPosition(mainDensityMap); + auto res2 = grid.getMaxZTPosition(mainDensityMap); BOOST_CHECK(res2.ok()); // Trk 2 is closer to z-axis and should yield higher density values // New maximum is therefore at z0Trk2 - BOOST_CHECK_EQUAL(*res2, z0Trk2); + BOOST_CHECK_EQUAL((*res2).first, z0Trk2); // Add small density values around the maximum of track 1 const float densityToAdd = 0.5; - mainDensityMap.at(4) += densityToAdd; - mainDensityMap.at(6) += densityToAdd; + mainDensityMap.at(std::make_pair(4, 0)) += densityToAdd; + mainDensityMap.at(std::make_pair(6, 0)) += densityToAdd; - auto res3 = grid.getMaxZPosition(mainDensityMap); + auto res3 = grid.getMaxZTPosition(mainDensityMap); BOOST_CHECK(res3.ok()); // Trk 2 still has the highest peak density value, however, the small // added densities for track 1 around its maximum should now lead to // a prediction at z0Trk1 again - BOOST_CHECK_EQUAL(*res3, z0Trk1); + BOOST_CHECK_EQUAL((*res3).first, z0Trk1); } -BOOST_AUTO_TEST_CASE(adaptive_gaussian_grid_density_track_removing_test) { - const int trkGridSize = 29; +BOOST_AUTO_TEST_CASE(track_removing) { + const int spatialTrkGridSize = 29; + const int temporalTrkGridSize = 29; + + int trkGridSize = spatialTrkGridSize * temporalTrkGridSize; - double binSize = 0.05; // mm + // bin extent in z and t direction + double binExtent = 0.05; - AdaptiveGridTrackDensity::Config cfg(binSize); - AdaptiveGridTrackDensity grid(cfg); + // 1D grid + AdaptiveGridTrackDensity::Config cfg1D(binExtent); + AdaptiveGridTrackDensity grid1D(cfg1D); + + // 2D grid + AdaptiveGridTrackDensity::Config + cfg2D(binExtent, binExtent); + AdaptiveGridTrackDensity grid2D( + cfg2D); // Create some test tracks - Covariance covMat(Covariance::Identity()); + Covariance covMat = makeRandomCovariance(); // Define z0 values for test tracks float z0Trk1 = -0.45; + float t0Trk1 = -0.15; float z0Trk2 = -0.25; + float t0Trk2 = t0Trk1; BoundVector paramVec0; - paramVec0 << 0.1, z0Trk1, 0, 0, 0, 0; + paramVec0 << 0.1, z0Trk1, 0, 0, 0, t0Trk1; BoundVector paramVec1; - paramVec1 << 0.1, z0Trk2, 0, 0, 0, 0; + paramVec1 << 0.1, z0Trk2, 0, 0, 0, t0Trk2; // Create perigee surface std::shared_ptr perigeeSurface = @@ -354,84 +516,119 @@ BOOST_AUTO_TEST_CASE(adaptive_gaussian_grid_density_track_removing_test) { BoundTrackParameters params1(perigeeSurface, paramVec1, covMat, ParticleHypothesis::pion()); - // Empty map - AdaptiveGridTrackDensity::DensityMap mainDensityMap; - - // Add track 0 - auto trackDensityMap0 = grid.addTrack(params0, mainDensityMap); - BOOST_CHECK(not mainDensityMap.empty()); - // Grid size should match trkGridSize - BOOST_CHECK_EQUAL(mainDensityMap.size(), trkGridSize); - - // Calculate total density - float densitySum0 = 0; - for (auto it = mainDensityMap.begin(); it != mainDensityMap.end(); it++) { - densitySum0 += it->second; - } + // Empty maps + AdaptiveGridTrackDensity::DensityMap mainDensityMap1D; + AdaptiveGridTrackDensity::DensityMap + mainDensityMap2D; + + // Lambda for calculating total density + auto densitySum = [](const auto& densityMap) { + float sum = 0.; + for (auto it = densityMap.begin(); it != densityMap.end(); it++) { + sum += it->second; + } + return sum; + }; - // Add track 0 again - trackDensityMap0 = grid.addTrack(params0, mainDensityMap); - BOOST_CHECK(not mainDensityMap.empty()); + // Add track 0 to 1D grid + auto firstTrackDensityMap1D = grid1D.addTrack(params0, mainDensityMap1D); + BOOST_CHECK(not mainDensityMap1D.empty()); + // Grid size should match spatialTrkGridSize + BOOST_CHECK_EQUAL(mainDensityMap1D.size(), spatialTrkGridSize); + float firstDensitySum1D = densitySum(mainDensityMap1D); + + // Add track 0 to 2D grid + auto firstTrackDensityMap2D = grid2D.addTrack(params0, mainDensityMap2D); + BOOST_CHECK(not mainDensityMap2D.empty()); + // Grid size should match spatialTrkGridSize + BOOST_CHECK_EQUAL(mainDensityMap2D.size(), trkGridSize); + float firstDensitySum2D = densitySum(mainDensityMap2D); + + // Add track 0 again to 1D grid + firstTrackDensityMap1D = grid1D.addTrack(params0, mainDensityMap1D); + BOOST_CHECK(not mainDensityMap1D.empty()); + // Grid size should still match spatialTrkGridSize + BOOST_CHECK_EQUAL(mainDensityMap1D.size(), spatialTrkGridSize); + // Calculate new total density ... + float secondDensitySum1D = densitySum(mainDensityMap1D); + // ... and check that it's twice as large as before + BOOST_CHECK(2 * firstDensitySum1D == secondDensitySum1D); + + // Add track 0 again to 2D grid + firstTrackDensityMap2D = grid2D.addTrack(params0, mainDensityMap2D); + BOOST_CHECK(not mainDensityMap2D.empty()); // Grid size should still match trkGridSize - BOOST_CHECK_EQUAL(mainDensityMap.size(), trkGridSize); - + BOOST_CHECK_EQUAL(mainDensityMap2D.size(), trkGridSize); + // Calculate new total density ... + float secondDensitySum2D = densitySum(mainDensityMap2D); + // ... and check that it's twice as large as before + BOOST_CHECK(2 * firstDensitySum2D == secondDensitySum2D); + + // Remove track 0 from 1D grid + grid1D.subtractTrack(firstTrackDensityMap1D, mainDensityMap1D); // Calculate new total density - float densitySum1 = 0; - for (auto it = mainDensityMap.begin(); it != mainDensityMap.end(); it++) { - densitySum1 += it->second; - } - - BOOST_CHECK(2 * densitySum0 == densitySum1); - - // Remove track 0 - grid.subtractTrack(trackDensityMap0, mainDensityMap); + float thirdDensitySum1D = densitySum(mainDensityMap1D); + // Density should be old one again + BOOST_CHECK(firstDensitySum1D == thirdDensitySum1D); + // Grid size should still match spatialTrkGridSize (removal does not change + // grid size) + BOOST_CHECK_EQUAL(mainDensityMap1D.size(), spatialTrkGridSize); + // Remove track 0 from 2D grid + grid2D.subtractTrack(firstTrackDensityMap2D, mainDensityMap2D); // Calculate new total density - float densitySum2 = 0; - for (auto it = mainDensityMap.begin(); it != mainDensityMap.end(); it++) { - densitySum2 += it->second; - } - + float thirdDensitySum2D = densitySum(mainDensityMap2D); // Density should be old one again - BOOST_CHECK(densitySum0 == densitySum2); - // Grid size should still match trkGridSize (removal does not touch grid size) - BOOST_CHECK_EQUAL(mainDensityMap.size(), trkGridSize); - - // Add track 1, overlapping track 0 - auto trackDensityMap1 = grid.addTrack(params1, mainDensityMap); - - int nNonOverlappingBins = int(std::abs(z0Trk1 - z0Trk2) / binSize + 1); - BOOST_CHECK_EQUAL(mainDensityMap.size(), trkGridSize + nNonOverlappingBins); - - float densitySum3 = 0; - for (auto it = mainDensityMap.begin(); it != mainDensityMap.end(); it++) { - densitySum3 += it->second; - } - - // Remove second track 0 - grid.subtractTrack(trackDensityMap0, mainDensityMap); - - float densitySum4 = 0; - for (auto it = mainDensityMap.begin(); it != mainDensityMap.end(); it++) { - densitySum4 += it->second; - } - + BOOST_CHECK(firstDensitySum2D == thirdDensitySum2D); + // Grid size should still match trkGridSize (removal does not change grid + // size) + BOOST_CHECK_EQUAL(mainDensityMap2D.size(), trkGridSize); + + // Add track 1 to 1D grid (overlaps with track 0!) + auto secondTrackDensityMap1D = grid1D.addTrack(params1, mainDensityMap1D); + int nNonOverlappingBins1D = int(std::abs(z0Trk1 - z0Trk2) / binExtent + 1); + BOOST_CHECK_EQUAL(mainDensityMap1D.size(), + spatialTrkGridSize + nNonOverlappingBins1D); + float fourthDensitySum1D = densitySum(mainDensityMap1D); + + // Add track 1 to 2D grid (overlaps with track 0!) + auto secondTrackDensityMap2D = grid2D.addTrack(params1, mainDensityMap2D); + int nNonOverlappingBins2D = nNonOverlappingBins1D * temporalTrkGridSize; + BOOST_CHECK_EQUAL(mainDensityMap2D.size(), + trkGridSize + nNonOverlappingBins2D); + float fourthDensitySum2D = densitySum(mainDensityMap2D); + + // Remove second track 0 from 1D grid + grid1D.subtractTrack(firstTrackDensityMap1D, mainDensityMap1D); + float fifthDensitySum1D = densitySum(mainDensityMap1D); // Density should match differences of removed tracks - CHECK_CLOSE_ABS(densitySum4, densitySum3 - densitySum0, 1e-5); + CHECK_CLOSE_REL(fifthDensitySum1D, fourthDensitySum1D - firstDensitySum1D, + 1e-5); - // Remove track 1 - grid.subtractTrack(trackDensityMap1, mainDensityMap); + // Remove second track 0 from 2D grid + grid2D.subtractTrack(firstTrackDensityMap2D, mainDensityMap2D); + float fifthDensitySum2D = densitySum(mainDensityMap2D); + // Density should match differences of removed tracks + CHECK_CLOSE_REL(fifthDensitySum2D, fourthDensitySum2D - firstDensitySum2D, + 1e-5); + // Remove track 1 from 1D grid + grid1D.subtractTrack(secondTrackDensityMap1D, mainDensityMap1D); // Size should not have changed - BOOST_CHECK_EQUAL(mainDensityMap.size(), trkGridSize + nNonOverlappingBins); - - float densitySum5 = 0; - for (auto it = mainDensityMap.begin(); it != mainDensityMap.end(); it++) { - densitySum5 += it->second; - } - - // Grid is now empty after all tracks were removed - CHECK_CLOSE_ABS(densitySum5, 0., 1e-5); + BOOST_CHECK_EQUAL(mainDensityMap1D.size(), + spatialTrkGridSize + nNonOverlappingBins1D); + float sixthDensitySum1D = densitySum(mainDensityMap1D); + // 1D grid is now empty after all tracks were removed + CHECK_CLOSE_ABS(sixthDensitySum1D, 0., 1e-4); + + // Remove track 1 from 2D grid + grid2D.subtractTrack(secondTrackDensityMap2D, mainDensityMap2D); + // Size should not have changed + BOOST_CHECK_EQUAL(mainDensityMap2D.size(), + trkGridSize + nNonOverlappingBins2D); + float sixthDensitySum2D = densitySum(mainDensityMap2D); + // 2D grid is now empty after all tracks were removed + CHECK_CLOSE_ABS(sixthDensitySum2D, 0., 1e-4); } } // namespace Test diff --git a/Tests/UnitTests/Core/Vertexing/ImpactPointEstimatorTests.cpp b/Tests/UnitTests/Core/Vertexing/ImpactPointEstimatorTests.cpp index d812db2f0f6..ada4d758d16 100644 --- a/Tests/UnitTests/Core/Vertexing/ImpactPointEstimatorTests.cpp +++ b/Tests/UnitTests/Core/Vertexing/ImpactPointEstimatorTests.cpp @@ -25,6 +25,7 @@ #include "Acts/MagneticField/NullBField.hpp" #include "Acts/Propagator/EigenStepper.hpp" #include "Acts/Propagator/Propagator.hpp" +#include "Acts/Propagator/StraightLineStepper.hpp" #include "Acts/Propagator/detail/VoidPropagatorComponents.hpp" #include "Acts/Surfaces/PerigeeSurface.hpp" #include "Acts/Surfaces/Surface.hpp" @@ -54,10 +55,13 @@ using namespace Acts::UnitLiterals; using Acts::VectorHelpers::makeVector4; using MagneticField = Acts::ConstantBField; +using StraightPropagator = Acts::Propagator; using Stepper = Acts::EigenStepper<>; using Propagator = Acts::Propagator; using Estimator = Acts::ImpactPointEstimator; +using StraightLineEstimator = + Acts::ImpactPointEstimator; const Acts::GeometryContext geoContext; const Acts::MagneticFieldContext magFieldContext; @@ -95,7 +99,7 @@ Estimator makeEstimator(double bZ) { Estimator::Config cfg(field, std::make_shared( std::move(stepper), detail::VoidNavigator(), - getDefaultLogger("Prop", Logging::Level::VERBOSE))); + getDefaultLogger("Prop", Logging::Level::WARNING))); return Estimator(cfg); } @@ -132,7 +136,7 @@ BOOST_AUTO_TEST_SUITE(VertexingImpactPointEstimator) // Check `calculateDistance`, `estimate3DImpactParameters`, and // `getVertexCompatibility`. -BOOST_DATA_TEST_CASE(SingleTrackDistanceParametersCompatibility3d, tracks, d0, +BOOST_DATA_TEST_CASE(SingleTrackDistanceParametersCompatibility3D, tracks, d0, l0, t0, phi, theta, p, q) { auto particleHypothesis = ParticleHypothesis::pion(); @@ -193,12 +197,24 @@ BOOST_DATA_TEST_CASE(SingleTrackDistanceParametersCompatibility3d, tracks, d0, BOOST_DATA_TEST_CASE(TimeAtPca, tracksWithoutIPs* vertices, t0, phi, theta, p, q, vx0, vy0, vz0, vt0) { using Propagator = Acts::Propagator; + using StraightPropagator = Acts::Propagator; + + // Set up quantities for constant B field auto field = std::make_shared(Vector3(0, 0, 2_T)); Stepper stepper(field); auto propagator = std::make_shared(std::move(stepper)); Estimator::Config cfg(field, propagator); Estimator ipEstimator(cfg); - Estimator::State state(magFieldCache()); + Estimator::State ipState(magFieldCache()); + + // Set up quantities for B = 0 + auto zeroField = std::make_shared(Vector3(0, 0, 0)); + StraightLineStepper straightLineStepper; + auto straightLinePropagator = + std::make_shared(straightLineStepper); + StraightLineEstimator::Config zeroFieldCfg(zeroField, straightLinePropagator); + StraightLineEstimator zeroFieldIPEstimator(zeroFieldCfg); + StraightLineEstimator::State zeroFieldIPState(magFieldCache()); // Vertex position and vertex object Vector4 vtxPos(vx0, vy0, vz0, vt0); @@ -240,7 +256,7 @@ BOOST_DATA_TEST_CASE(TimeAtPca, tracksWithoutIPs* vertices, t0, phi, theta, p, // Perigee surface at vertex position auto refPerigeeSurface = Surface::makeShared(refPoint); - // Propagate to the 2D PCA of the reference point + // Set up the propagator options (they are the same with and without B field) PropagatorOptions pOptions(geoContext, magFieldContext); auto intersection = refPerigeeSurface ->intersect(geoContext, params.position(geoContext), @@ -248,13 +264,21 @@ BOOST_DATA_TEST_CASE(TimeAtPca, tracksWithoutIPs* vertices, t0, phi, theta, p, .closest(); pOptions.direction = Direction::fromScalarZeroAsPositive(intersection.pathLength()); + + // Propagate to the 2D PCA of the reference point in a constant B field auto result = propagator->propagate(params, *refPerigeeSurface, pOptions); BOOST_CHECK(result.ok()); - const auto& refParams = *result->endParameters; + // Propagate to the 2D PCA of the reference point when B = 0 + auto zeroFieldResult = + straightLinePropagator->propagate(params, *refPerigeeSurface, pOptions); + BOOST_CHECK(zeroFieldResult.ok()); + const auto& zeroFieldRefParams = *zeroFieldResult->endParameters; + BOOST_TEST_CONTEXT( - "Check time at 2D PCA (i.e., function getImpactParameters)") { + "Check time at 2D PCA (i.e., function getImpactParameters) for helical " + "tracks") { // Calculate impact parameters auto ipParams = ipEstimator .getImpactParameters(refParams, vtx, geoContext, @@ -270,14 +294,14 @@ BOOST_DATA_TEST_CASE(TimeAtPca, tracksWithoutIPs* vertices, t0, phi, theta, p, 1e-3); } - BOOST_TEST_CONTEXT( - "Check time at 3D PCA (i.e., function getDistanceAndMomentum)") { + auto checkGetDistanceAndMomentum = [&vtxPos, &corrMomDir, corrTimeDiff]( + const auto& ipe, const auto& rParams, + auto& state) { // Find 4D distance and momentum of the track at the vertex starting from // the perigee representation at the reference position - auto distAndMom = - ipEstimator - .getDistanceAndMomentum<4>(geoContext, refParams, vtxPos, state) - .value(); + auto distAndMom = ipe.template getDistanceAndMomentum<4>( + geoContext, rParams, vtxPos, state) + .value(); ActsVector<4> distVec = distAndMom.first; Vector3 momDir = distAndMom.second; @@ -293,6 +317,18 @@ BOOST_DATA_TEST_CASE(TimeAtPca, tracksWithoutIPs* vertices, t0, phi, theta, p, // Momentum direction should correspond to the momentum direction at the // vertex CHECK_CLOSE_OR_SMALL(momDir, corrMomDir, 1e-5, 1e-4); + }; + + BOOST_TEST_CONTEXT( + "Check time at 3D PCA (i.e., function getDistanceAndMomentum) for " + "straight tracks") { + checkGetDistanceAndMomentum(zeroFieldIPEstimator, zeroFieldRefParams, + zeroFieldIPState); + } + BOOST_TEST_CONTEXT( + "Check time at 3D PCA (i.e., function getDistanceAndMomentum) for " + "helical tracks") { + checkGetDistanceAndMomentum(ipEstimator, refParams, ipState); } }