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 "Acts/Definitions/Units.hpp"
0013 #include "ActsExamples/EventData/SimParticle.hpp"
0014 #include "ActsExamples/Framework/DataHandle.hpp"
0015 #include "ActsExamples/Framework/WhiteBoard.hpp"
0016 #include "ActsExamples/Geant4/EventStore.hpp"
0017
0018 #include <ostream>
0019 #include <unordered_map>
0020 #include <utility>
0021
0022 #include <G4ChargedGeantino.hh>
0023 #include <G4Event.hh>
0024 #include <G4Geantino.hh>
0025 #include <G4ParticleDefinition.hh>
0026 #include <G4ParticleTable.hh>
0027 #include <G4PrimaryParticle.hh>
0028 #include <G4PrimaryVertex.hh>
0029 #include <G4UnitsTable.hh>
0030
0031 namespace ActsExamples::Geant4 {
0032
0033 SimParticleTranslation::SimParticleTranslation(
0034     const Config& cfg, std::unique_ptr<const Acts::Logger> logger)
0035     : G4VUserPrimaryGeneratorAction(),
0036       m_cfg(cfg),
0037       m_logger(std::move(logger)) {}
0038
0039 SimParticleTranslation::~SimParticleTranslation() = default;
0040
0041 void SimParticleTranslation::GeneratePrimaries(G4Event* anEvent) {
0042   anEvent->SetEventID(m_eventNr++);
0043   unsigned int eventID = anEvent->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   constexpr double convertLength = CLHEP::mm / Acts::UnitConstants::mm;
0073   constexpr double convertTime = CLHEP::ns / Acts::UnitConstants::ns;
0074   constexpr double convertEnergy = CLHEP::GeV / Acts::UnitConstants::GeV;
0075
0076   unsigned int pCounter = 0;
0077   unsigned int trackId = 1;
0078   // Loop over the input partilces and run
0079   for (const auto& part : inputParticles) {
0080     auto currentVertex = part.fourPosition();
0081     if (!lastVertex || !currentVertex.isApprox(*lastVertex)) {
0082       // Add the vertex to the event
0083       if (pVertex != nullptr) {
0084         anEvent->AddPrimaryVertex(pVertex);
0085         ACTS_DEBUG("Flushing " << pCounter
0086                                << " particles associated with vertex "
0087                                << lastVertex->transpose());
0088         pCounter = 0;
0089       }
0090       lastVertex = currentVertex;
0091       pVertex = new G4PrimaryVertex(
0092           currentVertex[0] * convertLength, currentVertex[1] * convertLength,
0093           currentVertex[2] * convertLength, currentVertex[3] * convertTime);
0094     }
0095
0096     // Add a new primary to the vertex
0097
0098     Acts::Vector4 mom4 = part.fourMomentum() * convertEnergy;
0099
0100     // Particle properties, may be forced to specific value
0101     G4int particlePdgCode = m_cfg.forcedPdgCode.value_or(part.pdg());
0102     G4double particleCharge = m_cfg.forcedCharge.value_or(part.charge());
0103     G4double particleMass =
0104         m_cfg.forcedMass.value_or(part.mass() * convertEnergy);
0105
0106     // Check if it is a Geantino / ChargedGeantino
0107     G4ParticleDefinition* particleDefinition =
0108         particleTable->FindParticle(particlePdgCode);
0109     if (particleDefinition == nullptr) {
0110       if (particlePdgCode == 0 && particleMass == 0 && particleCharge == 0) {
0111         particleDefinition = G4Geantino::Definition();
0112       }
0113       if (particlePdgCode == 0 && particleMass == 0 && particleCharge != 0) {
0114         if (particleCharge != 1) {
0115           ACTS_ERROR("invalid charged geantino charge " << particleCharge
0116                                                         << ". should be 1");
0117         }
0118         particleDefinition = G4ChargedGeantino::Definition();
0119       }
0120     }
0121
0122     // Skip if translation failed
0123     if (particleDefinition == nullptr) {
0124       ACTS_DEBUG(
0125           "Could not translate particle with PDG code : " << particlePdgCode);
0126       continue;
0127     }
0128
0129     ACTS_VERBOSE("Adding particle with name '"
0130                  << particleDefinition->GetParticleName()
0131                  << "' and properties:");
0132     ACTS_VERBOSE(" -> mass: " << particleMass);
0133     ACTS_VERBOSE(" -> charge: " << particleCharge);
0134     ACTS_VERBOSE(" -> momentum: " << mom4.transpose());
0135
0136     // G4 will delete this
0137     G4PrimaryParticle* particle = new G4PrimaryParticle(particleDefinition);
0138
0139     particle->SetMass(particleMass);
0140     particle->SetCharge(particleCharge);
0141     particle->Set4Momentum(mom4[0], mom4[1], mom4[2], mom4[3]);
0142     particle->SetTrackID(trackId++);
0143
0144     // Add the primary to the vertex
0145     pVertex->SetPrimary(particle);
0146
0147     eventStore().particlesInitial.insert(part);
0148     eventStore().trackIdMapping[particle->GetTrackID()] = part.particleId();
0149
0150     ++pCounter;
0151   }
0152   // Final vertex to be added
0153   if (pVertex != nullptr) {
0154     anEvent->AddPrimaryVertex(pVertex);
0155     ACTS_DEBUG("Flushing " << pCounter << " particles associated with vertex "
0156                            << lastVertex->transpose());
0157   }
0158 }
0159
0160 }  // namespace ActsExamples::Geant4