diff --git a/Core/include/Acts/Material/AccumulatedSurfaceMaterial.hpp b/Core/include/Acts/Material/AccumulatedSurfaceMaterial.hpp index 85b9310afd3..42828685d2d 100644 --- a/Core/include/Acts/Material/AccumulatedSurfaceMaterial.hpp +++ b/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 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 totalAverage( + bool recordBinCounts = false); /// Access to the accumulated material /// @return Reference to the matrix of accumulated material data diff --git a/Core/include/Acts/Material/BinnedSurfaceMaterial.hpp b/Core/include/Acts/Material/BinnedSurfaceMaterial.hpp index 65427a8426c..9bdaea53970 100644 --- a/Core/include/Acts/Material/BinnedSurfaceMaterial.hpp +++ b/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 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> 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>& 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> m_binCounts; }; } // namespace Acts diff --git a/Core/include/Acts/Material/BinnedSurfaceMaterialAccumulator.hpp b/Core/include/Acts/Material/BinnedSurfaceMaterialAccumulator.hpp index 0c5a35c8105..d3a4082fa6a 100644 --- a/Core/include/Acts/Material/BinnedSurfaceMaterialAccumulator.hpp +++ b/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 materialSurfaces = {}; }; diff --git a/Core/src/Material/AccumulatedSurfaceMaterial.cpp b/Core/src/Material/AccumulatedSurfaceMaterial.cpp index 9fb08edc903..f9b51a243b8 100644 --- a/Core/src/Material/AccumulatedSurfaceMaterial.cpp +++ b/Core/src/Material/AccumulatedSurfaceMaterial.cpp @@ -132,7 +132,7 @@ void Acts::AccumulatedSurfaceMaterial::trackAverage( /// Total average creates SurfaceMaterial std::unique_ptr -Acts::AccumulatedSurfaceMaterial::totalAverage() { +Acts::AccumulatedSurfaceMaterial::totalAverage(bool recordBinCounts) { if (m_binUtility.bins() == 1) { // Return HomogeneousSurfaceMaterial return std::make_unique( @@ -142,13 +142,25 @@ Acts::AccumulatedSurfaceMaterial::totalAverage() { MaterialSlabMatrix mpMatrix( m_binUtility.bins(1), MaterialSlabVector(m_binUtility.bins(0), MaterialSlab::Nothing())); + + std::vector> binCounts{}; + + if (recordBinCounts) { + binCounts.resize(m_binUtility.bins(1), + std::vector(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( - m_binUtility, std::move(mpMatrix), m_splitFactor); + m_binUtility, std::move(mpMatrix), m_splitFactor, binCounts); } diff --git a/Core/src/Material/BinnedSurfaceMaterial.cpp b/Core/src/Material/BinnedSurfaceMaterial.cpp index 303935973ad..d7dd4120beb 100644 --- a/Core/src/Material/BinnedSurfaceMaterial.cpp +++ b/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 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> 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 " diff --git a/Core/src/Material/BinnedSurfaceMaterialAccumulator.cpp b/Core/src/Material/BinnedSurfaceMaterialAccumulator.cpp index 096421bffcd..82735d5500a 100644 --- a/Core/src/Material/BinnedSurfaceMaterialAccumulator.cpp +++ b/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); } diff --git a/Examples/Scripts/MaterialMapping/material_map_stats.py b/Examples/Scripts/MaterialMapping/material_map_stats.py new file mode 100644 index 00000000000..219f71c149d --- /dev/null +++ b/Examples/Scripts/MaterialMapping/material_map_stats.py @@ -0,0 +1,69 @@ +import ROOT + + +def make_bin_counts_dist(filename="material-maps.root"): + + f = ROOT.TFile.Open(filename, "READ") + if not f or f.IsZombie(): + return + + for key in f.GetListOfKeys(): + + obj = key.ReadObj() + + if not obj.InheritsFrom("TDirectory"): + continue + + directory = obj + + h2 = directory.Get("binCounts") + + if not h2: + print("Histogram not found - continue") + continue + + print(f"\nProcessing: {directory.GetName()}") + + hdist = ROOT.TH1F( + "binCountsDist", + "BinCounts distribution;count;bins", + 20, + 0, + h2.GetMaximum() + 1, + ) + + n_empty_bins = 0 + + for ix in range(1, h2.GetNbinsX() + 1): + for iy in range(1, h2.GetNbinsY() + 1): + + val = h2.GetBinContent(ix, iy) + + if val == 0: + n_empty_bins += 1 + + hdist.Fill(val) + + print( + f"Mean = {hdist.GetMean():.3f}, " + f"StdDev = {hdist.GetStdDev():.3f}, " + f"Empty bins counted = {n_empty_bins}" + ) + + c = ROOT.TCanvas(f"{directory.GetName()}binCounts", "", 800, 600) + + hdist.Draw() + n_total_bins = h2.GetNbinsX() * h2.GetNbinsY() + empty_fraction = n_empty_bins / n_total_bins + pt = ROOT.TPaveText(0.60, 0.70, 0.90, 0.88, "NDC") + pt.AddText(f"Empty bins: {n_empty_bins}") + pt.AddText(f"Empty fraction: {100.0 * empty_fraction:.2f}%") + pt.AddText(f"Mean: {hdist.GetMean():.2f}") + pt.AddText(f"StdDev: {hdist.GetStdDev():.2f}") + pt.Draw() + + c.SaveAs(f"{directory.GetName()}_binCounts.png") + + +if __name__ == "__main__": + make_bin_counts_dist("material-maps.root") diff --git a/Plugins/Json/src/MaterialJsonConverter.cpp b/Plugins/Json/src/MaterialJsonConverter.cpp index c21b65fc3ce..1819a5ad269 100644 --- a/Plugins/Json/src/MaterialJsonConverter.cpp +++ b/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> 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); } } diff --git a/Plugins/Root/include/ActsPlugins/Root/RootMaterialMapIo.hpp b/Plugins/Root/include/ActsPlugins/Root/RootMaterialMapIo.hpp index d40324b3f60..d73de105cfa 100644 --- a/Plugins/Root/include/ActsPlugins/Root/RootMaterialMapIo.hpp +++ b/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 diff --git a/Plugins/Root/src/RootMaterialMapIo.cpp b/Plugins/Root/src/RootMaterialMapIo.cpp index 99f5f697940..a07072efc4c 100644 --- a/Plugins/Root/src/RootMaterialMapIo.cpp +++ b/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(b0) + 1, + static_cast(b1) + 1, bin); + } + } + binCounts.Write(); + } + t.Write(); x0.Write(); l0.Write();