feat: Add elemental fractions to material recording

Core/include/Acts/Material/AccumulatedSurfaceMaterial.hpp

@@ -92,22 +92,30 @@ class AccumulatedSurfaceMaterial {
   /// @param lp local position for the bin assignment
   /// @param mp material properties to be assigned
   /// @param pathCorrection Correction factor for the effective path length
+  /// @param elementZ to list all elements' Z value present (optional)
+  /// @param elementFrac to list the fraction that each element makes up (optional)
   ///
   /// @return the bin triple to which the material was assigned
   std::array<std::size_t, 3> accumulate(const Vector2& lp,
                                         const MaterialSlab& mp,
-                                        double pathCorrection = 1.);
+                                        double pathCorrection = 1.,
+					const std::vector<unsigned int>& elementZ = {},
+					const std::vector<float>& elementFrac = {});
 
   /// Assign a material properties object
   ///
   /// @param gp global position for the bin assignment
   /// @param mp material properties to be assigned
   /// @param pathCorrection Correction factor for the effective path length
+  /// @param elementZ to list all elements' Z value present (optional)
+  /// @param elementFrac to list the fraction that each element makes up (optional)
   ///
   /// @return the bin triple to which the material was assigned
   std::array<std::size_t, 3> accumulate(const Vector3& gp,
                                         const MaterialSlab& mp,
-                                        double pathCorrection = 1.);
+                                        double pathCorrection = 1.,
+					const std::vector<unsigned int>& elementZ = {},
+                                        const std::vector<float>& elementFrac = {});
 
   /// Use the accumulated material to update the material variance
   ///
