feat: generate measurements with muon space points

Examples/Algorithms/Digitization/include/ActsExamples/Digitization/MuonSpacePointDigitizer.hpp

@@ -11,6 +11,7 @@
 #include "Acts/Geometry/TrackingGeometry.hpp"
 #include "Acts/Utilities/Logger.hpp"
 #include "Acts/Utilities/RangeXD.hpp"
+#include "ActsExamples/EventData/Measurement.hpp"
 #include "ActsExamples/EventData/MuonSpacePoint.hpp"
 #include "ActsExamples/EventData/MuonSpacePointCalibrator.hpp"
 #include "ActsExamples/EventData/SimHit.hpp"
@@ -32,6 +33,12 @@ class MuonSpacePointDigitizer final : public IAlgorithm {
     std::string inputParticles{"particles_simulated"};
     /// @brief Name of the output space points collection
     std::string outputSpacePoints{"MuonSpacePoints"};
+    /// @brief Output measurements collection.
+    std::string outputMeasurements = "measurements";
+    /// @brief Output collection to map measured hits to simulated hits.
+    std::string outputMeasurementSimHitsMap = "measurement_simhits_map";
+    /// @brief Output collection to map particles to simulated hits.
+    std::string outputSimHitMeasurementsMap = "simhit_measurements_map";
     /// @brief Random number generator service
     std::shared_ptr<const RandomNumbers> randomNumbers{};
     /// @brief Pointer to the muon calibrator to fetch the smearing constants
@@ -99,5 +106,11 @@ class MuonSpacePointDigitizer final : public IAlgorithm {
   /// @brief Data handle for the output space points
   WriteDataHandle<MuonSpacePointContainer> m_outputSpacePoints{this,
                                                                "SpacePoints"};
+  WriteDataHandle<MeasurementContainer> m_outputMeasurements{
+      this, "OutputMeasurements"};
+  WriteDataHandle<IndexMultimap<Index>> m_outputMeasurementSimHitsMap{
+      this, "OutputMeasurementSimHitsMap"};
+  WriteDataHandle<InverseMultimap<Index>> m_outputSimHitMeasurementsMap{
+      this, "OutputSimHitMeasurementsMap"};
 };
 }  // namespace ActsExamples

Examples/Algorithms/Digitization/src/MuonSpacePointDigitizer.cpp

@@ -20,6 +20,7 @@
 #include "Acts/Utilities/Helpers.hpp"
 #include "Acts/Utilities/MathHelpers.hpp"
 #include "Acts/Utilities/StringHelpers.hpp"
+#include "ActsExamples/Digitization/MeasurementCreation.hpp"
 #include "ActsExamples/Digitization/Smearers.hpp"
 #include "ActsExamples/EventData/MuonSpacePoint.hpp"
 
