Acts/Vertexing/HoughVertexFinder.ipp

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 #pragma once
0010
0011 #include "Acts/Vertexing/HoughVertexFinder.hpp"
0012
0013 #include "Acts/Seeding/HoughTransformUtils.hpp"
0014
0015 #include <numeric>
0016
0017 template <typename spacepoint_t>
0018 Acts::HoughVertexFinder<spacepoint_t>::HoughVertexFinder(
0019     Config cfg, std::unique_ptr<const Logger> lgr)
0020     : m_cfg(std::move(cfg)), m_logger(std::move(lgr)) {
0021   if (m_cfg.absEtaFractions.size() != m_cfg.absEtaRanges.size()) {
0022     throw std::invalid_argument("size of the absEtaFractions is " +
0023                                 std::to_string(m_cfg.absEtaFractions.size()) +
0024                                 " but size of the absEtaRanges vector is " +
0025                                 std::to_string(m_cfg.absEtaRanges.size()) +
0026                                 "; these two have to be equal.");
0027   }
0028
0029   if (m_cfg.rangeIterZ.size() != m_cfg.nBinsZIterZ.size()) {
0030     throw std::invalid_argument("size of the rangeIterZ is " +
0031                                 std::to_string(m_cfg.rangeIterZ.size()) +
0032                                 " but size of the nBinsZIterZ vector is " +
0033                                 std::to_string(m_cfg.nBinsZIterZ.size()) +
0034                                 "; these two have to be equal.");
0035   }
0036
0037   if (m_cfg.rangeIterZ.size() != m_cfg.nBinsCotThetaIterZ.size()) {
0038     throw std::invalid_argument(
0039         "size of the rangeIterZ is " + std::to_string(m_cfg.rangeIterZ.size()) +
0040         " but size of the nBinsCotThetaIterZ vector is " +
0041         std::to_string(m_cfg.nBinsCotThetaIterZ.size()) +
0042         "; these two have to be equal.");
0043   }
0044 }
0045
0046 template <typename spacepoint_t>
0047 Acts::Result<Acts::Vector3> Acts::HoughVertexFinder<spacepoint_t>::find(
0048     const std::vector<spacepoint_t>& spacepoints) const {
0049   if (spacepoints.empty()) {
0050     return Acts::Result<Acts::Vector3>::failure(std::error_code());
0051   }
0052
0053   double absEtaRange = m_cfg.maxAbsEta;
0054   double totalFrac = 0.;
0055   for (unsigned int r = 0; r < m_cfg.absEtaRanges.size(); ++r) {
0056     double addToTotalFrac = m_cfg.absEtaFractions.at(r);
0057     if ((totalFrac + addToTotalFrac) * spacepoints.size() > m_cfg.targetSPs) {
0058       double needOnly = (m_cfg.targetSPs - totalFrac * spacepoints.size()) /
0059                         (addToTotalFrac * spacepoints.size());
0060       absEtaRange = (r != 0u ? m_cfg.absEtaRanges.at(r - 1) : 0.) +
0061                     (m_cfg.absEtaRanges.at(r) -
0062                      (r != 0u ? m_cfg.absEtaRanges.at(r - 1) : 0.)) *
0063                         needOnly;
0064       totalFrac += needOnly * addToTotalFrac;
0065       break;
0066     }
0067
0068     totalFrac += addToTotalFrac;
0069   }
0070   if (absEtaRange > m_cfg.maxAbsEta) {
0071     absEtaRange = m_cfg.maxAbsEta;
0072   }
0073   if (absEtaRange < m_cfg.minAbsEta) {
0074     absEtaRange = m_cfg.minAbsEta;
0075   }
0076
0077   const double maxCotTheta = std::sinh(absEtaRange);
0078   const double minCotTheta = -1. * maxCotTheta;
0079
0080   double binsNumDecrease = 1.;
0081   if (spacepoints.size() * totalFrac < m_cfg.targetSPs) {
0082     binsNumDecrease =
0083         std::pow(m_cfg.binsCotThetaDecrease,
0084                  std::log(m_cfg.targetSPs / (spacepoints.size() * totalFrac)));
0085   }
0086
0087   Acts::Vector3 vtx{m_cfg.defVtxPosition};
0088   for (std::uint32_t iter = 0; iter < m_cfg.rangeIterZ.size(); ++iter) {
0089     auto vtxNew = findHoughVertex(
0090         spacepoints, vtx, m_cfg.rangeIterZ.at(iter), m_cfg.nBinsZIterZ.at(iter),
0091         minCotTheta, maxCotTheta,
0092         static_cast<std::uint32_t>(m_cfg.nBinsCotThetaIterZ.at(iter) /
0093                                    binsNumDecrease));
0094
0095     if (!vtxNew.ok()) {
0096       // vertex not found
0097       return Acts::Result<Acts::Vector3>::failure(std::error_code());
0098     }
0099
0100     vtx = vtxNew.value();
0101   }
0102
0103   return Acts::Result<Acts::Vector3>::success(vtx);
0104 }
0105
0106 template <typename spacepoint_t>
0107 Acts::Result<Acts::Vector3>
0108 Acts::HoughVertexFinder<spacepoint_t>::findHoughVertex(
0109     const std::vector<spacepoint_t>& spacepoints, const Acts::Vector3& vtxOld,
0110     double rangeZ, std::uint32_t numZBins, double minCotTheta,
0111     double maxCotTheta, std::uint32_t numCotThetaBins) const {
0112   const double zBinSize = 2. * rangeZ / numZBins;
0113   const double invCotThetaBinSize =
0114       numCotThetaBins / (maxCotTheta - minCotTheta);
0115   const double minZ = vtxOld[2] - rangeZ;
0116   const double maxZ = vtxOld[2] + rangeZ;
0117   const double vtxOldX = vtxOld[0];
0118   const double vtxOldY = vtxOld[1];
0119
0120   HoughHist houghHist(HoughAxis(minZ, maxZ, numZBins),
0121                       HoughAxis(minCotTheta, maxCotTheta, numCotThetaBins));
0122
0123   std::vector<std::uint32_t> houghZProjection(numZBins, 0);
0124
0125   std::vector<double> vtxZPositions;
0126   for (std::uint32_t zBin = 0; zBin < numZBins; zBin++) {
0127     vtxZPositions.push_back(
0128         Acts::HoughTransformUtils::binCenter(minZ, maxZ, numZBins, zBin));
0129   }
0130
0131   for (const auto& sp : spacepoints) {
0132     if (sp.z() > maxZ) {
0133       if ((sp.z() - maxZ) / sp.r() > maxCotTheta) {
0134         continue;
0135       }
0136     } else if (sp.z() < minZ) {
0137       if ((sp.z() - minZ) / sp.r() < minCotTheta) {
0138         continue;
0139       }
0140     }
0141
0142     double sp_invr = 1. / std::hypot((sp.x() - vtxOldX), (sp.y() - vtxOldY));
0143
0144     std::uint32_t zFrom = static_cast<std::uint32_t>(
0145         std::max(((sp.z() - maxCotTheta / sp_invr) - minZ) / zBinSize + 1, 0.));
0146     std::uint32_t zTo = static_cast<std::uint32_t>(std::min(
0147         ((sp.z() - minCotTheta / sp_invr) - minZ) / zBinSize, 1. * numZBins));
0148
0149     for (std::uint32_t zBin = zFrom; zBin < zTo; zBin++) {
0150       double cotTheta = (sp.z() - vtxZPositions[zBin]) * sp_invr;
0151
0152       std::uint32_t cotThetaBin = static_cast<std::uint32_t>(
0153           (cotTheta - minCotTheta) * invCotThetaBinSize);
0154
0155       std::uint32_t cotThetaFrom =
0156           std::max<std::uint32_t>(cotThetaBin - m_cfg.fillNeighbours, 0);
0157       std::uint32_t cotThetaTo =
0158           std::min(cotThetaBin + m_cfg.fillNeighbours + 1, numCotThetaBins);
0159
0160       for (std::uint32_t cotBin = cotThetaFrom; cotBin < cotThetaTo; ++cotBin) {
0161         ++houghHist.atLocalBins({zBin, cotBin});
0162       }
0163     }
0164   }
0165
0166   for (std::uint32_t zBin = 0; zBin < numZBins; zBin++) {
0167     for (std::uint32_t cotBin = 0; cotBin < numCotThetaBins; ++cotBin) {
0168       auto rhs = houghHist.atLocalBins({zBin, cotBin});
0169       houghZProjection[zBin] +=
0170           static_cast<std::uint32_t>(rhs * (rhs >= m_cfg.minHits));
0171     }
0172   }
0173
0174   auto vtxNewZ = findHoughPeak(houghZProjection, vtxZPositions);
0175   if (vtxNewZ.ok()) {
0176     Acts::Vector3 newVertex{vtxOldX, vtxOldY, vtxNewZ.value()};
0177     return Acts::Result<Acts::Vector3>::success(newVertex);
0178   }
0179
0180   return Acts::Result<Acts::Vector3>::failure(std::error_code());
0181 }
0182
0183 template <typename spacepoint_t>
0184 Acts::Result<double> Acts::HoughVertexFinder<spacepoint_t>::findHoughPeak(
0185     const std::vector<std::uint32_t>& houghZProjection,
0186     const std::vector<double>& vtxZPositions) const {
0187   std::uint32_t numZBins = houghZProjection.size();
0188
0189   auto maxZElement =
0190       std::max_element(houghZProjection.begin(), houghZProjection.end());
0191   std::uint32_t maxZBin = std::distance(houghZProjection.begin(), maxZElement);
0192
0193   double avg =
0194       std::accumulate(houghZProjection.begin(), houghZProjection.end(), 0.) /
0195       houghZProjection.size();
0196   double sumEntries = 0;
0197   double meanZPeak = 0.;
0198
0199   for (std::uint32_t zBin =
0200            std::max<std::uint32_t>(maxZBin - m_cfg.peakWidth, 0);
0201        zBin <= std::min(numZBins - 1, maxZBin + m_cfg.peakWidth); ++zBin) {
0202     double countsInBin = std::max(houghZProjection.at(zBin) - avg, 0.);
0203     sumEntries += countsInBin;
0204     meanZPeak += vtxZPositions[zBin] * countsInBin;
0205   }
0206
0207   if (sumEntries != 0.) {
0208     meanZPeak /= sumEntries;
0209     return Acts::Result<double>::success(meanZPeak);
0210   }
0211
0212   // vertex not found; Hough image empty
0213   return Acts::Result<double>::failure(std::error_code());
0214 }