Geant4/src/Geant4Simulation.cpp

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