Skip to content

Commit

Permalink
GaussianGridTrackDensity and GridDensityVertexFinder move to cpp files
Browse files Browse the repository at this point in the history
  • Loading branch information
@paulgessinger
paulgessinger committed Feb 27, 2024
1 parent c0c26e8 commit 7ea5362
Show file tree
Hide file tree
Showing 8 changed files with 129 additions and 132 deletions.
43 changes: 27 additions & 16 deletions Core/include/Acts/Vertexing/GaussianGridTrackDensity.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,32 +21,43 @@ namespace Acts {
/// 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.
///
/// @tparam mainGridSize The size of the z-axis 1-dim main density 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 <int mainGridSize = 2000, int trkGridSize = 15>
class GaussianGridTrackDensity {
// Assert odd trkGridSize
static_assert(trkGridSize % 2);
// Assert bigger main grid than track grid
static_assert(mainGridSize > trkGridSize);

public:
using MainGridVector = Eigen::Matrix<float, mainGridSize, 1>;
using TrackGridVector = Eigen::Matrix<float, trkGridSize, 1>;
using MainGridVector = Eigen::Matrix<float, Eigen::Dynamic, 1>;
using TrackGridVector = Eigen::Matrix<float, Eigen::Dynamic, 1>;

/// The configuration struct
struct Config {
/// @param zMinMax_ The minimum and maximum z-values (in mm) that
/// should be covered by the main 1-dim density grid along
/// the z-axis
/// @tparam mainGridSize The size of the z-axis 1-dim main density 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.
/// @note The value of @p zMinMax_ together with @p mainGridSize determines the
/// overall bin size to be used as seen below
Config(float zMinMax_ = 100) : zMinMax(zMinMax_) {
Config(float zMinMax_ = 100, int mainGridSize_ = 2000,
int trkGridSize_ = 15)
: mainGridSize(mainGridSize_),
trkGridSize(trkGridSize_),
zMinMax(zMinMax_) {
binSize = 2. * zMinMax / mainGridSize;

if (trkGridSize % 2 == 0) {
throw std::runtime_error(
"GaussianGridTrackDensity: trkGridSize has to be an odd value!");
}
if (mainGridSize < trkGridSize) {
throw std::runtime_error(
"GaussianGridTrackDensity: mainGridSize has to be bigger than "
"trkGridSize!");
}
}

int mainGridSize;
int trkGridSize;

// Min and max z value of big grid
float zMinMax; // mm

Expand Down Expand Up @@ -105,6 +116,8 @@ class GaussianGridTrackDensity {
void removeTrackGridFromMainGrid(int zBin, const TrackGridVector& trkGrid,
MainGridVector& mainGrid) const;

const Config& config() const { return m_cfg; }

private:
/// @brief Helper function that actually adds the track to the
/// main density grid
Expand Down Expand Up @@ -176,5 +189,3 @@ class GaussianGridTrackDensity {
};

} // namespace Acts

#include "Acts/Vertexing/GaussianGridTrackDensity.ipp"
38 changes: 15 additions & 23 deletions Core/include/Acts/Vertexing/GridDensityVertexFinder.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,33 +29,23 @@ namespace Acts {
/// density vector) needs to be calculated. All track contributions along the
/// beam axis (main density grid) a superimposed and the z-value of the bin
/// with the highest track density is returned as a vertex candidate.
/// @tparam mainGridSize The size of the z-axis 1-dim main density 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 <int mainGridSize = 2000, int trkGridSize = 15>
class GridDensityVertexFinder final : public IVertexFinder {
// Assert odd trkGridSize
static_assert(trkGridSize % 2);
// Assert bigger main grid than track grid
static_assert(mainGridSize > trkGridSize);

using GridDensity = GaussianGridTrackDensity<mainGridSize, trkGridSize>;

public:
using MainGridVector = typename GridDensity::MainGridVector;
using TrackGridVector = typename GridDensity::TrackGridVector;
using MainGridVector = GaussianGridTrackDensity::MainGridVector;
using TrackGridVector = GaussianGridTrackDensity::TrackGridVector;

/// @brief The Config struct
struct Config {
///@param zMinMax min and max z value of big z-axis grid
Config(float zMinMax = 100)
: gridDensity(typename GridDensity::Config(zMinMax)) {}
Config(float zMinMax_ = 100, int mainGridSize_ = 2000,
int trkGridSize_ = 15)
: Config{GaussianGridTrackDensity::Config{zMinMax_, mainGridSize_,
trkGridSize_}} {}

///@param gDensity The grid density
Config(const GridDensity& gDensity) : gridDensity(gDensity) {}
Config(GaussianGridTrackDensity gDensity) : gridDensity(gDensity) {}

// The grid density object
GridDensity gridDensity;
GaussianGridTrackDensity gridDensity;

// Cache the main grid and the density contributions (trackGrid and z-bin)
// for every single track.
Expand All @@ -82,8 +72,10 @@ class GridDensityVertexFinder final : public IVertexFinder {
///
/// Only needed if cacheGridStateForTrackRemoval == true
struct State {
State(MainGridVector mainGrid_) : mainGrid(std::move(mainGrid_)) {}

// The main density grid
MainGridVector mainGrid = MainGridVector::Zero();
MainGridVector mainGrid;
// Map to store z-bin and track grid (i.e. the density contribution of
// a single track to the main grid) for every single track
std::map<InputTrack, std::pair<int, TrackGridVector>> binAndTrackGridMap;
Expand Down Expand Up @@ -115,7 +107,9 @@ class GridDensityVertexFinder final : public IVertexFinder {

IVertexFinder::State makeState(
const Acts::MagneticFieldContext& /*mctx*/) const override {
return IVertexFinder::State{State{}};
return IVertexFinder::State{
std::in_place_type<State>,
MainGridVector{m_cfg.gridDensity.config().mainGridSize}};
}

void setTracksToRemove(
Expand Down Expand Up @@ -151,5 +145,3 @@ class GridDensityVertexFinder final : public IVertexFinder {
};

} // namespace Acts

#include "GridDensityVertexFinder.ipp"
2 changes: 2 additions & 0 deletions Core/src/Vertexing/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,4 +17,6 @@ target_sources(
TrackDensityVertexFinder.cpp
GaussianTrackDensity.cpp
ImpactPointEstimator.cpp
GaussianGridTrackDensity.cpp
GridDensityVertexFinder.cpp
)
110 changes: 49 additions & 61 deletions ...ts/Vertexing/GaussianGridTrackDensity.ipp → ...rc/Vertexing/GaussianGridTrackDensity.cpp
View file
Original file line number Diff line number Diff line change
@@ -1,20 +1,22 @@
// This file is part of the Acts project.
//
// Copyright (C) 2020 CERN for the benefit of the Acts project
// Copyright (C) 2020-2023 CERN for the benefit of the Acts project
//
// This Source Code Form is subject to the terms of the Mozilla Public
// License, v. 2.0. If a copy of the MPL was not distributed with this
// file, You can obtain one at http://mozilla.org/MPL/2.0/.

#include "Acts/Vertexing/GaussianGridTrackDensity.hpp"

#include "Acts/Utilities/AlgebraHelpers.hpp"
#include "Acts/Vertexing/VertexingError.hpp"

#include <algorithm>
namespace Acts {

template <int mainGridSize, int trkGridSize>
Acts::Result<float>
Acts::GaussianGridTrackDensity<mainGridSize, trkGridSize>::getMaxZPosition(
Result<float> GaussianGridTrackDensity::getMaxZPosition(
MainGridVector& mainGrid) const {
if (mainGrid == MainGridVector::Zero()) {
if (mainGrid.isZero(0)) {
return VertexingError::EmptyInput;
}

Expand All @@ -29,13 +31,12 @@ Acts::GaussianGridTrackDensity<mainGridSize, trkGridSize>::getMaxZPosition(
}

// Derive corresponding z value
return (zbin - mainGridSize / 2.0f + 0.5f) * m_cfg.binSize;
return (zbin - m_cfg.mainGridSize / 2.0f + 0.5f) * m_cfg.binSize;
}

template <int mainGridSize, int trkGridSize>
Acts::Result<std::pair<float, float>> Acts::GaussianGridTrackDensity<
mainGridSize, trkGridSize>::getMaxZPositionAndWidth(MainGridVector&
mainGrid) const {
Result<std::pair<float, float>>
GaussianGridTrackDensity::getMaxZPositionAndWidth(
MainGridVector& mainGrid) const {
// Get z maximum value
auto maxZRes = getMaxZPosition(mainGrid);
if (!maxZRes.ok()) {
Expand All @@ -53,11 +54,9 @@ Acts::Result<std::pair<float, float>> Acts::GaussianGridTrackDensity<
return returnPair;
}

template <int mainGridSize, int trkGridSize>
std::pair<int, typename Acts::GaussianGridTrackDensity<
mainGridSize, trkGridSize>::TrackGridVector>
Acts::GaussianGridTrackDensity<mainGridSize, trkGridSize>::addTrack(
const Acts::BoundTrackParameters& trk, MainGridVector& mainGrid) const {
std::pair<int, GaussianGridTrackDensity::TrackGridVector>
GaussianGridTrackDensity::addTrack(const BoundTrackParameters& trk,
MainGridVector& mainGrid) const {
SquareMatrix2 cov = trk.spatialImpactParameterCovariance().value();
float d0 = trk.parameters()[0];
float z0 = trk.parameters()[1];
Expand All @@ -66,17 +65,17 @@ Acts::GaussianGridTrackDensity<mainGridSize, trkGridSize>::addTrack(
int dOffset = static_cast<int>(std::floor(d0 / m_cfg.binSize - 0.5) + 1);
// Check if current track does affect grid density
// in central bins at z-axis
if (std::abs(dOffset) > (trkGridSize - 1) / 2.) {
if (std::abs(dOffset) > (m_cfg.trkGridSize - 1) / 2.) {
// Current track is too far away to contribute
// to track density at z-axis bins
return {-1, TrackGridVector::Zero()};
return {-1, TrackGridVector::Zero(m_cfg.trkGridSize)};
}

// Calculate bin in z
int zBin = int(z0 / m_cfg.binSize + mainGridSize / 2.);
int zBin = int(z0 / m_cfg.binSize + m_cfg.mainGridSize / 2.);

if (zBin < 0 || zBin >= mainGridSize) {
return {-1, TrackGridVector::Zero()};
if (zBin < 0 || zBin >= m_cfg.mainGridSize) {
return {-1, TrackGridVector::Zero(m_cfg.trkGridSize)};
}
// Calculate the positions of the bin centers
float binCtrZ = (zBin + 0.5f) * m_cfg.binSize - m_cfg.zMinMax;
Expand All @@ -93,71 +92,61 @@ Acts::GaussianGridTrackDensity<mainGridSize, trkGridSize>::addTrack(
return {zBin, trackGrid};
}

template <int mainGridSize, int trkGridSize>
void Acts::GaussianGridTrackDensity<mainGridSize, trkGridSize>::
addTrackGridToMainGrid(int zBin, const TrackGridVector& trkGrid,
MainGridVector& mainGrid) const {
void GaussianGridTrackDensity::addTrackGridToMainGrid(
int zBin, const TrackGridVector& trkGrid, MainGridVector& mainGrid) const {
modifyMainGridWithTrackGrid(zBin, trkGrid, mainGrid, +1);
}

template <int mainGridSize, int trkGridSize>
void Acts::GaussianGridTrackDensity<mainGridSize, trkGridSize>::
removeTrackGridFromMainGrid(int zBin, const TrackGridVector& trkGrid,
MainGridVector& mainGrid) const {
void GaussianGridTrackDensity::removeTrackGridFromMainGrid(
int zBin, const TrackGridVector& trkGrid, MainGridVector& mainGrid) const {
modifyMainGridWithTrackGrid(zBin, trkGrid, mainGrid, -1);
}

template <int mainGridSize, int trkGridSize>
void Acts::GaussianGridTrackDensity<mainGridSize, trkGridSize>::
modifyMainGridWithTrackGrid(int zBin, const TrackGridVector& trkGrid,
MainGridVector& mainGrid,
int modifyModeSign) const {
int width = (trkGridSize - 1) / 2;
void GaussianGridTrackDensity::modifyMainGridWithTrackGrid(
int zBin, const TrackGridVector& trkGrid, MainGridVector& mainGrid,
int modifyModeSign) const {
int width = (m_cfg.trkGridSize - 1) / 2;
// Overlap left
int leftOL = zBin - width;
// Overlap right
int rightOL = zBin + width - mainGridSize + 1;
int rightOL = zBin + width - m_cfg.mainGridSize + 1;
if (leftOL < 0) {
int totalTrkSize = trkGridSize + leftOL;
int totalTrkSize = m_cfg.trkGridSize + leftOL;
mainGrid.segment(0, totalTrkSize) +=
modifyModeSign * trkGrid.segment(-leftOL, totalTrkSize);
return;
}
if (rightOL > 0) {
int totalTrkSize = trkGridSize - rightOL;
mainGrid.segment(mainGridSize - totalTrkSize, totalTrkSize) +=
int totalTrkSize = m_cfg.trkGridSize - rightOL;
mainGrid.segment(m_cfg.mainGridSize - totalTrkSize, totalTrkSize) +=
modifyModeSign * trkGrid.segment(0, totalTrkSize);
return;
}

mainGrid.segment(zBin - width, trkGridSize) += modifyModeSign * trkGrid;
mainGrid.segment(zBin - width, m_cfg.trkGridSize) += modifyModeSign * trkGrid;
}

template <int mainGridSize, int trkGridSize>
typename Acts::GaussianGridTrackDensity<mainGridSize,
trkGridSize>::TrackGridVector
Acts::GaussianGridTrackDensity<mainGridSize, trkGridSize>::createTrackGrid(
float d0, float distCtrZ, const Acts::SquareMatrix2& cov) const {
TrackGridVector trackGrid(TrackGridVector::Zero());
float floorHalfTrkGridSize = static_cast<float>(trkGridSize) / 2 - 0.5f;
GaussianGridTrackDensity::TrackGridVector
GaussianGridTrackDensity::createTrackGrid(float d0, float distCtrZ,
const SquareMatrix2& cov) const {
TrackGridVector trackGrid(TrackGridVector::Zero(m_cfg.trkGridSize));
float floorHalfTrkGridSize = static_cast<float>(m_cfg.trkGridSize) / 2 - 0.5f;

// Loop over columns
for (int j = 0; j < trkGridSize; j++) {
for (int j = 0; j < m_cfg.trkGridSize; j++) {
float z = (j - floorHalfTrkGridSize) * m_cfg.binSize;
trackGrid(j) = normal2D(-d0, z - distCtrZ, cov);
}
return trackGrid;
}

template <int mainGridSize, int trkGridSize>
Acts::Result<float>
Acts::GaussianGridTrackDensity<mainGridSize, trkGridSize>::estimateSeedWidth(
Result<float> GaussianGridTrackDensity::estimateSeedWidth(
MainGridVector& mainGrid, float maxZ) const {
if (mainGrid == MainGridVector::Zero()) {
if (mainGrid.isZero(0)) {
return VertexingError::EmptyInput;
}
// Get z bin of max density z value
int zBin = int(maxZ / m_cfg.binSize + mainGridSize / 2.);
int zBin = int(maxZ / m_cfg.binSize + m_cfg.mainGridSize / 2.);

const float maxValue = mainGrid(zBin);
float gridValue = mainGrid(zBin);
Expand Down Expand Up @@ -195,9 +184,8 @@ Acts::GaussianGridTrackDensity<mainGridSize, trkGridSize>::estimateSeedWidth(
return std::isnormal(width) ? width : 0.0f;
}

template <int mainGridSize, int trkGridSize>
float Acts::GaussianGridTrackDensity<mainGridSize, trkGridSize>::normal2D(
float d, float z, const Acts::SquareMatrix2& cov) const {
float GaussianGridTrackDensity::normal2D(float d, float z,
const SquareMatrix2& cov) const {
float det = cov.determinant();
float coef = 1 / std::sqrt(det);
float expo =
Expand All @@ -206,9 +194,8 @@ float Acts::GaussianGridTrackDensity<mainGridSize, trkGridSize>::normal2D(
return coef * safeExp(expo);
}

template <int mainGridSize, int trkGridSize>
int Acts::GaussianGridTrackDensity<mainGridSize, trkGridSize>::
getHighestSumZPosition(MainGridVector& mainGrid) const {
int GaussianGridTrackDensity::getHighestSumZPosition(
MainGridVector& mainGrid) 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)

Expand Down Expand Up @@ -255,9 +242,8 @@ int Acts::GaussianGridTrackDensity<mainGridSize, trkGridSize>::
return zFirstMax;
}

template <int mainGridSize, int trkGridSize>
double Acts::GaussianGridTrackDensity<mainGridSize, trkGridSize>::getDensitySum(
const MainGridVector& mainGrid, int pos) const {
double GaussianGridTrackDensity::getDensitySum(const MainGridVector& mainGrid,
int pos) const {
double sum = mainGrid(pos);
// Sum up only the density contributions from the
// neighboring bins if they are still within bounds
Expand All @@ -269,3 +255,5 @@ double Acts::GaussianGridTrackDensity<mainGridSize, trkGridSize>::getDensitySum(
}
return sum;
}

} // namespace Acts
Loading

0 comments on commit 7ea5362

Please sign in to comment.