Generic simple_calibration #108

LisaSchlueter · 2024-12-16T22:21:34Z

add a more generic simple_calibration function that can be accessed with keyword calib_type == :gamma which is useful for example for Co-60 data
add a minimalistic gamma peak-shape, that contains only Gaussian signal and flat background to fit low-statistics data

fhagemann · 2024-12-16T22:41:38Z

src/simple_calibration.jl

+    e_simple = e_uncal .* c
+    e_unit = u"keV"
+    # 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)


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

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.

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.

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:

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

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.

src/simple_calibration.jl

src/specfit.jl

LisaSchlueter · 2024-12-16T23:20:44Z

Is this for Co60 or Co56?

I wrote and tested this for my local Co60 data, but it should work for any other calibration spectra, because the anticipated gamma energies are passed as an argument.

codecov · 2024-12-16T23:33:12Z

Codecov Report

Attention: Patch coverage is 28.04878% with 59 lines in your changes missing coverage. Please review.

Project coverage is 21.24%. Comparing base (7ff7aac) to head (421e2c3).

Files with missing lines	Patch %	Lines
ext/LegendSpecFitsRecipesBaseExt.jl	0.00%	31 Missing ⚠️
src/simple_calibration.jl	12.00%	22 Missing ⚠️
src/pseudo_prior.jl	0.00%	2 Missing ⚠️
src/specfit.jl	90.90%	2 Missing ⚠️
src/specfit_functions.jl	0.00%	2 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #108      +/-   ##
==========================================
- Coverage   21.33%   21.24%   -0.09%     
==========================================
  Files          36       36              
  Lines        3305     3365      +60     
==========================================
+ Hits          705      715      +10     
- Misses       2600     2650      +50

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

theHenks

In general, the PR looks good. Thanks @LisaSchlueter .
My comments are mainly about the naming.
Do we wanna stick to the isotope naming or start having more generic names?

theHenks · 2024-12-17T14:50:46Z

ext/LegendSpecFitsRecipesBaseExt.jl

@@ -490,6 +490,41 @@ end
    end
 end

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


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

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.

theHenks · 2024-12-17T14:51:25Z

ext/LegendSpecFitsRecipesBaseExt.jl

@@ -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


Does this also work for the 228-Th data

theHenks · 2024-12-17T14:52:44Z

src/simple_calibration.jl

    # 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


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.

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

theHenks · 2024-12-17T14:54:18Z

src/simple_calibration.jl



-function simple_calibration_th228(e_uncal::Vector{<:Real}, th228_lines::Vector{<:Unitful.Energy{<:Real}}, window_sizes::Vector{<:Tuple{Unitful.Energy{<:Real}, Unitful.Energy{<:Real}}},; n_bins::Int=15000, quantile_perc::Float64=NaN, binning_peak_window::Unitful.Energy{<:Real}=10.0u"keV")
+function simple_calibration_th228(e_uncal::Vector{<:Real}, th228_lines::Vector{<:Unitful.Energy{<:Real}}, window_sizes::Vector{<:Tuple{Unitful.Energy{<:Real}, Unitful.Energy{<:Real}}},; n_bins::Int=15000, quantile_perc::Float64=NaN, binning_peak_window::Unitful.Energy{<:Real}=10.0u"keV", e_unit::Unitful.EnergyUnits=u"keV")


I always had the unit and therefore the scale fixed in the code since for some of the values I assumed that everything from there on is in keV. Is there a specific reason why to have this as a kwarg? In my opinion, it should be fine to define the "energy scale" in the code itself.

theHenks · 2024-12-17T14:54:37Z

src/simple_calibration.jl

