feat: Add the option of bin counting histograms for the surfaces material maps

Core/include/Acts/Material/AccumulatedSurfaceMaterial.hpp

@@ -142,8 +142,10 @@ class AccumulatedSurfaceMaterial {
   void trackAverage(const Vector3& gp, bool emptyHit = false);
 
   /// Total average creates SurfaceMaterial
-  /// @return Unique pointer to the averaged surface material
-  std::unique_ptr<const ISurfaceMaterial> totalAverage();
+  /// @param recordBinCounts flag to record the counts of the mapped bins from all the tracks
+  /// @return Unique pointer to the surface material with the bin counts information if enabled
+  std::unique_ptr<const ISurfaceMaterial> totalAverage(
+      bool recordBinCounts = false);
 
   /// Access to the accumulated material
   /// @return Reference to the matrix of accumulated material data

Core/include/Acts/Material/BinnedSurfaceMaterial.hpp

@@ -34,10 +34,12 @@ class BinnedSurfaceMaterial : public ISurfaceMaterial {
   /// @param binUtility defines the binning structure on the surface (copied)
   /// @param materialVector is the vector of material slabs as recorded (moved)
   /// @param splitFactor is the pre/post splitting directive
+  /// @param binCounts is the vector of the bin counts recorded
   /// @param mappingType is the type of surface mapping associated to the surface
   BinnedSurfaceMaterial(const BinUtility& binUtility,
                         MaterialSlabVector materialVector,
                         double splitFactor = 0.,
+                        std::vector<unsigned int> binCounts = {},
                         MappingType mappingType = MappingType::Default);
 
   /// Explicit constructor with only full MaterialSlab,
@@ -46,10 +48,12 @@ class BinnedSurfaceMaterial : public ISurfaceMaterial {
   /// @param binUtility defines the binning structure on the surface (copied)
   /// @param materialMatrix is the matrix of material slabs as recorded (moved)
   /// @param splitFactor is the pre/post splitting directive
+  /// @param binCounts is the matrix of the bin counts recorded
   /// @param mappingType is the type of surface mapping associated to the surface
   BinnedSurfaceMaterial(const BinUtility& binUtility,
                         MaterialSlabMatrix materialMatrix,
                         double splitFactor = 0.,
+                        std::vector<std::vector<unsigned int>> binCounts = {},
                         MappingType mappingType = MappingType::Default);
 
   /// Scale operation
@@ -66,6 +70,12 @@ class BinnedSurfaceMaterial : public ISurfaceMaterial {
   /// @return Reference to the complete matrix of material slabs
   const MaterialSlabMatrix& fullMaterial() const { return m_fullMaterial; }
 
+  /// @brief Retrieve the bin counts matrix of the mapped bins
+  /// @return Reference to the bin counts matrix of the map
+  const std::vector<std::vector<unsigned int>>& binCounts() const {
+    return m_binCounts;
+  }
+
   /// @copydoc ISurfaceMaterial::materialSlab(const Vector2&) const
   const MaterialSlab& materialSlab(const Vector2& lp) const final;
 
@@ -91,6 +101,9 @@ class BinnedSurfaceMaterial : public ISurfaceMaterial {
 
   /// The five different MaterialSlab
   MaterialSlabMatrix m_fullMaterial;
+
+  // The bin counts of the mapped bins
+  std::vector<std::vector<unsigned int>> m_binCounts;
 };
 
 }  // namespace Acts

Core/include/Acts/Material/BinnedSurfaceMaterialAccumulator.hpp

@@ -28,6 +28,9 @@ class BinnedSurfaceMaterialAccumulator final
     /// Correct for empty bins (recommended)
     bool emptyBinCorrection = true;
 
+    // Record the counts in the mapped bins (useful for statistics studies )
+    bool recordBinCounts = false;
+
     /// The surfaces to be used for the accumulation
     std::vector<const Surface*> materialSurfaces = {};
   };

Core/src/Material/AccumulatedSurfaceMaterial.cpp

@@ -132,7 +132,7 @@ void Acts::AccumulatedSurfaceMaterial::trackAverage(
 
 /// Total average creates SurfaceMaterial
 std::unique_ptr<const Acts::ISurfaceMaterial>
-Acts::AccumulatedSurfaceMaterial::totalAverage() {
+Acts::AccumulatedSurfaceMaterial::totalAverage(bool recordBinCounts) {
   if (m_binUtility.bins() == 1) {
     // Return HomogeneousSurfaceMaterial
     return std::make_unique<HomogeneousSurfaceMaterial>(
@@ -142,13 +142,25 @@ Acts::AccumulatedSurfaceMaterial::totalAverage() {
   MaterialSlabMatrix mpMatrix(
       m_binUtility.bins(1),
       MaterialSlabVector(m_binUtility.bins(0), MaterialSlab::Nothing()));
+
+  std::vector<std::vector<unsigned int>> binCounts{};
+
+  if (recordBinCounts) {
+    binCounts.resize(m_binUtility.bins(1),
+                     std::vector<unsigned int>(m_binUtility.bins(0), 0));
+  }
+
   // Loop over and fill
   for (std::size_t ib1 = 0; ib1 < m_binUtility.bins(1); ++ib1) {
     for (std::size_t ib0 = 0; ib0 < m_binUtility.bins(0); ++ib0) {
       mpMatrix[ib1][ib0] = m_accumulatedMaterial[ib1][ib0].totalAverage().first;
+      if (recordBinCounts) {
+        binCounts[ib1][ib0] =
+            m_accumulatedMaterial[ib1][ib0].totalAverage().second;
+      }
     }
   }
   // Now return the BinnedSurfaceMaterial
   return std::make_unique<const BinnedSurfaceMaterial>(
-      m_binUtility, std::move(mpMatrix), m_splitFactor);
+      m_binUtility, std::move(mpMatrix), m_splitFactor, binCounts);
 }

Core/src/Material/BinnedSurfaceMaterial.cpp

@@ -17,10 +17,10 @@
 
 namespace Acts {
 
-BinnedSurfaceMaterial::BinnedSurfaceMaterial(const BinUtility& binUtility,
-                                             MaterialSlabVector materialVector,
-                                             double splitFactor,
-                                             MappingType mappingType)
+BinnedSurfaceMaterial::BinnedSurfaceMaterial(
+    const BinUtility& binUtility, MaterialSlabVector materialVector,
+    double splitFactor, std::vector<unsigned int> binCounts,
+    MappingType mappingType)
     : ISurfaceMaterial(splitFactor, mappingType), m_binUtility(binUtility) {
   if (binUtility.dimensions() != 1) {
     throw std::invalid_argument(
@@ -32,15 +32,17 @@ BinnedSurfaceMaterial::BinnedSurfaceMaterial(const BinUtility& binUtility,
         "number of provided material slabs.");
   }
   m_fullMaterial.push_back(std::move(materialVector));
+  m_binCounts.push_back(std::move(binCounts));
 }
 
-BinnedSurfaceMaterial::BinnedSurfaceMaterial(const BinUtility& binUtility,
-                                             MaterialSlabMatrix materialMatrix,
-                                             double splitFactor,
-                                             MappingType mappingType)
+BinnedSurfaceMaterial::BinnedSurfaceMaterial(
+    const BinUtility& binUtility, MaterialSlabMatrix materialMatrix,
+    double splitFactor, std::vector<std::vector<unsigned int>> binCounts,
+    MappingType mappingType)
     : ISurfaceMaterial(splitFactor, mappingType),
       m_binUtility(binUtility),
-      m_fullMaterial(std::move(materialMatrix)) {
+      m_fullMaterial(std::move(materialMatrix)),
+      m_binCounts(std::move(binCounts)) {
   if (binUtility.dimensions() != 1 && binUtility.dimensions() != 2) {
     throw std::invalid_argument(
         "BinnedSurfaceMaterial with material matrix only supports 1D and 2D "

Core/src/Material/BinnedSurfaceMaterialAccumulator.cpp

@@ -160,7 +160,7 @@ Acts::BinnedSurfaceMaterialAccumulator::finalizeMaterial(
   // iterate over the map to call the total average
   for (auto& accMaterial : cState->accumulatedMaterial) {
     ACTS_DEBUG("Finalizing map for Surface " << accMaterial.first);
-    auto sMaterial = accMaterial.second.totalAverage();
+    auto sMaterial = accMaterial.second.totalAverage(m_cfg.recordBinCounts);
     sMaterials[accMaterial.first] = std::move(sMaterial);
   }
 

Examples/Scripts/MaterialMapping/Mat_map_stats.C

@@ -0,0 +1,79 @@
+// This file is part of the ACTS project.
+//
+// Copyright (C) 2026 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 https://mozilla.org/MPL/2.0/.
+
+
+#include <TROOT.h>
+#include <TFile.h>
+#include <TKey.h>
+#include <TDirectory.h>
+#include <TH2F.h>
+#include <TH1F.h>
+#include <TCanvas.h>
+
+#include <string>
+#include <fstream>
+#include <iostream>
+#include <sstream>
+
+void makeBinCountsDist(const char* filename = "material-maps.root")
+{
+  TFile* f = TFile::Open(filename, "READ");
+  if (!f || f->IsZombie()) return;
+
+  TIter nextDir(f->GetListOfKeys());
+  TKey* key;
+
+  while ((key = (TKey*)nextDir())) {
+
+    TObject* obj = key->ReadObj();
+
+    if (!obj->InheritsFrom(TDirectory::Class())){
+      continue;
+    }
+
+    TDirectory* dir = (TDirectory*)obj;
+
+    TH2F* h2 = nullptr;
+    dir->GetObject("binCounts", h2);
+
+    if (!h2){
+        std::cout<<"Histogram not found - continue"<<std::endl;
+        continue;
+    }
+
+    std::cout << "\nProcessing: " << dir->GetName() << std::endl;
+
+    TH1F hDist("binCountsDist",
+               "BinCounts distribution;count;bins",
+               100, 0, h2->GetMaximum() + 1);
+
+    std::size_t nEmptyBins{0};
+
+    for (int ix = 1; ix <= h2->GetNbinsX(); ++ix) {
+      for (int iy = 1; iy <= h2->GetNbinsY(); ++iy) {
+
+        double val = h2->GetBinContent(ix, iy);
+        if(val == 0){
+            nEmptyBins+=1;
+        }
+        hDist.Fill(val);
+
+      }
+    }
+
+    std::cout << "Mean = " << hDist.GetMean()
+              << "  ,StdDev = " << hDist.GetStdDev()
+              <<" ,Empty bins counted = "<< nEmptyBins
+              << std::endl;
+
+    TCanvas* c = new TCanvas();
+    hDist.Draw();
+
+    c->SaveAs(Form("%s_binCounts.png", dir->GetName()));
+  }
+}

Plugins/Json/src/MaterialJsonConverter.cpp

@@ -474,6 +474,7 @@ void Acts::from_json(const nlohmann::json& j,
   Acts::BinUtility bUtility;
   Acts::MaterialSlabMatrix mpMatrix;
   Acts::MappingType mapType = Acts::MappingType::Default;
+  std::vector<std::vector<unsigned int>> binCounts{};
   for (auto& [key, value] : jMaterial.items()) {
     if (key == Acts::jsonKey().binkey && !value.empty()) {
       from_json(value, bUtility);
@@ -491,7 +492,8 @@ void Acts::from_json(const nlohmann::json& j,
   } else if (bUtility.bins() == 1) {
     material = new Acts::HomogeneousSurfaceMaterial(mpMatrix[0][0], 1, mapType);
   } else {
-    material = new Acts::BinnedSurfaceMaterial(bUtility, mpMatrix, 1, mapType);
+    material = new Acts::BinnedSurfaceMaterial(bUtility, mpMatrix, 1, binCounts,
+                                               mapType);
   }
 }
 

Plugins/Root/include/ActsPlugins/Root/RootMaterialMapIo.hpp

@@ -61,6 +61,8 @@ class RootMaterialMapIo {
     std::string maxRangeHistName = "max";
     /// The thickness histogram name
     std::string thicknessHistName = "t";
+    // The bin count histogram name
+    std::string binCountsHistName = "binCounts";
     /// The x0 histogram name
     std::string x0HistName = "x0";
     /// The l0 histogram name

Plugins/Root/src/RootMaterialMapIo.cpp

@@ -227,6 +227,23 @@ void RootMaterialMapIo::fillBinnedSurfaceMaterial(
                         mat.material().massDensity());
     }
   }
+
+  // bin counts writing - if not empty and recorded from the binned surface
+  // material accumulator
+  const auto& binCountsMatrix = bsMaterial.binCounts();
+  if (!binCountsMatrix.empty()) {
+    TH2F binCounts(m_cfg.binCountsHistName.c_str(), "#binCounts ; b0; b1",
+                   bins0, -0.5, fBins0 - 0.5, bins1, -0.5, fBins1 - 0.5);
+
+    for (auto [b1, binVector] : enumerate(binCountsMatrix)) {
+      for (auto [b0, bin] : enumerate(binVector)) {
+        binCounts.SetBinContent(static_cast<int>(b0) + 1,
+                                static_cast<int>(b1) + 1, bin);
+      }
+    }
+    binCounts.Write();
+  }
+
   t.Write();
   x0.Write();
   l0.Write();