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