Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-07-09 07:49:42

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/Seeding/detail/CylindricalSpacePointGrid.hpp"
0012 
0013 #include <concepts>
0014 #include <numbers>
0015 
0016 namespace Acts {
0017 
0018 template <typename external_spacepoint_t>
0019 CylindricalSpacePointGrid<external_spacepoint_t>
0020 CylindricalSpacePointGridCreator::createGrid(
0021     const CylindricalSpacePointGridConfig& config,
0022     const CylindricalSpacePointGridOptions& options, const Logger& logger) {
0023   config.checkConfig();
0024 
0025   using AxisScalar = Vector3::Scalar;
0026   using namespace UnitLiterals;
0027 
0028   int phiBins = 0;
0029   // for no magnetic field, create 100 phi-bins
0030   if (options.bFieldInZ == 0) {
0031     phiBins = 100;
0032     ACTS_VERBOSE(
0033         "B-Field is 0 (z-coordinate), setting the number of bins in phi to "
0034         << phiBins);
0035   } else {
0036     // calculate circle intersections of helix and max detector radius in mm.
0037     // bFieldInZ is in (pT/radius) natively, no need for conversion
0038     float minHelixRadius = config.minPt / options.bFieldInZ;
0039 
0040     // sanity check: if yOuter takes the square root of a negative number
0041     if (minHelixRadius < config.rMax * 0.5) {
0042       throw std::domain_error(
0043           "The value of minHelixRadius cannot be smaller than rMax / 2. Please "
0044           "check the configuration of bFieldInZ and minPt");
0045     }
0046 
0047     float maxR2 = config.rMax * config.rMax;
0048     float xOuter = maxR2 / (2 * minHelixRadius);
0049     float yOuter = std::sqrt(maxR2 - xOuter * xOuter);
0050     float outerAngle = std::atan(xOuter / yOuter);
0051     // intersection of helix and max detector radius minus maximum R distance
0052     // from middle SP to top SP
0053     float innerAngle = 0;
0054     float rMin = config.rMax;
0055     if (config.rMax > config.deltaRMax) {
0056       rMin = config.rMax - config.deltaRMax;
0057       float innerCircleR2 =
0058           (config.rMax - config.deltaRMax) * (config.rMax - config.deltaRMax);
0059       float xInner = innerCircleR2 / (2 * minHelixRadius);
0060       float yInner = std::sqrt(innerCircleR2 - xInner * xInner);
0061       innerAngle = std::atan(xInner / yInner);
0062     }
0063 
0064     // evaluating the azimutal deflection including the maximum impact parameter
0065     float deltaAngleWithMaxD0 =
0066         std::abs(std::asin(config.impactMax / (rMin)) -
0067                  std::asin(config.impactMax / config.rMax));
0068 
0069     // evaluating delta Phi based on the inner and outer angle, and the azimutal
0070     // deflection including the maximum impact parameter
0071     // Divide by config.phiBinDeflectionCoverage since we combine
0072     // config.phiBinDeflectionCoverage number of consecutive phi bins in the
0073     // seed making step. So each individual bin should cover
0074     // 1/config.phiBinDeflectionCoverage of the maximum expected azimutal
0075     // deflection
0076     float deltaPhi = (outerAngle - innerAngle + deltaAngleWithMaxD0) /
0077                      config.phiBinDeflectionCoverage;
0078 
0079     // sanity check: if the delta phi is equal to or less than zero, we'll be
0080     // creating an infinite or a negative number of bins, which would be bad!
0081     if (deltaPhi <= 0.f) {
0082       throw std::domain_error(
0083           "Delta phi value is equal to or less than zero, leading to an "
0084           "impossible number of bins (negative or infinite)");
0085     }
0086 
0087     // divide 2pi by angle delta to get number of phi-bins
0088     // size is always 2pi even for regions of interest
0089     phiBins = static_cast<int>(std::ceil(2 * std::numbers::pi / deltaPhi));
0090     // need to scale the number of phi bins accordingly to the number of
0091     // consecutive phi bins in the seed making step.
0092     // Each individual bin should be approximately a fraction (depending on this
0093     // number) of the maximum expected azimutal deflection.
0094 
0095     // set protection for large number of bins, by default it is large
0096     phiBins = std::min(phiBins, config.maxPhiBins);
0097   }
0098 
0099   Axis<AxisType::Equidistant, AxisBoundaryType::Closed> phiAxis(
0100       config.phiMin, config.phiMax, phiBins);
0101 
0102   // vector that will store the edges of the bins of z
0103   std::vector<AxisScalar> zValues{};
0104 
0105   // If zBinEdges is not defined, calculate the edges as zMin + bin * zBinSize
0106   if (config.zBinEdges.empty()) {
0107     // TODO: can probably be optimized using smaller z bins
0108     // and returning (multiple) neighbors only in one z-direction for forward
0109     // seeds
0110     // FIXME: zBinSize must include scattering
0111     float zBinSize = config.cotThetaMax * config.deltaRMax;
0112     float zBins =
0113         std::max(1.f, std::floor((config.zMax - config.zMin) / zBinSize));
0114 
0115     zValues.reserve(static_cast<int>(zBins));
0116     for (int bin = 0; bin <= static_cast<int>(zBins); bin++) {
0117       AxisScalar edge =
0118           config.zMin + bin * ((config.zMax - config.zMin) / zBins);
0119       zValues.push_back(edge);
0120     }
0121 
0122   } else {
0123     // Use the zBinEdges defined in the config
0124     zValues.reserve(config.zBinEdges.size());
0125     for (float bin : config.zBinEdges) {
0126       zValues.push_back(bin);
0127     }
0128   }
0129 
0130   std::vector<AxisScalar> rValues{};
0131   rValues.reserve(std::max(2ul, config.rBinEdges.size()));
0132   if (config.rBinEdges.empty()) {
0133     rValues = {config.rMin, config.rMax};
0134   } else {
0135     rValues.insert(rValues.end(), config.rBinEdges.begin(),
0136                    config.rBinEdges.end());
0137   }
0138 
0139   Axis<AxisType::Variable, AxisBoundaryType::Open> zAxis(std::move(zValues));
0140   Axis<AxisType::Variable, AxisBoundaryType::Open> rAxis(std::move(rValues));
0141 
0142   ACTS_VERBOSE("Defining Grid:");
0143   ACTS_VERBOSE("- Phi Axis: " << phiAxis);
0144   ACTS_VERBOSE("- Z axis  : " << zAxis);
0145   ACTS_VERBOSE("- R axis  : " << rAxis);
0146 
0147   return CylindricalSpacePointGrid<external_spacepoint_t>(
0148       std::make_tuple(std::move(phiAxis), std::move(zAxis), std::move(rAxis)));
0149 }
0150 
0151 template <typename external_spacepoint_t,
0152           typename external_spacepoint_iterator_t>
0153 void CylindricalSpacePointGridCreator::fillGrid(
0154     const SeedFinderConfig<external_spacepoint_t>& config,
0155     const SeedFinderOptions& options,
0156     CylindricalSpacePointGrid<external_spacepoint_t>& grid,
0157     external_spacepoint_iterator_t spBegin,
0158     external_spacepoint_iterator_t spEnd, const Logger& logger) {
0159   (void)options;
0160 
0161   if (config.seedFilter == nullptr) {
0162     throw std::runtime_error("SeedFinderConfig has a null SeedFilter object");
0163   }
0164 
0165   // Space points are assumed to be ALREADY CORRECTED for beamspot position
0166   // phi, z and r space point selection comes naturally from the
0167   // grid axis definition. Calling `isInside` will let us know if we are
0168   // inside the grid range.
0169   // If a space point is outside the validity range of these quantities
0170   // it goes in an over- or under-flow bin. We want to avoid to consider those
0171   // and skip some computations.
0172   // Additional cuts can be applied by customizing the space point selector
0173   // in the config object.
0174 
0175   // keep track of changed bins while sorting
0176   std::vector<bool> usedBinIndex(grid.size(), false);
0177   std::vector<std::size_t> rBinsIndex;
0178   rBinsIndex.reserve(grid.size());
0179 
0180   ACTS_VERBOSE("Fetching " << std::distance(spBegin, spEnd)
0181                            << " space points to the grid");
0182   std::size_t counter = 0ul;
0183   for (external_spacepoint_iterator_t it = spBegin; it != spEnd; ++it) {
0184     const external_spacepoint_t& sp = *it;
0185 
0186     // remove SPs according to experiment specific cuts
0187     if (!config.spacePointSelector(sp)) {
0188       continue;
0189     }
0190 
0191     // fill rbins into grid
0192     Vector3 position(sp.phi(), sp.z(), sp.radius());
0193     if (!grid.isInside(position)) {
0194       continue;
0195     }
0196 
0197     std::size_t globIndex = grid.globalBinFromPosition(position);
0198     auto& rbin = grid.at(globIndex);
0199     rbin.push_back(&sp);
0200     ++counter;
0201 
0202     // keep track of the bins we modify so that we can later sort the SPs in
0203     // those bins only
0204     if (rbin.size() > 1 && !usedBinIndex[globIndex]) {
0205       usedBinIndex[globIndex] = true;
0206       rBinsIndex.push_back(globIndex);
0207     }
0208   }
0209 
0210   /// sort SPs in R for each filled bin
0211   for (std::size_t binIndex : rBinsIndex) {
0212     auto& rbin = grid.atPosition(binIndex);
0213     std::ranges::sort(rbin, {}, [](const auto& rb) { return rb->radius(); });
0214   }
0215 
0216   ACTS_VERBOSE(
0217       "Number of space points inserted (within grid range): " << counter);
0218 }
0219 
0220 template <typename external_spacepoint_t, typename external_collection_t>
0221   requires std::ranges::range<external_collection_t> &&
0222            std::same_as<typename external_collection_t::value_type,
0223                         external_spacepoint_t>
0224 void CylindricalSpacePointGridCreator::fillGrid(
0225     const SeedFinderConfig<external_spacepoint_t>& config,
0226     const SeedFinderOptions& options,
0227     CylindricalSpacePointGrid<external_spacepoint_t>& grid,
0228     const external_collection_t& collection, const Logger& logger) {
0229   CylindricalSpacePointGridCreator::fillGrid<external_spacepoint_t>(
0230       config, options, grid, std::ranges::begin(collection),
0231       std::ranges::end(collection), logger);
0232 }
0233 
0234 }  // namespace Acts