Merge pull request #600 from OpenGATE/new_kill_interacting_particle

New kill interacting particle

core/opengate_core/opengate_core.cpp

@@ -309,6 +309,8 @@ void init_GateARFTrainingDatasetActor(py::module &m);
 
 void init_GateKillActor(py::module &);
 
+void init_GateKillAccordingProcessesActor(py::module &);
+
 void init_itk_image(py::module &);
 
 void init_GateImageNestedParameterisation(py::module &);
@@ -583,6 +585,7 @@ PYBIND11_MODULE(opengate_core, m) {
   init_GateARFActor(m);
   init_GateARFTrainingDatasetActor(m);
   init_GateKillActor(m);
+  init_GateKillAccordingProcessesActor(m);
   init_GateDigiAttributeManager(m);
   init_GateVDigiAttribute(m);
   init_GateExceptionHandler(m);

core/opengate_core/opengate_lib/GateKillAccordingProcessesActor.cpp

@@ -0,0 +1,128 @@
+/* --------------------------------------------------
+   Copyright (C): OpenGATE Collaboration
+   This software is distributed under the terms
+   of the GNU Lesser General  Public Licence (LGPL)
+   See LICENSE.md for further details
+   ------------------------------------ -------------- */
+
+#include "GateKillAccordingProcessesActor.h"
+#include "G4LogicalVolumeStore.hh"
+#include "G4ParticleTable.hh"
+#include "G4PhysicalVolumeStore.hh"
+#include "G4ProcessManager.hh"
+#include "G4VProcess.hh"
+#include "G4ios.hh"
+#include "GateHelpers.h"
+#include "GateHelpersDict.h"
+
+G4Mutex SetNbKillAccordingProcessesMutex = G4MUTEX_INITIALIZER;
+
+GateKillAccordingProcessesActor::GateKillAccordingProcessesActor(
+    py::dict &user_info)
+    : GateVActor(user_info, true) {}
+
+std::vector<G4String>
+GateKillAccordingProcessesActor::GetListOfPhysicsListProcesses() {
+  std::vector<G4String> listOfAllProcesses = {};
+
+  G4ParticleTable *particleTable = G4ParticleTable::GetParticleTable();
+
+  for (G4int i = 0; i < particleTable->size(); ++i) {
+    G4ParticleDefinition *particle = particleTable->GetParticle(i);
+    G4String particleName = particle->GetParticleName();
+    G4ProcessManager *processManager = particle->GetProcessManager();
+    if (!processManager)
+      continue;
+    G4int numProcesses = processManager->GetProcessListLength();
+    for (G4int j = 0; j < numProcesses; ++j) {
+      const G4VProcess *process = (*(processManager->GetProcessList()))[j];
+      G4String processName = process->GetProcessName();
+      if (std::find(listOfAllProcesses.begin(), listOfAllProcesses.end(),
+                    processName) == listOfAllProcesses.end())
+        listOfAllProcesses.push_back(processName);
+    }
+  }
+  return listOfAllProcesses;
+}
+
+void GateKillAccordingProcessesActor::InitializeUserInfo(py::dict &user_info) {
+  GateVActor::InitializeUserInfo(user_info);
+  fProcessesToKill = DictGetVecStr(user_info, "processes_to_kill");
+  fIsRayleighAnInteraction =
+      DictGetBool(user_info, "is_rayleigh_an_interaction");
+}
+
+void GateKillAccordingProcessesActor::BeginOfRunAction(const G4Run *run) {
+  fNbOfKilledParticles = 0;
+  std::vector<G4String> listOfAllProcesses = GetListOfPhysicsListProcesses();
+  listOfAllProcesses.push_back("all");
+  for (auto process : fProcessesToKill) {
+    if (std::find(listOfAllProcesses.begin(), listOfAllProcesses.end(),
+                  process) == listOfAllProcesses.end()) {
+      G4String errorMessage =
+          "Process '" + process + "' not found. Existing processes are '";
+      for (auto aProcess : listOfAllProcesses) {
+        errorMessage = errorMessage + aProcess + "', ";
+      }
+      errorMessage.pop_back();
+      errorMessage.pop_back();
+      G4Exception("CheckProcessExistence", // Exception origin
+                  "ProcessNotFound.1",     // Exception code
+                  FatalException,          // Exception severity
+                  errorMessage);
+    }
+  }
+  if (fProcessesToKill[0] == "all") {
+    if (fProcessesToKill.size() == 1) {
+      fKillIfAnyInteraction = true;
+    }
+  }
+}
+
+void GateKillAccordingProcessesActor::PreUserTrackingAction(
+    const G4Track *track) {
+  fIsFirstStep = true;
+}
+
+void GateKillAccordingProcessesActor::SteppingAction(G4Step *step) {
+
+  G4String logNameMotherVolume = G4LogicalVolumeStore::GetInstance()
+                                     ->GetVolume(fAttachedToVolumeName)
+                                     ->GetName();
+  G4String physicalVolumeNamePreStep = "None";
+
+  G4String processName = "None";
+  const G4VProcess *process = step->GetPostStepPoint()->GetProcessDefinedStep();
+  if (process != 0)
+    processName = process->GetProcessName();
+
+  // Positron exception to retrieve the annihilation process, since it's an at
+  // rest process most of the time
+
+  if ((step->GetTrack()->GetParticleDefinition()->GetParticleName() == "e+") &&
+      (step->GetTrack()->GetTrackStatus() == 1))
+    processName = "annihil";
+
+  if (fKillIfAnyInteraction) {
+    if (processName != "Transportation") {
+      if (fIsRayleighAnInteraction == true) {
+        step->GetTrack()->SetTrackStatus(fKillTrackAndSecondaries);
+        G4AutoLock mutex(&SetNbKillAccordingProcessesMutex);
+        fNbOfKilledParticles++;
+      } else {
+        if (processName != "Rayl") {
+          step->GetTrack()->SetTrackStatus(fKillTrackAndSecondaries);
+          G4AutoLock mutex(&SetNbKillAccordingProcessesMutex);
+          fNbOfKilledParticles++;
+        }
+      }
+    }
+  } else {
+    if (std::find(fProcessesToKill.begin(), fProcessesToKill.end(),
+                  processName) != fProcessesToKill.end()) {
+      step->GetTrack()->SetTrackStatus(fKillTrackAndSecondaries);
+      G4AutoLock mutex(&SetNbKillAccordingProcessesMutex);
+      fNbOfKilledParticles++;
+    }
+  }
+}

core/opengate_core/opengate_lib/GateKillAccordingProcessesActor.h

@@ -0,0 +1,43 @@
+/* --------------------------------------------------
+   Copyright (C): OpenGATE Collaboration
+   This software is distributed under the terms
+   of the GNU Lesser General  Public Licence (LGPL)
+   See LICENSE.md for further details
+   -------------------------------------------------- */
+
+#ifndef GateKillAccordingProcessesActor_h
+#define GateKillAccordingProcessesActor_h
+
+#include "G4Cache.hh"
+#include "GateVActor.h"
+#include <pybind11/stl.h>
+
+namespace py = pybind11;
+
+class GateKillAccordingProcessesActor : public GateVActor {
+
+public:
+  // Constructor
+  GateKillAccordingProcessesActor(py::dict &user_info);
+  std::vector<G4String> GetListOfPhysicsListProcesses();
+
+  void InitializeUserInfo(py::dict &user_info) override;
+
+  void BeginOfRunAction(const G4Run *) override;
+
+  // Main function called every step in attached volume
+  void SteppingAction(G4Step *) override;
+
+  void PreUserTrackingAction(const G4Track *) override;
+
+  std::vector<G4String> fParticlesTypeToKill;
+  G4bool fIsFirstStep = true;
+  std::vector<std::string> fProcessesToKill;
+  std::vector<std::string> fListOfVolumeAncestor;
+  G4bool fKillIfAnyInteraction = false;
+  G4bool fIsRayleighAnInteraction = false;
+
+  long fNbOfKilledParticles{};
+};
+
+#endif

core/opengate_core/opengate_lib/GateKillActor.cpp

@@ -17,7 +17,10 @@ GateKillActor::GateKillActor(py::dict &user_info)
   fNbOfKilledParticles = 0;
 }
 
