Acts/Seeding/CombinatorialSeedSolver.ipp

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 #pragma once
0010
0011 #include "Acts/Seeding/CombinatorialSeedSolver.hpp"
0012
0013 namespace Acts::Experimental::detail {
0014
0015 template <Experimental::CompositeSpacePointPtr Point_t>
0016 std::array<std::size_t, 3> separateLayers(
0017     const std::array<Point_t, 3>& layerTriplet) {
0018   // make sure we have exacxtly one 2D and two 1D measurements available
0019   auto n2D = std::ranges::count_if(layerTriplet, [](const auto& sp) {
0020     unsigned int dimension = sp->measuresLoc0() + sp->measuresLoc1();
0021     return dimension == 2;
0022   });
0023
0024   auto n1D = std::ranges::count_if(layerTriplet, [](const auto& sp) {
0025     unsigned int dimension = sp->measuresLoc0() + sp->measuresLoc1();
0026     return dimension == 1;
0027   });
0028
0029   if (n2D != 1 || n1D != 2) {
0030     throw std::runtime_error(
0031         "Triplet must contain exactly one 2D and two 1D space points");
0032   }
0033
0034   std::array<std::size_t, 3> retIndices{
0035       std::numeric_limits<std::size_t>::max()};
0036
0037   for (std::size_t idx = 0; idx < layerTriplet.size(); ++idx) {
0038     const Point_t& spacePoint = layerTriplet[idx];
0039     unsigned int dimension =
0040         spacePoint->measuresLoc0() + spacePoint->measuresLoc1();
0041     if (dimension == 1) {
0042       retIndices[0] == std::numeric_limits<std::size_t>::max()
0043           ? retIndices[0] = idx
0044           : retIndices[1] = idx;
0045     } else {
0046       retIndices[2] = idx;
0047     }
0048   }
0049
0050   return retIndices;
0051 }
0052
0053 }  // namespace Acts::Experimental::detail
0054
0055 namespace Acts::Experimental::CombinatorialSeedSolver {
0056
0057 using Acts::Experimental::detail::separateLayers;
0058
0059 template <Experimental::CompositeSpacePointPtr Point_t>
0060 SquareMatrix2 betaMatrix(const std::array<Point_t, 4>& layerQuartett) {
0061   SquareMatrix2 bMatrix{SquareMatrix2::Identity()};
0062
0063   Vector3 b1Aterm = (layerQuartett[3]->sensorDirection().dot(
0064                         layerQuartett[1]->sensorDirection())) *
0065                         layerQuartett[1]->sensorDirection() -
0066                     layerQuartett[3]->sensorDirection();
0067   Vector3 b1Gterm = (layerQuartett[3]->sensorDirection().dot(
0068                         layerQuartett[0]->sensorDirection())) *
0069                         layerQuartett[0]->sensorDirection() -
0070                     (layerQuartett[3]->sensorDirection().dot(
0071                         layerQuartett[1]->sensorDirection())) *
0072                         layerQuartett[1]->sensorDirection();
0073
0074   Vector3 b2Aterm = layerQuartett[2]->sensorDirection() -
0075                     (layerQuartett[2]->sensorDirection().dot(
0076                         layerQuartett[1]->sensorDirection())) *
0077                         layerQuartett[1]->sensorDirection();
0078   Vector3 b2Kterm = (layerQuartett[2]->sensorDirection().dot(
0079                         layerQuartett[1]->sensorDirection())) *
0080                         layerQuartett[1]->sensorDirection() -
0081                     (layerQuartett[2]->sensorDirection().dot(
0082                         layerQuartett[0]->sensorDirection())) *
0083                         layerQuartett[0]->sensorDirection();
0084
0085   // get the distances of the layers along z direction
0086   double A = (layerQuartett[0]->localPosition().z() -
0087               layerQuartett[1]->localPosition().z());
0088   double G = (layerQuartett[0]->localPosition().z() -
0089               layerQuartett[2]->localPosition().z());
0090   double K = (layerQuartett[0]->localPosition().z() -
0091               layerQuartett[3]->localPosition().z());
0092
0093   // define B matrix
0094   Vector3 b1 = A * b1Aterm + G * b1Gterm;
0095   Vector3 b2 = K * b2Kterm + A * b2Aterm;
0096
0097   bMatrix.col(0) = Vector2(b1.x(), b1.y());
0098   bMatrix.col(1) = Vector2(b2.x(), b2.y());
0099
0100   return bMatrix;
0101 }
0102
0103 template <Experimental::CompositeSpacePointPtr Point_t>
0104 std::array<double, 4> defineParameters(
0105     const SquareMatrix2& betaMatrix,
0106     const std::array<Point_t, 4>& layerQuartett) {
0107   double A = (layerQuartett[0]->localPosition().z() -
0108               layerQuartett[1]->localPosition().z());
0109   double G = (layerQuartett[0]->localPosition().z() -
0110               layerQuartett[2]->localPosition().z());
0111   double K = (layerQuartett[0]->localPosition().z() -
0112               layerQuartett[3]->localPosition().z());
0113
0114   // Define y2 for the linear system
0115   Vector3 y0 = K * (layerQuartett[2]->localPosition() -
0116                     layerQuartett[0]->localPosition()) +
0117                G * (layerQuartett[0]->localPosition() -
0118                     layerQuartett[3]->localPosition());
0119   Vector3 y1 = A * (layerQuartett[3]->localPosition() -
0120                     layerQuartett[2]->localPosition()) +
0121                G * (layerQuartett[1]->localPosition() -
0122                     layerQuartett[3]->localPosition()) +
0123                K * (layerQuartett[2]->localPosition() -
0124                     layerQuartett[1]->localPosition());
0125   Vector3 y2 = (K - G) * (layerQuartett[0]->localPosition() -
0126                           layerQuartett[1]->localPosition()) -
0127                (y1.dot(layerQuartett[1]->sensorDirection())) *
0128                    layerQuartett[1]->sensorDirection() +
0129                (y0.dot(layerQuartett[0]->sensorDirection())) *
0130                    layerQuartett[0]->sensorDirection() +
0131                A * (layerQuartett[3]->localPosition() -
0132                     layerQuartett[2]->localPosition());
0133
0134   Vector2 solution = betaMatrix.inverse() * y2.block<2, 1>(0, 0);
0135   double kappa = solution.x();
0136   double gamma = solution.y();
0137
0138   double lambda = (y0.dot(layerQuartett[0]->sensorDirection()) +
0139                    K * gamma *
0140                        (layerQuartett[0]->sensorDirection().dot(
0141                            layerQuartett[2]->sensorDirection())) -
0142                    G * kappa *
0143                        (layerQuartett[0]->sensorDirection().dot(
0144                            layerQuartett[3]->sensorDirection()))) /
0145                   (K - G);
0146   double alpha = (y1.dot(layerQuartett[1]->sensorDirection()) +
0147                   (A - G) * kappa *
0148                       layerQuartett[3]->sensorDirection().dot(
0149                           layerQuartett[1]->sensorDirection()) +
0150                   (K - A) * gamma *
0151                       layerQuartett[2]->sensorDirection().dot(
0152                           layerQuartett[1]->sensorDirection())) /
0153                  (K - G);
0154
0155   return std::array<double, 4>({lambda, alpha, gamma, kappa});
0156 }
0157
0158 template <Experimental::CompositeSpacePointPtr Point_t>
0159 std::pair<Vector3, Vector3> seedSolution(
0160     const std::array<Point_t, 4>& layerQuartett,
0161     const std::array<double, 4>& parameters) {
0162   // estimate the seed positionInChamber from the 1st equation of the system of
0163   // the layer equations
0164   Vector3 seedPosition = layerQuartett[0]->localPosition() +
0165                          parameters[0] * layerQuartett[0]->sensorDirection();
0166
0167   // estimate the seed direction from the 2nd equation of the system of the
0168   // layer equations
0169   Vector3 seedDirection =
0170       ((layerQuartett[1]->localPosition() +
0171         parameters[1] * layerQuartett[1]->sensorDirection() - seedPosition))
0172           .normalized();
0173
0174   // calculate the position at z=0
0175   Intersection3D intersectionZ0 = PlanarHelper::intersectPlane(
0176       seedPosition, seedDirection, Vector3::UnitZ(), 0.0);
0177   Vector3 seedPositionZ0 = intersectionZ0.position();
0178
0179   return std::make_pair(seedPositionZ0,
0180                         copySign(seedDirection, seedDirection.z()));
0181 }
0182
0183 template <Experimental::CompositeSpacePointPtr Point_t>
0184 SquareMatrix2 betaMatrix(const std::array<Point_t, 3>& layerTriplet) {
0185   SquareMatrix2 bMatrix{SquareMatrix2::Identity()};
0186
0187   // keep the 1D spacepoints and the 2D spacepoint seperately- for invariance
0188   // under the order of the layers
0189   std::array<std::size_t, 3> indices = separateLayers(layerTriplet);
0190   Point_t spacePoint2D{layerTriplet[indices[2]]};
0191   std::array<Point_t, 2> spacePoints1D{layerTriplet[indices[0]],
0192                                        layerTriplet[indices[1]]};
0193
0194   double R =
0195       spacePoint2D->localPosition().z() - spacePoints1D[0]->localPosition().z();
0196   double L =
0197       spacePoint2D->localPosition().z() - spacePoints1D[1]->localPosition().z();
0198   double dotProd = spacePoints1D[0]->sensorDirection().dot(
0199       spacePoints1D[1]->sensorDirection());
0200
0201   // construct the betaMatrix
0202   bMatrix.col(0) = Vector2(L, L * dotProd);
0203   bMatrix.col(1) = Vector2(-R * dotProd, -R);
0204
0205   return bMatrix;
0206 }
0207
0208 template <Experimental::CompositeSpacePointPtr Point_t>
0209 std::array<double, 2> defineParameters(
0210     const SquareMatrix2& betaMatrix,
0211     const std::array<Point_t, 3>& layerTriplet) {
0212   std::array<std::size_t, 3> indices = separateLayers(layerTriplet);
0213   const Point_t& spacePoint2D{layerTriplet[indices[2]]};
0214   std::array<Point_t, 2> spacePoints1D{layerTriplet[indices[0]],
0215                                        layerTriplet[indices[1]]};
0216
0217   const Vector3& M = spacePoint2D->localPosition();
0218   double R =
0219       spacePoint2D->localPosition().z() - spacePoints1D[0]->localPosition().z();
0220   double L =
0221       spacePoint2D->localPosition().z() - spacePoints1D[1]->localPosition().z();
0222
0223   Vector3 y = R * spacePoints1D[1]->localPosition() -
0224               L * spacePoints1D[0]->localPosition() + (L - R) * M;
0225   Vector2 y2 = {y.dot(spacePoints1D[0]->sensorDirection()),
0226                 y.dot(spacePoints1D[1]->sensorDirection())};
0227
0228   Vector2 solution = betaMatrix.inverse() * y2.block<2, 1>(0, 0);
0229   double beta = solution.x();
0230   double delta = solution.y();
0231
0232   return std::array<double, 2>({beta, delta});
0233 }
0234
0235 template <Experimental::CompositeSpacePointPtr Point_t>
0236 std::pair<Vector3, Vector3> seedSolution(
0237     const std::array<Point_t, 3>& layerTriplet,
0238     const std::array<double, 2>& parameters) {
0239   // separate 2D and 1D spacepoints
0240   std::array<std::size_t, 3> indices = separateLayers(layerTriplet);
0241   Point_t spacePoint2D{layerTriplet[indices[2]]};
0242   std::array<Point_t, 2> spacePoints1D{layerTriplet[indices[0]],
0243                                        layerTriplet[indices[1]]};
0244
0245   // the position of the seed can be evaluated from the 2D measurement
0246   const Vector3& seedPosition = spacePoint2D->localPosition();
0247   // The direction of the seed can be estimated from the second layer equation
0248   Vector3 seedDirection =
0249       (spacePoints1D[0]->localPosition() +
0250        parameters[0] * spacePoints1D[0]->sensorDirection() - seedPosition)
0251           .normalized();
0252   // calculate the position at z=0
0253   Intersection3D intersectionZ0 = PlanarHelper::intersectPlane(
0254       seedPosition, seedDirection, Vector3::UnitZ(), 0.0);
0255   Vector3 seedPositionZ0 = intersectionZ0.position();
0256
0257   return std::make_pair(seedPositionZ0,
0258                         copySign(seedDirection, seedDirection.z()));
0259 }
0260
0261 }  // namespace Acts::Experimental::CombinatorialSeedSolver