Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Generic simple_calibration #108

Open
wants to merge 10 commits into
base: main
Choose a base branch
from
Open

Generic simple_calibration #108

Show file tree
Hide file tree
Changes from 3 commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
a0b628f
add more generic verison of simple calibration, e.g. for co60 data
LisaSchlueter Dec 12, 2024
fa21052
take most prominent peak not not largest energy
LisaSchlueter Dec 12, 2024
1639f38
add minimal peakshape without tail
LisaSchlueter Dec 14, 2024
888996b
update docstring and add fall-back solution for fwhm calculation
LisaSchlueter Dec 16, 2024
421e2c3
add e_unit as kwarg
LisaSchlueter Dec 16, 2024
692c1de
merge commong function for th228 and generic gamma
LisaSchlueter Jan 10, 2025
b8dda3f
clean up th228 remainders and docstrings
LisaSchlueter Jan 10, 2025
8e924d9
merge simple calibration plto recipes
LisaSchlueter Jan 10, 2025
13a85f9
add simple calibration test
LisaSchlueter Jan 11, 2025
4851690
bug fix tests
LisaSchlueter Jan 11, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
38 changes: 37 additions & 1 deletion ext/LegendSpecFitsRecipesBaseExt.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -490,6 +490,41 @@ end
end
end

@recipe function f(report::NamedTuple{(:h_calsimple, :h_uncal, :c, :peak_guess, :peakhists, :peakstats)}; cal=true)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't understand why this is a new recipe. Was this not here already before?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is almost a duplicate of the existing recipe. However, the old one has some stuff hardcoded, e.g. "FEP" in the plot label. I didn't want to mess with the existing one, since we're int he middle of the processing.

ylabel := "Counts"
legend := :topright
yscale := :log10
fill := false
if cal
h = LinearAlgebra.normalize(report.h_calsimple, mode = :density)
xlabel := "Energy (keV)"
xlims := (0, 3000)
xticks := (0:500:3000, ["$i" for i in 0:500:3000])
ylims := (0.2, maximum(h.weights)*1.1)
peak_guess = ustrip(report.c * report.peak_guess)
else
h = LinearAlgebra.normalize(report.h_uncal, mode = :density)
xlabel := "Energy (ADC)"
xlims := (0, 1.2*report.peak_guess)
xticks := (0:3000:1.2*report.peak_guess, ["$i" for i in 0:3000:1.2*report.peak_guess])
ylims := (0.2, maximum(h.weights)*1.1)
peak_guess = report.peak_guess
end
@series begin
seriestype := :stepbins
label := "Energy"
h
end
y_vline = 0.2:1:maximum(h.weights)*1.1
@series begin
seriestype := :line
label := "Peak guess"#: $(round(peak_guess, digits = 1))"
color := :red2
linewidth := 1.5
fill(peak_guess, length(y_vline)), y_vline
end
end

@recipe function f(report::NamedTuple{(:peakpos, :peakpos_cal, :h_uncal, :h_calsimple)}; cal=true)
legend := :topright
size := (1000, 600)
Expand Down Expand Up @@ -820,7 +855,8 @@ end
if plot_ribbon
ribbon := uncertainty.(report.f_fit.(0:1:1.2*value(maximum(report.x))))
end
0:1:1.2*value(maximum(report.x)), value.(report.f_fit.(0:1:1.2*value(maximum(report.x))))
xstep = value(maximum(report.x))/100
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Does this also work for the 228-Th data

