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}\frac{g_R}{g_L}\ln \frac{A_L}{A_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$ 的存在，数据点分布在一个很宽的范围内，无法直接看到清晰的直线。
In [2]:
TCanvas *c1 = new TCanvas;
TFile *f = new TFile("tree_hw.root");
TTree *tree = (TTree*)f->Get("tree");
In [3]:
tree->Draw("(tL+tR)/2.0 >> htof(200, 10, 90)", "pid==0");
c1->Draw();
In [4]:
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$¶

In [6]:
TCut center_cut = "pid==0 && abs(log(AL/AR)+0.5) < 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)$的关系¶

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

Left¶

In [8]:
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)", center_cut, "colz");
c2->cd(2);
TProfile *profL = walkL->ProfileX("profL");
profL->Fit(fFit, "RQ");
c2->Draw();
Double_t WL = fFit->GetParameter(1);
Double_t err_WL = fFit->GetParError(1);
Double_t CL = fFit->GetParameter(0);

Right¶

In [10]:
c2->Clear();
c2->Divide(2,1);
c2->cd(1);
tree->Draw("tR : 1.0/sqrt(AR) >> walkR(50, 0, 0.6)", center_cut, "colz");
c2->cd(2);
TProfile *profR = walkR->ProfileX("profR");
profR->Fit(fFit, "RQ");
c2->Draw();
Double_t WR = fFit->GetParameter(1);
Double_t err_WR = fFit->GetParError(1);
Double_t CR = fFit->GetParameter(0);
In [11]:
std::cout << Form(">>> 提取参数: WL=%.2f ± %.2f | WR=%.2f ± %.2f", WL, err_WL, WR, err_WR) << std::endl;

>>> 提取参数: WL=40.29 ± 0.19 | WR=45.15 ± 0.23

3. 应用修正¶

3.1 Walk系数修正¶

  • 蓝色峰： 修正前
  • 红色峰（修正后）： 峰宽显著变窄，形状接近对称的高斯分布。
In [13]:
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,30,70)", 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("h_old", "Uncorrected", "l");
leg->AddEntry("h_new", "Corrected", "l");
leg->Draw();

c3->Draw();

3.2 时间偏移修正¶

在完成 Time-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)$ 的分布理论上是一个矩形。但由于时间分辨率的存在，计数的分布在边缘处并非垂直下降，而是带有高斯卷积的阶梯。因此，导数的峰值位置才是物理边界最可靠的估计。

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

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

  • 刻度目标：识别左边界 $\Delta t_{left}$ 和右边界 $\Delta t_{right}$，建立线性映射：$x = k \cdot \Delta t + b$，使得重建出的 $x$ 恰好落在 $[-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)

3.2.3 位置重建验证¶

在完成 Time-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还需进一步修正。
In [16]:
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());
In [17]:
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();

3.3 TOF刻度¶

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

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

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

校准原理： $\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} $$

  • 图像分析（左图）： 绘制 $3.33 \cdot z - TOF_{meas}$ 的分布图并进行高斯拟合。拟合得到的均值 $mean \approx -26.27$ ns 即为该系统的绝对时间偏移量。
  • 图像分析（右图）： 应用偏移修正,归一化 $TOF$（$ns/m$）随位置 $x$ 的变化: 可以看到 $\gamma$ 射线带现在精固定在 $3.33$ ns/m 线上，且不再随位置 $x$ 发生倾斜或弯曲。
