From 624a8695eca87aa2e411a86889c0c0bc635ae7ca Mon Sep 17 00:00:00 2001
From: Claudia Alvarez Garcia <claudiaalvgar@hotmail.es>
Date: Mon, 10 Nov 2025 09:38:59 +0100
Subject: [PATCH 01/39] Function to define CylindricalPoint with units

---
 .../PointsAndVectors/Points.jl                        | 11 +++++++++++
 1 file changed, 11 insertions(+)

diff --git a/src/ConstructiveSolidGeometry/PointsAndVectors/Points.jl b/src/ConstructiveSolidGeometry/PointsAndVectors/Points.jl
index 442d5cecf..6cc3d2d3a 100644
--- a/src/ConstructiveSolidGeometry/PointsAndVectors/Points.jl
+++ b/src/ConstructiveSolidGeometry/PointsAndVectors/Points.jl
@@ -244,6 +244,17 @@ struct CylindricalPoint{T} <: AbstractCylindricalPoint{T}
     CylindricalPoint{T}(r::Real, φ::Real, z::Real) where {T} = new(T(r), mod(T(φ),T(2π)), T(z))
 end
 
+# Handles Unitful quantities or plain numbers 
+function CylindricalPoint(r, φ, z)
+    r_val = isa(r, Unitful.Quantity) ? float(ustrip(uconvert(u"m", r))) : float(r)
+    φ_val = isa(φ, Unitful.Quantity) ? float(ustrip(uconvert(u"rad", φ))) : float(φ)
+    z_val = isa(z, Unitful.Quantity) ? float(ustrip(uconvert(u"m", z))) : float(z)
+
+    T = promote_type(typeof(r_val), typeof(φ_val), typeof(z_val))
+    CylindricalPoint{T}(T(r_val), T(φ_val), T(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))

From cb2e60300c03868a6e1d60987dd1ae5468b065d7 Mon Sep 17 00:00:00 2001
From: Claudia Alvarez Garcia <claudiaalvgar@hotmail.es>
Date: Mon, 10 Nov 2025 09:40:23 +0100
Subject: [PATCH 02/39] Allow NBodyChargeCloud to receive CylindricalPoint in
 Floats and with units and CartesianPoint with units

---
 src/ChargeCloudModels/ChargeCloudModels.jl | 38 +++++++++++++++++++++-
 1 file changed, 37 insertions(+), 1 deletion(-)

