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