diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index 856cf8f93..bff6f8a1f 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -17,10 +17,10 @@ jobs: strategy: fail-fast: false matrix: - os: [ubuntu-latest, macos-latest, windows-latest] + os: [ubuntu-latest, macos-11, windows-latest] python-version: [3.8, 3.9, '3.10', '3.11'] exclude: - - os: macos-latest + - os: macos-11 python-version: '3.11' env: @@ -36,7 +36,7 @@ jobs: if [[ ${{ matrix.os }} == "ubuntu-latest" ]]; then export GIT_SSL_NO_VERIFY=1 git submodule update --init --recursive - elif [[ ${{ matrix.os }} == "macos-latest" ]]; then + elif [[ ${{ matrix.os }} == "macos-11" ]]; then export GIT_SSL_NO_VERIFY=1 git submodule update --init --recursive else @@ -56,7 +56,7 @@ jobs: varOS=`cat /etc/os-release | grep "VERSION=" | grep -oP '(?<=\").*?(?=\")'` varOS=($varOS) echo "release=${varOS[0]}" >> $GITHUB_OUTPUT - elif [[ ${{ matrix.os }} == "macos-latest" ]]; then + elif [[ ${{ matrix.os }} == "macos-11" ]]; then varOS=`sw_vers | grep "ProductVersion:"` varOS="${varOS#*:}" echo "release=${varOS:1}" >> $GITHUB_OUTPUT @@ -106,21 +106,21 @@ jobs: mv dist_opengate/* dist/ fi - uses: conda-incubator/setup-miniconda@v2 - if: (matrix.os == 'macos-latest') || (matrix.os == 'windows-latest') + if: (matrix.os == 'macos-11') || (matrix.os == 'windows-latest') with: auto-update-conda: true activate-environment: opengate_core python-version: ${{ matrix.python-version }} - name: Set up Homebrew - if: matrix.os == 'macos-latest' + if: matrix.os == 'macos-11' id: set-up-homebrew uses: Homebrew/actions/setup-homebrew@master - name: Create opengate_core Wheel Mac - if: matrix.os == 'macos-latest' + if: matrix.os == 'macos-11' shell: bash -l {0} run: | brew update - rm -rf /usr/local/bin/python3.11-config /usr/local/bin/2to3-3.11 /usr/local/bin/idle3.11 /usr/local/bin/pydoc3.11 /usr/local/bin/python3.11 + rm -rf /usr/local/bin/python3.1*-config /usr/local/bin/2to3-3.1* /usr/local/bin/idle3.1* /usr/local/bin/pydoc3.1* /usr/local/bin/python3.1* rm -rf /usr/local/bin/python3-config /usr/local/bin/2to3 /usr/local/bin/idle3 /usr/local/bin/pydoc3 /usr/local/bin/python3 brew install --force --verbose --overwrite \ ccache \ diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index c43f13afd..6abf66756 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -5,11 +5,11 @@ repos: - id: end-of-file-fixer - id: trailing-whitespace - repo: https://github.com/psf/black - rev: 23.11.0 + rev: 23.12.0 hooks: - id: black - repo: https://github.com/pre-commit/mirrors-clang-format - rev: v17.0.5 + rev: v17.0.6 hooks: - id: clang-format ci: diff --git a/core/opengate_core/opengate_lib/GateGANSource.cpp b/core/opengate_core/opengate_lib/GateGANSource.cpp index 8b451939f..25106edf1 100644 --- a/core/opengate_core/opengate_lib/GateGANSource.cpp +++ b/core/opengate_core/opengate_lib/GateGANSource.cpp @@ -88,7 +88,20 @@ void GateGANSource::PrepareNextRun() { // (no need to update th fSPS pos in GateGenericSource) // GateVSource::PrepareNextRun(); // FIXME remove this function ? - GateGenericSource::PrepareNextRun(); + // GateGenericSource::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 *pos = fSPS->GetPosDist(); + pos->SetCentreCoords(l.fGlobalTranslation); + + // orientation according to mother volume + auto rotation = l.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); } void GateGANSource::SetGeneratorFunction(ParticleGeneratorType &f) { diff --git a/core/opengate_core/opengate_lib/GateGenericSource.cpp b/core/opengate_core/opengate_lib/GateGenericSource.cpp index 374950c1b..6cd829c9d 100644 --- a/core/opengate_core/opengate_lib/GateGenericSource.cpp +++ b/core/opengate_core/opengate_lib/GateGenericSource.cpp @@ -13,9 +13,12 @@ #include GateGenericSource::GateGenericSource() : GateVSource() { - fNumberOfGeneratedEvents = 0; + /*fNumberOfGeneratedEvents = 0; fMaxN = 0; fActivity = 0; + fHalfLife = -1; + fLambda = -1; + */ fInitGenericIon = false; fA = 0; fZ = 0; @@ -23,8 +26,6 @@ GateGenericSource::GateGenericSource() : GateVSource() { fInitConfine = false; fWeight = -1; fWeightSigma = -1; - fHalfLife = -1; - fDecayConstant = -1; fTotalSkippedEvents = 0; fCurrentSkippedEvents = 0; fTotalZeroEvents = 0; @@ -34,6 +35,7 @@ GateGenericSource::GateGenericSource() : GateVSource() { fParticleDefinition = nullptr; fEffectiveEventTime = -1; fEffectiveEventTime = -1; + fforceRotation = false; } GateGenericSource::~GateGenericSource() { @@ -75,18 +77,20 @@ void GateGenericSource::InitializeUserInfo(py::dict &user_info) { CreateSPS(); // get user info about activity or nb of events + /* fMaxN = DictGetInt(user_info, "n"); fActivity = DictGetDouble(user_info, "activity"); fInitialActivity = fActivity; // half life ? fHalfLife = DictGetDouble(user_info, "half_life"); - fDecayConstant = log(2) / fHalfLife; - fUserParticleLifeTime = DictGetDouble(user_info, "user_particle_life_time"); + fLambda = log(2) / fHalfLife; + */ // weight fWeight = DictGetDouble(user_info, "weight"); fWeightSigma = DictGetDouble(user_info, "weight_sigma"); + fUserParticleLifeTime = DictGetDouble(user_info, "user_particle_life_time"); // get the user info for the particle InitializeParticle(user_info); @@ -103,14 +107,14 @@ void GateGenericSource::InitializeUserInfo(py::dict &user_info) { fCurrentSkippedEvents = 0; fTotalSkippedEvents = 0; fEffectiveEventTime = -1; + + fforceRotation = DictGetDouble(user_info, "force_rotation"); } void GateGenericSource::UpdateActivity(double time) { if (!fTAC_Times.empty()) return UpdateActivityWithTAC(time); - if (fHalfLife <= 0) - return; - fActivity = fInitialActivity * exp(-fDecayConstant * (time - fStartTime)); + GateVSource::UpdateActivity(time); } void GateGenericSource::UpdateActivityWithTAC(double time) { @@ -159,14 +163,15 @@ double GateGenericSource::PrepareNextTime(double current_simulation_time) { return -1; // get next time according to current fActivity - double next_time = - fEffectiveEventTime - log(G4UniformRand()) * (1.0 / fActivity); + double next_time = CalcNextTime(fEffectiveEventTime); if (next_time >= fEndTime) return -1; return next_time; } // check according to t MaxN + // std::cout<= fMaxN) { return -1; } @@ -180,15 +185,30 @@ void GateGenericSource::PrepareNextRun() { // This global transformation is given to the SPS that will // generate particles in the correct coordinate system auto &l = fThreadLocalData.Get(); + // auto user_info_pos = py::dict(puser_info["position"]); + // auto pos_init = DictGetG4ThreeVector(user_info_pos, "translation"); auto *pos = fSPS->GetPosDist(); pos->SetCentreCoords(l.fGlobalTranslation); // orientation according to mother volume auto rotation = l.fGlobalRotation; - G4ThreeVector r1(rotation(0, 0), rotation(0, 1), rotation(0, 2)); - G4ThreeVector r2(rotation(1, 0), rotation(1, 1), rotation(1, 2)); + 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 *ang = fSPS->GetAngDist(); + + if (fangType == "momentum" && fforceRotation) { + auto new_d = rotation * fInitializeMomentum; + ang->SetParticleMomentumDirection(new_d); + } + if (fangType == "focused" && fforceRotation) { + auto vec_f = fInitiliazeFocusPoint - fInitTranslation; + auto rot_f = rotation * vec_f; + auto new_f = rot_f + l.fGlobalTranslation; + ang->SetFocusPoint(new_f); + } } void GateGenericSource::UpdateEffectiveEventTime( @@ -226,7 +246,6 @@ void GateGenericSource::GeneratePrimaries(G4Event *event, // (acceptance angle is included) fSPS->SetParticleTime(current_simulation_time); fSPS->GeneratePrimaryVertex(event); - // update the time according to skipped events fEffectiveEventTime = current_simulation_time; auto &l = fThreadLocalDataAA.Get(); @@ -300,6 +319,7 @@ void GateGenericSource::InitializePosition(py::dict puser_info) { std::vector l = {"sphere", "point", "box", "disc", "cylinder"}; CheckIsIn(pos_type, l); auto translation = DictGetG4ThreeVector(user_info, "translation"); + fInitTranslation = translation; if (pos_type == "point") { pos->SetPosDisType("Point"); } @@ -365,6 +385,7 @@ void GateGenericSource::InitializeDirection(py::dict puser_info) { auto user_info = py::dict(puser_info["direction"]); auto *ang = fSPS->GetAngDist(); auto ang_type = DictGetStr(user_info, "type"); + fangType = ang_type; std::vector ll = {"iso", "momentum", "focused", "beam2d"}; // FIXME check on py side ? CheckIsIn(ang_type, ll); @@ -382,11 +403,13 @@ void GateGenericSource::InitializeDirection(py::dict puser_info) { if (ang_type == "momentum") { ang->SetAngDistType("planar"); // FIXME really ?? auto d = DictGetG4ThreeVector(user_info, "momentum"); + fInitializeMomentum = d; ang->SetParticleMomentumDirection(d); } if (ang_type == "focused") { ang->SetAngDistType("focused"); auto f = DictGetG4ThreeVector(user_info, "focus_point"); + fInitiliazeFocusPoint = f; ang->SetFocusPoint(f); } if (ang_type == "beam2d") { diff --git a/core/opengate_core/opengate_lib/GateGenericSource.h b/core/opengate_core/opengate_lib/GateGenericSource.h index 3d1e65b72..5034481d6 100644 --- a/core/opengate_core/opengate_lib/GateGenericSource.h +++ b/core/opengate_core/opengate_lib/GateGenericSource.h @@ -34,7 +34,7 @@ class GateGenericSource : public GateVSource { /// Current number of simulated events in this source /// (do not include skipped events) - unsigned long fNumberOfGeneratedEvents; + // unsigned long fNumberOfGeneratedEvents; /// Count the number of skipped events /// (e.g. Acceptance Angle or in GANSource) @@ -51,16 +51,22 @@ class GateGenericSource : public GateVSource { const std::vector &activities); protected: - unsigned long fMaxN; - // We cannot not use a std::unique_ptr - // (or maybe by controlling the deletion during the CleanWorkerThread ?) + // unsigned long fMaxN; + // We cannot not use a std::unique_ptr + // (or maybe by controlling the deletion during the CleanWorkerThread ?) GateSingleParticleSource *fSPS; - + /* double fActivity; double fInitialActivity; double fHalfLife; - double fDecayConstant; + double fLambda; + */ + G4ParticleDefinition *fParticleDefinition; + G4ThreeVector fInitializeMomentum; + G4ThreeVector fInitiliazeFocusPoint; + G4ThreeVector fInitTranslation; + G4String fangType; double fEffectiveEventTime; double fUserParticleLifeTime; @@ -78,6 +84,10 @@ class GateGenericSource : public GateVSource { double fWeight; double fWeightSigma; + // Force the rotation of momentum and focal point to follow rotation of the + // source, eg: needed for motion actor + bool fforceRotation; + // angular acceptance management struct threadLocalT { GateAcceptanceAngleTesterManager *fAAManager; @@ -106,7 +116,7 @@ class GateGenericSource : public GateVSource { virtual void InitializeEnergy(py::dict user_info); - virtual void UpdateActivity(double time); + virtual void UpdateActivity(double time) override; void UpdateEffectiveEventTime(double current_simulation_time, unsigned long skipped_particle); diff --git a/core/opengate_core/opengate_lib/GateMotionVolumeActor.cpp b/core/opengate_core/opengate_lib/GateMotionVolumeActor.cpp index 220f1c48a..bdd5e2363 100644 --- a/core/opengate_core/opengate_lib/GateMotionVolumeActor.cpp +++ b/core/opengate_core/opengate_lib/GateMotionVolumeActor.cpp @@ -63,7 +63,7 @@ void GateMotionVolumeActor::MoveGeometry(int run_id) { rot->set(r.rep3x3()); // close the geometry manager - gm->CloseGeometry(false, false, pv); + gm->CloseGeometry(true, false, pv); // G4RunManager::GetRunManager()->GeometryHasBeenModified(true); } diff --git a/core/opengate_core/opengate_lib/GatePhaseSpaceSource.cpp b/core/opengate_core/opengate_lib/GatePhaseSpaceSource.cpp index 3b0922c12..de2b90931 100644 --- a/core/opengate_core/opengate_lib/GatePhaseSpaceSource.cpp +++ b/core/opengate_core/opengate_lib/GatePhaseSpaceSource.cpp @@ -14,7 +14,8 @@ GatePhaseSpaceSource::GatePhaseSpaceSource() : GateVSource() { fCharge = 0; fMass = 0; - fMaxN = 0; + fCurrentBatchSize = 0; + // fMaxN = 0; fGlobalFag = false; } @@ -31,7 +32,7 @@ void GatePhaseSpaceSource::InitializeUserInfo(py::dict &user_info) { GateVSource::InitializeUserInfo(user_info); // Number of events to generate - fMaxN = DictGetInt(user_info, "n"); + // fMaxN = DictGetInt(user_info, "n"); // global (world) or local (mother volume) coordinate system fGlobalFag = DictGetBool(user_info, "global_flag"); @@ -54,6 +55,12 @@ void GatePhaseSpaceSource::InitializeUserInfo(py::dict &user_info) { l.fNumberOfGeneratedEvents = 0; l.fCurrentIndex = 0; l.fCurrentBatchSize = 0; + + l.fgenerate_until_next_primary = + DictGetBool(user_info, "generate_until_next_primary"); + l.fprimary_lower_energy_threshold = + DictGetDouble(user_info, "primary_lower_energy_threshold"); + l.fprimary_PDGCode = DictGetInt(user_info, "primary_PDGCode"); } void GatePhaseSpaceSource::PrepareNextRun() { @@ -63,6 +70,19 @@ void GatePhaseSpaceSource::PrepareNextRun() { double GatePhaseSpaceSource::PrepareNextTime(double current_simulation_time) { // check according to t MaxN + + UpdateActivity(current_simulation_time); + if (fMaxN <= 0) { + if (current_simulation_time < fStartTime) + return fStartTime; + if (current_simulation_time >= fEndTime) + return -1; + + double next_time = CalcNextTime(current_simulation_time); + if (next_time >= fEndTime) + return -1; + return next_time; + } auto &l = fThreadLocalDataPhsp.Get(); if (l.fNumberOfGeneratedEvents >= fMaxN) { return -1; @@ -93,19 +113,51 @@ void GatePhaseSpaceSource::GenerateBatchOfParticles() { void GatePhaseSpaceSource::GeneratePrimaries(G4Event *event, double current_simulation_time) { auto &l = 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 (l.fgenerate_until_next_primary) { + int num_primaries = 0; + while (num_primaries <= 2) { + // If batch is empty, we generate some particles + if (l.fCurrentIndex >= l.fCurrentBatchSize) + GenerateBatchOfParticles(); + + // check if next particle is primary + if (ParticleIsPrimary()) + num_primaries++; + + // don't generate the second primary + if (num_primaries < 2) { + // Go + GenerateOnePrimary(event, current_simulation_time); + + // update the index; + l.fCurrentIndex++; + + // // update the number of generated event + // fNumberOfGeneratedEvents++; + } else + break; + } + // update the number of generated event + l.fNumberOfGeneratedEvents++; + } - // If batch is empty, we generate some particles - if (l.fCurrentIndex >= l.fCurrentBatchSize) - GenerateBatchOfParticles(); + else { + // If batch is empty, we generate some particles + if (l.fCurrentIndex >= l.fCurrentBatchSize) + GenerateBatchOfParticles(); - // Go - GenerateOnePrimary(event, current_simulation_time); + // Go + GenerateOnePrimary(event, current_simulation_time); - // update the index; - l.fCurrentIndex++; + // update the index; + l.fCurrentIndex++; - // update the number of generated event - l.fNumberOfGeneratedEvents++; + // update the number of generated event + l.fNumberOfGeneratedEvents++; + } } void GatePhaseSpaceSource::GenerateOnePrimary(G4Event *event, @@ -221,3 +273,24 @@ void GatePhaseSpaceSource::SetDirectionZBatch( auto &l = fThreadLocalDataPhsp.Get(); l.fDirectionZ = PyBindGetVector(fDirectionZ); } + +bool GatePhaseSpaceSource::ParticleIsPrimary() { + auto &l = fThreadLocalDataPhsp.Get(); + // check if particle is primary + bool is_primary = false; + + // if PDGCode exists in file + if ((l.fPDGCode[l.fCurrentIndex] != 0) && (l.fprimary_PDGCode != 0)) { + if ((l.fprimary_PDGCode == l.fPDGCode[l.fCurrentIndex]) && + (l.fprimary_lower_energy_threshold <= l.fEnergy[l.fCurrentIndex])) { + is_primary = true; + } + } else { + G4Exception("GatePhaseSpaceSource::ParticleIsPrimary", "Error", + FatalException, "Particle type not defined in file"); + std::cout << "ERROR: PDGCode not defined in file. Aborting." << std::endl; + exit(1); + } + + return is_primary; +} diff --git a/core/opengate_core/opengate_lib/GatePhaseSpaceSource.h b/core/opengate_core/opengate_lib/GatePhaseSpaceSource.h index 4a4b0c0ac..24e6edb22 100644 --- a/core/opengate_core/opengate_lib/GatePhaseSpaceSource.h +++ b/core/opengate_core/opengate_lib/GatePhaseSpaceSource.h @@ -45,6 +45,10 @@ class GatePhaseSpaceSource : public GateVSource { void SetGeneratorFunction(ParticleGeneratorType &f) const; + bool ParticleIsPrimary(); + + // virtual void SetGeneratorInfo(py::dict &user_info); + void GenerateBatchOfParticles(); G4ParticleDefinition *fParticleDefinition; @@ -54,7 +58,9 @@ class GatePhaseSpaceSource : public GateVSource { bool fGlobalFag; bool fUseParticleTypeFromFile; - unsigned long fMaxN; + // unsigned long fMaxN; + long fNumberOfGeneratedEvents; + size_t fCurrentBatchSize; void SetPDGCodeBatch(const py::array_t &fPDGCode) const; @@ -77,6 +83,10 @@ class GatePhaseSpaceSource : public GateVSource { // For MT, all threads local variables are gathered here struct threadLocalTPhsp { + bool fgenerate_until_next_primary; + int fprimary_PDGCode; + double fprimary_lower_energy_threshold; + ParticleGeneratorType fGenerator; unsigned long fNumberOfGeneratedEvents; size_t fCurrentIndex; diff --git a/core/opengate_core/opengate_lib/GateVSource.cpp b/core/opengate_core/opengate_lib/GateVSource.cpp index ba139cb29..e62826075 100644 --- a/core/opengate_core/opengate_lib/GateVSource.cpp +++ b/core/opengate_core/opengate_lib/GateVSource.cpp @@ -7,6 +7,7 @@ #include "GateVSource.h" #include "G4PhysicalVolumeStore.hh" +#include "G4RandomTools.hh" #include "GateHelpers.h" #include "GateHelpersDict.h" #include "GateHelpersGeometry.h" @@ -20,6 +21,11 @@ GateVSource::GateVSource() { fMother = ""; fLocalTranslation = G4ThreeVector(); fLocalRotation = G4RotationMatrix(); + fNumberOfGeneratedEvents = 0; + fMaxN = 0; + fActivity = 0; + fHalfLife = -1; + fDecayConstant = -1; } GateVSource::~GateVSource() {} @@ -30,6 +36,27 @@ void GateVSource::InitializeUserInfo(py::dict &user_info) { fStartTime = DictGetDouble(user_info, "start_time"); fEndTime = DictGetDouble(user_info, "end_time"); fMother = DictGetStr(user_info, "mother"); + + // get user info about activity or nb of events + fMaxN = DictGetInt(user_info, "n"); + fActivity = DictGetDouble(user_info, "activity"); + fInitialActivity = fActivity; + + // half life ? + fHalfLife = DictGetDouble(user_info, "half_life"); + fDecayConstant = log(2) / fHalfLife; +} + +void GateVSource::UpdateActivity(double time) { + if (fHalfLife <= 0) + return; + fActivity = fInitialActivity * exp(-fDecayConstant * (time - fStartTime)); +} + +double GateVSource::CalcNextTime(double current_simulation_time) { + double next_time = + current_simulation_time - log(G4UniformRand()) * (1.0 / fActivity); + return next_time; } void GateVSource::PrepareNextRun() { SetOrientationAccordingToMotherVolume(); } diff --git a/core/opengate_core/opengate_lib/GateVSource.h b/core/opengate_core/opengate_lib/GateVSource.h index ae2631963..bc8a4bc92 100644 --- a/core/opengate_core/opengate_lib/GateVSource.h +++ b/core/opengate_core/opengate_lib/GateVSource.h @@ -29,6 +29,10 @@ class GateVSource { // Called at initialisation to set the source properties from a single dict virtual void InitializeUserInfo(py::dict &user_info); + virtual void UpdateActivity(double time); + + double CalcNextTime(double current_simulation_time); + virtual void PrepareNextRun(); virtual double PrepareNextTime(double current_simulation_time); @@ -40,6 +44,8 @@ class GateVSource { std::string fName; double fStartTime; double fEndTime; + unsigned long fNumberOfGeneratedEvents; + std::string fMother; std::vector fTranslations; std::vector fRotations; @@ -47,6 +53,16 @@ class GateVSource { G4ThreeVector fLocalTranslation; G4RotationMatrix fLocalRotation; + G4ThreeVector fGlobalTranslation; + G4RotationMatrix fGlobalRotation; + +protected: + unsigned long fMaxN; + double fActivity; + double fInitialActivity; + double fHalfLife; + double fDecayConstant; + struct threadLocalT { G4ThreeVector fGlobalTranslation; G4RotationMatrix fGlobalRotation; diff --git a/core/opengate_core/opengate_lib/GateVoxelsSource.cpp b/core/opengate_core/opengate_lib/GateVoxelsSource.cpp index 3551d27d2..8922e43f7 100644 --- a/core/opengate_core/opengate_lib/GateVoxelsSource.cpp +++ b/core/opengate_core/opengate_lib/GateVoxelsSource.cpp @@ -17,9 +17,23 @@ GateVoxelsSource::GateVoxelsSource() : GateGenericSource() { GateVoxelsSource::~GateVoxelsSource() {} void GateVoxelsSource::PrepareNextRun() { - GateGenericSource::PrepareNextRun(); - // rotation and translation to apply, according to mother volume + // GateGenericSource::PrepareNextRun(); + // rotation and translation to apply, according to mother volume + GateVSource::PrepareNextRun(); + // This global transformation is given to the SPS that will + // generate particles in the correct coordinate system auto &l = fThreadLocalData.Get(); + auto *pos = fSPS->GetPosDist(); + pos->SetCentreCoords(l.fGlobalTranslation); + + // orientation according to mother volume + auto rotation = l.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; // the direction is 'isotropic' so we don't care about rotating the direction. diff --git a/docs/source/developer_guide.md b/docs/source/developer_guide.md index 64e8ad65e..18b37dee2 100644 --- a/docs/source/developer_guide.md +++ b/docs/source/developer_guide.md @@ -475,6 +475,57 @@ Below are a list of hints (compared to boost-python). - Pure virtual need a trampoline class - Python debug: `python -q -X faulthandler` + +### Memory managment between Python and Geant4 - Why is there a segmentation fault? + +Code and memory (objects) can be on the python and/or C++(Geant4)-side. +To avoid memory leakage, the G4 objects need to be deleted at some point, either by the C++-side, meaning G4, or by the garbage collector on the python-side. This requires some understanding of how lifetime is managed in G4. GATE 10 uses the G4RunManager to initialize, run, and deconstruct a simulation. The G4RunManager's destructor triggers a nested sequence of calls to the destructors of many objects the run manager handles, e.g. geometry and physics. These objects, if they are created on the python-side, e.g. because the Construct() method is implemented in python, should never be deleted on the python-side because the G4RunManager is responsible for deletion and segfaults occur if the objects exist now longer at the time of supposed destruction. In this case, the no `py::nodelete` option in pybind11 is the correct way. + +​ +​There is also the python-side concerning objects which are deleted by the G4RunManager: If python stores a reference to the object, and it generally should (see below), this reference will not find the object anymore once the G4RunManager has been deleted. Therefore, we have implemented a mechanism which makes sure that references to G4 objects are unset before the G4RunManager is garbage collected (and thus its destructor is called). This is done via the close() and release_g4_references() methods combined with a `with SimulationEngine as se` context clause. +​ + +​The situation is somewhat different for objects that are not automatically deleted by the Geant4 session. This concerns objects which a G4 user would manually destruct and delete. The "user", from a G4 perspective, is the python-side of GATE, specifically the SimulationEngine, which creates and controls the G4RunManager. The objects in question should be destroyed on the python side, and in fact they (normally) are via garbage collection. To this end, they should **not** be bound with the `py::nodelete` option - which would prevent garbage collection. ​ + +In any case, G4 objects created on the python-side should not be destroyed (garbage collected) too early, i.e. not before the end of the G4 simulation, to avoid segfaults when G4 tries to a find the object to which it holds a pointer. It is important to note here that garbage collection in python is strongly tied to reference counting, meaning, objects may be garbage collected as soon as there are no more references pointing to them. Or in other words: you can prevent garbage collection by keeping a reference. In practice: If you create a G4 object that is required to persist for the duration of the simulation and you create this object in a local scope, e.g. in a function, you need to make sure to store a reference somewhere so that it lives beyond the local scope. Otherwise, when the function goes out of scope, the local reference no longer exists and the G4 object may be garbage collected. This is the reason why GATE 10 stores references to G4 objects in class attributes, either as plane reference, or via lists or dictionaries. A nasty detail: in case of a G4 object that is only referenced locally (implementation error), the moment when the segfault occurs might vary because garbage collection in python is typically scheduled (for performance reasons), meaning objects are collected any time after their last reference is released, but not necessarily at that moment. This can cause a seeming randomness in the segfaults. + +​ +​So, after a long explanation, the key points are: + +#. check the nature of your G4 object + +​#. if it is designed to be deleted by the G4RunManager, add the `py::no_delete` option in the pybind11 binding + +​#. In any case, make sure to store a non-local persistent reference on the python-side, ideally as a class attribute + +#. add a “release statement” in the release_g4_references() method for each G4 reference (strictly only to RunManager-handled objects, but it does not hurt for others) in the class, and make sure the release_g4_references() method is called by the class’s close() method. + +An example may be the implementation of the Tesselated Volume (STL). + +Pybindings in pyG4TriangularFacet.cpp: + +``py::class_(m, "G4TriangularFacet")`` +will cause a segmentation fault + +``py::class_>(m, "G4TriangularFacet")`` +will work as python is instructed not to delete the object. + +On the python-side in geometry/solids.py: + +```python +def __init__(self, *args, **kwargs): + super().__init__(*args, **kwargs) + + self.g4_solid = None + self.g4_tessellated_solid = None + self.facetArray = None +``` + +All data which is sent to Geant4 is included as a variable to the class. As such, it is only deleted at the end of the simulation, not when the function ends. +Which, would cause a segmentation fault. + + + ### ITK - The following [issue](https://github.com/OpenGATE/opengate/issues/216) occured in the past and was solved by updating ITK: diff --git a/docs/source/user_guide_2_2_sources.md b/docs/source/user_guide_2_2_sources.md index 899680513..c8fb7bd9d 100644 --- a/docs/source/user_guide_2_2_sources.md +++ b/docs/source/user_guide_2_2_sources.md @@ -133,7 +133,7 @@ Like all objects, by default, the source is located according to the coordinate ### Phase-Space sources -A phase-space source read particles properties (position, direction, energy, etc.) from a root file and use them as events. Here is an example: +A phase-space source reads particles properties (position, direction, energy, etc.) from a root file and use them as events. Here is an example: ```python source = sim.add_source("PhaseSpaceSource", "phsp_source") @@ -149,15 +149,56 @@ source.n = 20000 In that case, the key "PrePositionLocal" in the root tree file will be used to define the position of all generated particles. The flag "global_flag" is False so the position will be relative to the mother volume (the plane here) ; otherwise, position is considered as global (in the world coordinate system). + Limitation: the particle timestamps is NOT read from the phsp and not considered (yet) -The particle type can be set by ```source.particle = "proton"``` option (all generated particles will be for example proton), or read in the phsp file by using the PDGCode: +The particle type can be set by ```source.particle = "proton"``` option. Using this option all generated particles will be for example protons, overriding the particle type specified in the phase space. +Alternatively, by setting ```source.particle = None``` the particle type is read from the phase space file using the PDGCode. ```source.PDGCode_key = PDGCode``` specifies the name of the entry in the phase space file. +Full listing: ```python source.PDGCode_key = "PDGCode" source.particle = None ``` +The PDGCode is defined by the particle data group (see https://pdg.lbl.gov/2023/web/viewer.html?file=../reviews/rpp2022-rev-monte-carlo-numbering.pdf). +Here is a short overview of common particle types and its corresponding PDG Code +``` +proton: 2212 +neutron: 2211 +electron: 11 +gamma: 22 +carbon ion C12: 1000060120 +``` + +The naming of the phsp file entries generated by e.g. a GATE phase space actor changed over time, most notably from GATE v9 to GATE v10. +Setting ```source.position_key = "PrePositionLocal"``` will cause the phsp source to look for particle positions in ```PrePositionLocal_X, PrePositionLocal_Y, PrePositionLocal_Z```. +```source.direction_key = "PreDirectionLocal"``` will do the corresponding for the particle direction vector components in ```PreDirectionLocal_X, PreDirectionLocal_y, PreDirectionLocal_Z```. + +Alternatively, it is possible to directly set the individual components: +``` +source.position_key = None"PrePositionLocal" +source.position_key_x = Position_X" +source.position_key_y = Position_X +source.position_key_z = Position_X +source.direction_key = None +source.direction_key_x = Direction_X +source.direction_key_y = Direction_X +source.direction_key_z = Direction_X +``` + +In case of simulating a realistic particle mix, for example the output after a linac, a phsp file could contain a mixture of particles. +Typically, one would be interested in simulating a given number of primary particles (e.g. protons), simulating, but not counting as secondary particles (e.g. secondary electrons) in the number of particles to simulate. +This can be acchieved by setting ```generate_until_next_primary = True```. Furthermore, the PDG code of the primary particle needs to be specified, as well as a +lower energy threshold (in order to identify secondary particles of the same type as the primary particle). +The example below will consider protons above 90 MeV as primary particles. Every primary particle found in the phsp file will increase the number of particles simulated counter, while secondary particles (e.g. all other particles in the phsp file) +will be simulated, but not be considered in (e.g. not increasing) the number of particles simulated. +``` +source.generate_until_next_primary = True +source.primary_lower_energy_threshold = 90.0 * MeV +source.primary_PDGCode = 2212 +``` + For multithread: you need to indicate the ```entry_start``` for all threads, as an array, so that each thread starts in the phsp file at a different position. This done for example as follows (see ```test019_linac_phsp_source_MT.py```). Warning, if the phsp reach its end, it will cycle and start back at the beginning. ```python @@ -166,7 +207,7 @@ nb_of_threads = 4 source.entry_start = [total_nb_of_particle * p for p in range(nb_of_threads)] ``` -See all test019 as examples. +See all test019 and test060 as examples. ### GAN sources (Generative Adversarial Network) diff --git a/opengate/contrib/tps/treatmentPlanPhsSource.py b/opengate/contrib/tps/treatmentPlanPhsSource.py new file mode 100644 index 000000000..9f6feefc9 --- /dev/null +++ b/opengate/contrib/tps/treatmentPlanPhsSource.py @@ -0,0 +1,257 @@ +import numpy as np +from scipy.spatial.transform import Rotation +import opengate as gate + +from opengate.contrib.tps.ionbeamtherapy import * + +import os + +# TreatmentPlanSource + + +class TreatmentPlanPhsSource(TreatmentPlanSource): + def __init__(self, name, sim): + self.name = name + self.rotation = Rotation.identity() + self.translation = [0, 0, 0] + self.spots = None + self.phaseSpaceFolder = "" + self.phaseSpaceList_file_name = "" + self.phaseSpaceList = {} + self.position_key_x = "PrePositionLocal_X" + self.position_key_y = "PrePositionLocal_Y" + self.position_key_z = "PrePositionLocal_Z" + self.direction_key_x = "PreDirectionLocal_X" + self.direction_key_y = "PreDirectionLocal_Y" + self.direction_key_z = "PreDirectionLocal_Z" + self.energy_key = "KineticEnergy" + self.weight_key = "Weight" + self.PDGCode_key = "PDGCode" + self.generate_until_next_primary = False + self.primary_lower_energy_threshold = 0 + self.primary_PDGCode = 0 + self.n_sim = 0 + self.sim = sim # simulation obj to which we want to add the tpPhS source + self.distance_source_to_isocenter = None + # SMX to Isocenter distance + self.distance_stearmag_to_isocenter_x = None + # SMY to Isocenter distance + self.distance_stearmag_to_isocenter_y = None + self.batch_size = None + + def __del__(self): + pass + + def set_distance_source_to_isocenter(self, distance): + self.distance_source_to_isocenter = distance + self.d_nozzle_to_iso = distance + + def set_distance_stearmag_to_isocenter(self, distance_x, distance_y): + self.d_stearMag_to_iso_x = distance_x + self.d_stearMag_to_iso_y = distance_y + self.distance_stearmag_to_isocenter_x = distance_x + self.distance_stearmag_to_isocenter_y = distance_y + + def set_phaseSpaceList_file_name(self, file_name): + self.phaseSpaceList_file_name = str(file_name) + + def set_phaseSpaceFolder(self, folder_name): + self.phaseSpaceFolder = str(folder_name) + + def initialize_tpPhssource(self): + if self.batch_size is None: + self.batch_size = 30000 + # verify that all necessary parameters are set + self.verify_necessary_parameters_are_set() + # read in the phase space list + self.phaseSpaceList = self.read_list_of_Phs( + self.phaseSpaceList_file_name, self.phaseSpaceFolder + ) + # verify the phase space list + self.verify_phs_files_exist(self.phaseSpaceList) + spots_array = self.spots + sim = self.sim + nSim = self.n_sim + + # mapping factors between iso center plane and nozzle plane (due to steering magnets) + cal_proportion_factor = ( + lambda d_magnet_iso: 1 + if (d_magnet_iso == float("inf")) + else (d_magnet_iso - self.d_nozzle_to_iso) / d_magnet_iso + ) + self.proportion_factor_x = cal_proportion_factor(self.d_stearMag_to_iso_x) + self.proportion_factor_y = cal_proportion_factor(self.d_stearMag_to_iso_y) + + tot_sim_particles = 0 + # initialize a pencil beam for each spot + for i, spot in enumerate(spots_array): + # simulate a fraction of the beam particles for this spot + nspot = np.round(spot.beamFraction * nSim) + if nspot == 0: + continue + tot_sim_particles += nspot + # create a source + source = sim.add_source("PhaseSpaceSource", f"{self.name}_spot_{i}") + + # set energy + # find corresponding phase space file + if self.phaseSpaceList.get(spot.energy) is not None: + source.phsp_file = self.phaseSpaceList.get(spot.energy) + else: + print( + "ERROR in TreatmentPlanPhsSource: Energy requested from plan file does not exist. Aborting." + ) + print("Requested energy was: ", spot.energy) + exit(-1) + + # set keys of phase space file to use + source.position_key_x = self.position_key_x + source.position_key_y = self.position_key_y + source.position_key_z = self.position_key_z + source.direction_key_x = self.direction_key_x + source.direction_key_y = self.direction_key_y + source.direction_key_z = self.direction_key_z + source.energy_key = self.energy_key + source.weight_key = self.weight_key + source.PDGCode_key = self.PDGCode_key + + if self.batch_size is not None: + source.batch_size = self.batch_size + else: + source.batch_size = 30000 + + # POSITION: + source.translate_position = True + source.position.translation = self._get_pbs_position(spot) + print("source.position.translation: ", source.position.translation) + + # ROTATION: + source.rotate_direction = True + source.position.rotation = self._get_pbs_rotation(spot) + + # add weight + source.n = nspot + + # allow the possibility to count primaries + source.generate_until_next_primary = self.generate_until_next_primary + source.primary_lower_energy_threshold = self.primary_lower_energy_threshold + source.primary_PDGCode = self.primary_PDGCode + + self.actual_sim_particles = tot_sim_particles + + def verify_phs_files_exist(self, phs_dict): + """Check if all the files in the dictionary exist. + Returns True if all the files exist, False otherwise + If one file does not exist, the program is aborted.""" + + for key, phs_file in phs_dict.items(): + # print( + # "key: ", + # key, + # " phs_file: ", + # phs_file, + # " type: ", + # type(phs_file), + # ) + + if not os.path.isfile(phs_file): + print( + "ERROR in ThreatmenPlanPhsSource: File {} does not exist".format( + phs_file + ) + ) + print("Error: File in Phase space dictionary does not exist. Aborting.") + exit(-1) + return True + + def read_list_of_Phs(self, file_name: str, path_to_phsp=""): + """Read a list of Phs from a file. + + Parameters + ---------- + file_name : str + The name of the file to read from. + File needs to have at least two columns + First column is the energy in MeV which is the energy label + Second column is the phase space file name + Returns + ------- + dictionary of Phs""" + + # prepend the path to the phase space files + file_name = path_to_phsp + "/" + file_name + + phs_dict = {} + try: + input_arr = np.genfromtxt( + file_name, + delimiter="\t", + comments="#", + usecols=(0, 1), + dtype=None, + encoding=None, + ) + + # convert to dictionary + if input_arr.shape == (0,): + print( + "Error in TreatmentPlanPhsSource: No data found in file: ", + file_name, + " Aborting.", + ) + exit(-1) + if input_arr.ndim == 0: + # only single line read, convert to array + input_arr = np.array([input_arr]) + + # normal case, multiple lines + if path_to_phsp != "": + phs_dict = { + float(i[0]): str(path_to_phsp + "/" + (i[1])) for i in input_arr + } + else: + phs_dict = {float(i[0]): str(i[1]) for i in input_arr} + # print("phs_dict read: ", phs_dict) + except Exception as e: + print( + "Error in TreatmentPlanPhsSource: could not read the phase space file list. Aborting." + ) + print("The error was: ", e) + exit(-1) + if len(phs_dict) == 0: + print( + "Error in TreatmentPlanPhsSource: the phase space file list is empty. Aborting." + ) + exit(-1) + return phs_dict + + def verify_necessary_parameters_are_set(self): + """This function checks whether all necessary parameters were set. + E.g. not None + It does not check sensibility""" + if self.phaseSpaceList_file_name is None: + print( + "Error in TreatmentPlanPhsSource: phaseSpaceList_file_name is None. Aborting." + ) + exit(-1) + if self.distance_source_to_isocenter is None: + print( + "Error in TreatmentPlanPhsSource: distance_source_to_isocenter is None. Aborting." + ) + exit(-1) + if self.distance_stearmag_to_isocenter_x is None: + print( + "Error in TreatmentPlanPhsSource: distance_stearmag_to_isocenter_x is None. Aborting." + ) + exit(-1) + if self.distance_stearmag_to_isocenter_y is None: + print( + "Error in TreatmentPlanPhsSource: distance_stearmag_to_isocenter_y is None. Aborting." + ) + exit(-1) + if self.batch_size is None: + print("Error in TreatmentPlanPhsSource: batch_size is None. Aborting.") + exit(-1) + if self.spots is None: + print("Error in TreatmentPlanPhsSource: No spots have been set. Aborting.") + exit(-1) diff --git a/opengate/engines.py b/opengate/engines.py index 4a03947c7..7fbea3432 100644 --- a/opengate/engines.py +++ b/opengate/engines.py @@ -1245,7 +1245,7 @@ def initialize_g4_verbose(self): else: # no Geant4 output ui = UIsessionSilent() - if self.simulation.g4_verbose_level_tracking != -1: + if self.simulation.g4_verbose_level_tracking >= 0: self.simulation.add_g4_command_after_init( f"/tracking/verbose {self.simulation.g4_verbose_level_tracking}" ) diff --git a/opengate/sources/generic.py b/opengate/sources/generic.py index 0d384cbfd..71cbea957 100644 --- a/opengate/sources/generic.py +++ b/opengate/sources/generic.py @@ -219,6 +219,9 @@ def set_default_user_info(user_info): user_info.mother = __world_name__ user_info.start_time = None user_info.end_time = None + user_info.n = 0 + user_info.activity = 0 + user_info.half_life = -1 # negative value is no half_life def __init__(self, user_info): # type_name MUST be defined in class that inherit from SourceBase @@ -305,14 +308,12 @@ def set_default_user_info(user_info): # initial user info user_info.particle = "gamma" user_info.ion = Box() - user_info.n = 0 - user_info.activity = 0 user_info.weight = -1 user_info.weight_sigma = -1 - user_info.half_life = -1 # negative value is no half_life user_info.user_particle_life_time = -1 # negative means : by default user_info.tac_times = None user_info.tac_activities = None + user_info.force_rotation = False # ion user_info.ion = Box() user_info.ion.Z = 0 # Z: Atomic Number @@ -497,7 +498,6 @@ class TemplateSource(SourceBase): def set_default_user_info(user_info): SourceBase.set_default_user_info(user_info) # initial user info - user_info.n = 0 user_info.float_value = None user_info.vector_value = None diff --git a/opengate/sources/phspsources.py b/opengate/sources/phspsources.py index 50301d56c..2c1e49fde 100644 --- a/opengate/sources/phspsources.py +++ b/opengate/sources/phspsources.py @@ -38,6 +38,8 @@ def initialize(self, user_info): def read_phsp_and_keys(self): # convert str like 1e5 to int self.user_info.batch_size = int(float(self.user_info.batch_size)) + if self.user_info.batch_size < 1: + gate.fatal("PhaseSpaceSourceGenerator: Batch size should be > 0") if ( opengate_core.IsMultithreadedApplication() @@ -50,7 +52,14 @@ def read_phsp_and_keys(self): # FIXME could have an option to select the branch self.root_file = uproot.open(self.user_info.phsp_file) branches = self.root_file.keys() - self.root_file = self.root_file[branches[0]] + if len(branches) > 0: + self.root_file = self.root_file[branches[0]] + else: + gate.fatal( + f"PhaseSpaceSourceGenerator: No useable branches in the root file {self.user_info.phsp_file}. Aborting." + ) + exit() + self.num_entries = int(self.root_file.num_entries) # initialize the index to start @@ -144,24 +153,24 @@ def generate(self, source): f"in the phsp file and no source.particle" ) - # if override_position is set to True, the position + # if translate_position is set to True, the position # supplied will be added to the phsp file position - if ui.override_position: - x = batch[ui.position_key_x] + ui.position.translation[0] - y = batch[ui.position_key_y] + ui.position.translation[1] - z = batch[ui.position_key_z] + ui.position.translation[2] - source.SetPositionXBatch(x) - source.SetPositionYBatch(y) - source.SetPositionZBatch(z) + if ui.translate_position: + batch[ui.position_key_x] += ui.position.translation[0] + batch[ui.position_key_y] += ui.position.translation[1] + batch[ui.position_key_z] += ui.position.translation[2] + source.SetPositionXBatch(batch[ui.position_key_x]) + source.SetPositionYBatch(batch[ui.position_key_y]) + source.SetPositionZBatch(batch[ui.position_key_z]) else: source.SetPositionXBatch(batch[ui.position_key_x]) source.SetPositionYBatch(batch[ui.position_key_y]) source.SetPositionZBatch(batch[ui.position_key_z]) # direction is a rotation of the stored direction - # if override_direction is set to True, the direction + # if rotate_direction is set to True, the direction # in the root file will be rotated based on the supplied rotation matrix - if ui.override_direction: + if ui.rotate_direction: # create point vectors self.points = np.column_stack( ( @@ -224,7 +233,9 @@ def set_default_user_info(user_info): SourceBase.set_default_user_info(user_info) # initial user info user_info.phsp_file = None - user_info.n = 1 + user_info.n = 0 + user_info.activity = 0 + user_info.half_life = -1 # negative value is not half_life user_info.particle = "" # FIXME later as key user_info.entry_start = 0 # if a particle name is supplied, the particle type is set to it @@ -251,11 +262,14 @@ def set_default_user_info(user_info): # change position and direction of the source # position is relative to the stored coordinates # direction is a rotation of the stored direction - user_info.override_position = False - user_info.override_direction = False + user_info.translate_position = False + user_info.rotate_direction = False user_info.position = Box() user_info.position.translation = [0, 0, 0] user_info.position.rotation = Rotation.identity().as_matrix() + user_info.generate_until_next_primary = False + user_info.primary_lower_energy_threshold = 0 + user_info.primary_PDGCode = 0 # user_info.time_key = None # FIXME TODO later # for debug user_info.verbose_batch = False @@ -287,6 +301,18 @@ def initialize(self, run_timing_intervals): if ui.direction_key_z is None: ui.direction_key_z = f"{ui.direction_key}_Z" + # check if the source should generate particles until the second one + # which is identified as primary by name, PDGCode and above a threshold + if ui.generate_until_next_primary == True: + if ui.primary_PDGCode == 0: + gate.fatal( + f"PhaseSpaceSource: generate_until_next_primary is True but no primary particle is defined" + ) + if ui.primary_lower_energy_threshold <= 0: + gate.fatal( + f"PhaseSpaceSource: generate_until_next_primary is True but no primary_lower_energy_threshold is defined" + ) + # initialize the generator (read the phsp file) self.particle_generator.initialize(self.user_info) diff --git a/opengate/tests/data b/opengate/tests/data index 8cf78354d..a0e299e40 160000 --- a/opengate/tests/data +++ b/opengate/tests/data @@ -1 +1 @@ -Subproject commit 8cf78354d4655eb630f223e5097e32c3f7f46b73 +Subproject commit a0e299e40bb878b7bdeed76d0f9f1a6bfa993557 diff --git a/opengate/tests/src/test060_phsp_source_helpers.py b/opengate/tests/src/test060_phsp_source_helpers.py index ca2db9fb0..98a9dd651 100644 --- a/opengate/tests/src/test060_phsp_source_helpers.py +++ b/opengate/tests/src/test060_phsp_source_helpers.py @@ -332,7 +332,7 @@ def test_source_translation( source.particle = "proton" source.batch_size = 3000 source.n = number_of_particles - source.override_position = True + source.translate_position = True source.position.translation = [3 * cm, 0 * cm, 0 * cm] print(source) @@ -360,9 +360,9 @@ def test_source_rotation( source.particle = "proton" source.batch_size = 3000 source.n = number_of_particles - # source.override_position = True + # source.translate_position = True # source.position.translation = [3 * cm, 1 * cm, 0 * cm] - source.override_direction = True + source.rotate_direction = True # rotation = Rotation.from_euler("zyx", [30, 20, 10], degrees=True) rotation = Rotation.from_euler("x", [30], degrees=True) source.position.rotation = rotation.as_matrix() @@ -371,6 +371,36 @@ def test_source_rotation( sim.run() +def test_source_untilPrimary( + source_file_name="output/test_proton_offset.root", + phs_file_name_out="output/output/test_source_electron.root", +) -> None: + sim = create_phs_without_source( + phs_name=phs_file_name_out, + ) + number_of_particles = 2 + ########################################################################################## + # Source + ########################################################################################## + # phsp source + source = sim.add_source("PhaseSpaceSource", "phsp_source_global") + source.mother = "world" + source.phsp_file = source_file_name + source.position_key = "PrePosition" + source.direction_key = "PreDirection" + source.global_flag = True + source.particle = "" + source.batch_size = 3000 + source.n = number_of_particles + source.generate_until_next_primary = True + source.primary_lower_energy_threshold = 90.0 * MeV + source.primary_PDGCode = 2212 + print(source) + + sim.run() + output = sim.output + + def get_first_entry_of_key( file_name_root="output/test_source_electron.root", key="ParticleName" ) -> None: diff --git a/opengate/tests/src/test060_phsp_source_particlename_direct.py b/opengate/tests/src/test060_phsp_source_particletype_direct.py similarity index 100% rename from opengate/tests/src/test060_phsp_source_particlename_direct.py rename to opengate/tests/src/test060_phsp_source_particletype_direct.py diff --git a/opengate/tests/src/test060_phsp_source_particlename_fromphs_particlename_wip.py b/opengate/tests/src/test060_phsp_source_particletype_fromphs_particlename_wip.py similarity index 100% rename from opengate/tests/src/test060_phsp_source_particlename_fromphs_particlename_wip.py rename to opengate/tests/src/test060_phsp_source_particletype_fromphs_particlename_wip.py diff --git a/opengate/tests/src/test060_phsp_source_particlename_fromphs_pdgcode.py b/opengate/tests/src/test060_phsp_source_particletype_fromphs_pdgcode.py similarity index 100% rename from opengate/tests/src/test060_phsp_source_particlename_fromphs_pdgcode.py rename to opengate/tests/src/test060_phsp_source_particletype_fromphs_pdgcode.py diff --git a/opengate/tests/src/test060_phsp_source_untilNextPrimary.py b/opengate/tests/src/test060_phsp_source_untilNextPrimary.py new file mode 100644 index 000000000..c90d014c0 --- /dev/null +++ b/opengate/tests/src/test060_phsp_source_untilNextPrimary.py @@ -0,0 +1,95 @@ +import test060_phsp_source_helpers as t +import opengate as gate +import uproot +import numpy as np +import gatetools.phsp as phsp +from opengate.tests import utility + + +paths = utility.get_default_test_paths( + __file__, "test060_PhsSource_ParticleName_direct", output_folder="test060" +) + +# units +m = gate.g4_units.m +mm = gate.g4_units.mm +cm = gate.g4_units.cm +nm = gate.g4_units.nm +Bq = gate.g4_units.Bq +MeV = gate.g4_units.MeV +deg: float = gate.g4_units.deg + + +def createRootFile(file_name): + # create numpy arrays + kinetic_energy = np.array([100, 2.0, 3.0, 80, 2.0, 3.0, 2.0, 110.0, 3.0, 111.0, 2]) + pdg_code = np.array([2212, 11, 11, 2212, 11, 11, 11, 2212, 11, 2212, 11]) + pre_direction_local_x = np.zeros(11) + pre_direction_local_y = np.zeros(11) + pre_direction_local_z = np.ones(11) + pre_position_local_x = np.zeros(11) + pre_position_local_y = np.zeros(11) + pre_position_local_z = np.zeros(11) + pre_direction_x = np.zeros(11) + pre_direction_y = np.zeros(11) + pre_direction_z = np.ones(11) + pre_position_x = np.zeros(11) + pre_position_y = np.zeros(11) + pre_position_z = np.zeros(11) + weight = np.ones(11) + + # generate root file + tfile = uproot.recreate(file_name) + tfile["PhaseSpace1"] = { + "KineticEnergy": kinetic_energy, + "PDGCode": pdg_code, + "PreDirectionLocal_X": pre_direction_local_x, + "PreDirectionLocal_Y": pre_direction_local_y, + "PreDirectionLocal_Z": pre_direction_local_z, + "PrePositionLocal_X": pre_position_local_x, + "PrePositionLocal_Y": pre_position_local_y, + "PrePositionLocal_Z": pre_position_local_z, + "PreDirection_X": pre_direction_x, + "PreDirection_Y": pre_direction_y, + "PreDirection_Z": pre_direction_z, + "PrePosition_X": pre_position_x, + "PrePosition_Y": pre_position_y, + "PrePosition_Z": pre_position_z, + "Weight": weight, + } + tfile.close() + + +def main(): + print("create reference PhS file") + createRootFile(paths.output / "test_phs.root") + + print("testing until primary") + t.test_source_untilPrimary( + source_file_name=paths.output / "test_phs.root", + phs_file_name_out=paths.output / "test_source_untilPrimary.root", + ) + + # load data from root file + data_ref, keys_ref, m_ref = phsp.load( + paths.output / "test_source_untilPrimary.root" + ) + # the root file should contain 9 particles + # print(m_ref) + + # there should be 3 protons in the root file + index = keys_ref.index("PDGCode") + particle_list = data_ref[:, index] + # print(particle_list) + count = np.count_nonzero(particle_list == 2212) + is_ok = False + if m_ref == 9 and count == 3: + is_ok = True + print("test ok") + + # this is the end, my friend + utility.test_ok(is_ok) + + +if __name__ == "__main__": + main() diff --git a/opengate/tests/src/test061_TPPhsSource.py b/opengate/tests/src/test061_TPPhsSource.py new file mode 100644 index 000000000..dc2d3281a --- /dev/null +++ b/opengate/tests/src/test061_TPPhsSource.py @@ -0,0 +1,79 @@ +import test061_TPPhsSource_helpers as t +import opengate as gate +from scipy.spatial.transform import Rotation +from opengate.tests import utility + +from opengate.contrib.tps.treatmentPlanPhsSource import TreatmentPlanPhsSource +from opengate.contrib.tps.ionbeamtherapy import spots_info_from_txt, TreatmentPlanSource + + +paths = utility.get_default_test_paths(__file__, "test061_TPPhsSource") +paths.output_ref = paths.output_ref / "test061_ref" + +ref_path = paths.output_ref + +# units +m = gate.g4_units.m +mm = gate.g4_units.mm +cm = gate.g4_units.cm +nm = gate.g4_units.nm +Bq = gate.g4_units.Bq +MeV = gate.g4_units.MeV +deg: float = gate.g4_units.deg + + +def main(): + # print("create reference PhS file") + # t.create_test_Phs( + # particle="proton", + # phs_name=paths.output / "test_proton", + # number_of_particles=1, + # translation=[0 * cm, 0 * cm, 0 * mm], + # ) + print("Testing TPPhS source rotations") + + t.test_source_rotation_A( + plan_file_name=ref_path / "PlanSpot.txt", + phs_list_file_name="PhsList.txt", + phs_folder_name=ref_path, + phs_file_name_out=paths.output / "output.root", + ) + + a = t.check_value_from_root_file( + file_name_root=paths.output / "output.root", + key="KineticEnergy", + ref_value=150 * MeV, + ) + b = t.check_value_from_root_file( + file_name_root=paths.output / "output.root", + key="PrePositionLocal_X", + ref_value=-60, + ) + c = t.check_value_from_root_file( + file_name_root=paths.output / "output.root", + key="PrePositionLocal_Y", + ref_value=50, + ) + d = t.check_value_from_root_file( + file_name_root=paths.output / "output.root", + key="PrePositionLocal_Z", + ref_value=0, + ) + + e = t.check_value_from_root_file( + file_name_root=paths.output / "output.root", + key="PreDirectionLocal_X", + ref_value=-0.012, + ) + f = t.check_value_from_root_file( + file_name_root=paths.output / "output.root", + key="PreDirectionLocal_Y", + ref_value=0.01, + ) + + # this is the end, my friend + utility.test_ok(a and b and c and d and e and f) + + +if __name__ == "__main__": + main() diff --git a/opengate/tests/src/test061_TPPhsSource_helpers.py b/opengate/tests/src/test061_TPPhsSource_helpers.py new file mode 100644 index 000000000..d99e6d8d4 --- /dev/null +++ b/opengate/tests/src/test061_TPPhsSource_helpers.py @@ -0,0 +1,266 @@ +import opengate as gate + +# import numpy as np +# import matplotlib.pyplot as plot +from scipy.spatial.transform import Rotation +import gatetools.phsp as phsp +import os +from opengate.tests import utility +from opengate.contrib.tps.treatmentPlanPhsSource import TreatmentPlanPhsSource +from opengate.contrib.tps.ionbeamtherapy import spots_info_from_txt, TreatmentPlanSource + + +# paths = gate.get_default_test_paths( +# __file__, "gate_test019_linac_phsp", output_folder="test019" +# ) +# paths = {output: "output"} + +# units +m = gate.g4_units.m +mm = gate.g4_units.mm +cm = gate.g4_units.cm +nm = gate.g4_units.nm +Bq = gate.g4_units.Bq +MeV = gate.g4_units.MeV +deg: float = gate.g4_units.deg + + +def create_test_Phs( + particle="proton", + phs_name="output/test_proton.root", + number_of_particles=1, + translation=[0 * mm, 0 * mm, 0 * mm], +): + # create the simulation + sim = gate.Simulation() + + # main options + ui = sim.user_info + ui.g4_verbose = False + # ui.visu = True + ui.visu_type = "vrml" + ui.check_volumes_overlap = False + # ui.running_verbose_level = gate.EVENT + ui.number_of_threads = 1 + ui.random_seed = "auto" + + ########################################################################################## + # geometry + ########################################################################################## + # adapt world size + world = sim.world + world.size = [1 * m, 1 * m, 1 * m] + world.material = "G4_Galactic" + + # virtual plane for phase space + plane = sim.add_volume("Tubs", "phase_space_plane") + # plane.material = "G4_AIR" + plane.material = "G4_Galactic" + plane.rmin = 0 + plane.rmax = 30 * cm + plane.dz = 1 * nm # half height + # plane.rotation = Rotation.from_euler("xy", [180, 30], degrees=True).as_matrix() + # plane.translation = [0 * mm, 0 * mm, 0 * mm] + plane.translation = translation + + plane.color = [1, 0, 0, 1] # red + + ########################################################################################## + # Actors + ########################################################################################## + # Split the joined path into the directory path and filename + directory_path, filename = os.path.split(phs_name) + # Extract the base filename and extension + base_filename, extension = os.path.splitext(filename) + new_extension = ".root" + + # PhaseSpace Actor + ta1 = sim.add_actor("PhaseSpaceActor", "PhaseSpace1") + ta1.mother = plane.name + ta1.attributes = [ + "KineticEnergy", + "Weight", + "PrePosition", + "PrePositionLocal", + "ParticleName", + "PreDirection", + "PreDirectionLocal", + "PDGCode", + ] + new_joined_path = os.path.join(directory_path, base_filename + new_extension) + ta1.output = new_joined_path + ta1.debug = False + f = sim.add_filter("ParticleFilter", "f") + f.particle = particle + ta1.filters.append(f) + + sim.physics_manager.physics_list_name = "QGSP_BIC_EMZ" + + ########################################################################################## + # Source + ########################################################################################## + source = sim.add_source("GenericSource", "particle_source") + source.mother = "world" + # source.particle = "ion 6 12" # Carbon ions + source.particle = "proton" + source.position.type = "point" + source.direction.type = "momentum" + source.direction.momentum = [0, 0, 1] + source.position.translation = [0 * cm, 0 * cm, -0.1 * cm] + source.energy.type = "mono" + source.energy.mono = 150 * MeV + source.n = number_of_particles + + sim.run() + output = sim.output + + +def create_PhS_withoutSource( + phs_name="output/test_proton.root", +): + # create the simulation + sim = gate.Simulation() + + # main options + ui = sim.user_info + ui.g4_verbose = False + # ui.visu = True + ui.visu_type = "vrml" + ui.check_volumes_overlap = False + # ui.running_verbose_level = gate.EVENT + ui.number_of_threads = 1 + ui.random_seed = "auto" + + ########################################################################################## + # geometry + ########################################################################################## + # adapt world size + world = sim.world + world.size = [1 * m, 1 * m, 1 * m] + world.material = "G4_Galactic" + print(world.name) + + # virtual plane for phase space + plane = sim.add_volume("Tubs", "phase_space_plane") + # plane.material = "G4_AIR" + plane.material = "G4_Galactic" + plane.rmin = 0 + plane.rmax = 30 * cm + plane.dz = 1 * nm # half height + # plane.rotation = Rotation.from_euler("xy", [180, 30], degrees=True).as_matrix() + plane.translation = [0 * mm, 0 * mm, 0 * mm] + # plane.translation = translation + + plane.color = [1, 0, 0, 1] # red + + ########################################################################################## + # Actors + ########################################################################################## + + # PhaseSpace Actor + ta1 = sim.add_actor("PhaseSpaceActor", "PhaseSpace") + ta1.mother = plane.name + ta1.attributes = [ + "KineticEnergy", + "Weight", + "PrePosition", + "PrePositionLocal", + "ParticleName", + "PreDirection", + "PreDirectionLocal", + "PDGCode", + ] + ta1.output = phs_name + ta1.debug = False + + sim.physics_manager.physics_list_name = "QGSP_BIC_EMZ" + + # ########################################################################################## + # # Source + # ########################################################################################## + # # phsp source + # source = sim.add_source("PhaseSpaceSource", "phsp_source_global") + # source.mother = world.name + # source.phsp_file = source_name + # source.position_key = "PrePosition" + # source.direction_key = "PreDirection" + # source.global_flag = True + # source.particle = particle + # source.batch_size = 3000 + # source.n = number_of_particles / ui.number_of_threads + # # source.position.translation = [0 * cm, 0 * cm, -35 * cm] + + # output = sim.run() + # output = sim.start() + # output = sim.start(start_new_process=True) + return sim, plane + + +def test_source_rotation_A( + plan_file_name="output/test_proton_offset.root", + phs_list_file_name="PhsList.txt", + phs_folder_name="", + phs_file_name_out="output/output/test_source_electron.root", +) -> None: + sim, plane = create_PhS_withoutSource( + phs_name=phs_file_name_out, + ) + number_of_particles = 1 + ########################################################################################## + # Source + ########################################################################################## + # TreatmentPlanPhsSource source + tpPhSs = TreatmentPlanPhsSource("RT_plan", sim) + tpPhSs.set_phaseSpaceList_file_name(phs_list_file_name) + tpPhSs.set_phaseSpaceFolder(phs_folder_name) + spots, ntot, energies, G = spots_info_from_txt(plan_file_name, "") + # print(spots, ntot, energies, G) + tpPhSs.set_spots(spots) + tpPhSs.set_particles_to_simulate(number_of_particles) + tpPhSs.set_distance_source_to_isocenter(100 * cm) + tpPhSs.set_distance_stearmag_to_isocenter(5 * m, 5 * m) + tpPhSs.rotation = Rotation.from_euler("z", G, degrees=True) + tpPhSs.initialize_tpPhssource() + + # depending on the rotation of the gantry, the rotation of the phase space to catch the particles is different + plane.rotation = Rotation.from_euler("y", 90, degrees=True).as_matrix() + + sim.run() + output = sim.output + + +def get_first_entry_of_key( + file_name_root="output/test_source_electron.root", key="ParticleName" +) -> None: + # read root file + data_ref, keys_ref, m_ref = phsp.load(file_name_root) + # print(data_ref) + # print(keys_ref) + index = keys_ref.index(key) + # print(index) + # print(data_ref[index][0]) + return data_ref[0][index] + + +def check_value_from_root_file( + file_name_root="output/test_source_electron.root", + key="ParticleName", + ref_value="e-", +): + # read root file + value = get_first_entry_of_key(file_name_root=file_name_root, key=key) + if (type(ref_value) != str) and (type(value) != str): + is_ok = utility.check_diff_abs( + float(value), float(ref_value), tolerance=1e-3, txt=key + ) + # utility.check_diff_abs(float(value), float(ref_value), tolerance=1e-6, txt=key) + else: + if value == ref_value: + # print("Is correct") + is_ok = True + else: + # print("Is not correct") + is_ok = False + # print("ref_value: ", ref_value) + # print("value: ", value) + return is_ok diff --git a/opengate/tests/src/test066_phspsource_activity.py b/opengate/tests/src/test066_phspsource_activity.py new file mode 100644 index 000000000..c25495cd9 --- /dev/null +++ b/opengate/tests/src/test066_phspsource_activity.py @@ -0,0 +1,119 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import opengate as gate +import uproot +import numpy as np +from opengate.tests import utility + + +def validation_test(n, n_measured): + n_part_minus = n - 3 * np.sqrt(n) + n_part_max = n + 3 * np.sqrt(n) + if n_part_minus <= n_measured <= n_part_max: + return True + else: + return False + + +if __name__ == "__main__": + paths = utility.get_default_test_paths( + __file__, "test066_PhsSource_Activity", output_folder="test066" + ) + + # create the simulation + sim = gate.Simulation() + + # main options + ui = sim.user_info + ui.g4_verbose = False + # ui.visu = True + ui.visu_type = "vrml" + ui.check_volumes_overlap = False + # ui.running_verbose_level = gate.EVENT + ui.number_of_threads = 1 + ui.random_seed = "auto" + + # units + m = gate.g4_units.m + km = gate.g4_units.km + mm = gate.g4_units.mm + cm = gate.g4_units.cm + nm = gate.g4_units.nm + Bq = gate.g4_units.Bq + MeV = gate.g4_units.MeV + keV = gate.g4_units.keV + gcm3 = gate.g4_units.g / gate.g4_units.cm3 + + # adapt world size + world = sim.world + world.size = [1 * m, 1 * m, 2 * m] + world.material = "G4_Galactic" + + # source_1 + nb_part = 1000 + plane_1 = sim.add_volume("Box", "plan_1") + plane_1.material = "G4_Galactic" + plane_1.mother = world.name + plane_1.size = [1 * m, 1 * m, 1 * nm] + plane_1.color = [0.8, 0.2, 0.1, 1] + + source_1 = sim.add_source("PhaseSpaceSource", "phsp_source_global_1") + source_1.mother = plane_1.name + source_1.phsp_file = paths.output_ref / ".." / "test019" / "test019_hits.root" + source_1.position_key = "PrePosition" + source_1.direction_key = "PreDirection" + source_1.weight_key = "Weight" + source_1.global_flag = False + source_1.particle = "gamma" + source_1.batch_size = 100 + source_1.translate_position = True + source_1.position.translation = [0, 0, 300] + source_1.activity = nb_part * Bq + + # add phase space plan 1 + + phsp_plane_1 = sim.add_volume("Box", "phase_space_plane_1") + phsp_plane_1.mother = world.name + phsp_plane_1.material = "G4_Galactic" + phsp_plane_1.size = [1 * m, 1 * m, 1 * nm] + phsp_plane_1.translation = [0, 0, -1000 * nm] + phsp_plane_1.color = [1, 0, 0, 1] # red + + phsp_actor = sim.add_actor("PhaseSpaceActor", "PhaseSpace") + phsp_actor.mother = phsp_plane_1.name + phsp_actor.attributes = [ + "KineticEnergy", + ] + + phsp_actor.output = paths.output / "test066_output_data.root" + + s = sim.add_actor("SimulationStatisticsActor", "Stats") + s.track_types_flag = True + sim.physics_manager.physics_list_name = "G4EmStandardPhysics_option3" + sim.physics_manager.enable_decay = False + sim.physics_manager.global_production_cuts.gamma = 1 * mm + sim.physics_manager.global_production_cuts.electron = 1 * mm + sim.physics_manager.global_production_cuts.positron = 1 * mm + # + # # go ! + sim.run() + output = sim.output + + # + # # print results + stats = sim.output.get_actor("Stats") + h = output.get_actor("PhaseSpace") + print(stats) + # + f_phsp = uproot.open(paths.output / "test066_output_data.root") + arr = f_phsp["PhaseSpace"].arrays() + print("Number of detected events :", len(arr)) + print( + "Number of expected events :", + nb_part, + "+- 3*" + str(int(np.round(np.sqrt(nb_part)))), + ) + + is_ok = validation_test(nb_part, len(arr)) + utility.test_ok(is_ok) diff --git a/opengate/tests/src/test067_cbct_fluence_actor_mt.py b/opengate/tests/src/test067_cbct_fluence_actor_mt.py index 58f903a9e..d2161e999 100755 --- a/opengate/tests/src/test067_cbct_fluence_actor_mt.py +++ b/opengate/tests/src/test067_cbct_fluence_actor_mt.py @@ -85,9 +85,10 @@ # run sim.run() + output = sim.output # print output statistics - stats = sim.output.get_actor("stats") + stats = output.get_actor("stats") print(stats) # check images diff --git a/opengate/tests/src/test068_rotation_source.py b/opengate/tests/src/test068_rotation_source.py new file mode 100644 index 000000000..7c2fe7761 --- /dev/null +++ b/opengate/tests/src/test068_rotation_source.py @@ -0,0 +1,176 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import opengate as gate +from opengate.tests import utility +from scipy.spatial.transform import Rotation +import numpy as np +import uproot + + +# The test generates two differents generic sources, momentum and focused, which is attached to a plan. +# The plan is randomly rotated and we verify that the generated particles have a direction which is in accordance +# with the applied rotations and the transmissions. + + +def test068(tab, nb_run, nb_part, nb_source): + nb_event = len(tab) + nb_event_theo = nb_run * nb_part * nb_source + err = np.sqrt(nb_event_theo) + if nb_event_theo - 3 * err < nb_event < nb_event_theo + 3 * err: + return True + else: + return False + + +if __name__ == "__main__": + paths = utility.get_default_test_paths(__file__) + output_path = paths.output + + # create the simulation + sim = gate.Simulation() + + # main options + ui = sim.user_info + ui.g4_verbose = False + # ui.visu = True + ui.visu_type = "vrml" + ui.check_volumes_overlap = False + # ui.running_verbose_level = gate.logger.EVENT + ui.number_of_threads = 1 + ui.random_seed = "auto" + + # units + m = gate.g4_units.m + km = gate.g4_units.km + mm = gate.g4_units.mm + cm = gate.g4_units.cm + nm = gate.g4_units.nm + Bq = gate.g4_units.Bq + MeV = gate.g4_units.MeV + keV = gate.g4_units.keV + sec = gate.g4_units.s + gcm3 = gate.g4_units["g/cm3"] + + # adapt world size + world = sim.world + world.size = [1 * m, 1 * m, 1 * m] + world.material = "G4_Galactic" + + # nb of particles + n = 100 + + # Plan to attach the source + nb_source = 2 + plan = sim.add_volume("Box", "source_plan") + plan.mother = world.name + plan.material = "G4_Galactic" + plan.size = [10 * cm, 10 * cm, 1 * nm] + plan.translation = np.array([0 * cm, 0 * cm, 0 * cm]) + + source1 = sim.add_source("GenericSource", "photon_source") + source1.particle = "gamma" + source1.position.type = "box" + source1.mother = plan.name + source1.position.size = [10 * cm, 10 * cm, 1 * nm] + source1.direction.type = "momentum" + source1.force_rotation = True + # source1.direction.focus_point = [0*cm, 0*cm, -5 *cm] + source1.direction.momentum = [0, 0, -1] + source1.energy.type = "mono" + source1.energy.mono = 1 * MeV + source1.activity = n * Bq / ui.number_of_threads + + source2 = sim.add_source("GenericSource", "photon_source_2") + source2.particle = "gamma" + source2.position.type = "box" + source2.mother = plan.name + source2.position.size = [10 * cm, 10 * cm, 1 * nm] + source2.direction.type = "focused" + source2.force_rotation = True + source2.direction.focus_point = [0 * cm, 0 * cm, -5 * cm] + # source1.direction.momentum = [0, 0, -1] + source2.energy.type = "mono" + source2.energy.mono = 1 * MeV + source2.activity = n * Bq / ui.number_of_threads + + # Phase Space + + phsp_plan = sim.add_volume("Box", "phsp_plan") + phsp_plan.mother = world.name + phsp_plan.material = "G4_Galactic" + phsp_plan.size = [10 * cm, 10 * cm, 1 * nm] + phsp_plan.translation = np.array([0 * cm, 0 * cm, -5 * cm]) + + phsp = sim.add_actor("PhaseSpaceActor", "PhaseSpace") + phsp.mother = phsp_plan.name + phsp.attributes = ["EventID"] + phsp.output = output_path / "test068.root" + + # MOTION ACTOR for the source and actors + motion_source = sim.add_actor("MotionVolumeActor", "Move_source") + motion_source.mother = plan.name + motion_source.rotations = [] + motion_source.translations = [] + + motion_phsp = sim.add_actor("MotionVolumeActor", "Move_phsp") + motion_phsp.mother = phsp_plan.name + motion_phsp.rotations = [] + motion_phsp.translations = [] + sim.run_timing_intervals = [] + + for i in range(10): + x_translation = 20 * cm * np.random.random() + y_translation = 20 * cm * np.random.random() + z_translation = 20 * cm * np.random.random() + + x_rot = 360 * np.random.random() + y_rot = 360 * np.random.random() + z_rot = 360 * np.random.random() + + motion_source.translations.append([x_translation, y_translation, z_translation]) + rot = Rotation.from_euler( + "xyz", [x_rot, y_rot, z_rot], degrees=True + ).as_matrix() + vec_source_phsp = phsp_plan.translation - plan.translation + # vec_source_phsp = vec_source_phsp.reshape((len(vec_source_phsp),1)) + # print(vec_source_phsp) + translation_phsp = np.dot(rot, vec_source_phsp) + np.array( + [x_translation, y_translation, z_translation] + ) + motion_source.rotations.append(rot) + motion_phsp.rotations.append(rot) + motion_phsp.translations.append(translation_phsp) + sim.run_timing_intervals.append([i * sec, (i + 1) * sec]) + # stat actor + s = sim.add_actor("SimulationStatisticsActor", "Stats") + s.track_types_flag = True + + # Physic list and cuts + sim.physics_manager.physics_list_name = "G4EmStandardPhysics_option3" + sim.physics_manager.enable_decay = False + sim.physics_manager.global_production_cuts.gamma = 1 * mm + sim.physics_manager.global_production_cuts.electron = 1 * mm + sim.physics_manager.global_production_cuts.positron = 1 * mm + + # go ! + sim.run() + output = sim.output + + # Validation test + f_phsp = uproot.open(phsp.output) + arr = f_phsp["PhaseSpace"].arrays() + print("Number of detected events :", len(arr)) + nb_event = len(arr) + nb_part = n + nb_run = len(sim.run_timing_intervals) + err = np.sqrt(nb_part * nb_run * nb_source) + print( + "Number of expected events : [" + + str(int(-3 * err + nb_run * nb_part * nb_source)) + + "," + + str(int(+3 * err + nb_run * nb_part * nb_source)) + + "]" + ) + is_ok = test068(arr, len(sim.run_timing_intervals), nb_part, 2) + utility.test_ok(is_ok)