ユーザ用ツール

サイト用ツール


サイドバー

Menu

Latest

study:software:root:roofit

RooFit

RooFitとは

RooFitはFitツールの一種。

ヒストグラムやグラフに関数をfitするのは普通のROOTの機能で十分である。RooFitは以下のような特徴を持つ。

  • 最尤法(Most Likelihood method)なので、binnnigしていないdataにfitができる
    • もちろん最小二乗法も使える
    • pdfとdataからlikelihood functionを得ることもできるらしい
  • 激しくオブジェクト指向で構成されていて、変数やpdfなどがすべてc++オブジェクトになっている
    • ROOTの独特なシステムと違って、普通のオブジェクト指向プログラミングに感じる
    • でも正直、使いにくい
  • pdfの規格化ができる
  • 多次元fit
  • 多次元関数の積分
  • Toy Montecarlo Generationができる
  • 2つのpdfを合成して新たなpdfを作ることができる
    • 足して一次元pdfを作る
      • f(x)+g(x)
    • 掛け合わせて二次元pdfを作る
      • f(x)*g(y)
    • parameterを他のpdfで変化させて二次元pdfを作る
      • f(x,g(y))
    • convolutionを行う。高速フーリエ変換を使用できる。

ほとんど、ROOTとは独立した機能なので、特殊なpdfを使う以外の意味では使う必要はないと思う。

チュートリアル

$ROOTSYS/tutorials/roofit/rf101_basics.Cより

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooPlot.h"
using namespace RooFit ;


//RooFitでは変数もオブジェクト
//ROOTと同じで第1引数はname、第2引数はtitle。つまりなんでもいい
//第3引数はパラメータの下端、第4引数は上端
RooRealVar x("x","x",-10,10) ;


//引数が5つの場合は第3引数はパラメータの固定値
RooRealVar mean("mean","mean of gaussian",1,-10,10) ;
RooRealVar sigma("sigma","width of gaussian",1,0.1,10) ;


//予めRooFitではpdfがいくつか登録されている
RooGaussian gauss("gauss","gaussian PDF",x,mean,sigma) ;//RooRealVarはすべて参照渡しにする

//MC genration
RooDataSet* data = gauss.generate(x,10000) ;

//fit
gauss.fitTo(*data) ;//参照渡しに注意

  • x,mean,gaussはすべてRooRealVar
    • 違いはxだけ値が上限と下限を決めて、固定値をもっていないこと
    • つまりRooRealVarはパラメータと定数を両方扱えるクラス
      • なので、もちろんxの代わりにmeanを変数にしてグラフを描くこともできる
  • RooDataSetは任意の次元(たぶん上限はあると思う)のunbinned dataを保持できる
    • RooPlotにplotするときには自動的にbinngingされて描画されている
  • RooFitのメソッドはいわゆるlowerCamelCaseで書かれている
    • ROOTのメソッドはUpperCamelCaseで書かれているので区別できる

RooPlot* xframe = x.frame(Title("Gaussian p.d.f.")) ;


//pdfをxframeにplotする
gauss.plotOn(xframe) ;//ここではポインタ渡しをする。ややこしい。
//色をつける場合は、
//gauss.plotOn(xframe,LineColor(kRed)) ;

//dataのplot(binnningされた状態で表示される)
data->plotOn(xframe) ;
//統計情報の表示
data->statOn(xframe) ;

//TCanvasに描く
xframe->Draw() ;

  • RooFitではある変数のplotを作りCanvasに描くという手順を取る
  • RooPlotはある変数xを横軸にした時の縦軸の値をplotできる
    • RooRealVar::frameはRooPlot*を返り値にもつ
  • Title()というのは"RooGlobalFunc.h"で定義されている関数の一つで、RooFitのオブジェクトにtitleを与えることができる
    • もちろん、使わなくてもいい
  • TCanvasに描かれたオブジェクトはオブジェクトが死ぬと、キャンバスから消える
    • RooPlotの場合、一度plotすれば元のデータは失われても問題ない(RooPlot自体が情報を保持している)
      • なんなんだろう、ROOTの仕様が気に入らない人が設計したのか。なぜ、わざわざ…

