Bringing CMakeLists under svn maintenance
[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
a74e0f03 251void AliTMinuitToolkit::FitterFCN(int &/*npar*/, double */*dummy*/, double &fchisq, double *gin, int /*iflag*/){
4b23a517 252 //
253 // internal function which gives the specified function to the TMinuit function
254 //
255
4b23a517 256 //
257 AliTMinuitToolkit * fitter = (AliTMinuitToolkit*)TVirtualFitter::GetFitter()->GetObjectFit();
258 fchisq = 0;
259 Int_t nvar = fitter->GetPoints()->GetNcols()-1;
260 Int_t npoints = fitter->GetPoints()->GetNrows();
261
dc697472 262 // calculate mean deviation for normalization or use user-defined sigma
263 Double_t dev = 0.;
264 if (fitter->GetExpectedSigma() == 0 && fitter->GetStatus() == true) {
265 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
4b23a517 266 Double_t x[100];
dc697472 267 for (Int_t ivar=0; ivar<nvar; ivar++){
4b23a517 268 x[ivar] = (*fitter->GetPoints())(ipoint, ivar);
dc697472 269 }
270 Float_t funx = fitter->GetFormula()->EvalPar(x,gin);
271 Double_t delta = (*fitter->GetPoints())(ipoint, nvar) - funx;
272 dev += TMath::Sqrt(TMath::Abs(delta));
273 }
274 dev = dev/npoints;
275 } else {
276 dev = fitter->GetExpectedSigma();
4b23a517 277 }
dc697472 278 // calculate chisquare
4b23a517 279 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
280 Double_t x[100];
281 for (Int_t ivar=0; ivar<nvar; ivar++){
282 x[ivar] = (*fitter->GetPoints())(ipoint, ivar);
283 }
dc697472 284 Float_t funx = fitter->GetFormula()->EvalPar(x,gin);
285 Double_t delta = TMath::Abs((*fitter->GetPoints())(ipoint, nvar) - funx);
286 if (fitter->GetStatus() == true) {
287 delta = delta/dev; // normalization
288 if (delta <= 2.5) fchisq+= delta*delta; // new metric: Huber-k-estimator
289 if (delta > 2.5) fchisq+= 2*(2.5)*delta - (2.5*2.5);
290 } else {
291 Double_t weight = (*fitter->GetWeights())(ipoint);
292 fchisq+= delta*delta*weight; //old metric
293 }
4b23a517 294 }
dc697472 295 }
4b23a517 296
297
298void AliTMinuitToolkit::Fit() {
299 //
300 // internal function that calls the fitter
301 //
dc697472 302 Int_t nparam = fParam->GetNrows();
303 Int_t npoints = fPoints->GetNrows();
304 Int_t nvar = fPoints->GetNcols()-1;
4b23a517 305
306 // set all paramter limits to infinity as default
307 if (fParamLimits == 0) {
308 fParamLimits = new TMatrixD(nparam ,2);
309 for (Int_t iparam=0; iparam<nparam; iparam++){
310 (*fParamLimits)(iparam, 0) = 0;
311 (*fParamLimits)(iparam, 1) = 0;
312 }
313 }
314
315 // set all weights to 1 as default
dc697472 316 Bool_t weightFlag = false;
4b23a517 317 if (fWeightFunction == 0) {
318 fWeightFunction = new TFormula("constant", "1");
dc697472 319 } else {
320 weightFlag = true;
4b23a517 321 }
322
323 // migrad fit algorithm as default
bf7bd859 324 if (fFitAlgorithm == "") {
4b23a517 325 fFitAlgorithm = "migrad";
326 }
327
dc697472 328 // assign weights
329 if (fWeights == 0) {
330 fWeights = new TVectorD(npoints);
331 for (Int_t ipoint=0; ipoint<npoints; ipoint++) (*fWeights)(ipoint) = 1;
332 }
333
4b23a517 334 // set up the fitter
335 TVirtualFitter *minuit = TVirtualFitter::Fitter(0, nparam);
336 minuit->SetObjectFit(this);
dc697472 337 minuit->SetFCN(AliTMinuitToolkit::FitterFCN);
4b23a517 338
339 // initialize paramters (step size???)
340 for (Int_t iparam=0; iparam<nparam; iparam++){
341 minuit->SetParameter(iparam, Form("p[%d]",iparam), (*fParam)(iparam), (*fParam)(iparam)/10, (*fParamLimits)(iparam, 0), (*fParamLimits)(iparam, 1));
342 }
dc697472 343
344 //
4b23a517 345 Double_t argList[2];
346 argList[0] = fMaxCalls; //maximal number of calls
347 argList[1] = fPrecision; //tolerance normalized to 0.001
348 if (fMaxCalls == 500 && fPrecision == 1) minuit->ExecuteCommand(fFitAlgorithm, 0, 0);
349 if (fMaxCalls != 500 || fPrecision != 1) minuit->ExecuteCommand(fFitAlgorithm, argList, 2);
350 // two additional arguments can be specified ExecuteCommand("migrad", argList, 2) - use 0,0 for default
dc697472 351
352 // fill parameter vector
353 for (Int_t ivar=0; ivar<nparam; ivar++){
354 (*fParam)(ivar) = minuit->GetParameter(ivar);
355 fFormula->SetParameter(ivar, minuit->GetParameter(ivar));
356 }
4b23a517 357
dc697472 358 // if a weight function is specified -> enter 2nd run with weights
359 if (weightFlag == true && fUseRobust == true) {
360 // sort points for weighting
361 Double_t *sortList = new Double_t[npoints];
362 Int_t *indexList = new Int_t[npoints];
363 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
364 Double_t funx = fFormula->Eval((*fPoints)(ipoint, 0));
365 Double_t delta = TMath::Abs((*fPoints)[ipoint][nvar] - funx);
366 sortList[ipoint] = delta;
367 }
368 TMath::Sort(npoints, sortList, indexList, false);
369 for (Int_t ip=0; ip<npoints; ip++){
370 Double_t t = ip/(Double_t)npoints;
371 (*fWeights)(indexList[ip]) = fWeightFunction->Eval(t);
372 }
373
374 // set up the fitter
375 fUseRobust = false;
376 for (Int_t iparam=0; iparam<nparam; iparam++){
377 minuit->SetParameter(iparam, Form("p[%d]",iparam), (*fParam)(iparam), (*fParam)(iparam)/10, (*fParamLimits)(iparam, 0), (*fParamLimits)(iparam, 1));
378 }
379 // start fitting
380 if (fMaxCalls == 500 && fPrecision == 1) minuit->ExecuteCommand(fFitAlgorithm, 0, 0);
381 if (fMaxCalls != 500 || fPrecision != 1) minuit->ExecuteCommand(fFitAlgorithm, argList, 2);
382 fUseRobust = true;
383
4ce766eb 384 delete [] sortList;
385 delete [] indexList;
dc697472 386 }
387
4b23a517 388 // fill parameter vector
389 for (Int_t ivar=0; ivar<nparam; ivar++){
390 (*fParam)(ivar) = minuit->GetParameter(ivar);
dc697472 391 fFormula->SetParameter(ivar, minuit->GetParameter(ivar));
4b23a517 392 }
393
394 // fill covariance matrix
395 fCovar = new TMatrixD(nparam, nparam);
4b23a517 396 for(Int_t i=0; i < nparam; i++) {
397 for(Int_t j=0; j < nparam; j++) {
398 (*fCovar)(i,j) = minuit->GetCovarianceMatrixElement(i,j);
399 }
400 }
dc697472 401
402 if (weightFlag == false) fWeightFunction = 0;
4b23a517 403}
404
405
406
407void AliTMinuitToolkit::Test() {
408 //
409 // This test function shows the basic working principles of this class
410 // and illustrates how a robust fit can improve the results
411 //
305a4a95 412
413 // 1. provide some example histogram
dc697472 414 TH1F * hist = new TH1F("test", "with (red) and without (black) robust option", 20,0,4);
4b23a517 415 TRandom * rand = new TRandom();
416 for (Int_t i = 0; i < 10000; i++) {
417 hist->Fill(rand->Exp(1));
418 if (i < 1000) hist->Fill(3); //"outliers"
419 if (i < 1070) hist->Fill(3.5);
420 if (i < 670) hist->Fill(2);
421 if (i < 770) hist->Fill(1.5);//"outliers"
422 if (i < 740) hist->Fill(1);
423 }
424 TCanvas * canv = new TCanvas();
425 canv->cd(1);
426 hist->Draw();
dc697472 427
305a4a95 428 // 2. example fit without robust option
429 AliTMinuitToolkit * tool = new AliTMinuitToolkit();
b7d1d891 430 TFormula *aFormExp = new TFormula("formExp", "[0]*TMath::Exp(-[1]*x)");
431 tool->SetFitFunction(aFormExp);
305a4a95 432 TVectorD *vec1 = new TVectorD(2); // Set initial values
433 (*vec1)(0) = 1800;
434 (*vec1)(1) = 1;
435 tool->SetInitialParam(vec1);
436 tool->FitHistogram(hist);
437
4b23a517 438 // draw fit function
439 TF1 *func = new TF1("test", "[0]*TMath::Exp(-[1]*x)", 0, 6);
305a4a95 440 func->SetParameters((*tool->GetParameters())(0), (*tool->GetParameters())(1));
4b23a517 441 func->Draw("same");
dc697472 442
305a4a95 443 // 3 . robust fit
4b23a517 444 TVectorD *vec2 = new TVectorD(2);
445 (*vec2)(0) = 1800;
446 (*vec2)(1) = 1;
305a4a95 447 tool->SetInitialParam(vec2);
448 tool->EnableRobust(true, 10);
449 tool->SetWeightFunction("box", 0.75);
450 tool->FitHistogram(hist);
4b23a517 451 TF1 *func2 = new TF1("test2", "[0]*TMath::Exp(-[1]*x)", 0, 6);
305a4a95 452 func2->SetParameter(0, (*tool->GetParameters())(0));
453 func2->SetParameter(1, (*tool->GetParameters())(1));
4b23a517 454 func2->SetLineColor(kRed);
455 func2->Draw("same");
456
457}
458