- ROOT及びROOTに関する一般的事柄
・ROOTの起動は、コマンドライン上でrootと打つと使うことができる。(これを有効にするためにはちゃんとパスを通しておく必要がある。)
・ROOTの起動においては、オプションコマンドがあり、-lをつけることでスプラッシュを非表示で起動、-bをつけてバッチモードで起動できこれはsshなどで使う場合に便利、-qをつけておくとマクロの読み込みの終了と同時にrootが終了するようにできる。
・ROOTは、インタープリター(対話式、1行ずつ解釈して実行すること)であり、ROOTのコマンド(ROOTで予め用意されているクラス。例えばTH1Fなど)を打ち込んで使う。
・また、root []←この状態で{を打つと、複数行にかけてマクロを書くことができる。
・ROOTのコマンドラインは、CINTと呼ばれるCまたはC++用(で書かれた)のインタープリターになっているので、とても応用性が高い(しかもフリーのソフトウェア)。ここにCINTの機能などが分かりやすくまとめられている。
・CINTは、実行時間より迅速な開発が重要な状況で大変便利であり、このインタープリターを用いることで、コンパイルとリンクのサイクルがかなり短縮され、迅速な開発ができるようになる。
※今までCINTはCERNによってサポートされていたが、2013年にサポートを終了し、CINTのサポートは元の開発者である後藤さんに戻り、CERNでは新しいインタープリターであるClingに切り替わった。Clingに関してROOT公式ページに説明が書いてあったので、紹介程度に以下に記載。
・Clingとは、LLVMとClangライブラリ上に構築された対話型のC++インタープリターである(言い換えると、バックエンドでLLVM/Clangを用いたClingというC++インタープリターが動いている)。これは、CERNのROOTフレームワークを支援するグループにより、C/C++インタープリターに取って代わるものとして開発された。
他の標準的なインタープリターに対する利点は、コマンドラインプロンプトを持っていて、コンパイルにJust-In-Time(JIT)コンパイラを用いることである。JITコンパイラとは、ソフトウェアの実行時にコードのコンパイルを行い、実行速度の向上を図るコンパイラのこと。通常のコンパイラは、ソースコードから機械語への変換を実行前に事前に行う。このことから、このコンパイラをJITコンパイラに対してAhead-Of-Time(AOT)コンパイラと呼ぶ。
Clingの目標の一つが、ROOTにおける現在のC++コンパイラに替わるハイパフォーマンスを供給することらしい。CINTとの下位互換性は現在開発中とのこと。
Clingは、Clangができる全てのことができ、いくつかのインタープリター特有のC++拡張機能が解析できる。
Clingでは、以下のようなメタプロセッサコマンドが使える。
・.x filename.cpp - filename.cppというマクロをロードして、もし、filename()がきちんと定義されていれば、filename()が走る。
・.L library | filename.cpp - libraryまたはfilename.cppをロードする。
・.printAST - 各処理済み実体の後に抽象構文treeを表示する。
・.I path - includeパスを追加する。
補足的にLLVMとClangについて説明する。
・LLVMは、Low Level Virtual Machineの略で、C++で書かれたコンパイラ基盤で、任意のプログラミング言語で書かれたプログラムのコンパイルや実行時の最適化のために設計されたもの。
・Clangは、CまたC++などのプログラミング言語用のコンパイラフロントエンドである。バックエンドとしてLLVMが使われている。Clangのゴールは、GNU Compiler Collection (GCC)に代わるものを提供することで、Appleが開発を後援している。これは、フリーソフトウェアライセンスの下で利用可能。
調べてたら、めっちゃ色々書いてしまった、、笑
- マクロの実行について
・決まった処理の繰り返しや、ソースコードが複雑になってClingで1行ずつ書くのが大変になった場合は、マクロを使うと便利。
・最も簡単なマクロはClingでのコマンドを{}の中に入れたもの。この形式のマクロは、Clingでコマンドを打っていたのと同じように扱われる。注意点として、改行が文の終わりとは認識されないので、C++のコードを書くように文の終わりにはセミコロンをつける！
・もう一つの形式は、ソースコードを名前のついた関数の中に入れること。例えば、ファイルの中に関数が2つあれば、実行する時はファイル名と同名の関数が実行される。
以下にマクロの実行方法を記載する。方法としては3種類あり、例としてmacro.cppがあったとする。
1つ目は、ターミナルからrootを開くときに、$ root macro.cpp とする。
2つ目は、rootを開いて、root [0] .x macro.cpp とする。ここで、.xはマクロを実行するという意味。
3つ目は、rootを開いて、root [0] .L macro.cpp+ でマクロを読み込み、マクロ内で定義した関数を実行する。ここで、.Lはマクロの読み込みという意味。(.Lについての補足情報は(*)を参照)
ex)3つ目に対して、マクロ内で例えばmacro()を定義していれば、root [1] macro()のようにすると、macro()の中身が実行される。
注意点として、一つ目と二つ目は、ファイル名と実行する関数名は同じにしておかないといけない。
三つ目は、ファイル名と関数名は一致してなくてもよく、マクロ内で定義した関数であればどれでも実行できる。
・ROOTでマクロを実行する上で、大規模な解析をインタープリター形式で行うには限界がある。そこで、マクロを実行する別の方法として出てくるのが"ACLiC"を用いる方法である。
・ACLiCは、Automatic Compiler of Libraries for Cling??(Clingでなく、CINTの時はAutomatic Compiler of Libraries for CINTであった)の略で、マクロをコンパイルし共有ライブラリを作り、Cling上で実行することができる。
(*)ACLiCを呼び出すための方法として、あるマクロに対して、.L macro.cpp+のように.Lと+をつけることによって、ACLiCが呼び出され、オンザフライでコンパイルまた同時に共有ライブラリに変換する。
→この時、必要なヘッダファイルがincludeされてないと、コンパイルは失敗するので注意
・ROOT公式ページによると、インタープリターの代わりにコンパイル済みマクロとしてマクロを走らせることができ、ACLiCによりデバッグをしてくれたり、わずかだが高速で走らせることができるらしい。
「Clingと比較して」
- コードの信頼性が通常のマクロより高くなるので、重要な計算などを行う時はACLiCを使う方がいい。
- 文法がC++と同じで正しい(つまり、完全にC++のコンパイラに準拠した文法で書く必要がある)。
- コンパイル時に最適化することができるので、実行速度が高速
- Clingでは、宣言されていない変数が使われたとき自動的に変数を作るが、ACLiCではそういうことがない。
- コンパイル時に通常のコンパイラがWarningを出してくれる。
- Debugオプションをつけたとき、クラッシュした時などにスタックトレースを見ることができ、デバッグがしやすくなる。
「普通のコンパイルと比べて」
- コンパイルが簡単(マクロの読み込みに"+"を付け足すだけ。これをするだけでROOTが自動的にコンパイラを起動し、マクロを動的ライブラリとしてロードする。)
- 作成された共有ライブラリはClingで実行されるため、マクロ実行後シームレス(連続的に)にインタプリターの操作に移れる。
- ACLiCがコンパイル時に行うエラー解析は、ROOT用なので普通のコンパイラーよりも分かりやすい。
使い方は、例にmacro.cppを用意し、以下のように実行する。
root macro.cpp+ -->マクロをコンパイルして実行(ACLiCの利用。マクロが更新されてなかったら再コンパイルしない。)
root macro.cpp++ -->マクロを再コンパイルして実行
root macro.cpp+g -->マクロをデバッグ情報を付加してコンパイルして実行
root macro.cpp+o -->マクロを最適化してコンパイルして実行
もしくは、ROOTを起動し、root [] .x macro.cpp+ のようにしても良い。
ACLiCを用いる点での注意点として、
- Cling独自の機能は使えない
- 予め必要なヘッダーファイルをincludeする必要がある。
- rootfileについて
・通常、rootでroot fileを開きたかったら、root fileの名前がraw0001.rootとだったら、ターミナル上で$ root raw0001.rootとすれば、raw0001.rootが読み込める。
↑上記のやり方でもrootfileは読み込めるが、これは、rootを開いてTFile *_file0 = TFile::Open("raw0001.root")をしているのと同じ意味である。
例えば、あるTTreeがあり、そこからそのTTree内にあるBranchのヒストグラムを対話式で表示させようとする方法として、以下のように書けばできる。
[TTreeの名前]->Draw("[表示させたいBranchの名前]>>[適当にヒストグラムの名前定義](bin数, 最小値, 最大値)", "[条件式]", "[オプション]")
・ここで、上にも書いてある通り""の第一引数はヒストグラムの定義、第二引数は条件式で、例えばgammaEnergy>5000&&gammaEnergy<10000でgateをかけたり、ある特定のchannelだけ除いたhistogramを畫きたかったら、"channel!=1"のようにすればよく、第三引数はオプションで、例えばヒストグラムをエラーバ付きで表示させたかったら、"e"とか一つのキャンバスに複数のヒストグラムをDrawしたかったら、"same"とかすればいい。
以下に例を示す。
sampleTree->Draw("gammaEnergy>>h1(40000, 0, 40000)", "neutronEnergy>10&&neutronEnergy<20 & channel==1", "e")
↑上記の場合は、sampleTreeの中のBranchの一つである、gammaEnergyをDrawする操作。第一引数では、ヒストグラムの名前をh1とし、bin数を40000、最小値・最大値をそれぞれ0、40000とする。第二引数では、10~20の間でneutronEnergyでgateをかけると同時にchannel1のみを抽出する。第三引数では、eをつけることによりヒストグラムにエラーバー付きで表示させる。
・すでに一度root fileを読み込んだ状態で、もう一つ別のroot fileを新たに読み込むことができる。読み込みたいroot file名が"raw0002.root"だったら、
TFile *_file1 = TFile::Open("raw0002.root")
とすると、raw0001.rootに引き続き、raw0002.rootが読み込まれる。読み込まれただけでは良くなくて、raw0002.rootの中身をいじりたい場合は、_file1->cd()を行い、trueと返ってくることにより、raw0001.rootからraw0002.rootに移動することができる。
・中身が同じ形式のroot fileが複数あるとき、それを一つのroot fileにまとめる方法として、ターミナル上でhaddコマンドを実行する。
やり方の例を挙げる。まず、run1.root、run2.rootというroot fileがあり、これらをall.rootとして一つにまとめることを考える。以下のように実行すればいい。
$ hadd all.root run1.root run2.root
こうすることで、all.rootとして一つにまとめることができる。
- TChainについて
これは、同じ構造のTTreeが入った、複数のROOTファイルを読み込んで解析する場合に使用する。
読み込みたいTTreeの名前を今回はtreeとし、それぞれのrootfile名をtest1.root、test2.root、test3.rootとする。実際の書き方は以下の通り。
TChain *chain(適当に名前を定義) = new TChain("tree","tree title"); ここで、第一引数は、読み込むTTreeの名前、第二引数は、なくても良い。
chain->Add("~/rootfile/test1.root");
chain->Add("~/rootfile/test2.root");
chain->Add("~/rootfile/test3.root");
Addするときは、loopを回して読み込むこともできる。以下に使い方の例を示す。
TChain *chain = new TChain("chain", "chainname");
const Int_t fNFile=11;
Int_t iFile;
for (iFile=0; iFile
chain->Add(Form("~/rootfile/test%d.root", iFile+10));
}
- TBrowser
全てのROOTオブジェクトを見ることができる。使い方は以下の通り。
root [0] TBrowser [何か文字]
[何か文字]の部分は、文字だったらなんでもいい。例えば"a"とか"blue"とか。
開くとウィンドウの左側にリストができ、参照できる全てのROOTオブジェクトが表示される。どれかを選択すると、右側のボックス内に表示される。簡単にrootファイルの中身を見れるので、便利。
TBrowserにおいて縦軸のrangeを変更するときは、縦軸にカーソルを合わせ手のマークになったところで、右クリック→SetRangeUser→上から順に軸の下限と上限を入力したら変更可能。ex)下限を0、上限を10^5にしたかったら、上から順に0、1.0e+5と入力すればいい。
- TGraph
TGraphは、テキストファイルからグラフが書けたりする便利なやつ。
例えば、data.datというテキストファイルの中身が以下のようだとする。
100 1020
200 258
300 111
400 87
これをROOT上で描きたかったら、ROOTを開き、TGraph *g = new TGraph("data.dat")としてまず空の箱を作る。
そして、g->Draw()したら描ける。Drawのオプションとして、Draw("A*")としたらポインターがアスタリスクみたいな形になる。
- ヒストグラムのfitについて
TF1を使って、ROOTで作成したヒストグラムをガウシアンや任意の関数でfitすることができる。
ROOTにはいくつかデフォルトで入っている関数があり、それらでfitする場合は以下のように書いたらfitできる。
histo->Fit("gaus")
histo->Fit("landau")
histo->Fit("expo") //(expo=exp(p0+p1*x))
histo->Fit("polN") //(polN=p0+p1*x+p2*x^2+...+pN*x^N)
また、fitする範囲も指定することができる。
histo->Fit("gaus","","",-1,1)のように書くと、-1〜1でfitされる。
⭐︎ここからは、TF1コンストラクタを用いて任意のfit関数の定義の仕方とfitのやり方を紹介する。
やり方の順番としては、まずTF1でfit関数を定義して、それからfitting parameterを決めてfitする。
1.fit関数定義:TF1 *name = new TF1("fit関数名","関数の中身",最小値,最大値)
2.パラメータ設定:fit関数名->SetParameters(.....)※fitting parameterの数に合わせて変える。
3.fit実行:histo->Fit(fit関数名)※fit範囲を制限するときは、histo->Fit("fit関数名","","",min,max)のように書く。
※TF1はfitするときはもちろんのこと、関数を定義してDrawするときにも使える! 以下に例を示す。
root [] TF1 *f1 = new TF1("f1","cos(x)", -10, 10);
root [] f1->SetNpx(1000); //これはつけるかどうかは自由で、点数を増やしてより精度の高いグラフを描きたかったらつけたらいいと思う。デフォルトのSetNpxは100っぽい。
root [] f1->Draw();
統計ボックスについて
・ヒストグラムを書いた時などに、デフォルトでは右上にEntriesやStd Devなどの統計情報が記載されているボックスができる。これは、次のコマンドを打つことで、表示させる項目を変更することができる。
gStyle->SetOptStat()
()の中の数字を変えることで、変更可能。デフォルトは(000001111)。
(abcdefghi)で、
a = 1; kurtosis printed
a = 2; kurtosis and kurtosis error printed
b = 1; skewness printed
b = 2; skewness and skewness error printed
c = 1; integral of bins printed
c = 2; integral of bins with option "width" printed
d = 1; number of overflows printed
e = 1; number of underflows printed
f = 1; rms printed
f = 2; rms and rms error printed
g = 1; mean value printed
g = 2; mean and mean error printed
h = 1; number of entries printed
i = 1; name of histogram printed
・フィッティングを行った後、フィッティングの結果の表示と統計の表示の仕方は、gStyle->SetOptFit()とすれば、フィッティングの結果の表示が図にされる。非表示に戻すには、gStyle->SetOptFit(0)で全ての情報を表示するには、gStyle->SetOptFit(1111)。
gStyle->SetOptFit()
()の中の数字を変えることで、変更可能。デフォルトは、(0111)。
(abcd)で、
a = 1; print Probability
b = 1; print Chisquare(Chi2)/Number of degrees of freedom(Ndof)
c = 1; print errors (もしc=1なら, dは1になる)
d = 1; print name/values fo parameters
・論文などに載せる時、統計ボックスに記載する情報を、χ2/ndfとフィッティングパラメータのみにしたい時は、まずgStyle->SetOptStat(0)をして(なぜか2回やらないと消えなかった)、次にgStyle->SetOptFit(1)をすればできる。
- RooFitについて
・RooFitはFitツールの一種である。ヒストグラムやグラフをある関数でfitするときは通常のROOTのfitで十分であるが、Roofitには以下のような特徴がある。
・最尤法なので、binningしてないdataにfitできる。
・多次元fitや多次元関数の積分ができる。
・2つのpdfを合成して新たなpdfを作成できる。(足して1次元pdfを作る-->f(x)+g(x)、掛けて2次元pdfを作る-->f(x)*g(y)、parameterを他のpdfで変化させて2次元pdfを作る-->f(x,g(y))、convolution(畳み込み)ができる(高速フーリエ変換を使用できる))
・Convolutionをするとき(f(x)○g(x))は、RooFitの中のRooFFTConvPdfというクラスを使う。
・コマンドラインでRooFitクラスの一つを参照するとすぐに、ROOTは自動的にRooFitライブラリ(libRooFitCoreとlibRooFit)をロードする。マクロでは便宜上以下の1行を追加しておく。
using namespace RooFit;
・RooFitの重要な概念は、モデルがオブジェクト指向(複雑な処理などを無視して効率よく利用できるようにすること)で構築されることである。各RooFitクラスは、数学オブジェクトと1対1対応している。変数を表すクラスRooRealVar、関数を表す基本クラスRooAbsReal、確率密度関数を表す基本クラスRooAbsPdfなどがある。最も単純な数学関数でさえも複数のオブジェクト(関数自体とその変数)からなるので、全てのRooFitモデルも複数のオブジェクトからなる。
・ここで、ガウス確率密度関数の構築の例を以下に示す。
[ガウス確率密度関数の構築]
RooRealVar x("x","x", -10, 10);
RooRealVar mean("mean", "Mean of Gaussian", 0, -10, 10);
RooRealVar sigma("sigma", "Width of Gaussian" 3, -10, 10);
RooGaussian gauss("gauss", "gauss(x, mean, sigma)", x, mean, sigma);
1行目でgaussで使用される各変数はいくつかのプロパティで初期化される。順番にname、title、range。
最終行で、RooGaussianで実装されているように、ガウス確率密度関数(PDF)を作成している。RooGaussianは、抽象基本クラスRooAbsPdfの実装である。PDFガウスは、変数オブジェクトと同様にname、titleを持ち、コンストラクタで渡される参照を通じて変数x、mean、sigmaにリンクされている。
続いて、以下の記述を実行することにより、ガウシアンPDFを作成できる。
RooPlot* xframe = x.frame(Title("Gaussian PDF")); //ある変数(ここではx)を横軸にした時の縦軸の値をplotできる。Title()は"RooGlobalFunc.h"で定義されている関数の一つで、RooFitオブジェクトにtitleを与えることができる。
gauss.plotOn(frame); //pdfをxframeにplotする
frame->Draw(); //TCanvasに描く
上記を行う際は、"RooRealVar.h"、"RooDataSet.h"、"RooGaussian.h"、"TCanvas.h"、"TAxis.h"のincludeまた、using namespace RooFitを定義しておく必要がある。
また以下のような記述を加えると、複数のオブジェクトを一つのキャンバス内に書ける。例えば、パラメータsigmaの値を変えて、ガウシアンを2回書くことができる。
RooPlot* xframe = x.frame();
gauss.plotOn(frame);
sigma = 2;
gauss.plotOn(frame,LineColor(kRed));
frame->Draw();
この例では、代入演算子を用いて、plotOn()コマンドの後ろにRooRealVarシグマの値を変更する。LineColor(kRed)で曲線の色を赤に変更(デフォルトの色は青らしい)。
・ここからは、RooFitを用いて2種類の関数を畳み込んで、畳み込んだ関数でdataをfitするやり方を紹介する。
一般に、実験データ分布は、検出器応答関数によって補正される理論的分布の結果である。
最も一般的なケースは、ある理論モデルT(x,a)と検出器応答関数R(x,b)の畳み込みによって記述される。
[ 畳み込みの式をおいおい書く予定]
・畳み込み積分の明示的な計算を以下に示す。
- ROOTにおけるクラスの使用方法
・TMath::Pi()
これは、ROOT内で円周率πを定義するクラス。
・TMath::Factorial←()がないことに注意
これは、階乗を定義するクラス。使用方法は、TMath::Factorialと定義してから、求めたい階乗の値をTMath::Factorial(4)のように入力。この場合は4!なので、24と返ってくる。
・ROOTにおいて、例えばベッセル関数など特殊関数を自分で定義することなく、使うことができる。それが、MathCoreというものである。これらを使うには、まず"Math/SpecFunc.h"をincludeする必要がある。"Math/SpecFuncMathCore.h"も？
次に、例にwigner6j記号を取って具体的な使い方を書く。なお、使える特殊関数はここに載っている。
1. ターミナル上でROOTを開き、初めにgSystem->Load("libMathMore")を打つ。MathCoreを使うにはこれを最初にロードする必要がある。
2. ROOTにデフォルトで設定されているwigner6jの計算方法について説明する。ROOTでは、ROOT::Math::wigner_6j(two_j1, two_j2, two_j3, two_j4, two_j5, two_j6)の表式で定義されている。つまり、引数には元々の数値を2倍したものを入れる必要がある。また、引数にはDouble型は使えず、Int型のみしか使えないので注意が必要。
・TH1::SetDefaultSumw2()
これは、全てエラーバー付きで統計誤差の誤差伝播をしてくれる
- ROOT(C++)におけるマクロでのコンパイルエラーの解決法
・treating Unicode character as whitespace [-Wunicode-whitespace]
これはマクロの中でどこかに全角のスペースがある時に出る。
→全角スペースを半角に変えてやると解決する。
・unknown type name
これはROOTのコマンドを使おうとして必要なインクルードファイルがインクルードされてない時に出たりする。
→対応する型が定義されたヘッダファイルを、きちんとインクルードしたら大丈夫。
・use of undeclared identifier
これは宣言してないものを使ってると言ってる。
→マクロ内でちゃんと宣言したら直る。
・excess elements in array initializer
これは宣言している配列の数以上を、配列に入れた時に出るエラー。
→ちゃんと宣言してる配列の数に直せば直る。
・control may reach end of non-void function
これは、マクロをint main()のように書いてるのに、void()みたいに戻り値なしでreturnしてない時に出るエラー。
→returnを最後に付け足せば直る。※ returnを書いてても}の位置がおかしかったらこのエラーは出る。
・redefinition of "add"
これは、すでに一度定義してあるものをもう一度定義してる時に出るエラー。
→別のもので再定義してやれば直る。
・expected ';' at end of declaration
これは、終了を表す";"を付け忘れている時に出るエラー。
→付けてやれば直る。
・cannot initialize a variable of type 'int' with an lvalue of type 'const char [5]'
これは、クラスが衝突して発生するエラー。
→大量のエラーが列挙されるが、エラーの原因は一つだけだったりする。
・implicit instantiation of undefined template
これは、includeせずにマクロを書いた時に起こるエラー。
→必要なincludeファイルをincludeすれば解決する。
・extraneous ')' before ';'
これは、";"の前に余計な()が付いている時に吐くエラー。
→余計な()を消すと、解決する。
・expected unqualified-id
予約語の関数を再定義したり、再宣言してるときに出るエラー。
→予約語以外に定義や宣言を変えると解決する。
・*** Break *** segmentation violation
これがどういう時に出るか正確には理解できてないが、自分が経験した限りで大部分は、fileのpath指定の間違いや読み込むfile内のグラフなどの名前の間違いであることが多かった。
→きちんと名前などを確認して、正しくなかった場合は正確な名前に直す(合ってると思っても、案外細かいところで違ってることがある。例えば、"0129"としないといけないところを"129"として定義していたとか)。
・マクロを動かし終わった後に出る、(TFile::Append): Replacing existing TH1: hist (Potential memory leak)というエラー。
これはTH1Dなどでヒストグラムを定義した時に、TH1D *h = new TH1D()とTH1D *h2 = new TH1D()でhとh2で名前は分けてるが、()の中身が全く同じ時に出るエラー。
→()の中身を違うものに書き換えると解決する。
・ISO C++11 does not allow conversion from string literal to 'char *' [-Wwritable-strings]
これは、文字列リテラルから'char'型への変換は非推奨であるとの警告。C++では文字列リテラル(0文字以上の連続した文字列を示す定数)は、変更できないconst char配列として扱われるらしい。
→これを解決するためには、*char型をconst char*型に変更すればOK。
- 1次元ヒストグラムにおいて
・通常、1次元のヒストグラムを作るときは、以下のように定義して書く。
TH1D *h1 = new TH1D("name", "title;XXX;YYY", 100, -25., 25.);
"XXX"はx軸title、"YYY"はy軸titleになる。100はbinの数で、-25.と25.は横軸を-25から25まで取るということ。
以下に例を示す。
TH1D *h1 = new TH1D("tof_hist", "tof hist", 20000, 0, 40000);
上のヒストグラムの定義で、tof_histが右上のボックスの名前でtof histがヒストグラム自体の名前となる。
1binの幅は、第3引数÷第1引数で求めることができ、今回の場合は40000/20000=2より、1binの幅は2[横軸の単位]となる。
ここで、同じマクロ内でもう一つヒストグラムを定義する時に、上で定義したやつと同じnameで新たにnewしようとすると、"Warning in : Replacing existing TH1: name (Potential memory leak)."と出てしまう。-->これを消すためには、別の名前で定義してやらないといけない。
また、loop内ではnameをhprof[j] = new TH1F(Form("h%d", j), "", Bins[j], Xmin[j], Xmax[j]);などで変えてあげるといい。
・
- ヒストグラム及びグラフにおいての使い方色々
【ヒストグラム(TH1,TH2,TH3)編】
・Tree->Drawにおいて、ヒストグラムの上限と下限を定義する時は、演算過程ではなく演算結果を定義しないと、正しいヒストグラムが書けない。
ex) (20000,0.39219189658435688051,20000*1.0376873428237611829+0.39219189658435688051)→(20000,0.39219189658435688051,20754.13905)のようにする。
ここから、定義したヒストグラム名を"h"として話す。
・h->GetNbinsX()をすると、書いたヒストグラムのbin数が出力
・h->GetXaxis()->GetBinUpEdge(histo->GetNbinsX())で、ヒストグラムの最大値(最大のbinの右端のxの値)が知れる
※これが自分が定義したヒストグラムの最大の値と一致してなかったら、何かが間違っている。
・h->GetXaxis()->GetBinUpEdge(bin番号)で、指定したbinの右端のxの値を返す。
・h->GetEntries()で、総イベント数を出力。
・h->FindBin(横軸の値)で、横軸の値に対応したbin番号が知れる。
・h->GetBinContent(bin番号)で、bin番号のイベント数が知れる。
・h->GetBinError(bin番号)で、そのbin番号のイベント数のエラー(誤差)を知れる。
※補足として、をDrawする時にエラーバー付きでDrawしたければ、オプション"e"をつけて、Draw("e")のようにすればいい。これで出力されるエラーバーは、各binの統計誤差を誤差伝播している。(自分で確かめたので多分間違いない)。
これが気持ち悪く、きちんと確かめたければ、tmp = h->GetBinContent(bin番号), error = h->GetBinError(bin番号)として上で、h->SetBinContent(i, tmp), h->SetBinError(i, error)とし、h->Draw()とすれば、誤差付きでDrawされる。
これも面倒くさければ、main関数の中にTH1::SetDefaultSumw2();を入れておくと、自動的に誤差伝播をしてくれて、Drawした時にエラーバー(統計誤差、イベント数Nの√N)付きでDrawしてくれる
・h->GetMaximum()で、全てのbinの中での最大値を出力。GetMinimum()は逆に最小値を出力。
・h->GetMaximumBin()で、最大のイベント数が入ったbin番号が知れる。同様にGetMinimumBin()は逆のことをする。
つまり、h->GetBinContent(h->GetMaximumBin())とh->GetMaximum()は同じ意味。
・作った1次元ヒストグラムにおいて、binの数、binの最小値、binの最大値を得るためにはそれぞれ以下のようにしたらいい。
h->GetXaxis()->GetNbins()(h->GetNbinsX()でもいけそう)、h->GetXaxis()->GetXmin()、h->GetXaxis()->GetXmax()
・先ほどの続きでbinの下端、中心値、上端を得るにはそれぞれ以下のようにしたらいい。
h->GetXaxis()->GetBinLowEdge(i)(ここでiはbin番号を入力)、h->GetXaxis()->GetBinCenter(i)、h->GetXaxis()->GetBinUpEdge(i)
・h->SetBinContent(bin番号,代入したい値)とすると、入れたい値をそのbin番号に入れることができる。
・ヒストグラムをキャンバスに保存する時は、[キャンバスの名前]->SaveAs("保存したい名前と拡張子")と書く。
ex) c1->SaveAs("title.pdf")のような感じ。キャンバスを特に設定しなかったら、デフォルトのキャンバス名は"c1"となる。拡張子は.pdfじゃなくて、.jpgや.pngや.rootなども使える。
・h->SetXTitle("X軸名")やh->SetYTitle("Y軸名")でX,Y軸にタイトルが設定可能。これはグラフでは使えない。
h->GetXaxis()->SetTitle("X軸名")やh->GetYaxis()->SetYitle("Y軸名")であれば、ヒストグラムやグラフでも使える。
同様の使い方として、h->SetTitle("Title")とすると、グラフやヒストグラムの名前を定義できる。また、複数のグラフやヒストグラムを一つのキャンバスに描き、それぞれにタイトルをつけたかったら、h->SetTitle("h1")とh->SetTitle("h2")のようにするとできる。余談だが、キャンバスの右上の空白のところで右クリックをしてBuildLegendを選択すると、定義したタイトルでレジェンドを作ることができる。
・h->Sumw2()とすると、このヒストグラムをDrawした時に自動的に統計誤差の誤差伝播をして、Drawしてくれる。また、ヒストグラムを他のマクロからGet()などして取ってきた時は、エラーは継承されないので、Get()した後にh1->Sumw2()みたいにすると、統計誤差の誤差伝播をしてくれる。
また、TH1::SetDefaultSumw2()をmain関数の中に入れておくことで、わざわざSumw2()しなくても、ヒストグラムの演算をした場合にきちんと誤差伝播の式に則って各binの誤差を計算してくれる。
※デフォルトでは、bin同士の値の演算を終えてから演算後の値(N)に対して√Nを計算するため、誤差を過小評価することになる。注意!!!
(ex.)A ± δa, B ± δbというデータがあるとする。
・Sumw2()をした場合：A-Bの誤差は、√(δa)^2+(δb)^2で計算されるが、
・Sumw2()をしない場合：A-Bの誤差は、(δa)^2-(δb)^2=δcとすると、√δcが誤差として計算される。
【グラフ(TGraph)編】
ここから、定義したグラフ名を"g"として話す。
・TGraphにおいて、TH1などのFindBin()やGetBinContent()に対応するものが当然存在し、それは、GetPointX()やGetPointY()に対応する。例えば、x座標が500のところの横軸と縦軸の値が知りたければ、g->GetPointX(500), g->GetPointY(500)で知ることができる。
・複数のヒストグラムを一つのキャンバスに描いた時に、どれがどれを表しているか人に見せる時など分かりにくい場合がある。そういう時にはTLegendを使うと便利。TLegendの使い方は以下の通り。
まず、TLegendを載せる位置を決める。位置は、図全体から見てデフォルトで左下を(0, 0)、右上を(1, 1)となっていて、TLegend *l = new TLegend(左下x, 左下y, 右上x, 右上y)で表す。TLegend(0.8, 0.68, 0.99, 0.78)くらいがちょうどいいかも。
次に、実際にヒストグラムをAddしていく。Addの仕方は次の通り。l->AddEntry(TH1Dなどで定義しているヒストグラムの名前, "表示させたい名前", "オプション")となる。オプションは、"f"=box表示、"l"=line表示、"p"=marker表示となっている。
最後に、l->Draw()としてやれば、Legendがキャンバス上に表示されるはず。
・長方形を描きたかったら、[適当に名前を定義] = new TBox(下限のx値,下限のy値,上限のx値,上限のy値)と書けばいい。ここでFillのオプションとして、[適当に名前を定義]->SetFillStyle(3001)とすると、長方形の色が薄くなる。SetFillColor()でFillする色も決められる。
・線を書くにはTLineを使う。具体的な定義の仕方は、TLine *l = new TLine(x軸の下限,y軸の下限,x軸の上限,y軸の上限)のように書く。
線の種類も変更することができ、デフォルトは普通の直線だが、l->SetLineStyle()でかっこの中に1~10の数字を入れることで変更可能。線の種類はこのページに書かれてある。
また、SetLineColor(適当な数字)で線の色も変更可能。上のリンクにはマーカーの種類も書いてある。
・右肩に出てる統計情報のボックスは、gStyle->SetOptStat(0)で非表示、gStyle->SetOptStat(1111111)で全ての情報を表示、gStyle->SetOptStat(0001111)でデフォルトになる。デフォルトはSetOptStat()でも可。
・キャンバスにおいて、c1->SetGrid()とすると、キャンバスにGridが表示される。x,y軸それぞれでやりたかったら、c1->SetGridx()、c1->SetGridy()とする。
Gridの線のタイプは、gStyle->SetGridStyle(1)など数字を変えることで変更可能。(デフォルト値は3)
・TF1にて、ROOTで定義済みの関数として、expo,gaus,polなどが与えられている。これらの関数の和を描きたいならば、expo(0)+gaus(2)+pol(5)とする。()内の数字は、パラメータの開始番号を表す。
ちなみに、expo：exp([0]+[1]*x), gaus：[0]*exp(-0.5*( (x-[1])/[2])**2), polN：[0]+[1]*x+[2]*x^2+…+[N]*x^Nである。
【TF1の表示について】
f1->Eval(x) //xの値を代入してyの値を取得
f1->Derivative(x) // xの微分を得る
f1->SetParameters(p0,p1,...) // パラメータの初期値を設定
f1->SetParNames("A","B",...) // パラメータの名前を設定
f1->FixParameter(N) // N番目のパラメータを固定する
f1->SetParLimits(N, min, max) // N番目のパラメータの可動範囲をmin~mixに制限
// フィッティングパラメータ関係の値を取得
f1->GetParameter(N) // N番目のパラメータを取得(※Nは0から)
f1->GetParameters(par) // パラメータを配列に入れる。この時まず、double par[3]などとして配列を定義しておく。
f1->GetParError(N) // N番目のパラメータのerrorを取得
f1->GetChisquare() // chi squareを取得
f1->GetNDF() // 自由度を取得