20 #include "Acts/Seeding/CombinatorialSeedSolver.hpp"
25 using namespace Acts::Experimental::CombinatorialSeedSolver;
33 const auto*
re = prd->readoutElement();
34 switch(prd->channelType()) {
36 return re->stripDesign(prd->measurementHash());
38 return re->wireDesign(prd->measurementHash());
40 return re->padDesign(prd->measurementHash());
46 const auto& design = getDesign(sp);
49 return 0.5* design.stripLength(prd->channelNumber());
57 using namespace SegmentFit;
65 ATH_CHECK(m_visionTool.retrieve(DisableTool{m_visionTool.empty()}));
67 if (!(m_idHelperSvc->hasMM() || m_idHelperSvc->hasSTGC())) {
68 ATH_MSG_ERROR(
"MM or STGC not part of initialized detector layout");
69 return StatusCode::FAILURE;
72 return StatusCode::SUCCESS;
76 CombinatorialNSWSeedFinderAlg::emptyBookKeeper(
const HitLayVec& sortedSp)
const{
78 for (std::size_t
l = 0;
l < sortedSp.size(); ++
l) {
79 emptyKeeper[
l].resize(sortedSp[
l].
size(), 0);
85 CombinatorialNSWSeedFinderAlg::classifyStrip(
const SpacePoint& sp)
const{
88 const auto& design = getDesign(sp);
89 if (!design.hasStereoAngle()) {
92 return design.stereoAngle() > 0. ? StripOrient::U : StripOrient::V;
105 CombinatorialNSWSeedFinderAlg::hitFromIPCorridor(
const SpacePoint& testHit,
110 const Amg::Vector3D estPlaneArrivalUp = SeedingAux::extrapolateToPlane(beamSpotPos, dirEstUp, testHit);
111 const Amg::Vector3D estPlaneArrivalDn = SeedingAux::extrapolateToPlane(beamSpotPos, dirEstDn, testHit);
113 bool below{
true}, above{
true};
114 switch (classifyStrip(testHit)) {
118 const double halfLength = 0.5* stripHalfLength(testHit);
124 below = estPlaneArrivalDn.y() >
std::max(leftEdge.y(), rightEdge.y());
126 above = estPlaneArrivalUp.y() <
std::min(leftEdge.y(), rightEdge.y());
132 below = estPlaneArrivalDn.y() > hY;
134 above = estPlaneArrivalUp.y() < hY;
146 << (below || above ?
" is outside the window" :
" is inside the window"));
148 return HitWindow::tooLow;
151 return HitWindow::tooHigh;
157 #define TEST_HIT_CORRIDOR(LAYER, HIT_ITER, START_LAYER) \
159 const SpacePoint* testMe = combinatoricLayers[LAYER].get()[HIT_ITER]; \
160 if (usedHits[LAYER].get()[HIT_ITER]) { \
161 ATH_MSG_VERBOSE(__func__<<":"<<__LINE__<<" - " \
162 <<m_idHelperSvc->toString(testMe->identify()) \
163 <<" already used in good seed." ); \
166 const HitWindow inWindow = hitFromIPCorridor(*testMe, beamSpot, dirEstUp, dirEstDn); \
167 if(inWindow == HitWindow::tooHigh) { \
168 ATH_MSG_VERBOSE(__func__<<":"<<__LINE__<<" - Hit " \
169 <<m_idHelperSvc->toString(testMe->identify()) \
170 <<" is beyond the corridor. Break loop"); \
172 } else if (inWindow == HitWindow::tooLow) { \
173 START_LAYER = HIT_ITER + 1; \
174 ATH_MSG_VERBOSE(__func__<<":"<<__LINE__<<" - Hit " \
175 <<m_idHelperSvc->toString(testMe->identify()) \
176 <<" is still below the corridor. Update start to " \
187 seedHitsFromLayers.clear();
188 std::size_t maxSize{1};
189 for (
const HitVec& hitVec : combinatoricLayers) {
190 maxSize = maxSize * hitVec.size();
192 seedHitsFromLayers.reserve(maxSize);
195 unsigned int iterLay0{0}, iterLay1{0}, iterLay2{0}, iterLay3{0};
196 unsigned int startLay1{0}, startLay2{0}, startLay3{0};
198 for( ; iterLay0 < combinatoricLayers[0].get().
size() ; ++iterLay0){
200 if (usedHits[0].
get()[iterLay0]) {
203 const SpacePoint* hit0 = combinatoricLayers[0].get()[iterLay0];
214 for( iterLay1 = startLay1; iterLay1 < combinatoricLayers[1].get().
size() ; ++iterLay1){
216 for( iterLay2 = startLay2; iterLay2 < combinatoricLayers[2].get().
size() ; ++iterLay2){
218 for( iterLay3 = startLay3; iterLay3 < combinatoricLayers[3].get().
size(); ++iterLay3){
220 seedHitsFromLayers.emplace_back(
std::array{hit0, combinatoricLayers[1].get()[iterLay1],
221 combinatoricLayers[2].
get()[iterLay2],
222 combinatoricLayers[3].
get()[iterLay3]});
228 #undef TEST_HIT_CORRIDOR
240 for (std::size_t
i = 0;
i < extensionLayers.size(); ++
i) {
242 const Amg::Vector3D extrapPos = SeedingAux::extrapolateToPlane(startPos, direction, *
layer.front());
244 unsigned int indexOfHit =
layer.size() + 1;
245 unsigned int triedHit{0};
249 for (
unsigned int j = 0; j <
layer.size(); ++j) {
250 if (usedHits[
i].
get().at(j)) {
253 auto hit =
layer.at(j);
254 const double pull = std::sqrt(SeedingAux::chi2Term(extrapPos, direction, *hit));
255 ATH_MSG_VERBOSE(
"Trying extension with hit " << m_idHelperSvc->toString(hit->identify()));
258 if (
pull > minPull) {
272 if (minPull < m_minPullThreshold) {
273 const auto* bestCand =
layer.at(indexOfHit);
274 ATH_MSG_VERBOSE(
"Extension successfull - hit" << m_idHelperSvc->toString(bestCand->identify())
276 <<
", dir: "<<
Amg::toString(bestCand->sensorDirection())<<
" found with pull "<<minPull);
277 combinatoricHits.push_back(bestCand);
280 return combinatoricHits;
283 std::unique_ptr<SegmentSeed>
284 CombinatorialNSWSeedFinderAlg::buildSegmentSeed(
const InitialSeed_t& initialSeed,
290 std::array<double, 4>
params = defineParameters(bMatrix, initialSeed);
292 const auto [segPos, direction] = seedSolution(initialSeed,
params);
296 for (std::size_t
i = 0;
i < 4; ++
i) {
297 const double halfLength = stripHalfLength(*initialSeed[
i]);
299 if (std::abs(
params[
i]) > halfLength) {
300 ATH_MSG_VERBOSE(
"Seed Rejection: Invalid seed - outside of the strip's length");
307 double interceptX = segPos.x();
308 double interceptY = segPos.y();
312 auto extendedHits = extendHits(segPos, direction, extensionLayers, usedHits);
313 HitVec hits{initialSeed.begin(),initialSeed.end()};
314 hits.insert(
hits.end(), extendedHits.begin(), extendedHits.end());
315 return std::make_unique<SegmentSeed>(tanTheta, interceptY, tanPhi,
316 interceptX,
hits.size(),
317 std::move(
hits),
max.parentBucket());
321 std::vector<std::unique_ptr<SegmentSeed>>
326 const HitLayVec& stripHitsLayers{hitLayers.stripHits()};
327 const std::size_t layerSize = stripHitsLayers.size();
329 std::vector<std::unique_ptr<SegmentSeed>> seeds{};
338 if (m_visionTool.isEnabled()) {
341 constexpr
double legX{0.2};
353 const double pull = std::sqrt(SeedingAux::chi2Term(hitPos, hitDir, *sp));
354 const double pull2 = (mmClust->localPosition<1>().
x() - simHit->
localPosition().x()) / std::sqrt(mmClust->localCovariance<1>().x());
359 m_visionTool->visualizeBucket(Gaudi::Hive::currentContext(), *
max.parentBucket(),
360 "truth", std::move(primitives));
365 std::array<const SpacePoint*, 4> seedHits{};
369 for (std::size_t
i = 0;
i < layerSize - 3; ++
i) {
370 seedHits[0] = stripHitsLayers[
i].front();
371 for (std::size_t j =
i + 1; j < layerSize - 2; ++j) {
372 seedHits[1] = stripHitsLayers[j].front();
373 for (std::size_t
k = j + 1;
k < layerSize - 1; ++
k) {
374 seedHits[2] = stripHitsLayers[
k].front();
375 for (std::size_t
l =
k + 1;
l < layerSize; ++
l) {
376 seedHits[3] = stripHitsLayers[
l].front();
378 if (std::abs(bMatrix.determinant()) < 1.e-6) {
382 const HitLaySpan_t layers{stripHitsLayers[
i], stripHitsLayers[j], stripHitsLayers[
k], stripHitsLayers[
l]};
383 UsedHitSpan_t usedHits{allUsedHits[
i], allUsedHits[j], allUsedHits[
k], allUsedHits[
l]};
385 constructPrelimnarySeeds(globToLocal.translation(),
layers, usedHits, preLimSeeds);
390 usedExtensionHits.reserve(stripHitsLayers.size());
391 extensionLayers.reserve(stripHitsLayers.size());
392 for (std::size_t
e = 0 ;
e < stripHitsLayers.size(); ++
e) {
393 if (!(
e ==
i ||
e == j ||
e ==
k ||
e ==
l)){
394 extensionLayers.emplace_back(stripHitsLayers[
e]);
395 usedExtensionHits.emplace_back(allUsedHits[
e]);
400 for (
auto &combinatoricHits : preLimSeeds) {
401 auto seed = buildSegmentSeed(combinatoricHits, bMatrix,
max, extensionLayers, usedExtensionHits);
403 markHitsAsUsed(*seed,stripHitsLayers, allUsedHits);
404 seeds.push_back(std::move(seed));
414 void CombinatorialNSWSeedFinderAlg::markHitsAsUsed(
const SegmentSeed& seed,
418 for (
const auto* sp : seed.getHitsInMax()) {
420 for (std::size_t lIdx = 0; !
found && lIdx < allSortHits.size(); ++lIdx) {
421 const HitVec& hVec{allSortHits[lIdx]};
422 for (std::size_t hIdx = 0 ; hIdx < hVec.size(); ++hIdx) {
423 if (hVec[hIdx] == sp) {
424 usedHitMarker[lIdx][hIdx] =
true;
442 ATH_CHECK(writeMaxima.record(std::make_unique<SegmentSeedContainer>()));
447 std::vector<std::unique_ptr<SegmentSeed>> seeds = findSeedsFromMaximum(*
max, *gctx);
449 for(
const auto& hitMax :
max->getHitsInMax()){
450 ATH_MSG_VERBOSE(
"Hit "<<m_idHelperSvc->toString(hitMax->identify())<<
", "
455 for (
auto &seed : seeds) {
457 std::stringstream sstr{};
458 sstr<<
"Seed tanTheta = "<<seed->tanTheta()<<
", y0 = "<<seed->interceptY()
459 <<
", tanPhi = "<<seed->tanPhi()<<
", x0 = "<<seed->interceptX()<<
", hits in the seed "
460 <<seed->getHitsInMax().size()<<std::endl;
462 for(
const auto& hit : seed->getHitsInMax()){
463 sstr<<
" *** Hit "<<m_idHelperSvc->toString(hit->identify())<<
", "
468 if (m_visionTool.isEnabled()) {
469 m_visionTool->visualizeSeed(ctx, *seed,
"#phi-combinatorialSeed");
471 writeMaxima->push_back(std::move(seed));
474 return StatusCode::SUCCESS;