查看练习示例脚本
> 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
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]);
}
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.root
root [0] TBrowser a
# 在弹出窗口左侧文件列表中 "ROOT Files" 目录下找到 tree.root 文件,
# 双击 root 文件查看tree "tree",双击 tree 查看 Branch "X" 和 “Y”,
# 双击 Branch “X” 可以在右侧窗口中看到 “X” 的直方图。
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.C
root [0]
# 可以看到生成了一个一维直方图以及二维直方图
# 尝试改变画图选项,得到不同风格的直方图
单击菜单栏中 View 下拉菜单选中 Editor 选项可以调出编辑面版,可以手动编辑图片属性,如坐标轴标题,添加网格,直方图颜色等。尝试用鼠标单击图像不同位置,如直方图,坐标轴,标题处,观察编辑面板编辑选项的变化。尝试手动编辑图像属性。
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 里面找到函数参数的拟合值:
xxxxxxxxxx
NO. 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 Function
root[0] TF1 *fit = histo->GetFunction("fitFcn");
# value of the first parameter
root[1] Double_t p1 = fit->GetParameter(0);
# error of the first parameter
root[2] Double_t e1 = fit->GetParError(0);
# 同时我们还可以计算拟合的 chisquare 值
root[3] Double_t chi2 = fit->GetChisquare();
附加练习
尝试将本例中的朗道分布替换为高斯分布,完成 root 文件生成,root 文件读入,画直方图,拟合直方图。