先按需阅读 Python 基础。本页示例介绍作业 1–3 中用到的 ROOT 操作。
在 PyROOT notebook 中按顺序运行单元格。选择一种代码语言即可，两种版本采用相同的示例和文件格式。
0. 交互式使用 ROOT¶
用 import ROOT 加载 PyROOT，随后通过 ROOT 模块调用 ROOT 的类和函数。下面先输出版本号，便于排查不同系统上的运行差异。ROOT.gROOT 是 ROOT 的全局对象，这里只使用它的 GetVersion() 方法。
import ROOT
print(f"ROOT version: {ROOT.gROOT.GetVersion()}")
ROOT version: 6.40.02
JSROOT 可在 Jupyter 中交互式显示 ROOT 图。%jsroot on 是 Jupyter 命令，不是普通 Python 语句，通常只需运行一次。若当前安装不提供该命令，可跳过此单元格，使用该环境已配置的显示方式。
%jsroot on
ROOT 将函数、图和直方图绘制在画布 TCanvas 上。下面使用的构造函数为
ROOT.TCanvas(name, title, width, height)
name 是画布在 ROOT 中的内部名称，同一会话中应避免重复；title 是画布标题；width 和 height 是初始宽、高，单位为像素。本教程重复使用同一画布：Clear() 清除先前的绘图而不删除这里绘制的对象，Draw() 则在 notebook 中显示当前画布。
SetOptStat(0) 隐藏自动统计框；需要的数值在相应代码中输出。
ROOT.gStyle.SetOptStat(0)
c1 = ROOT.TCanvas("c1", "ROOT Tutorial I", 800, 520)
ROOT 的数学函数可通过 ROOT.TMath 调用。先用下面的单元格熟悉调用方式，再创建 ROOT 对象。
x = 2.0
print("sqrt(2) =", ROOT.TMath.Sqrt(x))
print("pi =", ROOT.TMath.Pi())
sqrt(2) = 1.4142135623730951 pi = 3.141592653589793
类（class）定义对象的类型，例如 ROOT.TCanvas 和 ROOT.TF1。调用构造函数可创建一个对象（object），如 c1 是一个具体画布。对象提供的操作称为方法（method），如 c1.Clear()。
PyROOT 用点号调用方法：
object.Method(arguments)
ROOT 文档常给出 C++ 写法 object->Method(arguments)。PyROOT 用 . 代替 ->，也不需要写 C++ 中的指针声明 Type* 和 new；调用的 ROOT 类和方法本身相同。
1. TF1：一维函数¶
TF1 表示单个自变量的数学函数，常用于能量刻度曲线、峰形和本底模型。它保存公式和参数值，不保存逐个测量事例。
本节使用的构造函数为
ROOT.TF1(name, formula, xmin, xmax)
name 是 ROOT 内部名称，formula 是包含数学表达式的字符串；xmin、xmax 给出函数计算和绘制的范围。自变量写作 x，可调参数依次写作 [0]、[1] 等。
1.1 线性刻度函数¶
探测器及其电子学通常给出 ADC channel，而不是以 keV 为单位的能量。在适当范围内，可以先用线性关系描述能量刻度：
$$ E = a + b \times \text{channel}. $$
ROOT 公式中的 [0] 表示零点偏移 a，[1] 表示增益 b。下面将函数定义在 channel 0–4000 范围内。
fcal = ROOT.TF1("fcal", "[0] + [1]*x", 0.0, 4000.0)
fcal.SetParNames("offset", "gain")
fcal.SetParameters(0.2, 0.5)
fcal.SetTitle("Linear calibration;Channel;Energy (keV)")
c1.Clear()
fcal.Draw()
c1.Draw()
SetParNames 为 [0]、[1] 设置便于辨认的参数名称，不设置数值。SetParameters(0.2, 0.5) 再按顺序设置零点偏移和增益，得到
$$ E = 0.2 + 0.5 \times \text{channel}. $$
ROOT 标题字符串用分号分隔图标题、x 轴标题和 y 轴标题：
main title;x-axis title;y-axis title
fcal.Draw() 将函数画在当前画布上，c1.Draw() 将画布显示在 Jupyter 中。这里的参数已经直接给定，没有进行拟合。
Eval(value) 计算指定自变量处的函数值，可用于按已知刻度关系将 channel 换算为能量。
channel = 1500.0
energy = fcal.Eval(channel)
print(f"Channel {channel:.0f} -> Energy = {energy:.1f} keV")
Channel 1500 -> Energy = 750.2 keV
1.2 Gaussian 峰形¶
辐射探测器的全能峰（full-energy peak）常先用 Gaussian 函数近似：
$$ f(E)=A\exp\left[-\frac{1}{2}\left(\frac{E-\mu}{\sigma}\right)^2\right]. $$
参数 A、mu、sigma 分别表示峰高、峰位和标准差，对应 ROOT 表达式中的 [0]、[1]、[2]。下面在 400–700 keV 范围内绘制函数。
fpeak = ROOT.TF1(
"fpeak",
"[0]*exp(-0.5*((x-[1])/[2])^2)",
400.0,
700.0,
)
fpeak.SetParNames("height", "centroid", "sigma")
fpeak.SetParameters(1000.0, 550.0, 20.0)
fpeak.SetTitle("Gaussian peak model;Energy (keV);Counts")
c1.Clear()
fpeak.Draw()
c1.Draw()
曲线在 550 keV 处的峰高为 1000，宽度参数为 \(\sigma=20\) keV。Gaussian 的半高全宽（FWHM）为
$$\mathrm{FWHM}=2\sqrt{2\ln 2}\,\sigma\approx2.355\sigma.$$
TF1 在这里仅定义曲线；第 3、4 节将生成逐个事例的能量，并填入直方图。
2. TGraph：数据点¶
TGraph 保存成对坐标 (x_i, y_i)，适合表示每个点都有明确坐标的数据，例如刻度峰的 channel 与已知能量。它保存传入的坐标，不划分 bin，也不自动统计重复测量的次数。
2.1 刻度数据点与误差
下面给出五组示意性的 channel–能量数据。实验中，每组数据来自一个已识别的 γ 峰：峰位拟合给出 channel 及其误差，核数据给出参考能量。这里假设参考能量误差可以忽略，并为 channel 给定一组示意误差，单位也是 channel。
数组中的同一索引对应同一个峰。误差应放在它所属的坐标轴上，不能将 channel 误差直接当作能量误差。
import numpy as np
ch = np.array([512.0, 1024.0, 1538.0, 2049.0, 2561.0], dtype=np.float64)
E = np.array([122.1, 245.4, 367.8, 489.0, 612.3], dtype=np.float64)
dch = np.array([1.0, 1.0, 2.0, 2.0, 3.0], dtype=np.float64)
dEref = np.zeros(len(ch), dtype=np.float64)
TGraph(n, x, y) 保存 n 对坐标；TGraphErrors(n, x, y, ex, ey) 进一步保存两个坐标的误差。下面横轴为 channel，纵轴为参考能量，所以第四、第五个数组分别为 channel 误差和能量误差。
SetName 设置保存到 ROOT 文件时使用的名称，SetMarkerStyle(20) 选择实心圆点。Draw("AP") 中 A 表示坐标轴，P 表示数据点；TGraphErrors 同时显示误差棒。
gcal = ROOT.TGraphErrors(len(ch), ch, E, dch, dEref)
gcal.SetName("gcal")
gcal.SetTitle("Energy calibration;Channel;Energy (keV)")
gcal.SetMarkerStyle(20)
c1.Clear()
gcal.Draw("AP")
c1.Draw()
2.2 拟合刻度数据
ROOT 内置公式 pol1 表示 \(f(x)=[0]+[1]x\)。从数据大致估计斜率约为 0.24 keV/channel，用它设置拟合初值。
Fit(flin, "RSF") 中：R 使用函数定义的范围，S 返回拟合结果，F 让多项式也使用一般数值 minimizer。这里 channel 误差在 x 轴上，F 很重要：ROOT 的专用 linear fitter 会忽略 x 误差。一般拟合按函数斜率将 x 误差折算到能量方向。
flin = ROOT.TF1("flin", "pol1", 0.0, 3000.0)
flin.SetParNames("offset", "slope")
flin.SetParameters(0.0, 0.24)
calibration_fit = gcal.Fit(flin, "RSF")
c1.Draw()
**************************************** Minimizer is Minuit2 / Migrad Chi2 = 10.8919 NDf = 3 Edm = 5.21828e-07 NCalls = 39 offset = 0.0914136 +/- 0.307312 slope = 0.239022 +/- 0.000258847
GetParameter(i) 读取拟合参数 i。参数 0 是零点偏移，单位为 keV；参数 1 是斜率，单位为 keV/channel。单位由坐标轴对应的物理量确定，ROOT 不会从数值中自动推断单位。
offset = flin.GetParameter(0)
slope = flin.GetParameter(1)
print("Fit status:", int(calibration_fit))
print("Offset:", offset, "+/-", flin.GetParError(0), "keV")
print("Slope:", slope, "+/-", flin.GetParError(1), "keV/channel")
Fit status: 0 Offset: 0.09141363127154742 +/- 0.3073116485098406 keV Slope: 0.23902166824139454 +/- 0.0002588474870812204 keV/channel
状态码 0 表示数值拟合正常完成；模型是否合适还要结合残差判断。GetParError(i) 返回第 i 个参数的拟合误差。由多个参数共同计算一个量时，还需要考虑参数之间的协方差，参见加权拟合与误差传播。
2.3 逐点构造残差图
在完整刻度范围内，小的偏离容易被坐标尺度掩盖。逐点计算 \(E_{\rm ref}-f(\mathrm{channel})\)，用残差图查看偏离的位置和大小。
SetPoint(i, x, y) 设置第 i 个点，SetPointError(i, ex, ey) 设置其误差，索引从 0 开始。这里的误差棒显示原数据的测量误差尺度，折算为 \(|b|\sigma_{\rm channel}\)。
cal_residual = ROOT.TGraphErrors()
cal_residual.SetName("cal_residual")
for i in range(len(ch)):
cal_residual.SetPoint(i, ch[i], E[i] - flin.Eval(ch[i]))
cal_residual.SetPointError(i, 0, abs(slope) * dch[i])
cal_residual.SetTitle("Calibration residual;Channel;Reference - fitted energy (keV)")
cal_residual.SetMarkerStyle(20)
c1.Clear()
cal_residual.Draw("AP")
zero_cal = ROOT.TLine(ch[0], 0, ch[-1], 0)
zero_cal.SetLineStyle(2)
zero_cal.Draw()
c1.Draw()
虚线对应零残差。本例有些点的偏离超过了给定的测量误差尺度，说明拟合收敛不等于数据与模型符合良好。实际刻度时，应检查峰的对应关系、峰位误差和线性刻度的适用范围；若残差随 channel 连续弯曲，尤其需要检查最后一项。
2.4 读入文本表格并插值
将calibration_points.txt 放在 notebook 所在目录。文件中第一列为 channel，第二列为能量（keV），与前面的五组坐标相同；该文件不含误差列。
TGraph(filename, "%lg %lg") 从每行读取两个浮点数，依次作为 x、y。先用 GetN() 检查读入的点数，再计算插值。
table = ROOT.TGraph("calibration_points.txt", "%lg %lg")
print("Read points:", table.GetN())
print("Interpolated energy:", table.Eval(1000.0), "keV")
print("Fitted energy:", flin.Eval(1000.0), "keV")
Read points: 5 Interpolated energy: 239.6203125 keV Fitted energy: 239.11308187266607 keV
TGraph::Eval 默认在相邻数据点之间做线性插值；TF1::Eval 则计算拟合函数值，所以两个结果不必相同。用于射程、阻止本领等表格时，同样先确认两列的物理量、单位和覆盖范围，再选择合理的插值。
3. TRandom3：随机采样
Gaus 生成 Gaussian 测量涨落，Uniform 生成均匀分布数值，Poisson 生成固定时长内的计数。这三种方法都可由同一个随机数发生器调用。
构造函数为
ROOT.TRandom3(seed)
随机种子（seed）确定发生器的初始状态。每次用同一个正整数种子重新创建发生器，会得到相同序列，便于复现教程结果；换一个种子，则得到来自同一概率分布的另一组随机样本。
先看 Rndm()：它不接收参数，每次返回一个 0–1 之间均匀分布的随机数，并将发生器的内部状态推进到下一步。
rng = ROOT.TRandom3(12345)
for i in range(5):
print(f"sample {i}: {rng.Rndm():.6f}")
sample 0: 0.929616 sample 1: 0.890155 sample 2: 0.316376 sample 3: 0.130707 sample 4: 0.183919
用相同的正整数种子重新创建发生器，可以复现该序列。
3.1 Gaussian 展宽：模拟探测器分辨率¶
假设单能辐射沉积的真实能量为 548 keV，有限的探测器分辨率使重复测量值在附近涨落。方法
rng.Gaus(mean, sigma)
返回一个 Gaussian 随机数。mean 是中心值，sigma 是标准差，单位与结果相同。因此 Gaus(548.0, 8.0) 可表示一次可能的能量测量，单位为 keV。
E_true = 548.0 # keV
sigma_E = 8.0 # keV
for i in range(10):
E_measured = rng.Gaus(E_true, sigma_E)
print(f"measurement {i}: {E_measured:.3f} keV")
measurement 0: 537.891 keV measurement 1: 553.223 keV measurement 2: 541.708 keV measurement 3: 540.530 keV measurement 4: 542.628 keV measurement 5: 554.556 keV measurement 6: 554.827 keV measurement 7: 561.549 keV measurement 8: 553.479 keV measurement 9: 547.698 keV
每次调用代表对同一真实能量的一次可能测量，十个结果并不是十种不同的 γ 射线能量。这里用 Gaussian 近似只考虑分辨率展宽，完整的探测器响应还可能包含其他成分。
3.2 均匀本底¶
方法
rng.Uniform(xmin, xmax)
返回 xmin 与 xmax 之间的实数。本例令 0–2000 keV 内的概率密度相同，用平坦本底演示不同分布的组合，不将其视为一般 γ 连续谱的物理模型。
for i in range(10):
E_background = rng.Uniform(0.0, 2000.0)
print(f"background sample {i}: {E_background:.3f} keV")
background sample 0: 747.842 keV background sample 1: 1307.140 keV background sample 2: 309.947 keV background sample 3: 1495.430 keV background sample 4: 1784.687 keV background sample 5: 1922.613 keV background sample 6: 53.579 keV background sample 7: 16.777 keV background sample 8: 583.005 keV background sample 9: 212.889 keV
3.3 Poisson 计数涨落¶
事例相互独立、平均发生率恒定时，固定时长内的计数通常可用 Poisson 分布描述。方法
rng.Poisson(mean_count)
接收该时长内的期望计数，返回一个整数计数；参数不是计数上限。Poisson 分布的方差等于均值，因此标准差为 sqrt(mean_count)。
mean_count = 10.0 # 平均计数率为 10 Hz，测量 1 s
for i in range(10):
print(rng.Poisson(mean_count))
18 13 10 12 11 9 6 9 15 11
Gaus 和 Uniform 返回能量等连续量，Poisson 返回离散的事例数。选哪种分布取决于所描述的实验量，而不是它们是否属于同一个随机数发生器。
4. TH1：一维直方图¶
直方图表示许多测量值或模拟值的分布。横轴划分为若干小区间，称为 bin。Fill(value) 找到 value 所在的 bin，并将其计数加一；直方图保存各 bin 的计数，不保存完整的原始数值列表。
TH1 是 ROOT 一维直方图的基类。本例使用 TH1F，其中 F 表示 bin 内容用单精度浮点数保存；TH1D 则使用双精度。对于本例的计数规模，两者都够用。
4.1 构造简单能谱¶
构造函数为
ROOT.TH1F(name, title, number_of_bins, xmin, xmax)
下面的直方图在 0–2000 keV 范围内有 400 个 bin，每个宽 5 keV。低于或超出范围的值进入专门的 underflow 或 overflow bin。
模型包含 30,000 个 548 keV 峰事例、18,000 个 1250 keV 峰事例，以及 30,000 个平坦本底事例。两个峰采用不同的 sigma，表示不同的测量宽度。这是一个简化能谱，用来练习组合随机样本并观察直方图。
hspec = ROOT.TH1F(
"hspec",
"Simulated energy spectrum;Energy (keV);Counts per 5 keV",
400,
0.0,
2000.0,
)
r1 = ROOT.TRandom3(24680)
for i in range(30000):
hspec.Fill(r1.Gaus(548.0, 10.0))
for i in range(18000):
hspec.Fill(r1.Gaus(1250.0, 18.0))
for i in range(30000):
hspec.Fill(r1.Uniform(0.0, 2000.0))
c1.Clear()
hspec.Draw()
c1.Draw()
每个循环加入一种成分，每次生成一个能量并传给 Fill。循环次数控制该成分的事例总数，不是峰的最大高度，因为峰事例会分布在多个 bin 中。不带选项的 Draw() 使用 ROOT 默认的一维直方图样式。
4.2 查看直方图信息¶
GetEntries() 返回 Fill 的调用次数，包括进入 underflow 和 overflow 的事例。GetMean()、GetStdDev() 描述直方图整体分布；GetNbinsX() 返回 x 轴普通 bin 的个数，GetBinWidth(1) 返回第一个 bin 的宽度。本例所有普通 bin 等宽。
print(f"Entries = {hspec.GetEntries():.0f}")
print(f"Number bins = {hspec.GetNbinsX()}")
print(f"Bin width = {hspec.GetXaxis().GetBinWidth(1):.1f} keV")
print(f"Mean = {hspec.GetMean():.3f} keV")
print(f"Std dev = {hspec.GetStdDev():.3f} keV")
Entries = 78000 Number bins = 400 Bin width = 5.0 keV Mean = 882.077 keV Std dev = 456.857 keV
总事例数应为 30000 + 18000 + 30000 = 78000。输出中的 Mean 和 Std dev 对应两个峰加本底的整体分布，不是某个峰的峰位和分辨率。单个峰的性质应通过选定峰区或适当拟合提取。
4.2.1 用 sidebands 拟合局部本底
先用 480–510 keV 和 590–620 keV 两段 sidebands 描述峰附近的局部本底。每个 bin 的中心和计数构成一个 TGraphErrors 点，计数误差取 \(\sqrt{n}\)。只用这些点拟合 pol1，避免 Gaussian 峰影响本底初值。
红点是参与本底拟合的 bin，绿色直线是拟合结果，并向峰下延伸。sidebands 应避开峰尾和邻峰；真实能谱还需要检查线性本底是否足够。
sideband = ROOT.TGraphErrors()
for b in range(hspec.FindBin(480), hspec.FindBin(620)):
x = hspec.GetBinCenter(b)
if x < 510 or x >= 590:
y = hspec.GetBinContent(b)
i = sideband.GetN()
sideband.SetPoint(i, x, y)
sideband.SetPointError(i, 0, ROOT.TMath.Sqrt(max(y, 1.0)))
background_line = ROOT.TF1("background_line", "pol1", 480, 620)
background_line.SetParameters(75, 0)
background_result = sideband.Fit(background_line, "RSQ0")
print("Background fit status:", int(background_result))
print("p0:", background_line.GetParameter(0))
print("p1:", background_line.GetParameter(1), "counts/bin/keV")
hspec.GetXaxis().SetRangeUser(470, 630)
c1.Clear()
hspec.Draw("E")
sideband.SetMarkerStyle(20)
sideband.SetMarkerColor(ROOT.kRed + 1)
sideband.Draw("P SAME")
background_line.SetLineColor(ROOT.kGreen + 2)
background_line.Draw("SAME")
c1.Draw()
Background fit status: 0 p0: 92.75668910048645 p1: -0.032710868543318346 counts/bin/keV
RSQ0 中 R 使用函数范围，S 返回拟合结果，Q 减少屏幕输出，0 表示暂不把拟合函数画在 graph 上。这里随后把直线叠加到完整峰区。本底拟合的两个参数将在下一步作为联合拟合的初值。
4.3 峰与本底的联合拟合
在 480–620 keV 内使用 gaus(0)+pol1(3)：参数 0、1、2 是 Gaussian 的 height、mean 和 sigma，参数 3、4 是线性本底的截距和斜率。
Gaussian 的 mean 取峰区最大 bin 的中心，height 取该 bin 计数减去本底估计，sigma 从图上估计为 10 keV；线性本底参数沿用 sideband 拟合结果。合理的初值可以避免 minimizer 从与数据相差很远的位置开始搜索。
peak_bin = hspec.FindBin(510)
for b in range(hspec.FindBin(510), hspec.FindBin(590)):
if hspec.GetBinContent(b) > hspec.GetBinContent(peak_bin):
peak_bin = b
mean0 = hspec.GetBinCenter(peak_bin)
height0 = hspec.GetBinContent(peak_bin) - background_line.Eval(mean0)
sigma0 = 10.0
fit_peak_548 = ROOT.TF1("fit_peak_548", "gaus(0)+pol1(3)", 480, 620)
fit_peak_548.SetParNames("height", "mean", "sigma", "background", "slope")
fit_peak_548.SetParameters(
height0, mean0, sigma0,
background_line.GetParameter(0), background_line.GetParameter(1)
)
fit_peak_548.SetParLimits(0, 0, 100000)
fit_peak_548.SetParLimits(1, 510, 590)
fit_peak_548.SetParLimits(2, 0.1, 40)
print("Initial height, mean, sigma:", height0, mean0, sigma0)
print("Initial background p0, p1:", background_line.GetParameter(0), background_line.GetParameter(1))
Initial height, mean, sigma: 5897.152511426981 547.5 10.0 Initial background p0, p1: 92.75668910048645 -0.032710868543318346
Fit(..., "LIRS") 中：L 使用 binned Poisson likelihood；I 在每个 bin 内对模型积分取平均；R 使用函数定义的拟合范围；S 返回拟合结果及 covariance matrix。联合拟合会同时调整峰和本底参数，sideband 结果只是初值。
下图红线是完整的 Gaussian + linear background 模型，绿色虚线是联合拟合得到的本底分量。两条曲线都覆盖 sidebands 和 peak 区域。
peak_fit_result = hspec.Fit(fit_peak_548, "LIRS")
fitted_background = ROOT.TF1("fitted_background", "pol1", 480, 620)
fitted_background.SetParameters(
fit_peak_548.GetParameter(3), fit_peak_548.GetParameter(4)
)
hspec.GetXaxis().SetRangeUser(470, 630)
c1.Clear()
hspec.Draw("E")
fit_peak_548.SetLineColor(ROOT.kRed + 1)
fit_peak_548.Draw("SAME")
fitted_background.SetLineColor(ROOT.kGreen + 2)
fitted_background.SetLineStyle(2)
fitted_background.Draw("SAME")
c1.Draw()
**************************************** Minimizer is Minuit2 / Migrad MinFCN = 14.8711 Chi2 = 29.7421 NDf = 23 Edm = 4.78666e-09 NCalls = 144 height = 5962.56 +/- 43.8631 (limited) mean = 547.913 +/- 0.0617975 (limited) sigma = 10.0305 +/- 0.0494087 (limited) background = 82.8929 +/- 23.2866 slope = -0.0128031 +/- 0.0419971
peak_height = fit_peak_548.GetParameter(0)
peak_mean = fit_peak_548.GetParameter(1)
peak_sigma = fit_peak_548.GetParameter(2)
bin_width = hspec.GetBinWidth(1)
area_scale = ROOT.TMath.Sqrt(2 * ROOT.TMath.Pi()) / bin_width
peak_area = area_scale * peak_height * peak_sigma
fit_result = peak_fit_result.Get()
var_height = fit_result.CovMatrix(0, 0)
var_sigma = fit_result.CovMatrix(2, 2)
cov_height_sigma = fit_result.CovMatrix(0, 2)
area_variance = area_scale**2 * (
peak_sigma**2 * var_height
+ peak_height**2 * var_sigma
+ 2 * peak_height * peak_sigma * cov_height_sigma
)
area_error = ROOT.TMath.Sqrt(max(area_variance, 0.0))
print("Fit status:", int(peak_fit_result))
print("Mean:", peak_mean, "+/-", fit_peak_548.GetParError(1), "keV")
print("Sigma:", peak_sigma, "+/-", fit_peak_548.GetParError(2), "keV")
print("FWHM:", 2.35482 * peak_sigma, "+/-", 2.35482 * fit_peak_548.GetParError(2), "keV")
print("Fitted background p0, p1:", fit_peak_548.GetParameter(3), fit_peak_548.GetParameter(4))
print("Gaussian area:", peak_area, "+/-", area_error, "counts")
print("Cov(height, sigma):", cov_height_sigma)
Fit status: 0 Mean: 547.912850335577 +/- 0.06179746177912193 keV Sigma: 10.030549768698704 +/- 0.04940868351727179 keV FWHM: 23.62013920632708 +/- 0.11634855612014196 keV Fitted background p0, p1: 82.89291995584088 -0.012803089802491724 Gaussian area: 29983.16847672527 +/- 178.54626023362655 counts Cov(height, sigma): -1.2834784256546472
对固定 bin 宽度 \(\Delta E\)，拟合函数的 Gaussian height 为每 bin 的计数尺度，因此 Gaussian 成分的总面积为
\[N_{\rm peak}=\frac{\sqrt{2\pi}}{\Delta E}H\sigma.\]
height 与 sigma 来自同一次拟合，通常相关。令 \(K=\sqrt{2\pi}/\Delta E\)，用拟合 covariance matrix 传播误差：
\[\sigma_N^2=K^2\left(\sigma^2 C_{HH}+H^2 C_{\sigma\sigma}+2H\sigma C_{H\sigma}\right).\]
若只把 height 和 sigma 的误差分别平方相加，会漏掉 covariance 项。这里的面积是拟合所得 Gaussian 成分的积分，不包括线性本底。FWHM 只与 sigma 相差固定系数，误差也按同一系数换算。
残差（residual）是同一 bin 中观测计数减去拟合期望计数。由于拟合使用 I 选项，下面用函数在每个 bin 内的积分除以 bin 宽度计算期望值。
residual = ROOT.TGraph()
for b in range(hspec.FindBin(480), hspec.FindBin(620)):
low = hspec.GetBinLowEdge(b)
width = hspec.GetBinWidth(b)
expected = fit_peak_548.Integral(low, low + width) / width
residual.SetPoint(
residual.GetN(), hspec.GetBinCenter(b),
hspec.GetBinContent(b) - expected
)
residual.SetTitle("Fit residual;Energy (keV);Observed - fitted counts")
residual.SetMarkerStyle(20)
c1.Clear()
residual.Draw("AP")
c1.Draw()
4.4 坐标范围与对数坐标
SetRangeUser 改变显示的 x 范围；SetLogy() 开启对数 y 轴，SetLogy(0) 恢复线性 y 轴。SetLogx() 同理，但 x 范围需为正值，布拉格曲线中会用到。下面恢复完整能谱显示，不改变各 bin 的内容。
直方图被设置轴范围后，GetMean() 等统计量也可能按选定范围重新计算；要查看整幅能谱的统计量，应先恢复完整轴范围。
hspec.GetXaxis().SetRange(0, 0) # 恢复完整的 x 轴范围
c1.Clear()
c1.SetLogy()
hspec.Draw()
c1.Draw()
5. TH2：二维直方图¶
二维直方图统计成对数值，而不是单个数值。在探测器分析中，这两个量应属于同一次测量或同一个事例，例如束斑坐标 (x, y)、望远镜中两层探测器的信号，或粒子鉴别图中的两个观测量。
TH2 是二维直方图的基类。本例采用 TH2F，其 bin 内容用单精度浮点数保存。
5.1 探测器平面上的束斑¶
假设位置灵敏探测器给出每个粒子的 x、y 坐标。构造函数为
ROOT.TH2F(name, title, nx, xmin, xmax, ny, ymin, ymax)
x 轴在 xmin 到 xmax 之间划分 nx 个 bin，y 轴同理。下面两个轴都在 −30 到 30 mm 之间划分 120 个 bin，每个方向的 bin 宽度均为 0.5 mm。
c1.SetLogy(0)
hxy = ROOT.TH2F(
"hxy",
"Beam spot on detector plane;x (mm);y (mm)",
120, -30.0, 30.0,
120, -30.0, 30.0,
)
r2 = ROOT.TRandom3(13579)
for i in range(50000):
x_position = r2.Gaus(2.0, 6.0)
y_position = r2.Gaus(-1.0, 4.0)
hxy.Fill(x_position, y_position)
c1.Clear()
hxy.Draw("COLZ")
c1.Draw()
Gaus(2.0, 6.0) 生成均值 2 mm、sigma 为 6 mm 的 x；Gaus(-1.0, 4.0) 生成均值 −1 mm、sigma 为 4 mm 的 y。一次循环中的坐标对代表一个被探测的粒子，通过 Fill(x, y) 填入直方图。
Draw("COLZ") 中，COL 用颜色表示 bin 内容，Z 添加颜色标尺。这里 x、y 独立采样；TH2 显示的是联合分布，但画成二维图本身并不意味着两个量存在物理关联。
print(f"Entries = {hxy.GetEntries():.0f}")
print(f"Mean x = {hxy.GetMean(1):.3f} mm")
print(f"Mean y = {hxy.GetMean(2):.3f} mm")
print(f"Std dev x = {hxy.GetStdDev(1):.3f} mm")
print(f"Std dev y = {hxy.GetStdDev(2):.3f} mm")
Entries = 50000 Mean x = 2.023 mm Mean y = -1.006 mm Std dev x = 6.025 mm Std dev y = 3.984 mm
对于 TH2，轴编号 1 表示 x，2 表示 y。计算出的均值和标准差应接近生成位置分布时使用的参数。真实束斑分析会用实测位置代替随机样本，填充和显示操作相同。
5.2 选定 y 区间后作 x 投影
假设只关心穿过 −3 ≤ y < 1 mm 水平区域的粒子，查看它们在 x 方向上的分布。ProjectionX 沿 y 方向累加，得到 x 分布；ProjectionY 则沿 x 方向累加。
ProjectionX(name,first_y_bin,last_y_bin) 指定要累加的 y 轴 bin，两端均包含。这里边界与 bin 边缘重合，最后一个 bin 用 FindBin(1)-1。不指定范围时，默认累加整个 y 轴，包括 underflow 和 overflow。
px = hxy.ProjectionX("px")
first_y_bin = hxy.GetYaxis().FindBin(-3)
last_y_bin = hxy.GetYaxis().FindBin(1) - 1
px_gate = hxy.ProjectionX("px_gate", first_y_bin, last_y_bin)
all_counts = px.Integral(0, px.GetNbinsX() + 1)
gate_counts = px_gate.Integral(0, px_gate.GetNbinsX() + 1)
print("All events:", all_counts)
print("Events in y gate:", gate_counts)
print("Accepted fraction:", gate_counts / all_counts)
All events: 50000.0 Events in y gate: 19236.0 Accepted fraction: 0.38472
二维图上的红线给出选择边界；右图叠加全部事例和所选事例的 x 分布。SetLineColor 只改变显示颜色，真正的选择由前面的 y 轴 bin 范围完成。
cProjection = ROOT.TCanvas("cProjection", "Position gate and projection", 1100, 460)
cProjection.Divide(2, 1)
cProjection.cd(1)
hxy.Draw("COLZ")
gate_low = ROOT.TLine(-30, -3, 30, -3)
gate_high = ROOT.TLine(-30, 1, 30, 1)
for line in [gate_low, gate_high]:
line.SetLineColor(ROOT.kRed + 1)
line.Draw()
cProjection.cd(2)
px.SetTitle("x projection; x (mm);Counts / 0.5 mm")
px.SetLineColor(ROOT.kBlack)
px.Draw("HIST")
px_gate.SetLineColor(ROOT.kRed + 1)
px_gate.Draw("HIST SAME")
legend = ROOT.TLegend(0.60, 0.72, 0.88, 0.88)
legend.AddEntry(px, "All events", "l")
legend.AddEntry(px_gate, "y gate", "l")
legend.Draw()
cProjection.Draw()
本例 x、y 独立采样，因此 y gate 主要减少计数，所选 x 分布的中心和宽度基本不变。如果两个量存在关联，条件投影的形状也会改变。
固定电荷区间内作 PSD 投影，以及望远镜中选定某类事例后查看其他测量量，都采用这种分析思路。已有二维直方图的范围选择按 bin 进行；TTree 中的逐事例条件和 graphical cut 见 Tutorial II。
6. TFile：保存和读回 ROOT 对象¶
TFile 按 ROOT 内部名称保存图、直方图、函数等对象。写入对象不同于保存图片：ROOT 文件保留数值对象，之后可以重新读取和分析。
6.1 将对象写入 ROOT 文件¶
构造函数需要文件名和打开模式：
ROOT.TFile(filename, mode)
RECREATE 创建文件，若同名文件已存在则覆盖；READ 只读打开已有文件；UPDATE 在保留已有内容的情况下允许读写。使用 RECREATE 前先核对文件名。
fout.cd() 将 fout 设为当前 ROOT 目录，各对象的 Write() 按内部名称写入该目录；Close() 完成写入并关闭文件。
fout = ROOT.TFile("root_tutorial1.root", "RECREATE")
fout.cd()
gcal.Write()
hspec.Write()
hxy.Write()
px_gate.Write()
fout.Close()
文件保存 gcal、hspec、hxy 和条件投影 px_gate。ROOT 按对象的内部名称保存，名称来自构造函数或 SetName，不由 C++ 或 Python 变量名决定。
6.2 打开文件并查看内容
Open 打开刚才写出的文件，ls() 列出其中的对象。若打不开课程数据文件，先检查文件名和位置。
fin = ROOT.TFile.Open("root_tutorial1.root", "READ")
fin.ls()
TFile** root_tutorial1.root TFile* root_tutorial1.root KEY: TGraphErrors gcal;1 Energy calibration KEY: TH1F hspec;1 Simulated energy spectrum KEY: TH2F hxy;1 Beam spot on detector plane KEY: TH1D px_gate;1 Beam spot on detector plane
6.3 读回一个对象¶
Get(name) 读取以 name 保存的对象。从 ROOT 文件中读出的直方图通常与该文件关联，因此使用它时保持输入文件打开。
hread = fin.Get("hspec")
c1.cd()
c1.Clear()
hread.Draw()
c1.Draw()
这三个操作的作用不同：
Draw()在画布上显示对象。Write()将对象写入当前输出文件。Get()按 ROOT 对象名从文件中读取对象。
不再使用读出的对象后，再关闭文件。关闭后不要继续使用 hread，除非重新读取，或已明确解除它与文件的关联。
fin.Close()
小结¶
TF1 表示一维函数模型，TGraph 保存坐标对，TRandom3 通过 Gaus、Uniform、Poisson 等方法生成相应分布的伪随机数。TH1F、TH2F 分别累积一维、二维分布，TFile 用于保存 ROOT 对象。
后续分析中要注意这些区别：
- 函数模型不是逐个事例的数据集合；
- TGraph 保存数据点，直方图保存各 bin 的内容；
- 随机数发生器本身不是某一种概率分布；
- 显示对象与将对象写入文件是不同操作。