ROOTのグラフを使う

概要

欧州原子核研究機構CERNが開発しているROOTのグラフについて記述しようと思います。
グラフは散布図を描くのにも使えるし、曲線を描くのにも使える便利なクラスです。
この記事で扱う内容は以下の通りです。

• グラフを定義する
• グラフを描く
• グラフをフィットする
• グラフでエラーバンドを描く

公式のクラスページを参考にしています。

グラフを定義する

グラフは変数の数が2つか3つか、そしてエラーがあるかで細かくクラスが分かれています。
まとめると次の表のようになります。

	エラーなし	エラーあり （対称）	エラーあり （非対称）
２変数(x, y)	TGraph	TGraphErrors	TGraphAsymErrors
3変数(x, y, z)	TGraph2D	TGraph2DErrors	（なし）

用途に応じて使うクラスを決めたあと、インスタンスを生成します。
このとき、定義するとき大きく分けて3つの方法があります。

1. データをあとから挿入していく方法
   g = TGraph()とまず定義したあとにg.SetPointで挿入していくことができます。この記事ではこの方法を採用しています。
2. データを配列で生成時に渡す方法
   g = TGraph(N, array_x, array_y)といった感じに配列で最初に渡せます。ベクトルでも構わないです。
3. データをファイルから渡す方法
   g = TGraph("data.txt", "%lg %lg", "")というようにデータファイルを指定してデータをグラフに詰めてくれます。自分でファイルを開いて、while文などを使ってループして...というコーディングが不要なので、めちゃ便利。

コンストラクタの詳細は表内のリンクを参照してください。

グラフを描く

早速以下のサンプルを使ってグラフを描いてみましょう。

test3-1.py

from ROOT import TGraph, TRandom, TCanvas

#
rdm = TRandom()
rdm.SetSeed(1)

#
c = TCanvas()

# グラフのインスタンスを生成
g = TGraph()
# グラフの表示スタイルなどを指定
g.SetMarkerStyle(20)
g.SetMarkerSize(1.1)
g.SetTitle("graph;Variable X;Variable Y")

# データを詰める
for x in range(10):
    y = 0.6 * x + 1.0
    g.SetPoint( g.GetN(), x, rdm.Gaus(y, 0.1*y) )

c.cd()
g.Draw("AP")
c.SaveAs("test3-1.png")

インスタンスを生成したあとに、g.SetMarkerStyle(20)などグラフの表示スタイルを変更しています。これは、デフォルトの表示スタイルだと、グラフのデータ点（マーカー）が小さすぎて見えにくいためです。
g.Draw("AP")で、グラフを描画しています。"AP"は表示オプションで、今回は軸（A）とマーカーで表示（P）するようにしています。その他のオプションについてはこちらを参照してください。
このコードを実行すると次のような図が出力されます。

グラフをフィットする

データを表示したあと、ある関数でフィットしたいときがあります。
そのような時は、TF1というクラスを使ってフィットします。

test3-2.py

from ROOT import TGraph, TRandom, TCanvas, TF1

#
rdm = TRandom()
rdm.SetSeed(1)

#
c = TCanvas()

# フィット用の関数の定義 y = p0 + p1 * x
f = TF1("f", "[0] + [1] * x", 0, 10)
f.SetParameter(0, 2)
f.SetParameter(1, 2)

#
g = TGraph()
g.SetMarkerStyle(20)
g.SetMarkerSize(1.1)
g.SetTitle("graph;Variable X;Variable Y")

for x in range(10):
    y = 0.6 * x + 1.0
    g.SetPoint( g.GetN(), x, rdm.Gaus(y, 0.1*y) )

c.cd()
g.Draw("AP")
g.Fit(f, "", "", 0, 10)
c.SaveAs("test3-2.png")

このコードを実行すると、ターミナルに次のようなフィット結果が表示されます。

$ python test3-2.py 

****************************************
Minimizer is Minuit / Migrad
Chi2                      =      1.17618
NDf                       =            8
Edm                       =  8.11236e-23
NCalls                    =           32
p0                        =      1.08717   +/-   0.225366    
p1                        =     0.566031   +/-   0.0422149   
Info in <TCanvas::Print>: png file test3-2.png has been created

