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 #pragma once
0010
0011 #include "Acts/Seeding/SeedFilter.hpp"
0012
0013 #include "Acts/EventData/Seed.hpp"
0014 #include "Acts/Seeding/detail/UtilityFunctions.hpp"
0015
0016 #include <algorithm>
0017 #include <utility>
0018
0019 namespace Acts {
0020
0021 template <typename external_spacepoint_t>
0022 SeedFilter<external_spacepoint_t>::SeedFilter(
0023     const SeedFilterConfig& config,
0024     IExperimentCuts<external_spacepoint_t>* expCuts /* = 0*/)
0025     : m_cfg(config), m_experimentCuts(expCuts) {}
0026
0027 template <typename external_spacepoint_t>
0028 SeedFilter<external_spacepoint_t>::SeedFilter(
0029     const SeedFilterConfig& config, std::unique_ptr<const Logger> logger,
0030     IExperimentCuts<external_spacepoint_t>* expCuts /* = 0*/)
0031     : m_cfg(config), m_logger(std::move(logger)), m_experimentCuts(expCuts) {}
0032
0033 template <typename external_spacepoint_t>
0034 void SeedFilter<external_spacepoint_t>::filterSeeds_2SpFixed(
0035     const SpacePointMutableData& mutableData,
0036     const external_spacepoint_t& bottomSp,
0037     const external_spacepoint_t& middleSp,
0038     const std::vector<const external_spacepoint_t*>& topSpVec,
0039     const std::vector<float>& invHelixDiameterVec,
0040     const std::vector<float>& impactParametersVec,
0041     SeedFilterState& seedFilterState,
0042     CandidatesForMiddleSp<const external_spacepoint_t>& candidatesCollector)
0043     const {
0044   // seed confirmation
0045   SeedConfirmationRangeConfig seedConfRange;
0046   if (m_cfg.seedConfirmation) {
0047     // check if bottom SP is in the central or forward region
0048     seedConfRange =
0049         (bottomSp.z() > m_cfg.centralSeedConfirmationRange.zMaxSeedConf ||
0050          bottomSp.z() < m_cfg.centralSeedConfirmationRange.zMinSeedConf)
0051             ? m_cfg.forwardSeedConfirmationRange
0052             : m_cfg.centralSeedConfirmationRange;
0053     // set the minimum number of top SP depending on whether the bottom SP is
0054     // in the central or forward region
0055     seedFilterState.nTopSeedConf =
0056         bottomSp.radius() > seedConfRange.rMaxSeedConf
0057             ? seedConfRange.nTopForLargeR
0058             : seedConfRange.nTopForSmallR;
0059   }
0060
0061   std::size_t maxWeightSeedIndex = 0;
0062   bool maxWeightSeed = false;
0063   float weightMax = std::numeric_limits<float>::lowest();
0064   float zOrigin = seedFilterState.zOrigin;
0065
0066   // initialize original index locations
0067   std::vector<std::size_t> topSpIndexVec(topSpVec.size());
0068   for (std::size_t i(0); i < topSpIndexVec.size(); ++i) {
0069     topSpIndexVec[i] = i;
0070   }
0071
0072   if (topSpVec.size() > 2) {
0073     // sort indexes based on comparing values in invHelixDiameterVec
0074     std::ranges::sort(topSpIndexVec, {},
0075                       [&invHelixDiameterVec](const std::size_t t) {
0076                         return invHelixDiameterVec[t];
0077                       });
0078   }
0079
0080   // vector containing the radius of all compatible seeds
0081   std::vector<float> compatibleSeedR;
0082   compatibleSeedR.reserve(m_cfg.compatSeedLimit);
0083
0084   std::size_t beginCompTopIndex = 0;
0085   // loop over top SPs and other compatible top SP candidates
0086   for (const std::size_t topSpIndex : topSpIndexVec) {
0087     // if two compatible seeds with high distance in r are found, compatible
0088     // seeds span 5 layers
0089     compatibleSeedR.clear();
0090
0091     float invHelixDiameter = invHelixDiameterVec[topSpIndex];
0092     float lowerLimitCurv = invHelixDiameter - m_cfg.deltaInvHelixDiameter;
0093     float upperLimitCurv = invHelixDiameter + m_cfg.deltaInvHelixDiameter;
0094     // use deltaR instead of top radius
0095     float currentTopR = m_cfg.useDeltaRorTopRadius
0096                             ? mutableData.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               ? mutableData.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           (candidatesCollector.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 < mutableData.quality(bottomSp.index()) &&
0194           weight < mutableData.quality(middleSp.index()) &&
0195           weight < mutableData.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         candidatesCollector.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       candidatesCollector.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       candidatesCollector.nHighQualityCandidates() == 0) {
0227     // if we have not yet reached our max number of seeds we add the new seed to
0228     // outCont
0229
0230     candidatesCollector.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     SpacePointMutableData& mutableData,
0241     CandidatesForMiddleSp<const external_spacepoint_t>& candidatesCollector,
0242     collection_t& outputCollection) const {
0243   // retrieve all candidates
0244   // this collection is already sorted
0245   // higher weights first
0246   std::size_t numQualitySeeds = candidatesCollector.nHighQualityCandidates();
0247   auto extended_collection = candidatesCollector.storage();
0248   filterSeeds_1SpFixed(mutableData, extended_collection, numQualitySeeds,
0249                        outputCollection);
0250 }
0251
0252 template <typename external_spacepoint_t>
0253 template <typename collection_t>
0254 void SeedFilter<external_spacepoint_t>::filterSeeds_1SpFixed(
0255     SpacePointMutableData& mutableData,
0256     std::vector<typename CandidatesForMiddleSp<
0257         const external_spacepoint_t>::value_type>& candidates,
0258     const std::size_t numQualitySeeds, collection_t& outputCollection) const {
0259   if (m_experimentCuts != nullptr) {
0260     candidates = m_experimentCuts->cutPerMiddleSP(std::move(candidates));
0261   }
0262
0263   unsigned int maxSeeds = candidates.size();
0264
0265   if (maxSeeds > m_cfg.maxSeedsPerSpM) {
0266     maxSeeds = m_cfg.maxSeedsPerSpM + 1;
0267   }
0268
0269   // default filter removes the last seeds if maximum amount exceeded
0270   // ordering by weight by filterSeeds_2SpFixed means these are the lowest
0271   // weight seeds
0272   unsigned int numTotalSeeds = 0;
0273   for (const auto& [bottom, medium, top, bestSeedQuality, zOrigin,
0274                     qualitySeed] : candidates) {
0275     // stop if we reach the maximum number of seeds
0276     if (numTotalSeeds >= maxSeeds) {
0277       break;
0278     }
0279
0280     if (m_cfg.seedConfirmation) {
0281       // continue if higher-quality seeds were found
0282       if (numQualitySeeds > 0 && !qualitySeed) {
0283         continue;
0284       }
0285       if (bestSeedQuality < mutableData.quality(bottom->index()) &&
0286           bestSeedQuality < mutableData.quality(medium->index()) &&
0287           bestSeedQuality < mutableData.quality(top->index())) {
0288         continue;
0289       }
0290     }
0291
0292     // set quality of seed components
0293     mutableData.setQuality(bottom->index(), bestSeedQuality);
0294     mutableData.setQuality(medium->index(), bestSeedQuality);
0295     mutableData.setQuality(top->index(), bestSeedQuality);
0296
0297     Seed<external_spacepoint_t> seed{*bottom, *medium, *top};
0298     seed.setVertexZ(zOrigin);
0299     seed.setQuality(bestSeedQuality);
0300
0301     ACTS_VERBOSE("Adding seed: [b=" << bottom->index() << ", m="
0302                                     << medium->index() << ", t=" << top->index()
0303                                     << "], quality=" << bestSeedQuality
0304                                     << ", vertexZ=" << zOrigin);
0305     detail::pushBackOrInsertAtEnd(outputCollection, std::move(seed));
0306     ++numTotalSeeds;
0307   }
0308   ACTS_VERBOSE("Identified " << numTotalSeeds << " seeds");
0309 }
0310
0311 }  // namespace Acts