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

3.2 利用可变长数组存储 DSSD Hit 结构

描述一个粒子与 DSSD 的相互作用，至少需要如下信息：

• 位置：探测器的 X/Y 面及对应条号
• 能量：原始幅度与刻度后能量
• 时间：原始读数与刻度后时间

这三类信息对应的是：粒子打在探测器哪一侧、哪一条；沉积了多少能量；以及对应的时间信号。后续的数据结构设计，目标就是把这些信息以适合逐事件分析的方式组织起来。

常用的固定数组表达方式

常用的数据表达方式是数组。假设有三个 DSSD（128×128），则 energy 可以写成：

Int_t xe[3][128], ye[3][128];            // raw
Double_t xec[3][128], yec[3][128];      // calibrated

time 可以写成：

Int_t xt[3][128], yt[3][128];            // raw
Double_t xtc[3][128], ytc[3][128];      // calibrated

这种写法的优点是结构直接，适合原始数据存储；后续若主要处理有效条，也可以转换为紧凑的 hit 列表。在本例中，一个事件通常只有少数条有信号。固定数组中重复的无效标记也可被 ROOT 压缩；hit 结构的主要优点是便于表达和处理有效条，而不是保证文件总会更小。固定数组的下标直接是 strip ID；hit 列表的下标则是本事件内的 hit 编号。两种结构都能使用，区别在于访问方式。

如果有多个同类型探测器

如果系统中有多个同类型探测器，更适合把事件整理成 hit structure。对于 x side，可以写成：

Int_t xhit;
Int_t xdet[384];      // dssd: 1-3
Int_t xid[384];       // 0-127
Double_t xe[384], xec[384];

对于 y side，可以写成：

Int_t yhit;
Int_t ydet[384];      // dssd: 1-3
Int_t yid[384];       // 0-127
Double_t ye[384], yec[384];

这时事件不再被表示成“每个探测器、每一条都留一个位置”，而是被表示成“这个事件一共有哪些有效 hit，每个 hit 属于哪一个 DSSD、哪一条、对应什么能量”。这些仍是 strip hits；一个粒子可能产生多个 hit，后面还需重建粒子。

上述数组容量 384 对应三个 128 条探测器；事件中的有效长度由 xhit、yhit 给出，写进 ROOT leaflist。下文改用当前 S4 数据的容量 48，以及 phit/rhit 演示。

%jsroot on

数据结构

-TTree 的Branch 格式为

Int_t pe[48];   // Pie side raw energy
Int_t re[48];   // Ring side raw energy

这种表示方式便于按通道保存原始数据；逐事件分析时，则经常需要先找出有信号的条。对一个具体事件来说，我们真正关心的是：Pie 一侧有几条有信号、Ring 一侧有几条有信号、它们的条号是什么、能量是多少。也就是说，我们关心的是 hit，而不是整个固定数组。

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

tree->SetScanField(0); // notebook 中关闭翻页，仍只查看指定的 1 个事件
tree->Scan("pe:re", "", "", 1, 1);//第4参数-事件数，第5参数-起始事件编号
//jupyter上的Scan必须指定第4个参数(事件数)

***********************************************
*    Row   * Instance *        pe *        re *
***********************************************
*        1 *        0 *       -10 *       -10 *
*        1 *        1 *       -10 *       -10 *
*        1 *        2 *       -10 *       -10 *
*        1 *        3 *       -10 *       -10 *
*        1 *        4 *       -10 *       -10 *
*        1 *        5 *       -10 *       -10 *
*        1 *        6 *       -10 *       -10 *
*        1 *        7 *       -10 *       -10 *
*        1 *        8 *       -10 *       -10 *
*        1 *        9 *       -10 *       -10 *
*        1 *       10 *       983 *       -10 *
*        1 *       11 *       -10 *       -10 *
*        1 *       12 *       -10 *       -10 *
*        1 *       13 *       -10 *       -10 *
*        1 *       14 *       -10 *       -10 *
*        1 *       15 *       -10 *       -10 *
*        1 *       16 *       -10 *       -10 *
*        1 *       17 *       -10 *       -10 *
*        1 *       18 *       -10 *       -10 *
*        1 *       19 *       -10 *       -10 *
*        1 *       20 *       -10 *       -10 *
*        1 *       21 *       -10 *       -10 *
*        1 *       22 *       -10 *       -10 *
*        1 *       23 *       -10 *       -10 *
*        1 *       24 *       -10 *       -10 *
*        1 *       25 *       -10 *       -10 *
*        1 *       26 *       -10 *       -10 *
*        1 *       27 *       -10 *       -10 *
*        1 *       28 *       -10 *       -10 *
*        1 *       29 *       -10 *       -10 *
*        1 *       30 *       -10 *       -10 *
*        1 *       31 *       -10 *       -10 *
*        1 *       32 *       -10 *       -10 *
*        1 *       33 *       -10 *       -10 *
*        1 *       34 *       -10 *       -10 *
*        1 *       35 *       -10 *       -10 *
*        1 *       36 *       -10 *       -10 *
*        1 *       37 *       -10 *       -10 *
*        1 *       38 *       -10 *       -10 *
*        1 *       39 *       -10 *       -10 *
*        1 *       40 *       -10 *       -10 *
*        1 *       41 *       -10 *       -10 *
*        1 *       42 *       -10 *       -10 *
*        1 *       43 *       -10 *       -10 *
*        1 *       44 *       -10 *       -10 *
*        1 *       45 *       -10 *       -10 *
*        1 *       46 *       -10 *       -10 *
*        1 *       47 *       -10 *       -10 *
***********************************************

