Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2026-05-12 08:00:24

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