19 #include "Acts/Seeding/CombinatorialSeedSolver.hpp"
24 using namespace Acts::Experimental::CombinatorialSeedSolver;
32 const auto*
re = prd->readoutElement();
33 switch(prd->channelType()) {
35 return re->stripDesign(prd->measurementHash());
37 return re->wireDesign(prd->measurementHash());
39 return re->padDesign(prd->measurementHash());
45 const auto& design = getDesign(sp);
48 return 0.5* design.stripLength(prd->channelNumber());
56 using namespace SegmentFit;
64 ATH_CHECK(m_visionTool.retrieve(DisableTool{m_visionTool.empty()}));
66 if (!(m_idHelperSvc->hasMM() || m_idHelperSvc->hasSTGC())) {
67 ATH_MSG_ERROR(
"MM or STGC not part of initialized detector layout");
68 return StatusCode::FAILURE;
71 return StatusCode::SUCCESS;
75 CombinatorialNSWSeedFinderAlg::emptyBookKeeper(
const HitLayVec& sortedSp)
const{
77 for (std::size_t
l = 0;
l < sortedSp.size(); ++
l) {
78 emptyKeeper[
l].resize(sortedSp[
l].
size(), 0);
84 CombinatorialNSWSeedFinderAlg::classifyStrip(
const SpacePoint& sp)
const{
87 const auto& design = getDesign(sp);
88 if (!design.hasStereoAngle()) {
91 return design.stereoAngle() > 0. ? StripOrient::U : StripOrient::V;
104 CombinatorialNSWSeedFinderAlg::hitFromIPCorridor(
const SpacePoint& testHit,
109 const Amg::Vector3D estPlaneArrivalUp = SeedingAux::extrapolateToPlane(beamSpotPos, dirEstUp, testHit);
110 const Amg::Vector3D estPlaneArrivalDn = SeedingAux::extrapolateToPlane(beamSpotPos, dirEstDn, testHit);
112 bool below{
true}, above{
true};
113 switch (classifyStrip(testHit)) {
117 const double halfLength = 0.5* stripHalfLength(testHit);
123 below = estPlaneArrivalDn.y() >
std::max(leftEdge.y(), rightEdge.y());
125 above = estPlaneArrivalUp.y() <
std::min(leftEdge.y(), rightEdge.y());
131 below = estPlaneArrivalDn.y() > hY;
133 above = estPlaneArrivalUp.y() < hY;
145 << (below || above ?
" is outside the window" :
" is inside the window"));
147 return HitWindow::tooLow;
150 return HitWindow::tooHigh;
156 #define TEST_HIT_CORRIDOR(LAYER, HIT_ITER, START_LAYER) \
158 const SpacePoint* testMe = combinatoricLayers[LAYER].get()[HIT_ITER]; \
159 if (usedHits[LAYER].get()[HIT_ITER]) { \
160 ATH_MSG_VERBOSE(__func__<<":"<<__LINE__<<" - " \
161 <<m_idHelperSvc->toString(testMe->identify()) \
162 <<" already used in good seed." ); \
165 const HitWindow inWindow = hitFromIPCorridor(*testMe, beamSpot, dirEstUp, dirEstDn); \
166 if(inWindow == HitWindow::tooHigh) { \
167 ATH_MSG_VERBOSE(__func__<<":"<<__LINE__<<" - Hit " \
168 <<m_idHelperSvc->toString(testMe->identify()) \
169 <<" is beyond the corridor. Break loop"); \
171 } else if (inWindow == HitWindow::tooLow) { \
172 START_LAYER = HIT_ITER + 1; \
173 ATH_MSG_VERBOSE(__func__<<":"<<__LINE__<<" - Hit " \
174 <<m_idHelperSvc->toString(testMe->identify()) \
175 <<" is still below the corridor. Update start to " \
186 seedHitsFromLayers.clear();
187 std::size_t maxSize{1};
188 for (
const HitVec& hitVec : combinatoricLayers) {
189 maxSize = maxSize * hitVec.size();
191 seedHitsFromLayers.reserve(maxSize);
194 unsigned int iterLay0{0}, iterLay1{0}, iterLay2{0}, iterLay3{0};
195 unsigned int startLay1{0}, startLay2{0}, startLay3{0};
197 for( ; iterLay0 < combinatoricLayers[0].get().
size() ; ++iterLay0){
199 if (usedHits[0].
get()[iterLay0]) {
202 const SpacePoint* hit0 = combinatoricLayers[0].get()[iterLay0];
213 for( iterLay1 = startLay1; iterLay1 < combinatoricLayers[1].get().
size() ; ++iterLay1){
215 for( iterLay2 = startLay2; iterLay2 < combinatoricLayers[2].get().
size() ; ++iterLay2){
217 for( iterLay3 = startLay3; iterLay3 < combinatoricLayers[3].get().
size(); ++iterLay3){
219 seedHitsFromLayers.emplace_back(
std::array{hit0, combinatoricLayers[1].get()[iterLay1],
220 combinatoricLayers[2].
get()[iterLay2],
221 combinatoricLayers[3].
get()[iterLay3]});
227 #undef TEST_HIT_CORRIDOR
239 for (std::size_t
i = 0;
i < extensionLayers.size(); ++
i) {
241 const Amg::Vector3D extrapPos = SeedingAux::extrapolateToPlane(startPos, direction, *
layer.front());
243 unsigned int indexOfHit =
layer.size() + 1;
244 unsigned int triedHit{0};
248 for (
unsigned int j = 0; j <
layer.size(); ++j) {
249 if (usedHits[
i].
get().at(j)) {
252 auto hit =
layer.at(j);
253 const double pull = std::sqrt(SeedingAux::chi2Term(extrapPos, direction, *hit));
254 ATH_MSG_VERBOSE(
"Trying extension with hit " << m_idHelperSvc->toString(hit->identify()));
257 if (
pull > minPull) {
271 if (minPull < m_minPullThreshold) {
272 const auto* bestCand =
layer.at(indexOfHit);
273 ATH_MSG_VERBOSE(
"Extension successfull - hit" << m_idHelperSvc->toString(bestCand->identify())
275 <<
", dir: "<<
Amg::toString(bestCand->sensorDirection())<<
" found with pull "<<minPull);
276 combinatoricHits.push_back(bestCand);
279 return combinatoricHits;
282 std::unique_ptr<SegmentSeed>
283 CombinatorialNSWSeedFinderAlg::buildSegmentSeed(
const InitialSeed_t& initialSeed,
289 std::array<double, 4>
params = defineParameters(bMatrix, initialSeed);
291 const auto [segPos, direction] = seedSolution(initialSeed,
params);
295 for (std::size_t
i = 0;
i < 4; ++
i) {
296 const double halfLength = stripHalfLength(*initialSeed[
i]);
298 if (std::abs(
params[
i]) > halfLength) {
299 ATH_MSG_VERBOSE(
"Seed Rejection: Invalid seed - outside of the strip's length");
306 double interceptX = segPos.x();
307 double interceptY = segPos.y();
311 auto extendedHits = extendHits(segPos, direction, extensionLayers, usedHits);
312 HitVec hits{initialSeed.begin(),initialSeed.end()};
313 hits.insert(
hits.end(), extendedHits.begin(), extendedHits.end());
314 return std::make_unique<SegmentSeed>(tanBeta, interceptY, tanAlpha,
315 interceptX,
hits.size(),
316 std::move(
hits),
max.parentBucket());
320 std::vector<std::unique_ptr<SegmentSeed>>
325 const HitLayVec& stripHitsLayers{hitLayers.stripHits()};
326 const std::size_t layerSize = stripHitsLayers.size();
328 std::vector<std::unique_ptr<SegmentSeed>> seeds{};
337 if (m_visionTool.isEnabled()) {
340 constexpr
double legX{0.2};
352 const double pull = std::sqrt(SeedingAux::chi2Term(hitPos, hitDir, *sp));
353 const double pull2 = (mmClust->localPosition<1>().
x() - simHit->
localPosition().x()) / std::sqrt(mmClust->localCovariance<1>().x());
358 m_visionTool->visualizeBucket(Gaudi::Hive::currentContext(), *
max.parentBucket(),
359 "truth", std::move(primitives));
364 std::array<const SpacePoint*, 4> seedHits{};
368 for (std::size_t
i = 0;
i < layerSize - 3; ++
i) {
369 seedHits[0] = stripHitsLayers[
i].front();
370 for (std::size_t j =
i + 1; j < layerSize - 2; ++j) {
371 seedHits[1] = stripHitsLayers[j].front();
372 for (std::size_t
k = j + 1;
k < layerSize - 1; ++
k) {
373 seedHits[2] = stripHitsLayers[
k].front();
374 for (std::size_t
l =
k + 1;
l < layerSize; ++
l) {
375 seedHits[3] = stripHitsLayers[
l].front();
377 if (std::abs(bMatrix.determinant()) < 1.e-6) {
381 const HitLaySpan_t layers{stripHitsLayers[
i], stripHitsLayers[j], stripHitsLayers[
k], stripHitsLayers[
l]};
382 UsedHitSpan_t usedHits{allUsedHits[
i], allUsedHits[j], allUsedHits[
k], allUsedHits[
l]};
384 constructPrelimnarySeeds(globToLocal.translation(),
layers, usedHits, preLimSeeds);
389 usedExtensionHits.reserve(stripHitsLayers.size());
390 extensionLayers.reserve(stripHitsLayers.size());
391 for (std::size_t
e = 0 ;
e < stripHitsLayers.size(); ++
e) {
392 if (!(
e ==
i ||
e == j ||
e ==
k ||
e ==
l)){
393 extensionLayers.emplace_back(stripHitsLayers[
e]);
394 usedExtensionHits.emplace_back(allUsedHits[
e]);
399 for (
auto &combinatoricHits : preLimSeeds) {
400 auto seed = buildSegmentSeed(combinatoricHits, bMatrix,
max, extensionLayers, usedExtensionHits);
402 markHitsAsUsed(*seed,stripHitsLayers, allUsedHits);
403 seeds.push_back(std::move(seed));
413 void CombinatorialNSWSeedFinderAlg::markHitsAsUsed(
const SegmentSeed& seed,
417 for (
const auto* sp : seed.getHitsInMax()) {
419 for (std::size_t lIdx = 0; !
found && lIdx < allSortHits.size(); ++lIdx) {
420 const HitVec& hVec{allSortHits[lIdx]};
421 for (std::size_t hIdx = 0 ; hIdx < hVec.size(); ++hIdx) {
422 if (hVec[hIdx] == sp) {
423 usedHitMarker[lIdx][hIdx] =
true;
441 ATH_CHECK(writeMaxima.record(std::make_unique<SegmentSeedContainer>()));
446 std::vector<std::unique_ptr<SegmentSeed>> seeds = findSeedsFromMaximum(*
max, *gctx);
448 for(
const auto& hitMax :
max->getHitsInMax()){
449 ATH_MSG_VERBOSE(
"Hit "<<m_idHelperSvc->toString(hitMax->identify())<<
", "
454 for (
auto &seed : seeds) {
456 std::stringstream sstr{};
457 sstr<<
"Seed tanBeta = "<<seed->tanBeta()<<
", y0 = "<<seed->interceptY()
458 <<
", tanAlpha = "<<seed->tanAlpha()<<
", x0 = "<<seed->interceptX()<<
", hits in the seed "
459 <<seed->getHitsInMax().size()<<std::endl;
461 for(
const auto& hit : seed->getHitsInMax()){
462 sstr<<
" *** Hit "<<m_idHelperSvc->toString(hit->identify())<<
", "
467 if (m_visionTool.isEnabled()) {
468 m_visionTool->visualizeSeed(ctx, *seed,
"#phi-combinatorialSeed");
470 writeMaxima->push_back(std::move(seed));
473 return StatusCode::SUCCESS;