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/Definitions/Units.hpp"
0013 #include "Acts/Geometry/GeometryIdentifier.hpp"
0014 #include "Acts/Propagator/detail/SteppingLogger.hpp"
0015 #include "Acts/Surfaces/Surface.hpp"
0016 #include "Acts/Utilities/MultiIndex.hpp"
0017 #include "ActsExamples/Geant4/AlgebraConverters.hpp"
0018 #include "ActsExamples/Geant4/EventStore.hpp"
0019 #include "ActsExamples/Geant4/SensitiveSurfaceMapper.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/version.hpp>
0037
0038 #if BOOST_VERSION >= 107800
0039 #include <boost/describe.hpp>
0040
0041 BOOST_DESCRIBE_ENUM(G4StepStatus, fWorldBoundary, fGeomBoundary,
0042                     fAtRestDoItProc, fAlongStepDoItProc, fPostStepDoItProc,
0043                     fUserDefinedLimit, fExclusivelyForcedProc, fUndefined);
0044
0045 BOOST_DESCRIBE_ENUM(G4ProcessType, fNotDefined, fTransportation,
0046                     fElectromagnetic, fOptical, fHadronic, fPhotolepton_hadron,
0047                     fDecay, fGeneral, fParameterisation, fUserDefined,
0048                     fParallel, fPhonon, fUCN);
0049
0050 BOOST_DESCRIBE_ENUM(G4TrackStatus, fAlive, fStopButAlive, fStopAndKill,
0051                     fKillTrackAndSecondaries, fSuspend, fPostponeToNextEvent);
0052 #endif
0053
0054 namespace {
0055
0056 std::array<Acts::Vector4, 4u> kinematicsOfStep(const G4Step* step) {
0057   const G4StepPoint* preStepPoint = step->GetPreStepPoint();
0058   const G4StepPoint* postStepPoint = step->GetPostStepPoint();
0059   using namespace ActsExamples::Geant4;
0060   Acts::Vector4 preStepPosition = convertPosition(
0061       preStepPoint->GetPosition(), preStepPoint->GetGlobalTime());
0062
0063   Acts::Vector4 preStepMomentum = convertMomentum(
0064       preStepPoint->GetMomentum(), preStepPoint->GetTotalEnergy());
0065   Acts::Vector4 postStepPosition = convertPosition(
0066       postStepPoint->GetPosition(), postStepPoint->GetGlobalTime());
0067
0068   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, ActsFatras::Barcode particleId,
0075                             Acts::GeometryIdentifier geoId,
0076                             std::int32_t index) {
0077   auto [preStepPosition, preStepMomentum, postStepPosition, postStepMomentum] =
0078       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   auto [preStepPosition, preStepMomentum, postStepPosition, postStepMomentum] =
0088       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 namespace ActsExamples::Geant4 {
0100
0101 SensitiveSteppingAction::SensitiveSteppingAction(
0102     const Config& cfg, std::unique_ptr<const Acts::Logger> logger)
0103     : G4UserSteppingAction(), m_cfg(cfg), m_logger(std::move(logger)) {}
0104
0105 void SensitiveSteppingAction::UserSteppingAction(const G4Step* step) {
0106   // Unit conversions G4->::ACTS
0107   static constexpr double convertLength = Acts::UnitConstants::mm / CLHEP::mm;
0108   static constexpr double convertEnergy = Acts::UnitConstants::GeV / CLHEP::GeV;
0109   static constexpr auto mappingPrefix = SensitiveSurfaceMapper::mappingPrefix;
0110
0111   // The particle after the step
0112   G4Track* track = step->GetTrack();
0113   G4PrimaryParticle* primaryParticle =
0114       track->GetDynamicParticle()->GetPrimaryParticle();
0115
0116   // Get PreStepPoint and PostStepPoint
0117   const G4StepPoint* preStepPoint = step->GetPreStepPoint();
0118   const G4StepPoint* postStepPoint = step->GetPostStepPoint();
0119
0120   // Bail out if charged & configured to do so
0121   G4double absCharge = 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   std::string volumeName = volume->GetName();
0147   ACTS_VERBOSE("Check whether volume " << volumeName << " is sensitive");
0148   if (!m_cfg.stepLogging &&
0149       volumeName.find(mappingPrefix) == std::string::npos) {
0150     return;
0151   }
0152
0153   // The G4Touchable for the matching
0154   const G4VTouchable* touchable = track->GetTouchable();
0155
0156   Acts::GeometryIdentifier geoId{};
0157   const Acts::Surface* surface = nullptr;
0158
0159   // Find the range of candidate surfaces for the current position in the map
0160   const auto surfaceMap_itr = m_surfaceMapping.find(volume);
0161   if (surfaceMap_itr != m_surfaceMapping.end()) {
0162     const auto& surfacesToG4Vol = surfaceMap_itr->second;
0163     ACTS_VERBOSE("Found " << surfacesToG4Vol.size()
0164                           << " candidate surfaces for volume " << volumeName);
0165     const Acts::Vector3 volumePos =
0166         convertPosition(touchable->GetTranslation());
0167     const auto lookUp_itr = surfacesToG4Vol.find(volumePos);
0168     if (lookUp_itr == surfacesToG4Vol.end() && !m_cfg.stepLogging) {
0169       ACTS_ERROR("No candidate surfaces found for volume " << volumeName);
0170       return;
0171     }
0172     surface = lookUp_itr->second;
0173     geoId = surface->geometryId();
0174     ACTS_VERBOSE("Replica assignment successful -> to surface " << geoId);
0175   } else if (!m_cfg.stepLogging) {
0176     ACTS_ERROR("No candidate surfaces found for volume " << volumeName);
0177     return;
0178   }
0179   // This is not the case if we have a particle-ID collision
0180   if (!eventStore().trackIdMapping.contains(track->GetTrackID())) {
0181     return;
0182   }
0183
0184   // Output is only strictly valid if step logging is not enabled
0185   const auto particleId = eventStore().trackIdMapping.at(track->GetTrackID());
0186   if (!m_cfg.stepLogging && surface != nullptr) {
0187     ACTS_VERBOSE("Step of " << particleId << " in sensitive volume " << geoId);
0188   } else if (m_cfg.stepLogging) {
0189     if (!eventStore().propagationRecords.contains(track->GetTrackID())) {
0190       // Create the propagation summary
0191       double absMomentum = track->GetMomentum().mag() * convertEnergy;
0192
0193       PropagationSummary iSummary(Acts::BoundTrackParameters::createCurvilinear(
0194           convertPosition(track->GetVertexPosition(), 0.),
0195           convertDirection(track->GetVertexMomentumDirection()),
0196           absCharge / absMomentum, std::nullopt,
0197           Acts::ParticleHypothesis::pion()));
0198
0199       eventStore().propagationRecords.insert({track->GetTrackID(), iSummary});
0200     }
0201     PropagationSummary& pSummary =
0202         eventStore().propagationRecords.at(track->GetTrackID());
0203
0204     // Increase the step counter
0205     pSummary.nSteps += 1;
0206
0207     double currentTrackLength = track->GetTrackLength() * convertLength;
0208     double currentStepLength = currentTrackLength - pSummary.pathLength;
0209     pSummary.pathLength = currentTrackLength;
0210
0211     // Create a new step for the step logging
0212     Acts::detail::Step pStep = stepFromG4Step(step);
0213     pStep.geoID = geoId;
0214     pStep.surface = surface != nullptr ? surface->getSharedPtr() : nullptr;
0215     // Check if last step was on same surface
0216     if (!pSummary.steps.empty() && pSummary.steps.back().geoID == geoId &&
0217         pSummary.steps.back().surface != nullptr) {
0218       auto& lastStep = pSummary.steps.back();
0219       lastStep.stepSize = Acts::ConstrainedStep(currentStepLength);
0220       lastStep.position = 0.5 * (pStep.position + lastStep.position);
0221       lastStep.momentum = 0.5 * (pStep.momentum + lastStep.momentum);
0222     } else {
0223       // Record the propagation state
0224       pStep.stepSize = Acts::ConstrainedStep(currentStepLength);
0225       pSummary.steps.emplace_back(std::move(pStep));
0226     }
0227     // You have nothing to do from here
0228     return;
0229   }
0230
0231   // Set particle hit count to zero, so we have this entry in the map later
0232   if (!eventStore().particleHitCount.contains(particleId)) {
0233     eventStore().particleHitCount[particleId] = 0;
0234   }
0235
0236   // Extract if we are at volume boundaries
0237   const bool preOnBoundary = preStepPoint->GetStepStatus() == fGeomBoundary;
0238   const bool postOnBoundary = postStepPoint->GetStepStatus() == fGeomBoundary ||
0239                               postStepPoint->GetStepStatus() == fWorldBoundary;
0240   const bool particleStopped = (postStepPoint->GetKineticEnergy() == 0.0);
0241   const bool particleDecayed =
0242       (postStepPoint->GetProcessDefinedStep()->GetProcessType() == fDecay);
0243
0244   auto print = [](auto s) {
0245 #if BOOST_VERSION >= 107800
0246     return boost::describe::enum_to_string(s, "unmatched");
0247 #else
0248     return s;
0249 #endif
0250   };
0251   ACTS_VERBOSE("status: pre="
0252                << print(preStepPoint->GetStepStatus())
0253                << ", post=" << print(postStepPoint->GetStepStatus())
0254                << ", post E_kin=" << std::boolalpha
0255                << postStepPoint->GetKineticEnergy() << ", process_type="
0256                << print(
0257                       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 auto 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