查看一个事例中两面的有效信号

tree->Scan("pe:re", "", "", 1, 2);

***********************************************
*    Row   * Instance *        pe *        re *
***********************************************
*        2 *        0 *       -10 *       -10 *
*        2 *        1 *       -10 *       -10 *
*        2 *        2 *       -10 *       -10 *
*        2 *        3 *       -10 *       -10 *
*        2 *        4 *       -10 *       -10 *
*        2 *        5 *       -10 *       -10 *
*        2 *        6 *       -10 *       -10 *
*        2 *        7 *       -10 *       -10 *
*        2 *        8 *       978 *       -10 *
*        2 *        9 *       -10 *       -10 *
*        2 *       10 *       -10 *       -10 *
*        2 *       11 *       -10 *       -10 *
*        2 *       12 *       -10 *       -10 *
*        2 *       13 *       -10 *       -10 *
*        2 *       14 *       -10 *       -10 *
*        2 *       15 *       -10 *       -10 *
*        2 *       16 *       -10 *       -10 *
*        2 *       17 *       -10 *       -10 *
*        2 *       18 *       -10 *       -10 *
*        2 *       19 *       -10 *       -10 *
*        2 *       20 *       -10 *       -10 *
*        2 *       21 *       -10 *       -10 *
*        2 *       22 *       -10 *       -10 *
*        2 *       23 *       -10 *       -10 *
*        2 *       24 *       -10 *       -10 *
*        2 *       25 *       -10 *       -10 *
*        2 *       26 *       -10 *       -10 *
*        2 *       27 *       -10 *       -10 *
*        2 *       28 *       -10 *       -10 *
*        2 *       29 *       -10 *       -10 *
*        2 *       30 *       -10 *       -10 *
*        2 *       31 *       -10 *       -10 *
*        2 *       32 *       -10 *       -10 *
*        2 *       33 *       -10 *       -10 *
*        2 *       34 *       -10 *       -10 *
*        2 *       35 *       -10 *       -10 *
*        2 *       36 *       -10 *       -10 *
*        2 *       37 *       -10 *       -10 *
*        2 *       38 *       -10 *       -10 *
*        2 *       39 *       -10 *       -10 *
*        2 *       40 *       -10 *       -10 *
*        2 *       41 *       -10 *       -10 *
*        2 *       42 *       -10 *       623 *
*        2 *       43 *       -10 *       287 *
*        2 *       44 *       -10 *       -10 *
*        2 *       45 *       -10 *       -10 *
*        2 *       46 *       -10 *       -10 *
*        2 *       47 *       -10 *       -10 *
***********************************************

• 将某一方向的信息记为 (id,energy),则上一个事件只有如下有效信号
• pie: (8,978)
• ring: (42,623), (43,287)

在固定数组结构中，"pe>0 && re>0" 只能用于判断下标相同的 pe[i] 与 re[i] 是否同时大于零；它并不能直接表达一个事件中 Pie 和 Ring 两侧所有有效 hit 之间的关联关系。

新的事件表达方式：Hit Structure

对每个事件，只保存真正有信号的条。

对当前 S4 数据，可以定义：

