ROOT/TMinuitでエネルギー較正のパラメーターを最適化する

TMinuitでエネルギー較正のパラメーターを最適化する†

TMinuit†

• https://root.cern.ch/doc/master/classTMinuit.html

参考サイト†

• https://qiita.com/pn11/items/36af3e3067f5d11f18eb
• https://qiita.com/kjyam23/items/90d3e7bb584090123c86
• https://llr.in2p3.fr/activites/physique/glast/workbook/pages/advanced_GRUG/GRUGminuit.pdf

やろうとしてること†

• 152Eu線源のガンマ線ピークのうち7本を使ってGe検出器のエネルギー較正を行う。生データでのチャンネル番号をx, ガンマ線エネルギーをyとして、

  y = ax + b

  の関係が成り立つとして、最適なa, bを決定する。

最小化する関数の定義†

• SetFCNで指定する。SetFCNは受け付ける関数の形式が決まっていて、

  void(*fcn)(Int_t &, Double_t *, Double_t &f, Double_t *, Int_t)

  という形式の静的外部関数(クラスに属するメンバー関数不可)を要求する。上記の参考サイトの見本を使うと

  void FCN(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t flag)

  になる。で、それぞれの変数の意味は

  npar : 最小化に含まれる自由パラメーターの数
  gin  : 計算された勾配値 (オプション)
  f    : 関数値そのもの
  par  : 定数および変数パラメーターのvector
  flag : FCN のいくつかのアクションを切り替える

  らしい。よくわからんけど、fとparだけ使えば最低限のことは出来る。

• で、最適なa, b を求める流れとしては、
  • まず適当にa,bを決める
  • 測定で得られた1本目のガンマ線ピークに対しy = ax + b を計算する
  • 計算で求めたyと真のyの差の二乗を計算する
  • 2本目以降も計算で求めたyと真のyの差の二乗を計算し、二乗値を加算していく
  • 最終的な二乗和をfに代入する
  • あとはTMinuitが勝手にfが最小になるa,bを決めてくれる(はず)
• FCNの例

  void FCN(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t flag){
  	double a,b;
  	float realpeak[7] = {121.7817, 244.6975, 344.2785, 778.9040, 867.378, 964.079, 1408.006};
  	float rawpeak[7]  = {259.4612, 521.1042, 733.0547, 1668.473, 1856.93, 2062.59, 3013.292};
  	double rss = 0;
  	a = par[0];
  	b = par[1];
  	for (int i=0; i<7; i++){
  		double y_calced = a * rawpeak[i] + b;
  		rss += (ealpeak[i] - y_calced) * (ealpeak[i] - y_calced);
  	}
  	f = rss;
  }

  • この例では、最終計算結果fが0に最も近くなるようにpar[0], par[1]を振って最適な値を決めてくれる。

• mainの例

  int main(){
  	TMinuit *min = new TMinuit(2); // パラメータ数
  	min->SetPrintLevel(-1); // 途中経過を非表示 (-1:非表示, 0:標準, 1:冗長表示)
  	min->SetFCN(FCN); // 関数を設定
  	min->DefineParameter(0, "a", 0.4, 0.0000001, 0.2, 0.5); // パラメータ初期値と刻み幅を設定
  	min->DefineParameter(1, "b", 4.4, 0.001, -10, +10);     // 第二引数の"a"とか"b"はnameを指定する。重複しないようにすれば適当でOK。
  	min->Migrad(); // 最小化処理を実行
  	double a,b,ea,eb;
  	min->GetParameter(0, a, dummy); // パラメータを取得
  	min->GetParameter(1, b, dummy);
  	printf("a : %f +- %f, b : %f +- %f \n", a, ea, b, eb);
  	return 0;
  }

• コンパイル

  g++ sample.cxx -o sample `root-config --cflags --libs` -lMinuit

• 注意
  • FCNには一般的な引数は設定できない。したがって、何らかの変数の引き渡しを行う場合はグローバル変数を使うしかない。

Topへ戻る