Actor for time of flight for neutron

.github/workflows/main.yml

@@ -10,7 +10,6 @@ on:
   pull_request:
     paths-ignore:
       - 'docs/**'
-    branches: [ master ]
   schedule:
     - cron: '0 0 * * 0,3'
   workflow_dispatch:

core/opengate_core/opengate_lib/GateFluenceActor.cpp

@@ -5,88 +5,242 @@
    See LICENSE.md for further details
    ------------------------------------ -------------- */
 
-#include "G4Navigator.hh"
+#include "G4EmCalculator.hh"
+#include "G4Gamma.hh"
+#include "G4ParticleDefinition.hh"
+#include "G4RandomTools.hh"
 #include "G4RunManager.hh"
 #include "G4Threading.hh"
+#include "G4Track.hh"
 
+#include "G4HadronInelasticProcess.hh"
 #include "GateFluenceActor.h"
+#include "GateHelpers.h"
 #include "GateHelpersDict.h"
 #include "GateHelpersImage.h"
+#include "GateMaterialMuHandler.h"
 
 #include <iostream>
 #include <itkAddImageFilter.h>
+#include <itkCastImageFilter.h>
+#include <itkImageFileWriter.h>
 #include <itkImageRegionIterator.h>
-
-// Mutex that will be used by thread to write the output image
-G4Mutex SetPixelFluenceMutex = G4MUTEX_INITIALIZER;
-G4Mutex SetNbEventMutexFluence = G4MUTEX_INITIALIZER;
+#include <vector>
 
 GateFluenceActor::GateFluenceActor(py::dict &user_info)
     : GateVActor(user_info, true) {
-
-  // Action for this actor: during stepping
-  fActions.insert("SteppingAction");
-  fActions.insert("BeginOfRunAction");
+  fMultiThreadReady = true;
 }
 
 void GateFluenceActor::InitializeUserInfo(py::dict &user_info) {
   // IMPORTANT: call the base class method
   GateVActor::InitializeUserInfo(user_info);
   fTranslation = DictGetG4ThreeVector(user_info, "translation");
+  fsize = DictGetG4ThreeVector(user_info, "size");
+  fspacing = DictGetG4ThreeVector(user_info, "spacing");
+  Nbbinstime = py::int_(user_info["timebins"]);
+  Nbbinsenergy = py::int_(user_info["energybins"]);
+  foutputname = std::string(py::str(user_info["output_filename"]));
 }
 
 void GateFluenceActor::InitializeCpp() {
   GateVActor::InitializeCpp();
-
-  // Create the image pointer
+  // Create the image pointers
   // (the size and allocation will be performed on the py side)
-  cpp_fluence_image = Image3DType::New();
+  tof_cpp_image = Image3DType::New();
+
+  Image3DType::RegionType region3;
+  Image3DType::SizeType size3;
+  size3[0] = fsize[0]; // Replace with actual size in X
+  size3[1] = fsize[1]; // Replace with actual size in Y
+  size3[2] = fsize[2]; // Replace with actual size in Z
+  region3.SetSize(size3);
+
+  Image3DType::SpacingType spacing3;
+  spacing3[0] = fspacing[0]; // Replace with actual spacing in X
+  spacing3[1] = fspacing[1]; // Replace with actual spacing in Y
+  spacing3[2] = fspacing[2]; // Replace with actual spacing in Z
+
+  Image3DType::PointType origin;
+  origin[0] = fTranslation[0]; // Replace with actual origin in X
+  origin[1] = fTranslation[1]; // Replace with actual origin in Y
+  origin[2] = fTranslation[2]; // Replace with actual origin in Z
+
+  tof_cpp_image->SetRegions(region3);
+  tof_cpp_image->SetSpacing(spacing3);
+  tof_cpp_image->SetOrigin(origin);
+  tof_cpp_image->Allocate();
+  tof_cpp_image->FillBuffer(0);
+
+  output_image = Image2DType::New();
+
+  Image2DType::RegionType region;
+  Image2DType::SizeType size;
+  size[0] = Nbbinsenergy; // Width
+  size[1] = Nbbinstime;   // Height
+  region.SetSize(size);
+
+  Image2DType::SpacingType spacing;
+  spacing[0] = 1.0; // Energy spacing
+  spacing[1] = 1.0; // Time spacing
+
+  output_image->SetRegions(region);
+  output_image->SetSpacing(spacing);
+  output_image->Allocate();
+  output_image->FillBuffer(0); // Initialize with zeros
+
+  norm = 0;
+  nb_inel = 0;
 }
 
