物理と数学とITのあれこれ

物理や数学、パソコンの設定とかプログラミングとかいろいろ書きなぐってます。

rootのTProfileで任意の幅でrebinしたりProfileどうしを足したり割ったりする方法

いくつかのTProfileor 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-1002分割したかったら{0,50,100}3つ値が必要ということ。

(2)の場合、TProfileRebinしても出力は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は使いものにならない。掛け算とかする関数もあるみたいだが、使ってないので使えるかわかんない。