PyROOT tips

# 環境

以下のtipsでは、

```python
import ROOT as R
```

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

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

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

# よく使うオプション一覧

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

# TTree を読む

```python
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を作りたい場合。

```python
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からヒストグラムを書く

- 一番シンプルな方法

```python
tree.Draw("pt")
```

- 自分の好きな様に定義したヒストグラムにリダイレクトする：

```python
hist = R.TH1D("hist_name", "title", 100, 0, 1)
tree.Draw("pt >> hist_name") # 「>> hist」 ではなく、ヒストグラムの名前を使用することに留意
```

- 条件を掛けてヒストグラムを作る

```python
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種類あって、

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

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

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

## 最大値・最小値

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

## エラーについて
ヒストグラムをスケールする際に注意しておく必要のある基本事項について。まず、以下の様な適当なヒストグラムを考える。適当に階段状にFillしたもの。
![before.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/137120/199f5b13-d9d9-12a7-f797-93797efa472f.png)
各ビンの中心に書かれた縦線はエラーバーを示していて、各ビンは

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

```python
hist.Scale(10) # 例えば10倍の補正を加える
```
とすると、次のような結果が得られる。
![after.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/137120/237e8231-1e1b-b386-89e4-e2e9440c80ab.png)
各ビンのエントリーは10倍になっていて狙い通りの補正が掛かっているが、各ビンのエラーバーは狙い通りの挙動を示していない。例えば右から2番目のビンはエントリー数が9→90となってエントリー自体は10倍されているが、エラーバーの大きさ10倍されている。つまり本来は

```math
10 \times N_{bin} \pm \sqrt{10 \times N_{bin}}
```
という振る舞いをしてほしいのに、単にスケールした場合は

```math
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$番目のイベントを示す）。

```math 
\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.

基本的に各ビンの

```math 
\sqrt{\sum_{j=0}^{N}w_{ij}^2}
```
を返すと考えておけば良い（weight=1のときには、$\sqrt{N}$と等しくなるだけ）。

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

```python
hist.Rebin(10)
```

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

```python

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`を使わずそのままのリストを渡すと、

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

って出ます。

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

```python
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)
```


## ビンの情報を取る

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


# TCanvas

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

```python
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`から触るのではなくて

```python
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/

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

```shell
>>> 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に必要なコマンドを書いて終わり

```shell
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

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

- https://sft.its.cern.ch/jira/si/jira.issueviews:issue-html/ROOT-9492/ROOT-9492.html

## root-6.18.04

```shell
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` が上書きされてしまったので再コンパイルする必要があった
- いざ再コンパイルすると、

```shell
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 からパッチの当たったバージョンを持ってきた。で、以下のコマンドを実行。

```shell
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にパスも通った。このときには

```shell
/usr/local/bin/python3 
```

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

## その他の手続き
### .bashrc に書くべき内容
ビルドディレクトリの名前は各自変更のこと。

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