Merge pull request #431 from Padniuk/main

Add IsotropicChargeDriftModel and add more materials

Project.toml

@@ -50,7 +50,7 @@ ArraysOfArrays = "0.5, 0.6"
 Clustering = "0.15"
 Distributions = "0.24.5, 0.25"
 FillArrays = "0.9, 0.10, 0.11, 0.12, 0.13, 1"
-Format = "1.2, 1.3"
+Format = "1.2"
 GPUArrays = "8, 9, 10"
 Geant4 = "0.1.13"
 Interpolations = "0.14, 0.15"
@@ -68,16 +68,16 @@ Polynomials = "2, 3, 4"
 ProgressMeter = "1.5"
 RadiationDetectorSignals = "0.3.5"
 Random = "<0.0.1, 1"
-RecipesBase = "1.0, 1.1, 1.2, 1.3"
-Requires = "1.1, 1.2, 1.3"
-Rotations = "1.0"
+RecipesBase = "1"
+Requires = "1.1"
+Rotations = "1"
 StaticArrays = "1"
 Statistics = "<0.0.1, 1"
 StatsBase = "0.33.7, 0.34"
 Suppressor = "0.2"
 Tables = "1.1"
 TypedTables = "1.2"
-Unitful = "1.11, 1.12, 1.13, 1.14, 1.15, 1.16, 1.17, 1.18, 1.19"
+Unitful = "1.11"
 UnitfulAtomic = "1"
 YAML = "0.4.2"
 julia = "1.10"

examples/example_config_files/IsotropicChargeDriftModel/mobilities.yaml

@@ -0,0 +1,4 @@
+model: IsotropicChargeDriftModel
+mobilities:
+  e: 1000cm^2/(V*s)
+  h: 1000cm^2/(V*s)

ext/Geant4/io_gdml.jl

@@ -463,6 +463,8 @@ end
     elseif material == "Polyethylene Naphthalate" return "G4_PEN"
     elseif material == "Polytetrafluorethylen" return "G4_PTFE"
     elseif material == "Cadmium zinc telluride" return "G4_CdZnTe"
+    elseif material == "CsPbBr3" return "G4_CsPbBr3"
+    elseif material == "Lead" return "G4_Pb"
     end
     throw("Material characteristics for \"$(material)\" not defined yet")
 end

ext/Geant4/materials.xml

@@ -64,4 +64,26 @@
         <fraction n="0.0278" ref="G4_Zn" />
         <fraction n="0.5422" ref="G4_Te" />
     </material>
+
+    <element name="G4_Cs" Z="55">
+        <atom value="132.9054" />
+    </element>
+    <element name="G4_Pb" Z="82">
+        <atom value="207.2" />
+    </element>
+    <element name="G4_Br" Z="35">
+        <atom value="79.904" />
+    </element>
+
+    <material name="G4_CsPbBr3" formula="CsPbBr3">
+        <D value="4.73" unit="g/cm3" />
+        <fraction n="0.2292" ref="G4_Cs" />
+        <fraction n="0.3574" ref="G4_Pb" />
+        <fraction n="0.4134" ref="G4_Br" />
+    </material>
+
+    <material name="G4_Pb" Z="82">
+        <D value="11.34" />
+        <atom value="207.2" />
+    </material>
 </materials>

src/ChargeDriftModels/ChargeDriftModels.jl

@@ -2,4 +2,5 @@ abstract type AbstractChargeDriftModel{T <: SSDFloat} end
 abstract type AbstractTemperatureModel{T <: SSDFloat} end
 
 include("ElectricFieldChargeDriftModel/ElectricFieldChargeDriftModel.jl")
-include("ADLChargeDriftModel/ADLChargeDriftModel.jl")
+include("ADLChargeDriftModel/ADLChargeDriftModel.jl")
+include("IsotropicChargeDriftModel/IsotropicChargeDriftModel.jl")

