Merge branch 'dev'

src/DetectorGeometries/Coax.jl

@@ -40,7 +40,7 @@ mutable struct Coax{T<:AbstractFloat} <: SolidStateDetector{T}
 
     ### Segmentaion
     n_total_contacts::Int
-    n_repetitive_segments::Int
+    # n_repetitive_segments::Int
     n_individual_segments::Int
 
     segmentation_r_ranges::Array{Tuple{T,T},1}
@@ -117,13 +117,13 @@ function Coax{T}(config_file::Dict)::Coax where T<:AbstractFloat
 
     ### Segmentation
     c.n_total_contacts = config_file["segmentation"]["n_contacts_total"]
-    c.n_repetitive_segments = config_file["segmentation"]["n_repetitive_segments"]
+    # c.n_repetitive_segments = config_file["segmentation"]["n_repetitive_segments"]
     c.n_individual_segments = config_file["segmentation"]["n_individual_segments"]
     c.custom_grouping = config_file["segmentation"]["custom_grouping"]
 
     # c.segmentation_boundaryWidth_vertical = round(c.geometry_unit_factor *config_file["segmentation"]["repetitive_segment"]["boundaryWidth"]["vertical"],digits=5)
     # c.segmentation_boundaryWidth_horizontal = round(config_file["segmentation"]["repetitive_segment"]["boundaryWidth"]["vertical"]*c.geometry_unit_factor *180/(π*c.crystal_radius), digits=5)
-    construct_segmentation_arrays_from_repetitive_segment(c,config_file)
+    # construct_segmentation_arrays_from_repetitive_segment(c, config_file)
     c.n_individual_segments > 0 ? construct_segmentation_arrays_for_individual_segments(c,config_file) : nothing
     construct_grouped_channels_array(c, config_file)
     construct_floating_boundary_arrays(c)

src/DetectorGeometries/DetectorGeometries.jl

