TTree 中使用 vector 类型的 Branch¶
在探测器数据分析中,一个事件中往往会出现多个 hit。若直接使用固定长度数组,或依赖 hit 数目的动态数组保存这些信息,虽然能够完成数据存储,但在后续分析中往往不便于按单个 hit 进行排序、筛选和配对。利用 vector,可以把同一个 hit 的条带编号、能量和时间统一组织起来,再将一个事件中的全部 hit 保存在同一个容器中,从而使数据结构更加清晰,也更便于后续分析处理。
以 data_16C.root 为例,原始 ROOT 文件中既包含固定长度数组,也包含由 hit 数目决定长度的动态数组。固定数组通常按探测器通道展开保存,动态数组则先记录 hit 数目,再分别保存条带编号、能量等信息。对于这类数据,可以进一步整理为结构体与 vector 的组合形式。
原始数据结构¶
固定数组常写成:
Double_t d1x[32], d2x[32], d3x[32]; // xe
Double_t d1y[32], d2y[32], d3y[32]; // ye
Double_t d1t[32], d2t[32], d3t[32]; // xt
动态数组则常写成:
```cpp
Int_t d1xhit, d2xhit, d3xhit; // hit
Int_t d1yhit, d2yhit, d3yhit;
Int_t d1xs[d1xhit], d2xs[d2xhit], d3xs[d3xhit]; // strip
Double_t d1xe[d1xhit], d2xe[d2xhit], d3xe[d3xhit]; // energy
Int_t d1ys[d1yhit], d2ys[d2yhit], d3ys[d3yhit];
Double_t d1ye[d1yhit], d2ye[d2yhit], d3ye[d3yhit];
如果改用 vector 来组织 hit 信息,可以先定义表示单个 hit 的结构体:
struct dssd
{
Int_t id;
Double_t e;
Double_t t;
};
再分别用 vector<dssd> 保存三层 DSSD 在 x、y 两侧的 hit:
vector<dssd> x1v, x2v, x3v;
vector<dssd> y1v, y2v, y3v;
这样,一个事件中某一侧探测器的全部 hit 就被统一保存在一个 vector 中。相比于将条带编号、能量、时间分别放在不同数组中,这种写法更加紧凑,也更符合后续分析的逻辑。
从普通 ROOT 文件生成含 vector 类型 Branch 的新文件¶
这一部分的目标,是将原始 ROOT 文件中以数组形式保存的 hit 信息,转换为以 vector 类型 Branch 保存的新 ROOT 文件。
通常的流程是:
- 先利用
MakeClass为原始树生成基础分析框架; - 在此基础上补充自己的分析代码;
- 将数组形式的输入数据整理为
vector,并写入新的 ROOT 文件。
程序中一般需要准备如下文件:
main.cppmakefileana.hana.cppLinkdef.h
其中,main.cpp 负责输入输出文件和树的管理;ana.h 与 ana.cpp 负责具体的数据处理;Linkdef.h 和 makefile 用于生成并编译 ROOT 字典。
main.cpp¶
main.cpp 负责打开输入文件 cal_16C.root,读取其中的树,创建输出文件 vec_16C.root 和新的输出树,然后调用分析类中的 Analysis() 完成逐事件处理,最后将输出树写入文件。
#include <iostream>
#include <TFile.h>
#include <TTree.h>
#include <TString.h>
#include "ana.h"
using namespace std;
int main(int argc, char** argv)
{
TString InputPath, OutputPath, infile, outfile;
InputPath = "./";
OutputPath = "./";
infile.Form("%scal_16C.root", InputPath.Data());
outfile.Form("%svec_16C.root", OutputPath.Data());
// input
TFile *ipf = new TFile(infile);
if(!ipf->IsOpen()) {
cout << "Cannot open input file: " << infile << endl;
return -1;
}
TTree *ipt = (TTree*)ipf->Get("tree");
// output
TFile *opf = new TFile(outfile, "RECREATE");
TTree *opt = new TTree("tree", "vector branch");
ana *tk = new ana(ipt, opt);
tk->Analysis();
opt->Write();
ipf->Close();
opf->Close();
return 1;
}
ana.h¶
在 ana.h 中,可以定义分析类 ana。它继承自 MakeClass 生成的基类,并声明六个 vector<dssd> 变量,用于分别保存三层 DSSD 在 x、y 两个方向上的 hit。同时,还定义设置输出树、处理单侧探测器数据以及主分析循环所需的成员函数。
#ifndef ana_h
#define ana_h
#include <TRandom3.h>
#include <iostream>
#include "test.h"
using namespace std;
struct dssd
{
Int_t id;
Double_t e;
Double_t t;
};
class ana : public test
{
public:
vector<dssd> x1v, x2v, x3v;
vector<dssd> y1v, y2v, y3v;
TTree *opt;
TRandom3 *gr;
ana(TTree* ipt_, TTree *opt_) : test(ipt_), opt(opt_) {}
virtual ~ana() {};
virtual void Analysis();
virtual void SetOutBranch();
virtual void ProcessDS(Double_t ee[32], vector<dssd> &vec);
};
#endif
ana.cpp¶
ana.cpp 中主要包括三个部分:
SetOutBranch():设置输出树中的 Branch;ProcessDS():将数组形式的数据整理为vector;Analysis():逐事件读取输入树,并完成转换后写入输出树。
#include "ana.h"
using namespace std;
void ana::SetOutBranch()
{
opt->Branch("x1v", &x1v);
opt->Branch("x2v", &x2v);
opt->Branch("x3v", &x3v);
opt->Branch("y1v", &y1v);
opt->Branch("y2v", &y2v);
opt->Branch("y3v", &y3v);
opt->Branch("sx1e", &sx1e, "sx1e");
opt->Branch("sx2e", &sx2e, "sx2e");
opt->Branch("sx3e", &sx3e, "sx3e");
opt->Branch("sy1e", &sy1e, "sy1e");
opt->Branch("sy2e", &sy2e, "sy2e");
opt->Branch("sy3e", &sy3e, "sy3e");
}
void ana::ProcessDS(Double_t ee[32], vector<dssd> &vec)
{
vec.clear();
dssd ds;
for(int i = 0; i < 32; i++) {
if(ee[i] < 1) continue;
ds.id = i;
ds.e = ee[i];
vec.push_back(ds);
}
}
void ana::Analysis()
{
SetOutBranch();
if (fChain == 0) return;
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;
ProcessDS(d1x, x1v);
ProcessDS(d1y, y1v);
ProcessDS(d2x, x2v);
ProcessDS(d2y, y2v);
ProcessDS(d3x, x3v);
ProcessDS(d3y, y3v);
bool b1 = x1v.size() > 0 || y1v.size() > 0;
bool b2 = x2v.size() > 0 || y2v.size() > 0;
bool b3 = x3v.size() > 0 || y3v.size() > 0;
if (b1 || b2 || b3) opt->Fill();
if ((jentry) % 1000 == 0) {
printf("Process %.2f %% %dk / %dk\r",
Double_t(jentry) / nentries * 100.,
int(jentry / 1000), int(nentries / 1000));
fflush(stdout);
}
}
}
这里的做法是:对每一个事件,分别把 d1x、d1y、d2x、d2y、d3x、d3y 中满足条件的通道转成 dssd 结构体,再压入对应的 vector 中。只要三层 DSSD 中至少有一层存在 hit,就把该事件写入输出树。
字典与编译¶
当 TTree 的 Branch 中保存的是自定义结构体,或 vector<自定义类型> 这类 STL 容器时,需要为相应类型生成 ROOT 字典。否则,ROOT 在读取文件时将无法正确识别这些类型,并可能给出找不到字典的警告信息,从而影响数据的正常读写与后续分析。
例如,在本节中,输出 Branch 中使用了自定义结构体 dssd 以及 vector<dssd>,因此需要在 Linkdef.h 中加入相应的字典声明:
#ifdef __CINT__
#pragma link off all globals;
#pragma link off all classes;
#pragma link off all functions;
#pragma link C++ nestedclasses;
#pragma link C++ class dssd+;
#pragma link C++ class vector<dssd>+;
#endif
在完成 Linkdef.h 的设置后,还需要在编译过程中调用 rootcling 生成字典源文件,并将该文件与主程序及其他源文件一起编译。为此,可以编写如下 makefile:
###############################################################################
OBJ = dssd
MainFile = main.cpp
###############################################################################
LinkFile = LinkDict.cc
###############################################################################
SourceFile := $(LinkFile) $(MainFile) \
$(wildcard $(shell pwd)/src/*.c $(shell pwd)/src/*.cc \
$(shell pwd)/src/*.C $(shell pwd)/src/*.cpp \
$(shell pwd)/src/*.cxx)
IncludeFile := $(wildcard $(shell pwd)/include/*.h \
$(shell pwd)/include/*.hh \
$(shell pwd)/include/*.hpp)
###############################################################################
ROOTCFLAGS = $(shell root-config --cflags)
ROOTLIBS = $(shell root-config --libs)
ROOTGLIBS = $(shell root-config --glibs)
GXX = g++
DIR_INC = -I$(ROOTSYS)/include -I$(shell pwd)/include
CFLAGS = -Wall -O2 $(DIR_INC) -I$(ROOTSYS)/include \
$(ROOTLIBS) -lSpectrum -lXMLParser \
-D_LARGEFILE_SOURCE -D_FILE_OFFSET_BITS=64
###############################################################################
all: $(OBJ)
# 在 Linkdef.h 之前,需要放置包含相关类型声明的头文件
$(LinkFile) : include/ana.h Linkdef.h
rootcling -f $(LinkFile) -c $^
$(OBJ): $(SourceFile)
$(GXX) $(CFLAGS) $(ROOTCFLAGS) $(ROOTLIBS) $(ROOTGLIBS) -o $@ $(SourceFile)
@echo "=============================================================="
@echo "$@ done !"
@echo "=============================================================="
clean:
rm -f *.o *.d $(OBJ) *Dict*
上述 makefile 中:
OBJ指定最终生成的可执行文件名;MainFile指定主程序文件;LinkFile为rootcling生成的字典源文件;SourceFile自动收集src目录下的源文件,并与main.cpp和字典文件一起参与编译;IncludeFile收集include目录中的头文件;ROOTCFLAGS、ROOTLIBS和ROOTGLIBS分别给出 ROOT 所需的编译与链接参数;all目标用于生成最终可执行文件;clean目标用于清除编译过程中产生的中间文件和可执行文件。
这里需要特别注意的是,在调用 rootcling 生成字典时,必须将包含相关类型定义的头文件放在 Linkdef.h 之前。例如本例中使用了 include/ana.h,其中定义了 dssd 以及相关 vector 类型,因此应写成:
$(LinkFile) : include/ana.h Linkdef.h
rootcling -f $(LinkFile) -c $^
如果顺序不正确,rootcling 在处理 Linkdef.h 时将无法识别这些类型,从而导致字典生成失败。也就是说,Linkdef.h 只负责声明需要生成字典的类型,而这些类型本身必须在它之前已经被包含并定义。
在 ROOT 命令行中使用含 vector 类型 Branch 的文件¶
生成含 vector Branch 的新文件后,可以直接在 ROOT 命令行中进行查看和分析。
例如:
TCanvas *c1 = new TCanvas;
TFile *ff = new TFile("vec_16C.root");
Warning in <TClass::Init>: no dictionary for class dssd is available
如果字典没有正确生成,ROOT 在打开文件时可能会给出类似如下的警告:
Warning in <TClass::Init>: no dictionary for class dssd is available
这说明 ROOT 虽然能够打开文件,但对自定义类型的支持并不完整,因此前面的字典生成步骤是必要的。
查看 vector 的长度¶
ROOT 提供了 @vec.size() 的写法,用于获得某个 vector 的长度。例如:
tree->Draw("@x1v.size()");//显示vector的大小
c1->Draw();
tree->Scan("x1v[0].id:x1v[0].e:x1v[0].t:x1v[1].id:x1v[1].e:x1v[1].t",
"@x1v.size()==2", "", 10, 1);
************************************************************************************ * Row * x1v[0].id * x1v[0].e * x1v[0].t * x1v[1].id * x1v[1].e * x1v[1].t * ************************************************************************************ * 1 * 17 * 4925.0811 * 1.988e-81 * 18 * 1007.1224 * 0 * * 2 * 23 * 3146.8572 * 1.988e-81 * 24 * 339.14570 * 0 * * 3 * 23 * 5465.0851 * 1.988e-81 * 24 * 41.435740 * 0 * * 4 * 21 * 5993.6417 * 1.988e-81 * 22 * 43.867701 * 0 * * 5 * 11 * 4196.0801 * 1.988e-81 * 12 * 1159.4884 * 0 * * 6 * 13 * 387.77402 * 1.988e-81 * 14 * 3889.7610 * 0 * ************************************************************************************ ==> 6 selected entries
tree->Scan("x1v.id:x1v.e:x1v.t", "@x1v.size()==2", "", 10, 1);
*********************************************************** * Row * Instance * x1v.id * x1v.e * x1v.t * *********************************************************** * 1 * 0 * 17 * 4925.0811 * 1.988e-81 * * 1 * 1 * 18 * 1007.1224 * 0 * * 2 * 0 * 23 * 3146.8572 * 1.988e-81 * * 2 * 1 * 24 * 339.14570 * 0 * * 3 * 0 * 23 * 5465.0851 * 1.988e-81 * * 3 * 1 * 24 * 41.435740 * 0 * * 4 * 0 * 21 * 5993.6417 * 1.988e-81 * * 4 * 1 * 22 * 43.867701 * 0 * * 5 * 0 * 11 * 4196.0801 * 1.988e-81 * * 5 * 1 * 12 * 1159.4884 * 0 * * 6 * 0 * 13 * 387.77402 * 1.988e-81 * * 6 * 1 * 14 * 3889.7610 * 0 * *********************************************************** ==> 12 selected entries
在这种情况下,ROOT 会自动把 vector 中的每个元素展开显示,并为每个元素给出对应的 Instance 编号。这样可以方便地查看一个事件中所有 hit 的详细信息。
在后续分析中继续保持 vector 类型 Branch¶
如果已经生成了含 vector Branch 的 ROOT 文件,后续分析时往往希望继续以 vector 的形式读取这些数据,而不是重新退回到数组形式。此时需要特别注意:不能直接沿用 MakeClass 的默认处理方式。
原因在于,MakeClass 对 vector Branch 的展开方式更接近于动态数组。也就是说,虽然树中原本保存的是 vector<dssd>,但 MakeClass 生成的代码会将其解释为长度变量和若干数组成员的组合。例如,执行:
tree->Print();
****************************************************************************** *Tree :tree : vector branch * *Entries : 2023076 : Total = 744983759 bytes File Size = 335786363 * * : : Tree compression factor = 2.22 * ****************************************************************************** *Br 0 :x1v : Int_t x1v_ * *Entries : 2023076 : Total Size= 16284062 bytes File Size = 3881644 * *Baskets : 764 : Basket Size= 32000 bytes Compression= 4.18 * *............................................................................* *Br 1 :x1v.id : Int_t id[x1v_] * *Entries : 2023076 : Total Size= 25352263 bytes File Size = 6320782 * *Baskets : 98 : Basket Size= 3794944 bytes Compression= 4.01 * *............................................................................* *Br 2 :x1v.e : Double_t e[x1v_] * *Entries : 2023076 : Total Size= 42606133 bytes File Size = 32524001 * *Baskets : 141 : Basket Size= 5379584 bytes Compression= 1.31 * *............................................................................* *Br 3 :x1v.t : Double_t t[x1v_] * *Entries : 2023076 : Total Size= 42606133 bytes File Size = 13246169 * *Baskets : 141 : Basket Size= 5379584 bytes Compression= 3.22 * *............................................................................* *Br 4 :x2v : Int_t x2v_ * *Entries : 2023076 : Total Size= 16285402 bytes File Size = 4196636 * *Baskets : 764 : Basket Size= 32000 bytes Compression= 3.87 * *............................................................................* *Br 5 :x2v.id : Int_t id[x2v_] * *Entries : 2023076 : Total Size= 28880153 bytes File Size = 8438087 * *Baskets : 112 : Basket Size= 4293120 bytes Compression= 3.42 * *............................................................................* *Br 6 :x2v.e : Double_t e[x2v_] * *Entries : 2023076 : Total Size= 49661789 bytes File Size = 38240269 * *Baskets : 168 : Basket Size= 6375936 bytes Compression= 1.30 * *............................................................................* *Br 7 :x2v.t : Double_t t[x2v_] * *Entries : 2023076 : Total Size= 49661789 bytes File Size = 17650173 * *Baskets : 168 : Basket Size= 6375936 bytes Compression= 2.81 * *............................................................................* *Br 8 :x3v : Int_t x3v_ * *Entries : 2023076 : Total Size= 16282902 bytes File Size = 3994156 * *Baskets : 764 : Basket Size= 32000 bytes Compression= 4.07 * *............................................................................* *Br 9 :x3v.id : Int_t id[x3v_] * *Entries : 2023076 : Total Size= 20228212 bytes File Size = 6121364 * *Baskets : 87 : Basket Size= 3377664 bytes Compression= 3.30 * *............................................................................* *Br 10 :x3v.e : Double_t e[x3v_] * *Entries : 2023076 : Total Size= 32357957 bytes File Size = 23241034 * *Baskets : 118 : Basket Size= 4545024 bytes Compression= 1.39 * *............................................................................* *Br 11 :x3v.t : Double_t t[x3v_] * *Entries : 2023076 : Total Size= 32357957 bytes File Size = 10010541 * *Baskets : 118 : Basket Size= 4545024 bytes Compression= 3.23 * *............................................................................* *Br 12 :x1v : Int_t x1v_ * *Entries : 2023076 : Total Size= 16284042 bytes File Size = 3881644 * *Baskets : 764 : Basket Size= 32000 bytes Compression= 4.18 * *............................................................................* *Br 13 :x1v.id : Int_t id[x1v_] * *Entries : 2023076 : Total Size= 25352259 bytes File Size = 6320782 * *Baskets : 98 : Basket Size= 3794944 bytes Compression= 4.01 * *............................................................................* *Br 14 :x1v.e : Double_t e[x1v_] * *Entries : 2023076 : Total Size= 42606129 bytes File Size = 32524001 * *Baskets : 141 : Basket Size= 5379584 bytes Compression= 1.31 * *............................................................................* *Br 15 :x1v.t : Double_t t[x1v_] * *Entries : 2023076 : Total Size= 42606129 bytes File Size = 13246169 * *Baskets : 141 : Basket Size= 5379584 bytes Compression= 3.22 * *............................................................................* *Br 16 :x2v : Int_t x2v_ * *Entries : 2023076 : Total Size= 16285402 bytes File Size = 4196636 * *Baskets : 764 : Basket Size= 32000 bytes Compression= 3.87 * *............................................................................* *Br 17 :x2v.id : Int_t id[x2v_] * *Entries : 2023076 : Total Size= 28880153 bytes File Size = 8438087 * *Baskets : 112 : Basket Size= 4293120 bytes Compression= 3.42 * *............................................................................* *Br 18 :x2v.e : Double_t e[x2v_] * *Entries : 2023076 : Total Size= 49661789 bytes File Size = 38240269 * *Baskets : 168 : Basket Size= 6375936 bytes Compression= 1.30 * *............................................................................* *Br 19 :x2v.t : Double_t t[x2v_] * *Entries : 2023076 : Total Size= 49661789 bytes File Size = 17650173 * *Baskets : 168 : Basket Size= 6375936 bytes Compression= 2.81 * *............................................................................* *Br 20 :x3v : Int_t x3v_ * *Entries : 2023076 : Total Size= 16282902 bytes File Size = 3994156 * *Baskets : 764 : Basket Size= 32000 bytes Compression= 4.07 * *............................................................................* *Br 21 :x3v.id : Int_t id[x3v_] * *Entries : 2023076 : Total Size= 20228212 bytes File Size = 6121364 * *Baskets : 87 : Basket Size= 3377664 bytes Compression= 3.30 * *............................................................................* *Br 22 :x3v.e : Double_t e[x3v_] * *Entries : 2023076 : Total Size= 32357957 bytes File Size = 23241034 * *Baskets : 118 : Basket Size= 4545024 bytes Compression= 1.39 * *............................................................................* *Br 23 :x3v.t : Double_t t[x3v_] * *Entries : 2023076 : Total Size= 32357957 bytes File Size = 10010541 * *Baskets : 118 : Basket Size= 4545024 bytes Compression= 3.23 * *............................................................................*
在后续分析中继续保持 vector 类型 Branch¶
如果已经生成了含 vector Branch 的 ROOT 文件,后续分析时往往希望继续以 vector 的形式读取这些数据,而不是重新退回到数组形式。此时需要特别注意:不能直接沿用 MakeClass 的默认处理方式。
原因在于,MakeClass 对 vector Branch 的展开方式更接近于动态数组。也就是说,虽然树中原本保存的是 vector<dssd>,但 MakeClass 生成的代码会将其解释为长度变量和若干数组成员的组合。例如,执行:
tree->Print();
可以看到 x1v 之类的 Branch 被展开为:
x1v_x1v.id[x1v_]x1v.e[x1v_]x1v.t[x1v_]
这种展开方式虽然可以访问数据,但不利于继续使用 STL 中针对 vector 的各种操作。
因此,若希望在分析代码中继续把输入量视为 vector<dssd>,就需要手动定义指针,并使用 SetBranchAddress() 建立关联。
输入 Branch 的定义¶
在新的分析类中,可以将输入 vector Branch 定义为指针,同时定义一个新的结构体用于保存 x-y 配对后的结果:
#ifndef ana_h
#define ana_h
#include <vector>
#include <algorithm>
#include <TFile.h>
#include <TTree.h>
using namespace std;
struct dssd
{
Int_t id;
Double_t e;
Double_t t;
};
struct DSSD
{
int xid;
int yid;
double xe;
double ye;
};
class ana
{
public:
vector<dssd> *br_x1v, *br_x2v, *br_x3v;
vector<dssd> *br_y1v, *br_y2v, *br_y3v;
Double_t sx1e, sx2e, sx3e;
Double_t sy1e, sy2e, sy3e;
TTree *ipt;
TTree *opt;
vector<DSSD> d1, d2, d3;
ana(TTree* ipt_, TTree *opt_) : ipt(ipt_), opt(opt_) {}
virtual ~ana() {};
virtual void SetBranchInput();
virtual void GetDSSD(vector<dssd> *x, vector<dssd> *y, vector<DSSD> &xy);
virtual void Analysis();
virtual void BranchOutput();
};
#endif
设置输入 Branch 地址¶
在使用 SetBranchAddress() 之前,应先将这些指针初始化为空:
#include "ana.h"
using namespace std;
void ana::SetBranchInput()
{
br_x1v = NULL;
br_x2v = NULL;
br_x3v = NULL;
br_y1v = NULL;
br_y2v = NULL;
br_y3v = NULL;
ipt->SetBranchAddress("x1v", &br_x1v);
ipt->SetBranchAddress("x2v", &br_x2v);
ipt->SetBranchAddress("x3v", &br_x3v);
ipt->SetBranchAddress("y1v", &br_y1v);
ipt->SetBranchAddress("y2v", &br_y2v);
ipt->SetBranchAddress("y3v", &br_y3v);
ipt->SetBranchAddress("sx1e", &sx1e);
ipt->SetBranchAddress("sx2e", &sx2e);
ipt->SetBranchAddress("sx3e", &sx3e);
ipt->SetBranchAddress("sy1e", &sy1e);
ipt->SetBranchAddress("sy2e", &sy2e);
ipt->SetBranchAddress("sy3e", &sy3e);
}
这一步非常重要。若不先初始化为空指针,程序在运行时可能会出现错误。
设置输出 Branch¶
输出时,可以将配对后的结果保存为新的 vector<DSSD> Branch:
void ana::BranchOutput()
{
opt->Branch("d1", &d1);
opt->Branch("d2", &d2);
opt->Branch("d3", &d3);
}
排序与配对¶
由于同一侧探测器可能存在多个 hit,在进行 x-y 关联之前,常常需要先对 hit 按能量进行排序。例如,可以定义一个比较函数:
bool SortDS(dssd &a, dssd &b)
{
return a.e > b.e;
}
随后对各 vector 分别排序:
sort(br_x1v->begin(), br_x1v->end(), SortDS);
sort(br_y1v->begin(), br_y1v->end(), SortDS);
在完成排序后,就可以进行 x-y 配对。一个简单的思路是:
- 取 x、y 两侧 hit 数目的较小值作为循环上限;
- 分别取出对应 hit 的条带编号和能量;
- 判断 x、y 两侧能量是否匹配;
- 若满足条件,则构造一个
DSSD结果并保存。
示意代码如下:
void ana::GetDSSD(vector<dssd> *x, vector<dssd> *y, vector<DSSD> &xy)
{
int hit = min(x->size(), y->size());
xy.clear();
DSSD dxy;
for (int i = 0; i < hit; i++) {
double xe = (*x)[i].e;
double ye = (*y)[i].e;
int ix = (*x)[i].id;
int iy = (*y)[i].id;
if (abs(xe - ye) < 50) {
dxy.xe = xe;
dxy.ye = ye;
dxy.xid = ix;
dxy.yid = iy;
xy.push_back(dxy);
}
}
}
这里的能量匹配条件、配对策略以及排序方式都可以根据具体实验数据进行调整。
后续分析主循环¶
在 Analysis() 中,完整流程通常包括:
- 设置输入 Branch;
- 设置输出 Branch;
- 逐事件读取数据;
- 对各层探测器的 x、y hit 分别排序;
- 进行 x-y 配对;
- 将得到的结果写入新的输出树。
void ana::Analysis()
{
SetBranchInput();
BranchOutput();
if (ipt == 0) return;
Long64_t nentries = ipt->GetEntriesFast();
for (Long64_t jentry = 0; jentry < nentries; jentry++) {
ipt->GetEntry(jentry);
sort(br_x1v->begin(), br_x1v->end(), SortDS);
sort(br_y1v->begin(), br_y1v->end(), SortDS);
sort(br_x2v->begin(), br_x2v->end(), SortDS);
sort(br_y2v->begin(), br_y2v->end(), SortDS);
sort(br_x3v->begin(), br_x3v->end(), SortDS);
sort(br_y3v->begin(), br_y3v->end(), SortDS);
GetDSSD(br_x1v, br_y1v, d1);
GetDSSD(br_x2v, br_y2v, d2);
GetDSSD(br_x3v, br_y3v, d3);
if (d1.size() > 0 || d2.size() > 0 || d3.size() > 0)
opt->Fill();
if ((jentry) % 1000 == 0) {
printf("Process %.2f %% %dk / %dk\r",
Double_t(jentry) / nentries * 100.,
int(jentry / 1000), int(nentries / 1000));
fflush(stdout);
}
}
}
后续分析时的 main.cpp 与前面类似,只是输入文件变为 vec_16C.root,输出文件变为 sort_16C.root:
#include <iostream>
#include <TFile.h>
#include <TTree.h>
#include <TString.h>
#include "ana.h"
using namespace std;
int main(int argc, char** argv)
{
TString InputPath, OutputPath, infile, outfile;
InputPath = "./";
OutputPath = "./";
infile.Form("%svec_16C.root", InputPath.Data());
outfile.Form("%ssort_16C.root", OutputPath.Data());
// input
TFile *ipf = new TFile(infile);
if(!ipf->IsOpen()) {
cout << "Cannot open input file: " << infile << endl;
return -1;
}
TTree *ipt = (TTree*)ipf->Get("tree");
// output
TFile *opf = new TFile(outfile, "RECREATE");
TTree *opt = new TTree("tree", "vector branch");
ana *tk = new ana(ipt, opt);
tk->Analysis();
opt->Write();
ipf->Close();
opf->Close();
return 1;
}
如果输出中又引入了新的结构体 DSSD 以及 vector<DSSD>,那么也需要继续在 Linkdef.h 中补充对应的字典声明:
#ifdef __CINT__
#pragma link off all globals;
#pragma link off all classes;
#pragma link off all functions;
#pragma link C++ nestedclasses;
#pragma link C++ class dssd+;
#pragma link C++ class vector<dssd>+;
#pragma link C++ class DSSD+;
#pragma link C++ class vector<DSSD>+;
#endif
实际代码¶
本节对应的完整可运行代码放在 GitHub 仓库的 chapt3/code 目录下,其中包含 code1 和 code2 两个子目录。code1 的主程序读取 cal_16C.root,输出 vec_16C.root;code2 的主程序读取 vec_16C.root,输出 sort_16C.root。这两个目录对应前后两个连续的分析步骤。
代码目录:
https://github.com/zhihuanli/Experimental-Data-Analysis-Course/tree/master/chapt3/code (GitHub)
code1 用于将原始 ROOT 文件中的数组形式数据整理成 vector 类型的 Branch,并写入新的 ROOT 文件。对应讲义中“从普通 ROOT 文件生成含 vector 类型 Branch 的新文件”这一部分。code1 中的 main.cpp 明确给出了输入文件 cal_16C.root 和输出文件 vec_16C.root,而 src/ana.cpp 中则完成了 x1v、x2v、x3v、y1v、y2v、y3v 等 vector Branch 的建立与填充。
code2 用于继续读取已经包含 vector Branch 的 ROOT 文件,在分析中保持 vector 类型读入,完成 x-y hit 的排序、配对,并将结果写入新的 ROOT 文件。对应讲义中“使用含 vector 类型 Branch 的 ROOT 文件继续进行分析”这一部分。code2 的 main.cpp 给出了输入文件 vec_16C.root 和输出文件 sort_16C.root,而 src/ana.cpp 中则通过 SetBranchAddress() 读取 vector<dssd>,并将配对结果输出为 d1、d2、d3 等新的 Branch。
vector 的常用操作¶
在实际分析中,vector 的优势不仅在于可以保存变长数据,更重要的是可以直接利用 STL 提供的各种操作完成排序、筛选和去重等任务。
常见的 vector 成员函数包括:
vector::size(); // 返回 vector 大小
vector::clear(); // 清空 vector
vector::push_back(value); // 在 vector 内填入元素
vector::begin(); // 返回 vector 的首地址
vector::end(); // 返回 vector 的尾地址
vector::assign(v1.begin(), v2.end()); // 将 v1 拷贝到 vector 中
vector::erase(it); // 删除 it 位置变量,返回下一个元素的指针
排序¶
若要按能量从高到低排序,可结合比较函数与 sort() 使用:
bool SortDS(dssd &a, dssd &b)
{
return a.e > b.e;
}
sort(x1v.begin(), x1v.end(), SortDS);
删除重复元素¶
如果需要删除重复 hit,可以先定义“重复”的判据,再结合 unique() 与 erase() 使用:
bool Equal(dssd &a, dssd &b)
{
return Long64_t(a.e - b.e) == 0 && a.id == b.id;
}
void Unique(vector<dssd> &a)
{
a.erase(unique(a.begin(), a.end(), Equal), a.end());
}
// remove identical events, probably caused by bugs in Digital modules
Unique(xvec);
按条件删除元素¶
在很多分析中,需要从 vector 中删除不满足条件的元素。例如,当某个 hit 的时间与参考时间相差过大时,可以将其剔除:
void tCut(vector<dssd> &a, vector<dssd> &b, double t1, double t2)
{
if (b.size() > 0) {
for (auto it = a.begin(); it != a.end(); ) {
int dt = Long64_t(it->t - b[0].t);
if (dt < t1 || dt > t2)
it = a.erase(it);
else
++it;
}
}
}
tCut(x1v, x1v, -20, 20);
这些例子说明,使用 vector 的优势不仅在于能够保存变长 hit 信息,还在于可以直接利用 STL 提供的排序、去重和筛选等操作,使后续分析实现起来更加自然。