Acts/Seeding/SeedFilter.ipp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2023 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 #include "Acts/Seeding/detail/UtilityFunctions.hpp"
0010
0011 #include <algorithm>
0012 #include <numeric>
0013 #include <utility>
0014
0015 namespace Acts {
0016 // constructor
0017 template <typename external_spacepoint_t>
0018 SeedFilter<external_spacepoint_t>::SeedFilter(
0019     SeedFilterConfig config,
0020     IExperimentCuts<external_spacepoint_t>* expCuts /* = 0*/)
0021     : m_cfg(config), m_experimentCuts(expCuts) {
0022   if (!config.isInInternalUnits) {
0023     throw std::runtime_error(
0024         "SeedFilterConfig not in ACTS internal units in SeedFilter");
0025   }
0026 }
0027 // function to filter seeds based on all seeds with same bottom- and
0028 // middle-spacepoint.
0029 // return vector must contain weight of each seed
0030 template <typename external_spacepoint_t>
0031 void SeedFilter<external_spacepoint_t>::filterSeeds_2SpFixed(
0032     Acts::SpacePointData& spacePointData,
0033     const InternalSpacePoint<external_spacepoint_t>& bottomSP,
0034     const InternalSpacePoint<external_spacepoint_t>& middleSP,
0035     const std::vector<const InternalSpacePoint<external_spacepoint_t>*>&
0036         topSpVec,
0037     const std::vector<float>& invHelixDiameterVec,
0038     const std::vector<float>& impactParametersVec,
0039     SeedFilterState& seedFilterState,
0040     CandidatesForMiddleSp<const InternalSpacePoint<external_spacepoint_t>>&
0041         candidates_collector) const {
0042   // seed confirmation
0043   SeedConfirmationRangeConfig seedConfRange;
0044   if (m_cfg.seedConfirmation) {
0045     // check if bottom SP is in the central or forward region
0046     seedConfRange =
0047         (bottomSP.z() > m_cfg.centralSeedConfirmationRange.zMaxSeedConf ||
0048          bottomSP.z() < m_cfg.centralSeedConfirmationRange.zMinSeedConf)
0049             ? m_cfg.forwardSeedConfirmationRange
0050             : m_cfg.centralSeedConfirmationRange;
0051     // set the minimum number of top SP depending on whether the bottom SP is
0052     // in the central or forward region
0053     seedFilterState.nTopSeedConf =
0054         bottomSP.radius() > seedConfRange.rMaxSeedConf
0055             ? seedConfRange.nTopForLargeR
0056             : seedConfRange.nTopForSmallR;
0057   }
0058
0059   std::size_t maxWeightSeedIndex = 0;
0060   bool maxWeightSeed = false;
0061   float weightMax = std::numeric_limits<float>::lowest();
0062   float zOrigin = seedFilterState.zOrigin;
0063
0064   // initialize original index locations
0065   std::vector<std::size_t> topSPIndexVec(topSpVec.size());
0066   for (std::size_t i(0); i < topSPIndexVec.size(); ++i) {
0067     topSPIndexVec[i] = i;
0068   }
0069
0070   if (topSpVec.size() > 2) {
0071     // sort indexes based on comparing values in invHelixDiameterVec
0072     std::sort(
0073         topSPIndexVec.begin(), topSPIndexVec.end(),
0074         [&invHelixDiameterVec](const std::size_t i1, const std::size_t i2) {
0075           return invHelixDiameterVec[i1] < invHelixDiameterVec[i2];
0076         });
0077   }
0078
0079   // vector containing the radius of all compatible seeds
0080   std::vector<float> compatibleSeedR;
0081   compatibleSeedR.reserve(m_cfg.compatSeedLimit);
0082
0083   std::size_t beginCompTopIndex = 0;
0084   // loop over top SPs and other compatible top SP candidates
0085   for (const std::size_t topSPIndex : topSPIndexVec) {
0086     // if two compatible seeds with high distance in r are found, compatible
0087     // seeds span 5 layers
0088     compatibleSeedR.clear();
0089
0090     float invHelixDiameter = invHelixDiameterVec[topSPIndex];
0091     float lowerLimitCurv = invHelixDiameter - m_cfg.deltaInvHelixDiameter;
0092     float upperLimitCurv = invHelixDiameter + m_cfg.deltaInvHelixDiameter;
0093     // use deltaR instead of top radius
0094     float currentTopR =
0095         m_cfg.useDeltaRorTopRadius
0096             ? spacePointData.deltaR(topSpVec[topSPIndex]->index())
0097             : topSpVec[topSPIndex]->radius();
0098     float impact = impactParametersVec[topSPIndex];
0099
0100     float weight = -(impact * m_cfg.impactWeightFactor);
0101
0102     // loop over compatible top SP candidates
0103     for (std::size_t variableCompTopIndex = beginCompTopIndex;
0104          variableCompTopIndex < topSPIndexVec.size(); variableCompTopIndex++) {
0105       std::size_t compatibleTopSPIndex = topSPIndexVec[variableCompTopIndex];
0106       if (compatibleTopSPIndex == topSPIndex) {
0107         continue;
0108       }
0109
0110       float otherTopR =
0111           m_cfg.useDeltaRorTopRadius
0112               ? spacePointData.deltaR(topSpVec[compatibleTopSPIndex]->index())
0113               : topSpVec[compatibleTopSPIndex]->radius();
0114
0115       // curvature difference within limits?
0116       if (invHelixDiameterVec[compatibleTopSPIndex] < lowerLimitCurv) {
0117         // the SPs are sorted in curvature so we skip unnecessary iterations
0118         beginCompTopIndex = variableCompTopIndex + 1;
0119         continue;
0120       }
0121       if (invHelixDiameterVec[compatibleTopSPIndex] > upperLimitCurv) {
0122         // the SPs are sorted in curvature so we skip unnecessary iterations
0123         break;
0124       }
0125       // compared top SP should have at least deltaRMin distance
0126       float deltaR = currentTopR - otherTopR;
0127       if (std::abs(deltaR) < m_cfg.deltaRMin) {
0128         continue;
0129       }
0130       bool newCompSeed = true;
0131       for (const float previousDiameter : compatibleSeedR) {
0132         // original ATLAS code uses higher min distance for 2nd found compatible
0133         // seed (20mm instead of 5mm)
0134         // add new compatible seed only if distance larger than rmin to all
0135         // other compatible seeds
0136         if (std::abs(previousDiameter - otherTopR) < m_cfg.deltaRMin) {
0137           newCompSeed = false;
0138           break;
0139         }
0140       }
0141       if (newCompSeed) {
0142         compatibleSeedR.push_back(otherTopR);
0143         weight += m_cfg.compatSeedWeight;
0144       }
0145       if (compatibleSeedR.size() >= m_cfg.compatSeedLimit) {
0146         break;
0147       }
0148     }
0149
0150     if (m_experimentCuts != nullptr) {
0151       // add detector specific considerations on the seed weight
0152       weight += m_experimentCuts->seedWeight(bottomSP, middleSP,
0153                                              *topSpVec[topSPIndex]);
0154       // discard seeds according to detector specific cuts (e.g.: weight)
0155       if (!m_experimentCuts->singleSeedCut(weight, bottomSP, middleSP,
0156                                            *topSpVec[topSPIndex])) {
0157         continue;
0158       }
0159     }
0160
0161     // increment in seed weight if number of compatible seeds is larger than
0162     // numSeedIncrement
0163     if (compatibleSeedR.size() > m_cfg.numSeedIncrement) {
0164       weight += m_cfg.seedWeightIncrement;
0165     }
0166
0167     if (m_cfg.seedConfirmation) {
0168       // seed confirmation cuts - keep seeds if they have specific values of
0169       // impact parameter, z-origin and number of compatible seeds inside a
0170       // pre-defined range that also depends on the region of the detector (i.e.
0171       // forward or central region) defined by SeedConfirmationRange
0172       int deltaSeedConf =
0173           compatibleSeedR.size() + 1 - seedFilterState.nTopSeedConf;
0174       if (deltaSeedConf < 0 ||
0175           (candidates_collector.nHighQualityCandidates() != 0 &&
0176            deltaSeedConf == 0)) {
0177         continue;
0178       }
0179       bool seedRangeCuts =
0180           bottomSP.radius() < seedConfRange.seedConfMinBottomRadius ||
0181           std::abs(zOrigin) > seedConfRange.seedConfMaxZOrigin;
0182       if (seedRangeCuts && deltaSeedConf == 0 &&
0183           impact > seedConfRange.minImpactSeedConf) {
0184         continue;
0185       }
0186
0187       // term on the weight that depends on the value of zOrigin
0188       weight += -(std::abs(zOrigin) * m_cfg.zOriginWeightFactor) +
0189                 m_cfg.compatSeedWeight;
0190
0191       // skip a bad quality seed if any of its constituents has a weight larger
0192       // than the seed weight
0193       if (weight < spacePointData.quality(bottomSP.index()) &&
0194           weight < spacePointData.quality(middleSP.index()) &&
0195           weight < spacePointData.quality(topSpVec[topSPIndex]->index())) {
0196         continue;
0197       }
0198
0199       if (deltaSeedConf > 0) {
0200         // if we have not yet reached our max number of quality seeds we add the
0201         // new seed to outCont
0202
0203         // Internally, "push" will also check the max number of quality seeds
0204         // for a middle sp.
0205         // If this is reached, we remove the seed with the lowest weight.
0206         candidates_collector.push(bottomSP, middleSP, *topSpVec[topSPIndex],
0207                                   weight, zOrigin, true);
0208       } else if (weight > weightMax) {
0209         // store weight and index of the best "lower quality" seed
0210         weightMax = weight;
0211         maxWeightSeedIndex = topSPIndex;
0212         maxWeightSeed = true;
0213       }
0214     } else {
0215       // keep the normal behavior without seed quality confirmation
0216       // if we have not yet reached our max number of seeds we add the new seed
0217       // to outCont
0218
0219       candidates_collector.push(bottomSP, middleSP, *topSpVec[topSPIndex],
0220                                 weight, zOrigin, false);
0221     }
0222   }  // loop on tops
0223   // if no high quality seed was found for a certain middle+bottom SP pair,
0224   // lower quality seeds can be accepted
0225   if (m_cfg.seedConfirmation && maxWeightSeed &&
0226       candidates_collector.nHighQualityCandidates() == 0) {
0227     // if we have not yet reached our max number of seeds we add the new seed to
0228     // outCont
0229
0230     candidates_collector.push(bottomSP, middleSP, *topSpVec[maxWeightSeedIndex],
0231                               weightMax, zOrigin, false);
0232   }
0233 }
0234
0235 // after creating all seeds with a common middle space point, filter again
0236
0237 template <typename external_spacepoint_t>
0238 template <typename collection_t>
0239 void SeedFilter<external_spacepoint_t>::filterSeeds_1SpFixed(
0240     Acts::SpacePointData& spacePointData,
0241     CandidatesForMiddleSp<const InternalSpacePoint<external_spacepoint_t>>&
0242         candidates_collector,
0243     collection_t& outputCollection) const {
0244   // retrieve all candidates
0245   // this collection is already sorted
0246   // higher weights first
0247   std::size_t numQualitySeeds = candidates_collector.nHighQualityCandidates();
0248   auto extended_collection = candidates_collector.storage();
0249   filterSeeds_1SpFixed(spacePointData, extended_collection, numQualitySeeds,
0250                        outputCollection);
0251 }
0252
0253 template <typename external_spacepoint_t>
0254 template <typename collection_t>
0255 void SeedFilter<external_spacepoint_t>::filterSeeds_1SpFixed(
0256     Acts::SpacePointData& spacePointData,
0257     std::vector<typename CandidatesForMiddleSp<
0258         const InternalSpacePoint<external_spacepoint_t>>::value_type>&
0259         candidates,
0260     const std::size_t numQualitySeeds, collection_t& outputCollection) const {
0261   if (m_experimentCuts != nullptr) {
0262     candidates = m_experimentCuts->cutPerMiddleSP(std::move(candidates));
0263   }
0264
0265   unsigned int maxSeeds = candidates.size();
0266
0267   if (maxSeeds > m_cfg.maxSeedsPerSpM) {
0268     maxSeeds = m_cfg.maxSeedsPerSpM + 1;
0269   }
0270
0271   // default filter removes the last seeds if maximum amount exceeded
0272   // ordering by weight by filterSeeds_2SpFixed means these are the lowest
0273   // weight seeds
0274   unsigned int numTotalSeeds = 0;
0275   for (const auto& [bottom, medium, top, bestSeedQuality, zOrigin,
0276                     qualitySeed] : candidates) {
0277     // stop if we reach the maximum number of seeds
0278     if (numTotalSeeds >= maxSeeds) {
0279       break;
0280     }
0281
0282     if (m_cfg.seedConfirmation) {
0283       // continue if higher-quality seeds were found
0284       if (numQualitySeeds > 0 && !qualitySeed) {
0285         continue;
0286       }
0287       if (bestSeedQuality < spacePointData.quality(bottom->index()) &&
0288           bestSeedQuality < spacePointData.quality(medium->index()) &&
0289           bestSeedQuality < spacePointData.quality(top->index())) {
0290         continue;
0291       }
0292     }
0293
0294     // set quality of seed components
0295     spacePointData.setQuality(bottom->index(), bestSeedQuality);
0296     spacePointData.setQuality(medium->index(), bestSeedQuality);
0297     spacePointData.setQuality(top->index(), bestSeedQuality);
0298
0299     Acts::detail::pushBackOrInsertAtEnd(
0300         outputCollection,
0301         Seed<external_spacepoint_t>{bottom->sp(), medium->sp(), top->sp(),
0302                                     zOrigin, bestSeedQuality});
0303
0304     ++numTotalSeeds;
0305   }
0306 }
0307
0308 }  // namespace Acts