feat: add elemental fractions to the material recording step. add a toggle to enable / disable

Author: Lauren Wills

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

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 recordElementFractions = 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 recordElementFractions = 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.recordElementFractions = m_cfg.recordElementFractions;
     // G4RunManager will take care of deletion
     auto steppingAction = new Geant4::MaterialSteppingAction(
         steppingCfg, this->logger().cloneWithSuffix("MaterialSteppingAction"));

Examples/Algorithms/Geant4/src/MaterialSteppingAction.cpp

@@ -90,6 +90,21 @@ void MaterialSteppingAction::UserSteppingAction(const G4Step* stepPtr) {
   mInteraction.materialSlab = slab;
   mInteraction.pathCorrection = step.GetStepLength() * convertLengthToActs;
 
+  if (m_cfg.recordElementFractions) {
+    if (nElements == 1) {
+	mInteraction.elementZ.push_back(
+      	static_cast<unsigned int>(material.GetZ()));
+    	mInteraction.elementFrac.push_back(1.0f);
+    } else {
+    	for (std::size_t i = 0; i < nElements; i++){
+      	mInteraction.elementZ.push_back(
+        static_cast<unsigned int>(elements->at(i).GetZ()));
+      	mInteraction.elementFrac.push_back(
+        static_cast<float>(fraction[i]));
+    	}
+    }
+  }
+
   assert(step.GetTrack() != nullptr);
   const G4Track& track = *step.GetTrack();
   const std::size_t trackId = track.GetTrackID();

Examples/Scripts/Python/material_recording.py

@@ -75,6 +75,7 @@ def runMaterialRecording(
         randomNumbers=rnd,
         inputParticles=hepmc3Converter.config.outputParticles,
         outputMaterialTracks=materialTrackCollectionName,
+        recordElementFractions= False,
     )
 
     s.addAlgorithm(g4Alg)

Plugins/Root/include/ActsPlugins/Root/RootMaterialTrackIo.hpp

@@ -155,6 +155,12 @@ class RootMaterialTrackIo {
     /// step material rho
     std::vector<float> stepMatRho;
     std::vector<float>* stepMatRhoPtr = &stepMatRho;
+    /// step material element
+    std::vector<std::vector<unsigned int>> stepElementZ;
+    std::vector<std::vector<unsigned int>>* stepElementZPtr = &stepElementZ;
+    /// step material fraction of elements
+    std::vector<std::vector<float>> stepFraction;
+    std::vector<std::vector<float>>* stepFractionPtr = &stepFraction;
   };
 
   struct MaterialSurfacePayload {

Plugins/Root/src/RootMaterialTrackIo.cpp

@@ -43,6 +43,8 @@ void ActsPlugins::RootMaterialTrackIo::connectForRead(TChain& materialChain) {
   materialChain.SetBranchAddress("mat_A", &m_stepPayload.stepMatAPtr);
   materialChain.SetBranchAddress("mat_Z", &m_stepPayload.stepMatZPtr);
   materialChain.SetBranchAddress("mat_rho", &m_stepPayload.stepMatRhoPtr);
+  materialChain.SetBranchAddress("elements", &m_stepPayload.stepElementZPtr);
+  materialChain.SetBranchAddress("fraction", &m_stepPayload.stepFractionPtr);
   if (m_cfg.surfaceInfo) {
     materialChain.SetBranchAddress("sur_id", &m_surfacePayload.surfaceIdPtr);
     materialChain.SetBranchAddress("sur_x", &m_surfacePayload.surfaceXPtr);
@@ -80,6 +82,8 @@ void ActsPlugins::RootMaterialTrackIo::connectForWrite(TTree& materialTree) {
   materialTree.Branch("mat_A", &m_stepPayload.stepMatA);
   materialTree.Branch("mat_Z", &m_stepPayload.stepMatZ);
   materialTree.Branch("mat_rho", &m_stepPayload.stepMatRho);
+  materialTree.Branch("elements", &m_stepPayload.stepElementZ);
+  materialTree.Branch("fraction", &m_stepPayload.stepFraction);
 
   if (m_cfg.prePostStepInfo) {
     materialTree.Branch("mat_sx", &m_stepPayload.stepXs);
@@ -128,6 +132,8 @@ void ActsPlugins::RootMaterialTrackIo::write(
   m_stepPayload.stepMatA.clear();
   m_stepPayload.stepMatZ.clear();
   m_stepPayload.stepMatRho.clear();
+  m_stepPayload.stepElementZ.clear();
+  m_stepPayload.stepFraction.clear();
 
   m_surfacePayload.surfaceId.clear();
   m_surfacePayload.surfaceX.clear();
@@ -162,6 +168,8 @@ void ActsPlugins::RootMaterialTrackIo::write(
   m_stepPayload.stepMatA.reserve(mints);
   m_stepPayload.stepMatZ.reserve(mints);
   m_stepPayload.stepMatRho.reserve(mints);
+  m_stepPayload.stepElementZ.reserve(mints);
+  m_stepPayload.stepFraction.reserve(mints);
 
   m_surfacePayload.surfaceId.reserve(mints);
   m_surfacePayload.surfaceX.reserve(mints);
@@ -274,6 +282,8 @@ void ActsPlugins::RootMaterialTrackIo::write(
     m_stepPayload.stepMatA.push_back(mprops.material().Ar());
     m_stepPayload.stepMatZ.push_back(mprops.material().Z());
     m_stepPayload.stepMatRho.push_back(mprops.material().massDensity());
+    m_stepPayload.stepElementZ.push_back(mint.elementZ);
+    m_stepPayload.stepFraction.push_back(mint.elementFrac);
     // re-calculate if defined to do so
     if (m_cfg.recalculateTotals) {
       m_summaryPayload.tX0 += mprops.thicknessInX0();
@@ -320,6 +330,11 @@ RecordedMaterialTrack ActsPlugins::RootMaterialTrackIo::read() const {
                                   m_stepPayload.stepMatZ[is],
                                   m_stepPayload.stepMatRho[is]),
         s);
+    /// adding the element vectors and fractions 
+    if (is < m_stepPayload.stepElementZ.size()) {
+      mInteraction.elementZ = m_stepPayload.stepElementZ[is];
+      mInteraction.elementFrac = m_stepPayload.stepFraction[is];
+    }
     if (m_cfg.surfaceInfo) {
       // add the surface information to the interaction this allows the
       // mapping to be speed up

Python/Examples/src/plugins/Geant4.cpp

@@ -168,7 +168,7 @@ PYBIND11_MODULE(ActsExamplesPythonBindingsGeant4, mod) {
     auto c = py::class_<Algorithm::Config, Geant4SimulationBase::Config,
                         std::shared_ptr<Algorithm::Config>>(alg, "Config")
                  .def(py::init<>());
-    ACTS_PYTHON_STRUCT(c, outputMaterialTracks, excludeMaterials);
+    ACTS_PYTHON_STRUCT(c, outputMaterialTracks, excludeMaterials, recordElementFractions);
   }
 
   {