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