【RubyROOT】ヒストグラムをフィットする例

RubyROOTで1次元ヒストグラムのテキストデータを読みこんで、TH1D(Double型の1次元ヒストグラム)に変換し、一次関数＋ガウシアンでフィットしてプロットする例を示します。RubyROOTのインストール方法は「RubyROOTをHomebrewでインストールする」を参照してください。

最終的に、以下のようなプロットが得られます。

コード

コードは以下のような構成です。コメントに書いてある通りで、シンプルな構成です。以下のコードをコピーペーストして保存してもよいですが、gistにもアップロードしたので、そちらからダウンロードする方が簡単です。ファイル名はread_and_fit_histogram.rbとしてください。

ヒストグラムのデータはhistogram_sampleからダウンロードしてください。ビンのデータが縦に並んだ超シンプルなデータファイルです。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40

require "RubyROOT"
include Root
include RootApp

inputFile="histogram_sample.data"
outputPDFFile="histogram_sample.pdf"

#create an instance of histogram
hist=TH1D.create("hist","Histogram",512,0,512)

#read text data that contains counts of each histogram bin
open(inputFile).each_with_index(){|line,i|
	binIndex=hist.GetXaxis().FindBin(i)
	hist.SetBinContent(binIndex,line.to_i)
}

#draw with statistical error
canvas=TCanvas.create("c","Canvas",600,600)
hist.Draw("e")

#fit the 40K line
f_pol1_gaus=TF1.new("f_pol1_gaus","gaus(0)+pol1(3)",150,200)
f_pol1_gaus.SetParameter(0,200)
f_pol1_gaus.SetParameter(1,180)
f_pol1_gaus.SetParameter(2,10)
f_pol1_gaus.SetParameter(3,200)
f_pol1_gaus.SetParameter(4,-1)
f_pol1_gaus.SetLineColor(KOrange)
hist.Fit("f_pol1_gaus","","",150,200)

#save as pdf
hist.SaveAs(outputPDFFile)

#print fit statistics
puts "NDF  = #{f_pol1_gaus.GetNDF()}"
puts "Chi2 = #{"%.2f"%f_pol1_gaus.GetChisquare()}"

#run ROOT GUI event loop
canvas.Update
run_app

スクリプトの解説

ヒストグラムオブジェクトの作成

1
2

#create an instance of histogram
hist=TH1D.create("hist","Histogram",512,0,512)

TH1Dという、1ビンがdoubleで表現されたROOTの標準的なヒストグラムクラスをもとに、ヒストグラムオブジェクト(インスタンス)を作成しています。

ROOTのTH1(1次元ヒストグラム)のAPIリファレンスはこちら→The Histogram classes

[browser-shot url=”https://root.cern.ch/root/html/TH1.html” width=”300”]

テキストファイルの読み込み

標準的なRubyスクリプトのやり方でテキストファイルを開いて、一行ごとに読み込みんでいます。lineに格納された1行の文字列を、.to_i()メソッドを使って整数に変換し、その結果をヒストグラムに格納しています。

もしも入力データが

1
2
3
4
5
6
7

10.0  2.0
20.0  3.5
30.0  5.0
40.0 7.9
50.0 12.0

のようにビンの値とカウントのペアになっているときは、読み込み部分は以下のようにString.split(delimitter)メソッドを用いて1行を空白文字で分割して配列に分けてから取り扱います。delimitterを”,”に変更すればCSVのファイルも読み込めます。

1
2
3
4
5
6
7
8
9

open(inputFile).each(){|line|
	#split with white space
	array=line.split(" ")
	#retrieve bin value and count
	binValue=array[0].to_f
	count=array[1].to_f
	binIndex=hist.GetXaxis().FindBin(binValue)
	hist.SetBinContent(binIndex,count)
}

この場合、ヒストグラムオブジェクトの作成のところでも、入力データの範囲に対応させて

1
2

#create an instance of histogram
hist=TH1D.create("hist","Histogram",10,0,100)

のようにして、ビンの分割数、xmin、xmaxを適切に調整する必要があります。

ヒストグラムへのfill

1
2

binIndex=hist.GetXaxis().FindBin(i)
hist.SetBinContent(binIndex,line.to_i)

データを詰めるべきビン番号を検索して、line.to_i()で整数に変換したビンの値を設定しています。

プロット

1
2
3

#draw with statistical error
canvas=TCanvas.create("c","Canvas",600,600)
hist.Draw("e")

TCanvasのインスタンスを作ってから、TH1のDraw()メソッドに、”e”で表される統計エラーオプションをつけて描画しています。簡単ですね。

フィット関数

1

f_pol1_gaus=TF1.new("f_pol1_gaus","gaus(0)+pol1(3)",150,200)

でフィットに使う関数を定義しています。

TFormulaのAPIリファレンスにあるようにgausは

gaus(0) is a substitute for [0]*exp(-0.5*((x-[1])/[2])**2) and (0) means start numbering parameters at 0

です。同様にpol1は

[0]+[1]*x

です。

関数の後ろにカッコ書きされた数字(gaus(0)とかpol1(3))は、パラメタ何番からをその関数のパラメタとして使うか、ということを表しています。その後ろは関数が定義される区間です。結局、上のインスタンス化は

1

[0]*exp(-0.5*((x-[1])/[2])**2) + [3] + [4]*x

という関数を作成していることになります。

fit statistics

カイ二乗や自由度を取得する時はTF1クラスのGetChisquare()メソッドやGetNDF()メソッドを使います。

1
2
3

#print fit statistics
puts "NDF  = #{f_pol1_gaus.GetNDF()}"
puts "Chi2 = #{"%.2f"%f_pol1_gaus.GetChisquare()}"

1
2

NDF  = 45
Chi2 = 119.20

実行結果

このスクリプトを実行すると、ターミナルには以下のように表示されます。 途中のp0〜p4がベストフィットパラメタです。

終了するときは、コマンドラインから「.q」と入力してリターンする(ROOTの標準的終了方法)か、GUIのウインドウのFileメニューからQuit ROOTを選択してください。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31

> ruby read_and_fit_histogram.rb
  *******************************************
  *                                         *
  *        W E L C O M E  to  R O O T       *
  *                                         *
  *   Version   5.34/26  20 February 2015   *
  *                                         *
  *  You are welcome to visit our Web site  *
  *          http://root.cern.ch            *
  *                                         *
  *******************************************

ROOT 5.34/26 (v5-34-26@v5-34-26, Feb 20 2015, 13:23:25 on macosx64)

CINT/ROOT C/C++ Interpreter version 5.18.00, July 2, 2010
Type ? for help. Commands must be C++ statements.
Enclose multiple statements between { }.
 FCN=119.203 FROM MIGRAD    STATUS=CONVERGED     219 CALLS         220 TOTAL
                     EDM=7.61759e-09    STRATEGY= 1      ERROR MATRIX ACCURATE
  EXT PARAMETER                                   STEP         FIRST
  NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE
   1  p0           1.13034e+02   4.30109e+00   1.75476e-02   1.88265e-06
   2  p1           1.72250e+02   3.15708e-01   1.38694e-03  -1.17718e-04
   3  p2           7.21296e+00   3.16813e-01   9.66173e-04  -4.12744e-04
   4  p3           1.85431e+02   1.85117e+01   6.20558e-03  -2.53449e-05
   5  p4          -7.83916e-01   9.76740e-02   3.42935e-05  -3.68753e-03
Info in <TH1D::SaveAs>: C++ Macro file: histogram_sample.pdf has been generated
NDF  = 45
Chi2 = 119.20
To quit: please select "Quit ROOT" from the pulldown menu "File"
root [0]