@@ -848,7 +848,7 @@ function is_boundary_point(d::SolidStateDetector, r::T, φ::T, z::T, rs::Vector{
     end
     digits::Int=6
     if  d.borehole_modulation == true
-        idx_r_closest_gridpoint_to_borehole = find_closest_gridpoint(rs, d.borehole_radius + d.borehole_ModulationFunction(φ))
+        idx_r_closest_gridpoint_to_borehole = searchsortednearest(rs, d.borehole_radius + d.borehole_ModulationFunction(φ))
         idx_r = findfirst(x->x==r,rs)
         bore_hole_r = rs[idx_r_closest_gridpoint_to_borehole]
     else
@@ -874,7 +874,7 @@ function is_boundary_point(d::SolidStateDetector, r::T, φ::T, z::T, rs::Vector{
                     if d.segmentation_types[iseg]=="Tubs"
                         rv = true
                     else
-                        if isapprox(rs[find_closest_gridpoint(rs,analytical_taper_r_from_φz(φ,z,
+                        if isapprox(rs[searchsortednearest(rs,analytical_taper_r_from_φz(φ,z,
                             d.segmentation_types[iseg],
                             d.segmentation_r_ranges[iseg],
                             d.segmentation_phi_ranges[iseg],
@@ -971,21 +971,19 @@ function analytical_taper_r_from_φz(φ,z,orientation,r_bounds,φ_bounds,z_bound
     r
 end
 
-# function modulated_borehole_r()
-
-function find_closest_gridpoint(array::Vector{<:T}, value::T)::Int where {T<:Real}
-    high::Int = findfirst(x -> x > value, array)
-    idx::Int = -1
-    if high==1
-        idx=1
-    elseif sum(array[high-1:high])/2 > value
-        idx = high - 1
+function searchsortednearest(a::AbstractArray{T, 1}, x::T)::Int where {T <: Real}
+    idx::Int = searchsortedfirst(a, x)
+    if (idx == 1) return idx end
+    if (idx > length(a)) return length(a) end
+    if (a[idx] == x) return idx end
+    if (abs(a[idx] - x) < abs(a[idx - 1] - x))
+        return idx
     else
-        idx = high
+        return idx - 1
     end
-    return idx
 end
 
+
 function get_segment_idx(d::SolidStateDetector,r::Real,φ::Real,z::Real)
     digits=6
     for iseg in 1:size(d.segment_bias_voltages,1)
@@ -1008,7 +1006,7 @@ function get_segment_idx(d::SolidStateDetector,r::Real,φ::Real,z::Real)
 end
 function get_segment_idx(d::SolidStateDetector,r::Real,φ::Real,z::Real,rs::Vector{<:Real})
     digits=6
-    idx_r_closest_gridpoint_to_borehole = find_closest_gridpoint(rs,d.borehole_radius+d.borehole_ModulationFunction(φ))
+    idx_r_closest_gridpoint_to_borehole = searchsortednearest(rs,d.borehole_radius+d.borehole_ModulationFunction(φ))
     idx_r = findfirst(x->x==r,rs)
     bore_hole_r = rs[idx_r_closest_gridpoint_to_borehole]
     for iseg in 1:size(d.segment_bias_voltages,1)
@@ -1057,7 +1055,9 @@ function get_charge_density(detector::SolidStateDetector, r::Real, ϕ::Real, z::
     top_net_charge_carrier_density::T = detector.charge_carrier_density_top * 1e10 * 1e6  #  1/cm^3 -> 1/m^3
     bot_net_charge_carrier_density::T = detector.charge_carrier_density_bot * 1e10 * 1e6  #  1/cm^3 -> 1/m^3
     slope::T = (top_net_charge_carrier_density - bot_net_charge_carrier_density) / detector.crystal_length
-    return bot_net_charge_carrier_density + z * slope
+    ρ::T = bot_net_charge_carrier_density + z * slope
+    # ρ = ρ - (r - detector.borehole_radius) * 0.2 * ρ / (detector.crystal_radius - detector.borehole_radius) 
+    return ρ 
 end
 
 function json_to_dict(inputfile::String)::Dict

src/Grids/CylindricalGrid.jl

@@ -144,26 +144,27 @@ function CylindricalGrid(detector::SolidStateDetector, r_arr::AbstractArray, φ_
     return CylindricalGrid{T}(cyclic, r_arr, φ_arr, z_arr, only_2d=only_2d)
 end
 
-function CylindricalGrid(detector::SolidStateDetector; init_grid_spacing::Real=2e-3, cyclic::Real=2π, only_2d::Bool = false)::CylindricalGrid
+function CylindricalGrid(detector::SolidStateDetector; init_grid_spacing::Vector{<:Real}=[2e-3, 2π / 72, 2e-3], cyclic::Real=2π, only_2d::Bool = false)::CylindricalGrid
     T = get_precision_type(detector)
 
     only_2d::Bool = cyclic == 0 ? true : false
 
-    grid_spacing::T = init_grid_spacing
-
     r_start::T = 0
-    r_stop::T = detector.crystal_radius + 4 * grid_spacing
-    r_step::T = grid_spacing # 1mm
+    r_step::T = init_grid_spacing[1] # 1mm
+    r_stop::T = detector.crystal_radius + 4 * r_step
     r_length::Int = div(r_stop - r_start, r_step)
     r_arr::Vector{T} = collect(range(r_start, stop=r_stop, length=r_length))
-
+
+    five_degree_in_rad::T = 2π / 72
+
     φ_start::T = 0
-    φ_step::T = grid_spacing / r_arr[2]
+    φ_step::T = init_grid_spacing[2]
 
     # φ_length::Int = div(φ_stop - φ_start, φ_step)
-    φ_length = 8
-    φ_step = cyclic / φ_length 
+    # φ_length = 8
+    # φ_step = cyclic / φ_length 
     φ_stop::T = cyclic - φ_step
+    φ_length::Int = div(φ_stop - φ_start, φ_step)
 
     if isodd(φ_length) φ_length += 1 end
     if φ_length == 0 φ_length = 2 end
@@ -173,9 +174,9 @@ function CylindricalGrid(detector::SolidStateDetector; init_grid_spacing::Real=2
         collect(range(φ_start, stop=φ_stop, length=φ_length) )
     end
 
-    z_start::T =  -4 * grid_spacing
-    z_step::T = grid_spacing # 1 mm
-    z_stop::T = detector.crystal_length + 4 * grid_spacing
+    z_step::T = init_grid_spacing[3] # 1 mm
+    z_start::T =  -4 * z_step
+    z_stop::T = detector.crystal_length + 4 * z_step
     z_length::Int = div(z_stop - z_start, z_step)
     if isodd(z_length)
         z_length + 1

src/MaterialProperties/MaterialProperties.jl

@@ -17,3 +17,21 @@ material_properties[:HPGe] = (
     ρ = 5.323Unitful.g / (Unitful.cm^3), # u"g/cm^3" throws warnings in precompilation
     name = "High Purity Germanium"
 )
+
+
+# These values might just be approximations
+material_properties[:Si] = (
+    E_ionisation = 3.62u"eV",
+    f_fano = 0.11,
+    ϵ_r = 11.7,
+    ρ = 2.3290Unitful.g / (Unitful.cm^3), # u"g/cm^3" throws warnings in precompilation
+    name = "Silicon"
+)
+
+material_properties[:LAr] = (
+    name = "Silicon",
+    E_ionisation = 0u"eV",
+    f_fano = 0.107,
+    ϵ_r = 1.505,
+    ρ = 1.396Unitful.g / (Unitful.ml), # u"g/cm^3" throws warnings in precompilation
+)

src/Potentials/CalculateElectricPotential.jl

@@ -19,6 +19,7 @@ There are serveral `<keyword arguments>` which can be used to tune the computati
 - `return_2π_grid::Bool`: If set to `true`, the grid is extended to 2π in `φ` after it is computated. Default is `true`. This keyword is ignored if the simulation is in 2D. Use `extent_2D_grid_to_3D()` function to extend a 2D grid into 3D.
 - `max_n_iterations::Int`: Set the maximum number of iterations which are performed after each grid refinement. Default is `10000`. If set to `-1` there will be no limit.
 - `verbose::Bool=true`: Boolean whether info output is produced or not.
+- `init_grid_spacing::Vector{<:Real}`: Initial spacing of the grid. Default is [2e-3, 2π / 72, 2e-3] <=> [2mm, 5 degree, 2mm ]
 
 # Additional Information
 - The function always turns the detector into a p-type detector to compute the potential. At the end it turns the signs to make it n-type again if it was n-type.
@@ -29,7 +30,7 @@ function calculate_electric_potential(  det::SolidStateDetector;
                                         max_refinements::Int=3,
                                         refinement_limits::Vector{<:Real}=[1e-4, 1e-4, 1e-4],
                                         min_grid_spacing::Vector{<:Real}=[1e-4, 1e-2, 1e-4],
-                                        init_grid_spacing::Real=2e-3, # 2 mm
+                                        init_grid_spacing::Vector{<:Real}=[2e-3, 2π / 72, 2e-3], # 2mm, 5 degree, 2 mm
                                         depletion_handling::Bool=false,
                                         nthreads::Int=Base.Threads.nthreads(),
                                         sor_consts::Vector{<:Real}=[1.4, 1.85],
@@ -41,7 +42,7 @@ function calculate_electric_potential(  det::SolidStateDetector;
     bias_voltage = T(det.segment_bias_voltages[1]);
     refinement_limits = T.(refinement_limits)
     min_grid_spacing = T.(min_grid_spacing)
-    init_grid_spacing = T(init_grid_spacing)
+    init_grid_spacing = T.(init_grid_spacing)
     sor_consts = T.(sor_consts)
     refine::Bool = max_refinements > 0 ? true : false
 
@@ -169,10 +170,13 @@ function calculate_electric_potential(  det::SolidStateDetector;
     abs_min::T = abs(minimum(cg.potential))
     if det_tmp.bulk_type == :ntype cg.potential *= -1 end
 
+    abs_max_bv::T = abs(maximum(det.segment_bias_voltages))
+    abs_min_bv::T = abs(minimum(det.segment_bias_voltages))
+
     # Checks
     if bias_voltage_abs < abs_max + abs_min
     # if abs_max > abs(det.segment_bias_voltages[1])
-        @warn "Detector not fully depleted at this bias voltage ($(det.segment_bias_voltages[1]) V). At least on grid point has a higher (or lower) potential value ($(sign(det.segment_bias_voltages[1]) * abs_max) V). This should not be."
+        @warn "Detector not fully depleted at this bias voltage ($( max(abs_max_bv, abs_min_bv) ) V). At least on grid point has a higher (or lower) potential value ($( max(abs_max, abs_min) ) V). This should not be. However, small overshoots might be due to over relaxation and not full convergence."
     end
     if !in(det.segment_bias_voltages[1], cg.potential)
         @warn "Something is wrong. Applied bias voltage ($(det.segment_bias_voltages[1]) V) does not occur in the grid. This should not be since it is a fixed value."

src/Potentials/CalculateWeightingPotential.jl

@@ -19,6 +19,7 @@ There are serveral `<keyword arguments>` which can be used to tune the computati
 - `sor_consts::Vector{<:Real}`: Two element array. First element contains the SOR constant for `r` = 0. Second contains the constant at the outer most grid point in `r`. A linear scaling is applied in between. First element should be smaller than the second one and both should be ∈ [1.0, 2.0].
 - `max_n_iterations::Int`: Set the maximum number of iterations which are performed after each grid refinement. Default is `10000`. If set to `-1` there will be no limit.
 - `verbose::Bool=true`: Boolean whether info output is produced or not.
+- `init_grid_spacing::Vector{<:Real}`: Initial spacing of the grid. Default is [2e-3, 2π / 72, 2e-3] <=> [2mm, 5 degree, 2mm ]
 
 # Additional Information
 
@@ -29,7 +30,7 @@ function calculate_weighting_potential( det::SolidStateDetector, channel::Int;#c
                                         max_refinements::Int=3,
                                         refinement_limits::Vector{<:Real}=[1e-4, 1e-4, 1e-4],
                                         min_grid_spacing::Vector{<:Real}=[1e-4, 1e-2, 1e-4],
-                                        init_grid_spacing::Real=2e-3, # 2 mm
+                                        init_grid_spacing::Vector{<:Real}=[2e-3, 2π / 72, 2e-3], # 2mm, 5 degree, 2 mm
                                         depletion_handling::Bool=false,
                                         nthreads::Int=Base.Threads.nthreads(),
                                         sor_consts::Vector{<:Real}=[1.4, 1.85],
@@ -41,7 +42,7 @@ function calculate_weighting_potential( det::SolidStateDetector, channel::Int;#c
     bias_voltage::T = 1
     refinement_limits = T.(refinement_limits)
     min_grid_spacing = T.(min_grid_spacing)
-    init_grid_spacing = T(init_grid_spacing)
+    init_grid_spacing = T.(init_grid_spacing)
     sor_consts = T.(sor_consts)
     refine::Bool = max_refinements > 0 ? true : false
 

src/Potentials/PrecalculatedWeightsCylindricalRedBlack.jl

@@ -72,8 +72,7 @@ function PrecalculatedWeightsCylindricalRedBlack(detector::SolidStateDetector, g
     odd::Int  = 1
     even::Int = 2
 
-    # @inbounds begin
-    begin
+    @inbounds begin
         nr::Int = length(grid.r)
         nφ::Int = length(grid.φ)
         nz::Int = length(grid.z)
@@ -289,15 +288,6 @@ function PrecalculatedWeightsCylindricalRedBlack(detector::SolidStateDetector, g
                     inside_detector::Bool = contains(detector, pos_r, pos_φ, pos_z)
                     if inside_detector point_types[pwinds...] += pn_junction_bit end
 
-                    # deprecated:
-                    # ρ = if inside_detector
-                    #     get_charge_density(detector, pos_r, pos_φ, pos_z) * effective_electron_charge_vacuum
-                    #         This charge is not fully correct.The charge should be a weighted sum of the neighbouring (8) cells.
-                    #         This will be correct later. The error should be small for now.
-                    # else
-                    #     zero(T)
-                    # end
-
                     charge[pwinds...] = dV * ρ
                 end
             end
@@ -414,11 +404,18 @@ function get_active_volume(grid::CylindricalGrid, pts::PointTypes{T}) where {T}
             end
         end
     end
+    if grid.cyclic > 0
+        active_volume *= 2π / grid.cyclic
+    end
 
     f::T = 10^6
     return active_volume * f * Unitful.cm * Unitful.cm * Unitful.cm
 end
 
+function is_fully_depleted(pts::PointTypes{T})::Bool where {T}
+    return !(undepleted_bit in (pts.pointtypes .& undepleted_bit))
+end
+
 
 function get_epsilon_r_and_charge(detector::SolidStateDetector, grid::CylindricalGrid; sor_consts::Vector{<:Real}=[1.4, 1.85], only_2d::Bool=false)
     T = get_precision_type(detector)
@@ -443,8 +440,7 @@ function get_epsilon_r_and_charge(detector::SolidStateDetector, grid::Cylindrica
     odd::Int  = 1
     even::Int = 2
 
-    # @inbounds begin
-    begin
+    @inbounds begin
         nr::Int = length(grid.r)
         nφ::Int = length(grid.φ)
         nz::Int = length(grid.z)
@@ -521,17 +517,9 @@ function get_epsilon_r_and_charge(detector::SolidStateDetector, grid::Cylindrica
                 while pos_φ < 0 
                     pos_φ += grid.cyclic
                 end
-                # if iz == 12
-                #     # @show pos_φ, pos_z
-                # end
-                # pos_φ += 0.1
                 for ir in 2:size(epsilon_r, 1)
                     pos_r = mpr[ir]
 
-                    if iz == 12 && iφ == 1 
-                        ct = contains(detector, pos_r, pos_φ, pos_z)
-                        @show pos_z, pos_φ, pos_r, ct
-                    end
                     if contains(detector, pos_r, pos_φ, pos_z)
                         epsilon_r[ir, iφ, iz]::T = detector_material_ϵ_r
                         charge_tmp[ir, iφ, iz]::T = get_charge_density(detector, pos_r, pos_φ, pos_z) * effective_electron_charge_vacuum

src/SolidStateDetectors.jl

@@ -57,8 +57,4 @@ include("IO/IO.jl")
 
 include("examples.jl")
 
-function check_borehole(x,testfunction)
-return 1 + testfunction(x)
-end
-
 end # module