Estimation of the geometrical volume span.
428 {
429
430 const CuboidVolumeBounds* box =
431 dynamic_cast<const CuboidVolumeBounds*>(&volBounds);
432 const TrapezoidVolumeBounds* trd =
433 dynamic_cast<const TrapezoidVolumeBounds*>(&volBounds);
434 const DoubleTrapezoidVolumeBounds* dtrd =
435 dynamic_cast<const DoubleTrapezoidVolumeBounds*>(&volBounds);
436 const BevelledCylinderVolumeBounds* bcyl =
437 dynamic_cast<const BevelledCylinderVolumeBounds*>(&volBounds);
438 const CylinderVolumeBounds* cyl =
439 dynamic_cast<const CylinderVolumeBounds*>(&volBounds);
440 const SubtractedVolumeBounds* sub =
441 dynamic_cast<const SubtractedVolumeBounds*>(&volBounds);
442 const CombinedVolumeBounds* comb =
443 dynamic_cast<const CombinedVolumeBounds*>(&volBounds);
444 const SimplePolygonBrepVolumeBounds* spb =
445 dynamic_cast<const SimplePolygonBrepVolumeBounds*>(&volBounds);
446 const PrismVolumeBounds* prism =
447 dynamic_cast<const PrismVolumeBounds*>(&volBounds);
448
450
451 if (sub) {
453 transform * sub->outer()->transform(), zTol,
454 phiTol);
455 }
456
457 if (comb) {
459 comb->first()->volumeBounds(),
460 transform * comb->first()->transform(), zTol, phiTol);
462 comb->second()->volumeBounds(),
463 transform * comb->second()->transform(), zTol, phiTol);
464
465 VolumeSpan scomb;
466 scomb.rMin = std::min((*s1).rMin, (*s2).rMin);
467 scomb.rMax = std::max((*s1).rMax, (*s2).rMax);
468 scomb.xMin = std::min((*s1).xMin, (*s2).xMin);
469 scomb.xMax = std::max((*s1).xMax, (*s2).xMax);
470 scomb.yMin = std::min((*s1).yMin, (*s2).yMin);
471 scomb.yMax = std::max((*s1).yMax, (*s2).yMax);
472 scomb.zMin = std::min((*s1).zMin, (*s2).zMin);
473 scomb.zMax = std::max((*s1).zMax, (*s2).zMax);
474 if ((*s1).phiMin < (*s1).phiMax && (*s2).phiMin < (*s2).phiMax) {
475 scomb.phiMin = std::min((*s1).phiMin, (*s2).phiMin);
476 scomb.phiMax = std::max((*s1).phiMax, (*s2).phiMax);
477 } else if ((*s1).phiMin < (*s1).phiMax && (*s2).phiMin > (*s2).phiMax) {
478 if ((*s1).phiMin > (*s2).phiMax) {
479 scomb.phiMin = std::min((*s1).phiMin, (*s2).phiMin);
480 scomb.phiMax = (*s2).phiMax;
481 } else if ((*s1).phiMax < (*s2).phiMin) {
482 scomb.phiMin = (*s2).phiMin;
483 scomb.phiMax = std::max((*s1).phiMax, (*s2).phiMax);
484 } else {
485 scomb.phiMin = 0.;
486 scomb.phiMax = 2 *
M_PI;
487 }
488 } else if ((*s1).phiMin > (*s1).phiMax && (*s2).phiMin < (*s2).phiMax) {
489 if ((*s2).phiMin > (*s1).phiMax) {
490 scomb.phiMin = std::min((*s1).phiMin, (*s2).phiMin);
491 scomb.phiMax = (*s1).phiMax;
492 } else if ((*s2).phiMax < (*s1).phiMin) {
493 scomb.phiMin = (*s1).phiMin;
494 scomb.phiMax = std::max((*s1).phiMax, (*s2).phiMax);
495 } else {
496 scomb.phiMin = 0.;
497 scomb.phiMax = 2 *
M_PI;
498 }
499 } else {
500 scomb.phiMin = std::min((*s1).phiMin, (*s2).phiMin);
501 scomb.phiMax = std::max((*s1).phiMax, (*s2).phiMax);
502 }
503 return std::make_unique<VolumeSpan>(scomb);
504 }
505
506
507 double minZ{1.e6};
508 double maxZ{-1.e6};
509 double minPhi{2 *
M_PI};
510 double maxPhi{0.};
511 double minR{1.e6};
512 double maxR{0.};
513 double minX{1.e6};
514 double maxX{-1.e6};
515 double minY{1.e6};
516 double maxY{-1.e6};
517
518
519 std::vector<Amg::Vector3D> vtx;
520 std::vector<std::pair<int, int>> edges;
522
523 if (box) {
524 vtx.emplace_back(box->halflengthX(), box->halflengthY(),
525 box->halflengthZ());
526 vtx.emplace_back(-box->halflengthX(), box->halflengthY(),
527 box->halflengthZ());
528 vtx.emplace_back(box->halflengthX(), -box->halflengthY(),
529 box->halflengthZ());
530 vtx.emplace_back(-box->halflengthX(), -box->halflengthY(),
531 box->halflengthZ());
532 vtx.emplace_back(box->halflengthX(), box->halflengthY(),
533 -box->halflengthZ());
534 vtx.emplace_back(-box->halflengthX(), box->halflengthY(),
535 -box->halflengthZ());
536 vtx.emplace_back(box->halflengthX(), -box->halflengthY(),
537 -box->halflengthZ());
538 vtx.emplace_back(-box->halflengthX(), -box->halflengthY(),
539 -box->halflengthZ());
540 edges.emplace_back(0, 1);
541 edges.emplace_back(0, 2);
542 edges.emplace_back(1, 3);
543 edges.emplace_back(2, 3);
544 edges.emplace_back(4, 5);
545 edges.emplace_back(4, 6);
546 edges.emplace_back(5, 7);
547 edges.emplace_back(6, 7);
548 edges.emplace_back(0, 4);
549 edges.emplace_back(1, 5);
550 edges.emplace_back(2, 6);
551 edges.emplace_back(3, 7);
552 }
553 if (trd) {
554 vtx.emplace_back(trd->maxHalflengthX(), trd->halflengthY(),
555 trd->halflengthZ());
556 vtx.emplace_back(-trd->maxHalflengthX(), trd->halflengthY(),
557 trd->halflengthZ());
558 vtx.emplace_back(trd->minHalflengthX(), -trd->halflengthY(),
559 trd->halflengthZ());
560 vtx.emplace_back(-trd->minHalflengthX(), -trd->halflengthY(),
561 trd->halflengthZ());
562 vtx.emplace_back(trd->maxHalflengthX(), trd->halflengthY(),
563 -trd->halflengthZ());
564 vtx.emplace_back(-trd->maxHalflengthX(), trd->halflengthY(),
565 -trd->halflengthZ());
566 vtx.emplace_back(trd->minHalflengthX(), -trd->halflengthY(),
567 -trd->halflengthZ());
568 vtx.emplace_back(-trd->minHalflengthX(), -trd->halflengthY(),
569 -trd->halflengthZ());
570 edges.emplace_back(0, 1);
571 edges.emplace_back(0, 2);
572 edges.emplace_back(1, 3);
573 edges.emplace_back(2, 3);
574 edges.emplace_back(4, 5);
575 edges.emplace_back(4, 6);
576 edges.emplace_back(5, 7);
577 edges.emplace_back(6, 7);
578 edges.emplace_back(0, 4);
579 edges.emplace_back(1, 5);
580 edges.emplace_back(2, 6);
581 edges.emplace_back(3, 7);
582 }
583 if (dtrd) {
584 vtx.emplace_back(dtrd->maxHalflengthX(), 2 * dtrd->halflengthY2(),
585 dtrd->halflengthZ());
586 vtx.emplace_back(-dtrd->maxHalflengthX(), 2 * dtrd->halflengthY2(),
587 dtrd->halflengthZ());
588 vtx.emplace_back(dtrd->medHalflengthX(), 0., dtrd->halflengthZ());
589 vtx.emplace_back(-dtrd->medHalflengthX(), 0., dtrd->halflengthZ());
590 vtx.emplace_back(dtrd->minHalflengthX(), -2 * dtrd->halflengthY1(),
591 dtrd->halflengthZ());
592 vtx.emplace_back(-dtrd->minHalflengthX(), -2 * dtrd->halflengthY1(),
593 dtrd->halflengthZ());
594 vtx.emplace_back(dtrd->maxHalflengthX(), 2 * dtrd->halflengthY2(),
595 -dtrd->halflengthZ());
596 vtx.emplace_back(-dtrd->maxHalflengthX(), 2 * dtrd->halflengthY2(),
597 -dtrd->halflengthZ());
598 vtx.emplace_back(dtrd->medHalflengthX(), 0., -dtrd->halflengthZ());
599 vtx.emplace_back(-dtrd->medHalflengthX(), 0., -dtrd->halflengthZ());
600 vtx.emplace_back(dtrd->minHalflengthX(), -2 * dtrd->halflengthY1(),
601 -dtrd->halflengthZ());
602 vtx.emplace_back(-dtrd->minHalflengthX(), -2 * dtrd->halflengthY1(),
603 -dtrd->halflengthZ());
604 edges.emplace_back(0, 1);
605 edges.emplace_back(0, 2);
606 edges.emplace_back(1, 3);
607 edges.emplace_back(2, 4);
608 edges.emplace_back(3, 5);
609 edges.emplace_back(4, 5);
610 edges.emplace_back(6, 7);
611 edges.emplace_back(6, 8);
612 edges.emplace_back(7, 9);
613 edges.emplace_back(8, 10);
614 edges.emplace_back(9, 11);
615 edges.emplace_back(10, 11);
616 edges.emplace_back(0, 6);
617 edges.emplace_back(1, 7);
618 edges.emplace_back(2, 8);
619 edges.emplace_back(3, 9);
620 edges.emplace_back(4, 10);
621 edges.emplace_back(5, 11);
622 }
623 if (bcyl) {
624 dPhi = bcyl->halfPhiSector();
625 vtx.emplace_back(0., 0., bcyl->halflengthZ());
626 vtx.emplace_back(0., 0., -bcyl->halflengthZ());
627 edges.emplace_back(0, 1);
629 const double cosDphi = std::cos(dPhi);
630 const double sinDphi = std::sin(dPhi);
631 vtx.emplace_back(bcyl->outerRadius() * cosDphi,
632 bcyl->outerRadius() * sinDphi,
633 bcyl->halflengthZ());
634 vtx.emplace_back(bcyl->innerRadius() * cosDphi,
635 bcyl->innerRadius() * sinDphi,
636 bcyl->halflengthZ());
637 vtx.emplace_back(bcyl->outerRadius() * cosDphi,
638 -bcyl->outerRadius() * sinDphi,
639 bcyl->halflengthZ());
640 vtx.emplace_back(bcyl->innerRadius() * cosDphi,
641 -bcyl->innerRadius() * sinDphi,
642 bcyl->halflengthZ());
643 vtx.emplace_back(bcyl->outerRadius() * cosDphi,
644 bcyl->outerRadius() * sinDphi,
645 -bcyl->halflengthZ());
646 vtx.emplace_back(bcyl->innerRadius() * cosDphi,
647 bcyl->innerRadius() * sinDphi,
648 -bcyl->halflengthZ());
649 vtx.emplace_back(bcyl->outerRadius() * cosDphi,
650 -bcyl->outerRadius() * sinDphi,
651 -bcyl->halflengthZ());
652 vtx.emplace_back(bcyl->innerRadius() * cosDphi,
653 -bcyl->innerRadius() * sinDphi,
654 -bcyl->halflengthZ());
655 vtx.emplace_back(bcyl->outerRadius(), 0.,
656 0.);
657
658 edges.emplace_back(2, 3);
659 edges.emplace_back(4, 5);
660 edges.emplace_back(6, 7);
661 edges.emplace_back(8, 9);
662 if (bcyl->type() == 1 || bcyl->type() == 3) {
663 edges.emplace_back(3, 5);
664 edges.emplace_back(7, 9);
665 }
666 if (bcyl->type() == 2 || bcyl->type() == 3) {
667 edges.emplace_back(2, 4);
668 edges.emplace_back(6, 8);
669 }
670 }
671 }
672 if (cyl) {
673 dPhi = cyl->halfPhiSector();
674 vtx.emplace_back(0., 0., cyl->halflengthZ());
675 vtx.emplace_back(0., 0., -cyl->halflengthZ());
676 edges.emplace_back(0, 1);
678 const double cosDphi = std::cos(dPhi);
679 const double sinDphi = std::sin(dPhi);
680 vtx.emplace_back(cyl->outerRadius() * cosDphi,
681 cyl->outerRadius() * sinDphi, cyl->halflengthZ());
682 vtx.emplace_back(cyl->innerRadius() * cosDphi,
683 cyl->innerRadius() * sinDphi, cyl->halflengthZ());
684 vtx.emplace_back(cyl->outerRadius() * cosDphi,
685 -cyl->outerRadius() * sinDphi, cyl->halflengthZ());
686 vtx.emplace_back(cyl->outerRadius() * cosDphi,
687 -cyl->outerRadius() * sinDphi, cyl->halflengthZ());
688 vtx.emplace_back(cyl->outerRadius() * cosDphi,
689 cyl->outerRadius() * sinDphi, -cyl->halflengthZ());
690 vtx.emplace_back(cyl->innerRadius() * cosDphi,
691 cyl->innerRadius() * sinDphi, -cyl->halflengthZ());
692 vtx.emplace_back(cyl->outerRadius() * cosDphi,
693 -cyl->outerRadius() * sinDphi,
694 -cyl->halflengthZ());
695 vtx.emplace_back(cyl->outerRadius() * cosDphi,
696 -cyl->outerRadius() * sinDphi,
697 -cyl->halflengthZ());
698 vtx.emplace_back(cyl->outerRadius(), 0.,
699 0.);
700
701 edges.emplace_back(2, 3);
702 edges.emplace_back(4, 5);
703 edges.emplace_back(6, 7);
704 edges.emplace_back(8, 9);
705 }
706 }
707
708 if (spb) {
709 const std::vector<std::pair<double, double>> vtcs = spb->xyVertices();
710 for (const auto& vtc : vtcs) {
711 vtx.emplace_back(vtc.first, vtc.second, spb->halflengthZ());
712 vtx.emplace_back(vtc.first, vtc.second, -spb->halflengthZ());
713 edges.emplace_back(vtx.size() - 2, vtx.size() - 1);
714 if (vtx.size() > 2) {
715 edges.emplace_back(
716 vtx.size() - 4, vtx.size() - 2);
717 edges.emplace_back(
718 vtx.size() - 3, vtx.size() - 1);
719 }
720 if (vtx.size() > 4) {
721 edges.emplace_back(vtx.size() - 2, 1);
722 edges.emplace_back(vtx.size() - 1, 0);
723 }
724 }
725 edges.emplace_back(0, vtx.size() - 2);
726 edges.emplace_back(1, vtx.size() - 1);
727 }
728
729 if (prism) {
730 const std::vector<std::pair<double, double>> vtcs = prism->xyVertices();
731 for (const auto& vtc : vtcs) {
732 vtx.emplace_back(vtc.first, vtc.second, prism->halflengthZ());
733 vtx.emplace_back(vtc.first, vtc.second, -prism->halflengthZ());
734 edges.emplace_back(vtx.size() - 2, vtx.size() - 1);
735 if (vtx.size() > 2) {
736 edges.emplace_back(
737 vtx.size() - 4, vtx.size() - 2);
738 edges.emplace_back(
739 vtx.size() - 3, vtx.size() - 1);
740 }
741 }
742 edges.emplace_back(0, vtx.size() - 2);
743 edges.emplace_back(1, vtx.size() - 1);
744 }
745
746 std::vector<Amg::Vector3D> vtxt;
747
748 for (
unsigned int ie = 0;
ie < vtx.size();
ie++) {
750 vtxt.push_back(gp);
751
753 double rad = gp.perp();
754
755
756 minX = std::min(minX, gp[0]);
757 maxX = std::max(maxX, gp[0]);
758 minY = std::min(minY, gp[1]);
759 maxY = std::max(maxY, gp[1]);
760 minZ = std::min(minZ, gp[2]);
761 maxZ = std::max(maxZ, gp[2]);
762 minR = std::min(minR, rad);
763 maxR = std::max(maxR, rad);
764 maxPhi = std::max(maxPhi,
phi);
765 minPhi = std::min(minPhi,
phi);
766 }
767
768 if (cyl || bcyl) {
769
770 double ro = cyl ? cyl->outerRadius() : bcyl->outerRadius();
771 double ri = cyl ? cyl->innerRadius() : bcyl->innerRadius();
772
774 (vtxt[edges[0].first] - vtxt[edges[0].second]).
unit();
775 maxZ += ro *
sin(
dir.theta());
776 minZ += -ro *
sin(
dir.theta());
777
778 if (ro < minR) {
779
780 PerigeeSurface peri;
781 Intersection closest = peri.straightLineIntersection(vtxt[1], dir);
782 double le = (vtxt[0] - vtxt[1]).
norm();
783 if ((closest.position - vtxt[0]).norm() < le &&
784 (closest.position - vtxt[1]).norm() < le) {
785 if (minR > closest.position.perp() - ro)
786 minR = std::max(0., closest.position.perp() - ro);
787
788 double phiClosest = closest.position.phi() +
M_PI;
789 if (phiClosest < minPhi || phiClosest > maxPhi) {
790 double phiTmp = minPhi;
791 minPhi = maxPhi;
792 maxPhi = phiTmp;
793 }
794 } else
795 minR = std::max(0., minR - ro * std::abs(
dir.z()));
796
797 const double aTan = std::atan2(ro, minR);
798 minPhi += -aTan;
799 maxPhi += aTan;
800 if (minPhi < 0)
802 if (maxPhi > 2 *
M_PI)
804
805 maxR += ro * std::abs(
cos(
dir.theta()));
806 } else {
807
808 double rAx = std::max(vtxt[0].
perp(), vtxt[1].
perp());
809 if (rAx < ri)
810 minR = ri - rAx;
811 else
812 minR = std::max(0., minR - ro * std::abs(
cos(
dir.theta())));
813
814
815 PerigeeSurface peri;
816 for (
unsigned int ie = 0;
ie < edges.size();
ie++) {
818 (vtxt[edges[
ie].first] - vtxt[edges[
ie].second]).
unit();
819 Intersection closest =
820 peri.straightLineIntersection(vtxt[edges[ie].second], dir);
821 double le =
822 (vtxt[edges[
ie].first] - vtxt[edges[
ie].second]).
norm();
823 if ((closest.position - vtxt[edges[ie].first]).norm() < le &&
824 (closest.position - vtxt[edges[ie].second]).norm() < le)
825 if (minR > closest.position.perp())
826 minR = closest.position.perp();
827 }
828
829 if (vtxt.size() > 10) {
830
831
832 double phiSecLmin = std::min(
835 double phiSecLmax = std::max(
838
839
840 double phiSecUmin = std::min(
843 double phiSecUmax = std::max(
846 minPhi = std::min(std::min(phiSecLmin, phiSecLmax),
847 std::min(phiSecUmin, phiSecUmax));
848 maxPhi = std::max(std::max(phiSecLmin, phiSecLmax),
849 std::max(phiSecUmin, phiSecUmax));
850 if (vtxt[10].
phi() +
M_PI < minPhi ||
851 vtxt[10].
phi() +
M_PI > maxPhi) {
853 maxPhi = 0.;
854 double phiTmp;
855 for (unsigned int iv = 2; iv < vtxt.size(); iv++) {
856 phiTmp = vtxt[iv].phi() +
M_PI;
859 minPhi = phiTmp < minPhi ? phiTmp : minPhi;
860 maxPhi = phiTmp > maxPhi ? phiTmp : maxPhi;
861 }
862 if (minPhi > 2 *
M_PI)
864 if (maxPhi > 2 *
M_PI)
866 }
867 } else {
868 minPhi = 0.;
870 maxR += ro * std::abs(std::cos(
dir.theta()));
871 }
872 if (minPhi >= maxPhi && (minPhi - maxPhi) <
M_PI) {
873 minPhi = 0.;
875 }
876 }
877 }
878
879 if (!cyl && !bcyl) {
880
881 PerigeeSurface peri;
882 for (
unsigned int ie = 0;
ie < edges.size();
ie++) {
884 (vtxt[edges[
ie].first] - vtxt[edges[
ie].second]).
unit();
885 Intersection closest =
886 peri.straightLineIntersection(vtxt[edges[ie].second], dir);
887 double le = (vtxt[edges[
ie].first] - vtxt[edges[
ie].second]).
norm();
888 if ((closest.position - vtxt[edges[ie].first]).norm() < le &&
889 (closest.position - vtxt[edges[ie].second]).norm() < le)
890 if (minR > closest.position.perp())
891 minR = closest.position.perp();
892 }
893
894 if (std::abs(maxPhi - minPhi) >
M_PI) {
895 double phiTmp = minPhi;
897 maxPhi = 0.;
898 for (unsigned int iv = 0; iv < vtxt.size(); iv++) {
899 phiTmp = vtxt[iv].phi() +
M_PI;
902 minPhi = phiTmp < minPhi ? phiTmp : minPhi;
903 maxPhi = phiTmp > maxPhi ? phiTmp : maxPhi;
904 }
905 if (minPhi > 2 *
M_PI)
907 if (maxPhi > 2 *
M_PI)
909 if (minPhi >= maxPhi && (minPhi - maxPhi) <
M_PI) {
910 minPhi = 0.;
912 }
913 }
914 }
915
916 if (cyl || bcyl || box || trd || dtrd || spb || prism) {
917 span.zMin = minZ - zTol;
918 span.zMax = maxZ - +zTol;
919 minPhi = (minPhi - phiTol) < 0 ? minPhi - phiTol + 2 *
M_PI
920 : minPhi - phiTol;
921 span.phiMin = minPhi;
922 maxPhi = (maxPhi + phiTol) > 2 *
M_PI ? maxPhi + phiTol - 2 *
M_PI
923 : maxPhi + phiTol;
924 span.phiMax = maxPhi;
925 span.rMin = std::max(70.001, minR - zTol);
926 span.rMax = maxR + zTol;
927 span.xMin = minX - zTol;
928 span.xMax = maxX - +zTol;
929 span.yMin = minY - zTol;
930 span.yMax = maxY - +zTol;
931 } else {
933 }
934 return std::make_unique<VolumeSpan>(span);
935}
Scalar perp() const
perp method - perpendicular length
const PlainObject unit() const
This is a plugin that makes Eigen look like CLHEP & defines some convenience methods.
#define ATH_MSG_WARNING(x)
std::unique_ptr< VolumeSpan > findVolumeSpan(const VolumeBounds &volBounds, const Amg::Transform3D &transform, double zTol, double phiTol) const
Estimation of the geometrical volume span.
Amg::Vector3D transform(Amg::Vector3D &v, Amg::Transform3D &tr)
Transform a point from a Trasformation3D.
Eigen::Matrix< double, 3, 1 > Vector3D
span(T *ptr, std::size_t sz) -> span< T >
A couple needed deduction guides.
bool dPhi(const xAOD::TauJet &tau, const xAOD::CaloVertexedTopoCluster &cluster, float &out)