]>
Commit | Line | Data |
---|---|---|
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 | ||
124 | ClassImp(AliTMinuitToolkit) | |
125 | ||
126 | AliTMinuitToolkit::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 | 152 | AliTMinuitToolkit::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 | ||
173 | AliTMinuitToolkit& AliTMinuitToolkit::operator=(const AliTMinuitToolkit&) { | |
174 | ||
175 | return *this; | |
176 | } | |
177 | ||
178 | ||
4b23a517 | 179 | |
180 | AliTMinuitToolkit::~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 | 194 | void 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 | 212 | void 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 | 237 | void 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 | ||
251 | void 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 | ||
302 | void 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 | ||
411 | void 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 |