Update simulate_waveforms to allow for max_interaction_distance

Author: fhagemann

ext/SolidStateDetectorsLegendHDF5IOExt.jl

@@ -5,7 +5,7 @@ module SolidStateDetectorsLegendHDF5IOExt
 import ..LegendHDF5IO
 
 using SolidStateDetectors
-using SolidStateDetectors: RealQuantity, SSDFloat, to_internal_units, chunked_ranges
+using SolidStateDetectors: RealQuantity, SSDFloat, to_internal_units, chunked_ranges, LengthQuantity
 using TypedTables, Unitful
 using Format
 
@@ -34,6 +34,7 @@ function SolidStateDetectors.simulate_waveforms( mcevents::TypedTables.Table, si
                              self_repulsion::Bool = false,
                              number_of_carriers::Int = 1,
                              number_of_shells::Int = 1,
+                             max_interaction_distance::Union{<:Real, <:LengthQuantity} = NaN,
                              verbose = false) where {T <: SSDFloat}
     n_total_physics_events = length(mcevents)
     Δtime = T(to_internal_units(Δt)) 
@@ -48,7 +49,7 @@ function SolidStateDetectors.simulate_waveforms( mcevents::TypedTables.Table, si
     for evtrange in evt_ranges
         ofn = joinpath(output_dir, "$(output_base_name)_evts_$(nfmt(first(evtrange)))-$(nfmt(last(evtrange))).h5")
         @info "Now simulating $(evtrange) and storing it in\n\t \"$ofn\""
-        mcevents_sub = simulate_waveforms(mcevents[evtrange], sim; Δt, max_nsteps, diffusion, self_repulsion, number_of_carriers, number_of_shells, verbose)
+        mcevents_sub = simulate_waveforms(mcevents[evtrange], sim; Δt, max_nsteps, diffusion, self_repulsion, number_of_carriers, number_of_shells, max_interaction_distance, verbose)
 
         LegendHDF5IO.lh5open(ofn, "w") do h5f
             LegendHDF5IO.writedata(h5f.data_store, "generated_waveforms", mcevents_sub)

src/MCEventsProcessing/MCEventsProcessing.jl

@@ -27,6 +27,7 @@ If [LegendHDF5IO.jl](https://github.com/legend-exp/LegendHDF5IO.jl) is loaded, t
 * `self_repulsion::Bool = false`: Activate or deactive self-repulsion of charge carriers of the same type.
 * `number_of_carriers::Int = 1`: Number of charge carriers to be used in the N-Body simulation of an energy deposition. 
 * `number_of_shells::Int = 1`: Number of shells around the `center` point of the energy deposition.
+* `max_interaction_distance = NaN`: Maximum distance for which charge clouds will interact with each other (if `NaN`, then all charge clouds drift independently).
 * `verbose = false`: Activate or deactivate additional info output.
 * `chunk_n_physics_events::Int = 1000` (`LegendHDF5IO` only): Number of events that should be saved in a single HDF5 output file.
 
@@ -53,6 +54,7 @@ function simulate_waveforms( mcevents::TypedTables.Table, sim::Simulation{T};
                              self_repulsion::Bool = false,
                              number_of_carriers::Int = 1,
                              number_of_shells::Int = 1,
+                             max_interaction_distance::Union{<:Real, <:LengthQuantity} = NaN,
                              verbose::Bool = false ) where {T <: SSDFloat}
     n_total_physics_events = length(mcevents)
     Δtime = T(to_internal_units(Δt)) 
@@ -71,7 +73,7 @@ function simulate_waveforms( mcevents::TypedTables.Table, sim::Simulation{T};
 
     # First simulate drift paths
     drift_paths_and_edeps = _simulate_charge_drifts(mcevents, sim, Δt, max_nsteps, electric_field, 
-        diffusion, self_repulsion, number_of_carriers, number_of_shells, verbose)
+        diffusion, self_repulsion, number_of_carriers, number_of_shells, max_interaction_distance, verbose)
     drift_paths = map(x -> vcat([vcat(ed...) for ed in x[1]]...), drift_paths_and_edeps)
     edeps = map(x -> vcat([vcat(ed...) for ed in x[2]]...), drift_paths_and_edeps)
     # now iterate over contacts and generate the waveform for each contact
@@ -94,19 +96,40 @@ function simulate_waveforms( mcevents::TypedTables.Table, sim::Simulation{T};
     return vcat(mcevents_chns...)  
 end
 
-_convertEnergyDepsToChargeDeps(pos::AbstractVector{<:SVector{3}}, edep::AbstractVector{<:Quantity}, det::SolidStateDetector; kwargs...) = 
-    _convertEnergyDepsToChargeDeps([[p] for p in pos], [[e] for e in edep], det; kwargs...)
+function _convertEnergyDepsToChargeDeps(
+        pos::AbstractVector{<:SVector{3}}, edep::AbstractVector{<:Quantity}, det::SolidStateDetector{T}; 
+        particle_type::Type{PT} = Gamma, 
+        radius::Vector{<:Union{<:Real, <:LengthQuantity}} = radius_guess.(T.(to_internal_units.(edep)), particle_type),
+        max_interaction_distance::Union{<:Real, <:LengthQuantity} = NaN,
+        kwargs...
+    ) where {T <: SSDFloat, PT <: ParticleType} 
+
+    @assert eachindex(pos) == eachindex(edep) == eachindex(radius)
+
+    d::T = T(to_internal_units(max_interaction_distance))
+    @assert isnan(d) || d >= 0 "Max. interaction distance must be positive or NaN (no grouping), but $(max_interaction_distance) was given."
+
+    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)
+    _convertEnergyDepsToChargeDeps(loc, ene, det; radius = rad, kwargs...)
+end
+
+function _convertEnergyDepsToChargeDeps(
+        pos::AbstractVector{<:AbstractVector}, edep::AbstractVector{<:AbstractVector}, det::SolidStateDetector{T}; 
+        particle_type::Type{PT} = Gamma,
+        radius::AbstractVector{<:AbstractVector{<:Union{<:Real, <:LengthQuantity}}} = map(e -> radius_guess.(to_internal_units.(e), particle_type), edep),
+        number_of_carriers::Int = 1, number_of_shells::Int = 1, 
+        max_interaction_distance::Union{<:Real, <:LengthQuantity} = NaN # is ignored here
+    ) where {T <: SSDFloat, PT <: ParticleType}
 
-function _convertEnergyDepsToChargeDeps(pos::AbstractVector{<:AbstractVector}, edep::AbstractVector{<:AbstractVector}, det::SolidStateDetector; 
-        number_of_carriers::Int = 1, number_of_shells::Int = 1)
     charge_clouds = broadcast(
         iEdep_indep -> broadcast(
             i_together -> begin 
                 nbcc = NBodyChargeCloud(
-                    CartesianPoint(to_internal_units.((pos[iEdep_indep][i_together]))), 
-                    edep[iEdep_indep][i_together],
+                    CartesianPoint(T.(to_internal_units.(pos[iEdep_indep][i_together]))), 
+                    T.(to_internal_units(edep[iEdep_indep][i_together])),
                     number_of_carriers,
-                    number_of_shells = number_of_shells
+                    number_of_shells = number_of_shells,
+                    radius = T(to_internal_units(edep[iEdep_indep][i_together]))
                 )
                 move_charges_inside_semiconductor!([nbcc.locations], [nbcc.energies], det)
                 nbcc
@@ -115,8 +138,8 @@ function _convertEnergyDepsToChargeDeps(pos::AbstractVector{<:AbstractVector}, e
         ), 
         eachindex(edep)
     )
-    locations = [map(cc -> cc.locations, ccs) for ccs in charge_clouds]
-    edeps = [map(cc -> cc.energies, ccs) for ccs in charge_clouds]
+    locations = map.(cc -> cc.locations, charge_clouds)
+    edeps = map.(cc -> cc.energies, charge_clouds)
     locations, edeps
 end
 
@@ -127,9 +150,10 @@ function _simulate_charge_drifts( mcevents::TypedTables.Table, sim::Simulation{T
                                   self_repulsion::Bool,
                                   number_of_carriers::Int, 
                                   number_of_shells::Int,
+                                  max_interaction_distance::Union{<:Real, <:LengthQuantity},
                                   verbose::Bool ) where {T <: SSDFloat}
     @showprogress map(mcevents) do phyevt
-        locations, edeps = _convertEnergyDepsToChargeDeps(phyevt.pos, phyevt.edep, sim.detector; number_of_carriers, number_of_shells)
+        locations, edeps = _convertEnergyDepsToChargeDeps(phyevt.pos, phyevt.edep, sim.detector; number_of_carriers, number_of_shells, max_interaction_distance)
         drift_paths = map( i -> _drift_charges(sim.detector, sim.electric_field.grid, sim.point_types, 
                 VectorOfArrays(locations[i]), VectorOfArrays(edeps[i]),
                 electric_field, T(Δt.val) * unit(Δt), max_nsteps = max_nsteps,