ThOr 函子

摘要

• 在基础 DaVinci 脚本中引入更复杂的函子

• 理解这些复杂函子的基本原理

• 掌握使用函子进行基本计算的方法

最小化脚本

在上一节中，我们学习了 lbexec 的用法以及如何运行最小化的 DaVinci作业。这需要一个包含 main 函数的 Python 文件和一个用于配置参数的 YAML文件。下面是我们将在本节中构建的最小化脚本示例。您可以直接尝试此示例：将Python 和 YAML 文件复制到 lxplus 的目录中，并运行 Python脚本顶部的命令。

此脚本不再只是将衰变树输出到终端，而是生成包含实际物理测量结果的 ROOT格式N元组（ntuple）,这正是物理导论中所承诺的物理测量!

我们将使用的触发线（trigger line）同时选择\(B_s^0\to (D_s^- \to K^- K^+ \pi^-) \pi^+\)和\(B^0\to (D_s^- \to K^- K^+ \pi^-) \pi^+\)。因此需在 DaVinci 文件中包含衰变描述符，且必须与编写触发线时使用的描述符匹配。

DST文件包含海量信息：除了大量事件（其中仅有部分符合我们的衰变候选），还包含采集（或模拟）数据的丰富物理信息。为将这些信息写入ROOT文件，我们使用称为函子（Functors）的对象。函子还能基于顶点拟合结果计算衍生信息（例如粒子寿命）。

下面的 DaVinci脚本将遍历模拟的\(B_s^0\to (D_s^- \to K^- K^+ \pi^-) \pi^+\)的DST，从子粒子及其各自的母粒子中选择信息（具体来说，是质量、总动量和横向动量），并使用FunTuple 类将一个名为Hlt2B2OC_BdToDsmPi_DsmToKpKmPim/DecayTree的TTree 写入一个名为 tuple.root 的 ROOT 文件中。请看：

运行此脚本

请记住，在命令中更改你的 python 文件的名称（此处为DV_basic.py），同样也要更改你的 yaml 文件的名称（此处为tuple_options.yaml）

detdesc 与 DD4hep

运行 MC 时，必须通过 -c 参数指定包含 +detdesc 的平台。这是因为 Run 3 中探测器的几何结构和数据采集条件存储在两种不同系统中：

• detdesc：传统方法，仍用于 MC 数据（模拟组正在更新迁移）
• DD4hep：现代方法，用于真实数据，若未指定平台则为 LHCb 软件默认选项

'''
lb-run -c x86_64_v3-el9-gcc13+detdesc-opt+g DaVinci/v65r0 lbexec DV_basic:main tuple_options.yaml
'''
import Functors as F
import FunTuple.functorcollections as FC
from DaVinci import Options, make_config
from DaVinci.algorithms import create_lines_filter
from FunTuple import FunctorCollection
from FunTuple import FunTuple_Particles as Funtuple
from PyConf.reading import get_particles, get_pvs
from RecoConf.event_filters import require_pvs


def main(options: Options):
    line = "Hlt2B2OC_BdToDsmPi_DsmToKpKmPim"
    data = get_particles(f"/Event/HLT2/{line}/Particles")
    line_prefilter = create_lines_filter(name=f"PreFilter_{line}", lines=[line])
    pvs = get_pvs()

    fields = {
        "Bs": "[B0 -> (D_s- -> K- K+ pi-) pi+]CC", #Note here, we call the parent a B0 even though it's a Bs, because it needs to match the reconstruction in the trigger line
        "Ds": "[B0 -> ^(D_s- -> K- K+ pi-) pi+]CC",
        "Bs_pi": "[B0 -> (D_s- -> K- K+ pi-) ^pi+]CC",
        "Ds_pi": "[B0 -> (D_s- -> K- K+ ^pi-) pi+]CC",
        "Km": "[B0 -> (D_s- -> ^K- K+ pi-) pi+]CC",
        "Kp": "[B0 -> (D_s- -> K- ^K+ pi-) pi+]CC",
    }

    all_vars = FunctorCollection({
        "M": F.MASS,
        "P": F.P,
        "PT": F.PT
    })

    variables = {"ALL": all_vars}

    funtuple = Funtuple(
        name=line,
        tuple_name="DecayTree",
        fields=fields,
        variables=variables,
        inputs=data,
    )

    algs = {line: [line_prefilter, require_pvs(pvs), funtuple]}

    return make_config(options, algs)

