0. 交互式使用 ROOT¶
ROOT C++ kernel 通常已加载常用 ROOT 库。下面仍列出头文件，便于看清使用了哪些类，也方便将代码移到 macro 或编译程序中。标准 C++ 头文件提供终端输出、格式控制和数学运算。
#include <iostream>
#include <iomanip>
#include <TROOT.h>
#include <TCanvas.h>
#include <TMath.h>
#include <TF1.h>
#include <TGraph.h>
#include <TGraphErrors.h>
#include <TLine.h>
#include <TLegend.h>
#include <TStyle.h>
#include <TRandom3.h>
#include <TH1.h>
#include <TH2.h>
#include <TFile.h>
#include <TFitResultPtr.h>
std::cout << "ROOT version: " << gROOT->GetVersion() << std::endl;
ROOT version: 6.40.02
#include <TF1.h> 使编译器能够识别 TF1 的声明，其他 ROOT 头文件作用相同。gROOT 是 ROOT 的全局对象，这里只调用 GetVersion() 查看版本。
JSROOT 可在 Jupyter 中交互式显示 ROOT 图。%jsroot on 是 Jupyter 命令，不是 C++ 语句，通常只需运行一次。若当前安装不提供此命令，可跳过该单元格，使用环境已配置的显示方式。
%jsroot on
ROOT 将函数、图和直方图绘制在画布 TCanvas 上。下面使用的构造函数为
TCanvas(name, title, width, height)
name 是 ROOT 内部名称，同一会话中应避免重复；title 是画布标题；width 和 height 给出初始宽、高，单位为像素。
new TCanvas(...) 动态创建对象并返回地址。因此 c1 声明为 TCanvas*，即指向 TCanvas 的指针。通过指针调用方法用 ->，如 c1->Clear()。本教程重复使用同一画布。
SetOptStat(0) 隐藏自动统计框；需要的数值在相应代码中输出。
gStyle->SetOptStat(0);
TCanvas *c1 = new TCanvas("c1", "ROOT Tutorial I", 800, 520);
ROOT 的数学函数可通过 TMath 调用。TMath::Sqrt(x) 返回 x 的平方根，TMath::Pi() 返回 π。std::cout 将文字和数值写到单元格输出区，std::endl 用于换行。
double x0 = 2.0;
std::cout << "sqrt(2) = " << TMath::Sqrt(x0) << std::endl;
std::cout << "pi = " << TMath::Pi() << std::endl;
sqrt(2) = 1.41421 pi = 3.14159
类（class）定义对象类型，例如 TCanvas 和 TF1；c1 指向一个具体画布。对象提供的操作称为方法（method）。
C++ 常见的两种方法调用方式为：
object.Method(arguments); // object itself
pointer->Method(arguments); // pointer to an object
许多 ROOT 示例用 new 创建需要持续使用的分析对象，并保存其指针，因此经常出现 ->。普通局部对象使用 .，如后面创建的随机数发生器。C++ 语句以分号结束。
1. TF1：一维函数¶
TF1 表示单个自变量的数学函数，常用于能量刻度曲线、峰形和本底模型。它保存公式和参数值，不保存逐个测量事例。
本节使用的构造函数为
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 范围内。
TF1 *fcal = new 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 轴标题。fcal->Draw() 将函数画在当前画布上，c1->Draw() 在 Jupyter 中显示画布。这里直接给定参数，没有进行拟合。
Eval(value) 计算指定自变量处的函数值，可用于按已知刻度关系将 channel 换算为能量。
double channel_test = 1500.0;
double energy_test = fcal->Eval(channel_test);
std::cout << "Channel " << channel_test
<< " -> Energy = " << energy_test << " keV" << std::endl;
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 分别对应 [0]、[1]、[2]。下面在 400–700 keV 范围内绘制函数。
TF1 *fpeak = new 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 误差直接当作能量误差。
const int n_cal = 5;
double ch[n_cal] = {512.0, 1024.0, 1538.0, 2049.0, 2561.0};
double Eref[n_cal] = {122.1, 245.4, 367.8, 489.0, 612.3};
double dch[n_cal] = {1.0, 1.0, 2.0, 2.0, 3.0};
double dEref[n_cal] = {0, 0, 0, 0, 0};
TGraph(n, x, y) 保存 n 对坐标；TGraphErrors(n, x, y, ex, ey) 进一步保存两个坐标的误差。下面横轴为 channel，纵轴为参考能量，所以第四、第五个数组分别为 channel 误差和能量误差。
SetName 设置保存到 ROOT 文件时使用的名称，SetMarkerStyle(20) 选择实心圆点。Draw("AP") 中 A 表示坐标轴，P 表示数据点；TGraphErrors 同时显示误差棒。
TGraphErrors *gcal = new TGraphErrors(n_cal, ch, Eref, 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 误差折算到能量方向。
TF1 *flin = new TF1("flin", "pol1", 0.0, 3000.0);
flin->SetParNames("offset", "slope");
flin->SetParameters(0.0, 0.24);
TFitResultPtr 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 不会从数值中自动推断单位。
double fitted_offset = flin->GetParameter(0);
double fitted_slope = flin->GetParameter(1);
std::cout << "Fit status: " << int(calibration_fit) << std::endl;
std::cout << "Offset: " << fitted_offset << " +/- " << flin->GetParError(0) << " keV" << std::endl;
std::cout << "Slope: " << fitted_slope << " +/- " << flin->GetParError(1) << " keV/channel" << std::endl;
Fit status: 0 Offset: 0.0914136 +/- 0.307312 keV Slope: 0.239022 +/- 0.000258847 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}\)。
TGraphErrors *cal_residual = new TGraphErrors();
cal_residual->SetName("cal_residual");
for (int i = 0; i < n_cal; ++i) {
cal_residual->SetPoint(i, ch[i], Eref[i] - flin->Eval(ch[i]));
cal_residual->SetPointError(i, 0, TMath::Abs(fitted_slope) * dch[i]);
}
cal_residual->SetTitle("Calibration residual;Channel;Reference - fitted energy (keV)");
cal_residual->SetMarkerStyle(20);
c1->Clear();
cal_residual->Draw("AP");
TLine *zero_cal = new TLine(ch[0], 0, ch[n_cal-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() 检查读入的点数，再计算插值。
TGraph *table = new TGraph("calibration_points.txt", "%lg %lg");
std::cout << "Read points: " << table->GetN() << std::endl;
std::cout << "Interpolated energy: " << table->Eval(1000.0) << " keV" << std::endl;
std::cout << "Fitted energy: " << flin->Eval(1000.0) << " keV" << std::endl;
Read points: 5 Interpolated energy: 239.62 keV Fitted energy: 239.113 keV
TGraph::Eval 默认在相邻数据点之间做线性插值；TF1::Eval 则计算拟合函数值，所以两个结果不必相同。用于射程、阻止本领等表格时，同样先确认两列的物理量、单位和覆盖范围，再选择合理的插值。
3. TRandom3：随机采样
Gaus 生成 Gaussian 测量涨落，Uniform 生成均匀分布数值，Poisson 生成固定时长内的计数。这三种方法都可由同一个随机数发生器调用。
构造函数为
TRandom3(seed)
随机种子（seed）确定发生器的初始状态。每次用同一个正整数种子重新创建对象，会得到相同序列；换一个种子，则得到来自同一分布的另一组样本。
这里 rng 是普通 C++ 对象，不是指针，因此用 . 调用方法。Rndm() 不接收参数，每次返回一个 0–1 之间均匀分布的随机数。
TRandom3 rng(12345);
for (int i = 0; i < 5; ++i) {
std::cout << "sample " << i << ": " << rng.Rndm() << std::endl;
}
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。
double E_true = 548.0; // keV
double sigma_E = 8.0; // keV
for (int i = 0; i < 10; ++i) {
double E_measured = rng.Gaus(E_true, sigma_E);
std::cout << "measurement " << i << ": "
<< E_measured << " keV" << std::endl;
}
measurement 0: 537.891 keV measurement 1: 553.223 keV measurement 2: 541.708 keV measurement 3: 540.53 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 (int i = 0; i < 10; ++i) {
double E_background = rng.Uniform(0.0, 2000.0);
std::cout << "background sample " << i << ": "
<< E_background << " keV" << std::endl;
}
background sample 0: 747.842 keV background sample 1: 1307.14 keV background sample 2: 309.947 keV background sample 3: 1495.43 keV background sample 4: 1784.69 keV background sample 5: 1922.61 keV background sample 6: 53.5795 keV background sample 7: 16.7766 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)。
double mean_count = 10.0; // 平均计数率为 10 Hz，测量 1 s
for (int i = 0; i < 10; ++i) {
std::cout << rng.Poisson(mean_count) << std::endl;
}
18 13 10 12 11 9 6 9 15 11
Gaus 和 Uniform 返回能量等连续量，Poisson 返回离散的事例数。选哪种分布取决于实验观测量，而不是它们是否属于同一个发生器类。
4. TH1：一维直方图¶
直方图表示许多测量值或模拟值的分布。横轴划分为若干小区间，称为 bin。Fill(value) 找到 value 所在的 bin 并将计数加一；直方图保存各 bin 的计数，不保存完整的原始数值列表。
TH1 是一维直方图的基类。本例使用 TH1F，其中 F 表示 bin 内容用单精度浮点数保存；TH1D 则使用双精度。对于本例的计数规模，两者都够用。
4.1 构造简单能谱¶
构造函数为
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。这是一个简化能谱，用来练习随机采样与直方图操作。
TH1F *hspec = new TH1F(
"hspec",
"Simulated energy spectrum;Energy (keV);Counts per 5 keV",
400, 0.0, 2000.0
);
TRandom3 r1(24680);
for (int i = 0; i < 30000; ++i) {
hspec->Fill(r1.Gaus(548.0, 10.0));
}
for (int i = 0; i < 18000; ++i) {
hspec->Fill(r1.Gaus(1250.0, 18.0));
}
for (int i = 0; i < 30000; ++i) {
hspec->Fill(r1.Uniform(0.0, 2000.0));
}
c1->Clear();
hspec->Draw();
c1->Draw();
每个循环加入一种谱成分，每次生成一个能量并传给 Fill。循环次数控制该成分的事例总数，不是最高 bin 的计数，因为峰事例会分布在多个 bin 中。不带选项的 Draw() 使用 ROOT 默认的一维直方图样式。
4.2 查看直方图信息¶
GetEntries() 返回 Fill 的调用次数，包括 underflow 和 overflow。GetMean()、GetStdDev() 描述整体分布；GetNbinsX() 返回 x 轴普通 bin 的个数，GetBinWidth(1) 返回第一个 bin 的宽度。
std::cout << "Entries = " << hspec->GetEntries() << std::endl;
std::cout << "Number bins = " << hspec->GetNbinsX() << std::endl;
std::cout << "Bin width = " << hspec->GetXaxis()->GetBinWidth(1)
<< " keV" << std::endl;
std::cout << "Mean = " << hspec->GetMean() << " keV" << std::endl;
std::cout << "Std dev = " << hspec->GetStdDev() << " keV" << std::endl;
Entries = 78000 Number bins = 400 Bin width = 5 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 应避开峰尾和邻峰；真实能谱还需要检查线性本底是否足够。
TGraphErrors *sideband = new TGraphErrors();
for (int b = hspec->FindBin(480); b < hspec->FindBin(620); ++b) {
double x = hspec->GetBinCenter(b);
if (x < 510 || x >= 590) {
double y = hspec->GetBinContent(b);
int i = sideband->GetN();
sideband->SetPoint(i, x, y);
sideband->SetPointError(i, 0, TMath::Sqrt(y > 0 ? y : 1.0));
}
}
TF1 *background_line = new TF1("background_line", "pol1", 480, 620);
background_line->SetParameters(75, 0);
TFitResultPtr background_result = sideband->Fit(background_line, "RSQ0");
std::cout << "Background fit status: " << int(background_result) << std::endl;
std::cout << "p0: " << background_line->GetParameter(0) << std::endl;
std::cout << "p1: " << background_line->GetParameter(1) << " counts/bin/keV" << std::endl;
hspec->GetXaxis()->SetRangeUser(470, 630);
c1->Clear();
hspec->Draw("E");
sideband->SetMarkerStyle(20);
sideband->SetMarkerColor(kRed + 1);
sideband->Draw("P SAME");
background_line->SetLineColor(kGreen + 2);
background_line->Draw("SAME");
c1->Draw();
Background fit status: 0 p0: 92.7567 p1: -0.0327109 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 从与数据相差很远的位置开始搜索。
int peak_bin = hspec->FindBin(510);
for (int b = hspec->FindBin(510); b < hspec->FindBin(590); ++b) {
if (hspec->GetBinContent(b) > hspec->GetBinContent(peak_bin)) peak_bin = b;
}
double mean0 = hspec->GetBinCenter(peak_bin);
double height0 = hspec->GetBinContent(peak_bin) - background_line->Eval(mean0);
double sigma0 = 10.0;
TF1 *fit_peak_548 = new 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);
std::cout << "Initial height, mean, sigma: " << height0 << " " << mean0 << " " << sigma0 << std::endl;
std::cout << "Initial background p0, p1: " << background_line->GetParameter(0) << " " << background_line->GetParameter(1) << std::endl;
Initial height, mean, sigma: 5897.15 547.5 10 Initial background p0, p1: 92.7567 -0.0327109
Fit(..., "LIRS") 中：L 使用 binned Poisson likelihood；I 在每个 bin 内对模型积分取平均；R 使用函数定义的拟合范围；S 返回拟合结果及 covariance matrix。联合拟合会同时调整峰和本底参数，sideband 结果只是初值。
下图红线是完整的 Gaussian + linear background 模型，绿色虚线是联合拟合得到的本底分量。两条曲线都覆盖 sidebands 和 peak 区域。
TFitResultPtr peak_fit_result = hspec->Fit(fit_peak_548, "LIRS");
TF1 *fitted_background = new 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(kRed + 1);
fit_peak_548->Draw("SAME");
fitted_background->SetLineColor(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
double peak_height = fit_peak_548->GetParameter(0);
double peak_mean = fit_peak_548->GetParameter(1);
double peak_sigma = fit_peak_548->GetParameter(2);
double bin_width = hspec->GetBinWidth(1);
double area_scale = TMath::Sqrt(2*TMath::Pi()) / bin_width;
double peak_area = area_scale * peak_height * peak_sigma;
double var_height = peak_fit_result->CovMatrix(0, 0);
double var_sigma = peak_fit_result->CovMatrix(2, 2);
double cov_height_sigma = peak_fit_result->CovMatrix(0, 2);
double area_variance = area_scale*area_scale * (
peak_sigma*peak_sigma*var_height
+ peak_height*peak_height*var_sigma
+ 2*peak_height*peak_sigma*cov_height_sigma
);
double area_error = TMath::Sqrt(area_variance > 0 ? area_variance : 0.0);
std::cout << "Fit status: " << int(peak_fit_result) << std::endl;
std::cout << "Mean: " << peak_mean << " +/- " << fit_peak_548->GetParError(1) << " keV" << std::endl;
std::cout << "Sigma: " << peak_sigma << " +/- " << fit_peak_548->GetParError(2) << " keV" << std::endl;
std::cout << "FWHM: " << 2.35482*peak_sigma << " +/- " << 2.35482*fit_peak_548->GetParError(2) << " keV" << std::endl;
std::cout << "Fitted background p0, p1: " << fit_peak_548->GetParameter(3) << " " << fit_peak_548->GetParameter(4) << std::endl;
std::cout << "Gaussian area: " << peak_area << " +/- " << area_error << " counts" << std::endl;
std::cout << "Cov(height, sigma): " << cov_height_sigma << std::endl;
Fit status: 0 Mean: 547.913 +/- 0.0617975 keV Sigma: 10.0305 +/- 0.0494087 keV FWHM: 23.6201 +/- 0.116349 keV Fitted background p0, p1: 82.8929 -0.0128031 Gaussian area: 29983.2 +/- 178.546 counts Cov(height, sigma): -1.28348
对固定 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 宽度计算期望值。
TGraph *residual = new TGraph();
for (int b = hspec->FindBin(480); b < hspec->FindBin(620); ++b) {
double low = hspec->GetBinLowEdge(b);
double width = hspec->GetBinWidth(b);
double 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：二维直方图¶
二维直方图统计成对数值。在探测器分析中，两个量应属于同一次测量或同一个事例，例如束斑坐标、两个探测器的信号或粒子鉴别图中的观测量。
TH2 是二维直方图的基类。本例使用 TH2F，其 bin 内容用单精度浮点数保存。
5.1 探测器平面上的束斑¶
假设位置灵敏探测器给出每个粒子的 x、y 坐标。构造函数为
TH2F(name, title, nx, xmin, xmax, ny, ymin, ymax)
x 轴在 xmin 到 xmax 之间划分 nx 个 bin，y 轴同理。下面两个轴都在 −30 到 30 mm 之间划分 120 个 bin，宽度均为 0.5 mm。
c1->SetLogy(0);
TH2F *hxy = new TH2F(
"hxy",
"Beam spot on detector plane;x (mm);y (mm)",
120, -30.0, 30.0,
120, -30.0, 30.0
);
TRandom3 r2(13579);
for (int i = 0; i < 50000; ++i) {
double x_position = r2.Gaus(2.0, 6.0);
double 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。一次循环产生一个 (x, y) 坐标对，并用 Fill(x, y) 填入。
Draw("COLZ") 中，COL 用颜色表示 bin 内容，Z 添加颜色标尺。这里 x、y 独立采样；TH2 显示联合分布，使用二维图本身并不意味着两个量存在物理关联。
std::cout << "Entries = " << hxy->GetEntries() << std::endl;
std::cout << "Mean x = " << hxy->GetMean(1) << " mm" << std::endl;
std::cout << "Mean y = " << hxy->GetMean(2) << " mm" << std::endl;
std::cout << "Std dev x = " << hxy->GetStdDev(1) << " mm" << std::endl;
std::cout << "Std dev y = " << hxy->GetStdDev(2) << " mm" << std::endl;
Entries = 50000 Mean x = 2.02314 mm Mean y = -1.00589 mm Std dev x = 6.02503 mm Std dev y = 3.98422 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。
TH1D *px = hxy->ProjectionX("px");
int first_y_bin = hxy->GetYaxis()->FindBin(-3);
int last_y_bin = hxy->GetYaxis()->FindBin(1) - 1;
TH1D *px_gate = hxy->ProjectionX("px_gate", first_y_bin, last_y_bin);
double all_counts = px->Integral(0, px->GetNbinsX() + 1);
double gate_counts = px_gate->Integral(0, px_gate->GetNbinsX() + 1);
std::cout << "All events: " << all_counts << std::endl;
std::cout << "Events in y gate: " << gate_counts << std::endl;
std::cout << "Accepted fraction: " << gate_counts / all_counts << std::endl;
All events: 50000 Events in y gate: 19236 Accepted fraction: 0.38472
二维图上的红线给出选择边界；右图叠加全部事例和所选事例的 x 分布。SetLineColor 只改变显示颜色，真正的选择由前面的 y 轴 bin 范围完成。
TCanvas *cProjection = new TCanvas("cProjection", "Position gate and projection", 1100, 460);
cProjection->Divide(2, 1);
cProjection->cd(1);
hxy->Draw("COLZ");
TLine *gate_low = new TLine(-30, -3, 30, -3);
TLine *gate_high = new TLine(-30, 1, 30, 1);
gate_low->SetLineColor(kRed + 1);
gate_high->SetLineColor(kRed + 1);
gate_low->Draw();
gate_high->Draw();
cProjection->cd(2);
px->SetTitle("x projection;x (mm);Counts / 0.5 mm");
px->SetLineColor(kBlack);
px->Draw("HIST");
px_gate->SetLineColor(kRed + 1);
px_gate->Draw("HIST SAME");
TLegend *legend = new 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 文件¶
构造函数需要文件名和打开模式：
TFile(filename, mode)
RECREATE 创建文件并覆盖已有同名文件；READ 只读打开；UPDATE 在保留已有内容的情况下允许读写。
fout->cd() 将 fout 设为当前 ROOT 目录，各对象的 Write() 按内部名称写入该目录；Close() 完成写入并关闭文件。
TFile *fout = new 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() 列出其中的对象。若打不开课程数据文件，先检查文件名和位置。
TFile *fin = 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 读回一个对象
GetObject("hspec", hread) 按保存的 ROOT 名称读取直方图，并将指针赋给 hread。
TH1F *hread = nullptr;
fin->GetObject("hspec", hread);
c1->cd();
c1->Clear();
hread->Draw();
c1->Draw();
这三个操作的作用不同：
Draw()显示对象。Write()将对象写入当前输出文件。GetObject()按 ROOT 对象名读取对象。
不再使用读出的对象后，再关闭输入文件。关闭后不要继续使用 hread，除非重新读取，或已解除直方图与文件的关联。
if (fin) fin->Close();
小结¶
TF1 表示一维模型，TGraph 保存坐标对，TRandom3 通过 Gaus、Uniform、Poisson 生成相应分布的伪随机数。TH1F、TH2F 分别累积一维、二维分布，TFile 保存 ROOT 对象。
C++ 版本还用到了以下语法：
- 通过对象调用方法用
.，通过对象指针调用用->； new动态创建对象并返回指针；::访问命名空间或类的静态方法，如TMath::Pi()和TFile::Open()；- 分号结束语句，花括号划定循环、条件判断或函数的代码块。