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   Acts::GeometryIdentifier geoId;
0064   geoId.setVolume(vol);
0065   geoId.setLayer(lyr);
0066   geoId.setSensitive(mod);
0067
0068   std::string var(matches[4].str());
0069
0070   return {geoId, var};
0071 }
0072
0073 std::map<Acts::GeometryIdentifier, ActsExamples::ScalingCalibrator::MapTuple>
0074 readMaps(const std::filesystem::path& path) {
0075   std::map<Acts::GeometryIdentifier, ActsExamples::ScalingCalibrator::MapTuple>
0076       maps;
0077
0078   TFile ifile(path.c_str(), "READ");
0079   if (ifile.IsZombie()) {
0080     throw std::runtime_error("Unable to open TFile: " + path.string());
0081   }
0082
0083   TList* lst = ifile.GetListOfKeys();
0084   assert(lst != nullptr);
0085
0086   for (auto it = lst->begin(); it != lst->end(); ++it) {
0087     TKey* key = static_cast<TKey*>(*it);
0088     if (key != nullptr && std::strcmp(key->GetClassName(), "TH2D") == 0) {
0089       auto [geoId, var] = parseMapKey(key->GetName());
0090
0091       TH2D hist;
0092       key->Read(&hist);
0093
0094       if (var == "x_offset") {
0095         maps[geoId].x_offset = hist;
0096       } else if (var == "x_scale") {
0097         maps[geoId].x_scale = hist;
0098       } else if (var == "y_offset") {
0099         maps[geoId].y_offset = hist;
0100       } else if (var == "y_scale") {
0101         maps[geoId].y_scale = hist;
0102       } else {
0103         throw std::runtime_error("Unrecognized var: " + var);
0104       }
0105     }
0106   }
0107   return maps;
0108 }
0109
0110 std::bitset<3> readMask(const std::filesystem::path& path) {
0111   TFile ifile(path.c_str(), "READ");
0112   if (ifile.IsZombie()) {
0113     throw std::runtime_error("Unable to open TFile: " + path.string());
0114   }
0115
0116   TString* tstr = ifile.Get<TString>("v_mask");
0117   if (tstr == nullptr) {
0118     throw std::runtime_error("Unable to read mask");
0119   }
0120
0121   return std::bitset<3>(std::string{*tstr});
0122 }
0123
0124 }  // namespace detail
0125
0126 ActsExamples::ScalingCalibrator::ScalingCalibrator(
0127     const std::filesystem::path& path)
0128     : m_calib_maps{::detail::readMaps(path)},
0129       m_mask{::detail::readMask(path)} {}
0130
0131 void ActsExamples::ScalingCalibrator::calibrate(
0132     const MeasurementContainer& measurements, const ClusterContainer* clusters,
0133     const Acts::GeometryContext& /*gctx*/,
0134     const Acts::CalibrationContext& /*cctx*/,
0135     const Acts::SourceLink& sourceLink,
0136     Acts::VectorMultiTrajectory::TrackStateProxy& trackState) const {
0137   trackState.setUncalibratedSourceLink(Acts::SourceLink{sourceLink});
0138   const IndexSourceLink& idxSourceLink = sourceLink.get<IndexSourceLink>();
0139
0140   assert((idxSourceLink.index() < measurements.size()) &&
0141          "Source link index is outside the container bounds");
0142
0143   auto geoId = trackState.referenceSurface().geometryId();
0144   Acts::GeometryIdentifier mgid;
0145   mgid.setVolume(geoId.volume() *
0146                  static_cast<Acts::GeometryIdentifier::Value>(m_mask[2]));
0147   mgid.setLayer(geoId.layer() *
0148                 static_cast<Acts::GeometryIdentifier::Value>(m_mask[1]));
0149   mgid.setSensitive(geoId.sensitive() *
0150                     static_cast<Acts::GeometryIdentifier::Value>(m_mask[0]));
0151   const Cluster& cl = clusters->at(idxSourceLink.index());
0152   ConstantTuple ct = m_calib_maps.at(mgid).at(cl.sizeLoc0, cl.sizeLoc1);
0153
0154   const ConstVariableBoundMeasurementProxy measurement =
0155       measurements.getMeasurement(idxSourceLink.index());
0156
0157   assert(measurement.contains(Acts::eBoundLoc0) &&
0158          "Measurement does not contain the required bound loc0");
0159   assert(measurement.contains(Acts::eBoundLoc1) &&
0160          "Measurement does not contain the required bound loc1");
0161
0162   auto boundLoc0 = measurement.indexOf(Acts::eBoundLoc0);
0163   auto boundLoc1 = measurement.indexOf(Acts::eBoundLoc1);
0164
0165   Acts::visit_measurement(measurement.size(), [&](auto N) -> void {
0166     constexpr std::size_t kMeasurementSize = decltype(N)::value;
0167     const ConstFixedBoundMeasurementProxy<kMeasurementSize> fixedMeasurement =
0168         static_cast<ConstFixedBoundMeasurementProxy<kMeasurementSize>>(
0169             measurement);
0170
0171     Acts::ActsVector<kMeasurementSize> calibratedParameters =
0172         fixedMeasurement.parameters();
0173     Acts::ActsSquareMatrix<kMeasurementSize> calibratedCovariance =
0174         fixedMeasurement.covariance();
0175
0176     calibratedParameters[boundLoc0] += ct.x_offset;
0177     calibratedParameters[boundLoc1] += ct.y_offset;
0178     calibratedCovariance(boundLoc0, boundLoc0) *= ct.x_scale;
0179     calibratedCovariance(boundLoc1, boundLoc1) *= ct.y_scale;
0180
0181     trackState.allocateCalibrated(calibratedParameters, calibratedCovariance);
0182     trackState.setProjectorSubspaceIndices(fixedMeasurement.subspaceIndices());
0183   });
0184 }