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