Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Support tori with partial θ ranges in Geant4 extension #459

Merged
merged 2 commits into from
Feb 13, 2025
Merged

Support tori with partial θ ranges in Geant4 extension #459

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
e0e2451
Bug fix in `parse_origin` and `parse_rotation` in Geant4 extension
fhagemann Feb 12, 2025
69104fe
Support tori with partial θ ranges in Geant4 extension
fhagemann Feb 12, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
78 changes: 62 additions & 16 deletions ext/Geant4/io_gdml.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,13 +14,13 @@ using LightXML
@inline parse_origin(e::AbstractVolumePrimitive) = origin(e)

# Returns position vector for leftmost Primitive in geometry tree
@inline parse_origin(e::AbstractConstructiveGeometry) = parse_origin(e.a)
@inline parse_origin(e::AbstractConstructiveGeometry) = parse_origin(has_volume(e.a) ? e.a : e.b)

# Returns rotation matrix for Primitive relative to standard basis
@inline parse_rotation(e::AbstractVolumePrimitive) = rotation(e)

# Returns rotation matrix for leftmost Primitive in geometry tree
@inline parse_rotation(e::AbstractConstructiveGeometry) = parse_rotation(e.a)
@inline parse_rotation(e::AbstractConstructiveGeometry) = parse_rotation(has_volume(e.a) ? e.a : e.b)


# Add <position> to <define> section, referenced in the geometry definition (in <solids>) via the name
Expand Down Expand Up @@ -322,21 +322,26 @@ function parse_geometry(e::Ellipsoid{T,<:Any, TR, Nothing, Nothing}, x_solids::X
nothing
end

# @inline parse_θ(::Type{T}, θ::Nothing) where {T} = (T(0), T(360), "deg")
# @inline parse_θ(::Type{T}, θ::Tuple{T, T}) where {T} = throw(AssertionError("Torus: θ values are currently not supported"))
function parse_geometry(e::Torus, x_solids::XMLElement, x_define::XMLElement, id::Integer, pf::AbstractString, v::Bool)::Nothing
if has_volume(e, v)
function parse_geometry(e::Torus{T,<:Any,<:Any,Nothing,Nothing}, x_solids::XMLElement, x_define::XMLElement, id::Integer, pf::AbstractString, v::Bool)::Nothing where {T}

# if e.r_torus == 0, the torus is basically a sphere => convert it to a sphere
if iszero(e.r_torus)
parse_geometry(
if isa(e.r_tube, Tuple)
CSGDifference(
Ellipsoid{T}(r = e.r_tube[2], origin = e.origin, rotation = e.rotation),
Ellipsoid{T}(r = e.r_tube[1], origin = e.origin, rotation = e.rotation)
)
else
Ellipsoid{T}(r = e.r_tube, origin = e.origin, rotation = e.rotation)
end,
x_solids, x_define, id, pf, v
)

elseif has_volume(e, v)
y = new_child(x_solids, "torus")
# phi, aunit = parse_φ(T, e.φ)
# theta1, theta2, aunit = parse_θ(T, e.θ)

rmin, rmax = if isa(e.r_tube, Tuple)
e.r_tube[1], e.r_tube[2]
else
0, e.r_tube
end

# rmin = 0
rmin, rmax = isa(e.r_tube, Tuple) ? e.r_tube : (0, e.r_tube)

set_attributes(y, OrderedDict(
"name" => pf * string(id),
Expand All @@ -348,7 +353,48 @@ function parse_geometry(e::Torus, x_solids::XMLElement, x_define::XMLElement, id
"lunit" => SolidStateDetectors.internal_length_unit,
"aunit" => "deg"
))

end
end

function _rectangle_point(x::T)::Tuple{T,T} where {T <: Real}
return (
clamp(tand(abs(mod(x-90,360)-180)-90),-1,1),
clamp(tand(abs(mod( x,360)-180)-90),-1,1)
)
end
Comment on lines +359 to +364
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Some back info: For tori with partial θ ranges, I intersect a full torus with a PolyCone.
The function _rectangle_point returns a point on a rectangle at a given θ (in degrees):

using SolidStateDetectors
using Geant4
using Plots

const G4 = Base.get_extension(SolidStateDetectors, :SolidStateDetectorsGeant4Ext)

base_plot = plot([1,1,-1,-1,1],[1,-1,-1,1,1], color = :black, ratio = 1, size = (400,400), label = "", xlims = (-2,2), ylims = (-2,2))
@gif for x in 0:5.:360-5
    p = G4._rectangle_point(x)
    plot(deepcopy(base_plot))
    scatter!([p[1]], [p[2]], label = "", color = :black)
end

tmp


function parse_geometry(e::Torus{T,<:Any,<:Any,Nothing,Tuple{T,T}}, x_solids::XMLElement, x_define::XMLElement, id::Integer, pf::AbstractString, v::Bool)::Nothing where {T}
if has_volume(e, v)

theta1::T, theta2::T = rad2deg.(e.θ)
rmax = isa(e.r_tube, Tuple) ? e.r_tube[2] : e.r_tube

parse_geometry(
CSGIntersection(
# create full-θ Torus
Torus(r_torus = e.r_torus, r_tube = e.r_tube, φ = e.φ, θ = nothing, origin = e.origin, rotation = e.rotation),

# create Polycone to intersect with
begin
points = Tuple{T,T}[(zero(T), zero(T))]
push!(points, _rectangle_point(theta2))
tmp::T = theta2 - mod(theta2 - 45, 90)
while theta1 < tmp && tmp > 0
push!(points, _rectangle_point(tmp))
tmp -= 90
end
push!(points, _rectangle_point(theta1))
push!(points, (zero(T), zero(T)))
Comment on lines +379 to +387
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I can use this to create a shape that spans over a rectangle with given theta1 and theta2:

T = Float64
# generate random values for theta1 and theta2
theta2 = 360 * rand(T)
theta1 = theta2 * rand(T)

points = Tuple{T,T}[(zero(T), zero(T))]
push!(points, G4._rectangle_point(theta2))
tmp::T = theta2 - mod(theta2 - 45, 90)
while theta1 < tmp && tmp > 0
    push!(points, G4._rectangle_point(tmp))
    tmp -= 90
end
push!(points, G4._rectangle_point(theta1))
push!(points, (zero(T), zero(T)))

plot(deepcopy(base_plot), legend = :topleft)
plot!(getindex.(points,1), getindex.(points,2), color = :red, lw = 2, label = "θ1 = $(theta1),\nθ2 = $(theta2)")

image

Polycone(
r = getindex.(points,1) .* rmax .+ e.r_torus,
z = getindex.(points,2) .* rmax,
origin = e.origin,
rotation = e.rotation
)
Comment on lines +388 to +393
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

... to then intersect the full-θ Torus with the newly defined Polycone

end
),
x_solids, x_define, id, pf, v
)
end
end

Expand Down
Loading