Core/SpacePointFormation/PixelSpacePointBuilderTests.cpp

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 #include <boost/test/unit_test.hpp>
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Definitions/Units.hpp"
0013 #include "Acts/Geometry/GeometryContext.hpp"
0014 #include "Acts/SpacePointFormation/PixelSpacePointBuilder.hpp"
0015 #include "Acts/Surfaces/PlaneSurface.hpp"
0016 #include "Acts/Surfaces/RectangleBounds.hpp"
0017 #include "Acts/Utilities/MathHelpers.hpp"
0018
0019 #include <cmath>
0020 #include <memory>
0021 #include <numbers>
0022
0023 using namespace Acts;
0024 using namespace Acts::UnitLiterals;
0025
0026 namespace ActsTests {
0027
0028 namespace {
0029
0030 const GeometryContext gctx = GeometryContext::dangerouslyDefaultConstruct();
0031
0032 constexpr double sigmaLoc0 = 15_um;
0033 constexpr double sigmaLoc1 = 40_um;
0034
0035 /// A plane whose local axes are given explicitly, centred at @p center
0036 std::shared_ptr<PlaneSurface> makePlane(const Vector3& center,
0037                                         const Vector3& loc0,
0038                                         const Vector3& loc1) {
0039   Transform3 transform = Transform3::Identity();
0040   transform.matrix().block<3, 1>(0, 0) = loc0.normalized();
0041   transform.matrix().block<3, 1>(0, 1) = loc1.normalized();
0042   transform.matrix().block<3, 1>(0, 2) = loc0.cross(loc1).normalized();
0043   transform.matrix().block<3, 1>(0, 3) = center;
0044   return Surface::makeShared<PlaneSurface>(
0045       transform, std::make_shared<RectangleBounds>(50_mm, 50_mm));
0046 }
0047
0048 SquareMatrix2 localCov() {
0049   SquareMatrix2 cov = SquareMatrix2::Zero();
0050   cov(0, 0) = sigmaLoc0 * sigmaLoc0;
0051   cov(1, 1) = sigmaLoc1 * sigmaLoc1;
0052   return cov;
0053 }
0054
0055 /// The same variance from a numerically differentiated Jacobian
0056 Vector2 numericVarianceZR(const PlaneSurface& surface, const Vector3& position,
0057                           const SquareMatrix2& cov) {
0058   const SquareMatrix3 rot =
0059       surface.referenceFrame(gctx, position, Vector3::UnitZ());
0060   constexpr double h = 1e-6;
0061
0062   SquareMatrix2 jac = SquareMatrix2::Zero();
0063   for (int i = 0; i < 2; ++i) {
0064     const Vector3 step = h * rot.col(i);
0065     const Vector3 plus = position + step;
0066     const Vector3 minus = position - step;
0067     jac(0, i) = (plus.z() - minus.z()) / (2 * h);
0068     jac(1, i) =
0069         (fastHypot(plus.x(), plus.y()) - fastHypot(minus.x(), minus.y())) /
0070         (2 * h);
0071   }
0072   return (jac * cov * jac.transpose()).diagonal();
0073 }
0074
0075 }  // namespace
0076
0077 BOOST_AUTO_TEST_SUITE(PixelSpacePointBuilderSuite)
0078
0079 /// A barrel module measures r-phi and z, so all of the variance lands in z
0080 BOOST_AUTO_TEST_CASE(BarrelModule) {
0081   const double phi = 0.3;
0082   const double r = 120_mm;
0083   const Vector3 position(r * std::cos(phi), r * std::sin(phi), 250_mm);
0084   const Vector3 rPhi(-std::sin(phi), std::cos(phi), 0);
0085   const Vector3 z = Vector3::UnitZ();
0086
0087   auto surface = makePlane(position, rPhi, z);
0088   const Vector2 varZR = PixelSpacePointBuilder::computeVarianceZR(
0089       gctx, *surface, position, localCov());
0090
0091   BOOST_CHECK_CLOSE(varZR[0], sigmaLoc1 * sigmaLoc1, 1e-6);
0092   BOOST_CHECK_SMALL(varZR[1], 1e-12);
0093 }
0094
0095 /// An endcap module measures r and r-phi, so all of the variance lands in r.
0096 /// This pins the scale of the r Jacobian: a factor f in dr/d{x,y} shows as f²
0097 BOOST_AUTO_TEST_CASE(EndcapModule) {
0098   const double phi = -1.1;
0099   const double r = 250_mm;
0100   const Vector3 position(r * std::cos(phi), r * std::sin(phi), 1500_mm);
0101   const Vector3 radial(std::cos(phi), std::sin(phi), 0);
0102   const Vector3 rPhi(-std::sin(phi), std::cos(phi), 0);
0103
0104   auto surface = makePlane(position, radial, rPhi);
0105   const Vector2 varZR = PixelSpacePointBuilder::computeVarianceZR(
0106       gctx, *surface, position, localCov());
0107
0108   BOOST_CHECK_SMALL(varZR[0], 1e-12);
0109   BOOST_CHECK_CLOSE(varZR[1], sigmaLoc0 * sigmaLoc0, 1e-6);
0110 }
0111
0112 /// A module inclined against both, where the local directions mix into r
0113 BOOST_AUTO_TEST_CASE(InclinedModuleAgainstNumericJacobian) {
0114   const double phi = 0.7;
0115   const double r = 180_mm;
0116   const Vector3 position(r * std::cos(phi), r * std::sin(phi), -800_mm);
0117   const Vector3 radial(std::cos(phi), std::sin(phi), 0);
0118   const Vector3 rPhi(-std::sin(phi), std::cos(phi), 0);
0119
0120   // tilt out of the transverse plane, then rotate within the module plane
0121   const double tilt = 30 * std::numbers::pi / 180;
0122   const double skew = 20 * std::numbers::pi / 180;
0123   const Vector3 tilted =
0124       std::cos(tilt) * radial + std::sin(tilt) * Vector3::UnitZ();
0125   const Vector3 loc0 = std::cos(skew) * tilted + std::sin(skew) * rPhi;
0126   const Vector3 loc1 = -std::sin(skew) * tilted + std::cos(skew) * rPhi;
0127
0128   auto surface = makePlane(position, loc0, loc1);
0129   const SquareMatrix2 cov = localCov();
0130   const Vector2 varZR =
0131       PixelSpacePointBuilder::computeVarianceZR(gctx, *surface, position, cov);
0132   const Vector2 expected = numericVarianceZR(*surface, position, cov);
0133
0134   BOOST_CHECK_CLOSE(varZR[0], expected[0], 1e-3);
0135   BOOST_CHECK_CLOSE(varZR[1], expected[1], 1e-3);
0136 }
0137
0138 BOOST_AUTO_TEST_SUITE_END()
0139
0140 }  // namespace ActsTests