Geant4/src/SensitiveSteppingAction.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/SensitiveSteppingAction.hpp"
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Definitions/Units.hpp"
0013 #include "Acts/Geometry/GeometryIdentifier.hpp"
0014 #include "Acts/Propagator/detail/SteppingLogger.hpp"
0015 #include "Acts/Surfaces/Surface.hpp"
0016 #include "Acts/Utilities/MultiIndex.hpp"
0017 #include "ActsExamples/Geant4/EventStore.hpp"
0018 #include "ActsExamples/Geant4/SensitiveSurfaceMapper.hpp"
0019 #include "ActsFatras/EventData/Barcode.hpp"
0020
0021 #include <algorithm>
0022 #include <array>
0023 #include <cstddef>
0024 #include <string>
0025 #include <unordered_map>
0026 #include <utility>
0027
0028 #include <G4RunManager.hh>
0029 #include <G4Step.hh>
0030 #include <G4StepPoint.hh>
0031 #include <G4Track.hh>
0032 #include <G4UnitsTable.hh>
0033 #include <G4VPhysicalVolume.hh>
0034 #include <G4VTouchable.hh>
0035 #include <boost/version.hpp>
0036
0037 #if BOOST_VERSION >= 107800
0038 #include <boost/describe.hpp>
0039
0040 BOOST_DESCRIBE_ENUM(G4StepStatus, fWorldBoundary, fGeomBoundary,
0041                     fAtRestDoItProc, fAlongStepDoItProc, fPostStepDoItProc,
0042                     fUserDefinedLimit, fExclusivelyForcedProc, fUndefined);
0043
0044 BOOST_DESCRIBE_ENUM(G4ProcessType, fNotDefined, fTransportation,
0045                     fElectromagnetic, fOptical, fHadronic, fPhotolepton_hadron,
0046                     fDecay, fGeneral, fParameterisation, fUserDefined,
0047                     fParallel, fPhonon, fUCN);
0048
0049 BOOST_DESCRIBE_ENUM(G4TrackStatus, fAlive, fStopButAlive, fStopAndKill,
0050                     fKillTrackAndSecondaries, fSuspend, fPostponeToNextEvent);
0051 #endif
0052
0053 namespace {
0054
0055 std::array<Acts::Vector4, 4u> kinematicsOfStep(const G4Step* step) {
0056   static constexpr double convertLength = Acts::UnitConstants::mm / CLHEP::mm;
0057   static constexpr double convertEnergy = Acts::UnitConstants::GeV / CLHEP::GeV;
0058   static constexpr double convertTime = Acts::UnitConstants::ns / CLHEP::ns;
0059
0060   const G4StepPoint* preStepPoint = step->GetPreStepPoint();
0061   const G4StepPoint* postStepPoint = step->GetPostStepPoint();
0062
0063   Acts::Vector4 preStepPosition(convertLength * preStepPoint->GetPosition().x(),
0064                                 convertLength * preStepPoint->GetPosition().y(),
0065                                 convertLength * preStepPoint->GetPosition().z(),
0066                                 convertTime * preStepPoint->GetGlobalTime());
0067   Acts::Vector4 preStepMomentum(convertEnergy * preStepPoint->GetMomentum().x(),
0068                                 convertEnergy * preStepPoint->GetMomentum().y(),
0069                                 convertEnergy * preStepPoint->GetMomentum().z(),
0070                                 convertEnergy * preStepPoint->GetTotalEnergy());
0071   Acts::Vector4 postStepPosition(
0072       convertLength * postStepPoint->GetPosition().x(),
0073       convertLength * postStepPoint->GetPosition().y(),
0074       convertLength * postStepPoint->GetPosition().z(),
0075       convertTime * postStepPoint->GetGlobalTime());
0076   Acts::Vector4 postStepMomentum(
0077       convertEnergy * postStepPoint->GetMomentum().x(),
0078       convertEnergy * postStepPoint->GetMomentum().y(),
0079       convertEnergy * postStepPoint->GetMomentum().z(),
0080       convertEnergy * postStepPoint->GetTotalEnergy());
0081
0082   return {preStepPosition, preStepMomentum, postStepPosition, postStepMomentum};
0083 }
0084
0085 ActsFatras::Hit hitFromStep(const G4Step* step, ActsFatras::Barcode particleId,
0086                             Acts::GeometryIdentifier geoId,
0087                             std::int32_t index) {
0088   auto [preStepPosition, preStepMomentum, postStepPosition, postStepMomentum] =
0089       kinematicsOfStep(step);
0090
0091   return ActsFatras::Hit(geoId, particleId,
0092                          0.5 * (preStepPosition + postStepPosition),
0093                          preStepMomentum, postStepMomentum, index);
0094 }
0095
0096 Acts::detail::Step stepFromG4Step(const G4Step* step) {
0097   Acts::detail::Step pStep;
0098   auto [preStepPosition, preStepMomentum, postStepPosition, postStepMomentum] =
0099       kinematicsOfStep(step);
0100
0101   pStep.navDir = Acts::Direction::Forward();
0102   pStep.position = 0.5 * (preStepPosition + postStepPosition).block<3, 1>(0, 0);
0103   pStep.momentum = 0.5 * (preStepMomentum + postStepMomentum).block<3, 1>(0, 0);
0104   pStep.nTotalTrials = 1;
0105   return pStep;
0106 }
0107
0108 }  // namespace
0109
0110 namespace ActsExamples::Geant4 {
0111
0112 SensitiveSteppingAction::SensitiveSteppingAction(
0113     const Config& cfg, std::unique_ptr<const Acts::Logger> logger)
0114     : G4UserSteppingAction(), m_cfg(cfg), m_logger(std::move(logger)) {}
0115
0116 void SensitiveSteppingAction::UserSteppingAction(const G4Step* step) {
0117   // Unit conversions G4->::ACTS
0118   static constexpr double convertLength = Acts::UnitConstants::mm / CLHEP::mm;
0119   static constexpr double convertEnergy = Acts::UnitConstants::GeV / CLHEP::GeV;
0120
0121   // The particle after the step
0122   G4Track* track = step->GetTrack();
0123   G4PrimaryParticle* primaryParticle =
0124       track->GetDynamicParticle()->GetPrimaryParticle();
0125
0126   // Get PreStepPoint and PostStepPoint
0127   const G4StepPoint* preStepPoint = step->GetPreStepPoint();
0128   const G4StepPoint* postStepPoint = step->GetPostStepPoint();
0129
0130   // Bail out if charged & configured to do so
0131   G4double absCharge = std::abs(track->GetParticleDefinition()->GetPDGCharge());
0132   if (!m_cfg.charged && absCharge > 0.) {
0133     return;
0134   }
0135
0136   // Bail out if neutral & configured to do so
0137   if (!m_cfg.neutral && absCharge == 0.) {
0138     return;
0139   }
0140
0141   // Bail out if it is a primary & configured to be ignored
0142   if (!m_cfg.primary && primaryParticle != nullptr) {
0143     return;
0144   }
0145
0146   // Bail out if it is a secondary & configured to be ignored
0147   if (!m_cfg.secondary && primaryParticle == nullptr) {
0148     return;
0149   }
0150
0151   // Get the physical volume & check if it has the sensitive string name
0152   const G4VPhysicalVolume* volume = track->GetVolume();
0153   if (volume == nullptr) {
0154     throw std::runtime_error("No volume found, terminate simulation");
0155   }
0156   std::string volumeName = volume->GetName();
0157
0158   if (volumeName.find(SensitiveSurfaceMapper::mappingPrefix) ==
0159           std::string::npos &&
0160       !m_cfg.stepLogging) {
0161     return;
0162   }
0163
0164   // The G4Touchable for the matching
0165   const G4VTouchable* touchable = track->GetTouchable();
0166
0167   Acts::GeometryIdentifier geoId{};
0168
0169   // Find the range of candidate surfaces for the current position in the
0170   // mapping multimap
0171   auto [bsf, esf] = m_surfaceMapping.equal_range(volume);
0172   std::size_t nSurfaces = std::distance(bsf, esf);
0173
0174   ACTS_VERBOSE("Found " << nSurfaces << " candidate surfaces for volume "
0175                         << volumeName);
0176
0177   const Acts::Surface* surface = nullptr;
0178   if (nSurfaces == 0 && !m_cfg.stepLogging) {
0179     ACTS_ERROR("No candidate surfaces found for volume " << volumeName);
0180     return;
0181   } else if (nSurfaces == 1u) {
0182     geoId = bsf->second->geometryId();
0183     ACTS_VERBOSE("Unique assignment successful -> to surface " << geoId);
0184   } else {
0185     // Find the closest surface to the current position
0186     Acts::GeometryContext gctx;
0187     for (; bsf != esf; ++bsf) {
0188       surface = bsf->second;
0189       const G4ThreeVector& translation = touchable->GetTranslation();
0190       Acts::Vector3 g4VolumePosition(convertLength * translation.x(),
0191                                      convertLength * translation.y(),
0192                                      convertLength * translation.z());
0193       if (surface->center(gctx).isApprox(g4VolumePosition)) {
0194         geoId = surface->geometryId();
0195         break;
0196       }
0197     }
0198     ACTS_VERBOSE("Replica assignment successful -> to surface " << geoId);
0199   }
0200
0201   // This is not the case if we have a particle-ID collision
0202   if (!eventStore().trackIdMapping.contains(track->GetTrackID())) {
0203     return;
0204   }
0205
0206   // Output is only strictly valid if step logging is not enabled
0207   const auto particleId = eventStore().trackIdMapping.at(track->GetTrackID());
0208   if (!m_cfg.stepLogging && surface != nullptr) {
0209     ACTS_VERBOSE("Step of " << particleId << " in sensitive volume " << geoId);
0210   } else if (m_cfg.stepLogging) {
0211     if (!eventStore().propagationRecords.contains(track->GetTrackID())) {
0212       // Create the propagation summary
0213       double xVtx = track->GetVertexPosition().x() * convertLength;
0214       double yVtx = track->GetVertexPosition().y() * convertLength;
0215       double zVtx = track->GetVertexPosition().z() * convertLength;
0216       double xDirVtx = track->GetVertexMomentumDirection().x();
0217       double yDirVtx = track->GetVertexMomentumDirection().y();
0218       double zDirVtx = track->GetVertexMomentumDirection().z();
0219       double absMomentum = track->GetMomentum().mag() * convertEnergy;
0220
0221       PropagationSummary iSummary(Acts::CurvilinearTrackParameters(
0222           Acts::Vector4(xVtx, yVtx, zVtx, 0.),
0223           Acts::Vector3(xDirVtx, yDirVtx, zDirVtx), absCharge / absMomentum,
0224           std::nullopt, Acts::ParticleHypothesis::pion()));
0225
0226       eventStore().propagationRecords.insert({track->GetTrackID(), iSummary});
0227     }
0228     PropagationSummary& pSummary =
0229         eventStore().propagationRecords.at(track->GetTrackID());
0230
0231     // Increase the step counter
0232     pSummary.nSteps += 1;
0233
0234     double currentTrackLength = track->GetTrackLength() * convertLength;
0235     double currentStepLength = currentTrackLength - pSummary.pathLength;
0236     pSummary.pathLength = currentTrackLength;
0237
0238     // Create a new step for the step logging
0239     Acts::detail::Step pStep = stepFromG4Step(step);
0240     pStep.geoID = geoId;
0241     pStep.surface = surface != nullptr ? surface->getSharedPtr() : nullptr;
0242     // Check if last step was on same surface
0243     if (!pSummary.steps.empty() && pSummary.steps.back().geoID == geoId &&
0244         pSummary.steps.back().surface != nullptr) {
0245       auto& lastStep = pSummary.steps.back();
0246       lastStep.stepSize = Acts::ConstrainedStep(currentStepLength);
0247       lastStep.position = 0.5 * (pStep.position + lastStep.position);
0248       lastStep.momentum = 0.5 * (pStep.momentum + lastStep.momentum);
0249     } else {
0250       // Record the propagation state
0251       pStep.stepSize = Acts::ConstrainedStep(currentStepLength);
0252       pSummary.steps.emplace_back(std::move(pStep));
0253     }
0254     // You have nothing to do from here
0255     return;
0256   }
0257
0258   // Set particle hit count to zero, so we have this entry in the map later
0259   if (!eventStore().particleHitCount.contains(particleId)) {
0260     eventStore().particleHitCount[particleId] = 0;
0261   }
0262
0263   // Extract if we are at volume boundaries
0264   const bool preOnBoundary = preStepPoint->GetStepStatus() == fGeomBoundary;
0265   const bool postOnBoundary = postStepPoint->GetStepStatus() == fGeomBoundary ||
0266                               postStepPoint->GetStepStatus() == fWorldBoundary;
0267   const bool particleStopped = (postStepPoint->GetKineticEnergy() == 0.0);
0268   const bool particleDecayed =
0269       (postStepPoint->GetProcessDefinedStep()->GetProcessType() == fDecay);
0270
0271   auto print = [](auto s) {
0272 #if BOOST_VERSION >= 107800
0273     return boost::describe::enum_to_string(s, "unmatched");
0274 #else
0275     return s;
0276 #endif
0277   };
0278   ACTS_VERBOSE("status: pre="
0279                << print(preStepPoint->GetStepStatus())
0280                << ", post=" << print(postStepPoint->GetStepStatus())
0281                << ", post E_kin=" << std::boolalpha
0282                << postStepPoint->GetKineticEnergy() << ", process_type="
0283                << print(
0284                       postStepPoint->GetProcessDefinedStep()->GetProcessType())
0285                << ", particle="
0286                << track->GetParticleDefinition()->GetParticleName()
0287                << ", process_name="
0288                << postStepPoint->GetProcessDefinedStep()->GetProcessName()
0289                << ", track status=" << print(track->GetTrackStatus()));
0290
0291   // Case A: The step starts at the entry of the volume and ends at the exit.
0292   // Add hit to collection.
0293   if (preOnBoundary && postOnBoundary) {
0294     ACTS_VERBOSE("-> merge single step to hit");
0295     ++eventStore().particleHitCount[particleId];
0296     eventStore().hits.push_back(
0297         hitFromStep(step, particleId, geoId,
0298                     eventStore().particleHitCount.at(particleId) - 1));
0299
0300     eventStore().numberGeantSteps += 1ul;
0301     eventStore().maxStepsForHit = std::max(eventStore().maxStepsForHit, 1ul);
0302     return;
0303   }
0304
0305   // Case B: The step ends at the exit of the volume. Add hit to hit buffer, and
0306   // then combine hit buffer
0307   if (postOnBoundary || particleStopped || particleDecayed) {
0308     ACTS_VERBOSE("-> merge buffer to hit");
0309     auto& buffer = eventStore().hitBuffer;
0310     buffer.push_back(hitFromStep(step, particleId, geoId, -1));
0311
0312     const auto pos4 =
0313         0.5 * (buffer.front().fourPosition() + buffer.back().fourPosition());
0314
0315     ++eventStore().particleHitCount[particleId];
0316     eventStore().hits.emplace_back(
0317         geoId, particleId, pos4, buffer.front().momentum4Before(),
0318         buffer.back().momentum4After(),
0319         eventStore().particleHitCount.at(particleId) - 1);
0320
0321     assert(std::ranges::all_of(
0322         buffer, [&](const auto& h) { return h.geometryId() == geoId; }));
0323     assert(std::ranges::all_of(
0324         buffer, [&](const auto& h) { return h.particleId() == particleId; }));
0325
0326     eventStore().numberGeantSteps += buffer.size();
0327     eventStore().maxStepsForHit =
0328         std::max(eventStore().maxStepsForHit, buffer.size());
0329
0330     buffer.clear();
0331     return;
0332   }
0333
0334   // Case C: The step doesn't end at the exit of the volume. Add the hit to the
0335   // hit buffer.
0336   if (!postOnBoundary) {
0337     // ACTS_VERBOSE("-> add hit to buffer");
0338     eventStore().hitBuffer.push_back(hitFromStep(step, particleId, geoId, -1));
0339     return;
0340   }
0341
0342   assert(false && "should never reach this");
0343 }
0344
0345 }  // namespace ActsExamples::Geant4