Add RadfordImpurityDensity model to the SSD extension

Project.toml

@@ -11,6 +11,7 @@ IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 LRUCache = "8ac3fa9e-de4c-5943-b1dc-09c6b5f20637"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+LsqFit = "2fda8390-95c7-5789-9bda-21331edee243"
 MIMEs = "6c6e2e6c-3030-632d-7369-2d6c69616d65"
 Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a"
 Measurements = "eff96d63-e80a-5855-80a2-b1b0885c5ab7"
@@ -49,6 +50,7 @@ Glob = "1.3"
 IntervalSets = "0.6, 0.7"
 JSON = "0.21.2, 1"
 LRUCache = "1.5"
+LsqFit = "0.13, 0.14, 0.15"
 LegendDataTypes = "0.1.13"
 LegendHDF5IO = "0.1.14"
 LinearAlgebra = "<0.0.1, 1"

ext/LegendDataManagementSolidStateDetectorsExt.jl

@@ -6,9 +6,35 @@
 using LegendDataManagement
 using Unitful
 using PropDicts
+using LsqFit
 
 const _SSDDefaultNumtype = Float32
 
+struct RadfordImpurityDensity{T} <: SolidStateDetectors.AbstractImpurityDensity{T}
+    # a + b*z + c*exp((z-L)/tau) -> needs at least 4 points
+    a::T
+    b::T 
+    c::T 
+    tau::T
+    L::T
+    det_z0::T
+end
+
+function SolidStateDetectors.get_impurity_density(
+    idm::RadfordImpurityDensity, pt::SolidStateDetectors.AbstractCoordinatePoint{T}
+    )::T where {T}
+    cpt = CartesianPoint(pt)
+    z = cpt[3]
+
+    # the function parameters are in crystal axis coordinates i.e. z = 0 is seed end, z = L crystal length 
+    # -> convert to detector coordiantes where z = 0 corresponds to p+ contact i.e. z -> det_z0 - z
+    idm.a .+ idm.b * (idm.det_z0 .- z) .+ idm.c * exp.((idm.det_z0 .- z .- idm.L)/idm.tau)
+end
+
+function SolidStateDetectors.ImpurityDensity(T::DataType, t::Val{:radford}, dict::AbstractDict, input_units::NamedTuple)
+    RadfordImpurityDensity{T}(dict["parameters"]..., )
+end
+
 
 """
     SolidStateDetector[{T<:AbstractFloat}](data::LegendData, detector::DetectorIdLike)
@@ -19,34 +45,34 @@
 `SolidStateDetector` can be constructed from LEGEND metadata using the
 methods above.
 """
-function SolidStateDetectors.SolidStateDetector(data::LegendData, detector::DetectorIdLike)
-    SolidStateDetectors.SolidStateDetector{_SSDDefaultNumtype}(data, detector)
+function SolidStateDetectors.SolidStateDetector(data::LegendData, detector::DetectorIdLike; kwargs...)
+    SolidStateDetectors.SolidStateDetector{_SSDDefaultNumtype}(data, detector; kwargs...)
 end
 
-function SolidStateDetectors.SolidStateDetector{T}(data::LegendData, detector::DetectorIdLike) where {T<:AbstractFloat}
+function SolidStateDetectors.SolidStateDetector{T}(data::LegendData, detector::DetectorIdLike; kwargs...) where {T<:AbstractFloat}
     detector_props = getproperty(data.metadata.hardware.detectors.germanium.diodes, Symbol(detector))
     xtal_props = getproperty(data.metadata.hardware.detectors.germanium.crystals, Symbol(string(detector)[1:end-1]))
-    SolidStateDetector{T}(LegendData, detector_props, xtal_props)
+    SolidStateDetector{T}(LegendData, detector_props, xtal_props; kwargs...)
 end
 
-function SolidStateDetectors.SolidStateDetector{T}(::Type{LegendData}, filename::String) where {T<:AbstractFloat}
-    SolidStateDetector{T}(LegendData, readprops(filename, subst_pathvar = false, subst_env = false, trim_null = false))
+function SolidStateDetectors.SolidStateDetector{T}(::Type{LegendData}, filename::String; kwargs...) where {T<:AbstractFloat}
+    SolidStateDetector{T}(LegendData, readprops(filename, subst_pathvar = false, subst_env = false, trim_null = false); kwargs...)
 end
 
-function SolidStateDetectors.SolidStateDetector(::Type{LegendData}, filename::String)
-    SolidStateDetector{_SSDDefaultNumtype}(LegendData, filename)
+function SolidStateDetectors.SolidStateDetector(::Type{LegendData}, filename::String; kwargs...)
+    SolidStateDetector{_SSDDefaultNumtype}(LegendData, filename; kwargs...)
 end
 
