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

3.8 用 vector 保存 TTree 数据

在探测器数据分析中，一个事件中往往会出现多个 hit。若直接使用固定长度数组，或依赖 hit 数目的动态数组保存这些信息，虽然能够完成数据存储，但在后续分析中往往不便于按单个 hit 进行排序、筛选和配对。利用 vector，可以把同一个 hit 的条带编号、能量和时间统一组织起来，再将一个事件中的全部 hit 保存在同一个容器中，从而使数据结构更加清晰，也更便于后续分析处理。

以 data_16C.root 为例，原始 ROOT 文件中既包含固定长度数组，也包含由 hit 数目决定长度的动态数组。固定数组通常按探测器通道展开保存，动态数组则先记录 hit 数目，再分别保存条带编号、能量等信息。对于这类数据，可以进一步整理为结构体与 vector 的组合形式。

原始数据结构

固定数组常写成：

Double_t d1x[32], d2x[32], d3x[32];
Double_t d1y[32], d2y[32], d3y[32];

可变长 Branch 则由 hit 数控制写入长度；内存中仍保留足够的固定容量，例如一层的 X 面：

Int_t d1xhit;
Int_t d1xs[32];
Double_t d1xe[32];
// tree->Branch("d1xe", d1xe, "d1xe[d1xhit]/D");

如果改用 vector 来组织 hit 信息，可以先定义表示单个 hit 的结构体：

struct dssd
{
  Int_t id;
  Double_t e;
  Double_t t = std::numeric_limits<double>::quiet_NaN();
};

再分别用 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.cpp
• makefile
• ana.h
• ana.cpp
• Linkdef.h

其中，main.cpp 负责输入输出文件和树的管理；ana.h 与 ana.cpp 负责具体的数据处理；Linkdef.h 和 makefile 用于生成并编译 ROOT 字典。

main.cpp

main.cpp 负责打开输入文件 strip_arrays_16C.root，读取其中的树，创建输出文件 vec_16C.root 和新的输出树，然后调用分析类中的 Analysis() 完成逐事件处理，最后将输出树写入文件。

#include <TFile.h>
#include <TTree.h>
#include <iostream>
#include "ana.h"
int main(int argc, char** argv) {
    const char* inputName = argc>1 ? argv[1] : "../../data/strip_arrays_16C.root";
    const char* outputName = argc>2 ? argv[2] : "../../vec_16C.root";
    TFile* input = TFile::Open(inputName);
    if (!input || input->IsZombie()) return 1;
    TTree* tin = input->Get<TTree>("tree");
    if (!tin) return 1;
    TFile output(outputName,"RECREATE");
    TTree* tout = new TTree("tree","vector branch");
    {
        ana analysis(tin,tout);
        analysis.Analysis();
        std::cout << "Input=" << tin->GetEntries() << ", output=" << tout->GetEntries() << '\n';
        output.cd();
        tout->Write();
    }
    // MakeClass 基类的析构函数已释放 input。
    return 0;
}

ana.h

在 ana.h 中，可以定义分析类 ana。它继承自 MakeClass 生成的基类，并声明六个 vector<dssd> 变量，用于分别保存三层 DSSD 在 x、y 两个方向上的 hit。同时，还定义设置输出树、处理单侧探测器数据以及主分析循环所需的成员函数。

#ifndef ana_h
#define ana_h