しかし、ROOTなら、Draw()の一行で済ませられたことを、RooFitでは3-4行は絶対かかってしまう。さすがにやってられないので、一つの関数でまとめた方がいいと思う。

継承関係、RooAbsArgとRooAbsData

RooFitで重要なクラスの継承関係を示す。

ROOTに例えるなら上で囲んだカテゴリーはTF1で、下で囲んだカテゴリーはTTree、TH1に対応していると思えば理解しやすいと思う。

            RooAbsArg
                |
            RooAbsReal
            /        \                    
RooAbsPdf(pdf)     RooAbsRealLValue(変数)    (抽象クラス)
     |                      |               --------------------
RooGaussian            RooRealVar                 (実体)
...                    RooLinearVar    
                 RooAbsData
          /           |         \                    
RooDataHist      RooDataSet     RooTreeData

AbsというのはAbstruct(抽象クラス)を表している。

RooAbsArgやRooAbsDataはTNamedを継承している。つまり普通のRootのオブジェクトと同等の性質を兼ね備えている。

多次元の場合の使い方、RooArgSetとRooArgList

RooArgSetとRooArgListは共にRooAbsCollectionを継承していて、RooAbsArg(つまり変数や関数)のコンテナ。RooArgListは順序がある点だけ異なるが、どちらもほぼ同じように使うことができる。

例えば、

RooDataSet data("name","title",x);//一次元
RooDataSet data("name","title",RooArgSet(x,y));//二次元、RooArgListでも(たぶん)同じ
RooDataSet data("name","title",RooArgSet(x,y,z));//三次元

のように、一次元ではRooRealVarを入れていたところにRooArgSetを使うことで二次元データにすることができる。詳細は各項目に譲る。

多次元のdataに1次元pdfでfitした場合は、その変数への射影にfitを行う。多次元pdfでfitした場合は多次元fitも可能。

(特に知る必要はないが)RooDataSetのコンストラクタの挙動を説明する。

コンストラクタは

	RooDataSet(const char* name, const char* title, const RooArgSet& vars, const char* wgtVarName = 0)

と定義されている(他にも定義されているが、これを例にする)。

第3引数はRooArgSetだけど、RooArgSetのコンストラクタには

	RooArgSet(const RooArgList& list)
	RooArgSet(const RooAbsArg& var1, const char* name = "")

というのが定義されていて、RooArgList、RooAbsArgからRooArgSetへのキャスト(変換コンストラクタかな)が暗黙に可能(c++では引数が一つの場合はこういうことができる)。

なので、結果としてRooArgSet、RooArgList、RooAbsArgのいずれも代入できる。

ちなみにRooDataSetの場合はRooAbsPdf系を代入してもerrorになる(コンパイルはできる)。

RooCmdArg

RooCmgArgというのは(よくわかっていないが)オブジェクトの色とか線のサイズといったオプションを表すクラスらしい。

RooFitのメソッドはRooCmgArgを引数にすることができて、順不同にこれらをいれていけばオブジェクトの細かい設定ができる。

RooGlobalFunc.hで定義されているRooFit::TitleやRooFit::LineColorの返り値として生成される。

下に幾つか例をあげる。

frame

RooRealVar::frameの引数はすべてRooCmdArg。これらは好きな順番で入れられる。

RooPlot* xframe = x.frame(Name("..."),Title("..."),Bins(10),Range(-10.,10.));

plotOn

RooAbsPdf::plotOnの第二引数以降はRooCmdArg。

gauss.plotOn(xframe,LineColor(kRed),Binning(100));

fitTo

gauss.fitTo(data,Range(-10.,10.));

定数