@@ -162,6 +170,23 @@ class AccumulatedSurfaceMaterial {
 
   /// The stored accumulated material matrix
   AccumulatedMatrix m_accumulatedMaterial;
+
+  /// The running sum of thickness per bin
+  /// Used to normalize the vector of fractions or element
+  std::vector<std::vector<double>> m_totalThickness;
+
+  /// The list of atomic numbers of elements present during mapping, per bin
+  std::vector<std::vector<std::vector<unsigned int>>> m_elementZ;
+
+  /// The fraction of each element present, normalized
+  /// Weighted by m_totalThickness to get final fraction
+  std::vector<std::vector<std::vector<float>>> m_weightedFrac;
+
+  void accumulateElementData(
+    std::size_t bin0, std::size_t bin1,
+    const MaterialSlab& mp, double pathCorrection,
+    const std::vector<unsigned int>& elementZ,
+    const std::vector<float>& elementFrac);
 };
 
 inline const BinUtility& AccumulatedSurfaceMaterial::binUtility() const {

Core/include/Acts/Material/BinnedSurfaceMaterial.hpp

@@ -27,6 +27,17 @@ class BinnedSurfaceMaterial : public ISurfaceMaterial {
   /// Default Constructor - deleted
   BinnedSurfaceMaterial() = delete;
 
+
+  ///Alias for Grid Level values of elementZ vector
+  /// [bin0][bin1] -> Z values per bin
+  using ElementZMatrix = 
+      std::vector<std::vector<std::vector<unsigned int>>>;
+
+  ///Alias for Grid Level values of ElementFrac vetcor
+  using ElementFracMatrix = 
+      std::vector<std::vector<std::vector<float>>>;
+
+
   /// Explicit constructor with only full MaterialSlab,
   /// for one-dimensional binning.
   ///
@@ -42,7 +53,9 @@ class BinnedSurfaceMaterial : public ISurfaceMaterial {
   BinnedSurfaceMaterial(const BinUtility& binUtility,
                         MaterialSlabVector fullProperties,
                         double splitFactor = 0.,
-                        MappingType mappingType = MappingType::Default);
+                        MappingType mappingType = MappingType::Default,
+			ElementZMatrix elementZ = {},
+			ElementFracMatrix elementFrac = {});
 
   /// Explicit constructor with only full MaterialSlab,
   /// for two-dimensional binning.
@@ -59,7 +72,9 @@ class BinnedSurfaceMaterial : public ISurfaceMaterial {
   BinnedSurfaceMaterial(const BinUtility& binUtility,
                         MaterialSlabMatrix fullProperties,
                         double splitFactor = 0.,
-                        MappingType mappingType = MappingType::Default);
+                        MappingType mappingType = MappingType::Default,
+			ElementZMatrix elementZ = {},
+			ElementFracMatrix elementFrac = {});
 
   /// Copy Move Constructor
   ///
@@ -106,6 +121,25 @@ class BinnedSurfaceMaterial : public ISurfaceMaterial {
 
   using ISurfaceMaterial::materialSlab;
 
+  ///Return 3D Grid of element Z values
+  ///@return Reference to the ElementZMatrix
+  const ElementZMatrix& elementZMatrix() const;
+
+  ///Return 3D Grid of fractions for each element Z
+  ///@return Reference to the ElementFracMatrix
+  const ElementFracMatrix& elementFracMatrix() const;
+
+  ///Return elementZ for a specific bin
+  ///@return Reference to the elementZ vector for that bin
+  const std::vector<unsigned int>& elementZAt(std::size_t bin0,
+		 std::size_t bin1) const;
+
+  ///Return elementFrac for a specific bin
+  ///@return Reference to the elementFrac vector for that bin
+  const std::vector<float>& elementFracAt(std::size_t bin0,
+		std::size_t bin1) const; 
+
+
   /// Output Method for std::ostream, to be overloaded by child classes
   /// @param sl The output stream to write to
   /// @return Reference to the output stream after writing
@@ -117,6 +151,14 @@ class BinnedSurfaceMaterial : public ISurfaceMaterial {
 
   /// The five different MaterialSlab
   MaterialSlabMatrix m_fullMaterial;
+
+  ///The Per-bin element list - indexed[bin1][bin0]
+  ///empty if using older material recording (element list was not recorded)
+  ElementZMatrix m_elementZ;
+
+  ///The Per-bin list of fractions weighted by thickness- sum to 1 per bin
+  /// Empty if using older material recording (fraction was not recorded)
+  ElementFracMatrix m_elementFrac;
 };
 
 inline const BinUtility& BinnedSurfaceMaterial::binUtility() const {
@@ -127,4 +169,24 @@ inline const MaterialSlabMatrix& BinnedSurfaceMaterial::fullMaterial() const {
   return m_fullMaterial;
 }
 
+inline const BinnedSurfaceMaterial::ElementZMatrix&
+BinnedSurfaceMaterial::elementZMatrix() const {
+  return m_elementZ;
+}
+
+inline const BinnedSurfaceMaterial::ElementFracMatrix&
+BinnedSurfaceMaterial::elementFracMatrix() const {
+  return m_elementFrac;
+}
+
+inline const std::vector<unsigned int>& BinnedSurfaceMaterial::elementZAt(
+    std::size_t bin0, std::size_t bin1) const{
+  return m_elementZ[bin0][bin1];
+}
+
+inline const std::vector<float>& BinnedSurfaceMaterial::elementFracAt(
+    std::size_t bin0, std::size_t bin1) const{
+  return m_elementFrac[bin0][bin1];
+}
+
 }  // namespace Acts

Core/include/Acts/Material/MaterialInteraction.hpp

@@ -81,6 +81,10 @@ struct MaterialInteraction {
   double pathCorrection = 1.;
   /// The effective, passed material properties including the path correction.
   MaterialSlab materialSlab = MaterialSlab::Nothing();
+  /// Vectors of the individual elements present
+  std::vector<unsigned int> elementZ = {};
+  /// How much of each element are present
+  std::vector<float> elementFrac = {};
 };
 
 /// Simple result struct to be returned

Core/src/Material/AccumulatedSurfaceMaterial.cpp

@@ -19,6 +19,13 @@ Acts::AccumulatedSurfaceMaterial::AccumulatedSurfaceMaterial(double splitFactor)
     : m_splitFactor(splitFactor) {
   AccumulatedVector accMat = {{AccumulatedMaterialSlab()}};
   m_accumulatedMaterial = {{accMat}};
+  /// initializing the three vectors needed for Detailed Material description:
+  m_totalThickness = std::vector<std::vector<double>>(
+	1, std::vector<double>(1, 0.0));
+  m_elementZ = std::vector<std::vector<std::vector<unsigned int>>>(
+	1, std::vector<std::vector<unsigned int>>(1));
+  m_weightedFrac = std::vector<std::vector<std::vector<float>>>(
+	1, std::vector<std::vector<float>>(1));
 }
 
 // Binned Material constructor with split factor
@@ -29,33 +36,87 @@ Acts::AccumulatedSurfaceMaterial::AccumulatedSurfaceMaterial(
   std::size_t bins1 = m_binUtility.bins(1);
   AccumulatedVector accVec(bins0, AccumulatedMaterialSlab());
   m_accumulatedMaterial = AccumulatedMatrix(bins1, accVec);
+  m_totalThickness = std::vector<std::vector<double>>(
+	bins1, std::vector<double>(bins0, 0.0));
+  m_elementZ = std::vector<std::vector<std::vector<unsigned int>>>(
+	bins1, std::vector<std::vector<unsigned int>>(bins0));
+  m_weightedFrac = std::vector<std::vector<std::vector<float>>>(
+	bins1, std::vector<std::vector<float>>(bins0));
 }
 
 // Assign a material properties object
 std::array<std::size_t, 3> Acts::AccumulatedSurfaceMaterial::accumulate(
-    const Vector2& lp, const MaterialSlab& mp, double pathCorrection) {
+    const Vector2& lp, const MaterialSlab& mp, double pathCorrection,
+    const std::vector<unsigned int>& elementZ,
+    const std::vector<float>& elementFrac) {
   if (m_binUtility.dimensions() == 0) {
     m_accumulatedMaterial[0][0].accumulate(mp, pathCorrection);
+    accumulateElementData(0, 0, mp, pathCorrection, elementZ, elementFrac);
     return {0, 0, 0};
   }
   std::size_t bin0 = m_binUtility.bin(lp, 0);
   std::size_t bin1 = m_binUtility.bin(lp, 1);
   m_accumulatedMaterial[bin1][bin0].accumulate(mp, pathCorrection);
+  accumulateElementData(bin0, bin1, mp, pathCorrection, elementZ, elementFrac);
   return {bin0, bin1, 0};
 }
 
+
 // Assign a material properties object
 std::array<std::size_t, 3> Acts::AccumulatedSurfaceMaterial::accumulate(
-    const Vector3& gp, const MaterialSlab& mp, double pathCorrection) {
+    const Vector3& gp, const MaterialSlab& mp, double pathCorrection,
+    const std::vector<unsigned int>& elementZ,
+    const std::vector<float>& elementFrac) {
   if (m_binUtility.dimensions() == 0) {
     m_accumulatedMaterial[0][0].accumulate(mp, pathCorrection);
+    accumulateElementData(0, 0, mp, pathCorrection, elementZ, elementFrac);
     return {0, 0, 0};
   }
   std::array<std::size_t, 3> bTriple = m_binUtility.binTriple(gp);
   m_accumulatedMaterial[bTriple[1]][bTriple[0]].accumulate(mp, pathCorrection);
+  accumulateElementData(bTriple[0],bTriple[1], mp, pathCorrection, elementZ, elementFrac);
   return bTriple;
 }
 
+// Helper- new element list and fraction accumulation calculation
+void Acts::AccumulatedSurfaceMaterial::accumulateElementData(
+    std::size_t bin0, size_t bin1, 
+    const MaterialSlab& mp, double pathCorrection,
+    const std::vector<unsigned int>& elementZ,
+    const std::vector<float>& elementFrac) {
+  if (elementZ.empty() || elementZ.size() != elementFrac.size()) {
+    return;
+  }
+
+  double stepThickness = 
+	  static_cast<double>(mp.thickness()) / pathCorrection;
+
+  if (stepThickness <= 0.) {
+	  return;
+  }
+
+  m_totalThickness[bin1][bin0] += stepThickness;
+  for (std::size_t i = 0; i <elementZ.size(); ++i) {
+    unsigned int z = elementZ[i];
+    float frac = elementFrac[i];
+    /// check to see if this element is already in the accumulated list
+    auto it = std::find(m_elementZ[bin1][bin0].begin(),
+			  m_elementZ[bin1][bin0].end(), z);
+    if (it != m_elementZ[bin1][bin0].end()) {
+      std::size_t idx =
+	  std::distance(m_elementZ[bin1][bin0].begin(), it);
+      m_weightedFrac[bin1][bin0][idx] += 
+	      frac * static_cast<float>(stepThickness);
+    } else{
+      m_elementZ[bin1][bin0].push_back(z);
+      m_weightedFrac[bin1][bin0].push_back(
+		      frac * static_cast<float>(stepThickness));
+    }
+  }
+}
+
+
+
 // Void average for vacuum assignment
 void Acts::AccumulatedSurfaceMaterial::trackVariance(const Vector3& gp,
                                                      MaterialSlab slabReference,
@@ -148,7 +209,31 @@ Acts::AccumulatedSurfaceMaterial::totalAverage() {
       mpMatrix[ib1][ib0] = m_accumulatedMaterial[ib1][ib0].totalAverage().first;
     }
   }
+  // Build the element composition grids
+  BinnedSurfaceMaterial::ElementZMatrix finalElementZ(
+	m_binUtility.bins(1),
+	std::vector<std::vector<unsigned int>>(m_binUtility.bins(0)));
+  BinnedSurfaceMaterial::ElementFracMatrix finalElementFrac(
+	m_binUtility.bins(1),
+	std::vector<std::vector<float>>(m_binUtility.bins(0)));
+  for (std::size_t ib1 = 0; ib1 < m_binUtility.bins(1); ++ib1){
+    for (std::size_t ib0 = 0; ib0 < m_binUtility.bins(0); ++ib0){
+      finalElementZ[ib1][ib0] = m_elementZ[ib1][ib0];
+      double totalThickness = m_totalThickness[ib1][ib0];
+      if (totalThickness > 0. && !m_weightedFrac[ib1][ib0].empty()) {
+        finalElementFrac[ib1][ib0].resize(m_weightedFrac[ib1][ib0].size());
+	for (std::size_t i = 0; i < m_weightedFrac[ib1][ib0].size(); ++i){
+	  finalElementFrac[ib1][ib0][i] = 
+		m_weightedFrac[ib1][ib0][i] / static_cast<float>(totalThickness);
+	}
+      }
+    }
+  }
+
+
   // Now return the BinnedSurfaceMaterial
   return std::make_unique<const BinnedSurfaceMaterial>(
-      m_binUtility, std::move(mpMatrix), m_splitFactor);
+      m_binUtility, std::move(mpMatrix), m_splitFactor,
+      MappingType::Default,
+      std::move(finalElementZ), std::move(finalElementFrac));
 }

Core/src/Material/BinnedSurfaceMaterial.cpp

@@ -16,18 +16,25 @@
 
 Acts::BinnedSurfaceMaterial::BinnedSurfaceMaterial(
     const BinUtility& binUtility, MaterialSlabVector fullProperties,
-    double splitFactor, Acts::MappingType mappingType)
-    : ISurfaceMaterial(splitFactor, mappingType), m_binUtility(binUtility) {
-  // fill the material with deep copy
+    double splitFactor, Acts::MappingType mappingType,
+    ElementZMatrix elementZ, ElementFracMatrix elementFrac)
+    : ISurfaceMaterial(splitFactor, mappingType), 
+	m_binUtility(binUtility),
+	m_elementZ(std::move(elementZ)),
+	m_elementFrac(std::move(elementFrac)) {
   m_fullMaterial.push_back(std::move(fullProperties));
 }
 
 Acts::BinnedSurfaceMaterial::BinnedSurfaceMaterial(
     const BinUtility& binUtility, MaterialSlabMatrix fullProperties,
-    double splitFactor, Acts::MappingType mappingType)
+    double splitFactor, Acts::MappingType mappingType,
+    ElementZMatrix elementZ, ElementFracMatrix elementFrac)
     : ISurfaceMaterial(splitFactor, mappingType),
       m_binUtility(binUtility),
