核物理实验数据处理 Experimental Data Analysis in Nuclear Physics

1.2：Time-Walk 效应的校准与修正

1. 目的

1. timewalk修正方法，时间刻度方法：
2. 掌握TTree 逐事件处理

2. 原理

2.1 Time-Walk 的数学描述

信号幅度 $A$ 对触发时间 $t$ 的附加延迟在本例中用下式作经验描述： $$t_{measured} = t_{true} + \frac{W}{\sqrt{A}}$$ 当粒子击中长条闪烁体时，测量到的时间包含飞行时间、光传导延迟和 Walk 效应： $$t_L = TOF + \frac{L+x}{v_{sc}} + t_{0L} + \frac{W_{L}}{\sqrt{A_L}}$$ $$t_R = TOF + \frac{L-x}{v_{sc}} + t_{0R} + \frac{W_{R}}{\sqrt{A_R}}$$

$t$ 与入射位置 $x$ 、TOF 以及 A 三者同时相关。

2.2 切片法

利用两个物理约束来“冻结”无关变量：

1. 利用位置切片（Position Cut）： 限制入射位置在探测器的某一位置附近（利用$\log{A_L/A_R}$）。
2. 利用 Gamma 射线： 选出 $\gamma$ 的事件，此时在窄位置切片内，真实 TOF 近似相同。

• $ x_q = \frac{\lambda}{2}\ln \frac{A_R g_L}{A_L g_R} $

在上述两个条件下，测量时间简化为：$t_{L/R} \approx \text{Constant} + \frac{W_{{L/R}}}{\sqrt{A_{L/R}}}$。

---

2.3 观察 Walk 效应

• 上图： Gamma 峰非常弥散，且由于 Walk 效应存在严重的长尾（低能端延迟更大）。
• 下图： tL 随 1/sqrt(AL) 呈现弥散的带状。由于位置 $x$ 的存在，数据点分布在一个很宽的范围内，无法直接看到清晰的直线。

%jsroot on

TCanvas *c1 = new TCanvas;
TFile *f = new TFile("tree_hw_demo.root");
TTree *tree = (TTree*)f->Get("tree");

tree->Draw("(tL+tR)/2.0 >> htof(200, 10, 90)", "pid==0");
c1->Draw();

tree->Draw("tL : 1/sqrt(AL) >> h2(100,0,0.6, 100, 20, 80)", "pid==0", "colz");
c1->Draw();

3. 切片法 提取walk系数 $W_L, W_R$

TCut center_cut = "pid==0 && abs(log(AL/AR)-log(10.0/15.0)) < 0.05";
tree->Draw("log(AL/AR)>>hqxcut(100,-1.2,0.4)", center_cut, "");
((TH1F*)gPad->GetPrimitive("hqxcut"))->SetFillColor(kGreen);
tree->Draw("log(AL/AR)>>hqx(100,-1.2,0.4)", "pid==0", "same");
c1->Draw();

3.1 gamma 打在固定位置时 t～ 1/$\sqrt(A)$的关系

• 将两端的 $t : 1/\sqrt{A}$ 二维图保存到TProfile，利用线性方程 $Y = W_A \cdot X + C$ 拟合，得到walk系数 $W_A$。
• 可以把探测器切成多个独立的“切片（Slices）”，对每一个切片分别拟合求出斜率，最后用加权平均合并结果，从而利用更高统计量，以减少$W_A $的误差。

Left

TProfile 直接保存每个横轴 bin 内的事件平均时间及其误差，拟合斜率给出 W。左图用 TH2 显示分布，右图用相同事件直接填充 TProfile，避免把纵轴 bin 中心当作原始时间。位置切片有有限宽度，幅度比也有噪声，因此这是近似校准，可比较不同切片的结果。示例中的 pid 用于展示已知粒子的响应，作业改用图形选择。

拟合使用 Profile 的均值和均值误差。只有一个事件的 bin 无法从样本估计方差，故不把它作为加权拟合点；它仍保留在原始二维图和 Profile 中。下面用 TGraphErrors 明确列出参与拟合的均值点。

