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