Core/TrackFitting/Gx2fTests.cpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2023 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/Units.hpp"
0013 #include "Acts/EventData/VectorMultiTrajectory.hpp"
0014 #include "Acts/EventData/VectorTrackContainer.hpp"
0015 #include "Acts/EventData/detail/TestSourceLink.hpp"
0016 #include "Acts/Geometry/CuboidVolumeBuilder.hpp"
0017 #include "Acts/Geometry/Layer.hpp"
0018 #include "Acts/Geometry/TrackingGeometry.hpp"
0019 #include "Acts/Geometry/TrackingGeometryBuilder.hpp"
0020 #include "Acts/MagneticField/ConstantBField.hpp"
0021 #include "Acts/Material/HomogeneousSurfaceMaterial.hpp"
0022 #include "Acts/Material/HomogeneousVolumeMaterial.hpp"
0023 #include "Acts/Material/MaterialSlab.hpp"
0024 #include "Acts/Propagator/EigenStepper.hpp"
0025 #include "Acts/Propagator/Navigator.hpp"
0026 #include "Acts/Propagator/Propagator.hpp"
0027 #include "Acts/Propagator/StraightLineStepper.hpp"
0028 #include "Acts/Surfaces/RectangleBounds.hpp"
0029 #include "Acts/Tests/CommonHelpers/DetectorElementStub.hpp"
0030 #include "Acts/Tests/CommonHelpers/MeasurementsCreator.hpp"
0031 #include "Acts/Tests/CommonHelpers/PredefinedMaterials.hpp"
0032 #include "Acts/TrackFitting/GlobalChiSquareFitter.hpp"
0033 #include "Acts/Utilities/Logger.hpp"
0034
0035 #include <vector>
0036
0037 #include "FitterTestsCommon.hpp"
0038
0039 using namespace Acts::UnitLiterals;
0040 using namespace Acts::detail::Test;
0041
0042 Acts::Logging::Level logLevel = Acts::Logging::VERBOSE;
0043 const auto gx2fLogger = Acts::getDefaultLogger("Gx2f", logLevel);
0044
0045 namespace Acts::Test {
0046
0047 //// Construct initial track parameters.
0048 Acts::CurvilinearTrackParameters makeParameters(
0049     const ActsScalar x = 0.0_m, const ActsScalar y = 0.0_m,
0050     const ActsScalar z = 0.0_m, const ActsScalar w = 42_ns,
0051     const ActsScalar phi = 0_degree, const ActsScalar theta = 90_degree,
0052     const ActsScalar p = 2_GeV, const ActsScalar q = 1_e) {
0053   // create covariance matrix from reasonable standard deviations
0054   Acts::BoundVector stddev;
0055   stddev[Acts::eBoundLoc0] = 100_um;
0056   stddev[Acts::eBoundLoc1] = 100_um;
0057   stddev[Acts::eBoundTime] = 25_ns;
0058   stddev[Acts::eBoundPhi] = 2_degree;
0059   stddev[Acts::eBoundTheta] = 2_degree;
0060   stddev[Acts::eBoundQOverP] = 1 / 100_GeV;
0061   const Acts::BoundSquareMatrix cov = stddev.cwiseProduct(stddev).asDiagonal();
0062   // define a track in the transverse plane along x
0063   const Acts::Vector4 mPos4(x, y, z, w);
0064   return Acts::CurvilinearTrackParameters(mPos4, phi, theta, q / p, cov,
0065                                           Acts::ParticleHypothesis::pion());
0066 }
0067
0068 static std::vector<Acts::SourceLink> prepareSourceLinks(
0069     const std::vector<TestSourceLink>& sourceLinks) {
0070   std::vector<Acts::SourceLink> result;
0071   std::transform(sourceLinks.begin(), sourceLinks.end(),
0072                  std::back_inserter(result),
0073                  [](const auto& sl) { return Acts::SourceLink{sl}; });
0074   return result;
0075 }
0076
0077 /// @brief Create a simple telescope detector
0078 ///
0079 /// @param geoCtx
0080 /// @param nSurfaces Number of surfaces
0081 std::shared_ptr<const TrackingGeometry> makeToyDetector(
0082     const Acts::GeometryContext& geoCtx, const std::size_t nSurfaces = 5) {
0083   if (nSurfaces < 1) {
0084     throw std::invalid_argument("At least 1 surfaces needs to be created.");
0085   }
0086
0087   // Define the dimensions of the square surfaces
0088   const double halfSizeSurface = 1_m;
0089
0090   // Rotation of the surfaces
0091   const double rotationAngle = M_PI * 0.5;
0092   const Vector3 xPos(cos(rotationAngle), 0., sin(rotationAngle));
0093   const Vector3 yPos(0., 1., 0.);
0094   const Vector3 zPos(-sin(rotationAngle), 0., cos(rotationAngle));
0095
0096   // Construct builder
0097   CuboidVolumeBuilder cvb;
0098
0099   // Create configurations for surfaces
0100   std::vector<CuboidVolumeBuilder::SurfaceConfig> surfaceConfig;
0101   for (unsigned int surfPos = 1; surfPos <= nSurfaces; surfPos++) {
0102     // Position of the surfaces
0103     CuboidVolumeBuilder::SurfaceConfig cfg;
0104     cfg.position = {surfPos * UnitConstants::m, 0., 0.};
0105
0106     // Rotation of the surfaces
0107     cfg.rotation.col(0) = xPos;
0108     cfg.rotation.col(1) = yPos;
0109     cfg.rotation.col(2) = zPos;
0110
0111     // Boundaries of the surfaces (shape)
0112     cfg.rBounds = std::make_shared<const RectangleBounds>(
0113         RectangleBounds(halfSizeSurface, halfSizeSurface));
0114
0115     // Material of the surfaces
0116     MaterialSlab matProp(makeBeryllium(), 0.5_mm);
0117     cfg.surMat = std::make_shared<HomogeneousSurfaceMaterial>(matProp);
0118
0119     // Thickness of the detector element
0120     cfg.thickness = 1_um;
0121
0122     cfg.detElementConstructor =
0123         [](const Transform3& trans,
0124            const std::shared_ptr<const RectangleBounds>& bounds,
0125            double thickness) {
0126           return new DetectorElementStub(trans, bounds, thickness);
0127         };
0128     surfaceConfig.push_back(cfg);
0129   }
0130
0131   // Build layer configurations
0132   std::vector<CuboidVolumeBuilder::LayerConfig> layerConfig;
0133   for (auto& sCfg : surfaceConfig) {
0134     CuboidVolumeBuilder::LayerConfig cfg;
0135     cfg.surfaceCfg = {sCfg};
0136     cfg.active = true;
0137     cfg.envelopeX = {-0.1_mm, 0.1_mm};
0138     cfg.envelopeY = {-0.1_mm, 0.1_mm};
0139     cfg.envelopeZ = {-0.1_mm, 0.1_mm};
0140     layerConfig.push_back(cfg);
0141   }
0142
0143   // Inner Volume - Build volume configuration
0144   CuboidVolumeBuilder::VolumeConfig volumeConfig;
0145   volumeConfig.length = {(nSurfaces + 1) * 1_m, 2 * halfSizeSurface,
0146                          2 * halfSizeSurface};
0147   volumeConfig.position = {volumeConfig.length.x() / 2, 0., 0.};
0148   volumeConfig.layerCfg = layerConfig;
0149   volumeConfig.name = "Test volume";
0150
0151   // Outer volume - Build TrackingGeometry configuration
0152   CuboidVolumeBuilder::Config config;
0153   config.length = {(nSurfaces + 1) * 1_m, 2 * halfSizeSurface,
0154                    2 * halfSizeSurface};
0155   config.position = {volumeConfig.length.x() / 2, 0., 0.};
0156   config.volumeCfg = {volumeConfig};
0157
0158   cvb.setConfig(config);
0159
0160   TrackingGeometryBuilder::Config tgbCfg;
0161
0162   tgbCfg.trackingVolumeBuilders.push_back(
0163       [=](const auto& context, const auto& inner, const auto&) {
0164         return cvb.trackingVolume(context, inner, nullptr);
0165       });
0166
0167   TrackingGeometryBuilder tgb(tgbCfg);
0168
0169   std::unique_ptr<const TrackingGeometry> detector =
0170       tgb.trackingGeometry(geoCtx);
0171   return detector;
0172 }
0173
0174 struct Detector {
0175   // geometry
0176   std::shared_ptr<const TrackingGeometry> geometry;
0177 };
0178
0179 BOOST_AUTO_TEST_SUITE(Gx2fTest)
0180 ACTS_LOCAL_LOGGER(Acts::getDefaultLogger("Gx2fTests", logLevel))
0181
0182 // Context objects
0183 const Acts::GeometryContext geoCtx;
0184 const Acts::MagneticFieldContext magCtx;
0185 const Acts::CalibrationContext calCtx;
0186
0187 // Measurement resolutions
0188 const MeasurementResolution resPixel = {MeasurementType::eLoc01,
0189                                         {25_um, 50_um}};
0190 const MeasurementResolution resStrip0 = {MeasurementType::eLoc0, {25_um}};
0191 const MeasurementResolution resStrip1 = {MeasurementType::eLoc1, {50_um}};
0192 const MeasurementResolutionMap resMapAllPixel = {
0193     {Acts::GeometryIdentifier().setVolume(0), resPixel}};
0194
0195 // This test checks if the call to the fitter works and no errors occur in the
0196 // framework, without fitting and updating any parameters
0197 BOOST_AUTO_TEST_CASE(NoFit) {
0198   ACTS_INFO("*** Test: NoFit -- Start");
0199
0200   std::default_random_engine rng(42);
0201
0202   ACTS_DEBUG("Create the detector");
0203   const std::size_t nSurfaces = 5;
0204   Detector detector;
0205   detector.geometry = makeToyDetector(geoCtx, nSurfaces);
0206
0207   ACTS_DEBUG("Set the start parameters for measurement creation and fit");
0208   const auto parametersMeasurements = makeParameters();
0209   const auto startParametersFit = makeParameters(
0210       7_mm, 11_mm, 15_mm, 42_ns, 10_degree, 80_degree, 1_GeV, 1_e);
0211
0212   ACTS_DEBUG("Create the measurements");
0213   using SimPropagator =
0214       Acts::Propagator<Acts::StraightLineStepper, Acts::Navigator>;
0215   const SimPropagator simPropagator = makeStraightPropagator(detector.geometry);
0216   const auto measurements =
0217       createMeasurements(simPropagator, geoCtx, magCtx, parametersMeasurements,
0218                          resMapAllPixel, rng);
0219   const auto sourceLinks = prepareSourceLinks(measurements.sourceLinks);
0220   ACTS_VERBOSE("sourceLinks.size() = " << sourceLinks.size());
0221
0222   ACTS_DEBUG("Set up the fitter");
0223   // Reuse the SimPropagator, since we will not actually use it
0224   using Gx2Fitter =
0225       Experimental::Gx2Fitter<SimPropagator, VectorMultiTrajectory>;
0226   const Gx2Fitter fitter(simPropagator, gx2fLogger->clone());
0227
0228   const Surface* rSurface = &parametersMeasurements.referenceSurface();
0229
0230   Experimental::Gx2FitterExtensions<VectorMultiTrajectory> extensions;
0231   extensions.calibrator
0232       .connect<&testSourceLinkCalibrator<VectorMultiTrajectory>>();
0233   TestSourceLink::SurfaceAccessor surfaceAccessor{*detector.geometry};
0234   extensions.surfaceAccessor
0235       .connect<&TestSourceLink::SurfaceAccessor::operator()>(&surfaceAccessor);
0236
0237   Experimental::Gx2FitterOptions gx2fOptions(
0238       geoCtx, magCtx, calCtx, extensions, PropagatorPlainOptions(), rSurface,
0239       false, false, FreeToBoundCorrection(false), 0, true, 0);
0240
0241   Acts::TrackContainer tracks{Acts::VectorTrackContainer{},
0242                               Acts::VectorMultiTrajectory{}};
0243
0244   ACTS_DEBUG("Fit the track");
0245   ACTS_VERBOSE("startParameter unsmeared:\n" << parametersMeasurements);
0246   ACTS_VERBOSE("startParameter fit:\n" << startParametersFit);
0247   const auto res = fitter.fit(sourceLinks.begin(), sourceLinks.end(),
0248                               startParametersFit, gx2fOptions, tracks);
0249
0250   BOOST_REQUIRE(res.ok());
0251
0252   const auto& track = *res;
0253   BOOST_CHECK_EQUAL(track.tipIndex(), Acts::MultiTrajectoryTraits::kInvalid);
0254   BOOST_CHECK(track.hasReferenceSurface());
0255
0256   // Track quantities
0257   BOOST_CHECK_EQUAL(track.chi2(), 0.);
0258   BOOST_CHECK_EQUAL(track.nDoF(), 0u);
0259   BOOST_CHECK_EQUAL(track.nHoles(), 0u);
0260   BOOST_CHECK_EQUAL(track.nMeasurements(), 0u);
0261   BOOST_CHECK_EQUAL(track.nSharedHits(), 0u);
0262   BOOST_CHECK_EQUAL(track.nOutliers(), 0u);
0263
0264   // Parameters
0265   BOOST_CHECK_EQUAL(track.parameters(), startParametersFit.parameters());
0266   BOOST_CHECK_EQUAL(track.covariance(), BoundMatrix::Identity());
0267
0268   // Convergence
0269   BOOST_CHECK_EQUAL(
0270       (track.template component<
0271           std::size_t,
0272           hashString(Experimental::Gx2fConstants::gx2fnUpdateColumn)>()),
0273       0);
0274
0275   ACTS_INFO("*** Test: NoFit -- Finish");
0276 }
0277
0278 BOOST_AUTO_TEST_CASE(Fit5Iterations) {
0279   ACTS_INFO("*** Test: Fit5Iterations -- Start");
0280
0281   std::default_random_engine rng(42);
0282
0283   ACTS_DEBUG("Create the detector");
0284   const std::size_t nSurfaces = 5;
0285   Detector detector;
0286   detector.geometry = makeToyDetector(geoCtx, nSurfaces);
0287
0288   ACTS_DEBUG("Set the start parameters for measurement creation and fit");
0289   const auto parametersMeasurements = makeParameters();
0290   const auto startParametersFit = makeParameters(
0291       7_mm, 11_mm, 15_mm, 42_ns, 10_degree, 80_degree, 1_GeV, 1_e);
0292
0293   ACTS_DEBUG("Create the measurements");
0294   using SimPropagator =
0295       Acts::Propagator<Acts::StraightLineStepper, Acts::Navigator>;
0296   const SimPropagator simPropagator = makeStraightPropagator(detector.geometry);
0297   const auto measurements =
0298       createMeasurements(simPropagator, geoCtx, magCtx, parametersMeasurements,
0299                          resMapAllPixel, rng);
0300   const auto sourceLinks = prepareSourceLinks(measurements.sourceLinks);
0301   ACTS_VERBOSE("sourceLinks.size() = " << sourceLinks.size());
0302
0303   BOOST_REQUIRE_EQUAL(sourceLinks.size(), nSurfaces);
0304
0305   ACTS_DEBUG("Set up the fitter");
0306   const Surface* rSurface = &parametersMeasurements.referenceSurface();
0307
0308   using RecoStepper = EigenStepper<>;
0309   const auto recoPropagator =
0310       makeConstantFieldPropagator<RecoStepper>(detector.geometry, 0_T);
0311
0312   using RecoPropagator = decltype(recoPropagator);
0313   using Gx2Fitter =
0314       Experimental::Gx2Fitter<RecoPropagator, VectorMultiTrajectory>;
0315   const Gx2Fitter fitter(recoPropagator, gx2fLogger->clone());
0316
0317   Experimental::Gx2FitterExtensions<VectorMultiTrajectory> extensions;
0318   extensions.calibrator
0319       .connect<&testSourceLinkCalibrator<VectorMultiTrajectory>>();
0320   TestSourceLink::SurfaceAccessor surfaceAccessor{*detector.geometry};
0321   extensions.surfaceAccessor
0322       .connect<&TestSourceLink::SurfaceAccessor::operator()>(&surfaceAccessor);
0323
0324   const Experimental::Gx2FitterOptions gx2fOptions(
0325       geoCtx, magCtx, calCtx, extensions, PropagatorPlainOptions(), rSurface,
0326       false, false, FreeToBoundCorrection(false), 5, true, 0);
0327
0328   Acts::TrackContainer tracks{Acts::VectorTrackContainer{},
0329                               Acts::VectorMultiTrajectory{}};
0330
0331   ACTS_DEBUG("Fit the track");
0332   ACTS_VERBOSE("startParameter unsmeared:\n" << parametersMeasurements);
0333   ACTS_VERBOSE("startParameter fit:\n" << startParametersFit);
0334   const auto res = fitter.fit(sourceLinks.begin(), sourceLinks.end(),
0335                               startParametersFit, gx2fOptions, tracks);
0336
0337   BOOST_REQUIRE(res.ok());
0338
0339   const auto& track = *res;
0340
0341   BOOST_CHECK_EQUAL(track.tipIndex(), nSurfaces - 1);
0342   BOOST_CHECK(track.hasReferenceSurface());
0343
0344   // Track quantities
0345   CHECK_CLOSE_ABS(track.chi2(), 8., 2.);
0346   BOOST_CHECK_EQUAL(track.nDoF(), 10u);
0347   BOOST_CHECK_EQUAL(track.nHoles(), 0u);
0348   BOOST_CHECK_EQUAL(track.nMeasurements(), nSurfaces);
0349   BOOST_CHECK_EQUAL(track.nSharedHits(), 0u);
0350   BOOST_CHECK_EQUAL(track.nOutliers(), 0u);
0351
0352   // Parameters
0353   // We need quite coarse checks here, since on different builds
0354   // the created measurements differ in the randomness
0355   BOOST_CHECK_CLOSE(track.parameters()[eBoundLoc0], -11., 7e0);
0356   BOOST_CHECK_CLOSE(track.parameters()[eBoundLoc1], -15., 6e0);
0357   BOOST_CHECK_CLOSE(track.parameters()[eBoundPhi], 1e-5, 1e3);
0358   BOOST_CHECK_CLOSE(track.parameters()[eBoundTheta], M_PI / 2, 1e-3);
0359   BOOST_CHECK_EQUAL(track.parameters()[eBoundQOverP], 1);
0360   BOOST_CHECK_CLOSE(track.parameters()[eBoundTime], 12591.2832360000, 1e-6);
0361   BOOST_CHECK_CLOSE(track.covariance().determinant(), 1e-27, 4e0);
0362
0363   // Convergence
0364   BOOST_CHECK_EQUAL(
0365       (track.template component<
0366           std::size_t,
0367           hashString(Experimental::Gx2fConstants::gx2fnUpdateColumn)>()),
0368       5);
0369
0370   ACTS_INFO("*** Test: Fit5Iterations -- Finish");
0371 }
0372
0373 BOOST_AUTO_TEST_CASE(MixedDetector) {
0374   ACTS_INFO("*** Test: MixedDetector -- Start");
0375
0376   std::default_random_engine rng(42);
0377
0378   ACTS_DEBUG("Create the detector");
0379   const std::size_t nSurfaces = 7;
0380   Detector detector;
0381   detector.geometry = makeToyDetector(geoCtx, nSurfaces);
0382
0383   ACTS_DEBUG("Set the start parameters for measurement creation and fit");
0384   const auto parametersMeasurements = makeParameters();
0385   const auto startParametersFit = makeParameters(
0386       7_mm, 11_mm, 15_mm, 42_ns, 10_degree, 80_degree, 1_GeV, 1_e);
0387
0388   ACTS_DEBUG("Create the measurements");
0389   const MeasurementResolutionMap resMap = {
0390       {Acts::GeometryIdentifier().setVolume(2).setLayer(2), resPixel},
0391       {Acts::GeometryIdentifier().setVolume(2).setLayer(4), resStrip0},
0392       {Acts::GeometryIdentifier().setVolume(2).setLayer(6), resStrip1},
0393       {Acts::GeometryIdentifier().setVolume(2).setLayer(8), resPixel},
0394       {Acts::GeometryIdentifier().setVolume(2).setLayer(10), resStrip0},
0395       {Acts::GeometryIdentifier().setVolume(2).setLayer(12), resStrip1},
0396       {Acts::GeometryIdentifier().setVolume(2).setLayer(14), resPixel},
0397   };
0398
0399   using SimPropagator =
0400       Acts::Propagator<Acts::StraightLineStepper, Acts::Navigator>;
0401   const SimPropagator simPropagator = makeStraightPropagator(detector.geometry);
0402   const auto measurements = createMeasurements(
0403       simPropagator, geoCtx, magCtx, parametersMeasurements, resMap, rng);
0404   const auto sourceLinks = prepareSourceLinks(measurements.sourceLinks);
0405   ACTS_VERBOSE("sourceLinks.size() = " << sourceLinks.size());
0406
0407   BOOST_REQUIRE_EQUAL(sourceLinks.size(), nSurfaces);
0408
0409   ACTS_DEBUG("Set up the fitter");
0410   const Surface* rSurface = &parametersMeasurements.referenceSurface();
0411
0412   using RecoStepper = EigenStepper<>;
0413   const auto recoPropagator =
0414       makeConstantFieldPropagator<RecoStepper>(detector.geometry, 0_T);
0415
0416   using RecoPropagator = decltype(recoPropagator);
0417   using Gx2Fitter =
0418       Experimental::Gx2Fitter<RecoPropagator, VectorMultiTrajectory>;
0419   const Gx2Fitter fitter(recoPropagator, gx2fLogger->clone());
0420
0421   Experimental::Gx2FitterExtensions<VectorMultiTrajectory> extensions;
0422   extensions.calibrator
0423       .connect<&testSourceLinkCalibrator<VectorMultiTrajectory>>();
0424   TestSourceLink::SurfaceAccessor surfaceAccessor{*detector.geometry};
0425   extensions.surfaceAccessor
0426       .connect<&TestSourceLink::SurfaceAccessor::operator()>(&surfaceAccessor);
0427
0428   const Experimental::Gx2FitterOptions gx2fOptions(
0429       geoCtx, magCtx, calCtx, extensions, PropagatorPlainOptions(), rSurface,
0430       false, false, FreeToBoundCorrection(false), 5, true, 0);
0431
0432   Acts::TrackContainer tracks{Acts::VectorTrackContainer{},
0433                               Acts::VectorMultiTrajectory{}};
0434
0435   ACTS_DEBUG("Fit the track");
0436   ACTS_VERBOSE("startParameter unsmeared:\n" << parametersMeasurements);
0437   ACTS_VERBOSE("startParameter fit:\n" << startParametersFit);
0438   const auto res = fitter.fit(sourceLinks.begin(), sourceLinks.end(),
0439                               startParametersFit, gx2fOptions, tracks);
0440
0441   BOOST_REQUIRE(res.ok());
0442
0443   const auto& track = *res;
0444   BOOST_CHECK_EQUAL(track.tipIndex(), nSurfaces - 1);
0445   BOOST_CHECK(track.hasReferenceSurface());
0446
0447   // Track quantities
0448   CHECK_CLOSE_ABS(track.chi2(), 8.5, 4.);
0449   BOOST_CHECK_EQUAL(track.nDoF(), 10u);
0450   BOOST_CHECK_EQUAL(track.nHoles(), 0u);
0451   BOOST_CHECK_EQUAL(track.nMeasurements(), nSurfaces);
0452   BOOST_CHECK_EQUAL(track.nSharedHits(), 0u);
0453   BOOST_CHECK_EQUAL(track.nOutliers(), 0u);
0454
0455   // Parameters
0456   // We need quite coarse checks here, since on different builds
0457   // the created measurements differ in the randomness
0458   BOOST_CHECK_CLOSE(track.parameters()[eBoundLoc0], -11., 7e0);
0459   BOOST_CHECK_CLOSE(track.parameters()[eBoundLoc1], -15., 6e0);
0460   BOOST_CHECK_CLOSE(track.parameters()[eBoundPhi], 1e-5, 1e3);
0461   BOOST_CHECK_CLOSE(track.parameters()[eBoundTheta], M_PI / 2, 1e-3);
0462   BOOST_CHECK_EQUAL(track.parameters()[eBoundQOverP], 1);
0463   BOOST_CHECK_CLOSE(track.parameters()[eBoundTime], 12591.2832360000, 1e-6);
0464   BOOST_CHECK_CLOSE(track.covariance().determinant(), 2e-28, 1e0);
0465
0466   // Convergence
0467   BOOST_CHECK_EQUAL(
0468       (track.template component<
0469           std::size_t,
0470           hashString(Experimental::Gx2fConstants::gx2fnUpdateColumn)>()),
0471       5);
0472
0473   ACTS_INFO("*** Test: MixedDetector -- Finish");
0474 }
0475
0476 // This test checks if we can fit QOverP, when a magnetic field is introduced
0477 BOOST_AUTO_TEST_CASE(FitWithBfield) {
0478   ACTS_INFO("*** Test: FitWithBfield -- Start");
0479
0480   std::default_random_engine rng(42);
0481
0482   ACTS_DEBUG("Create the detector");
0483   const std::size_t nSurfaces = 5;
0484   Detector detector;
0485   detector.geometry = makeToyDetector(geoCtx, nSurfaces);
0486
0487   ACTS_DEBUG("Set the start parameters for measurement creation and fit");
0488   const auto parametersMeasurements = makeParameters();
0489   const auto startParametersFit = makeParameters(
0490       7_mm, 11_mm, 15_mm, 42_ns, 10_degree, 80_degree, 1_GeV, 1_e);
0491
0492   ACTS_DEBUG("Create the measurements");
0493   using SimStepper = EigenStepper<>;
0494   const auto simPropagator =
0495       makeConstantFieldPropagator<SimStepper>(detector.geometry, 0.3_T);
0496
0497   const auto measurements =
0498       createMeasurements(simPropagator, geoCtx, magCtx, parametersMeasurements,
0499                          resMapAllPixel, rng);
0500
0501   const auto sourceLinks = prepareSourceLinks(measurements.sourceLinks);
0502   ACTS_VERBOSE("sourceLinks.size() = " << sourceLinks.size());
0503
0504   BOOST_REQUIRE_EQUAL(sourceLinks.size(), nSurfaces);
0505
0506   ACTS_DEBUG("Set up the fitter");
0507   const Surface* rSurface = &parametersMeasurements.referenceSurface();
0508
0509   // Reuse the SimPropagator, since it already uses the EigenStepper<>
0510   using SimPropagator = decltype(simPropagator);
0511   using Gx2Fitter =
0512       Experimental::Gx2Fitter<SimPropagator, VectorMultiTrajectory>;
0513   const Gx2Fitter fitter(simPropagator, gx2fLogger->clone());
0514
0515   Experimental::Gx2FitterExtensions<VectorMultiTrajectory> extensions;
0516   extensions.calibrator
0517       .connect<&testSourceLinkCalibrator<VectorMultiTrajectory>>();
0518   TestSourceLink::SurfaceAccessor surfaceAccessor{*detector.geometry};
0519   extensions.surfaceAccessor
0520       .connect<&TestSourceLink::SurfaceAccessor::operator()>(&surfaceAccessor);
0521
0522   const Experimental::Gx2FitterOptions gx2fOptions(
0523       geoCtx, magCtx, calCtx, extensions, PropagatorPlainOptions(), rSurface,
0524       false, false, FreeToBoundCorrection(false), 5, false, 0);
0525
0526   Acts::TrackContainer tracks{Acts::VectorTrackContainer{},
0527                               Acts::VectorMultiTrajectory{}};
0528
0529   ACTS_DEBUG("Fit the track");
0530   ACTS_VERBOSE("startParameter unsmeared:\n" << parametersMeasurements);
0531   ACTS_VERBOSE("startParameter fit:\n" << startParametersFit);
0532   const auto res = fitter.fit(sourceLinks.begin(), sourceLinks.end(),
0533                               startParametersFit, gx2fOptions, tracks);
0534
0535   BOOST_REQUIRE(res.ok());
0536
0537   const auto& track = *res;
0538   BOOST_CHECK_EQUAL(track.tipIndex(), nSurfaces - 1);
0539   BOOST_CHECK(track.hasReferenceSurface());
0540
0541   // Track quantities
0542   CHECK_CLOSE_ABS(track.chi2(), 7.5, 1.5);
0543   BOOST_CHECK_EQUAL(track.nDoF(), 10u);
0544   BOOST_CHECK_EQUAL(track.nHoles(), 0u);
0545   BOOST_CHECK_EQUAL(track.nMeasurements(), nSurfaces);
0546   BOOST_CHECK_EQUAL(track.nSharedHits(), 0u);
0547   BOOST_CHECK_EQUAL(track.nOutliers(), 0u);
0548
0549   // Parameters
0550   // We need quite coarse checks here, since on different builds
0551   // the created measurements differ in the randomness
0552   // TODO investigate further the reference values for eBoundPhi and det(cov)
0553   BOOST_CHECK_CLOSE(track.parameters()[eBoundLoc0], -11., 8e0);
0554   BOOST_CHECK_CLOSE(track.parameters()[eBoundLoc1], -15., 6e0);
0555   BOOST_CHECK_CLOSE(track.parameters()[eBoundPhi], 1e-4, 1e3);
0556   BOOST_CHECK_CLOSE(track.parameters()[eBoundTheta], M_PI / 2, 1e-3);
0557   BOOST_CHECK_CLOSE(track.parameters()[eBoundQOverP], 0.5, 2e-1);
0558   BOOST_CHECK_CLOSE(track.parameters()[eBoundTime], 12591.2832360000, 1e-6);
0559   BOOST_CHECK_CLOSE(track.covariance().determinant(), 8e-35, 4e0);
0560
0561   // Convergence
0562   BOOST_CHECK_EQUAL(
0563       (track.template component<
0564           std::size_t,
0565           hashString(Experimental::Gx2fConstants::gx2fnUpdateColumn)>()),
0566       5);
0567
0568   ACTS_INFO("*** Test: FitWithBfield -- Finish");
0569 }
0570
0571 BOOST_AUTO_TEST_CASE(relChi2changeCutOff) {
0572   ACTS_INFO("*** Test: relChi2changeCutOff -- Start");
0573
0574   std::default_random_engine rng(42);
0575
0576   ACTS_DEBUG("Create the detector");
0577   const std::size_t nSurfaces = 5;
0578   Detector detector;
0579   detector.geometry = makeToyDetector(geoCtx, nSurfaces);
0580
0581   ACTS_DEBUG("Set the start parameters for measurement creation and fit");
0582   const auto parametersMeasurements = makeParameters();
0583   const auto startParametersFit = makeParameters(
0584       7_mm, 11_mm, 15_mm, 42_ns, 10_degree, 80_degree, 1_GeV, 1_e);
0585
0586   ACTS_DEBUG("Create the measurements");
0587   // simulation propagator
0588   using SimPropagator =
0589       Acts::Propagator<Acts::StraightLineStepper, Acts::Navigator>;
0590   const SimPropagator simPropagator = makeStraightPropagator(detector.geometry);
0591   const auto measurements =
0592       createMeasurements(simPropagator, geoCtx, magCtx, parametersMeasurements,
0593                          resMapAllPixel, rng);
0594   const auto sourceLinks = prepareSourceLinks(measurements.sourceLinks);
0595   ACTS_VERBOSE("sourceLinks.size() = " << sourceLinks.size());
0596
0597   BOOST_REQUIRE_EQUAL(sourceLinks.size(), nSurfaces);
0598
0599   ACTS_DEBUG("Set up the fitter");
0600   const Surface* rSurface = &parametersMeasurements.referenceSurface();
0601
0602   using RecoStepper = EigenStepper<>;
0603   const auto recoPropagator =
0604       makeConstantFieldPropagator<RecoStepper>(detector.geometry, 0_T);
0605
0606   using RecoPropagator = decltype(recoPropagator);
0607   using Gx2Fitter =
0608       Experimental::Gx2Fitter<RecoPropagator, VectorMultiTrajectory>;
0609   const Gx2Fitter fitter(recoPropagator, gx2fLogger->clone());
0610
0611   Experimental::Gx2FitterExtensions<VectorMultiTrajectory> extensions;
0612   extensions.calibrator
0613       .connect<&testSourceLinkCalibrator<VectorMultiTrajectory>>();
0614   TestSourceLink::SurfaceAccessor surfaceAccessor{*detector.geometry};
0615   extensions.surfaceAccessor
0616       .connect<&TestSourceLink::SurfaceAccessor::operator()>(&surfaceAccessor);
0617
0618   const Experimental::Gx2FitterOptions gx2fOptions(
0619       geoCtx, magCtx, calCtx, extensions, PropagatorPlainOptions(), rSurface,
0620       false, false, FreeToBoundCorrection(false), 500, true, 1e-5);
0621
0622   Acts::TrackContainer tracks{Acts::VectorTrackContainer{},
0623                               Acts::VectorMultiTrajectory{}};
0624
0625   ACTS_DEBUG("Fit the track");
0626   ACTS_VERBOSE("startParameter unsmeared:\n" << parametersMeasurements);
0627   ACTS_VERBOSE("startParameter fit:\n" << startParametersFit);
0628   const auto res = fitter.fit(sourceLinks.begin(), sourceLinks.end(),
0629                               startParametersFit, gx2fOptions, tracks);
0630
0631   BOOST_REQUIRE(res.ok());
0632
0633   const auto& track = *res;
0634   BOOST_CHECK_EQUAL(track.tipIndex(), nSurfaces - 1);
0635   BOOST_CHECK(track.hasReferenceSurface());
0636
0637   // Track quantities
0638   CHECK_CLOSE_ABS(track.chi2(), 8., 2.);
0639   BOOST_CHECK_EQUAL(track.nDoF(), 10u);
0640   BOOST_CHECK_EQUAL(track.nHoles(), 0u);
0641   BOOST_CHECK_EQUAL(track.nMeasurements(), nSurfaces);
0642   BOOST_CHECK_EQUAL(track.nSharedHits(), 0u);
0643   BOOST_CHECK_EQUAL(track.nOutliers(), 0u);
0644
0645   // Parameters
0646   // We need quite coarse checks here, since on different builds
0647   // the created measurements differ in the randomness
0648   BOOST_CHECK_CLOSE(track.parameters()[eBoundLoc0], -11., 7e0);
0649   BOOST_CHECK_CLOSE(track.parameters()[eBoundLoc1], -15., 6e0);
0650   BOOST_CHECK_CLOSE(track.parameters()[eBoundPhi], 1e-5, 1e3);
0651   BOOST_CHECK_CLOSE(track.parameters()[eBoundTheta], M_PI / 2, 1e-3);
0652   BOOST_CHECK_EQUAL(track.parameters()[eBoundQOverP], 1);
0653   BOOST_CHECK_CLOSE(track.parameters()[eBoundTime], 12591.2832360000, 1e-6);
0654   BOOST_CHECK_CLOSE(track.covariance().determinant(), 1e-27, 4e0);
0655
0656   // Convergence
0657   BOOST_CHECK_LT(
0658       (track.template component<
0659           std::size_t,
0660           hashString(Experimental::Gx2fConstants::gx2fnUpdateColumn)>()),
0661       10);
0662
0663   ACTS_INFO("*** Test: relChi2changeCutOff -- Finish");
0664 }
0665
0666 BOOST_AUTO_TEST_CASE(DidNotConverge) {
0667   ACTS_INFO("*** Test: DidNotConverge -- Start");
0668
0669   std::default_random_engine rng(42);
0670
0671   ACTS_DEBUG("Create the detector");
0672   const std::size_t nSurfaces = 5;
0673   Detector detector;
0674   detector.geometry = makeToyDetector(geoCtx, nSurfaces);
0675
0676   ACTS_DEBUG("Set the start parameters for measurement creation and fit");
0677   const auto parametersMeasurements = makeParameters();
0678   const auto startParametersFit = makeParameters(
0679       7_mm, 11_mm, 15_mm, 42_ns, 10_degree, 80_degree, 1_GeV, 1_e);
0680
0681   ACTS_DEBUG("Create the measurements");
0682   // simulation propagator
0683   using SimPropagator =
0684       Acts::Propagator<Acts::StraightLineStepper, Acts::Navigator>;
0685   const SimPropagator simPropagator = makeStraightPropagator(detector.geometry);
0686   const auto measurements =
0687       createMeasurements(simPropagator, geoCtx, magCtx, parametersMeasurements,
0688                          resMapAllPixel, rng);
0689   const auto sourceLinks = prepareSourceLinks(measurements.sourceLinks);
0690   ACTS_VERBOSE("sourceLinks.size() = " << sourceLinks.size());
0691
0692   BOOST_REQUIRE_EQUAL(sourceLinks.size(), nSurfaces);
0693
0694   ACTS_DEBUG("Set up the fitter");
0695   const Surface* rSurface = &parametersMeasurements.referenceSurface();
0696
0697   using RecoStepper = EigenStepper<>;
0698   const auto recoPropagator =
0699       makeConstantFieldPropagator<RecoStepper>(detector.geometry, 0_T);
0700
0701   using RecoPropagator = decltype(recoPropagator);
0702   using Gx2Fitter =
0703       Experimental::Gx2Fitter<RecoPropagator, VectorMultiTrajectory>;
0704   const Gx2Fitter fitter(recoPropagator, gx2fLogger->clone());
0705
0706   Experimental::Gx2FitterExtensions<VectorMultiTrajectory> extensions;
0707   extensions.calibrator
0708       .connect<&testSourceLinkCalibrator<VectorMultiTrajectory>>();
0709   TestSourceLink::SurfaceAccessor surfaceAccessor{*detector.geometry};
0710   extensions.surfaceAccessor
0711       .connect<&TestSourceLink::SurfaceAccessor::operator()>(&surfaceAccessor);
0712
0713   // The relChi2changeCutOff = 0 prevents to stop the fitter after convergence,
0714   // therefore all updates will be done (even if the result does not change).
0715   // Since we didn't break due to convergence, we reach nUpdatesMax and
0716   // therefore fail the fit.
0717   const Experimental::Gx2FitterOptions gx2fOptions(
0718       geoCtx, magCtx, calCtx, extensions, PropagatorPlainOptions(), rSurface,
0719       false, false, FreeToBoundCorrection(false), 6, true, 0);
0720
0721   Acts::TrackContainer tracks{Acts::VectorTrackContainer{},
0722                               Acts::VectorMultiTrajectory{}};
0723
0724   ACTS_DEBUG("Fit the track");
0725   ACTS_VERBOSE("startParameter unsmeared:\n" << parametersMeasurements);
0726   ACTS_VERBOSE("startParameter fit:\n" << startParametersFit);
0727   const auto res = fitter.fit(sourceLinks.begin(), sourceLinks.end(),
0728                               startParametersFit, gx2fOptions, tracks);
0729
0730   BOOST_REQUIRE(!res.ok());
0731   BOOST_CHECK_EQUAL(
0732       res.error(),
0733       Acts::Experimental::GlobalChiSquareFitterError::DidNotConverge);
0734
0735   ACTS_INFO("*** Test: DidNotConverge -- Finish");
0736 }
0737
0738 BOOST_AUTO_TEST_CASE(NotEnoughMeasurements) {
0739   ACTS_INFO("*** Test: NotEnoughMeasurements -- Start");
0740
0741   std::default_random_engine rng(42);
0742
0743   ACTS_DEBUG("Create the detector");
0744   const std::size_t nSurfaces = 2;
0745   Detector detector;
0746   detector.geometry = makeToyDetector(geoCtx, nSurfaces);
0747
0748   ACTS_DEBUG("Set the start parameters for measurement creation and fit");
0749   const auto parametersMeasurements = makeParameters();
0750   const auto startParametersFit = makeParameters(
0751       7_mm, 11_mm, 15_mm, 42_ns, 10_degree, 80_degree, 1_GeV, 1_e);
0752
0753   ACTS_DEBUG("Create the measurements");
0754   // simulation propagator
0755   using SimPropagator =
0756       Acts::Propagator<Acts::StraightLineStepper, Acts::Navigator>;
0757   const SimPropagator simPropagator = makeStraightPropagator(detector.geometry);
0758   const auto measurements =
0759       createMeasurements(simPropagator, geoCtx, magCtx, parametersMeasurements,
0760                          resMapAllPixel, rng);
0761   const auto sourceLinks = prepareSourceLinks(measurements.sourceLinks);
0762   ACTS_VERBOSE("sourceLinks.size() = " << sourceLinks.size());
0763
0764   BOOST_REQUIRE_EQUAL(sourceLinks.size(), nSurfaces);
0765
0766   ACTS_DEBUG("Set up the fitter");
0767   const Surface* rSurface = &parametersMeasurements.referenceSurface();
0768
0769   using RecoStepper = EigenStepper<>;
0770   const auto recoPropagator =
0771       makeConstantFieldPropagator<RecoStepper>(detector.geometry, 0_T);
0772
0773   using RecoPropagator = decltype(recoPropagator);
0774   using Gx2Fitter =
0775       Experimental::Gx2Fitter<RecoPropagator, VectorMultiTrajectory>;
0776   const Gx2Fitter fitter(recoPropagator, gx2fLogger->clone());
0777
0778   Experimental::Gx2FitterExtensions<VectorMultiTrajectory> extensions;
0779   extensions.calibrator
0780       .connect<&testSourceLinkCalibrator<VectorMultiTrajectory>>();
0781   TestSourceLink::SurfaceAccessor surfaceAccessor{*detector.geometry};
0782   extensions.surfaceAccessor
0783       .connect<&TestSourceLink::SurfaceAccessor::operator()>(&surfaceAccessor);
0784
0785   const Experimental::Gx2FitterOptions gx2fOptions(
0786       geoCtx, magCtx, calCtx, extensions, PropagatorPlainOptions(), rSurface,
0787       false, false, FreeToBoundCorrection(false), 6, true, 0);
0788
0789   Acts::TrackContainer tracks{Acts::VectorTrackContainer{},
0790                               Acts::VectorMultiTrajectory{}};
0791
0792   ACTS_DEBUG("Fit the track");
0793   ACTS_VERBOSE("startParameter unsmeared:\n" << parametersMeasurements);
0794   ACTS_VERBOSE("startParameter fit:\n" << startParametersFit);
0795   const auto res = fitter.fit(sourceLinks.begin(), sourceLinks.end(),
0796                               startParametersFit, gx2fOptions, tracks);
0797
0798   BOOST_REQUIRE(!res.ok());
0799   BOOST_CHECK_EQUAL(
0800       res.error(),
0801       Acts::Experimental::GlobalChiSquareFitterError::NotEnoughMeasurements);
0802
0803   ACTS_INFO("*** Test: NotEnoughMeasurements -- Finish");
0804 }
0805
0806 BOOST_AUTO_TEST_CASE(FindHoles) {
0807   ACTS_INFO("*** Test: FindHoles -- Start");
0808
0809   std::default_random_engine rng(42);
0810
0811   ACTS_DEBUG("Create the detector");
0812   //  const std::size_t nSurfaces = 7;
0813   const std::size_t nSurfaces = 8;
0814   Detector detector;
0815   detector.geometry = makeToyDetector(geoCtx, nSurfaces);
0816
0817   ACTS_DEBUG("Set the start parameters for measurement creation and fit");
0818   const auto parametersMeasurements = makeParameters();
0819   const auto startParametersFit = makeParameters(
0820       7_mm, 11_mm, 15_mm, 42_ns, 10_degree, 80_degree, 1_GeV, 1_e);
0821
0822   ACTS_DEBUG("Create the measurements");
0823   using SimPropagator =
0824       Acts::Propagator<Acts::StraightLineStepper, Acts::Navigator>;
0825   const SimPropagator simPropagator = makeStraightPropagator(detector.geometry);
0826   const auto measurements =
0827       createMeasurements(simPropagator, geoCtx, magCtx, parametersMeasurements,
0828                          resMapAllPixel, rng);
0829   auto sourceLinks = prepareSourceLinks(measurements.sourceLinks);
0830   ACTS_VERBOSE("sourceLinks.size() [before] = " << sourceLinks.size());
0831
0832   // We remove the first measurement in the list. This does not create a hole.
0833   sourceLinks.erase(std::next(sourceLinks.begin(), 0));
0834   ACTS_VERBOSE(
0835       "sourceLinks.size() [after first erase] = " << sourceLinks.size());
0836
0837   // We remove the last measurement in the list. This does not create a hole.
0838   sourceLinks.pop_back();
0839   ACTS_VERBOSE("sourceLinks.size() [after pop] = " << sourceLinks.size());
0840
0841   // We remove the second to last measurement in the list. This effectivly
0842   // creates a hole on that surface.
0843   const std::size_t indexHole = sourceLinks.size() - 2;
0844   ACTS_VERBOSE("Remove measurement " << indexHole);
0845   sourceLinks.erase(std::next(sourceLinks.begin(), indexHole));
0846   ACTS_VERBOSE("sourceLinks.size() [after second-to-last erase]= "
0847                << sourceLinks.size());
0848
0849   // We removed 3 measurements
0850   //  const std::size_t nMeasurements = nSurfaces - 2;
0851   const std::size_t nMeasurements = nSurfaces - 3;
0852   BOOST_REQUIRE_EQUAL(sourceLinks.size(), nMeasurements);
0853
0854   ACTS_DEBUG("Set up the fitter");
0855   const Surface* rSurface = &parametersMeasurements.referenceSurface();
0856
0857   using RecoStepper = EigenStepper<>;
0858   const auto recoPropagator =
0859       makeConstantFieldPropagator<RecoStepper>(detector.geometry, 0_T);
0860
0861   using RecoPropagator = decltype(recoPropagator);
0862   using Gx2Fitter =
0863       Experimental::Gx2Fitter<RecoPropagator, VectorMultiTrajectory>;
0864   const Gx2Fitter fitter(recoPropagator, gx2fLogger->clone());
0865
0866   Experimental::Gx2FitterExtensions<VectorMultiTrajectory> extensions;
0867   extensions.calibrator
0868       .connect<&testSourceLinkCalibrator<VectorMultiTrajectory>>();
0869   TestSourceLink::SurfaceAccessor surfaceAccessor{*detector.geometry};
0870   extensions.surfaceAccessor
0871       .connect<&TestSourceLink::SurfaceAccessor::operator()>(&surfaceAccessor);
0872
0873   const Experimental::Gx2FitterOptions gx2fOptions(
0874       geoCtx, magCtx, calCtx, extensions, PropagatorPlainOptions(), rSurface,
0875       false, false, FreeToBoundCorrection(false), 20, true, 1e-5);
0876
0877   Acts::TrackContainer tracks{Acts::VectorTrackContainer{},
0878                               Acts::VectorMultiTrajectory{}};
0879
0880   ACTS_DEBUG("Fit the track");
0881   ACTS_VERBOSE("startParameter unsmeared:\n" << parametersMeasurements);
0882   ACTS_VERBOSE("startParameter fit:\n" << startParametersFit);
0883   const auto res = fitter.fit(sourceLinks.begin(), sourceLinks.end(),
0884                               startParametersFit, gx2fOptions, tracks);
0885
0886   BOOST_REQUIRE(res.ok());
0887
0888   const auto& track = *res;
0889
0890   // -1, because the index starts at 0
0891   // -2, because the first and the last surface are not part of the track
0892   BOOST_CHECK_EQUAL(track.tipIndex(), nSurfaces - 1 - 2);
0893   BOOST_CHECK(track.hasReferenceSurface());
0894
0895   // Track quantities
0896   CHECK_CLOSE_ABS(track.chi2(), 6.5, 2.);
0897   BOOST_CHECK_EQUAL(track.nDoF(), 10u);
0898   BOOST_CHECK_EQUAL(track.nHoles(), 1u);
0899   BOOST_CHECK_EQUAL(track.nMeasurements(), nMeasurements);
0900   BOOST_CHECK_EQUAL(track.nSharedHits(), 0u);
0901   BOOST_CHECK_EQUAL(track.nOutliers(), 0u);
0902
0903   // Parameters
0904   // We need quite coarse checks here, since on different builds
0905   // the created measurements differ in the randomness
0906   BOOST_CHECK_CLOSE(track.parameters()[eBoundLoc0], -11., 7e0);
0907   BOOST_CHECK_CLOSE(track.parameters()[eBoundLoc1], -15., 6e0);
0908   BOOST_CHECK_CLOSE(track.parameters()[eBoundPhi], 1e-5, 1e3);
0909   BOOST_CHECK_CLOSE(track.parameters()[eBoundTheta], M_PI / 2, 1e-3);
0910   BOOST_CHECK_EQUAL(track.parameters()[eBoundQOverP], 1);
0911   BOOST_CHECK_CLOSE(track.parameters()[eBoundTime], 12591.2832360000, 1e-6);
0912   BOOST_CHECK_CLOSE(track.covariance().determinant(), 4.7e-28, 2e0);
0913
0914   ACTS_INFO("*** Test: FindHoles -- Finish");
0915 }
0916 BOOST_AUTO_TEST_SUITE_END()
0917 }  // namespace Acts::Test