Remove obsolete drawing methods and macros
[u/mrichter/AliRoot.git] / STAT / AliTMinuitToolkit.cxx
CommitLineData
b7d1d891 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
4b23a517 15
b7d1d891 16/* $Id$ */
17
18#include <TCanvas.h>
4b23a517 19#include <TF1.h>
20#include <TFormula.h>
b7d1d891 21#include <TH1F.h>
22#include <TH2F.h>
4b23a517 23#include <TMath.h>
b7d1d891 24#include <TMatrix.h>
4b23a517 25#include <TRandom.h>
b7d1d891 26#include <TString.h>
27#include <TVector.h>
28#include <TVectorD.h>
29#include <TVirtualFitter.h>
4b23a517 30
b7d1d891 31#include "AliTMinuitToolkit.h"
4b23a517 32
dc697472 33//--------------------------------------------------------------------------------------
4b23a517 34//
35// The AliTMinuitToolkit serves as an easy to use interface for the TMinuit
36// package:
37//
38// - It allows to fit a curve to one and two dimensional histograms
39// (TH2F::Fit() only allows to fit a hyperplane).
40// - Or n points can be specified directly via a n x 2 matrix.
dc697472 41// - An option for robust fitting with non-linear functions is implemented.
305a4a95 42//
43// A small example illustrating the usage of AliTMinuitToolkit is given in the function
44// "AliTMinuitToolkit::Test()".
4b23a517 45//
46//
47// 1. Setting the formula:
48//
49// The formula is simply set via "void SetFitFunction(TFormula * formula)".
50//
51//
52// 2. Adding the data points
53//
54// - In order to fit a histogram, use "void FitHistogram(TH1F * his)" or
55// "void FitHistogram(TH2F * his)". The fitter is started automatically
56// - Alternatively, the direct specification of the points is possible via
57// "void SetPoints(TMatrixD * points)". Note, that the each point
58// corresponds to one row in the matrix. The fitter is then started with
dc697472 59// the command "void Fit()". The weight of each point can be specified
60// with an n-dimensional vector using "void SetWeights(TVectorD * weights)".
4b23a517 61//
62//
63// 3. Accessing the fit results
64//
65// The N parameters of the formula are stored in a N-dimensional vector which
dc697472 66// is returned by "TVectorD * GetParameters()". In a similar way the covariance
4b23a517 67// matrix of the fit is returned via "TMatrixD * GetCovarianceMatrix()" which
68// is of the type N x N.
69//
70//
71// 4. Non-linear robust fitting:
72//
73// Even a few outliers can lead to wrong results of a least-squares fitting
74// procedure. In this case the use of robust(resistant) methods can be
75// helpful, but a stronger dependence on starting values or convergence to
76// local minima can occur.
77//
dc697472 78// The robust option becomes active if EnableRobust(true, sigma) is called. It is
79// very much recommended that a normalization value (scale variable) corresponding
80// to an expected deviation (sigma) is specified via
81// "EnableRobust(Bool_t b, Double_t sigma)".
82//
83// Performing the fit without knowledge of sigma is also possible if only
84// "EnableRobust(true)" is activated, but this is NOT RECOMMENDED.
85//
86// The method is based on another estimator instead of chi^2. For small deviations
87// the function behaves like x^2 and for larger deviations like |x| - the so
88// called Huber estimator:
89//
90// h(x) = x^2 , for x < 2.5*sigma
91// h(x) = 2*(2.5*sigma)*x - (2.5*sigma)^2 , for x > 2.5*sigma
92//
93// If a weighting function is specified in addition, a second run with the ordinary
94// metric is started, but before entering the iteration every point is weighted
95// according to its distance to the outcoming function of the first run. The weighting
96// function w(x) must be defined on the intervall x in [0,1]. w(0) then
97// corresponds to the weight of the closest point and w(1) to the point with the
98// largest distance.
4b23a517 99//
100// Some standard weighting functions are predefined in
101// "SetWeightFunction(Char_t * name, Float_t param1, Float_t param2 = 0)":
102// - "BOX" equals to 1 if x < param1 and to 0 if x > param1.
103// - "EXPONENTIAL" corresponds to "Math::Exp(-TMath::Log(param1)*x)"
104// - "ERRORFUNCTION" corresponds to "TMath::Erfc((x-param1)/param2)"
105//
dc697472 106//
107// REFERENCE for non-linear robust fitting:
108// Ekblom H. and Madsen K. (1988), Alogrithms for non-linear Huber estimation,
109// BIT Numerical Mathematics 29 (1989) 60-76.
110// internet: http://www.springerlink.com/content/m277218542988344/
111//
112//
113// 5. examples:
114//
115// A small example illustrating the working principles of AliTMinuitToolkit is given
116// in the function "AliTMinuitToolkit::Test()".
117//
118//
119//
120// Comments and questions are always welcome: A.Kalweit@gsi.de
121//--------------------------------------------------------------------------------------
4b23a517 122
123
124ClassImp(AliTMinuitToolkit)
125
126AliTMinuitToolkit::AliTMinuitToolkit() :
127 TNamed(),
128 fFormula(0),
129 fWeightFunction(0),
bf7bd859 130 fFitAlgorithm(""),
4b23a517 131 fPoints(0),
dc697472 132 fWeights(0),
4b23a517 133 fParam(0),
134 fParamLimits(0),
135 fCovar(0),
da7c9ae3 136 fChi2(0),
137 fMaxCalls(0),
dc697472 138 fPrecision(0),
139 fUseRobust(0),
140 fExpectedSigma(0)
4b23a517 141{
142 //
143 // standard constructor
144 //
145 fMaxCalls = 500;
146 fPrecision = 1;
dc697472 147 fUseRobust = false;
148 fExpectedSigma = 0;
4b23a517 149}
150
151
da7c9ae3 152AliTMinuitToolkit::AliTMinuitToolkit(const AliTMinuitToolkit&) :
153 TNamed(),
154 fFormula(0),
155 fWeightFunction(0),
bf7bd859 156 fFitAlgorithm(""),
da7c9ae3 157 fPoints(0),
dc697472 158 fWeights(0),
da7c9ae3 159 fParam(0),
160 fParamLimits(0),
161 fCovar(0),
162 fChi2(0),
163 fMaxCalls(0),
dc697472 164 fPrecision(0),
165 fUseRobust(0),
166 fExpectedSigma(0)
da7c9ae3 167{
168
169
170}
171
172
173AliTMinuitToolkit& AliTMinuitToolkit::operator=(const AliTMinuitToolkit&) {
174
175 return *this;
176}
177
178
4b23a517 179
180AliTMinuitToolkit::~AliTMinuitToolkit(){
181 //
182 // destructor
183 //
184 delete fPoints;
dc697472 185 delete fWeights;
4b23a517 186 delete fWeightFunction;
187 delete fParamLimits;
188 delete fFormula;
189 delete fParam;
190 delete fCovar;
191 delete fChi2;
192}
193
b7d1d891 194void AliTMinuitToolkit::FitHistogram(TH1F *const his) {
4b23a517 195 //
196 // Fit a one dimensional histogram
197 //
198 fPoints = new TMatrixD(his->GetNbinsX(), 2);
199
200 for(Int_t ibin=0; ibin < his->GetNbinsX(); ibin++) {
201 Double_t x = his->GetXaxis()->GetBinCenter(ibin+1);
202 Double_t y = his->GetBinContent(ibin+1);
203
204 (*fPoints)(ibin, 0) = x;
205 (*fPoints)(ibin, 1) = y;
206 }
207
208 Fit();
209}
210
211
b7d1d891 212void AliTMinuitToolkit::FitHistogram(TH2F *const his) {
4b23a517 213 //
dc697472 214 // Fit a curve to a two dimensional histogram
4b23a517 215 //
216 fPoints = new TMatrixD((Long64_t)his->GetEntries(), 2);
217 Long64_t entry = 0;
218
219 for(Int_t ibin=0; ibin < his->GetNbinsX(); ibin++) {
220 Double_t x = his->GetXaxis()->GetBinCenter(ibin);
221 for(Int_t jbin=0; jbin < his->GetNbinsY(); jbin++) {
222 Long64_t n = his->GetBin(ibin, jbin);
223 Double_t y = his->GetYaxis()->GetBinCenter(jbin);
224 for(Int_t ientries=0; ientries < his->GetBinContent(n); ientries++) {
225 (*fPoints)(entry,0) = x;
226 (*fPoints)(entry,1) = y;
227 entry++;
228 }
229
230 }
231 }
232
233 Fit();
234}
235
236
b7d1d891 237void AliTMinuitToolkit::SetWeightFunction(const Char_t *name, Float_t param1, Float_t param2) {
4b23a517 238 //
239 // Set the weight function which must be defined on the interval [0,1].
240 //
241 TString FuncType(name);
242 FuncType.ToUpper();
243
244 if (FuncType == "EXPONENTIAL") fWeightFunction = new TFormula("exp", Form("TMath::Exp(-TMath::Log(%f)*x)", param1));
245 if (FuncType == "BOX") fWeightFunction = new TFormula("box", Form("TMath::Erfc((x-%f)/0.0001)", param1));
dc697472 246 if (FuncType == "ERRORFUNCTION") fWeightFunction = new TFormula("err", Form("TMath::Erfc((x-%f)/%f)", param1, param2));
4b23a517 247
248}
249
250
251void AliTMinuitToolkit::FitterFCN(int &npar, double *dummy, double &fchisq, double *gin, int iflag){
252 //
253 // internal function which gives the specified function to the TMinuit function
254 //
255
256 // suppress warnings for unused variables:
257 dummy = dummy;
258 iflag = iflag;
259 npar = npar;
260 //
261 AliTMinuitToolkit * fitter = (AliTMinuitToolkit*)TVirtualFitter::GetFitter()->GetObjectFit();
262 fchisq = 0;
263 Int_t nvar = fitter->GetPoints()->GetNcols()-1;
264 Int_t npoints = fitter->GetPoints()->GetNrows();
265
dc697472 266 // calculate mean deviation for normalization or use user-defined sigma
267 Double_t dev = 0.;
268 if (fitter->GetExpectedSigma() == 0 && fitter->GetStatus() == true) {
269 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
4b23a517 270 Double_t x[100];
dc697472 271 for (Int_t ivar=0; ivar<nvar; ivar++){
4b23a517 272 x[ivar] = (*fitter->GetPoints())(ipoint, ivar);
dc697472 273 }
274 Float_t funx = fitter->GetFormula()->EvalPar(x,gin);
275 Double_t delta = (*fitter->GetPoints())(ipoint, nvar) - funx;
276 dev += TMath::Sqrt(TMath::Abs(delta));
277 }
278 dev = dev/npoints;
279 } else {
280 dev = fitter->GetExpectedSigma();
4b23a517 281 }
dc697472 282 // calculate chisquare
4b23a517 283 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
284 Double_t x[100];
285 for (Int_t ivar=0; ivar<nvar; ivar++){
286 x[ivar] = (*fitter->GetPoints())(ipoint, ivar);
287 }
dc697472 288 Float_t funx = fitter->GetFormula()->EvalPar(x,gin);
289 Double_t delta = TMath::Abs((*fitter->GetPoints())(ipoint, nvar) - funx);
290 if (fitter->GetStatus() == true) {
291 delta = delta/dev; // normalization
292 if (delta <= 2.5) fchisq+= delta*delta; // new metric: Huber-k-estimator
293 if (delta > 2.5) fchisq+= 2*(2.5)*delta - (2.5*2.5);
294 } else {
295 Double_t weight = (*fitter->GetWeights())(ipoint);
296 fchisq+= delta*delta*weight; //old metric
297 }
4b23a517 298 }
dc697472 299 }
4b23a517 300
301
302void AliTMinuitToolkit::Fit() {
303 //
304 // internal function that calls the fitter
305 //
dc697472 306 Int_t nparam = fParam->GetNrows();
307 Int_t npoints = fPoints->GetNrows();
308 Int_t nvar = fPoints->GetNcols()-1;
4b23a517 309
310 // set all paramter limits to infinity as default
311 if (fParamLimits == 0) {
312 fParamLimits = new TMatrixD(nparam ,2);
313 for (Int_t iparam=0; iparam<nparam; iparam++){
314 (*fParamLimits)(iparam, 0) = 0;
315 (*fParamLimits)(iparam, 1) = 0;
316 }
317 }
318
319 // set all weights to 1 as default
dc697472 320 Bool_t weightFlag = false;
4b23a517 321 if (fWeightFunction == 0) {
322 fWeightFunction = new TFormula("constant", "1");
dc697472 323 } else {
324 weightFlag = true;
4b23a517 325 }
326
327 // migrad fit algorithm as default
bf7bd859 328 if (fFitAlgorithm == "") {
4b23a517 329 fFitAlgorithm = "migrad";
330 }
331
dc697472 332 // assign weights
333 if (fWeights == 0) {
334 fWeights = new TVectorD(npoints);
335 for (Int_t ipoint=0; ipoint<npoints; ipoint++) (*fWeights)(ipoint) = 1;
336 }
337
4b23a517 338 // set up the fitter
339 TVirtualFitter *minuit = TVirtualFitter::Fitter(0, nparam);
340 minuit->SetObjectFit(this);
dc697472 341 minuit->SetFCN(AliTMinuitToolkit::FitterFCN);
4b23a517 342
343 // initialize paramters (step size???)
344 for (Int_t iparam=0; iparam<nparam; iparam++){
345 minuit->SetParameter(iparam, Form("p[%d]",iparam), (*fParam)(iparam), (*fParam)(iparam)/10, (*fParamLimits)(iparam, 0), (*fParamLimits)(iparam, 1));
346 }
dc697472 347
348 //
4b23a517 349 Double_t argList[2];
350 argList[0] = fMaxCalls; //maximal number of calls
351 argList[1] = fPrecision; //tolerance normalized to 0.001
352 if (fMaxCalls == 500 && fPrecision == 1) minuit->ExecuteCommand(fFitAlgorithm, 0, 0);
353 if (fMaxCalls != 500 || fPrecision != 1) minuit->ExecuteCommand(fFitAlgorithm, argList, 2);
354 // two additional arguments can be specified ExecuteCommand("migrad", argList, 2) - use 0,0 for default
dc697472 355
356 // fill parameter vector
357 for (Int_t ivar=0; ivar<nparam; ivar++){
358 (*fParam)(ivar) = minuit->GetParameter(ivar);
359 fFormula->SetParameter(ivar, minuit->GetParameter(ivar));
360 }
4b23a517 361
dc697472 362 // if a weight function is specified -> enter 2nd run with weights
363 if (weightFlag == true && fUseRobust == true) {
364 // sort points for weighting
365 Double_t *sortList = new Double_t[npoints];
366 Int_t *indexList = new Int_t[npoints];
367 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
368 Double_t funx = fFormula->Eval((*fPoints)(ipoint, 0));
369 Double_t delta = TMath::Abs((*fPoints)[ipoint][nvar] - funx);
370 sortList[ipoint] = delta;
371 }
372 TMath::Sort(npoints, sortList, indexList, false);
373 for (Int_t ip=0; ip<npoints; ip++){
374 Double_t t = ip/(Double_t)npoints;
375 (*fWeights)(indexList[ip]) = fWeightFunction->Eval(t);
376 }
377
378 // set up the fitter
379 fUseRobust = false;
380 for (Int_t iparam=0; iparam<nparam; iparam++){
381 minuit->SetParameter(iparam, Form("p[%d]",iparam), (*fParam)(iparam), (*fParam)(iparam)/10, (*fParamLimits)(iparam, 0), (*fParamLimits)(iparam, 1));
382 }
383 // start fitting
384 if (fMaxCalls == 500 && fPrecision == 1) minuit->ExecuteCommand(fFitAlgorithm, 0, 0);
385 if (fMaxCalls != 500 || fPrecision != 1) minuit->ExecuteCommand(fFitAlgorithm, argList, 2);
386 fUseRobust = true;
387
4ce766eb 388 delete [] sortList;
389 delete [] indexList;
dc697472 390 }
391
4b23a517 392 // fill parameter vector
393 for (Int_t ivar=0; ivar<nparam; ivar++){
394 (*fParam)(ivar) = minuit->GetParameter(ivar);
dc697472 395 fFormula->SetParameter(ivar, minuit->GetParameter(ivar));
4b23a517 396 }
397
398 // fill covariance matrix
399 fCovar = new TMatrixD(nparam, nparam);
4b23a517 400 for(Int_t i=0; i < nparam; i++) {
401 for(Int_t j=0; j < nparam; j++) {
402 (*fCovar)(i,j) = minuit->GetCovarianceMatrixElement(i,j);
403 }
404 }
dc697472 405
406 if (weightFlag == false) fWeightFunction = 0;
4b23a517 407}
408
409
410
411void AliTMinuitToolkit::Test() {
412 //
413 // This test function shows the basic working principles of this class
414 // and illustrates how a robust fit can improve the results
415 //
305a4a95 416
417 // 1. provide some example histogram
dc697472 418 TH1F * hist = new TH1F("test", "with (red) and without (black) robust option", 20,0,4);
4b23a517 419 TRandom * rand = new TRandom();
420 for (Int_t i = 0; i < 10000; i++) {
421 hist->Fill(rand->Exp(1));
422 if (i < 1000) hist->Fill(3); //"outliers"
423 if (i < 1070) hist->Fill(3.5);
424 if (i < 670) hist->Fill(2);
425 if (i < 770) hist->Fill(1.5);//"outliers"
426 if (i < 740) hist->Fill(1);
427 }
428 TCanvas * canv = new TCanvas();
429 canv->cd(1);
430 hist->Draw();
dc697472 431
305a4a95 432 // 2. example fit without robust option
433 AliTMinuitToolkit * tool = new AliTMinuitToolkit();
b7d1d891 434 TFormula *aFormExp = new TFormula("formExp", "[0]*TMath::Exp(-[1]*x)");
435 tool->SetFitFunction(aFormExp);
305a4a95 436 TVectorD *vec1 = new TVectorD(2); // Set initial values
437 (*vec1)(0) = 1800;
438 (*vec1)(1) = 1;
439 tool->SetInitialParam(vec1);
440 tool->FitHistogram(hist);
441
4b23a517 442 // draw fit function
443 TF1 *func = new TF1("test", "[0]*TMath::Exp(-[1]*x)", 0, 6);
305a4a95 444 func->SetParameters((*tool->GetParameters())(0), (*tool->GetParameters())(1));
4b23a517 445 func->Draw("same");
dc697472 446
305a4a95 447 // 3 . robust fit
4b23a517 448 TVectorD *vec2 = new TVectorD(2);
449 (*vec2)(0) = 1800;
450 (*vec2)(1) = 1;
305a4a95 451 tool->SetInitialParam(vec2);
452 tool->EnableRobust(true, 10);
453 tool->SetWeightFunction("box", 0.75);
454 tool->FitHistogram(hist);
4b23a517 455 TF1 *func2 = new TF1("test2", "[0]*TMath::Exp(-[1]*x)", 0, 6);
305a4a95 456 func2->SetParameter(0, (*tool->GetParameters())(0));
457 func2->SetParameter(1, (*tool->GetParameters())(1));
4b23a517 458 func2->SetLineColor(kRed);
459 func2->Draw("same");
460
461}
462