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/Material/Material.hpp"
0013 #include "Acts/Material/MaterialInteraction.hpp"
0014 #include "Acts/Material/MaterialSlab.hpp"
0015 #include "ActsExamples/Geant4/AlgebraConverters.hpp"
0016 #include "ActsExamples/Geant4/EventStore.hpp"
0017 #include "ActsExamples/Geant4/UnitConversion.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 void MaterialSteppingAction::UserSteppingAction(const G4Step* stepPtr) {
0035   assert(stepPtr != nullptr);
0036   const G4Step& step = *stepPtr;
0037
0038   // Get the material & check if it is present
0039   const G4Material* materialPtr = step.GetPreStepPoint()->GetMaterial();
0040   if (materialPtr == nullptr) {
0041     return;
0042   }
0043   const G4Material& material = *materialPtr;
0044
0045   // First check for exclusion
0046   const std::string materialName = material.GetName();
0047   for (const std::string& emat : m_cfg.excludeMaterials) {
0048     if (emat == materialName) {
0049       ACTS_VERBOSE("Exclude step in material '" << materialName << "'.");
0050       return;
0051     }
0052   }
0053
0054   ACTS_VERBOSE("Performing a step with step size = "
0055                << convertLengthToActs * step.GetStepLength());
0056
0057   // Quantities valid for elemental materials and mixtures
0058   const double X0 = convertLengthToActs * material.GetRadlen();
0059   const double L0 = convertLengthToActs * material.GetNuclearInterLength();
0060   const double rho = convertDensityToActs * material.GetDensity();
0061
0062   // Get{A,Z} is only meaningful for single-element materials (according to
0063   // the Geant4 docs). Need to compute average manually.
0064   const G4ElementVector* elements = material.GetElementVector();
0065   const G4double* fraction = material.GetFractionVector();
0066   const std::size_t nElements = material.GetNumberOfElements();
0067   double Ar = 0.;
0068   double Z = 0.;
0069   if (nElements == 1) {
0070     Ar = material.GetA() * (CLHEP::mole / CLHEP::gram);
0071     Z = material.GetZ();
0072   } else {
0073     for (std::size_t i = 0; i < nElements; i++) {
0074       Ar += elements->at(i)->GetA() * fraction[i] * (CLHEP::mole / CLHEP::gram);
0075       Z += elements->at(i)->GetZ() * fraction[i];
0076     }
0077   }
0078
0079   // Construct passed material slab for the step
0080   const Acts::MaterialSlab slab(
0081       Acts::Material::fromMassDensity(X0, L0, Ar, Z, rho),
0082       convertLengthToActs * step.GetStepLength());
0083
0084   // Create the RecordedMaterialSlab
0085   Acts::MaterialInteraction mInteraction;
0086   mInteraction.position =
0087       convertPosition(step.GetPreStepPoint()->GetPosition());
0088   mInteraction.direction =
0089       convertDirection(step.GetPreStepPoint()->GetMomentum()).normalized();
0090   mInteraction.materialSlab = slab;
0091   mInteraction.pathCorrection = step.GetStepLength() * convertLengthToActs;
0092
0093   if (m_cfg.recordElementFractions) {
0094     if (nElements == 1) {
0095       mInteraction.elementZ.push_back(
0096           static_cast<unsigned int>(material.GetZ()));
0097       mInteraction.elementFrac.push_back(1.0f);
0098     } else {
0099       for (std::size_t i = 0; i < nElements; i++) {
0100         mInteraction.elementZ.push_back(
0101             static_cast<unsigned int>(elements->at(i)->GetZ()));
0102         mInteraction.elementFrac.push_back(static_cast<float>(fraction[i]));
0103       }
0104     }
0105   }
0106
0107   assert(step.GetTrack() != nullptr);
0108   const G4Track& track = *step.GetTrack();
0109   const std::size_t trackId = track.GetTrackID();
0110   if (!eventStore().materialTracks.contains(trackId - 1)) {
0111     Acts::RecordedMaterialTrack rmTrack;
0112     const Acts::Vector3 vertex = convertPosition(track.GetVertexPosition());
0113     const Acts::Vector3 direction =
0114         convertDirection(track.GetMomentumDirection());
0115     rmTrack.first = {vertex, direction};
0116     rmTrack.second.materialInteractions.push_back(mInteraction);
0117     eventStore().materialTracks[trackId - 1] = rmTrack;
0118   } else {
0119     eventStore()
0120         .materialTracks[trackId - 1]
0121         .second.materialInteractions.push_back(mInteraction);
0122   }
0123 }
0124
0125 }  // namespace ActsExamples::Geant4