Skip to content

Commit 32851dc

Browse files
fhagemannfhagemann
authored
Merge pull request #461 from JuliaPhysics/group_hits
Allow automatic grouping of hits based on a `group_distance` when creating an `Event`
2 parents 03320b4 + cacc71f commit 32851dc

File tree

9 files changed

+183
-25
lines changed

9 files changed

+183
-25
lines changed

‎Project.toml‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -6,6 +6,7 @@ version = "0.10.12"
66
Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
77
ArraysOfArrays = "65a8f2f4-9b39-5baf-92e2-a9cc46fdf018"
88
Clustering = "aaaa29a8-35af-508c-8bc3-b662a17a0fe5"
9+
DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
910
Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
1011
FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
1112
Format = "1fa38f19-a742-5d3f-a2b9-30dd87b9d5f8"
@@ -48,6 +49,7 @@ SolidStateDetectorsGeant4Ext = "Geant4"
4849
Adapt = "3, 4"
4950
ArraysOfArrays = "0.5, 0.6"
5051
Clustering = "0.15"
52+
DataStructures = "0.17, 0.18"
5153
Distributions = "0.24.5, 0.25"
5254
FillArrays = "0.9, 0.10, 0.11, 0.12, 0.13, 1"
5355
Format = "1.2"

‎ext/SolidStateDetectorsLegendHDF5IOExt.jl‎

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -5,7 +5,7 @@ module SolidStateDetectorsLegendHDF5IOExt
55
import ..LegendHDF5IO
66

77
using SolidStateDetectors
8-
using SolidStateDetectors: RealQuantity, SSDFloat, to_internal_units, chunked_ranges
8+
using SolidStateDetectors: RealQuantity, SSDFloat, to_internal_units, chunked_ranges, LengthQuantity
99
using TypedTables, Unitful
1010
using Format
1111

@@ -34,6 +34,7 @@ function SolidStateDetectors.simulate_waveforms( mcevents::TypedTables.Table, si
3434
self_repulsion::Bool = false,
3535
number_of_carriers::Int = 1,
3636
number_of_shells::Int = 1,
37+
max_interaction_distance::Union{<:Real, <:LengthQuantity} = NaN,
3738
verbose = false) where {T <: SSDFloat}
3839
n_total_physics_events = length(mcevents)
3940
Δtime = T(to_internal_units(Δt))
@@ -48,7 +49,7 @@ function SolidStateDetectors.simulate_waveforms( mcevents::TypedTables.Table, si
4849
for evtrange in evt_ranges
4950
ofn = joinpath(output_dir, "$(output_base_name)_evts_$(nfmt(first(evtrange)))-$(nfmt(last(evtrange))).h5")
5051
@info "Now simulating $(evtrange) and storing it in\n\t \"$ofn\""
51-
mcevents_sub = simulate_waveforms(mcevents[evtrange], sim; Δt, max_nsteps, diffusion, self_repulsion, number_of_carriers, number_of_shells, verbose)
52+
mcevents_sub = simulate_waveforms(mcevents[evtrange], sim; Δt, max_nsteps, diffusion, self_repulsion, number_of_carriers, number_of_shells, max_interaction_distance, verbose)
5253

5354
LegendHDF5IO.lh5open(ofn, "w") do h5f
5455
LegendHDF5IO.writedata(h5f.data_store, "generated_waveforms", mcevents_sub)

‎src/ChargeClustering/ChargeClustering.jl‎

Lines changed: 54 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -72,3 +72,57 @@ function cluster_detector_hits(table::TypedTables.Table, cluster_radius::RealQua
7272
)
7373
))
7474
end
75+
76+
77+
function _group_points_by_distance(pts::AbstractVector{CartesianPoint{T}}, group_distance::T)::Tuple{Vector{Int}, Vector{Int}} where {T}
78+
79+
n::Int = length(pts)
80+
81+
# use BFS to find connected components
82+
visited = falses(n)
83+
clustersidx = similar(eachindex(pts))
84+
elem_ptr = similar(eachindex(pts), n+1)
85+
queue = DataStructures.CircularBuffer{Int}(n)
86+
87+
counter = firstindex(pts)
88+
Cidx = firstindex(elem_ptr)
89+
elem_ptr[Cidx] = counter
90+
91+
@inbounds for start in eachindex(pts)
92+
if !visited[start]
93+
push!(queue, start)
94+
visited[start] = true
95+
while !isempty(queue)
96+
node = popfirst!(queue)
97+
clustersidx[counter] = node
98+
counter += 1
99+
for neighbor in eachindex(pts)
100+
if !visited[neighbor] && distance_squared(pts[node], pts[neighbor]) <= group_distance^2
101+
push!(queue, neighbor)
102+
visited[neighbor] = true
103+
end
104+
end
105+
end
106+
Cidx += 1
107+
elem_ptr[Cidx] = counter
108+
end
109+
end
110+
clustersidx, elem_ptr[begin:Cidx]
111+
end
112+
113+
function group_points_by_distance(pts::AbstractVector{CartesianPoint{T}}, group_distance::T)::Tuple{VectorOfVectors{CartesianPoint{T}}, VectorOfVectors{T}} where {T <: SSDFloat}
114+
clustersidx, elem_ptr = _group_points_by_distance(pts, group_distance)
115+
VectorOfVectors(pts[clustersidx], elem_ptr)
116+
end
117+
118+
function group_points_by_distance(pts::AbstractVector{CartesianPoint{T}}, energies::AbstractVector{T}, group_distance::T)::Tuple{VectorOfVectors{CartesianPoint{T}}, VectorOfVectors{T}} where {T <: SSDFloat}
119+
@assert eachindex(pts) == eachindex(energies)
120+
clustersidx, elem_ptr = _group_points_by_distance(pts, group_distance)
121+
VectorOfVectors(pts[clustersidx], elem_ptr), VectorOfVectors(energies[clustersidx], elem_ptr)
122+
end
123+
124+
function group_points_by_distance(pts::AbstractVector{CartesianPoint{T}}, energies::AbstractVector{T}, radius::AbstractVector{T}, group_distance::T)::Tuple{VectorOfVectors{CartesianPoint{T}}, VectorOfVectors{T}, VectorOfVectors{T}} where {T <: SSDFloat}
125+
@assert eachindex(pts) == eachindex(energies) == eachindex(radius)
126+
clustersidx, elem_ptr = _group_points_by_distance(pts, group_distance)
127+
VectorOfVectors(pts[clustersidx], elem_ptr), VectorOfVectors(energies[clustersidx], elem_ptr), VectorOfVectors(radius[clustersidx], elem_ptr)
128+
end