RooFitでは変数もRooAbsRealクラスを継承している必要があるわけだけど、RooRealVar(変数)ではなく、定数を代入した場合もあるはず。

RooConst(3);

RooConstはRooConstArgクラスを返り値に持っている。RooConstArgはRooAbsRealを継承している。

てゆうか、これこそ変換コンストラクタつかって、定数をそのまま代入するようにすればいいのに。全く不親切な設計。

plotOn

RooRealVarが失われた場合

繰り返すが、RooFitではDataの可視化をヒストグラムを使わずに、RooPlotを使用する。

RooRealVar x("x","x",-10,10);
RooPlot *framex = x.frame();

ところで、ここでスコープから外れて、RooRealVar xが失われた場合にどうやって可視化するのか疑問を起こすと思う。

RooRealVar x2("x","hoge",-10,10);
RooPlot *framex = x2.frame();

として、元の変数と同じオブジェクト名でRooRealVarを宣言して、RooPlotを生成すれば、同様にできる。

RooDataSetからRooRealVarを取り出すこともできるが、あまり想定していない気がする。どちらにしても、オブジェクト名が必要。

ようするにplotOnするときに本当に必要なのはRooRealVarではなく、オブジェクト名(const char*)である。この辺りにRooFitの設計の悪さがうかがえる。

というか無理やり違う思想の体型をROOTにいれているから、こういう変な現象が起きている気もする。

多次元の場合の挙動について、Rangeの制限

RooDataSetは多次元データを保持できる。持っている情報量としてはRooDataSetとTTreeはあまり違いはない(TTreeはクラスをBranchにできるが)。

RooPlotは何故かは知らないが、1次元しか可視化することができない。なので、plotした場合は1変数への射影を描く。

RooPlot *xframe = x.frame();
data->plotOn(xframe);

射影する場合、他の変数の範囲を制限したい場合があると思う。

例えば、dataにはx,y,zの値が入っていた場合は

x.setRange("selection",xmin,xmax);
y.setRange("selection",ymin,ymax);
z.setRange("selection",zmin,zmax);
data->plotOn(xframe,CutRange("selection"));
//data->plotOn(xframe,CutRange("selection,selection2")); ,でor条件で連結できるみたい

などとする。ここで注意するのは変数の値そのもののRangeは変わらないということ。

変数のRangeを変えたい場合は、

x.setRange(xmin,xmax);

とする。

Binning

Binning(nbin,low,high);

const int nbin = n;
double bin[nbin+1] = {...};
Binning(RooBinning(nbin,bin));

dataの変換

RooAbsData->TH1

RooPlotにしなくてもTH1を作って可視化することもできる。

  //dataはRooAbsDataを継承しているクラス
  TH1 *h = data.createHistogram("h",x,Binning(100));
  TH2 *h2 = (TH2*)data.createHistogram( "h2",x,Binning(100), YVar(y,Binning(100)) );
  TH3 *h3 = (TH3*)data.createHistogram( "h3",x,Binning(100), YVar(y,Binning(100)) , ZVar(z,Binning(100)) );

  • 一次元ではRooPlotを作るのとTH1を作るのにはそれほど差はないが、RooPlotはなぜか二次元以上を可視化することはできないのでTH1を作る必要がある
  • 4次元以降はサポートされていないと思われる。まあ可視化目的以外ではTH1に詰め直す必要性はないのでとりあえず無問題だと思う
  • createHistogramの第3引数以降はRooCmdArg。たぶん順不同。

TH1->RooDataHist

RooDataHistはRooAbsDataを継承したクラスでbinned dataを保存できる。

RooRealVar x("x","x",-10,10) ; 
TH1* h = ...; 
RooDataHist data("data","dataset with x",x,h) ; 

TTree->RooDataSet

  TTree* tree = (TTree*)gDirectory->Get("tree");
  //treeにはx,yというBranchがあるとする

  RooRealVar x("x","x",xmin,xmax) ;  
  RooDataSet data1("data1","dataset with x",tree,x);

  //二次元以上ではRooArgSetを使う
  RooRealVar x("y","y",ymin,ymax) ;  
  RooDataSet data2("data2","dataset with x,y",tree,RooArgSet(x,y));

直接値をいれてRooDataSetを作る

  //一次元
  RooRealVar x("x","x",xmin,xmax) ;  
  RooDataSet data1("data1","dataset with x",x);
  for(int i=0;i<1000;++i){
  
    x.setVal(gRandom->Rndm());
    data1.add(x);
    //data1.add(x,weight=1,weighterror=0);
  }

  //二次元
  RooRealVar y("y","y",ymin,ymax) ; 
  RooDataSet data2("data2","dataset with x,y",RooArgSet(x,y));
  for(int i=0;i<1000;++i){
  
    x.setVal(gRandom->Rndm());
    y.setVal(gRandom->Rndm());
    data2.add(RooArgSet(x,y));

  }

RooAbsArg(関数、データ)から条件を満たすデータをRooAbsArgを生成する

RooDataSet* data2 = (RooDataSet*)data->reduce("条件式");

条件式は、x>-2 && y < 0という感じ。pdfの場合も同様にreduce可能。

RooDataSetからRooArgSetを取り出す

RooArgSet *as = data->get(number); //RooDataSet* data

number番目の値を返す。引数なしだと、現在の値を返す。

RooArgSetから値を取り出す

"name"という名前の変数の値を取り出す

as->getRealValue("name") //RooArgSet *as

RooArgSetからRooRealVarを取り出す

"name"という名前の変数を取り出す場合。

//asはRooArgSetの実体とする
RooRealVar *x = (RooRealVar*)&(as["name"]); // オペレータ[]はアドレスではなく参照返し。

RooRealVarから値、エラーを取り出す

Double_t num = x->getValV();
Double_t error = x->getError();

2つのpdfを組み合わせる

addition

convolution

  • convolution(畳み込み)を行う
  • デフォルトでは高速フーリエ変換で計算
  • 1e-7の精度(あまり詳しいことは知らない)

    RooRealVar x("x","x",-10,10) ;

    RooRealVar meanl("meanl","mean of Landau",2) ;
    RooRealVar sigmal("sigmal","sigma of Landau",1) ;
    RooLandau landau("landau","landau",x,meanl,sigmal) ;

    RooRealVar meang("meang","mean of Gaussian",0) ;
    RooRealVar sigmag("sigmag","sigma of Gaussian",2) ;
    RooGaussian gauss("gauss","gauss",x,meang,sigmag) ;

    RooNumConvPdf model("model","model",x,landau,gauss) ;

    RooPlot* frame = x.frame() ;
    model.plotOn(frame) ;
    landau.plotOn(frame,LineStyle(kDashed)) ;
    frame->Draw() ;

convolutionの計算は[-inf, +inf]で積分するのでリソースをかなり使う。

例えばgaussianならtailの裾で0に近づくので、数sigmaの範囲で計算すれば十分なので範囲を制限したほうがいい。

下のようにすると[meang-5*sigmag, meang+5*sigmag]で計算できる。

model.setConvolutionWindow(meang,sigmag,5); 

このメソッドは実体を定義してからでも意味あるのだろうか。

よくわからないけど、積分計算自体を実行するのはplotOnとかfitToとかするときなのかなあ。積分計算を解析的にするわけではないので、積分結果がfitするようにパラメータを動かすってことか。めちゃくちゃ計算量多い気がする、

production

    RooRealVar x("x","x",-10,10) ;
    RooRealVar meanx("meanx","meanx",0,-10,10) ;
    RooRealVar sigmax("sigmax","sigmax",3,0.,10.) ;
    RooGaussian gaussx("gaussx","gaussx",x,meanx,sigmax) ;

    RooRealVar y("y","y",-10,10) ;
    RooRealVar meany("meany","meany",0,-10,10) ;
    RooRealVar sigmay("sigmay","sigmay",2,0.,10.) ;
    RooGaussian gaussy("gaussy","gaussy",y,meany,sigmay) ;

    RooProdPdf gaussxy("gaussxy","gaussxy",RooArgSet(gaussx,gaussy)) ;

f(x,m)->f(x;m(y,a))

あるpdfのparameterをyの関数として動かして二次元分布を作る


 //gaussianの用のparameter定義
  RooRealVar x("x","x",-10,10) ;
  RooRealVar sigma("sigma","width of gaussian",3,0.1,10) ;

  //linearのparameter定義
  RooRealVar y("y","y",-10,10) ;
  RooRealVar a0("a0","a0",0,-10,10) ;
  RooRealVar a1("a1","a1",3,-10,10) ;

  //RooPolyVarは任意の次元のpolynominalを作れる
  //この場合はm=a0+a1*y
  RooPolyVar m("m","m",y,RooArgList(a0,a1)) ;


  //mはRooPolyVarであることに注意。つまりmが変数yに依存している
  RooGaussian gauss("gauss","gaussian PDF",x,m,sigma) ;

  //二次元の場合はRooPlotで可視化することはできない(RooPlotは射影に使う)
  //代わりにTH2を使うことにする
  TH2 *h = gauss.createHistogram( "h",x,Binning(100), YVar(y,Binning(30)) );

  h->Draw("lego2 z") ;

  • RooGaussianの第4引数以降はRooRealVarである必要はなく、RooAbsRealを継承していればいい。
    • 噛み砕いていうと、pdfのパラメータは変数だけでなく、関数を使ってもいい

general function

RooGenericPdf

    //double func(double,double,double);
    RooRealVar x("x","x",0,-5,5);
    RooRealVar a("name1","title",0,-5,5);
    RooRealVar b("name2","title",0,-5,5);
    RooGenericPdf g("g","func(x,name1,name2)",RooArgSet(x,a,b));

ここで作られるpdfはnormalizationされない。fitする上ではそれほど問題ないと思うかもしれないが、normalizationできないとヒストグラムに関数型をあわすことができない(つまり見た目でfitがうまくいったかどうか判断しにくくなる)。

normalizationするためには解析的な積分関数を自分で定義しておく必要があるが、RooGenericPdfではたぶん無理。

RooClassFactory

c++/TF1 binding

#include "RooTFnBinding.h"

  TF1 *f = new TF1("f","sin(x)/x",0,10);
  RooRealVar x("x","x",0.01,20) ;
  RooAbsReal* rf = bindFunction(f,x) ;

  • 多変数の場合はTF2、TF3を使う
  • パラメータ([0],[1],…)を定義式に使ってもいいが、固定しておかないと意味ない

//double f(double,double,double);
  RooRealVar x(...) ;
  RooRealVar a(...) ;
  RooRealVar b(...) ;
  RooAbsReal* rf = bindFunction("name",f,x,a,b) ;

使用する関数の引数の数に応じて、ヘッダーをincludeする必要がある。

#ifndef __CINT__
#include "RooCFunction1Binding.h"
#include "RooCFunction2Binding.h"
#include "RooCFunction3Binding.h"
#include "RooCFunction4Binding.h"
#endif

引数5個以上は今のところ用意されていない。

#include "RooTFnPdfBinding.h"

  RooAbsPdf* rf = bindPdf(f,x) ;

TF1のパラメータ(i.e.定義式中の[0],[1],…のこと)はたぶん定数扱い。

  RooAbsPdf* rf = bindPdf("name",f,x,a,b) ;

  • pdfは正の値しかとらない
  • 規格化されるわけではない
study/software/root/roofit.txt · 最終更新: 2014/02/18 08:58 by kamo

ページ用ツール