ROOT
本讲要点 什么是ROOT 登录ROOT环境和体验中心 ROOT的语法简介 ROOT的函数,直方图,随机数,文件,散点图 TTree的定义、填充、保存、读取及查看 直方图加减乘除运算、归一及拟合
什么是 ROOT ? ROOT: Executive Summary ... provides a set of OO frameworks with all the functionality needed to handle and analyse large amounts of data in a very efficient way.... (摘自http://root.cern.ch/root/Mission.html) 关键字:面向对象的框架、所有功能、海量数据、非常有效
安装ROOT ubuntu 14.04虚拟机已经安装并设置好了ROOT ! 如果没有使用虚拟机,且尚未安装ROOT: 设置ROOT的环境变量(bash) export ROOTSYS=/projects/soft/ext/root export PATH=$ROOTSYS/bin:$PATH export LD_LIBRARY_PATH=$ROOTSYS/lib:$LD_LIBRARY_PATH 可以把上面这3行放到~/.login或者~/.bashrc文件中,这样每次登录到Linux系统,系统就自动设置ROOT的环境变量
登录ROOT环境 运行 >root [-l] 退出 root[0].q 键入 help 指令,如 root[0]? root[1].ls root[2].!ls ROOT环境其它常用指令: .L macro.C Load文件macro.C .x macro.C 执行文件macro.C .ls 显示ROOT当前环境的所有信息 .! ls 显示Linux系统当前目录的所有信息 注:ROOT环境中,ROOT指令都以“.”开头 系统指令都以“.!”开头
ROOT体验中心(1) 在$ROOTSYS/tutorials目录下,有五花八门的例子。 以后会经常与这个目录打交道。先尝试一下吧。 尝试方法: >cd /workdir/root >cp -r $ROOTSYS/tutorials . (注意不要把这个"."漏掉了) >cd tutorials 然后找个感兴趣的目录/文件, 执行ROOT脚本,比如 >cd roofit >root -l rf607_fitresult.C 替换掉这个例子 小技巧提示: 根据关键字"xxxx"从tuotorials的例子中寻找线索 grep -sirn "xxxx" $ROOTSYS/tutorials 比如找随机数用法:grep -sirn "random" $ROOTSYS/tutorials
ROOT体验中心(2) 还可以在ROOT网站上看到一些ROOT图片: http://root.cern.ch/drupal/image 事例产生、探测器模拟、事例重建、数据采集、数据分析 Check this link!
ROOT参考资料 ROOT官方网站http://root.cern.ch ROOT手册(Users Guide) (ver 5.34):http://root.cern.ch/drupal/content/root-users-guide-534 ROOT Primer (ver 5.34): http://root.cern.ch/drupal/content/root-primer-534 ROOT Reference guide (ver 5.34): http://root.cern.ch/root/html534/ClassIndex.html $ROOTSYS/tutorials http://root.cern.ch/root/html/tutorials/ 搜索引擎
下载本讲例子到本地计算机 cd workdir/root tar zxvf rootexamples.tgz cd examples wget hep.tsinghua.edu.cn/~zhuxl/weihai/ppt/rootexamples.tgz tar zxvf rootexamples.tgz cd examples export EXAMPLEDIR=`pwd` 本讲所有例题及练习题相关文件都在$EXAMPLEDIR目录下。
ROOT语法(1)—基本信息 ROOT使用C++语法 一段C++程序可以直接在ROOT环境运行 数据类型重定义 int Int_t float Float_t double Double_t ...... ROOT的类都以T开头 如TFile, TH1F, TTree, ... 详细规定参阅ROOT手册第二章关于Convention和Global Variables部分。 可以直接在ROOT环境中运行macro文件(自动调用cint编译器),也可以在makefile中设置好相关参数用g++编译得到可执行文件运行。 也许需要更新ROOT手册。
ROOT语法(2)—直方图类 ROOT中有众多已经定义好的类可供使用,比如直方图家族
ROOT语法(2)—其它类 其它常用类 数学函数:TF1, TF2, TF3... 图 形:TGraph, TGraphErrors, TGraph2D,... 文 件:TFile 画 布:TCanvas, TPad, ... 随 机 数:TRandom,TRandom1,TRandom2,TRandom3 周期 109 10171 1026 106000 速度(ns/call) 34 242 37 45 比如跟数据结构和分析有关的: TTree, TChain, ... 参见 http://root.cern.ch/root/html534/ClassIndex.html 还有很多全局变量,多数以g开头,如: gRandom, gROOT, gStyle, gPad, gEnv, gFile... 速度与CPU和编译器有关 可以直接查看类的定义得到用法。注释文档非常详尽。
ROOT语法(3)—随机数 使用前可根据需要改变随机数种子和机制 gRandom是指向当前随机数产生子的指针,该产生子默认为TRandom3对象。 http://root.cern.ch/root/html534/TRandom.html (为什么看TRandom?因为TRandom1/2/3都继承自TRandom) gRandom->Binomial(ntot, p): 二项分布 gRandom->BreiWigner(mean, gamma) Breit-Wigner分布 gRandom->Exp(tau) 指数分布 gRandom->Gaus(mean,sigma) 高斯分布 gRandom->Integer(imax) (0,imax-1)随机整数 gRandom->Landau(mean,sigma) Landau分布 gRandom->Poisson(mean) 泊松分布(返回int) gRandom->PoissonD(mean) 泊松分布(返回double) gRandom->Rndm() (0,1]均匀分布 gRandom->Uniform(x1,x2) (x1,x2]均匀分布 .... 使用前可根据需要改变随机数种子和机制
ROOT脚本文件示例(1):Macro文件 $EXAMPLEDIR/ex31.C { cout << "Hello ROOT" << endl; int Num=5; for (int i=0;i<Num;i++) { cout << "i=" << i << endl; } Keep it. 纯粹C++语法,执行的时候只需要在命令提示行: cd $EXAMPLEDIR root -l ex31.C
ROOT中的数学函数 画图时采用 root[1]fun_name.Draw(); 制作一维函数曲线图 制作二维函数曲线图 制作三维函数曲线图 root[0]TF1 *f1 = new TF1("f1","x*sin(x)",-5,5); 制作一维函数曲线图 TF1 *fun_name = new TF1("fun_name","expression", x_low,x_high); root[0]TF2 *f2 = new TF2("f2","x*sin(x)+y*cos(y)", -5,5,-10,10); 制作二维函数曲线图 TF2 *fun_name = new TF2("fun_name","expression", x_low,x_high, y_low,y_high); 画图时采用 root[1]fun_name.Draw(); root[0]TF3 *f3 = new TF2("f3","x*sin(x)+y*cos(y) +z*exp(z)",-5,5,-10,10,-20,20); 制作三维函数曲线图 TF3 *fun_name = new TF3("fun_name","expression", x_low,x_high,y_low,y_high,z_low,z_high);
数学函数的定义方式(1) ROOT中定义数学函数的方式多种多样 以上函数都不含参数,但在数据拟合时,我们往往需要定义含未知参数的函数 利用c++数学表达式 TF1* f1 = new TF1("f1","sin(x)/x",0,10); 利用TMath定义的函数 TF1 *f1 = new TF1("f1","TMath::DiLog(x)",0,10); 利用自定义c++数学函数 Double_t myFun(x) { return x+sqrt(x); } TF1* f1 = new TF1("f1","myFun(x)",0,10); 以上函数都不含参数,但在数据拟合时,我们往往需要定义含未知参数的函数
数学函数的定义方式(2) ROOT中定义含未知参数的数学函数 ROOT已经预定义了几种常用的含参函数 gaus:3个参数 expo:2个参数 f(x)=p0*exp(-0.5*((x-p1)/p2)^2)) expo:2个参数 f(x)=exp(p0+p1*x) polN:N+1个参数 f(x)=p0+p1*x+p2*x^2+... 其中N=0,1,2,...,使用时根据需要用pol0,pol1,pol2... landau:3个参数 朗道分布,没有解析表达式 这些预定义函数可直接使用,比如 histogram->Fit("gaus"); //对直方图进行高斯拟合 TF1 *f1=new TF1("f1","gaus",-5,5);
数学函数的定义方式(3) ROOT中自定义含未知参数的数学函数 利用c++数学表达式 利用c++数学表达式以及ROOT预定义函数 TF1* f1 = new TF1("f1","[0]*sin([1]*x)/x",0,10); 利用c++数学表达式以及ROOT预定义函数 TF1* f1 = new TF1("f1","gaus(0)+[3]*x",0,3); 利用自定义的c++数学函数 Double_t myFun(Double_t *x, Double_t *par) { Double_t xx=x[0]; Double_t f=par[0]*exp(-xx/par[1]); return f; } TF1* f1 = new TF1("f1","myFun",0,10,2); 指定参数数目 定义了含参的TF1对象f1之后,可以设定参数初值,比如 f1->SetParameter(0,value); //为第0个参数设初值为value
ROOT脚本文件示例(2):数学函数定义 $EXAMPLEDIR/ex32.C //a simple ROOT macro, ex32.C //说明ROOT中数学函数的使用,如TF1 void ex32() { //定义函数 TF1 *f1 = new TF1("func1","sin(x)/x",0,10); f1->Draw();//画出函数图像 TF1 *f2 = new TF1("func1",“TMath::Gaus(x,0,1)",0,10); f2->SetLineColor(2);//设置颜色为红色 f2->Draw(“same”);//用参数”same”,把f1,f2画在同一个画布上 } 函数名称 函数表达式 函数区间 运行:在命令提示行下 root -l ex32.C 或在ROOT环境下 .x ex32.C 提示:1)脚本中void函数的名字必须与文件名相同(如ex32) 2)ROOT环境中定义类指针之后,如TF1 *f1,之后 输入“f1->”,然后按一下Tab键,可以自动列出 该类对象的成员函数和成员变量
ROOT脚本文件示例(3): 画布,保存图片 $EXAMPLEDIR/ex33.C //说明ROOT画布的使用,TCanvas,保存图形 void ex33() { //define a function sin(x)/x TF1 *f1 = new TF1("func1","sin(x)/x",0,10); //define a Gaussian function, mean=0, sigma=1 TF1 *f2 = new TF1("func2","Gaus(x,0,1)",-3,3); //定义一个画布, TCanvas TCanvas *myC1 = new TCanvas("myC1","A Canvas",10,10,800,600); //将画布分成两部分 myC1->Divide(2,1); myC1->cd(1); //进入第一部分 f1->Draw(); myC1->cd(2); //进入第二部分 f2->Draw(); myC1->SaveAs(“myex33.gif”); myC1->SaveAs(“myex33.eps”); } 像素坐标 (10,10):左上角 (800,600):右下角 描述 名称 运行:在命令提示行下 > root -l ex33.C 或在ROOT环境下 root[0] .x ex33.C
ROOT中统计直方图 定制一维直方图 定制二维图 定制三维图 填充统计图 绘图:root[0]hist_name.Draw(); TH1F *hist_name = new TH1F(“hist_name”,”hist_title”, num_bins,x_low,x_high); 定制二维图 TH2F *hist_name = new TH2F(“hist_name”,”hist_title”, num_bins_x,x_low,x_high,num_bins_y,y_low,y_high); 定制三维图 TH3F *hist_name = new TH3F(“hist_name”,”hist_title”, num_bins_x,x_low,x_high,num_bins_y,y_low,y_high, num_bins_z,z_low,z_high); 填充统计图 hist_name.Fill(x); hist_name.Fill(x,y); Hist_name.Fill(x,y,z); 绘图:root[0]hist_name.Draw();
ROOT脚本文件示例(4): 直方图,TFile,随机数 $EXAMPLEDIR/ex34.C //说明ROOT直方图、随机数的使用,如TH1F, gRandom void ex34 () { const Int_t NEntry = 10000 ; //创建一个root文件 TFile *file = new TFile(“hist1.root”,”RECREATE”); TH1F *h1 = new TH1F("h1","A simple histo",100,0,1); //填充直方图10000次,用(0,1)均匀分布 for (int i=0;i<NEntry;i++) h1->Fill( gRandom->Uniform() ); h1->Draw(); h1->GetYaxis()->SetRangeUser(0,150); h1->GetXaxis()->SetTitle("x"); h1->GetXaxis()->CenterTitle(); file->cd(); //进入文件file h1->Write();//将h1写入文件 } 名称 描述 No. of Bin 区间 调用均匀分布Uniform(),其它: Landau(mean,sigma); Binomial(ntot,prob); Poisson(mean); Exp(tau); BreitWigner (mean,sigma); 执行的时候只需要在命令提示行 root -l ex34.C 或者进入ROOT环境之后,运行 .x ex34.C
直方图、打开root文件 打开已有的root文件,如hist1.root: 终端提示行下: root -l hist1.root 直方图的Title 直方图统计信息 事例数:Entries 均 值:Mean 方 差:RMS 参见ROOT手册第3章 “Statistics Display” X轴的名称 打开已有的root文件,如hist1.root: 终端提示行下: root -l hist1.root ROOT环境下: TFile f1(“hist1.root”); .ls h1->Draw();
ROOT脚本文件示例(5): 散点图 $EXAMPLEDIR/ex35.C //2维直方图TH2F,散点图,散点图的协方差 void ex35() { const Int_t NEntry = 10000 ; TH2F *hXY = new TH2F("hXY","2d histo",100,0,1,100,-3,3); for (int i=0;i<NEntry;i++) { float x = gRandom->Rndm() ; float y = gRandom->Gaus(0,1) ; hXY->Fill(x,y) ; //填充2维直方图 } hXY->Draw(); //2维直方图的散点图 hXY->GetXaxis()->SetTitle("X: Uniform" ); hXY->GetYaxis()->SetTitle("Y: Gaussian"); Float_t covar = hXY->GetCovariance(); //协方差 cout << "Covariance = " << covar << endl; 二维直方图的Draw()函数 有很多选项,请自行选择 运行:在命令提示行下 root -l ex35.C 或在ROOT环境下 .x ex35.C
ROOT脚本文件示例(5): colz图 $EXAMPLEDIR/ex35.C hXY->Draw(“colz”);
TTree要点 ROOT中tree的概念(类TTree) 什么是tree,为什么tree存取数据 如何定义、填充TTree并写入文件 TChain
root 文件与它的 tree 概念 一个 root 文件就像 UNIX 中的一个目录,它可以 包含目录和没有层次限制的目标模块。 目录中存放不同的类对象或普通数据。 如要求存储大量的同类目标模块,需要引入概念: TTree: 减小磁盘空间和增加读取速度方面被优化 TNtuple:只能存储浮点数的TTree。尽量避免使用。 TTree 采用了 branch 的体系,每个 branch 的读取 可以与别的 branch 无关。
为什么使用TTree 适用于大量的类型相同的对象 可以存储包括类对象、数组等各种类型数据 一般情况下,tree的Branch,Leaf信息就是一个事例的完整信息 有了tree之后,可以很方便对事例进行循环处理。 占用空间少,读取速度快 TTree是ROOT最强大的概念之一
TTree的定义 参见 http://root.cern.ch/root/html534/TTree.html 构造函数: TTree(const char* name, const char* title, Int_t splitlevel = 99); Branch成员函数: virtual TBranch*Branch(const char* name, void* address, const char* leaflist, Int_t bufsize = 32000); 名称 描述 创建TTree,并设置Branch,比如: Int_t RunID; TTree *t1 = new TTree("t1","test tree"); TBranch *br = t1->Branch("RunID",&RunID,"RunID/I"); Branch可以是单独的变量,也可以是一串变量,也可以是定长或不定长数组,也可以是C结构体,或者类对象(继承自TObject,如TH1F对象)。
如何写一个简单的TTree $EXAMPLEDIR/ex41.C void ex41() { TFile *f = new TFile("tree1.root","recreate"); TTree *t1 = new TTree("t1","test tree"); gRandom->SetSeed(0); Float_t px,py,pz; Double_t random; Int_t i; //Set the Branches of tree t1->Branch("px",&px,"px/F"); t1->Branch("py",&py,"py/F"); t1->Branch("pz",&pz,"pz/F"); t1->Branch("random",&random,"random/D"); t1->Branch("i",&i,"i/I"); for (i=0;i<5000;i++) { gRandom->Rannor(px,py); pz = px*px + py*py; random = gRandom->Rndm(); t1->Fill();//Fill tree } t1->Write(); 定义tree,参数分别为tree的名称和描述 设置Branch,参数分别为Branch的“名称”、“地址”以及“leaf列表和类型”。这里只有一个leaf,如果多个则用冒号分开。 常用类型:C,I,F,D分别表示字符串、整型、浮点型和双精度型,参见ROOT手册第12章 为每个leaf赋值,每个事例结束时填充一次。这里一共填充5000事例。 好的做法是实验一个事例填充一次!!! 将tree写入root文件中存盘 运行:root -l ex41.C 或ROOT环境中: .x ex41.C
查看Tree的信息 >root -l tree1.root 打开root文件 root[1].ls 查看文件信息, 发现TTree t1 root[2]t1->Show(0); 显示第0个event的信息 root[3]t1->GetEntries() 总事例数 root[4]t1->Scan(); root[5]t1->Print(); root[6]t1->Draw("px");
查看Tree的信息(续) 也可以 >root -l 进入root root[0]TFile *f1=new TFile("tree1.root"); root[1]t1->Draw( "sqrt(px*px+py*py)" ); root[2]TH1F *h1; root[3]t1->Draw("px>>h1"); root[4]t1->Draw("py","px>0","sames"); root[5]t1->Draw("py","","sames"); root[5]t1->StartViewer()
这是当前目录,双击进入并选择要打开的root文件,以及文件中的tree,最后可以看到tree的各个leaf TBrowser b TBrowser打开一个浏览器,从中可以选择root文件,并一层层进入其中的tree,branch以及leaf。类似于Windows下的Explorer。 这是当前目录,双击进入并选择要打开的root文件,以及文件中的tree,最后可以看到tree的各个leaf 双击leaf可以查看leaf的直方图 适合初步浏览,但不适合具体的数据分析处理。并不推荐使用。
如何读写含有不定长数组的tree(1) $EXAMPLEDIR/ex42.C ... const Int_t kMaxTrack = 50; Int_t ntrack; Float_t px[kMaxTrack]; Float_t py[kMaxTrack]; Float_t zv[kMaxTrack]; Double_t pv[3]; TFile f(rootfile,"recreate"); TTree *t3 = new TTree("t3", "Reconst events"); t3->Branch("ntrack",&ntrack,"ntrack/I"); t3->Branch("px",px,"px[ntrack]/F"); t3->Branch("py",py,"py[ntrack]/F"); t3->Branch("zv",zv,"zv[ntrack]/F"); t3->Branch("pv",pv,"pv[3]/D"); void ex42r() {//读取数据,适用于简单分析 TFile *f = new TFile(rootfile); TTree *t3 = (TTree*)f->Get("t3"); t3->Draw("sqrt(px*px+py*py)"); htemp->SetLineColor(2); t3->Draw("sqrt(px*px+py*py)", "zv>100","sames"); } !!如何获取root文件中的tree指针 直接画出Branch/Leaf, 可以加很多条件。 1)估计不定长数组的最大维数,以该维数定义数组;如float zv[kMaxTrack] 2)定义某变量,用于存放数组的实际维数。如int ntrack,表示一个事例中实际的径迹数。 3)定义tree,设置Branch。第三个参数给出数组的实际维数。如”zv[ntrack]/F” 运行:进入ROOT环境后 .L ex42.C ex42w() ex42r() 很多时候不定长数组是必要的,比如正负电子对撞,记录末态粒子的信息,末态粒子数目是不固定的。
如何读写含有不定长数组的tree(2) $EXAMPLEDIR/ex42.C void ex42r2() { TFile *f = new TFile(rootfile); TTree *t3 = (TTree*)f->Get("t3"); //步骤1:定义好必要的变量 const Int_t kMaxTrack = 100; Int_t ntrack; Float_t px[kMaxTrack]; //[ntrack] Float_t py[kMaxTrack]; //[ntrack] Float_t zv[kMaxTrack]; //[ntrack] Double_t pv[3]; //步骤2: 用SetBranchAddress函数 //将tree的Branch与定义好的变量 //地址联系起来。 t3->SetBranchAddress("ntrack", &ntrack); t3->SetBranchAddress("px", px); t3->SetBranchAddress("py", py); t3->SetBranchAddress("zv", zv); t3->SetBranchAddress("pv", pv); //获取事例总数 Int_t nentries = t3->GetEntries(); TH1I *hntrack = new TH1I("hntrack",“trk n",25,0,50); TH1F *hpt = new TH1F("hpt" ,“trk pt" ,100,0,10); //步骤3:对所有事例循环 for (int i=0;i<nentries;i++) { t3->GetEntry(i); //获取第i个事例 hntrack->Fill( ntrack ); for (int j=0;j<ntrack;j++) { Float_t pt = sqrt(px[j]*px[j]+py[j]*py[j]); hpt->Fill( pt ); } TCanvas *myC = new TCanvas("myC","",10,10,600,400); hntrack->Draw("e"); hntrack->GetYAxis()->SetRangeUser(0,60); TCanvas *myC1 = new TCanvas("myC1","",10,10,600,400); hpt->Draw(); 问题:对ntrack和px的处理有什么差别?为什么? 运行: 进入ROOT环境后 .L ex42.C ex42w() ex42r2() 每获取一个事例,这些Branch直接赋值给指定的变量。可以在循环中设定选取条件,选择分析数据。进行复杂细致的分析推荐使用这种方法。 注:程序中//[ntrack]的意义参见ROOT手册Streamer
如何将类对象设定为tree的Branch $EXAMPLEDIR/ex43.C ... TTree *t3 = new TTree("t3","Reconst events"); //Evt_t是已经定义好的类(详见ex43.C)。这里的myevt一定要用new的方式定义, //然后就可以直接将该对象设为tree的一个Branch。 //第1个参数为Branch的名字,第2个为类的名字(可省略),第3个为对象指针的地址(!!), //第4个为缓存大小,第5个为分割级别(split-level)。 //参见手册“Adding a Branch to Hold an Object” Evt_t *myevt = new Evt_t(); t3->Branch(“evt”,“Evt_t”,&myevt,32000,1); //读取tree时,可以直接将对象指针的地址的赋给相应的Branch //需要提醒的是,第1个参数为Branch的名称,必须与写入时指定的名称相同。 Evt_t *myevt = 0 ; t3->SetBranchAddress("evt",&myevt); //GetEntry(i)获得第i个entry后,通过myevt->ntrack(或其它成员变量)可以获得 //相应的数据。 t3->GetEntry(i); hntrack->Fill( myevt->ntrack ); 略过,但并非不重要 语法更严格,必须include需要的头文件 详见ex43.C。运行时必须通过外部编译器如:root -l ex43.C+ 或者在ROOT环境中 .x ex43.C+ 这里的”+”是必须的。
如何从ASCII文件中读取数据转为ROOT中的TTree /home/yangzw/examples/Lec4/ex44.C 有时候记录的实验数据是ASCII格式,或者二进制格式。我们可以读取这些数据转换成ROOT中的tree,方便进行数据分析。ex44.C读取basic.dat(ASCII码),转成最简单的TTree。basic.dat中的数据分3列,分别为x,y,z坐标。 void ex44() { ifstream in; // 定义文件流对象, i表示in,f表示file,即从某文件中读取数据 in.open(“basic.dat”); //打开该文件 Float_t x,y,z; Int_t nlines = 0; TFile *f = new TFile("ex44.root","RECREATE"); TH1F *h1 = new TH1F("h1","x distribution",100,-4,4); //TNtuple可以看成特殊的的TTree,只可以存放浮点型数据。 TNtuple *ntuple = new TNtuple("ntuple","data aus ascii","x:y:z"); while (1) { in >> x >> y >> z; //从文件中读取一行,分别赋值给x,y,z。 if (!in.good()) break; if (nlines < 5) printf("x=%8f, y=%8f, z=%8f\n",x,y,z); h1->Fill(x); ntuple->Fill(x,y,z); //填充TNtuple nlines++; } printf(" found %d points\n",nlines); in.close(); f->Write(); 注:3个Branch 分别为x,y,z 每个事例填充一次 运行:root -l ex44.C 或者在root环境中: .x ex44.C
改成TTree方式: mytree->Branch(“y”,&y,”y/F”); void ex44() { Float_t x,y,z; Int_t nlines = 0; TFile *f = new TFile("ex44.root","RECREATE"); TH1F *h1 = new TH1F("h1","x distribution",100,-4,4); //TNtuple可以看成特殊的的TTree,只可以存放浮点型数据。 TTree *mytree = new TTree(“mytree”,”aaaaaaaaaaaa”); mytree->Branch(“x”,&x,”x/F”); mytree->Branch(“y”,&y,”y/F”); mytree->Branch(“z”,&z,”z/F”); mytree->Branch(“w”,&w,”w/C”); while (1) { in >> x >> y >> z; //从文件中读取一行,分别赋值给x,y,z。 if (!in.good()) break; if (nlines < 5) printf("x=%8f, y=%8f, z=%8f\n",x,y,z); h1->Fill(x); mytree->Fill();//填充TNtuple nlines++; } printf(" found %d points\n",nlines); in.close(); f->Write(); 略过,但并非不重要
如何从ASCII文件中读取数据转为ROOT中的TTree $EXAMPLEDIR/ex44a.C 读取ASCII格式文件还有个更简便的方式,即用TTree的ReadFile函数: tree->ReadFile(parameter1,parameter2); 参见ROOT手册TTree那一章 void ex44a() { TFile *f = new TFile("ex44.root","RECREATE"); TTree *T = new TTree("ntuple","data from ascii file"); //第1个参数为要打开的文件名称 //第2个参数是Branch的描述,即设定3个Branch x,y,z Long64_t nlines = T->ReadFile("basic.dat","x:y:z"); printf(" found %lld points\n",nlines); T->Write(); } 运行:root -l ex44a.C 或者在root环境中: .x ex44a.C
TChain: 分析多个root文件的利器(1) TChain对象是包含相同tree的ROOT文件的列表。 参见User’s Guide中TChain相关章节以及 http://root.cern.ch/root/html534/TChain.html void ex45() { //定义TChain,t3为root文件中tree的名称!!!!!!! TChain* fChain= new TChain(“t3”); //添加所有文件至fChain,或根据需要添加部分root文件 fChain->Add(“rootfiles/*.root”); //画出t3的某个leaf,如ntrack fChain->Draw(“ntrack”); } 注意:此时,fChain等同于一个大root文件中的一个类"t3",该文件包含的事例数为所有文件中事例数之和(1012以内) 问题:root文件中多个tree怎么办?
直方图的操作 直方图的运算 加减乘除: Add,Divide,... 归一化:Scale ROOT中直方图拟合 h1->Fit();
直方图归一化(1) 直方图的归一化 void TH1::Scale(Double_t c1, Option_t *option) http://root.cern.ch/root/html534/TH1.html#TH1:Scale 直方图的归一化 void TH1::Scale(Double_t c1, Option_t *option) 默认c1=1,把直方图每个区间的值(BinContent)乘以c1 假设 sum=h1->Integral() h1->Scale(c1)之后, h1->Integral() = c1*sum 不加参数时,h1->Scale() 没有变化(默认c1=1) "归一化"后,不仅BinContent变化了,BinError也变化了
直方图的归一化(2) 归一化常用于比较两种分布,找出区别。 所以,将2个直方图归一化到积分相同进行比较才直观。 root[0]TH1F *h1=new TH1F("h1","",100,-5,5); root[1]TH1F *h2=new TH1F("h2","",100,-5,5); root[2]h1->FillRandom("gaus",5000); root[3]h2->FillRandom("gaus",10000); root[4]float norm=1000; root[5]h1->Scale(norm/h1->Integral()); root[6]h2->Scale(norm/h2->Integral()); root[7]h1->Draw("e"); root[8]h2->Draw(“esames”) ; 注意Draw()函数的选项 "归一化"之后,h1或h2->Integral()=norm 在同一张图上可以看出比较2个分布的差别。
直方图四则运算(1) 重要提示: 对直方图进行四则运算操作,一定要想明白运算的意义 比如两个直方图的相加与两个随机变量的卷积有什么区别 2. 两个直方图的四则运算,区间大小和区间数相同才有意义 四则运算"加减乘除"分别对应 统计量(BinContent)的相加、相减、相乘、相除 3. 如果需要正确处理统计误差,需要在对ROOT脚本中 调用TH1的某个静态成员函数,即 TH1::SetDefaultSumw2(); void SetDefaultSumw2(Bool_t sumw2 = kTRUE) //static function. When this static function is called with sumw2=kTRUE, all new histograms will automatically activate the storage of the sum of squares of errors, ie TH1::Sumw2 is automatically called.
直方图的四则运算(2) 相加:常用于相同实验的数据叠加,增加统计量。 ...... root[1]TH1::SetDefaultSumw2(); root[1]TH1F *h3=new TH1F(*h1); root[2]h3->Add(h1,h2,a,b); 结果:h3的BinContent被a*h1+b*h2替换,一般a=b=1 相加:常用于相同实验的数据叠加,增加统计量。 ...... root[1]TH1F *h3=new TH1F(*h1); root[2]h3->Sumw2();//也可在定义h3前TH1::SetDefaultSumw2(); root[3]h3->Add(h1,h2,a,-b); 结果:h3的BinContent被a*h1+b*h2替换,一般a=-b=1 相减:常用于从实验测量的分布中扣除本底。
直方图的四则运算(3) 相除 root[1]TH1F *h3=new TH1F(*h1); root[2]h3->Sumw2(); 常用于效率的计算。 root[1]TH1F *h3=new TH1F(*h1); root[2]h3->Sumw2(); root[3]h3->Divide(h1,h2,a,b); root[4]h3->Divide(h1,h2,a,b); 思考:如果h1和h2不独立,怎么办?比如h1包含于h2 root[4]h3->Divide(h1,h2,a,b,"B"); 相乘 常用于对分布进行诸如效率等的修正。 root>TH1F *h3=new TH1F(*h1); root>h3->Sumw2(); root>h3->Multiply(h1,h2,a,b);
直方图四则运算的误差处理 虽然ROOT都提供了较完善的一维直方图运算功能,但对最终结果的误差一定要仔细检查。很多情况下,用户需要从图中读出各频数数值与误差值,并确认运算无误。 包括归一化和加减乘除在内, 如果希望使用直方图的误差,都需要调用 TH1::SetDefaultSumw2(); 或者,对每个直方图(如hist)调用 hist->Sumw2();
拟合直方图(1) 将鼠标放到直方图上,右键,出现直方图操作选项,选择FitPanel,可以在FitPanel中选择拟合的各个选项,比如用什么函数拟合,拟合的区间,等等。
拟合直方图(2) 并不推荐这种拟合方式: 1)不适合自定义函数拟合 2)不适合批处理 用默认的高斯拟合,并在Options菜单中选上Fit Parameters选项,可以看到拟合的结果。 拟合结果给出了高斯分布的3个参数: 常系数、均值、均方差,以及拟合的好坏chi2/ndf 并不推荐这种拟合方式: 1)不适合自定义函数拟合 2)不适合批处理
拟合直方图(3) $EXAMPLEDIR/ex51.C hpx->Fit("gaus"); hpx->Fit("gaus","","",-3,3); 自定义拟合函数 TF1 *fcn = new TF1("fcn","gaus",-3,3); hpx->Fit(fcn,”R”); gStyle->SetOptFit();//设置拟合选项 拟合前往往需要给出合理的参数初值 fcn->SetParameters(500,mean,sigma); 拟合后取出拟合得到的参数 Double_t mypar[3]; fcn->GetParameters(&mypar[0]); 运行:root -l root [0] .L ex51.C root [1] ex51r() root [2] ex51r2() 用自定义的函数拟合直方图
拟合直方图(3) $EXAMPLEDIR/ex52.C 共振峰(Breit-Wigner分布)加上二次函数本底的拟合(一共6个参数) 先自定义本底函数(background)和共振峰函数(lorentianPeak),再定义这两个函数的和为拟合函数:fitFunction 利用fitFunction定义TF1函数 TF1 *fitFcn = new TF1("fitFcn",fitFunction,0,3,6); 这里指定函数区间为0-3,6个参数 fitFcn->SetParameter(4, 0.2); 为某个参数设初值(width) fitFcn->SetParLimits(5, 0.6,1.4); 为某参数设置取值范围 运行:root -l root [0] .L ex52.C 注意TLegend的使用
ROOT小结 设定ROOT环境变量: ROOTSYS,PATH,LD_LIBRARY_PATH 绘制各种直方图,散点图,数学函数 TH1F,TH2F,TF1,... 随机数产生子,各种分布 gRandom->Rndm,Uniform,Gaus,Exp,... 创建、保存root文件 TFile *f = new TFile("myfile.root”,”recreate"); f->Write(); TTree, TChain的使用 TTree *mytree = new TTree("mytree”,”my tree"); mytree->Branch(........); 用TChain分析相同格式的数据文件。 直方图的运算,拟合 h1->Fit("function_name");
ROOT的重要功能或用法(1) 带参数的macro;ACLIC模式,参考ROOT手册 ROOT手册13、14章分别是数学库和线性代数,提供很多数学功能,比如Lorentz矢量的操作,特殊函数,矩阵求解运算,求极值等等 ROOT手册第4章介绍Graphs,适用于不等距数据的图形分析(当然也可以构造不等bin的直方图) RooFit,最大似然法拟合等 神经网络分析方法,TMVA(多元数据分析) ROOT中使用PYTHIA、Geant3/4 图形接口...
ROOT的重要功能或用法(2) MakeClass,MakeSelector的运用 比如$EXAMPLEDIR目录下有文件ex51.root,其中含有复杂的tree。可以用MakeClass或MakeSelector自动产生分析文件和头文件: root [0] TFile f("ex51.root"); root [1] .ls TFile** ex51.root TFile* ex51.root KEY: TTree t4;1 Reconst events root [2] t4->MakeClass(); 或: t4->MakeClass(“MyClass”); 自动产生以t4.h和t4.C文件, 或MyClass.h和MyClass.C文件。 类的定义以及Branch地址设定、分析框架都已经自动完成。 MakeSelector的用法类似: root [0] TFile f("ex51.root"); root [1] t4->MakeSelector(); 或 t4->MakeSelector(“MySelector”);
直方图画图练习 1. 写一个ROOT脚本,ex3_gaus.C, 调用随机数产生子产生高斯分布,区间(-6,6),分30个bin,画出直方图,比较不同的参数的分布。 参数组合为:(mean,sigma)=(0,1), (0,2), (1,1), (1,2), 把这4个分布的直方图画在同一个图中进行比较。 hint:高斯分布用gRandom->Gaus(mean,sigma)产生。 使用Draw()函数的"same"参数可以在一个画板上画多个图。 2. 写一个ROOT脚本,ex3_pdf.C,作4个直方图,分别产生10000事例的Gauss,Poisson,Binomial,Landau分布。创建画布,分成2*2块,将4个直方图画在画布的1-4部分。注意不同分布的参数选择合理性,比如Binomial(ntot,p), ntot>0, 0<p<1. 定义一个二维直方图(TH2F), 将随机产生的1000个坐标(x,y)填充到直方图中,其中x和y都是(0,1)之间的均匀分布。画出散点图,查看x和y的关联。用hint:用gRandom->Rndm()产生均匀分布。 3. 将练习2中产生的直方图储存到mypdf.root文件中。 将所画直方图的x/y轴添加上名称,不同分布用不同颜色。 将画布存成eps文件和gif文件 4. cp –r $ROOTSYS/tutorials ~/workdir/root 运行tutorials/hist/以下几个文件,查看ROOT直方图的常用功能如何实现 twoscales.C, transpad.C, multicolor.C, logscales.C, hstack.C
TTree练习 1. $EXAMPLEDIR/ascii_data/random1.dat 该文本文件有5列数据,前2列为整型,之后两列为浮点型,最后一列为字符串型。读取该文件,把这些数据写入tree,并存到random1.root文件中。为tree设置5个Branch:npart, ntrack,px,py,ch,分别为整型、整型、浮点型、浮点型、字符串 提示:参照ex41.C和ex44.C,将二者综合起来。 设置字符串的Branch可以如下: Char_t ch[4]; tree1->Branch(“ch”,ch,”ch/C”); 2. 用练习1的程序将$EXAMPLEDIR/ascii_data/random2.dat, random3.dat分别读取为random2.root, random3.root。加上习题1的root文件,共3个。 用这3个root文件数据,画出npart,ntrack,px以及py的分布。 画出npart:ntrack的散点图。 提示:参照ex45.C,用TChain将这3个root文件连起来。 tree->Draw(“npart:ntrack”);可以看2个变量的关联
直方图运算和拟合练习 练习需要的root文件都存放在下面目录里: $EXAMPLEDIR/exercise/ 1. 查看该double_gaus.root文件。其中存储了名为tree1的TTree。画出tree1中pz的分布,并对该分布进行拟合,在图上显示出拟合的结果,并在屏幕上打印出拟合结果。 (提示:该分布为两个高斯的叠加,可以自定义一个包含6个参数的TF1进行拟合,分布比较复杂的时候,需要先估计参数的大概值,为拟合函数预设估计值。) 思考:假设函数fun=p0*exp(-(x-p1)^2/2/p2^2) +p3*exp(-(x-p4)^2/2/p5^2) 由拟合得到的结果,比较两个高斯的份额 2. hist.root中有两个直方图,对这两个直方图进行加减乘除运算。除法时,查看用”B”选项和不用“B”选项时误差的不同。 (提示:h1的事例包含于h2的事例,计算误差需要用”B”选项) 3. 利用1.root和2.root,将其中的px分别画到两个直方图h1,h2中。对h1,h2进行加减乘除运算,查看误差情况。比较调用与不调用Sumw2()的差别。