TF1 *fFit = new TF1("fFit", "pol1", 0.05, 0.6); 

TCanvas *c2 = new TCanvas("c2", "Calibration", 800, 400);
c2->Clear();
c2->Divide(2,1);
c2->cd(1);
tree->Draw("tL : 1.0/sqrt(AL) >> walkL(50,0,0.6,240,20,100)", center_cut, "colz");
c2->cd(2);
TProfile *profL = new TProfile("profL","Walk calibration;1/sqrt(AL);Mean time (ns)",50,0,0.6);
tree->Draw("tL:1/sqrt(AL)>>profL", center_cut, "prof");
TGraphErrors *meansL = new TGraphErrors;
for (int b=1; b<=profL->GetNbinsX(); ++b) {
    // 单个事件不能估计该 bin 的时间方差，因此不作为拟合点。
    if (profL->GetBinEntries(b)<2 || profL->GetBinError(b)<=0) continue;
    int n = meansL->GetN();
    meansL->SetPoint(n, profL->GetBinCenter(b), profL->GetBinContent(b));
    meansL->SetPointError(n, 0, profL->GetBinError(b));
}
meansL->Fit(fFit, "RQ");
fFit->Draw("same");
c2->Draw();
Double_t WL = fFit->GetParameter(1);
Double_t err_WL = fFit->GetParError(1);
Double_t CL = fFit->GetParameter(0);

Right

c2->Clear();
c2->Divide(2,1);
c2->cd(1);
tree->Draw("tR : 1.0/sqrt(AR) >> walkR(50,0,0.6,240,20,100)", center_cut, "colz");
c2->cd(2);
TProfile *profR = new TProfile("profR","Walk calibration;1/sqrt(AR);Mean time (ns)",50,0,0.6);
tree->Draw("tR:1/sqrt(AR)>>profR", center_cut, "prof");
TGraphErrors *meansR = new TGraphErrors;
for (int b=1; b<=profR->GetNbinsX(); ++b) {
    // 单个事件不能估计该 bin 的时间方差，因此不作为拟合点。
    if (profR->GetBinEntries(b)<2 || profR->GetBinError(b)<=0) continue;
    int n = meansR->GetN();
    meansR->SetPoint(n, profR->GetBinCenter(b), profR->GetBinContent(b));
    meansR->SetPointError(n, 0, profR->GetBinError(b));
}
meansR->Fit(fFit, "RQ");
fFit->Draw("same");
c2->Draw();
Double_t WR = fFit->GetParameter(1);
Double_t err_WR = fFit->GetParError(1);
Double_t CR = fFit->GetParameter(0);

std::cout << Form(">>> 提取参数: WL=%.2f ± %.2f | WR=%.2f ± %.2f", WL, err_WL, WR, err_WR) << std::endl;

>>> 提取参数: WL=39.85 ± 0.19 | WR=45.49 ± 0.23

3. 应用修正

3.1 Walk系数修正

• 蓝色峰： 修正前
• 红色峰（修正后）： 峰宽显著变窄，形状接近对称的高斯分布。

TCanvas *c3 = new TCanvas("c3", "Verification");
tree->SetLineColor(kBlue);
// 定义 Walk 修正后的时间
TString tL_corr = Form("(tL - %f/sqrt(AL))", WL);//"tL-45.2/sqrt(AL)"
TString tR_corr = Form("(tR - %f/sqrt(AR))", WR);
TString tof = Form("(%s + %s)/2.", tL_corr.Data(),tR_corr.Data());
TString tx = Form("(%s - %s)/2.", tL_corr.Data(),tR_corr.Data());

TString draw_tof = Form("%s >> h_new(200,10,90)", tof.Data());
tree->SetLineColor(kRed);
tree->Draw(draw_tof, "pid==0", "");

htof->SetLineColor(kBlue);
htof->Draw("same");

TLegend *leg = new TLegend(0.6, 0.7, 0.9, 0.9);
leg->AddEntry("htof", "Uncorrected", "l");
leg->AddEntry("h_new", "Corrected", "l");
leg->Draw();

c3->Draw();

3.2 时间偏移修正