‎src/ChargeDrift/ChargeDrift.jl‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -130,7 +130,7 @@ function _add_fieldvector_selfrepulsion!(step_vectors::Vector{CartesianVector{T}
130130
if iszero(direction) # if the two charges are at the exact same position
131131
continue # don't let them self-repel each other but treat them as same change
132132
end # if diffusion is simulated, they will very likely be separated in the next step
133-
tmp::T = elementary_charge * inv(4 * pi * ϵ0 * ϵ_r * max(sum(direction.^2), T(1e-10))) # minimum distance is 10µm
133+
tmp::T = elementary_charge * inv(4π * ϵ0 * ϵ_r * max(distance_squared(direction), T(1e-10))) # minimum distance is 10µm
134134
step_vectors[n] += charges[m] * tmp * normalize(direction)
135135
step_vectors[m] -= charges[n] * tmp * normalize(direction)
136136
end

‎src/ConstructiveSolidGeometry/PointsAndVectors/PointsAndVectors.jl‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -3,3 +3,7 @@ abstract type AbstractCoordinateVector{T, S} <: StaticArrays.FieldVector{3, T} e
33

44
include("Points.jl")
55
include("Vectors.jl")
6+
7+
@inline distance_squared(v::CartesianVector) = v.x * v.x + v.y * v.y + v.z * v.z
8+
@inline distance_squared(p1::CartesianPoint{T}, p2::CartesianPoint{T}) where {T <: Real} = distance_squared(CartesianVector(p1 .- p2))
9+
export distance_squared

‎src/Event/Event.jl‎

Lines changed: 42 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -26,10 +26,16 @@ function Event(location::AbstractCoordinatePoint{T}, energy::RealQuantity = one(
2626
return evt
2727
end
2828

29-
function Event(locations::Vector{<:AbstractCoordinatePoint{T}}, energies::Vector{<:RealQuantity} = ones(T, length(locations)))::Event{T} where {T <: SSDFloat}
29+
function Event(locations::Vector{<:AbstractCoordinatePoint{T}}, energies::Vector{<:RealQuantity} = ones(T, length(locations)); max_interaction_distance::Union{<:Real, <:LengthQuantity} = NaN)::Event{T} where {T <: SSDFloat}
30+
d::T = T(to_internal_units(max_interaction_distance))
31+
@assert isnan(d) || d >= 0 "Max. interaction distance must be positive or NaN (no grouping), but $(max_interaction_distance) was given."
3032
evt = Event{T}()
31-
evt.locations = VectorOfArrays(broadcast(pt -> [CartesianPoint(pt)], locations))
32-
evt.energies = VectorOfArrays(broadcast(E -> [T(to_internal_units(E))], energies))
33+
if isnan(d) # default: no grouping, the charges from different hits drift independently
34+
evt.locations = VectorOfArrays(broadcast(pt -> [CartesianPoint(pt)], locations))
35+
evt.energies = VectorOfArrays(broadcast(E -> [T(to_internal_units(E))], energies))
36+
else
37+
evt.locations, evt.energies = group_points_by_distance(CartesianPoint.(locations), T.(to_internal_units.(energies)), d)
38+
end
3339
return evt
3440
end
3541

@@ -42,6 +48,7 @@ end
4248

4349
function Event(nbcc::NBodyChargeCloud{T})::Event{T} where {T <: SSDFloat}
4450
evt = Event{T}()
51+
# charges in one NBodyChargeCloud should see each other by default
4552
evt.locations = VectorOfArrays([nbcc.locations])
4653
evt.energies = VectorOfArrays([nbcc.energies])
4754
return evt
@@ -55,15 +62,40 @@ function Event(nbccs::Vector{<:NBodyChargeCloud{T}})::Event{T} where {T <: SSDFl
5562
end
5663

5764
function Event(locations::Vector{<:AbstractCoordinatePoint{T}}, energies::Vector{<:RealQuantity}, N::Int;
58-
particle_type::Type{PT} = Gamma, number_of_shells::Int = 2,
59-
radius::Vector{<:RealQuantity} = radius_guess.(T.(to_internal_units.(energies)), particle_type)
60-
)::Event{T} where {T <: SSDFloat, PT <: ParticleType}
61-
65+
particle_type::Type{PT} = Gamma, number_of_shells::Int = 2,
66+
radius::Vector{<:Union{<:Real, <:LengthQuantity}} = radius_guess.(T.(to_internal_units.(energies)), particle_type),
67+
max_interaction_distance::Union{<:Real, <:LengthQuantity} = NaN
68+
)::Event{T} where {T <: SSDFloat, PT <: ParticleType}
69+
70+
6271
@assert eachindex(locations) == eachindex(energies) == eachindex(radius)
63-
return Event(broadcast(i ->
64-
NBodyChargeCloud(locations[i], energies[i], N, particle_type,
72+
73+
d::T = T(to_internal_units(max_interaction_distance))
74+
@assert isnan(d) || d >= 0 "Max. interaction distance must be positive or NaN (no grouping), but $(max_interaction_distance) was given."
75+
76+
if isnan(d) # default: no grouping, the charges from different hits drift independently
77+
Event(
78+
broadcast(i ->
79+
NBodyChargeCloud(locations[i], energies[i], N,
6580
radius = T(to_internal_units(radius[i])), number_of_shells = number_of_shells),
66-
eachindex(locations)))
81+
eachindex(locations))
82+
)
83+
else # otherwise: group the CENTERS by distance and compose the NBodyChargeClouds around them
84+
loc, ene, rad = group_points_by_distance(CartesianPoint.(locations), T.(to_internal_units.(energies)), T.(to_internal_units.(radius)), d)
85+
nbccs = broadcast(i ->
86+
broadcast(j ->
87+
let nbcc = NBodyChargeCloud(loc[i][j], ene[i][j], N,
88+
radius = rad[i][j], number_of_shells = number_of_shells)
89+
nbcc.locations, nbcc.energies
90+
end,
91+
eachindex(loc[i])),
92+
eachindex(loc))
93+
94+
Event(
95+
broadcast(nbcc -> vcat(getindex.(nbcc, 1)...), nbccs),
96+
broadcast(nbcc -> vcat(getindex.(nbcc, 2)...), nbccs)
97+
)
98+
end
6799
end
68100

69101
function Event(

‎src/MCEventsProcessing/MCEventsProcessing.jl‎

Lines changed: 35 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -27,6 +27,7 @@ If [LegendHDF5IO.jl](https://github.com/legend-exp/LegendHDF5IO.jl) is loaded, t
2727
* `self_repulsion::Bool = false`: Activate or deactive self-repulsion of charge carriers of the same type.
2828
* `number_of_carriers::Int = 1`: Number of charge carriers to be used in the N-Body simulation of an energy deposition.
2929
* `number_of_shells::Int = 1`: Number of shells around the `center` point of the energy deposition.
30+
* `max_interaction_distance = NaN`: Maximum distance for which charge clouds will interact with each other (if `NaN`, then all charge clouds drift independently).
3031
* `verbose = false`: Activate or deactivate additional info output.
3132
* `chunk_n_physics_events::Int = 1000` (`LegendHDF5IO` only): Number of events that should be saved in a single HDF5 output file.
3233
@@ -53,6 +54,7 @@ function simulate_waveforms( mcevents::TypedTables.Table, sim::Simulation{T};
5354
self_repulsion::Bool = false,
5455
number_of_carriers::Int = 1,
5556
number_of_shells::Int = 1,
57+
max_interaction_distance::Union{<:Real, <:LengthQuantity} = NaN,
5658
verbose::Bool = false ) where {T <: SSDFloat}
5759
n_total_physics_events = length(mcevents)
5860
Δtime = T(to_internal_units(Δt))
@@ -71,7 +73,7 @@ function simulate_waveforms( mcevents::TypedTables.Table, sim::Simulation{T};
7173

7274
# First simulate drift paths
7375
drift_paths_and_edeps = _simulate_charge_drifts(mcevents, sim, Δt, max_nsteps, electric_field,
74-
diffusion, self_repulsion, number_of_carriers, number_of_shells, verbose)
76+
diffusion, self_repulsion, number_of_carriers, number_of_shells, max_interaction_distance, verbose)
7577
drift_paths = map(x -> vcat([vcat(ed...) for ed in x[1]]...), drift_paths_and_edeps)
7678
edeps = map(x -> vcat([vcat(ed...) for ed in x[2]]...), drift_paths_and_edeps)
7779
# now iterate over contacts and generate the waveform for each contact
@@ -94,19 +96,40 @@ function simulate_waveforms( mcevents::TypedTables.Table, sim::Simulation{T};
9496
return vcat(mcevents_chns...)
9597
end
9698

97-
_convertEnergyDepsToChargeDeps(pos::AbstractVector{<:SVector{3}}, edep::AbstractVector{<:Quantity}, det::SolidStateDetector; kwargs...) =
98-
_convertEnergyDepsToChargeDeps([[p] for p in pos], [[e] for e in edep], det; kwargs...)
99+
function _convertEnergyDepsToChargeDeps(
100+
pos::AbstractVector{<:SVector{3}}, edep::AbstractVector{<:Quantity}, det::SolidStateDetector{T};
101+
particle_type::Type{PT} = Gamma,
102+
radius::Vector{<:Union{<:Real, <:LengthQuantity}} = radius_guess.(T.(to_internal_units.(edep)), particle_type),
103+
max_interaction_distance::Union{<:Real, <:LengthQuantity} = NaN,
104+
kwargs...
105+
) where {T <: SSDFloat, PT <: ParticleType}
106+
107+
@assert eachindex(pos) == eachindex(edep) == eachindex(radius)
108+
109+
d::T = T(to_internal_units(max_interaction_distance))
110+
@assert isnan(d) || d >= 0 "Max. interaction distance must be positive or NaN (no grouping), but $(max_interaction_distance) was given."
111+
112+
loc, ene, rad = group_points_by_distance(CartesianPoint.(map(x -> to_internal_units.(x), pos)), T.(to_internal_units.(edep)), T.(to_internal_units.(radius)), d)
113+
_convertEnergyDepsToChargeDeps(loc, ene, det; radius = rad, kwargs...)
114+
end
115+
116+
function _convertEnergyDepsToChargeDeps(
117+
pos::AbstractVector{<:AbstractVector}, edep::AbstractVector{<:AbstractVector}, det::SolidStateDetector{T};
118+
particle_type::Type{PT} = Gamma,
119+
radius::AbstractVector{<:AbstractVector{<:Union{<:Real, <:LengthQuantity}}} = map(e -> radius_guess.(to_internal_units.(e), particle_type), edep),
120+
number_of_carriers::Int = 1, number_of_shells::Int = 1,
121+
max_interaction_distance::Union{<:Real, <:LengthQuantity} = NaN # is ignored here
122+
) where {T <: SSDFloat, PT <: ParticleType}
99123

100-
function _convertEnergyDepsToChargeDeps(pos::AbstractVector{<:AbstractVector}, edep::AbstractVector{<:AbstractVector}, det::SolidStateDetector;
101-
number_of_carriers::Int = 1, number_of_shells::Int = 1)
102124
charge_clouds = broadcast(
103125
iEdep_indep -> broadcast(
104126
i_together -> begin
105127
nbcc = NBodyChargeCloud(
106-
CartesianPoint(to_internal_units.((pos[iEdep_indep][i_together]))),
107-
edep[iEdep_indep][i_together],
128+
CartesianPoint(T.(to_internal_units.(pos[iEdep_indep][i_together]))),
129+
T.(to_internal_units(edep[iEdep_indep][i_together])),
108130
number_of_carriers,
109-
number_of_shells = number_of_shells
131+
number_of_shells = number_of_shells,
132+
radius = T(to_internal_units(edep[iEdep_indep][i_together]))
110133
)
111134
move_charges_inside_semiconductor!([nbcc.locations], [nbcc.energies], det)
112135
nbcc
@@ -115,8 +138,8 @@ function _convertEnergyDepsToChargeDeps(pos::AbstractVector{<:AbstractVector}, e
115138
),
116139
eachindex(edep)
117140
)
118-
locations = [map(cc -> cc.locations, ccs) for ccs in charge_clouds]
119-
edeps = [map(cc -> cc.energies, ccs) for ccs in charge_clouds]
141+
locations = map.(cc -> cc.locations, charge_clouds)
142+
edeps = map.(cc -> cc.energies, charge_clouds)
120143
locations, edeps
121144
end
122145

@@ -127,9 +150,10 @@ function _simulate_charge_drifts( mcevents::TypedTables.Table, sim::Simulation{T
127150
self_repulsion::Bool,
128151
number_of_carriers::Int,
129152
number_of_shells::Int,
153+
max_interaction_distance::Union{<:Real, <:LengthQuantity},
130154
verbose::Bool ) where {T <: SSDFloat}
131155
@showprogress map(mcevents) do phyevt
132-
locations, edeps = _convertEnergyDepsToChargeDeps(phyevt.pos, phyevt.edep, sim.detector; number_of_carriers, number_of_shells)
156+
locations, edeps = _convertEnergyDepsToChargeDeps(phyevt.pos, phyevt.edep, sim.detector; number_of_carriers, number_of_shells, max_interaction_distance)
133157
drift_paths = map( i -> _drift_charges(sim.detector, sim.electric_field.grid, sim.point_types,
134158
VectorOfArrays(locations[i]), VectorOfArrays(edeps[i]),
135159
electric_field, T(Δt.val) * unit(Δt), max_nsteps = max_nsteps,

‎src/SolidStateDetectors.jl‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -48,6 +48,7 @@ import .ConstructiveSolidGeometry: sample, to_internal_units, from_internal_unit
4848
export CartesianPoint, CartesianVector, CylindricalPoint
4949

5050
import Clustering
51+
import DataStructures
5152
import Distributions
5253
import GPUArrays
5354
import IntervalSets

‎test/test_charge_drift_models.jl‎

Lines changed: 41 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -3,7 +3,8 @@
33
using Test
44

55
using SolidStateDetectors
6-
using SolidStateDetectors: getVe, getVh, Vl, get_path_to_example_config_files, AbstractChargeDriftModel
6+
using SolidStateDetectors: getVe, getVh, Vl, get_path_to_example_config_files, AbstractChargeDriftModel, group_points_by_distance, distance_squared
7+
using ArraysOfArrays
78
using InteractiveUtils
89
using StaticArrays
910
using LinearAlgebra
@@ -324,4 +325,43 @@ end
324325
@test cdm0.holes == cdmdict.holes
325326
@test cdm0.crystal_orientation == cdmdict.crystal_orientation
326327
end
328+
end
329+
330+
@timed_testset "Test grouping of charges" begin
331+
for N in 1:100
332+
@timed_testset "Test for length $(N)" begin
333+
334+
# Generate random data
335+
pts = [CartesianPoint{T}(rand(3)...) for _ in Base.OneTo(N)]
336+
E = rand(T, N)
337+
d = rand(T)
338+
339+
# Evaluate method
340+
ptsg, Eg = group_points_by_distance(pts, E, d)
341+
342+
# Test correctness
343+
s0 = Set(pts)
344+
345+
# result should be a VectorOfVectors
346+
@test ptsg isa VectorOfVectors{<:CartesianPoint{T}}
347+
348+
# all points should appear in the result
349+
@test Set(flatview(ptsg)) == Set(pts)
350+
@test length(flatview(ptsg)) == length(pts)
351+
@test Set(flatview(Eg)) == Set(E)
352+
@test length(flatview(Eg)) == length(E)
353+
354+
for (i,group) in enumerate(ptsg)
355+
@testset "Length $(N), subgroup $(i)" begin
356+
sg = Set(group)
357+
sc = setdiff(s0, group)
358+
for pt in group
359+
swo = setdiff(sg, Set([pt]))
360+
@test isempty(swo) || any(distance_squared.(Ref(pt), swo) .<= d^2)
361+
@test all(distance_squared.(Ref(pt), sc) .> d^2)
362+
end
363+
end
364+
end
365+
end
366+
end
327367
end

0 commit comments

Comments
 (0)