粒子物理与核物理实验中的数据分析课件.ppt

ROOT,2,本讲要点,什么是ROOT登录ROOT环境和体验中心ROOT的语法简介ROOT的函数,直方图,随机数,文件,散点图TTree的定义、填充、保存、读取及查看直方图加减乘除运算、归一及拟合,3,什么是ROOT？,ROOT:ExecutiveSvidesasetofOOframeworkswithallthefunctionalityneededtohandleandanalyselargeamountsofdatainaveryefficientway.(摘自http:/root.cern.ch/root/Mission.html)关键字：面向对象的框架、所有功能、海量数据、非常有效,4,安装ROOT,ubuntu14.04虚拟机已经安装并设置好了ROOT!如果没有使用虚拟机，且尚未安装ROOT：首先查看ROOT网站http:/root.cern.ch上是否有相应操作系统（或Linux发行版）及编译器的预编译好的程序包，有的话，直接下载解压缩到你想要的安装目录下，比如/projects/soft/ext/root如果没有适合你的操作系统的预编译包，就需要到官网下载ROOT的源代码，按照安装指南编译安装设置ROOT的环境变量（bash）exportROOTSYS=/projects/soft/ext/rootexportPATH=$ROOTSYS/bin:$PATHexportLD_LIBRARY_PATH=$ROOTSYS/lib:$LD_LIBRARY_PATH可以把上面这3行放到/.login或者/.bashrc文件中，这样每次登录到Linux系统，系统就自动设置ROOT的环境变量,5,登录ROOT环境,运行root-l退出root0.q键入help指令，如root0?root1.lsroot2.!ls,ROOT环境其它常用指令：.Lmacro.CLoad文件macro.C.xmacro.C执行文件macro.C.ls显示ROOT当前环境的所有信息.!ls显示Linux系统当前目录的所有信息注：ROOT环境中，ROOT指令都以“.”开头系统指令都以“.!”开头,6,ROOT体验中心(1),在$ROOTSYS/tutorials目录下，有五花八门的例子。以后会经常与这个目录打交道。先尝试一下吧。尝试方法：cd/workdir/rootcp-r$ROOTSYS/tutorials.(注意不要把这个.漏掉了)cdtutorials然后找个感兴趣的目录/文件，执行ROOT脚本，比如cdroofitroot-lrf607_fitresult.C,小技巧提示：根据关键字xxxx从tuotorials的例子中寻找线索grep-sirnxxxx$ROOTSYS/tutorials比如找随机数用法：grep-sirnrandom$ROOTSYS/tutorials,7,ROOT体验中心(2),还可以在ROOT网站上看到一些ROOT图片：http:/root.cern.ch/drupal/image当然，ROOT的功能不只是做图，它不是一个作图工具。跟数据分析有关的东西，基本都是ROOT的擅长；跟物理有关的很多东西，ROOT基本都可以做得很好：事例产生、探测器模拟、事例重建、数据采集、数据分析,ROOT参考资料,8,ROOT官方网站http:/root.cern.chROOT手册（UsersGuide）(ver5.34)：http:/root.cern.ch/drupal/content/root-users-guide-534ROOTPrimer(ver5.34):http:/root.cern.ch/drupal/content/root-primer-534ROOTReferenceguide(ver5.34):http:/root.cern.ch/root/html534/ClassIndex.html$ROOTSYS/tutorialshttp:/root.cern.ch/root/html/tutorials/搜索引擎,9,下载本讲例子到本地计算机,cdworkdir/,10,ROOT语法(1)基本信息,ROOT使用C+语法一段C+程序可以直接在ROOT环境运行数据类型重定义intInt_tfloatFloat_tdoubleDouble_t.ROOT的类都以T开头如TFile,TH1F,TTree,.详细规定参阅ROOT手册第二章关于Convention和GlobalVariables部分。可以直接在ROOT环境中运行macro文件(自动调用cint编译器)，也可以在makefile中设置好相关参数用g+编译得到可执行文件运行。,11,ROOT语法(2)直方图类,ROOT中有众多已经定义好的类可供使用，比如直方图家族,12,ROOT语法(2)其它类,其它常用类数学函数：TF1,TF2,TF3.图形：TGraph,TGraphErrors,TGraph2D,.文件：TFile画布：TCanvas,TPad,.随机数：TRandom,TRandom1,TRandom2,TRandom3周期109101711026106000速度(ns/call)342423745比如跟数据结构和分析有关的：TTree,TChain,.参见http:/root.cern.ch/root/html534/ClassIndex.html还有很多全局变量，多数以g开头，如：gRandom,gROOT,gStyle,gPad,gEnv,gFile.,速度与CPU和编译器有关,13,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均匀分布.,使用前可根据需要改变随机数种子和机制,14,ROOT脚本文件示例(1)：Macro文件$EXAMPLEDIR/ex31.C,用花括号括起来，后缀名一般用”.C”coutHelloROOTendl;intNum=5;for(inti=0;iNum;i+)couti=iSetParameter(0,value);/为第0个参数设初值为value,19,ROOT脚本文件示例(2):数学函数定义$EXAMPLEDIR/ex32.C,/asimpleROOTmacro,ex32.C/说明ROOT中数学函数的使用，如TF1voidex32()/定义函数TF1*f1=newTF1(func1,sin(x)/x,0,10);f1-Draw();/画出函数图像TF1*f2=newTF1(func1,“TMath:Gaus(x,0,1),0,10);f2-SetLineColor(2);/设置颜色为红色f2-Draw(“same”);/用参数”same”，把f1,f2画在同一个画布上,函数名称,函数表达式,函数区间,提示：1)脚本中void函数的名字必须与文件名相同(如ex32)2)ROOT环境中定义类指针之后，如TF1*f1，之后输入“f1-”，然后按一下Tab键，可以自动列出该类对象的成员函数和成员变量,运行：在命令提示行下root-lex32.C或在ROOT环境下.xex32.C,20,ROOT脚本文件示例(3):画布，保存图片$EXAMPLEDIR/ex33.C,/说明ROOT画布的使用，TCanvas，保存图形voidex33()/defineafunctionsin(x)/xTF1*f1=newTF1(func1,sin(x)/x,0,10);/defineaGaussianfunction,mean=0,sigma=1TF1*f2=newTF1(func2,Gaus(x,0,1),-3,3);/定义一个画布,TCanvasTCanvas*myC1=newTCanvas(myC1,ACanvas,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”);,运行：在命令提示行下root-lex33.C或在ROOT环境下root0.xex33.C,名称,描述,像素坐标(10,10):左上角(800,600):右下角,21,ROOT中统计直方图,绘图：root0hist_name.Draw();,22,ROOT脚本文件示例(4):直方图，TFile，随机数$EXAMPLEDIR/ex34.C,/说明ROOT直方图、随机数的使用，如TH1F,gRandomvoidex34()constInt_tNEntry=10000;/创建一个root文件TFile*file=newTFile(“hist1.root”,”RECREATE”);TH1F*h1=newTH1F(h1,Asimplehisto,100,0,1);/填充直方图10000次，用(0,1)均匀分布for(inti=0;iFill(gRandom-Uniform();h1-Draw();h1-GetYaxis()-SetRangeUser(0,150);h1-GetXaxis()-SetTitle(x);h1-GetXaxis()-CenterTitle();file-cd();/进入文件fileh1-Write();/将h1写入文件,执行的时候只需要在命令提示行root-lex34.C或者进入ROOT环境之后，运行.xex34.C,名称,调用均匀分布Uniform()，其它：Landau(mean,sigma);Binomial(ntot,prob);Poisson(mean);Exp(tau);BreitWigner(mean,sigma);,描述,No.ofBin,区间,23,X轴的名称,直方图统计信息事例数：Entries均值：Mean方差：RMS参见ROOT手册第3章“StatisticsDisplay”,直方图的Title,打开已有的root文件，如hist1.root:终端提示行下:root-lhist1.rootROOT环境下：TFilef1(“hist1.root”);.lsh1-Draw();,直方图、打开root文件,24,/2维直方图TH2F，散点图，散点图的协方差voidex35()constInt_tNEntry=10000;TH2F*hXY=newTH2F(hXY,2dhisto,100,0,1,100,-3,3);for(inti=0;iRndm();floaty=gRandom-Gaus(0,1);hXY-Fill(x,y);/填充2维直方图hXY-Draw();/2维直方图的散点图hXY-GetXaxis()-SetTitle(X:Uniform);hXY-GetYaxis()-SetTitle(Y:Gaussian);Float_tcovar=hXY-GetCovariance();/协方差coutCovariance=covarSetSeed(0);Float_tpx,py,pz;Double_trandom;Int_ti;/SettheBranchesoftreet1-Branch(px,如何写一个简单的TTree$EXAMPLEDIR/ex41.C,运行：root-lex41.C或ROOT环境中:.xex41.C,设置Branch，参数分别为Branch的“名称”、“地址”以及“leaf列表和类型”。这里只有一个leaf，如果多个则用冒号分开。常用类型:C,I,F,D分别表示字符串、整型、浮点型和双精度型，参见ROOT手册第12章,为每个leaf赋值，每个事例结束时填充一次。这里一共填充5000事例。好的做法是实验一个事例填充一次！,将tree写入root文件中存盘,31,查看Tree的信息,root-ltree1.root打开root文件root1.ls查看文件信息，发现TTreet1root2t1-Show(0);显示第0个event的信息root3t1-GetEntries()总事例数root4t1-Scan();root5t1-Print();root6t1-Draw(px);,32,查看Tree的信息(续),也可以root-l进入rootroot0TFile*f1=newTFile(tree1.root);root1t1-Draw(sqrt(px*px+py*py);root2TH1F*h1;root3t1-Draw(pxh1);root4t1-Draw(py,px0,sames);root5t1-Draw(py,sames);root5t1-StartViewer(),33,查看Tree的信息(2),TBrowserbTBrowser打开一个浏览器，从中可以选择root文件，并一层层进入其中的tree，branch以及leaf。类似于Windows下的Explorer。,这是当前目录，双击进入并选择要打开的root文件，以及文件中的tree，最后可以看到tree的各个leaf,双击leaf可以查看leaf的直方图,适合初步浏览，但不适合具体的数据分析处理。并不推荐使用。,34,.constInt_tkMaxTrack=50;Int_tntrack;Float_tpxkMaxTrack;Float_tpykMaxTrack;Float_tzvkMaxTrack;Double_tpv3;TFilef(rootfile,recreate);TTree*t3=newTTree(t3,Reconstevents);t3-Branch(ntrack,.,运行:进入ROOT环境后.Lex42.Cex42w()ex42r(),voidex42r()/读取数据，适用于简单分析TFile*f=newTFile(rootfile);TTree*t3=(TTree*)f-Get(t3);t3-Draw(sqrt(px*px+py*py);htemp-SetLineColor(2);t3-Draw(sqrt(px*px+py*py),zv100,sames);,如何读写含有不定长数组的tree(1)$EXAMPLEDIR/ex42.C,1)估计不定长数组的最大维数，以该维数定义数组；如floatzvkMaxTrack2)定义某变量，用于存放数组的实际维数。如intntrack,表示一个事例中实际的径迹数。3)定义tree，设置Branch。第三个参数给出数组的实际维数。如”zvntrack/F”,很多时候不定长数组是必要的，比如正负电子对撞，记录末态粒子的信息，末态粒子数目是不固定的。,35,voidex42r2()TFile*f=newTFile(rootfile);TTree*t3=(TTree*)f-Get(t3);/步骤1:定义好必要的变量constInt_tkMaxTrack=100;Int_tntrack;Float_tpxkMaxTrack;/ntrackFloat_tpykMaxTrack;/ntrackFloat_tzvkMaxTrack;/ntrackDouble_tpv3;/步骤2:用SetBranchAddress函数/将tree的Branch与定义好的变量/地址联系起来。t3-SetBranchAddress(ntrack,运行:进入ROOT环境后.Lex42.Cex42w()ex42r2(),如何读写含有不定长数组的tree(2)$EXAMPLEDIR/ex42.C,每获取一个事例，这些Branch直接赋值给指定的变量。可以在循环中设定选取条件，选择分析数据。进行复杂细致的分析推荐使用这种方法。注：程序中/ntrack的意义参见ROOT手册Streamer,/获取事例总数Int_tnentries=t3-GetEntries();TH1I*hntrack=newTH1I(hntrack,“trkn,25,0,50);TH1F*hpt=newTH1F(hpt,“trkpt,100,0,10);/步骤3:对所有事例循环for(inti=0;iGetEntry(i);/获取第i个事例hntrack-Fill(ntrack);for(intj=0;jFill(pt);TCanvas*myC=newTCanvas(myC,10,10,600,400);hntrack-Draw(e);hntrack-GetYAxis()-SetRangeUser(0,60);TCanvas*myC1=newTCanvas(myC1,10,10,600,400);hpt-Draw();,36,.TTree*t3=newTTree(t3,Reconstevents);/Evt_t是已经定义好的类(详见ex43.C)。这里的myevt一定要用new的方式定义，/然后就可以直接将该对象设为tree的一个Branch。/第1个参数为Branch的名字，第2个为类的名字(可省略)，第3个为对象指针的地址(!)，/第4个为缓存大小，第5个为分割级别(split-level)。/参见手册“AddingaBranchtoHoldanObject”Evt_t*myevt=newEvt_t();t3-Branch(“evt”,“Evt_t”,如何将类对象设定为tree的Branch$EXAMPLEDIR/ex43.C,详见ex43.C。运行时必须通过外部编译器如：root-lex43.C+或者在ROOT环境中.xex43.C+这里的”+”是必须的。,语法更严格，必须include需要的头文件,略过，但并非不重要,37,有时候记录的实验数据是ASCII格式，或者二进制格式。我们可以读取这些数据转换成ROOT中的tree，方便进行数据分析。ex44.C读取basic.dat(ASCII码)，转成最简单的TTree。basic.dat中的数据分3列，分别为x,y,z坐标。,如何从ASCII文件中读取数据转为ROOT中的TTree/home/yangzw/examples/Lec4/ex44.C,运行：root-lex44.C或者在root环境中:.xex44.C,voidex44()ifstreamin;/定义文件流对象,i表示in，f表示file，即从某文件中读取数据in.open(“basic.dat”);/打开该文件Float_tx,y,z;Int_tnlines=0;TFile*f=newTFile(ex44.root,RECREATE);TH1F*h1=newTH1F(h1,xdistribution,100,-4,4);/TNtuple可以看成特殊的的TTree，只可以存放浮点型数据。TNtuple*ntuple=newTNtuple(ntuple,dataausascii,x:y:z);while(1)inxyz;/从文件中读取一行，分别赋值给x,y,z。if(!in.good()break;if(nlinesFill(x);ntuple-Fill(x,y,z);/填充TNtuplenlines+;printf(found%dpointsn,nlines);in.close();f-Write();,38,改成TTree方式：,voidex44()Float_tx,y,z;Int_tnlines=0;TFile*f=newTFile(ex44.root,RECREATE);TH1F*h1=newTH1F(h1,xdistribution,100,-4,4);/TNtuple可以看成特殊的的TTree，只可以存放浮点型数据。TTree*mytree=newTTree(“mytree”,”aaaaaaaaaaaa”);mytree-Branch(“x”,略过，但并非不重要,39,读取ASCII格式文件还有个更简便的方式，即用TTree的ReadFile函数：tree-ReadFile(parameter1,parameter2);参见ROOT手册TTree那一章,如何从ASCII文件中读取数据转为ROOT中的TTree$EXAMPLEDIR/ex44a.C,运行：root-lex44a.C或者在root环境中:.xex44a.C,voidex44a()TFile*f=newTFile(ex44.root,RECREATE);TTree*T=newTTree(ntuple,datafromasciifile);/第1个参数为要打开的文件名称/第2个参数是Branch的描述，即设定3个Branchx,y,zLong64_tnlines=T-ReadFile(basic.dat,x:y:z);printf(found%lldpointsn,nlines);T-Write();,40,TChain对象是包含相同tree的ROOT文件的列表。参见UsersGuide中TChain相关章节以及http:/root.cern.ch/root/html534/TChain.html,TChain:分析多个root文件的利器(1),voidex45()/定义TChain，t3为root文件中tree的名称!TChain*fChain=newTChain(“t3”);/添加所有文件至fChain，或根据需要添加部分root文件fChain-Add(“rootfiles/*.root”);/画出t3的某个leaf，如ntrackfChain-Draw(“ntrack”);,注意：此时，fChain等同于一个大root文件中的一个类t3,该文件包含的事例数为所有文件中事例数之和(1012以内),问题：root文件中多个tree怎么办？,41,直方图的操作,直方图的运算加减乘除：Add,Divide,.归一化：ScaleROOT中直方图拟合h1-Fit();,42,直方图归一化(1),归一化后，不仅BinContent变化了，BinError也变化了,直方图的归一化voidTH1:Scale(Double_tc1,Option_t*option),http:/root.cern.ch/root/html534/TH1.html#TH1:Scale,默认c1=1，把直方图每个区间的值(BinContent)乘以c1假设sum=h1-Integral()h1-Scale(c1)之后，h1-Integral()=c1*sum不加参数时，h1-Scale()没有变化(默认c1=1),43,直方图的归一化(2),root0TH1F*h1=newTH1F(h1,100,-5,5);root1TH1F*h2=newTH1F(h2,100,-5,5);root2h1-FillRandom(gaus,5000);root3h2-FillRandom(gaus,10000);root4floatnorm=1000;root5h1-Scale(norm/h1-Integral();root6h2-Scale(norm/h2-Integral();root7h1-Draw(e);root8h2-Draw(“esames”);,归一化之后，h1或h2-Integral()=norm在同一张图上可以看出比较2个分布的差别。,归一化常用于比较两种分布，找出区别。所以，将2个直方图归一化到积分相同进行比较才直观。,注意Draw()函数的选项,44,直方图四则运算(1),重要提示：对直方图进行四则运算操作，一定要想明白运算的意义比如两个直方图的相加与两个随机变量的卷积有什么区别2.两个直方图的四则运算，区间大小和区间数相同才有意义四则运算加减乘除分别对应统计量(BinContent)的相加、相减、相乘、相除3.如果需要正确处理统计误差，需要在对ROOT脚本中调用TH1的某个静态成员函数，即TH1:SetDefaultSumw2();,voidSetDefaultSumw2(Bool_tsumw2=kTRUE)/staticfunction.Whenthisstaticfunctioniscalledwithsumw2=kTRUE,allnewhistogramswillautomaticallyactivatethestorageofthesumofsquaresoferrors,ieTH1:Sumw2isautomaticallycalled.,45,直方图的四则运算(2),46,直方图的四则运算(3),相除,root1TH1F*h3=newTH1F(*h1);root2h3-Sumw2();root3h3-Divide(h1,h2,a,b);root4h3-Divide(h1,h2,a,b);,常用于效率的计算。,相乘,rootTH1F*h3=newTH1F(*h1);rooth3-Sumw2();rooth3-Multiply(h1,h2,a,b);,常用于对分布进行诸如效率等的修正。,思考：如果h1和h2不独立，怎么办？比如h1包含于h2,root4h3-Divide(h1,h2,a,b,B);,47,直方图四则运算的误差处理,虽然ROOT都提供了较完善的一维直方图运算功能，但对最终结果的误差一定要仔细检查。很多情况下，用户需要从图中读出各频数数值与误差值，并确认运算无误。包括归一化和加减乘除在内，如果希望使用直方图的误差，都需要调用TH1:SetDefaultSumw2();或者，对每个直方图(如hist)调用hist-Sumw2();,48,拟合直方图(1),将鼠标放到直方图上，右键，出现直方图操作选项，选择FitPanel，可以在FitPanel中选择拟合的各个选项，比如用什么函数拟合，拟合的区间，等等。,49,用默认的高斯拟合，并在Options菜单中选上FitParameters选项，可以看到拟合的结果。,拟合直方图(2),并不推荐这种拟合方式:1)不适合自定义函数拟合2)不适合批处理,50,hpx-Fit(gaus);hpx-Fit(gaus,-3,3);,拟合直方图(3)$EXAMPLEDIR/ex51.C,自定义拟合函数TF1*fcn=newTF1(fcn,gaus,-3,3);hpx-Fit(fcn,”R”);,gStyle-SetOptFit();/设置拟合选项拟合前往往需要给出合理的参数初值fcn-SetParameters(500,mean,sigma);拟合后取出拟合得到的参数Double_tmypar3;fcn-GetParameters(,运行：root-lroot0.Lex51.Croot1ex51r()root2ex51r2(),用自定义的函数拟合直方图,51,共振峰(Breit-Wigner分布)加上二次函数本底的拟合(一共6个参数)先自定义本底函数(background)和共振峰函数(lorentianPeak)，再定义这两个函数的和为拟合函数:fitFunction利用fitFunction定义TF1函数,拟合直方图(3)$EXAMPLEDIR/ex52.C,这里指定函数区间为0-3，6个参数,运行：root-lroot0.Lex52.C,注意TLegend的使用,TF1*fitFcn=newTF1(fitFcn,fitFunction,0,3,6);,fitFcn-SetParameter(4,0.2);为某个参数设初值(width)fitFcn-SetParLimits(5,0.6,1.4);为某参数设置取值范围,52,ROOT小结,设定ROOT环境变量:ROOTSYS，PATH，LD_LIBRARY_PATH绘制各种直方图，散点图，数学函数TH1F，TH2F，TF1，.随机数产生子，各种分布gRandom-Rndm,Uniform,Gaus,Exp,.创建、保存root文件TFile*f=newTFile(myfile.root”,”recreate);f-Write();TTree,TChain的使用TTree*mytree=newTTree(mytree”,”mytree);mytree-Branch(.);用TChain分析相同格式的数据文件。直方图的运算，拟合h1-Fit(function_name);,53,ROOT的重要功能或用法(1),带参数的macro；ACLIC模式，参考ROOT手册ROOT手册13、14章分别是数学库和线性代数，提供很多数学功能，比如Lorentz矢量的操作，特殊函数，矩阵求解运算，求极值等等ROOT手册第4章介绍Graphs，适用于不等距数据的图形分析(当然也可以构造不等bin的直方图)RooFit，最大似然法拟合等神经网络分析方法,TMVA(多元数据分析)ROOT中使用PYTHIA、Geant3/4图形接口.,54,ROOT的重要功能或用法(2),MakeClass,MakeSelector的运用,比如$EXAMPLEDIR目录下有文件ex51.root，其中含有复杂的tree。可以用MakeClass或MakeSelector自动产生分析文件和头文件：,root0TFilef(ex51.root);root1.lsTFile*ex51.rootTFile*ex51.rootKEY:TTreet4;1Reconsteventsroot2t4-MakeClass();或：t4-MakeClass(“MyClass”);,MakeSelector的用法类似：root0TFilef(ex51.root);root1t4-MakeSelector();或t4-MakeSelector(“MySelector”);,55,直方图画图练习,1.写一个ROOT脚本，ex3_gaus.C,调用随机数产生子产生高斯分布，区间(-6,6)，分30个bin，画出直方图，比较不同的参数的分布。参数组合为：(mean,sigma)=(0,1),(0,2),(1,1