input_files:
- root://eoslhcb.cern.ch//eos/lhcb/wg/dpa/wp7/Run3SK/exampleDST/00257716_00000011_1.hlt2.dst
input_type: ROOT
output_type: ROOT
input_raw_format: 0.5
simulation: true
input_process: Hlt2
input_stream: b2oc
conddb_tag: sim10-2024.W37.39-v00.00-md100
dddb_tag: dddb-20240427
lumi: false
data_type: '2024'
evt_max: 1000
print_freq: 100
ntuple_file: tuple.root

如先前所述，现在我们将重点介绍分析中可能用到的更多实用函子

关于 ThOr 函子的注意点

在旧版文档中您可能会看到 LoKi 函子，这是 DaVinci Run 1 和 Run 2 时期使用的旧版函子系统。如今我们已升级为 ThOr 函子（Throughput Oriented functors，面向吞吐量优化的函子）。其详细工作原理可参阅官方文档。

函子

这里有一份所有已记录函子的详细列表，还有一个非常棒的术语表资源。下面介绍一些在LHCb 分析中常用且至关重要的函子。

冲击参数（Impact Parameter , IP）

IP 是一条径迹与某个顶点（如主顶点 PV）之间的最短距离。

冲击参数卡方（IPCHI2）

它定量衡量一条径迹来自该顶点的可能性，同时考虑了 IP测量的不确定性。例如，IPCHI2 = 0 表示该径迹与顶点完全兼容；IPCHI2越大，则径迹与顶点之间的不一致性越强。

IPCHI2 的说明细节

IPCHI2 是包含和不包含所考虑粒子时顶点拟合的\(\chi^2\)增量。它的表现近似于撞击参数（IP）除以其不确定性的平方，但更加稳定。考虑一条被纳入初级顶点拟合的轨迹的撞击参数及其不确定性，此时撞击参数与初级顶点（PV）位置之间存在反相关性，轨迹会将初级顶点拉向自身。再看一个未被纳入拟合的粒子（可能因为距离较远，或者因为它是像B 粒子这样的复合粒子），这种情况下就不存在这种相关性。因此，"撞击参数/不确定性" 这一量的含义会因轨迹的使用方式不同而有所差异。而 IPCHI2则不存在这个问题：它始终是在比较包含和不包含该粒子时的初级顶点。

最近距离（DOCA）

它表示一对轨迹之间的最短距离。例如，在我们研究的\(B_{s}^0 \to (D_{s}^- \to K^- K^+ \pi^-) \pi^+\)衰变中，我们可能希望知道构成候选粒子\(D_{s}\)的子粒子之间的最近接近距离，或者其他粒子\(\pi\)之间的最近接近距离。理论上，若它们来自同一个衰变顶点，这个距离应为0，但由于外推误差、分辨率和噪声的影响，可能会出现较小的最近接近距离值。

方向角（DIRA）

方向角是从初级顶点（PV）到衰变顶点的连线与衰变产物的总 4 动量方向之间的夹角。函子 F.DIRA 实际返回的是方向角的余弦值（cos(DIRA)），因为当 cos (DIRA)≈1时，动量方向更接近指向初级顶点，这意味着该候选粒子更有可能是真实的 B 粒子。

粒子鉴别变量（PID Variables）

我们主要使用两类粒子鉴别变量，其中较基础的函子如下：

• F.PID_P → 质子识别
• F.PID_K → K 介子识别
• F.PID_E → 电子识别
• F.PID_MU → μ 子识别
• F.PID_PI → π 介子识别

这些 PID变量的详细解释可参考此处。其基本原理是：三个提供PID 信息的子探测器（RICH、CALO、MUON系统）各自给出轨迹属于某类粒子的对数似然值。将这些值相乘得到一个总对数似然值。随后，与假设为π 介子时的对数似然求差，这些差值被称为 DLL（对数似然差值，ΔLog Likelihood）。因此，PID_K 就是把"假设为 K 介子"与"假设为 π 介子"的对数似然之差，以此类推。

还有一类更高级的 PID 函子称为 PROBNN。为了超越 PID变量这种较简单的对数似然方法，PROBNN 应运而生。它们还考虑了探测器子系统与追踪之间的相关性，通过训练包含这些特征的机器学习模型（MVA）来生成每种粒子的概率。其形式如下：

• F.PROBNN_K
• F.PROBNN_PI
• F.PROBNN_P
• F.PROBNN_E
• F.PROBNN_MU
• F.PROBNN_GHOST

这些是比 DLL 更强大的 PID 变量，在分析中广泛使用。

鬼径迹（GHOSTs）