Int_t phit;      // number of valid Pie hits
Int_t pid[48];   // Pie strip id
Int_t per[48];   // Pie raw energy

Int_t rhit;      // number of valid Ring hits
Int_t rid[48];   // Ring strip id
Int_t rer[48];   // Ring raw energy

如果加入刻度后的能量，也可以进一步定义：

Double_t pee[48];   // Pie calibrated energy
Double_t ree[48];   // Ring calibrated energy

这样，上面的事件就可以写成：

phit = 1;
pid[0] = 8;
per[0] = 978;

rhit = 2;
rid[0] = 42; rer[0] = 623;
rid[1] = 43; rer[1] = 287;

这个结构直接回答了事件分析中的关键问题：有几条、是哪几条、各自能量是多少。注意 rid[0]=42 表示“第 0 个 hit 来自第 42 条”，不是“第 0 条”。

在 TTree 中建立可变长数组分支

在 ROOT 里，可变长数组分支的基本写法是：先定义一个“长度变量”，再定义一个足够大的本地数组，随后在 Branch 的 leaflist 里用这个长度变量控制每个事件的实际数组长度。例如：

Int_t phit;
Int_t pid[48];
Int_t per[48];

tree->Branch("phit", &phit, "phit/I");
tree->Branch("pid",  pid,   "pid[phit]/I");
tree->Branch("per",  per,   "per[phit]/I");

Ring 一侧同理：

Int_t rhit;
Int_t rid[48];
Int_t rer[48];

tree->Branch("rhit", &rhit, "rhit/I");
tree->Branch("rid",  rid,   "rid[rhit]/I");
tree->Branch("rer",  rer,   "rer[rhit]/I");

这里要注意两点。

第一，内存中的数组仍然要预先给出一个足够大的最大长度，例如 48；
第二，真正写入每个事件的数据长度由 phit 或 rhit 决定。也就是说，ROOT 允许“每个事件实际存多少个元素”随事件而变化，而 branch 的声明方式就是 array_name[length_var]/type。

生成hit结构数据

• 逐事件扫描 48 条，只把真正有信号的条重新组织成 hit 列表。这样，新的树就不再只是“通道数组”，而是“事件中的有效 hit 集合”。

Long64_t nentries = tin->GetEntries();
for(Long64_t jentry = 0; jentry < nentries; jentry++) {
    tin->GetEntry(jentry);

    phit = 0;
    rhit = 0;
    for(int i = 0; i < 48; i++) {
        if(pe[i] > 0) {
            pid[phit] = i;
            per[phit] = pe[i];
            phit++;
        }
        if(re[i] > 0) {
            rid[rhit] = i;
            rer[rhit] = re[i];
            rhit++;
        }
    }
    tout->Fill();
}

Hit 结构数据

TFile *fhit=new TFile("s4hit.root");
if (!fhit || fhit->IsZombie()) throw std::runtime_error("无法打开输入 ROOT 文件");
TTree *tree=(TTree*) fhit->Get("tree");
if (!tree) throw std::runtime_error("输入文件中缺少 tree");

tree->Print();//ree，pee为刻度后数据。

