Geant4/src/SensitiveSteppingAction.cpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 #include "ActsExamples/Geant4/SensitiveSteppingAction.hpp"
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Geometry/GeometryIdentifier.hpp"
0013 #include "Acts/Propagator/detail/SteppingLogger.hpp"
0014 #include "Acts/Surfaces/Surface.hpp"
0015 #include "ActsExamples/EventData/SimParticle.hpp"
0016 #include "ActsExamples/Geant4/AlgebraConverters.hpp"
0017 #include "ActsExamples/Geant4/EventStore.hpp"
0018 #include "ActsExamples/Geant4/SensitiveSurfaceMapper.hpp"
0019 #include "ActsExamples/Geant4/UnitConversion.hpp"
0020 #include "ActsFatras/EventData/Barcode.hpp"
0021
0022 #include <algorithm>
0023 #include <array>
0024 #include <cstddef>
0025 #include <string>
0026 #include <unordered_map>
0027 #include <utility>
0028
0029 #include <G4RunManager.hh>
0030 #include <G4Step.hh>
0031 #include <G4StepPoint.hh>
0032 #include <G4Track.hh>
0033 #include <G4UnitsTable.hh>
0034 #include <G4VPhysicalVolume.hh>
0035 #include <G4VTouchable.hh>
0036 #include <boost/describe.hpp>
0037
0038 BOOST_DESCRIBE_ENUM(G4StepStatus, fWorldBoundary, fGeomBoundary,
0039                     fAtRestDoItProc, fAlongStepDoItProc, fPostStepDoItProc,
0040                     fUserDefinedLimit, fExclusivelyForcedProc, fUndefined);
0041
0042 BOOST_DESCRIBE_ENUM(G4ProcessType, fNotDefined, fTransportation,
0043                     fElectromagnetic, fOptical, fHadronic, fPhotolepton_hadron,
0044                     fDecay, fGeneral, fParameterisation, fUserDefined,
0045                     fParallel, fPhonon, fUCN);
0046
0047 BOOST_DESCRIBE_ENUM(G4TrackStatus, fAlive, fStopButAlive, fStopAndKill,
0048                     fKillTrackAndSecondaries, fSuspend, fPostponeToNextEvent);
0049
0050 namespace ActsExamples::Geant4 {
0051
0052 namespace {
0053
0054 std::array<Acts::Vector4, 4u> kinematicsOfStep(const G4Step& step) {
0055   assert(step.GetPreStepPoint() != nullptr);
0056   assert(step.GetPostStepPoint() != nullptr);
0057   const G4StepPoint& preStepPoint = *step.GetPreStepPoint();
0058   const G4StepPoint& postStepPoint = *step.GetPostStepPoint();
0059
0060   const Acts::Vector4 preStepPosition =
0061       convertPosition(preStepPoint.GetPosition(), preStepPoint.GetGlobalTime());
0062
0063   const Acts::Vector4 preStepMomentum = convertMomentum(
0064       preStepPoint.GetMomentum(), preStepPoint.GetTotalEnergy());
0065   const Acts::Vector4 postStepPosition = convertPosition(
0066       postStepPoint.GetPosition(), postStepPoint.GetGlobalTime());
0067
0068   const Acts::Vector4 postStepMomentum = convertMomentum(
0069       postStepPoint.GetMomentum(), postStepPoint.GetTotalEnergy());
0070
0071   return {preStepPosition, preStepMomentum, postStepPosition, postStepMomentum};
0072 }
0073
0074 ActsFatras::Hit hitFromStep(const G4Step& step, SimBarcode particleId,
0075                             Acts::GeometryIdentifier geoId,
0076                             std::int32_t index) {
0077   const auto [preStepPosition, preStepMomentum, postStepPosition,
0078               postStepMomentum] = kinematicsOfStep(step);
0079
0080   return ActsFatras::Hit(geoId, particleId,
0081                          0.5 * (preStepPosition + postStepPosition),
0082                          preStepMomentum, postStepMomentum, index);
0083 }
0084
0085 Acts::detail::Step stepFromG4Step(const G4Step& step) {
0086   Acts::detail::Step pStep;
0087   const auto [preStepPosition, preStepMomentum, postStepPosition,
0088               postStepMomentum] = kinematicsOfStep(step);
0089
0090   pStep.navDir = Acts::Direction::Forward();
0091   pStep.position = 0.5 * (preStepPosition + postStepPosition).block<3, 1>(0, 0);
0092   pStep.momentum = 0.5 * (preStepMomentum + postStepMomentum).block<3, 1>(0, 0);
0093   pStep.nTotalTrials = 1;
0094   return pStep;
0095 }
0096
0097 }  // namespace
0098
0099 SensitiveSteppingAction::SensitiveSteppingAction(
0100     const Config& cfg, std::unique_ptr<const Acts::Logger> logger)
0101     : G4UserSteppingAction(), m_cfg(cfg), m_logger(std::move(logger)) {}
0102
0103 void SensitiveSteppingAction::UserSteppingAction(const G4Step* stepPtr) {
0104   assert(stepPtr != nullptr);
0105   const G4Step& step = *stepPtr;
0106
0107   // The particle after the step
0108   assert(step.GetTrack() != nullptr);
0109   const G4Track& track = *step.GetTrack();
0110   G4PrimaryParticle* primaryParticle =
0111       track.GetDynamicParticle()->GetPrimaryParticle();
0112
0113   // Get PreStepPoint and PostStepPoint
0114   assert(step.GetPreStepPoint() != nullptr);
0115   assert(step.GetPostStepPoint() != nullptr);
0116   const G4StepPoint& preStepPoint = *step.GetPreStepPoint();
0117   const G4StepPoint& postStepPoint = *step.GetPostStepPoint();
0118
0119   // Bail out if charged & configured to do so
0120   const double absCharge =
0121       std::abs(track.GetParticleDefinition()->GetPDGCharge());
0122   if (!m_cfg.charged && absCharge > 0.) {
0123     return;
0124   }
0125
0126   // Bail out if neutral & configured to do so
0127   if (!m_cfg.neutral && absCharge == 0.) {
0128     return;
0129   }
0130
0131   // Bail out if it is a primary & configured to be ignored
0132   if (!m_cfg.primary && primaryParticle != nullptr) {
0133     return;
0134   }
0135
0136   // Bail out if it is a secondary & configured to be ignored
0137   if (!m_cfg.secondary && primaryParticle == nullptr) {
0138     return;
0139   }
0140
0141   // Get the physical volume & check if it has the sensitive string name
0142   const G4VPhysicalVolume* volume = track.GetVolume();
0143   if (volume == nullptr) {
0144     throw std::runtime_error("No volume found, terminate simulation");
0145   }
0146   const std::string volumeName = volume->GetName();
0147   ACTS_VERBOSE("Check whether volume " << volumeName << " is sensitive");
0148   if (!m_cfg.stepLogging &&
0149       volumeName.find(SensitiveSurfaceMapper::mappingPrefix) ==
0150           std::string::npos) {
0151     return;
0152   }
0153
0154   // The G4Touchable for the matching
0155   const G4VTouchable* touchable = track.GetTouchable();
0156
0157   Acts::GeometryIdentifier geoId{};
0158   const Acts::Surface* surface = nullptr;
0159
0160   // Find the range of candidate surfaces for the current position in the map
0161   const auto surfaceMap_itr = m_surfaceMapping.find(volume);
0162   if (surfaceMap_itr != m_surfaceMapping.end()) {
0163     const auto& surfacesToG4Vol = surfaceMap_itr->second;
0164     ACTS_VERBOSE("Found " << surfacesToG4Vol.size()
0165                           << " candidate surfaces for volume " << volumeName);
0166     const Acts::Vector3 volumePos =
0167         convertPosition(touchable->GetTranslation());
0168     const auto lookUp_itr = surfacesToG4Vol.find(volumePos);
0169     if (lookUp_itr == surfacesToG4Vol.end() && !m_cfg.stepLogging) {
0170       ACTS_ERROR("No candidate surfaces found for volume " << volumeName);
0171       return;
0172     }
0173     surface = lookUp_itr->second;
0174     geoId = surface->geometryId();
0175     ACTS_VERBOSE("Replica assignment successful -> to surface " << geoId);
0176   } else if (!m_cfg.stepLogging) {
0177     ACTS_ERROR("No candidate surfaces found for volume " << volumeName);
0178     return;
0179   }
0180   // This is not the case if we have a particle-ID collision
0181   if (!eventStore().trackIdMapping.contains(track.GetTrackID())) {
0182     return;
0183   }
0184
0185   // Output is only strictly valid if step logging is not enabled
0186   const SimBarcode particleId =
0187       eventStore().trackIdMapping.at(track.GetTrackID());
0188   if (!m_cfg.stepLogging && surface != nullptr) {
0189     ACTS_VERBOSE("Step of " << particleId << " in sensitive volume " << geoId);
0190   } else if (m_cfg.stepLogging) {
0191     if (!eventStore().propagationRecords.contains(track.GetTrackID())) {
0192       // Create the propagation summary
0193       const double absMomentum =
0194           track.GetMomentum().mag() * convertEnergyToActs;
0195
0196       PropagationSummary iSummary(Acts::BoundTrackParameters::createCurvilinear(
0197           convertPosition(track.GetVertexPosition(), 0.),
0198           convertDirection(track.GetVertexMomentumDirection()),
0199           absCharge / absMomentum, std::nullopt,
0200           Acts::ParticleHypothesis::pion()));
0201
0202       eventStore().propagationRecords.insert({track.GetTrackID(), iSummary});
0203     }
0204     PropagationSummary& pSummary =
0205         eventStore().propagationRecords.at(track.GetTrackID());
0206
0207     // Increase the step counter
0208     pSummary.nSteps += 1;
0209
0210     const double currentTrackLength =
0211         track.GetTrackLength() * convertLengthToActs;
0212     const double currentStepLength = currentTrackLength - pSummary.pathLength;
0213     pSummary.pathLength = currentTrackLength;
0214
0215     // Create a new step for the step logging
0216     Acts::detail::Step pStep = stepFromG4Step(step);
0217     pStep.geoID = geoId;
0218     pStep.surface = surface != nullptr ? surface->getSharedPtr() : nullptr;
0219     // Check if last step was on same surface
0220     if (!pSummary.steps.empty() && pSummary.steps.back().geoID == geoId &&
0221         pSummary.steps.back().surface != nullptr) {
0222       auto& lastStep = pSummary.steps.back();
0223       lastStep.stepSize = Acts::ConstrainedStep(currentStepLength);
0224       lastStep.position = 0.5 * (pStep.position + lastStep.position);
0225       lastStep.momentum = 0.5 * (pStep.momentum + lastStep.momentum);
0226     } else {
0227       // Record the propagation state
0228       pStep.stepSize = Acts::ConstrainedStep(currentStepLength);
0229       pSummary.steps.emplace_back(std::move(pStep));
0230     }
0231     // You have nothing to do from here
0232     return;
0233   }
0234
0235   // Set particle hit count to zero, so we have this entry in the map later
0236   if (!eventStore().particleHitCount.contains(particleId)) {
0237     eventStore().particleHitCount[particleId] = 0;
0238   }
0239
0240   // Extract if we are at volume boundaries
0241   const bool preOnBoundary = preStepPoint.GetStepStatus() == fGeomBoundary;
0242   const bool postOnBoundary = postStepPoint.GetStepStatus() == fGeomBoundary ||
0243                               postStepPoint.GetStepStatus() == fWorldBoundary;
0244   const bool particleStopped = (postStepPoint.GetKineticEnergy() == 0.0);
0245   const bool particleDecayed =
0246       (postStepPoint.GetProcessDefinedStep()->GetProcessType() == fDecay);
0247
0248   const auto print = [](auto s) {
0249     return boost::describe::enum_to_string(s, "unmatched");
0250   };
0251
0252   ACTS_VERBOSE("status: pre="
0253                << print(preStepPoint.GetStepStatus())
0254                << ", post=" << print(postStepPoint.GetStepStatus())
0255                << ", post E_kin=" << std::boolalpha
0256                << postStepPoint.GetKineticEnergy() << ", process_type="
0257                << print(postStepPoint.GetProcessDefinedStep()->GetProcessType())
0258                << ", particle="
0259                << track.GetParticleDefinition()->GetParticleName()
0260                << ", process_name="
0261                << postStepPoint.GetProcessDefinedStep()->GetProcessName()
0262                << ", track status=" << print(track.GetTrackStatus()));
0263
0264   // Case A: The step starts at the entry of the volume and ends at the exit.
0265   // Add hit to collection.
0266   if (preOnBoundary && postOnBoundary) {
0267     ACTS_VERBOSE("-> merge single step to hit");
0268     ++eventStore().particleHitCount[particleId];
0269     eventStore().hits.push_back(
0270         hitFromStep(step, particleId, geoId,
0271                     eventStore().particleHitCount.at(particleId) - 1));
0272
0273     eventStore().numberGeantSteps += 1ul;
0274     eventStore().maxStepsForHit = std::max(eventStore().maxStepsForHit, 1ul);
0275     return;
0276   }
0277
0278   // Case B: The step ends at the exit of the volume. Add hit to hit buffer, and
0279   // then combine hit buffer
0280   if (postOnBoundary || particleStopped || particleDecayed) {
0281     ACTS_VERBOSE("-> merge buffer to hit");
0282     auto& buffer = eventStore().hitBuffer;
0283     buffer.push_back(hitFromStep(step, particleId, geoId, -1));
0284
0285     const Acts::Vector4 pos4 =
0286         0.5 * (buffer.front().fourPosition() + buffer.back().fourPosition());
0287
0288     ++eventStore().particleHitCount[particleId];
0289     eventStore().hits.emplace_back(
0290         geoId, particleId, pos4, buffer.front().momentum4Before(),
0291         buffer.back().momentum4After(),
0292         eventStore().particleHitCount.at(particleId) - 1);
0293
0294     assert(std::ranges::all_of(
0295         buffer, [&](const auto& h) { return h.geometryId() == geoId; }));
0296     assert(std::ranges::all_of(
0297         buffer, [&](const auto& h) { return h.particleId() == particleId; }));
0298
0299     eventStore().numberGeantSteps += buffer.size();
0300     eventStore().maxStepsForHit =
0301         std::max(eventStore().maxStepsForHit, buffer.size());
0302
0303     buffer.clear();
0304     return;
0305   }
0306
0307   // Case C: The step doesn't end at the exit of the volume. Add the hit to the
0308   // hit buffer.
0309   if (!postOnBoundary) {
0310     // ACTS_VERBOSE("-> add hit to buffer");
0311     eventStore().hitBuffer.push_back(hitFromStep(step, particleId, geoId, -1));
0312     return;
0313   }
0314
0315   assert(false && "should never reach this");
0316 }
0317
0318 }  // namespace ActsExamples::Geant4