Added AliTMinuitToolkit (A.Kalweit)
authormarian <marian@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 17 Sep 2007 12:27:51 +0000 (12:27 +0000)
committermarian <marian@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 17 Sep 2007 12:27:51 +0000 (12:27 +0000)
// The AliTMinuitToolkit serves as an easy to use interface for the TMinuit
// package:
//
// - It allows to fit a curve to one and two dimensional histograms
//   (TH2F::Fit() only allows to fit a hyperplane).
// - Or n points can be specified directly via a n x 2 matrix.
// - An option for robust fitting of non-linear functions is implemented.

--------------------------------------------------------------------

STAT/AliTMinuitToolkit.cxx [new file with mode: 0644]
STAT/AliTMinuitToolkit.h [new file with mode: 0644]
STAT/Makefile
STAT/kDStatLinkDef.h
STAT/libKDStat.pkg

diff --git a/STAT/AliTMinuitToolkit.cxx b/STAT/AliTMinuitToolkit.cxx
new file mode 100644 (file)
index 0000000..23c6f60
--- /dev/null
@@ -0,0 +1,326 @@
+
+#include "AliTMinuitToolkit.h"
+#include <TNamed.h>
+#include <TVirtualFitter.h>
+#include <TH1F.h>
+#include <TH2F.h>
+#include <TF1.h>
+#include <TFormula.h>
+#include <TVectorD.h>
+#include <TMatrixD.h>
+#include <TMath.h>
+#include <TString.h>
+#include <TROOT.h>
+#include <TCanvas.h>
+#include <TRandom.h>
+
+
+//--------------------------------------------------------------------------
+//
+// The AliTMinuitToolkit serves as an easy to use interface for the TMinuit 
+// package: 
+// 
+// - It allows to fit a curve to one and two dimensional histograms 
+//   (TH2F::Fit() only allows to fit a hyperplane).
+// - Or n points can be specified directly via a n x 2 matrix.
+// - An option for robust fitting of non-linear functions is implemented.
+// 
+// 
+// 1. Setting the formula:
+//
+//  The formula is simply set via "void SetFitFunction(TFormula * formula)".
+//
+//
+// 2. Adding the data points
+//
+//  - In order to fit a histogram, use "void FitHistogram(TH1F * his)" or 
+//    "void FitHistogram(TH2F * his)". The fitter is started automatically
+//  - Alternatively, the direct specification of the points is possible via
+//    "void SetPoints(TMatrixD * points)". Note, that the each point 
+//    corresponds to one row in the matrix. The fitter is then started with 
+//    the command "void Fit()".
+//
+//
+// 3. Accessing the fit results
+//
+//  The N parameters of the formula are stored in a N-dimensional vector which
+//  is returned by "TVectorD * GetParameters()". In a similar the covariance 
+//  matrix of the fit is returned via "TMatrixD * GetCovarianceMatrix()" which
+//  is of the type N x N.
+//
+//
+// 4. Non-linear robust fitting:
+//
+//  Even a few outliers can lead to wrong results of a least-squares fitting 
+//  procedure. In this case the use of robust(resistant) methods can be 
+//  helpful, but a stronger dependence on starting values or convergence to
+//  local minima can occur.
+//
+//  The robust option becomes active if a weighting function is specified. 
+//  All points are sorted according to their distance to the curve and
+//  weighted. The weighting function must be defined on the interval [0,1].
+//
+//  Some standard weighting functions are predefined in 
+//  "SetWeightFunction(Char_t * name, Float_t param1, Float_t param2 = 0)":
+//   - "BOX" equals to 1 if x < param1 and to 0 if x > param1.
+//   - "EXPONENTIAL" corresponds to "Math::Exp(-TMath::Log(param1)*x)"
+//   - "ERRORFUNCTION" corresponds to "TMath::Erfc((x-param1)/param2)"
+//
+//-------------------------------------------------------------------------
+
+
+ClassImp(AliTMinuitToolkit)
+
+AliTMinuitToolkit::AliTMinuitToolkit() : 
+   TNamed(),
+   fFormula(0),
+   fWeightFunction(0),
+   fFitAlgorithm(0),
+   fPoints(0),
+   fParam(0),
+   fParamLimits(0),
+   fCovar(0),
+   fChi2(0)
+{
+ //
+ // standard constructor
+ //
+ fMaxCalls = 500;
+ fPrecision = 1;
+}
+
+
+
+AliTMinuitToolkit::~AliTMinuitToolkit(){
+  //
+  // destructor
+  //
+  delete fPoints;
+  delete fWeightFunction;
+  delete fParamLimits;
+  delete fFormula;
+  delete fParam;
+  delete fCovar;
+  delete fChi2;
+}
+
+void AliTMinuitToolkit::FitHistogram(TH1F * his) {
+ //
+ // Fit a one dimensional histogram
+ //
+ fPoints = new TMatrixD(his->GetNbinsX(), 2);
+ for(Int_t ibin=0; ibin < his->GetNbinsX(); ibin++) {
+  Double_t x = his->GetXaxis()->GetBinCenter(ibin+1);
+  Double_t y = his->GetBinContent(ibin+1);
+  
+  (*fPoints)(ibin, 0) = x;
+  (*fPoints)(ibin, 1) = y;
+ }
+ Fit();
+}
+
+
+void AliTMinuitToolkit::FitHistogram(TH2F * his) {
+ //
+ // Fit a two dimensional histogram
+ //
+ fPoints = new TMatrixD((Long64_t)his->GetEntries(), 2);
+ Long64_t entry = 0;
+ for(Int_t ibin=0; ibin < his->GetNbinsX(); ibin++) {
+  Double_t x = his->GetXaxis()->GetBinCenter(ibin);
+  for(Int_t jbin=0; jbin < his->GetNbinsY(); jbin++) {   
+   Long64_t n = his->GetBin(ibin, jbin);
+   Double_t y = his->GetYaxis()->GetBinCenter(jbin);
+   for(Int_t ientries=0; ientries < his->GetBinContent(n); ientries++) {
+    (*fPoints)(entry,0) = x;
+    (*fPoints)(entry,1) = y;
+    entry++;
+   }
+   
+  }
+ }
+
+ Fit();
+}
+
+
+void AliTMinuitToolkit::SetWeightFunction(Char_t * name, Float_t param1, Float_t param2) {
+ //
+ // Set the weight function which must be defined on the interval [0,1].
+ //
+ TString FuncType(name);
+ FuncType.ToUpper();
+ if (FuncType == "EXPONENTIAL") fWeightFunction = new TFormula("exp", Form("TMath::Exp(-TMath::Log(%f)*x)", param1));
+ if (FuncType == "BOX") fWeightFunction = new TFormula("box", Form("TMath::Erfc((x-%f)/0.0001)", param1));
+ if (FuncType == "ERRORFUNCTION") fWeightFunction = new TFormula("err", Form("TMath::Erfc((x-%f)/%f)", param1, param2));// !!!!!!!!!!!!!!!!!
+}
+
+
+void AliTMinuitToolkit::FitterFCN(int &npar, double *dummy, double &fchisq, double *gin, int iflag){
+  //
+  // internal function which gives the specified function to the TMinuit function
+  //  
+
+  // suppress warnings for unused variables:
+  dummy = dummy;
+  iflag = iflag;
+  npar = npar;
+  //
+  AliTMinuitToolkit * fitter = (AliTMinuitToolkit*)TVirtualFitter::GetFitter()->GetObjectFit();
+  fchisq = 0;
+  Int_t nvar       = fitter->GetPoints()->GetNcols()-1;
+  Int_t npoints    = fitter->GetPoints()->GetNrows();
+  
+  // sort points for weighting
+  Double_t sortList[npoints];
+  Int_t indexList[npoints];
+  
+  TVectorD *fWeight = new TVectorD(npoints);
+  
+  for (Int_t ipoint=0; ipoint<npoints; ipoint++){
+    Double_t x[100];
+    for (Int_t ivar=0; ivar<nvar; ivar++){
+      x[ivar] = (*fitter->GetPoints())(ipoint, ivar);      
+    }
+    Float_t funx   = fitter->GetFormula()->EvalPar(x,gin);
+    sortList[ipoint] = TMath::Abs((*fitter->GetPoints())(ipoint, nvar) - funx);
+  }
+  
+  TMath::Sort(npoints, sortList, indexList, false);
+
+  Double_t t;
+  for (Int_t ipoint=0; ipoint<npoints; ipoint++){
+   t = indexList[ipoint]/(Double_t)npoints;
+   (*fWeight)(ipoint) = fitter->GetWeightFunction()->Eval(t);
+  }
+  //
+  // calculate chisquare
+  for (Int_t ipoint=0; ipoint<npoints; ipoint++){
+    Double_t x[100];
+    for (Int_t ivar=0; ivar<nvar; ivar++){
+      x[ivar] = (*fitter->GetPoints())(ipoint, ivar);      
+    }
+    Float_t funx   = fitter->GetFormula()->EvalPar(x,gin);
+    
+    Double_t delta = (*fitter->GetPoints())(ipoint, nvar) - funx;
+    fchisq+= delta*delta*(*fWeight)(ipoint);
+
+  }
+  delete fWeight;
+}
+
+void AliTMinuitToolkit::Fit() {
+  //
+  // internal function that calls the fitter
+  //
+  Int_t nparam = fParam->GetNrows();
+  
+  // set all paramter limits to infinity as default
+  if (fParamLimits == 0) {
+   fParamLimits = new TMatrixD(nparam ,2);
+   for (Int_t iparam=0; iparam<nparam; iparam++){
+    (*fParamLimits)(iparam, 0) = 0;
+    (*fParamLimits)(iparam, 1) = 0;
+   }
+  }
+  
+  // set all weights to 1 as default
+  if (fWeightFunction == 0) {
+   fWeightFunction = new TFormula("constant", "1");
+  }
+  
+  // migrad fit algorithm as default
+  if (fFitAlgorithm == 0) {
+   fFitAlgorithm = "migrad";
+  }
+  
+  // set up the fitter
+  TVirtualFitter *minuit = TVirtualFitter::Fitter(0, nparam);
+  minuit->SetObjectFit(this); 
+  minuit->SetFCN((void*)(AliTMinuitToolkit::FitterFCN));
+  
+  // initialize paramters (step size???)
+  for (Int_t iparam=0; iparam<nparam; iparam++){
+   minuit->SetParameter(iparam, Form("p[%d]",iparam), (*fParam)(iparam), (*fParam)(iparam)/10, (*fParamLimits)(iparam, 0), (*fParamLimits)(iparam, 1));
+  }
+  
+  Double_t argList[2];
+  argList[0] = fMaxCalls; //maximal number of calls 
+  argList[1] = fPrecision; //tolerance normalized to 0.001 
+  if (fMaxCalls == 500 && fPrecision == 1) minuit->ExecuteCommand(fFitAlgorithm, 0, 0); 
+  if (fMaxCalls != 500 || fPrecision != 1) minuit->ExecuteCommand(fFitAlgorithm, argList, 2);
+  // two additional arguments can be specified ExecuteCommand("migrad", argList, 2) - use 0,0 for default
+
+  // fill parameter vector
+  for (Int_t ivar=0; ivar<nparam; ivar++){
+   (*fParam)(ivar) = minuit->GetParameter(ivar);
+  }
+  
+  // fill covariance matrix
+  fCovar = new TMatrixD(nparam, nparam);
+  //TVirtualFitter *fitCov = TVirtualFitter::GetFitter();
+  for(Int_t i=0; i < nparam; i++) {
+   for(Int_t j=0; j < nparam; j++) {
+    (*fCovar)(i,j) = minuit->GetCovarianceMatrixElement(i,j);
+   }
+  }
+}
+
+
+
+void AliTMinuitToolkit::Test() {
+ //
+ // This test function shows the basic working principles of this class 
+ // and illustrates how a robust fit can improve the results
+ //
+ TFormula *FormExp = new TFormula("formExp", "[0]*TMath::Exp(-[1]*x)");
+ SetFitFunction(FormExp);
+ SetFitAlgorithm("simplex");
+ // Set initial values
+ TVectorD *vec1 = new TVectorD(2);
+ (*vec1)(0) = 1800;
+ (*vec1)(1) = 1;
+ SetInitialParam(vec1);
+ //provide some example histogram
+ TH1F * hist = new TH1F("bla", "with (red) and without (black) robust option", 20,0,4);
+ TRandom * rand = new TRandom();
+ for (Int_t i = 0; i < 10000; i++) {
+  hist->Fill(rand->Exp(1));
+  if (i < 1000) hist->Fill(3); //"outliers"
+  if (i < 1070) hist->Fill(3.5);
+  if (i < 670) hist->Fill(2);
+  if (i < 770) hist->Fill(1.5);//"outliers"
+  if (i < 740) hist->Fill(1);
+ }
+ TCanvas * canv = new TCanvas();
+ canv->cd(1);
+ hist->Draw();
+ // fit it with the exponential decay
+ FitHistogram(hist);
+ // draw fit function
+ TF1 *func = new TF1("test", "[0]*TMath::Exp(-[1]*x)", 0, 6);
+ func->SetParameter(0, (*GetParameters())(0));
+ func->SetParameter(1, (*GetParameters())(1));
+ func->Draw("same");
+ // robust fit
+ TVectorD *vec2 = new TVectorD(2);
+ (*vec2)(0) = 1800;
+ (*vec2)(1) = 1;
+ SetInitialParam(vec2);
+ SetWeightFunction("Box", 0.7);
+ FitHistogram(hist);
+ TF1 *func2 = new TF1("test2", "[0]*TMath::Exp(-[1]*x)", 0, 6);
+ func2->SetParameter(0, (*GetParameters())(0));
+ func2->SetParameter(1, (*GetParameters())(1));
+ func2->SetLineColor(kRed);
+ func2->Draw("same");
+}
+
diff --git a/STAT/AliTMinuitToolkit.h b/STAT/AliTMinuitToolkit.h
new file mode 100644 (file)
index 0000000..3986331
--- /dev/null
@@ -0,0 +1,67 @@
+#ifndef AliTMINUITTOOLKIT_H
+#define AliTMINUITTOOLKIT_H
+
+
+#include <TNamed.h>
+#include <TVirtualFitter.h>
+#include <TH1F.h>
+#include <TH2F.h>
+#include <TF1.h>
+#include <TFormula.h>
+#include <TVectorD.h>
+#include <TMatrixD.h>
+#include <TMath.h>
+#include <TString.h>
+
+
+class AliTMinuitToolkit: public TNamed{
+  //
+  //
+  //
+public:
+  AliTMinuitToolkit();
+  virtual ~AliTMinuitToolkit();
+  
+  void FitHistogram(TH1F * his);
+  void FitHistogram(TH2F * his);
+  void SetInitialParam(TVectorD *param) { fParam=param;};
+  void SetParamLimits(TMatrixD *paramLimits) { fParamLimits=paramLimits;};
+  //
+  void SetFitFunction(TFormula * formula) {fFormula=formula;};
+  void SetWeightFunction(TFormula * weightFunction) {fWeightFunction=weightFunction;};
+  void SetWeightFunction(Char_t * name, Float_t param1, Float_t param2 = 0.1);
+  void SetFitAlgorithm(Char_t * name) {fFitAlgorithm=name;};
+  void SetMaxCalls(Int_t calls) {fMaxCalls = calls;};
+  void SetTolerance(Double_t tol) {fPrecision = tol;};
+  void SetPoints(TMatrixD * points) {fPoints = points;};
+  void Test();
+  static void FitterFCN(int &npar, double  *dummy, double &fchisq, double *gin, int iflag);
+  void Fit();
+  //
+  TMatrixD * GetPoints() {return fPoints;};
+  TFormula * GetFormula() {return fFormula;};
+  TVectorD * GetParameters() {return fParam;};
+  TFormula * GetWeightFunction() {return fWeightFunction;};
+  TMatrixD * GetCovarianceMatrix() {return fCovar;};
+private:
+  //
+  TFormula        * fFormula;            // formula of the fitted function
+  TFormula        * fWeightFunction;     // weight function, must be defined between 0 and 1
+  Char_t          * fFitAlgorithm;       // fit algorithm for TMinuit: migrad, simplex, ...
+  //
+  //
+  TMatrixD        * fPoints;             // fitted points
+  TVectorD        * fParam;              // parameter values
+  TMatrixD        * fParamLimits;        // limits of the parameters (up, down)
+  TMatrixD        * fCovar;              // covariance matrix
+  Double_t        * fChi2;               // chi-square
+  Int_t             fMaxCalls;           // maximum number of calls, default value depends on fit algorithm
+  Double_t          fPrecision;          // normalized tolerance, default value depends on fit algorithm
+  //  
+  ClassDef(AliTMinuitToolkit,1);
+};
+
+
+
+#endif
index 43ebeb8..cdcd380 100644 (file)
@@ -17,9 +17,9 @@ OutPutOpt = -o
 
 CXX       = g++
 LD        = g++
-CXXFLAGS  = -O -g -pg
-LDFLAGS   = -O -g -pg
-SOFLAGS   = -shared -g -pg
+CXXFLAGS  = -O -g -pg -fPIC
+LDFLAGS   = -O -g -pg -fPIC
+SOFLAGS   = -shared -g -pg -fPIC
 
 ROOTCXX  = $(shell root-config --cflags)
 ROOTLIBS = $(shell root-config --libs) -lMinuit
index 49e2c0c..cff2356 100644 (file)
@@ -14,6 +14,7 @@
 
 #pragma link C++ class TKDInterpolator+;
 #pragma link C++ class TKDSpline+;
+#pragma link C++ class AliTMinuitToolkit+;
 
 
 #endif
index 6467df6..1fa4156 100644 (file)
@@ -1,4 +1,4 @@
-SRCS=  TKDTree.cxx TKDInterpolator.cxx TKDSpline.cxx
+SRCS=  TKDTree.cxx TKDInterpolator.cxx TKDSpline.cxx AliTMinuitToolkit.cxx
 
 HDRS= $(SRCS:.cxx=.h)