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