src/ChargeDriftModels/IsotropicChargeDriftModel/IsotropicChargeDriftModel.jl

@@ -0,0 +1,52 @@
+"""
+    struct IsotropicChargeDriftModel{T <: SSDFloat} <: AbstractChargeDriftModel{T}
+        
+Charge drift model in which the electrons and holes drift along the electric field with a constant mobility in m²/Vs
+
+## Fields
+- `μ_e::Union{RealQuantity, AbstractString}`: Mobility of the electrons in m²/Vs.
+- `μ_h::Union{RealQuantity, AbstractString}`: Mobility of the holes in m²/Vs.
+
+"""
+struct IsotropicChargeDriftModel{T <: SSDFloat} <: AbstractChargeDriftModel{T}
+    μ_e::T
+    μ_h::T
+    IsotropicChargeDriftModel{T}(
+        μ_e::Union{RealQuantity, AbstractString}, 
+        μ_h::Union{RealQuantity, AbstractString}
+    ) where {T <: SSDFloat} = new{T}(_parse_value.(T, (μ_e, μ_h), internal_mobility_unit)...)
+end
+
+IsotropicChargeDriftModel(args...; T::Type = Float32, kwargs...) = IsotropicChargeDriftModel{T}(args...; kwargs...)
+
+const default_icd_config_file = joinpath(get_path_to_example_config_files(), "IsotropicChargeDriftModel/mobilities.yaml")
+
+function IsotropicChargeDriftModel{T}(config_filename::AbstractString = default_icd_config_file) where {T <: SSDFloat}
+    IsotropicChargeDriftModel(parse_config_file(config_filename))
+end
+
+function IsotropicChargeDriftModel{T}(config::AbstractDict) where {T <: SSDFloat}
+    if !haskey(config, "mobilities")
+        throw(ConfigFileError("IsotropicChargeDriftModel config file needs entry 'mobilities'."))
+    end
+
+    for axis in ("e", "h")
+        if !haskey(config["mobilities"], axis)
+            throw(ConfigFileError("IsotropicChargeDriftModel config file needs entry 'mobilities/$axis'."))
+        end
+    end
+
+    IsotropicChargeDriftModel{T}(config["mobilities"]["e"], config["mobilities"]["h"])
+end
+
+@fastmath function getVe(fv::SVector{3, T}, cdm::IsotropicChargeDriftModel{T}) where {T <: SSDFloat}
+    @inbounds begin 
+        -cdm.μ_e*fv
+    end
+end
+
+@fastmath function getVh(fv::SVector{3, T}, cdm::IsotropicChargeDriftModel{T}) where {T <: SSDFloat}
+    @inbounds begin 
+        cdm.μ_h*fv
+    end
+end

src/MaterialProperties/MaterialProperties.jl

@@ -89,6 +89,21 @@ material_properties[:CdZnTe] = (
     ϵ_r = 10.9,
     ρ = 5.78u"g*cm^-3"
 )
+
+abstract type CsPbBr3 <: AbstractDriftMaterial end
+Symbol(::Type{CsPbBr3}) = :CsPbBr3
+material_properties[:CsPbBr3] = (
+    name = "CsPbBr3",
+    E_ionisation = 6.9u"eV",
+    ϵ_r = 16.46,
+    ρ = 4.73u"g*cm^-3"
+)
+
+material_properties[:Pb] = (
+    name = "Lead",
+    ϵ_r = 1e6, # high value for lead as for conductor. Usage of lead is shielding.
+    ρ = 11.35u"g*cm^-3"
+)
 # Add new materials above this line
 # and just put different spellings into the dict `materials` below
 
@@ -101,7 +116,8 @@ materials = Dict{String, Symbol}(
     "Al"  => :Al,
     "LAr" => :LAr,
     "CZT" => :CdZnTe,
-    "Si" => :Si
+    "Si" => :Si,
+    "Lead" => :Pb
 )
 
 

src/SolidStateDetectors.jl

