Geant4/src/Geant4MaterialStepConverter.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 "ActsPlugins/Geant4/Geant4MaterialStepConverter.hpp"
0010
0011 #include "Acts/Definitions/Units.hpp"
0012 #include "Acts/Material/Material.hpp"
0013 #include "Acts/Material/MaterialSlab.hpp"
0014
0015 #include <algorithm>
0016 #include <cstddef>
0017 #include <optional>
0018 #include <utility>
0019
0020 #include <CLHEP/Units/SystemOfUnits.h>
0021 #include <G4Material.hh>
0022 #include <G4Step.hh>
0023 #include <G4Track.hh>
0024
0025 namespace {
0026
0027 constexpr double convertLengthToActs = Acts::UnitConstants::mm / CLHEP::mm;
0028 constexpr double convertDensityToActs =
0029     (Acts::UnitConstants::g / Acts::UnitConstants::mm3) /
0030     (CLHEP::gram / CLHEP::mm3);
0031
0032 Acts::Vector3 convertPosition(const G4ThreeVector& g4vec) {
0033   return Acts::Vector3(g4vec[0], g4vec[1], g4vec[2]) * convertLengthToActs;
0034 }
0035
0036 Acts::Vector3 convertDirection(const G4ThreeVector& g4vec) {
0037   return Acts::Vector3(g4vec[0], g4vec[1], g4vec[2]);
0038 }
0039
0040 }  // namespace
0041
0042 namespace ActsPlugins {
0043
0044 Acts::RecordedMaterialTrack Geant4MaterialStepConverter::startTrack(
0045     const G4Track& track) const {
0046   Acts::RecordedMaterialTrack rmTrack;
0047   rmTrack.first = {convertPosition(track.GetVertexPosition()),
0048                    convertDirection(track.GetVertexMomentumDirection())};
0049   return rmTrack;
0050 }
0051
0052 bool Geant4MaterialStepConverter::addStep(
0053     const G4Step& step, Acts::RecordedMaterialTrack& track) const {
0054   std::optional<Acts::MaterialInteraction> mInteraction =
0055       materialInteraction(step);
0056   if (!mInteraction.has_value()) {
0057     return false;
0058   }
0059   Acts::RecordedMaterial& recorded = track.second;
0060   recorded.materialInX0 += mInteraction->materialSlab.thicknessInX0();
0061   recorded.materialInL0 += mInteraction->materialSlab.thicknessInL0();
0062   recorded.materialInteractions.push_back(std::move(*mInteraction));
0063   return true;
0064 }
0065
0066 std::optional<Acts::MaterialInteraction>
0067 Geant4MaterialStepConverter::materialInteraction(const G4Step& step) const {
0068   const G4Material* materialPtr = step.GetPreStepPoint()->GetMaterial();
0069   if (materialPtr == nullptr) {
0070     return std::nullopt;
0071   }
0072   const G4Material& material = *materialPtr;
0073
0074   if (std::ranges::find(m_cfg.excludeMaterials, material.GetName()) !=
0075       m_cfg.excludeMaterials.end()) {
0076     return std::nullopt;
0077   }
0078
0079   // Quantities valid for elemental materials and mixtures
0080   const double X0 = convertLengthToActs * material.GetRadlen();
0081   const double L0 = convertLengthToActs * material.GetNuclearInterLength();
0082   const double rho = convertDensityToActs * material.GetDensity();
0083
0084   // Get{A,Z} is only meaningful for single-element materials (according to
0085   // the Geant4 docs). Need to compute average manually.
0086   const G4ElementVector* elements = material.GetElementVector();
0087   const G4double* fraction = material.GetFractionVector();
0088   const std::size_t nElements = material.GetNumberOfElements();
0089   double Ar = 0.;
0090   double Z = 0.;
0091   if (nElements == 1) {
0092     Ar = material.GetA() * (CLHEP::mole / CLHEP::gram);
0093     Z = material.GetZ();
0094   } else {
0095     for (std::size_t i = 0; i < nElements; i++) {
0096       Ar += elements->at(i)->GetA() * fraction[i] * (CLHEP::mole / CLHEP::gram);
0097       Z += elements->at(i)->GetZ() * fraction[i];
0098     }
0099   }
0100
0101   // Construct passed material slab for the step, ACTS stores the material in
0102   // single precision
0103   const Acts::MaterialSlab slab(
0104       Acts::Material::fromMassDensity(
0105           static_cast<float>(X0), static_cast<float>(L0),
0106           static_cast<float>(Ar), static_cast<float>(Z),
0107           static_cast<float>(rho)),
0108       static_cast<float>(convertLengthToActs * step.GetStepLength()));
0109
0110   Acts::MaterialInteraction mInteraction;
0111   mInteraction.position =
0112       convertPosition(step.GetPreStepPoint()->GetPosition());
0113   mInteraction.direction =
0114       convertDirection(step.GetPreStepPoint()->GetMomentum()).normalized();
0115   mInteraction.materialSlab = slab;
0116   mInteraction.pathCorrection = step.GetStepLength() * convertLengthToActs;
0117
0118   if (m_cfg.recordElementFractions) {
0119     if (nElements == 1) {
0120       mInteraction.elementZ.push_back(
0121           static_cast<unsigned int>(material.GetZ()));
0122       mInteraction.elementFrac.push_back(1.0f);
0123     } else {
0124       for (std::size_t i = 0; i < nElements; i++) {
0125         mInteraction.elementZ.push_back(
0126             static_cast<unsigned int>(elements->at(i)->GetZ()));
0127         mInteraction.elementFrac.push_back(static_cast<float>(fraction[i]));
0128       }
0129     }
0130   }
0131
0132   return mInteraction;
0133 }
0134
0135 }  // namespace ActsPlugins