Merge pull request #489 from JuliaPhysics/point-types

Change AbstractCoordinatePoint to not be a FieldVector

Author: fhagemann

.github/workflows/Downgrade.yml

@@ -19,7 +19,7 @@ jobs:
     runs-on: ubuntu-latest
     strategy:
       matrix:
-        version: ['1.10', '1']
+        version: ['1.10', '1.12']
     steps:
       - uses: actions/checkout@v5
       - uses: julia-actions/setup-julia@v2

Project.toml

@@ -14,6 +14,7 @@ Format = "1fa38f19-a742-5d3f-a2b9-30dd87b9d5f8"
 GPUArrays = "0c68f7d7-f131-5f86-a1c3-88cf8149b2d7"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953"
+InverseFunctions = "3587e190-3f89-42d0-90ee-14403ec27112"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
 LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
@@ -58,6 +59,7 @@ GPUArrays = "8, 9, 10, 11"
 Geant4 = "0.1.13, 0.2"
 Interpolations = "0.14, 0.15, 0.16"
 IntervalSets = "0.6, 0.7"
+InverseFunctions = "0.1"
 JSON = "0.21.2, 1.1"
 KernelAbstractions = "0.8, 0.9"
 LaTeXStrings = "1.1"
@@ -69,7 +71,7 @@ Parameters = "0.12"
 PolygonOps = "0.1.1"
 Polynomials = "2, 3, 4"
 ProgressMeter = "1.5"
-RadiationDetectorSignals = "0.3.5"
+RadiationDetectorSignals = "0.3.10"
 Random = "<0.0.1, 1"
 RecipesBase = "1"
 Requires = "1.1.3"

devenv/README.md

@@ -3,7 +3,8 @@
 This directory contains the Julia project environments `cpu` and `cuda`
 that can be used when developing SolidStateDetectors with Julia >= v1.11.
 They contain all direct, test and doc-gen dependencies of SolidStateDetectors,
-plus BenchmarkTools and Cthulhu.
+plus BenchmarkTools, Cthulhu and JLD2 (for snapshots and regression tests
+during development).
 
 Note: These environments can't be used with Julia versions <= v1.10, as they
 use a `[sources]` section in the `Project.toml` to ensure SolidStateDetectors

devenv/cuda/Project.toml

@@ -18,6 +18,7 @@ InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953"
 InverseFunctions = "3587e190-3f89-42d0-90ee-14403ec27112"
+JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
 LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"

docs/Project.toml

@@ -1,4 +1,5 @@
 [deps]
+ArraysOfArrays = "65a8f2f4-9b39-5baf-92e2-a9cc46fdf018"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 Geant4 = "559df036-b7a0-42fd-85df-7d5dd9d70f44"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"

docs/src/tutorials/dead_layer_simulation_lit.jl

@@ -11,14 +11,14 @@ using SolidStateDetectors
 T = Float64
 
 # the geometry parameters of the model for following display
-mm=T(1/1000)
-det_rin=1*mm
-det_z=det_r=10*mm
-z_draw=det_z/2
-pn_r=8.957282*mm # this one was calculated by searching the zero impurity point (displayed in the following section)
+mm = T(1/1000)
+det_rin = 1*mm
+det_z = det_r = 10*mm
+z_draw = det_z/2
+pn_r = 8.957282*mm # this one was calculated by searching the zero impurity point (displayed in the following section)
 
 sim = Simulation{T}(SSD_examples[:TrueCoaxial])
-cfn=SSD_examples[:TrueCoaxial]
+cfn = SSD_examples[:TrueCoaxial]
 print(open(f -> read(f, String), cfn))
 plot(sim.detector, xunit = u"mm", yunit = u"mm", zunit = u"mm")
 #jl savefig("tutorial_det_dl.pdf") # hide
@@ -28,10 +28,10 @@ plot(sim.detector, xunit = u"mm", yunit = u"mm", zunit = u"mm")
 
 # ## Display the impurity profile defined
 # > Above position of the PN junction boundary `pn_r` is calculate by searching the point on the curve where the density is zero.
-r_list=(0*mm):(0.01*mm):det_r
+r_list = (0*mm):(0.01*mm):det_r
 imp_list = T[]
 for r in r_list
-    pt=CylindricalPoint{T}(r, 0, z_draw)
+    pt = CylindricalPoint{T}(r, 0, z_draw)
     push!(imp_list, SolidStateDetectors.get_impurity_density(sim.detector.semiconductor.impurity_density_model, pt))
 end
 imp_list = imp_list ./ 1e6 # in cm^-3
