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

3.4 DSSD 正反面关联（I）：pixel-pixel 相对刻度

正反面能量一致性与相对刻度

同一个粒子在一个 pixel 内沉积能量，两面收集到的电荷来自同一次能量沉积。对电荷收集完整的 single-pixel 事例， $$E_{\rm true}=E_{x,i}=E_{y,j}.$$

电子学实际记录的是 ADC 幅度 $A_{x,i},A_{y,j}$。由于各条增益和 offset 不同，原始道值不必相等。选择 X 面的一条 $X_r$ 定义参考尺度： $$E^{(\rm rel)}\equiv A_{x,r},\qquad k_{x,r}=1,\quad b_{x,r}=0.$$ 其他条通过线性变换归一到这个尺度： $$E_{x,i}^{(\rm rel)}=k_{x,i}A_{x,i}+b_{x,i},\qquad E_{y,j}^{(\rm rel)}=k_{y,j}A_{y,j}+b_{y,j}.$$ 正确配对的 single-pixel 事例应满足 $E_{x,i}^{(\rm rel)}\approx E_{y,j}^{(\rm rel)}$。本节求的是这些相对系数，还不是 MeV 刻度。

从一个 pixel 向另一面传递参考尺度

固定参考条 $X_{i^*}$，选取与 $Y_j$ 对应的事例，拟合 $$A_{x,i^*}=k_{y,j}A_{y,j}+b_{y,j}.$$ 得到 $E_{y,j}^{(\rm rel)}=k_{y,j}A_{y,j}+b_{y,j}$ 后，可以选择已刻度的 $Y_{j^*}$，再拟合 $$E_{y,j^*}^{(\rm rel)}=k_{x,k}A_{x,k}+b_{x,k},$$ 从而得到 $E_{x,k}^{(\rm rel)}=k_{x,k}A_{x,k}+b_{x,k}$。这种 pixel-pixel 方法直观，但每个交叉 pixel 都需要足够统计量。3.5 再把多个 pixel 合并到 strip-wise 拟合中。

与后续多粒子配对的关系

两粒子分别产生 $X_1,X_2$ 和 $Y_1,Y_2$ 时，仅凭条号不能区分直接与交换配对。归一后的 $$\Delta E=|E_x^{(\rm rel)}-E_y^{(\rm rel)}|$$ 可用于检验候选。能量相近时仍可能多解，因此“选最小差”不能保证正确；3.6 会把候选数和配对状态明确保存。

例如，两个粒子分别打在 $(X_1,Y_1)$ 和 $(X_2,Y_2)$，电子学只分别读出两条 X、两条 Y。仅凭这些条号，真正的组合 $\{(X_1,Y_1),(X_2,Y_2)\}$ 与交叉的 ghost 组合 $\{(X_1,Y_2),(X_2,Y_1)\}$ 都可能成立。若两粒子的能量明显不同，正确组合的两面能量应分别一致，错误组合则产生较大差值。

所以相对刻度不仅是让二维条带变窄，也是为后续事例重建提供可比较的物理量。本节先用较简单的 single-pixel 样本求系数；多击样本中的组合选择在 3.6 中展开。

%jsroot on

实验与输入文件

本例使用 25 MeV/u 的 $^{16}$C 束流轰击 $^9$Be 靶，测量前向带电碎片。靶后零度方向的望远镜由三层 32×32 DSSD 和 CsI 组成，条宽约 2 mm、条间隔约 0.1 mm。

data/data_16C.root 已去除 pedestal 事例，硬件触发要求 D1、D2 的 X 面 multiplicity≥2。下列 branch 同时提供按 strip ID 排列的原始数组与紧凑 hit 列表；本节从原始幅度的关联出发。

// xenergy, yenergy, xtime (Array of 32 for each detector)
d1x[32], d1y[32], d1t[32];
d2x[32], d2y[32], d2t[32];
d3x[32], d3y[32], d3t[32];

// hit multiplicity, energy, time, strip ID (Compact arrays based on hits)
d1xhit, d1xe[d1xhit], d1xt[d1xhit], d1xs[d1xhit];
d2xhit, d2xe[d2xhit], d2xt[d2xhit], d2xs[d2xhit];
d3xhit, d3xe[d3xhit], d3xt[d3xhit], d3xs[d3xhit];
d1yhit, d1ye[d1yhit], d1yt[d1yhit], d1ys[d1yhit];
d2yhit, d2ye[d2yhit], d2yt[d2yhit], d2ys[d2yhit];
d3yhit, d3ye[d3yhit], d3yt[d3yhit], d3ys[d3yhit];

TFile *ipf = new TFile("./data/data_16C.root");
if (!ipf || ipf->IsZombie()) throw std::runtime_error("无法打开输入 ROOT 文件");
TTree *tree = (TTree*)ipf->Get("tree");
if (!tree) throw std::runtime_error("输入文件中缺少 tree");
TCanvas *c1 = new TCanvas("c1","c1");

