Geant4/src/SimParticleTranslation.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/SimParticleTranslation.hpp"
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "ActsExamples/EventData/SimParticle.hpp"
0013 #include "ActsExamples/Framework/DataHandle.hpp"
0014 #include "ActsExamples/Geant4/EventStore.hpp"
0015 #include "ActsExamples/Geant4/UnitConversion.hpp"
0016
0017 #include <ostream>
0018 #include <unordered_map>
0019 #include <utility>
0020
0021 #include <G4ChargedGeantino.hh>
0022 #include <G4Event.hh>
0023 #include <G4Geantino.hh>
0024 #include <G4ParticleDefinition.hh>
0025 #include <G4ParticleTable.hh>
0026 #include <G4PrimaryParticle.hh>
0027 #include <G4PrimaryVertex.hh>
0028 #include <G4UnitsTable.hh>
0029
0030 namespace ActsExamples::Geant4 {
0031
0032 SimParticleTranslation::SimParticleTranslation(
0033     const Config& cfg, std::unique_ptr<const Acts::Logger> logger)
0034     : G4VUserPrimaryGeneratorAction(),
0035       m_cfg(cfg),
0036       m_logger(std::move(logger)) {}
0037
0038 void SimParticleTranslation::GeneratePrimaries(G4Event* eventPtr) {
0039   assert(eventPtr != nullptr);
0040   G4Event& event = *eventPtr;
0041
0042   event.SetEventID(m_eventNr++);
0043   std::uint32_t eventID = event.GetEventID();
0044
0045   ACTS_DEBUG("Primary Generator Action for Event: " << eventID);
0046
0047   if (eventStore().store == nullptr) {
0048     ACTS_WARNING("No WhiteBoard instance could be found for this event!");
0049     return;
0050   }
0051
0052   if (eventStore().inputParticles == nullptr) {
0053     ACTS_WARNING("No input particle handle found");
0054     return;
0055   }
0056
0057   // Get the number of input particles
0058   const auto inputParticles =
0059       (*eventStore().inputParticles)(*eventStore().store);
0060
0061   // Reserve hopefully enough hit space
0062   eventStore().hits.reserve(inputParticles.size() *
0063                             m_cfg.reserveHitsPerParticle);
0064
0065   // Default particle kinematic
0066   G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
0067   G4PrimaryVertex* pVertex = nullptr;
0068
0069   // We are looping through the particles and flush per vertex
0070   std::optional<Acts::Vector4> lastVertex;
0071
0072   std::uint32_t pCounter = 0;
0073   std::uint32_t trackId = 1;
0074   // Loop over the input partilces and run
0075   for (const auto& part : inputParticles) {
0076     const Acts::Vector4 currentVertex = part.fourPosition();
0077     if (!lastVertex || !currentVertex.isApprox(*lastVertex)) {
0078       // Add the vertex to the event
0079       if (pVertex != nullptr) {
0080         event.AddPrimaryVertex(pVertex);
0081         ACTS_DEBUG("Flushing " << pCounter
0082                                << " particles associated with vertex "
0083                                << lastVertex->transpose());
0084         pCounter = 0;
0085       }
0086       lastVertex = currentVertex;
0087       pVertex = new G4PrimaryVertex(currentVertex[0] * convertLengthToGeant4,
0088                                     currentVertex[1] * convertLengthToGeant4,
0089                                     currentVertex[2] * convertLengthToGeant4,
0090                                     currentVertex[3] * convertTimeToGeant4);
0091     }
0092
0093     // Add a new primary to the vertex
0094
0095     const Acts::Vector4 mom4 = part.fourMomentum() * convertEnergyToGeant4;
0096
0097     // Particle properties, may be forced to specific value
0098     const G4int particlePdgCode = m_cfg.forcedPdgCode.value_or(part.pdg());
0099     const G4double particleCharge = m_cfg.forcedCharge.value_or(part.charge());
0100     const G4double particleMass =
0101         m_cfg.forcedMass.value_or(part.mass() * convertEnergyToGeant4);
0102
0103     // Check if it is a Geantino / ChargedGeantino
0104     const G4ParticleDefinition* particleDefinition =
0105         particleTable->FindParticle(particlePdgCode);
0106     if (particleDefinition == nullptr) {
0107       if (particlePdgCode == 0 && particleMass == 0 && particleCharge == 0) {
0108         particleDefinition = G4Geantino::Definition();
0109       }
0110       if (particlePdgCode == 0 && particleMass == 0 && particleCharge != 0) {
0111         if (particleCharge != 1) {
0112           ACTS_ERROR("invalid charged geantino charge " << particleCharge
0113                                                         << ". should be 1");
0114         }
0115         particleDefinition = G4ChargedGeantino::Definition();
0116       }
0117     }
0118
0119     // Skip if translation failed
0120     if (particleDefinition == nullptr) {
0121       ACTS_DEBUG(
0122           "Could not translate particle with PDG code : " << particlePdgCode);
0123       continue;
0124     }
0125
0126     ACTS_VERBOSE("Adding particle with name '"
0127                  << particleDefinition->GetParticleName()
0128                  << "' and properties:");
0129     ACTS_VERBOSE(" -> mass: " << particleMass);
0130     ACTS_VERBOSE(" -> charge: " << particleCharge);
0131     ACTS_VERBOSE(" -> momentum: " << mom4.transpose());
0132
0133     // G4 will delete this
0134     G4PrimaryParticle* particle = new G4PrimaryParticle(particleDefinition);
0135
0136     particle->SetMass(particleMass);
0137     particle->SetCharge(particleCharge);
0138     particle->Set4Momentum(mom4[0], mom4[1], mom4[2], mom4[3]);
0139     particle->SetTrackID(trackId++);
0140
0141     // Add the primary to the vertex
0142     pVertex->SetPrimary(particle);
0143
0144     eventStore().particlesInitial.insert(part);
0145     eventStore().trackIdMapping[particle->GetTrackID()] = part.particleId();
0146
0147     ++pCounter;
0148   }
0149   // Final vertex to be added
0150   if (pVertex != nullptr) {
0151     event.AddPrimaryVertex(pVertex);
0152     ACTS_DEBUG("Flushing " << pCounter << " particles associated with vertex "
0153                            << lastVertex->transpose());
0154   }
0155 }
0156
0157 }  // namespace ActsExamples::Geant4