From 3d59079e716e3d250d7b8473f3b5898af4d3bb41 Mon Sep 17 00:00:00 2001 From: LuisFelipeCoelho Date: Mon, 13 Mar 2023 10:03:08 +0100 Subject: [PATCH] split SSS and PPP --- Core/include/Acts/Seeding/SeedFilter.hpp | 45 ++ Core/include/Acts/Seeding/SeedFilter.ipp | 428 +++++++++++ Core/include/Acts/Seeding/SeedFinder.hpp | 35 +- Core/include/Acts/Seeding/SeedFinder.ipp | 722 ++++++++++++++++-- .../TrackFinding/src/SeedingAlgorithm.cpp | 2 +- Examples/Python/python/acts/examples/itk.py | 1 + .../python/acts/examples/reconstruction.py | 4 +- 7 files changed, 1175 insertions(+), 62 deletions(-) diff --git a/Core/include/Acts/Seeding/SeedFilter.hpp b/Core/include/Acts/Seeding/SeedFilter.hpp index 276e2e93059..e5281b82f3c 100644 --- a/Core/include/Acts/Seeding/SeedFilter.hpp +++ b/Core/include/Acts/Seeding/SeedFilter.hpp @@ -68,6 +68,24 @@ class SeedFilter { CandidatesForMiddleSp>& candidates_collector) const; + virtual void filterSeeds_2SpFixedPPP( + InternalSpacePoint& bottomSP, + InternalSpacePoint& middleSP, + std::vector*>& topSpVec, + std::vector& invHelixDiameterVec, + std::vector& impactParametersVec, SeedFilterState& seedFilterState, + CandidatesForMiddleSp>& + candidates_collector) const; + + virtual void filterSeeds_2SpFixedSSS( + InternalSpacePoint& bottomSP, + InternalSpacePoint& middleSP, + std::vector*>& topSpVec, + std::vector& invHelixDiameterVec, + std::vector& impactParametersVec, SeedFilterState& seedFilterState, + CandidatesForMiddleSp>& + candidates_collector) const; + /// Filter seeds once all seeds for one middle space point have been created /// @param candidates_collector collection of seed candidates /// @param numQualitySeeds number of high quality seeds in seed confirmation @@ -80,6 +98,19 @@ class SeedFilter { std::back_insert_iterator>> outIt) const; + virtual void filterSeeds_1SpFixedPPP( + CandidatesForMiddleSp>& + candidates_collector, + std::size_t& numQualitySeeds, + std::back_insert_iterator>> outIt) + const; + + virtual void filterSeeds_1SpFixedSSS( + CandidatesForMiddleSp>& + candidates_collector, + std::back_insert_iterator>> outIt) + const; + /// Filter seeds once all seeds for one middle space point have been created /// @param candidates collection of seed candidates /// @param numQualitySeeds number of high quality seeds in seed confirmation @@ -92,6 +123,20 @@ class SeedFilter { std::back_insert_iterator>> outIt) const; + virtual void filterSeeds_1SpFixedPPP( + std::vector>::value_type>& candidates, + std::size_t& numQualitySeeds, + std::back_insert_iterator>> outIt) + const; + + virtual void filterSeeds_1SpFixedSSS( + std::vector>::value_type>& candidates, + std::back_insert_iterator>> outIt) + const; + + const SeedFilterConfig getSeedFilterConfig() const { return m_cfg; } const IExperimentCuts* getExperimentCuts() const { return m_experimentCuts; diff --git a/Core/include/Acts/Seeding/SeedFilter.ipp b/Core/include/Acts/Seeding/SeedFilter.ipp index 271a785c20d..d308c15963e 100644 --- a/Core/include/Acts/Seeding/SeedFilter.ipp +++ b/Core/include/Acts/Seeding/SeedFilter.ipp @@ -233,6 +233,320 @@ void SeedFilter::filterSeeds_2SpFixed( } } +// function to filter seeds based on all seeds with same bottom- and +// middle-spacepoint. +// return vector must contain weight of each seed +template +void SeedFilter::filterSeeds_2SpFixedPPP( + InternalSpacePoint& bottomSP, + InternalSpacePoint& middleSP, + std::vector*>& topSpVec, + std::vector& invHelixDiameterVec, + std::vector& impactParametersVec, SeedFilterState& seedFilterState, + CandidatesForMiddleSp>& + candidates_collector) const { + // seed confirmation + SeedConfirmationRangeConfig seedConfRange; + + // check if bottom SP is in the central or forward region + seedConfRange = + (bottomSP.z() > m_cfg.centralSeedConfirmationRange.zMaxSeedConf || + bottomSP.z() < m_cfg.centralSeedConfirmationRange.zMinSeedConf) + ? m_cfg.forwardSeedConfirmationRange + : m_cfg.centralSeedConfirmationRange; + // set the minimum number of top SP depending on whether the bottom SP is + // in the central or forward region + seedFilterState.nTopSeedConf = bottomSP.radius() > seedConfRange.rMaxSeedConf + ? seedConfRange.nTopForLargeR + : seedConfRange.nTopForSmallR; + + size_t maxWeightSeedIndex = 0; + bool maxWeightSeed = false; + float weightMax = -std::numeric_limits::max(); + float zOrigin = seedFilterState.zOrigin; + + // initialize original index locations + std::vector topSPIndexVec(topSpVec.size()); + for (std::size_t i(0); i < topSPIndexVec.size(); ++i) { + topSPIndexVec[i] = i; + } + + if (topSpVec.size() > 2) { + // sort indexes based on comparing values in invHelixDiameterVec + std::sort(topSPIndexVec.begin(), topSPIndexVec.end(), + [&invHelixDiameterVec](const size_t& i1, const size_t& i2) { + return invHelixDiameterVec[i1] < invHelixDiameterVec[i2]; + }); + } + + size_t beginCompTopIndex = 0; + // loop over top SPs and other compatible top SP candidates + for (auto& topSPIndex : topSPIndexVec) { + // if two compatible seeds with high distance in r are found, compatible + // seeds span 5 layers + // -> weaker requirement for a good seed + std::vector compatibleSeedR; + + float invHelixDiameter = invHelixDiameterVec[topSPIndex]; + float lowerLimitCurv = invHelixDiameter - m_cfg.deltaInvHelixDiameter; + float upperLimitCurv = invHelixDiameter + m_cfg.deltaInvHelixDiameter; + // use deltaR instead of top radius + float currentTopR = topSpVec[topSPIndex]->deltaR(); + float impact = impactParametersVec[topSPIndex]; + + float weight = -(impact * m_cfg.impactWeightFactor); + + // loop over compatible top SP candidates + for (size_t variableCompTopIndex = beginCompTopIndex; + variableCompTopIndex < topSPIndexVec.size(); variableCompTopIndex++) { + size_t compatibleTopSPIndex = topSPIndexVec[variableCompTopIndex]; + // continue if the two SPs are the same + if (compatibleTopSPIndex == topSPIndex) { + continue; + } + + float otherTopR = + topSpVec[compatibleTopSPIndex]->deltaR(); + + // curvature difference within limits? + if (invHelixDiameterVec[compatibleTopSPIndex] < lowerLimitCurv) { + // if SPs are sorted in curvature we skip unnecessary iterations + beginCompTopIndex = variableCompTopIndex + 1; + continue; + } + if (invHelixDiameterVec[compatibleTopSPIndex] > upperLimitCurv) { + // if SPs are sorted in curvature we skip unnecessary iterations + break; + } + // compared top SP should have at least deltaRMin distance + float deltaR = currentTopR - otherTopR; + if (std::abs(deltaR) < m_cfg.deltaRMin) { + continue; + } + bool newCompSeed = true; + for (float previousDiameter : compatibleSeedR) { + // original ATLAS code uses higher min distance for 2nd found compatible + // seed (20mm instead of 5mm) + // add new compatible seed only if distance larger than rmin to all + // other compatible seeds + if (std::abs(previousDiameter - otherTopR) < m_cfg.deltaRMin) { + newCompSeed = false; + break; + } + } + if (newCompSeed) { + compatibleSeedR.push_back(otherTopR); + weight += m_cfg.compatSeedWeight; + } + if (compatibleSeedR.size() >= m_cfg.compatSeedLimit) { + break; + } + } + + if (m_experimentCuts != nullptr) { + // add detector specific considerations on the seed weight + weight += m_experimentCuts->seedWeight(bottomSP, middleSP, + *topSpVec[topSPIndex]); + // discard seeds according to detector specific cuts (e.g.: weight) + if (!m_experimentCuts->singleSeedCut(weight, bottomSP, middleSP, + *topSpVec[topSPIndex])) { + continue; + } + } + + // increment in seed weight if number of compatible seeds is larger than + // numSeedIncrement + if (compatibleSeedR.size() > m_cfg.numSeedIncrement) { + weight += m_cfg.seedWeightIncrement; + } + + // seed confirmation cuts - keep seeds if they have specific values of + // impact parameter, z-origin and number of compatible seeds inside a + // pre-defined range that also depends on the region of the detector (i.e. + // forward or central region) defined by SeedConfirmationRange + int deltaSeedConf = + compatibleSeedR.size() + 1 - seedFilterState.nTopSeedConf; + if (deltaSeedConf < 0 || + (seedFilterState.numQualitySeeds != 0 and deltaSeedConf == 0)) { + continue; + } + bool seedRangeCuts = + bottomSP.radius() < seedConfRange.seedConfMinBottomRadius || + std::abs(zOrigin) > seedConfRange.seedConfMaxZOrigin; + if (seedRangeCuts and deltaSeedConf == 0 and + impact > seedConfRange.minImpactSeedConf) { + continue; + } + + // term on the weight that depends on the value of zOrigin + weight += -(std::abs(zOrigin) * m_cfg.zOriginWeightFactor) + + m_cfg.compatSeedWeight; + + // skip a bad quality seed if any of its constituents has a weight larger + // than the seed weight + if (weight < bottomSP.quality() and weight < middleSP.quality() and + weight < topSpVec[topSPIndex]->quality()) { + continue; + } + + if (deltaSeedConf > 0) { + // if we have not yet reached our max number of quality seeds we add the + // new seed to outCont + + // Internally, "push" will also checks the max number of quality seeds + // for a middle sp. + // If this is reached, we remove the seed with lowest weight. + candidates_collector.push(bottomSP, middleSP, *topSpVec[topSPIndex], + weight, zOrigin, true); + if (seedFilterState.numQualitySeeds < m_cfg.maxQualitySeedsPerSpMConf) { + // fill high quality seed + seedFilterState.numQualitySeeds++; + } + + } else if (weight > weightMax) { + // store weight and index of the best "lower quality" seed + weightMax = weight; + maxWeightSeedIndex = topSPIndex; + maxWeightSeed = true; + } + } // loop on tops + // if no high quality seed was found for a certain middle+bottom SP pair, + // lower quality seeds can be accepted + if (maxWeightSeed and seedFilterState.numQualitySeeds == 0) { + // if we have not yet reached our max number of seeds we add the new seed to + // outCont + + candidates_collector.push(bottomSP, middleSP, *topSpVec[maxWeightSeedIndex], + weightMax, zOrigin, false); + if (seedFilterState.numSeeds < m_cfg.maxSeedsPerSpMConf) { + // fill seed + seedFilterState.numSeeds++; + } + } +} // namespace Acts + + +template +void SeedFilter::filterSeeds_2SpFixedSSS( + InternalSpacePoint& bottomSP, + InternalSpacePoint& middleSP, + std::vector*>& topSpVec, + std::vector& invHelixDiameterVec, + std::vector& impactParametersVec, SeedFilterState& seedFilterState, + CandidatesForMiddleSp>& + candidates_collector) const { + float zOrigin = seedFilterState.zOrigin; + + // initialize original index locations + std::vector topSPIndexVec(topSpVec.size()); + for (std::size_t i(0); i < topSPIndexVec.size(); ++i) { + topSPIndexVec[i] = i; + } + + if (topSpVec.size() > 2) { + // sort indexes based on comparing values in invHelixDiameterVec + std::sort(topSPIndexVec.begin(), topSPIndexVec.end(), + [&invHelixDiameterVec](const size_t& i1, const size_t& i2) { + return invHelixDiameterVec[i1] < invHelixDiameterVec[i2]; + }); + } + + size_t beginCompTopIndex = 0; + // loop over top SPs and other compatible top SP candidates + for (auto& topSPIndex : topSPIndexVec) { + // if two compatible seeds with high distance in r are found, compatible + // seeds span 5 layers + // -> weaker requirement for a good seed + std::vector compatibleSeedR; + + float invHelixDiameter = invHelixDiameterVec[topSPIndex]; + float lowerLimitCurv = invHelixDiameter - m_cfg.deltaInvHelixDiameter; + float upperLimitCurv = invHelixDiameter + m_cfg.deltaInvHelixDiameter; + // use deltaR instead of top radius + float currentTopR = topSpVec[topSPIndex]->radius(); + float impact = impactParametersVec[topSPIndex]; + + float weight = -(impact * m_cfg.impactWeightFactor); + + // loop over compatible top SP candidates + for (size_t variableCompTopIndex = beginCompTopIndex; + variableCompTopIndex < topSPIndexVec.size(); variableCompTopIndex++) { + size_t compatibleTopSPIndex = topSPIndexVec[variableCompTopIndex]; + // continue if the two SPs are the same + if (compatibleTopSPIndex == topSPIndex) { + continue; + } + + float otherTopR = m_cfg.useDeltaRorTopRadius + ? topSpVec[compatibleTopSPIndex]->deltaR() + : topSpVec[compatibleTopSPIndex]->radius(); + + // curvature difference within limits? + if (invHelixDiameterVec[compatibleTopSPIndex] < lowerLimitCurv) { + // if SPs are sorted in curvature we skip unnecessary iterations + beginCompTopIndex = variableCompTopIndex + 1; + continue; + } + if (invHelixDiameterVec[compatibleTopSPIndex] > upperLimitCurv) { + // if SPs are sorted in curvature we skip unnecessary iterations + break; + } + // compared top SP should have at least deltaRMin distance + float deltaR = currentTopR - otherTopR; + if (std::abs(deltaR) < m_cfg.deltaRMin) { + continue; + } + bool newCompSeed = true; + for (float previousDiameter : compatibleSeedR) { + // original ATLAS code uses higher min distance for 2nd found compatible + // seed (20mm instead of 5mm) + // add new compatible seed only if distance larger than rmin to all + // other compatible seeds + if (std::abs(previousDiameter - otherTopR) < m_cfg.deltaRMin) { + newCompSeed = false; + break; + } + } + if (newCompSeed) { + compatibleSeedR.push_back(otherTopR); + weight += m_cfg.compatSeedWeight; + } + if (compatibleSeedR.size() >= m_cfg.compatSeedLimit) { + break; + } + } + + if (m_experimentCuts != nullptr) { + // add detector specific considerations on the seed weight + weight += m_experimentCuts->seedWeight(bottomSP, middleSP, + *topSpVec[topSPIndex]); + // discard seeds according to detector specific cuts (e.g.: weight) + if (!m_experimentCuts->singleSeedCut(weight, bottomSP, middleSP, + *topSpVec[topSPIndex])) { + continue; + } + } + + // increment in seed weight if number of compatible seeds is larger than + // numSeedIncrement + if (compatibleSeedR.size() > m_cfg.numSeedIncrement) { + weight += m_cfg.seedWeightIncrement; + } + + // keep the normal behavior without seed quality confirmation + // if we have not yet reached our max number of seeds we add the new seed + // to outCont + + candidates_collector.push(bottomSP, middleSP, *topSpVec[topSPIndex], weight, + zOrigin, false); + if (seedFilterState.numSeeds < m_cfg.maxSeedsPerSpMConf) { + // fill seed + seedFilterState.numSeeds++; + } + } +} + + // after creating all seeds with a common middle space point, filter again template @@ -249,6 +563,33 @@ void SeedFilter::filterSeeds_1SpFixed( filterSeeds_1SpFixed(extended_collection, numQualitySeeds, outIt); } +template +void SeedFilter::filterSeeds_1SpFixedPPP( + CandidatesForMiddleSp>& + candidates_collector, + std::size_t& numQualitySeeds, + std::back_insert_iterator>> outIt) +const { + // retrieve all candidates + // this collection is alredy sorted + // higher weights first + auto extended_collection = candidates_collector.storage(); + filterSeeds_1SpFixedPPP(extended_collection, numQualitySeeds, outIt); +} + +template +void SeedFilter::filterSeeds_1SpFixedSSS( + CandidatesForMiddleSp>& + candidates_collector, + std::back_insert_iterator>> outIt) +const { + // retrieve all candidates + // this collection is alredy sorted + // higher weights first + auto extended_collection = candidates_collector.storage(); + filterSeeds_1SpFixedSSS(extended_collection, outIt); +} + template void SeedFilter::filterSeeds_1SpFixed( std::vector::filterSeeds_1SpFixed( } } +template +void SeedFilter::filterSeeds_1SpFixedPPP( + std::vector>::value_type>& candidates, + std::size_t& numQualitySeeds, + std::back_insert_iterator>> outIt) +const { + if (m_experimentCuts != nullptr) { + candidates = m_experimentCuts->cutPerMiddleSP(std::move(candidates)); + } + + unsigned int maxSeeds = candidates.size(); + + if (maxSeeds > m_cfg.maxSeedsPerSpM) { + maxSeeds = m_cfg.maxSeedsPerSpM + 1; + } + + // default filter removes the last seeds if maximum amount exceeded + // ordering by weight by filterSeeds_2SpFixed means these are the lowest + // weight seeds + unsigned int numTotalSeeds = 0; + for (auto& [bottom, medium, top, bestSeedQuality, zOrigin, qualitySeed] : + candidates) { + // stop if we reach the maximum number of seeds + if (numTotalSeeds >= maxSeeds) { + break; + } + + // continue if higher-quality seeds were found + if (numQualitySeeds > 0 and not qualitySeed) { + continue; + } + if (bestSeedQuality < bottom->quality() and + bestSeedQuality < medium->quality() and + bestSeedQuality < top->quality()) { + continue; + } + + // set quality of seed components + bottom->setQuality(bestSeedQuality); + medium->setQuality(bestSeedQuality); + top->setQuality(bestSeedQuality); + + outIt = Seed{bottom->sp(), medium->sp(), top->sp(), + zOrigin, bestSeedQuality}; + ++numTotalSeeds; + } +} + +template +void SeedFilter::filterSeeds_1SpFixedSSS( + std::vector>::value_type>& candidates, + std::back_insert_iterator>> outIt) +const { + if (m_experimentCuts != nullptr) { + candidates = m_experimentCuts->cutPerMiddleSP(std::move(candidates)); + } + + unsigned int maxSeeds = candidates.size(); + + if (maxSeeds > m_cfg.maxSeedsPerSpM) { + maxSeeds = m_cfg.maxSeedsPerSpM + 1; + } + + // default filter removes the last seeds if maximum amount exceeded + // ordering by weight by filterSeeds_2SpFixed means these are the lowest + // weight seeds + unsigned int numTotalSeeds = 0; + for (auto& [bottom, medium, top, bestSeedQuality, zOrigin, qualitySeed] : + candidates) { + // stop if we reach the maximum number of seeds + if (numTotalSeeds >= maxSeeds) { + break; + } + + // set quality of seed components + bottom->setQuality(bestSeedQuality); + medium->setQuality(bestSeedQuality); + top->setQuality(bestSeedQuality); + + outIt = Seed{bottom->sp(), medium->sp(), top->sp(), + zOrigin, bestSeedQuality}; + ++numTotalSeeds; + } +} + } // namespace Acts diff --git a/Core/include/Acts/Seeding/SeedFinder.hpp b/Core/include/Acts/Seeding/SeedFinder.hpp index 79c9d4c2143..4d325382d72 100644 --- a/Core/include/Acts/Seeding/SeedFinder.hpp +++ b/Core/include/Acts/Seeding/SeedFinder.hpp @@ -86,6 +86,21 @@ class SeedFinder { sp_range_t& bottomSPs, sp_range_t& middleSPs, sp_range_t& topSPs, const Acts::Range1D& rMiddleSPRange) const; + template