dose rate as a proper CLI

core/opengate_core/opengate_lib/GateMaterialMuHandler.cpp

@@ -9,21 +9,23 @@
 #include "GateHelpers.h"
 #include "GateMuDatabase.h"
 #include "GateMuTables.h"
+#include <G4AutoLock.hh>
 #include <G4Gamma.hh>
 #include <G4ProcessManager.hh>
 #include <G4ProcessVector.hh>
 #include <G4ProductionCutsTable.hh>
 #include <G4VAtomDeexcitation.hh>
 #include <G4VEmProcess.hh>
 
-// GateMaterialMuHandler *GateMaterialMuHandler::fSingletonMaterialMuHandler =
-// nullptr;
+G4Mutex sMaterialMuHandlerMutex = G4MUTEX_INITIALIZER;
+
 std::map<std::tuple<std::string, double>,
          std::shared_ptr<GateMaterialMuHandler>>
     GateMaterialMuHandler::fInstances;
 
 std::shared_ptr<GateMaterialMuHandler>
 GateMaterialMuHandler::GetInstance(std::string database, double energy_max) {
+  G4AutoLock mutex(&sMaterialMuHandlerMutex);
   // Create a key based on database and energy_max
   auto key = std::make_tuple(database, energy_max);
 
@@ -49,55 +51,64 @@ GateMaterialMuHandler::GateMaterialMuHandler() {
   fEnergyNumber = 40;
   fAtomicShellEnergyMin = 1. * CLHEP::keV;
   fPrecision = 0.01;
-  fLastCouple = nullptr;
-  fLastMuTable = nullptr;
 }
 
 GateMaterialMuHandler::~GateMaterialMuHandler() { delete[] fElementsTable; }
 
-void GateMaterialMuHandler::CheckLastCall(const G4MaterialCutsCouple *couple) {
+void GateMaterialMuHandler::CheckLastCall(const G4MaterialCutsCouple *) {
   if (!fIsInitialized) {
     Initialize();
   }
-  if (couple != fLastCouple) {
-    fLastCouple = couple;
-    fLastMuTable = fCoupleTable[fLastCouple];
+}
+
+GateMuTable *
+GateMaterialMuHandler::GetMuTable(const G4MaterialCutsCouple *couple) {
+  if (!fIsInitialized) {
+    Initialize();
   }
+  thread_local const GateMaterialMuHandler *lastHandler = nullptr;
+  thread_local const G4MaterialCutsCouple *lastCouple = nullptr;
+  thread_local GateMuTable *lastTable = nullptr;
+  if (this == lastHandler && couple == lastCouple && lastTable != nullptr) {
+    return lastTable;
+  }
+  auto it = fCoupleTable.find(couple);
+  if (it != fCoupleTable.end()) {
+    lastHandler = this;
+    lastCouple = couple;
+    lastTable = it->second;
+    return lastTable;
+  }
+  return nullptr;
 }
 
 double GateMaterialMuHandler::GetDensity(const G4MaterialCutsCouple *couple) {
-  CheckLastCall(couple);
-  return fLastMuTable->GetDensity();
+  auto *muTable = GetMuTable(couple);
+  return muTable ? muTable->GetDensity() : 0.0;
 }
 
 double GateMaterialMuHandler::GetMuEnOverRho(const G4MaterialCutsCouple *couple,
                                              double energy) {
-  CheckLastCall(couple);
-  return fLastMuTable->GetMuEnOverRho(energy);
+  auto *muTable = GetMuTable(couple);
+  return muTable ? muTable->GetMuEnOverRho(energy) : 0.0;
 }
 
 double GateMaterialMuHandler::GetMuEn(const G4MaterialCutsCouple *couple,
                                       double energy) {
-  CheckLastCall(couple);
-  return fLastMuTable->GetMuEn(energy);
+  auto *muTable = GetMuTable(couple);
+  return muTable ? muTable->GetMuEn(energy) : 0.0;
 }
 
 double GateMaterialMuHandler::GetMuOverRho(const G4MaterialCutsCouple *couple,
                                            double energy) {
-  CheckLastCall(couple);
-  return fLastMuTable->GetMuOverRho(energy);
+  auto *muTable = GetMuTable(couple);
+  return muTable ? muTable->GetMuOverRho(energy) : 0.0;
 }
 
 double GateMaterialMuHandler::GetMu(const G4MaterialCutsCouple *couple,
                                     double energy) {
-  CheckLastCall(couple);
-  return fLastMuTable->GetMu(energy);
-}
-
-GateMuTable *
-GateMaterialMuHandler::GetMuTable(const G4MaterialCutsCouple *couple) {
-  CheckLastCall(couple);
-  return fLastMuTable;
+  auto *muTable = GetMuTable(couple);
+  return muTable ? muTable->GetMu(energy) : 0.0;
 }
 
 inline double interpolation(double Xa, double Xb, double Ya, double Yb,
@@ -106,6 +117,7 @@ inline double interpolation(double Xa, double Xb, double Ya, double Yb,
 }
 
 void GateMaterialMuHandler::Initialize() {
+  G4AutoLock mutex(&sMaterialMuHandlerMutex);
   if (fIsInitialized)
     return;
 

core/opengate_core/opengate_lib/GateMaterialMuHandler.h

@@ -109,8 +109,6 @@ class GateMaterialMuHandler {
   int fEnergyNumber;
   double fAtomicShellEnergyMin;
   double fPrecision;
-  const G4MaterialCutsCouple *fLastCouple;
-  GateMuTable *fLastMuTable;
 };
 
 #endif

core/opengate_core/opengate_lib/GateMuTables.cpp

@@ -12,10 +12,6 @@ GateMuTable::GateMuTable(const G4MaterialCutsCouple *couple, G4int size) {
   fMu = new double[size];
   fMuEn = new double[size];
   fSize = size;
-  fLastMu = -1.0;
-  fLastMuEn = -1.0;
-  fLastEnergyMu = -1.0;
-  fLastEnergyMuEn = -1.0;
 
   mCouple = couple;
   mDensity = -1;
@@ -44,69 +40,57 @@ inline double interpol(double x1, double x, double x2, double y1, double y2) {
   //   storage
 }
 
-double GateMuTable::GetMuEnOverRho(double energy) {
-  if (energy != fLastEnergyMuEn) {
-    fLastEnergyMuEn = energy;
-
-    energy = log(energy);
-
-    int inf = 0;
-    int sup = fSize - 1;
-    while (sup - inf > 1) {
-      int tmp_bound = (inf + sup) / 2;
-      if (fEnergy[tmp_bound] > energy) {
-        sup = tmp_bound;
-      } else {
-        inf = tmp_bound;
-      }
-    }
-    double e_inf = fEnergy[inf];
-    double e_sup = fEnergy[sup];
+double GateMuTable::GetMuEnOverRho(double energy) const {
+  const double log_energy = log(energy);
 
-    if (energy > e_inf && energy < e_sup) {
-      fLastMuEn = exp(interpol(e_inf, energy, e_sup, fMuEn[inf], fMuEn[sup]));
+  int inf = 0;
+  int sup = fSize - 1;
+  while (sup - inf > 1) {
+    int tmp_bound = (inf + sup) / 2;
+    if (fEnergy[tmp_bound] > log_energy) {
+      sup = tmp_bound;
     } else {
-      fLastMuEn = exp(fMuEn[inf]);
+      inf = tmp_bound;
     }
   }
+  double e_inf = fEnergy[inf];
+  double e_sup = fEnergy[sup];
 
-  return fLastMuEn;
+  if (log_energy > e_inf && log_energy < e_sup) {
+    return exp(interpol(e_inf, log_energy, e_sup, fMuEn[inf], fMuEn[sup]));
+  } else {
+    return exp(fMuEn[inf]);
+  }
 }
 
-double GateMuTable::GetMuEn(double energy) {
+double GateMuTable::GetMuEn(double energy) const {
   return (GetMuEnOverRho(energy) * mDensity);
 }
 
-double GateMuTable::GetMuOverRho(double energy) {
-  if (energy != fLastEnergyMu) {
-    fLastEnergyMu = energy;
-
-    energy = log(energy);
-
-    int inf = 0;
-    int sup = fSize - 1;
-    while (sup - inf > 1) {
-      int tmp_bound = (inf + sup) / 2;
-      if (fEnergy[tmp_bound] > energy) {
-        sup = tmp_bound;
-      } else {
-        inf = tmp_bound;
-      }
-    }
-    double e_inf = fEnergy[inf];
-    double e_sup = fEnergy[sup];
+double GateMuTable::GetMuOverRho(double energy) const {
+  const double log_energy = log(energy);
 
-    if (energy > e_inf && energy < e_sup) {
-      fLastMu = exp(interpol(e_inf, energy, e_sup, fMu[inf], fMu[sup]));
+  int inf = 0;
+  int sup = fSize - 1;
+  while (sup - inf > 1) {
+    int tmp_bound = (inf + sup) / 2;
+    if (fEnergy[tmp_bound] > log_energy) {
+      sup = tmp_bound;
     } else {
-      fLastMu = exp(fMu[inf]);
+      inf = tmp_bound;
     }
   }
+  double e_inf = fEnergy[inf];
+  double e_sup = fEnergy[sup];
 
-  return fLastMu;
+  if (log_energy > e_inf && log_energy < e_sup) {
+    return exp(interpol(e_inf, log_energy, e_sup, fMu[inf], fMu[sup]));
+  } else {
+    return exp(fMu[inf]);
+  }
 }
 
-double GateMuTable::GetMu(double energy) {
+double GateMuTable::GetMu(double energy) const {
   return (GetMuOverRho(energy) * mDensity);
 }
 

core/opengate_core/opengate_lib/GateMuTables.h

@@ -20,13 +20,13 @@ class GateMuTable {
 
   void PutValue(int index, double energy, double mu, double mu_en) const;
 
-  double GetMuEn(double energy);
+  double GetMuEn(double energy) const;
 
-  double GetMuEnOverRho(double energy);
+  double GetMuEnOverRho(double energy) const;
 
-  double GetMu(double energy);
+  double GetMu(double energy) const;
 
-  double GetMuOverRho(double energy);
+  double GetMuOverRho(double energy) const;
 
   const G4MaterialCutsCouple *GetMaterialCutsCouple() const;
 
@@ -48,10 +48,6 @@ class GateMuTable {
   double *fEnergy;
   double *fMu;
   double *fMuEn;
-  double fLastEnergyMu;
-  double fLastEnergyMuEn;
-  double fLastMu;
-  double fLastMuEn;
   G4int fSize;
 };
 

core/opengate_core/opengate_lib/GateTLEDoseActor.cpp

@@ -129,6 +129,7 @@ void GateTLEDoseActor::BeginOfEventAction(const G4Event *event) {
         fMaterialMuHandler =
             GateMaterialMuHandler::GetInstance(fDatabase, 5 * CLHEP::MeV);
       }
+      fMaterialMuHandler->Initialize();
     }
   }