removinf effc++ warnings (Marian)
[u/mrichter/AliRoot.git] / STAT / AliTMinuitToolkit.cxx
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
72 ClassImp(AliTMinuitToolkit)
73
74 AliTMinuitToolkit::AliTMinuitToolkit() : 
75    TNamed(),
76    fFormula(0),
77    fWeightFunction(0),
78    fFitAlgorithm(0),
79    fPoints(0),
80    fParam(0),
81    fParamLimits(0),
82    fCovar(0),
83    fChi2(0),
84    fMaxCalls(0),
85    fPrecision(0)
86 {
87  //
88  // standard constructor
89  //
90  fMaxCalls = 500;
91  fPrecision = 1;
92 }
93
94
95 AliTMinuitToolkit::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
113 AliTMinuitToolkit& AliTMinuitToolkit::operator=(const AliTMinuitToolkit&) {
114
115  return *this;
116 }
117
118
119
120 AliTMinuitToolkit::~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
133 void 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
151 void 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
176 void 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
190 void 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
206   Double_t *sortList = new Double_t[npoints];
207   Int_t  *indexList = new Int_t[npoints];
208   
209   TVectorD *Weight = new TVectorD(npoints);
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;
225    (*Weight)(ipoint) = fitter->GetWeightFunction()->Eval(t);
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;
237     fchisq+= delta*delta*(*Weight)(ipoint);
238
239   }
240   delete Weight;
241   delete sortList;
242   delete indexList;
243 }
244  
245
246 void 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
306 void 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