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