Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-08-05 08:09:48 

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/Units.hpp"
0012 #include "Acts/Seeding/detail/StrawLineFitAuxiliaries.hpp"
0013 #include "Acts/Surfaces/detail/LineHelper.hpp"
0014 #include "Acts/Surfaces/detail/PlanarHelper.hpp"
0015 #include "Acts/Utilities/StringHelpers.hpp"
0016 
0017 using namespace Acts;
0018 using namespace Acts::detail;
0019 using namespace Acts::Experimental::detail;
0020 using namespace Acts::UnitLiterals;
0021 using namespace Acts::PlanarHelper;
0022 
0023 using Line_t = StrawLineFitAuxiliaries::Line_t;
0024 using Config_t = StrawLineFitAuxiliaries::Config;
0025 using ParIdx = StrawLineFitAuxiliaries::FitParIndex;
0026 
0027 using Vector = Line_t::Vector;
0028 using Pars_t = Line_t::ParamVector;
0029 
0030 namespace Acts::Test {
0031 class TestSpacePoint {
0032  public:
0033   /// @brief Mask to express which track directions are measured
0034   ///        by the space point
0035   enum ProjectorMask : std::size_t {
0036     bendingDir = 1 << 0,
0037     nonBendingDir = 1 << 1,
0038     timeOfArrival = 1 << 2,
0039     bothDirections = bendingDir | nonBendingDir
0040   };
0041 
0042   TestSpacePoint(
0043       const Vector3& pos, const Vector3& wireDir, const double radius,
0044       bool measNonPrec = false,
0045       const std::array<double, 3>& cov = Acts::filledArray<double, 3>(0.))
0046       : m_pos{pos},
0047         m_dir{wireDir},
0048         m_radius{radius},
0049         m_dirMask{measNonPrec ? bothDirections : bendingDir},
0050         m_cov{cov} {}
0051 
0052   TestSpacePoint(
0053       const Vector3& pos, const Vector3& stripDir, const Vector3& stripNorm,
0054       ProjectorMask mask,
0055       const std::array<double, 3>& cov = Acts::filledArray<double, 3>(0.))
0056       : m_pos{pos},
0057         m_dir{stripDir},
0058         m_toNext{stripNorm},
0059         m_dirMask{mask},
0060         m_cov{cov},
0061         m_isStraw{false} {}
0062 
0063   const Vector3& localPosition() const { return m_pos; }
0064   const Vector3& sensorDirection() const { return m_dir; }
0065   const Vector3& toNextSensor() const { return m_toNext; }
0066   const Vector3& planeNormal() const { return m_planeNorm; }
0067 
0068   bool isStraw() const { return m_isStraw; }
0069   bool hasTime() const {
0070     return (m_dirMask & (ProjectorMask::timeOfArrival)) != 0u;
0071   }
0072   bool measuresLoc1() const {
0073     return (m_dirMask & (ProjectorMask::bendingDir)) != 0u;
0074   }
0075   bool measuresLoc0() const {
0076     return (m_dirMask & (ProjectorMask::nonBendingDir)) != 0u;
0077   }
0078 
0079   double driftRadius() const { return m_radius; }
0080   double time() const { return m_time; }
0081   const std::array<double, 3>& covariance() const { return m_cov; }
0082 
0083  private:
0084   Vector3 m_pos{Vector3::Zero()};
0085   Vector3 m_dir{Vector3::Zero()};
0086   Vector3 m_toNext{Vector3::Zero()};
0087   Vector3 m_planeNorm{m_dir.cross(m_toNext).normalized()};
0088   double m_radius{0.};
0089   double m_time{0.};
0090   ProjectorMask m_dirMask{ProjectorMask::bothDirections};
0091   std::array<double, 3> m_cov{Acts::filledArray<double, 3>(0.)};
0092   bool m_isStraw{true};
0093 };
0094 
0095 BOOST_AUTO_TEST_SUITE(StrawLineSeederTest)
0096 
0097 void testResidual(const Pars_t& linePars, const TestSpacePoint& testPoint) {
0098   constexpr auto logLvl = Logging::Level::INFO;
0099   using namespace Acts::detail::LineHelper;
0100   using ResidualIdx = StrawLineFitAuxiliaries::ResidualIdx;
0101   Config_t resCfg{};
0102   resCfg.useHessian = true;
0103   resCfg.calcAlongStrip = false;
0104   resCfg.parsToUse = {ParIdx::x0, ParIdx::y0, ParIdx::phi, ParIdx::theta};
0105   Line_t line{};
0106   line.updateParameters(linePars);
0107 
0108   std::cout << "\n\n\nResidual test - Test line: " << toString(line.position())
0109             << ", " << toString(line.direction()) << std::endl;
0110 
0111   StrawLineFitAuxiliaries resCalc{resCfg,
0112                                   Acts::getDefaultLogger("testRes", logLvl)};
0113   resCalc.updateSpatialResidual(line, testPoint);
0114   if (testPoint.isStraw()) {
0115     const double lineDist =
0116         signedDistance(testPoint.localPosition(), testPoint.sensorDirection(),
0117                        line.position(), line.direction());
0118     std::cout << "Test residual w.r.t. sp with position: "
0119               << toString(testPoint.localPosition())
0120               << ", orientation: " << toString(testPoint.sensorDirection())
0121               << ", r: " << testPoint.driftRadius() << std::endl;
0122 
0123     std::cout << "Residual: " << toString(resCalc.residual())
0124               << ", line distance: " << lineDist << std::endl;
0125     ///
0126     BOOST_CHECK_CLOSE(resCalc.residual()[ResidualIdx::bending],
0127                       lineDist - testPoint.driftRadius(), 1.e-12);
0128   } else {
0129     std::cout << "Test residual w.r.t strip space point -- position: "
0130               << toString(testPoint.localPosition())
0131               << ", normal: " << toString(testPoint.planeNormal())
0132               << ", strip: " << toString(testPoint.sensorDirection())
0133               << ", to-next: " << toString(testPoint.toNextSensor())
0134               << std::endl;
0135     const Vector3& n{testPoint.planeNormal()};
0136     const Vector3 planeIsect = intersectPlane(line.position(), line.direction(),
0137                                               n, testPoint.localPosition())
0138                                    .position();
0139     const Vector3 delta = (planeIsect - testPoint.localPosition());
0140 
0141     std::cout << "Residual: " << toString(resCalc.residual())
0142               << ", delta: " << toString(delta)
0143               << ", <delta, n> =" << delta.dot(n) << std::endl;
0144 
0145     Acts::ActsSquareMatrix<3> coordTrf{Acts::ActsSquareMatrix<3>::Zero()};
0146     coordTrf.col(ResidualIdx::bending) = testPoint.measuresLoc1()
0147                                              ? testPoint.toNextSensor()
0148                                              : testPoint.sensorDirection();
0149     coordTrf.col(ResidualIdx::nonBending) = testPoint.measuresLoc1()
0150                                                 ? testPoint.sensorDirection()
0151                                                 : testPoint.toNextSensor();
0152     coordTrf.col(2) = n;
0153     std::cout << "Coordinate trf: \n" << coordTrf << std::endl;
0154 
0155     const Vector3 trfDelta = coordTrf.inverse() * delta;
0156     std::cout << "Transformed delta: " << toString(trfDelta) << std::endl;
0157 
0158     if (testPoint.measuresLoc1()) {
0159       BOOST_CHECK_CLOSE(trfDelta[ResidualIdx::bending],
0160                         resCalc.residual()[ResidualIdx::bending], 1.e-10);
0161     } else {
0162       BOOST_CHECK_CLOSE(trfDelta[ResidualIdx::bending] *
0163                             static_cast<double>(resCfg.calcAlongStrip),
0164                         resCalc.residual()[ResidualIdx::bending], 1.e-10);
0165     }
0166 
0167     if (testPoint.measuresLoc0()) {
0168       BOOST_CHECK_CLOSE(trfDelta[ResidualIdx::nonBending],
0169                         resCalc.residual()[ResidualIdx::nonBending], 1.e-10);
0170     } else {
0171       BOOST_CHECK_CLOSE(trfDelta[ResidualIdx::nonBending] *
0172                             static_cast<double>(resCfg.calcAlongStrip),
0173                         resCalc.residual()[ResidualIdx::nonBending], 1.e-10);
0174     }
0175   }
0176 
0177   constexpr double h = 5.e-9;
0178   constexpr double tolerance = 1.e-3;
0179   for (auto par : resCfg.parsToUse) {
0180     Pars_t lineParsUp{linePars}, lineParsDn{linePars};
0181     lineParsUp[static_cast<std::size_t>(par)] += h;
0182     lineParsDn[static_cast<std::size_t>(par)] -= h;
0183     Line_t lineUp{}, lineDn{};
0184     lineUp.updateParameters(lineParsUp);
0185     lineDn.updateParameters(lineParsDn);
0186 
0187     StrawLineFitAuxiliaries resCalcUp{
0188         resCfg, Acts::getDefaultLogger("testResUp", logLvl)};
0189     StrawLineFitAuxiliaries resCalcDn{
0190         resCfg, Acts::getDefaultLogger("testResDn", logLvl)};
0191 
0192     resCalcUp.updateSpatialResidual(lineUp, testPoint);
0193     resCalcDn.updateSpatialResidual(lineDn, testPoint);
0194 
0195     const Vector numDeriv =
0196         (resCalcUp.residual() - resCalcDn.residual()) / (2. * h);
0197 
0198     std::cout << "Derivative test: " << StrawLineFitAuxiliaries::parName(par)
0199               << ", derivative: " << toString(resCalc.gradient(par))
0200               << ",  numerical: " << toString(numDeriv) << std::endl;
0201     BOOST_CHECK_LE((numDeriv - resCalc.gradient(par)).norm() /
0202                        std::max(numDeriv.norm(), 1.),
0203                    tolerance);
0204     /// Next attempt to calculate the second derivative
0205     for (auto par1 : resCfg.parsToUse) {
0206       lineParsUp = linePars;
0207       lineParsDn = linePars;
0208       lineParsUp[static_cast<std::size_t>(par1)] += h;
0209       lineParsDn[static_cast<std::size_t>(par1)] -= h;
0210 
0211       lineUp.updateParameters(lineParsUp);
0212       lineDn.updateParameters(lineParsDn);
0213 
0214       resCalcUp.updateSpatialResidual(lineUp, testPoint);
0215       resCalcDn.updateSpatialResidual(lineDn, testPoint);
0216 
0217       const Vector numDeriv1{
0218           (resCalcUp.gradient(par) - resCalcDn.gradient(par)) / (2. * h)};
0219       const Vector& analyticDeriv = resCalc.hessian(par, par1);
0220       std::cout << "Second deriv (" << StrawLineFitAuxiliaries::parName(par)
0221                 << ", " << StrawLineFitAuxiliaries::parName(par1)
0222                 << ") -- numerical: " << toString(numDeriv1)
0223                 << ", analytic: " << toString(analyticDeriv) << std::endl;
0224       BOOST_CHECK_LE((numDeriv1 - analyticDeriv).norm(), tolerance);
0225     }
0226   }
0227 }
0228 
0229 BOOST_AUTO_TEST_CASE(WireResidualTest) {
0230   /// Set the line to be 45 degrees
0231   using Pars_t = Line_t::ParamVector;
0232   using ParIdx = Line_t::ParIndex;
0233   Pars_t linePars{};
0234   linePars[static_cast<std::size_t>(ParIdx::phi)] = 90. * 1_degree;
0235   linePars[static_cast<std::size_t>(ParIdx::theta)] = 45. * 1_degree;
0236 
0237   const Vector wireDir = Vector::UnitX();
0238   // /// Generate the first test measurement
0239   testResidual(linePars,
0240                TestSpacePoint{Vector{100. * 1_cm, 50. * 1_cm, 30. * 1_cm},
0241                               wireDir, 10. * 1_cm});
0242   testResidual(
0243       linePars,
0244       TestSpacePoint{Vector{100. * 1_cm, 50. * 1_cm, 30. * 1_cm},
0245                      makeDirectionFromPhiTheta(10. * 1_degree, 30. * 1_degree),
0246                      10. * 1_cm});
0247 
0248   testResidual(linePars,
0249                TestSpacePoint{Vector{100. * 1_cm, 50. * 1_cm, 30. * 1_cm},
0250                               wireDir, 10. * 1_cm, true});
0251 
0252   linePars[static_cast<std::size_t>(ParIdx::phi)] = 30. * 1_degree;
0253   linePars[static_cast<std::size_t>(ParIdx::theta)] = 60. * 1_degree;
0254 
0255   testResidual(linePars,
0256                TestSpacePoint{Vector{100. * 1_cm, 50. * 1_cm, 30. * 1_cm},
0257                               wireDir, 10. * 1_cm});
0258 
0259   testResidual(linePars,
0260                TestSpacePoint{Vector{100. * 1_cm, 50. * 1_cm, 30. * 1_cm},
0261                               wireDir, 10. * 1_cm, true});
0262 
0263   linePars[static_cast<std::size_t>(ParIdx::phi)] = 60. * 1_degree;
0264   linePars[static_cast<std::size_t>(ParIdx::theta)] = 30. * 1_degree;
0265 
0266   testResidual(linePars,
0267                TestSpacePoint{Vector{100. * 1_cm, 50. * 1_cm, 30. * 1_cm},
0268                               wireDir, 10. * 1_cm});
0269   testResidual(linePars,
0270                TestSpacePoint{Vector{100. * 1_cm, 50. * 1_cm, 30. * 1_cm},
0271                               wireDir, 10. * 1_cm, true});
0272 
0273   testResidual(
0274       linePars,
0275       TestSpacePoint{Vector{100. * 1_cm, 50. * 1_cm, 30. * 1_cm},
0276                      makeDirectionFromPhiTheta(20. * 1_degree, 30 * 1_degree),
0277                      10. * 1_cm});
0278   testResidual(
0279       linePars,
0280       TestSpacePoint{Vector{100. * 1_cm, 50. * 1_cm, 30. * 1_cm},
0281                      makeDirectionFromPhiTheta(40. * 1_degree, 50 * 1_degree),
0282                      10. * 1_cm, true});
0283 }
0284 
0285 BOOST_AUTO_TEST_CASE(StripResidual) {
0286   Pars_t linePars{};
0287   linePars[static_cast<std::size_t>(ParIdx::phi)] = 60. * 1_degree;
0288   linePars[static_cast<std::size_t>(ParIdx::theta)] = 45 * 1_degree;
0289 
0290   testResidual(
0291       linePars,
0292       TestSpacePoint{Vector{75. * 1_cm, -75. * 1_cm, 100. * 1_cm},
0293                      makeDirectionFromPhiTheta(90. * 1_degree, 90. * 1_degree),
0294                      makeDirectionFromPhiTheta(0. * 1_degree, 90 * 1_degree),
0295                      TestSpacePoint::bendingDir});
0296 
0297   testResidual(
0298       linePars,
0299       TestSpacePoint{Vector{75. * 1_cm, -75. * 1_cm, 100. * 1_cm},
0300                      makeDirectionFromPhiTheta(0. * 1_degree, 90. * 1_degree),
0301                      makeDirectionFromPhiTheta(90. * 1_degree, 90 * 1_degree),
0302                      TestSpacePoint::nonBendingDir});
0303 
0304   testResidual(
0305       linePars,
0306       TestSpacePoint{Vector{75. * 1_cm, -75. * 1_cm, 100. * 1_cm},
0307                      makeDirectionFromPhiTheta(90. * 1_degree, 90. * 1_degree),
0308                      makeDirectionFromPhiTheta(0. * 1_degree, 90 * 1_degree),
0309                      TestSpacePoint::bothDirections});
0310 
0311   testResidual(
0312       linePars,
0313       TestSpacePoint{Vector{75. * 1_cm, -75. * 1_cm, 100. * 1_cm},
0314                      makeDirectionFromPhiTheta(30. * 1_degree, 90. * 1_degree),
0315                      makeDirectionFromPhiTheta(60. * 1_degree, 90 * 1_degree),
0316                      TestSpacePoint::bothDirections});
0317 }
0318 
0319 }  // namespace Acts::Test
0320 BOOST_AUTO_TEST_SUITE_END()
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!