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