LoginSignup
0
3

More than 3 years have passed since last update.

@kabupen

PyROOT tips

Last updated at Posted at 2019-11-05

環境

以下のtipsでは、

import ROOT as R

の宣言をしていることを前提にしています。また適宜 tree、histだったり変数が初期化なしに使用している箇所があると思いますが、

tree = input_file.Get("my_tree")
hist = input_file.Get("hist")
hist = R.TH1D(...)

とかをしている前提なんだなと思っておいてください（あまりにもspecificな変数であればその都度初期化も明記します）。

よく使うオプション一覧

R.gStyle.SetOptStat(0) # 統計ボックスを消す

TTree を読む

input_file = R.TFile("input.root")
tree = input_file.Get("mytree")

for e in tree:
    print(e.pt) # 例えば pt というブランチがあれば、これでアクセスできる

TTree から TTreeを書き出す

大元のntupleをskimして（slimして）、新しいTTreeを作りたい場合。

input_file  = R.TFile("orig.root")
output_file = R.TFile("new.root", "recreate")

for tree in input_file.GetListOfKeys():                                                                                                                                                                                                                       

    tree = tree.ReadObj()

    R.gROOT.ProcessLine('\
            struct MyStruct{\
            Double_t pT;\
            };')
    struct = R.MyStruct()

    new_tree = R.TTree(tree.GetName(), tree.GetName())
    new_tree.Branch("pT", R.AddressOf(struct, "pT"), "pT/D")

    for e in tree:
        if e.pT < 100 : continue

        struct.pT = p.sT
        new_tree.Fill()

    new_tree.Write()

上記のスクリプトを使うことで、入力ファイルのTTreeからとある条件を掛けた場合（例えば運動量pTが100GeV以上のものを取り出す、とか）のエントリーだけを新しいTTrreとして定義してnew.rootに書き出すことができる。この他にも（試してはいないけど）

  • CopyTree で skimする

という方法もあるので、もしかしたらそちらのほうが簡単かも？検討ください。

TTreeからヒストグラムを書く

  • 一番シンプルな方法
tree.Draw("pt")
  • 自分の好きな様に定義したヒストグラムにリダイレクトする：
hist = R.TH1D("hist_name", "title", 100, 0, 1)
tree.Draw("pt >> hist_name") # 「>> hist」 ではなく、ヒストグラムの名前を使用することに留意
  • 条件を掛けてヒストグラムを作る
tree.Draw("pt", "pt > 100")
tree.Draw("pt", "eta > 1")  # 別ブランチの変数も条件に使用できる
tree.Draw("pt", "eta > 1 && pt > 100 || energy > 100")  # AND や OR も取れる

TH1, TH2

ヒストグラムについて

Fillには2種類あって、

hist         .Fill(x)
hist_weighted.Fill(x, w)

Fillを呼び出した回数がそのままイベント数になる場合と、そうはならない場合。前者は実データをFillするときに使用する。後者の例として（例えば）シミュレーションサンプルを使う場合、実データを再現するように何らかの補正（weight, w）を掛けるために使用する。hist.Fill(x)も weight 1の場合と考えると良い。weight付きヒストグラムを扱う際には、注意が必要な点がいくつかある。

GetEntries()        # エントリー数を返す。単にFillを呼び出した回数に相当する
GetBinContent(ibin) # i番目のビンの sum of weight を返す

最大値・最小値

hist.GetMaximum() # y軸最大値
hist.GetMinimum()　# y軸最大

エラーについて

ヒストグラムをスケールする際に注意しておく必要のある基本事項について。まず、以下の様な適当なヒストグラムを考える。適当に階段状にFillしたもの。
before.png
各ビンの中心に書かれた縦線はエラーバーを示していて、各ビンは

N_{bin} \pm \sqrt{N_{bin}}

のふらつきを持って分布していると考える（エラーの大きさの$\sqrt{N}$はポワソン分布から出てくる）。例えばこれに対してoverallに何らかの補正を加える場合を考える（全体に対して何らかのweightを掛け忘れていた、とか）。そのときに何も考えず、

hist.Scale(10) # 例えば10倍の補正を加える

とすると、次のような結果が得られる。
after.png
各ビンのエントリーは10倍になっていて狙い通りの補正が掛かっているが、各ビンのエラーバーは狙い通りの挙動を示していない。例えば右から2番目のビンはエントリー数が9→90となってエントリー自体は10倍されているが、エラーバーの大きさ10倍されている。つまり本来は