-void GateFluenceActor::BeginOfEventAction(const G4Event *event) {
-  G4AutoLock mutex(&SetNbEventMutexFluence);
-  NbOfEvent++;
-}
+void GateFluenceActor::BeginOfEventAction(const G4Event *event) { norm++; }
 
 void GateFluenceActor::BeginOfRunActionMasterThread(int run_id) {
   // Important ! The volume may have moved, so we (re-)attach each run
-  AttachImageToVolume<Image3DType>(cpp_fluence_image, fPhysicalVolumeName,
+  AttachImageToVolume<Image3DType>(tof_cpp_image, fPhysicalVolumeName,
                                    fTranslation);
-  NbOfEvent = 0;
 }
 
 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();
+  // Get the voxel index
+  if (step->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName() ==
+          "neutronInelastic" ||
+      step->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName() ==
+          "protonInelastic") {
+  }
+  if (step->GetTrack()->GetParticleDefinition()->GetParticleName() !=
+      "proton") {
+    return;
+  }
+  auto secondaries = step->GetSecondary();
+  for (size_t i = 0; i < secondaries->size(); i++) {
+    auto secondary = secondaries->at(i);
+    auto secondary_def = secondary->GetParticleDefinition();
+
+    if (secondary_def != G4Gamma::Gamma()) {
+      continue;
+    }
+    if (secondary->GetCreatorProcess()->GetProcessName() != "protonInelastic") {
+      continue;
+    }
+
+    // Get the position of the interaction
+    // auto preposition = step->GetPreStepPoint()->GetPosition();
+    auto position = step->GetPostStepPoint()->GetPosition();
     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;
+    // Get the voxel index
     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) {
-      G4AutoLock FluenceMutex(&SetPixelFluenceMutex);
-      ImageAddValue<Image3DType>(cpp_fluence_image, index, w);
-    } // else : outside the image
+    G4bool isInside =
+        tof_cpp_image->TransformPhysicalPointToIndex(point, index);
+
+    if (!isInside) {
+      return;
+    }
+
+    auto energyPG =
+        secondary->GetKineticEnergy(); // Get the energy of the secondary
+
+    if (energyPG < 0) {
+      std::cerr << "Warning: Negative energy detected." << std::endl;
+      continue;
+    }
+
+    G4double energy_range_EP =
+        10.0 * CLHEP::MeV; // Replace with the actual range of energy values
+
+    Image2DType::IndexType ind;
+    G4double widthbin = (energy_range_EP / Nbbinsenergy);
+
+    // Assign bin for time (ind[1])
+    ind[0] = static_cast<int>(energyPG / widthbin);
+    if (ind[0] >= Nbbinsenergy) {
+      ind[0] = Nbbinsenergy - 1; // Clamp to the maximum bin index
+    }
+
+    // Clamp energyPG to the valid range
+    if (energyPG < 0)
+      energyPG = 0;
+    if (energyPG > energy_range_EP)
+      energyPG = energy_range_EP;
+
+    // time of the fragmentation :
+    creationtime = step->GetPostStepPoint()->GetGlobalTime();
+
+    // time of the gamma emission
+    G4double currenttime = secondary->GetGlobalTime();
+
+    // time of emission
+    G4double time = (currenttime - creationtime) * 1e6; // fs
+
+    if (time == 0) {
+      continue;
+    }
+
+    std::cout << "Current time: " << currenttime << std::endl;
+    std::cout << "Creation time: " << creationtime << std::endl;
+    std::cout << "Time: " << time << std::endl;
+
+    G4double time_max = 1e6;
+    G4double time_min = 1; // to avoid log(0))
+    G4double log_time_min = std::log(time_min);
+    G4double log_time_max = std::log(time_max);
+    G4double log_time_width = (log_time_max - log_time_min) / Nbbinstime;
+
+    // Assign bin for time (log scale)
+    if (time < time_min)
+      time = time_min;
+    if (time > time_max)
+      time = time_max;
+    ind[1] = static_cast<int>((std::log(time) - log_time_min) / log_time_width);
+    if (ind[1] >= Nbbinstime)
+      ind[1] = Nbbinstime - 1; // Sécurité
+
+    std::cout << "Voxel index: " << ind[0] << ", " << ind[1] << std::endl;
+    ImageAddValue<Image2DType>(output_image, ind, 1);
+    nb_inel = nb_inel + 1;
   }
 }
