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