diff --git a/core/opengate_core/opengate_lib/GateFluenceActor.cpp b/core/opengate_core/opengate_lib/GateFluenceActor.cpp index da0bfc7999..52f2604508 100644 --- a/core/opengate_core/opengate_lib/GateFluenceActor.cpp +++ b/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(cpp_fluence_image, fPhysicalVolumeName, - fTranslation); + + if (fFluenceScoringMode == "fluence") { + AttachImageToVolume(cpp_fluence_image, fPhysicalVolumeName, + fTranslation); + auto sp = cpp_fluence_image->GetSpacing(); + fVoxelVolume = sp[0] * sp[1] * sp[2]; + } else if (fFluenceScoringMode == "sum_tracks") { + AttachImageToVolume(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(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(cpp_fluence_image, index, w); + } + } else if (fFluenceScoringMode == "sum_tracks") { + // Score track length density for all steps + ImageAddValue(cpp_fluence_sum_tracks_image, index, + step_length / fVoxelVolume); + } + } } } diff --git a/core/opengate_core/opengate_lib/GateFluenceActor.h b/core/opengate_core/opengate_lib/GateFluenceActor.h index 00298f3982..f1223a16ba 100644 --- a/core/opengate_core/opengate_lib/GateFluenceActor.h +++ b/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 diff --git a/core/opengate_core/opengate_lib/pyGateFluenceActor.cpp b/core/opengate_core/opengate_lib/pyGateFluenceActor.cpp index fd90006be7..f8af2a276b 100644 --- a/core/opengate_core/opengate_lib/pyGateFluenceActor.cpp +++ b/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); } diff --git a/opengate/actors/doseactors.py b/opengate/actors/doseactors.py index 61a2b209be..97eed5ef2e 100644 --- a/opengate/actors/doseactors.py +++ b/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