在完成 Time-Walk 修正后，获得了反映粒子入射位置的已扣除主要 walk 的时间差信息。本节的目标是确定时间差与物理空间坐标 $x$ 的映射关系。

3.2.1 利用微分法进行位置刻度

探测器位置x与两侧时间差(修正后)之间有如下关系：

$$ x_{t} = \frac{v_{sc}}{2} (t_L - t_R - (t_{0L}-t_{0R})) $$

• 通过下列步骤得到刻度系数：$v_{sc}/{2}$ 和 $t_{0L}-t_{0R}$.

1. 物理思想：寻找“消失的边缘” 由于探测器在空间上是有限的（长度为 $2L$），当粒子均匀入射时，时间差 $\Delta t = (t_L - t_R)$ 的分布理论上是一个矩形。但由于时间分辨率的存在，计数的分布在边缘处并非垂直下降，而是带有高斯卷积的阶梯。因此，可用导数极值定位边缘，也可直接拟合展宽后的阶梯形状。

• 建立时间差的直方图(TH1) $h1$

• 微分法的逻辑：对 $h1$ 分布求导（数值微分），计数下降最快的地方（导数的极值点）即对应探测器的物理边界。

• 刻度目标：识别左边界 $\Delta t_{left}$ 和右边界 $\Delta t_{right}$，建立线性映射：$x = k \cdot \Delta t + b$，使两处边缘的中心映射到 −L 和 L；分辨会使少量重建值落在边界外。

• 获取 X 轴的 Bin 总数（不含 Underflow 和 Overflow）：

  Int_t nBins = h1->GetNbinsX();
  

  注：ROOT 的直方图索引从 1 到 nBins。0 号 Bin 是 Underflow（低于范围），nBins+1 号 Bin 是 Overflow（高于范围）。

• 在循环中遍历所有 Bin：

  for (int i = 1; i <= h1->GetNbinsX(); i++) {
      Double_t content = h1->GetBinContent(i);
      // ... 处理代码
  }
  

  • 数值微分：最简单的实现是 diff = h1->GetBinContent(i+1) - h1->GetBinContent(i)，将微分结果存入一个新的直方图。

  • 寻找“向下的峰”：在右边缘，导数是负的（计数在减少）。当你使用 gaus 拟合这个负峰时，必须手动给定初值，否则拟合会失败： , 需设定拟合参数的合理初值。

    TF1 *f1 = new TF1("f1","[0]*TMath::Exp(-0.5*((x-[1])/[2])^2)",xbin1,xbin2);//定义函数的方法 TF1
    //gaus:f(x) = p0*exp(-0.5*((x-p1)/p2)^2)
    f1->SetParameters(-A0_guess, x_guess, sigma_guess); // 振幅设为负(-A0_guess)，给定位置初值(x_guess)和宽度初值(sigma_guess)
    

有限差分应止于最后一个相邻的普通 bin，不把 overflow 混入。邻近差分共享计数，彼此相关；用导数图定位边缘时，不宜把各点简单当作独立计数来求精确误差。

3.2.3 位置重建验证

在完成 Time-Walk 修正后，获得了反映粒子入射位置的已扣除主要 walk 的时间差信息。

• 物理公式： 利用两侧光传导的时间差重建位置： $$x = \frac{V_{sc}}{2} \cdot \left( (t_{L,corr} - t_{R,corr}) - \Delta t_{0} \right)$$ 其中，$\Delta t_{0}$ 是通过微分法或几何中心确定的时间偏移。

• 飞行路径： 对于距离靶心垂直距离为 $D$ 的探测器，打在位置 $x$ 处的事件，其理论飞行距离为 $z = \sqrt{D^2 + x^2}$。

下面的分析直接采用分析代码中的参数设定：

const Double_t Vsc = 0.075; // m/ns
    const Double_t t0L_t0R = 5.5 - 20.4;   //ns