-      m_fullMaterial(std::move(fullProperties)) {}
+      m_fullMaterial(std::move(fullProperties)),
+      m_elementZ(std::move(elementZ)),
+      m_elementFrac(std::move(elementFrac)) {
+}
 
 Acts::BinnedSurfaceMaterial& Acts::BinnedSurfaceMaterial::scale(double factor) {
   for (auto& materialVector : m_fullMaterial) {

Core/src/Material/SurfaceMaterialMapper.cpp

@@ -390,7 +390,8 @@ Result<void> SurfaceMaterialMapper::mapInteraction(
     }
     // Now assign the material for the accumulation process
     auto tBin = currentAccMaterial->second.accumulate(
-        currentPos, rmIter->materialSlab, currentPathCorrection);
+        currentPos, rmIter->materialSlab, currentPathCorrection,
+	rmIter->elementZ, rmIter->elementFrac);
     if (!touchedMapBins.contains(&(currentAccMaterial->second))) {
       touchedMapBins.insert(MapBin(&(currentAccMaterial->second), tBin));
     }
@@ -482,7 +483,8 @@ void SurfaceMaterialMapper::mapSurfaceInteraction(
 
     // Now assign the material for the accumulation process
     auto tBin = currentAccMaterial->second.accumulate(
-        currentPos, rmIter->materialSlab, rmIter->pathCorrection);
+        currentPos, rmIter->materialSlab, rmIter->pathCorrection,
+	rmIter->elementZ, rmIter->elementFrac);
     if (!touchedMapBins.contains(&(currentAccMaterial->second))) {
       touchedMapBins.insert(MapBin(&(currentAccMaterial->second), tBin));
     }

Examples/Algorithms/Geant4/include/ActsExamples/Geant4/Geant4Simulation.hpp

@@ -191,6 +191,7 @@ class Geant4MaterialRecording final : public Geant4SimulationBase {
 
     /// Materials to exclude from the recording.
     std::vector<std::string> excludeMaterials = {"Air", "Vacuum"};
+    bool DetailedMaterial = false;
   };
 
   /// Material recording constructor

Examples/Algorithms/Geant4/include/ActsExamples/Geant4/MaterialSteppingAction.hpp

@@ -35,6 +35,7 @@ class MaterialSteppingAction final : public G4UserSteppingAction {
     std::shared_ptr<EventStore> eventStore;
 
     std::vector<std::string> excludeMaterials = {};
+    bool DetailedMaterial = false;
   };
 
   /// Construct the action

Examples/Algorithms/Geant4/src/Geant4Simulation.cpp

@@ -392,6 +392,7 @@ Geant4MaterialRecording::Geant4MaterialRecording(
     Geant4::MaterialSteppingAction::Config steppingCfg;
     steppingCfg.eventStore = m_eventStore;
     steppingCfg.excludeMaterials = m_cfg.excludeMaterials;
+    steppingCfg.DetailedMaterial = m_cfg.DetailedMaterial;
     // G4RunManager will take care of deletion
     auto steppingAction = new Geant4::MaterialSteppingAction(
         steppingCfg, this->logger().cloneWithSuffix("MaterialSteppingAction"));