Add geometric validation for electrode termination in charge drift

[zql2021]

Problem
Charge carriers incorrectly terminated when point_types flagged positions as CD_ELECTRODE due to grid discretization, even though they hadn't reached actual electrode geometry.(see
#474)

Solution
Added crossing_pos in det.contacts validation in _check_and_update_position!():
Only terminate when charge actually reaches electrode geometry
Reclassify failed cases as CD_FLOATING_BOUNDARY to continue surface drift

[fhagemann]

Thank you for this initiative!
Do you have any experience how much this would slow down large-scale simulations? Because the geometry checks can become quite resourceful if done too often, especially if the detector has a lot of (complex-geometry) contacts.

[fhagemann]

@zql2021 You had a sketch in #474, can you show a before/after with this PR for the exact same event. And also what the execution time would be for the two?

[zql2021]

@zql2021 You had a sketch in #474, can you show a before/after with this PR for the exact same event. And also what the execution time would be for the two?

Sorry for the late reply. Below are some test code and results.

Test single event ：

Initial position of an event is set at the center of the gap (0, 1.5, 2.8), with 500 keV energy deposition. Considering the charge cloud, it is set to 500 charges using nbcc = NBodyChargeCloud(pos, Edep, 500, number_of_shells = 2).

Theoretically, the collected charge should be 500*1e3/2.95 = 169491.525
Before code modification, the charge collection was incomplete (real collected charge was 165574.144 ).
After code modification, the collected charge is 169491.525, which matches the theoretical value.

begin
    T=Float64
    position = (0.0, 1.5, 2.8)  # mm
    pos = CartesianPoint{T}(position[1] / 1000.0, position[2] / 1000.0, position[3] / 1000.0)# convert mm to m
    Edep = T(500) * u"keV"
    
    nbcc = NBodyChargeCloud(pos, Edep, 500, number_of_shells = 2)
    evt = Event(nbcc)
   
    simulate!(evt, sim, max_nsteps = 10000, Δt = 1u"ns", diffusion=true, self_repulsion=true, verbose = true)
    
    contact5end=ustrip(evt.waveforms[5].signal[end]) 
    contact6end=ustrip(evt.waveforms[6].signal[end])  
    
    theory_signal=500*1e3/2.95
    real_signal=contact5end+contact6end

    println("theory_end", theory_signal)
    println("real_end",real_signal)
end

The following figures show the histogram distribution of the x-coordinates at the end positions of charges (holes here) before and after code modification. The red dashed lines (x=0.25, x=-0.25) represent the electrode boundaries.
Andthe waveform comparison before and after modification.
It can be seen that before modification, due to charge carriers not reaching the electrodes, there was incomplete charge collection.

h_end_x = [path.h_path[end][1] for path in evt.drift_paths]
histogram(h_end_x*1000, bins = 20,  
    xlabel = "PosX (mm)",ylabel = "Counts",
    title = "PosX @ Hole End",
    legend = :top  ,fmt = :png, 
    label="Hole EndPosX", unit=u"mm")
vline!([0.25],  color = :red, 
    linestyle = :dash,  linewidth = 2,label = "X = 0.25 mm")
vline!([-0.25],color = :red, 
    linestyle = :dash, linewidth = 2,label = "X = -0.25 mm" )

begin
using Plots, CSV, DataFrames
df_before = CSV.read("waveforms_before_modify.csv", DataFrame)
df_after = CSV.read("waveforms_after_modify.csv", DataFrame)

p = plot(layout=(2,1), size=(800, 600))

# Contact 5
plot!(p[1], df_before.time_step, df_before.contact_5, 
      label="Before modify", linewidth=2, color=:blue)
plot!(p[1], df_after.time_step, df_after.contact_5, 
      label="After modify", linewidth=2, color=:red, linestyle=:dash)
plot!(p[1], title="Contact 5", 
      xlabel="Time (ns)", ylabel="Amplitude (e)")

# Contact 6 
plot!(p[2], df_before.time_step, df_before.contact_6, 
      label="Before modify", linewidth=2, color=:blue)
plot!(p[2], df_after.time_step, df_after.contact_6, 
      label="After modify", linewidth=2, color=:red, linestyle=:dash)
plot!(p[2], title="Contact 6", 
      xlabel="Time (ns)", ylabel="Amplitude  (e)")

plot!(p, plot_title="Waveform Comparison: Before vs After Modification")
display(p)
end

Speed Test for 1000 events：

Set 1000 events in the gap center (x=0.0) and directly under the contact (x=1.5) respectively , with y fixed at 1.5 and z uniformly sampled between 0.5 and 13.5.

begin
    using Random, Statistics
    T = Float64
    # Random.seed!(123)
    start_time = time()
    for i in 1:1000
        x = 0.   # gap center
        #x = 1.5  # under contact area 
        y= 1.5
        z = 0.5 + rand() * (13.5 - 0.5)  

        position = (x, y, z)  # mm
        pos = CartesianPoint{T}(position[1] / 1000.0, position[2] / 1000.0, position[3] / 1000.0)  
        Edep = T(500) * u"keV"
        
        nbcc = NBodyChargeCloud(pos, Edep, 500, number_of_shells = 2)
        evt = Event(nbcc)
        simulate!(evt, sim, max_nsteps = 10000, Δt = 1u"ns", diffusion=true, self_repulsion=true, verbose = false)
    end
    total_time = time() - start_time
    println("total_time: $(round(total_time, digits=2)) s")
end

In the gap center,
before code modification, simulating 1000 events cost 1244.92 seconds,
after code modification, simulating 1000 events cost 1848.94 seconds.

Under the contact,
before code modification, simulating 1000 events cost 1204.57 seconds,
after code modification, simulating 1000 events cost 1173.17 seconds.

Therefore, after code modification, the simulation time for events under the contact remains nearly the same, while in the gap region, it increased about 48.5% (604.02 seconds). However, this increase is related to the number of charges in the charge cloud (which is 500 in this case).

[fhagemann]

I see, so for events that might be affected by this issue, there seems to be a significant increase in runtime.
Maybe it makes sense to make this extra geometry check optional (we could talk about what the default should be) and allow the user to choose when this geometry check should be applied (?)

[fhagemann]

I added a keyword geometry_check now to all functions that drift charges.
The default is false to obtain the old behavior, but using simulate!(evt, sim, geometry_check = true) or drift_charges!(evt, sim, geometry_check = true) will now additionally perform the geometry check to determine if charge carriers have reached the contact or not,.

[fhagemann]

If that's fine for everyone, I would merge and release this still as 0.10.24 before merging the rest and going for 0.11
@oschulz @claudiaalvgar