Introduce CustomIDArray to access Arrays with contact IDs

docs/src/man/plotting.md

@@ -190,7 +190,7 @@ The waveforms can be extended by calling [`add_baseline_and_extend_tail`](@ref)
 Again, the units of the axes can be set by calling a `plot` command with units before plotting the waveforms.
 ````@example tutorial
 plot(u"µs", u"fC")
-plot!(add_baseline_and_extend_tail.(evt.waveforms,0,400), 
+plot!(add_baseline_and_extend_tail(evt.waveforms,0,400), 
       linewidth = 4, linestyle = :dash, 
       label = "", unitformat = :slash)
 ````

examples/example_config_files/ivc_splitted_config/channels/01_PointContact.yaml

@@ -1,5 +1,5 @@
 material: HPGe
-id: 1
+id: 100
 potential: 0
 geometry:
   translate:

examples/example_config_files/ivc_splitted_config/channels/02_Mantle.yaml

@@ -1,5 +1,5 @@
 material: HPGe
-id: 2
+id: 200
 potential: 3500
 geometry:
   union:

src/CustomIDVector.jl

@@ -0,0 +1,63 @@
+struct CustomIDArray{T, N} <: AbstractArray{T, N}
+    data::Array{T, N}
+    idx::Vector{Int}
+end
+
+const CustomIDVector{T} = CustomIDArray{T, 1}
+CustomIDVector(data::Vector{T}, idx::Vector{Int}) where {T} = CustomIDVector{T}(data, idx)
+
+const CustomIDMatrix{T} = CustomIDArray{T, 2}
+CustomIDMatrix(data::Matrix{T}, idx::Vector{Int}) where {T} = CustomIDMatrix{T}(data, idx)
+
+@inline size(cia::CustomIDArray) = size(cia.data)
+@inline length(cia::CustomIDArray) = length(cia.data)
+@inline iterate(cia::CustomIDArray, args...) = iterate(cia.data, args...)
+@inline axes(cia::CustomIDVector) = (cia.idx,)
+
+@inline function getindex(cia::CustomIDVector, i::Int)
+    !(i in cia.idx) && throw(ArgumentError("No entry for contact with ID "*string(i)))
+    idx::Int = Int(findfirst(ix -> ix == i, cia.idx))
+    getindex(cia.data, idx)
+end
+
+@inline function getindex(cia::CustomIDArray{T,N}, I::Vararg{Int,N}) where {T, N}
+    !all(in.(I, Ref(cia.idx))) && throw(ArgumentError("No entry for contact with ID "*string(I)))
+    idx = broadcast(i -> findfirst(idx -> idx == i, cia.idx), I)
+    getindex(cia.data, idx...)
+end
+
+@inline function push!(cia::CustomIDVector{T}, data::T, i::Int) where {T}
+    push!(cia.data, data)
+    push!(cia.idx, i)
+end
+
+@inline function setindex!(cia::CustomIDVector{T}, data::T, i::Int) where {T}
+    if i in cia.idx
+        idx::Int = Int(findfirst(idx -> idx == i, cia.idx))
+        setindex!(cia.data, data, idx)
+    else
+        push!(cia, data, i)
+    end
+end
+
+@inline function setindex!(cia::CustomIDArray{T, N}, data::T, I::Vararg{Int,N}) where {T, N}
+    !all(in.(I, Ref(cia.idx))) && throw(ArgumentError("No entry for contact with ID "*string(I)))
+    idx = broadcast(i -> findfirst(idx -> idx == i, cia.idx), I)
+    setindex!(cia.data, data, idx...)
+end
+
+@inline function setindex!(cia::CustomIDArray{T, N}, data, I::Vararg{Int, N}) where {T, N}
+    setindex!(cia, convert(T, data), I...)
+end
+
+@inline function Base.isequal(cia1::CIA1, cia2::CIA2) where {CIA1 <: CustomIDArray, CIA2 <: CustomIDArray}
+    CIA1 == CIA2 && cia1.idx == cia2.idx && cia1.data == cia2.data
+end
+
+@inline function Base.:(==)(cia1::CIA1, cia2::CIA2) where {CIA1 <: CustomIDArray, CIA2 <: CustomIDArray}
+    CIA1 == CIA2 && cia1.idx == cia2.idx && cia1.data == cia2.data
+end
+
+@inline print(io::IO, cia::CustomIDArray) = print(io, cia.data)
+@inline show(io::IO, cia::CustomIDArray) = print(io, cia)
+@inline show(io::IO, ::MIME"text/plain", cia::CustomIDArray) = show(io, cia)

src/Event/Event.jl

@@ -15,8 +15,8 @@ mutable struct Event{T <: SSDFloat}
     locations::VectorOfArrays{CartesianPoint{T}}
     energies::VectorOfArrays{T}
     drift_paths::Union{Vector{EHDriftPath{T}}, Missing}
-    waveforms::Union{Vector{<:Any}, Missing}
-    Event{T}() where {T <: SSDFloat} = new{T}(VectorOfArrays(Vector{CartesianPoint{T}}[]), VectorOfArrays(Vector{T}[]), missing, missing)
+    waveforms::CustomIDVector{RadiationDetectorSignals.RDWaveform}
+    Event{T}() where {T <: SSDFloat} = new{T}(VectorOfArrays(Vector{CartesianPoint{T}}[]), VectorOfArrays(Vector{T}[]), missing, CustomIDVector(RadiationDetectorSignals.RDWaveform[], Int[]))
 end
 
 function Event(location::AbstractCoordinatePoint{T}, energy::RealQuantity = one(T))::Event{T} where {T <: SSDFloat}
@@ -35,7 +35,7 @@ end
 
 function Event(locations::Vector{<:Vector{<:AbstractCoordinatePoint{T}}}, energies::Vector{<:Vector{<:RealQuantity}} = [[one(T) for i in j] for j in locations])::Event{T} where {T <: SSDFloat}
     evt = Event{T}()
-    evt.locations = VectorOfArrays(broadcast(pts -> CartesianPoint{T}.(pts), locations))
+    evt.locations = VectorOfArrays(broadcast(pts -> CartesianPoint.(pts), locations))
     evt.energies = VectorOfArrays(Vector{T}.(to_internal_units.(energies)))
     return evt
 end
@@ -155,10 +155,7 @@ function drift_charges!(evt::Event{T}, sim::Simulation{T}; max_nsteps::Int = 100
 end
 function get_signal!(evt::Event{T}, sim::Simulation{T}, contact_id::Int; Δt::RealQuantity = 5u"ns")::Nothing where {T <: SSDFloat}
     @assert !ismissing(evt.drift_paths) "No drift path for this event. Use `drift_charges(evt::Event, sim::Simulation)` first."
-    @assert !ismissing(sim.weighting_potentials[contact_id]) "No weighting potential for contact $(contact_id). Use `calculate_weighting_potential!(sim::Simulation, contact_id::Int)` first."
-    if ismissing(evt.waveforms)
-        evt.waveforms = Union{Missing, RadiationDetectorSignals.RDWaveform}[missing for i in eachindex(sim.detector.contacts)]
-    end
+    @assert contact_id in sim.weighting_potentials.idx "No weighting potential for contact $(contact_id). Use `calculate_weighting_potential!(sim::Simulation, contact_id::Int)` first."
     evt.waveforms[contact_id] = get_signal(sim, evt.drift_paths, flatview(evt.energies), contact_id, Δt = Δt)
     nothing
 end
@@ -189,12 +186,8 @@ get_signals!(evt, sim, Δt = 1u"ns") # if evt.drift_paths were calculated in tim
 """
 function get_signals!(evt::Event{T}, sim::Simulation{T}; Δt::RealQuantity = 5u"ns")::Nothing where {T <: SSDFloat}
     @assert !ismissing(evt.drift_paths) "No drift path for this event. Use `drift_charges!(evt::Event, sim::Simulation)` first."
-    if ismissing(evt.waveforms)
-        evt.waveforms = Union{Missing, RadiationDetectorSignals.RDWaveform}[missing for i in eachindex(sim.detector.contacts)]
-    end
     for contact in sim.detector.contacts
-        if any(ismissing, sim.weighting_potentials) "No weighting potential(s) for some contact(s).." end
-        if !ismissing(sim.weighting_potentials[contact.id])
+        if contact.id in sim.weighting_potentials.idx
             evt.waveforms[contact.id] = get_signal(sim, evt.drift_paths, flatview(evt.energies), contact.id, Δt = Δt)
         end
     end
@@ -271,6 +264,10 @@ function add_baseline_and_extend_tail(wv::RadiationDetectorSignals.RDWaveform{T,
     return RDWaveform( new_times, new_signal )
 end
 
+function add_baseline_and_extend_tail(wvs::CustomIDVector{RadiationDetectorSignals.RDWaveform}, args...)
+    CustomIDVector(add_baseline_and_extend_tail.(wvs.data, args...), wvs.idx)
+end
+
 
 """
     get_electron_and_hole_contribution(evt::Event{T}, sim::Simulation{T}, contact_id::Int)

src/PlotRecipes/Waveform.jl

@@ -1,6 +1,8 @@
-# This should be in RadiationDetectorSignals.jl
-@recipe function f(wvs::Vector{Union{Missing, RadiationDetectorSignals.RDWaveform}})
+@recipe function f(wvs::CustomIDVector{<:RadiationDetectorSignals.RDWaveform})
     @series begin
-        RadiationDetectorSignals.RDWaveform[wv for wv in skipmissing(wvs)]
+        label --> wvs.idx'
+        linecolor --> wvs.idx'
+        unitformat --> :slash
+        [wv for wv in wvs.data]
     end
 end

src/Simulation/Capacitance.jl

@@ -45,7 +45,9 @@ function _calculate_mutual_capacitance(
 
     phi_2D && (c_ij *= 2)
     consider_multiplicity && (c_ij *= multiplicity(grid))
-    return uconvert(u"pF", c_ij*u"m" * ϵ0*u"F/m" )
+
+    # ϵ0 is stored as Float64, so convert to T here
+    return T(uconvert(u"pF", c_ij*u"m" * ϵ0*u"F/m"))
 end
 
 function _calculate_mutual_capacitance(grid::Grid{T, 3, CS}, grid_mps, int_ϵ_r, int_p1, int_p2) where {T, CS}
@@ -123,11 +125,13 @@ to numerical precision in the integration due to the different grids of the two
 function calculate_capacitance_matrix(sim::Simulation{T}; consider_multiplicity::Bool = true) where {T}
     @assert !ismissing(sim.ϵ_r) "The electric potential needs to be calculated first."
     @assert !ismissing(sim.weighting_potentials) "The weighting_potentials needs to be calculated first."
-    n = length(sim.weighting_potentials)
-    C = zeros(typeof(one(T) * u"pF"), (n, n))
-    for i in 1:n
-        for j in 1:n
-            C[j, i] = if !ismissing(sim.weighting_potentials[i]) && !ismissing(sim.weighting_potentials[j]) 
+    n::Int = length(sim.detector.contacts)
+    contact_ids::Vector{Int} = map(c -> c.id, sim.detector.contacts)
+    TT = all(in.(contact_ids, Ref(sim.weighting_potentials.idx))) ? typeof(one(T) * u"pF") : Union{Missing,typeof(one(T) * u"pF")}
+    C = CustomIDMatrix{TT}(zeros(typeof(one(T) * u"pF"), (n, n)), contact_ids)
+    for i in contact_ids
+        for j in contact_ids
+            C[j, i] = if (i in sim.weighting_potentials.idx) && (j in sim.weighting_potentials.idx)
                 calculate_mutual_capacitance(sim, (i, j); consider_multiplicity)
             else
                 missing

src/Simulation/Simulation.jl

@@ -36,7 +36,7 @@ mutable struct Simulation{T <: SSDFloat, CS <: AbstractCoordinateSystem} <: Abst
     ϵ_r::Union{DielectricDistribution{T}, Missing}
     point_types::Union{PointTypes{T}, Missing}
     electric_potential::Union{ElectricPotential{T}, Missing}
-    weighting_potentials::Vector{Any}
+    weighting_potentials::CustomIDVector{<:WeightingPotential}
     electric_field::Union{ElectricField{T}, Missing}
 end
 
@@ -53,7 +53,7 @@ function Simulation{T,CS}() where {T <: SSDFloat, CS <: AbstractCoordinateSystem
         missing,
         missing,
         missing,
-        [missing],
+        CustomIDVector(WeightingPotential{T,3,CS}[], Int[]),
         missing
     )
 end
@@ -72,7 +72,9 @@ function NamedTuple(sim::Simulation{T}) where {T <: SSDFloat}
         ϵ_r = NamedTuple(sim.ϵ_r),
         point_types = NamedTuple(sim.point_types),
         electric_field = NamedTuple(sim.electric_field),
-        weighting_potentials = NamedTuple{Tuple(Symbol.(wpots_strings))}(NamedTuple.(sim.weighting_potentials))
+        weighting_potentials = let wp = sim.weighting_potentials
+            NamedTuple{Tuple(Symbol.("WeightingPotential_".*string.(wp.idx)))}(NamedTuple.(wp.data))
+        end
     )
     return nt
 end
@@ -101,12 +103,16 @@ function Simulation(nt::NamedTuple)
         ImpurityScale(nt.imp_scale)
     end
     sim.electric_field = haskey(nt, :electric_field) && nt.electric_field !== missing_tuple ? ElectricField(nt.electric_field) : missing
-    sim.weighting_potentials = if haskey(nt, :weighting_potentials) 
-        [let wp = Symbol("WeightingPotential_$(contact.id)")
-            haskey(nt.weighting_potentials, wp) && getfield(nt.weighting_potentials, wp) !== missing_tuple ? WeightingPotential(getfield(nt.weighting_potentials, wp)) : missing 
-        end for contact in sim.detector.contacts]
-    else
-        [missing for contact in sim.detector.contacts]
+    grid = Grid(sim, for_weighting_potential = true)
+    sim.weighting_potentials = CustomIDVector(WeightingPotential{T,3,get_coordinate_system(grid),typeof(grid.axes)}[], Int[])
+    if haskey(nt, :weighting_potentials) 
+        for contact in sim.detector.contacts
+            let wp = Symbol("WeightingPotential_$(contact.id)")
+                if haskey(nt.weighting_potentials, wp) && getfield(nt.weighting_potentials, wp) !== missing_tuple 
+                    sim.weighting_potentials[contact.id] = WeightingPotential(getfield(nt.weighting_potentials, wp))
+                end
+            end
+        end
     end
     return sim
 end
@@ -131,7 +137,7 @@ function println(io::IO, sim::Simulation{T}) where {T <: SSDFloat}
     println("  Weighting potentials: ")
     for contact in sim.detector.contacts
         print("    Contact $(contact.id): ")
-        println(!ismissing(sim.weighting_potentials[contact.id]) ? size(sim.weighting_potentials[contact.id]) : missing)
+        println(contact.id in sim.weighting_potentials.idx ? size(sim.weighting_potentials[contact.id]) : missing)
     end
 end
 
@@ -179,7 +185,8 @@ function Simulation{T}(dict::Dict)::Simulation{T} where {T <: SSDFloat}
             World(CS, world_limits)
         end
     end
-    sim.weighting_potentials = Missing[ missing for i in 1:length(sim.detector.contacts)]
+    grid = Grid(sim, for_weighting_potential = true)
+    sim.weighting_potentials = CustomIDVector(WeightingPotential{T,3,CS,typeof(grid.axes)}[], Int[])
     return sim
 end
 

src/SolidStateDetectors.jl

@@ -53,7 +53,7 @@ import IntervalSets
 import Tables
 import TypedTables
 
-import Base: size, sizeof, length, getindex, setindex!, axes, getproperty, broadcast,
+import Base: size, sizeof, length, getindex, setindex!, push!, axes, getproperty, broadcast, isequal,
              range, ndims, eachindex, enumerate, iterate, IndexStyle, eltype, in, convert,
              show, print, println, display, +, -, &
 
@@ -78,6 +78,7 @@ const SSDFloat = Union{Float16, Float32, Float64}
 
 
 include("Units.jl")
+include("CustomIDVector.jl")
 include("Axes/DiscreteAxis.jl")
 include("World/World.jl")
 include("Grids/Grids.jl")

test/comparison_to_analytic_solutions.jl

@@ -20,7 +20,7 @@ struct DummyImpurityDensity{T} <: SolidStateDetectors.AbstractImpurityDensity{T}
     W_true = uconvert(u"J", (ϵr * ϵ0 * E_true^2 / 2) * A * Δd)
     C_true = uconvert(u"pF", 2 * W_true / (BV_true^2))
     C_Analytical = [C_true -C_true; -C_true C_true]
-    C_ssd = calculate_capacitance_matrix(sim)
+    C_ssd = calculate_capacitance_matrix(sim).data
     @testset "Capacity" begin
         @test all(isapprox.(C_ssd, C_Analytical, rtol = 0.001))
     end