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