10 \times N_{bin} \pm \sqrt{10 \times N_{bin}}

という振る舞いをしてほしいのに、単にスケールした場合は

10 \times \left( N_{bin} \pm \sqrt{N_{bin}} \right)

とエラーの大きさも含めた全体が10倍になる挙動となってしまう。そのため今注目したビンでは、$\sqrt{9}\times 10=30$の大きさのエラーが付いていることになる。これを防ぐには、Scaleする前にSumw2()を呼んでおきましょう、という基本的な事項。

Sumw2

Create structure to store sum of squares of weights.
If histogram is already filled, the sum of squares of weights is filled with the existing bin contents.
The error per bin will be computed as sqrt(sum of squares of weight) for each bin.

i番目のビンにおける sumw2 とは、そのビンにFillされるイベント（$j=0〜N$）に対して掛けられるweightの二乗を足し合わせたものを表す（下式の添字$ij$は$i$番目ビンの$j$番目のイベントを示す）。

\sum_{j=0}^{N}w_{ij}^2

この定義式から、weightを与えずにFillした場合($w=1$)のsumw2はそのビンのエントリー数と等しいことが分かる。

GetBinError

If the sum of squares of weights has been defined (via Sumw2), this function returns the sqrt(sum of w2). otherwise it returns the sqrt(contents) for this bin.

基本的に各ビンの

\sqrt{\sum_{j=0}^{N}w_{ij}^2}

を返すと考えておけば良い（weight=1のときには、$\sqrt{N}$と等しくなるだけ）。

ビン幅の調整

今、[0,1]の範囲でビン数１００のヒストグラムがあるとする（＝0から1を0.01刻みで分割したヒストグラム）。これを0.1刻み（より粗いビンニングに変更する）に変更する場合は、

hist.Rebin(10)

とすればよい。100個のビンを等間隔に10個まとめて1ビンとする操作を意味する。また、ユーザー定義のビン幅、特に非等幅なビンニングにするときには以下の操作が有効（ちなみにimport arrayを使っているのはRebinの第三引数が配列の先頭ポインタを取ることに由来します。これが無いと実行時クラッシュします）。同じRebin関数を使って、


bins = [0, 0.1, 0.5, 1] # 各ビンの左端のxの値を指定

import array
bins = array.array('d', bins)
hist.Rebin(len(rebins)-1, "new_hist_name", rebins)

とすれば、ユーザー定義のビン間隔で簡単に調整することができる。ちなみに、arrayを使わずそのままのリストを渡すと、

TypeError: TH1* TH1::Rebin(int ngroup = 2, const char* newname = "", const double* xbins = 0) =>
    could not convert argument 3

って出ます。

TH2の場合

TH2にはTH1と違ってヒストグラムを作成してから適当にrebinするための関数が無い（等幅ビンであれば用意されている）ので、コンストラクタの段階で可変ビンとして宣言しておく必要がある。

rebins_x = array.array('d', [0, 10, 20, 100])   # テキトーな例
rebins_y = array.array('d', [0, 50, 100, 1000]) # テキトーな例

hist = R.TH2D("","", len(rebins_x)-1, rebins_x, len(rebins_y)-1, rebins_y)

ビンの情報を取る

hist.GetBinContent(ibin) # i-th bin のエントリー数
hist.GetBinLowEdge(ibin) # i-th bin の左端のx軸

TCanvas

軸の範囲は様々なメソッドから変更できますが、下記のやり方が一番簡単かと思います（もちろん慣れている方法が一番ですが）。

canvas = R.TCanvas("","")
canvas.DrawFrame(0,0,1,1) # x軸[0,1]、y軸[0,1]のフレームを描画する
hist.Draw("same")         # same オプションが必要

こうすることで、たとえhistのx軸の範囲が0~100とかであれ、ビニングを変更することなく [0.1]の範囲だけを描画することができる。ただしコメントにも示したとおり、sameオプションを付けておかないとDrawを呼んだときに新たに軸が描画されてしまうので DrawFrameで書いたフレームは上書きされて消えます。

これはDrawFrameの返り値がTH1Fであることに由来します。DrawFrameを呼ぶことで空のヒストグラムを書いていて、そのキャンバスに自分のヒストグラムを書くという流れです。また、軸のタイトルを触りたいときはhistから触るのではなくて

