]>
Commit | Line | Data |
---|---|---|
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), | |
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 | |
d99c0447 | 180 | Double_t *sortList = new Double_t[npoints]; |
181 | Int_t *indexList = new Int_t[npoints]; | |
4b23a517 | 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; | |
d99c0447 | 215 | delete [] sortList; |
216 | delete [] indexList; | |
217 | ||
4b23a517 | 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 |