Core/Seeding/EstimateTrackParamsFromSeedTest.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/TrackParametrization.hpp"
0013 #include "Acts/Definitions/Units.hpp"
0014 #include "Acts/EventData/TrackParameters.hpp"
0015 #include "Acts/EventData/detail/TestSourceLink.hpp"
0016 #include "Acts/Geometry/GeometryContext.hpp"
0017 #include "Acts/Geometry/GeometryIdentifier.hpp"
0018 #include "Acts/Geometry/LayerCreator.hpp"
0019 #include "Acts/Geometry/TrackingGeometry.hpp"
0020 #include "Acts/MagneticField/ConstantBField.hpp"
0021 #include "Acts/MagneticField/MagneticFieldContext.hpp"
0022 #include "Acts/Propagator/EigenStepper.hpp"
0023 #include "Acts/Propagator/Navigator.hpp"
0024 #include "Acts/Propagator/Propagator.hpp"
0025 #include "Acts/Seeding/EstimateTrackParamsFromSeed.hpp"
0026 #include "Acts/Surfaces/Surface.hpp"
0027 #include "Acts/Tests/CommonHelpers/CylindricalTrackingGeometry.hpp"
0028 #include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
0029 #include "Acts/Tests/CommonHelpers/MeasurementsCreator.hpp"
0030 #include "Acts/Utilities/Logger.hpp"
0031
0032 #include <algorithm>
0033 #include <array>
0034 #include <cmath>
0035 #include <iterator>
0036 #include <map>
0037 #include <memory>
0038 #include <optional>
0039 #include <random>
0040 #include <utility>
0041 #include <vector>
0042
0043 #include "SpacePoint.hpp"
0044
0045 namespace {
0046
0047 using namespace Acts;
0048 using namespace Acts::Test;
0049 using namespace Acts::UnitLiterals;
0050
0051 using ConstantFieldStepper = Acts::EigenStepper<>;
0052 using ConstantFieldPropagator =
0053     Acts::Propagator<ConstantFieldStepper, Acts::Navigator>;
0054
0055 const GeometryContext geoCtx;
0056 const MagneticFieldContext magCtx;
0057
0058 // detector geometry
0059 CylindricalTrackingGeometry geometryStore(geoCtx);
0060 const auto geometry = geometryStore();
0061
0062 // Two dimensional measurement with zero resolution
0063 const MeasurementResolutionMap resolutions = {
0064     {GeometryIdentifier(),
0065      MeasurementResolution{MeasurementType::eLoc01, {0, 0}}}};
0066
0067 // Construct initial track parameters.
0068 CurvilinearTrackParameters makeParameters(double phi, double theta, double p,
0069                                           double q) {
0070   // create covariance matrix from reasonable standard deviations
0071   Acts::BoundVector stddev;
0072   stddev[Acts::eBoundLoc0] = 100_um;
0073   stddev[Acts::eBoundLoc1] = 100_um;
0074   stddev[Acts::eBoundTime] = 25_ns;
0075   stddev[Acts::eBoundPhi] = 2_degree;
0076   stddev[Acts::eBoundTheta] = 2_degree;
0077   stddev[Acts::eBoundQOverP] = 1 / 100_GeV;
0078   BoundSquareMatrix cov = stddev.cwiseProduct(stddev).asDiagonal();
0079   // Let the particle starts from the origin
0080   Vector4 mPos4(0., 0., 0., 0.);
0081   return CurvilinearTrackParameters(mPos4, phi, theta, q / p, cov,
0082                                     ParticleHypothesis::pionLike(std::abs(q)));
0083 }
0084
0085 std::default_random_engine rng(42);
0086
0087 }  // namespace
0088
0089 BOOST_AUTO_TEST_CASE(trackparameters_estimation_test) {
0090   // Construct a propagator with the cylinderal geometry and a constant magnetic
0091   // field along z
0092   Acts::Navigator navigator({
0093       geometry,
0094       true,  // sensitive
0095       true,  // material
0096       false  // passive
0097   });
0098   const Vector3 bField(0, 0, 2._T);
0099   auto field = std::make_shared<Acts::ConstantBField>(bField);
0100   ConstantFieldStepper stepper(std::move(field));
0101
0102   ConstantFieldPropagator propagator(std::move(stepper), std::move(navigator));
0103
0104   std::array<double, 2> pArray = {0.5_GeV, 1.0_GeV};
0105   std::array<double, 3> phiArray = {20._degree, 0._degree - 20._degree};
0106   std::array<double, 3> thetaArray = {80._degree, 90.0_degree, 100._degree};
0107   std::array<double, 2> qArray = {1, -1};
0108
0109   auto logger = Acts::getDefaultLogger("estimateTrackParamsFromSeed",
0110                                        Acts::Logging::INFO);
0111
0112   for (const auto& p : pArray) {
0113     for (const auto& phi : phiArray) {
0114       for (const auto& theta : thetaArray) {
0115         for (const auto& q : qArray) {
0116           BOOST_TEST_INFO("Test track with p = " << p << ", phi = " << phi
0117                                                  << ", theta = " << theta
0118                                                  << ", q = " << q);
0119           auto start = makeParameters(phi, theta, p, q);
0120           auto measurements = createMeasurements(propagator, geoCtx, magCtx,
0121                                                  start, resolutions, rng);
0122
0123           // Create space points from different detector layers
0124           std::map<GeometryIdentifier::Value, SpacePoint> spacePoints;
0125           const Surface* bottomSurface = nullptr;
0126           for (const auto& sl : measurements.sourceLinks) {
0127             const auto geoId = sl.m_geometryId;
0128             const auto& layer = geoId.layer();
0129             auto it = spacePoints.find(layer);
0130             // Avoid to use space point from the same layers
0131             if (it != spacePoints.end()) {
0132               continue;
0133             }
0134             const auto surface = geometry->findSurface(geoId);
0135             const auto& localPos = sl.parameters;
0136             Vector3 globalFakeMom(1, 1, 1);
0137             Vector3 globalPos =
0138                 surface->localToGlobal(geoCtx, localPos, globalFakeMom);
0139             // Create a space point (varianceR and varianceZ are lazily set to
0140             // zero since they are not important for the test)
0141             float r = std::hypot(globalPos.x(), globalPos.y());
0142             spacePoints.emplace(
0143                 layer, SpacePoint{static_cast<float>(globalPos.x()),
0144                                   static_cast<float>(globalPos.y()),
0145                                   static_cast<float>(globalPos.z()), r,
0146                                   static_cast<int>(geoId.layer()), 0., 0.,
0147                                   std::nullopt, std::nullopt});
0148             if (spacePoints.size() == 1) {
0149               bottomSurface = surface;
0150             }
0151           }
0152
0153           // Check if there are at least 3 space points
0154           if (spacePoints.size() < 3) {
0155             BOOST_TEST_WARN("Number of space points less than 3.");
0156             continue;
0157           }
0158
0159           // The truth track parameters at the bottom space point
0160           const auto& expParams = measurements.truthParameters[0];
0161           BOOST_TEST_INFO(
0162               "The truth track parameters at the bottom space point: \n"
0163               << expParams.transpose());
0164
0165           // The space point pointers
0166           std::array<const SpacePoint*, 3> spacePointPtrs{};
0167           std::transform(spacePoints.begin(), std::next(spacePoints.begin(), 3),
0168                          spacePointPtrs.begin(),
0169                          [](const auto& sp) { return &sp.second; });
0170
0171           // Test the free track parameters estimator
0172           FreeVector estFreeParams =
0173               estimateTrackParamsFromSeed(spacePointPtrs, bField);
0174           BOOST_CHECK(!estFreeParams.hasNaN());
0175
0176           // Test the bound track parameters estimator
0177           auto estFullParamsResult = estimateTrackParamsFromSeed(
0178               geoCtx, spacePointPtrs, *bottomSurface, bField);
0179           BOOST_CHECK(estFullParamsResult.ok());
0180           const auto& estFullParams = estFullParamsResult.value();
0181           BOOST_TEST_INFO(
0182               "The estimated full track parameters at the bottom space point: "
0183               "\n"
0184               << estFullParams.transpose());
0185
0186           CHECK_CLOSE_ABS(estFullParams[eBoundLoc0], expParams[eBoundLoc0],
0187                           1e-5);
0188           CHECK_CLOSE_ABS(estFullParams[eBoundLoc1], expParams[eBoundLoc1],
0189                           1e-5);
0190           // @todo Understand why the estimated phi has a limited precision
0191           CHECK_CLOSE_ABS(estFullParams[eBoundPhi], expParams[eBoundPhi], 1e-1);
0192           CHECK_CLOSE_ABS(estFullParams[eBoundTheta], expParams[eBoundTheta],
0193                           1e-2);
0194           CHECK_CLOSE_ABS(estFullParams[eBoundQOverP], expParams[eBoundQOverP],
0195                           1e-2);
0196           // time is not estimated so we check if it is default zero
0197           CHECK_CLOSE_ABS(estFullParams[eBoundTime], 0, 1e-6);
0198         }
0199       }
0200     }
0201   }
0202 }