• 图像分析（左图）： 绘制重建后的位置 $x$ 的一维分布。可以看到在 $[-1.0, 1.0]$ 米范围内，计数分布接近平直的矩形。这说明探测器对均匀入射的粒子具有良好的空间响应，且物理边界（条长 2m）在刻度后得到了正确还原。
• 图像分析（右图）： 绘制单位距离的 TOF（tof/z）随位置 $x$ 的二维关联图。此时尚未进行绝对时间校准，可以看到 $\gamma$ 射线带（$pid==0$）分布在 $y \approx 8.5$ ns/m 附近。数值不对，而且进行飞行距离归一化的tof与 $x$ 存在依赖关系，说明tof还需进一步修正。

const Double_t D = 5.0;      // m, 靶到探测器距离
const Double_t L = 1.0;      // m, 探测器半长
const Double_t Vsc = 0.075;  // m/ns
const Double_t t0L_t0R = 5.5 - 20.4;    //ns  

//TString tof = Form("(%s + %s)/2.", tL_corr.Data(),tR_corr.Data());
//TString tx = Form("(%s - %s)/2.", tL_corr.Data(),tR_corr.Data());

TString x = Form("%f * ((%s - %s) - %f)/2.", Vsc, tL_corr.Data(), tR_corr.Data(), t0L_t0R);
TString z = Form("sqrt(pow(%f,2) + pow(%s,2))", D, x.Data());
TString tof_norm = Form("%s / %s", tof.Data(), z.Data());
TString draw_x_tof_norm = Form("%s:%s >> htofx(100,-1,1, 100,7,10)", tof_norm.Data(), x.Data());

TCanvas *c4 = new TCanvas("c4", "Calibrate", 1000, 500);
c4->Divide(2,1);
c4->cd(1);
tree->Draw(x, "pid==1", "colz");
c4->cd(2);
tree->Draw(draw_x_tof_norm, "pid==0 ", "colz");
c4->Draw();

位置刻度后的物理边界与 cut

本例塑料闪烁体长 2 m，真实入射位置在 $-L<x<L$，其中 $L=1$ m。有限时间分辨使重建位置 $x_{\rm rec}=x+\delta x$ 在边缘外也有少量分布。刻度应让边缘响应的中心对应 ±L，而不是把最外侧的非零 bin 拉回 ±L。

若边缘附近照明均匀、位置误差近似 Gaussian，则位置分布是矩形与 Gaussian 的卷积。下面用 $\sigma_x=0.03$ m 的假设值说明形状，不是对上面模拟结果的分辨拟合。左图画右边缘附近的测量分布；右图画给定真实位置时，通过 $|x_{\rm rec}|<L$ 的概率。

TCanvas *cBoundary = new TCanvas("cBoundary","Scintillator edge",1000,380);
cBoundary->Divide(2,1);
TF1 *edge = new TF1("edge",
    "0.5*(TMath::Erf((x+[0])/(sqrt(2)*[1]))-TMath::Erf((x-[0])/(sqrt(2)*[1])))",
    0.85,1.15);
edge->SetParameters(L,0.03); // L=1 m；sigma=0.03 m 为示意参数
cBoundary->cd(1);
edge->SetTitle("Uniform illumination, Gaussian response;x_{rec} (m);Relative density");
edge->Draw();
TLine *boundary = new TLine(L,0,L,1);
boundary->SetLineStyle(2); boundary->Draw();
cBoundary->cd(2);
TF1 *survival = (TF1*)edge->Clone("survival");
survival->SetRange(0.85,L);
survival->SetTitle("Acceptance of |x_{rec}|<L;True x (m);Cut survival probability");
survival->Draw();
cout << "At x=L: " << survival->Eval(L)
     << "; at x=L-sigma: " << survival->Eval(L-0.03) << endl;
cBoundary->Draw();

At x=L: 0.5; at x=L-sigma: 0.841345

在位置误差无偏、对称且另一侧边界足够远的条件下，真实位置恰在边缘处的事例约有一半通过这个 cut；这不表示整个探测器损失一半。总损失还取决于入射位置分布，有多少粒子靠近边缘。

怎样处理边缘外的重建位置

不必因为 $x_{\rm rec}$ 略超出 ±L 就判定事例没有物理意义。若目的是保留有效探测事例，可以结合已知分辨与本底保留相容的尾部；若分析明确要求 $|x_{\rm rec}|<L$，则把被 cut 去的有效事例计作位置选择损失，而不是探测器无响应。

