Line data Source code
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 : **************************************************************************/
15 :
16 : /* $Id$ */
17 :
18 : #include <TCanvas.h>
19 : #include <TF1.h>
20 : #include <TFormula.h>
21 : #include <TH1F.h>
22 : #include <TH2F.h>
23 : #include <TMath.h>
24 : #include <TMatrix.h>
25 : #include <TRandom.h>
26 : #include <TString.h>
27 : #include <TVector.h>
28 : #include <TVectorD.h>
29 : #include <TVirtualFitter.h>
30 :
31 : #include "AliTMinuitToolkit.h"
32 :
33 : //--------------------------------------------------------------------------------------
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.
41 : // - An option for robust fitting with non-linear functions is implemented.
42 : //
43 : // A small example illustrating the usage of AliTMinuitToolkit is given in the function
44 : // "AliTMinuitToolkit::Test()".
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
59 : // the command "void Fit()". The weight of each point can be specified
60 : // with an n-dimensional vector using "void SetWeights(TVectorD * weights)".
61 : //
62 : //
63 : // 3. Accessing the fit results
64 : //
65 : // The N parameters of the formula are stored in a N-dimensional vector which
66 : // is returned by "TVectorD * GetParameters()". In a similar way the covariance
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 : //
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.
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 : //
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 : //--------------------------------------------------------------------------------------
122 :
123 :
124 128 : ClassImp(AliTMinuitToolkit)
125 :
126 : AliTMinuitToolkit::AliTMinuitToolkit() :
127 0 : TNamed(),
128 0 : fFormula(0),
129 0 : fWeightFunction(0),
130 0 : fFitAlgorithm(""),
131 0 : fPoints(0),
132 0 : fWeights(0),
133 0 : fParam(0),
134 0 : fParamLimits(0),
135 0 : fCovar(0),
136 0 : fChi2(0),
137 0 : fMaxCalls(0),
138 0 : fPrecision(0),
139 0 : fUseRobust(0),
140 0 : fExpectedSigma(0)
141 0 : {
142 : //
143 : // standard constructor
144 : //
145 0 : fMaxCalls = 500;
146 0 : fPrecision = 1;
147 0 : fUseRobust = false;
148 0 : fExpectedSigma = 0;
149 0 : }
150 :
151 :
152 : AliTMinuitToolkit::AliTMinuitToolkit(const AliTMinuitToolkit&) :
153 0 : TNamed(),
154 0 : fFormula(0),
155 0 : fWeightFunction(0),
156 0 : fFitAlgorithm(""),
157 0 : fPoints(0),
158 0 : fWeights(0),
159 0 : fParam(0),
160 0 : fParamLimits(0),
161 0 : fCovar(0),
162 0 : fChi2(0),
163 0 : fMaxCalls(0),
164 0 : fPrecision(0),
165 0 : fUseRobust(0),
166 0 : fExpectedSigma(0)
167 0 : {
168 :
169 :
170 0 : }
171 :
172 :
173 : AliTMinuitToolkit& AliTMinuitToolkit::operator=(const AliTMinuitToolkit&) {
174 :
175 0 : return *this;
176 : }
177 :
178 :
179 :
180 0 : AliTMinuitToolkit::~AliTMinuitToolkit(){
181 : //
182 : // destructor
183 : //
184 0 : delete fPoints;
185 0 : delete fWeights;
186 0 : delete fWeightFunction;
187 0 : delete fParamLimits;
188 0 : delete fFormula;
189 0 : delete fParam;
190 0 : delete fCovar;
191 0 : delete fChi2;
192 0 : }
193 :
194 : void AliTMinuitToolkit::FitHistogram(TH1F *const his) {
195 : //
196 : // Fit a one dimensional histogram
197 : //
198 0 : fPoints = new TMatrixD(his->GetNbinsX(), 2);
199 :
200 0 : for(Int_t ibin=0; ibin < his->GetNbinsX(); ibin++) {
201 0 : Double_t x = his->GetXaxis()->GetBinCenter(ibin+1);
202 0 : Double_t y = his->GetBinContent(ibin+1);
203 :
204 0 : (*fPoints)(ibin, 0) = x;
205 0 : (*fPoints)(ibin, 1) = y;
206 : }
207 :
208 0 : Fit();
209 0 : }
210 :
211 :
212 : void AliTMinuitToolkit::FitHistogram(TH2F *const his) {
213 : //
214 : // Fit a curve to a two dimensional histogram
215 : //
216 0 : fPoints = new TMatrixD((Long64_t)his->GetEntries(), 2);
217 : Long64_t entry = 0;
218 :
219 0 : for(Int_t ibin=0; ibin < his->GetNbinsX(); ibin++) {
220 0 : Double_t x = his->GetXaxis()->GetBinCenter(ibin);
221 0 : for(Int_t jbin=0; jbin < his->GetNbinsY(); jbin++) {
222 0 : Long64_t n = his->GetBin(ibin, jbin);
223 0 : Double_t y = his->GetYaxis()->GetBinCenter(jbin);
224 0 : for(Int_t ientries=0; ientries < his->GetBinContent(n); ientries++) {
225 0 : (*fPoints)(entry,0) = x;
226 0 : (*fPoints)(entry,1) = y;
227 0 : entry++;
228 : }
229 :
230 : }
231 : }
232 :
233 0 : Fit();
234 0 : }
235 :
236 :
237 : void AliTMinuitToolkit::SetWeightFunction(const Char_t *name, Float_t param1, Float_t param2) {
238 : //
239 : // Set the weight function which must be defined on the interval [0,1].
240 : //
241 0 : TString FuncType(name);
242 0 : FuncType.ToUpper();
243 :
244 0 : if (FuncType == "EXPONENTIAL") fWeightFunction = new TFormula("exp", Form("TMath::Exp(-TMath::Log(%f)*x)", param1));
245 0 : if (FuncType == "BOX") fWeightFunction = new TFormula("box", Form("TMath::Erfc((x-%f)/0.0001)", param1));
246 0 : if (FuncType == "ERRORFUNCTION") fWeightFunction = new TFormula("err", Form("TMath::Erfc((x-%f)/%f)", param1, param2));
247 :
248 0 : }
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 : //
257 0 : AliTMinuitToolkit * fitter = (AliTMinuitToolkit*)TVirtualFitter::GetFitter()->GetObjectFit();
258 0 : fchisq = 0;
259 0 : Int_t nvar = fitter->GetPoints()->GetNcols()-1;
260 0 : Int_t npoints = fitter->GetPoints()->GetNrows();
261 :
262 : // calculate mean deviation for normalization or use user-defined sigma
263 : Double_t dev = 0.;
264 0 : if (fitter->GetExpectedSigma() == 0 && fitter->GetStatus() == true) {
265 0 : for (Int_t ipoint=0; ipoint<npoints; ipoint++){
266 0 : Double_t x[100];
267 0 : for (Int_t ivar=0; ivar<nvar; ivar++){
268 0 : x[ivar] = (*fitter->GetPoints())(ipoint, ivar);
269 : }
270 0 : Float_t funx = fitter->GetFormula()->EvalPar(x,gin);
271 0 : Double_t delta = (*fitter->GetPoints())(ipoint, nvar) - funx;
272 0 : dev += TMath::Sqrt(TMath::Abs(delta));
273 0 : }
274 0 : dev = dev/npoints;
275 0 : } else {
276 0 : dev = fitter->GetExpectedSigma();
277 : }
278 : // calculate chisquare
279 0 : for (Int_t ipoint=0; ipoint<npoints; ipoint++){
280 0 : Double_t x[100];
281 0 : for (Int_t ivar=0; ivar<nvar; ivar++){
282 0 : x[ivar] = (*fitter->GetPoints())(ipoint, ivar);
283 : }
284 0 : Float_t funx = fitter->GetFormula()->EvalPar(x,gin);
285 0 : Double_t delta = TMath::Abs((*fitter->GetPoints())(ipoint, nvar) - funx);
286 0 : if (fitter->GetStatus() == true) {
287 0 : delta = delta/dev; // normalization
288 0 : if (delta <= 2.5) fchisq+= delta*delta; // new metric: Huber-k-estimator
289 0 : if (delta > 2.5) fchisq+= 2*(2.5)*delta - (2.5*2.5);
290 : } else {
291 0 : Double_t weight = (*fitter->GetWeights())(ipoint);
292 0 : fchisq+= delta*delta*weight; //old metric
293 : }
294 0 : }
295 0 : }
296 :
297 :
298 : void AliTMinuitToolkit::Fit() {
299 : //
300 : // internal function that calls the fitter
301 : //
302 0 : Int_t nparam = fParam->GetNrows();
303 0 : Int_t npoints = fPoints->GetNrows();
304 0 : Int_t nvar = fPoints->GetNcols()-1;
305 :
306 : // set all paramter limits to infinity as default
307 0 : if (fParamLimits == 0) {
308 0 : fParamLimits = new TMatrixD(nparam ,2);
309 0 : for (Int_t iparam=0; iparam<nparam; iparam++){
310 0 : (*fParamLimits)(iparam, 0) = 0;
311 0 : (*fParamLimits)(iparam, 1) = 0;
312 : }
313 0 : }
314 :
315 : // set all weights to 1 as default
316 : Bool_t weightFlag = false;
317 0 : if (fWeightFunction == 0) {
318 0 : fWeightFunction = new TFormula("constant", "1");
319 0 : } else {
320 : weightFlag = true;
321 : }
322 :
323 : // migrad fit algorithm as default
324 0 : if (fFitAlgorithm == "") {
325 0 : fFitAlgorithm = "migrad";
326 0 : }
327 :
328 : // assign weights
329 0 : if (fWeights == 0) {
330 0 : fWeights = new TVectorD(npoints);
331 0 : for (Int_t ipoint=0; ipoint<npoints; ipoint++) (*fWeights)(ipoint) = 1;
332 0 : }
333 :
334 : // set up the fitter
335 0 : TVirtualFitter *minuit = TVirtualFitter::Fitter(0, nparam);
336 0 : minuit->SetObjectFit(this);
337 0 : minuit->SetFCN(AliTMinuitToolkit::FitterFCN);
338 :
339 : // initialize paramters (step size???)
340 0 : for (Int_t iparam=0; iparam<nparam; iparam++){
341 0 : minuit->SetParameter(iparam, Form("p[%d]",iparam), (*fParam)(iparam), (*fParam)(iparam)/10, (*fParamLimits)(iparam, 0), (*fParamLimits)(iparam, 1));
342 : }
343 :
344 : //
345 0 : Double_t argList[2];
346 0 : argList[0] = fMaxCalls; //maximal number of calls
347 0 : argList[1] = fPrecision; //tolerance normalized to 0.001
348 0 : if (fMaxCalls == 500 && fPrecision == 1) minuit->ExecuteCommand(fFitAlgorithm, 0, 0);
349 0 : 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
351 :
352 : // fill parameter vector
353 0 : for (Int_t ivar=0; ivar<nparam; ivar++){
354 0 : (*fParam)(ivar) = minuit->GetParameter(ivar);
355 0 : fFormula->SetParameter(ivar, minuit->GetParameter(ivar));
356 : }
357 :
358 : // if a weight function is specified -> enter 2nd run with weights
359 0 : if (weightFlag == true && fUseRobust == true) {
360 : // sort points for weighting
361 0 : Double_t *sortList = new Double_t[npoints];
362 0 : Int_t *indexList = new Int_t[npoints];
363 0 : for (Int_t ipoint=0; ipoint<npoints; ipoint++){
364 0 : Double_t funx = fFormula->Eval((*fPoints)(ipoint, 0));
365 0 : Double_t delta = TMath::Abs((*fPoints)[ipoint][nvar] - funx);
366 0 : sortList[ipoint] = delta;
367 : }
368 0 : TMath::Sort(npoints, sortList, indexList, false);
369 0 : for (Int_t ip=0; ip<npoints; ip++){
370 0 : Double_t t = ip/(Double_t)npoints;
371 0 : (*fWeights)(indexList[ip]) = fWeightFunction->Eval(t);
372 : }
373 :
374 : // set up the fitter
375 0 : fUseRobust = false;
376 0 : for (Int_t iparam=0; iparam<nparam; iparam++){
377 0 : minuit->SetParameter(iparam, Form("p[%d]",iparam), (*fParam)(iparam), (*fParam)(iparam)/10, (*fParamLimits)(iparam, 0), (*fParamLimits)(iparam, 1));
378 : }
379 : // start fitting
380 0 : if (fMaxCalls == 500 && fPrecision == 1) minuit->ExecuteCommand(fFitAlgorithm, 0, 0);
381 0 : if (fMaxCalls != 500 || fPrecision != 1) minuit->ExecuteCommand(fFitAlgorithm, argList, 2);
382 0 : fUseRobust = true;
383 :
384 0 : delete [] sortList;
385 0 : delete [] indexList;
386 0 : }
387 :
388 : // fill parameter vector
389 0 : for (Int_t ivar=0; ivar<nparam; ivar++){
390 0 : (*fParam)(ivar) = minuit->GetParameter(ivar);
391 0 : fFormula->SetParameter(ivar, minuit->GetParameter(ivar));
392 : }
393 :
394 : // fill covariance matrix
395 0 : fCovar = new TMatrixD(nparam, nparam);
396 0 : for(Int_t i=0; i < nparam; i++) {
397 0 : for(Int_t j=0; j < nparam; j++) {
398 0 : (*fCovar)(i,j) = minuit->GetCovarianceMatrixElement(i,j);
399 : }
400 : }
401 :
402 0 : if (weightFlag == false) fWeightFunction = 0;
403 0 : }
404 :
405 :
406 :
407 : void 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 : //
412 :
413 : // 1. provide some example histogram
414 0 : TH1F * hist = new TH1F("test", "with (red) and without (black) robust option", 20,0,4);
415 0 : TRandom * rand = new TRandom();
416 0 : for (Int_t i = 0; i < 10000; i++) {
417 0 : hist->Fill(rand->Exp(1));
418 0 : if (i < 1000) hist->Fill(3); //"outliers"
419 0 : if (i < 1070) hist->Fill(3.5);
420 0 : if (i < 670) hist->Fill(2);
421 0 : if (i < 770) hist->Fill(1.5);//"outliers"
422 0 : if (i < 740) hist->Fill(1);
423 : }
424 0 : TCanvas * canv = new TCanvas();
425 0 : canv->cd(1);
426 0 : hist->Draw();
427 :
428 : // 2. example fit without robust option
429 0 : AliTMinuitToolkit * tool = new AliTMinuitToolkit();
430 0 : TFormula *aFormExp = new TFormula("formExp", "[0]*TMath::Exp(-[1]*x)");
431 0 : tool->SetFitFunction(aFormExp);
432 0 : TVectorD *vec1 = new TVectorD(2); // Set initial values
433 0 : (*vec1)(0) = 1800;
434 0 : (*vec1)(1) = 1;
435 0 : tool->SetInitialParam(vec1);
436 0 : tool->FitHistogram(hist);
437 :
438 : // draw fit function
439 0 : TF1 *func = new TF1("test", "[0]*TMath::Exp(-[1]*x)", 0, 6);
440 0 : func->SetParameters((*tool->GetParameters())(0), (*tool->GetParameters())(1));
441 0 : func->Draw("same");
442 :
443 : // 3 . robust fit
444 0 : TVectorD *vec2 = new TVectorD(2);
445 0 : (*vec2)(0) = 1800;
446 0 : (*vec2)(1) = 1;
447 0 : tool->SetInitialParam(vec2);
448 0 : tool->EnableRobust(true, 10);
449 0 : tool->SetWeightFunction("box", 0.75);
450 0 : tool->FitHistogram(hist);
451 0 : TF1 *func2 = new TF1("test2", "[0]*TMath::Exp(-[1]*x)", 0, 6);
452 0 : func2->SetParameter(0, (*tool->GetParameters())(0));
453 0 : func2->SetParameter(1, (*tool->GetParameters())(1));
454 0 : func2->SetLineColor(kRed);
455 0 : func2->Draw("same");
456 :
457 0 : }
458 :
|