LHCb中的一条真实径迹应对应一个真实粒子，其各层探测器中的命中点应良好对齐并指向物理顶点。然而，实际情况并非总是如此，即使经过触发和数据整理，我们仍会得到所谓"鬼径迹"。它们由子探测器中的随机命中组合产生，可能来自多个粒子的共同命中或噪声，并不对应真实粒子。因此，使用诸如 F.PROBNN_GHOST 的函子可以极大地帮助去除样本中的背景。

函子计算

我们不仅可以直接应用这些函子，还可以对它们进行计算。例如，我们可以使用以下函子：

F.MAX
F.MIN
F.SUM

在 DaVinci中对其他变量进行计算。这意味着最终在你的元组中会生成一个额外的分支来存储计算结果。例如，如果你想获取某个复合粒子的子粒子中的最大横向动量（PT），可以这样传递：

"MAX_PT": F.MAX(F.PT)

并将其应用于\(B_{s}\)或\(D_{s}\)粒子。

你还会看到一种重要的结构，例如：

"NHITS": F.VALUE_OR(-1) @ F.NHITS @ F.TRACK

这个结构应该从右到左解读。它的含义是，我们使用函子F.TRACK获取轨迹对象，然后对该轨迹对象应用F.NHITS以获取每条轨迹的击中数。最后，F.VALUE_OR(-1)表示如果对该轨迹应用F.NHITS，若没有返回值，则赋予其值为 -1。

函子运算符

你还可以使用诸如&（与）之类的运算符来要求多个条件同时成立。此外，你可以使用比较运算符如<（小于）和 >（大于）来指定数值阈值。下面是一个使用F.NINTREE来计算衰变树中满足以下条件的基本粒子子粒子数量的示例：

• 具有关联的轨迹（F.TRACK）。

• 横向动量（PT）大于 1500 MeV。

"N_HIGHPT_TRACKS": F.NINTREE(F.ISBASICPARTICLE & (F.HAS_VALUE @ F.TRACK) & (F.PT > 1500)),

DaVinci 脚本示例

现在我们将演示其中一些新函子的实际应用。我们可以使用之前的 yaml文件，下面是一个包含这些新函子的 Python 脚本示例，并在每个函子旁边做了简要定义。

关于单位的说明

在 LHCb软件中，默认使用的单位是MeV（兆电子伏特）、纳秒（nanoseconds）和毫米（millimetres）！

下面这个脚本的唯一其他变化是我们如何为衰变中的每个粒子分配不同的函子。之前我们只使用all_vars集合，它会将F.MASS、F.P和F.PT应用于衰变中的所有粒子。但我们可能有一些函子只想应用于复合粒子，例如\(B_{s}\)或\(D_{s}\)。你可以理解为什么要对这些粒子应用 DOCA变量，因为它返回的是其子粒子的最近接近距离。同样，我们可能有一些变量只想应用于轨迹，例如 PROBNN 变量，因为实际上是kaon、pion、质子等粒子在探测器中留下击中信息。因此，我们现在有 4 个函子集合：all_vars,base_composite_variables，Ds_extra_doca和track_variables。然后我们在变量集合中这样应用它们：

variables = {
    "ALL": all_vars,                                    # Variables for all particles
    "Bs": base_composite_variables,                     # Variables specific to Bs
    "Ds": base_composite_variables + Ds_extra_doca,     # Variables specific to Ds 
    "Bs_pi": track_variables,                           # Variables for pion
    "Ds_pi": track_variables,                           # Variables for pion 
    "Kp":    track_variables,                           # Variables for kaon 
    "Km":    track_variables,                           # Variables for kaon 
}

我们需要指定衰变中的粒子以及希望它们拥有哪些函子集合。记住，我们仍然有ALL类别，它会将该集合应用于候选粒子中的所有粒子！

轮到你了！

现在，轮到你用之前的 yaml 文件运行这个更新后的 Python脚本了。请确保你理解所使用的函子，并尝试从函子文档中添加更多函子。你可以使用与之前相同的yaml 文件，只需像下面这样更新 python 文件即可。

OWNPV 与 BPV 的区别

你可能会在其他示例中看到一些函子，例如F.BPVFD(pvs)与F.OWNPVFD。BPV代表"最佳初级顶点（Best Primary Vertex）"，使用该函子时需要传入初级顶点（pvs）的列表。OWNPV与之功能等效，但经证实可以节省CPU 时间，因为它无需传入初级顶点列表。

通常认为BPV已过时，在一般情况下应避免使用，而优先选择OWNPV。

