-- ShionChen - 2015-01-05

参考にしたもの:

* KamonoWiki: 京都の人のtipまとめ。かなりオススメ

* ATLAS japan root tutorial: 田中先生が作ったマニュアル。ちょっと古いがまとまっててよい。

Quick draw using the interactive mode of root

TTreeから情報を抜き取って結果(ヒストグラムなど)にする方法はおおまかに4つ

MakeClass してひな形を作ってMakefile書いてまじめにコンパイル

② C++のコードは書くけどg++じゃなくてrootのコンパイラ(rootCINT)でコンパイル

③ C++っぽいマクロ書いてrootにそのまま食わせる

④ rootのコマンドライン ("interactive mode") で強行

上に行くほど格調高いがめんどい。下に行くほど難しいことはやりづらいがお手軽といった感じ。

ここでは④ (と最後に③) の技を紹介します。

基本的にDrawをすればよいのですが, オプションとの合わせ技によって意外と複雑なことができます。



1次元ヒストを極めよう

例が欲しいので, まずは適当なsmall treeを開いて下さい。

$ cd /home/chen/gpfs/samples/BHole_slimmed
$ root -l slimmed_BM_c1_n4_Mth6000_MD3000_13TeV.root  
 
第1引数で描きたい変数を指定

- 対象としているbranchの型がただの数 (int , float ,double) の場合, Drawするとただ毎イベントの値がFillされる

e.g. 各イベントのjet数

   $ smtruth -> Draw("jet_n") 

- classがbranchに入ってる場合はそのメンバ・メソッドが普通にアクセスできる

e.g. truth BH mass

   $ smtruth -> Draw("tr_pBH.M()")


e.g. MET vectorの方向

   $ smtruth -> Draw("pMis_Int.Phi()")

この辺はまあよいでしょう。


- vectorや配列の場合は, 何も指定しないと全部の要素がFillされる。

e.g. (全イベントの) 各jetのpt

   $ smtruth -> Draw("v_jet.Pt()")


- .at(i)もしくは[i]をつけてやると各イベントのi番目の要素だけFillされる。

e.g. (全イベントの) 0番目のjetのpt

   $ smtruth -> Draw("v_jet[0].Pt()")

 ※small treeではobjectはpTの大きい順に入ってるので, この場合はイベント中でpTが最大のjet (leading jet) を見てることになる.


- 演算もok. 普通の算数はもちろん, branch間の演算も

e.g. 再構成したBH mass と truth BH mass の差(再構成精度)

    $ smtruth -> Draw("pvisObj_tot.M() - tr_pBH.M()")


- tlv同士の足し算もいける

e.g. pTの最も大きい2 jetで組んだmass

    $ smtruth -> Draw("(v_jet[0]+v_jet[1]).M()")

 ※ " "で囲まれた中では基本的にpythonの書式が有効.

  c++の関数は一部使えない(例えばfabs()はダメ。abs()しか使えない)



第2引数: cut

- 第2引数に条件式を書き込むと, それが真となるイベントだけFillされる

e.g. MET>50GeVのイベントのleading jet pT分布

   $ smtruth -> Draw("v_jet[0]", "pMis_Int.Pt()>50");

しかも頭いいことにこいつはイベント単位の判定(条件を満たしたイベントは残して, そうじゃないイベントは捨てるというタイプの判定)だけじゃなくて, もっと細かい単位でcutもかけられる

e.g. pT>100GeVのjetの eta分布

   $ smtruth -> Draw("v_jet.Eta()", "v_jet.Pt>100");

 ※この書式は, 第一引数にいる変数と第二引数にいる変数で, 毎イベントvector(or 配列)のサイズが同じときのみ有効

  今の例でいうとv_jet.Eta()とv_jet.Pt()の要素数は毎イベント同じなので大丈夫



第2引数-2: weight

- ヒストグラムの縦のスケールを変えたい or イベントごとにFillのweightを変えたい場合は対応するweightを第二引数に書いてやればよい

