feat: Simple scaling calibrator (#2085)

This PR adds a new measurement calibrator, which offsets the measurement positions and scales the error based on a set of 2-D maps of cluster sizes (loc0 / loc1). The granularity of the maps can be different based on Volume, Layer, or Sensitive surface.

This PR does not add documentation of the method -- that will follow later once I fully develop a good setup to derive some good maps (I have a WIP implementation that I used for testing).

Lastly, note that this is not meant to be a state of the art method, it's basically to showcase the simplest non-trivial use of the measurement calibration interface. 

Ping @benjaminhuth if you can look at it, or suggest someone else!

Core/include/Acts/EventData/Measurement.hpp

@@ -106,7 +106,7 @@ class Measurement {
   /// The measurement indices
   constexpr std::array<indices_t, kSize> indices() const {
     std::array<uint8_t, kSize> subInds = m_subspace.indices();
-    std::array<indices_t, kSize> inds;
+    std::array<indices_t, kSize> inds{};
     for (size_t i = 0; i < kSize; i++) {
       inds[i] = static_cast<indices_t>(subInds[i]);
     }

Examples/Framework/CMakeLists.txt

@@ -3,6 +3,7 @@ include(ActsTargetLinkLibrariesSystem)
 add_library(
   ActsExamplesFramework SHARED
   src/EventData/MeasurementCalibration.cpp
+  src/EventData/ScalingCalibrator.cpp
   src/Framework/IAlgorithm.cpp
   src/Framework/SequenceElement.cpp
   src/Framework/WhiteBoard.cpp

Examples/Framework/include/ActsExamples/EventData/MeasurementCalibration.hpp

@@ -29,7 +29,7 @@ class MeasurementCalibrator {
           trackState) const = 0;
 
   virtual ~MeasurementCalibrator() = default;
-  virtual bool needsClusters() { return false; }
+  virtual bool needsClusters() const { return false; }
 };
 
 // Calibrator to convert an index source link to a measurement as-is

Examples/Framework/include/ActsExamples/EventData/ScalingCalibrator.hpp

@@ -0,0 +1,65 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2023 CERN for the benefit of the Acts project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#pragma once
+
+#include <Acts/Geometry/GeometryIdentifier.hpp>
+#include <ActsExamples/EventData/MeasurementCalibration.hpp>
+
+#include <filesystem>
+
+#include <TFile.h>
+#include <TH2D.h>
+
+namespace ActsExamples {
+
+class ScalingCalibrator : public MeasurementCalibrator {
+ public:
+  struct ConstantTuple {
+    double x_offset{0};
+    double x_scale{1};
+    double y_offset{0};
+    double y_scale{1};
+  };
+
+  struct MapTuple {
+    TH2D x_offset;
+    TH2D x_scale;
+    TH2D y_offset;
+    TH2D y_scale;
+
+    ConstantTuple at(size_t sizeLoc0, size_t sizeLoc1) const {
+      ConstantTuple ct;
+      ct.x_offset =
+          x_offset.GetBinContent(x_offset.FindFixBin(sizeLoc0, sizeLoc1));
+      ct.x_scale =
+          x_scale.GetBinContent(x_scale.FindFixBin(sizeLoc0, sizeLoc1));
+      ct.y_offset =
+          y_offset.GetBinContent(y_offset.FindFixBin(sizeLoc0, sizeLoc1));
+      ct.y_scale =
+          y_scale.GetBinContent(y_scale.FindFixBin(sizeLoc0, sizeLoc1));
+      return ct;
+    }
+  };
+
+  ScalingCalibrator(const std::filesystem::path& path);
+
+  void calibrate(
+      const MeasurementContainer& measurements,
+      const ClusterContainer* clusters, const Acts::GeometryContext& /*gctx*/,
+      Acts::MultiTrajectory<Acts::VectorMultiTrajectory>::TrackStateProxy&
+          trackState) const override;
+
+  bool needsClusters() const override { return true; }
+
+ private:
+  std::map<Acts::GeometryIdentifier, MapTuple> m_calib_maps;
+  std::bitset<3> m_mask;
+};
+
+}  // namespace ActsExamples

Examples/Framework/src/EventData/MeasurementCalibration.cpp

@@ -13,8 +13,9 @@ void ActsExamples::PassThroughCalibrator::calibrate(
     const ClusterContainer* /*clusters*/, const Acts::GeometryContext& /*gctx*/,
     Acts::MultiTrajectory<Acts::VectorMultiTrajectory>::TrackStateProxy&
         trackState) const {
-  const IndexSourceLink& sourceLink =
-      trackState.getUncalibratedSourceLink().get<IndexSourceLink>();
+  Acts::SourceLink usl = trackState.getUncalibratedSourceLink();
+  const IndexSourceLink& sourceLink = usl.get<IndexSourceLink>();
+
   assert((sourceLink.index() < measurements.size()) and
          "Source link index is outside the container bounds");
 

Examples/Framework/src/EventData/ScalingCalibrator.cpp

@@ -0,0 +1,150 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2023 CERN for the benefit of the Acts project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#include <Acts/Definitions/TrackParametrization.hpp>
+#include <ActsExamples/EventData/ScalingCalibrator.hpp>
+
+#include <regex>
+
+#include <TKey.h>
+
+namespace detail {
+
+std::pair<Acts::GeometryIdentifier, std::string> parseMapKey(
+    const std::string& mapkey) {
+  std::regex reg("^map_([0-9]+)-([0-9]+)-([0-9]+)_([xy]_.*)$");
+  std::smatch matches;
+
+  if (std::regex_search(mapkey, matches, reg) && matches.size() == 5) {
+    size_t vol = std::stoull(matches[1].str());
+    size_t lyr = std::stoull(matches[2].str());
+    size_t mod = std::stoull(matches[3].str());
+
+    Acts::GeometryIdentifier geoId;
+    geoId.setVolume(vol);
+    geoId.setLayer(lyr);
+    geoId.setSensitive(mod);
+
+    std::string var(matches[4].str());
+
+    return std::make_pair(geoId, var);
+  } else {
+    throw std::runtime_error("Invalid map key: " + mapkey);
+  }
+}
+
+std::map<Acts::GeometryIdentifier, ActsExamples::ScalingCalibrator::MapTuple>
+readMaps(const std::filesystem::path& path) {
+  std::map<Acts::GeometryIdentifier, ActsExamples::ScalingCalibrator::MapTuple>
+      maps;
+
+  TFile ifile(path.c_str(), "READ");
+  if (ifile.IsZombie()) {
+    throw std::runtime_error("Unable to open TFile: " + path.string());
+  }
+
+  TList* lst = ifile.GetListOfKeys();
+  assert(lst != nullptr);
+
+  for (auto it = lst->begin(); it != lst->end(); ++it) {
+    TKey* key = static_cast<TKey*>(*it);
+    if (std::strcmp(key->GetClassName(), "TH2D") == 0) {
+      auto [geoId, var] = parseMapKey(key->GetName());
+
+      TH2D hist;
+      key->Read(&hist);
+
+      if (var == "x_offset") {
+        maps[geoId].x_offset = hist;
+      } else if (var == "x_scale") {
+        maps[geoId].x_scale = hist;
+      } else if (var == "y_offset") {
+        maps[geoId].y_offset = hist;
+      } else if (var == "y_scale") {
+        maps[geoId].y_scale = hist;
+      } else {
+        throw std::runtime_error("Unrecognized var: " + var);
+      }
+    }
+  }
+  return maps;
+}
+
+std::bitset<3> readMask(const std::filesystem::path& path) {
+  TFile ifile(path.c_str(), "READ");
+  if (ifile.IsZombie()) {
+    throw std::runtime_error("Unable to open TFile: " + path.string());
+  }
+
+  TString* tstr = ifile.Get<TString>("v_mask");
+  if (tstr == nullptr) {
+    throw std::runtime_error("Unable to read mask");
+  }
+
+  return std::bitset<3>(std::string(*tstr));
+}
+
+}  // namespace detail
+
+ActsExamples::ScalingCalibrator::ScalingCalibrator(
+    const std::filesystem::path& path)
+    : m_calib_maps{::detail::readMaps(path)},
+      m_mask{::detail::readMask(path)} {}
+
+void ActsExamples::ScalingCalibrator::calibrate(
+    const MeasurementContainer& measurements, const ClusterContainer* clusters,
+    const Acts::GeometryContext& /*gctx*/,
+    Acts::MultiTrajectory<Acts::VectorMultiTrajectory>::TrackStateProxy&
+        trackState) const {
+  Acts::SourceLink usl = trackState.getUncalibratedSourceLink();
+  const IndexSourceLink& sourceLink = usl.get<IndexSourceLink>();
+
+  assert((sourceLink.index() < measurements.size()) and
+         "Source link index is outside the container bounds");
+
+  Acts::GeometryIdentifier mgid;
+  mgid.setVolume(sourceLink.geometryId().volume() *
+                 static_cast<Acts::GeometryIdentifier::Value>(m_mask[2]));
+  mgid.setLayer(sourceLink.geometryId().layer() *
+                static_cast<Acts::GeometryIdentifier::Value>(m_mask[1]));
+  mgid.setSensitive(sourceLink.geometryId().sensitive() *
+                    static_cast<Acts::GeometryIdentifier::Value>(m_mask[0]));
+  const Cluster& cl = clusters->at(sourceLink.index());
+  ConstantTuple ct = m_calib_maps.at(mgid).at(cl.sizeLoc0, cl.sizeLoc1);
+
+  std::visit(
+      [&](const auto& meas) {
+        auto E = meas.expander();
+        auto P = meas.projector();
+
+        Acts::ActsVector<Acts::eBoundSize> fpar = E * meas.parameters();
+
+        Acts::ActsSymMatrix<Acts::eBoundSize> fcov =
+            E * meas.covariance() * E.transpose();
+
+        fpar[Acts::eBoundLoc0] += ct.x_offset;
+        fpar[Acts::eBoundLoc1] += ct.y_offset;
+        fcov(Acts::eBoundLoc0, Acts::eBoundLoc0) *= ct.x_scale;
+        fcov(Acts::eBoundLoc1, Acts::eBoundLoc1) *= ct.y_scale;
+
+        constexpr size_t kSize =
+            std::remove_reference_t<decltype(meas)>::size();
+        std::array<Acts::BoundIndices, kSize> indices = meas.indices();
+        Acts::ActsVector<kSize> cpar = P * fpar;
+        Acts::ActsSymMatrix<kSize> ccov = P * fcov * P.transpose();
+
+        Acts::SourceLink sl{sourceLink.geometryId(), sourceLink};
+
+        Acts::Measurement<Acts::BoundIndices, kSize> cmeas(std::move(sl),
+                                                           indices, cpar, ccov);
+
+        trackState.allocateCalibrated(cmeas.size());
+        trackState.setCalibrated(cmeas);
+      },
+      (measurements)[sourceLink.index()]);
+}

Examples/Python/python/acts/examples/reconstruction.py

@@ -819,6 +819,7 @@ def addKalmanTracks(
     inputProtoTracks: str = "truth_particle_tracks",
     multipleScattering: bool = True,
     energyLoss: bool = True,
+    calibrationConfigFile: str = None,
     clusters: str = None,
     logLevel: Optional[acts.logging.Level] = None,
 ) -> None:
@@ -843,6 +844,11 @@ def addKalmanTracks(
         "level": customLogLevel(),
     }
 
+    if calibrationConfigFile is None:
+        calibrator = acts.examples.makePassThroughCalibrator()
+    else:
+        calibrator = acts.examples.makeScalingCalibrator(calibrationConfigFile)
+
     fitAlg = acts.examples.TrackFittingAlgorithm(
         level=customLogLevel(),
         inputMeasurements="measurements",
@@ -855,7 +861,7 @@ def addKalmanTracks(
         fit=acts.examples.makeKalmanFitterFunction(
             trackingGeometry, field, **kalmanOptions
         ),
-        calibrator=acts.examples.makePassThroughCalibrator(),
+        calibrator=calibrator,
     )
     s.addAlgorithm(fitAlg)
     s.addWhiteboardAlias("tracks", fitAlg.config.outputTracks)

Examples/Python/src/TrackFitting.cpp

@@ -9,6 +9,7 @@
 #include "Acts/MagneticField/MagneticFieldProvider.hpp"
 #include "Acts/Plugins/Python/Utilities.hpp"
 #include "Acts/Utilities/Logger.hpp"
+#include "ActsExamples/EventData/ScalingCalibrator.hpp"
 #include "ActsExamples/TrackFitting/RefittingAlgorithm.hpp"
 #include "ActsExamples/TrackFitting/SurfaceSortingAlgorithm.hpp"
 #include "ActsExamples/TrackFitting/TrackFitterFunction.hpp"
@@ -75,6 +76,13 @@ void addTrackFitting(Context& ctx) {
               return std::make_shared<PassThroughCalibrator>();
             });
 
+    mex.def(
+        "makeScalingCalibrator",
+        [](const char* path) -> std::shared_ptr<MeasurementCalibrator> {
+          return std::make_shared<ActsExamples::ScalingCalibrator>(path);
+        },
+        py::arg("path"));
+
     py::enum_<Acts::MixtureReductionMethod>(mex, "FinalReductionMethod")
         .value("mean", Acts::MixtureReductionMethod::eMean)
         .value("maxWeight", Acts::MixtureReductionMethod::eMaxWeight);