Create new datatype when saving CartesianPoint to HDF5 using LegendHDF5IO

docs/src/man/charge_drift.md

@@ -416,7 +416,7 @@ using SolidStateDetectors #hide
 using Unitful #hide
 using Plots #hide
 T = Float64 #hide
-center = CartesianPoint{T}([0,0,0])
+center = CartesianPoint{T}(0,0,0)
 energy = 1460u"keV"
 nbcc = NBodyChargeCloud(center, energy)
 plot(nbcc)
@@ -428,7 +428,7 @@ plot(nbcc, color = :red, size = (500,500), xlims = (-0.0012, 0.0012), ylims = (-
 For an [`NBodyChargeCloud`](@ref) consisting of more than around 50 charges, the shells should consist of more than 20 point charges and the approach with using [Platonic Solids](@ref) for the shell structure might not be favored anymore. For this, a [second algorithm](https://www.cmu.edu/biolphys/deserno/pdf/sphere_equi.pdf) was implemented that generates point charges equally distributed on the surface of a regular sphere. The approximate number of charges needs to be passed to the constructor of [`NBodyChargeCloud`](@ref) to use this method.
 
 ````@example NBodyChargeCloud
-center = CartesianPoint{T}([0,0,0])
+center = CartesianPoint{T}(0,0,0)
 energy = 1460u"keV"
 nbcc = NBodyChargeCloud(center, energy, 100)
 plot(nbcc)
@@ -441,7 +441,7 @@ plot(nbcc, color = :red, size = (500,500), xlims = (-0.0012, 0.0012), ylims = (-
 Diffusion describes the random thermal motion of charge carriers. In SolidStateDetectors.jl, diffusion is simulated using a random walk algorithm. Diffusion is simulated using fixed-step vectors where the magnitude of the step vectors depends on the diffusion constant $D$ and the time steps $\Delta t$.
 
 ```julia
-center = CartesianPoint{T}([0,0,0])
+center = CartesianPoint{T}(0,0,0)
 energy = 1460u"keV"
 nbcc = NBodyChargeCloud(center, energy, 100)
 evt = Event(nbcc)
@@ -478,7 +478,7 @@ After the creation electron-hole pairs, both the electron and the hole clouds re
 SolidStateDetectors.jl does not account for attraction of electron and holes but only for repulsion of charge carriers of the same type. The determination of the electric field vector is calculated pair-wise for each pair of charge carriers.
 
 ```julia
-center = CartesianPoint{T}([0,0,0])
+center = CartesianPoint{T}(0,0,0)
 energy = 1460u"keV"
 nbcc = NBodyChargeCloud(center, energy, 100)
 evt = Event(nbcc)
@@ -491,7 +491,7 @@ simulate!(evt, sim, self_repulsion = true)
 
 [Diffusion](@ref) and [Self-Repulsion](@ref) can be simulated both at once to get the most realistic picture:
 ```julia
-center = CartesianPoint{T}([0,0,0])
+center = CartesianPoint{T}(0,0,0)
 energy = 1460u"keV"
 nbcc = NBodyChargeCloud(center, energy, 100)
 evt = Event(nbcc)

ext/SolidStateDetectorsLegendHDF5IOExt.jl

@@ -6,6 +6,8 @@ import ..LegendHDF5IO
 
 using ..SolidStateDetectors
 using ..SolidStateDetectors: RealQuantity, SSDFloat, to_internal_units, chunked_ranges, LengthQuantity
+using ..SolidStateDetectors.ConstructiveSolidGeometry: AbstractCoordinatePoint
+
 using ArraysOfArrays
 import Tables
 using TypedTables, Unitful, UnitfulAtomic
@@ -65,4 +67,97 @@ function SolidStateDetectors.simulate_waveforms( mcevents::AbstractVector{<:Name
     end    
 end
 
+
+# Add support for reading and writing CartesianPoint and CylindricalPoint using LegendHDF5IO
+function __init__()
+    LegendHDF5IO._datatype_dict["cartesianpoint"] = CartesianPoint
+    LegendHDF5IO._datatype_dict["cylindricalpoint"] = CylindricalPoint
+end
+
+LegendHDF5IO.datatype_to_string(::Type{<:CartesianPoint}) = "cartesianpoint"
+LegendHDF5IO.datatype_to_string(::Type{<:CylindricalPoint}) = "cylindricalpoint"
+
+
+
+# support writing points to HDF5 file using LH5Array
+function LegendHDF5IO.create_entry(parent::LegendHDF5IO.LHDataStore, name::AbstractString, 
+    pt::T; kwargs...) where {T <: AbstractCoordinatePoint}
+    LegendHDF5IO.create_entry(parent, name, getindex.(Ref(pt), 1:3); kwargs...)
+    LegendHDF5IO.setdatatype!(parent.data_store[name], T)
+    nothing
+end
+
+function LegendHDF5IO.LH5Array(ds::LegendHDF5IO.HDF5.Dataset, T::Type{<:AbstractCoordinatePoint})
+    T(read(ds) .* LegendHDF5IO.getunits(ds)...)
+end
+
+# support writing arrays of points to HDF5 file using LH5Array
+function LegendHDF5IO.create_entry(parent::LegendHDF5IO.LHDataStore, name::AbstractString, 
+    pts::T; kwargs...) where {T <: AbstractVector{<:AbstractCoordinatePoint}}
+    LegendHDF5IO.create_entry(parent, name, VectorOfVectors([getindex.(Ref(pt), 1:3) for pt in pts]); kwargs...)
+    LegendHDF5IO.setdatatype!(parent.data_store[name], T)
+    nothing
+end
+
+function LegendHDF5IO.LH5Array(ds::LegendHDF5IO.HDF5.Group, T::Type{<:AbstractVector{<:CartesianPoint}})
+    data = LegendHDF5IO.LH5Array(ds, VectorOfVectors)
+    Vector([CartesianPoint(d...) for d in data])
+end
+
+function LegendHDF5IO.LH5Array(ds::LegendHDF5IO.HDF5.Group, T::Type{<:AbstractVector{<:CylindricalPoint}})
+    data = LegendHDF5IO.LH5Array(ds, VectorOfVectors)
+    Vector([CylindricalPoint(d...) for d in data])
+end
+
+# support writing points to HDF5 file using readdata/writedata
+function LegendHDF5IO.writedata(
+    output::LegendHDF5IO.HDF5.H5DataStore, name::AbstractString,
+    pt::AbstractCoordinatePoint,
+    fulldatatype::DataType = typeof(pt)
+)
+    data = getindex.(Ref(pt), 1:3)
+    LegendHDF5IO.writedata(output, name, data, fulldatatype)
+end
+
+function LegendHDF5IO.readdata(
+    input::LegendHDF5IO.HDF5.H5DataStore, name::AbstractString,
+    fulldatatype::Type{T}
+) where {T <: AbstractCoordinatePoint}
+    data = LegendHDF5IO.readdata(input, name, Vector)
+    T(data...)
+end
+
+
+# support writing arrays of points to HDF5 file using readdata/writedata
+function LegendHDF5IO.writedata(
+    output::LegendHDF5IO.HDF5.H5DataStore, name::AbstractString,
+    pts::AbstractVector{<:AbstractCoordinatePoint},
+    fulldatatype::DataType = typeof(pts)
+)
+    data = VectorOfVectors([getindex.(Ref(pt), 1:3) for pt in pts])
+    LegendHDF5IO.writedata(output, name, data, fulldatatype)
+    LegendHDF5IO.setdatatype!(output[name], fulldatatype)
+    nothing
+end
+
+function LegendHDF5IO.readdata(
+    input::LegendHDF5IO.HDF5.H5DataStore, name::AbstractString,
+    fulldatatype::Type{T}
+) where {T <:AbstractVector{<:CartesianPoint}}
+    data = LegendHDF5IO.readdata(input, name, VectorOfVectors)
+    output = [CartesianPoint(d...) for d in data]
+    @assert output isa T
+    output
+end
+
+function LegendHDF5IO.readdata(
+    input::LegendHDF5IO.HDF5.H5DataStore, name::AbstractString,
+    fulldatatype::Type{T}
+) where {T <:AbstractVector{<:CylindricalPoint}}
+    data = LegendHDF5IO.readdata(input, name, VectorOfVectors)
+    output = [CylindricalPoint(d...) for d in data]
+    @assert output isa T
+    output
+end
+
 end # module LegendHDF5IO

src/ChargeCloudModels/ChargeCloudModels.jl

@@ -52,6 +52,12 @@ NBodyChargeCloud(center, energy, number_of_shells = 3, shell_structure = SolidSt
 !!! note
     Using values with units for `energy` requires the package [Unitful.jl](https://github.com/JuliaPhysics/Unitful.jl).
 """
+function NBodyChargeCloud(center::Union{<:CartesianPoint{<:Unitful.Quantity},<:CylindricalPoint{<:Union{SSDFloat, Unitful.Quantity}}},
+    energy::RealQuantity; kwargs...)
+    internal_center = to_internal_units(center)
+    return NBodyChargeCloud(internal_center, energy; kwargs...)
+end
+
 function NBodyChargeCloud(center::CartesianPoint{T}, energy::RealQuantity, particle_type::Type{PT} = Gamma;
         radius::T = radius_guess(T(to_internal_units(energy)), particle_type), number_of_shells::Int = 2, shell_structure = Dodecahedron
     )::NBodyChargeCloud{T, number_of_shells, typeof(shell_structure{T})} where {T, PT <: ParticleType}
@@ -121,6 +127,12 @@ NBodyChargeCloud(center, energy, 200, number_of_shells = 3)
 !!! note
     Using values with units for `energy` requires the package [Unitful.jl](https://github.com/JuliaPhysics/Unitful.jl).
 """
+function NBodyChargeCloud(center::Union{<:CartesianPoint{<:Unitful.Quantity},<:CylindricalPoint{<:Union{SSDFloat, Unitful.Quantity}}},
+    energy::RealQuantity, N::Integer; kwargs...)
+    internal_center = to_internal_units(center)
+    return NBodyChargeCloud(internal_center, energy, N; kwargs...)
+end
+
 function NBodyChargeCloud(center::CartesianPoint{T}, energy::RealQuantity, N::Integer, particle_type::Type{PT} = Gamma;
         radius::T = radius_guess(T(to_internal_units(energy)), particle_type), number_of_shells::Int = 2,
     )::NBodyChargeCloud{T, number_of_shells, NTuple{N > 1 ? number_of_shells : 0, Int}} where {T, PT <: ParticleType}
@@ -142,4 +154,4 @@ function NBodyChargeCloud(center::CartesianPoint{T}, energy::RealQuantity, N::In
     end
 
     return NBodyChargeCloud{T, number_of_shells, NTuple{length(n), Int}}( locations, energies./sum(energies) * T(to_internal_units(energy)), Tuple(n))
-end
+end

src/ChargeClustering/ChargeClustering.jl

@@ -16,8 +16,7 @@ function cluster_detector_hits(
     r_edep = similar(edep, 0)
     r_pos = similar(pos, 0)
 
-    posunit = unit(PT)
-    ustripped_cradius = ustrip(posunit, cluster_radius)
+    ustripped_cradius = to_internal_units(cluster_radius)
 
     for d_hits_nt in grouped
         d_hits = TypedTables.Table(d_hits_nt)

src/ConstructiveSolidGeometry/PointsAndVectors/Points.jl

@@ -13,12 +13,6 @@ function Base.isapprox(a::AbstractCartesianPoint, b::AbstractCartesianPoint; kwa
            isapprox(get_z(a), get_z(b); kwargs...)
 end
 
-function to_internal_units(x::AbstractCartesianPoint)
-    CartesianPoint(to_internal_units(get_x(x)), to_internal_units(get_y(x)), to_internal_units(get_z(x)))
-end
-
-
-
 # Unitful uses the `*` and `/` operators to combine values with units. But mathematically that's really not an algebraic
 # product (not defined for affine points), but a cartesian product, so supporting this should be fine:
 Base.:(*)(pt::AbstractCartesianPoint{<:Real}, u::Unitful.Units{<:Any,Unitful.𝐋}) = CartesianPoint(get_x(pt) * u, get_y(pt) * u, get_z(pt) * u)
@@ -72,6 +66,22 @@ function CartesianPoint(x::TX, y::TY, z::TZ) where {TX<:Real,TY<:Real,TZ<:Real}
     CartesianPoint{T}(T(x),T(y),T(z))
 end
 
+#Unit support
+function CartesianPoint(x::Unitful.Quantity, y::Unitful.Quantity, z::Unitful.Quantity)
+    for (i, pt) in enumerate((x, y, z))
+        !(pt isa Unitful.Length) && throw(ArgumentError(
+            "Coordinate $i must be a length, got $(pt) with unit $(Unitful.unit(pt))"
+        ))
+    end
+    CartesianPoint(x::Unitful.Length, y::Unitful.Length, z::Unitful.Length)
+end
+
+function CartesianPoint(x::Unitful.Length, y::Unitful.Length, z::Unitful.Length)
+    vals = to_internal_units.( (x, y, z) )
+    T = all(v -> v isa Float32, vals) ? Float32 : Float64
+    CartesianPoint{T}(T.(vals)...)
+end
+
 CartesianPoint(; x = 0, y = 0, z = 0) = CartesianPoint(x, y, z)
 
 CartesianPoint{T}(;x = 0, y = 0, z = 0) where {T} = CartesianPoint{T}(T(x),T(y),T(z))
@@ -204,7 +214,8 @@ end
 @inline Base.zero(::CartesianZero{T}) where {T} = CartesianZero{T}()
 @inline Base.iszero(::CartesianZero) = true
 
-Base.:(*)(::CartesianZero{T}, u::Unitful.Units{<:Any,Unitful.𝐋}) where {T<:Real} = CartesianZero{Quantity{T, Unitful.𝐋, typeof(u)}}()
+#Base.:(*)(::CartesianZero{T}, u::Unitful.Units{<:Any,Unitful.𝐋}) where {T<:Real} = CartesianZero{Quantity{T, Unitful.𝐋, typeof(u)}}()
+Base.:*(z::CartesianZero, u::Unitful.Quantity) = z
 
 function Unitful.uconvert(u::Unitful.Units{T,Unitful.NoDims}, pt::CartesianZero{<:Quantity{U, Unitful.NoDims}}) where {T,U}
     CartesianZero{promote_type(T,U)}()
@@ -244,6 +255,27 @@ struct CylindricalPoint{T} <: AbstractCylindricalPoint{T}
     CylindricalPoint{T}(r::Real, φ::Real, z::Real) where {T} = new(T(r), mod(T(φ),T(2π)), T(z))
 end
 
+#Units support
+function CylindricalPoint(r, φ, z)
+    if !(r isa Real || r isa Unitful.Length)
+        throw(ArgumentError("Expected `r` to be a length or Real, got unit $(Unitful.unit(r))"))
+    end
+
+    if !(φ isa Real || φ isa Unitful.Quantity{<:Real, NoDims})
+        throw(ArgumentError("Expected `φ` to be an angle or Real, got unit $(Unitful.unit(φ))"))
+    end
+
+    if !(z isa Real || z isa Unitful.Length)
+        throw(ArgumentError("Expected `z` to be a length or Real, got unit $(Unitful.unit(z))"))
+    end
+
+    r_val = to_internal_units(r)
+    φ_val = to_internal_units(φ)
+    z_val = to_internal_units(z)
+
+    return CylindricalPoint(r_val, φ_val, z_val)
+end
+
 function CylindricalPoint(r::TR, φ::TP, z::TZ) where {TR<:Real,TP<:Real,TZ<:Real}
     # ToDo: Simplify this:
     eltypes = _csg_get_promoted_eltype.((TR,TP,TZ))
@@ -303,4 +335,21 @@ 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))
+
+@inline Base.adjoint(pt::CylindricalPoint) = CylindricalPoint(adjoint(pt.r), adjoint(pt.φ), adjoint(pt.z))
+
+function to_internal_units(pt::CartesianPoint)
+    CartesianPoint(to_internal_units(pt.x), to_internal_units(pt.y), to_internal_units(pt.z))
+end
+
+function to_internal_units(pt::CylindricalPoint)
+    _r = to_internal_units(pt.r)
+    _φ = to_internal_units(pt.φ)
+    _z = to_internal_units(pt.z)
+
+    return CartesianPoint(_r * cos(_φ), _r * sin(_φ), _z)
+end
+
+function to_internal_units(pt::AbstractCoordinatePoint)
+    error("Unsupported point type $(typeof(pt)). Expected CartesianPoint or CylindricalPoint.")
+end

src/Event/Event.jl

@@ -1,6 +1,6 @@
 """
-    mutable struct Event{T <: SSDFloat}
-
+  mutable struct Event{T <: SSDFloat}
+         
 Collection struct for individual events. 
 This (mutable) struct is meant to be used to look at individual events,
 not to process a huge amount of events.
@@ -19,30 +19,40 @@ mutable struct Event{T <: SSDFloat}
     Event{T}() where {T <: SSDFloat} = new{T}(VectorOfArrays(Vector{CartesianPoint{T}}[]), VectorOfArrays(Vector{T}[]), missing, missing)
 end
 
-function Event(location::AbstractCoordinatePoint{T}, energy::RealQuantity = one(T))::Event{T} where {T <: SSDFloat}
+function Event(location::AbstractCoordinatePoint, energy::RealQuantity = 1)
+    pt_internal = to_internal_units(location)
+    T = eltype(pt_internal)
     evt = Event{T}()
-    evt.locations = VectorOfArrays([[CartesianPoint(location)]])
-    evt.energies = VectorOfArrays([[T(to_internal_units(energy))]])
+    evt.locations = VectorOfArrays([[pt_internal]])
+    evt.energies  = VectorOfArrays([[T(to_internal_units(energy))]])
     return evt
 end
 
-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}
+function Event(locations::Vector{<:AbstractCoordinatePoint}, energies::Vector{<:RealQuantity} = fill(1, length(locations)); max_interaction_distance::Union{<:Real, <:LengthQuantity} = NaN)
+    pts_internal = to_internal_units.(locations)
+    T = eltype(first(pts_internal))
     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."
     evt = Event{T}()
+
     if isnan(d) # default: no grouping, the charges from different hits drift independently
-        evt.locations = VectorOfArrays(broadcast(pt -> [CartesianPoint(pt)], locations))
-        evt.energies = VectorOfArrays(broadcast(E -> [T(to_internal_units(E))], energies))
+        evt.locations = VectorOfArrays(broadcast(pt -> [pt], pts_internal))
+        evt.energies  = VectorOfArrays(broadcast(E -> [T(to_internal_units(E))], energies))
     else
-        evt.locations, evt.energies = group_points_by_distance(CartesianPoint.(locations), T.(to_internal_units.(energies)), d)
+        evt.locations, evt.energies = group_points_by_distance(pts_internal, T.(to_internal_units.(energies)), d)
     end
     return evt
 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}
+function Event(locations::Vector{<:Vector{<:AbstractCoordinatePoint}}, energies::Vector{<:Vector{<:RealQuantity}} = [[1 for _ in group] for group in locations])
+
+    T = eltype(to_internal_units(first(first(locations))))
+    locs = [[convert(CartesianPoint{T}, to_internal_units(pt)) for pt in group] for group in locations]
+    ens  = [[T(to_internal_units(E)) for E in group] for group in energies]
+
     evt = Event{T}()
-    evt.locations = VectorOfArrays(broadcast(pts -> CartesianPoint{T}.(pts), locations))
-    evt.energies = VectorOfArrays(Vector{T}.(to_internal_units.(energies)))
+    evt.locations = VectorOfArrays(locs)
+    evt.energies  = VectorOfArrays(ens)
     return evt
 end
 
@@ -67,7 +77,6 @@ function Event(locations::Vector{<:AbstractCoordinatePoint{T}}, energies::Vector
         max_interaction_distance::Union{<:Real, <:LengthQuantity} = NaN
     )::Event{T} where {T <: SSDFloat, PT <: ParticleType}
 
-
     @assert eachindex(locations) == eachindex(energies) == eachindex(radius)
 
     d::T = T(to_internal_units(max_interaction_distance))
@@ -99,14 +108,16 @@ function Event(locations::Vector{<:AbstractCoordinatePoint{T}}, energies::Vector
 end
 
 function Event(
-        locations::Vector{<:Vector{<:AbstractCoordinatePoint{T}}}, 
-        energies::Vector{<:Vector{<:RealQuantity}}, N::Int;
-        particle_type::Type{PT} = Gamma, number_of_shells::Int = 2,
-        radius::Vector{<:Vector{<:RealQuantity}} = map(e -> radius_guess.(T.(to_internal_units.(e)), particle_type), energies)
-    ) where {T <: SSDFloat, PT <: ParticleType}
-
-    @assert eachindex(locations) == eachindex(energies) == eachindex(radius)
-    events = map(i -> Event(locations[i], energies[i], N; particle_type, number_of_shells, radius = radius[i]), eachindex(locations))
+    locations::Vector{<:Vector{<:AbstractCoordinatePoint}},
+    energies::Vector{<:Vector{<:RealQuantity}}, N::Int;
+    particle_type::Type = Gamma, number_of_shells::Int = 2,
+    radius::Vector{<:Vector{<:RealQuantity}} = map(e -> radius_guess.(to_internal_units.(e), particle_type), energies)
+)
+    locs_internal = [to_internal_units.(loc) for loc in locations]
+
+    @assert eachindex(locs_internal) == eachindex(energies) == eachindex(radius)
+    events = map(i -> Event(locs_internal[i], energies[i], N; particle_type, number_of_shells, radius = radius[i]), eachindex(locations))
+
     Event(flatview.(getfield.(events, :locations)), flatview.(getfield.(events, :energies)))
 end
 
@@ -375,4 +386,4 @@ function get_electron_and_hole_contribution(evt::Event{T}, sim::Simulation{T, S}
             hole_contribution = RDWaveform(range(zero(T) * u"ns", step = dt * u"ns", length = length(signal_h)), signal_h * calibration_factor))
 end
 
-export get_electron_and_hole_contribution
+export get_electron_and_hole_contribution