Merge pull request #1005 from OpenGATE/ffd

Author: dsarrut

.gitignore

@@ -94,6 +94,8 @@ opengate_temporaire
 /opengate/tests/src/test030_save
 /opengate/tests/src/test_log*
 /opengate/tests/src/save_tests
+/opengate/tests/src/external/coresi/logs/*
+/opengate/tests/src/external/coresi/checkpoi*
 
 /opengate/data/isomeric_transition/save
 /opengate/data/isomeric_transition/save2

core/opengate_core/opengate_lib/GateFluenceActor.cpp

@@ -9,8 +9,10 @@
 #define GateFluenceActor_h
 
 #include "G4Cache.hh"
+#include "G4EmCalculator.hh"
 #include "G4VPrimitiveScorer.hh"
 #include "GateVActor.h"
+#include "digitizer/GateDigiAttributeLastProcessDefinedStepInVolumeActor.h"
 #include "itkImage.h"
 #include <iostream>
 #include <pybind11/stl.h>
@@ -29,16 +31,25 @@ class GateFluenceActor : public GateVActor {
 
   // Function called every step in attached volume
   // This where the scoring takes place
+
   void SteppingAction(G4Step *) override;
 
+  void StartSimulationAction() override;
+
   void BeginOfEventAction(const G4Event *event) override;
 
+  void EndOfEventAction(const G4Event *event) override;
+
   void BeginOfRunActionMasterThread(int run_id) override;
 
+  void BeginOfRunAction(const G4Run *run) override;
+
   inline std::string GetPhysicalVolumeName() { return fPhysicalVolumeName; }
 
   inline void SetPhysicalVolumeName(std::string s) { fPhysicalVolumeName = s; }
 
+  int fNbOfEvent;
+
   int NbOfEvent = 0;
 
   // Image type is 3D float by default
@@ -49,7 +60,77 @@ class GateFluenceActor : public GateVActor {
   Size4DType size_4D;
 
   // The image is accessible on py side (shared by all threads)
-  Image3DType::Pointer cpp_fluence_image;
+  Image3DType::Pointer cpp_counts_image;
+  Image3DType::Pointer cpp_counts_compton_image;
+  Image3DType::Pointer cpp_counts_rayleigh_image;
+  Image3DType::Pointer cpp_counts_secondaries_image;
+  Image3DType::Pointer cpp_counts_primaries_image;
+
+  Image3DType::Pointer cpp_energy_image;
+  Image3DType::Pointer cpp_energy_compton_image;
+  Image3DType::Pointer cpp_energy_rayleigh_image;
+  Image3DType::Pointer cpp_energy_secondaries_image;
+  Image3DType::Pointer cpp_energy_primaries_image;
+
+  Image3DType::Pointer cpp_counts_squared_image;
+  Image3DType::Pointer cpp_counts_squared_compton_image;
+  Image3DType::Pointer cpp_counts_squared_rayleigh_image;
+  Image3DType::Pointer cpp_counts_squared_secondaries_image;
+  Image3DType::Pointer cpp_counts_squared_primaries_image;
+
+  Image3DType::Pointer cpp_energy_squared_image;
+  Image3DType::Pointer cpp_energy_squared_compton_image;
+  Image3DType::Pointer cpp_energy_squared_rayleigh_image;
+  Image3DType::Pointer cpp_energy_squared_secondaries_image;
+  Image3DType::Pointer cpp_energy_squared_primaries_image;
+
+  Image3DType::SizeType size_region{};
+
+  struct threadLocalT {
+    G4EmCalculator emcalc;
+    std::vector<double> squared_worker_flatimg;
+    std::vector<int> lastid_worker_flatimg;
+  };
+
+  G4Cache<threadLocalT> fThreadLocalDataCounts;
+  G4Cache<threadLocalT> fThreadLocalDataComptCounts;
+  G4Cache<threadLocalT> fThreadLocalDataRaylCounts;
+  G4Cache<threadLocalT> fThreadLocalDataSecCounts;
+  G4Cache<threadLocalT> fThreadLocalDataPrimCounts;
+  G4Cache<threadLocalT> fThreadLocalDataEnergy;
+  G4Cache<threadLocalT> fThreadLocalDataComptEnergy;
+  G4Cache<threadLocalT> fThreadLocalDataRaylEnergy;
+  G4Cache<threadLocalT> fThreadLocalDataSecEnergy;
+  G4Cache<threadLocalT> fThreadLocalDataPrimEnergy;
+
+  G4bool fCountsSquaredFlag;
+  G4bool fEnergyFlag;
+  G4bool fEnergySquaredFlag;
+  G4bool fSecondaries;
+  GateDigiAttributeLastProcessDefinedStepInVolumeActor *fLastProcessActor;
+  // GateVActor* fLastProcessActor;
+
+  void FlushSquaredValues(threadLocalT &data,
+                          const Image3DType::Pointer &cpp_image);
+  void ScoreSquaredValue(threadLocalT &data,
+                         const Image3DType::Pointer &cpp_image,
+                         const double value, const int event_id,
+                         const Image3DType::IndexType &index);
+  int sub2ind(Image3DType::IndexType index3D);
+  void PrepareLocalDataForRun(threadLocalT &data,
+                              const unsigned int numberOfVoxels);
+
+  bool GetEnergySquaredFlag() const { return fEnergySquaredFlag; }
+
+  void SetEnergySquaredFlag(const bool b) { fEnergySquaredFlag = b; }
+
+  void SetEnergyFlag(const bool b) { fEnergyFlag = b; }
+
+  bool GetEnergyFlag() const { return fEnergyFlag; }
+
+  void SetCountsSquaredFlag(const bool b) { fCountsSquaredFlag = b; }
+
+  bool GetCountsSquaredFlag() const { return fCountsSquaredFlag; }
 
 private:
   std::string fPhysicalVolumeName;

core/opengate_core/opengate_lib/GateFluenceActor.h

@@ -60,14 +60,13 @@ G4VParticleChange *GateScatterSplittingFreeFlightOptn::ApplyFinalStateBiasing(
     const G4BiasingProcessInterface *callingProcess, const G4Track *track,
     const G4Step *step, G4bool &b) {
   // This is the reference pure G4 version
-  // return ApplyFinalStateBiasing_V1_PostStepDoIt(callingProcess, track, step,
-  // b);
+  return ApplyFinalStateBiasing_V1_PostStepDoIt(callingProcess, track, step, b);
 
   // This is a faster (1x5 speedup) version with direct Compton sampling
   // return ApplyFinalStateBiasing_V3_SampleScatter(callingProcess, track, step,
   // b);
-  return ApplyFinalStateBiasing_V4_SampleComptonOnly(callingProcess, track,
-                                                     step, b);
+  // return ApplyFinalStateBiasing_V4_SampleComptonOnly(callingProcess, track,
+  //                                                 step, b);
 
   // Those are tests, not conclusive
   // return ApplyFinalStateBiasing_V2_SampleSecondaries(callingProcess, track,

core/opengate_core/opengate_lib/biasing/GateScatterSplittingFreeFlightOptn.cpp

@@ -39,6 +39,49 @@ void init_GateFluenceActor(py::module &m) {
            &GateFluenceActor::EndOfRunActionMasterThread)
       .def("GetPhysicalVolumeName", &GateFluenceActor::GetPhysicalVolumeName)
       .def("SetPhysicalVolumeName", &GateFluenceActor::SetPhysicalVolumeName)
+      .def("GetEnergySquaredFlag", &GateFluenceActor::GetEnergySquaredFlag)
+      .def("SetEnergySquaredFlag", &GateFluenceActor::SetEnergySquaredFlag)
+      .def("GetEnergyFlag", &GateFluenceActor::GetEnergyFlag)
+      .def("SetEnergyFlag", &GateFluenceActor::SetEnergyFlag)
+      .def("GetCountsSquaredFlag", &GateFluenceActor::GetCountsSquaredFlag)
+      .def("SetCountsSquaredFlag", &GateFluenceActor::SetCountsSquaredFlag)
       .def_readwrite("NbOfEvent", &GateFluenceActor::NbOfEvent)
-      .def_readwrite("cpp_fluence_image", &GateFluenceActor::cpp_fluence_image);
+      .def_readwrite("cpp_counts_image", &GateFluenceActor::cpp_counts_image)
+      .def_readwrite("cpp_energy_image", &GateFluenceActor::cpp_energy_image)
+      .def_readwrite("cpp_counts_squared_image",
+                     &GateFluenceActor::cpp_counts_squared_image)
+      .def_readwrite("cpp_energy_squared_image",
+                     &GateFluenceActor::cpp_energy_squared_image)
+      .def_readwrite("cpp_counts_compton_image",
+                     &GateFluenceActor::cpp_counts_compton_image)
+      .def_readwrite("cpp_energy_compton_image",
+                     &GateFluenceActor::cpp_energy_compton_image)
+      .def_readwrite("cpp_counts_squared_compton_image",
+                     &GateFluenceActor::cpp_counts_squared_compton_image)
+      .def_readwrite("cpp_energy_squared_compton_image",
+                     &GateFluenceActor::cpp_energy_squared_compton_image)
+      .def_readwrite("cpp_counts_rayleigh_image",
+                     &GateFluenceActor::cpp_counts_rayleigh_image)
+      .def_readwrite("cpp_energy_rayleigh_image",
+                     &GateFluenceActor::cpp_energy_rayleigh_image)
+      .def_readwrite("cpp_counts_squared_rayleigh_image",
+                     &GateFluenceActor::cpp_counts_squared_rayleigh_image)
+      .def_readwrite("cpp_energy_squared_rayleigh_image",
+                     &GateFluenceActor::cpp_energy_squared_rayleigh_image)
+      .def_readwrite("cpp_counts_secondaries_image",
+                     &GateFluenceActor::cpp_counts_secondaries_image)
+      .def_readwrite("cpp_energy_secondaries_image",
+                     &GateFluenceActor::cpp_energy_secondaries_image)
+      .def_readwrite("cpp_counts_squared_secondaries_image",
+                     &GateFluenceActor::cpp_counts_squared_secondaries_image)
+      .def_readwrite("cpp_energy_squared_secondaries_image",
+                     &GateFluenceActor::cpp_energy_squared_secondaries_image)
+      .def_readwrite("cpp_counts_primaries_image",
+                     &GateFluenceActor::cpp_counts_primaries_image)
+      .def_readwrite("cpp_energy_primaries_image",
+                     &GateFluenceActor::cpp_energy_primaries_image)
+      .def_readwrite("cpp_counts_squared_primaries_image",
+                     &GateFluenceActor::cpp_counts_squared_primaries_image)
+      .def_readwrite("cpp_energy_squared_primaries_image",
+                     &GateFluenceActor::cpp_energy_squared_primaries_image);
 }

core/opengate_core/opengate_lib/pyGateFluenceActor.cpp

@@ -305,8 +305,31 @@ FluenceActor
 Description
 ~~~~~~~~~~~
 
-This actor scores the particle fluence on a voxel grid, essentially by counting the number of particles passing through each voxel. The FluenceActor will be extended in the future with features to handle scattered radiation, e.g. in cone beam CT imaging.
+This actor scores the particle fluence on a voxel grid, essentially by counting the number of particles passing through each voxel. When a particle enters a voxel, it records either the counts or the kinetic energy of the incoming particle. If a particle is created within a voxel of the geometry to which the Fluence Actor is attached, this particle is not counted as part of the incoming flux for that voxel. A basic example of its usage is provided below:
 
+.. code-block:: python
+
+    fluence_actor = sim.add_actor("FluenceActor", "fluence_actor")
+    fluence_actor.counts_uncertainty.active = True
+    fluence_actor.counts_squared.active = True
+    fluence_actor.energy.active = True
+    fluence_actor.energy_uncertainty.active = True
+    fluence_actor.energy_squared.active = True
+    fluence_actor.output_filename = "test099.mhd"
+    fluence_actor.attached_to = fluence_plane
+    fluence_actor.size = [10, 10, 1]
+    ts = [10 * cm, 10 * cm, 1 * nm]
+    fluence_actor.spacing = [x / y for x, y in zip(ts, fluence_actor.size)]
+    fluence_actor.hit_type = "random"
+
+
+In addition, it is possible to generate separate fluence maps resolving the particle's tracking state and underlying physics processes (`primaries`, `secondaries`, `compton`, and `rayleigh`). These maps are created as follows: particles originating directly from the source without interacting are recorded as `primaries`, while all others are recorded as `secondaries`. Furthermore, if the incoming particle is a gamma photon and its last interaction was either Compton or Rayleigh scattering, the counts (and optionally the photon energy) are also recorded in the corresponding scattering-process maps. At present, pair production is not yet included as a recordable process. This actor is also compatible with the FreeFlightAngularAcceptance variance reduction technique. To enable the recording of these additional maps, simply set the following boolean to True:
+
+.. code-block:: python
+
+    fluence_actor.score_by_process = True
+
+The activation of the squared counts (and energies) and their associated uncertainty maps is handled automatically, according to the global settings defined by the user for counts and energy scoring.
 
 Reference
 ~~~~~~~~~