From 82fb686dcea63d5e317dbd41627204966a95377a Mon Sep 17 00:00:00 2001 From: Felix Hagemann Date: Tue, 4 Nov 2025 14:46:12 -0800 Subject: [PATCH 1/9] Fix plotting methods for `Box` --- src/ConstructiveSolidGeometry/SurfacePrimitives/Polygon.jl | 2 +- src/ConstructiveSolidGeometry/VolumePrimitives/Box.jl | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/ConstructiveSolidGeometry/SurfacePrimitives/Polygon.jl b/src/ConstructiveSolidGeometry/SurfacePrimitives/Polygon.jl index 8fdee1446..39a34b1f5 100644 --- a/src/ConstructiveSolidGeometry/SurfacePrimitives/Polygon.jl +++ b/src/ConstructiveSolidGeometry/SurfacePrimitives/Polygon.jl @@ -35,7 +35,7 @@ function _sample_excluding_border(t::Triangle{T}, spacing::T)::Vector{CartesianP p2 = t.points[2] + (push/sφ)*(-u/nu + w/nw) if (p2-p1) ⋅ u > 0 su = norm(p2-p1)/nu - [(su*a*u + su*b*v) + p1 + [p1 + (su*a*u + su*b*v) for a in range(0, stop = 1, length = max(2, 1 + Int(ceil(su*nu*sθ/spacing)))) for b in range(0, stop = 1 - a, length = max(2, 1 + Int(ceil(su*nv*(1 - a)/spacing))))] else diff --git a/src/ConstructiveSolidGeometry/VolumePrimitives/Box.jl b/src/ConstructiveSolidGeometry/VolumePrimitives/Box.jl index 58a5b5feb..35bce8dd8 100644 --- a/src/ConstructiveSolidGeometry/VolumePrimitives/Box.jl +++ b/src/ConstructiveSolidGeometry/VolumePrimitives/Box.jl @@ -145,12 +145,12 @@ function vertices(b::Box{T}) where {T} return _transform_into_global_coordinate_system.(( CartesianPoint{T}(-b.hX, -b.hY, -b.hZ), CartesianPoint{T}(+b.hX, -b.hY, -b.hZ), - CartesianPoint{T}(-b.hX, +b.hY, -b.hZ), CartesianPoint{T}(+b.hX, +b.hY, -b.hZ), + CartesianPoint{T}(-b.hX, +b.hY, -b.hZ), CartesianPoint{T}(-b.hX, -b.hY, +b.hZ), CartesianPoint{T}(+b.hX, -b.hY, +b.hZ), - CartesianPoint{T}(-b.hX, +b.hY, +b.hZ), - CartesianPoint{T}(+b.hX, +b.hY, +b.hZ) + CartesianPoint{T}(+b.hX, +b.hY, +b.hZ), + CartesianPoint{T}(-b.hX, +b.hY, +b.hZ) ), Ref(b)) end From 3436c02a2263338ba36fa163e45f8ad333dd6fb6 Mon Sep 17 00:00:00 2001 From: Felix Hagemann Date: Tue, 4 Nov 2025 17:27:06 -0800 Subject: [PATCH 2/9] Add default `markersize` and `markerstrokewidth` for `samplesurface` plots --- src/ConstructiveSolidGeometry/plotting/CSG/CSG.jl | 9 +++++++++ .../plotting/SurfacePrimitives/SurfacePrimitives.jl | 9 ++++++++- .../plotting/VolumePrimitives/VolumePrimitives.jl | 9 +++++++++ 3 files changed, 26 insertions(+), 1 deletion(-) diff --git a/src/ConstructiveSolidGeometry/plotting/CSG/CSG.jl b/src/ConstructiveSolidGeometry/plotting/CSG/CSG.jl index 206a9b20f..039a0c271 100644 --- a/src/ConstructiveSolidGeometry/plotting/CSG/CSG.jl +++ b/src/ConstructiveSolidGeometry/plotting/CSG/CSG.jl @@ -8,6 +8,15 @@ end @recipe function f(csg::AbstractConstructiveGeometry{T}, primitive::AbstractPrimitive{T}, axis::Symbol, slice::T, st::Symbol; n_samples = 40, CSG_scale = missing, projection = :none, full_det = false, linewidth = :auto) where {T} if st == :samplesurface + label --> l + seriesalpha --> 0.2 + seriescolor --> 1 + markerstrokewidth --> 0 + markersize --> 3 + if occursin("GRBackend", string(typeof(plotattributes[:plot_object].backend))) + aspect_ratio --> 1.0 + end + seriestype := :scatter CSG_scale = ismissing(CSG_scale) ? get_scale(csg) : CSG_scale spacing = T(CSG_scale/n_samples) pts_nounits = filter(p -> in(p,csg), sample(primitive, spacing)) diff --git a/src/ConstructiveSolidGeometry/plotting/SurfacePrimitives/SurfacePrimitives.jl b/src/ConstructiveSolidGeometry/plotting/SurfacePrimitives/SurfacePrimitives.jl index d8c3a86cc..61f527632 100644 --- a/src/ConstructiveSolidGeometry/plotting/SurfacePrimitives/SurfacePrimitives.jl +++ b/src/ConstructiveSolidGeometry/plotting/SurfacePrimitives/SurfacePrimitives.jl @@ -39,7 +39,14 @@ end mesh(s, n_arc, n_vert_lines) elseif st == :samplesurface label --> l - seriesalpha --> 0.4 + seriesalpha --> 0.2 + seriescolor --> 1 + markerstrokewidth --> 0 + markersize --> 3 + if occursin("GRBackend", string(typeof(plotattributes[:plot_object].backend))) + aspect_ratio --> 1.0 + end + seriestype := :scatter sample(s, extremum(s)/n_samples) elseif st == :slice @warn ":slice is not supported for AbstractSurfacePrimitive. Use :csg, :wireframe, :mesh3d, or :samplesurface." diff --git a/src/ConstructiveSolidGeometry/plotting/VolumePrimitives/VolumePrimitives.jl b/src/ConstructiveSolidGeometry/plotting/VolumePrimitives/VolumePrimitives.jl index a7a638a81..894f59ac5 100644 --- a/src/ConstructiveSolidGeometry/plotting/VolumePrimitives/VolumePrimitives.jl +++ b/src/ConstructiveSolidGeometry/plotting/VolumePrimitives/VolumePrimitives.jl @@ -31,6 +31,15 @@ end @recipe function f(p::AbstractPrimitive{T}, axis::Symbol, slice::T, st::Symbol; n_samples = 40, CSG_scale = missing, projection = :none, full_det = false, linewidth = :auto) where {T} if st == :samplesurface + label --> l + seriesalpha --> 0.2 + seriescolor --> 1 + markerstrokewidth --> 0 + markersize --> 3 + if occursin("GRBackend", string(typeof(plotattributes[:plot_object].backend))) + aspect_ratio --> 1.0 + end + seriestype := :scatter CSG_scale = ismissing(CSG_scale) ? extremum(p) : CSG_scale spacing = T(CSG_scale/n_samples) sample(p, spacing) From c85a13c8690c394c7dee2b77d9cbc7b33b08acf6 Mon Sep 17 00:00:00 2001 From: Felix Hagemann Date: Wed, 5 Nov 2025 07:32:14 -0800 Subject: [PATCH 3/9] Add `barycenter` for `Vector{<:AbstractCoordinatePoints}` --- .../PointsAndVectors/Points.jl | 35 +++++++++++-------- 1 file changed, 20 insertions(+), 15 deletions(-) diff --git a/src/ConstructiveSolidGeometry/PointsAndVectors/Points.jl b/src/ConstructiveSolidGeometry/PointsAndVectors/Points.jl index 7426c7298..d23256296 100644 --- a/src/ConstructiveSolidGeometry/PointsAndVectors/Points.jl +++ b/src/ConstructiveSolidGeometry/PointsAndVectors/Points.jl @@ -125,10 +125,8 @@ end @inline Base.adjoint(pt::CartesianPoint) = CartesianPoint(adjoint(pt.x), adjoint(pt.y), adjoint(pt.z)) # Barycentric combination -function Statistics.mean(A::AbstractArray{<:CartesianPoint{T}}; dims = :) where T - # Equivalent to A .- Ref(cartesian_zero), but allocation-free: - B = reinterpret(CartesianVector{T}, A) # ToDo: Avoid this memory allocation. - return cartesian_zero + mean(B; dims = dims) +function Statistics.mean(A::AbstractArray{<:CartesianPoint}; dims = :) + barycenter(A; dims) end function Statistics.mean(A::AbstractArray{<:CartesianPoint{T}}, weights::StatsBase.AbstractWeights; dims = :) where T @@ -137,6 +135,12 @@ function Statistics.mean(A::AbstractArray{<:CartesianPoint{T}}, weights::StatsBa return cartesian_zero + mean(B, weights; dims = dims) end +barycenter(X; kwargs...) = _barycenter_impl(eltype(X), X; kwargs...) +function _barycenter_impl(::Type{T}, points; kwargs...) where {T<:AbstractCoordinatePoint} + mean_vector = mean(Base.Fix2(-, cartesian_zero), CartesianPoint.(points); kwargs...) + return convert(T, cartesian_zero + mean_vector)::T +end + # Base.broadcastable(p::CartesianPoint) = (p.x, p.y, p.z) @@ -275,6 +279,14 @@ function Base.convert(::Type{CylindricalPoint{T}}, pt::CylindricalPoint{U}) wher return CylindricalPoint{T}(convert(T, pt.r), convert(T, pt.φ), convert(T, pt.z)) end +function Base.convert(::Type{CylindricalPoint{T}}, pt::AbstractCoordinatePoint)::CylindricalPoint{T} where {T} + return CylindricalPoint(pt) +end + +function Base.convert(::Type{CartesianPoint{T}}, pt::AbstractCoordinatePoint)::CartesianPoint{T} where {T} + return CartesianPoint(pt) +end + AbstractCoordinatePoint{T, Cylindrical}(r::Real, φ::Real, z::Real) where T = CylindricalPoint{T}(r, φ, z) Base.:(==)(a::CylindricalPoint, b::CylindricalPoint) = a.r == b.r && a.φ == b.φ && a.z == b.z @@ -298,19 +310,12 @@ Base.iszero(pt::CylindricalPoint) = iszero(pt.r) && iszero(pt.z) @inline Base.adjoint(pt::CylindricalPoint) = CylindricalPoint(adjoint(pt.r), adjoint(pt.φ), adjoint(pt.z)) # Barycentric combination -function Statistics.mean(A::AbstractArray{<:CylindricalPoint}; dims = :) - CylindricalPoint(mean(CartesianPoint, A; dims = dims)) +function Statistics.mean(A::AbstractArray{<:CylindricalPoint}; dims = :) + B = CartesianPoint.(A) + CylindricalPoint(mean(B; dims = dims)) end function Statistics.mean(A::AbstractArray{<:CylindricalPoint}, weights::StatsBase.AbstractWeights; dims = :) B = CartesianPoint.(A) CylindricalPoint(mean(B, weights; dims = dims)) -end - - -# function _Δφ(φ1::T, φ2::T)::T where {T} -# δφ = mod(φ2 - φ1, T(2π)) -# min(δφ, T(2π) - δφ) -# end -# -# _φNear(φ::Real, φMin::T, φMax::T) where {T} = _Δφ(T(φ),φMin) ≤ _Δφ(T(φ),φMax) ? φMin : φMax +end \ No newline at end of file From 7512e9b7135157f171156537d64d366208d4822e Mon Sep 17 00:00:00 2001 From: Felix Hagemann Date: Wed, 5 Nov 2025 07:32:31 -0800 Subject: [PATCH 4/9] Add tests for `barycenter` and `mean` Tests for `mean` are broken when using `SVectors` --- .../CSG_coordinates.jl | 20 ++++++++++++++++++- 1 file changed, 19 insertions(+), 1 deletion(-) diff --git a/test/ConstructiveSolidGeometry/CSG_coordinates.jl b/test/ConstructiveSolidGeometry/CSG_coordinates.jl index cb8104953..bcfeb42db 100644 --- a/test/ConstructiveSolidGeometry/CSG_coordinates.jl +++ b/test/ConstructiveSolidGeometry/CSG_coordinates.jl @@ -3,9 +3,11 @@ using Test using SolidStateDetectors.ConstructiveSolidGeometry: CartesianPoint, CartesianVector, - CartesianZero, cartesian_zero, CylindricalPoint, LocalAffineFrame, global_frame, frame_transformation + CartesianZero, cartesian_zero, CylindricalPoint, LocalAffineFrame, global_frame, frame_transformation, barycenter using StaticArrays: Size, SVector, SMatrix using InverseFunctions: inverse +using Statistics: mean + import Unitful @testset "points_and_vectors" begin @@ -44,6 +46,14 @@ import Unitful @test @inferred(CartesianZero{Float32}() * Unitful.u"mm") === CartesianZero{typeof(zero(Float32) * Unitful.u"mm")}() + + A = [CartesianPoint{Float32}(x,0,0) for x in -2:2] + @test isapprox(mean(A), CartesianPoint{Float32}(0,0,0)) + @test isapprox(barycenter(A), CartesianPoint{Float32}(0,0,0)) + + S = SVector{length(A)}(A) + @test_broken isapprox(mean(S), CartesianPoint{Float32}(0,0,0)) + @test isapprox(barycenter(S), CartesianPoint{Float32}(0,0,0)) # ToDo: Add more tests! end @@ -68,5 +78,13 @@ import Unitful @test @inferred(CylindricalPoint(a[1], a[2], a[3])) === a @test @inferred(CylindricalPoint(a[1], a[2], a[3])) == a @test @inferred(CylindricalPoint(a[1], a[2], a[3])) ≈ a + + A = [CylindricalPoint{Float32}(x,0,0) for x in -2:2] + @test isapprox(mean(A), CylindricalPoint{Float32}(0,0,0)) + @test isapprox(barycenter(A), CylindricalPoint{Float32}(0,0,0)) + + S = SVector{length(A)}(A) + @test_broken isapprox(mean(S), CylindricalPoint{Float32}(0,0,0)) + @test isapprox(barycenter(S), CylindricalPoint{Float32}(0,0,0)) end end From eee740728a27035f029cb08135dc6d4a2567ceb0 Mon Sep 17 00:00:00 2001 From: Felix Hagemann Date: Tue, 4 Nov 2025 17:30:37 -0800 Subject: [PATCH 5/9] Fix `extremum` for `Polygon` --- src/ConstructiveSolidGeometry/SurfacePrimitives/Polygon.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ConstructiveSolidGeometry/SurfacePrimitives/Polygon.jl b/src/ConstructiveSolidGeometry/SurfacePrimitives/Polygon.jl index 39a34b1f5..40c3b04fd 100644 --- a/src/ConstructiveSolidGeometry/SurfacePrimitives/Polygon.jl +++ b/src/ConstructiveSolidGeometry/SurfacePrimitives/Polygon.jl @@ -55,7 +55,7 @@ function sample(p::Polygon{N,T}, spacing::T)::Vector{CartesianPoint{T}} where {N end function extremum(p::Polygon{N,T})::T where {N,T} - c = sum(p.points)/N + c = barycenter(p.points) m = maximum([norm(point - c) for point in p.points]) end From 050aeccbd09b6510f8b764db15289d8f68a3a884 Mon Sep 17 00:00:00 2001 From: Felix Hagemann Date: Fri, 7 Nov 2025 09:10:32 -0800 Subject: [PATCH 6/9] Remove support for `Statistics.mean(::Vector{<:AbstractCoordinatePoint})` --- src/ChargeClustering/ChargeClustering.jl | 4 +-- .../PointsAndVectors/Points.jl | 35 ++++++------------- .../plotting/SurfacePrimitives/Polygon.jl | 2 +- src/Event/Event.jl | 4 +-- .../CSG_coordinates.jl | 6 ---- 5 files changed, 15 insertions(+), 36 deletions(-) diff --git a/src/ChargeClustering/ChargeClustering.jl b/src/ChargeClustering/ChargeClustering.jl index da37e161b..823c8edb6 100644 --- a/src/ChargeClustering/ChargeClustering.jl +++ b/src/ChargeClustering/ChargeClustering.jl @@ -39,10 +39,10 @@ function cluster_detector_hits( push!(r_edep, esum_u) esum = ustrip(esum_u) if esum ≈ 0 - push!(r_pos, mean(c_hits.pos)) + push!(r_pos, barycenter(c_hits.pos)) else weights = ustrip.(c_hits.edep) .* inv(esum) - push!(r_pos, sum(c_hits.pos .* weights)) + push!(r_pos, barycenter(c_hits.pos, weights)) end end else diff --git a/src/ConstructiveSolidGeometry/PointsAndVectors/Points.jl b/src/ConstructiveSolidGeometry/PointsAndVectors/Points.jl index d23256296..442d5cecf 100644 --- a/src/ConstructiveSolidGeometry/PointsAndVectors/Points.jl +++ b/src/ConstructiveSolidGeometry/PointsAndVectors/Points.jl @@ -124,25 +124,21 @@ end @inline Base.transpose(pt::CartesianPoint) = CartesianPoint(transpose(pt.x), transpose(pt.y), transpose(pt.z)) @inline Base.adjoint(pt::CartesianPoint) = CartesianPoint(adjoint(pt.x), adjoint(pt.y), adjoint(pt.z)) -# Barycentric combination -function Statistics.mean(A::AbstractArray{<:CartesianPoint}; dims = :) - barycenter(A; dims) -end - -function Statistics.mean(A::AbstractArray{<:CartesianPoint{T}}, weights::StatsBase.AbstractWeights; dims = :) where T - # Equivalent to A .- Ref(cartesian_zero), but allocation-free: - B = reinterpret(CartesianVector{T}, A) # ToDo: Avoid this memory allocation. - return cartesian_zero + mean(B, weights; dims = dims) -end - -barycenter(X; kwargs...) = _barycenter_impl(eltype(X), X; kwargs...) +barycenter(X, args...; kwargs...) = _barycenter_impl(eltype(X), X, args...; kwargs...) function _barycenter_impl(::Type{T}, points; kwargs...) where {T<:AbstractCoordinatePoint} mean_vector = mean(Base.Fix2(-, cartesian_zero), CartesianPoint.(points); kwargs...) return convert(T, cartesian_zero + mean_vector)::T end -# Base.broadcastable(p::CartesianPoint) = (p.x, p.y, p.z) +function _barycenter_impl(::Type{T}, points, weights::StatsBase.AbstractWeights; kwargs...) where {T<:AbstractCoordinatePoint} + mean_vector = mean(CartesianPoint.(points) .- Ref(cartesian_zero), weights; kwargs...) + return convert(T, cartesian_zero + mean_vector)::T +end +# Fallback if a Vector of Vectors (e.g. StaticVectors) is passed to barycenter +function _barycenter_impl(::Type{T}, args...; kwargs...) where {T <: AbstractVector} + convert(T, mean(args...; kwargs...)) +end """ @@ -307,15 +303,4 @@ Base.iszero(pt::CylindricalPoint) = iszero(pt.r) && iszero(pt.z) @inline Base.:(-)(a::CylindricalPoint, b::CylindricalPoint) = CartesianPoint(a) - CartesianPoint(b) @inline Base.transpose(pt::CylindricalPoint) = CylindricalPoint(transpose(pt.r), transpose(pt.φ), transpose(pt.z)) -@inline Base.adjoint(pt::CylindricalPoint) = CylindricalPoint(adjoint(pt.r), adjoint(pt.φ), adjoint(pt.z)) - -# Barycentric combination -function Statistics.mean(A::AbstractArray{<:CylindricalPoint}; dims = :) - B = CartesianPoint.(A) - CylindricalPoint(mean(B; dims = dims)) -end - -function Statistics.mean(A::AbstractArray{<:CylindricalPoint}, weights::StatsBase.AbstractWeights; dims = :) - B = CartesianPoint.(A) - CylindricalPoint(mean(B, weights; dims = dims)) -end \ No newline at end of file +@inline Base.adjoint(pt::CylindricalPoint) = CylindricalPoint(adjoint(pt.r), adjoint(pt.φ), adjoint(pt.z)) \ No newline at end of file diff --git a/src/ConstructiveSolidGeometry/plotting/SurfacePrimitives/Polygon.jl b/src/ConstructiveSolidGeometry/plotting/SurfacePrimitives/Polygon.jl index b3a45ad7d..e85409f60 100644 --- a/src/ConstructiveSolidGeometry/plotting/SurfacePrimitives/Polygon.jl +++ b/src/ConstructiveSolidGeometry/plotting/SurfacePrimitives/Polygon.jl @@ -19,7 +19,7 @@ @series begin label := nothing seriestype := :vector - mean(p.points), Plane(p).normal / 5 + barycenter(p.points), Plane(p).normal / 5 end end end \ No newline at end of file diff --git a/src/Event/Event.jl b/src/Event/Event.jl index c0338565a..0ec73403d 100644 --- a/src/Event/Event.jl +++ b/src/Event/Event.jl @@ -144,7 +144,7 @@ function move_charges_inside_semiconductor!( idx_in = broadcast( pt -> pt in det.semiconductor, locations[n]); if !all(idx_in) edep_weights = ustrip.(energies[n]) - charge_center = mean(locations[n], Weights(edep_weights)) + charge_center = barycenter(locations[n], Weights(edep_weights)) @assert charge_center in det.semiconductor "The center of the charge cloud ($(charge_center)) is not inside the semiconductor." surf = ConstructiveSolidGeometry.surfaces(det.semiconductor.geometry) for (k,m) in enumerate(findall(.!idx_in)) @@ -162,7 +162,7 @@ function move_charges_inside_semiconductor!( end end end - charge_center_new = mean(locations[n], Weights(edep_weights)) + charge_center_new = barycenter(locations[n], Weights(edep_weights)) if verbose @warn "$(sum(.!idx_in)) charges of the charge cloud at $(round.(charge_center - cartesian_zero, digits = (T == Float64 ? 12 : 6)))"* " are outside. Moving them inside...\nThe new charge center is at $(round.(charge_center_new - cartesian_zero, digits = (T == Float64 ? 12 : 6))).\n" diff --git a/test/ConstructiveSolidGeometry/CSG_coordinates.jl b/test/ConstructiveSolidGeometry/CSG_coordinates.jl index bcfeb42db..38078c7df 100644 --- a/test/ConstructiveSolidGeometry/CSG_coordinates.jl +++ b/test/ConstructiveSolidGeometry/CSG_coordinates.jl @@ -6,7 +6,6 @@ using SolidStateDetectors.ConstructiveSolidGeometry: CartesianPoint, CartesianVe CartesianZero, cartesian_zero, CylindricalPoint, LocalAffineFrame, global_frame, frame_transformation, barycenter using StaticArrays: Size, SVector, SMatrix using InverseFunctions: inverse -using Statistics: mean import Unitful @@ -48,13 +47,10 @@ import Unitful @test @inferred(CartesianZero{Float32}() * Unitful.u"mm") === CartesianZero{typeof(zero(Float32) * Unitful.u"mm")}() A = [CartesianPoint{Float32}(x,0,0) for x in -2:2] - @test isapprox(mean(A), CartesianPoint{Float32}(0,0,0)) @test isapprox(barycenter(A), CartesianPoint{Float32}(0,0,0)) S = SVector{length(A)}(A) - @test_broken isapprox(mean(S), CartesianPoint{Float32}(0,0,0)) @test isapprox(barycenter(S), CartesianPoint{Float32}(0,0,0)) - # ToDo: Add more tests! end @testset "cylindrical" begin @@ -80,11 +76,9 @@ import Unitful @test @inferred(CylindricalPoint(a[1], a[2], a[3])) ≈ a A = [CylindricalPoint{Float32}(x,0,0) for x in -2:2] - @test isapprox(mean(A), CylindricalPoint{Float32}(0,0,0)) @test isapprox(barycenter(A), CylindricalPoint{Float32}(0,0,0)) S = SVector{length(A)}(A) - @test_broken isapprox(mean(S), CylindricalPoint{Float32}(0,0,0)) @test isapprox(barycenter(S), CylindricalPoint{Float32}(0,0,0)) end end From 1fb4479aa6f26ee6bfb7181150d690b29d8122fa Mon Sep 17 00:00:00 2001 From: Felix Hagemann Date: Fri, 7 Nov 2025 09:11:12 -0800 Subject: [PATCH 7/9] Import `barycenter` and `distance_squared` into SSD --- .../PointsAndVectors/PointsAndVectors.jl | 1 - src/SolidStateDetectors.jl | 11 ++++++----- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/ConstructiveSolidGeometry/PointsAndVectors/PointsAndVectors.jl b/src/ConstructiveSolidGeometry/PointsAndVectors/PointsAndVectors.jl index 6f57f1d2e..ad5697c5a 100644 --- a/src/ConstructiveSolidGeometry/PointsAndVectors/PointsAndVectors.jl +++ b/src/ConstructiveSolidGeometry/PointsAndVectors/PointsAndVectors.jl @@ -18,4 +18,3 @@ include("AffineFrames.jl") @inline distance_squared(v::CartesianVector) = v.x * v.x + v.y * v.y + v.z * v.z @inline distance_squared(p1::CartesianPoint{T}, p2::CartesianPoint{T}) where {T <: Real} = distance_squared(p1 - p2) -export distance_squared diff --git a/src/SolidStateDetectors.jl b/src/SolidStateDetectors.jl index 80052f9b9..4841a2a81 100644 --- a/src/SolidStateDetectors.jl +++ b/src/SolidStateDetectors.jl @@ -41,12 +41,13 @@ using .ConstructiveSolidGeometry: Cartesian, Cylindrical, AbstractCoordinateSystem, CoordinateSystemType, Geometry, AbstractGeometry, AbstractSurfacePrimitive, parse_rotation_matrix, parse_translate_vector, parse_CSG_transformation, - transform, CSG_dict, Transformations, combine_transformations, - ConfigFileError, _parse_value, - LengthQuantity, AngleQuantity, get_scale, - LocalAffineFrame, cartesian_zero, global_frame, frame_transformation - + transform, CSG_dict, Transformations, combine_transformations, barycenter, + LocalAffineFrame, cartesian_zero, global_frame, frame_transformation, + ConfigFileError, _parse_value, distance_squared, + LengthQuantity, AngleQuantity, get_scale + import .ConstructiveSolidGeometry: sample, to_internal_units, from_internal_units + export CartesianPoint, CartesianVector, CartesianZero, cartesian_zero, CylindricalPoint import Clustering From d6a62f381b3e134060519d44b05a54be2a50fe72 Mon Sep 17 00:00:00 2001 From: Felix Hagemann Date: Tue, 11 Nov 2025 02:41:26 -0800 Subject: [PATCH 8/9] Add test for `cluster_detector_hits` --- src/ChargeClustering/ChargeClustering.jl | 39 +++++++++++------------- test/test_geant4.jl | 21 +++++++++++++ 2 files changed, 39 insertions(+), 21 deletions(-) diff --git a/src/ChargeClustering/ChargeClustering.jl b/src/ChargeClustering/ChargeClustering.jl index 823c8edb6..4559db170 100644 --- a/src/ChargeClustering/ChargeClustering.jl +++ b/src/ChargeClustering/ChargeClustering.jl @@ -1,20 +1,16 @@ # # This file is a part of SolidStateDetectors.jl, licensed under the MIT License (MIT). -# breaking changes in Clustering v0.14 -> v0.15 -@inline _get_clusters(clusters::Clustering.DbscanResult) = clusters.clusters -@inline _get_clusters(clusters::Vector{Clustering.DbscanCluster}) = clusters - function cluster_detector_hits( - detno::AbstractVector{<:Integer}, - edep::AbstractVector{TT}, - pos::AbstractVector{<:StaticVector{3,PT}}, - cluster_radius::RealQuantity -) where {TT<:RealQuantity, PT <: RealQuantity} - Table = TypedTables.Table - unsorted = Table(detno = detno, edep = edep, pos = pos) + detno::AbstractVector{<:Integer}, + edep::AbstractVector{TT}, + pos::AbstractVector{<:Union{<:StaticVector{3,PT}, <:CartesianPoint{PT}}}, + cluster_radius::RealQuantity + ) where {TT <: RealQuantity, PT <: RealQuantity} + + unsorted = TypedTables.Table(detno = detno, edep = edep, pos = pos) sorting_idxs = sortperm(unsorted.detno) sorted = unsorted[sorting_idxs] - grouped = Table(consgroupedview(sorted.detno, TypedTables.columns(sorted))) + grouped = TypedTables.Table(consgroupedview(sorted.detno, TypedTables.columns(sorted))) r_detno = similar(detno, 0) r_edep = similar(edep, 0) @@ -24,25 +20,26 @@ function cluster_detector_hits( ustripped_cradius = ustrip(posunit, cluster_radius) for d_hits_nt in grouped - d_hits = Table(d_hits_nt) + d_hits = TypedTables.Table(d_hits_nt) d_detno = first(d_hits.detno) @assert all(isequal(d_detno), d_hits.detno) if length(d_hits) > 3 - clusters = Clustering.dbscan(ustrip.(flatview(d_hits.pos)), ustripped_cradius, leafsize = 20, min_neighbors = 1, min_cluster_size = 1) - for c in _get_clusters(clusters) + + clusters = Clustering.dbscan(hcat((ustrip.(getindex.(d_hits.pos,i)) for i in 1:3)...)', + ustripped_cradius, leafsize = 20, min_neighbors = 1, min_cluster_size = 1).clusters + for c in clusters idxs = vcat(c.boundary_indices, c.core_indices) @assert length(idxs) == c.size c_hits = view(d_hits, idxs) push!(r_detno, d_detno) - esum_u = sum(c_hits.edep) - push!(r_edep, esum_u) - esum = ustrip(esum_u) - if esum ≈ 0 + esum = sum(c_hits.edep) + push!(r_edep, esum) + if esum ≈ zero(TT) push!(r_pos, barycenter(c_hits.pos)) else - weights = ustrip.(c_hits.edep) .* inv(esum) - push!(r_pos, barycenter(c_hits.pos, weights)) + weights = ustrip.(Unitful.NoUnits, c_hits.edep .* inv(esum)) + push!(r_pos, barycenter(c_hits.pos, StatsBase.Weights(weights))) end end else diff --git a/test/test_geant4.jl b/test/test_geant4.jl index 22fcef0f3..6444c6a76 100644 --- a/test/test_geant4.jl +++ b/test/test_geant4.jl @@ -3,7 +3,9 @@ using Test using SolidStateDetectors +using ArraysOfArrays using RadiationDetectorSignals +using StaticArrays using TypedTables using Unitful @@ -66,6 +68,12 @@ end @test evts isa Table @test length(evts) == 100 + # Cluster events by radius + clustered_evts = @test_nowarn SolidStateDetectors.cluster_detector_hits(evts, 10u"µm") + @test length(clustered_evts) == length(evts) + @test length(flatview(clustered_evts.pos)) <= length(flatview(evts.pos)) + @test eltype(first(clustered_evts.pos)) <: CartesianPoint + # Generate waveforms simulate!(sim, refinement_limits = [0.2,0.1,0.05,0.03,0.02]) wf = simulate_waveforms(evts, sim, Δt = 1u"ns", max_nsteps = 2000) @@ -73,4 +81,17 @@ end @test :waveform in columnnames(wf) @test length(wf) == length(evts) * sum(.!ismissing.(sim.weighting_potentials)) + # Try the same method using StaticVectors as eltype of pos + evts_static = Table(evts; pos = VectorOfVectors(broadcast.(p -> SVector{3}(p.x, p.y, p.z), evts.pos))) + clustered_evts_static = @test_nowarn SolidStateDetectors.cluster_detector_hits(evts_static, 10u"µm") + @test length(clustered_evts_static) == length(evts_static) + @test length(flatview(clustered_evts_static.pos)) <= length(flatview(evts_static.pos)) + @test eltype(first(clustered_evts_static.pos)) <: StaticVector{3} + + # Generate waveforms using StaticVectors as eltype of pos + wf_static = simulate_waveforms(evts_static, sim, Δt = 1u"ns", max_nsteps = 2000) + @test wf_static isa Table + @test :waveform in columnnames(wf_static) + @test length(wf_static) == length(evts_static) * sum(.!ismissing.(sim.weighting_potentials)) + end \ No newline at end of file From 99334b195bb8243b12b3218ee70fe920edf3387e Mon Sep 17 00:00:00 2001 From: Felix Hagemann Date: Tue, 11 Nov 2025 10:44:27 -0800 Subject: [PATCH 9/9] Fix downgrade tests on julia-1.10 --- .github/workflows/Downgrade.yml | 4 ++++ test/test_geant4.jl | 4 ++-- 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/.github/workflows/Downgrade.yml b/.github/workflows/Downgrade.yml index c32013fb9..573d82e9a 100644 --- a/.github/workflows/Downgrade.yml +++ b/.github/workflows/Downgrade.yml @@ -28,6 +28,10 @@ jobs: - name: Add LegendJuliaRegistry run: julia -e 'using Pkg; Pkg.Registry.add("General"); Pkg.Registry.add(Pkg.RegistrySpec(url = "https://github.com/legend-exp/LegendJuliaRegistry"))' shell: bash + - name: Downgrade NearestNeighbors # Temporary fix: this is needed to work with older versions of StaticArrays + run: julia --project=. -e 'using Pkg; Pkg.add(name = "NearestNeighbors", version = "0.4.13"); Pkg.compat("NearestNeighbors", "0.4.13")' + shell: bash + if: matrix.version == '1.10' - uses: julia-actions/julia-downgrade-compat@v2 with: julia_version: ${{ matrix.version }} diff --git a/test/test_geant4.jl b/test/test_geant4.jl index 6444c6a76..a717e1ad7 100644 --- a/test/test_geant4.jl +++ b/test/test_geant4.jl @@ -69,7 +69,7 @@ end @test length(evts) == 100 # Cluster events by radius - clustered_evts = @test_nowarn SolidStateDetectors.cluster_detector_hits(evts, 10u"µm") + clustered_evts = SolidStateDetectors.cluster_detector_hits(evts, 10u"µm") @test length(clustered_evts) == length(evts) @test length(flatview(clustered_evts.pos)) <= length(flatview(evts.pos)) @test eltype(first(clustered_evts.pos)) <: CartesianPoint @@ -83,7 +83,7 @@ end # Try the same method using StaticVectors as eltype of pos evts_static = Table(evts; pos = VectorOfVectors(broadcast.(p -> SVector{3}(p.x, p.y, p.z), evts.pos))) - clustered_evts_static = @test_nowarn SolidStateDetectors.cluster_detector_hits(evts_static, 10u"µm") + clustered_evts_static = SolidStateDetectors.cluster_detector_hits(evts_static, 10u"µm") @test length(clustered_evts_static) == length(evts_static) @test length(flatview(clustered_evts_static.pos)) <= length(flatview(evts_static.pos)) @test eltype(first(clustered_evts_static.pos)) <: StaticVector{3}