Merge pull request #553 from hervasa2/main

Use vectors instead of points to `move_charges_inside_semiconductor`

Author: fhagemann

Project.toml

@@ -79,7 +79,7 @@ Rotations = "1"
 SpecialFunctions = "1.6.1, 2"
 StaticArrays = "1"
 Statistics = "<0.0.1, 1"
-StatsBase = "0.33.7, 0.34"
+StatsBase = "0.33.14, 0.34"
 Suppressor = "0.2"
 Tables = "1.1"
 TypedTables = "1.2"

src/Event/Event.jl

@@ -156,16 +156,16 @@ function move_charges_inside_semiconductor!(
                         proj::T = (i - charge_center) ⋅ (l - charge_center) / sum((l - charge_center).^2)
                         if 0 ≤ proj ≤ 1
                             if i in det.semiconductor
-                                locations[n][m] = (1 - fraction * proj) * i + fraction * proj * charge_center
+                                locations[n][m] = cartesian_zero + (1 - fraction * proj) * (i - cartesian_zero) + fraction * proj * (charge_center - cartesian_zero)
                             end
                         end
                     end
                 end
             end
             charge_center_new = mean(locations[n], Weights(edep_weights))
             if verbose
-                @warn "$(sum(.!idx_in)) charges of the charge cloud at $(round.(charge_center, digits = (T == Float64 ? 12 : 6)))"*
-                " are outside. Moving them inside...\nThe new charge center is at $(round.(charge_center_new, digits = (T == Float64 ? 12 : 6))).\n"
+                @warn "$(sum(.!idx_in)) charges of the charge cloud at $(round.(charge_center - cartesian_zero, digits = (T == Float64 ? 12 : 6)))"*
+                " are outside. Moving them inside...\nThe new charge center is at $(round.(charge_center_new - cartesian_zero, digits = (T == Float64 ? 12 : 6))).\n"
             end
         end
     end

src/MCEventsProcessing/MCEventsProcessing.jl

@@ -135,7 +135,8 @@ function _convertEnergyDepsToChargeDeps(
         particle_type::Type{PT} = Gamma,
         radius::AbstractVector{<:AbstractVector{<:Union{<:Real, <:LengthQuantity}}} = map(e -> radius_guess.(to_internal_units.(e), particle_type), edep),
         number_of_carriers::Int = 1, number_of_shells::Int = 1, 
-        max_interaction_distance::Union{<:Real, <:LengthQuantity} = NaN # is ignored here
+        max_interaction_distance::Union{<:Real, <:LengthQuantity} = NaN, # is ignored here
+        kwargs...
     ) where {T <: SSDFloat, PT <: ParticleType}
 
     charge_clouds = broadcast(
@@ -148,7 +149,7 @@ function _convertEnergyDepsToChargeDeps(
                     number_of_shells = number_of_shells,
                     radius = T(to_internal_units(radius[iEdep_indep][i_together]))
                 )
-                move_charges_inside_semiconductor!([nbcc.locations], [nbcc.energies], det)
+                move_charges_inside_semiconductor!([nbcc.locations], [nbcc.energies], det; kwargs...)
                 nbcc
             end,
             eachindex(edep[iEdep_indep])
@@ -172,7 +173,7 @@ function _simulate_charge_drifts( mcevents::TypedTables.Table, sim::Simulation{T
                                   geometry_check::Bool,
                                   verbose::Bool ) where {T <: SSDFloat}
     @showprogress map(mcevents) do phyevt
-        locations, edeps = _convertEnergyDepsToChargeDeps(phyevt.pos, phyevt.edep, sim.detector; number_of_carriers, number_of_shells, max_interaction_distance)
+        locations, edeps = _convertEnergyDepsToChargeDeps(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

test/test_charge_drift_models.jl

@@ -31,6 +31,27 @@ nbcc = NBodyChargeCloud(pos, Edep, 40, radius = T(0.0005), number_of_shells = 2)
     @test isapprox( signalsum, T(2), atol = 5e-3 )
 end
 
+@timed_testset "Move charges inside semiconductor" begin
+
+    # define a charge cloud very outside of the detector => should throw an AssertionError
+    nbcc_top = NBodyChargeCloud(CartesianPoint{T}(0.01,0,0.081), Edep, 10, radius = T(0.002), number_of_shells = 1)
+    evt = Event(nbcc_top)
+    @test_throws AssertionError timed_simulate!(evt, sim, self_repulsion = true, diffusion = true)
+
+    # define a charge cloud very close to the top surface of the example inverted coax detector (top = 80mm) => should throw a warning
+    nbcc_top = NBodyChargeCloud(CartesianPoint{T}(0.01,0,0.079), Edep, 10, radius = T(0.002), number_of_shells = 1)
+    evt = Event(nbcc_top)
+    @test_logs (:warn, r".*Moving them inside.*") timed_simulate!(evt, sim, self_repulsion = true, diffusion = true)
+
+    signalsum = T(0)
+    for i in 1:length(evt.waveforms)
+        signalsum += abs(ustrip(evt.waveforms[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
+
+
 @timed_testset "Charge Trapping: BoggsChargeTrappingModel" begin
     sim.detector = SolidStateDetector(sim.detector, BoggsChargeTrappingModel{T}())
     evt = Event(pos, Edep)