Generic simple_calibration

ext/LegendSpecFitsRecipesBaseExt.jl

@@ -456,24 +456,25 @@ end
     end
 end
 
-@recipe function f(report::NamedTuple{(:h_calsimple, :h_uncal, :c, :fep_guess, :peakhists, :peakstats)}; cal=true)
+@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:200:3000, ["$i" for i in 0:200:3000])
+        xticks := (0:500:3000, ["$i" for i in 0:500:3000])
         ylims := (0.2, maximum(h.weights)*1.1)
-        fep_guess = 2614.5
+        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.fep_guess)
-        xticks := (0:3000:1.2*report.fep_guess, ["$i" for i in 0:3000:1.2*report.fep_guess])
+        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)
-        fep_guess = report.fep_guess
+        peak_guess = report.peak_guess
     end
     @series begin
         seriestype := :stepbins
@@ -483,10 +484,10 @@ end
     y_vline = 0.2:1:maximum(h.weights)*1.1
     @series begin
         seriestype := :line
-        label := "FEP Guess"
-        color := :red
+        label := "Peak estimate"#: $(round(peak_guess, digits = 1))"
+        color := :red2
         linewidth := 1.5
-        fill(fep_guess, length(y_vline)), y_vline
+        fill(peak_guess, length(y_vline)), y_vline
     end
 end
 
@@ -820,7 +821,8 @@ end
         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/chi2fit.jl

