diff --git a/models/__init__.py b/models/__init__.py index 63e7476566a3bededd86feaa540e5bf721594cf5..92e327126ae088b14bb9631d2a67fd35fc3c0b98 100644 --- a/models/__init__.py +++ b/models/__init__.py @@ -6,12 +6,14 @@ from .dayabay_v0 import model_dayabay_v0 from .dayabay_v0b import model_dayabay_v0b from .dayabay_v0c import model_dayabay_v0c from .dayabay_v0d import model_dayabay_v0d +# from .dayabay_v0e import model_dayabay_v0e _dayabay_models = { "v0": model_dayabay_v0, "v0b": model_dayabay_v0b, "v0c": model_dayabay_v0c, "v0d": model_dayabay_v0d, + # "v0e": model_dayabay_v0e, } _available_sources = ("tsv", "hdf5", "root", "npz") diff --git a/models/dayabay_v0d.py b/models/dayabay_v0d.py index 9b692459837cac6428e6d0282cddef9c0e0081a7..7defefed7a6611c55f4ad10d81816713b4a2cdc8 100644 --- a/models/dayabay_v0d.py +++ b/models/dayabay_v0d.py @@ -59,6 +59,10 @@ class model_dayabay_v0d: Purpose: - introduce alternative data inputs - provide configuration to switch between inputs and compare them + - introduce dataset A: + - livetimes, efficiencies, accidentals rates + - muon decay background + - proper correlations between background rate parameters Attributes ---------- @@ -717,7 +721,7 @@ class model_dayabay_v0d: ) load_parameters( path="bkg.uncertainty_scale_by_site", - load=path_parameters / "bkg_rate_uncertainty_scale_site.yaml", + load=path_parameters / "bkg_rate_uncertainty_scale_site_dataset_a.yaml", replicate=combinations["site.period"], ignore_keys=inactive_backgrounds, ) @@ -2671,12 +2675,23 @@ class model_dayabay_v0d: unused_keys = list(labels_mk.walkjoinedkeys()) may_ignore = { + # future "bkg.spectrum_per_day", "statistic.nuisance.parts.bkg.rate.acc", "statistic.nuisance.parts.bkg.rate.amc", "statistic.nuisance.parts.bkg.rate.lihe", "statistic.nuisance.parts.bkg.rate.fastn", "statistic.nuisance.parts.bkg.rate.muonx", + # past + "daily_data.detector.rate_acc", + "daily_data.detector.rate_acc_s_day", + "daily_data.detector.num_acc_s_day", + "bkg.count_fixed", + "bkg.count", + "bkg.spectrum.muonx", + "bkg.spectrum_shape.muonx", + "statistic.nuisance.parts.bkg", + "summary", } for key_may_ignore in may_ignore: for i, key_unused in reversed(tuple(enumerate(unused_keys))): diff --git a/models/dayabay_v0d_labels.yaml b/models/dayabay_v0d_labels.yaml index bf874c943351388c2a8b1d9e55dadfceec38650e..e237077da58759722bf8cb2f971a89fdf34b258a 100644 --- a/models/dayabay_v0d_labels.yaml +++ b/models/dayabay_v0d_labels.yaml @@ -114,6 +114,12 @@ daily_data: livetime: group: text: "Daily detector wall clock livetime {key} [s]" + rate_acc: + group: + text: "Daily rate of accidental background events [#/day]" + num_acc_s_day: + group: + text: "Daily # of accidental background events 脳 day/s [#]" reactor: power: group: @@ -221,22 +227,66 @@ bkg: latex: "Reconstructed energy $E_{\\rm rec}$ bin edges [MeV] (for background events)" axis: "$E_{\\rm rec}$ [MeV]" plotmethod: 'none' + count_acc_fixed_s_day: + group: + text: "# of accidental background events in {key} 脳 day/s [#]" + count_fixed: + acc: + group: + text: "# of accidental background events in {index[0]} during {index[1]} [#]\\nbefore nuisance" + alphan: + group: + text: "# of 鹿鲁C(伪,n)鹿鈦禣 background events in {index[0]} during {index[1]} [#]\\nbefore nuisance" + amc: + group: + text: "# of 虏鈦绰笰m鹿鲁C related background events in {index[0]} during {index[1]} [#]\\nbefore nuisance" + fastn: + group: + text: "# of fast neutron background events in {index[0]} during {index[1]} [#]\\nbefore nuisance" + lihe: + group: + text: "# of 鈦筁i/鈦窰e background events in {index[0]} during {index[1]} [#]\\nbefore nuisance" + muonx: + group: + text: "# of 渭 decay background events in {index[0]} during {index[1]} [#]\\nbefore nuisance" + count: + acc: + group: + text: "# of accidental background events in {index[0]} during {index[1]} [#]" + alphan: + group: + text: "# of 鹿鲁C(伪,n)鹿鈦禣 background events in {index[0]} during {index[1]} [#]" + amc: + group: + text: "# of 虏鈦绰笰m鹿鲁C related background events in {index[0]} during {index[1]} [#]" + fastn: + group: + text: "# of fast neutron background events in {index[0]} during {index[1]} [#]" + lihe: + group: + text: "# of 鈦筁i/鈦窰e background events in {index[0]} during {index[1]} [#]" + muonx: + group: + text: "# of 渭 decay background events in {index[0]} during {index[1]} [#]" spectrum: acc: group: - text: "# of accidental background events in {key} [#]" + text: "# of accidental background events in {index[0]} during {index[1]} [#]" alphan: group: - text: "# of 鹿鲁C(伪,n)鹿鈦禣 background events in {key} [#]" + text: "# of 鹿鲁C(伪,n)鹿鈦禣 background events in {index[0]} during {index[1]} [#]" amc: group: - text: "# of 虏鈦绰笰m鹿鲁C related background events in {key} [#]" + text: "# of 虏鈦绰笰m鹿鲁C related background events in {index[0]} during {index[1]} [#]" fastn: group: - text: "# of fast neutron background events in {key} [#]" + text: "# of fast neutron background events in {index[0]} during {index[1]} [#]" lihe: group: - text: "# of 鈦筁i/鈦窰e background events in {key} [#]" + text: "# of 鈦筁i/鈦窰e background events in {index[0]} during {index[1]} [#]" + muonx: + group: + text: "# of 渭 decay background events in {index[0]} during {index[1]} [#]" spectrum_per_day: acc: group: @@ -269,6 +319,9 @@ bkg: lihe: group: text: "Shape of spectrum of 鈦筁i/鈦窰e background events in {key}" + muonx: + group: + text: "Shape of spectrum of 渭 decay background events in {key}" oscprob: group: text: "谓虆(e) survival probability {key}" @@ -642,3 +695,56 @@ statistic: text: "蠂虏 nuisance term for parameters for background rate\\nfast neutrons" lihe: text: "蠂虏 nuisance term for parameters for background rate\\n鈦筁i/鈦窰e" + rate_scale: + acc: + text: "蠂虏 nuisance term for parameters for background rate\\naccidentals" + uncertainty_scale: + amc: + text: "蠂虏 nuisance term for parameters for background rate\\n虏鈦绰笰m鹿鲁C" + uncertainty_scale_by_site: + fastn: + text: "蠂虏 nuisance term for parameters for background rate\\nfast neutrons" + lihe: + text: "蠂虏 nuisance term for parameters for background rate\\n鈦筁i/鈦窰e" + muonx: + text: "蠂虏 nuisance term for parameters for background rate\\n渭 decay" +summary: + total: + bkg_count: + group: + text: "Total # of {index[0]} background events in {index[1]}" + bkg_rate_s: + group: + text: "Total rate of {index[0]} background events in {index[1]}" + bkg_rate: + group: + text: "Total rate {index[0]} background events in {index[1]}" + eff: + group: + text: "Total detector {index[0]} efficiency" + efflivetime: + group: + text: "Total detector {index[0]} effective livetime [s]" + livetime: + group: + text: "Total detector {index[0]} wall clock livetime [s]" + periods: + bkg_count: + group: + text: "Total # of {index[0]} background events in {index[2]} during {index[1]}" + bkg_rate_s: + group: + text: "Total rate of {index[0]} background events in {index[2]} during {index[1]}" + bkg_rate: + group: + text: "Total rate {index[0]} background events in {index[2]} during {index[1]}" + eff: + group: + text: "Total detector {index[1]} efficiency during {index[0]}" + efflivetime: + group: + text: "Total detector {index[1]} effective livetime during {index[0]} [s]" + livetime: + group: + text: "Total detector {index[1]} wall clock livetime during {index[0]} [s]" + diff --git a/scripts/covmatrix_plot.py b/scripts/covmatrix_plot.py index 79d7ff5658a0792b862b7fe27e01d3781e853fa8..f306aa15ff8737c43ee70810dac683b88ef5daaf 100755 --- a/scripts/covmatrix_plot.py +++ b/scripts/covmatrix_plot.py @@ -8,6 +8,7 @@ from typing import TYPE_CHECKING from h5py import File, Group from matplotlib import pyplot as plt from matplotlib.backends.backend_pdf import PdfPages +from matplotlib.colors import SymLogNorm from numpy import arange, diagonal, isnan, quantile from numpy.typing import NDArray @@ -17,6 +18,8 @@ from dagflow.tools.logger import logger if TYPE_CHECKING: from typing import Literal, Mapping, Sequence + from matplotlib.collections import QuadMesh + def main(opts: Namespace) -> None: cmap = "RdBu_r" @@ -26,7 +29,6 @@ def main(opts: Namespace) -> None: model = group["model"][:] edges = group["edges"][:] - blocksize = edges.size - 1 sfile, sgroup, smodel = None, None, None if opts.subtract: @@ -48,6 +50,13 @@ def main(opts: Namespace) -> None: except KeyError: elements = None + blocksize = edges.size - 1 + block_selected = 1 + block_selection = slice( + blocksize * block_selected, blocksize * (block_selected + 1) + ) + block_name = elements[block_selected] + if opts.output is not None: ofile = opts.output and opts.output[0] or opts.input.replace(".hdf5", ".pdf") assert ofile != opts.input and ofile.endswith(".pdf") @@ -156,8 +165,8 @@ def main(opts: Namespace) -> None: ylabel="bin", title=f"Covariance matrix {name}{title_suffix}", ) - pcolor_with_blocks(matrix_cov, blocks=elements, cmap=cmap) - if pdf: + hm = pcolor_with_blocks(matrix_cov, blocks=elements, cmap=cmap) + if pdf and hm is not None: pdf.savefig() if vmax is not None: @@ -166,12 +175,23 @@ def main(opts: Namespace) -> None: 111, xlabel="", ylabel="bin", - title=f"Covariance matrix (truncated) {name}{title_suffix}", + title=f"Covariance matrix {name} (truncated){title_suffix}", ) - pcolor_with_blocks(matrix_cov, blocks=elements, cmap=cmap, vmax=vmax) - if pdf: + hm = pcolor_with_blocks(matrix_cov, blocks=elements, cmap=cmap, vmax=vmax) + if pdf and hm is not None: pdf.savefig() + plt.figure(figsize=figsize_2d) + ax = plt.subplot( + 111, + xlabel="", + ylabel="bin", + title=f"Covariance matrix {name} (symlog){title_suffix}", + ) + hm = pcolor_with_blocks(matrix_cov, blocks=elements, cmap=cmap, symlog=True) + if pdf and hm is not None: + pdf.savefig() + plt.figure(figsize=figsize_2d) ax = plt.subplot( 111, @@ -179,8 +199,19 @@ def main(opts: Namespace) -> None: ylabel="bin", title=f"Covariance matrix {name} (blocks){title_suffix}", ) - pcolor_with_blocks(bmatrix_cov, blocks=elements, cmap=cmap) - if pdf: + hm = pcolor_with_blocks(bmatrix_cov, blocks=elements, cmap=cmap) + if pdf and hm is not None: + pdf.savefig() + + plt.figure(figsize=figsize_2d) + ax = plt.subplot( + 111, + xlabel="", + ylabel="bin", + title=f"Covariance matrix {name} (blocks, symlog){title_suffix}", + ) + hm = pcolor_with_blocks(bmatrix_cov, blocks=elements, cmap=cmap, symlog=True) + if pdf and hm is not None: pdf.savefig() plt.figure(figsize=figsize_2d) @@ -188,12 +219,28 @@ def main(opts: Namespace) -> None: 111, xlabel="", ylabel="bin", - title=f"Covariance matrix {name} ({elements[0]}){title_suffix}", + title=f"Covariance matrix {name} (block {block_name}){title_suffix}", ) - pcolor_with_blocks( - matrix_cov[:blocksize, :blocksize], blocks=elements[:1], cmap=cmap + hm = pcolor_with_blocks( + matrix_cov[block_selection, block_selection], blocks=elements[:1], cmap=cmap ) - if pdf: + if pdf and hm is not None: + pdf.savefig() + + plt.figure(figsize=figsize_2d) + ax = plt.subplot( + 111, + xlabel="", + ylabel="bin", + title=f"Covariance matrix {name} (block {block_name}, symlog){title_suffix}", + ) + hm = pcolor_with_blocks( + matrix_cov[block_selection, block_selection], + blocks=elements[:1], + cmap=cmap, + symlog=True, + ) + if pdf and hm is not None: pdf.savefig() plt.figure(figsize=figsize_2d) @@ -203,8 +250,8 @@ def main(opts: Namespace) -> None: ylabel="bin", title=rf"Relative covariance matrix {name}, %虏{title_suffix}", ) - pcolor_with_blocks(matrix_cov_rel, blocks=elements, cmap=cmap) - if pdf: + hm = pcolor_with_blocks(matrix_cov_rel, blocks=elements, cmap=cmap) + if pdf and hm is not None: pdf.savefig() if vmax_rel is not None: @@ -213,12 +260,12 @@ def main(opts: Namespace) -> None: 111, xlabel="", ylabel="bin", - title=rf"Relative covariance matrix {name}, %虏{title_suffix}", + title=rf"Relative covariance matrix {name} (truncated), %虏{title_suffix}", ) - pcolor_with_blocks( + hm = pcolor_with_blocks( matrix_cov_rel, blocks=elements, cmap=cmap, vmax=vmax_rel ) - if pdf: + if pdf and hm is not None: pdf.savefig() plt.figure(figsize=figsize_2d) @@ -226,12 +273,41 @@ def main(opts: Namespace) -> None: 111, xlabel="", ylabel="bin", - title=f"Relative covariance matrix {name} ({elements[0]}){title_suffix}", + title=rf"Relative covariance matrix {name} (symlog), %虏{title_suffix}", + ) + hm = pcolor_with_blocks(matrix_cov_rel, blocks=elements, cmap=cmap, symlog=True) + if pdf and hm is not None: + pdf.savefig() + + plt.figure(figsize=figsize_2d) + ax = plt.subplot( + 111, + xlabel="", + ylabel="bin", + title=f"Relative covariance matrix {name} (block {block_name}){title_suffix}", ) - pcolor_with_blocks( - matrix_cov_rel[:blocksize, :blocksize], blocks=elements[:1], cmap=cmap + hm = pcolor_with_blocks( + matrix_cov_rel[block_selection, block_selection], + blocks=elements[:1], + cmap=cmap, ) - if pdf: + if pdf and hm is not None: + pdf.savefig() + + plt.figure(figsize=figsize_2d) + ax = plt.subplot( + 111, + xlabel="", + ylabel="bin", + title=f"Relative covariance matrix {name} (block {block_name}, symlog){title_suffix}", + ) + hm = pcolor_with_blocks( + matrix_cov_rel[block_selection, block_selection], + blocks=elements[:1], + cmap=cmap, + symlog=True, + ) + if pdf and hm is not None: pdf.savefig() plt.figure(figsize=figsize_2d) @@ -241,8 +317,8 @@ def main(opts: Namespace) -> None: ylabel="bin", title=f"Correlation matrix {name}{title_suffix}", ) - pcolor_with_blocks(matrix_cor, blocks=elements, cmap=cmap) - if pdf: + hm = pcolor_with_blocks(matrix_cor, blocks=elements, cmap=cmap) + if pdf and hm is not None: pdf.savefig() plt.figure(figsize=figsize_2d) @@ -250,16 +326,16 @@ def main(opts: Namespace) -> None: 111, xlabel="", ylabel="bin", - title=f"Correlation matrix {name} ({elements[0]}){title_suffix}", + title=f"Correlation matrix {name} (block {block_name}){title_suffix}", ) hm = pcolor_with_blocks( - matrix_cor[:blocksize, :blocksize], + matrix_cor[block_selection, block_selection], blocks=elements[:1], pcolormesh=True, cmap=cmap, ) # heatmap_show_values(hm, lower_triangle=True) - if pdf: + if pdf and hm is not None: pdf.savefig() plt.figure(figsize=figsize_2d) @@ -276,7 +352,7 @@ def main(opts: Namespace) -> None: logger.info(f"Plot {name}") - if pdf: + if pdf and hm is not None: pdf.savefig() if opts.show: @@ -285,7 +361,7 @@ def main(opts: Namespace) -> None: if not opts.show: plt.close("all") - if pdf: + if pdf and hm is not None: pdf.__exit__(None, None, None) logger.info(f"Write output file: {ofile}") @@ -327,7 +403,7 @@ def covariance_get_matrices( def heatmap_show_values( - pc: "QuadMesh", fmt: str = "%.2f", lower_triangle: bool = False, **kwargs + pc: QuadMesh, fmt: str = "%.2f", lower_triangle: bool = False, **kwargs ): from numpy import mean, unravel_index @@ -364,9 +440,9 @@ def matrix_sum_blocks(matrix: NDArray, blocksize: int) -> NDArray: nblocks = nr // blocksize ret = empty((nblocks, nblocks), dtype=matrix.dtype) for row in range(nblocks): + row1 = row * blocksize + row2 = row1 + blocksize for col in range(nblocks): - row1 = row * blocksize - row2 = row1 + blocksize col1 = col * blocksize col2 = col1 + blocksize @@ -394,28 +470,41 @@ def pcolor_with_blocks( pcolormesh: bool = False, sep_kwargs: Mapping = {}, vmax: float | None = None, + symlog: bool = False, **kwargs, ): from numpy import fabs - from numpy.ma import array dmin = data.min() dmax = data.max() + if dmin == dmax == 0: + return None + bound = max(fabs(dmin), dmax) if vmax is not None: vmin, vmax = -vmax, vmax else: vmin, vmax = -bound, bound + if symlog: + pcolor_kwargs = dict( + norm=SymLogNorm( + linthresh=bound * 0.01, linscale=0.05, vmin=vmin, vmax=vmax + ), + **kwargs, + ) + else: + pcolor_kwargs = dict(vmin=vmin, vmax=vmax, **kwargs) + # fdata = fabs(data) # data = array(data, mask=(fdata < 1.0e-9)) ax = plt.gca() ax.set_aspect("equal") if pcolormesh: - hm = ax.pcolormesh(data, *args, vmin=vmin, vmax=vmax, **kwargs) + hm = ax.pcolormesh(data, *args, **pcolor_kwargs) else: - hm = ax.pcolorfast(data, *args, vmin=vmin, vmax=vmax, **kwargs) + hm = ax.pcolorfast(data, *args, **pcolor_kwargs) hm.set_rasterized(True) if colorbar: add_colorbar(hm, rasterized=True) @@ -436,9 +525,7 @@ def pcolor_with_blocks( def _get_blocks_data(size: int, blocks: Sequence[str]) -> NDArray: n_blocks = len(blocks) bins_in_block = size // n_blocks - xs = arange(0, n_blocks + 1) * bins_in_block - - return xs + return arange(0, n_blocks + 1) * bins_in_block def _plot_separators( diff --git a/scripts/run_dayabay_v0.py b/scripts/run_dayabay_v0.py index 05e07f805c961eb63340fec4830f82bfda7222f5..4223a7dff73791d41395f45b5760b0dd1b048608 100755 --- a/scripts/run_dayabay_v0.py +++ b/scripts/run_dayabay_v0.py @@ -119,15 +119,16 @@ def main(opts: Namespace) -> None: accept_index=accept_index, filter=accept_index, ) - storage["parameters.sigma"].savegraphs( - opts.graphs / "parameters" / "sigma", - mindepth=mindepth, - maxdepth=maxdepth, - keep_direction=True, - show="all", - accept_index=accept_index, - filter=accept_index, - ) + with suppress(KeyError): + storage["parameters.sigma"].savegraphs( + opts.graphs / "parameters" / "sigma", + mindepth=mindepth, + maxdepth=maxdepth, + keep_direction=True, + show="all", + accept_index=accept_index, + filter=accept_index, + ) storage["nodes"].savegraphs( opts.graphs, mindepth=mindepth, diff --git a/tests/models/test_dayabay_v0.py b/tests/models/test_dayabay_v0.py index da4683c5c27af0d2f2792b69a519e5de692a9e93..a02a950870c51a5c2e8b5833ed487de4a535829c 100644 --- a/tests/models/test_dayabay_v0.py +++ b/tests/models/test_dayabay_v0.py @@ -7,9 +7,6 @@ from models import available_models, load_model @mark.parametrize("modelname", available_models()) def test_dayabay_v0(modelname: str): - # TODO: remove when the model is done - if modelname=="v0d": - return model = load_model(modelname, close=True, strict=True) graph = model.graph