Merge pull request #187 from nkrah/activate_additional_physics_constr…

…uctors

core/opengate_core/g4_bindings/pyG4VModularPhysicsList.cpp

@@ -74,6 +74,7 @@ void init_G4VModularPhysicsList(py::module &m) {
                                     py::const_),
            py::return_value_policy::reference)
       .def("RegisterPhysics", &G4VModularPhysicsList::RegisterPhysics)
+      .def("ReplacePhysics", &G4VModularPhysicsList::ReplacePhysics)
       .def("RemovePhysics", py::overload_cast<G4VPhysicsConstructor *>(
                                 &G4VModularPhysicsList::RemovePhysics));
 }

core/opengate_core/g4_bindings/pyPhysicsLists.cpp

@@ -55,26 +55,24 @@ namespace py = pybind11;
 #include "G4VUserPhysicsList.hh"
 
 // macro for adding physics lists: no parameter
-#define ADD_PHYSICS_LIST0(m, plname)                                           \
-  py::class_<plname, G4VModularPhysicsList>(m, #plname).def(py::init<>());     \
-  AddPhysicsList(#plname);
+// #define ADD_PHYSICS_LIST0(m, plname) \
+//   py::class_<plname, G4VModularPhysicsList>(m, #plname).def(py::init<>()); \
+//   AddPhysicsList(#plname);
 
 // macro for adding physics lists: one int parameter
-#define ADD_PHYSICS_LIST1(m, plname)                                           \
-  py::class_<plname, G4VUserPhysicsList>(m, #plname).def(py::init<G4int>());   \
-  AddPhysicsList(#plname);
+// #define ADD_PHYSICS_LIST1(m, plname) \
+//   py::class_<plname, G4VUserPhysicsList>(m, #plname).def(py::init<G4int>());
+//   \ AddPhysicsList(#plname);
 
 // macro for adding physics lists: int+str parameter
-#define ADD_PHYSICS_LIST2(m, plname)                                           \
-  py::class_<plname, G4VUserPhysicsList>(m, #plname)                           \
-      .def(py::init<G4int, G4String>());                                       \
-  AddPhysicsList(#plname);
-
-// macro for adding physics constructor: one int parameter
-// (nodelete is needed because it is deleted in cpp side (runmanager?)
-// NK: Yes, the RunManager destructor calls the destructors of all
-// G4VPhysicsConstructor objects in a physics list
-// then also on py side, so seg fault during garbage collection)
+// #define ADD_PHYSICS_LIST2(m, plname) \
+//   py::class_<plname, G4VUserPhysicsList>(m, #plname) \
+//       .def(py::init<G4int, G4String>()); \
+//   AddPhysicsList(#plname);
+
+// macro for adding physics constructor: one int parameter (verbosity),
+// nodelete is needed because the G4 run manager deletes the physics list
+// on the cpp side, python should not delete the object
 #define ADD_PHYSICS_CONSTRUCTOR(plname)                                        \
   py::class_<plname, G4VPhysicsConstructor,                                    \
              std::unique_ptr<plname, py::nodelete>>(m, #plname)                \

opengate/physics/PhysicsEngine.py

@@ -81,7 +81,6 @@ def initialize_before_runmanager(self):
 
         """
         self.initialize_physics_list()
-        self.initialize_decay()
         self.initialize_em_options()
         self.initialize_user_limits_physics()
         self.initialize_parallel_world_physics()
@@ -94,7 +93,6 @@ def initialize_after_runmanager(self):
         # the global cuts with the physics list defaults.
         self.initialize_global_cuts()
         self.initialize_regions()
-
         self.initialize_g4_em_parameters()
 
     def initialize_parallel_world_physics(self):
@@ -116,28 +114,6 @@ def initialize_physics_list(self):
             )
         )
 
-    def initialize_decay(self):
-        """
-        G4DecayPhysics - defines all particles and their decay processes
-        G4RadioactiveDecayPhysics - defines radioactiveDecay for GenericIon
-        """
-        ui = self.physics_manager.user_info
-        if not ui.enable_decay:
-            return
-        # check if decay/radDecay already exist in the physics list
-        # (keep p and pp in self to prevent destruction)
-        # NOTE: can we use Replace() method from G4VModularPhysicsList?
-        self.g4_decay = self.g4_physics_list.GetPhysics("Decay")
-        if not self.g4_decay:
-            self.g4_decay = g4.G4DecayPhysics(1)
-            self.g4_physics_list.RegisterPhysics(self.g4_decay)
-        self.g4_radioactive_decay = self.g4_physics_list.GetPhysics(
-            "G4RadioactiveDecay"
-        )
-        if not self.g4_radioactive_decay:
-            self.g4_radioactive_decay = g4.G4RadioactiveDecayPhysics(1)
-            self.g4_physics_list.RegisterPhysics(self.g4_radioactive_decay)
-
     def initialize_em_options(self):
         # later
         pass

opengate/physics/PhysicsListManager.py

@@ -1,6 +1,7 @@
 import sys
 
 from opengate_core import G4PhysListFactory, G4VModularPhysicsList
+import opengate_core as g4
 
 from ..Decorators import requires_fatal
 from ..helpers import fatal
@@ -24,6 +25,14 @@ class PhysicsListManager(GateObjectSingleton):
         "G4OpticalPhysics",
     ]
 
+    special_physics_constructor_classes = {}
+    special_physics_constructor_classes["G4DecayPhysics"] = g4.G4DecayPhysics
+    special_physics_constructor_classes[
+        "G4RadioactiveDecayPhysics"
+    ] = g4.G4RadioactiveDecayPhysics
+    special_physics_constructor_classes["G4OpticalPhysics"] = g4.G4OpticalPhysics
+    special_physics_constructor_classes["G4EmDNAPhysics"] = g4.G4EmDNAPhysics
+
     def __init__(self, physics_manager, *args, **kwargs):
         super().__init__(*args, **kwargs)
         self.physics_manager = physics_manager
@@ -38,11 +47,13 @@ def create_physics_list_classes(self):
 
     def get_physics_list(self, physics_list_name):
         if physics_list_name in self.created_physics_list_classes:
-            return self.created_physics_list_classes[physics_list_name]()
+            physics_list = self.created_physics_list_classes[physics_list_name](
+                self.physics_manager.simulation.user_info.g4_verbose_level
+            )
         else:
             g4_factory = G4PhysListFactory()
             if g4_factory.IsReferencePhysList(physics_list_name):
-                return g4_factory.GetReferencePhysList(physics_list_name)
+                physics_list = g4_factory.GetReferencePhysList(physics_list_name)
             else:
                 s = (
                     f"Cannot find the physic list: {physics_list_name}\n"
@@ -51,6 +62,24 @@ def get_physics_list(self, physics_list_name):
                     f"Help : https://geant4-userdoc.web.cern.ch/UsersGuides/PhysicsListGuide/html/physicslistguide.html"
                 )
                 fatal(s)
+        # add special physics constructors
+        for (
+            spc,
+            switch,
+        ) in self.physics_manager.user_info.special_physics_constructors.items():
+            print(f"{spc} : {switch}")
+            if switch is True:
+                try:
+                    physics_list.ReplacePhysics(
+                        self.special_physics_constructor_classes[spc](
+                            self.physics_manager.simulation.user_info.g4_verbose_level
+                        )
+                    )
+                except KeyError:
+                    fatal(
+                        f"Special physics constructor named '{spc}' not found. Available constructors are: {self.special_physics_constructor_classes.keys()}."
+                    )
+        return physics_list
 
     def dump_info_physics_lists(self):
         g4_factory = G4PhysListFactory()
@@ -103,14 +132,12 @@ def create_modular_physics_list_class(g4_physics_constructor_class_name):
     return cls
 
 
-def init_method(self):
+def init_method(self, verbosity):
     """
     Init method of the dynamically created physics list class.
     - call the init method of the super class (G4VModularPhysicsList)
     - Create and register the physics constructor (G4VPhysicsConstructor)
     """
     G4VModularPhysicsList.__init__(self)
-    self.g4_physics_constructor = self.g4_physics_constructor_class(
-        1
-    )  # argument 1 means verbose=1
+    self.g4_physics_constructor = self.g4_physics_constructor_class(verbosity)
     self.RegisterPhysics(self.g4_physics_constructor)

opengate/physics/PhysicsManager.py

@@ -65,7 +65,6 @@ def _default_parameters(self):
         # keep the name to be able to come back to default
         self.default_physic_list = "QGSP_BERT_EMV"
         ui.physics_list_name = self.default_physic_list
-        ui.enable_decay = False
         """
         FIXME Energy range not clear : does not work in mono-thread mode
         Ignored for the moment (keep them to None)

opengate/physics/PhysicsUserInfo.py

@@ -1,7 +1,8 @@
-import opengate as gate
 import opengate_core as g4
 from box import Box
 from .PhysicsManager import PhysicsManager
+from .PhysicsListManager import PhysicsListManager
+from ..helpers import fatal
 
 
 class PhysicsUserInfo:
@@ -14,9 +15,14 @@ def __init__(self, simulation):
         # keep pointer to ref
         self.simulation = simulation
 
-        # physics list and decay
+        # physics list and optional physics constructors
         self.physics_list_name = None
-        self.enable_decay = False
+        # dictionary with names of additional physics constructors
+        # to be added to physics list (False = off by default)
+        # content is maintained by the PhysicsListManager
+        self.special_physics_constructors = Box()
+        for spc in PhysicsListManager.special_physics_constructor_classes:
+            self.special_physics_constructors[spc] = False
 
         # options related to the cuts and user limits
         # self.production_cuts = Box()
@@ -59,3 +65,28 @@ def __str__(self):
             f"user limits particles : {self.user_limits_particles}"
         )
         return s
+
+    # properties to quickly enable decay
+    # makes tests backwards compatible
+    # To be discussed whether an enable_decay switch makes sense
+    # Issue: if a physics list has already G4Decay in it, Gate will not
+    # deactivate it even if enable_decay = False
+    # Either enhance the logic, or leave it to the user to understand
+    # what is in the physics list they use
+    @property
+    def enable_decay(self):
+        switch1 = self.special_physics_constructors["G4DecayPhysics"]
+        switch2 = self.special_physics_constructors["G4RadioactiveDecayPhysics"]
+        if switch1 is True and switch2 is True:
+            return True
+        elif switch1 is False and switch2 is False:
+            return False
+        else:
+            fatal(
+                f"Inconsistent G4Decay constructors: G4DecayPhysics = {switch1}, G4RadioactiveDecayPhysics = {switch2}."
+            )
+
+    @enable_decay.setter
+    def enable_decay(self, value):
+        self.special_physics_constructors["G4DecayPhysics"] = value
+        self.special_physics_constructors["G4RadioactiveDecayPhysics"] = value

opengate/tests/src/test013_phys_lists_1.py

@@ -18,7 +18,11 @@
 # change physics
 p = sim.get_physics_user_info()
 p.physics_list_name = "G4EmStandardPhysics_option4"
-p.enable_decay = True
+# enable decay via switch:
+# p.enable_decay = True
+# or by activating the physics constructors:
+p.special_physics_constructors.G4DecayPhysics = True
+p.special_physics_constructors.G4RadioactiveDecayPhysics = True
 
 um = gate.g4_units("um")
 global_cut = 7 * um