diff --git a/.gitignore b/.gitignore index f285a8920..c3aab4dc6 100644 --- a/.gitignore +++ b/.gitignore @@ -76,6 +76,8 @@ data/ # --- Tests opengate/tests/output_dashboard/ +opengate/tests/src/external/pytomography/output +opengate/tests/src/external/pytomography/data opengate/tests/log .clang-format a.txt diff --git a/core/opengate_core/opengate_core.cpp b/core/opengate_core/opengate_core.cpp index ab8e2c503..669857fa7 100644 --- a/core/opengate_core/opengate_core.cpp +++ b/core/opengate_core/opengate_core.cpp @@ -403,6 +403,8 @@ void init_GateAttributeComparisonFilter(py::module &); // Gate actors void init_GateDoseActor(py::module &m); +void init_GateDebugActor(py::module &m); + void init_GateVoxelizedPromptGammaTLEActor(py::module &m); void init_GateVoxelizedPromptGammaAnalogActor(py::module &m); @@ -530,6 +532,8 @@ void init_GateSourceManager(py::module &); void init_GateGenericSource(py::module &); +void init_GateDebugSource(py::module &); + void init_GateTreatmentPlanPBSource(py::module &); void init_GateTemplateSource(py::module &); @@ -781,6 +785,7 @@ PYBIND11_MODULE(opengate_core, m) { init_GateLastVertexSource(m); init_GateSourceManager(m); init_GateGenericSource(m); + init_GateDebugSource(m); init_GateTreatmentPlanPBSource(m); init_GateTemplateSource(m); init_GatePencilBeamSource(m); @@ -807,6 +812,7 @@ PYBIND11_MODULE(opengate_core, m) { init_GateProcessDefinedStepInVolumeAttribute(m); init_GateDoseActor(m); + init_GateDebugActor(m); init_GateTLEDoseActor(m); init_GateVoxelizedPromptGammaTLEActor(m); init_GateVoxelizedPromptGammaAnalogActor(m); diff --git a/core/opengate_core/opengate_lib/GateDebugActor.cpp b/core/opengate_core/opengate_lib/GateDebugActor.cpp new file mode 100644 index 000000000..1cf5c2347 --- /dev/null +++ b/core/opengate_core/opengate_lib/GateDebugActor.cpp @@ -0,0 +1,64 @@ +/* -------------------------------------------------- + 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 "GateDebugActor.h" +#include "G4RunManager.hh" +#include "GateHelpers.h" + +GateDebugActor::GateDebugActor(py::dict &user_info) + : GateVActor(user_info, true) { + DDD("(cpp) DebugActor constructor") +} + +void GateDebugActor::InitializeUserInfo(py::dict &user_info) { + DDD("(cpp) DebugActor::InitializeUserInfo") + // IMPORTANT: call the base class method + GateVActor::InitializeUserInfo(user_info); +} + +void GateDebugActor::InitializeCpp() { + DDD("(cpp) DebugActor::InitializeCpp") + GateVActor::InitializeCpp(); +} + +void GateDebugActor::BeginOfRunActionMasterThread(int run_id) { + DDD("(cpp) GateDebugActor::BeginOfRunActionMasterThread") +} + +void GateDebugActor::BeginOfRunAction(const G4Run *run) { + DDD("(cpp) GateDebugActor::BeginOfRunAction", " ", run->GetRunID()); +} + +void GateDebugActor::BeginOfEventAction(const G4Event *event) { + DDD("(cpp) GateDebugActor::BeginOfEventAction", event->GetEventID()); +} + +void GateDebugActor::PreUserTrackingAction(const G4Track *track) { + DDD("(cpp) GateDebugActor::PreUserTrackingAction", " ", track->GetTrackID()); +} + +void GateDebugActor::PostUserTrackingAction(const G4Track *track) { + DDD("(cpp) GateDebugActor::PostUserTrackingAction", " ", track->GetTrackID()); +} + +void GateDebugActor::SteppingAction(G4Step *step) { + DDD("(cpp) GateDebugActor::SteppingAction", " ", + step->GetTrack()->GetCurrentStepNumber()); +} + +void GateDebugActor::EndOfEventAction(const G4Event *event) { + DDD("(cpp) GateDebugActor::EndOfEventAction", " ", event->GetEventID()); +} + +int GateDebugActor::EndOfRunActionMasterThread(int run_id) { + DDD("(cpp) GateDebugActor::EndOfRunActionMasterThread") + return 0; +} + +void GateDebugActor::EndOfRunAction(const G4Run *run) { + DDD("(cpp) GateDebugActor::EndOfRunAction", " ", run->GetRunID()); +} diff --git a/core/opengate_core/opengate_lib/GateDebugActor.h b/core/opengate_core/opengate_lib/GateDebugActor.h new file mode 100644 index 000000000..f15dc6831 --- /dev/null +++ b/core/opengate_core/opengate_lib/GateDebugActor.h @@ -0,0 +1,43 @@ +/* -------------------------------------------------- + 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 GateDebugActor_h +#define GateDebugActor_h + +#include "GateVActor.h" + +namespace py = pybind11; + +class GateDebugActor : public GateVActor { + +public: + explicit GateDebugActor(py::dict &user_info); + + void InitializeUserInfo(py::dict &user_info) override; + + void InitializeCpp() override; + + void BeginOfRunAction(const G4Run *run) override; + + void BeginOfRunActionMasterThread(int run_id) override; + + void BeginOfEventAction(const G4Event *event) override; + + void PreUserTrackingAction(const G4Track *track) override; + + void PostUserTrackingAction(const G4Track *track) override; + + void SteppingAction(G4Step *) override; + + void EndOfEventAction(const G4Event *event) override; + + void EndOfRunAction(const G4Run *run) override; + + int EndOfRunActionMasterThread(int run_id) override; +}; + +#endif // GateDebugActor_h diff --git a/core/opengate_core/opengate_lib/GateDebugSource.cpp b/core/opengate_core/opengate_lib/GateDebugSource.cpp new file mode 100644 index 000000000..668b42c6b --- /dev/null +++ b/core/opengate_core/opengate_lib/GateDebugSource.cpp @@ -0,0 +1,64 @@ +/* -------------------------------------------------- + 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 "GateDebugSource.h" +#include "GateHelpers.h" +#include "GateHelpersDict.h" +#include "fmt/core.h" +#include +#include +#include + +GateDebugSource::GateDebugSource() : GateVSource() { + DDD("GateDebugSource constructor"); + fDebugValue = 0.0; +} + +GateDebugSource::~GateDebugSource() { DDD("GateDebugSource destructor"); } + +void GateDebugSource::CleanWorkerThread() { + DDD("GateDebugSource::CleanWorkerThread"); +} + +void GateDebugSource::InitializeUserInfo(py::dict &user_info) { + DDD("GateDebugSource::InitializeUserInfo"); + GateVSource::InitializeUserInfo(user_info); + fDebugValue = DictGetDouble(user_info, "debug_value"); +} + +void GateDebugSource::UpdateActivity(const double time) { + DDD("GateDebugSource::UpdateActivity"); + GateVSource::UpdateActivity(time); +} + +double GateDebugSource::PrepareNextTime(const double current_simulation_time, + double NumberOfGeneratedEvents) { + DDD("GateDebugSource::PrepareNextTime", " ", + G4BestUnit(current_simulation_time, "Time"), " ", + NumberOfGeneratedEvents); + return GateVSource::PrepareNextTime(current_simulation_time, + NumberOfGeneratedEvents); +} + +void GateDebugSource::PrepareNextRun() { + DDD("GateDebugSource::PrepareNextRun"); + GateVSource::PrepareNextRun(); +} + +void GateDebugSource::GeneratePrimaries(G4Event *event, + const double current_simulation_time) { + DDD("GateDebugSource::GeneratePrimaries", " ", event->GetEventID(), " ", + fMaxN); + fNumberOfGeneratedEvents++; + fDebugValue += 1; + DDD("debug value = ", fDebugValue); +} + +double GateDebugSource::GetDebugValue() { + DDD("GateDebugSource::GetDebugValue ", fDebugValue) + return fDebugValue; +} diff --git a/core/opengate_core/opengate_lib/GateDebugSource.h b/core/opengate_core/opengate_lib/GateDebugSource.h new file mode 100644 index 000000000..2f46ffd6d --- /dev/null +++ b/core/opengate_core/opengate_lib/GateDebugSource.h @@ -0,0 +1,42 @@ +/* -------------------------------------------------- + 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 GateDebugSource_h +#define GateDebugSource_h + +#include "GateVSource.h" +#include +#include + +namespace py = pybind11; + +class GateDebugSource : public GateVSource { + +public: + GateDebugSource(); + ~GateDebugSource() override; + + void CleanWorkerThread() override; + + void InitializeUserInfo(py::dict &user_info) override; + + void UpdateActivity(const double time) override; + + double PrepareNextTime(double current_simulation_time, + double NumberOfGeneratedEvents) override; + + void PrepareNextRun() override; + + void GeneratePrimaries(G4Event *event, + double current_simulation_time) override; + + double GetDebugValue(); + + double fDebugValue = 0.0; +}; + +#endif // GateDebugSource_h diff --git a/core/opengate_core/opengate_lib/GateGANPairSource.cpp b/core/opengate_core/opengate_lib/GateGANPairSource.cpp index 3ffca6e7b..e9811bfbb 100644 --- a/core/opengate_core/opengate_lib/GateGANPairSource.cpp +++ b/core/opengate_core/opengate_lib/GateGANPairSource.cpp @@ -14,8 +14,7 @@ GateGANPairSource::~GateGANPairSource() = default; void GateGANPairSource::InitializeUserInfo(py::dict &user_info) { GateGANSource::InitializeUserInfo(user_info); - auto &l = GetThreadLocalDataGenericSource(); - if (l.fAAManager->IsEnabled()) { + if (fAAManager->IsEnabled()) { std::ostringstream oss; oss << "Error, cannot use AngularAcceptance with GAN pairs (yet), for the " "source '" @@ -71,8 +70,7 @@ void GateGANPairSource::GeneratePrimaries(G4Event *event, fCurrentIndex++; // update the number of generated event - auto &l = fThreadLocalData.Get(); - l.fNumberOfGeneratedEvents++; + fNumberOfGeneratedEvents++; } void GateGANPairSource::GeneratePrimariesPair(G4Event *event, @@ -89,12 +87,11 @@ void GateGANPairSource::GeneratePrimariesPair(G4Event *event, fDirectionZ2[fCurrentIndex]); // move position according to mother volume - auto &l = fThreadLocalData.Get(); - position = l.fGlobalRotation * position + l.fGlobalTranslation; + position = fGlobalRotation * position + fGlobalTranslation; // normalize (needed) direction = direction / direction.mag(); // move according to mother volume - direction = l.fGlobalRotation * direction; + direction = fGlobalRotation * direction; // energy of the second particle double energy = fEnergy2[fCurrentIndex]; @@ -105,22 +102,20 @@ void GateGANPairSource::GeneratePrimariesPair(G4Event *event, if (!accept_energy) { energy = 0; // at least one of the two vertices has been skipped with zeroE - auto &ll = GetThreadLocalDataGenericSource(); - ll.fCurrentZeroEvents = 1; + fCurrentZeroEvents = 1; } - auto &ll = fThreadLocalDataGenericSource.Get(); if (fTime_is_set_by_GAN && accept_energy) { // time - double time = ll.fEffectiveEventTime; + double time = fEffectiveEventTime; if (fRelativeTiming) time += fTime2[fCurrentIndex]; else time = fTime2[fCurrentIndex]; // consider the earliest one - ll.fEffectiveEventTime = std::min(time, ll.fEffectiveEventTime); + fEffectiveEventTime = std::min(time, fEffectiveEventTime); } else { - ll.fEffectiveEventTime = current_simulation_time; + fEffectiveEventTime = current_simulation_time; } // weights @@ -130,6 +125,6 @@ void GateGANPairSource::GeneratePrimariesPair(G4Event *event, } // Vertex - AddOnePrimaryVertex(event, position, direction, energy, - ll.fEffectiveEventTime, w); + AddOnePrimaryVertex(event, position, direction, energy, fEffectiveEventTime, + w); } diff --git a/core/opengate_core/opengate_lib/GateGANSource.cpp b/core/opengate_core/opengate_lib/GateGANSource.cpp index 448afa6db..8996edac2 100644 --- a/core/opengate_core/opengate_lib/GateGANSource.cpp +++ b/core/opengate_core/opengate_lib/GateGANSource.cpp @@ -48,9 +48,8 @@ void GateGANSource::InitializeUserInfo(py::dict &user_info) { // AAManager is already set in GenericSource BUT MUST be iso direction here? auto d = py::dict(user_info["direction"]); auto dd = DictToMap(d["angular_acceptance"]); - auto &ll = GetThreadLocalDataGenericSource(); - ll.fAAManager->Initialize(dd, true); - ll.fSPS->SetAAManager(ll.fAAManager); + fAAManager->Initialize(dd, true); + fSPS->SetAAManager(fAAManager); // energy threshold mode auto s = DictGetStr(user_info, "skip_policy"); @@ -68,12 +67,12 @@ void GateGANSource::InitializeUserInfo(py::dict &user_info) { } // FIXME generic ion ? Should be possible, not tested - if (ll.fInitGenericIon) { + if (fInitGenericIon) { Fatal("Sorry, generic ion is not implemented with GAN source"); } // FIXME confine ? - if (ll.fInitConfine) { + if (fInitConfine) { Fatal("Sorry, confine is not implemented with GAN source"); } } @@ -87,13 +86,11 @@ void GateGANSource::PrepareNextRun() { GateVSource::PrepareNextRun(); // This global transformation is given to the SPS that will // generate particles in the correct coordinate system - auto &l = fThreadLocalData.Get(); - auto &lll = GetThreadLocalDataGenericSource(); - auto *pos = lll.fSPS->GetPosDist(); - pos->SetCentreCoords(l.fGlobalTranslation); + auto *pos = fSPS->GetPosDist(); + pos->SetCentreCoords(fGlobalTranslation); // orientation according to mother volume - auto rotation = l.fGlobalRotation; + auto rotation = fGlobalRotation; G4ThreeVector r1(rotation(0, 0), rotation(1, 0), rotation(2, 0)); G4ThreeVector r2(rotation(0, 1), rotation(1, 1), rotation(2, 1)); pos->SetPosRot1(r1); @@ -161,15 +158,13 @@ void GateGANSource::GeneratePrimaries(G4Event *event, fCurrentIndex++; // update the number of generated event - auto &l = fThreadLocalData.Get(); - l.fNumberOfGeneratedEvents++; + fNumberOfGeneratedEvents++; } void GateGANSource::GenerateOnePrimary(G4Event *event, double current_simulation_time) { // If AA (Angular Acceptance) is enabled, we perform rejection - auto &l = GetThreadLocalDataGenericSource(); - if (l.fAAManager->IsEnabled()) + if (fAAManager->IsEnabled()) return GenerateOnePrimaryWithAA(event, current_simulation_time); // If no AA, we loop until energy is acceptable, @@ -177,19 +172,19 @@ void GateGANSource::GenerateOnePrimary(G4Event *event, G4ThreeVector position; G4ThreeVector direction; double energy = 0; - l.fCurrentZeroEvents = 0; - l.fCurrentSkippedEvents = 0; + fCurrentZeroEvents = 0; + fCurrentSkippedEvents = 0; while (energy == 0 && fCurrentIndex < fCurrentBatchSize) { // position // (if it is not set by the GAN, we may avoid to sample at each iteration) - if (fPosition_is_set_by_GAN || l.fCurrentSkippedEvents == 0) + if (fPosition_is_set_by_GAN || fCurrentSkippedEvents == 0) position = GeneratePrimariesPosition(); // FIXME change position is not set by GAN // direction // (if it is not set by the GAN, we may avoid to sample at each iteration) - if (fDirection_is_set_by_GAN || l.fCurrentSkippedEvents == 0) + if (fDirection_is_set_by_GAN || fCurrentSkippedEvents == 0) direction = GeneratePrimariesDirection(); // energy @@ -202,11 +197,11 @@ void GateGANSource::GenerateOnePrimary(G4Event *event, if (fSkipEnergyPolicy == SEPolicyType::AAZeroEnergy || fCurrentIndex >= fCurrentBatchSize - 1) { energy = -1; - l.fCurrentZeroEvents = 1; + fCurrentZeroEvents = 1; } else { // energy is not ok, we skip the event and try the next one energy = 0; - l.fCurrentSkippedEvents++; + fCurrentSkippedEvents++; fCurrentIndex++; } } @@ -216,7 +211,7 @@ void GateGANSource::GenerateOnePrimary(G4Event *event, // we continue with a zeroE event if (energy == -1 || fCurrentIndex >= fCurrentBatchSize - 1) { energy = 0; - l.fCurrentZeroEvents = 1; + fCurrentZeroEvents = 1; } // timing @@ -230,8 +225,6 @@ void GateGANSource::GenerateOnePrimary(G4Event *event, } G4ThreeVector GateGANSource::GeneratePrimariesPosition() { - auto &l = GetThreadLocalData(); - auto &ll = GetThreadLocalDataGenericSource(); G4ThreeVector position; if (fPosition_is_set_by_GAN) { position = @@ -239,15 +232,14 @@ G4ThreeVector GateGANSource::GeneratePrimariesPosition() { fPositionZ[fCurrentIndex]); position = fLocalRotation * position + fLocalTranslation; // FIXME // move position according to mother volume - position = l.fGlobalRotation * position + l.fGlobalTranslation; + position = fGlobalRotation * position + fGlobalTranslation; } else - position = ll.fSPS->GetPosDist()->VGenerateOne(); + position = fSPS->GetPosDist()->VGenerateOne(); return position; } G4ThreeVector GateGANSource::GeneratePrimariesDirection() { G4ThreeVector direction; - auto &ll = GetThreadLocalDataGenericSource(); if (fDirection_is_set_by_GAN) { direction = G4ParticleMomentum(fDirectionX[fCurrentIndex], fDirectionY[fCurrentIndex], @@ -255,54 +247,50 @@ G4ThreeVector GateGANSource::GeneratePrimariesDirection() { // normalize (needed) direction = direction / direction.mag(); // move according to mother volume - auto &l = fThreadLocalData.Get(); - direction = l.fGlobalRotation * direction; + direction = fGlobalRotation * direction; } else - direction = ll.fSPS->GetAngDist()->GenerateOne(); + direction = fSPS->GetAngDist()->GenerateOne(); return direction; } double GateGANSource::GeneratePrimariesEnergy() { double energy; - auto &ll = GetThreadLocalDataGenericSource(); if (fEnergy_is_set_by_GAN) energy = fEnergy[fCurrentIndex]; else - energy = ll.fSPS->GetEneDist()->VGenerateOne(fParticleDefinition); + energy = fSPS->GetEneDist()->VGenerateOne(fParticleDefinition); return energy; } double GateGANSource::GeneratePrimariesTime(double current_simulation_time) { - auto &ll = fThreadLocalDataGenericSource.Get(); - if (!fTime_is_set_by_GAN) { - ll.fEffectiveEventTime = current_simulation_time; - return ll.fEffectiveEventTime; + fEffectiveEventTime = current_simulation_time; + return fEffectiveEventTime; } - if (ll.fCurrentZeroEvents > 0) { - ll.fEffectiveEventTime = current_simulation_time; - return ll.fEffectiveEventTime; + if (fCurrentZeroEvents > 0) { + fEffectiveEventTime = current_simulation_time; + return fEffectiveEventTime; } // if the time is managed by the GAN, it can be relative or absolute. if (fRelativeTiming) { // update the real time (important as the event is in the // future according to the current_simulation_time) - UpdateEffectiveEventTime(current_simulation_time, ll.fCurrentSkippedEvents); + UpdateEffectiveEventTime(current_simulation_time, fCurrentSkippedEvents); } // Get the time from the GAN except if it is a zeroE - if (ll.fCurrentZeroEvents > 0) - ll.fEffectiveEventTime = current_simulation_time; + if (fCurrentZeroEvents > 0) + fEffectiveEventTime = current_simulation_time; else { if (fRelativeTiming) - ll.fEffectiveEventTime += fTime[fCurrentIndex]; + fEffectiveEventTime += fTime[fCurrentIndex]; else - ll.fEffectiveEventTime = fTime[fCurrentIndex]; + fEffectiveEventTime = fTime[fCurrentIndex]; } - return ll.fEffectiveEventTime; + return fEffectiveEventTime; } double GateGANSource::GeneratePrimariesWeight() { @@ -317,10 +305,9 @@ void GateGANSource::GenerateOnePrimaryWithAA(G4Event *event, G4ThreeVector direction; double energy = 0; bool cont = true; - auto &l = GetThreadLocalDataGenericSource(); - l.fCurrentZeroEvents = 0; - l.fCurrentSkippedEvents = 0; - l.fAAManager->StartAcceptLoop(); + fCurrentZeroEvents = 0; + fCurrentSkippedEvents = 0; + fAAManager->StartAcceptLoop(); while (cont) { // position @@ -330,17 +317,17 @@ void GateGANSource::GenerateOnePrimaryWithAA(G4Event *event, direction = GeneratePrimariesDirection(); // check AA - bool accept_angle = l.fAAManager->TestIfAccept(position, direction); + bool accept_angle = fAAManager->TestIfAccept(position, direction); if (!accept_angle && - l.fAAManager->GetPolicy() == GateAcceptanceAngleManager::AAZeroEnergy) { + fAAManager->GetPolicy() == GateAcceptanceAngleManager::AAZeroEnergy) { energy = 0; cont = false; continue; // stop here } if (!accept_angle && - l.fAAManager->GetPolicy() == GateAcceptanceAngleManager::AASkipEvent) { - l.fCurrentSkippedEvents++; + fAAManager->GetPolicy() == GateAcceptanceAngleManager::AASkipEvent) { + fCurrentSkippedEvents++; fCurrentIndex++; continue; // no need to check energy now } @@ -355,11 +342,11 @@ void GateGANSource::GenerateOnePrimaryWithAA(G4Event *event, if (fSkipEnergyPolicy == SEPolicyType::AAZeroEnergy) { cont = false; energy = 0; - l.fCurrentZeroEvents = 1; + fCurrentZeroEvents = 1; continue; // stop here } else { // energy is not ok, we skip the event and try the next one - l.fCurrentSkippedEvents++; + fCurrentSkippedEvents++; fCurrentIndex++; } } @@ -368,7 +355,7 @@ void GateGANSource::GenerateOnePrimaryWithAA(G4Event *event, if (fCurrentIndex >= fCurrentBatchSize) { cont = false; energy = 0; - l.fCurrentZeroEvents = 1; + fCurrentZeroEvents = 1; continue; // stop here } else { cont = false; diff --git a/core/opengate_core/opengate_lib/GateGenericSource.cpp b/core/opengate_core/opengate_lib/GateGenericSource.cpp index 194f006d0..b2785b682 100644 --- a/core/opengate_core/opengate_lib/GateGenericSource.cpp +++ b/core/opengate_core/opengate_lib/GateGenericSource.cpp @@ -32,6 +32,14 @@ GateGenericSource::GateGenericSource() : GateVSource() { fDirectionRelativeToAttachedVolume = false; fUserParticleLifeTime = -1; fBackToBackMode = false; + fSPS = nullptr; + fAAManager = nullptr; + fFDManager = nullptr; + fInitConfine = false; + fInitGenericIon = false; + fEffectiveEventTime = -1; + fCurrentSkippedEvents = 0; + fCurrentZeroEvents = 0; } GateGenericSource::~GateGenericSource() { @@ -39,25 +47,18 @@ GateGenericSource::~GateGenericSource() { // I dont know exactly why. // Maybe because it has been created in a thread which // can be different from the thread that delete. - auto &l = fThreadLocalDataGenericSource.Get(); - if (l.fAAManager != nullptr) { - // delete l.fAAManager; + if (fAAManager != nullptr) { + // delete fAAManager; } // delete fSPS; } -GateGenericSource::threadLocalGenericSource & -GateGenericSource::GetThreadLocalDataGenericSource() const { - return fThreadLocalDataGenericSource.Get(); -} - void GateGenericSource::CleanWorkerThread() { // Not used yet. Maybe later to clean local data in a thread. } void GateGenericSource::CreateSPS() { - auto &l = fThreadLocalDataGenericSource.Get(); - l.fSPS = new GateSingleParticleSource(fAttachedToVolumeName); + fSPS = new GateSingleParticleSource(fAttachedToVolumeName); } void GateGenericSource::SetEnergyCDF(const std::vector &cdf) { @@ -140,19 +141,18 @@ void GateGenericSource::UpdateActivityWithTAC(const double time) { double GateGenericSource::PrepareNextTime(const double current_simulation_time, double NumberOfGeneratedEvents) { - auto &ll = GetThreadLocalDataGenericSource(); // initialization of the effective event time (it can be in the // future according to the current_simulation_time) - if (ll.fEffectiveEventTime < current_simulation_time) { - ll.fEffectiveEventTime = current_simulation_time; + if (fEffectiveEventTime < current_simulation_time) { + fEffectiveEventTime = current_simulation_time; } - fTotalSkippedEvents += ll.fCurrentSkippedEvents; // FIXME lock ? - fTotalZeroEvents += ll.fCurrentZeroEvents; - ll.fCurrentZeroEvents = 0; - const auto cse = ll.fCurrentSkippedEvents; - ll.fCurrentSkippedEvents = 0; + fTotalSkippedEvents += fCurrentSkippedEvents; + fTotalZeroEvents += fCurrentZeroEvents; + fCurrentZeroEvents = 0; + const auto cse = fCurrentSkippedEvents; + fCurrentSkippedEvents = 0; - return GateVSource::PrepareNextTime(ll.fEffectiveEventTime, + return GateVSource::PrepareNextTime(fEffectiveEventTime, NumberOfGeneratedEvents + cse); } @@ -163,13 +163,11 @@ void GateGenericSource::PrepareNextRun() { // This global transformation is given to the SPS that will // generate particles in the correct coordinate system - auto &l = GetThreadLocalData(); - auto &ll = GetThreadLocalDataGenericSource(); - auto *pos = ll.fSPS->GetPosDist(); - pos->SetCentreCoords(l.fGlobalTranslation); + auto *pos = fSPS->GetPosDist(); + pos->SetCentreCoords(fGlobalTranslation); // orientation according to mother volume - const auto rotation = l.fGlobalRotation; + const auto rotation = fGlobalRotation; const G4ThreeVector r1(rotation(0, 0), rotation(1, 0), rotation(2, 0)); const G4ThreeVector r2(rotation(0, 1), rotation(1, 1), rotation(2, 1)); pos->SetPosRot1(r1); @@ -177,10 +175,10 @@ void GateGenericSource::PrepareNextRun() { // For the direction, the orientation may or may not be // relative to the volume according to the user option - auto *ang = ll.fSPS->GetAngDist(); + auto *ang = fSPS->GetAngDist(); ang->fDirectionRelativeToAttachedVolume = fDirectionRelativeToAttachedVolume; - ang->fGlobalRotation = l.fGlobalRotation; - ang->fGlobalTranslation = l.fGlobalTranslation; + ang->fGlobalRotation = fGlobalRotation; + ang->fGlobalTranslation = fGlobalTranslation; if (fangType == "momentum" && fDirectionRelativeToAttachedVolume) { const auto new_d = rotation * fInitializeMomentum; ang->SetParticleMomentumDirection(new_d); @@ -189,7 +187,7 @@ void GateGenericSource::PrepareNextRun() { if (fangType == "focused" && fDirectionRelativeToAttachedVolume) { const auto vec_f = fInitializeFocusPoint - fInitTranslation; const auto rot_f = rotation * vec_f; - const auto new_f = rot_f + l.fGlobalTranslation; + const auto new_f = rot_f + fGlobalTranslation; ang->SetFocusPoint(new_f); ang->fDirectionRelativeToAttachedVolume = false; } @@ -197,54 +195,51 @@ void GateGenericSource::PrepareNextRun() { void GateGenericSource::UpdateEffectiveEventTime( const double current_simulation_time, - const unsigned long skipped_particle) const { - auto &ll = GetThreadLocalDataGenericSource(); + const unsigned long skipped_particle) { unsigned long n = 0; - ll.fEffectiveEventTime = current_simulation_time; - while (n < skipped_particle && ll.fEffectiveEventTime < fEndTime) { - ll.fEffectiveEventTime = - ll.fEffectiveEventTime - log(G4UniformRand()) * (1.0 / fActivity); + fEffectiveEventTime = current_simulation_time; + while (n < skipped_particle && fEffectiveEventTime < fEndTime) { + fEffectiveEventTime = + fEffectiveEventTime - log(G4UniformRand()) * (1.0 / fActivity); n++; } } void GateGenericSource::GeneratePrimaries( G4Event *event, const double current_simulation_time) { - auto &ll = GetThreadLocalDataGenericSource(); // Generic ion cannot be created at initialization. // It must be created the first time we get there - if (ll.fInitGenericIon) { + if (fInitGenericIon) { auto *ion_table = G4IonTable::GetIonTable(); auto *ion = ion_table->GetIon(fZ, fA, fE); - ll.fSPS->SetParticleDefinition(ion); + fSPS->SetParticleDefinition(ion); SetLifeTime(ion); - ll.fInitGenericIon = false; // only the first time + fInitGenericIon = false; // only the first time } // Confine cannot be initialized at initialization (because need all volumes // to be created) It must be set here, the first time we get there - if (ll.fInitConfine) { - auto *pos = ll.fSPS->GetPosDist(); + if (fInitConfine) { + auto *pos = fSPS->GetPosDist(); pos->ConfineSourceToVolume(fConfineVolume); - ll.fInitConfine = false; + fInitConfine = false; } // sample the particle properties with SingleParticleSource // (the acceptance angle or forced direction is included) - ll.fSPS->SetParticleTime(current_simulation_time); - ll.fSPS->GeneratePrimaryVertex(event); + fSPS->SetParticleTime(current_simulation_time); + fSPS->GeneratePrimaryVertex(event); // update the time according to skipped events - ll.fEffectiveEventTime = current_simulation_time; - if (ll.fAAManager->IsEnabled()) { - if (ll.fAAManager->GetPolicy() == GateAcceptanceAngleManager::AASkipEvent) { + fEffectiveEventTime = current_simulation_time; + if (fAAManager->IsEnabled()) { + if (fAAManager->GetPolicy() == GateAcceptanceAngleManager::AASkipEvent) { UpdateEffectiveEventTime(current_simulation_time, - ll.fAAManager->GetNumberOfNotAcceptedEvents()); - ll.fCurrentSkippedEvents = ll.fAAManager->GetNumberOfNotAcceptedEvents(); - event->GetPrimaryVertex(0)->SetT0(ll.fEffectiveEventTime); + fAAManager->GetNumberOfNotAcceptedEvents()); + fCurrentSkippedEvents = fAAManager->GetNumberOfNotAcceptedEvents(); + event->GetPrimaryVertex(0)->SetT0(fEffectiveEventTime); } else { - ll.fCurrentZeroEvents = - ll.fAAManager->GetNumberOfNotAcceptedEvents(); // 1 or 0 + fCurrentZeroEvents = fAAManager->GetNumberOfNotAcceptedEvents(); // 1 or 0 } } @@ -262,19 +257,17 @@ void GateGenericSource::GeneratePrimaries( } } - auto &l = GetThreadLocalData(); - l.fNumberOfGeneratedEvents++; + fNumberOfGeneratedEvents++; } void GateGenericSource::InitializeParticle(py::dict &user_info) { - auto &ll = fThreadLocalDataGenericSource.Get(); std::string pname = DictGetStr(user_info, "particle"); // Is the particle an ion (name start with ion) ? if (pname.rfind("ion", 0) == 0) { InitializeIon(user_info); return; } - ll.fInitGenericIon = false; + fInitGenericIon = false; // Is the particle a back to back ? if (pname.rfind("back_to_back") == 0) { InitializeBackToBackMode(user_info); @@ -287,7 +280,7 @@ void GateGenericSource::InitializeParticle(py::dict &user_info) { if (fParticleDefinition == nullptr) { Fatal("Cannot find the particle '" + pname + "'."); } - ll.fSPS->SetParticleDefinition(fParticleDefinition); + fSPS->SetParticleDefinition(fParticleDefinition); SetLifeTime(fParticleDefinition); } @@ -296,24 +289,22 @@ void GateGenericSource::InitializeIon(py::dict &user_info) { fA = DictGetInt(u, "A"); fZ = DictGetInt(u, "Z"); fE = DictGetDouble(u, "E"); - auto &ll = fThreadLocalDataGenericSource.Get(); - ll.fInitGenericIon = true; + fInitGenericIon = true; } void GateGenericSource::InitializeBackToBackMode(py::dict &user_info) { - auto &ll = fThreadLocalDataGenericSource.Get(); auto u = py::dict(user_info["direction"]); bool accolinearityFlag = DictGetBool(u, "accolinearity_flag"); - ll.fSPS->SetBackToBackMode(true, accolinearityFlag); + fSPS->SetBackToBackMode(true, accolinearityFlag); if (accolinearityFlag == true) { // Change the value if user provided one. double accolinearityFWHM = DictGetDouble(u, "accolinearity_fwhm"); - ll.fSPS->SetAccolinearityFWHM(accolinearityFWHM); + fSPS->SetAccolinearityFWHM(accolinearityFWHM); } // this is photon auto *particle_table = G4ParticleTable::GetParticleTable(); fParticleDefinition = particle_table->FindParticle("gamma"); - ll.fSPS->SetParticleDefinition(fParticleDefinition); + fSPS->SetParticleDefinition(fParticleDefinition); // The energy is fixed to 511 keV in the python side } @@ -325,9 +316,8 @@ void GateGenericSource::InitializePosition(py::dict puser_info) { * New interface -> point box sphere disc (later: ellipse) * translation rotation size radius */ - auto &ll = fThreadLocalDataGenericSource.Get(); auto user_info = py::dict(puser_info["position"]); - auto *pos = ll.fSPS->GetPosDist(); + auto *pos = fSPS->GetPosDist(); auto pos_type = DictGetStr(user_info, "type"); std::vector l = {"sphere", "point", "box", "disc", "cylinder"}; CheckIsIn(pos_type, l); @@ -382,7 +372,7 @@ void GateGenericSource::InitializePosition(py::dict puser_info) { auto v = DictGetStr(user_info, "confine"); if (v != "None") { fConfineVolume = v; - ll.fInitConfine = true; + fInitConfine = true; } } } @@ -395,9 +385,8 @@ void GateGenericSource::InitializeDirection(py::dict puser_info) { * New ones: iso, focus, direction * (Later: beam, user defined) */ - auto &ll = fThreadLocalDataGenericSource.Get(); auto user_info = py::dict(puser_info["direction"]); - auto *ang = ll.fSPS->GetAngDist(); + auto *ang = fSPS->GetAngDist(); auto ang_type = DictGetStr(user_info, "type"); fangType = ang_type; std::vector llt = {"iso", "histogram", "momentum", "focused", @@ -474,24 +463,23 @@ void GateGenericSource::InitializeDirection(py::dict puser_info) { auto dd = DictToMap(d["angular_acceptance"]); const auto is_valid_type = ang->GetDistType() == "iso" || ang->GetDistType() == "user"; - ll.fAAManager = new GateAcceptanceAngleManager; - ll.fAAManager->Initialize(dd, is_valid_type); - ll.fSPS->SetAAManager(ll.fAAManager); + fAAManager = new GateAcceptanceAngleManager; + fAAManager->Initialize(dd, is_valid_type); + fSPS->SetAAManager(fAAManager); // set Forced Direction - ll.fFDManager = new GateForcedDirectionManager; - ll.fFDManager->Initialize(dd, ang->GetDistType() == "iso"); - ll.fSPS->SetFDManager(ll.fFDManager); + fFDManager = new GateForcedDirectionManager; + fFDManager->Initialize(dd, ang->GetDistType() == "iso"); + fSPS->SetFDManager(fFDManager); } void GateGenericSource::InitializePolarization(py::dict puser_info) { // Set the polarization - auto &ll = fThreadLocalDataGenericSource.Get(); auto polarization = DictGetVecDouble(puser_info, "polarization"); if (polarization.size() == 3) { auto polarisation_tree_vector = G4ThreeVector(polarization[0], polarization[1], polarization[2]); - ll.fSPS->SetPolarization(polarisation_tree_vector); + fSPS->SetPolarization(polarisation_tree_vector); } } @@ -505,9 +493,8 @@ void GateGenericSource::InitializeEnergy(py::dict puser_info) { * New interface: mono gauss // later 'user' * */ - auto &ll = fThreadLocalDataGenericSource.Get(); auto user_info = py::dict(puser_info["energy"]); - auto *ene = ll.fSPS->GetEneDist(); + auto *ene = fSPS->GetEneDist(); auto ene_type = DictGetStr(user_info, "type"); auto is_cdf = DictGetBool(user_info, "is_cdf"); @@ -725,8 +712,7 @@ void GateGenericSource::Visualize() const { G4VisManager *vis_manager = G4VisManager::GetInstance(); auto *scene = vis_manager->GetCurrentScene(); - auto &ll = GetThreadLocalDataGenericSource(); - auto *pos = ll.fSPS->GetPosDist(); + auto *pos = fSPS->GetPosDist(); auto *point_cloud = new PosPointCloud(fVisColour, fVisSize); for (int k = 0; k < fVisCount; ++k) { diff --git a/core/opengate_core/opengate_lib/GateGenericSource.h b/core/opengate_core/opengate_lib/GateGenericSource.h index cfd3a3b90..2984d1d49 100644 --- a/core/opengate_core/opengate_lib/GateGenericSource.h +++ b/core/opengate_core/opengate_lib/GateGenericSource.h @@ -94,25 +94,17 @@ class GateGenericSource : public GateVSource { // source, eg: needed for motion actor bool fDirectionRelativeToAttachedVolume; - // thread local structure - struct threadLocalGenericSource { - GateSingleParticleSource *fSPS = nullptr; - GateAcceptanceAngleManager *fAAManager = nullptr; - GateForcedDirectionManager *fFDManager = nullptr; - bool fInitConfine = false; - bool fInitGenericIon = false; - double fEffectiveEventTime = -1; - unsigned long fCurrentSkippedEvents = 0; - unsigned long fCurrentZeroEvents = 0; - }; - G4Cache fThreadLocalDataGenericSource; - - // sum of all threads + GateSingleParticleSource *fSPS = nullptr; + GateAcceptanceAngleManager *fAAManager = nullptr; + GateForcedDirectionManager *fFDManager = nullptr; + bool fInitConfine = false; + bool fInitGenericIon = false; + double fEffectiveEventTime = -1; + unsigned long fCurrentSkippedEvents = 0; + unsigned long fCurrentZeroEvents = 0; unsigned long fTotalSkippedEvents = 0; unsigned long fTotalZeroEvents = 0; - threadLocalGenericSource &GetThreadLocalDataGenericSource() const; - // if confine is used, must be defined after the initialization // bool fInitConfine; std::string fConfineVolume; @@ -142,7 +134,7 @@ class GateGenericSource : public GateVSource { void UpdateActivity(double time) override; void UpdateEffectiveEventTime(double current_simulation_time, - unsigned long skipped_particle) const; + unsigned long skipped_particle); }; #endif // GateGenericSource_h diff --git a/core/opengate_core/opengate_lib/GateLastVertexSource.cpp b/core/opengate_core/opengate_lib/GateLastVertexSource.cpp index bec27eb1a..ab61b05c6 100644 --- a/core/opengate_core/opengate_lib/GateLastVertexSource.cpp +++ b/core/opengate_core/opengate_lib/GateLastVertexSource.cpp @@ -85,7 +85,6 @@ void GateLastVertexSource::GenerateOnePrimary(G4Event *event, G4double weight = containerToSplit.GetWeight(); fProcessToSplit = containerToSplit.GetProcessNameToSplit(); - auto &l = fThreadLocalData.Get(); auto *particle_table = G4ParticleTable::GetParticleTable(); auto *fParticleDefinition = particle_table->FindParticle(particleName); auto *particle = new G4PrimaryParticle(fParticleDefinition); diff --git a/core/opengate_core/opengate_lib/GatePencilBeamSource.cpp b/core/opengate_core/opengate_lib/GatePencilBeamSource.cpp index ee6a094f7..3ea940e8d 100644 --- a/core/opengate_core/opengate_lib/GatePencilBeamSource.cpp +++ b/core/opengate_core/opengate_lib/GatePencilBeamSource.cpp @@ -10,21 +10,16 @@ #include "GateHelpersDict.h" #include -GatePencilBeamSource::GatePencilBeamSource() : GateGenericSource() {} +GatePencilBeamSource::GatePencilBeamSource() : GateGenericSource() { + fSPS_PB = nullptr; +} GatePencilBeamSource::~GatePencilBeamSource() = default; -GatePencilBeamSource::threadLocalPencilBeamSource & -GatePencilBeamSource::GetThreadLocalDataPencilBeamSource() { - return fThreadLocalDataPencilBeamSource.Get(); -} - void GatePencilBeamSource::CreateSPS() { - auto &lll = GetThreadLocalDataPencilBeamSource(); - lll.fSPS_PB = new GateSingleParticleSourcePencilBeam(std::string(), - fAttachedToVolumeName); - auto &ll = GetThreadLocalDataGenericSource(); - ll.fSPS = lll.fSPS_PB; + fSPS_PB = new GateSingleParticleSourcePencilBeam(std::string(), + fAttachedToVolumeName); + fSPS = fSPS_PB; } void GatePencilBeamSource::PrepareNextRun() { @@ -33,9 +28,7 @@ void GatePencilBeamSource::PrepareNextRun() { GateVSource::PrepareNextRun(); // This global transformation is given to the SPS that will // generate particles in the correct coordinate system - auto &l = fThreadLocalData.Get(); - auto &lll = GetThreadLocalDataPencilBeamSource(); - lll.fSPS_PB->SetSourceRotTransl(l.fGlobalTranslation, l.fGlobalRotation); + fSPS_PB->SetSourceRotTransl(fGlobalTranslation, fGlobalRotation); } void GatePencilBeamSource::InitializeDirection(py::dict puser_info) { @@ -44,16 +37,14 @@ void GatePencilBeamSource::InitializeDirection(py::dict puser_info) { auto dir_info = py::dict(puser_info["direction"]); auto x_param = DictGetVecDouble(dir_info, "partPhSp_x"); auto y_param = DictGetVecDouble(dir_info, "partPhSp_y"); - auto &lll = GetThreadLocalDataPencilBeamSource(); - lll.fSPS_PB->SetPBSourceParam(x_param, y_param); + fSPS_PB->SetPBSourceParam(x_param, y_param); // angle acceptance ? auto d = py::dict(puser_info["direction"]); auto dd = DictToMap(d["angular_acceptance"]); - auto &l = fThreadLocalDataGenericSource.Get(); - l.fAAManager = new GateAcceptanceAngleManager; - l.fAAManager->Initialize(dd, false); - if (l.fAAManager->IsEnabled()) { + fAAManager = new GateAcceptanceAngleManager; + fAAManager->Initialize(dd, false); + if (fAAManager->IsEnabled()) { Fatal("Sorry, cannot use Acceptance Angle with Pencil Beam source (yet)."); } } diff --git a/core/opengate_core/opengate_lib/GatePencilBeamSource.h b/core/opengate_core/opengate_lib/GatePencilBeamSource.h index 99523e5b3..53d27d1ed 100644 --- a/core/opengate_core/opengate_lib/GatePencilBeamSource.h +++ b/core/opengate_core/opengate_lib/GatePencilBeamSource.h @@ -26,15 +26,7 @@ class GatePencilBeamSource : public GateGenericSource { void PrepareNextRun() override; protected: - // thread local structure - struct threadLocalPencilBeamSource { - // the fSPS will be a GateSingleParticleSourcePencilBeam - // we store the two pointers fSPS and fSPS_PB to the same object - GateSingleParticleSourcePencilBeam *fSPS_PB = nullptr; - }; - G4Cache fThreadLocalDataPencilBeamSource; - - threadLocalPencilBeamSource &GetThreadLocalDataPencilBeamSource(); + GateSingleParticleSourcePencilBeam *fSPS_PB = nullptr; void CreateSPS() override; }; diff --git a/core/opengate_core/opengate_lib/GatePhaseSpaceSource.cpp b/core/opengate_core/opengate_lib/GatePhaseSpaceSource.cpp index 63fa819b7..f29d9b795 100644 --- a/core/opengate_core/opengate_lib/GatePhaseSpaceSource.cpp +++ b/core/opengate_core/opengate_lib/GatePhaseSpaceSource.cpp @@ -20,23 +20,12 @@ GatePhaseSpaceSource::GatePhaseSpaceSource() : GateVSource() { fVerbose = false; fParticleTable = nullptr; fUseParticleTypeFromFile = false; - auto &l = fThreadLocalData.Get(); - auto &ll = fThreadLocalDataPhsp.Get(); - ll.fParticleDefinition = nullptr; + fParticleDefinition = nullptr; } -GatePhaseSpaceSource::~GatePhaseSpaceSource() { - // It seems that this is required to prevent seg fault at the end - // I don't understand why - auto &l = fThreadLocalData.Get(); - auto &ll = fThreadLocalDataPhsp.Get(); -} +GatePhaseSpaceSource::~GatePhaseSpaceSource() = default; void GatePhaseSpaceSource::InitializeUserInfo(py::dict &user_info) { - auto &l = fThreadLocalData.Get(); - auto &ll = fThreadLocalDataPhsp.Get(); - // the following initialize all GenericSource options - // and the SPS GateSingleParticleSource GateVSource::InitializeUserInfo(user_info); // global (world) or local (mother volume) coordinate system @@ -54,24 +43,24 @@ void GatePhaseSpaceSource::InitializeUserInfo(py::dict &user_info) { fUseParticleTypeFromFile = true; else { fUseParticleTypeFromFile = false; - ll.fParticleDefinition = fParticleTable->FindParticle(pname); - if (ll.fParticleDefinition == nullptr) { + fParticleDefinition = fParticleTable->FindParticle(pname); + if (fParticleDefinition == nullptr) { Fatal("GatePhaseSpaceSource: PDGCode not found. Aborting."); } - fCharge = static_cast(ll.fParticleDefinition->GetPDGCharge()); - fMass = static_cast(ll.fParticleDefinition->GetPDGMass()); + fCharge = static_cast(fParticleDefinition->GetPDGCharge()); + fMass = static_cast(fParticleDefinition->GetPDGMass()); } // Init - l.fNumberOfGeneratedEvents = 0; - ll.fCurrentIndex = 0; - ll.fCurrentBatchSize = 0; + fNumberOfGeneratedEvents = 0; + fCurrentIndex = 0; + fCurrentBatchSize = 0; - ll.fGenerateUntilNextPrimary = + fGenerateUntilNextPrimary = DictGetBool(user_info, "generate_until_next_primary"); - ll.fPrimaryLowerEnergyThreshold = static_cast( - DictGetDouble(user_info, "primary_lower_energy_threshold")); - ll.fPrimaryPDGCode = DictGetInt(user_info, "primary_PDGCode"); + fPrimaryLowerEnergyThreshold = + DictGetDouble(user_info, "primary_lower_energy_threshold"); + fPrimaryPDGCode = DictGetInt(user_info, "primary_PDGCode"); } void GatePhaseSpaceSource::PrepareNextRun() { @@ -79,10 +68,8 @@ void GatePhaseSpaceSource::PrepareNextRun() { GateVSource::PrepareNextRun(); } -void GatePhaseSpaceSource::SetGeneratorFunction( - ParticleGeneratorType &f) const { - auto &ll = fThreadLocalDataPhsp.Get(); - ll.fGenerator = f; +void GatePhaseSpaceSource::SetGeneratorFunction(ParticleGeneratorType &f) { + fGenerator = f; } void GatePhaseSpaceSource::GenerateBatchOfParticles() { @@ -90,34 +77,31 @@ void GatePhaseSpaceSource::GenerateBatchOfParticles() { // (does not seem needed) // py::gil_scoped_acquire acquire; - // This function (ll.fGenerator) is defined on Python side + // This function (fGenerator) is defined on Python side // It fills all values needed for the particles (position, dir, energy, etc) // Alternative: build vector of G4ThreeVector in GenerateBatchOfParticles ? // (unsure if it is faster) - auto &ll = fThreadLocalDataPhsp.Get(); - ll.fCurrentBatchSize = ll.fGenerator(this, G4Threading::G4GetThreadId()); - ll.fCurrentIndex = 0; + fCurrentBatchSize = fGenerator(this, G4Threading::G4GetThreadId()); + fCurrentIndex = 0; } void GatePhaseSpaceSource::GeneratePrimaries(G4Event *event, double current_simulation_time) { - auto &l = fThreadLocalData.Get(); - auto &ll = fThreadLocalDataPhsp.Get(); // check if we should simulate until next primary // in this case, generate until a second primary is in the list, excluding the // second primary // If batch is empty, we generate some particles, could happen at the first // execution - if (ll.fCurrentBatchSize == 0) { + if (fCurrentBatchSize == 0) { GenerateBatchOfParticles(); } - if (ll.fGenerateUntilNextPrimary) { + if (fGenerateUntilNextPrimary) { int num_primaries = 0; while (num_primaries <= 2) { // If batch is empty, we generate some particles - if (ll.fCurrentIndex >= ll.fCurrentBatchSize) + if (fCurrentIndex >= fCurrentBatchSize) GenerateBatchOfParticles(); // check if next particle is primary @@ -132,26 +116,26 @@ void GatePhaseSpaceSource::GeneratePrimaries(G4Event *event, GenerateOnePrimary(event, current_simulation_time); // update the root file index; - ll.fCurrentIndex++; + fCurrentIndex++; } else break; } // update the number of generated event - l.fNumberOfGeneratedEvents++; + fNumberOfGeneratedEvents++; } else { // If batch is empty, we generate some particles - if (ll.fCurrentIndex >= ll.fCurrentBatchSize) + if (fCurrentIndex >= fCurrentBatchSize) GenerateBatchOfParticles(); // Go GenerateOnePrimary(event, current_simulation_time); // update the root file index; - ll.fCurrentIndex++; + fCurrentIndex++; // update the number of generated event - l.fNumberOfGeneratedEvents++; + fNumberOfGeneratedEvents++; } } @@ -168,33 +152,29 @@ G4ParticleMomentum GatePhaseSpaceSource::GenerateRandomDirection() { void GatePhaseSpaceSource::GenerateOnePrimary(G4Event *event, double current_simulation_time) { - auto &ll = fThreadLocalDataPhsp.Get(); - - auto position = G4ThreeVector(ll.fPositionX[ll.fCurrentIndex], - ll.fPositionY[ll.fCurrentIndex], - ll.fPositionZ[ll.fCurrentIndex]); + G4ThreeVector position(fPositionX[fCurrentIndex], fPositionY[fCurrentIndex], + fPositionZ[fCurrentIndex]); G4ParticleMomentum direction; if (fIsotropicMomentum == false) { - direction = G4ParticleMomentum(ll.fDirectionX[ll.fCurrentIndex], - ll.fDirectionY[ll.fCurrentIndex], - ll.fDirectionZ[ll.fCurrentIndex]); + direction = G4ParticleMomentum(fDirectionX[fCurrentIndex], + fDirectionY[fCurrentIndex], + fDirectionZ[fCurrentIndex]); } else { direction = GenerateRandomDirection(); } - auto energy = ll.fEnergy[ll.fCurrentIndex]; - auto weight = ll.fWeight[ll.fCurrentIndex]; + auto energy = fEnergy[fCurrentIndex]; + auto weight = fWeight[fCurrentIndex]; - // FIXME auto time = fTime[ll.fCurrentIndex]; + // FIXME auto time = fTime[fCurrentIndex]; // transform according to mother if (!fGlobalFag) { - auto &ls = fThreadLocalData.Get(); - position = ls.fGlobalRotation * position + ls.fGlobalTranslation; + position = fGlobalRotation * position + fGlobalTranslation; direction = direction / direction.mag(); - direction = ls.fGlobalRotation * direction; + direction = fGlobalRotation * direction; } // Create the final vertex AddOnePrimaryVertex(event, position, direction, energy, @@ -207,26 +187,25 @@ void GatePhaseSpaceSource::AddOnePrimaryVertex(G4Event *event, double energy, double time, double w) { auto *particle = new G4PrimaryParticle(); - auto &ll = fThreadLocalDataPhsp.Get(); if (fUseParticleTypeFromFile) { - auto pdg = ll.fPDGCode[ll.fCurrentIndex]; - if (ll.fPDGCode[ll.fCurrentIndex] != 0) { + auto pdg = fPDGCode[fCurrentIndex]; + if (fPDGCode[fCurrentIndex] != 0) { // find if particle exists - ll.fParticleDefinition = fParticleTable->FindParticle(pdg); + fParticleDefinition = fParticleTable->FindParticle(pdg); // if not, find if it is an ion - if (ll.fParticleDefinition == nullptr) { + if (fParticleDefinition == nullptr) { G4IonTable *ionTable = fParticleTable->GetIonTable(); - ll.fParticleDefinition = ionTable->GetIon(pdg); + fParticleDefinition = ionTable->GetIon(pdg); } - if (ll.fParticleDefinition == nullptr) { + if (fParticleDefinition == nullptr) { Fatal("GatePhaseSpaceSource: PDGCode not found. Aborting."); } - particle->SetParticleDefinition(ll.fParticleDefinition); + particle->SetParticleDefinition(fParticleDefinition); } else { Fatal("GatePhaseSpaceSource: PDGCode not available. Aborting."); } } else { - particle->SetParticleDefinition(ll.fParticleDefinition); + particle->SetParticleDefinition(fParticleDefinition); particle->SetMass(fMass); particle->SetCharge(fCharge); } @@ -240,8 +219,7 @@ void GatePhaseSpaceSource::AddOnePrimaryVertex(G4Event *event, // weights event->GetPrimaryVertex(0)->SetWeight(w); if (fVerbose) { - std::cout << "Particle PDGCode: " - << ll.fParticleDefinition->GetPDGEncoding() + std::cout << "Particle PDGCode: " << fParticleDefinition->GetPDGEncoding() << " Energy: " << energy << " Weight: " << w << " Position: " << position << " Direction: " << direction << " Time: " << time << " EventID: " << event->GetEventID() @@ -250,68 +228,58 @@ void GatePhaseSpaceSource::AddOnePrimaryVertex(G4Event *event, } void GatePhaseSpaceSource::SetPDGCodeBatch( - const py::array_t &fPDGCode) const { - auto &ll = fThreadLocalDataPhsp.Get(); - ll.fPDGCode = PyBindGetVector(fPDGCode); + const py::array_t &fPDGCode) { + this->fPDGCode = PyBindGetVector(fPDGCode); } void GatePhaseSpaceSource::SetEnergyBatch( - const py::array_t &fEnergy) const { - auto &ll = fThreadLocalDataPhsp.Get(); - ll.fEnergy = PyBindGetVector(fEnergy); + const py::array_t &fEnergy) { + this->fEnergy = PyBindGetVector(fEnergy); } void GatePhaseSpaceSource::SetWeightBatch( - const py::array_t &fWeight) const { - auto &ll = fThreadLocalDataPhsp.Get(); - ll.fWeight = PyBindGetVector(fWeight); + const py::array_t &fWeight) { + this->fWeight = PyBindGetVector(fWeight); } void GatePhaseSpaceSource::SetPositionXBatch( - const py::array_t &fPositionX) const { - auto &ll = fThreadLocalDataPhsp.Get(); - ll.fPositionX = PyBindGetVector(fPositionX); + const py::array_t &fPositionX) { + this->fPositionX = PyBindGetVector(fPositionX); } void GatePhaseSpaceSource::SetPositionYBatch( - const py::array_t &fPositionY) const { - auto &ll = fThreadLocalDataPhsp.Get(); - ll.fPositionY = PyBindGetVector(fPositionY); + const py::array_t &fPositionY) { + this->fPositionY = PyBindGetVector(fPositionY); } void GatePhaseSpaceSource::SetPositionZBatch( - const py::array_t &fPositionZ) const { - auto &ll = fThreadLocalDataPhsp.Get(); - ll.fPositionZ = PyBindGetVector(fPositionZ); + const py::array_t &fPositionZ) { + this->fPositionZ = PyBindGetVector(fPositionZ); } void GatePhaseSpaceSource::SetDirectionXBatch( - const py::array_t &fDirectionX) const { - auto &ll = fThreadLocalDataPhsp.Get(); - ll.fDirectionX = PyBindGetVector(fDirectionX); + const py::array_t &fDirectionX) { + this->fDirectionX = PyBindGetVector(fDirectionX); } void GatePhaseSpaceSource::SetDirectionYBatch( - const py::array_t &fDirectionY) const { - auto &ll = fThreadLocalDataPhsp.Get(); - ll.fDirectionY = PyBindGetVector(fDirectionY); + const py::array_t &fDirectionY) { + this->fDirectionY = PyBindGetVector(fDirectionY); } void GatePhaseSpaceSource::SetDirectionZBatch( - const py::array_t &fDirectionZ) const { - auto &ll = fThreadLocalDataPhsp.Get(); - ll.fDirectionZ = PyBindGetVector(fDirectionZ); + const py::array_t &fDirectionZ) { + this->fDirectionZ = PyBindGetVector(fDirectionZ); } bool GatePhaseSpaceSource::ParticleIsPrimary() const { - auto &ll = fThreadLocalDataPhsp.Get(); // check if particle is primary bool is_primary = false; // if PDGCode exists in file - if ((ll.fPDGCode[ll.fCurrentIndex] != 0) && (ll.fPrimaryPDGCode != 0)) { - if ((ll.fPrimaryPDGCode == ll.fPDGCode[ll.fCurrentIndex]) && - (ll.fPrimaryLowerEnergyThreshold <= ll.fEnergy[ll.fCurrentIndex])) { + if ((fPDGCode[fCurrentIndex] != 0) && (fPrimaryPDGCode != 0)) { + if ((fPrimaryPDGCode == fPDGCode[fCurrentIndex]) && + (fPrimaryLowerEnergyThreshold <= fEnergy[fCurrentIndex])) { is_primary = true; } } else { diff --git a/core/opengate_core/opengate_lib/GatePhaseSpaceSource.h b/core/opengate_core/opengate_lib/GatePhaseSpaceSource.h index 42cf52f80..47466d53e 100644 --- a/core/opengate_core/opengate_lib/GatePhaseSpaceSource.h +++ b/core/opengate_core/opengate_lib/GatePhaseSpaceSource.h @@ -41,7 +41,7 @@ class GatePhaseSpaceSource : public GateVSource { const G4ThreeVector &momentum_direction, double energy, double time, double w); - void SetGeneratorFunction(ParticleGeneratorType &f) const; + void SetGeneratorFunction(ParticleGeneratorType &f); bool ParticleIsPrimary() const; @@ -57,52 +57,47 @@ class GatePhaseSpaceSource : public GateVSource { bool fVerbose; G4bool fIsotropicMomentum; - void SetPDGCodeBatch(const py::array_t &fPDGCode) const; + void SetPDGCodeBatch(const py::array_t &fPDGCode); - void SetEnergyBatch(const py::array_t &fEnergy) const; + void SetEnergyBatch(const py::array_t &fEnergy); - void SetWeightBatch(const py::array_t &fWeight) const; + void SetWeightBatch(const py::array_t &fWeight); - void SetPositionXBatch(const py::array_t &fPositionX) const; + void SetPositionXBatch(const py::array_t &fPositionX); - void SetPositionYBatch(const py::array_t &fPositionY) const; + void SetPositionYBatch(const py::array_t &fPositionY); - void SetPositionZBatch(const py::array_t &fPositionZ) const; + void SetPositionZBatch(const py::array_t &fPositionZ); - void SetDirectionXBatch(const py::array_t &fDirectionX) const; + void SetDirectionXBatch(const py::array_t &fDirectionX); - void SetDirectionYBatch(const py::array_t &fDirectionY) const; + void SetDirectionYBatch(const py::array_t &fDirectionY); - void SetDirectionZBatch(const py::array_t &fDirectionZ) const; + void SetDirectionZBatch(const py::array_t &fDirectionZ); - // For MT, all threads local variables are gathered here - struct threadLocalTPhsp { - G4ParticleDefinition *fParticleDefinition; +protected: + G4ParticleDefinition *fParticleDefinition = nullptr; - bool fGenerateUntilNextPrimary; - std::int32_t fPrimaryPDGCode; - std::float_t fPrimaryLowerEnergyThreshold; + bool fGenerateUntilNextPrimary = false; + std::int32_t fPrimaryPDGCode = 0; + std::float_t fPrimaryLowerEnergyThreshold = 0.0; - ParticleGeneratorType fGenerator; - unsigned long fNumberOfGeneratedEvents; - size_t fCurrentIndex; - size_t fCurrentBatchSize; + ParticleGeneratorType fGenerator; + size_t fCurrentIndex = 0; + size_t fCurrentBatchSize = 0; - std::int32_t *fPDGCode; + std::int32_t *fPDGCode = nullptr; - std::float_t *fPositionX; - std::float_t *fPositionY; - std::float_t *fPositionZ; + std::float_t *fPositionX = nullptr; + std::float_t *fPositionY = nullptr; + std::float_t *fPositionZ = nullptr; - std::float_t *fDirectionX; - std::float_t *fDirectionY; - std::float_t *fDirectionZ; + std::float_t *fDirectionX = nullptr; + std::float_t *fDirectionY = nullptr; + std::float_t *fDirectionZ = nullptr; - std::float_t *fEnergy; - std::float_t *fWeight; - // double * fTime; // FIXME todo - }; - G4Cache fThreadLocalDataPhsp; + std::float_t *fEnergy = nullptr; + std::float_t *fWeight = nullptr; }; #endif // GatePhaseSpaceSource_h diff --git a/core/opengate_core/opengate_lib/GateSourceManager.cpp b/core/opengate_core/opengate_lib/GateSourceManager.cpp index 40fb46ca0..38d7bd59b 100644 --- a/core/opengate_core/opengate_lib/GateSourceManager.cpp +++ b/core/opengate_core/opengate_lib/GateSourceManager.cpp @@ -284,7 +284,6 @@ void GateSourceManager::CheckForNextRun() const { G4RunManager::GetRunManager()->AbortRun(true); // FIXME true or false ? l.fStartNewRun = true; l.fNextRunId++; - /* if (l.fNextRunId >= fSimulationTimes.size()) { // Sometimes, the source must clean some data in its own thread, not by // the master thread (for example, with a G4SingleParticleSource object) @@ -293,7 +292,6 @@ void GateSourceManager::CheckForNextRun() const { source->CleanWorkerThread(); } } - */ } } diff --git a/core/opengate_core/opengate_lib/GateTemplateSource.cpp b/core/opengate_core/opengate_lib/GateTemplateSource.cpp index a715085fa..aab5bcee6 100644 --- a/core/opengate_core/opengate_lib/GateTemplateSource.cpp +++ b/core/opengate_core/opengate_lib/GateTemplateSource.cpp @@ -59,8 +59,7 @@ void GateTemplateSource::GeneratePrimaries(G4Event *event, // the position is changed according to fGlobalTranslation and fGlobalRotation auto pos = G4ThreeVector(fVectorValue[0], fVectorValue[1], fVectorValue[2]); - auto &l = fThreadLocalData.Get(); - pos = l.fGlobalRotation * pos + l.fGlobalTranslation; + pos = fGlobalRotation * pos + fGlobalTranslation; // create the vertex auto *vertex = new G4PrimaryVertex(pos, current_simulation_time); diff --git a/core/opengate_core/opengate_lib/GateTreatmentPlanPBSource.cpp b/core/opengate_core/opengate_lib/GateTreatmentPlanPBSource.cpp index 190a5c6e9..3a07c8e91 100644 --- a/core/opengate_core/opengate_lib/GateTreatmentPlanPBSource.cpp +++ b/core/opengate_core/opengate_lib/GateTreatmentPlanPBSource.cpp @@ -25,22 +25,21 @@ GateTreatmentPlanPBSource::GateTreatmentPlanPBSource() : GateVSource() { fSortedSpotGenerationFlag = false; fPDF = nullptr; fTotalNumberOfSpots = 0; + + fSPS_PB = nullptr; + fCurrentSpot = 0; + fPreviousSpot = -1; + fInitGenericIon = false; } GateTreatmentPlanPBSource::~GateTreatmentPlanPBSource() = default; -GateTreatmentPlanPBSource::threadLocalTPSource & -GateTreatmentPlanPBSource::GetThreadLocalDataTPSource() { - return fThreadLocalDataTPSource.Get(); -} - void GateTreatmentPlanPBSource::InitializeUserInfo(py::dict &user_info) { GateVSource::InitializeUserInfo(user_info); - auto &ll = GetThreadLocalDataTPSource(); // Create single particle source only once. Parameters are then updated for // the different spots. - ll.fSPS_PB = new GateSingleParticleSourcePencilBeam(std::string(), - fAttachedToVolumeName); + fSPS_PB = new GateSingleParticleSourcePencilBeam(std::string(), + fAttachedToVolumeName); // common to all spots InitializeParticle(user_info); @@ -58,9 +57,9 @@ void GateTreatmentPlanPBSource::InitializeUserInfo(py::dict &user_info) { fSpotPosition = DictGetVecG4ThreeVector(user_info, "positions"); fSpotRotation = DictGetVecG4RotationMatrix(user_info, "rotations"); - fTotalNumberOfSpots = static_cast(fSpotWeight.size()); - ll.fNbGeneratedSpots.resize(fTotalNumberOfSpots, - 0); // keep track for debug + fTotalNumberOfSpots = fSpotWeight.size(); + fNbGeneratedSpots.resize(fTotalNumberOfSpots, + 0); // keep track for debug // Init the random fEngine InitRandomEngine(); @@ -71,12 +70,11 @@ void GateTreatmentPlanPBSource::InitializeUserInfo(py::dict &user_info) { } void GateTreatmentPlanPBSource::InitNbPrimariesVec() { - auto &ll = GetThreadLocalDataTPSource(); // Initialize all spots to zero particles - ll.fNbIonsToGenerate.resize(fTotalNumberOfSpots, 0); - for (unsigned long i = 0; i < fMaxN; i++) { - int bin = static_cast(fTotalNumberOfSpots * fDistriGeneral->fire()); - ++ll.fNbIonsToGenerate[bin]; + fNbIonsToGenerate.resize(fTotalNumberOfSpots, 0); + for (long int i = 0; i < fMaxN; i++) { + int bin = fTotalNumberOfSpots * fDistriGeneral->fire(); + ++fNbIonsToGenerate[bin]; } } void GateTreatmentPlanPBSource::InitRandomEngine() { @@ -121,117 +119,106 @@ void GateTreatmentPlanPBSource::PrepareNextRun() { void GateTreatmentPlanPBSource::GeneratePrimaries( G4Event *event, double current_simulation_time) { - auto &ll = GetThreadLocalDataTPSource(); // Find next spot to initialize FindNextSpot(); // if we moved to a new spot, we need to update the SPS parameters - if (ll.fCurrentSpot != ll.fPreviousSpot) { + if (fCurrentSpot != fPreviousSpot) { ConfigureSingleSpot(); } // Generate vertex - ll.fSPS_PB->SetParticleTime(current_simulation_time); - ll.fSPS_PB->GeneratePrimaryVertex(event); + fSPS_PB->SetParticleTime(current_simulation_time); + fSPS_PB->GeneratePrimaryVertex(event); // weight - double w = fSpotWeight[ll.fCurrentSpot]; + double w = fSpotWeight[fCurrentSpot]; for (auto i = 0; i < event->GetNumberOfPrimaryVertex(); i++) { event->GetPrimaryVertex(i)->SetWeight(w); } // update number of generated events // fNumberOfGeneratedEvents++; - ll.fNbGeneratedSpots[ll.fCurrentSpot]++; + fNbGeneratedSpots[fCurrentSpot]++; if (fSortedSpotGenerationFlag) { // we generated an ion from this spot, so we remove it from the vector - --ll.fNbIonsToGenerate[ll.fCurrentSpot]; + --fNbIonsToGenerate[fCurrentSpot]; } // update previous spot - ll.fPreviousSpot = ll.fCurrentSpot; + fPreviousSpot = fCurrentSpot; } void GateTreatmentPlanPBSource::FindNextSpot() { - auto &ll = GetThreadLocalDataTPSource(); if (fSortedSpotGenerationFlag) { // move to next spot if there are no more particles to generate in the // current one - while ((ll.fCurrentSpot < fTotalNumberOfSpots) && - (ll.fNbIonsToGenerate[ll.fCurrentSpot] <= 0)) { - ll.fCurrentSpot++; + while ((fCurrentSpot < fTotalNumberOfSpots) && + (fNbIonsToGenerate[fCurrentSpot] <= 0)) { + fCurrentSpot++; } } else { // select random spot according to PDF - int bin = static_cast(fTotalNumberOfSpots * fDistriGeneral->fire()); - ll.fCurrentSpot = bin; + int bin = fTotalNumberOfSpots * fDistriGeneral->fire(); + fCurrentSpot = bin; } } void GateTreatmentPlanPBSource::ConfigureSingleSpot() { - auto &ll = GetThreadLocalDataTPSource(); // Particle definition if ion - if (ll.fInitGenericIon) { + if (fInitGenericIon) { auto *ion_table = G4IonTable::GetIonTable(); auto *ion = ion_table->GetIon(fZ, fA, fE); - ll.fSPS_PB->SetParticleDefinition(ion); - ll.fInitGenericIon = false; // only the first time + fSPS_PB->SetParticleDefinition(ion); + fInitGenericIon = false; // only the first time } // Energy - double energy = fSpotEnergy[ll.fCurrentSpot]; - double sigmaE = fSigmaEnergy[ll.fCurrentSpot]; + double energy = fSpotEnergy[fCurrentSpot]; + double sigmaE = fSigmaEnergy[fCurrentSpot]; UpdateEnergySPS(energy, sigmaE); // rotation and translation - G4ThreeVector translation = fSpotPosition[ll.fCurrentSpot]; - G4RotationMatrix rotation = fSpotRotation[ll.fCurrentSpot]; + G4ThreeVector translation = fSpotPosition[fCurrentSpot]; + G4RotationMatrix rotation = fSpotRotation[fCurrentSpot]; UpdatePositionSPS(translation, rotation); // Phase space parameters - std::vector x_param = fPhSpaceX[ll.fCurrentSpot]; - std::vector y_param = fPhSpaceY[ll.fCurrentSpot]; - ll.fSPS_PB->SetPBSourceParam(x_param, y_param); + std::vector x_param = fPhSpaceX[fCurrentSpot]; + std::vector y_param = fPhSpaceY[fCurrentSpot]; + fSPS_PB->SetPBSourceParam(x_param, y_param); } void GateTreatmentPlanPBSource::UpdatePositionSPS( const G4ThreeVector &localTransl, const G4RotationMatrix &localRot) { // update local translation and rotation - auto &ll = GetThreadLocalDataTPSource(); - auto &l = fThreadLocalData.Get(); - - l.fGlobalTranslation = localTransl; - l.fGlobalRotation = localRot; // ConvertToG4RotationMatrix(localRot); - - // // update global rotation - // GateVSource::SetOrientationAccordingToAttachedVolume(); + fGlobalTranslation = localTransl; + fGlobalRotation = localRot; // ConvertToG4RotationMatrix(localRot); // No change in the translation rotation if mother is the world if (fAttachedToVolumeName == "world") { // set it to the vertex - ll.fSPS_PB->SetSourceRotTransl(l.fGlobalTranslation, l.fGlobalRotation); + fSPS_PB->SetSourceRotTransl(fGlobalTranslation, fGlobalRotation); return; } // compute global translation rotation. - // l.fGlobalTranslation and l.fGlobalRotation values are updated here. + // fGlobalTranslation and fGlobalRotation values are updated here. ComputeTransformationFromVolumeToWorld( - fAttachedToVolumeName, l.fGlobalTranslation, l.fGlobalRotation, false); + fAttachedToVolumeName, fGlobalTranslation, fGlobalRotation, false); // set it to the vertex - ll.fSPS_PB->SetSourceRotTransl(l.fGlobalTranslation, l.fGlobalRotation); + fSPS_PB->SetSourceRotTransl(fGlobalTranslation, fGlobalRotation); } void GateTreatmentPlanPBSource::UpdateEnergySPS(double energy, double sigma) { - auto &ll = GetThreadLocalDataTPSource(); - auto *ene = ll.fSPS_PB->GetEneDist(); + auto *ene = fSPS_PB->GetEneDist(); ene->SetEnergyDisType("Gauss"); ene->SetMonoEnergy(energy); ene->SetBeamSigmaInE(sigma); } void GateTreatmentPlanPBSource::InitializeParticle(py::dict &user_info) { - auto &ll = GetThreadLocalDataTPSource(); std::string pname = DictGetStr(user_info, "particle"); // If the particle is an ion (name start with ion) if (pname.rfind("ion", 0) == 0) { @@ -245,23 +232,22 @@ void GateTreatmentPlanPBSource::InitializeParticle(py::dict &user_info) { if (fParticleDefinition == nullptr) { Fatal("Cannot find the particle '" + pname + "'."); } - ll.fSPS_PB->SetParticleDefinition(fParticleDefinition); + fSPS_PB->SetParticleDefinition(fParticleDefinition); } void GateTreatmentPlanPBSource::InitializeIon(py::dict &user_info) { - auto &ll = GetThreadLocalDataTPSource(); auto u = py::dict(user_info["ion"]); fA = DictGetInt(u, "A"); fZ = DictGetInt(u, "Z"); fE = DictGetDouble(u, "E"); - ll.fInitGenericIon = true; + fInitGenericIon = true; } py::list GateTreatmentPlanPBSource::GetGeneratedPrimaries() { py::list n_spot_vec; - // for (const auto &item : fNbGeneratedSpots) { - // n_spot_vec.append(item); - // } + for (const auto &item : fNbGeneratedSpots) { + n_spot_vec.append(item); + } return n_spot_vec; } diff --git a/core/opengate_core/opengate_lib/GateTreatmentPlanPBSource.h b/core/opengate_core/opengate_lib/GateTreatmentPlanPBSource.h index dd10eb88f..ceee46dc8 100644 --- a/core/opengate_core/opengate_lib/GateTreatmentPlanPBSource.h +++ b/core/opengate_core/opengate_lib/GateTreatmentPlanPBSource.h @@ -34,18 +34,12 @@ class GateTreatmentPlanPBSource : public GateVSource { py::list GetGeneratedPrimaries(); protected: - // thread local structure - struct threadLocalTPSource { - GateSingleParticleSourcePencilBeam *fSPS_PB = nullptr; - std::vector fNbIonsToGenerate; - std::vector fNbGeneratedSpots; - int fCurrentSpot = 0; - int fPreviousSpot = -1; - bool fInitGenericIon = false; - }; - G4Cache fThreadLocalDataTPSource; - - threadLocalTPSource &GetThreadLocalDataTPSource(); + GateSingleParticleSourcePencilBeam *fSPS_PB = nullptr; + std::vector fNbIonsToGenerate; + std::vector fNbGeneratedSpots; + int fCurrentSpot = 0; + int fPreviousSpot = -1; + bool fInitGenericIon = false; // variables common to all spots CLHEP::HepRandomEngine *fEngine; diff --git a/core/opengate_core/opengate_lib/GateVSource.cpp b/core/opengate_core/opengate_lib/GateVSource.cpp index 5e98a62b2..55d5ee488 100644 --- a/core/opengate_core/opengate_lib/GateVSource.cpp +++ b/core/opengate_core/opengate_lib/GateVSource.cpp @@ -27,10 +27,6 @@ GateVSource::GateVSource() { GateVSource::~GateVSource() = default; -GateVSource::threadLocalT &GateVSource::GetThreadLocalData() { - return fThreadLocalData.Get(); -} - void GateVSource::InitializeUserInfo(py::dict &user_info) { // get info from the dict fName = DictGetStr(user_info, "name"); @@ -66,10 +62,9 @@ double GateVSource::CalcNextTime(double current_simulation_time) { } void GateVSource::PrepareNextRun() { - auto &l = GetThreadLocalData(); - l.fNumberOfGeneratedEvents = 0; - fMaxN = fVectorOfMaxN[l.fRunID]; - l.fRunID++; + fNumberOfGeneratedEvents = 0; + fMaxN = fVectorOfMaxN[fRunID]; + fRunID++; SetOrientationAccordingToAttachedVolume(); } @@ -97,9 +92,8 @@ void GateVSource::GeneratePrimaries(G4Event * /*event*/, double /*time*/) { } void GateVSource::SetOrientationAccordingToAttachedVolume() { - auto &l = GetThreadLocalData(); - l.fGlobalRotation = fLocalRotation; - l.fGlobalTranslation = fLocalTranslation; + fGlobalRotation = fLocalRotation; + fGlobalTranslation = fLocalTranslation; // No change in the translation rotation if mother is the world if (fAttachedToVolumeName == "world") @@ -108,7 +102,7 @@ void GateVSource::SetOrientationAccordingToAttachedVolume() { // compute global translation rotation and keep it. // Will be used, for example, in GenericSource to change position ComputeTransformationFromVolumeToWorld( - fAttachedToVolumeName, l.fGlobalTranslation, l.fGlobalRotation, false); + fAttachedToVolumeName, fGlobalTranslation, fGlobalRotation, false); } unsigned long diff --git a/core/opengate_core/opengate_lib/GateVSource.h b/core/opengate_core/opengate_lib/GateVSource.h index 687b40030..0f6996c69 100644 --- a/core/opengate_core/opengate_lib/GateVSource.h +++ b/core/opengate_core/opengate_lib/GateVSource.h @@ -52,10 +52,7 @@ class GateVSource { virtual unsigned long GetExpectedNumberOfEvents(const TimeInterval &time_interval); - G4int GetNumberOfSimulatedEvents() { - auto &l = fThreadLocalData.Get(); - return l.fNumberOfGeneratedEvents; - } + G4int GetNumberOfSimulatedEvents() { return fNumberOfGeneratedEvents; } std::vector GetVectorOfSimulatedEvents() { return fVectorOfMaxN; } @@ -82,15 +79,8 @@ class GateVSource { double fHalfLife; double fDecayConstant; - struct threadLocalT { - unsigned long fNumberOfGeneratedEvents = 0; - G4ThreeVector fGlobalTranslation; - G4RotationMatrix fGlobalRotation; - G4int fRunID = 0; - }; - G4Cache fThreadLocalData; - - virtual threadLocalT &GetThreadLocalData(); + unsigned long fNumberOfGeneratedEvents = 0; + G4int fRunID = 0; }; #endif // GateVSource_h diff --git a/core/opengate_core/opengate_lib/GateVoxelSource.cpp b/core/opengate_core/opengate_lib/GateVoxelSource.cpp index 7070253dc..cf78d513b 100644 --- a/core/opengate_core/opengate_lib/GateVoxelSource.cpp +++ b/core/opengate_core/opengate_lib/GateVoxelSource.cpp @@ -14,32 +14,28 @@ GateVoxelSource::GateVoxelSource() : GateGenericSource() { GateVoxelSource::~GateVoxelSource() = default; void GateVoxelSource::PrepareNextRun() { - // GateGenericSource::PrepareNextRun(); - // rotation and translation to apply, according to mother volume + // The following compute the global transformation from + // the local volume (mother) to the world GateVSource::PrepareNextRun(); // This global transformation is given to the SPS that will // generate particles in the correct coordinate system - auto &l = GetThreadLocalData(); - auto &ll = GetThreadLocalDataGenericSource(); - auto *pos = ll.fSPS->GetPosDist(); - pos->SetCentreCoords(l.fGlobalTranslation); + auto *pos = fSPS->GetPosDist(); + pos->SetCentreCoords(fGlobalTranslation); // orientation according to mother volume - auto rotation = l.fGlobalRotation; + const auto rotation = fGlobalRotation; G4ThreeVector r1(rotation(0, 0), rotation(1, 0), rotation(2, 0)); G4ThreeVector r2(rotation(0, 1), rotation(1, 1), rotation(2, 1)); pos->SetPosRot1(r1); pos->SetPosRot2(r2); - // auto &l = fThreadLocalData.Get(); - fVoxelPositionGenerator->fGlobalRotation = l.fGlobalRotation; - fVoxelPositionGenerator->fGlobalTranslation = l.fGlobalTranslation; + fVoxelPositionGenerator->fGlobalRotation = fGlobalRotation; + fVoxelPositionGenerator->fGlobalTranslation = fGlobalTranslation; // the direction is 'isotropic' so we don't care about rotating the direction. } void GateVoxelSource::InitializePosition(py::dict) { - auto &ll = GetThreadLocalDataGenericSource(); - ll.fSPS->SetPosGenerator(fVoxelPositionGenerator); + fSPS->SetPosGenerator(fVoxelPositionGenerator); // we set a fake value (not used) fVoxelPositionGenerator->SetPosDisType("Point"); } diff --git a/core/opengate_core/opengate_lib/GateVoxelizedPromptGammaTLESource.cpp b/core/opengate_core/opengate_lib/GateVoxelizedPromptGammaTLESource.cpp index 999492d9f..e78677e67 100644 --- a/core/opengate_core/opengate_lib/GateVoxelizedPromptGammaTLESource.cpp +++ b/core/opengate_core/opengate_lib/GateVoxelizedPromptGammaTLESource.cpp @@ -23,27 +23,23 @@ void GateVoxelizedPromptGammaTLESource::PrepareNextRun() { GateVSource::PrepareNextRun(); // This global transformation is given to the SPS that will // generate particles in the correct coordinate system - auto &l = GetThreadLocalData(); - auto &ll = GetThreadLocalDataGenericSource(); - auto *pos = ll.fSPS->GetPosDist(); - pos->SetCentreCoords(l.fGlobalTranslation); + auto *pos = fSPS->GetPosDist(); + pos->SetCentreCoords(fGlobalTranslation); // orientation according to mother volume - auto rotation = l.fGlobalRotation; + auto rotation = fGlobalRotation; G4ThreeVector r1(rotation(0, 0), rotation(1, 0), rotation(2, 0)); G4ThreeVector r2(rotation(0, 1), rotation(1, 1), rotation(2, 1)); pos->SetPosRot1(r1); pos->SetPosRot2(r2); - // auto &l = fThreadLocalData.Get(); - fVoxelPositionGenerator->fGlobalRotation = l.fGlobalRotation; - fVoxelPositionGenerator->fGlobalTranslation = l.fGlobalTranslation; + fVoxelPositionGenerator->fGlobalRotation = fGlobalRotation; + fVoxelPositionGenerator->fGlobalTranslation = fGlobalTranslation; // the direction is 'isotropic' so we don't care about rotating the direction. } void GateVoxelizedPromptGammaTLESource::InitializePosition(py::dict) { - auto &ll = GetThreadLocalDataGenericSource(); - ll.fSPS->SetPosGenerator(fVoxelPositionGenerator); + fSPS->SetPosGenerator(fVoxelPositionGenerator); // we set a fake value (not used) fVoxelPositionGenerator->SetPosDisType("Point"); } diff --git a/core/opengate_core/opengate_lib/pyGateDebugActor.cpp b/core/opengate_core/opengate_lib/pyGateDebugActor.cpp new file mode 100644 index 000000000..12655fa47 --- /dev/null +++ b/core/opengate_core/opengate_lib/pyGateDebugActor.cpp @@ -0,0 +1,20 @@ +/* -------------------------------------------------- + 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 "GateDebugActor.h" + +void init_GateDebugActor(py::module &m) { + py::class_, + GateVActor>(m, "GateDebugActor") + .def(py::init()) + .def("InitializeUserInfo", &GateDebugActor::InitializeUserInfo); +} diff --git a/core/opengate_core/opengate_lib/pyGateDebugSource.cpp b/core/opengate_core/opengate_lib/pyGateDebugSource.cpp new file mode 100644 index 000000000..76365e02f --- /dev/null +++ b/core/opengate_core/opengate_lib/pyGateDebugSource.cpp @@ -0,0 +1,17 @@ +/* -------------------------------------------------- + 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 "GateDebugSource.h" + +void init_GateDebugSource(py::module &m) { + + py::class_(m, "GateDebugSource") + .def(py::init()) + .def("GetDebugValue", &GateDebugSource::GetDebugValue); +} diff --git a/docs/source/user_guide/user_guide_chemistry.rst b/docs/source/user_guide/user_guide_chemistry.rst index e68bb9454..238093710 100644 --- a/docs/source/user_guide/user_guide_chemistry.rst +++ b/docs/source/user_guide/user_guide_chemistry.rst @@ -282,7 +282,7 @@ Adding a chemistry actor ------------------------ The current reference chemistry actor is ``ChemicalCountingActor``. -Further chemistry actors can be developed and made available in the future. +Further chemistry actors can be developed and made available in the future. It is intended as a passive chemistry-scoring actor. Simulation-wide chemistry control stays on ``ChemistryManager``. diff --git a/opengate/actors/miscactors.py b/opengate/actors/miscactors.py index e2411a1c6..9ddde9ece 100644 --- a/opengate/actors/miscactors.py +++ b/opengate/actors/miscactors.py @@ -636,6 +636,63 @@ def BeginOfRunAction(self, run): self.user_output.attenuation_image.end_of_simulation() +class DebugActor(ActorBase, g4.GateDebugActor): + """ + Process tracking for debugging and education purposes. + + Example usage in Python: + debug = sim.add_actor("DebugActor", "debug") + debug.debug_flag = True + """ + + user_info_defaults = {"debug_flag": (False, {"doc": "Test option"})} + + def __init__(self, *args, **kwargs): + print(f"(python) DebugActor: __init__") + ActorBase.__init__(self, *args, **kwargs) + self.__initcpp__() + + def __initcpp__(self): + print(f"(python) DebugActor ({self.name}) : __initcpp__") + g4.GateDebugActor.__init__(self, self.user_info) + print(f"(python) DebugActor ({self.name}) : AddActions") + self.AddActions( + { + "BeginOfSimulationAction", + "BeginOfRunAction", + "PreUserTrackingAction", + "PostUserTrackingAction", + "BeginOfEventAction", + "EndOfEventAction", + "SteppingAction", + "EndOfRunAction", + "EndOfSimulationAction", + } + ) + + def __getstate__(self): + print(f"(python) DebugActor ({self.name}) : __getstate__") + return ActorBase.__getstate__(self) + + def __setstate__(self, state): + print(f"(python) DebugActor ({self.name}) : __setstate__") + ActorBase.__setstate__(self, state) + + def initialize(self): + print(f"(python) DebugActor ({self.name}) : initialize") + ActorBase.initialize(self) + self.InitializeUserInfo(self.user_info) + self.InitializeCpp() + + def BeginOfSimulationAction(self): + print(f"(python) DebugActor ({self.name}) : BeginOfSimulationAction") + g4.GateDebugActor.BeginOfSimulationAction(self) + + def EndOfSimulationAction(self): + print(f"(python) DebugActor ({self.name}) : EndOfSimulationAction") + g4.GateDebugActor.EndOfSimulationAction(self) + + process_cls(ActorOutputStatisticsActor) process_cls(SimulationStatisticsActor) process_cls(KillActor) @@ -646,3 +703,4 @@ def BeginOfRunAction(self, run): process_cls(ActorOutputKillNonInteractingParticleActor) process_cls(KillNonInteractingParticleActor) process_cls(AttenuationImageActor) +process_cls(DebugActor) diff --git a/opengate/engines.py b/opengate/engines.py index 4f9f67730..1434cf7ee 100644 --- a/opengate/engines.py +++ b/opengate/engines.py @@ -87,7 +87,7 @@ def __init__(self, simulation_engine): self.progress_bar = False self.expected_number_of_events = "unknown" - # List of sources (GateSource), for all threads + # List of sources (GateSource) self.sources = [] # The source manager will be constructed at build (during ActionManager) @@ -109,6 +109,8 @@ def source_manager(self): def close(self): if self.verbose_close: warning("Closing SourceEngine") + for source in self.source_manager.sources.values(): + source.close() self.release_g4_references() super().close() @@ -128,6 +130,11 @@ def initialize(self, run_timing_intervals, progress_bar=False): self.progress_bar = progress_bar self.initialize_dynamic_parametrisations() + # Pre-create C++ sources for all threads (1 master + N workers) on the main thread + num_instances = 1 + self.simulation_engine.simulation.number_of_threads + for source in self.source_manager.sources.values(): + source.pre_create_g4_sources(num_instances) + def initialize_actors(self): """ Set all actors to the master source manager @@ -153,6 +160,12 @@ def create_master_source_manager(self): self.g4_master_source_manager = self.create_g4_thread_source_manager( append=False ) + # Initialize all worker thread-local C++ sources on the main thread to size G4Cache properly. + # This is safe because physics and geometry are fully initialized now. + for source in self.source_manager.sources.values(): + if source.g4_thread_sources: + for g4_src in source.g4_thread_sources[1:]: + source.initialize_g4_source(g4_src, self.run_timing_intervals) return self.g4_master_source_manager def create_g4_thread_source_manager(self, append=True): @@ -166,8 +179,13 @@ def create_g4_thread_source_manager(self, append=True): # create all sources for this source manager (for all threads) source_manager = self.simulation_engine.simulation.source_manager for source in source_manager.sources.values(): - source.initialize(self.run_timing_intervals) - source.add_to_source_manager(ms) + g4_source = source.get_next_g4_source() + if g4_source is not None: + source.initialize_g4_source(g4_source, self.run_timing_intervals) + ms.AddSource(g4_source) + else: + source.initialize(self.run_timing_intervals) + source.add_to_source_manager(ms) # store all the sources (will be used later by SimulationOutput) self.sources.append(source) @@ -217,6 +235,11 @@ def start(self): # once terminated, packup the sources (if needed) for source in self.sources: source.prepare_output() + if source.g4_thread_sources: + source.gather_outputs(source.g4_thread_sources) + # Clear C++ sources now on the main thread while G4Cache is active + source.g4_thread_sources = [] + source.g4_thread_sources_index = 0 def can_predict_expected_number_of_event(self): # can_predict = True @@ -1499,6 +1522,9 @@ def close(self): self.g4_RunManager.SetVerboseLevel(0) self._is_closed = True self.g4_RunManager = None + if self.run_manager_finalizer: + self.run_manager_finalizer.detach() + self.run_manager_finalizer = None def __enter__(self): return self diff --git a/opengate/managers.py b/opengate/managers.py index 7dc7aa65b..48b6128c8 100644 --- a/opengate/managers.py +++ b/opengate/managers.py @@ -1,3 +1,4 @@ +from PIL import DcxImagePlugin import copy import io import shutil @@ -34,6 +35,7 @@ ) from .processing import dispatch_to_subprocess from .serialization import dump_json, dumps_json, load_json, loads_json +from .sources.base import DebugSource from .sources.beamsources import IonPencilBeamSource, TreatmentPlanPBSource from .sources.gansources import GANPairsSource, GANSource from .sources.generic import GenericSource, SourceBase @@ -52,6 +54,7 @@ source_types = { "GenericSource": GenericSource, + "DebugSource": DebugSource, "LastVertexSource": LastVertexSource, "PhaseSpaceSource": PhaseSpaceSource, "VoxelSource": VoxelSource, @@ -110,6 +113,7 @@ KillActor, KillNonInteractingParticleActor, SimulationStatisticsActor, + DebugActor, ) from .actors.pgactors import ( VoxelizedPromptGammaAnalogActor, @@ -169,6 +173,7 @@ "ChemicalCountingActor": ChemicalCountingActor, "ARFActor": ARFActor, "ARFTrainingDatasetActor": ARFTrainingDatasetActor, + "DebugActor": DebugActor, # digit "PhaseSpaceActor": PhaseSpaceActor, "DigitizerAdderActor": DigitizerAdderActor, @@ -376,6 +381,18 @@ def initialize_before_g4_engine(self): if source.initialize_source_before_g4_engine: source.initialize_source_before_g4_engine(source) + def to_dictionary(self): + d = super().to_dictionary() + d["sources"] = dict([(k, v.to_dictionary()) for k, v in self.sources.items()]) + return d + + def from_dictionary(self, d): + self.sources = {} + super().from_dictionary(d) + for k, v in d["sources"].items(): + s = self.add_source(v["object_type"], name=v["user_info"]["name"]) + s.from_dictionary(v) + class ActorManager(GateObject): """ @@ -1929,6 +1946,7 @@ def to_dictionary(self): d["physics_manager"] = self.physics_manager.to_dictionary() d["chemistry_manager"] = self.chemistry_manager.to_dictionary() d["actor_manager"] = self.actor_manager.to_dictionary() + d["source_manager"] = self.source_manager.to_dictionary() d["auxiliary_attributes"] = dict( [(k, v.to_dictionary()) for k, v in self.auxiliary_attributes.items()] ) @@ -1941,6 +1959,8 @@ def from_dictionary(self, d): if "chemistry_manager" in d: self.chemistry_manager.from_dictionary(d["chemistry_manager"]) self.actor_manager.from_dictionary(d["actor_manager"]) + if "source_manager" in d: + self.source_manager.from_dictionary(d["source_manager"]) self.auxiliary_attributes = {} for _, v in d.get("auxiliary_attributes", {}).items(): a = self.activate_auxiliary_attribute( @@ -1949,11 +1969,9 @@ def from_dictionary(self, d): a.from_dictionary(v) def to_json_string(self): - warning("Only parts of the simulation can currently be dumped as JSON") return dumps_json(self.to_dictionary()) def to_json_file(self, directory=None, filename=None): - warning("Only parts of the simulation can currently be dumped as JSON.") d = self.to_dictionary() if filename is None: filename = self.json_archive_filename @@ -1965,11 +1983,9 @@ def to_json_file(self, directory=None, filename=None): self.copy_input_files(directory, dct=d) def from_json_string(self, json_string): - warning("Only parts of the simulation can currently be reloaded from JSON.") self.from_dictionary(loads_json(json_string)) def from_json_file(self, path): - warning("Only parts of the simulation can currently be reloaded from JSON.") with open(path, "r") as f: self.from_dictionary(load_json(f)) @@ -2209,7 +2225,6 @@ def run(self, start_new_process=False): # FIXME: temporary workaround to copy from output the additional # information of the source (such as fTotalSkippedEvents) - s = {} for source in self.source_manager.sources.values(): # WARNING: when multithread, the sources are stored in # simulation_output.sources_by_thread @@ -2219,12 +2234,7 @@ def run(self, start_new_process=False): s = output.get_source(source.name) except: continue - if "total_zero_events" in s.__dict__: - source.total_zero_events = s.__dict__["total_zero_events"] - source.total_skipped_events = s.__dict__["total_skipped_events"] - if "particle_generators" in s.__dict__: - source.particle_generators = s.__dict__["particle_generators"] - source.num_entries = s.__dict__["num_entries"] + source.recover_user_output(s) else: # Nothing special to do if the simulation engine ran in the native python process diff --git a/opengate/sources/base.py b/opengate/sources/base.py index 0119d5f6c..f279e7438 100644 --- a/opengate/sources/base.py +++ b/opengate/sources/base.py @@ -1,4 +1,6 @@ +import os import numpy as np +import opengate_core as g4 from ..actors.base import _setter_hook_attached_to from ..base import GateObject, DynamicGateObject, process_cls @@ -71,6 +73,8 @@ def __init__(self, *args, **kwargs): GateObject.__init__(self, *args, **kwargs) # all times intervals self.run_timing_intervals = None + self.g4_thread_sources = [] + self.g4_thread_sources_index = 0 def __initcpp__(self): """Nothing to do in the base class.""" @@ -115,11 +119,53 @@ def initialize(self, run_timing_intervals): self.InitializeUserInfo(self.user_info) def add_to_source_manager(self, source_manager): - source_manager.AddSource(self) + if hasattr(self, "g4_source") and self.g4_source is not None: + source_manager.AddSource(self.g4_source) + else: + source_manager.AddSource(self) + + def close(self): + # remove the g4 objects + for v in list(self.__dict__.keys()): + if "g4_" in v: + self.__dict__[v] = None + # close the base GateObject + GateObject.close(self) def prepare_output(self): pass + def pre_create_g4_sources(self, num_instances): + self.g4_thread_sources = [] + self.g4_thread_sources_index = 0 + for _ in range(num_instances): + g4_src = self.create_g4_source() + if g4_src is not None: + self.g4_thread_sources.append(g4_src) + + def get_next_g4_source(self): + if self.g4_thread_sources: + tid = g4.G4GetThreadId() + idx = tid + 1 if tid >= 0 else 0 + if idx < len(self.g4_thread_sources): + return self.g4_thread_sources[idx] + return None + + def create_g4_source(self): + return None + + def initialize_g4_source(self, g4_source, run_timing_intervals): + pass + + def gather_outputs(self, thread_sources): + pass + + def recover_user_output(self, s): + pid = os.getpid() + print(f"(python) recover_user_output {self.name} pid={pid}") + for k, v in s.user_info.items(): + self.user_info[k] = v + def can_predict_number_of_events(self): return True @@ -145,4 +191,46 @@ def check_ui_activity(self, ui): ui.activity = 0 +class DebugSource(SourceBase): + + user_info_defaults = { + "debug_flag": (False, {"doc": "Fake parameter."}), + "debug_value": (0.0, {"doc": "Fake parameter."}), + } + + def __init__(self, *args, **kwargs): + pid = os.getpid() + print(f"(python) DebugSource::__init__ pid={pid}") + SourceBase.__init__(self, *args, **kwargs) + + def create_g4_source(self): + pid = os.getpid() + print(f"(python) DebugSource::create_g4_source pid={pid}") + return g4.GateDebugSource() + + def initialize_g4_source(self, g4_source, run_timing_intervals): + pid = os.getpid() + print(f"(python) DebugSource::initialize_g4_source pid={pid}") + self.initialize_start_end_time(run_timing_intervals) + self.check_ui_activity(self.user_info) + g4_source.InitializeUserInfo(self.user_info) + + def initialize_start_end_time(self, run_timing_intervals): + pid = os.getpid() + print(f"(python) DebugSource::initialize_start_end_time {self.name} pid={pid}") + SourceBase.initialize_start_end_time(self, run_timing_intervals) + + def gather_outputs(self, thread_sources): + values = [ + g4_src.GetDebugValue() for g4_src in thread_sources if g4_src is not None + ] + print(f"(python) DebugSource::gather_outputs values = {values}") + if values: + self.debug_value = np.sum(np.array(values)) + print( + f"(python) DebugSource::gather_outputs selected max value = {self.debug_value}" + ) + + process_cls(SourceBase) +process_cls(DebugSource) diff --git a/opengate/sources/beamsources.py b/opengate/sources/beamsources.py index 13ea054e4..25cdef89b 100644 --- a/opengate/sources/beamsources.py +++ b/opengate/sources/beamsources.py @@ -97,7 +97,7 @@ def validate(self, parent_obj, attr_name: str, parent_context: str = None): _pbs_dir_validator = PBSDirectionValidator() -class IonPencilBeamSource(GenericSource, g4.GatePencilBeamSource): +class IonPencilBeamSource(GenericSource): """ Pencil Beam source """ @@ -110,19 +110,15 @@ class IonPencilBeamSource(GenericSource, g4.GatePencilBeamSource): } def __init__(self, *args, **kwargs): - super().__init__(self, *args, **kwargs) - self.__initcpp__() + GenericSource.__init__(self, *args, **kwargs) self.position.type = "disc" self._dir_validator = PBSDirectionValidator() - def __initcpp__(self): - g4.GatePencilBeamSource.__init__(self) + def create_g4_source(self): + return g4.GatePencilBeamSource() - def initialize(self, run_timing_intervals): - GenericSource.initialize(self, run_timing_intervals) - -class TreatmentPlanPBSource(SourceBase, g4.GateTreatmentPlanPBSource): +class TreatmentPlanPBSource(SourceBase): """ Treatment Plan source Pencil Beam """ @@ -234,8 +230,7 @@ class TreatmentPlanPBSource(SourceBase, g4.GateTreatmentPlanPBSource): } def __init__(self, *args, **kwargs): - super().__init__(self, *args, **kwargs) - self.__initcpp__() + SourceBase.__init__(self, *args, **kwargs) # initialize internal members self.spots = None self.rotation = None @@ -246,19 +241,16 @@ def __init__(self, *args, **kwargs): self.proportion_factor_x = None self.proportion_factor_y = None - def __initcpp__(self): - g4.GateTreatmentPlanPBSource.__init__(self) + def create_g4_source(self): + return g4.GateTreatmentPlanPBSource() - def initialize(self, run_timing_intervals): + def initialize_g4_source(self, g4_source, run_timing_intervals): if not self.beam_data_dict and not self.plan_path: raise ValueError( "User must provide either a treatment plan file path or a beam data dictionary " "with spots and gantry angle." ) - # if len(self.user_info.n_primaries_vector) != len(self.user_info.run_timing_intervals): - # raise ValueError("Particles per run must have the same length of the number of runs") - # set pbs param self._set_pbs_param_all_spots() @@ -272,11 +264,14 @@ def initialize(self, run_timing_intervals): if len(words) > 3: self.ion.E = words[3] - # initialize - SourceBase.initialize(self, run_timing_intervals) + self.initialize_start_end_time(run_timing_intervals) + self.check_ui_activity(self.user_info) + g4_source.InitializeUserInfo(self.user_info) def get_generated_primaries(self): - return self.g4_source.GetGeneratedPrimaries() + if self.g4_thread_sources: + return self.g4_thread_sources[0].GetGeneratedPrimaries() + return [] def _set_pbs_param_all_spots(self): # initialize vectors every time, otherwise issues with MT diff --git a/opengate/sources/gansources.py b/opengate/sources/gansources.py index 560f1e662..74145865b 100644 --- a/opengate/sources/gansources.py +++ b/opengate/sources/gansources.py @@ -267,7 +267,7 @@ def generate_condition(self, n): return np.column_stack((p, v)) -class GANSource(GenericSource, g4.GateGANSource): +class GANSource(GenericSource): """ GAN source: the Generator produces particles Input is a neural network Generator trained with a GAN @@ -385,31 +385,27 @@ class GANSource(GenericSource, g4.GateGANSource): } def __init__(self, *args, **kwargs): - super().__init__(self, *args, **kwargs) - self.__initcpp__() + GenericSource.__init__(self, *args, **kwargs) - def __initcpp__(self): - g4.GateGANSource.__init__(self) - - def initialize(self, run_timing_intervals): - # FIXME -> check input user_info - # initialize the mother class generic source - GenericSource.initialize(self, run_timing_intervals) + def create_g4_source(self): + return g4.GateGANSource() + def initialize_g4_source(self, g4_source, run_timing_intervals): # default generator or set by the user if self.user_info.generator is None: self.set_default_generator() gen = self.user_info.generator # initialize the generator (read the GAN) - # this function must have 1) the generator function 2) the associated info gen.initialize() # set the function pointer to the cpp side - self.SetGeneratorFunction(gen.generator) + g4_source.SetGeneratorFunction(gen.generator) # set the parameters to the cpp side - self.SetGeneratorInfo(gen.gan_info) + g4_source.SetGeneratorInfo(gen.gan_info) + + GenericSource.initialize_g4_source(self, g4_source, run_timing_intervals) def set_default_generator(self): # non-conditional generator @@ -426,18 +422,17 @@ def set_default_generator(self): ) -class GANPairsSource(GANSource, g4.GateGANPairSource): +class GANPairsSource(GANSource): """ GAN source: the Generator produces pairs of particles (for PET) Input is a neural network Generator trained with a GAN """ def __init__(self, *args, **kwargs): - super().__init__(self, *args, **kwargs) - self.__initcpp__() + GANSource.__init__(self, *args, **kwargs) - def __initcpp__(self): - g4.GateGANPairSource.__init__(self) + def create_g4_source(self): + return g4.GateGANPairSource() def set_default_generator(self): # non-conditional generator diff --git a/opengate/sources/generic.py b/opengate/sources/generic.py index d768df4f4..818c6874d 100644 --- a/opengate/sources/generic.py +++ b/opengate/sources/generic.py @@ -17,6 +17,7 @@ all_beta_plus_radionuclides, compute_cdf_and_total_yield, get_spectrum, + _setter_hook_generic_source_particle, ) @@ -86,24 +87,6 @@ def visualization_parameters(): ) -def _setter_hook_generic_source_particle(self, particle): - # The particle parameter must be a str - if not isinstance(particle, str): - fatal(f"the .particle user info must be a str, while it is {type(str)}") - # if it does not start with ion, we consider this is a simple particle (gamma, e+, etc.) - if not particle.startswith("ion"): - return particle - # if start with ion, it is like 'ion 9 18' with Z A E - words = particle.split(" ") - if len(words) > 1: - self.ion.Z = int(words[1]) - if len(words) > 2: - self.ion.A = int(words[2]) - if len(words) > 3: - self.ion.E = int(words[3]) - return particle - - class PositionValidator(UserInfoValidatorBase): """Validates the 'position' Box.""" @@ -347,7 +330,7 @@ def validate(self, parent_obj, attr_name: str, parent_context: str = None): return context_name -class GenericSource(SourceBase, g4.GateGenericSource): +class GenericSource(SourceBase): """ GenericSource close to the G4 SPS, but a bit simpler. The G4 source created by this class is GateGenericSource. @@ -439,8 +422,7 @@ class GenericSource(SourceBase, g4.GateGenericSource): } def __init__(self, *args, **kwargs): - super().__init__(self, *args, **kwargs) - self.__initcpp__() + SourceBase.__init__(self, *args, **kwargs) # to validate the parameters self._pos_validator = PositionValidator() self._dir_validator = DirectionValidator() @@ -458,10 +440,10 @@ def __init__(self, *args, **kwargs): if len(words) > 3: self.user_info.ion.E = words[3] - def __initcpp__(self): - g4.GateGenericSource.__init__(self) + def create_g4_source(self): + return g4.GateGenericSource() - def initialize(self, run_timing_intervals): + def initialize_g4_source(self, g4_source, run_timing_intervals): # Check the sub-parameters self._pos_validator.validate(self, "position") self._ene_validator.validate(self, "energy") @@ -484,10 +466,10 @@ def initialize(self, run_timing_intervals): # total = total * 1000 # (because was in MeV) # self.user_info.activity *= total self.energy.is_cdf = True - self.SetEnergyCDF(ene) - self.SetProbabilityCDF(cdf) + g4_source.SetEnergyCDF(ene) + g4_source.SetProbabilityCDF(cdf) - self.update_tac_activity() + self.update_tac_activity(g4_source) # histogram parameters: histogram_weight, histogram_energy" ene = self.energy @@ -505,8 +487,9 @@ def initialize(self, run_timing_intervals): if self.user_particle_life_time < 0: self.user_particle_life_time = 0 - # initialize - SourceBase.initialize(self, run_timing_intervals) + self.initialize_start_end_time(run_timing_intervals) + self.check_ui_activity(self.user_info) + g4_source.InitializeUserInfo(self.user_info) # warning for non-used ? def check_confine(self, ui): @@ -519,13 +502,19 @@ def check_confine(self, ui): f"confine is used, while position.type is point ... really ?" ) - def prepare_output(self): - SourceBase.prepare_output(self) - # store the output from G4 objects - self.total_zero_events = self.GetTotalZeroEvents() - self.total_skipped_events = self.GetTotalSkippedEvents() - - def update_tac_activity(self): + def gather_outputs(self, thread_sources): + self.total_zero_events = sum( + g4_src.GetTotalZeroEvents() + for g4_src in thread_sources + if g4_src is not None + ) + self.total_skipped_events = sum( + g4_src.GetTotalSkippedEvents() + for g4_src in thread_sources + if g4_src is not None + ) + + def update_tac_activity(self, g4_source): if self.tac_times is None and self.tac_activities is None: return if len(self.tac_times) != len(self.tac_activities): @@ -536,7 +525,7 @@ def update_tac_activity(self): # may start later than the simulation timing self.start_time = self.tac_times[0] self.activity = self.tac_activities[0] - self.SetTAC(self.tac_times, self.tac_activities) + g4_source.SetTAC(self.tac_times, self.tac_activities) def can_predict_number_of_events(self): aa = self.direction.angular_acceptance diff --git a/opengate/sources/lastvertexsources.py b/opengate/sources/lastvertexsources.py index b90db6852..6fa120ed4 100644 --- a/opengate/sources/lastvertexsources.py +++ b/opengate/sources/lastvertexsources.py @@ -11,21 +11,22 @@ from ..exception import fatal, warning -class LastVertexSource(SourceBase, g4.GateLastVertexSource): +class LastVertexSource(SourceBase): """ The source used to replay position, energy, direction and weight of last vertex particles actor """ def __init__(self, *args, **kwargs): - super().__init__(self, *args, **kwargs) - self.__initcpp__() + SourceBase.__init__(self, *args, **kwargs) - def __initcpp__(self): - g4.GateLastVertexSource.__init__(self) + def create_g4_source(self): + return g4.GateLastVertexSource() - def initialize(self, run_timing_intervals): + def initialize_g4_source(self, g4_source, run_timing_intervals): self.user_info.n = np.zeros(len(run_timing_intervals)) + 1 - SourceBase.initialize(self, run_timing_intervals) + self.initialize_start_end_time(run_timing_intervals) + self.check_ui_activity(self.user_info) + g4_source.InitializeUserInfo(self.user_info) process_cls(LastVertexSource) diff --git a/opengate/sources/phidsources.py b/opengate/sources/phidsources.py index 7ea5f00fa..1a2ea4945 100644 --- a/opengate/sources/phidsources.py +++ b/opengate/sources/phidsources.py @@ -50,7 +50,7 @@ class PhotonFromIonDecaySource(GenericSource): } def __init__(self, *args, **kwargs): - super().__init__(self, *args, **kwargs) + GenericSource.__init__(self, *args, **kwargs) # list of all sub_sources self.sub_sources = [] # False if this is the main first source @@ -61,36 +61,36 @@ def __init__(self, *args, **kwargs): # used during sub_sources creation self.daughters = None + # metadata for daughter/mother ions + self.ion_gamma_mother = None + self.ion_gamma_daughter = None + # debug self.log = "" self.debug_first_daughter_only = False - def __initcpp__(self): - g4.GateGenericSource.__init__(self) - - def initialize(self, run_timing_intervals): + def pre_create_g4_sources(self, num_instances): if not self.is_a_sub_source: - # this is called only one single time to build the list of sub_sources - if (g4.IsMultithreadedApplication() and g4.G4GetThreadId() == -1) or ( - not g4.IsMultithreadedApplication() - ): - phid_build_all_sub_sources(self) + phid_build_all_sub_sources(self) + for sub_source in self.sub_sources: + sub_source.pre_create_g4_sources(num_instances) + else: + GenericSource.pre_create_g4_sources(self, num_instances) - # super class conventional initialization + def initialize(self, run_timing_intervals): + if not self.is_a_sub_source and not self.sub_sources: + phid_build_all_sub_sources(self) - GenericSource.initialize(self, run_timing_intervals) self.initialize_start_end_time(run_timing_intervals) - for sub_source in self.sub_sources: + def add_to_source_manager(self, source_manager): + for sub_source in self.sub_sources: sub_source.start_time = self.start_time sub_source.end_time = self.end_time - # this will initialize and set user_info to the cpp side - # for all sub_sources - for sub_source in self.sub_sources: # get tac from decay p = Box(sub_source.tac_from_decay_parameters) - sub_source.tac_times, sub_source.tac_activities = get_tac_from_decay( + tac_times, tac_activities = get_tac_from_decay( p.ion_name, p.daughter, sub_source.activity, @@ -98,32 +98,52 @@ def initialize(self, run_timing_intervals): sub_source.end_time, p.bins, ) - # update the tac - update_sub_source_tac_activity(sub_source) + # Get the thread-local C++ source instance for the sub_source + g4_src = sub_source.get_next_g4_source() + + # update the tac on the C++ source + update_sub_source_tac_activity( + sub_source, g4_src, tac_times, tac_activities + ) # check self.check_ui_activity(sub_source) self.check_confine(sub_source) - # final initialize to cpp side - # sub_source.n = np.array([sub_source.n],dtype = int) - sub_source.InitializeUserInfo(sub_source.user_info) + # Initialize the thread-local C++ source instance + if g4_src is not None: + sub_source.initialize_g4_source(g4_src, self.run_timing_intervals) + source_manager.AddSource(g4_src) # dump log if self.user_info.dump_log is not None: with open(self.user_info.dump_log, "w") as outfile: outfile.write(self.log) - def add_to_source_manager(self, source_manager): - for g4_source in self.sub_sources: - source_manager.AddSource(g4_source) + def prepare_output(self): + for sub_source in self.sub_sources: + sub_source.prepare_output() + if sub_source.g4_thread_sources: + sub_source.gather_outputs(sub_source.g4_thread_sources) + sub_source.g4_thread_sources = [] + sub_source.g4_thread_sources_index = 0 + self.total_zero_events = sum( + s.total_zero_events + for s in self.sub_sources + if hasattr(s, "total_zero_events") + ) + self.total_skipped_events = sum( + s.total_skipped_events + for s in self.sub_sources + if hasattr(s, "total_skipped_events") + ) -def update_sub_source_tac_activity(sub_source): - if sub_source.tac_times is None and sub_source.tac_activities is None: +def update_sub_source_tac_activity(sub_source, g4_source, tac_times, tac_activities): + if tac_times is None and tac_activities is None: return - n = len(sub_source.tac_times) - if n != len(sub_source.tac_activities): + n = len(tac_times) + if n != len(tac_activities): fatal( f"option tac_activities must have the same size than tac_times in source '{sub_source.name}'" ) @@ -136,26 +156,25 @@ def update_sub_source_tac_activity(sub_source): # it is important to set the starting time for this source as the tac # may start later than the simulation timing i = 0 - while i < len(sub_source.tac_activities) and sub_source.tac_activities[i] <= 0: + while i < len(tac_activities) and tac_activities[i] <= 0: i += 1 - if i >= len(sub_source.tac_activities): + if i >= len(tac_activities): # gate.warning(f"Source '{sub_source.name}' TAC with zero activity.") sub_source.start_time = sub_source.end_time + 1 * sec else: - sub_source.start_time = sub_source.tac_times[i] + sub_source.start_time = tac_times[i] # IMPORTANT : activities must be x by total here # (not before, because it can be called several times in MT mode) - sub_source.SetTAC( - sub_source.tac_times, np.array(sub_source.tac_activities) * total - ) + if g4_source is not None: + g4_source.SetTAC(tac_times, np.array(tac_activities) * total) if sub_source.verbose: print( f"GammaFromIon source {sub_source.name} total = {total * 100:8.2f}% " f" gammas lines = {len(sub_source.energy.spectrum_weights)} " - f" total activity = {sum(sub_source.tac_activities) / Bq:10.3f}" - f" first activity = {sub_source.tac_activities[0] / Bq:5.2f}" - f" last activity = {sub_source.tac_activities[-1] / Bq:5.2f}" + f" total activity = {sum(tac_activities) / Bq:10.3f}" + f" first activity = {tac_activities[0] / Bq:5.2f}" + f" last activity = {tac_activities[-1] / Bq:5.2f}" ) @@ -503,7 +522,7 @@ def isomeric_transition_load_from_iaea_website(a, rad_name): df2 = lc_read_csv(url) if not df2.empty: # Identify overlapping columns - overlapping_columns = df.columns.plane_intersection_torch(df2.columns) + overlapping_columns = df.columns.intersection(df2.columns) # Convert columns in df2 to the same type as df1 for col in overlapping_columns: @@ -747,8 +766,8 @@ def phid_build_one_sub_source(stype, source, daughter, ene, w, first_nuclide): s.verbose = source.verbose s.particle = "gamma" s.energy.type = "spectrum_discrete" - s.energy.ion_gamma_mother = Box({"z": first_nuclide.Z, "a": first_nuclide.A}) - s.energy.ion_gamma_daughter = ion_gamma_daughter + s.ion_gamma_mother = Box({"z": first_nuclide.Z, "a": first_nuclide.A}) + s.ion_gamma_daughter = ion_gamma_daughter s.energy.spectrum_weights = w s.energy.spectrum_energies = ene s.activity = source.activity diff --git a/opengate/sources/phspsources.py b/opengate/sources/phspsources.py index 7bf225dd3..cecaa0e0b 100644 --- a/opengate/sources/phspsources.py +++ b/opengate/sources/phspsources.py @@ -114,7 +114,7 @@ def get_entry_start(self, entry_start): ) return n - def generate(self, source, pid): + def generate(self, g4_source, pid): """ Main function called from C++ to generate a batch of particles. """ @@ -130,7 +130,7 @@ def generate(self, source, pid): self.cycle_changed_flag = False # --- 2. Calculate Batch Size --- - requested_batch_size = source.batch_size + requested_batch_size = self.phsp_source.batch_size if self.current_index + requested_batch_size >= self.num_entries: requested_batch_size = self.num_entries - self.current_index @@ -139,11 +139,11 @@ def generate(self, source, pid): if requested_batch_size == 0: self.current_index = 0 - requested_batch_size = min(source.batch_size, self.num_entries) + requested_batch_size = min(self.phsp_source.batch_size, self.num_entries) self.cycle_count += 1 self.cycle_changed_flag = True - if source.verbose_batch: + if self.phsp_source.verbose_batch: print( f"Thread {g4.G4GetThreadId()} reading {requested_batch_size} events from index {self.current_index}" ) @@ -190,39 +190,41 @@ def get_data(key, dtype, must_exist=True): # We assign to 'self.X' to ensure the Python object survives # as long as the C++ side needs it (until the next batch overwrites it). - self.pos_x = get_data(source.position_key_x, np.float32) + self.pos_x = get_data(self.phsp_source.position_key_x, np.float32) actual_size = len(self.pos_x) - self.pos_y = get_data(source.position_key_y, np.float32) - self.pos_z = get_data(source.position_key_z, np.float32) + self.pos_y = get_data(self.phsp_source.position_key_y, np.float32) + self.pos_z = get_data(self.phsp_source.position_key_z, np.float32) - self.dir_x = get_data(source.direction_key_x, np.float32) - self.dir_y = get_data(source.direction_key_y, np.float32) - self.dir_z = get_data(source.direction_key_z, np.float32) + self.dir_x = get_data(self.phsp_source.direction_key_x, np.float32) + self.dir_y = get_data(self.phsp_source.direction_key_y, np.float32) + self.dir_z = get_data(self.phsp_source.direction_key_z, np.float32) - self.energy = get_data(source.energy_key, np.float32) + self.energy = get_data(self.phsp_source.energy_key, np.float32) # Weights self.weight = None - if source.weight_key: - self.weight = get_data(source.weight_key, np.float32, must_exist=False) + if self.phsp_source.weight_key: + self.weight = get_data( + self.phsp_source.weight_key, np.float32, must_exist=False + ) if self.weight is None: self.weight = np.ones(actual_size, dtype=np.float32) # PDG Code self.pdg = None - if not source.particle: - self.pdg = get_data(source.PDGCode_key, np.int32) + if not self.phsp_source.particle: + self.pdg = get_data(self.phsp_source.PDGCode_key, np.int32) # --- 6. Transforms (Modify SELF) --- - if source.translate_position: - self.pos_x += float(source.position.translation[0]) - self.pos_y += float(source.position.translation[1]) - self.pos_z += float(source.position.translation[2]) + if self.phsp_source.translate_position: + self.pos_x += float(self.phsp_source.position.translation[0]) + self.pos_y += float(self.phsp_source.position.translation[1]) + self.pos_z += float(self.phsp_source.position.translation[2]) - if source.rotate_direction: + if self.phsp_source.rotate_direction: points = np.column_stack((self.dir_x, self.dir_y, self.dir_z)) - r = Rotation.from_matrix(source.position.rotation) + r = Rotation.from_matrix(self.phsp_source.position.rotation) rotated = r.apply(points) self.dir_x = np.ascontiguousarray(rotated[:, 0], dtype=np.float32) self.dir_y = np.ascontiguousarray(rotated[:, 1], dtype=np.float32) @@ -233,22 +235,22 @@ def get_data(key, dtype, must_exist=True): fatal(f"Size mismatch: Pos {actual_size} vs Energy {len(self.energy)}") # Pass the SELF variables to C++ - source.SetPositionXBatch(self.pos_x) - source.SetPositionYBatch(self.pos_y) - source.SetPositionZBatch(self.pos_z) - source.SetDirectionXBatch(self.dir_x) - source.SetDirectionYBatch(self.dir_y) - source.SetDirectionZBatch(self.dir_z) - source.SetEnergyBatch(self.energy) - source.SetWeightBatch(self.weight) + g4_source.SetPositionXBatch(self.pos_x) + g4_source.SetPositionYBatch(self.pos_y) + g4_source.SetPositionZBatch(self.pos_z) + g4_source.SetDirectionXBatch(self.dir_x) + g4_source.SetDirectionYBatch(self.dir_y) + g4_source.SetDirectionZBatch(self.dir_z) + g4_source.SetEnergyBatch(self.energy) + g4_source.SetWeightBatch(self.weight) if self.pdg is not None: - source.SetPDGCodeBatch(self.pdg) + g4_source.SetPDGCodeBatch(self.pdg) return actual_size -class PhaseSpaceSource(SourceBase, g4.GatePhaseSpaceSource): +class PhaseSpaceSource(SourceBase): """ Source of particles from a (root) phase space. Read position + direction + energy + weight from the root and use them as event. @@ -441,16 +443,12 @@ class PhaseSpaceSource(SourceBase, g4.GatePhaseSpaceSource): } def __init__(self, *args, **kwargs): - super().__init__(self, *args, **kwargs) - self.__initcpp__() + SourceBase.__init__(self, *args, **kwargs) # there will be one particle generator per thread self.particle_generators = {} # number of entries in the phsp root file self.num_entries = None - def __initcpp__(self): - g4.GatePhaseSpaceSource.__init__(self) - def __getstate__(self): # the particle generator cannot (?) being pickled # we convert in a dict with the cycle count values @@ -461,28 +459,27 @@ def __getstate__(self): state_dict = super().__getstate__() return state_dict - def initialize(self, run_timing_intervals): - # create a generator for each thread - tid = g4.G4GetThreadId() - self.particle_generators[tid] = PhaseSpaceSourceGenerator(tid) + def create_g4_source(self): + return g4.GatePhaseSpaceSource() - # initialize source - SourceBase.initialize(self, run_timing_intervals) + def initialize_g4_source(self, g4_source, run_timing_intervals): + # Calculate the target thread ID for this g4_source + if g4.IsMultithreadedApplication(): + try: + thread_idx = self.g4_thread_sources.index(g4_source) + except ValueError: + thread_idx = 0 + tid = thread_idx - 1 + else: + tid = 0 - # check user info - if self.position_key_x is None: - self.position_key_x = f"{self.position_key}_X" - if self.position_key_y is None: - self.position_key_y = f"{self.position_key}_Y" - if self.position_key_z is None: - self.position_key_z = f"{self.position_key}_Z" + # Master source does not need generator + if tid < 0: + return - if self.direction_key_x is None: - self.direction_key_x = f"{self.direction_key}_X" - if self.direction_key_y is None: - self.direction_key_y = f"{self.direction_key}_Y" - if self.direction_key_z is None: - self.direction_key_z = f"{self.direction_key}_Z" + # Create/initialize the generator for this thread + if tid not in self.particle_generators: + self.particle_generators[tid] = PhaseSpaceSourceGenerator(tid) # check if the source should generate particles until the second one # which is identified as primary by name, PDGCode and above a threshold @@ -502,14 +499,30 @@ def initialize(self, run_timing_intervals): if not g4.IsMultithreadedApplication(): self.entry_start = 0 else: - # create an entry_start array with the correct number of start entries - # all entries are spaced by the number of particles/thread - # FIXME: check this line. I corrected it because it seemed like a typo (NK) n_threads = self.simulation.number_of_threads - # n_threads = self.simulation.phsp_source.number_of_threads - step = np.ceil(self.n / n_threads) + 1 # Specify the increment value + step = np.ceil(self.n / n_threads) + 1 self.entry_start = [i * step for i in range(n_threads)] + # check user info + if self.position_key_x is None: + self.position_key_x = f"{self.position_key}_X" + if self.position_key_y is None: + self.position_key_y = f"{self.position_key}_Y" + if self.position_key_z is None: + self.position_key_z = f"{self.position_key}_Z" + + if self.direction_key_x is None: + self.direction_key_x = f"{self.direction_key}_X" + if self.direction_key_y is None: + self.direction_key_y = f"{self.direction_key}_Y" + if self.direction_key_z is None: + self.direction_key_z = f"{self.direction_key}_Z" + + # Initialize base start/end times and check activity on the python side + self.initialize_start_end_time(run_timing_intervals) + self.check_ui_activity(self.user_info) + g4_source.InitializeUserInfo(self.user_info) + # initialize the generator (read the phsp file) self.particle_generators[tid].initialize(self) @@ -517,17 +530,19 @@ def initialize(self, run_timing_intervals): self.num_entries = self.particle_generators[tid].num_entries # set the function pointer to the cpp side - self.SetGeneratorFunction(self.particle_generators[tid].generate) + g4_source.SetGeneratorFunction(self.particle_generators[tid].generate) @property def cycle_count(self): if not g4.IsMultithreadedApplication(): - tid = g4.G4GetThreadId() - return self.particle_generators[tid].cycle_count + if not self.particle_generators: + return 0 + key = list(self.particle_generators.keys())[0] + return self.particle_generators[key].cycle_count else: s = " ".join( str(self.particle_generators[tid].cycle_count) - for tid in self.particle_generators.keys() + for tid in sorted(self.particle_generators.keys()) ) return s diff --git a/opengate/sources/utility.py b/opengate/sources/utility.py index 4c503c844..de54ae01b 100644 --- a/opengate/sources/utility.py +++ b/opengate/sources/utility.py @@ -320,3 +320,21 @@ def __get_rad_beta_spectrum(rad: str): data.energies = np.array(energies) return data + + +def _setter_hook_generic_source_particle(self, particle): + # The particle parameter must be a str + if not isinstance(particle, str): + fatal(f"the .particle user info must be a str, while it is {type(str)}") + # if it does not start with ion, we consider this is a simple particle (gamma, e+, etc.) + if not particle.startswith("ion"): + return particle + # if start with ion, it is like 'ion 9 18' with Z A E + words = particle.split(" ") + if len(words) > 1: + self.ion.Z = int(words[1]) + if len(words) > 2: + self.ion.A = int(words[2]) + if len(words) > 3: + self.ion.E = int(words[3]) + return particle diff --git a/opengate/sources/voxelsources.py b/opengate/sources/voxelsources.py index 0b2e8132c..571312c74 100644 --- a/opengate/sources/voxelsources.py +++ b/opengate/sources/voxelsources.py @@ -12,7 +12,7 @@ from ..actors.dynamicactors import SourceActivityImageChanger -class VoxelSource(GenericSource, g4.GateVoxelSource): +class VoxelSource(GenericSource): """ VoxelSource = 3D activity distribution. Sampled with cumulative distribution functions. @@ -34,13 +34,13 @@ class VoxelSource(GenericSource, g4.GateVoxelSource): } def __init__(self, *args, **kwargs): - self.__initcpp__() - super().__init__(self, *args, **kwargs) + GenericSource.__init__(self, *args, **kwargs) # the loaded image self._current_itk_image = None - - def __initcpp__(self): - g4.GateVoxelSource.__init__(self) + # cached CDFs + self._cdf_x = None + self._cdf_y = None + self._cdf_z = None def create_changers(self): changers = super().create_changers() @@ -61,10 +61,10 @@ def create_changers(self): ) return changers - def set_transform_from_user_info(self): + def set_transform_from_user_info(self, g4_source): src_info = get_info_from_image(self._current_itk_image) # get the pointer to SPSVoxelPosDistribution - pg = self.GetSPSVoxelPosDistribution() + pg = g4_source.GetSPSVoxelPosDistribution() # update cpp image info (no need to allocate) update_image_py_to_cpp(self._current_itk_image, pg.cpp_edep_image, False) # set spacing @@ -77,44 +77,57 @@ def set_transform_from_user_info(self): ) pg.cpp_edep_image.set_origin(c) - def cumulative_distribution_functions(self): + def cumulative_distribution_functions(self, g4_source): """ Compute the Cumulative Distribution Function of the image Composed of: CDF_Z = 1D, CDF_Y = 2D, CDF_X = 3D """ - cdf_x, cdf_y, cdf_z = compute_image_3D_CDF(self._current_itk_image) + if self._cdf_x is None or self._cdf_y is None or self._cdf_z is None: + self._cdf_x, self._cdf_y, self._cdf_z = compute_image_3D_CDF( + self._current_itk_image + ) # set CDF to the position generator - pg = self.GetSPSVoxelPosDistribution() - pg.SetCumulativeDistributionFunction(cdf_z, cdf_y, cdf_x) + pg = g4_source.GetSPSVoxelPosDistribution() + pg.SetCumulativeDistributionFunction(self._cdf_z, self._cdf_y, self._cdf_x) def update_activity_image(self, filename): # read source image self._current_itk_image = itk.imread(ensure_filename_is_str(filename)) - # compute position - self.set_transform_from_user_info() - - # create Cumulative Distribution Function - self.cumulative_distribution_functions() - - def initialize(self, run_timing_intervals): - self.update_activity_image(self.image) - - # FIXME -> check other option in position not used here - + # Reset CDF cache + self._cdf_x = None + self._cdf_y = None + self._cdf_z = None + + # update all thread-local sources + for g4_source in self.g4_thread_sources: + # compute position + self.set_transform_from_user_info(g4_source) + # create Cumulative Distribution Function + self.cumulative_distribution_functions(g4_source) + + def create_g4_source(self): + return g4.GateVoxelSource() + + def initialize_g4_source(self, g4_source, run_timing_intervals): + if self._current_itk_image is None: + self._current_itk_image = itk.imread(ensure_filename_is_str(self.image)) + self.set_transform_from_user_info(g4_source) + self.cumulative_distribution_functions(g4_source) # initialise standard options (particle energy, etc.) - # we temporarily set the position attribute to reuse - # the GenericSource verification - GenericSource.initialize(self, run_timing_intervals) + GenericSource.initialize_g4_source(self, g4_source, run_timing_intervals) -class VoxelizedPromptGammaTLESource(GenericSource, g4.GateVoxelSource): +class VoxelizedPromptGammaTLESource(VoxelSource): """ VoxelizedPromptGammaTLESource = 3D PG distribution. Sampled with cumulative distribution functions. """ + def create_g4_source(self): + return g4.GateVoxelizedPromptGammaTLESource() + process_cls(VoxelSource) process_cls(VoxelizedPromptGammaTLESource) diff --git a/opengate/tests/src/actors/test023_filters_auxiliary_attribute.py b/opengate/tests/src/actors/test023_filters_auxiliary_attribute.py old mode 100644 new mode 100755 diff --git a/opengate/tests/src/actors/test023_interaction_counter_auxiliary_attribute.py b/opengate/tests/src/actors/test023_interaction_counter_auxiliary_attribute.py old mode 100644 new mode 100755 diff --git a/opengate/tests/src/actors/test023_last_interaction_position_auxiliary_attribute.py b/opengate/tests/src/actors/test023_last_interaction_position_auxiliary_attribute.py old mode 100644 new mode 100755 diff --git a/opengate/tests/src/actors/test023_last_interaction_position_auxiliary_attribute_propagation.py b/opengate/tests/src/actors/test023_last_interaction_position_auxiliary_attribute_propagation.py old mode 100644 new mode 100755 diff --git a/opengate/tests/src/actors/test102_debug_actor.py b/opengate/tests/src/actors/test102_debug_actor.py new file mode 100755 index 000000000..6c1a5df94 --- /dev/null +++ b/opengate/tests/src/actors/test102_debug_actor.py @@ -0,0 +1,73 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import os +from scipy.spatial.transform import Rotation +import opengate as gate +from opengate.tests import utility + +if __name__ == "__main__": + paths = utility.get_default_test_paths(__file__, "", "test102_debug_actor") + + # create the simulation + sim = gate.Simulation() + sim.visu = False + sim.random_seed = "auto" + sim.output_dir = paths.output + sim.number_of_threads = 2 + if os.name == "nt": + sim.number_of_threads = 1 + sim.store_json_archive = True + print(paths) + + # units + m = gate.g4_units.m + cm = gate.g4_units.cm + mm = gate.g4_units.mm + MeV = gate.g4_units.MeV + + # world size + world = sim.world + world.size = [1 * m, 1 * m, 1 * m] + + # fake volume + fake = sim.add_volume("Box", "fake") + fake.size = [40 * cm, 40 * cm, 40 * cm] + fake.translation = [1 * cm, 2 * cm, 3 * cm] + fake.rotation = Rotation.from_euler("x", 10, degrees=True).as_matrix() + fake.material = "G4_AIR" + fake.color = [1, 0, 1, 1] + + # source + source = sim.add_source("GenericSource", "mysource") + source.energy.mono = 1 * MeV + source.particle = "gamma" + source.position.type = "sphere" + source.position.radius = 10 * mm + source.direction.type = "iso" + source.n = 2 + + # add debug actor + debug = sim.add_actor("DebugActor", "debug") + debug.attached_to = "fake" + debug.debug_flag = True + + # add stat actor + stat = sim.add_actor("SimulationStatisticsActor", "Stats") + stat.track_types_flag = True + + # start simulation in another process + sim.run(start_new_process=True) + + # print results at the end + print(stat) + print(debug) + + print(f"Simulation json saved in {sim.output_dir / sim.json_archive_filename}") + + # assertions to verify MT execution and output recovery + assert stat.counts.events == 2 * sim.number_of_threads + assert stat.counts.runs == sim.number_of_threads + + is_ok = True + utility.test_ok(is_ok) diff --git a/opengate/tests/src/chemistry/test102_chemical_counting_actor_aggregate.py b/opengate/tests/src/chemistry/test102_chemical_counting_actor_aggregate_mt.py similarity index 100% rename from opengate/tests/src/chemistry/test102_chemical_counting_actor_aggregate.py rename to opengate/tests/src/chemistry/test102_chemical_counting_actor_aggregate_mt.py diff --git a/opengate/tests/src/geometry/test007_volumes.py b/opengate/tests/src/geometry/test007_volumes.py index 305f6e0b3..c6c3a78d2 100755 --- a/opengate/tests/src/geometry/test007_volumes.py +++ b/opengate/tests/src/geometry/test007_volumes.py @@ -316,6 +316,6 @@ def check_mat(se): sec = gate.g4_units.second c.duration = 2.5267 * sec print("-" * 80) - is_ok = utility.assert_stats(stats, stats_ref, 0.15) + is_ok = utility.assert_stats(stats, stats_ref, 0.16) utility.test_ok(is_ok) diff --git a/opengate/tests/src/misc/test069_sim_as_dict.py b/opengate/tests/src/misc/test069_sim_as_dict.py index e3eb0a869..8c0c22977 100755 --- a/opengate/tests/src/misc/test069_sim_as_dict.py +++ b/opengate/tests/src/misc/test069_sim_as_dict.py @@ -99,18 +99,31 @@ # run sim.run() - # test the file content -> NO, there are some abs filenames ... + # test the file content fn1 = paths.output / "test069" / sim.json_archive_filename print(fn1) - """fn2 = paths.output_ref / "test069" / sim.json_archive_filename - f1 = open(fn1) - j1 = json.load(f1) - f2 = open(fn2) - j2 = json.load(f2) - is_ok = j1 == j2 - print(fn1) - print(fn2) - utility.print_test(is_ok, f"Compare json gate output with reference")""" - is_ok = fn1.exists() + from opengate.serialization import load_json + + with open(fn1) as f: + dct = load_json(f) + + is_ok = "source_manager" in dct + if is_ok: + sources_dct = dct["source_manager"]["sources"] + is_ok = is_ok and ("mysource" in sources_dct) + if "mysource" in sources_dct: + mysource_dct = sources_dct["mysource"] + is_ok = is_ok and (mysource_dct["user_info"]["particle"] == "proton") + is_ok = is_ok and (mysource_dct["user_info"]["energy"]["mono"] == 230 * MeV) + + # Test reading back + sim2 = gate.Simulation() + sim2.from_dictionary(dct) + is_ok = is_ok and (sim2.source_manager.get_source("mysource").particle == "proton") + is_ok = is_ok and ( + sim2.source_manager.get_source("mysource").energy.mono == 230 * MeV + ) + + is_ok = is_ok and fn1.exists() utility.test_ok(is_ok) diff --git a/opengate/tests/src/misc/test069_sim_from_dict.py b/opengate/tests/src/misc/test069_sim_from_dict.py index 21b0647de..ee57faf68 100755 --- a/opengate/tests/src/misc/test069_sim_from_dict.py +++ b/opengate/tests/src/misc/test069_sim_from_dict.py @@ -34,6 +34,12 @@ def main(a, b, c, dependency="test069_sim_as_dict.py"): sim.volume_manager.get_volume("waterbox_with_hole").creator_volumes[0].name == "Waterbox" ) + # Assert that sources have been read back correctly + assert "mysource" in sim.source_manager.sources + assert sim.source_manager.get_source("mysource").particle == "proton" + assert ( + sim.source_manager.get_source("mysource").energy.mono == 230 * gate.g4_units.MeV + ) # If we make it until here without exception, the test is passed utility.test_ok(True) diff --git a/opengate/tests/src/source/test034_gan_phsp_linac.py b/opengate/tests/src/source/test034_gan_phsp_linac.py index 9595c7c80..df1b629f7 100755 --- a/opengate/tests/src/source/test034_gan_phsp_linac.py +++ b/opengate/tests/src/source/test034_gan_phsp_linac.py @@ -109,7 +109,7 @@ sim.run() s = sim.source_manager.get_source("gaga") - print(f"Source, nb of E<=0: {s.GetTotalSkippedEvents()}") + print(f"Source, nb of E<=0: {s.total_skipped_events}") # print results gate.exception.warning(f"Check stats") diff --git a/opengate/tests/src/source/test040_gan_phsp_pet_gan.py b/opengate/tests/src/source/test040_gan_phsp_pet_gan.py index b2532cc03..72d2d1a42 100755 --- a/opengate/tests/src/source/test040_gan_phsp_pet_gan.py +++ b/opengate/tests/src/source/test040_gan_phsp_pet_gan.py @@ -207,7 +207,7 @@ def gen_cond(n): s = sim.source_manager.get_source("gaga") else: s = sim.source_manager.get_source_mt("gaga", 0) - print(f"Source, nb of skipped particles : {s.GetTotalSkippedEvents()}") + print(f"Source, nb of skipped particles : {s.total_skipped_events}") b = gsource.total_skipped_events print(f"Source, nb of skipped particles (check) : {b}") diff --git a/opengate/tests/src/source/test053_phid_06_it_model_mt.py b/opengate/tests/src/source/test053_phid_06_it_model_mt.py index 7ceb619f0..e6e825422 100755 --- a/opengate/tests/src/source/test053_phid_06_it_model_mt.py +++ b/opengate/tests/src/source/test053_phid_06_it_model_mt.py @@ -1,5 +1,6 @@ #!/usr/bin/env python3 # -*- coding: utf-8 -*- + from test053_phid_helpers2 import * import opengate as gate diff --git a/opengate/tests/src/source/test059_tpsource_LUT_beammodel.py b/opengate/tests/src/source/test059_tpsource_LUT_beammodel.py old mode 100644 new mode 100755 index f3e9d93cc..21e6c3e41 --- a/opengate/tests/src/source/test059_tpsource_LUT_beammodel.py +++ b/opengate/tests/src/source/test059_tpsource_LUT_beammodel.py @@ -1,3 +1,6 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + import gatetools import numpy as np diff --git a/opengate/tests/src/source/test059_tpsource_MU_to_N.py b/opengate/tests/src/source/test059_tpsource_MU_to_N.py old mode 100644 new mode 100755 diff --git a/opengate/tests/src/source/test097_dynamic_voxel_source.py b/opengate/tests/src/source/test097_dynamic_voxel_source.py old mode 100644 new mode 100755 diff --git a/opengate/tests/src/source/test097_dynamic_voxel_source_half_life.py b/opengate/tests/src/source/test097_dynamic_voxel_source_half_life.py old mode 100644 new mode 100755 diff --git a/opengate/tests/src/source/test097_dynamic_voxel_source_mt.py b/opengate/tests/src/source/test097_dynamic_voxel_source_mt.py old mode 100644 new mode 100755 diff --git a/opengate/tests/src/source/test097_dynamic_voxel_source_tac.py b/opengate/tests/src/source/test097_dynamic_voxel_source_tac.py old mode 100644 new mode 100755 diff --git a/opengate/tests/src/source/test103_debug_source.py b/opengate/tests/src/source/test103_debug_source.py new file mode 100755 index 000000000..6b0764165 --- /dev/null +++ b/opengate/tests/src/source/test103_debug_source.py @@ -0,0 +1,85 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import os +from scipy.spatial.transform import Rotation +import opengate as gate +from opengate.tests import utility + +if __name__ == "__main__": + paths = utility.get_default_test_paths(__file__, "", "test103_debug_source") + + # create the simulation + sim = gate.Simulation() + sim.visu = False + sim.random_seed = "auto" + sim.output_dir = paths.output + sim.number_of_threads = 3 + if os.name == "nt": + sim.number_of_threads = 1 + sim.store_json_archive = True + + # units + m = gate.g4_units.m + cm = gate.g4_units.cm + mm = gate.g4_units.mm + MeV = gate.g4_units.MeV + + # world size + world = sim.world + world.size = [1 * m, 1 * m, 1 * m] + + # fake volume + fake = sim.add_volume("Box", "fake") + fake.size = [40 * cm, 40 * cm, 40 * cm] + fake.translation = [1 * cm, 2 * cm, 3 * cm] + fake.rotation = Rotation.from_euler("x", 10, degrees=True).as_matrix() + fake.material = "G4_AIR" + fake.color = [1, 0, 1, 1] + + # debug source + debug_source = sim.add_source("DebugSource", "debug_source") + debug_source.n = 3 + debug_source.debug_flag = True + + # add stat actor + stat = sim.add_actor("SimulationStatisticsActor", "Stats") + stat.track_types_flag = True + + # start simulation in another process + sim.run(start_new_process=True) + + # print results at the end + print(stat) + print(debug_source) + + print(f"Simulation json saved in {sim.output_dir / sim.json_archive_filename}") + + # assertions to verify MT execution and output recovery + assert stat.counts.events == debug_source.n * sim.number_of_threads + assert stat.counts.runs == sim.number_of_threads + + assert hasattr(debug_source, "debug_flag") + assert debug_source.debug_flag == True + + assert hasattr(debug_source, "debug_value") + assert debug_source.debug_value == debug_source.n * sim.number_of_threads + + # check simulation json source_manager entry + json_path = sim.output_dir / sim.json_archive_filename + assert json_path.exists() + + from opengate.serialization import load_json + + with open(json_path) as f: + dct = load_json(f) + + assert "source_manager" in dct + sources_dct = dct["source_manager"]["sources"] + assert "debug_source" in sources_dct + debug_src_dct = sources_dct["debug_source"] + assert debug_src_dct["user_info"]["n"] == [debug_source.n] + assert debug_src_dct["user_info"]["debug_flag"] == True + + is_ok = True + utility.test_ok(is_ok)