2 #include "AliTMinuitToolkit.h"
4 #include <TVirtualFitter.h>
18 //--------------------------------------------------------------------------
20 // The AliTMinuitToolkit serves as an easy to use interface for the TMinuit
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.
29 // 1. Setting the formula:
31 // The formula is simply set via "void SetFitFunction(TFormula * formula)".
34 // 2. Adding the data points
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()".
44 // 3. Accessing the fit results
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.
52 // 4. Non-linear robust fitting:
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.
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].
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)"
69 //-------------------------------------------------------------------------
72 ClassImp(AliTMinuitToolkit)
74 AliTMinuitToolkit::AliTMinuitToolkit() :
86 // standard constructor
94 AliTMinuitToolkit::~AliTMinuitToolkit(){
99 delete fWeightFunction;
107 void AliTMinuitToolkit::FitHistogram(TH1F * his) {
109 // Fit a one dimensional histogram
111 fPoints = new TMatrixD(his->GetNbinsX(), 2);
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);
117 (*fPoints)(ibin, 0) = x;
118 (*fPoints)(ibin, 1) = y;
125 void AliTMinuitToolkit::FitHistogram(TH2F * his) {
127 // Fit a two dimensional histogram
129 fPoints = new TMatrixD((Long64_t)his->GetEntries(), 2);
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;
150 void AliTMinuitToolkit::SetWeightFunction(Char_t * name, Float_t param1, Float_t param2) {
152 // Set the weight function which must be defined on the interval [0,1].
154 TString FuncType(name);
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));// !!!!!!!!!!!!!!!!!
164 void AliTMinuitToolkit::FitterFCN(int &npar, double *dummy, double &fchisq, double *gin, int iflag){
166 // internal function which gives the specified function to the TMinuit function
169 // suppress warnings for unused variables:
174 AliTMinuitToolkit * fitter = (AliTMinuitToolkit*)TVirtualFitter::GetFitter()->GetObjectFit();
176 Int_t nvar = fitter->GetPoints()->GetNcols()-1;
177 Int_t npoints = fitter->GetPoints()->GetNrows();
179 // sort points for weighting
180 Double_t *sortList = new Double_t[npoints];
181 Int_t *indexList = new Int_t[npoints];
183 TVectorD *fWeight = new TVectorD(npoints);
185 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
187 for (Int_t ivar=0; ivar<nvar; ivar++){
188 x[ivar] = (*fitter->GetPoints())(ipoint, ivar);
190 Float_t funx = fitter->GetFormula()->EvalPar(x,gin);
191 sortList[ipoint] = TMath::Abs((*fitter->GetPoints())(ipoint, nvar) - funx);
194 TMath::Sort(npoints, sortList, indexList, false);
197 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
198 t = indexList[ipoint]/(Double_t)npoints;
199 (*fWeight)(ipoint) = fitter->GetWeightFunction()->Eval(t);
202 // calculate chisquare
203 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
205 for (Int_t ivar=0; ivar<nvar; ivar++){
206 x[ivar] = (*fitter->GetPoints())(ipoint, ivar);
208 Float_t funx = fitter->GetFormula()->EvalPar(x,gin);
210 Double_t delta = (*fitter->GetPoints())(ipoint, nvar) - funx;
211 fchisq+= delta*delta*(*fWeight)(ipoint);
221 void AliTMinuitToolkit::Fit() {
223 // internal function that calls the fitter
225 Int_t nparam = fParam->GetNrows();
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;
236 // set all weights to 1 as default
237 if (fWeightFunction == 0) {
238 fWeightFunction = new TFormula("constant", "1");
241 // migrad fit algorithm as default
242 if (fFitAlgorithm == 0) {
243 fFitAlgorithm = "migrad";
247 TVirtualFitter *minuit = TVirtualFitter::Fitter(0, nparam);
248 minuit->SetObjectFit(this);
249 minuit->SetFCN((void*)(AliTMinuitToolkit::FitterFCN));
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));
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
263 // fill parameter vector
264 for (Int_t ivar=0; ivar<nparam; ivar++){
265 (*fParam)(ivar) = minuit->GetParameter(ivar);
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);
281 void AliTMinuitToolkit::Test() {
283 // This test function shows the basic working principles of this class
284 // and illustrates how a robust fit can improve the results
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);
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);
305 TCanvas * canv = new TCanvas();
308 // fit it with the exponential decay
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));
316 TVectorD *vec2 = new TVectorD(2);
319 SetInitialParam(vec2);
320 SetWeightFunction("Box", 0.7);
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);