揭秘振幅分析

目录

• 介绍
• 安装
• 使用方法
• 待办事项

介绍

本项目始于希望创建一个快速且易于使用的界面，用于将振幅拟合到粒子物理数据中。话虽如此，也有一些性能良好的方法，如用C++编写的 AmpTools，但在我的个人经验中，它可能有点难以使用和扩展，并且通常需要大量模板代码来生成新的振幅或绘图脚本。另一方面，还有像 PyPWA（用Python编写）这样的库，它们看起来可能易于使用，但通常由于Python的限制性语法、速度问题和普遍缺乏文档以及持续开发，在这方面常常失败。已经尝试在AmpTools和Python之间架起桥梁，最近（并且成功）的是 PyAmpTools。这种方法的问题在于它依赖于PyROOT，这意味着您还需要安装ROOT（并且使用您版本的Python构建）。到目前为止，我将省略对ROOT的反感，只说ROOT应该是一个可选的，而不是必需的。那么，这又让 rustitude 处于什么位置呢？

正如其名所示，rustitude是用Rust编写的，所以让我们先谈谈明显的缺点：知道如何编写Rust代码的粒子物理学家并不多。希望这种情况会在接下来的十年内有所改变（而且美国政府在所有地方都提供了一些支持）。虽然与C++相比，Rust在知名度上相对较低，但它也有很多好处。没有null，就没有空引用（Tony Hoare的"十亿美元的错误"）。在Rust中，引用总是有效的，这是由一个非常有帮助且偶尔令人沮丧的借用检查器保证的。Rust "crate"的设置鼓励文档编写（参见rustitude-core的文档），而且对于使用Python的可选类型检查的人来说，其基本语法相当容易学习。Rust最大的好处之一可能就是它很容易实现并行化，但我最喜欢它的两个原因之一是它非常容易编写Python绑定（毕竟，这就是这个库的目的所在），并且它有一个包管理器。这一点非常重要——与C/C++不同，开发者会被一些Makefile、CMakeLists.txt或一些可能只在"X"系统上工作且只有在安装和使用make、cmake、g++等工具的情况下才能工作的怪物项目所淹没，而Rust只需在Cargo.toml中添加一行即可添加一个包（或者使用cargo add命令）。在许多方面，Rust的包管理实际上比Python更简单，因为只有一个首选的方法来创建和管理项目、格式化、代码审查和编译。

现在我已经解释了为什么我不喜欢一些现有的解决方案，以及为什么我选择使用Rust，但这个项目有什么特别之处？以下是一些吸引你的原因：

• rustitude将自动将振幅并行化到数据集中的事件上。开发者没有必要自己编写并行化代码。
• 在结构体上实现Node即可将其用作振幅。这意味着开发者只需编写两个到三个方法来描述振幅的全部功能，其中一个方法只是给出了振幅输入参数的名称和顺序。
• rustitude的一个主要目标是牺牲内存来提高处理速度。这是通过预先计算那些当自由参数输入改变时不会改变的振幅部分来实现的。《AmpTools》支持这种版本，但只是在每个通用振幅的层面上，而不是在个体层面上。最简单的例子是Ylm振幅（球谐函数），给定l和m的值，可以完全预先计算。在《AmpTools》中，不同实例的Ylm与不同的l和m共享预先计算的数据，而在《rustitude》中则不是。在《AmpTools》中实现的Ylm需要在每次函数调用时计算每个事件的球谐函数，而《rustitude》只需在数组中查找值即可！
• 库的大多数部分（公共接口）都有 Python 绑定，因此如果没有自定义振幅的需求，开发者实际上永远不需要编写任何 Rust 代码，生成的计算性能将与用 Rust 编写的性能相当。

安装

将 rustitude 添加到现有的 Rust 项目就像

cargo add rustitude

Python 安装同样简单

pip install rustitude

使用方法

有关在 Rust 中编写自定义振幅的更深入教程，请参阅 rustitude-core 包。此包主要关注 Python 方面。以下是示例分析的设置：

import rustitude as rt
from rustitude import gluex
import numpy as np

# Start by creating some amplitudes:
f0p = gluex.resonances.KMatrixF0('f0+', channel=2)
f0n = gluex.resonances.KMatrixF0('f0-', channel=2)
f2 = gluex.resonances.KMatrixF2('f2', channel=2)
a0p = gluex.resonances.KMatrixA0('a0+', channel=1)
a0n = gluex.resonances.KMatrixA0('a0-', channel=1)
a2 = gluex.resonances.KMatrixA2('a2', channel=1)
s0p = gluex.harmonics.Zlm('Z00+', l=0, m=0, reflectivity='+')
s0n = gluex.harmonics.Zlm('Z00-', 0, 0, '-')
d2p = gluex.harmonics.Zlm('Z22+', 2, 2) # positive reflectivity is the default

# Next, let's put them together into a model
# The API supports addition and multiplication and has additional methods for the real part (`real`) and imaginary part (`imag`).
pos_re_sum = (f0p + a0p) * s0p.real() + (f2 + a2) * d2p.real()
pos_im_sum = (f0p + a0p) * s0p.imag() + (f2 + a2) * d2p.imag()
neg_re_sum = (f0n + a0n) * s0n.real()
neg_im_sum = (f0n + a0n) * s0n.imag()

mod = rt.Model([pos_re_sum, pos_im_sum, neg_re_sum, neg_im_sum])

# There is no need to constrain amplitudes, since each named amplitude is only ever evaluated once and a cached value gets pulled if we run across an amplitude by the same name!
# We should, however, fix some of the values to make the fit less ambiguous. For instance, suppose we are above the threshold for the f_0(500) which is included in the F0 K-Matrix:
mod.fix("f0+", "f0_500 re", 0.0)
mod.fix("f0+", "f0_500 im", 0.0)
mod.fix("f0-", "f0_500 re", 0.0)
mod.fix("f0-", "f0_500 im", 0.0)

# As mentioned, we should also fix at least one complex phase per coherent sum:
mod.fix("f0+", "f0_980 im", 0.0)
mod.fix("f0-", "f0_980 im", 0.0)

# All done! Now let's load our model into a Manager, which helps coordinate the model with datasets.
# First, load up some datasets. rustitude provides an open operation that uses uproot under the hood:
ds = rt.open('data.root')
ds_mc = rt.open('mc.root')

m_data = rt.Manager(mod, ds)
m_mc = rt.Manager(mod, ds_mc)

# We could stop here and evaluate just one dataset at a time:

# res = m_data([1.0, 3.4, 2.8, -3.6, ... ])
# -> [5.3, 0.2, ...], a list of intensities from the amplitude calculation

# Or, what we probably want to do is find the negative log-likelihood for some choice of parameters:

nll = rt.ExtendedLogLikelihood(m_data, m_mc)

res = nll([10.0] * mod.n_free) # automatic CPU parallelism without GIL
print(res) # prints some value for the NLL

mi = rt.minimizer(nll, method='Minuit') # use iminuit to create a Minuit minimizer
mi.migrad() # run the Migrad algorithm
print(mi) # print the fit result

可以通过支持该功能的函数调用禁用 CPU 上的自动并行处理（例如，nll([10.0] * mod.n_free, parallel=False) 将在没有并行处理的情况下运行），并且可以通过 RAYON_NUM_THREADS 环境变量来控制使用的 CPU 数量，可以在代码运行之前设置或修改代码中的值（例如，os.environ['RAYON_NUM_THREADS] = '5' 将确保代码中只使用五个线程）。默认情况下，未设置值或 '0' 的值将使用所有可用核心。