X、Y 两面的有效条数

硬件 multiplicity 条件影响样本组成。软件阈值、读出效率和能量共享会使离线 hit 数与触发 hit 数不同；X 侧触发也不保证 Y 侧具有相同 multiplicity。

tree->Draw("d1xhit: d1yhit>>(15, 0, 15, 15, 0, 15)", "", "colz");
gPad->SetLogz();
c1->Draw();

多条信号可能来自多个粒子，也可能来自一个粒子的电荷共享或串扰。仅凭这张 multiplicity 图还不能确定各类的比例，下面继续检查条号和能量关联。

tree->Scan("d1xe: d1xs: d1ye: d1ys", "", "", 10, 1);

***********************************************************************
*    Row   * Instance *      d1xe *      d1xs *      d1ye *      d1ys *
***********************************************************************
*        1 *        0 *      4800 *        17 *      5775 *        18 *
*        1 *        1 *       982 *        18 *           *           *
*        2 *        0 *      3069 *        23 *      3352 *        11 *
*        2 *        1 *       341 *        24 *       111 *        10 *
*        3 *        0 *      5325 *        23 *      4522 *        15 *
*        3 *        1 *           *           *       922 *        16 *
*        4 *        0 *      5822 *        21 *      5688 *        15 *
*        4 *        1 *           *           *       271 *        14 *
*        5 *        0 *      4048 *        11 *      5213 *        12 *
*        5 *        1 *      1149 *        12 *           *           *
*        6 *        0 *      3719 *        14 *      4255 *        16 *
*        6 *        1 *       374 *        13 *           *           *
*        7 *        0 *      4546 *        12 *      4961 *         9 *
*        7 *        1 *       297 *        11 *           *           *
*        8 *        0 *      4051 *        22 *      3877 *         9 *
*        8 *        1 *           *           *       329 *        10 *
*        9 *        0 *      6720 *        13 *      6513 *        12 *
*        9 *        1 *       142 *        10 *       418 *        13 *
*        9 *        2 *       109 *        14 *           *           *
*       10 *        0 *      3873 *        20 *      3007 *         8 *
*       10 *        1 *           *           *       974 *         9 *
***********************************************************************

正反面幅度关联

先看 X[12] 与 Y[13] 的关联。主条带是建立相对刻度的依据；条带外既可能有 charge sharing，也可能有多个粒子形成的错误组合、阈值效应或其他本底。不能把全部离群点都归为 sharing。

gPad->SetLogy(0);
tree->Draw("d1x[12]: d1y[13]>>(1000, 0, 8000, 1000, 0, 8000)", "", "colz");
c1->Draw();//front-back correlation

1. 先试两面各一个 hit

xhit==1 && yhit==1 可以减少组合歧义。这里硬件偏向多击事件，这种选择的统计量有限，因此还要考察局部条件。

c1->Clear();
TCut chit1 = "d1xs==12 && d1ys==13 && d1xhit==1 && d1yhit==1";
tree->Draw("d1xe: d1ye>>(1000, 0, 8000, 1000, 0, 8000)", chit1, "colz");
c1->Draw();

2. 邻条无有效信号：local isolation

对于 X[12] 检查 X[11]、X[13]，对于 Y[13] 检查 Y[12]、Y[14]。下面用邻条幅度小于 50 的条件作示例；这个软件阈值应结合各通道 pedestal 和噪声确定，不等同于所有通道已知的硬件阈值。

它排除了可见的邻条信号，但阈值以下的 sharing 仍可能存在。两条都 isolated，也不能证明它们一定来自同一个粒子；继续检查二维条带和 residual。

c1->Clear();
TCut cveto = "d1x[11]<50 && d1x[13]<50 && d1y[12]<50 && d1y[14]<50"; // 邻条无超过所选阈值的信号
TCut c1213 = "d1x[12]>200 && d1y[13]>200 && d1y[13]<8000" && cveto;
tree->Draw("d1x[12]: d1y[13]>>h2(1000, 0, 8000, 1000, 0, 8000)", c1213, "colz");
c1->Draw();

1. 取出逐事例关联数据

将所选事例直接交给 TGraph 做 unbinned 直线拟合，避免把 TH2 的 bin 中心当成测量值。TH2 用来显示同一批事例；并不是所有关联拟合都必须采用这种做法。

tree->SetEstimate(tree->GetEntries()+1);
tree->Draw("d1x[12]:d1y[13]", c1213, "goff");
TGraph *gr = new TGraph(tree->GetSelectedRows(), tree->GetV2(), tree->GetV1());
gr->SetMarkerSize(0.2);
gr->Draw("A*");
c1->Draw();

普通 least squares 与 ROB 拟合的比较

