Simple geometrical splitting actor

core/opengate_core/opengate_core.cpp

@@ -259,6 +259,8 @@ void init_GateARFTrainingDatasetActor(py::module &m);
 
 void init_GateKillActor(py::module &);
 
+void init_GateSurfaceSplittingActor(py::module &);
+
 void init_itk_image(py::module &);
 
 void init_GateImageNestedParameterisation(py::module &);
@@ -505,6 +507,7 @@ PYBIND11_MODULE(opengate_core, m) {
   init_GateARFActor(m);
   init_GateARFTrainingDatasetActor(m);
   init_GateKillActor(m);
+  init_GateSurfaceSplittingActor(m);
   init_GateDigiAttributeManager(m);
   init_GateVDigiAttribute(m);
   init_GateExceptionHandler(m);

core/opengate_core/opengate_lib/GateSurfaceSplittingActor.cpp

@@ -0,0 +1,96 @@
+/* --------------------------------------------------
+   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 "GateSurfaceSplittingActor.h"
+#include "G4LogicalVolumeStore.hh"
+#include "G4ios.hh"
+#include "GateHelpers.h"
+#include "GateHelpersDict.h"
+
+GateSurfaceSplittingActor::GateSurfaceSplittingActor(py::dict &user_info)
+    : GateVActor(user_info, false) {
+  fActions.insert("StartSimulationAction");
+  fActions.insert("SteppingAction");
+  fActions.insert("PreUserTrackingAction");
+  fMotherVolumeName = DictGetStr(user_info, "mother");
+  fWeightThreshold = DictGetBool(user_info, "weight_threshold");
+  fSplittingFactor = DictGetInt(user_info, "splitting_factor");
+  fSplitEnteringParticles = DictGetBool(user_info, "split_entering_particles");
+  fSplitExitingParticles = DictGetBool(user_info, "split_exiting_particles");
+}
+
+void GateSurfaceSplittingActor::ActorInitialize() {}
+
+void GateSurfaceSplittingActor::StartSimulationAction() {
+  fNbOfKilledParticles = 0;
+}
+
+void GateSurfaceSplittingActor::PreUserTrackingAction(const G4Track *track) {
+
+  fIsFirstStep = true;
+}
+
+void GateSurfaceSplittingActor::SteppingAction(G4Step *step) {
+  auto track = step->GetTrack();
+  auto weight = track->GetWeight();
+
+  if (weight >= fWeightThreshold) {
+    if (fSplitEnteringParticles) {
+      G4String logicalVolumeNamePreStep = step->GetPreStepPoint()
+                                              ->GetPhysicalVolume()
+                                              ->GetLogicalVolume()
+                                              ->GetName();
+      if (((fIsFirstStep) && (step->GetPreStepPoint()->GetStepStatus() == 1) &&
+           (logicalVolumeNamePreStep == fMotherVolumeName)) ||
+          ((fIsFirstStep) &&
+           (track->GetLogicalVolumeAtVertex()->GetName() !=
+            logicalVolumeNamePreStep) &&
+           (track->GetLogicalVolumeAtVertex()->GetName() !=
+            fMotherVolumeName))) {
+        G4ThreeVector position = step->GetPreStepPoint()->GetPosition();
+        G4ThreeVector momentum = step->GetPreStepPoint()->GetMomentum();
+        G4double ekin = step->GetPreStepPoint()->GetKineticEnergy();
+
+        const G4DynamicParticle *particleType = track->GetDynamicParticle();
+        G4double time = step->GetPreStepPoint()->GetGlobalTime();
+        G4TrackVector *trackVector = step->GetfSecondary();
+
+        for (int i = 0; i < fSplittingFactor - 1; i++) {
+          G4DynamicParticle *particleTypeToAdd =
+              new G4DynamicParticle(*particleType);
+          G4Track *clone = new G4Track(particleTypeToAdd, time, position);
+          clone->SetKineticEnergy(ekin);
+          clone->SetMomentumDirection(momentum);
+          clone->SetWeight(weight / fSplittingFactor);
+          trackVector->push_back(clone);
+        }
+        step->GetPreStepPoint()->SetWeight(weight / fSplittingFactor);
+        step->GetPostStepPoint()->SetWeight(weight / fSplittingFactor);
+        track->SetWeight(weight / fSplittingFactor);
+      }
+    }
+
+    if (fSplitExitingParticles) {
+      G4String logicalVolumeNamePostStep = step->GetPostStepPoint()
+                                               ->GetPhysicalVolume()
+                                               ->GetLogicalVolume()
+                                               ->GetName();
+      if (std::find(fListOfVolumeAncestor.begin(), fListOfVolumeAncestor.end(),
+                    logicalVolumeNamePostStep) != fListOfVolumeAncestor.end()) {
+        G4TrackVector *trackVector = step->GetfSecondary();
+        for (int i = 0; i < fSplittingFactor - 1; i++) {
+          G4Track *clone = new G4Track(*track);
+          clone->SetWeight(weight / fSplittingFactor);
+          trackVector->push_back(clone);
+        }
+        step->GetPostStepPoint()->SetWeight(weight / fSplittingFactor);
+        track->SetWeight(weight / fSplittingFactor);
+      }
+    }
+  }
+  fIsFirstStep = false;
+}

core/opengate_core/opengate_lib/GateSurfaceSplittingActor.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 GateSurfaceSplittingActor_h
+#define GateSurfaceSplittingActor_h
+
+#include "G4Cache.hh"
+#include "GateVActor.h"
+#include <pybind11/stl.h>
+
+namespace py = pybind11;
+
+class GateSurfaceSplittingActor : public GateVActor {
+
+public:
+  // Constructor
+  GateSurfaceSplittingActor(py::dict &user_info);
+
+  void ActorInitialize() override;
+
+  void StartSimulationAction() override;
+
+  // Main function called every step in attached volume
+  void SteppingAction(G4Step *) override;
+
+  void PreUserTrackingAction(const G4Track *) override;
+
+  G4bool fSplitEnteringParticles = false;
+  G4bool fSplitExitingParticles = false;
+  G4int fSplittingFactor;
+  G4bool fIsFirstStep;
+  G4bool fWeightThreshold;
+  G4String fMotherVolumeName;
+  std::vector<std::string> fListOfVolumeAncestor;
+
+  long fNbOfKilledParticles{};
+};
+
+#endif

core/opengate_core/opengate_lib/pyGateSurfaceSplittingActor.cpp

@@ -0,0 +1,22 @@
+/* --------------------------------------------------
+   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 "GateSurfaceSplittingActor.h"
+
+void init_GateSurfaceSplittingActor(py::module &m) {
+  py::class_<GateSurfaceSplittingActor,
+             std::unique_ptr<GateSurfaceSplittingActor, py::nodelete>,
+             GateVActor>(m, "GateSurfaceSplittingActor")
+      .def(py::init<py::dict &>())
+      .def_readwrite("fListOfVolumeAncestor",
+                     &GateSurfaceSplittingActor::fListOfVolumeAncestor)
+      .def(py::init<py::dict &>());
+}

opengate/actors/actorbuilders.py

@@ -17,6 +17,7 @@
     SourceInfoActor,
     TestActor,
     KillActor,
+    SurfaceSplittingActor,
     BremSplittingActor,
     ComptSplittingActor,
 )
@@ -44,6 +45,7 @@
     ARFTrainingDatasetActor,
     TestActor,
     KillActor,
+    SurfaceSplittingActor,
     BremSplittingActor,
     ComptSplittingActor,
     DynamicGeometryActor,

opengate/actors/miscactors.py

@@ -5,6 +5,7 @@
 import numpy as np
 import time
 import platform
+from anytree import Node, RenderTree
 import opengate_core as g4
 from .base import ActorBase
 from ..exception import fatal
@@ -429,6 +430,38 @@ def __init__(self, user_info):
         g4.GateKillActor.__init__(self, user_info.__dict__)
 
 
+class SurfaceSplittingActor(g4.GateSurfaceSplittingActor, ActorBase):
+    type_name = "SurfaceSplittingActor"
+
+    def set_default_user_info(user_info):
+        ActorBase.set_default_user_info(user_info)
+        user_info.list_of_volume_name = []
+        user_info.splitting_factor = 1
+        user_info.split_entering_particles = False
+        user_info.split_exiting_particles = False
+        user_info.weight_threshold = 0
+
+    def __init__(self, user_info):
+        ActorBase.__init__(self, user_info)
+        g4.GateSurfaceSplittingActor.__init__(self, user_info.__dict__)
+        self.list_of_volume_name = user_info.list_of_volume_name
+        self.user_info.mother = user_info.mother
+
+    def initialize(self, volume_engine=None):
+
+        super().initialize(volume_engine)
+        volume_tree = self.simulation.volume_manager.get_volume_tree()
+        dico_of_volume_tree = {}
+        for pre, _, node in RenderTree(volume_tree):
+            dico_of_volume_tree[str(node.name)] = node
+        volume_name = self.user_info.mother
+        while volume_name != "world":
+            node = dico_of_volume_tree[volume_name]
+            volume_name = node.mother
+            self.list_of_volume_name.append(volume_name)
+        self.fListOfVolumeAncestor = self.list_of_volume_name
+
+
 """
 class ComptonSplittingActor(g4.GateComptonSplittingActor,ActorBase):
     type_name = "ComptonSplittingActor"

opengate/managers.py

@@ -984,6 +984,9 @@ def dump_volume_types(self):
     def dump_material_database_names(self):
         return list(self.material_database.filenames)
 
+    def get_volume_tree(self):
+        return self.volume_tree_root
+
 
 def setter_hook_verbose_level(self, verbose_level):
     log.setLevel(verbose_level)