Compton splitting + pseudo transportation

core/opengate_core/opengate_core.cpp

@@ -289,6 +289,10 @@ void init_GateARFTrainingDatasetActor(py::module &m);
 
 void init_GateKillActor(py::module &);
 
+void init_GateKillNonInteractingParticleActor(py::module &);
+
+void init_GateSurfaceSplittingActor(py::module &);
+
 void init_itk_image(py::module &);
 
 void init_GateImageNestedParameterisation(py::module &);
@@ -321,10 +325,10 @@ void init_GateSimulationStatisticsActor(py::module &);
 
 void init_GatePhaseSpaceActor(py::module &);
 
-// void init_GateComptonSplittingActor(py::module &);
-
 void init_GateOptrComptSplittingActor(py::module &m);
 
+void init_GateOptrComptPseudoTransportationActor(py::module &m);
+
 void init_GateBOptrBremSplittingActor(py::module &m);
 
 void init_G4VBiasingOperator(py::module &m);
@@ -535,8 +539,8 @@ PYBIND11_MODULE(opengate_core, m) {
   init_GateLETActor(m);
   init_GateSimulationStatisticsActor(m);
   init_GatePhaseSpaceActor(m);
-  // init_GateComptonSplittingActor(m);
   init_GateBOptrBremSplittingActor(m);
+  init_GateOptrComptPseudoTransportationActor(m);
   init_GateOptrComptSplittingActor(m);
   init_GateHitsCollectionActor(m);
   init_GateMotionVolumeActor(m);
@@ -550,6 +554,8 @@ PYBIND11_MODULE(opengate_core, m) {
   init_GateARFActor(m);
   init_GateARFTrainingDatasetActor(m);
   init_GateKillActor(m);
+  init_GateKillNonInteractingParticleActor(m);
+  init_GateSurfaceSplittingActor(m);
   init_GateDigiAttributeManager(m);
   init_GateVDigiAttribute(m);
   init_GateExceptionHandler(m);

core/opengate_core/opengate_lib/GateKillNonInteractingParticleActor.cpp

@@ -0,0 +1,75 @@
+/* --------------------------------------------------
+   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 "GateKillNonInteractingParticleActor.h"
+#include "G4LogicalVolumeStore.hh"
+#include "G4PhysicalVolumeStore.hh"
+#include "G4ios.hh"
+#include "GateHelpers.h"
+#include "GateHelpersDict.h"
+
+GateKillNonInteractingParticleActor::GateKillNonInteractingParticleActor(
+    py::dict &user_info)
+    : GateVActor(user_info, false) {
+  fActions.insert("StartSimulationAction");
+  fActions.insert("SteppingAction");
+  fActions.insert("PreUserTrackingAction");
+}
+
+void GateKillNonInteractingParticleActor::ActorInitialize() {}
+
+void GateKillNonInteractingParticleActor::StartSimulationAction() {
+  fNbOfKilledParticles = 0;
+}
+
+void GateKillNonInteractingParticleActor::PreUserTrackingAction(
+    const G4Track *track) {
+  fIsFirstStep = true;
+  fKineticEnergyAtTheEntrance = 0;
+  ftrackIDAtTheEntrance = 0;
+  fPassedByTheMotherVolume = false;
+}
+
+void GateKillNonInteractingParticleActor::SteppingAction(G4Step *step) {
+
+  G4String logNameMotherVolume = G4LogicalVolumeStore::GetInstance()
+                                     ->GetVolume(fMotherVolumeName)
+                                     ->GetName();
+  G4String physicalVolumeNamePreStep =
+      step->GetPreStepPoint()->GetPhysicalVolume()->GetName();
+  if ((step->GetTrack()->GetLogicalVolumeAtVertex()->GetName() !=
+       logNameMotherVolume) &&
+      (fIsFirstStep)) {
+    if ((fPassedByTheMotherVolume == false) &&
+        (physicalVolumeNamePreStep == fMotherVolumeName) &&
+        (step->GetPreStepPoint()->GetStepStatus() == 1)) {
+      fPassedByTheMotherVolume = true;
+      fKineticEnergyAtTheEntrance = step->GetPreStepPoint()->GetKineticEnergy();
+      ftrackIDAtTheEntrance = step->GetTrack()->GetTrackID();
+    }
+  }
+
+  G4String logicalVolumeNamePostStep = step->GetPostStepPoint()
+                                           ->GetPhysicalVolume()
+                                           ->GetLogicalVolume()
+                                           ->GetName();
+  if ((fPassedByTheMotherVolume) &&
+      (step->GetPostStepPoint()->GetStepStatus() == 1)) {
+    if (std::find(fListOfVolumeAncestor.begin(), fListOfVolumeAncestor.end(),
+                  logicalVolumeNamePostStep) != fListOfVolumeAncestor.end()) {
+      if ((step->GetTrack()->GetTrackID() == ftrackIDAtTheEntrance) &&
+          (step->GetPostStepPoint()->GetKineticEnergy() ==
+           fKineticEnergyAtTheEntrance)) {
+        auto track = step->GetTrack();
+        track->SetTrackStatus(fStopAndKill);
+        fNbOfKilledParticles++;
+      }
+    }
+  }
+
+  fIsFirstStep = false;
+}

core/opengate_core/opengate_lib/GateKillNonInteractingParticleActor.h

@@ -0,0 +1,42 @@
+/* --------------------------------------------------
+   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 GateKillNonInteractingParticleActor_h
+#define GateKillNonInteractingParticleActor_h
+
+#include "G4Cache.hh"
+#include "GateVActor.h"
+#include <pybind11/stl.h>
+
+namespace py = pybind11;
+
+class GateKillNonInteractingParticleActor : public GateVActor {
+
+public:
+  // Constructor
+  GateKillNonInteractingParticleActor(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;
+
+  std::vector<G4String> fParticlesTypeToKill;
+  G4bool fPassedByTheMotherVolume = false;
+  G4double fKineticEnergyAtTheEntrance = 0;
+  G4int ftrackIDAtTheEntrance = 0;
+  G4bool fIsFirstStep = true;
+  std::vector<std::string> fListOfVolumeAncestor;
+
+  long fNbOfKilledParticles{};
+};
+
+#endif

core/opengate_core/opengate_lib/GateOptnForceFreeFlight.cpp

@@ -0,0 +1,127 @@
+//
+// ********************************************************************
+// * License and Disclaimer                                           *
+// *                                                                  *
+// * The  Geant4 software  is  copyright of the Copyright Holders  of *
+// * the Geant4 Collaboration.  It is provided  under  the terms  and *
+// * conditions of the Geant4 Software License,  included in the file *
+// * LICENSE and available at  http://cern.ch/geant4/license .  These *
+// * include a list of copyright holders.                             *
+// *                                                                  *
+// * Neither the authors of this software system, nor their employing *
+// * institutes,nor the agencies providing financial support for this *
+// * work  make  any representation or  warranty, express or implied, *
+// * regarding  this  software system or assume any liability for its *
+// * use.  Please see the license in the file  LICENSE  and URL above *
+// * for the full disclaimer and the limitation of liability.         *
+// *                                                                  *
+// * This  code  implementation is the result of  the  scientific and *
+// * technical work of the GEANT4 collaboration.                      *
+// * By using,  copying,  modifying or  distributing the software (or *
+// * any work based  on the software)  you  agree  to acknowledge its *
+// * use  in  resulting  scientific  publications,  and indicate your *
+// * acceptance of all terms of the Geant4 Software license.          *
+// ********************************************************************
+//
+#include "GateOptnForceFreeFlight.h"
+#include "G4BiasingProcessInterface.hh"
+#include "G4ILawForceFreeFlight.hh"
+#include "G4Step.hh"
+
+// This operator is used to transport the particle without interaction, and then
+// correct the weight of the particle
+// according to the probablity for the photon to not interact within the matter.
+// To do that, we use a GEANT4 modifed Interaction Law (G4ILawForceFreeFlight),
+// which modify, for the biased process the probability for the interaction to
+// occur :  Never. This occurence is called during the tracking for each step.
+// Here the step is the largest possible, and one step correspond to the path of
+// particle in the media.
+
+GateOptnForceFreeFlight ::GateOptnForceFreeFlight(G4String name)
+    : G4VBiasingOperation(name), fOperationComplete(true) {
+  fForceFreeFlightInteractionLaw =
+      new G4ILawForceFreeFlight("LawForOperation" + name);
+}
+
+GateOptnForceFreeFlight ::~GateOptnForceFreeFlight() {
+  if (fForceFreeFlightInteractionLaw)
+    delete fForceFreeFlightInteractionLaw;
+}
+
+const G4VBiasingInteractionLaw *
+GateOptnForceFreeFlight ::ProvideOccurenceBiasingInteractionLaw(
+    const G4BiasingProcessInterface *,
+    G4ForceCondition &proposeForceCondition) {
+  fOperationComplete = false;
+  proposeForceCondition = Forced;
+  return fForceFreeFlightInteractionLaw;
+}
+
+G4VParticleChange *GateOptnForceFreeFlight ::ApplyFinalStateBiasing(
+    const G4BiasingProcessInterface *callingProcess, const G4Track *track,
+    const G4Step *step, G4bool &forceFinalState) {
+
+  fParticleChange.Initialize(*track);
+  forceFinalState = true;
+  fCountProcess++;
+
+  fProposedWeight *=
+      fWeightChange[callingProcess->GetWrappedProcess()->GetProcessName()];
+
+  if (fRussianRouletteForWeights) {
+    G4int nbOfOrderOfMagnitude = std::log10(fInitialWeight / fMinWeight);
+
+    if ((fProposedWeight < fMinWeight) ||
+        ((fProposedWeight < 0.1 * fInitialWeight) &&
+         (fNbOfRussianRoulette == nbOfOrderOfMagnitude - 1))) {
+      fParticleChange.ProposeTrackStatus(G4TrackStatus::fStopAndKill);
+      fNbOfRussianRoulette = 0;
+      return &fParticleChange;
+    }
+
+    if ((fProposedWeight < 0.1 * fInitialWeight) && (fCountProcess == 4)) {
+      G4double probability = G4UniformRand();
+      for (int i = 2; i <= nbOfOrderOfMagnitude; i++) {
+        G4double RRprobability = 1 / std::pow(10, i);
+        if (fProposedWeight * 1 / std::pow(10, fNbOfRussianRoulette) <
+            10 * RRprobability * fMinWeight) {
+          fParticleChange.ProposeTrackStatus(G4TrackStatus::fStopAndKill);
+          fNbOfRussianRoulette = 0;
+          return &fParticleChange;
+        }
+        if ((fProposedWeight >= RRprobability * fInitialWeight) &&
+            (fProposedWeight < 10 * RRprobability * fInitialWeight)) {
+          if (probability > 10 * RRprobability) {
+            fParticleChange.ProposeTrackStatus(G4TrackStatus::fStopAndKill);
+            fNbOfRussianRoulette = 0;
+            return &fParticleChange;
+          } else {
+            fProposedWeight = fProposedWeight / (10 * RRprobability);
+            fNbOfRussianRoulette = fNbOfRussianRoulette + i - 1;
+          }
+        }
+      }
+    }
+  } else {
+    if (fProposedWeight < fMinWeight) {
+      fParticleChange.ProposeTrackStatus(G4TrackStatus::fStopAndKill);
+      return &fParticleChange;
+    }
+  }
+  fParticleChange.ProposeWeight(fProposedWeight);
+  fOperationComplete = true;
+  return &fParticleChange;
+}
+
+void GateOptnForceFreeFlight ::AlongMoveBy(
+    const G4BiasingProcessInterface *callingProcess, const G4Step *,
+    G4double weightChange)
+
+{
+  G4String processName = callingProcess->GetWrappedProcess()->GetProcessName();
+  if (processName != "Rayl") {
+    fWeightChange[processName] = weightChange;
+  } else {
+    fWeightChange[processName] = 1;
+  }
+}