Refactor boule impurity densities (#477)

* refactor boule impurity densities

* default impurity param vals

Author: hervasa2

src/ConstructiveSolidGeometry/Units.jl

@@ -11,6 +11,7 @@ to_internal_units(x::AngleQuantity)  = ustrip(internal_angle_unit, x)
 to_internal_units(x::Quantity{<:Real, Unitful.𝐋^(-1)}) = ustrip(internal_length_unit^(-1), x) # charge trapping
 to_internal_units(x::Quantity{<:Real, Unitful.𝐋^(-3)}) = ustrip(internal_length_unit^(-3), x) # densities
 to_internal_units(x::Quantity{<:Real, Unitful.𝐋^(-4)}) = ustrip(internal_length_unit^(-4), x) # density gradients
+to_internal_units(x::Quantity{<:Real, Unitful.𝐋^(-5)}) = ustrip(internal_length_unit^(-5), x) # density quadratic components
 to_internal_units(x::AbstractArray) = to_internal_units.(x)
 
 from_internal_units(x::Real, unit::Unitful.Units{<:Any, Unitful.𝐋}) = uconvert(unit, x * internal_length_unit)

src/ImpurityDensities/BouleImpurityDensities/BouleImpurityDensities.jl

@@ -1,2 +1,4 @@
+include("ParBouleImpurityDensity.jl")
 include("LinBouleImpurityDensity.jl")
-include("LinExpBouleImpurityDensity.jl")
+include("LinExpBouleImpurityDensity.jl")
+include("ParExpBouleImpurityDensity.jl")

src/ImpurityDensities/BouleImpurityDensities/LinBouleImpurityDensity.jl

@@ -16,6 +16,13 @@ struct LinBouleImpurityDensity{T <: SSDFloat} <: AbstractImpurityDensity{T}
     det_z0::T
 end
 
+function ImpurityDensity(T::DataType, t::Val{:linear_boule}, dict::AbstractDict, input_units::NamedTuple)
+    a::T = haskey(dict, "a") ? _parse_value(T, dict["a"], input_units.length^(-3)) : T(0)
+    b::T = haskey(dict, "b") ? _parse_value(T, dict["b"], input_units.length^(-4)) : T(0)
+    det_z0::T = haskey(dict, "det_z0") ? _parse_value(T, dict["det_z0"], input_units.length) : T(0)
+    LinBouleImpurityDensity{T}(a, b, det_z0)
+end
+
 function get_impurity_density(idm::LinBouleImpurityDensity, pt::AbstractCoordinatePoint{T})::T where {T}
     cpt = CartesianPoint(pt)
     z = cpt[3]

src/ImpurityDensities/BouleImpurityDensities/LinExpBouleImpurityDensity.jl

@@ -1,29 +1,39 @@
 """
     struct LinExpBouleImpurityDensity{T <: SSDFloat} <: AbstractImpurityDensity{T}
 
-a + b*z + c*exp((z-L)/tau)
+a + b*z + n*exp((z-l)/m)
 Impurity density model which assumes a linear + exponential component in impurity density in z.
 
 ## Fields
 * `a:T`: impurity density values at the boule origin.
 * `b::T`: linear slope in `z` direction.
-* `c::T`: exponential coefficient.
-* `L::T`: exponential offset.
-* `tau::T`: exponential scale factor.
+* `n::T`: exponential coefficient.
+* `l::T`: exponential offset.
+* `m::T`: exponential scale factor.
 * `det_z0::T`: z coordinate of the detector origin in boule coordinates. The z-direction of the detector is opposite to the z-direction of the boule coordinates.
 """
 
 struct LinExpBouleImpurityDensity{T <: SSDFloat} <: AbstractImpurityDensity{T}
     a::T
     b::T 
-    c::T 
-    L::T
-    tau::T
+    n::T 
+    l::T
+    m::T
     det_z0::T
 end
 
+function ImpurityDensity(T::DataType, t::Val{:linear_exponential_boule}, dict::AbstractDict, input_units::NamedTuple)
+    a::T = haskey(dict, "a") ? _parse_value(T, dict["a"], input_units.length^(-3)) : T(0)
+    b::T = haskey(dict, "b") ? _parse_value(T, dict["b"], input_units.length^(-4)) : T(0)
+    n::T = haskey(dict, "n") ? _parse_value(T, dict["n"], input_units.length^(-3)) : T(0)
+    l::T = haskey(dict, "l") ? _parse_value(T, dict["l"], input_units.length) : T(0)
+    m::T = haskey(dict, "m") ? _parse_value(T, dict["m"], input_units.length) : T(1)
+    det_z0::T = haskey(dict, "det_z0") ? _parse_value(T, dict["det_z0"], input_units.length) : T(0)
+    LinExpBouleImpurityDensity{T}(a, b, n, l, m, det_z0)
+end
+
 function get_impurity_density(idm::LinExpBouleImpurityDensity, pt::AbstractCoordinatePoint{T})::T where {T}
     cpt = CartesianPoint(pt)
     z = cpt[3]
-    idm.a + idm.b * (idm.det_z0 - z) + idm.c * exp((idm.det_z0 - z - idm.L)/idm.tau)
+    idm.a + idm.b * (idm.det_z0 - z) + idm.n * exp((idm.det_z0 - z - idm.l)/idm.m)
 end

src/ImpurityDensities/BouleImpurityDensities/ParBouleImpurityDensity.jl

@@ -0,0 +1,33 @@
+"""
+    struct ParBouleImpurityDensity{T <: SSDFloat} <: AbstractImpurityDensity{T}
+
+a + b*z + c*z^2
+Impurity density model which assumes a linear + exponential component in impurity density in z.
+
+## Fields
+* `a:T`: impurity density values at the boule origin.
+* `b::T`: linear slope in `z` direction.
+* `c::T`: quadratic coefficient.
+* `det_z0::T`: z coordinate of the detector origin in boule coordinates. The z-direction of the detector is opposite to the z-direction of the boule coordinates.
+"""
+
+struct ParBouleImpurityDensity{T <: SSDFloat} <: AbstractImpurityDensity{T}
+    a::T
+    b::T 
+    c::T 
+    det_z0::T
+end
+
+function ImpurityDensity(T::DataType, t::Val{:parabolic_boule}, dict::AbstractDict, input_units::NamedTuple)
+    a::T = haskey(dict, "a") ? _parse_value(T, dict["a"], input_units.length^(-3)) : T(0)
+    b::T = haskey(dict, "b") ? _parse_value(T, dict["b"], input_units.length^(-4)) : T(0)
+    c::T = haskey(dict, "c") ? _parse_value(T, dict["c"], input_units.length^(-5)) : T(0)
+    det_z0::T = haskey(dict, "det_z0") ? _parse_value(T, dict["det_z0"], input_units.length) : T(0)
+    ParExpBouleImpurityDensity{T}(a, b, c, det_z0)
+end
+
+function get_impurity_density(idm::ParBouleImpurityDensity, pt::AbstractCoordinatePoint{T})::T where {T}
+    cpt = CartesianPoint(pt)
+    z = cpt[3]
+    idm.a + idm.b * (idm.det_z0 - z) + idm.c * (idm.det_z0 - z)^2
+end

src/ImpurityDensities/BouleImpurityDensities/ParExpBouleImpurityDensity.jl

@@ -0,0 +1,42 @@
+"""
+    struct ParExpBouleImpurityDensity{T <: SSDFloat} <: AbstractImpurityDensity{T}
+
+a + b*z + c*z^2 + n*exp((z-l)/m)
+Impurity density model which assumes a linear + exponential component in impurity density in z.
+
+## Fields
+* `a:T`: impurity density values at the boule origin.
+* `b::T`: linear slope in `z` direction.
+* `c::T`: quadratic coefficient.
+* `n::T`: exponential coefficient.
+* `l::T`: exponential offset.
+* `m::T`: exponential scale factor.
+* `det_z0::T`: z coordinate of the detector origin in boule coordinates. The z-direction of the detector is opposite to the z-direction of the boule coordinates.
+"""
+
+struct ParExpBouleImpurityDensity{T <: SSDFloat} <: AbstractImpurityDensity{T}
+    a::T
+    b::T 
+    c::T 
+    n::T
+    l::T
+    m::T
+    det_z0::T
+end
+
+function ImpurityDensity(T::DataType, t::Val{:parabolic_exponential_boule}, dict::AbstractDict, input_units::NamedTuple)
+    a::T = haskey(dict, "a") ? _parse_value(T, dict["a"], input_units.length^(-3)) : T(0)
+    b::T = haskey(dict, "b") ? _parse_value(T, dict["b"], input_units.length^(-4)) : T(0)
+    c::T = haskey(dict, "c") ? _parse_value(T, dict["c"], input_units.length^(-5)) : T(0)
+    n::T = haskey(dict, "n") ? _parse_value(T, dict["n"], input_units.length^(-3)) : T(0)
+    l::T = haskey(dict, "l") ? _parse_value(T, dict["l"], input_units.length) : T(0)
+    m::T = haskey(dict, "m") ? _parse_value(T, dict["m"], input_units.length) : T(1)
+    det_z0::T = haskey(dict, "det_z0") ? _parse_value(T, dict["det_z0"], input_units.length) : T(0)
+    ParExpBouleImpurityDensity{T}(a, b, c, n, l, m, det_z0)
+end
+
+function get_impurity_density(idm::ParExpBouleImpurityDensity, pt::AbstractCoordinatePoint{T})::T where {T}
+    cpt = CartesianPoint(pt)
+    z = cpt[3]
+    idm.a + idm.b * (idm.det_z0 - z) + idm.c * (idm.det_z0 - z)^2 + idm.n * exp((idm.det_z0 - z - idm.l)/idm.m)
+end

src/SolidStateDetectors.jl

@@ -78,7 +78,7 @@ export Simulation, simulate!
 export Event, drift_charges!
 export add_baseline_and_extend_tail
 export NBodyChargeCloud
-export LinearImpurityDensity, LinBouleImpurityDensity, LinExpBouleImpurityDensity
+export LinBouleImpurityDensity, ParBouleImpurityDensity, LinExpBouleImpurityDensity, ParExpBouleImpurityDensity
 
 using Unitful: RealOrRealQuantity as RealQuantity
 const SSDFloat = Union{Float16, Float32, Float64}

test/BEGe_02.yaml

@@ -0,0 +1,76 @@
+name: Point-contact detector
+units:
+  length: mm
+  angle: deg
+  potential: V
+  temperature: K
+grid:
+  coordinates: cylindrical
+  axes:
+    r:
+      to: 60
+      boundaries: inf
+    phi:
+      from: 0
+      to: 0
+      boundaries:
+        left: periodic
+        right: periodic
+    z:
+      from: -20
+      to: 60
+      boundaries:
+        left: inf
+        right: inf
+medium: vacuum
+
+detectors:
+- semiconductor:
+    material: HPGe
+    impurity_density:
+      name: linear_boule
+      a: -1e7
+      b: -1e5
+      det_z0: 120
+    geometry:
+      translate:
+        tube:
+          r: 40
+          h: 40
+        z: 20
+  contacts:
+  - name: Core
+    material: HPGe
+    id: 1
+    potential: -3000
+    geometry:
+      tube:
+        r: 3
+        h: 0.3
+        origin:
+          z: 39.85
+  - name: Mantle
+    material: HPGe
+    id: 2
+    potential: 0
+    geometry:
+      union:
+      - tube:
+          r:
+            from: 0
+            to: 40
+          h: 0
+      - tube:
+          r:
+            from: 40
+            to: 40
+          h: 40
+          origin:
+            z: 20
+      - tube:
+          r:
+            from: 10
+            to: 40
+          h: 0
+          origin:
+            z: 40

test/test_boule_impurity_densities.jl

@@ -15,13 +15,17 @@ T = Float32
 
     @test det_ρ0 == boule_ρ0 + boule_gradient * det_z0
 
-    boule_c = T(-2e15)
-    boule_L = T(0.05)
-    boule_tau = T(0.03)
-    idm = LinExpBouleImpurityDensity{T}(boule_ρ0, boule_gradient, boule_c, boule_L, boule_tau, det_z0)
+    sim2 = Simulation{T}("BEGe_02.yaml")
+
+    @test sim.detector.semiconductor.impurity_density_model == sim2.detector.semiconductor.impurity_density_model
+
+    boule_n = T(-2e15)
+    boule_l = T(0.05)
+    boule_m = T(0.03)
+    idm = LinExpBouleImpurityDensity{T}(boule_ρ0, boule_gradient, boule_n, boule_l, boule_m, det_z0)
     sim.detector = SolidStateDetector(sim.detector, idm)
 
     det_ρ0 = SolidStateDetectors.get_impurity_density(sim.detector.semiconductor.impurity_density_model, CylindricalPoint{T}(0,0,0))
 
-    @test det_ρ0 == boule_ρ0 + boule_gradient * det_z0 + boule_c * exp((det_z0 - boule_L)/boule_tau)
+    @test det_ρ0 == boule_ρ0 + boule_gradient * det_z0 + boule_n * exp((det_z0 - boule_l)/boule_m)
 end