@@ -47,12 +47,11 @@ vline!([pn_r/mm], lw = 2, ls = :dash, color = :darkred, label = "PN junction bou
 using SolidStateDetectors: Electron, Hole
 cdm = sim.detector.semiconductor.charge_drift_model
 depth_list = 0:(0.01*mm):(det_r-pn_r)
-hole_mobility_list=[]
-electron_mobility_list=[]
+hole_mobility_list = []
+electron_mobility_list = []
 for depth in depth_list
-    r=det_r-depth
-    pt = CylindricalPoint{T}([r, 0, z_draw])
-    pt = CartesianPoint{T}(pt)
+    r = det_r-depth
+    pt = CartesianPoint{T}(r, 0, z_draw)
     push!(hole_mobility_list, SolidStateDetectors.calculate_mobility(cdm, pt, Hole))
     push!(electron_mobility_list, SolidStateDetectors.calculate_mobility(cdm, pt, Electron))
 end
@@ -70,9 +69,9 @@ plot!(ylabel = "Mobility [cm\$^2\$/V/s]", xlabel = "Depth to surface [mm]",
 calculate_electric_potential!(sim, max_n_iterations = 10, grid = Grid(sim), verbose = false, depletion_handling = true)
 g = sim.electric_potential.grid
 ax1, ax2, ax3 = g.axes
-bulk_tick_dis=0.05*mm
-dl_tick_dis=0.01*mm
-user_additional_ticks_ax1=sort(vcat(ax1.interval.left:bulk_tick_dis:pn_r, pn_r:dl_tick_dis:ax1.interval.right))
+bulk_tick_dis = 0.05*mm
+dl_tick_dis = 0.01*mm
+user_additional_ticks_ax1 = sort(vcat(ax1.interval.left:bulk_tick_dis:pn_r, pn_r:dl_tick_dis:ax1.interval.right))
 user_ax1 = typeof(ax1)(ax1.interval, user_additional_ticks_ax1)
 user_g = typeof(g)((user_ax1, ax2, ax3))
 calculate_electric_potential!(sim, refinement_limits = 0.1, use_nthreads = 8, grid = user_g, depletion_handling = true)
@@ -81,7 +80,7 @@ calculate_weighting_potential!(sim, 1, use_nthreads = 8, depletion_handling = tr
 calculate_weighting_potential!(sim, 2, use_nthreads = 8, depletion_handling = true);
 plot(
     begin
-        imp=plot(sim.imp_scale, φ = 0, xunit = u"mm", yunit = u"mm", title = "impurity scale")
+        imp = plot(sim.imp_scale, φ = 0, xunit = u"mm", yunit = u"mm", title = "impurity scale")
         vline!([pn_r/mm], lw = 2, ls = :dash, color = :darkred, label = "PN junction boundary", legendfontsize = 6)
     end,
     begin
@@ -109,25 +108,25 @@ pulse_list = []
 totTime = 5 # us
 totEnergy = 1 # 1 keV --> simulating ~339 carrier pairs
 max_nsteps = Int(totTime * 1000)
-N=Int(totEnergy*1000÷2.95)
+N = Int(totEnergy*1000÷2.95)
 for depth in depth_list
-    r=det_r-depth
+    r = det_r-depth
     energy_depos = fill(T(2.95/1000), N) * u"keV"
-    starting_positions = repeat([CylindricalPoint{T}(r, deg2rad(0), z_draw)], N)
+    starting_positions = repeat([CartesianPoint{T}(r, deg2rad(0), z_draw)], N)
     evt = Event(starting_positions, energy_depos);
     simulate!(evt, sim, Δt = 1u"ns", max_nsteps = max_nsteps, diffusion = true, end_drift_when_no_field = false, self_repulsion = false)
-    charge=ustrip(evt.waveforms[1].signal)
+    charge = ustrip(evt.waveforms[1].signal)
     push!(pulse_list, charge)
 end
-pulse_plot=plot()
+pulse_plot = plot()
 eff_list = []
 for (i, depth) in enumerate(depth_list)
-    depth=round(depth/mm, digits = 1)
+    depth = round(depth/mm, digits = 1)
     plot!(pulse_list[i], label = "Depth: $(depth) mm", lw = 2, xlabel = "Time [ns]", ylabel = "Amplitude [e]",
         legend = :topright, grid = :on, minorgrid = :on, frame = :box)
     push!(eff_list, maximum(pulse_list[i])/N)
 end
-cce_plot=plot(depth_list/mm, eff_list, xlabel = "Depth to surface [mm]", lw = 2, ylabel = "Charge collection efficiency",
+cce_plot = plot(depth_list/mm, eff_list, xlabel = "Depth to surface [mm]", lw = 2, ylabel = "Charge collection efficiency",
     frame = :box, grid = :on, minorgrid = :on, xticks = 0:0.2:1.2, yticks = 0:0.1:1, dpi = 500, color = :black, label = "")
 plot(pulse_plot, cce_plot, layout = (1, 2), size = (1000, 400), margin = 5Plots.mm)
 #jl savefig("tutorial_pulse_cce_dl.pdf") # hide

docs/src/tutorials/geant4_ssd_lit.jl

@@ -12,7 +12,7 @@ using Geant4
 # The extension features a function that creates a `Geant4.G4JLApplcation` object from a SSD `Simulation` object and a particle source.
 
 using Plots
-using Unitful
+using ArraysOfArrays, Unitful
 
 # Two types of particle source are pre-defined in `SolidStateDetectors`:
 # 
@@ -89,7 +89,7 @@ events = run_geant4_simulation(app, N_events)
 
 plot(sim.detector, show_passives = false, size = (500,500), fmt = :png)
 plot!(source_1)
-plot!(CartesianPoint.(broadcast(p -> ustrip.(u"m", p), events[1:1000].pos.data)), ms = 0.5, msw = 0, color=:black, label = "")
+plot!(ustrip.(u"m", 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

ext/Geant4/io_gdml.jl

@@ -26,7 +26,7 @@ using LightXML
 # Add <position> to <define> section, referenced in the geometry definition (in <solids>) via the name
 function create_position(e::AbstractConstructiveGeometry, x_define::XMLElement, name::String, v::Bool)
     # Calculate relative position of the volumes contained in the parent volume
-    x, y, z = parse_origin(e.b) .- parse_origin(e.a)
+    x, y, z = parse_origin(e.b) - parse_origin(e.a)
 
     xpos = new_child(x_define, "position")
     set_attributes(xpos, OrderedDict(
@@ -515,9 +515,10 @@ function add_to_world(x_world::XMLElement, x_define::XMLElement, e::AbstractGeom
 
     x_pos_ref = new_child(x_world, "positionref")
     set_attribute(x_pos_ref, "ref", "pos_" * name)
-        
+    
     x_pos_ref = new_child(x_define, "position")
-    x, y, z = parse_origin(e)
+    parsed_e = parse_origin(e)
+    x, y, z = parsed_e.x, parsed_e.y, parsed_e.z
     set_attributes(x_pos_ref, OrderedDict(
         "name" => "pos_" * name,
         "x" => x,

ext/SolidStateDetectorsGeant4Ext.jl

@@ -219,7 +219,10 @@ function SolidStateDetectors.run_geant4_simulation(app::G4JLApplication, number_
         next!(p)
         evtno += 1
     end
-    return RadiationDetectorSignals.group_by_evtno(Table(evts))
+
+    evts_with_points = map(evt -> merge(evt, (pos = cartesian_zero + evt.pos,)), evts)
+
+    return RadiationDetectorSignals.group_by_evtno(Table(evts_with_points))
 end 
 
 

ext/SolidStateDetectorsLegendHDF5IOExt.jl

@@ -6,6 +6,8 @@ import ..LegendHDF5IO
 
 using ..SolidStateDetectors
 using ..SolidStateDetectors: RealQuantity, SSDFloat, to_internal_units, chunked_ranges, LengthQuantity
+using ArraysOfArrays
+import Tables
 using TypedTables, Unitful
 using Format
 
@@ -24,7 +26,7 @@ end
 
 
 
-function SolidStateDetectors.simulate_waveforms( mcevents::TypedTables.Table, sim::Simulation{T},
+function SolidStateDetectors.simulate_waveforms( mcevents::AbstractVector{<:NamedTuple}, sim::Simulation{T},
                              output_dir::AbstractString, 
                              output_base_name::AbstractString = "generated_waveforms";
                              chunk_n_physics_events::Int = 1000, 
@@ -42,19 +44,22 @@ function SolidStateDetectors.simulate_waveforms( mcevents::TypedTables.Table, si
     Δtime = T(to_internal_units(Δt)) 
     n_contacts = length(sim.detector.contacts)
     @info "Detector has $(n_contacts) contact(s)"
-    contacts = sim.detector.contacts
     if !ispath(output_dir) mkpath(output_dir) end
     nfmt(i::Int) = format(i, zeropadding = true, width = length(digits(n_total_physics_events)))
     evt_ranges = chunked_ranges(n_total_physics_events, chunk_n_physics_events)
-    @info "-> $(sum(length.(mcevents.edep))) energy depositions to simulate."
+    @info "-> $(length(mcevents)) events to simulate."
 
     for evtrange in evt_ranges
         ofn = joinpath(output_dir, "$(output_base_name)_evts_$(nfmt(first(evtrange)))-$(nfmt(last(evtrange))).h5")
         @info "Now simulating $(evtrange) and storing it in\n\t \"$ofn\""
         mcevents_sub = simulate_waveforms(mcevents[evtrange], sim; Δt, max_nsteps, diffusion, self_repulsion, number_of_carriers, number_of_shells, 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,)))
+
         LegendHDF5IO.lh5open(ofn, "w") do h5f
-            LegendHDF5IO.writedata(h5f.data_store, "generated_waveforms", mcevents_sub)
+            LegendHDF5IO.writedata(h5f.data_store, "generated_waveforms", new_mcevents_sub)
         end        
     end    
 end