feat: Compute material variance in material mapping (#1318)

This PR is in preparation to the one that will add the automatic optimisation of the material mapping.

In this PR we add the option of computing the material variance of each bin. For this, the material mapping need to be run a second time with the material map used as an input. This is control by an option and will thus not affect the performance if not turned on. 

This PR also include a new material decorator that can be used as an input for the mapping, this decorator takes in input a map of ID and binning size and use it to create protomaterial for all the mapping surface (this will allow us to change the binning of the material without having to change the input json file).

Core/include/Acts/Material/AccumulatedMaterialSlab.hpp

@@ -48,6 +48,15 @@ class AccumulatedMaterialSlab {
   ///  in material structures.
   void accumulate(MaterialSlab slabAlongTrack, float pathCorrection = 1);
 
+  /// Use the accumulated material to update the material variance
+  ///
+  /// @param slabReference reference slab (from the map) used to compute the variance
+  /// @param useEmptyTrack indicate whether to consider an empty track store
+  ///
+  /// The material variance can be used to optimised the mapping process as it
+  /// should be inversly proportionnal to the map quality
+  void trackVariance(MaterialSlab slabReference, bool useEmptyTrack = false);
+
   /// Add the accumulated material for the current track to the total average.
   ///
   /// @param useEmptyTrack indicate whether to consider an empty track store
@@ -72,11 +81,23 @@ class AccumulatedMaterialSlab {
   /// the average thickness seen by the tracks.
   std::pair<MaterialSlab, unsigned int> totalAverage() const;
 
+  /// Return the material variance from all accumulated tracks.
+  ///
+  /// @returns Average material properties and the number of contributing tracks
+  ///
+  /// Only contains the information up to the last `.trackVariance(...)` call.
+  /// If there have been additional calls to `.accumulate(...)` afterwards, the
+  /// information is not part of the total average. The number of tracks is only
+  /// opdated on the call of `.trackAverage(...)`
+  std::pair<float, unsigned int> totalVariance() const;
+
  private:
   /// Averaged properties for a single track.
   MaterialSlab m_trackAverage;
   /// Averaged properties over multiple tracks.
   MaterialSlab m_totalAverage;
+  /// Averaged variance over multiple tracks.
+  float m_totalVariance = 0.0;
   // Number of tracks contributing to the total average.
   unsigned int m_totalCount = 0u;
 };

Core/include/Acts/Material/AccumulatedSurfaceMaterial.hpp

@@ -99,6 +99,23 @@ class AccumulatedSurfaceMaterial {
   std::array<size_t, 3> accumulate(const Vector3& gp, const MaterialSlab& mp,
                                    double pathCorrection = 1.);
 
+  /// Use the accumulated material to update the material variance
+  ///
+  /// @param trackBins The bins that were touched by this event
+  /// @param emptyHit indicator if this is an empty assignment
+  /// @param slabReference reference slab (from the map) used to compute the variance
+  /// If none is given, the average runs over all bins in the surface map
+  void trackVariance(const std::vector<std::array<size_t, 3>>& trackBins,
+                     MaterialSlab slabReference, bool emptyHit = false);
+
+  /// Use the accumulated material to update the material variance
+  ///
+  /// @param gp global position for the bin assignment
+  /// @param emptyHit indicator if this is an empty assignment
+  /// @param slabReference indicator if this is an empty assignment
+  void trackVariance(const Vector3& gp, MaterialSlab slabReference,
+                     bool emptyHit = false);
+
   /// Average the information accumulated from one mapped track
   ///
   /// @param trackBins The bins that were touched by this event

Core/include/Acts/Material/SurfaceMaterialMapper.hpp

@@ -90,6 +90,8 @@ class SurfaceMaterialMapper {
     bool emptyBinCorrection = true;
     /// Mapping output to debug stream
     bool mapperDebugOutput = false;
+    /// Compute the variance of each material slab (only if using an input map)
+    bool computeVariance = false;
   };
 
   /// @struct State
@@ -109,6 +111,10 @@ class SurfaceMaterialMapper {
     std::map<GeometryIdentifier, std::unique_ptr<const ISurfaceMaterial>>
         surfaceMaterial;
 
+    /// The surface material of the input tracking geometry
+    std::map<GeometryIdentifier, std::shared_ptr<const ISurfaceMaterial>>
+        inputSurfaceMaterial;
+
     /// The volume material of the input tracking geometry
     std::map<GeometryIdentifier, std::shared_ptr<const IVolumeMaterial>>
         volumeMaterial;

Core/src/Material/AccumulatedMaterialSlab.cpp

@@ -18,6 +18,24 @@ void Acts::AccumulatedMaterialSlab::accumulate(MaterialSlab slab,
   m_trackAverage = detail::combineSlabs(m_trackAverage, slab);
 }
 
+void Acts::AccumulatedMaterialSlab::trackVariance(MaterialSlab slabReference,
+                                                  bool useEmptyTrack) {
+  // Only use real tracks or if empty tracks are allowed.
+  if (useEmptyTrack or (0 < m_trackAverage.thickness())) {
+    float variance = ((1 / m_trackAverage.material().X0()) -
+                      (1 / slabReference.material().X0())) *
+                     ((1 / m_trackAverage.material().X0()) -
+                      (1 / slabReference.material().X0()));
+    if (m_totalCount == 0u) {
+      m_totalVariance = variance;
+    } else {
+      double weightTotal = m_totalCount / (m_totalCount + 1.0);
+      double weightTrack = 1 / (m_totalCount + 1.0);
+      m_totalVariance = weightTotal * m_totalVariance + weightTrack * variance;
+    }
+  }
+}
+
 void Acts::AccumulatedMaterialSlab::trackAverage(bool useEmptyTrack) {
   // average only real tracks or if empty tracks are allowed.
   if (useEmptyTrack or (0 < m_trackAverage.thickness())) {
@@ -43,3 +61,8 @@ std::pair<Acts::MaterialSlab, unsigned int>
 Acts::AccumulatedMaterialSlab::totalAverage() const {
   return {m_totalAverage, m_totalCount};
 }
+
+std::pair<float, unsigned int> Acts::AccumulatedMaterialSlab::totalVariance()
+    const {
+  return {m_totalVariance, m_totalCount};
+}

Core/src/Material/AccumulatedSurfaceMaterial.cpp

@@ -56,6 +56,43 @@ std::array<size_t, 3> Acts::AccumulatedSurfaceMaterial::accumulate(
   return bTriple;
 }
 
+// Void average for vacuum assignment
+void Acts::AccumulatedSurfaceMaterial::trackVariance(const Vector3& gp,
+                                                     MaterialSlab slabReference,
+                                                     bool emptyHit) {
+  if (m_binUtility.dimensions() == 0) {
+    m_accumulatedMaterial[0][0].trackVariance(slabReference, emptyHit);
+    return;
+  }
+  std::array<size_t, 3> bTriple = m_binUtility.binTriple(gp);
+  std::vector<std::array<size_t, 3>> trackBins = {bTriple};
+  trackVariance(trackBins, slabReference);
+}
+
+// Average the information accumulated during one event
+void Acts::AccumulatedSurfaceMaterial::trackVariance(
+    const std::vector<std::array<size_t, 3>>& trackBins,
+    MaterialSlab slabReference, bool emptyHit) {
+  // the homogeneous material case
+  if (m_binUtility.dimensions() == 0) {
+    m_accumulatedMaterial[0][0].trackVariance(slabReference, emptyHit);
+    return;
+  }
+  // The touched bins are known, so you can access them directly
+  if (not trackBins.empty()) {
+    for (auto bin : trackBins) {
+      m_accumulatedMaterial[bin[1]][bin[0]].trackVariance(slabReference);
+    }
+  } else {
+    // Touched bins are not known: Run over all bins
+    for (auto& matVec : m_accumulatedMaterial) {
+      for (auto& mat : matVec) {
+        mat.trackVariance(slabReference);
+      }
+    }
+  }
+}
+
 // Void average for vacuum assignment
 void Acts::AccumulatedSurfaceMaterial::trackAverage(const Vector3& gp,
                                                     bool emptyHit) {

Core/src/Material/SurfaceMaterialMapper.cpp

@@ -119,6 +119,10 @@ void Acts::SurfaceMaterialMapper::checkAndInsert(State& mState,
   auto surfaceMaterial = surface.surfaceMaterial();
   // Check if the surface has a proxy
   if (surfaceMaterial != nullptr) {
+    if (m_cfg.computeVariance) {
+      mState.inputSurfaceMaterial[surface.geometryId()] =
+          surface.surfaceMaterialSharedPtr();
+    }
     auto geoID = surface.geometryId();
     size_t volumeID = geoID.volume();
     ACTS_DEBUG("Material surface found with volumeID " << volumeID);
@@ -257,9 +261,12 @@ void Acts::SurfaceMaterialMapper::mapMaterialTrack(
 
   // To remember the bins of this event
   using MapBin = std::pair<AccumulatedSurfaceMaterial*, std::array<size_t, 3>>;
+  using MaterialBin = std::pair<AccumulatedSurfaceMaterial*,
+                                std::shared_ptr<const ISurfaceMaterial>>;
   std::multimap<AccumulatedSurfaceMaterial*, std::array<size_t, 3>>
       touchedMapBins;
-
+  std::map<AccumulatedSurfaceMaterial*, std::shared_ptr<const ISurfaceMaterial>>
+      touchedMaterialBin;
   if (sfIter != mappingSurfaces.end() &&
       sfIter->surface->surfaceMaterial()->mappingType() ==
           Acts::MappingType::PostMapping) {
@@ -365,7 +372,11 @@ void Acts::SurfaceMaterialMapper::mapMaterialTrack(
     // Now assign the material for the accumulation process
     auto tBin = currentAccMaterial->second.accumulate(
         currentPos, rmIter->materialSlab, currentPathCorrection);
-    touchedMapBins.insert(MapBin(&(currentAccMaterial->second), tBin));
+    auto itBin =
+        touchedMapBins.insert(MapBin(&(currentAccMaterial->second), tBin));
+    if (m_cfg.computeVariance) {
+      touchedMaterialBin[itBin->first] = mState.inputSurfaceMaterial[currentID];
+    }
     ++assignedMaterial[currentID];
     // Update the material interaction with the associated surface
     rmIter->surface = sfIter->surface;
@@ -381,6 +392,11 @@ void Acts::SurfaceMaterialMapper::mapMaterialTrack(
   // After mapping this track, average the touched bins
   for (auto tmapBin : touchedMapBins) {
     std::vector<std::array<size_t, 3>> trackBins = {tmapBin.second};
+    if (m_cfg.computeVariance) {
+      tmapBin.first->trackVariance(
+          trackBins, touchedMaterialBin[tmapBin.first]->materialSlab(
+                         trackBins[0][0], trackBins[0][1]));
+    }
     tmapBin.first->trackAverage(trackBins);
   }
 
@@ -395,6 +411,12 @@ void Acts::SurfaceMaterialMapper::mapMaterialTrack(
       // list of assigned surfaces
       if (assignedMaterial[mgID] == 0) {
         auto missedMaterial = mState.accumulatedMaterial.find(mgID);
+        if (m_cfg.computeVariance) {
+          missedMaterial->second.trackVariance(
+              mSurface.position,
+              mState.inputSurfaceMaterial[currentID]->materialSlab(
+                  mSurface.position));
+        }
         missedMaterial->second.trackAverage(mSurface.position, true);
       }
     }