Add test for cluster_detector_hits

Author: fhagemann

src/ChargeClustering/ChargeClustering.jl

@@ -7,14 +7,14 @@
 function cluster_detector_hits(
     detno::AbstractVector{<:Integer},
     edep::AbstractVector{TT},
-    pos::AbstractVector{<:StaticVector{3,PT}},
+    pos::AbstractVector{<:Union{<:StaticVector{PT}, <:CartesianPoint{PT}}},
     cluster_radius::RealQuantity
 ) where {TT<:RealQuantity, PT <: RealQuantity}
-    Table = TypedTables.Table
-    unsorted = Table(detno = detno, edep = edep, pos = pos)
+
+    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 +24,25 @@ 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)
+            
+            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)
             for c in _get_clusters(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

test/test_geant4.jl

@@ -3,6 +3,7 @@
 using Test
 
 using SolidStateDetectors
+using ArraysOfArrays
 using RadiationDetectorSignals
 using TypedTables
 using Unitful
@@ -66,6 +67,11 @@ 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))
+
     # 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)