From 26d8ee60712af82906b63238ff016cd841bbdf96 Mon Sep 17 00:00:00 2001 From: Sofia Calgaro Date: Thu, 20 Feb 2025 14:49:09 +0100 Subject: [PATCH 1/5] fixed user naming of nuis pars --- docs/src/config.md | 4 ++-- src/utils.jl | 27 ++++++++++++++------------- test/io/test_get_partitions.jl | 22 +++++++++++++++++++++- 3 files changed, 37 insertions(+), 16 deletions(-) diff --git a/docs/src/config.md b/docs/src/config.md index d019908..08362db 100644 --- a/docs/src/config.md +++ b/docs/src/config.md @@ -70,9 +70,9 @@ Moreover, the config requires the following block for nuisance parameters, ie en In particular, you can set `"correlated": true` if you want to use one variable to correlate the nuisance parameters (eg to speed up the computation times), and `"fixed": false` if you want to include a prior for nuisance parameters (otherwise these parameters they will be fixed to their partition value and not constrained). If a variable is correlated (either `energy_bias` or `energy_res` or `efficiency`), the code will search for a field in the `fit_groups` block of the partitions JSON file to use a correlated variable per each fit group. -In particular, the field has to be specified as: +In particular, the field can be specified as: - `"energy_bias_group_name": "..."` -- `"energy_res_group_name": "..."` +- `"energy_reso_group_name": "..."` - `"efficiency_group_name": "..."` Parameters are then added to the model called `αr_\$name` (for resolution), `αe_\$name` for efficiency and `αb_\$name` for bias. diff --git a/src/utils.jl b/src/utils.jl index a9a1c5d..e88c56c 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -218,11 +218,12 @@ function get_partitions_new(part_path::String) end ## defaults to 'all' - if (haskey("efficiency_group_name", part_data_json["fit_groups"][fit_group])) + if haskey(part_data_json["fit_groups"][fit_group], "efficiency_group_name") append!( arrays["eff_par_name"], [ - "αe_" * Symbol( + Symbol( + "αe_" * part_data_json["fit_groups"][fit_group]["efficiency_group_name"], ), ], @@ -230,32 +231,32 @@ function get_partitions_new(part_path::String) else append!(arrays["eff_par_name"], [:αe_all]) end - - ## defaults to 'all' - if (haskey("energy_scale_group_name", part_data_json["fit_groups"][fit_group])) + if haskey(part_data_json["fit_groups"][fit_group], "energy_bias_group_name") append!( - arrays["energy_reso_name"], + arrays["energy_bias_name"], [ Symbol( - "αr_" * - part_data_json["fit_groups"][fit_group]["energy_scale_group_name"], + "αb_" * + part_data_json["fit_groups"][fit_group]["energy_bias_group_name"], ), ], ) + else + append!(arrays["energy_bias_name"], [:αb_all]) + end + if haskey(part_data_json["fit_groups"][fit_group], "energy_reso_group_name") append!( - arrays["energy_bias_name"], + arrays["energy_reso_name"], [ Symbol( - "αb_" * - part_data_json["fit_groups"][fit_group]["energy_scale_group_name"], + "αr_" * + part_data_json["fit_groups"][fit_group]["energy_reso_group_name"], ), ], ) else append!(arrays["energy_reso_name"], [:αr_all]) - append!(arrays["energy_bias_name"], [:αb_all]) - end end diff --git a/test/io/test_get_partitions.jl b/test/io/test_get_partitions.jl index 81db366..5a96b1d 100644 --- a/test/io/test_get_partitions.jl +++ b/test/io/test_get_partitions.jl @@ -4,6 +4,9 @@ Pkg.instantiate() # instantiate the environment include("../../src/ZeroNuFit.jl") using .ZeroNuFit using TypedTables +using JSON + +Base.exit(code::Int) = throw(ArgumentError("exit code $code")) @testset "test_get_partitions" begin @@ -129,8 +132,25 @@ using TypedTables [2209.1, 2350.0], ], ) - #@test partitions == expected_partitions @test partitions == expected_partitions @test fit_ranges == expected_fit_ranges + # not existing path + config["partitions"] = [joinpath(present_dir, "../inputs/partitions_fake.json")] + @test_throws ArgumentError ZeroNuFit.Utils.get_partitions(config) + + # test handling of correlated parameters with specified key names + partitions_file = joinpath(present_dir, "../inputs/partitions_test_3.json") + partitions, fit_groups, fit_ranges = ZeroNuFit.Utils.get_partitions_new(partitions_file) + partitions_json = JSON.parsefile(partitions_file) + @test partitions[1].energy_bias_name == Symbol( + "αb_" * partitions_json["fit_groups"]["all_phase_II"]["energy_bias_group_name"], + ) + @test partitions[1].energy_reso_name == Symbol( + "αr_" * partitions_json["fit_groups"]["all_phase_II"]["energy_reso_group_name"], + ) + @test partitions[1].eff_name == Symbol( + "αe_" * partitions_json["fit_groups"]["all_phase_II"]["efficiency_group_name"], + ) + end From 1dfd2bcb8203afd9e99c58dffff9e4ebf9306da8 Mon Sep 17 00:00:00 2001 From: Sofia Calgaro Date: Fri, 21 Feb 2025 14:37:06 +0100 Subject: [PATCH 2/5] added build prior test --- src/likelihood.jl | 3 +- test/inputs/partitions_test_3.json | 32 +++++ test/statistics/test_all.jl | 1 + test/statistics/test_build_prior.jl | 179 ++++++++++++++++++++++++++++ 4 files changed, 213 insertions(+), 2 deletions(-) create mode 100644 test/inputs/partitions_test_3.json create mode 100644 test/statistics/test_build_prior.jl diff --git a/src/likelihood.jl b/src/likelihood.jl index a86080e..1a0f2f6 100644 --- a/src/likelihood.jl +++ b/src/likelihood.jl @@ -852,8 +852,7 @@ function build_prior( for name in unique_list priors[Symbol(name)] = Truncated(Normal(0, 1), -Inf, Inf) pretty_names[Symbol(name)] = L"\alpha_{b} (" * split(String(name), "_")[2] * ")" - nuisance_info[string(name)] = - [["combined", "", "", part.detector, 0, 1, -Inf, Inf]] + nuisance_info[string(name)] = [["combined", "", "", 0, 1, -Inf, Inf]] end end diff --git a/test/inputs/partitions_test_3.json b/test/inputs/partitions_test_3.json new file mode 100644 index 0000000..495259c --- /dev/null +++ b/test/inputs/partitions_test_3.json @@ -0,0 +1,32 @@ +{ + "fit_groups": { + "all_phase_II": { + "range": [ + [1930.0, 2098.511], [2108.511, 2113.513], [2123.513, 2190.0] + ], + "model": "uniform", + "bkg_name": "B_gerda_all_pII", + "energy_bias_group_name": "gerdaII_bias", + "energy_reso_group_name": "gerdaII_reso", + "efficiency_group_name": "gerdaII_eff" + } + }, + "partitions": { + "all_phase_II": [ + { + "experiment": "GERDA", + "detector": "ANG4", + "part_name": "part00", + "start_ts": 1450622572, + "end_ts": 1469119346, + "eff_tot": 0.476981, + "eff_tot_sigma": 0.0395483, + "exposure": 0.987039, + "bias": -0.33885135, + "bias_sigma": 0.07638651, + "width": 1.328561380042463, + "width_sigma": 0.08067940552016985 + } + ] + } +} \ No newline at end of file diff --git a/test/statistics/test_all.jl b/test/statistics/test_all.jl index ce6e6e4..1afc0e9 100644 --- a/test/statistics/test_all.jl +++ b/test/statistics/test_all.jl @@ -13,4 +13,5 @@ Test.@testset "statistics" begin include("test_exp_stable.jl") include("test_norm_exponential.jl") include("test_get_signal_bkg_priors.jl") + include("test_build_prior.jl") end diff --git a/test/statistics/test_build_prior.jl b/test/statistics/test_build_prior.jl new file mode 100644 index 0000000..9c3b27d --- /dev/null +++ b/test/statistics/test_build_prior.jl @@ -0,0 +1,179 @@ +using Pkg +Pkg.activate(".") +Pkg.instantiate() +using Random +include("../../src/ZeroNuFit.jl") +using .ZeroNuFit +using IntervalSets +using Distributions +using JSON +using LaTeXStrings + +@testset "test_build_prior" begin + + @info "Testing retrieval of signal & background pdfs (function 'build_prior' in src/likelihood.jl)" + present_dir = @__DIR__ + + config = Dict( + "bkg_only" => false, + "bkg" => Dict( + "correlated" => Dict("range" => "none", "mode" => "none"), + "prior" => "uniform", + "upper_bound" => 0.1, + ), + "signal" => Dict("prior" => "uniform", "upper_bound" => 1000), + "events" => [joinpath(present_dir, "../inputs/events_test.json")], + "partitions" => [joinpath(present_dir, "../inputs/partitions_test.json")], + "output_path" => "tests", + "bat_fit" => Dict("nsteps" => 10000.0, "nchains" => 4), + "nuisance" => Dict( + "efficiency" => Dict("fixed" => true, "correlated" => true), + "energy_bias" => Dict("fixed" => true, "correlated" => false), + "energy_res" => Dict("fixed" => true, "correlated" => false), + ), + "plot" => Dict( + "bandfit_and_data" => false, + "alpha" => 0.3, + "fit_and_data" => false, + "scheme" => "blue", + ), + ) + + # fixed efficiency & energy scale parameters + partitions, fit_ranges = ZeroNuFit.Utils.get_partitions(config) + events = ZeroNuFit.Utils.get_events(config["events"][1], partitions) + part_event_index = ZeroNuFit.Utils.get_partition_event_index(events, partitions) + settings = ZeroNuFit.Utils.get_settings(config) + bkg_name = + JSON.parsefile(config["partitions"][1])["fit_groups"]["all_phase_II"]["bkg_name"] + priors = nothing + pretty_names = nothing + nuisance_info = nothing + try + priors, pretty_names, nuisance_info = + ZeroNuFit.Likelihood.build_prior(partitions, part_event_index, config, settings) + catch e + @error "Error in build_prior: $e" + throw(e) + end + # S & B priors + @test Uniform(0, 1000) == priors.S + @test Uniform(0, 0.1) == priors.B_gerda_all_pII + # pretty names (just check it once, these are unique) + @test string("S [") * L"10^{-27}" * string("yr") * L"^{-1}" * string("]") == + pretty_names[:S] + @test string(bkg_name) * " [cts/keV/kg/yr]" == pretty_names[Symbol(bkg_name)] + @test L"\alpha_{\varepsilon}" == pretty_names[:α] + @test L"\alpha_{r}" == pretty_names[:αr] + @test L"\alpha_{b}" == pretty_names[:αb] + # nuisance_info + @test [] == nuisance_info["α"] + @test [] == nuisance_info["αr"] + @test [] == nuisance_info["αb"] + @test [] == nuisance_info["γ"] + @test [] == nuisance_info["ε"] + @test [] == nuisance_info["ω"] + @test [] == nuisance_info["𝛥"] + + # varying efficiency & energy scale parameters (uncorrelated) + config["nuisance"]["efficiency"]["fixed"] = false + config["nuisance"]["energy_bias"]["fixed"] = false + config["nuisance"]["energy_res"]["fixed"] = false + partitions, fit_ranges = ZeroNuFit.Utils.get_partitions(config) + events = ZeroNuFit.Utils.get_events(config["events"][1], partitions) + part_event_index = ZeroNuFit.Utils.get_partition_event_index(events, partitions) + settings = ZeroNuFit.Utils.get_settings(config) + bkg_name = + JSON.parsefile(config["partitions"][1])["fit_groups"]["all_phase_II"]["bkg_name"] + priors = nothing + pretty_names = nothing + nuisance_info = nothing + priors, pretty_names, nuisance_info = + ZeroNuFit.Likelihood.build_prior(partitions, part_event_index, config, settings) + + # priors + @test Uniform(0, 1000) == priors.S + @test Uniform(0, 0.1) == priors.B_gerda_all_pII + α_min = maximum(-0.476981 ./ 0.0395483) + @test Truncated(Normal(0, 1), α_min, Inf) == priors.αe_all + prior_ω = + Vector{Truncated{Normal{Float64},Continuous,Float64,Float64,Float64}}(undef, 1) + prior_ω[1] = Truncated(Normal(1.328561380042463, 0.08067940552016985), 0, Inf) + @test prior_ω == priors.ω.v + prior_𝛥 = + Vector{Truncated{Normal{Float64},Continuous,Float64,Float64,Float64}}(undef, 1) + prior_𝛥[1] = Truncated(Normal(-0.33885135, 0.07638651), -Inf, Inf) + @test prior_𝛥 == priors.𝛥.v + # pretty names + long_name = "GERDA part00 ANG4" + @test "Energy Resolution " * L"(\omega)" * " " * long_name * " [keV]" == + pretty_names[:ω][1] + @test "Energy Scale Bias " * L"(\Delta)" * " - " * long_name * " [keV]" == + pretty_names[:𝛥][1] + @test L"\alpha_{\varepsilon} (" * "all)" == pretty_names[:αe_all] + # nuisance_info + @test [] == nuisance_info["α"] + @test [] == nuisance_info["αr"] + @test [] == nuisance_info["αb"] + @test [] == nuisance_info["γ"] + @test [] == nuisance_info["ε"] + @test ["GERDA", "part00", "ANG4", 1.328561380042463, 0.08067940552016985, 0, Inf] == + nuisance_info["ω"][1] + @test ["GERDA", "part00", "ANG4", -0.33885135, 0.07638651, -Inf, Inf] == + nuisance_info["𝛥"][1] + @test ["combined", "", "", 0, 1, α_min, Inf] == nuisance_info["αe_all"][1] + + + # correlated energy scale parameters + config["nuisance"]["energy_bias"]["correlated"] = true + config["nuisance"]["energy_res"]["correlated"] = true + partitions, fit_ranges = ZeroNuFit.Utils.get_partitions(config) + events = ZeroNuFit.Utils.get_events(config["events"][1], partitions) + part_event_index = ZeroNuFit.Utils.get_partition_event_index(events, partitions) + settings = ZeroNuFit.Utils.get_settings(config) + bkg_name = + JSON.parsefile(config["partitions"][1])["fit_groups"]["all_phase_II"]["bkg_name"] + priors = nothing + pretty_names = nothing + nuisance_info = nothing + priors, pretty_names, nuisance_info = + ZeroNuFit.Likelihood.build_prior(partitions, part_event_index, config, settings) + + # priors + @test Uniform(0, 1000) == priors.S + @test Uniform(0, 0.1) == priors.B_gerda_all_pII + α_min = maximum(-0.476981 ./ 0.0395483) + @test Truncated(Normal(0, 1), α_min, Inf) == priors.αe_all + @test Truncated(Normal(0, 1), -Inf, Inf) == priors.αb_all + αr_min = maximum(-1.328561380042463 ./ 0.08067940552016985) + @test Truncated(Normal(0, 1), αr_min, Inf) == priors.αr_all + # pretty names + long_name = "GERDA part00 ANG4" + @test L"\alpha_{\varepsilon} (" * "all)" == pretty_names[:αe_all] + @test L"\alpha_{b} (" * "all)" == pretty_names[:αb_all] + @test L"\alpha_{r} (" * "all)" == pretty_names[:αr_all] + # nuisance_info + @test [] == nuisance_info["α"] + @test [] == nuisance_info["αr"] + @test [] == nuisance_info["αb"] + @test [] == nuisance_info["γ"] + @test [] == nuisance_info["ε"] + @test ["combined", "", "", 0, 1, -12.060720688373456, Inf] == nuisance_info["αe_all"][1] + @test ["combined", "", "", 0, 1, -16.467168684210527, Inf] == nuisance_info["αr_all"][1] + @test ["combined", "", "", 0, 1, -Inf, Inf] == nuisance_info["αb_all"][1] + + # no signal = no S prior + config["bkg_only"] = true + partitions, fit_ranges = ZeroNuFit.Utils.get_partitions(config) + events = ZeroNuFit.Utils.get_events(config["events"][1], partitions) + part_event_index = ZeroNuFit.Utils.get_partition_event_index(events, partitions) + settings = ZeroNuFit.Utils.get_settings(config) + bkg_name = + JSON.parsefile(config["partitions"][1])["fit_groups"]["all_phase_II"]["bkg_name"] + priors = nothing + pretty_names = nothing + nuisance_info = nothing + priors, pretty_names, nuisance_info = + ZeroNuFit.Likelihood.build_prior(partitions, part_event_index, config, settings) + @test !haskey(priors, :S) +end From 7383846ac90f9de8aee6bc4035ef6a76767ba9f7 Mon Sep 17 00:00:00 2001 From: Sofia Calgaro Date: Fri, 21 Feb 2025 16:48:25 +0100 Subject: [PATCH 3/5] extensive build prior and get posterior tests --- src/likelihood.jl | 14 +-- src/utils.jl | 2 +- test/inputs/events_test_2.json | 10 ++ test/inputs/partitions_test_4.json | 48 +++++++++ test/statistics/test_build_prior.jl | 147 ++++++++++++++++++++++++++- test/utils/test_all.jl | 1 + test/utils/test_get_par_posterior.jl | 61 +++++++++++ 7 files changed, 268 insertions(+), 15 deletions(-) create mode 100644 test/inputs/events_test_2.json create mode 100644 test/inputs/partitions_test_4.json create mode 100644 test/utils/test_get_par_posterior.jl diff --git a/src/likelihood.jl b/src/likelihood.jl index 1a0f2f6..cb751dc 100644 --- a/src/likelihood.jl +++ b/src/likelihood.jl @@ -755,7 +755,6 @@ function build_prior( :α => L"\alpha_{\varepsilon}", :αr => L"\alpha_{r}", :αb => L"\alpha_{b}", - :γ => [], :ε => [], :ω => [], :𝛥 => [], @@ -776,15 +775,8 @@ function build_prior( end # dictionary with info on the prior parameters - nuisance_info = OrderedDict( - "α" => [], - "αr" => [], - "αb" => [], - "γ" => [], - "ε" => [], - "ω" => [], - "𝛥" => [], - ) + nuisance_info = + OrderedDict("α" => [], "αr" => [], "αb" => [], "ε" => [], "ω" => [], "𝛥" => []) ### EFFICIENCY prior @@ -889,7 +881,7 @@ function build_prior( ) elseif part.signal_name == :gaussian_plus_lowEtail - # let's define some intervals in +-5σ (always with res>0) + # let's define some intervals in +-5σ bias_min = part.bias - 5 * part.bias_sigma bias_max = part.bias + 5 * part.bias_sigma bias[i_new] = diff --git a/src/utils.jl b/src/utils.jl index e88c56c..77b9190 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -55,7 +55,7 @@ Function that retrieves the parameter posterior. - `par`: name of the parameter for which we want to extract the marginalized mode. - `idx`: index for multiparameters """ -function get_par_posterior(samples, par; idx) +function get_par_posterior(samples, par; idx = nothing) pars = [] diff --git a/test/inputs/events_test_2.json b/test/inputs/events_test_2.json new file mode 100644 index 0000000..e0e5b00 --- /dev/null +++ b/test/inputs/events_test_2.json @@ -0,0 +1,10 @@ +{ + "events": [ + { + "experiment": "MJD", + "energy": 2195.2452, + "timestamp": 2.5, + "detector": "DS0" + } + ] +} diff --git a/test/inputs/partitions_test_4.json b/test/inputs/partitions_test_4.json new file mode 100644 index 0000000..00894e6 --- /dev/null +++ b/test/inputs/partitions_test_4.json @@ -0,0 +1,48 @@ +{ + "fit_groups": { + "mjd-DS0": { + "range": [ + [ + 1950.0, + 2098.511 + ], + [ + 2108.511, + 2113.513 + ], + [ + 2123.513, + 2199.1 + ], + [ + 2209.1, + 2350.0 + ] + ], + "signal_name": "fancy_gaussian", + "model": "uniform", + "bkg_name": "mjd-DS0" + } + }, + "partitions": { + "mjd-DS0": [ + { + "experiment": "MJD", + "detector": "DS0", + "part_name": "ds0", + "start_ts": 0, + "end_ts": 5, + "eff_tot": 0.2, + "eff_tot_sigma": 0.01, + "frac": 0.10, + "sigma": 1.10, + "tau": 1.05, + "width": 1.05, + "width_sigma": 0.040, + "exposure": 1.08, + "bias": 0.02, + "bias_sigma": 0.08 + } + ] + } +} diff --git a/test/statistics/test_build_prior.jl b/test/statistics/test_build_prior.jl index 9c3b27d..fedc078 100644 --- a/test/statistics/test_build_prior.jl +++ b/test/statistics/test_build_prior.jl @@ -9,6 +9,8 @@ using Distributions using JSON using LaTeXStrings +Base.exit(code::Int) = throw(ArgumentError("exit code $code")) + @testset "test_build_prior" begin @info "Testing retrieval of signal & background pdfs (function 'build_prior' in src/likelihood.jl)" @@ -70,7 +72,6 @@ using LaTeXStrings @test [] == nuisance_info["α"] @test [] == nuisance_info["αr"] @test [] == nuisance_info["αb"] - @test [] == nuisance_info["γ"] @test [] == nuisance_info["ε"] @test [] == nuisance_info["ω"] @test [] == nuisance_info["𝛥"] @@ -115,7 +116,6 @@ using LaTeXStrings @test [] == nuisance_info["α"] @test [] == nuisance_info["αr"] @test [] == nuisance_info["αb"] - @test [] == nuisance_info["γ"] @test [] == nuisance_info["ε"] @test ["GERDA", "part00", "ANG4", 1.328561380042463, 0.08067940552016985, 0, Inf] == nuisance_info["ω"][1] @@ -156,7 +156,6 @@ using LaTeXStrings @test [] == nuisance_info["α"] @test [] == nuisance_info["αr"] @test [] == nuisance_info["αb"] - @test [] == nuisance_info["γ"] @test [] == nuisance_info["ε"] @test ["combined", "", "", 0, 1, -12.060720688373456, Inf] == nuisance_info["αe_all"][1] @test ["combined", "", "", 0, 1, -16.467168684210527, Inf] == nuisance_info["αr_all"][1] @@ -176,4 +175,146 @@ using LaTeXStrings priors, pretty_names, nuisance_info = ZeroNuFit.Likelihood.build_prior(partitions, part_event_index, config, settings) @test !haskey(priors, :S) + + # eff + config["nuisance"]["efficiency"]["fixed"] = false + config["nuisance"]["efficiency"]["correlated"] = false + partitions, fit_ranges = ZeroNuFit.Utils.get_partitions(config) + events = ZeroNuFit.Utils.get_events(config["events"][1], partitions) + part_event_index = ZeroNuFit.Utils.get_partition_event_index(events, partitions) + settings = ZeroNuFit.Utils.get_settings(config) + priors = nothing + pretty_names = nothing + nuisance_info = nothing + priors, pretty_names, nuisance_info = + ZeroNuFit.Likelihood.build_prior(partitions, part_event_index, config, settings) + # priors + prior_ε = + Vector{Truncated{Normal{Float64},Continuous,Float64,Float64,Float64}}(undef, 1) + prior_ε[1] = Truncated(Normal(0.476981, 0.0395483), 0, 1) + @test prior_ε == priors.ε.v + # pretty names + @test "Efficiency " * L"(\varepsilon)" * " - " * long_name * "" == pretty_names[:ε][1] + # nuisance info + @test ["GERDA", "part00", "ANG4", 0.476981, 0.0395483, 0, 1] == nuisance_info["ε"][1] + + + # bkg shapes parameters (eg linear/exponential) + config["bkg"] = Dict( + "correlated" => Dict("range" => "none", "mode" => "none"), + "prior" => "uniform", + "upper_bound" => 0.1, + "shape" => Dict("name" => "linear", "pars" => Dict("slope" => [-10, 10])), + ) + partitions, fit_ranges = ZeroNuFit.Utils.get_partitions(config) + events = ZeroNuFit.Utils.get_events(config["events"][1], partitions) + part_event_index = ZeroNuFit.Utils.get_partition_event_index(events, partitions) + settings = ZeroNuFit.Utils.get_settings(config) + bkg_shape, bkg_shape_pars = ZeroNuFit.Utils.get_bkg_info(config) + priors = nothing + pretty_names = nothing + nuisance_info = nothing + priors, pretty_names, nuisance_info = ZeroNuFit.Likelihood.build_prior( + partitions, + part_event_index, + config, + settings, + shape_pars = bkg_shape_pars, + ) + # priors + @test Uniform(-10, 10) == priors.B_gerda_all_pII_slope + # pretty names + long_name = "GERDA part00 ANG4" + @test bkg_name * "_slope" == pretty_names[:B_gerda_all_pII_slope] + + # gaussian + low E tail (MJD-like event) + config["bkg"] = Dict( + "correlated" => Dict("range" => "none", "mode" => "none"), + "prior" => "uniform", + "upper_bound" => 0.1, + ) + config["nuisance"]["energy_bias"]["fixed"] = false + config["nuisance"]["energy_bias"]["correlated"] = false + config["nuisance"]["energy_res"]["fixed"] = false + config["nuisance"]["energy_res"]["correlated"] = false + config["events"] = [joinpath(present_dir, "../inputs/events_test_2.json")] + config["partitions"] = [joinpath(present_dir, "../inputs/partitions_test_2.json")] + partitions, fit_ranges = ZeroNuFit.Utils.get_partitions(config) + events = ZeroNuFit.Utils.get_events(config["events"][1], partitions) + part_event_index = ZeroNuFit.Utils.get_partition_event_index(events, partitions) + settings = ZeroNuFit.Utils.get_settings(config) + priors = nothing + pretty_names = nothing + nuisance_info = nothing + priors, pretty_names, nuisance_info = + ZeroNuFit.Likelihood.build_prior(partitions, part_event_index, config, settings) + # priors + prior_𝛥 = + Vector{Truncated{Normal{Float64},Continuous,Float64,Float64,Float64}}(undef, 1) + prior_𝛥[1] = Truncated(Normal(0.02, 0.08), 0.02 - 5 * 0.08, 0.02 + 5 * 0.08) + @test prior_𝛥 == priors.𝛥.v + prior_ω = + Vector{Truncated{Normal{Float64},Continuous,Float64,Float64,Float64}}(undef, 1) + prior_ω[1] = Truncated(Normal(1.05, 0.040), 1.05 - 5 * 0.040, 1.05 + 5 * 0.040) + @test prior_ω == priors.ω.v + # nuisance info + @test ["MJD", "ds0", "DS0", 0.02, 0.08, 0.02 - 5 * 0.08, 0.02 + 5 * 0.08] == + nuisance_info["𝛥"][1] + @test ["MJD", "ds0", "DS0", 1.05, 0.04, 1.05 - 5 * 0.040, 1.05 + 5 * 0.040] == + nuisance_info["ω"][1] + + # same, but with a not-existing signal shape + config["partitions"] = [joinpath(present_dir, "../inputs/partitions_test_4.json")] + partitions, fit_ranges = ZeroNuFit.Utils.get_partitions(config) + events = ZeroNuFit.Utils.get_events(config["events"][1], partitions) + part_event_index = ZeroNuFit.Utils.get_partition_event_index(events, partitions) + settings = ZeroNuFit.Utils.get_settings(config) + @test_throws ArgumentError ZeroNuFit.Likelihood.build_prior( + partitions, + part_event_index, + config, + settings, + ) + + # hierarchical (back to GERDA-like event) + config["bkg"]["correlated"] = Dict("mode" => "lognormal", "range" => [0, 1]) + config["events"] = [joinpath(present_dir, "../inputs/events_test.json")] + config["partitions"] = [joinpath(present_dir, "../inputs/partitions_test.json")] + partitions, fit_ranges = ZeroNuFit.Utils.get_partitions(config) + events = ZeroNuFit.Utils.get_events(config["events"][1], partitions) + part_event_index = ZeroNuFit.Utils.get_partition_event_index(events, partitions) + settings = ZeroNuFit.Utils.get_settings(config) + corr, hier_mode, hier_range = ZeroNuFit.Utils.get_corr_info(config) + priors = nothing + pretty_names = nothing + nuisance_info = nothing + priors, pretty_names, nuisance_info = ZeroNuFit.Likelihood.build_prior( + partitions, + part_event_index, + config, + settings, + hierachical = corr, + hierachical_mode = hier_mode, + hierachical_range = hier_range, + ) + # pretty names + @test "B [cts/keV/kg/yr]" == pretty_names[:B] + @test L"\sigma_B" * string("[cts/keV/kg/yr]") == pretty_names[:σB] + + # unknown hierarchical mode + config["bkg"]["correlated"] = Dict("mode" => "fancynormal", "range" => [0, 1]) + partitions, fit_ranges = ZeroNuFit.Utils.get_partitions(config) + events = ZeroNuFit.Utils.get_events(config["events"][1], partitions) + part_event_index = ZeroNuFit.Utils.get_partition_event_index(events, partitions) + settings = ZeroNuFit.Utils.get_settings(config) + corr, hier_mode, hier_range = ZeroNuFit.Utils.get_corr_info(config) + @test_throws ArgumentError ZeroNuFit.Likelihood.build_prior( + partitions, + part_event_index, + config, + settings, + hierachical = corr, + hierachical_mode = hier_mode, + hierachical_range = hier_range, + ) end diff --git a/test/utils/test_all.jl b/test/utils/test_all.jl index 1ed7459..7862551 100644 --- a/test/utils/test_all.jl +++ b/test/utils/test_all.jl @@ -2,4 +2,5 @@ using Test Test.@testset "likelihood" begin include("test_get_range.jl") + include("test_get_par_posterior.jl") end diff --git a/test/utils/test_get_par_posterior.jl b/test/utils/test_get_par_posterior.jl new file mode 100644 index 0000000..eaaba4d --- /dev/null +++ b/test/utils/test_get_par_posterior.jl @@ -0,0 +1,61 @@ +using Pkg +Pkg.activate(".") # activate the environment +Pkg.instantiate() # instantiate the environment +include("../../src/ZeroNuFit.jl") +using .ZeroNuFit + +struct Sample + v::Dict{Symbol,Any} + weight::Int +end + +@testset "test_get_par_posterior " begin + + @info "Testing function to retrieve samples of a posterior pdf (function 'get_par_posterior ' in src/utils.jl)" + + # testing when idx is nothing + @testset "Test get_par_posterior function" begin + sample1 = Sample(Dict(:alpha => [1, 2, 3], :beta => [4, 5]), 2) + sample2 = Sample(Dict(:alpha => [6, 7, 8], :beta => [9, 10]), 1) + samples = [sample1, sample2] + + # test when idx is nothing (should return the entire value of par) + result = ZeroNuFit.Utils.get_par_posterior(samples, :alpha) + @test result == [1, 2, 3, 1, 2, 3, 6, 7, 8] + end + + #testing when idx is a specific index (should return the indexed value) + @testset "Test get_par_posterior function with idx" begin + sample1 = Sample(Dict(:alpha => [1, 2, 3], :beta => [4, 5]), 2) + sample2 = Sample(Dict(:alpha => [6, 7, 8], :beta => [9, 10]), 1) + samples = [sample1, sample2] + + # test when idx is 2 (should return the 2nd element of :alpha) + result = ZeroNuFit.Utils.get_par_posterior(samples, :alpha, idx = 2) + @test result == [2, 2, 7] + end + + # empty sample list + @testset "Edge case with empty samples" begin + samples = [] + + # test with empty sample list (should return an empty result) + result = ZeroNuFit.Utils.get_par_posterior(samples, :alpha) + @test result == [] + end + + # case where a parameter doesn't exist in the sample + @testset "Test missing parameter" begin + sample1 = Sample(Dict(:alpha => [1, 2, 3], :beta => [4, 5]), 1) + samples = [sample1] + + # test when the parameter doesn't exist (should raise an error or return an empty list) + try + result = ZeroNuFit.Utils.get_par_posterior(samples, :gamma) + @test result == [] + catch e + @test true + end + end + +end From a772c1d0fc10c1dc837ead6d688a24e695820602 Mon Sep 17 00:00:00 2001 From: Sofia Calgaro Date: Fri, 21 Feb 2025 17:02:11 +0100 Subject: [PATCH 4/5] not existing peak shape for resolutions --- test/statistics/test_build_prior.jl | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/test/statistics/test_build_prior.jl b/test/statistics/test_build_prior.jl index fedc078..777a58f 100644 --- a/test/statistics/test_build_prior.jl +++ b/test/statistics/test_build_prior.jl @@ -275,6 +275,18 @@ Base.exit(code::Int) = throw(ArgumentError("exit code $code")) config, settings, ) + config["nuisance"]["energy_bias"]["fixed"] = true + config["nuisance"]["energy_bias"]["correlated"] = false + partitions, fit_ranges = ZeroNuFit.Utils.get_partitions(config) + events = ZeroNuFit.Utils.get_events(config["events"][1], partitions) + part_event_index = ZeroNuFit.Utils.get_partition_event_index(events, partitions) + settings = ZeroNuFit.Utils.get_settings(config) + @test_throws ArgumentError ZeroNuFit.Likelihood.build_prior( + partitions, + part_event_index, + config, + settings, + ) # hierarchical (back to GERDA-like event) config["bkg"]["correlated"] = Dict("mode" => "lognormal", "range" => [0, 1]) From 29d79f5b4450c509a85e81de438c195188747aa4 Mon Sep 17 00:00:00 2001 From: Sofia Calgaro Date: Fri, 21 Feb 2025 17:09:55 +0100 Subject: [PATCH 5/5] added script to plot surviving events --- attic/py_scripts/plot_events.py | 133 ++++++++++++++++++++++++++++++++ 1 file changed, 133 insertions(+) create mode 100644 attic/py_scripts/plot_events.py diff --git a/attic/py_scripts/plot_events.py b/attic/py_scripts/plot_events.py new file mode 100644 index 0000000..88b5ae8 --- /dev/null +++ b/attic/py_scripts/plot_events.py @@ -0,0 +1,133 @@ +""" +Script to plot energy events for different experiments given entry JSON files. +Main Authors: Sofia Calgaro +""" +import os,json,h5py,math +import numpy as np +import matplotlib +import matplotlib.patches as mpatches +import matplotlib.pyplot as plt +from matplotlib.backends.backend_pdf import PdfPages + +### some colours +c_pistachio = (0.58, 0.87, 0.45) +c_columbiablue = (0.61, 0.87, 1.0) +c_dark_columbiablue = (0.61, 0.87*0.85, 1) +c_frenchblue = (0.0, 0.45, 0.73) +c_dark_gray = (0.33, 0.33, 0.33) +c_dark_pistachio = (0.58*0.85, 0.87*0.85, 0.45*0.85) +c_lava = (0.94, 0.01, 0.05) + +### function to read output h5 files with saved posteriors +def read_samples(path,par): + out=None + if os.path.isfile(path): + with h5py.File(path, 'r') as hdf: + + v = hdf[f"v/{par}"][:] + w = hdf[f"weight"][:] + + out= np.repeat(v, w) + return out + else: + print(f"File {path} not found!") + exit() + + +bin_width = 1 +xmin = np.min(1900) +xmax = np.max(3000) +nbins = np.arange(xmin, xmax + bin_width, bin_width) +gm_type = "marginalized_modes" +center = 1930 +expo = 42.25 +fit_range = [[1930, 2099], [2109, 2114], [2124, 2190]] +Qbb = 2039.06 + + +with PdfPages(f"distribution_of_events.pdf") as pdf: + gerdaI_events_json = json.load(open("../legend-0vbb-config/gerda/events_gerda_pI.json")) + gerdaII_events_json = json.load(open("../legend-0vbb-config/gerda/events_gerda.json")) + l200_events_json = json.load(open("../legend-0vbb-config/legend/events_l200_all_05022025.json")) + mjd_events_json = json.load(open("../legend-0vbb-config/mjd/events_mjd_og.json")) + custom_colors = [ + "deepskyblue", + "dodgerblue", + "blue", + "navy", + "darkviolet", + "teal", + "#228b22", + "#ff7f0e", + "chocolate", + "maroon", + ] + + json_dict = { + "ph1_golden": gerdaI_events_json, + "ph1_bege": gerdaI_events_json, + "ph1_silver": gerdaI_events_json, + "ph1_extra": gerdaI_events_json, + "B_gerda_all_pII": gerdaII_events_json, + "B_l200a_Nu24": l200_events_json, + "B_l200a_extra": l200_events_json, + "mjd-DS0": mjd_events_json, + "mjd-mod1": mjd_events_json, + "mjd-mod2": mjd_events_json, + } + + title_names = { + "ph1_golden": "GERDA PhI golden", + "ph1_bege": "GERDA PhI BEGe", + "ph1_silver": "GERDA PhI silver", + "ph1_extra": "GERDA PhI extra", + "B_gerda_all_pII": "GERDA PhII", + "B_l200a_Nu24": "LEGEND Nu24", + "B_l200a_extra": "LEGEND extra", + "mjd-DS0": "MJD DS0", + "mjd-mod1": "MJD module 1", + "mjd-mod2": "MJD module 2", + } + + for idx,bkg_name in enumerate(["ph1_golden", "ph1_bege", "ph1_silver", "ph1_extra", "B_gerda_all_pII", "B_l200a_Nu24", "B_l200a_extra", "mjd-DS0", "mjd-mod1", "mjd-mod2"]): + fig, ax = plt.subplots(figsize=(5,2)) + evt_dict = json_dict[bkg_name] + + events = [] + for entry in evt_dict["events"]: + if bkg_name=="ph1_golden" and entry["experiment"]=="GERDA_pI" and entry["detector"]=="coax" and (entry["timestamp"] > 1320849782-1 and entry["timestamp"] <1369143782+1): events.append(entry["energy"]) + if bkg_name=="ph1_bege" and entry["experiment"]=="GERDA_pI" and entry["detector"]=="bege" and (entry["timestamp"] > 1320849782-1 and entry["timestamp"] <1369143782+1): events.append(entry["energy"]) + if bkg_name=="ph1_silver" and entry["experiment"]=="GERDA_pI" and (entry["timestamp"] > 339945251-1 and entry["timestamp"] <342589688+1): events.append(entry["energy"]) + if bkg_name=="ph1_extra" and entry["experiment"]=="GERDA_pI" and (entry["timestamp"] > 1370007782-1 and entry["timestamp"] <1380548582+1): events.append(entry["energy"]) + + if bkg_name=="B_gerda_all_pII" and entry["experiment"]=="GERDA": events.append(entry["energy"]) + + if bkg_name=="B_l200a_Nu24" and entry["experiment"]=="L200" and entry["detector"][0]!="C" and entry["detector"][0:2]!="V05": events.append(entry["energy"]) + if bkg_name=="B_l200a_extra" and entry["experiment"]=="L200" and (entry["detector"][0]=="C" or entry["detector"][0:2]=="V05"): events.append(entry["energy"]) + + if bkg_name=="mjd-DS0" and entry["experiment"]=="MJD" and "DS0" in entry["detector"]: events.append(entry["energy"]) + if bkg_name=="mjd-mod1" and entry["experiment"]=="MJD": + if entry["detector"] in ["DS1","DS2", "DS3", "DS5a", "DS5b", "DS5c", "DS6a", "DS6b", "DS6c","DS7", "DS8P"] and entry["timestamp"] < 140: events.append(entry["energy"]) + if bkg_name=="mjd-mod2" and entry["experiment"]=="MJD": + if entry["detector"] in ["DS4", "DS5a", "DS5b", "DS5c", "DS6a", "DS6b", "DS6c", "DS8P", "DS8I"] and (entry["timestamp"] >= 140 and entry["timestamp"] <280): events.append(entry["energy"]) + + plt.hist(events, bins=nbins, histtype='stepfilled', density=False, label=F"{title_names[bkg_name]}: {len(events)} events", color=custom_colors[idx]) + plt.xlabel("Energy (keV)") + plt.ylabel('Counts / keV') + plt.axvspan(2099, 2109, facecolor="black", alpha=0.4) + plt.axvspan(2114, 2124, facecolor="black", alpha=0.4) + if "mjd" in bkg_name: + plt.xlim(1950,2350) + plt.axvspan(2199.1, 2209.1, facecolor="black", alpha=0.4) + plt.xticks(np.linspace(1950,2350,4)) + else: + plt.xlim(1930,2190) + plt.xticks(np.linspace(1930,2190,4)) + plt.ylim(0,4.1) + plt.yticks([0,1,2,3,4]) + ax.legend(loc=(0,1), frameon=False, ncol=3) + pdf.savefig(bbox_inches='tight') + plt.close() + + +