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