今年のGWはMinuit2で遊んでいました。自分なりに理解してきちんと動いたので、まとめてみます。

Minuit2とは?

Minuit2とは、元々FORTRANで開発されていたMinuitをC++用に開発しなおしたものです(ROOTに入っている従来のMinuitは単にFORTRANコードの移植であり、オブジェクト化されていない)。 Fitterとして普通に使う分には問題ない(TH1::FitはデフォルトとしてTMinuitを使っています)ですが、自分でオブジェクト指向プログラミングで書く場合(自分はARTEMISに組み込むために使いました)や、並列計算をする時にはMinuit2を使う必要があります。 また、tutorialに入っているtutorials/fit/minuit2FitBench.Cを走らせてみると分かりますが、Minuit2の方がコンパイラがより効率良く最適化してくれているからか、若干早いです。

因みにMinuit2をROOTで使用するのには、Buildするときにオプション--enable-minuit2をつけなければなりません。ご注意ください。

使い方

Minuit2には1つのページになっているIntroductionがあります。ここにサンプルコードもあります。また、MinuitのページにはMinuit2単体も配布されていて、この中にもサンプルコードが入っています。このpostでは、これを日本語で噛み砕いてまとめようと思います。

前提

今回は、横軸が波形で縦軸が波高値となるFlash ADCの波形データ(オシロスコープの出力と思ってください)に対してFittingするということをやります。波形は簡単のため正規分布に従っているとします。 そして、最後にFitしたあとの関数と波形データを重ねてDrawするようにします。

Fitting functionの定義

まず、Fitする関数のクラスを作ります。この関数は横軸の値を引数にし、縦軸の値を返す関数を作ります。

FittingFunction.h

#ifndef FITTING_FUNCTION_H
#define FITTING_FUNCTION_H

class FittingFunction {
public:
	FittingFunction();
	FittingFunction(Double_t,Double_t,Double_t);
	~FittingFunction();
	Double_t operator()(Double_t) const;
	Double_t operator()(Double_t*, Double_t*) const;
	
	Double_t GetConstant() const { return fConstant; }
	Double_t GetMean()     const { return fMean; }
	Double_t GetSigma()    const { return fSigma; }

private:
	Double_t fConstant;
	Double_t fMean;
	Double_t fSigma;

};

#endif

FittingFunction.cc

#include "FittingFunction.h"

FittingFunction::FittingFunction()
	: fConstant(kInvalidD), fMean(kInvalidD), fSigma(kInvalidD)
{
}

FittingFunction::FittingFunction(Double_t constant, Double_t mean, Double_t sigma)
	: fConstant(constant), fMean(mean), fSigma(sigma)
{
}

FittingFunction::~FittingFunction()
{
}

Double_t FittingFunction::operator()(Double_t timing) const
{
	// この関数オブジェクトはMinuit2で関数を使う時に必要
	// 引数が横軸の値(ここではtiming)に対応します
	Double_t val;
	val = TMath::Exp(-TMath::Power(timing - GetMean(),2.)/2./TMath::Power(GetSigma(),2.))/TMath::Sqrt(2.*TMath::Pi())/GetSigma();
	return val;
}