グラフでエラーバンドを描く

あるモデルがデータと合致しているかみたいときに、モデルの予想とそのエラー、そしてデータを同時にみたいときがあります。このようなときにTGraphを使うと便利です。
サンプルコードは以下の通り。

test3-3.py

import ROOT
from ROOT import TGraph, TRandom, TCanvas, TF1, TMath, TGraphErrors

#
rdm = TRandom()
rdm.SetSeed(1)

#
c = TCanvas()

# データ用のグラフ
g = TGraph()
g.SetMarkerStyle(20)
g.SetMarkerSize(1.1)
g.SetTitle("graph;Variable X;Variable Y")

# データを詰める
for x in range(10):
    y = TMath.Sin(0.5*x) + 0.6 * x + 1.0
    g.SetPoint( g.GetN(), x, rdm.Gaus(y, 0.1*y) )

# 表示フレーム用のヒストグラムを取得
frame = g.GetHistogram()

# フィット関数を定義 y = p0 + p1 * x
f = TF1("f", "[0] + [1] * x", 0, 10)
f.SetParameter(0, 2)
f.SetParameter(1, 2)
# フィット(N0のオプションで自動的にフィット結果を図に表示しないようにしている)
g.Fit(f, "N0", "", 0, 10)

# フィット後の関数をモデルとして、
# モデルの予想とエラーを表示するためのグラフを用意
model = TGraphErrors()
model_1s = TGraphErrors()
model.SetLineColor(ROOT.kBlack)
model_1s.SetFillColor(ROOT.kGray)

# モデル用のグラフに予想を詰める
# エラーは平方和の平方根でフィットの誤差を伝搬させて計算
for i in range(50):
        x = i * 10. / 50.
        y = f.Eval(x)
        delta_p0 = f.GetParError(0)
        delta_p1 = f.GetParError(1)
        y_err = TMath.Hypot( (x * delta_p1), delta_p0 )
        model.SetPoint( i, x, y)
        model_1s.SetPoint( i, x, y)
        model_1s.SetPointError( i, 0., y_err)
# 結果を描画
c.cd()
frame.Draw("AXIS")
model_1s.Draw("C3 SAME")
model.Draw("LSAME")
g.Draw("P SAME")
frame.Draw("AXIS SAME")
c.SaveAs("test3-3.png")

今回はモデル予想用のグラフにTGraphErrorを使用しました。
modelの方には誤差を入力しないので、TGraphでも構わないです。
modelとmodel_1sには、ある横軸の値(x)とその値に対してEval関数を使って算出した縦軸の値(y)を挿入しています。
model_1sにはフィット結果から計算した誤差をmodel_1s.SetPointError( i, 0., y_err)のように、各データ点に対して挿入しています。

描画部分は少し複雑に見えると思います。
まず、frame.Draw("AXIS")でグラフを描画するための外枠を描画しています。

model_1sの表示オプションである"C3"`によって、各データ点間をスムーズにつなぎ（C）、エラーをバンドとして表示（3）させています。今回は直線でフィットしているので、"L3"でも大丈夫だと思います。

ここでのポイントは、model_1s.Draw("C3 SAME")、model.Draw("LSAME")、g.Draw("P SAME")をこの順番で描画しています。
ROOTではコードの上から順番にオブジェクトを上書きしていくので、データがエラーバンドで見えなくなることを防ぐためにこうしています。
ちょっとやっかいなのが、この上書きが軸に対しても行われてしまいます。
つまり、エラーバンドが軸の上に描かれてしまうので、軸や軸のメモリがバンドで隠れて見えなくなってしまいます。
軸を見えるようにするために最後にframe.Draw("AXIS SAME")と軸だけ上書きしています。

このコードを実行すると、下のフィット結果と図が表示されます。

$ python test3-3.py

****************************************
Minimizer is Minuit / Migrad
Chi2                      =      3.04961
NDf                       =            8
Edm                       =  4.98288e-23
NCalls                    =           32
p0                        =      1.99379   +/-   0.362888    
p1                        =     0.406821   +/-   0.0679751   
Info in <TCanvas::Print>: png file test3-3.png has been created

注意点はいくつかありますが、このように綺麗な図が作れました。