diff --git a/core/opengate_core/g4_bindings/chemistryadaptator.h b/core/opengate_core/g4_bindings/chemistryadaptator.h new file mode 100644 index 000000000..fa0d422d4 --- /dev/null +++ b/core/opengate_core/g4_bindings/chemistryadaptator.h @@ -0,0 +1,58 @@ +/* -------------------------------------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + -------------------------------------------------- */ + +#ifndef chemistryadaptator_h +#define chemistryadaptator_h + +#include + +class G4DNAMolecularReactionTable; +#include + +template class ChemistryAdaptator : public C { +public: + using ConstructReactionTableHook = + std::function; + +public: + ChemistryAdaptator(int verbosity) { + C::SetVerboseLevel(verbosity); + _chemistryLists.push_back(this); + } + + void + ConstructTimeStepModel(G4DNAMolecularReactionTable *reactionTable) override { + C::ConstructTimeStepModel(reactionTable); + reactionTable->PrintTable(); + } + + void + ConstructReactionTable(G4DNAMolecularReactionTable *reactionTable) override { + C::ConstructReactionTable(reactionTable); + if (_constructReactionTableHook) + _constructReactionTableHook(reactionTable); + } + + static C *getChemistryList() { + for (auto *chemistryList : _chemistryLists) { + auto *ptr = dynamic_cast(chemistryList); + if (ptr != nullptr) + return ptr; + } + return nullptr; + } + + template static void setConstructReactionTableHook(T fn) { + _constructReactionTableHook = fn; + } + +private: + inline static ConstructReactionTableHook _constructReactionTableHook; + inline static std::vector _chemistryLists; +}; + +#endif diff --git a/core/opengate_core/g4_bindings/pyG4UserStackingAction.cpp b/core/opengate_core/g4_bindings/pyG4UserStackingAction.cpp new file mode 100644 index 000000000..9d4d4f30c --- /dev/null +++ b/core/opengate_core/g4_bindings/pyG4UserStackingAction.cpp @@ -0,0 +1,30 @@ +/* -------------------------------------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + -------------------------------------------------- */ + +#include + +namespace py = pybind11; + +#include "G4UserStackingAction.hh" + +class PyG4UserStackingAction : public G4UserStackingAction { +public: + /* Inherit the constructors */ + using G4UserStackingAction::G4UserStackingAction; + + void NewStage() override { + PYBIND11_OVERLOAD(void, G4UserStackingAction, NewStage); + } +}; + +void init_G4UserStackingAction(py::module &m) { + + py::class_( + m, "G4UserStackingAction") + .def(py::init()) + .def("NewStage", &G4UserStackingAction::NewStage); +} diff --git a/core/opengate_core/g4_bindings/pyG4VUserActionInitialization.cpp b/core/opengate_core/g4_bindings/pyG4VUserActionInitialization.cpp index 121ce6d20..3b3f96bfd 100644 --- a/core/opengate_core/g4_bindings/pyG4VUserActionInitialization.cpp +++ b/core/opengate_core/g4_bindings/pyG4VUserActionInitialization.cpp @@ -12,6 +12,7 @@ namespace py = pybind11; #include "G4SteppingVerbose.hh" #include "G4UserEventAction.hh" #include "G4UserRunAction.hh" +#include "G4UserStackingAction.hh" #include "G4UserSteppingAction.hh" #include "G4UserTrackingAction.hh" #include "G4VUserActionInitialization.hh" @@ -67,6 +68,11 @@ class PyG4VUserActionInitialization : public G4VUserActionInitialization { G4VUserActionInitialization::SetUserAction(e); } + void SetUserAction(G4UserStackingAction *e) { + // std::cout << "PyG4VUserActionInitialization::SetUserAction" << std::endl; + G4VUserActionInitialization::SetUserAction(e); + } + void SetUserAction(G4UserSteppingAction *e) { // std::cout << "PyG4VUserActionInitialization::SetUserAction" << std::endl; G4VUserActionInitialization::SetUserAction(e); @@ -113,6 +119,10 @@ void init_G4VUserActionInitialization(py::module &m) { (void(G4VUserActionInitialization::*)(G4UserTrackingAction *)) & PyG4VUserActionInitialization::SetUserAction) + .def("SetUserAction", + (void(G4VUserActionInitialization::*)(G4UserStackingAction *)) & + PyG4VUserActionInitialization::SetUserAction) + .def("SetUserAction", (void(G4VUserActionInitialization::*)(G4UserEventAction *)) & PyG4VUserActionInitialization::SetUserAction) diff --git a/core/opengate_core/g4_bindings/pyPhysicsLists.cpp b/core/opengate_core/g4_bindings/pyPhysicsLists.cpp index 578b7691a..472241382 100644 --- a/core/opengate_core/g4_bindings/pyPhysicsLists.cpp +++ b/core/opengate_core/g4_bindings/pyPhysicsLists.cpp @@ -4,6 +4,7 @@ of the GNU Lesser General Public Licence (LGPL) See LICENSE.md for further details -------------------------------------------------- */ +#include #include #include @@ -41,6 +42,14 @@ namespace py = pybind11; #include "G4EmStandardPhysics_option4.hh" #include "G4EmDNAPhysics.hh" +#include "G4EmDNAPhysics_option1.hh" +#include "G4EmDNAPhysics_option2.hh" +#include "G4EmDNAPhysics_option3.hh" +#include "G4EmDNAPhysics_option4.hh" +#include "G4EmDNAPhysics_option5.hh" +#include "G4EmDNAPhysics_option6.hh" +#include "G4EmDNAPhysics_option7.hh" +#include "G4EmDNAPhysics_option8.hh" #include "G4EmLivermorePhysics.hh" #include "G4EmLivermorePolarizedPhysics.hh" #include "G4EmLowEPPhysics.hh" @@ -52,10 +61,17 @@ namespace py = pybind11; #include "G4DecayPhysics.hh" #include "G4RadioactiveDecayPhysics.hh" +#include "G4EmDNAChemistry.hh" +#include "G4EmDNAChemistry_option1.hh" +#include "G4EmDNAChemistry_option2.hh" +#include "G4EmDNAChemistry_option3.hh" + #include "G4VModularPhysicsList.hh" #include "G4VPhysicsConstructor.hh" #include "G4VUserPhysicsList.hh" +#include "chemistryadaptator.h" + // macro for adding physics lists: no parameter // #define ADD_PHYSICS_LIST0(m, plname) \ // py::class_(m, #plname).def(py::init<>()); \ @@ -80,6 +96,13 @@ namespace py = pybind11; std::unique_ptr>(m, #plname) \ .def(py::init()); +// copied from ADD_PHYSICS_CONSTRUCTOR, adapted to chemistry +#define ADD_CHEMISTRY_CONSTRUCTOR(clname) \ + py::class_, G4VPhysicsConstructor, \ + std::unique_ptr, py::nodelete>>( \ + m, #clname) \ + .def(py::init()); + // FIXME ? A bit different for the biasing classe which do not take as argument // a int. Moreover, we need at least the function PhysicsBias to put a bias, so // this constructor needs probably its own function ?. @@ -163,12 +186,25 @@ void init_G4PhysicsLists(py::module &m) { ADD_PHYSICS_CONSTRUCTOR(G4EmLivermorePolarizedPhysics) ADD_PHYSICS_CONSTRUCTOR(G4EmPenelopePhysics) ADD_PHYSICS_CONSTRUCTOR(G4EmDNAPhysics) + ADD_PHYSICS_CONSTRUCTOR(G4EmDNAPhysics_option1) + ADD_PHYSICS_CONSTRUCTOR(G4EmDNAPhysics_option2) + ADD_PHYSICS_CONSTRUCTOR(G4EmDNAPhysics_option3) + ADD_PHYSICS_CONSTRUCTOR(G4EmDNAPhysics_option4) + ADD_PHYSICS_CONSTRUCTOR(G4EmDNAPhysics_option5) + ADD_PHYSICS_CONSTRUCTOR(G4EmDNAPhysics_option6) + ADD_PHYSICS_CONSTRUCTOR(G4EmDNAPhysics_option7) + ADD_PHYSICS_CONSTRUCTOR(G4EmDNAPhysics_option8) ADD_PHYSICS_CONSTRUCTOR(G4OpticalPhysics) ADD_PHYSICS_CONSTRUCTOR_BIASING(G4GenericBiasingPhysics) ADD_PHYSICS_CONSTRUCTOR(G4DecayPhysics) ADD_PHYSICS_CONSTRUCTOR(G4RadioactiveDecayPhysics) + ADD_CHEMISTRY_CONSTRUCTOR(G4EmDNAChemistry) + ADD_CHEMISTRY_CONSTRUCTOR(G4EmDNAChemistry_option1) + ADD_CHEMISTRY_CONSTRUCTOR(G4EmDNAChemistry_option2) + ADD_CHEMISTRY_CONSTRUCTOR(G4EmDNAChemistry_option3) + // sort PL vector std::sort(plList.begin(), plList.end()); } diff --git a/core/opengate_core/opengate_core.cpp b/core/opengate_core/opengate_core.cpp index bff200f36..b3ca95e25 100644 --- a/core/opengate_core/opengate_core.cpp +++ b/core/opengate_core/opengate_core.cpp @@ -87,6 +87,8 @@ void init_G4UserEventAction(py::module &); void init_G4UserTrackingAction(py::module &); +void init_G4UserStackingAction(py::module &); + void init_G4UserSteppingAction(py::module &); void init_G4Track(py::module &); @@ -296,6 +298,10 @@ void init_GateFluenceActor(py::module &m); void init_GateLETActor(py::module &m); +void init_GateChemistryActor(py::module &m); + +void init_GateChemistryLongTimeActor(py::module &m); + void init_GateARFActor(py::module &m); void init_GateARFTrainingDatasetActor(py::module &m); @@ -332,6 +338,8 @@ void init_GateEventAction(py::module &); void init_GateTrackingAction(py::module &); +void init_GateStackingAction(py::module &); + void init_GateSimulationStatisticsActor(py::module &); void init_GatePhaseSpaceActor(py::module &); @@ -426,6 +434,7 @@ PYBIND11_MODULE(opengate_core, m) { init_G4PrimaryVertex(m); init_G4UserEventAction(m); init_G4UserTrackingAction(m); + init_G4UserStackingAction(m); init_G4StepPoint(m); init_G4Track(m); init_G4Step(m); @@ -553,9 +562,12 @@ PYBIND11_MODULE(opengate_core, m) { init_GateRunAction(m); init_GateEventAction(m); init_GateTrackingAction(m); + init_GateStackingAction(m); init_GateDoseActor(m); init_GateFluenceActor(m); init_GateLETActor(m); + init_GateChemistryActor(m); + init_GateChemistryLongTimeActor(m); init_GateSimulationStatisticsActor(m); init_GatePhaseSpaceActor(m); // init_GateComptonSplittingActor(m); diff --git a/core/opengate_core/opengate_lib/GateChemistryActor.cpp b/core/opengate_core/opengate_lib/GateChemistryActor.cpp new file mode 100644 index 000000000..f16e5ec53 --- /dev/null +++ b/core/opengate_core/opengate_lib/GateChemistryActor.cpp @@ -0,0 +1,244 @@ +/* -------------------------------------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + ------------------------------------ -------------- */ + +#include "GateChemistryActor.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "GateHelpers.h" +#include "GateHelpersDict.h" + +#include "../g4_bindings/chemistryadaptator.h" +#include "GateVActor.h" + +GateChemistryActor::GateChemistryActor(pybind11::dict &user_info) + : GateVActor(user_info, true) { + fActions.insert("NewStage"); + fActions.insert("EndOfRunAction"); + fActions.insert("EndOfEventAction"); + fActions.insert("SteppingAction"); + fActions.insert("EndSimulationAction"); + + _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")); + + G4MoleculeCounter::Instance()->SetVerbose(_moleculeCounterVerbose); + G4MoleculeCounter::Instance()->Use(); + G4MoleculeCounter::Instance()->DontRegister(G4H2O::Definition()); + G4MolecularConfiguration::PrintAll(); + + auto timeStepModel = fIRT; + if (_timeStepModelStr == "IRT") + timeStepModel = fIRT; + else if (_timeStepModelStr == "SBS") + timeStepModel = fSBS; + else if (_timeStepModelStr == "IRT_syn") + timeStepModel = fIRT_syn; + else /* TODO error; detect Python-side? */ + ; + + G4Scheduler::Instance()->SetEndTime(_endTime); + + if (!_reactions.empty()) + setupConstructReactionTableHook(); + else if (!_keepDefaultReactions) + Fatal("Disabling default reactions requires to provide reactions"); + + auto *chemistryList = + ChemistryAdaptator::getChemistryList(); + if (chemistryList != nullptr) + chemistryList->SetTimeStepModel(timeStepModel); +} + +void GateChemistryActor::Initialize(G4HCofThisEvent *hce) { + GateVActor::Initialize(hce); + + G4MoleculeCounter::Instance()->ResetCounter(); + G4DNAChemistryManager::Instance()->ResetCounterWhenRunEnds(false); +} + +void GateChemistryActor::EndSimulationAction() {} + +void GateChemistryActor::EndOfRunAction(G4Run const *) { + for (auto &[molecule, map] : _speciesInfoPerTime) { + for (auto &[time, data] : map) { + data.g /= _nbEvents; // mean value of g + data.sqG /= _nbEvents; // mean value of g²(TODO?) + } + } +} + +void GateChemistryActor::EndOfEventAction(G4Event const *) { + auto *moleculeCounter = G4MoleculeCounter::Instance(); + + if (not G4EventManager::GetEventManager() + ->GetConstCurrentEvent() + ->IsAborted()) { + auto species = moleculeCounter->GetRecordedMolecules(); + if (species && !species->empty()) { + for (auto const *molecule : *species) { + auto &speciesInfo = _speciesInfoPerTime[molecule]; + + for (auto time : _timesToRecord) { + auto nbMol = moleculeCounter->GetNMoleculesAtTime(molecule, time); + + if (nbMol < 0) { + G4cerr << "Invalid molecule count: " << nbMol << " < 0 " << G4endl; + G4Exception("", "N < 0", FatalException, ""); + } + + double gValue = (nbMol / (_edepSum / CLHEP::eV)) * 100.; + + auto &molInfo = speciesInfo[time]; + molInfo.count += nbMol; + molInfo.g += gValue; + molInfo.sqG += gValue * gValue; + } + } + } else + G4cout << "[GateChemistryActor] No molecule recorded, edep = " + << G4BestUnit(_edepSum, "Energy") << G4endl; + } + + ++_nbEvents; + _edepSumRun += _edepSum; + _edepSum = 0.; + moleculeCounter->ResetCounter(); +} + +void GateChemistryActor::SteppingAction(G4Step *step) { + auto edep = step->GetTotalEnergyDeposit(); + if (edep <= 0.) + return; + + edep *= step->GetPreStepPoint()->GetWeight(); + _edepSum += edep; +} + +void GateChemistryActor::NewStage() { + auto *stackManager = G4EventManager::GetEventManager()->GetStackManager(); + if (stackManager != nullptr && stackManager->GetNTotalTrack() == 0) { + G4DNAChemistryManager::Instance()->Run(); + } +} + +void GateChemistryActor::setTimeBinsCount(int n) { + double timeMin = 1 * CLHEP::ps; + double timeMax = G4Scheduler::Instance()->GetEndTime() - 1 * CLHEP::ps; + double timeMinLog = std::log10(timeMin); + double timeMaxLog = std::log10(timeMax); + double timeStepLog = (timeMaxLog - timeMinLog) / (n - 1); + + _timesToRecord.clear(); + for (int i = 0; i < n; ++i) + _timesToRecord.insert(std::pow(10, timeMinLog + i * timeStepLog)); +} + +pybind11::list GateChemistryActor::getTimes() const { + pybind11::list o; + std::for_each(std::begin(_timesToRecord), std::end(_timesToRecord), + [&o](auto const &v) { o.append(v); }); + return o; +} + +pybind11::dict GateChemistryActor::getData() const { + pybind11::dict o; + for (auto const &[molecule, map] : _speciesInfoPerTime) { + pybind11::dict dMolecule; + auto const *name = molecule->GetName().c_str(); + + py::list count; + py::list g; + py::list sqG; + + for (auto const &[time, data] : + map) { // time order is guaranteed by std::map + count.append(data.count); + g.append(data.g); + sqG.append(data.sqG); + } + + if (not o.contains(name)) { + dMolecule["count"] = count; + dMolecule["g"] = g; + dMolecule["sqG"] = sqG; + + o[name] = dMolecule; + } + } + + return o; +} + +GateChemistryActor::ReactionInputs +GateChemistryActor::getReactionInputs(pybind11::dict &user_info, + std::string const &key) { + ReactionInputs reactionInputs; + + auto reactions = DictGetVecList(user_info, key); + for (auto const &reaction : reactions) { + ReactionInput reactionInput; + + reactionInput.reactants = reaction[0].cast>(); + reactionInput.products = reaction[1].cast>(); + reactionInput.fix = reaction[2].cast(); + reactionInput.rate = reaction[3].cast(); + reactionInput.type = reaction[4].cast(); + + reactionInputs.push_back(reactionInput); + } + + 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 new file mode 100644 index 000000000..15dd8e62c --- /dev/null +++ b/core/opengate_core/opengate_lib/GateChemistryActor.h @@ -0,0 +1,80 @@ +/* -------------------------------------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + -------------------------------------------------- */ + +#ifndef GateChemistryActor_h +#define GateChemistryActor_h + +#include +#include +#include +#include + +#include "GateVActor.h" + +class G4EmCalculator; +class G4HCofThisEvent; + +class GateChemistryActor : public GateVActor { +public: + // Constructor + GateChemistryActor(pybind11::dict &user_info); + + void Initialize(G4HCofThisEvent *hce) override; + + void EndSimulationAction() override; + void EndOfRunAction(const G4Run *event) override; + void EndOfEventAction(const G4Event *event) override; + void SteppingAction(G4Step *step) override; + void NewStage() override; + + void setTimeBinsCount(int); + + [[nodiscard]] pybind11::list getTimes() const; + [[nodiscard]] pybind11::dict getData() const; + +public: + struct ReactionInput { + std::vector reactants; + std::vector products; + std::string fix; + double rate; + int type; + }; + + using ReactionInputs = std::vector; + + struct SpeciesInfo { + int count = 0; + double g = 0.; + double sqG = 0.; + }; + + using SpeciesPtr = G4MolecularConfiguration const *; + using InnerSpeciesMap = std::map; + using SpeciesMap = std::map; + +protected: + static ReactionInputs getReactionInputs(pybind11::dict &user_info, + std::string const &key); + void setupConstructReactionTableHook(); + +private: + SpeciesMap _speciesInfoPerTime; + + double _edepSum = 0; + double _edepSumRun = 0; + unsigned _nbEvents = 0; + std::set _timesToRecord; + + int _moleculeCounterVerbose = 0; + std::string _timeStepModelStr = "IRT"; + double _endTime; + bool _keepDefaultReactions = true; + std::vector _reactions; +}; + +#endif diff --git a/core/opengate_core/opengate_lib/GateChemistryLongTimeActor.cpp b/core/opengate_core/opengate_lib/GateChemistryLongTimeActor.cpp new file mode 100644 index 000000000..3cadbe100 --- /dev/null +++ b/core/opengate_core/opengate_lib/GateChemistryLongTimeActor.cpp @@ -0,0 +1,379 @@ +/* -------------------------------------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + ------------------------------------ -------------- */ + +#include "GateChemistryLongTimeActor.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "GateHelpersDict.h" +#include +#include + +#include "../g4_bindings/chemistryadaptator.h" +#include "GateVActor.h" + +GateChemistryLongTimeActor::ChemistryWorld::ChemistryWorld( + GateChemistryLongTimeActor *actor) + : _actor(actor) {} + +void GateChemistryLongTimeActor::ChemistryWorld::ConstructChemistryBoundary() { + auto const &size = _actor->_boundarySize; + std::initializer_list boundingBox{ + // high, low + size[0], -size[0], // x + size[1], -size[1], // y + size[2], -size[2] // z + }; + fpChemistryBoundary = std::make_unique(boundingBox); +} + +void GateChemistryLongTimeActor::ChemistryWorld:: + ConstructChemistryComponents() { + constexpr auto water = 55.3; // from NIST material database + constexpr auto moleLiter = CLHEP::mole * CLHEP::liter; + + constexpr double pKw = 14; // at 25°C pK of water is 14 + auto const &pH = _actor->_pH; + + auto *moleculeTable = G4MoleculeTable::Instance(); + + auto *mH2O = moleculeTable->GetConfiguration("H2O"); + fpChemicalComponent[mH2O] = water / moleLiter; + + auto *mH3Op = moleculeTable->GetConfiguration("H3Op(B)"); + fpChemicalComponent[mH3Op] = std::pow(10, -pH) / moleLiter; // pH = 7 + + auto *mOHm = moleculeTable->GetConfiguration("OHm(B)"); + fpChemicalComponent[mOHm] = std::pow(10, -(pKw - pH)) / moleLiter; // pH = 7 + + auto *mO2 = moleculeTable->GetConfiguration("O2"); + fpChemicalComponent[mO2] = (0. / 100) * 0.0013 / moleLiter; + + // apply scavanger configurations + for (auto const &scavengerConfig : _actor->_scavengerConfigs) { + auto *m = moleculeTable->GetConfiguration(scavengerConfig.species); + auto const &concentration = scavengerConfig.concentration; + auto const &unitWithPrefix = scavengerConfig.unit; + if (not unitWithPrefix.empty()) { + char unit = *unitWithPrefix.rbegin(); + if (unit == 'M') { + double factor = 1; + if (unitWithPrefix.size() > 1) { + auto prefix = unitWithPrefix.substr(0, unitWithPrefix.size() - 1); + if (prefix == "u") + factor = 1e-6; + else if (prefix == "m") + factor = 1e-3; + else // TODO error + ; + } + + fpChemicalComponent[m] = factor * concentration / moleLiter; + } else if (unit == '%') { + constexpr auto o2 = 0.0013; + double concentrationInM = fpChemicalComponent[m] = + (concentration / 100) * o2 / moleLiter; + } + } + } +} + +/* *** */ + +GateChemistryLongTimeActor::GateChemistryLongTimeActor( + pybind11::dict &user_info) + : GateVActor(user_info, true) { + fActions.insert("NewStage"); + fActions.insert("EndOfRunAction"); + fActions.insert("EndOfEventAction"); + fActions.insert("SteppingAction"); + fActions.insert("EndSimulationAction"); + + _timeStepModelStr = DictGetStr(user_info, "timestep_model"); + _endTime = DictGetDouble(user_info, "end_time"); + _reactions = getReactionInputs(user_info, "reactions"); + _moleculeCounterVerbose = DictGetInt(user_info, "molecule_counter_verbose"); + + setTimeBinsCount(DictGetInt(user_info, "time_bins_count")); + + _pH = DictGetDouble(user_info, "pH"); + _scavengerConfigs = getScavengerConfigs(user_info, "scavengers"); + + _doseCutOff = DictGetDouble(user_info, "dose_cutoff"); + + _resetScavengerForEachBeam = + DictGetBool(user_info, "reset_scavenger_for_each_beam"); + + // TODO remove + _boundarySize = DictGetVecDouble(user_info, "boundary_size"); +} + +void GateChemistryLongTimeActor::Initialize(G4HCofThisEvent *hce) { + GateVActor::Initialize(hce); + + _chemistryWorld = std::make_unique(this); + _chemistryWorld->ConstructChemistryComponents(); + + G4Scheduler::Instance()->ResetScavenger(_resetScavengerForEachBeam); + + G4MoleculeCounter::Instance()->SetVerbose(_moleculeCounterVerbose); + G4MoleculeCounter::Instance()->Use(); + G4MoleculeCounter::Instance()->DontRegister(G4H2O::Definition()); + G4MolecularConfiguration::PrintAll(); + + auto timeStepModel = fIRT; + if (_timeStepModelStr == "IRT") + timeStepModel = fIRT; + else if (_timeStepModelStr == "SBS") + timeStepModel = fSBS; + else if (_timeStepModelStr == "IRT_syn") + timeStepModel = fIRT_syn; + else /* TODO error; detect Python-side? */ + ; + + 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); + } + + auto *chemistryList = + ChemistryAdaptator::getChemistryList(); + if (chemistryList != nullptr) + chemistryList->SetTimeStepModel(timeStepModel); + + G4MoleculeCounter::Instance()->ResetCounter(); + G4DNAChemistryManager::Instance()->ResetCounterWhenRunEnds(false); + + G4MoleculeCounter::Instance()->ResetCounter(); +} + +void GateChemistryLongTimeActor::EndSimulationAction() {} + +void GateChemistryLongTimeActor::EndOfRunAction(G4Run const *) { + for (auto &[molecule, map] : _speciesInfoPerTime) { + for (auto &[time, data] : map) { + data.g /= _nbEvents; // mean value of g + data.sqG /= _nbEvents; // mean value of g²(TODO?) + } + } +} + +void GateChemistryLongTimeActor::EndOfEventAction(G4Event const *) { + auto *moleculeCounter = G4MoleculeCounter::Instance(); + + if (not G4EventManager::GetEventManager() + ->GetConstCurrentEvent() + ->IsAborted()) { + auto species = moleculeCounter->GetRecordedMolecules(); + if (species && !species->empty()) { + for (auto const *molecule : *species) { + auto &speciesInfo = _speciesInfoPerTime[molecule]; + + for (auto time : _timesToRecord) { + auto nbMol = moleculeCounter->GetNMoleculesAtTime(molecule, time); + + if (nbMol < 0) { + G4cerr << "Invalid molecule count: " << nbMol << " < 0 " << G4endl; + G4Exception("", "N < 0", FatalException, ""); + } + + double gValue = (nbMol / (_edepSum / CLHEP::eV)) * 100.; + + auto &molInfo = speciesInfo[time]; + molInfo.count += nbMol; + molInfo.g += gValue; + molInfo.sqG += gValue * gValue; + } + } + } else + G4cout << "[GateChemistryLongTimeActor] No molecule recorded, edep = " + << G4BestUnit(_edepSum, "Energy") << G4endl; + } + + ++_nbEvents; + _edepSumRun += _edepSum; + _edepSum = 0.; + moleculeCounter->ResetCounter(); +} + +void GateChemistryLongTimeActor::SteppingAction(G4Step *step) { + auto edep = step->GetTotalEnergyDeposit(); + if (edep <= 0.) + return; + + edep *= step->GetPreStepPoint()->GetWeight(); + _edepSum += edep; + + auto *track = step->GetTrack(); + + if (track->GetParentID() == 0 && track->GetCurrentStepNumber() == 1) { + constexpr auto density = .001; + auto volume = _chemistryWorld->GetChemistryBoundary()->Volume(); + double dose = (_edepSum / CLHEP::eV) / (density * volume * 6.242e+18); + if (dose > _doseCutOff) { + auto *eventManager = G4EventManager::GetEventManager(); + auto *primary = eventManager->GetConstCurrentEvent() + ->GetPrimaryVertex() + ->GetPrimary(); + auto const &name = primary->GetParticleDefinition()->GetParticleName(); + auto energy = primary->GetKineticEnergy(); + G4cout << "[GateChemistryLongTimeActor] stop beam line '" << name << "' (" + << energy << " MeV) at dose: " << dose << " Gy" << G4endl; + + track->SetTrackStatus(fStopAndKill); + auto const *secondaries = track->GetStep()->GetSecondaryInCurrentStep(); + for (auto const *secondary : *secondaries) { + if (secondary != nullptr) { + // FIXME + // from UHDR example + // must find how to get non-const access to secondaries + auto *secondaryMut = const_cast(secondary); + secondaryMut->SetTrackStatus(fStopAndKill); + } + } + eventManager->GetStackManager()->ClearUrgentStack(); + } + } +} + +void GateChemistryLongTimeActor::NewStage() { + auto *stackManager = G4EventManager::GetEventManager()->GetStackManager(); + if (stackManager != nullptr && stackManager->GetNTotalTrack() == 0) { + G4DNAChemistryManager::Instance()->Run(); + } +} + +void GateChemistryLongTimeActor::setTimeBinsCount(int n) { + double timeMin = 1 * CLHEP::ps; + double timeMax = G4Scheduler::Instance()->GetEndTime() - 1 * CLHEP::ps; + double timeMinLog = std::log10(timeMin); + double timeMaxLog = std::log10(timeMax); + double timeStepLog = (timeMaxLog - timeMinLog) / (n - 1); + + _timesToRecord.clear(); + for (int i = 0; i < n; ++i) + _timesToRecord.insert(std::pow(10, timeMinLog + i * timeStepLog)); +} + +pybind11::list GateChemistryLongTimeActor::getTimes() const { + pybind11::list o; + std::for_each(std::begin(_timesToRecord), std::end(_timesToRecord), + [&o](auto const &v) { o.append(v); }); + return o; +} + +pybind11::dict GateChemistryLongTimeActor::getData() const { + pybind11::dict o; + for (auto const &[molecule, map] : _speciesInfoPerTime) { + pybind11::dict dMolecule; + auto const *name = molecule->GetName().c_str(); + + py::list count; + py::list g; + py::list sqG; + + for (auto const &[time, data] : + map) { // time order is guaranteed by std::map + count.append(data.count); + g.append(data.g); + sqG.append(data.sqG); + } + + if (not o.contains(name)) { + dMolecule["count"] = count; + dMolecule["g"] = g; + dMolecule["sqG"] = sqG; + + o[name] = dMolecule; + } + } + + return o; +} + +GateChemistryLongTimeActor::ScavengerConfigs +GateChemistryLongTimeActor::getScavengerConfigs(pybind11::dict &user_info, + std::string const &key) { + ScavengerConfigs scavengerConfigs; + + auto configs = DictGetVecList(user_info, key); + for (auto const &config : configs) { + ScavengerConfig scavengerConfig; + + scavengerConfig.species = config[0].cast(); + scavengerConfig.concentration = config[1].cast(); + scavengerConfig.unit = config[2].cast(); + } + + return scavengerConfigs; +} + +GateChemistryLongTimeActor::ReactionInputs +GateChemistryLongTimeActor::getReactionInputs(pybind11::dict &user_info, + std::string const &key) { + ReactionInputs reactionInputs; + + auto reactions = DictGetVecList(user_info, key); + for (auto const &reaction : reactions) { + ReactionInput reactionInput; + + reactionInput.reactants = reaction[0].cast>(); + reactionInput.products = reaction[1].cast>(); + reactionInput.fix = reaction[2].cast(); + reactionInput.rate = reaction[3].cast(); + reactionInput.type = reaction[4].cast(); + + reactionInputs.push_back(reactionInput); + } + + return reactionInputs; +} diff --git a/core/opengate_core/opengate_lib/GateChemistryLongTimeActor.h b/core/opengate_core/opengate_lib/GateChemistryLongTimeActor.h new file mode 100644 index 000000000..8f79f3a71 --- /dev/null +++ b/core/opengate_core/opengate_lib/GateChemistryLongTimeActor.h @@ -0,0 +1,114 @@ +/* -------------------------------------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + -------------------------------------------------- */ + +#ifndef GateChemistryLongTimeActor_h +#define GateChemistryLongTimeActor_h + +#include +#include +#include +#include +#include +#include + +#include "GateVActor.h" + +class G4EmCalculator; +class G4HCofThisEvent; + +class GateChemistryLongTimeActor : public GateVActor { +public: + GateChemistryLongTimeActor(pybind11::dict &user_info); + + void Initialize(G4HCofThisEvent *hce) override; + + void EndSimulationAction() override; + void EndOfRunAction(const G4Run *event) override; + void EndOfEventAction(const G4Event *event) override; + void SteppingAction(G4Step *step) override; + void NewStage() override; + + void setTimeBinsCount(int); + + [[nodiscard]] pybind11::list getTimes() const; + [[nodiscard]] pybind11::dict getData() const; + +public: + struct ScavengerConfig { + std::string species; + double concentration; + std::string unit; + }; + + using ScavengerConfigs = std::vector; + + struct ReactionInput { + std::vector reactants; + std::vector products; + std::string fix; + double rate; + int type; + }; + + using ReactionInputs = std::vector; + + struct SpeciesInfo { + int count = 0; + double g = 0.; + double sqG = 0.; + }; + + using SpeciesPtr = G4MolecularConfiguration const *; + using InnerSpeciesMap = std::map; + using SpeciesMap = std::map; + +private: + class ChemistryWorld : public G4VChemistryWorld { + public: + ChemistryWorld(GateChemistryLongTimeActor *actor); + + void ConstructChemistryBoundary() override; + void ConstructChemistryComponents() override; + + private: + GateChemistryLongTimeActor *_actor; + }; + +protected: + static ScavengerConfigs getScavengerConfigs(pybind11::dict &user_info, + std::string const &key); + + static ReactionInputs getReactionInputs(pybind11::dict &user_info, + std::string const &key); + +private: + SpeciesMap _speciesInfoPerTime; + + double _edepSum = 0; + double _edepSumRun = 0; + unsigned _nbEvents = 0; + std::set _timesToRecord; + + int _moleculeCounterVerbose = 0; + std::string _timeStepModelStr = "SBS"; + double _endTime; + std::vector _reactions; + + std::unique_ptr _chemistryWorld; + + double _pH = 7; + ScavengerConfigs _scavengerConfigs; + + double _doseCutOff = 0; + + bool _resetScavengerForEachBeam = false; + + // TODO remove + std::vector _boundarySize; +}; + +#endif diff --git a/core/opengate_core/opengate_lib/GateStackingAction.cpp b/core/opengate_core/opengate_lib/GateStackingAction.cpp new file mode 100644 index 000000000..e0e104871 --- /dev/null +++ b/core/opengate_core/opengate_lib/GateStackingAction.cpp @@ -0,0 +1,22 @@ +/* -------------------------------------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + -------------------------------------------------- */ + +#include "GateStackingAction.h" + +void GateStackingAction::RegisterActor(GateVActor *actor) { + auto actions = actor->fActions; + auto beg = std::find(std::begin(actions), std::end(actions), "NewStage"); + if (beg != actions.end()) { + fNewStageActors.push_back(actor); + } +} + +void GateStackingAction::NewStage() { + for (auto *actor : fNewStageActors) { + actor->NewStage(); + } +} diff --git a/core/opengate_core/opengate_lib/GateStackingAction.h b/core/opengate_core/opengate_lib/GateStackingAction.h new file mode 100644 index 000000000..44ada6bbf --- /dev/null +++ b/core/opengate_core/opengate_lib/GateStackingAction.h @@ -0,0 +1,31 @@ +/* -------------------------------------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + -------------------------------------------------- */ + +#ifndef GateStackingAction_h +#define GateStackingAction_h + +#include "G4UserStackingAction.hh" +#include "GateVActor.h" + +class GateStackingAction : public G4UserStackingAction { + +public: + GateStackingAction() = default; + + ~GateStackingAction() override = default; + + void RegisterActor(GateVActor *actor); + + void NewStage() override; + + auto *stackManager() { return G4UserStackingAction::stackManager; } + +protected: + std::vector fNewStageActors; +}; + +#endif // GateTrackingAction_h diff --git a/core/opengate_core/opengate_lib/GateVActor.h b/core/opengate_core/opengate_lib/GateVActor.h index 09731b83c..a8bdded28 100644 --- a/core/opengate_core/opengate_lib/GateVActor.h +++ b/core/opengate_core/opengate_lib/GateVActor.h @@ -96,6 +96,9 @@ class GateVActor : public G4VPrimitiveScorer { // Called every FillHits, should be overloaded virtual void SteppingAction(G4Step *) {} + // TODO + virtual void NewStage() {} + // List of actions (set to trigger some actions) // Can be set either on cpp or py side std::set fActions; diff --git a/core/opengate_core/opengate_lib/pyGateChemistryActor.cpp b/core/opengate_core/opengate_lib/pyGateChemistryActor.cpp new file mode 100644 index 000000000..f5d7c1130 --- /dev/null +++ b/core/opengate_core/opengate_lib/pyGateChemistryActor.cpp @@ -0,0 +1,22 @@ +/* -------------------------------------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + -------------------------------------------------- */ + +#include +#include + +namespace py = pybind11; + +#include "GateChemistryActor.h" + +void init_GateChemistryActor(py::module &m) { + py::class_, GateVActor>( + m, "GateChemistryActor") + .def(py::init()) + .def("get_times", &GateChemistryActor::getTimes) + .def("get_data", &GateChemistryActor::getData); +} diff --git a/core/opengate_core/opengate_lib/pyGateChemistryLongTimeActor.cpp b/core/opengate_core/opengate_lib/pyGateChemistryLongTimeActor.cpp new file mode 100644 index 000000000..2a37901ef --- /dev/null +++ b/core/opengate_core/opengate_lib/pyGateChemistryLongTimeActor.cpp @@ -0,0 +1,22 @@ +/* -------------------------------------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + -------------------------------------------------- */ + +#include +#include + +namespace py = pybind11; + +#include "GateChemistryLongTimeActor.h" + +void init_GateChemistryLongTimeActor(py::module &m) { + py::class_, + GateVActor>(m, "GateChemistryLongTimeActor") + .def(py::init()) + .def("get_times", &GateChemistryLongTimeActor::getTimes) + .def("get_data", &GateChemistryLongTimeActor::getData); +} diff --git a/core/opengate_core/opengate_lib/pyGateStackingAction.cpp b/core/opengate_core/opengate_lib/pyGateStackingAction.cpp new file mode 100644 index 000000000..34f396871 --- /dev/null +++ b/core/opengate_core/opengate_lib/pyGateStackingAction.cpp @@ -0,0 +1,24 @@ +/* -------------------------------------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + -------------------------------------------------- */ + +#include +#include + +namespace py = pybind11; + +#include "G4UserStackingAction.hh" +#include "GateSourceManager.h" +#include "GateStackingAction.h" + +void init_GateStackingAction(py::module &m) { + + py::class_>( + m, "GateStackingAction") + .def(py::init<>()) + .def("RegisterActor", &GateStackingAction::RegisterActor); +} diff --git a/docs/source/user_guide_2_3_physics.md b/docs/source/user_guide_2_3_physics.md index b14ed19ad..23739d37a 100644 --- a/docs/source/user_guide_2_3_physics.md +++ b/docs/source/user_guide_2_3_physics.md @@ -72,6 +72,14 @@ G4EmLivermorePhysics G4EmLivermorePolarizedPhysics G4EmPenelopePhysics G4EmDNAPhysics +G4EmDNAPhysics_option1 +G4EmDNAPhysics_option2 +G4EmDNAPhysics_option3 +G4EmDNAPhysics_option4 +G4EmDNAPhysics_option5 +G4EmDNAPhysics_option6 +G4EmDNAPhysics_option7 +G4EmDNAPhysics_option8 G4OpticalPhysics ``` @@ -92,6 +100,13 @@ Under the hood, this will add two processed to the Geant4 list of processes, G4D - - +**WARNING** The chemistry modules, if needed, must be added explicitly. This is done with: + +``` +sim.enable_dnachemistry(???) +``` + + ### Optical Physics Processes #### G4OpticalPhysics physics list diff --git a/opengate/actors/actorbuilders.py b/opengate/actors/actorbuilders.py index 8b6c2f697..41656e79c 100644 --- a/opengate/actors/actorbuilders.py +++ b/opengate/actors/actorbuilders.py @@ -1,5 +1,6 @@ from .arfactors import ARFActor, ARFTrainingDatasetActor from .doseactors import DoseActor, LETActor, FluenceActor +from .chemistryactors import ChemistryActor, ChemistryLongTimeActor from .digitizers import ( DigitizerAdderActor, DigitizerReadoutActor, @@ -29,6 +30,8 @@ DoseActor, FluenceActor, LETActor, + ChemistryActor, + ChemistryLongTimeActor, SourceInfoActor, PhaseSpaceActor, DigitizerHitsCollectionActor, diff --git a/opengate/actors/chemistryactors.py b/opengate/actors/chemistryactors.py new file mode 100644 index 000000000..d9a46a583 --- /dev/null +++ b/opengate/actors/chemistryactors.py @@ -0,0 +1,96 @@ +import opengate_core as g4 +from .base import ActorBase +from ..utility import g4_units + + +class ChemistryActor(g4.GateChemistryActor, ActorBase): + """ """ + + type_name = "ChemistryActor" + + def set_default_user_info(user_info): + ActorBase.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 + + def __init__(self, user_info): + ActorBase.__init__(self, user_info) + g4.GateChemistryActor.__init__(self, user_info.__dict__) + + def __str__(self): + u = self.user_info + s = f'ChemistryActor "{u.name}": dim={u.size} spacing={u.spacing} {u.output} tr={u.translation}' + return s + + def __getstate__(self): + # superclass getstate + ActorBase.__getstate__(self) + return self.__dict__ + + def initialize(self, volume_engine=None): + super().initialize(volume_engine) + + def StartSimulationAction(self): + pass + + def EndSimulationAction(self): + pass + + +class ChemistryLongTimeActor(g4.GateChemistryLongTimeActor, ActorBase): + """ """ + + type_name = "ChemistryLongTimeActor" + + def set_default_user_info(user_info): + ActorBase.set_default_user_info(user_info) + + user_info.output = "output.root" + user_info.timestep_model = "IRT" + user_info.reactions = [] + user_info.end_time = 1 * g4_units.s + user_info.time_bins_count = 10 + user_info.molecule_counter_verbose = 0 + + user_info.pH = 7 + user_info.scavengers = [] + + user_info.dose_cutoff = 0 + + user_info.reset_scavenger_for_each_beam = False + + # TODO remove + user_info.boundary_size = [ + 0.8 * g4_units.um, + 0.8 * g4_units.um, + 0.8 * g4_units.um, + ] + + def __init__(self, user_info): + ActorBase.__init__(self, user_info) + g4.GateChemistryLongTimeActor.__init__(self, user_info.__dict__) + + def __str__(self): + u = self.user_info + s = f'ChemistryLongTimeActor "{u.name}": dim={u.size} spacing={u.spacing} {u.output} tr={u.translation}' + return s + + def __getstate__(self): + # superclass getstate + ActorBase.__getstate__(self) + return self.__dict__ + + def initialize(self, volume_engine=None): + super().initialize(volume_engine) + + def StartSimulationAction(self): + pass + + def EndSimulationAction(self): + pass diff --git a/opengate/engines.py b/opengate/engines.py index 9d9255adc..a34c6cacd 100644 --- a/opengate/engines.py +++ b/opengate/engines.py @@ -437,6 +437,7 @@ def __init__(self, simulation_engine): self.g4_RunAction = [] self.g4_EventAction = [] self.g4_TrackingAction = [] + self.g4_StackingAction = [] def close(self): if self.verbose_close: @@ -450,6 +451,7 @@ def release_g4_references(self): self.g4_RunAction = None self.g4_EventAction = None self.g4_TrackingAction = None + self.g4_StackingAction = None def BuildForMaster(self): # This function is call only in MT mode, for the master thread @@ -492,6 +494,11 @@ def Build(self): self.SetUserAction(ta) self.g4_TrackingAction.append(ta) + # add chemistry stage + sa = g4.GateStackingAction() + self.SetUserAction(sa) + self.g4_StackingAction.append(sa) + class ActorEngine(EngineBase): """ @@ -566,6 +573,9 @@ def register_all_actions(self, actor): # Track for ta in self.simulation_engine_wr().action_engine.g4_TrackingAction: ta.RegisterActor(actor) + # Stacking + for sa in self.simulation_engine_wr().action_engine.g4_StackingAction: + sa.RegisterActor(actor) # initialization actor.ActorInitialize() diff --git a/opengate/managers.py b/opengate/managers.py index d2f8620b4..2f3b44b34 100644 --- a/opengate/managers.py +++ b/opengate/managers.py @@ -343,7 +343,19 @@ class PhysicsListManager(GateObject): "G4EmLivermorePolarizedPhysics", "G4EmPenelopePhysics", "G4EmDNAPhysics", + "G4EmDNAPhysics_option1", + "G4EmDNAPhysics_option2", + "G4EmDNAPhysics_option3", + "G4EmDNAPhysics_option4", + "G4EmDNAPhysics_option5", + "G4EmDNAPhysics_option6", + "G4EmDNAPhysics_option7", + "G4EmDNAPhysics_option8", "G4OpticalPhysics", + "G4EmDNAChemistry", + "G4EmDNAChemistry_option1", + "G4EmDNAChemistry_option2", + "G4EmDNAChemistry_option3", "G4GenericBiasingPhysics", ] @@ -358,6 +370,22 @@ class PhysicsListManager(GateObject): g4.G4GenericBiasingPhysics ) + special_chemistry_constructor_classes = {} + special_chemistry_constructor_classes["G4EmDNAChemistry_option3"] = ( + g4.G4EmDNAChemistry_option3 + ) + + special_physics_constructor_classes["G4EmDNAChemistry"] = g4.G4EmDNAChemistry + special_physics_constructor_classes["G4EmDNAChemistry_option1"] = ( + g4.G4EmDNAChemistry_option1 + ) + special_physics_constructor_classes["G4EmDNAChemistry_option2"] = ( + g4.G4EmDNAChemistry_option2 + ) + special_physics_constructor_classes["G4EmDNAChemistry_option3"] = ( + g4.G4EmDNAChemistry_option3 + ) + def __init__(self, physics_manager, *args, **kwargs): super().__init__(*args, **kwargs) self.physics_manager = physics_manager @@ -550,6 +578,17 @@ class PhysicsManager(GateObject): "doc": "Define the process to bias (if wanted) on the different particle types." }, ), + "special_chemistry_constructors": ( + Box( + [ + (spc, False) + for spc in PhysicsListManager.special_chemistry_constructor_classes + ] + ), + { + "doc": "Special chemistry constructors to be added to the physics list, e.g. G4EmDNAChemistry, G4EmDNAChemistry_option1, …, G4EmDNAChemistry_option4. " + }, + ), } def __init__(self, simulation, *args, **kwargs):