Merge pull request JuliaPhysics#526 from JuliaPhysics/IsotopeSource

Fix bug in source emission in `Geant4` extension

Author: fhagemann

ext/Geant4/g4jl_application.jl

@@ -13,17 +13,17 @@ struct SSDPhysics <: G4VUserPhysicsList
 end
 
 struct SDData{T} <: G4JLSDData
-  detectorHits::DetectorHits
-  SDData{T}() where {T} = new(DetectorHits((evtno = Int32[], detno = Int32[], thit = typeof(zero(T)*u"s")[], edep = typeof(zero(T)*u"keV")[], pos = SVector{3,typeof(zero(T)*u"mm")}[])))
+    detectorHits::DetectorHits
+    SDData{T}() where {T} = new(DetectorHits((evtno = Int32[], detno = Int32[], thit = typeof(zero(T)*u"s")[], edep = typeof(zero(T)*u"keV")[], pos = SVector{3,typeof(zero(T)*u"mm")}[])))
 end
 
 function _initialize(::G4HCofThisEvent, data::SDData)::Nothing
-  empty!(data.detectorHits)
-  return
+    empty!(data.detectorHits)
+    return
 end
 
 function _endOfEvent(::G4HCofThisEvent, data::SDData)::Nothing
-  return
+    return
 end
 
 function _processHits(step::G4Step, ::G4TouchableHistory, data::SDData{T})::Bool where {T}
@@ -48,12 +48,19 @@ function _processHits(step::G4Step, ::G4TouchableHistory, data::SDData{T})::Bool
 end
 
 function endeventaction(evt::G4Event, app::G4JLApplication)
+    # direction = evt |> (evt -> GetPrimaryVertex(evt, 0)) |> (vertex -> GetPrimary(vertex, 0)) |> GetMomentumDirection
+    # pos = evt |> (evt -> GetPrimaryVertex(evt, 0)) |> GetPosition
+    # d = [x(direction), y(direction), z(direction)]
+    # @info [x(pos), y(pos), z(pos)], normalize(d)
     return
 end
 
 
+# We are using G4SingleParticleSource instead of G4ParticleGun for uniformity 
+# across MonoenergeticSource and IsotopeSource, see discussion
+# https://github.com/JuliaPhysics/SolidStateDetectors.jl/pull/526
 @with_kw mutable struct GeneratorData <: G4JLGeneratorData
-    gun::Union{Nothing, CxxPtr{G4ParticleGun}, CxxPtr{G4SingleParticleSource}} = nothing
+    gun::Union{Nothing, CxxPtr{G4SingleParticleSource}} = nothing
     particle::Union{Nothing, CxxPtr{G4ParticleDefinition}} = nothing
     position::G4ThreeVector = G4ThreeVector(0, 0, 0)
     direction::G4ThreeVector = G4ThreeVector(0, 0, 0)
@@ -66,31 +73,39 @@ function SSDGenerator(source::SolidStateDetectors.MonoenergeticSource;  kwargs..
     data = GeneratorData(;kwargs...)
 
     function _init(data::GeneratorData, ::Any)
-        gun = data.gun = move!(G4ParticleGun())
+        data.gun = move!(G4SingleParticleSource())
         particle = data.particle = FindParticle(source.particle_type)
+        SetParticleDefinition(data.gun, particle)
+
+        # Place the SingleParticleSource as point source
         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)
-        data.direction = G4ThreeVector(source.direction.x, source.direction.y, source.direction.z)
-        SetParticleMomentumDirection(gun, data.direction)
-        SetParticleEnergy(gun, ustrip(u"MeV", source.energy))
-        SetParticleDefinition(gun, particle)
+        posdist = GetPosDist(data.gun)
+        SetPosDisType(posdist, "Point")
+        SetCentreCoords(posdist, data.position)
+
+        # Define the source emission
+        d::CartesianVector = normalize(source.direction)
+        a::CartesianVector = normalize(d × (abs(d.x) == 1 ? CartesianVector(0,1,0) : CartesianVector(1,0,0)))
+        b::CartesianVector = normalize(a × d)
+        data.direction = G4ThreeVector(d.x, d.y, d.z)
+
+        angdist = GetAngDist(data.gun)
+        SetAngDistType(angdist, "iso")
+        DefineAngRefAxes(angdist, "angref1", G4ThreeVector(a...))
+        DefineAngRefAxes(angdist, "angref2", G4ThreeVector(b...))
+        SetMinTheta(angdist, 0.0)
+        SetMaxTheta(angdist, ustrip(NoUnits, source.opening_angle))  # convert to radians
+        # SetParticleMomentumDirection(angdist, data.direction)
+
+        # set energy
+        SetMonoEnergy(GetEneDist(data.gun), ustrip(u"MeV", source.energy))
     end
 
     function _gen(evt::G4Event, data::GeneratorData)::Nothing
