[u/mrichter/AliRoot.git] / STAT / AliTMinuitToolkit.cxx


#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),
   fMaxCalls(0),
   fPrecision(0)
{
 //
 // standard constructor
 //
 fMaxCalls = 500;
 fPrecision = 1;
}


AliTMinuitToolkit::AliTMinuitToolkit(const AliTMinuitToolkit&) :
   TNamed(),
   fFormula(0),
   fWeightFunction(0),
   fFitAlgorithm(0),
   fPoints(0),
   fParam(0),
   fParamLimits(0),
   fCovar(0),
   fChi2(0),
   fMaxCalls(0),
   fPrecision(0)
{


}


AliTMinuitToolkit& AliTMinuitToolkit::operator=(const AliTMinuitToolkit&) {

 return *this;
}


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 = new Double_t[npoints];
  Int_t  *indexList = new Int_t[npoints];
  
  TVectorD *Weight = 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;
   (*Weight)(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*(*Weight)(ipoint);

  }
  delete Weight;
  delete sortList;
  delete indexList;
}
 

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");
 
}
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),
da7c9ae3	83	fChi2(0),
	84	fMaxCalls(0),
	85	fPrecision(0)
4b23a517	86	{
	87	//
	88	// standard constructor
	89	//
	90	fMaxCalls = 500;
	91	fPrecision = 1;
	92	}
	93
	94
da7c9ae3	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
4b23a517	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
d99c0447	206	Double_t *sortList = new Double_t[npoints];
da7c9ae3	207	Int_t *indexList = new Int_t[npoints];
4b23a517	208
da7c9ae3	209	TVectorD *Weight = new TVectorD(npoints);
4b23a517	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;
da7c9ae3	225	(*Weight)(ipoint) = fitter->GetWeightFunction()->Eval(t);
4b23a517	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;
da7c9ae3	237	fchisq+= deltadelta(*Weight)(ipoint);
4b23a517	238
4b23a517	239	}
da7c9ae3	240	delete Weight;
	241	delete sortList;
	242	delete indexList;
4b23a517	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