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