Core/Seeding/EstimateTrackParamsFromSeedTest.cpp

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