diff --git a/docs/src/man/charge_drift.md b/docs/src/man/charge_drift.md index 7f6ca2c40..b02660ba7 100644 --- a/docs/src/man/charge_drift.md +++ b/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 diff --git a/examples/example_config_files/public_TrueCoaxial.yaml b/examples/example_config_files/public_TrueCoaxial.yaml index 2550881cf..4976e1018 100644 --- a/examples/example_config_files/public_TrueCoaxial.yaml +++ b/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 \ No newline at end of file + z: 5.0 \ No newline at end of file diff --git a/src/ChargeTrapping/ChargeTrapping.jl b/src/ChargeTrapping/ChargeTrapping.jl index b26513692..ee914d9a5 100644 --- a/src/ChargeTrapping/ChargeTrapping.jl +++ b/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 \ No newline at end of file diff --git a/src/SolidStateDetector/Semiconductor.jl b/src/SolidStateDetector/Semiconductor.jl index adb536089..a90932e16 100644 --- a/src/SolidStateDetector/Semiconductor.jl +++ b/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 diff --git a/src/SolidStateDetectors.jl b/src/SolidStateDetectors.jl index 092971631..0e9bf5128 100644 --- a/src/SolidStateDetectors.jl +++ b/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 diff --git a/test/test_charge_drift_models.jl b/test/test_charge_drift_models.jl index 58d66806b..c2907e06d 100644 --- a/test/test_charge_drift_models.jl +++ b/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)