環境
以下の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したもの。
各ビンの中心に書かれた縦線はエラーバーを示していて、各ビンは
N_{bin} \pm \sqrt{N_{bin}}
のふらつきを持って分布していると考える(エラーの大きさの$\sqrt{N}$はポワソン分布から出てくる)。例えばこれに対してoverallに何らかの補正を加える場合を考える(全体に対して何らかのweightを掛け忘れていた、とか)。そのときに何も考えず、
hist.Scale(10) # 例えば10倍の補正を加える
とすると、次のような結果が得られる。
各ビンのエントリーは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]の範囲でビン数100のヒストグラムがあるとする(=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