src/Surfaces/AnnulusBounds.cpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 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 https://mozilla.org/MPL/2.0/.
0008
0009 #include "Acts/Surfaces/AnnulusBounds.hpp"
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Surfaces/detail/VerticesHelper.hpp"
0013 #include "Acts/Utilities/VectorHelpers.hpp"
0014 #include "Acts/Utilities/detail/periodic.hpp"
0015
0016 #include <algorithm>
0017 #include <cmath>
0018 #include <iomanip>
0019 #include <iostream>
0020 #include <limits>
0021
0022 namespace Acts {
0023
0024 namespace {
0025
0026 Vector2 closestOnSegment(const Vector2& a, const Vector2& b, const Vector2& p,
0027                          const SquareMatrix2& metric) {
0028   // connecting vector
0029   auto n = b - a;
0030   // squared norm of line
0031   auto f = (n.transpose() * metric * n).value();
0032   // weighted scalar product of line to point and segment line
0033   auto u = ((p - a).transpose() * metric * n).value() / f;
0034   // clamp to [0, 1], convert to point
0035   return std::clamp(u, 0., 1.) * n + a;
0036 }
0037
0038 double squaredNorm(const Vector2& v, const SquareMatrix2& metric) {
0039   return (v.transpose() * metric * v).value();
0040 }
0041
0042 }  // namespace
0043
0044 AnnulusBounds::AnnulusBounds(const std::array<double, eSize>& values) noexcept(
0045     false)
0046     : m_values(values), m_moduleOrigin({values[eOriginX], values[eOriginY]}) {
0047   checkConsistency();
0048   m_rotationStripPC = Translation2(Vector2(0, -get(eAveragePhi)));
0049   m_translation = Translation2(m_moduleOrigin);
0050
0051   m_shiftXY = m_moduleOrigin * -1;
0052   m_shiftPC =
0053       Vector2(VectorHelpers::perp(m_shiftXY), VectorHelpers::phi(m_shiftXY));
0054
0055   // we need the corner points of the module to do the inside
0056   // checking, calculate them here once, they don't change
0057
0058   // find inner outer radius at edges in STRIP PC
0059   auto circIx = [](double O_x, double O_y, double r, double phi) -> Vector2 {
0060     //                      _____________________________________________
0061     //                     /      2  2                    2    2  2    2
0062     //     O_x + O_y*m - \/  - O_x *m  + 2*O_x*O_y*m - O_y  + m *r  + r
0063     // x = --------------------------------------------------------------
0064     //                                  2
0065     //                                 m  + 1
0066     //
0067     // y = m*x
0068     //
0069     double m = std::tan(phi);
0070     Vector2 dir(std::cos(phi), std::sin(phi));
0071     double x1 = (O_x + O_y * m -
0072                  std::sqrt(-std::pow(O_x, 2) * std::pow(m, 2) +
0073                            2 * O_x * O_y * m - std::pow(O_y, 2) +
0074                            std::pow(m, 2) * std::pow(r, 2) + std::pow(r, 2))) /
0075                 (std::pow(m, 2) + 1);
0076     double x2 = (O_x + O_y * m +
0077                  std::sqrt(-std::pow(O_x, 2) * std::pow(m, 2) +
0078                            2 * O_x * O_y * m - std::pow(O_y, 2) +
0079                            std::pow(m, 2) * std::pow(r, 2) + std::pow(r, 2))) /
0080                 (std::pow(m, 2) + 1);
0081
0082     Vector2 v1(x1, m * x1);
0083     if (v1.dot(dir) > 0) {
0084       return v1;
0085     }
0086     return {x2, m * x2};
0087   };
0088
0089   // calculate corners in STRIP XY, keep them we need them for minDistance()
0090   m_outLeftStripXY =
0091       circIx(m_moduleOrigin[0], m_moduleOrigin[1], get(eMaxR), get(eMaxPhiRel));
0092   m_inLeftStripXY =
0093       circIx(m_moduleOrigin[0], m_moduleOrigin[1], get(eMinR), get(eMaxPhiRel));
0094   m_outRightStripXY =
0095       circIx(m_moduleOrigin[0], m_moduleOrigin[1], get(eMaxR), get(eMinPhiRel));
0096   m_inRightStripXY =
0097       circIx(m_moduleOrigin[0], m_moduleOrigin[1], get(eMinR), get(eMinPhiRel));
0098
0099   m_outLeftStripPC = {m_outLeftStripXY.norm(),
0100                       VectorHelpers::phi(m_outLeftStripXY)};
0101   m_inLeftStripPC = {m_inLeftStripXY.norm(),
0102                      VectorHelpers::phi(m_inLeftStripXY)};
0103   m_outRightStripPC = {m_outRightStripXY.norm(),
0104                        VectorHelpers::phi(m_outRightStripXY)};
0105   m_inRightStripPC = {m_inRightStripXY.norm(),
0106                       VectorHelpers::phi(m_inRightStripXY)};
0107
0108   m_outLeftModulePC = stripXYToModulePC(m_outLeftStripXY);
0109   m_inLeftModulePC = stripXYToModulePC(m_inLeftStripXY);
0110   m_outRightModulePC = stripXYToModulePC(m_outRightStripXY);
0111   m_inRightModulePC = stripXYToModulePC(m_inRightStripXY);
0112
0113   // Pre-calculate center as average of corner vertices
0114   // Use simple 4-corner approximation for efficiency
0115   Vector2 centerXY = 0.25 * (m_outLeftStripXY + m_inLeftStripXY +
0116                              m_outRightStripXY + m_inRightStripXY);
0117
0118   // Convert center from cartesian (strip XY) to polar (strip PC) coordinates
0119   // Apply inverse rotation to account for average phi, matching other methods
0120   Vector2 centerPC = Vector2(centerXY.norm(), VectorHelpers::phi(centerXY));
0121   m_center = m_rotationStripPC.inverse() * centerPC;
0122 }
0123
0124 std::vector<double> AnnulusBounds::values() const {
0125   return {m_values.begin(), m_values.end()};
0126 }
0127
0128 void AnnulusBounds::checkConsistency() noexcept(false) {
0129   if (get(eMinR) < 0. || get(eMaxR) < 0. || get(eMinR) > get(eMaxR) ||
0130       std::abs(get(eMinR) - get(eMaxR)) < s_epsilon) {
0131     throw std::invalid_argument("AnnulusBounds: invalid radial setup.");
0132   }
0133   if (get(eMinPhiRel) != detail::radian_sym(get(eMinPhiRel)) ||
0134       get(eMaxPhiRel) != detail::radian_sym(get(eMaxPhiRel)) ||
0135       get(eMinPhiRel) > get(eMaxPhiRel)) {
0136     throw std::invalid_argument("AnnulusBounds: invalid phi boundary setup.");
0137   }
0138   if (get(eAveragePhi) != detail::radian_sym(get(eAveragePhi))) {
0139     throw std::invalid_argument("AnnulusBounds: invalid phi positioning.");
0140   }
0141 }
0142
0143 std::vector<Vector2> AnnulusBounds::corners() const {
0144   auto rot = m_rotationStripPC.inverse();
0145
0146   return {rot * m_outRightStripPC, rot * m_outLeftStripPC,
0147           rot * m_inLeftStripPC, rot * m_inRightStripPC};
0148 }
0149
0150 std::vector<Vector2> AnnulusBounds::vertices(
0151     unsigned int quarterSegments) const {
0152   if (quarterSegments > 0u) {
0153     using VectorHelpers::phi;
0154
0155     double phiMinInner = phi(m_inRightStripXY - m_moduleOrigin);
0156     double phiMaxInner = phi(m_inLeftStripXY - m_moduleOrigin);
0157
0158     double phiMinOuter = phi(m_outRightStripXY - m_moduleOrigin);
0159     double phiMaxOuter = phi(m_outLeftStripXY - m_moduleOrigin);
0160
0161     // Inner bow from phi_min -> phi_max (needs to be reversed)
0162     std::vector<Vector2> rvertices =
0163         detail::VerticesHelper::segmentVertices<Vector2, Transform2>(
0164             {get(eMinR), get(eMinR)}, phiMinInner, phiMaxInner, {},
0165             quarterSegments);
0166     std::reverse(rvertices.begin(), rvertices.end());
0167
0168     // Outer bow from phi_min -> phi_max
0169     auto overtices =
0170         detail::VerticesHelper::segmentVertices<Vector2, Transform2>(
0171             {get(eMaxR), get(eMaxR)}, phiMinOuter, phiMaxOuter, {},
0172             quarterSegments);
0173     rvertices.insert(rvertices.end(), overtices.begin(), overtices.end());
0174
0175     std::ranges::for_each(rvertices,
0176                           [&](Vector2& rv) { rv += m_moduleOrigin; });
0177     return rvertices;
0178   }
0179   return {m_inLeftStripXY, m_inRightStripXY, m_outRightStripXY,
0180           m_outLeftStripXY};
0181 }
0182
0183 SquareMatrix2 AnnulusBounds::boundToCartesianJacobian(
0184     const Vector2& lposition) const {
0185   SquareMatrix2 j;
0186   j(0, 0) = std::cos(lposition[1] - get(eAveragePhi));
0187   j(0, 1) = -lposition[0] * std::sin(lposition[1] - get(eAveragePhi));
0188   j(1, 0) = std::sin(lposition[1] - get(eAveragePhi));
0189   j(1, 1) = lposition[0] * std::cos(lposition[1] - get(eAveragePhi));
0190   return j;
0191 }
0192
0193 SquareMatrix2 AnnulusBounds::boundToCartesianMetric(
0194     const Vector2& lposition) const {
0195   SquareMatrix2 m;
0196   m(0, 0) = 1;
0197   m(0, 1) = 0;
0198   m(1, 0) = 0;
0199   m(1, 1) = lposition[0] * lposition[0];
0200   return m;
0201 }
0202
0203 bool AnnulusBounds::inside(const Vector2& lposition) const {
0204   // locpo is PC in STRIP SYSTEM
0205   // need to perform internal rotation induced by average phi
0206   Vector2 locpo_rotated = m_rotationStripPC * lposition;
0207   double phiLoc = locpo_rotated[1];
0208   double rLoc = locpo_rotated[0];
0209
0210   if (phiLoc < get(eMinPhiRel) || phiLoc > get(eMaxPhiRel)) {
0211     return false;
0212   }
0213
0214   // calculate R in MODULE SYSTEM to evaluate R-bounds
0215   // don't need R, can use R^2
0216   double r_mod2 = m_shiftPC[0] * m_shiftPC[0] + rLoc * rLoc +
0217                   2 * m_shiftPC[0] * rLoc * cos(phiLoc - m_shiftPC[1]);
0218
0219   if (r_mod2 < get(eMinR) * get(eMinR) || r_mod2 > get(eMaxR) * get(eMaxR)) {
0220     return false;
0221   }
0222
0223   return true;
0224 }
0225
0226 SquareMatrix2 AnnulusBounds::stripPCToModulePCJacobian(
0227     const Vector2& lpositionRotated) const {
0228   double dphi = get(eAveragePhi);
0229   double phi_strip = lpositionRotated[1];
0230   double r_strip = lpositionRotated[0];
0231   double O_x = m_shiftXY[0];
0232   double O_y = m_shiftXY[1];
0233
0234   // For a transformation from cartesian into polar coordinates
0235   //
0236   //              [         _________      ]
0237   //              [        /  2    2       ]
0238   //              [      \/  x  + y        ]
0239   //     [ r' ]   [                        ]
0240   // v = [    ] = [      /       y        \]
0241   //     [phi']   [2*atan|----------------|]
0242   //              [      |       _________|]
0243   //              [      |      /  2    2 |]
0244   //              [      \x + \/  x  + y  /]
0245   //
0246   // Where x, y are polar coordinates that can be rotated by dPhi
0247   //
0248   // [x]   [O_x + r*cos(dPhi - phi)]
0249   // [ ] = [                       ]
0250   // [y]   [O_y - r*sin(dPhi - phi)]
0251   //
0252   // The general jacobian is:
0253   //
0254   //        [d        d      ]
0255   //        [--(f_x)  --(f_x)]
0256   //        [dx       dy     ]
0257   // Jgen = [                ]
0258   //        [d        d      ]
0259   //        [--(f_y)  --(f_y)]
0260   //        [dx       dy     ]
0261   //
0262   // which means in this case:
0263   //
0264   //     [     d                   d           ]
0265   //     [ ----------(rMod)    ---------(rMod) ]
0266   //     [ dr_{strip}          dphiStrip       ]
0267   // J = [                                     ]
0268   //     [    d                   d            ]
0269   //     [----------(phiMod)  ---------(phiMod)]
0270   //     [dr_{strip}          dphiStrip        ]
0271   //
0272   // Performing the derivative one gets:
0273   //
0274   //     [B*O_x + C*O_y + rStrip  rStrip*(B*O_y + O_x*sin(dPhi - phiStrip))]
0275   //     [----------------------  -----------------------------------------]
0276   //     [          ___                               ___                  ]
0277   //     [        \/ A                              \/ A                   ]
0278   // J = [                                                                 ]
0279   //     [  -(B*O_y - C*O_x)           rStrip*(B*O_x + C*O_y + rStrip)     ]
0280   //     [  -----------------          -------------------------------     ]
0281   //     [          A                                 A                    ]
0282   //
0283   // where
0284   //        2                                          2
0285   // A = O_x  + 2*O_x*rStrip*cos(dPhi - phiStrip) + O_y
0286   //                                                 2
0287   //     - 2*O_y*rStrip*sin(dPhi - phiStrip) + rStrip
0288   // B = cos(dPhi - phiStrip)
0289   // C = -sin(dPhi - phiStrip)
0290
0291   double cosDPhiPhiStrip = std::cos(dphi - phi_strip);
0292   double sinDPhiPhiStrip = std::sin(dphi - phi_strip);
0293
0294   double A = O_x * O_x + 2 * O_x * r_strip * cosDPhiPhiStrip + O_y * O_y -
0295              2 * O_y * r_strip * sinDPhiPhiStrip + r_strip * r_strip;
0296   double sqrtA = std::sqrt(A);
0297
0298   double B = cosDPhiPhiStrip;
0299   double C = -sinDPhiPhiStrip;
0300
0301   SquareMatrix2 j;
0302   j(0, 0) = (B * O_x + C * O_y + r_strip) / sqrtA;
0303   j(0, 1) = r_strip * (B * O_y + O_x * sinDPhiPhiStrip) / sqrtA;
0304   j(1, 0) = -(B * O_y - C * O_x) / A;
0305   j(1, 1) = r_strip * (B * O_x + C * O_y + r_strip) / A;
0306   return j;
0307 }
0308
0309 Vector2 AnnulusBounds::closestPoint(const Vector2& lposition,
0310                                     const SquareMatrix2& metric) const {
0311   // lposition is PC in STRIP SYSTEM
0312   // we need to rotate the local position
0313   Vector2 lpositionRotated = m_rotationStripPC * lposition;
0314
0315   SquareMatrix2 jacobianStripPCToModulePC =
0316       stripPCToModulePCJacobian(lpositionRotated);
0317
0318   // calculate the metrics for STRIP PC and MODULE PC
0319   SquareMatrix2 metricStripPC = metric;
0320   SquareMatrix2 metricModulePC = jacobianStripPCToModulePC.transpose() *
0321                                  metricStripPC * jacobianStripPCToModulePC;
0322
0323   // minimum distance and associated point
0324   double minDist = std::numeric_limits<double>::max();
0325   Vector2 closest = Vector2::Zero();
0326
0327   // first: STRIP system. lpositionRotated is in STRIP PC already
0328
0329   {
0330     Vector2 currentClosest = closestOnSegment(m_inLeftStripPC, m_outLeftStripPC,
0331                                               lpositionRotated, metricStripPC);
0332     double currentDist =
0333         squaredNorm(lpositionRotated - currentClosest, metricStripPC);
0334     if (currentDist < minDist) {
0335       minDist = currentDist;
0336       closest = m_rotationStripPC.inverse() * currentClosest;
0337     }
0338   }
0339
0340   {
0341     Vector2 currentClosest = closestOnSegment(
0342         m_inRightStripPC, m_outRightStripPC, lpositionRotated, metricStripPC);
0343     double currentDist =
0344         squaredNorm(lpositionRotated - currentClosest, metricStripPC);
0345     if (currentDist < minDist) {
0346       minDist = currentDist;
0347       closest = m_rotationStripPC.inverse() * currentClosest;
0348     }
0349   }
0350
0351   // now: MODULE system. Need to transform lposition to MODULE PC
0352   //  transform is STRIP PC -> STRIP XY -> MODULE XY -> MODULE PC
0353   Vector2 lpositionModulePC = stripPCToModulePC(lpositionRotated);
0354
0355   // now check edges in MODULE PC (inner and outer circle)
0356
0357   {
0358     Vector2 currentClosest = closestOnSegment(
0359         m_inLeftModulePC, m_inRightModulePC, lpositionModulePC, metricModulePC);
0360     double currentDist =
0361         squaredNorm(lpositionModulePC - currentClosest, metricModulePC);
0362     if (currentDist < minDist) {
0363       minDist = currentDist;
0364       closest = m_rotationStripPC.inverse() * modulePCToStripPC(currentClosest);
0365     }
0366   }
0367
0368   {
0369     Vector2 currentClosest =
0370         closestOnSegment(m_outLeftModulePC, m_outRightModulePC,
0371                          lpositionModulePC, metricModulePC);
0372     double currentDist =
0373         squaredNorm(lpositionModulePC - currentClosest, metricModulePC);
0374     if (currentDist < minDist) {
0375       minDist = currentDist;
0376       closest = m_rotationStripPC.inverse() * modulePCToStripPC(currentClosest);
0377     }
0378   }
0379
0380   return closest;
0381 }
0382
0383 Vector2 AnnulusBounds::stripXYToModulePC(const Vector2& vStripXY) const {
0384   Vector2 vModuleXY = vStripXY + m_shiftXY;
0385   return {vModuleXY.norm(), VectorHelpers::phi(vModuleXY)};
0386 }
0387
0388 Vector2 AnnulusBounds::stripPCToModulePC(const Vector2& vStripPC) const {
0389   Vector2 vStripXY = {vStripPC[0] * std::cos(vStripPC[1]),
0390                       vStripPC[0] * std::sin(vStripPC[1])};
0391   return stripXYToModulePC(vStripXY);
0392 }
0393
0394 Vector2 AnnulusBounds::modulePCToStripPC(const Vector2& vModulePC) const {
0395   Vector2 vModuleXY = {vModulePC[0] * std::cos(vModulePC[1]),
0396                        vModulePC[0] * std::sin(vModulePC[1])};
0397   Vector2 vStripXY = vModuleXY - m_shiftXY;
0398   return {vStripXY.norm(), VectorHelpers::phi(vStripXY)};
0399 }
0400
0401 Vector2 AnnulusBounds::moduleOrigin() const {
0402   return Eigen::Rotation2D<double>(get(eAveragePhi)) * m_moduleOrigin;
0403 }
0404
0405 Vector2 AnnulusBounds::center() const {
0406   // Return pre-calculated center that accounts for the complex coordinate
0407   // transformations between radial and angular coordinate systems
0408   return m_center;
0409 }
0410
0411 std::ostream& AnnulusBounds::toStream(std::ostream& sl) const {
0412   sl << std::setiosflags(std::ios::fixed);
0413   sl << std::setprecision(7);
0414   sl << "Acts::AnnulusBounds:  (innerRadius, outerRadius, minPhi, maxPhi) = ";
0415   sl << "(" << get(eMinR) << ", " << get(eMaxR) << ", " << phiMin() << ", "
0416      << phiMax() << ")" << '\n';
0417   sl << " - shift xy = " << m_shiftXY.x() << ", " << m_shiftXY.y() << '\n';
0418   sl << " - shift pc = " << m_shiftPC.x() << ", " << m_shiftPC.y() << '\n';
0419   sl << std::setprecision(-1);
0420   return sl;
0421 }
0422
0423 }  // namespace Acts