Add option to disable default reactions table

OpenGATE · May 17, 2024 · edee2ec · edee2ec
1 parent 299e608
commit edee2ec
Show file tree

Hide file tree

Showing 3 changed files with 36 additions and 27 deletions.
diff --git a/core/opengate_core/opengate_lib/GateChemistryActor.cpp b/core/opengate_core/opengate_lib/GateChemistryActor.cpp
@@ -24,6 +24,7 @@
 #include <pybind11/pybind11.h>
 #include <pybind11/pytypes.h>
 
+#include "GateHelpers.h"
 #include "GateHelpersDict.h"
 
 #include "../g4_bindings/chemistryadaptator.h"
@@ -40,6 +41,7 @@ GateChemistryActor::GateChemistryActor(pybind11::dict &user_info)
   _timeStepModelStr = DictGetStr(user_info, "timestep_model");
   _endTime = DictGetDouble(user_info, "end_time");
   _reactions = getReactionInputs(user_info, "reactions");
+	_keepDefaultReactions = DictGetBool(user_info, "default_reactions");
 	_moleculeCounterVerbose = DictGetInt(user_info, "molecule_counter_verbose");
 
   setTimeBinsCount(DictGetInt(user_info, "time_bins_count"));
@@ -61,31 +63,10 @@ GateChemistryActor::GateChemistryActor(pybind11::dict &user_info)
 
   G4Scheduler::Instance()->SetEndTime(_endTime);
 
-  {
-    auto constructReactionTable =
-        [&reactions = _reactions](G4DNAMolecularReactionTable *reactionTable) {
-          reactionTable->Reset();
-
-          for (auto const &reaction : reactions) {
-            double rate =
-                reaction.rate * (1e-3 * CLHEP::m3 / (CLHEP::mole * CLHEP::s));
-            auto *reactionData = new G4DNAMolecularReactionData(
-                rate, reaction.reactants[0], reaction.reactants[1]);
-            for (auto const &product : reaction.products)
-              if (product != "H2O")
-                reactionData->AddProduct(product);
-            reactionData->ComputeEffectiveRadius();
-            reactionData->SetReactionType(reaction.type);
-
-            reactionTable->SetReaction(reactionData);
-          }
-
-          reactionTable->PrintTable();
-        };
-
-    ChemistryAdaptator<G4EmDNAChemistry_option3>::setConstructReactionTableHook(
-        constructReactionTable);
-  }
+  if (!_reactions.empty())
+    setupConstructReactionTableHook();
+  else if (!_keepDefaultReactions)
+    Fatal("Disabling default reactions requires to provide reactions");
 
   auto *chemistryList =
       ChemistryAdaptator<G4EmDNAChemistry_option3>::getChemistryList();
@@ -103,7 +84,6 @@ void GateChemistryActor::Initialize(G4HCofThisEvent *hce) {
 void GateChemistryActor::EndSimulationAction() {}
 
 void GateChemistryActor::EndOfRunAction(G4Run const *) {
-	G4cerr << "~~ DEBUG EndOfRunAction" << G4endl;
   for (auto &[molecule, map] : _speciesInfoPerTime) {
     for (auto &[time, data] : map) {
       data.g /= _nbEvents;   // mean value of g
@@ -113,7 +93,6 @@ void GateChemistryActor::EndOfRunAction(G4Run const *) {
 }
 
 void GateChemistryActor::EndOfEventAction(G4Event const *) {
-	G4cerr << "~~ DEBUG EndOfEventAction" << G4endl;
   auto *moleculeCounter = G4MoleculeCounter::Instance();
 
   if (not G4EventManager::GetEventManager()
@@ -235,3 +214,30 @@ GateChemistryActor::getReactionInputs(pybind11::dict &user_info,
 
   return reactionInputs;
 }
+
+void GateChemistryActor::setupConstructReactionTableHook() {
+  auto constructReactionTable =
+    [&reactions = _reactions, &keepDefaultReactions = _keepDefaultReactions]
+    (G4DNAMolecularReactionTable *reactionTable) {
+      if (!keepDefaultReactions) reactionTable->Reset();
+
+      for (auto const &reaction : reactions) {
+        double rate =
+          reaction.rate * (1e-3 * CLHEP::m3 / (CLHEP::mole * CLHEP::s));
+        auto *reactionData = new G4DNAMolecularReactionData(
+          rate, reaction.reactants[0], reaction.reactants[1]);
+        for (auto const &product : reaction.products)
+          if (product != "H2O")
+            reactionData->AddProduct(product);
+        reactionData->ComputeEffectiveRadius();
+        reactionData->SetReactionType(reaction.type);
+
+        reactionTable->SetReaction(reactionData);
+      }
+
+      reactionTable->PrintTable();
+    };
+
+  ChemistryAdaptator<G4EmDNAChemistry_option3>::setConstructReactionTableHook(
+    constructReactionTable);
+}
diff --git a/core/opengate_core/opengate_lib/GateChemistryActor.h b/core/opengate_core/opengate_lib/GateChemistryActor.h
@@ -60,6 +60,7 @@ class GateChemistryActor : public GateVActor {
 protected:
   static ReactionInputs getReactionInputs(pybind11::dict &user_info,
                                           std::string const &key);
+	void setupConstructReactionTableHook();
 
 private:
   SpeciesMap _speciesInfoPerTime;
@@ -72,6 +73,7 @@ class GateChemistryActor : public GateVActor {
   int _moleculeCounterVerbose = 0;
   std::string _timeStepModelStr = "IRT";
   double _endTime;
+	bool _keepDefaultReactions = true;
   std::vector<ReactionInput> _reactions;
 };
 

diff --git a/opengate/actors/chemistryactors.py b/opengate/actors/chemistryactors.py
@@ -14,6 +14,7 @@ def set_default_user_info(user_info):
         user_info.output = "output.root"
         user_info.timestep_model = "IRT"
         user_info.reactions = []
+        user_info.default_reactions = True
         user_info.end_time = 1 * g4_units.s
         user_info.time_bins_count = 10
         user_info.molecule_counter_verbose = 0