Geant4/src/SimParticleTranslation.cpp

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