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