Plugins/Mille/ActsToMillePseudoResiduals.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/Algebra.hpp"
0012 #include "Acts/Definitions/Alignment.hpp"
0013 #include "Acts/Definitions/TrackParametrization.hpp"
0014 #include "Acts/Surfaces/PlaneSurface.hpp"
0015 #include "Acts/Surfaces/RectangleBounds.hpp"
0016 #include "Acts/Surfaces/Surface.hpp"
0017 #include "Acts/Utilities/Logger.hpp"
0018 #include "ActsAlignment/Kernel/detail/AlignmentEngine.hpp"
0019 #include "ActsPlugins/Mille/ActsToMille.hpp"
0020
0021 #include <cmath>
0022 #include <memory>
0023 #include <string>
0024 #include <unordered_map>
0025 #include <vector>
0026
0027 #include "Mille/MilleDecoder.h"
0028 #include "Mille/MilleFactory.h"
0029
0030 // The Kalman correlations are written as pseudo-measurements along the
0031 // eigenvectors of the correlation term. The record is expanded around the
0032 // smoothed track, whose measurement chi2 gradient H^T V^-1 r is in general not
0033 // zero (multiple scattering kinks). The pseudo-measurements must carry
0034 // residuals that cancel it, so that the local fit in Millepede is at its
0035 // minimum on the smoothed track. Otherwise the local fit moves the track and
0036 // the alignment derivatives are biased.
0037
0038 BOOST_AUTO_TEST_SUITE(ActsToMillePseudoResidualsTests)
0039
0040 BOOST_AUTO_TEST_CASE(SmoothedTrackIsLocalFitMinimum) {
0041   // three track states, loc0 and loc1 measured on each
0042   constexpr std::size_t nStates = 3;
0043   constexpr std::size_t nPar = nStates * Acts::eBoundSize;
0044   constexpr std::size_t nMeas = 2 * nStates;
0045   // the last measurement is on a surface that is not aligned
0046   constexpr std::size_t iMeasNotAligned = nMeas - 1;
0047
0048   auto surface = Acts::Surface::makeShared<Acts::PlaneSurface>(
0049       Acts::Transform3::Identity(),
0050       std::make_shared<Acts::RectangleBounds>(10., 10.));
0051
0052   ActsAlignment::detail::TrackAlignmentState state;
0053   state.measurementDim = nMeas;
0054   state.trackParametersDim = nPar;
0055   state.alignmentDof = Acts::eAlignmentSize;
0056   state.alignedSurfaces.emplace(surface.get(), std::make_pair(0u, 0u));
0057
0058   state.measurementCovariance =
0059       Acts::DynamicMatrix::Identity(nMeas, nMeas) * 0.05 * 0.05;
0060   state.projectionMatrix = Acts::DynamicMatrix::Zero(nMeas, nPar);
0061   for (std::size_t s = 0; s < nStates; ++s) {
0062     const std::size_t o = s * Acts::eBoundSize;
0063     state.projectionMatrix(2 * s, o + Acts::eBoundLoc0) = 1.;
0064     state.projectionMatrix(2 * s + 1, o + Acts::eBoundLoc1) = 1.;
0065   }
0066
0067   // Track parameter covariance of a Kalman fit: the inverse of the measurement
0068   // information plus a dense, positive definite correlation term K, which
0069   // couples the parameters on different surfaces.
0070   Acts::DynamicMatrix b(nPar, nPar);
0071   for (std::size_t i = 0; i < nPar; ++i) {
0072     for (std::size_t j = 0; j < nPar; ++j) {
0073       b(i, j) = std::sin(1. + i + 2. * j);
0074     }
0075   }
0076   const Acts::DynamicMatrix correlationTerm =
0077       b.transpose() * b + Acts::DynamicMatrix::Identity(nPar, nPar);
0078   const Acts::DynamicMatrix weightMatMeasurements =
0079       state.projectionMatrix.transpose() *
0080       state.measurementCovariance.inverse() * state.projectionMatrix;
0081   state.trackParametersCovariance =
0082       (weightMatMeasurements + correlationTerm).inverse();
0083
0084   // residuals of a track with kinks: the measurement chi2 gradient is not zero
0085   state.residual = Acts::DynamicVector::Zero(nMeas);
0086   state.alignmentToResidualDerivative =
0087       Acts::DynamicMatrix::Zero(nMeas, state.alignmentDof);
0088   for (std::size_t i = 0; i < nMeas; ++i) {
0089     state.residual(i) = 0.01 * (i % 2 == 0 ? 1. : -2.) * (i + 1);
0090     if (i == iMeasNotAligned) {
0091       continue;
0092     }
0093     for (std::size_t a = 0; a < Acts::eAlignmentSize; ++a) {
0094       state.alignmentToResidualDerivative(i, a) = 0.1 * (a + 1) + 0.01 * i;
0095     }
0096   }
0097   const Acts::DynamicVector measGradient =
0098       state.projectionMatrix.transpose() *
0099       state.measurementCovariance.inverse() * state.residual;
0100   BOOST_REQUIRE(measGradient.norm() > 1.);
0101
0102   auto logger = Acts::getDefaultLogger("ActsToMillePseudoResiduals",
0103                                        Acts::Logging::WARNING);
0104
0105   const std::string fname = "ActsToMillePseudoResiduals.dat";
0106   {
0107     std::unique_ptr<Mille::MilleRecord> out = Mille::spawnMilleRecord(fname);
0108     BOOST_REQUIRE(out != nullptr);
0109     ActsPlugins::ActsToMille::dumpToMille(state, *out, false, *logger);
0110   }
0111
0112   auto reader = Mille::spawnMilleReader(fname);
0113   BOOST_REQUIRE(reader != nullptr);
0114   BOOST_REQUIRE(reader->open(fname));
0115   Mille::MilleDecoder decoder;
0116   std::vector<Mille::MilleMeasurement> measurements;
0117   BOOST_REQUIRE(decoder.decode(*reader, measurements) ==
0118                 Mille::MilleDecoder::ReadResult::OK);
0119   BOOST_REQUIRE_EQUAL(measurements.size(), nMeas + nPar);
0120
0121   // Gradient of the local fit chi2 at the expansion point (the smoothed
0122   // track), summed over the measurements and the pseudo-measurements. Each
0123   // term is compared to the scale of the measurement contribution; the record
0124   // is written in single precision.
0125   Acts::DynamicVector localGradient = Acts::DynamicVector::Zero(nPar);
0126   std::size_t nPseudoNonZero = 0;
0127   for (const auto& measurement : measurements) {
0128     const double weight =
0129         1. / (measurement.uncertainty * measurement.uncertainty);
0130     for (std::size_t iLoc = 0; iLoc < measurement.localLabels.size(); ++iLoc) {
0131       localGradient(measurement.localLabels[iLoc] - 1) +=
0132           weight * measurement.measurement * measurement.localDerivatives[iLoc];
0133     }
0134     if (measurement.globalLabels.empty() &&
0135         measurement.localLabels.size() > 1 && measurement.measurement != 0) {
0136       ++nPseudoNonZero;
0137     }
0138   }
0139   BOOST_CHECK_EQUAL(nPseudoNonZero, nPar);
0140   BOOST_CHECK_SMALL(localGradient.norm() / measGradient.norm(), 1e-5);
0141
0142   // The reader still separates the surface measurements from the
0143   // pseudo-measurements, including the one without alignment derivatives.
0144   auto reader2 = Mille::spawnMilleReader(fname);
0145   BOOST_REQUIRE(reader2 != nullptr);
0146   BOOST_REQUIRE(reader2->open(fname));
0147   const std::unordered_map<const Acts::Surface*, std::size_t>
0148       idxedAlignSurfaces{{surface.get(), 0}};
0149   ActsAlignment::detail::TrackAlignmentState readState;
0150   BOOST_REQUIRE(ActsPlugins::ActsToMille::unpackMilleRecord(
0151                     *reader2, readState, idxedAlignSurfaces, *logger) ==
0152                 Mille::MilleDecoder::ReadResult::OK);
0153   BOOST_CHECK_EQUAL(readState.measurementDim, nMeas);
0154   BOOST_CHECK_EQUAL(readState.trackParametersDim, nPar);
0155   BOOST_REQUIRE_EQUAL(readState.residual.size(), nMeas);
0156   for (std::size_t i = 0; i < nMeas; ++i) {
0157     BOOST_CHECK_CLOSE(readState.residual(i), state.residual(i), 1e-4);
0158   }
0159   BOOST_CHECK(readState.projectionMatrix.isApprox(state.projectionMatrix));
0160 }
0161
0162 BOOST_AUTO_TEST_SUITE_END()