diff --git a/core/opengate_core/opengate_lib/GateChemistryActor.cpp b/core/opengate_core/opengate_lib/GateChemistryActor.cpp index 64fd7584e..3cdd9654c 100644 --- a/core/opengate_core/opengate_lib/GateChemistryActor.cpp +++ b/core/opengate_core/opengate_lib/GateChemistryActor.cpp @@ -8,15 +8,19 @@ #include "GateChemistryActor.h" #include -#include #include -#include #include +#include #include -#include -#include +#include +#include +#include #include +#include #include +#include +#include +#include #include #include @@ -25,204 +29,193 @@ #include "../g4_bindings/chemistryadaptator.h" #include "GateVActor.h" -GateChemistryActor::GateChemistryActor(pybind11::dict &user_info) - : GateVActor(user_info, true) { +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"); - G4MoleculeCounter::Instance()->Use(); - G4MolecularConfiguration::PrintAll(); - - auto timeStepModelStr = DictGetStr(user_info, "timestep_model"); - 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? */ - ; - - // TODO user defined - G4Scheduler::Instance()->SetEndTime(DictGetDouble(user_info, "end_time")); - setTimeBinsCount(DictGetInt(user_info, "time_bins_count")); - - { - std::vector reactions = - getReactionInputs(user_info, "reactions"); - - auto constructReactionTable = - [reactions = - std::move(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); + _timeStepModelStr = DictGetStr(user_info, "timestep_model"); + _endTime = DictGetDouble(user_info, "end_time"); + _reactions = getReactionInputs(user_info, "reactions"); + + setTimeBinsCount(DictGetInt(user_info, "time_bins_count")); } -void GateChemistryActor::Initialize(G4HCofThisEvent *hce) { - GateVActor::Initialize(hce); +void GateChemistryActor::Initialize(G4HCofThisEvent* hce) { + GateVActor::Initialize(hce); + + 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(); - G4MoleculeCounter::Instance()->ResetCounter(); + 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 GateChemistryActor::EndSimulationAction() {} +void GateChemistryActor::EndSimulationAction() { +} -void GateChemistryActor::EndOfRunAction(G4Run const *) { - G4cout << "************ EndOfRun *******" << G4endl; - for (auto [species, map] : _speciesInfoPerTime) { - for (auto [time, data] : map) { - G4cout << species->GetName() << " @ " << time << ": " << data.number - << ", " << data.g << ", " << data.sqG << G4endl; - } - } +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 *molecularCounter = G4MoleculeCounter::Instance(); - - if (not G4EventManager::GetEventManager() - ->GetConstCurrentEvent() - ->IsAborted()) { - auto species = molecularCounter->GetRecordedMolecules(); - if (species && !species->empty()) { - for (auto const *molecule : *species) { - auto &speciesInfo = _speciesInfoPerTime[molecule]; - - for (auto time : _timesToRecord) { - auto nbMol = molecularCounter->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.number += nbMol; - molInfo.g += gValue; - molInfo.sqG += gValue * gValue; - } - } - } else - G4cout << "************* No molecule recorded, edep = " - << G4BestUnit(_edepSum, "Energy") << G4endl; - } - - ++_nbEvents; +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.; - molecularCounter->ResetCounter(); + moleculeCounter->ResetCounter(); } -void GateChemistryActor::SteppingAction(G4Step *step) { - auto edep = step->GetTotalEnergyDeposit(); - if (edep <= 0.) - return; +void GateChemistryActor::SteppingAction(G4Step* step) { + auto edep = step->GetTotalEnergyDeposit(); + if(edep <= 0.) return; - edep *= step->GetPreStepPoint()->GetWeight(); - _edepSum += edep; + edep *= step->GetPreStepPoint()->GetWeight(); + _edepSum += edep; } void GateChemistryActor::NewStage() { - auto *stackManager = G4EventManager::GetEventManager()->GetStackManager(); - if (stackManager != nullptr && stackManager->GetNTotalTrack() == 0) { - G4DNAChemistryManager::Instance()->Run(); - } + 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)); + 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; + 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 &[species, map] : _speciesInfoPerTime) { - pybind11::dict dMolecule; - - 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.number); - g.append(data.g); - sqG.append(data.sqG); - } - - dMolecule["count"] = count; - dMolecule["g"] = g; - dMolecule["sqG"] = sqG; - - o[species->GetName()] = dMolecule; - } - - return 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; +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; + 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(); + 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); - } + reactionInputs.push_back(reactionInput); + } - return reactionInputs; + return reactionInputs; } diff --git a/core/opengate_core/opengate_lib/GateChemistryActor.h b/core/opengate_core/opengate_lib/GateChemistryActor.h index 2a29101a0..cc35e20a4 100644 --- a/core/opengate_core/opengate_lib/GateChemistryActor.h +++ b/core/opengate_core/opengate_lib/GateChemistryActor.h @@ -48,7 +48,7 @@ class GateChemistryActor : public GateVActor { using ReactionInputs = std::vector; struct SpeciesInfo { - int number = 0; + int count = 0; double g = 0.; double sqG = 0.; }; @@ -65,8 +65,14 @@ class GateChemistryActor : public GateVActor { SpeciesMap _speciesInfoPerTime; double _edepSum = 0; + double _edepSumRun = 0; unsigned _nbEvents = 0; std::set _timesToRecord; + + int _moleculeCounterVerbose = 0; + std::string _timeStepModelStr = "IRT"; + double _endTime; + std::vector _reactions; }; #endif