******************************************************************************
*Tree    :tree      : alpha                                                  *
*Entries :  5790825 : Total =      2541888463 bytes  File  Size =  188942263 *
*        :          : Tree compression factor =  13.46                       *
******************************************************************************
*Br    0 :pe        : pe[48]/I                                               *
*Entries :  5790825 : Total  Size= 1112292512 bytes  File Size  =   29248218 *
*Baskets :     5099 : Basket Size=   25600000 bytes  Compression=  38.03     *
*............................................................................*
*Br    1 :re        : re[48]/I                                               *
*Entries :  5790825 : Total  Size= 1112292512 bytes  File Size  =   41772927 *
*Baskets :     5099 : Basket Size=   25600000 bytes  Compression=  26.62     *
*............................................................................*
*Br    2 :rhit      : rhit/I                                                 *
*Entries :  5790825 : Total  Size=   23179708 bytes  File Size  =    2833849 *
*Baskets :      177 : Basket Size=     291840 bytes  Compression=   8.18     *
*............................................................................*
*Br    3 :rid       : rid[rhit]/I                                            *
*Entries :  5790825 : Total  Size=   52791843 bytes  File Size  =   14089701 *
*Baskets :      357 : Basket Size=    7382528 bytes  Compression=   3.75     *
*............................................................................*
*Br    4 :ree       : ree[rhit]/F                                            *
*Entries :  5790825 : Total  Size=   52791850 bytes  File Size  =   34140607 *
*Baskets :      357 : Basket Size=    7382528 bytes  Compression=   1.55     *
*............................................................................*
*Br    5 :rer       : rer[rhit]/I                                            *
*Entries :  5790825 : Total  Size=   52791843 bytes  File Size  =   22519855 *
*Baskets :      357 : Basket Size=    7382528 bytes  Compression=   2.34     *
*............................................................................*
*Br    6 :phit      : phit/I                                                 *
*Entries :  5790825 : Total  Size=   23179708 bytes  File Size  =    2488050 *
*Baskets :      177 : Basket Size=     291840 bytes  Compression=   9.31     *
*............................................................................*
*Br    7 :pid       : pid[phit]/I                                            *
*Entries :  5790825 : Total  Size=   37522837 bytes  File Size  =   10440561 *
*Baskets :      288 : Basket Size=    7190016 bytes  Compression=   3.59     *
*............................................................................*
*Br    8 :pee       : pee[phit]/F                                            *
*Entries :  5790825 : Total  Size=   37522844 bytes  File Size  =   18283102 *
*Baskets :      288 : Basket Size=    7190016 bytes  Compression=   2.05     *
*............................................................................*
*Br    9 :per       : per[phit]/I                                            *
*Entries :  5790825 : Total  Size=   37522837 bytes  File Size  =   13024115 *
*Baskets :      288 : Basket Size=    7190016 bytes  Compression=   2.88     *
*............................................................................*

tree->Scan("phit:pid:per:rhit:rid:rer","","",10,1);

***********************************************************************************************
*    Row   * Instance *      phit *       pid *       per *      rhit *       rid *       rer *
***********************************************************************************************
*        1 *        0 *         1 *        10 *       983 *         0 *           *           *
*        2 *        0 *         1 *         8 *       978 *         2 *        42 *       623 *
*        2 *        1 *         1 *           *           *         2 *        43 *       287 *
*        3 *        0 *         1 *        27 *       928 *         2 *         4 *       740 *
*        3 *        1 *         1 *           *           *         2 *         5 *       169 *
*        4 *        0 *         0 *           *           *         2 *         1 *       670 *
*        4 *        1 *         0 *           *           *         2 *         2 *       379 *
*        5 *        0 *         1 *        44 *      1402 *         0 *           *           *
*        6 *        0 *         0 *           *           *         2 *        27 *        58 *
*        6 *        1 *         0 *           *           *         2 *        28 *       854 *
*        7 *        0 *         0 *           *           *         2 *        20 *       240 *
*        7 *        1 *         0 *           *           *         2 *        21 *       807 *
*        8 *        0 *         0 *           *           *         2 *        45 *      1236 *
*        8 *        1 *         0 *           *           *         2 *        46 *       183 *
*        9 *        0 *         0 *           *           *         1 *        11 *        52 *
*       10 *        0 *         0 *           *           *         2 *        35 *       811 *
*       10 *        1 *         0 *           *           *         2 *        36 *       624 *
***********************************************************************************************

每个事件的结构：

• phit、rhit 给出 Pie 和 Ring 两侧的 multiplicity；
• pid、rid 给出有效条号；
• per、rer 给出对应的能量。

与原始固定数组相比，现在可以直接判断一个事件是单击、相邻双击，还是更复杂的多击事件

排序不能拆开 hit

若按能量排序，同一个 hit 的 strip ID、raw energy、calibrated energy 和 time 要一起移动。不能只排序能量数组，再用原来的条号。判断相邻条用 abs(rid[i]-rid[j])==1，不是 j==i+1。

最小的交换例子如下，std::swap 逐项保持对应关系。多个 hit 可先对索引排序，或采用 3.8 中的 hit 对象。

if (rhit==2 && ree[0]<ree[1]) {
    std::swap(rid[0],rid[1]);
    std::swap(rer[0],rer[1]);
    std::swap(ree[0],ree[1]);
    // 若记录了时间，也交换两个 hit 的时间。
}

利用 hit 结构分析事例

有效条数（multiplicity）

这两个分布给出 Pie 和 Ring 两侧的 multiplicity，是后续事件分类的起点。

• pie方向：phit=1的事件为主
• ring方向：rhit=2的事件为主

