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/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 {
0051
0052 std::array<Acts::Vector4, 4u> kinematicsOfStep(const G4Step* step) {
0053   const G4StepPoint* preStepPoint = step->GetPreStepPoint();
0054   const G4StepPoint* postStepPoint = step->GetPostStepPoint();
0055   using namespace ActsExamples::Geant4;
0056   Acts::Vector4 preStepPosition = convertPosition(
0057       preStepPoint->GetPosition(), preStepPoint->GetGlobalTime());
0058
0059   Acts::Vector4 preStepMomentum = convertMomentum(
0060       preStepPoint->GetMomentum(), preStepPoint->GetTotalEnergy());
0061   Acts::Vector4 postStepPosition = convertPosition(
0062       postStepPoint->GetPosition(), postStepPoint->GetGlobalTime());
0063
0064   Acts::Vector4 postStepMomentum = convertMomentum(
0065       postStepPoint->GetMomentum(), postStepPoint->GetTotalEnergy());
0066
0067   return {preStepPosition, preStepMomentum, postStepPosition, postStepMomentum};
0068 }
0069
0070 ActsFatras::Hit hitFromStep(const G4Step* step, ActsFatras::Barcode particleId,
0071                             Acts::GeometryIdentifier geoId,
0072                             std::int32_t index) {
0073   auto [preStepPosition, preStepMomentum, postStepPosition, postStepMomentum] =
0074       kinematicsOfStep(step);
0075
0076   return ActsFatras::Hit(geoId, particleId,
0077                          0.5 * (preStepPosition + postStepPosition),
0078                          preStepMomentum, postStepMomentum, index);
0079 }
0080
0081 Acts::detail::Step stepFromG4Step(const G4Step* step) {
0082   Acts::detail::Step pStep;
0083   auto [preStepPosition, preStepMomentum, postStepPosition, postStepMomentum] =
0084       kinematicsOfStep(step);
0085
0086   pStep.navDir = Acts::Direction::Forward();
0087   pStep.position = 0.5 * (preStepPosition + postStepPosition).block<3, 1>(0, 0);
0088   pStep.momentum = 0.5 * (preStepMomentum + postStepMomentum).block<3, 1>(0, 0);
0089   pStep.nTotalTrials = 1;
0090   return pStep;
0091 }
0092
0093 }  // namespace
0094
0095 namespace ActsExamples::Geant4 {
0096
0097 SensitiveSteppingAction::SensitiveSteppingAction(
0098     const Config& cfg, std::unique_ptr<const Acts::Logger> logger)
0099     : G4UserSteppingAction(), m_cfg(cfg), m_logger(std::move(logger)) {}
0100
0101 void SensitiveSteppingAction::UserSteppingAction(const G4Step* step) {
0102   // Unit conversions G4->::ACTS
0103   static constexpr double convertLength = Acts::UnitConstants::mm / CLHEP::mm;
0104   static constexpr double convertEnergy = Acts::UnitConstants::GeV / CLHEP::GeV;
0105   static constexpr auto mappingPrefix = SensitiveSurfaceMapper::mappingPrefix;
0106
0107   // The particle after the step
0108   G4Track* track = step->GetTrack();
0109   G4PrimaryParticle* primaryParticle =
0110       track->GetDynamicParticle()->GetPrimaryParticle();
0111
0112   // Get PreStepPoint and PostStepPoint
0113   const G4StepPoint* preStepPoint = step->GetPreStepPoint();
0114   const G4StepPoint* postStepPoint = step->GetPostStepPoint();
0115
0116   // Bail out if charged & configured to do so
0117   G4double absCharge = std::abs(track->GetParticleDefinition()->GetPDGCharge());
0118   if (!m_cfg.charged && absCharge > 0.) {
0119     return;
0120   }
0121
0122   // Bail out if neutral & configured to do so
0123   if (!m_cfg.neutral && absCharge == 0.) {
0124     return;
0125   }
0126
0127   // Bail out if it is a primary & configured to be ignored
0128   if (!m_cfg.primary && primaryParticle != nullptr) {
0129     return;
0130   }
0131
0132   // Bail out if it is a secondary & configured to be ignored
0133   if (!m_cfg.secondary && primaryParticle == nullptr) {
0134     return;
0135   }
0136
0137   // Get the physical volume & check if it has the sensitive string name
0138   const G4VPhysicalVolume* volume = track->GetVolume();
0139   if (volume == nullptr) {
0140     throw std::runtime_error("No volume found, terminate simulation");
0141   }
0142   std::string volumeName = volume->GetName();
0143   ACTS_VERBOSE("Check whether volume " << volumeName << " is sensitive");
0144   if (!m_cfg.stepLogging &&
0145       volumeName.find(mappingPrefix) == std::string::npos) {
0146     return;
0147   }
0148
0149   // The G4Touchable for the matching
0150   const G4VTouchable* touchable = track->GetTouchable();
0151
0152   Acts::GeometryIdentifier geoId{};
0153   const Acts::Surface* surface = nullptr;
0154
0155   // Find the range of candidate surfaces for the current position in the map
0156   const auto surfaceMap_itr = m_surfaceMapping.find(volume);
0157   if (surfaceMap_itr != m_surfaceMapping.end()) {
0158     const auto& surfacesToG4Vol = surfaceMap_itr->second;
0159     ACTS_VERBOSE("Found " << surfacesToG4Vol.size()
0160                           << " candidate surfaces for volume " << volumeName);
0161     const Acts::Vector3 volumePos =
0162         convertPosition(touchable->GetTranslation());
0163     const auto lookUp_itr = surfacesToG4Vol.find(volumePos);
0164     if (lookUp_itr == surfacesToG4Vol.end() && !m_cfg.stepLogging) {
0165       ACTS_ERROR("No candidate surfaces found for volume " << volumeName);
0166       return;
0167     }
0168     surface = lookUp_itr->second;
0169     geoId = surface->geometryId();
0170     ACTS_VERBOSE("Replica assignment successful -> to surface " << geoId);
0171   } else if (!m_cfg.stepLogging) {
0172     ACTS_ERROR("No candidate surfaces found for volume " << volumeName);
0173     return;
0174   }
0175   // This is not the case if we have a particle-ID collision
0176   if (!eventStore().trackIdMapping.contains(track->GetTrackID())) {
0177     return;
0178   }
0179
0180   // Output is only strictly valid if step logging is not enabled
0181   const auto particleId = eventStore().trackIdMapping.at(track->GetTrackID());
0182   if (!m_cfg.stepLogging && surface != nullptr) {
0183     ACTS_VERBOSE("Step of " << particleId << " in sensitive volume " << geoId);
0184   } else if (m_cfg.stepLogging) {
0185     if (!eventStore().propagationRecords.contains(track->GetTrackID())) {
0186       // Create the propagation summary
0187       double absMomentum = track->GetMomentum().mag() * convertEnergy;
0188
0189       PropagationSummary iSummary(Acts::BoundTrackParameters::createCurvilinear(
0190           convertPosition(track->GetVertexPosition(), 0.),
0191           convertDirection(track->GetVertexMomentumDirection()),
0192           absCharge / absMomentum, std::nullopt,
0193           Acts::ParticleHypothesis::pion()));
0194
0195       eventStore().propagationRecords.insert({track->GetTrackID(), iSummary});
0196     }
0197     PropagationSummary& pSummary =
0198         eventStore().propagationRecords.at(track->GetTrackID());
0199
0200     // Increase the step counter
0201     pSummary.nSteps += 1;
0202
0203     double currentTrackLength = track->GetTrackLength() * convertLength;
0204     double currentStepLength = currentTrackLength - pSummary.pathLength;
0205     pSummary.pathLength = currentTrackLength;
0206
0207     // Create a new step for the step logging
0208     Acts::detail::Step pStep = stepFromG4Step(step);
0209     pStep.geoID = geoId;
0210     pStep.surface = surface != nullptr ? surface->getSharedPtr() : nullptr;
0211     // Check if last step was on same surface
0212     if (!pSummary.steps.empty() && pSummary.steps.back().geoID == geoId &&
0213         pSummary.steps.back().surface != nullptr) {
0214       auto& lastStep = pSummary.steps.back();
0215       lastStep.stepSize = Acts::ConstrainedStep(currentStepLength);
0216       lastStep.position = 0.5 * (pStep.position + lastStep.position);
0217       lastStep.momentum = 0.5 * (pStep.momentum + lastStep.momentum);
0218     } else {
0219       // Record the propagation state
0220       pStep.stepSize = Acts::ConstrainedStep(currentStepLength);
0221       pSummary.steps.emplace_back(std::move(pStep));
0222     }
0223     // You have nothing to do from here
0224     return;
0225   }
0226
0227   // Set particle hit count to zero, so we have this entry in the map later
0228   if (!eventStore().particleHitCount.contains(particleId)) {
0229     eventStore().particleHitCount[particleId] = 0;
0230   }
0231
0232   // Extract if we are at volume boundaries
0233   const bool preOnBoundary = preStepPoint->GetStepStatus() == fGeomBoundary;
0234   const bool postOnBoundary = postStepPoint->GetStepStatus() == fGeomBoundary ||
0235                               postStepPoint->GetStepStatus() == fWorldBoundary;
0236   const bool particleStopped = (postStepPoint->GetKineticEnergy() == 0.0);
0237   const bool particleDecayed =
0238       (postStepPoint->GetProcessDefinedStep()->GetProcessType() == fDecay);
0239
0240   auto print = [](auto s) {
0241     return boost::describe::enum_to_string(s, "unmatched");
0242   };
0243
0244   ACTS_VERBOSE("status: pre="
0245                << print(preStepPoint->GetStepStatus())
0246                << ", post=" << print(postStepPoint->GetStepStatus())
0247                << ", post E_kin=" << std::boolalpha
0248                << postStepPoint->GetKineticEnergy() << ", process_type="
0249                << print(
0250                       postStepPoint->GetProcessDefinedStep()->GetProcessType())
0251                << ", particle="
0252                << track->GetParticleDefinition()->GetParticleName()
0253                << ", process_name="
0254                << postStepPoint->GetProcessDefinedStep()->GetProcessName()
0255                << ", track status=" << print(track->GetTrackStatus()));
0256
0257   // Case A: The step starts at the entry of the volume and ends at the exit.
0258   // Add hit to collection.
0259   if (preOnBoundary && postOnBoundary) {
0260     ACTS_VERBOSE("-> merge single step to hit");
0261     ++eventStore().particleHitCount[particleId];
0262     eventStore().hits.push_back(
0263         hitFromStep(step, particleId, geoId,
0264                     eventStore().particleHitCount.at(particleId) - 1));
0265
0266     eventStore().numberGeantSteps += 1ul;
0267     eventStore().maxStepsForHit = std::max(eventStore().maxStepsForHit, 1ul);
0268     return;
0269   }
0270
0271   // Case B: The step ends at the exit of the volume. Add hit to hit buffer, and
0272   // then combine hit buffer
0273   if (postOnBoundary || particleStopped || particleDecayed) {
0274     ACTS_VERBOSE("-> merge buffer to hit");
0275     auto& buffer = eventStore().hitBuffer;
0276     buffer.push_back(hitFromStep(step, particleId, geoId, -1));
0277
0278     const auto pos4 =
0279         0.5 * (buffer.front().fourPosition() + buffer.back().fourPosition());
0280
0281     ++eventStore().particleHitCount[particleId];
0282     eventStore().hits.emplace_back(
0283         geoId, particleId, pos4, buffer.front().momentum4Before(),
0284         buffer.back().momentum4After(),
0285         eventStore().particleHitCount.at(particleId) - 1);
0286
0287     assert(std::ranges::all_of(
0288         buffer, [&](const auto& h) { return h.geometryId() == geoId; }));
0289     assert(std::ranges::all_of(
0290         buffer, [&](const auto& h) { return h.particleId() == particleId; }));
0291
0292     eventStore().numberGeantSteps += buffer.size();
0293     eventStore().maxStepsForHit =
0294         std::max(eventStore().maxStepsForHit, buffer.size());
0295
0296     buffer.clear();
0297     return;
0298   }
0299
0300   // Case C: The step doesn't end at the exit of the volume. Add the hit to the
0301   // hit buffer.
0302   if (!postOnBoundary) {
0303     // ACTS_VERBOSE("-> add hit to buffer");
0304     eventStore().hitBuffer.push_back(hitFromStep(step, particleId, geoId, -1));
0305     return;
0306   }
0307
0308   assert(false && "should never reach this");
0309 }
0310
0311 }  // namespace ActsExamples::Geant4