//========================================================================== // AIDA Detector description implementation //-------------------------------------------------------------------------- // Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN) // All rights reserved. // // For the licensing terms see $DD4hepINSTALL/LICENSE. // For the list of contributors see $DD4hepINSTALL/doc/CREDITS. // // Author : M.Frank // //========================================================================== // Framework include files #include <DDG4/Geant4PhysicsList.h> #include <DDG4/Geant4UIMessenger.h> #include <DDG4/Geant4Particle.h> #include <DDG4/Geant4Kernel.h> #include <DD4hep/InstanceCount.h> #include <DD4hep/Printout.h> #include <DD4hep/Plugins.h> // Geant4 include files #include <G4VPhysicsConstructor.hh> #include <G4PhysListFactory.hh> #include <G4ProcessManager.hh> #include <G4ParticleTable.hh> #include <G4RunManager.hh> #include <G4VProcess.hh> #include <G4Decay.hh> // C/C++ include files #include <stdexcept> #include <regex.h> using namespace std; using namespace dd4hep; using namespace dd4hep::sim; namespace { struct MyPhysics : G4VUserPhysicsList { void AddTransportation() { this->G4VUserPhysicsList::AddTransportation(); } }; struct EmptyPhysics : public G4VModularPhysicsList { EmptyPhysics() {} virtual ~EmptyPhysics() {} }; struct ParticlePhysics : public G4VPhysicsConstructor { Geant4PhysicsListActionSequence* seq; G4VUserPhysicsList* phys; ParticlePhysics(Geant4PhysicsListActionSequence* s, G4VUserPhysicsList* p) : seq(s), phys(p) {} virtual ~ParticlePhysics() {} virtual void ConstructProcess() { seq->constructProcesses(phys); if ( seq->transportation() ) { MyPhysics* ph = (MyPhysics*)phys; ph->AddTransportation(); } } virtual void ConstructParticle() { seq->constructParticles(phys); } }; } /// Default constructor Geant4PhysicsList::Process::Process() : ordAtRestDoIt(-1), ordAlongSteptDoIt(-1), ordPostStepDoIt(-1) { } /// Copy constructor Geant4PhysicsList::Process::Process(const Process& p) : name(p.name), ordAtRestDoIt(p.ordAtRestDoIt), ordAlongSteptDoIt(p.ordAlongSteptDoIt), ordPostStepDoIt(p.ordPostStepDoIt) { } /// Assignment operator Geant4PhysicsList::Process& Geant4PhysicsList::Process::operator=(const Process& p) { if ( this != &p ) { name = p.name; ordAtRestDoIt = p.ordAtRestDoIt; ordAlongSteptDoIt = p.ordAlongSteptDoIt; ordPostStepDoIt = p.ordPostStepDoIt; } return *this; } /// Standard constructor Geant4PhysicsList::Geant4PhysicsList(Geant4Context* ctxt, const string& nam) : Geant4Action(ctxt, nam) { InstanceCount::increment(this); } /// Default destructor Geant4PhysicsList::~Geant4PhysicsList() { InstanceCount::decrement(this); } /// Install command control messenger if wanted void Geant4PhysicsList::installCommandMessenger() { control()->addCall("dump", "Dump content of " + name(), Callback(this).make(&Geant4PhysicsList::dump)); } /// Dump content to stdout void Geant4PhysicsList::dump() { if ( !m_physics.empty() && !m_particles.empty() && !m_processes.empty() ) { printout(ALWAYS,name(),"+++ Geant4PhysicsList Dump"); } for ( const auto& p : m_physics) printout(ALWAYS,name(),"+++ PhysicsConstructor: %s",p.c_str()); for ( const auto& p : m_particles ) printout(ALWAYS,name(),"+++ ParticleConstructor: %s",p.c_str()); for ( const auto& p : m_particlegroups ) printout(ALWAYS,name(),"+++ ParticleGroupConstructor: %s",p.c_str()); for ( const auto& [part_name,procs] : m_discreteProcesses ) { printout(ALWAYS,name(),"+++ DiscretePhysicsProcesses of particle %s",part_name.c_str()); for (ParticleProcesses::const_iterator ip = procs.begin(); ip != procs.end(); ++ip) { printout(ALWAYS,name(),"+++ Process %s", (*ip).name.c_str()); } } for ( const auto& [part_name, procs] : m_processes ) { printout(ALWAYS,name(),"+++ PhysicsProcesses of particle %s",part_name.c_str()); for ( const Process& p : procs ) { printout(ALWAYS,name(),"+++ Process %s ordAtRestDoIt=%d ordAlongSteptDoIt=%d ordPostStepDoIt=%d", p.name.c_str(),p.ordAtRestDoIt,p.ordAlongSteptDoIt,p.ordPostStepDoIt); } } } /// Add physics particle constructor by name void Geant4PhysicsList::addParticleConstructor(const std::string& part_name) { particles().emplace_back(part_name); } /// Add physics particle constructor by name void Geant4PhysicsList::addParticleGroup(const std::string& part_name) { particlegroups().emplace_back(part_name); } /// Add particle process by name with arguments void Geant4PhysicsList::addParticleProcess(const std::string& part_name, const std::string& proc_name, int ordAtRestDoIt, int ordAlongSteptDoIt, int ordPostStepDoIt) { Process p; p.name = proc_name; p.ordAtRestDoIt = ordAtRestDoIt; p.ordAlongSteptDoIt = ordAlongSteptDoIt; p.ordPostStepDoIt = ordPostStepDoIt; processes(part_name).emplace_back(p); } /// Add discrete particle process by name with arguments void Geant4PhysicsList::addDiscreteParticleProcess(const std::string& part_name, const std::string& proc_name) { Process p; p.name = proc_name; discreteProcesses(part_name).emplace_back(p); } /// Add PhysicsConstructor by name void Geant4PhysicsList::addPhysicsConstructor(const std::string& phys_name) { physics().emplace_back(phys_name); } /// Access processes for one particle type Geant4PhysicsList::ParticleProcesses& Geant4PhysicsList::processes(const string& nam) { if (auto i = m_processes.find(nam); i != m_processes.end()) return (*i).second; auto ret = m_processes.emplace(nam, ParticleProcesses()); return (*(ret.first)).second; } /// Access processes for one particle type (CONST) const Geant4PhysicsList::ParticleProcesses& Geant4PhysicsList::processes(const string& nam) const { if (auto i = m_processes.find(nam); i != m_processes.end()) return (*i).second; except("Failed to access the physics process '%s' [Unknown-Process]", nam.c_str()); throw runtime_error("Failed to access the physics process"); // never called anyway } /// Access discrete processes for one particle type Geant4PhysicsList::ParticleProcesses& Geant4PhysicsList::discreteProcesses(const string& nam) { if (auto i = m_discreteProcesses.find(nam); i != m_discreteProcesses.end()) return (*i).second; pair<PhysicsProcesses::iterator, bool> ret = m_discreteProcesses.emplace(nam, ParticleProcesses()); return (*(ret.first)).second; } /// Access discrete processes for one particle type (CONST) const Geant4PhysicsList::ParticleProcesses& Geant4PhysicsList::discreteProcesses(const string& nam) const { if (auto i = m_discreteProcesses.find(nam); i != m_discreteProcesses.end()) return (*i).second; except("Failed to access the physics process '%s' [Unknown-Process]", nam.c_str()); throw runtime_error("Failed to access the physics process"); // never called anyway } /// Access physics constructor by name (CONST) Geant4PhysicsList::PhysicsConstructor Geant4PhysicsList::physics(const std::string& nam) const { for ( const auto& ctor : m_physics ) { if ( ctor == nam ) { if ( nullptr == ctor.pointer ) except("Failed to instaniate the physics for constructor '%s'", nam.c_str()); return ctor; } } except("Failed to access the physics for constructor '%s' [Unknown physics]", nam.c_str()); throw runtime_error("Failed to access the physics process"); // never called anyway } /// Add PhysicsConstructor by name void Geant4PhysicsList::adoptPhysicsConstructor(Geant4Action* action) { if ( 0 != action ) { if ( G4VPhysicsConstructor* p = dynamic_cast<G4VPhysicsConstructor*>(action) ) { PhysicsConstructor ctor(action->name()); ctor.pointer = p; action->addRef(); m_physics.emplace_back(ctor); return; } except("Failed to adopt action object %s as physics constructor. [Invalid-Base]",action->c_name()); } except("Failed to adopt invalid Geant4Action as PhysicsConstructor. [Invalid-object]"); } /// Callback to construct particle decays void Geant4PhysicsList::constructPhysics(G4VModularPhysicsList* physics_pointer) { debug("constructPhysics %p", physics_pointer); for ( auto& ctor : m_physics ) { if ( 0 == ctor.pointer ) { if ( G4VPhysicsConstructor* p = PluginService::Create<G4VPhysicsConstructor*>(ctor) ) ctor.pointer = p; else except("Failed to create the physics for G4VPhysicsConstructor '%s'", ctor.c_str()); } physics_pointer->RegisterPhysics(ctor.pointer); info("Registered Geant4 physics constructor %s to physics list", ctor.c_str()); } } /// constructParticle callback void Geant4PhysicsList::constructParticles(G4VUserPhysicsList* physics_pointer) { debug("constructParticles %p", physics_pointer); /// Now define all particles for ( const auto& ctor : m_particles ) { G4ParticleDefinition* def = PluginService::Create<G4ParticleDefinition*>(ctor); if ( !def ) { /// Check if we have here a particle group constructor long* result = (long*) PluginService::Create<long>(ctor); if ( !result || *result != 1L ) { except("Failed to create particle type '%s' result=%d", ctor.c_str(), result); } info("Constructed Geant4 particle %s [using signature long (*)()]",ctor.c_str()); } } /// Now define all particles for ( const auto& ctor : m_particlegroups ) { /// Check if we have here a particle group constructor long* result = (long*) PluginService::Create<long>(ctor); if ( !result || *result != 1L ) { except("Failed to create particle type '%s' result=%d", ctor.c_str(), result); } info("Constructed Geant4 particle group %s [using signature long (*)()]",ctor.c_str()); } } /// Callback to construct processes (uses the G4 particle table) void Geant4PhysicsList::constructProcesses(G4VUserPhysicsList* physics_pointer) { debug("constructProcesses %p", physics_pointer); for ( const auto& [part_name, procs] : m_discreteProcesses ) { vector<G4ParticleDefinition*> defs(Geant4ParticleHandle::g4DefinitionsRegEx(part_name)); if ( defs.empty() ) { except("Particle:%s Cannot find the corresponding entry in the particle table.", part_name.c_str()); } for ( const Process& p : procs ) { if ( G4VProcess* g4 = PluginService::Create<G4VProcess*>(p.name) ) { for ( G4ParticleDefinition* particle : defs ) { G4ProcessManager* mgr = particle->GetProcessManager(); mgr->AddDiscreteProcess(g4); info("Particle:%s -> [%s] added discrete process %s", part_name.c_str(), particle->GetParticleName().c_str(), p.name.c_str()); } continue; } except("Cannot create discrete physics process %s", p.name.c_str()); } } for ( const auto& [part_name, procs] : m_processes ) { vector<G4ParticleDefinition*> defs(Geant4ParticleHandle::g4DefinitionsRegEx(part_name)); if (defs.empty()) { except("Particle:%s Cannot find the corresponding entry in the particle table.", part_name.c_str()); } for ( const Process& p : procs ) { if ( G4VProcess* g4 = PluginService::Create<G4VProcess*>(p.name) ) { for ( G4ParticleDefinition* particle : defs ) { G4ProcessManager* mgr = particle->GetProcessManager(); mgr->AddProcess(g4, p.ordAtRestDoIt, p.ordAlongSteptDoIt, p.ordPostStepDoIt); info("Particle:%s -> [%s] added process %s with flags (%d,%d,%d)", part_name.c_str(), particle->GetParticleName().c_str(), p.name.c_str(), p.ordAtRestDoIt, p.ordAlongSteptDoIt, p.ordPostStepDoIt); } continue; } except("Cannot create physics process %s", p.name.c_str()); } } } /// Enable physics list: actions necessary to be propagated to Geant4. void Geant4PhysicsList::enable(G4VUserPhysicsList* /* physics */) { } /// Standard constructor Geant4PhysicsListActionSequence::Geant4PhysicsListActionSequence(Geant4Context* ctxt, const string& nam) : Geant4Action(ctxt, nam), m_transportation(false), m_decays(false), m_rangecut(0.7*CLHEP::mm) { declareProperty("transportation", m_transportation); declareProperty("extends", m_extends); declareProperty("decays", m_decays); declareProperty("rangecut", m_rangecut); m_needsControl = true; InstanceCount::increment(this); } /// Default destructor Geant4PhysicsListActionSequence::~Geant4PhysicsListActionSequence() { m_actors(&Geant4Action::release); m_actors.clear(); m_process.clear(); InstanceCount::decrement(this); } #include <G4FastSimulationPhysics.hh> /// Extend physics list from factory: G4VUserPhysicsList* Geant4PhysicsListActionSequence::extensionList() { G4VModularPhysicsList* physics = ( m_extends.empty() ) ? new EmptyPhysics() : G4PhysListFactory().GetReferencePhysList(m_extends); #if 0 G4FastSimulationPhysics* fastSimulationPhysics = new G4FastSimulationPhysics(); // -- We now configure the fastSimulationPhysics object. // -- The gflash model (GFlashShowerModel, see ExGflashDetectorConstruction.cc) // -- is applicable to e+ and e- : we augment the physics list for these // -- particles (by adding a G4FastSimulationManagerProcess with below's // -- calls), this will make the fast simulation to be activated: fastSimulationPhysics->ActivateFastSimulation("e-"); fastSimulationPhysics->ActivateFastSimulation("e+"); // -- Register this fastSimulationPhysics to the physicsList, // -- when the physics list will be called by the run manager // -- (will happen at initialization of the run manager) // -- for physics process construction, the fast simulation // -- configuration will be applied as well. physics->RegisterPhysics( fastSimulationPhysics ); #endif // Register all physics constructors with the physics list constructPhysics(physics); // Ensure the particles and processes declared are also invoked. // Hence: Create a special physics constructor doing so. // Ownership is transferred to the physics list. // Do not delete this pointer afterwards.... physics->RegisterPhysics(new ParticlePhysics(this,physics)); return physics; } /// Install command control messenger if wanted void Geant4PhysicsListActionSequence::installCommandMessenger() { control()->addCall("dump", "Dump content of " + name(), Callback(this).make(&Geant4PhysicsListActionSequence::dump)); } /// Dump content to stdout void Geant4PhysicsListActionSequence::dump() { printout(ALWAYS,name(),"+++ Dump of physics list component(s)"); printout(ALWAYS,name(),"+++ Extension name %s",m_extends.c_str()); printout(ALWAYS,name(),"+++ Transportation flag: %d",m_transportation); printout(ALWAYS,name(),"+++ Program decays: %d",m_decays); printout(ALWAYS,name(),"+++ RangeCut: %f",m_rangecut); m_actors(&Geant4PhysicsList::dump); } /// Add an actor responding to all callbacks. Sequence takes ownership. void Geant4PhysicsListActionSequence::adopt(Geant4PhysicsList* action) { if (action) { action->addRef(); m_actors.add(action); return; } except("Geant4EventActionSequence: Attempt to add invalid actor!"); } /// Callback to construct particles void Geant4PhysicsListActionSequence::constructParticles(G4VUserPhysicsList* physics_pointer) { m_particle(physics_pointer); m_actors(&Geant4PhysicsList::constructParticles, physics_pointer); } /// Callback to execute physics constructors void Geant4PhysicsListActionSequence::constructPhysics(G4VModularPhysicsList* physics_pointer) { m_physics(physics_pointer); m_actors(&Geant4PhysicsList::constructPhysics, physics_pointer); } /// constructProcess callback void Geant4PhysicsListActionSequence::constructProcesses(G4VUserPhysicsList* physics_pointer) { m_actors(&Geant4PhysicsList::constructProcesses, physics_pointer); m_process(physics_pointer); if (m_decays) { constructDecays(physics_pointer); } } /// Callback to construct particle decays void Geant4PhysicsListActionSequence::constructDecays(G4VUserPhysicsList* physics_pointer) { G4ParticleTable* pt = G4ParticleTable::GetParticleTable(); G4ParticleTable::G4PTblDicIterator* iter = pt->GetIterator(); // Add Decay Process G4Decay* decay = new G4Decay(); info("ConstructDecays %p",physics_pointer); iter->reset(); while ((*iter)()) { G4ParticleDefinition* p = iter->value(); G4ProcessManager* mgr = p->GetProcessManager(); if (decay->IsApplicable(*p)) { mgr->AddProcess(decay); // set ordering for PostStepDoIt and AtRestDoIt mgr->SetProcessOrdering(decay, idxPostStep); mgr->SetProcessOrdering(decay, idxAtRest); } } } /// Enable physics list: actions necessary to be propagated to Geant4. void Geant4PhysicsListActionSequence::enable(G4VUserPhysicsList* physics_pointer) { m_actors(&Geant4PhysicsList::enable, physics_pointer); }