AlignmentMillePede/src/ActsSolverFromMille.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 "ActsExamples/AlignmentMillePede/ActsSolverFromMille.hpp"
0010
0011 #include "Acts/Definitions/Alignment.hpp"
0012 #include "Acts/Geometry/GeometryContext.hpp"
0013 #include "Acts/Geometry/GeometryIdentifier.hpp"
0014 #include "Acts/Material/MaterialInteraction.hpp"
0015 #include "Acts/Propagator/Navigator.hpp"
0016 #include "Acts/Propagator/Propagator.hpp"
0017 #include "Acts/Surfaces/Surface.hpp"
0018 #include "Acts/Utilities/Logger.hpp"
0019 #include "ActsAlignment/Kernel/Alignment.hpp"
0020 #include "ActsAlignment/Kernel/detail/AlignmentEngine.hpp"
0021 #include "ActsExamples/Framework/ProcessCode.hpp"
0022 #include "ActsPlugins/Mille/ActsToMille.hpp"
0023
0024 #include <memory>
0025
0026 #include <Mille/MilleFactory.h>
0027
0028 namespace ActsExamples {
0029
0030 ActsSolverFromMille::ActsSolverFromMille(
0031     Config cfg, std::unique_ptr<const Acts::Logger> logger)
0032     : IAlgorithm("ActsSolverFromMille", std::move(logger)),
0033       m_cfg(std::move(cfg)) {
0034   // retrieve tracking geo
0035   m_trackingGeometry = m_cfg.trackingGeometry;
0036
0037   // instantiate the alignment tool instance
0038   Acts::Navigator::Config navcfg{m_cfg.trackingGeometry};
0039   navcfg.resolvePassive = false;
0040   navcfg.resolveMaterial = true;
0041   navcfg.resolveSensitive = true;
0042   Acts::Navigator navigator(navcfg,
0043                             this->logger().cloneWithSuffix("Navigator"));
0044   Stepper stepper{m_cfg.magneticField};
0045   m_align = std::make_shared<Alignment>(
0046       Fitter(Propagator(stepper, Acts::Navigator(navcfg))));
0047 }
0048
0049 ProcessCode ActsSolverFromMille::execute(
0050     const AlgorithmContext& /*ClangShutUp*/) const {
0051   // this algorithm will not do anything at event-time
0052   return ProcessCode::SUCCESS;
0053 }
0054
0055 ProcessCode ActsSolverFromMille::finalize() {
0056   auto milleReader = Mille::spawnMilleReader(m_cfg.milleInput);
0057   auto openStat = milleReader->open(m_cfg.milleInput);
0058   if (!openStat) {
0059     ACTS_FATAL("Failed to read the mille binary " << m_cfg.milleInput);
0060     return ProcessCode::ABORT;
0061   }
0062
0063   // Assign indices to the alignable surfaces
0064
0065   // We wish to have a relation between alignment parameter indices and real
0066   // geometry. The unordered_map does not give us this - so perform a manual
0067   // sorting.
0068   std::vector<std::pair<Acts::GeometryIdentifier, const Acts::Surface*>>
0069       sortedGeo;
0070   ActsAlignment::AlignmentResult alignResult;
0071
0072   sortedGeo.insert(sortedGeo.end(),
0073                    m_trackingGeometry->geoIdSurfaceMap().begin(),
0074                    m_trackingGeometry->geoIdSurfaceMap().end());
0075   std::sort(
0076       sortedGeo.begin(), sortedGeo.end(),
0077       [&](const std::pair<Acts::GeometryIdentifier, const Acts::Surface*>& lhs,
0078           const std::pair<Acts::GeometryIdentifier, const Acts::Surface*>&
0079               rhs) { return (lhs.first.layer() < rhs.first.layer()); });
0080
0081   const Acts::Surface* firstSurf = nullptr;
0082   unsigned int iSurface = 0;
0083   for (auto& [geoID, surface] : sortedGeo) {
0084     // only consider sensitive surfaces
0085     if (!surface->isSensitive()) {
0086       continue;
0087     }
0088     // use the first sensitive surface as trajectory reference in the kalman
0089     if (firstSurf == nullptr) {
0090       firstSurf = surface;
0091     }
0092     if (!m_cfg.fixModules.contains(geoID)) {
0093       alignResult.idxedAlignSurfaces.emplace(surface, iSurface);
0094       iSurface++;
0095     }
0096   }
0097
0098   ACTS_INFO("Performing alignment fit on collected Mille records");
0099   std::vector<ActsAlignment::detail::TrackAlignmentState> alignmentStates;
0100   ActsAlignment::detail::TrackAlignmentState state;
0101   std::size_t iRec = 0;
0102   while (ActsPlugins::ActsToMille::unpackMilleRecord(
0103              *milleReader, state, alignResult.idxedAlignSurfaces) ==
0104          Mille::MilleDecoder::ReadResult::OK) {
0105     if (++iRec % 10000 == 0) {
0106       ACTS_INFO("     Reading input record " << iRec);
0107     }
0108     alignmentStates.push_back(state);
0109   }
0110
0111   /// TODO: Should try a local iteration without track state info.
0112   /// Can use the linearised info in the Track Alignment State
0113   /// to calculate approximate track parameter & residual updates
0114   /// and then repeat the solution. As in Millepede, probably
0115   /// safe to keep the "big matrix" and only update the right hand side.
0116   m_align->calculateAlignmentParameters(alignmentStates, alignResult);
0117
0118   /// in a real experiment, the results would be written out
0119   /// and stored e.g. in a DB file for further use / validation.
0120   /// For this initial demo, we just print them out.
0121
0122   ACTS_INFO("Performed internal alignment. ");
0123   ACTS_INFO(std::setw(16) << "  Tracks used: " << alignmentStates.size());
0124   ACTS_INFO(std::setw(16) << "  avg Chi2/NDF = "
0125                           << alignResult.averageChi2ONdf);
0126   ACTS_INFO(std::setw(16) << "  Chi2   = " << alignResult.chi2);
0127   ACTS_INFO(std::setw(16) << "  delta Chi2   = " << alignResult.deltaChi2);
0128   ACTS_INFO(std::setw(16) << "  Alignment parameter updates: ");
0129   std::vector<std::string> parLabels{"dx", "dy", "dz", "rx", "ry", "rz"};
0130   for (auto [surface, index] : alignResult.idxedAlignSurfaces) {
0131     ACTS_INFO(std::setw(20)
0132               << " Surface with geo ID " << surface->geometryId() << ": ");
0133     for (std::size_t i = 0; i < Acts::eAlignmentSize; ++i) {
0134       std::size_t row = Acts::eAlignmentSize * index + i;
0135       ACTS_INFO(std::setw(20)
0136                 << parLabels[i] << " = " << std::setw(10)
0137                 << alignResult.deltaAlignmentParameters(row) << std::setw(6)
0138                 << " +/- " << std::setw(10)
0139                 << std::sqrt(alignResult.alignmentCovariance(row, row)));
0140     }
0141   }
0142
0143   return ProcessCode::SUCCESS;
0144 }
0145
0146 }  // namespace ActsExamples