diff --git a/src/ChargeCloudModels/ChargeCloudModels.jl b/src/ChargeCloudModels/ChargeCloudModels.jl
index 1f0a9e542..07950fadc 100644
--- a/src/ChargeCloudModels/ChargeCloudModels.jl
+++ b/src/ChargeCloudModels/ChargeCloudModels.jl
@@ -52,6 +52,24 @@ 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::CartesianPoint{<:Unitful.Quantity}, energy::RealQuantity;
+                          kwargs...)
+    internal_center = to_internal_point(center)
+    return NBodyChargeCloud(internal_center, energy; kwargs...)
+end
+
+function NBodyChargeCloud(center::CylindricalPoint{<:Unitful.Quantity}, energy::RealQuantity;
+                          kwargs...)
+    internal_center = to_internal_point(center)
+    return NBodyChargeCloud(internal_center, energy; kwargs...)
+end
+
+function NBodyChargeCloud(center::CylindricalPoint{<:SSDFloat}, energy::RealQuantity;
+                          kwargs...)
+    internal_center = to_internal_point(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 +139,24 @@ 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::CartesianPoint{<:Unitful.Quantity}, energy::RealQuantity, N::Integer;
+                          kwargs...)
+    internal_center = to_internal_point(center)
+    return NBodyChargeCloud(internal_center, energy; kwargs...)
+end
+
+function NBodyChargeCloud(center::CylindricalPoint{<:Unitful.Quantity}, energy::RealQuantity, N::Integer;
+                          kwargs...)
+    internal_center = to_internal_point(center)
+    return NBodyChargeCloud(internal_center, energy; kwargs...)
+end
+
+function NBodyChargeCloud(center::CylindricalPoint{<:SSDFloat}, energy::RealQuantity, N::Integer;
+                          kwargs...)
+    internal_center = to_internal_point(center)
+    return NBodyChargeCloud(internal_center, energy; 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 +178,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

From fa00d52dea9bf2de9791347d1143f8da3fe269d7 Mon Sep 17 00:00:00 2001
From: Claudia Alvarez Garcia <claudiaalvgar@hotmail.es>
Date: Mon, 10 Nov 2025 09:41:26 +0100
Subject: [PATCH 03/39] Use the defined function to_internal_point to convert
 Points to Cartesian when passed to Events

---
 src/Event/Event.jl | 85 +++++++++++++++++++++++++++++++++-------------
 1 file changed, 62 insertions(+), 23 deletions(-)

diff --git a/src/Event/Event.jl b/src/Event/Event.jl
index eed4f9c3d..e25c369ee 100644
--- a/src/Event/Event.jl
+++ b/src/Event/Event.jl
@@ -1,6 +1,30 @@
 """
-    mutable struct Event{T <: SSDFloat}
-
+  to_internal_point(pt::AbstractCoordinatePoint)
+        
+Converts any `AbstractCoordinatePoint` (Cartesian or Cylindrical, with or without units) to a `CartesianPoint`.
+ """
+function to_internal_point(pt::AbstractCoordinatePoint)
+    if pt isa CartesianPoint
+        T = pt.x isa Unitful.Quantity ? float(typeof(ustrip(pt.x))) : typeof(pt.x)
+        x = pt.x isa Unitful.Quantity ? float(to_internal_units(pt.x)) : pt.x
+        y = pt.y isa Unitful.Quantity ? float(to_internal_units(pt.y)) : pt.y
+        z = pt.z isa Unitful.Quantity ? float(to_internal_units(pt.z)) : pt.z
+        return CartesianPoint{T}(x, y, z)
+        
+    elseif pt isa CylindricalPoint
+        T = pt.r isa Unitful.Quantity ? float(typeof(ustrip(pt.r))) : float(typeof(pt.r))
+        r = pt.r isa Unitful.Quantity ? float(to_internal_units(pt.r)) : float(pt.r)
+        φ = pt.φ isa Unitful.Quantity ? float(ustrip(uconvert(u"rad", pt.φ))) : float(pt.φ)
+        z = pt.z isa Unitful.Quantity ? float(to_internal_units(pt.z)) : float(pt.z)
+        return CartesianPoint{T}(r * cos(φ), r * sin(φ), z)
+        
+     else
+        error("Unsupported point type $(typeof(pt)). Expected CartesianPoint or CylindricalPoint.")
+     end
+end
+"""
+  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 +43,43 @@ 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}
-    evt = Event{T}()
-    evt.locations = VectorOfArrays([[CartesianPoint(location)]])
-    evt.energies = VectorOfArrays([[T(to_internal_units(energy))]])
+function Event(location::AbstractCoordinatePoint{T}, energy::RealQuantity = 1.0) where {T <: RealQuantity}
+    pt_internal = to_internal_point(location)
+    T_ = typeof(pt_internal.x)
+    evt = Event{T_}()
+    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}
-    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))
+function Event(locations::Vector{<:AbstractCoordinatePoint{T}}, energies::Vector{<:RealQuantity} = fill(1.0, length(locations)); max_interaction_distance::Union{<:Real, <:LengthQuantity} = NaN) where {T <: RealQuantity}
+    pts_internal = [to_internal_point(pt) for pt in locations]
+    T_ = typeof(first(pts_internal).x)
+    d::T_ = T_(to_internal_units(max_interaction_distance))
+    @assert isnan(d) || d >= 0
+
+    evt = Event{T_}()
+    
+    if isnan(d)
+        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}
-    evt = Event{T}()
-    evt.locations = VectorOfArrays(broadcast(pts -> CartesianPoint{T}.(pts), locations))
-    evt.energies = VectorOfArrays(Vector{T}.(to_internal_units.(energies)))
+function Event(locations::Vector{<:Vector{<:AbstractCoordinatePoint{T}}}, energies::Vector{<:Vector{<:RealQuantity}} = [[1.0 for _ in group] for group in locations]) where {T <: RealQuantity}
+    locs_internal = [[to_internal_point(pt) for pt in group] for group in locations]
+    ens_internal  = [[to_internal_units(E) for E in group] for group in energies]
+    T_ = typeof(first(first(locs_internal)).x)
+    
+    locs_T = [[convert(CartesianPoint{T_}, pt) for pt in group] for group in locs_internal]
+    ens_T  = [[T_(E) for E in group] for group in ens_internal]
+
+    evt = Event{T_}()
+    evt.locations = VectorOfArrays(locs_T)
+    evt.energies  = VectorOfArrays(ens_T)
     return evt
 end
 
@@ -67,7 +104,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))
@@ -103,10 +139,13 @@ function Event(
         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}
+    ) where {T <: RealQuantity, PT <: ParticleType}
+
+    locs_internal = [to_internal_point.(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))
     
-    @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))
     Event(flatview.(getfield.(events, :locations)), flatview.(getfield.(events, :energies)))
 end
 
@@ -375,4 +414,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

From 1f999ac234df405c210fa217cc794585b2f0c815 Mon Sep 17 00:00:00 2001
From: Claudia Alvarez Garcia <claudiaalvgar@hotmail.es>
Date: Mon, 10 Nov 2025 10:51:23 +0100
Subject: [PATCH 04/39] Fix manual with new behaviour

---
 docs/src/man/charge_drift.md | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

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)

From 9482e5d608c947f9b5303c5fb9a14ff001ff958e Mon Sep 17 00:00:00 2001
From: Claudia Alvarez Garcia <claudiaalvgar@hotmail.es>
Date: Wed, 12 Nov 2025 12:55:42 +0100
Subject: [PATCH 05/39] convert CartesianPoint to SVector to be returned from
 run_geant4_simulation fixes bug

---
 ext/SolidStateDetectorsGeant4Ext.jl | 7 +++++--
 1 file changed, 5 insertions(+), 2 deletions(-)

diff --git a/ext/SolidStateDetectorsGeant4Ext.jl b/ext/SolidStateDetectorsGeant4Ext.jl
index feabfb2c1..b07a4ea7b 100644
--- a/ext/SolidStateDetectorsGeant4Ext.jl
+++ b/ext/SolidStateDetectorsGeant4Ext.jl
@@ -220,9 +220,12 @@ function SolidStateDetectors.run_geant4_simulation(app::G4JLApplication, number_
         evtno += 1
     end
 
-    evts_with_points = map(evt -> merge(evt, (pos = cartesian_zero + evt.pos,)), evts)
+    evts_with_points = map(evts) do evt
+        merge(evt, (pos = SolidStateDetectors.to_internal_point(evt.pos),))
+    end
 
-    return RadiationDetectorSignals.group_by_evtno(Table(evts_with_points))
+    tbl = Table(evts_with_points)
+    return RadiationDetectorSignals.group_by_evtno(tbl)
 end 
 
 

From 80be3b45c940e70a730d82e85c978ce43b7582ac Mon Sep 17 00:00:00 2001
From: Claudia Alvarez Garcia <claudiaalvgar@hotmail.es>
Date: Wed, 12 Nov 2025 12:55:58 +0100
Subject: [PATCH 06/39] add test for bug found

---
 test/test_geant4.jl | 11 ++++++++++-
 1 file changed, 10 insertions(+), 1 deletion(-)

diff --git a/test/test_geant4.jl b/test/test_geant4.jl
index 23913f5ba..213bbdd62 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,14 @@ 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)
+    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)
@@ -108,4 +117,4 @@ end
     @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
+end

From 5b3db1b869af24acf090efd867ae3434cc07f73d Mon Sep 17 00:00:00 2001
From: Claudia Alvarez Garcia <claudiaalvgar@hotmail.es>
Date: Wed, 12 Nov 2025 12:56:21 +0100
Subject: [PATCH 07/39] pass N and unify functions

---
 src/ChargeCloudModels/ChargeCloudModels.jl | 22 ++++------------------
 1 file changed, 4 insertions(+), 18 deletions(-)

diff --git a/src/ChargeCloudModels/ChargeCloudModels.jl b/src/ChargeCloudModels/ChargeCloudModels.jl
index 07950fadc..310385764 100644
--- a/src/ChargeCloudModels/ChargeCloudModels.jl
+++ b/src/ChargeCloudModels/ChargeCloudModels.jl
@@ -58,14 +58,7 @@ function NBodyChargeCloud(center::CartesianPoint{<:Unitful.Quantity}, energy::Re
     return NBodyChargeCloud(internal_center, energy; kwargs...)
 end
 
-function NBodyChargeCloud(center::CylindricalPoint{<:Unitful.Quantity}, energy::RealQuantity;
-                          kwargs...)
-    internal_center = to_internal_point(center)
-    return NBodyChargeCloud(internal_center, energy; kwargs...)
-end
-
-function NBodyChargeCloud(center::CylindricalPoint{<:SSDFloat}, energy::RealQuantity;
-                          kwargs...)
+function NBodyChargeCloud(center::CylindricalPoint{T}, energy::RealQuantity; kwargs...) where T<:Union{SSDFloat, Unitful.Quantity}
     internal_center = to_internal_point(center)
     return NBodyChargeCloud(internal_center, energy; kwargs...)
 end
@@ -142,19 +135,12 @@ NBodyChargeCloud(center, energy, 200, number_of_shells = 3)
 function NBodyChargeCloud(center::CartesianPoint{<:Unitful.Quantity}, energy::RealQuantity, N::Integer;
                           kwargs...)
     internal_center = to_internal_point(center)
-    return NBodyChargeCloud(internal_center, energy; kwargs...)
-end
-
-function NBodyChargeCloud(center::CylindricalPoint{<:Unitful.Quantity}, energy::RealQuantity, N::Integer;
-                          kwargs...)
-    internal_center = to_internal_point(center)
-    return NBodyChargeCloud(internal_center, energy; kwargs...)
+    return NBodyChargeCloud(internal_center, energy, N; kwargs...)
 end
 
-function NBodyChargeCloud(center::CylindricalPoint{<:SSDFloat}, energy::RealQuantity, N::Integer;
-                          kwargs...)
+function NBodyChargeCloud(center::CylindricalPoint{T}, energy::RealQuantity, N::Integer; kwargs...) where T<:Union{SSDFloat, Unitful.Quantity}
     internal_center = to_internal_point(center)
-    return NBodyChargeCloud(internal_center, energy; kwargs...)
+    return NBodyChargeCloud(internal_center, energy, N; kwargs...)
 end
 
 function NBodyChargeCloud(center::CartesianPoint{T}, energy::RealQuantity, N::Integer, particle_type::Type{PT} = Gamma;

From 9bfe583bf8ac82a6ea4716a63a7498c0078a01a9 Mon Sep 17 00:00:00 2001
From: Claudia Alvarez Garcia <claudiaalvgar@hotmail.es>
Date: Wed, 12 Nov 2025 12:56:41 +0100
Subject: [PATCH 08/39] restrict Unitful.Quantity

---
 src/ConstructiveSolidGeometry/PointsAndVectors/Points.jl | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/ConstructiveSolidGeometry/PointsAndVectors/Points.jl b/src/ConstructiveSolidGeometry/PointsAndVectors/Points.jl
index 6cc3d2d3a..2475dbc48 100644
--- a/src/ConstructiveSolidGeometry/PointsAndVectors/Points.jl
+++ b/src/ConstructiveSolidGeometry/PointsAndVectors/Points.jl
@@ -246,9 +246,9 @@ end
 
 # Handles Unitful quantities or plain numbers 
 function CylindricalPoint(r, φ, z)
-    r_val = isa(r, Unitful.Quantity) ? float(ustrip(uconvert(u"m", r))) : float(r)
-    φ_val = isa(φ, Unitful.Quantity) ? float(ustrip(uconvert(u"rad", φ))) : float(φ)
-    z_val = isa(z, Unitful.Quantity) ? float(ustrip(uconvert(u"m", z))) : float(z)
+    r_val = isa(r, Unitful.Length) ? float(ustrip(internal_length_unit, r)) : float(r)
+    φ_val = isa(φ, Unitful.Quantity) ? float(ustrip(internal_angle_unit, φ)) : float(φ)
+    z_val = isa(z, Unitful.Length) ? float(ustrip(internal_length_unit, z)) : float(z)
 
     T = promote_type(typeof(r_val), typeof(φ_val), typeof(z_val))
     CylindricalPoint{T}(T(r_val), T(φ_val), T(z_val))

From 035398f32fecfbde9a05529f336dbbf8e127cd9e Mon Sep 17 00:00:00 2001
From: Claudia Alvarez Garcia <claudiaalvgar@hotmail.es>
Date: Wed, 12 Nov 2025 12:56:56 +0100
Subject: [PATCH 09/39] implemented suggestions

---
 src/Event/Event.jl | 92 ++++++++++++++++++++++++----------------------
 1 file changed, 49 insertions(+), 43 deletions(-)

diff --git a/src/Event/Event.jl b/src/Event/Event.jl
index e25c369ee..d658410e7 100644
--- a/src/Event/Event.jl
+++ b/src/Event/Event.jl
@@ -2,26 +2,32 @@
   to_internal_point(pt::AbstractCoordinatePoint)
         
 Converts any `AbstractCoordinatePoint` (Cartesian or Cylindrical, with or without units) to a `CartesianPoint`.
- """
+"""
+function to_internal_point(pt::CartesianPoint)
+    T = pt.x isa Unitful.Quantity ? float(typeof(ustrip(pt.x))) : typeof(pt.x)
+    x = pt.x isa Unitful.Quantity ? float(to_internal_units(pt.x)) : pt.x
+    y = pt.y isa Unitful.Quantity ? float(to_internal_units(pt.y)) : pt.y
+    z = pt.z isa Unitful.Quantity ? float(to_internal_units(pt.z)) : pt.z
+    return CartesianPoint{T}(x, y, z)
+end
+
+function to_internal_point(pt::CylindricalPoint)
+    T = pt.r isa Unitful.Quantity ? float(typeof(ustrip(pt.r))) : float(typeof(pt.r))
+    r = pt.r isa Unitful.Quantity ? float(to_internal_units(pt.r)) : float(pt.r)
+    φ = pt.φ isa Unitful.Quantity ? float(ustrip(uconvert(u"rad", pt.φ))) : float(pt.φ)
+    z = pt.z isa Unitful.Quantity ? float(to_internal_units(pt.z)) : float(pt.z)
+    return CartesianPoint{T}(r * cos(φ), r * sin(φ), z)
+end
+
 function to_internal_point(pt::AbstractCoordinatePoint)
-    if pt isa CartesianPoint
-        T = pt.x isa Unitful.Quantity ? float(typeof(ustrip(pt.x))) : typeof(pt.x)
-        x = pt.x isa Unitful.Quantity ? float(to_internal_units(pt.x)) : pt.x
-        y = pt.y isa Unitful.Quantity ? float(to_internal_units(pt.y)) : pt.y
-        z = pt.z isa Unitful.Quantity ? float(to_internal_units(pt.z)) : pt.z
-        return CartesianPoint{T}(x, y, z)
-        
-    elseif pt isa CylindricalPoint
-        T = pt.r isa Unitful.Quantity ? float(typeof(ustrip(pt.r))) : float(typeof(pt.r))
-        r = pt.r isa Unitful.Quantity ? float(to_internal_units(pt.r)) : float(pt.r)
-        φ = pt.φ isa Unitful.Quantity ? float(ustrip(uconvert(u"rad", pt.φ))) : float(pt.φ)
-        z = pt.z isa Unitful.Quantity ? float(to_internal_units(pt.z)) : float(pt.z)
-        return CartesianPoint{T}(r * cos(φ), r * sin(φ), z)
-        
-     else
-        error("Unsupported point type $(typeof(pt)). Expected CartesianPoint or CylindricalPoint.")
-     end
+    error("Unsupported point type $(typeof(pt)). Expected CartesianPoint or CylindricalPoint.")
+end
+
+function to_internal_point(p::SVector{3, <:Unitful.Quantity})
+    vals = ustrip.(uconvert.(u"m", p))   # convert to meters and strip units
+    return SVector{3, typeof(vals[1])}(vals)  # preserves type
 end
+
 """
   mutable struct Event{T <: SSDFloat}
          
@@ -43,41 +49,42 @@ 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 = 1.0) where {T <: RealQuantity}
+function Event(location::AbstractCoordinatePoint, energy::RealQuantity = 1.0)
     pt_internal = to_internal_point(location)
