Geant4/src/Geant4Decay.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 "ActsFatras/Geant4/Geant4Decay.hpp"
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Definitions/Common.hpp"
0013 #include "ActsFatras/EventData/Barcode.hpp"
0014 #include "ActsFatras/EventData/ProcessType.hpp"
0015 #include "ActsFatras/Geant4/DummyDetectorConstruction.hpp"
0016
0017 #include <cstdint>
0018 #include <utility>
0019
0020 #include "G4DecayProducts.hh"
0021 #include "G4DecayTable.hh"
0022
0023 ActsFatras::Geant4Decay::Geant4Decay()
0024     : m_g4RunManager(ensureGeant4RunManager()) {}
0025
0026 std::vector<ActsFatras::Particle> ActsFatras::Geant4Decay::decayParticle(
0027     const ActsFatras::Particle& parent) const {
0028   std::vector<Particle> children;
0029
0030   // Find the particle type that will decay
0031   G4ParticleDefinition* pDef =
0032       m_pdgToG4Conv.getParticleDefinition(parent.pdg());
0033   if (pDef == nullptr) {
0034     return children;
0035   }
0036   // Get the particles decay table
0037   G4DecayTable* dt = pDef->GetDecayTable();
0038   if (dt == nullptr) {
0039     return children;
0040   }
0041   // Select a decay channel
0042   G4VDecayChannel* channel = dt->SelectADecayChannel();
0043   if (channel == nullptr) {
0044     return children;
0045   }
0046   // Get the decay products from the selected channel
0047   G4DecayProducts* products = channel->DecayIt();
0048   if (products == nullptr) {
0049     return children;
0050   }
0051
0052   // Boost the decay products using the parents four-momentum
0053   const Acts::Vector4 mom4 = parent.fourMomentum();
0054   products->Boost(mom4[Acts::eMom0] / mom4[Acts::eEnergy],
0055                   mom4[Acts::eMom1] / mom4[Acts::eEnergy],
0056                   mom4[Acts::eMom2] / mom4[Acts::eEnergy]);
0057
0058   G4int nProducts = products->entries();
0059   for (G4int i = 0; i < nProducts; i++) {
0060     G4DynamicParticle* prod = products->PopProducts();
0061     if (prod == nullptr) {
0062       continue;
0063     }
0064
0065     // Convert the decay product from Geant4 to Acts
0066     const G4ThreeVector& mom = prod->GetMomentum();
0067     constexpr double convertEnergy = Acts::UnitConstants::GeV / CLHEP::GeV;
0068     Acts::Vector3 amgMom(mom.x(), mom.y(), mom.z());
0069     amgMom *= convertEnergy;
0070     const std::int32_t pdg = prod->GetPDGcode();
0071
0072     Particle childParticle(parent.particleId().makeDescendant(i),
0073                            static_cast<Acts::PdgParticle>(pdg));
0074     childParticle.setPosition4(parent.fourPosition())
0075         .setAbsoluteMomentum(amgMom.norm())
0076         .setDirection(amgMom)
0077         .setProcess(ProcessType::eDecay);
0078
0079     // Store the particle
0080     children.push_back(std::move(childParticle));
0081   }
0082   return children;
0083 }