@@ -110,12 +111,26 @@ MuonSpacePointDigitizer::MuonSpacePointDigitizer(const Config& cfg,
   if (m_cfg.outputSpacePoints.empty()) {
     throw std::invalid_argument("No output space points were defined");
   }
+  if (m_cfg.outputMeasurements.empty()) {
+    throw std::invalid_argument("Missing measurements output collection");
+  }
+  if (m_cfg.outputMeasurementSimHitsMap.empty()) {
+    throw std::invalid_argument(
+        "Missing hit-to-simulated-hits map output collection");
+  }
+  if (m_cfg.outputSimHitMeasurementsMap.empty()) {
+    throw std::invalid_argument(
+        "Missing particle-to-simulated-hits map output collection");
+  }
   ACTS_DEBUG("Retrieve sim hits and particles from "
              << m_cfg.inputSimHits << " & " << m_cfg.inputParticles);
   ACTS_DEBUG("Write produced space points to " << m_cfg.outputSpacePoints);
   m_inputSimHits.initialize(m_cfg.inputSimHits);
   m_inputParticles.initialize(m_cfg.inputParticles);
   m_outputSpacePoints.initialize(m_cfg.outputSpacePoints);
+  m_outputMeasurements.initialize(m_cfg.outputMeasurements);
+  m_outputMeasurementSimHitsMap.initialize(m_cfg.outputMeasurementSimHitsMap);
+  m_outputSimHitMeasurementsMap.initialize(m_cfg.outputSimHitMeasurementsMap);
 }
 
 ProcessCode MuonSpacePointDigitizer::initialize() {
@@ -174,6 +189,13 @@ ProcessCode MuonSpacePointDigitizer::execute(
 
   GeometryContext gctx{};
 
+  // Prepare output containers
+  // need list here for stable addresses
+  MeasurementContainer measurements;
+  IndexMultimap<Index> measurementSimHitsMap;
+  measurements.reserve(gotSimHits.size());
+  measurementSimHitsMap.reserve(gotSimHits.size());
+
   using MuonId_t = MuonSpacePoint::MuonId;
   auto rndEngine = m_cfg.randomNumbers->spawnGenerator(ctx);
   /// temporary output container to group the hits per chamber volume
@@ -203,6 +225,7 @@ ProcessCode MuonSpacePointDigitizer::execute(
     // Iterate over all simHits in a single module
     for (auto h = moduleSimHits.begin(); h != moduleSimHits.end(); ++h) {
       const auto& simHit = *h;
+      const auto simHitIdx = gotSimHits.index_of(h);
 
       // Convert the hit trajectory into local coordinates
       const Vector3 locPos = surfLocToGlob.inverse() * simHit.position();
@@ -281,6 +304,24 @@ ProcessCode MuonSpacePointDigitizer::execute(
                   calibCfg.rpcPhiStripPitch, calibCfg.rpcEtaStripPitch,
                   m_cfg.digitizeTime ? calibCfg.rpcTimeResolution : 0.);
 
+              // Set measurement
+              DigitizedParameters dParameters;
+
+              auto cov = newSp.covariance();
+
+              dParameters.indices.push_back(Acts::eBoundLoc0);
+              dParameters.values.push_back(smearedHit[ePos0]);
+              dParameters.variances.push_back(cov[0]);
+
+              dParameters.indices.push_back(Acts::eBoundLoc1);
+              dParameters.values.push_back(smearedHit[ePos1]);
+              dParameters.variances.push_back(cov[1]);
+
+              auto measurement =
+                  createMeasurement(measurements, moduleGeoId, dParameters);
+              measurementSimHitsMap.emplace_hint(
+                  measurementSimHitsMap.end(), measurement.index(), simHitIdx);
+
               break;
             }
             /// Endcap strips not yet available
@@ -321,15 +362,16 @@ ProcessCode MuonSpacePointDigitizer::execute(
             convertSp = false;
             break;
           }
-          double driftR =
+          double smearedDriftR =
               (*Digitization::Gauss{uncert}(unsmearedR, rndEngine)).first;
 
           // bounds
           const auto& lBounds = static_cast<const LineBounds&>(bounds);
           const double maxR = lBounds.get(LineBounds::eR);
           const double maxZ = lBounds.get(LineBounds::eHalfLengthZ);
           /// The generated hit is unphysical
-          if (driftR < 0. || driftR > maxR || std::abs(nominalPos.z()) > maxZ) {
+          if (smearedDriftR < 0. || smearedDriftR > maxR ||
+              std::abs(nominalPos.z()) > maxZ) {
             convertSp = false;
             break;
           }
@@ -346,7 +388,10 @@ ProcessCode MuonSpacePointDigitizer::execute(
 
           const double sigmaZ = 0.5 * maxZ;
 
-          newSp.setRadius(driftR);
+          const double smearedZ =
+              (*Digitization::Gauss{sigmaZ}(nominalPos.z(), rndEngine)).first;
+
+          newSp.setRadius(smearedDriftR);
           newSp.setCovariance(square(uncert), square(sigmaZ), 0.);
 
           newSp.defineCoordinates(
@@ -362,6 +407,23 @@ ProcessCode MuonSpacePointDigitizer::execute(
           id.setCoordFlags(true, false);
           newSp.setId(id);
 
+          // Set measurement
+          DigitizedParameters dParameters;
+
+          auto cov = newSp.covariance();
+
+          dParameters.indices.push_back(Acts::eBoundLoc0);
+          dParameters.values.push_back(smearedDriftR);
+          dParameters.variances.push_back(cov[0]);
+
+          dParameters.indices.push_back(Acts::eBoundLoc1);
+          dParameters.values.push_back(smearedZ);
+          dParameters.variances.push_back(cov[1]);
+
+          auto measurement =
+              createMeasurement(measurements, moduleGeoId, dParameters);
+          measurementSimHitsMap.emplace_hint(measurementSimHitsMap.end(),
+                                             measurement.index(), simHitIdx);
           break;
         }
         default:
@@ -396,6 +458,10 @@ ProcessCode MuonSpacePointDigitizer::execute(
   }
 
   m_outputSpacePoints(ctx, std::move(outSpacePoints));
+  m_outputMeasurements(ctx, std::move(measurements));
+  m_outputSimHitMeasurementsMap(ctx,
+                                invertIndexMultimap(measurementSimHitsMap));
+  m_outputMeasurementSimHitsMap(ctx, std::move(measurementSimHitsMap));
 
   return ProcessCode::SUCCESS;
 }