Acts/Seeding/SeedFinderUtils.ipp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2023 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 namespace Acts {
0010 template <typename external_spacepoint_t>
0011 inline LinCircle transformCoordinates(
0012     const InternalSpacePoint<external_spacepoint_t>& sp,
0013     const InternalSpacePoint<external_spacepoint_t>& spM, bool bottom) {
0014   auto extractFunction =
0015       [](const InternalSpacePoint<external_spacepoint_t>& obj)
0016       -> std::array<float, 6> {
0017     std::array<float, 6> output{obj.x(),      obj.y(),         obj.z(),
0018                                 obj.radius(), obj.varianceR(), obj.varianceZ()};
0019     return output;
0020   };
0021
0022   return transformCoordinates<InternalSpacePoint<external_spacepoint_t>>(
0023       sp, spM, bottom, std::move(extractFunction));
0024 }
0025
0026 template <typename external_spacepoint_t, typename callable_t>
0027 inline LinCircle transformCoordinates(const external_spacepoint_t& sp,
0028                                       const external_spacepoint_t& spM,
0029                                       bool bottom,
0030                                       callable_t&& extractFunction) {
0031   // The computation inside this function is exactly identical to that in the
0032   // vectorized version of this function, except that it operates on a single
0033   // spacepoint. Please see the other version of this function for more
0034   // detailed comments.
0035
0036   auto [xM, yM, zM, rM, varianceRM, varianceZM] = extractFunction(spM);
0037   auto [xSP, ySP, zSP, rSP, varianceRSP, varianceZSP] = extractFunction(sp);
0038
0039   float cosPhiM = xM / rM;
0040   float sinPhiM = yM / rM;
0041   float deltaX = xSP - xM;
0042   float deltaY = ySP - yM;
0043   float deltaZ = zSP - zM;
0044   float xNewFrame = deltaX * cosPhiM + deltaY * sinPhiM;
0045   float yNewFrame = deltaY * cosPhiM - deltaX * sinPhiM;
0046   float deltaR2 = (xNewFrame * xNewFrame + yNewFrame * yNewFrame);
0047   float iDeltaR2 = 1. / (deltaX * deltaX + deltaY * deltaY);
0048   float iDeltaR = std::sqrt(iDeltaR2);
0049   int bottomFactor = bottom ? -1 : 1;
0050   float cotTheta = deltaZ * iDeltaR * bottomFactor;
0051
0052   // conformal transformation u=x/(x²+y²) v=y/(x²+y²) transform the
0053   // circle into straight lines in the u/v plane the line equation can
0054   // be described in terms of aCoef and bCoef, where v = aCoef * u +
0055   // bCoef
0056   const float U = xNewFrame * iDeltaR2;
0057   const float V = yNewFrame * iDeltaR2;
0058
0059   // error term for sp-pair without correlation of middle space point
0060   const float Er = ((varianceZM + varianceZSP) +
0061                     (cotTheta * cotTheta) * (varianceRM + varianceRSP)) *
0062                    iDeltaR2;
0063
0064   sp.setDeltaR(std::sqrt(deltaR2 + (deltaZ * deltaZ)));
0065   return LinCircle(cotTheta, iDeltaR, Er, U, V, xNewFrame, yNewFrame);
0066 }
0067
0068 template <typename external_spacepoint_t>
0069 inline void transformCoordinates(
0070     Acts::SpacePointData& spacePointData,
0071     const std::vector<InternalSpacePoint<external_spacepoint_t>*>& vec,
0072     const InternalSpacePoint<external_spacepoint_t>& spM, bool bottom,
0073     std::vector<LinCircle>& linCircleVec) {
0074   auto extractFunction =
0075       [](const InternalSpacePoint<external_spacepoint_t>& obj)
0076       -> std::array<float, 6> {
0077     std::array<float, 6> output{obj.x(),      obj.y(),         obj.z(),
0078                                 obj.radius(), obj.varianceR(), obj.varianceZ()};
0079     return output;
0080   };
0081
0082   transformCoordinates<InternalSpacePoint<external_spacepoint_t>>(
0083       spacePointData, vec, spM, bottom, linCircleVec,
0084       std::move(extractFunction));
0085 }
0086
0087 template <typename external_spacepoint_t, typename callable_t>
0088 inline void transformCoordinates(Acts::SpacePointData& spacePointData,
0089                                  const std::vector<external_spacepoint_t*>& vec,
0090                                  const external_spacepoint_t& spM, bool bottom,
0091                                  std::vector<LinCircle>& linCircleVec,
0092                                  callable_t&& extractFunction) {
0093   auto [xM, yM, zM, rM, varianceRM, varianceZM] = extractFunction(spM);
0094
0095   // resize + operator[] is faster than reserve and push_back
0096   linCircleVec.resize(vec.size());
0097
0098   float cosPhiM = xM / rM;
0099   float sinPhiM = yM / rM;
0100
0101   int bottomFactor = bottom ? -1 : 1;
0102
0103   for (std::size_t idx(0); idx < vec.size(); ++idx) {
0104     auto& sp = vec[idx];
0105     auto [xSP, ySP, zSP, rSP, varianceRSP, varianceZSP] = extractFunction(*sp);
0106
0107     float deltaX = xSP - xM;
0108     float deltaY = ySP - yM;
0109     float deltaZ = zSP - zM;
0110     // calculate projection fraction of spM->sp vector pointing in same
0111     // direction as
0112     // vector origin->spM (x) and projection fraction of spM->sp vector pointing
0113     // orthogonal to origin->spM (y)
0114     float xNewFrame = deltaX * cosPhiM + deltaY * sinPhiM;
0115     float yNewFrame = deltaY * cosPhiM - deltaX * sinPhiM;
0116     // 1/(length of M -> SP)
0117     float deltaR2 = (xNewFrame * xNewFrame + yNewFrame * yNewFrame);
0118     float iDeltaR2 = 1. / deltaR2;
0119     float iDeltaR = std::sqrt(iDeltaR2);
0120     //
0121     // cot_theta = (deltaZ/deltaR)
0122     float cotTheta = deltaZ * iDeltaR * bottomFactor;
0123     // transformation of circle equation (x,y) into linear equation (u,v)
0124     // x^2 + y^2 - 2x_0*x - 2y_0*y = 0
0125     // is transformed into
0126     // 1 - 2x_0*u - 2y_0*v = 0
0127     // using the following m_U and m_V
0128     // (u = A + B*v); A and B are created later on
0129     float U = xNewFrame * iDeltaR2;
0130     float V = yNewFrame * iDeltaR2;
0131     // error term for sp-pair without correlation of middle space point
0132     float Er = ((varianceZM + varianceZSP) +
0133                 (cotTheta * cotTheta) * (varianceRM + varianceRSP)) *
0134                iDeltaR2;
0135
0136     // Fill Line Circle
0137     linCircleVec[idx].cotTheta = cotTheta;
0138     linCircleVec[idx].iDeltaR = iDeltaR;
0139     linCircleVec[idx].Er = Er;
0140     linCircleVec[idx].U = U;
0141     linCircleVec[idx].V = V;
0142     linCircleVec[idx].x = xNewFrame;
0143     linCircleVec[idx].y = yNewFrame;
0144
0145     spacePointData.setDeltaR(sp->index(),
0146                              std::sqrt(deltaR2 + (deltaZ * deltaZ)));
0147   }
0148 }
0149
0150 template <typename external_spacepoint_t>
0151 inline bool xyzCoordinateCheck(
0152     Acts::SpacePointData& spacePointData,
0153     const Acts::SeedFinderConfig<external_spacepoint_t>& m_config,
0154     const Acts::InternalSpacePoint<external_spacepoint_t>& sp,
0155     const double* spacepointPosition, double* outputCoordinates) {
0156   // check the compatibility of SPs coordinates in xyz assuming the
0157   // Bottom-Middle direction with the strip measurement details
0158   bool hasValueStored = spacePointData.hasDynamicVariable();
0159   if (!hasValueStored) {
0160     return false;
0161   }
0162
0163   std::size_t index = sp.index();
0164
0165   // prepare variables
0166   const Acts::Vector3& topStripVector = spacePointData.getTopStripVector(index);
0167   const Acts::Vector3& bottomStripVector =
0168       spacePointData.getBottomStripVector(index);
0169   const Acts::Vector3& stripCenterDistance =
0170       spacePointData.getStripCenterDistance(index);
0171
0172   const double xTopStripVector = topStripVector[0];
0173   const double yTopStripVector = topStripVector[1];
0174   const double zTopStripVector = topStripVector[2];
0175   const double xBottomStripVector = bottomStripVector[0];
0176   const double yBottomStripVector = bottomStripVector[1];
0177   const double zBottomStripVector = bottomStripVector[2];
0178
0179   // cross product between top strip vector and spacepointPosition
0180   double d1[3] = {yTopStripVector * spacepointPosition[2] -
0181                       zTopStripVector * spacepointPosition[1],
0182                   zTopStripVector * spacepointPosition[0] -
0183                       xTopStripVector * spacepointPosition[2],
0184                   xTopStripVector * spacepointPosition[1] -
0185                       yTopStripVector * spacepointPosition[0]};
0186
0187   // scalar product between bottom strip vector and d1
0188   double bd1 = xBottomStripVector * d1[0] + yBottomStripVector * d1[1] +
0189                zBottomStripVector * d1[2];
0190
0191   // compatibility check using distance between strips to evaluate if
0192   // spacepointPosition is inside the bottom detector element
0193   double s1 = (stripCenterDistance[0] * d1[0] + stripCenterDistance[1] * d1[1] +
0194                stripCenterDistance[2] * d1[2]);
0195   if (std::abs(s1) > std::abs(bd1) * m_config.toleranceParam) {
0196     return false;
0197   }
0198
0199   // cross product between bottom strip vector and spacepointPosition
0200   double d0[3] = {yBottomStripVector * spacepointPosition[2] -
0201                       zBottomStripVector * spacepointPosition[1],
0202                   zBottomStripVector * spacepointPosition[0] -
0203                       xBottomStripVector * spacepointPosition[2],
0204                   xBottomStripVector * spacepointPosition[1] -
0205                       yBottomStripVector * spacepointPosition[0]};
0206
0207   // compatibility check using distance between strips to evaluate if
0208   // spacepointPosition is inside the top detector element
0209   double s0 = (stripCenterDistance[0] * d0[0] + stripCenterDistance[1] * d0[1] +
0210                stripCenterDistance[2] * d0[2]);
0211   if (std::abs(s0) > std::abs(bd1) * m_config.toleranceParam) {
0212     return false;
0213   }
0214
0215   // if arrive here spacepointPosition is compatible with strip directions and
0216   // detector elements
0217
0218   const Acts::Vector3& topStripCenterPosition =
0219       spacePointData.getTopStripCenterPosition(index);
0220
0221   // spacepointPosition corrected with respect to the top strip position and
0222   // direction and the distance between the strips
0223   s0 = s0 / bd1;
0224   outputCoordinates[0] = topStripCenterPosition[0] + xTopStripVector * s0;
0225   outputCoordinates[1] = topStripCenterPosition[1] + yTopStripVector * s0;
0226   outputCoordinates[2] = topStripCenterPosition[2] + zTopStripVector * s0;
0227   return true;
0228 }
0229 }  // namespace Acts