-    T_ = typeof(pt_internal.x)
-    evt = Event{T_}()
+    T = eltype(pt_internal)
+    evt = Event{T}()
     evt.locations = VectorOfArrays([[pt_internal]])
-    evt.energies  = VectorOfArrays([[T_(to_internal_units(energy))]])
+    evt.energies  = VectorOfArrays([[T(to_internal_units(energy))]])
     return evt
 end
 
-function Event(locations::Vector{<:AbstractCoordinatePoint{T}}, energies::Vector{<:RealQuantity} = fill(1.0, length(locations)); max_interaction_distance::Union{<:Real, <:LengthQuantity} = NaN) where {T <: RealQuantity}
-    pts_internal = [to_internal_point(pt) for pt in locations]
-    T_ = typeof(first(pts_internal).x)
-    d::T_ = T_(to_internal_units(max_interaction_distance))
+function Event(locations::Vector{<:AbstractCoordinatePoint}, energies::Vector{<:RealQuantity} = fill(1.0, length(locations)); max_interaction_distance::Union{<:Real, <:LengthQuantity} = NaN)
+    pts_internal = to_internal_point.(locations)
+    T = eltype(first(pts_internal))
+    d::T = T(to_internal_units(max_interaction_distance))
     @assert isnan(d) || d >= 0
-
-    evt = Event{T_}()
+    evt = Event{T}()
     
     if isnan(d)
         evt.locations = VectorOfArrays(broadcast(pt -> [pt], pts_internal))
-        evt.energies  = VectorOfArrays(broadcast(E -> [T_(to_internal_units(E))], energies))
+        evt.energies  = VectorOfArrays(broadcast(E -> [T(to_internal_units(E))], energies))
     else
-        evt.locations, evt.energies = group_points_by_distance(pts_internal, 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}} = [[1.0 for _ in group] for group in locations]) where {T <: RealQuantity}
+function Event(locations::Vector{<:Vector{<:AbstractCoordinatePoint}}, energies::Vector{<:Vector{<:RealQuantity}} = [[1.0 for _ in group] for group in locations])
+
     locs_internal = [[to_internal_point(pt) for pt in group] for group in locations]
     ens_internal  = [[to_internal_units(E) for E in group] for group in energies]
-    T_ = typeof(first(first(locs_internal)).x)
+
+    T = eltype(first(first(locs_internal)))
     
-    locs_T = [[convert(CartesianPoint{T_}, pt) for pt in group] for group in locs_internal]
-    ens_T  = [[T_(E) for E in group] for group in ens_internal]
+    locs_T = [[convert(CartesianPoint{T}, pt) for pt in group] for group in locs_internal]
+    ens_T  = [[T(E) for E in group] for group in ens_internal]
 
-    evt = Event{T_}()
+    evt = Event{T}()
     evt.locations = VectorOfArrays(locs_T)
     evt.energies  = VectorOfArrays(ens_T)
     return evt