-void GateKillActor::StartSimulationAction() { fNbOfKilledParticles = 0; }
+void GateKillActor::StartSimulationAction() {
+  fNbOfKilledParticles = 0;
+  std::cout << "lol" << std::endl;
+}
 
 void GateKillActor::SteppingAction(G4Step *step) {
   auto track = step->GetTrack();

core/opengate_core/opengate_lib/pyGateKillAccordingProcessesActor.cpp

@@ -0,0 +1,21 @@
+/* --------------------------------------------------
+   Copyright (C): OpenGATE Collaboration
+   This software is distributed under the terms
+   of the GNU Lesser General  Public Licence (LGPL)
+   See LICENSE.md for further details
+   -------------------------------------------------- */
+
+#include <pybind11/pybind11.h>
+
+namespace py = pybind11;
+
+#include "GateKillAccordingProcessesActor.h"
+
+void init_GateKillAccordingProcessesActor(py::module &m) {
+  py::class_<GateKillAccordingProcessesActor,
+             std::unique_ptr<GateKillAccordingProcessesActor, py::nodelete>,
+             GateVActor>(m, "GateKillAccordingProcessesActor")
+      .def_readwrite("fListOfVolumeAncestor",
+                     &GateKillAccordingProcessesActor::fListOfVolumeAncestor)
+      .def(py::init<py::dict &>());
+}

opengate/actors/miscactors.py

@@ -5,8 +5,9 @@
 from ..utility import g4_units, g4_best_unit_tuple
 from .actoroutput import ActorOutputBase
 from ..serialization import dump_json
-from ..exception import warning
+from ..exception import fatal, warning
 from ..base import process_cls
+from anytree import Node, RenderTree
 
 
 def _setter_hook_stats_actor_output_filename(self, output_filename):
@@ -262,6 +263,106 @@ def EndSimulationAction(self):
         self.user_output.stats.write_data_if_requested()
 
 
+"""
+    It is feasible to get callback every Run, Event, Track, Step in the python side.
+    However, it is VERY time consuming. For SteppingAction, expect large performance drop.
+    It could be however useful for prototyping or tests.
+
+    it requires "trampoline functions" on the cpp side.
+
+    # it is feasible but very slow !
+    def SteppingAction(self, step, touchable):
+        g4.GateSimulationStatisticsActor.SteppingAction(self, step, touchable)
+        do_something()
+"""
+
+
+class ActorOutputKillAccordingProcessesActor(ActorOutputBase):
+
+    def __init__(self, *args, **kwargs):
+        super().__init__(*args, **kwargs)
+        self.number_of_killed_particles = 0
+
+    def get_processed_output(self):
+        d = {}
+        d["particles killed"] = self.number_of_killed_particles
+        return d
+
+    def __str__(self):
+        s = ""
+        for k, v in self.get_processed_output().items():
+            s = k + ": " + str(v)
+            s += "\n"
+        return s
+
+
+class KillAccordingProcessesActor(ActorBase, g4.GateKillAccordingProcessesActor):
+    # hints for IDE
+    processes_to_kill: list
+    is_rayleigh_an_interaction: bool
+
+    """
+    This actor enables the user to kill particles according to one or more processes which occur in a volume. If the user
+    wants to kill a particle whenever a proces occurs (except transportation), the "all" option is available.
+    """
+
+    user_info_defaults = {
+        "processes_to_kill": (
+            [],
+            {
+                "doc": "If a processes belonging to this list occured, the particle and its potential secondaries are killed. the variable all can be set up to kill a particle if an interaction occured."
+            },
+        ),
+        "is_rayleigh_an_interaction": (
+            True,
+            {
+                "doc": "Specific case to be faster. If a user wants to kill all interactions which implies an energy loss, this boolean enables to not account Rayleigh process as an interaction"
+            },
+        ),
+    }
+
+    """
+    If a particle, not generated or generated within the volume at which our actor is attached, crosses the volume
+    without interaction, the particle is killed.
+    """
+
+    def __init__(self, *args, **kwargs):
+        ActorBase.__init__(self, *args, **kwargs)
+        self._add_user_output(
+            ActorOutputKillAccordingProcessesActor, "kill_interacting_particles"
+        )
+        self.__initcpp__()
+        self.number_of_killed_particles = 0
+
+    def __initcpp__(self):
+        g4.GateKillAccordingProcessesActor.__init__(self, self.user_info)
+        self.AddActions(
+            {
+                "BeginOfRunAction",
+                "BeginOfEventAction",
+                "PreUserTrackingAction",
+                "SteppingAction",
+                "EndSimulationAction",
+            }
+        )
+
+    def initialize(self):
+        ActorBase.initialize(self)
+        self.InitializeUserInfo(self.user_info)
+        self.InitializeCpp()
+        if len(self.user_info.processes_to_kill) == 0:
+            fatal("You have to select at least one process ! ")
+
+    def EndSimulationAction(self):
+        self.user_output.kill_interacting_particles.number_of_killed_particles = (
+            self.number_of_killed_particles
+        )
+
+    def __str__(self):
+        s = self.user_output["kill_non_interacting_particles"].__str__()
+        return s
+
+
 class KillActor(ActorBase, g4.GateKillActor):
     """
     Actor which kill a particle entering in a volume with the following attached actor.
@@ -442,6 +543,8 @@ def initialize(self):
 process_cls(ActorOutputStatisticsActor)
 process_cls(SimulationStatisticsActor)
 process_cls(KillActor)
+process_cls(ActorOutputKillAccordingProcessesActor)
+process_cls(KillAccordingProcessesActor)
 process_cls(SplittingActorBase)
 process_cls(ComptSplittingActor)
 process_cls(BremSplittingActor)

opengate/managers.py

@@ -90,6 +90,7 @@
 from .actors.miscactors import (
     SimulationStatisticsActor,
     KillActor,
+    KillAccordingProcessesActor,
     SplittingActorBase,
     ComptSplittingActor,
     BremSplittingActor,
@@ -124,6 +125,7 @@
     "ARFTrainingDatasetActor": ARFTrainingDatasetActor,
     "SimulationStatisticsActor": SimulationStatisticsActor,
     "KillActor": KillActor,
+    "KillAccordingProcessesActor": KillAccordingProcessesActor,
     "BremSplittingActor": BremSplittingActor,
     "ComptSplittingActor": ComptSplittingActor,
     "DigitizerAdderActor": DigitizerAdderActor,
@@ -1195,6 +1197,9 @@ def dump_volume_tree(self):
             s += len(pre) * " " + f"{node.name}\n"
         return s
 
+    def get_volume_tree(self):
+        return self.volume_tree_root
+
     def print_volume_tree(self):
         print(self.dump_volume_tree())