+
+void GateFluenceActor::EndOfRunAction(const G4Run *run) {
+  itk::ImageRegionIterator<Image2DType> it(
+      output_image, output_image->GetLargestPossibleRegion());
+  for (it.GoToBegin(); !it.IsAtEnd(); ++it) {
+    it.Set(it.Get() / norm);
+  }
+  std::cout << "incidentpart" << norm << std::endl;
+  std::cout << "nb gammas" << nb_inel << std::endl;
+}
+
+std::string GateFluenceActor::GetOutputImage() {
+
+  using InputImageType = itk::Image<double, 2>;
+  using OutputImageType = itk::Image<float, 2>;
+  using CastFilterType = itk::CastImageFilter<InputImageType, OutputImageType>;
+
+  try {
+    auto castFilter = CastFilterType::New();
+    castFilter->SetInput(output_image);
+    castFilter->Update();
+
+    using WriterType = itk::ImageFileWriter<OutputImageType>;
+    auto writer = WriterType::New();
+    writer->SetFileName(foutputname);
+    writer->SetInput(castFilter->GetOutput());
+    writer->Update();
+  } catch (const itk::ExceptionObject &e) {
+    std::cerr << "Error during ITK operations: " << e << std::endl;
+    throw;
+  }
+  return foutputname;
+}
+
+int GateFluenceActor::EndOfRunActionMasterThread(int run_id) { return 0; }

core/opengate_core/opengate_lib/GateFluenceActor.h

@@ -8,13 +8,24 @@
 #ifndef GateFluenceActor_h
 #define GateFluenceActor_h
 
+#include "GateDoseActor.h"
+#include "GateMaterialMuHandler.h"
+
 #include "G4Cache.hh"
+#include "G4EmCalculator.hh"
+#include "G4NistManager.hh"
 #include "G4VPrimitiveScorer.hh"
-#include "GateVActor.h"
-#include "itkImage.h"
+
 #include <iostream>
 #include <pybind11/stl.h>
 
+#include "GateVActor.h"
+#include "itkImage.h"
+#include <G4Threading.hh>
+
+#include <itkCastImageFilter.h>
+#include <itkImageFileWriter.h>
+
 namespace py = pybind11;
 
 class GateFluenceActor : public GateVActor {
@@ -23,38 +34,47 @@ class GateFluenceActor : public GateVActor {
   // Constructor
   GateFluenceActor(py::dict &user_info);
 
-  void InitializeCpp() override;
-
   void InitializeUserInfo(py::dict &user_info) override;
 
-  // Function called every step in attached volume
-  // This where the scoring takes place
-  void SteppingAction(G4Step *) override;
+  void InitializeCpp() override;
+
+  void BeginOfRunActionMasterThread(int run_id) override;
 
   void BeginOfEventAction(const G4Event *event) override;
 
-  void BeginOfRunActionMasterThread(int run_id) override;
+  int EndOfRunActionMasterThread(int run_id) override;
 
-  inline std::string GetPhysicalVolumeName() { return fPhysicalVolumeName; }
+  void EndOfRunAction(const G4Run *run) override;
 
-  inline void SetPhysicalVolumeName(std::string s) { fPhysicalVolumeName = s; }
+  std::string GetOutputImage();
 
-  int NbOfEvent = 0;
+  // Main function called every step in attached volume
+  void SteppingAction(G4Step *) override;
 
-  // Image type is 3D float by default
-  typedef itk::Image<float, 3> Image3DType;
-  typedef itk::Image<float, 4> Image4DType;
-  typedef itk::Image<int, 4> ImageInt4DType;
-  using Size4DType = Image4DType::SizeType;
-  Size4DType size_4D;
+  inline std::string GetPhysicalVolumeName() const {
+    return fPhysicalVolumeName;
+  }
 
-  // The image is accessible on py side (shared by all threads)
-  Image3DType::Pointer cpp_fluence_image;
+  inline void SetPhysicalVolumeName(std::string s) { fPhysicalVolumeName = s; }
 
-private:
   std::string fPhysicalVolumeName;
+
+  typedef itk::Image<double, 3> Image3DType;
+  Image3DType::Pointer tof_cpp_image;
+  typedef itk::Image<double, 2> Image2DType;
+  Image2DType::Pointer output_image;
+
   G4ThreeVector fTranslation;
-  std::string fHitType;
+  G4ThreeVector fsize;
+  G4ThreeVector fspacing;
+
+private:
+  G4int norm;
+  G4int nb_inel;
+  G4int Nbbinstime;
+  G4int Nbbinsenergy;
+  G4double creationtime;
+  std::string foutputname;
 };
 
 #endif // GateFluenceActor_h