include <vector>
include <limits>
include <iostream>
include "test.h"  //包含基类头文件
using namespace std;
struct dssd
{
Int_t id;
Double_t e;
Double_t t = std::numeric_limits<double>::quiet_NaN();
};
class ana : public test //从test类中继承其成员变量和成员函数
{
public:
vector<dssd> x1v,x2v,x3v;
vector<dssd> y1v,y2v,y3v;
Long64_t source_entry = 0;
TTree *opt;
ana(TTree* ipt_,TTree *opt_): test(ipt_),opt(opt_) {}
virtual ~ana() {};
virtual void     Analysis();//分析函数，作用等价于原Loop函数
virtual void     SetOutBranch();
virtual void     ProcessDS(const 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("source_entry", &source_entry, "source_entry/L");
  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/D");
  opt->Branch("sx2e",&sx2e,"sx2e/D");
  opt->Branch("sx3e",&sx3e,"sx3e/D");
  opt->Branch("sy1e",&sy1e,"sy1e/D");
  opt->Branch("sy2e",&sy2e,"sy2e/D");
  opt->Branch("sy3e",&sy3e,"sy3e/D");

}
void ana::ProcessDS(const Double_t ee[32], vector<dssd> &vec)
{
vec.clear(); // 每个事件重新建立 hit 列表
for(int i=0; i<32; ++i) {
if(ee[i]<1) continue;
dssd hit;
hit.id=i;
hit.e=ee[i];
vec.push_back(hit); // t 保持 NaN，输入没有时间信息
}
}
void ana::Analysis()
{
if (fChain == 0) return;
SetOutBranch();
Long64_t nentries = fChain->GetEntriesFast();
for (Long64_t jentry=0; jentry<nentries;jentry++) {
Long64_t ientry = LoadTree(jentry);
if (ientry < 0) break;
fChain->GetEntry(jentry);
ProcessDS(d1x,x1v);
ProcessDS(d1y,y1v);
ProcessDS(d2x,x2v);
ProcessDS(d2y,y2v);
ProcessDS(d3x,x3v);
ProcessDS(d3y,y3v);
source_entry = jentry;
opt->Fill(); // 保留空事件及事件对应关系
}
}

这里的做法是：对每一个事件，分别把 d1x、d1y、d2x、d2y、d3x、d3y 中满足条件的通道转成 dssd 结构体，再压入对应的 vector 中。每个输入事件都写入输出树，包括空 vector；这样后续仍能追溯原事件。

---

字典与编译

当 TTree 的 Branch 中保存的是自定义结构体，或 vector<自定义类型> 这类 STL 容器时，需要为相应类型生成 ROOT 字典。否则，ROOT 在读取文件时将无法正确识别这些类型，并可能给出找不到字典的警告信息，从而影响数据的正常读写与后续分析。

例如，在本节中，输出 Branch 中使用了自定义结构体 dssd 以及 vector<dssd>，因此需要在 Linkdef.h 中加入相应的字典声明：

#ifdef __CLING__

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：

CXX = c++
CPPFLAGS = -Iinclude $(shell root-config --cflags)
CXXFLAGS = -O2 -Wall
LDLIBS = $(shell root-config --libs)
SOURCES = main.cpp $(wildcard src/*.cpp src/*.C) LinkDict.cc
HEADERS = $(wildcard include/*.h)

all: dssd libhits.so
LinkDict.cc: $(HEADERS) Linkdef.h
	rootcling -f $@ -Iinclude include/ana.h Linkdef.h
dssd: $(SOURCES) $(HEADERS)
$(CXX) $(CPPFLAGS) $(CXXFLAGS) $(SOURCES) $(LDLIBS) -o $@
libhits.so: LinkDict.cc $(HEADERS)
	$(CXX) $(CPPFLAGS) $(CXXFLAGS) -fPIC -shared LinkDict.cc $(LDLIBS) -o $@
clean:
rm -f libhits.so dssd LinkDict.cc LinkDict_rdict.pcm

Makefile 中，SOURCES 收集主程序、分析源文件和生成的 LinkDict.cc；HEADERS 作为依赖。CPPFLAGS、CXXFLAGS、LDLIBS 分别提供头文件路径、编译选项与链接库。它同时构建独立程序 dssd 和供 ROOT 会话载入的 libhits.so。

• all 构建 dssd 和 libhits.so；clean 只清理本工程的构建产物。
• clean 目标用于清除编译过程中产生的中间文件和可执行文件。

这里需要特别注意的是，在调用 rootcling 生成字典时，必须将包含相关类型定义的头文件放在 Linkdef.h 之前。例如本例中使用了 include/ana.h，其中定义了 dssd 以及相关 vector 类型，因此应写成：

LinkDict.cc: $(HEADERS) Linkdef.h
	rootcling -f $@ -Iinclude include/ana.h Linkdef.h

如果顺序不正确，rootcling 在处理 Linkdef.h 时将无法识别这些类型，从而导致字典生成失败。也就是说，Linkdef.h 只负责声明需要生成字典的类型，而这些类型本身必须在它之前已经被包含并定义。

---

在 ROOT 命令行中使用含 vector 类型 Branch 的文件

生成含 vector Branch 的新文件后，可以直接在 ROOT 命令行中进行查看和分析。

例如：

本节输入为 data/strip_arrays_16C.root，保留刻度后的固定数组，区别于 3.6 的 compact-hit 文件。它没有时间数据，结构体中的 t 显式设为 NaN，不生成虚构时间。两个转换程序均保留原事件编号和空事件。

%jsroot on

在 ROOT 命令行中使用含 vector 类型 Branch 的文件

生成含 vector Branch 的新文件后，可以直接在 ROOT 命令行中进行查看和分析。

例如：

gSystem->Load("code/code1/libhits.so"); // 本节 make 同时生成读取字典
TCanvas *c1 = new TCanvas;
TFile *ff = new TFile("vec_16C.root");
if (!ff || ff->IsZombie()) throw std::runtime_error("无法打开输入 ROOT 文件");
TTree *tree = (TTree*)ff->Get("tree");
if (!tree) throw std::runtime_error("输入文件中缺少 tree");

如果字典没有正确生成，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();

该语句可以直接画出 x1v 在各事件中的 hit 数分布。

访问 vector 中单个元素的成员

可以通过

vec[i].member

的形式访问 vector 中第 i 个元素的成员变量。例如：

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 *       nan *        18 * 1007.1224 *       nan *
*        2 *        23 * 3146.8572 *       nan *        24 * 339.14570 *       nan *
*        3 *        23 * 5465.0851 *       nan *        24 * 41.435740 *       nan *
*        4 *        21 * 5993.6417 *       nan *        22 * 43.867701 *       nan *
*        5 *        11 * 4196.0801 *       nan *        12 * 1159.4884 *       nan *
*        6 *        13 * 387.77402 *       nan *        14 * 3889.7610 *       nan *
************************************************************************************
==> 6 selected entries

表示查看 x1v 中前两个 hit 的条带编号、能量和时间。

展开查看全部 hit 信息

除了按下标逐项访问之外，也可以直接写成：

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 *       nan *
*        1 *        1 *        18 * 1007.1224 *       nan *
*        2 *        0 *        23 * 3146.8572 *       nan *
*        2 *        1 *        24 * 339.14570 *       nan *
*        3 *        0 *        23 * 5465.0851 *       nan *
*        3 *        1 *        24 * 41.435740 *       nan *
*        4 *        0 *        21 * 5993.6417 *       nan *
*        4 *        1 *        22 * 43.867701 *       nan *
*        5 *        0 *        11 * 4196.0801 *       nan *
*        5 *        1 *        12 * 1159.4884 *       nan *
*        6 *        0 *        13 * 387.77402 *       nan *
*        6 *        1 *        14 * 3889.7610 *       nan *
***********************************************************
==> 12 selected entries

在这种情况下，ROOT 会自动把 vector 中的每个元素展开显示，并为每个元素给出对应的 Instance 编号。这样可以方便地查看一个事件中所有 hit 的详细信息。

在后续分析中继续保持 vector 类型 Branch

如果已经生成了含 vector Branch 的 ROOT 文件，后续分析时往往希望继续以 vector 的形式读取这些数据，而不是重新退回到数组形式。此时需要特别注意：不能直接沿用 MakeClass 的默认处理方式。

对本例中 split 的 vector<dssd> Branch，MakeClass 会按长度与成员数组生成读取代码。这种方式可以访问数据，但不直接提供整个 vector 的容器操作。需要继续按 vector 读入时，可用下面的对象指针绑定方式。

tree->Print();

******************************************************************************
*Tree    :tree      : vector branch                                          *
*Entries :  1926502 : Total =       825074564 bytes  File  Size =  342577290 *
*        :          : Tree compression factor =   2.41                       *
******************************************************************************
*Br    0 :source_entry : source_entry/L                                      *
*Entries :  1926502 : Total  Size=   15417703 bytes  File Size  =    2963374 *
*Baskets :       54 : Basket Size=    1520128 bytes  Compression=   5.20     *
*............................................................................*
*Br    1 :x1v       : Int_t x1v_                                             *
*Entries :  1926502 : Total  Size=   15507888 bytes  File Size  =    3614131 *
*Baskets :      735 : Basket Size=      32000 bytes  Compression=   4.28     *
*............................................................................*
*Br    2 :x1v.id    : Int_t id[x1v_]                                         *
*Entries :  1926502 : Total  Size=   24661757 bytes  File Size  =    6135902 *
*Baskets :       96 : Basket Size=    3631104 bytes  Compression=   4.02     *
*............................................................................*
*Br    3 :x1v.e     : Double_t e[x1v_]                                       *
*Entries :  1926502 : Total  Size=   41611521 bytes  File Size  =   35194103 *
*Baskets :      138 : Basket Size=    5157376 bytes  Compression=   1.18     *
*............................................................................*
*Br    4 :x1v.t     : Double_t t[x1v_]                                       *
*Entries :  1926502 : Total  Size=   41611521 bytes  File Size  =    3061385 *
*Baskets :      138 : Basket Size=    5157376 bytes  Compression=  13.59     *
*............................................................................*
*Br    5 :x2v       : Int_t x2v_                                             *
*Entries :  1926502 : Total  Size=   15509108 bytes  File Size  =    3988457 *
*Baskets :      735 : Basket Size=      32000 bytes  Compression=   3.88     *
*............................................................................*
*Br    6 :x2v.id    : Int_t id[x2v_]                                         *
*Entries :  1926502 : Total  Size=   27790829 bytes  File Size  =    8149989 *
*Baskets :      108 : Basket Size=    4084736 bytes  Compression=   3.41     *
*............................................................................*
*Br    7 :x2v.e     : Double_t e[x2v_]                                       *
*Entries :  1926502 : Total  Size=   47869745 bytes  File Size  =   41377037 *
*Baskets :      163 : Basket Size=    6064640 bytes  Compression=   1.16     *
*............................................................................*
*Br    8 :x2v.t     : Double_t t[x2v_]                                       *
*Entries :  1926502 : Total  Size=   47869745 bytes  File Size  =    3326076 *
*Baskets :      163 : Basket Size=    6064640 bytes  Compression=  14.39     *
*............................................................................*
*Br    9 :x3v       : Int_t x3v_                                             *
*Entries :  1926502 : Total  Size=   15506668 bytes  File Size  =    3801045 *
*Baskets :      735 : Basket Size=      32000 bytes  Compression=   4.07     *
*............................................................................*
*Br   10 :x3v.id    : Int_t id[x3v_]                                         *
*Entries :  1926502 : Total  Size=   19228061 bytes  File Size  =    5806670 *
*Baskets :       84 : Basket Size=    3197952 bytes  Compression=   3.31     *
*............................................................................*
*Br   11 :x3v.e     : Double_t e[x3v_]                                       *
*Entries :  1926502 : Total  Size=   30744157 bytes  File Size  =   24490887 *
*Baskets :      114 : Basket Size=    4291584 bytes  Compression=   1.26     *
*............................................................................*
*Br   12 :x3v.t     : Double_t t[x3v_]                                       *
*Entries :  1926502 : Total  Size=   30744157 bytes  File Size  =    2681812 *
*Baskets :      114 : Basket Size=    4291584 bytes  Compression=  11.46     *
*............................................................................*
*Br   13 :y1v       : Int_t y1v_                                             *
*Entries :  1926502 : Total  Size=   15507908 bytes  File Size  =    4073264 *
*Baskets :      735 : Basket Size=      32000 bytes  Compression=   3.80     *
*............................................................................*
*Br   14 :y1v.id    : Int_t id[y1v_]                                         *
*Entries :  1926502 : Total  Size=   24196297 bytes  File Size  =    6776205 *
*Baskets :       96 : Basket Size=    3644416 bytes  Compression=   3.57     *
*............................................................................*
*Br   15 :y1v.e     : Double_t e[y1v_]                                       *
*Entries :  1926502 : Total  Size=   40680705 bytes  File Size  =   34401533 *
*Baskets :      139 : Basket Size=    5184512 bytes  Compression=   1.18     *
*............................................................................*
*Br   16 :y1v.t     : Double_t t[y1v_]                                       *
*Entries :  1926502 : Total  Size=   40680705 bytes  File Size  =    3161987 *
*Baskets :      139 : Basket Size=    5184512 bytes  Compression=  12.86     *
*............................................................................*
*Br   17 :y2v       : Int_t y2v_                                             *
*Entries :  1926502 : Total  Size=   15509068 bytes  File Size  =    4135907 *
*Baskets :      735 : Basket Size=      32000 bytes  Compression=   3.74     *
*............................................................................*
*Br   18 :y2v.id    : Int_t id[y2v_]                                         *
*Entries :  1926502 : Total  Size=   27906997 bytes  File Size  =    8156518 *
*Baskets :      108 : Basket Size=    4072960 bytes  Compression=   3.42     *
*............................................................................*
*Br   19 :y2v.e     : Double_t e[y2v_]                                       *
*Entries :  1926502 : Total  Size=   48101981 bytes  File Size  =   41621527 *
*Baskets :      162 : Basket Size=    6041600 bytes  Compression=   1.16     *
*............................................................................*
*Br   20 :y2v.t     : Double_t t[y2v_]                                       *
*Entries :  1926502 : Total  Size=   48101981 bytes  File Size  =    3336585 *
*Baskets :      162 : Basket Size=    6041600 bytes  Compression=  14.42     *
*............................................................................*
*Br   21 :y3v       : Int_t y3v_                                             *
*Entries :  1926502 : Total  Size=   15506768 bytes  File Size  =    3829699 *
*Baskets :      735 : Basket Size=      32000 bytes  Compression=   4.04     *
*............................................................................*
*Br   22 :y3v.id    : Int_t id[y3v_]                                         *
*Entries :  1926502 : Total  Size=   19574150 bytes  File Size  =    5895363 *
*Baskets :       85 : Basket Size=    3237888 bytes  Compression=   3.32     *
*............................................................................*
*Br   23 :y3v.e     : Double_t e[y3v_]                                       *
*Entries :  1926502 : Total  Size=   31436333 bytes  File Size  =   25193096 *
*Baskets :      116 : Basket Size=    4371456 bytes  Compression=   1.25     *
*............................................................................*
*Br   24 :y3v.t     : Double_t t[y3v_]                                       *
*Entries :  1926502 : Total  Size=   31436333 bytes  File Size  =    2718790 *
*Baskets :      116 : Basket Size=    4371456 bytes  Compression=  11.56     *
*............................................................................*
*Br   25 :sx1e      : sx1e/D                                                 *
*Entries :  1926502 : Total  Size=   15417239 bytes  File Size  =    9161955 *
*Baskets :       54 : Basket Size=    1519616 bytes  Compression=   1.68     *
*............................................................................*
*Br   26 :sx2e      : sx2e/D                                                 *
*Entries :  1926502 : Total  Size=   15417239 bytes  File Size  =    9313360 *
*Baskets :       54 : Basket Size=    1519616 bytes  Compression=   1.66     *
*............................................................................*
*Br   27 :sx3e      : sx3e/D                                                 *
*Entries :  1926502 : Total  Size=   15417239 bytes  File Size  =    8957602 *
*Baskets :       54 : Basket Size=    1519616 bytes  Compression=   1.72     *
*............................................................................*
*Br   28 :sy1e      : sy1e/D                                                 *
*Entries :  1926502 : Total  Size=   15417239 bytes  File Size  =    8916249 *
*Baskets :       54 : Basket Size=    1519616 bytes  Compression=   1.73     *
*............................................................................*
*Br   29 :sy2e      : sy2e/D                                                 *
*Entries :  1926502 : Total  Size=   15417239 bytes  File Size  =    9295908 *
*Baskets :       54 : Basket Size=    1519616 bytes  Compression=   1.66     *
*............................................................................*
*Br   30 :sy3e      : sy3e/D                                                 *
*Entries :  1926502 : Total  Size=   15417239 bytes  File Size  =    8982079 *
*Baskets :       54 : Basket Size=    1519616 bytes  Compression=   1.72     *
*............................................................................*

---

在后续分析中继续保持 vector 类型 Branch

如果已经生成了含 vector Branch 的 ROOT 文件，后续分析时往往希望继续以 vector 的形式读取这些数据，而不是重新退回到数组形式。此时需要特别注意：不能直接沿用 MakeClass 的默认处理方式。

对本例中 split 的 vector<dssd> Branch，MakeClass 会按长度与成员数组生成读取代码。这种方式可以访问数据，但不直接提供整个 vector 的容器操作。需要继续按 vector 读入时，可用下面的对象指针绑定方式。

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 <limits>
include <algorithm>
include <TFile.h>
include <TTree.h>
using namespace std;
struct dssd//aside
{
Int_t id;
Double_t e;
Double_t t = std::numeric_limits<double>::quiet_NaN();
};
struct DSSD//x-y side
{
int xid;
int yid;
double xe; // X 面幅度
double ye; // Y 面幅度
};
class ana
{
public:
vector<dssd> *br_x1v, *br_x2v, *br_x3v;  //声明vector指针
vector<dssd> *br_y1v, *br_y2v, *br_y3v;
Double_t sx1e,sx2e,sx3e;//sum，与输入 tree 的 /D 一致
Double_t sy1e,sy2e,sy3e;
TTree *ipt;
Long64_t source_entry = 0;
TTree *opt;
vector<DSSD> d1,d2,d3; //output
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() 之前，应先将这些指针初始化为空：

void ana::SetBranchInput()
{
  ipt->SetBranchAddress("source_entry", &source_entry);
  br_x1v = nullptr; // ROOT 读入后使指针指向对应的 vector
  br_x2v = nullptr;
  br_x3v = nullptr;
  br_y1v = nullptr;
  br_y2v = nullptr;
  br_y3v = nullptr;
  ipt->SetBranchAddress("x1v", &br_x1v); //将变量指向对应Branch的地址
  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("source_entry", &source_entry, "source_entry/L");
  opt->Branch("d1",&d1);
  opt->Branch("d2",&d2);
  opt->Branch("d3",&d3);
}

排序与配对

由于同一侧探测器可能存在多个 hit，在进行 x-y 关联之前，常常需要先对 hit 按能量进行排序。例如，可以定义一个比较函数：

bool SortDS(const dssd &a, const 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)
{
    xy.clear();
    const size_t nPairs=std::min(x->size(),y->size());
    for(size_t i=0; i<nPairs; ++i) {
        const dssd &xhit=(*x)[i]; // 引用完整 hit，条号和幅度保持对应
        const dssd &yhit=(*y)[i];
        if(std::abs(xhit.e-yhit.e)<50) {
            DSSD pair;
            pair.xid=xhit.id;
            pair.yid=yhit.id;
            pair.xe=xhit.e;
            pair.ye=yhit.e;
            xy.push_back(pair);
        }
    }
}

这里的能量匹配条件、配对策略以及排序方式都可以根据具体实验数据进行调整。

后续分析主循环

在 Analysis() 中，完整流程通常包括：

• 设置输入 Branch；
• 设置输出 Branch；
• 逐事件读取数据；
• 对各层探测器的 x、y hit 分别排序；
• 进行 x-y 配对；
• 将得到的结果写入新的输出树。

void ana::Analysis()
{
  if (ipt == 0) return;
  SetBranchInput();
  BranchOutput();
  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);
    opt->Fill(); // 无候选时保存空 vector，不改变事件顺序

}
}

后续分析时的 main.cpp 与前面类似，只是输入文件变为 vec_16C.root，输出文件变为 sort_16C.root：

#include <TFile.h>
#include <TTree.h>
#include <iostream>
#include "ana.h"
int main(int argc, char** argv) {
    const char* inputName = argc>1 ? argv[1] : "../../vec_16C.root";
    const char* outputName = argc>2 ? argv[2] : "../../sort_16C.root";
    TFile* input = TFile::Open(inputName);
    if (!input || input->IsZombie()) return 1;
    TTree* tin = input->Get<TTree>("tree");
    if (!tin) return 1;
    TFile output(outputName,"RECREATE");
    TTree* tout = new TTree("tree","vector branch");
    {
        ana analysis(tin,tout);
        analysis.Analysis();
        std::cout << "Input=" << tin->GetEntries() << ", output=" << tout->GetEntries() << '\n';
        output.cd();
        tout->Write();
    }
    delete input;
    return 0;
}

如果输出中又引入了新的结构体 DSSD 以及 vector<DSSD>，那么也需要继续在 Linkdef.h 中补充对应的字典声明：

#ifdef __CLING__

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 的主程序读取 strip_arrays_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 明确给出了输入文件 strip_arrays_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 成员函数包括：

x1v.size();             // 当前 hit 数
x1v.clear();            // 清空当前事件的 hit
x1v.push_back(ds);       // 加入一个 dssd 元素
x1v.begin();             // 指向第一个元素的迭代器
x1v.end();               // 指向末尾之后的位置，不能解引用
x1v.assign(x2v.begin(), x2v.end()); // 复制另一个 vector 的内容
// 若 it 指向有效元素，erase 返回删除位置之后的迭代器：
// it = x1v.erase(it);

排序

若要按能量从高到低排序，可结合比较函数与 sort() 使用：

bool SortDS(const dssd &a, const dssd &b)
{
  return a.e > b.e;
}
sort(x1v.begin(), x1v.end(), SortDS);

删除重复元素

如果需要删除重复 hit，可以先定义“重复”的判据，再结合 unique() 与 erase() 使用：

bool Equal(dssd &a, dssd &b)
{
  return a.id == b.id && a.e == b.e && a.t == b.t;
}

void Unique(vector<dssd> &a)
{
a.erase(unique(a.begin(), a.end(), Equal), a.end());
}
// 仅演示删除相邻的等价元素；先独立确认确为重复记录
Unique(xvec);

按条件删除元素

在很多分析中，需要从 vector 中删除不满足条件的元素。例如，当某个 hit 的时间与参考时间相差过大时，可以将其剔除：

void tCut(vector<dssd> &a, vector<dssd> &b, double t1, double t2)
{
  if (b.size() > 0) {
    const double referenceTime = b[0].t; // 删除元素前保存参考时间
    for (auto it = a.begin(); it != a.end(); ) {
      double dt = it->t - referenceTime;
      if (dt < t1 || dt > t2)
        it = a.erase(it);
      else
        ++it;
    }
  }
}

tCut(x1v, x1v, -20, 20);

这些例子说明，使用 vector 的优势不仅在于能够保存变长 hit 信息，还在于可以直接利用 STL 提供的排序、去重和筛选等操作，使后续分析实现起来更加自然。

unique 只去掉相邻的等价元素；相同条号和近似相等的能量不能证明是重复读出，真实 pileup 也可能如此。数据去重应有原始事件标识、timestamp 或电子学重复记录的证据。本节参考文件没有时间，不能运行基于 t 的物理 cut；时间选择代码仅说明有有效时间输入时的容器用法。

按能量排序后同下标配对，是单条响应和清晰能量分离下的初步候选算法，并未覆盖 3.6 的 sharing、共用条及多解事件。保留未匹配结果，再按需要回到完整重建。