Generic simple_calibration

ext/LegendSpecFitsRecipesBaseExt.jl

@@ -486,24 +486,25 @@
     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
@@ -513,10 +514,10 @@
     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
 
@@ -850,7 +851,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/chi2fit.jl

@@ -94,7 +94,13 @@
             # gof 
             chi2min = res.objective
             dof = length(x) - length(v_chi2)
-            pvalue = ccdf(Chisq(dof), chi2min)
+            pvalue = if iszero(dof)
+                @warn "The number of fit parameters is equivalent to number of data points --> dof = 0 ; p-value not meaningful, set to NaN"
+                NaN
+            else 
+                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/filter_optimization.jl

@@ -132,11 +132,16 @@
     modes = Vector{Float64}(undef, length(e_grid_ft))
     fts_success  = Bool.(zeros(length(e_grid_ft)))
 
-    Threads.@threads for f in eachindex(e_grid_ft)
+    # Threads.@threads 
+    for f in eachindex(e_grid_ft)
         # get ft
         ft = e_grid_ft[f]
 
         # get e values for this ft
+        if all(ismissing.(e_grid[f,:]))
+            @debug "All values are missing, skipping"
+            continue
+        end
         e_ft = Array{Float64}(flatview(e_grid)[f, :])
         e_ft = e_ft[isfinite.(e_ft)]
 

src/fit_calibration.jl

@@ -51,8 +51,8 @@
     result_fit = merge(result_fit, (par = par_unit,))
 
     # built function in string 
-    func = join(["$(mvalue(par[i]))$e_unit .* ($(e_expression)).^$(i-1)" for i in eachindex(par)], " .+ ")
-    func_err = join(["($(par[i]))$e_unit .* ($(e_expression)).^$(i-1)" for i in eachindex(par)], " .+ ")
+    func = join(["$(mvalue(par[i]))$e_unit * ($(e_expression))^$(i-1)" for i in eachindex(par)], " + ")
+    func_err = join(["($(par[i]))$e_unit * ($(e_expression))^$(i-1)" for i in eachindex(par)], " + ")
 
     result = merge(result_fit, (func = func, func_err = func_err, µ = μ, peaks = peaks))
     report = merge(report_fit, (e_unit = e_unit, par = result.par, type = :cal))

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