rootのTProfileで任意の幅でrebinしたりProfileどうしを足したり割ったりする方法
いくつかのTProfile(or histogram)同士でrebinしたりProfile同士を足したり割ったりしたい。
今まではGetBinContent()でビンの中身持って来て足したりとか脳筋っぽいことをしてたんだけど、
histogramと違ってProfileの場合はエラーの計算とかがめんどくさい。
ヒストはエントリー数足してルートとればエラーになるけど、profileは誤差伝搬とか面倒だからね。
そんなわけで便利な関数がないか調べてみた。
1、TProfile同士を足す。
http://www-he.scphys.kyoto-u.ac.jp/member/n.kamo/wiki/doku.php?id=study:software:root:main
このページに書いてあった。
hist->Add(hist1,x);//各binの値=現在の値+x*hist1の値
重みとか気にしないならx=1がデフォルトなので引数一個だけ入れれば大丈夫。
例えば BeforeProfile[N] というTProfileのオブジェクトを10個足したい場合
TProfile *profile;
for(int i=0;i<N;i++){
if(i==0)profile = BeforeProfile[i];
else profile->Add(BeforeProfile[i]);
}
みたいな感じで書けばよい。最初のTProfileには何も入ってないのでfor文の一回目は普通に代入して
2回目からAddしてるのがミソ。最初にこれを忘れて見事にbreakした。初期化って大事だね。
2、Rebinする
次にRebinする。今回は等間隔じゃなくて自分で任意の幅を決めてrebinしたい。
これもさっきのページに書いてあった。
・等分割の場合
TH1D h_rebin* = h1->Rebin(n); //(1)
TH1D h_rebin* = h1->Rebin(n,"h_rebin"); //第二引数は新しく作るhistoの名前
・好きな幅で分割する
double xbins[n+1] = {};
TH1D h_rebin* = h1->Rebin(n,"h_rebin",xbins);
第3引数に配列をぶっこむことですきな場所でrebinできるらしい。
(1)はn個で1つのbinにするという意味。全体でn個のbinに分けるという意味ではないので注意。
逆に(2)のnは何個に分けるかという意味なのでややこしい。配列の中身がn+1になってるのは最後の値の上限値があるから。
つまり0-100を2分割したかったら{0,50,100}と3つ値が必要ということ。
(2)の場合、TProfileをRebinしても出力はTH1の形式になるので注意。
3、TProfile同士で割る
これもさっきのAddの代わりにDivide()という関数を使ったら値は出たのだが、エラーの値が超でかい。。
よくみるとWarningが出ていて
"WARNING!!: The algorithm in TProfile::Divide computing the errors is not accurate
Instead of Divide(TProfile *h1, TProfile *h2, do:
TH1D *p1 = h1->ProjectionX();
TH1D *p2 = h2->ProjectionX();
p1->Divide(p2);
とか出てる。Profile::Divideだと正しくエラー計算できないってこと?どういうこっちゃ。
もともとのソースコード見るとコメントに
//THE ALGORITHM COMPUTING THE ERRORS IS WRONG. HELP REQUIRED
help requiredって…
どうやらDivide()は使い物にならなそう。GetBinContentとかでビンの値とエラーを持って来て
誤差伝搬とかで計算するしかなさそう。
結論、TProfileでもAddしたりRebinしたりする関数はあるが、Divideは使いものにならない。掛け算とかする関数もあるみたいだが、使ってないので使えるかわかんない。