Acts/Seeding/CombinatorialSeedSolver.hpp

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/Definitions/Algebra.hpp"
0012 #include "Acts/EventData/CompositeSpacePoint.hpp"
0013 #include "Acts/Surfaces/detail/PlanarHelper.hpp"
0014 #include "Acts/Utilities/Intersection.hpp"
0015 #include "Acts/Utilities/MathHelpers.hpp"
0016
0017 #include <concepts>
0018 #include <iostream>
0019 #include <ranges>
0020 #include <vector>
0021
0022 namespace Acts::Experimental::CombinatorialSeedSolver {
0023
0024 /// @brief A Combinatorial Seed Solver for seed estimation from combinatoric hits from four layers (e.g Muon NSW seeding)
0025
0026 ///***The layers equations for the four layers (Si,Di) can be */
0027 /// S1 + lambda*D1 = M
0028 /// S2 + alpha*D2 = M + A*Dm
0029 /// S3 + gamma*D3 = M + G*Dm
0030 /// S4 + kappa*D4 = M + K*Dm
0031 /// where {Si,Di}  are the position and direction of the strip
0032 /// {M,Dm} the position and direction of the muon's trajectory on the 1st plane
0033 /// solving the system we can reduce it to a 2x2 system for kappa and gamma --->
0034 /// https://gitlab.cern.ch/atlas-nextgen/work-package-2.5/analyticalsegment
0035 /// Bx1*kappa + Bx2*gamma = Y2x
0036 /// By1*kappa + By2*gamma = Y2y
0037
0038 /// @brief defines the betaMatrix calculated from the combinatoric hits
0039 /// @tparam spacePointContainer the space point container
0040 /// @param layerQuartett the space points of the combinatorics
0041 /// @return the 2x2 beta matrix
0042 template <Experimental::CompositeSpacePointPtr Point_t>
0043 SquareMatrix2 betaMatrix(const std::array<Point_t, 4>& layerQuartett) {
0044   SquareMatrix2 bMatrix{SquareMatrix2::Identity()};
0045
0046   Vector3 b1Aterm = (layerQuartett[3]->sensorDirection().dot(
0047                         layerQuartett[1]->sensorDirection())) *
0048                         layerQuartett[1]->sensorDirection() -
0049                     layerQuartett[3]->sensorDirection();
0050   Vector3 b1Gterm = (layerQuartett[3]->sensorDirection().dot(
0051                         layerQuartett[0]->sensorDirection())) *
0052                         layerQuartett[0]->sensorDirection() -
0053                     (layerQuartett[3]->sensorDirection().dot(
0054                         layerQuartett[1]->sensorDirection())) *
0055                         layerQuartett[1]->sensorDirection();
0056
0057   Vector3 b2Aterm = layerQuartett[2]->sensorDirection() -
0058                     (layerQuartett[2]->sensorDirection().dot(
0059                         layerQuartett[1]->sensorDirection())) *
0060                         layerQuartett[1]->sensorDirection();
0061   Vector3 b2Kterm = (layerQuartett[2]->sensorDirection().dot(
0062                         layerQuartett[1]->sensorDirection())) *
0063                         layerQuartett[1]->sensorDirection() -
0064                     (layerQuartett[2]->sensorDirection().dot(
0065                         layerQuartett[0]->sensorDirection())) *
0066                         layerQuartett[0]->sensorDirection();
0067
0068   // get the distances of the layers along z direction
0069   double A = (layerQuartett[0]->localPosition().z() -
0070               layerQuartett[1]->localPosition().z());
0071   double G = (layerQuartett[0]->localPosition().z() -
0072               layerQuartett[2]->localPosition().z());
0073   double K = (layerQuartett[0]->localPosition().z() -
0074               layerQuartett[3]->localPosition().z());
0075
0076   // define B matrix
0077   Vector3 b1 = A * b1Aterm + G * b1Gterm;
0078   Vector3 b2 = K * b2Kterm + A * b2Aterm;
0079
0080   bMatrix.col(0) = Vector2(b1.x(), b1.y());
0081   bMatrix.col(1) = Vector2(b2.x(), b2.y());
0082
0083   return bMatrix;
0084 }
0085
0086 /// @brief calculates the parameters lambda,alpha,gamma,kappa of the system
0087 /// @tparam spacePointContainer the space point container
0088 /// @param betaMatrix the betaMatrix for the system
0089 /// @param layerQuartett the space points of the combinatorics
0090 /// @return an array of the calculated parameters
0091 template <Experimental::CompositeSpacePointPtr Point_t>
0092 std::array<double, 4> defineParameters(
0093     const SquareMatrix2& betaMatrix,
0094     const std::array<Point_t, 4>& layerQuartett) {
0095   double A = (layerQuartett[0]->localPosition().z() -
0096               layerQuartett[1]->localPosition().z());
0097   double G = (layerQuartett[0]->localPosition().z() -
0098               layerQuartett[2]->localPosition().z());
0099   double K = (layerQuartett[0]->localPosition().z() -
0100               layerQuartett[3]->localPosition().z());
0101
0102   // Define y2 for the linear system
0103   Vector3 y0 = K * (layerQuartett[2]->localPosition() -
0104                     layerQuartett[0]->localPosition()) +
0105                G * (layerQuartett[0]->localPosition() -
0106                     layerQuartett[3]->localPosition());
0107   Vector3 y1 = A * (layerQuartett[3]->localPosition() -
0108                     layerQuartett[2]->localPosition()) +
0109                G * (layerQuartett[1]->localPosition() -
0110                     layerQuartett[3]->localPosition()) +
0111                K * (layerQuartett[2]->localPosition() -
0112                     layerQuartett[1]->localPosition());
0113   Vector3 y2 = (K - G) * (layerQuartett[0]->localPosition() -
0114                           layerQuartett[1]->localPosition()) -
0115                (y1.dot(layerQuartett[1]->sensorDirection())) *
0116                    layerQuartett[1]->sensorDirection() +
0117                (y0.dot(layerQuartett[0]->sensorDirection())) *
0118                    layerQuartett[0]->sensorDirection() +
0119                A * (layerQuartett[3]->localPosition() -
0120                     layerQuartett[2]->localPosition());
0121
0122   Vector2 solution = betaMatrix.inverse() * y2.block<2, 1>(0, 0);
0123   double kappa = solution.x();
0124   double gamma = solution.y();
0125
0126   double lambda = (y0.dot(layerQuartett[0]->sensorDirection()) +
0127                    K * gamma *
0128                        (layerQuartett[0]->sensorDirection().dot(
0129                            layerQuartett[2]->sensorDirection())) -
0130                    G * kappa *
0131                        (layerQuartett[0]->sensorDirection().dot(
0132                            layerQuartett[3]->sensorDirection()))) /
0133                   (K - G);
0134   double alpha = (y1.dot(layerQuartett[1]->sensorDirection()) +
0135                   (A - G) * kappa *
0136                       layerQuartett[3]->sensorDirection().dot(
0137                           layerQuartett[1]->sensorDirection()) +
0138                   (K - A) * gamma *
0139                       layerQuartett[2]->sensorDirection().dot(
0140                           layerQuartett[1]->sensorDirection())) /
0141                  (K - G);
0142
0143   return std::array<double, 4>({lambda, alpha, gamma, kappa});
0144 }
0145
0146 // Solve the equations for the seed position and direction (M,DM)
0147 /// @brief solves the equation system to calculate the seed
0148 /// @tparam spacePointContainr the space point container
0149 /// @param layerQuartett the space points of the combinatorics
0150 /// @param parameters the lambda,alpha,gamma,kappa parameters of the four layers
0151 /// @return the pair of the seed position and direction
0152 template <Experimental::CompositeSpacePointPtr Point_t>
0153 std::pair<Vector3, Vector3> seedSolution(
0154     const std::array<Point_t, 4>& layerQuartett,
0155     const std::array<double, 4>& parameters) {
0156   // estimate the seed positionInChamber from the 1st equation of the system of
0157   // the layer equations
0158   Vector3 seedPosition = layerQuartett[0]->localPosition() +
0159                          parameters[0] * layerQuartett[0]->sensorDirection();
0160
0161   // estimate the seed direction from the 2nd equation of the system of the
0162   // layer equations
0163   Vector3 seedDirection =
0164       ((layerQuartett[1]->localPosition() +
0165         parameters[1] * layerQuartett[1]->sensorDirection() - seedPosition))
0166           .normalized();
0167
0168   // calculate the position at z=0
0169   Intersection3D intersectionZ0 = PlanarHelper::intersectPlane(
0170       seedPosition, seedDirection, Vector3::UnitZ(), 0.0);
0171   Vector3 seedPositionZ0 = intersectionZ0.position();
0172
0173   return std::make_pair(seedPositionZ0,
0174                         copySign(seedDirection, seedDirection.z()));
0175 };
0176
0177 }  // namespace Acts::Experimental::CombinatorialSeedSolver