src/Surfaces/AnnulusBounds.cpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2020 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/Surfaces/AnnulusBounds.hpp"
0010
0011 #include "Acts/Definitions/TrackParametrization.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 Acts::AnnulusBounds::AnnulusBounds(
0023     const std::array<double, eSize>& values) noexcept(false)
0024     : m_values(values), m_moduleOrigin({values[eOriginX], values[eOriginY]}) {
0025   checkConsistency();
0026   m_rotationStripPC = Translation2(Vector2(0, -get(eAveragePhi)));
0027   m_translation = Translation2(m_moduleOrigin);
0028
0029   m_shiftXY = m_moduleOrigin * -1;
0030   m_shiftPC =
0031       Vector2(VectorHelpers::perp(m_shiftXY), VectorHelpers::phi(m_shiftXY));
0032
0033   // we need the corner points of the module to do the inside
0034   // checking, calculate them here once, they don't change
0035
0036   // find inner outer radius at edges in STRIP PC
0037   auto circIx = [](double O_x, double O_y, double r, double phi) -> Vector2 {
0038     //                      _____________________________________________
0039     //                     /      2  2                    2    2  2    2
0040     //     O_x + O_y*m - \/  - O_x *m  + 2*O_x*O_y*m - O_y  + m *r  + r
0041     // x = --------------------------------------------------------------
0042     //                                  2
0043     //                                 m  + 1
0044     //
0045     // y = m*x
0046     //
0047     double m = std::tan(phi);
0048     Vector2 dir(std::cos(phi), std::sin(phi));
0049     double x1 = (O_x + O_y * m -
0050                  std::sqrt(-std::pow(O_x, 2) * std::pow(m, 2) +
0051                            2 * O_x * O_y * m - std::pow(O_y, 2) +
0052                            std::pow(m, 2) * std::pow(r, 2) + std::pow(r, 2))) /
0053                 (std::pow(m, 2) + 1);
0054     double x2 = (O_x + O_y * m +
0055                  std::sqrt(-std::pow(O_x, 2) * std::pow(m, 2) +
0056                            2 * O_x * O_y * m - std::pow(O_y, 2) +
0057                            std::pow(m, 2) * std::pow(r, 2) + std::pow(r, 2))) /
0058                 (std::pow(m, 2) + 1);
0059
0060     Vector2 v1(x1, m * x1);
0061     if (v1.dot(dir) > 0) {
0062       return v1;
0063     }
0064     return {x2, m * x2};
0065   };
0066
0067   // calculate corners in STRIP XY, keep them we need them for minDistance()
0068   m_outLeftStripXY =
0069       circIx(m_moduleOrigin[eBoundLoc0], m_moduleOrigin[eBoundLoc1], get(eMaxR),
0070              get(eMaxPhiRel));
0071   m_inLeftStripXY =
0072       circIx(m_moduleOrigin[eBoundLoc0], m_moduleOrigin[eBoundLoc1], get(eMinR),
0073              get(eMaxPhiRel));
0074   m_outRightStripXY =
0075       circIx(m_moduleOrigin[eBoundLoc0], m_moduleOrigin[eBoundLoc1], get(eMaxR),
0076              get(eMinPhiRel));
0077   m_inRightStripXY =
0078       circIx(m_moduleOrigin[eBoundLoc0], m_moduleOrigin[eBoundLoc1], get(eMinR),
0079              get(eMinPhiRel));
0080
0081   m_outLeftStripPC = {m_outLeftStripXY.norm(),
0082                       VectorHelpers::phi(m_outLeftStripXY)};
0083   m_inLeftStripPC = {m_inLeftStripXY.norm(),
0084                      VectorHelpers::phi(m_inLeftStripXY)};
0085   m_outRightStripPC = {m_outRightStripXY.norm(),
0086                        VectorHelpers::phi(m_outRightStripXY)};
0087   m_inRightStripPC = {m_inRightStripXY.norm(),
0088                       VectorHelpers::phi(m_inRightStripXY)};
0089
0090   m_outLeftModulePC = stripXYToModulePC(m_outLeftStripXY);
0091   m_inLeftModulePC = stripXYToModulePC(m_inLeftStripXY);
0092   m_outRightModulePC = stripXYToModulePC(m_outRightStripXY);
0093   m_inRightModulePC = stripXYToModulePC(m_inRightStripXY);
0094 }
0095
0096 std::vector<Acts::Vector2> Acts::AnnulusBounds::corners() const {
0097   auto rot = m_rotationStripPC.inverse();
0098
0099   return {rot * m_outRightStripPC, rot * m_outLeftStripPC,
0100           rot * m_inLeftStripPC, rot * m_inRightStripPC};
0101 }
0102
0103 std::vector<Acts::Vector2> Acts::AnnulusBounds::vertices(
0104     unsigned int lseg) const {
0105   if (lseg > 0) {
0106     // List of vertices counter-clockwise starting with left inner
0107     std::vector<Acts::Vector2> rvertices;
0108
0109     using VectorHelpers::phi;
0110     auto phisInner = detail::VerticesHelper::phiSegments(
0111         phi(m_inRightStripXY - m_moduleOrigin),
0112         phi(m_inLeftStripXY - m_moduleOrigin));
0113     auto phisOuter = detail::VerticesHelper::phiSegments(
0114         phi(m_outLeftStripXY - m_moduleOrigin),
0115         phi(m_outRightStripXY - m_moduleOrigin));
0116
0117     // Inner bow from phi_min -> phi_max
0118     for (unsigned int iseg = 0; iseg < phisInner.size() - 1; ++iseg) {
0119       int addon = (iseg == phisInner.size() - 2) ? 1 : 0;
0120       detail::VerticesHelper::createSegment<Vector2, Transform2>(
0121           rvertices, {get(eMinR), get(eMinR)}, phisInner[iseg],
0122           phisInner[iseg + 1], lseg, addon);
0123     }
0124     // Upper bow from phi_max -> phi_min
0125     for (unsigned int iseg = 0; iseg < phisOuter.size() - 1; ++iseg) {
0126       int addon = (iseg == phisOuter.size() - 2) ? 1 : 0;
0127       detail::VerticesHelper::createSegment<Vector2, Transform2>(
0128           rvertices, {get(eMaxR), get(eMaxR)}, phisOuter[iseg],
0129           phisOuter[iseg + 1], lseg, addon);
0130     }
0131     std::for_each(rvertices.begin(), rvertices.end(),
0132                   [&](Acts::Vector2& rv) { rv += m_moduleOrigin; });
0133     return rvertices;
0134   }
0135   return {m_inLeftStripXY, m_inRightStripXY, m_outRightStripXY,
0136           m_outLeftStripXY};
0137 }
0138
0139 bool Acts::AnnulusBounds::inside(const Vector2& lposition, double tolR,
0140                                  double tolPhi) const {
0141   // locpo is PC in STRIP SYSTEM
0142   // need to perform internal rotation induced by average phi
0143   Vector2 locpo_rotated = m_rotationStripPC * lposition;
0144   double phiLoc = locpo_rotated[eBoundLoc1];
0145   double rLoc = locpo_rotated[eBoundLoc0];
0146
0147   if (phiLoc < (get(eMinPhiRel) - tolPhi) ||
0148       phiLoc > (get(eMaxPhiRel) + tolPhi)) {
0149     return false;
0150   }
0151
0152   // calculate R in MODULE SYSTEM to evaluate R-bounds
0153   if (tolR == 0.) {
0154     // don't need R, can use R^2
0155     double r_mod2 =
0156         m_shiftPC[eBoundLoc0] * m_shiftPC[eBoundLoc0] + rLoc * rLoc +
0157         2 * m_shiftPC[eBoundLoc0] * rLoc * cos(phiLoc - m_shiftPC[eBoundLoc1]);
0158
0159     if (r_mod2 < get(eMinR) * get(eMinR) || r_mod2 > get(eMaxR) * get(eMaxR)) {
0160       return false;
0161     }
0162   } else {
0163     // use R
0164     double r_mod = sqrt(
0165         m_shiftPC[eBoundLoc0] * m_shiftPC[eBoundLoc0] + rLoc * rLoc +
0166         2 * m_shiftPC[eBoundLoc0] * rLoc * cos(phiLoc - m_shiftPC[eBoundLoc1]));
0167
0168     if (r_mod < (get(eMinR) - tolR) || r_mod > (get(eMaxR) + tolR)) {
0169       return false;
0170     }
0171   }
0172   return true;
0173 }
0174
0175 bool Acts::AnnulusBounds::inside(const Vector2& lposition,
0176                                  const BoundaryCheck& bcheck) const {
0177   // locpo is PC in STRIP SYSTEM
0178   if (bcheck.type() == BoundaryCheck::Type::eAbsolute) {
0179     return inside(lposition, bcheck.tolerance()[eBoundLoc0],
0180                   bcheck.tolerance()[eBoundLoc1]);
0181   } else {
0182     // first check if inside. We don't need to look into the covariance if
0183     // inside
0184     if (inside(lposition, 0., 0.)) {
0185       return true;
0186     }
0187
0188     // we need to rotate the locpo
0189     Vector2 locpo_rotated = m_rotationStripPC * lposition;
0190
0191     // covariance is given in STRIP SYSTEM in PC
0192     // we need to convert the covariance to the MODULE SYSTEM in PC
0193     // via jacobian.
0194     // The following transforms into STRIP XY, does the shift into MODULE XY,
0195     // and then transforms into MODULE PC
0196     double dphi = get(eAveragePhi);
0197     double phi_strip = locpo_rotated[eBoundLoc1];
0198     double r_strip = locpo_rotated[eBoundLoc0];
0199     double O_x = m_shiftXY[eBoundLoc0];
0200     double O_y = m_shiftXY[eBoundLoc1];
0201
0202     // For a transformation from cartesian into polar coordinates
0203     //
0204     //              [         _________      ]
0205     //              [        /  2    2       ]
0206     //              [      \/  x  + y        ]
0207     //     [ r' ]   [                        ]
0208     // v = [    ] = [      /       y        \]
0209     //     [phi']   [2*atan|----------------|]
0210     //              [      |       _________|]
0211     //              [      |      /  2    2 |]
0212     //              [      \x + \/  x  + y  /]
0213     //
0214     // Where x, y are polar coordinates that can be rotated by dPhi
0215     //
0216     // [x]   [O_x + r*cos(dPhi - phi)]
0217     // [ ] = [                       ]
0218     // [y]   [O_y - r*sin(dPhi - phi)]
0219     //
0220     // The general jacobian is:
0221     //
0222     //        [d        d      ]
0223     //        [--(f_x)  --(f_x)]
0224     //        [dx       dy     ]
0225     // Jgen = [                ]
0226     //        [d        d      ]
0227     //        [--(f_y)  --(f_y)]
0228     //        [dx       dy     ]
0229     //
0230     // which means in this case:
0231     //
0232     //     [     d                   d           ]
0233     //     [ ----------(rMod)    ---------(rMod) ]
0234     //     [ dr_{strip}          dphiStrip       ]
0235     // J = [                                     ]
0236     //     [    d                   d            ]
0237     //     [----------(phiMod)  ---------(phiMod)]
0238     //     [dr_{strip}          dphiStrip        ]
0239     //
0240     // Performing the derivative one gets:
0241     //
0242     //     [B*O_x + C*O_y + rStrip  rStrip*(B*O_y + O_x*sin(dPhi - phiStrip))]
0243     //     [----------------------  -----------------------------------------]
0244     //     [          ___                               ___                  ]
0245     //     [        \/ A                              \/ A                   ]
0246     // J = [                                                                 ]
0247     //     [  -(B*O_y - C*O_x)           rStrip*(B*O_x + C*O_y + rStrip)     ]
0248     //     [  -----------------          -------------------------------     ]
0249     //     [          A                                 A                    ]
0250     //
0251     // where
0252     //        2                                          2 2
0253     // A = O_x  + 2*O_x*rStrip*cos(dPhi - phiStrip) + O_y  -
0254     // 2*O_y*rStrip*sin(dPhi - phiStrip) + rStrip B = cos(dPhi - phiStrip) C =
0255     // -sin(dPhi - phiStrip)
0256
0257     double cosDPhiPhiStrip = std::cos(dphi - phi_strip);
0258     double sinDPhiPhiStrip = std::sin(dphi - phi_strip);
0259
0260     double A = O_x * O_x + 2 * O_x * r_strip * cosDPhiPhiStrip + O_y * O_y -
0261                2 * O_y * r_strip * sinDPhiPhiStrip + r_strip * r_strip;
0262     double sqrtA = std::sqrt(A);
0263
0264     double B = cosDPhiPhiStrip;
0265     double C = -sinDPhiPhiStrip;
0266     ActsMatrix<2, 2> jacobianStripPCToModulePC;
0267     jacobianStripPCToModulePC(0, 0) = (B * O_x + C * O_y + r_strip) / sqrtA;
0268     jacobianStripPCToModulePC(0, 1) =
0269         r_strip * (B * O_y + O_x * sinDPhiPhiStrip) / sqrtA;
0270     jacobianStripPCToModulePC(1, 0) = -(B * O_y - C * O_x) / A;
0271     jacobianStripPCToModulePC(1, 1) =
0272         r_strip * (B * O_x + C * O_y + r_strip) / A;
0273
0274     // covariance is given in STRIP PC
0275     auto covStripPC = bcheck.covariance();
0276     // calculate covariance in MODULE PC using jacobian from above
0277     auto covModulePC = jacobianStripPCToModulePC * covStripPC *
0278                        jacobianStripPCToModulePC.transpose();
0279
0280     // Mahalanobis distance uses inverse covariance as weights
0281     auto weightStripPC = covStripPC.inverse();
0282     auto weightModulePC = covModulePC.inverse();
0283
0284     double minDist = std::numeric_limits<double>::max();
0285
0286     Vector2 currentClosest;
0287     double currentDist = 0;
0288
0289     // do projection in STRIP PC
0290
0291     // first: STRIP system. locpo is in STRIP PC already
0292     currentClosest = closestOnSegment(m_inLeftStripPC, m_outLeftStripPC,
0293                                       locpo_rotated, weightStripPC);
0294     currentDist = squaredNorm(locpo_rotated - currentClosest, weightStripPC);
0295     minDist = currentDist;
0296
0297     currentClosest = closestOnSegment(m_inRightStripPC, m_outRightStripPC,
0298                                       locpo_rotated, weightStripPC);
0299     currentDist = squaredNorm(locpo_rotated - currentClosest, weightStripPC);
0300     if (currentDist < minDist) {
0301       minDist = currentDist;
0302     }
0303
0304     // now: MODULE system. Need to transform locpo to MODULE PC
0305     //  transform is STRIP PC -> STRIP XY -> MODULE XY -> MODULE PC
0306     Vector2 locpoStripXY(
0307         locpo_rotated[eBoundLoc0] * std::cos(locpo_rotated[eBoundLoc1]),
0308         locpo_rotated[eBoundLoc0] * std::sin(locpo_rotated[eBoundLoc1]));
0309     Vector2 locpoModulePC = stripXYToModulePC(locpoStripXY);
0310
0311     // now check edges in MODULE PC (inner and outer circle)
0312     // assuming Mahalanobis distances are of same unit if covariance
0313     // is correctly transformed
0314     currentClosest = closestOnSegment(m_inLeftModulePC, m_inRightModulePC,
0315                                       locpoModulePC, weightModulePC);
0316     currentDist = squaredNorm(locpoModulePC - currentClosest, weightModulePC);
0317     if (currentDist < minDist) {
0318       minDist = currentDist;
0319     }
0320
0321     currentClosest = closestOnSegment(m_outLeftModulePC, m_outRightModulePC,
0322                                       locpoModulePC, weightModulePC);
0323     currentDist = squaredNorm(locpoModulePC - currentClosest, weightModulePC);
0324     if (currentDist < minDist) {
0325       minDist = currentDist;
0326     }
0327
0328     // compare resulting Mahalanobis distance to configured
0329     // "number of sigmas"
0330     // we square it b/c we never took the square root of the distance
0331     return minDist < bcheck.tolerance()[0] * bcheck.tolerance()[0];
0332   }
0333 }
0334
0335 Acts::Vector2 Acts::AnnulusBounds::stripXYToModulePC(
0336     const Vector2& vStripXY) const {
0337   Vector2 vecModuleXY = vStripXY + m_shiftXY;
0338   return {vecModuleXY.norm(), VectorHelpers::phi(vecModuleXY)};
0339 }
0340
0341 Acts::Vector2 Acts::AnnulusBounds::closestOnSegment(
0342     const Vector2& a, const Vector2& b, const Vector2& p,
0343     const SquareMatrix2& weight) const {
0344   // connecting vector
0345   auto n = b - a;
0346   // squared norm of line
0347   auto f = (n.transpose() * weight * n).value();
0348   // weighted scalar product of line to point and segment line
0349   auto u = ((p - a).transpose() * weight * n).value() / f;
0350   // clamp to [0, 1], convert to point
0351   return std::min(std::max(u, static_cast<ActsScalar>(0)),
0352                   static_cast<ActsScalar>(1)) *
0353              n +
0354          a;
0355 }
0356
0357 double Acts::AnnulusBounds::squaredNorm(const Vector2& v,
0358                                         const SquareMatrix2& weight) const {
0359   return (v.transpose() * weight * v).value();
0360 }
0361
0362 Acts::Vector2 Acts::AnnulusBounds::moduleOrigin() const {
0363   return Eigen::Rotation2D<ActsScalar>(get(eAveragePhi)) * m_moduleOrigin;
0364 }
0365
0366 // Ostream operator overload
0367 std::ostream& Acts::AnnulusBounds::toStream(std::ostream& sl) const {
0368   sl << std::setiosflags(std::ios::fixed);
0369   sl << std::setprecision(7);
0370   sl << "Acts::AnnulusBounds:  (innerRadius, outerRadius, minPhi, maxPhi) = ";
0371   sl << "(" << get(eMinR) << ", " << get(eMaxR) << ", " << phiMin() << ", "
0372      << phiMax() << ")" << '\n';
0373   sl << " - shift xy = " << m_shiftXY.x() << ", " << m_shiftXY.y() << '\n';
0374   sl << " - shift pc = " << m_shiftPC.x() << ", " << m_shiftPC.y() << '\n';
0375   sl << std::setprecision(-1);
0376   return sl;
0377 }