0:xstep:1.2*value(maximum(report.x)), value.(report.f_fit.(0:xstep:1.2*value(maximum(report.x))))
end
@series begin
seriestype := :scatter
Expand Down
2 changes: 2 additions & 0 deletions src/pseudo_prior.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -33,6 +33,8 @@ function get_standard_pseudo_prior(h::Histogram, ps::NamedTuple{(:peak_pos, :pea
NamedTupleDist(; μ, σ, n, skew_fraction, skew_width, background, step_amplitude, background_slope)
elseif fit_func == :gamma_bckExp
NamedTupleDist(; μ, σ, n, skew_fraction, skew_width, background, step_amplitude, background_exp)
elseif fit_func == :gamma_minimal
NamedTupleDist(; μ, σ, n, background)
else
throw(ArgumentError("Unknown fit function: $fit_func"))
end
Expand Down
63 changes: 59 additions & 4 deletions src/simple_calibration.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,20 +18,22 @@ Returns
function simple_calibration end
export simple_calibration

function simple_calibration(e_uncal::Vector{<:Real}, th228_lines::Vector{<:Unitful.Energy{<:Real}}, window_sizes::Vector{<:Tuple{Unitful.Energy{<:Real}, Unitful.Energy{<:Real}}},; kwargs...)
function simple_calibration(e_uncal::Vector{<:Real}, gamma_lines::Vector{<:Unitful.Energy{<:Real}}, window_sizes::Vector{<:Tuple{Unitful.Energy{<:Real}, Unitful.Energy{<:Real}}},; kwargs...)
# remove calib type from kwargs
@assert haskey(kwargs, :calib_type) "Calibration type not specified"
calib_type = kwargs[:calib_type]
# remove :calib_type from kwargs
kwargs = pairs(NamedTuple(filter(k -> !(:calib_type in k), kwargs)))
if calib_type == :th228
@debug "Use simple calibration for Th228 lines"
return simple_calibration_th228(e_uncal, th228_lines, window_sizes,; kwargs...)
return simple_calibration_th228(e_uncal, gamma_lines, window_sizes,; kwargs...)
elseif calib_type == :gamma
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why not call this :co60 if it is for the 60-Co data? My original idea was here to have something for different isotopes since you always look for different features in the data.
But I am actually open for other suggestions.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This works also for any calibration data. I think we should have only 1 generic simple_calibration function in the future and only special one if we can't make it work with the generic one

return simple_calibration_gamma(e_uncal, gamma_lines, window_sizes,; kwargs...)
else
error("Calibration type not supported")
end
end
simple_calibration(e_uncal::Vector{<:Real}, th228_lines::Vector{<:Unitful.Energy{<:Real}}, left_window_sizes::Vector{<:Unitful.Energy{<:Real}}, right_window_sizes::Vector{<:Unitful.Energy{<:Real}}; kwargs...) = simple_calibration(e_uncal, th228_lines, [(l,r) for (l,r) in zip(left_window_sizes, right_window_sizes)],; kwargs...)
simple_calibration(e_uncal::Vector{<:Real}, gamma_lines::Vector{<:Unitful.Energy{<:Real}}, left_window_sizes::Vector{<:Unitful.Energy{<:Real}}, right_window_sizes::Vector{<:Unitful.Energy{<:Real}}; kwargs...) = simple_calibration(e_uncal, gamma_lines, [(l,r) for (l,r) in zip(left_window_sizes, right_window_sizes)],; kwargs...)


function simple_calibration_th228(e_uncal::Vector{<:Real}, th228_lines::Vector{<:Unitful.Energy{<:Real}}, window_sizes::Vector{<:Tuple{Unitful.Energy{<:Real}, Unitful.Energy{<:Real}}},; n_bins::Int=15000, quantile_perc::Float64=NaN, binning_peak_window::Unitful.Energy{<:Real}=10.0u"keV")
Expand Down Expand Up @@ -81,7 +83,60 @@ function simple_calibration_th228(e_uncal::Vector{<:Real}, th228_lines::Vector{<
return result, report
end


"""
simple_calibration_gamma(e_uncal::Array, gamma_lines::Array, window_size::Float64=25.0, n_bins::Int=15000, quantile_perc::Float64=NaN, binning_peak_window::Float64=10.0, peak_quantile::ClosedInterval{<:Real} = 0.5..1.0)
Perform a simple calibration for the uncalibrated energy array `e_uncal`
* Find peaks within peak_quantile::ClosedInterval{<:Real} = 0.5..1.0 of the data
* Estimate the calibration constant `c` by dividing the maximum gamma line energy by the peak guess
* Return peakstats and peakhists of gamma lines for further calibration
"""
fhagemann marked this conversation as resolved.
Show resolved Hide resolved
function simple_calibration_gamma(e_uncal::Vector{<:Real}, gamma_lines::Vector{<:Unitful.Energy{<:Real}}, window_sizes::Vector{<:Tuple{Unitful.Energy{<:Real}, Unitful.Energy{<:Real}}},; n_bins::Int=15000, quantile_perc::Float64=NaN, binning_peak_window::Unitful.Energy{<:Real}=10.0u"keV", peak_quantile::ClosedInterval{<:Real} = 0.5..1.0)
e_min = quantile(e_uncal, leftendpoint(peak_quantile))
e_max = quantile(e_uncal, rightendpoint(peak_quantile))
# initial binning
bin_width = get_friedman_diaconis_bin_width(filter(in(e_min..e_max), e_uncal))
# create initial peak search histogram
h_uncal = fit(Histogram, e_uncal, 0:bin_width:maximum(e_uncal))
peak_guess, peak_idx = if isnan(quantile_perc)
h_peaksearch = fit(Histogram, e_uncal, e_min:bin_width:e_max)
# search all possible peak candidates
h_decon, peakpos = RadiationSpectra.peakfinder(h_peaksearch, σ=2.0, backgroundRemove=true, threshold=10)
# find the most prominent peak in the deconvoluted histogram
peakpos_idxs = StatsBase.binindex.(Ref(h_decon), peakpos)
cts_peakpos = h_decon.weights[peakpos_idxs]
peakpos[argmax(cts_peakpos)], argmax(cts_peakpos)
else
quantile(e_uncal, quantile_perc), length(gamma_lines)
end
@info "Identified most prominent peak (peak $(peak_idx), $(gamma_lines[peak_idx])). Peak guess: $peak_guess"
# get calibration constant for simple calibration
c = gamma_lines[peak_idx] / peak_guess
e_simple = e_uncal .* c
e_unit = u"keV"
fhagemann marked this conversation as resolved.
Show resolved Hide resolved
# get peakhists and peakstats
peakhists, peakstats, h_calsimple, bin_widths = get_peakhists_th228(e_simple, gamma_lines, window_sizes; binning_peak_window=binning_peak_window, e_unit=e_unit)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

get_peakhists_th228 also works for generic simple_calibration_gamma?
Or would that also need a get_peakhists_gamma?

Copy link
Contributor Author

@LisaSchlueter LisaSchlueter Dec 16, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes it does.
In principle we could consider re-writing/re-naming all ..._th228 calibration functions to be more generic for other gamma spectra. Some don't even need any modification at all and work out of the box for other sources.
I didn't do that to keep the changes as minimal as possible.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It makes sense to rename it, but maybe this could go along then with a PR onto the dataflow where we do a search and replace for these functions.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I agree. I can take care of this beginning of next year.

In that context, I'd suggest to change the julia energy-calibration metadata as well:

  1. Add a field "source" in the metadata, which is "th228" for LEGEND-200, but can also be something else. Like that the dataflow can pick the correct gamma_lines and gamma_names from the energy calibration config.

OR

  1. Rename "th228_lines" etc in the metadata to "gamma_lines".

I prefer option 1, because this allows you to switch easily back and fourth between different calibration sources.


result = (
h_calsimple = h_calsimple,
h_uncal = h_uncal,
c = c,
unit = e_unit,
bin_width = median(bin_widths),
peak_guess = peak_guess,
peakbinwidths = bin_widths,
peakhists = peakhists,
peakstats = peakstats
)
report = (
h_calsimple = result.h_calsimple,
h_uncal = result.h_uncal,
c = result.c,
peak_guess = result.peak_guess,
peakhists = result.peakhists,
peakstats = result.peakstats
)
return result, report
end
"""
get_peakhists_th228(e::Array, th228_lines::Array, window_sizes::Array, e_unit::String="keV", proxy_binning_peak::Float64=2103.5, proxy_binning_peak_window::Float64=10.0)

Expand Down
55 changes: 31 additions & 24 deletions src/specfit.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -194,7 +194,11 @@ export fit_single_peak_th228
calculate centroid of gamma peak from fit parameters
"""
function peak_centroid(v::NamedTuple)
centroid = v.μ - v.skew_fraction * (v.µ * v.skew_width)
if haskey(v, :skew_fraction)
centroid = v.μ - v.skew_fraction * (v.µ * v.skew_width)
else
centroid = v.μ
end
if haskey(v, :skew_fraction_highE)
centroid += v.skew_fraction_highE * (v.µ * v.skew_width_highE)
end
Expand All @@ -208,31 +212,34 @@ Get the FWHM of a peak from the fit parameters.
# Returns
* `fwhm`: the FWHM of the peak
"""
function estimate_fwhm(v::NamedTuple)

f_sigWithTail = Base.Fix2(get_th228_fit_functions().gamma_sigWithTail,v)
try
e_low, e_high = v.skew_fraction <= 0.5 ? (v.μ - v.σ, v.μ + v.σ) : (v.μ * (1 - v.skew_width), v.μ * (1 + v.skew_width))

max_sig = -Inf
for e in e_low:0.001:e_high
fe = f_sigWithTail(e)
if fe > max_sig
max_sig = fe
else
# if the maximum is reached,
# no need to further continue
break
function estimate_fwhm(v::NamedTuple)
if haskey(v, :skew_fraction)
f_sigWithTail = Base.Fix2(get_th228_fit_functions().gamma_sigWithTail,v)
try
e_low, e_high = v.skew_fraction <= 0.5 ? (v.μ - v.σ, v.μ + v.σ) : (v.μ * (1 - v.skew_width), v.μ * (1 + v.skew_width))

max_sig = -Inf
for e in e_low:0.001:e_high
fe = f_sigWithTail(e)
if fe > max_sig
max_sig = fe
else
# if the maximum is reached,
# no need to further continue
break
end
end
half_max_sig = max_sig/2

tmp = x -> f_sigWithTail(x) - half_max_sig
roots_low = find_zero(tmp, e_low, maxiter=100)
roots_high = find_zero(tmp, e_high, maxiter=100)
return roots_high - roots_low
catch
return NaN
end
half_max_sig = max_sig/2

tmp = x -> f_sigWithTail(x) - half_max_sig
roots_low = find_zero(tmp, e_low, maxiter=100)
roots_high = find_zero(tmp, e_high, maxiter=100)
return roots_high - roots_low
catch
return NaN
else
return 2 * sqrt(2 * log(2)) * v.σ
end
fhagemann marked this conversation as resolved.
Show resolved Hide resolved
end
"""
Expand Down
6 changes: 6 additions & 0 deletions src/specfit_functions.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,7 @@ This function defines the gamma peakshape fit functions used in the calibration
* gamma_bckExp: default gamma peakshape + exponential background
* gamma_bckFlat: default gamma peakshape - step background (only flat component!)
* gamma_tails_bckFlat: default gamma peakshape + high-energy tail - step background (only flat component!)
* gamma_minimal: only Gaussian signal and flat background
"""
function get_th228_fit_functions(; background_center::Union{Real,Nothing} = nothing)
merge(
Expand All @@ -16,6 +17,7 @@ function get_th228_fit_functions(; background_center::Union{Real,Nothing} = noth
gamma_sigWithTail = (x, v) -> signal_peakshape(x, v.μ, v.σ, v.n, v.skew_fraction) + lowEtail_peakshape(x, v.μ, v.σ, v.n, v.skew_fraction, v.skew_width),
gamma_bckFlat = (x, v) -> gamma_peakshape(x, v.μ, v.σ, v.n, 0.0, v.skew_fraction, v.skew_width, v.background),
gamma_tails_bckFlat = (x, v) -> gamma_peakshape(x, v.μ, v.σ, v.n, 0.0, v.skew_fraction, v.skew_width , v.background; skew_fraction_highE = v.skew_fraction_highE, skew_width_highE = v.skew_width_highE),
gamma_minimal = (x, v) -> gamma_peakshape(x, v.μ, v.σ, v.n, 0.0, 0.0, 0.0, v.background),
),
if isnothing(background_center)
(gamma_bckSlope = (x, v) -> gamma_peakshape(x, v.μ, v.σ, v.n, v.step_amplitude, v.skew_fraction, v.skew_width, v.background; background_slope = v.background_slope, background_center = v.μ),
Expand Down Expand Up @@ -61,6 +63,10 @@ function peakshape_components(fit_func::Symbol; background_center::Union{Real,No
f_lowEtail = (x, v) -> lowEtail_peakshape(x, v.μ, v.σ, v.n, v.skew_fraction, v.skew_width),
f_highEtail = (x, v) -> highEtail_peakshape(x, v.μ, v.σ, v.n, v.skew_fraction_highE, v.skew_width_highE),
f_bck = (x, v) -> background_peakshape(x, v.μ, v.σ, 0.0, v.background))
elseif fit_func == :gamma_minimal
funcs = (f_sig = (x, v) -> signal_peakshape(x, v.μ, v.σ, v.n, 0.0),
f_lowEtail = (x, v) -> lowEtail_peakshape(x, v.μ, v.σ, v.n, 0.0, 0.0),
f_bck = (x, v) -> background_peakshape(x, v.μ, v.σ, 0.0, v.background))
end
labels = (f_sig = "Signal", f_lowEtail = "Low-energy tail", f_bck = "Background", f_highEtail = "High-energy tail")
colors = (f_sig = :orangered1, f_lowEtail = :orange, f_bck = :dodgerblue2, f_highEtail = :forestgreen)
Expand Down
Loading