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