Add and apply RadfordImpurityDensity in 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"
@@ -30,8 +31,8 @@ UUIDs = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [weakdeps]
-LegendHDF5IO = "c9265ca6-b027-5446-b1a4-febfa8dd10b0"
 LegendDataTypes = "99e09c13-5545-5ee2-bfa2-77f358fb75d8"
+LegendHDF5IO = "c9265ca6-b027-5446-b1a4-febfa8dd10b0"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 SolidStateDetectors = "71e43887-2bd9-5f77-aebd-47f656f0a3f0"
@@ -48,10 +49,11 @@ Format = "1.2"
 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"
-LRUCache = "1.5"
 MIMEs = "0.1"
 Markdown = "<0.0.1, 1"
 Measurements = "2.2.1"

ext/LegendDataManagementSolidStateDetectorsExt.jl

@@ -6,9 +6,36 @@ using SolidStateDetectors
 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)
@@ -527,11 +554,27 @@ function create_SSD_config_dict_from_LEGEND_metadata(meta::PropDict, xtal_meta::
         mantle_contact_parts
     end
 
-
-    config_dict["detectors"][1]["semiconductor"]["impurity_density"] = dicttype(
-        "name" => "constant", 
-        "value" => "-1e9cm^-3"
-    )
+
+    if X == PropDict
+
+        if !haskey(xtal_meta, :impurity_measurements)
+            @warn "No information regarding impurity density for $(xtal_meta.name)"
+        end
+        # TODO: check units of impurity density against the values in the config file
+        # 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                 # units: e/m^-3
+        fit_result = curve_fit(fit_model, pos, val, ones(Float64,5))
+
+        config_dict["detectors"][1]["semiconductor"]["impurity_density"] =  dicttype(
+            "name" => "radford", 
+            "parameters" => vcat(fit_result.param..., xtal_meta.slices[Symbol(meta.name[end])].detector_offset_in_mm / 1000)
+        )
+    else
+        @warn "No crystal metadata found for detector $(meta.name)"
+        # TODO: Implement a default impurity density for cases without crystal metadata
+    end
 
     # evaluate "include" statements - needed for the charge drift model
     SolidStateDetectors.scan_and_merge_included_json_files!(config_dict, "")