Double_t FittingFunction::operator()(Double_t *x, Double_t *p) const
{
	// この関数オブジェクトはTF1のオブジェクトを作るのに必要
	// x[0]がここで言うtiming, pはFitting parameters (constant, mean, sigma)に対応します
	Double_t val;
	val = TMath::Exp(-TMath::Power(x[0] - p[1],2.)/2./TMath::Power(p[2],2.))/TMath::Sqrt(2.*TMath::Pi())/p[2];
	return val;

Fitting functionの定義でのポイントはDouble_t型の引数を持つoperatorを作ることです。これは横軸の値を引数とし、パラメータはメンバ変数にアクセスすることにより戻り値を返す関数オブジェクトになります。 もう一つの関数オブジェクトに関しては後ほど述べます。

最小化させたい関数(FCN)の定義

続いて、Minuitで最小化させたい関数であるFCNの定義するクラスを作ります。ここで大体はminimum chi2やらLog-likelihoodを定義することになります。 ここでは、最小二乗法を計算するためのchi2関数を作ることにします。

FCNはROOT::Minuit2::FCNBaseを継承します。

FittingFCN.h

#ifndef FITTING_FCN_H
#define FITTING_FCN_H

#include "Minuit2/FCNBase.h"
#include <vector>

class FittingFCN : public ROOT::Minuit2::FCNBase {
public:
	FittingFCN();
	~FittingFCN();
	FittingFCN(const std::vector<Double_t>&, const std::vector<Double_t>&, const std::vector<Double_t>&);
	
	virtual Double_t operator()(const std::vector<Double_t>&) const;
	virtual Double_t Up() const { return fErrorDef; }

	std::vector<Double_t> GetSample()    const { return fSample; }
	std::vector<Double_t> GetClock()     const { return fClock; }
	std::vector<Double_t> GetMVariance() const { return fMVariance; }

	void SetErrorDef(Double_t val) { fErrorDef = val; }
	
private:

	std::vector<Double_t> fSample;
	std::vector<Double_t> fClock;
	std::vector<Double_t> fMVariance;
	Double_t fErrorDef;
};

#endif

FittingFCN.cc

#include "FittingFCN.h"
#include "FittingFunction.h"
#include "TMath.h"
#include <cassert>


FittingFCN::FittingFCN()
{
}

FittingFCN::~FittingFCN()
{
}

FittingFCN::FittingFCN(const std::vector<Double_t>& sample,
                       const std::vector<Double_t>& clock,
					   const std::vector<Double_t>& mvar)
	: fSample(sample), fClock(clock), fMVariance(mvar), fErrorDef(1.)
{
}

Double_t FittingFCN::operator()(const std::vector<Double_t>& par) const
{
	assert(par.size() == 3); // fitting parameterの数が3以外だとプログラムを落とす
	Double_t chi2 = 0.;
	FittingFunction func(par[0], par[1], par[2]);
	for (UInt_t iData = 0; iData != fSample.size(); ++iData) {
		chi2 += TMath::Power(func(fClock[iData]) - fSample[iData],2.)/fMVariance[iData];
	}
	return chi2;
}

今回扱うデータは波形データを想定しているので、変数名は分かりやすく横軸がclock, 縦軸はsampleとしました。mvarは1点での誤差を表しています。fErrorDefは何シグマでの誤差を評価するかを指定するのに使用する変数です。

ここではFittingFunction(Double_t constant, Double_t mean, Double_t sigma)のコンストラクタを用いてそれぞれFitting parameterを関数に入れて、そのパラメータを取ったときの、ある時刻での波高値の計算値$f(t_i)$とデータでの波高値$y_i$との残渣二乗和を取っていることが分かるかと思います。FCNの定義で重要なのは、Fitting parameterを引数に取る関数オブジェクトを作ることです。このオペレーター内に残渣二乗和を計算するコードを書きます。

オペレーションするコード

そして最後に、mainとなる関数です。

FittingProcessor.cc

#include "FittingFunction.h"
#include "FittingFCN.h"
#include "Minuit2/FuncitonMinimum.h"
#include "Minuit2/MnMigrad.h"
#include "Minuit2/MnPrint.h"
#include "Minuit2/MnUserParameterState.h"
#include <iostream>
#include <vector>
#include "TH1D.h"
#include "TRandom.h"

int main(void) 
{
	// 疑似波形データを生成
	TH1D hist("pulse","pulse shape",500,0,500);
	const UInt_t nData = 1e5;
	for (UInt_t iData = 0; iData != nData; ++iData) {
		hist.Fill(gRandom->Gaus(250,100));
	}
	// 横軸(clock),縦軸(sample),縦軸誤差(mvar)を定義し、波形データをvectorに詰める
	std::vector<Double_t> fSample;
	std::vector<Double_t> fClock;
	std::vector<Double_t> fSampleError;
	for (UInt_t iBin = 0; iBin != GetNbinsX(); ++iBin) {
		fClock.push_back(hist.GetBinCenter(iBin));
		fSample.push_back(hist.GetBinContent(iBin));
		fSampleError.push_back(10.);
	}
	// 初期値を入れておく
	Double_t peakY = hist.GetMaximum();
	Double_t peakT = hist.GetBinCenter(hist.GetMaximumBin());
	Double_t minT  = hist.GetBinCenter(hist.FindFirstBinAbove());
	Double_t maxT  = hist.GetBinCenter(hist.FindLastBinAbove());

	// Minuit2計算し、出力する
	FittingFCN fFCN(fSample, fClock, fSampleError);
	ROOT::Minuit2::MnUserParameters upar;
	upar.Add("Constant",peakY,10,0,1e6); // 引数は(変数名, 初期値, ステップ幅, 最小値, 最大値)
	upar.Add("Mean",peakT,0.1,minT,maxT);
	upar.Add("Sigma",10,0.1,0.001,1000);
	ROOT::Minuit2::MnMigrad migrad(fFCN, upar);
	ROOT::Minuit2::FunctionMinimum min = migrad(); // このタイミングで最小化を実行
	ROOT::Minuit2::MnUserParameterState state = min.UserState(); // Fitting resultを拾うためのオブジェクト
	std::cout << min << std::endl; // MnPrint.hをincludeすれば演算子 '<<' をoverloadしてくれるので、Printできるようになる
	// TF1型のFunctionを作り、得られたFitting parameterを取得し、関数のparameterに代入する
	const UInt_t nPrm = 3;
	FittingFunction fobj;
	TF1 func("func",fobj,minT,maxT,nPrm);

	Double_t fit_prm[nPrm], fit_prm_err[nPrm];
	for (UInt_t iPrm = 0; iPrm != nPrm; ++iPrm) {
		fit_prm[iPrm] = state.Value(iPrm);
		fit_prm_err[iPrm] = state.Error(iPrm);
		func.SetParameter(iPrm,fit_prm[iPrm]);
		}
	hist.SetTitle(TString::Format("chi2 = %f",min.Fval())); // 残渣(最小化後のFCNの値)はmin.Fval()とすると取れる
	hist.Draw();
	func.Draw("samel");
	
}

Minuit2のnamespaceはROOT::Minuit2となっております。Migradを使う場合はMnMigradを定義し、FunctionMinimumを呼びます。 振るParameterはMnUserParameters::Add(char*, double, double, double, double)を使用します。MnUserParametersは変数名で変数を引けるように作られています。Parameter番号はAddを呼ぶ順番によって決まっております。

また、FunctionのDrawingに関して補足しておきます。これはTF1のClass Refferenceのgeneral C++ function objectのところに書かれている方法を応用したものです。 TF1のFunction objectにはDouble_t operator()(Double_t *x, Double_t *p)のoperator(関数オブジェクト)を作れば、これを使ってくれるコンストラクタが用意されています。そのため、このオペレーターを予めFittingFunction.ccにて定義しておきました。

Minos errorを評価する、ParameterをFixするなど色々と出来ることがありますが、一先ずMinimumとしてはこれくらいで使えると思います。 詳しくはReferenceを見てください。