ROOT 分析作业
查看练习示例脚本
> cd root_exercise/> ls-rw-r--r-- 1 user group 1792 Aug 16 19:45 exercise0.C-rw-r--r-- 1 user group 969 Aug 16 19:49 exercise1.C-rw-r--r-- 1 user group 2848 Aug 16 19:40 exercise2.C
1，生成数据样本
1.1 定义分布函数
Lorentzian Peak function:
x // Quadratic background function Double_t background(Double_t *x, Double_t *par) { return par[0] + par[1]*x[0] + par[2]*x[0]*x[0]; } // Lorentzian Peak function Double_t lorentzianPeak(Double_t *x, Double_t *par) { return (0.5*par[0]*par[1]/TMath::Pi()) / TMath::Max(1.e-10, (x[0]-par[2])*(x[0]-par[2])+ .25*par[1]*par[1]); } // Sum of background and peak function Double_t Function(Double_t *x, Double_t *par) { return background(x,par) + lorentzianPeak(x,&par[3]); }1.2 生成数据样本，填tree，写入root文件
x
// create the file, the Tree and a few branches *f = ::Open("tree.root","recreate"); TTree *tree = new TTree("tree","treelibrated tree"); Float_t x, y; // create a branch with energy tree->Branch("X",&x); tree->Branch("Y",&y); // create a TF1 with the range from 0 to 3 and 6 parameters TF1 *Fcnx = new TF1("Fcnx", Function, 0,3,6); TF1 *Fcny = new TF1("Fcny", Function, 0,3,6); // Fix the Parameters of function Fcnx->SetParameters(-1, 45, -13.3, 13.8, 0.2, 1); Fcny->SetParameters(-1, 45, -13.3, 13.8, 0.1, 1.5); // fill some events with random numbers Int_t nevent=1000; for (Int_t iev=0;iev<nevent;iev++) { x = Fcnx->GetRandom(); y = Fcny->GetRandom(); tree->Fill(); // fill the tree with the current event } // save the Tree header. The file will be automatically closed // when going out of the function scope tree.Write();操作练习
xxxxxxxxxx> root -l exercise0.C root [0] Processing exercise0.C...root [1] .q# 查看生成的 root 文件> ls -l -rw-r--r-- 1 user group 1792 Aug 16 19:45 exercise0.C-rw-r--r-- 1 user group 969 Aug 16 19:49 exercise1.C-rw-r--r-- 1 user group 2848 Aug 16 19:40 exercise2.C-rw-r--r-- 1 user group 76796 Aug 16 20:03 tree.root# 打开 root 文件, 通过 TBrowser 查看内容> root -l tree.rootroot [0] TBrowser a # 在弹出窗口左侧文件列表中 "ROOT Files" 目录下找到 tree.root 文件，# 双击 root 文件查看tree "tree"，双击 tree 查看 Branch "X" 和 “Y”，# 双击 Branch “X” 可以在右侧窗口中看到 “X” 的直方图。
2. 读取root文件，填直方图，画出直方图
x
*f = ::Open("tree.root"); TTree *t1 = (TTree*)f->Get("tree"); Float_t x, y; t1->SetBranchAddress("px",&px); t1->SetBranchAddress("py",&py); // create two histograms TH1F *hx = new TH1F("hx","x distribution",100,-3,3); TH2F *hxy = new TH2F("hxy","y vs x",30,-3,3,30,-3,3); // read all entries and fill the histograms Long64_t nentries = t1->GetEntries(); for (Long64_t i=0;i<nentries;i++) { t1->GetEntry(i); hx->Fill(x); hxy->Fill(x,y); } // draw the histograms TCanvas *c1 = new TCanvas(); hx->Draw("E0"); c1->Print("hx.png") // eps, ps, jpg TCanvas *c2 = new TCanvas(); hxy->Draw("colz"); // BOX, SCAT, ARR c1->Print("hxy.png") 操作练习
xxxxxxxxxx> root -l exercise1.Croot [0]# 可以看到生成了一个一维直方图以及二维直方图# 尝试改变画图选项，得到不同风格的直方图
单击菜单栏中 View 下拉菜单选中 Editor 选项可以调出编辑面版，可以手动编辑图片属性，如坐标轴标题，添加网格，直方图颜色等。尝试用鼠标单击图像不同位置，如直方图，坐标轴，标题处，观察编辑面板编辑选项的变化。尝试手动编辑图像属性。
3. 拟合直方图
x
// Read the Tree *f = ::Open("tree.root"); TTree *t1 = (TTree*)f->Get("tree"); TCanvas *c1 = new TCanvas("c1","Fitting Demo",10,10,700,500); c1->SetFillColor(33); c1->SetFrameFillColor(41); c1->SetGrid(); // Create one histogram TH1F *histo = new TH1F("histo", "Lorentzian Peak on Quadratic Background",60,0,3); histo->SetMarkerStyle(21); histo->SetMarkerSize(0.8); histo->SetStats(0); // Fill the histograms t1->Project("histo", "X"); // create a TF1 with the range from 0 to 3 and 6 parameters TF1 *fitFcn = new TF1("fitFcn",fitFunction,0,3,6); fitFcn->SetNpx(500); fitFcn->SetLineWidth(4); fitFcn->SetLineColor(kMagenta); // first try without starting values for the parameters // This defaults to 1 for each param. // this results in an ok fit for the polynomial function // however the non-linear part (lorenzian) does not // respond well. fitFcn->SetParameters(1,1,1,1,1,1); histo->Fit("fitFcn","0"); // second try: set start values for some parameters fitFcn->SetParameter(4,0.2); // width fitFcn->SetParameter(5,1); // peak histo->Fit("fitFcn","V+","ep"); // improve the picture: TF1 *backFcn = new TF1("backFcn",background,0,3,3); backFcn->SetLineColor(kRed); TF1 *signalFcn = new TF1("signalFcn",lorentzianPeak,0,3,3); signalFcn->SetLineColor(kBlue); signalFcn->SetNpx(500); Double_t par[6]; // writes the fit results into the par array fitFcn->GetParameters(par); backFcn->SetParameters(par); backFcn->Draw("same"); signalFcn->SetParameters(&par[3]); signalFcn->Draw("same"); // draw the legend TLegend *legend=new TLegend(0.6,0.65,0.88,0.85); legend->SetTextFont(72); legend->SetTextSize(0.04); legend->AddEntry(histo,"Data","lpe"); legend->AddEntry(backFcn,"Background fit","l"); legend->AddEntry(signalFcn,"Signal fit","l"); legend->AddEntry(fitFcn,"Global Fit","l"); legend->Draw();操作练习
xxxxxxxxxx> root -l exercise1.C得到拟合结果：
可以在输出 log 里面找到函数参数的拟合值：
xxxxxxxxxxNO. NAME VALUE ERROR NEGATIVE POSITIVE 1 p0 -6.98862e+00 1.03331e+00 2 p1 2.36951e+02 5.46818e+00 3 p2 -6.93593e+01 2.12246e+00 4 p3 8.09132e+01 5.59809e+00 5 p4 2.05657e-01 1.77811e-02 6 p5 1.01006e+00 5.72858e-03 也可以通过以下方式拿到单个参数值及其误差：
xxxxxxxxxx# Get Associated Functionroot[0] TF1 *fit = histo->GetFunction("fitFcn");# value of the first parameterroot[1] Double_t p1 = fit->GetParameter(0);# error of the first parameterroot[2] Double_t e1 = fit->GetParError(0);# 同时我们还可以计算拟合的 chisquare 值root[3] Double_t chi2 = fit->GetChisquare();
附加练习
尝试将本例中的朗道分布替换为高斯分布，完成 root 文件生成，root 文件读入，画直方图，拟合直方图。