'''
lb-run -c x86_64_v3-el9-gcc13+detdesc-opt+g DaVinci/v65r0 lbexec DV_intermediate:main tuple_options.yaml
'''
import Functors as F
import FunTuple.functorcollections as FC
from DaVinci import Options, make_config
from DaVinci.algorithms import create_lines_filter
from FunTuple import FunctorCollection
from FunTuple import FunTuple_Particles as Funtuple
from PyConf.reading import get_particles, get_pvs
from RecoConf.event_filters import require_pvs


def main(options: Options):
    line = "Hlt2B2OC_BdToDsmPi_DsmToKpKmPim"
    data = get_particles(f"/Event/HLT2/{line}/Particles")
    line_prefilter = create_lines_filter(name=f"PreFilter_{line}", lines=[line])
    pvs = get_pvs()

    fields = {
        "Bs": "[B0 -> (D_s- -> K- K+ pi-) pi+]CC",
        "Ds": "[B0 -> ^(D_s- -> K- K+ pi-) pi+]CC",
        "Bs_pi": "[B0 -> (D_s- -> K- K+ pi-) ^pi+]CC",
        "Ds_pi": "[B0 -> (D_s- -> K- K+ ^pi-) pi+]CC",
        "Km": "[B0 -> (D_s- -> ^K- K+ pi-) pi+]CC",
        "Kp": "[B0 -> (D_s- -> K- ^K+ pi-) pi+]CC",
    }

    all_vars = FunctorCollection({
        "M": F.MASS,
        "P": F.P,
        "PT": F.PT,
        "ENERGY": F.ENERGY,
        "PX": F.PX,
        "PY": F.PY,
        "PZ": F.PZ,
        "ID": F.PARTICLE_ID,            # PDG ID of the particle
        "Q": F.CHARGE,                  # Electric charge
        "ETA": F.ETA,                   # Pseudorapidity
        "PHI": F.PHI,                   # Azimuthal angle
        "CHI2": F.CHI2,                 # χ²
        "CHI2DOF": F.CHI2DOF,           # χ² of degrees of freedom
        "OWNPVIP": F.OWNPVIP,           # Impact parameter wrt own PV
        "OWNPVIPCHI2": F.OWNPVIPCHI2,   # Impact parameter χ² wrt own PV

    })

    base_composite_variables = FunctorCollection({
        "VTXCHI2NDOF": F.CHI2DOF,         # Vertex fit χ²/ndf
        "END_VX": F.END_VX,               # x-coordinate of decay vertex
        "END_VY": F.END_VY,               # y-coordinate of decay vertex
        "END_VZ": F.END_VZ,               # z-coordinate of decay vertex
        # OWNPV values
        "OWNPV_X": F.OWNPVX,              # x-coordinate of best PV
        "OWNPV_Y": F.OWNPVY,              # y-coordinate of best PV
        "OWNPV_Z": F.OWNPVZ,              # z-coordinate of best PV
        "OWNPV_DIRA": F.OWNPVDIRA,        # Direction angle cosine wrt own PV
        "OWNPV_FD": F.OWNPVFD,            # Flight distance wrt own PV
        "OWNPV_FDCHI2": F.OWNPVFDCHI2,    # Flight distance χ² wrt own PV
        "OWNPV_VDRHO": F.OWNPVVDRHO,      # Radial flight distance wrt own PV
        "OWNPV_VDZ": F.OWNPVVDZ,          # z-direction flight distance
        "OWNPV_LTIME": F.OWNPVLTIME,      # Proper lifetime
        "OWNPV_DLS": F.OWNPVDLS,          # Decay length significance
        # DOCA
        "DOCA12": F.DOCA(1, 2),           # DOCA between first and second daughter
        "DOCA12CHI2": F.DOCACHI2(1, 2),   # DOCA χ² between first and second daughter
        # Daughter Max, Min and Sums
        "MAX_PT": F.MAX(F.PT),            # Maximum PT of daughters
        "MIN_PT": F.MIN(F.PT),            # Minimum PT of daughters
        "SUM_PT": F.SUM(F.PT),            # Sum of daughters' PT
        "MAX_P": F.MAX(F.P),              # Maximum momentum of daughters
        "MIN_P": F.MIN(F.P),              # Minimum momentum of daughters
        "SUM_P": F.SUM(F.P),              # Sum of daughters' momentum
        "MAX_OWNPVIPCHI2": F.MAX(F.OWNPVIPCHI2),  # Max IP χ² of daughters
        "MIN_OWNPVIPCHI2": F.MIN(F.OWNPVIPCHI2),  # Min IP χ² of daughters
        "SUM_OWNPVIPCHI2": F.SUM(F.OWNPVIPCHI2),  # Sum of daughters' IP χ²
        "MAXDOCACHI2": F.MAXDOCACHI2,      # Maximum DOCA χ² between any daughters
        "MAXDOCA": F.MAXDOCA,              # Maximum DOCA between any daughters
        "MAXSDOCACHI2": F.MAXSDOCACHI2,    # Maximum signed DOCA χ²
        "MAXSDOCA": F.MAXSDOCA,            # Maximum signed DOCA

    })

    # Need extra DOCA combinations since Ds decays to 3 final state particles
    Ds_extra_doca = FunctorCollection({
        "DOCA13": F.DOCA(1, 3),
        "DOCA23": F.DOCA(2, 3),
        "DOCA13CHI2": F.DOCACHI2(1, 3),
        "DOCA23CHI2": F.DOCACHI2(2, 3),
    })

    track_variables = FunctorCollection({
        # Standard PID
        "PIDp": F.PID_P,              # Proton PID likelihood
        "PIDK": F.PID_K,              # Kaon PID likelihood
        "PIDPi": F.PID_PI,            # Pion PID likelihood
        "PIDe": F.PID_E,              # Electron PID likelihood
        "PIDmu": F.PID_MU,            # Muon PID likelihood
        # PROBNNs
        "PROBNN_pi": F.PROBNN_PI,        # Neural net probability of being a pion
        "PROBNN_p": F.PROBNN_P,          # Neural net probability of being a proton
        "PROBNN_K": F.PROBNN_K,          # Neural net probability of being a kaon
        "PROBNN_e": F.PROBNN_E,          # Neural net probability of being an electron
        "PROBNN_mu": F.PROBNN_MU,        # Neural net probability of being a muon
        "PROBNN_GHOST": F.PROBNN_GHOST,  # Neural net probability of being a ghost track
        "ISMUON": F.ISMUON,              # Boolean: is it identified as a muon 0 or 1?
        # Additional track related info
        # F.TRACK gets the track object
        # F.NDOF gets the NDOF for that track
        # F.VALUE_OR(-1) means if no value exists return -1 instead of failing
        "TRNDOF": F.VALUE_OR(-1) @ F.NDOF @ F.TRACK,                # NDOF in track fit, if this is higher then more hits used in fit
        "NHITS": F.VALUE_OR(-1) @ F.NHITS @ F.TRACK,                # Total number of hits in all detectors
        "NVPHITS": F.VALUE_OR(-1) @ F.NVPHITS @ F.TRACK,            # Total number of hits in VELO phi sensors
        "NUTHITS": F.VALUE_OR(-1) @ F.NUTHITS @ F.TRACK,            # Total number of hits in UT
        "NFTHITS": F.VALUE_OR(-1) @ F.NFTHITS @ F.TRACK,            # Total number of hits in Fibre Tracker (SciFi)
        "TRACKHISTORY": F.VALUE_OR(-1) @ F.TRACKHISTORY @ F.TRACK,  # Track reconstruction history
    })

    variables = {
        "ALL": all_vars,                                    # Variables for all particles
        "Bs": base_composite_variables,                     # Variables specific to Bs
        "Ds": base_composite_variables + Ds_extra_doca,     # Variables specific to Ds 
        "Bs_pi": track_variables,                           # Variables for pion
        "Ds_pi": track_variables,                           # Variables for pion 
        "Kp":    track_variables,                           # Variables for kaon 
        "Km":    track_variables,                           # Variables for kaon 
    }

    funtuple = Funtuple(
        name=line,
        tuple_name="DecayTree",
        fields=fields,
        variables=variables,
        inputs=data,
    )

    algs = {line: [line_prefilter, require_pvs(pvs), funtuple]}
    return make_config(options, algs)

这只是运行 DaVinci 的一个示例，完整的 DaVinci教程部分可参考此处。接下来我们可以

使用 root 命令打开 root 文件：

root -l tuple.root

进入我们的树所在目录：

Hlt2B2OC_BdToDsmPi_DsmToKpKmPim->cd()

列出树内容：

.ls

打印叶子的信息

DecayTree->GetListOfLeaves()->Print()

绘制一些变量

DecayTree->Draw("Bs_M")

VS Code 扩展

如果你使用 VS Code，还可以安装名为 "ROOT File Viewer"的扩展，它能很好地用于查看小型 root 文件！

运行上述命令后，应该会生成类似这样的图像（在本示例中，我们已在 optionsyaml 文件中增加了候选粒子的数量）

太棒了！如果能成功运行并得到这样的结果，你就可以继续学习下一部分关于函子集合的内容了！