-function SolidStateDetectors.SolidStateDetector(::Type{LegendData}, meta::AbstractDict)
-    SolidStateDetectors.SolidStateDetector{_SSDDefaultNumtype}(LegendData, meta)
+function SolidStateDetectors.SolidStateDetector(::Type{LegendData}, meta::AbstractDict; kwargs...)
+    SolidStateDetectors.SolidStateDetector{_SSDDefaultNumtype}(LegendData, meta; kwargs...)
 end
 
-function SolidStateDetectors.SolidStateDetector{T}(::Type{LegendData}, meta::AbstractDict) where {T<:AbstractFloat}
-    SolidStateDetectors.SolidStateDetector{T}(LegendData, convert(PropDict, meta), LegendDataManagement.NoSuchPropsDBEntry("",[]))
+function SolidStateDetectors.SolidStateDetector{T}(::Type{LegendData}, meta::AbstractDict; kwargs...) where {T<:AbstractFloat}
+    SolidStateDetectors.SolidStateDetector{T}(LegendData, convert(PropDict, meta), LegendDataManagement.NoSuchPropsDBEntry("",[]); kwargs...)
 end
 
-function SolidStateDetectors.SolidStateDetector{T}(::Type{LegendData}, meta::PropDict, xtal_meta::Union{PropDict, LegendDataManagement.NoSuchPropsDBEntry}) where {T<:AbstractFloat}
-    config_dict = create_SSD_config_dict_from_LEGEND_metadata(meta, xtal_meta)
+function SolidStateDetectors.SolidStateDetector{T}(::Type{LegendData}, meta::PropDict, xtal_meta::Union{PropDict, LegendDataManagement.NoSuchPropsDBEntry}; kwargs...) where {T<:AbstractFloat}
+    config_dict = create_SSD_config_dict_from_LEGEND_metadata(meta, xtal_meta; kwargs...)
     return SolidStateDetector{T}(config_dict, SolidStateDetectors.construct_units(config_dict))
 end
 
