diff --git a/core/opengate_core/opengate_core.cpp b/core/opengate_core/opengate_core.cpp index e9ab29218..07e526ce7 100644 --- a/core/opengate_core/opengate_core.cpp +++ b/core/opengate_core/opengate_core.cpp @@ -320,6 +320,10 @@ void init_GateKillActor(py::module &); void init_GateKillAccordingProcessesActor(py::module &); +void init_GateKillNonInteractingParticleActor(py::module &); + +void init_GateKillAccordingParticleNameActor(py::module &); + void init_GateAttenuationImageActor(py::module &); void init_itk_image(py::module &); @@ -340,6 +344,8 @@ void init_GatePhaseSpaceActor(py::module &); void init_GateOptrComptSplittingActor(py::module &m); +void init_GateLastVertexInteractionSplittingActor(py::module &m); + void init_GateBOptrBremSplittingActor(py::module &m); void init_GateOptrFreeFlightActor(py::module &m); @@ -378,6 +384,8 @@ void init_GateVolumeDepthID(py::module &m); // Gate source void init_GateVSource(py::module &); +void init_GateLastVertexSource(py::module &); + void init_GateSourceManager(py::module &); void init_GateGenericSource(py::module &); @@ -569,6 +577,7 @@ PYBIND11_MODULE(opengate_core, m) { init_GateImageNestedParameterisation(m); init_GateRepeatParameterisation(m); init_GateVSource(m); + init_GateLastVertexSource(m); init_GateSourceManager(m); init_GateGenericSource(m); init_GateTreatmentPlanPBSource(m); @@ -592,6 +601,7 @@ PYBIND11_MODULE(opengate_core, m) { init_GatePhaseSpaceActor(m); init_GateBOptrBremSplittingActor(m); init_GateOptrComptSplittingActor(m); + init_GateLastVertexInteractionSplittingActor(m); init_GateOptrFreeFlightActor(m); init_GateHitsCollectionActor(m); init_GateVDigitizerWithOutputActor(m); @@ -606,6 +616,8 @@ PYBIND11_MODULE(opengate_core, m) { init_GateARFTrainingDatasetActor(m); init_GateKillActor(m); init_GateKillAccordingProcessesActor(m); + init_GateKillNonInteractingParticleActor(m); + init_GateKillAccordingParticleNameActor(m); init_GateAttenuationImageActor(m); init_GateDigiAttributeManager(m); init_GateVDigiAttribute(m); diff --git a/core/opengate_core/opengate_lib/GateGenericSource.cpp b/core/opengate_core/opengate_lib/GateGenericSource.cpp index 1f30c32ea..3c23ac7d4 100644 --- a/core/opengate_core/opengate_lib/GateGenericSource.cpp +++ b/core/opengate_core/opengate_lib/GateGenericSource.cpp @@ -145,7 +145,6 @@ double GateGenericSource::PrepareNextTime(double current_simulation_time) { return fStartTime; if (ll.fEffectiveEventTime >= fEndTime) return -1; - // get next time according to current fActivity double next_time = CalcNextTime(ll.fEffectiveEventTime); if (next_time >= fEndTime) diff --git a/core/opengate_core/opengate_lib/GateKillAccordingParticleNameActor.cpp b/core/opengate_core/opengate_lib/GateKillAccordingParticleNameActor.cpp new file mode 100644 index 000000000..f4c295245 --- /dev/null +++ b/core/opengate_core/opengate_lib/GateKillAccordingParticleNameActor.cpp @@ -0,0 +1,57 @@ +/* -------------------------------------------------- + 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 "GateKillAccordingParticleNameActor.h" +#include "G4LogicalVolumeStore.hh" +#include "G4PhysicalVolumeStore.hh" +#include "G4TransportationManager.hh" +#include "G4ios.hh" +#include "GateHelpers.h" +#include "GateHelpersDict.h" + +G4Mutex SetNbKillAcordingParticleMutex = G4MUTEX_INITIALIZER; + +GateKillAccordingParticleNameActor::GateKillAccordingParticleNameActor( + py::dict &user_info) + : GateVActor(user_info, false) {} + +void GateKillAccordingParticleNameActor::InitializeUserInfo( + py::dict &user_info) { + GateVActor::InitializeUserInfo(user_info); + fParticlesNameToKill = DictGetVecStr(user_info, "particles_name_to_kill"); +} + +void GateKillAccordingParticleNameActor::PreUserTrackingAction( + const G4Track *track) { + auto &l = fThreadLocalData.Get(); + l.fIsAParticleToKill = false; + + G4String particleName = track->GetParticleDefinition()->GetParticleName(); + if (std::find(fParticlesNameToKill.begin(), fParticlesNameToKill.end(), + particleName) != fParticlesNameToKill.end()) { + l.fIsAParticleToKill = true; + } +} + +void GateKillAccordingParticleNameActor::SteppingAction(G4Step *step) { + + if (step->GetPostStepPoint()->GetStepStatus() == 1) { + G4String logicalVolumeNamePostStep = step->GetPostStepPoint() + ->GetPhysicalVolume() + ->GetLogicalVolume() + ->GetName(); + if (std::find(fListOfVolumeAncestor.begin(), fListOfVolumeAncestor.end(), + logicalVolumeNamePostStep) != fListOfVolumeAncestor.end()) { + auto &l = fThreadLocalData.Get(); + if (l.fIsAParticleToKill) { + step->GetTrack()->SetTrackStatus(fStopAndKill); + G4AutoLock mutex(&SetNbKillAcordingParticleMutex); + fNbOfKilledParticles++; + } + } + } +} \ No newline at end of file diff --git a/core/opengate_core/opengate_lib/GateKillAccordingParticleNameActor.h b/core/opengate_core/opengate_lib/GateKillAccordingParticleNameActor.h new file mode 100644 index 000000000..4da5a246a --- /dev/null +++ b/core/opengate_core/opengate_lib/GateKillAccordingParticleNameActor.h @@ -0,0 +1,40 @@ +/* -------------------------------------------------- + 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 GateKillAccordingParticleNameActor_h +#define GateKillAccordingParticleNameActor_h + +#include "G4Cache.hh" +#include "GateVActor.h" +#include + +namespace py = pybind11; + +class GateKillAccordingParticleNameActor : public GateVActor { + +public: + // Constructor + GateKillAccordingParticleNameActor(py::dict &user_info); + struct threadLocalT { + G4bool fIsAParticleToKill = false; + }; + G4Cache fThreadLocalData; + std::vector fParticlesNameToKill; + std::vector fListOfVolumeAncestor; + + // Main function called every step in attached volume + void PreUserTrackingAction(const G4Track *) override; + void SteppingAction(G4Step *) override; + void InitializeUserInfo(py::dict &user_info) override; + + inline long GetNumberOfKilledParticles() { return fNbOfKilledParticles; } + +private: + long fNbOfKilledParticles = 0; +}; + +#endif diff --git a/core/opengate_core/opengate_lib/GateKillNonInteractingParticleActor.cpp b/core/opengate_core/opengate_lib/GateKillNonInteractingParticleActor.cpp new file mode 100644 index 000000000..0a0a1ea22 --- /dev/null +++ b/core/opengate_core/opengate_lib/GateKillNonInteractingParticleActor.cpp @@ -0,0 +1,81 @@ +/* -------------------------------------------------- + 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 "GateKillNonInteractingParticleActor.h" +#include "G4LogicalVolumeStore.hh" +#include "G4PhysicalVolumeStore.hh" +#include "G4TransportationManager.hh" +#include "G4ios.hh" +#include "GateHelpers.h" +#include "GateHelpersDict.h" + +GateKillNonInteractingParticleActor::GateKillNonInteractingParticleActor( + py::dict &user_info) + : GateVActor(user_info, false) {} + +void GateKillNonInteractingParticleActor::StartSimulationAction() { + fNbOfKilledParticles = 0; +} + +void GateKillNonInteractingParticleActor::PreUserTrackingAction( + const G4Track *track) { + fIsFirstStep = true; + fKineticEnergyAtTheEntrance = 0; + ftrackIDAtTheEntrance = 0; + fPassedByTheMotherVolume = false; +} + +void GateKillNonInteractingParticleActor::SteppingAction(G4Step *step) { + G4Navigator *navigator = G4TransportationManager::GetTransportationManager() + ->GetNavigatorForTracking(); + + G4String logNameMotherVolume = G4LogicalVolumeStore::GetInstance() + ->GetVolume(fAttachedToVolumeName) + ->GetName(); + G4String physicalVolumeNamePreStep = "None"; + if (step->GetPreStepPoint()->GetPhysicalVolume() != 0) + physicalVolumeNamePreStep = + step->GetPreStepPoint()->GetPhysicalVolume()->GetName(); + // if (((step->GetTrack()->GetLogicalVolumeAtVertex()->GetName() != + // logNameMotherVolume) && (fIsFirstStep)) || ((fIsFirstStep) && + // (step->GetTrack()->GetParentID() == 0))) { + if ((step->GetTrack()->GetLogicalVolumeAtVertex()->GetName() != + logNameMotherVolume) && + (fIsFirstStep)) { + // if ((fPassedByTheMotherVolume == false) && + // (((step->GetPreStepPoint()->GetStepStatus() == 1) && + // (physicalVolumeNamePreStep == fAttachedToVolumeName)) || ((fIsFirstStep) + // && (step->GetTrack()->GetParentID() == 0)))) { + if ((fPassedByTheMotherVolume == false) && + (((step->GetPreStepPoint()->GetStepStatus() == 1) && + (physicalVolumeNamePreStep == fAttachedToVolumeName)))) { + fPassedByTheMotherVolume = true; + fKineticEnergyAtTheEntrance = step->GetPreStepPoint()->GetKineticEnergy(); + ftrackIDAtTheEntrance = step->GetTrack()->GetTrackID(); + } + } + + G4String logicalVolumeNamePostStep = step->GetPostStepPoint() + ->GetPhysicalVolume() + ->GetLogicalVolume() + ->GetName(); + if ((fPassedByTheMotherVolume) && + (step->GetPostStepPoint()->GetStepStatus() == 1)) { + if (std::find(fListOfVolumeAncestor.begin(), fListOfVolumeAncestor.end(), + logicalVolumeNamePostStep) != fListOfVolumeAncestor.end()) { + if ((step->GetTrack()->GetTrackID() == ftrackIDAtTheEntrance) && + (step->GetPostStepPoint()->GetKineticEnergy() == + fKineticEnergyAtTheEntrance)) { + auto track = step->GetTrack(); + track->SetTrackStatus(fStopAndKill); + fNbOfKilledParticles++; + } + } + } + + fIsFirstStep = false; +} diff --git a/core/opengate_core/opengate_lib/GateKillNonInteractingParticleActor.h b/core/opengate_core/opengate_lib/GateKillNonInteractingParticleActor.h new file mode 100644 index 000000000..b0e0066d1 --- /dev/null +++ b/core/opengate_core/opengate_lib/GateKillNonInteractingParticleActor.h @@ -0,0 +1,40 @@ +/* -------------------------------------------------- + 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 GateKillNonInteractingParticleActor_h +#define GateKillNonInteractingParticleActor_h + +#include "G4Cache.hh" +#include "GateVActor.h" +#include + +namespace py = pybind11; + +class GateKillNonInteractingParticleActor : public GateVActor { + +public: + // Constructor + GateKillNonInteractingParticleActor(py::dict &user_info); + + void StartSimulationAction() override; + + // Main function called every step in attached volume + void SteppingAction(G4Step *) override; + + void PreUserTrackingAction(const G4Track *) override; + + std::vector fParticlesTypeToKill; + G4bool fPassedByTheMotherVolume = false; + G4double fKineticEnergyAtTheEntrance = 0; + G4int ftrackIDAtTheEntrance = 0; + G4bool fIsFirstStep = true; + std::vector fListOfVolumeAncestor; + + long fNbOfKilledParticles{}; +}; + +#endif diff --git a/core/opengate_core/opengate_lib/GateLastVertexInteractionSplittingActor.cpp b/core/opengate_core/opengate_lib/GateLastVertexInteractionSplittingActor.cpp new file mode 100644 index 000000000..26f0a8a3c --- /dev/null +++ b/core/opengate_core/opengate_lib/GateLastVertexInteractionSplittingActor.cpp @@ -0,0 +1,757 @@ +// +// ******************************************************************** +// * License and Disclaimer * +// * * +// * The Geant4 software is copyright of the Copyright Holders of * +// * the Geant4 Collaboration. It is provided under the terms and * +// * conditions of the Geant4 Software License, included in the file * +// * LICENSE and available at http://cern.ch/geant4/license . These * +// * include a list of copyright holders. * +// * * +// * Neither the authors of this software system, nor their employing * +// * institutes,nor the agencies providing financial support for this * +// * work make any representation or warranty, express or implied, * +// * regarding this software system or assume any liability for its * +// * use. Please see the license in the file LICENSE and URL above * +// * for the full disclaimer and the limitation of liability. * +// * * +// * This code implementation is the result of the scientific and * +// * technical work of the GEANT4 collaboration. * +// * By using, copying, modifying or distributing the software (or * +// * any work based on the software) you agree to acknowledge its * +// * use in resulting scientific publications, and indicate your * +// * acceptance of all terms of the Geant4 Software license. * +// ******************************************************************** +// +// +/// \file GateLastVertexInteractionSplittingActor.cc +/// \brief Implementation of the GateLastVertexInteractionSplittingActor class + +#include "GateHelpersDict.h" +#include "GateHelpersImage.h" + +#include "CLHEP/Units/SystemOfUnits.h" +#include "G4BiasingProcessInterface.hh" +#include "G4Gamma.hh" +#include "G4LogicalVolumeStore.hh" +#include "G4ParticleTable.hh" +#include "G4PhysicalVolumeStore.hh" +#include "G4Positron.hh" +#include "G4ProcessManager.hh" +#include "G4ProcessVector.hh" +#include "G4RunManager.hh" +#include "G4TrackStatus.hh" +#include "G4TrackingManager.hh" +#include "G4VParticleChange.hh" +#include "G4eplusAnnihilation.hh" +#include "GateLastVertexInteractionSplittingActor.h" +#include "GateLastVertexSource.h" +#include "GateLastVertexSplittingPostStepDoIt.h" +#include "GateOptnComptSplitting.h" +#include + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +GateLastVertexInteractionSplittingActor:: + GateLastVertexInteractionSplittingActor(py::dict &user_info) + : GateVActor(user_info, false) {} + +void GateLastVertexInteractionSplittingActor::InitializeUserInfo( + py::dict &user_info) { + GateVActor::InitializeUserInfo(user_info); + fAttachedToVolumeName = DictGetStr(user_info, "attached_to"); + fSplittingFactor = DictGetDouble(user_info, "splitting_factor"); + fRotationVectorDirector = DictGetBool(user_info, "rotation_vector_director"); + fAngularKill = DictGetBool(user_info, "angular_kill"); + fVectorDirector = DictGetG4ThreeVector(user_info, "vector_director"); + fMaxTheta = DictGetDouble(user_info, "max_theta"); + fBatchSize = DictGetDouble(user_info, "batch_size"); + fNbOfMaxBatchPerEvent = DictGetInt(user_info, "nb_of_max_batch_per_event"); +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +void GateLastVertexInteractionSplittingActor::print_tree( + const tree &tr, + tree::pre_order_iterator it, + tree::pre_order_iterator end) { + if (!tr.is_valid(it)) + return; + int rootdepth = tr.depth(it); + std::cout << "-----" << std::endl; + while (it != end) { + for (int i = 0; i < tr.depth(it) - rootdepth; ++i) + std::cout << " "; + std::cout << (*it) << std::endl << std::flush; + ++it; + } + std::cout << "-----" << std::endl; +} + +G4bool GateLastVertexInteractionSplittingActor::DoesParticleEmittedInSolidAngle( + G4ThreeVector dir, G4ThreeVector vectorDirector) { + G4double cosTheta = vectorDirector * dir; + if (cosTheta < fCosMaxTheta) + return false; + return true; +} + +G4VProcess *GateLastVertexInteractionSplittingActor::GetProcessFromProcessName( + G4String particleName, G4String pName) { + auto *particle_table = G4ParticleTable::GetParticleTable(); + G4ParticleDefinition *particleDefinition = + particle_table->FindParticle(particleName); + G4ProcessManager *processManager = particleDefinition->GetProcessManager(); + G4ProcessVector *processList = processManager->GetProcessList(); + G4VProcess *nullProcess = nullptr; + for (size_t i = 0; i < processList->size(); ++i) { + auto process = (*processList)[i]; + if (process->GetProcessName() == pName) { + return process; + } + } + return nullProcess; +} + +G4Track *GateLastVertexInteractionSplittingActor::CreateATrackFromContainer( + LastVertexDataContainer theContainer) { + + auto *particle_table = G4ParticleTable::GetParticleTable(); + SimpleContainer container = theContainer.GetContainerToSplit(); + if (container.GetParticleNameToSplit() != "None") { + G4ParticleDefinition *particleDefinition = + particle_table->FindParticle(container.GetParticleNameToSplit()); + G4ThreeVector momentum = container.GetMomentum(); + G4double energy = container.GetEnergy(); + if (energy < 0) { + energy = 0; + momentum = {0, 0, 0}; + } + G4int trackStatus = container.GetTrackStatus(); + G4ThreeVector position = container.GetVertexPosition(); + G4ThreeVector polarization = container.GetPolarization(); + G4DynamicParticle *dynamicParticle = + new G4DynamicParticle(particleDefinition, momentum, energy); + G4double time = 0; + G4Track *aTrack = new G4Track(dynamicParticle, time, position); + aTrack->SetTouchableHandle(container.GetTouchableHandle()); + aTrack->SetPolarization(polarization); + if (trackStatus == 0) { + aTrack->SetTrackStatus(fAlive); + } + if (trackStatus == 1) { + aTrack->SetTrackStatus(fStopButAlive); + } + if ((trackStatus == 2) || (trackStatus == 3)) { + aTrack->SetTrackStatus(fAlive); + } + aTrack->SetWeight(container.GetWeight()); + return aTrack; + } + + return nullptr; +} + +G4Track *GateLastVertexInteractionSplittingActor::CreateComptonTrack( + G4ParticleChangeForGamma *gammaProcess, G4Track track, G4double weight) { + + G4double energy = gammaProcess->GetProposedKineticEnergy(); + G4double globalTime = track.GetGlobalTime(); + G4ThreeVector polarization = gammaProcess->GetProposedPolarization(); + const G4ThreeVector momentum = gammaProcess->GetProposedMomentumDirection(); + const G4ThreeVector position = track.GetPosition(); + G4Track *newTrack = new G4Track(track); + + newTrack->SetKineticEnergy(energy); + newTrack->SetMomentumDirection(momentum); + newTrack->SetPosition(position); + newTrack->SetPolarization(polarization); + newTrack->SetWeight(weight); + return newTrack; +} + +void GateLastVertexInteractionSplittingActor::ComptonSplitting( + G4Step *initStep, G4Step *CurrentStep, G4VProcess *process, + LastVertexDataContainer container, G4double batchSize) { + + // G4TrackVector *trackVector = CurrentStep->GetfSecondary(); + GateGammaEmPostStepDoIt *emProcess = (GateGammaEmPostStepDoIt *)process; + for (int i = 0; i < batchSize; i++) { + G4VParticleChange *processFinalState = + emProcess->PostStepDoIt(*fTrackToSplit, *initStep); + G4ParticleChangeForGamma *gammaProcessFinalState = + (G4ParticleChangeForGamma *)processFinalState; + + G4ThreeVector momentum = + gammaProcessFinalState->GetProposedMomentumDirection(); + + G4Track *newTrack = + CreateComptonTrack(gammaProcessFinalState, *fTrackToSplit, fWeight); + + if ((fAngularKill) && + (DoesParticleEmittedInSolidAngle(newTrack->GetMomentumDirection(), + fVectorDirector) == false)) { + delete newTrack; + } else { + fStackManager->PushOneTrack(newTrack); + } + + // Special case here, since we generate independently each particle, we will + // not attach an electron to exiting compton photon, but we will the + // secondaries. + + if (processFinalState->GetNumberOfSecondaries() > 0) { + delete processFinalState->GetSecondary(0); + } + + processFinalState->Clear(); + gammaProcessFinalState->Clear(); + } +} + +G4Track GateLastVertexInteractionSplittingActor::eBremProcessFinalState( + G4Track *track, G4Step *step, G4VProcess *process) { + // It seem's that the the along step method apply only to brem results to no + // deposited energy but a change in momentum direction according to the + // process Whereas the along step method applied to the ionisation well change + // the deposited energy but not the momentum. Then I apply both to have a + // correct momentum and deposited energy before the brem effect. + G4String particleName = track->GetDefinition()->GetParticleName(); + G4VProcess *eIoniProcess = GetProcessFromProcessName(particleName, "eIoni"); + G4VProcess *eBremProcess = GetProcessFromProcessName(particleName, "eBrem"); + G4VParticleChange *eIoniProcessAlongState = + eIoniProcess->AlongStepDoIt(*track, *step); + G4VParticleChange *eBremProcessAlongState = + eBremProcess->AlongStepDoIt(*track, *step); + G4ParticleChangeForLoss *eIoniProcessAlongStateForLoss = + (G4ParticleChangeForLoss *)eIoniProcessAlongState; + G4ParticleChangeForLoss *eBremProcessAlongStateForLoss = + (G4ParticleChangeForLoss *)eBremProcessAlongState; + G4double LossEnergy = eIoniProcessAlongStateForLoss->GetLocalEnergyDeposit(); + G4ThreeVector momentum = + eBremProcessAlongStateForLoss->GetProposedMomentumDirection(); + G4ThreeVector polarization = + eBremProcessAlongStateForLoss->GetProposedPolarization(); + G4Track aTrack = G4Track(*track); + + aTrack.SetKineticEnergy(track->GetKineticEnergy() - LossEnergy); + aTrack.SetMomentumDirection(momentum); + aTrack.SetPolarization(polarization); + eIoniProcessAlongState->Clear(); + eBremProcessAlongState->Clear(); + return aTrack; +} + +void GateLastVertexInteractionSplittingActor::SecondariesSplitting( + G4Step *initStep, G4Step *CurrentStep, G4VProcess *process, + LastVertexDataContainer theContainer, G4double batchSize) { + SimpleContainer container = theContainer.GetContainerToSplit(); + G4String particleName = + fTrackToSplit->GetParticleDefinition()->GetParticleName(); + // G4TrackVector *trackVector = CurrentStep->GetfSecondary(); + + G4VParticleChange *processFinalState = nullptr; + GateBremPostStepDoIt *bremProcess = nullptr; + GateGammaEmPostStepDoIt *emProcess = nullptr; + GateplusannihilAtRestDoIt *eplusAnnihilProcess = nullptr; + + if ((container.GetAnnihilationFlag() == "PostStep") && + (fTrackToSplit->GetKineticEnergy() > 0)) { + emProcess = (GateGammaEmPostStepDoIt *)process; + } + if ((container.GetAnnihilationFlag() == "AtRest") || + (fTrackToSplit->GetKineticEnergy() == 0)) { + eplusAnnihilProcess = (GateplusannihilAtRestDoIt *)process; + } + for (int j = 0; j < batchSize; j++) { + G4int NbOfSecondaries = 0; + + G4int count = 0; + while (NbOfSecondaries == 0) { + if (process->GetProcessName() == "eBrem") { + G4Track aTrack = + eBremProcessFinalState(fTrackToSplit, initStep, process); + bremProcess = (GateBremPostStepDoIt *)process; + processFinalState = + bremProcess->GateBremPostStepDoIt::PostStepDoIt(aTrack, *initStep); + } else { + if ((container.GetAnnihilationFlag() == "PostStep") && + (fTrackToSplit->GetKineticEnergy() > 0)) { + processFinalState = + emProcess->PostStepDoIt(*fTrackToSplit, *initStep); + } + if ((container.GetAnnihilationFlag() == "AtRest") || + (fTrackToSplit->GetKineticEnergy() == 0)) { + processFinalState = + eplusAnnihilProcess->GateplusannihilAtRestDoIt::AtRestDoIt( + *fTrackToSplit, *initStep); + } + } + NbOfSecondaries = processFinalState->GetNumberOfSecondaries(); + if (NbOfSecondaries == 0) { + processFinalState->Clear(); + } + count++; + // Security break, in case of infinite loop + if (count > 10000) { + G4ExceptionDescription ed; + ed << " infinite loop detected during the track creation for the " + << process->GetProcessName() << " process" << G4endl; + G4Exception( + "GateLastVertexInteractionSplittingActor::SecondariesSplitting", + "BIAS.LV1", JustWarning, ed); + G4RunManager::GetRunManager()->AbortEvent(); + break; + } + } + + G4int idx = 0; + G4bool IsPushBack = false; + for (int i = 0; i < NbOfSecondaries; i++) { + G4Track *newTrack = processFinalState->GetSecondary(i); + G4ThreeVector momentum = newTrack->GetMomentumDirection(); + + if (!(isnan(momentum[0]))) { + if ((fAngularKill) && (DoesParticleEmittedInSolidAngle( + momentum, fVectorDirector) == false)) { + delete newTrack; + } else if (IsPushBack == true) { + delete newTrack; + } else { + newTrack->SetWeight(fWeight); + newTrack->SetCreatorProcess(process); + // trackVector->emplace_back(newTrack); + fStackManager->PushOneTrack(newTrack); + // delete newTrack; + IsPushBack = true; + } + } else { + delete newTrack; + } + } + processFinalState->Clear(); + } +} + +void GateLastVertexInteractionSplittingActor::CreateNewParticleAtTheLastVertex( + G4Step *initStep, G4Step *step, LastVertexDataContainer theContainer, + G4double batchSize) { + // We retrieve the process associated to the process name to split and we + // split according the process. Since for compton scattering, the gamma is not + // a secondary particles, this one need to have his own splitting function. + + G4String processName = fProcessNameToSplit; + G4int nbOfTrackAlreadyInStack = fStackManager->GetNTotalTrack(); + if ((fProcessToSplit == 0) || (fProcessToSplit == nullptr)) { + SimpleContainer container = theContainer.GetContainerToSplit(); + fProcessToSplit = GetProcessFromProcessName( + container.GetParticleNameToSplit(), processName); + } + + if (processName == "compt") { + ComptonSplitting(initStep, step, fProcessToSplit, theContainer, batchSize); + } + + else if ((processName != "msc") && (processName != "conv")) { + SecondariesSplitting(initStep, step, fProcessToSplit, theContainer, + batchSize); + } + fNumberOfTrackToSimulate = + fStackManager->GetNTotalTrack() - nbOfTrackAlreadyInStack; + fNbOfBatchForExitingParticle++; + if (fNbOfBatchForExitingParticle > fNbOfMaxBatchPerEvent) { + fStackManager->clear(); + fRemovedParticle++; + } + // stackManager->clear(); +} + +void GateLastVertexInteractionSplittingActor::CreateListOfbiasedVolume( + G4LogicalVolume *volume) { + G4int nbOfDaughters = volume->GetNoDaughters(); + if (nbOfDaughters > 0) { + for (int i = 0; i < nbOfDaughters; i++) { + G4String LogicalVolumeName = + volume->GetDaughter(i)->GetLogicalVolume()->GetName(); + G4LogicalVolume *logicalDaughtersVolume = + volume->GetDaughter(i)->GetLogicalVolume(); + if (!(std::find(fListOfBiasedVolume.begin(), fListOfBiasedVolume.end(), + LogicalVolumeName) != fListOfBiasedVolume.end())) + fListOfBiasedVolume.push_back( + volume->GetDaughter(i)->GetLogicalVolume()->GetName()); + CreateListOfbiasedVolume(logicalDaughtersVolume); + } + } +} + +void GateLastVertexInteractionSplittingActor::FillOfDataTree(G4Step *step) { + + G4String processName = "None"; + if (step->GetPostStepPoint()->GetProcessDefinedStep() != 0) + processName = + step->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName(); + + G4String creatorProcessName = "None"; + if (step->GetTrack()->GetCreatorProcess() != 0) + creatorProcessName = + step->GetTrack()->GetCreatorProcess()->GetProcessName(); + + if ((step->GetTrack()->GetParticleDefinition()->GetParticleName() == "e+") && + ((step->GetTrack()->GetTrackStatus() == 1) || + (step->GetTrack()->GetTrackStatus() == 2))) { + processName = "annihil"; + } + + G4String annihilFlag = "None"; + if (processName == "annihil") { + if (step->GetPostStepPoint()->GetProcessDefinedStep() != 0) { + if (processName == + step->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName()) { + annihilFlag = "PostStep"; + } else if (processName != step->GetPostStepPoint() + ->GetProcessDefinedStep() + ->GetProcessName()) { + annihilFlag = "AtRest"; + } + } + } + + if (fIsFirstStep) { + LastVertexDataContainer newContainer = LastVertexDataContainer(); + newContainer.SetTrackID(step->GetTrack()->GetTrackID()); + newContainer.SetParticleName( + step->GetTrack()->GetDefinition()->GetParticleName()); + newContainer.SetCreationProcessName(creatorProcessName); + + if (fTree.empty()) { + fTree.set_head(newContainer); + } + + for (auto it = fTree.begin_post(); it != fTree.end_post(); ++it) { + LastVertexDataContainer container = *it; + G4int trackID = container.GetTrackID(); + + if (step->GetTrack()->GetParentID() == trackID) { + newContainer = container.ContainerFromParentInformation(step); + fTree.append_child(it, newContainer); + break; + } + } + + for (auto it = fTree.begin_post(); it != fTree.end_post(); ++it) { + LastVertexDataContainer container = *it; + G4int trackID = container.GetTrackID(); + if (step->GetTrack()->GetTrackID() == trackID) { + fIterator = it; + break; + } + } + } + + LastVertexDataContainer *container = &(*fIterator); + G4int trackID = container->GetTrackID(); + if ((processName != "Transportation") && (processName != "None") && + (processName != "Rayl")) { + if (step->GetTrack()->GetTrackID() == trackID) { + G4ThreeVector position = step->GetTrack()->GetPosition(); + G4ThreeVector prePosition = step->GetPreStepPoint()->GetPosition(); + G4ThreeVector momentum; + if ((processName == "annihil")) + momentum = step->GetPostStepPoint()->GetMomentumDirection(); + else { + momentum = step->GetPreStepPoint()->GetMomentumDirection(); + } + G4ThreeVector polarization = step->GetPreStepPoint()->GetPolarization(); + G4String particleName = + step->GetTrack()->GetDefinition()->GetParticleName(); + G4double energy = step->GetPreStepPoint()->GetKineticEnergy(); + G4double weight = step->GetTrack()->GetWeight(); + G4int trackStatus = step->GetTrack()->GetTrackStatus(); + G4int nbOfSecondaries = step->GetfSecondary()->size(); + G4double stepLength = step->GetStepLength(); + G4TouchableHandle tch = step->GetTrack()->GetTouchableHandle(); + if (((processName == "annihil"))) { + energy -= (step->GetTotalEnergyDeposit()); + } + SimpleContainer containerToSplit = + SimpleContainer(processName, energy, momentum, position, polarization, + particleName, weight, trackStatus, nbOfSecondaries, + annihilFlag, stepLength, prePosition, tch); + container->SetContainerToSplit(containerToSplit); + container->PushListOfSplittingParameters(containerToSplit); + } + } +} + +G4bool GateLastVertexInteractionSplittingActor::IsParticleExitTheBiasedVolume( + G4Step *step) { + + if ((step->GetPostStepPoint()->GetStepStatus() == 1)) { + G4String logicalVolumeNamePostStep = "None"; + if (step->GetPostStepPoint()->GetPhysicalVolume() != 0) + logicalVolumeNamePostStep = step->GetPostStepPoint() + ->GetPhysicalVolume() + ->GetLogicalVolume() + ->GetName(); + + if (std::find(fListOfVolumeAncestor.begin(), fListOfVolumeAncestor.end(), + logicalVolumeNamePostStep) != fListOfVolumeAncestor.end()) { + return true; + } + /* + else if (std::find(fListOfBiasedVolume.begin(), fListOfBiasedVolume.end(), + logicalVolumeNamePostStep) != fListOfBiasedVolume.end()) { return false; + } + */ + } + if (step->GetPostStepPoint()->GetStepStatus() == 0) + return true; + return false; +} + +G4bool GateLastVertexInteractionSplittingActor::IsTheParticleUndergoesAProcess( + G4Step *step) { + G4String processName = "None"; + G4String particleName = + step->GetTrack()->GetParticleDefinition()->GetParticleName(); + if (step->GetPostStepPoint()->GetProcessDefinedStep() != 0) + processName = + step->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName(); + + if (std::find(fListOfProcessesAccordingParticles[particleName].begin(), + fListOfProcessesAccordingParticles[particleName].end(), + processName) != + fListOfProcessesAccordingParticles[particleName].end()) { + return true; + } + return false; +} + +G4bool GateLastVertexInteractionSplittingActor:: + IsTheParticleUndergoesALossEnergyProcess(G4Step *step) { + if (step->GetPostStepPoint()->GetKineticEnergy() - + step->GetPreStepPoint()->GetKineticEnergy() != + 0) + return true; + return false; +} + +void GateLastVertexInteractionSplittingActor::BeginOfRunAction( + const G4Run *run) { + + fListOfProcessesAccordingParticles["gamma"] = {"compt", "phot", "conv"}; + fListOfProcessesAccordingParticles["e-"] = {"eBrem", "eIoni", "msc"}; + fListOfProcessesAccordingParticles["e+"] = {"eBrem", "eIoni", "msc", + "annihil"}; + + G4LogicalVolume *biasingVolume = + G4LogicalVolumeStore::GetInstance()->GetVolume(fAttachedToVolumeName); + fListOfBiasedVolume.push_back(biasingVolume->GetName()); + CreateListOfbiasedVolume(biasingVolume); + + auto *source = fSourceManager->FindSourceByName("source_vertex"); + fVertexSource = (GateLastVertexSource *)source; + + fCosMaxTheta = std::cos(fMaxTheta); + fStackManager = G4EventManager::GetEventManager()->GetStackManager(); + + if (fRotationVectorDirector) { + G4VPhysicalVolume *physBiasingVolume = + G4PhysicalVolumeStore::GetInstance()->GetVolume(fAttachedToVolumeName); + auto rot = physBiasingVolume->GetObjectRotationValue(); + fVectorDirector = rot * fVectorDirector; + } +} + +void GateLastVertexInteractionSplittingActor::BeginOfEventAction( + const G4Event *event) { + fEventID = event->GetEventID(); + fIsAnnihilAlreadySplit = false; + fNbOfBatchForExitingParticle = 0; + if (fEventID % 50000 == 0) + std::cout << "event ID : " << fEventID << std::endl; + if (fCopyInitStep != 0) { + delete fCopyInitStep; + fCopyInitStep = nullptr; + } + fSplitCounter = 0; + fNumberOfTrackToSimulate = 0; + fKilledBecauseOfProcess = false; + + if (fActiveSource == "source_vertex") { + auto *source = fSourceManager->FindSourceByName(fActiveSource); + GateLastVertexSource *vertexSource = (GateLastVertexSource *)source; + fContainer = vertexSource->GetLastVertexContainer(); + fProcessNameToSplit = vertexSource->GetProcessToSplit(); + if (fProcessToSplit != 0) { + fProcessToSplit = nullptr; + } + if (fTrackToSplit != 0) { + delete fTrackToSplit; + fTrackToSplit = nullptr; + } + fTrackToSplit = CreateATrackFromContainer(fContainer); + if (fTrackToSplit != 0) + fWeight = fTrackToSplit->GetWeight() / fSplittingFactor; + } +} + +void GateLastVertexInteractionSplittingActor::PreUserTrackingAction( + const G4Track *track) { + fToSplit = true; + fIsFirstStep = true; +} + +void GateLastVertexInteractionSplittingActor::SteppingAction(G4Step *step) { + + if (fActiveSource != "source_vertex") { + FillOfDataTree(step); + + if (IsParticleExitTheBiasedVolume(step)) { + if ((fAngularKill == false) || + ((fAngularKill == true) && + (DoesParticleEmittedInSolidAngle( + step->GetTrack()->GetMomentumDirection(), fVectorDirector) == + true))) { + if ((*fIterator).GetContainerToSplit().GetProcessNameToSplit() != + "None") { + fListOfContainer.push_back((*fIterator)); + fNumberOfReplayedParticle++; + } + } + + step->GetTrack()->SetTrackStatus(fStopAndKill); + } + } + + if (fOnlyTree == false) { + if (fActiveSource == "source_vertex") { + + if (fIsFirstStep) { + fTrackID = step->GetTrack()->GetTrackID(); + fEkin = step->GetPostStepPoint()->GetKineticEnergy(); + } else { + if ((fTrackID == step->GetTrack()->GetTrackID()) && + (fEkin != step->GetPreStepPoint()->GetKineticEnergy())) { + fToSplit = false; + } else { + fEkin = step->GetPostStepPoint()->GetKineticEnergy(); + } + } + if (fToSplit) { + G4String creatorProcessName = "None"; + if (step->GetTrack()->GetCreatorProcess() != 0) + creatorProcessName = + step->GetTrack()->GetCreatorProcess()->GetProcessName(); + if (((step->GetTrack()->GetParentID() == 0) && + (step->GetTrack()->GetTrackID() == 1)) || + ((creatorProcessName == "annihil") && + (step->GetTrack()->GetParentID() == 1))) { + + if ((fProcessNameToSplit != "annihil") || + ((fProcessNameToSplit == "annihil") && + (fIsAnnihilAlreadySplit == false))) { + + // FIXME : list of process which are not splitable yet + + if ((fProcessNameToSplit != "msc") && + (fProcessNameToSplit != "conv") && + (fProcessNameToSplit != "eIoni") && + (((fProcessNameToSplit != "phot") || + ((fProcessNameToSplit == "phot") && + (step->GetTrack() + ->GetParticleDefinition() + ->GetParticleName() != "gamma"))))) { + fCopyInitStep = new G4Step(*step); + if (fProcessNameToSplit == "eBrem") { + fCopyInitStep->SetStepLength( + fContainer.GetContainerToSplit().GetStepLength()); + fCopyInitStep->GetPreStepPoint()->SetKineticEnergy( + fContainer.GetContainerToSplit().GetEnergy()); + } + CreateNewParticleAtTheLastVertex(fCopyInitStep, step, fContainer, + fBatchSize); + } + + step->GetTrack()->SetTrackStatus(fKillTrackAndSecondaries); + + if (fProcessNameToSplit == "annihil") { + fIsAnnihilAlreadySplit = true; + } + } + + else if ((fProcessNameToSplit == "annihil") && + (fIsAnnihilAlreadySplit == true)) { + step->GetTrack()->SetTrackStatus(fKillTrackAndSecondaries); + } + + } + + else { + if (fIsFirstStep) { + fNumberOfTrackToSimulate--; + if (fKilledBecauseOfProcess == false) { + fSplitCounter += 1; + } else { + fKilledBecauseOfProcess = false; + } + + if (fSplitCounter > fSplittingFactor) { + step->GetTrack()->SetTrackStatus(fKillTrackAndSecondaries); + fStackManager->clear(); + } + } + + if (IsTheParticleUndergoesALossEnergyProcess(step)) { + step->GetTrack()->SetTrackStatus(fKillTrackAndSecondaries); + fKilledBecauseOfProcess = true; + } + + if (fIsFirstStep) { + if (fSplitCounter <= fSplittingFactor) { + if (fNumberOfTrackToSimulate == 0) { + // CreateNewParticleAtTheLastVertex( + // fCopyInitStep, step, fContainer, + //(fSplittingFactor - fSplitCounter + 1) / fSplittingFactor * + // fBatchSize); + CreateNewParticleAtTheLastVertex(fCopyInitStep, step, + fContainer, fBatchSize); + } + } + } + } + } + } + } + + fIsFirstStep = false; +} + +void GateLastVertexInteractionSplittingActor::EndOfEventAction( + const G4Event *event) { + + if (fActiveSource != "source_vertex") { + + // print_tree(fTree,fTree.begin(),fTree.end()); + fVertexSource->SetNumberOfEventToSimulate(fListOfContainer.size()); + fVertexSource->SetNumberOfGeneratedEvent(0); + fVertexSource->SetListOfVertexToSimulate(fListOfContainer); + fTree.clear(); + fListOfContainer.clear(); + } + + if (fOnlyTree == false) { + + auto *source = fSourceManager->FindSourceByName("source_vertex"); + GateLastVertexSource *vertexSource = (GateLastVertexSource *)source; + if (vertexSource->GetNumberOfGeneratedEvent() < + vertexSource->GetNumberOfEventToSimulate()) { + fSourceManager->SetActiveSourcebyName("source_vertex"); + } + fActiveSource = fSourceManager->GetActiveSourceName(); + } +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... diff --git a/core/opengate_core/opengate_lib/GateLastVertexInteractionSplittingActor.h b/core/opengate_core/opengate_lib/GateLastVertexInteractionSplittingActor.h new file mode 100644 index 000000000..d1a551dcb --- /dev/null +++ b/core/opengate_core/opengate_lib/GateLastVertexInteractionSplittingActor.h @@ -0,0 +1,149 @@ +// +// ******************************************************************** +// * License and Disclaimer * +// * * +// * The Geant4 software is copyright of the Copyright Holders of * +// * the Geant4 Collaboration. It is provided under the terms and * +// * conditions of the Geant4 Software License, included in the file * +// * LICENSE and available at http://cern.ch/geant4/license . These * +// * include a list of copyright holders. * +// * * +// * Neither the authors of this software system, nor their employing * +// * institutes,nor the agencies providing financial support for this * +// * work make any representation or warranty, express or implied, * +// * regarding this software system or assume any liability for its * +// * use. Please see the license in the file LICENSE and URL above * +// * for the full disclaimer and the limitation of liability. * +// * * +// * This code implementation is the result of the scientific and * +// * technical work of the GEANT4 collaboration. * +// * By using, copying, modifying or distributing the software (or * +// * any work based on the software) you agree to acknowledge its * +// * use in resulting scientific publications, and indicate your * +// * acceptance of all terms of the Geant4 Software license. * +// ******************************************************************** +// +/// +/// \file GateLastVertexInteractionSplittingActor.h +/// \brief Definition of the GateLastVertexInteractionSplittingActor class +#ifndef GateLastVertexInteractionSplittingActor_h +#define GateLastVertexInteractionSplittingActor_h 1 + +#include "CLHEP/Vector/ThreeVector.h" +#include "G4ParticleChangeForGamma.hh" +#include "G4StackManager.hh" +#include "G4VEnergyLossProcess.hh" +#include "GateLastVertexSource.h" +#include "GateLastVertexSplittingDataContainer.h" +#include "GateVActor.h" +#include "tree.hh" +#include "tree_util.hh" +#include +#include +using CLHEP::Hep3Vector; + +namespace py = pybind11; + +class GateLastVertexInteractionSplittingActor : public GateVActor { +public: + GateLastVertexInteractionSplittingActor(py::dict &user_info); + virtual ~GateLastVertexInteractionSplittingActor() {} + + G4double fSplittingFactor; + G4bool fAngularKill; + G4bool fRotationVectorDirector; + G4ThreeVector fVectorDirector; + G4double fMaxTheta; + G4double fCosMaxTheta; + G4int fTrackIDOfSplittedTrack = 0; + G4int fParentID = -1; + G4int fEventID; + G4int fEventIDOfSplittedTrack; + G4int fEventIDOfInitialSplittedTrack; + G4int fTrackID; + G4double fEkin; + G4int fTrackIDOfInitialSplittedTrack = 0; + G4int ftmpTrackID; + G4bool fIsFirstStep = true; + G4bool fSuspendForAnnihil = false; + G4double fWeightOfEnteringParticle = 0; + G4double fSplitCounter = 0; + G4bool fToSplit = true; + G4String fActiveSource = "None"; + G4bool fIsAnnihilAlreadySplit = false; + G4int fCounter; + G4bool fKilledBecauseOfProcess = false; + G4bool fFirstSplittedPart = true; + G4bool fOnlyTree = false; + G4double fWeight; + G4double fBatchSize; + G4int fNumberOfTrackToSimulate = 0; + G4int fNbOfBatchForExitingParticle = 0; + G4int fTracksCounts = 0; + GateLastVertexSource *fVertexSource = nullptr; + tree fTree; + tree::post_order_iterator fIterator; + std::vector fListOfContainer; + G4StackManager *fStackManager = nullptr; + G4int fNbOfMaxBatchPerEvent; + G4int fRemovedParticle = 0; + G4int fNumberOfReplayedParticle = 0; + + G4Track *fTrackToSplit = nullptr; + G4Step *fCopyInitStep = nullptr; + G4String fProcessNameToSplit; + G4VProcess *fProcessToSplit; + LastVertexDataContainer fContainer; + + std::vector fTracksToPostpone; + std::map> fListOfProcessesAccordingParticles; + std::map fDataMap; + + std::vector fListOfVolumeAncestor; + std::vector fListOfBiasedVolume; + + std::vector fListOfProcesses = {"compt", "annihil", "eBrem", "conv", + "phot"}; + + void InitializeUserInfo(py::dict &user_info) override; + virtual void SteppingAction(G4Step *) override; + virtual void BeginOfEventAction(const G4Event *) override; + virtual void EndOfEventAction(const G4Event *) override; + virtual void BeginOfRunAction(const G4Run *run) override; + virtual void PreUserTrackingAction(const G4Track *track) override; + + // Pure splitting functions + G4bool DoesParticleEmittedInSolidAngle(G4ThreeVector dir, + G4ThreeVector vectorDirector); + G4Track *CreateComptonTrack(G4ParticleChangeForGamma *, G4Track, G4double); + void ComptonSplitting(G4Step *initStep, G4Step *CurrentStep, + G4VProcess *process, LastVertexDataContainer container, + G4double batchSize); + void SecondariesSplitting(G4Step *initStep, G4Step *CurrentStep, + G4VProcess *process, + LastVertexDataContainer container, + G4double batchSize); + + void CreateNewParticleAtTheLastVertex(G4Step *init, G4Step *current, + LastVertexDataContainer, + G4double batchSize); + G4Track *CreateATrackFromContainer(LastVertexDataContainer container); + G4bool IsTheParticleUndergoesAProcess(G4Step *step); + G4bool IsTheParticleUndergoesALossEnergyProcess(G4Step *step); + G4VProcess *GetProcessFromProcessName(G4String particleName, G4String pName); + G4Track eBremProcessFinalState(G4Track *track, G4Step *step, + G4VProcess *process); + + void FillOfDataTree(G4Step *step); + G4bool IsParticleExitTheBiasedVolume(G4Step *step); + void CreateListOfbiasedVolume(G4LogicalVolume *volume); + void print_tree(const tree &tr, + tree::pre_order_iterator it, + tree::pre_order_iterator end); + inline long GetNumberOfKilledParticles() { return fRemovedParticle; } + inline long GetNumberOfReplayedParticles() { + return fNumberOfReplayedParticle; + } +}; + +#endif diff --git a/core/opengate_core/opengate_lib/GateLastVertexSource.cpp b/core/opengate_core/opengate_lib/GateLastVertexSource.cpp new file mode 100644 index 000000000..c870637a4 --- /dev/null +++ b/core/opengate_core/opengate_lib/GateLastVertexSource.cpp @@ -0,0 +1,109 @@ +/* -------------------------------------------------- + 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 "GateLastVertexSource.h" +#include "G4ParticleTable.hh" +#include "GateHelpersDict.h" +#include + +GateLastVertexSource::GateLastVertexSource() : GateVSource() {} + +GateLastVertexSource::~GateLastVertexSource() {} + +void GateLastVertexSource::InitializeUserInfo(py::dict &user_info) { + GateVSource::InitializeUserInfo(user_info); + // get user info about activity or nb of events + fN = DictGetInt(user_info, "n"); +} + +double GateLastVertexSource::PrepareNextTime(double current_simulation_time) { + + /* + // If all N events have been generated, we stop (negative time) + if (fNumberOfGeneratedEvents >= fN){ + std::cout << "LV: "<< -1<<" "<< fNumberOfGeneratedEvents <<" "<< + fN<= fN) + return -1; + + return fStartTime + 1; +} + +void GateLastVertexSource::PrepareNextRun() { + // The following compute the global transformation from + // the local volume (mother) to the world + GateVSource::PrepareNextRun(); + + // The global transformation to apply for the current RUN is known in : + // fGlobalTranslation & fGlobalRotation + + // init the number of generated events (here, for each run) + fNumberOfGeneratedEvents = 0; +} + +void GateLastVertexSource::GenerateOnePrimary(G4Event *event, + double current_simulation_time, + G4int idx) { + + if (fNumberOfGeneratedEvents >= fN) { + auto *particle_table = G4ParticleTable::GetParticleTable(); + auto *fParticleDefinition = particle_table->FindParticle("geantino"); + auto *particle = new G4PrimaryParticle(fParticleDefinition); + particle->SetKineticEnergy(0); + particle->SetMomentumDirection({1, 0, 0}); + particle->SetWeight(1); + auto *vertex = new G4PrimaryVertex({0, 0, 0}, current_simulation_time); + vertex->SetPrimary(particle); + event->AddPrimaryVertex(vertex); + } else { + + SimpleContainer containerToSplit = + fListOfContainer[idx].GetContainerToSplit(); + G4double energy = containerToSplit.GetEnergy(); + if (energy < 0) { + energy = 0; + } + fContainer = fListOfContainer[idx]; + G4ThreeVector position = containerToSplit.GetVertexPosition(); + G4ThreeVector momentum = containerToSplit.GetMomentum(); + G4String particleName = containerToSplit.GetParticleNameToSplit(); + 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); + particle->SetKineticEnergy(energy); + particle->SetMomentumDirection(momentum); + particle->SetWeight(weight); + auto *vertex = new G4PrimaryVertex(position, current_simulation_time); + vertex->SetPrimary(particle); + event->AddPrimaryVertex(vertex); + } +} + +void GateLastVertexSource::GeneratePrimaries(G4Event *event, + double current_simulation_time) { + + GenerateOnePrimary(event, current_simulation_time, fNumberOfGeneratedEvents); + fNumberOfGeneratedEvents++; + if (fNumberOfGeneratedEvents == fListOfContainer.size()) { + fListOfContainer.clear(); + } +} diff --git a/core/opengate_core/opengate_lib/GateLastVertexSource.h b/core/opengate_core/opengate_lib/GateLastVertexSource.h new file mode 100644 index 000000000..eebf43e90 --- /dev/null +++ b/core/opengate_core/opengate_lib/GateLastVertexSource.h @@ -0,0 +1,72 @@ +/* -------------------------------------------------- + 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 GateLastVertexSource_h +#define GateLastVertexSource_h + +#include "GateAcceptanceAngleTesterManager.h" +#include "GateLastVertexSplittingDataContainer.h" +#include "GateSingleParticleSource.h" +#include "GateVSource.h" +#include + +namespace py = pybind11; + +/* + This is NOT a real source type but a template to help writing your own source + type. Copy-paste this file with a different name ("MyNewSource.hh") and start + building. You also need to copy : GateLastVertexSource.hh + GateLastVertexSource.cpp pyGateLastVertexSource.cpp And add the source + declaration in opengate_core.cpp + */ + +class GateLastVertexSource : public GateVSource { + +public: + GateLastVertexSource(); + + ~GateLastVertexSource() override; + + void InitializeUserInfo(py::dict &user_info) override; + + double PrepareNextTime(double current_simulation_time) override; + + void PrepareNextRun() override; + + void GeneratePrimaries(G4Event *event, double time) override; + + void GenerateOnePrimary(G4Event *event, double time, G4int idx); + + void SetListOfVertexToSimulate(std::vector list) { + fListOfContainer = list; + } + + void SetNumberOfGeneratedEvent(G4int nbEvent) { + fNumberOfGeneratedEvents = nbEvent; + } + + void SetNumberOfEventToSimulate(G4int N) { fN = N; } + + G4int GetNumberOfEventToSimulate() { return fN; } + + G4int GetNumberOfGeneratedEvent() { return fNumberOfGeneratedEvents; } + + G4String GetProcessToSplit() { return fProcessToSplit; } + + LastVertexDataContainer GetLastVertexContainer() { return fContainer; } + +protected: + G4int fNumberOfGeneratedEvents = 0; + G4int fN = 0; + double fFloatValue; + std::vector fVectorValue; + std::vector fListOfContainer; + G4String fProcessToSplit = "None"; + LastVertexDataContainer fContainer; +}; + +#endif // GateLastVertexSource_h diff --git a/core/opengate_core/opengate_lib/GateLastVertexSplittingDataContainer.h b/core/opengate_core/opengate_lib/GateLastVertexSplittingDataContainer.h new file mode 100644 index 000000000..4264766d8 --- /dev/null +++ b/core/opengate_core/opengate_lib/GateLastVertexSplittingDataContainer.h @@ -0,0 +1,146 @@ +// +// ******************************************************************** +// * License and Disclaimer * +// * * +// * The Geant4 software is copyright of the Copyright Holders of * +// * the Geant4 Collaboration. It is provided under the terms and * +// * conditions of the Geant4 Software License, included in the file * +// * LICENSE and available at http://cern.ch/geant4/license . These * +// * include a list of copyright holders. * +// * * +// * Neither the authors of this software system, nor their employing * +// * institutes,nor the agencies providing financial support for this * +// * work make any representation or warranty, express or implied, * +// * regarding this software system or assume any liability for its * +// * use. Please see the license in the file LICENSE and URL above * +// * for the full disclaimer and the limitation of liability. * +// * * +// * This code implementation is the result of the scientific and * +// * technical work of the GEANT4 collaboration. * +// * By using, copying, modifying or distributing the software (or * +// * any work based on the software) you agree to acknowledge its * +// * use in resulting scientific publications, and indicate your * +// * acceptance of all terms of the Geant4 Software license. * +// ******************************************************************** +// +/// +#ifndef LastVertexDataContainer_h +#define LastVertexDataContainer_h + +#include "G4Electron.hh" +#include "G4EmBiasingManager.hh" +#include "G4EmParameters.hh" +#include "G4EntanglementAuxInfo.hh" +#include "G4Gamma.hh" +#include "G4MaterialCutsCouple.hh" +#include "G4PhysicalConstants.hh" +#include "G4PhysicsModelCatalog.hh" +#include "G4Positron.hh" +#include "G4Track.hh" +#include "G4VEmProcess.hh" +#include "G4VEnergyLossProcess.hh" +#include "G4VParticleChange.hh" +#include "G4eeToTwoGammaModel.hh" +#include "G4eplusAnnihilation.hh" +#include "G4eplusAnnihilationEntanglementClipBoard.hh" +#include "GateLastVertexSplittingSimpleContainer.h" +#include + +class LastVertexDataContainer { + +public: + LastVertexDataContainer() {} + + ~LastVertexDataContainer() {} + + void SetTrackID(G4int trackID) { fTrackID = trackID; } + + G4int GetTrackID() { return fTrackID; } + + void SetParticleName(G4String name) { fParticleName = name; } + + G4String GetParticleName() { return fParticleName; } + + void SetCreationProcessName(G4String creationProcessName) { + fCreationProcessName = creationProcessName; + } + + G4String GetCreationProcessName() { return fCreationProcessName; } + + void SetContainerToSplit(SimpleContainer container) { + fContainerToSplit = container; + } + + SimpleContainer GetContainerToSplit() { return fContainerToSplit; } + + void PushListOfSplittingParameters(SimpleContainer container) { + fVectorOfContainerToSplit.emplace_back(container); + } + + LastVertexDataContainer ContainerFromParentInformation(G4Step *step) { + LastVertexDataContainer aContainer = LastVertexDataContainer(); + + aContainer.fTrackID = step->GetTrack()->GetTrackID(); + aContainer.fParticleName = + step->GetTrack()->GetDefinition()->GetParticleName(); + if (this->fContainerToSplit.GetProcessNameToSplit() != "None") { + if (this->fVectorOfContainerToSplit.size() != 0) { + G4ThreeVector vertexPosition = step->GetTrack()->GetVertexPosition(); + for (int i = 0; i < this->fVectorOfContainerToSplit.size(); i++) { + if (vertexPosition == + this->fVectorOfContainerToSplit[i].GetVertexPosition()) { + SimpleContainer tmpContainer = this->fVectorOfContainerToSplit[i]; + // std::cout<<"1 "<fContainerToSplit; + // std::cout<<"2 "< fVectorOfContainerToSplit; +}; + +#endif diff --git a/core/opengate_core/opengate_lib/GateLastVertexSplittingPostStepDoIt.h b/core/opengate_core/opengate_lib/GateLastVertexSplittingPostStepDoIt.h new file mode 100644 index 000000000..ba7b9b366 --- /dev/null +++ b/core/opengate_core/opengate_lib/GateLastVertexSplittingPostStepDoIt.h @@ -0,0 +1,103 @@ +// +// ******************************************************************** +// * License and Disclaimer * +// * * +// * The Geant4 software is copyright of the Copyright Holders of * +// * the Geant4 Collaboration. It is provided under the terms and * +// * conditions of the Geant4 Software License, included in the file * +// * LICENSE and available at http://cern.ch/geant4/license . These * +// * include a list of copyright holders. * +// * * +// * Neither the authors of this software system, nor their employing * +// * institutes,nor the agencies providing financial support for this * +// * work make any representation or warranty, express or implied, * +// * regarding this software system or assume any liability for its * +// * use. Please see the license in the file LICENSE and URL above * +// * for the full disclaimer and the limitation of liability. * +// * * +// * This code implementation is the result of the scientific and * +// * technical work of the GEANT4 collaboration. * +// * By using, copying, modifying or distributing the software (or * +// * any work based on the software) you agree to acknowledge its * +// * use in resulting scientific publications, and indicate your * +// * acceptance of all terms of the Geant4 Software license. * +// ******************************************************************** +// +/// +#ifndef GateLastVertexSplittingPostStepDoIt_h +#define GateLastVertexSplittingPostStepDoIt_h + +#include "G4Electron.hh" +#include "G4EmBiasingManager.hh" +#include "G4EmParameters.hh" +#include "G4EntanglementAuxInfo.hh" +#include "G4Gamma.hh" +#include "G4MaterialCutsCouple.hh" +#include "G4PhysicalConstants.hh" +#include "G4PhysicsModelCatalog.hh" +#include "G4Positron.hh" +#include "G4VEmProcess.hh" +#include "G4VEnergyLossProcess.hh" +#include "G4VParticleChange.hh" +#include "G4eeToTwoGammaModel.hh" +#include "G4eplusAnnihilation.hh" +#include "G4eplusAnnihilationEntanglementClipBoard.hh" +#include + +class GateBremPostStepDoIt : public G4VEnergyLossProcess { +public: + GateBremPostStepDoIt(); + + ~GateBremPostStepDoIt(); + + virtual G4VParticleChange *PostStepDoIt(const G4Track &track, + const G4Step &step) override { + const G4MaterialCutsCouple *couple = + step.GetPreStepPoint()->GetMaterialCutsCouple(); + currentCouple = couple; + G4VParticleChange *particleChange = + G4VEnergyLossProcess::PostStepDoIt(track, step); + return particleChange; + } + + virtual G4VParticleChange *AlongStepDoIt(const G4Track &track, + const G4Step &step) override { + G4VParticleChange *particleChange = + G4VEnergyLossProcess::AlongStepDoIt(track, step); + return particleChange; + } +}; + +class GateGammaEmPostStepDoIt : public G4VEmProcess { +public: + GateGammaEmPostStepDoIt(); + + ~GateGammaEmPostStepDoIt(); + + virtual G4VParticleChange *PostStepDoIt(const G4Track &track, + const G4Step &step) override { + const G4MaterialCutsCouple *couple = + step.GetPreStepPoint()->GetMaterialCutsCouple(); + currentCouple = couple; + G4VParticleChange *particleChange = G4VEmProcess::PostStepDoIt(track, step); + return particleChange; + } +}; + +class GateplusannihilAtRestDoIt : public G4eplusAnnihilation { +public: + GateplusannihilAtRestDoIt(); + ~GateplusannihilAtRestDoIt(); + + virtual G4VParticleChange *AtRestDoIt(const G4Track &track, + const G4Step &step) override + // Performs the e+ e- annihilation when both particles are assumed at rest. + { + G4Track copyTrack = G4Track(track); + copyTrack.SetStep(&step); + G4VParticleChange *particleChange = + G4eplusAnnihilation::AtRestDoIt(copyTrack, step); + return particleChange; + } +}; +#endif diff --git a/core/opengate_core/opengate_lib/GateLastVertexSplittingSimpleContainer.h b/core/opengate_core/opengate_lib/GateLastVertexSplittingSimpleContainer.h new file mode 100644 index 000000000..3980916de --- /dev/null +++ b/core/opengate_core/opengate_lib/GateLastVertexSplittingSimpleContainer.h @@ -0,0 +1,163 @@ +// +// ******************************************************************** +// * License and Disclaimer * +// * * +// * The Geant4 software is copyright of the Copyright Holders of * +// * the Geant4 Collaboration. It is provided under the terms and * +// * conditions of the Geant4 Software License, included in the file * +// * LICENSE and available at http://cern.ch/geant4/license . These * +// * include a list of copyright holders. * +// * * +// * Neither the authors of this software system, nor their employing * +// * institutes,nor the agencies providing financial support for this * +// * work make any representation or warranty, express or implied, * +// * regarding this software system or assume any liability for its * +// * use. Please see the license in the file LICENSE and URL above * +// * for the full disclaimer and the limitation of liability. * +// * * +// * This code implementation is the result of the scientific and * +// * technical work of the GEANT4 collaboration. * +// * By using, copying, modifying or distributing the software (or * +// * any work based on the software) you agree to acknowledge its * +// * use in resulting scientific publications, and indicate your * +// * acceptance of all terms of the Geant4 Software license. * +// ******************************************************************** +// +/// +#ifndef SimpleContainer_h +#define SimpleContainer_h + +#include "G4Electron.hh" +#include "G4EmBiasingManager.hh" +#include "G4EmParameters.hh" +#include "G4EntanglementAuxInfo.hh" +#include "G4Gamma.hh" +#include "G4MaterialCutsCouple.hh" +#include "G4PhysicalConstants.hh" +#include "G4PhysicsModelCatalog.hh" +#include "G4Positron.hh" +#include "G4ReferenceCountedHandle.hh" +#include "G4TouchableHandle.hh" +#include "G4Track.hh" +#include "G4VEmProcess.hh" +#include "G4VEnergyLossProcess.hh" +#include "G4VParticleChange.hh" +#include "G4eeToTwoGammaModel.hh" +#include "G4eplusAnnihilation.hh" + +class SimpleContainer { + +public: + SimpleContainer(G4String processName, G4double energy, G4ThreeVector momentum, + G4ThreeVector position, G4ThreeVector polarization, + G4String name, G4double weight, G4int trackStatus, + G4int nbSec, G4String flag, G4double length, + G4ThreeVector prePos, G4TouchableHandle aTouchable) { + + fProcessNameToSplit = processName; + fEnergyToSplit = energy; + fMomentumToSplit = momentum; + fPositionToSplit = position; + fPolarizationToSplit = polarization; + fParticleNameToSplit = name; + fWeightToSplit = weight; + fTrackStatusToSplit = trackStatus; + fNumberOfSecondariesToSplit = nbSec; + fAnnihilProcessFlag = flag; + fStepLength = length; + fPrePosition = prePos; + fpTouchable.push_back(aTouchable); + } + + SimpleContainer() {} + + ~SimpleContainer() {} + + void SetProcessNameToSplit(G4String processName) { + fProcessNameToSplit = processName; + } + + G4String GetProcessNameToSplit() { return fProcessNameToSplit; } + + void SetEnergy(G4double energy) { fEnergyToSplit = energy; } + + G4double GetEnergy() { return fEnergyToSplit; } + + void SetWeight(G4double weight) { fWeightToSplit = weight; } + + G4double GetWeight() { return fWeightToSplit; } + + void SetPolarization(G4ThreeVector polarization) { + fPolarizationToSplit = polarization; + } + + G4ThreeVector GetPolarization() { return fPolarizationToSplit; } + + void SetMomentum(G4ThreeVector momentum) { fMomentumToSplit = momentum; } + + G4ThreeVector GetMomentum() { return fMomentumToSplit; } + + void SetVertexPosition(G4ThreeVector position) { + fPositionToSplit = position; + } + + G4ThreeVector GetVertexPosition() { return fPositionToSplit; } + + void SetParticleNameToSplit(G4String name) { fParticleNameToSplit = name; } + + G4String GetParticleNameToSplit() { return fParticleNameToSplit; } + + void SetTrackStatus(G4int trackStatus) { fTrackStatusToSplit = trackStatus; } + + G4int GetTrackStatus() { return fTrackStatusToSplit; } + + void SetNbOfSecondaries(G4int nbSec) { fNumberOfSecondariesToSplit = nbSec; } + + G4int GetNbOfSecondaries() { return fNumberOfSecondariesToSplit; } + + void SetAnnihilationFlag(G4String flag) { fAnnihilProcessFlag = flag; } + + G4String GetAnnihilationFlag() { return fAnnihilProcessFlag; } + + void SetStepLength(G4double length) { fStepLength = length; } + + G4double GetStepLength() { return fStepLength; } + + void SetPrePositionToSplit(G4ThreeVector prePos) { fPrePosition = prePos; } + + G4ThreeVector GetPrePositionToSplit() { return fPrePosition; } + + G4TouchableHandle GetTouchableHandle() { return fpTouchable[0]; } + + void DumpInfoToSplit() { + std::cout << "Particle name of the particle to split: " + << fParticleNameToSplit << std::endl; + std::cout << "Kinetic Energy of the particle to split: " << fEnergyToSplit + << std::endl; + std::cout << "Momentum of the particle to split: " << fMomentumToSplit + << std::endl; + std::cout << "Initial position of the particle to split: " + << fPositionToSplit << std::endl; + std::cout << "ProcessNameToSplit: " << fProcessNameToSplit << std::endl; + std::cout << " " << std::endl; + } + + // G4TouchableHandle fpTouchable; + +private: + G4String fParticleNameToSplit = "None"; + G4String fProcessNameToSplit = "None"; + G4double fEnergyToSplit = 0; + G4ThreeVector fMomentumToSplit; + G4ThreeVector fPositionToSplit; + G4ThreeVector fPolarizationToSplit; + G4double fWeightToSplit; + G4int fTrackStatusToSplit; + G4int fNumberOfSecondariesToSplit; + G4String fAnnihilProcessFlag; + G4double fStepLength; + G4ThreeVector fPrePosition; + std::vector fpTouchable; +}; + +#endif diff --git a/core/opengate_core/opengate_lib/GateSourceManager.cpp b/core/opengate_core/opengate_lib/GateSourceManager.cpp index d183567e9..61dfb8361 100644 --- a/core/opengate_core/opengate_lib/GateSourceManager.cpp +++ b/core/opengate_core/opengate_lib/GateSourceManager.cpp @@ -255,6 +255,7 @@ void GateSourceManager::PrepareNextSource() { l.fNextSimulationTime = t; } } + // If no next time in the current interval, active source is NULL } diff --git a/core/opengate_core/opengate_lib/GateSourceManager.h b/core/opengate_core/opengate_lib/GateSourceManager.h index 3c8bc00e1..55e75ff8a 100644 --- a/core/opengate_core/opengate_lib/GateSourceManager.h +++ b/core/opengate_core/opengate_lib/GateSourceManager.h @@ -65,6 +65,21 @@ class GateSourceManager : public G4VUserPrimaryGeneratorAction { // Return a source GateVSource *FindSourceByName(std::string name) const; + G4String GetActiveSourceName() { + auto &l = fThreadLocalData.Get(); + if (l.fNextActiveSource != 0) { + G4String name = l.fNextActiveSource->fName; + return name; + } + return "None"; + } + + void SetActiveSourcebyName(G4String sourceName) { + auto &l = fThreadLocalData.Get(); + auto *source = FindSourceByName(sourceName); + l.fNextActiveSource = source; + } + // [available on py side] start the simulation, master thread only void StartMasterThread(); diff --git a/core/opengate_core/opengate_lib/pyGateKillAccordingParticleNameActor.cpp b/core/opengate_core/opengate_lib/pyGateKillAccordingParticleNameActor.cpp new file mode 100644 index 000000000..364261cec --- /dev/null +++ b/core/opengate_core/opengate_lib/pyGateKillAccordingParticleNameActor.cpp @@ -0,0 +1,23 @@ +/* -------------------------------------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + -------------------------------------------------- */ + +#include + +namespace py = pybind11; + +#include "GateKillAccordingParticleNameActor.h" + +void init_GateKillAccordingParticleNameActor(py::module &m) { + py::class_, + GateVActor>(m, "GateKillAccordingParticleNameActor") + .def(py::init()) + .def_readwrite("fListOfVolumeAncestor", + &GateKillAccordingParticleNameActor::fListOfVolumeAncestor) + .def("GetNumberOfKilledParticles", + &GateKillAccordingParticleNameActor::GetNumberOfKilledParticles); +} diff --git a/core/opengate_core/opengate_lib/pyGateKillNonInteractingParticleActor.cpp b/core/opengate_core/opengate_lib/pyGateKillNonInteractingParticleActor.cpp new file mode 100644 index 000000000..d454f60a4 --- /dev/null +++ b/core/opengate_core/opengate_lib/pyGateKillNonInteractingParticleActor.cpp @@ -0,0 +1,24 @@ +/* -------------------------------------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + -------------------------------------------------- */ + +#include + +namespace py = pybind11; + +#include "GateKillNonInteractingParticleActor.h" + +void init_GateKillNonInteractingParticleActor(py::module &m) { + py::class_, + GateVActor>(m, "GateKillNonInteractingParticleActor") + .def_readwrite( + "fListOfVolumeAncestor", + &GateKillNonInteractingParticleActor::fListOfVolumeAncestor) + .def_readwrite("number_of_killed_particles", + &GateKillNonInteractingParticleActor::fNbOfKilledParticles) + .def(py::init()); +} diff --git a/core/opengate_core/opengate_lib/pyGateLastVertexInteractionSplittingActor.cpp b/core/opengate_core/opengate_lib/pyGateLastVertexInteractionSplittingActor.cpp new file mode 100644 index 000000000..9ac022e2a --- /dev/null +++ b/core/opengate_core/opengate_lib/pyGateLastVertexInteractionSplittingActor.cpp @@ -0,0 +1,27 @@ +/* -------------------------------------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + -------------------------------------------------- */ +#include + +namespace py = pybind11; +#include "GateLastVertexInteractionSplittingActor.h" + +void init_GateLastVertexInteractionSplittingActor(py::module &m) { + + py::class_< + GateLastVertexInteractionSplittingActor, GateVActor, + std::unique_ptr>( + m, "GateLastVertexInteractionSplittingActor") + .def_readwrite( + "fListOfVolumeAncestor", + &GateLastVertexInteractionSplittingActor::fListOfVolumeAncestor) + .def(py::init()) + .def("GetNumberOfKilledParticles", + &GateLastVertexInteractionSplittingActor::GetNumberOfKilledParticles) + .def("GetNumberOfReplayedParticles", + &GateLastVertexInteractionSplittingActor:: + GetNumberOfReplayedParticles); +} diff --git a/core/opengate_core/opengate_lib/pyGateLastVertexSource.cpp b/core/opengate_core/opengate_lib/pyGateLastVertexSource.cpp new file mode 100644 index 000000000..0cb59f146 --- /dev/null +++ b/core/opengate_core/opengate_lib/pyGateLastVertexSource.cpp @@ -0,0 +1,22 @@ +/* -------------------------------------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + -------------------------------------------------- */ + +#include + +namespace py = pybind11; + +#include "GateLastVertexSource.h" + +void init_GateLastVertexSource(py::module &m) { + + py::class_(m, "GateLastVertexSource") + .def(py::init()) + .def("InitializeUserInfo", &GateLastVertexSource::InitializeUserInfo) + // If needed: add your own class functions that will be accessible from + // python side. + ; +} diff --git a/core/opengate_core/opengate_lib/tree.hh b/core/opengate_core/opengate_lib/tree.hh new file mode 100644 index 000000000..17e501f48 --- /dev/null +++ b/core/opengate_core/opengate_lib/tree.hh @@ -0,0 +1,3490 @@ + +// STL-like templated tree class. +// +// Copyright (C) 2001-2024 Kasper Peeters +// Distributed under the GNU General Public License version 3. +// +// Special permission to use tree.hh under the conditions of a +// different license can be requested from the author. + +/** \mainpage tree.hh + \author Kasper Peeters + \version 3.20 + \date 2024-04-12 + \see http://github.com/kpeeters/tree.hh/ + + The tree.hh library for C++ provides an STL-like container class + for n-ary trees, templated over the data stored at the + nodes. Various types of iterators are provided (post-order, + pre-order, and others). Where possible the access methods are + compatible with the STL or alternative algorithms are + available. +*/ + +#ifndef tree_hh_ +#define tree_hh_ + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +/// A node in the tree, combining links to other nodes as well as the actual +/// data. +template +class tree_node_ { // size: 5*4=20 bytes (on 32 bit arch), can be reduced by 8. +public: + tree_node_(); + tree_node_(const T &); + tree_node_(T &&); + + tree_node_ *parent; + tree_node_ *first_child, *last_child; + tree_node_ *prev_sibling, *next_sibling; + T data; +}; + +template +tree_node_::tree_node_() + : parent(0), first_child(0), last_child(0), prev_sibling(0), + next_sibling(0) {} + +template +tree_node_::tree_node_(const T &val) + : parent(0), first_child(0), last_child(0), prev_sibling(0), + next_sibling(0), data(val) {} + +template +tree_node_::tree_node_(T &&val) + : parent(0), first_child(0), last_child(0), prev_sibling(0), + next_sibling(0), data(val) {} + +// Throw an exception with a stacktrace. + +// template +// void throw_with_trace(const E& e) +// { +// throw boost::enable_error_info(e) +// << traced(boost::stacktrace::stacktrace()); +// } + +class navigation_error : public std::logic_error { +public: + navigation_error(const std::string &s) : std::logic_error(s) { + // assert(1==0); + // std::ostringstream str; + // std::cerr << boost::stacktrace::stacktrace() << + // std::endl; str << boost::stacktrace::stacktrace(); + // stacktrace=str.str(); + } + + // virtual const char *what() const noexcept override + // { + // return (std::logic_error::what()+std::string("; + //")+stacktrace).c_str(); + // } + // + // std::string stacktrace; +}; + +template >> +class tree { +protected: + typedef tree_node_ tree_node; + +public: + /// Value of the data stored at a node. + typedef T value_type; + + class iterator_base; + class pre_order_iterator; + class post_order_iterator; + class sibling_iterator; + class leaf_iterator; + + tree(); // empty constructor + tree(const T &); // constructor setting given element as head + tree(const iterator_base &); + tree(const tree &); // copy constructor + tree(tree &&); // move constructor + ~tree(); + tree & + operator=(const tree &); // copy assignment + tree & + operator=(tree &&); // move assignment + + /// Base class for iterators, only pointers stored, no traversal logic. +#ifdef __SGI_STL_PORT + class iterator_base : public stlport::bidirectional_iterator { +#else + class iterator_base { +#endif + public: + typedef T value_type; + typedef T *pointer; + typedef T &reference; + typedef size_t size_type; + typedef ptrdiff_t difference_type; + typedef std::bidirectional_iterator_tag iterator_category; + + iterator_base(); + iterator_base(tree_node *); + + T &operator*() const; + T *operator->() const; + + /// When called, the next increment/decrement skips children of this node. + void skip_children(); + void skip_children(bool skip); + /// Number of children of the node pointed to by the iterator. + unsigned int number_of_children() const; + + sibling_iterator begin() const; + sibling_iterator end() const; + + tree_node *node; + + protected: + bool skip_current_children_; + }; + + /// Depth-first iterator, first accessing the node, then its children. + class pre_order_iterator : public iterator_base { + public: + pre_order_iterator(); + pre_order_iterator(tree_node *); + pre_order_iterator(const iterator_base &); + pre_order_iterator(const sibling_iterator &); + + bool operator==(const pre_order_iterator &) const; + bool operator!=(const pre_order_iterator &) const; + pre_order_iterator &operator++(); + pre_order_iterator &operator--(); + pre_order_iterator operator++(int); + pre_order_iterator operator--(int); + pre_order_iterator &operator+=(unsigned int); + pre_order_iterator &operator-=(unsigned int); + + pre_order_iterator &next_skip_children(); + }; + + /// Depth-first iterator, first accessing the children, then the node itself. + class post_order_iterator : public iterator_base { + public: + post_order_iterator(); + post_order_iterator(tree_node *); + post_order_iterator(const iterator_base &); + post_order_iterator(const sibling_iterator &); + + bool operator==(const post_order_iterator &) const; + bool operator!=(const post_order_iterator &) const; + post_order_iterator &operator++(); + post_order_iterator &operator--(); + post_order_iterator operator++(int); + post_order_iterator operator--(int); + post_order_iterator &operator+=(unsigned int); + post_order_iterator &operator-=(unsigned int); + + /// Set iterator to the first child as deep as possible down the tree. + void descend_all(); + }; + + /// Breadth-first iterator, using a queue + class breadth_first_queued_iterator : public iterator_base { + public: + breadth_first_queued_iterator(); + breadth_first_queued_iterator(tree_node *); + breadth_first_queued_iterator(const iterator_base &); + + bool operator==(const breadth_first_queued_iterator &) const; + bool operator!=(const breadth_first_queued_iterator &) const; + breadth_first_queued_iterator &operator++(); + breadth_first_queued_iterator operator++(int); + breadth_first_queued_iterator &operator+=(unsigned int); + + private: + std::queue traversal_queue; + }; + + /// The default iterator types throughout the tree class. + typedef pre_order_iterator iterator; + typedef breadth_first_queued_iterator breadth_first_iterator; + + /// Iterator which traverses only the nodes at a given depth from the root. + class fixed_depth_iterator : public iterator_base { + public: + fixed_depth_iterator(); + fixed_depth_iterator(tree_node *); + fixed_depth_iterator(const iterator_base &); + fixed_depth_iterator(const sibling_iterator &); + fixed_depth_iterator(const fixed_depth_iterator &); + + void swap(fixed_depth_iterator &, fixed_depth_iterator &); + fixed_depth_iterator &operator=(fixed_depth_iterator); + + bool operator==(const fixed_depth_iterator &) const; + bool operator!=(const fixed_depth_iterator &) const; + fixed_depth_iterator &operator++(); + fixed_depth_iterator &operator--(); + fixed_depth_iterator operator++(int); + fixed_depth_iterator operator--(int); + fixed_depth_iterator &operator+=(unsigned int); + fixed_depth_iterator &operator-=(unsigned int); + + tree_node *top_node; + }; + + /// Iterator which traverses only the nodes which are siblings of each other. + class sibling_iterator : public iterator_base { + public: + sibling_iterator(); + sibling_iterator(tree_node *); + sibling_iterator(const sibling_iterator &); + sibling_iterator(const iterator_base &); + + void swap(sibling_iterator &, sibling_iterator &); + sibling_iterator &operator=(sibling_iterator); + + bool operator==(const sibling_iterator &) const; + bool operator!=(const sibling_iterator &) const; + sibling_iterator &operator++(); + sibling_iterator &operator--(); + sibling_iterator operator++(int); + sibling_iterator operator--(int); + sibling_iterator &operator+=(unsigned int); + sibling_iterator &operator-=(unsigned int); + + tree_node *range_first() const; + tree_node *range_last() const; + tree_node *parent_; + + private: + void set_parent_(); + }; + + /// Iterator which traverses only the leaves. + class leaf_iterator : public iterator_base { + public: + leaf_iterator(); + leaf_iterator(tree_node *, tree_node *top = 0); + leaf_iterator(const sibling_iterator &); + leaf_iterator(const iterator_base &); + + bool operator==(const leaf_iterator &) const; + bool operator!=(const leaf_iterator &) const; + leaf_iterator &operator++(); + leaf_iterator &operator--(); + leaf_iterator operator++(int); + leaf_iterator operator--(int); + leaf_iterator &operator+=(unsigned int); + leaf_iterator &operator-=(unsigned int); + + private: + tree_node *top_node; + }; + + /// Return iterator to the beginning of the tree. + inline pre_order_iterator begin() const; + /// Return iterator to the end of the tree. + inline pre_order_iterator end() const; + /// Return post-order iterator to the beginning of the tree. + post_order_iterator begin_post() const; + /// Return post-order end iterator of the tree. + post_order_iterator end_post() const; + /// Return fixed-depth iterator to the first node at a given depth from the + /// given iterator. If 'walk_back=true', a depth=0 iterator will be taken from + /// the beginning of the sibling range, not the current node. + fixed_depth_iterator begin_fixed(const iterator_base &, unsigned int, + bool walk_back = true) const; + /// Return fixed-depth end iterator. + fixed_depth_iterator end_fixed(const iterator_base &, unsigned int) const; + /// Return breadth-first iterator to the first node at a given depth. + breadth_first_queued_iterator begin_breadth_first() const; + /// Return breadth-first end iterator. + breadth_first_queued_iterator end_breadth_first() const; + /// Return sibling iterator to the first child of given node. + static sibling_iterator begin(const iterator_base &); + /// Return sibling end iterator for children of given node. + static sibling_iterator end(const iterator_base &); + /// Return leaf iterator to the first leaf of the tree. + leaf_iterator begin_leaf() const; + /// Return leaf end iterator for entire tree. + leaf_iterator end_leaf() const; + /// Return leaf iterator to the first leaf of the subtree at the given node. + leaf_iterator begin_leaf(const iterator_base &top) const; + /// Return leaf end iterator for the subtree at the given node. + leaf_iterator end_leaf(const iterator_base &top) const; + + typedef std::vector path_t; + /// Return a path (to be taken from the 'top' node) corresponding to a node in + /// the tree. The first integer in path_t is the number of steps you need to + /// go 'right' in the sibling chain (so 0 if we go straight to the children). + path_t path_from_iterator(const iterator_base &iter, + const iterator_base &top) const; + /// Return an iterator given a path from the 'top' node. + iterator iterator_from_path(const path_t &, const iterator_base &top) const; + + /// Return iterator to the parent of a node. Throws a `navigation_error` if + /// the node does not have a parent. + template static iter parent(iter); + /// Return iterator to the previous sibling of a node. + template static iter previous_sibling(iter); + /// Return iterator to the next sibling of a node. + template static iter next_sibling(iter); + /// Return iterator to the next node at a given depth. + template iter next_at_same_depth(iter) const; + + /// Erase all nodes of the tree. + void clear(); + /// Erase element at position pointed to by iterator, return incremented + /// iterator. + template iter erase(iter); + /// Erase all children of the node pointed to by iterator. + void erase_children(const iterator_base &); + /// Erase all siblings to the right of the iterator. + void erase_right_siblings(const iterator_base &); + /// Erase all siblings to the left of the iterator. + void erase_left_siblings(const iterator_base &); + + /// Insert empty node as last/first child of node pointed to by position. + template iter append_child(iter position); + template iter prepend_child(iter position); + /// Insert node as last/first child of node pointed to by position. + template iter append_child(iter position, const T &x); + template iter append_child(iter position, T &&x); + template iter prepend_child(iter position, const T &x); + template iter prepend_child(iter position, T &&x); + /// Append the node (plus its children) at other_position as last/first child + /// of position. + template + iter append_child(iter position, iter other_position); + template + iter prepend_child(iter position, iter other_position); + /// Append the nodes in the from-to range (plus their children) as last/first + /// children of position. + template + iter append_children(iter position, sibling_iterator from, + sibling_iterator to); + template + iter prepend_children(iter position, sibling_iterator from, + sibling_iterator to); + + /// Short-hand to insert topmost node in otherwise empty tree. + pre_order_iterator set_head(const T &x); + pre_order_iterator set_head(T &&x); + /// Insert node as previous sibling of node pointed to by position. + template iter insert(iter position, const T &x); + template iter insert(iter position, T &&x); + /// Specialisation of previous member. + sibling_iterator insert(sibling_iterator position, const T &x); + sibling_iterator insert(sibling_iterator position, T &&x); + /// Insert node (with children) pointed to by subtree as previous sibling of + /// node pointed to by position. Does not change the subtree itself (use + /// move_in or move_in_below for that). + template + iter insert_subtree(iter position, const iterator_base &subtree); + /// Insert node as next sibling of node pointed to by position. + template iter insert_after(iter position, const T &x); + template iter insert_after(iter position, T &&x); + /// Insert node (with children) pointed to by subtree as next sibling of node + /// pointed to by position. + template + iter insert_subtree_after(iter position, const iterator_base &subtree); + + /// Replace node at 'position' with other node (keeping same children); + /// 'position' becomes invalid. + template iter replace(iter position, const T &x); + /// Replace node at 'position' with subtree starting at 'from' (do not erase + /// subtree at 'from'); see above. + template + iter replace(iter position, const iterator_base &from); + /// Replace string of siblings (plus their children) with copy of a new string + /// (with children); see above + sibling_iterator replace(sibling_iterator orig_begin, + sibling_iterator orig_end, + sibling_iterator new_begin, + sibling_iterator new_end); + + /// Move all children of node at 'position' to be siblings, returns position. + template iter flatten(iter position); + /// Move nodes in range to be children of 'position'. + template + iter reparent(iter position, sibling_iterator begin, sibling_iterator end); + /// Move all child nodes of 'from' to be children of 'position'. + template iter reparent(iter position, iter from); + + /// Replace node with a new node, making the old node (plus subtree) a child + /// of the new node. + template iter wrap(iter position, const T &x); + /// Replace the range of sibling nodes (plus subtrees), making these children + /// of the new node. + template iter wrap(iter from, iter to, const T &x); + + /// Move 'source' node (plus its children) to become the next sibling of + /// 'target'. + template iter move_after(iter target, iter source); + /// Move 'source' node (plus its children) to become the previous sibling of + /// 'target'. + template iter move_before(iter target, iter source); + sibling_iterator move_before(sibling_iterator target, + sibling_iterator source); + /// Move 'source' node (plus its children) to become the node at 'target' + /// (erasing the node at 'target'). + template iter move_ontop(iter target, iter source); + + /// Extract the subtree starting at the indicated node, removing it from the + /// original tree. + tree move_out(iterator); + /// Inverse of take_out: inserts the given tree as previous sibling of + /// indicated node by a move operation, that is, the given tree becomes empty. + /// Returns iterator to the top node. + template iter move_in(iter, tree &); + /// As above, but now make the tree the last child of the indicated node. + template iter move_in_below(iter, tree &); + /// As above, but now make the tree the nth child of the indicated node (if + /// possible). + template iter move_in_as_nth_child(iter, size_t, tree &); + + /// Merge with other tree, creating new branches and leaves only if they are + /// not already present. + void merge(sibling_iterator, sibling_iterator, sibling_iterator, + sibling_iterator, bool duplicate_leaves = false); + /// As above, but using two trees with a single top node at the 'to' and + /// 'from' positions. + void merge(iterator to, iterator from, bool duplicate_leaves); + /// Sort (std::sort only moves values of nodes, this one moves children as + /// well). + void sort(sibling_iterator from, sibling_iterator to, bool deep = false); + template + void sort(sibling_iterator from, sibling_iterator to, StrictWeakOrdering comp, + bool deep = false); + /// Compare two ranges of nodes (compares nodes as well as tree structure). + template + bool equal(const iter &one, const iter &two, const iter &three) const; + template + bool equal(const iter &one, const iter &two, const iter &three, + BinaryPredicate) const; + template + bool equal_subtree(const iter &one, const iter &two) const; + template + bool equal_subtree(const iter &one, const iter &two, BinaryPredicate) const; + /// Extract a new tree formed by the range of siblings plus all their + /// children. + tree subtree(sibling_iterator from, sibling_iterator to) const; + void subtree(tree &, sibling_iterator from, sibling_iterator to) const; + /// Exchange the node (plus subtree) with its sibling node (do nothing if no + /// sibling present). + void swap(sibling_iterator it); + /// Exchange two nodes (plus subtrees). The iterators will remain valid and + /// keep pointing to the same nodes, which now sit at different locations in + /// the tree. + void swap(iterator, iterator); + + /// Count the total number of nodes. + size_t size() const; + /// Count the total number of nodes below the indicated node (plus one). + size_t size(const iterator_base &) const; + /// Check if tree is empty. + bool empty() const; + /// Compute the depth to the root or to a fixed other iterator. + static int depth(const iterator_base &); + static int depth(const iterator_base &, const iterator_base &); + /// Compute the depth to the root, counting all levels for which predicate + /// returns true. + template + static int depth(const iterator_base &, Predicate p); + /// Compute the depth distance between two nodes, counting all levels for + /// which predicate returns true. + template + static int distance(const iterator_base &top, const iterator_base &bottom, + Predicate p); + /// Determine the maximal depth of the tree. An empty tree has max_depth=-1. + int max_depth() const; + /// Determine the maximal depth of the tree with top node at the given + /// position. + int max_depth(const iterator_base &) const; + /// Count the number of children of node at position. + static unsigned int number_of_children(const iterator_base &); + /// Count the number of siblings (left and right) of node at iterator. Total + /// nodes at this level is +1. + unsigned int number_of_siblings(const iterator_base &) const; + /// Determine whether node at position is in the subtrees with indicated top + /// node. + bool is_in_subtree(const iterator_base &position, + const iterator_base &top) const; + /// Determine whether node at position is in the subtrees with root in the + /// range. + bool is_in_subtree(const iterator_base &position, const iterator_base &begin, + const iterator_base &end) const; + /// Determine whether the iterator is an 'end' iterator and thus not actually + /// pointing to a node. + bool is_valid(const iterator_base &) const; + /// Determine whether the iterator is one of the 'head' nodes at the top + /// level, i.e. has no parent. + static bool is_head(const iterator_base &); + /// Find the lowest common ancestor of two nodes, that is, the deepest node + /// such that both nodes are descendants of it. + iterator lowest_common_ancestor(const iterator_base &, + const iterator_base &) const; + + /// Determine the index of a node in the range of siblings to which it + /// belongs. + unsigned int index(sibling_iterator it) const; + /// Inverse of 'index': return the n-th child of the node at position. + static sibling_iterator child(const iterator_base &position, unsigned int); + /// Return iterator to the sibling indicated by index + sibling_iterator sibling(const iterator_base &position, unsigned int) const; + + /// For debugging only: verify internal consistency by inspecting all pointers + /// in the tree (which will also trigger a valgrind error in case something + /// got corrupted). + void debug_verify_consistency() const; + + /// Comparator class for iterators (compares pointer values; why doesn't this + /// work automatically?) + class iterator_base_less { + public: + bool operator()( + const typename tree::iterator_base &one, + const typename tree::iterator_base &two) const { + return one.node < two.node; + } + }; + tree_node *head, *feet; // head/feet are always dummy; if an iterator points + // to them it is invalid +private: + tree_node_allocator alloc_; + void head_initialise_(); + void copy_(const tree &other); + + /// Comparator class for two nodes of a tree (used for sorting and searching). + template class compare_nodes { + public: + compare_nodes(StrictWeakOrdering comp) : comp_(comp) {} + + bool operator()(const tree_node *a, const tree_node *b) const { + return comp_(a->data, b->data); + } + + private: + StrictWeakOrdering comp_; + }; +}; + +// template +// class iterator_base_less { +// public: +// bool operator()(const typename tree::iterator_base& one, +// const typename tree::iterator_base& two) const +// { +// txtout << "operatorclass<" << one.node < two.node << +// std::endl; return one.node < two.node; +// } +// }; + +// template +// bool operator<(const typename tree::iterator& one, +// const typename tree::iterator& two) +// { +// txtout << "operator< " << one.node < two.node << std::endl; +// if(one.node < two.node) return true; +// return false; +// } +// +// template +// bool operator==(const typename tree::iterator& one, +// const typename tree::iterator& two) +// { +// txtout << "operator== " << one.node == two.node << std::endl; +// if(one.node == two.node) return true; +// return false; +// } +// +// template +// bool operator>(const typename tree::iterator_base& +// one, const typename tree::iterator_base& two) +// { +// txtout << "operator> " << one.node < two.node << std::endl; +// if(one.node > two.node) return true; +// return false; +// } + +// Tree + +template +tree::tree() { + head_initialise_(); +} + +template +tree::tree(const T &x) { + head_initialise_(); + set_head(x); +} + +template +tree::tree(tree &&x) { + head_initialise_(); + if (x.head->next_sibling != x.feet) { // move tree if non-empty only + head->next_sibling = x.head->next_sibling; + feet->prev_sibling = x.feet->prev_sibling; + x.head->next_sibling->prev_sibling = head; + x.feet->prev_sibling->next_sibling = feet; + x.head->next_sibling = x.feet; + x.feet->prev_sibling = x.head; + } +} + +template +tree::tree(const iterator_base &other) { + head_initialise_(); + set_head((*other)); + replace(begin(), other); +} + +template +tree::~tree() { + clear(); + std::allocator_traits::destroy(alloc_, head); + std::allocator_traits::destroy(alloc_, feet); + std::allocator_traits::deallocate(alloc_, head, 1); + std::allocator_traits::deallocate(alloc_, feet, 1); +} + +template +void tree::head_initialise_() { + head = std::allocator_traits::allocate(alloc_, 1, 0); + feet = std::allocator_traits::allocate(alloc_, 1, 0); + std::allocator_traits::construct(alloc_, head, + tree_node_()); + std::allocator_traits::construct(alloc_, feet, + tree_node_()); + + head->parent = 0; + head->first_child = 0; + head->last_child = 0; + head->prev_sibling = 0; // head; + head->next_sibling = feet; // head; + + feet->parent = 0; + feet->first_child = 0; + feet->last_child = 0; + feet->prev_sibling = head; + feet->next_sibling = 0; +} + +template +tree &tree::operator=( + const tree &other) { + if (this != &other) + copy_(other); + return *this; +} + +template +tree & +tree::operator=(tree &&x) { + if (this != &x) { + clear(); // clear any existing data. + + head->next_sibling = x.head->next_sibling; + feet->prev_sibling = x.feet->prev_sibling; + x.head->next_sibling->prev_sibling = head; + x.feet->prev_sibling->next_sibling = feet; + x.head->next_sibling = x.feet; + x.feet->prev_sibling = x.head; + } + return *this; +} + +template +tree::tree(const tree &other) { + head_initialise_(); + copy_(other); +} + +template +void tree::copy_( + const tree &other) { + clear(); + pre_order_iterator it = other.begin(), to = begin(); + while (it != other.end()) { + to = insert(to, (*it)); + it.skip_children(); + ++it; + } + to = begin(); + it = other.begin(); + while (it != other.end()) { + to = replace(to, it); + to.skip_children(); + it.skip_children(); + ++to; + ++it; + } +} + +template +void tree::clear() { + if (head) + while (head->next_sibling != feet) + erase(pre_order_iterator(head->next_sibling)); +} + +template +void tree::erase_children(const iterator_base &it) { + // std::cout << "erase_children " << it.node << std::endl; + if (it.node == 0) + return; + + tree_node *cur = it.node->first_child; + tree_node *prev = 0; + + while (cur != 0) { + prev = cur; + cur = cur->next_sibling; + erase_children(pre_order_iterator(prev)); + std::allocator_traits::destroy(alloc_, prev); + std::allocator_traits::deallocate(alloc_, prev, 1); + } + it.node->first_child = 0; + it.node->last_child = 0; + // std::cout << "exit" << std::endl; +} + +template +void tree::erase_right_siblings( + const iterator_base &it) { + if (it.node == 0) + return; + + tree_node *cur = it.node->next_sibling; + tree_node *prev = 0; + + while (cur != 0) { + prev = cur; + cur = cur->next_sibling; + erase_children(pre_order_iterator(prev)); + std::allocator_traits::destroy(alloc_, prev); + std::allocator_traits::deallocate(alloc_, prev, 1); + } + it.node->next_sibling = 0; + if (it.node->parent != 0) + it.node->parent->last_child = it.node; +} + +template +void tree::erase_left_siblings( + const iterator_base &it) { + if (it.node == 0) + return; + + tree_node *cur = it.node->prev_sibling; + tree_node *prev = 0; + + while (cur != 0) { + prev = cur; + cur = cur->prev_sibling; + erase_children(pre_order_iterator(prev)); + std::allocator_traits::destroy(alloc_, prev); + std::allocator_traits::deallocate(alloc_, prev, 1); + } + it.node->prev_sibling = 0; + if (it.node->parent != 0) + it.node->parent->first_child = it.node; +} + +template +template +iter tree::erase(iter it) { + tree_node *cur = it.node; + assert(cur != head); + iter ret = it; + ret.skip_children(); + ++ret; + erase_children(it); + if (cur->prev_sibling == 0) { + cur->parent->first_child = cur->next_sibling; + } else { + cur->prev_sibling->next_sibling = cur->next_sibling; + } + if (cur->next_sibling == 0) { + cur->parent->last_child = cur->prev_sibling; + } else { + cur->next_sibling->prev_sibling = cur->prev_sibling; + } + + std::allocator_traits::destroy(alloc_, cur); + std::allocator_traits::deallocate(alloc_, cur, 1); + return ret; +} + +template +typename tree::pre_order_iterator +tree::begin() const { + return pre_order_iterator(head->next_sibling); +} + +template +typename tree::pre_order_iterator +tree::end() const { + return pre_order_iterator(feet); +} + +template +typename tree::breadth_first_queued_iterator +tree::begin_breadth_first() const { + return breadth_first_queued_iterator(head->next_sibling); +} + +template +typename tree::breadth_first_queued_iterator +tree::end_breadth_first() const { + return breadth_first_queued_iterator(); +} + +template +typename tree::post_order_iterator +tree::begin_post() const { + tree_node *tmp = head->next_sibling; + if (tmp != feet) { + while (tmp->first_child) + tmp = tmp->first_child; + } + return post_order_iterator(tmp); +} + +template +typename tree::post_order_iterator +tree::end_post() const { + return post_order_iterator(feet); +} + +template +typename tree::fixed_depth_iterator +tree::begin_fixed(const iterator_base &pos, + unsigned int dp, + bool walk_back) const { + typename tree::fixed_depth_iterator ret; + ret.top_node = pos.node; + + tree_node *tmp = pos.node; + unsigned int curdepth = 0; + while (curdepth < dp) { // go down one level + while (tmp->first_child == 0) { + if (tmp->next_sibling == 0) { + // try to walk up and then right again + do { + if (tmp == ret.top_node) + throw std::range_error("tree: begin_fixed out of range"); + tmp = tmp->parent; + if (tmp == 0) + throw std::range_error("tree: begin_fixed out of range"); + --curdepth; + } while (tmp->next_sibling == 0); + } + tmp = tmp->next_sibling; + } + tmp = tmp->first_child; + ++curdepth; + } + + // Now walk back to the first sibling in this range. + if (walk_back) + while (tmp->prev_sibling != 0) + tmp = tmp->prev_sibling; + + ret.node = tmp; + return ret; +} + +template +typename tree::fixed_depth_iterator +tree::end_fixed(const iterator_base &pos, + unsigned int dp) const { + assert(1 == + 0); // FIXME: not correct yet: use is_valid() as a temporary workaround + tree_node *tmp = pos.node; + unsigned int curdepth = 1; + while (curdepth < dp) { // go down one level + while (tmp->first_child == 0) { + tmp = tmp->next_sibling; + if (tmp == 0) + throw std::range_error("tree: end_fixed out of range"); + } + tmp = tmp->first_child; + ++curdepth; + } + return tmp; +} + +template +typename tree::sibling_iterator +tree::begin(const iterator_base &pos) { + assert(pos.node != 0); + if (pos.node->first_child == 0) { + return end(pos); + } + return pos.node->first_child; +} + +template +typename tree::sibling_iterator +tree::end(const iterator_base &pos) { + sibling_iterator ret(0); + ret.parent_ = pos.node; + return ret; +} + +template +typename tree::leaf_iterator +tree::begin_leaf() const { + tree_node *tmp = head->next_sibling; + if (tmp != feet) { + while (tmp->first_child) + tmp = tmp->first_child; + } + return leaf_iterator(tmp); +} + +template +typename tree::leaf_iterator +tree::end_leaf() const { + return leaf_iterator(feet); +} + +template +typename tree::path_t +tree::path_from_iterator( + const iterator_base &iter, const iterator_base &top) const { + path_t path; + tree_node *walk = iter.node; + + do { + if (path.size() > 0) + walk = walk->parent; + int num = 0; + while (walk != top.node && walk->prev_sibling != 0 && + walk->prev_sibling != head) { + ++num; + walk = walk->prev_sibling; + } + path.push_back(num); + } while (walk->parent != 0 && walk != top.node); + + std::reverse(path.begin(), path.end()); + return path; +} + +template +typename tree::iterator +tree::iterator_from_path( + const path_t &path, const iterator_base &top) const { + iterator it = top; + tree_node *walk = it.node; + + for (size_t step = 0; step < path.size(); ++step) { + if (step > 0) + walk = walk->first_child; + if (walk == 0) + throw std::range_error( + "tree::iterator_from_path: no more nodes at step " + + std::to_string(step)); + + for (int i = 0; i < path[step]; ++i) { + walk = walk->next_sibling; + if (walk == 0) + throw std::range_error( + "tree::iterator_from_path: out of siblings at step " + + std::to_string(step)); + } + } + it.node = walk; + return it; +} + +template +typename tree::leaf_iterator +tree::begin_leaf(const iterator_base &top) const { + tree_node *tmp = top.node; + while (tmp->first_child) + tmp = tmp->first_child; + return leaf_iterator(tmp, top.node); +} + +template +typename tree::leaf_iterator +tree::end_leaf(const iterator_base &top) const { + return leaf_iterator(top.node, top.node); +} + +template +template +iter tree::parent(iter position) { + if (position.node == 0) + throw navigation_error("tree: attempt to navigate from null iterator."); + + if (position.node->parent == 0) + throw navigation_error("tree: attempt to navigate up past head node."); + + return iter(position.node->parent); +} + +template +template +iter tree::previous_sibling(iter position) { + assert(position.node != 0); + iter ret(position); + ret.node = position.node->prev_sibling; + return ret; +} + +template +template +iter tree::next_sibling(iter position) { + assert(position.node != 0); + iter ret(position); + ret.node = position.node->next_sibling; + return ret; +} + +template +template +iter tree::next_at_same_depth(iter position) const { + // We make use of a temporary fixed_depth iterator to implement this. + + typename tree::fixed_depth_iterator tmp( + position.node); + + ++tmp; + return iter(tmp); + + // assert(position.node!=0); + // iter ret(position); + // + // if(position.node->next_sibling) { + // ret.node=position.node->next_sibling; + // } + // else { + // int relative_depth=0; + // upper: + // do { + // ret.node=ret.node->parent; + // if(ret.node==0) return ret; + // --relative_depth; + // } while(ret.node->next_sibling==0); + // lower: + // ret.node=ret.node->next_sibling; + // while(ret.node->first_child==0) { + // if(ret.node->next_sibling==0) + // goto upper; + // ret.node=ret.node->next_sibling; + // if(ret.node==0) return ret; + // } + // while(relative_depth<0 && ret.node->first_child!=0) { + // ret.node=ret.node->first_child; + // ++relative_depth; + // } + // if(relative_depth<0) { + // if(ret.node->next_sibling==0) goto upper; + // else goto lower; + // } + // } + // return ret; +} + +template +template +iter tree::append_child(iter position) { + assert(position.node != head); + assert(position.node != feet); + assert(position.node); + + tree_node *tmp = + std::allocator_traits::allocate(alloc_, 1, 0); + std::allocator_traits::construct(alloc_, tmp, + tree_node_()); + tmp->first_child = 0; + tmp->last_child = 0; + + tmp->parent = position.node; + if (position.node->last_child != 0) { + position.node->last_child->next_sibling = tmp; + } else { + position.node->first_child = tmp; + } + tmp->prev_sibling = position.node->last_child; + position.node->last_child = tmp; + tmp->next_sibling = 0; + return tmp; +} + +template +template +iter tree::prepend_child(iter position) { + assert(position.node != head); + assert(position.node != feet); + assert(position.node); + + tree_node *tmp = + std::allocator_traits::allocate(alloc_, 1, 0); + std::allocator_traits::construct(alloc_, tmp, + tree_node_()); + tmp->first_child = 0; + tmp->last_child = 0; + + tmp->parent = position.node; + if (position.node->first_child != 0) { + position.node->first_child->prev_sibling = tmp; + } else { + position.node->last_child = tmp; + } + tmp->next_sibling = position.node->first_child; + position.node->prev_child = tmp; + tmp->prev_sibling = 0; + return tmp; +} + +template +template +iter tree::append_child(iter position, const T &x) { + // If your program fails here you probably used 'append_child' to add the top + // node to an empty tree. From version 1.45 the top element should be added + // using 'insert'. See the documentation for further information, and sorry + // about the API change. + assert(position.node != head); + assert(position.node != feet); + assert(position.node); + + tree_node *tmp = + std::allocator_traits::allocate(alloc_, 1, 0); + std::allocator_traits::construct(alloc_, tmp, x); + tmp->first_child = 0; + tmp->last_child = 0; + + tmp->parent = position.node; + if (position.node->last_child != 0) { + position.node->last_child->next_sibling = tmp; + } else { + position.node->first_child = tmp; + } + tmp->prev_sibling = position.node->last_child; + position.node->last_child = tmp; + tmp->next_sibling = 0; + return tmp; +} + +template +template +iter tree::append_child(iter position, T &&x) { + assert(position.node != head); + assert(position.node != feet); + assert(position.node); + + tree_node *tmp = + std::allocator_traits::allocate(alloc_, 1, 0); + std::allocator_traits::construct( + alloc_, tmp); // Here is where the move semantics kick in + std::swap(tmp->data, x); + + tmp->first_child = 0; + tmp->last_child = 0; + + tmp->parent = position.node; + if (position.node->last_child != 0) { + position.node->last_child->next_sibling = tmp; + } else { + position.node->first_child = tmp; + } + tmp->prev_sibling = position.node->last_child; + position.node->last_child = tmp; + tmp->next_sibling = 0; + return tmp; +} + +template +template +iter tree::prepend_child(iter position, const T &x) { + assert(position.node != head); + assert(position.node != feet); + assert(position.node); + + tree_node *tmp = + std::allocator_traits::allocate(alloc_, 1, 0); + std::allocator_traits::construct(alloc_, tmp, x); + tmp->first_child = 0; + tmp->last_child = 0; + + tmp->parent = position.node; + if (position.node->first_child != 0) { + position.node->first_child->prev_sibling = tmp; + } else { + position.node->last_child = tmp; + } + tmp->next_sibling = position.node->first_child; + position.node->first_child = tmp; + tmp->prev_sibling = 0; + return tmp; +} + +template +template +iter tree::prepend_child(iter position, T &&x) { + assert(position.node != head); + assert(position.node != feet); + assert(position.node); + + tree_node *tmp = + std::allocator_traits::allocate(alloc_, 1, 0); + std::allocator_traits::construct(alloc_, tmp); + std::swap(tmp->data, x); + + tmp->first_child = 0; + tmp->last_child = 0; + + tmp->parent = position.node; + if (position.node->first_child != 0) { + position.node->first_child->prev_sibling = tmp; + } else { + position.node->last_child = tmp; + } + tmp->next_sibling = position.node->first_child; + position.node->first_child = tmp; + tmp->prev_sibling = 0; + return tmp; +} + +template +template +iter tree::append_child(iter position, iter other) { + assert(position.node != head); + assert(position.node != feet); + assert(position.node); + + sibling_iterator aargh = append_child(position, value_type()); + return replace(aargh, other); +} + +template +template +iter tree::prepend_child(iter position, iter other) { + assert(position.node != head); + assert(position.node != feet); + assert(position.node); + + sibling_iterator aargh = prepend_child(position, value_type()); + return replace(aargh, other); +} + +template +template +iter tree::append_children(iter position, + sibling_iterator from, + sibling_iterator to) { + assert(position.node != head); + assert(position.node != feet); + assert(position.node); + + iter ret = from; + + while (from != to) { + insert_subtree(position.end(), from); + ++from; + } + return ret; +} + +template +template +iter tree::prepend_children(iter position, + sibling_iterator from, + sibling_iterator to) { + assert(position.node != head); + assert(position.node != feet); + assert(position.node); + + if (from == to) + return from; // should return end of tree? + + iter ret; + do { + --to; + ret = insert_subtree(position.begin(), to); + } while (to != from); + + return ret; +} + +template +typename tree::pre_order_iterator +tree::set_head(const T &x) { + assert(head->next_sibling == feet); + return insert(iterator(feet), x); +} + +template +typename tree::pre_order_iterator +tree::set_head(T &&x) { + assert(head->next_sibling == feet); + return insert(iterator(feet), x); +} + +template +template +iter tree::insert(iter position, const T &x) { + if (position.node == 0) { + position.node = feet; // Backward compatibility: when calling insert on a + // null node, insert before the feet. + } + assert(position.node != head); // Cannot insert before head. + + tree_node *tmp = + std::allocator_traits::allocate(alloc_, 1, 0); + std::allocator_traits::construct(alloc_, tmp, x); + tmp->first_child = 0; + tmp->last_child = 0; + + tmp->parent = position.node->parent; + tmp->next_sibling = position.node; + tmp->prev_sibling = position.node->prev_sibling; + position.node->prev_sibling = tmp; + + if (tmp->prev_sibling == 0) { + if (tmp->parent) // when inserting nodes at the head, there is no parent + tmp->parent->first_child = tmp; + } else + tmp->prev_sibling->next_sibling = tmp; + return tmp; +} + +template +template +iter tree::insert(iter position, T &&x) { + if (position.node == 0) { + position.node = feet; // Backward compatibility: when calling insert on a + // null node, insert before the feet. + } + tree_node *tmp = + std::allocator_traits::allocate(alloc_, 1, 0); + std::allocator_traits::construct(alloc_, tmp); + std::swap(tmp->data, x); // Move semantics + tmp->first_child = 0; + tmp->last_child = 0; + + tmp->parent = position.node->parent; + tmp->next_sibling = position.node; + tmp->prev_sibling = position.node->prev_sibling; + position.node->prev_sibling = tmp; + + if (tmp->prev_sibling == 0) { + if (tmp->parent) // when inserting nodes at the head, there is no parent + tmp->parent->first_child = tmp; + } else + tmp->prev_sibling->next_sibling = tmp; + return tmp; +} + +template +typename tree::sibling_iterator +tree::insert(sibling_iterator position, const T &x) { + tree_node *tmp = + std::allocator_traits::allocate(alloc_, 1, 0); + std::allocator_traits::construct(alloc_, tmp, x); + tmp->first_child = 0; + tmp->last_child = 0; + + tmp->next_sibling = position.node; + if (position.node == 0) { // iterator points to end of a subtree + tmp->parent = position.parent_; + tmp->prev_sibling = position.range_last(); + tmp->parent->last_child = tmp; + } else { + tmp->parent = position.node->parent; + tmp->prev_sibling = position.node->prev_sibling; + position.node->prev_sibling = tmp; + } + + if (tmp->prev_sibling == 0) { + if (tmp->parent) // when inserting nodes at the head, there is no parent + tmp->parent->first_child = tmp; + } else + tmp->prev_sibling->next_sibling = tmp; + return tmp; +} + +template +typename tree::sibling_iterator +tree::insert(sibling_iterator position, T &&x) { + tree_node *tmp = + std::allocator_traits::allocate(alloc_, 1, 0); + std::allocator_traits::construct(alloc_, tmp); + std::swap(tmp->data, x); // Move semantics + + tmp->first_child = 0; + tmp->last_child = 0; + + tmp->next_sibling = position.node; + if (position.node == 0) { // iterator points to end of a subtree + tmp->parent = position.parent_; + tmp->prev_sibling = position.range_last(); + tmp->parent->last_child = tmp; + } else { + tmp->parent = position.node->parent; + tmp->prev_sibling = position.node->prev_sibling; + position.node->prev_sibling = tmp; + } + + if (tmp->prev_sibling == 0) { + if (tmp->parent) // when inserting nodes at the head, there is no parent + tmp->parent->first_child = tmp; + } else + tmp->prev_sibling->next_sibling = tmp; + return tmp; +} + +template +template +iter tree::insert_after(iter position, const T &x) { + tree_node *tmp = + std::allocator_traits::allocate(alloc_, 1, 0); + std::allocator_traits::construct(alloc_, tmp, x); + tmp->first_child = 0; + tmp->last_child = 0; + + tmp->parent = position.node->parent; + tmp->prev_sibling = position.node; + tmp->next_sibling = position.node->next_sibling; + position.node->next_sibling = tmp; + + if (tmp->next_sibling == 0) { + if (tmp->parent) // when inserting nodes at the head, there is no parent + tmp->parent->last_child = tmp; + } else { + tmp->next_sibling->prev_sibling = tmp; + } + return tmp; +} + +template +template +iter tree::insert_after(iter position, T &&x) { + tree_node *tmp = + std::allocator_traits::allocate(alloc_, 1, 0); + std::allocator_traits::construct(alloc_, tmp); + std::swap(tmp->data, x); // move semantics + tmp->first_child = 0; + tmp->last_child = 0; + + tmp->parent = position.node->parent; + tmp->prev_sibling = position.node; + tmp->next_sibling = position.node->next_sibling; + position.node->next_sibling = tmp; + + if (tmp->next_sibling == 0) { + if (tmp->parent) // when inserting nodes at the head, there is no parent + tmp->parent->last_child = tmp; + } else { + tmp->next_sibling->prev_sibling = tmp; + } + return tmp; +} + +template +template +iter tree::insert_subtree( + iter position, const iterator_base &subtree) { + // insert dummy + iter it = insert(position, value_type()); + // replace dummy with subtree + return replace(it, subtree); +} + +template +template +iter tree::insert_subtree_after( + iter position, const iterator_base &subtree) { + // insert dummy + iter it = insert_after(position, value_type()); + // replace dummy with subtree + return replace(it, subtree); +} + +// template +// template +// iter tree::insert_subtree(sibling_iterator position, +// iter subtree) +// { +// // insert dummy +// iter it(insert(position, value_type())); +// // replace dummy with subtree +// return replace(it, subtree); +// } + +template +template +iter tree::replace(iter position, const T &x) { + // kp::destructor(&position.node->data); + // kp::constructor(&position.node->data, x); + position.node->data = x; + // alloc_.destroy(position.node); + // alloc_.construct(position.node, x); + return position; +} + +template +template +iter tree::replace(iter position, + const iterator_base &from) { + assert(position.node != head); + tree_node *current_from = from.node; + tree_node *start_from = from.node; + tree_node *current_to = position.node; + + // replace the node at position with head of the replacement tree at from + // std::cout << "warning!" << position.node << std::endl; + erase_children(position); + // std::cout << "no warning!" << std::endl; + tree_node *tmp = + std::allocator_traits::allocate(alloc_, 1, 0); + std::allocator_traits::construct(alloc_, tmp, (*from)); + tmp->first_child = 0; + tmp->last_child = 0; + if (current_to->prev_sibling == 0) { + if (current_to->parent != 0) + current_to->parent->first_child = tmp; + } else { + current_to->prev_sibling->next_sibling = tmp; + } + tmp->prev_sibling = current_to->prev_sibling; + if (current_to->next_sibling == 0) { + if (current_to->parent != 0) + current_to->parent->last_child = tmp; + } else { + current_to->next_sibling->prev_sibling = tmp; + } + tmp->next_sibling = current_to->next_sibling; + tmp->parent = current_to->parent; + // kp::destructor(¤t_to->data); + std::allocator_traits::destroy(alloc_, current_to); + std::allocator_traits::deallocate(alloc_, current_to, 1); + current_to = tmp; + + // only at this stage can we fix 'last' + tree_node *last = from.node->next_sibling; + + pre_order_iterator toit = tmp; + // copy all children + do { + assert(current_from != 0); + if (current_from->first_child != 0) { + current_from = current_from->first_child; + toit = append_child(toit, current_from->data); + } else { + while (current_from->next_sibling == 0 && current_from != start_from) { + current_from = current_from->parent; + toit = parent(toit); + assert(current_from != 0); + } + current_from = current_from->next_sibling; + if (current_from != last) { + toit = append_child(parent(toit), current_from->data); + } + } + } while (current_from != last); + + return current_to; +} + +template +typename tree::sibling_iterator +tree::replace(sibling_iterator orig_begin, + sibling_iterator orig_end, + sibling_iterator new_begin, + sibling_iterator new_end) { + tree_node *orig_first = orig_begin.node; + tree_node *new_first = new_begin.node; + tree_node *orig_last = orig_first; + while ((++orig_begin) != orig_end) + orig_last = orig_last->next_sibling; + tree_node *new_last = new_first; + while ((++new_begin) != new_end) + new_last = new_last->next_sibling; + + // insert all siblings in new_first..new_last before orig_first + bool first = true; + pre_order_iterator ret; + while (1 == 1) { + pre_order_iterator tt = insert_subtree(pre_order_iterator(orig_first), + pre_order_iterator(new_first)); + if (first) { + ret = tt; + first = false; + } + if (new_first == new_last) + break; + new_first = new_first->next_sibling; + } + + // erase old range of siblings + bool last = false; + tree_node *next = orig_first; + while (1 == 1) { + if (next == orig_last) + last = true; + next = next->next_sibling; + erase((pre_order_iterator)orig_first); + if (last) + break; + orig_first = next; + } + return ret; +} + +template +template +iter tree::flatten(iter position) { + if (position.node->first_child == 0) + return position; + + tree_node *tmp = position.node->first_child; + while (tmp) { + tmp->parent = position.node->parent; + tmp = tmp->next_sibling; + } + if (position.node->next_sibling) { + position.node->last_child->next_sibling = position.node->next_sibling; + position.node->next_sibling->prev_sibling = position.node->last_child; + } else { + position.node->parent->last_child = position.node->last_child; + } + position.node->next_sibling = position.node->first_child; + position.node->next_sibling->prev_sibling = position.node; + position.node->first_child = 0; + position.node->last_child = 0; + + return position; +} + +template +template +iter tree::reparent(iter position, + sibling_iterator begin, + sibling_iterator end) { + tree_node *first = begin.node; + tree_node *last = first; + + assert(first != position.node); + + if (begin == end) + return begin; + // determine last node + while ((++begin) != end) { + last = last->next_sibling; + } + // move subtree + if (first->prev_sibling == 0) { + first->parent->first_child = last->next_sibling; + } else { + first->prev_sibling->next_sibling = last->next_sibling; + } + if (last->next_sibling == 0) { + last->parent->last_child = first->prev_sibling; + } else { + last->next_sibling->prev_sibling = first->prev_sibling; + } + if (position.node->first_child == 0) { + position.node->first_child = first; + position.node->last_child = last; + first->prev_sibling = 0; + } else { + position.node->last_child->next_sibling = first; + first->prev_sibling = position.node->last_child; + position.node->last_child = last; + } + last->next_sibling = 0; + + tree_node *pos = first; + for (;;) { + pos->parent = position.node; + if (pos == last) + break; + pos = pos->next_sibling; + } + + return first; +} + +template +template +iter tree::reparent(iter position, iter from) { + if (from.node->first_child == 0) + return position; + return reparent(position, from.node->first_child, end(from)); +} + +template +template +iter tree::wrap(iter position, const T &x) { + assert(position.node != 0); + sibling_iterator fr = position, to = position; + ++to; + iter ret = insert(position, x); + reparent(ret, fr, to); + return ret; +} + +template +template +iter tree::wrap(iter from, iter to, const T &x) { + assert(from.node != 0); + iter ret = insert(from, x); + reparent(ret, from, to); + return ret; +} + +template +template +iter tree::move_after(iter target, iter source) { + tree_node *dst = target.node; + tree_node *src = source.node; + assert(dst); + assert(src); + + if (dst == src) + return source; + if (dst->next_sibling) + if (dst->next_sibling == src) // already in the right spot + return source; + + // take src out of the tree + if (src->prev_sibling != 0) + src->prev_sibling->next_sibling = src->next_sibling; + else + src->parent->first_child = src->next_sibling; + if (src->next_sibling != 0) + src->next_sibling->prev_sibling = src->prev_sibling; + else + src->parent->last_child = src->prev_sibling; + + // connect it to the new point + if (dst->next_sibling != 0) + dst->next_sibling->prev_sibling = src; + else + dst->parent->last_child = src; + src->next_sibling = dst->next_sibling; + dst->next_sibling = src; + src->prev_sibling = dst; + src->parent = dst->parent; + return src; +} + +template +template +iter tree::move_before(iter target, iter source) { + tree_node *dst = target.node; + tree_node *src = source.node; + assert(dst); + assert(src); + + if (dst == src) + return source; + if (dst->prev_sibling) + if (dst->prev_sibling == src) // already in the right spot + return source; + + // take src out of the tree + if (src->prev_sibling != 0) + src->prev_sibling->next_sibling = src->next_sibling; + else + src->parent->first_child = src->next_sibling; + if (src->next_sibling != 0) + src->next_sibling->prev_sibling = src->prev_sibling; + else + src->parent->last_child = src->prev_sibling; + + // connect it to the new point + if (dst->prev_sibling != 0) + dst->prev_sibling->next_sibling = src; + else + dst->parent->first_child = src; + src->prev_sibling = dst->prev_sibling; + dst->prev_sibling = src; + src->next_sibling = dst; + src->parent = dst->parent; + return src; +} + +// specialisation for sibling_iterators +template +typename tree::sibling_iterator +tree::move_before(sibling_iterator target, + sibling_iterator source) { + tree_node *dst = target.node; + tree_node *src = source.node; + tree_node *dst_prev_sibling; + if (dst == 0) { // must then be an end iterator + dst_prev_sibling = target.parent_->last_child; + assert(dst_prev_sibling); + } else + dst_prev_sibling = dst->prev_sibling; + assert(src); + + if (dst == src) + return source; + if (dst_prev_sibling) + if (dst_prev_sibling == src) // already in the right spot + return source; + + // take src out of the tree + if (src->prev_sibling != 0) + src->prev_sibling->next_sibling = src->next_sibling; + else + src->parent->first_child = src->next_sibling; + if (src->next_sibling != 0) + src->next_sibling->prev_sibling = src->prev_sibling; + else + src->parent->last_child = src->prev_sibling; + + // connect it to the new point + if (dst_prev_sibling != 0) + dst_prev_sibling->next_sibling = src; + else + target.parent_->first_child = src; + src->prev_sibling = dst_prev_sibling; + if (dst) { + dst->prev_sibling = src; + src->parent = dst->parent; + } else { + src->parent = dst_prev_sibling->parent; + } + src->next_sibling = dst; + return src; +} + +template +template +iter tree::move_ontop(iter target, iter source) { + tree_node *dst = target.node; + tree_node *src = source.node; + assert(dst); + assert(src); + + if (dst == src) + return source; + + // if(dst==src->prev_sibling) { + // + // } + + // remember connection points + tree_node *b_prev_sibling = dst->prev_sibling; + tree_node *b_next_sibling = dst->next_sibling; + tree_node *b_parent = dst->parent; + + // remove target + erase(target); + + // take src out of the tree + if (src->prev_sibling != 0) + src->prev_sibling->next_sibling = src->next_sibling; + else { + assert(src->parent != 0); + src->parent->first_child = src->next_sibling; + } + if (src->next_sibling != 0) + src->next_sibling->prev_sibling = src->prev_sibling; + else { + assert(src->parent != 0); + src->parent->last_child = src->prev_sibling; + } + + // connect it to the new point + if (b_prev_sibling != 0) + b_prev_sibling->next_sibling = src; + else { + assert(b_parent != 0); + b_parent->first_child = src; + } + if (b_next_sibling != 0) + b_next_sibling->prev_sibling = src; + else { + assert(b_parent != 0); + b_parent->last_child = src; + } + src->prev_sibling = b_prev_sibling; + src->next_sibling = b_next_sibling; + src->parent = b_parent; + return src; +} + +template +tree +tree::move_out(iterator source) { + tree ret; + + // Move source node into the 'ret' tree. + ret.head->next_sibling = source.node; + ret.feet->prev_sibling = source.node; + + // Close the links in the current tree. + if (source.node->prev_sibling != 0) + source.node->prev_sibling->next_sibling = source.node->next_sibling; + + if (source.node->next_sibling != 0) + source.node->next_sibling->prev_sibling = source.node->prev_sibling; + + // If the moved-out node was a first or last child of + // the parent, adjust those links. + if (source.node->parent->first_child == source.node) { + if (source.node->next_sibling != 0) + source.node->parent->first_child = source.node->next_sibling; + else + source.node->parent->first_child = 0; + } + if (source.node->parent->last_child == source.node) { + if (source.node->prev_sibling != 0) + source.node->parent->last_child = source.node->prev_sibling; + else + source.node->parent->last_child = 0; + } + source.node->parent = 0; + + // Fix source prev/next links. + source.node->prev_sibling = ret.head; + source.node->next_sibling = ret.feet; + + return ret; // A good compiler will move this, not copy. +} + +template +template +iter tree::move_in(iter loc, tree &other) { + if (other.head->next_sibling == other.feet) + return loc; // other tree is empty + + tree_node *other_first_head = other.head->next_sibling; + tree_node *other_last_head = other.feet->prev_sibling; + + sibling_iterator prev(loc); + --prev; + + prev.node->next_sibling = other_first_head; + loc.node->prev_sibling = other_last_head; + other_first_head->prev_sibling = prev.node; + other_last_head->next_sibling = loc.node; + + // Adjust parent pointers. + tree_node *walk = other_first_head; + while (true) { + walk->parent = loc.node->parent; + if (walk == other_last_head) + break; + walk = walk->next_sibling; + } + + // Close other tree. + other.head->next_sibling = other.feet; + other.feet->prev_sibling = other.head; + + return other_first_head; +} + +template +template +iter tree::move_in_below(iter loc, tree &other) { + if (other.head->next_sibling == other.feet) + return loc; // other tree is empty + + auto n = other.number_of_children(loc); + return move_in_as_nth_child(loc, n, other); +} + +template +template +iter tree::move_in_as_nth_child(iter loc, size_t n, + tree &other) { + if (other.head->next_sibling == other.feet) + return loc; // other tree is empty + + tree_node *other_first_head = other.head->next_sibling; + tree_node *other_last_head = other.feet->prev_sibling; + + if (n == 0) { + if (loc.node->first_child == 0) { + loc.node->first_child = other_first_head; + loc.node->last_child = other_last_head; + other_last_head->next_sibling = 0; + other_first_head->prev_sibling = 0; + } else { + loc.node->first_child->prev_sibling = other_last_head; + other_last_head->next_sibling = loc.node->first_child; + loc.node->first_child = other_first_head; + other_first_head->prev_sibling = 0; + } + } else { + --n; + tree_node *walk = loc.node->first_child; + while (true) { + if (walk == 0) + throw std::range_error( + "tree: move_in_as_nth_child position out of range"); + if (n == 0) + break; + --n; + walk = walk->next_sibling; + } + if (walk->next_sibling == 0) + loc.node->last_child = other_last_head; + else + walk->next_sibling->prev_sibling = other_last_head; + other_last_head->next_sibling = walk->next_sibling; + walk->next_sibling = other_first_head; + other_first_head->prev_sibling = walk; + } + + // Adjust parent pointers. + tree_node *walk = other_first_head; + while (true) { + walk->parent = loc.node; + if (walk == other_last_head) + break; + walk = walk->next_sibling; + } + + // Close other tree. + other.head->next_sibling = other.feet; + other.feet->prev_sibling = other.head; + + return other_first_head; +} + +template +void tree::merge(sibling_iterator to1, + sibling_iterator to2, + sibling_iterator from1, + sibling_iterator from2, + bool duplicate_leaves) { + sibling_iterator fnd; + while (from1 != from2) { + if ((fnd = std::find(to1, to2, (*from1))) != to2) { // element found + if (from1.begin() == from1.end()) { // full depth reached + if (duplicate_leaves) + append_child(parent(to1), (*from1)); + } else { // descend further + merge(fnd.begin(), fnd.end(), from1.begin(), from1.end(), + duplicate_leaves); + } + } else { // element missing + insert_subtree(to2, from1); + } + ++from1; + } +} + +template +void tree::merge(iterator to, iterator from, + bool duplicate_leaves) { + sibling_iterator to1(to); + sibling_iterator to2 = to1; + ++to2; + sibling_iterator from1(from); + sibling_iterator from2 = from1; + ++from2; + + merge(to1, to2, from1, from2, duplicate_leaves); +} + +template +void tree::sort(sibling_iterator from, + sibling_iterator to, bool deep) { + std::less comp; + sort(from, to, comp, deep); +} + +template +template +void tree::sort(sibling_iterator from, + sibling_iterator to, + StrictWeakOrdering comp, bool deep) { + if (from == to) + return; + // make list of sorted nodes + // CHECK: if multiset stores equivalent nodes in the order in which they + // are inserted, then this routine should be called 'stable_sort'. + std::multiset> nodes(comp); + sibling_iterator it = from, it2 = to; + while (it != to) { + nodes.insert(it.node); + ++it; + } + // reassemble + --it2; + + // prev and next are the nodes before and after the sorted range + tree_node *prev = from.node->prev_sibling; + tree_node *next = it2.node->next_sibling; + typename std::multiset>::iterator + nit = nodes.begin(), + eit = nodes.end(); + if (prev == 0) { + if ((*nit)->parent != + 0) // to catch "sorting the head" situations, when there is no parent + (*nit)->parent->first_child = (*nit); + } else + prev->next_sibling = (*nit); + + --eit; + while (nit != eit) { + (*nit)->prev_sibling = prev; + if (prev) + prev->next_sibling = (*nit); + prev = (*nit); + ++nit; + } + // prev now points to the last-but-one node in the sorted range + if (prev) + prev->next_sibling = (*eit); + + // eit points to the last node in the sorted range. + (*eit)->next_sibling = next; + (*eit)->prev_sibling = prev; // missed in the loop above + if (next == 0) { + if ((*eit)->parent != + 0) // to catch "sorting the head" situations, when there is no parent + (*eit)->parent->last_child = (*eit); + } else + next->prev_sibling = (*eit); + + if (deep) { // sort the children of each node too + sibling_iterator bcs(*nodes.begin()); + sibling_iterator ecs(*eit); + ++ecs; + while (bcs != ecs) { + sort(begin(bcs), end(bcs), comp, deep); + ++bcs; + } + } +} + +template +template +bool tree::equal(const iter &one_, const iter &two, + const iter &three_) const { + std::equal_to comp; + return equal(one_, two, three_, comp); +} + +template +template +bool tree::equal_subtree(const iter &one_, + const iter &two_) const { + std::equal_to comp; + return equal_subtree(one_, two_, comp); +} + +template +template +bool tree::equal(const iter &one_, const iter &two, + const iter &three_, + BinaryPredicate fun) const { + pre_order_iterator one(one_), three(three_); + + // if(one==two && is_valid(three) && three.number_of_children()!=0) + // return false; + while (one != two && is_valid(three)) { + if (!fun(*one, *three)) + return false; + if (one.number_of_children() != three.number_of_children()) + return false; + ++one; + ++three; + } + return true; +} + +template +template +bool tree::equal_subtree(const iter &one_, + const iter &two_, + BinaryPredicate fun) const { + pre_order_iterator one(one_), two(two_); + + if (!fun(*one, *two)) + return false; + if (number_of_children(one) != number_of_children(two)) + return false; + return equal(begin(one), end(one), begin(two), fun); +} + +template +tree +tree::subtree(sibling_iterator from, + sibling_iterator to) const { + assert(from != + to); // if from==to, the range is empty, hence no tree to return. + + tree tmp; + tmp.set_head(value_type()); + tmp.replace(tmp.begin(), tmp.end(), from, to); + return tmp; +} + +template +void tree::subtree(tree &tmp, sibling_iterator from, + sibling_iterator to) const { + assert(from != + to); // if from==to, the range is empty, hence no tree to return. + + tmp.set_head(value_type()); + tmp.replace(tmp.begin(), tmp.end(), from, to); +} + +template +size_t tree::size() const { + size_t i = 0; + pre_order_iterator it = begin(), eit = end(); + while (it != eit) { + ++i; + ++it; + } + return i; +} + +template +size_t tree::size(const iterator_base &top) const { + size_t i = 0; + pre_order_iterator it = top, eit = top; + eit.skip_children(); + ++eit; + while (it != eit) { + ++i; + ++it; + } + return i; +} + +template +bool tree::empty() const { + pre_order_iterator it = begin(), eit = end(); + return (it == eit); +} + +template +int tree::depth(const iterator_base &it) { + tree_node *pos = it.node; + assert(pos != 0); + int ret = 0; + while (pos->parent != 0) { + pos = pos->parent; + ++ret; + } + return ret; +} + +template +int tree::depth(const iterator_base &it, + const iterator_base &root) { + tree_node *pos = it.node; + assert(pos != 0); + int ret = 0; + while (pos->parent != 0 && pos != root.node) { + pos = pos->parent; + ++ret; + } + return ret; +} + +template +template +int tree::depth(const iterator_base &it, Predicate p) { + tree_node *pos = it.node; + assert(pos != 0); + int ret = 0; + while (pos->parent != 0) { + pos = pos->parent; + if (p(pos)) + ++ret; + } + return ret; +} + +template +template +int tree::distance(const iterator_base &top, + const iterator_base &bottom, + Predicate p) { + tree_node *pos = bottom.node; + assert(pos != 0); + int ret = 0; + while (pos->parent != 0 && pos != top.node) { + pos = pos->parent; + if (p(pos)) + ++ret; + } + return ret; +} + +template +int tree::max_depth() const { + int maxd = -1; + for (tree_node *it = head->next_sibling; it != feet; it = it->next_sibling) + maxd = std::max(maxd, max_depth(it)); + + return maxd; +} + +template +int tree::max_depth(const iterator_base &pos) const { + tree_node *tmp = pos.node; + + if (tmp == 0 || tmp == head || tmp == feet) + return -1; + + int curdepth = 0, maxdepth = 0; + while (true) { // try to walk the bottom of the tree + while (tmp->first_child == 0) { + if (tmp == pos.node) + return maxdepth; + if (tmp->next_sibling == 0) { + // try to walk up and then right again + do { + tmp = tmp->parent; + if (tmp == 0) + return maxdepth; + --curdepth; + } while (tmp->next_sibling == 0); + } + if (tmp == pos.node) + return maxdepth; + tmp = tmp->next_sibling; + } + tmp = tmp->first_child; + ++curdepth; + maxdepth = std::max(curdepth, maxdepth); + } +} + +template +unsigned int +tree::number_of_children(const iterator_base &it) { + tree_node *pos = it.node->first_child; + if (pos == 0) + return 0; + + unsigned int ret = 1; + // while(pos!=it.node->last_child) { + // ++ret; + // pos=pos->next_sibling; + // } + while ((pos = pos->next_sibling)) + ++ret; + return ret; +} + +template +unsigned int tree::number_of_siblings( + const iterator_base &it) const { + tree_node *pos = it.node; + unsigned int ret = 0; + // count forward + while (pos->next_sibling && pos->next_sibling != head && + pos->next_sibling != feet) { + ++ret; + pos = pos->next_sibling; + } + // count backward + pos = it.node; + while (pos->prev_sibling && pos->prev_sibling != head && + pos->prev_sibling != feet) { + ++ret; + pos = pos->prev_sibling; + } + + return ret; +} + +template +void tree::swap(sibling_iterator it) { + tree_node *nxt = it.node->next_sibling; + if (nxt) { + if (it.node->prev_sibling) + it.node->prev_sibling->next_sibling = nxt; + else + it.node->parent->first_child = nxt; + nxt->prev_sibling = it.node->prev_sibling; + tree_node *nxtnxt = nxt->next_sibling; + if (nxtnxt) + nxtnxt->prev_sibling = it.node; + else + it.node->parent->last_child = it.node; + nxt->next_sibling = it.node; + it.node->prev_sibling = nxt; + it.node->next_sibling = nxtnxt; + } +} + +template +void tree::swap(iterator one, iterator two) { + // if one and two are adjacent siblings, use the sibling swap + if (one.node->next_sibling == two.node) + swap(one); + else if (two.node->next_sibling == one.node) + swap(two); + else { + tree_node *nxt1 = one.node->next_sibling; + tree_node *nxt2 = two.node->next_sibling; + tree_node *pre1 = one.node->prev_sibling; + tree_node *pre2 = two.node->prev_sibling; + tree_node *par1 = one.node->parent; + tree_node *par2 = two.node->parent; + + // reconnect + one.node->parent = par2; + one.node->next_sibling = nxt2; + if (nxt2) + nxt2->prev_sibling = one.node; + else + par2->last_child = one.node; + one.node->prev_sibling = pre2; + if (pre2) + pre2->next_sibling = one.node; + else + par2->first_child = one.node; + + two.node->parent = par1; + two.node->next_sibling = nxt1; + if (nxt1) + nxt1->prev_sibling = two.node; + else + par1->last_child = two.node; + two.node->prev_sibling = pre1; + if (pre1) + pre1->next_sibling = two.node; + else + par1->first_child = two.node; + } +} + +// template +// tree::iterator tree::find_subtree( sibling_iterator subfrom, +// sibling_iterator subto, iterator from, iterator to, BinaryPredicate fun) +// const +// { +// assert(1==0); // this routine is not finished yet. +// while(from!=to) { +// if(fun(*subfrom, *from)) { +// +// } +// } +// return to; +// } + +template +bool tree::is_in_subtree( + const iterator_base &it, const iterator_base &top) const { + sibling_iterator first = top; + sibling_iterator last = first; + ++last; + return is_in_subtree(it, first, last); +} + +template +bool tree::is_in_subtree( + const iterator_base &it, const iterator_base &begin, + const iterator_base &end) const { + // FIXME: this should be optimised. + pre_order_iterator tmp = begin; + while (tmp != end) { + if (tmp == it) + return true; + ++tmp; + } + return false; +} + +template +bool tree::is_valid(const iterator_base &it) const { + if (it.node == 0 || it.node == feet || it.node == head) + return false; + else + return true; +} + +template +bool tree::is_head(const iterator_base &it) { + if (it.node->parent == 0) + return true; + return false; +} + +template +typename tree::iterator +tree::lowest_common_ancestor( + const iterator_base &one, const iterator_base &two) const { + std::set parents; + + // Walk up from 'one' storing all parents. + iterator walk = one; + do { + walk = parent(walk); + parents.insert(walk); + } while (walk.node->parent); + + // Walk up from 'two' until we encounter a node in parents. + walk = two; + do { + walk = parent(walk); + if (parents.find(walk) != parents.end()) + break; + } while (walk.node->parent); + + return walk; +} + +template +unsigned int tree::index(sibling_iterator it) const { + unsigned int ind = 0; + if (it.node->parent == 0) { + while (it.node->prev_sibling != head) { + it.node = it.node->prev_sibling; + ++ind; + } + } else { + while (it.node->prev_sibling != 0) { + it.node = it.node->prev_sibling; + ++ind; + } + } + return ind; +} + +template +typename tree::sibling_iterator +tree::sibling(const iterator_base &it, + unsigned int num) const { + tree_node *tmp; + if (it.node->parent == 0) { + tmp = head->next_sibling; + while (num) { + tmp = tmp->next_sibling; + --num; + } + } else { + tmp = it.node->parent->first_child; + while (num) { + assert(tmp != 0); + tmp = tmp->next_sibling; + --num; + } + } + return tmp; +} + +template +void tree::debug_verify_consistency() const { + iterator it = begin(); + while (it != end()) { + // std::cerr << *it << " (" << it.node << ")" << std::endl; + if (it.node->parent != 0) { + if (it.node->prev_sibling == 0) + assert(it.node->parent->first_child == it.node); + else + assert(it.node->prev_sibling->next_sibling == it.node); + if (it.node->next_sibling == 0) + assert(it.node->parent->last_child == it.node); + else + assert(it.node->next_sibling->prev_sibling == it.node); + } + ++it; + } +} + +template +typename tree::sibling_iterator +tree::child(const iterator_base &it, unsigned int num) { + tree_node *tmp = it.node->first_child; + while (num--) { + assert(tmp != 0); + tmp = tmp->next_sibling; + } + return tmp; +} + +// Iterator base + +template +tree::iterator_base::iterator_base() + : node(0), skip_current_children_(false) {} + +template +tree::iterator_base::iterator_base(tree_node *tn) + : node(tn), skip_current_children_(false) {} + +template +T &tree::iterator_base::operator*() const { + return node->data; +} + +template +T *tree::iterator_base::operator->() const { + return &(node->data); +} + +template +bool tree::post_order_iterator::operator!=( + const post_order_iterator &other) const { + if (other.node != this->node) + return true; + else + return false; +} + +template +bool tree::post_order_iterator::operator==( + const post_order_iterator &other) const { + if (other.node == this->node) + return true; + else + return false; +} + +template +bool tree::pre_order_iterator::operator!=( + const pre_order_iterator &other) const { + if (other.node != this->node) + return true; + else + return false; +} + +template +bool tree::pre_order_iterator::operator==( + const pre_order_iterator &other) const { + if (other.node == this->node) + return true; + else + return false; +} + +template +bool tree::sibling_iterator::operator!=( + const sibling_iterator &other) const { + if (other.node != this->node) + return true; + else + return false; +} + +template +bool tree::sibling_iterator::operator==( + const sibling_iterator &other) const { + if (other.node == this->node) + return true; + else + return false; +} + +template +bool tree::leaf_iterator::operator!=( + const leaf_iterator &other) const { + if (other.node != this->node) + return true; + else + return false; +} + +template +bool tree::leaf_iterator::operator==( + const leaf_iterator &other) const { + if (other.node == this->node && other.top_node == this->top_node) + return true; + else + return false; +} + +template +typename tree::sibling_iterator +tree::iterator_base::begin() const { + if (node->first_child == 0) + return end(); + + sibling_iterator ret(node->first_child); + ret.parent_ = this->node; + return ret; +} + +template +typename tree::sibling_iterator +tree::iterator_base::end() const { + sibling_iterator ret(0); + ret.parent_ = node; + return ret; +} + +template +void tree::iterator_base::skip_children() { + skip_current_children_ = true; +} + +template +void tree::iterator_base::skip_children(bool skip) { + skip_current_children_ = skip; +} + +template +unsigned int +tree::iterator_base::number_of_children() const { + tree_node *pos = node->first_child; + if (pos == 0) + return 0; + + unsigned int ret = 1; + while (pos != node->last_child) { + ++ret; + pos = pos->next_sibling; + } + return ret; +} + +// Pre-order iterator + +template +tree::pre_order_iterator::pre_order_iterator() + : iterator_base(0) {} + +template +tree::pre_order_iterator::pre_order_iterator( + tree_node *tn) + : iterator_base(tn) {} + +template +tree::pre_order_iterator::pre_order_iterator( + const iterator_base &other) + : iterator_base(other.node) {} + +template +tree::pre_order_iterator::pre_order_iterator( + const sibling_iterator &other) + : iterator_base(other.node) { + if (this->node == 0) { + if (other.range_last() != 0) + this->node = other.range_last(); + else + this->node = other.parent_; + this->skip_children(); + ++(*this); + } +} + +template +typename tree::pre_order_iterator & +tree::pre_order_iterator::operator++() { + assert(this->node != 0); + if (!this->skip_current_children_ && this->node->first_child != 0) { + this->node = this->node->first_child; + } else { + this->skip_current_children_ = false; + while (this->node->next_sibling == 0) { + this->node = this->node->parent; + if (this->node == 0) + return *this; + } + this->node = this->node->next_sibling; + } + return *this; +} + +template +typename tree::pre_order_iterator & +tree::pre_order_iterator::operator--() { + assert(this->node != 0); + if (this->node->prev_sibling) { + this->node = this->node->prev_sibling; + while (this->node->last_child) + this->node = this->node->last_child; + } else { + this->node = this->node->parent; + if (this->node == 0) + return *this; + } + return *this; +} + +template +typename tree::pre_order_iterator +tree::pre_order_iterator::operator++(int) { + pre_order_iterator copy = *this; + ++(*this); + return copy; +} + +template +typename tree::pre_order_iterator & +tree::pre_order_iterator::next_skip_children() { + (*this).skip_children(); + (*this)++; + return *this; +} + +template +typename tree::pre_order_iterator +tree::pre_order_iterator::operator--(int) { + pre_order_iterator copy = *this; + --(*this); + return copy; +} + +template +typename tree::pre_order_iterator & +tree::pre_order_iterator::operator+=(unsigned int num) { + while (num > 0) { + ++(*this); + --num; + } + return (*this); +} + +template +typename tree::pre_order_iterator & +tree::pre_order_iterator::operator-=(unsigned int num) { + while (num > 0) { + --(*this); + --num; + } + return (*this); +} + +// Post-order iterator + +template +tree::post_order_iterator::post_order_iterator() + : iterator_base(0) {} + +template +tree::post_order_iterator::post_order_iterator( + tree_node *tn) + : iterator_base(tn) {} + +template +tree::post_order_iterator::post_order_iterator( + const iterator_base &other) + : iterator_base(other.node) {} + +template +tree::post_order_iterator::post_order_iterator( + const sibling_iterator &other) + : iterator_base(other.node) { + if (this->node == 0) { + if (other.range_last() != 0) + this->node = other.range_last(); + else + this->node = other.parent_; + this->skip_children(); + ++(*this); + } +} + +template +typename tree::post_order_iterator & +tree::post_order_iterator::operator++() { + assert(this->node != 0); + if (this->node->next_sibling == 0) { + this->node = this->node->parent; + this->skip_current_children_ = false; + } else { + this->node = this->node->next_sibling; + if (this->skip_current_children_) { + this->skip_current_children_ = false; + } else { + while (this->node->first_child) + this->node = this->node->first_child; + } + } + return *this; +} + +template +typename tree::post_order_iterator & +tree::post_order_iterator::operator--() { + assert(this->node != 0); + if (this->skip_current_children_ || this->node->last_child == 0) { + this->skip_current_children_ = false; + while (this->node->prev_sibling == 0) + this->node = this->node->parent; + this->node = this->node->prev_sibling; + } else { + this->node = this->node->last_child; + } + return *this; +} + +template +typename tree::post_order_iterator +tree::post_order_iterator::operator++(int) { + post_order_iterator copy = *this; + ++(*this); + return copy; +} + +template +typename tree::post_order_iterator +tree::post_order_iterator::operator--(int) { + post_order_iterator copy = *this; + --(*this); + return copy; +} + +template +typename tree::post_order_iterator & +tree::post_order_iterator::operator+=( + unsigned int num) { + while (num > 0) { + ++(*this); + --num; + } + return (*this); +} + +template +typename tree::post_order_iterator & +tree::post_order_iterator::operator-=( + unsigned int num) { + while (num > 0) { + --(*this); + --num; + } + return (*this); +} + +template +void tree::post_order_iterator::descend_all() { + assert(this->node != 0); + while (this->node->first_child) + this->node = this->node->first_child; +} + +// Breadth-first iterator + +template +tree::breadth_first_queued_iterator:: + breadth_first_queued_iterator() + : iterator_base() {} + +template +tree::breadth_first_queued_iterator:: + breadth_first_queued_iterator(tree_node *tn) + : iterator_base(tn) { + traversal_queue.push(tn); +} + +template +tree::breadth_first_queued_iterator:: + breadth_first_queued_iterator(const iterator_base &other) + : iterator_base(other.node) { + traversal_queue.push(other.node); +} + +template +bool tree::breadth_first_queued_iterator::operator!=( + const breadth_first_queued_iterator &other) const { + if (other.node != this->node) + return true; + else + return false; +} + +template +bool tree::breadth_first_queued_iterator::operator==( + const breadth_first_queued_iterator &other) const { + if (other.node == this->node) + return true; + else + return false; +} + +template +typename tree::breadth_first_queued_iterator & +tree::breadth_first_queued_iterator::operator++() { + assert(this->node != 0); + + // Add child nodes and pop current node + sibling_iterator sib = this->begin(); + while (sib != this->end()) { + traversal_queue.push(sib.node); + ++sib; + } + traversal_queue.pop(); + if (traversal_queue.size() > 0) + this->node = traversal_queue.front(); + else + this->node = 0; + return (*this); +} + +template +typename tree::breadth_first_queued_iterator +tree::breadth_first_queued_iterator::operator++(int) { + breadth_first_queued_iterator copy = *this; + ++(*this); + return copy; +} + +template +typename tree::breadth_first_queued_iterator & +tree::breadth_first_queued_iterator::operator+=( + unsigned int num) { + while (num > 0) { + ++(*this); + --num; + } + return (*this); +} + +// Fixed depth iterator + +template +tree::fixed_depth_iterator::fixed_depth_iterator() + : iterator_base() {} + +template +tree::fixed_depth_iterator::fixed_depth_iterator( + tree_node *tn) + : iterator_base(tn), top_node(0) {} + +template +tree::fixed_depth_iterator::fixed_depth_iterator( + const iterator_base &other) + : iterator_base(other.node), top_node(0) {} + +template +tree::fixed_depth_iterator::fixed_depth_iterator( + const sibling_iterator &other) + : iterator_base(other.node), top_node(0) {} + +template +tree::fixed_depth_iterator::fixed_depth_iterator( + const fixed_depth_iterator &other) + : iterator_base(other.node), top_node(other.top_node) {} + +template +void tree::fixed_depth_iterator::swap( + fixed_depth_iterator &first, fixed_depth_iterator &second) { + std::swap(first.node, second.node); + std::swap(first.skip_current_children_, second.skip_current_children_); + std::swap(first.top_node, second.top_node); +} + +template +typename tree::fixed_depth_iterator & +tree::fixed_depth_iterator::operator=( + fixed_depth_iterator other) { + swap(*this, other); + return *this; +} + +template +bool tree::fixed_depth_iterator::operator==( + const fixed_depth_iterator &other) const { + if (other.node == this->node && other.top_node == top_node) + return true; + else + return false; +} + +template +bool tree::fixed_depth_iterator::operator!=( + const fixed_depth_iterator &other) const { + if (other.node != this->node || other.top_node != top_node) + return true; + else + return false; +} + +template +typename tree::fixed_depth_iterator & +tree::fixed_depth_iterator::operator++() { + assert(this->node != 0); + + if (this->node->next_sibling) { + this->node = this->node->next_sibling; + } else { + int relative_depth = 0; + upper: + do { + if (this->node == this->top_node) { + this->node = 0; // FIXME: return a proper fixed_depth end iterator once + // implemented + return *this; + } + this->node = this->node->parent; + if (this->node == 0) + return *this; + --relative_depth; + } while (this->node->next_sibling == 0); + lower: + this->node = this->node->next_sibling; + while (this->node->first_child == 0) { + if (this->node->next_sibling == 0) + goto upper; + this->node = this->node->next_sibling; + if (this->node == 0) + return *this; + } + while (relative_depth < 0 && this->node->first_child != 0) { + this->node = this->node->first_child; + ++relative_depth; + } + if (relative_depth < 0) { + if (this->node->next_sibling == 0) + goto upper; + else + goto lower; + } + } + return *this; +} + +template +typename tree::fixed_depth_iterator & +tree::fixed_depth_iterator::operator--() { + assert(this->node != 0); + + if (this->node->prev_sibling) { + this->node = this->node->prev_sibling; + } else { + int relative_depth = 0; + upper: + do { + if (this->node == this->top_node) { + this->node = 0; + return *this; + } + this->node = this->node->parent; + if (this->node == 0) + return *this; + --relative_depth; + } while (this->node->prev_sibling == 0); + lower: + this->node = this->node->prev_sibling; + while (this->node->last_child == 0) { + if (this->node->prev_sibling == 0) + goto upper; + this->node = this->node->prev_sibling; + if (this->node == 0) + return *this; + } + while (relative_depth < 0 && this->node->last_child != 0) { + this->node = this->node->last_child; + ++relative_depth; + } + if (relative_depth < 0) { + if (this->node->prev_sibling == 0) + goto upper; + else + goto lower; + } + } + return *this; + + // + // + // assert(this->node!=0); + // if(this->node->prev_sibling!=0) { + // this->node=this->node->prev_sibling; + // assert(this->node!=0); + // if(this->node->parent==0 && this->node->prev_sibling==0) // head + // element this->node=0; + // } + // else { + // tree_node *par=this->node->parent; + // do { + // par=par->prev_sibling; + // if(par==0) { // FIXME: need to keep track of this! + // this->node=0; + // return *this; + // } + // } while(par->last_child==0); + // this->node=par->last_child; + // } + // return *this; +} + +template +typename tree::fixed_depth_iterator +tree::fixed_depth_iterator::operator++(int) { + fixed_depth_iterator copy = *this; + ++(*this); + return copy; +} + +template +typename tree::fixed_depth_iterator +tree::fixed_depth_iterator::operator--(int) { + fixed_depth_iterator copy = *this; + --(*this); + return copy; +} + +template +typename tree::fixed_depth_iterator & +tree::fixed_depth_iterator::operator-=( + unsigned int num) { + while (num > 0) { + --(*this); + --(num); + } + return (*this); +} + +template +typename tree::fixed_depth_iterator & +tree::fixed_depth_iterator::operator+=( + unsigned int num) { + while (num > 0) { + ++(*this); + --(num); + } + return *this; +} + +// Sibling iterator + +template +tree::sibling_iterator::sibling_iterator() + : iterator_base() { + set_parent_(); +} + +template +tree::sibling_iterator::sibling_iterator(tree_node *tn) + : iterator_base(tn) { + set_parent_(); +} + +template +tree::sibling_iterator::sibling_iterator( + const iterator_base &other) + : iterator_base(other.node) { + set_parent_(); +} + +template +tree::sibling_iterator::sibling_iterator( + const sibling_iterator &other) + : iterator_base(other), parent_(other.parent_) {} + +template +void tree::sibling_iterator::swap( + sibling_iterator &first, sibling_iterator &second) { + std::swap(first.node, second.node); + std::swap(first.skip_current_children_, second.skip_current_children_); + std::swap(first.parent_, second.parent_); +} + +template +typename tree::sibling_iterator & +tree::sibling_iterator::operator=( + sibling_iterator other) { + swap(*this, other); + return *this; +} + +template +void tree::sibling_iterator::set_parent_() { + parent_ = 0; + if (this->node == 0) + return; + if (this->node->parent != 0) + parent_ = this->node->parent; +} + +template +typename tree::sibling_iterator & +tree::sibling_iterator::operator++() { + if (this->node) + this->node = this->node->next_sibling; + return *this; +} + +template +typename tree::sibling_iterator & +tree::sibling_iterator::operator--() { + if (this->node) + this->node = this->node->prev_sibling; + else { + assert(parent_); + this->node = parent_->last_child; + } + return *this; +} + +template +typename tree::sibling_iterator +tree::sibling_iterator::operator++(int) { + sibling_iterator copy = *this; + ++(*this); + return copy; +} + +template +typename tree::sibling_iterator +tree::sibling_iterator::operator--(int) { + sibling_iterator copy = *this; + --(*this); + return copy; +} + +template +typename tree::sibling_iterator & +tree::sibling_iterator::operator+=(unsigned int num) { + while (num > 0) { + ++(*this); + --num; + } + return (*this); +} + +template +typename tree::sibling_iterator & +tree::sibling_iterator::operator-=(unsigned int num) { + while (num > 0) { + --(*this); + --num; + } + return (*this); +} + +template +typename tree::tree_node * +tree::sibling_iterator::range_first() const { + tree_node *tmp = parent_->first_child; + return tmp; +} + +template +typename tree::tree_node * +tree::sibling_iterator::range_last() const { + return parent_->last_child; +} + +// Leaf iterator + +template +tree::leaf_iterator::leaf_iterator() + : iterator_base(0), top_node(0) {} + +template +tree::leaf_iterator::leaf_iterator(tree_node *tn, + tree_node *top) + : iterator_base(tn), top_node(top) {} + +template +tree::leaf_iterator::leaf_iterator( + const iterator_base &other) + : iterator_base(other.node), top_node(0) {} + +template +tree::leaf_iterator::leaf_iterator( + const sibling_iterator &other) + : iterator_base(other.node), top_node(0) { + if (this->node == 0) { + if (other.range_last() != 0) + this->node = other.range_last(); + else + this->node = other.parent_; + ++(*this); + } +} + +template +typename tree::leaf_iterator & +tree::leaf_iterator::operator++() { + assert(this->node != 0); + if (this->node->first_child != + 0) { // current node is no longer leaf (children got added) + while (this->node->first_child) + this->node = this->node->first_child; + } else { + while (this->node->next_sibling == 0) { + if (this->node->parent == 0) + return *this; + this->node = this->node->parent; + if (top_node != 0 && this->node == top_node) + return *this; + } + this->node = this->node->next_sibling; + while (this->node->first_child) + this->node = this->node->first_child; + } + return *this; +} + +template +typename tree::leaf_iterator & +tree::leaf_iterator::operator--() { + assert(this->node != 0); + while (this->node->prev_sibling == 0) { + if (this->node->parent == 0) + return *this; + this->node = this->node->parent; + if (top_node != 0 && this->node == top_node) + return *this; + } + this->node = this->node->prev_sibling; + while (this->node->last_child) + this->node = this->node->last_child; + return *this; +} + +template +typename tree::leaf_iterator +tree::leaf_iterator::operator++(int) { + leaf_iterator copy = *this; + ++(*this); + return copy; +} + +template +typename tree::leaf_iterator +tree::leaf_iterator::operator--(int) { + leaf_iterator copy = *this; + --(*this); + return copy; +} + +template +typename tree::leaf_iterator & +tree::leaf_iterator::operator+=(unsigned int num) { + while (num > 0) { + ++(*this); + --num; + } + return (*this); +} + +template +typename tree::leaf_iterator & +tree::leaf_iterator::operator-=(unsigned int num) { + while (num > 0) { + --(*this); + --num; + } + return (*this); +} + +#endif + +// Local variables: +// tab-width: 3 +// End: diff --git a/core/opengate_core/opengate_lib/tree_util.hh b/core/opengate_core/opengate_lib/tree_util.hh new file mode 100644 index 000000000..2706b3e2a --- /dev/null +++ b/core/opengate_core/opengate_lib/tree_util.hh @@ -0,0 +1,90 @@ +/* + + A collection of miscellaneous utilities that operate on the templated + tree.hh class. + + + Copyright (C) 2001-2009 Kasper Peeters + + (At the moment this only contains a printing utility, thanks to Linda + Buisman ) + + This program is free software: you can redistribute it and/or + modify it under the terms of the GNU General Public License as + published by the Free Software Foundation, either version 3 of the + License, or (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . + +*/ + +#ifndef tree_util_hh_ +#define tree_util_hh_ + +#include "tree.hh" +#include + +namespace kptree { + +template +void print_tree_bracketed(const tree &t, std::ostream &str = std::cout); + +template +void print_subtree_bracketed(const tree &t, typename tree::iterator iRoot, + std::ostream &str = std::cout); + +// Iterate over all roots (the head) and print each one on a new line +// by calling printSingleRoot. + +template +void print_tree_bracketed(const tree &t, std::ostream &str) { + int headCount = t.number_of_siblings(t.begin()); + int headNum = 0; + for (typename tree::sibling_iterator iRoots = t.begin(); iRoots != t.end(); + ++iRoots, ++headNum) { + print_subtree_bracketed(t, iRoots, str); + if (headNum != headCount) { + str << std::endl; + } + } +} + +// Print everything under this root in a flat, bracketed structure. + +template +void print_subtree_bracketed(const tree &t, typename tree::iterator iRoot, + std::ostream &str) { + if (t.empty()) + return; + if (t.number_of_children(iRoot) == 0) { + str << *iRoot; + } else { + // parent + str << *iRoot; + str << "("; + // child1, ..., childn + int siblingCount = t.number_of_siblings(t.begin(iRoot)); + int siblingNum; + typename tree::sibling_iterator iChildren; + for (iChildren = t.begin(iRoot), siblingNum = 0; iChildren != t.end(iRoot); + ++iChildren, ++siblingNum) { + // recursively print child + print_subtree_bracketed(t, iChildren, str); + // comma after every child except the last one + if (siblingNum != siblingCount) { + str << ", "; + } + } + str << ")"; + } +} + +} // namespace kptree + +#endif diff --git a/opengate/actors/miscactors.py b/opengate/actors/miscactors.py index 1edc93fb8..e2fe2c16e 100644 --- a/opengate/actors/miscactors.py +++ b/opengate/actors/miscactors.py @@ -1,12 +1,15 @@ from box import Box import platform +from anytree import Node, RenderTree import opengate_core as g4 +from anytree import Node, RenderTree from .base import ActorBase from ..utility import g4_units, g4_best_unit_tuple from .actoroutput import ActorOutputBase, ActorOutputSingleImage from ..serialization import dump_json from ..exception import fatal, warning from ..base import process_cls +from anytree import RenderTree """ It is feasible to get callback every Run, Event, Track, Step in the python side. @@ -326,17 +329,13 @@ class KillAccordingProcessesActor(ActorBase, g4.GateKillAccordingProcessesActor) }, ), } + user_output_config = { "kill_according_processes": { "actor_output_class": ActorOutputKillAccordingProcessesActor, }, } - """ - If a particle, not generated or generated within the volume at which our actor is attached, crosses the volume - without interaction, the particle is killed. - """ - def __init__(self, *args, **kwargs): ActorBase.__init__(self, *args, **kwargs) self.__initcpp__() @@ -371,6 +370,52 @@ def __str__(self): return s +class KillAccordingParticleNameActor(ActorBase, g4.GateKillAccordingParticleNameActor): + """Actor which kills a particle according the particle name provied by the user at the exit of the + actorified volume.""" + + particles_name_to_kill: list + + user_info_defaults = { + "particles_name_to_kill": ( + [], + { + "doc": "Put particles name the user wants to kill at the exit of the volume" + }, + ), + } + + def __init__(self, *args, **kwargs): + ActorBase.__init__(self, *args, **kwargs) + self.number_of_killed_particles = 0 + self.__initcpp__() + self.list_of_volume_name = [] + + def __initcpp__(self): + g4.GateKillAccordingParticleNameActor.__init__(self, self.user_info) + self.AddActions( + {"PreUserTrackingAction", "SteppingAction", "EndSimulationAction"} + ) + + def initialize(self): + ActorBase.initialize(self) + self.InitializeUserInfo(self.user_info) + self.InitializeCpp() + volume_tree = self.simulation.volume_manager.get_volume_tree() + dico_of_volume_tree = {} + for pre, _, node in RenderTree(volume_tree): + dico_of_volume_tree[str(node.name)] = node + volume_name = self.user_info.attached_to + while volume_name != "world": + node = dico_of_volume_tree[volume_name] + volume_name = node.mother + self.list_of_volume_name.append(volume_name) + self.fListOfVolumeAncestor = self.list_of_volume_name + + def EndSimulationAction(self): + self.number_of_killed_particles = self.GetNumberOfKilledParticles() + + class KillActor(ActorBase, g4.GateKillActor): """Actor which kills a particle entering a volume.""" @@ -394,6 +439,331 @@ def EndSimulationAction(self): self.number_of_killed_particles = self.GetNumberOfKilledParticles() +class ActorOutputKillNonInteractingParticleActor(ActorOutputBase): + + def __init__(self, *args, **kwargs): + super().__init__(*args, **kwargs) + self.number_of_killed_particles = 0 + + def get_processed_output(self): + d = {} + d["particles killed"] = self.number_of_killed_particles + return d + + def __str__(self): + s = "" + for k, v in self.get_processed_output().items(): + s = k + ": " + str(v) + s += "\n" + return s + + +class KillNonInteractingParticleActor( + ActorBase, g4.GateKillNonInteractingParticleActor +): + """ + If a particle, not generated or generated within the volume at which our actor is attached, crosses the volume + without interaction, the particle is killed. Warning : this actor being based on energy measurement, Rayleigh photon + may not be killed. + """ + + def __init__(self, *args, **kwargs): + ActorBase.__init__(self, *args, **kwargs) + self._add_user_output( + ActorOutputKillNonInteractingParticleActor, "kill_non_interacting_particles" + ) + self.__initcpp__() + self.list_of_volume_name = [] + self.number_of_killed_particles = 0 + + def __initcpp__(self): + g4.GateKillNonInteractingParticleActor.__init__(self, self.user_info) + self.AddActions( + { + "StartSimulationAction", + "PreUserTrackingAction", + "SteppingAction", + "EndOfSimulationAction", + } + ) + + def initialize(self): + ActorBase.initialize(self) + self.InitializeUserInfo(self.user_info) + self.InitializeCpp() + volume_tree = self.simulation.volume_manager.get_volume_tree() + dico_of_volume_tree = {} + for pre, _, node in RenderTree(volume_tree): + dico_of_volume_tree[str(node.name)] = node + volume_name = self.user_info.attached_to + while volume_name != "world": + node = dico_of_volume_tree[volume_name] + volume_name = node.mother + self.list_of_volume_name.append(volume_name) + self.fListOfVolumeAncestor = self.list_of_volume_name + + def EndSimulationAction(self): + self.user_output.kill_non_interacting_particles.number_of_killed_particles = ( + self.number_of_killed_particles + ) + + def __str__(self): + s = self.user_output["kill_non_interacting_particles"].__str__() + return s + + +def _setter_hook_particles(self, value): + if isinstance(value, str): + return [value] + else: + return list(value) + + +class SplittingActorBase(ActorBase): + """ + Actors based on the G4GenericBiasing class of GEANT4. This class provides tools to interact with GEANT4 processes + during a simulation, allowing direct modification of process properties. Additionally, it enables non-physics-based + particle splitting (e.g., pure geometrical splitting) to introduce biasing into simulations. SplittingActorBase + serves as a foundational class for particle splitting operations, with parameters for configuring the splitting + behavior based on various conditions. + """ + + # hints for IDE + splitting_factor: int + bias_primary_only: bool + bias_only_once: bool + particles: list + + user_info_defaults = { + "splitting_factor": ( + 1, + { + "doc": "Specifies the number of particles to generate each time the splitting mechanism is applied", + }, + ), + "bias_primary_only": ( + True, + { + "doc": "If true, the splitting mechanism is applied only to particles with a ParentID of 1", + }, + ), + "bias_only_once": ( + True, + { + "doc": "If true, the splitting mechanism is applied only once per particle history", + }, + ), + "particles": ( + [ + "all", + ], + { + "doc": "Specifies the particles to split. The default value, all, includes all particles", + "setter_hook": _setter_hook_particles, + }, + ), + } + + +class ComptSplittingActor(SplittingActorBase, g4.GateOptrComptSplittingActor): + """ + This splitting actor enables process-based splitting specifically for Compton interactions. Each time a Compton + process occurs, its behavior is modified by generating multiple Compton scattering tracks + (splitting factor - 1 additional tracks plus the original) associated with the initial particle. + Compton electrons produced in the interaction are also included, in accordance with the secondary cut settings + provided by the user. + """ + + # hints for IDE + min_weight_of_particle: float + russian_roulette: bool + rotation_vector_director: bool + vector_director: list + max_theta: float + + user_info_defaults = { + "min_weight_of_particle": ( + 0, + { + "doc": "Defines a minimum weight for particles. Particles with weights below this threshold will not be split, limiting the splitting cascade of low-weight particles generated during Compton interactions.", + }, + ), + "russian_roulette": ( + False, + { + "doc": "If enabled (True), applies a Russian roulette mechanism. Particles emitted in undesired directions are discarded if a random number exceeds 1 / splitting_factor", + }, + ), + "vector_director": ( + [0, 0, 1], + { + "doc": "Specifies the particle’s direction of interest for the Russian roulette. In this direction, the Russian roulette is not applied", + }, + ), + "rotation_vector_director": ( + False, + { + "doc": "If enabled, allows the vector_director to rotate based on any rotation applied to a volume to which this actor is attached", + }, + ), + "max_theta": ( + 90 * g4_units.deg, + { + "doc": "Sets the angular range (in degrees) around vector_director within which the Russian roulette mechanism is not applied.", + }, + ), + } + + processes = ("compt",) + + def __init__(self, *args, **kwargs): + SplittingActorBase.__init__(self, *args, **kwargs) + self.__initcpp__() + + def __initcpp__(self): + g4.GateOptrComptSplittingActor.__init__(self, {"name": self.name}) + + def initialize(self): + SplittingActorBase.initialize(self) + self.InitializeUserInfo(self.user_info) + self.InitializeCpp() + + +class LastVertexInteractionSplittingActor( + ActorBase, g4.GateLastVertexInteractionSplittingActor +): + """This splitting actor proposes an interaction splitting at the last particle vertex before the exit + of the biased volume. This actor can be usefull for application where collimation are important, + such as in medical LINAC (Linear Accelerator) simulations or radiation shielding. + """ + + # hints for IDE + splitting_factor: int + angular_kill: bool + max_theta: float + vector_director: list + rotation_vector_director: bool + batch_size: int + nb_of_max_batch_per_event: int + + user_info_defaults = { + "splitting_factor": ( + 1, + { + "doc": "Defines the number of particles exiting at each split process. Unlike other split actors, this splitting factor counts particles that actually exit, not just those generated.", + }, + ), + "angular_kill": ( + False, + { + "doc": "If enabled, particles with momentum outside a specified angular range are killed.", + }, + ), + "max_theta": ( + 90 * g4_units.deg, + { + "doc": "Defines the maximum angle (in degrees) from a central axis within which particles are retained. Particles with momentum beyond this angle are removed. The angular range spans from 0 to max_theta, measured from the vector_director", + }, + ), + "vector_director": ( + [0, 0, 1], + { + "doc": "Specifies the reference direction vector from which max_theta is measured. Particles’ angular range is calculated based on this direction.", + }, + ), + "rotation_vector_director": ( + False, + { + "doc": "If enabled, the vector_director rotates in alignment with any rotation applied to the biased volume attached to this actor.", + }, + ), + "batch_size": ( + 1, + { + "doc": "Defines a batch of number of processes to regenerate. The optimal value depends on the collimation setup; for example, a batch_size of 10 works well for LINAC head configurations.", + }, + ), + "nb_of_max_batch_per_event": ( + 500, + { + "doc": "Defines a maximum number of attempt to enable the particles to exit. Useful to avoid an important loss of time for extremely rare events", + }, + ), + } + + def __init__(self, *args, **kwargs): + ActorBase.__init__(self, *args, **kwargs) + self.__initcpp__() + self.list_of_volume_name = [] + + def __initcpp__(self): + g4.GateLastVertexInteractionSplittingActor.__init__(self, {"name": self.name}) + self.AddActions( + { + "BeginOfRunAction", + "BeginOfEventAction", + "PreUserTrackingAction", + "SteppingAction", + "PostUserTrackingAction", + "EndOfEventAction", + "EndSimulationAction", + } + ) + + def initialize(self): + ActorBase.initialize(self) + self.InitializeUserInfo(self.user_info) + self.InitializeCpp() + volume_tree = self.simulation.volume_manager.get_volume_tree() + dico_of_volume_tree = {} + for pre, _, node in RenderTree(volume_tree): + dico_of_volume_tree[str(node.name)] = node + volume_name = self.user_info.attached_to + while volume_name != "world": + node = dico_of_volume_tree[volume_name] + volume_name = node.mother + self.list_of_volume_name.append(volume_name) + self.fListOfVolumeAncestor = self.list_of_volume_name + + def EndSimulationAction(self): + print("Number of replayed particles: ", self.GetNumberOfReplayedParticles()) + print("Number of killed particle:", self.GetNumberOfKilledParticles()) + + +class BremSplittingActor(SplittingActorBase, g4.GateBOptrBremSplittingActor): + """ + This splitting actor enables process-based splitting specifically for bremsstrahlung process. Each time a Brem + process occurs, its behavior is modified by generating multiple secondary Brem scattering tracks + (splitting factor) attached to the initial charged particle. + """ + + # hints for IDE + processes: list + + user_info_defaults = { + "processes": ( + ["eBrem"], + { + "doc": "Specifies the process split by this actor. This parameter is set to eBrem, as the actor is specifically developed for this process. It is recommended not to modify this setting.", + }, + ), + } + + processes = ("eBrem",) + + def __init__(self, *args, **kwargs): + SplittingActorBase.__init__(self, *args, **kwargs) + self.__initcpp__() + + def __initcpp__(self): + g4.GateBOptrBremSplittingActor.__init__(self, {"name": self.name}) + + def initialize(self): + SplittingActorBase.initialize(self) + self.InitializeUserInfo(self.user_info) + self.InitializeCpp() + + class AttenuationImageActor(ActorBase, g4.GateAttenuationImageActor): """ This actor generates an attenuation image for a simulation run. @@ -463,4 +833,9 @@ def BeginOfRunAction(self, run): process_cls(KillActor) process_cls(ActorOutputKillAccordingProcessesActor) process_cls(KillAccordingProcessesActor) +process_cls(LastVertexInteractionSplittingActor) +process_cls(KillNonInteractingParticleActor) +process_cls(SplittingActorBase) +process_cls(ComptSplittingActor) +process_cls(BremSplittingActor) process_cls(AttenuationImageActor) diff --git a/opengate/contrib/linacs/elektaversa.py b/opengate/contrib/linacs/elektaversa.py index bcaf1a2db..f9066c3a5 100644 --- a/opengate/contrib/linacs/elektaversa.py +++ b/opengate/contrib/linacs/elektaversa.py @@ -39,7 +39,7 @@ def add_empty_linac_box(sim, linac_name, sad=1000): linac.material = "G4_AIR" linac.size = [1 * m, 1 * m, 0.52 * m] translation_linac_box = np.array([0 * mm, 0, sad - linac.size[2] / 2 + 3.5 * mm]) - # Isocenter begin at the end of the target, That's why 3.5 mm is added. + # Isocenter start from the end of the target, That's why 3.5 mm is added. linac.translation = translation_linac_box linac.color = [1, 1, 1, 0] return linac @@ -539,7 +539,8 @@ def mlc_leaf(linac_name): interleaf_gap = 0.09 * mm leaf_length = 155 * mm leaf_height = 90 * mm - leaf_mean_width = 1.76 * mm + # leaf_mean_width = 1.76 * mm + leaf_mean_width = 1.69 * mm tongues_length = 0.8 * mm cyl = volumes.TubsVolume(name=f"{linac_name}_cylinder_leaf") @@ -552,8 +553,10 @@ def mlc_leaf(linac_name): trap_leaf = volumes.TrapVolume(name=f"{linac_name}_trap_leaf") dz = leaf_height / 2 dy1 = leaf_length / 2 - dx1 = 1.94 * mm / 2 - dx3 = 1.58 * mm / 2 + # dx1 = 1.94 * mm / 2 + # dx3 = 1.58 * mm / 2 + dx1 = 1.91 * mm / 2 + dx3 = 1.47 * mm / 2 theta = 0 alpha1 = 0 @@ -583,11 +586,13 @@ def mlc_leaf(linac_name): dy1 = leaf_length / 2 ##FIXME I need to remove 2 um to the tongues to avoid an overleap between leaves - dx1 = (interleaf_gap - 2.2 * 10 ** (-3) * mm) / 2 + # dx1 = (interleaf_gap - 2.2 * 10 ** (-3) * mm) / 2 + dx1 = (interleaf_gap - 3.3 * 10 ** (-3) * mm) / 2 dx3 = dx1 alpha1 = 0 alpha2 = alpha1 - theta = np.arctan((1.58 * mm - 1.94 * mm) * 0.5 / leaf_height) + # theta = np.arctan((1.58 * mm - 1.94 * mm) * 0.5 / leaf_height) + theta = np.arctan((1.47 * mm - 1.91 * mm) * 0.5 / leaf_height) phi = 0 dy2 = dy1 dx2 = dx1 @@ -624,17 +629,19 @@ def add_mlc(sim, linac_name): z_linac = linac.size[2] center_mlc = 349.3 * mm + 3.5 * mm interleaf_gap = 0.09 * mm - leaf_width = 1.76 * mm + # leaf_width = 1.76 * mm + leaf_width = 1.69 * mm leaf_lenght = 155 * mm nb_leaf = 160 - rotation_angle = np.arctan((1.94 * mm - 1.58 * mm) * 0.5 / leaf_height) + # rotation_angle = np.arctan((1.94 * mm - 1.58 * mm) * 0.5 / leaf_height) + rotation_angle = np.arctan((1.91 * mm - 1.47 * mm) * 0.5 / leaf_height) mlc = sim.add_volume("Box", f"{linac_name}_mlc") mlc_bank_rotation = Rotation.from_euler( - "X", np.arctan(3.25 / 349.3), degrees=False + "X", -np.arctan(3.25 / 349.3), degrees=False ).as_matrix() mlc.rotation = mlc_bank_rotation - mlc.size = [linac.size[0] - 2 * cm, linac.size[1] - 2 * cm, 100 * mm] + mlc.size = [linac.size[0] - 2 * cm, linac.size[1] - 2 * cm, 95 * mm] mlc.translation = np.array([0, 0, z_linac / 2 - center_mlc]) mlc.mother = linac_name mlc.color = [0, 0, 0, 0] diff --git a/opengate/managers.py b/opengate/managers.py index e30f6f6a0..e3f2c8c74 100644 --- a/opengate/managers.py +++ b/opengate/managers.py @@ -39,6 +39,7 @@ from .processing import dispatch_to_subprocess from .sources.generic import SourceBase, GenericSource +from .sources.lastvertexsources import LastVertexSource from .sources.phspsources import PhaseSpaceSource from .sources.voxelsources import VoxelSource from .sources.gansources import GANSource, GANPairsSource @@ -48,6 +49,7 @@ source_types = { "GenericSource": GenericSource, + "LastVertexSource": LastVertexSource, "PhaseSpaceSource": PhaseSpaceSource, "VoxelSource": VoxelSource, "GANSource": GANSource, @@ -89,6 +91,10 @@ SimulationStatisticsActor, KillActor, KillAccordingProcessesActor, + LastVertexInteractionSplittingActor, + KillNonInteractingParticleActor, + KillAccordingParticleNameActor, + SplittingActorBase, AttenuationImageActor, ) from .actors.biasingactors import ( @@ -129,6 +135,10 @@ "SimulationStatisticsActor": SimulationStatisticsActor, "KillActor": KillActor, "KillAccordingProcessesActor": KillAccordingProcessesActor, + "KillNonInteractingParticleActor": KillNonInteractingParticleActor, + "KillAccordingParticleNameActor": KillAccordingParticleNameActor, + "BremSplittingActor": BremSplittingActor, + "ComptSplittingActor": ComptSplittingActor, "DynamicGeometryActor": DynamicGeometryActor, "ARFActor": ARFActor, "ARFTrainingDatasetActor": ARFTrainingDatasetActor, @@ -142,6 +152,9 @@ "DigitizerProjectionActor": DigitizerProjectionActor, "DigitizerEnergyWindowsActor": DigitizerEnergyWindowsActor, "DigitizerHitsCollectionActor": DigitizerHitsCollectionActor, + "PhaseSpaceActor": PhaseSpaceActor, + "LastVertexInteractionSplittingActor": LastVertexInteractionSplittingActor, + "AttenuationImageActor": AttenuationImageActor, # biasing "BremSplittingActor": BremSplittingActor, "ComptSplittingActor": ComptSplittingActor, @@ -1224,12 +1237,24 @@ def dump_volume_types(self): s += f"{vt} " return s + def get_volume_tree(self): + return self.volume_tree_root + def print_volume_types(self): print(self.dump_volume_types()) def dump_material_database_names(self): return list(self.material_database.filenames) + def get_volume_tree(self): + return self.volume_tree_root + + def get_volume_tree(self): + return self.volume_tree_root + + def get_volume_tree(self): + return self.volume_tree_root + def print_material_database_names(self): print(self.dump_material_database_names()) diff --git a/opengate/sources/builders.py b/opengate/sources/builders.py new file mode 100644 index 000000000..6d14ecfe8 --- /dev/null +++ b/opengate/sources/builders.py @@ -0,0 +1,28 @@ +from .generic import GenericSource, TemplateSource, LastVertexSource +from .voxelsources import VoxelsSource +from .gansources import GANSource, GANPairsSource +from .beamsources import IonPencilBeamSource, TreatmentPlanPBSource +from .phspsources import PhaseSpaceSource +from .phidsources import PhotonFromIonDecaySource +from ..utility import make_builders + + +""" + List of source types: Generic, Voxels etc + + Energy spectra for beta+ emitters +""" + +source_type_names = { + GenericSource, + VoxelsSource, + GANSource, + GANPairsSource, + IonPencilBeamSource, + TemplateSource, + LastVertexSource, + PhaseSpaceSource, + TreatmentPlanPBSource, + PhotonFromIonDecaySource, +} +source_builders = make_builders(source_type_names) diff --git a/opengate/sources/generic.py b/opengate/sources/generic.py index b9b3e5f74..f7f115de5 100644 --- a/opengate/sources/generic.py +++ b/opengate/sources/generic.py @@ -360,4 +360,41 @@ def can_predict_number_of_events(self): return True +class TemplateSource(SourceBase): + """ + Source template: to create a new type of source, copy-paste + this file and adapt to your needs. + Also declare the source type in the file helpers_source.py + """ + + type_name = "TemplateSource" + + @staticmethod + def set_default_user_info(user_info): + SourceBase.set_default_user_info(user_info) + # initial user info + user_info.float_value = None + user_info.vector_value = None + + def create_g4_source(self): + return opengate_core.GateTemplateSource() + + def __init__(self, user_info): + super().__init__(user_info) + + def initialize(self, run_timing_intervals): + # Check user_info type + if self.user_info.float_value is None: + fatal( + f"Error for source {self.user_info.name}, float_value must be a float" + ) + if self.user_info.vector_value is None: + fatal( + f"Error for source {self.user_info.name}, vector_value must be a vector" + ) + + # initialize + SourceBase.initialize(self, run_timing_intervals) + + process_cls(GenericSource) diff --git a/opengate/sources/lastvertexsources.py b/opengate/sources/lastvertexsources.py new file mode 100644 index 000000000..6712cac08 --- /dev/null +++ b/opengate/sources/lastvertexsources.py @@ -0,0 +1,32 @@ +from box import Box +from scipy.spatial.transform import Rotation + +import opengate_core as g4 +from .base import ( + SourceBase, + all_beta_plus_radionuclides, + read_beta_plus_spectra, + compute_cdf_and_total_yield, +) +from ..base import process_cls +from ..utility import g4_units +from ..exception import fatal, warning + + +class LastVertexSource(SourceBase, g4.GateLastVertexSource): + """ + 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__() + + def __initcpp__(self): + g4.GateLastVertexSource.__init__(self) + + def initialize(self, run_timing_intervals): + SourceBase.initialize(self, run_timing_intervals) + + +process_cls(LastVertexSource) diff --git a/opengate/tests/src/test072_pseudo_transport_RR.py b/opengate/tests/src/test072_pseudo_transport_RR.py new file mode 100644 index 000000000..3ed0038e4 --- /dev/null +++ b/opengate/tests/src/test072_pseudo_transport_RR.py @@ -0,0 +1,231 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import opengate as gate +import uproot +import numpy as np +import matplotlib.pyplot as plt +import matplotlib as mpl +import scipy +from scipy.spatial.transform import Rotation +from opengate.tests import utility + + +def Ekin_per_event(weight, EventID, Ekin): + l_edep = [] + Edep_per_event = 0 + last_EventID = 0 + for i in range(len(weight)): + if i > 0: + last_EventID = EventID[i - 1] + if EventID[i] != last_EventID and i > 0: + l_edep.append(Edep_per_event) + Edep_per_event = 0 + Edep_per_event += Ekin[i] * weight[i] + return np.asarray(l_edep) + + +def validation_test(arr_no_RR, arr_RR, nb_event, splitting_factor, tol=0.15): + + Tracks_no_RR = arr_no_RR[ + (arr_no_RR["TrackCreatorProcess"] == "biasWrapper(compt)") + & (arr_no_RR["ParticleName"] == "gamma") + ] + + Tracks_RR = arr_RR[ + (arr_RR["TrackCreatorProcess"] == "biasWrapper(compt)") + & (arr_RR["ParticleName"] == "gamma") + ] + + Ekin_no_RR = Tracks_no_RR["KineticEnergy"] + w_no_RR = Tracks_no_RR["Weight"] + Event_ID_no_RR = Tracks_no_RR["EventID"] + + Ekin_per_event_no_RR = Ekin_per_event(w_no_RR, Event_ID_no_RR, Ekin_no_RR) + + Ekin_RR = Tracks_RR["KineticEnergy"] + w_RR = Tracks_RR["Weight"] + Event_ID_RR = Tracks_RR["EventID"] + + bool_RR = (np.min(w_RR) > 0.1 * 1 / splitting_factor) & ( + np.max(w_RR) < 1 / splitting_factor + ) + + Ekin_per_event_RR = Ekin_per_event(w_RR, Event_ID_RR, Ekin_RR) + + mean_ekin_event_no_RR = np.sum(Ekin_per_event_no_RR) / nb_event + mean_ekin_event_RR = np.sum(Ekin_per_event_RR) / nb_event + diff = (mean_ekin_event_RR - mean_ekin_event_no_RR) / mean_ekin_event_RR + print( + "Mean kinetic energy per event without russian roulette :", + np.round(1000 * mean_ekin_event_no_RR, 1), + "keV", + ) + print( + "Mean kinetic energy per event with russian roulette :", + np.round(1000 * mean_ekin_event_RR, 1), + "keV", + ) + print("Percentage of difference :", diff * 100, "%") + + return (abs(diff) < tol) & (bool_RR) + # mean_energy_per_history_no_RR = np.sum(Tracks_no_RR["KineticEnergy"]*Tracks_no_RR["Weight"])/nb_event + # mean_energy_per_history_RR = np.sum(Tracks_RR["KineticEnergy"] * Tracks_RR["Weight"]) / nb_event + + # print(mean_energy_per_history_no_RR,mean_energy_per_history_RR) + + +if __name__ == "__main__": + + for j in range(2): + paths = utility.get_default_test_paths( + __file__, "test072_pseudo_transportation", output_folder="test072" + ) + + # 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 + um = gate.g4_units.um + 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 + deg = gate.g4_units.deg + + # adapt world size + world = sim.world + world.size = [1.2 * m, 1.2 * m, 2 * m] + world.material = "G4_Galactic" + + ####### GEOMETRY TO IRRADIATE ############# + sim.volume_manager.material_database.add_material_weights( + "Tungsten", + ["W"], + [1], + 19.3 * gcm3, + ) + simple_collimation = sim.add_volume("Box", "colli_box") + simple_collimation.material = "G4_Galactic" + simple_collimation.mother = world.name + simple_collimation.size = [1 * m, 1 * m, 40 * cm] + simple_collimation.color = [0.3, 0.1, 0.8, 1] + + W_leaf = sim.add_volume("Box", "W_leaf") + W_leaf.mother = simple_collimation.name + W_leaf.size = [1 * m, 1 * m, 0.5 * cm] + W_leaf.material = "Tungsten" + leaf_translation = [] + for i in range(10): + leaf_translation.append( + [ + 0, + 0, + -0.5 * simple_collimation.size[2] + + 0.5 * W_leaf.size[2] + + i * W_leaf.size[2], + ] + ) + W_leaf.translation = leaf_translation + W_leaf.color = [0.8, 0.2, 0.1, 1] + + ######## pseudo_transportation ACTOR ######### + nb_split = 5 + pseudo_transportation_actor = sim.add_actor( + "ComptPseudoTransportationActor", "pseudo_transportation_actor" + ) + pseudo_transportation_actor.mother = simple_collimation.name + pseudo_transportation_actor.splitting_factor = nb_split + pseudo_transportation_actor.relative_min_weight_of_particle = 10000 + if j == 0: + pseudo_transportation_actor.russian_roulette_for_weights = False + if j == 1: + pseudo_transportation_actor.russian_roulette_for_weights = True + list_processes_to_bias = pseudo_transportation_actor.gamma_processes + + ##### PHASE SPACE plan ######" + plan = sim.add_volume("Box", "phsp") + plan.material = "G4_Galactic" + plan.mother = world.name + plan.size = [1 * m, 1 * m, 1 * nm] + plan.color = [0.2, 1, 0.8, 1] + plan.translation = [0, 0, -20 * cm - 1 * nm] + + ####### gamma source ########### + nb_event = 20000 + source = sim.add_source("GenericSource", "source1") + source.particle = "gamma" + source.n = nb_event + source.position.type = "sphere" + source.position.radius = 1 * nm + source.direction.type = "momentum" + source.direction.momentum = [0, 0, -1] + source.energy.type = "mono" + source.energy.mono = 6 * MeV + source.position.translation = [0, 0, 18 * cm] + + ####### PHASE SPACE ACTOR ############## + + phsp_actor = sim.add_actor("PhaseSpaceActor", "PhaseSpace") + phsp_actor.mother = plan.name + phsp_actor.attributes = [ + "EventID", + "TrackCreatorProcess", + "TrackVertexPosition", + "PrePosition", + "Weight", + "KineticEnergy", + "ParentID", + "ParticleName", + ] + data_name = "test072_output_data_RR_" + str(j) + ".root" + phsp_actor.output = paths.output / data_name + + ##### MODIFIED PHYSICS LIST ############### + + s = sim.add_actor("SimulationStatisticsActor", "Stats") + s.track_types_flag = True + sim.physics_manager.physics_list_name = "G4EmStandardPhysics_option3" + ### Perhaps avoid the user to call the below boolean function ? ### + sim.physics_manager.special_physics_constructors.G4GenericBiasingPhysics = True + sim.physics_manager.processes_to_bias.gamma = list_processes_to_bias + s = f"/process/em/UseGeneralProcess false" + sim.add_g4_command_before_init(s) + + sim.physics_manager.global_production_cuts.gamma = 1 * mm + sim.physics_manager.global_production_cuts.electron = 1 * km + sim.physics_manager.global_production_cuts.positron = 1 * km + + output = sim.run(True) + + # + # # print results + stats = sim.output.get_actor("Stats") + h = sim.output.get_actor("PhaseSpace") + print(stats) + # + + f_1 = uproot.open(paths.output / "test072_output_data_RR_0.root") + f_2 = uproot.open(paths.output / "test072_output_data_RR_1.root") + + arr_no_RR = f_1["PhaseSpace"].arrays() + arr_RR = f_2["PhaseSpace"].arrays() + + is_ok = validation_test(arr_no_RR, arr_RR, nb_event, nb_split) + utility.test_ok(is_ok) diff --git a/opengate/tests/src/test072_pseudo_transportation.py b/opengate/tests/src/test072_pseudo_transportation.py new file mode 100644 index 000000000..64f8e822f --- /dev/null +++ b/opengate/tests/src/test072_pseudo_transportation.py @@ -0,0 +1,194 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import opengate as gate +import uproot +import numpy as np +import matplotlib.pyplot as plt +import matplotlib as mpl +import scipy +from scipy.spatial.transform import Rotation +from opengate.tests import utility + + +def validation_test(arr, NIST_data, nb_split, tol=0.01): + tab_ekin = NIST_data[:, 0] + mu_att = NIST_data[:, -2] + + log_ekin = np.log(tab_ekin) + log_mu = np.log(mu_att) + + f_mu = scipy.interpolate.interp1d(log_ekin, log_mu, kind="cubic") + # print(arr["TrackCreatorProcess"]) + Tracks = arr[ + (arr["TrackCreatorProcess"] == "biasWrapper(compt)") + & (arr["KineticEnergy"] > 0.1) + & (arr["ParticleName"] == "gamma") + & (arr["Weight"] > 10 ** (-20)) + ] + ekin_tracks = Tracks["KineticEnergy"] + x_vertex = Tracks["TrackVertexPosition_X"] + y_vertex = Tracks["TrackVertexPosition_Y"] + z_vertex = Tracks["TrackVertexPosition_Z"] + x = Tracks["PrePosition_X"] + y = Tracks["PrePosition_Y"] + z = Tracks["PrePosition_Z"] + weights = Tracks["Weight"] * nb_split + dist = np.sqrt((x - x_vertex) ** 2 + (y - y_vertex) ** 2 + (z - z_vertex) ** 2) + + G4_mu = -np.log(weights) / (0.1 * dist * 19.3) + X_com_mu = np.exp(f_mu(np.log(ekin_tracks))) + diff = (G4_mu - X_com_mu) / G4_mu + print( + "Median difference between mu calculated from XCOM database and from GEANT4 free flight operation:", + np.round(100 * np.median(diff), 1), + "%", + ) + return np.median(diff) < tol + + +if __name__ == "__main__": + paths = utility.get_default_test_paths( + __file__, "test072_pseudo_transportation", output_folder="test072" + ) + + # 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 + um = gate.g4_units.um + 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 + deg = gate.g4_units.deg + + # adapt world size + world = sim.world + world.size = [1.2 * m, 1.2 * m, 2 * m] + world.material = "G4_Galactic" + + ####### GEOMETRY TO IRRADIATE ############# + sim.volume_manager.material_database.add_material_weights( + "Tungsten", + ["W"], + [1], + 19.3 * gcm3, + ) + simple_collimation = sim.add_volume("Box", "colli_box") + simple_collimation.material = "G4_Galactic" + simple_collimation.mother = world.name + simple_collimation.size = [1 * m, 1 * m, 40 * cm] + simple_collimation.color = [0.3, 0.1, 0.8, 1] + + W_leaf = sim.add_volume("Box", "W_leaf") + W_leaf.mother = simple_collimation.name + W_leaf.size = [1 * m, 1 * m, 2 * cm] + W_leaf.material = "Tungsten" + leaf_translation = [] + for i in range(10): + leaf_translation.append( + [ + 0, + 0, + -0.5 * simple_collimation.size[2] + + 0.5 * W_leaf.size[2] + + i * W_leaf.size[2], + ] + ) + print(leaf_translation) + W_leaf.translation = leaf_translation + W_leaf.color = [0.8, 0.2, 0.1, 1] + + ######## pseudo_transportation ACTOR ######### + nb_split = 5 + pseudo_transportation_actor = sim.add_actor( + "ComptPseudoTransportationActor", "pseudo_transportation_actor" + ) + pseudo_transportation_actor.mother = simple_collimation.name + pseudo_transportation_actor.splitting_factor = nb_split + pseudo_transportation_actor.relative_min_weight_of_particle = np.inf + list_processes_to_bias = pseudo_transportation_actor.gamma_processes + + ##### PHASE SPACE plan ######" + plan = sim.add_volume("Box", "phsp") + plan.material = "G4_Galactic" + plan.mother = world.name + plan.size = [1 * m, 1 * m, 1 * nm] + plan.color = [0.2, 1, 0.8, 1] + plan.translation = [0, 0, -20 * cm - 1 * nm] + + ####### gamma source ########### + source = sim.add_source("GenericSource", "source1") + source.particle = "gamma" + source.n = 1000 + source.position.type = "sphere" + source.position.radius = 1 * nm + source.direction.type = "momentum" + source.direction.momentum = [0, 0, -1] + source.energy.type = "mono" + source.energy.mono = 6 * MeV + source.position.translation = [0, 0, 18 * cm] + + ####### PHASE SPACE ACTOR ############## + + phsp_actor = sim.add_actor("PhaseSpaceActor", "PhaseSpace") + phsp_actor.mother = plan.name + phsp_actor.attributes = [ + "EventID", + "TrackCreatorProcess", + "TrackVertexPosition", + "PrePosition", + "Weight", + "KineticEnergy", + "ParentID", + "ParticleName", + ] + + phsp_actor.output = paths.output / "test072_output_data.root" + + ##### MODIFIED PHYSICS LIST ############### + + s = sim.add_actor("SimulationStatisticsActor", "Stats") + s.track_types_flag = True + sim.physics_manager.physics_list_name = "G4EmStandardPhysics_option3" + ### Perhaps avoid the user to call the below boolean function ? ### + sim.physics_manager.special_physics_constructors.G4GenericBiasingPhysics = True + sim.physics_manager.processes_to_bias.gamma = list_processes_to_bias + s = f"/process/em/UseGeneralProcess false" + sim.add_g4_command_before_init(s) + + sim.physics_manager.global_production_cuts.gamma = 1 * mm + sim.physics_manager.global_production_cuts.electron = 1 * km + sim.physics_manager.global_production_cuts.positron = 1 * km + + output = sim.run() + + # + # # print results + stats = sim.output.get_actor("Stats") + h = sim.output.get_actor("PhaseSpace") + print(stats) + # + f_phsp = uproot.open(paths.output / "test072_output_data.root") + data_NIST_W = np.loadtxt(paths.data / "NIST_W.txt", delimiter="|") + arr = f_phsp["PhaseSpace"].arrays() + # + is_ok = validation_test(arr, data_NIST_W, nb_split) + utility.test_ok(is_ok) diff --git a/opengate/tests/src/test074_kill_non_interacting_particles.py b/opengate/tests/src/test074_kill_non_interacting_particles.py new file mode 100644 index 000000000..36616f210 --- /dev/null +++ b/opengate/tests/src/test074_kill_non_interacting_particles.py @@ -0,0 +1,212 @@ +#!/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 +from anytree import Node, RenderTree +import uproot + + +def test074_test(entry_data, exit_data_1, exit_data_2): + liste_ekin = [] + liste_evtID = [] + liste_trackID = [] + evt_ID_entry_data = entry_data["EventID"] + j = 0 + i = 0 + while i < len(evt_ID_entry_data): + if ( + j < len(exit_data_1["EventID"]) + and evt_ID_entry_data[i] == exit_data_1["EventID"][j] + ): + TID_entry = entry_data["TrackID"][i] + TID_exit = exit_data_1["TrackID"][j] + Ekin_entry = entry_data["KineticEnergy"][i] + Ekin_exit = exit_data_1["KineticEnergy"][j] + + if (TID_entry == TID_exit) and (Ekin_exit == Ekin_entry): + liste_ekin.append(exit_data_1["KineticEnergy"][j]) + liste_evtID.append(exit_data_1["EventID"][j]) + liste_trackID.append(exit_data_1["TrackID"][j]) + if (j < len(exit_data_1["EventID"]) - 1) and ( + exit_data_1["EventID"][j] == exit_data_1["EventID"][j + 1] + ): + i = i - 1 + j += 1 + i += 1 + liste_ekin = np.asarray(liste_ekin) + print("Number of tracks to kill =", len(liste_ekin)) + print( + "Number of killed tracks =", + (len(exit_data_1["EventID"]) - len(exit_data_2["EventID"])), + ) + + return len(liste_ekin) == ( + len(exit_data_1["EventID"]) - len(exit_data_2["EventID"]) + ) + + +if __name__ == "__main__": + paths = utility.get_default_test_paths(__file__) + output_path = paths.output + + print(output_path) + # 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"] + + sim.volume_manager.material_database.add_material_weights( + "Tungsten", + ["W"], + [1], + 19.3 * gcm3, + ) + + # adapt world size + world = sim.world + world.size = [1 * m, 1 * m, 1 * m] + world.material = "G4_AIR" + + big_box = sim.add_volume("Box", "big_box") + big_box.mother = world.name + big_box.material = "G4_AIR" + big_box.size = [0.8 * m, 0.8 * m, 0.8 * m] + + actor_box = sim.add_volume("Box", "actor_box") + actor_box.mother = big_box.name + actor_box.material = "G4_AIR" + actor_box.size = [0.6 * m, 0.6 * m, 0.6 * m] + actor_box.translation = [0, 0, -0.1 * m] + + source = sim.add_source("GenericSource", "photon_source") + source.particle = "gamma" + source.position.type = "box" + source.mother = world.name + source.position.size = [6 * cm, 6 * cm, 6 * cm] + source.position.translation = [0, 0, 0.3 * m] + source.direction.type = "momentum" + source.force_rotation = True + # source1.direction.focus_point = [0*cm, 0*cm, -5 *cm] + source.direction.momentum = [0, 0, -1] + source.energy.type = "mono" + source.energy.mono = 6 * MeV + source.n = 1000 + + tungsten_leaves = sim.add_volume("Box", "tungsten_leaves") + tungsten_leaves.mother = actor_box + tungsten_leaves.size = [0.6 * m, 0.6 * m, 0.3 * cm] + tungsten_leaves.material = "Tungsten" + liste_translation_W = [] + for i in range(7): + liste_translation_W.append([0, 0, 0.25 * m - i * 6 * cm]) + tungsten_leaves.translation = liste_translation_W + tungsten_leaves.color = [0.9, 0.0, 0.4, 0.8] + + kill_No_int_act = sim.add_actor("KillNonInteractingParticleActor", "killact") + kill_No_int_act.mother = actor_box.name + + entry_phase_space = sim.add_volume("Box", "entry_phase_space") + entry_phase_space.mother = big_box + entry_phase_space.size = [0.8 * m, 0.8 * m, 1 * nm] + entry_phase_space.material = "G4_AIR" + entry_phase_space.translation = [0, 0, 0.21 * m] + entry_phase_space.color = [0.5, 0.9, 0.3, 1] + + exit_phase_space_1 = sim.add_volume("Box", "exit_phase_space_1") + exit_phase_space_1.mother = actor_box + exit_phase_space_1.size = [0.6 * m, 0.6 * m, 1 * nm] + exit_phase_space_1.material = "G4_AIR" + exit_phase_space_1.translation = [0, 0, -0.3 * m + 1 * nm] + exit_phase_space_1.color = [0.5, 0.9, 0.3, 1] + + exit_phase_space_2 = sim.add_volume("Box", "exit_phase_space_2") + exit_phase_space_2.mother = world.name + exit_phase_space_2.size = [0.6 * m, 0.6 * m, 1 * nm] + exit_phase_space_2.material = "G4_AIR" + exit_phase_space_2.translation = [0, 0, -0.4 * m - 1 * nm] + exit_phase_space_2.color = [0.5, 0.9, 0.3, 1] + + # print(sim.volume_manager.dump_volume_tree()) + liste_phase_space_name = [ + entry_phase_space.name, + exit_phase_space_1.name, + exit_phase_space_2.name, + ] + for name in liste_phase_space_name: + + phsp = sim.add_actor("PhaseSpaceActor", "PhaseSpace_" + name) + phsp.mother = name + phsp.attributes = ["EventID", "TrackID", "KineticEnergy"] + name_phsp = "test074_" + name + ".root" + phsp.output = output_path / name_phsp + + 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 + + s = sim.add_actor("SimulationStatisticsActor", "Stats") + s.track_types_flag = True + + # go ! + sim.run() + output = sim.output + stats = sim.output.get_actor("Stats") + print(stats) + + entry_phsp = uproot.open( + str(output_path) + + "/test074_" + + liste_phase_space_name[0] + + ".root" + + ":PhaseSpace_" + + liste_phase_space_name[0] + ) + exit_phase_space_1 = uproot.open( + str(output_path) + + "/test074_" + + liste_phase_space_name[1] + + ".root" + + ":PhaseSpace_" + + liste_phase_space_name[1] + ) + exit_phase_space_2 = uproot.open( + str(output_path) + + "/test074_" + + liste_phase_space_name[2] + + ".root" + + ":PhaseSpace_" + + liste_phase_space_name[2] + ) + + df_entry = entry_phsp.arrays() + df_exit_1 = exit_phase_space_1.arrays() + df_exit_2 = exit_phase_space_2.arrays() + + is_ok = test074_test(df_entry, df_exit_1, df_exit_2) + + utility.test_ok(is_ok) diff --git a/opengate/tests/src/test075_geometrical_splitting_volume_surface.py b/opengate/tests/src/test075_geometrical_splitting_volume_surface.py new file mode 100644 index 000000000..58af57f28 --- /dev/null +++ b/opengate/tests/src/test075_geometrical_splitting_volume_surface.py @@ -0,0 +1,185 @@ +#!/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 +from anytree import Node, RenderTree +import uproot + + +def test075(entry_data, exit_data, splitting_factor): + splitted_particle_data_entry = entry_data[ + entry_data["TrackCreatorProcess"] == "none" + ] + + splitted_particle_data_exit = exit_data[exit_data["TrackCreatorProcess"] == "none"] + + array_weight_1 = splitted_particle_data_entry["Weight"] + array_weight_2 = splitted_particle_data_exit["Weight"] + + if (np.round(np.sum(array_weight_1), 3) == 1) and len( + array_weight_1 + ) == splitting_factor: + if (np.round(np.sum(array_weight_2), 3) == 1) and len( + array_weight_2 + ) == splitting_factor: + return True + else: + return False + else: + return False + + +if __name__ == "__main__": + paths = utility.get_default_test_paths(__file__) + output_path = paths.output + + print(output_path) + # 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" + + big_box = sim.add_volume("Box", "big_box") + big_box.mother = world.name + big_box.material = "G4_Galactic" + big_box.size = [0.8 * m, 0.8 * m, 0.8 * m] + + actor_box = sim.add_volume("Box", "actor_box") + actor_box.mother = big_box.name + actor_box.material = "G4_Galactic" + actor_box.size = [0.6 * m, 0.6 * m, 0.6 * m] + actor_box.translation = [0, 0, -0.1 * m] + + source_1 = sim.add_source("GenericSource", "elec_source_1") + source_1.particle = "e-" + source_1.position.type = "box" + source_1.mother = big_box.name + source_1.position.size = [1 * cm, 1 * cm, 1 * cm] + source_1.position.translation = [0, 0.35 * m, 0] + source_1.direction.type = "momentum" + source_1.direction.momentum = [0, -1, 0] + source_1.energy.type = "mono" + source_1.energy.mono = 10 * MeV + source_1.n = 1 + + source_2 = sim.add_source("GenericSource", "elec_source_2") + source_2.particle = "e-" + source_2.position.type = "box" + source_2.mother = big_box.name + source_2.position.size = [1 * cm, 1 * cm, 1 * cm] + source_2.position.translation = [0, 0, -0.39 * m] + source_2.direction.type = "momentum" + source_2.direction.momentum = [0, 0, -1] + source_2.energy.type = "mono" + source_2.energy.mono = 10 * MeV + source_2.n = 1 + + geom_splitting = sim.add_actor("SurfaceSplittingActor", "splitting_act") + geom_splitting.mother = actor_box.name + geom_splitting.splitting_factor = 10 + geom_splitting.weight_threshold = 1 + geom_splitting.split_entering_particles = True + geom_splitting.split_exiting_particles = True + + entry_phase_space = sim.add_volume("Box", "entry_phase_space") + entry_phase_space.mother = actor_box + entry_phase_space.size = [0.6 * m, 1 * nm, 0.6 * m] + entry_phase_space.material = "G4_Galactic" + entry_phase_space.translation = [0, 0.3 * m - 0.5 * nm, 0] + entry_phase_space.color = [0.5, 0.9, 0.3, 1] + + exit_phase_space = sim.add_volume("Box", "exit_phase_space") + exit_phase_space.mother = world.name + exit_phase_space.size = [0.6 * m, 0.6 * m, 1 * nm] + exit_phase_space.material = "G4_Galactic" + exit_phase_space.translation = [0, 0, -0.4 * m - 1 * nm] + exit_phase_space.color = [0.5, 0.9, 0.3, 1] + + # # print(sim.volume_manager.dump_volume_tree()) + liste_phase_space_name = [ + entry_phase_space.name, + exit_phase_space.name, + ] + for name in liste_phase_space_name: + + phsp = sim.add_actor("PhaseSpaceActor", "PhaseSpace_" + name) + phsp.mother = name + phsp.attributes = [ + "EventID", + "TrackID", + "Weight", + "PDGCode", + "TrackCreatorProcess", + ] + name_phsp = "test075_" + name + ".root" + phsp.output = output_path / name_phsp + + 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 + + s = sim.add_actor("SimulationStatisticsActor", "Stats") + s.track_types_flag = True + + # go ! + sim.run() + output = sim.output + stats = sim.output.get_actor("Stats") + print(stats) + + # + entry_phsp = uproot.open( + str(output_path) + + "/test075_" + + liste_phase_space_name[0] + + ".root" + + ":PhaseSpace_" + + liste_phase_space_name[0] + ) + exit_phase_space = uproot.open( + str(output_path) + + "/test075_" + + liste_phase_space_name[1] + + ".root" + + ":PhaseSpace_" + + liste_phase_space_name[1] + ) + + # + df_entry = entry_phsp.arrays() + df_exit = exit_phase_space.arrays() + # + + is_ok = test075(df_entry, df_exit, geom_splitting.splitting_factor) + + utility.test_ok(is_ok) diff --git a/opengate/tests/src/test077_kill_interacting_particles.py b/opengate/tests/src/test077_kill_interacting_particles.py new file mode 100644 index 000000000..fc342cdf5 --- /dev/null +++ b/opengate/tests/src/test077_kill_interacting_particles.py @@ -0,0 +1,213 @@ +#!/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 +from anytree import Node, RenderTree +import uproot + + +def test077_test(entry_data, exit_data_1, exit_data_2): + + liste_ekin = [] + liste_evtID = [] + liste_trackID = [] + evt_ID_entry_data = entry_data["EventID"] + j = 0 + i = 0 + while i < len(evt_ID_entry_data): + if ( + j < len(exit_data_1["EventID"]) + and evt_ID_entry_data[i] == exit_data_1["EventID"][j] + ): + TID_entry = entry_data["TrackID"][i] + TID_exit = exit_data_1["TrackID"][j] + Ekin_entry = entry_data["KineticEnergy"][i] + Ekin_exit = exit_data_1["KineticEnergy"][j] + + if (TID_entry != TID_exit) or (Ekin_exit != Ekin_entry): + liste_ekin.append(exit_data_1["KineticEnergy"][j]) + liste_evtID.append(exit_data_1["EventID"][j]) + liste_trackID.append(exit_data_1["TrackID"][j]) + if (j < len(exit_data_1["EventID"]) - 1) and ( + exit_data_1["EventID"][j] == exit_data_1["EventID"][j + 1] + ): + i = i - 1 + j += 1 + i += 1 + liste_ekin = np.asarray(liste_ekin) + print("Number of tracks to kill =", len(liste_ekin)) + print( + "Number of killed tracks =", + (len(exit_data_1["EventID"]) - len(exit_data_2["EventID"])), + ) + + return len(liste_ekin) == ( + len(exit_data_1["EventID"]) - len(exit_data_2["EventID"]) + ) + + +if __name__ == "__main__": + paths = utility.get_default_test_paths(__file__) + output_path = paths.output + + print(output_path) + # 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"] + + sim.volume_manager.material_database.add_material_weights( + "Tungsten", + ["W"], + [1], + 19.3 * gcm3, + ) + + # adapt world size + world = sim.world + world.size = [1 * m, 1 * m, 1 * m] + world.material = "G4_AIR" + + big_box = sim.add_volume("Box", "big_box") + big_box.mother = world.name + big_box.material = "G4_AIR" + big_box.size = [0.8 * m, 0.8 * m, 0.8 * m] + + actor_box = sim.add_volume("Box", "actor_box") + actor_box.mother = big_box.name + actor_box.material = "G4_AIR" + actor_box.size = [0.6 * m, 0.6 * m, 0.6 * m] + actor_box.translation = [0, 0, -0.1 * m] + + source = sim.add_source("GenericSource", "photon_source") + source.particle = "gamma" + source.position.type = "box" + source.mother = world.name + source.position.size = [6 * cm, 6 * cm, 6 * cm] + source.position.translation = [0, 0, 0.3 * m] + source.direction.type = "momentum" + source.force_rotation = True + # source1.direction.focus_point = [0*cm, 0*cm, -5 *cm] + source.direction.momentum = [0, 0, -1] + source.energy.type = "mono" + source.energy.mono = 6 * MeV + source.n = 1000 + + tungsten_leaves = sim.add_volume("Box", "tungsten_leaves") + tungsten_leaves.mother = actor_box + tungsten_leaves.size = [0.6 * m, 0.6 * m, 0.3 * cm] + tungsten_leaves.material = "Tungsten" + liste_translation_W = [] + for i in range(7): + liste_translation_W.append([0, 0, 0.25 * m - i * 6 * cm]) + tungsten_leaves.translation = liste_translation_W + tungsten_leaves.color = [0.9, 0.0, 0.4, 0.8] + + kill_int_act = sim.add_actor("KillInteractingParticleActor", "killact") + kill_int_act.mother = actor_box.name + + entry_phase_space = sim.add_volume("Box", "entry_phase_space") + entry_phase_space.mother = big_box + entry_phase_space.size = [0.8 * m, 0.8 * m, 1 * nm] + entry_phase_space.material = "G4_AIR" + entry_phase_space.translation = [0, 0, 0.21 * m] + entry_phase_space.color = [0.5, 0.9, 0.3, 1] + + exit_phase_space_1 = sim.add_volume("Box", "exit_phase_space_1") + exit_phase_space_1.mother = actor_box + exit_phase_space_1.size = [0.6 * m, 0.6 * m, 1 * nm] + exit_phase_space_1.material = "G4_AIR" + exit_phase_space_1.translation = [0, 0, -0.3 * m + 1 * nm] + exit_phase_space_1.color = [0.5, 0.9, 0.3, 1] + + exit_phase_space_2 = sim.add_volume("Box", "exit_phase_space_2") + exit_phase_space_2.mother = world.name + exit_phase_space_2.size = [0.6 * m, 0.6 * m, 1 * nm] + exit_phase_space_2.material = "G4_AIR" + exit_phase_space_2.translation = [0, 0, -0.4 * m - 1 * nm] + exit_phase_space_2.color = [0.5, 0.9, 0.3, 1] + + # print(sim.volume_manager.dump_volume_tree()) + liste_phase_space_name = [ + entry_phase_space.name, + exit_phase_space_1.name, + exit_phase_space_2.name, + ] + for name in liste_phase_space_name: + + phsp = sim.add_actor("PhaseSpaceActor", "PhaseSpace_" + name) + phsp.mother = name + phsp.attributes = ["EventID", "TrackID", "KineticEnergy"] + name_phsp = "test077_" + name + ".root" + phsp.output = output_path / name_phsp + + 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 + + s = sim.add_actor("SimulationStatisticsActor", "Stats") + s.track_types_flag = True + + # go ! + sim.run() + output = sim.output + stats = sim.output.get_actor("Stats") + print(stats) + + entry_phsp = uproot.open( + str(output_path) + + "/test077_" + + liste_phase_space_name[0] + + ".root" + + ":PhaseSpace_" + + liste_phase_space_name[0] + ) + exit_phase_space_1 = uproot.open( + str(output_path) + + "/test077_" + + liste_phase_space_name[1] + + ".root" + + ":PhaseSpace_" + + liste_phase_space_name[1] + ) + exit_phase_space_2 = uproot.open( + str(output_path) + + "/test077_" + + liste_phase_space_name[2] + + ".root" + + ":PhaseSpace_" + + liste_phase_space_name[2] + ) + + df_entry = entry_phsp.arrays() + df_exit_1 = exit_phase_space_1.arrays() + df_exit_2 = exit_phase_space_2.arrays() + + is_ok = test077_test(df_entry, df_exit_1, df_exit_2) + + utility.test_ok(is_ok) diff --git a/opengate/tests/src/test082_kill_according_all_processes.py b/opengate/tests/src/test082_kill_according_all_processes.py index 145555e77..15a9a0749 100755 --- a/opengate/tests/src/test082_kill_according_all_processes.py +++ b/opengate/tests/src/test082_kill_according_all_processes.py @@ -10,6 +10,7 @@ def test082_test(df): + df = df[df["PDGCode"] == 22] nb_event = len(df["ParentID"]) nb_event_to_interest = len(df["ParentID"][df["ParentID"] == 0]) @@ -133,6 +134,7 @@ def test082_test(df): sim.run() print(kill_proc_act) # + phsp = uproot.open( str(output_path) + "/test082_PhaseSpace_all.root" + ":PhaseSpace" ) diff --git a/opengate/tests/src/test082_kill_according_processes.py b/opengate/tests/src/test082_kill_according_processes.py index 974cc1955..7bc4d45e0 100755 --- a/opengate/tests/src/test082_kill_according_processes.py +++ b/opengate/tests/src/test082_kill_according_processes.py @@ -115,6 +115,7 @@ def test082_test(df): "PreDirection", "PDGCode", ] + name_phsp = "test082_" + phsp.name + ".root" phsp.output_filename = name_phsp @@ -135,8 +136,8 @@ def test082_test(df): # # go ! sim.run() # - phsp = uproot.open(str(output_path) + "/test082_PhaseSpace.root" + ":PhaseSpace") + phsp = uproot.open(str(output_path) + "/test082_PhaseSpace.root" + ":PhaseSpace") df = phsp.arrays() is_ok = test082_test(df) # diff --git a/opengate/tests/src/test082_kill_non_interacting_particles.py b/opengate/tests/src/test082_kill_non_interacting_particles.py new file mode 100644 index 000000000..73c439179 --- /dev/null +++ b/opengate/tests/src/test082_kill_non_interacting_particles.py @@ -0,0 +1,211 @@ +#!/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 +from anytree import Node, RenderTree +import uproot + + +def test082_test(entry_data, exit_data_1, exit_data_2): + liste_ekin = [] + liste_evtID = [] + liste_trackID = [] + evt_ID_entry_data = entry_data["EventID"] + j = 0 + i = 0 + while i < len(evt_ID_entry_data): + if ( + j < len(exit_data_1["EventID"]) + and evt_ID_entry_data[i] == exit_data_1["EventID"][j] + ): + TID_entry = entry_data["TrackID"][i] + TID_exit = exit_data_1["TrackID"][j] + Ekin_entry = entry_data["KineticEnergy"][i] + Ekin_exit = exit_data_1["KineticEnergy"][j] + + if (TID_entry == TID_exit) and (Ekin_exit == Ekin_entry): + liste_ekin.append(exit_data_1["KineticEnergy"][j]) + liste_evtID.append(exit_data_1["EventID"][j]) + liste_trackID.append(exit_data_1["TrackID"][j]) + if (j < len(exit_data_1["EventID"]) - 1) and ( + exit_data_1["EventID"][j] == exit_data_1["EventID"][j + 1] + ): + i = i - 1 + j += 1 + i += 1 + liste_ekin = np.asarray(liste_ekin) + print("Number of tracks to kill =", len(liste_ekin)) + print( + "Number of killed tracks =", + (len(exit_data_1["EventID"]) - len(exit_data_2["EventID"])), + ) + + return len(liste_ekin) == ( + len(exit_data_1["EventID"]) - len(exit_data_2["EventID"]) + ) + + +if __name__ == "__main__": + paths = utility.get_default_test_paths(__file__) + output_path = paths.output + # create the simulation + sim = gate.Simulation() + sim.output_dir = output_path + + # 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"] + + sim.volume_manager.material_database.add_material_weights( + "Tungsten", + ["W"], + [1], + 19.3 * gcm3, + ) + + # adapt world size + world = sim.world + world.size = [1 * m, 1 * m, 1 * m] + world.material = "G4_Galactic" + + big_box = sim.add_volume("Box", "big_box") + big_box.mother = world.name + big_box.material = "G4_Galactic" + big_box.size = [0.8 * m, 0.8 * m, 0.8 * m] + + actor_box = sim.add_volume("Box", "actor_box") + actor_box.mother = big_box.name + actor_box.material = "G4_Galactic" + actor_box.size = [0.6 * m, 0.6 * m, 0.6 * m] + actor_box.translation = [0, 0, -0.1 * m] + + source = sim.add_source("GenericSource", "photon_source") + source.particle = "gamma" + source.position.type = "box" + source.mother = world.name + source.position.size = [6 * cm, 6 * cm, 4 * cm] + source.position.translation = [0, 0, 0.205 * m] + source.direction.type = "momentum" + source.direction_relative_to_attached_volume = True + # source1.direction.focus_point = [0*cm, 0*cm, -5 *cm] + source.direction.momentum = [0, 0, -1] + source.energy.type = "mono" + source.energy.mono = 6 * MeV + source.n = 5000 + + tungsten_leaves = sim.add_volume("Box", "tungsten_leaves") + tungsten_leaves.mother = actor_box + tungsten_leaves.size = [0.6 * m, 0.6 * m, 0.3 * cm] + tungsten_leaves.material = "Tungsten" + liste_translation_W = [] + for i in range(7): + liste_translation_W.append([0, 0, 0.25 * m - i * 6 * cm]) + tungsten_leaves.translation = liste_translation_W + tungsten_leaves.color = [0.9, 0.0, 0.4, 0.8] + + kill_No_int_act = sim.add_actor("KillNonInteractingParticleActor", "killact") + kill_No_int_act.attached_to = actor_box.name + + entry_phase_space = sim.add_volume("Box", "entry_phase_space") + entry_phase_space.mother = actor_box.name + entry_phase_space.size = [0.6 * m, 0.6 * m, 1 * nm] + entry_phase_space.material = "G4_Galactic" + entry_phase_space.translation = [0, 0, 0.255 * m] + entry_phase_space.color = [0.5, 0.9, 0.3, 1] + + exit_phase_space_1 = sim.add_volume("Box", "exit_phase_space_1") + exit_phase_space_1.mother = actor_box + exit_phase_space_1.size = [0.6 * m, 0.6 * m, 1 * nm] + exit_phase_space_1.material = "G4_Galactic" + exit_phase_space_1.translation = [0, 0, -0.3 * m + 1 * nm] + exit_phase_space_1.color = [0.5, 0.9, 0.3, 1] + + exit_phase_space_2 = sim.add_volume("Box", "exit_phase_space_2") + exit_phase_space_2.mother = world.name + exit_phase_space_2.size = [0.6 * m, 0.6 * m, 1 * nm] + exit_phase_space_2.material = "G4_Galactic" + exit_phase_space_2.translation = [0, 0, -0.4 * m - 1 * nm] + exit_phase_space_2.color = [0.5, 0.9, 0.3, 1] + + # print(sim.volume_manager.dump_volume_tree()) + liste_phase_space_name = [ + entry_phase_space.name, + exit_phase_space_1.name, + exit_phase_space_2.name, + ] + for name in liste_phase_space_name: + print(name) + phsp = sim.add_actor("PhaseSpaceActor", "PhaseSpace_" + name) + phsp.attached_to = name + phsp.attributes = ["EventID", "TrackID", "KineticEnergy"] + name_phsp = "test082_" + name + ".root" + phsp.output_filename = name_phsp + + 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 + + s = sim.add_actor("SimulationStatisticsActor", "Stats") + s.track_types_flag = True + + # go ! + sim.run() + print(s) + print(output_path) + entry_phsp = uproot.open( + str(output_path) + + "/test082_" + + liste_phase_space_name[0] + + ".root" + + ":PhaseSpace_" + + liste_phase_space_name[0] + ) + exit_phase_space_1 = uproot.open( + str(output_path) + + "/test082_" + + liste_phase_space_name[1] + + ".root" + + ":PhaseSpace_" + + liste_phase_space_name[1] + ) + exit_phase_space_2 = uproot.open( + str(output_path) + + "/test082_" + + liste_phase_space_name[2] + + ".root" + + ":PhaseSpace_" + + liste_phase_space_name[2] + ) + + df_entry = entry_phsp.arrays() + df_exit_1 = exit_phase_space_1.arrays() + df_exit_2 = exit_phase_space_2.arrays() + + is_ok = test082_test(df_entry, df_exit_1, df_exit_2) + + print(kill_No_int_act.user_output.kill_non_interacting_particles) + + utility.test_ok(is_ok) diff --git a/opengate/tests/src/test082_last_vertex_splittting.py b/opengate/tests/src/test082_last_vertex_splittting.py new file mode 100755 index 000000000..6a40856a6 --- /dev/null +++ b/opengate/tests/src/test082_last_vertex_splittting.py @@ -0,0 +1,223 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import opengate as gate +import uproot +import numpy as np +from scipy.spatial.transform import Rotation +from opengate.tests import utility + + +def validation_test(arr_ref, arr_data, nb_split): + arr_ref = arr_ref[arr_ref["ParticleName"] == "gamma"] + arr_data = arr_data[arr_data["ParticleName"] == "gamma"] + + weight_data = np.round(np.mean(arr_data["Weight"]), 4) + bool_weight = False + print( + "Weight mean is equal to", weight_data, "and need to be equal to", 1 / nb_split + ) + if weight_data == 1 / nb_split: + bool_weight = True + + bool_events = False + sigma = np.sqrt((len(arr_data["KineticEnergy"]) / nb_split)) * nb_split + nb_events_ref = len(arr_ref["KineticEnergy"]) + nb_events_data = len(arr_data["KineticEnergy"]) + print("Reference counts number:", nb_events_ref) + print("Biased counts number:", nb_events_data) + if nb_events_data - 4 * sigma <= nb_events_ref <= nb_events_data + 4 * sigma: + bool_events = True + + keys = [ + "KineticEnergy", + "PreDirection_X", + "PreDirection_Y", + "PreDirection_Z", + "PrePosition_X", + "PrePosition_Y", + ] + + bool_distrib = True + for key in keys: + ref = arr_ref[key] + data = arr_data[key] + + mean_ref = np.mean(ref) + mean_data = np.mean(data) + + std_dev_ref = np.std(ref, ddof=1) + std_dev_data = np.std(data, ddof=1) + + std_err_ref = std_dev_ref / np.sqrt(len(ref)) + std_err_data = std_dev_data / (nb_split * np.sqrt(len(data) / nb_split)) + + print( + key, + "mean ref value:", + np.round(mean_ref, 3), + "+-", + np.round(std_err_ref, 3), + ) + print( + key, + "mean data value:", + np.round(mean_data, 3), + "+-", + np.round(std_err_data, 3), + ) + + if ( + mean_data - 4 * np.sqrt(std_err_data**2 + std_err_ref**2) > mean_ref + ) or (mean_data + 4 * np.sqrt(std_err_data**2 + std_err_ref**2) < mean_ref): + bool_distrib = False + + if bool_distrib and bool_events and bool_weight: + return True + else: + return False + + +if __name__ == "__main__": + for i in range(2): + if i == 0: + bias = False + else: + bias = True + paths = utility.get_default_test_paths( + __file__, "test084_last_vertex_splitting", output_folder="test084" + ) + + # create the simulation + sim = gate.Simulation() + + # main options + ui = sim.user_info + ui.g4_verbose = False + ui.visu = False + ui.visu_type = "vrml" + ui.check_volumes_overlap = False + ui.number_of_threads = 1 + ui.random_seed = 123456789 + + # units + m = gate.g4_units.m + km = gate.g4_units.km + mm = gate.g4_units.mm + um = gate.g4_units.um + 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 + deg = gate.g4_units.deg + gcm3 = gate.g4_units.g / gate.g4_units.cm3 + + # adapt world size + world = sim.world + world.size = [0.25 * m, 0.25 * m, 0.25 * m] + world.material = "G4_Galactic" + + ####### GEOMETRY TO IRRADIATE ############# + sim.volume_manager.material_database.add_material_weights( + "Tungsten", + ["W"], + [1], + 19.3 * gcm3, + ) + + W_tubs = sim.add_volume("Tubs", "W_box") + W_tubs.material = "Tungsten" + W_tubs.mother = world.name + + W_tubs.rmin = 0 + W_tubs.rmax = 0.4 * cm + W_tubs.dz = 0.05 * m + W_tubs.color = [0.8, 0.2, 0.1, 1] + angle_x = 45 + angle_y = 70 + angle_z = 80 + + rotation = Rotation.from_euler( + "xyz", [angle_y, angle_y, angle_z], degrees=True + ).as_matrix() + W_tubs.rotation = rotation + + if bias: + ###### Last vertex Splitting ACTOR ######### + nb_split = 10 + vertex_splitting_actor = sim.add_actor( + "LastVertexInteractionSplittingActor", "vertexSplittingW" + ) + vertex_splitting_actor.attached_to = W_tubs.name + vertex_splitting_actor.splitting_factor = nb_split + vertex_splitting_actor.angular_kill = True + vertex_splitting_actor.vector_director = [0, 0, -1] + vertex_splitting_actor.max_theta = 90 * deg + vertex_splitting_actor.batch_size = 10 * nb_split + vertex_splitting_actor.nb_of_max_batch_per_event = 500 + + plan = sim.add_volume("Box", "plan_phsp") + plan.material = "G4_Galactic" + plan.size = [5 * cm, 5 * cm, 1 * nm] + plan.translation = [0, 0, -1 * cm] + + ####### gamma source ########### + source = sim.add_source("GenericSource", "source1") + source.particle = "gamma" + source.n = 100000 + if bias: + source.n = source.n / nb_split + + source.position.type = "sphere" + source.position.radius = 1 * nm + source.direction.type = "momentum" + # source.direction.momentum = [0,0,-1] + source.direction.momentum = np.dot(rotation, np.array([0, 0, -1])) + source.energy.type = "mono" + source.energy.mono = 4 * MeV + # + ###### LastVertexSource ############# + if bias: + source_0 = sim.add_source("LastVertexSource", "source_vertex") + source_0.n = 1 + + ####### PHASE SPACE ACTOR ############## + sim.output_dir = paths.output + phsp_actor = sim.add_actor("PhaseSpaceActor", "PhaseSpace") + phsp_actor.attached_to = plan.name + phsp_actor.attributes = [ + "EventID", + "TrackID", + "Weight", + "ParticleName", + "KineticEnergy", + "PreDirection", + "PrePosition", + "TrackCreatorProcess", + ] + if bias: + phsp_actor.output_filename = "test084_output_data_last_vertex_biased.root" + else: + phsp_actor.output_filename = "test084_output_data_last_vertex_ref.root" + + s = sim.add_actor("SimulationStatisticsActor", "Stats") + s.track_types_flag = True + sim.physics_manager.physics_list_name = "G4EmStandardPhysics_option3" + s = f"/process/em/UseGeneralProcess false" + sim.g4_commands_before_init.append(s) + + sim.physics_manager.global_production_cuts.gamma = 1 * mm + sim.physics_manager.global_production_cuts.electron = 1000 * km + sim.physics_manager.global_production_cuts.positron = 1000 * km + + output = sim.run(True) + print(s) + + f_data = uproot.open(paths.output / "test084_output_data_last_vertex_biased.root") + f_ref_data = uproot.open(paths.output / "test084_output_data_last_vertex_ref.root") + arr_data = f_data["PhaseSpace"].arrays() + arr_ref_data = f_ref_data["PhaseSpace"].arrays() + # # + is_ok = validation_test(arr_ref_data, arr_data, nb_split) + utility.test_ok(is_ok) diff --git a/opengate/tests/src/test082_last_vertex_splittting_angular_kill.py b/opengate/tests/src/test082_last_vertex_splittting_angular_kill.py new file mode 100755 index 000000000..d1fa56cf2 --- /dev/null +++ b/opengate/tests/src/test082_last_vertex_splittting_angular_kill.py @@ -0,0 +1,171 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import opengate as gate +import uproot +import numpy as np +from scipy.spatial.transform import Rotation +from opengate.tests import utility + + +def validation_test(arr_data, vector_director, max_theta): + arr_data = arr_data[arr_data["ParticleName"] == "gamma"] + qt_mvt_data = arr_data[["PreDirection_X", "PreDirection_Y", "PreDirection_Z"]] + mom_data = np.zeros((len(qt_mvt_data["PreDirection_X"]), 3)) + mom_data[:, 0] += np.array(qt_mvt_data["PreDirection_X"]) + mom_data[:, 1] += np.array(qt_mvt_data["PreDirection_Y"]) + mom_data[:, 2] += np.array(qt_mvt_data["PreDirection_Z"]) + + l_theta = np.zeros(len(mom_data[:, 0])) + + for i in range(len(l_theta)): + theta = np.arccos(mom_data[i].dot(vector_director)) + l_theta[i] = theta + print( + "Number of particles with an angle higher than max_theta:", + len(l_theta[l_theta > max_theta]), + ) + if len(l_theta[l_theta > max_theta]) == 0: + return True + else: + return False + + +if __name__ == "__main__": + bias = True + paths = utility.get_default_test_paths( + __file__, "test084_last_vertex_splitting", output_folder="test084" + ) + + # create the simulation + sim = gate.Simulation() + + # main options + ui = sim.user_info + ui.g4_verbose = False + ui.visu = False + ui.visu_type = "vrml" + ui.check_volumes_overlap = False + ui.number_of_threads = 1 + ui.random_seed = 123456789 + + # units + m = gate.g4_units.m + km = gate.g4_units.km + mm = gate.g4_units.mm + um = gate.g4_units.um + 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 + deg = gate.g4_units.deg + gcm3 = gate.g4_units.g / gate.g4_units.cm3 + + # adapt world size + world = sim.world + world.size = [0.25 * m, 0.25 * m, 0.25 * m] + world.material = "G4_Galactic" + + ####### GEOMETRY TO IRRADIATE ############# + sim.volume_manager.material_database.add_material_weights( + "Tungsten", + ["W"], + [1], + 19.3 * gcm3, + ) + + W_tubs = sim.add_volume("Tubs", "W_box") + W_tubs.material = "Tungsten" + W_tubs.mother = world.name + + W_tubs.rmin = 0 + W_tubs.rmax = 0.4 * cm + W_tubs.dz = 0.05 * m + W_tubs.color = [0.8, 0.2, 0.1, 1] + angle_x = 45 + angle_y = 70 + angle_z = 80 + + rotation = Rotation.from_euler( + "xyz", [angle_y, angle_y, angle_z], degrees=True + ).as_matrix() + W_tubs.rotation = rotation + + if bias: + ###### Last vertex Splitting ACTOR ######### + nb_split = 10 + vertex_splitting_actor = sim.add_actor( + "LastVertexInteractionSplittingActor", "vertexSplittingW" + ) + vertex_splitting_actor.attached_to = W_tubs.name + vertex_splitting_actor.splitting_factor = nb_split + vertex_splitting_actor.angular_kill = True + vertex_splitting_actor.vector_director = [0, 0, -1] + vertex_splitting_actor.max_theta = 15 * deg + vertex_splitting_actor.batch_size = 10 + + plan = sim.add_volume("Box", "plan_phsp") + plan.material = "G4_Galactic" + plan.size = [5 * cm, 5 * cm, 1 * nm] + plan.translation = [0, 0, -1 * cm] + + if bias: + vector_director = np.array(vertex_splitting_actor.vector_director) + + ####### gamma source ########### + source = sim.add_source("GenericSource", "source1") + source.particle = "gamma" + source.n = 100000 + if bias: + source.n = source.n / nb_split + + source.position.type = "sphere" + source.position.radius = 1 * nm + source.direction.type = "momentum" + # source.direction.momentum = [0,0,-1] + source.direction.momentum = np.dot(rotation, np.array([0, 0, -1])) + source.energy.type = "mono" + source.energy.mono = 4 * MeV + + ###### LastVertexSource ############# + if bias: + source_0 = sim.add_source("LastVertexSource", "source_vertex") + source_0.n = 1 + + ####### PHASE SPACE ACTOR ############## + sim.output_dir = paths.output + phsp_actor = sim.add_actor("PhaseSpaceActor", "PhaseSpace") + phsp_actor.attached_to = plan.name + phsp_actor.attributes = [ + "EventID", + "TrackID", + "Weight", + "ParticleName", + "KineticEnergy", + "PreDirection", + "PrePosition", + "TrackCreatorProcess", + ] + if bias: + phsp_actor.output_filename = "test084_output_data_last_vertex_angular_kill.root" + + s = sim.add_actor("SimulationStatisticsActor", "Stats") + s.track_types_flag = True + sim.physics_manager.physics_list_name = "G4EmStandardPhysics_option3" + s = f"/process/em/UseGeneralProcess false" + sim.g4_commands_before_init.append(s) + + sim.physics_manager.global_production_cuts.gamma = 1 * mm + sim.physics_manager.global_production_cuts.electron = 1000 * km + sim.physics_manager.global_production_cuts.positron = 1000 * km + + output = sim.run() + print(s) + + f_data = uproot.open( + paths.output / "test084_output_data_last_vertex_angular_kill.root" + ) + arr_data = f_data["PhaseSpace"].arrays() + is_ok = validation_test(arr_data, vector_director, vertex_splitting_actor.max_theta) + utility.test_ok(is_ok) diff --git a/opengate/tests/src/test084_last_vertex_splittting.py b/opengate/tests/src/test084_last_vertex_splittting.py new file mode 100755 index 000000000..390c5d1ae --- /dev/null +++ b/opengate/tests/src/test084_last_vertex_splittting.py @@ -0,0 +1,222 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import opengate as gate +import uproot +import numpy as np +from scipy.spatial.transform import Rotation +from opengate.tests import utility + + +def validation_test(arr_ref, arr_data, nb_split): + arr_ref = arr_ref[arr_ref["ParticleName"] == "gamma"] + arr_data = arr_data[arr_data["ParticleName"] == "gamma"] + + weight_data = np.round(np.mean(arr_data["Weight"]), 4) + bool_weight = False + print( + "Weight mean is equal to", weight_data, "and need to be equal to", 1 / nb_split + ) + if weight_data == 1 / nb_split: + bool_weight = True + + bool_events = False + sigma = np.sqrt((len(arr_data["KineticEnergy"]) / nb_split)) * nb_split + nb_events_ref = len(arr_ref["KineticEnergy"]) + nb_events_data = len(arr_data["KineticEnergy"]) + print("Reference counts number:", nb_events_ref) + print("Biased counts number:", nb_events_data) + if nb_events_data - 4 * sigma <= nb_events_ref <= nb_events_data + 4 * sigma: + bool_events = True + + keys = [ + "KineticEnergy", + "PreDirection_X", + "PreDirection_Y", + "PreDirection_Z", + "PrePosition_X", + "PrePosition_Y", + ] + + bool_distrib = True + for key in keys: + ref = arr_ref[key] + data = arr_data[key] + + mean_ref = np.mean(ref) + mean_data = np.mean(data) + + std_dev_ref = np.std(ref, ddof=1) + std_dev_data = np.std(data, ddof=1) + + std_err_ref = std_dev_ref / np.sqrt(len(ref)) + std_err_data = std_dev_data / (nb_split * np.sqrt(len(data) / nb_split)) + + print( + key, + "mean ref value:", + np.round(mean_ref, 3), + "+-", + np.round(std_err_ref, 3), + ) + print( + key, + "mean data value:", + np.round(mean_data, 3), + "+-", + np.round(std_err_data, 3), + ) + + if ( + mean_data - 4 * np.sqrt(std_err_data**2 + std_err_ref**2) > mean_ref + ) or (mean_data + 4 * np.sqrt(std_err_data**2 + std_err_ref**2) < mean_ref): + bool_distrib = False + + if bool_distrib and bool_events and bool_weight: + return True + else: + return False + + +if __name__ == "__main__": + for i in range(2): + if i == 0: + bias = False + else: + bias = True + paths = utility.get_default_test_paths( + __file__, "test084_last_vertex_splitting", output_folder="test084" + ) + + # create the simulation + sim = gate.Simulation() + + # main options + ui = sim.user_info + ui.g4_verbose = False + ui.visu = False + ui.visu_type = "vrml" + ui.check_volumes_overlap = False + ui.number_of_threads = 1 + ui.random_seed = 123456789 + + # units + m = gate.g4_units.m + km = gate.g4_units.km + mm = gate.g4_units.mm + um = gate.g4_units.um + 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 + deg = gate.g4_units.deg + gcm3 = gate.g4_units.g / gate.g4_units.cm3 + + # adapt world size + world = sim.world + world.size = [0.25 * m, 0.25 * m, 0.25 * m] + world.material = "G4_Galactic" + + ####### GEOMETRY TO IRRADIATE ############# + sim.volume_manager.material_database.add_material_weights( + "Tungsten", + ["W"], + [1], + 19.3 * gcm3, + ) + + W_tubs = sim.add_volume("Tubs", "W_box") + W_tubs.material = "Tungsten" + W_tubs.mother = world.name + + W_tubs.rmin = 0 + W_tubs.rmax = 0.4 * cm + W_tubs.dz = 0.05 * m + W_tubs.color = [0.8, 0.2, 0.1, 1] + angle_x = 45 + angle_y = 70 + angle_z = 80 + + rotation = Rotation.from_euler( + "xyz", [angle_y, angle_y, angle_z], degrees=True + ).as_matrix() + W_tubs.rotation = rotation + + if bias: + ###### Last vertex Splitting ACTOR ######### + nb_split = 10 + vertex_splitting_actor = sim.add_actor( + "LastVertexInteractionSplittingActor", "vertexSplittingW" + ) + vertex_splitting_actor.attached_to = W_tubs.name + vertex_splitting_actor.splitting_factor = nb_split + vertex_splitting_actor.angular_kill = True + vertex_splitting_actor.vector_director = [0, 0, -1] + vertex_splitting_actor.max_theta = 90 * deg + vertex_splitting_actor.batch_size = 10 + + plan = sim.add_volume("Box", "plan_phsp") + plan.material = "G4_Galactic" + plan.size = [5 * cm, 5 * cm, 1 * nm] + plan.translation = [0, 0, -1 * cm] + + ####### gamma source ########### + source = sim.add_source("GenericSource", "source1") + source.particle = "gamma" + source.n = 100000 + if bias: + source.n = source.n / nb_split + + source.position.type = "sphere" + source.position.radius = 1 * nm + source.direction.type = "momentum" + # source.direction.momentum = [0,0,-1] + source.direction.momentum = np.dot(rotation, np.array([0, 0, -1])) + source.energy.type = "mono" + source.energy.mono = 4 * MeV + # + ###### LastVertexSource ############# + if bias: + source_0 = sim.add_source("LastVertexSource", "source_vertex") + source_0.n = 1 + + ####### PHASE SPACE ACTOR ############## + sim.output_dir = paths.output + phsp_actor = sim.add_actor("PhaseSpaceActor", "PhaseSpace") + phsp_actor.attached_to = plan.name + phsp_actor.attributes = [ + "EventID", + "TrackID", + "Weight", + "ParticleName", + "KineticEnergy", + "PreDirection", + "PrePosition", + "TrackCreatorProcess", + ] + if bias: + phsp_actor.output_filename = "test084_output_data_last_vertex_biased.root" + else: + phsp_actor.output_filename = "test084_output_data_last_vertex_ref.root" + + s = sim.add_actor("SimulationStatisticsActor", "Stats") + s.track_types_flag = True + sim.physics_manager.physics_list_name = "G4EmStandardPhysics_option3" + s = f"/process/em/UseGeneralProcess false" + sim.g4_commands_before_init.append(s) + + sim.physics_manager.global_production_cuts.gamma = 1 * mm + sim.physics_manager.global_production_cuts.electron = 1000 * km + sim.physics_manager.global_production_cuts.positron = 1000 * km + + output = sim.run(True) + print(s) + + f_data = uproot.open(paths.output / "test084_output_data_last_vertex_biased.root") + f_ref_data = uproot.open(paths.output / "test084_output_data_last_vertex_ref.root") + arr_data = f_data["PhaseSpace"].arrays() + arr_ref_data = f_ref_data["PhaseSpace"].arrays() + # # + is_ok = validation_test(arr_ref_data, arr_data, nb_split) + utility.test_ok(is_ok) diff --git a/opengate/tests/src/test084_last_vertex_splittting_angular_kill.py b/opengate/tests/src/test084_last_vertex_splittting_angular_kill.py new file mode 100755 index 000000000..d1fa56cf2 --- /dev/null +++ b/opengate/tests/src/test084_last_vertex_splittting_angular_kill.py @@ -0,0 +1,171 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import opengate as gate +import uproot +import numpy as np +from scipy.spatial.transform import Rotation +from opengate.tests import utility + + +def validation_test(arr_data, vector_director, max_theta): + arr_data = arr_data[arr_data["ParticleName"] == "gamma"] + qt_mvt_data = arr_data[["PreDirection_X", "PreDirection_Y", "PreDirection_Z"]] + mom_data = np.zeros((len(qt_mvt_data["PreDirection_X"]), 3)) + mom_data[:, 0] += np.array(qt_mvt_data["PreDirection_X"]) + mom_data[:, 1] += np.array(qt_mvt_data["PreDirection_Y"]) + mom_data[:, 2] += np.array(qt_mvt_data["PreDirection_Z"]) + + l_theta = np.zeros(len(mom_data[:, 0])) + + for i in range(len(l_theta)): + theta = np.arccos(mom_data[i].dot(vector_director)) + l_theta[i] = theta + print( + "Number of particles with an angle higher than max_theta:", + len(l_theta[l_theta > max_theta]), + ) + if len(l_theta[l_theta > max_theta]) == 0: + return True + else: + return False + + +if __name__ == "__main__": + bias = True + paths = utility.get_default_test_paths( + __file__, "test084_last_vertex_splitting", output_folder="test084" + ) + + # create the simulation + sim = gate.Simulation() + + # main options + ui = sim.user_info + ui.g4_verbose = False + ui.visu = False + ui.visu_type = "vrml" + ui.check_volumes_overlap = False + ui.number_of_threads = 1 + ui.random_seed = 123456789 + + # units + m = gate.g4_units.m + km = gate.g4_units.km + mm = gate.g4_units.mm + um = gate.g4_units.um + 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 + deg = gate.g4_units.deg + gcm3 = gate.g4_units.g / gate.g4_units.cm3 + + # adapt world size + world = sim.world + world.size = [0.25 * m, 0.25 * m, 0.25 * m] + world.material = "G4_Galactic" + + ####### GEOMETRY TO IRRADIATE ############# + sim.volume_manager.material_database.add_material_weights( + "Tungsten", + ["W"], + [1], + 19.3 * gcm3, + ) + + W_tubs = sim.add_volume("Tubs", "W_box") + W_tubs.material = "Tungsten" + W_tubs.mother = world.name + + W_tubs.rmin = 0 + W_tubs.rmax = 0.4 * cm + W_tubs.dz = 0.05 * m + W_tubs.color = [0.8, 0.2, 0.1, 1] + angle_x = 45 + angle_y = 70 + angle_z = 80 + + rotation = Rotation.from_euler( + "xyz", [angle_y, angle_y, angle_z], degrees=True + ).as_matrix() + W_tubs.rotation = rotation + + if bias: + ###### Last vertex Splitting ACTOR ######### + nb_split = 10 + vertex_splitting_actor = sim.add_actor( + "LastVertexInteractionSplittingActor", "vertexSplittingW" + ) + vertex_splitting_actor.attached_to = W_tubs.name + vertex_splitting_actor.splitting_factor = nb_split + vertex_splitting_actor.angular_kill = True + vertex_splitting_actor.vector_director = [0, 0, -1] + vertex_splitting_actor.max_theta = 15 * deg + vertex_splitting_actor.batch_size = 10 + + plan = sim.add_volume("Box", "plan_phsp") + plan.material = "G4_Galactic" + plan.size = [5 * cm, 5 * cm, 1 * nm] + plan.translation = [0, 0, -1 * cm] + + if bias: + vector_director = np.array(vertex_splitting_actor.vector_director) + + ####### gamma source ########### + source = sim.add_source("GenericSource", "source1") + source.particle = "gamma" + source.n = 100000 + if bias: + source.n = source.n / nb_split + + source.position.type = "sphere" + source.position.radius = 1 * nm + source.direction.type = "momentum" + # source.direction.momentum = [0,0,-1] + source.direction.momentum = np.dot(rotation, np.array([0, 0, -1])) + source.energy.type = "mono" + source.energy.mono = 4 * MeV + + ###### LastVertexSource ############# + if bias: + source_0 = sim.add_source("LastVertexSource", "source_vertex") + source_0.n = 1 + + ####### PHASE SPACE ACTOR ############## + sim.output_dir = paths.output + phsp_actor = sim.add_actor("PhaseSpaceActor", "PhaseSpace") + phsp_actor.attached_to = plan.name + phsp_actor.attributes = [ + "EventID", + "TrackID", + "Weight", + "ParticleName", + "KineticEnergy", + "PreDirection", + "PrePosition", + "TrackCreatorProcess", + ] + if bias: + phsp_actor.output_filename = "test084_output_data_last_vertex_angular_kill.root" + + s = sim.add_actor("SimulationStatisticsActor", "Stats") + s.track_types_flag = True + sim.physics_manager.physics_list_name = "G4EmStandardPhysics_option3" + s = f"/process/em/UseGeneralProcess false" + sim.g4_commands_before_init.append(s) + + sim.physics_manager.global_production_cuts.gamma = 1 * mm + sim.physics_manager.global_production_cuts.electron = 1000 * km + sim.physics_manager.global_production_cuts.positron = 1000 * km + + output = sim.run() + print(s) + + f_data = uproot.open( + paths.output / "test084_output_data_last_vertex_angular_kill.root" + ) + arr_data = f_data["PhaseSpace"].arrays() + is_ok = validation_test(arr_data, vector_director, vertex_splitting_actor.max_theta) + utility.test_ok(is_ok)