src/Vertexing/HoughVertexFinder2.cpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 #include "Acts/Vertexing/HoughVertexFinder2.hpp"
0010
0011 #include "Acts/Seeding/HoughTransformUtils.hpp"
0012 #include "Acts/Utilities/AxisDefinitions.hpp"
0013 #include "Acts/Utilities/Grid.hpp"
0014
0015 #include <numeric>
0016
0017 namespace Acts {
0018
0019 namespace {
0020
0021 /// Type alias for count values in Hough histogram bins
0022 using HoughCount = std::uint16_t;
0023 /// Type alias for equidistant axis used in Hough transform
0024 using HoughAxis = Axis<AxisType::Equidistant, AxisBoundaryType::Open>;
0025 /// Type alias for 2D histogram used in Hough transform
0026 using HoughHist = Grid<HoughCount, HoughAxis, HoughAxis>;
0027
0028 }  // namespace
0029
0030 HoughVertexFinder2::HoughVertexFinder2(Config cfg,
0031                                        std::unique_ptr<const Logger> lgr)
0032     : m_cfg(std::move(cfg)), m_logger(std::move(lgr)) {
0033   if (m_cfg.absEtaFractions.size() != m_cfg.absEtaRanges.size()) {
0034     throw std::invalid_argument("size of the absEtaFractions is " +
0035                                 std::to_string(m_cfg.absEtaFractions.size()) +
0036                                 " but size of the absEtaRanges vector is " +
0037                                 std::to_string(m_cfg.absEtaRanges.size()) +
0038                                 "; these two have to be equal.");
0039   }
0040
0041   if (m_cfg.rangeIterZ.size() != m_cfg.nBinsZIterZ.size()) {
0042     throw std::invalid_argument("size of the rangeIterZ is " +
0043                                 std::to_string(m_cfg.rangeIterZ.size()) +
0044                                 " but size of the nBinsZIterZ vector is " +
0045                                 std::to_string(m_cfg.nBinsZIterZ.size()) +
0046                                 "; these two have to be equal.");
0047   }
0048
0049   if (m_cfg.rangeIterZ.size() != m_cfg.nBinsCotThetaIterZ.size()) {
0050     throw std::invalid_argument(
0051         "size of the rangeIterZ is " + std::to_string(m_cfg.rangeIterZ.size()) +
0052         " but size of the nBinsCotThetaIterZ vector is " +
0053         std::to_string(m_cfg.nBinsCotThetaIterZ.size()) +
0054         "; these two have to be equal.");
0055   }
0056 }
0057
0058 Result<Vector3> HoughVertexFinder2::find(
0059     const SpacePointContainer2& spacePoints) const {
0060   if (spacePoints.empty()) {
0061     return Result<Vector3>::failure(std::error_code());
0062   }
0063
0064   double absEtaRange = m_cfg.maxAbsEta;
0065   double totalFrac = 0.;
0066   for (unsigned int r = 0; r < m_cfg.absEtaRanges.size(); ++r) {
0067     double addToTotalFrac = m_cfg.absEtaFractions.at(r);
0068     if ((totalFrac + addToTotalFrac) * spacePoints.size() > m_cfg.targetSPs) {
0069       double needOnly = (m_cfg.targetSPs - totalFrac * spacePoints.size()) /
0070                         (addToTotalFrac * spacePoints.size());
0071       absEtaRange = (r != 0u ? m_cfg.absEtaRanges.at(r - 1) : 0.) +
0072                     (m_cfg.absEtaRanges.at(r) -
0073                      (r != 0u ? m_cfg.absEtaRanges.at(r - 1) : 0.)) *
0074                         needOnly;
0075       totalFrac += needOnly * addToTotalFrac;
0076       break;
0077     }
0078
0079     totalFrac += addToTotalFrac;
0080   }
0081   if (absEtaRange > m_cfg.maxAbsEta) {
0082     absEtaRange = m_cfg.maxAbsEta;
0083   }
0084   if (absEtaRange < m_cfg.minAbsEta) {
0085     absEtaRange = m_cfg.minAbsEta;
0086   }
0087
0088   const double maxCotTheta = std::sinh(absEtaRange);
0089   const double minCotTheta = -1. * maxCotTheta;
0090
0091   double binsNumDecrease = 1.;
0092   if (spacePoints.size() * totalFrac < m_cfg.targetSPs) {
0093     binsNumDecrease =
0094         std::pow(m_cfg.binsCotThetaDecrease,
0095                  std::log(m_cfg.targetSPs / (spacePoints.size() * totalFrac)));
0096   }
0097
0098   Vector3 vtx{m_cfg.defVtxPosition};
0099   for (std::uint32_t iter = 0; iter < m_cfg.rangeIterZ.size(); ++iter) {
0100     auto vtxNew = findHoughVertex(
0101         spacePoints, vtx, m_cfg.rangeIterZ.at(iter), m_cfg.nBinsZIterZ.at(iter),
0102         minCotTheta, maxCotTheta,
0103         static_cast<std::uint32_t>(m_cfg.nBinsCotThetaIterZ.at(iter) /
0104                                    binsNumDecrease));
0105
0106     if (!vtxNew.ok()) {
0107       // vertex not found
0108       return Result<Vector3>::failure(std::error_code());
0109     }
0110
0111     vtx = vtxNew.value();
0112   }
0113
0114   return Result<Vector3>::success(vtx);
0115 }
0116
0117 Result<Vector3> HoughVertexFinder2::findHoughVertex(
0118     const SpacePointContainer2& spacePoints, const Vector3& vtxOld,
0119     double rangeZ, std::uint32_t numZBins, double minCotTheta,
0120     double maxCotTheta, std::uint32_t numCotThetaBins) const {
0121   const double zBinSize = 2. * rangeZ / numZBins;
0122   const double invCotThetaBinSize =
0123       numCotThetaBins / (maxCotTheta - minCotTheta);
0124   const double minZ = vtxOld[2] - rangeZ;
0125   const double maxZ = vtxOld[2] + rangeZ;
0126   const double vtxOldX = vtxOld[0];
0127   const double vtxOldY = vtxOld[1];
0128
0129   HoughHist houghHist(HoughAxis(minZ, maxZ, numZBins),
0130                       HoughAxis(minCotTheta, maxCotTheta, numCotThetaBins));
0131
0132   std::vector<std::uint32_t> houghZProjection(numZBins, 0);
0133
0134   std::vector<double> vtxZPositions;
0135   for (std::uint32_t zBin = 0; zBin < numZBins; zBin++) {
0136     vtxZPositions.push_back(
0137         HoughTransformUtils::binCenter(minZ, maxZ, numZBins, zBin));
0138   }
0139
0140   for (const auto& sp : spacePoints) {
0141     double sp_invr = 1. / std::hypot((sp.x() - vtxOldX), (sp.y() - vtxOldY));
0142     if (sp.z() > maxZ) {
0143       if ((sp.z() - maxZ + 0.5 * zBinSize) * sp_invr > maxCotTheta) {
0144         continue;
0145       }
0146     } else if (sp.z() < minZ) {
0147       if ((sp.z() - minZ - 0.5 * zBinSize) * sp_invr < minCotTheta) {
0148         continue;
0149       }
0150     }
0151
0152     std::uint32_t zFrom = static_cast<std::uint32_t>(
0153         std::max(((sp.z() - maxCotTheta / sp_invr) - minZ) / zBinSize + 1, 0.));
0154     std::uint32_t zTo = static_cast<std::uint32_t>(std::min(
0155         ((sp.z() - minCotTheta / sp_invr) - minZ) / zBinSize, 1. * numZBins));
0156
0157     for (std::uint32_t zBin = zFrom; zBin < zTo; zBin++) {
0158       double cotTheta = (sp.z() - vtxZPositions[zBin]) * sp_invr;
0159
0160       std::uint32_t cotThetaBin = static_cast<std::uint32_t>(
0161           (cotTheta - minCotTheta) * invCotThetaBinSize);
0162
0163       std::uint32_t cotThetaFrom = static_cast<std::uint32_t>(
0164           std::max<int>(cotThetaBin - m_cfg.fillNeighbours, 0));
0165
0166       std::uint32_t cotThetaTo =
0167           std::min(cotThetaBin + m_cfg.fillNeighbours + 1u, numCotThetaBins);
0168
0169       for (std::uint32_t cotBin = cotThetaFrom; cotBin < cotThetaTo; ++cotBin) {
0170         ++houghHist.atLocalBins({zBin, cotBin});
0171       }
0172     }
0173   }
0174
0175   for (std::uint32_t zBin = 0; zBin < numZBins; zBin++) {
0176     for (std::uint32_t cotBin = 0; cotBin < numCotThetaBins; ++cotBin) {
0177       auto rhs = houghHist.atLocalBins({zBin, cotBin});
0178       if (rhs >= m_cfg.minHits) {
0179         houghZProjection[zBin] += static_cast<std::uint32_t>(rhs);
0180       }
0181     }
0182   }
0183
0184   auto vtxNewZ = findHoughPeak(houghZProjection, vtxZPositions);
0185   if (vtxNewZ.ok()) {
0186     Vector3 newVertex{vtxOldX, vtxOldY, vtxNewZ.value()};
0187     return Result<Vector3>::success(newVertex);
0188   }
0189
0190   return Result<Vector3>::failure(std::error_code());
0191 }
0192
0193 Result<double> HoughVertexFinder2::findHoughPeak(
0194     const std::vector<std::uint32_t>& houghZProjection,
0195     const std::vector<double>& vtxZPositions) const {
0196   std::uint32_t numZBins = houghZProjection.size();
0197
0198   auto maxZElement =
0199       std::max_element(houghZProjection.begin(), houghZProjection.end());
0200   std::uint32_t maxZBin = std::distance(houghZProjection.begin(), maxZElement);
0201
0202   double avg =
0203       std::accumulate(houghZProjection.begin(), houghZProjection.end(), 0.) /
0204       houghZProjection.size();
0205   double sumEntries = 0;
0206   double meanZPeak = 0.;
0207
0208   for (std::uint32_t zBin =
0209            std::max<std::uint32_t>(maxZBin - m_cfg.peakWidth, 0);
0210        zBin <= std::min(numZBins - 1, maxZBin + m_cfg.peakWidth); ++zBin) {
0211     double countsInBin = std::max(houghZProjection.at(zBin) - avg, 0.);
0212     sumEntries += countsInBin;
0213     meanZPeak += vtxZPositions[zBin] * countsInBin;
0214   }
0215
0216   if (sumEntries != 0.) {
0217     meanZPeak /= sumEntries;
0218     return Result<double>::success(meanZPeak);
0219   }
0220
0221   // vertex not found; Hough image empty
0222   return Result<double>::failure(std::error_code());
0223 }
0224
0225 }  // namespace Acts