In [19]:
TString c_tof = Form("3.33* %s - %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<<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.2577

3.4. TOF刻度验证¶

  • 将所有入射粒子的测量 TOF 按照其实际飞行路径 $z$ 进行归一化。
  • 图像分析：
    • $\gamma$ 射线峰（左侧尖峰）： 位于 $3.33$ ns/m 处。修正后的 $\gamma$ 峰（红色）比修正前（蓝色）更宽，峰宽的大小直接反映了探测器的时间分辨率。
    • 中子分布（右侧宽峰）： 位于 $6$ ns/m 之后。由于中子质量大、速度慢，其飞行时间远大于 $\gamma$ 射线。
  • 结论： 经过 Walk 修正、位置刻度、偏移校准后的数据，与模拟中的飞行时间一致。
In [21]:
TCanvas *c6 = new TCanvas("tof","tof",1000,400);
c6->Divide(2,1);
c6->cd(1);
tree->Draw("tof_true/sqrt(x_true*x_true+5*5)>>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 的逐事件分析、筛选与二次存储¶

实验数据分析的一般流程：从原始数据（Raw Data）中提取信息，通过计算生成新的物理量，并根据特定的物理判据（如能量阈值）筛选有效事件，最终将结果保存到新的文件中。

语法¶

  • 1. 数据读取绑定 (SetBranchAddress)：

    tree->SetBranchAddress("BranchName", &Variable); 这是写入操作 Branch() 的逆过程。它将现有的数据内TTree分支(BranchName)与本地内存变量(Variable)的地址绑定，相当于在数据源与程序内存之间建立了一根专门的“数据管道”。

  • 2. 按需提取事件 (GetEntry)：

    tree->GetEntry(i); 的作用是加载第 i 个事件的数据。GetEntry 解压并读取那些通过 SetBranchAddress 挂载过的分支。执行后，绑定的本地变量会同步更新为该事件的真实数值。

示例代码：Analyze.C¶

#include <iostream>
#include <TFile.h>
#include <TTree.h>
#include <TMath.h>

void Analyze() {
    // 1. 输入配置：打开原始文件并绑定地址映射
    TFile *fIn = new TFile("tree_hw.root", "READ");
    TTree *tIn = (TTree*)fIn->Get("tree");

    // 内存工作区：声明本地变量接收原始数据
    Int_t raw_pid; Double_t raw_tL, raw_tR, raw_qL, raw_qR;
    tIn->SetBranchAddress("pid", &raw_pid);
    tIn->SetBranchAddress("tL",  &raw_tL);
    tIn->SetBranchAddress("tR",  &raw_tR);
    tIn->SetBranchAddress("qL",  &raw_qL);
    tIn->SetBranchAddress("qR",  &raw_qR);

    // 2. 输出配置：创建新文件与定标树
    TFile *fOut = new TFile("calibrated.root", "RECREATE");
    TTree *tOut = new TTree("tree_cal", "Filtered & Calibrated Data");

    // 定义定标物理量并建立输出分支绑定
    Double_t xt, xq;
    tOut->Branch("xt",  &xt,  "xt/D");      // 计算后的时间差位置
    tOut->Branch("xq",  &xq,  "xq/D");      // 计算后的电荷比对数
    tOut->Branch("pid", &raw_pid, "pid/I"); // 同时转存原始标识

    // 3. 事件循环：逐事件进行读取、定标与筛选
    Long64_t nEvents = tIn->GetEntries();
    std::cout << ">>> 开始处理并筛选事件，总计: " << nEvents << std::endl;

    for (Long64_t i = 0; i < nEvents; i++) {
        // 将当前事件的所有原始分支数据提取到本地变量中
        tIn->GetEntry(i); 

        // --- 物理筛选条件：仅记录两端电荷响应均大于 10 的有效事件 ---
        if (raw_qL > 10 && raw_qR > 10) {
            
            // 执行物理校准算法
            xt = raw_tL - raw_tR;
            xq = TMath::Log(raw_qL / raw_qR);

            // 【条件填充】：只有满足 if 条件的事件才会写入新树的缓冲区
            tOut->Fill(); 
        }

        // 实时进度反馈
        if (i % (nEvents / 20) == 0 || i == nEvents - 1) {
            std::cout << "\rProcessing: " << (int)((double)i/(nEvents-1)*100) << "%" << std::flush;
        }
    }

    // 4. 持久化存储与资源清理
    fOut->cd();      // 确保进入输出文件目录
    tOut->Write();   // 
    fOut->Close();   // 关闭并保存输出文件
    fIn->Close();    // 关闭输入文件
    std::cout << "\n>>> 处理完成，校准数据已保存至 calibrated.root" << std::endl;
}

作业：¶

  • 利用Graphic cut选出gamma和中子条件，保存为cutGamma.cc, cutNeutron.cc文件。
  • 在下面的分析中调用保存的Graphic cut来选择粒子种类，不得使用pid。
  • 下列分析中所需变量通过逐事件遍历生成，不得使用讲义的做法。
    • example：
TH1F *htL_corr = new TH1F("htl_corr","htl_corr", 200, 10, 200);
for (Long64_t i = 0; i < t->GetEntries(); i++) {
    t->GetEntry(i);
    if (qL <= 0 || qR <= 0) continue;
    Double_t tL_cor = tL - A_L / TMath::Sqrt(qL) - B_L;
    htl_corr->Fill(tL_corr);
    ...
}
  1. 用位置的多个切片，对时间信号进行 Time Walk 修正，得到两侧的平均walk系数，给出系数的误差。
  2. 利用微分法进行 位置刻度，刻度系数与模拟代码中的数值进行比较。
  3. 计算AL,AR两端信号刻度系数比值$g_L/g_R$ 以及闪烁体中光的衰减系数$\lambda$.
  4. 计算中子能量(En)与模拟数据进行比较。
  5. 将刻度后的事件信息(tL_cal,tR_cal,x_cal,tof_cal, En_cal) 新文件。
In [ ]: