--- /dev/null
+
+#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");
+
+}
+
--- /dev/null
+#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
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff