Geant4/src/ParticleTrackingAction.cpp

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