alocate memory dynamically not on stack
[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 {
85  //
86  // standard constructor
87  //
88  fMaxCalls = 500;
89  fPrecision = 1;
90 }
91
92
93
94 AliTMinuitToolkit::~AliTMinuitToolkit(){
95   //
96   // destructor
97   //
98   delete fPoints;
99   delete fWeightFunction;
100   delete fParamLimits;
101   delete fFormula;
102   delete fParam;
103   delete fCovar;
104   delete fChi2;
105 }
106
107 void AliTMinuitToolkit::FitHistogram(TH1F * his) {
108  //
109  // Fit a one dimensional histogram
110  //
111  fPoints = new TMatrixD(his->GetNbinsX(), 2);
112  
113  for(Int_t ibin=0; ibin < his->GetNbinsX(); ibin++) {
114   Double_t x = his->GetXaxis()->GetBinCenter(ibin+1);
115   Double_t y = his->GetBinContent(ibin+1);
116   
117   (*fPoints)(ibin, 0) = x;
118   (*fPoints)(ibin, 1) = y;
119  }
120  
121  Fit();
122 }
123
124
125 void AliTMinuitToolkit::FitHistogram(TH2F * his) {
126  //
127  // Fit a two dimensional histogram
128  //
129  fPoints = new TMatrixD((Long64_t)his->GetEntries(), 2);
130  Long64_t entry = 0;
131  
132  for(Int_t ibin=0; ibin < his->GetNbinsX(); ibin++) {
133   Double_t x = his->GetXaxis()->GetBinCenter(ibin);
134   for(Int_t jbin=0; jbin < his->GetNbinsY(); jbin++) {   
135    Long64_t n = his->GetBin(ibin, jbin);
136    Double_t y = his->GetYaxis()->GetBinCenter(jbin);
137    for(Int_t ientries=0; ientries < his->GetBinContent(n); ientries++) {
138     (*fPoints)(entry,0) = x;
139     (*fPoints)(entry,1) = y;
140     entry++;
141    }
142    
143   }
144  }
145
146  Fit();
147 }
148
149
150 void AliTMinuitToolkit::SetWeightFunction(Char_t * name, Float_t param1, Float_t param2) {
151  //
152  // Set the weight function which must be defined on the interval [0,1].
153  //
154  TString FuncType(name);
155  FuncType.ToUpper();
156  
157  if (FuncType == "EXPONENTIAL") fWeightFunction = new TFormula("exp", Form("TMath::Exp(-TMath::Log(%f)*x)", param1));
158  if (FuncType == "BOX") fWeightFunction = new TFormula("box", Form("TMath::Erfc((x-%f)/0.0001)", param1));
159  if (FuncType == "ERRORFUNCTION") fWeightFunction = new TFormula("err", Form("TMath::Erfc((x-%f)/%f)", param1, param2));// !!!!!!!!!!!!!!!!!
160  
161 }
162
163
164 void AliTMinuitToolkit::FitterFCN(int &npar, double *dummy, double &fchisq, double *gin, int iflag){
165   //
166   // internal function which gives the specified function to the TMinuit function
167   //  
168
169   // suppress warnings for unused variables:
170   dummy = dummy;
171   iflag = iflag;
172   npar = npar;
173   //
174   AliTMinuitToolkit * fitter = (AliTMinuitToolkit*)TVirtualFitter::GetFitter()->GetObjectFit();
175   fchisq = 0;
176   Int_t nvar       = fitter->GetPoints()->GetNcols()-1;
177   Int_t npoints    = fitter->GetPoints()->GetNrows();
178   
179   // sort points for weighting
180   Double_t *sortList = new Double_t[npoints];
181   Int_t *indexList   = new Int_t[npoints];
182   
183   TVectorD *fWeight = new TVectorD(npoints);
184   
185   for (Int_t ipoint=0; ipoint<npoints; ipoint++){
186     Double_t x[100];
187     for (Int_t ivar=0; ivar<nvar; ivar++){
188       x[ivar] = (*fitter->GetPoints())(ipoint, ivar);      
189     }
190     Float_t funx   = fitter->GetFormula()->EvalPar(x,gin);
191     sortList[ipoint] = TMath::Abs((*fitter->GetPoints())(ipoint, nvar) - funx);
192   }
193   
194   TMath::Sort(npoints, sortList, indexList, false);
195
196   Double_t t;
197   for (Int_t ipoint=0; ipoint<npoints; ipoint++){
198    t = indexList[ipoint]/(Double_t)npoints;
199    (*fWeight)(ipoint) = fitter->GetWeightFunction()->Eval(t);
200   }
201   //
202   // calculate chisquare
203   for (Int_t ipoint=0; ipoint<npoints; ipoint++){
204     Double_t x[100];
205     for (Int_t ivar=0; ivar<nvar; ivar++){
206       x[ivar] = (*fitter->GetPoints())(ipoint, ivar);      
207     }
208     Float_t funx   = fitter->GetFormula()->EvalPar(x,gin);
209     
210     Double_t delta = (*fitter->GetPoints())(ipoint, nvar) - funx;
211     fchisq+= delta*delta*(*fWeight)(ipoint);
212
213   }
214   delete fWeight;
215   delete [] sortList;
216   delete [] indexList;
217
218 }
219  
220
221 void AliTMinuitToolkit::Fit() {
222   //
223   // internal function that calls the fitter
224   //
225   Int_t nparam = fParam->GetNrows();
226   
227   // set all paramter limits to infinity as default
228   if (fParamLimits == 0) {
229    fParamLimits = new TMatrixD(nparam ,2);
230    for (Int_t iparam=0; iparam<nparam; iparam++){
231     (*fParamLimits)(iparam, 0) = 0;
232     (*fParamLimits)(iparam, 1) = 0;
233    }
234   }
235   
236   // set all weights to 1 as default
237   if (fWeightFunction == 0) {
238    fWeightFunction = new TFormula("constant", "1");
239   }
240   
241   // migrad fit algorithm as default
242   if (fFitAlgorithm == 0) {
243    fFitAlgorithm = "migrad";
244   }
245   
246   // set up the fitter
247   TVirtualFitter *minuit = TVirtualFitter::Fitter(0, nparam);
248   minuit->SetObjectFit(this); 
249   minuit->SetFCN((void*)(AliTMinuitToolkit::FitterFCN));
250   
251   // initialize paramters (step size???)
252   for (Int_t iparam=0; iparam<nparam; iparam++){
253    minuit->SetParameter(iparam, Form("p[%d]",iparam), (*fParam)(iparam), (*fParam)(iparam)/10, (*fParamLimits)(iparam, 0), (*fParamLimits)(iparam, 1));
254   }
255   
256   Double_t argList[2];
257   argList[0] = fMaxCalls; //maximal number of calls 
258   argList[1] = fPrecision; //tolerance normalized to 0.001 
259   if (fMaxCalls == 500 && fPrecision == 1) minuit->ExecuteCommand(fFitAlgorithm, 0, 0); 
260   if (fMaxCalls != 500 || fPrecision != 1) minuit->ExecuteCommand(fFitAlgorithm, argList, 2);
261   // two additional arguments can be specified ExecuteCommand("migrad", argList, 2) - use 0,0 for default
262
263   // fill parameter vector
264   for (Int_t ivar=0; ivar<nparam; ivar++){
265    (*fParam)(ivar) = minuit->GetParameter(ivar);
266   }
267   
268   // fill covariance matrix
269   fCovar = new TMatrixD(nparam, nparam);
270   //TVirtualFitter *fitCov = TVirtualFitter::GetFitter();
271   for(Int_t i=0; i < nparam; i++) {
272    for(Int_t j=0; j < nparam; j++) {
273     (*fCovar)(i,j) = minuit->GetCovarianceMatrixElement(i,j);
274    }
275   }
276  
277 }
278
279
280
281 void AliTMinuitToolkit::Test() {
282  //
283  // This test function shows the basic working principles of this class 
284  // and illustrates how a robust fit can improve the results
285  //
286  TFormula *FormExp = new TFormula("formExp", "[0]*TMath::Exp(-[1]*x)");
287  SetFitFunction(FormExp);
288  SetFitAlgorithm("simplex");
289  // Set initial values
290  TVectorD *vec1 = new TVectorD(2);
291  (*vec1)(0) = 1800;
292  (*vec1)(1) = 1;
293  SetInitialParam(vec1);
294  //provide some example histogram
295  TH1F * hist = new TH1F("bla", "with (red) and without (black) robust option", 20,0,4);
296  TRandom * rand = new TRandom();
297  for (Int_t i = 0; i < 10000; i++) {
298   hist->Fill(rand->Exp(1));
299   if (i < 1000) hist->Fill(3); //"outliers"
300   if (i < 1070) hist->Fill(3.5);
301   if (i < 670) hist->Fill(2);
302   if (i < 770) hist->Fill(1.5);//"outliers"
303   if (i < 740) hist->Fill(1);
304  }
305  TCanvas * canv = new TCanvas();
306  canv->cd(1);
307  hist->Draw();
308  // fit it with the exponential decay
309  FitHistogram(hist);
310  // draw fit function
311  TF1 *func = new TF1("test", "[0]*TMath::Exp(-[1]*x)", 0, 6);
312  func->SetParameter(0, (*GetParameters())(0));
313  func->SetParameter(1, (*GetParameters())(1));
314  func->Draw("same");
315  // robust fit
316  TVectorD *vec2 = new TVectorD(2);
317  (*vec2)(0) = 1800;
318  (*vec2)(1) = 1;
319  SetInitialParam(vec2);
320  SetWeightFunction("Box", 0.7);
321  FitHistogram(hist);
322  TF1 *func2 = new TF1("test2", "[0]*TMath::Exp(-[1]*x)", 0, 6);
323  func2->SetParameter(0, (*GetParameters())(0));
324  func2->SetParameter(1, (*GetParameters())(1));
325  func2->SetLineColor(kRed);
326  func2->Draw("same");
327  
328 }
329