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