Acts/Seeding/SeedFilter.ipp

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