Defining an extended non-point-like source for event generation

[rajxabc]

Hi Felix ,
could you provide some example of "virtual_drift_volumes" usage.
I tried like this ,
virtual_drift_volume: - name: srcVOL model: dead geometry: box: hX: 250 hY: 250 hZ: 0.001 origin: x: 0 y: 0 z: 530.0005

in sim it shows "missing,missing"

what i am trying is : to define virtual thin volume and use it as initial 2D source ..something like
pt (rand(x), rand(y), z)
if pt in sim.detector.virtual_drift_volumes)
source will be set to pt .
is this possible ?
Thanks

[fhagemann]

Two things:

1. In the config file, you would need to use virtual_drift_volumes with an s at the end, see here, so that is probably why it is listed as "missing", because it's not actually read in with the wrong key word.
2. Virtual drift volumes should be used to modify the drift, not to define sources. Here, it would probably make most sense to define a new source type that also draws the initial position of the source to be random within a volume.

[fhagemann]

SolidStateDetectors.jl/ext/Geant4/g4jl_application.jl

Lines 62 to 100 in e4daebe

	function SSDGenerator(source::SolidStateDetectors.MonoenergeticSource; kwargs...)
	iszero(source.direction) && throw(ArgumentError("The direction of the source cannot be zero."))
	data = GeneratorData(;kwargs...)
	function _init(data::GeneratorData, ::Any)
	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
	)
	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)
	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
	G4JLPrimaryGenerator("SSDGenerator", data; init_method=_init, generate_method=_gen)
	end

This is how it is done for the MonoenergeticSource.
You would then, instead of definiting data.position once at the beginning of the application in _init, have to define some algorithm to draw a random number in _gen, e.g. something like (if you were to use the MonoenergeticSource as input.

function SSDGenerator(source::SolidStateDetectors.MonoenergeticSource;  kwargs...)

    iszero(source.direction) && throw(ArgumentError("The direction of the source cannot be zero."))

    data = GeneratorData(;kwargs...)

    function _init(data::GeneratorData, ::Any)
        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
        )
        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)
    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
        
        # here comes the new code to generate a random position 
        # (pseudo-code, please adjust to your needs)
        position = G4ThreeVector( 
             rand() * Geant4.SystemOfUnits.meter, 
             rand() * Geant4.SystemOfUnits.meter, 
             rand() * Geant4.SystemOfUnits.meter 
         ) 
         SetParticlePosition(data.gun, position) 
         # new code end

        GeneratePrimaryVertex(data.gun, CxxPtr(evt))
    end

    G4JLPrimaryGenerator("SSDGenerator", data; init_method=_init, generate_method=_gen)
end

[fhagemann]

For a source that looks like a Box, you can try this:

using SolidStateDetectors, Geant4, Unitful
const G4Ext = Base.get_extension(SolidStateDetectors, :SolidStateDetectorsGeant4Ext)

using SolidStateDetectors: AbstractParticleSource, RealQuantity
using SolidStateDetectors.ConstructiveSolidGeometry: Box, _transform_into_global_coordinate_system

struct MonoenergeticBoxSource <: AbstractParticleSource
    particle_type::String
    energy::RealQuantity
    position::Box
end

function G4Ext.SSDGenerator(source::MonoenergeticBoxSource;  kwargs...)

    data = G4Ext.GeneratorData(; kwargs...)

    function _init(data::G4Ext.GeneratorData, ::Any)
        gun = data.gun = move!(G4ParticleGun())
        particle = data.particle = FindParticle(source.particle_type)
        data.position = G4ThreeVector(source.position.origin .* Geant4.SystemOfUnits.meter...)
        SetParticlePosition(gun, data.position)
        data.direction = G4ThreeVector(1,0,0) # arbitrary initial direction
        SetParticleMomentumDirection(gun, data.direction)
        SetParticleEnergy(gun, ustrip(u"MeV", source.energy))
        SetParticleDefinition(gun, particle)
    end

    function _gen(evt::G4Event, data::G4Ext.GeneratorData, box::Box{T} = source.position)::Nothing where {T}
        
        # generate random emission direction
        ϕ = rand()*2π
        θ = acos(1 - 2*rand())
        direction = G4ThreeVector(sin(θ)*cos(ϕ), sin(θ)*sin(ϕ), cos(θ)) 
        SetParticleMomentumDirection(data.gun, direction)

        # generate random point in the Box
        pos =  _transform_into_global_coordinate_system(
            CartesianPoint{T}(
                box.hX * (2 * rand() - 1),
                box.hY * (2 * rand() - 1),
                box.hZ * (2 * rand() - 1)
            ), box)
        position = G4ThreeVector(Float64.(pos) .* Geant4.SystemOfUnits.meter...) 
        SetParticlePosition(data.gun, position)
        
        GeneratePrimaryVertex(data.gun, CxxPtr(evt))
    end

    G4JLPrimaryGenerator("SSDGenerator", data; init_method=_init, generate_method=_gen)
end

sim = Simulation(SSD_examples[:InvertedCoax])
source = MonoenergeticBoxSource("gamma", 2.615u"MeV", Box(hX = 0.04, hY = 0.04, hZ = 1e-6, origin = CartesianPoint(0,0,0.1)))

app = G4JLApplication(sim, source)

run_geant4_simulation(app, 10_000)

Maybe add a passive (!! not a virtual_drift_volume) to your simulation, so that you can also capture effects like self-absorption of the source material, or X-ray activation, if that might be expected.

Note that you might still need to add a plot recipe to plot the source in a 3D plot.

[fhagemann]

An example plot recipe would look like this:

using RecipesBase
@recipe function f(m::MonoenergeticBoxSource)
    @series begin
        seriescolor := :gray
        markerstrokewidth --> 1
        linecolor := :black
        label --> string(typeof(m).name.name)
        m.position
    end
end

using Plots
plot(sim.detector)
plot!(source)

[rajxabc]

Thanks Felix ,
i will modify it for hollow rectangle and hollow disc source
thanks again

[fhagemann]

With #492, the custom plot recipe is also no longer needed.

[rajxabc]

One related question , if i want to check
"CatesianPoint in Box" or "CatesianPoint in CSGdifference"
How to do it ?
is there some function to convert them to geometry ?

[fhagemann]

Maybe also have a look at #526, where instead of a point-like source, you might be able to implement another position distribution
(see page 29 in the Geant4 Book for Application Developers): Source position and structure.pdf

SolidStateDetectors.jl/ext/Geant4/g4jl_application.jl

Line 84 in a28411a

SetPosDisType(posdist, "Point")

[rajxabc]

ok thanks Felix , for a time being custom box method was working fine for me.