Merge a4eb3bf into fea5c06

Project.toml

@@ -25,7 +25,7 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 [extensions]
 LegendMakieExt = ["Makie", "Dates", "FileIO", "Format", "KernelDensity", "LaTeXStrings", "LinearAlgebra", "MathTeXEngine", "StatsBase", "Unitful"]
 LegendMakieLegendDataManagementExt = ["Makie", "LegendDataManagement", "Format", "Measurements", "PropDicts", "TypedTables", "Unitful"]
-LegendMakieMeasurementsExt = ["Makie", "LaTeXStrings", "Measurements"]
+LegendMakieMeasurementsExt = ["Makie", "LaTeXStrings", "Measurements", "Unitful"]
 LegendMakieRadiationDetectorSignalsExt = "RadiationDetectorSignals"
 
 [compat]

ext/LegendMakieMeasurementsExt.jl

@@ -7,8 +7,18 @@ module LegendMakieMeasurementsExt
     import LaTeXStrings
     import Makie
     import Measurements
+    import Unitful
+
+    import LegendMakie: pt, aoecorrectionplot!, energycalibrationplot!
+
+    # function to compose labels when errorbars are scaled
+    function label_errscaling(xerrscaling::Real, yerrscaling::Real)
+        result = String[]
+        xerrscaling != 1 && push!(result, "x-Error ×$(xerrscaling)")
+        yerrscaling != 1 && push!(result, "y-Error ×$(yerrscaling)")
+        isempty(result) ? "" : " ($(join(result,", ")))"
+    end
 