e.g. electronの再構成効率を考慮した electron pT分布

    $ smtruth -> Draw("v_el.Pt()","el_reco_eff")

e.g. ヒストグラムをcross-sectionとluminosity (10fb-1) で規格化

    $ smtruth -> Draw("v_el.Pt()","xsec*10/nMC")

  ※ この場合はTH1の"Scale()"でヒストグラムを縦軸方向に丸ごと拡大・縮小してやってもよい

   cut optionと同時に使いたい場合は, 条件文をカッコで囲んで, weightはかけ算してあげればok.

e.g.

    $ smtruth -> Draw("jet_n","(pMis_Int.Pt()>50)*el_reco_eff")



第3引数: draw option

- "same"で重ね書き

e.g. 無回転BHと回転してるBH (c3) でjetの角度分布を比べてみる

   $ root -l
   $ TFile *_file0 = new TFile("slimmed_BM_c1_n4_Mth6000_MD3000_13TeV.root");
   $ TFile *_file1 = new TFile("slimmed_BM_c3_n4_Mth6000_MD3000_13TeV.root");

// 2つのTFileを読み込む. 単純に $ root -l slimmed_BM_c*_n4_Mth6000_MD3000_13TeV.root と打つだけでも同じことが起こる

   $ TTree *t0 = (TTree*) _file0 -> Get("smtruth");
   $ TTree *t1 = (TTree*) _file1 -> Get("smtruth");

// TFileからTTreeを読み込む

   $t0 -> Draw("v_jet.Eta()");
   $t1 -> Draw("v_jet.Eta()", "", "same");

  ※"same"を"sames"にすると統計boxも2つ出てくる

- "E": ヒストにエラーバーつける

 ※ デフォルトだとPoissonエラー。エラーを自分設定したい場合はSetBinError()で1binずつ変えるしかない。

 ※ weightかけたりヒスト全体をスケールしたりするとエラー計算が複雑になるときがある。entryに対する相対的なエラーの大きさを固定したい場合は h -> Sumw2() というコマンドを最初に打っとくこと。



ヒストの設定(色とか)を変えたい → TH1インスタンスを取得しよう

第一引数の最後に >>をつけるとヒストをTH1のインスタンスとして引っぱってこれる.

e.g.

    $ smtruth -> Draw("jet_n>>h");

以降このインスタンス "h" はコマンドラインで有効. 普通のc++のように関数呼んで色々できる。

e.g.

 ヒストの線の色を赤にする

    $ h -> SetLineColor(2)   

ヒストを斜線で塗りつぶす

    $ h -> SetFillStyle(3004)

ヒストの塗りつぶしの色を青にする

    $ h -> SetFillColor(4)

ヒストの上(下/左/右)の余白を調節(0.3: 余白のスペースの全体に対する割合)

   $ h -> SetHistTop (Bottom/Left/Right) Margin(0.3) 

ヒストのx軸の範囲を1~5に変更

   $ h -> GetXaxis() -> SetRangeUser(1,5);

などなど

 ※ 描画設定に関する関数は基本的に Set + ("Line" / "Marker" / "Fill") + ("Size" / "Color" / "Style" / "Width") みたいな構造.

 ※ root color/style code: https://root.cern.ch/root/html/TAttFill.html

 ※ TH1でできる設定: http://root.cern.ch/root/html/TH1.html, "Set~"で始まる関数が大体そうです。

 ※ "h"を宣言する際にrange & binning指定したい場合

     $ smtruth -> Draw("jet_n>>h(10, 0, 10)");

  (nbins, xmin, xmaxの順)



DrawしないでTH1のインスタンスだけ取ってきたい

最初にTH1を宣言しとく

   $ TH1F *h = new TH1F("h","h",10,0,10);

TTree::Projectを使う. (e.g. MET>50GeVのイベントのjet_n分布)

   $ smtruth -> Project("h", jet_n, "pMis_Int.Pt()>50")

※ Projectの第2・第3引数は, Drawの第1・第2引数にそれぞれ完全対応



