Generic simple_calibration

ext/LegendSpecFitsRecipesBaseExt.jl

@@ -490,6 +490,41 @@
     end
 end
 
+@recipe function f(report::NamedTuple{(:h_calsimple, :h_uncal, :c, :peak_guess, :peakhists, :peakstats)}; cal=true)
+    ylabel := "Counts"
+    legend := :topright
+    yscale := :log10
+    fill := false
+    if cal
+        h = LinearAlgebra.normalize(report.h_calsimple, mode = :density)
+        xlabel := "Energy (keV)"
+        xlims := (0, 3000)
+        xticks := (0:500:3000, ["$i" for i in 0:500:3000])
+        ylims := (0.2, maximum(h.weights)*1.1)
+        peak_guess = ustrip(report.c * report.peak_guess)
+    else
+        h = LinearAlgebra.normalize(report.h_uncal, mode = :density)
+        xlabel := "Energy (ADC)"
+        xlims := (0, 1.2*report.peak_guess)
+        xticks := (0:3000:1.2*report.peak_guess, ["$i" for i in 0:3000:1.2*report.peak_guess])
+        ylims := (0.2, maximum(h.weights)*1.1)
+        peak_guess = report.peak_guess
+    end
+    @series begin
+        seriestype := :stepbins
+        label := "Energy"
+        h
+    end
+    y_vline = 0.2:1:maximum(h.weights)*1.1
+    @series begin
+        seriestype := :line
+        label := "Peak guess"#: $(round(peak_guess, digits = 1))"
+        color := :red2
+        linewidth := 1.5
+        fill(peak_guess, length(y_vline)), y_vline
+    end
+end
+
 @recipe function f(report::NamedTuple{(:peakpos, :peakpos_cal, :h_uncal, :h_calsimple)}; cal=true)
     legend := :topright
     size := (1000, 600)
@@ -820,7 +855,8 @@
         if plot_ribbon
             ribbon := uncertainty.(report.f_fit.(0:1:1.2*value(maximum(report.x))))
         end
-        0:1:1.2*value(maximum(report.x)), value.(report.f_fit.(0:1:1.2*value(maximum(report.x))))
+        xstep = value(maximum(report.x))/100
+        0:xstep:1.2*value(maximum(report.x)), value.(report.f_fit.(0:xstep:1.2*value(maximum(report.x))))
     end
     @series begin
         seriestype := :scatter

src/pseudo_prior.jl

@@ -33,6 +33,8 @@
         NamedTupleDist(; μ, σ, n, skew_fraction, skew_width, background, step_amplitude, background_slope)
     elseif fit_func == :gamma_bckExp
         NamedTupleDist(; μ, σ, n, skew_fraction, skew_width, background, step_amplitude, background_exp)
+    elseif fit_func == :gamma_minimal
+        NamedTupleDist(; μ, σ, n, background)
     else
         throw(ArgumentError("Unknown fit function: $fit_func"))
     end

src/simple_calibration.jl

@@ -18,23 +18,25 @@
 function simple_calibration end
 export simple_calibration
 
-function simple_calibration(e_uncal::Vector{<:Real}, th228_lines::Vector{<:Unitful.Energy{<:Real}}, window_sizes::Vector{<:Tuple{Unitful.Energy{<:Real}, Unitful.Energy{<:Real}}},; kwargs...)
+function simple_calibration(e_uncal::Vector{<:Real}, gamma_lines::Vector{<:Unitful.Energy{<:Real}}, window_sizes::Vector{<:Tuple{Unitful.Energy{<:Real}, Unitful.Energy{<:Real}}},; kwargs...)
     # remove calib type from kwargs
     @assert haskey(kwargs, :calib_type) "Calibration type not specified"
     calib_type = kwargs[:calib_type]
     # remove :calib_type from kwargs
     kwargs = pairs(NamedTuple(filter(k -> !(:calib_type in k), kwargs)))
     if calib_type == :th228
         @debug "Use simple calibration for Th228 lines"
-        return simple_calibration_th228(e_uncal, th228_lines, window_sizes,; kwargs...)
+        return simple_calibration_th228(e_uncal, gamma_lines, window_sizes,; kwargs...)
+    elseif calib_type == :gamma
+        return simple_calibration_gamma(e_uncal, gamma_lines, window_sizes,; kwargs...)
     else
         error("Calibration type not supported")
     end
 end
