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/SimulationOutcome.hpp"
0017
0018 #include <cassert>
0019 #include <ostream>
0020 #include <unordered_map>
0021 #include <utility>
0022
0023 #include <G4ParticleDefinition.hh>
0024 #include <G4RunManager.hh>
0025 #include <G4Track.hh>
0026 #include <G4UnitsTable.hh>
0027
0028 namespace ActsExamples::Geant4 {
0029
0030 ParticleTrackingAction::ParticleTrackingAction(
0031     const Config& cfg, std::unique_ptr<const Acts::Logger> logger)
0032     : G4UserTrackingAction(), m_cfg(cfg), m_logger(std::move(logger)) {}
0033
0034 void ParticleTrackingAction::PreUserTrackingAction(const G4Track* trackPtr) {
0035   assert(trackPtr != nullptr);
0036   const G4Track& track = *trackPtr;
0037
0038   // If this is not the case, there are unhandled cases of particle stopping in
0039   // the SensitiveSteppingAction
0040   // TODO We could also merge the remaining hits to a hit here, but it would be
0041   // nicer to investigate, if we can handle all particle stop conditions in the
0042   // SensitiveSteppingAction... This seems to happen O(1) times in a ttbar
0043   // event, so seems not to be too problematic
0044   if (!eventStore().hitBuffer.empty()) {
0045     eventStore().hitBuffer.clear();
0046     ACTS_WARNING("Hit buffer not empty after track");
0047   }
0048
0049   const std::optional<SimBarcode> barcode =
0050       makeParticleId(track.GetTrackID(), track.GetParentID());
0051
0052   // There is already a warning printed in the makeParticleId function if this
0053   // indicates a failure
0054   if (!barcode.has_value()) {
0055     return;
0056   }
0057
0058   const SimParticleState fatrasParticle = convert(track, *barcode);
0059   const SimParticle particle(fatrasParticle, fatrasParticle);
0060   const auto [it, success] = eventStore().particlesInitial.insert(particle);
0061
0062   // Only register particle at the initial state AND if there is no particle ID
0063   // collision
0064   if (success) {
0065     eventStore().trackIdMapping[track.GetTrackID()] = particle.particleId();
0066   } else {
0067     eventStore().particleIdCollisionsInitial++;
0068     ACTS_WARNING("Particle ID collision with "
0069                  << particle.particleId()
0070                  << " detected for initial particles. Skip particle");
0071   }
0072 }
0073
0074 void ParticleTrackingAction::PostUserTrackingAction(const G4Track* trackPtr) {
0075   assert(trackPtr != nullptr);
0076   const G4Track& track = *trackPtr;
0077
0078   // The initial particle maybe was not registered because of a particle ID
0079   // collision
0080   if (!eventStore().trackIdMapping.contains(track.GetTrackID())) {
0081     ACTS_WARNING("Particle ID for track ID " << track.GetTrackID()
0082                                              << " not registered. Skip");
0083     return;
0084   }
0085
0086   const SimBarcode barcode = eventStore().trackIdMapping.at(track.GetTrackID());
0087
0088   const bool hasHits = eventStore().particleHitCount.contains(barcode) &&
0089                        eventStore().particleHitCount.at(barcode) > 0;
0090   if (!m_cfg.keepParticlesWithoutHits && !hasHits) {
0091     [[maybe_unused]] const std::size_t n =
0092         eventStore().particlesSimulated.erase(
0093             SimParticle(barcode, Acts::PdgParticle::eInvalid));
0094     assert(n == 1);
0095     return;
0096   }
0097
0098   const auto particleIt = eventStore().particlesInitial.find(barcode);
0099   if (particleIt == eventStore().particlesInitial.end()) {
0100     ACTS_WARNING("Particle ID " << barcode
0101                                 << " not found in initial particles");
0102     return;
0103   }
0104   SimParticle particle = *particleIt;
0105   particle.finalState() = convert(track, barcode);
0106
0107   const auto [it, success] = eventStore().particlesSimulated.insert(particle);
0108
0109   if (!success) {
0110     eventStore().particleIdCollisionsFinal++;
0111     ACTS_WARNING("Particle ID collision with "
0112                  << particle.particleId()
0113                  << " detected for final particles. Skip particle");
0114   }
0115 }
0116
0117 SimParticleState ParticleTrackingAction::convert(const G4Track& track,
0118                                                  SimBarcode particleId) const {
0119   // Get all the information from the Track
0120   const G4ParticleDefinition* particleDef = track.GetParticleDefinition();
0121   const G4int pdg = particleDef->GetPDGEncoding();
0122   const G4double charge = particleDef->GetPDGCharge();
0123   const G4double mass = convertEnergyToActs * particleDef->GetPDGMass();
0124   const G4ThreeVector pPosition = convertLengthToActs * track.GetPosition();
0125   const G4double pTime = convertTimeToActs * track.GetGlobalTime();
0126   const G4ThreeVector pDirection = track.GetMomentumDirection();
0127   const G4double p = convertEnergyToActs * track.GetKineticEnergy();
0128
0129   std::uint32_t numberOfHits = 0;
0130   if (const auto it = eventStore().particleHitCount.find(particleId);
0131       it != eventStore().particleHitCount.end()) {
0132     numberOfHits = it->second;
0133   }
0134
0135   ActsFatras::SimulationOutcome particleOutcome =
0136       ActsFatras::SimulationOutcome::Alive;
0137   if (const auto it = eventStore().particleOutcome.find(particleId);
0138       it != eventStore().particleOutcome.end()) {
0139     particleOutcome = it->second;
0140   }
0141
0142   // Now create the Particle
0143   SimParticleState aParticle(particleId, Acts::PdgParticle{pdg}, charge, mass);
0144   aParticle.setPosition4(pPosition[0], pPosition[1], pPosition[2], pTime);
0145   aParticle.setDirection(pDirection[0], pDirection[1], pDirection[2]);
0146   aParticle.setAbsoluteMomentum(p);
0147   aParticle.setNumberOfHits(numberOfHits);
0148   aParticle.setOutcome(particleOutcome);
0149   return aParticle;
0150 }
0151
0152 std::optional<SimBarcode> ParticleTrackingAction::makeParticleId(
0153     G4int trackId, G4int parentId) const {
0154   // We already have this particle registered (it is one of the input particles
0155   // or we are making a final particle state)
0156   if (eventStore().trackIdMapping.contains(trackId)) {
0157     return std::nullopt;
0158   }
0159
0160   if (!eventStore().trackIdMapping.contains(parentId)) {
0161     ACTS_DEBUG("Parent particle " << parentId
0162                                   << " not registered, cannot build barcode");
0163     eventStore().parentIdNotFound++;
0164     return std::nullopt;
0165   }
0166
0167   SimBarcode pid = eventStore().trackIdMapping.at(parentId).makeDescendant();
0168   const SimBarcode key = pid.withoutSubparticle();
0169   ++eventStore().subparticleMap[key];
0170   pid = pid.withSubParticle(eventStore().subparticleMap[key]);
0171
0172   return pid;
0173 }
0174
0175 }  // namespace ActsExamples::Geant4