# 面向 BESIII 用户的 TF-PWA 使用手册 * 本手册仅面向使用 BESIII 集群的用户 * 原作者——蒋艺的 **Github** 地址：[https://github.com/jiangyi15/tf-pwa](https://github.com/jiangyi15/tf-pwa) * 我在 **Gitee** 克隆的地址：[https://gitee.com/wwdws1/tf-pwa](https://gitee.com/wwdws1/tf-pwa) * 我在所内 **GitLab** 克隆的地址(仅 **dev** 分支)：[https://code.ihep.ac.cn/wangshi/tf-pwa](https://code.ihep.ac.cn/wangshi/tf-pwa) * 原作者——蒋艺的完整的用户手册：[https://tf-pwa.readthedocs.io/en/latest/index.html](https://tf-pwa.readthedocs.io/en/latest/index.html) * 有任何问题和意见，欢迎给我发邮件：[wangshi@ihep.ac.cn](mailto:wangshi@ihep.ac.cn) * TF-PWA 作者蒋艺的邮箱：[jiangyi15@mails.ucas.ac.cn](mailto:jiangyi15@mails.ucas.ac.cn) # 目录 [TOC] # 1. 使用前的准备工作 ## 1.1 申请 GPU 账号 在已有集群账号的前提下，申请 Slum GPU 集群账号，详情见 Offline Software 组用户手册的3.2.2节：[链接](http://afsapply.ihep.ac.cn/cchelp/zh/local-cluster/jobs/slurm/#3222-slum-gpu%E9%9B%86%E7%BE%A4%E4%BD%BF%E7%94%A8%E6%96%B9%E6%B3%95)。分波分析申请的是 GPUPWA 组。 ## 1.2 准备末态粒子的四动量 * 数据文件每一行都是一个粒子的四动量 `E px py pz`，每m行为一组，表示m个最终粒子。 * 支持纯文本和numpy的数据文件，即 `.txt` 和 `.npy`。纯文本文件推荐精度为12位小数以上。 * 至少需要 data 和带效率的相空间 MC 的四动量文件。 # 2. 安装 ## 2.1 无安装使用方法 原作者已经将装好环境的 python 进行了打包，路径为： ``` /hpcfs/others/ucas/jiangy/python3.10.8-cp310-cp310-manylinux2014_x86_64.AppImage ``` 直接将其作为安装好的 python3 使用即可, 可以使用环境变量PYTHONPATH指定包的路径。 ## 2.2 使用 miniconda 和 pip 安装 ### 2.2.0 安装前的准备 把 `$HOME` 路径下的 `.cache` 和 `.local` 两个目录移动到其他的盘再[软链接](https://www.runoob.com/note/29134)回来，以免 `$HOME` 目录超出限额。 ### 2.2.1 在集群上安装 miniconda * ~~mambaforge 已被弃用~~。[miniconda 下载链接](https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh) * 下载完成之后，使用 bash 运行（如果你是 tcsh，那就换 bash），即： ``` bash Miniconda3-latest-Linux-x86_64.sh ``` 最开始是同意协议，yes 就行了。然后是指定安装位置，需要在 /hpcfs 盘，因为 gpu 作业节点没有接其他的盘。例（这里的 `<user>` 是让你填自己的用户名的）： ``` /hpcfs/bes/gpupwa/<user>/miniconda3 ``` 再往后就是安装过程，然后最后会问你是否初始化（`init`），选择初始化（`yes`），然后把初始化的那部分从 `.bashrc` 里拿出来单独放一个文件 `tf.sh`，每次需要这个环境的时候再 `source tf.sh`。因为这一步会比较卡（仅限在集群上，因为磁盘性能差），如果写在 `.bashrc` 里每次登录都会卡一会儿。然后新建一个新的环境： ``` conda create -n tf-pwa ``` 进入新环境： ``` conda activate tf-pwa ``` ### 2.2.2 下载源代码 找一个~~风水好的~~地方，然后 ``` git clone https://github.com/jiangyi15/tf-pwa.git ``` ~~然后这时候会有大概率失败，因为集群和 GitHub 的连接时有时无~~。现在集群已经稳定地完全无法连接 GitHub, 请务必使用在最开始写的几个克隆到国内的地址 ``` git clone https://code.ihep.ac.cn/wangshi/tf-pwa.git ``` 或者 ``` git clone https://gitee.com/wwdws1/tf-pwa.git ``` 更新方法：在没有作业正在跑的情况下，到克隆下来的 tf-pwa 的路径下，运行 ``` git pull ``` ### 2.2.3 安装 TensorFlow 你已经建好新的环境并且进了。官方教程在[这里](https://tensorflow.google.cn/install/pip#linux)。 **注：此节教程编写时间为2023年12月12日。** 1. 安装 python3.10 ``` conda install python=3.10 -c conda-forge ``` 2. 更新 pip（注：这里 `pip` 使用了清华的镜像加速下载） ``` python -m pip install -i https://pypi.tuna.tsinghua.edu.cn/simple --upgrade pip ``` 3. 安装 TensorFlow 和 CUDA ``` python -m pip install -i https://pypi.tuna.tsinghua.edu.cn/simple tensorflow[and-cuda]==2.14.1 ``` > tensorflow 2.15.0开始依赖cu12，但是集群的显卡驱动太老。最后一个使用cu11的是2.14.1 ### 2.2.4 安装其他依赖 然后安装 `uproot`，其作用为： > Uproot 是一个 Python 库，用于读取和处理 ROOT 文件格式的数据。 安装的命令为： ``` python -m pip install -i https://pypi.tuna.tsinghua.edu.cn/simple uproot ``` `cd` 到你克隆下来的 `tf-pwa` 里，运行如下命令安装依赖： ``` python -m pip install -i https://pypi.tuna.tsinghua.edu.cn/simple -e . ``` 对于此句的解释，原作者的解释为： > (Using `-e`, so it can be updated by `git pull` directly.) 到这里，所有的安装都已经结束，可以进行分波分析了。 # 3. 修改配置文件 配置文件的样板为 `tf-pwa` 路径下的 `config.sample.yml`。下面将分区域进行解释。特别注意，`yaml` 文件中的从属关系是通过缩进确定的。 ## 3.1 data ``` ## data entry is the data file name list used for fit data: ## 4-momentum order of final particles in dat files [option] ## the `finals` in `particle` will be used by default. ## it is necessary when dat files has momentum out of final particles dat_order: [B, C, D] ## basic data files, every line is `E px py pz` as 4-momentum of a particle, ## and every m lines group as m final particls data: ["data/data4600_new.dat"] ## data weight file, each line for a weight # data_weight: ["data/data_weight.dat"] ## phase space data for normalize amplitude phsp: ["data/PHSP4600_new.dat"] ## phase space weight file, each line for a weight # phsp_weight: ["data/phsp_weight.dat"] ## phase space without efficiencies, used for calculating fit fraction. # phsp_noeff: ["data/PHSP4600_noeff.dat"] ## much more phase space, used for plotting # phsp_plot: ["data/PHSP4600_plot.dat"] ## optional weight background data, need bg_weight bg: ["data/bg4600_new.dat"] ## background weight bg_weight: 0.731 ## inject MC in data # inmc: ["data/inMC.dat"] # inject_ratio: 0.01 # float_inmc_ratio_in_pdf: False ## using total monumtum direction as initial axis (or labrary axis) # random_z: True ## whether boost data to center mass frame first # center_mass: True # cached data file # cached_data: "data/all_data.npy" # data_charge: ["data/data4600_cc.dat"] # charge conjugation condition same as weight # cp_trans: True # when used charge conjugation as above, this do p -> -p for charge conjugation process. ``` 我们对该部分进行逐行讲解 ### 3.1.1 `dat_order: [B, C, D]` 这里指定的是你输入文件中对应的的粒子的顺序，有几个写几个，可以用更明显的名字代替，需要和decay和particle部分的配置一致。比如 $J/\psi\to\pi^+\pi^-\pi^0$，这里你可以写成 `dat_order: [pip, pim, piz]` 用以明显区分。 ### 3.1.2 `data: ["data/data4600_new.dat"]` 这里指定的是输入的 data 的四动量文件的路径，`.dat` 文件是纯文本文件，同样支持输入 `.npy` 文件。可以多个 config 文件使用同一个四动量文件用以节省磁盘。建议写成绝对路径避免出错。 ### 3.1.3 `phsp: ["data/PHSP4600_new.dat"]` 这里指定的是输入的相空间 MC 的四动量文件的路径。使用 truth 或 reconstructed 均可，但是需要带效率。 ### 3.1.4 `bg: ["data/bg4600_new.dat"]` (optional) 这里指定的是输入的本底的四动量文件的路径。目前只支持输入一个文件，所以有多个本底的需合成一个文件。**可以没有本底。** ### 3.1.5 `bg_weight: 0.731` (optional) 这里指定的是输入的本底的权重。如果输入本底的每个事例的 weight 一样，这里使用浮点数即可。如果每个事例有不同的 weight，则需要准备一个文件存入 weight，并与本底事例一一对应，并且输入的是一个文件。纯文本文件和 `.npy` 文件均可。 ### 3.1.6 `random_z: True` 该选项默认为 `True` 且注释，但是对于 BESIII 探测器的数据需要将该项取消注释并且修改为 `False`。（这是因为初态粒子直接从正负电子出来，极化的方向固定，如果不是正负电子直接产生的，不需要`False`） ### 3.1.7 `center_mass: True` 分波分析需要输入的四动量在母粒子静止系下，但是对于我们正常取得的数据则为实验室系下，存在一个沿着 $x$ 轴的分量。你可以选择输出四动量的时候就进行 boost，也可以输出实验室系下的四动量，在这里打开这个选项。建议取消注释。 ### 3.1.8 其他额外选项 TBD ## 3.2 decay ``` ## `decay` describe the decay structure, each node can be a list of particle name decay: ## each entry is a list of direct decay particle, or a list or the list. ## `dec_file` entry for loading `*.dec` files will be supported in the feature. ## `-` is often used for list in yaml A: ## for example A -> R_CD + B, R_CD -> C + D ## the optional `key: value` mean option parameters of the decay ## note there, R1 is overrated by `particle` entry as R1 and R11 - [R_CD, B, model: default] - - R_BD - C - model: "default" - - R_BC - D ## set l_list to allow only l=0 # - [ NR(BC)SP, D, curve_style: "g-", {has_barrier_factor: False, l_list: [0]}] # - [R_BC2, D, p_break: True] # allow P parity violate # R_BC2: [B, C, c_break: False] # enable C parity select C=(-1)^(l+s) R_CD: - C - D R_BD: [B, D] R_BC: [B, C] ## set l_list to allow only l=1 # NR(BC)SP: [B, C, {has_barrier_factor: False, l_list: [1]}] ``` ### 3.2.1 综述 在这一节内，你需要指定从母粒子经过多级衰变到你的输入末态的衰变链。并且这里仅支持两体衰变，相空间过程需要通过共振态模拟来实现。 ### 3.2.2 第一级衰变 以这里的写法，衰变为 $A\to B+C+D$，这里有三种可能的情况： 1. $A\to R\_BC+D, R\_BC\to B+C$ 2. $A\to R\_BD+C, R\_BD\to B+D$ 3. $A\to R\_CD+B, R\_CD\to C+D$ 这里的 $R\_BC$，$R\_BD$ 和 $R\_CD$ 指的是衰变到 $B+C$，$B+D$ 和 $C+D$的共振态（Resonance）。$R\_BC$，$R\_BD$ 和 $R\_CD$ 可以为各种粒子，这里只是每个衰变链上的总称，每个共振态的具体粒子在下面指定。 两种写法均可： ``` - [R_CD, B, model: default] ``` 和 ``` - - R_CD - B - model: "default" ``` 两种写法是一个意思。这里的 `default` 带不带双引号都可以，使用默认 model 的时候 `model: default` 可以不写。 ### 3.2.3 第二级衰变 由于这里的例子是三体末态，最多只有两级衰变。这里以 $R\_BC$ 的衰变举例：$R\_BC\to B+C$ 对应的写法为： ``` R_BC: [B, C] ``` 或 ``` R_BC: - B - C ``` 均可。同样的，使用默认的 model 的时候可以不写。 ### 3.2.4 其他选项 TBD ## 3.3 particle ``` ## `particle` describe the particle included ## top, finals, include are used as keywords particle: ## top decay particle $top: A: J: 1 # spin, use 1/2 or 0.5 for half spin P: -1 # parity, alias: `Par` mass: 4.6 # spins: [-1, 1] # list of possible of spin-projected quantum numbers # polarization: "vector" # possible polarization option ## final particles $finals: B: J: 1 P: -1 mass: 2.1 C: J: 1 P: -1 mass: 1.8 D: J: 0 P: -1 mass: 0.1 ## map for possible resonances in the same position ## each entry is a list of resonances id [option] ## include detailed parameters defined in `Resonances.yml` $include: Resonances.sample.yml R_BD: - D2_2460p - D1_2430p - D1_2420p R_BC: [Zc_4025, Zc_4160] ## set `[]` for particle which is not in real decays, ## just used for calculating angle data. R_CD: [] ## particle can be defined in here or included by `Resonances.yml` ## if both have the defination, defination in here will overtite defination in `Resonances.yml` # NR(BC)SP: {J: 1, P: -1, model: one} ``` ### 3.3.1 初、末态粒子定义 这里可以看到两类：`$top` 对应的是母粒子，`$final` 对应的是末态粒子。其中几个重要的参数为： 1. `J`：对应粒子的自旋，半整数可以使用分数（`1/2`）或者小数（`0.5`）来表示。 2. `P`：对应粒子的 P 宇称，正的就是 `+1`，负的就是 `-1`。 3. `mass`：对应粒子的质量，单位为 $\mathrm{GeV}/c^2$。 4. `spins` (optional)：对应粒子自旋可用的投影，对于 BESIII 实验，$J/\psi$ 粒子是通过虚光子而来，所以自旋的投影不能为 $0$。因此需要在这里指定其自旋的投影只有 $-1$ 和 $1$。 ### 3.3.2 中间共振态粒子选择 对于中间共振态，即 $R\_BC$，$R\_BD$ 和 $R\_CD$，他们所对应的粒子的定义在 `$include` 这一部分。这些中间共振态的具体定义可以写在这里，也可以引用写在 `Resonance.sample.yml` 中的粒子。这里有几个注意事项： 1. 如果在某一个质量谱上你没有放共振态，但是又想画它的质量分布，这里就需要写成 `[]`。例如： ``` R_BC: [] ``` 2. 对于在这里和 `Resonance.sample.yml` 中都定义过的粒子，这里的优先级是高于 `Resonance.sample.yml` 的。 ### 3.3.3 中间共振态粒子的定义 这里以 `Resonance.sample.yml` 中的粒子举例： ``` ## Resonance configure list ## yaml format version , see more: https://yaml.org/ ## required pyyaml module ## head resonance ## resonance name as id ## quotation marks are option ## head resonance will has head fixed to 1.0 D2_2460: ## spin [required] J: 2 ## parity (alias `Par`) [required] P: 1 ## resonance mass (alias `mass`) : m_0 [required] m0: 2.4607 m_max: 2.4611 m_min: 2.4603 ## resonance width (alias `width`): \Gamma_0 [required] g0: 0.0475 g_max: 0.0486 g_min: 0.0464 ## set mass and width (gamma0) as variables, and float in fit # float: mg ## add gaussian constraint term (e.g. -ln L -> -lm [L * Gauss(m0; mean, sigma)] # gauss_constr: ## the sigma of mass, mean is m0 # mass: 0.001 # width: 0.002 ## Breit Wigner formula [option] model: default ## whether used in decay disable: False Ds1_2700_2860: J: 1 P: -1 params: mass1: 2.7083 mass2: 2.859 width1: 0.12 width2: 0.159 mass1_range: [2.7, 2.72] mass1_sigma: 0.008 mass1_free: True mass1_constr: False m0: 2.83 model: Kmatrix ## other resonances D2_2460p: J: 2 P: 1 ## coef_head make D2_2460p to use the params of D2_2460 by the same name, (only works in simple cases) ## only the total phase angle different coef_head: "D2_2460" g0: 0.0467 g_max: 0.0479 g_min: 0.0455 m0: 2.4654 m_max: 2.4667 m_min: 2.4644 ## or you can used `{}` to avoid indent, use `,` to distinguish items. D1_2420: { J: 1, Par: 1, g0: 0.0317, g_max: 0.0342, g_min: 0.0292, m0: 2.4208, m_max: 2.4213, m_min: 2.4203, } D1_2420p: { J: 1, Par: 1, coef_head: "D1_2420", g0: 0.025, g_max: 0.021, g_min: 0.019, m0: 2.4232, m_max: 2.4256, m_min: 2.4208, } D1_2430: { J: 1, Par: 1, g0: 0.384, g_max: 0.514, g_min: 0.274, m0: 2.427, m_max: 2.467, m_min: 2.387, } D1_2430p: { J: 1, Par: 1, coef_head: "D1_2430", g0: 0.384, g_max: 0.514, g_min: 0.274, m0: 2.427, m_max: 2.467, m_min: 2.387, } # D0_2550: { # J: 0, Par: -1, # m0: 2.518, # g0: 0.199, # } # D0_2550p: { # J: 0, Par: -1, # coef_head: "D0_2550", # m0: 2.518, # g0: 0.199, # } # D1_2600: { # J: 1, Par: -1, # m0: 2.6419, # g0: 0.149, # } # D1_2600p: { # J: 1, Par: -1, # coef_head: "D1_2600", # m0: 2.6419, # g0: 0.149, # } Zc_4025: { J: 1, Par: 1, g0: 0.0248, m0: 4.0263 } Zc_4160: { J: 1, Par: 1, g0: 0.0921485, m0: 4.17329 } Zc_4210: J: 1 P: 1 m0: 4.21 g0: 0.11 ## trans model examples model: expr: a * BWR where: a: -1.0 Kmatrix_S: J: 1 P: 1 ## most model use m0 to calulate `q0` in barrier factor, m0: 4.21 ## KmatrixSingleChannel paramters model: KMatrixSingleChannel mass_list: [4.21, 4.31] width_list: [0.03, 0.06] ``` 这里我们分条讲解： 1. `J`：自旋。 2. `P`：P 宇称。 3. `m0`：质量初值，亦可写作 `mass`，单位 $\mathrm{GeV}/c^2$。 4. `m_max` 和 `m_min` (optional)：质量变化范围的上限和下限。 5. `g0`：宽度初值，亦可写作 `width`，单位 $\mathrm{GeV}$。 6. `g_max` 和 `g_min` (optional)：宽度变化范围的上限和下限。 7. `float: mg` (optional)：放开质量和宽度作为参数拟合，可以只写 `m` 或 `g`，用以只放开质量或宽度。 8. `gauss_constr` (optional)：高斯约束，指定的是质量或宽度的误差。 9. `coef_head` (optional)：固定参数与另一个共振态相同，让这两个共振态只有相角不同。 ## 3.4 约束 ``` ## constrains for params [WIP] constrains: # particle: # equal: # mass: [[D1_2430p, D1_2430]] ## fix the first decay chain total to 1.0 decay: { fix_chain_idx: 0, fix_chain_val: 1 } # fix_var: # "A->Zc_4025.DZc_4025->B.C_total_0r": 1.0 # var_range: # "A->Zc_4025.DZc_4025->B.C_total_0r": [1.0, null] # var_equal: # - ["A->D2_2460p.CD2_2460p->B.D_total_0r", "A->D2_2460.BD2_2460->C.D_total_0r"] ``` 1. 指定两个共振态的质量相同。 2. 指定第 0 个衰变链的 `total` 为 $1+0i$，这里是为了方便其他衰变链向其归一。 3. 固定某个参数为定值。 4. 指定某个参数的范围。 5. 指定某两个参数相同。 ## 3.5 画图 ``` ## plot describe the configuration of plotting 1-d data distribution plot: ## plot configuration TODO config: bins: 50 # legend_outside: True ## Make the legend outside the plot, default is False (inside) ## invariant mass mass: R_CD: ## label of display ## use `$$` for latex math display: "$M_{CD}$" # upper_ylim: 200 # trans: x*x R_BD: display: "$M_{BD}$" R_BC: display: "$M_{BC}$" ## helicity angle angle: ## If it named as `[.../]A/R_BC`, it refer to the decay A -> R_BC ... in the decay chain has A and R_BC, ## the angle is determined by the last 2 particles, the other order is random. ## for example, D->E+F in [A->R_BC+D, D->E+F, R_BC->B+C] is `R_BC/D/E` or `A/R_BC/D/E` or `R_BC/A/D/E` A/R_CD: ## include "alpha, cos(beta), gamma", # add `cos` for cos(beta) alpha: display: "$\\phi_{1}$" cos(beta): display: "$\\cos(\\theta_1)$" bins: 50 # number of bins is 50 legend: False # do not plot the legend range: [-1, 1] # range for plot A/R_BD: alpha: display: "$\\phi_{2}$" cos(beta): display: "$\\cos(\\theta_2)$" A/R_BC: alpha: display: "$\\phi_{3}$" cos(beta): display: "$\\cos(\\theta_3)$" R_BC/B: alpha: display: "$\\phi_{11}$" cos(beta): display: "$\\cos(\\theta_{11})$" R_CD/C: alpha: display: "$\\phi_{21}$" cos(beta): display: "$\\cos(\\theta_{21})$" R_BD/B: alpha: display: "$\\phi_{31}$" cos(beta): display: "$\\cos(\\theta_{31})$" ## 2D plot 2Dplot: m_R_CD & m_R_BC: display: "$M_{CD}$ vs $M_{BC}$" plot_figs: ["data", "sideband", "fitted"] dalitz_12: x: m_R_CD**2 y: m_R_BC**2 display: "$M_{CD}^2$ vs $M_{BC}^2$" plot_figs: ["data", "sideband", "fitted"] ``` ### 3.5.1 config 这里指定的参数适用于各个一维分布。优先级为：每个图单独的设置 > 总的设置 > 默认。常用参数如下： * `bins` (type=`int`)：设置分 bin 数，默认是 50。 * `legend` (type=`bool`)：是否画图例，质量谱默认画，角分布默认不画。 * `legend_outside` (type=`bool`)：画图例的话，是在图内还是在图外，默认是图内。示意图如下： | <center>  </center> | <center>  </center> | | -------- | -------- | | <center>`legend_outside=False` </center> | <center> `legend_outside=True` </center> | ### 3.5.2 mass 这部分是画各个质量谱的，即你前面定义个各种 `R_*`。选项有： * `display` (type=`"str"`)：图和 x 轴的 title。支持 latex 语法，需要用两个英文的 `$`，并且其中的反斜线命令需要两个反斜线。比如 $\Lambda$ 在这里需要写成 `$\\Lambda$`。 * `upper_ylim` (type=`int`)：y 轴的上限。 * `range` (type=`[double, double]`)：指定画图的范围，如果不指定则自动获取。 ### 3.5.3 angle 画角分布的时候比质量谱多了一点就是指定参考系。例如 `A/R_BD` 指在 $A$ 的静止系下画 $R\_BD$ 的角分布，这里的 alpha 和 beta 对应的是球坐标系下的 $\phi$ 角和 $\theta$ 角。还支持画角分布的余弦 `cos()`。角分布默认不画 legend。 ### 3.5.4 2Dplot 支持两种画图方式： * `x & y`：此种方法画的是散点图。x 和 y 是一维图的文件名，对于质量谱来说就是 `m_XXX`，角分布是 `XX_YY_ZZ`。如 `m_R_BC` 和 `A_R_BC_cos(beta)`。 * `XXXXX`：这里 XXXXX 是保存的文件名，下面的 x 和 y 是指定每一维画的什么。图的名字的规则和前面的一样，但是支持平方，即可以画 dalitz 图。 这里的 `plot_figs` 支持的如下： * `"data"`, `"sideband"`：对应 root 中二维图的 `SCAT`。 * `"data_hist"`, `"sideband_hist"`, `"fitted"`：对应 root 中二维图的 `COLZ`。 | <center>  </center> | <center>  </center> | | -------- | -------- | | <center> **data** </center> | <center> **data_hist** </center> | # 4. 提交作业 ## 4.1 准备作业运行脚本 这里使用计算中心提供的模板修改而来： ```#! /bin/bash ######## Part 1 ######### # Script parameters # ######################### # Specify the partition name from which resources will be allocated, mandatory option #SBATCH --partition=gpu # Specify the QOS, mandatory option #SBATCH --qos=pwanormal # Specify which group you belong to, mandatory option # This is for the accounting, so if you belong to many group, # write the experiment which will pay for your resource consumption #SBATCH --account=gpupwa # Specify your job name, optional option, but strongly recommand to specify some name #SBATCH --job-name=<name> # Specify how many cores you will need, default is one if not specified #SBATCH --ntasks=1 # Specify the output file path of your job # Attention!! Your afs account must have write access to the path # Or the job will be FAILED! #SBATCH --output=/hpcfs/bes/gpupwa/<user>/<dir>/<script_name>-%j.out # Specify memory to use, or slurm will allocate all available memory in MB #SBATCH --mem-per-cpu=16384 # Specify how many GPU cards to use #SBATCH --gres=gpu:v100:1 ######## Part 2 ###### # Script workload # ###################### # Replace the following lines with your real workload # list the allocated hosts srun -l hostname # list the GPU cards of the host /usr/bin/nvidia-smi -L python fit.py -c config.yml ``` 第一部分为作业配置，每一个选项分别对应： 1. `#SBATCH --partition=gpu`：资源分区。gpupwa 组可用的资源区有两个，分别为 `gpu` 和 `gpupwa`。 2. `#SBATCH --qos=pwanormal`：队列。gpupwa 组 gpu 分区可用的 QOS 有 `pwanormal` 和 `debug`，gpupwa 组 gpupwa 分区可用的 QOS 有 `pwadedicate` 和 `pwadebug`。具体区别详见[集群资源、队列和组别说明](http://afsapply.ihep.ac.cn/cchelp/zh/local-cluster/jobs/slurm/#%E9%9B%86%E7%BE%A4%E8%B5%84%E6%BA%90%E3%80%81%E9%98%9F%E5%88%97%E5%92%8C%E7%BB%84%E5%88%AB%E8%AF%B4%E6%98%8E)。 3. `#SBATCH --account=gpupwa`：组别。你就是 gpupwa 组的，这个不改。 4. `#SBATCH --job-name=<name>`：作业名。自己起个名改了就行了。 5. `#SBATCH --ntasks=1`：申请的 CPU 核心数。1 个够了，并且建议用 1 个。 6. `#SBATCH --output=/hpcfs/bes/gpupwa/<user>/<dir>/<script_name>-%j.out`：输出 log 文件绝对路径。这里的路径一定要是已经存在的，否则运行失败。其中的 `%j` 会被替换成作业的号。 7. `#SBATCH --mem-per-cpu=16384`：每个 CPU 核心数对应的内存大小，单位 MB。 8. `#SBATCH --gres=gpu:v100:1`：申请的 GPU 名称和数量。每个资源分区对应的 GPU 名称和数量详见[集群资源、队列和组别说明](http://afsapply.ihep.ac.cn/cchelp/zh/local-cluster/jobs/slurm/#%E9%9B%86%E7%BE%A4%E8%B5%84%E6%BA%90%E3%80%81%E9%98%9F%E5%88%97%E5%92%8C%E7%BB%84%E5%88%AB%E8%AF%B4%E6%98%8E)。 第二部分是运行区域，可以理解为 `bash` 脚本。 ## 4.2 提交作业命令 ``` # command to submit a job $ sbatch <job_script.sh> # <job_script.sh> is the name of the script,e.g: v100_test.sh, then the command is: $ sbatch v100_test.sh # There will be a jobid returned as a message if the job is submitted successfully ``` ## 4.3 查询作业 ``` # command to check the job queue $ squeue # command to check the jobs submitted by user $ sacct ``` ## 4.4 取消作业 ``` # command to cancel the job $ scancel <jobid> # <jobid> can be found using the command sacct ``` # 5. 查看拟合结果 ## 5.1 log 很多的信息都会保存在 log 中，可以自行查看。 ## 5.2 图 程序会根据 `config.yml` 里的设置进行画图，保存在提交作业所在目录下的 `figure` 目录内。 ## 5.3 参数 拟合中的各种参数，保存在 `final_params.json` 文件中。如果你放开了共振态的质量宽度，则结果也在里面。第一部分为中心值，第二部分为误差，第三部分为拟合状态，NLL(negative log likelihood) 和自由度。 ## 5.4 共振态占比 (fit fraction) 中心值保存在 `fit_frac.csv` 中。其中对角元为某个共振态自己的占比，非对角元为行列对应共振态之间的干涉。误差则保存在 `fit_frac_err.csv` 中。 # 6. 进阶 本节将按照不同的目标进行讲解。 ## 6.1 产生无效率的 PHSP 样本 ``` import tensorflow as tf from tf_pwa.config_loader import ConfigLoader def main(): config = ConfigLoader("/hpcfs/bes/mlgpu/wangshi/tfpwa/pksp/pwa/all/config.yml") N = 2000000 with tf.device("/device:GPU:0"): phsp = config.generate_phsp(N) config.data.savetxt("./data/phsp_noeff.npy", phsp) if __name__ == "__main__": main() ``` 其中 `config.yml`、事例数和 `.npy` 文件的名字和位置可以修改。 ## 6.2 基于分波结果产生 toy MC 样本 下面以 $\psi(3686)\to pK^0_S\bar{\Sigma}^-$ 为例说明如何产生 toy MC 样本 ``` import csv import json import time from pprint import pprint import matplotlib # avoid using Xwindow matplotlib.use("agg") import numpy as np import tensorflow as tf # examples of custom particle model from tf_pwa.amp import simple_resonance from tf_pwa.config_loader import ConfigLoader, MultiConfig from tf_pwa.experimental import extra_amp, extra_data from tf_pwa.utils import error_print, tuple_table devices = "/device:GPU:0" num_data = 10000000 final_paricle_num = 3 lb = ["pp", "pm"] resons = ["all", "R1", "R2", "R3", "R4", "R5"] particles_num = np.array([2212, 310, 3222], dtype=int) # p, K_S0, Sigma+ anti_params = np.array([1, 0, -1], dtype=int) # p, K_S0, anti-Sigma- config_file = "config.yml" def main(): for reson in resons: for i in range(len(lb)): final_paras = "final_params_" + str(reson) + ".json" data_dat = "data/data_" + str(lb[i]) + "_" + str(reson) + ".npy" pdata_dat = "data/pdata_" + str(lb[i]) + "_" + str(reson) + ".dat" config = ConfigLoader(config_file) print("Loaded parameter file: " + final_paras) config.set_params(final_paras) with tf.device(devices): data = config.generate_toy(num_data) config.data.savetxt(data_dat, data) print("Saved dat file: " + data_dat) p4 = np.load(data_dat) print("Loaded dat file: " + data_dat) f1 = open(pdata_dat, "w") for j in range(num_data): f1.write("%d\n" % final_paricle_num) for k in range(final_paricle_num): if anti_params[k] == 0: particle_num = int(particles_num[k]) else: particle_num = int( (-1) ** (i + 2) * anti_params[k] * particles_num[k] ) f1.write( str(particle_num) + " " + str("%.12f" % p4[j][k][1]) + " " + str("%.12f" % p4[j][k][2]) + " " + str("%.12f" % p4[j][k][3]) + " " + str("%.12f" % p4[j][k][0]) + "\n" ) f1.close() print("Wirted pdat file: " + pdata_dat) if __name__ == "__main__": main() ``` 这里是基于你拟合的 `config.yml` 文件，以及拟合结果的 `final_params.json` 文件，产生无效率的 toy MC truth，`num_data` 是你要产生的事例数。 修改好了之后交 GPU 作业，会产生 `.npy` 和 `.dat` 文件。其中 `.npy` 文件是纯四动量，三维矩阵：事例数，粒子数，$E,p_x,p_y,p_z$。注意这里产生的四动量是母粒子静止系下的，BesEvtGen 在模拟的时候会自动 boost。`.dat` 文件是用于 BesEvtGen 模拟是使用的，文件内的格式为： ``` ······ 3 2212 -0.657680440696 0.252044740774 0.505645335494 1.277537408514 -321 0.143385615443 -0.327185941150 0.480461257571 0.775996182438 -3212 0.514294825253 0.075141200376 -0.986106593065 1.632466391790 ······ ``` 这里具体的格式详见 `$BESEVTGENROOT/src/EvtGen/EvtGenModels/EvtTrackGen.cc` 文件开头的注释部分。对于如何使用四动量产生 toy MC，BesEvtGen 会在提供的四动量中抽样。boss 作业的修改方法如下： * 修改 decay card 产生相空间 MC 的是这样 ``` Decay psi(2S) 1.0000 p+ K- anti-Sigma0 PHSP; Enddecay ······ ``` 修改后是这样 ``` Decay psi(2S) 1.0000 p+ K- anti-Sigma0 TrackGen 0; Enddecay ······ ``` * 修改模拟的 jobOption 默认的里面关于 EvtGen 的是这样 ``` //*************job options for EvtGen*************** #include "$BESEVTGENROOT/share/BesEvtGen.txt" EvtDecay.userDecayTableName = "decay.dec"; ``` 修改之后是这样 ``` //*************job options for EvtGen*************** #include "$BESEVTGENROOT/share/BesEvtGen.txt" EvtDecay.userDecayTableName = "decay.dec"; EvtDecay.FileForTrackGen={"pdata.dat"}; ``` **注意：HTCondor 节点没有挂载 `/hpcfs`，所以四动量文件需要复制到非 `/hpcfs` 盘后方可调用。** ## 6.3 质量分辨 关于考虑质量分辨的原理，详见蒋艺的手册：[Resolution in Partial Wave Analysis](https://tf-pwa.readthedocs.io/en/latest/resolution.html)。 关于具体实现，这里分为几个部分： 1. 将 data 扩大 n 倍 2. 将 background 同样扩大 n 倍 3. 画图的时候，MC 需要带分辨的样本，即重建值 这里用到三个新的脚本产生上述 1 和 2 步的扩大后的样本： 1. `detector.yml`：这里面写的是你的分辨和偏移，单位 $\mathrm{GeV}/c^2$ ``` sigma: 0.001951045 bias: 0.0 ``` 2. `detector.py`：这个是产生新样本的时候调用的。这里需要修改三个 `Define` 的部分，分别是共振态的名字，相空间上限和下限的约束以及 `config.yml` 文件的路径。**注意**，其中trans_function(m1,m2)需要自己定义，这里是简单高斯加上一个线性效率以及范围截断。只考虑高斯分辨可以把eff处的改为 `eff=np.where((m2<mass_max_cut)&(m2>mass_min_cut),1.,0.)*phsp_factor(m1)`. 对于一般情况trans_function返回值是相空间MC中truth(m1)和rec(m2的二维分布概率密度。一般的做法是按照m2分bin，每个bin里得到一个m1的分布，这样总的分布就是m2本身的分布乘以每个bin里的m1的分布。 ``` """ Example of detector model """ import numpy as np import tensorflow as tf import yaml from tf_pwa.config_loader import ConfigLoader from tf_pwa.data import data_mask from tf_pwa.data_trans.helicity_angle import HelicityAngle # Define decay chain to be smeared resonance = "R_CD" # Define the mass of each particle, in GeV, order by decay chain mass_mother = 3.097 mass_one = 1.192642 mass_two = 0.134977 mass_another = 1.115683 mass_min_cut = mass_one + mass_two mass_max_cut = mass_mother - mass_another # Define the config.yml file config_file = "/hpcfs/bes/mlgpu/wangshi/tfpwa/jlsp/pwa/all/config.yml" config = ConfigLoader(config_file) decay_chain = config.get_decay(False).get_decay_chain(resonance) ha = HelicityAngle(decay_chain) # load detector model parameters with open("detector.yml") as f: detector_config = yaml.safe_load(f) def detector(data): """ run detector for data """ ms, costheta, phi = ha.find_variable(data) p_map = {str(i): i for i in ms.keys()} m = ms[p_map[resonance]] # smear mass new_m = ( m + tf.random.normal(shape=m.shape, dtype=m.dtype) * detector_config["sigma"] + detector_config["bias"] ) # selected with effecency, eff(m) = m # weight = m # max_w = 2.0 # cut = max_w * tf.random.uniform(shape=m.shape, dtype=m.dtype) < weight # selected in reconstructed value, a mass cut # cut = (cut) & (new_m < mass_max_cut) & (new_m > mass_min_cut) cut = (new_m < mass_max_cut) & (new_m > mass_min_cut) toy_truth = data_mask(data, cut) ms[p_map[resonance]] = new_m toy_rec = ha.build_data(*data_mask((ms, costheta, phi), cut)) toy_rec = config.data.cal_angle(toy_rec) return toy_truth, toy_rec def relative_p(m0, m1, m2): """breakup monmentum""" s2 = (m0**2 - (m1 + m2) ** 2) * (m0**2 - (m1 - m2) ** 2) s2 = np.where(s2 < 0, 0, s2) return np.sqrt(s2) / 2 / m0 def phsp_factor(m): """phase space factor""" p1 = relative_p(mass_mother, m, mass_another) # A -> R_CD + B p2 = relative_p(m, mass_one, mass_two) # R_CD -> C + D return p1 * p2 def gauss(delta, sigma): """simple gauss function without normalisation""" return np.exp(-(delta**2) / 2 / sigma**2) def trans_function(m1, m2): """ transisation function """ delta = m2 - m1 - detector_config["bias"] r1 = gauss(delta, sigma=detector_config["sigma"]) eff = np.where((m2 > mass_min_cut) & (m2 < mass_max_cut), 1.0, 0.0) * phsp_factor( m1 ) return r1 * eff def rec_function(m2): """pdf of reconstructed value""" ret = 0 N = 10000 m1 = np.linspace(mass_min_cut, mass_max_cut, N) ret = trans_function(m1[:, None], m2) return np.mean(ret, axis=0) def truth_function(m1): """pdf of truth value""" ret = 0 N = 10000 m2 = np.linspace(mass_min_cut, mass_max_cut, N) ret = trans_function(m1[:, None], m2) return np.mean(ret, axis=1) def delta_function(delta): """pdf of reconstructed value - truth value""" N = 10000 ms_min = mass_min_cut * 2 + np.abs(delta) # np.where(delta<0, 0.2+delta, 0.2-delta) ms_max = mass_max_cut * 2 - np.abs(delta) # np.where(delta<0, 1.9+delta, 1.9-delta) # print(delta, ms_min, ms_max) ms = ( np.linspace(0.0 + 1 / 2 / N, 1.0 - 1 / 2 / N, N)[:, None] * (ms_max - ms_min) + ms_min ) # print(ms, (ms-delta)/2, (ms+delta)/2) ret = trans_function((ms - delta) / 2, (ms + delta) / 2) return np.mean(ret, axis=0) * (ms_max - ms_min) def run(name): """run detector for data sample""" toy = config.data.load_data("data/" + name + ".dat") toy_truth, toy_rec = detector(toy) toy_rec.savetxt("data/" + name + "_rec.dat", config.get_dat_order()) toy_truth.savetxt("data/" + name + "_truth.dat", config.get_dat_order()) if __name__ == "__main__": # run("toy") run("phsp") # run("phsp_plot") ``` 3. `sample.py`：产生扩大后样本的脚本。其中需要修改的部分仍为三个 `Define` 的部分，分别为本底的 weight，共振态的名字（和 `detector.py` 中一致）以及 `config.yml` 文件的路径（和 `detector.py` 中一致）。这里是通过trans_function均匀取点，也可以改成不均匀的取点，提高精度。 ``` import numpy as np import tensorflow as tf import yaml from detector import trans_function from scipy.stats import norm from tf_pwa.config_loader import ConfigLoader from tf_pwa.data_trans.helicity_angle import ( HelicityAngle, HelicityAngle1, generate_p, ) # loda detector model parameters with open("detector.yml") as f: detector_config = yaml.safe_load(f) # Define the total weight of background bg_weight_total = 1.0 # Define decay chain to be smeared resonance = "R_CD" # Define the config.yml file config_file = "/hpcfs/bes/mlgpu/wangshi/tfpwa/jlsp/pwa/all/config.yml" def smear_function(m, m_min, m_max, i, N): """generate mass of truth sample and weight""" delta_min = m_min - m delta_max = m_max - m sigma = detector_config["sigma"] bias = detector_config["bias"] delta_min = np.max([delta_min, -5 * sigma * np.ones_like(delta_min) - bias], axis=0) delta_max = np.min([delta_max, 5 * sigma * np.ones_like(delta_max) - bias], axis=0) delta = ( np.linspace(1 / N / 2, 1 - 1 / N / 2, N)[i] * (delta_max - delta_min) + delta_min ) w = trans_function(m + delta, m) return m + delta, w def smear(toy, decay_chain, name, function, idx, N): """generate full truth sample and weight""" ha = HelicityAngle(decay_chain) ms, costheta, phi = ha.find_variable(toy) m_min = 0 m_max = np.inf par = None for i in decay_chain: if str(i.core) == str(name): m_min = ms[i.outs[0]] + ms[i.outs[1]] par = i.core if str(name) in [str(j) for j in i.outs]: m_max = ms[i.core] for j in i.outs: if str(j) == str(name): continue m_max = m_max - ms[j] m_smear, w = function(ms[par], m_min, m_max, idx, N) # print("smear", idx, m_smear, ms[par], m_smear-ms[par], w) # exit() ms[par] = m_smear return ha, ha.build_data(ms, costheta, phi), w def random_sample(config, decay_chain, toy): """generate total truth sample and weight""" all_p4 = [] ws = [] for i in range(config.resolution_size): # decay_chain = [i for i in toy["decay"].keys() if "(p+, p-, pi+, pi-)" in str(i)][0] # config.get_decay(False).get_decay_chain("(B, C)") ha, toy_smear, w = smear( toy, decay_chain, resonance, smear_function, i, config.resolution_size ) ws.append(w) # print(toy_smear) # toy_smear = {i: j for i,j in zip(ha.par, toy_smear)} # print(toy_smear) p4 = [ toy_smear[i].numpy() for i in config.get_dat_order() ] # (particle, idx, mu) all_p4.append(np.stack(p4)) # (resolution, particle, idx, mu) ws = np.stack(ws) ws = ws / np.sum(ws, axis=0) return np.stack(all_p4).transpose((2, 0, 1, 3)), ws def random_sample_bg(config, decay_chain, toy): """generate total truth sample and weight""" all_p4 = [] ws = [] for i in range(config.resolution_size): # decay_chain = [i for i in toy["decay"].keys() if "(p+, p-, pi+, pi-)" in str(i)][0] # config.get_decay(False).get_decay_chain("(B, C)") ha, toy_smear, w = smear( toy, decay_chain, resonance, smear_function, i, config.resolution_size ) ws.append(w) # print(toy_smear) # toy_smear = {i: j for i,j in zip(ha.par, toy_smear)} # print(toy_smear) p4 = [ toy_smear[i].numpy() for i in config.get_dat_order() ] # (particle, idx, mu) all_p4.append(np.stack(p4)) # (resolution, particle, idx, mu) ws = np.stack(ws) ws = ws / np.sum(ws, axis=0) * bg_weight_total return np.stack(all_p4).transpose((2, 0, 1, 3)), ws def main(): config = ConfigLoader(config_file) decay_chain = config.get_decay(False).get_decay_chain(resonance) toy = config.get_data("data_rec")[0] # print(toy) # exit() # ha = HelicityAngle(decay_chain) # ms, costheta, phi = ha.find_variable(toy) # dat = ha.build_data(ms, costheta, phi) # print(var) # print(ha.build_data(*var)) p4, w = random_sample(config, decay_chain, toy) # exit() np.save("data/data.npy", np.stack(p4).reshape((-1, 4))) np.savetxt( "data/data_w.dat", np.transpose(w).reshape((-1,)) ) # np.ones(p4.reshape((-1,4)).shape[0])) toy = config.get_data("bg_rec")[0] p4, w = random_sample_bg(config, decay_chain, toy) np.save("data/bg.npy", np.stack(p4).reshape((-1, 4))) np.savetxt("data/bg_w.dat", np.transpose(w).reshape((-1,))) # toy = config.get_data("phsp_rec")[0] # p4, w = random_sample(config, decay_chain, toy) # np.savetxt("data/phsp_re.dat", p4.reshape((-1,4))) # np.savetxt("data/phsp_re_w.dat", np.transpose(w * toy.get_weight()).reshape((-1,))) # np.repeat(toy.get_weight(), config.resolution_size)) if __name__ == "__main__": main() ``` 然后需要修改 `config.yml` 文件中的 `data` 部分： ``` ······ # data entry is the data file name list used for fit data: dat_order: [B, C, D] resolution_size: 25 data: ["resolution/data/data.npy"] data_weight: ["resolution/data/data_w.dat"] data_rec: ["data.dat"] phsp: ["phsp_truth.dat"] phsp_rec: ["phsp_rec.dat"] phsp_noeff: ["resolution/data/phsp_noeff.npy"] bg: ["resolution/data/bg.npy"] bg_weight: ["resolution/data/bg_w.dat"] bg_rec: ["bg.dat"] bg_rec_weight: 1.0 ······ ``` 这里**多余的部分和注释就不写了**，只需要注意我下面列出的几条： * `resolution_size`：扩大的倍数，请量力而为。 * `data` 和 `data_weight`：这里放的是扩大后的 data 和其一一对应的 weight。把 data 存成 `.npy` 是为了节省空间，下同。 * `data_rec`：这是扩大前，原本的 data，这里用来画图。 * `phsp`，`phsp_rec` 和 `phsp_noeff`：第一个是在拟合中做 MC 积分的，第二个是画图的时候为了和有质量分辨的 data 对应的有分辨的 MC 样本，第三个是用来计算 fit ftaction 时用的无效率的 PHSP MC 样本。 * `bg` 和 `bg_weight`：扩大后的本底和其一一对应的 weight。 * `bg_rec` 和 `bg_rec_weight`：这是扩大前的本底事例，以及其归一化因子。注意这里的归一化因子应该和 `sample.py` 中要修改的权重一致。 搞好这四个文件之后，把 `sample.py` 交到 GPU 作业上，等跑完就将数据扩完了。但是这里需要检查两个 `weight` 文件中是否有 nan。如果有，说明这个事例由于分辨原因超出了相空间允许的范围，需要从原文件中删除。这里有一个写好的检查并删除对应事例的脚本 `re_cal_data.py`： ``` # Read the .txt file to get the line number, then print the line number per 10 lines import numpy as np line_number_file = "resolution/data/nan_check.txt" four_momentum_file = "data.dat" save_four_momentum_file = "resolution/data/data_new.dat" def read_line_number(file_name): line_number = np.loadtxt(file_name) len_line_number = line_number.shape[0] return line_number, len_line_number def read_four_momentum(file_name): four_momentum = np.loadtxt(file_name) return four_momentum def save_four_momentum(file_name, four_momentum): np.savetxt(file_name, four_momentum, delimiter=" ", fmt="%.12f") def main(): line_number, len_line_number = read_line_number(line_number_file) four_momentum = read_four_momentum(four_momentum_file) for i in range(int(len_line_number / 10) - 1, -1, -1): print(int(line_number[i * 10 + 9] / 10)) four_momentum = np.delete( four_momentum, (int(line_number[i * 10 + 9] / 10) * 3), 0 ) four_momentum = np.delete( four_momentum, (int(line_number[i * 10 + 9] / 10) * 3 - 1), 0 ) four_momentum = np.delete( four_momentum, (int(line_number[i * 10 + 9] / 10) * 3 - 2), 0 ) save_four_momentum(save_four_momentum_file, four_momentum) if __name__ == "__main__": main() ``` 用法为： * 先产生10倍扩大的样本 * 使用 `grep -n "nan" resolution/data/data_w.dat >resolution/data/nan_check.txt` 将 nan 对应的行数输出到 `nan_check.txt` 中，然后将该文件中的除数字外的所有字符删除 * 运行 `re_cal_data.py`，将新产生的 data 覆盖原有的文件 * 再跑一遍 `sample.py` 确保其正常工作。 * bg 和 data 一样。 当所有的都确定没有问题之后，就可以基于改好的 `config.yml` 跑拟合了。 ## 6.4 大样本优化选项 ### 6.4.1 第一组 ``` ······ # data entry is the data file name list used for fit data: ······ use_tf_function: True jit_compile: True ······ ``` TensorFlow提供的编译选项，将动态图转化为静态图。 不支持TensorFlow以外的运算（会固定到第一次的值）。会造成第一步编译耗费大量时间，在不拟合只是计算一下结果时（或者迭代次数较少时）建议关闭。 * `use_tf_function`：tensorflow 优化项 * `jit_compile`：实时（just in time）编译，，tensorflow 优化项。会需要额外的内存。 * **注意**: 1. 如果是使用 `conda` 安装的，会有大概率启用此选项后无法正常拟合，则只需注释这一行即可。 2. 目前集群上的显卡的驱动所对应的 CUDA 版本低于要求，所以使用最新方法安装的 TensorFlow 会有问题，请勿使用该选项。 ### 6.4.2 第二组 ``` ······ # data entry is the data file name list used for fit data: ······ lazy_call: True no_id_cached: True cached_lazy_call: "cached_data/" ······ ``` 这一组一般是一起使用的，其中每一项为： * `lazy_call`：启用惰性运行。主要针对角度计算，启用后会在每次计算时计算角度，避免角度占用大量内存。 * `no_id_cached`：在使用`lazy_call: True`时使用`use_tf_function: True`需要设为`True`。 * `cached_lazy_call`：指定缓存文件存储位置，如果指定位置为空（`""`），则缓存文件会存储在内存中。不指定则不会缓存，而是每次都重新计算。**注意**，每次修改 `config.yml` 文件(以及更新代码）后需要删除旧的缓存文件重新生成，未改变则无需重新生成。 ### 6.4.3 第三组 ``` ······ # data entry is the data file name list used for fit data: ······ amp_model: cached_amp preprocessor: cached_amp ······ ``` 新的 `amp_model` 和 `preprocessor`。改变amplitude不同部分的计算顺序，以及需要缓存的数据内容。支持的选项有`default`,`cached_amp`,`cached_shape`,`p4_directly`. 不同选项的具体表现取决于具体的过程。一般当初末态粒子的自旋皆为0时，cached_amp或cached_shape更快。对于初末态自旋较大时，p4_directly更快。 | | preprocessor | amp_model | | ---------- | ------------------ | --------------- | | default | p4 -> m + angle | $a_i f_i(m) D_i(angle)$ | |cached_amp | p4 -> m + angle + $D_i$ | $a_i f_i(m) D_i$ | |cached_shape| p4 -> m + angle + $f_i D_i + D_j$ | $a_i f_i D_i + a_j f_j(m) D_j$ | |p4_directly | p4 | $a_i f_i(m(p4))D_i(angle(p4))$ | $a_i$ 是拟合系数， $f_i(m)$ 是质量依赖部分， $D_i(angle)$ 是角度部分。 存储需求如下表（参考值：$10^7 \times 16\approx1.2$G） | | | | -- | -- | | p4 | 4 N(event) N(particles) | | m | N(event) N(particles) | | angle | N(event) N(topo chains)(9 N(final particles)-6) | | $f_iD_i, D_i$ | N(event) N(helicity combination) N(partial waves) | `N(topo chains)`取决与`config.yml`的`decay`里定义的数目，去掉没用到的decay chain，可以减少角度计算和缓存需求。 **注意** 画图是用的preprocessor 的结果，没有的将不能画图（比如p4_directly只有p4,没法画m和angle)。 ### 6.4.4 第四组 ``` ······ # data entry is the data file name list used for fit data: ······ no_left_angle: True no_angle: True no_p4: True ······ ``` 可以去掉部分不需要的角度，缓解存储压力。 * `no_left_angle`: 去掉反向粒子的helicity角。 * `no_angle`: 不画角分布。 * `no_p4`: 去掉原始的四动量。 ### 6.4.5 第五组 ``` ······ # data entry is the data file name list used for fit data: ······ lazy_call: True lazy_file: True # bool lazy_prefetch: 2 # int , -1 for auto, 0 for no, >0 for exact value ······ ``` 这些选项是为了超级大的 PHSP 样本，例如一亿的 PHSP。其中： * `lazy_file`：读取 `.npy` 文件会与内存对应，但是非常大的 `.npy` 会占用大量内存，使用memmap直接把读取内存映射到文件，减小内存使用。并且其需要与 `lazy_call` 一起使用。 * `lazy_prefetch`：prefetch会在GPU计算的同时，预先读取下面几步需要的数据，会消耗额外的内存。可以手动设置读取次数。 上述选项可以组合使用。 ## 6.5 其他脚本 ### 6.5.1 提前计算缓存文件 在使用 `cached_data` 时，可以基于 `config.yml` 文件提前计算缓存文件，提前释放内存，避免拟合时超限。路径：[create_lazy_cached_data.py](https://github.com/jiangyi15/tf-pwa/blob/dev/tutorials/examples/create_lazy_cached_data.py)