乱数の Fill 3通り(CERN ROOT)

乱数を振ってヒストグラムに積むやり方についてのメモを書いておく。

元々正規分布ガウス分布)に従う乱数を、C++ の標準ライブラリを用いて生成して Fill してというのを for 文で回すみたいなことをしていた。どうやら ROOT の機能だけを用いて簡単に乱数生成と Fill ができるらしいと知ったので自分のコードを書き直した。ROOTで3通り、C++標準ライブラリで1通り下にまとめたつもり。

その際、まだクラスと言う概念を理解しておらず苦戦したので、未来の自分と誰か用のメモとして、ROOT での乱数生成と Fill について書いた。
*1

TH1::FillRandom()

    hist->FillRandom("funcname", number);

とすれば良いらしい。ということで、使い方。

"funcname" としては、自分で定義した TF1 クラスのROOT内の名前を入力する。"gaus" や "landau" などのもともとパラメタ込みで定義されているものも使える。
https://www-he.scphys.kyoto-u.ac.jp/member/n.kamo/wiki/doku.php?id=study:software:root:main#fillrandom

{
    TH1D *hist = new TH1D("", "", 100, 0, 5);
    TF1* func = new TF1("gaussdayo", "TMath::Gaus(x, 2.0, 0.5, true)", 0, 5);  // ここの "gaussdayo" というROOTのための名前が、 TF1 インスタンスに与えられる。ポインタ名 func で呼び出せる。
    hist->FillRandom("gaussdayo", 1000000);  // ROOTのためのTF1インスタンスの名前 "gaussdayo" を入れる。1,000,000個Fillする。
    hist->Draw();
}

ちなみに、"" なしではTF1を定義できなかった。

f:id:nyaaaaoooon:20180922024053p:plain
gaus で画面に出力した画像

TF1::GetRandom()

TF1::GetRandom は、TF1 クラスのメンバ関数なので、TF1 class のポインタにアロー演算子(->)をくっつけて GetRandom() としてやると乱数の値が得られる。

{
    TH1D *histo = new TH1D("", "", 100, 0, 5);
    TF1* funct = new TF1("gausdayo", "TMath::Gaus(x, 2, 0.5, true)", 0., 5.);
    for (Int_t i = 0; i<1000000; ++i){
        Double_t r = funct->GetRandom();  // TF1のポインタ funct に->演算子を付け、GetRandom()を呼び出して乱数の値を得る。その値を r に代入する。 
        histo->Fill(r);
    }
    histo->Draw();
}

とすると同様に以下のようなヒストグラムの画像が得られる。

f:id:nyaaaaoooon:20180922024309p:plain
gausdayoで画面に出力した画像


TRandom::分布名

分布一覧(メンバ関数一覧)はこれ。Sphere とか Circle とかを簡単に作れるのが良いと思った。
ROOT: TRandom Class Reference

TRandom::Gaus() で乱数を得る。

    TRandom3 r; // TRandom を継承するTRandom3 classを用いる。
    
    // generate a gaussian distributed number with:
    double x1 = r.Gaus(); // average = 0, sigma = 1 (初期値)
    double x2 = r.Gaus(10,3); // use average = 10, sigma = 3 のガウス分布

ROOTUsersGuide

となっていた。(乱数の「シード」がどこでどう定義されているかは理解していない。)
これを用いて、

{
    TH1D *histog = new TH1D("", "", 100, 0, 5);
    TRandom3 *rnd = new TRandom3();
    for (Int_t i = 0; i<1000000; ++i){
        Double_t r = rnd->Gaus(2, 0.5);
        histog->Fill(r);
    }
    histog->Draw();
}

とできた。

f:id:nyaaaaoooon:20180922024532p:plain
gausus で画面に出力した画像
(ただ使っただけとは言っても初心者な私にはまだ簡単ではないので参考にしたサイトがあり、これ。)
ROOTでランダムな数を生成するTRandommorephotons.wordpress.com


C++標準ライブラリの を用いる方法。)

ROOT のクラス群を使わずに書いていた時はこうしていた。
にも normal 以外の分布はたくさんあるが、自分の分野でROOTで解析する限り ROOT での乱数生成のほうが楽に書けると思う。どうでしょう。

// normaldist.cxx
#include <random>
using namespace std;    // 毎回std::と打たなくてよくなる

int normaldist(void){

    TH1D *hist = new TH1D("<random>", "", 100, 0, 5);
    
    int n=1000000; // 出力数

    // <random> ライブラリから
    mt19937_64 mt(3);  // mersenne twister の 64 bit ver. のオブジェクトに、シードとして3を入力
    normal_distribution<> norm(2, 0.5); // 平均2, 標準偏差0.5(元にしたプログラムでは分散と書かれていたが、Std Devが0.5程度になったので標準偏差のはず)の正規分布に整えてもらう
    
    for (int i = 0; i<n; ++i){
        hist->Fill(norm(mt));
    }    
    hist->Draw();
    
    return 0;
}

f:id:nyaaaaoooon:20180922023322p:plain
normaldist で出力した画像



Fit 用のTF1::Gausのメモ。

関数のパラメタを後から代入してFillしないだけのものも置いておく。2行増えてFillしないだけ。

void gaus(){
    TH1D *hist = new TH1D("", "", 100, 0, 5);
    TF1* func = new TF1("gaussda", "TMath::Gaus(x, [0], [1], true)", 0, 5);  // "gaussda"と言う名前を TF1 のオブジェクトに与える。"" なしでは定義できなかった。
    func->SetParameter(0, 2);   // TMath::Gaus の [0] の所にはガウス分布の平均値が入る。
    func->SetParameter(1, 0.5);   // [1] の所にはガウス分布の標準偏差が入る。
}

*1:注意としてたまに書かれていること:TRandom クラスが出力する擬似乱数列は、目的によっては性質がよくないので使えない場合もある。以下の記事のTH1やTF1の持つ乱数生成器はTRandom から来ているので、そのままでは用いることが出来ないことがある。(if(gRandom) delete gRandom; gRandom = new TRandom3(0); とせよというQ&Aを見かけたのでそうすることにしている)