3.2 利用可变长数组存储 DSSD Hit 结构¶
描述一个粒子与 DSSD 的相互作用，至少需要如下信息：
- position: x or y side, strip number in x/y side
- energy: raw, calibrated
- time: raw, calibrated
这三类信息对应的是：粒子打在探测器哪一侧、哪一条；沉积了多少能量；以及对应的时间信号。后续的数据结构设计，目标就是把这些信息以适合逐事件分析的方式组织起来。
常用的固定数组表达方式¶
常用的数据表达方式是数组。假设有三个 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
这种写法的优点是结构直接，适合原始数据存储；但对 DSSD 这种稀疏事件来说，它并不适合后续分析。因为在绝大多数情况下，一个 DSSD 只有少数几条包含有用信息，数组中大部分元素都是无用的。固定数组的记录方式不仅造成内存和存储空间浪费，也不利于表达一个 DSSD 事件的整体特性，例如不方便一次检查一个 DSSD 中多条的能量信息。
如果有多个同类型探测器¶
如果系统中有多个同类型探测器，更适合把事件整理成 hit structure。对于 x side，可以写成：
Int_t xhit;
Int_t xdet[xhit]; // dssd: 1-3
Int_t xid[xhit]; // 1-128
Double_t xe[xhit], xec[xhit];
对于 y side，可以写成：
Int_t yhit;
Int_t ydet[yhit]; // dssd: 1-3
Int_t yid[yhit]; // 1-128
Double_t ye[yhit], yec[yhit];
这时事件不再被表示成“每个探测器、每一条都留一个位置”，而是被表示成“这个事件一共有哪些有效 hit，每个 hit 属于哪一个 DSSD、哪一条、对应什么能量”。这样，数据结构和物理事件就是一致的。
数据结构¶
-TTree 的Branch 格式为
Int_t pe[48]; // Pie side raw energy
Int_t re[48]; // Ring side raw energy
这种表示方式适合原始存储，但不适合逐事件分析。对一个具体事件来说，我们真正关心的是：Pie 一侧有几条有信号、Ring 一侧有几条有信号、它们的条号是什么、能量是多少。也就是说，我们关心的是 hit，而不是整个固定数组。
In [3]:
//%jsroot on
TFile *f = new TFile("s4.root");
TTree *tree = (TTree*)f->Get("tree");
TCanvas *c1 = new TCanvas;
In [4]:
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 * ***********************************************
Type <CR> to continue or q to quit ==>
观察事件中x，y有效信息¶
In [6]:
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 * ***********************************************
Type <CR> to continue or q to quit ==>
- 将某一方向的信息记为 (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 pec[48]; // Pie calibrated energy
Double_t rec[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;
这个结构直接回答了事件分析中的关键问题：有几条、是哪几条、各自能量是多少。它比固定数组更适合后续所有相关分析。
ROOT - Adding a Branch with a Variable Length Array¶
在 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 结构数据¶
In [11]:
//%jsroot on
TFile *fhit=new TFile("s4hit.root");
TTree *tree=(TTree*) fhit->Get("tree");
In [12]:
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 * *............................................................................*
In [13]:
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给出对应的能量。
与原始固定数组相比，现在可以直接判断一个事件是单击、相邻双击，还是更复杂的多击事件
In [16]:
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();
In [18]:
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
相邻条能量关联¶
- 出现清楚的k=-1关联，主要对应相邻条能量共享(见3.3节)。
In [20]:
c2->Clear();
c2->Divide(2,1);
c2->cd(1);
tree->Draw("rer[0]:rer[1]>>h1(400,0,1600,400,0,1600)","abs(rid[0]-rid[1])==1","colz");
c2->cd(2);
tree->Draw("rer[0]:rer[1]>>h2(400,0,1600,400,0,1600)","abs(rid[0]-rid[1])!=1","colz");
c2->Draw();
ring的相邻条能量之和与pie的能量关联¶
- 出现k=1的关联(见3.4节)。
In [22]:
c2->Clear();
c2->Divide(1,1);
tree->Draw("per[0]:rer[0]+rer[1]>>(400,0,1600,400,0,1600)","rhit ==2 && abs(rid[0]-rid[1])==1","colz");
c2->Draw();
选择特定条之间关联¶
In [24]:
c2->Clear();
c2->Divide(2,1);
c2->cd(1);
tree->Draw("rer[0]:rer[1]>>h1(400,0,1600,400,0,1600)","rid[0] == 13 && rid[1] == 14 ","colz");
c2->cd(2);
tree->Draw("per[0]:rer[0]+rer[1]>>(400,0,1600,400,0,1600)","pid[0] == 13 && rid[0] == 13 && rid[1] == 14 ","colz");
c2->Draw();
In [26]:
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();
能量已刻度¶
In [28]:
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();
作业¶
- 对
pe参数连续化，利用一般策略（选出拟合 $\chi^2/\text{NDF}$ 最小的一组作为实际的峰位-能量组合）完成 Pie 侧所有 48 条的能量刻度。分析每个条的 $\sigma$ 和 $\chi^2/\text{NDF}$，根据结果调整拟合策略。 - 将带有高斯拟合曲线的所有单条原始能谱、刻度后的单条能谱保存到一个 ROOT 文件中。将刻度后所有条的能谱累积在一个 TH1 上，给出 4 个参与刻度的 $\alpha$ 峰的分辨。
- 提取各条的分辨率，画出 FWHM-条id 的关系图并放在一个
TGraph上（参照下图）。同时，对比刻度前后 Pie 侧所有条的能量分布。 - 将刻度前和刻度后的数据转换成 hit 结构，对比并展示刻度前后 Pie 侧所有条的整体能量分布情况。其中，Ring 侧的条事件需按照能量由大到小的顺序进行排序后，再填入 hit 结构。
In [ ]: