ROOT Tutorial II：TTree、关联分析与事件选择
建议学完第一章后阅读。编程基础见 Tutorial I。
以 ΔE–E 望远镜为例，将同一粒子的探测器信号保存在 TTree 中，再通过逻辑 cut 和 graphical cut 选择粒子、构建能谱与 PID 投影。
按顺序运行 notebook。两种语言使用相同数据模型和筛选条件，选择一种即可。
0. 在 Python 中启动 ROOT¶
按顺序从头运行 notebook，后面的部分单元格会读取前面生成的文件。重启 kernel 会清除 Python 变量，但不会删除已写到磁盘的 ROOT 文件。
import numpy as np
import math
import ROOT
print("ROOT version:", ROOT.gROOT.GetVersion())
ROOT version: 6.40.02
ROOT 提供 PyROOT 接口，NumPy 数组用作 branch 的内存 buffer，math 用于数值计算。
%jsroot on 是 Jupyter 命令，用于开启交互画布，不是普通 Python 语句。在 Python 脚本中运行时去掉这一行。
%jsroot on
ROOT.gStyle.SetOptStat(0)
c1 = ROOT.TCanvas("c1", "Telescope analysis", 800, 520)
各图沿用这张画布。SetOptStat(0) 隐藏默认统计框，计数在需要时直接输出。
1. 探测器背景与事例记录
ΔE–E 望远镜测量什么
带电粒子进入硅后，通过电离和激发损失动能。硅探测器收集电离产生的电荷，信号经能量刻度后给出粒子在该层的能量沉积。将两层探测器前后排列，前层测量能量损失 ΔE，后层测量剩余能量 E，这就是 ΔE–E 望远镜。
本例中的粒子穿透前层、停止在后层：穿透表示离开该层时仍有动能，停止表示在该层耗尽剩余动能。忽略层间损失时，入射动能满足 \(E_0=\Delta E+E\)。如果粒子也穿出后层，两层信号之和就不再给出完整的入射能量。
为什么不同同位素形成条带
在非相对论能区，忽略缓慢变化的项时，阻止本领有 \(-dE/dx\propto Z^2/v^2\) 的变化趋势。这里 Z 和 A 表示入射粒子的电荷数和质量数。由 \(v^2\propto E_0/A\)，薄层中的能损近似随 \(Z^2A/E_0\) 变化。因此同一种粒子的 ΔE 和 E 彼此关联，不同种类形成不同的弯曲条带。相同 Z、相同速度的同位素阻止本领近似相同；在相同总动能下，它们的速度不同，能损也不同。
只有 ΔE 的一维分布会把不同能量、不同粒子的信号叠在一起。保留同一粒子的 E 后，才能在二维图中沿条带选择某种粒子。条带有宽度，来自能损涨落、测量分辨等因素，邻近条带也可能重叠。
event、entry 与 branch
一次入射粒子的完整记录称为事例（event）。TTree 通常用一个 entry 保存一个事例，用不同 branch 保存字段。两层信号放在同一个数组 branch 的两个元素中；读取某个 entry 时，两者仍属于同一粒子。直方图只累积 bin 计数，TTree 则保留逐事例数据，方便以后改变 cut 或重新计算观测量。
本例中的望远镜事例¶
设 p、d、t、³He、⁴He、⁶He 混合轻离子入射到薄硅 ΔE 探测器，后接足够厚、能够阻止粒子的 E 探测器。第一层测得 signal[0]，第二层测得剩余能量 signal[1]。两个值来自同一个粒子，因此保存在同一个 entry 中。
每条模拟记录包含：
eventID：按顺序生成的事例编号；A、Z：生成粒子的质量数和电荷数；E0：生成的入射动能，单位为 MeV；signal[0]：薄探测器测得的能量损失，单位为 MeV；signal[1]：厚探测器测得的能量，单位为 MeV。
生成的 (A, Z) = (1, 1)、(2, 1)、(3, 1)、(3, 2)、(4, 2)、(6, 2) 依次对应 p、d、t、³He、⁴He、⁶He。A、Z 和 E0 是模拟真值（simulation truth），可用于核对示例；普通实验数据不会自动给出这些真值。后面先根据测量信号筛选，再用真值标签检查。
生成示意性的望远镜数据
为练习 TTree，先生成一组简单的双层探测器数据。根据非相对论近似下 \(dE/dx\propto Z^2/v^2\) 和 \(v^2\propto E_0/A\) 的变化趋势，采用示意响应
$$\Delta E_{\rm true}=K\frac{Z^2 A}{E_0},\qquad E_{\rm true}=E_0-\Delta E_{\rm true},$$
其中 \(K=8\) MeV²，\(20\le E_0\le80\) MeV。假设粒子停止在第二层，再分别给两个信号加入 Gaussian 测量涨落。下面的参数仅供示例使用，不是实际探测器的刻度结果。这个解析模型只用于生成便于练习的条带，不代替第一章基于射程数据的能损计算；作业 4.1 使用作业 1.2 的三层探测器能量沉积。
2. 创建并填充 TTree
先打开输出文件，定义 branch，再逐事例填充。每调用一次 Fill()，就把当前全部字段保存为一条记录。
fout = ROOT.TFile("tree_telescope_python.root", "RECREATE")
tree_out = ROOT.TTree("telescope", "two-detector charged-particle telescope events")
eventID = np.zeros(1, dtype=np.int32)
A = np.zeros(1, dtype=np.int32)
Z = np.zeros(1, dtype=np.int32)
E0 = np.zeros(1, dtype=np.float32)
signal = np.zeros(2, dtype=np.float32)
TFile(filename, mode) 创建或打开 ROOT 文件。RECREATE 创建新文件，并覆盖已有同名文件，因此仅用于允许覆盖的输出文件。
TTree(name, title) 创建空 tree。telescope 是之后用 Get() 读取时使用的内部名称，title 用于说明内容。
基本类型的 branch 从固定内存位置读取待保存的值。np.zeros(1, dtype=np.int32) 创建长度为 1 的 32 位有符号整数数组，np.zeros(1, dtype=np.float32) 创建长度为 1 的 32 位浮点数组；第一个参数是长度，初值均为零。标量 branch 的值放在索引 0 处，signal 则用两个元素保存两层信号。调用 Branch() 之后保留这些数组对象，每次 Fill() 前只更新其元素；Fill() 将当前值复制到新的 entry 中。
tree_out.Branch("eventID", eventID, "eventID/I")
tree_out.Branch("A", A, "A/I")
tree_out.Branch("Z", Z, "Z/I")
tree_out.Branch("E0", E0, "E0/F")
tree_out.Branch("signal", signal, "signal[2]/F")
<cppyy.gbl.TBranch object at 0xc1f98ed00>
上面创建基本类型 branch 的形式为
tree.Branch(branch_name, memory_buffer, leaf_description)
第一个参数是文件中的 branch 名称，第二个参数是缓冲数组（buffer），Fill() 会复制其中的当前值。第三个参数是 leaf list：/ 前为 leaf 名称，后面的字母指定存储类型；/I 表示 32 位有符号整数，/F 表示 32 位浮点数。signal[2]/F 中的 [2] 表示每条记录含两个浮点元素。
buffer 的类型应与 leaf list 一致。例如将 np.float32 配成 /I，会使同一段内存被按错误类型解释。ROOT 不会自动推断物理单位，需在名称、标题或说明中注明。
rng = ROOT.TRandom3(2026)
n_events = 60000
K_loss = 8.0 # 示意模型中的常数，单位为 MeV^2
species_A = (1, 2, 3, 3, 4, 6)
species_Z = (1, 1, 1, 2, 2, 2)
对每个粒子，先更新粒子标签、入射能量和两个信号，再调用一次 Fill()，将这些值共同保存为一条记录。
for i in range(n_events):
eventID[0] = i
species_index = int(rng.Integer(6))
A[0] = species_A[species_index]
Z[0] = species_Z[species_index]
E0[0] = rng.Uniform(20.0, 80.0)
dE_true = K_loss * Z[0]**2 * A[0] / E0[0]
E_true = E0[0] - dE_true
sigma_dE = 0.03 + 0.04 * math.sqrt(dE_true)
sigma_E = 0.08 + 0.015 * math.sqrt(E_true)
signal[0] = max(0.0, rng.Gaus(dE_true, sigma_dE))
signal[1] = max(0.0, rng.Gaus(E_true, sigma_E))
tree_out.Fill()
print("Entries held in memory:", tree_out.GetEntries())
Entries held in memory: 60000
Integer(6) 等概率选取六种同位素之一；Gaus(mean, sigma) 生成测量信号，其中两种 sigma 表达式都是本示意模型的假设。Fill() 将当前 branch buffer 复制到下一条记录；不要重新创建 buffer，只更新其中的值。
fout.cd()
tree_out.Write()
fout.Close()
print("Created tree_telescope_python.root")
Created tree_telescope_python.root
fout.cd() 将输出文件设为当前 ROOT 目录，Write() 保存 TTree 的 branch 结构和已填入的记录，Close() 完成文件写入并关闭。仅调用 Fill() 不能保证 TTree 已完整保存到磁盘，还需要完成写入步骤。
3. 重新打开并查看 TTree
先检查文件对象、branch 名称与类型、数组长度和 entry 数，再进行分析。
fin = ROOT.TFile.Open("tree_telescope_python.root", "READ")
fin.ls()
tree_in = fin.Get("telescope")
TFile** tree_telescope_python.root TFile* tree_telescope_python.root KEY: TTree telescope;1 two-detector charged-particle telescope events
ls() 列出文件中的对象。按保存的名称 telescope 取出 TTree，并在分析期间保持输入文件打开。若文件或 tree 不存在，先核对文件名及对象列表。
tree_in.Print()
****************************************************************************** *Tree :telescope : two-detector charged-particle telescope events * *Entries : 60000 : Total = 1446637 bytes File Size = 806254 * * : : Tree compression factor = 1.79 * ****************************************************************************** *Br 0 :eventID : eventID/I * *Entries : 60000 : Total Size= 241125 bytes File Size = 84525 * *Baskets : 8 : Basket Size= 32000 bytes Compression= 2.85 * *............................................................................* *Br 1 :A : A/I * *Entries : 60000 : Total Size= 241053 bytes File Size = 39081 * *Baskets : 8 : Basket Size= 32000 bytes Compression= 6.16 * *............................................................................* *Br 2 :Z : Z/I * *Entries : 60000 : Total Size= 241053 bytes File Size = 25163 * *Baskets : 8 : Basket Size= 32000 bytes Compression= 9.56 * *............................................................................* *Br 3 :E0 : E0/F * *Entries : 60000 : Total Size= 241065 bytes File Size = 212170 * *Baskets : 8 : Basket Size= 32000 bytes Compression= 1.13 * *............................................................................* *Br 4 :signal : signal[2]/F * *Entries : 60000 : Total Size= 481883 bytes File Size = 444124 * *Baskets : 16 : Basket Size= 32000 bytes Compression= 1.08 * *............................................................................*
Print() 显示 entry 数、branch 名称、类型和维度。这里 signal[2]/F 表示每条记录包含两个 Float_t 值，GetEntries() 返回 entry 总数。
tree_in.Show(10)
======> EVENT:10
eventID = 10
A = 2
Z = 1
E0 = 71.276
signal = 0.266728,
71.1437
Show(entry_number) 输出某条记录的所有 branch。entry 编号从 0 开始，因此 Show(10) 显示的是第十一条记录。它用于检查数据，不是提取分析结果的方法。
tree_in.Scan(
"eventID:A:Z:E0:signal[0]:signal[1]",
"",
"",
10,
0,
)
************************************************************************************ * Row * eventID * A * Z * E0 * signal[0] * signal[1] * ************************************************************************************ * 0 * 0 * 2 * 1 * 65.996543 * 0.2959895 * 65.936271 * * 1 * 1 * 6 * 2 * 25.333940 * 7.5492363 * 17.606323 * * 2 * 2 * 6 * 2 * 57.311428 * 3.4103038 * 54.082107 * * 3 * 3 * 2 * 1 * 53.495006 * 0.3692665 * 53.215637 * * 4 * 4 * 4 * 2 * 38.148593 * 3.3014807 * 35.070316 * * 5 * 5 * 2 * 1 * 52.400577 * 0.3346374 * 52.181098 * * 6 * 6 * 1 * 1 * 69.937141 * 0.1737575 * 69.656311 * * 7 * 7 * 2 * 1 * 41.822879 * 0.4038591 * 41.643936 * * 8 * 8 * 1 * 1 * 21.835584 * 0.2760921 * 21.359891 * * 9 * 9 * 2 * 1 * 24.489458 * 0.6236841 * 23.949422 * ************************************************************************************
10
这里使用 Scan(expression, selection, option, number, first_entry)。用冒号分隔要显示的表达式；selection 和 option 留空表示不筛选、采用默认格式。最后两个参数指定从编号 0 开始显示十条记录。Scan() 适合抽查少量数据，不宜靠打印成千上万条记录完成分析或验证。
4. 二维关联、逻辑 cut 与 graphical cut
4.1 从二维关联图开始
Draw(expression, selection, option) 遍历 TTree。二维表达式写作 y:x，所以 signal[0]:signal[1] 表示横轴 E、纵轴 ΔE；TH2F 构造函数则先给 x 轴、再给 y 轴的 bin 设置。
hPID = ROOT.TH2F("hPID", "Telescope;E (MeV);#DeltaE (MeV)",
850, 0, 85, 1100, 0, 11)
tree_in.Draw("signal[0]:signal[1]>>hPID", "", "COLZ")
c1.Draw()
>>hPID 将计算值填入指定直方图，COLZ 用颜色表示计数。这里 bin 宽度为 0.1 MeV × 0.01 MeV。直接重复 Draw 会先清空同名直方图；>>+hPID 才是累加。
氢同位素的 ΔE 较小，在全范围图中挤在底部。下面只放大坐标显示，不重新生成或筛选数据。
hPID.GetXaxis().SetRangeUser(18, 38)
hPID.GetYaxis().SetRangeUser(0, 1.5)
hPID.Draw("COLZ")
labels = ROOT.TLatex()
labels.SetTextSize(0.035)
labels.DrawLatex(22, 0.22, "p")
labels.DrawLatex(22, 0.57, "d")
labels.DrawLatex(22, 0.92, "t")
c1.Draw()
局部图中由下到上是 p、d、t。增加 bin、放大坐标有助于看清原有结构，但不会改善探测器分辨；能量较高处的氢条带仍可能重叠。先用 SetRange(0, 0) 恢复全轴范围，后面继续分析全部事例。
hPID.GetXaxis().SetRange(0, 0)
hPID.GetYaxis().SetRange(0, 0)
tree_in.SetAlias("Etot", "signal[0]+signal[1]")
True
SetAlias 为表达式取别名。Etot 可用于后面的 Draw 和 cut，但不会增加一个 branch。
4.2 用 AND、OR、NOT 组合条件
cut 判断每条记录是否被接受。下面选择总能量在 20–40 MeV 或 60–80 MeV，并且两层信号均为正的事例。
| 写法 | 含义 | 例子 |
|---|---|---|
&&（AND） | 两项都成立 | Etot>=20 && Etot<40 |
||（OR） | 至少一项成立 | (Etot<40) || (Etot>=60) |
!（NOT） | 条件取反 | !(Etot>=40 && Etot<60) |
==、!= | 相等、不相等 | A==3 && Z==2（模拟中的 ³He 标签） |
无论外层代码是 PyROOT 还是 C++，传给 ROOT 的字符串表达式都采用上表写法。区间要写成两个比较，不写 20<Etot<40；混合 AND 和 OR 时，用括号明确组合。
valid = "signal[0]>0 && signal[1]>0"
low = "Etot>=20 && Etot<40"
high = "Etot>=60 && Etot<80"
windows = f"({valid}) && (({low}) || ({high}))"
print("Low window:", tree_in.GetEntries(f"({valid}) && ({low})"))
print("High window:", tree_in.GetEntries(f"({valid}) && ({high})"))
print("Either window:", tree_in.GetEntries(windows))
print("Outside both:", tree_in.GetEntries(f"({valid}) && !(({low}) || ({high}))"))
Low window: 20133 High window: 19885 Either window: 40018 Outside both: 19959
PyROOT 中可用字符串变量组合表达式；C++ 中的 TCut 是对条件字符串的封装，可用 &&、|| 组合。GetEntries(cut) 只统计满足条件的 entry，不画图。这两个能区不重叠，因此 OR 的计数等于两者相加；一般的 OR 不会重复计数同时满足两项的事例。
hWindows = ROOT.TH2F("hWindows", "Two total-energy windows;E (MeV);#DeltaE (MeV)",
850, 0, 85, 1100, 0, 11)
tree_in.Draw("signal[0]:signal[1]>>hWindows", windows, "COLZ")
c1.Draw()
图中保留两段总能量区间，边界由 ΔE + E 决定，因此不是竖直线。筛选字段可以不同于绘图字段：例如三层探测器中，可用第三层 e[2] 的信号筛选，再画前两层 e[0]:e[1] 的关联。
4.3 用 TCutG 圈选 ³He 条带
能量区间适合用逻辑 cut 表达；沿弯曲条带选择粒子时，封闭多边形更直接，称为 graphical cut。下面沿 ³He 条带画一个有限范围的区域，再将区域内的事例用于其他分析。
TCutG(name, n) 创建 n 个顶点的 cut，SetPoint(i, x, y) 按轮廓顺序设置顶点，最后一个点与第一个点相同以闭合。SetVarX、SetVarY 说明横纵坐标对应哪个 tree 表达式；这里分别为后层 E 和前层 ΔE，不要与 Draw 的 y:x 顺序混淆。
xgate = [22, 30, 40, 50, 60, 70, 70, 60, 50, 40, 30, 22, 22]
ygate = [4.03, 3.24, 2.56, 2.12, 1.80, 1.55, 1.12, 1.31, 1.57, 1.98, 2.56, 3.43, 4.03]
he3 = ROOT.TCutG("he3_cut", len(xgate))
he3.SetVarX("signal[1]")
he3.SetVarY("signal[0]")
for i in range(len(xgate)):
he3.SetPoint(i, xgate[i], ygate[i])
hPID.Draw("COLZ")
he3.SetLineColor(ROOT.kRed + 1)
he3.SetLineWidth(2)
he3.Draw("L SAME")
labels.DrawLatex(34, 3.7, "^{3}He gate")
c1.Draw()
顶点坐标来自这张图中的条带，只适用于本例。在 ROOT 桌面画布中，也可用图形编辑器的 CutG 工具逐点画出轮廓，再将默认名称 CUTG 改成有意义的名称。不同 notebook 前端的编辑功能不完全相同，直接设置顶点的写法便于复现和修改。
将 cut 名称 "he3_cut" 放入 selection，即可选择区域内的事例；也能继续写 "he3_cut && Etot<50" 或 "!he3_cut"。仅把轮廓画在图上不会筛选数据。
hHe3 = ROOT.TH2F("hHe3", "Selected by the graphical cut;E (MeV);#DeltaE (MeV)",
850, 0, 85, 1100, 0, 11)
tree_in.Draw("signal[0]:signal[1]>>hHe3", "he3_cut", "COLZ")
he3.Draw("L SAME")
c1.Draw()
n_gate = tree_in.GetEntries("he3_cut")
n_true = tree_in.GetEntries("he3_cut && A==3 && Z==2")
print("Selected entries:", n_gate)
print("True 3He among selected:", n_true)
print("Other isotopes among selected:", n_gate - n_true)
Selected entries: 7422 True 3He among selected: 7416 Other isotopes among selected: 6
这次选择只用了测量信号，随后才用模拟的 A、Z 核对成分。区域只覆盖一段条带，边缘也可能漏选或混入其他粒子，因此所选计数不等于全部 ³He 的产额。
4.4 从条带到 PID 投影
第一章用 PID 参数将弯曲条带变成近似水平带，再作一维投影。对本页的示意模型，由 \(\Delta E_{\rm true}=KZ^2A/E_0\) 和 \(E_0=\Delta E_{\rm true}+E_{\rm true}\)，可定义
\[P=\sqrt{\Delta E(\Delta E+E)}\approx\sqrt{KZ^2A}.\]
无测量涨落时，同一核素的 P 在这个模型中是常数；加入测量涨落后形成有宽度的峰。这个表达式与本例模型相配，并不是实际望远镜数据的通用 PID 公式。
tree_in.SetAlias("pid", "sqrt(signal[0]*Etot)")
hPidE = ROOT.TH2F("hPidE", "PID transformation;E (MeV);P (MeV)",
425, 0, 85, 320, 0, 16)
tree_in.Draw("pid:signal[1]>>hPidE", "", "goff")
hAllP = hPidE.ProjectionY("hAllP")
hGateP = ROOT.TH1D("hGateP", "", 320, 0, 16)
tree_in.Draw("pid>>hGateP", "he3_cut", "goff")
7422
goff 表示先填充、不立即画图。ProjectionY 沿 x 方向累加二维 bin，得到纵轴 P 的分布。cut 仍使用原来的 ΔE、E 坐标，但绘图量已经换成 P。
cPID = ROOT.TCanvas("cPID", "PID and projection", 1100, 460)
cPID.Divide(2, 1)
cPID.cd(1)
hPidE.Draw("COLZ")
cPID.cd(2)
hAllP.SetTitle("PID projection;P (MeV);Counts / 0.05 MeV")
hAllP.SetLineColor(ROOT.kBlack)
hAllP.SetMaximum(1.35 * hAllP.GetMaximum())
hAllP.Draw("HIST")
hGateP.SetLineColor(ROOT.kRed + 1)
hGateP.Draw("HIST SAME")
legendPID = ROOT.TLegend(0.55, 0.76, 0.88, 0.89)
legendPID.AddEntry(hAllP, "All events", "l")
legendPID.AddEntry(hGateP, "3He graphical cut", "l")
legendPID.Draw()
cPID.Draw()
左图显示变换后的条带，右图比较全部事例与 ³He cut 的投影。p、d、t、³He、⁴He、⁶He 的中心依次约为 2.83、4.00、4.90、9.80、11.31、13.86 MeV，来自 \(\sqrt{KZ^2A}\)。红色投影落在 ³He 峰附近；氢同位素投影有明显重叠，不能把每个局部起伏都解释成独立核素。
4.5 数组选择与事例选择
signal[0] 每个 entry 只取一个值；直接使用 signal 则展开数组元素，一个 entry 可能贡献多个值。若要判断整个事例，可用 Sum$、Max$ 等在同一 entry 内汇总数组。
print("Both positive (explicit):", tree_in.GetEntries("signal[0]>0 && signal[1]>0"))
print("Both positive (array):", tree_in.GetEntries("Sum$(signal>0)==2"))
print("At least one above 40 MeV:", tree_in.GetEntries("Max$(signal)>40"))
Both positive (explicit): 59977 Both positive (array): 59977 At least one above 40 MeV: 37990
Sum$(signal>0) 统计本事例中正信号的个数，Max$(signal) 取本事例的最大信号，不是在所有 entry 之间求和或取最大值。三层探测器和波形数组采用同样的规则。
Draw 的 selection 实际是填充权重：逻辑条件得到 0 或 1。若直接写 "signal[0]"，含义是用能量作为权重，而不是“选择这个 branch”。本例均用逻辑条件，直方图内容就是所选事例的计数。
5. 在逐事例循环中复用 cut
复杂的逐事例修正、波形积分等计算更适合显式循环。这里保留同一个 ³He 区域，用 IsInside(x, y) 判断当前事例是否在多边形内，并计算所选粒子的总能量谱。它的参数顺序是 x、y，即 E、ΔE。
for i, event in enumerate(tree_in):
print(
int(event.eventID),
int(event.A),
int(event.Z),
float(event.signal[0]),
float(event.signal[1]),
)
if i == 4:
break
0 2 1 0.2959895730018616 65.93627166748047 1 6 2 7.549236297607422 17.6063232421875 2 6 2 3.410303831100464 54.08210754394531 3 2 1 0.3692665994167328 53.21563720703125 4 4 2 3.301480770111084 35.070316314697266
这里只读五条记录，用来查看访问方式。int()、float() 将 ROOT 标量代理转成普通 Python 数值，便于输出；若直接将数值传给 ROOT 直方图方法，通常不需要这些转换。
先用 Draw 得到对照能谱，再建立相同 bin 的空直方图用于循环。这里两幅能谱都只填入 cut 内的事例，不混用不同粒子的信号。
hEnergyDraw = ROOT.TH1F("hEnergyDraw", "3He gate;#DeltaE + E (MeV);Counts / MeV", 85, 0, 85)
hEnergyLoop = ROOT.TH1F("hEnergyLoop", "", 85, 0, 85)
tree_in.Draw("Etot>>hEnergyDraw", "he3_cut", "goff")
7422
selected_entries = 0
for event in tree_in:
dE = float(event.signal[0])
E = float(event.signal[1])
if he3.IsInside(E, dE):
hEnergyLoop.Fill(dE + E)
selected_entries += 1
print("Selected by IsInside:", selected_entries)
Selected by IsInside: 7422
这个循环与 Draw(..., "he3_cut") 应选择完全相同的事例。循环中的条件是普通程序语句：C++ 使用 &&、||、!，Python 使用 and、or、not，区别于前面的 ROOT 字符串表达式。
c1.cd()
hEnergyDraw.SetLineColor(ROOT.kBlack)
hEnergyDraw.SetMaximum(1.35 * hEnergyDraw.GetMaximum())
hEnergyDraw.Draw("HIST")
hEnergyLoop.SetMarkerStyle(20)
hEnergyLoop.SetMarkerSize(0.6)
hEnergyLoop.SetMarkerColor(ROOT.kRed + 1)
hEnergyLoop.Draw("P SAME")
legendE = ROOT.TLegend(0.60, 0.73, 0.88, 0.88)
legendE.AddEntry(hEnergyDraw, "TTree.Draw", "l")
legendE.AddEntry(hEnergyLoop, "Loop + IsInside", "p")
legendE.Draw()
c1.Draw()
max_difference = 0.0
for b in range(hEnergyDraw.GetNbinsX() + 2):
difference = abs(hEnergyDraw.GetBinContent(b) - hEnergyLoop.GetBinContent(b))
max_difference = max(max_difference, difference)
print("Largest bin-content difference:", max_difference)
Largest bin-content difference: 0.0
两种实现应逐 bin 一致，比较也包含 underflow 和 overflow。能谱的两端受 graphical cut 覆盖范围限制；这里展示的是所选区域内的总能量分布，不是完整的 ³He 入射能谱。
用 GetEntry() 按索引读取¶
直接遍历通常是最清楚的 PyROOT 写法，也可以按索引读取。GetEntry(i) 将编号为 i 的记录加载到 tree 对象中，之后用 tree.branch 读取各 branch。需要指定 entry 编号，或改写已有 C++ 循环时，可采用此方式。
for i in range(3):
bytes_read = tree_in.GetEntry(i)
print(
"entry", i,
"bytes", bytes_read,
"eventID", int(tree_in.eventID),
"dE", float(tree_in.signal[0]),
"E", float(tree_in.signal[1]),
)
entry 0 bytes 24 eventID 0 dE 0.2959895730018616 E 65.93627166748047 entry 1 bytes 24 eventID 1 dE 7.549236297607422 E 17.6063232421875 entry 2 bytes 24 eventID 2 dE 3.410303831100464 E 54.08210754394531
GetEntry() 返回读取的字节数，不返回事例对象；返回值非正表示没有读到 entry 数据。普通 PyROOT 标量分析通常不必显式调用 SetBranchAddress()，因为 branch 已可作为属性访问。C++ 版本将演示编译代码中常用的 buffer 与 SetBranchAddress() 写法。
6. 保存 cut、直方图与所选事例
将所选直方图和 cut 一起保存，下次可以复查选择区域。若后续只分析这些事例，CopyTree(cut) 可复制满足条件的完整 entry，包括所有原有 branch。
fanalysis = ROOT.TFile("telescope_analysis_python.root", "RECREATE")
fanalysis.cd()
he3.Write()
hPID.Write()
hAllP.Write()
hGateP.Write()
hEnergyDraw.Write()
selected_tree = tree_in.CopyTree("he3_cut")
selected_tree.Write("he3")
print("Saved entries:", selected_tree.GetEntries())
fanalysis.Close()
Saved entries: 7422
fanalysis.cd() 将新 tree 建立在输出文件中；写出的内部名称为 he3。原始输入 TTree 没有改变。he3_cut 则保存了多边形顶点和坐标表达式，重新读取时仍应核对它与待分析数据的 branch、单位和刻度是否一致。
7. 作业中的定长数组
同样的 branch 写法也可保存探测器能量或波形采样点：
| 数据来源 | leaf 描述 | 一条 entry 的内容 |
|---|---|---|
| 本例 | signal[2]/F | 两层探测器的信号 |
| 作业 4.1 | e[3]/F, A/I, Z/I | D1、D2、D3 的能量沉积及粒子标签 |
| 作业 5.1 | tree wave 中的 adc[250]/D | 一个脉冲的 250 个采样点 |
/F 对应 Float_t，/D 对应 Double_t，/I 对应 Int_t。内存中数组的类型和长度应与文件中的定义一致。
对于波形，GetEntry(i) 一次加载整个脉冲，再通过内层循环处理采样数组，计算基线或电荷积分。下一次 GetEntry 会用下一个脉冲替换 buffer 中的内容。
PyROOT 也支持 C++ 式的 buffer 读取。以上面的 signal branch 为例：
signal_read = np.zeros(2, dtype=np.float32)
tree_in.SetBranchAddress("signal", signal_read)
tree_in.GetEntry(0)
print("First event signals:", signal_read[0], signal_read[1])
tree_in.ResetBranchAddresses()
First event signals: 0.29598957 65.93627
SetBranchAddress 将数组与 branch 连接，GetEntry 将数据读入数组，结束后用 ResetBranchAddresses 解除连接。读取 adc[250]/D 时，则准备含 250 个 np.float64 元素的数组。
fin.Close()