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.cpp
  • makefile
  • ana.h
  • ana.cpp
  • Linkdef.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 命令行中进行查看和分析。

例如:

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

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

例如:

In [4]:
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 的长度。例如:

In [8]:
tree->Draw("@x1v.size()");//显示vector的大小
c1->Draw();

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

访问 vector 中单个元素的成员¶

可以通过

vec[i].member

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

In [11]:
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

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

展开查看全部 hit 信息¶

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

In [14]:
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 生成的代码会将其解释为长度变量和若干数组成员的组合。例如,执行:

In [9]:
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 提供的排序、去重和筛选等操作,使后续分析实现起来更加自然。

In [ ]: