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