445 {
446
447
448 const float MAX_R = 1030.0;
449 const float MAX_Z = 3000.0;
450 const float maxCornerDist = 15.0;
451
452
453
454 MagField::AtlasFieldCache fieldCache;
455
459 return -1;
460 }
461
462 fieldCondObj->getInitializedCache (fieldCache);
463
464 road.clear();
465
466 unsigned int nSP =
seed.size();
467
468 if(nSP < 3) return -2;
469
470 std::vector<unsigned int> seedHashes;
471
472 for(unsigned int spIdx=0;spIdx<nSP;spIdx++) {
473 const auto&
sp =
seed.at(spIdx);
474 const Trk::PrepRawData* prd =
sp->clusterList().first;
475 const InDet::PixelCluster* pPixelHit = static_cast<const InDet::PixelCluster*>(prd);
477 seedHashes.push_back(hash);
478 }
479
480 std::vector<std::array<float,2> > zr;
481 zr.resize(nSP+2);
482
483 for(unsigned int spIdx=0;spIdx<nSP;spIdx++) {
484 const auto&
sp =
seed.at(spIdx);
485 zr[spIdx+1][0] =
sp->globalPosition().z();
486 zr[spIdx+1][1] =
sp->globalPosition().perp();
487 }
488
489
490
491 zr[0][0] = zr[1][0] - zr[1][1]*(zr[2][0]-zr[1][0])/(zr[2][1]-zr[1][1]);
492 zr[0][1] = 0;
493
494
495 float zlast = zr[nSP-1][0] + (MAX_R - zr[nSP-1][1])*(zr[nSP][0]-zr[nSP-1][0])/(zr[nSP][1]-zr[nSP-1][1]);
496 float rlast = MAX_R;
497
498 if (std::fabs(zlast) > MAX_Z) {
499 if(zlast > 0) zlast = MAX_Z;
500 else zlast = -MAX_Z;
501 rlast = zr[nSP-1][1] + (zlast - zr[nSP-1][0])*(zr[nSP][1]-zr[nSP-1][1])/(zr[nSP][0]-zr[nSP-1][0]);
502 }
503 zr[nSP+1][0] = zlast;
504 zr[nSP+1][1] = rlast;
505
506 std::map<unsigned int, SearchInterval> rzIntervals;
507
508 for(unsigned int k1 = 0;k1<zr.size()-1;k1++) {
509
510 unsigned int k2 = k1+1;
511
512 float z1 = zr[k1][0];
513 float r1 = zr[k1][1];
514 float z2 = zr[k2][0];
515 float r2 = zr[k2][1];
516 float dz21 = z2-z1;
518 float L = std::sqrt(dz21*dz21 + dr21*dr21);
520 float sinF = dr21*invL;
521 float cosF = dz21*invL;
522
524
525 if (z1 < vb.m_zr[0] && z2 < vb.m_zr[0]) continue;
526 if (z1 > vb.m_zr[1] && z2 > vb.m_zr[1]) continue;
527 if (r1 < vb.m_zr[2] && r2 < vb.m_zr[2]) continue;
528 if (r1 > vb.m_zr[3] && r2 > vb.m_zr[3]) continue;
529
530
531
532 float zc[4] = {vb.m_zr[0], vb.m_zr[1], vb.m_zr[1], vb.m_zr[0]};
533 float rc[4] = {vb.m_zr[2], vb.m_zr[2], vb.m_zr[3], vb.m_zr[3]};
534
535 int nUp(0), nDn(0);
536
537 float minDistances[4];
538
539 for (
int ic=0;
ic<4;
ic++) {
540 minDistances[
ic] = (
rc[
ic] -
r1)*cosF - (zc[ic] - z1)*sinF;
541 }
542
543 float minH = std::abs(minDistances[0]);
544
545 for (
int ic=0;
ic<4;
ic++) {
546 float h = minDistances[
ic];
547 if (h <=0) nDn += 1;
548 else nUp += 1;
549 if(std::abs(h) < minH) minH = std::abs(h);
550 }
551
552 if (nUp == 4 || nDn == 4) {
553 if(minH > maxCornerDist) {
554 continue;
555 }
556 }
557
558
559
560 for(auto lIdx : vb.m_layers) {
562 for(
int s1=0;
s1<
lb.m_nVertices-1;
s1++) {
564
565 float h1 = (
lb.m_r[
s1] -
r1)*cosF - (
lb.m_z[s1] - z1)*sinF;
566 float h2 = (
lb.m_r[
s2] -
r1)*cosF - (
lb.m_z[s2] - z1)*sinF;
567 if (h1*h2 > 0) continue;
568 float l1 = (
lb.m_z[
s1] - z1)*cosF + (
lb.m_r[s1] - r1)*sinF;
569 float l2 = (
lb.m_z[
s2] - z1)*cosF + (
lb.m_r[s2] - r1)*sinF;
570 if (l1 < 0 && l2 < 0) continue;
571 if (l1 > L && l2 > L) continue;
573
574 if (lx < 0 || lx > 1) continue;
575 float zx = z2*lx + (1-lx)*z1;
576 float rx =
r2*lx + (1-lx)*r1;
577
578 auto radItr = rzIntervals.find(
lb.m_lay_id);
579
580 if(radItr == rzIntervals.end()) {
582 }
583 else {
584 (*radItr).second.addPoint(zx, rx);
585 }
586 }
587 }
588 }
589 }
590
591
592
593 unsigned int sp1Idx = 0;
594 unsigned int sp3Idx = nSP-1;
595 unsigned int sp2Idx = (sp3Idx+sp1Idx)/2;
596
597
598
599 float uv_coords[2][2];
600
601 float dx =
seed.at(sp3Idx)->globalPosition().x() -
seed.at(sp2Idx)->globalPosition().x();
602 float dy =
seed.at(sp3Idx)->globalPosition().y() -
seed.at(sp2Idx)->globalPosition().y();
603
604 uv_coords[1][0] = -std::sqrt(dx*dx + dy*dy);
605 uv_coords[1][1] = 0.0;
606
607 float cos_theta =
dx/(-uv_coords[1][0]);
608 float sin_theta =
dy/(-uv_coords[1][0]);
609
610 float rot_matrix[2][2];
611 float inv_rot_matrix[2][2];
612
613 rot_matrix[0][0] = cos_theta;
614 rot_matrix[0][1] = sin_theta;
615 rot_matrix[1][0] = -sin_theta;
616 rot_matrix[1][1] = cos_theta;
617
618 inv_rot_matrix[0][0] = cos_theta;
619 inv_rot_matrix[0][1] = -sin_theta;
620 inv_rot_matrix[1][0] = sin_theta;
621 inv_rot_matrix[1][1] = cos_theta;
622
623 float sp3_coords[3];
624 sp3_coords[0] =
seed.at(sp3Idx)->globalPosition().x();
625 sp3_coords[1] =
seed.at(sp3Idx)->globalPosition().y();
626 sp3_coords[2] =
seed.at(sp3Idx)->globalPosition().z();
627
628
629
630 float u_c = rot_matrix[0][0]*(0-sp3_coords[0]) + rot_matrix[0][1]*(0-sp3_coords[1]);
631 float v_c = rot_matrix[1][0]*(0-sp3_coords[0]) + rot_matrix[1][1]*(0-sp3_coords[1]);
632
633 int sign_up = u_c < 0 ? 1 : -1;
634
635
636
637 float dR1[2];
638
639 dR1[0] =
seed.at(sp1Idx)->globalPosition().x() - sp3_coords[0];
640 dR1[1] =
seed.at(sp1Idx)->globalPosition().y() - sp3_coords[1];
641
642 uv_coords[0][0] = rot_matrix[0][0]*dR1[0] + rot_matrix[0][1]*dR1[1];
643 uv_coords[0][1] = rot_matrix[1][0]*dR1[0] + rot_matrix[1][1]*dR1[1];
644
645
646
647 float a = uv_coords[0][1]/(uv_coords[0][0]*(uv_coords[0][0]-uv_coords[1][0]));
648 float b = -
a*uv_coords[1][0];
649
650
651
652 std::vector<unsigned int> pixelHashIds;
653 std::vector<unsigned int> stripHashIds;
654
655 for(auto& ip : rzIntervals) {
656
657 float R =
ip.second.getAverageRadius();
658
660
661 float dRv = R2-v_c*v_c;
662
663 if(dRv < 0) continue;
664
665 float u_p = u_c + sign_up*std::sqrt(dRv);
666
667 float v_p =
b*u_p +
a*u_p*u_p;
668 float dv2 = (v_p-v_c)*(v_p-v_c);
669
670 if (R2-dv2<0) continue;
671
672 float u_star = u_c + sign_up*std::sqrt(R2-dv2);
673 float v_star =
b*u_star +
a*u_star*u_star;
674
675 float x_star = inv_rot_matrix[0][0]*u_star + inv_rot_matrix[0][1]*v_star + sp3_coords[0];
676 float y_star = inv_rot_matrix[1][0]*u_star + inv_rot_matrix[1][1]*v_star + sp3_coords[1];
677
678 ip.second.m_r = std::sqrt(x_star*x_star + y_star*y_star);
679 ip.second.m_phi = std::atan2(y_star, x_star);
680
681 if(
ip.first > 15000) {
682
683 bool hasHit = false;
684
685 for(
unsigned int i=1;
i<nSP-1;
i++) {
686
687 float dpr =
ip.second.m_r - zr[
i][1];
688 float dpz =
ip.second.m_z - zr[
i][0];
689 float dist = std::sqrt(dpr*dpr + dpz*dpz);
690 if(dist < maxCornerDist) {
691 hasHit = true;
692 break;
693 }
694 }
696 }
697 else {
699 }
700 }
701
702 std::set<unsigned int> pixelHashSet(pixelHashIds.begin(), pixelHashIds.end());
703
704 for(auto id : seedHashes) {
705 if(pixelHashSet.find(id) == pixelHashSet.end()) pixelHashIds.push_back(id);
706 }
707
708 std::vector<std::pair<float, const InDetDD::SiDetectorElement*> > theRoad;
709
710 for(auto hash_id : pixelHashIds) {
711 const InDetDD::SiDetectorElement *
p =
m_pixelManager->getDetectorElement(hash_id);
712 if(p == nullptr) continue;
714 float dist = std::sqrt(
C(0)*
C(0) +
C(1)*
C(1) +
C(2)*
C(2));
715 theRoad.push_back(std::make_pair(dist,p));
716 }
717
718 for(auto hash_id : stripHashIds) {
719 const InDetDD::SiDetectorElement *
p =
m_stripManager->getDetectorElement(hash_id);
720 if(p == nullptr) continue;
722 float dist = std::sqrt(
C(0)*
C(0) +
C(1)*
C(1) +
C(2)*
C(2));
723 theRoad.push_back(std::make_pair(dist,p));
724 }
725
726 std::sort(theRoad.begin(), theRoad.end());
727
728 for(const auto & dp : theRoad) {
729 road.push_back(
dp.second);
730 }
731
732 return (int)theRoad.size();
733}
#define ATH_MSG_ERROR(x,...)
virtual IdentifierHash identifyHash() const override final
identifier hash (inline)
virtual const InDetDD::SiDetectorElement * detectorElement() const override final
return the detector element corresponding to this PRD The pointer will be zero if the det el is not d...
double R(const INavigable4Momentum *p1, const double v_eta, const double v_phi)
void sort(typename DataModel_detail::iterator< DVL > beg, typename DataModel_detail::iterator< DVL > end)
Specialization of sort for DataVector/List.
1.17.0