ヒストグラムの統計情報を取得したい

描いたヒストのMeanとかRMSは統計boxに表示されますが, プログラムにそのまま渡したくなった場合は

   $ h -> GetMean()
   $ h -> GetRMS()

とかで返り値として持ってこれます。他にもGetMeanError (=RMS/sqrt(n)) や GetRMSError, GetOverflow など色々あります。詳しくは TH1のページを参照。

ちなみに統計boxを消したい場合はTStyleクラスのグローバル変数 "gStyle"を使って

   $gStyle -> SetOptStat(000000);

逆にoverflow/underflowなどの情報まで表示させたい場合は

   $gStyle -> SetOptStat(111111);

※ デフォルトのflagは"111100"です



TTreeの組み込み変数・関数たち(Sum$は超便利!!)

rootのinteractive modeにはいくつか便利な変数・関数が最初から用意されてる。例えば

Entry$: TTreeのentry index

e.g. Treeの3イベント目のjet PtをFill

   $ smtruth -> Draw("v_jet.Pt()","Entry$==3")

Sum$(): 引数はvector(と条件式)。引数の取り方で挙動が変わる

① Sum$(vector A): 各イベントにおけるAの全要素の和を返す

e.g. 各イベントにおけるpT 100GeV 以上の jetのpT sum

     $ smtruth -> Draw("Sum$(v_jet.Pt())","v_jet.Pt()>100")

② Sum$(vector Aに関する条件式): 条件式を満たすvector Aの要素の数を返す

e.g. 各イベントのpT 100GeV 以上の jetの数

     $ smtruth -> Draw("Sum$(v_jet.Pt()>100)")

他にもMaxIf$とかMinIf$とか便利なものがたくさんあるので, 興味があったら https://root.cern.ch/root/html532/TTree.html の真ん中らへんを見てみてください。

「こんな機能あったら便利なんだけど」って思ったやつは大概あります。



2Dヒストを書こう

基本的に1Dと同じ。

第1引数で描きたい2変数を指定

Drawの第一変数(Projectの第二変数)で変数の間にコロンを打つだけ。

e.g. BH massの再構成精度とMETの大きさの2Dヒスト (cut: jet_n>=4)

   $ smtruth -> Draw("tr_pBH.M()-tr_pBH.M():pMis_Int.Pt()","jet_n>=4")

なぜかコロンの前の変数がy軸, コロンの後の変数がx軸に描かれます。



TH2インスタンスを持ってくる

1Dのときと同じく第1引数の最後に>>つければok. bin数と範囲指定したい場合は以下の通り.(e.g. tr_pBH.M()-tr_pBH.M(): 100bin, -1000GeV~1000GeV, pMis_Int.Pt(): 100bin, 0GeV~5000GeV)

   $ smtruth -> Draw("tr_pBH.M()-tr_pBH.M():pMis_Int.Pt() >> h(100,0,5000, 100, -1000,1000)","jet_n>=4")

h(x_nbins, xmin, xmax, y_nbins, ymin, ymax)という順序です。第1引数のx-y notaionと逆なのでよく混乱します。

第2引数 - cut & weight

1Dヒストのときと全く同じ



第3引数 - Draw option

デフォルトだと黒の汚い散布図ですが, 点の数が増えてくるとかなりしんどいので, 色でentryの密度を表す"colz" (color-z)オプションをオススメします。

実質こいつはz軸みたいな扱いなので, 色の範囲を指定したかったら1Dのときと同様に ->GetZaxis() -> SetRangeUser (zmin, zmax) みたいにすればokです。

他にも"lego"とか"box"とか意味わかんないオプションがいっぱいあるので興味あったら こちらを.

"same"系は残念ながら使えませんが, 後述の「画面分割」やヒストグラム同士の除法( link, page 21)などで直接比較はできます。



画面分割

TCanvasというクラスのインスタンスをまず宣言しよう。

   $ TCanvas *c = new TCanvas("c","c");

