揭示振幅分析的秘密

目录

• 简介
• 安装
• 使用
• 待办事项

简介

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

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

现在我解释了为什么我不喜欢一些现有的解决方案，以及为什么我选择使用Rust，但这个项目有什么使其脱颖而出的地方？以下是一些吸引你的原因

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

安装

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

cargo add rustitude

Python 安装同样简单

pip install rustitude

使用

有关在 Rust 中编写自定义振幅的更深入教程，请参阅 rustitude-core crate。这个包主要关注 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' 的值将使用所有可用核心。