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