Plugins/Mille/ActsMilleReadRecordSequence.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/Surfaces/PlaneSurface.hpp"
0014 #include "Acts/Surfaces/RectangleBounds.hpp"
0015 #include "Acts/Surfaces/Surface.hpp"
0016 #include "Acts/Utilities/Logger.hpp"
0017 #include "ActsAlignment/Kernel/detail/AlignmentEngine.hpp"
0018 #include "ActsPlugins/Mille/ActsToMille.hpp"
0019
0020 #include <algorithm>
0021 #include <map>
0022 #include <memory>
0023 #include <set>
0024 #include <string>
0025 #include <unordered_map>
0026 #include <vector>
0027
0028 #include "Mille/MilleFactory.h"
0029
0030 // Reading several records into the same TrackAlignmentState, as
0031 // ActsSolverFromMille does. Tracks do not all cross the same surfaces
0032 // (holes, tracks leaving the detector): every record must be unpacked on its
0033 // own, with the derivative matrix sized and indexed for its own surfaces.
0034
0035 namespace {
0036
0037 constexpr std::size_t nSurfaces = 6;
0038
0039 /// one measurement: residual, the surface it is on, the alignment labels
0040 /// (dof 0..5) with non-zero derivative
0041 struct Hit {
0042   float residual;
0043   std::size_t surface;
0044   std::vector<std::size_t> dofs = {0, 1, 2, 3, 4, 5};
0045 };
0046
0047 int label(std::size_t surface, std::size_t dof) {
0048   return static_cast<int>(surface * Acts::eAlignmentSize + dof + 1);
0049 }
0050 double derivative(std::size_t surface, std::size_t dof) {
0051   return 0.25 + 0.5 * surface + 0.125 * dof;
0052 }
0053
0054 void writeRecord(Mille::MilleRecord& record, const std::vector<Hit>& hits) {
0055   unsigned int local = 1;
0056   for (const Hit& hit : hits) {
0057     std::vector<int> labels;
0058     std::vector<double> derivatives;
0059     for (std::size_t dof : hit.dofs) {
0060       labels.push_back(label(hit.surface, dof));
0061       derivatives.push_back(derivative(hit.surface, dof));
0062     }
0063     // one local parameter per measurement, so that the local system is
0064     // well defined
0065     record.addData(hit.residual, 0.1f, std::vector<unsigned int>{local},
0066                    std::vector<double>{1.0}, labels, derivatives);
0067     ++local;
0068   }
0069   record.writeRecord();
0070 }
0071
0072 }  // namespace
0073
0074 BOOST_AUTO_TEST_SUITE(ActsMilleReadRecordSequenceTests)
0075
0076 BOOST_AUTO_TEST_CASE(ReuseStateAcrossRecords) {
0077   ACTS_LOCAL_LOGGER(
0078       Acts::getDefaultLogger("ReuseStateAcrossRecords", Acts::Logging::INFO));
0079   // the alignable surfaces and their indices in the geometry
0080   std::vector<std::shared_ptr<Acts::Surface>> surfaces;
0081   std::unordered_map<const Acts::Surface*, std::size_t> idxedAlignSurfaces;
0082   for (std::size_t s = 0; s < nSurfaces; ++s) {
0083     surfaces.push_back(Acts::Surface::makeShared<Acts::PlaneSurface>(
0084         Acts::Transform3::Identity(),
0085         std::make_shared<Acts::RectangleBounds>(10., 10.)));
0086     idxedAlignSurfaces.emplace(surfaces.back().get(), s);
0087   }
0088
0089   // record sequence: a full track first, then tracks with fewer surfaces
0090   const std::vector<std::vector<Hit>> records = {
0091       // full track
0092       {{0.01f, 0}, {0.02f, 1}, {0.03f, 2}, {0.04f, 3}, {0.05f, 4}, {0.06f, 5}},
0093       // track on the last surface only
0094       {{0.07f, 5}},
0095       // track with a hole on surface 2
0096       {{0.01f, 0}, {0.02f, 1}, {0.04f, 3}, {0.05f, 4}, {0.06f, 5}},
0097       // track leaving after surface 1
0098       {{0.01f, 0}, {0.02f, 1}},
0099       // surface 4 seen only through dof 2..5: its first label is missing
0100       {{0.03f, 3}, {0.05f, 4, {2, 3, 4, 5}}},
0101       // a hit with a residual of exactly zero is still a hit
0102       {{0.f, 1}, {0.02f, 2}},
0103   };
0104
0105   const std::string fname = "ActsMilleReadRecordSequence.dat";
0106   {
0107     std::unique_ptr<Mille::MilleRecord> out = Mille::spawnMilleRecord(fname);
0108     BOOST_REQUIRE(out != nullptr);
0109     for (const auto& hits : records) {
0110       writeRecord(*out, hits);
0111     }
0112   }
0113
0114   auto reader = Mille::spawnMilleReader(fname);
0115   BOOST_REQUIRE(reader != nullptr);
0116   BOOST_REQUIRE(reader->open(fname));
0117
0118   // one state for all records, as in ActsSolverFromMille
0119   ActsAlignment::detail::TrackAlignmentState state;
0120   for (const auto& hits : records) {
0121     BOOST_REQUIRE(ActsPlugins::ActsToMille::unpackMilleRecord(
0122                       *reader, state, idxedAlignSurfaces, logger()) ==
0123                   Mille::MilleDecoder::ReadResult::OK);
0124
0125     std::set<std::size_t> onTrack;
0126     for (const Hit& hit : hits) {
0127       onTrack.insert(hit.surface);
0128     }
0129     // exactly the surfaces of this record, numbered 0..k-1 in label order
0130     BOOST_CHECK_EQUAL(state.alignedSurfaces.size(), onTrack.size());
0131     BOOST_CHECK_EQUAL(state.alignmentDof,
0132                       Acts::eAlignmentSize * onTrack.size());
0133     BOOST_CHECK_EQUAL(state.measurementDim, hits.size());
0134     BOOST_CHECK_EQUAL(state.alignmentToResidualDerivative.rows(),
0135                       static_cast<long>(hits.size()));
0136     BOOST_CHECK_EQUAL(state.alignmentToResidualDerivative.cols(),
0137                       static_cast<long>(state.alignmentDof));
0138     std::map<std::size_t, std::size_t> internalOf;
0139     std::size_t expectedInternal = 0;
0140     for (std::size_t s : onTrack) {
0141       internalOf[s] = expectedInternal++;
0142     }
0143     for (const auto& [surface, indices] : state.alignedSurfaces) {
0144       const auto [globalIndex, internalIndex] = indices;
0145       BOOST_CHECK_EQUAL(globalIndex, idxedAlignSurfaces.at(surface));
0146       BOOST_REQUIRE(internalOf.contains(globalIndex));
0147       BOOST_CHECK_EQUAL(internalIndex, internalOf.at(globalIndex));
0148     }
0149     // every derivative in the column of its label, zero elsewhere
0150     for (std::size_t iMeas = 0; iMeas < hits.size(); ++iMeas) {
0151       const Hit& hit = hits[iMeas];
0152       for (std::size_t s : onTrack) {
0153         for (std::size_t dof = 0; dof < Acts::eAlignmentSize; ++dof) {
0154           const bool written =
0155               s == hit.surface &&
0156               std::ranges::find(hit.dofs, dof) != hit.dofs.end();
0157           const double value = state.alignmentToResidualDerivative(
0158               iMeas, internalOf.at(s) * Acts::eAlignmentSize + dof);
0159           if (written) {
0160             BOOST_CHECK_CLOSE(value, derivative(s, dof), 1e-4);
0161           } else {
0162             BOOST_CHECK_EQUAL(value, 0.);
0163           }
0164         }
0165       }
0166     }
0167   }
0168   BOOST_CHECK(ActsPlugins::ActsToMille::unpackMilleRecord(
0169                   *reader, state, idxedAlignSurfaces, logger()) ==
0170               Mille::MilleDecoder::ReadResult::atEof);
0171 }
0172
0173 BOOST_AUTO_TEST_CASE(UnknownSurfaceIsReadError) {
0174   ACTS_LOCAL_LOGGER(
0175       Acts::getDefaultLogger("UnknownSurfaceIsReadError", Acts::Logging::INFO));
0176   std::vector<std::shared_ptr<Acts::Surface>> surfaces;
0177   std::unordered_map<const Acts::Surface*, std::size_t> idxedAlignSurfaces;
0178   for (std::size_t s = 0; s < 2; ++s) {
0179     surfaces.push_back(Acts::Surface::makeShared<Acts::PlaneSurface>(
0180         Acts::Transform3::Identity(),
0181         std::make_shared<Acts::RectangleBounds>(10., 10.)));
0182     idxedAlignSurfaces.emplace(surfaces.back().get(), s);
0183   }
0184
0185   const std::string fname = "ActsMilleReadUnknownSurface.dat";
0186   {
0187     std::unique_ptr<Mille::MilleRecord> out = Mille::spawnMilleRecord(fname);
0188     BOOST_REQUIRE(out != nullptr);
0189     writeRecord(*out, {{0.01f, 0}, {0.02f, 1}});
0190     // surface index 4 is not in the indexed list
0191     writeRecord(*out, {{0.01f, 0}, {0.02f, 4}});
0192   }
0193   auto reader = Mille::spawnMilleReader(fname);
0194   BOOST_REQUIRE(reader != nullptr);
0195   BOOST_REQUIRE(reader->open(fname));
0196
0197 #ifdef ACTS_ENABLE_LOG_FAILURE_THRESHOLD
0198   // the read error below is expected: do not fail on it
0199   Acts::Logging::ScopedFailureThreshold threshold{Acts::Logging::FATAL};
0200 #endif
0201   ActsAlignment::detail::TrackAlignmentState state;
0202   BOOST_REQUIRE(ActsPlugins::ActsToMille::unpackMilleRecord(
0203                     *reader, state, idxedAlignSurfaces, logger()) ==
0204                 Mille::MilleDecoder::ReadResult::OK);
0205   BOOST_CHECK_EQUAL(state.alignmentDof, 2 * Acts::eAlignmentSize);
0206   BOOST_CHECK(ActsPlugins::ActsToMille::unpackMilleRecord(
0207                   *reader, state, idxedAlignSurfaces, logger()) ==
0208               Mille::MilleDecoder::ReadResult::error);
0209   // a failed read leaves the state untouched
0210   BOOST_CHECK_EQUAL(state.alignmentDof, 2 * Acts::eAlignmentSize);
0211   BOOST_CHECK_EQUAL(state.alignedSurfaces.size(), 2u);
0212 }
0213
0214 BOOST_AUTO_TEST_CASE(GlobalsWithoutLocalsIsReadError) {
0215   ACTS_LOCAL_LOGGER(Acts::getDefaultLogger("GlobalsWithoutLocalsIsReadError",
0216                                            Acts::Logging::INFO));
0217   std::vector<std::shared_ptr<Acts::Surface>> surfaces;
0218   std::unordered_map<const Acts::Surface*, std::size_t> idxedAlignSurfaces;
0219   for (std::size_t s = 0; s < 2; ++s) {
0220     surfaces.push_back(Acts::Surface::makeShared<Acts::PlaneSurface>(
0221         Acts::Transform3::Identity(),
0222         std::make_shared<Acts::RectangleBounds>(10., 10.)));
0223     idxedAlignSurfaces.emplace(surfaces.back().get(), s);
0224   }
0225
0226   const std::string fname = "ActsMilleReadGlobalsWithoutLocals.dat";
0227   {
0228     std::unique_ptr<Mille::MilleRecord> out = Mille::spawnMilleRecord(fname);
0229     BOOST_REQUIRE(out != nullptr);
0230     writeRecord(*out, {{0.01f, 0}, {0.02f, 1}});
0231     // second record: the hit on surface 1 has alignment derivatives but its
0232     // only local derivative is zero, which Mille does not write
0233     out->addData(0.01f, 0.1f, std::vector<unsigned int>{1u},
0234                  std::vector<double>{1.0}, std::vector<int>{label(0, 0)},
0235                  std::vector<double>{derivative(0, 0)});
0236     out->addData(0.02f, 0.1f, std::vector<unsigned int>{2u},
0237                  std::vector<double>{0.0}, std::vector<int>{label(1, 0)},
0238                  std::vector<double>{derivative(1, 0)});
0239     out->writeRecord();
0240   }
0241   auto reader = Mille::spawnMilleReader(fname);
0242   BOOST_REQUIRE(reader != nullptr);
0243   BOOST_REQUIRE(reader->open(fname));
0244
0245 #ifdef ACTS_ENABLE_LOG_FAILURE_THRESHOLD
0246   // the read error below is expected: do not fail on it
0247   Acts::Logging::ScopedFailureThreshold threshold{Acts::Logging::FATAL};
0248 #endif
0249   ActsAlignment::detail::TrackAlignmentState state;
0250   BOOST_REQUIRE(ActsPlugins::ActsToMille::unpackMilleRecord(
0251                     *reader, state, idxedAlignSurfaces, logger()) ==
0252                 Mille::MilleDecoder::ReadResult::OK);
0253   BOOST_CHECK(ActsPlugins::ActsToMille::unpackMilleRecord(
0254                   *reader, state, idxedAlignSurfaces, logger()) ==
0255               Mille::MilleDecoder::ReadResult::error);
0256   // a failed read leaves the state untouched
0257   BOOST_CHECK_EQUAL(state.alignmentDof, 2 * Acts::eAlignmentSize);
0258 }
0259
0260 BOOST_AUTO_TEST_SUITE_END()