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