src/EventData/ScalingCalibrator.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/EventData/ScalingCalibrator.hpp"
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Definitions/TrackParametrization.hpp"
0013 #include "Acts/EventData/MultiTrajectory.hpp"
0014 #include "Acts/EventData/SourceLink.hpp"
0015 #include "Acts/Geometry/GeometryContext.hpp"
0016 #include "Acts/Geometry/GeometryIdentifier.hpp"
0017 #include "Acts/Utilities/CalibrationContext.hpp"
0018 #include "ActsExamples/EventData/Cluster.hpp"
0019 #include "ActsExamples/EventData/IndexSourceLink.hpp"
0020 #include "ActsExamples/EventData/Measurement.hpp"
0021
0022 #include <algorithm>
0023 #include <array>
0024 #include <bitset>
0025 #include <cassert>
0026 #include <cstring>
0027 #include <filesystem>
0028 #include <map>
0029 #include <regex>
0030 #include <stdexcept>
0031 #include <string>
0032 #include <type_traits>
0033 #include <utility>
0034 #include <variant>
0035 #include <vector>
0036
0037 #include <TCollection.h>
0038 #include <TFile.h>
0039 #include <TH2.h>
0040 #include <TKey.h>
0041 #include <TList.h>
0042 #include <TString.h>
0043
0044 namespace Acts {
0045 class VectorMultiTrajectory;
0046 }  // namespace Acts
0047
0048 namespace detail {
0049
0050 std::pair<Acts::GeometryIdentifier, std::string> parseMapKey(
0051     const std::string& mapkey) {
0052   std::regex reg("^map_([0-9]+)-([0-9]+)-([0-9]+)_([xy]_.*)$");
0053   std::smatch matches;
0054
0055   if (!std::regex_search(mapkey, matches, reg) || matches.size() != 5) {
0056     throw std::runtime_error("Invalid map key: " + mapkey);
0057   }
0058
0059   std::size_t vol = std::stoull(matches[1].str());
0060   std::size_t lyr = std::stoull(matches[2].str());
0061   std::size_t mod = std::stoull(matches[3].str());
0062
0063   auto geoId =
0064       Acts::GeometryIdentifier().withVolume(vol).withLayer(lyr).withSensitive(
0065           mod);
0066
0067   std::string var(matches[4].str());
0068
0069   return {geoId, var};
0070 }
0071
0072 std::map<Acts::GeometryIdentifier, ActsExamples::ScalingCalibrator::MapTuple>
0073 readMaps(const std::filesystem::path& path) {
0074   std::map<Acts::GeometryIdentifier, ActsExamples::ScalingCalibrator::MapTuple>
0075       maps;
0076
0077   TFile ifile(path.c_str(), "READ");
0078   if (ifile.IsZombie()) {
0079     throw std::runtime_error("Unable to open TFile: " + path.string());
0080   }
0081
0082   TList* lst = ifile.GetListOfKeys();
0083   assert(lst != nullptr);
0084
0085   for (auto it = lst->begin(); it != lst->end(); ++it) {
0086     TKey* key = static_cast<TKey*>(*it);
0087     if (key != nullptr && std::strcmp(key->GetClassName(), "TH2D") == 0) {
0088       auto [geoId, var] = parseMapKey(key->GetName());
0089
0090       TH2D hist;
0091       key->Read(&hist);
0092
0093       if (var == "x_offset") {
0094         maps[geoId].x_offset = hist;
0095       } else if (var == "x_scale") {
0096         maps[geoId].x_scale = hist;
0097       } else if (var == "y_offset") {
0098         maps[geoId].y_offset = hist;
0099       } else if (var == "y_scale") {
0100         maps[geoId].y_scale = hist;
0101       } else {
0102         throw std::runtime_error("Unrecognized var: " + var);
0103       }
0104     }
0105   }
0106   return maps;
0107 }
0108
0109 std::bitset<3> readMask(const std::filesystem::path& path) {
0110   TFile ifile(path.c_str(), "READ");
0111   if (ifile.IsZombie()) {
0112     throw std::runtime_error("Unable to open TFile: " + path.string());
0113   }
0114
0115   TString* tstr = ifile.Get<TString>("v_mask");
0116   if (tstr == nullptr) {
0117     throw std::runtime_error("Unable to read mask");
0118   }
0119
0120   return std::bitset<3>(std::string{*tstr});
0121 }
0122
0123 }  // namespace detail
0124
0125 ActsExamples::ScalingCalibrator::ScalingCalibrator(
0126     const std::filesystem::path& path)
0127     : m_calib_maps{::detail::readMaps(path)},
0128       m_mask{::detail::readMask(path)} {}
0129
0130 void ActsExamples::ScalingCalibrator::calibrate(
0131     const MeasurementContainer& measurements, const ClusterContainer* clusters,
0132     const Acts::GeometryContext& /*gctx*/,
0133     const Acts::CalibrationContext& /*cctx*/,
0134     const Acts::SourceLink& sourceLink,
0135     Acts::VectorMultiTrajectory::TrackStateProxy& trackState) const {
0136   trackState.setUncalibratedSourceLink(Acts::SourceLink{sourceLink});
0137   const IndexSourceLink& idxSourceLink = sourceLink.get<IndexSourceLink>();
0138
0139   assert((idxSourceLink.index() < measurements.size()) &&
0140          "Source link index is outside the container bounds");
0141
0142   auto geoId = trackState.referenceSurface().geometryId();
0143   auto mgid =
0144       Acts::GeometryIdentifier()
0145           .withVolume(geoId.volume() *
0146                       static_cast<Acts::GeometryIdentifier::Value>(m_mask[2]))
0147           .withLayer(geoId.layer() *
0148                      static_cast<Acts::GeometryIdentifier::Value>(m_mask[1]))
0149           .withSensitive(
0150               geoId.sensitive() *
0151               static_cast<Acts::GeometryIdentifier::Value>(m_mask[0]));
0152   const Cluster& cl = clusters->at(idxSourceLink.index());
0153   ConstantTuple ct = m_calib_maps.at(mgid).at(cl.sizeLoc0, cl.sizeLoc1);
0154
0155   const ConstVariableBoundMeasurementProxy measurement =
0156       measurements.getMeasurement(idxSourceLink.index());
0157
0158   assert(measurement.contains(Acts::eBoundLoc0) &&
0159          "Measurement does not contain the required bound loc0");
0160   assert(measurement.contains(Acts::eBoundLoc1) &&
0161          "Measurement does not contain the required bound loc1");
0162
0163   auto boundLoc0 = measurement.indexOf(Acts::eBoundLoc0);
0164   auto boundLoc1 = measurement.indexOf(Acts::eBoundLoc1);
0165
0166   Acts::visit_measurement(measurement.size(), [&](auto N) -> void {
0167     constexpr std::size_t kMeasurementSize = decltype(N)::value;
0168     const ConstFixedBoundMeasurementProxy<kMeasurementSize> fixedMeasurement =
0169         static_cast<ConstFixedBoundMeasurementProxy<kMeasurementSize>>(
0170             measurement);
0171
0172     Acts::ActsVector<kMeasurementSize> calibratedParameters =
0173         fixedMeasurement.parameters();
0174     Acts::ActsSquareMatrix<kMeasurementSize> calibratedCovariance =
0175         fixedMeasurement.covariance();
0176
0177     calibratedParameters[boundLoc0] += ct.x_offset;
0178     calibratedParameters[boundLoc1] += ct.y_offset;
0179     calibratedCovariance(boundLoc0, boundLoc0) *= ct.x_scale;
0180     calibratedCovariance(boundLoc1, boundLoc1) *= ct.y_scale;
0181
0182     trackState.allocateCalibrated(calibratedParameters, calibratedCovariance);
0183     trackState.setProjectorSubspaceIndices(fixedMeasurement.subspaceIndices());
0184   });
0185 }