diff --git a/scripts/make_event_rate_tables.py b/scripts/make_event_rate_tables.py new file mode 100644 index 0000000..071b38c --- /dev/null +++ b/scripts/make_event_rate_tables.py @@ -0,0 +1,175 @@ +# %% +# %load_ext autoreload +# %autoreload 2 +# %% +import sys +import os + +from matplotlib import pyplot as plt +import numpy as np +import toml +# %% +sys.path.append("../") +import logging +from microfit.fileio import from_json, to_json +from microfit.analysis import MultibandAnalysis +from microfit.parameters import Parameter, ParameterSet +from unblinding_functions import ( + plot_chi2_distribution, +) +from microfit.histogram import MultiChannelHistogram +from microfit.category_definitions import get_category_label +logging.basicConfig(level=logging.INFO) + +# %% +config_file = "../config_files/old_model_ana_with_detvars.toml" +output_dir = "../old_model_ana_output_with_3a/" + +analysis = MultibandAnalysis.from_toml( + config_file, + logging_level=logging.DEBUG, + overwrite_cached_df=False, + overwrite_cached_df_detvars=False, + output_dir=output_dir, +) +# %% +h0_params = ParameterSet([Parameter("signal_strength", 0.0)]) +h1_params = ParameterSet([Parameter("signal_strength", 1.0)]) +analysis.set_parameters(h0_params) +data_hist = analysis.get_data_hist() +ext_hist = analysis.get_data_hist(type="ext") +assert ext_hist is not None +assert data_hist is not None +total_mc_hist_h0_constrained = analysis.get_mc_hist( + include_multisim_errors=True, + use_sideband=True, + add_precomputed_detsys=True +) +total_pred_hist_h0_constrained = total_mc_hist_h0_constrained + ext_hist +total_mc_hist_unconstrained = analysis.get_mc_hist( + include_multisim_errors=True, + use_sideband=False, + add_precomputed_detsys=True +) +total_pred_hist_h0_unconstrained = total_mc_hist_unconstrained + ext_hist + +analysis.set_parameters(h1_params) + +mc_hists = analysis.get_mc_hists(category_column="paper_category") + +total_mc_hist_h1 = analysis.get_mc_hist( + include_multisim_errors=True, + use_sideband=True, + add_precomputed_detsys=True, +) +total_pred_hist_h1 = total_mc_hist_h1 + ext_hist +# %% +print("Total H0 pred (constrained):") +assert isinstance(total_pred_hist_h0_constrained, MultiChannelHistogram) +for channel_hist in total_pred_hist_h0_constrained: + print(f"{channel_hist.binning.label}: {channel_hist.sum():.3f} +/- {channel_hist.sum_std():.3f}") +print("Total H0 pred (unconstrained):") +assert isinstance(total_pred_hist_h0_unconstrained, MultiChannelHistogram) +for channel_hist in total_pred_hist_h0_unconstrained: + print(f"{channel_hist.binning.label}: {channel_hist.sum():.3f} +/- {channel_hist.sum_std():.3f}") +print("Split by category:") +for category, category_mc_hist in mc_hists.items(): + category_label = get_category_label(category_column="paper_category", category=category) + for channel_hist in category_mc_hist: + print(f"{category_label} {channel_hist.binning.label}: {channel_hist.sum():.3f}") +# treat EXT also as a category +for channel_hist in ext_hist: + print(f"EXT {channel_hist.binning.label}: {channel_hist.sum():.3f}") +print("Total H1 pred:") +assert isinstance(total_pred_hist_h1, MultiChannelHistogram) +for channel_hist in total_pred_hist_h1: + print(f"{channel_hist.binning.label}: {channel_hist.sum():.3f} +/- {channel_hist.sum_std():.3f}") +print("Data counts:") +assert isinstance(data_hist, MultiChannelHistogram) +for channel_hist in data_hist: + # No errors because this is data + print(f"{channel_hist.binning.label}: {channel_hist.sum():.0f}") +# %% +# Turn the printout into a table, with one column for each channel +# and one row for each of the following: +# - Data counts +# - Total pred. H0 counts (constrained), with uncertainties +# - Total pred. H0 counts (unconstrained), with uncertainties +# - Total pred. H1 counts, with uncertainties +# - Total pred. H0 counts split by category (no uncertainties) +# - Data counts (no uncertainties) + +# First, get the channel labels +channel_labels = [channel_hist.binning.label for channel_hist in data_hist] +print("&" + "&".join(channel_labels) + r"\\") # type: ignore +# Data counts +print("Data counts", end="") +for channel_hist in data_hist: + print(f"&{channel_hist.sum():.0f}", end="") +print(r"\\") +# MC H0 counts (constrained) +print("Total $H_0$ pred. (constrained)", end="") +for channel_hist in total_pred_hist_h0_constrained: + print(f"&{channel_hist.sum():.1f} $\\pm$ {channel_hist.sum_std():.1f}", end="") +print(r"\\") +# MC H0 counts (unconstrained) +print("Total $H_0$ pred. (unconstrained)", end="") +for channel_hist in total_pred_hist_h0_unconstrained: + print(f"&{channel_hist.sum():.1f} $\\pm$ {channel_hist.sum_std():.1f}", end="") +print(r"\\") +# MC H1 counts +print("Total $H_1$ pred.", end="") +for channel_hist in total_pred_hist_h1: + print(f"&{channel_hist.sum():.1f} $\\pm$ {channel_hist.sum_std():.1f}", end="") +print(r"\\") +# Predicted H0 counts split by category. +print("$H_0$ pred. split by category&&\\\\") + +# treat EXT also as a category +print("EXT", end="") +for channel_hist in ext_hist: + print(f"&{channel_hist.sum():.1f}", end="") +print(r"\\") +for category, category_mc_hist in mc_hists.items(): + category_label = get_category_label(category_column="paper_category", category=category) + print(category_label, end="") + for channel_hist in category_mc_hist: + print(f"&{channel_hist.sum():.1f}", end="") + print(r"\\") + +# %% +# Now we slighly hack the configuration to use the signal model 2 +ana_config = toml.load(config_file) +for gen_config in ana_config["generator"]: + if gen_config["load_runs"].get("load_lee", False): + gen_config["load_runs"]["use_new_signal_model"] = True + +analysis = MultibandAnalysis( + configuration=ana_config, + output_dir=output_dir, + logging_level=logging.DEBUG, + overwrite_cached_df=False, + overwrite_cached_df_detvars=False, +) +# %% +analysis.plot_signals() +# %% +analysis.set_parameters(h1_params) +total_mc_hist_h1_new_model = analysis.get_mc_hist( + include_multisim_errors=True, + use_sideband=True, + add_precomputed_detsys=True, +) +total_pred_hist_h1_new_model = total_mc_hist_h1_new_model + ext_hist +mc_hists_new_model = analysis.get_mc_hists(category_column="paper_category") +# %% +print("Total H1 pred (new model, constrained):") +assert isinstance(total_pred_hist_h1_new_model, MultiChannelHistogram) +for channel_hist in total_pred_hist_h1_new_model: + print(f"{channel_hist.binning.label}: {channel_hist.sum():.3f} +/- {channel_hist.sum_std():.3f}") +print("Split by category:") +for category, category_mc_hist in mc_hists_new_model.items(): + category_label = get_category_label(category_column="paper_category", category=category) + for channel_hist in category_mc_hist: + print(f"{category_label} {channel_hist.binning.label}: {channel_hist.sum():.3f}") +# %% diff --git a/scripts/old_model_normalization_test.py b/scripts/old_model_normalization_test.py index ad21bba..9e9230f 100644 --- a/scripts/old_model_normalization_test.py +++ b/scripts/old_model_normalization_test.py @@ -67,6 +67,18 @@ def run_unblinding(signal_channels, constraint_channels, plot_suffix, plot_title figsize=hist_plot_figsize, override_channel_titles=override_channel_titles, ) + analysis.set_parameters(h0_params) + data_hist = analysis.get_data_hist() + assert data_hist is not None + total_mc_hist = analysis.get_mc_hist( + include_multisim_errors=True, use_sideband=True, add_precomputed_detsys=True + ) + print("Data histogram") + print(data_hist.bin_counts) + print(data_hist.std_devs) + print("Total MC histogram") + print(total_mc_hist.bin_counts) + print(total_mc_hist.std_devs) chi2_at_h0 = analysis.get_chi2_at_hypothesis(h0_params)