研究内部响应时，也可借独立的位置参考选取 $|x_{\rm ref}|<L-d$，避开边缘。$d$ 应结合参考位置误差及边缘响应选择；不能用待测探测器自身的位置条件定义效率分母。模拟中可用真实位置核验，真实实验则需要独立径迹、准直器或位置扫描。

对同一入射参考样本，分别数出 $N_{\rm ref}$、有有效读数的 $N_{\rm valid}$ 和通过位置 cut 的 $N_{\rm pass}$：

$$\epsilon_{\rm read}=\frac{N_{\rm valid}}{N_{\rm ref}},\qquad \epsilon_{\rm cut}=\frac{N_{\rm pass}}{N_{\rm valid}},\qquad \epsilon_{\rm selected}=\frac{N_{\rm pass}}{N_{\rm ref}} =\epsilon_{\rm read}\epsilon_{\rm cut}.$$

这是分阶段计数的恒等关系，并不要求两个过程独立。如果总效率已经包含此 cut，就不再单独重复修正边缘损失。

3.3 TOF刻度

位置刻度完成后，由于电子学系统的固有延迟（如电缆长度、触发延迟），测量到的 $TOF$ 并非绝对物理时间。

$$ \text{TOF}_{rec} = \frac{t_L + t_R}{2} + c_{tof}$$

• $c_{tof} = -L/v_{sc} - (t_{0L}+t_{0R})/2$ 为待求刻度系数。

校准原理： $\gamma$ 射线以光速 $c$ 飞行。已知靶心到入射点的几何距离为 $z$，则理论飞行时间应为 $z/c$（约 $3.33 \times z$ ns）。测量值与理论值的差值即为全局时间偏移（Time Offset）。 $$ TOF_{\gamma}= 3.33 * z $$

$$ c_{tof} = 3.33*z - \frac{t_L + t_R}{2} $$

• 左图：对几何飞行时间减去已扣除 walk 的平均时间作图，拟合中心值作为时间偏移；数值见代码输出。
• 图像分析（右图）： 应用偏移修正,归一化 $TOF$（$ns/m$）随位置 $x$ 的变化: 可以看到 $\gamma$ 射线带现在位于在 $3.33$ ns/m 线上，且不再随位置 $x$ 发生倾斜或弯曲。

TString c_tof = Form("%s/0.299792458 - %s >> hconst(100, -30, -22)", z.Data(),tof.Data());
TCanvas *c5 = new TCanvas("c5", "tof Calibration", 800, 400);
c5->Divide(2,1);
c5->cd(1);
tree->Draw(c_tof,"pid==0");
TH1F *hconst = (TH1F*)gROOT->FindObject("hconst");
hconst->Fit("gaus","Q");
TF1 *fitFunc = hconst->GetFunction("gaus");
double mean = fitFunc->GetParameter(1);
cout << "time offset of tof = " << mean << " +/- " << fitFunc->GetParError(1) << " ns" << endl;

c5->cd(2);
TString tof_norm_corr = Form("(%s+%f)/%s", tof.Data(), mean, z.Data());
TString draw_x_tof_norm_corr = Form("%s:%s >> htofx(100,-1,1,100,1.5,5)", tof_norm_corr.Data(), x.Data());
tree->Draw(draw_x_tof_norm_corr, "pid==0 ", "colz");
c5->Draw();

time offset of tof = -26.2583 +/- 0.00151439 ns

3.4 TOF 刻度验证

左图蓝色为模拟真值，红色为重建后按估计路径归一化的 TOF。gamma 真值位于 \(1/c=3.33564\) ns/m，重建分布还含时间涨落、位置误差和未知反应深度的影响；不能把红色峰宽直接等同于单端时间分辨。中子速度较小，分布位于更长 TOF 处。

右图用两端幅度的几何平均观察 TOF—幅度关联。

