RubyROOTで1次元ヒストグラムのテキストデータを読みこんで、TH1D(Double型の1次元ヒストグラム)に変換し、一次関数＋ガウシアンでフィットしてプロットする例を示します。RubyROOTのインストール方法は「RubyROOTをHomebrewでインストールする」を参照してください。
最終的に、以下のようなプロットが得られます。
コード
コードは以下のような構成です。コメントに書いてある通りで、シンプルな構成です。以下のコードをコピーペーストして保存してもよいですが、gistにもアップロードしたので、そちらからダウンロードする方が簡単です。ファイル名はread_and_fit_histogram.rbとしてください。
ヒストグラムのデータはhistogram_sampleからダウンロードしてください。ビンのデータが縦に並んだ超シンプルなデータファイルです。
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
スクリプトの解説
ヒストグラムオブジェクトの作成
#create an instance of histogram
hist=TH1D.create("hist","Histogram",512,0,512)
TH1Dという、1ビンがdoubleで表現されたROOTの標準的なヒストグラムクラスをもとに、ヒストグラムオブジェクト(インスタンス)を作成しています。
ROOTのTH1(1次元ヒストグラム)のAPIリファレンスはこちら→The Histogram classes
テキストファイルの読み込み
標準的なRubyスクリプトのやり方でテキストファイルを開いて、一行ごとに読み込みんでいます。lineに格納された1行の文字列を、.to_i()メソッドを使って整数に変換し、その結果をヒストグラムに格納しています。
もしも入力データが
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のファイルも読み込めます。
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)
}
この場合、ヒストグラムオブジェクトの作成のところでも、入力データの範囲に対応させて
#create an instance of histogram
hist=TH1D.create("hist","Histogram",10,0,100)
のようにして、ビンの分割数、xmin、xmaxを適切に調整する必要があります。
ヒストグラムへのfill
binIndex=hist.GetXaxis().FindBin(i) hist.SetBinContent(binIndex,line.to_i)
データを詰めるべきビン番号を検索して、line.to_i()で整数に変換したビンの値を設定しています。
プロット
#draw with statistical error
canvas=TCanvas.create("c","Canvas",600,600)
hist.Draw("e")
TCanvasのインスタンスを作ってから、TH1のDraw()メソッドに、”e”で表される統計エラーオプションをつけて描画しています。簡単ですね。
フィット関数
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))は、パラメタ何番からをその関数のパラメタとして使うか、ということを表しています。その後ろは関数が定義される区間です。結局、上のインスタンス化は
[0]*exp(-0.5*((x-[1])/[2])**2) + [3] + [4]*x
という関数を作成していることになります。
fit statistics
カイ二乗や自由度を取得する時はTF1クラスのGetChisquare()メソッドやGetNDF()メソッドを使います。
#print fit statistics
puts "NDF = #{f_pol1_gaus.GetNDF()}"
puts "Chi2 = #{"%.2f"%f_pol1_gaus.GetChisquare()}"
NDF = 45 Chi2 = 119.20
実行結果
このスクリプトを実行すると、ターミナルには以下のように表示されます。
途中のp0〜p4がベストフィットパラメタです。
終了するときは、コマンドラインから「.q」と入力してリターンする(ROOTの標準的終了方法)か、GUIのウインドウのFileメニューからQuit ROOTを選択してください。
> 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]
Leave a Reply