DDG4/src/Geant4PhysicsList.cpp

0001 //==========================================================================
0002 //  AIDA Detector description implementation
0003 //--------------------------------------------------------------------------
0004 // Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
0005 // All rights reserved.
0006 //
0007 // For the licensing terms see $DD4hepINSTALL/LICENSE.
0008 // For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
0009 //
0010 // Author     : M.Frank
0011 //
0012 //==========================================================================
0013
0014 // Framework include files
0015 #include <DDG4/Geant4PhysicsList.h>
0016 #include <DDG4/Geant4UIMessenger.h>
0017 #include <DDG4/Geant4Particle.h>
0018 #include <DDG4/Geant4Kernel.h>
0019 #include <DD4hep/InstanceCount.h>
0020 #include <DD4hep/Printout.h>
0021 #include <DD4hep/Plugins.h>
0022
0023 // Geant4 include files
0024 #include <G4VPhysicsConstructor.hh>
0025 #include <G4PhysListFactory.hh>
0026 #include <G4ProcessManager.hh>
0027 #include <G4ParticleTable.hh>
0028 #include <G4RunManager.hh>
0029 #include <G4VProcess.hh>
0030 #include <G4Decay.hh>
0031 #include <G4EmParameters.hh>
0032 #include <G4HadronicParameters.hh>
0033
0034 // C/C++ include files
0035 #include <stdexcept>
0036 #include <regex.h>
0037
0038 using namespace dd4hep::sim;
0039
0040 namespace {
0041
0042   struct MyPhysics : G4VUserPhysicsList  {
0043     void AddTransportation()   {  this->G4VUserPhysicsList::AddTransportation(); }
0044   };
0045
0046   struct EmptyPhysics : public G4VModularPhysicsList {
0047     EmptyPhysics() = default;
0048     virtual ~EmptyPhysics() = default;
0049   };
0050   struct ParticlePhysics : public G4VPhysicsConstructor {
0051     Geant4PhysicsListActionSequence* seq;
0052     G4VUserPhysicsList*              phys;
0053     ParticlePhysics(Geant4PhysicsListActionSequence* s, G4VUserPhysicsList* p) : seq(s), phys(p)  { }
0054     virtual ~ParticlePhysics() = default;
0055     virtual void ConstructProcess()  {
0056       seq->constructProcesses(phys);
0057       if ( seq->transportation() )   {
0058         MyPhysics* ph = (MyPhysics*)phys;
0059         ph->AddTransportation();
0060       }
0061     }
0062     virtual void ConstructParticle()  {
0063       seq->constructParticles(phys);
0064     }
0065   };
0066 }
0067
0068 /// Default constructor
0069 Geant4PhysicsList::Process::Process()
0070   : ordAtRestDoIt(-1), ordAlongSteptDoIt(-1), ordPostStepDoIt(-1)  {
0071 }
0072 /// Copy constructor
0073 Geant4PhysicsList::Process::Process(const Process& p)
0074   : name(p.name), ordAtRestDoIt(p.ordAtRestDoIt), ordAlongSteptDoIt(p.ordAlongSteptDoIt), ordPostStepDoIt(p.ordPostStepDoIt)  {
0075 }
0076
0077 /// Assignment operator
0078 Geant4PhysicsList::Process& Geant4PhysicsList::Process::operator=(const Process& p)  {
0079   if ( this != &p )  {
0080     name = p.name;
0081     ordAtRestDoIt     = p.ordAtRestDoIt;
0082     ordAlongSteptDoIt = p.ordAlongSteptDoIt;
0083     ordPostStepDoIt   = p.ordPostStepDoIt;
0084   }
0085   return *this;
0086 }
0087
0088 /// Standard constructor
0089 Geant4PhysicsList::Geant4PhysicsList(Geant4Context* ctxt, const std::string& nam)
0090   : Geant4Action(ctxt, nam)  {
0091   InstanceCount::increment(this);
0092 }
0093
0094 /// Default destructor
0095 Geant4PhysicsList::~Geant4PhysicsList()  {
0096   InstanceCount::decrement(this);
0097 }
0098
0099 /// Install command control messenger if wanted
0100 void Geant4PhysicsList::installCommandMessenger()   {
0101   control()->addCall("dump", "Dump content of " + name(), Callback(this).make(&Geant4PhysicsList::dump));
0102 }
0103
0104 /// Dump content to stdout
0105 void Geant4PhysicsList::dump()    {
0106   if ( !m_physics.empty() && !m_particles.empty() && !m_processes.empty() )  {
0107     printout(ALWAYS,name(),"+++ Geant4PhysicsList Dump");
0108   }
0109   for ( const auto& p : m_physics)
0110     printout(ALWAYS,name(),"+++ PhysicsConstructor:           %s",p.c_str());
0111   for ( const auto& p : m_particles )
0112     printout(ALWAYS,name(),"+++ ParticleConstructor:          %s",p.c_str());
0113   for ( const auto& p : m_particlegroups )
0114     printout(ALWAYS,name(),"+++ ParticleGroupConstructor:     %s",p.c_str());
0115   for ( const auto& [part_name,procs] : m_discreteProcesses )   {
0116     printout(ALWAYS,name(),"+++ DiscretePhysicsProcesses of particle  %s",part_name.c_str());
0117     for (ParticleProcesses::const_iterator ip = procs.begin(); ip != procs.end(); ++ip)  {
0118       printout(ALWAYS,name(),"+++        Process    %s", (*ip).name.c_str());
0119     }
0120   }
0121   for ( const auto& [part_name, procs] : m_processes )  {
0122     printout(ALWAYS,name(),"+++ PhysicsProcesses of particle  %s",part_name.c_str());
0123     for ( const Process& p : procs )    {
0124       printout(ALWAYS,name(),"+++        Process    %s  ordAtRestDoIt=%d ordAlongSteptDoIt=%d ordPostStepDoIt=%d",
0125                p.name.c_str(),p.ordAtRestDoIt,p.ordAlongSteptDoIt,p.ordPostStepDoIt);
0126     }
0127   }
0128 }
0129
0130 /// Add physics particle constructor by name
0131 void Geant4PhysicsList::addParticleConstructor(const std::string& part_name)   {
0132   particles().emplace_back(part_name);
0133 }
0134
0135 /// Add physics particle constructor by name
0136 void Geant4PhysicsList::addParticleGroup(const std::string& part_name)   {
0137   particlegroups().emplace_back(part_name);
0138 }
0139
0140 /// Add particle process by name with arguments
0141 void Geant4PhysicsList::addParticleProcess(const std::string& part_name,
0142                                            const std::string& proc_name,
0143                                            int ordAtRestDoIt,
0144                                            int ordAlongSteptDoIt,
0145                                            int ordPostStepDoIt)
0146 {
0147   Process p;
0148   p.name = proc_name;
0149   p.ordAtRestDoIt     = ordAtRestDoIt;
0150   p.ordAlongSteptDoIt = ordAlongSteptDoIt;
0151   p.ordPostStepDoIt   = ordPostStepDoIt;
0152   processes(part_name).emplace_back(p);
0153 }
0154
0155 /// Add discrete particle process by name with arguments
0156 void Geant4PhysicsList::addDiscreteParticleProcess(const std::string& part_name,
0157                                                    const std::string& proc_name)
0158 {
0159   Process p;
0160   p.name = proc_name;
0161   discreteProcesses(part_name).emplace_back(p);
0162 }
0163
0164 /// Add PhysicsConstructor by name
0165 void Geant4PhysicsList::addPhysicsConstructor(const std::string& phys_name)  {
0166   physics().emplace_back(phys_name);
0167 }
0168
0169 /// Access processes for one particle type
0170 Geant4PhysicsList::ParticleProcesses& Geant4PhysicsList::processes(const std::string& nam)  {
0171   if (auto i = m_processes.find(nam); i != m_processes.end())
0172     return (*i).second;
0173   auto ret = m_processes.emplace(nam, ParticleProcesses());
0174   return (*(ret.first)).second;
0175 }
0176
0177 /// Access processes for one particle type (CONST)
0178 const Geant4PhysicsList::ParticleProcesses& Geant4PhysicsList::processes(const std::string& nam) const {
0179   if (auto i = m_processes.find(nam); i != m_processes.end())
0180     return (*i).second;
0181   except("Failed to access the physics process '%s' [Unknown-Process]", nam.c_str());
0182   throw std::runtime_error("Failed to access the physics process"); // never called anyway
0183 }
0184
0185 /// Access discrete processes for one particle type
0186 Geant4PhysicsList::ParticleProcesses& Geant4PhysicsList::discreteProcesses(const std::string& nam)  {
0187   if (auto i = m_discreteProcesses.find(nam); i != m_discreteProcesses.end())
0188     return (*i).second;
0189   auto ret = m_discreteProcesses.emplace(nam, ParticleProcesses());
0190   return (*(ret.first)).second;
0191 }
0192
0193 /// Access discrete processes for one particle type (CONST)
0194 const Geant4PhysicsList::ParticleProcesses& Geant4PhysicsList::discreteProcesses(const std::string& nam) const {
0195   if (auto i = m_discreteProcesses.find(nam); i != m_discreteProcesses.end())
0196     return (*i).second;
0197   except("Failed to access the physics process '%s' [Unknown-Process]", nam.c_str());
0198   throw std::runtime_error("Failed to access the physics process"); // never called anyway
0199 }
0200
0201 /// Access physics constructor by name (CONST)
0202 Geant4PhysicsList::PhysicsConstructor Geant4PhysicsList::physics(const std::string& nam)  const    {
0203   for ( const auto& ctor : m_physics )   {
0204     if ( ctor == nam )  {
0205       if ( nullptr == ctor.pointer )
0206         except("Failed to instaniate the physics for constructor '%s'", nam.c_str());
0207       return ctor;
0208     }
0209   }
0210   except("Failed to access the physics for constructor '%s' [Unknown physics]", nam.c_str());
0211   throw std::runtime_error("Failed to access the physics process"); // never called anyway
0212 }
0213
0214 /// Add PhysicsConstructor by name
0215 void Geant4PhysicsList::adoptPhysicsConstructor(Geant4Action* action)  {
0216   if ( 0 != action )   {
0217     if ( G4VPhysicsConstructor* p = dynamic_cast<G4VPhysicsConstructor*>(action) )  {
0218       PhysicsConstructor ctor(action->name());
0219       ctor.physics_list = nullptr;
0220       ctor.pointer = p;
0221       action->addRef();
0222       m_physics.emplace_back(ctor);
0223       return;
0224     }
0225     except("Failed to adopt action object %s as physics constructor. [Invalid-Base]",action->c_name());
0226   }
0227   except("Failed to adopt invalid Geant4Action as PhysicsConstructor. [Invalid-object]");
0228 }
0229
0230 /// Callback to construct particle decays
0231 void Geant4PhysicsList::constructPhysics(G4VModularPhysicsList* physics_pointer)  {
0232   debug("constructPhysics %p", physics_pointer);
0233   for ( auto& ctor : m_physics )   {
0234     if ( 0 == ctor.pointer )   {
0235       if ( G4VPhysicsConstructor* p = PluginService::Create<G4VPhysicsConstructor*>(ctor) )  {
0236         ctor.pointer = p;
0237       }
0238       else  {
0239         except("Failed to create the physics for G4VPhysicsConstructor '%s'", ctor.c_str());
0240       }
0241     }
0242     if( !ctor.physics_list )  {
0243       ctor.physics_list = physics_pointer;
0244       physics_pointer->RegisterPhysics(ctor.pointer);
0245       info("+++ Registered Geant4 physics constructor %s to modular physics list id:%d",
0246            ctor.c_str(), physics_pointer->GetInstanceID());
0247     }
0248   }
0249 }
0250
0251 /// Callback to add a physics type to the physics list
0252 G4VPhysicsConstructor* Geant4PhysicsList::addPhysicsConstructorType(const std::string& physics_type)  {
0253   debug("addPhysics %s", physics_type.c_str());
0254   for ( auto& ctor : m_physics )   {
0255     if( physics_type == ctor )  {
0256       warning("+++ Physics type %s is already present. [using existing, creation denied]",
0257               physics_type.c_str());
0258       return ctor.pointer;
0259     }
0260   }
0261   PhysicsConstructor ctor(physics_type);
0262   ctor.physics_list = nullptr;
0263   ctor.pointer = PluginService::Create<G4VPhysicsConstructor*>(physics_type);
0264   if( ctor.pointer == nullptr )  {
0265     return nullptr;
0266   }
0267   m_physics.push_back(ctor);
0268   info("+++ Added Geant4 physics constructor %s to physics list", ctor.c_str());
0269   return ctor.pointer;
0270 }
0271
0272 /// constructParticle callback
0273 void Geant4PhysicsList::constructParticles(G4VUserPhysicsList* physics_pointer)   {
0274   debug("constructParticles %p", physics_pointer);
0275   /// Now define all particles
0276   for ( const auto& ctor : m_particles )   {
0277     G4ParticleDefinition* def = PluginService::Create<G4ParticleDefinition*>(ctor);
0278     if ( !def )  {
0279       /// Check if we have here a particle group constructor
0280       long* result = (long*) PluginService::Create<long>(ctor);
0281       if ( !result || *result != 1L )   {
0282         except("Failed to create particle type '%s' result=%d", ctor.c_str(), result);
0283       }
0284       info("Constructed Geant4 particle %s [using signature long (*)()]",ctor.c_str());
0285     }
0286   }
0287   /// Now define all particles
0288   for ( const auto& ctor : m_particlegroups )  {
0289     /// Check if we have here a particle group constructor
0290     long* result = (long*) PluginService::Create<long>(ctor);
0291     if ( !result || *result != 1L )  {
0292       except("Failed to create particle type '%s' result=%d", ctor.c_str(), result);
0293     }
0294     info("Constructed Geant4 particle group %s [using signature long (*)()]",ctor.c_str());
0295   }
0296 }
0297
0298 /// Callback to construct processes (uses the G4 particle table)
0299 void Geant4PhysicsList::constructProcesses(G4VUserPhysicsList* physics_pointer)   {
0300   debug("constructProcesses %p", physics_pointer);
0301   for ( const auto& [part_name, procs] : m_discreteProcesses )  {
0302     std::vector<G4ParticleDefinition*> defs(Geant4ParticleHandle::g4DefinitionsRegEx(part_name));
0303     if ( defs.empty() )  {
0304       except("Particle:%s Cannot find the corresponding entry in the particle table.", part_name.c_str());
0305     }
0306     for ( const Process& p : procs )  {
0307       if ( G4VProcess* g4 = PluginService::Create<G4VProcess*>(p.name) )   {
0308         for ( G4ParticleDefinition* particle : defs )   {
0309           G4ProcessManager* mgr = particle->GetProcessManager();
0310           mgr->AddDiscreteProcess(g4);
0311           info("Particle:%s -> [%s] added discrete process %s",
0312                part_name.c_str(), particle->GetParticleName().c_str(), p.name.c_str());
0313         }
0314         continue;
0315       }
0316       except("Cannot create discrete physics process %s", p.name.c_str());
0317     }
0318   }
0319   for ( const auto& [part_name, procs] : m_processes )   {
0320     std::vector<G4ParticleDefinition*> defs(Geant4ParticleHandle::g4DefinitionsRegEx(part_name));
0321     if (defs.empty())  {
0322       except("Particle:%s Cannot find the corresponding entry in the particle table.", part_name.c_str());
0323     }
0324     for ( const Process& p : procs )  {
0325       if ( G4VProcess* g4 = PluginService::Create<G4VProcess*>(p.name) )   {
0326         for ( G4ParticleDefinition* particle : defs )   {
0327           G4ProcessManager* mgr = particle->GetProcessManager();
0328           mgr->AddProcess(g4, p.ordAtRestDoIt, p.ordAlongSteptDoIt, p.ordPostStepDoIt);
0329           info("Particle:%s -> [%s] added process %s with flags (%d,%d,%d)",
0330                part_name.c_str(), particle->GetParticleName().c_str(), p.name.c_str(),
0331                p.ordAtRestDoIt, p.ordAlongSteptDoIt, p.ordPostStepDoIt);
0332         }
0333         continue;
0334       }
0335       except("Cannot create physics process %s", p.name.c_str());
0336     }
0337   }
0338 }
0339
0340 /// Enable physics list: actions necessary to be propagated to Geant4.
0341 void Geant4PhysicsList::enable(G4VUserPhysicsList* /* physics */)  {
0342 }
0343
0344 /// Standard constructor
0345 Geant4PhysicsListActionSequence::Geant4PhysicsListActionSequence(Geant4Context* ctxt, const std::string& nam)
0346   : Geant4Action(ctxt, nam), m_rangecut(0.7*CLHEP::mm)
0347 {
0348   declareProperty("transportation", m_transportation);
0349   declareProperty("extends",  m_extends);
0350   declareProperty("decays",   m_decays);
0351   declareProperty("rangecut", m_rangecut);
0352   declareProperty("verbosity", m_verbosity);
0353   declareProperty("modular_physics_verbosity", m_physics_verbosity);
0354   m_needsControl = true;
0355   InstanceCount::increment(this);
0356 }
0357
0358 /// Default destructor
0359 Geant4PhysicsListActionSequence::~Geant4PhysicsListActionSequence()  {
0360   m_actors(&Geant4Action::release);
0361   m_actors.clear();
0362   m_process.clear();
0363   InstanceCount::decrement(this);
0364 }
0365
0366 #include <G4FastSimulationPhysics.hh>
0367
0368 /// Extend physics list from factory:
0369 G4VUserPhysicsList* Geant4PhysicsListActionSequence::extensionList()    {
0370   G4VModularPhysicsList* physics = ( m_extends.empty() )
0371     ? new EmptyPhysics()
0372     : G4PhysListFactory().GetReferencePhysList(m_extends);
0373
0374 #if 0
0375   G4FastSimulationPhysics* fastSimulationPhysics = new G4FastSimulationPhysics();
0376   // -- We now configure the fastSimulationPhysics object.
0377   // -- The gflash model (GFlashShowerModel, see ExGflashDetectorConstruction.cc)
0378   // -- is applicable to e+ and e- : we augment the physics list for these
0379   // -- particles (by adding a G4FastSimulationManagerProcess with below's
0380   // -- calls), this will make the fast simulation to be activated:
0381   fastSimulationPhysics->ActivateFastSimulation("e-");
0382   fastSimulationPhysics->ActivateFastSimulation("e+");
0383   // -- Register this fastSimulationPhysics to the physicsList,
0384   // -- when the physics list will be called by the run manager
0385   // -- (will happen at initialization of the run manager)
0386   // -- for physics process construction, the fast simulation
0387   // -- configuration will be applied as well.
0388   physics->RegisterPhysics( fastSimulationPhysics );
0389 #endif
0390   // Register all physics constructors with the physics list
0391   constructPhysics(physics);
0392   // Ensure the particles and processes declared are also invoked.
0393   // Hence: Create a special physics constructor doing so.
0394   // Ownership is transferred to the physics list.
0395   // Do not delete this pointer afterwards....
0396   physics->RegisterPhysics(new ParticlePhysics(this,physics));
0397
0398   //Setting verbosity for pieces of the physics
0399   physics->SetVerboseLevel(m_verbosity);
0400   G4EmParameters::Instance()->SetVerbose(m_verbosity);
0401   G4HadronicParameters::Instance()->SetVerboseLevel(m_verbosity);
0402
0403   return physics;
0404 }
0405
0406 /// Install command control messenger if wanted
0407 void Geant4PhysicsListActionSequence::installCommandMessenger()   {
0408   control()->addCall("dump", "Dump content of " + name(), Callback(this).make(&Geant4PhysicsListActionSequence::dump));
0409 }
0410
0411 /// Dump content to stdout
0412 void Geant4PhysicsListActionSequence::dump()    {
0413   printout(ALWAYS,name(),"+++ Dump of physics list component(s)");
0414   printout(ALWAYS,name(),"+++ Extension name       %s",m_extends.c_str());
0415   printout(ALWAYS,name(),"+++ Transportation flag: %d",m_transportation);
0416   printout(ALWAYS,name(),"+++ Program decays:      %d",m_decays);
0417   printout(ALWAYS,name(),"+++ RangeCut:            %f",m_rangecut);
0418   printout(ALWAYS,name(),"+++ Verbosity:           %i",m_verbosity);
0419   m_actors(&Geant4PhysicsList::dump);
0420 }
0421
0422 /// Add an actor responding to all callbacks. Sequence takes ownership.
0423 void Geant4PhysicsListActionSequence::adopt(Geant4PhysicsList* action)  {
0424   if (action)  {
0425     action->addRef();
0426     m_actors.add(action);
0427     return;
0428   }
0429   except("Geant4EventActionSequence: Attempt to add invalid actor!");
0430 }
0431
0432 /// Callback to construct particles
0433 void Geant4PhysicsListActionSequence::constructParticles(G4VUserPhysicsList* physics_pointer)  {
0434   m_particle(physics_pointer);
0435   m_actors(&Geant4PhysicsList::constructParticles, physics_pointer);
0436 }
0437
0438 /// Callback to execute physics constructors
0439 void Geant4PhysicsListActionSequence::constructPhysics(G4VModularPhysicsList* physics_pointer)  {
0440   m_physics(physics_pointer);
0441   m_actors(&Geant4PhysicsList::constructPhysics, physics_pointer);
0442 }
0443
0444 /// constructProcess callback
0445 void Geant4PhysicsListActionSequence::constructProcesses(G4VUserPhysicsList* physics_pointer)  {
0446   m_actors(&Geant4PhysicsList::constructProcesses, physics_pointer);
0447   m_process(physics_pointer);
0448   if (m_decays)  {
0449     constructDecays(physics_pointer);
0450   }
0451 }
0452
0453 /// Callback to construct particle decays
0454 void Geant4PhysicsListActionSequence::constructDecays(G4VUserPhysicsList* physics_pointer)  {
0455   G4ParticleTable* pt = G4ParticleTable::GetParticleTable();
0456   G4ParticleTable::G4PTblDicIterator* iter = pt->GetIterator();
0457   // Add Decay Process
0458   G4Decay* decay = new G4Decay();
0459   info("ConstructDecays %p",physics_pointer);
0460   iter->reset();
0461   while ((*iter)())  {
0462     G4ParticleDefinition* p = iter->value();
0463     G4ProcessManager* mgr = p->GetProcessManager();
0464     if (decay->IsApplicable(*p))  {
0465       mgr->AddProcess(decay);
0466       // set ordering for PostStepDoIt and AtRestDoIt
0467       mgr->SetProcessOrdering(decay, idxPostStep);
0468       mgr->SetProcessOrdering(decay, idxAtRest);
0469     }
0470   }
0471 }
0472
0473 /// Enable physics list: actions necessary to be propagated to Geant4.
0474 void Geant4PhysicsListActionSequence::enable(G4VUserPhysicsList* physics_pointer)   {
0475   if( m_physics_verbosity >= 0 )  {
0476     physics_pointer->SetVerboseLevel(m_physics_verbosity);
0477   }
0478   m_actors(&Geant4PhysicsList::enable, physics_pointer);
0479 }
0480