Acts/Seeding/EstimateTrackParamsFromSeed.hpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2021 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 #pragma once
0010
0011 #include "Acts/Definitions/TrackParametrization.hpp"
0012 #include "Acts/Definitions/Units.hpp"
0013 #include "Acts/Geometry/GeometryContext.hpp"
0014 #include "Acts/Seeding/Seed.hpp"
0015 #include "Acts/Surfaces/Surface.hpp"
0016 #include "Acts/Utilities/Logger.hpp"
0017
0018 #include <array>
0019 #include <cmath>
0020 #include <iostream>
0021 #include <iterator>
0022 #include <optional>
0023 #include <vector>
0024
0025 namespace Acts {
0026 /// @todo:
0027 /// 1) Implement the simple Line and Circle fit based on Taubin Circle fit
0028 /// 2) Implement the simple Line and Parabola fit (from HPS reconstruction by
0029 /// Robert Johnson)
0030
0031 /// Estimate the track parameters on the xy plane from at least three space
0032 /// points. It assumes the trajectory projection on the xy plane is a circle,
0033 /// i.e. the magnetic field is along global z-axis.
0034 ///
0035 /// The method is based on V. Karimaki NIM A305 (1991) 187-191:
0036 /// https://doi.org/10.1016/0168-9002(91)90533-V
0037 /// - no weights are used in Karimaki's fit; d0 is the distance of the point of
0038 /// closest approach to the origin, 1/R is the curvature, phi is the angle of
0039 /// the direction propagation (counter clockwise as positive) at the point of
0040 /// cloest approach.
0041 ///
0042 /// @tparam spacepoint_iterator_t The type of space point iterator
0043 ///
0044 /// @param spBegin is the begin iterator for the space points
0045 /// @param spEnd is the end iterator for the space points
0046 /// @param logger A logger instance
0047 ///
0048 /// @return optional bound track parameters with the estimated d0, phi and 1/R
0049 /// stored with the indices, eBoundLoc0, eBoundPhi and eBoundQOverP,
0050 /// respectively. The bound parameters with other indices are set to zero.
0051 template <typename spacepoint_iterator_t>
0052 std::optional<BoundVector> estimateTrackParamsFromSeed(
0053     spacepoint_iterator_t spBegin, spacepoint_iterator_t spEnd,
0054     const Logger& logger = getDummyLogger()) {
0055   // Check the number of provided space points
0056   std::size_t numSP = std::distance(spBegin, spEnd);
0057   if (numSP < 3) {
0058     ACTS_ERROR("At least three space points are required.")
0059     return std::nullopt;
0060   }
0061
0062   ActsScalar x2m = 0., xm = 0.;
0063   ActsScalar xym = 0.;
0064   ActsScalar y2m = 0., ym = 0.;
0065   ActsScalar r2m = 0., r4m = 0.;
0066   ActsScalar xr2m = 0., yr2m = 0.;
0067
0068   for (spacepoint_iterator_t it = spBegin; it != spEnd; it++) {
0069     if (*it == nullptr) {
0070       ACTS_ERROR("Empty space point found. This should not happen.")
0071       return std::nullopt;
0072     }
0073     const auto& sp = *it;
0074
0075     ActsScalar x = sp->x();
0076     ActsScalar y = sp->y();
0077     ActsScalar r2 = x * x + y * y;
0078     x2m += x * x;
0079     xm += x;
0080     xym += x * y;
0081     y2m += y * y;
0082     ym += y;
0083     r2m += r2;
0084     r4m += r2 * r2;
0085     xr2m += x * r2;
0086     yr2m += y * r2;
0087     numSP++;
0088   }
0089   x2m = x2m / numSP;
0090   xm = xm / numSP;
0091   xym = xym / numSP;
0092   y2m = y2m / numSP;
0093   ym = ym / numSP;
0094   r2m = r2m / numSP;
0095   r4m = r4m / numSP;
0096   xr2m = xr2m / numSP;
0097   yr2m = yr2m / numSP;
0098
0099   ActsScalar Cxx = x2m - xm * xm;
0100   ActsScalar Cxy = xym - xm * ym;
0101   ActsScalar Cyy = y2m - ym * ym;
0102   ActsScalar Cxr2 = xr2m - xm * r2m;
0103   ActsScalar Cyr2 = yr2m - ym * r2m;
0104   ActsScalar Cr2r2 = r4m - r2m * r2m;
0105
0106   ActsScalar q1 = Cr2r2 * Cxy - Cxr2 * Cyr2;
0107   ActsScalar q2 = Cr2r2 * (Cxx - Cyy) - Cxr2 * Cxr2 + Cyr2 * Cyr2;
0108
0109   ActsScalar phi = 0.5 * std::atan(2 * q1 / q2);
0110   ActsScalar k = (std::sin(phi) * Cxr2 - std::cos(phi) * Cyr2) * (1. / Cr2r2);
0111   ActsScalar delta = -k * r2m + std::sin(phi) * xm - std::cos(phi) * ym;
0112
0113   ActsScalar rho = (2 * k) / (std::sqrt(1 - 4 * delta * k));
0114   ActsScalar d0 = (2 * delta) / (1 + std::sqrt(1 - 4 * delta * k));
0115
0116   // Initialize the bound parameters vector
0117   BoundVector params = BoundVector::Zero();
0118   params[eBoundLoc0] = d0;
0119   params[eBoundPhi] = phi;
0120   params[eBoundQOverP] = rho;
0121
0122   return params;
0123 }
0124
0125 /// Estimate the full track parameters from three space points
0126 ///
0127 /// This method is based on the conformal map transformation. It estimates the
0128 /// full bound track parameters, i.e. (loc0, loc1, phi, theta, q/p, t) at the
0129 /// bottom space point. The bottom space is assumed to be the first element
0130 /// in the range defined by the iterators. It must lie on the surface
0131 /// provided for the representation of the bound track parameters. The magnetic
0132 /// field (which might be along any direction) is also necessary for the
0133 /// momentum estimation.
0134 ///
0135 /// It resembles the method used in ATLAS for the track parameters
0136 /// estimated from seed, i.e. the function InDet::SiTrackMaker_xk::getAtaPlane
0137 /// here:
0138 /// https://acode-browser.usatlas.bnl.gov/lxr/source/athena/InnerDetector/InDetRecTools/SiTrackMakerTool_xk/src/SiTrackMaker_xk.cxx
0139 ///
0140 /// @tparam spacepoint_iterator_t  The type of space point iterator
0141 ///
0142 /// @param gctx is the geometry context
0143 /// @param spBegin is the begin iterator for the space points
0144 /// @param spEnd is the end iterator for the space points
0145 /// @param surface is the surface of the bottom space point. The estimated bound
0146 /// track parameters will be represented also at this surface
0147 /// @param bField is the magnetic field vector
0148 /// @param bFieldMin is the minimum magnetic field required to trigger the
0149 /// estimation of q/pt
0150 /// @param logger A logger instance
0151 ///
0152 /// @return optional bound parameters
0153 template <typename spacepoint_iterator_t>
0154 std::optional<BoundVector> estimateTrackParamsFromSeed(
0155     const GeometryContext& gctx, spacepoint_iterator_t spBegin,
0156     spacepoint_iterator_t spEnd, const Surface& surface, const Vector3& bField,
0157     ActsScalar bFieldMin, const Acts::Logger& logger = getDummyLogger()) {
0158   // Check the number of provided space points
0159   std::size_t numSP = std::distance(spBegin, spEnd);
0160   if (numSP != 3) {
0161     ACTS_ERROR("There should be exactly three space points provided.")
0162     return std::nullopt;
0163   }
0164
0165   // Convert bField to Tesla
0166   ActsScalar bFieldInTesla = bField.norm() / UnitConstants::T;
0167   ActsScalar bFieldMinInTesla = bFieldMin / UnitConstants::T;
0168   // Check if magnetic field is too small
0169   if (bFieldInTesla < bFieldMinInTesla) {
0170     // @todo shall we use straight-line estimation and use default q/pt in such
0171     // case?
0172     ACTS_WARNING("The magnetic field at the bottom space point: B = "
0173                  << bFieldInTesla << " T is smaller than |B|_min = "
0174                  << bFieldMinInTesla << " T. Estimation is not performed.")
0175     return std::nullopt;
0176   }
0177
0178   // The global positions of the bottom, middle and space points
0179   std::array<Vector3, 3> spGlobalPositions = {Vector3::Zero(), Vector3::Zero(),
0180                                               Vector3::Zero()};
0181   std::array<std::optional<float>, 3> spGlobalTimes = {
0182       std::nullopt, std::nullopt, std::nullopt};
0183   // The first, second and third space point are assumed to be bottom, middle
0184   // and top space point, respectively
0185   for (std::size_t isp = 0; isp < 3; ++isp) {
0186     spacepoint_iterator_t it = std::next(spBegin, isp);
0187     if (*it == nullptr) {
0188       ACTS_ERROR("Empty space point found. This should not happen.")
0189       return std::nullopt;
0190     }
0191     const auto& sp = *it;
0192     spGlobalPositions[isp] = Vector3(sp->x(), sp->y(), sp->z());
0193     spGlobalTimes[isp] = sp->t();
0194   }
0195
0196   // Define a new coordinate frame with its origin at the bottom space point, z
0197   // axis long the magnetic field direction and y axis perpendicular to vector
0198   // from the bottom to middle space point. Hence, the projection of the middle
0199   // space point on the transverse plane will be located at the x axis of the
0200   // new frame.
0201   Vector3 relVec = spGlobalPositions[1] - spGlobalPositions[0];
0202   Vector3 newZAxis = bField.normalized();
0203   Vector3 newYAxis = newZAxis.cross(relVec).normalized();
0204   Vector3 newXAxis = newYAxis.cross(newZAxis);
0205   RotationMatrix3 rotation;
0206   rotation.col(0) = newXAxis;
0207   rotation.col(1) = newYAxis;
0208   rotation.col(2) = newZAxis;
0209   // The center of the new frame is at the bottom space point
0210   Translation3 trans(spGlobalPositions[0]);
0211   // The transform which constructs the new frame
0212   Transform3 transform(trans * rotation);
0213
0214   // The coordinate of the middle and top space point in the new frame
0215   Vector3 local1 = transform.inverse() * spGlobalPositions[1];
0216   Vector3 local2 = transform.inverse() * spGlobalPositions[2];
0217
0218   // In the new frame the bottom sp is at the origin, while the middle
0219   // sp in along the x axis. As such, the x-coordinate of the circle is
0220   // at: x-middle / 2.
0221   // The y coordinate can be found by using the straight line passing
0222   // between the mid point between the middle and top sp and perpendicular to
0223   // the line connecting them
0224   Vector2 circleCenter;
0225   circleCenter(0) = 0.5 * local1(0);
0226
0227   ActsScalar deltaX21 = local2(0) - local1(0);
0228   ActsScalar sumX21 = local2(0) + local1(0);
0229   // straight line connecting the two points
0230   // y = a * x + c (we don't care about c right now)
0231   // we simply need the slope
0232   // we compute 1./a since this is what we need for the following computation
0233   ActsScalar ia = deltaX21 / local2(1);
0234   // Perpendicular line is then y = -1/a *x + b
0235   // we can evaluate b given we know a already by imposing
0236   // the line passes through P = (0.5 * (x2 + x1), 0.5 * y2)
0237   ActsScalar b = 0.5 * (local2(1) + ia * sumX21);
0238   circleCenter(1) = -ia * circleCenter(0) + b;
0239   // Radius is a signed distance between circleCenter and first sp, which is at
0240   // (0, 0) in the new frame. Sign depends on the slope a (positive vs negative)
0241   int sign = ia > 0 ? -1 : 1;
0242   const ActsScalar R = circleCenter.norm();
0243   ActsScalar invTanTheta =
0244       local2.z() /
0245       (2.f * R * std::asin(std::hypot(local2.x(), local2.y()) / (2.f * R)));
0246   // The momentum direction in the new frame (the center of the circle has the
0247   // coordinate (-1.*A/(2*B), 1./(2*B)))
0248   ActsScalar A = -circleCenter(0) / circleCenter(1);
0249   Vector3 transDirection(1., A, std::hypot(1, A) * invTanTheta);
0250   // Transform it back to the original frame
0251   Vector3 direction = rotation * transDirection.normalized();
0252
0253   // Initialize the bound parameters vector
0254   BoundVector params = BoundVector::Zero();
0255
0256   // The estimated phi and theta
0257   params[eBoundPhi] = VectorHelpers::phi(direction);
0258   params[eBoundTheta] = VectorHelpers::theta(direction);
0259
0260   // Transform the bottom space point to local coordinates of the provided
0261   // surface
0262   auto lpResult = surface.globalToLocal(gctx, spGlobalPositions[0], direction);
0263   if (!lpResult.ok()) {
0264     ACTS_ERROR(
0265         "Global to local transformation did not succeed. Please make sure the "
0266         "bottom space point lies on the provided surface.");
0267     return std::nullopt;
0268   }
0269   Vector2 bottomLocalPos = lpResult.value();
0270   // The estimated loc0 and loc1
0271   params[eBoundLoc0] = bottomLocalPos.x();
0272   params[eBoundLoc1] = bottomLocalPos.y();
0273   params[eBoundTime] = spGlobalTimes[0].value_or(0.);
0274
0275   // The estimated q/pt in [GeV/c]^-1 (note that the pt is the projection of
0276   // momentum on the transverse plane of the new frame)
0277   ActsScalar qOverPt = sign * (UnitConstants::m) / (0.3 * bFieldInTesla * R);
0278   // The estimated q/p in [GeV/c]^-1
0279   params[eBoundQOverP] = qOverPt / std::hypot(1., invTanTheta);
0280
0281   if (params.hasNaN()) {
0282     ACTS_ERROR(
0283         "The NaN value exists at the estimated track parameters from seed with"
0284         << "\nbottom sp: " << spGlobalPositions[0] << "\nmiddle sp: "
0285         << spGlobalPositions[1] << "\ntop sp: " << spGlobalPositions[2]);
0286     return std::nullopt;
0287   }
0288   return params;
0289 }
0290
0291 }  // namespace Acts