"c"というタイトルがついたキャンバスが出てくると思います。

※ interactive modeの場合, 宣言しなくてもDrawをすると, デフォルトで"c1"というTCanvasのインスタンスが勝手に作られる。

※ キャンバスの寸法は宣言時にいじれます

e.g. 四隅の座標をピクセルで指定

       $ TCanvas *c = new TCanvas("c","c", xmin, ymin, xmax, ymax);

e.g. 横×縦の大きさを指定

      $ TCanvas *c = new TCanvas("c","c", xwidth, ywidth);

画面分割は TCanvas::Divide(columns, rows)を使います.

$ c -> Divide(2,1);

※ この場合, 横2×縦1になる.

分割した画面 ("pad")にはTCanvas::cd() でアクセスできます。

e.g. 2つの2Dヒストを同じキャンバスに描く

   $ c -> cd(1);
$ smtruth -> Draw("tr_pBH.M()-tr_pBH.M():pMis_Int.Pt()","jet_n>=4")
   $ c -> cd(2);
$ smtruth -> Draw("tr_pBH.M()-tr_pBH.M():pMis_Int.Pt()","jet_n<4")

※ padの番号は "左上→右上→左下→右下" の順番

※ 絵の保存はTCanvas::Print(char** filename, char** filetype)です.

e.g. c -> Print("home.eps","eps"); // eps形式で保存




マクロを書こう(冒頭の③の方法)

まじめにやろうとすると1枚の絵を描くのに数十行のコマンドを打ち込むハメになることもざらにありますが、一行ずつ手打ちだとやっぱしんどくなってくるときがあります。(保存してなかった絵をもう一回描きたい, ループ回して同じようなスタイルの絵をいっぱい描きたいetc.)

interactive modeではこれらのコマンドをただ羅列しただけのソースファイルを普通にマクロとして読んでくれるので、ちゃんとした絵を描きたくなったときはマクロが非常に強力です。

書式はクソ簡単。

{

(command1);

(command2);

(command3);

....

}

各コマンドにセミコロンつけて、中カッコで囲むだけ。関数名もヘッダーもreturnもいらん。ファイルの名前は"~.C"とか"~.cc"とかにしとけばok。

※ 基本的にC/C++の関数・書式は全部通用する.(ただしコンパイラを通してないので若干仕様が変. 多少バグってても空気読んで勝手に解釈して動く. なのであんま複雑なC/C++コードは実装しない方が無難です. )

後はrootを起動するときに

$ root -l hoge.C

のように呼ぶか, interactive mode に入った後に

$ .x hoge.C

みたいな感じで呼んであげればokです。

Styleを定義するマクロ

Drawのデフォルトの設定(LineColor, FillStyle, タイトルの大きさ・書式 etc.)はTH1とかを通じてじゃなくても, TStyleというクラスを通じて設定できます。

なので多用する設定はマクロにまとめて、起動時に(or 他のマクロの中から)呼ぶようにすると効率よく書式をいじれます。

例をいくつか:

- Chen Style: (login icepp) /home/chen/rootlogon.C

- ATLAS Style ("root atlasっぽいplotを描く"):

実行のしかた

interactive modeのコマンドラインから実行する場合は上と同様

$ .x hoge.C

他のマクロから呼ぶ場合は, そのマクロの先頭かなんかに

gROOT -> ProcessLine(".x hoge.C");

という行を加えればokです.

Edit | Attach | Watch | Print version | History: r3 < r2 < r1 | Backlinks | Raw View | WYSIWYG | More topic actions
Topic revision: r3 - 2015-08-11 - ShionChen
 
    • Cern Search Icon Cern Search
    • TWiki Search Icon TWiki Search
    • Google Search Icon Google Search

    Main All webs login

This site is powered by the TWiki collaboration platform Powered by PerlCopyright &© 2008-2024 by the contributing authors. All material on this collaboration platform is the property of the contributing authors.
or Ideas, requests, problems regarding TWiki? use Discourse or Send feedback