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 | |
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), |
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 |
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 | |
4b23a517 |
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 |
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 | |
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 | |