-simple_calibration(e_uncal::Vector{<:Real}, th228_lines::Vector{<:Unitful.Energy{<:Real}}, left_window_sizes::Vector{<:Unitful.Energy{<:Real}}, right_window_sizes::Vector{<:Unitful.Energy{<:Real}}; kwargs...) = simple_calibration(e_uncal, th228_lines, [(l,r) for (l,r) in zip(left_window_sizes, right_window_sizes)],; kwargs...)
+simple_calibration(e_uncal::Vector{<:Real}, gamma_lines::Vector{<:Unitful.Energy{<:Real}}, left_window_sizes::Vector{<:Unitful.Energy{<:Real}}, right_window_sizes::Vector{<:Unitful.Energy{<:Real}}; kwargs...) = simple_calibration(e_uncal, gamma_lines, [(l,r) for (l,r) in zip(left_window_sizes, right_window_sizes)],; kwargs...)
 
 
-function simple_calibration_th228(e_uncal::Vector{<:Real}, th228_lines::Vector{<:Unitful.Energy{<:Real}}, window_sizes::Vector{<:Tuple{Unitful.Energy{<:Real}, Unitful.Energy{<:Real}}},; n_bins::Int=15000, quantile_perc::Float64=NaN, binning_peak_window::Unitful.Energy{<:Real}=10.0u"keV")
+function simple_calibration_th228(e_uncal::Vector{<:Real}, th228_lines::Vector{<:Unitful.Energy{<:Real}}, window_sizes::Vector{<:Tuple{Unitful.Energy{<:Real}, Unitful.Energy{<:Real}}},; n_bins::Int=15000, quantile_perc::Float64=NaN, binning_peak_window::Unitful.Energy{<:Real}=10.0u"keV", e_unit::Unitful.EnergyUnits=u"keV")
     # initial binning
     bin_width = get_friedman_diaconis_bin_width(filter(in(quantile(e_uncal, 0.05)..quantile(e_uncal, 0.5)), e_uncal))
     # create initial peak search histogram
@@ -55,7 +57,6 @@
     # get calibration constant for simple calibration
     c = 2614.5*u"keV" / fep_guess
     e_simple = e_uncal .* c
-    e_unit = u"keV"
     # get peakhists and peakstats
     peakhists, peakstats, h_calsimple, bin_widths = get_peakhists_th228(e_simple, th228_lines, window_sizes; binning_peak_window=binning_peak_window, e_unit=e_unit)
 
@@ -81,7 +82,68 @@
     return result, report
 end
 
