Add support for simulation for dead layer-Step3: BulkSurfaceConstantLifetimeChargeTrappingModel

docs/src/man/charge_drift.md

@@ -219,7 +219,7 @@ In this case, the signal is calculated using the Schockley-Ramo theorem, i.e. by
 
 ### `BoggsChargeTrappingModel`
 
-In SolidStateDetectors.jl, the only charge trapping model implemented so far is the one presented in [Boggs _et al._ (2023)](https://doi.org/10.1016/j.nima.2023.168756).
+In SolidStateDetectors.jl, the Boggs charge trapping model implemented is the one presented in [Boggs _et al._ (2023)](https://doi.org/10.1016/j.nima.2023.168756).
 In this model, the charge cloud loses part of its charge at every point `path[i]` of the charge drift, depending on its drift and thermal velocity, as well as the trapping product $[n\sigma_{e/h}]^{-1}$ for electrons and holes.
 The charge signal is then given by the charge-decreased charge cloud reaching the contacts and the charges trapped on the way.
 
@@ -261,6 +261,59 @@ parameters = Dict("parameters" =>
 sim.detector = SolidStateDetector(sim.detector, BoggsChargeTrappingModel{T}(parameters))
 ```
 
+### `ConstantLifetimeChargeTrappingModel`
+This constant-lifetime-based charge trapping model assumes electrons and holes to have a constant free lifetime throughout the crystal.
+In this model, the charge cloud loses part of its charge at every point `path[i]` of the charge drift, depending on the lifetime.
+The charge signal is then given by the charge-decreased charge cloud reaching the contacts and the charges trapped on the way.
+
+This idea is presented in [Dai _et al._ (2023)](https://doi.org/10.1016/j.apradiso.2022.110638).
+
+Besides, the model implemented has two sets of parameters for the sensitive (bulk) and inactive (dead layer/surface layer) volume respectively.
+The inactive layer effect will only be considered if the corresponding `virtual_drift_volume` is defined in the configuration file.
+
+The `ConstantLifetimeChargeTrappingModel` can be applied in the configuration file by adding a field `charge_trapping_model` to the `semiconductor` with `model: ConstantLifetime`, `parameters` defining the constant lifetime and `inactive_layer_geometry`(optional) defining the inactive volume:
+```yaml 
+detectors:
+  - semiconductor:
+      material: #...
+      geometry: #...
+      charge_trapping_model:
+        model: ConstantLifetime
+        parameters:
+          τh: 1ms
+          τe: 1ms
+          τh_inactive: 80ns
+          τe_inactive: 80ns
+        inactive_layer_geometry:
+          tube:
+            r:
+              from: 9.0
+              to: 10.0
+            h: 10.0
+            origin:
+              z: 5.0
+```
+
+The inactive layer could only be modeled with the help of the configuration file.
+But for the sensitive volume, the parameters can always be applied to an already read-in `SolidStateDetector` using for example:
+```julia
+using SolidStateDetectors, Unitful
+T = Float64
+sim = Simulation{T}(SSD_examples[:TrueCoaxial])
+simulate!(sim)
+
+# only model the sensitive volume
+τh = τe = 1u"ms"
+parameters = Dict("parameters" => Dict("τh" => τh, "τe" => τe))
+
+# model both the sensitive and inactive volumes based on the inactive volume geometry defined in the configuration file
+τh = τe = 1u"ms"
+τh_inactive = τe_inactive = 80u"ns"
+parameters = Dict("parameters" => Dict("τh" => τh, "τe" => τe, "τh_inactive" => τh_inactive, "τe_inactive" => τe_inactive), "inactive_layer_geometry" => sim.detector.semiconductor.charge_trapping_model.inactive_layer_geometry)
+
+sim.detector = SolidStateDetector(sim.detector, ConstantLifetimeChargeTrappingModel{T}(parameters))
+```
+
 
 ## Group Effects
 

examples/example_config_files/public_TrueCoaxial.yaml

@@ -32,14 +32,29 @@ detectors:
       value: 1e8cm^-3
     charge_drift_model:
       include: ADLChargeDriftModel/drift_velocity_config.yaml
+    charge_trapping_model:
+      model: ConstantLifetime
+      parameters:
+        τh: 1ms
+        τe: 1ms
+        τh_inactive: 80ns
+        τe_inactive: 80ns
+      inactive_layer_geometry:
+        tube:
+          r:
+            from: 9.0
+            to: 10.0
+          h: 10.0
+          origin:
+            z: 5.0
     geometry:
       cone:
         r:
           from: 1.0
           to: 10.0
-        h: 10 
+        h: 10.0
         origin:
-          z: 5
+          z: 5.0
   contacts:
     - name: "P+"
       id: 1
@@ -49,10 +64,10 @@ detectors:
         tube:
           r:
             from: 0.5
-            to: 1
-          h: 10
+            to: 1.0
+          h: 10.0
           origin:
-            z: 5
+            z: 5.0
     - name: "N+"
       id: 2
       material: HPGe
@@ -61,7 +76,7 @@ detectors:
         tube:
           r:
             from: 10.0
-            to: 10.5
-          h: 10.
+            to: 10.0
+          h: 10.0
           origin:
-            z: 5
+            z: 5.0

src/ChargeTrapping/ChargeTrapping.jl

@@ -121,3 +121,86 @@ end
 
 
 
+
+"""
+    struct ConstantLifetimeChargeTrappingModel{T <: SSDFloat} <: AbstractChargeTrappingModel{T}
+
+This constant-lifetime-based charge trapping model is similar to the Boggs model, which is constant-mean-free-path based.
+
+The model implemented has two sets of parameters for the sensitive (bulk) and inactive (dead layer/surface layer) volume respectively.
+The inactive layer effect will only be considered if the corresponding `inactive_layer_geometry` is defined in the configuration file.
+
+## Fields
+* `τh::T`: Lifetime for holes in bulk volume (sensitive region) (default: `τh = 1ms`).
+* `τe::T`: Lifetime for electrons in bulk volume (default: `τe = 1ms`).
+* `τh_inactive::T`: Lifetime for holes in surface layer (dead layer/inactive layer) (default: `τh_inactive = 1ms`).
+* `τe_inactive::T`: Lifetime for electrons in surface layer (default: `τe_inactive = 1ms`).
+* `inactive_layer_geometry`: The geometry of the inactive layer. When this is defined, the model will consider the inactive effect.
+> Note: All τ must be much bigger than `Δt`.
+
+See also [Charge Trapping Models](@ref).
+"""
+struct ConstantLifetimeChargeTrappingModel{T <: SSDFloat, G <: Union{<:AbstractGeometry, Nothing}} <: AbstractChargeTrappingModel{T}
+    τh::T
+    τe::T
+    τh_inactive::T
+    τe_inactive::T
+    inactive_layer_geometry::G
+end
+
+function _calculate_signal( 
+        ctm::ConstantLifetimeChargeTrappingModel{T}, 
+        path::AbstractVector{CartesianPoint{T}}, 
+        pathtimestamps::AbstractVector{T}, 
+        charge::T,          
+        wpot::Interpolations.Extrapolation{T, 3}, 
+        S::CoordinateSystemType
+    )::Vector{T} where {T <: SSDFloat}
+
+    tmp_signal::Vector{T} = Vector{T}(undef, length(pathtimestamps))
+
+    running_sum::T = zero(T)
+    q::T = charge
+
+    inactive_layer_exist = !isnothing(ctm.inactive_layer_geometry)
+    τ::T = ifelse(charge > 0, ctm.τh, ctm.τe)
+    τ_inactive::T = ifelse(charge > 0, ctm.τh_inactive, ctm.τe_inactive)
+
+    @inbounds for i in eachindex(tmp_signal)
+        Δt::T = (i > 1) ? (pathtimestamps[i] - pathtimestamps[i-1]) : zero(T)
+        τi::T = (inactive_layer_exist && in(path[i], ctm.inactive_layer_geometry)) ? τ_inactive : τ
+        Δq::T = q * Δt/τi
+        q -= Δq
+        w::T = i > 1 ? get_interpolation(wpot, path[i], S) : zero(T)
+        running_sum += w * Δq
+        tmp_signal[i] = running_sum + w * q
+    end
+
+    tmp_signal
+end
+
+
+ConstantLifetimeChargeTrappingModel(args...; T::Type{<:SSDFloat}, kwargs...) = ConstantLifetimeChargeTrappingModel{T}(args...; kwargs...)
+function ConstantLifetimeChargeTrappingModel{T}(config_dict::AbstractDict = Dict()) where {T <: SSDFloat}
+    τh::T = ustrip(u"s", 1u"ms")
+    τe::T = ustrip(u"s", 1u"ms")
+    τh_inactive::T = ustrip(u"s", 80u"ns")
+    τe_inactive::T = ustrip(u"s", 80u"ns")
+
+    if haskey(config_dict, "model") && !haskey(config_dict, "parameters")
+        throw(ConfigFileError("`ConstantLifetimeChargeTrappingModel` does not have `parameters`"))
+    end
+
+    parameters = haskey(config_dict, "parameters") ? config_dict["parameters"] : config_dict
+    inactive_layer_geometry = haskey(config_dict, "inactive_layer_geometry") ? config_dict["inactive_layer_geometry"] : nothing
+
+    allowed_keys = ("τh", "τe", "τh_inactive", "τe_inactive")
+    k = filter(k -> !(k in allowed_keys), keys(parameters))
+    !isempty(k) && @warn "The following keys will be ignored: $(k).\nAllowed keys are: $(allowed_keys)"
+
+    if haskey(parameters, "τh")   τh =     _parse_value(T, parameters["τh"], internal_time_unit) end
+    if haskey(parameters, "τe")   τe =     _parse_value(T, parameters["τe"], internal_time_unit) end
+    if haskey(parameters, "τh_inactive")   τh_inactive =     _parse_value(T, parameters["τh_inactive"], internal_time_unit) end
+    if haskey(parameters, "τe_inactive")   τe_inactive =     _parse_value(T, parameters["τe_inactive"], internal_time_unit) end
+    ConstantLifetimeChargeTrappingModel{T, typeof(inactive_layer_geometry)}(τh, τe, τh_inactive, τe_inactive, inactive_layer_geometry)
+end

src/SolidStateDetector/Semiconductor.jl

@@ -81,17 +81,26 @@ function Semiconductor{T}(dict::AbstractDict, input_units::NamedTuple, outer_tra
         T(293)
     end
 
+    inner_transformations = parse_CSG_transformation(T, dict, input_units)
+    transformations = combine_transformations(inner_transformations, outer_transformations)
+    geometry = Geometry(T, dict["geometry"], input_units, transformations)
+
     charge_trapping_model = if haskey(dict, "charge_trapping_model") &&
                                haskey(dict["charge_trapping_model"], "model") &&
                                dict["charge_trapping_model"]["model"] == "Boggs"
         BoggsChargeTrappingModel{T}(dict["charge_trapping_model"], temperature = temperature)
+    elseif haskey(dict, "charge_trapping_model") &&
+                               haskey(dict["charge_trapping_model"], "model") &&
+                               dict["charge_trapping_model"]["model"] == "ConstantLifetime"
+        constant_lifetime_ctm_dict = deepcopy(dict["charge_trapping_model"])
+        if haskey(constant_lifetime_ctm_dict, "inactive_layer_geometry")
+            constant_lifetime_ctm_dict["inactive_layer_geometry"] = Geometry(T, dict["charge_trapping_model"]["inactive_layer_geometry"], input_units, transformations)
+        end
+        ConstantLifetimeChargeTrappingModel{T}(constant_lifetime_ctm_dict)
     else
         NoChargeTrappingModel{T}()
     end
 
-    inner_transformations = parse_CSG_transformation(T, dict, input_units)
-    transformations = combine_transformations(inner_transformations, outer_transformations)
-    geometry = Geometry(T, dict["geometry"], input_units, transformations)
     return Semiconductor(temperature, material, impurity_density_model, charge_drift_model, charge_trapping_model, geometry)
 end
 

src/SolidStateDetectors.jl

@@ -68,7 +68,7 @@ export ElectricPotential, PointTypes, EffectiveChargeDensity, DielectricDistribu
 export apply_initial_state!
 export calculate_electric_potential!, calculate_weighting_potential!, calculate_electric_field!, calculate_drift_fields!
 export ElectricFieldChargeDriftModel, ADLChargeDriftModel, IsotropicChargeDriftModel
-export NoChargeTrappingModel, BoggsChargeTrappingModel
+export NoChargeTrappingModel, BoggsChargeTrappingModel, ConstantLifetimeChargeTrappingModel
 export get_active_volume, is_depleted, estimate_depletion_voltage
 export calculate_stored_energy, calculate_mutual_capacitance, calculate_capacitance_matrix
 export simulate_waveforms

test/test_charge_drift_models.jl

@@ -31,7 +31,7 @@ nbcc = NBodyChargeCloud(pos, Edep, 40, radius = T(0.0005), number_of_shells = 2)
     @test isapprox( signalsum, T(2), atol = 5e-3 )
 end
 
-@timed_testset "Charge Trapping" begin
+@timed_testset "Charge Trapping: BoggsChargeTrappingModel" begin
     sim.detector = SolidStateDetector(sim.detector, BoggsChargeTrappingModel{T}())
     evt = Event(pos, Edep)
     timed_simulate!(evt, sim)
@@ -81,6 +81,112 @@ end
     end
 end
 
+@timed_testset "Charge Trapping: ConstantLifetimeChargeTrappingModel" begin
+    # test 1: parse the lifetimes and the inactive layer geometry
+    config_dict = SolidStateDetectors.parse_config_file(SSD_examples[:TrueCoaxial])
+    simA = @test_nowarn Simulation{T}(config_dict)
+    @testset "Parse the TrueCoaxial config file" begin
+        @test simA.detector.semiconductor.charge_trapping_model isa ConstantLifetimeChargeTrappingModel{T}
+        @test simA.detector.semiconductor.charge_trapping_model.τh == T(1e-3)
+        @test simA.detector.semiconductor.charge_trapping_model.τe == T(1e-3)
+        @test simA.detector.semiconductor.charge_trapping_model.τh_inactive == T(8e-8)
+        @test simA.detector.semiconductor.charge_trapping_model.τe_inactive  == T(8e-8)
+        @test simA.detector.semiconductor.charge_trapping_model.inactive_layer_geometry.origin == CartesianPoint{T}(0.0, 0.0, 0.005)
+        r0, r1 = T.((0.009, 0.01))
+        @test simA.detector.semiconductor.charge_trapping_model.inactive_layer_geometry.r == tuple((r0, r1), (r0, r1))
+        @test simA.detector.semiconductor.charge_trapping_model.inactive_layer_geometry.hZ == T(0.005)
+    end
+
+    simA = Simulation{T}(SSD_examples[:TrueCoaxial])
+    timed_simulate!(simA, convergence_limit = 1e-5, device_array_type = device_array_type, refinement_limits = [0.2, 0.1, 0.05], verbose = false)
+    simA_inactive_layer_geometry=deepcopy(simA.detector.semiconductor.charge_trapping_model.inactive_layer_geometry)
+
+    # test2: only model the bulk volume, test the bulk signals while varying the lifetime: τ
+    pos = CartesianPoint{T}(0.005,0,0.005); Edep = 1u"eV"
+    signalsum_list = []
+    τ_list = (10 .^ (-3:0.3:0))*u"ms"
+    for τ in τ_list
+        parameters = Dict("parameters" => Dict("τh" => τ, "τe" => τ))
+        trapping_model=ConstantLifetimeChargeTrappingModel{T}(parameters)
+        simA.detector = SolidStateDetector(simA.detector, trapping_model)
+        evt = Event(pos, Edep)
+        timed_simulate!(evt, simA)
+        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(simA.detector.semiconductor.material)))
+        push!(signalsum_list, signalsum)
+    end
+    @info signalsum_list
+    @test all(signalsum_list .< T(2.0))
+    @test all(diff(signalsum_list) .> 0)
+
+    # test3: model both the bulk/inactive volumes
+    ## test 3.1: the bulk/inactive signals while varying the lifetimes: τ, τ_inactive
+    pos_bulk = CartesianPoint{T}(0.0085,0,0.005); Edep = 1u"eV"
+    @test !in(pos_bulk, simA_inactive_layer_geometry)
+    signalsum_list_bulk = []
+    τ_list = (10 .^ (-2:0.2:0))*u"ms"
+    pos_inactive = CartesianPoint{T}(0.0095,0,0.005); Edep = 1u"eV"
+    @test in(pos_inactive, simA_inactive_layer_geometry)
+    signalsum_list_inactive = []
+    τ_inactive_list = τ_list/100
+    for (τ,τ_inactive) in zip(τ_list, τ_inactive_list)
+        parameters = Dict("parameters" => Dict("τh" => τ, "τe" => τ, "τh_inactive" => τ_inactive, "τe_inactive" => τ_inactive),
+            "inactive_layer_geometry" => simA_inactive_layer_geometry)
+        trapping_model=ConstantLifetimeChargeTrappingModel{T}(parameters)
+        simA.detector = SolidStateDetector(simA.detector, trapping_model)
+        evt_bulk = Event(pos_bulk , Edep)
+        timed_simulate!(evt_bulk, simA)
+        signalsum_bulk  = T(0)
+        for i in 1:length(evt_bulk.waveforms)
+            signalsum_bulk  += abs(ustrip(evt_bulk.waveforms[i].signal[end]))
+        end
+        signalsum_bulk  *= inv(ustrip(SolidStateDetectors._convert_internal_energy_to_external_charge(simA.detector.semiconductor.material)))
+        push!(signalsum_list_bulk, signalsum_bulk)
+
+        evt_inactive = Event(pos_inactive , Edep)
+        timed_simulate!(evt_inactive, simA)
+        signalsum_inactive  = T(0)
+        for i in 1:length(evt_inactive.waveforms)
+            signalsum_inactive  += abs(ustrip(evt_inactive.waveforms[i].signal[end]))
+        end
+        signalsum_inactive  *= inv(ustrip(SolidStateDetectors._convert_internal_energy_to_external_charge(simA.detector.semiconductor.material)))
+        push!(signalsum_list_inactive, signalsum_inactive)
+    end
+    @info signalsum_list_bulk
+    @test all(signalsum_list_bulk .< T(2.0))
+    @test all(diff(signalsum_list_bulk) .> 0)
+    @info signalsum_list_inactive
+    @test all(signalsum_list_inactive .< T(2.0))
+    @test all(diff(signalsum_list_inactive) .> 0)
+    @test all(signalsum_list_bulk .> signalsum_list_inactive)
+
+    ## test 3.2: the inactive layer signals while varying the depth
+    τ, τ_inactive = 1u"ms", 100u"ns"
+    parameters = Dict("parameters" => Dict("τh" => τ, "τe" => τ, "τh_inactive" => τ_inactive, "τe_inactive" => τ_inactive),
+        "inactive_layer_geometry" => simA_inactive_layer_geometry)
+    trapping_model=ConstantLifetimeChargeTrappingModel{T}(parameters)
+    simA.detector = SolidStateDetector(simA.detector, trapping_model)
+    signalsum_list_inactive = []
+    for depth in (0.1:0.1:0.9)/1000
+        pos_inactive = CartesianPoint{T}(0.01-depth,0,0.005); Edep = 1u"eV"
+        @test in(pos_inactive, simA_inactive_layer_geometry)
+        evt_inactive = Event(pos_inactive , Edep)
+        timed_simulate!(evt_inactive, simA, diffusion=false, Δt=1u"ns")
+        signalsum_inactive  = T(0)
+        for i in 1:length(evt_inactive.waveforms)
+            signalsum_inactive  += abs(ustrip(evt_inactive.waveforms[i].signal[end]))
+        end
+        signalsum_inactive  *= inv(ustrip(SolidStateDetectors._convert_internal_energy_to_external_charge(simA.detector.semiconductor.material)))
+        push!(signalsum_list_inactive, signalsum_inactive)
+    end
+    @info signalsum_list_inactive
+    @test all(signalsum_list_inactive .< T(2.0))
+    @test all(diff(signalsum_list_inactive) .> 0)
+end
+
 @timed_testset "Test completeness of charge drift models" begin
     for c in InteractiveUtils.subtypes(AbstractChargeDriftModel)
         if isstructtype(c)