+* `binning_peak_window`: window size for the peak search histogram (default: 10 keV)
+* `peak_quantile`: quantile range for the peak search (default: 0.5..1.0)
+"""
+function simple_calibration_gamma(e_uncal::Vector{<:Real}, gamma_lines::Vector{<:Unitful.Energy{<:Real}}, window_sizes::Vector{<:Tuple{Unitful.Energy{<:Real}, Unitful.Energy{<:Real}}},; n_bins::Int=15000, quantile_perc::Float64=NaN, binning_peak_window::Unitful.Energy{<:Real}=10.0u"keV", peak_quantile::ClosedInterval{<:Real} = 0.5..1.0,  e_unit::Unitful.EnergyUnits=u"keV")


Maybe also rename to simpel_calibration_gamma

theHenks · 2024-12-17T14:54:54Z

src/simple_calibration.jl

+    e_simple = e_uncal .* c
+    e_unit = u"keV"
+    # 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)


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.

LisaSchlueter · 2024-12-17T16:21:06Z

Thanks for your input! I really appreciate it. I suggest the following before merging this PR:

LegendSpecFits:

rename everything th228 related into gamma
merge the two simple_calibration functions into one combined function. Same for plot recipe.
test everything on th228 and co60 data

julia-legend-dataflow:

adapt processor to new naming

metadata / legend-jldataflow-config:

add source field (see above)

fhagemann · 2024-12-20T10:28:49Z

That sounds like a great plan.
We can still define functions called simple_calibration_th228 that will just call simple_calibration_gamma, something along

simple_calibration_th228(args...; kwargs...) = simple_calibration_gamma(args...; source = :th228, kwargs...)

where the keyword argument is source, calib_type, ... or whatever.

LisaSchlueter added 3 commits December 12, 2024 15:23

add more generic verison of simple calibration, e.g. for co60 data

a0b628f

take most prominent peak not not largest energy

fa21052

add minimal peakshape without tail

1639f38

LisaSchlueter added the enhancement New feature or request label Dec 16, 2024

LisaSchlueter requested review from fhagemann and theHenks December 16, 2024 22:21

LisaSchlueter self-assigned this Dec 16, 2024

fhagemann reviewed Dec 16, 2024

View reviewed changes

LisaSchlueter added 2 commits December 16, 2024 15:22

update docstring and add fall-back solution for fwhm calculation

888996b

add e_unit as kwarg

421e2c3

theHenks requested changes Dec 17, 2024

View reviewed changes

LisaSchlueter added 5 commits January 10, 2025 14:04

merge commong function for th228 and generic gamma

692c1de

clean up th228 remainders and docstrings

b8dda3f

merge simple calibration plto recipes

8e924d9

add simple calibration test

13a85f9

bug fix tests

4851690

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Generic simple_calibration #108

Generic simple_calibration #108

LisaSchlueter commented Dec 16, 2024 •

edited

Loading

fhagemann Dec 16, 2024

LisaSchlueter Dec 16, 2024 •

edited

Loading

theHenks Dec 17, 2024

LisaSchlueter Dec 17, 2024

LisaSchlueter commented Dec 16, 2024

codecov bot commented Dec 16, 2024

theHenks left a comment

theHenks Dec 17, 2024

LisaSchlueter Dec 17, 2024

theHenks Dec 17, 2024

theHenks Dec 17, 2024

LisaSchlueter Dec 17, 2024

theHenks Dec 17, 2024

theHenks Dec 17, 2024

theHenks Dec 17, 2024

LisaSchlueter commented Dec 17, 2024

fhagemann commented Dec 20, 2024



		function simple_calibration_th228(e_uncal::Vector{<:Real}, th228_lines::Vector{<:Unitful.Energy{<:Real}}, window_sizes::Vector{<:Tuple{Unitful.Energy{<:Real}, Unitful.Energy{<:Real}}},; n_bins::Int=15000, quantile_perc::Float64=NaN, binning_peak_window::Unitful.Energy{<:Real}=10.0u"keV")
		function simple_calibration_th228(e_uncal::Vector{<:Real}, th228_lines::Vector{<:Unitful.Energy{<:Real}}, window_sizes::Vector{<:Tuple{Unitful.Energy{<:Real}, Unitful.Energy{<:Real}}},; n_bins::Int=15000, quantile_perc::Float64=NaN, binning_peak_window::Unitful.Energy{<:Real}=10.0u"keV", e_unit::Unitful.EnergyUnits=u"keV")

Generic simple_calibration #108

Are you sure you want to change the base?

Generic simple_calibration #108

Conversation

LisaSchlueter commented Dec 16, 2024 • edited Loading

Choose a reason for hiding this comment

LisaSchlueter Dec 16, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

LisaSchlueter commented Dec 16, 2024

codecov bot commented Dec 16, 2024

Codecov Report

theHenks left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

LisaSchlueter commented Dec 17, 2024

fhagemann commented Dec 20, 2024

LisaSchlueter commented Dec 16, 2024 •

edited

Loading

LisaSchlueter Dec 16, 2024 •

edited

Loading