Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2026-09-19 08:28:56

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/Geant4Simulation.hpp"
0010 
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Definitions/Units.hpp"
0013 #include "Acts/Utilities/Logger.hpp"
0014 #include "ActsExamples/Framework/AlgorithmContext.hpp"
0015 #include "ActsExamples/Framework/IAlgorithm.hpp"
0016 #include "ActsExamples/Framework/RandomNumbers.hpp"
0017 #include "ActsExamples/Framework/WhiteBoard.hpp"
0018 #include "ActsExamples/Geant4/EventStore.hpp"
0019 #include "ActsExamples/Geant4/Geant4Manager.hpp"
0020 #include "ActsExamples/Geant4/MagneticFieldWrapper.hpp"
0021 #include "ActsExamples/Geant4/MaterialPhysicsList.hpp"
0022 #include "ActsExamples/Geant4/MaterialSteppingAction.hpp"
0023 #include "ActsExamples/Geant4/ParticleKillAction.hpp"
0024 #include "ActsExamples/Geant4/ParticleTrackingAction.hpp"
0025 #include "ActsExamples/Geant4/SensitiveSteppingAction.hpp"
0026 #include "ActsExamples/Geant4/SensitiveSurfaceMapper.hpp"
0027 #include "ActsExamples/Geant4/SimParticleTranslation.hpp"
0028 #include "ActsExamples/Geant4/SteppingActionList.hpp"
0029 #include "ActsPlugins/FpeMonitoring/FpeMonitor.hpp"
0030 
0031 #include <stdexcept>
0032 #include <utility>
0033 
0034 #include <G4FieldManager.hh>
0035 #include <G4PropagatorInField.hh>
0036 #include <G4RunManager.hh>
0037 #include <G4TransportationManager.hh>
0038 #include <G4UniformMagField.hh>
0039 #include <G4UserEventAction.hh>
0040 #include <G4UserLimits.hh>
0041 #include <G4UserRunAction.hh>
0042 #include <G4UserSteppingAction.hh>
0043 #include <G4UserTrackingAction.hh>
0044 #include <G4VUserDetectorConstruction.hh>
0045 #include <G4VUserPhysicsList.hh>
0046 #include <G4Version.hh>
0047 #include <Randomize.hh>
0048 
0049 namespace ActsExamples {
0050 
0051 Geant4SimulationBase::Geant4SimulationBase(
0052     const Config& cfg, const std::string& name,
0053     std::unique_ptr<const Acts::Logger> logger)
0054     : IAlgorithm(name, std::move(logger)) {
0055   if (cfg.inputParticles.empty()) {
0056     throw std::invalid_argument("Missing input particle collection");
0057   }
0058   if (cfg.detector == nullptr) {
0059     throw std::invalid_argument("Missing detector construction factory");
0060   }
0061   if (cfg.randomNumbers == nullptr) {
0062     throw std::invalid_argument("Missing random numbers");
0063   }
0064 
0065   m_eventStore = std::make_shared<Geant4::EventStore>();
0066 
0067   // tweak logging
0068   // If we are in VERBOSE mode, set the verbose level in Geant4 to 2.
0069   // 3 would be also possible, but that produces infinite amount of output.
0070   m_geant4Level = this->logger().level() == Acts::Logging::VERBOSE ? 2 : 0;
0071 }
0072 
0073 Geant4SimulationBase::~Geant4SimulationBase() = default;
0074 
0075 void Geant4SimulationBase::commonInitialization() {
0076   // Set the detector construction
0077   {
0078     // Clear detector construction if it exists
0079     if (runManager().GetUserDetectorConstruction() != nullptr) {
0080       delete runManager().GetUserDetectorConstruction();
0081     }
0082     // G4RunManager will take care of deletion
0083     m_detectorConstruction =
0084         config()
0085             .detector
0086             ->buildGeant4DetectorConstruction(config().constructionOptions)
0087             .release();
0088     runManager().SetUserInitialization(m_detectorConstruction);
0089     runManager().InitializeGeometry();
0090   }
0091 
0092   m_geant4Instance->tweakLogging(m_geant4Level);
0093 }
0094 
0095 G4RunManager& Geant4SimulationBase::runManager() const {
0096   return *m_geant4Instance->runManager;
0097 }
0098 
0099 Geant4::EventStore& Geant4SimulationBase::eventStore() const {
0100   return *m_eventStore;
0101 }
0102 
0103 ProcessCode Geant4SimulationBase::initialize() {
0104   // Initialize the Geant4 run manager
0105   runManager().Initialize();
0106 
0107   return ProcessCode::SUCCESS;
0108 }
0109 
0110 ProcessCode Geant4SimulationBase::execute(const AlgorithmContext& ctx) const {
0111   // Ensure exclusive access to the Geant4 run manager
0112   std::lock_guard<std::mutex> guard(m_geant4Instance->mutex);
0113 
0114   // If requested cap the max propagator in field to avoid initial stepping
0115   // errors in case volumes are larger than field map
0116 
0117   if (std::isfinite(config().propagatorLargestAcceptableStep)) {
0118     G4PropagatorInField* propagator =
0119         G4TransportationManager::GetTransportationManager()
0120             ->GetPropagatorInField();
0121     propagator->SetLargestAcceptableStep(
0122         config().propagatorLargestAcceptableStep / Acts::UnitConstants::mm *
0123         CLHEP::mm);
0124   }
0125 
0126   // Set the seed new per event, so that we get reproducible results
0127   G4Random::setTheSeed(config().randomNumbers->generateSeed(ctx));
0128 
0129   // Get and reset event registry state
0130   eventStore() = Geant4::EventStore{};
0131 
0132   // Register the current event store to the registry
0133   // this will allow access from the User*Actions
0134   eventStore().store = &(ctx.eventStore);
0135 
0136   // Register the input particle read handle
0137   eventStore().inputParticles = &m_inputParticles;
0138 
0139   // Only reaches the user actions. The G4 geometry itself is not context aware,
0140   // so use Fatras to simulate a misaligned detector.
0141   eventStore().geoContext = ctx.simGeoContext;
0142 
0143   ACTS_DEBUG("Sending Geant RunManager the BeamOn() command.");
0144   {
0145     ActsPlugins::FpeMonitor mon{0};  // disable all FPEs while we're in Geant4
0146     // Start simulation. each track is simulated as a separate Geant4 event.
0147     runManager().BeamOn(1);
0148   }
0149 
0150   // Print out warnings about possible particle collision if happened
0151   if (eventStore().particleIdCollisionsInitial > 0 ||
0152       eventStore().particleIdCollisionsFinal > 0 ||
0153       eventStore().parentIdNotFound > 0) {
0154     ACTS_WARNING(
0155         "Particle ID collisions detected, don't trust the particle "
0156         "identification!");
0157     ACTS_WARNING(
0158         "- initial states: " << eventStore().particleIdCollisionsInitial);
0159     ACTS_WARNING("- final states: " << eventStore().particleIdCollisionsFinal);
0160     ACTS_WARNING("- parent ID not found: " << eventStore().parentIdNotFound);
0161   }
0162 
0163   if (eventStore().hits.empty()) {
0164     ACTS_DEBUG("Step merging: No steps recorded");
0165   } else {
0166     ACTS_DEBUG("Step merging: mean hits per hit: "
0167                << static_cast<double>(eventStore().numberGeantSteps) /
0168                       eventStore().hits.size());
0169     ACTS_DEBUG(
0170         "Step merging: max hits per hit: " << eventStore().maxStepsForHit);
0171   }
0172 
0173   return ProcessCode::SUCCESS;
0174 }
0175 
0176 std::shared_ptr<Geant4Handle> Geant4SimulationBase::geant4Handle() const {
0177   return m_geant4Instance;
0178 }
0179 
0180 Geant4Simulation::Geant4Simulation(const Config& cfg,
0181                                    std::unique_ptr<const Acts::Logger> logger)
0182     : Geant4SimulationBase(cfg, "Geant4Simulation", std::move(logger)),
0183       m_cfg(cfg) {
0184   m_geant4Instance =
0185       m_cfg.geant4Handle
0186           ? m_cfg.geant4Handle
0187           : Geant4Manager::instance().createHandle(m_cfg.physicsList);
0188   if (m_geant4Instance->physicsListName != m_cfg.physicsList) {
0189     throw std::runtime_error("inconsistent physics list");
0190   }
0191 
0192   commonInitialization();
0193 
0194   // Set the primarty generator
0195   {
0196     // Clear primary generation action if it exists
0197     if (runManager().GetUserPrimaryGeneratorAction() != nullptr) {
0198       delete runManager().GetUserPrimaryGeneratorAction();
0199     }
0200     Geant4::SimParticleTranslation::Config prCfg;
0201     prCfg.eventStore = m_eventStore;
0202     // G4RunManager will take care of deletion
0203     auto primaryGeneratorAction = new Geant4::SimParticleTranslation(
0204         prCfg, this->logger().cloneWithSuffix("SimParticleTranslation"));
0205     // Set the primary generator action
0206     runManager().SetUserAction(primaryGeneratorAction);
0207   }
0208 
0209   // Particle action
0210   {
0211     // Clear tracking action if it exists
0212     if (runManager().GetUserTrackingAction() != nullptr) {
0213       delete runManager().GetUserTrackingAction();
0214     }
0215     Geant4::ParticleTrackingAction::Config trackingCfg;
0216     trackingCfg.eventStore = m_eventStore;
0217     trackingCfg.keepParticlesWithoutHits = cfg.keepParticlesWithoutHits;
0218     // G4RunManager will take care of deletion
0219     auto trackingAction = new Geant4::ParticleTrackingAction(
0220         trackingCfg, this->logger().cloneWithSuffix("ParticleTracking"));
0221     runManager().SetUserAction(trackingAction);
0222   }
0223 
0224   // Stepping actions
0225   Geant4::SensitiveSteppingAction* sensitiveSteppingActionAccess = nullptr;
0226   {
0227     // Clear stepping action if it exists
0228     if (runManager().GetUserSteppingAction() != nullptr) {
0229       delete runManager().GetUserSteppingAction();
0230     }
0231 
0232     Geant4::ParticleKillAction::Config particleKillCfg;
0233     particleKillCfg.eventStore = m_eventStore;
0234     particleKillCfg.volume = cfg.killVolume;
0235     particleKillCfg.maxTime = cfg.killAfterTime;
0236     particleKillCfg.secondaries = cfg.killSecondaries;
0237 
0238     Geant4::SensitiveSteppingAction::Config stepCfg;
0239     stepCfg.eventStore = m_eventStore;
0240     stepCfg.charged = cfg.recordHitsOfCharged;
0241     stepCfg.neutral = cfg.recordHitsOfNeutrals;
0242     stepCfg.primary = cfg.recordHitsOfPrimaries;
0243     stepCfg.secondary = cfg.recordHitsOfSecondaries;
0244     stepCfg.stepLogging = cfg.recordPropagationSummaries;
0245 
0246     Geant4::SteppingActionList::Config steppingCfg;
0247     steppingCfg.actions.push_back(std::make_unique<Geant4::ParticleKillAction>(
0248         particleKillCfg, this->logger().cloneWithSuffix("Killer")));
0249 
0250     auto sensitiveSteppingAction =
0251         std::make_unique<Geant4::SensitiveSteppingAction>(
0252             stepCfg, this->logger().cloneWithSuffix("SensitiveStepping"));
0253     sensitiveSteppingActionAccess = sensitiveSteppingAction.get();
0254 
0255     steppingCfg.actions.push_back(std::move(sensitiveSteppingAction));
0256 
0257     // G4RunManager will take care of deletion
0258     auto steppingAction = new Geant4::SteppingActionList(steppingCfg);
0259     runManager().SetUserAction(steppingAction);
0260   }
0261 
0262   // Get the g4World cache
0263   G4VPhysicalVolume* g4World = m_detectorConstruction->Construct();
0264 
0265   // Please note:
0266   // The following two blocks rely on the fact that the Acts
0267   // detector constructions cache the world volume
0268 
0269   // Set the magnetic field
0270   if (cfg.magneticField) {
0271     ACTS_LOG_WITH_LOGGER(this->logger(), Acts::Logging::INFO,
0272                          "Setting ACTS configured field to Geant4.");
0273 
0274     Geant4::MagneticFieldWrapper::Config g4FieldCfg;
0275     g4FieldCfg.magneticField = cfg.magneticField;
0276     m_magneticField =
0277         std::make_unique<Geant4::MagneticFieldWrapper>(g4FieldCfg);
0278 
0279     // Set the field or the G4Field manager
0280     m_fieldManager = std::make_unique<G4FieldManager>();
0281     m_fieldManager->SetDetectorField(m_magneticField.get());
0282     m_fieldManager->CreateChordFinder(m_magneticField.get());
0283 
0284     // Propagate down to all childrend
0285     g4World->GetLogicalVolume()->SetFieldManager(m_fieldManager.get(), true);
0286   }
0287 
0288   // ACTS sensitive surfaces are provided, so hit creation is turned on
0289   if (cfg.sensitiveSurfaceMapper != nullptr) {
0290     Geant4::SensitiveSurfaceMapper::State sState;
0291     ACTS_LOG_WITH_LOGGER(this->logger(), Acts::Logging::INFO,
0292                          "Remapping selected volumes from Geant4 to "
0293                          "Acts::Surface::GeometryID");
0294     cfg.sensitiveSurfaceMapper->remapSensitiveNames(
0295         sState, Acts::GeometryContext::dangerouslyDefaultConstruct(), g4World,
0296         Acts::Transform3::Identity());
0297 
0298     auto allSurfacesMapped = cfg.sensitiveSurfaceMapper->checkMapping(
0299         sState, Acts::GeometryContext::dangerouslyDefaultConstruct(), false,
0300         false);
0301     if (!allSurfacesMapped) {
0302       ACTS_LOG_WITH_LOGGER(this->logger(), Acts::Logging::WARNING,
0303                            "Not all sensitive surfaces have been mapped to "
0304                            "Geant4 volumes!");
0305     }
0306 
0307     sensitiveSteppingActionAccess->assignSurfaceMapping(
0308         sState.g4VolumeToSurfaces);
0309   }
0310 
0311   m_inputParticles.initialize(cfg.inputParticles);
0312   m_outputSimHits.initialize(cfg.outputSimHits);
0313   m_outputParticles.initialize(cfg.outputParticles);
0314 
0315   if (cfg.recordPropagationSummaries) {
0316     m_outputPropagationSummaries.initialize(cfg.outputPropagationSummaries);
0317   }
0318 }
0319 
0320 Geant4Simulation::~Geant4Simulation() = default;
0321 
0322 ProcessCode Geant4Simulation::execute(const AlgorithmContext& ctx) const {
0323   auto ret = Geant4SimulationBase::execute(ctx);
0324   if (ret != ProcessCode::SUCCESS) {
0325     return ret;
0326   }
0327 
0328   // Output handling: Simulation
0329   m_outputParticles(
0330       ctx, SimParticleContainer(eventStore().particlesSimulated.begin(),
0331                                 eventStore().particlesSimulated.end()));
0332 
0333   m_outputSimHits(
0334       ctx, SimHitContainer(eventStore().hits.begin(), eventStore().hits.end()));
0335 
0336   // Output the propagation summaries if requested
0337   if (m_cfg.recordPropagationSummaries) {
0338     PropagationSummaries summaries;
0339     summaries.reserve(eventStore().propagationRecords.size());
0340     for (auto& [trackId, summary] : eventStore().propagationRecords) {
0341       summaries.push_back(std::move(summary));
0342     }
0343     m_outputPropagationSummaries(ctx, std::move(summaries));
0344   }
0345 
0346   return ProcessCode::SUCCESS;
0347 }
0348 
0349 Geant4MaterialRecording::Geant4MaterialRecording(
0350     const Config& cfg, std::unique_ptr<const Acts::Logger> logger)
0351     : Geant4SimulationBase(cfg, "Geant4Simulation", std::move(logger)),
0352       m_cfg(cfg) {
0353   auto physicsListName = "MaterialPhysicsList";
0354   m_geant4Instance =
0355       m_cfg.geant4Handle
0356           ? m_cfg.geant4Handle
0357           : Geant4Manager::instance().createHandle(
0358                 std::make_unique<Geant4::MaterialPhysicsList>(
0359                     this->logger().cloneWithSuffix("MaterialPhysicsList")),
0360                 physicsListName);
0361   if (m_geant4Instance->physicsListName != physicsListName) {
0362     throw std::runtime_error("inconsistent physics list");
0363   }
0364 
0365   commonInitialization();
0366 
0367   // Set the primarty generator
0368   {
0369     // Clear primary generation action if it exists
0370     if (runManager().GetUserPrimaryGeneratorAction() != nullptr) {
0371       delete runManager().GetUserPrimaryGeneratorAction();
0372     }
0373 
0374     Geant4::SimParticleTranslation::Config prCfg;
0375     prCfg.eventStore = m_eventStore;
0376     prCfg.forcedPdgCode = 0;
0377     prCfg.forcedCharge = 0.;
0378     prCfg.forcedMass = 0.;
0379 
0380     // G4RunManager will take care of deletion
0381     auto primaryGeneratorAction = new Geant4::SimParticleTranslation(
0382         prCfg, this->logger().cloneWithSuffix("SimParticleTranslation"));
0383     // Set the primary generator action
0384     runManager().SetUserAction(primaryGeneratorAction);
0385   }
0386 
0387   // Particle action
0388   {
0389     // Clear tracking action if it exists
0390     if (runManager().GetUserTrackingAction() != nullptr) {
0391       delete runManager().GetUserTrackingAction();
0392     }
0393     Geant4::ParticleTrackingAction::Config trackingCfg;
0394     trackingCfg.eventStore = m_eventStore;
0395     trackingCfg.keepParticlesWithoutHits = true;
0396     // G4RunManager will take care of deletion
0397     auto trackingAction = new Geant4::ParticleTrackingAction(
0398         trackingCfg, this->logger().cloneWithSuffix("ParticleTracking"));
0399     runManager().SetUserAction(trackingAction);
0400   }
0401 
0402   // Stepping action
0403   {
0404     // Clear stepping action if it exists
0405     if (runManager().GetUserSteppingAction() != nullptr) {
0406       delete runManager().GetUserSteppingAction();
0407     }
0408     Geant4::MaterialSteppingAction::Config steppingCfg;
0409     steppingCfg.eventStore = m_eventStore;
0410     steppingCfg.excludeMaterials = m_cfg.excludeMaterials;
0411     steppingCfg.recordElementFractions = m_cfg.recordElementFractions;
0412     // G4RunManager will take care of deletion
0413     auto steppingAction = new Geant4::MaterialSteppingAction(
0414         steppingCfg, this->logger().cloneWithSuffix("MaterialSteppingAction"));
0415     runManager().SetUserAction(steppingAction);
0416   }
0417 
0418   runManager().Initialize();
0419 
0420   m_inputParticles.initialize(cfg.inputParticles);
0421   m_outputMaterialTracks.initialize(cfg.outputMaterialTracks);
0422 }
0423 
0424 Geant4MaterialRecording::~Geant4MaterialRecording() = default;
0425 
0426 ProcessCode Geant4MaterialRecording::execute(
0427     const AlgorithmContext& ctx) const {
0428   const auto ret = Geant4SimulationBase::execute(ctx);
0429   if (ret != ProcessCode::SUCCESS) {
0430     return ret;
0431   }
0432 
0433   // Output handling: Material tracks
0434   m_outputMaterialTracks(
0435       ctx, decltype(eventStore().materialTracks)(eventStore().materialTracks));
0436 
0437   return ProcessCode::SUCCESS;
0438 }
0439 
0440 }  // namespace ActsExamples