Alignment/Kernel/AlignmentTests.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/Units.hpp"
0012 #include "Acts/EventData/TrackParameters.hpp"
0013 #include "Acts/EventData/detail/TestSourceLink.hpp"
0014 #include "Acts/Geometry/CuboidVolumeBounds.hpp"
0015 #include "Acts/Geometry/CuboidVolumeBuilder.hpp"
0016 #include "Acts/Geometry/GeometryContext.hpp"
0017 #include "Acts/Geometry/LayerArrayCreator.hpp"
0018 #include "Acts/Geometry/LayerCreator.hpp"
0019 #include "Acts/Geometry/PlaneLayer.hpp"
0020 #include "Acts/Geometry/SurfaceArrayCreator.hpp"
0021 #include "Acts/Geometry/TrackingGeometry.hpp"
0022 #include "Acts/Geometry/TrackingGeometryBuilder.hpp"
0023 #include "Acts/Geometry/TrackingVolume.hpp"
0024 #include "Acts/MagneticField/ConstantBField.hpp"
0025 #include "Acts/MagneticField/MagneticFieldContext.hpp"
0026 #include "Acts/Material/HomogeneousSurfaceMaterial.hpp"
0027 #include "Acts/Material/ISurfaceMaterial.hpp"
0028 #include "Acts/Propagator/EigenStepper.hpp"
0029 #include "Acts/Propagator/Navigator.hpp"
0030 #include "Acts/Propagator/Propagator.hpp"
0031 #include "Acts/Propagator/StraightLineStepper.hpp"
0032 #include "Acts/Surfaces/PlaneSurface.hpp"
0033 #include "Acts/Surfaces/RectangleBounds.hpp"
0034 #include "Acts/Tests/CommonHelpers/DetectorElementStub.hpp"
0035 #include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
0036 #include "Acts/Tests/CommonHelpers/MeasurementsCreator.hpp"
0037 #include "Acts/Tests/CommonHelpers/PredefinedMaterials.hpp"
0038 #include "Acts/TrackFitting/GainMatrixSmoother.hpp"
0039 #include "Acts/TrackFitting/GainMatrixUpdater.hpp"
0040 #include "Acts/TrackFitting/KalmanFitter.hpp"
0041 #include "Acts/TrackFitting/detail/KalmanGlobalCovariance.hpp"
0042 #include "Acts/Utilities/CalibrationContext.hpp"
0043 #include "ActsAlignment/Kernel/Alignment.hpp"
0044
0045 #include <cmath>
0046 #include <random>
0047 #include <string>
0048
0049 namespace {
0050 using namespace Acts;
0051 using namespace ActsAlignment;
0052 using namespace Acts::Test;
0053 using namespace Acts::detail::Test;
0054 using namespace Acts::UnitLiterals;
0055
0056 using StraightPropagator =
0057     Acts::Propagator<Acts::StraightLineStepper, Acts::Navigator>;
0058 using ConstantFieldStepper = Acts::EigenStepper<>;
0059 using ConstantFieldPropagator =
0060     Acts::Propagator<ConstantFieldStepper, Acts::Navigator>;
0061
0062 using KalmanUpdater = Acts::GainMatrixUpdater;
0063 using KalmanSmoother = Acts::GainMatrixSmoother;
0064 using KalmanFitterType =
0065     Acts::KalmanFitter<ConstantFieldPropagator, VectorMultiTrajectory>;
0066
0067 KalmanUpdater kfUpdater;
0068 KalmanSmoother kfSmoother;
0069
0070 // Create a test context
0071 const GeometryContext geoCtx;
0072 const MagneticFieldContext magCtx;
0073 const CalibrationContext calCtx;
0074
0075 std::normal_distribution<double> normalDist(0., 1.);
0076 std::default_random_engine rng(42);
0077
0078 KalmanFitterExtensions<VectorMultiTrajectory> getExtensions() {
0079   KalmanFitterExtensions<VectorMultiTrajectory> extensions;
0080   extensions.calibrator
0081       .connect<&testSourceLinkCalibrator<VectorMultiTrajectory>>();
0082   extensions.updater.connect<&KalmanUpdater::operator()<VectorMultiTrajectory>>(
0083       &kfUpdater);
0084   extensions.smoother
0085       .connect<&KalmanSmoother::operator()<VectorMultiTrajectory>>(&kfSmoother);
0086   return extensions;
0087 }
0088
0089 ///
0090 /// @brief Construct a telescope-like detector
0091 ///
0092 struct TelescopeDetector {
0093   /// Default constructor for the Cubic tracking geometry
0094   ///
0095   /// @param gctx the geometry context for this geometry at building time
0096   TelescopeDetector(std::reference_wrapper<const GeometryContext> gctx)
0097       : geoContext(gctx) {
0098     // Construct the rotation
0099     rotation.col(0) = Acts::Vector3(0, 0, -1);
0100     rotation.col(1) = Acts::Vector3(0, 1, 0);
0101     rotation.col(2) = Acts::Vector3(1, 0, 0);
0102
0103     // Boundaries of the surfaces
0104     rBounds = std::make_shared<const RectangleBounds>(0.1_m, 0.1_m);
0105
0106     // Material of the surfaces
0107     MaterialSlab matProp(makeSilicon(), 80_um);
0108
0109     surfaceMaterial = std::make_shared<HomogeneousSurfaceMaterial>(matProp);
0110   }
0111
0112   ///
0113   /// Call operator to build the standard cubic tracking geometry
0114   ///
0115   std::shared_ptr<const TrackingGeometry> operator()() {
0116     using namespace UnitLiterals;
0117
0118     unsigned int nLayers = 6;
0119     std::vector<double> positions = {-500_mm, -300_mm, -100_mm,
0120                                      100_mm,  300_mm,  500_mm};
0121     auto length = positions.back() - positions.front();
0122
0123     std::vector<LayerPtr> layers(nLayers);
0124     for (unsigned int i = 0; i < nLayers; ++i) {
0125       // The transform
0126       Translation3 trans(0., 0., positions[i]);
0127       Transform3 trafo(rotation * trans);
0128       auto detElement = std::make_shared<DetectorElementStub>(
0129           trafo, rBounds, 1._um, surfaceMaterial);
0130       // The surface is not right!!!
0131       auto surface = detElement->surface().getSharedPtr();
0132       // Add it to the event store
0133       detectorStore.push_back(std::move(detElement));
0134       std::unique_ptr<SurfaceArray> surArray(new SurfaceArray(surface));
0135       // The layer thickness should not be too large
0136       layers[i] =
0137           PlaneLayer::create(trafo, rBounds, std::move(surArray),
0138                              1._mm);  // Associate the layer to the surface
0139       auto mutableSurface = const_cast<Surface*>(surface.get());
0140       mutableSurface->associateLayer(*layers[i]);
0141     }
0142
0143     // The volume transform
0144     Translation3 transVol(0, 0, 0);
0145     Transform3 trafoVol(rotation * transVol);
0146     auto boundsVol = std::make_shared<CuboidVolumeBounds>(
0147         rBounds->halfLengthX() + 10._mm, rBounds->halfLengthY() + 10._mm,
0148         length + 10._mm);
0149
0150     LayerArrayCreator::Config lacConfig;
0151     LayerArrayCreator layArrCreator(
0152         lacConfig, getDefaultLogger("LayerArrayCreator", Logging::INFO));
0153     LayerVector layVec;
0154     for (unsigned int i = 0; i < nLayers; i++) {
0155       layVec.push_back(layers[i]);
0156     }
0157
0158     // Create the layer array
0159     std::unique_ptr<const LayerArray> layArr(layArrCreator.layerArray(
0160         geoContext, layVec, positions.front() - 2._mm, positions.back() + 2._mm,
0161         BinningType::arbitrary, AxisDirection::AxisX));
0162
0163     // Build the tracking volume
0164     auto trackVolume = std::make_shared<TrackingVolume>(
0165         trafoVol, boundsVol, nullptr, std::move(layArr), nullptr,
0166         MutableTrackingVolumeVector{}, "Telescope");
0167
0168     return std::make_shared<const TrackingGeometry>(trackVolume);
0169   }
0170
0171   RotationMatrix3 rotation = RotationMatrix3::Identity();
0172   std::shared_ptr<const RectangleBounds> rBounds = nullptr;
0173   std::shared_ptr<const ISurfaceMaterial> surfaceMaterial = nullptr;
0174
0175   std::vector<std::shared_ptr<DetectorElementStub>> detectorStore;
0176
0177   std::reference_wrapper<const GeometryContext> geoContext;
0178 };
0179
0180 // Construct a straight-line propagator.
0181 StraightPropagator makeStraightPropagator(
0182     std::shared_ptr<const TrackingGeometry> geo) {
0183   Navigator::Config cfg{std::move(geo)};
0184   cfg.resolvePassive = false;
0185   cfg.resolveMaterial = true;
0186   cfg.resolveSensitive = true;
0187   Navigator navigator(cfg);
0188   StraightLineStepper stepper;
0189   return StraightPropagator(stepper, std::move(navigator));
0190 }
0191
0192 // Construct a propagator using a constant magnetic field along z.
0193 ConstantFieldPropagator makeConstantFieldPropagator(
0194     std::shared_ptr<const TrackingGeometry> geo, double bz,
0195     std::unique_ptr<const Logger> logger) {
0196   Navigator::Config cfg{std::move(geo)};
0197   cfg.resolvePassive = false;
0198   cfg.resolveMaterial = true;
0199   cfg.resolveSensitive = true;
0200   Navigator navigator(cfg, logger->cloneWithSuffix("Nav"));
0201   auto field = std::make_shared<ConstantBField>(Vector3(0.0, 0.0, bz));
0202   ConstantFieldStepper stepper(std::move(field));
0203   return ConstantFieldPropagator(std::move(stepper), std::move(navigator),
0204                                  logger->cloneWithSuffix("Prop"));
0205 }
0206
0207 // Construct initial track parameters.
0208 CurvilinearTrackParameters makeParameters() {
0209   // create covariance matrix from reasonable standard deviations
0210   BoundVector stddev;
0211   stddev[eBoundLoc0] = 100_um;
0212   stddev[eBoundLoc1] = 100_um;
0213   stddev[eBoundTime] = 25_ns;
0214   stddev[eBoundPhi] = 0.5_degree;
0215   stddev[eBoundTheta] = 0.5_degree;
0216   stddev[eBoundQOverP] = 1 / 100_GeV;
0217   BoundSquareMatrix cov = stddev.cwiseProduct(stddev).asDiagonal();
0218
0219   auto loc0 = 0. + stddev[eBoundLoc0] * normalDist(rng);
0220   auto loc1 = 0. + stddev[eBoundLoc1] * normalDist(rng);
0221   auto t = 42_ns + stddev[eBoundTime] * normalDist(rng);
0222   auto phi = 0_degree + stddev[eBoundPhi] * normalDist(rng);
0223   auto theta = 90_degree + stddev[eBoundTheta] * normalDist(rng);
0224   auto qOverP = 1_e / 1_GeV + stddev[eBoundQOverP] * normalDist(rng);
0225
0226   // define a track in the transverse plane along x
0227   Vector4 mPos4(-1_m, loc0, loc1, t);
0228
0229   return CurvilinearTrackParameters(mPos4, phi, theta, qOverP, cov,
0230                                     ParticleHypothesis::pion());
0231 }
0232
0233 // detector resolutions
0234 const MeasurementResolution resPixel = {MeasurementType::eLoc01,
0235                                         {30_um, 50_um}};
0236 const MeasurementResolutionMap resolutions = {
0237     {GeometryIdentifier(), resPixel},
0238 };
0239
0240 struct KalmanFitterInputTrajectory {
0241   // The source links
0242   std::vector<TestSourceLink> sourceLinks;
0243   // The start parameters
0244   std::optional<CurvilinearTrackParameters> startParameters;
0245 };
0246
0247 ///
0248 /// Function to create trajectories for kalman fitter
0249 ///
0250 std::vector<KalmanFitterInputTrajectory> createTrajectories(
0251     std::shared_ptr<const TrackingGeometry> geo, std::size_t nTrajectories) {
0252   // simulation propagator
0253   const auto simPropagator = makeStraightPropagator(std::move(geo));
0254
0255   std::vector<KalmanFitterInputTrajectory> trajectories;
0256   trajectories.reserve(nTrajectories);
0257
0258   for (unsigned int iTrack = 0; iTrack < nTrajectories; iTrack++) {
0259     auto start = makeParameters();
0260     // Launch and collect - the measurements
0261     auto measurements = createMeasurements(simPropagator, geoCtx, magCtx, start,
0262                                            resolutions, rng);
0263
0264     // Extract measurements from result of propagation.
0265     KalmanFitterInputTrajectory traj;
0266     traj.startParameters = start;
0267     traj.sourceLinks = measurements.sourceLinks;
0268
0269     trajectories.push_back(std::move(traj));
0270   }
0271   return trajectories;
0272 }
0273 }  // namespace
0274
0275 ///
0276 /// @brief Unit test for KF-based alignment algorithm
0277 ///
0278 BOOST_AUTO_TEST_CASE(ZeroFieldKalmanAlignment) {
0279   // Build detector
0280   TelescopeDetector detector(geoCtx);
0281   const auto geometry = detector();
0282
0283   // reconstruction propagator and fitter
0284   auto kfLogger = getDefaultLogger("KalmanFilter", Logging::INFO);
0285   const auto kfZeroPropagator =
0286       makeConstantFieldPropagator(geometry, 0_T, std::move(kfLogger));
0287   auto kfZero = KalmanFitterType(kfZeroPropagator);
0288
0289   // alignment
0290   auto alignLogger = getDefaultLogger("Alignment", Logging::INFO);
0291   const auto alignZero = Alignment(std::move(kfZero), std::move(alignLogger));
0292
0293   // Create 10 trajectories
0294   const auto& trajectories = createTrajectories(geometry, 10);
0295
0296   // Construct the KalmanFitter options
0297
0298   auto extensions = getExtensions();
0299   TestSourceLink::SurfaceAccessor surfaceAccessor{*geometry};
0300   extensions.surfaceAccessor
0301       .connect<&TestSourceLink::SurfaceAccessor::operator()>(&surfaceAccessor);
0302   KalmanFitterOptions kfOptions(geoCtx, magCtx, calCtx, extensions,
0303                                 PropagatorPlainOptions(geoCtx, magCtx));
0304
0305   // Construct a non-updating alignment updater
0306   AlignedTransformUpdater voidAlignUpdater =
0307       [](DetectorElementBase* /*element*/, const GeometryContext& /*gctx*/,
0308          const Transform3& /*transform*/) { return true; };
0309
0310   // Construct the alignment options
0311   AlignmentOptions<KalmanFitterOptions<VectorMultiTrajectory>> alignOptions(
0312       kfOptions, voidAlignUpdater);
0313   alignOptions.maxIterations = 1;
0314
0315   // Set the surfaces to be aligned (fix the layer 8)
0316   unsigned int iSurface = 0;
0317   std::unordered_map<const Surface*, std::size_t> idxedAlignSurfaces;
0318   // Loop over the detector elements
0319   for (auto& det : detector.detectorStore) {
0320     const auto& surface = det->surface();
0321     if (surface.geometryId().layer() != 8) {
0322       alignOptions.alignedDetElements.push_back(det.get());
0323       idxedAlignSurfaces.emplace(&surface, iSurface);
0324       iSurface++;
0325     }
0326   }
0327
0328   // Test the method to evaluate alignment state for a single track
0329   const auto& inputTraj = trajectories.front();
0330   kfOptions.referenceSurface = &(*inputTraj.startParameters).referenceSurface();
0331
0332   auto evaluateRes = alignZero.evaluateTrackAlignmentState(
0333       kfOptions.geoContext, inputTraj.sourceLinks, *inputTraj.startParameters,
0334       kfOptions, idxedAlignSurfaces, AlignmentMask::All);
0335   BOOST_CHECK(evaluateRes.ok());
0336
0337   const auto& alignState = evaluateRes.value();
0338   CHECK_CLOSE_ABS(alignState.chi2 / alignState.alignmentDof, 0.5, 1);
0339
0340   // Check the dimensions
0341   BOOST_CHECK_EQUAL(alignState.measurementDim, 12);
0342   BOOST_CHECK_EQUAL(alignState.trackParametersDim, 36);
0343   // Check the alignment dof
0344   BOOST_CHECK_EQUAL(alignState.alignmentDof, 30);
0345   BOOST_CHECK_EQUAL(alignState.alignedSurfaces.size(), 5);
0346   // Check the measurements covariance
0347   BOOST_CHECK_EQUAL(alignState.measurementCovariance.rows(), 12);
0348   const SquareMatrix2 measCov =
0349       alignState.measurementCovariance.block<2, 2>(2, 2);
0350   SquareMatrix2 cov2D;
0351   cov2D << 30_um * 30_um, 0, 0, 50_um * 50_um;
0352   CHECK_CLOSE_ABS(measCov, cov2D, 1e-10);
0353   // Check the track parameters covariance matrix. Its rows/columns scales
0354   // with the number of measurement states
0355   BOOST_CHECK_EQUAL(alignState.trackParametersCovariance.rows(), 36);
0356   // Check the projection matrix
0357   BOOST_CHECK_EQUAL(alignState.projectionMatrix.rows(), 12);
0358   BOOST_CHECK_EQUAL(alignState.projectionMatrix.cols(), 36);
0359   const ActsMatrix<2, 6> proj = alignState.projectionMatrix.block<2, 6>(0, 0);
0360   const ActsMatrix<2, 6> refProj = ActsMatrix<2, 6>::Identity();
0361   CHECK_CLOSE_ABS(proj, refProj, 1e-10);
0362   // Check the residual
0363   BOOST_CHECK_EQUAL(alignState.residual.size(), 12);
0364   // Check the residual covariance
0365   BOOST_CHECK_EQUAL(alignState.residualCovariance.rows(), 12);
0366   // Check the alignment to residual derivative
0367   BOOST_CHECK_EQUAL(alignState.alignmentToResidualDerivative.rows(), 12);
0368   BOOST_CHECK_EQUAL(alignState.alignmentToResidualDerivative.cols(), 30);
0369   // Check the chi2 derivative
0370   BOOST_CHECK_EQUAL(alignState.alignmentToChi2Derivative.size(), 30);
0371   BOOST_CHECK_EQUAL(alignState.alignmentToChi2SecondDerivative.rows(), 30);
0372
0373   // Test the align method
0374   std::vector<std::vector<TestSourceLink>> trajCollection;
0375   trajCollection.reserve(10);
0376   std::vector<CurvilinearTrackParameters> sParametersCollection;
0377   sParametersCollection.reserve(10);
0378   for (const auto& traj : trajectories) {
0379     trajCollection.push_back(traj.sourceLinks);
0380     sParametersCollection.push_back(*traj.startParameters);
0381   }
0382   auto alignRes =
0383       alignZero.align(trajCollection, sParametersCollection, alignOptions);
0384
0385   // BOOST_CHECK(alignRes.ok());
0386 }