Geant4/src/MaterialSteppingAction.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/MaterialSteppingAction.hpp"
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Definitions/Units.hpp"
0013 #include "Acts/Material/Material.hpp"
0014 #include "Acts/Material/MaterialInteraction.hpp"
0015 #include "Acts/Material/MaterialSlab.hpp"
0016 #include "ActsExamples/Geant4/AlgebraConverters.hpp"
0017 #include "ActsExamples/Geant4/EventStore.hpp"
0018
0019 #include <cstddef>
0020 #include <ostream>
0021 #include <unordered_map>
0022 #include <utility>
0023
0024 #include <G4Material.hh>
0025 #include <G4RunManager.hh>
0026 #include <G4Step.hh>
0027
0028 namespace ActsExamples::Geant4 {
0029
0030 MaterialSteppingAction::MaterialSteppingAction(
0031     const Config& cfg, std::unique_ptr<const Acts::Logger> logger)
0032     : G4UserSteppingAction(), m_cfg(cfg), m_logger(std::move(logger)) {}
0033
0034 MaterialSteppingAction::~MaterialSteppingAction() = default;
0035
0036 void MaterialSteppingAction::UserSteppingAction(const G4Step* step) {
0037   // Get the material & check if it is present
0038   G4Material* material = step->GetPreStepPoint()->GetMaterial();
0039   if (material == nullptr) {
0040     return;
0041   }
0042
0043   // First check for exclusion
0044   std::string materialName = material->GetName();
0045   for (const auto& emat : m_cfg.excludeMaterials) {
0046     if (emat == materialName) {
0047       ACTS_VERBOSE("Exclude step in material '" << materialName << "'.");
0048       return;
0049     }
0050   }
0051
0052   constexpr double convertLength = Acts::UnitConstants::mm / CLHEP::mm;
0053   constexpr double convertDensity =
0054       (Acts::UnitConstants::g / Acts::UnitConstants::mm3) /
0055       (CLHEP::gram / CLHEP::mm3);
0056
0057   ACTS_VERBOSE("Performing a step with step size = "
0058                << convertLength * step->GetStepLength());
0059
0060   // Quantities valid for elemental materials and mixtures
0061   double X0 = convertLength * material->GetRadlen();
0062   double L0 = convertLength * material->GetNuclearInterLength();
0063   double rho = convertDensity * material->GetDensity();
0064
0065   // Get{A,Z} is only meaningful for single-element materials (according to
0066   // the Geant4 docs). Need to compute average manually.
0067   const G4ElementVector* elements = material->GetElementVector();
0068   const G4double* fraction = material->GetFractionVector();
0069   std::size_t nElements = material->GetNumberOfElements();
0070   double Ar = 0.;
0071   double Z = 0.;
0072   if (nElements == 1) {
0073     Ar = material->GetA() / (CLHEP::gram / CLHEP::mole);
0074     Z = material->GetZ();
0075   } else {
0076     for (std::size_t i = 0; i < nElements; i++) {
0077       Ar += elements->at(i)->GetA() * fraction[i] / (CLHEP::gram / CLHEP::mole);
0078       Z += elements->at(i)->GetZ() * fraction[i];
0079     }
0080   }
0081
0082   // Construct passed material slab for the step
0083   const auto slab =
0084       Acts::MaterialSlab(Acts::Material::fromMassDensity(X0, L0, Ar, Z, rho),
0085                          convertLength * step->GetStepLength());
0086
0087   // Create the RecordedMaterialSlab
0088   Acts::MaterialInteraction mInteraction;
0089   mInteraction.position =
0090       convertPosition(step->GetPreStepPoint()->GetPosition());
0091   mInteraction.direction =
0092       convertDirection(step->GetPreStepPoint()->GetMomentum()).normalized();
0093   mInteraction.materialSlab = slab;
0094   mInteraction.pathCorrection = (step->GetStepLength() / CLHEP::mm);
0095
0096   G4Track* g4Track = step->GetTrack();
0097   std::size_t trackID = g4Track->GetTrackID();
0098   auto& materialTracks = eventStore().materialTracks;
0099   if (!materialTracks.contains(trackID - 1)) {
0100     Acts::RecordedMaterialTrack rmTrack;
0101     Acts::Vector3 vertex = convertPosition(g4Track->GetVertexPosition());
0102     Acts::Vector3 direction = convertDirection(g4Track->GetMomentumDirection());
0103     rmTrack.first = {vertex, direction};
0104     rmTrack.second.materialInteractions.push_back(mInteraction);
0105     materialTracks[trackID - 1] = rmTrack;
0106   } else {
0107     materialTracks[trackID - 1].second.materialInteractions.push_back(
0108         mInteraction);
0109   }
0110 }
0111
0112 }  // namespace ActsExamples::Geant4