feat: adding G4 stepping

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

@@ -9,6 +9,7 @@
 #pragma once
 
 #include "Acts/Material/MaterialInteraction.hpp"
+#include "ActsExamples/EventData/PropagationSummary.hpp"
 #include "ActsExamples/EventData/SimHit.hpp"
 #include "ActsExamples/EventData/SimParticle.hpp"
 #include "ActsExamples/Framework/DataHandle.hpp"
@@ -52,6 +53,9 @@ struct EventStore {
   /// Tracks recorded in material mapping
   std::unordered_map<std::size_t, Acts::RecordedMaterialTrack> materialTracks;
 
+  /// Stepping records for step plotting
+  std::unordered_map<G4int, PropagationSummary> propagationRecords;
+
   /// Particle hit count (for hit indexing)
   std::unordered_map<SimBarcode, std::size_t> particleHitCount;
   /// Particle status

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

@@ -10,6 +10,7 @@
 
 #include "Acts/Material/MaterialInteraction.hpp"
 #include "Acts/Utilities/Logger.hpp"
+#include "ActsExamples/EventData/PropagationSummary.hpp"
 #include "ActsExamples/EventData/SimHit.hpp"
 #include "ActsExamples/EventData/SimParticle.hpp"
 #include "ActsExamples/Framework/DataHandle.hpp"
@@ -124,6 +125,9 @@ class Geant4Simulation final : public Geant4SimulationBase {
     /// Name of the output collection : final particles
     std::string outputParticlesFinal = "particles_final";
 
+    /// Name of the output collection : propagation records (debugging)
+    std::string outputPropagationSummaries = "propagation_summaries";
+
     /// The ACTS sensitive surfaces in a mapper, used for hit creation
     std::shared_ptr<const SensitiveSurfaceMapper> sensitiveSurfaceMapper =
         nullptr;
@@ -138,11 +142,17 @@ class Geant4Simulation final : public Geant4SimulationBase {
     double killAfterTime = std::numeric_limits<double>::infinity();
     bool killSecondaries = false;
 
+    bool recordHitsOfCharged = true;
+
     bool recordHitsOfNeutrals = false;
 
+    bool recordHitsOfPrimaries = true;
+
     bool recordHitsOfSecondaries = true;
 
     bool keepParticlesWithoutHits = true;
+
+    bool recordPropagationSummaries = false;
   };
 
   /// Simulation constructor
@@ -175,6 +185,9 @@ class Geant4Simulation final : public Geant4SimulationBase {
   WriteDataHandle<SimParticleContainer> m_outputParticlesFinal{
       this, "OutputParticlesFinal"};
   WriteDataHandle<SimHitContainer> m_outputSimHits{this, "OutputSimHIts"};
+
+  WriteDataHandle<PropagationSummaries> m_outputPropagationSummaries{
+      this, "OutputPropagationSummaries"};
 };
 
 class Geant4MaterialRecording final : public Geant4SimulationBase {

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

@@ -41,6 +41,9 @@ class SensitiveSteppingAction : public G4UserSteppingAction {
     bool neutral = false;
     bool primary = true;
     bool secondary = true;
+
+    /// step logging mode
+    bool stepLogging = false;
   };
 
   /// Construct the stepping action

Examples/Algorithms/Geant4/src/Geant4Simulation.cpp

@@ -233,10 +233,11 @@ ActsExamples::Geant4Simulation::Geant4Simulation(const Config& cfg,
 
     SensitiveSteppingAction::Config stepCfg;
     stepCfg.eventStore = m_eventStore;
-    stepCfg.charged = true;
+    stepCfg.charged = cfg.recordHitsOfCharged;
     stepCfg.neutral = cfg.recordHitsOfNeutrals;
-    stepCfg.primary = true;
+    stepCfg.primary = cfg.recordHitsOfPrimaries;
     stepCfg.secondary = cfg.recordHitsOfSecondaries;
+    stepCfg.stepLogging = cfg.recordPropagationSummaries;
 
     SteppingActionList::Config steppingCfg;
     steppingCfg.actions.push_back(std::make_unique<ParticleKillAction>(
@@ -300,6 +301,10 @@ ActsExamples::Geant4Simulation::Geant4Simulation(const Config& cfg,
   m_outputSimHits.initialize(cfg.outputSimHits);
   m_outputParticlesInitial.initialize(cfg.outputParticlesInitial);
   m_outputParticlesFinal.initialize(cfg.outputParticlesFinal);
+
+  if (cfg.recordPropagationSummaries) {
+    m_outputPropagationSummaries.initialize(cfg.outputPropagationSummaries);
+  }
 }
 
 ActsExamples::Geant4Simulation::~Geant4Simulation() = default;
@@ -330,6 +335,16 @@ ActsExamples::ProcessCode ActsExamples::Geant4Simulation::execute(
       ctx, SimHitContainer(eventStore().hits.begin(), eventStore().hits.end()));
 #endif
 
+  // Output the propagation summaries if requested
+  if (m_cfg.recordPropagationSummaries) {
+    PropagationSummaries summaries;
+    summaries.reserve(eventStore().propagationRecords.size());
+    for (auto& [trackId, summary] : eventStore().propagationRecords) {
+      summaries.push_back(std::move(summary));
+    }
+    m_outputPropagationSummaries(ctx, std::move(summaries));
+  }
+
   return ActsExamples::ProcessCode::SUCCESS;
 }
 

Examples/Algorithms/Geant4/src/SensitiveSteppingAction.cpp

@@ -11,6 +11,7 @@
 #include "Acts/Definitions/Algebra.hpp"
 #include "Acts/Definitions/Units.hpp"
 #include "Acts/Geometry/GeometryIdentifier.hpp"
+#include "Acts/Propagator/detail/SteppingLogger.hpp"
 #include "Acts/Surfaces/Surface.hpp"
 #include "Acts/Utilities/MultiIndex.hpp"
 #include "ActsExamples/EventData/SimHit.hpp"
@@ -19,6 +20,7 @@
 #include "ActsFatras/EventData/Barcode.hpp"
 
 #include <algorithm>
+#include <array>
 #include <cstddef>
 #include <string>
 #include <unordered_map>
@@ -53,46 +55,59 @@ BOOST_DESCRIBE_ENUM(G4TrackStatus, fAlive, fStopButAlive, fStopAndKill,
 
 namespace {
 
-ActsFatras::Hit hitFromStep(const G4StepPoint* preStepPoint,
-                            const G4StepPoint* postStepPoint,
-                            ActsFatras::Barcode particleId,
-                            Acts::GeometryIdentifier geoId,
-                            std::int32_t index) {
-  static constexpr double convertTime = Acts::UnitConstants::s / CLHEP::s;
+std::array<Acts::Vector4, 4u> kinematicsOfStep(const G4Step* step) {
   static constexpr double convertLength = Acts::UnitConstants::mm / CLHEP::mm;
   static constexpr double convertEnergy = Acts::UnitConstants::GeV / CLHEP::GeV;
+  static constexpr double convertTime = Acts::UnitConstants::ns / CLHEP::ns;
+
+  const G4StepPoint* preStepPoint = step->GetPreStepPoint();
+  const G4StepPoint* postStepPoint = step->GetPostStepPoint();
+
+  Acts::Vector4 preStepPosition(convertLength * preStepPoint->GetPosition().x(),
+                                convertLength * preStepPoint->GetPosition().y(),
+                                convertLength * preStepPoint->GetPosition().z(),
+                                convertTime * preStepPoint->GetGlobalTime());
+  Acts::Vector4 preStepMomentum(convertEnergy * preStepPoint->GetMomentum().x(),
+                                convertEnergy * preStepPoint->GetMomentum().y(),
+                                convertEnergy * preStepPoint->GetMomentum().z(),
+                                convertEnergy * preStepPoint->GetTotalEnergy());
+  Acts::Vector4 postStepPosition(
+      convertLength * postStepPoint->GetPosition().x(),
+      convertLength * postStepPoint->GetPosition().y(),
+      convertLength * postStepPoint->GetPosition().z(),
+      convertTime * postStepPoint->GetGlobalTime());
+  Acts::Vector4 postStepMomentum(
+      convertEnergy * postStepPoint->GetMomentum().x(),
+      convertEnergy * postStepPoint->GetMomentum().y(),
+      convertEnergy * postStepPoint->GetMomentum().z(),
+      convertEnergy * postStepPoint->GetTotalEnergy());
+
+  return {preStepPosition, preStepMomentum, postStepPosition, postStepMomentum};
+}
 
-  G4ThreeVector preStepPosition = convertLength * preStepPoint->GetPosition();
-  G4double preStepTime = convertTime * preStepPoint->GetGlobalTime();
-  G4ThreeVector postStepPosition = convertLength * postStepPoint->GetPosition();
-  G4double postStepTime = convertTime * postStepPoint->GetGlobalTime();
-
-  G4ThreeVector preStepMomentum = convertEnergy * preStepPoint->GetMomentum();
-  G4double preStepEnergy = convertEnergy * preStepPoint->GetTotalEnergy();
-  G4ThreeVector postStepMomentum = convertEnergy * postStepPoint->GetMomentum();
-  G4double postStepEnergy = convertEnergy * postStepPoint->GetTotalEnergy();
-
-  Acts::ActsScalar hX = 0.5 * (preStepPosition[0] + postStepPosition[0]);
-  Acts::ActsScalar hY = 0.5 * (preStepPosition[1] + postStepPosition[1]);
-  Acts::ActsScalar hZ = 0.5 * (preStepPosition[2] + postStepPosition[2]);
-  Acts::ActsScalar hT = 0.5 * (preStepTime + postStepTime);
-
-  Acts::ActsScalar mXpre = preStepMomentum[0];
-  Acts::ActsScalar mYpre = preStepMomentum[1];
-  Acts::ActsScalar mZpre = preStepMomentum[2];
-  Acts::ActsScalar mEpre = preStepEnergy;
-  Acts::ActsScalar mXpost = postStepMomentum[0];
-  Acts::ActsScalar mYpost = postStepMomentum[1];
-  Acts::ActsScalar mZpost = postStepMomentum[2];
-  Acts::ActsScalar mEpost = postStepEnergy;
-
-  Acts::Vector4 particlePosition(hX, hY, hZ, hT);
-  Acts::Vector4 beforeMomentum(mXpre, mYpre, mZpre, mEpre);
-  Acts::Vector4 afterMomentum(mXpost, mYpost, mZpost, mEpost);
-
-  return ActsFatras::Hit(geoId, particleId, particlePosition, beforeMomentum,
-                         afterMomentum, index);
+ActsFatras::Hit hitFromStep(const G4Step* step, ActsFatras::Barcode particleId,
+                            Acts::GeometryIdentifier geoId,
+                            std::int32_t index) {
+  auto [preStepPosition, preStepMomentum, postStepPosition, postStepMomentum] =
+      kinematicsOfStep(step);
+
+  return ActsFatras::Hit(geoId, particleId,
+                         0.5 * (preStepPosition + postStepPosition),
+                         preStepMomentum, postStepMomentum, index);
 }
+
+Acts::detail::Step stepFromG4Step(const G4Step* step) {
+  Acts::detail::Step pStep;
+  auto [preStepPosition, preStepMomentum, postStepPosition, postStepMomentum] =
+      kinematicsOfStep(step);
+
+  pStep.navDir = Acts::Direction::Forward;
+  pStep.position = 0.5 * (preStepPosition + postStepPosition).block<3, 1>(0, 0);
+  pStep.momentum = 0.5 * (preStepMomentum + postStepMomentum).block<3, 1>(0, 0);
+  pStep.nTotalTrials = 1;
+  return pStep;
+}
+
 }  // namespace
 
 ActsExamples::SensitiveSteppingAction::SensitiveSteppingAction(
@@ -109,6 +124,10 @@ void ActsExamples::SensitiveSteppingAction::UserSteppingAction(
   G4PrimaryParticle* primaryParticle =
       track->GetDynamicParticle()->GetPrimaryParticle();
 
+  // Get PreStepPoint and PostStepPoint
+  const G4StepPoint* preStepPoint = step->GetPreStepPoint();
+  const G4StepPoint* postStepPoint = step->GetPostStepPoint();
+
   // Bail out if charged & configured to do so
   G4double absCharge = std::abs(track->GetParticleDefinition()->GetPDGCharge());
   if (!m_cfg.charged && absCharge > 0.) {
@@ -138,14 +157,11 @@ void ActsExamples::SensitiveSteppingAction::UserSteppingAction(
   std::string volumeName = volume->GetName();
 
   if (volumeName.find(SensitiveSurfaceMapper::mappingPrefix) ==
-      std::string_view::npos) {
+          std::string::npos &&
+      !m_cfg.stepLogging) {
     return;
   }
 
-  // Get PreStepPoint and PostStepPoint
-  const G4StepPoint* preStepPoint = step->GetPreStepPoint();
-  const G4StepPoint* postStepPoint = step->GetPostStepPoint();
-
   // The G4Touchable for the matching
   const G4VTouchable* touchable = track->GetTouchable();
 
@@ -159,7 +175,8 @@ void ActsExamples::SensitiveSteppingAction::UserSteppingAction(
   ACTS_VERBOSE("Found " << nSurfaces << " candidate surfaces for volume "
                         << volumeName);
 
-  if (nSurfaces == 0) {
+  const Acts::Surface* surface = nullptr;
+  if (nSurfaces == 0 && !m_cfg.stepLogging) {
     ACTS_ERROR("No candidate surfaces found for volume " << volumeName);
     return;
   } else if (nSurfaces == 1u) {
@@ -169,7 +186,7 @@ void ActsExamples::SensitiveSteppingAction::UserSteppingAction(
     // Find the closest surface to the current position
     Acts::GeometryContext gctx;
     for (; bsf != esf; ++bsf) {
-      const Acts::Surface* surface = bsf->second;
+      surface = bsf->second;
       const G4ThreeVector& translation = touchable->GetTranslation();
       Acts::Vector3 g4VolumePosition(convertLength * translation.x(),
                                      convertLength * translation.y(),
@@ -187,9 +204,59 @@ void ActsExamples::SensitiveSteppingAction::UserSteppingAction(
     return;
   }
 
+  // Output is only strictly valid if step logging is not enabled
   const auto particleId = eventStore().trackIdMapping.at(track->GetTrackID());
-
-  ACTS_VERBOSE("Step of " << particleId << " in sensitive volume " << geoId);
+  if (!m_cfg.stepLogging && surface != nullptr) {
+    ACTS_VERBOSE("Step of " << particleId << " in sensitive volume " << geoId);
+  } else if (m_cfg.stepLogging) {
+    if (!eventStore().propagationRecords.contains(track->GetTrackID())) {
+      // Create the propagation summary
+      Acts::ActsScalar xVtx = track->GetVertexPosition().x();
+      Acts::ActsScalar yVtx = track->GetVertexPosition().y();
+      Acts::ActsScalar zVtx = track->GetVertexPosition().z();
+      Acts::ActsScalar xDirVtx = track->GetVertexMomentumDirection().x();
+      Acts::ActsScalar yDirVtx = track->GetVertexMomentumDirection().y();
+      Acts::ActsScalar zDirVtx = track->GetVertexMomentumDirection().z();
+      Acts::ActsScalar absMomentum = track->GetMomentum().mag();
+
+      PropagationSummary iSummary(Acts::CurvilinearTrackParameters(
+          Acts::Vector4(xVtx, yVtx, zVtx, 0.),
+          Acts::Vector3(xDirVtx, yDirVtx, zDirVtx), absCharge / absMomentum,
+          std::nullopt, Acts::ParticleHypothesis::pion()));
+
+      eventStore().propagationRecords.insert({track->GetTrackID(), iSummary});
+    }
+    PropagationSummary& pSummary =
+        eventStore().propagationRecords.at(track->GetTrackID());
+
+    // Increase the step counter
+    pSummary.nSteps += 1;
+
+    Acts::ActsScalar currentTrackLength =
+        track->GetTrackLength() * convertLength;
+    Acts::ActsScalar currentStepLength =
+        currentTrackLength - pSummary.pathLength;
+    pSummary.pathLength = currentTrackLength;
+
+    // Create a new step for the step logging
+    Acts::detail::Step pStep = stepFromG4Step(step);
+    pStep.geoID = geoId;
+    pStep.surface = surface != nullptr ? surface->getSharedPtr() : nullptr;
+    // Check if last step was on same surface
+    if (!pSummary.steps.empty() && pSummary.steps.back().geoID == geoId &&
+        pSummary.steps.back().surface != nullptr) {
+      auto& lastStep = pSummary.steps.back();
+      lastStep.stepSize = Acts::ConstrainedStep(currentStepLength);
+      lastStep.position = 0.5 * (pStep.position + lastStep.position);
+      lastStep.momentum = 0.5 * (pStep.momentum + lastStep.momentum);
+    } else {
+      // Record the propagation state
+      pStep.stepSize = Acts::ConstrainedStep(currentStepLength);
+      pSummary.steps.emplace_back(std::move(pStep));
+    }
+    // You have nothing to do from here
+    return;
+  }
 
   // Set particle hit count to zero, so we have this entry in the map later
   if (!eventStore().particleHitCount.contains(particleId)) {
@@ -230,7 +297,7 @@ void ActsExamples::SensitiveSteppingAction::UserSteppingAction(
     ACTS_VERBOSE("-> merge single step to hit");
     ++eventStore().particleHitCount[particleId];
     eventStore().hits.push_back(
-        hitFromStep(preStepPoint, postStepPoint, particleId, geoId,
+        hitFromStep(step, particleId, geoId,
                     eventStore().particleHitCount.at(particleId) - 1));
 
     eventStore().numberGeantSteps += 1ul;
@@ -243,8 +310,7 @@ void ActsExamples::SensitiveSteppingAction::UserSteppingAction(
   if (postOnBoundary || particleStopped || particleDecayed) {
     ACTS_VERBOSE("-> merge buffer to hit");
     auto& buffer = eventStore().hitBuffer;
-    buffer.push_back(
-        hitFromStep(preStepPoint, postStepPoint, particleId, geoId, -1));
+    buffer.push_back(hitFromStep(step, particleId, geoId, -1));
 
     const auto pos4 =
         0.5 * (buffer.front().fourPosition() + buffer.back().fourPosition());
@@ -272,8 +338,7 @@ void ActsExamples::SensitiveSteppingAction::UserSteppingAction(
   // hit buffer.
   if (!postOnBoundary) {
     // ACTS_VERBOSE("-> add hit to buffer");
-    eventStore().hitBuffer.push_back(
-        hitFromStep(preStepPoint, postStepPoint, particleId, geoId, -1));
+    eventStore().hitBuffer.push_back(hitFromStep(step, particleId, geoId, -1));
     return;
   }
 

Examples/Algorithms/Propagation/CMakeLists.txt

@@ -1,9 +1,4 @@
-add_library(
-    ActsExamplesPropagation
-    SHARED
-    src/PropagationAlgorithm.cpp
-    src/SimHitToSummaryConversion.cpp
-)
+add_library(ActsExamplesPropagation SHARED src/PropagationAlgorithm.cpp)
 
 target_include_directories(
     ActsExamplesPropagation