@@ -59,39 +85,39 @@
 `Simulation` can be constructed from LEGEND metadata using the
 methods above.
 """
-function SolidStateDetectors.Simulation(data::LegendData, detector::DetectorIdLike)
-    SolidStateDetectors.Simulation{_SSDDefaultNumtype}(data, detector)
+function SolidStateDetectors.Simulation(data::LegendData, detector::DetectorIdLike; kwargs...)
+    SolidStateDetectors.Simulation{_SSDDefaultNumtype}(data, detector; kwargs...)
 end
 
-function SolidStateDetectors.Simulation{T}(data::LegendData, detector::DetectorIdLike) where {T<:AbstractFloat}
+function SolidStateDetectors.Simulation{T}(data::LegendData, detector::DetectorIdLike; kwargs...) where {T<:AbstractFloat}
     detector_props = getproperty(data.metadata.hardware.detectors.germanium.diodes, Symbol(detector))
     xtal_props = getproperty(data.metadata.hardware.detectors.germanium.crystals, Symbol(string(detector)[1:end-1]))
-    Simulation{T}(LegendData, detector_props, xtal_props)
+    Simulation{T}(LegendData, detector_props, xtal_props; kwargs...)
 end
 
-function SolidStateDetectors.Simulation{T}(::Type{LegendData}, filename::String) where {T<:AbstractFloat}
-    Simulation{T}(LegendData, readprops(filename, subst_pathvar = false, subst_env = false, trim_null = false))
+function SolidStateDetectors.Simulation{T}(::Type{LegendData}, filename::String; kwargs...) where {T<:AbstractFloat}
+    Simulation{T}(LegendData, readprops(filename, subst_pathvar = false, subst_env = false, trim_null = false); kwargs...)
 end
 
-function SolidStateDetectors.Simulation(::Type{LegendData}, filename::String)
-    Simulation{_SSDDefaultNumtype}(LegendData, filename)
+function SolidStateDetectors.Simulation(::Type{LegendData}, filename::String; kwargs...)
+    Simulation{_SSDDefaultNumtype}(LegendData, filename; kwargs...)
 end
 
-function SolidStateDetectors.Simulation(::Type{LegendData}, meta::AbstractDict)
-    SolidStateDetectors.Simulation{_SSDDefaultNumtype}(LegendData, meta)
+function SolidStateDetectors.Simulation(::Type{LegendData}, meta::AbstractDict; kwargs...)
+    SolidStateDetectors.Simulation{_SSDDefaultNumtype}(LegendData, meta; kwargs...)
 end
 
-function SolidStateDetectors.Simulation{T}(::Type{LegendData}, meta::AbstractDict) where {T<:AbstractFloat}
-    SolidStateDetectors.Simulation{T}(LegendData, convert(PropDict, meta), LegendDataManagement.NoSuchPropsDBEntry("", []))
+function SolidStateDetectors.Simulation{T}(::Type{LegendData}, meta::AbstractDict; kwargs...) where {T<:AbstractFloat}
+    SolidStateDetectors.Simulation{T}(LegendData, convert(PropDict, meta), LegendDataManagement.NoSuchPropsDBEntry("", []); kwargs...)
 end
 
-function SolidStateDetectors.Simulation{T}(::Type{LegendData}, meta::PropDict, xtal_meta::Union{PropDict, LegendDataManagement.NoSuchPropsDBEntry}) where {T<:AbstractFloat}
-    config_dict = create_SSD_config_dict_from_LEGEND_metadata(meta, xtal_meta)
+function SolidStateDetectors.Simulation{T}(::Type{LegendData}, meta::PropDict, xtal_meta::Union{PropDict, LegendDataManagement.NoSuchPropsDBEntry}; kwargs...) where {T<:AbstractFloat}
+    config_dict = create_SSD_config_dict_from_LEGEND_metadata(meta, xtal_meta; kwargs...)
     return Simulation{T}(config_dict)
 end
 
 
-function create_SSD_config_dict_from_LEGEND_metadata(meta::PropDict, xtal_meta::X; dicttype = Dict{String,Any}) where {X <: Union{PropDict, LegendDataManagement.NoSuchPropsDBEntry}}
+function create_SSD_config_dict_from_LEGEND_metadata(meta::PropDict, xtal_meta::X; dicttype = Dict{String,Any}, crystal_impurity::Bool = false) where {X <: Union{PropDict, LegendDataManagement.NoSuchPropsDBEntry}}
 
     # Not all possible configurations are yet implemented!
     # https://github.com/legend-exp/legend-metadata/blob/main/hardware/detectors/detector-metadata_1.pdf
@@ -527,11 +553,43 @@
         mantle_contact_parts
     end
 
-
-    config_dict["detectors"][1]["semiconductor"]["impurity_density"] = dicttype(
-        "name" => "constant", 
-        "value" => "-1e9cm^-3"
-    )
+    ### IMPURITY DENSITY ###
+
+    if X != PropDict
+        @warn "No crystal metadata found for detector $(meta.name)"
+    end
+
+    if X == PropDict && !haskey(xtal_meta, :impurity_measurements)
+        @warn "No information regarding impurity density for $(xtal_meta.name)"
+    end
+
+    if !crystal_impurity
+        @warn """
+        Reading the impurity density from the crystal metadata will be ignored.
+        Set `crystal_impurity=true` to load the impurity density from the crystal metadata
+        """
+    end
+
+    config_dict["detectors"][1]["semiconductor"]["impurity_density"] = if crystal_impurity && X == PropDict && haskey(xtal_meta, :impurity_measurements) && haskey(xtal_meta.impurity_measurements, :value_in_1e9e_cm3) && !isempty(xtal_meta.impurity_measurements.value_in_1e9e_cm3)
+        @info "Reading impurity density values from crystal metadata $(xtal_meta.order)$(xtal_meta.name)"
+        # Fit the impurity measurement data to a Radford model
+        @. fit_model(z, p) = p[1] + p[2]*z + p[3]*exp((z-p[5])/p[4])
+        pos =  xtal_meta.impurity_measurements.distance_from_seed_end_mm * 1e-3  # units: m
+        val = -xtal_meta.impurity_measurements.value_in_1e9e_cm3 * 1e15          # units: e/m^-3 (p-type => negative values)
+        fit_result = curve_fit(fit_model, pos, val, [-1e15, -1e15, -1e15, 1., 1.])
+        dicttype(
+            "name" => "radford", 
+            "parameters" => vcat(fit_result.param..., xtal_meta.slices[Symbol(meta.name[end])].detector_offset_in_mm * 1e-3)
+        )
+    else
+        # default impurity density for cases without crystal metadata
+        default_impurity_value = "-0.9e10cm^-3"
+        @info "Set impurity density to constant default value of $(default_impurity_value)"
+        dicttype(
+            "name" => "constant", 
+            "value" => default_impurity_value
+        )
+    end
 
     # evaluate "include" statements - needed for the charge drift model
     SolidStateDetectors.scan_and_merge_included_json_files!(config_dict, "")

test/test_ext_ssd.jl

@@ -16,11 +16,15 @@ include("testing_utils.jl")
         @testset "$(detname)" begin
             det = SolidStateDetector{Float64}(l200, detname) 
             @test det isa SolidStateDetector
-            sim = Simulation{Float64}(l200, detname)
+            sim = Simulation{Float64}(l200, detname, crystal_impurity = true)
             @test sim isa Simulation
 
-            # Compare active volume from SSD to active volume from LegendDataManagement
             SolidStateDetectors.apply_initial_state!(sim, ElectricPotential, Grid(sim, max_tick_distance = 0.1u"mm"))
+
+            # Check that all crystals are p-type
+            @test all(sim.q_eff_imp.data .<= 0)
+
+            # Compare active volume from SSD to active volume from LegendDataManagement
             active_volume_ssd = SolidStateDetectors.get_active_volume(sim.point_types)
             active_volume_ldm = LegendDataManagement.get_active_volume(l200.metadata.hardware.detectors.germanium.diodes[Symbol(detname)], 0.0)
             @test isapprox(active_volume_ssd, active_volume_ldm, rtol = 0.01)