Changing once more (hopefully we get it correct this time...) the logic to trig the...
[u/mrichter/AliRoot.git] / STAT / AliTMinuitToolkit.cxx
CommitLineData
4b23a517 1
2#include "AliTMinuitToolkit.h"
3#include <TNamed.h>
4#include <TVirtualFitter.h>
5#include <TH1F.h>
6#include <TH2F.h>
7#include <TF1.h>
8#include <TFormula.h>
dc697472 9#include <TVector.h>
10#include <TMatrix.h>
4b23a517 11#include <TMath.h>
12#include <TString.h>
13#include <TROOT.h>
14#include <TCanvas.h>
15#include <TRandom.h>
16
17
dc697472 18//--------------------------------------------------------------------------------------
4b23a517 19//
20// The AliTMinuitToolkit serves as an easy to use interface for the TMinuit
21// package:
22//
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.
dc697472 26// - An option for robust fitting with non-linear functions is implemented.
305a4a95 27//
28// A small example illustrating the usage of AliTMinuitToolkit is given in the function
29// "AliTMinuitToolkit::Test()".
4b23a517 30//
31//
32// 1. Setting the formula:
33//
34// The formula is simply set via "void SetFitFunction(TFormula * formula)".
35//
36//
37// 2. Adding the data points
38//
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
dc697472 44// the command "void Fit()". The weight of each point can be specified
45// with an n-dimensional vector using "void SetWeights(TVectorD * weights)".
4b23a517 46//
47//
48// 3. Accessing the fit results
49//
50// The N parameters of the formula are stored in a N-dimensional vector which
dc697472 51// is returned by "TVectorD * GetParameters()". In a similar way the covariance
4b23a517 52// matrix of the fit is returned via "TMatrixD * GetCovarianceMatrix()" which
53// is of the type N x N.
54//
55//
56// 4. Non-linear robust fitting:
57//
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.
62//
dc697472 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)".
67//
68// Performing the fit without knowledge of sigma is also possible if only
69// "EnableRobust(true)" is activated, but this is NOT RECOMMENDED.
70//
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:
74//
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
77//
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
83// largest distance.
4b23a517 84//
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)"
90//
dc697472 91//
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/
96//
97//
98// 5. examples:
99//
100// A small example illustrating the working principles of AliTMinuitToolkit is given
101// in the function "AliTMinuitToolkit::Test()".
102//
103//
104//
105// Comments and questions are always welcome: A.Kalweit@gsi.de
106//--------------------------------------------------------------------------------------
4b23a517 107
108
109ClassImp(AliTMinuitToolkit)
110
111AliTMinuitToolkit::AliTMinuitToolkit() :
112 TNamed(),
113 fFormula(0),
114 fWeightFunction(0),
115 fFitAlgorithm(0),
116 fPoints(0),
dc697472 117 fWeights(0),
4b23a517 118 fParam(0),
119 fParamLimits(0),
120 fCovar(0),
da7c9ae3 121 fChi2(0),
122 fMaxCalls(0),
dc697472 123 fPrecision(0),
124 fUseRobust(0),
125 fExpectedSigma(0)
4b23a517 126{
127 //
128 // standard constructor
129 //
130 fMaxCalls = 500;
131 fPrecision = 1;
dc697472 132 fUseRobust = false;
133 fExpectedSigma = 0;
4b23a517 134}
135
136
da7c9ae3 137AliTMinuitToolkit::AliTMinuitToolkit(const AliTMinuitToolkit&) :
138 TNamed(),
139 fFormula(0),
140 fWeightFunction(0),
141 fFitAlgorithm(0),
142 fPoints(0),
dc697472 143 fWeights(0),
da7c9ae3 144 fParam(0),
145 fParamLimits(0),
146 fCovar(0),
147 fChi2(0),
148 fMaxCalls(0),
dc697472 149 fPrecision(0),
150 fUseRobust(0),
151 fExpectedSigma(0)
da7c9ae3 152{
153
154
155}
156
157
158AliTMinuitToolkit& AliTMinuitToolkit::operator=(const AliTMinuitToolkit&) {
159
160 return *this;
161}
162
163
4b23a517 164
165AliTMinuitToolkit::~AliTMinuitToolkit(){
166 //
167 // destructor
168 //
169 delete fPoints;
dc697472 170 delete fWeights;
4b23a517 171 delete fWeightFunction;
172 delete fParamLimits;
173 delete fFormula;
174 delete fParam;
175 delete fCovar;
176 delete fChi2;
177}
178
179void AliTMinuitToolkit::FitHistogram(TH1F * his) {
180 //
181 // Fit a one dimensional histogram
182 //
183 fPoints = new TMatrixD(his->GetNbinsX(), 2);
184
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);
188
189 (*fPoints)(ibin, 0) = x;
190 (*fPoints)(ibin, 1) = y;
191 }
192
193 Fit();
194}
195
196
197void AliTMinuitToolkit::FitHistogram(TH2F * his) {
198 //
dc697472 199 // Fit a curve to a two dimensional histogram
4b23a517 200 //
201 fPoints = new TMatrixD((Long64_t)his->GetEntries(), 2);
202 Long64_t entry = 0;
203
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;
212 entry++;
213 }
214
215 }
216 }
217
218 Fit();
219}
220
221
222void AliTMinuitToolkit::SetWeightFunction(Char_t * name, Float_t param1, Float_t param2) {
223 //
224 // Set the weight function which must be defined on the interval [0,1].
225 //
226 TString FuncType(name);
227 FuncType.ToUpper();
228
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));
dc697472 231 if (FuncType == "ERRORFUNCTION") fWeightFunction = new TFormula("err", Form("TMath::Erfc((x-%f)/%f)", param1, param2));
4b23a517 232
233}
234
235
236void AliTMinuitToolkit::FitterFCN(int &npar, double *dummy, double &fchisq, double *gin, int iflag){
237 //
238 // internal function which gives the specified function to the TMinuit function
239 //
240
241 // suppress warnings for unused variables:
242 dummy = dummy;
243 iflag = iflag;
244 npar = npar;
245 //
246 AliTMinuitToolkit * fitter = (AliTMinuitToolkit*)TVirtualFitter::GetFitter()->GetObjectFit();
247 fchisq = 0;
248 Int_t nvar = fitter->GetPoints()->GetNcols()-1;
249 Int_t npoints = fitter->GetPoints()->GetNrows();
250
dc697472 251 // calculate mean deviation for normalization or use user-defined sigma
252 Double_t dev = 0.;
253 if (fitter->GetExpectedSigma() == 0 && fitter->GetStatus() == true) {
254 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
4b23a517 255 Double_t x[100];
dc697472 256 for (Int_t ivar=0; ivar<nvar; ivar++){
4b23a517 257 x[ivar] = (*fitter->GetPoints())(ipoint, ivar);
dc697472 258 }
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));
262 }
263 dev = dev/npoints;
264 } else {
265 dev = fitter->GetExpectedSigma();
4b23a517 266 }
dc697472 267 // calculate chisquare
4b23a517 268 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
269 Double_t x[100];
270 for (Int_t ivar=0; ivar<nvar; ivar++){
271 x[ivar] = (*fitter->GetPoints())(ipoint, ivar);
272 }
dc697472 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);
279 } else {
280 Double_t weight = (*fitter->GetWeights())(ipoint);
281 fchisq+= delta*delta*weight; //old metric
282 }
4b23a517 283 }
dc697472 284 }
4b23a517 285
286
287void AliTMinuitToolkit::Fit() {
288 //
289 // internal function that calls the fitter
290 //
dc697472 291 Int_t nparam = fParam->GetNrows();
292 Int_t npoints = fPoints->GetNrows();
293 Int_t nvar = fPoints->GetNcols()-1;
4b23a517 294
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;
301 }
302 }
303
304 // set all weights to 1 as default
dc697472 305 Bool_t weightFlag = false;
4b23a517 306 if (fWeightFunction == 0) {
307 fWeightFunction = new TFormula("constant", "1");
dc697472 308 } else {
309 weightFlag = true;
4b23a517 310 }
311
312 // migrad fit algorithm as default
313 if (fFitAlgorithm == 0) {
314 fFitAlgorithm = "migrad";
315 }
316
dc697472 317 // assign weights
318 if (fWeights == 0) {
319 fWeights = new TVectorD(npoints);
320 for (Int_t ipoint=0; ipoint<npoints; ipoint++) (*fWeights)(ipoint) = 1;
321 }
322
4b23a517 323 // set up the fitter
324 TVirtualFitter *minuit = TVirtualFitter::Fitter(0, nparam);
325 minuit->SetObjectFit(this);
dc697472 326 minuit->SetFCN(AliTMinuitToolkit::FitterFCN);
4b23a517 327
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));
331 }
dc697472 332
333 //
4b23a517 334 Double_t argList[2];
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
dc697472 340
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));
345 }
4b23a517 346
dc697472 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;
356 }
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);
361 }
362
363 // set up the fitter
364 fUseRobust = false;
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));
367 }
368 // start fitting
369 if (fMaxCalls == 500 && fPrecision == 1) minuit->ExecuteCommand(fFitAlgorithm, 0, 0);
370 if (fMaxCalls != 500 || fPrecision != 1) minuit->ExecuteCommand(fFitAlgorithm, argList, 2);
371 fUseRobust = true;
372
373 delete sortList;
374 delete indexList;
375 }
376
4b23a517 377 // fill parameter vector
378 for (Int_t ivar=0; ivar<nparam; ivar++){
379 (*fParam)(ivar) = minuit->GetParameter(ivar);
dc697472 380 fFormula->SetParameter(ivar, minuit->GetParameter(ivar));
4b23a517 381 }
382
383 // fill covariance matrix
384 fCovar = new TMatrixD(nparam, nparam);
4b23a517 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);
388 }
389 }
dc697472 390
391 if (weightFlag == false) fWeightFunction = 0;
4b23a517 392}
393
394
395
396void AliTMinuitToolkit::Test() {
397 //
398 // This test function shows the basic working principles of this class
399 // and illustrates how a robust fit can improve the results
400 //
305a4a95 401
402 // 1. provide some example histogram
dc697472 403 TH1F * hist = new TH1F("test", "with (red) and without (black) robust option", 20,0,4);
4b23a517 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);
412 }
413 TCanvas * canv = new TCanvas();
414 canv->cd(1);
415 hist->Draw();
dc697472 416
305a4a95 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
422 (*vec1)(0) = 1800;
423 (*vec1)(1) = 1;
424 tool->SetInitialParam(vec1);
425 tool->FitHistogram(hist);
426
4b23a517 427 // draw fit function
428 TF1 *func = new TF1("test", "[0]*TMath::Exp(-[1]*x)", 0, 6);
305a4a95 429 func->SetParameters((*tool->GetParameters())(0), (*tool->GetParameters())(1));
4b23a517 430 func->Draw("same");
dc697472 431
305a4a95 432 // 3 . robust fit
4b23a517 433 TVectorD *vec2 = new TVectorD(2);
434 (*vec2)(0) = 1800;
435 (*vec2)(1) = 1;
305a4a95 436 tool->SetInitialParam(vec2);
437 tool->EnableRobust(true, 10);
438 tool->SetWeightFunction("box", 0.75);
439 tool->FitHistogram(hist);
4b23a517 440 TF1 *func2 = new TF1("test2", "[0]*TMath::Exp(-[1]*x)", 0, 6);
305a4a95 441 func2->SetParameter(0, (*tool->GetParameters())(0));
442 func2->SetParameter(1, (*tool->GetParameters())(1));
4b23a517 443 func2->SetLineColor(kRed);
444 func2->Draw("same");
445
446}
447