TCanvas *c2 = new TCanvas("c2","c2",800,400);
c2->Clear();
c2->Divide(2,1);
c2->cd(1);
tree->Draw("rhit");
c2->cd(2);
tree->Draw("phit");
gStyle->SetOptStat(0);
c2->Draw();

Ring 面的双击事例

对于 rhit==2 的事件，可以直接画出两个 Ring hit 的条号关系：

• 说明主要为相邻条事件

TCanvas *c1 =new TCanvas("c1","c1",600,600);
c1->Clear();
c1->Divide(1,1);
tree->Draw("abs(rid[0]-rid[1])>>(10,0,10)", "rhit==2", "colz");
c1->Draw();

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

相邻条能量关联

• 图中有清楚的负斜率条带，对应同一粒子的能量在相邻条间分配；两条响应接近时斜率约为 −1，增益不同时则不必如此，见 3.3 节。

c2->Clear();
c2->Divide(2,1);
c2->cd(1);
tree->Draw("rer[0]:rer[1]>>h1(400,0,1600,400,0,1600)","rhit==2 && abs(rid[0]-rid[1])==1","colz");
c2->cd(2);
tree->Draw("rer[0]:rer[1]>>h2(400,0,1600,400,0,1600)","rhit==2 && abs(rid[0]-rid[1])!=1","colz");
c2->Draw();

ring的相邻条能量之和与pie的能量关联

• 两面的幅度呈正相关。统一各条的响应后，对完整收集的单粒子事例，两面能量应一致，见 3.4 节。

c2->Clear();
c2->Divide(1,1);
tree->Draw("per[0]:rer[0]+rer[1]>>(400,0,1600,400,0,1600)","phit==1 && rhit==2 && abs(rid[0]-rid[1])==1","colz");
c2->Draw();

选择特定条之间关联

c2->Clear();
c2->Divide(2,1);
c2->cd(1);
tree->Draw("rer[0]:rer[1]>>h1(400,0,1600,400,0,1600)","rhit==2 && rid[0] == 13 && rid[1] == 14 ","colz");
c2->cd(2);
tree->Draw("per[0]:rer[0]+rer[1]>>(400,0,1600,400,0,1600)","phit==1 && rhit==2 && pid[0] == 13 && rid[0] == 13 && rid[1] == 14 ","colz");
c2->Draw();

观察pie和ring的所有条的能量分布

利用 hit 结构可以非常方便地检查所有有效条的能量分布。通过这种方式，我们可以直接分析阈值特性、单条响应以及整体数据质量，而无需再回到原有的固定数组表示。

能量未刻度

c2->Clear();
c2->Divide(2,1);
c2->cd(1);
tree->Draw("per:pid>>hh1(48,0,48,800,0,1600)","","colz");
c2->cd(2);
tree->Draw("rer:rid>>hh2(48,0,48,800,0,1600)","","colz");
c2->Draw();

能量已刻度

c2->Clear();
c2->Divide(2,1);
c2->cd(1);
tree->Draw("pee:pid>>ha1(48,0,48,800,0,10000)","","colz");
c2->cd(2);
tree->Draw("ree:rid>>ha2(48,0,48,800,0,10000)","","colz");
c2->Draw();

作业

1. 参照 3.1 节完成 Pie 侧所有 48 条的能量刻度。先对有效的整数道值进行 dithering，再应用各条刻度。按已知谱线的能量顺序枚举候选峰组合，比较拟合残差与峰身份，确定对应关系；再比较各条的峰宽和拟合质量。
2. 将带有高斯拟合曲线的所有单条原始能谱、刻度后的单条能谱保存到一个 ROOT 文件中。将刻度后所有条的能谱累积在一个 TH1 上，给出 4 个参与刻度的 $\alpha$ 峰的分辨。
3. 提取各条的分辨率，画出 FWHM-条id 的关系图并放在一个 TGraph 上（参照下图）。同时，对比刻度前后 Pie 侧所有条的能量分布。
4. 将刻度前和刻度后的数据转换成 hit 结构，对比并展示刻度前后 Pie 侧所有条的整体能量分布情况。其中，Ring 侧的条事件需按照能量由大到小的顺序进行排序后，再填入 hit 结构。

固定容量数组在每次 Fill() 前应满足 0 ≤ phit,rhit ≤ 48。示例文件的 pee/ree 单位为 keV，与 3.1 中的 MeV 相差 1000 倍。