Geant4/src/ParticleTrackingAction.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/ParticleTrackingAction.hpp"
0010
0011 #include "Acts/Definitions/PdgParticle.hpp"
0012 #include "ActsExamples/EventData/SimParticle.hpp"
0013 #include "ActsExamples/Geant4/EventStore.hpp"
0014 #include "ActsExamples/Geant4/UnitConversion.hpp"
0015 #include "ActsFatras/EventData/Barcode.hpp"
0016 #include "ActsFatras/EventData/GenerationProcess.hpp"
0017 #include "ActsFatras/EventData/SimulationOutcome.hpp"
0018
0019 #include <cassert>
0020 #include <ostream>
0021 #include <unordered_map>
0022 #include <utility>
0023
0024 #include <G4EmProcessSubType.hh>
0025 #include <G4ParticleDefinition.hh>
0026 #include <G4ProcessType.hh>
0027 #include <G4RunManager.hh>
0028 #include <G4Track.hh>
0029 #include <G4UnitsTable.hh>
0030 #include <G4VProcess.hh>
0031
0032 namespace ActsExamples::Geant4 {
0033
0034 namespace {
0035
0036 /// Map the Geant4 process that created a particle to an ActsFatras
0037 /// GenerationProcess enum value, so it can be stored on the SimParticle.
0038 /// Primaries have no creator process.
0039 ///
0040 /// We translate the G4 creator process into the finer GenerationProcess codes:
0041 /// - decay                          -> eDecay        (G4ProcessType fDecay)
0042 /// - true nuclear / photo-nuclear   -> eNuclearInteraction
0043 ///                                         (fHadronic / fPhotolepton_hadron)
0044 /// EM processes defined in G4EmProcessSubType.hh:
0045 /// - bremsstrahlung                -> eBremsstrahlung  (fBremsstrahlung in G4)
0046 /// - photon conversion (pair prod.)-> ePhotonConversion(fGammaConversion in G4)
0047 /// - ionisation (delta-ray)        -> eIonisation      (fIonisation in G4)
0048 /// - anything else (Compton, photo-
0049 ///   electric, annihilation, ...)  -> eOther (usually rare)
0050 /// Primaries have no creator process and stay eUndefined.
0051 ActsFatras::GenerationProcess g4CreatorToGenerationProcess(
0052     const G4VProcess* proc) {
0053   using ActsFatras::GenerationProcess;
0054   if (proc == nullptr) {
0055     return GenerationProcess::eUndefined;  // primary / no creator process
0056   }
0057   const G4ProcessType type = proc->GetProcessType();
0058   if (type == fDecay) {
0059     return GenerationProcess::eDecay;
0060   }
0061   if (type == fHadronic || type == fPhotolepton_hadron) {
0062     return GenerationProcess::eNuclearInteraction;  // genuinely nuclear
0063   }
0064   switch (proc->GetProcessSubType()) {
0065     case fBremsstrahlung:
0066       return GenerationProcess::eBremsstrahlung;
0067     case fGammaConversion:
0068       return GenerationProcess::ePhotonConversion;
0069     case fIonisation:
0070       return GenerationProcess::eIonisation;  // delta-ray
0071     default:
0072       return GenerationProcess::eOther;  // residual EM (Compton, photoelectric,
0073                                          // ...)
0074   }
0075 }
0076
0077 }  // namespace
0078
0079 ParticleTrackingAction::ParticleTrackingAction(
0080     const Config& cfg, std::unique_ptr<const Acts::Logger> logger)
0081     : G4UserTrackingAction(), m_cfg(cfg), m_logger(std::move(logger)) {}
0082
0083 void ParticleTrackingAction::PreUserTrackingAction(const G4Track* trackPtr) {
0084   assert(trackPtr != nullptr);
0085   const G4Track& track = *trackPtr;
0086
0087   // If this is not the case, there are unhandled cases of particle stopping in
0088   // the SensitiveSteppingAction
0089   // TODO We could also merge the remaining hits to a hit here, but it would be
0090   // nicer to investigate, if we can handle all particle stop conditions in the
0091   // SensitiveSteppingAction... This seems to happen O(1) times in a ttbar
0092   // event, so seems not to be too problematic
0093   if (!eventStore().hitBuffer.empty()) {
0094     eventStore().hitBuffer.clear();
0095     ACTS_WARNING("Hit buffer not empty after track");
0096   }
0097
0098   const std::optional<SimBarcode> barcode =
0099       makeParticleId(track.GetTrackID(), track.GetParentID());
0100
0101   // There is already a warning printed in the makeParticleId function if this
0102   // indicates a failure
0103   if (!barcode.has_value()) {
0104     return;
0105   }
0106
0107   const SimParticleState fatrasParticle = convert(track, *barcode);
0108   SimParticle particle(fatrasParticle, fatrasParticle);
0109
0110   // Record the parent particle so parentParticleId() is available downstream.
0111   // Secondaries have their G4 parent registered in trackIdMapping; input
0112   // primaries have parentId == 0 (not registered) and keep the default parent.
0113   if (const auto pit = eventStore().trackIdMapping.find(track.GetParentID());
0114       pit != eventStore().trackIdMapping.end()) {
0115     particle.setParentParticleId(pit->second);
0116   }
0117   const auto [it, success] = eventStore().particlesInitial.insert(particle);
0118
0119   // Only register particle at the initial state AND if there is no particle ID
0120   // collision
0121   if (success) {
0122     eventStore().trackIdMapping[track.GetTrackID()] = particle.particleId();
0123   } else {
0124     eventStore().particleIdCollisionsInitial++;
0125     ACTS_WARNING("Particle ID collision with "
0126                  << particle.particleId()
0127                  << " detected for initial particles. Skip particle");
0128   }
0129 }
0130
0131 void ParticleTrackingAction::PostUserTrackingAction(const G4Track* trackPtr) {
0132   assert(trackPtr != nullptr);
0133   const G4Track& track = *trackPtr;
0134
0135   // The initial particle maybe was not registered because of a particle ID
0136   // collision
0137   if (!eventStore().trackIdMapping.contains(track.GetTrackID())) {
0138     ACTS_WARNING("Particle ID for track ID " << track.GetTrackID()
0139                                              << " not registered. Skip");
0140     return;
0141   }
0142
0143   const SimBarcode barcode = eventStore().trackIdMapping.at(track.GetTrackID());
0144
0145   const bool hasHits = eventStore().particleHitCount.contains(barcode) &&
0146                        eventStore().particleHitCount.at(barcode) > 0;
0147   if (!m_cfg.keepParticlesWithoutHits && !hasHits) {
0148     [[maybe_unused]] const std::size_t n =
0149         eventStore().particlesSimulated.erase(
0150             SimParticle(barcode, Acts::PdgParticle::eInvalid));
0151     assert(n == 1);
0152     return;
0153   }
0154
0155   const auto particleIt = eventStore().particlesInitial.find(barcode);
0156   if (particleIt == eventStore().particlesInitial.end()) {
0157     ACTS_WARNING("Particle ID " << barcode
0158                                 << " not found in initial particles");
0159     return;
0160   }
0161   SimParticle particle = *particleIt;
0162   particle.finalState() = convert(track, barcode);
0163   // set parent id from the initial state so both states carry it consistently
0164   particle.setParentParticleId(particle.parentParticleId());
0165
0166   const auto [it, success] = eventStore().particlesSimulated.insert(particle);
0167
0168   if (!success) {
0169     eventStore().particleIdCollisionsFinal++;
0170     ACTS_WARNING("Particle ID collision with "
0171                  << particle.particleId()
0172                  << " detected for final particles. Skip particle");
0173   }
0174 }
0175
0176 SimParticleState ParticleTrackingAction::convert(const G4Track& track,
0177                                                  SimBarcode particleId) const {
0178   // Get all the information from the Track
0179   const G4ParticleDefinition* particleDef = track.GetParticleDefinition();
0180   const G4int pdg = particleDef->GetPDGEncoding();
0181   const G4double charge = particleDef->GetPDGCharge();
0182   const G4double mass = convertEnergyToActs * particleDef->GetPDGMass();
0183   const G4ThreeVector pPosition = convertLengthToActs * track.GetPosition();
0184   const G4double pTime = convertTimeToActs * track.GetGlobalTime();
0185   const G4ThreeVector pDirection = track.GetMomentumDirection();
0186   const G4double p = convertEnergyToActs * track.GetMomentum().mag();
0187
0188   std::uint32_t numberOfHits = 0;
0189   if (const auto it = eventStore().particleHitCount.find(particleId);
0190       it != eventStore().particleHitCount.end()) {
0191     numberOfHits = it->second;
0192   }
0193
0194   ActsFatras::SimulationOutcome particleOutcome =
0195       ActsFatras::SimulationOutcome::Alive;
0196   if (const auto it = eventStore().particleOutcome.find(particleId);
0197       it != eventStore().particleOutcome.end()) {
0198     particleOutcome = it->second;
0199   }
0200
0201   // Now create the Particle
0202   SimParticleState aParticle(particleId, Acts::PdgParticle{pdg}, charge, mass);
0203   aParticle.setPosition4(pPosition[0], pPosition[1], pPosition[2], pTime);
0204   aParticle.setDirection(pDirection[0], pDirection[1], pDirection[2]);
0205   aParticle.setAbsoluteMomentum(p);
0206   aParticle.setNumberOfHits(numberOfHits);
0207   aParticle.setOutcome(particleOutcome);
0208   // Record which Geant4 process created this particle (decay / material / ...).
0209   // Primaries have no creator process and stay eUndefined.
0210   aParticle.setProcess(g4CreatorToGenerationProcess(track.GetCreatorProcess()));
0211   return aParticle;
0212 }
0213
0214 std::optional<SimBarcode> ParticleTrackingAction::makeParticleId(
0215     G4int trackId, G4int parentId) const {
0216   // We already have this particle registered (it is one of the input particles
0217   // or we are making a final particle state)
0218   if (eventStore().trackIdMapping.contains(trackId)) {
0219     return std::nullopt;
0220   }
0221
0222   if (!eventStore().trackIdMapping.contains(parentId)) {
0223     ACTS_DEBUG("Parent particle " << parentId
0224                                   << " not registered, cannot build barcode");
0225     eventStore().parentIdNotFound++;
0226     return std::nullopt;
0227   }
0228
0229   SimBarcode pid = eventStore().trackIdMapping.at(parentId).makeDescendant();
0230   const SimBarcode key = pid.withoutSubparticle();
0231   ++eventStore().subparticleMap[key];
0232   pid = pid.withSubParticle(eventStore().subparticleMap[key]);
0233
0234   return pid;
0235 }
0236
0237 }  // namespace ActsExamples::Geant4