-        if !iszero(source.opening_angle)
-            d::CartesianVector = normalize(source.direction)
-            a::CartesianVector = normalize(d × (abs(d.x) == 1 ? CartesianVector(0,1,0) : CartesianVector(1,0,0)))
-            b::CartesianVector = normalize(a × d)
-            ϕ = rand()*2π
-            θ = acos(1 - (1 - cos(source.opening_angle))*rand())
-            v = (cos(θ) * d + sin(θ) * (cos(ϕ) * a + sin(ϕ) * b))
-            direction = G4ThreeVector(v.x, v.y, v.z)
-            SetParticleMomentumDirection(data.gun, direction)
-        end
         GeneratePrimaryVertex(data.gun, CxxPtr(evt))
     end
 
@@ -103,16 +118,33 @@ function SSDGenerator(source::SolidStateDetectors.IsotopeSource; kwargs...)
 
     data = GeneratorData(;kwargs...)
     function _init(data::GeneratorData, ::Any)
-        gun = data.gun = move!(G4SingleParticleSource())
+        data.gun = move!(G4SingleParticleSource())
+
+        # Place the SingleParticleSource as point source
         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)
-        data.direction = G4ThreeVector(source.direction.x, source.direction.y, source.direction.z)
-        SetParticleMomentumDirection(GetAngDist(gun), data.direction)
-        SetMonoEnergy(gun |> GetEneDist, 0.0 * Geant4.SystemOfUnits.MeV)
+        posdist = GetPosDist(data.gun)
+        SetPosDisType(posdist, "Point")
+        SetCentreCoords(posdist, data.position)
+
+        # Define the source emission
+        d::CartesianVector = normalize(source.direction)
+        a::CartesianVector = normalize(d × (abs(d.x) == 1 ? CartesianVector(0,1,0) : CartesianVector(1,0,0)))
+        b::CartesianVector = normalize(a × d)
+        data.direction = G4ThreeVector(d.x, d.y, d.z)
+
+        angdist = GetAngDist(data.gun)
+        SetAngDistType(angdist, "iso")
+        DefineAngRefAxes(angdist, "angref1", G4ThreeVector(a...))
+        DefineAngRefAxes(angdist, "angref2", G4ThreeVector(b...))
+        SetMinTheta(angdist, 0.0)
+        SetMaxTheta(angdist, ustrip(NoUnits, source.opening_angle))  # convert to radians
+
+        # Set kinetic energy of the decaying isotope to zero ("decay at rest")
+        SetMonoEnergy(GetEneDist(data.gun), 0.0 * Geant4.SystemOfUnits.MeV)
     end
 
     function _gen(evt::G4Event, data::GeneratorData)::Nothing
@@ -121,16 +153,6 @@ function SSDGenerator(source::SolidStateDetectors.IsotopeSource; kwargs...)
             SetParticleDefinition(data.gun, data.particle)
             SetParticleCharge(data.gun, source.ionCharge)
         end
-        if !iszero(source.opening_angle)
-            d::CartesianVector = normalize(source.direction)
-            a::CartesianVector = normalize(d × (abs(d.x) == 1 ? CartesianVector(0,1,0) : CartesianVector(1,0,0)))
-            b::CartesianVector = normalize(a × d)
-            ϕ = rand()*2π
-            θ = acos(1 - (1 - cos(source.opening_angle))*rand())
-            v = (cos(θ) * d + sin(θ) * (cos(ϕ) * a + sin(ϕ) * b))
-            direction = G4ThreeVector(v.x, v.y, v.z)
-            SetParticleMomentumDirection(GetAngDist(data.gun), direction)
-        end
         GeneratePrimaryVertex(data.gun, CxxPtr(evt))
     end