feat: Step merger for Geant4 and digitization energy threshold (#2184)

This PR does tow things:

* Add a step merger to the Geant4 sim. This fixes the way how we record hits: Geant4 can do multiple steps through the material, which need to be combined to one single hit. Bevor, we recorded each step as hit.

* Add a energy threshold for the hit acceptance to the digitization: This allows the digitization to reject hits with a very small energy deposit. The default value for this is `1000 * 3.65 eV` (1000 electron-hole pairs for activation).

Together, this reduces the all-over hits, that the `Geant4Simulation` produces on a ttbar event drastically. However, the number of "target hits" (hits of particles with pT > 500 MeV and nhits > 3 in this case) stays roughly the same.

![image](https://github.com/acts-project/acts/assets/37871400/5524cf19-ed0d-4cf1-adab-01478749a862)

The energy threshold seems reasonable: In a particle-gun muon-only run in the ODD we see that almost all muons pass the threshold, but a lot low energy deposit electron hits are rejected by the digitization (energy deposit is negative, since it is stored as energy loss of the particle):

![image](https://github.com/acts-project/acts/assets/37871400/c6ae11a2-af5f-48b1-90cc-e76bffd9558f)

Examples/Algorithms/Digitization/include/ActsExamples/Digitization/DigitizationConfig.hpp

@@ -9,6 +9,7 @@
 #pragma once
 
 #include "Acts/Definitions/Algebra.hpp"
+#include "Acts/Definitions/Units.hpp"
 #include "Acts/Geometry/GeometryHierarchyMap.hpp"
 #include "Acts/Geometry/TrackingGeometry.hpp"
 #include "Acts/Utilities/BinUtility.hpp"
@@ -130,6 +131,11 @@ class DigitizationConfig {
   const double mergeNsigma;
   /// Consider clusters that share a corner as merged (8-cell connectivity)
   const bool mergeCommonCorner;
+  /// Energy deposit threshold for accepting a hit
+  /// For a generic readout frontend we assume 1000 e/h pairs, in Si each
+  /// e/h-pair requiers on average an energy of 3.65 eV (PDG  review 2023,
+  /// Table 35.10)
+  double minEnergyDeposit = 1000 * 3.65 * Acts::UnitConstants::eV;
   /// The digitizers per GeometryIdentifiers
   Acts::GeometryHierarchyMap<DigiComponentsConfig> digitizationConfigs;
 

Examples/Algorithms/Geant4/CMakeLists.txt

@@ -29,7 +29,7 @@ target_include_directories(
   PUBLIC $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/include>)
 target_link_libraries(
   ActsExamplesGeant4
-  PUBLIC ActsCore ActsExamplesFramework ActsExamplesDetectorTelescope ${Geant4_LIBRARIES})
+  PUBLIC ActsCore ActsExamplesFramework ActsExamplesDetectorTelescope Boost::headers ${Geant4_LIBRARIES})
 
 if (ACTS_BUILD_EXAMPLES_DD4HEP)
   if(${DD4hep_VERSION} VERSION_LESS 1.11)

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

@@ -50,6 +50,13 @@ class EventStoreRegistry {
     /// The hits in sensitive detectors
     SimHitContainer::sequence_type hits;
 
+    /// Hit buffer for step merging (multiple steps in senstive volume)
+    std::vector<ActsFatras::Hit> hitBuffer;
+
+    /// Some statistics for the step merging
+    std::size_t numberGeantSteps = 0;
+    std::size_t maxStepsForHit = 0;
+
     /// Tracks recorded in material mapping
     std::unordered_map<size_t, Acts::RecordedMaterialTrack> materialTracks;
 

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

@@ -34,7 +34,7 @@ namespace ActsExamples {
 /// elements of the Acts::TrackingGeoemtry w/o map lookup.
 class SensitiveSurfaceMapper {
  public:
-  static constexpr const char* mappingPrefix = "ActsGeoID#";
+  constexpr static std::string_view mappingPrefix = "ActsGeoID#";
 
   /// Configuration struct for the surface mapper
   struct Config {

Examples/Algorithms/Geant4/src/Geant4Simulation.cpp

@@ -83,11 +83,13 @@ ActsExamples::Geant4Simulation::Geant4Simulation(
       geantVerboseLevel);
 
   // Suppress the printing of physics information.
+  if (logger().level() > Acts::Logging::DEBUG) {
 #if G4VERSION_NUMBER >= 1100
-  G4HadronicParameters::Instance()->SetVerboseLevel(0);
-  G4HadronicProcessStore::Instance()->SetVerbose(0);
-  G4EmParameters::Instance()->SetIsPrintedFlag(true);
+    G4HadronicParameters::Instance()->SetVerboseLevel(geantVerboseLevel);
+    G4HadronicProcessStore::Instance()->SetVerbose(geantVerboseLevel);
+    G4EmParameters::Instance()->SetIsPrintedFlag(true);
 #endif
+  }
 
   // Set the detector construction
   m_cfg.runManager->SetUserInitialization(m_cfg.detectorConstruction);
@@ -210,5 +212,10 @@ ActsExamples::ProcessCode ActsExamples::Geant4Simulation::execute(
         ctx, decltype(eventData.materialTracks)(eventData.materialTracks));
   }
 
+  ACTS_INFO("Step merging: mean hits per hit: "
+            << static_cast<double>(eventData.numberGeantSteps) /
+                   eventData.hits.size());
+  ACTS_INFO("Step merging: max hits per hit: " << eventData.maxStepsForHit);
+
   return ActsExamples::ProcessCode::SUCCESS;
 }

Examples/Algorithms/Geant4/src/ParticleTrackingAction.cpp

@@ -26,6 +26,17 @@ void ActsExamples::ParticleTrackingAction::PreUserTrackingAction(
     const G4Track* aTrack) {
   auto& eventData = EventStoreRegistry::eventData();
 
+  // If this is not the case, there are unhandled cases of particle stopping in
+  // the SensitiveSteppingAction
+  // TODO We could also merge the remaining hits to a hit here, but it would be
+  // nicer to investigate, if we can handle all particle stop conditions in the
+  // SensitiveSteppingAction... This seems to happen O(1) times in a ttbar
+  // event, so seems not to be too problematic
+  if (not eventData.hitBuffer.empty()) {
+    eventData.hitBuffer.clear();
+    ACTS_WARNING("Hit buffer not empty after track");
+  }
+
   auto particleId = makeParticleId(aTrack->GetTrackID(), aTrack->GetParentID());
 
   // There is already a warning printed in the makeParticleId function

Examples/Algorithms/Geant4/src/SensitiveSteppingAction.cpp

@@ -19,16 +19,72 @@
 #include <G4UnitsTable.hh>
 #include <G4VPhysicalVolume.hh>
 
+#if BOOST_VERSION >= 107800
+#include <boost/describe.hpp>
+
+BOOST_DESCRIBE_ENUM(G4StepStatus, fWorldBoundary, fGeomBoundary,
+                    fAtRestDoItProc, fAlongStepDoItProc, fPostStepDoItProc,
+                    fUserDefinedLimit, fExclusivelyForcedProc, fUndefined);
+
+BOOST_DESCRIBE_ENUM(G4ProcessType, fNotDefined, fTransportation,
+                    fElectromagnetic, fOptical, fHadronic, fPhotolepton_hadron,
+                    fDecay, fGeneral, fParameterisation, fUserDefined,
+                    fParallel, fPhonon, fUCN);
+
+BOOST_DESCRIBE_ENUM(G4TrackStatus, fAlive, fStopButAlive, fStopAndKill,
+                    fKillTrackAndSecondaries, fSuspend, fPostponeToNextEvent);
+#endif
+
+namespace {
+
+ActsFatras::Hit hitFromStep(const G4StepPoint* preStepPoint,
+                            const G4StepPoint* postStepPoint,
+                            ActsFatras::Barcode particleId,
+                            Acts::GeometryIdentifier geoId, int32_t index) {
+  static constexpr double convertTime = Acts::UnitConstants::s / CLHEP::s;
+  static constexpr double convertLength = Acts::UnitConstants::mm / CLHEP::mm;
+  static constexpr double convertEnergy = Acts::UnitConstants::GeV / CLHEP::GeV;
+
+  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);
+}
+}  // namespace
+
 ActsExamples::SensitiveSteppingAction::SensitiveSteppingAction(
     const Config& cfg, std::unique_ptr<const Acts::Logger> logger)
     : G4UserSteppingAction(), m_cfg(cfg), m_logger(std::move(logger)) {}
 
 void ActsExamples::SensitiveSteppingAction::UserSteppingAction(
     const G4Step* step) {
   // Unit conversions G4->::ACTS
-  constexpr double convertTime = Acts::UnitConstants::s / CLHEP::s;
-  constexpr double convertLength = Acts::UnitConstants::mm / CLHEP::mm;
-  constexpr double convertEnergy = Acts::UnitConstants::GeV / CLHEP::GeV;
 
   // Retrieve the event data registry
   auto& eventData = EventStoreRegistry::eventData();
@@ -60,68 +116,121 @@ void ActsExamples::SensitiveSteppingAction::UserSteppingAction(
   }
 
   // Get the physical volume & check if it has the sensitive string name
-  G4VPhysicalVolume* volume = track->GetVolume();
-  std::string volumeName = volume->GetName();
-
-  std::string mappingPfx(SensitiveSurfaceMapper::mappingPrefix);
+  std::string_view volumeName = track->GetVolume()->GetName();
 
-  if (volumeName.find(mappingPfx) == std::string::npos) {
+  if (volumeName.find(SensitiveSurfaceMapper::mappingPrefix) ==
+      std::string_view::npos) {
     return;
   }
 
-  ACTS_VERBOSE("Step in senstive volume " << volumeName);
-
-  // Get PreStepPoint and PostStepPoint
-  G4StepPoint* preStepPoint = step->GetPreStepPoint();
-  G4StepPoint* postStepPoint = step->GetPostStepPoint();
-
-  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;
-
   // Cast out the GeometryIdentifier
-  volumeName.erase(0, mappingPfx.size());
-
-  Acts::GeometryIdentifier::Value sGeoVal = std::stoul(volumeName);
-  Acts::GeometryIdentifier geoID(sGeoVal);
+  volumeName = volumeName.substr(SensitiveSurfaceMapper::mappingPrefix.size(),
+                                 std::string_view::npos);
+  char* end = nullptr;
+  const Acts::GeometryIdentifier geoId(
+      std::strtoul(volumeName.data(), &end, 10));
 
   // This is not the case if we have a particle-ID collision
   if (eventData.trackIdMapping.find(track->GetTrackID()) ==
       eventData.trackIdMapping.end()) {
     return;
   }
 
-  auto particleID = eventData.trackIdMapping.at(track->GetTrackID());
+  const auto particleId = eventData.trackIdMapping.at(track->GetTrackID());
 
-  Acts::Vector4 particlePosition(hX, hY, hZ, hT);
-  Acts::Vector4 beforeMomentum(mXpre, mYpre, mZpre, mEpre);
-  Acts::Vector4 afterMomentum(mXpost, mYpost, mZpost, mEpost);
+  ACTS_VERBOSE("Step of " << particleId << " in senstive volume " << geoId);
+
+  // Get PreStepPoint and PostStepPoint
+  const G4StepPoint* preStepPoint = step->GetPreStepPoint();
+  const G4StepPoint* postStepPoint = step->GetPostStepPoint();
 
-  // Increase counter (starts at 1 because of ++)
-  ++eventData.particleHitCount[particleID];
+  // Set particle hit count to zero, so we have this entry in the map later
+  if (eventData.particleHitCount.find(particleId) ==
+      eventData.particleHitCount.end()) {
+    eventData.particleHitCount[particleId] = 0;
+  }
+
+  // Extract if we are at volume boundaries
+  const bool preOnBoundary = preStepPoint->GetStepStatus() == fGeomBoundary;
+  const bool postOnBoundary = postStepPoint->GetStepStatus() == fGeomBoundary or
+                              postStepPoint->GetStepStatus() == fWorldBoundary;
+  const bool particleStopped = (postStepPoint->GetKineticEnergy() == 0.0);
+  const bool particleDecayed =
+      (postStepPoint->GetProcessDefinedStep()->GetProcessType() == fDecay);
+
+  auto print = [](auto s) {
+#if BOOST_VERSION >= 107800
+    return boost::describe::enum_to_string(s, "unmatched");
+#else
+    return s;
+#endif
+  };
+  ACTS_VERBOSE("status: pre="
+               << print(preStepPoint->GetStepStatus())
+               << ", post=" << print(postStepPoint->GetStepStatus())
+               << ", post E_kin=" << std::boolalpha
+               << postStepPoint->GetKineticEnergy() << ", process_type="
+               << print(
+                      postStepPoint->GetProcessDefinedStep()->GetProcessType())
+               << ", particle="
+               << track->GetParticleDefinition()->GetParticleName()
+               << ", process_name="
+               << postStepPoint->GetProcessDefinedStep()->GetProcessName()
+               << ", track status=" << print(track->GetTrackStatus()));
+
+  // Case A: The step starts at the entry of the volume and ends at the exit.
+  // Add hit to collection.
+  if (preOnBoundary and postOnBoundary) {
+    ACTS_VERBOSE("-> merge single step to hit");
+    ++eventData.particleHitCount[particleId];
+    eventData.hits.push_back(
+        hitFromStep(preStepPoint, postStepPoint, particleId, geoId,
+                    eventData.particleHitCount.at(particleId) - 1));
+
+    eventData.numberGeantSteps += 1ul;
+    eventData.maxStepsForHit = std::max(eventData.maxStepsForHit, 1ul);
+    return;
+  }
+
+  // Case B: The step ends at the exit of the volume. Add hit to hit buffer, and
+  // then combine hit buffer
+  if (postOnBoundary or particleStopped or particleDecayed) {
+    ACTS_VERBOSE("-> merge buffer to hit");
+    auto& buffer = eventData.hitBuffer;
+    buffer.push_back(
+        hitFromStep(preStepPoint, postStepPoint, particleId, geoId, -1));
+
+    const auto pos4 =
+        0.5 * (buffer.front().fourPosition() + buffer.back().fourPosition());
+
+    ++eventData.particleHitCount[particleId];
+    eventData.hits.emplace_back(geoId, particleId, pos4,
+                                buffer.front().momentum4Before(),
+                                buffer.back().momentum4After(),
+                                eventData.particleHitCount.at(particleId) - 1);
+
+    assert(std::all_of(buffer.begin(), buffer.end(),
+                       [&](const auto& h) { return h.geometryId() == geoId; }));
+    assert(std::all_of(buffer.begin(), buffer.end(), [&](const auto& h) {
+      return h.particleId() == particleId;
+    }));
+
+    eventData.numberGeantSteps += buffer.size();
+    eventData.maxStepsForHit =
+        std::max(eventData.maxStepsForHit, buffer.size());
+
+    buffer.clear();
+    return;
+  }
+
+  // Case C: The step doesn't end at the exit of the volume. Add the hit to the
+  // hit buffer.
+  if (not postOnBoundary) {
+    // ACTS_VERBOSE("-> add hit to buffer");
+    eventData.hitBuffer.push_back(
+        hitFromStep(preStepPoint, postStepPoint, particleId, geoId, -1));
+    return;
+  }
 
-  // Fill into the registry (subtract 1 from hit-count to get 0-based index)
-  eventData.hits.emplace_back(geoID, particleID, particlePosition,
-                              beforeMomentum, afterMomentum,
-                              eventData.particleHitCount.at(particleID) - 1);
+  assert(false && "should never reach this");
 }

Examples/Algorithms/Geant4/src/SensitiveSurfaceMapper.cpp

@@ -105,7 +105,7 @@ void ActsExamples::SensitiveSurfaceMapper::remapSensitiveNames(
     // contains the GeometryID/
     if (mappedSurface != nullptr) {
       ++sCounter;
-      std::string mappedVolumeName = SensitiveSurfaceMapper::mappingPrefix;
+      std::string mappedVolumeName(SensitiveSurfaceMapper::mappingPrefix);
       mappedVolumeName += std::to_string(mappedSurface->geometryId().value());
       ACTS_VERBOSE("Found matching surface " << mappedSurface->geometryId()
                                              << " at position "

Examples/Python/python/acts/examples/simulation.py

@@ -699,6 +699,7 @@ def addDigitization(
     outputDirRoot: Optional[Union[Path, str]] = None,
     rnd: Optional[acts.examples.RandomNumbers] = None,
     doMerge: Optional[bool] = None,
+    minEnergyDeposit: Optional[float] = None,
     logLevel: Optional[acts.logging.Level] = None,
 ) -> acts.examples.Sequencer:
     """This function steers the digitization step
@@ -738,6 +739,11 @@ def addDigitization(
         outputMeasurementSimHitsMap="measurement_simhits_map",
         doMerge=doMerge,
     )
+
+    # Not sure how to do this in our style
+    if minEnergyDeposit is not None:
+        digiCfg.minEnergyDeposit = minEnergyDeposit
+
     digiAlg = acts.examples.DigitizationAlgorithm(digiCfg, customLogLevel())
 
     s.addAlgorithm(digiAlg)

Examples/Python/src/Digitization.cpp

@@ -57,6 +57,7 @@ void addDigitization(Context& ctx) {
     ACTS_PYTHON_MEMBER(trackingGeometry);
     ACTS_PYTHON_MEMBER(randomNumbers);
     ACTS_PYTHON_MEMBER(doMerge);
+    ACTS_PYTHON_MEMBER(minEnergyDeposit);
     ACTS_PYTHON_MEMBER(digitizationConfigs);
     ACTS_PYTHON_STRUCT_END();
 

Examples/Python/tests/root_file_hashes.txt

@@ -4,7 +4,7 @@ test_fatras__particles_initial.root: 0e2d9974fdd4aa5549c5c05817ec54fe6c7d9cbe796
 test_fatras__hits.root: d46b760e1647ab94c49d6976a78ee5a6178d41c45f2238f86f3987938a7901fb
 test_geant4__particles_final.root: 94069908d2c77d5dd54193fbaa0d22d401b98df17c148fff6b48dbfa9e0b8a55
 test_geant4__particles_initial.root: 0c70278caa0adebcae8ce64e8f262b99daab0c905bc4aefb9600bcb22106e7c6
-test_geant4__hits.root: 0ec94c05c6d657fbf56ac9407a66c6a19514c3cd421cd47de905952824902575
+test_geant4__hits.root: 70b592a546fd362c9341d87c9068b1f8fed657e0036e62e7262c01fa5eb2d469
 test_seeding__estimatedparams.root: fe5ce98fb6b6886e05f5e6c17e037c898017a48201534b07bc56ccbf672c7a33
 test_seeding__performance_seeding.root: 992f9c611d30dde0d3f3ab676bab19ada61ab6a4442828e27b65ec5e5b7a2880
 test_seeding__particles.root: 74e08ee12bdaf9f7d273369c71c742df345475c8b73bb7ce223619d0dd1f87ee