-    import LegendMakie: pt, aoecorrectionplot!
 
     Makie.@recipe(AoECorrectionPlot, report) do scene
         Makie.Attributes(
@@ -31,4 +41,172 @@ module LegendMakieMeasurementsExt
         p
     end
 
+
+    Makie.@recipe(EnergyCalibrationPlot, report, additional_pts) do scene
+        Makie.Attributes(
+            color = LegendMakie.AchatBlue,
+            plot_ribbon = true,
+            xerrscaling = 1,
+            yerrscaling = 1
+        )
+    end
+
+    # Needed for creatings legend using Makie recipes
+    # https://discourse.julialang.org/t/makie-defining-legend-output-for-a-makie-recipe/121567
+    function Makie.get_plots(p::EnergyCalibrationPlot)
+        return p.plots
+    end
+
+    function Makie.plot!(p::EnergyCalibrationPlot{<:Tuple{<:NamedTuple{(:par, :f_fit, :x, :y, :gof)}, <:NamedTuple}})
+
+        report = p.report[]
+        additional_pts = p.additional_pts[]
+        xerrscaling = p.xerrscaling[]
+        yerrscaling = p.yerrscaling[]
+
+        # plot fit
+        xfit = 0:1:1.2*Measurements.value(maximum(report.x))
+        yfit = report.f_fit.(xfit)
+        yfit_m = Measurements.value.(yfit)
+        Makie.lines!(xfit, yfit_m, label = "Best Fit" * (!isempty(report.gof) ? " (p = $(round(report.gof.pvalue, digits=2))| χ²/ndf = $(round(report.gof.chi2min, digits=2)) / $(report.gof.dof))" : ""), color = p.color)
+        if p.plot_ribbon[]
+            Δyfit = Measurements.uncertainty.(yfit)
+            Makie.band!(xfit, yfit_m .- Δyfit, yfit_m .+ Δyfit, color = (p.color[], 0.2))
+        end
+
+        # scatter points with error bars
+        xvalues = Measurements.value.(report.x)
+        yvalues = Measurements.value.(report.y)
+        Makie.errorbars!(p, xvalues, yvalues, Measurements.uncertainty.(report.x) .* xerrscaling, direction = :x, color = :black)
+        Makie.errorbars!(p, xvalues, yvalues, Measurements.uncertainty.(report.y) .* yerrscaling, color = :black)
+        Makie.scatter!(p, xvalues, yvalues, marker = :circle, color = :black, label = "Data" * label_errscaling(xerrscaling, yerrscaling))
+
+        # plot additional points
+        if !isempty(additional_pts)
+            xvalues = Measurements.value.(additional_pts.x)
+            yvalues = Measurements.value.(additional_pts.y)
+            Makie.errorbars!(p, xvalues, yvalues, Measurements.uncertainty.(additional_pts.x) .* xerrscaling, direction = :x, color = :black)
+            Makie.errorbars!(p, xvalues, yvalues, Measurements.uncertainty.(additional_pts.y) .* yerrscaling, color = :black)
+            Makie.scatter!(p, xvalues, yvalues, marker = :circle, color = :silver, strokewidth = 1, strokecolor = :black, label = "Data not used for fit" * label_errscaling(xerrscaling, yerrscaling))
+        end
+
+        p
+    end
+
+
+    function LegendMakie.lplot!(
+            report::NamedTuple{(:par, :f_fit, :x, :y, :gof)};
+            additional_pts::NamedTuple = NamedTuple(),
+            xlims = (0, 1.2*Measurements.value(maximum(report.x))), ylims = nothing,
+            xlabel = "Energy (ADC)", ylabel = "Energy (calibrated)", title::AbstractString = "",
+            show_residuals::Bool = true, plot_ribbon::Bool = true, legend_position = :lt,
+            xerrscaling::Real = 1, yerrscaling::Real = 1, row::Int = 1, col::Int = 1,
+            watermark::Bool = true, final::Bool = (title != ""), kwargs...
+        )
+
+        fig = Makie.current_figure()
+
+        g = Makie.GridLayout(fig[row,col])
+        ax = Makie.Axis(g[1,1],
+            title = title,
+            limits = (xlims, ylims),
+            xlabel = xlabel, ylabel = ylabel,
+        )
+
+        LegendMakie.energycalibrationplot!(ax, report, additional_pts; plot_ribbon, xerrscaling, yerrscaling)
+        legend_position != :none && Makie.axislegend(ax, position = legend_position)
+
+        if !isempty(report.gof) && show_residuals
+
+            ax.xticklabelsize = 0
+            ax.xticksize = 0
+            ax.xlabel = ""
+
+            ax2 = Makie.Axis(g[2,1], yticks = -3:3:3, limits = (xlims,(-5,5)), xlabel = xlabel, ylabel = "Residuals (σ)")
+            LegendMakie.residualplot!(ax2, (x = Measurements.value.(report.x), residuals_norm = report.gof.residuals_norm))
+            # add the additional points
+            if !isempty(additional_pts)
+                Makie.scatter!(ax2, Measurements.value.(additional_pts.x), additional_pts.residuals_norm, 
+                        marker = :circle, color = :silver, strokewidth = 1, strokecolor = :black)
+            end
+
+            # link axis and put plots together
+            Makie.linkxaxes!(ax, ax2)
+            Makie.rowgap!(g, 0)
+            Makie.rowsize!(g, 1, Makie.Auto(4))
+
+            # align ylabels
+            yspace = maximum(Makie.tight_yticklabel_spacing!, (ax, ax2))
+            ax.yticklabelspace = yspace
+            ax2.yticklabelspace = yspace
+        end
+
+        # add watermarks
+        Makie.current_axis!(ax)
+        watermark && LegendMakie.add_watermarks!(; final, kwargs...)
+
+        fig
+    end
+
+
+    function LegendMakie.lplot!(
+            report::NamedTuple{(:par, :f_fit, :x, :y, :gof, :e_unit, :type)}; 
+            additional_pts::NamedTuple = NamedTuple(), kwargs...
+        )
+
+        if report.type != :cal
+            @warn "Unknown calibration type $(report.type), no plot generated."
+            return Makie.current_figure()
+        end
+
+        LegendMakie.lplot!(
+            report[(:par, :f_fit, :x, :y, :gof)],
+            additional_pts = if !isempty(additional_pts)
+                # strip the units from the additional points
+                μ_strip = Unitful.unit(first(additional_pts.μ)) != Unitful.NoUnits ? Unitful.ustrip.(report.e_unit, additional_pts.μ) : additional_pts.μ
+                p_strip = Unitful.unit(first(additional_pts.peaks)) != Unitful.NoUnits ? Unitful.ustrip.(report.e_unit, additional_pts.peaks) : additional_pts.peaks    
+                μ_cal = report.f_fit.(μ_strip)
+                (x = μ_strip, y = p_strip, residuals_norm = (Measurements.value.(μ_cal) .- Measurements.value.(p_strip))./ Measurements.uncertainty.(μ_cal))
+            else
+                NamedTuple()
+            end,
+            xlabel = "Energy (ADC)", ylabel = "Energy ($(report.e_unit))", xlims = (0, 1.1*Measurements.value(maximum(report.x))); kwargs...
+        )
+    end
+
+    function LegendMakie.lplot!(
+            report::NamedTuple{(:par, :f_fit, :x, :y, :gof, :e_unit, :qbb, :type)};
+            additional_pts::NamedTuple = NamedTuple(), xlims = (0,3000), title::AbstractString = "",
+            kwargs...
+        )
+
+        if report.type != :fwhm
+            @warn "Unknown calibration type $(report.type), no plot generated."
+            return Makie.current_figure()
+        end
+
+        fig = LegendMakie.lplot!(
+            report[(:par, :f_fit, :x, :y, :gof)],
+            additional_pts = if !isempty(additional_pts)
+                fwhm_cal = report.f_fit.(Unitful.ustrip.(additional_pts.peaks))
+                (x = Unitful.ustrip.(report.e_unit, additional_pts.peaks), y = Unitful.ustrip.(report.e_unit, additional_pts.fwhm),
+                    residuals_norm = (Measurements.value.(fwhm_cal .- Unitful.ustrip.(report.e_unit, additional_pts.fwhm))) ./ Measurements.uncertainty.(fwhm_cal))
+            else
+                NamedTuple()
+            end,
+            xlabel = "Energy ($(report.e_unit))", ylabel = "FWHM ($(report.e_unit))",
+            legend_position = :none; xlims, title, kwargs...
+        )
+
+        ax = first(fig.content)
+        Makie.current_axis!(ax)
+        qbb = Unitful.ustrip(report.e_unit, Measurements.value(report.qbb))
+        Δqbb = Unitful.ustrip(report.e_unit, Measurements.uncertainty(report.qbb))
+        Makie.hlines!(ax, [qbb], color = LegendMakie.CoaxGreen, label = LaTeXStrings.latexstring("\\fontfamily{Roboto}Q_{\\beta \\beta}:" * " $(round(Unitful.ustrip(report.e_unit, report.qbb), digits=2))\\;\\text{$(report.e_unit)}"))
+        Makie.band!(ax, range(xlims..., length = 2), fill(qbb - Δqbb, 2), fill(qbb + Δqbb, 2), color = (LegendMakie.CoaxGreen, 0.2))
+        Makie.axislegend(ax, position = :lt)
+
+        fig
+    end
+
 end

ext/recipes/lplot.jl

@@ -367,6 +367,8 @@ end
 # fallback method: use Makie.plot!
 function LegendMakie.lplot!(args...; watermark::Bool = false, kwargs...)
 
+    @info "No `LegendMakie` plot recipe found for this set of arguments. Using `Makie.plot!`"
+
     fig = Makie.current_figure()
     ax = isnothing(Makie.current_axis()) ? Makie.Axis(fig[1,1]) : Makie.current_axis()
 

src/lplot.jl

@@ -25,6 +25,7 @@ export lsavefig
 
 # recipes
 function residualplot! end
+function energycalibrationplot! end
 function aoecorrectionplot! end
 function parameterplot! end
 

test/test_lplot.jl

@@ -23,7 +23,7 @@ using Test
 
         # test default watermark
         ax = Axis(fig[1,1])
-        @test_nowarn lplot!(StatsBase.fit(StatsBase.Histogram, randn(10000)))
+        lplot!(StatsBase.fit(StatsBase.Histogram, randn(10000)))
         @test_nowarn LegendMakie.add_watermarks!()
 
         # test alternative watermark
@@ -39,16 +39,29 @@ using Test
         end
 
         @testset "Energy calibration" begin
-            ecal = vcat(rand(Distributions.Exponential(50_000),97_500), 261_450 .+ 200 .* randn(2_000), 210_350 .+ 150 .* randn(300), 159_300 .+ 100 .* randn(200))
-            result_simple, report_simple = LegendSpecFits.simple_calibration(ecal, [1592.513, 2103.512, 2614.511]u"keV", [25, 25, 35]u"keV", [25, 25, 30]u"keV", calib_type = :th228)
+            ecal = filter(x -> x <= 262_000, vcat(rand(Distributions.Exponential(70_000),97_000), 261_450 .+ 200 .* randn(2_000), 210_350 .+ 185 .* randn(500), 159_300 .+ 170 .* randn(500)))
+            lines = [:Tl208DEP, :Tl208SEP, :Tl208FEP]
+            energies = [1592.513, 2103.512, 2614.511]u"keV"
+            result_simple, report_simple = LegendSpecFits.simple_calibration(ecal, energies, [25, 25, 35]u"keV", [25, 25, 30]u"keV", calib_type = :th228)
             @testset "Simple energy calibration" begin
                 @test_nowarn lplot(report_simple)
                 @test_nowarn lplot(report_simple, cal = false)
             end
             m_cal_simple = result_simple.c
-            result_fit, report_fit = LegendSpecFits.fit_peaks(result_simple.peakhists, result_simple.peakstats, [:Tl208DEP, :Tl208SEP, :Tl208FEP]; e_unit=result_simple.unit, calib_type=:th228, m_cal_simple=m_cal_simple)
+            result_fit, report_fit = LegendSpecFits.fit_peaks(result_simple.peakhists, result_simple.peakstats, lines; e_unit=result_simple.unit, calib_type=:th228, m_cal_simple=m_cal_simple)
+            @testset "Fit peaks for energy calibration" begin
+                @test_nowarn lplot(report_fit, figsize = (600, 400*length(report_fit)), watermark = false)
+            end
+            μ_fit = getfield.(getindex.(Ref(result_fit), lines), :centroid)
+            result_calib, report_calib = LegendSpecFits.fit_calibration(1, μ_fit, energies)
             @testset "Fit energy calibration" begin
-                @test_nowarn lplot(report_fit, figsize = (600, 400*length(report_fit)), title = "Test title")
+                @test_nowarn lplot(report_calib, xerrscaling=10, yerrscaling=10, additional_pts=(μ = [100_000], peaks = [1000u"keV"]), title = "Test")
+            end
+            f_cal_widths(x) = report_calib.f_fit(x) .* report_calib.e_unit .- first(report_calib.par)
+            fwhm_fit = f_cal_widths.(getfield.(getindex.(Ref(result_fit), lines), :fwhm))
+            result_fwhm, report_fwhm = LegendSpecFits.fit_fwhm(1, energies, fwhm_fit, uncertainty=true)
+            @testset "FWHM energy calibration" begin
+                @test_nowarn lplot(report_fwhm, additional_pts=(peaks = [1000u"keV"], fwhm = [3.5u"keV"]), title = "Test")
             end
         end