@@ -67,7 +67,7 @@ export Grid, GridPoint
 export ElectricPotential, PointTypes, EffectiveChargeDensity, DielectricDistribution, WeightingPotential, ElectricField
 export apply_initial_state!
 export calculate_electric_potential!, calculate_weighting_potential!, calculate_electric_field!, calculate_drift_fields!
-export ElectricFieldChargeDriftModel, ADLChargeDriftModel
+export ElectricFieldChargeDriftModel, ADLChargeDriftModel, IsotropicChargeDriftModel
 export NoChargeTrappingModel, BoggsChargeTrappingModel
 export get_active_volume, is_depleted, estimate_depletion_voltage
 export calculate_stored_energy, calculate_mutual_capacitance, calculate_capacitance_matrix

test/Project.toml

@@ -1,6 +1,7 @@
 [deps]
 ArraysOfArrays = "65a8f2f4-9b39-5baf-92e2-a9cc46fdf018"
 Geant4 = "559df036-b7a0-42fd-85df-7d5dd9d70f44"
+InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 LegendDataTypes = "99e09c13-5545-5ee2-bfa2-77f358fb75d8"
 LegendHDF5IO = "c9265ca6-b027-5446-b1a4-febfa8dd10b0"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -17,5 +18,5 @@ TypedTables = "9d95f2ec-7b3d-5a63-8d20-e2491e220bb9"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [compat]
-LegendHDF5IO = "0.1"
 LegendDataTypes = "0.1.3"
+LegendHDF5IO = "0.1"

test/test_charge_drift_models.jl

@@ -3,7 +3,8 @@
 using Test
 
 using SolidStateDetectors
-using SolidStateDetectors: getVe, getVh, Vl, get_path_to_example_config_files
+using SolidStateDetectors: getVe, getVh, Vl, get_path_to_example_config_files, AbstractChargeDriftModel
+using InteractiveUtils
 using StaticArrays
 using LinearAlgebra
 using Unitful
@@ -43,6 +44,52 @@ end
     @test signalsum < T(1.99)
 end
 
+@timed_testset "Test completeness of charge drift models" begin
+    for c in InteractiveUtils.subtypes(AbstractChargeDriftModel)
+        if isstructtype(c)
+            @testset "$(c)" begin 
+                cdm = c{T}
+                fv = @SVector T[1,0,0]
+
+                # default constructor
+                @test cdm() isa c
+                # electron-drift velocity
+                @test hasmethod(getVe, Tuple{SVector{3,T}, cdm})
+                @test getVe(fv, cdm()) isa SVector
+                # hole-drift velocity
+                @test hasmethod(getVh, Tuple{SVector{3,T}, cdm})
+                @test getVh(fv, cdm()) isa SVector
+            end
+        end
+    end
+end
+
+
+@timed_testset "Test IsotropicChargeDriftModel" begin
+
+    @testset "Test constructors of IsotropicChargeDriftModel" begin
+        cdm0 = IsotropicChargeDriftModel{T}() # default charge drift model
+        @test cdm0.μ_e  == 0.1f0
+        @test cdm0.μ_h == 0.1f0
+
+        cdm1 = IsotropicChargeDriftModel{T}(1000u"cm^2/(V*s)", 1000u"cm^2/(V*s)")
+        @test cdm1 == cdm0
+
+        cdm2 = IsotropicChargeDriftModel{T}(0.1, 0.1)
+        @test cdm2 == cdm0
+
+        config_dict = Dict(
+            "model" => "IsotropicChargeDriftModel",
+            "mobilities" => Dict(
+                "e" => "1000cm^2/(V*s)",
+                "h" => "1000cm^2/(V*s)"
+            )
+        )
+        cdm3 = IsotropicChargeDriftModel{T}(config_dict)
+        @test cdm3 == cdm0
+    end
+end
+
 @timed_testset "Test ADLChargeDriftModel" begin
 
     config_dict = Dict(