TCanvas *c6 = new TCanvas("tof","tof",1000,400);
c6->Divide(2,1);
c6->cd(1);
tree->Draw("tof_true/path_true>>htof_all(200,0,15)","","");
TH1F *htof_all =(TH1F*) gROOT->FindObject("htof_all");
htof_all->SetLineColor(kBlue);

TString draw_tof_corr = Form("%s>>htof_all_corr(200,0,15)",tof_norm_corr.Data());
tree->Draw(draw_tof_corr,"","same");

c6->cd(2);
TString draw_tof_corr_q = Form("sqrt(AL*AR):%s>>htof_all_corr_q(200,0,15,200,0,700)",tof_norm_corr.Data());
tree->Draw(draw_tof_corr_q,"","colz");
c6->Draw();

4. TTree 的逐事件读取与二次存储

SetBranchAddress 将已有分支绑定到同类型变量；GetEntry(i) 读取第 i 个事件的活动分支。需要只读部分分支时，再用 SetBranchStatus 设置。

示例：Analyze.C

把前面拟合得到的 WL、WR 和 TOF offset 传给事件循环。这里的几何与增益取模拟设定，作业中再从数据确定。CloneTree(0) 先复制结构，Fill 时同时保存原量和新增物理量。

#include <TFile.h>
#include <TTree.h>
#include <cmath>
#include <iostream>
void Analyze(double WL, double WR, double timeOffset) {
    TFile input("tree_hw_demo.root");
    TTree *tin = input.Get<TTree>("tree");
    double AL, AR, tL, tR;
    tin->SetBranchAddress("AL", &AL);
    tin->SetBranchAddress("AR", &AR);
    tin->SetBranchAddress("tL", &tL);
    tin->SetBranchAddress("tR", &tR);

TFile output("calibrated.root", "RECREATE");
TTree *tout = tin-&gt;CloneTree(0); // 保留原有分支，只复制结构
double xt, xq, tof_cal;
Long64_t entry;
tout-&gt;Branch("source_entry", &amp;entry, "source_entry/L");
tout-&gt;Branch("xt", &amp;xt, "xt/D");
tout-&gt;Branch("xq", &amp;xq, "xq/D");
tout-&gt;Branch("tof_cal", &amp;tof_cal, "tof_cal/D");
for (entry=0; entry&lt;tin-&gt;GetEntries(); ++entry) {
    tin-&gt;GetEntry(entry);
    if (AL&lt;=0 || AR&lt;=0) continue;
    double tl = tL-WL/std::sqrt(AL);
    double tr = tR-WR/std::sqrt(AR);
    xt = 0.075/2 * (tl-tr-(5.5-20.4)); // m
    xq = 3.8/2 * std::log(AR*10/(AL*15)); // m
    tof_cal = (tl+tr)/2 + timeOffset; // ns
    tout-&gt;Fill();
}
std::cout &lt;&lt; "Input=" &lt;&lt; tin-&gt;GetEntries() &lt;&lt; ", output=" &lt;&lt; tout-&gt;GetEntries() &lt;&lt; '\n';
tout-&gt;Write();
}

载入后调用 Analyze(WL, WR, mean)。原始输入不修改，结果另存为 calibrated.root。

作业

1. 在 TOF—幅度图上选择 gamma 和中子区域，保存并应用 TCutG，不使用 pid 作分析条件。
2. 比较多个位置切片的 walk 系数及误差，修正两端时间。用探测器边缘确定位置刻度，与模拟设定比较。
3. 从两端幅度比求增益比与衰减长度，比较时间差位置和电荷比位置。改变相对电荷分辨，观察两种位置估计的变化。
4. 用 gamma 峰完成 TOF 绝对刻度。按估计路径归一到 1 m，再由 \(\beta=d/(ct)\)、\(E_n=(1/\sqrt{1-\beta^2}-1)m_nc^2\) 求中子动能，与真值比较。将原始量和新增的时间、位置、TOF、能量逐事件保存到新树中。

重建时只使用测量量；真值用于最后核验。几何路径的单位为 m，时间为 ns，能量为 MeV。

gROOT->ProcessLine(".L Analyze.C");
Analyze(WL, WR, mean);

Input=449201, output=449201