@@ -94,7 +94,13 @@ function chi2fit(f_fit::Function, x::AbstractVector{<:Union{Real,Measurement{<:R
             # gof 
             chi2min = res.objective
             dof = length(x) - length(v_chi2)
-            pvalue = ccdf(Chisq(dof), chi2min)
+            if dof == 0
+                pvalue = NaN
+                  @warn "The number of fit parameters is equivalent to number of data points --> dof = 0 ; p-value not meaningful, set to NaN"
+            else 
+                pvalue = ccdf(Chisq(dof), chi2min)
+            end 
+
             function get_y_pred_err(f_fit, x_val, x_err, pars)  # get final uncertainties for normalized residuals
                 dual_x = ForwardDiff.Dual{UncertTag}(x_val, x_err)
                 dual_y = f_fit(dual_x, pars...)

src/pseudo_prior.jl

@@ -33,6 +33,8 @@ function get_standard_pseudo_prior(h::Histogram, ps::NamedTuple{(:peak_pos, :pea
         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/specfit.jl

@@ -194,7 +194,11 @@ export fit_single_peak_th228
 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 @@ Get the FWHM of a peak from the fit parameters.
 # 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 @@ This function defines the gamma peakshape fit functions used in the calibration
 * 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 @@ function get_th228_fit_functions(; background_center::Union{Real,Nothing} = noth
         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 @@ function peakshape_components(fit_func::Symbol; background_center::Union{Real,No
             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)

test/Project.toml

@@ -8,6 +8,7 @@ Measurements = "eff96d63-e80a-5855-80a2-b1b0885c5ab7"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
+IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953"
 
 [compat]
 Distributions = "0.24, 0.25"
@@ -16,4 +17,5 @@ Interpolations = "0.15"
 LegendDataTypes = "0.1"
 Measurements = "2"
 StatsBase = "0.32, 0.33, 0.34"
-Unitful = "1"
+Unitful = "1"
+IntervalSets = "0.7"

test/runtests.jl

@@ -5,6 +5,7 @@ import Test
 Test.@testset "Package LegendSpecFits" begin
     #include("test_aqua.jl")
     include("test_specfit.jl")
+    include("test_simplecal.jl")
     include("test_fit_chisq.jl")
     include("test_singlefit.jl")
     include("test_docs.jl")

test/test_fit_chisq.jl

@@ -6,18 +6,24 @@ using Test
 @testset "fit_chisq" begin
     par_true = [5, 2]
     f_lin(x,p1,p2)  = p1 +  p2 * x 
-    x = [1,2,3,4,5,6,7,8,9,10] #.± ones(10)
+    x = [1,2,3,4,5,6,7,8,9,10]
     y = f_lin.(x,par_true...) .+ 0.5.*randn(10)
     @info "chisq fit without uncertainties on x and y "
     result, report       = chi2fit(1, x, y; uncertainty=true) 
-
-    x = [1,2,3,4,5,6,7,8,9,10] .± ones(10)
+    @test isapprox(result.par[1], par_true[1], atol = 0.2*par_true[1])
+    @test isapprox(result.par[2], par_true[2], atol = 0.2*par_true[2])
+
+    x = measurement.([1,2,3,4,5,6,7,8,9,10], ones(10))
     y = f_lin.(x,par_true...) .+ 0.5.*randn(10)
     @info "chisq fit with uncertainties on x and y"
     result, report       = chi2fit(1, x, y; uncertainty=true) 
-
-    x = [1,2,3,4,5,6,7,8,9,10] .± ones(10)
+    @test isapprox(result.par[1], par_true[1], atol = 0.2*par_true[1])
+    @test isapprox(result.par[2], par_true[2], atol = 0.2*par_true[2])
+
+    x = measurement.([1,2,3,4,5,6,7,8,9,10], ones(10))
     y = f_lin.(x,par_true...) .+ 0.5.*randn(10)
     @info "chisq fit with uncertainties on x and y"
     result, report       = chi2fit(1, x, y; pull_t = [(mean = par_true[1], std= 0.1),(mean = par_true[2],std= 0.1)], uncertainty=true) 
+    @test isapprox(result.par[1], par_true[1], atol = 0.2*par_true[1])
+    @test isapprox(result.par[2], par_true[2], atol = 0.2*par_true[2])
 end

test/test_simplecal.jl

@@ -0,0 +1,32 @@
+# This file is a part of LegendSpecFits.jl, licensed under the MIT License (MIT).
+using LegendSpecFits
+using Test
+using LegendDataTypes: fast_flatten
+using Measurements: value as mvalue
+using Distributions
+using StatsBase
+using Interpolations
+using Unitful 
+using IntervalSets
+include("test_utils.jl")
+
+@testset "simplecal" begin
+    # generate fake data - very simplified
+    energy_test, th228_lines = generate_mc_spectrum(200000)
+    m_cal_simple = 0.123u"keV"
+    e_uncal = energy_test ./ m_cal_simple
+
+    # binning and peak-finder settings that should work for this fake data set 
+    window_sizes =  vcat([(25.0u"keV",25.0u"keV") for _ in 1:6], (30.0u"keV",30.0u"keV"))
+    quantile_perc=NaN
+    peakfinder_σ = 2.0
+    peakfinder_threshold = 6.0
+    peak_quantile = 0.7..1.0
+    bin_quantile =  0.05..0.5
+    quantile_perc = NaN 
+    kwargs = (peakfinder_σ=peakfinder_σ, peakfinder_threshold=peakfinder_threshold, peak_quantile=peak_quantile, bin_quantile=bin_quantile, quantile_perc = quantile_perc)
+
+    # simple calibration
+    result_simple, report_simple = simple_calibration(e_uncal, th228_lines, window_sizes,; calib_type= :th228, kwargs...);
+    @test isapprox(result_simple.c, m_cal_simple, atol = 0.05.*(m_cal_simple))
+end

test/test_specfit.jl

@@ -15,7 +15,7 @@ include("test_utils.jl")
 
     # simple calibration fit 
     window_sizes =  vcat([(25.0u"keV",25.0u"keV") for _ in 1:6], (30.0u"keV",30.0u"keV"))
-    result_simple, report_simple = simple_calibration(ustrip.(energy_test), th228_lines, window_sizes, n_bins=10000,; calib_type=:th228, quantile_perc=0.995)
+    result_simple, report_simple = simple_calibration(ustrip.(energy_test), th228_lines, window_sizes,; calib_type=:th228, quantile_perc=0.995)
 
     # fit 
     result, report = fit_peaks(result_simple.peakhists, result_simple.peakstats, ustrip.(th228_lines),; uncertainty=true,calib_type = :th228);