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