Allow radioactive isotope as source and arbitrary source energy

JuliaPhysics · Jul 30, 2024 · f0a19c7 · f0a19c7
1 parent 7af1497
commit f0a19c7
Show file tree

Hide file tree

Showing 3 changed files with 54 additions and 44 deletions.
diff --git a/ext/Geant4/g4jl_application.jl b/ext/Geant4/g4jl_application.jl
@@ -66,7 +66,7 @@ end
     direction::G4ThreeVector = G4ThreeVector(0, 0, 0)
 end
 
-function SSDGenerator(source::SolidStateDetectors.SSDSource;  kwargs...)
+function SSDGenerator(source::SolidStateDetectors.MonoenergeticSource;  kwargs...)
 
     @assert source.direction isa SolidStateDetectors.CartesianVector || source.direction == :isotropic
 
@@ -75,16 +75,16 @@ function SSDGenerator(source::SolidStateDetectors.SSDSource;  kwargs...)
         gun = data.gun = move!(G4ParticleGun())
         particle = data.particle = FindParticle(source.particle_type)
         data.position = G4ThreeVector(
-		source.position.x * Geant4.SystemOfUnits.meter,
-		source.position.y * Geant4.SystemOfUnits.meter,
-		source.position.z * Geant4.SystemOfUnits.meter
-	)
+            source.position.x * Geant4.SystemOfUnits.meter,
+            source.position.y * Geant4.SystemOfUnits.meter,
+            source.position.z * Geant4.SystemOfUnits.meter
+        )
         SetParticlePosition(gun, data.position)
         if source.direction isa SolidStateDetectors.CartesianVector
           data.direction = G4ThreeVector(source.direction.x, source.direction.y, source.direction.z)
         end
         SetParticleMomentumDirection(gun, data.direction)
-        SetParticleEnergy(gun, 2.615GeV)
+        SetParticleEnergy(gun, ustrip(u"MeV", source.energy))
         SetParticleDefinition(gun, particle)
     end
 
@@ -102,43 +102,41 @@ function SSDGenerator(source::SolidStateDetectors.SSDSource;  kwargs...)
     G4JLPrimaryGenerator("SSDGenerator", data; init_method=_init, generate_method=_gen)
 end
 
-#=
-@with_kw mutable struct GeneratorC3aData <: G4JLGeneratorData
-    gun::Union{Nothing, CxxPtr{G4ParticleGun}} = nothing
-    ion::Union{Nothing, CxxPtr{G4ParticleDefinition}} = nothing
-    Z::Int64 = 82
-    A::Int64 = 212
-    ionCharge::Float64 = 0eplus
-    excitEnergy::Float64 = 0keV
-    position::G4ThreeVector = G4ThreeVector(5cm, 0cm, 5cm)
-    direction::G4ThreeVector = G4ThreeVector(-1, 0, 0)
-end
-function GeneratorC3a(;kwargs...)
-    data = GeneratorC3aData(;kwargs...)
-    function _init(data::GeneratorC3aData, ::Any)
+function SSDGenerator(source::SolidStateDetectors.IsotopeSource;  kwargs...)
+
+    @assert source.direction isa SolidStateDetectors.CartesianVector || source.direction == :isotropic
+
+    data = GeneratorData(;kwargs...)
+    function _init(data::GeneratorData, ::Any)
         gun = data.gun = move!(G4ParticleGun())
-        data.direction = direction = G4ThreeVector(-10,rand()*2 - 1,rand()*2 - 1)
+        data.position = G4ThreeVector(
+            source.position.x * Geant4.SystemOfUnits.meter,
+            source.position.y * Geant4.SystemOfUnits.meter,
+            source.position.z * Geant4.SystemOfUnits.meter
+        )
         SetParticlePosition(gun, data.position)
+        if source.direction isa SolidStateDetectors.CartesianVector
+            data.direction = G4ThreeVector(source.direction.x, source.direction.y, source.direction.z)
+        end
         SetParticleMomentumDirection(gun, data.direction)
-        SetParticleEnergy(gun, 0eV)
+        # SetParticleEnergy(gun, source.excitEnergy)
     end
-    function _gen(evt::G4Event, data::GeneratorC3aData)::Nothing
-        if isnothing(data.ion)  # late initialize (after physics processes)
-            ion = data.ion = GetIon(data.Z, data.A, data.excitEnergy)
-            SetParticleDefinition(data.gun, ion)
-            SetParticleCharge(data.gun, data.ionCharge)
+
+    function _gen(evt::G4Event, data::GeneratorData)::Nothing
+        if isnothing(data.particle)  # late initialize (after physics processes)
+            data.particle = GetIon(source.Z, source.A, source.excitEnergy)
+            SetParticleDefinition(data.gun, data.particle)
+            SetParticleCharge(data.gun, source.ionCharge)
         end
-        #position = data.position + G4ThreeVector((rand()-0.5)*1cm, (rand()-0.5)*1cm, (rand()-0.5)*1cm)
-        #SetParticlePosition(data.gun, position)
-        data.direction = direction = G4ThreeVector(-10,rand()*2 - 1,rand()*2 - 1)
-        SetParticleMomentumDirection(data.gun, direction)
-        #data.direction = direction
-        SetParticlePosition(data.gun, data.position)
+        if source.direction == :isotropic
+            ϕ = rand()*2π 
+            θ = acos(1 - 2*rand())
+            direction = G4ThreeVector(cos(ϕ)*sin(θ), sin(ϕ)*sin(θ), cos(θ))
+            data.direction = direction
+        end
+        SetParticleMomentumDirection(data.gun, data.direction)
+        SetParticlePosition(data.gun, data.position) # needed ?
         GeneratePrimaryVertex(data.gun, CxxPtr(evt))
     end
-    G4JLPrimaryGenerator("GeneratorC3a", data; init_method=_init, generate_method=_gen)
-end
-=#
-
-
-
+    G4JLPrimaryGenerator("SSDGenerator", data; init_method=_init, generate_method=_gen)
+end
diff --git a/ext/SolidStateDetectorsGeant4Ext.jl b/ext/SolidStateDetectorsGeant4Ext.jl
@@ -109,7 +109,7 @@ end
 
 function Geant4.G4JLApplication(
     sim::Union{SolidStateDetectors.Simulation, String},
-    source::SSDSource;
+    source::SolidStateDetectors.AbstractParticleSource;
     physics_type = SSDPhysics,
     endeventaction_method = endeventaction,
     verbose = true,

diff --git a/src/ChargeCloudModels/ParticleTypes.jl b/src/ChargeCloudModels/ParticleTypes.jl
@@ -18,10 +18,22 @@ radius_guess(charge::T, ::Type{Beta}) where {T} = T(0.0005)
 radius_guess(charge::T, ::Type{Gamma}) where {T} = T(0.0005) 
 
 
-struct SSDSource
+abstract type AbstractParticleSource end
+
+struct MonoenergeticSource{T} <: AbstractParticleSource
   particle_type::String
-  position::SolidStateDetectors.CartesianPoint
-  direction::Union{Symbol, SolidStateDetectors.CartesianVector}
+  energy::RealQuantity{T}
+  position::CartesianPoint
+  direction::Union{Symbol, CartesianVector}
+end
+
+struct IsotopeSource{T} <: AbstractParticleSource
+  Z::Integer
+  A::Integer
+  ionCharge::T
+  excitEnergy::T
+  position::CartesianPoint
+  direction::Union{Symbol, CartesianVector}
 end
 
-export SSDSource
+export AbstractParticleSource, MonoenergeticSource, IsotopeSource