diff --git a/docs/src/man/charge_drift.md b/docs/src/man/charge_drift.md index 2b4c2392e..94fbb4dcc 100644 --- a/docs/src/man/charge_drift.md +++ b/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) diff --git a/docs/src/tutorials/geant4_ssd_lit.jl b/docs/src/tutorials/geant4_ssd_lit.jl index 57f2cc54e..3a7f913d3 100644 --- a/docs/src/tutorials/geant4_ssd_lit.jl +++ b/docs/src/tutorials/geant4_ssd_lit.jl @@ -89,7 +89,7 @@ events = run_geant4_simulation(app, N_events) plot(sim.detector, show_passives = false, size = (500,500), fmt = :png) plot!(source_1) -plot!(ustrip.(u"m", flatview(events[1:1000].pos)), ms = 0.5, msw = 0, color=:black, label = "") +plot!(flatview(events[1:1000].pos), ms = 0.5, msw = 0, color=:black, label = "") #jl savefig("events.pdf") # hide #md savefig("events.pdf") # hide #md savefig("events.svg"); nothing # hide @@ -150,4 +150,4 @@ plot(w[1:20], label = "") #jl savefig("wf_and_amplitude.pdf") # hide #md savefig("wf_and_amplitude.pdf") # hide #md savefig("wf_and_amplitude.svg"); nothing # hide -#md # [![wf_and_amplitude](wf_and_amplitude.svg)](wf_and_amplitude.pdf) \ No newline at end of file +#md # [![wf_and_amplitude](wf_and_amplitude.svg)](wf_and_amplitude.pdf) diff --git a/ext/SolidStateDetectorsLegendHDF5IOExt.jl b/ext/SolidStateDetectorsLegendHDF5IOExt.jl index 7147e71ee..bad97ae46 100644 --- a/ext/SolidStateDetectorsLegendHDF5IOExt.jl +++ b/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 @@ -55,14 +57,107 @@ function SolidStateDetectors.simulate_waveforms( mcevents::AbstractVector{<:Name @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, signal_unit, max_interaction_distance, end_drift_when_no_field, geometry_check, verbose) - # LH5 can't handle CartesianPoint, turn positions into CartesianVectors which will be saved as SVectors - pos_vec = VectorOfVectors([[CartesianPoint(p...) - cartesian_zero for p in ps] for ps in mcevents_sub.pos]) - new_mcevents_sub = TypedTables.Table(merge(Tables.columns(mcevents_sub), (pos = pos_vec,))) + # # LH5 couldn't handle CartesianPoint, turn positions into CartesianVectors which will be saved as SVectors + # pos_vec = VectorOfVectors([[CartesianPoint(p...) - cartesian_zero for p in ps] for ps in mcevents_sub.pos]) + # new_mcevents_sub = TypedTables.Table(merge(Tables.columns(mcevents_sub), (pos = pos_vec,))) LegendHDF5IO.lh5open(ofn, "w") do h5f - LegendHDF5IO.writedata(h5f.data_store, "generated_waveforms", new_mcevents_sub) + LegendHDF5IO.writedata(h5f.data_store, "generated_waveforms", mcevents_sub) end 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 diff --git a/src/ChargeCloudModels/ChargeCloudModels.jl b/src/ChargeCloudModels/ChargeCloudModels.jl index 1f0a9e542..7f555fb66 100644 --- a/src/ChargeCloudModels/ChargeCloudModels.jl +++ b/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 \ No newline at end of file +end diff --git a/src/ChargeClustering/ChargeClustering.jl b/src/ChargeClustering/ChargeClustering.jl index 8b03f78e7..8b8bbb68d 100644 --- a/src/ChargeClustering/ChargeClustering.jl +++ b/src/ChargeClustering/ChargeClustering.jl @@ -5,7 +5,7 @@ function cluster_detector_hits( edep::AbstractVector{TT}, pos::AbstractVector{<:Union{<:StaticVector{3,PT}, <:CartesianPoint{PT}}}, cluster_radius::RealQuantity - ) where {TT <: RealQuantity, PT <: RealQuantity} + ) where {TT <: RealQuantity, T <: Real, PT <: Union{T, <:Unitful.Length{T}}} unsorted = TypedTables.Table(detno = detno, edep = edep, pos = pos) sorting_idxs = sortperm(unsorted.detno) @@ -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 = _parse_value(float(T), cluster_radius, internal_length_unit) for d_hits_nt in grouped d_hits = TypedTables.Table(d_hits_nt) @@ -25,7 +24,7 @@ function cluster_detector_hits( @assert all(isequal(d_detno), d_hits.detno) if length(d_hits) > 3 - clusters = Clustering.dbscan(hcat((ustrip.(getindex.(d_hits.pos,i)) for i in 1:3)...)', + clusters = Clustering.dbscan(hcat((to_internal_units.(getindex.(d_hits.pos,i)) for i in 1:3)...)', ustripped_cradius, leafsize = 20, min_neighbors = 1, min_cluster_size = 1).clusters for c in clusters idxs = vcat(c.boundary_indices, c.core_indices) diff --git a/src/ConstructiveSolidGeometry/PointsAndVectors/Points.jl b/src/ConstructiveSolidGeometry/PointsAndVectors/Points.jl index 442d5cecf..153a073aa 100644 --- a/src/ConstructiveSolidGeometry/PointsAndVectors/Points.jl +++ b/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) @@ -62,6 +56,21 @@ struct CartesianPoint{T} <: AbstractCartesianPoint{T} x::T y::T z::T + CartesianPoint{T}(x::T, y::T, z::T) where {T<:AbstractFloat} = new(x, y, z) + CartesianPoint{T}(x::Real, y::Real, z::Real) where {T} = new(T(x), T(y), T(z)) +end + +#Units support +function CartesianPoint(x, y, z) + for (name, pt) in zip((:x, :y, :z), (x, y, z)) + (pt isa Real || pt isa Unitful.Length) || + throw(ArgumentError("Expected $(name) to be a length or Real, got unit $(Unitful.unit(pt))")) + end + x_val = to_internal_units(x) + y_val = to_internal_units(y) + z_val = to_internal_units(z) + + return CartesianPoint(x_val, y_val, z_val) end #Type promotion happens here @@ -204,7 +213,9 @@ 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)}}() +# ToDo: Revert this once we support units internally. +#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)) \ No newline at end of file + +@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 diff --git a/src/Event/Event.jl b/src/Event/Event.jl index eed4f9c3d..182e01381 100644 --- a/src/Event/Event.jl +++ b/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 \ No newline at end of file +export get_electron_and_hole_contribution diff --git a/src/MCEventsProcessing/MCEventsProcessing.jl b/src/MCEventsProcessing/MCEventsProcessing.jl index 6e309599a..51c435bf6 100644 --- a/src/MCEventsProcessing/MCEventsProcessing.jl +++ b/src/MCEventsProcessing/MCEventsProcessing.jl @@ -116,7 +116,6 @@ function simulate_waveforms( mcevents_table::AbstractVector{<:NamedTuple}, sim:: return vcat(mcevents_chns...) end - # Support detectors hits with positions given as vectors: function _get_unitless_positions(point_vectors::AbstractVector{<:StaticVector{3}}) Ref(cartesian_zero) .+ to_internal_units.(point_vectors) @@ -190,7 +189,7 @@ function _simulate_charge_drifts( mcevents::TypedTables.Table, sim::Simulation{T @assert is_detector_hits_table(mcevents) "Table does not have the correct format" @showprogress map(mcevents) do phyevt - locations, edeps = _convertEnergyDepsToChargeDeps(phyevt.pos, phyevt.edep, sim.detector; number_of_carriers, number_of_shells, max_interaction_distance, verbose) + locations, edeps = _convertEnergyDepsToChargeDeps(to_internal_units.(phyevt.pos), phyevt.edep, sim.detector; number_of_carriers, number_of_shells, max_interaction_distance, verbose) 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, diffusion, self_repulsion, end_drift_when_no_field, geometry_check, verbose diff --git a/test/ConstructiveSolidGeometry/CSG_IO.jl b/test/ConstructiveSolidGeometry/CSG_IO.jl index 17da0a73b..c6e341a84 100644 --- a/test/ConstructiveSolidGeometry/CSG_IO.jl +++ b/test/ConstructiveSolidGeometry/CSG_IO.jl @@ -4,8 +4,11 @@ using SolidStateDetectors import SolidStateDetectors.ConstructiveSolidGeometry as CSG import SolidStateDetectors.ConstructiveSolidGeometry: Geometry -T = Float64 +import LegendHDF5IO +import Unitful: @u_str +import ArraysOfArrays: VectorOfVectors +T = Float64 example_primitive_dir = joinpath(@__DIR__, "../../examples/example_primitive_files") @testset "Test primitive read-in" begin @@ -44,4 +47,64 @@ example_primitive_dir = joinpath(@__DIR__, "../../examples/example_primitive_fil torus = Geometry(T, joinpath(example_primitive_dir, "Torus.yaml")) @test torus isa CSG.Torus{T} end +end + +@testset "Test LegendHDF5IO with AbstractCoordinatePoint" begin + x = rand(T, 3) + x_unit = x * u"mm" + pt = CartesianPoint(x...) + pt_unit = CartesianPoint(x_unit...) + pt_cyl = CylindricalPoint(x...) + pt_cyl_unit = CylindricalPoint(x .* [u"mm", u"ยฐ", u"mm"]...) + + @testset begin + mktemp() do tmpfile, io + for x in (x, x_unit, pt, pt_unit, pt_cyl) + LegendHDF5IO.lh5open(tmpfile, "w") do h + + @testset "I/O for AbstractCoordinatePoints" begin + # write to file using writedata and LH5Array + @test_nowarn LegendHDF5IO.writedata(h, "x1", x) + @test_nowarn h["x2"] = x + + # read using readddata + x1 = @test_nowarn LegendHDF5IO.readdata(h, "x1") + @test x1 isa typeof(x) + @test x1 == x + + # read using LH5Array + x2 = @test_nowarn h["x2"] + @test x2 == x + + # cross-compatibility between the two methods + x1_ = @test_nowarn h["x1"] + x2_ = @test_nowarn LegendHDF5IO.readdata(h, "x2") + @test x1_ == x2_ == x + end + + @testset "I/O for Vector of AbstractCoordinatePoints" begin + + y = x isa Vector ? VectorOfVectors([x]) : [x] + + # write to file using writedata and LH5Array + @test_nowarn LegendHDF5IO.writedata(h, "y1", y) + @test_nowarn h["y2"] = y + + # read using readddata + y1 = @test_nowarn LegendHDF5IO.readdata(h, "y1") + @test y1 isa typeof(y) + @test y == y1 + + # read using LH5Array + @test y == h["y2"] + + # cross-compatibility between the two methods + y1_ = @test_nowarn h["y1"] + y2_ = @test_nowarn LegendHDF5IO.readdata(h, "y2") + @test y1_ == y2_ == y + end + end + end + end + end end \ No newline at end of file diff --git a/test/ConstructiveSolidGeometry/CSG_coordinates.jl b/test/ConstructiveSolidGeometry/CSG_coordinates.jl index 38078c7df..9a155eac4 100644 --- a/test/ConstructiveSolidGeometry/CSG_coordinates.jl +++ b/test/ConstructiveSolidGeometry/CSG_coordinates.jl @@ -7,7 +7,7 @@ using SolidStateDetectors.ConstructiveSolidGeometry: CartesianPoint, CartesianVe using StaticArrays: Size, SVector, SMatrix using InverseFunctions: inverse -import Unitful +import Unitful: @u_str @testset "points_and_vectors" begin @testset "cartesian" begin @@ -36,7 +36,6 @@ import Unitful @test @inferred(CartesianPoint(a[1], a[2], a[3])) == a @test @inferred(CartesianPoint(a[1], a[2], a[3])) โ‰ˆ a - frame = LocalAffineFrame(b, A) f = frame_transformation(frame, global_frame) @@ -44,13 +43,36 @@ import Unitful @test @inferred(inverse(f)(f(a))) โ‰ˆ a - @test @inferred(CartesianZero{Float32}() * Unitful.u"mm") === CartesianZero{typeof(zero(Float32) * Unitful.u"mm")}() + #@test @inferred(CartesianZero{Float32}() * u"mm") === CartesianZero{typeof(zero(Float32) * u"mm")}() A = [CartesianPoint{Float32}(x,0,0) for x in -2:2] @test isapprox(barycenter(A), CartesianPoint{Float32}(0,0,0)) S = SVector{length(A)}(A) @test isapprox(barycenter(S), CartesianPoint{Float32}(0,0,0)) + + # test types and units + @test CartesianPoint(1, 2, 3) isa CartesianPoint{Float64} + @test CartesianPoint(1, 2, 3f0) isa CartesianPoint{Float32} + @test CartesianPoint(1, 2.0f0, Float16(3)) isa CartesianPoint{Float32} + @test CartesianPoint(1.0, 2.0f0, Float16(3)) isa CartesianPoint{Float64} + @test CartesianPoint(1, 2, Float16(3)) isa CartesianPoint{Float16} + @test CartesianPoint(1.0u"m", 2.0u"m", 3.0f0u"m") isa CartesianPoint{Float64} + @test CartesianPoint(1.0u"m", 2.0f0u"m", 3.0f0u"m") isa CartesianPoint{Float64} + @test CartesianPoint(1.0f0u"m", 2.0f0u"m", 3.0f0u"m") isa CartesianPoint{Float32} + @test CartesianPoint(1u"m", 2u"m", 3u"m") isa CartesianPoint{Float64} + @test CartesianPoint(1.0u"mm", 2.0u"cm", 3.0f0u"m") isa CartesianPoint{Float64} + @test CartesianPoint(1.0f0u"mm", 2.0f0u"mm", 3.0u"m") isa CartesianPoint{Float64} + @test CartesianPoint(1.0u"mm", 2.0f0u"cm", 3.0f0u"m") isa CartesianPoint{Float64} + @test CartesianPoint(1.0f0u"mm", 2.0f0u"cm", 3.0f0u"m") isa CartesianPoint{Float32} + @test CartesianPoint(1u"mm", 2u"cm", 3u"m") isa CartesianPoint{Float64} + @test CartesianPoint(1u"mm", 0, 0) isa CartesianPoint{Float64} + @test CartesianPoint(1u"m", 2u"m", 3f0u"m") isa CartesianPoint{Float32} + + # test throwing errors with wrong units + @test_throws ArgumentError CartesianPoint(1u"m", 2u"rad", 3u"m") + @test_throws ArgumentError CartesianPoint(1u"s", 2u"m", 3u"m") + @test_throws ArgumentError CartesianPoint(1u"m", 2u"m", 3u"kg") end @testset "cylindrical" begin @@ -80,5 +102,28 @@ import Unitful S = SVector{length(A)}(A) @test isapprox(barycenter(S), CylindricalPoint{Float32}(0,0,0)) + + # test types and units + @test CylindricalPoint(1, 2, 3) isa CylindricalPoint{Float64} + @test CylindricalPoint(1, 2, 3f0) isa CylindricalPoint{Float32} + @test CylindricalPoint(1, 2.0f0, Float16(3)) isa CylindricalPoint{Float32} + @test CylindricalPoint(1.0, 2.0f0, Float16(3)) isa CylindricalPoint{Float64} + @test CylindricalPoint(1, 2, Float16(3)) isa CylindricalPoint{Float16} + @test CylindricalPoint(1.0u"m", 2.0u"rad", 3.0f0u"m") isa CylindricalPoint{Float64} + @test CylindricalPoint(1.0u"m", 2.0f0u"rad", 3.0f0u"m") isa CylindricalPoint{Float64} + @test CylindricalPoint(1.0f0u"m", 2.0f0u"rad", 3.0f0u"m") isa CylindricalPoint{Float32} + @test CylindricalPoint(1u"m", 2u"rad", 3u"m") isa CylindricalPoint{Float64} + @test CylindricalPoint(1.0u"mm", 2.0u"rad", 3.0f0u"m") isa CylindricalPoint{Float64} + @test CylindricalPoint(1.0f0u"mm", 2.0f0u"rad", 3.0u"m") isa CylindricalPoint{Float64} + @test CylindricalPoint(1.0u"mm", 2.0f0u"rad", 3.0f0u"m") isa CylindricalPoint{Float64} + @test CylindricalPoint(1.0f0u"mm", 2.0f0u"rad", 3.0f0u"m") isa CylindricalPoint{Float32} + @test CylindricalPoint(1u"mm", 2u"rad", 3u"m") isa CylindricalPoint{Float64} + @test CylindricalPoint(1u"mm", 0, 0) isa CylindricalPoint{Float64} + @test CylindricalPoint(1u"m", 2u"rad", 3f0u"m") isa CylindricalPoint{Float32} + + # test throwing errors with wrong units + @test_throws ArgumentError CylindricalPoint(1u"rad", 2u"rad", 3u"m") + @test_throws ArgumentError CylindricalPoint(1u"m", 2u"m", 3u"m") + @test_throws ArgumentError CylindricalPoint(1u"m", 2u"rad", 3u"rad") end end diff --git a/test/test_charge_drift_models.jl b/test/test_charge_drift_models.jl index 385197838..c234af0ea 100644 --- a/test/test_charge_drift_models.jl +++ b/test/test_charge_drift_models.jl @@ -565,3 +565,48 @@ end end end end + +@timed_testset "NBodyChargeCloud Units" begin + + cartesian_unit = CartesianPoint(0.01u"m", 0.0u"m", 0.05u"m") + cartesian_float = CartesianPoint(0.01, 0.0, 0.05) + cylindrical_unit = CylindricalPoint(0.01u"m", (ฯ€/4)u"rad", 0.05u"m") + cylindrical_float = CylindricalPoint(0.01, ฯ€/4, 0.05) + Edep = 1u"MeV" + + # Variant 1: CartesianPoint with Units, energy + nb1 = NBodyChargeCloud(cartesian_unit, Edep) + @test all(x -> x isa CartesianPoint, nb1.locations) + @test isapprox(sum(nb1.energies), ustrip.(u"eV", Edep)) + + # Variant 2: CylindricalPoint with Units, energy + nb2 = NBodyChargeCloud(cylindrical_unit, Edep) + @test all(x -> x isa CartesianPoint, nb2.locations) + @test isapprox(sum(nb2.energies), ustrip.(u"eV", Edep)) + + # Variant 3: CartesianPoint with Units, energy, N + nb3 = NBodyChargeCloud(cartesian_unit, Edep, 10, radius = 0.001, number_of_shells = 1) + @test all(x -> x isa CartesianPoint, nb3.locations) + @test length(nb3.locations) >= 10 + @test isapprox(sum(nb3.energies), ustrip.(u"eV", Edep)) + + # Variant 4: CylindricalPoint with Units, energy, N + nb4 = NBodyChargeCloud(cylindrical_unit, Edep, 10, radius = 0.001, number_of_shells = 1) + @test all(x -> x isa CartesianPoint, nb4.locations) + @test length(nb4.locations) >= 10 + @test isapprox(sum(nb4.energies), ustrip.(u"eV", Edep)) + + # Edge case, zero radius + nb_zero = NBodyChargeCloud(cartesian_float, Edep, 5, radius = 0.0, number_of_shells = 1) + @test all(x -> x isa CartesianPoint, nb_zero.locations) + @test all(x -> x == nb_zero.locations[1], nb_zero.locations) # all points coincide at the center + @test isapprox(sum(nb_zero.energies), ustrip.(u"eV", Edep)) + @test all(x -> all(isfinite, (x.x, x.y, x.z)), nb_zero.locations) + + nb_zero_cyl = NBodyChargeCloud(cylindrical_float, Edep, 5, radius = 0.0, number_of_shells = 1) + @test all(x -> x isa CartesianPoint, nb_zero_cyl.locations) + @test all(x -> x == nb_zero_cyl.locations[1], nb_zero_cyl.locations) + @test isapprox(sum(nb_zero_cyl.energies), ustrip.(u"eV", Edep)) + @test all(x -> all(isfinite, (x.x, x.y, x.z)), nb_zero_cyl.locations) + +end diff --git a/test/test_geant4.jl b/test/test_geant4.jl index 23913f5ba..bcbb6fe5e 100644 --- a/test/test_geant4.jl +++ b/test/test_geant4.jl @@ -8,6 +8,7 @@ using RadiationDetectorSignals using StaticArrays using TypedTables using Unitful +using LegendHDF5IO import Geant4 @@ -82,6 +83,18 @@ end @test length(flatview(clustered_evts.pos)) <= length(flatview(evts.pos)) @test eltype(first(clustered_evts.pos)) <: CartesianPoint + # Temporary HDF5 file + mktemp() do tmpfile, io + LegendHDF5IO.lh5open(tmpfile, "w") do h5 + LegendHDF5IO.writedata(h5, "geant4", evts) + end + @test isfile(tmpfile) + + evts == LegendHDF5IO.lh5open(tmpfile) do h5 + LegendHDF5IO.readdata(h5, "geant4") + end + end + # 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) @@ -90,22 +103,29 @@ end @test length(wf) == length(evts) * sum(.!ismissing.(sim.weighting_potentials)) # Use the table with added Fano noise - wf = simulate_waveforms(evts_fano, sim, ฮ”t = 1u"ns", max_nsteps = 2000) - @test wf isa Table - @test :waveform in columnnames(wf) - @test length(wf) == length(evts_fano) * sum(.!ismissing.(sim.weighting_potentials)) + wf_fano = simulate_waveforms(evts_fano, sim, ฮ”t = 1u"ns", max_nsteps = 2000) + @test wf_fano isa Table + @test :waveform in columnnames(wf_fano) + @test length(wf_fano) == length(evts_fano) * sum(.!ismissing.(sim.weighting_potentials)) - # Try the same method using StaticVectors as eltype of pos - evts_static = Table(evts; pos = VectorOfVectors(broadcast.(p -> SVector{3}(p.x, p.y, p.z), evts.pos))) + # Try the same method using StaticVectors as eltype of pos (with units of mm) + evts_static = Table(evts; pos = VectorOfVectors(broadcast.(p -> SVector{3}(p.x, p.y, p.z) * 1000u"mm", SolidStateDetectors.to_internal_units(evts.pos)))) clustered_evts_static = SolidStateDetectors.cluster_detector_hits(evts_static, 10u"ยตm") @test length(clustered_evts_static) == length(evts_static) @test length(flatview(clustered_evts_static.pos)) <= length(flatview(evts_static.pos)) @test eltype(first(clustered_evts_static.pos)) <: StaticVector{3} + # Check clustering with cluster radius being unitless or with the wrong unit + clustered_evts_static_unitless = SolidStateDetectors.cluster_detector_hits(evts_static, 1e-5) + @test clustered_evts_static == clustered_evts_static_unitless + @test_throws Exception SolidStateDetectors.cluster_detector_hits(evts_static, 2u"kg") + # Generate waveforms using StaticVectors as eltype of pos wf_static = simulate_waveforms(evts_static, sim, ฮ”t = 1u"ns", max_nsteps = 2000) @test wf_static isa Table @test :waveform in columnnames(wf_static) @test length(wf_static) == length(evts_static) * sum(.!ismissing.(sim.weighting_potentials)) -end \ No newline at end of file + # Result should not depend whether input positions are SVector or CartesianPoint + @test wf_static.waveform โ‰ˆ wf.waveform +end diff --git a/test/test_io.jl b/test/test_io.jl index 939851249..17284934a 100644 --- a/test/test_io.jl +++ b/test/test_io.jl @@ -64,51 +64,55 @@ end sim = Simulation{T}(SSD_examples[:InvertedCoax]) timed_simulate!(sim, convergence_limit = 1e-5, device_array_type = device_array_type, refinement_limits = [0.2, 0.1], verbose = false) - evt_table = Table( - evtno = Int32[1], - detno = Int32[1], - thit = VectorOfVectors([T[0] * u"s"]), - edep = VectorOfVectors([T[1] * u"eV"]), - pos = VectorOfVectors([[SVector{3, T}.(0.01, 0.01, 0.01) * u"m"]]) - ) - - contact_charge_signals = timed_simulate_waveforms( - evt_table, - sim, - max_nsteps = 4000, - ฮ”t = 1u"ns", - number_of_carriers = 20, - number_of_shells = 2, - geometry_check = true, - verbose = false); - - @timed_testset "Regular simulate_waveforms" begin - signalsum = T(0) - for i in 1:length(contact_charge_signals.waveform) - signalsum += abs(ustrip(contact_charge_signals.waveform[i].signal[end])) + # test I/O for both SVector (with and without units) and CartesianPoint + for pos in ([SVector{3, T}.(10, 10, 10) * u"mm"], [SVector{3, T}.(0.01, 0.01, 0.01)], [CartesianPoint{T}(0.01, 0.01, 0.01)]) + + evt_table = Table( + evtno = Int32[1], + detno = Int32[1], + thit = VectorOfVectors([T[0] * u"s"]), + edep = VectorOfVectors([T[1] * u"eV"]), + pos = VectorOfVectors([pos]) + ) + + contact_charge_signals = timed_simulate_waveforms( + evt_table, + sim, + max_nsteps = 4000, + ฮ”t = 1u"ns", + number_of_carriers = 20, + number_of_shells = 2, + geometry_check = true, + verbose = false); + + @timed_testset "Regular simulate_waveforms" begin + signalsum = T(0) + for i in 1:length(contact_charge_signals.waveform) + signalsum += abs(ustrip(contact_charge_signals.waveform[i].signal[end])) + end + signalsum *= inv(ustrip(SolidStateDetectors._convert_internal_energy_to_external_charge(sim.detector.semiconductor.material))) + @test isapprox( signalsum, T(2), atol = 5e-3 ) end - signalsum *= inv(ustrip(SolidStateDetectors._convert_internal_energy_to_external_charge(sim.detector.semiconductor.material))) - @test isapprox( signalsum, T(2), atol = 5e-3 ) - end - @timed_testset "LegendHDF5IO simulate_waveforms" begin - timed_simulate_waveforms( - evt_table, - sim, - ".", - chunk_n_physics_events = 1, - max_nsteps = 4000, - ฮ”t = 1u"ns", - number_of_carriers = 20, - number_of_shells = 2, - geometry_check = true, - verbose = false - ) - @info isfile("generated_waveforms_evts_1-1.h5") - @test contact_charge_signals == LegendHDF5IO.lh5open("generated_waveforms_evts_1-1.h5") do h5f - @test haskey(h5f, "generated_waveforms") - LegendHDF5IO.readdata(h5f.data_store, "generated_waveforms") + @timed_testset "LegendHDF5IO simulate_waveforms" begin + timed_simulate_waveforms( + evt_table, + sim, + ".", + chunk_n_physics_events = 1, + max_nsteps = 4000, + ฮ”t = 1u"ns", + number_of_carriers = 20, + number_of_shells = 2, + geometry_check = true, + verbose = false + ) + @info isfile("generated_waveforms_evts_1-1.h5") + @test contact_charge_signals == LegendHDF5IO.lh5open("generated_waveforms_evts_1-1.h5") do h5f + @test haskey(h5f, "generated_waveforms") + LegendHDF5IO.readdata(h5f.data_store, "generated_waveforms") + end + rm("generated_waveforms_evts_1-1.h5") end - rm("generated_waveforms_evts_1-1.h5") end -end \ No newline at end of file +end