frame = canvas.DrawFrame(0,0,1,1) # 返り値はTH1Fなので「フレーム」ではないですが、気分はフレームのつもりです
frame.GetXaxis().SetTitle("好きなように")
frame.GetYaxis().SetTitle("好きなように")

といった形で触るのがベストです。

軸がかぶった

canvas.RedrawAxis()

THStack

凡例（レジェンド）

TLegend

legend = R.TLegend(x_min, y_min, x_max, y_max)
legend.SetBorderSize(0) // 枠線で囲まない
legend.SetFillStyle(0) // デフォルト色で背景を塗りつぶさない
legend.AddEntry( hist, "legend title", "l") // lオプション：線、p：ポイント、f：四角

TLatex

legend = R.TLatex()
legend.SetNDC(1) 
legend.DrawLatex(0.1,0.5, "p_T")

ROOTのインストール (mac)

URL

基本知識

https://root.cern/install/
そもそも ROOT、PyROOTをローカルPCにインストールするためのレシピ。
pythonのバージョンを上げたりすると、その都度コンパイルが必要となるので注意。

>>> import ROOT as R
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/usr/local/root/root-6.18.04/cmk_py3/lib/ROOT.py", line 24, in <module>
    import cppyy
  File "/usr/local/root/root-6.18.04/cmk_py3/lib/cppyy.py", line 61, in <module>
    import libPyROOT as _backend
ImportError: dlopen(/usr/local/root/root-6.18.04/cmk_py3/lib/libPyROOT.so, 2): Library not loaded: /usr/local/opt/python/Frameworks/Python.framework/Versions/3.7/Python
  Referenced from: /usr/local/root/root-6.18.04/cmk_py3/lib/libPyROOT.so
  Reason: image not found

のエラーはpythonのバージョンを上げただけで、再コンパイルしていなかったりすると出る。
再コンパイルする際には、新しくcmakeディレクトリを作り直すのが一番簡単。

共通のこと

  • ROOTのソースコードをダウンロードしてくる
  • バージョンを複数準備したいなら、それに対応する（分かりやすい任意の名前の）ディレクトリを用意する
  • それぞれ cmake ビルド用のディレクトリを用意して、そこで各々ビルドする
  • ビルドが終わったら、bashrcに必要なコマンドを書いて終わり
mkdir /usr/local/root/root-6.06.02/
cd /usr/local/root/root-6.06.02/
mkdir cmake_build_python3
cd cmake_build_python3
cmake .. -Dall=ON -DPYTHON_EXECUTABLE=/usr/local/bin/python3 -Dsoversion="ON"

root-6.06.02

cmake .. -Dall=ON -DPYTHON_EXECUTABLE=/usr/local/bin/python3 -Dsoversion="ON"
make -j10

root-6.18.04

mkdir cmake_python3.8
cd cmake_python3.8
cmake .. -Dall=ON -DPYTHON_EXECUTABLE=/usr/local/bin/python3 -Dbuiltin_xrootd=ON
make -j10
  • pythonのバージョン毎にディレクトリを作り置きしている

root-6.22.04

  • 別件でbrew install imagemagickを実行したら、/usr/local/bin/python3 が上書きされてしまったので再コンパイルする必要があった
  • いざ再コンパイルすると、
error:
downloading 'http://grid-deployment.web.cern.ch/grid-deployment/dms/lcgutil/tar/davix/0.6.7/davix-embedded-0.6.7.tar.gz' failed
status_code: 22
status_string: "HTTP response code said error"

っていうエラーが出たので、ひとまず ROOT のバージョンを上げてみた（root-6.18.14->root-6.22.04）

  • すると、error: "std::basic_string_view" is ambiguous というエラーが出てcompileがコケたので GitHub からパッチの当たったバージョンを持ってきた。で、以下のコマンドを実行。
git clone https://github.com/root-project/root.git
git checkout -b v6-22-00-patches remotes/origin/v6-22-00-patches
mkdir build; cd build
cmake ..; make

でコンパイルが通り、ROOTにパスも通った。このときには

/usr/local/bin/python3

からでないと、自分の環境ではimport ROOTできなかったので注意。

その他の手続き

.bashrc に書くべき内容

ビルドディレクトリの名前は各自変更のこと。

cd /usr/local/root/root-6.06.02/cmake_build
source bin/thisroot.sh
cd - > /dev/null
0
3
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
0
3