diff --git a/core/opengate_core/opengate_lib/GateChemistryActor.cpp b/core/opengate_core/opengate_lib/GateChemistryActor.cpp index da2c4d572..ab0e66a97 100644 --- a/core/opengate_core/opengate_lib/GateChemistryActor.cpp +++ b/core/opengate_core/opengate_lib/GateChemistryActor.cpp @@ -24,6 +24,7 @@ #include #include +#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::setConstructReactionTableHook( - constructReactionTable); - } + if (!_reactions.empty()) + setupConstructReactionTableHook(); + else if (!_keepDefaultReactions) + Fatal("Disabling default reactions requires to provide reactions"); auto *chemistryList = ChemistryAdaptator::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::setConstructReactionTableHook( + constructReactionTable); +} diff --git a/core/opengate_core/opengate_lib/GateChemistryActor.h b/core/opengate_core/opengate_lib/GateChemistryActor.h index 49eb73f5b..fe05deff6 100644 --- 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 _reactions; }; diff --git a/opengate/actors/chemistryactors.py b/opengate/actors/chemistryactors.py index 42ae95cf0..75ec085d7 100644 --- 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