Dev fluence actor

core/opengate_core/opengate_lib/GateFluenceActor.cpp

@@ -33,6 +33,7 @@ void GateFluenceActor::InitializeUserInfo(py::dict &user_info) {
   // IMPORTANT: call the base class method
   GateVActor::InitializeUserInfo(user_info);
   fTranslation = DictGetG4ThreeVector(user_info, "translation");
+  fHitType = DictGetStr(user_info, "hit_type");
 }
 
 void GateFluenceActor::InitializeCpp() {
@@ -41,6 +42,11 @@ void GateFluenceActor::InitializeCpp() {
   // Create the image pointer
   // (the size and allocation will be performed on the py side)
   cpp_fluence_image = Image3DType::New();
+
+  // Create sum_tracks image if needed based on scoring mode
+  if (fFluenceScoringMode == "sum_tracks") {
+    cpp_fluence_sum_tracks_image = Image3DType::New();
+  }
 }
 
 void GateFluenceActor::BeginOfEventAction(const G4Event *event) {
@@ -50,43 +56,97 @@ void GateFluenceActor::BeginOfEventAction(const G4Event *event) {
 
 void GateFluenceActor::BeginOfRunActionMasterThread(int run_id) {
   // Important ! The volume may have moved, so we (re-)attach each run
-  AttachImageToVolume<Image3DType>(cpp_fluence_image, fPhysicalVolumeName,
-                                   fTranslation);
+
+  if (fFluenceScoringMode == "fluence") {
+    AttachImageToVolume<Image3DType>(cpp_fluence_image, fPhysicalVolumeName,
+                                     fTranslation);
+    auto sp = cpp_fluence_image->GetSpacing();
+    fVoxelVolume = sp[0] * sp[1] * sp[2];
+  } else if (fFluenceScoringMode == "sum_tracks") {
+    AttachImageToVolume<Image3DType>(cpp_fluence_sum_tracks_image,
+                                     fPhysicalVolumeName, fTranslation);
+    auto sp = cpp_fluence_sum_tracks_image->GetSpacing();
+    fVoxelVolume = sp[0] * sp[1] * sp[2];
+  }
+
   NbOfEvent = 0;
 }
 
+void GateFluenceActor::GetVoxelPosition(G4Step *step, G4ThreeVector &position,
+                                        bool &isInside,
+                                        Image3DType::IndexType &index,
+                                        Image3DType::Pointer &image) const {
+  auto preGlobal = step->GetPreStepPoint()->GetPosition();
+  auto postGlobal = step->GetPostStepPoint()->GetPosition();
+  auto touchable = step->GetPreStepPoint()->GetTouchable();
+
+  // consider random position between pre and post
+  if (fHitType == "pre") {
+    position = preGlobal;
+  }
+  if (fHitType == "post") {
+    position = postGlobal;
+  }
+  if (fHitType == "random") {
+    auto x = G4UniformRand();
+    auto direction = postGlobal - preGlobal;
+    position = preGlobal + x * direction;
+  }
+  if (fHitType == "middle") {
+    auto direction = postGlobal - preGlobal;
+    position = preGlobal + 0.5 * direction;
+  }
+
+  auto localPosition =
+      touchable->GetHistory()->GetTransform(0).TransformPoint(position);
+
+  // convert G4ThreeVector to itk PointType
+  Image3DType::PointType point;
+  point[0] = localPosition[0];
+  point[1] = localPosition[1];
+  point[2] = localPosition[2];
+
+  isInside = image->TransformPhysicalPointToIndex(point, index);
+}
+
 void GateFluenceActor::SteppingAction(G4Step *step) {
-  // same method to consider only entering tracks
-  if (step->GetPreStepPoint()->GetStepStatus() == fGeomBoundary) {
-    // the pre-position is at the edge
-    auto preGlobal = step->GetPreStepPoint()->GetPosition();
-    auto dir = step->GetPreStepPoint()->GetMomentumDirection();
-    auto touchable = step->GetPreStepPoint()->GetTouchable();
-
-    // consider position in the local volume, slightly shifted by 0.1 nm because
-    // otherwise, it can be considered as outside the volume by isInside.
-    auto position = preGlobal + 0.1 * CLHEP::nm * dir;
-    auto localPosition =
-        touchable->GetHistory()->GetTransform(0).TransformPoint(position);
-
-    // convert G4ThreeVector to itk PointType
-    Image3DType::PointType point;
-    point[0] = localPosition[0];
-    point[1] = localPosition[1];
-    point[2] = localPosition[2];
-
-    // get weight
-    auto w = step->GetTrack()->GetWeight();
-
-    // get pixel index
-    Image3DType::IndexType index;
-    bool isInside =
-        cpp_fluence_image->TransformPhysicalPointToIndex(point, index);
-
-    // set value
-    if (isInside) {
+
+  const auto event_id =
+      G4RunManager::GetRunManager()->GetCurrentEvent()->GetEventID();
+  auto preGlobal = step->GetPreStepPoint()->GetPosition();
+  auto postGlobal = step->GetPostStepPoint()->GetPosition();
+
+  // FIXME If the volume has multiple copy, touchable->GetCopyNumber(0) ?
+
+  // Get the voxel index
+  G4ThreeVector position;
+  bool isInside;
+  Image3DType::IndexType index;
+  if (fFluenceScoringMode == "fluence") {
+    GetVoxelPosition(step, position, isInside, index, cpp_fluence_image);
+  } else if (fFluenceScoringMode == "sum_tracks") {
+    GetVoxelPosition(step, position, isInside, index,
+                     cpp_fluence_sum_tracks_image);
+  }
+
+  if (isInside) {
+
+    auto w = step->GetPreStepPoint()->GetWeight();
+    double step_length = step->GetStepLength() / CLHEP::mm * w;
+
+    {
       G4AutoLock FluenceMutex(&SetPixelFluenceMutex);
-      ImageAddValue<Image3DType>(cpp_fluence_image, index, w);
-    } // else : outside the image
+      // Score based on selected mode
+      if (fFluenceScoringMode == "fluence") {
+        // Score fluence only at geometric boundaries
+        if (step->GetPreStepPoint()->GetStepStatus() == fGeomBoundary) {
+          ImageAddValue<Image3DType>(cpp_fluence_image, index, w);
+        }
+      } else if (fFluenceScoringMode == "sum_tracks") {
+        // Score track length density for all steps
+        ImageAddValue<Image3DType>(cpp_fluence_sum_tracks_image, index,
+                                   step_length / fVoxelVolume);
+      }
+    }
   }
 }

core/opengate_core/opengate_lib/GateFluenceActor.h

@@ -39,6 +39,20 @@ class GateFluenceActor : public GateVActor {
 
   inline void SetPhysicalVolumeName(std::string s) { fPhysicalVolumeName = s; }
 
+  // Image type is 3D float by default
+  void SetSumTracksFlag(const bool b) { fSumTracksFlag = b; }
+
+  bool GetSumTracksFlag() const { return fSumTracksFlag; }
+
+  // Set the fluence scoring mode
+  inline void SetFluenceScoringMode(std::string mode) {
+    fFluenceScoringMode = mode;
+  }
+
+  inline std::string GetFluenceScoringMode() const {
+    return fFluenceScoringMode;
+  }
+
   int NbOfEvent = 0;
 
   // Image type is 3D float by default
@@ -48,13 +62,26 @@ class GateFluenceActor : public GateVActor {
   using Size4DType = Image4DType::SizeType;
   Size4DType size_4D;
 
+  void GetVoxelPosition(G4Step *step, G4ThreeVector &position, bool &isInside,
+                        Image3DType::IndexType &index,
+                        Image3DType::Pointer &image) const;
+
   // The image is accessible on py side (shared by all threads)
   Image3DType::Pointer cpp_fluence_image;
+  Image3DType::Pointer cpp_fluence_sum_tracks_image;
+
+  // Option: Is the fluence as sum of the tracks to be scored?
+  bool fSumTracksFlag{};
+
+  // store the voexl volume for later use (for example to compute dose from
+  // fluence)
+  double fVoxelVolume{};
 
 private:
   std::string fPhysicalVolumeName;
   G4ThreeVector fTranslation;
   std::string fHitType;
+  std::string fFluenceScoringMode = "sum_tracks";
 };
 
 #endif // GateFluenceActor_h

core/opengate_core/opengate_lib/pyGateFluenceActor.cpp

@@ -39,6 +39,12 @@ void init_GateFluenceActor(py::module &m) {
            &GateFluenceActor::EndOfRunActionMasterThread)
       .def("GetPhysicalVolumeName", &GateFluenceActor::GetPhysicalVolumeName)
       .def("SetPhysicalVolumeName", &GateFluenceActor::SetPhysicalVolumeName)
+      .def("SetSumTracksFlag", &GateFluenceActor::SetSumTracksFlag)
+      .def("GetSumTracksFlag", &GateFluenceActor::GetSumTracksFlag)
+      .def("SetFluenceScoringMode", &GateFluenceActor::SetFluenceScoringMode)
+      .def("GetFluenceScoringMode", &GateFluenceActor::GetFluenceScoringMode)
       .def_readwrite("NbOfEvent", &GateFluenceActor::NbOfEvent)
-      .def_readwrite("cpp_fluence_image", &GateFluenceActor::cpp_fluence_image);
+      .def_readwrite("cpp_fluence_image", &GateFluenceActor::cpp_fluence_image)
+      .def_readwrite("cpp_fluence_sum_tracks_image",
+                     &GateFluenceActor::cpp_fluence_sum_tracks_image);
 }

opengate/actors/doseactors.py

@@ -1599,13 +1599,14 @@ def EndSimulationAction(self):
 
 class FluenceActor(VoxelDepositActor, g4.GateFluenceActor):
     """
-    FluenceActor: compute a 3D map of fluence
+    FluenceActor: compute a 3D map of fluence or track length density
     FIXME: add scatter order and uncertainty
     """
 
     # hints for IDE
     uncertainty: bool
     scatter: bool
+    fluence_scoring_mode: str
 
     user_info_defaults = {
         "uncertainty": (
@@ -1620,11 +1621,23 @@ class FluenceActor(VoxelDepositActor, g4.GateFluenceActor):
                 "doc": "FIXME",
             },
         ),
+        "fluence_scoring_mode": (
+            "sum_tracks",
+            {
+                "doc": "Choose between 'fluence' (particle count at voxel boundaries) or 'sum_tracks' (track length density per voxel)",
+                "allowed_values": ("fluence", "sum_tracks"),
+            },
+        ),
     }
 
     user_output_config = {
         "fluence": {
             "actor_output_class": ActorOutputSingleImage,
+            "active": False,
+        },
+        "sum_tracks": {
+            "actor_output_class": ActorOutputSingleImage,
+            "active": True,
         },
     }
 
@@ -1652,21 +1665,54 @@ def initialize(self):
             fatal("FluenceActor : uncertainty and scatter not implemented yet")
 
         self.InitializeUserInfo(self.user_info)
+
+        # Set the fluence scoring mode
+        if self.fluence_scoring_mode == "sum_tracks":
+            self.user_output.sum_tracks.set_active(True)
+            self.user_output.fluence.set_active(False)
+        elif self.fluence_scoring_mode == "fluence":
+            self.user_output.fluence.set_active(True)
+            self.user_output.sum_tracks.set_active(False)
+        else:
+            fatal(
+                f"FluenceActor: unknown fluence_scoring_mode '{self.fluence_scoring_mode}'"
+            )
+
+        self.SetFluenceScoringMode(self.fluence_scoring_mode)
         # Set the physical volume name on the C++ side
         self.SetPhysicalVolumeName(self.get_physical_volume_name())
         self.InitializeCpp()
 
     def BeginOfRunActionMasterThread(self, run_index):
-        self.prepare_output_for_run("fluence", run_index)
-        self.push_to_cpp_image("fluence", run_index, self.cpp_fluence_image)
+        if self.user_output.fluence.get_active():
+            self.prepare_output_for_run("fluence", run_index)
+            self.push_to_cpp_image("fluence", run_index, self.cpp_fluence_image)
+
+        if self.user_output.sum_tracks.get_active():
+            self.prepare_output_for_run("sum_tracks", run_index)
+            self.push_to_cpp_image(
+                "sum_tracks", run_index, self.cpp_fluence_sum_tracks_image
+            )
+
         g4.GateFluenceActor.BeginOfRunActionMasterThread(self, run_index)
 
     def EndOfRunActionMasterThread(self, run_index):
-        self.fetch_from_cpp_image("fluence", run_index, self.cpp_fluence_image)
-        self._update_output_coordinate_system("fluence", run_index)
-        self.user_output.fluence.store_meta_data(
-            run_index, number_of_samples=self.NbOfEvent
-        )
+        if self.user_output.fluence.get_active():
+            self.fetch_from_cpp_image("fluence", run_index, self.cpp_fluence_image)
+            self._update_output_coordinate_system("fluence", run_index)
+            self.user_output.fluence.store_meta_data(
+                run_index, number_of_samples=self.NbOfEvent
+            )
+
+        if self.user_output.sum_tracks.get_active():
+            self.fetch_from_cpp_image(
+                "sum_tracks", run_index, self.cpp_fluence_sum_tracks_image
+            )
+            self._update_output_coordinate_system("sum_tracks", run_index)
+            self.user_output.sum_tracks.store_meta_data(
+                run_index, number_of_samples=self.NbOfEvent
+            )
+
         VoxelDepositActor.EndOfRunActionMasterThread(self, run_index)
         return 0