rootで作ったヒストグラムはガウス関数など、任意の関数でFitすることができる。
ROOTにはMINUITという最小化ツールが含まれていて、Fitをするときににはこれを使っている。MINUITはもともとFortranで書かれていたが、ROOT用にc++に書きなおされているのでどちらでも使うことができる。よく知らないがpython版やJava版もあるらしい。
rootには、もともと定義済みの関数が入っており、それらでfitする場合は下のように書くだけでいい。
//histはfitしたいヒストグラム hist->Fit("gaus"); hist->Fit("landau"); hist->Fit("expo");//exp(p0+p1*x) hist->Fit("pol1");//polNでp0 + p1*x + ... + pN*x^N hist->Fit("pol2");
これはTH1だけではなくTGraphやTTreeでもほとんど同じようにFitできる。範囲の指定をしたい場合は下のようにする。
hist->Fit("gaus","","",-1,1);
試しにFit("gaus")とした後、コマンドラインでgausと打ってみるといい。TF1型のオブジェクトが生成されているのが確認できる。 当然、gROOT→FindObject("gaus")でも、このオブジェクトを得ることができる。
一応、Fitのチュートリアルでもかいておきますが、初心者のかたは他のページで勉強して下さい。
[ ] TH1F *hist = new TH1F("name","title",100,-10.,10.); [ ] hist->FillRandom("gaus",1000);//FillRandomは与えられた関数の分布に合わせてFillするメソッド [ ] hist->Fit("gaus");
TF1 *name = new TF1( "const char* name" , "const char* formula" , Double_t xmin , Doble_t xmax);
rootで予め定義されていない関数でFitしたい場合は自分でTF1を定義する必要がある。
例えば下のようにする。
TF1 *f = new TF1("f","[0]*pow(x,[1])",-20,20); hist->Fit(f);
ただし、上のようにするとパラメータが初期化されていないので、きちんとFitされないはず。
自分で定義した関数でfitする場合は、ある程度値を荒くセットしておかないとうまくfitできない。
TF1 *f1 = new TF1("f1","[0]*x+[1]",0,10); f1->SetParameter(0,5); //第1引数が変数の番号、第2引数がその値 f1->SetParameter(1,2);
SetParametersを使えば一度にパラメータをセットできる。
f1->SetParameters(25,10);
f->SetParLimits(1,0.0,1.0);
第一引数はパラメータの番号、第二引数が下限、第三引数が上限。
パラメータの固定も出来る。
f->FixParLimits(1,0.0);
fitする範囲を狭めるとFitがうまくいくことがある。
範囲を制限するには
hist->Fit("hoge","","",min,max);
とする。TF1の定義域は無関係らしい。
gnuplotみたいに関数を描く目的で使うことも当然できる。例えばy=3x+2のような関数を描きたいとする。
TF1 *f = new TF1("name","[0]*x+[1]",-5,5);//-5から+5の範囲で描く。 f->SetParameter(0,3); f->SetParameter(1,2); f->Draw();
ちなみに
TF1 *f = new TF1("name","[0]*x+[1];Xtitle;Ytitle");
などとすれば軸のtitleをつけることができる。そう、TF1の式の記述は実はtitleだったのだー(な、なんだってー!?)。
定義したTF1同士を演算して、TF1を定義することもできる。
TF1 *f1 = new TF1("func1","[0]*x",0,5); TF1 *f2 = new TF1("func2","[0]",0,5); TF1 *f3 = new TF1("func3","func1+func2",0,5);//=[0]*x + [1]になる
TF1 *f1 = new TF1("f2","[0]*TMath::Gaus(x,[1],[2])"); f1->SetParameters(3,25,10);//TMass::Gausは第2引数がmean、第3引数がsigma
#include "TF1.h" double f(double x,double c){ double y; y = pow(x,3)+c; return y; } void file_name(){ TF1* g = new TF1("g","f(x,[0])",0,5); g->SetParameter(0,2); g->Draw(); }
上をrootに読み込ませると、x^3+2がDrawされる。
TH1::Fitの第2引数は下に引用したフィッティングのオプション、第3引数はTH1::Drawのグラフィックオプションである。
"W" 空でないすべての階級に1の重みを設定する。ここではエラーバーは無視する
"WW" 空の階級も含めてすべての階級に1の重みを設定する。エラーバーは無視する
"I" 階級の中央の値を使用する代わりに、関数の積分値を使用する
"L" 対数尤度メソッドを使用する (既定では、カイ2乗を使用する)
"U" ユーザの指定するフィッティングアルゴリズムを使用する
"Q" 静穏モード(最低限の表示のモード)
"V" 冗長モード(既定のモードはQとVの間である)
"E" Minos の技術を使用してより良いエラー推定を行う
"M" フィッティングの結果を改良する
"R" 関数レンジで指定される値域を使用する
"N" グラフィック関数を保存しない、表示をしない
"0" フィッティングの結果を表示しない。既定の設定では、先述の"N"が指定されなければ、フィッティングされた関数が表示される。
"+" フィッティングされた関数のリストに、新しいフィッティングされた関数を追加する (既定では、最新の関数のみが保持されて、それ以前のものは消去される)
"B" 1つ以上のパラメータを修正したくて、フィッティング関数が、polN、expo 、landau、 gausと同じようであれば、このオプションを使用する。
"LL" とても統計量が低くて、内容が整数ではない場合のために改良した、対数尤度によるフィッティング。ビンの内容が(100よりも)大きければ、このオプションは使用しない。
"C" 線形のフィッティングの場合は、カイ二乗を計算しない(時間の節約である)。
"F" polNでフィッティングを行う場合、Minuitフィッターに切り替える。(既定の設定では、polN関数は線形フィッターでフィッティングされる)
TF1はTAttLineを継承しているので、下のように設定を変えることができる。
func->SetLineColor(kRed)
func->SetLineWidth(2)
func->SetLineStyle(2)
あるいはgStyleを使えば、それ以降に生成された関数の設定はすべて変更できる。
[ ] gStyle->SetFuncColor(kRed) [ ] gStyle->SetFuncWidth(2) [ ] gStyle->SetFuncstyle(2)
gausやlandauは下のように色をかえればいい(ただし、一度何らかのメソッドでこれらの関数をつかっていないとできない)。
[ ] gaus->SetLineColor(kRed)
これはオブジェクト名で呼び出しているので、c++の文法では正式には下のようにする。
(TF1*)(gROOT->FindObject("gaus"))->SetLineColor(kRed);
あるいはnameという名前のTF1でfitしたとすると、
hist->GetFunction("name")->SetLineColor(kRed);
などとする方法もある。ただし、Canvasに一度描いてから色をかえる場合は
c1->Modified();
などとしてcanvasの変更を伝える必要がある。
統計ボックスで表示される名前。デフォルトだとp1,p2,…になる。
f->SetParNames("par1","par2","...")
gStyle->SetOptFit(); gStyle->SetOptFit(1111);
引数は10進数の4桁の整数で表され、左から、確率、χ2/Ndof、エラー、パラメータの名前と値を表示する。
例えば、0111でχ2/Ndof、エラー、パラメータの名前と値を表示する(デフォルト)。
fitした関数のパラメータを知りたい場合は統計情報を表示すればディスプレイで確認できるが、
それらの値をマクロで得たい場合は、関数にアクセスして下のようにする。
TF1 *f = new TF1("name","[0]*exp(-0.5*((x-[1])/[2])^2)",-10,10); hist->Fit("name"); Double_t p0 = f->GetParameter(0);//0番目のパラメータ Double_t p0e = f->GetParError(0);//0番目のパラメータのエラー Double_t chi2 = f->GetChisquare();//χ二乗値 Int_t Ndof = f->GetNDF();//自由度の数(フィット点の数 - パラメータ数)
TF1 *fit = hist->GetFunction("name");