2 #include "AliTMinuitToolkit.h"
4 #include <TVirtualFitter.h>
18 //--------------------------------------------------------------------------------------
20 // The AliTMinuitToolkit serves as an easy to use interface for the TMinuit
23 // - It allows to fit a curve to one and two dimensional histograms
24 // (TH2F::Fit() only allows to fit a hyperplane).
25 // - Or n points can be specified directly via a n x 2 matrix.
26 // - An option for robust fitting with non-linear functions is implemented.
28 // A small example illustrating the usage of AliTMinuitToolkit is given in the function
29 // "AliTMinuitToolkit::Test()".
32 // 1. Setting the formula:
34 // The formula is simply set via "void SetFitFunction(TFormula * formula)".
37 // 2. Adding the data points
39 // - In order to fit a histogram, use "void FitHistogram(TH1F * his)" or
40 // "void FitHistogram(TH2F * his)". The fitter is started automatically
41 // - Alternatively, the direct specification of the points is possible via
42 // "void SetPoints(TMatrixD * points)". Note, that the each point
43 // corresponds to one row in the matrix. The fitter is then started with
44 // the command "void Fit()". The weight of each point can be specified
45 // with an n-dimensional vector using "void SetWeights(TVectorD * weights)".
48 // 3. Accessing the fit results
50 // The N parameters of the formula are stored in a N-dimensional vector which
51 // is returned by "TVectorD * GetParameters()". In a similar way the covariance
52 // matrix of the fit is returned via "TMatrixD * GetCovarianceMatrix()" which
53 // is of the type N x N.
56 // 4. Non-linear robust fitting:
58 // Even a few outliers can lead to wrong results of a least-squares fitting
59 // procedure. In this case the use of robust(resistant) methods can be
60 // helpful, but a stronger dependence on starting values or convergence to
61 // local minima can occur.
63 // The robust option becomes active if EnableRobust(true, sigma) is called. It is
64 // very much recommended that a normalization value (scale variable) corresponding
65 // to an expected deviation (sigma) is specified via
66 // "EnableRobust(Bool_t b, Double_t sigma)".
68 // Performing the fit without knowledge of sigma is also possible if only
69 // "EnableRobust(true)" is activated, but this is NOT RECOMMENDED.
71 // The method is based on another estimator instead of chi^2. For small deviations
72 // the function behaves like x^2 and for larger deviations like |x| - the so
73 // called Huber estimator:
75 // h(x) = x^2 , for x < 2.5*sigma
76 // h(x) = 2*(2.5*sigma)*x - (2.5*sigma)^2 , for x > 2.5*sigma
78 // If a weighting function is specified in addition, a second run with the ordinary
79 // metric is started, but before entering the iteration every point is weighted
80 // according to its distance to the outcoming function of the first run. The weighting
81 // function w(x) must be defined on the intervall x in [0,1]. w(0) then
82 // corresponds to the weight of the closest point and w(1) to the point with the
85 // Some standard weighting functions are predefined in
86 // "SetWeightFunction(Char_t * name, Float_t param1, Float_t param2 = 0)":
87 // - "BOX" equals to 1 if x < param1 and to 0 if x > param1.
88 // - "EXPONENTIAL" corresponds to "Math::Exp(-TMath::Log(param1)*x)"
89 // - "ERRORFUNCTION" corresponds to "TMath::Erfc((x-param1)/param2)"
92 // REFERENCE for non-linear robust fitting:
93 // Ekblom H. and Madsen K. (1988), Alogrithms for non-linear Huber estimation,
94 // BIT Numerical Mathematics 29 (1989) 60-76.
95 // internet: http://www.springerlink.com/content/m277218542988344/
100 // A small example illustrating the working principles of AliTMinuitToolkit is given
101 // in the function "AliTMinuitToolkit::Test()".
105 // Comments and questions are always welcome: A.Kalweit@gsi.de
106 //--------------------------------------------------------------------------------------
109 ClassImp(AliTMinuitToolkit)
111 AliTMinuitToolkit::AliTMinuitToolkit() :
128 // standard constructor
137 AliTMinuitToolkit::AliTMinuitToolkit(const AliTMinuitToolkit&) :
158 AliTMinuitToolkit& AliTMinuitToolkit::operator=(const AliTMinuitToolkit&) {
165 AliTMinuitToolkit::~AliTMinuitToolkit(){
171 delete fWeightFunction;
179 void AliTMinuitToolkit::FitHistogram(TH1F * his) {
181 // Fit a one dimensional histogram
183 fPoints = new TMatrixD(his->GetNbinsX(), 2);
185 for(Int_t ibin=0; ibin < his->GetNbinsX(); ibin++) {
186 Double_t x = his->GetXaxis()->GetBinCenter(ibin+1);
187 Double_t y = his->GetBinContent(ibin+1);
189 (*fPoints)(ibin, 0) = x;
190 (*fPoints)(ibin, 1) = y;
197 void AliTMinuitToolkit::FitHistogram(TH2F * his) {
199 // Fit a curve to a two dimensional histogram
201 fPoints = new TMatrixD((Long64_t)his->GetEntries(), 2);
204 for(Int_t ibin=0; ibin < his->GetNbinsX(); ibin++) {
205 Double_t x = his->GetXaxis()->GetBinCenter(ibin);
206 for(Int_t jbin=0; jbin < his->GetNbinsY(); jbin++) {
207 Long64_t n = his->GetBin(ibin, jbin);
208 Double_t y = his->GetYaxis()->GetBinCenter(jbin);
209 for(Int_t ientries=0; ientries < his->GetBinContent(n); ientries++) {
210 (*fPoints)(entry,0) = x;
211 (*fPoints)(entry,1) = y;
222 void AliTMinuitToolkit::SetWeightFunction(Char_t * name, Float_t param1, Float_t param2) {
224 // Set the weight function which must be defined on the interval [0,1].
226 TString FuncType(name);
229 if (FuncType == "EXPONENTIAL") fWeightFunction = new TFormula("exp", Form("TMath::Exp(-TMath::Log(%f)*x)", param1));
230 if (FuncType == "BOX") fWeightFunction = new TFormula("box", Form("TMath::Erfc((x-%f)/0.0001)", param1));
231 if (FuncType == "ERRORFUNCTION") fWeightFunction = new TFormula("err", Form("TMath::Erfc((x-%f)/%f)", param1, param2));
236 void AliTMinuitToolkit::FitterFCN(int &npar, double *dummy, double &fchisq, double *gin, int iflag){
238 // internal function which gives the specified function to the TMinuit function
241 // suppress warnings for unused variables:
246 AliTMinuitToolkit * fitter = (AliTMinuitToolkit*)TVirtualFitter::GetFitter()->GetObjectFit();
248 Int_t nvar = fitter->GetPoints()->GetNcols()-1;
249 Int_t npoints = fitter->GetPoints()->GetNrows();
251 // calculate mean deviation for normalization or use user-defined sigma
253 if (fitter->GetExpectedSigma() == 0 && fitter->GetStatus() == true) {
254 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
256 for (Int_t ivar=0; ivar<nvar; ivar++){
257 x[ivar] = (*fitter->GetPoints())(ipoint, ivar);
259 Float_t funx = fitter->GetFormula()->EvalPar(x,gin);
260 Double_t delta = (*fitter->GetPoints())(ipoint, nvar) - funx;
261 dev += TMath::Sqrt(TMath::Abs(delta));
265 dev = fitter->GetExpectedSigma();
267 // calculate chisquare
268 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
270 for (Int_t ivar=0; ivar<nvar; ivar++){
271 x[ivar] = (*fitter->GetPoints())(ipoint, ivar);
273 Float_t funx = fitter->GetFormula()->EvalPar(x,gin);
274 Double_t delta = TMath::Abs((*fitter->GetPoints())(ipoint, nvar) - funx);
275 if (fitter->GetStatus() == true) {
276 delta = delta/dev; // normalization
277 if (delta <= 2.5) fchisq+= delta*delta; // new metric: Huber-k-estimator
278 if (delta > 2.5) fchisq+= 2*(2.5)*delta - (2.5*2.5);
280 Double_t weight = (*fitter->GetWeights())(ipoint);
281 fchisq+= delta*delta*weight; //old metric
287 void AliTMinuitToolkit::Fit() {
289 // internal function that calls the fitter
291 Int_t nparam = fParam->GetNrows();
292 Int_t npoints = fPoints->GetNrows();
293 Int_t nvar = fPoints->GetNcols()-1;
295 // set all paramter limits to infinity as default
296 if (fParamLimits == 0) {
297 fParamLimits = new TMatrixD(nparam ,2);
298 for (Int_t iparam=0; iparam<nparam; iparam++){
299 (*fParamLimits)(iparam, 0) = 0;
300 (*fParamLimits)(iparam, 1) = 0;
304 // set all weights to 1 as default
305 Bool_t weightFlag = false;
306 if (fWeightFunction == 0) {
307 fWeightFunction = new TFormula("constant", "1");
312 // migrad fit algorithm as default
313 if (fFitAlgorithm == 0) {
314 fFitAlgorithm = "migrad";
319 fWeights = new TVectorD(npoints);
320 for (Int_t ipoint=0; ipoint<npoints; ipoint++) (*fWeights)(ipoint) = 1;
324 TVirtualFitter *minuit = TVirtualFitter::Fitter(0, nparam);
325 minuit->SetObjectFit(this);
326 minuit->SetFCN(AliTMinuitToolkit::FitterFCN);
328 // initialize paramters (step size???)
329 for (Int_t iparam=0; iparam<nparam; iparam++){
330 minuit->SetParameter(iparam, Form("p[%d]",iparam), (*fParam)(iparam), (*fParam)(iparam)/10, (*fParamLimits)(iparam, 0), (*fParamLimits)(iparam, 1));
335 argList[0] = fMaxCalls; //maximal number of calls
336 argList[1] = fPrecision; //tolerance normalized to 0.001
337 if (fMaxCalls == 500 && fPrecision == 1) minuit->ExecuteCommand(fFitAlgorithm, 0, 0);
338 if (fMaxCalls != 500 || fPrecision != 1) minuit->ExecuteCommand(fFitAlgorithm, argList, 2);
339 // two additional arguments can be specified ExecuteCommand("migrad", argList, 2) - use 0,0 for default
341 // fill parameter vector
342 for (Int_t ivar=0; ivar<nparam; ivar++){
343 (*fParam)(ivar) = minuit->GetParameter(ivar);
344 fFormula->SetParameter(ivar, minuit->GetParameter(ivar));
347 // if a weight function is specified -> enter 2nd run with weights
348 if (weightFlag == true && fUseRobust == true) {
349 // sort points for weighting
350 Double_t *sortList = new Double_t[npoints];
351 Int_t *indexList = new Int_t[npoints];
352 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
353 Double_t funx = fFormula->Eval((*fPoints)(ipoint, 0));
354 Double_t delta = TMath::Abs((*fPoints)[ipoint][nvar] - funx);
355 sortList[ipoint] = delta;
357 TMath::Sort(npoints, sortList, indexList, false);
358 for (Int_t ip=0; ip<npoints; ip++){
359 Double_t t = ip/(Double_t)npoints;
360 (*fWeights)(indexList[ip]) = fWeightFunction->Eval(t);
365 for (Int_t iparam=0; iparam<nparam; iparam++){
366 minuit->SetParameter(iparam, Form("p[%d]",iparam), (*fParam)(iparam), (*fParam)(iparam)/10, (*fParamLimits)(iparam, 0), (*fParamLimits)(iparam, 1));
369 if (fMaxCalls == 500 && fPrecision == 1) minuit->ExecuteCommand(fFitAlgorithm, 0, 0);
370 if (fMaxCalls != 500 || fPrecision != 1) minuit->ExecuteCommand(fFitAlgorithm, argList, 2);
377 // fill parameter vector
378 for (Int_t ivar=0; ivar<nparam; ivar++){
379 (*fParam)(ivar) = minuit->GetParameter(ivar);
380 fFormula->SetParameter(ivar, minuit->GetParameter(ivar));
383 // fill covariance matrix
384 fCovar = new TMatrixD(nparam, nparam);
385 for(Int_t i=0; i < nparam; i++) {
386 for(Int_t j=0; j < nparam; j++) {
387 (*fCovar)(i,j) = minuit->GetCovarianceMatrixElement(i,j);
391 if (weightFlag == false) fWeightFunction = 0;
396 void AliTMinuitToolkit::Test() {
398 // This test function shows the basic working principles of this class
399 // and illustrates how a robust fit can improve the results
402 // 1. provide some example histogram
403 TH1F * hist = new TH1F("test", "with (red) and without (black) robust option", 20,0,4);
404 TRandom * rand = new TRandom();
405 for (Int_t i = 0; i < 10000; i++) {
406 hist->Fill(rand->Exp(1));
407 if (i < 1000) hist->Fill(3); //"outliers"
408 if (i < 1070) hist->Fill(3.5);
409 if (i < 670) hist->Fill(2);
410 if (i < 770) hist->Fill(1.5);//"outliers"
411 if (i < 740) hist->Fill(1);
413 TCanvas * canv = new TCanvas();
417 // 2. example fit without robust option
418 AliTMinuitToolkit * tool = new AliTMinuitToolkit();
419 TFormula *FormExp = new TFormula("formExp", "[0]*TMath::Exp(-[1]*x)");
420 tool->SetFitFunction(FormExp);
421 TVectorD *vec1 = new TVectorD(2); // Set initial values
424 tool->SetInitialParam(vec1);
425 tool->FitHistogram(hist);
428 TF1 *func = new TF1("test", "[0]*TMath::Exp(-[1]*x)", 0, 6);
429 func->SetParameters((*tool->GetParameters())(0), (*tool->GetParameters())(1));
433 TVectorD *vec2 = new TVectorD(2);
436 tool->SetInitialParam(vec2);
437 tool->EnableRobust(true, 10);
438 tool->SetWeightFunction("box", 0.75);
439 tool->FitHistogram(hist);
440 TF1 *func2 = new TF1("test2", "[0]*TMath::Exp(-[1]*x)", 0, 6);
441 func2->SetParameter(0, (*tool->GetParameters())(0));
442 func2->SetParameter(1, (*tool->GetParameters())(1));
443 func2->SetLineColor(kRed);