2.2 PPAC 信号处理-II
1. Tracking
利用多个 PPAC 的位置 (x,y,z) 信息进行 tracking 的原理如下图所示。由于 PPAC 的探测效率通常小于 100%，因此对于每个入射粒子，往往只有部分 PPAC 能够在 x 或 y 方向给出有效的位置信息。
假设对于某一入射粒子，有两个或两个以上探测器给出了有效的测量位置信息（如图中的红点所示，例如 1Ax, 2Ax 等），入射粒子在 x-z 或 y-z 平面上的飞行径迹即可由这些有效点进行线性拟合得到。
得到径迹的线性方程后，束流线上其他任意位置 (z) 处的 (x,y) 坐标，均可利用上述方程通过内插或外推计算得出。图中黑点即代表由线性方程计算/拟合得到的位置（它代表了我们所能推测出的粒子实际飞行路径）。
一般重建策略与本例的选择
在实际实验中，若强制要求单层探测器必须同时给出 $(x,y)$ 二维有效信号才参与拟合，会导致严重的整体探测效率损失。因此可以分别收集有效的 $x$ 和 $y$ 测量点，重建两个投影。下面的演示固定选择 1A、2A、3 在两个方向均有效的样本；作业再实现按方向选择有效点。即：只要某层探测器测得了有效的 $x$ 坐标，即便其 $y$ 信号缺失，该点仍独立参与 $x-z$ 平面的径迹拟合。
2. 径迹有效条件 对于 $x$ 或 $y$ 径迹，必须满足以下逻辑条件之一才被视为有效重建：
情况一： 最靠近下游的 F8PPAC3 具有有效信号
- 同时要求上游的 F8PPAC1 和 F8PPAC2 中，至少有 1 层 PPAC 提供有效信号。
情况二： 最靠近下游的 F8PPAC3 无有效信号
- 则完全依赖上游探测器，需满足以下 a) 或 b)：
- a) F8PPAC1 和 F8PPAC2 中，合计有 3 层或 3 层以上的 PPAC 提供有效信号。
- b) 若合计仅有 2 层 PPAC 提供有效信号，则要求这两层必须分别来自 F8PPAC1 和 F8PPAC2。 (若仅有的2层信号来自同一位置（如都在 F8PPAC1 内），因两点间距过近，向靶点外推时会产生巨大的角度误差)
- 则完全依赖上游探测器，需满足以下 a) 或 b)：
束流径迹重建结果（所有PPAC的信息都有x/y位置信息）
2. ROOT 文件中 TTree 的 Branch 说明
本实验使用的ROOT文件中，已完成位置刻度的 PPAC 数据存储在二维数组 PPACF8[i][j] 中。
- $z$ 坐标：代表各探测器平面在束流线上的绝对 $z$ 方向位置，由实验室几何测量事先给定。
- $x/y/z$坐标：单位mm。
- 当出现位置计算超界、信号堆积 (Pile-up) 或缺失时，数据统一用
-999或-1000填充。
数据映射关系表（从原始数据到物理分析变量）：
| 探测器层级 | 原始 Branch | 物理意义 (已刻度) | 讲义代码对应变量 |
|---|---|---|---|
| PPAC 1 Layer A | PPACF8[0][0]PPACF8[0][1]PPACF8[0][2]PPACF8[0][3] |
$X$ 坐标 (mm) $Y$ 坐标 (mm) $X$ 平面对应的 $Z$ 坐标 $Y$ 平面对应的 $Z$ 坐标 |
xx[0]yy[0]xz[0]yz[0] |
| PPAC 1 Layer B | PPACF8[1][0~4] |
(同上结构，含 Anode 时间) | (本示例暂未使用) |
| PPAC 2 Layer A | PPACF8[2][0~4] |
(同上结构) | xx[1], yy[1]xz[1], yz[1] |
| PPAC 2 Layer B (作为待测层 DUT) |
PPACF8[3][0~4] |
(同上结构) | xx2b[0], yy2b[0]xz2b, yz2banode2b |
| PPAC 3 | PPACF8[4][0~4] |
(同上结构) | xx[2], yy[2]xz[2], yz[2] |
3. 径迹重建 (Tracking) 事例代码
下面示例代码 假设 1A, 2A, 3 这三个探测器均给出了有效位置信息，利用这三点进行直线拟合，并将径迹外推，用来评估 2B 探测器的各项探测效率及系统分辨率。
算法步骤如下：
- 1. 空间独立直线拟合：
- $x-z$ 平面拟合直线方程 $x = f_x(z)$，使用已知数据点：
(xx[0],xz[0]), (xx[1],xz[1]), (xx[2],xz[2]) - $y-z$ 平面拟合直线方程 $y = f_y(z)$，使用已知数据点：
(yy[0],yz[0]), (yy[1],yz[1]), (yy[2],yz[2])
- $x-z$ 平面拟合直线方程 $x = f_x(z)$，使用已知数据点：
- 2. included residual 计算（先检查径迹与测量的一致性）：
- 测量值与拟合值的差：$dx[i] = xx[i] - f_x(xz[i])$ ； $dy[i] = yy[i] - f_y(yz[i])$ $(i=0,1,2)$
- 3. 待测层 F8PPAC2_B (2B) 的对比：
- 实验测量坐标：
(xx2b[0], xz2b)和(yy2b[0], yz2b) - 外推预期坐标： 利用拟合直线算得
(xx2b[1], xz2b)和(yy2b[1], yz2b) - 阳极触发信号：
anode2b（用于统计参考样本内的阳极响应效率）
- 实验测量坐标：
- 4. 物理靶位坐标外推：
- 利用拟合直线外推至靶点平面（假设靶位于 $z=0$），求得靶上坐标
(tx, ty)。
- 利用拟合直线外推至靶点平面（假设靶位于 $z=0$），求得靶上坐标
%jsroot on
靶点位置与入射方向的不确定度
外推精度由单层测量误差、探测器间距以及外推距离共同决定。三个测量点不一定比相距较远的两个点有更好的角度精度。这里先假定每层位置误差为 1 mm，说明 ROOT 如何把测量误差传播到靶点与方向；这个数值不是本数据的分辨测量结果。
2. 协方差矩阵 (Covariance Matrix)
假设在 $X-Z$ 平面内用线性方程 $x(z) = p_0 + p_1 \cdot z$ 拟合粒子径迹。 最小二乘法在给出截距 $p_0$ 和斜率 $p_1$ 的同时，必然会给出一个 $2 \times 2$ 的协方差矩阵 $C$：
$$ C = \begin{pmatrix} \sigma_{p_0}^2 & Cov(p_0, p_1) \\ Cov(p_1, p_0) & \sigma_{p_1}^2 \end{pmatrix} $$
- 对角线元素： 分别是截距和斜率自身误差的平方（方差）。
slope 与 intercept 的 covariance 取决于 z 原点及探测器布局，符号不总为负。外推误差使用完整协方差矩阵，不能预先假定忽略交叉项一定使误差偏大或偏小。
3. 原理推导：外推靶点 $(t_x)$ 与入射方向的投影角 $(\theta)$ 的误差传递
假设物理靶位于 $z = z_{target}$ 处。
A. 靶上位置 $t_x$ 的误差推导（线性传递）
根据直线方程，外推点的计算公式为：$t_x = p_0 + p_1 \cdot z_{target}$。 利用多元函数误差传递公式： $$ \sigma_{f(p_0, p_1)}^2 = \left( \frac{\partial f}{\partial p_0} \right)^2 \sigma_{p_0}^2 + \left( \frac{\partial f}{\partial p_1} \right)^2 \sigma_{p_1}^2 + 2 \left( \frac{\partial f}{\partial p_0} \right) \left( \frac{\partial f}{\partial p_1} \right) Cov(p_0, p_1) $$
代入偏导数（$\frac{\partial t_x}{\partial p_0} = 1$, $\frac{\partial t_x}{\partial p_1} = z_{target}$），得到模型内的靶点位置方差： $$ \mathbf{\sigma_{t_x}^2 = \sigma_{p_0}^2 + (z_{target})^2 \cdot \sigma_{p_1}^2 + 2 \cdot z_{target} \cdot Cov(p_0, p_1)} $$ (注：Y 方向的 $t_y$ 误差推导与 X 方向完全相同，只需将 X 平面的拟合参数代入即可。)
B. 物理入射方向的投影角 $\theta$ 的误差推导（非线性传递）
拟合出的斜率 $p_1$ 等于角度的正切值，真正的物理角度为：$\theta = \arctan(p_1)$。 利用一元非线性误差传递公式： $$ \sigma_\theta^2 = \left( \frac{d\theta}{dp_1} \right)^2 \sigma_{p_1}^2 $$ 根据微积分导数公式 $\frac{d}{dx}\arctan(x) = \frac{1}{1+x^2}$，一阶误差传播得到投影角的标准不确定度： $$ \mathbf{\sigma_\theta = \frac{\sigma_{p_1}}{1 + p_1^2}} $$
\(\arctan k_x\)、\(\arctan k_y\) 是两个投影角；相对 z 轴的极角为 \(\theta=\arctan\sqrt{k_x^2+k_y^2}\)。两点也可确定直线并在已知测量误差下传播参数误差，但 NDF=0，不能用 χ²/NDF 检验拟合。
束斑中不同事例的真实位置本来就有分布，束斑中心不应仅按各事例的 tracking 误差加权；inverse-variance 平均要求它们测量同一个参数，或在模型中同时纳入束斑的本征宽度。
用外推误差选择径迹
逐事件保存外推位置和方向的不确定度，是为了让后续分析按所需精度选择径迹，而不是只按参与拟合的探测器数量分类。例如，若分析允许靶点两方向各有 1.5 mm 的标准不确定度，可对输出树使用：
TCut goodTrack = "sigma_tx>0 && sigma_ty>0 && sigma_tx<1.5 && sigma_ty<1.5";
tree->Draw("ty:tx", goodTrack, "colz");这里的 1.5 mm 是精度要求的示例，不是固定的 PPAC cut；当前误差还依赖前面假定的单层分辨。比较选择前后的靶点与角度分布，可以检查精度选择是否同时改变了接受度。
4. 位置分辨：先区分两种 residual
准直源或 mask 可用于独立测量。若测量位置是真实照明位置与独立读出误差之和，则方差相加： $\sigma_{\rm measured}^2=\sigma_{\rm illumination}^2+\sigma_{\rm det}^2$。 均匀照明宽度为 $w$ 的狭缝，其位置方差是 $w^2/12$，不是 Gaussian 峰宽。拟合时可把已知照明分布与分辨函数卷积。源、束流的电离密度及边缘散射不同，所得分辨也可能不同。
束流实验可用其他 PPAC 预测待测层（DUT）的位置。定义 $r_i=x_i-\hat x_{-i}(z_i)$，其中下标 $-i$ 表示拟合中没有使用第 i 层。若 DUT 与参考测量误差独立，且直线模型足够描述径迹，
$$\sigma_{r_i}^2=\sigma_i^2+\sigma_{\rm pred}^2(z_i),\qquad \sigma_{\rm pred}^2(z_i)=C_{00}+2z_i C_{01}+z_i^2 C_{11}.$$
因此，residual 的宽度还要扣除 reference track 的预测方差，才能估计单层分辨。参考层分辨若未知，可结合多层的 excluded residual 建立方差方程；不能把演示中假设的 1 mm 当作已经测出的输入。多重散射、对准误差和读出相关性也会影响这个关系。
相反，若第 i 层已经参与拟合，$x_i$ 与预测位置相关，included residual 通常较窄，不能套用上述方差相加公式。下面先画 included residual 检查径迹一致性，再画未参与拟合的 2B residual。
5. 探测效率：分母由谁提供？
用阳极有效的事例作分母，得到的是给定阳极响应后的位置读出效率： $\epsilon_{x|a}=N_{x\cap a}/N_a$。它不计入阳极自身漏掉的事例。
用独立 reference tracks 作分母，可估计该触发及照明样本内的响应效率。先要求参考径迹通过 DUT 的内部有效区域，再统计 DUT 的 anode、x、y 与 x-y 是否有效。参考径迹及分母的 cut 不使用 DUT 信号；x-y 联合效率直接计数，不假定 x 与 y 独立。
本文件中的位置已在上游处理阶段作过有效性选择。因此，下文“有效读数”指文件中仍保留的有效坐标，不能据此拆分上游已丢失的电子学响应与 time-sum 选择损失。若要分开测量这些效率，应从原始信号逐级计数。
6. 分析代码框架
利用 ROOT 的 MakeClass 功能，生成用于遍历 TTree 的 C++ 框架代码：
# 在终端中启动 ROOT 并加载数据文件
root -l f8ppac001.root
root [1] tree->MakeClass("tracking");
root [2] .q
提示：执行后会生成 tracking.h 和 tracking.C，按照下文修改这两个文件。
修改后运行分析的命令：
root -l
root [1] .L tracking.C
root [2] tracking t;
root [3] t.Loop();
root [4] .q
运行结束后，生成 tracking_demo.root 。
tracking.h
先对本节的 f8ppac001.root 运行 MakeClass，再在 public: 后补充这些成员。生成的分支绑定、构造函数和 LoadTree 等内容保持原样。
Double_t xx[3], xz[3], yy[3], yz[3], dx[3], dy[3];
Double_t xx2b[2], yy2b[2], xz2b, yz2b, anode2b;
Double_t tx,ty,theta_x,theta_y,sigma_tx,sigma_ty,sigma_thetax,sigma_thetay,c2nx,c2ny;
Long64_t source_entry;
void SetBranch(TTree *tree);
void TrackInit();
void SetTrace(TH2D *h, Double_t k, Double_t b, Int_t min, Int_t max);头文件需要 #include <TH2.h>。完整文件与下面的 tracking.C 放在同一目录。
tracking.C：逐事件计算
#define tracking_cxx
#include "tracking.h"
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TF1.h>
#include <TGraphErrors.h> // 输入测量点及其误差
#include <TFitResult.h>
#include <TMatrixDSym.h> // 拟合参数的协方差矩阵
#include <iostream>
#include <cmath>
using namespace std;
void tracking::SetBranch(TTree *tree)
{
tree->Branch("source_entry", &source_entry, "source_entry/L");
tree->Branch("xx", xx, "xx[3]/D");
tree->Branch("xz", xz, "xz[3]/D");
tree->Branch("yy", yy, "yy[3]/D");
tree->Branch("yz", yz, "yz[3]/D");
tree->Branch("dx", dx, "dx[3]/D");
tree->Branch("dy", dy, "dy[3]/D");
tree->Branch("xx2b", xx2b, "xx2b[2]/D");
tree->Branch("yy2b", yy2b, "yy2b[2]/D");
tree->Branch("anode2b", &anode2b, "anode2b/D");
// 保存位置、投影角与模型内传播的不确定度
tree->Branch("tx", &tx, "tx/D");
tree->Branch("ty", &ty, "ty/D");
tree->Branch("theta_x", &theta_x, "theta_x/D");
tree->Branch("theta_y", &theta_y, "theta_y/D");
tree->Branch("sigma_tx", &sigma_tx, "sigma_tx/D");
tree->Branch("sigma_ty", &sigma_ty, "sigma_ty/D");
tree->Branch("sigma_thetax", &sigma_thetax, "sigma_thetax/D");
tree->Branch("sigma_thetay", &sigma_thetay, "sigma_thetay/D");
tree->Branch("c2nx", &c2nx, "c2nx/D");
tree->Branch("c2ny", &c2ny, "c2ny/D");
tree->Branch("beamTrig", &beamTrig, "beamTrig/I");
tree->Branch("must2Trig", &must2Trig, "must2Trig/I");
tree->Branch("targetX", &targetX, "targetX/F");
tree->Branch("targetY", &targetY, "targetY/F");
}
void tracking::TrackInit()
{
// 初始化所有计算变量为无效值，防止上一个事件的数据污染
tx = -999; ty = -999;
c2nx = -1; c2ny = -1;
for (int i=0;i<3;++i) { dx[i]=-999; dy[i]=-999; }
theta_x = -999; theta_y = -999;
sigma_tx = -1; sigma_ty = -1; // 误差初始化为负数代表无效
sigma_thetax = -1; sigma_thetay = -1;
xx[0] = PPACF8[0][0]; yy[0] = PPACF8[0][1]; xz[0] = PPACF8[0][2]; yz[0] = PPACF8[0][3];
xx[1] = PPACF8[2][0]; yy[1] = PPACF8[2][1]; xz[1] = PPACF8[2][2]; yz[1] = PPACF8[2][3];
xx[2] = PPACF8[4][0]; yy[2] = PPACF8[4][1]; xz[2] = PPACF8[4][2]; yz[2] = PPACF8[4][3];
xx2b[0] = PPACF8[3][0]; yy2b[0] = PPACF8[3][1];
xz2b = PPACF8[3][2]; yz2b = PPACF8[3][3];
anode2b = PPACF8[3][4];
xx2b[1] = -1000; yy2b[1] = -1000;
}
void tracking::SetTrace(TH2D *h, Double_t k, Double_t b, Int_t min, Int_t max){
if(h == 0 || min >= max) return;
for(int i = min; i < max; i++){
h->Fill(i, i * k + b);
}
}
void tracking::Loop()
{
if (fChain == 0) return;
TFile *opf = new TFile("tracking_demo.root", "RECREATE");
TTree *tree = new TTree("tree", "PPAC Tracking Data");
SetBranch(tree);
TH2D *htf8xz = new TH2D("htf8xz", "X-Z Plane Trace; Z (mm); X (mm)", 2200, -2000, 200, 300, -150, 150);
TH2D *htf8yz = new TH2D("htf8yz", "Y-Z Plane Trace; Z (mm); Y (mm)", 2200, -2000, 200, 300, -150, 150);
// 输入假设的单层位置误差，演示误差传播
TGraphErrors *grx = new TGraphErrors(3);
TGraphErrors *gry = new TGraphErrors(3);
TF1 *fx = new TF1("fx", "pol1", -2000, 0);
TF1 *fy = new TF1("fy", "pol1", -2000, 0);
// 假设：所有PPAC每层的本征位置分辨率为 1.0 mm
const double det_resolution = 1.0;
const double z_target = 0.0; // 物理靶所在Z坐标位置
Long64_t nentries = fChain->GetEntriesFast();
Long64_t nbytes = 0, nb = 0;
for (Long64_t jentry = 0; jentry < nentries; jentry++) {
Long64_t ientry = LoadTree(jentry);
if (ientry < 0) break;
nb = fChain->GetEntry(jentry); nbytes += nb;
source_entry = jentry;
TrackInit();
bool b1a = abs(xx[0]) < 150 && abs(yy[0]) < 150;
bool b2a = abs(xx[1]) < 150 && abs(yy[1]) < 150;
bool b3 = abs(xx[2]) < 100 && abs(yy[2]) < 100;
if(!b1a || !b2a || !b3) continue;
// ================= X-Z 平面径迹拟合与误差计算 =================
for(int i=0; i<3; i++) {
grx->SetPoint(i, xz[i], xx[i]);
grx->SetPointError(i, 0.0, det_resolution); // 关键：输入Z和X的误差
}
// S 保存结果及协方差；Q 减少日志；N 不附加逐事件拟合函数
TFitResultPtr rx = grx->Fit(fx, "SQN");
if (int(rx)==0 && rx.Get() && rx->IsValid()) {
double p0_x = fx->GetParameter(0);
double p1_x = fx->GetParameter(1);
// 提取中心值
xx2b[1] = fx->Eval(xz2b);
tx = p0_x + p1_x * z_target;
theta_x = atan(p1_x); // 物理出射角 (rad)
// 提取协方差矩阵并计算严谨物理误差
TMatrixDSym cov_x = rx->GetCovarianceMatrix();
double var_p0 = cov_x(0, 0);
double var_p1 = cov_x(1, 1);
double cov_p0_p1 = cov_x(0, 1);
// 外推位置误差传递公式
double err2_tx = var_p0 + (z_target * z_target * var_p1) + (2.0 * z_target * cov_p0_p1);
sigma_tx = sqrt(err2_tx);
// 角度非线性误差传递公式
sigma_thetax = sqrt(var_p1) / (1.0 + p1_x * p1_x);
c2nx = rx->Chi2() / rx->Ndf();
if (jentry < 10000) SetTrace(htf8xz, p1_x, p0_x, -1800, 0);
for(int i=0; i<3; i++) dx[i] = xx[i] - fx->Eval(xz[i]);
}
// ================= Y-Z 平面径迹拟合与误差计算 =================
for(int i=0; i<3; i++) {
gry->SetPoint(i, yz[i], yy[i]);
gry->SetPointError(i, 0.0, det_resolution);
}
TFitResultPtr ry = gry->Fit(fy, "SQN");
if (int(ry)==0 && ry.Get() && ry->IsValid()) {
double p0_y = fy->GetParameter(0);
double p1_y = fy->GetParameter(1);
yy2b[1] = fy->Eval(yz2b);
ty = p0_y + p1_y * z_target;
theta_y = atan(p1_y);
TMatrixDSym cov_y = ry->GetCovarianceMatrix();
double var_p0 = cov_y(0, 0);
double var_p1 = cov_y(1, 1);
double cov_p0_p1 = cov_y(0, 1);
double err2_ty = var_p0 + (z_target * z_target * var_p1) + (2.0 * z_target * cov_p0_p1);
sigma_ty = sqrt(err2_ty);
sigma_thetay = sqrt(var_p1) / (1.0 + p1_y * p1_y);
c2ny = ry->Chi2() / ry->Ndf();
if (jentry < 10000) SetTrace(htf8yz, p1_y, p0_y, -1800, 0);
for(int i=0; i<3; i++) dy[i] = yy[i] - fy->Eval(yz[i]);
}
// 将本事件结果写入 Tree (包括新算出的误差)
if (c2nx>=0 && c2ny>=0) tree->Fill();
if(jentry % 10000 == 0) cout << "Processing Event: " << jentry << " / " << nentries << endl;
}
// 释放内存并保存结果
delete grx; delete gry;
delete fx; delete fy;
htf8xz->Write();
htf8yz->Write();
tree->Write();
cout << "Input events=" << nentries << ", accepted reference tracks=" << tree->GetEntries() << endl;
opf->Close();
cout << "Tracking finished. Output saved to tracking_demo.root" << endl;
}
本例固定使用 1A、2A、3 的有效测量，2B 不参与拟合。SetPointError(i,0,1) 输入假设的 1 mm 位置误差；S 返回结果和协方差，Q 减少逐事件日志，N 不向图附加每次拟合的函数。轨迹图只累积输入前 10000 个事例中的合格径迹，结果树保留全部合格事例。
gROOT->ProcessLine(".L tracking.C");
gROOT->ProcessLine("{ TFile *input=new TFile(\"f8ppac001.root\"); tracking tr(input->Get<TTree>(\"tree\")); tr.Loop(); }");
TFile *f = new TFile("tracking_demo.root");
TTree *tree = (TTree*)f->Get("tree");
TCanvas *c1 = new TCanvas("c1","Tracking");
Processing Event: 30000 / 739685 Processing Event: 50000 / 739685 Processing Event: 110000 / 739685 Processing Event: 130000 / 739685 Processing Event: 140000 / 739685 Processing Event: 150000 / 739685 Processing Event: 190000 / 739685 Processing Event: 240000 / 739685 Processing Event: 270000 / 739685 Processing Event: 360000 / 739685 Processing Event: 390000 / 739685 Processing Event: 400000 / 739685 Processing Event: 410000 / 739685 Processing Event: 440000 / 739685 Processing Event: 450000 / 739685 Processing Event: 490000 / 739685 Processing Event: 550000 / 739685 Processing Event: 560000 / 739685 Processing Event: 590000 / 739685 Processing Event: 600000 / 739685 Processing Event: 620000 / 739685 Processing Event: 630000 / 739685 Processing Event: 640000 / 739685 Processing Event: 650000 / 739685 Input events=739685, accepted reference tracks=232180 Tracking finished. Output saved to tracking_demo.root
束流径迹
TH2 *hxz=(TH2*)f->Get("htf8xz");
hxz->Draw("colz");
c1->Draw();
TH2 *hyz=(TH2*)f->Get("htf8yz");
hyz->Draw("colz");
c1->Draw();
束流在靶上投影
- 假设靶与束流线垂直。
tree->Draw("ty:tx>>htx(120,-60,60,120,-60,60)","must2Trig","colz");
c1->Draw();
tree->Draw("ty:tx>>htx_beam(120,-60,60,120,-60,60)","beamTrig","colz");
c1->Draw();
$\chi^2/Ndf : tx$
Residual 与 χ²/NDF
先看参与拟合的 1A residual：主体较窄，同时带有较宽的尾部。气体探测器中需要考虑以下来源：
- 粒子在电极膜和气体中散射，使实际轨迹偏离理想直线；较少的大角散射可以形成非 Gaussian 尾部。
- δ 电子在主粒子轨迹以外产生电离，也可能使 delay-line 提前读出信号，从而影响位置重建。
- Pileup、误配及电子学噪声也会增加远离主峰的事例。
用 double-Gaussian 分别描述窄核和宽尾，可以比较主要事例的 residual 宽度及尾部变化。两项是经验描述，不各自对应唯一物理过程；提取单层位置分辨仍采用前面的 excluded residual 方法，并扣除参考径迹的预测方差。
χ²/NDF 检查直线与测量点的一致性。下面比较不同 χ²/NDF 条件下的 residual、靶点位置和方向分布，观察限制尾部时同时损失了哪些事例。本例的位置误差统一取 1 mm，图中的 10、20 用于演示选择效果；实际阈值结合测得的分辨与参考样本确定。两点拟合时 NDF=0，不使用这个比值。
关于散射尾部和 delay-line 的 δ-ray 响应，参见 PDG, Passage of Particles Through Matter, §34.3；H. Kumagai et al., Development of Parallel Plate Avalanche Counter PPAC for BigRIPS fragment separator。
// 先给窄核初值，再拟合包含尾部的分布。
TF1 *g1 = new TF1("g1", "gaus");
TF1 *g2 = new TF1("g2", "gaus");
TF1 *total = new TF1("total", "gaus(0) + gaus(3)");
tree->Draw("dx[0]>>hdx(200,-5,5)");
TH1 *hdx = (TH1*)gROOT->FindObject("hdx");
g1->SetParameters(hdx->GetMaximum(),0,0.3);
hdx->Fit(g1,"Q","",-0.5,0.5);
double sigma = fabs(g1->GetParameter(2));
total->SetParameters(g1->GetParameter(0),g1->GetParameter(1),sigma,hdx->GetMaximum()*0.1,0,3*sigma);
total->SetParLimits(0,0,2*hdx->GetMaximum());
total->SetParLimits(3,0,2*hdx->GetMaximum());
total->SetParLimits(2,0.05,5);
total->SetParLimits(5,0.1,15);
TFitResultPtr residualFit=hdx->Fit(total,"S");
cout << "Core sigma=" << total->GetParameter(2) << ", tail sigma=" << total->GetParameter(5) << " mm" << endl;
gPad->SetLogy(0);
c1->Draw();
**************************************** Minimizer is Minuit2 / Migrad Chi2 = 4053.25 NDf = 194 Edm = 1.69637e-06 NCalls = 205 p0 = 15535.9 +/- 54.2825 (limited) p1 = -0.0690366 +/- 0.000593738 p2 = 0.216938 +/- 0.000651831 (limited) p3 = 852.798 +/- 8.09227 (limited) p4 = -0.0870692 +/- 0.00572479 p5 = 1.36868 +/- 0.00718109 (limited) Core sigma=0.216938, tail sigma=1.36868 mm
tree->Draw("c2ny:dy[0]>>hh(40,-10,10,200,0,1000)","","colz");
c1->SetLogy(0);
c1->Draw();
tree->Draw("c2ny>>hh(200,0,1000)","","");
gPad->SetLogy();
c1->Draw();
gPad->SetLogy(0);
tree->Draw("ty:tx>>hbeam(120,-60,60,120,-60,60)","c2nx<10 && c2ny<10 && beamTrig ","colz");
c1->Draw();//
tree->Draw("ty:tx>>(120,-60,60,120,-60,60)","(c2nx>20 || c2ny>20) && beamTrig ","colz");
c1->Draw();//
待测层 2B 的 excluded residual
2B 没有参与 1A、2A、3 的拟合，因此可直接比较 xx2b[0]（测量值）和 xx2b[1]（预测值）。以下拟合中心 ±1.5 mm 的主峰，打印的是 residual 的 core σ。先检查中心是否偏离零以及尾部；分辨的进一步提取再使用上面的预测方差公式。
TCanvas *cDUT=new TCanvas("cDUT","Excluded residual of PPAC 2B",1000,380);
cDUT->Divide(2,1);
cDUT->cd(1);
tree->Draw("xx2b[0]-xx2b[1]>>hDutX(240,-6,6)","xx2b[0]>-900","hist");
TH1 *hDutX=(TH1*)gROOT->FindObject("hDutX");
TF1 *dutX=new TF1("dutX","gaus",-1.5,1.5);
hDutX->Fit(dutX,"RQ"); dutX->Draw("same");
cDUT->cd(2);
tree->Draw("yy2b[0]-yy2b[1]>>hDutY(240,-6,6)","yy2b[0]>-900","hist");
TH1 *hDutY=(TH1*)gROOT->FindObject("hDutY");
TF1 *dutY=new TF1("dutY","gaus",-1.5,1.5);
hDutY->Fit(dutY,"RQ"); dutY->Draw("same");
cout << "Excluded core sigma X=" << abs(dutX->GetParameter(2))
<< ", Y=" << abs(dutY->GetParameter(2)) << " mm (not intrinsic resolution)\n";
cDUT->Draw();
Excluded core sigma X=0.597244, Y=0.605476 mm (not intrinsic resolution)
PPAC2B x,y,x-y的探测效率
在 ROOT 中，利用 TCut（条件切割）和 TTree::GetEntries(selection) 来快速统计满足特定物理条件的事件数，从而计算效率。
0. 定义物理筛选条件 (TCut)
首先，将复杂的逻辑判断定义为直观的 TCut 变量：
TCut c2btrack = "abs(xx2b[1])<100 && abs(yy2b[1])<60"; // reference position
TCut c2ba = "anode2b>-900";
TCut c2bx = "xx2b[0]>-900"; // 文件中的有效读数，与边界 cut 分开
TCut c2by = "yy2b[0]>-900";
TCut c2bInside = "abs(xx2b[0])<120 && abs(yy2b[0])<75";
1. 选择 reference tracks
分母只由 reference track 定义。内部区域避开物理边缘，不使用 2B 的测量位置。以下统计合并触发样本；比较 beamTrig 与 must2Trig 时，分子、分母应同时加上相同触发条件。
// 绘制外推位置的 2D 分布 (预期击中位置)
tree->Draw("yy2b[1]:xx2b[1]>>h_expected(200,-100,100,200,-100,100)", c2btrack, "colz");
c1->Draw();
// 统计预期穿过灵敏区域的总数
Long64_t N_track = tree->GetEntries(c2btrack);
cout << "预期穿过 PPAC2B 灵敏区的粒子总数 N_track = " << N_track << endl;
预期穿过 PPAC2B 灵敏区的粒子总数 N_track = 232178
2. 计算并对比两类探测效率
实际气态探测器的效率与入射粒子的种类 $(A, Z)$、能量以及击中位置均有关，示例代码求的是全局平均探测效率。
// 下面的分母为同一触发和 fiducial 区域内的 reference tracks。
Long64_t Ntrack = tree->GetEntries(c2btrack);
Long64_t Na = tree->GetEntries(c2btrack && c2ba);
Long64_t Nx = tree->GetEntries(c2btrack && c2bx);
Long64_t Ny = tree->GetEntries(c2btrack && c2by);
Long64_t Nxy = tree->GetEntries(c2btrack && c2bx && c2by);
if (Ntrack==0 || Na==0) throw std::runtime_error("empty efficiency denominator");
cout << "Reference tracks=" << Ntrack << ", anode=" << Na << endl;
for (Long64_t n : {Nx,Ny,Nxy}) {
double eff = double(n)/Ntrack;
cout << "eff=" << eff << " +/- " << sqrt(eff*(1-eff)/Ntrack) << endl;
}
Long64_t Nxa=tree->GetEntries(c2btrack && c2ba && c2bx);
Long64_t Nya=tree->GetEntries(c2btrack && c2ba && c2by);
Long64_t Nxya=tree->GetEntries(c2btrack && c2ba && c2bx && c2by);
cout << "Given anode: x=" << double(Nxa)/Na << ", y=" << double(Nya)/Na
<< ", xy=" << double(Nxya)/Na << endl;
Long64_t Npass=tree->GetEntries(c2btrack && c2bx && c2by && c2bInside);
if (Nxy==0) throw std::runtime_error("no valid xy readings");
cout << "Nref=" << Ntrack << ", Nvalid(xy)=" << Nxy << ", Npass=" << Npass << '\n';
cout << "xy read efficiency=" << double(Nxy)/Ntrack
<< ", boundary-cut survival=" << double(Npass)/Nxy
<< ", selected efficiency=" << double(Npass)/Ntrack << '\n';
Reference tracks=232178, anode=230382 eff=0.931156 +/- 0.000525452 eff=0.929791 +/- 0.000530248 eff=0.868683 +/- 0.00070094 Given anode: x=0.938415, y=0.937039, xy=0.875455 Nref=232178, Nvalid(xy)=201689, Npass=201689 xy read efficiency=0.868683, boundary-cut survival=1, selected efficiency=0.868683
效率的含义
这里得到的是 reference tracking 样本和指定面积内的条件效率。阳极触发条件下的效率用与阳极同时有效的 numerator，不能把未要求 anode 的 Nx 直接除以 Na。给出的误差是大样本 binomial 近似；x-y 联合效率直接计数，不假设 x、y 独立。
示例的 fiducial 区域是 |x|<100 mm、|y|<60 mm，位于探测器有效面内。比较不同触发或 reference 组合时使用相同区域，并同时报告分母。
本次运行中，内部 reference 区域的有效 x-y 读数全部通过物理边界 cut，故 εcut=1。这与刻意避开边缘的选择一致，不能据此推断边缘也没有损失；上游已删除的坐标还需回到 raw data 才能检查。
是否要求阳极符合
不要求 anode 时，可以保留阴极位置有效而阳极漏读的事例，用其他探测器的径迹检查这些位置是否可靠。要求 anode 并应用 time-sum 条件，则能抑制部分噪声、δ-ray 和 pileup 造成的误读，但也会损失阳极未响应或未通过时间选择的事例。
比较两种处理时，以同一批独立 reference tracks 为分母，分别统计位置有效，以及 anode、位置和 time-sum 同时有效的比例；这样才能直接看出新增条件的保留效率。若改用 anode 计数作分母，得到的是另一种条件效率，不能直接当作整体效率。本节文件已做过上游选择，研究 time-sum 的单独影响需要回到原始信号。
作业
Task A: 物理靶点外推与重建信息记录
示例代码演示了固定的 3 层探测器（1A, 2A, 3）拟合。但在实际物理分析中，为了最大化统计量，我们需要利用 F8PPAC1a, 1b, 2a, 2b, 3 这 5 层 的所有有效位置信息，进行动态径迹重建。
修改示例代码：
1. 径迹拟合与状态记录
- 筛选有效点： 对每一个事件，逐一检查 5 层 PPAC 是否给出了有效的位置坐标（排除 -999 或越界噪声）。
- 状态记录： 引入变量
Ndet_x和Ndet_y记录实际参与 $X$ 和 $Y$ 平面拟合的探测器层数（需满足 $N \ge 2$ 才可拟合）；引入数组DetHitX[5], DetHitY[5]记录具体是哪几层参与了拟合（1代表参与，0代表未参与）。
2. 靶位坐标
实验中，物理靶并非垂直放置，而是朝向望远镜倾斜了 $45^\circ$（靶的平面方程x+z=0 ）。
求解靶上坐标： 请将空间直线方程 $x = f_x(z)$ 和 $y = f_y(z)$ 与倾斜的靶平面方程联立，解出束流打在靶面上的真实三维交点坐标
targetX, targetY, targetZ。计算出靶上交点坐标的物理误差 ，以及束流入射方向的投影角误差 。
将上述所有拟合卡方值与 NDF（
NDF>0 时另存 chi2/ndf）、残差（Residual）、交点坐标及物理误差作为新 Branch 存入 TTree。
3. 数据分析
- 靶区接受度： 在
beamTrig（束流触发）条件下，统计“重建的束流打在给定物理靶尺寸范围内”的事件数比例（说明分母采用全部束流触发，还是成功重建的束流触发；两者分别还包含或不包含重建损失）。 - 分别画出不同策略下的事件在靶位上的位置以及入射方向的投影角分布，以及各项的误差分布。
Task B: PPAC 探测效率的多重验证
计算位于最上游的 PPAC1a 和最下游的 PPAC3 的 $x$ 效率、$y$ 效率以及 $(x,y)$ 二维联合效率。
要求进行交叉验证：
- 对比效率类型： 分别计算它们的reference 样本内的响应效率（方法一：Tracking 外推作为分母）和位置读出效率（方法二：自身阳极信号作为分母）。
- 在使用 Tracking 方法时，尝试改变“参考探测器”的组合（例如：测 PPAC3 时，分别尝试用
1a+2a拟合，或用1b+2b拟合），验证在不同基准下求出的效率值是否具有自洽性；差异可能来自几何接受度、粒子组成或参考选择的偏差。
使用两层时记录 NDF=0，不计算 χ²/NDF。对不同 reference 组合比较效率时，检查共同的几何范围和样本条件，不预先要求结果相等。