-
+"""
+    simple_calibration_gamma(e_uncal::Vector{<:Real}, gamma_lines::Vector{<:Unitful.Energy{<:Real}}, window_sizes::Vector{<:Tuple{Unitful.Energy{<:Real}, Unitful.Energy{<:Real}}},; n_bins::Int=15000, quantile_perc::Float64=NaN, binning_peak_window::Unitful.Energy{<:Real}=10.0u"keV", peak_quantile::ClosedInterval{<:Real} = 0.5..1.0)
+Perform a simple calibration for the uncalibrated energy array `e_uncal` for any gamma source 
+* 1. Run a peak-finding algorithm. Only the energy window that is defined by `peak_quantile` is considered for the peak search.
+* 2. Estimate the calibration constant `c` by dividing the energy of the `gamma_lines` by the energy of the most prominent peak from peak-finder 
+* Return peakstats and peakhists of gamma lines for further calibration 
+INPUTS:
+* `e_uncal`: uncalibrated energy array
+* `gamma_lines`: array of gamma lines
+* `window_sizes`: array of tuples with left and right window sizes around the gamma lines
+KEYWORD ARGUMENTS:
+* `n_bins`: number of bins for the histogram (not used, legacy from simple_calibration)
+* `quantile_perc`: quantile percentage for the most prominent peak  (not used, legacy from simple_calibration)
+* `binning_peak_window`: window size for the peak search histogram (default: 10 keV)
+* `peak_quantile`: quantile range for the peak search (default: 0.5..1.0)
+"""
+function simple_calibration_gamma(e_uncal::Vector{<:Real}, gamma_lines::Vector{<:Unitful.Energy{<:Real}}, window_sizes::Vector{<:Tuple{Unitful.Energy{<:Real}, Unitful.Energy{<:Real}}},; n_bins::Int=15000, quantile_perc::Float64=NaN, binning_peak_window::Unitful.Energy{<:Real}=10.0u"keV", peak_quantile::ClosedInterval{<:Real} = 0.5..1.0,  e_unit::Unitful.EnergyUnits=u"keV")
+    e_min = quantile(e_uncal, leftendpoint(peak_quantile))
+    e_max = quantile(e_uncal, rightendpoint(peak_quantile))
+    # initial binning
+    bin_width = get_friedman_diaconis_bin_width(filter(in(e_min..e_max), e_uncal))
+    # create initial peak search histogram
+    h_uncal = fit(Histogram, e_uncal, 0:bin_width:maximum(e_uncal))
+    peak_guess, peak_idx = if isnan(quantile_perc)
+        h_peaksearch = fit(Histogram, e_uncal, e_min:bin_width:e_max)
+        # search all possible peak candidates
+        h_decon, peakpos = RadiationSpectra.peakfinder(h_peaksearch, σ=2.0, backgroundRemove=true, threshold=10)
+        # find the most prominent peak in the deconvoluted histogram
+        peakpos_idxs = StatsBase.binindex.(Ref(h_decon), peakpos)
+        cts_peakpos = h_decon.weights[peakpos_idxs]
+        peakpos[argmax(cts_peakpos)], argmax(cts_peakpos)
+    else
+        quantile(e_uncal, quantile_perc), length(gamma_lines)
+    end
+    @info "Identified most prominent peak (peak $(peak_idx), $(gamma_lines[peak_idx])). Peak guess: $peak_guess"
+    # get calibration constant for simple calibration
+    c = gamma_lines[peak_idx] / peak_guess
+    e_simple = e_uncal .* c
+    # get peakhists and peakstats
+    peakhists, peakstats, h_calsimple, bin_widths = get_peakhists_th228(e_simple, gamma_lines, window_sizes; binning_peak_window=binning_peak_window, e_unit=e_unit)
+
+    result = (
+        h_calsimple = h_calsimple, 
+        h_uncal = h_uncal, 
+        c = c,
+        unit = e_unit,
+        bin_width = median(bin_widths),
+        peak_guess = peak_guess,
+        peakbinwidths = bin_widths,
+        peakhists = peakhists,
+        peakstats = peakstats
+        )
+    report = (
+        h_calsimple = result.h_calsimple,
+        h_uncal = result.h_uncal,
+        c = result.c,
+        peak_guess = result.peak_guess,
+        peakhists = result.peakhists,
+        peakstats = result.peakstats
+    )
+    return result, report
+end
 """
     get_peakhists_th228(e::Array, th228_lines::Array, window_sizes::Array, e_unit::String="keV", proxy_binning_peak::Float64=2103.5, proxy_binning_peak_window::Float64=10.0)
 

src/specfit.jl

@@ -194,7 +194,11 @@
 calculate centroid of gamma peak from fit parameters
 """
 function peak_centroid(v::NamedTuple)
-    centroid = v.μ - v.skew_fraction * (v.µ * v.skew_width)
+    if haskey(v, :skew_fraction)
+        centroid = v.μ - v.skew_fraction * (v.µ * v.skew_width)
+    else
+        centroid = v.μ
+    end
     if haskey(v, :skew_fraction_highE)
         centroid += v.skew_fraction_highE * (v.µ * v.skew_width_highE)
     end
@@ -208,30 +212,35 @@
 # Returns
     * `fwhm`: the FWHM of the peak
 """
-function estimate_fwhm(v::NamedTuple)
-
-    f_sigWithTail = Base.Fix2(get_th228_fit_functions().gamma_sigWithTail,v)
-    try
-        e_low, e_high = v.skew_fraction <= 0.5 ? (v.μ - v.σ, v.μ + v.σ) : (v.μ * (1 - v.skew_width), v.μ * (1 + v.skew_width))
-
-        max_sig = -Inf
-        for e in e_low:0.001:e_high
-            fe = f_sigWithTail(e)
-            if fe > max_sig
-                max_sig = fe
-            else
-                # if the maximum is reached,
-                # no need to further continue
-                break
+function estimate_fwhm(v::NamedTuple) 
+    if haskey(v, :skew_fraction)
+        f_sigWithTail = Base.Fix2(get_th228_fit_functions().gamma_sigWithTail,v)
+        try
+            e_low, e_high = v.skew_fraction <= 0.5 ? (v.μ - v.σ, v.μ + v.σ) : (v.μ * (1 - v.skew_width), v.μ * (1 + v.skew_width))
+
+            max_sig = -Inf
+            for e in e_low:0.001:e_high
+                fe = f_sigWithTail(e)
+                if fe > max_sig
+                    max_sig = fe
+                else
+                    # if the maximum is reached,
+                    # no need to further continue
+                    break
+                end
             end
+            half_max_sig = max_sig/2
+
+            tmp = x -> f_sigWithTail(x) - half_max_sig
+            roots_low = find_zero(tmp, e_low, maxiter=100)
+            roots_high = find_zero(tmp, e_high, maxiter=100)
+            return roots_high - roots_low
+        catch
+            return NaN
         end
-        half_max_sig = max_sig/2
-
-        tmp = x -> f_sigWithTail(x) - half_max_sig
-        roots_low = find_zero(tmp, e_low, maxiter=100)
-        roots_high = find_zero(tmp, e_high, maxiter=100)
-        return roots_high - roots_low
-    catch
+    elseif haskey(v, :σ)
+        return 2 * sqrt(2 * log(2)) * v.σ
+    else
         return NaN
     end
 end

src/specfit_functions.jl

@@ -8,6 +8,7 @@
 * gamma_bckExp: default gamma peakshape + exponential background 
 * gamma_bckFlat: default gamma peakshape - step background (only flat component!)
 * gamma_tails_bckFlat: default gamma peakshape + high-energy tail - step background (only flat component!)
+* gamma_minimal: only Gaussian signal and flat background 
 """
 function get_th228_fit_functions(; background_center::Union{Real,Nothing} = nothing)
     merge( 
@@ -16,6 +17,7 @@
         gamma_sigWithTail = (x, v) -> signal_peakshape(x, v.μ, v.σ, v.n, v.skew_fraction) + lowEtail_peakshape(x, v.μ, v.σ, v.n, v.skew_fraction, v.skew_width),
         gamma_bckFlat = (x, v) -> gamma_peakshape(x, v.μ, v.σ, v.n, 0.0, v.skew_fraction, v.skew_width, v.background),
         gamma_tails_bckFlat = (x, v) -> gamma_peakshape(x, v.μ, v.σ, v.n, 0.0, v.skew_fraction, v.skew_width , v.background; skew_fraction_highE = v.skew_fraction_highE, skew_width_highE = v.skew_width_highE),
+        gamma_minimal = (x, v) -> gamma_peakshape(x, v.μ, v.σ, v.n, 0.0, 0.0, 0.0, v.background),
         ),
     if isnothing(background_center)
         (gamma_bckSlope = (x, v) -> gamma_peakshape(x, v.μ, v.σ, v.n, v.step_amplitude, v.skew_fraction, v.skew_width, v.background; background_slope =  v.background_slope, background_center = v.μ),
@@ -61,6 +63,10 @@
             f_lowEtail = (x, v) -> lowEtail_peakshape(x, v.μ, v.σ, v.n, v.skew_fraction, v.skew_width),
             f_highEtail = (x, v) -> highEtail_peakshape(x, v.μ, v.σ, v.n, v.skew_fraction_highE, v.skew_width_highE), 
             f_bck = (x, v) -> background_peakshape(x, v.μ, v.σ, 0.0, v.background))
+    elseif fit_func == :gamma_minimal
+            funcs = (f_sig = (x, v) -> signal_peakshape(x, v.μ, v.σ, v.n, 0.0),
+            f_lowEtail = (x, v) -> lowEtail_peakshape(x, v.μ, v.σ, v.n, 0.0, 0.0),
+            f_bck = (x, v) -> background_peakshape(x, v.μ, v.σ, 0.0, v.background))
     end 
     labels = (f_sig = "Signal", f_lowEtail = "Low-energy tail", f_bck = "Background", f_highEtail = "High-energy tail")
     colors = (f_sig = :orangered1, f_lowEtail = :orange, f_bck = :dodgerblue2, f_highEtail = :forestgreen)