From fbf6c763c7a56eeccedf31bc883c5de80c39f889 Mon Sep 17 00:00:00 2001 From: LisaSchlueter Date: Tue, 12 Mar 2024 22:58:22 +0100 Subject: [PATCH 01/68] generic fit routine based on least-squares. will be used to fit calibration and resolution curves --- src/LegendSpecFits.jl | 1 + src/fit_chisq.jl | 54 +++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 55 insertions(+) create mode 100644 src/fit_chisq.jl diff --git a/src/LegendSpecFits.jl b/src/LegendSpecFits.jl index c20fdb81..cd6b0861 100644 --- a/src/LegendSpecFits.jl +++ b/src/LegendSpecFits.jl @@ -47,6 +47,7 @@ include("aoefit.jl") include("filter_optimization.jl") include("singlefit.jl") include("specfit.jl") +include("fit_chisq.jl") include("fwhm.jl") include("simple_calibration.jl") include("auto_calibration.jl") diff --git a/src/fit_chisq.jl b/src/fit_chisq.jl new file mode 100644 index 00000000..daa32501 --- /dev/null +++ b/src/fit_chisq.jl @@ -0,0 +1,54 @@ + +""" +fit_chisq(x::AbstractVector{<:Real},y::AbstractVector{<:Real},yerr::AbstractVector{<:Real}, f_fit::Function;pull_t::Vector{<:NamedTuple} = fill(NamedTuple(), first(methods(f_fit)).nargs - 2), v_init::Vector = []) +- least square fit with chi2 minimization +# input: +- x : x-values +- y : y-values +- yerr : 1 sigma uncertainty on y +- f_fit : fit/model function. e.g. for a linear function: f_lin(x,p1,p2) = p1 .* x .+ p2 +The numer of fit parameter is determined with `first(methods(f_fit)).nargs - 2`. That's why it's important that f_fit has the synthax f(x,arg1,arg2,arg3,...) +pull_t : pull term, a vector of NamedTuple with fields `mean` and `std`. A Gaussian pull term is added to the chi2 function to account for systematic uncertainties. If left blank, no pull term is used. +v_init : initial value for fit parameter optimization. If left blank, the initial value is set to 1 or guessed roughly for all fit parameters +""" +function fit_chisq(x::AbstractVector{<:Real},y::AbstractVector{<:Real},yerr::AbstractVector{<:Real}, f_fit::Function;pull_t::Vector{<:NamedTuple} = fill(NamedTuple(), first(methods(f_fit)).nargs - 2), v_init::Vector = []) + # prepare pull terms + f_pull(v::Number,pull_t::NamedTuple) = isempty(pull_t) ? zero(v) : (v .- pull_t.mean) .^2 ./ pull_t.std.^2 # pull term is zero if pull_t is zero + f_pull(v::Vector,pull_t::Vector) = sum(f_pull.(v,pull_t)) + pull_t_sum = Base.Fix2(f_pull, pull_t) + + # chi2 function + f_chi2 = let x = x, y = y, yerr = yerr, f_fit = f_fit, f_pull = pull_t_sum + v -> sum((y - f_fit(x, v...) ).^2 ./ yerr.^2) + f_pull(v) + end + + # init guess for fit parameter: this could be improved. + npar = first(methods(f_fit)).nargs - 2 # number of fit parameter (including nuisance parameters) + if isempty(v_init) + if npar==2 + v_init = [y[1]/x[1],1.0] # linear fit : guess slope + else + v_init = ones(npar) + end + end + + # minimization and error estimation + opt_r = optimize(f_chi2, v_init) + v_chi2 = Optim.minimizer(opt_r) + covmat = inv(ForwardDiff.hessian(f_chi2, v_chi2)) + v_chi2_err = sqrt.(diag(abs.(covmat))) + + # gof + chi2min = minimum(opt_r) + dof = length(x)-length(v_chi2) + pvalue = ccdf(Chisq(dof),chi2min) + + # result + result = (par = v_chi2, + err = v_chi2_err, + gof = (pvalue = pvalue, chi2min = chi2min, dof = dof, covmat = covmat)) + + return result +end + +export fit_chisq \ No newline at end of file From 6efee749e41456d54d8a306210b2da4b4c5be63a Mon Sep 17 00:00:00 2001 From: LisaSchlueter Date: Tue, 12 Mar 2024 22:58:56 +0100 Subject: [PATCH 02/68] test for fit_chisq fit routine. do linear and quadratic fit on dummydata --- test/runtests.jl | 1 + test/test_fit_chisq.jl | 27 +++++++++++++++++++++++++++ 2 files changed, 28 insertions(+) create mode 100644 test/test_fit_chisq.jl diff --git a/test/runtests.jl b/test/runtests.jl index 58172ab6..839ceaad 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -5,6 +5,7 @@ import Test Test.@testset "Package LegendSpecFits" begin #include("test_aqua.jl") include("test_specfit.jl") + include("test_fit_chisq.jl") include("test_docs.jl") isempty(Test.detect_ambiguities(LegendSpecFits)) end # testset diff --git a/test/test_fit_chisq.jl b/test/test_fit_chisq.jl new file mode 100644 index 00000000..e0b17f64 --- /dev/null +++ b/test/test_fit_chisq.jl @@ -0,0 +1,27 @@ + +using LegendSpecFits +using Test + +@testset "fit_chisq" begin + # linear fit : simple case with and without pull term + f_lin(x,p1,p2) = p1 .* x .+ p2 + + x = [1.0,2.0,3.0,4.0,5.0] + y = f_lin(x,5.0,10.0) + yerr = sqrt.(y) + @info "linear chisq fit" + result = fit_chisq(x,y,yerr,f_lin) + @info "linear chisq fit with pull term" + result = fit_chisq(x,y,yerr,f_lin;pull_t = [NamedTuple(), (mean = 10.0, std = 0.1)]) + + + # quadratic fit : simple case with and without pull term + f_quad(x,p1,p2,p3) = p1 .* x.^2 .+ p2 .* x .+ p3 + x = [1.0,2.0,3.0,4.0,5.0] + y = f_quad(x,2.0,5.0,10.0) + yerr = sqrt.(y) + @info "quadratic chisq fit" + result = fit_chisq(x,y,yerr,f_quad) + @info "quadratic chisq fit with pull term" + result = fit_chisq(x,y,yerr,f_quad;pull_t = [(mean = 2.0, std = 0.1),NamedTuple() ,NamedTuple(), ]) +end \ No newline at end of file From ec6459d95ebc6e298fb73b4004c8f2ac444ed96d Mon Sep 17 00:00:00 2001 From: LisaSchlueter Date: Fri, 15 Mar 2024 15:54:16 +0100 Subject: [PATCH 03/68] chi2fit for Measurements --- src/LegendSpecFits.jl | 7 +++- src/chi2fit.jl | 92 ++++++++++++++++++++++++++++++++++++++++++ src/fit_chisq.jl | 54 ------------------------- test/Project.toml | 2 + test/test_fit_chisq.jl | 40 +++++++++--------- 5 files changed, 118 insertions(+), 77 deletions(-) create mode 100644 src/chi2fit.jl delete mode 100644 src/fit_chisq.jl diff --git a/src/LegendSpecFits.jl b/src/LegendSpecFits.jl index cd6b0861..f91a657b 100644 --- a/src/LegendSpecFits.jl +++ b/src/LegendSpecFits.jl @@ -47,7 +47,7 @@ include("aoefit.jl") include("filter_optimization.jl") include("singlefit.jl") include("specfit.jl") -include("fit_chisq.jl") +include("chi2fit.jl") include("fwhm.jl") include("simple_calibration.jl") include("auto_calibration.jl") @@ -58,4 +58,9 @@ include("qc.jl") include("gof.jl") include("precompile.jl") +abstract type UncertTag end +ForwardDiff.:(≺)(::Type{<:ForwardDiff.Tag}, ::Type{UncertTag}) = true +ForwardDiff.:(≺)(::Type{UncertTag}, ::Type{<:ForwardDiff.Tag}) = false +ForwardDiff.:(≺)(::Type{UncertTag}, ::Type{UncertTag}) = false + end # module diff --git a/src/chi2fit.jl b/src/chi2fit.jl new file mode 100644 index 00000000..e52cb0cb --- /dev/null +++ b/src/chi2fit.jl @@ -0,0 +1,92 @@ + +""" +fit_chisq(x::AbstractVector{<:Real},y::AbstractVector{<:Real},yerr::AbstractVector{<:Real}, f_fit::Function;pull_t::Vector{<:NamedTuple} = fill(NamedTuple(), first(methods(f_fit)).nargs - 2), v_init::Vector = []) +- least square fit with chi2 minimization +# input: +- x : x-values +- y : y-values +- yerr : 1 sigma uncertainty on y +- f_fit : fit/model function. e.g. for a linear function: f_lin(x,p1,p2) = p1 .* x .+ p2 +The numer of fit parameter is determined with `first(methods(f_fit)).nargs - 2`. That's why it's important that f_fit has the synthax f(x,arg1,arg2,arg3,...) +pull_t : pull term, a vector of NamedTuple with fields `mean` and `std`. A Gaussian pull term is added to the chi2 function to account for systematic uncertainties. If left blank, no pull term is used. +v_init : initial value for fit parameter optimization. If left blank, the initial value is set to 1 or guessed roughly for all fit parameters +""" +function chi2fit(f_fit::Function, x::AbstractVector{<:Union{Real,Measurement{<:Real}}}, y::AbstractVector{<:Union{Real,Measurement{<:Real}}}; pull_t::Vector{<:NamedTuple}=fill(NamedTuple(), first(methods(f_fit)).nargs - 2), v_init::Vector = [], uncertainty::Bool=true ) + # prepare pull terms + f_pull(v::Number,pull_t::NamedTuple) = isempty(pull_t) ? zero(v) : (v .- pull_t.mean) .^2 ./ pull_t.std.^2 # pull term is zero if pull_t is zero + f_pull(v::Vector,pull_t::Vector) = sum(f_pull.(v,pull_t)) + pull_t_sum = Base.Fix2(f_pull, pull_t) + + # get rid of measurements + X_val = mvalue.(x) + Y_val = mvalue.(y) + X_err= muncert.(x) + Y_err = muncert.(y) + if all(X_err .== 0) && all(Y_err .== 0) + Y_err = ones(length(Y_val)) + end + # muncert_res(x) = x==0 ? 1.0 : muncert(x) # if uncertainties are zero use 1.0 instead (equal weights) + + # chi2-function + function chi2xy(f_fit, x_val, x_err, y_val, y_err, pars) + dual_x = ForwardDiff.Dual{UncertTag}(x_val, x_err) + dual_y = f_fit(dual_x, pars...) + y_pred_val = ForwardDiff.value(UncertTag, dual_y) + y_pred_err = only(ForwardDiff.partials(UncertTag, dual_y)) + chi2 = (y_val - y_pred_val)^2 / (y_err^2 + y_pred_err^2) + return chi2 + end + + # function that is minimized -> chi-squared with pull terms + f_opt = let X_val = X_val, Y_val = Y_val, X_err = X_err, Y_err = Y_err, f_pull = pull_t_sum, f_fit = f_fit + v -> sum(chi2xy.(Ref(f_fit), X_val, X_err, Y_val, Y_err, Ref(v))) + f_pull(v) + end + + # init guess for fit parameter: this could be improved. + npar = length(pull_t) # number of fit parameter (including nuisance parameters) + if isempty(v_init) + if npar==2 + v_init = [Y_val[1]/X_val[1], 1.0] # linear fit : guess slope + else + v_init = ones(npar) + end + end + + # minimization and error estimation + opt_r = optimize(f_opt, v_init) + v_chi2 = Optim.minimizer(opt_r) + f_fit_v(x) = f_fit(x, v_chi2...) # fit function with optimized parameters + par = measurement.(v_chi2,Ref(NaN)) # if ucnertainty is not calculated, return NaN + result = (par = par, f_fit_v = f_fit_v) + report = (par = result.par, f_fit = f_fit, x= x, y = y) + + if uncertainty + covmat = inv(ForwardDiff.hessian(f_opt, v_chi2)) + v_chi2_err = sqrt.(diag(abs.(covmat)))#mvalue.(sqrt.(diag(abs.(covmat)))) + par = measurement.(v_chi2,v_chi2_err) + f_fit_v(x) = f_fit(x, par...) # fit function with optimized parameters + + # gof + chi2min = minimum(opt_r) + dof = length(x) - length(v_chi2) + pvalue = ccdf(Chisq(dof),chi2min) + result = (par = par, f_fit_v = f_fit_v, gof = (pvalue = pvalue, chi2min = chi2min, dof = dof, covmat = covmat)) + report = merge(report, (gof = result.gof,)) + end + + return result, report +end + +function chi2fit(n_poly::Int, x, y; pull_t::Vector{<:NamedTuple}=fill(NamedTuple(), n_poly+1), kwargs...) + @assert length(pull_t) == n_poly+1 "Length of pull_t does not match the order of the polynomial" + chi2fit((x, a...) -> PolCalFunc.(a...).(x), x, y; pull_t = pull_t, kwargs...) +end + +chi2fit(f_fit, x,y::AbstractVector{<:Real},yerr::AbstractVector{<:Real};kwargs...) = chi2fit(f_fit,x,measurement.(y,yerr);kwargs...) +chi2fit(f_fit,x::AbstractVector{<:Real}, y::AbstractVector{<:Real}, yerr::AbstractVector{<:Real}, xerr::AbstractVector{<:Real} ;kwargs...) = chi2fit(f_fit, measurement.(x,xerr), measurement.(y,yerr); kwargs...) + +export chi2fit + + + +#f5(x,p1,p2,p3) = PolCalFun(p1,p2,p3) \ No newline at end of file diff --git a/src/fit_chisq.jl b/src/fit_chisq.jl deleted file mode 100644 index daa32501..00000000 --- a/src/fit_chisq.jl +++ /dev/null @@ -1,54 +0,0 @@ - -""" -fit_chisq(x::AbstractVector{<:Real},y::AbstractVector{<:Real},yerr::AbstractVector{<:Real}, f_fit::Function;pull_t::Vector{<:NamedTuple} = fill(NamedTuple(), first(methods(f_fit)).nargs - 2), v_init::Vector = []) -- least square fit with chi2 minimization -# input: -- x : x-values -- y : y-values -- yerr : 1 sigma uncertainty on y -- f_fit : fit/model function. e.g. for a linear function: f_lin(x,p1,p2) = p1 .* x .+ p2 -The numer of fit parameter is determined with `first(methods(f_fit)).nargs - 2`. That's why it's important that f_fit has the synthax f(x,arg1,arg2,arg3,...) -pull_t : pull term, a vector of NamedTuple with fields `mean` and `std`. A Gaussian pull term is added to the chi2 function to account for systematic uncertainties. If left blank, no pull term is used. -v_init : initial value for fit parameter optimization. If left blank, the initial value is set to 1 or guessed roughly for all fit parameters -""" -function fit_chisq(x::AbstractVector{<:Real},y::AbstractVector{<:Real},yerr::AbstractVector{<:Real}, f_fit::Function;pull_t::Vector{<:NamedTuple} = fill(NamedTuple(), first(methods(f_fit)).nargs - 2), v_init::Vector = []) - # prepare pull terms - f_pull(v::Number,pull_t::NamedTuple) = isempty(pull_t) ? zero(v) : (v .- pull_t.mean) .^2 ./ pull_t.std.^2 # pull term is zero if pull_t is zero - f_pull(v::Vector,pull_t::Vector) = sum(f_pull.(v,pull_t)) - pull_t_sum = Base.Fix2(f_pull, pull_t) - - # chi2 function - f_chi2 = let x = x, y = y, yerr = yerr, f_fit = f_fit, f_pull = pull_t_sum - v -> sum((y - f_fit(x, v...) ).^2 ./ yerr.^2) + f_pull(v) - end - - # init guess for fit parameter: this could be improved. - npar = first(methods(f_fit)).nargs - 2 # number of fit parameter (including nuisance parameters) - if isempty(v_init) - if npar==2 - v_init = [y[1]/x[1],1.0] # linear fit : guess slope - else - v_init = ones(npar) - end - end - - # minimization and error estimation - opt_r = optimize(f_chi2, v_init) - v_chi2 = Optim.minimizer(opt_r) - covmat = inv(ForwardDiff.hessian(f_chi2, v_chi2)) - v_chi2_err = sqrt.(diag(abs.(covmat))) - - # gof - chi2min = minimum(opt_r) - dof = length(x)-length(v_chi2) - pvalue = ccdf(Chisq(dof),chi2min) - - # result - result = (par = v_chi2, - err = v_chi2_err, - gof = (pvalue = pvalue, chi2min = chi2min, dof = dof, covmat = covmat)) - - return result -end - -export fit_chisq \ No newline at end of file diff --git a/test/Project.toml b/test/Project.toml index 973f7ea5..5ffde185 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -5,6 +5,7 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" LegendDataTypes = "99e09c13-5545-5ee2-bfa2-77f358fb75d8" LegendHDF5IO = "c9265ca6-b027-5446-b1a4-febfa8dd10b0" +Measurements = "eff96d63-e80a-5855-80a2-b1b0885c5ab7" Optim = "429524aa-4258-5aef-a3af-852621145aeb" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" @@ -15,4 +16,5 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" Documenter = "1" Interpolations = "0.15" LegendDataTypes = "0.1" +Measurements = "2" Unitful = "1" \ No newline at end of file diff --git a/test/test_fit_chisq.jl b/test/test_fit_chisq.jl index e0b17f64..18cff90f 100644 --- a/test/test_fit_chisq.jl +++ b/test/test_fit_chisq.jl @@ -1,27 +1,23 @@ -using LegendSpecFits +using LegendSpecFits +using Measurements using Test @testset "fit_chisq" begin - # linear fit : simple case with and without pull term - f_lin(x,p1,p2) = p1 .* x .+ p2 - - x = [1.0,2.0,3.0,4.0,5.0] - y = f_lin(x,5.0,10.0) - yerr = sqrt.(y) - @info "linear chisq fit" - result = fit_chisq(x,y,yerr,f_lin) - @info "linear chisq fit with pull term" - result = fit_chisq(x,y,yerr,f_lin;pull_t = [NamedTuple(), (mean = 10.0, std = 0.1)]) - + 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) + 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) - # quadratic fit : simple case with and without pull term - f_quad(x,p1,p2,p3) = p1 .* x.^2 .+ p2 .* x .+ p3 - x = [1.0,2.0,3.0,4.0,5.0] - y = f_quad(x,2.0,5.0,10.0) - yerr = sqrt.(y) - @info "quadratic chisq fit" - result = fit_chisq(x,y,yerr,f_quad) - @info "quadratic chisq fit with pull term" - result = fit_chisq(x,y,yerr,f_quad;pull_t = [(mean = 2.0, std = 0.1),NamedTuple() ,NamedTuple(), ]) -end \ No newline at end of file + x = [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) + 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) +end From 0e4ee9bbdfa5d35f4909426ab5183cbdd52a6525 Mon Sep 17 00:00:00 2001 From: Florian Henkes Date: Fri, 15 Mar 2024 18:18:02 +0100 Subject: [PATCH 04/68] Updated bin_width calcualtion with binning per peak --- src/simple_calibration.jl | 24 +++++++++++++----------- 1 file changed, 13 insertions(+), 11 deletions(-) diff --git a/src/simple_calibration.jl b/src/simple_calibration.jl index 2a3b782e..ca63a67d 100644 --- a/src/simple_calibration.jl +++ b/src/simple_calibration.jl @@ -34,7 +34,7 @@ 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...) -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, proxy_binning_peak::Unitful.Energy{<:Real}=2103.5u"keV", proxy_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") # create initial peak search histogram h_uncal = fit(Histogram, e_uncal, nbins=n_bins) # search all possible peak candidates @@ -50,16 +50,17 @@ function simple_calibration_th228(e_uncal::Vector{<:Real}, th228_lines::Vector{< e_simple = e_uncal .* c e_unit = u"keV" # get peakhists and peakstats - peakhists, peakstats, h_calsimple, bin_width = get_peakhists_th228(e_simple, th228_lines, window_sizes; proxy_binning_peak=proxy_binning_peak, proxy_binning_peak_window=proxy_binning_peak_window, e_unit=e_unit) + 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) result = ( h_calsimple = h_calsimple, h_uncal = h_uncal, c = c, unit = e_unit, - bin_width = bin_width, - fep_guess = fep_guess, - peakhists = peakhists, + bin_width = median(bin_widths), + fep_guess = fep_guess, + peakbinwidths = bin_widths, + peakhists = peakhists, peakstats = peakstats ) report = ( @@ -82,18 +83,19 @@ Create histograms around the calibration lines and return the histograms and the * `peakhists`: array of histograms around the calibration lines * `peakstats`: array of statistics for the calibration line fits """ -function get_peakhists_th228(e::Vector{<:Unitful.Energy{<:Real}}, th228_lines::Vector{<:Unitful.Energy{<:Real}}, window_sizes::Vector{<:Tuple{Unitful.Energy{<:Real}, Unitful.Energy{<:Real}}},; e_unit::Unitful.EnergyUnits=u"keV", proxy_binning_peak::Unitful.Energy{<:Real}=2103.5u"keV", proxy_binning_peak_window::Unitful.Energy{<:Real}=10.0u"keV") - bin_window_cut = proxy_binning_peak - proxy_binning_peak_window .< e .< proxy_binning_peak + proxy_binning_peak_window +function get_peakhists_th228(e::Vector{<:Unitful.Energy{<:Real}}, th228_lines::Vector{<:Unitful.Energy{<:Real}}, window_sizes::Vector{<:Tuple{Unitful.Energy{<:Real}, Unitful.Energy{<:Real}}},; e_unit::Unitful.EnergyUnits=u"keV", binning_peak_window::Unitful.Energy{<:Real}=10.0u"keV") # get optimal binning for simple calibration - bin_width = get_friedman_diaconis_bin_width(e[bin_window_cut]) + bin_widths = [get_friedman_diaconis_bin_width(filter(in(peak - binning_peak_window..peak + binning_peak_window), e)) for peak in th228_lines] + # create histogram for simple calibration e_min, e_max = 0u"keV", 3000u"keV" - h = fit(Histogram, ustrip.(e_unit, e), ustrip(e_unit, e_min):ustrip(e_unit, bin_width):ustrip(e_unit, e_max)) + h = fit(Histogram, ustrip.(e_unit, e), ustrip(e_unit, e_min):ustrip(e_unit, median(bin_widths)):ustrip(e_unit, e_max)) # get histograms around calibration lines and peakstats - peakhists = LegendSpecFits.subhist.(Ref(h), [ustrip.(e_unit, (peak-first(window), peak+last(window))) for (peak, window) in zip(th228_lines, window_sizes)]) + peakhists = [fit(Histogram, ustrip.(e_unit, e), ustrip(e_unit, peak - first(window)):ustrip(e_unit, bin_width):ustrip(e_unit, peak + last(window))) for (peak, window, bin_width) in zip(th228_lines, window_sizes, bin_widths)] + # peakhists = LegendSpecFits.subhist.(Ref(h), [ustrip.(e_unit, (peak-first(window), peak+last(window))) for (peak, window) in zip(th228_lines, window_sizes)]) peakstats = StructArray(estimate_single_peak_stats.(peakhists)) - peakhists, peakstats, h, bin_width + peakhists, peakstats, h, bin_widths end export get_peakhists_th228 From f58fbc99121297d80082f02472f2c4a900511789 Mon Sep 17 00:00:00 2001 From: Florian Henkes Date: Fri, 15 Mar 2024 18:18:27 +0100 Subject: [PATCH 05/68] New peakstats estimation for background and steps, updated priors --- src/specfit.jl | 38 ++++++++++++++++++++++++++------------ 1 file changed, 26 insertions(+), 12 deletions(-) diff --git a/src/specfit.jl b/src/specfit.jl index 72cfe6c0..472c7230 100644 --- a/src/specfit.jl +++ b/src/specfit.jl @@ -40,6 +40,7 @@ export estimate_single_peak_stats function estimate_single_peak_stats_th228(h::Histogram{T}) where T<:Real W = h.weights E = first(h.edges) + bin_width = step(E) peak_amplitude, peak_idx = findmax(W) fwhm_idx_left = findfirst(w -> w >= (first(W) + peak_amplitude) / 2, W) fwhm_idx_right = findlast(w -> w >= (last(W) + peak_amplitude) / 2, W) @@ -61,9 +62,20 @@ function estimate_single_peak_stats_th228(h::Histogram{T}) where T<:Real peak_sigma = 1.0 peak_fwhm = 2.0 end - #peak_area = peak_amplitude * peak_sigma * sqrt(2*π) - mean_background = (first(W) + last(W)) / 2 + # peak_area = peak_amplitude * peak_sigma * sqrt(2*π) + # calculate mean background and step + idx_bkg_left = something(findfirst(x -> x >= peak_pos - 15*peak_sigma, E[2:end]), 7) + idx_bkg_right = something(findfirst(x -> x >= peak_pos + 15*peak_sigma, E[2:end]), length(W) - 7) + mean_background_left, mean_background_right = mean(view(W, 1:idx_bkg_left)), mean(view(W, idx_bkg_right:length(W))) + + mean_background_step = (mean_background_left - mean_background_right) / bin_width + mean_background = (mean_background_left + mean_background_right) / 2 / bin_width + mean_background_std = std(view(W, 1:idx_bkg_left)) / bin_width + + # sanity checks mean_background = ifelse(mean_background == 0, 0.01, mean_background) + mean_background_step = ifelse(mean_background_step < 1e-2, 1e-2, mean_background_step) + peak_counts = inv(0.761) * (sum(view(W,fwhm_idx_left:fwhm_idx_right)) - mean_background * peak_fwhm) peak_counts = ifelse(peak_counts < 0.0, inv(0.761) * sum(view(W,fwhm_idx_left:fwhm_idx_right)), peak_counts) if !isnan(peak_fwqm) @@ -75,7 +87,9 @@ function estimate_single_peak_stats_th228(h::Histogram{T}) where T<:Real peak_fwhm = peak_fwhm, peak_sigma = peak_sigma, peak_counts = peak_counts, - mean_background = mean_background + mean_background = mean_background, + mean_background_step = mean_background_step, + mean_background_std = mean_background_std ) end @@ -125,8 +139,8 @@ function fit_peaks_th228(peakhists::Array, peakstats::StructArray, th228_lines:: end # save results if !isnothing(e_unit) - keys_to_keep = [k for k in keys(result_peak) if k != :gof] - result_peak = merge(NamedTuple{Tuple(keys_to_keep)}([result_peak[k] .* e_unit for k in keys_to_keep]...), (gof = result_peak.gof, n = result_peak.n, background = n = result_peak.background)) + keys_with_unit = [:μ, :σ, :fwhm] + result_peak = merge(result_peak, NamedTuple{Tuple(keys_with_unit)}([result_peak[k] .* e_unit for k in keys_with_unit]...)) end result[peak] = result_peak report[peak] = report_peak @@ -143,7 +157,7 @@ Also, FWHM is calculated from the fitted peakshape with MC error propagation. Th * `result`: NamedTuple of the fit results containing values and errors * `report`: NamedTuple of the fit report which can be plotted """ -function fit_single_peak_th228(h::Histogram, ps::NamedTuple{(:peak_pos, :peak_fwhm, :peak_sigma, :peak_counts, :mean_background), NTuple{5, T}}; +function fit_single_peak_th228(h::Histogram, ps::NamedTuple{(:peak_pos, :peak_fwhm, :peak_sigma, :peak_counts, :mean_background, :mean_background_step, :mean_background_std), NTuple{7, T}}; uncertainty::Bool=true, low_e_tail::Bool=true, fixed_position::Bool=false, pseudo_prior::NamedTupleDist=NamedTupleDist(empty = true), fit_func::Symbol=:f_fit) where T<:Real # create standard pseudo priors @@ -151,10 +165,10 @@ function fit_single_peak_th228(h::Histogram, ps::NamedTuple{(:peak_pos, :peak_fw μ = ifelse(fixed_position, ConstValueDist(ps.peak_pos), Uniform(ps.peak_pos-10, ps.peak_pos+10)), σ = weibull_from_mx(ps.peak_sigma, 2*ps.peak_sigma), n = weibull_from_mx(ps.peak_counts, 2*ps.peak_counts), - step_amplitude = weibull_from_mx(ps.mean_background, 2*ps.mean_background), - skew_fraction = ifelse(low_e_tail, Uniform(0.005, 0.25), ConstValueDist(0.0)), - skew_width = ifelse(low_e_tail, LogUniform(0.001, 0.1), ConstValueDist(1.0)), - background = weibull_from_mx(ps.mean_background, 2*ps.mean_background), + step_amplitude = weibull_from_mx(ps.mean_background_step, ps.mean_background_step + 5*ps.mean_background_std), + skew_fraction = ifelse(low_e_tail, truncated(weibull_from_mx(0.01, 0.05), 0.0, 0.25), ConstValueDist(0.0)), + skew_width = ifelse(low_e_tail, weibull_from_mx(0.001, 1e-3), ConstValueDist(1.0)), + background = weibull_from_mx(ps.mean_background, ps.mean_background + 5*ps.mean_background_std), ) # use standard priors in case of no overwrites given if !(:empty in keys(pseudo_prior)) @@ -171,7 +185,7 @@ function fit_single_peak_th228(h::Histogram, ps::NamedTuple{(:peak_pos, :peak_fw f_trafo = BAT.DistributionTransform(Normal, pseudo_prior) # start values for MLE - v_init = mean(pseudo_prior) + v_init = Vector(mean(f_trafo.target_dist)) # create loglikehood function: f_loglike(v) that can be evaluated for any set of v (fit parameter) f_loglike = let f_fit=th228_fit_functions[fit_func], h=h @@ -179,7 +193,7 @@ function fit_single_peak_th228(h::Histogram, ps::NamedTuple{(:peak_pos, :peak_fw end # MLE - opt_r = optimize((-) ∘ f_loglike ∘ inverse(f_trafo), f_trafo(v_init), Optim.Options(time_limit = 60, iterations = 500)) + opt_r = optimize((-) ∘ f_loglike ∘ inverse(f_trafo), v_init, Optim.Options(time_limit = 60, iterations = 500)) # best fit results v_ml = inverse(f_trafo)(Optim.minimizer(opt_r)) From 92eeb5c1d3f0d2358474b39a6c26de6d34413c55 Mon Sep 17 00:00:00 2001 From: Florian Henkes Date: Mon, 18 Mar 2024 11:06:22 +0100 Subject: [PATCH 06/68] Fixed precompile error due to function overwrite --- src/chi2fit.jl | 26 ++++++++++++++------------ 1 file changed, 14 insertions(+), 12 deletions(-) diff --git a/src/chi2fit.jl b/src/chi2fit.jl index e52cb0cb..ae83d99c 100644 --- a/src/chi2fit.jl +++ b/src/chi2fit.jl @@ -1,8 +1,8 @@ """ -fit_chisq(x::AbstractVector{<:Real},y::AbstractVector{<:Real},yerr::AbstractVector{<:Real}, f_fit::Function;pull_t::Vector{<:NamedTuple} = fill(NamedTuple(), first(methods(f_fit)).nargs - 2), v_init::Vector = []) -- least square fit with chi2 minimization -# input: + fit_chisq(x::AbstractVector{<:Real},y::AbstractVector{<:Real},yerr::AbstractVector{<:Real}, f_fit::Function;pull_t::Vector{<:NamedTuple} = fill(NamedTuple(), first(methods(f_fit)).nargs - 2), v_init::Vector = []) +Least square fit with chi2 minimization +# Input: - x : x-values - y : y-values - yerr : 1 sigma uncertainty on y @@ -10,6 +10,10 @@ fit_chisq(x::AbstractVector{<:Real},y::AbstractVector{<:Real},yerr::AbstractVect The numer of fit parameter is determined with `first(methods(f_fit)).nargs - 2`. That's why it's important that f_fit has the synthax f(x,arg1,arg2,arg3,...) pull_t : pull term, a vector of NamedTuple with fields `mean` and `std`. A Gaussian pull term is added to the chi2 function to account for systematic uncertainties. If left blank, no pull term is used. v_init : initial value for fit parameter optimization. If left blank, the initial value is set to 1 or guessed roughly for all fit parameters +# Return: +- result : NamedTuple with the optimized fit parameter and the fit function +- report: + """ function chi2fit(f_fit::Function, x::AbstractVector{<:Union{Real,Measurement{<:Real}}}, y::AbstractVector{<:Union{Real,Measurement{<:Real}}}; pull_t::Vector{<:NamedTuple}=fill(NamedTuple(), first(methods(f_fit)).nargs - 2), v_init::Vector = [], uncertainty::Bool=true ) # prepare pull terms @@ -37,40 +41,38 @@ function chi2fit(f_fit::Function, x::AbstractVector{<:Union{Real,Measurement{<:R return chi2 end - # function that is minimized -> chi-squared with pull terms + # function that is minimized -> chi-squared with pull terms f_opt = let X_val = X_val, Y_val = Y_val, X_err = X_err, Y_err = Y_err, f_pull = pull_t_sum, f_fit = f_fit v -> sum(chi2xy.(Ref(f_fit), X_val, X_err, Y_val, Y_err, Ref(v))) + f_pull(v) end - # init guess for fit parameter: this could be improved. + # init guess for fit parameter: this could be improved. npar = length(pull_t) # number of fit parameter (including nuisance parameters) if isempty(v_init) if npar==2 v_init = [Y_val[1]/X_val[1], 1.0] # linear fit : guess slope else - v_init = ones(npar) + v_init = ones(npar) end end # minimization and error estimation opt_r = optimize(f_opt, v_init) v_chi2 = Optim.minimizer(opt_r) - f_fit_v(x) = f_fit(x, v_chi2...) # fit function with optimized parameters par = measurement.(v_chi2,Ref(NaN)) # if ucnertainty is not calculated, return NaN - result = (par = par, f_fit_v = f_fit_v) - report = (par = result.par, f_fit = f_fit, x= x, y = y) + result = (par = par, f_fit_v = x -> f_fit(x, v_chi2...)) # fit function with optimized parameters + report = (par = result.par, f_fit = f_fit, x = x, y = y) if uncertainty covmat = inv(ForwardDiff.hessian(f_opt, v_chi2)) v_chi2_err = sqrt.(diag(abs.(covmat)))#mvalue.(sqrt.(diag(abs.(covmat)))) - par = measurement.(v_chi2,v_chi2_err) - f_fit_v(x) = f_fit(x, par...) # fit function with optimized parameters + par = measurement.(v_chi2, v_chi2_err) # gof chi2min = minimum(opt_r) dof = length(x) - length(v_chi2) pvalue = ccdf(Chisq(dof),chi2min) - result = (par = par, f_fit_v = f_fit_v, gof = (pvalue = pvalue, chi2min = chi2min, dof = dof, covmat = covmat)) + result = (par = par, f_fit_v = x -> f_fit(x, par...), gof = (pvalue = pvalue, chi2min = chi2min, dof = dof, covmat = covmat)) report = merge(report, (gof = result.gof,)) end From c5e151ba7ac04fc55a63b3010d4c02c2d17989d6 Mon Sep 17 00:00:00 2001 From: Florian Henkes Date: Mon, 18 Mar 2024 11:24:40 +0100 Subject: [PATCH 07/68] Fixed NaN for mean_background_std in estimate peakstats --- src/specfit.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/specfit.jl b/src/specfit.jl index 472c7230..565c2a1a 100644 --- a/src/specfit.jl +++ b/src/specfit.jl @@ -75,6 +75,7 @@ function estimate_single_peak_stats_th228(h::Histogram{T}) where T<:Real # sanity checks mean_background = ifelse(mean_background == 0, 0.01, mean_background) mean_background_step = ifelse(mean_background_step < 1e-2, 1e-2, mean_background_step) + mean_background_std = ifelse(!isfinite(mean_background_std) || mean_background_std == 0, sqrt(mean_background), mean_background_std) peak_counts = inv(0.761) * (sum(view(W,fwhm_idx_left:fwhm_idx_right)) - mean_background * peak_fwhm) peak_counts = ifelse(peak_counts < 0.0, inv(0.761) * sum(view(W,fwhm_idx_left:fwhm_idx_right)), peak_counts) From 7b66057c1e95abef0cdfd0a1942c6cfff9aa9180 Mon Sep 17 00:00:00 2001 From: LisaSchlueter Date: Mon, 18 Mar 2024 17:30:17 +0100 Subject: [PATCH 08/68] fit calibration curve with chi2fit --- src/LegendSpecFits.jl | 1 + src/chi2fit.jl | 8 ++++---- src/fit_calibration.jl | 30 ++++++++++++++++++++++++++++++ 3 files changed, 35 insertions(+), 4 deletions(-) create mode 100644 src/fit_calibration.jl diff --git a/src/LegendSpecFits.jl b/src/LegendSpecFits.jl index f91a657b..07a95046 100644 --- a/src/LegendSpecFits.jl +++ b/src/LegendSpecFits.jl @@ -48,6 +48,7 @@ include("filter_optimization.jl") include("singlefit.jl") include("specfit.jl") include("chi2fit.jl") +include("fit_calibration.jl") include("fwhm.jl") include("simple_calibration.jl") include("auto_calibration.jl") diff --git a/src/chi2fit.jl b/src/chi2fit.jl index ae83d99c..e963a026 100644 --- a/src/chi2fit.jl +++ b/src/chi2fit.jl @@ -60,8 +60,8 @@ function chi2fit(f_fit::Function, x::AbstractVector{<:Union{Real,Measurement{<:R opt_r = optimize(f_opt, v_init) v_chi2 = Optim.minimizer(opt_r) par = measurement.(v_chi2,Ref(NaN)) # if ucnertainty is not calculated, return NaN - result = (par = par, f_fit_v = x -> f_fit(x, v_chi2...)) # fit function with optimized parameters - report = (par = result.par, f_fit = f_fit, x = x, y = y) + result = (par = par,) # fit function with optimized parameters + report = (par = result.par, f_fit = x -> f_fit(x, v_chi2...), x = x, y = y) if uncertainty covmat = inv(ForwardDiff.hessian(f_opt, v_chi2)) @@ -72,8 +72,8 @@ function chi2fit(f_fit::Function, x::AbstractVector{<:Union{Real,Measurement{<:R chi2min = minimum(opt_r) dof = length(x) - length(v_chi2) pvalue = ccdf(Chisq(dof),chi2min) - result = (par = par, f_fit_v = x -> f_fit(x, par...), gof = (pvalue = pvalue, chi2min = chi2min, dof = dof, covmat = covmat)) - report = merge(report, (gof = result.gof,)) + result = (par = par, gof = (pvalue = pvalue, chi2min = chi2min, dof = dof, covmat = covmat)) + report = merge(report, (par = result.par,), (gof = result.gof,), (f_fit = x -> f_fit(x, v_chi2...),)) end return result, report diff --git a/src/fit_calibration.jl b/src/fit_calibration.jl new file mode 100644 index 00000000..ed55627c --- /dev/null +++ b/src/fit_calibration.jl @@ -0,0 +1,30 @@ +""" + fit_calibration(n_poly::Int, µ::AbstractVector{<:Union{Real,Measurement{<:Real}}}, peaks::AbstractVector{<:Union{Real,Measurement{<:Real}}}; pull_t::Vector{<:NamedTuple}=fill(NamedTuple(), n_poly+1), v_init::Vector = [], uncertainty::Bool=true ) +Fit the calibration lines with polynomial function of n_poly order + n_poly == 1 -> linear function + n_poly == 2 -> quadratic function +# Returns + * `result`: NamedTuple with the following fields + * `par`: best-fit parameters + * `gof`: godness of fit + * `report`: +""" +function fit_calibration(n_poly::Int, µ::AbstractVector{<:Union{Real,Measurement{<:Real}}}, peaks::AbstractVector{<:Quantity};e_type::Union{Symbol,String}="x") + @assert length(peaks) == length(μ) + e_unit = u"keV" + @debug "Fit calibration curve with $(n_poly)-order polynominal function" + result_fit, report = chi2fit(n_poly, µ, ustrip.(e_unit,peaks)) + + par = result_fit.par + par_unit = par .* e_unit + + result = merge(result_fit, (par = par_unit,)) + + # built function in string + func = join(["$(mvalue(par[i])) * $(e_type)^$(i-1)" for i=1:length(par)], " + ") + func_generic = join(["p$(i-1) * $(e_type)^$(i-1)" for i=1:length(par)], " + ") + + result = merge(result, (func = func, func_generic = func_generic)) + return result, report +end +export fit_ecal \ No newline at end of file From de611eed3dc88649d16b0cd2ce983873f385a4e1 Mon Sep 17 00:00:00 2001 From: Florian Henkes Date: Tue, 19 Mar 2024 16:20:32 +0100 Subject: [PATCH 09/68] Gof output includes normalized residuals --- src/chi2fit.jl | 22 ++++++++++++++-------- 1 file changed, 14 insertions(+), 8 deletions(-) diff --git a/src/chi2fit.jl b/src/chi2fit.jl index ae83d99c..da166a00 100644 --- a/src/chi2fit.jl +++ b/src/chi2fit.jl @@ -60,8 +60,8 @@ function chi2fit(f_fit::Function, x::AbstractVector{<:Union{Real,Measurement{<:R opt_r = optimize(f_opt, v_init) v_chi2 = Optim.minimizer(opt_r) par = measurement.(v_chi2,Ref(NaN)) # if ucnertainty is not calculated, return NaN - result = (par = par, f_fit_v = x -> f_fit(x, v_chi2...)) # fit function with optimized parameters - report = (par = result.par, f_fit = f_fit, x = x, y = y) + result = (par = par, ) # fit function with optimized parameters + report = (par = result.par, f_fit = x -> f_fit(x, v_chi2...), x = x, y = y, gof = NamedTuple()) if uncertainty covmat = inv(ForwardDiff.hessian(f_opt, v_chi2)) @@ -71,9 +71,10 @@ function chi2fit(f_fit::Function, x::AbstractVector{<:Union{Real,Measurement{<:R # gof chi2min = minimum(opt_r) dof = length(x) - length(v_chi2) - pvalue = ccdf(Chisq(dof),chi2min) - result = (par = par, f_fit_v = x -> f_fit(x, par...), gof = (pvalue = pvalue, chi2min = chi2min, dof = dof, covmat = covmat)) - report = merge(report, (gof = result.gof,)) + pvalue = ccdf(Chisq(dof), chi2min) + residuals_norm = (Y_val - f_fit(X_val, v_chi2...)) ./ Y_err + result = (par = par, gof = (pvalue = pvalue, chi2min = chi2min, dof = dof, covmat = covmat, residuals_norm = residuals_norm)) + report = merge(report, (par = par, gof = result.gof, f_fit = x -> f_fit(x, par...))) end return result, report @@ -81,11 +82,16 @@ end function chi2fit(n_poly::Int, x, y; pull_t::Vector{<:NamedTuple}=fill(NamedTuple(), n_poly+1), kwargs...) @assert length(pull_t) == n_poly+1 "Length of pull_t does not match the order of the polynomial" - chi2fit((x, a...) -> PolCalFunc.(a...).(x), x, y; pull_t = pull_t, kwargs...) + chi2fit((x, a...) -> PolCalFunc(a...).(x), x, y; pull_t = pull_t, kwargs...) end -chi2fit(f_fit, x,y::AbstractVector{<:Real},yerr::AbstractVector{<:Real};kwargs...) = chi2fit(f_fit,x,measurement.(y,yerr);kwargs...) -chi2fit(f_fit,x::AbstractVector{<:Real}, y::AbstractVector{<:Real}, yerr::AbstractVector{<:Real}, xerr::AbstractVector{<:Real} ;kwargs...) = chi2fit(f_fit, measurement.(x,xerr), measurement.(y,yerr); kwargs...) +function chi2fit(f_outer::Function, n_poly::Int, x, y; pull_t::Vector{<:NamedTuple}=fill(NamedTuple(), n_poly + 1), kwargs...) + @assert length(pull_t) == n_poly + 1 "Length of pull_t does not match the order of the polynomial" + chi2fit((x, a...) -> f_outer.(PolCalFunc(a...).(x)), x, y; pull_t = pull_t, kwargs...) +end + +chi2fit(f_fit, x::AbstractVector, y::AbstractVector{<:Real}, yerr::AbstractVector{<:Real}; kwargs...) = chi2fit(f_fit, x, measurement.(y,yerr); kwargs...) +chi2fit(f_fit, x::AbstractVector{<:Real}, y::AbstractVector{<:Real}, yerr::AbstractVector{<:Real}, xerr::AbstractVector{<:Real}; kwargs...) = chi2fit(f_fit, measurement.(x, xerr), measurement.(y, yerr); kwargs...) export chi2fit From d084cfca77fc25772e13e0649ae786c9f43425ac Mon Sep 17 00:00:00 2001 From: Florian Henkes Date: Tue, 19 Mar 2024 16:21:59 +0100 Subject: [PATCH 10/68] Recipes for fwhm and peak fits with residuals --- ext/LegendSpecFitsRecipesBaseExt.jl | 209 ++++++++++++++++++++++++++-- src/fwhm.jl | 47 +++---- 2 files changed, 219 insertions(+), 37 deletions(-) diff --git a/ext/LegendSpecFitsRecipesBaseExt.jl b/ext/LegendSpecFitsRecipesBaseExt.jl index a8a5a460..d4db3132 100644 --- a/ext/LegendSpecFitsRecipesBaseExt.jl +++ b/ext/LegendSpecFitsRecipesBaseExt.jl @@ -3,7 +3,7 @@ module LegendSpecFitsRecipesBaseExt using RecipesBase -using Unitful, Formatting, Measurements +using Unitful, Formatting, Measurements, LaTeXStrings using Measurements: value, uncertainty using StatsBase, LinearAlgebra @@ -104,48 +104,146 @@ end end end -@recipe function f(report::NamedTuple{(:v, :h, :f_fit, :f_sig, :f_lowEtail, :f_bck)}; show_label=true, show_fit=true) - xlabel := "Energy (keV)" - ylabel := "Counts" - legend := :bottomright - yscale := :log10 - ylim_max = max(5*report.f_sig(report.v.μ), 5*maximum(report.h.weights)) +@recipe function f(report::NamedTuple{(:v, :h, :f_fit, :f_sig, :f_lowEtail, :f_bck, :gof)}; show_label=true, show_fit=true, _subplot=1) + # thickness_scaling := 2.0 + legend := :topright + ylim_max = max(3*value(report.f_sig(report.v.μ)), 3*maximum(report.h.weights)) ylim_max = ifelse(ylim_max == 0.0, 1e5, ylim_max) - ylims := (1, ylim_max) + ylim_min = 0.1*minimum(filter(x -> x > 0, report.h.weights)) + framestyle --> :box @series begin seriestype := :stepbins label := ifelse(show_label, "Data", "") bins --> :sqrt + yscale := :log10 + ylims := (ylim_min, ylim_max) + xlabel := "Energy (keV)" + xlims := (minimum(report.h.edges[1]), maximum(report.h.edges[1])) + ylabel := "Counts / $(round(step(report.h.edges[1]), digits=2)) keV" + subplot --> _subplot LinearAlgebra.normalize(report.h, mode = :density) end if show_fit @series begin seriestype := :line - label := ifelse(show_label, "Best Fit", "") + if !isempty(report.gof) + label := ifelse(show_label, "Best Fit (p = $(round(report.gof.pvalue, digits=2)))", "") + else + label := ifelse(show_label, "Best Fit", "") + end color := :red - minimum(report.h.edges[1]):0.1:maximum(report.h.edges[1]), report.f_fit + subplot --> _subplot + # ribbon := uncertainty.(report.f_fit.(minimum(report.h.edges[1]):0.1:maximum(report.h.edges[1]))) + minimum(report.h.edges[1]):0.1:maximum(report.h.edges[1]), value.(report.f_fit.(minimum(report.h.edges[1]):0.1:maximum(report.h.edges[1]))) end @series begin seriestype := :line label := ifelse(show_label, "Signal", "") + subplot --> _subplot color := :green minimum(report.h.edges[1]):0.1:maximum(report.h.edges[1]), report.f_sig end @series begin seriestype := :line label := ifelse(show_label, "Low Energy Tail", "") + subplot --> _subplot color := :blue minimum(report.h.edges[1]):0.1:maximum(report.h.edges[1]), report.f_lowEtail end @series begin seriestype := :line label := ifelse(show_label, "Background", "") + subplot --> _subplot + margins --> (0, :mm) + bottom_margin --> (-4, :mm) + xlabel := "" + xticks --> ([]) + ylabel := "Counts / $(round(step(report.h.edges[1]), digits=2)) keV" + ylims := (ylim_min, ylim_max) + xlims := (minimum(report.h.edges[1]), maximum(report.h.edges[1])) color := :black minimum(report.h.edges[1]):0.1:maximum(report.h.edges[1]), report.f_bck end + if !isempty(report.gof) + ylims_res_max, ylims_res_min = 5, -5 + if any(report.gof.residuals_norm .> 5) || any(report.gof.residuals_norm .< -5) + abs_max = 1.2*maximum(abs.(report.gof.residuals_norm)) + ylims_res_max, ylims_res_min = abs_max, -abs_max + end + layout --> @layout([a{0.8h}; b{0.2h}]) + margins --> (0, :mm) + link --> :x + @series begin + seriestype := :hline + ribbon := 3 + subplot --> _subplot + 1 + fillalpha := 0.5 + label := "" + fillcolor := :lightgrey + linecolor := :darkgrey + [0.0] + end + @series begin + seriestype := :hline + ribbon := 1 + subplot --> _subplot + 1 + fillalpha := 0.5 + label := "" + fillcolor := :grey + linecolor := :darkgrey + [0.0] + end + @series begin + seriestype := :scatter + subplot --> _subplot + 1 + label := "" + title := "" + markercolor --> :black + ylabel --> "Residuals (σ)" + xlabel := "Energy (keV)" + link --> :x + top_margin --> (0, :mm) + ylims --> (ylims_res_min, ylims_res_max) + xlims := (minimum(report.h.edges[1]), maximum(report.h.edges[1])) + yscale --> :identity + if ylims_res_max == 5 + yticks --> ([-3, 0, 3]) + end + report.gof.bin_centers, report.gof.residuals_norm + end + else + @series begin + seriestype := :line + label := ifelse(show_label, "Background", "") + subplot --> _subplot + margins --> (0, :mm) + bottom_margin --> (-4, :mm) + xlabel := "Energy (keV)" + xticks --> ([]) + ylabel := "Counts / $(round(step(report.h.edges[1]), digits=2)) keV" + ylims := (ylim_min, ylim_max) + xlims := (minimum(report.h.edges[1]), maximum(report.h.edges[1])) + color := :black + minimum(report.h.edges[1]):0.1:maximum(report.h.edges[1]), report.f_bck + end + end end end +# TODO: Add a recipe for the report_dict --> Feeeeeeeeeeeeelix :* +# @recipe function f(report_dict::Dict{Symbol, NamedTuple}) +# # layout := grid(length(report_dict)*2, 1, heights=vcat(fill([0.8, 0.2], length(report_dict))...)) +# size := (1800, 1000*length(report_dict)) +# bottom_margin := (0, :mm) +# for (i, k) in enumerate(string.(keys(report_dict))) +# @series begin +# title := string(k) +# _subplot --> 2*i-1 +# report_dict[Symbol(k)] +# end +# end +# end + @recipe function f(report::NamedTuple{((:v, :h, :f_fit, :f_sig, :f_bck))}) xlabel := "A/E (a.u.)" ylabel := "Counts" @@ -370,6 +468,97 @@ end end +@recipe function f(report::NamedTuple{(:par, :f_fit, :x, :y, :gof)}; plot_ribbon=true) + thickness_scaling := 2.0 + xlims := (0, 1.2*maximum(report.x)) + framestyle := :box + xformatter := :plain + yformatter := :plain + layout --> @layout([a{0.8h}; b{0.2h}]) + margins --> (-15, :mm) + link --> :x + @series begin + seriestype := :line + subplot := 1 + xticks = :none + if !isempty(report.gof) + label := "Best Fit (p = $(round(report.gof.pvalue, digits=2)))" + else + label := "Best Fit" + end + color := :orange + linewidth := 2 + fillalpha := 0.2 + if plot_ribbon + ribbon := uncertainty.(report.f_fit.(0:1:1.2*maximum(report.x))) + end + 0:1:1.2*maximum(report.x), value.(report.f_fit.(0:1:1.2*maximum(report.x))) + end + @series begin + seriestype := :scatter + subplot := 1 + label := "Data" + markercolor --> :black + report.x, report.y + end + @series begin + seriestype := :hline + ribbon := 3 + subplot := 2 + fillalpha := 0.5 + label := "" + fillcolor := :lightgrey + linecolor := :darkgrey + [0.0] + end + @series begin + seriestype := :hline + ribbon := 1 + subplot := 2 + fillalpha := 0.5 + label := "" + fillcolor := :grey + linecolor := :darkgrey + [0.0] + end + @series begin + seriestype := :scatter + subplot := 2 + label := "" + markercolor --> :black + ylabel --> "\n\nResiduals (σ)" + ylims --> (-5, 5) + yticks --> ([-3, 0, 3]) + report.x, report.gof.residuals_norm + end +end + +@recipe function f(report::NamedTuple{(:par, :f_fit, :x, :y, :gof, :qbb, :type)}) + bottom_margin --> (0, :mm) + if report.type == :fwhm + xlabel := "Energy (keV)" + legend := :bottomright + framestyle := :box + xlims := (0, 3000) + xticks := (0:500:3000, ["$i" for i in 0:500:3000]) + @series begin + grid --> :all + (par = report.par, f_fit = report.f_fit, x = report.x, y = report.y, gof = report.gof) + end + @series begin + seriestype := :hline + label := L"Q_{\beta \beta}:" * " $(round(u"keV", report.qbb, digits=2))" + color := :green + fillalpha := 0.2 + linewidth := 2.5 + xticks := :none + ylabel := "FWHM" + subplot := 1 + ribbon := uncertainty(report.qbb) + [value(report.qbb)] + end + end +end # @recipe function f(x::AbstractArray, y::AbstractArray{<:Measurement}; plot_ribbon = false) # if plot_ribbon diff --git a/src/fwhm.jl b/src/fwhm.jl index 20122e49..a24706c8 100644 --- a/src/fwhm.jl +++ b/src/fwhm.jl @@ -1,6 +1,7 @@ # f_fwhm(x, p) = sqrt.((x .* x .* p[3] .+ x .* p[2] .+ p[1]) .* heaviside.(x .* x .* p[3] .+ x .* p[2] .+ p[1])) f_fwhm(x::T, p::AbstractArray{<:T}) where T<:Unitful.RealOrRealQuantity = sqrt((x * x * p[3] + x * p[2] + p[1]) * heaviside(x^2 * p[3] + x * p[2] + p[1])) f_fwhm(x::Array{<:T}, p::AbstractArray{<:T}) where T<:Unitful.RealOrRealQuantity = Base.Fix2(f_fwhm, p).(x) +f_fwhm(x, p1, p2, p3) = f_fwhm(x, [p1, p2, p3]) """ fitFWHM(fit_fwhm(peaks::Vector{T}, fwhm::Vector{T}) where T<:Real @@ -11,37 +12,29 @@ Fit the FWHM of the peaks to a quadratic function. * `v`: the fit result parameters * `f_fit`: the fitted function """ -function fit_fwhm(peaks::Vector{<:Unitful.Energy{<:Real}}, fwhm::Vector{<:Unitful.Energy{<:Real}}) +function fit_fwhm(peaks::Vector{<:Unitful.Energy{<:Real}}, fwhm::Vector{<:Unitful.Energy{<:Real}}; pol_order::Int=1, e_type_cal::Symbol=:e_cal, e_expression::Union{Symbol, String}="e", uncertainty::Bool=true) + @assert length(peaks) == length(fwhm) "Peaks and FWHM must have the same length" # fit FWHM fit function e_unit = u"keV" - p_start = [1*e_unit^2, 0.001*e_unit, 0.0*e_unit] - lower_bound = [0.0*e_unit^2, 0.0*e_unit, 0.0*e_unit] - fwhm_fit_result = curve_fit(f_fwhm, ustrip.(e_unit, peaks), ustrip.(e_unit, fwhm), ustrip.(p_start), lower=ustrip.(lower_bound)) + p_start = append!([1, 2.96e-3*0.11], fill(0.0, pol_order-1)) .* [e_unit^i for i in pol_order:-1:0] + + # fit FWHM fit function as a square root of a polynomial + result_chi2, report_chi2 = chi2fit(x -> LegendSpecFits.heaviside(x)*sqrt(abs(x)), pol_order, ustrip.(e_unit, peaks), ustrip.(e_unit, fwhm); v_init=ustrip.(p_start), uncertainty=uncertainty) + + # get pars and apply unit + par = result_chi2.par + par_unit = par .* [e_unit^i for i in pol_order:-1:0] - # get FWHM at Qbb with error - err = standard_errors(fwhm_fit_result) - enc, fano, ct = fwhm_fit_result.param[1]*e_unit^2, fwhm_fit_result.param[2]*e_unit, fwhm_fit_result.param[3] - fwhm_qbb = f_fwhm(2039u"keV", [enc, fano, ct]) - fwhm_pars_rand = rand(MvNormal(fwhm_fit_result.param, estimate_covar(fwhm_fit_result)), 10000) - # prevent negative parameter values in fit - # TODO: check if this is the right way to do this - fwhm_pars_rand_cleaned = fwhm_pars_rand[:, findall(x -> all(x .> 0), eachcol(fwhm_pars_rand))] - fwhm_pars_rand_cleaned = [[p1, p2, p3] for (p1, p2, p3) in zip(fwhm_pars_rand_cleaned[1,:] .*e_unit^2, fwhm_pars_rand_cleaned[2,:] .* e_unit, fwhm_pars_rand_cleaned[3,:])] - # get all FWHM at Qbb with error - fwhm_qbb_rand = f_fwhm.(2039u"keV", fwhm_pars_rand_cleaned) - fwhm_qbb_err = std(fwhm_qbb_rand) + # built function in string + func = "sqrt($(join(["$(mvalue(par[i])) * ($(e_expression))^$(i-1)" for i=1:length(par)], " + ")))$e_unit" + func_cal = "sqrt($(join(["$(mvalue(par[i])) * $(e_type_cal)^$(i-1) * keV^$(length(par)+1-i)" for i=1:length(par)], " + ")))" + func_generic = "sqrt($(join(["par[$(i-1)] * $(e_type_cal)^$(i-1)" for i=1:length(par)], " + ")))" + + # get fwhm at Qbb + qbb = report_chi2.f_fit(2039) * e_unit + result = merge(result_chi2, (par = par_unit , qbb = qbb, func = func, func_cal = func_cal, func_generic = func_generic)) + report = merge(report_chi2, (par = result.par, qbb = result.qbb, type = :fwhm)) - result = ( - qbb = measurement(fwhm_qbb, fwhm_qbb_err), - enc = measurement(fwhm_fit_result.param[1]*e_unit^2, err[1]*e_unit^2), - fano = measurement(fwhm_fit_result.param[2]*e_unit, err[2]*e_unit), - ct = measurement(fwhm_fit_result.param[3], err[3]) - ) - report = ( - qbb = result.qbb, - v = [enc, fano, ct], - f_fit = x -> Base.Fix2(f_fwhm, [enc, fano, ct])(x) - ) return result, report end export fit_fwhm \ No newline at end of file From e31cc46474ff923a18a51b71bb387e66107f3af7 Mon Sep 17 00:00:00 2001 From: Florian Henkes Date: Tue, 19 Mar 2024 16:22:28 +0100 Subject: [PATCH 11/68] report f_fit with uncertainties --- src/specfit.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/specfit.jl b/src/specfit.jl index 565c2a1a..4d646925 100644 --- a/src/specfit.jl +++ b/src/specfit.jl @@ -242,10 +242,11 @@ function fit_single_peak_th228(h::Histogram, ps::NamedTuple{(:peak_pos, :peak_fw report = ( v = v_ml, h = h, - f_fit = x -> Base.Fix2(th228_fit_functions.f_fit, v_ml)(x), + f_fit = x -> Base.Fix2(th228_fit_functions.f_fit, result)(x), f_sig = x -> Base.Fix2(th228_fit_functions.f_sig, v_ml)(x), f_lowEtail = x -> Base.Fix2(th228_fit_functions.f_lowEtail, v_ml)(x), f_bck = x -> Base.Fix2(th228_fit_functions.f_bck, v_ml)(x), + gof = result.gof ) else # get fwhm of peak @@ -266,7 +267,8 @@ function fit_single_peak_th228(h::Histogram, ps::NamedTuple{(:peak_pos, :peak_fw f_fit = x -> Base.Fix2(th228_fit_functions.f_fit, v_ml)(x), f_sig = x -> Base.Fix2(th228_fit_functions.f_sig, v_ml)(x), f_lowEtail = x -> Base.Fix2(th228_fit_functions.f_lowEtail, v_ml)(x), - f_bck = x -> Base.Fix2(th228_fit_functions.f_bck, v_ml)(x) + f_bck = x -> Base.Fix2(th228_fit_functions.f_bck, v_ml)(x), + gof = NamedTuple() ) end return result, report From 30d870755719c09a79552feb4e11b11330a28bb8 Mon Sep 17 00:00:00 2001 From: Florian Henkes Date: Tue, 19 Mar 2024 16:22:50 +0100 Subject: [PATCH 12/68] Require LatexStrings and Measures for plotting --- Project.toml | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/Project.toml b/Project.toml index b6a65bc2..83250c76 100644 --- a/Project.toml +++ b/Project.toml @@ -15,10 +15,12 @@ ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953" InverseFunctions = "3587e190-3f89-42d0-90ee-14403ec27112" IrrationalConstants = "92d709cd-6900-40b7-9082-c6be49f344b6" +LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f" LegendDataManagement = "9feedd95-f0e0-423f-a8dc-de0970eae6b3" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" LinearRegression = "92481ed7-9fb7-40fd-80f2-46fd0f076581" LsqFit = "2fda8390-95c7-5789-9bda-21331edee243" +Measures = "442fdcdd-2543-5da2-b0f3-8c86c306513e" Measurements = "eff96d63-e80a-5855-80a2-b1b0885c5ab7" Optim = "429524aa-4258-5aef-a3af-852621145aeb" PropDicts = "4dc08600-4268-439e-8673-d706fafbb426" @@ -55,10 +57,12 @@ ForwardDiff = "0.10" IntervalSets = "0.7" InverseFunctions = "0.1" IrrationalConstants = "0.1, 0.2" +LaTeXStrings = "1.3" LegendDataManagement = "0.2.7" LinearAlgebra = "1" LinearRegression = "0.2" LsqFit = "0.14, 0.15" +Measures = "0.3" Measurements = "2" Optim = "1" PropDicts = "0.2" From 2a3e3d77fe9b87826e7b8e4e1ccb711d95443f55 Mon Sep 17 00:00:00 2001 From: LisaSchlueter Date: Wed, 20 Mar 2024 11:26:24 +0100 Subject: [PATCH 13/68] add pvalue in ml-fit of compton bands and add chi2fit --- ext/LegendSpecFitsRecipesBaseExt.jl | 18 ++++++++++++++++++ src/aoe_calibration.jl | 19 +++++++++++-------- src/aoefit.jl | 29 ++++++++++------------------- src/gof.jl | 4 ++-- 4 files changed, 41 insertions(+), 29 deletions(-) diff --git a/ext/LegendSpecFitsRecipesBaseExt.jl b/ext/LegendSpecFitsRecipesBaseExt.jl index a8a5a460..b9857ab0 100644 --- a/ext/LegendSpecFitsRecipesBaseExt.jl +++ b/ext/LegendSpecFitsRecipesBaseExt.jl @@ -378,5 +378,23 @@ end # x, value.(y) # end +# @recipe function f(covmat::Matrix{<:Real};mode::Symbol=:cov) +# if mode==:cov +# colormap = :viridis +# title_str = "covariance" +# elseif mode==:corr +# cm = cor(cm); +# colormap = :greys +# title_str = "correlation" +# end +# @series begin +# seriestype := :heatmap +# yflip := true +# colorbar_title_location := :right +# guidefontsize := 16 +# tickfontsize := 14 +# size := (1000, 600) +# end +# end end # module LegendSpecFitsRecipesBaseExt \ No newline at end of file diff --git a/src/aoe_calibration.jl b/src/aoe_calibration.jl index b7a95466..2fd6290c 100644 --- a/src/aoe_calibration.jl +++ b/src/aoe_calibration.jl @@ -20,16 +20,19 @@ function fit_aoe_corrections(e::Array{<:Unitful.Energy{<:Real}}, μ::Array{<:Rea # fit compton band mus with linear function μ_cut = (mean(μ) - 2*std(μ) .< μ .< mean(μ) + 2*std(μ)) e, μ, σ = e[μ_cut], μ[μ_cut], σ[μ_cut] - # μ_scs = linregress(e[μ_cut], μ[μ_cut]) e_unit = unit(first(e)) e = ustrip.(e_unit, e) - # fit compton band mus with linear function - μ_scs = linregress(e, μ) - μ_scs_slope, μ_scs_intercept = LinearRegression.slope(μ_scs)[1], LinearRegression.bias(μ_scs)[1] - μ_scs = curve_fit(f_aoe_mu, e, μ, [μ_scs_intercept, μ_scs_slope]) - μ_scs_param = [μ_scs.param[1], μ_scs.param[2] / e_unit] - @debug "μ_scs_slope : $μ_scs_slope" - @debug "μ_scs_intercept: $μ_scs_intercept" + + # fit compton band with linear function + result_µ, report_µ = LegendSpecFits.chi2fit(1, e, µ; uncertainty=true) +result_µ + + # μ_scs = linregress(e, μ) + # μ_scs_slope, μ_scs_intercept = LinearRegression.slope(μ_scs)[1], LinearRegression.bias(μ_scs)[1] + # μ_scs = curve_fit(f_aoe_mu, e, μ, [μ_scs_intercept, μ_scs_slope]) + # μ_scs_param = [μ_scs.param[1], μ_scs.param[2] / e_unit] + # @debug "μ_scs_slope : $μ_scs_slope" + # @debug "μ_scs_intercept: $μ_scs_intercept" σ_cut = (mean(σ) - std(σ) .< σ .< mean(σ) + std(σ)) # fit compton band sigmas with exponential function diff --git a/src/aoefit.jl b/src/aoefit.jl index 57d984fc..61e4a9fc 100644 --- a/src/aoefit.jl +++ b/src/aoefit.jl @@ -276,22 +276,6 @@ function fit_single_aoe_compton(h::Histogram, ps::NamedTuple; uncertainty::Bool= v -> - hist_loglike(x -> x in Interval(extrema(h.edges[1])...) ? f_fit(x, v...) : 0, h) end - # get p-value - mle = f_loglike(v_ml) - - weights_rand = rand(Product(Poisson.(h.weights)), 10000) - - f_loglike_h = let f_fit=f_aoe_compton, v=v_ml - w -> hist_loglike.( - x -> x in Interval(extrema(h.edges[1])...) ? f_fit(x, v) : 0, - fit.(Histogram, Ref(midpoints(h.edges[1])), weights.(w), Ref(h.edges[1])) - ) - end - - mle_rand = f_loglike_h(eachcol(weights_rand)) - - p_value = count(mle_rand .> mle) / length(mle_rand) - if uncertainty # Calculate the Hessian matrix using ForwardDiff H = ForwardDiff.hessian(f_loglike_array, tuple_to_array(v_ml)) @@ -304,10 +288,16 @@ function fit_single_aoe_compton(h::Histogram, ps::NamedTuple; uncertainty::Bool= # Calculate the parameter covariance matrix param_covariance = inv(H) end - + if ~isposdef(param_covariance) + param_covariance = nearestSPD(param_covariance) + end # Extract the parameter uncertainties v_ml_err = array_to_tuple(sqrt.(abs.(diag(param_covariance))), v_ml) + # get p-value + pval, chi2, dof = p_value(f_aoe_compton, h, v_ml) + # calculate normalized residuals + residuals, residuals_norm, _, bin_centers = get_residuals(f_aoe_compton, h, v_ml) @debug "Best Fit values" @debug "μ: $(v_ml.μ) ± $(v_ml_err.μ)" @@ -316,7 +306,8 @@ function fit_single_aoe_compton(h::Histogram, ps::NamedTuple; uncertainty::Bool= @debug "B: $(v_ml.B) ± $(v_ml_err.B)" result = merge(NamedTuple{keys(v_ml)}([measurement(v_ml[k], v_ml_err[k]) for k in keys(v_ml)]...), - (p_value = p_value, )) + (gof = (pvalue = pval, chi2 = chi2, dof = dof, covmat = param_covariance, + residuals = residuals, residuals_norm = residuals_norm, bin_centers = bin_centers),)) else @debug "Best Fit values" @debug "μ: $(v_ml.μ)" @@ -324,7 +315,7 @@ function fit_single_aoe_compton(h::Histogram, ps::NamedTuple; uncertainty::Bool= @debug "n: $(v_ml.n)" @debug "B: $(v_ml.B)" - result = merge(v_ml, (p_value = p_value, )) + result = merge(v_ml, ) end report = ( v = v_ml, diff --git a/src/gof.jl b/src/gof.jl index 9a1be491..d69e7b49 100644 --- a/src/gof.jl +++ b/src/gof.jl @@ -29,7 +29,7 @@ end """ p_value(f_fit::Base.Callable, h::Histogram{<:Real,1},v_ml::NamedTuple) -calculate p-value based on least-squares +calculate p-value based on least-squares, assuming poisson uncertainty baseline method to get goodness-of-fit (gof) # input: * `f_fit`function handle of fit function (peakshape) @@ -57,7 +57,7 @@ function p_value(f_fit::Base.Callable, h::Histogram{<:Real,1}, v_ml::NamedTuple) pval = ccdf(Chisq(dof),chi2) end if any(model_counts .<= 5) - @warn "Bin width <= $(round(minimum(model_counts), digits=0)) counts - Chi2 test might be not valid" + @debug "Bin width <= $(round(minimum(model_counts), digits=0)) counts - Chi2 test might be not valid" else @debug "p-value = $(round(pval, digits=2))" end From 5ad93bafd0b37781a14c830f85ce63065a788911 Mon Sep 17 00:00:00 2001 From: Florian Henkes Date: Wed, 20 Mar 2024 15:07:03 +0100 Subject: [PATCH 14/68] Fixed recipes for calib and fwhm --- ext/LegendSpecFitsRecipesBaseExt.jl | 52 ++++++++++++++++++++++++----- 1 file changed, 43 insertions(+), 9 deletions(-) diff --git a/ext/LegendSpecFitsRecipesBaseExt.jl b/ext/LegendSpecFitsRecipesBaseExt.jl index d4db3132..013edfd2 100644 --- a/ext/LegendSpecFitsRecipesBaseExt.jl +++ b/ext/LegendSpecFitsRecipesBaseExt.jl @@ -468,9 +468,9 @@ end end -@recipe function f(report::NamedTuple{(:par, :f_fit, :x, :y, :gof)}; plot_ribbon=true) +@recipe function f(report::NamedTuple{(:par, :f_fit, :x, :y, :gof)}; plot_ribbon=true, xerrscaling=1) thickness_scaling := 2.0 - xlims := (0, 1.2*maximum(report.x)) + xlims := (0, 1.2*value(maximum(report.x))) framestyle := :box xformatter := :plain yformatter := :plain @@ -480,7 +480,7 @@ end @series begin seriestype := :line subplot := 1 - xticks = :none + xticks --> :none if !isempty(report.gof) label := "Best Fit (p = $(round(report.gof.pvalue, digits=2)))" else @@ -490,16 +490,21 @@ end linewidth := 2 fillalpha := 0.2 if plot_ribbon - ribbon := uncertainty.(report.f_fit.(0:1:1.2*maximum(report.x))) + ribbon := uncertainty.(report.f_fit.(0:1:1.2*value(maximum(report.x)))) end - 0:1:1.2*maximum(report.x), value.(report.f_fit.(0:1:1.2*maximum(report.x))) + 0:1:1.2*value(maximum(report.x)), value.(report.f_fit.(0:1:1.2*value(maximum(report.x)))) end @series begin seriestype := :scatter subplot := 1 - label := "Data" + if xerrscaling == 1 + label := "Data" + else + label := "Data (Error x$(xerrscaling))" + end markercolor --> :black - report.x, report.y + xerror := uncertainty.(report.x) .* xerrscaling + value.(report.x), report.y end @series begin seriestype := :hline @@ -526,10 +531,10 @@ end subplot := 2 label := "" markercolor --> :black - ylabel --> "\n\nResiduals (σ)" + ylabel --> "Residuals (σ)" ylims --> (-5, 5) yticks --> ([-3, 0, 3]) - report.x, report.gof.residuals_norm + value.(report.x), report.gof.residuals_norm end end @@ -560,6 +565,35 @@ end end end +@recipe function f(report::NamedTuple{(:par, :f_fit, :x, :y, :gof, :e_unit, :type)}; xerrscaling=1) + bottom_margin --> (0, :mm) + if report.type == :cal + xlabel := "Energy (ADC)" + legend := :bottomright + framestyle := :box + xlims := (0, 21000) + xticks := (0:2000:22000) + @series begin + grid --> :all + xerrscaling := xerrscaling + (par = report.par, f_fit = report.f_fit, x = report.x, y = report.y, gof = report.gof) + end + @series begin + seriestype := :hline + label := L"Q_{\beta \beta}" + color := :green + fillalpha := 0.2 + linewidth := 2.5 + xticks := :none + ylabel := "Energy ($(report.e_unit))" + ylims := (0, 1.2*value(maximum(report.y))) + yticks := (500:500:3000) + subplot := 1 + [2039] + end + end +end + # @recipe function f(x::AbstractArray, y::AbstractArray{<:Measurement}; plot_ribbon = false) # if plot_ribbon # ribbon := uncertainty.(y) From 10a9468474affe6426aac8a8bf5eba0ed899e712 Mon Sep 17 00:00:00 2001 From: Florian Henkes Date: Wed, 20 Mar 2024 15:07:35 +0100 Subject: [PATCH 15/68] Chi2fit returns x and y values in results --- src/chi2fit.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/chi2fit.jl b/src/chi2fit.jl index da166a00..5d388225 100644 --- a/src/chi2fit.jl +++ b/src/chi2fit.jl @@ -72,7 +72,7 @@ function chi2fit(f_fit::Function, x::AbstractVector{<:Union{Real,Measurement{<:R chi2min = minimum(opt_r) dof = length(x) - length(v_chi2) pvalue = ccdf(Chisq(dof), chi2min) - residuals_norm = (Y_val - f_fit(X_val, v_chi2...)) ./ Y_err + residuals_norm = (Y_val - f_fit(X_val, v_chi2...)) ./ ifelse(all(Y_err .== 0), ones(length(Y_val)), Y_err) result = (par = par, gof = (pvalue = pvalue, chi2min = chi2min, dof = dof, covmat = covmat, residuals_norm = residuals_norm)) report = merge(report, (par = par, gof = result.gof, f_fit = x -> f_fit(x, par...))) end From d4e2d77bfb9ae2881ebd086c346d5d72ad64a832 Mon Sep 17 00:00:00 2001 From: Florian Henkes Date: Wed, 20 Mar 2024 15:07:55 +0100 Subject: [PATCH 16/68] ctc result returns propfunc string for calibration --- src/ctc.jl | 23 +++++++++++++++++------ 1 file changed, 17 insertions(+), 6 deletions(-) diff --git a/src/ctc.jl b/src/ctc.jl index de2105b5..0febbb1a 100644 --- a/src/ctc.jl +++ b/src/ctc.jl @@ -34,9 +34,12 @@ of `window`. The drift time dependence is given by `qdrift`. * `peak`: peak position * `window`: window size * `fct`: correction factor - * `bin_width`: optimal bin width + * `fwhm_before`: FWHM before correction + * `fwhm_after`: FWHM after correction + * `func`: function to correct energy + * `func_generic`: generic function to correct energy """ -function ctc_energy(e::Array{<:Unitful.Energy{<:Real}}, qdrift::Array{<:Real}, peak::Unitful.Energy{<:Real}, window::Tuple{<:Unitful.Energy{<:Real}, <:Unitful.Energy{<:Real}}) +function ctc_energy(e::Array{<:Unitful.Energy{<:Real}}, qdrift::Array{<:Real}, peak::Unitful.Energy{<:Real}, window::Tuple{<:Unitful.Energy{<:Real}, <:Unitful.Energy{<:Real}}, m_cal_simple::Unitful.Energy{<:Real}=1.0u"keV"; e_expression::Union{Symbol, String}="e") # create cut window around peak cut = peak - first(window) .< e .< peak + last(window) e_cut, qdrift_cut = e[cut], qdrift[cut] @@ -65,17 +68,25 @@ function ctc_energy(e::Array{<:Unitful.Energy{<:Real}}, qdrift::Array{<:Real}, p # calculate drift time corrected energy e_ctc = e_cut .+ fct .* qdrift_cut + # get FWHM after correction # fit peak h_after = fit(Histogram, ustrip.(e_unit, e_ctc), ustrip(e_unit, minimum(e_ctc)):ustrip(e_unit, bin_width):ustrip(e_unit, maximum(e_ctc))) ps_after = estimate_single_peak_stats(h_after) result_after, report_after = fit_single_peak_th228(h_after, ps_after; uncertainty=true) + + # get ADC fct factor and cal PropertyFunctions + fct = fct / m_cal_simple + e_ctc_func = "$e_expression + $(mvalue(fct)) * qdrift" + e_ctc_func_generic = "$e_expression + fct * qdrift" + result = ( peak = peak, window = window, + m_cal_simple = m_cal_simple, + func = e_ctc_func, + func_generic = e_ctc_func_generic, fct = fct, - bin_width = bin_width, - bin_width_qdrift = bin_width_qdrift, fwhm_before = result_before.fwhm*e_unit, fwhm_after = result_after.fwhm*e_unit, ) @@ -83,8 +94,8 @@ function ctc_energy(e::Array{<:Unitful.Energy{<:Real}}, qdrift::Array{<:Real}, p peak = result.peak, window = result.window, fct = result.fct, - bin_width = result.bin_width, - bin_width_qdrift = result.bin_width_qdrift, + bin_width = bin_width, + bin_width_qdrift = bin_width_qdrift, e_peak = e_cut, e_ctc = e_ctc, qdrift_peak = qdrift_cut, From d54c07ff0c11b404faddef90595f9084bfc5b862 Mon Sep 17 00:00:00 2001 From: Florian Henkes Date: Wed, 20 Mar 2024 15:08:13 +0100 Subject: [PATCH 17/68] Fixed result issues and report for recipes --- src/fit_calibration.jl | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/src/fit_calibration.jl b/src/fit_calibration.jl index ed55627c..dd28887b 100644 --- a/src/fit_calibration.jl +++ b/src/fit_calibration.jl @@ -9,22 +9,23 @@ Fit the calibration lines with polynomial function of n_poly order * `gof`: godness of fit * `report`: """ -function fit_calibration(n_poly::Int, µ::AbstractVector{<:Union{Real,Measurement{<:Real}}}, peaks::AbstractVector{<:Quantity};e_type::Union{Symbol,String}="x") +function fit_calibration(n_poly::Int, µ::AbstractVector{<:Union{Real,Measurement{<:Real}}}, peaks::AbstractVector{<:Quantity}; e_expression::Union{Symbol, String}="e") @assert length(peaks) == length(μ) e_unit = u"keV" @debug "Fit calibration curve with $(n_poly)-order polynominal function" - result_fit, report = chi2fit(n_poly, µ, ustrip.(e_unit,peaks)) + result_fit, report = chi2fit(n_poly, µ, ustrip.(e_unit, peaks)) par = result_fit.par - par_unit = par .* e_unit - + par_unit = par .* e_unit + result = merge(result_fit, (par = par_unit,)) # built function in string - func = join(["$(mvalue(par[i])) * $(e_type)^$(i-1)" for i=1:length(par)], " + ") - func_generic = join(["p$(i-1) * $(e_type)^$(i-1)" for i=1:length(par)], " + ") - - result = merge(result, (func = func, func_generic = func_generic)) + func = join(["$(mvalue(par[i]))$e_unit * ($(e_expression))^$(i-1)" for i in eachindex(par)], " + ") + func_generic = join(["par[$(i-1)] * ($(e_expression))^$(i-1)" for i in eachindex(par)], " + ") + + result = merge(result, (func = func, func_generic = func_generic, µ = µ, peaks = peaks)) + report = merge(report, (e_unit = e_unit, par = result.par, type = :cal)) return result, report end -export fit_ecal \ No newline at end of file +export fit_calibration \ No newline at end of file From 89cc31c2732e25422d59691b5e8daca3543fa67a Mon Sep 17 00:00:00 2001 From: Florian Henkes Date: Wed, 20 Mar 2024 15:08:42 +0100 Subject: [PATCH 18/68] Fixed fwhm plot report for recipes, adapted results object --- src/fwhm.jl | 8 ++++---- src/specfit.jl | 19 +------------------ 2 files changed, 5 insertions(+), 22 deletions(-) diff --git a/src/fwhm.jl b/src/fwhm.jl index a24706c8..8454ab60 100644 --- a/src/fwhm.jl +++ b/src/fwhm.jl @@ -26,13 +26,13 @@ function fit_fwhm(peaks::Vector{<:Unitful.Energy{<:Real}}, fwhm::Vector{<:Unitfu par_unit = par .* [e_unit^i for i in pol_order:-1:0] # built function in string - func = "sqrt($(join(["$(mvalue(par[i])) * ($(e_expression))^$(i-1)" for i=1:length(par)], " + ")))$e_unit" - func_cal = "sqrt($(join(["$(mvalue(par[i])) * $(e_type_cal)^$(i-1) * keV^$(length(par)+1-i)" for i=1:length(par)], " + ")))" - func_generic = "sqrt($(join(["par[$(i-1)] * $(e_type_cal)^$(i-1)" for i=1:length(par)], " + ")))" + func = "sqrt($(join(["$(mvalue(par[i])) * ($(e_expression))^$(i-1)" for i in eachindex(length(par))], " + ")))$e_unit" + func_cal = "sqrt($(join(["$(mvalue(par[i])) * $(e_type_cal)^$(i-1) * keV^$(length(par)+1-i)" for i in eachindex(length(par))], " + ")))" + func_generic = "sqrt($(join(["par[$(i-1)] * $(e_type_cal)^$(i-1)" for i in eachindex(length(par))], " + ")))" # get fwhm at Qbb qbb = report_chi2.f_fit(2039) * e_unit - result = merge(result_chi2, (par = par_unit , qbb = qbb, func = func, func_cal = func_cal, func_generic = func_generic)) + result = merge(result_chi2, (par = par_unit , qbb = qbb, func = func, func_cal = func_cal, func_generic = func_generic, peaks = peaks, fwhm = fwhm)) report = merge(report_chi2, (par = result.par, qbb = result.qbb, type = :fwhm)) return result, report diff --git a/src/specfit.jl b/src/specfit.jl index 4d646925..ca5b110c 100644 --- a/src/specfit.jl +++ b/src/specfit.jl @@ -321,21 +321,4 @@ function get_peak_fwhm_th228(v_ml::NamedTuple, v_ml_err::Union{Matrix,NamedTuple fwhm_err = std(fwhm_mc[isfinite.(fwhm_mc)]) return fwhm, fwhm_err end -export get_peak_fwhm_th228 - -""" - fit_calibration(μ::Vector{<:Real}, peaks::Vector{<:Unitful.Energy{<:Real}}) -Fit the calibration lines to a linear function. -# Returns - * `slope`: the slope of the linear fit - * `intercept`: the intercept of the linear fit -""" -function fit_calibration(μ::Vector{<:Real}, peaks::Vector{<:Unitful.Energy{<:Real}}) - @assert length(peaks) == length(μ) - @debug "Fit calibration curve with linear function" - peaks_unit = unit(first(peaks)) - calib_fit_result = linregress(μ, ustrip.(peaks_unit, peaks)) - return LinearRegression.slope(calib_fit_result)[1]*peaks_unit, LinearRegression.bias(calib_fit_result)[1]*peaks_unit -end -export fit_calibration -fit_calibration(μ::Vector{<:Unitful.Energy{<:Real}}, peaks::Vector{<:Unitful.Energy{<:Real}}) = fit_calibration(ustrip.(μ), peaks) \ No newline at end of file +export get_peak_fwhm_th228 \ No newline at end of file From 3d16d0214cea134df8e45b69c5bfbf512eba3fe8 Mon Sep 17 00:00:00 2001 From: Florian Henkes Date: Wed, 20 Mar 2024 17:04:43 +0100 Subject: [PATCH 19/68] CTC result improvements in output --- src/ctc.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/ctc.jl b/src/ctc.jl index 0febbb1a..dce7bdb8 100644 --- a/src/ctc.jl +++ b/src/ctc.jl @@ -101,8 +101,8 @@ function ctc_energy(e::Array{<:Unitful.Energy{<:Real}}, qdrift::Array{<:Real}, p qdrift_peak = qdrift_cut, h_before = h_before, h_after = h_after, - fwhm_before = result_before.fwhm, - fwhm_after = result_after.fwhm, + fwhm_before = result.fwhm_before, + fwhm_after = result.fwhm_after, report_before = report_before, report_after = report_after ) From b67167f0379d270e670b785dec76320582aca57c Mon Sep 17 00:00:00 2001 From: Florian Henkes Date: Wed, 20 Mar 2024 17:04:59 +0100 Subject: [PATCH 20/68] Fixed ctc plot recipe --- ext/LegendSpecFitsRecipesBaseExt.jl | 61 ++++++++++++++++++----------- 1 file changed, 39 insertions(+), 22 deletions(-) diff --git a/ext/LegendSpecFitsRecipesBaseExt.jl b/ext/LegendSpecFitsRecipesBaseExt.jl index 013edfd2..ed5e8687 100644 --- a/ext/LegendSpecFitsRecipesBaseExt.jl +++ b/ext/LegendSpecFitsRecipesBaseExt.jl @@ -337,17 +337,26 @@ end end @recipe function f(report_ctc::NamedTuple{(:peak, :window, :fct, :bin_width, :bin_width_qdrift, :e_peak, :e_ctc, :qdrift_peak, :h_before, :h_after, :fwhm_before, :fwhm_after, :report_before, :report_after)}) - layout := (2,2) - thickness_scaling := 2.0 - size := (2400, 1600) + layout := (1, 3) + thickness_scaling := 1.0 + xtickfontsize := 12 + xlabelfontsize := 14 + ylabelfontsize := 14 + ytickfontsize := 12 + legendfontsize := 10 + size := (1000, 300) + margin := (8, :mm) @series begin seriestype := :histogram2d bins := (ustrip.(unit(first(report_ctc.e_peak)), minimum(report_ctc.e_peak):report_ctc.bin_width:maximum(report_ctc.e_peak)), quantile(report_ctc.qdrift_peak, 0.01):report_ctc.bin_width_qdrift:quantile(report_ctc.qdrift_peak, 0.99)) color := :inferno - xlabel := "Energy (keV)" - ylabel := "Drift Time" + xlabel := "Energy" + ylabel := "QDrift" title := "Before Correction" + titlelocation := (0.5, 1.1) xlims := (2600, 2630) + ylims := (0, quantile(report_ctc.qdrift_peak, 0.99)) + yformatter := :plain legend := :none colorbar_scale := :log10 subplot := 1 @@ -357,27 +366,31 @@ end seriestype := :histogram2d bins := (ustrip.(unit(first(report_ctc.e_peak)), minimum(report_ctc.e_peak):report_ctc.bin_width:maximum(report_ctc.e_peak)), quantile(report_ctc.qdrift_peak, 0.01):report_ctc.bin_width_qdrift:quantile(report_ctc.qdrift_peak, 0.99)) color := :magma - xlabel := "Energy (keV)" - ylabel := "Drift Time" + xlabel := "Energy" + ylabel := "QDrift" title := "After Correction" xlims := (2600, 2630) + titlelocation := (0.5, 1.1) + xlims := (2600, 2630) + ylims := (0, quantile(report_ctc.qdrift_peak, 0.99)) + yformatter := :plain legend := :none colorbar_scale := :log10 - subplot := 2 - report_ctc.e_ctc, report_ctc.qdrift_peak - end - @series begin - seriestype := :stepbins - color := :red - label := "Before CTC" - xlabel := "Energy (keV)" - ylabel := "Counts" - title := "FWHM $(round(report_ctc.fwhm_before, digits=2))" - yscale := :log10 subplot := 3 - report_ctc.h_before + report_ctc.e_ctc, report_ctc.qdrift_peak end # @series begin + # seriestype := :stepbins + # color := :red + # label := "Before CTC" + # xlabel := "Energy (keV)" + # ylabel := "Counts" + # title := "FWHM $(round(report_ctc.fwhm_before, digits=2))" + # yscale := :log10 + # subplot := 3 + # report_ctc.h_before + # end + # @series begin # # seriestype := :stepbins # color := :red # # label := "Before CTC" @@ -408,7 +421,7 @@ end xlabel := "Energy (keV)" ylabel := "Counts" yscale := :log10 - subplot := 4 + subplot := 2 report_ctc.h_before end @series begin @@ -417,9 +430,13 @@ end label := "After CTC" xlabel := "Energy (keV)" ylabel := "Counts" - title := "FWHM $(round(report_ctc.fwhm_after, digits=2))" + title := "FWHM $(round(u"keV", report_ctc.fwhm_after, digits=2))" + titlelocation := (0.5, 1.1) + xlims := (2600, 2630) + xticks := (2600:10:2630) + legend := :bottomright yscale := :log10 - subplot := 4 + subplot := 2 report_ctc.h_after end end From b1c73a903eea62e72d949256fba107a197e25c3e Mon Sep 17 00:00:00 2001 From: LisaSchlueter Date: Thu, 21 Mar 2024 19:09:12 +0100 Subject: [PATCH 21/68] update aoe correction/normalization fit routines with chi2fit and ljl_probfunc --- src/aoe_calibration.jl | 59 ++++++++++++++++++++++-------------------- 1 file changed, 31 insertions(+), 28 deletions(-) diff --git a/src/aoe_calibration.jl b/src/aoe_calibration.jl index 2fd6290c..1ceae0a8 100644 --- a/src/aoe_calibration.jl +++ b/src/aoe_calibration.jl @@ -1,8 +1,6 @@ # f_aoe_sigma(x, p) = p[1] .+ p[2]*exp.(-p[3]./x) @. f_aoe_sigma(x, p) = sqrt(abs(p[1]) + abs(p[2])/x^2) f_aoe_mu(x, p) = p[1] .+ p[2].*x - - """ fit_aoe_corrections(e::Array{<:Unitful.Energy{<:Real}}, μ::Array{<:Real}, σ::Array{<:Real}) @@ -16,40 +14,45 @@ Fit the corrections for the AoE value of the detector. - `σ_scs`: Fit result for the sigma values - `f_σ_scs`: Fit function for the sigma values """ -function fit_aoe_corrections(e::Array{<:Unitful.Energy{<:Real}}, μ::Array{<:Real}, σ::Array{<:Real}) +function fit_aoe_corrections(e::Array{<:Unitful.Energy{<:Real}}, μ::Array{<:Real}, σ::Array{<:Real}; e_expression::Union{String,Symbol}="e") # fit compton band mus with linear function μ_cut = (mean(μ) - 2*std(μ) .< μ .< mean(μ) + 2*std(μ)) e, μ, σ = e[μ_cut], μ[μ_cut], σ[μ_cut] e_unit = unit(first(e)) e = ustrip.(e_unit, e) - # fit compton band with linear function - result_µ, report_µ = LegendSpecFits.chi2fit(1, e, µ; uncertainty=true) -result_µ - - # μ_scs = linregress(e, μ) - # μ_scs_slope, μ_scs_intercept = LinearRegression.slope(μ_scs)[1], LinearRegression.bias(μ_scs)[1] - # μ_scs = curve_fit(f_aoe_mu, e, μ, [μ_scs_intercept, μ_scs_slope]) - # μ_scs_param = [μ_scs.param[1], μ_scs.param[2] / e_unit] - # @debug "μ_scs_slope : $μ_scs_slope" - # @debug "μ_scs_intercept: $μ_scs_intercept" + # fit compton band µ with linear function + result_µ, report_µ = chi2fit(1, e, µ; uncertainty=true) + func_µ = "$(mvalue(result_µ.par[1])) + ($e_expression) * $(mvalue(result_µ.par[2]))$e_unit^-1" + func_generic_µ = "p[1] + ($e_expression) * p[2]" + par_µ = [result_µ.par[i] ./ e_unit^(i-1) for i=1:length(result_µ.par)] # add units + result_µ = merge(result_µ, (par = par_µ, func = func_µ, func_generic = func_generic_µ, µ = µ)) + @debug "Compton band µ correction: $(result_µ.func)" + # fit compton band σ with sqrt function σ_cut = (mean(σ) - std(σ) .< σ .< mean(σ) + std(σ)) - # fit compton band sigmas with exponential function - σ_scs = curve_fit(f_aoe_sigma, e[σ_cut], σ[σ_cut], [median(σ)^2, 1]) - σ_scs_param = [σ_scs.param[1], σ_scs.param[2] * e_unit^2] - @debug "σ_scs offset: $(σ_scs_param[1])" - @debug "σ_scs shift : $(σ_scs_param[2])" - - ( - e = e .* e_unit, - μ = μ, - μ_scs = μ_scs_param, - f_μ_scs = x -> Base.Fix2(f_aoe_mu, μ_scs_param).(x), - σ = σ, - σ_scs = σ_scs_param, - f_σ_scs = x -> Base.Fix2(f_aoe_sigma, σ_scs_param).(x), - ) + f_fit_σ = f_aoe_sigma # fit function + result_σ, report_σ = chi2fit((x, p1, p2) -> f_fit_σ(x,[p1,p2]), e[σ_cut], µ[σ_cut]; uncertainty=true) + par_σ = [result_σ.par[1], result_σ.par[2] * e_unit^2] # add unit + func_σ = nothing + func_generic_σ = nothing + if string(f_fit_σ) == "f_aoe_sigma" + func_σ = "sqrt( abs($(mvalue(result_σ.par[1]))) + abs($(mvalue(result_σ.par[2]))$(e_unit)^2) / ($e_expression)^2)" # add units!! + func_generic_σ = "sqrt(abs(p[1]) + abs(p[2]) / ($e_expression)^2)" + end + result_σ = merge(result_σ, (par = par_σ, func = func_σ, func_generic = func_generic_σ, σ = σ)) + @debug "Compton band σ normalization: $(result_σ.func)" + + # put everything together into A/E correction/normalization function + aoe_str = "(a / (($e_expression)$e_unit^-1))" # get aoe, but without unit. + func_aoe_corr = "($aoe_str - ($(result_µ.func))) / ($(result_σ.func))" + func_generic_aoe_corr = "(aoe - $(result_µ.func_generic)) / $(result_σ.func_generic)" + func_aoe_corr_ecal = replace(func_aoe_corr,e_expression => "e_cal") # function that can be used for already calibrated energies + + result = (µ_compton = result_µ, σ_compton = result_σ, compton_bands = (e = e,), func = func_aoe_corr, func_generic = func_generic_aoe_corr, func_ecal = func_aoe_corr_ecal) + report = (report_µ = report_µ, report_σ = report_σ) + + return result, report end export fit_aoe_corrections From 33bf45951c79a299be2b280e6bcb4b2b771ed13e Mon Sep 17 00:00:00 2001 From: LisaSchlueter Date: Thu, 21 Mar 2024 19:45:57 +0100 Subject: [PATCH 22/68] deactivate chi2fit test until new LegendDataManagement release --- test/runtests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index 839ceaad..822649c6 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -5,7 +5,7 @@ import Test Test.@testset "Package LegendSpecFits" begin #include("test_aqua.jl") include("test_specfit.jl") - include("test_fit_chisq.jl") + # include("test_fit_chisq.jl") include("test_docs.jl") isempty(Test.detect_ambiguities(LegendSpecFits)) end # testset From 85e0c085e982ef1bc9f8dcf4bcb6810ed348115d Mon Sep 17 00:00:00 2001 From: LisaSchlueter Date: Mon, 25 Mar 2024 12:11:57 +0100 Subject: [PATCH 23/68] bug fix: fit mu instead of sigma in compton band fits --- src/aoe_calibration.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/aoe_calibration.jl b/src/aoe_calibration.jl index 1ceae0a8..8d36a1c1 100644 --- a/src/aoe_calibration.jl +++ b/src/aoe_calibration.jl @@ -24,15 +24,16 @@ function fit_aoe_corrections(e::Array{<:Unitful.Energy{<:Real}}, μ::Array{<:Rea # fit compton band µ with linear function result_µ, report_µ = chi2fit(1, e, µ; uncertainty=true) func_µ = "$(mvalue(result_µ.par[1])) + ($e_expression) * $(mvalue(result_µ.par[2]))$e_unit^-1" - func_generic_µ = "p[1] + ($e_expression) * p[2]" + func_generic_µ = "p1 + ($e_expression) * p2" # "p[1] + ($e_expression) * p[2]" par_µ = [result_µ.par[i] ./ e_unit^(i-1) for i=1:length(result_µ.par)] # add units result_µ = merge(result_µ, (par = par_µ, func = func_µ, func_generic = func_generic_µ, µ = µ)) + # report_µ = merge(report_µ, (µ_fit = ustrip.(result_µ.par[1]) .+ ustrip.(result_µ.par[2] .* e), e_fit = e*e_unit)) @debug "Compton band µ correction: $(result_µ.func)" - + # fit compton band σ with sqrt function σ_cut = (mean(σ) - std(σ) .< σ .< mean(σ) + std(σ)) f_fit_σ = f_aoe_sigma # fit function - result_σ, report_σ = chi2fit((x, p1, p2) -> f_fit_σ(x,[p1,p2]), e[σ_cut], µ[σ_cut]; uncertainty=true) + result_σ, report_σ = chi2fit((x, p1, p2) -> f_fit_σ(x,[p1,p2]), e[σ_cut], σ[σ_cut]; uncertainty=true) par_σ = [result_σ.par[1], result_σ.par[2] * e_unit^2] # add unit func_σ = nothing func_generic_σ = nothing From d00e6c9e8b9f36df59307eeed3644eadda01c02c Mon Sep 17 00:00:00 2001 From: LisaSchlueter Date: Mon, 25 Mar 2024 13:52:57 +0100 Subject: [PATCH 24/68] add plot receipe for compton band fits --- ext/LegendSpecFitsRecipesBaseExt.jl | 95 +++++++++++++++++++++-------- src/aoe_calibration.jl | 5 +- 2 files changed, 74 insertions(+), 26 deletions(-) diff --git a/ext/LegendSpecFitsRecipesBaseExt.jl b/ext/LegendSpecFitsRecipesBaseExt.jl index 3d7bc840..e172df72 100644 --- a/ext/LegendSpecFitsRecipesBaseExt.jl +++ b/ext/LegendSpecFitsRecipesBaseExt.jl @@ -611,30 +611,77 @@ end end end -# @recipe function f(x::AbstractArray, y::AbstractArray{<:Measurement}; plot_ribbon = false) -# if plot_ribbon -# ribbon := uncertainty.(y) -# end -# x, value.(y) -# end +@recipe function f(report::NamedTuple{(:par, :f_fit, :x, :y, :gof, :e_unit, :label_y, :label_fit)}) + layout --> @layout([a{0.8h}; b{0.2h}]) + margins --> (0, :mm) + link --> :x + size := (1200, 700) + layout --> @layout([a{0.7h}; b{0.3h}]) + xmin = floor(Int, minimum(report.x)/100)*100 + xmax = ceil(Int, maximum(report.x)/100)*100 + @series begin + seriestype := :line + subplot --> 1 + color := :orange + ms := 3 + linewidth := 3 + label := report.label_fit + report.x, report.f_fit(report.x) + end + @series begin + ylabel := "$(report.label_y) (a.u.)" + seriestype := :scatter + subplot --> 1 + color := :black + ms := 3 + yguidefontsize := 18 + xguidefontsize := 18 + ytickfontsize := 12 + xtickfontsize := 12 + legendfontsize := 14 + foreground_color_legend := :silver + background_color_legend := :white + xlims := (xmin,xmax) + xticks := (xmin:250:xmax, fill(" ", length(xmin:250:xmax) )) + ylims := (0.98 * (Measurements.value(minimum(report.y)) - Measurements.uncertainty(median(report.y))), 1.02 * (Measurements.value(maximum(report.y)) + Measurements.uncertainty(median(report.y)) ) ) + margin := (10, :mm) + bottom_margin := (-7, :mm) + framestyle := :box + grid := :false + label := "Compton band fits: Gaussian $(report.label_y)(A/E)" + report.x, report.y + end -# @recipe function f(covmat::Matrix{<:Real};mode::Symbol=:cov) -# if mode==:cov -# colormap = :viridis -# title_str = "covariance" -# elseif mode==:corr -# cm = cor(cm); -# colormap = :greys -# title_str = "correlation" -# end -# @series begin -# seriestype := :heatmap -# yflip := true -# colorbar_title_location := :right -# guidefontsize := 16 -# tickfontsize := 14 -# size := (1000, 600) -# end -# end + @series begin + seriestype := :hline + ribbon := 3 + subplot --> 2 + fillalpha := 0.5 + label := "" + fillcolor := :lightgrey + linecolor := :darkgrey + [0.0] + end + @series begin + xlabel := "Energy ($(report.e_unit))" + ylabel := "Residuals (σ) \n" + seriestype := :scatter + subplot --> 2 + color := :black + ms := 3 + label := false + framestyle := :box + grid := :false + xlims := (xmin,xmax) + xticks := xmin:250:xmax + ylims := (floor(minimum(report.gof.residuals_norm)-1), ceil(maximum(report.gof.residuals_norm))+1) + yticks := [-5,0,5] + yguidefontsize := 18 + xguidefontsize := 18 + ytickfontsize := 12 + xtickfontsize := 12 + report.x, report.gof.residuals_norm + end +end end # module LegendSpecFitsRecipesBaseExt \ No newline at end of file diff --git a/src/aoe_calibration.jl b/src/aoe_calibration.jl index 8d36a1c1..21818f37 100644 --- a/src/aoe_calibration.jl +++ b/src/aoe_calibration.jl @@ -27,7 +27,7 @@ function fit_aoe_corrections(e::Array{<:Unitful.Energy{<:Real}}, μ::Array{<:Rea func_generic_µ = "p1 + ($e_expression) * p2" # "p[1] + ($e_expression) * p[2]" par_µ = [result_µ.par[i] ./ e_unit^(i-1) for i=1:length(result_µ.par)] # add units result_µ = merge(result_µ, (par = par_µ, func = func_µ, func_generic = func_generic_µ, µ = µ)) - # report_µ = merge(report_µ, (µ_fit = ustrip.(result_µ.par[1]) .+ ustrip.(result_µ.par[2] .* e), e_fit = e*e_unit)) + report_µ = merge(report_µ, (e_unit = e_unit, label_y = "µ", label_fit = "Best Fit: $(mvalue(round(result_µ.par[1], digits=2))) + E * $(mvalue(round(ustrip(result_µ.par[2]) * 1e6, digits=2)))1e-6")) @debug "Compton band µ correction: $(result_µ.func)" # fit compton band σ with sqrt function @@ -42,7 +42,8 @@ function fit_aoe_corrections(e::Array{<:Unitful.Energy{<:Real}}, μ::Array{<:Rea func_generic_σ = "sqrt(abs(p[1]) + abs(p[2]) / ($e_expression)^2)" end result_σ = merge(result_σ, (par = par_σ, func = func_σ, func_generic = func_generic_σ, σ = σ)) - @debug "Compton band σ normalization: $(result_σ.func)" + report_σ = merge(report_σ, (e_unit = e_unit, label_y = "σ", label_fit = "Best fit: sqrt($(round(mvalue(result_σ.par[1])*1e6, digits=1))e-6 + $(round(ustrip(mvalue(result_σ.par[2])), digits=2)) / E^2)")) + @debug "Compton band σ normalization: $(result_σ.func)" # put everything together into A/E correction/normalization function aoe_str = "(a / (($e_expression)$e_unit^-1))" # get aoe, but without unit. From e4da68fd6f832e8f99ecd6c412e0d06b01873b63 Mon Sep 17 00:00:00 2001 From: Florian Henkes Date: Mon, 25 Mar 2024 21:12:54 +0100 Subject: [PATCH 25/68] Removed correct_aoe --> Utils --- src/aoe_calibration.jl | 25 +++++-------------------- 1 file changed, 5 insertions(+), 20 deletions(-) diff --git a/src/aoe_calibration.jl b/src/aoe_calibration.jl index 21818f37..ed3315ca 100644 --- a/src/aoe_calibration.jl +++ b/src/aoe_calibration.jl @@ -29,7 +29,7 @@ function fit_aoe_corrections(e::Array{<:Unitful.Energy{<:Real}}, μ::Array{<:Rea result_µ = merge(result_µ, (par = par_µ, func = func_µ, func_generic = func_generic_µ, µ = µ)) report_µ = merge(report_µ, (e_unit = e_unit, label_y = "µ", label_fit = "Best Fit: $(mvalue(round(result_µ.par[1], digits=2))) + E * $(mvalue(round(ustrip(result_µ.par[2]) * 1e6, digits=2)))1e-6")) @debug "Compton band µ correction: $(result_µ.func)" - + # fit compton band σ with sqrt function σ_cut = (mean(σ) - std(σ) .< σ .< mean(σ) + std(σ)) f_fit_σ = f_aoe_sigma # fit function @@ -40,10 +40,10 @@ function fit_aoe_corrections(e::Array{<:Unitful.Energy{<:Real}}, μ::Array{<:Rea if string(f_fit_σ) == "f_aoe_sigma" func_σ = "sqrt( abs($(mvalue(result_σ.par[1]))) + abs($(mvalue(result_σ.par[2]))$(e_unit)^2) / ($e_expression)^2)" # add units!! func_generic_σ = "sqrt(abs(p[1]) + abs(p[2]) / ($e_expression)^2)" - end - result_σ = merge(result_σ, (par = par_σ, func = func_σ, func_generic = func_generic_σ, σ = σ)) - report_σ = merge(report_σ, (e_unit = e_unit, label_y = "σ", label_fit = "Best fit: sqrt($(round(mvalue(result_σ.par[1])*1e6, digits=1))e-6 + $(round(ustrip(mvalue(result_σ.par[2])), digits=2)) / E^2)")) - @debug "Compton band σ normalization: $(result_σ.func)" + end + result_σ = merge(result_σ, (par = par_σ, func = func_σ, func_generic = func_generic_σ, σ = σ)) + report_σ = merge(report_σ, (e_unit = e_unit, label_y = "σ", label_fit = "Best fit: sqrt($(round(mvalue(result_σ.par[1])*1e6, digits=1))e-6 + $(round(ustrip(mvalue(result_σ.par[2])), digits=2)) / E^2)")) + @debug "Compton band σ normalization: $(result_σ.func)" # put everything together into A/E correction/normalization function aoe_str = "(a / (($e_expression)$e_unit^-1))" # get aoe, but without unit. @@ -58,21 +58,6 @@ function fit_aoe_corrections(e::Array{<:Unitful.Energy{<:Real}}, μ::Array{<:Rea end export fit_aoe_corrections -""" - correctaoe!(aoe::Array{T}, e::Array{T}, aoe_corrections::NamedTuple) where T<:Real - -Correct the AoE values in the `aoe` array using the corrections in `aoe_corrections`. -""" -function correct_aoe!(aoe::Array{<:Unitful.RealOrRealQuantity}, e::Array{<:Unitful.Energy{<:Real}}, aoe_corrections::NamedTuple) - aoe .-= Base.Fix2(f_aoe_mu, mvalue.(aoe_corrections.μ_scs)).(e) - # aoe .-= 1.0*unit(aoe[1]) - aoe ./= Base.Fix2(f_aoe_sigma, mvalue.(aoe_corrections.σ_scs)).(e) - ustrip.(aoe) -end -export correct_aoe! - -correct_aoe!(aoe, e, aoe_corrections::PropDict) = correct_aoe!(aoe, e, NamedTuple(aoe_corrections)) - """ prepare_dep_peakhist(e::Array{T}, dep::T,; relative_cut::T=0.5, n_bins_cut::Int=500) where T<:Real From d4be2b9fea922364b235e413d9d4dcd6956f5e63 Mon Sep 17 00:00:00 2001 From: Florian Henkes Date: Mon, 25 Mar 2024 21:13:45 +0100 Subject: [PATCH 26/68] Update types, include assert to check that aoe and e have same length --- src/aoefit.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/aoefit.jl b/src/aoefit.jl index 61e4a9fc..41579639 100644 --- a/src/aoefit.jl +++ b/src/aoefit.jl @@ -65,7 +65,7 @@ function estimate_single_peak_stats_psd(h::Histogram{T}) where T<:Real end """ - generate_aoe_compton_bands(aoe::Array{<:Real}, e::Array{<:Real}, compton_bands::Array{<:Real}, compton_window::T) where T<:Real + generate_aoe_compton_bands(aoe::Vector{<:Real}, e::Vector{<:T}, compton_bands::Vector{<:T}, compton_window::T) where T<:Unitful.Energy{<:Real} Generate histograms for the A/E Compton bands and estimate peak parameters. The compton bands are cutted out of the A/E spectrum and then binned using the Freedman-Diaconis Rule. For better performance @@ -83,7 +83,8 @@ the binning is only done in the area around the peak. The peak parameters are es * `simple_pars_aoe_σ`: Simple curve fit parameters for the peak sigma energy depencence * `simple_pars_error_aoe_σ`: Simple curve fit parameter errors for the peak sigma energy depencence """ -function generate_aoe_compton_bands(aoe::Array{<:Real}, e::Array{<:T}, compton_bands::Array{<:T}, compton_window::T) where T<:Unitful.Energy{<:Real} +function generate_aoe_compton_bands(aoe::Vector{<:Real}, e::Vector{<:T}, compton_bands::Vector{<:T}, compton_window::T) where T<:Unitful.Energy{<:Real} + @assert length(aoe) == length(e) "A/E and Energy arrays must have the same length" e_unit = u"keV" # get aoe values in compton bands aoe_compton_bands = [aoe[c .< e .< c + compton_window .&& aoe .> 0.0] for c in compton_bands] From a881c8ad5d268a77472e942714cb154e0860fb5f Mon Sep 17 00:00:00 2001 From: Florian Henkes Date: Tue, 26 Mar 2024 01:37:01 +0100 Subject: [PATCH 27/68] Bug Fix sigma and mu recipe in compton band fits --- ext/LegendSpecFitsRecipesBaseExt.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/ext/LegendSpecFitsRecipesBaseExt.jl b/ext/LegendSpecFitsRecipesBaseExt.jl index e172df72..c4403b9e 100644 --- a/ext/LegendSpecFitsRecipesBaseExt.jl +++ b/ext/LegendSpecFitsRecipesBaseExt.jl @@ -626,10 +626,11 @@ end ms := 3 linewidth := 3 label := report.label_fit - report.x, report.f_fit(report.x) + ribbon := uncertainty.(report.f_fit(report.x)) + report.x, value.(report.f_fit(report.x)) end @series begin - ylabel := "$(report.label_y) (a.u.)" + ylabel := "$(report.label_y) (a.u.)" seriestype := :scatter subplot --> 1 color := :black From b394e43c68def177f324e6704e4df82c7672a2d6 Mon Sep 17 00:00:00 2001 From: Florian Henkes Date: Tue, 26 Mar 2024 11:10:01 +0100 Subject: [PATCH 28/68] Fix problem in case of compton corrections chi2sqit with 0.0 uncerts --- src/aoe_calibration.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/aoe_calibration.jl b/src/aoe_calibration.jl index ed3315ca..c9389405 100644 --- a/src/aoe_calibration.jl +++ b/src/aoe_calibration.jl @@ -16,7 +16,7 @@ Fit the corrections for the AoE value of the detector. """ function fit_aoe_corrections(e::Array{<:Unitful.Energy{<:Real}}, μ::Array{<:Real}, σ::Array{<:Real}; e_expression::Union{String,Symbol}="e") # fit compton band mus with linear function - μ_cut = (mean(μ) - 2*std(μ) .< μ .< mean(μ) + 2*std(μ)) + μ_cut = (mean(μ) - 2*std(μ) .< μ .< mean(μ) + 2*std(μ)) .&& muncert.(μ) .> 0.0 e, μ, σ = e[μ_cut], μ[μ_cut], σ[μ_cut] e_unit = unit(first(e)) e = ustrip.(e_unit, e) @@ -31,7 +31,7 @@ function fit_aoe_corrections(e::Array{<:Unitful.Energy{<:Real}}, μ::Array{<:Rea @debug "Compton band µ correction: $(result_µ.func)" # fit compton band σ with sqrt function - σ_cut = (mean(σ) - std(σ) .< σ .< mean(σ) + std(σ)) + σ_cut = (mean(σ) - std(σ) .< σ .< mean(σ) + std(σ)) .&& muncert.(σ) .> 0.0 f_fit_σ = f_aoe_sigma # fit function result_σ, report_σ = chi2fit((x, p1, p2) -> f_fit_σ(x,[p1,p2]), e[σ_cut], σ[σ_cut]; uncertainty=true) par_σ = [result_σ.par[1], result_σ.par[2] * e_unit^2] # add unit From 4ec1619380872454aa0dad4b41672d7ea45d8356 Mon Sep 17 00:00:00 2001 From: LisaSchlueter Date: Wed, 3 Apr 2024 19:51:57 +0200 Subject: [PATCH 29/68] bug fix chi2fit/normalized residuals did not xerr into account. --- src/chi2fit.jl | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/src/chi2fit.jl b/src/chi2fit.jl index 5d388225..ce5342e9 100644 --- a/src/chi2fit.jl +++ b/src/chi2fit.jl @@ -72,7 +72,16 @@ function chi2fit(f_fit::Function, x::AbstractVector{<:Union{Real,Measurement{<:R chi2min = minimum(opt_r) dof = length(x) - length(v_chi2) pvalue = ccdf(Chisq(dof), chi2min) - residuals_norm = (Y_val - f_fit(X_val, v_chi2...)) ./ ifelse(all(Y_err .== 0), ones(length(Y_val)), Y_err) + 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...) + y_pred_err = only(ForwardDiff.partials(UncertTag, dual_y)) + return y_pred_err + end + Y_pred_err = get_y_pred_err.(Ref(f_fit), X_val, X_err, Ref(v_chi2)) + Y_err_tot = sqrt.(Y_err.^2 .+ Y_pred_err.^2) + Y_err_tot = ifelse(all(Y_err_tot .== 0), ones(length(Y_val)), Y_err_tot) + residuals_norm = (Y_val - f_fit(X_val, v_chi2...)) ./ Y_err_tot result = (par = par, gof = (pvalue = pvalue, chi2min = chi2min, dof = dof, covmat = covmat, residuals_norm = residuals_norm)) report = merge(report, (par = par, gof = result.gof, f_fit = x -> f_fit(x, par...))) end From dc01c52b8c50c1bb12903144801f5b0bbaa52827 Mon Sep 17 00:00:00 2001 From: Florian Henkes Date: Thu, 4 Apr 2024 11:25:42 +0200 Subject: [PATCH 30/68] Fixed plot errors for f_fit plots with fixed steprangelens steps --- ext/LegendSpecFitsRecipesBaseExt.jl | 26 +++++++++++++++++--------- 1 file changed, 17 insertions(+), 9 deletions(-) diff --git a/ext/LegendSpecFitsRecipesBaseExt.jl b/ext/LegendSpecFitsRecipesBaseExt.jl index c4403b9e..e2fdb7e0 100644 --- a/ext/LegendSpecFitsRecipesBaseExt.jl +++ b/ext/LegendSpecFitsRecipesBaseExt.jl @@ -7,6 +7,14 @@ using Unitful, Formatting, Measurements, LaTeXStrings using Measurements: value, uncertainty using StatsBase, LinearAlgebra +function round_wo_units(x::Unitful.RealOrRealQuantity, digits::Integer=2) + if unit(x) == NoUnits + round(x, digits=digits) + else + round(unit(x), x, digits=2) + end +end + # @recipe function f(x::Vector{T}, cuts::NamedTuple{(:low, :high, :max), Tuple{T, T, T}}) where T<:Unitful.RealOrRealQuantity @recipe function f(report::NamedTuple{(:f_fit, :μ, :σ, :n)}, x::Vector{T}, cuts::NamedTuple{(:low, :high, :max), Tuple{T, T, T}}) where {T <: Unitful.RealOrRealQuantity} ylabel := "Normalized Counts" @@ -20,7 +28,7 @@ using StatsBase, LinearAlgebra end @series begin color := :red - label := "Normal Fit (μ = $(round(unit(report.μ), report.μ, digits=2)), σ = $(round(unit(report.σ), report.σ, digits=2)))" + label := "Normal Fit (μ = $(round_wo_units(report.μ, digits=2)), σ = $(round_wo_units(report.σ, digits=2)))" lw := 3 ustrip(cuts.low):ustrip(Measurements.value(report.σ / 1000)):ustrip(cuts.high), t -> report.f_fit(t) end @@ -134,21 +142,21 @@ end color := :red subplot --> _subplot # ribbon := uncertainty.(report.f_fit.(minimum(report.h.edges[1]):0.1:maximum(report.h.edges[1]))) - minimum(report.h.edges[1]):0.1:maximum(report.h.edges[1]), value.(report.f_fit.(minimum(report.h.edges[1]):0.1:maximum(report.h.edges[1]))) + minimum(report.h.edges[1]):ustrip(value(report.v.σ))/1000:maximum(report.h.edges[1]), value.(report.f_fit.(minimum(report.h.edges[1]):0.1:maximum(report.h.edges[1]))) end @series begin seriestype := :line label := ifelse(show_label, "Signal", "") subplot --> _subplot color := :green - minimum(report.h.edges[1]):0.1:maximum(report.h.edges[1]), report.f_sig + minimum(report.h.edges[1]):ustrip(value(report.v.σ))/1000:maximum(report.h.edges[1]), report.f_sig end @series begin seriestype := :line label := ifelse(show_label, "Low Energy Tail", "") subplot --> _subplot color := :blue - minimum(report.h.edges[1]):0.1:maximum(report.h.edges[1]), report.f_lowEtail + minimum(report.h.edges[1]):ustrip(value(report.v.σ))/1000:maximum(report.h.edges[1]), report.f_lowEtail end @series begin seriestype := :line @@ -162,7 +170,7 @@ end ylims := (ylim_min, ylim_max) xlims := (minimum(report.h.edges[1]), maximum(report.h.edges[1])) color := :black - minimum(report.h.edges[1]):0.1:maximum(report.h.edges[1]), report.f_bck + minimum(report.h.edges[1]):ustrip(value(report.v.σ))/1000:maximum(report.h.edges[1]), report.f_bck end if !isempty(report.gof) ylims_res_max, ylims_res_min = 5, -5 @@ -224,7 +232,7 @@ end ylims := (ylim_min, ylim_max) xlims := (minimum(report.h.edges[1]), maximum(report.h.edges[1])) color := :black - minimum(report.h.edges[1]):0.1:maximum(report.h.edges[1]), report.f_bck + minimum(report.h.edges[1]):ustrip(value(report.v.σ))/1000:maximum(report.h.edges[1]), report.f_bck end end end @@ -259,19 +267,19 @@ end seriestype := :line label := "Best Fit" color := :red - minimum(report.h.edges[1]):1e-4:maximum(report.h.edges[1]), report.f_fit + minimum(report.h.edges[1]):ustrip(value(report.v.σ))/1000:maximum(report.h.edges[1]), report.f_fit end @series begin seriestype := :line label := "Signal" color := :green - minimum(report.h.edges[1]):1e-4:maximum(report.h.edges[1]), report.f_sig + minimum(report.h.edges[1]):ustrip(value(report.v.σ))/1000:maximum(report.h.edges[1]), report.f_sig end @series begin seriestype := :line label := "Background" color := :black - minimum(report.h.edges[1]):1e-4:maximum(report.h.edges[1]), report.f_bck + minimum(report.h.edges[1]):ustrip(value(report.v.σ))/1000:maximum(report.h.edges[1]), report.f_bck end end From 09d592d20372d962b2173cc73e540eea07732648 Mon Sep 17 00:00:00 2001 From: Florian Henkes Date: Thu, 4 Apr 2024 15:29:01 +0200 Subject: [PATCH 31/68] Include e_unit in report --- src/fwhm.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/fwhm.jl b/src/fwhm.jl index 8454ab60..8759ca2e 100644 --- a/src/fwhm.jl +++ b/src/fwhm.jl @@ -33,7 +33,7 @@ function fit_fwhm(peaks::Vector{<:Unitful.Energy{<:Real}}, fwhm::Vector{<:Unitfu # get fwhm at Qbb qbb = report_chi2.f_fit(2039) * e_unit result = merge(result_chi2, (par = par_unit , qbb = qbb, func = func, func_cal = func_cal, func_generic = func_generic, peaks = peaks, fwhm = fwhm)) - report = merge(report_chi2, (par = result.par, qbb = result.qbb, type = :fwhm)) + report = merge(report_chi2, (e_unit = e_unit, par = result.par, qbb = result.qbb, type = :fwhm)) return result, report end From f3bdf0743465293322b8f4b4a2c37d4e04c45906 Mon Sep 17 00:00:00 2001 From: Florian Henkes Date: Thu, 4 Apr 2024 15:29:27 +0200 Subject: [PATCH 32/68] Fixed recipes for A/E, included possibility to plot additional points in cal plots --- ext/LegendSpecFitsRecipesBaseExt.jl | 171 ++++++++++++++-------------- 1 file changed, 83 insertions(+), 88 deletions(-) diff --git a/ext/LegendSpecFitsRecipesBaseExt.jl b/ext/LegendSpecFitsRecipesBaseExt.jl index e2fdb7e0..9fba77f5 100644 --- a/ext/LegendSpecFitsRecipesBaseExt.jl +++ b/ext/LegendSpecFitsRecipesBaseExt.jl @@ -112,8 +112,8 @@ end end end -@recipe function f(report::NamedTuple{(:v, :h, :f_fit, :f_sig, :f_lowEtail, :f_bck, :gof)}; show_label=true, show_fit=true, _subplot=1) - # thickness_scaling := 2.0 +@recipe function f(report::NamedTuple{(:v, :h, :f_fit, :f_sig, :f_lowEtail, :f_bck, :gof)}; show_label=true, show_fit=true, f_fit_x_step_scaling=1/100, _subplot=1) + f_fit_x_step = ustrip(value(report.v.σ)) * f_fit_x_step_scaling legend := :topright ylim_max = max(3*value(report.f_sig(report.v.μ)), 3*maximum(report.h.weights)) ylim_max = ifelse(ylim_max == 0.0, 1e5, ylim_max) @@ -142,21 +142,21 @@ end color := :red subplot --> _subplot # ribbon := uncertainty.(report.f_fit.(minimum(report.h.edges[1]):0.1:maximum(report.h.edges[1]))) - minimum(report.h.edges[1]):ustrip(value(report.v.σ))/1000:maximum(report.h.edges[1]), value.(report.f_fit.(minimum(report.h.edges[1]):0.1:maximum(report.h.edges[1]))) + minimum(report.h.edges[1]):f_fit_x_step:maximum(report.h.edges[1]), value.(report.f_fit.(minimum(report.h.edges[1]):f_fit_x_step:maximum(report.h.edges[1]))) end @series begin seriestype := :line label := ifelse(show_label, "Signal", "") subplot --> _subplot color := :green - minimum(report.h.edges[1]):ustrip(value(report.v.σ))/1000:maximum(report.h.edges[1]), report.f_sig + minimum(report.h.edges[1]):f_fit_x_step:maximum(report.h.edges[1]), report.f_sig end @series begin seriestype := :line label := ifelse(show_label, "Low Energy Tail", "") subplot --> _subplot color := :blue - minimum(report.h.edges[1]):ustrip(value(report.v.σ))/1000:maximum(report.h.edges[1]), report.f_lowEtail + minimum(report.h.edges[1]):f_fit_x_step:maximum(report.h.edges[1]), report.f_lowEtail end @series begin seriestype := :line @@ -170,7 +170,7 @@ end ylims := (ylim_min, ylim_max) xlims := (minimum(report.h.edges[1]), maximum(report.h.edges[1])) color := :black - minimum(report.h.edges[1]):ustrip(value(report.v.σ))/1000:maximum(report.h.edges[1]), report.f_bck + minimum(report.h.edges[1]):f_fit_x_step:maximum(report.h.edges[1]), report.f_bck end if !isempty(report.gof) ylims_res_max, ylims_res_min = 5, -5 @@ -232,27 +232,14 @@ end ylims := (ylim_min, ylim_max) xlims := (minimum(report.h.edges[1]), maximum(report.h.edges[1])) color := :black - minimum(report.h.edges[1]):ustrip(value(report.v.σ))/1000:maximum(report.h.edges[1]), report.f_bck + minimum(report.h.edges[1]):f_fit_x_step:maximum(report.h.edges[1]), report.f_bck end end end end -# TODO: Add a recipe for the report_dict --> Feeeeeeeeeeeeelix :* -# @recipe function f(report_dict::Dict{Symbol, NamedTuple}) -# # layout := grid(length(report_dict)*2, 1, heights=vcat(fill([0.8, 0.2], length(report_dict))...)) -# size := (1800, 1000*length(report_dict)) -# bottom_margin := (0, :mm) -# for (i, k) in enumerate(string.(keys(report_dict))) -# @series begin -# title := string(k) -# _subplot --> 2*i-1 -# report_dict[Symbol(k)] -# end -# end -# end - -@recipe function f(report::NamedTuple{((:v, :h, :f_fit, :f_sig, :f_bck))}) +@recipe function f(report::NamedTuple{((:v, :h, :f_fit, :f_sig, :f_bck))}; f_fit_x_step_scaling=1/100) + f_fit_x_step = ustrip(value(report.v.σ)) * f_fit_x_step_scaling xlabel := "A/E (a.u.)" ylabel := "Counts" legend := :topleft @@ -267,19 +254,19 @@ end seriestype := :line label := "Best Fit" color := :red - minimum(report.h.edges[1]):ustrip(value(report.v.σ))/1000:maximum(report.h.edges[1]), report.f_fit + minimum(report.h.edges[1]):f_fit_x_step:maximum(report.h.edges[1]), report.f_fit end @series begin seriestype := :line label := "Signal" color := :green - minimum(report.h.edges[1]):ustrip(value(report.v.σ))/1000:maximum(report.h.edges[1]), report.f_sig + minimum(report.h.edges[1]):f_fit_x_step:maximum(report.h.edges[1]), report.f_sig end @series begin seriestype := :line label := "Background" color := :black - minimum(report.h.edges[1]):ustrip(value(report.v.σ))/1000:maximum(report.h.edges[1]), report.f_bck + minimum(report.h.edges[1]):f_fit_x_step:maximum(report.h.edges[1]), report.f_bck end end @@ -387,41 +374,6 @@ end subplot := 3 report_ctc.e_ctc, report_ctc.qdrift_peak end - # @series begin - # seriestype := :stepbins - # color := :red - # label := "Before CTC" - # xlabel := "Energy (keV)" - # ylabel := "Counts" - # title := "FWHM $(round(report_ctc.fwhm_before, digits=2))" - # yscale := :log10 - # subplot := 3 - # report_ctc.h_before - # end - # @series begin - # # seriestype := :stepbins - # color := :red - # # label := "Before CTC" - # # xlabel := "Energy (keV)" - # # ylabel := "Counts" - # # yscale := :log10 - # subplot := 3 - # # report_ctc.h_before - # minimum(report_ctc.e_peak):0.001:maximum(report_ctc.e_peak), t -> report_ctc.report_before.f_fit(t) - # end - # @series begin - # # seriestype := :stepbins - # color := :red - # # label := "Before CTC" - # # xlabel := "Energy (keV)" - # # ylabel := "Counts" - # # yscale := :log10 - # subplot := 4 - # # report_ctc.h_before - # minimum(report_ctc.e_peak):0.001:maximum(report_ctc.e_peak), t -> report_ctc.report_after.f_fit(t) - # end - - @series begin seriestype := :stepbins color := :red @@ -469,13 +421,6 @@ end linewidth := 3 ustrip.([report_window_cut.low_cut, report_window_cut.high_cut]) end - # @series begin - # seriestype := :vline - # label := "Center" - # color := :blue - # linewidth := 3 - # ustrip.([report_window_cut.center]) - # end @series begin seriestype := :vline label := "Fit Window" @@ -493,14 +438,14 @@ end end -@recipe function f(report::NamedTuple{(:par, :f_fit, :x, :y, :gof)}; plot_ribbon=true, xerrscaling=1) +@recipe function f(report::NamedTuple{(:par, :f_fit, :x, :y, :gof)}; plot_ribbon=true, xerrscaling=1, additional_pts=NamedTuple()) thickness_scaling := 2.0 xlims := (0, 1.2*value(maximum(report.x))) framestyle := :box xformatter := :plain yformatter := :plain layout --> @layout([a{0.8h}; b{0.2h}]) - margins --> (-15, :mm) + margins --> (-11.5, :mm) link --> :x @series begin seriestype := :line @@ -528,33 +473,63 @@ end label := "Data (Error x$(xerrscaling))" end markercolor --> :black - xerror := uncertainty.(report.x) .* xerrscaling value.(report.x), report.y end + if !isempty(additional_pts) + @series begin + seriestype := :scatter + subplot --> 1 + if xerrscaling == 1 + label := "Data not used for fit" + else + label := "Data not used for fit (Error x$(xerrscaling))" + end + ms --> 3 + markershape --> :circle + markerstrokecolor --> :black + linewidth --> 0.5 + markercolor --> :silver + xerror := uncertainty.(additional_pts.x) .* xerrscaling + value.(additional_pts.x), additional_pts.y + end + end @series begin seriestype := :hline - ribbon := 3 - subplot := 2 - fillalpha := 0.5 - label := "" - fillcolor := :lightgrey - linecolor := :darkgrey + ribbon --> 3 + subplot --> 2 + fillalpha --> 0.5 + label --> "" + fillcolor --> :lightgrey + linecolor --> :darkgrey [0.0] end @series begin - seriestype := :hline - ribbon := 1 - subplot := 2 - fillalpha := 0.5 - label := "" - fillcolor := :grey - linecolor := :darkgrey + seriestype --> :hline + ribbon --> 1 + subplot --> 2 + fillalpha --> 0.5 + label --> "" + fillcolor --> :grey + linecolor --> :darkgrey [0.0] end + if !isempty(additional_pts) + @series begin + seriestype := :scatter + label --> :none + ms --> 3 + markershape --> :circle + markerstrokecolor --> :black + linewidth --> 0.5 + markercolor --> :silver + subplot --> 2 + value.(additional_pts.x), additional_pts.residuals_norm + end + end @series begin - seriestype := :scatter - subplot := 2 - label := "" + seriestype --> :scatter + subplot --> 2 + label --> "" markercolor --> :black ylabel --> "Residuals (σ)" ylims --> (-5, 5) @@ -563,9 +538,18 @@ end end end -@recipe function f(report::NamedTuple{(:par, :f_fit, :x, :y, :gof, :qbb, :type)}) +@recipe function f(report::NamedTuple{(:par, :f_fit, :x, :y, :gof, :e_unit, :qbb, :type)}; additional_pts=NamedTuple()) bottom_margin --> (0, :mm) if report.type == :fwhm + y_max = value(maximum(report.y)) + additional_pts = if !isempty(additional_pts) + fwhm_cal = report.f_fit.(ustrip.(additional_pts.peaks)) .* report.e_unit + y_max = max(y_max, value(maximum(ustrip.(report.e_unit, additional_pts.fwhm)))) + (x = additional_pts.peaks, y = additional_pts.fwhm, + residuals_norm = (value.(fwhm_cal) .- additional_pts.fwhm)./ uncertainty.(fwhm_cal)) + else + NamedTuple() + end xlabel := "Energy (keV)" legend := :bottomright framestyle := :box @@ -573,6 +557,8 @@ end xticks := (0:500:3000, ["$i" for i in 0:500:3000]) @series begin grid --> :all + xerrscaling --> 1 + additional_pts --> additional_pts (par = report.par, f_fit = report.f_fit, x = report.x, y = report.y, gof = report.gof) end @series begin @@ -583,6 +569,7 @@ end linewidth := 2.5 xticks := :none ylabel := "FWHM" + ylims := (0, 1.2*y_max) subplot := 1 ribbon := uncertainty(report.qbb) [value(report.qbb)] @@ -590,9 +577,16 @@ end end end -@recipe function f(report::NamedTuple{(:par, :f_fit, :x, :y, :gof, :e_unit, :type)}; xerrscaling=1) +@recipe function f(report::NamedTuple{(:par, :f_fit, :x, :y, :gof, :e_unit, :type)}; xerrscaling=1, additional_pts=NamedTuple()) bottom_margin --> (0, :mm) if report.type == :cal + additional_pts = if !isempty(additional_pts) + μ_cal = report.f_fit.(additional_pts.μ) .* report.e_unit + (x = additional_pts.μ, y = additional_pts.peaks, + residuals_norm = (value.(μ_cal) .- additional_pts.peaks)./ uncertainty.(μ_cal)) + else + NamedTuple() + end xlabel := "Energy (ADC)" legend := :bottomright framestyle := :box @@ -601,6 +595,7 @@ end @series begin grid --> :all xerrscaling := xerrscaling + additional_pts := additional_pts (par = report.par, f_fit = report.f_fit, x = report.x, y = report.y, gof = report.gof) end @series begin From 9d0777e009dd7ee6a9746605d215725f8eabe7b7 Mon Sep 17 00:00:00 2001 From: Florian Henkes Date: Thu, 4 Apr 2024 15:29:37 +0200 Subject: [PATCH 33/68] Bug Fixes --- src/aoefit.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/aoefit.jl b/src/aoefit.jl index 41579639..61782d92 100644 --- a/src/aoefit.jl +++ b/src/aoefit.jl @@ -307,8 +307,8 @@ function fit_single_aoe_compton(h::Histogram, ps::NamedTuple; uncertainty::Bool= @debug "B: $(v_ml.B) ± $(v_ml_err.B)" result = merge(NamedTuple{keys(v_ml)}([measurement(v_ml[k], v_ml_err[k]) for k in keys(v_ml)]...), - (gof = (pvalue = pval, chi2 = chi2, dof = dof, covmat = param_covariance, - residuals = residuals, residuals_norm = residuals_norm, bin_centers = bin_centers),)) + (gof = (pvalue = pval, chi2 = chi2, dof = dof, covmat = param_covariance, + residuals = residuals, residuals_norm = residuals_norm, bin_centers = bin_centers),)) else @debug "Best Fit values" @debug "μ: $(v_ml.μ)" From 8394d78fbed7f1dfc45d040f60b2713e6d3938f0 Mon Sep 17 00:00:00 2001 From: Florian Henkes Date: Thu, 4 Apr 2024 19:34:08 +0200 Subject: [PATCH 34/68] Fixed QC fit reports --- ext/LegendSpecFitsRecipesBaseExt.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/ext/LegendSpecFitsRecipesBaseExt.jl b/ext/LegendSpecFitsRecipesBaseExt.jl index 9fba77f5..7f7908dc 100644 --- a/ext/LegendSpecFitsRecipesBaseExt.jl +++ b/ext/LegendSpecFitsRecipesBaseExt.jl @@ -403,7 +403,6 @@ end @recipe function f(report_window_cut::NamedTuple{(:h, :f_fit, :x_fit, :low_cut, :high_cut, :low_cut_fit, :high_cut_fit, :center, :σ)}) xlims := (value(ustrip(report_window_cut.center - 5*report_window_cut.σ)), value(ustrip(report_window_cut.center + 5*report_window_cut.σ))) - xlims := (ustrip(report_window_cut.center - 5*report_window_cut.σ), ustrip(report_window_cut.center + 5*report_window_cut.σ)) @series begin report_window_cut.h end From 676e728bd0356ca0663ba948fb574e8114b601ef Mon Sep 17 00:00:00 2001 From: Florian Henkes Date: Thu, 4 Apr 2024 19:34:31 +0200 Subject: [PATCH 35/68] Update f_aoe_sigma func for more stability --- src/aoe_calibration.jl | 21 ++++++++++++--------- 1 file changed, 12 insertions(+), 9 deletions(-) diff --git a/src/aoe_calibration.jl b/src/aoe_calibration.jl index c9389405..9d4a59c6 100644 --- a/src/aoe_calibration.jl +++ b/src/aoe_calibration.jl @@ -1,5 +1,6 @@ # f_aoe_sigma(x, p) = p[1] .+ p[2]*exp.(-p[3]./x) -@. f_aoe_sigma(x, p) = sqrt(abs(p[1]) + abs(p[2])/x^2) +# @. f_aoe_sigma(x, p) = sqrt(abs(p[1]) + abs(p[2])/x^2) +@. f_aoe_sigma(x, p) = sqrt(p[1]^2 + p[2]^2/x^2) f_aoe_mu(x, p) = p[1] .+ p[2].*x """ fit_aoe_corrections(e::Array{<:Unitful.Energy{<:Real}}, μ::Array{<:Real}, σ::Array{<:Real}) @@ -34,27 +35,29 @@ function fit_aoe_corrections(e::Array{<:Unitful.Energy{<:Real}}, μ::Array{<:Rea σ_cut = (mean(σ) - std(σ) .< σ .< mean(σ) + std(σ)) .&& muncert.(σ) .> 0.0 f_fit_σ = f_aoe_sigma # fit function result_σ, report_σ = chi2fit((x, p1, p2) -> f_fit_σ(x,[p1,p2]), e[σ_cut], σ[σ_cut]; uncertainty=true) - par_σ = [result_σ.par[1], result_σ.par[2] * e_unit^2] # add unit + par_σ = [result_σ.par[1], result_σ.par[2] * e_unit^2] func_σ = nothing func_generic_σ = nothing if string(f_fit_σ) == "f_aoe_sigma" - func_σ = "sqrt( abs($(mvalue(result_σ.par[1]))) + abs($(mvalue(result_σ.par[2]))$(e_unit)^2) / ($e_expression)^2)" # add units!! - func_generic_σ = "sqrt(abs(p[1]) + abs(p[2]) / ($e_expression)^2)" + # func_σ = "sqrt( abs($(mvalue(result_σ.par[1]))) + abs($(mvalue(result_σ.par[2]))$(e_unit)^2) / ($e_expression)^2)" + func_σ = "sqrt( ($(mvalue(result_σ.par[1])))^2 + ($(mvalue(result_σ.par[2]))$(e_unit))^2 / ($e_expression)^2 )" + # func_generic_σ = "sqrt(abs(p[1]) + abs(p[2]) / ($e_expression)^2)" + func_generic_σ = "sqrt( (p[1])^2 + (p[2])^2 / ($e_expression)^2 )" end - result_σ = merge(result_σ, (par = par_σ, func = func_σ, func_generic = func_generic_σ, σ = σ)) + result_σ = merge(result_σ, (par = par_σ, func = func_σ, func_generic = func_generic_σ, σ = σ)) report_σ = merge(report_σ, (e_unit = e_unit, label_y = "σ", label_fit = "Best fit: sqrt($(round(mvalue(result_σ.par[1])*1e6, digits=1))e-6 + $(round(ustrip(mvalue(result_σ.par[2])), digits=2)) / E^2)")) @debug "Compton band σ normalization: $(result_σ.func)" # put everything together into A/E correction/normalization function - aoe_str = "(a / (($e_expression)$e_unit^-1))" # get aoe, but without unit. - func_aoe_corr = "($aoe_str - ($(result_µ.func))) / ($(result_σ.func))" + aoe_str = "(a / ( ($e_expression) $e_unit^-1) )" # get aoe, but without unit. + func_aoe_corr = "($aoe_str - ($(result_µ.func)) ) / ($(result_σ.func))" func_generic_aoe_corr = "(aoe - $(result_µ.func_generic)) / $(result_σ.func_generic)" - func_aoe_corr_ecal = replace(func_aoe_corr,e_expression => "e_cal") # function that can be used for already calibrated energies + func_aoe_corr_ecal = replace(func_aoe_corr, e_expression => "e_cal") # function that can be used for already calibrated energies result = (µ_compton = result_µ, σ_compton = result_σ, compton_bands = (e = e,), func = func_aoe_corr, func_generic = func_generic_aoe_corr, func_ecal = func_aoe_corr_ecal) report = (report_µ = report_µ, report_σ = report_σ) - return result, report + return result, report end export fit_aoe_corrections From 6e0f206968e55c798f8e68eb6a9a5a01b965f73f Mon Sep 17 00:00:00 2001 From: Florian Henkes Date: Fri, 5 Apr 2024 11:12:00 +0200 Subject: [PATCH 36/68] Changed `psd` to `aoe` --- src/aoe_calibration.jl | 74 +++++++++++++++++++++--------------------- 1 file changed, 37 insertions(+), 37 deletions(-) diff --git a/src/aoe_calibration.jl b/src/aoe_calibration.jl index 9d4a59c6..4e2687f5 100644 --- a/src/aoe_calibration.jl +++ b/src/aoe_calibration.jl @@ -88,69 +88,69 @@ export prepare_dep_peakhist """ - get_n_after_psd_cut(psd_cut::T, aoe::Array{T}, e::Array{T}, peak::T, window::Array{T}, bin_width::T, result_before::NamedTuple, peakstats::NamedTuple; uncertainty=true) where T<:Real + get_n_after_aoe_cut(aoe_cut::T, aoe::Array{T}, e::Array{T}, peak::T, window::Array{T}, bin_width::T, result_before::NamedTuple, peakstats::NamedTuple; uncertainty=true) where T<:Real -Get the number of counts after a PSD cut value `psd_cut` for a given `peak` and `window` size whiile performing a peak fit with fixed position. The number of counts is determined by fitting the peak with a pseudo prior for the peak position. +Get the number of counts after a AoE cut value `aoe_cut` for a given `peak` and `window` size whiile performing a peak fit with fixed position. The number of counts is determined by fitting the peak with a pseudo prior for the peak position. # Returns - `n`: Number of counts after the cut - `n_err`: Uncertainty of the number of counts after the cut """ -function get_n_after_psd_cut(psd_cut::Unitful.RealOrRealQuantity, aoe::Vector{<:Unitful.RealOrRealQuantity}, e::Vector{<:T}, peak::T, window::Vector{T}, bin_width::T, result_before::NamedTuple, peakstats::NamedTuple; uncertainty::Bool=true, fixed_position::Bool=true) where T<:Unitful.Energy{<:Real} +function get_n_after_aoe_cut(aoe_cut::Unitful.RealOrRealQuantity, aoe::Vector{<:Unitful.RealOrRealQuantity}, e::Vector{<:T}, peak::T, window::Vector{T}, bin_width::T, result_before::NamedTuple, peakstats::NamedTuple; uncertainty::Bool=true, fixed_position::Bool=true) where T<:Unitful.Energy{<:Real} # get energy after cut and create histogram - peakhist = fit(Histogram, ustrip.(e[aoe .> psd_cut]), ustrip(peak-first(window):bin_width:peak+last(window))) + peakhist = fit(Histogram, ustrip.(e[aoe .> aoe_cut]), ustrip(peak-first(window):bin_width:peak+last(window))) # create pseudo_prior with known peak sigma in signal for more stable fit - pseudo_prior = NamedTupleDist(σ = Normal(result_before.σ, 0.3), ) + pseudo_prior = NamedTupleDist(σ = Normal(result_before.σ, 0.1), ) # fit peak and return number of signal counts result, _ = fit_single_peak_th228(peakhist, peakstats,; uncertainty=uncertainty, fixed_position=fixed_position, low_e_tail=false, pseudo_prior=pseudo_prior) return result.n end -export get_n_after_psd_cut +export get_n_after_aoe_cut """ - get_psd_cut(aoe::Vector{<:Unitful.RealOrRealQuantity}, e::Vector{<:T},; dep::T=1592.53u"keV", window::Vector{T}=[12.0, 10.0]u"keV", dep_sf::Float64=0.9, cut_search_interval::Tuple{Quantity{<:Real}, Quantity{<:Real}}=(-25.0u"keV^-1", 0.0u"keV^-1"), rtol::Float64=0.001, bin_width_window::T=3.0u"keV", fixed_position::Bool=true, sigma_high_sided::Float64=NaN, uncertainty::Bool=true, maxiters::Int=200) where T<:Unitful.Energy{<:Real} + get_aoe_cut(aoe::Vector{<:Unitful.RealOrRealQuantity}, e::Vector{<:T},; dep::T=1592.53u"keV", window::Vector{T}=[12.0, 10.0]u"keV", dep_sf::Float64=0.9, cut_search_interval::Tuple{Quantity{<:Real}, Quantity{<:Real}}=(-25.0u"keV^-1", 0.0u"keV^-1"), rtol::Float64=0.001, bin_width_window::T=3.0u"keV", fixed_position::Bool=true, sigma_high_sided::Float64=NaN, uncertainty::Bool=true, maxiters::Int=200) where T<:Unitful.Energy{<:Real} -Get the PSD cut value for a given `dep` and `window` size while performing a peak fit with fixed position. The PSD cut value is determined by finding the cut value for which the number of counts after the cut is equal to `dep_sf` times the number of counts before the cut. +Get the AoE cut value for a given `dep` and `window` size while performing a peak fit with fixed position. The AoE cut value is determined by finding the cut value for which the number of counts after the cut is equal to `dep_sf` times the number of counts before the cut. The algorhithm utilizes a root search algorithm to find the cut value with a relative tolerance of `rtol`. # Returns -- `cut`: PSD cut value +- `cut`: AoE cut value - `n0`: Number of counts before the cut - `nsf`: Number of counts after the cut """ -function get_psd_cut(aoe::Vector{<:Unitful.RealOrRealQuantity}, e::Vector{<:T},; dep::T=1592.53u"keV", window::Vector{<:T}=[12.0, 10.0]u"keV", dep_sf::Float64=0.9, cut_search_interval::Tuple{<:Unitful.RealOrRealQuantity, <:Unitful.RealOrRealQuantity}=(-25.0u"keV^-1", 0.0u"keV^-1"), rtol::Float64=0.001, bin_width_window::T=3.0u"keV", fixed_position::Bool=true, sigma_high_sided::Float64=NaN, uncertainty::Bool=true, maxiters::Int=200) where T<:Unitful.Energy{<:Real} - # check if a high sided AoE cut should be applied before the PSD cut is generated +function get_aoe_cut(aoe::Vector{<:Unitful.RealOrRealQuantity}, e::Vector{<:T},; dep::T=1592.53u"keV", window::Vector{<:T}=[12.0, 10.0]u"keV", dep_sf::Float64=0.9, cut_search_interval::Tuple{<:Unitful.RealOrRealQuantity, <:Unitful.RealOrRealQuantity}=(-25.0*unit(first(aoe)), 1.0*unit(first(aoe))), rtol::Float64=0.001, bin_width_window::T=3.0u"keV", fixed_position::Bool=true, sigma_high_sided::Float64=NaN, uncertainty::Bool=true, maxiters::Int=200) where T<:Unitful.Energy{<:Real} + # check if a high sided AoE cut should be applied before the AoE cut is generated if !isnan(sigma_high_sided) e = e[aoe .< sigma_high_sided] aoe = aoe[aoe .< sigma_high_sided] end + # cut window around peak + aoe = aoe[dep-first(window) .< e .< dep+last(window)] + e = e[dep-first(window) .< e .< dep+last(window)] # estimate bin width - bin_width = get_friedman_diaconis_bin_width(e[e .> dep - bin_width_window .&& e .< dep + bin_width_window]) + bin_width = get_friedman_diaconis_bin_width(e[dep - bin_width_window .< e .< dep + bin_width_window]) # create histogram dephist = fit(Histogram, ustrip.(e), ustrip(dep-first(window):bin_width:dep+last(window))) # get peakstats depstats = estimate_single_peak_stats(dephist) - # cut window around peak - aoe = aoe[dep-first(window) .< e .< dep+last(window)] - e = e[dep-first(window) .< e .< dep+last(window)] # fit before cut result_before, _ = fit_single_peak_th228(dephist, depstats,; uncertainty=uncertainty, fixed_position=fixed_position, low_e_tail=false) # get n0 before cut nsf = result_before.n * dep_sf - # get psd cut - n_surrival_dep_f = cut -> get_n_after_psd_cut(cut, aoe, e, dep, window, bin_width, mvalue(result_before), depstats; uncertainty=false, fixed_position=fixed_position) - nsf - psd_cut = find_zero(n_surrival_dep_f, cut_search_interval, Bisection(), rtol=rtol, maxiters=maxiters) + # get aoe cut + n_surrival_dep_f = cut -> get_n_after_aoe_cut(cut, aoe, e, dep, window, bin_width, mvalue(result_before), depstats; uncertainty=false, fixed_position=fixed_position) - nsf + aoe_cut = find_zero(n_surrival_dep_f, cut_search_interval, Bisection(), rtol=rtol, maxiters=maxiters) # return n_surrival_dep_f.(0.25:0.001:0.5) # get nsf after cut - nsf = get_n_after_psd_cut(psd_cut, aoe, e, dep, window, bin_width, mvalue(result_before), depstats; uncertainty=uncertainty, fixed_position=fixed_position) - return (lowcut = measurement(psd_cut, psd_cut * rtol), highcut = sigma_high_sided, n0 = result_before.n, nsf = nsf, sf = nsf / result_before.n * 100*u"percent") + nsf = get_n_after_aoe_cut(aoe_cut, aoe, e, dep, window, bin_width, mvalue(result_before), depstats; uncertainty=uncertainty, fixed_position=fixed_position) + return (lowcut = measurement(aoe_cut, aoe_cut * rtol), highcut = sigma_high_sided, n0 = result_before.n, nsf = nsf, sf = nsf / result_before.n * 100*u"percent") end -export get_psd_cut +export get_aoe_cut """ - get_peak_surrival_fraction(aoe::Array{T}, e::Array{T}, peak::T, window::Array{T}, psd_cut::T,; uncertainty::Bool=true, low_e_tail::Bool=true) where T<:Real + get_peak_surrival_fraction(aoe::Array{T}, e::Array{T}, peak::T, window::Array{T}, aoe_cut::T,; uncertainty::Bool=true, low_e_tail::Bool=true) where T<:Real -Get the surrival fraction of a peak after a PSD cut value `psd_cut` for a given `peak` and `window` size whiile performing a peak fit with fixed position. +Get the surrival fraction of a peak after a AoE cut value `aoe_cut` for a given `peak` and `window` size whiile performing a peak fit with fixed position. # Returns - `peak`: Peak position @@ -159,7 +159,7 @@ Get the surrival fraction of a peak after a PSD cut value `psd_cut` for a given - `sf`: Surrival fraction - `err`: Uncertainties """ -function get_peak_surrival_fraction(aoe::Vector{<:Unitful.RealOrRealQuantity}, e::Vector{<:T}, peak::T, window::Vector{T}, psd_cut::Unitful.RealOrRealQuantity,; uncertainty::Bool=true, low_e_tail::Bool=true, bin_width_window::T=2.0u"keV", sigma_high_sided::Unitful.RealOrRealQuantity=NaN) where T<:Unitful.Energy{<:Real} +function get_peak_surrival_fraction(aoe::Vector{<:Unitful.RealOrRealQuantity}, e::Vector{<:T}, peak::T, window::Vector{T}, aoe_cut::Unitful.RealOrRealQuantity,; uncertainty::Bool=true, low_e_tail::Bool=true, bin_width_window::T=2.0u"keV", sigma_high_sided::Unitful.RealOrRealQuantity=NaN) where T<:Unitful.Energy{<:Real} # estimate bin width bin_width = get_friedman_diaconis_bin_width(e[e .> peak - bin_width_window .&& e .< peak + bin_width_window]) # get energy before cut and create histogram @@ -171,9 +171,9 @@ function get_peak_surrival_fraction(aoe::Vector{<:Unitful.RealOrRealQuantity}, e # get e after cut if !isnan(sigma_high_sided) - e = e[psd_cut .< aoe .< sigma_high_sided] + e = e[aoe_cut .< aoe .< sigma_high_sided] else - e = e[psd_cut .< aoe] + e = e[aoe_cut .< aoe] end # estimate bin width @@ -210,21 +210,21 @@ export get_peak_surrival_fraction """ - get_peaks_surrival_fractions(aoe::Array{T}, e::Array{T}, peaks::Array{T}, peak_names::Array{Symbol}, windows::Array{T}, psd_cut::T,; uncertainty=true) where T<:Real + get_peaks_surrival_fractions(aoe::Array{T}, e::Array{T}, peaks::Array{T}, peak_names::Array{Symbol}, windows::Array{T}, aoe_cut::T,; uncertainty=true) where T<:Real -Get the surrival fraction of a peak after a PSD cut value `psd_cut` for a given `peak` and `window` size while performing a peak fit with fixed position. +Get the surrival fraction of a peak after a AoE cut value `aoe_cut` for a given `peak` and `window` size while performing a peak fit with fixed position. # Return - `result`: Dict of results for each peak - `report`: Dict of reports for each peak """ -function get_peaks_surrival_fractions(aoe::Vector{<:Unitful.RealOrRealQuantity}, e::Vector{<:T}, peaks::Vector{<:T}, peak_names::Vector{Symbol}, windows::Vector{<:Tuple{T, T}}, psd_cut::Unitful.RealOrRealQuantity,; uncertainty::Bool=true, bin_width_window::T=2.0u"keV", low_e_tail::Bool=true, sigma_high_sided::Unitful.RealOrRealQuantity=NaN) where T<:Unitful.Energy{<:Real} +function get_peaks_surrival_fractions(aoe::Vector{<:Unitful.RealOrRealQuantity}, e::Vector{<:T}, peaks::Vector{<:T}, peak_names::Vector{Symbol}, windows::Vector{<:Tuple{T, T}}, aoe_cut::Unitful.RealOrRealQuantity,; uncertainty::Bool=true, bin_width_window::T=2.0u"keV", low_e_tail::Bool=true, sigma_high_sided::Unitful.RealOrRealQuantity=NaN) where T<:Unitful.Energy{<:Real} # create return and result dicts result = Dict{Symbol, NamedTuple}() report = Dict{Symbol, NamedTuple}() # iterate throuh all peaks for (peak, name, window) in zip(peaks, peak_names, windows) # fit peak - result_peak, report_peak = get_peak_surrival_fraction(aoe, e, peak, collect(window), psd_cut; uncertainty=uncertainty, bin_width_window=bin_width_window, low_e_tail=low_e_tail, sigma_high_sided=sigma_high_sided) + result_peak, report_peak = get_peak_surrival_fraction(aoe, e, peak, collect(window), aoe_cut; uncertainty=uncertainty, bin_width_window=bin_width_window, low_e_tail=low_e_tail, sigma_high_sided=sigma_high_sided) # save results result[name] = result_peak report[name] = report_peak @@ -232,12 +232,12 @@ function get_peaks_surrival_fractions(aoe::Vector{<:Unitful.RealOrRealQuantity}, return result, report end export get_peaks_surrival_fractions -get_peaks_surrival_fractions(aoe, e, peaks, peak_names, left_window_sizes::Vector{<:Unitful.Energy{<:Real}}, right_window_sizes::Vector{<:Unitful.Energy{<:Real}}, psd_cut; kwargs...) = get_peaks_surrival_fractions(aoe, e, peaks, peak_names, [(l,r) for (l,r) in zip(left_window_sizes, right_window_sizes)], psd_cut; kwargs...) +get_peaks_surrival_fractions(aoe, e, peaks, peak_names, left_window_sizes::Vector{<:Unitful.Energy{<:Real}}, right_window_sizes::Vector{<:Unitful.Energy{<:Real}}, aoe_cut; kwargs...) = get_peaks_surrival_fractions(aoe, e, peaks, peak_names, [(l,r) for (l,r) in zip(left_window_sizes, right_window_sizes)], aoe_cut; kwargs...) """ - get_continuum_surrival_fraction(aoe:::Vector{<:Unitful.RealOrRealQuantity}, e::Vector{<:T}, center::T, window::T, psd_cut::Unitful.RealOrRealQuantity,; sigma_high_sided::Unitful.RealOrRealQuantity=NaN) where T<:Unitful.Energy{<:Real} + get_continuum_surrival_fraction(aoe:::Vector{<:Unitful.RealOrRealQuantity}, e::Vector{<:T}, center::T, window::T, aoe_cut::Unitful.RealOrRealQuantity,; sigma_high_sided::Unitful.RealOrRealQuantity=NaN) where T<:Unitful.Energy{<:Real} -Get the surrival fraction of a continuum after a PSD cut value `psd_cut` for a given `center` and `window` size. +Get the surrival fraction of a continuum after a AoE cut value `aoe_cut` for a given `center` and `window` size. # Returns - `center`: Center of the continuum @@ -246,15 +246,15 @@ Get the surrival fraction of a continuum after a PSD cut value `psd_cut` for a g - `n_after`: Number of counts after the cut - `sf`: Surrival fraction """ -function get_continuum_surrival_fraction(aoe::Vector{<:Unitful.RealOrRealQuantity}, e::Vector{<:T}, center::T, window::T, psd_cut::Unitful.RealOrRealQuantity,; sigma_high_sided::Unitful.RealOrRealQuantity=NaN) where T<:Unitful.Energy{<:Real} +function get_continuum_surrival_fraction(aoe::Vector{<:Unitful.RealOrRealQuantity}, e::Vector{<:T}, center::T, window::T, aoe_cut::Unitful.RealOrRealQuantity,; sigma_high_sided::Unitful.RealOrRealQuantity=NaN) where T<:Unitful.Energy{<:Real} # get number of events in window before cut n_before = length(e[center - window .< e .< center + window]) # get number of events after cut - n_after = length(e[aoe .> psd_cut .&& center - window .< e .< center + window]) + n_after = length(e[aoe .> aoe_cut .&& center - window .< e .< center + window]) if !isnan(sigma_high_sided) - n_after = length(e[psd_cut .< aoe .< sigma_high_sided .&& center - window .< e .< center + window]) + n_after = length(e[aoe_cut .< aoe .< sigma_high_sided .&& center - window .< e .< center + window]) end - n_after = length(e[aoe .> psd_cut .&& center - window .< e .< center + window]) + n_after = length(e[aoe .> aoe_cut .&& center - window .< e .< center + window]) # calculate surrival fraction sf = n_after / n_before result = ( From b0727e22ac526bd63fde4a7c1aa3dec3047fb3bf Mon Sep 17 00:00:00 2001 From: Florian Henkes Date: Fri, 5 Apr 2024 11:36:43 +0200 Subject: [PATCH 37/68] Bug Fix expression string in AoE calibrations --- src/aoe_calibration.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/aoe_calibration.jl b/src/aoe_calibration.jl index 4e2687f5..aa01dfd6 100644 --- a/src/aoe_calibration.jl +++ b/src/aoe_calibration.jl @@ -49,7 +49,7 @@ function fit_aoe_corrections(e::Array{<:Unitful.Energy{<:Real}}, μ::Array{<:Rea @debug "Compton band σ normalization: $(result_σ.func)" # put everything together into A/E correction/normalization function - aoe_str = "(a / ( ($e_expression) $e_unit^-1) )" # get aoe, but without unit. + aoe_str = "(a / ( ($e_expression)$e_unit^-1) )" # get aoe, but without unit. func_aoe_corr = "($aoe_str - ($(result_µ.func)) ) / ($(result_σ.func))" func_generic_aoe_corr = "(aoe - $(result_µ.func_generic)) / $(result_σ.func_generic)" func_aoe_corr_ecal = replace(func_aoe_corr, e_expression => "e_cal") # function that can be used for already calibrated energies From 9832f84e49ff86147da3b5c5a51be6bf42820d49 Mon Sep 17 00:00:00 2001 From: Florian Henkes Date: Fri, 12 Apr 2024 14:14:32 +0200 Subject: [PATCH 38/68] Bug Fix func_ecal generation --- src/aoe_calibration.jl | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/src/aoe_calibration.jl b/src/aoe_calibration.jl index aa01dfd6..e4f2b3de 100644 --- a/src/aoe_calibration.jl +++ b/src/aoe_calibration.jl @@ -52,7 +52,7 @@ function fit_aoe_corrections(e::Array{<:Unitful.Energy{<:Real}}, μ::Array{<:Rea aoe_str = "(a / ( ($e_expression)$e_unit^-1) )" # get aoe, but without unit. func_aoe_corr = "($aoe_str - ($(result_µ.func)) ) / ($(result_σ.func))" func_generic_aoe_corr = "(aoe - $(result_µ.func_generic)) / $(result_σ.func_generic)" - func_aoe_corr_ecal = replace(func_aoe_corr, e_expression => "e_cal") # function that can be used for already calibrated energies + func_aoe_corr_ecal = replace(func_aoe_corr, string(e_expression) => "e_cal") # function that can be used for already calibrated energies result = (µ_compton = result_µ, σ_compton = result_σ, compton_bands = (e = e,), func = func_aoe_corr, func_generic = func_generic_aoe_corr, func_ecal = func_aoe_corr_ecal) report = (report_µ = report_µ, report_σ = report_σ) @@ -99,9 +99,13 @@ function get_n_after_aoe_cut(aoe_cut::Unitful.RealOrRealQuantity, aoe::Vector{<: # get energy after cut and create histogram peakhist = fit(Histogram, ustrip.(e[aoe .> aoe_cut]), ustrip(peak-first(window):bin_width:peak+last(window))) # create pseudo_prior with known peak sigma in signal for more stable fit - pseudo_prior = NamedTupleDist(σ = Normal(result_before.σ, 0.1), ) + pseudo_prior = if fixed_position + NamedTupleDist(σ = Normal(result_before.σ, 0.1), μ = ConstValueDist(result_before.μ)) + else + NamedTupleDist(σ = Normal(result_before.σ, 0.1), ) + end # fit peak and return number of signal counts - result, _ = fit_single_peak_th228(peakhist, peakstats,; uncertainty=uncertainty, fixed_position=fixed_position, low_e_tail=false, pseudo_prior=pseudo_prior) + result, _ = fit_single_peak_th228(peakhist, peakstats,; uncertainty=uncertainty, low_e_tail=false, pseudo_prior=pseudo_prior) return result.n end export get_n_after_aoe_cut @@ -181,7 +185,7 @@ function get_peak_surrival_fraction(aoe::Vector{<:Unitful.RealOrRealQuantity}, e # get energy after cut and create histogram peakhist = fit(Histogram, ustrip(e), ustrip(peak-first(window):bin_width:peak+last(window))) # create pseudo_prior with known peak sigma in signal for more stable fit - pseudo_prior = NamedTupleDist(μ = ConstValueDist(result_before.μ), σ = Normal(result_before.σ, 0.1)) + # pseudo_prior = NamedTupleDist(μ = ConstValueDist(result_before.μ), σ = Normal(result_before.σ, 0.1)) pseudo_prior = NamedTupleDist(σ = Normal(result_before.σ, 0.1), ) # estimate peak stats peakstats = estimate_single_peak_stats(peakhist) From 65fcf32f93728bdd005e1c21114d4e165c95d500 Mon Sep 17 00:00:00 2001 From: LisaSchlueter Date: Sat, 13 Apr 2024 05:59:49 +0200 Subject: [PATCH 39/68] add plot receipe for covariance matrix --- ext/LegendSpecFitsRecipesBaseExt.jl | 55 +++++++++++++++++++++++++++++ src/specfit.jl | 3 +- 2 files changed, 57 insertions(+), 1 deletion(-) diff --git a/ext/LegendSpecFitsRecipesBaseExt.jl b/ext/LegendSpecFitsRecipesBaseExt.jl index 7f7908dc..07309eef 100644 --- a/ext/LegendSpecFitsRecipesBaseExt.jl +++ b/ext/LegendSpecFitsRecipesBaseExt.jl @@ -238,6 +238,61 @@ end end end +@recipe function f(report::NamedTuple{(:v, :h, :f_fit, :f_sig, :f_lowEtail, :f_bck, :gof)}, cormat::Bool = true) + cm = cor(report.gof.covmat) + cm_plt = NaN .* cm + for i in range(1, stop = size(cm)[1]) + cm_plt[i:end,i] = cm[i:end,i] + end + + # prepare labels + par_names = fieldnames(report.v) + tick_names = String.(collect(par_names)) + tick_names[tick_names .== "step_amplitude"] .= L"\textrm{bkg}_\textrm{step}" + tick_names[tick_names .== "background"] .= "bkg" + tick_names[tick_names .== "skew_width"] .= L"\textrm{tail}_\textrm{skew}" + tick_names[tick_names .== "skew_fraction"] .= L"\textrm{tail}_\textrm{frac}" + + @series begin + seriestype := :heatmap + # c := cgrad(:grays, rev = true) + c := :RdBu_3#:bam + cm_plt + end + + # annotation for correlation coefficient + cm_vec = round.(vec(cm_plt), digits = 2) + xvec = vec(hcat([fill(i, 7) for i in 1:7]...)) + yvec = vec(hcat([fill(i, 7) for i in 1:7]...)') + xvec = xvec[isfinite.(cm_vec)] + yvec = yvec[isfinite.(cm_vec)] + cm_vec = cm_vec[isfinite.(cm_vec)] + @series begin + seriestype := :scatter + xticks := (1:length(par_names),tick_names) + yticks := (1:length(par_names),tick_names) + yflip := true + markeralpha := 0 + colorbar := false + label := false + legend := false + thickness_scaling := 1.0 + xtickfontsize := 12 + xlabelfontsize := 14 + ylabelfontsize := 14 + tickfontsize:= 12 + legendfontsize := 10 + xlims := (0.5, 7.5) + ylims := (0.5, 7.5) + size := (600, 575) + grid := false + title := "Correlation Matrix" + series_annotations := [("$(cm_vec[i])", :center, :center, :black, 12, "Helvetica Bold") for i in eachindex(xvec)] + xvec, yvec + end + +end + @recipe function f(report::NamedTuple{((:v, :h, :f_fit, :f_sig, :f_bck))}; f_fit_x_step_scaling=1/100) f_fit_x_step = ustrip(value(report.v.σ)) * f_fit_x_step_scaling xlabel := "A/E (a.u.)" diff --git a/src/specfit.jl b/src/specfit.jl index ca5b110c..8ef7735e 100644 --- a/src/specfit.jl +++ b/src/specfit.jl @@ -168,9 +168,10 @@ function fit_single_peak_th228(h::Histogram, ps::NamedTuple{(:peak_pos, :peak_fw n = weibull_from_mx(ps.peak_counts, 2*ps.peak_counts), step_amplitude = weibull_from_mx(ps.mean_background_step, ps.mean_background_step + 5*ps.mean_background_std), skew_fraction = ifelse(low_e_tail, truncated(weibull_from_mx(0.01, 0.05), 0.0, 0.25), ConstValueDist(0.0)), - skew_width = ifelse(low_e_tail, weibull_from_mx(0.001, 1e-3), ConstValueDist(1.0)), + skew_width = ifelse(low_e_tail, weibull_from_mx(0.001, 1e-2), ConstValueDist(1.0)), background = weibull_from_mx(ps.mean_background, ps.mean_background + 5*ps.mean_background_std), ) + # use standard priors in case of no overwrites given if !(:empty in keys(pseudo_prior)) # check if input overwrite prior has the same fields as the standard prior set From 97f88bec59f1fead1c300f6080091677183af317 Mon Sep 17 00:00:00 2001 From: LisaSchlueter Date: Mon, 15 Apr 2024 19:35:44 +0200 Subject: [PATCH 40/68] update cormat and covmat receipe --- ext/LegendSpecFitsRecipesBaseExt.jl | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/ext/LegendSpecFitsRecipesBaseExt.jl b/ext/LegendSpecFitsRecipesBaseExt.jl index 07309eef..42c8380c 100644 --- a/ext/LegendSpecFitsRecipesBaseExt.jl +++ b/ext/LegendSpecFitsRecipesBaseExt.jl @@ -238,8 +238,16 @@ end end end -@recipe function f(report::NamedTuple{(:v, :h, :f_fit, :f_sig, :f_lowEtail, :f_bck, :gof)}, cormat::Bool = true) - cm = cor(report.gof.covmat) +@recipe function f(report::NamedTuple{(:v, :h, :f_fit, :f_sig, :f_lowEtail, :f_bck, :gof)}, mode::Symbol = :cormat) + if mode == :cormat + cm = cor(report.gof.covmat) + elseif mode == :covmat + cm = report.gof.covmat + else + @debug "mode $mode not supported" + return + end + cm_plt = NaN .* cm for i in range(1, stop = size(cm)[1]) cm_plt[i:end,i] = cm[i:end,i] From 6027b1d06eab97326607be21e9ee4f317287893e Mon Sep 17 00:00:00 2001 From: LisaSchlueter Date: Mon, 15 Apr 2024 19:53:11 +0200 Subject: [PATCH 41/68] bug fix plot receipe --- ext/LegendSpecFitsRecipesBaseExt.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/ext/LegendSpecFitsRecipesBaseExt.jl b/ext/LegendSpecFitsRecipesBaseExt.jl index 42c8380c..2ba550c3 100644 --- a/ext/LegendSpecFitsRecipesBaseExt.jl +++ b/ext/LegendSpecFitsRecipesBaseExt.jl @@ -238,13 +238,13 @@ end end end -@recipe function f(report::NamedTuple{(:v, :h, :f_fit, :f_sig, :f_lowEtail, :f_bck, :gof)}, mode::Symbol = :cormat) +@recipe function f(report::NamedTuple{(:v, :h, :f_fit, :f_sig, :f_lowEtail, :f_bck, :gof)}, mode::Symbol) if mode == :cormat cm = cor(report.gof.covmat) elseif mode == :covmat cm = report.gof.covmat else - @debug "mode $mode not supported" + @debug "mode $mode not supported - has to be :cormat or :covmat" return end From e3f9df2a237a26feba2e9a1052d68e0d603300ec Mon Sep 17 00:00:00 2001 From: LisaSchlueter Date: Mon, 15 Apr 2024 20:30:39 +0200 Subject: [PATCH 42/68] implement Felix suggestion to make arguments optional --- ext/LegendSpecFitsRecipesBaseExt.jl | 28 ++++++++++++++-------------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/ext/LegendSpecFitsRecipesBaseExt.jl b/ext/LegendSpecFitsRecipesBaseExt.jl index 2ba550c3..37d13f55 100644 --- a/ext/LegendSpecFitsRecipesBaseExt.jl +++ b/ext/LegendSpecFitsRecipesBaseExt.jl @@ -279,22 +279,22 @@ end seriestype := :scatter xticks := (1:length(par_names),tick_names) yticks := (1:length(par_names),tick_names) - yflip := true + yflip --> true markeralpha := 0 - colorbar := false - label := false - legend := false + colorbar --> false + label --> false + legend --> false thickness_scaling := 1.0 - xtickfontsize := 12 - xlabelfontsize := 14 - ylabelfontsize := 14 - tickfontsize:= 12 - legendfontsize := 10 - xlims := (0.5, 7.5) - ylims := (0.5, 7.5) - size := (600, 575) - grid := false - title := "Correlation Matrix" + xtickfontsize --> 12 + xlabelfontsize --> 14 + ylabelfontsize --> 14 + tickfontsize --> 12 + legendfontsize --> 10 + xlims --> (0.5, length(par_names)+0.5) + ylims --> (0.5, length(par_names)+0.5) + size --> (600, 575) + grid --> false + title --> "Correlation Matrix" series_annotations := [("$(cm_vec[i])", :center, :center, :black, 12, "Helvetica Bold") for i in eachindex(xvec)] xvec, yvec end From ea3c963339437158e59e1d44e74d8e10423951d7 Mon Sep 17 00:00:00 2001 From: Florian Henkes Date: Sat, 20 Apr 2024 19:06:55 +0200 Subject: [PATCH 43/68] Fix problem with decay time fitting --- ext/LegendSpecFitsRecipesBaseExt.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/ext/LegendSpecFitsRecipesBaseExt.jl b/ext/LegendSpecFitsRecipesBaseExt.jl index 37d13f55..4845547b 100644 --- a/ext/LegendSpecFitsRecipesBaseExt.jl +++ b/ext/LegendSpecFitsRecipesBaseExt.jl @@ -3,11 +3,11 @@ module LegendSpecFitsRecipesBaseExt using RecipesBase -using Unitful, Formatting, Measurements, LaTeXStrings +using Unitful, Format, Measurements, LaTeXStrings using Measurements: value, uncertainty using StatsBase, LinearAlgebra -function round_wo_units(x::Unitful.RealOrRealQuantity, digits::Integer=2) +function round_wo_units(x::Unitful.RealOrRealQuantity; digits::Integer=2) if unit(x) == NoUnits round(x, digits=digits) else @@ -19,12 +19,13 @@ end @recipe function f(report::NamedTuple{(:f_fit, :μ, :σ, :n)}, x::Vector{T}, cuts::NamedTuple{(:low, :high, :max), Tuple{T, T, T}}) where {T <: Unitful.RealOrRealQuantity} ylabel := "Normalized Counts" legend := :bottomright + framestyle := :box @series begin seriestype := :histogram bins --> :fd normalize --> :pdf label := "Data" - ustrip(x[x .> cuts.low .&& x .< cuts.high]) + ustrip.(x[x .> cuts.low .&& x .< cuts.high]) end @series begin color := :red From b5313095ec8aa93d9c049953721d004d0094e342 Mon Sep 17 00:00:00 2001 From: Florian Henkes Date: Sat, 20 Apr 2024 19:07:10 +0200 Subject: [PATCH 44/68] Switch from Formatting to Format.jl --- Project.toml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index 83250c76..a7111766 100644 --- a/Project.toml +++ b/Project.toml @@ -10,7 +10,7 @@ ChangesOfVariables = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0" DensityInterface = "b429d917-457f-4dbc-8f4c-0cc954292b1d" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b" -Formatting = "59287772-0a20-5a39-b81b-1366585eb4c0" +Format = "1fa38f19-a742-5d3f-a2b9-30dd87b9d5f8" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953" InverseFunctions = "3587e190-3f89-42d0-90ee-14403ec27112" @@ -52,7 +52,7 @@ ChangesOfVariables = "0.1.1" DensityInterface = "0.4" Distributions = "0.24, 0.25" FillArrays = "0.7,0.8, 0.9, 0.10, 0.11, 0.12, 0.13, 1" -Formatting = "0.4" +Format = "1.2, 1.3" ForwardDiff = "0.10" IntervalSets = "0.7" InverseFunctions = "0.1" From 32f57a759734f0203d2cb367518fdab6a890c81f Mon Sep 17 00:00:00 2001 From: Felix Hagemann Date: Thu, 25 Apr 2024 16:12:34 +0200 Subject: [PATCH 45/68] Add chisq fit tests --- test/runtests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index 822649c6..839ceaad 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -5,7 +5,7 @@ import Test Test.@testset "Package LegendSpecFits" begin #include("test_aqua.jl") include("test_specfit.jl") - # include("test_fit_chisq.jl") + include("test_fit_chisq.jl") include("test_docs.jl") isempty(Test.detect_ambiguities(LegendSpecFits)) end # testset From c753d8291e3230163cfac19c9bba5e7d23d8148d Mon Sep 17 00:00:00 2001 From: Felix Hagemann Date: Sun, 28 Apr 2024 21:10:45 +0200 Subject: [PATCH 46/68] Fix ylims of residuals in plot recipe --- ext/LegendSpecFitsRecipesBaseExt.jl | 31 +++++++++-------------------- 1 file changed, 9 insertions(+), 22 deletions(-) diff --git a/ext/LegendSpecFitsRecipesBaseExt.jl b/ext/LegendSpecFitsRecipesBaseExt.jl index 4845547b..9c66d582 100644 --- a/ext/LegendSpecFitsRecipesBaseExt.jl +++ b/ext/LegendSpecFitsRecipesBaseExt.jl @@ -125,7 +125,7 @@ end label := ifelse(show_label, "Data", "") bins --> :sqrt yscale := :log10 - ylims := (ylim_min, ylim_max) + ylims --> (ylim_min, ylim_max) xlabel := "Energy (keV)" xlims := (minimum(report.h.edges[1]), maximum(report.h.edges[1])) ylabel := "Counts / $(round(step(report.h.edges[1]), digits=2)) keV" @@ -165,10 +165,12 @@ end subplot --> _subplot margins --> (0, :mm) bottom_margin --> (-4, :mm) - xlabel := "" - xticks --> ([]) + if !isempty(report.gof) + xlabel := "" + xticks --> ([]) + end ylabel := "Counts / $(round(step(report.h.edges[1]), digits=2)) keV" - ylims := (ylim_min, ylim_max) + ylims --> (ylim_min, ylim_max) xlims := (minimum(report.h.edges[1]), maximum(report.h.edges[1])) color := :black minimum(report.h.edges[1]):f_fit_x_step:maximum(report.h.edges[1]), report.f_bck @@ -208,33 +210,18 @@ end label := "" title := "" markercolor --> :black - ylabel --> "Residuals (σ)" + ylabel := "Residuals (σ)" xlabel := "Energy (keV)" link --> :x top_margin --> (0, :mm) - ylims --> (ylims_res_min, ylims_res_max) + ylims := (ylims_res_min, ylims_res_max) xlims := (minimum(report.h.edges[1]), maximum(report.h.edges[1])) yscale --> :identity if ylims_res_max == 5 - yticks --> ([-3, 0, 3]) + yticks := ([-3, 0, 3]) end report.gof.bin_centers, report.gof.residuals_norm end - else - @series begin - seriestype := :line - label := ifelse(show_label, "Background", "") - subplot --> _subplot - margins --> (0, :mm) - bottom_margin --> (-4, :mm) - xlabel := "Energy (keV)" - xticks --> ([]) - ylabel := "Counts / $(round(step(report.h.edges[1]), digits=2)) keV" - ylims := (ylim_min, ylim_max) - xlims := (minimum(report.h.edges[1]), maximum(report.h.edges[1])) - color := :black - minimum(report.h.edges[1]):f_fit_x_step:maximum(report.h.edges[1]), report.f_bck - end end end end From 7bd52c0c2fd8bac806e5c95d5bfe5add32c02398 Mon Sep 17 00:00:00 2001 From: Felix Hagemann Date: Sun, 28 Apr 2024 23:17:50 +0200 Subject: [PATCH 47/68] Add function to estimate survival fraction from subpeaks --- src/specfit.jl | 197 ++++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 196 insertions(+), 1 deletion(-) diff --git a/src/specfit.jl b/src/specfit.jl index 8ef7735e..176786a3 100644 --- a/src/specfit.jl +++ b/src/specfit.jl @@ -322,4 +322,199 @@ function get_peak_fwhm_th228(v_ml::NamedTuple, v_ml_err::Union{Matrix,NamedTuple fwhm_err = std(fwhm_mc[isfinite.(fwhm_mc)]) return fwhm, fwhm_err end -export get_peak_fwhm_th228 \ No newline at end of file +export get_peak_fwhm_th228 + + + +""" + fit_single_peak_th228(h::Histogram, ps::NamedTuple{(:peak_pos, :peak_fwhm, :peak_sigma, :peak_counts, :mean_background), NTuple{5, T}};, uncertainty::Bool=true, fixed_position::Bool=false, low_e_tail::Bool=true) where T<:Real + +Perform a simultaneous fit of two peaks (`h_survived` and `h_cut`) that together would form a histogram `h`, from which the result `h_result` was already determined using `fit_single_peak_th228`. +Also, FWHM is calculated from the fitted peakshape with MC error propagation. The peak position can be fixed to the value in `ps` by setting `fixed_position=true`. If the low-E tail should not be fitted, it can be disabled by setting `low_e_tail=false`. +# Returns + * `result`: NamedTuple of the fit results containing values and errors, in particular the signal survival fraction `sf` and the background survival frachtion `bsf`. + * `report`: NamedTuple of the fit report which can be plotted +""" +function fit_subpeaks_th228( + h_survived::Histogram, h_cut::Histogram, h_result; + uncertainty::Bool=false, low_e_tail::Bool=true, pseudo_prior::NamedTupleDist=NamedTupleDist(empty = true), fit_func::Symbol=:f_fit +) + + # create standard pseudo priors + standard_pseudo_prior = let ps = h_result, ps_cut = estimate_single_peak_stats(h_cut), ps_survived = estimate_single_peak_stats(h_survived) + NamedTupleDist( + μ = ConstValueDist(Measurements.value(ps.μ)), + σ = ConstValueDist(Measurements.value(ps.σ)), + n = ConstValueDist(Measurements.value(ps.n)), + sf = Uniform(0,1), + bsf = Uniform(0,1), + step_amplitude_survived = weibull_from_mx(ps_survived.mean_background_step, ps_survived.mean_background_step + 5*ps_survived.mean_background_std), + skew_fraction_survived = ifelse(low_e_tail, truncated(weibull_from_mx(0.01, 0.05), 0.0, 0.1), ConstValueDist(0.0)), + skew_width_survived = ifelse(low_e_tail, weibull_from_mx(0.001, 1e-2), ConstValueDist(1.0)), + step_amplitude_cut = weibull_from_mx(ps_cut.mean_background_step, ps_cut.mean_background_step + 5*ps_cut.mean_background_std), + skew_fraction_cut = ifelse(low_e_tail, truncated(weibull_from_mx(0.01, 0.05), 0.0, 0.1), ConstValueDist(0.0)), + skew_width_cut = ifelse(low_e_tail, weibull_from_mx(0.001, 1e-2), ConstValueDist(1.0)), + background = ConstValueDist(Measurements.value(ps.background)) + ) + end + + # use standard priors in case of no overwrites given + if !(:empty in keys(pseudo_prior)) + # check if input overwrite prior has the same fields as the standard prior set + @assert all(f -> f in keys(standard_pseudo_prior), keys(standard_pseudo_prior)) "Pseudo priors can only have $(keys(standard_pseudo_prior)) as fields." + # replace standard priors with overwrites + pseudo_prior = merge(standard_pseudo_prior, pseudo_prior) + else + # take standard priors as pseudo priors with overwrites + pseudo_prior = standard_pseudo_prior + end + + # transform back to frequency space + f_trafo = BAT.DistributionTransform(Normal, pseudo_prior) + + # start values for MLE + v_init = Vector(mean(f_trafo.target_dist)) + + # create loglikehood function: f_loglike(v) that can be evaluated for any set of v (fit parameter) + f_loglike = let f_fit=th228_fit_functions[fit_func], h_cut=h_cut, h_survived=h_survived + v -> begin + v_survived = (μ = v.μ, σ = v.σ, n = v.n * v.sf, + step_amplitude = v.step_amplitude_survived, + skew_fraction = v.skew_fraction_survived, + skew_width = v.skew_width_survived, + background = v.background * v.bsf + ) + v_cut = (μ = v.μ, σ = v.σ, n = v.n * (1 - v.sf), + step_amplitude = v.step_amplitude_cut, + skew_fraction = v.skew_fraction_cut, + skew_width = v.skew_width_cut, + background = v.background * (1 - v.bsf) + ) + hist_loglike(Base.Fix2(f_fit, v_survived), h_survived) + hist_loglike(Base.Fix2(f_fit, v_cut), h_cut) + end + end + + # MLE + opt_r = optimize((-) ∘ f_loglike ∘ inverse(f_trafo), v_init, Optim.Options(time_limit = 60, iterations = 500)) + + # best fit results + v_ml = inverse(f_trafo)(Optim.minimizer(opt_r)) + + v_ml_survived = ( + μ = v_ml.μ, + σ = v_ml.σ, + n = v_ml.n * v_ml.sf, + step_amplitude = v_ml.step_amplitude_survived, + skew_fraction = v_ml.skew_fraction_survived, + skew_width = v_ml.skew_width_survived, + background = v_ml.background * v_ml.bsf + ) + + v_ml_cut = ( + μ = v_ml.μ, + σ = v_ml.σ, + n = v_ml.n * (1 - v_ml.sf), + step_amplitude = v_ml.step_amplitude_cut, + skew_fraction = v_ml.skew_fraction_cut, + skew_width = v_ml.skew_width_cut, + background = v_ml.background * (1 - v_ml.bsf) + ) + + gof_survived = NamedTuple() + gof_cut = NamedTuple() + + if uncertainty + + f_loglike_array = let v_keys = keys(pseudo_prior) + v -> -f_loglike(NamedTuple{v_keys}(v)) + end + + # Calculate the Hessian matrix using ForwardDiff + H = ForwardDiff.hessian(f_loglike_array, tuple_to_array(v_ml)) + + # Calculate the parameter covariance matrix + param_covariance_raw = inv(H) + param_covariance = nearestSPD(param_covariance_raw) + + # Extract the parameter uncertainties + v_ml_err = array_to_tuple(sqrt.(abs.(diag(param_covariance))), v_ml) + + # calculate all of this for each histogram individually + gofs = [ + begin + + h_part = Dict("survived" => h_survived, "cut" => h_cut)[part] + v_ml_part = Dict("survived" => v_ml_survived, "cut" => v_ml_cut)[part] + + # calculate p-value + pval, chi2, dof = p_value(th228_fit_functions[fit_func], h_part, v_ml_part) + + # calculate normalized residuals + residuals, residuals_norm, p_value_binwise, bin_centers = get_residuals(th228_fit_functions[fit_func], h_part, v_ml_part) + + gof = ( + pvalue = pval, chi2 = chi2, dof = dof, + covmat = param_covariance, covmat_raw = param_covariance_raw, + residuals = residuals, residuals_norm = residuals_norm, + pvalue_binwise = p_value_binwise, bin_centers = bin_centers + ) + + end for part in ("survived", "cut") + ] + + # get fwhm of peak + fwhm, fwhm_err = try + get_peak_fwhm_th228(v_ml, param_covariance) + catch e + get_peak_fwhm_th228(v_ml, v_ml_err) + end + + @debug "Best Fit values for $(part)" + @debug "μ: $(v_ml.μ) ± $(v_ml_err.μ)" + @debug "σ: $(v_ml.σ) ± $(v_ml_err.σ)" + @debug "n: $(v_ml.n) ± $(v_ml_err.n)" + @debug "p: $pval , chi2 = $(chi2) with $(dof) dof" + @debug "FWHM: $(fwhm) ± $(fwhm_err)" + + result = merge(NamedTuple{keys(v_ml)}([measurement(v_ml[k], v_ml_err[k]) for k in keys(v_ml)]...), + (fwhm = measurement(fwhm, fwhm_err),), NamedTuple{(:gof_survived, :gof_cut)}(gofs)) + + else + # get fwhm of peak + fwhm, fwhm_err = get_peak_fwhm_th228(v_ml, v_ml, false) + + @debug "Best Fit values" + @debug "μ: $(v_ml.μ)" + @debug "σ: $(v_ml.σ)" + @debug "n: $(v_ml.n)" + @debug "FWHM: $(fwhm)" + + result = merge(NamedTuple{keys(v_ml)}([measurement(v_ml[k], NaN) for k in keys(v_ml)]...), + (fwhm = measurement(fwhm, NaN), gof_survived = NamedTuple(), gof_cut = NamedTuple())) + end + + report = ( + survived = ( + v = v_ml_survived, + h = h_survived, + f_fit = x -> Base.Fix2(th228_fit_functions.f_fit, v_ml_survived)(x), + f_sig = x -> Base.Fix2(th228_fit_functions.f_sig, v_ml_survived)(x), + f_lowEtail = x -> Base.Fix2(th228_fit_functions.f_lowEtail, v_ml_survived)(x), + f_bck = x -> Base.Fix2(th228_fit_functions.f_bck, v_ml_survived)(x), + gof = result.gof_survived + ), + cut = ( + v = v_ml_cut, + h = h_cut, + f_fit = x -> Base.Fix2(th228_fit_functions.f_fit, v_ml_cut)(x), + f_sig = x -> Base.Fix2(th228_fit_functions.f_sig, v_ml_cut)(x), + f_lowEtail = x -> Base.Fix2(th228_fit_functions.f_lowEtail, v_ml_cut)(x), + f_bck = x -> Base.Fix2(th228_fit_functions.f_bck, v_ml_cut)(x), + gof = result.gof_cut + ), + sf = v_ml.sf, + bsf = v_ml.bsf + ) + + return result, report +end \ No newline at end of file From ae1b9058fb298020ca9a985260598ba8d9ada963 Mon Sep 17 00:00:00 2001 From: Felix Hagemann Date: Sun, 28 Apr 2024 23:18:02 +0200 Subject: [PATCH 48/68] Add plot recipe for survival fraction report --- ext/LegendSpecFitsRecipesBaseExt.jl | 25 +++++++++++++++++++++++++ 1 file changed, 25 insertions(+) diff --git a/ext/LegendSpecFitsRecipesBaseExt.jl b/ext/LegendSpecFitsRecipesBaseExt.jl index 9c66d582..df33c87e 100644 --- a/ext/LegendSpecFitsRecipesBaseExt.jl +++ b/ext/LegendSpecFitsRecipesBaseExt.jl @@ -321,6 +321,31 @@ end end end +@recipe function f(report::NamedTuple{(:survived, :cut, :sf, :bsf)}; cut_value = missing) + size --> (1200,400) + left_margin --> (10, :mm) + layout := @layout ([A{0.01h}; [B{0.75h}; C{0.175h}] [D{0.7h}; E{0.175h}]]) + @series begin + title := (ismissing(cut_value) ? "" : "Cut value: $(cut_value) ") * "Survival fraction: $(round(report.sf * 100, digits = 2))%" + grid := false + showaxis := false + label := nothing + bottom_margin := (-20, :px) + [] + end + @series begin + title := "Survived\n\n" + _subplot := 2 + report.survived + end + @series begin + title := "Cut\n\n" + _subplot := 4 + report.cut + end +end + + @recipe function f(x::Vector{T}, cut::NamedTuple{(:low, :high, :max), Tuple{T, T, T}}) where T<:Unitful.RealOrRealQuantity ylabel := "Counts" legend := :topright From b2c138d7a25b1b39a61b4f44ecbb88ab64ba4571 Mon Sep 17 00:00:00 2001 From: Felix Hagemann Date: Mon, 29 Apr 2024 13:13:17 +0200 Subject: [PATCH 49/68] Reduce max skew_fraction from 0.25 to 0.1 --- src/specfit.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/specfit.jl b/src/specfit.jl index 176786a3..12c92ba9 100644 --- a/src/specfit.jl +++ b/src/specfit.jl @@ -167,7 +167,7 @@ function fit_single_peak_th228(h::Histogram, ps::NamedTuple{(:peak_pos, :peak_fw σ = weibull_from_mx(ps.peak_sigma, 2*ps.peak_sigma), n = weibull_from_mx(ps.peak_counts, 2*ps.peak_counts), step_amplitude = weibull_from_mx(ps.mean_background_step, ps.mean_background_step + 5*ps.mean_background_std), - skew_fraction = ifelse(low_e_tail, truncated(weibull_from_mx(0.01, 0.05), 0.0, 0.25), ConstValueDist(0.0)), + skew_fraction = ifelse(low_e_tail, truncated(weibull_from_mx(0.01, 0.05), 0.0, 0.1), ConstValueDist(0.0)), skew_width = ifelse(low_e_tail, weibull_from_mx(0.001, 1e-2), ConstValueDist(1.0)), background = weibull_from_mx(ps.mean_background, ps.mean_background + 5*ps.mean_background_std), ) From 13045bd117194642b449efc4b1b241f06b2b03dd Mon Sep 17 00:00:00 2001 From: Felix Hagemann Date: Mon, 29 Apr 2024 22:43:21 +0200 Subject: [PATCH 50/68] Update `fit_subpeaks_th228` --- src/specfit.jl | 66 ++++++++++++++++++++++++++------------------------ 1 file changed, 34 insertions(+), 32 deletions(-) diff --git a/src/specfit.jl b/src/specfit.jl index 12c92ba9..0294420b 100644 --- a/src/specfit.jl +++ b/src/specfit.jl @@ -337,24 +337,26 @@ Also, FWHM is calculated from the fitted peakshape with MC error propagation. Th """ function fit_subpeaks_th228( h_survived::Histogram, h_cut::Histogram, h_result; - uncertainty::Bool=false, low_e_tail::Bool=true, pseudo_prior::NamedTupleDist=NamedTupleDist(empty = true), fit_func::Symbol=:f_fit + uncertainty::Bool=false, low_e_tail::Bool=true, fix_σ::Bool = true, fix_skew_fraction::Bool = true, fix_skew_width::Bool = true, + pseudo_prior::NamedTupleDist=NamedTupleDist(empty = true), fit_func::Symbol=:f_fit ) # create standard pseudo priors standard_pseudo_prior = let ps = h_result, ps_cut = estimate_single_peak_stats(h_cut), ps_survived = estimate_single_peak_stats(h_survived) NamedTupleDist( - μ = ConstValueDist(Measurements.value(ps.μ)), - σ = ConstValueDist(Measurements.value(ps.σ)), - n = ConstValueDist(Measurements.value(ps.n)), - sf = Uniform(0,1), - bsf = Uniform(0,1), - step_amplitude_survived = weibull_from_mx(ps_survived.mean_background_step, ps_survived.mean_background_step + 5*ps_survived.mean_background_std), - skew_fraction_survived = ifelse(low_e_tail, truncated(weibull_from_mx(0.01, 0.05), 0.0, 0.1), ConstValueDist(0.0)), - skew_width_survived = ifelse(low_e_tail, weibull_from_mx(0.001, 1e-2), ConstValueDist(1.0)), - step_amplitude_cut = weibull_from_mx(ps_cut.mean_background_step, ps_cut.mean_background_step + 5*ps_cut.mean_background_std), - skew_fraction_cut = ifelse(low_e_tail, truncated(weibull_from_mx(0.01, 0.05), 0.0, 0.1), ConstValueDist(0.0)), - skew_width_cut = ifelse(low_e_tail, weibull_from_mx(0.001, 1e-2), ConstValueDist(1.0)), - background = ConstValueDist(Measurements.value(ps.background)) + μ = ConstValueDist(mvalue(ps.μ)), + σ_survived = ifelse(fix_σ, ConstValueDist(mvalue(ps.σ)), weibull_from_mx(mvalue(ps.σ), 2*mvalue(ps.σ))), + σ_cut = ifelse(fix_σ, ConstValueDist(mvalue(ps.σ)), weibull_from_mx(mvalue(ps.σ), 2*mvalue(ps.σ))), + n = ConstValueDist(mvalue(ps.n)), + sf = Uniform(0,1), # signal survival fraction + bsf = Uniform(0,1), # background survival fraction + sasf = Uniform(0,1), # step amplitude survival fraction + step_amplitude = ConstValueDist(mvalue(ps.step_amplitude)), + skew_fraction_survived = ifelse(low_e_tail, ifelse(fix_skew_fraction, ConstValueDist(mvalue(ps.skew_fraction)), truncated(weibull_from_mx(0.01, 0.05), 0.0, 0.1)), ConstValueDist(0.0)), + skew_fraction_cut = ifelse(low_e_tail, ifelse(fix_skew_fraction, ConstValueDist(mvalue(ps.skew_fraction)), truncated(weibull_from_mx(0.01, 0.05), 0.0, 0.1)), ConstValueDist(0.0)), + skew_width_survived = ifelse(low_e_tail, ifelse(fix_skew_width, mvalue(ps.skew_width), weibull_from_mx(0.001, 1e-2)), ConstValueDist(1.0)), + skew_width_cut = ifelse(low_e_tail, ifelse(fix_skew_width, mvalue(ps.skew_width), weibull_from_mx(0.001, 1e-2)), ConstValueDist(1.0)), + background = ConstValueDist(mvalue(ps.background)) ) end @@ -378,14 +380,14 @@ function fit_subpeaks_th228( # create loglikehood function: f_loglike(v) that can be evaluated for any set of v (fit parameter) f_loglike = let f_fit=th228_fit_functions[fit_func], h_cut=h_cut, h_survived=h_survived v -> begin - v_survived = (μ = v.μ, σ = v.σ, n = v.n * v.sf, - step_amplitude = v.step_amplitude_survived, + v_survived = (μ = v.μ, σ = v.σ_survived, n = v.n * v.sf, + step_amplitude = v.step_amplitude * v.sasf, skew_fraction = v.skew_fraction_survived, skew_width = v.skew_width_survived, background = v.background * v.bsf ) - v_cut = (μ = v.μ, σ = v.σ, n = v.n * (1 - v.sf), - step_amplitude = v.step_amplitude_cut, + v_cut = (μ = v.μ, σ = v.σ_cut, n = v.n * (1 - v.sf), + step_amplitude = v.step_amplitude * (1 - v.sasf), skew_fraction = v.skew_fraction_cut, skew_width = v.skew_width_cut, background = v.background * (1 - v.bsf) @@ -402,9 +404,9 @@ function fit_subpeaks_th228( v_ml_survived = ( μ = v_ml.μ, - σ = v_ml.σ, + σ = v_ml.σ_survived, n = v_ml.n * v_ml.sf, - step_amplitude = v_ml.step_amplitude_survived, + step_amplitude = v_ml.step_amplitude * v_ml.sasf, skew_fraction = v_ml.skew_fraction_survived, skew_width = v_ml.skew_width_survived, background = v_ml.background * v_ml.bsf @@ -412,9 +414,9 @@ function fit_subpeaks_th228( v_ml_cut = ( μ = v_ml.μ, - σ = v_ml.σ, + σ = v_ml.σ_cut, n = v_ml.n * (1 - v_ml.sf), - step_amplitude = v_ml.step_amplitude_cut, + step_amplitude = v_ml.step_amplitude * (1 - v_ml.sasf), skew_fraction = v_ml.skew_fraction_cut, skew_width = v_ml.skew_width_cut, background = v_ml.background * (1 - v_ml.bsf) @@ -469,12 +471,12 @@ function fit_subpeaks_th228( get_peak_fwhm_th228(v_ml, v_ml_err) end - @debug "Best Fit values for $(part)" - @debug "μ: $(v_ml.μ) ± $(v_ml_err.μ)" - @debug "σ: $(v_ml.σ) ± $(v_ml_err.σ)" - @debug "n: $(v_ml.n) ± $(v_ml_err.n)" - @debug "p: $pval , chi2 = $(chi2) with $(dof) dof" - @debug "FWHM: $(fwhm) ± $(fwhm_err)" + # @debug "Best Fit values for $(part)" + # @debug "μ: $(v_ml.μ) ± $(v_ml_err.μ)" + # @debug "σ: $(v_ml.σ) ± $(v_ml_err.σ)" + # @debug "n: $(v_ml.n) ± $(v_ml_err.n)" + # @debug "p: $pval , chi2 = $(chi2) with $(dof) dof" + # @debug "FWHM: $(fwhm) ± $(fwhm_err)" result = merge(NamedTuple{keys(v_ml)}([measurement(v_ml[k], v_ml_err[k]) for k in keys(v_ml)]...), (fwhm = measurement(fwhm, fwhm_err),), NamedTuple{(:gof_survived, :gof_cut)}(gofs)) @@ -483,11 +485,11 @@ function fit_subpeaks_th228( # get fwhm of peak fwhm, fwhm_err = get_peak_fwhm_th228(v_ml, v_ml, false) - @debug "Best Fit values" - @debug "μ: $(v_ml.μ)" - @debug "σ: $(v_ml.σ)" - @debug "n: $(v_ml.n)" - @debug "FWHM: $(fwhm)" + # @debug "Best Fit values" + # @debug "μ: $(v_ml.μ)" + # @debug "σ: $(v_ml.σ)" + # @debug "n: $(v_ml.n)" + # @debug "FWHM: $(fwhm)" result = merge(NamedTuple{keys(v_ml)}([measurement(v_ml[k], NaN) for k in keys(v_ml)]...), (fwhm = measurement(fwhm, NaN), gof_survived = NamedTuple(), gof_cut = NamedTuple())) From d2f12615f63e14848cb48d79d6e8111b292018b8 Mon Sep 17 00:00:00 2001 From: Felix Hagemann Date: Mon, 29 Apr 2024 22:44:12 +0200 Subject: [PATCH 51/68] Fix bug when passing pseudo priors --- src/specfit.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/specfit.jl b/src/specfit.jl index 0294420b..5137c701 100644 --- a/src/specfit.jl +++ b/src/specfit.jl @@ -175,7 +175,7 @@ function fit_single_peak_th228(h::Histogram, ps::NamedTuple{(:peak_pos, :peak_fw # use standard priors in case of no overwrites given if !(:empty in keys(pseudo_prior)) # check if input overwrite prior has the same fields as the standard prior set - @assert all(f -> f in keys(standard_pseudo_prior), keys(standard_pseudo_prior)) "Pseudo priors can only have $(keys(standard_pseudo_prior)) as fields." + @assert all(f -> f in keys(standard_pseudo_prior), keys(pseudo_prior)) "Pseudo priors can only have $(keys(standard_pseudo_prior)) as fields." # replace standard priors with overwrites pseudo_prior = merge(standard_pseudo_prior, pseudo_prior) else @@ -363,7 +363,7 @@ function fit_subpeaks_th228( # use standard priors in case of no overwrites given if !(:empty in keys(pseudo_prior)) # check if input overwrite prior has the same fields as the standard prior set - @assert all(f -> f in keys(standard_pseudo_prior), keys(standard_pseudo_prior)) "Pseudo priors can only have $(keys(standard_pseudo_prior)) as fields." + @assert all(f -> f in keys(standard_pseudo_prior), keys(pseudo_prior)) "Pseudo priors can only have $(keys(standard_pseudo_prior)) as fields." # replace standard priors with overwrites pseudo_prior = merge(standard_pseudo_prior, pseudo_prior) else From d63f6ea7508f350da64cb4fda9f869a2e718135b Mon Sep 17 00:00:00 2001 From: Felix Hagemann Date: Mon, 29 Apr 2024 22:59:15 +0200 Subject: [PATCH 52/68] Make plots of `fit_subpeaks_th228` report symmetric --- ext/LegendSpecFitsRecipesBaseExt.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ext/LegendSpecFitsRecipesBaseExt.jl b/ext/LegendSpecFitsRecipesBaseExt.jl index df33c87e..0faba9ee 100644 --- a/ext/LegendSpecFitsRecipesBaseExt.jl +++ b/ext/LegendSpecFitsRecipesBaseExt.jl @@ -324,7 +324,7 @@ end @recipe function f(report::NamedTuple{(:survived, :cut, :sf, :bsf)}; cut_value = missing) size --> (1200,400) left_margin --> (10, :mm) - layout := @layout ([A{0.01h}; [B{0.75h}; C{0.175h}] [D{0.7h}; E{0.175h}]]) + layout := @layout ([A{0.01h}; [B{0.75h}; C{0.175h}] [D{0.75h}; E{0.175h}]]) @series begin title := (ismissing(cut_value) ? "" : "Cut value: $(cut_value) ") * "Survival fraction: $(round(report.sf * 100, digits = 2))%" grid := false From d4c936b9e99b882763bd9709f4b4940d17e95400 Mon Sep 17 00:00:00 2001 From: Felix Hagemann Date: Mon, 29 Apr 2024 23:28:27 +0200 Subject: [PATCH 53/68] Add dummy test --- test/Project.toml | 9 ++++++--- test/test_specfit.jl | 43 +++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 49 insertions(+), 3 deletions(-) diff --git a/test/Project.toml b/test/Project.toml index 5ffde185..efe6f101 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,6 +1,6 @@ [deps] Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" -# BAT = "c0cd4b16-88b7-57fa-983b-ab80aecada7e" +Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" LegendDataTypes = "99e09c13-5545-5ee2-bfa2-77f358fb75d8" @@ -9,12 +9,15 @@ Measurements = "eff96d63-e80a-5855-80a2-b1b0885c5ab7" Optim = "429524aa-4258-5aef-a3af-852621145aeb" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" +StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" [compat] +Distributions = "0.24, 0.25" Documenter = "1" -Interpolations = "0.15" -LegendDataTypes = "0.1" +Interpolations = "0.15" +LegendDataTypes = "0.1" Measurements = "2" +StatsBase = "0.32, 0.33, 0.34" Unitful = "1" \ No newline at end of file diff --git a/test/test_specfit.jl b/test/test_specfit.jl index f20d7fd8..64bfa557 100644 --- a/test/test_specfit.jl +++ b/test/test_specfit.jl @@ -2,6 +2,9 @@ using LegendSpecFits using Test using LegendDataTypes: fast_flatten +using Measurements: value as mvalue +using Distributions +using StatsBase using Interpolations using Unitful include("test_utils.jl") @@ -17,3 +20,43 @@ include("test_utils.jl") # fit result, report = fit_peaks(result_simple.peakhists, result_simple.peakstats, ustrip.(th228_lines),; uncertainty=true,calib_type = :th228); end + +@testset "fit_subpeaks_th228" begin + + # Compose a dummy energy peak (E) + n = 100000 + b = 20000 + μ = 1000.0 + σ = 20.0 + E = vcat(μ .+ σ .* randn(n), rand(Uniform(μ - 10*σ, μ + 10*σ), b)) + + # Fit the initial histogram and check that the determined quantities match the given ones + bins = range(μ - 10*σ, stop = μ + 10*σ, length = round(Int, sqrt(b))) + h = fit(Histogram, E, bins) + h_result, h_report = fit_single_peak_th228(h, estimate_single_peak_stats(h), low_e_tail = false) + @testset "fit single_peak_th228" begin + @test isapprox(mvalue(h_result.μ), μ, rtol = 0.01) + @test isapprox(mvalue(h_result.σ), σ, rtol = 0.01) + @test isapprox(mvalue(h_result.n), n, rtol = 0.01) + @test isapprox(mvalue(h_result.background), b / (20σ), rtol = 0.1) + end + + # Introduce a dummy quantity (Q) which has two peaks (around 0 and around 1) + # that should have 90% of the signal and 50% of the background around 0 + # and 10% of the signal and 50% of the background around 1 + sf = 0.9 + bsf = 0.5 + Q = vcat(0.2 * randn(round(Int, n*sf)), 0.2 * randn(round(Int, n*(1-sf))) .+ 1, 0.2 * randn(round(Int, b*bsf)), 0.2 * randn(round(Int, b*(1-bsf))) .+ 1) + + # Perform the cut on this dummy Q by selecting events with Q < 0.5 (the peak around 0) + # and check that the survival fractions for signal and background agree with what was initially given + cut_value = 0.5 + h_survived = fit(Histogram, E[findall(Q .<= cut_value)], bins) + h_cut = fit(Histogram, E[findall(Q .> cut_value)], bins) + result, report = LegendSpecFits.fit_subpeaks_th228(h_survived, h_cut, h_result, uncertainty = true) + @testset "fit_subpeaks_th228" begin + @test isapprox(mvalue(result.sf), sf, rtol = 0.01) + @test isapprox(mvalue(result.bsf), bsf, rtol = 0.01) + end + +end From a3d28b3a69590d6142fdb892fa5f72e9c94d1d96 Mon Sep 17 00:00:00 2001 From: Felix Hagemann Date: Tue, 30 Apr 2024 09:18:33 +0200 Subject: [PATCH 54/68] Relax constraints on `Int64` in `get_mc_value_shapes` --- src/utils.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/utils.jl b/src/utils.jl index fc6b0b7c..065d2888 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -67,11 +67,11 @@ end """ - get_mc_value_shapes(v::NamedTuple, v_err::NamedTuple, n::Int64) + get_mc_value_shapes(v::NamedTuple, v_err::NamedTuple, n::Union{Int64,Int32}) Return a `NamedTuple` with the same fields as `v` and `v_err` but with `Normal` distributions for each field. """ -function get_mc_value_shapes(v::NamedTuple, v_err::NamedTuple, n::Int64) +function get_mc_value_shapes(v::NamedTuple, v_err::NamedTuple, n::Union{Int64,Int32}) vs = BAT.distprod(map(Normal, v, v_err)) NamedTuple.(rand(vs, n)) end From ad8ebc56ec523af56c2922b8329a3ac862c3f256 Mon Sep 17 00:00:00 2001 From: Florian Henkes Date: Wed, 8 May 2024 23:39:05 +0200 Subject: [PATCH 55/68] single git uncertainties use abs values to avoid numerical issues --- src/singlefit.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/singlefit.jl b/src/singlefit.jl index 3de23950..5e9a90cb 100644 --- a/src/singlefit.jl +++ b/src/singlefit.jl @@ -47,8 +47,8 @@ function fit_single_trunc_gauss(x::Vector{T}, cuts::NamedTuple{(:low, :high, :ma param_covariance = inv(H) # Extract the parameter uncertainties - μ_uncertainty = sqrt(param_covariance[1, 1]) - σ_uncertainty = sqrt(param_covariance[2, 2]) + μ_uncertainty = sqrt(abs(param_covariance[1, 1])) + σ_uncertainty = sqrt(abs(param_covariance[2, 2])) @debug "μ: $μ ± $μ_uncertainty" @debug "σ: $σ ± $σ_uncertainty" From 7c4ef594dda2ead9afe013e61a0463e61156b33f Mon Sep 17 00:00:00 2001 From: Felix Hagemann Date: Thu, 9 May 2024 17:13:03 +0200 Subject: [PATCH 56/68] Update determination of A/E cut survival fraction --- src/aoe_calibration.jl | 76 ++++++++++++++++++++++++++++++-------- src/filter_optimization.jl | 11 +++--- 2 files changed, 67 insertions(+), 20 deletions(-) diff --git a/src/aoe_calibration.jl b/src/aoe_calibration.jl index e4f2b3de..6a289960 100644 --- a/src/aoe_calibration.jl +++ b/src/aoe_calibration.jl @@ -90,10 +90,9 @@ export prepare_dep_peakhist """ get_n_after_aoe_cut(aoe_cut::T, aoe::Array{T}, e::Array{T}, peak::T, window::Array{T}, bin_width::T, result_before::NamedTuple, peakstats::NamedTuple; uncertainty=true) where T<:Real -Get the number of counts after a AoE cut value `aoe_cut` for a given `peak` and `window` size whiile performing a peak fit with fixed position. The number of counts is determined by fitting the peak with a pseudo prior for the peak position. +Get the number of counts after a AoE cut value `aoe_cut` for a given `peak` and `window` size while performing a peak fit with fixed position. The number of counts is determined by fitting the peak with a pseudo prior for the peak position. # Returns - `n`: Number of counts after the cut -- `n_err`: Uncertainty of the number of counts after the cut """ function get_n_after_aoe_cut(aoe_cut::Unitful.RealOrRealQuantity, aoe::Vector{<:Unitful.RealOrRealQuantity}, e::Vector{<:T}, peak::T, window::Vector{T}, bin_width::T, result_before::NamedTuple, peakstats::NamedTuple; uncertainty::Bool=true, fixed_position::Bool=true) where T<:Unitful.Energy{<:Real} # get energy after cut and create histogram @@ -111,6 +110,26 @@ end export get_n_after_aoe_cut +""" + get_sf_after_aoe_cut(aoe_cut::T, aoe::Array{T}, e::Array{T}, peak::T, window::Array{T}, bin_width::T, result_before::NamedTuple; uncertainty=true) where T<:Real + +Get the survival fraction after a AoE cut value `aoe_cut` for a given `peak` and `window` size from a combined fit to the survived and cut histograms. + +# Returns +- `sf`: Survival fraction after the cut +""" +function get_sf_after_aoe_cut(aoe_cut::Unitful.RealOrRealQuantity, aoe::Vector{<:Unitful.RealOrRealQuantity}, e::Vector{<:T}, peak::T, window::Vector{T}, bin_width::T, result_before::NamedTuple; uncertainty::Bool=true) where T<:Unitful.Energy{<:Real} + # get energy after cut and create histogram + survived = fit(Histogram, ustrip.(e[aoe .>= aoe_cut]), ustrip(peak-first(window):bin_width:peak+last(window))) + cut = fit(Histogram, ustrip.(e[aoe .< aoe_cut]), ustrip(peak-first(window):bin_width:peak+last(window))) + # fit peak and return number of signal counts + result, _ = fit_subpeaks_th228(survived, cut, result_before; uncertainty=uncertainty, low_e_tail=false) + return result.sf +end +export get_sf_after_aoe_cut + + + """ get_aoe_cut(aoe::Vector{<:Unitful.RealOrRealQuantity}, e::Vector{<:T},; dep::T=1592.53u"keV", window::Vector{T}=[12.0, 10.0]u"keV", dep_sf::Float64=0.9, cut_search_interval::Tuple{Quantity{<:Real}, Quantity{<:Real}}=(-25.0u"keV^-1", 0.0u"keV^-1"), rtol::Float64=0.001, bin_width_window::T=3.0u"keV", fixed_position::Bool=true, sigma_high_sided::Float64=NaN, uncertainty::Bool=true, maxiters::Int=200) where T<:Unitful.Energy{<:Real} @@ -121,6 +140,7 @@ The algorhithm utilizes a root search algorithm to find the cut value with a rel - `n0`: Number of counts before the cut - `nsf`: Number of counts after the cut """ +#= function get_aoe_cut(aoe::Vector{<:Unitful.RealOrRealQuantity}, e::Vector{<:T},; dep::T=1592.53u"keV", window::Vector{<:T}=[12.0, 10.0]u"keV", dep_sf::Float64=0.9, cut_search_interval::Tuple{<:Unitful.RealOrRealQuantity, <:Unitful.RealOrRealQuantity}=(-25.0*unit(first(aoe)), 1.0*unit(first(aoe))), rtol::Float64=0.001, bin_width_window::T=3.0u"keV", fixed_position::Bool=true, sigma_high_sided::Float64=NaN, uncertainty::Bool=true, maxiters::Int=200) where T<:Unitful.Energy{<:Real} # check if a high sided AoE cut should be applied before the AoE cut is generated if !isnan(sigma_high_sided) @@ -149,6 +169,35 @@ function get_aoe_cut(aoe::Vector{<:Unitful.RealOrRealQuantity}, e::Vector{<:T},; return (lowcut = measurement(aoe_cut, aoe_cut * rtol), highcut = sigma_high_sided, n0 = result_before.n, nsf = nsf, sf = nsf / result_before.n * 100*u"percent") end export get_aoe_cut +=# + +function get_aoe_cut(aoe::Vector{<:Unitful.RealOrRealQuantity}, e::Vector{<:T},; dep::T=1592.53u"keV", window::Vector{<:T}=[12.0, 10.0]u"keV", dep_sf::Float64=0.9, cut_search_interval::Tuple{<:Unitful.RealOrRealQuantity, <:Unitful.RealOrRealQuantity}=(-25.0*unit(first(aoe)), 1.0*unit(first(aoe))), rtol::Float64=0.001, bin_width_window::T=3.0u"keV", fixed_position::Bool=true, sigma_high_sided::Float64=NaN, uncertainty::Bool=true, maxiters::Int=200) where T<:Unitful.Energy{<:Real} + # check if a high sided AoE cut should be applied before the AoE cut is generated + if !isnan(sigma_high_sided) + e = e[aoe .< sigma_high_sided] + aoe = aoe[aoe .< sigma_high_sided] + end + # cut window around peak + aoe = aoe[dep-first(window) .< e .< dep+last(window)] + e = e[dep-first(window) .< e .< dep+last(window)] + # estimate bin width + bin_width = get_friedman_diaconis_bin_width(e[dep - bin_width_window .< e .< dep + bin_width_window]) + # create histogram + dephist = fit(Histogram, ustrip.(e), ustrip(dep-first(window):bin_width:dep+last(window))) + # get peakstats + depstats = estimate_single_peak_stats(dephist) + # fit before cut + result_before, _ = fit_single_peak_th228(dephist, depstats; uncertainty=uncertainty, fixed_position=fixed_position, low_e_tail=false) + # get aoe cut + sf_dep_f = cut -> get_sf_after_aoe_cut(cut, aoe, e, dep, window, bin_width, mvalue(result_before); uncertainty=false) - dep_sf + aoe_cut = find_zero(sf_dep_f, cut_search_interval, Bisection(), rtol=rtol, maxiters=maxiters) + # get sf after cut + sf = get_sf_after_aoe_cut(aoe_cut, aoe, e, dep, window, bin_width, mvalue(result_before); uncertainty=uncertainty) + @info sf + return (lowcut = measurement(aoe_cut, aoe_cut * rtol), highcut = sigma_high_sided, n0 = result_before.n, nsf = result_before.n * sf, sf = sf * 100*u"percent") +end +export get_aoe_cut + """ @@ -175,40 +224,37 @@ function get_peak_surrival_fraction(aoe::Vector{<:Unitful.RealOrRealQuantity}, e # get e after cut if !isnan(sigma_high_sided) - e = e[aoe_cut .< aoe .< sigma_high_sided] + e = e[aoe .< sigma_high_sided] else - e = e[aoe_cut .< aoe] + e_survived = e[aoe_cut .<= aoe] + e_cut = e[aoe_cut .> aoe] end # estimate bin width bin_width = get_friedman_diaconis_bin_width(e[e .> peak - bin_width_window .&& e .< peak + bin_width_window]) # get energy after cut and create histogram - peakhist = fit(Histogram, ustrip(e), ustrip(peak-first(window):bin_width:peak+last(window))) - # create pseudo_prior with known peak sigma in signal for more stable fit - # pseudo_prior = NamedTupleDist(μ = ConstValueDist(result_before.μ), σ = Normal(result_before.σ, 0.1)) - pseudo_prior = NamedTupleDist(σ = Normal(result_before.σ, 0.1), ) - # estimate peak stats - peakstats = estimate_single_peak_stats(peakhist) + survived = fit(Histogram, ustrip(e_survived), ustrip(peak-first(window):bin_width:peak+last(window))) + cut = fit(Histogram, ustrip(e_cut), ustrip(peak-first(window):bin_width:peak+last(window))) # fit peak and return number of signal counts - result_after, report_after = fit_single_peak_th228(peakhist, peakstats,; uncertainty=uncertainty, low_e_tail=low_e_tail) + result_after, report_after = fit_subpeaks_th228(survived, cut, result_before; uncertainty=uncertainty, low_e_tail=low_e_tail) # result_after, report_after = fit_single_peak_th228(peakhist, peakstats,; uncertainty=uncertainty, low_e_tail=low_e_tail, pseudo_prior=pseudo_prior) # calculate surrival fraction - sf = result_after.n / result_before.n * 100.0 * u"percent" + sf = result_after.sf * 100u"percent" result = ( peak = peak, n_before = result_before.n, - n_after = result_after.n, + n_after = result_before.n * result_after.sf, sf = sf ) report = ( peak = result.peak, n_before = result.n_before, - n_after = result.n_after, + n_after = result.n_before * result_after.sf, sf = result.sf, before = report_before, after = report_after, ) - return result, report + return result, report_after end export get_peak_surrival_fraction diff --git a/src/filter_optimization.jl b/src/filter_optimization.jl index cffb1477..6746e0fe 100644 --- a/src/filter_optimization.jl +++ b/src/filter_optimization.jl @@ -166,7 +166,7 @@ Fit the SG window length for the SEP data and return the optimal window length a - `result`: optimal window length and corresponding survival fraction - `report`: report with all window lengths and survival fractions """ -function fit_sg_wl(dep_sep_data::NamedTuple{(:dep, :sep)}, a_grid_wl_sg::StepRangeLen, optimization_config::PropDict) +function fit_sg_wl(dep_sep_data::NamedTuple{(:dep, :sep)}, a_grid_wl_sg::StepRangeLen, optimization_config::PropDict; uncertainty::Bool = false) # unpack config dep, dep_window = optimization_config.dep, Float64.(optimization_config.dep_window) sep, sep_window = optimization_config.sep, Float64.(optimization_config.sep_window) @@ -176,7 +176,7 @@ function fit_sg_wl(dep_sep_data::NamedTuple{(:dep, :sep)}, a_grid_wl_sg::StepRan aoe_dep, aoe_sep = dep_sep_data.dep.aoe, dep_sep_data.sep.aoe # prepare peakhist - result_dep, _ = prepare_dep_peakhist(e_dep, dep; n_bins_cut=optimization_config.nbins_dep_cut, relative_cut=optimization_config.dep_rel_cut, uncertainty=false) + result_dep, _ = prepare_dep_peakhist(e_dep, dep; n_bins_cut=optimization_config.nbins_dep_cut, relative_cut=optimization_config.dep_rel_cut, uncertainty=uncertainty) yield() @@ -186,6 +186,7 @@ function fit_sg_wl(dep_sep_data::NamedTuple{(:dep, :sep)}, a_grid_wl_sg::StepRan # create empty arrays for sf and sf_err sep_sfs = Quantity{Measurement}[] + sep_reports = [] wls = Quantity{<:Real}[] # for each window lenght, calculate the survival fraction in the SEP @@ -199,15 +200,15 @@ function fit_sg_wl(dep_sep_data::NamedTuple{(:dep, :sep)}, a_grid_wl_sg::StepRan min_aoe_dep_i = quantile(aoe_dep_i, optimization_config.min_aoe_quantile) + optimization_config.min_aoe_offset aoe_dep_i_hist = fit(Histogram, ustrip.(aoe_dep_i), ustrip(min_aoe_dep_i):ustrip(get_friedman_diaconis_bin_width(aoe_dep_i[min_aoe_dep_i .< aoe_dep_i .< max_aoe_dep_i])):ustrip(max_aoe_dep_i)) - max_aoe_dep_i = first(aoe_dep_i_hist.edges)[argmax(aoe_dep_i_hist.weights)]*unit(max_aoe_dep_i) + max_aoe_dep_i = first(aoe_dep_i_hist.edges)[min(end, argmax(aoe_dep_i_hist.weights)+1)] * unit(max_aoe_dep_i) try - psd_cut = get_psd_cut(aoe_dep_i, e_dep_i; window=dep_window, cut_search_interval=(min_aoe_dep_i, max_aoe_dep_i), uncertainty=false) + psd_cut = get_aoe_cut(aoe_dep_i, e_dep_i; window=dep_window, cut_search_interval=(min_aoe_dep_i, max_aoe_dep_i), uncertainty=uncertainty) aoe_sep_i = aoe_sep[i_aoe, :][isfinite.(aoe_sep[i_aoe, :])] ./ result_dep.m_calib e_sep_i = e_sep_calib[isfinite.(aoe_sep[i_aoe, :])] - result_sep, _ = get_peak_surrival_fraction(aoe_sep_i, e_sep_i, sep, sep_window, psd_cut.cut; uncertainty=false, low_e_tail=false) + result_sep, _ = get_peak_surrival_fraction(aoe_sep_i, e_sep_i, sep, sep_window, psd_cut.lowcut; uncertainty=uncertainty, low_e_tail=false) push!(sep_sfs, result_sep.sf) push!(wls, wl) From cc878e40ed82a803304da7167b0a628453622e50 Mon Sep 17 00:00:00 2001 From: Felix Hagemann Date: Thu, 9 May 2024 17:19:15 +0200 Subject: [PATCH 57/68] f --- src/aoe_calibration.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/aoe_calibration.jl b/src/aoe_calibration.jl index 6a289960..8540f957 100644 --- a/src/aoe_calibration.jl +++ b/src/aoe_calibration.jl @@ -224,6 +224,7 @@ function get_peak_surrival_fraction(aoe::Vector{<:Unitful.RealOrRealQuantity}, e # get e after cut if !isnan(sigma_high_sided) + # TODO: decide how to deal with the high sided cut! e = e[aoe .< sigma_high_sided] else e_survived = e[aoe_cut .<= aoe] @@ -254,7 +255,7 @@ function get_peak_surrival_fraction(aoe::Vector{<:Unitful.RealOrRealQuantity}, e before = report_before, after = report_after, ) - return result, report_after + return result, report end export get_peak_surrival_fraction From 3c5218af3e4ecd659a4614adf81721b4000ed56d Mon Sep 17 00:00:00 2001 From: Felix Hagemann Date: Thu, 9 May 2024 17:34:36 +0200 Subject: [PATCH 58/68] Remove printing survival fractions as @info --- src/aoe_calibration.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/aoe_calibration.jl b/src/aoe_calibration.jl index 8540f957..ecc1f528 100644 --- a/src/aoe_calibration.jl +++ b/src/aoe_calibration.jl @@ -193,7 +193,6 @@ function get_aoe_cut(aoe::Vector{<:Unitful.RealOrRealQuantity}, e::Vector{<:T},; aoe_cut = find_zero(sf_dep_f, cut_search_interval, Bisection(), rtol=rtol, maxiters=maxiters) # get sf after cut sf = get_sf_after_aoe_cut(aoe_cut, aoe, e, dep, window, bin_width, mvalue(result_before); uncertainty=uncertainty) - @info sf return (lowcut = measurement(aoe_cut, aoe_cut * rtol), highcut = sigma_high_sided, n0 = result_before.n, nsf = result_before.n * sf, sf = sf * 100*u"percent") end export get_aoe_cut From cf2540bed12ac6145d76ace2482fac8eade97e35 Mon Sep 17 00:00:00 2001 From: Florian Henkes Date: Thu, 9 May 2024 19:25:13 +0200 Subject: [PATCH 59/68] Fixed sg window length optimization --- src/filter_optimization.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/filter_optimization.jl b/src/filter_optimization.jl index 6746e0fe..6b1f4756 100644 --- a/src/filter_optimization.jl +++ b/src/filter_optimization.jl @@ -168,8 +168,8 @@ Fit the SG window length for the SEP data and return the optimal window length a """ function fit_sg_wl(dep_sep_data::NamedTuple{(:dep, :sep)}, a_grid_wl_sg::StepRangeLen, optimization_config::PropDict; uncertainty::Bool = false) # unpack config - dep, dep_window = optimization_config.dep, Float64.(optimization_config.dep_window) - sep, sep_window = optimization_config.sep, Float64.(optimization_config.sep_window) + dep, dep_window = optimization_config.dep, optimization_config.dep_window + sep, sep_window = optimization_config.sep, optimization_config.sep_window # unpack data e_dep, e_sep = dep_sep_data.dep.e, dep_sep_data.sep.e @@ -192,8 +192,8 @@ function fit_sg_wl(dep_sep_data::NamedTuple{(:dep, :sep)}, a_grid_wl_sg::StepRan # for each window lenght, calculate the survival fraction in the SEP for (i_aoe, wl) in enumerate(a_grid_wl_sg) - aoe_dep_i = aoe_dep[i_aoe, :][isfinite.(aoe_dep[i_aoe, :])] ./ mvalue(result_dep.m_calib) - e_dep_i = e_dep_calib[isfinite.(aoe_dep[i_aoe, :])] + aoe_dep_i = flatview(aoe_dep)[i_aoe, :][isfinite.(flatview(aoe_dep)[i_aoe, :])] ./ mvalue(result_dep.m_calib) + e_dep_i = e_dep_calib[isfinite.(flatview(aoe_dep)[i_aoe, :])] # prepare AoE max_aoe_dep_i = quantile(aoe_dep_i, optimization_config.max_aoe_quantile) + optimization_config.max_aoe_offset @@ -205,8 +205,8 @@ function fit_sg_wl(dep_sep_data::NamedTuple{(:dep, :sep)}, a_grid_wl_sg::StepRan try psd_cut = get_aoe_cut(aoe_dep_i, e_dep_i; window=dep_window, cut_search_interval=(min_aoe_dep_i, max_aoe_dep_i), uncertainty=uncertainty) - aoe_sep_i = aoe_sep[i_aoe, :][isfinite.(aoe_sep[i_aoe, :])] ./ result_dep.m_calib - e_sep_i = e_sep_calib[isfinite.(aoe_sep[i_aoe, :])] + aoe_sep_i = flatview(aoe_sep)[i_aoe, :][isfinite.(flatview(aoe_sep)[i_aoe, :])] ./ result_dep.m_calib + e_sep_i = e_sep_calib[isfinite.(flatview(aoe_sep)[i_aoe, :])] result_sep, _ = get_peak_surrival_fraction(aoe_sep_i, e_sep_i, sep, sep_window, psd_cut.lowcut; uncertainty=uncertainty, low_e_tail=false) From daa515aec015a0cbe543e8cfe4c91e3d373fe6dc Mon Sep 17 00:00:00 2001 From: Florian Henkes Date: Fri, 10 May 2024 17:38:37 +0200 Subject: [PATCH 60/68] remove static qc cut generation --> Only fits to blmean, sigma slope --- src/qc.jl | 69 ++++++++++-------------------------------------- src/singlefit.jl | 10 ++++--- 2 files changed, 20 insertions(+), 59 deletions(-) diff --git a/src/qc.jl b/src/qc.jl index 5c573f2f..60310ec4 100644 --- a/src/qc.jl +++ b/src/qc.jl @@ -1,72 +1,31 @@ - - -""" - qc_sg_optimization(dsp_dep, dsp_sep, optimization_config) - -Perform simple QC cuts on the DEP and SEP data and return the data for the optimization of the SG window length. -""" -function qc_sg_optimization(dsp_dep::NamedTuple{(:aoe, :e, :blmean, :blslope, :t50)}, dsp_sep::NamedTuple{(:aoe, :e, :blmean, :blslope, :t50)}, optimization_config::PropDict) - ### DEP - # Load DEP data and prepare Pile-up cut - blslope_dep, t50_dep = dsp_dep.blslope[isfinite.(dsp_dep.e)], dsp_dep.t50[isfinite.(dsp_dep.e)] - aoe_dep, e_dep = dsp_dep.aoe[:, isfinite.(dsp_dep.e)], dsp_dep.e[isfinite.(dsp_dep.e)] - # get half truncated centered cut on blslope for pile-up rejection - result_dep_slope_cut, report_dep_slope_cut = get_centered_gaussian_window_cut(blslope_dep, -0.1u"ns^-1", 0.1u"ns^-1", optimization_config.cuts.dep.blslope_sigma, ; n_bins_cut=optimization_config.cuts.dep.nbins_blslope_cut, relative_cut=optimization_config.cuts.dep.rel_cut_blslope_cut) - # Cut on blslope, energy and t0 for simple QC - qc_cut_dep = blslope_dep .> result_dep_slope_cut.low_cut .&& blslope_dep .< result_dep_slope_cut.high_cut .&& e_dep .> optimization_config.cuts.dep.min_e .&& quantile(e_dep, first(optimization_config.cuts.dep.e_quantile)) .< e_dep .< quantile(e_dep, last(optimization_config.cuts.dep.e_quantile)) .&& first(optimization_config.cuts.dep.t50) .< t50_dep .< last(optimization_config.cuts.dep.t50) - aoe_dep, e_dep = aoe_dep[:, qc_cut_dep], e_dep[qc_cut_dep] - - ### SEP - # Load SEP data and prepare Pile-up cut - blslope_sep, t50_sep = dsp_sep.blslope[isfinite.(dsp_sep.e)], dsp_sep.t50[isfinite.(dsp_sep.e)] - aoe_sep, e_sep = dsp_sep.aoe[:, isfinite.(dsp_sep.e)], dsp_sep.e[isfinite.(dsp_sep.e)] - - # get half truncated centered cut on blslope for pile-up rejection - result_sep_slope_cut, report_sep_slope_cut = get_centered_gaussian_window_cut(blslope_sep, -0.1u"ns^-1", 0.1u"ns^-1", optimization_config.cuts.sep.blslope_sigma, ; n_bins_cut=optimization_config.cuts.sep.nbins_blslope_cut, relative_cut=optimization_config.cuts.sep.rel_cut_blslope_cut) - - # Cut on blslope, energy and t0 for simple QC - qc_cut_sep = blslope_sep .> result_sep_slope_cut.low_cut .&& blslope_sep .< result_sep_slope_cut.high_cut .&& e_sep .> optimization_config.cuts.sep.min_e .&& quantile(e_sep, first(optimization_config.cuts.sep.e_quantile)) .< e_sep .< quantile(e_sep, last(optimization_config.cuts.sep.e_quantile)) .&& first(optimization_config.cuts.sep.t50) .< t50_sep .< last(optimization_config.cuts.sep.t50) - aoe_sep, e_sep = aoe_sep[:, qc_cut_sep], e_sep[qc_cut_sep] - - return (dep=(aoe=aoe_dep, e=e_dep), sep=(aoe=aoe_sep, e=e_sep)) -end -export qc_sg_optimization - - """ - qc_cal_energy(data, qc_config) + baseline_qc(data, qc_config) -Perform simple QC cuts on the data and return the data for energy calibration. +Perform simple Gaussian fits on the baseline disitrbutions for a given data set. +# Returns +- `result`: Namedtuple containing the cut and fit results for each baseline variable +- `report`: Namedtuple containing plotable objects. """ -function qc_cal_energy(data::Q, qc_config::PropDict) where Q<:Table +function baseline_qc(data::Q, qc_config::PropDict) where Q<:Table # get bl mean cut - result_blmean, _ = get_centered_gaussian_window_cut(data.blmean, qc_config.blmean.min, qc_config.blmean.max, qc_config.blmean.sigma, ; n_bins_cut=convert(Int64, round(length(data) * qc_config.blmean.n_bins_fraction)), relative_cut=qc_config.blmean.relative_cut, fixed_center=false, left=true) + result_blmean, report_blmean = get_centered_gaussian_window_cut(data.blmean, qc_config.blmean.min, qc_config.blmean.max, qc_config.blmean.sigma, ; n_bins_cut=convert(Int64, round(length(data) * qc_config.blmean.n_bins_fraction)), relative_cut=qc_config.blmean.relative_cut, fixed_center=false, left=true) blmean_qc = result_blmean.low_cut .< data.blmean .< result_blmean.high_cut @debug format("Baseline Mean cut surrival fraction {:.2f}%", count(blmean_qc) / length(data) * 100) # get bl slope cut - result_blslope, _ = get_centered_gaussian_window_cut(data.blslope, qc_config.blslope.min, qc_config.blslope.max, qc_config.blslope.sigma, ; n_bins_cut=convert(Int64, round(length(data) * qc_config.blslope.n_bins_fraction)), relative_cut=qc_config.blslope.relative_cut, fixed_center=true, left=false, center=zero(data.blslope[1])) + result_blslope, report_blslope = get_centered_gaussian_window_cut(data.blslope, qc_config.blslope.min, qc_config.blslope.max, qc_config.blslope.sigma, ; n_bins_cut=convert(Int64, round(length(data) * qc_config.blslope.n_bins_fraction)), relative_cut=qc_config.blslope.relative_cut, fixed_center=true, left=false, center=zero(data.blslope[1])) blslope_qc = result_blslope.low_cut .< data.blslope .< result_blslope.high_cut @debug format("Baseline Slope cut surrival fraction {:.2f}%", count(blslope_qc) / length(data) * 100) # get blsigma cut - result_blsigma, _ = get_centered_gaussian_window_cut(data.blsigma, qc_config.blsigma.min, qc_config.blsigma.max, qc_config.blsigma.sigma, ; n_bins_cut=convert(Int64, round(length(data) * qc_config.blsigma.n_bins_fraction)), relative_cut=qc_config.blsigma.relative_cut, fixed_center=false, left=true) + result_blsigma, report_blsigma = get_centered_gaussian_window_cut(data.blsigma, qc_config.blsigma.min, qc_config.blsigma.max, qc_config.blsigma.sigma, ; n_bins_cut=convert(Int64, round(length(data) * qc_config.blsigma.n_bins_fraction)), relative_cut=qc_config.blsigma.relative_cut, fixed_center=false, left=true) blsigma_qc = result_blsigma.low_cut .< data.blsigma .< result_blsigma.high_cut @debug format("Baseline Sigma cut surrival fraction {:.2f}%", count(blsigma_qc) / length(data) * 100) - # get t0 cut - t0_qc = qc_config.t0.min .< data.t0 .< qc_config.t0.max - @debug format("t0 cut surrival fraction {:.2f}%", count(t0_qc) / length(data) * 100) - # get intrace pile-up cut - inTrace_qc = .!(data.inTrace_intersect .> data.t0 .+ 2 .* data.drift_time .&& data.inTrace_n .> 1) - @debug format("Intrace pile-up cut surrival fraction {:.2f}%", count(inTrace_qc) / length(data) * 100) - # get energy cut - energy_qc = qc_config.e_trap.min .< data.e_trap .&& isfinite.(data.e_trap) .&& isfinite.(data.e_zac) .&& isfinite.(data.e_cusp) - @debug format("Energy cut surrival fraction {:.2f}%", count(energy_qc) / length(data) * 100) - + + result = (blmean = result_blmean, blslope = result_blslope, blsigma = result_blsigma) + report = (blmean = report_blmean, blslope = report_blslope, blsigma = report_blsigma) # combine all cuts - qc_tab = TypedTables.Table(blmean = blmean_qc, blslope = blslope_qc, blsigma = blsigma_qc, t0 = t0_qc, inTrace = inTrace_qc, energy = energy_qc, qc = blmean_qc .&& blslope_qc .&& blsigma_qc .&& t0_qc .&& inTrace_qc .&& energy_qc) - @debug format("Total QC cut surrival fraction {:.2f}%", count(qc) / length(data) * 100) - return qc_tab, (blmean = result_blmean, blslope = result_blslope, blsigma = result_blsigma) + return result, report end -export qc_cal_energy +export baseline_qc function get_pulser_identified_idx(data::Table, key::Symbol, σ_search::Vector{Int}, pulser_config::PropDict; min_identified::Int=10) diff --git a/src/singlefit.jl b/src/singlefit.jl index 5e9a90cb..6ce80311 100644 --- a/src/singlefit.jl +++ b/src/singlefit.jl @@ -24,8 +24,9 @@ function fit_single_trunc_gauss(x::Vector{T}, cuts::NamedTuple{(:low, :high, :ma pseudo_prior = NamedTupleDist( μ = Uniform(ps.peak_pos-20, ps.peak_pos+20), σ = weibull_from_mx(ps.peak_sigma, 3*ps.peak_sigma), - n = Uniform(ps.peak_counts-100, ps.peak_counts+100) + n = weibull_from_mx(ps.peak_counts, 3*ps.peak_counts) ) + # create fit model f_trafo = BAT.DistributionTransform(Normal, pseudo_prior) # f_trafo = LegendSpecFitsBATExt.get_distribution_transform(Normal, pseudo_prior) @@ -94,8 +95,9 @@ function fit_half_centered_trunc_gauss(x::Vector{T}, μ::T, cuts::NamedTuple{(:l pseudo_prior = NamedTupleDist( μ = ConstValueDist(μ), σ = weibull_from_mx(ps.peak_sigma, 3*ps.peak_sigma), - n = Uniform(ps.peak_counts-100, ps.peak_counts+100) + n = weibull_from_mx(ps.peak_counts, 3*ps.peak_counts) ) + # create fit model f_trafo = BAT.DistributionTransform(Normal, pseudo_prior) # f_trafo = LegendSpecFitsBATExt.get_distribution_transform(Normal, pseudo_prior) @@ -144,7 +146,7 @@ export fit_half_centered_trunc_gauss """ fit_half_centered_trunc_gauss(x::Array, cuts::NamedTuple{(:low, :high, :max), Tuple{Float64, Float64, Float64}}) Fit a single truncated Gaussian to the data `x` between `cut.low` and `cut.high`. The peak center is fixed at `μ` and the peak is cut in half either in the left or right half. -# Returns `report` and `result`` with: +# Returns `report` and `result` with: * `f_fit`: fitted function * `μ`: mean of the Gaussian * `μ_err`: error of the mean @@ -166,7 +168,7 @@ function fit_half_trunc_gauss(x::Vector{T}, cuts::NamedTuple{(:low, :high, :max) pseudo_prior = NamedTupleDist( μ = Uniform(ps.peak_pos-2*ps.peak_sigma, ps.peak_pos+2*ps.peak_sigma), σ = weibull_from_mx(ps.peak_sigma, 3*ps.peak_sigma), - n = Uniform(ps.peak_counts-100, ps.peak_counts+100) + n = weibull_from_mx(ps.peak_counts, 3*ps.peak_counts) ) # create fit model f_trafo = BAT.DistributionTransform(Normal, pseudo_prior) From 0d799cab58e3a47d1719c465feb8afa8e0ae5f28 Mon Sep 17 00:00:00 2001 From: Florian Henkes Date: Fri, 10 May 2024 17:38:54 +0200 Subject: [PATCH 61/68] Bug fix baseline fit recipes --- ext/LegendSpecFitsRecipesBaseExt.jl | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/ext/LegendSpecFitsRecipesBaseExt.jl b/ext/LegendSpecFitsRecipesBaseExt.jl index 0faba9ee..feea2587 100644 --- a/ext/LegendSpecFitsRecipesBaseExt.jl +++ b/ext/LegendSpecFitsRecipesBaseExt.jl @@ -478,10 +478,14 @@ end end @recipe function f(report_window_cut::NamedTuple{(:h, :f_fit, :x_fit, :low_cut, :high_cut, :low_cut_fit, :high_cut_fit, :center, :σ)}) - xlims := (value(ustrip(report_window_cut.center - 5*report_window_cut.σ)), value(ustrip(report_window_cut.center + 5*report_window_cut.σ))) + xlims := (value(ustrip(report_window_cut.center - 7*report_window_cut.σ)), value(ustrip(report_window_cut.center + 7*report_window_cut.σ))) + ylabel := "Density" + framestyle := :box @series begin - report_window_cut.h + label := "Data" + report_window_cut.h end + ylims --> (0, 1.1*maximum(report_window_cut.h.weights)) @series begin seriestype := :line label := "Best Fit" From bdc638bfa0961fb2712fb8a91ccddd6d9f06ef70 Mon Sep 17 00:00:00 2001 From: Florian Henkes Date: Sat, 11 May 2024 16:44:54 +0200 Subject: [PATCH 62/68] New simple calibration routines exclusivley on peakfinder --- src/simple_calibration.jl | 19 +++++++++++++------ 1 file changed, 13 insertions(+), 6 deletions(-) diff --git a/src/simple_calibration.jl b/src/simple_calibration.jl index ca63a67d..c8c4aec9 100644 --- a/src/simple_calibration.jl +++ b/src/simple_calibration.jl @@ -25,7 +25,7 @@ function simple_calibration(e_uncal::Vector{<:Real}, th228_lines::Vector{<:Unitf # remove :calib_type from kwargs kwargs = pairs(NamedTuple(filter(k -> !(:calib_type in k), kwargs))) if calib_type == :th228 - @info "Use simple calibration for Th228 lines" + @debug "Use simple calibration for Th228 lines" return simple_calibration_th228(e_uncal, th228_lines, window_sizes,; kwargs...) else error("Calibration type not supported") @@ -35,13 +35,20 @@ simple_calibration(e_uncal::Vector{<:Real}, th228_lines::Vector{<:Unitful.Energy 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") + # 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 - h_uncal = fit(Histogram, e_uncal, nbins=n_bins) - # search all possible peak candidates - _, peakpos = RadiationSpectra.peakfinder(h_uncal, σ=5.0, backgroundRemove=true, threshold=10) - # the FEP ist the last peak in the list + h_uncal = fit(Histogram, e_uncal, 0:bin_width:maximum(e_uncal)) fep_guess = if isnan(quantile_perc) - sort(peakpos)[end] + # expect FEP in the last 10% of the data + min_e_fep = quantile(e_uncal, 0.9) + h_fepsearch = fit(Histogram, e_uncal, min_e_fep:bin_width:maximum(e_uncal)) + # search all possible peak candidates + h_decon, peakpos = RadiationSpectra.peakfinder(h_fepsearch, σ=5.0, backgroundRemove=true, threshold=10) + # the FEP is 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)] else quantile(e_uncal, quantile_perc) end From 3c59eb03d488ff224590ebbecf6923a0514630e2 Mon Sep 17 00:00:00 2001 From: Florian Henkes Date: Sat, 11 May 2024 17:09:53 +0200 Subject: [PATCH 63/68] Remove info message --- src/specfit.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/specfit.jl b/src/specfit.jl index 5137c701..a8202bfa 100644 --- a/src/specfit.jl +++ b/src/specfit.jl @@ -108,7 +108,6 @@ function fit_peaks(peakhists::Array, peakstats::StructArray, th228_lines::Vector # remove :calib_type from kwargs kwargs = pairs(NamedTuple(filter(k -> !(:calib_type in k), kwargs))) if calib_type == :th228 - @info "Fit peaks for Th228 lines" return fit_peaks_th228(peakhists, peakstats, th228_lines,; kwargs...) else error("Calibration type not supported") From f76f6eedf02b93aa67d8c26355b14be3a0610ba1 Mon Sep 17 00:00:00 2001 From: Felix Hagemann Date: Sat, 11 May 2024 17:06:27 +0200 Subject: [PATCH 64/68] Bugfix in `get_peak_survival_fraction` --- src/aoe_calibration.jl | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/aoe_calibration.jl b/src/aoe_calibration.jl index ecc1f528..2a598672 100644 --- a/src/aoe_calibration.jl +++ b/src/aoe_calibration.jl @@ -225,10 +225,9 @@ function get_peak_surrival_fraction(aoe::Vector{<:Unitful.RealOrRealQuantity}, e if !isnan(sigma_high_sided) # TODO: decide how to deal with the high sided cut! e = e[aoe .< sigma_high_sided] - else - e_survived = e[aoe_cut .<= aoe] - e_cut = e[aoe_cut .> aoe] end + e_survived = e[aoe_cut .<= aoe] + e_cut = e[aoe_cut .> aoe] # estimate bin width bin_width = get_friedman_diaconis_bin_width(e[e .> peak - bin_width_window .&& e .< peak + bin_width_window]) From d6889b9e40d6aee782ec4aafde3700c290e1a731 Mon Sep 17 00:00:00 2001 From: Florian Henkes Date: Wed, 15 May 2024 00:45:12 +0200 Subject: [PATCH 65/68] Delete legacy functions --- src/aoe_calibration.jl | 33 +-------------------------------- 1 file changed, 1 insertion(+), 32 deletions(-) diff --git a/src/aoe_calibration.jl b/src/aoe_calibration.jl index 2a598672..d162ecf1 100644 --- a/src/aoe_calibration.jl +++ b/src/aoe_calibration.jl @@ -140,37 +140,6 @@ The algorhithm utilizes a root search algorithm to find the cut value with a rel - `n0`: Number of counts before the cut - `nsf`: Number of counts after the cut """ -#= -function get_aoe_cut(aoe::Vector{<:Unitful.RealOrRealQuantity}, e::Vector{<:T},; dep::T=1592.53u"keV", window::Vector{<:T}=[12.0, 10.0]u"keV", dep_sf::Float64=0.9, cut_search_interval::Tuple{<:Unitful.RealOrRealQuantity, <:Unitful.RealOrRealQuantity}=(-25.0*unit(first(aoe)), 1.0*unit(first(aoe))), rtol::Float64=0.001, bin_width_window::T=3.0u"keV", fixed_position::Bool=true, sigma_high_sided::Float64=NaN, uncertainty::Bool=true, maxiters::Int=200) where T<:Unitful.Energy{<:Real} - # check if a high sided AoE cut should be applied before the AoE cut is generated - if !isnan(sigma_high_sided) - e = e[aoe .< sigma_high_sided] - aoe = aoe[aoe .< sigma_high_sided] - end - # cut window around peak - aoe = aoe[dep-first(window) .< e .< dep+last(window)] - e = e[dep-first(window) .< e .< dep+last(window)] - # estimate bin width - bin_width = get_friedman_diaconis_bin_width(e[dep - bin_width_window .< e .< dep + bin_width_window]) - # create histogram - dephist = fit(Histogram, ustrip.(e), ustrip(dep-first(window):bin_width:dep+last(window))) - # get peakstats - depstats = estimate_single_peak_stats(dephist) - # fit before cut - result_before, _ = fit_single_peak_th228(dephist, depstats,; uncertainty=uncertainty, fixed_position=fixed_position, low_e_tail=false) - # get n0 before cut - nsf = result_before.n * dep_sf - # get aoe cut - n_surrival_dep_f = cut -> get_n_after_aoe_cut(cut, aoe, e, dep, window, bin_width, mvalue(result_before), depstats; uncertainty=false, fixed_position=fixed_position) - nsf - aoe_cut = find_zero(n_surrival_dep_f, cut_search_interval, Bisection(), rtol=rtol, maxiters=maxiters) - # return n_surrival_dep_f.(0.25:0.001:0.5) - # get nsf after cut - nsf = get_n_after_aoe_cut(aoe_cut, aoe, e, dep, window, bin_width, mvalue(result_before), depstats; uncertainty=uncertainty, fixed_position=fixed_position) - return (lowcut = measurement(aoe_cut, aoe_cut * rtol), highcut = sigma_high_sided, n0 = result_before.n, nsf = nsf, sf = nsf / result_before.n * 100*u"percent") -end -export get_aoe_cut -=# - function get_aoe_cut(aoe::Vector{<:Unitful.RealOrRealQuantity}, e::Vector{<:T},; dep::T=1592.53u"keV", window::Vector{<:T}=[12.0, 10.0]u"keV", dep_sf::Float64=0.9, cut_search_interval::Tuple{<:Unitful.RealOrRealQuantity, <:Unitful.RealOrRealQuantity}=(-25.0*unit(first(aoe)), 1.0*unit(first(aoe))), rtol::Float64=0.001, bin_width_window::T=3.0u"keV", fixed_position::Bool=true, sigma_high_sided::Float64=NaN, uncertainty::Bool=true, maxiters::Int=200) where T<:Unitful.Energy{<:Real} # check if a high sided AoE cut should be applied before the AoE cut is generated if !isnan(sigma_high_sided) @@ -225,6 +194,7 @@ function get_peak_surrival_fraction(aoe::Vector{<:Unitful.RealOrRealQuantity}, e if !isnan(sigma_high_sided) # TODO: decide how to deal with the high sided cut! e = e[aoe .< sigma_high_sided] + aoe = aoe[aoe .< sigma_high_sided] end e_survived = e[aoe_cut .<= aoe] e_cut = e[aoe_cut .> aoe] @@ -236,7 +206,6 @@ function get_peak_surrival_fraction(aoe::Vector{<:Unitful.RealOrRealQuantity}, e cut = fit(Histogram, ustrip(e_cut), ustrip(peak-first(window):bin_width:peak+last(window))) # fit peak and return number of signal counts result_after, report_after = fit_subpeaks_th228(survived, cut, result_before; uncertainty=uncertainty, low_e_tail=low_e_tail) - # result_after, report_after = fit_single_peak_th228(peakhist, peakstats,; uncertainty=uncertainty, low_e_tail=low_e_tail, pseudo_prior=pseudo_prior) # calculate surrival fraction sf = result_after.sf * 100u"percent" result = ( From 3227a18bc05be86d159b5a6dac298a400137bdc0 Mon Sep 17 00:00:00 2001 From: Florian Henkes Date: Wed, 15 May 2024 00:46:32 +0200 Subject: [PATCH 66/68] Allow fit_calibration to fit mu with units. Strips unit to not break recipes --- src/fit_calibration.jl | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/fit_calibration.jl b/src/fit_calibration.jl index dd28887b..563a5fb4 100644 --- a/src/fit_calibration.jl +++ b/src/fit_calibration.jl @@ -9,9 +9,13 @@ Fit the calibration lines with polynomial function of n_poly order * `gof`: godness of fit * `report`: """ -function fit_calibration(n_poly::Int, µ::AbstractVector{<:Union{Real,Measurement{<:Real}}}, peaks::AbstractVector{<:Quantity}; e_expression::Union{Symbol, String}="e") +function fit_calibration(n_poly::Int, µ::AbstractVector{<:Union{Unitful.RealOrRealQuantity,Measurement{<:Unitful.RealOrRealQuantity}}}, peaks::AbstractVector{<:Quantity}; e_expression::Union{Symbol, String}="e") @assert length(peaks) == length(μ) e_unit = u"keV" + if unit(first(μ)) != NoUnits + @warn "The unit of µ is not $(e_unit), it will be converted to $(e_unit) and stripped." + µ = ustrip.(e_unit, µ) + end @debug "Fit calibration curve with $(n_poly)-order polynominal function" result_fit, report = chi2fit(n_poly, µ, ustrip.(e_unit, peaks)) From f7a67c36e3687e6fa9bae26f25778d4d3bf62fa6 Mon Sep 17 00:00:00 2001 From: Florian Henkes Date: Wed, 15 May 2024 13:31:47 +0200 Subject: [PATCH 67/68] Fixed recipes to adapt to new values of absolute ADC --- ext/LegendSpecFitsRecipesBaseExt.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/ext/LegendSpecFitsRecipesBaseExt.jl b/ext/LegendSpecFitsRecipesBaseExt.jl index feea2587..b7a2d55c 100644 --- a/ext/LegendSpecFitsRecipesBaseExt.jl +++ b/ext/LegendSpecFitsRecipesBaseExt.jl @@ -669,8 +669,8 @@ end xlabel := "Energy (ADC)" legend := :bottomright framestyle := :box - xlims := (0, 21000) - xticks := (0:2000:22000) + xlims := (0, 168000) + xticks := (0:16000:176000) @series begin grid --> :all xerrscaling := xerrscaling From c84abb3496870c13422e675bfc9ba1605974e8cb Mon Sep 17 00:00:00 2001 From: Oliver Schulz Date: Fri, 24 May 2024 10:08:32 +0200 Subject: [PATCH 68/68] Bump LegendDataManagement version bound --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index a7111766..237a091c 100644 --- a/Project.toml +++ b/Project.toml @@ -58,7 +58,7 @@ IntervalSets = "0.7" InverseFunctions = "0.1" IrrationalConstants = "0.1, 0.2" LaTeXStrings = "1.3" -LegendDataManagement = "0.2.7" +LegendDataManagement = "0.2.7, 0.3" LinearAlgebra = "1" LinearRegression = "0.2" LsqFit = "0.14, 0.15"