removinf effc++ warnings (Marian)
[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>
9#include <TVectorD.h>
10#include <TMatrixD.h>
11#include <TMath.h>
12#include <TString.h>
13#include <TROOT.h>
14#include <TCanvas.h>
15#include <TRandom.h>
16
17
18//--------------------------------------------------------------------------
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.
26// - An option for robust fitting of non-linear functions is implemented.
27//
28//
29// 1. Setting the formula:
30//
31// The formula is simply set via "void SetFitFunction(TFormula * formula)".
32//
33//
34// 2. Adding the data points
35//
36// - In order to fit a histogram, use "void FitHistogram(TH1F * his)" or
37// "void FitHistogram(TH2F * his)". The fitter is started automatically
38// - Alternatively, the direct specification of the points is possible via
39// "void SetPoints(TMatrixD * points)". Note, that the each point
40// corresponds to one row in the matrix. The fitter is then started with
41// the command "void Fit()".
42//
43//
44// 3. Accessing the fit results
45//
46// The N parameters of the formula are stored in a N-dimensional vector which
47// is returned by "TVectorD * GetParameters()". In a similar the covariance
48// matrix of the fit is returned via "TMatrixD * GetCovarianceMatrix()" which
49// is of the type N x N.
50//
51//
52// 4. Non-linear robust fitting:
53//
54// Even a few outliers can lead to wrong results of a least-squares fitting
55// procedure. In this case the use of robust(resistant) methods can be
56// helpful, but a stronger dependence on starting values or convergence to
57// local minima can occur.
58//
59// The robust option becomes active if a weighting function is specified.
60// All points are sorted according to their distance to the curve and
61// weighted. The weighting function must be defined on the interval [0,1].
62//
63// Some standard weighting functions are predefined in
64// "SetWeightFunction(Char_t * name, Float_t param1, Float_t param2 = 0)":
65// - "BOX" equals to 1 if x < param1 and to 0 if x > param1.
66// - "EXPONENTIAL" corresponds to "Math::Exp(-TMath::Log(param1)*x)"
67// - "ERRORFUNCTION" corresponds to "TMath::Erfc((x-param1)/param2)"
68//
69//-------------------------------------------------------------------------
70
71
72ClassImp(AliTMinuitToolkit)
73
74AliTMinuitToolkit::AliTMinuitToolkit() :
75 TNamed(),
76 fFormula(0),
77 fWeightFunction(0),
78 fFitAlgorithm(0),
79 fPoints(0),
80 fParam(0),
81 fParamLimits(0),
82 fCovar(0),
da7c9ae3 83 fChi2(0),
84 fMaxCalls(0),
85 fPrecision(0)
4b23a517 86{
87 //
88 // standard constructor
89 //
90 fMaxCalls = 500;
91 fPrecision = 1;
92}
93
94
da7c9ae3 95AliTMinuitToolkit::AliTMinuitToolkit(const AliTMinuitToolkit&) :
96 TNamed(),
97 fFormula(0),
98 fWeightFunction(0),
99 fFitAlgorithm(0),
100 fPoints(0),
101 fParam(0),
102 fParamLimits(0),
103 fCovar(0),
104 fChi2(0),
105 fMaxCalls(0),
106 fPrecision(0)
107{
108
109
110}
111
112
113AliTMinuitToolkit& AliTMinuitToolkit::operator=(const AliTMinuitToolkit&) {
114
115 return *this;
116}
117
118
4b23a517 119
120AliTMinuitToolkit::~AliTMinuitToolkit(){
121 //
122 // destructor
123 //
124 delete fPoints;
125 delete fWeightFunction;
126 delete fParamLimits;
127 delete fFormula;
128 delete fParam;
129 delete fCovar;
130 delete fChi2;
131}
132
133void AliTMinuitToolkit::FitHistogram(TH1F * his) {
134 //
135 // Fit a one dimensional histogram
136 //
137 fPoints = new TMatrixD(his->GetNbinsX(), 2);
138
139 for(Int_t ibin=0; ibin < his->GetNbinsX(); ibin++) {
140 Double_t x = his->GetXaxis()->GetBinCenter(ibin+1);
141 Double_t y = his->GetBinContent(ibin+1);
142
143 (*fPoints)(ibin, 0) = x;
144 (*fPoints)(ibin, 1) = y;
145 }
146
147 Fit();
148}
149
150
151void AliTMinuitToolkit::FitHistogram(TH2F * his) {
152 //
153 // Fit a two dimensional histogram
154 //
155 fPoints = new TMatrixD((Long64_t)his->GetEntries(), 2);
156 Long64_t entry = 0;
157
158 for(Int_t ibin=0; ibin < his->GetNbinsX(); ibin++) {
159 Double_t x = his->GetXaxis()->GetBinCenter(ibin);
160 for(Int_t jbin=0; jbin < his->GetNbinsY(); jbin++) {
161 Long64_t n = his->GetBin(ibin, jbin);
162 Double_t y = his->GetYaxis()->GetBinCenter(jbin);
163 for(Int_t ientries=0; ientries < his->GetBinContent(n); ientries++) {
164 (*fPoints)(entry,0) = x;
165 (*fPoints)(entry,1) = y;
166 entry++;
167 }
168
169 }
170 }
171
172 Fit();
173}
174
175
176void AliTMinuitToolkit::SetWeightFunction(Char_t * name, Float_t param1, Float_t param2) {
177 //
178 // Set the weight function which must be defined on the interval [0,1].
179 //
180 TString FuncType(name);
181 FuncType.ToUpper();
182
183 if (FuncType == "EXPONENTIAL") fWeightFunction = new TFormula("exp", Form("TMath::Exp(-TMath::Log(%f)*x)", param1));
184 if (FuncType == "BOX") fWeightFunction = new TFormula("box", Form("TMath::Erfc((x-%f)/0.0001)", param1));
185 if (FuncType == "ERRORFUNCTION") fWeightFunction = new TFormula("err", Form("TMath::Erfc((x-%f)/%f)", param1, param2));// !!!!!!!!!!!!!!!!!
186
187}
188
189
190void AliTMinuitToolkit::FitterFCN(int &npar, double *dummy, double &fchisq, double *gin, int iflag){
191 //
192 // internal function which gives the specified function to the TMinuit function
193 //
194
195 // suppress warnings for unused variables:
196 dummy = dummy;
197 iflag = iflag;
198 npar = npar;
199 //
200 AliTMinuitToolkit * fitter = (AliTMinuitToolkit*)TVirtualFitter::GetFitter()->GetObjectFit();
201 fchisq = 0;
202 Int_t nvar = fitter->GetPoints()->GetNcols()-1;
203 Int_t npoints = fitter->GetPoints()->GetNrows();
204
205 // sort points for weighting
d99c0447 206 Double_t *sortList = new Double_t[npoints];
da7c9ae3 207 Int_t *indexList = new Int_t[npoints];
4b23a517 208
da7c9ae3 209 TVectorD *Weight = new TVectorD(npoints);
4b23a517 210
211 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
212 Double_t x[100];
213 for (Int_t ivar=0; ivar<nvar; ivar++){
214 x[ivar] = (*fitter->GetPoints())(ipoint, ivar);
215 }
216 Float_t funx = fitter->GetFormula()->EvalPar(x,gin);
217 sortList[ipoint] = TMath::Abs((*fitter->GetPoints())(ipoint, nvar) - funx);
218 }
219
220 TMath::Sort(npoints, sortList, indexList, false);
221
222 Double_t t;
223 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
224 t = indexList[ipoint]/(Double_t)npoints;
da7c9ae3 225 (*Weight)(ipoint) = fitter->GetWeightFunction()->Eval(t);
4b23a517 226 }
227 //
228 // calculate chisquare
229 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
230 Double_t x[100];
231 for (Int_t ivar=0; ivar<nvar; ivar++){
232 x[ivar] = (*fitter->GetPoints())(ipoint, ivar);
233 }
234 Float_t funx = fitter->GetFormula()->EvalPar(x,gin);
235
236 Double_t delta = (*fitter->GetPoints())(ipoint, nvar) - funx;
da7c9ae3 237 fchisq+= delta*delta*(*Weight)(ipoint);
4b23a517 238
239 }
da7c9ae3 240 delete Weight;
241 delete sortList;
242 delete indexList;
4b23a517 243}
244
245
246void AliTMinuitToolkit::Fit() {
247 //
248 // internal function that calls the fitter
249 //
250 Int_t nparam = fParam->GetNrows();
251
252 // set all paramter limits to infinity as default
253 if (fParamLimits == 0) {
254 fParamLimits = new TMatrixD(nparam ,2);
255 for (Int_t iparam=0; iparam<nparam; iparam++){
256 (*fParamLimits)(iparam, 0) = 0;
257 (*fParamLimits)(iparam, 1) = 0;
258 }
259 }
260
261 // set all weights to 1 as default
262 if (fWeightFunction == 0) {
263 fWeightFunction = new TFormula("constant", "1");
264 }
265
266 // migrad fit algorithm as default
267 if (fFitAlgorithm == 0) {
268 fFitAlgorithm = "migrad";
269 }
270
271 // set up the fitter
272 TVirtualFitter *minuit = TVirtualFitter::Fitter(0, nparam);
273 minuit->SetObjectFit(this);
274 minuit->SetFCN((void*)(AliTMinuitToolkit::FitterFCN));
275
276 // initialize paramters (step size???)
277 for (Int_t iparam=0; iparam<nparam; iparam++){
278 minuit->SetParameter(iparam, Form("p[%d]",iparam), (*fParam)(iparam), (*fParam)(iparam)/10, (*fParamLimits)(iparam, 0), (*fParamLimits)(iparam, 1));
279 }
280
281 Double_t argList[2];
282 argList[0] = fMaxCalls; //maximal number of calls
283 argList[1] = fPrecision; //tolerance normalized to 0.001
284 if (fMaxCalls == 500 && fPrecision == 1) minuit->ExecuteCommand(fFitAlgorithm, 0, 0);
285 if (fMaxCalls != 500 || fPrecision != 1) minuit->ExecuteCommand(fFitAlgorithm, argList, 2);
286 // two additional arguments can be specified ExecuteCommand("migrad", argList, 2) - use 0,0 for default
287
288 // fill parameter vector
289 for (Int_t ivar=0; ivar<nparam; ivar++){
290 (*fParam)(ivar) = minuit->GetParameter(ivar);
291 }
292
293 // fill covariance matrix
294 fCovar = new TMatrixD(nparam, nparam);
295 //TVirtualFitter *fitCov = TVirtualFitter::GetFitter();
296 for(Int_t i=0; i < nparam; i++) {
297 for(Int_t j=0; j < nparam; j++) {
298 (*fCovar)(i,j) = minuit->GetCovarianceMatrixElement(i,j);
299 }
300 }
301
302}
303
304
305
306void AliTMinuitToolkit::Test() {
307 //
308 // This test function shows the basic working principles of this class
309 // and illustrates how a robust fit can improve the results
310 //
311 TFormula *FormExp = new TFormula("formExp", "[0]*TMath::Exp(-[1]*x)");
312 SetFitFunction(FormExp);
313 SetFitAlgorithm("simplex");
314 // Set initial values
315 TVectorD *vec1 = new TVectorD(2);
316 (*vec1)(0) = 1800;
317 (*vec1)(1) = 1;
318 SetInitialParam(vec1);
319 //provide some example histogram
320 TH1F * hist = new TH1F("bla", "with (red) and without (black) robust option", 20,0,4);
321 TRandom * rand = new TRandom();
322 for (Int_t i = 0; i < 10000; i++) {
323 hist->Fill(rand->Exp(1));
324 if (i < 1000) hist->Fill(3); //"outliers"
325 if (i < 1070) hist->Fill(3.5);
326 if (i < 670) hist->Fill(2);
327 if (i < 770) hist->Fill(1.5);//"outliers"
328 if (i < 740) hist->Fill(1);
329 }
330 TCanvas * canv = new TCanvas();
331 canv->cd(1);
332 hist->Draw();
333 // fit it with the exponential decay
334 FitHistogram(hist);
335 // draw fit function
336 TF1 *func = new TF1("test", "[0]*TMath::Exp(-[1]*x)", 0, 6);
337 func->SetParameter(0, (*GetParameters())(0));
338 func->SetParameter(1, (*GetParameters())(1));
339 func->Draw("same");
340 // robust fit
341 TVectorD *vec2 = new TVectorD(2);
342 (*vec2)(0) = 1800;
343 (*vec2)(1) = 1;
344 SetInitialParam(vec2);
345 SetWeightFunction("Box", 0.7);
346 FitHistogram(hist);
347 TF1 *func2 = new TF1("test2", "[0]*TMath::Exp(-[1]*x)", 0, 6);
348 func2->SetParameter(0, (*GetParameters())(0));
349 func2->SetParameter(1, (*GetParameters())(1));
350 func2->SetLineColor(kRed);
351 func2->Draw("same");
352
353}
354