Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-07-05 08:26:27

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2024 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 http://mozilla.org/MPL/2.0/.
0008 
0009 template <typename SpacePoint>
0010 Acts::CylindricalSpacePointGrid<SpacePoint>
0011 Acts::CylindricalSpacePointGridCreator::createGrid(
0012     const Acts::CylindricalSpacePointGridConfig& config,
0013     const Acts::CylindricalSpacePointGridOptions& options) {
0014   if (!config.isInInternalUnits) {
0015     throw std::runtime_error(
0016         "CylindricalSpacePointGridConfig not in ACTS internal units in "
0017         "CylindricalSpacePointGridCreator::createGrid");
0018   }
0019   if (!options.isInInternalUnits) {
0020     throw std::runtime_error(
0021         "CylindricalSpacePointGridOptions not in ACTS internal units in "
0022         "CylindricalSpacePointGridCreator::createGrid");
0023   }
0024   using AxisScalar = Acts::Vector3::Scalar;
0025   using namespace Acts::UnitLiterals;
0026 
0027   int phiBins = 0;
0028   // for no magnetic field, create 100 phi-bins
0029   if (options.bFieldInZ == 0) {
0030     phiBins = 100;
0031   } else {
0032     // calculate circle intersections of helix and max detector radius
0033     float minHelixRadius =
0034         config.minPt /
0035         (1_T * 1e6 *
0036          options.bFieldInZ);  // in mm -> R[mm] =pT[GeV] / (3·10−4×B[T])
0037                               // = pT[MeV] / (300 *Bz[kT])
0038 
0039     // sanity check: if yOuter takes the square root of a negative number
0040     if (minHelixRadius < config.rMax / 2) {
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 * M_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     if (phiBins > config.maxPhiBins) {
0096       phiBins = config.maxPhiBins;
0097     }
0098   }
0099 
0100   Acts::Axis<AxisType::Equidistant, AxisBoundaryType::Closed> phiAxis(
0101       config.phiMin, config.phiMax, phiBins);
0102 
0103   // vector that will store the edges of the bins of z
0104   std::vector<AxisScalar> zValues;
0105 
0106   // If zBinEdges is not defined, calculate the edges as zMin + bin * zBinSize
0107   if (config.zBinEdges.empty()) {
0108     // TODO: can probably be optimized using smaller z bins
0109     // and returning (multiple) neighbors only in one z-direction for forward
0110     // seeds
0111     // FIXME: zBinSize must include scattering
0112     float zBinSize = config.cotThetaMax * config.deltaRMax;
0113     float zBins =
0114         std::max(1.f, std::floor((config.zMax - config.zMin) / zBinSize));
0115 
0116     zValues.reserve(static_cast<int>(zBins));
0117     for (int bin = 0; bin <= static_cast<int>(zBins); bin++) {
0118       AxisScalar edge =
0119           config.zMin + bin * ((config.zMax - config.zMin) / zBins);
0120       zValues.push_back(edge);
0121     }
0122 
0123   } else {
0124     // Use the zBinEdges defined in the config
0125     zValues.reserve(config.zBinEdges.size());
0126     for (float bin : config.zBinEdges) {
0127       zValues.push_back(bin);
0128     }
0129   }
0130 
0131   Axis<AxisType::Variable, AxisBoundaryType::Bound> zAxis(std::move(zValues));
0132   return Acts::CylindricalSpacePointGrid<SpacePoint>(
0133       std::make_tuple(std::move(phiAxis), std::move(zAxis)));
0134 }
0135 
0136 template <typename external_spacepoint_t,
0137           typename external_spacepoint_iterator_t, typename callable_t>
0138 void Acts::CylindricalSpacePointGridCreator::fillGrid(
0139     const Acts::SeedFinderConfig<external_spacepoint_t>& config,
0140     const Acts::SeedFinderOptions& options,
0141     Acts::CylindricalSpacePointGrid<external_spacepoint_t>& grid,
0142     external_spacepoint_iterator_t spBegin,
0143     external_spacepoint_iterator_t spEnd, callable_t&& toGlobal,
0144     Acts::Extent& rRangeSPExtent) {
0145   using iterated_value_t =
0146       typename std::iterator_traits<external_spacepoint_iterator_t>::value_type;
0147   using iterated_t = typename std::remove_const<
0148       typename std::remove_pointer<iterated_value_t>::type>::type;
0149   static_assert(std::is_pointer<iterated_value_t>::value,
0150                 "Iterator must contain pointers to space points");
0151   static_assert(std::is_same<iterated_t, external_spacepoint_t>::value,
0152                 "Iterator does not contain type this class was templated with");
0153 
0154   if (!config.isInInternalUnits) {
0155     throw std::runtime_error(
0156         "SeedFinderConfig not in ACTS internal units in BinnedSPGroup");
0157   }
0158   if (config.seedFilter == nullptr) {
0159     throw std::runtime_error("SeedFinderConfig has a null SeedFilter object");
0160   }
0161   if (!options.isInInternalUnits) {
0162     throw std::runtime_error(
0163         "SeedFinderOptions not in ACTS internal units in BinnedSPGroup");
0164   }
0165 
0166   // sort by radius
0167   // add magnitude of beamPos to rMax to avoid excluding measurements
0168   // create number of bins equal to number of millimeters rMax
0169   // (worst case minR: configured minR + 1mm)
0170   // binSizeR allows to increase or reduce numRBins if needed
0171   std::size_t numRBins = static_cast<std::size_t>(
0172       (config.rMax + options.beamPos.norm()) / config.binSizeR);
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   std::size_t counter = 0ul;
0180   for (external_spacepoint_iterator_t it = spBegin; it != spEnd;
0181        it++, ++counter) {
0182     if (*it == nullptr) {
0183       continue;
0184     }
0185     const external_spacepoint_t& sp = **it;
0186     const auto& [spPosition, variance, spTime] =
0187         toGlobal(sp, config.zAlign, config.rAlign, config.sigmaError);
0188 
0189     float spX = spPosition[0];
0190     float spY = spPosition[1];
0191     float spZ = spPosition[2];
0192 
0193     // store x,y,z values in extent
0194     rRangeSPExtent.extend({spX, spY, spZ});
0195 
0196     // remove SPs according to experiment specific cuts
0197     if (!config.spacePointSelector(sp)) {
0198       continue;
0199     }
0200 
0201     // remove SPs outside z and phi region
0202     if (spZ > config.zMax || spZ < config.zMin) {
0203       continue;
0204     }
0205 
0206     float spPhi = std::atan2(spY, spX);
0207     if (spPhi > config.phiMax || spPhi < config.phiMin) {
0208       continue;
0209     }
0210 
0211     auto isp = std::make_unique<InternalSpacePoint<external_spacepoint_t>>(
0212         counter, sp, spPosition, options.beamPos, variance, spTime);
0213     // calculate r-Bin index and protect against overflow (underflow not
0214     // possible)
0215     std::size_t rIndex =
0216         static_cast<std::size_t>(isp->radius() / config.binSizeR);
0217     // if index out of bounds, the SP is outside the region of interest
0218     if (rIndex >= numRBins) {
0219       continue;
0220     }
0221 
0222     // fill rbins into grid
0223     std::size_t globIndex =
0224         grid.globalBinFromPosition(Acts::Vector2{isp->phi(), isp->z()});
0225     auto& rbin = grid.at(globIndex);
0226     rbin.push_back(std::move(isp));
0227 
0228     // keep track of the bins we modify so that we can later sort the SPs in
0229     // those bins only
0230     if (rbin.size() > 1 && !usedBinIndex[globIndex]) {
0231       usedBinIndex[globIndex] = true;
0232       rBinsIndex.push_back(globIndex);
0233     }
0234   }
0235 
0236   /// sort SPs in R for each filled bin
0237   for (std::size_t binIndex : rBinsIndex) {
0238     auto& rbin = grid.atPosition(binIndex);
0239     std::sort(rbin.begin(), rbin.end(),
0240               [](const auto& a, const auto& b) -> bool {
0241                 return a->radius() < b->radius();
0242               });
0243   }
0244 }