这里横轴为 Y[13] 的原始幅度，纵轴为 X[12] 的原始幅度，拟合 X[12]=b+k*Y[13]。普通 least squares 对远离主条带的点较敏感；ROOT 的 ROB 使用 least trimmed squares，以残差较小的子集估计直线。它适用于参数线性模型，不能把这个选项直接用于任意非线性 TF1。

下面比较两条线及各自的 residual，不预先认定某个拟合更好。ROB 也不能代替事件鉴别。参见 ROOT robust fit 示例。

// 1. 普通最小二乘拟合
TF1 *fStd = new TF1("fStd", "pol1", 200, 8000);
fStd->SetLineColor(kRed);
gr->Fit(fStd, "RQ"); 

// 2. ROB：Least Trimmed Squares，降低离群点的影响
TF1 *fRob = new TF1("fRob", "pol1", 200, 8000);
fRob->SetLineColor(kBlue);
fRob->SetLineWidth(1);
// R：使用 TF1 的范围；Q：简短输出；+：保留前一个拟合函数
gr->Fit(fRob, "R+ ROB Q"); 

// 在同一组散点上显示两种拟合
gr->Draw("A*");

// 3. 图例
TLegend *leg = new TLegend(0.15, 0.75, 0.45, 0.88);
leg->SetBorderSize(0);
leg->AddEntry(fStd, "Least squares", "l");
leg->AddEntry(fRob, "ROB", "l");
leg->Draw();

c1->Draw();

3. 残差与幅度的关系

定义 residual 为 d1x[12] − (p0 + p1*d1y[13])。除检查主条带是否在零附近，还要看低、中、高幅度区是否存在不同的偏移或趋势。

本例低幅度事例密集，逐事例拟合时这一段对目标函数的累计贡献较大。即使使用 ROB，拟合也可能更侧重低幅度区，高幅度区的符合程度却未必同样好；若刻度近似或残留事件成分随幅度变化，这种影响会反映在 residual 中。

3.5 逐条归一采用按待刻度条幅度抽样的处理，降低密集区的主导作用，再用完整候选样本检查归一结果。若抽样后仍有系统趋势，还需检查电子学非线性、两轴测量误差、阈值及残留混合事件；不能仅凭事例密度就断定趋势一定是统计效应。

// 提取 ROB 拟合的截距与斜率
double p0 = fRob->GetParameter(0);
double p1 = fRob->GetParameter(1);

// 横轴 residual = X-(p0+p1*Y)，纵轴为 Y 的原始幅度
TString stree;
stree.Form("d1y[13]:d1x[12]-(%f*d1y[13]+%f)>>ha(100,-40,40,1000,0,8000)", p1, p0);

// 检查 residual 是否随幅度变化
c1->Clear();
tree->Draw(stree.Data(), c1213, "colz");
c1->Draw();

作业：保存 front-back 候选

对三层 DSSD，参照上述邻条条件和二维主条带，保存可用于相对刻度的 x-y 候选组合。先在单击样本上观察主条带，再决定多击样本可采用的范围。

输出记录 source_entry、ix、iy、xe、ye。一个输入事例可以产生多个候选组合，所以不要把输出行数当成粒子数。TCutG::IsInside(x,y) 的参数顺序对应作图的横、纵轴。

事件循环的核心如下；输入分支的绑定和输出树建立沿用 1.2 的方法。

// 在输入的每个事例中，枚举满足局部条件的 X、Y 条。
for(int i=0;i<32;++i) {
    if(d1x[i]<200) continue;
    if(i>0 && d1x[i-1]>=50) continue;
    if(i<31 && d1x[i+1]>=50) continue;
    for(int j=0;j<32;++j) {
        if(d1y[j]<200) continue;
        if(j>0 && d1y[j-1]>=50) continue;
        if(j<31 && d1y[j+1]>=50) continue;
        // 本节画的是 X:Y，故横轴 Y、纵轴 X。
        if(!cutXY->IsInside(d1y[j],d1x[i])) continue;
        ix=i; iy=j; xe=d1x[i]; ye=d1y[j];
        source_entry=jentry;
        tout->Fill();
    }
}

另存分析输出，不覆盖提供的原始数据。下面的 data/d1xy.root 是供 3.5 使用的参考候选文件。

DSSD1 参考候选样本

TCanvas *c1 = new TCanvas("c1","c1");
TFile *fin = new TFile("./data/d1xy.root");
if (!fin || fin->IsZombie()) throw std::runtime_error("无法打开输入 ROOT 文件");
TTree *tree =(TTree *)fin->Get("tree");
if (!tree) throw std::runtime_error("输入文件中缺少 tree");
tree->Draw("ye:xe>>(1000,0,8000,1000,0,8000)","","colz");
c1->SetLogz();
c1->Draw();

Warning in <TCanvas::Constructor>: Deleting canvas with same name: c1