Skip to content

重新封装了 GMM 生成 ser/1.h5 的流程

lynn requested to merge fix-feature into master

暂时手扶了一下，把生成规则写在 Makefile 里留个记录，之后放到 production 仓库里的 Make 来调可能更好。

用 laser run 4056 刻度结果生成 ser，已经跑上了，明早能出来 1.h5，测试了一个子文件的结构没什么问题

with h5py.File("ser/0.h5", "r") as ipt:
    gain = ipt["gain"][:]
    ser = ipt["ser"][:]

gain
> array([( 1500, 5, [0.20387835, 0.43091337, 0.2276655 , 0.06557146, 0.07197132, 0.        ], [0.47503709, 0.91206184, 1.4828905 , 2.39981589, 4.20885362, 0.        ], [9.83889727e-02, 2.64838799e-01, 6.02569372e-01, 1.86670519e+00, 8.31123535e+00, 0.00000000e+00]),
       ( 2500, 5, [0.07651152, 0.52028664, 0.21107156, 0.15179516, 0.04033513, 0.        ], [0.33625932, 0.95616685, 1.65116006, 2.40332857, 3.38371998, 0.        ], [4.13855389e-01, 1.57090463e+00, 3.74975429e+00, 9.19732713e+00, 2.62184121e+01, 0.00000000e+00]),...
dtype=[('ch', '<u2'), ('nmean', '<u2'), ('ratio', '<f8', (6,)), ('mean', '<f8', (6,)), ('sig2', '<f8', (6,))])

ser
> array([( 1500, 1113.62219181, 868.18419774, 2.12841252e+05, [...], 7.1100527 ), ...
 dtype=[('ch', '<i8'), ('Gm', '<f8'), ('mus', '<f8'), ('sig2s', '<f8'), ('ser', '<f8', (77,)), ('sig2w', '<f8')])

对 MCP 的刻度结果做了 3，4，5，6 高斯分解然后选了 BIC 最小的那个。因为 laser run 的电荷数据量有限，刻度的结果除非对翁书记的模型再进行改进，不然暂时就是现阶段这样了，中位的卡方在 2.85 左右，均值 3.27，可以参考 https://juno.ihep.ac.cn/Dev_DocDB/0137/013708/002/20250417_mcp_progress.pdf 的附录。我明天也随机抽一些补一下高斯分解和 pdf 对照的图 @Berrysoft

/cc @xuchuang

Merge request reports