@@ -135,17 +142,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 <: RealQuantity, PT <: ParticleType}
-
+    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_point.(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
 

From 652d3242b75162910e58c41e8aa6f009dc894c43 Mon Sep 17 00:00:00 2001
From: Claudia Alvarez Garcia <claudiaalvgar@hotmail.es>
Date: Wed, 12 Nov 2025 12:57:09 +0100
Subject: [PATCH 10/39] export to_internal_point function

---
 src/SolidStateDetectors.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/SolidStateDetectors.jl b/src/SolidStateDetectors.jl
index 4841a2a81..39df2a6cc 100644
--- a/src/SolidStateDetectors.jl
+++ b/src/SolidStateDetectors.jl
@@ -81,7 +81,7 @@ export calculate_stored_energy, calculate_mutual_capacitance, calculate_capacita
 export simulate_waveforms
 export run_geant4_simulation
 export Simulation, simulate!
-export Event, drift_charges!
+export Event, drift_charges!, to_internal_point
 export add_baseline_and_extend_tail
 export NBodyChargeCloud
 

From 2f0981dc0a47d081f71592fb09011d5f14c96ea2 Mon Sep 17 00:00:00 2001
From: Claudia Alvarez Garcia <claudiaalvgar@hotmail.es>
Date: Wed, 12 Nov 2025 16:01:13 +0100
Subject: [PATCH 11/39] Add tests for new NBodyChargeCloud functions

---
 test/test_charge_drift_models.jl | 45 ++++++++++++++++++++++++++++++++
 1 file changed, 45 insertions(+)

diff --git a/test/test_charge_drift_models.jl b/test/test_charge_drift_models.jl
index 385197838..bf7f54383 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", π/4u"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.(uconvert.(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.(uconvert.(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.(uconvert.(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.(uconvert.(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.(uconvert.(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.(uconvert.(u"eV", Edep)))
+    @test all(x -> all(isfinite, (x.x, x.y, x.z)), nb_zero_cyl.locations)
+
+end

From 67343ad917e7d8966c4fc33e7e21bf6a07c156d9 Mon Sep 17 00:00:00 2001
From: Claudia Alvarez Garcia <claudiaalvgar@hotmail.es>
Date: Thu, 13 Nov 2025 16:48:44 +0100
Subject: [PATCH 12/39] comment erro line

---
 ext/SolidStateDetectorsGeant4Ext.jl | 9 +++------
 1 file changed, 3 insertions(+), 6 deletions(-)

diff --git a/ext/SolidStateDetectorsGeant4Ext.jl b/ext/SolidStateDetectorsGeant4Ext.jl
index b07a4ea7b..458a4c185 100644
--- a/ext/SolidStateDetectorsGeant4Ext.jl
+++ b/ext/SolidStateDetectorsGeant4Ext.jl
@@ -220,12 +220,9 @@ function SolidStateDetectors.run_geant4_simulation(app::G4JLApplication, number_
         evtno += 1
     end
 
-    evts_with_points = map(evts) do evt
-        merge(evt, (pos = SolidStateDetectors.to_internal_point(evt.pos),))
-    end
-
-    tbl = Table(evts_with_points)
-    return RadiationDetectorSignals.group_by_evtno(tbl)
+    #evts_with_points = map(evt -> merge(evt, (pos = cartesian_zero + evt.pos,)), evts) #This line triggers the error when using `writedata` function
+    
+    return RadiationDetectorSignals.group_by_evtno(Table(evts))
 end 
 
 

From 6d6c10372c0678b7f6affd4dfacb69519b0e8b76 Mon Sep 17 00:00:00 2001
From: Claudia Alvarez Garcia <claudiaalvgar@hotmail.es>
Date: Thu, 13 Nov 2025 16:48:58 +0100
Subject: [PATCH 13/39] add test

---
 test/test_geant4.jl | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/test/test_geant4.jl b/test/test_geant4.jl
index 213bbdd62..62a247e79 100644
--- a/test/test_geant4.jl
+++ b/test/test_geant4.jl
@@ -89,6 +89,10 @@ end
             LegendHDF5IO.writedata(h5, "geant4", evts)
         end
         @test isfile(tmpfile)
+        
+        evts == LegendHDF5IO.lh5open(tmpfile) do h5
+            LegendHDF5IO.readdata(h5, "geant4")
+        end
     end
     
     # Generate waveforms

From 1df5e980f800ff4738a16b2bf4d77085a110cb95 Mon Sep 17 00:00:00 2001
From: Claudia Alvarez Garcia <claudiaalvgar@hotmail.es>
Date: Thu, 13 Nov 2025 16:49:18 +0100
Subject: [PATCH 14/39] unify functions, just 1 function failed

---
 src/ChargeCloudModels/ChargeCloudModels.jl | 22 ++++++----------------
 1 file changed, 6 insertions(+), 16 deletions(-)

diff --git a/src/ChargeCloudModels/ChargeCloudModels.jl b/src/ChargeCloudModels/ChargeCloudModels.jl
index 310385764..7f555fb66 100644
--- a/src/ChargeCloudModels/ChargeCloudModels.jl
+++ b/src/ChargeCloudModels/ChargeCloudModels.jl
@@ -52,14 +52,9 @@ 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::CartesianPoint{<:Unitful.Quantity}, energy::RealQuantity;
-                          kwargs...)
-    internal_center = to_internal_point(center)
-    return NBodyChargeCloud(internal_center, energy; kwargs...)
-end
-
-function NBodyChargeCloud(center::CylindricalPoint{T}, energy::RealQuantity; kwargs...) where T<:Union{SSDFloat, Unitful.Quantity}
-    internal_center = to_internal_point(center)
+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
 
@@ -132,14 +127,9 @@ 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::CartesianPoint{<:Unitful.Quantity}, energy::RealQuantity, N::Integer;
-                          kwargs...)
-    internal_center = to_internal_point(center)
-    return NBodyChargeCloud(internal_center, energy, N; kwargs...)
-end
-
-function NBodyChargeCloud(center::CylindricalPoint{T}, energy::RealQuantity, N::Integer; kwargs...) where T<:Union{SSDFloat, Unitful.Quantity}
-    internal_center = to_internal_point(center)
+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
 

From c6a34e25086e907dc2cd2221e07322a36ef993f1 Mon Sep 17 00:00:00 2001
From: Claudia Alvarez Garcia <claudiaalvgar@hotmail.es>
Date: Thu, 13 Nov 2025 16:49:38 +0100
Subject: [PATCH 15/39] add functions to to_internal_units

---
 .../PointsAndVectors/Points.jl                | 37 ++++++++++++++-----
 1 file changed, 27 insertions(+), 10 deletions(-)

diff --git a/src/ConstructiveSolidGeometry/PointsAndVectors/Points.jl b/src/ConstructiveSolidGeometry/PointsAndVectors/Points.jl
index 2475dbc48..dc29e893c 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)
@@ -246,9 +240,9 @@ end
 
 # Handles Unitful quantities or plain numbers 
 function CylindricalPoint(r, φ, z)
-    r_val = isa(r, Unitful.Length) ? float(ustrip(internal_length_unit, r)) : float(r)
-    φ_val = isa(φ, Unitful.Quantity) ? float(ustrip(internal_angle_unit, φ)) : float(φ)
-    z_val = isa(z, Unitful.Length) ? float(ustrip(internal_length_unit, z)) : float(z)
+    r_val = isa(r, Unitful.Length) ? float(to_internal_units(r)) : float(r)
+    φ_val = isa(φ, Unitful.Quantity) ? float(to_internal_units(φ)) : float(φ)
+    z_val = isa(z, Unitful.Length) ? float(to_internal_units(z)) : float(z)
 
     T = promote_type(typeof(r_val), typeof(φ_val), typeof(z_val))
     CylindricalPoint{T}(T(r_val), T(φ_val), T(z_val))
@@ -314,4 +308,27 @@ 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
+
+function to_internal_units(p::SVector{3, <:Unitful.Quantity})
+    vals = ustrip.(internal_length_unit, p)
+    return SVector{3, typeof(vals[1])}(vals)
+end
+

From 09a64e3c280f68d61be81dbf7e2ffe53711438be Mon Sep 17 00:00:00 2001
From: Claudia Alvarez Garcia <claudiaalvgar@hotmail.es>
Date: Thu, 13 Nov 2025 16:49:52 +0100
Subject: [PATCH 16/39] added suggested changes

---
 src/Event/Event.jl | 60 ++++++++++------------------------------------
 1 file changed, 13 insertions(+), 47 deletions(-)

diff --git a/src/Event/Event.jl b/src/Event/Event.jl
index d658410e7..182e01381 100644
--- a/src/Event/Event.jl
+++ b/src/Event/Event.jl
@@ -1,33 +1,3 @@
-"""
-  to_internal_point(pt::AbstractCoordinatePoint)
-        
-Converts any `AbstractCoordinatePoint` (Cartesian or Cylindrical, with or without units) to a `CartesianPoint`.
-"""
-function to_internal_point(pt::CartesianPoint)
-    T = pt.x isa Unitful.Quantity ? float(typeof(ustrip(pt.x))) : typeof(pt.x)
-    x = pt.x isa Unitful.Quantity ? float(to_internal_units(pt.x)) : pt.x
-    y = pt.y isa Unitful.Quantity ? float(to_internal_units(pt.y)) : pt.y
-    z = pt.z isa Unitful.Quantity ? float(to_internal_units(pt.z)) : pt.z
-    return CartesianPoint{T}(x, y, z)
-end
-
-function to_internal_point(pt::CylindricalPoint)
-    T = pt.r isa Unitful.Quantity ? float(typeof(ustrip(pt.r))) : float(typeof(pt.r))
-    r = pt.r isa Unitful.Quantity ? float(to_internal_units(pt.r)) : float(pt.r)
-    φ = pt.φ isa Unitful.Quantity ? float(ustrip(uconvert(u"rad", pt.φ))) : float(pt.φ)
-    z = pt.z isa Unitful.Quantity ? float(to_internal_units(pt.z)) : float(pt.z)
-    return CartesianPoint{T}(r * cos(φ), r * sin(φ), z)
-end
-
-function to_internal_point(pt::AbstractCoordinatePoint)
-    error("Unsupported point type $(typeof(pt)). Expected CartesianPoint or CylindricalPoint.")
-end
-
-function to_internal_point(p::SVector{3, <:Unitful.Quantity})
-    vals = ustrip.(uconvert.(u"m", p))   # convert to meters and strip units
-    return SVector{3, typeof(vals[1])}(vals)  # preserves type
-end
-
 """
   mutable struct Event{T <: SSDFloat}
          
@@ -49,8 +19,8 @@ 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, energy::RealQuantity = 1.0)
-    pt_internal = to_internal_point(location)
+function Event(location::AbstractCoordinatePoint, energy::RealQuantity = 1)
+    pt_internal = to_internal_units(location)
     T = eltype(pt_internal)
     evt = Event{T}()
     evt.locations = VectorOfArrays([[pt_internal]])
@@ -58,14 +28,14 @@ function Event(location::AbstractCoordinatePoint, energy::RealQuantity = 1.0)
     return evt
 end
 
-function Event(locations::Vector{<:AbstractCoordinatePoint}, energies::Vector{<:RealQuantity} = fill(1.0, length(locations)); max_interaction_distance::Union{<:Real, <:LengthQuantity} = NaN)
-    pts_internal = to_internal_point.(locations)
+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
+    @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)
+    if isnan(d) # default: no grouping, the charges from different hits drift independently
         evt.locations = VectorOfArrays(broadcast(pt -> [pt], pts_internal))
         evt.energies  = VectorOfArrays(broadcast(E -> [T(to_internal_units(E))], energies))
     else
@@ -74,19 +44,15 @@ function Event(locations::Vector{<:AbstractCoordinatePoint}, energies::Vector{<:
     return evt
 end
 
-function Event(locations::Vector{<:Vector{<:AbstractCoordinatePoint}}, energies::Vector{<:Vector{<:RealQuantity}} = [[1.0 for _ in group] for group in locations])
-
-    locs_internal = [[to_internal_point(pt) for pt in group] for group in locations]
-    ens_internal  = [[to_internal_units(E) for E in group] for group in energies]
+function Event(locations::Vector{<:Vector{<:AbstractCoordinatePoint}}, energies::Vector{<:Vector{<:RealQuantity}} = [[1 for _ in group] for group in locations])
 
-    T = eltype(first(first(locs_internal)))
+    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]
     
-    locs_T = [[convert(CartesianPoint{T}, pt) for pt in group] for group in locs_internal]
-    ens_T  = [[T(E) for E in group] for group in ens_internal]
-
     evt = Event{T}()
-    evt.locations = VectorOfArrays(locs_T)
-    evt.energies  = VectorOfArrays(ens_T)
+    evt.locations = VectorOfArrays(locs)
+    evt.energies  = VectorOfArrays(ens)
     return evt
 end
 
@@ -147,7 +113,7 @@ function Event(
     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_point.(loc) for loc in locations]
+    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))

From 8d6268d35dbb0a5df69464375381d1c508052ba1 Mon Sep 17 00:00:00 2001
From: Claudia Alvarez Garcia <claudiaalvgar@hotmail.es>
Date: Thu, 13 Nov 2025 16:50:10 +0100
Subject: [PATCH 17/39] delete export to_internal_point

---
 src/SolidStateDetectors.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/SolidStateDetectors.jl b/src/SolidStateDetectors.jl
index 39df2a6cc..4841a2a81 100644
--- a/src/SolidStateDetectors.jl
+++ b/src/SolidStateDetectors.jl
@@ -81,7 +81,7 @@ export calculate_stored_energy, calculate_mutual_capacitance, calculate_capacita
 export simulate_waveforms
 export run_geant4_simulation
 export Simulation, simulate!
-export Event, drift_charges!, to_internal_point
+export Event, drift_charges!
 export add_baseline_and_extend_tail
 export NBodyChargeCloud
 

From a266d8d216fa425b1debdebce8e72c809b0cea58 Mon Sep 17 00:00:00 2001
From: Claudia Alvarez Garcia <claudiaalvgar@hotmail.es>
Date: Thu, 13 Nov 2025 16:50:58 +0100
Subject: [PATCH 18/39] delete ucovert

---
 test/test_charge_drift_models.jl | 15 ++++++++-------
 1 file changed, 8 insertions(+), 7 deletions(-)

diff --git a/test/test_charge_drift_models.jl b/test/test_charge_drift_models.jl
index bf7f54383..c54d98af5 100644
--- a/test/test_charge_drift_models.jl
+++ b/test/test_charge_drift_models.jl
@@ -570,43 +570,44 @@ end
 
     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", π/4u"rad", 0.05u"m")
+    #cylindrical_unit = CylindricalPoint(0.01u"m", π/4u"rad", 0.05u"m")
+    cylindrical_unit = CylindricalPoint(0.01, π/4, 0.05)
     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.(uconvert.(u"eV", Edep)))
+    @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.(uconvert.(u"eV", Edep)))
+    @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.(uconvert.(u"eV", Edep)))
+    @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.(uconvert.(u"eV", Edep)))
+    @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.(uconvert.(u"eV", Edep)))
+    @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.(uconvert.(u"eV", Edep)))
+    @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

From ebe8db6a60bbf88bf224cdd331b06b1bbd334f52 Mon Sep 17 00:00:00 2001
From: Felix Hagemann <hagemann@berkeley.edu>
Date: Thu, 13 Nov 2025 17:37:24 -0800
Subject: [PATCH 19/39] Add `LegendHDF5IO` support for
 `Cartesian/CylindricalPoint`

---
 ext/SolidStateDetectorsLegendHDF5IOExt.jl | 95 +++++++++++++++++++++++
 test/ConstructiveSolidGeometry/CSG_IO.jl  | 65 +++++++++++++++-
 2 files changed, 159 insertions(+), 1 deletion(-)

diff --git a/ext/SolidStateDetectorsLegendHDF5IOExt.jl b/ext/SolidStateDetectorsLegendHDF5IOExt.jl
index 7147e71ee..1bb23f377 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
@@ -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
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

From a6a74a10300a8716cfb269f572103cd85c00a8d3 Mon Sep 17 00:00:00 2001
From: Felix Hagemann <hagemann@berkeley.edu>
Date: Thu, 13 Nov 2025 18:18:32 -0800
Subject: [PATCH 20/39] Revert "comment erro line"

This reverts commit 178caea4be60ec96cf394c349e7eedd32860bb66.
---
 ext/SolidStateDetectorsGeant4Ext.jl | 9 ++++++---
 1 file changed, 6 insertions(+), 3 deletions(-)

diff --git a/ext/SolidStateDetectorsGeant4Ext.jl b/ext/SolidStateDetectorsGeant4Ext.jl
index 458a4c185..b07a4ea7b 100644
--- a/ext/SolidStateDetectorsGeant4Ext.jl
+++ b/ext/SolidStateDetectorsGeant4Ext.jl
@@ -220,9 +220,12 @@ function SolidStateDetectors.run_geant4_simulation(app::G4JLApplication, number_
         evtno += 1
     end
 
-    #evts_with_points = map(evt -> merge(evt, (pos = cartesian_zero + evt.pos,)), evts) #This line triggers the error when using `writedata` function
-    
-    return RadiationDetectorSignals.group_by_evtno(Table(evts))
+    evts_with_points = map(evts) do evt
+        merge(evt, (pos = SolidStateDetectors.to_internal_point(evt.pos),))
+    end
+
+    tbl = Table(evts_with_points)
+    return RadiationDetectorSignals.group_by_evtno(tbl)
 end 
 
 

From 023d2cf5bf273b6fab930df4e0774a159ebbd2a8 Mon Sep 17 00:00:00 2001
From: Felix Hagemann <hagemann@berkeley.edu>
Date: Thu, 13 Nov 2025 18:23:36 -0800
Subject: [PATCH 21/39] Revert "convert CartesianPoint to SVector to be
 returned from run_geant4_simulation fixes bug"

This reverts commit cc73fea1d44f2a588029b4a8acb22faa21d11443.
---
 ext/SolidStateDetectorsGeant4Ext.jl | 7 ++-----
 1 file changed, 2 insertions(+), 5 deletions(-)

diff --git a/ext/SolidStateDetectorsGeant4Ext.jl b/ext/SolidStateDetectorsGeant4Ext.jl
index b07a4ea7b..feabfb2c1 100644
--- a/ext/SolidStateDetectorsGeant4Ext.jl
+++ b/ext/SolidStateDetectorsGeant4Ext.jl
@@ -220,12 +220,9 @@ function SolidStateDetectors.run_geant4_simulation(app::G4JLApplication, number_
         evtno += 1
     end
 
-    evts_with_points = map(evts) do evt
-        merge(evt, (pos = SolidStateDetectors.to_internal_point(evt.pos),))
-    end
+    evts_with_points = map(evt -> merge(evt, (pos = cartesian_zero + evt.pos,)), evts)
 
-    tbl = Table(evts_with_points)
-    return RadiationDetectorSignals.group_by_evtno(tbl)
+    return RadiationDetectorSignals.group_by_evtno(Table(evts_with_points))
 end 
 
 

From a4dd5febfb89fd7316e1a48091bd8cb1a1bb9e26 Mon Sep 17 00:00:00 2001
From: Claudia Alvarez Garcia <claudiaalvgar@hotmail.es>
Date: Fri, 14 Nov 2025 11:05:35 +0100
Subject: [PATCH 22/39] put back Cylindrical test with units

---
 test/test_charge_drift_models.jl | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/test/test_charge_drift_models.jl b/test/test_charge_drift_models.jl
index c54d98af5..c234af0ea 100644
--- a/test/test_charge_drift_models.jl
+++ b/test/test_charge_drift_models.jl
@@ -570,8 +570,7 @@ end
 
     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", π/4u"rad", 0.05u"m")
-    cylindrical_unit = CylindricalPoint(0.01, π/4, 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"
     

From f000ab8b1e8088cf1ff6dfe1e0f2fa743949566d Mon Sep 17 00:00:00 2001
From: Claudia Alvarez Garcia <claudiaalvgar@hotmail.es>
Date: Fri, 14 Nov 2025 11:06:06 +0100
Subject: [PATCH 23/39] change CylindricalPoint function and delete function
 to_internal_units(p::SVector{3, <:Unitful.Quantity})

---
 .../PointsAndVectors/Points.jl                | 37 ++++++++++++-------
 1 file changed, 24 insertions(+), 13 deletions(-)

diff --git a/src/ConstructiveSolidGeometry/PointsAndVectors/Points.jl b/src/ConstructiveSolidGeometry/PointsAndVectors/Points.jl
index dc29e893c..884f95a6f 100644
--- a/src/ConstructiveSolidGeometry/PointsAndVectors/Points.jl
+++ b/src/ConstructiveSolidGeometry/PointsAndVectors/Points.jl
@@ -238,17 +238,34 @@ struct CylindricalPoint{T} <: AbstractCylindricalPoint{T}
     CylindricalPoint{T}(r::Real, φ::Real, z::Real) where {T} = new(T(r), mod(T(φ),T(2π)), T(z))
 end
 
-# Handles Unitful quantities or plain numbers 
 function CylindricalPoint(r, φ, z)
-    r_val = isa(r, Unitful.Length) ? float(to_internal_units(r)) : float(r)
-    φ_val = isa(φ, Unitful.Quantity) ? float(to_internal_units(φ)) : float(φ)
-    z_val = isa(z, Unitful.Length) ? float(to_internal_units(z)) : float(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
 
-    T = promote_type(typeof(r_val), typeof(φ_val), typeof(z_val))
-    CylindricalPoint{T}(T(r_val), T(φ_val), T(z_val))
+    if !(φ isa Real || φ isa Unitful.Quantity)
+        throw(ArgumentError("Expected `φ` to be an angle or Real, got unit $(Unitful.unit(φ))"))
+    end
 
-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)
 
+    vals = (r_val, φ_val, z_val)
+    has_int   = any(v -> v isa Integer, vals)
+    has_float = any(v -> v isa AbstractFloat, vals)
+
+    if has_int && has_float
+        throw(ArgumentError("Mixed integer and float inputs, use all integers or all floats."))
+    end
+    
+    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))
@@ -326,9 +343,3 @@ end
 function to_internal_units(pt::AbstractCoordinatePoint)
     error("Unsupported point type $(typeof(pt)). Expected CartesianPoint or CylindricalPoint.")
 end
-
-function to_internal_units(p::SVector{3, <:Unitful.Quantity})
-    vals = ustrip.(internal_length_unit, p)
-    return SVector{3, typeof(vals[1])}(vals)
-end
-

From 5fe64182fe1a3c68631ad9da7a20a111571ed324 Mon Sep 17 00:00:00 2001
From: Felix Hagemann <hagemann@berkeley.edu>
Date: Sat, 15 Nov 2025 10:47:59 -0800
Subject: [PATCH 24/39] Update stripping units of `CylindricalPoint`

---
 .../PointsAndVectors/Points.jl                         | 10 +---------
 1 file changed, 1 insertion(+), 9 deletions(-)

diff --git a/src/ConstructiveSolidGeometry/PointsAndVectors/Points.jl b/src/ConstructiveSolidGeometry/PointsAndVectors/Points.jl
index 884f95a6f..34799df3c 100644
--- a/src/ConstructiveSolidGeometry/PointsAndVectors/Points.jl
+++ b/src/ConstructiveSolidGeometry/PointsAndVectors/Points.jl
@@ -243,7 +243,7 @@ function CylindricalPoint(r, φ, z)
         throw(ArgumentError("Expected `r` to be a length or Real, got unit $(Unitful.unit(r))"))
     end
 
-    if !(φ isa Real || φ isa Unitful.Quantity)
+    if !(φ isa Real || φ isa Unitful.Quantity{<:Real, NoDims})
         throw(ArgumentError("Expected `φ` to be an angle or Real, got unit $(Unitful.unit(φ))"))
     end
 
@@ -255,14 +255,6 @@ function CylindricalPoint(r, φ, z)
     φ_val = to_internal_units(φ)
     z_val = to_internal_units(z)
 
-    vals = (r_val, φ_val, z_val)
-    has_int   = any(v -> v isa Integer, vals)
-    has_float = any(v -> v isa AbstractFloat, vals)
-
-    if has_int && has_float
-        throw(ArgumentError("Mixed integer and float inputs, use all integers or all floats."))
-    end
-    
     return CylindricalPoint(r_val, φ_val, z_val)
 end
                                                      

From 253ec305b063484b31c4e13dd0a4697320842abc Mon Sep 17 00:00:00 2001
From: Felix Hagemann <hagemann@berkeley.edu>
Date: Sat, 15 Nov 2025 10:48:23 -0800
Subject: [PATCH 25/39] Add tests for units/type promotion of
 `CylindricalPoint`

---
 .../CSG_coordinates.jl                        | 47 +++++++++++++++++--
 1 file changed, 44 insertions(+), 3 deletions(-)

diff --git a/test/ConstructiveSolidGeometry/CSG_coordinates.jl b/test/ConstructiveSolidGeometry/CSG_coordinates.jl
index 38078c7df..7d6bbec8a 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,35 @@ 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.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 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 +101,25 @@ 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.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 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

From 70b277bf310e7367e1ede8bb16a209050f919dd5 Mon Sep 17 00:00:00 2001
From: Claudia Alvarez Garcia <claudiaalvgar@hotmail.es>
Date: Tue, 25 Nov 2025 10:45:15 +0100
Subject: [PATCH 26/39] convert cluster_radius to internal units

---
 src/ChargeClustering/ChargeClustering.jl | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/src/ChargeClustering/ChargeClustering.jl b/src/ChargeClustering/ChargeClustering.jl
index 8b03f78e7..17bc6437d 100644
--- a/src/ChargeClustering/ChargeClustering.jl
+++ b/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)

From 4416ac539b1a5a85dc1a3b5649e7f5a621c953af Mon Sep 17 00:00:00 2001
From: Claudia Alvarez Garcia <claudiaalvgar@hotmail.es>
Date: Tue, 25 Nov 2025 10:45:45 +0100
Subject: [PATCH 27/39] Add CartesianPoint constructor to avoid units

---
 .../PointsAndVectors/Points.jl                | 22 +++++++++++++++++--
 1 file changed, 20 insertions(+), 2 deletions(-)

diff --git a/src/ConstructiveSolidGeometry/PointsAndVectors/Points.jl b/src/ConstructiveSolidGeometry/PointsAndVectors/Points.jl
index 34799df3c..a421b8337 100644
--- a/src/ConstructiveSolidGeometry/PointsAndVectors/Points.jl
+++ b/src/ConstructiveSolidGeometry/PointsAndVectors/Points.jl
@@ -66,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))
@@ -198,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)}()
@@ -238,6 +255,7 @@ 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))"))
@@ -257,7 +275,7 @@ function CylindricalPoint(r, φ, 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))

From 8e25960eb42e7cd92d92623f5a7f7c8a0b72ff42 Mon Sep 17 00:00:00 2001
From: Claudia Alvarez Garcia <claudiaalvgar@hotmail.es>
Date: Tue, 25 Nov 2025 10:46:02 +0100
Subject: [PATCH 28/39] uncomment tests

---
 test/ConstructiveSolidGeometry/CSG_coordinates.jl | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/test/ConstructiveSolidGeometry/CSG_coordinates.jl b/test/ConstructiveSolidGeometry/CSG_coordinates.jl
index 7d6bbec8a..2e270ec40 100644
--- a/test/ConstructiveSolidGeometry/CSG_coordinates.jl
+++ b/test/ConstructiveSolidGeometry/CSG_coordinates.jl
@@ -43,7 +43,7 @@ import Unitful: @u_str
         @test @inferred(inverse(f)(f(a))) ≈ a
 
 
-        @test @inferred(CartesianZero{Float32}() * u"mm") === CartesianZero{typeof(zero(Float32) * 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))
@@ -51,7 +51,6 @@ import Unitful: @u_str
         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}
@@ -63,6 +62,7 @@ import Unitful: @u_str
         @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}
@@ -71,7 +71,6 @@ import Unitful: @u_str
         @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
@@ -113,6 +112,7 @@ import Unitful: @u_str
         @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}

From ddaa78c7ae8ed0a293aed9d5f183e68537eec658 Mon Sep 17 00:00:00 2001
From: Claudia Alvarez Garcia <claudiaalvgar@hotmail.es>
Date: Tue, 25 Nov 2025 10:46:15 +0100
Subject: [PATCH 29/39] remove units from test

---
 test/test_io.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/test/test_io.jl b/test/test_io.jl
index 939851249..51d677909 100644
--- a/test/test_io.jl
+++ b/test/test_io.jl
@@ -69,7 +69,7 @@ end
         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"]])
+        pos = VectorOfVectors([[SVector{3, T}.(0.01, 0.01, 0.01)]])
     )
 
     contact_charge_signals = timed_simulate_waveforms(      
@@ -111,4 +111,4 @@ end
         end
         rm("generated_waveforms_evts_1-1.h5")
     end
-end
\ No newline at end of file
+end

From 366581b9efd948db6bba3e7d53397eff7f9d088b Mon Sep 17 00:00:00 2001
From: Claudia Alvarez Garcia <claudiaalvgar@hotmail.es>
Date: Thu, 27 Nov 2025 17:35:55 +0100
Subject: [PATCH 30/39] add back units to pos

---
 test/test_io.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/test_io.jl b/test/test_io.jl
index 51d677909..a1105b78c 100644
--- a/test/test_io.jl
+++ b/test/test_io.jl
@@ -69,7 +69,7 @@ end
         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)]])
+        pos = VectorOfVectors([[SVector{3, T}.(0.01, 0.01, 0.01) * u"m"]])
     )
 
     contact_charge_signals = timed_simulate_waveforms(      

From 1f2eb3450bb25fc2dfc1804db8704fbc2de95b77 Mon Sep 17 00:00:00 2001
From: Claudia Alvarez Garcia <claudiaalvgar@hotmail.es>
Date: Thu, 27 Nov 2025 17:36:14 +0100
Subject: [PATCH 31/39] add tests

---
 test/ConstructiveSolidGeometry/CSG_coordinates.jl | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/test/ConstructiveSolidGeometry/CSG_coordinates.jl b/test/ConstructiveSolidGeometry/CSG_coordinates.jl
index 2e270ec40..9a155eac4 100644
--- a/test/ConstructiveSolidGeometry/CSG_coordinates.jl
+++ b/test/ConstructiveSolidGeometry/CSG_coordinates.jl
@@ -66,6 +66,8 @@ import Unitful: @u_str
         @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")
@@ -116,6 +118,8 @@ import Unitful: @u_str
         @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")

From f1e85c62a51d2903b61fbb7329ef996e59ca8d02 Mon Sep 17 00:00:00 2001
From: Claudia Alvarez Garcia <claudiaalvgar@hotmail.es>
Date: Thu, 27 Nov 2025 17:36:47 +0100
Subject: [PATCH 32/39] remove units

---
 src/MCEventsProcessing/MCEventsProcessing.jl | 12 +++++++++++-
 1 file changed, 11 insertions(+), 1 deletion(-)

diff --git a/src/MCEventsProcessing/MCEventsProcessing.jl b/src/MCEventsProcessing/MCEventsProcessing.jl
index 6e309599a..1108aab92 100644
--- a/src/MCEventsProcessing/MCEventsProcessing.jl
+++ b/src/MCEventsProcessing/MCEventsProcessing.jl
@@ -63,6 +63,13 @@ simulate_waveforms(mcevents, sim, "output_dir", "my_basename", Δt = 1u"ns", ver
 !!! note
     Using values with units for `Δt` requires the package [Unitful.jl](https://github.com/JuliaPhysics/Unitful.jl).
 """
+function _get_unitless_positions_svec(point_vectors::AbstractVector{<:SVector{3, <:Quantity}})
+    converted_type = typeof(to_internal_units(point_vectors[1][1]))
+    [SVector{3, converted_type}(to_internal_units.(v)) for v in point_vectors]
+end
+_get_unitless_positions_svec(point_vectors::AbstractVector{<:SVector{3, <:Real}}) = point_vectors
+_get_unitless_positions_svec(point_vectors::AbstractVector{<:CartesianPoint{<:Real}}) = point_vectors
+
 function simulate_waveforms( mcevents_table::AbstractVector{<:NamedTuple}, sim::Simulation{T};
                              Δt::RealQuantity = 4u"ns",
                              max_nsteps::Int = 1000,
@@ -76,6 +83,10 @@ function simulate_waveforms( mcevents_table::AbstractVector{<:NamedTuple}, sim::
                              geometry_check::Bool = false,
                              verbose::Bool = false ) where {T <: SSDFloat}
     mcevents = TypedTables.Table(mcevents_table)
+    #ToDo: Provisional until we support units in CartesianPoints
+    unitless_pos = [collect(_get_unitless_positions_svec(v)) for v in mcevents.pos]
+    mcevents = TypedTables.Table((evtno = mcevents.evtno, detno = mcevents.detno, thit  = mcevents.thit, edep  = mcevents.edep, pos   = unitless_pos))
+
     n_total_physics_events = length(mcevents)
     Δtime = T(to_internal_units(Δt)) 
     n_contacts = length(sim.detector.contacts)
@@ -116,7 +127,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)

From a700a075555591339a107e08dbea1d62772336de Mon Sep 17 00:00:00 2001
From: Claudia Alvarez Garcia <claudiaalvgar@hotmail.es>
Date: Thu, 27 Nov 2025 17:37:00 +0100
Subject: [PATCH 33/39] account for suggestions

---
 .../PointsAndVectors/Points.jl                | 38 +++++++++++--------
 1 file changed, 23 insertions(+), 15 deletions(-)

diff --git a/src/ConstructiveSolidGeometry/PointsAndVectors/Points.jl b/src/ConstructiveSolidGeometry/PointsAndVectors/Points.jl
index a421b8337..86f26af77 100644
--- a/src/ConstructiveSolidGeometry/PointsAndVectors/Points.jl
+++ b/src/ConstructiveSolidGeometry/PointsAndVectors/Points.jl
@@ -58,30 +58,37 @@ struct CartesianPoint{T} <: AbstractCartesianPoint{T}
     z::T
 end
 
-#Type promotion happens here
-function CartesianPoint(x::TX, y::TY, z::TZ) where {TX<:Real,TY<:Real,TZ<:Real}
-    # ToDo: Simplify this:
-    eltypes = _csg_get_promoted_eltype.((TX,TY,TZ))
-    T = float(promote_type(eltypes...))
-    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))"
+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
-    CartesianPoint(x::Unitful.Length, y::Unitful.Length, z::Unitful.Length)
+    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
 
 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
+    vals = to_internal_units.((x, y, z))
+    tys = map(typeof, vals)
+    T = any(T -> T === Float64, tys) ? Float64 :
+        any(T -> T === Float32, tys) ? Float32 :
+        Float64
+
     CartesianPoint{T}(T.(vals)...)
 end
 
+#Type promotion happens here
+function CartesianPoint(x::TX, y::TY, z::TZ) where {TX<:Real,TY<:Real,TZ<:Real}
+    # ToDo: Simplify this:
+    eltypes = _csg_get_promoted_eltype.((TX,TY,TZ))
+    T = float(promote_type(eltypes...))
+    CartesianPoint{T}(T(x),T(y),T(z))
+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))
@@ -214,6 +221,7 @@ end
 @inline Base.zero(::CartesianZero{T}) where {T} = CartesianZero{T}()
 @inline Base.iszero(::CartesianZero) = true
 
+# 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
 

From 49c900c882f8608ab364f9e143906ba898d1a95a Mon Sep 17 00:00:00 2001
From: Claudia Alvarez Garcia <claudiaalvgar@hotmail.es>
Date: Fri, 28 Nov 2025 10:54:56 +0100
Subject: [PATCH 34/39] Suggested changes and fix doc without units

---
 docs/src/tutorials/geant4_ssd_lit.jl          |  6 +++--
 .../PointsAndVectors/Points.jl                | 22 ++++++-------------
 src/MCEventsProcessing/MCEventsProcessing.jl  |  9 +-------
 3 files changed, 12 insertions(+), 25 deletions(-)

diff --git a/docs/src/tutorials/geant4_ssd_lit.jl b/docs/src/tutorials/geant4_ssd_lit.jl
index 57f2cc54e..d47375651 100644
--- a/docs/src/tutorials/geant4_ssd_lit.jl
+++ b/docs/src/tutorials/geant4_ssd_lit.jl
@@ -89,7 +89,9 @@ 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 = "")
+#ToDo: Put this back when units are supported
+#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 +152,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/src/ConstructiveSolidGeometry/PointsAndVectors/Points.jl b/src/ConstructiveSolidGeometry/PointsAndVectors/Points.jl
index 86f26af77..153a073aa 100644
--- a/src/ConstructiveSolidGeometry/PointsAndVectors/Points.jl
+++ b/src/ConstructiveSolidGeometry/PointsAndVectors/Points.jl
@@ -56,29 +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
 
-#Unit support
+#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))"
-        ))
+    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
 
-function CartesianPoint(x::Unitful.Length, y::Unitful.Length, z::Unitful.Length)
-    vals = to_internal_units.((x, y, z))
-    tys = map(typeof, vals)
-    T = any(T -> T === Float64, tys) ? Float64 :
-        any(T -> T === Float32, tys) ? Float32 :
-        Float64
-
-    CartesianPoint{T}(T.(vals)...)
+    return CartesianPoint(x_val, y_val, z_val)
 end
 
 #Type promotion happens here
diff --git a/src/MCEventsProcessing/MCEventsProcessing.jl b/src/MCEventsProcessing/MCEventsProcessing.jl
index 1108aab92..8e15fd8c3 100644
--- a/src/MCEventsProcessing/MCEventsProcessing.jl
+++ b/src/MCEventsProcessing/MCEventsProcessing.jl
@@ -63,13 +63,6 @@ simulate_waveforms(mcevents, sim, "output_dir", "my_basename", Δt = 1u"ns", ver
 !!! note
     Using values with units for `Δt` requires the package [Unitful.jl](https://github.com/JuliaPhysics/Unitful.jl).
 """
-function _get_unitless_positions_svec(point_vectors::AbstractVector{<:SVector{3, <:Quantity}})
-    converted_type = typeof(to_internal_units(point_vectors[1][1]))
-    [SVector{3, converted_type}(to_internal_units.(v)) for v in point_vectors]
-end
-_get_unitless_positions_svec(point_vectors::AbstractVector{<:SVector{3, <:Real}}) = point_vectors
-_get_unitless_positions_svec(point_vectors::AbstractVector{<:CartesianPoint{<:Real}}) = point_vectors
-
 function simulate_waveforms( mcevents_table::AbstractVector{<:NamedTuple}, sim::Simulation{T};
                              Δt::RealQuantity = 4u"ns",
                              max_nsteps::Int = 1000,
@@ -84,7 +77,7 @@ function simulate_waveforms( mcevents_table::AbstractVector{<:NamedTuple}, sim::
                              verbose::Bool = false ) where {T <: SSDFloat}
     mcevents = TypedTables.Table(mcevents_table)
     #ToDo: Provisional until we support units in CartesianPoints
-    unitless_pos = [collect(_get_unitless_positions_svec(v)) for v in mcevents.pos]
+    unitless_pos = [collect(to_internal_units.(v)) for v in mcevents.pos]
     mcevents = TypedTables.Table((evtno = mcevents.evtno, detno = mcevents.detno, thit  = mcevents.thit, edep  = mcevents.edep, pos   = unitless_pos))
 
     n_total_physics_events = length(mcevents)

From efce34823cbb30383904d11d21ccbefee61d47bc Mon Sep 17 00:00:00 2001
From: Felix Hagemann <hagemann@berkeley.edu>
Date: Fri, 28 Nov 2025 11:08:02 -0800
Subject: [PATCH 35/39] No need to deal with units, the plot recipe will deal
 with this

---
 docs/src/tutorials/geant4_ssd_lit.jl | 2 --
 1 file changed, 2 deletions(-)

diff --git a/docs/src/tutorials/geant4_ssd_lit.jl b/docs/src/tutorials/geant4_ssd_lit.jl
index d47375651..3a7f913d3 100644
--- a/docs/src/tutorials/geant4_ssd_lit.jl
+++ b/docs/src/tutorials/geant4_ssd_lit.jl
@@ -89,8 +89,6 @@ events = run_geant4_simulation(app, N_events)
 
 plot(sim.detector, show_passives = false, size = (500,500), fmt = :png)
 plot!(source_1)
-#ToDo: Put this back when units are supported
-#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

From f6261cac2cde98d8484df416d162b903c25ecd9e Mon Sep 17 00:00:00 2001
From: Felix Hagemann <hagemann@berkeley.edu>
Date: Fri, 28 Nov 2025 11:10:06 -0800
Subject: [PATCH 36/39] Strip position units in `simulate_waveforms` further
 down

---
 src/MCEventsProcessing/MCEventsProcessing.jl | 6 +-----
 1 file changed, 1 insertion(+), 5 deletions(-)

diff --git a/src/MCEventsProcessing/MCEventsProcessing.jl b/src/MCEventsProcessing/MCEventsProcessing.jl
index 8e15fd8c3..51c435bf6 100644
--- a/src/MCEventsProcessing/MCEventsProcessing.jl
+++ b/src/MCEventsProcessing/MCEventsProcessing.jl
@@ -76,10 +76,6 @@ function simulate_waveforms( mcevents_table::AbstractVector{<:NamedTuple}, sim::
                              geometry_check::Bool = false,
                              verbose::Bool = false ) where {T <: SSDFloat}
     mcevents = TypedTables.Table(mcevents_table)
-    #ToDo: Provisional until we support units in CartesianPoints
-    unitless_pos = [collect(to_internal_units.(v)) for v in mcevents.pos]
-    mcevents = TypedTables.Table((evtno = mcevents.evtno, detno = mcevents.detno, thit  = mcevents.thit, edep  = mcevents.edep, pos   = unitless_pos))
-
     n_total_physics_events = length(mcevents)
     Δtime = T(to_internal_units(Δt)) 
     n_contacts = length(sim.detector.contacts)
@@ -193,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

From 861e5dbe35773a412b2c00dd4ee68fa4201a588e Mon Sep 17 00:00:00 2001
From: Felix Hagemann <hagemann@berkeley.edu>
Date: Fri, 28 Nov 2025 11:11:12 -0800
Subject: [PATCH 37/39] Add test for different types of `pos` in
 `simulate_waveforms`

---
 test/test_geant4.jl | 14 ++++++++------
 1 file changed, 8 insertions(+), 6 deletions(-)

diff --git a/test/test_geant4.jl b/test/test_geant4.jl
index 62a247e79..f9a57760d 100644
--- a/test/test_geant4.jl
+++ b/test/test_geant4.jl
@@ -103,13 +103,13 @@ 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))
@@ -121,4 +121,6 @@ end
     @test :waveform in columnnames(wf_static)
     @test length(wf_static) == length(evts_static) * sum(.!ismissing.(sim.weighting_potentials))
 
+    # Result should not depend whether input positions are SVector or CartesianPoint
+    @test wf_static.waveform ≈ wf.waveform
 end

From cfba8d486c3c5cb2fc9aa1114520a06db7ed1fa5 Mon Sep 17 00:00:00 2001
From: Felix Hagemann <hagemann@berkeley.edu>
Date: Fri, 28 Nov 2025 11:42:18 -0800
Subject: [PATCH 38/39] Ensure cluster radius is charge clustering has correct
 units

---
 src/ChargeClustering/ChargeClustering.jl | 6 +++---
 test/test_geant4.jl                      | 5 +++++
 2 files changed, 8 insertions(+), 3 deletions(-)

diff --git a/src/ChargeClustering/ChargeClustering.jl b/src/ChargeClustering/ChargeClustering.jl
index 17bc6437d..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,7 +16,7 @@ function cluster_detector_hits(
     r_edep = similar(edep, 0)
     r_pos = similar(pos, 0)
 
-    ustripped_cradius = to_internal_units(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)
@@ -24,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/test/test_geant4.jl b/test/test_geant4.jl
index f9a57760d..bcbb6fe5e 100644
--- a/test/test_geant4.jl
+++ b/test/test_geant4.jl
@@ -115,6 +115,11 @@ end
     @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

From 92fedb0a4476c2956b47addfd32fcc4bb6ad0e15 Mon Sep 17 00:00:00 2001
From: Felix Hagemann <hagemann@berkeley.edu>
Date: Fri, 28 Nov 2025 11:58:08 -0800
Subject: [PATCH 39/39] Avoid converting `CartesianPoint` to `SVector` when
 writing to HDF5

---
 ext/SolidStateDetectorsLegendHDF5IOExt.jl |  8 +-
 test/test_io.jl                           | 90 ++++++++++++-----------
 2 files changed, 51 insertions(+), 47 deletions(-)

diff --git a/ext/SolidStateDetectorsLegendHDF5IOExt.jl b/ext/SolidStateDetectorsLegendHDF5IOExt.jl
index 1bb23f377..bad97ae46 100644
--- a/ext/SolidStateDetectorsLegendHDF5IOExt.jl
+++ b/ext/SolidStateDetectorsLegendHDF5IOExt.jl
@@ -57,12 +57,12 @@ 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
diff --git a/test/test_io.jl b/test/test_io.jl
index a1105b78c..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