Back to home page

sPhenix code displayed by LXR

 
 

    


File indexing completed on 2025-08-05 08:09:28

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2018-2019 CERN for the benefit of the Acts project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
0008 
0009 #include "Acts/Utilities/VectorHelpers.hpp"
0010 
0011 template <typename value_t, std::size_t DIM, std::size_t SIDES>
0012 template <std::size_t D, std::enable_if_t<D == 2, int>>
0013 Acts::Frustum<value_t, DIM, SIDES>::Frustum(const VertexType& origin,
0014                                             const VertexType& dir,
0015                                             value_type opening_angle)
0016     : m_origin(origin) {
0017   using rotation_t = Eigen::Rotation2D<value_type>;
0018 
0019   static_assert(SIDES == 2, "2D frustum can only have 2 sides");
0020   assert(opening_angle < M_PI);
0021 
0022   translation_t translation(origin);
0023   value_type angle = VectorHelpers::phi(dir);
0024   Eigen::Rotation2D<value_type> rot(angle);
0025 
0026   value_type normal_angle = 0.5 * M_PI - 0.5 * opening_angle;
0027   VertexType normal1 = rotation_t(normal_angle) * VertexType::UnitX();
0028   VertexType normal2 = rotation_t(-normal_angle) * VertexType::UnitX();
0029 
0030   m_normals = {rot * VertexType::UnitX(), rot * normal1, rot * normal2};
0031 }
0032 
0033 template <typename value_t, std::size_t DIM, std::size_t SIDES>
0034 template <std::size_t D, std::enable_if_t<D == 3, int>>
0035 Acts::Frustum<value_t, DIM, SIDES>::Frustum(const VertexType& origin,
0036                                             const VertexType& dir,
0037                                             value_type opening_angle)
0038     : m_origin(origin) {
0039   static_assert(SIDES > 2, "3D frustum must have 3 or more sides");
0040   assert(opening_angle < M_PI);
0041   using angle_axis_t = Eigen::AngleAxis<value_type>;
0042 
0043   const VertexType ldir = VertexType::UnitZ();
0044   const VertexType lup = VertexType::UnitX();
0045 
0046   transform_type transform;
0047   transform = (Eigen::Quaternion<value_type>().setFromTwoVectors(ldir, dir));
0048 
0049   m_normals[0] = ldir;
0050 
0051   const value_type phi_sep = 2 * M_PI / sides;
0052   transform_type rot;
0053   rot = angle_axis_t(phi_sep, ldir);
0054 
0055   value_type half_opening_angle = opening_angle / 2.;
0056   auto calculate_normal =
0057       [&ldir, &half_opening_angle](const VertexType& out) -> VertexType {
0058     const VertexType tilt_axis = -1 * out.cross(ldir);
0059     return (-1 * (angle_axis_t(half_opening_angle, tilt_axis) * out))
0060         .normalized();
0061   };
0062 
0063   VertexType current_outward = lup;
0064   m_normals[1] = calculate_normal(current_outward);
0065 
0066   for (std::size_t i = 1; i < sides; i++) {
0067     current_outward = rot * current_outward;
0068     m_normals[i + 1] = calculate_normal(current_outward);
0069   }
0070 
0071   for (auto& normal : m_normals) {
0072     normal = transform * normal;
0073   }
0074 }
0075 
0076 template <typename value_t, std::size_t DIM, std::size_t SIDES>
0077 template <std::size_t D, std::enable_if_t<D == 3, int>>
0078 void Acts::Frustum<value_t, DIM, SIDES>::draw(IVisualization3D& helper,
0079                                               value_type far_distance) const {
0080   static_assert(DIM == 3, "Drawing is only supported in 3D");
0081 
0082   // Iterate around normals, calculate cross with "far" plane
0083   // to get intersection lines.
0084   // Work in local reference frame of the frustum, and only convert to global
0085   // right before drawing.
0086   VertexType far_normal = m_normals[0];  // far has same normal as pseudo-near
0087   VertexType far_center = m_normals[0] * far_distance;
0088   std::array<std::pair<VertexType, VertexType>, SIDES> planeFarIXs;
0089 
0090   auto ixPlanePlane = [](const auto& n1, const auto& p1, const auto& n2,
0091                          const auto& p2) -> std::pair<VertexType, VertexType> {
0092     const VertexType m = n1.cross(n2).normalized();
0093     const double j = (n2.dot(p2 - p1)) / (n2.dot(n1.cross(m)));
0094     const VertexType q = p1 + j * n1.cross(m);
0095     return {m, q};
0096   };
0097 
0098   auto ixLineLine = [](const auto& p1, const auto& d1, const auto& p2,
0099                        const auto& d2) -> VertexType {
0100     return p1 + (((p2 - p1).cross(d2)).norm() / (d1.cross(d2)).norm()) * d1;
0101   };
0102 
0103   // skip i=0 <=> pseudo-near
0104   for (std::size_t i = 1; i < n_normals; i++) {
0105     const auto ixLine =
0106         ixPlanePlane(far_normal, far_center, m_normals[i], VertexType::Zero());
0107     planeFarIXs.at(i - 1) = ixLine;
0108   }
0109 
0110   std::array<VertexType, SIDES> points;
0111 
0112   for (std::size_t i = 0; i < std::size(planeFarIXs); i++) {
0113     std::size_t j = (i + 1) % std::size(planeFarIXs);
0114     const auto& l1 = planeFarIXs.at(i);
0115     const auto& l2 = planeFarIXs.at(j);
0116     const VertexType ix =
0117         m_origin + ixLineLine(l1.second, l1.first, l2.second, l2.first);
0118     points.at(i) = ix;
0119   }
0120 
0121   for (std::size_t i = 0; i < std::size(points); i++) {
0122     std::size_t j = (i + 1) % std::size(points);
0123     helper.face(
0124         std::vector<VertexType>({m_origin, points.at(i), points.at(j)}));
0125   }
0126 }
0127 
0128 template <typename value_t, std::size_t DIM, std::size_t SIDES>
0129 template <std::size_t D, std::enable_if_t<D == 2, int>>
0130 std::ostream& Acts::Frustum<value_t, DIM, SIDES>::svg(std::ostream& os,
0131                                                       value_type w,
0132                                                       value_type h,
0133                                                       value_type far_distance,
0134                                                       value_type unit) const {
0135   static_assert(DIM == 2, "SVG is only supported in 2D");
0136 
0137   VertexType mid(w / 2., h / 2.);
0138 
0139   // set up transform for svg. +y is down, normally, and unit is pixels.
0140   // We flip the y axis, and scale up by `unit`.
0141   transform_type trf = transform_type::Identity();
0142   trf.translate(mid);
0143   trf = trf * Eigen::Scaling(VertexType(1, -1));
0144   trf.scale(unit);
0145 
0146   std::array<std::string, 3> colors({"orange", "blue", "red"});
0147 
0148   auto draw_line = [&](const VertexType& left_, const VertexType& right_,
0149                        const std::string& color, std::size_t width) {
0150     VertexType left = trf * left_;
0151     VertexType right = trf * right_;
0152     os << "<line ";
0153 
0154     os << "x1=\"" << left.x() << "\" ";
0155     os << "y1=\"" << left.y() << "\" ";
0156     os << "x2=\"" << right.x() << "\" ";
0157     os << "y2=\"" << right.y() << "\" ";
0158 
0159     os << " stroke=\"" << color << "\" stroke-width=\"" << width << "\"/>\n";
0160   };
0161 
0162   auto draw_point = [&](const VertexType& p_, const std::string& color,
0163                         std::size_t r) {
0164     VertexType p = trf * p_;
0165     os << "<circle ";
0166     os << "cx=\"" << p.x() << "\" cy=\"" << p.y() << "\" r=\"" << r << "\"";
0167     os << " fill=\"" << color << "\"";
0168     os << "/>\n";
0169   };
0170 
0171   using vec3 = Eigen::Matrix<value_type, 3, 1>;
0172   auto ixLineLine = [](const VertexType& p1_2, const VertexType& d1_2,
0173                        const VertexType& p2_2,
0174                        const VertexType& d2_2) -> VertexType {
0175     const vec3 p1(p1_2.x(), p1_2.y(), 0);
0176     const vec3 p2(p2_2.x(), p2_2.y(), 0);
0177     const vec3 d1(d1_2.x(), d1_2.y(), 0);
0178     const vec3 d2(d2_2.x(), d2_2.y(), 0);
0179 
0180     vec3 num = (p2 - p1).cross(d2);
0181     vec3 den = d1.cross(d2);
0182 
0183     value_type f = 1.;
0184     value_type dot = num.normalized().dot(den.normalized());
0185     if (std::abs(dot) - 1 < 1e-9 && dot < 0) {
0186       f = -1.;
0187     }
0188 
0189     const vec3 p = p1 + f * (num.norm() / den.norm()) * d1;
0190     assert(std::abs(p.z()) < 1e-9);
0191     return {p.x(), p.y()};
0192   };
0193 
0194   const VertexType far_dir = {m_normals[0].y(), -m_normals[0].x()};
0195   const VertexType far_point = m_normals[0] * far_distance;
0196 
0197   std::array<VertexType, 2> points;
0198 
0199   for (std::size_t i = 1; i < n_normals; i++) {
0200     VertexType plane_dir(m_normals[i].y(), -m_normals[i].x());
0201 
0202     const VertexType ix = ixLineLine(far_point, far_dir, {0, 0}, plane_dir);
0203     draw_point(m_origin + ix, "black", 3);
0204     draw_line(m_origin, m_origin + ix, "black", 2);
0205     points.at(i - 1) = ix;
0206   }
0207 
0208   draw_line(m_origin + points.at(0), m_origin + points.at(1), "black", 2);
0209 
0210   draw_line(m_origin, m_origin + m_normals[0] * 2, colors[0], 3);
0211 
0212   draw_point({0, 0}, "yellow", 5);
0213   draw_point(m_origin, "green", 5);
0214 
0215   return os;
0216 }
0217 
0218 template <typename value_t, std::size_t DIM, std::size_t SIDES>
0219 Acts::Frustum<value_t, DIM, SIDES>
0220 Acts::Frustum<value_t, DIM, SIDES>::transformed(
0221     const transform_type& trf) const {
0222   const auto& rot = trf.rotation();
0223 
0224   std::array<VertexType, n_normals> new_normals;
0225   for (std::size_t i = 0; i < n_normals; i++) {
0226     new_normals[i] = rot * m_normals[i];
0227   }
0228 
0229   return Frustum<value_t, DIM, SIDES>(trf * m_origin, std::move(new_normals));
0230 }