Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-12-16 09:23:18

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/Seeding2/CylindricalSpacePointGrid2.hpp"
0010 
0011 namespace Acts {
0012 
0013 CylindricalSpacePointGrid2::CylindricalSpacePointGrid2(
0014     const Config& config, std::unique_ptr<const Logger> _logger)
0015     : m_cfg(config), m_logger(std::move(_logger)) {
0016   if (m_cfg.phiMin < -std::numbers::pi_v<float> ||
0017       m_cfg.phiMax > std::numbers::pi_v<float>) {
0018     throw std::runtime_error(
0019         "CylindricalSpacePointGrid2: phiMin (" + std::to_string(m_cfg.phiMin) +
0020         ") and/or phiMax (" + std::to_string(m_cfg.phiMax) +
0021         ") are outside the allowed phi range, defined as "
0022         "[-std::numbers::pi_v<float>, std::numbers::pi_v<float>]");
0023   }
0024   if (m_cfg.phiMin > m_cfg.phiMax) {
0025     throw std::runtime_error(
0026         "CylindricalSpacePointGrid2: phiMin is bigger then phiMax");
0027   }
0028   if (m_cfg.rMin > m_cfg.rMax) {
0029     throw std::runtime_error(
0030         "CylindricalSpacePointGrid2: rMin is bigger then rMax");
0031   }
0032   if (m_cfg.zMin > m_cfg.zMax) {
0033     throw std::runtime_error(
0034         "CylindricalSpacePointGrid2: zMin is bigger than zMax");
0035   }
0036 
0037   int phiBins = 0;
0038   if (m_cfg.bFieldInZ == 0) {
0039     // for no magnetic field, use the maximum number of phi bins
0040     phiBins = m_cfg.maxPhiBins;
0041     ACTS_VERBOSE(
0042         "B-Field is 0 (z-coordinate), setting the number of bins in phi to "
0043         << phiBins);
0044   } else {
0045     // calculate circle intersections of helix and max detector radius in mm.
0046     // bFieldInZ is in (pT/radius) natively, no need for conversion
0047     const float minHelixRadius = m_cfg.minPt / m_cfg.bFieldInZ;
0048 
0049     // sanity check: if yOuter takes the square root of a negative number
0050     if (minHelixRadius < m_cfg.rMax * 0.5) {
0051       throw std::domain_error(
0052           "The value of minHelixRadius cannot be smaller than rMax / 2. Please "
0053           "check the m_cfguration of bFieldInZ and minPt");
0054     }
0055 
0056     const float maxR2 = m_cfg.rMax * m_cfg.rMax;
0057     const float xOuter = maxR2 / (2 * minHelixRadius);
0058     const float yOuter = std::sqrt(maxR2 - xOuter * xOuter);
0059     const float outerAngle = std::atan(xOuter / yOuter);
0060     // intersection of helix and max detector radius minus maximum R distance
0061     // from middle SP to top SP
0062     float innerAngle = 0;
0063     float rMin = m_cfg.rMax;
0064     if (m_cfg.rMax > m_cfg.deltaRMax) {
0065       rMin = m_cfg.rMax - m_cfg.deltaRMax;
0066       const float innerCircleR2 =
0067           (m_cfg.rMax - m_cfg.deltaRMax) * (m_cfg.rMax - m_cfg.deltaRMax);
0068       const float xInner = innerCircleR2 / (2 * minHelixRadius);
0069       const float yInner = std::sqrt(innerCircleR2 - xInner * xInner);
0070       innerAngle = std::atan(xInner / yInner);
0071     }
0072 
0073     // evaluating the azimutal deflection including the maximum impact parameter
0074     const float deltaAngleWithMaxD0 =
0075         std::abs(std::asin(m_cfg.impactMax / rMin) -
0076                  std::asin(m_cfg.impactMax / m_cfg.rMax));
0077 
0078     // evaluating delta Phi based on the inner and outer angle, and the azimutal
0079     // deflection including the maximum impact parameter
0080     // Divide by m_cfg.phiBinDeflectionCoverage since we combine
0081     // m_cfg.phiBinDeflectionCoverage number of consecutive phi bins in the
0082     // seed making step. So each individual bin should cover
0083     // 1/m_cfg.phiBinDeflectionCoverage of the maximum expected azimutal
0084     // deflection
0085     const float deltaPhi = (outerAngle - innerAngle + deltaAngleWithMaxD0) /
0086                            m_cfg.phiBinDeflectionCoverage;
0087 
0088     // sanity check: if the delta phi is equal to or less than zero, we'll be
0089     // creating an infinite or a negative number of bins, which would be bad!
0090     if (deltaPhi <= 0.f) {
0091       throw std::domain_error(
0092           "Delta phi value is equal to or less than zero, leading to an "
0093           "impossible number of bins (negative or infinite)");
0094     }
0095 
0096     // divide 2pi by angle delta to get number of phi-bins
0097     // size is always 2pi even for regions of interest
0098     phiBins = static_cast<int>(std::ceil(2 * std::numbers::pi / deltaPhi));
0099     // need to scale the number of phi bins accordingly to the number of
0100     // consecutive phi bins in the seed making step.
0101     // Each individual bin should be approximately a fraction (depending on this
0102     // number) of the maximum expected azimutal deflection.
0103 
0104     // set protection for large number of bins, by default it is large
0105     phiBins = std::min(phiBins, m_cfg.maxPhiBins);
0106   }
0107 
0108   PhiAxisType phiAxis(AxisClosed, m_cfg.phiMin, m_cfg.phiMax, phiBins);
0109 
0110   // vector that will store the edges of the bins of z
0111   std::vector<double> zValues;
0112 
0113   // If zBinEdges is not defined, calculate the edges as zMin + bin * zBinSize
0114   if (m_cfg.zBinEdges.empty()) {
0115     // TODO: can probably be optimized using smaller z bins
0116     // and returning (multiple) neighbors only in one z-direction for forward
0117     // seeds
0118     // FIXME: zBinSize must include scattering
0119     const float zBinSize = m_cfg.cotThetaMax * m_cfg.deltaRMax;
0120     const float zBins =
0121         std::max(1.f, std::floor((m_cfg.zMax - m_cfg.zMin) / zBinSize));
0122 
0123     zValues.reserve(static_cast<int>(zBins));
0124     for (int bin = 0; bin <= static_cast<int>(zBins); bin++) {
0125       const double edge =
0126           m_cfg.zMin + bin * ((m_cfg.zMax - m_cfg.zMin) / zBins);
0127       zValues.push_back(edge);
0128     }
0129   } else {
0130     // Use the zBinEdges defined in the m_cfg
0131     zValues.reserve(m_cfg.zBinEdges.size());
0132     for (float bin : m_cfg.zBinEdges) {
0133       zValues.push_back(bin);
0134     }
0135   }
0136 
0137   std::vector<double> rValues;
0138   rValues.reserve(std::max(2ul, m_cfg.rBinEdges.size()));
0139   if (m_cfg.rBinEdges.empty()) {
0140     rValues = {m_cfg.rMin, m_cfg.rMax};
0141   } else {
0142     rValues.insert(rValues.end(), m_cfg.rBinEdges.begin(),
0143                    m_cfg.rBinEdges.end());
0144   }
0145 
0146   ZAxisType zAxis(AxisOpen, std::move(zValues));
0147   RAxisType rAxis(AxisOpen, std::move(rValues));
0148 
0149   ACTS_VERBOSE("Defining Grid:");
0150   ACTS_VERBOSE("- Phi Axis: " << phiAxis);
0151   ACTS_VERBOSE("- Z axis  : " << zAxis);
0152   ACTS_VERBOSE("- R axis  : " << rAxis);
0153 
0154   GridType grid(
0155       std::make_tuple(std::move(phiAxis), std::move(zAxis), std::move(rAxis)));
0156   m_binnedGroup.emplace(std::move(grid), m_cfg.bottomBinFinder.value(),
0157                         m_cfg.topBinFinder.value(), m_cfg.navigation);
0158   m_grid = &m_binnedGroup->grid();
0159 }
0160 
0161 void CylindricalSpacePointGrid2::clear() {
0162   for (std::size_t i = 0; i < grid().size(); ++i) {
0163     BinType& bin = grid().at(i);
0164     bin.clear();
0165   }
0166   m_counter = 0;
0167 }
0168 
0169 std::optional<std::size_t> CylindricalSpacePointGrid2::insert(
0170     SpacePointIndex index, float phi, float z, float r) {
0171   const std::optional<std::size_t> gridIndex = binIndex(phi, z, r);
0172   if (gridIndex.has_value()) {
0173     BinType& bin = grid().at(*gridIndex);
0174     bin.push_back(index);
0175     ++m_counter;
0176   }
0177   return gridIndex;
0178 }
0179 
0180 void CylindricalSpacePointGrid2::extend(
0181     const SpacePointContainer2::ConstRange& spacePoints) {
0182   ACTS_VERBOSE("Inserting " << spacePoints.size()
0183                             << " space points to the grid");
0184 
0185   for (const ConstSpacePointProxy2& sp : spacePoints) {
0186     insert(sp);
0187   }
0188 }
0189 
0190 void CylindricalSpacePointGrid2::sortBinsByR(
0191     const SpacePointContainer2& spacePoints) {
0192   ACTS_VERBOSE("Sorting the grid");
0193 
0194   for (std::size_t i = 0; i < grid().size(); ++i) {
0195     BinType& bin = grid().at(i);
0196     std::ranges::sort(bin, {}, [&](SpacePointIndex2 spIndex) {
0197       return spacePoints[spIndex].zr()[1];
0198     });
0199   }
0200 
0201   ACTS_VERBOSE(
0202       "Number of space points inserted (within grid range): " << m_counter);
0203 }
0204 
0205 Range1D<float> CylindricalSpacePointGrid2::computeRadiusRange(
0206     const SpacePointContainer2& spacePoints) const {
0207   float minRange = std::numeric_limits<float>::max();
0208   float maxRange = std::numeric_limits<float>::lowest();
0209   for (const BinType& bin : grid()) {
0210     if (bin.empty()) {
0211       continue;
0212     }
0213     auto first = spacePoints[bin.front()];
0214     auto last = spacePoints[bin.back()];
0215     minRange = std::min(first.zr()[1], minRange);
0216     maxRange = std::max(last.zr()[1], maxRange);
0217   }
0218   return {minRange, maxRange};
0219 }
0220 
0221 }  // namespace Acts