-- Main.RyuSawada - 2016-12-06 ---+!! <nop>%TOPIC% %TOC% <!-- this line is optional --> %STARTINCLUDE% ---+ イントロダクション ROOTとは簡単にいうと、データ解析をするための道具集。 ROOTの中に、多くのライブラリ(用途ごとにクラスをまとめたもの)があり、これらのクラスを呼び出して自分の解析で使う。 ROOTのクラス名は普通大文字のTから始まるので、そのような名前を見たらそれはクラス名だと思って大丈夫です。 使い方は大きく分けて二つあり、 * ROOTのパッケージに含まれるインタープリター(Cling)上で使う * 自分のマクロやプログラムの中でROOTのクラスを使う どちらでも基本的な使い方は同じです。 Clingは簡単にイベントの分布を見たいときなどに使います。もう少し込み入った解析をしたかったり、同じ手続きを再利用したいときにはマクロを使い、もっと大掛かりな解析(eventの再構成や、物理フィット等)の場合はC++のプログラムにします。この講習では簡単のためClingを主に使いますが、同じ文法でマクロやC++でも解析に行うことができます。 ---++ データ解析でよく使うクラス データ解析の際に行う処理は単純にいうと、データを読み込み、何らかの計算を行い、選択してからデータ構造として格納するプロセスの繰り返しと言えます。 最終結果は分岐比や質量などの物理量や図として論文に載せます。 データ保存時のデータ構造としてROOTでは ツリー( =TTree= )を使います。 また、図を書く際は基本的には、ヒストグラム(TH1, TH2)かグラフ(TGraph, TGraphErrors)を使います。 ATLASでは生データをセントラルプロセスでDAODという形式のROOTファイルを作り、イベント再構成の結果(観測されたジェットや、レプトン、横消失エネルギー等)が書かれます。各実験グループがここから必要な情報を抜き出し、ツリーに詰め直したファイルを作ります。このファイルをntupleと呼んでいます。皆さんがATLASの解析を始めて最初に扱うのはこのntupleになる場合が多いと思います。 ---++ 環境設定 今回の講習はlogin.icepp.jpで行います。sshでloginしてからsetupATLASとlocalsetupROOTを実行してください。 ATLASではlocalSetupROOT等でROOTの使用に必要な環境設定が自動的にできますが、もしATLASのソフトウェアがない環境でROOTを使いたい場合は、ROOTSYS, PATH, LD_LIBRARY_PATH (MacではDYLD_LIBRARY_PATH)の三つを設定してください。ROOTSYSはROOTがインストールされているディレクトリ、PATHは通常は$ROOTSYS/bin/を追加、LD_LIBRARY_PATHには$ROOTSYS/libを追加します。 ---++ インタラクティブセッションの起動と終了 ROOTのインタラクティブセッション(Cling)は通常は以下のように =root= とタイプして始めます。 以下の例にあるように、ほとんどのC++のコードが実行可能です。またROOTに中のクラスや定数は#includeすることなく使えます。 それに加えて累乗のような拡張機能とClingの特殊コマンドがあります。 特殊コマンドはドットで始まります。特殊コマンドの詳細は.helpで見られます。 よく使うのは、メモリ上にあるオブジェクトのリストを表示する.lsと、セッションを終了する.qです。 通常は.qで終了しますが、プログラムがクラッシュした時など.qで終了しない場合があります。その場合は.qqqや.qqqqq等のようにqの数を増やすことで強制終了できます。 %SYNTAX{lang="c++" numbered="off"}% $ root root [] cout<<"Hello world"<<endl; root [] 3 * 2 root [] cos(TMath::Pi()) root [] 3^2 root [] for (int i = 0; i < 10; i++) { root (cont'ed, cancel with .@) [] cout<<i<<endl; root (cont'ed, cancel with .@) [] } root [] TTree tmp("tmp", "tmp") root [] .ls root [] .help root [] .q %ENDSYNTAX% ---++ 便利な機能 Clingにはヒストリー機能があるので、矢印やEmacs風のキーバインドでコマンドの履歴を遡れます。その上でreturnキーを押すとそのコマンドが実行されますが、その代わりにCTRL+oで実行するとその後に履歴の次のコマンドが表示されます。複数のコマンドを繰り返す場合にこの機能を使うと便利です。 また、コマンドの補完機能もあります。例えば、 %SYNTAX{lang="c++" numbered="off"}% root [] TGraphEr %ENDSYNTAX% の後でTabを打つと、Clingが知っているクラス名で補完してくれます。 さらに、下のような状態でTabを打つと、 %SYNTAX{lang="c++" numbered="off"}% root [] TGraphErrors a( %ENDSYNTAX% このクラスのコンストラクタの引数を知ることができます。 ---+ ヒストグラムとグラフ ヒストグラムとは一次元以上の変数の範囲を複数に区切ったもので、各領域(ビン)のイベント数が格納されます。 物理解析ではイベント数だけではなく誤差も重要なので、ROOTのヒストグラムクラスでは誤差も保存されます。 それに対して、グラフは二次元上の位置を持った点の配列で二つの変数の空間でのイベントの分布や、横軸に時間をとってある変数の時間変化のプロットを作る際等によく使われます。 ヒストグラムの例 %ATTACHURL%/histogram.png グラフの例 %ATTACHURL%/graph.png ---++ ヒストグラム ROOTの基本的なヒストグラムのクラスはTH1F, TH1S, TH2D等の名前がついています。3文字目が次元の数で通常は1か2です。 4文字目はビンの中身(普通はイベント数)の型を表していて、S: short, I: int, F:float, D:double 等があります。イベントに重みをつけたり、全体をスケールする時に浮動小数点を使う場合もあるので、通常はFを使っておくのがいいと思います。Dだと精度は高くなりますが、そこまでの精度が必要な場合は多くありません。Dにするとメモリー消費が増えます。 それでは一次元のヒストグラムを作ってみます。 %SYNTAX{lang="c++" numbered="off"}% $ root root [] TH1F hist("hist", "Histogram;x;number of events", 100, -10, 10) root [] for (int i = 0; i < 10000; i++) { hist.Fill(gRandom->Gaus(0, 7)); } root [] TCanvas c1("c1", "c1") root [] hist.Draw() %ENDSYNTAX% ここでは、一次元のヒストグラムを作り、ガウス分布で乱数を作ってヒストグラムに格納しました。 TH1Fの引数は順に、このオブジェクトの名前、タイトル、ビンの数、範囲の下限と上限です。 タイトルの中で =;= で区切られた二番目はx軸、三番目はy軸のタイトルとして使われます。 gRandomはROOTの中でグローバルに定義されているオブジェクトで乱数の生成に使われます。 (他にもgROOT, gSystem, gPad等のグローバルな変数があり、小文字のgで始まる慣習になっています。) ヒストグラムを表示するためのキャンバスをc1という名前で作成し、そのあと実際にヒストグラムを表示しました。 ヒストグラムの情報を得てみます。 %SYNTAX{lang="c++" numbered="off"}% root [] hist.GetEntries() root [] hist.GetMean() root [] hist.GetRMS() root [] hist.GetBinContent(50) root [] hist.GetBinError(50) root [] hist.GetBinContent(0) root [] hist.GetBinContent(101) %ENDSYNTAX% ここでは、histが持っている情報(エントリー数、平均、標準偏差)をとりました。 その次に例として50番目のビンの高さとそのエラーを取っています。ROOTでは何も指定しなければ、高さの二乗根がエラーとなります。 それで正しくない場合は自分で設定できます。 ヒストグラムの最初と最後のビンは特殊でアンダーフローとオーバーフローの数がはいっています。これらはDrawした時には表示されません。 横道にそれますが、ROOTではコードをポータブル(OS等の環境に依存しない)ようにするために、Int_t, UInt_t, Float_t等の基本型を定義しています。 基本型のフルリスト: https://root.cern.ch/root/html528/ListOfTypes.html では、ヒストグラムの見た目を変えてみます。 %SYNTAX{lang="c++" numbered="off"}% root [] hist.SetLineColor(kRed) root [] hist.SetLineWidth(2) root [] hist.SetLineStyle(1) root [] hist.SetFillColor(kRed) root [] hist.Draw("") root [] hist.Draw("e") %ENDSYNTAX% 軸のタイトルは後から変えることもできます。 %SYNTAX{lang="c++" numbered="off"}% root [] hist.GetXaxis()->SetTitle("#theta") %ENDSYNTAX% 上の例のように =#= を使うとLaTeXのコマンドの一部を使うことができます。 今回はDrawする時にeオプションをつけることでエラーバーも同時に表示しました。 同じようなことは、グラフィカルな操作も可能です。 ヒストグラムや軸上を右クリックすることで出てくるメニューから色々な操作ができます。 例えば、x軸を右クリックし、SetRangeUserを選ぶと表示する範囲をだいたい変えることができます。 元に戻すにはUnZoomを選びます。 ヒストグラムが選ばれている時はポインタが矢印、軸が選ばれている時は手の形に変わります。 では、ヒストグラムをフィットしてみます。 %SYNTAX{lang="c++" numbered="off"}% root [] hist.Fit("gausn", "", "", -5, 5) root [] hist.GetFunction("gausn")->GetParameter(0) root [] hist.GetFunction("gausn")->GetParameter(1) root [] hist.GetFunction("gausn")->GetParameter(2) %ENDSYNTAX% ここでは、正規化されたガウス分布を-5から5の範囲だけでフィットしています。 =gausn= は正規化されたガウス分布で、パラメータは順に、全区間の積分、mean、 sigmaです。 正規化する必要がない場合は =gaus= にすると、最初のパラメータがmeanでのy軸方向の高さになります。 ATLASでは通常は統計情報はプロットに表示しないようにしていますが、これを表示したい場合には、 %SYNTAX{lang="c++" numbered="off"}% root [] gStyle->SetOptStat(1) root [] gStyle->SetOptFit(1) root [] hist.Draw() %ENDSYNTAX% y軸を対数表示にする場合は %SYNTAX{lang="c++" numbered="off"}% root [] c1.SetLogy(kTRUE) root [] hist.Draw() %ENDSYNTAX% とします。 横道にそれますが、ROOTでは定数はkから始まります。ブーリアンの真と偽はkTRUE, kFALSEと書かれます。(true/falseを使っても大丈夫です。) キャンバスからView->Toolbarをクリックして出るメニューからグラフィカルに操作してプロットの中に線や文字を書くことができます。 同じことは、コマンドラインでTLineやTLatexのオブジェクトを作ることでも可能です。 キャンバスはSaveAs関数で保存することができます。例えば、 %SYNTAX{lang="c++" numbered="off"}% root [] c1.SaveAs("canvas.eps") root [] c1.SaveAs("canvas.pdf") root [] c1.SaveAs("canvas.png") root [] c1.SaveAs("canvas.C") %ENDSYNTAX% ファイルの形式は拡張子によって自動的に変化します。論文やプレゼンテーションで使う場合にはイメージとして保存します。 あとで、プロットの変更(色や線の太さなど)を変えたい時には、マクロとして保存しておくと便利です。 キャンバスは複数のパッドに分割することができます。 ここでは、縦横1x2に分割してみます。 他のオブジェクトをDrawする前にcdすることで描画するパッドを選択することができます。 パッドのインデックスは1始まりです。(0はキャンバス全体を指します。) %SYNTAX{lang="c++" numbered="off"}% root [] c1.Clear() root [] c1.Divide(2,1) root [] c1.cd(1) %ENDSYNTAX% 解析をする時には、イベントに重みを付けたい場合があります。 例えばMCのイベントをフィルする時にはそのイベントが起こる確率を考慮して重み付けます。 その場合には、Fill関数の二つ目の引数に重みを書きます。例えば、 %SYNTAX{lang="c++" numbered="off"}% root [] hist.Reset() root [] hist.Fill(2.5, 1.5) root [] hist.Draw() %ENDSYNTAX% とすると、一つのビンに1.5が入っていると思います。 また、Fillした後に、例えば積分ルミノシティーをかけて全体をスケールしたりします。その場合は、 %SYNTAX{lang="c++" numbered="off"}% root [] hist.Scale(36500) root [] hist.Draw() %ENDSYNTAX% とします。 次に二次元のヒストグラムを作ってみます。 %SYNTAX{lang="c++" numbered="off"}% root [] c1.cd(2) root [] TH2F hist2D("hist2D", "2D histogram", 100, -10, 10, 100, -10, 10) root [] Double_t x,y; root [] for (int i = 0; i < 100000; i++) { x = gRandom->Uniform(-10,10); y = x/10; hist2D.Fill(x, gRandom->Gaus(y, 3)); } root [] hist2D.Draw("lego") root [] hist2D.Draw("colz") %ENDSYNTAX% コンストラクタでは1次元の場合に加えて、y軸のビン数と範囲を指定します。また、Fillにはxとyの値を与えます。 重み付けをする場合は三番目の引数として与えてください。 一番目のDrawはレゴプロットで、z軸を含めて3次元で表示します。マウスのクリックで回転することができます。 二番目のDrawではz軸をカラーコードで表現しています。 2次元分布をフィットするのにはTProfileが使えます。 (他にもFitSlicesYを使う方法もありますが割愛します) %SYNTAX{lang="c++" numbered="off"}% root [] TProfile *prof = hist2D.ProfileX("prof") root [] prof->SetLineColor(kBlue) root [] prof->Draw("same") root [] prof->Fit("pol1") %ENDSYNTAX% ここでは、二次元の分布のプロファイルを作り、一次関数でフィットしました。 ---++ グラフ グラフの基本的なクラスは =TGraph= ですが、ほとんど場合はエラーバー付きのプロットが求められるのでここでは =TGraphErrors= を使います。 %SYNTAX{lang="c++" numbered="off"}% root [] c1.Clear() root [] Double_t xm[10] = {0.29, 0.86, 2.23, 2.99, 4.24, 4.98, 5.72, 6.97, 8.00, 8.87} root [] Double_t xe[10] = {0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3} root [] Double_t ym[10] = {1.30, 0.27, 4.08, 8.82, 14.82, 25.87, 36.40, 50.48, 64.52, 81.53} root [] Double_t ye[10] = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1} root [] TGraphErrors g(10, xm, ym, xe, ye) root [] g.SetMarkerStyle(7) root [] g.Draw("apl") %ENDSYNTAX% この時のオプションの意味は * a: グラフの内容に応じて自動的にプロットのxとyの範囲を設定 (二個目以後のグラフをDrawする時にはaは外す) * p: 点を表示 * l: 点と点の間を直線でつなぐ です。 グラフはそのままフィットできます。 %SYNTAX{lang="c++" numbered="off"}% root [] g.Fit("pol2") %ENDSYNTAX% ---+ ツリー =TTree= はいわゆるn-tuple (複数の値の組み) の列を格納するためのクラスで、データの保存のために広く使われています。例えば、アトラスのAODファイルには、再構成された複数の粒子の情報が各イベントで組みとして、複数のイベントが保存されています。 =TTree= には基本的な型の変数だけではなく、 =vector= やクラス等のオブジェクトも格納することができます。 (ROOTには =TTree= とは別に =TNtuple= というクラスもありますが、私が知る限り使うことはありません。) それでは、実際にツリーを作ってみますが、その前に注意があります。 実際の解析ではツリーの大きさが何百GBになることがあります。 その場合はツリーをRAM上に保持することができなくなり、プログラムが止まったり、他のプロセスに影響を及ぼす場合があります。 それを避けるために、ツリーを作る前に =TFile= を作ってください。これによりTreeは一定の大きさになる度にディスクに書かれるようになります。 (あらかじめサイズが大きくならないことを知っている場合はその必要はありません。) %SYNTAX{lang="c++" numbered="off"}% root [] Int_t EventNumber root [] Double_t METPt root [] Double_t METEta root [] vector<Double_t> muonPt root [] vector<Double_t> muonEta root [] vector<Double_t> MSTrackPt root [] vector<Double_t> MSTrackEta root [] TFile *file = TFile::Open("tree.root", "RECREATE") root [] TTree tree("tree", "tree") root [] tree.Branch("EventNumber", &EventNumber, "EventNumber/I") root [] tree.Branch("METPt", &METPt, "METPt/D") root [] tree.Branch("METEta", &METEta, "METEta/D") root [] tree.Branch("muonPt", &muonPt) root [] tree.Branch("muonEta", &muonEta) root [] tree.Branch("MSTrackPt", &MSTrackPt) root [] tree.Branch("MSTrackEta", &MSTrackEta) root [] for(Int_t i=0;i<100;i++){ root (cont'ed, cancel with .@) [] EventNumber=i; root (cont'ed, cancel with .@) [] METPt=gRandom->Uniform(20,500); root (cont'ed, cancel with .@) [] METEta=gRandom->Uniform(-2,2); root (cont'ed, cancel with .@) [] muonPt.clear(); root (cont'ed, cancel with .@) [] muonEta.clear(); root (cont'ed, cancel with .@) [] MSTrackPt.clear(); root (cont'ed, cancel with .@) [] MSTrackEta.clear(); root (cont'ed, cancel with .@) [] Int_t nmuon = gRandom->Poisson(10); root (cont'ed, cancel with .@) [] for (Int_t imuon = 0; imuon < nmuon; imuon++) { muonPt.push_back(gRandom->Uniform(20,500)); } root (cont'ed, cancel with .@) [] for (Int_t imuon = 0; imuon < nmuon; imuon++) { muonEta.push_back(gRandom->Uniform(-2,2)); } root (cont'ed, cancel with .@) [] Int_t nMSTrack = gRandom->Poisson(10); root (cont'ed, cancel with .@) [] for (Int_t iMSTrack = 0; iMSTrack < nMSTrack; iMSTrack++) { MSTrackPt.push_back(gRandom->Uniform(20,500)); } root (cont'ed, cancel with .@) [] for (Int_t iMSTrack = 0; iMSTrack < nMSTrack; iMSTrack++) { MSTrackEta.push_back(gRandom->Uniform(-2,2)); } root (cont'ed, cancel with .@) [] tree.Fill(); root (cont'ed, cancel with .@) [] } root [] tree.Write() root [] hist.Write() root [] hist2D.Write() root [] g.Write("graph") root [] file->Close() root [] delete file root [] .q %ENDSYNTAX% 上の例では、EventNumber, MET, muon, MSTrackを100イベント分ツリーに書いて保存しています。 METは各イベントに一つしかありませんが、muonは複数あるのでvectorに書き込んでいます。 ツリーがある変数を格納する場所はブランチと呼ばれます。Branch関数でブランチを作成し、それに対応する変数のアドレスをセットします。 基本的な型の場合とvectorの場合ではブランチの作り方が少し違います。 色々な型でブランチを作ることができますが、この例にあるような基本的な型とvectorが使えればほとんど場合は問題ないと思います。 ツリーのFillが呼ばれると、その時の変数の値がブランチにコピーされて、ツリーのエントリーが増えていきます。 最後にツリーをファイルに書いて、ファイルを閉じています。 ついでにヒストグラムとグラフもファイルに書きました。グラフは作った時に名前が設定されていないので、ファイルを書く際に指定しました。 TFileを開く際に"RECREATE"オプションを使うと、同じ名前のファイルがあった時に、上書きされます。 上書きせずに新しいオブジェクトを追加したい場合は"UPDATE"で開きます。 オプションなしの場合は"READ"オプションで開かれるので読み込み専用になります。 ファイルはCloseするまでは完全にディスクに書かれないので気をつけましょう。 次に、今書いたツリーを読み込んでプロットを作りましょう。 もし、どこかで失敗してうまくファイルが作れなかった人は、 /home/sawada/public/tutorial16/tree.root を使ってください。 %SYNTAX{lang="c++" numbered="off"}% $ root tree.root root [] .ls root [] tree->Print() root [] tree->Show(0) root [] tree->Draw("METPt") %ENDSYNTAX% rootコマンドの後にファイル名を指定することで、ファイルを開いた状態でClingが起動します。Clingの中では、ファイルに書かれているものは、その名前をポインタとして使うことができます。 マクロやプログラムでは、このように自動的にポインタが作られるわけではないので、明示的に %SYNTAX{lang="c++" numbered="off"}% TTree *tree = dynamic_cast<TTree*>(_file0->Get("tree")); %ENDSYNTAX% のようにする必要があります。 ツリーのプリントはツリーの情報をブランチのリストを表示します。 ツリーのShow関数は一つのイベントの内容を表示します。 ツリーでプロットを書く際には上の例のようにDraw関数を使うと簡単です。 この時表示されたヒストグラムはhtempという名前で自動的に作られたものですが、Drawを呼ぶ度に作り直されます。 このヒストグラムを後で再利用したい場合や、ビンニングや範囲を指定したいには、以下のようにヒストグラムを作ってからツリーのデータを流しこむようにしてください。 %SYNTAX{lang="c++" numbered="off"}% root [] TH1F frame("frame", "frame", 50, 20, 500) root [] tree->Draw("METPt>>frame") %ENDSYNTAX% 2次元(や3次元)のプロットは変数をコロン:でつなげることでDrawします。 %SYNTAX{lang="c++" numbered="off"}% root [] tree->SetMarkerStyle(8) root [] tree->Draw("METPt:EventNumber") %ENDSYNTAX% 2次元の場合は最初の変数がy軸方向、二番目の変数がx軸方向に来るのに注意してください。 なお3次元の場合は、左から順にx,y,zになります。 全イベントをプロットせずに、イベントを選ぶ場合は、 %SYNTAX{lang="c++" numbered="off"}% root [] tree->Draw("METEta:EventNumber", "") root [] tree->SetMarkerColor(kRed) root [] tree->Draw("METEta:EventNumber", "fabs(METEta)<1&&EventNumber>50", "same") %ENDSYNTAX% のように引数にカットを書きます。カットはC++のの条件式と同じように書きます。 同様にミューオンについてもPtの分布をDrawすることができます。 この場合は条件を満たすミューオン全てがプロットされますので、一つのイベントが複数のエントリーを持っています。 %SYNTAX{lang="c++" numbered="off"}% root [] tree->Draw("muonPt", "fabs(muonEta)<1&&EventNumber>50") %ENDSYNTAX% あえて最初のミューオンだけでプロットを作りたい場合は %SYNTAX{lang="c++" numbered="off"}% root [] tree->Draw("muonPt[0]", "fabs(muonEta[0])<1&&EventNumber>50") %ENDSYNTAX% のようにできます。 次のような間違えはよくあるものなので注意してください。 MSTrack (ミューオンスペクトロメータ)のトラックのPtが大きなミューオンのPtを書きたいと思い次のようなコマンドでプロットを作ったしまった場合、その分布は正しくありません。 (実際はもっと複雑なカットの中で間違える場合が多いですが、ここではその間違えの部分だけ抜き出しています。) %SYNTAX{lang="c++" numbered="off"}% root [] tree->Draw("muonPt", "MSTrackPt>200") %ENDSYNTAX% ROOTで配列やvector型をインデックスなしでDrawした場合は共通のインデックス同士で組み合わせられます。例えば =MSTrackPT[5]>200= であれば、 =muonPt[5]= がDrawされます。 しかし、実際には内部検出器も含めて再構成されるミューオンとMSトラックは一対一に対応していません。(MSトラックにはハドロンのパンチスルーによるトラックが含まれていたり、内部検出器内で散乱してセレクションに通らないミューオンのトラック等も含まれます。)このような状況でROOTは =muonPt= と =MSTrackPt= を対応させる方法を知らないので、Drawで簡単に行うことができません。 少し複雑なカットを伴うプロットを作る際には、Drawは使わずにイベント毎に変数を読みながら解析することをお勧めします。 %SYNTAX{lang="c++" numbered="off"}% root [] Int_t EventNumber root [] Double_t METPt root [] Double_t METEta root [] vector<Double_t> *muonPt = 0 root [] vector<Double_t> *muonEta = 0 root [] vector<Double_t> *MSTrackPt = 0 root [] vector<Double_t> *MSTrackEta = 0 root [] tree->SetBranchAddress("EventNumber", &EventNumber) root [] tree->SetBranchAddress("METPt", &METPt) root [] tree->SetBranchAddress("METEta", &METEta) root [] tree->SetBranchAddress("muonPt", &muonPt) root [] tree->SetBranchAddress("muonEta", &muonEta) root [] tree->SetBranchAddress("MSTrackPt", &MSTrackPt) root [] tree->SetBranchAddress("MSTrackEta", &MSTrackEta) root [] TH1F hmupt("hmupt", "Muon Pt", 100, 0, 500) root [] for(Int_t i=0;i<tree->GetEntries();i++){ root (cont'ed, cancel with .@) [] tree->GetEntry(i); root (cont'ed, cancel with .@) [] Int_t nmuon = muonPt->size(); root (cont'ed, cancel with .@) [] for (Int_t imuon = 0; imuon < nmuon; imuon++) { hmupt.Fill((*muonPt)[imuon]); } root (cont'ed, cancel with .@) [] } %ENDSYNTAX% 上のコードではブランチに対応するローカルな変数のアドレスをセットした後に、ツリーの全イベントでループを回して1イベントずつ処理をしています。 例では単純にミューオンのPtをヒストグラムにFillしていますが、C++でコードを書いてより複雑なカットや計算をすることができます。 例えば、実際の物理解析では、ミューオンとMSトラックのEtaとPhi両方の差を計算して、近いペアを探すことによってマッチングを行ったりします。 参考: たくさんのブランチを読む時に自分でコードを書くのが大変な場合はMakeClassを使うと便利です。 %SYNTAX{lang="c++" numbered="off"}% root [] tree.MakeClass("readtree") %ENDSYNTAX% ---+ マクロ 決まった処理を繰り返し行ったり、コードが複雑になりClingで1行ずつ書くのが大変になった場合はマクロを使うと便利です。 もっとも簡単なマクロはClingでのコマンドを ={}= の中に入れたものです。 以下のコードを =macro.C= というファイルに書いてください。 %SYNTAX{lang="c++" numbered="off"}% { TFile *file = TFile::Open("tree.root"); TTree *tree = dynamic_cast<TTree*>(file->Get("tree")); Int_t EventNumber; Double_t METPt; Double_t METEta; vector<Double_t> *muonPt = 0; vector<Double_t> *muonEta = 0; vector<Double_t> *MSTrackPt = 0; vector<Double_t> *MSTrackEta = 0; tree->SetBranchAddress("EventNumber", &EventNumber); tree->SetBranchAddress("METPt", &METPt); tree->SetBranchAddress("METEta", &METEta); tree->SetBranchAddress("muonPt", &muonPt); tree->SetBranchAddress("muonEta", &muonEta); tree->SetBranchAddress("MSTrackPt", &MSTrackPt); tree->SetBranchAddress("MSTrackEta", &MSTrackEta); TH1F hmupt("hmupt", "Muon Pt", 100, 0, 500); for(Int_t i=0;i<tree->GetEntries();i++){ tree->GetEntry(i); Int_t nmuon = muonPt->size(); for (Int_t imuon = 0; imuon < nmuon; imuon++) { hmupt.Fill((*muonPt)[imuon]); } } hmupt.Draw(); } %ENDSYNTAX% その後、以下のコマンドによってマクロが実行できます。 %SYNTAX{lang="c++" numbered="off"}% $ root macro.C %ENDSYNTAX% 同じことは、次のように、一回rootを立ち上げてから =.x= コマンドで実行することもできます。 %SYNTAX{lang="c++" numbered="off"}% $ root root [] .x macro.C %ENDSYNTAX% この形式のマクロの場合は、今までClingでコマンドを打っていたのと同じように扱われます。 マクロ内で定義した変数(例えばヒストグラム)をそのままCling内で継続して使用することができます。 ただし、改行が文の終わりとは認識されないので、C++のコードを書くように文の終わりにはセミコロンを付けてください。 もう一つの形式はコードを名前のついた関数内に入れます。 例えば、次のようになります。 =macro2.C= というファイルに保存してください。 %SYNTAX{lang="c++" numbered="off"}% Bool_t CheckDEta(Double_t eta1, Double_t eta2, Double_t threshold) { return (fabs(eta1 - eta2) < threshold); } Int_t macro2(Double_t threshold = 0.01) { TFile *file = TFile::Open("tree.root"); TTree *tree = dynamic_cast<TTree*>(file->Get("tree")); Int_t EventNumber; Double_t METPt; Double_t METEta; vector<Double_t> *muonPt = 0; vector<Double_t> *muonEta = 0; vector<Double_t> *MSTrackPt = 0; vector<Double_t> *MSTrackEta = 0; tree->SetBranchAddress("EventNumber", &EventNumber); tree->SetBranchAddress("METPt", &METPt); tree->SetBranchAddress("METEta", &METEta); tree->SetBranchAddress("muonPt", &muonPt); tree->SetBranchAddress("muonEta", &muonEta); tree->SetBranchAddress("MSTrackPt", &MSTrackPt); tree->SetBranchAddress("MSTrackEta", &MSTrackEta); TH1F *hmupt = new TH1F("hmupt", "Muon Pt", 100, 0, 500); for(Int_t i=0;i<tree->GetEntries();i++){ tree->GetEntry(i); Int_t nmuon = muonPt->size(); Int_t ntrack = MSTrackPt->size(); for (Int_t imuon = 0; imuon < nmuon; imuon++) { for (Int_t itrack = 0; itrack < ntrack; itrack++) { if (CheckDEta((*muonEta)[imuon], (*MSTrackEta)[itrack], threshold)) { hmupt->Fill((*muonPt)[imuon]); } } } } hmupt->Draw(); return 0; } %ENDSYNTAX% ここでは一つのファイルに関数が二つ書かれています。この形式のマクロを実行すると、ファイル名と同じ名前の関数(macro2)が実行されます。 実行時に関数に引数を渡すこともできます。 %SYNTAX{lang="c++" numbered="off"}% $ root 'macro2.C(0.05)' %ENDSYNTAX% この場合は最初に紹介した名前なしの場合と違い、関数内にスコープを持つ変数をマクロの実行が終了後にClingから使うことはできません。 (この例では、hmuptはnewで作っているため、関数が終了してもオブジェクトは消されていません。このため、関数から受けた後でもClingが名前に対応するポインタを自動的に用意されてくれるので、hmuptが利用できます。) ここでは、説明のため関数の中でnewしたヒストグラムをそのままにしましたが、呼び出し側でdeleteしないとメモリーリークになるので理由がない限りこのようなコードは書かないでください。 (Cloneを目的とする関数で意図的にこのようにするのはokです。) ---++ ACLiC ROOTのマクロをコンパイルしてから実行する仕組みをACLiCといいます。コードの信頼性が通常のマクロよりも高くなるので、特に重要な計算を行うときには(少し面倒ではありますが)ACLiCを使うことを勧めます。より具体的には、 * 完全にC++のコンパイラーに準拠した文法で書く必要がある * Clingでは宣言されていない変数が使われた時に(パイソンのように)自動的に変数を作ります。便利な反面タイプミスをしてもそのままマクロが実行されてしまいます。ACLiCではそういうことがありません。 * コンパイル時に通常のコンパイラーが(エラーではない部分でも)ワーニングを出してくれます。 * Debugオプションをつけた場合には、クラッシュした時にスタックトレースを見ることができ、デバッグが容易になります。 * コンパイラする際に最適化することができるので、実行速度が速くなります。 先ほどのマクロをACLiCで使えるようには次のように変更します。 違いは、ファイルの先頭でヘッダーをインクルードした点です。(using namespace std;を入れるかどうかは好みの問題だと思います。) %SYNTAX{lang="c++" numbered="off"}% #include <cmath> #include <TFile.h> #include <TTree.h> #include <TH1.h> using namespace std; Bool_t CheckDEta(Double_t eta1, Double_t eta2, Double_t threshold) { return (fabs(eta1 - eta2) < threshold); } Int_t macro2(Double_t threshold = 0.01) { TFile *file = TFile::Open("tree.root"); TTree *tree = dynamic_cast<TTree*>(file->Get("tree")); Int_t EventNumber; Double_t METPt; Double_t METEta; vector<Double_t> *muonPt = 0; vector<Double_t> *muonEta = 0; vector<Double_t> *MSTrackPt = 0; vector<Double_t> *MSTrackEta = 0; tree->SetBranchAddress("EventNumber", &EventNumber); tree->SetBranchAddress("METPt", &METPt); tree->SetBranchAddress("METEta", &METEta); tree->SetBranchAddress("muonPt", &muonPt); tree->SetBranchAddress("muonEta", &muonEta); tree->SetBranchAddress("MSTrackPt", &MSTrackPt); tree->SetBranchAddress("MSTrackEta", &MSTrackEta); TH1F hmupt("hmupt", "Muon Pt", 100, 0, 500); for(Int_t i=0;i<tree->GetEntries();i++){ tree->GetEntry(i); Int_t nmuon = muonPt->size(); Int_t ntrack = MSTrackPt->size(); for (Int_t imuon = 0; imuon < nmuon; imuon++) { for (Int_t itrack = 0; itrack < ntrack; itrack++) { if (CheckDEta((*muonEta)[imuon], (*MSTrackEta)[itrack], threshold)) { hmupt.Fill((*muonPt)[imuon]); } } } } return 0; } %ENDSYNTAX% 実行時は次のように、ファイル名の後に =+= を追加します。 %SYNTAX{lang="c++" numbered="off"}% $ root macro2.C+ %ENDSYNTAX% さらにその後に =g= を追加するとデバッグオプション付きになり、 =O= を追加すると最適化されます。 スタンドアローンのプログラムを書くよりかはずっと簡単に、速くて信頼性の高い処理が行えるACLiCがお勧めです。 ---+ C++プログラムからの利用 ここではあまり深入りしませんがC++のプログラムからROOTのクラスを利用することも可能です。 基本的には、コード内でのROOTのクラスの使い方はACLiCの例で示したようなものになります。 ROOT関連の初期化処理は自動的に行われますので、通常特別なコードを入れる必要はありませんが、グラフィカルな使い方をする場合にはTApplicationというクラスを使う必要があります。 コンパイルとリンクのフラグは =root-config --cflags= と =root-config --glibs= で得ることができます。 独自のクラスをツリーに保存したい時にはディクショナリーを作る必要がありますが、これについてはこの講習の範疇を超えるので割愛します。 ---+ 複雑な関数でのフィット 実際の解析でピークのフィットを行う時には、単純なガウス分布だけでフィットできる場合は稀で、実際は複数の分解能成分があるダブルガウシアン、非対称な分布だったり、バックグラウンドの分布の上に信号のピークがのっていたりします。 ---++ 定義済み関数を使う場合 その場合でもフィットの手続きは変わりませんが、関数の定義が複雑になります。 ダブルガウス分布は次のように定義します。 %SYNTAX{lang="c++" numbered="off"}% root [] TF1 fdgaus("fdgaus", "gausn(0)+gausn(3)", -10, 10) root [] fdgaus.SetParameters(100, 0, 1, 20, 0, 2) root [] fdgaus.Draw() %ENDSYNTAX% TF1のコンストラクタの二つ目の引数で関数系を定義していて、その次の二つで定義範囲を設定しています。 SetParametersでフィットをする前にパラメータの初期値を設定しています。 パラメータの意味は、最初の三つが一つ目のガウス分布の積分値、ミーン、シグマで次の三つが二つ目のガウス分布のパラメータです。 %SYNTAX{lang="c++" numbered="off"}% root [] TH1F hist("hist", "Histogram;x;number of events", 100, -10, 10) root [] for (int i = 0; i < 10000; i++) { hist.Fill(fdgaus.GetRandom()); } root [] hist.Draw() root [] hist.Fit("fdgaus") %ENDSYNTAX% のようにフィットすることができます。 (ここでは、ヒストグラムに関数を使ってイベントを詰めていますが、実際は実験やMCのイベントが入ったヒストグラムをフィットします。) 例えば事前に分布のミーンが0に来ると分かっていてパラメータを固定してフィットしたい場合は、 %SYNTAX{lang="c++" numbered="off"}% root [] fdgaus.FixParameter(1, 0) root [] fdgaus.FixParameter(4, 0) root [] hist.Fit("fdgaus") %ENDSYNTAX% とします。FixParameterの最初の引数はパラメータのインデックス。二つ目は固定する値です。 また、フィットでパラメータが動ける範囲を制限したい時には、 %SYNTAX{lang="c++" numbered="off"}% root [] fdgaus.SetParLimits(2, 0, 3) root [] fdgaus.SetParLimits(5, 0, 7) root [] hist.Fit("fdgaus") %ENDSYNTAX% のようにします。引数の意味は、パラメータのインデックスと、制限する範囲です。 ---++ 式の形を指定する場合 解析的な式を書いて関数を定義することもできます。 %SYNTAX{lang="c++" numbered="off"}% root [] TF1 fdcos("fdcos", "[0]*exp(-x/[1])*cos(x)", 0, 100) root [] fdcos.SetParameters(20, 10) root [] fdcos.Draw() %ENDSYNTAX% =[0]= と =[1]= が関数のパラメータを表していています。今は =[1]= までなので、自動的に二つのパラメータを持つ関数が作られます。 ---++ C++の関数を使うばあい さらに複雑になり、1行で書くのが難しい場合は次の例のように関数を定義します。 = func.C =というファイルを開いて以下の内容を書いて保存してください。 %SYNTAX{lang="c++" numbered="off"}% Double_t ExpGaus(Double_t *x, Double_t *par) { //http://pibeta.phys.virginia.edu/~pibeta/docs/publications/penny_diss/node35.html //par[0] : height of Gaussian //par[1] : peak position //par[2] : sigma of Gaussian //par[3] : transition point between gaussian and exponential Double_t fitval; if(x[0] > par[1] + par[3]) { if(par[2] != 0) { fitval = par[0] * TMath::Exp(-1 * (x[0] - par[1]) * (x[0] - par[1]) / 2 / par[2] / par[2]); } else { fitval = 0; } } else { if(par[2] != 0) { fitval = par[0] * TMath::Exp(par[3] / par[2] / par[2] * (par[3] / 2 - (x[0] - par[1]))); } else { fitval = 0; } } return fitval; } %ENDSYNTAX% これは、あるxを境にガウス分布と指数分布が切り替わる関数で、エネルギーのリークが無視できない場合のカロリーメータのレスポンス関数で使われます。 パラメータは4つあり、関数の先頭で説明されている通りです。 それでは、このC++の関数を使って、ROOTの関数を作ってみます。 %SYNTAX{lang="c++" numbered="off"}% root [] .L func.C root [] TF1 fegaus("fegaus", ExpGaus, 0, 100, 4) root [] fegaus.SetParameters(100, 50, 2, -1) root [] fegaus.Draw() %ENDSYNTAX% 1行目で =ExpGaus= 関数を読み込んでいます。この形式でTF1を作る場合にはROOTにパラメータの数を教えてあげる必要があるので、TF1のコンストラクタに4を与えています。 このようにTF1を使う場合のC++の関数の引数と返り値の型は決まっているので必ず例にある通りにしてください。 ---++ さらに難しいFit 今まで紹介した方法で、1次元の分布に対するフィットはほとんど場合こなせると思いますが、他変数のフィットや、変数に制限を加える場合には、自分でカイ二乗の計算を定義してMINUITパッケージを利用してフィットしたり、尤度関数を定義してRooFitパッケージを使ってライクリフッドフィットを行ったりします。 これらについては、この講習の範疇を超えますので、割愛します。 ---+ 発展 この講習ではROOTのほんの一部だけ紹介しました。 実際には、Graphical User interface、ソフトウェア開発のためのシステムツール(マルチスレッド、ソケット通信...)、統計解析のためのパッケージ(RooFit)、ピークサーチ、フーリエ変換、多変量解析等、実験物理研究者がソフトウェアを書く際に便利な色々なツールが入っています。 ---+ 練習 練習として、実際にATLASの解析で使われているntupleを使ってプロットを作ってみましょう。 /home/sawada/public/tutorial16/sample.rootはZがeeに崩壊するMCから作られたntupleです。 元々の目的は、Tag-and-probe法によって、電子がPixel検出器とSVTの間で散乱されてトラックとして再構成されない確率をPtの関数として計算することです。 Zeeのうちの一方の再構成された(陽)電子をTagとして、Probe側の(陽)電子がどのように見えるかを調べます。 Probe側の(陽)電子は必ずしもトラックとして再構成されないので、カロリメータの情報からPtや方向を再構成します。 練習として、このファイルを使って、 * Tag側とProbe側のファイ方向の散布図 * 電子と陽電子で組んだ不変質量分布 を作ってください。 以下のブランチを使ってください。 タグ側Pt : ZeeProbClusterTagPt タグ側Eta : ZeeProbClusterTagEta タグ側Phi : ZeeProbClusterTagPhi プローブ側Pt : ZeeProbClusterProbeClusterPt プローブ側Eta : ZeeProbClusterProbeClusterEta プローブ側Phi : ZeeProbClusterProbeClusterPhi 不変質量の答え: ZeeProbClusterM これらのブランチの型はvector<float>です。 すでにZeeらしいTagとProbeのペアになっているので、同じインデックスでペアとして使って大丈夫です。 不変質量の計算は自分で書いてもよいし、TVector3, TLorentzVector等の便利なクラスを使ってもよいです。 これらを使う場合は、ROOTのClass referenceを参照してください。 https://root.cern.ch/doc/master/classTVector3.html https://root.cern.ch/doc/master/classTLorentzVector.html %STOPINCLUDE%
Attachments
Attachments
Topic attachments
I
Attachment
History
Action
Size
Date
Who
Comment
png
graph.png
r1
manage
11.1 K
2016-12-26 - 02:43
RyuSawada
png
histogram.png
r1
manage
22.3 K
2016-12-26 - 02:43
RyuSawada
This topic: Main
>
TWikiUsers
>
RyuSawada
>
AtlasJapanSoftwareTutorial16DecROOT
Topic revision: r9 - 2016-12-27 - RyuSawada
Copyright &© 2008-2023 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