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