Keep track of missing DCS points in DDL maps (flagged by 'x')
[u/mrichter/AliRoot.git] / STAT / AliTMinuitToolkit.cxx
index e63edae..6feb618 100644 (file)
@@ -1,21 +1,36 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/* $Id$ */
 
-#include "AliTMinuitToolkit.h"
-#include <TNamed.h>
-#include <TVirtualFitter.h>
-#include <TH1F.h>
-#include <TH2F.h>
+#include <TCanvas.h>
 #include <TF1.h>
 #include <TFormula.h>
-#include <TVectorD.h>
-#include <TMatrixD.h>
+#include <TH1F.h>
+#include <TH2F.h>
 #include <TMath.h>
-#include <TString.h>
-#include <TROOT.h>
-#include <TCanvas.h>
+#include <TMatrix.h>
 #include <TRandom.h>
+#include <TString.h>
+#include <TVector.h>
+#include <TVectorD.h>
+#include <TVirtualFitter.h>
 
+#include "AliTMinuitToolkit.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.
+// - An option for robust fitting with non-linear functions is implemented.
+//
+// A small example illustrating the usage of AliTMinuitToolkit is given in the function 
+// "AliTMinuitToolkit::Test()".
 // 
 // 
 // 1. Setting the formula:
 //  - 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()".
+//    the command "void Fit()". The weight of each point can be specified
+//    with an n-dimensional vector using "void SetWeights(TVectorD * weights)".
 //
 //
 // 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 
+//  is returned by "TVectorD * GetParameters()". In a similar way the covariance 
 //  matrix of the fit is returned via "TMatrixD * GetCovarianceMatrix()" which
 //  is of the type N x N.
 //
 //  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].
+//  The robust option becomes active if EnableRobust(true, sigma) is called. It is
+//  very much recommended that a normalization value (scale variable) corresponding 
+//  to an expected deviation (sigma) is specified via 
+//  "EnableRobust(Bool_t b, Double_t sigma)".
+//
+//  Performing the fit without knowledge of sigma is also possible if only
+//  "EnableRobust(true)" is activated, but this is NOT RECOMMENDED.
+//
+//  The method is based on another estimator instead of chi^2. For small deviations 
+//  the function behaves like x^2 and for larger deviations like |x| - the so 
+//  called Huber estimator:
+//
+//   h(x) = x^2                              , for x < 2.5*sigma
+//   h(x) = 2*(2.5*sigma)*x - (2.5*sigma)^2  , for x > 2.5*sigma
+//
+//  If a weighting function is specified in addition, a second run with the ordinary 
+//  metric is started, but before entering the iteration every point is weighted 
+//  according to its distance to the outcoming function of the first run. The weighting
+//  function w(x) must be defined on the intervall x in [0,1]. w(0) then 
+//  corresponds to the weight of the closest point and w(1) to the point with the
+//  largest distance.
 //
 //  Some standard weighting functions are predefined in 
 //  "SetWeightFunction(Char_t * name, Float_t param1, Float_t param2 = 0)":
 //   - "EXPONENTIAL" corresponds to "Math::Exp(-TMath::Log(param1)*x)"
 //   - "ERRORFUNCTION" corresponds to "TMath::Erfc((x-param1)/param2)"
 //
-//-------------------------------------------------------------------------
+//
+//  REFERENCE for non-linear robust fitting:
+//  Ekblom H. and Madsen K. (1988), Alogrithms for non-linear Huber estimation,
+//  BIT Numerical Mathematics 29 (1989) 60-76.
+//  internet: http://www.springerlink.com/content/m277218542988344/
+//
+//
+// 5. examples:
+//
+//  A small example illustrating the working principles of AliTMinuitToolkit is given
+//  in the function "AliTMinuitToolkit::Test()".
+//
+//
+//
+// Comments and questions are always welcome: A.Kalweit@gsi.de
+//--------------------------------------------------------------------------------------
 
 
 ClassImp(AliTMinuitToolkit)
@@ -77,18 +129,23 @@ AliTMinuitToolkit::AliTMinuitToolkit() :
    fWeightFunction(0),
    fFitAlgorithm(0),
    fPoints(0),
+   fWeights(0),
    fParam(0),
    fParamLimits(0),
    fCovar(0),
    fChi2(0),
    fMaxCalls(0),
-   fPrecision(0)
+   fPrecision(0),
+   fUseRobust(0),
+   fExpectedSigma(0)
 {
  //
  // standard constructor
  //
  fMaxCalls = 500;
  fPrecision = 1;
+ fUseRobust = false;
+ fExpectedSigma = 0;
 }
 
 
@@ -98,12 +155,15 @@ AliTMinuitToolkit::AliTMinuitToolkit(const AliTMinuitToolkit&) :
    fWeightFunction(0),
    fFitAlgorithm(0),
    fPoints(0),
+   fWeights(0),
    fParam(0),
    fParamLimits(0),
    fCovar(0),
    fChi2(0),
    fMaxCalls(0),
-   fPrecision(0)
+   fPrecision(0),
+   fUseRobust(0),
+   fExpectedSigma(0)
 {
 
 
@@ -122,6 +182,7 @@ AliTMinuitToolkit::~AliTMinuitToolkit(){
   // destructor
   //
   delete fPoints;
+  delete fWeights;
   delete fWeightFunction;
   delete fParamLimits;
   delete fFormula;
@@ -130,7 +191,7 @@ AliTMinuitToolkit::~AliTMinuitToolkit(){
   delete fChi2;
 }
 
-void AliTMinuitToolkit::FitHistogram(TH1F * his) {
+void AliTMinuitToolkit::FitHistogram(TH1F *const his) {
  //
  // Fit a one dimensional histogram
  //
@@ -148,9 +209,9 @@ void AliTMinuitToolkit::FitHistogram(TH1F * his) {
 }
 
 
-void AliTMinuitToolkit::FitHistogram(TH2F * his) {
+void AliTMinuitToolkit::FitHistogram(TH2F *const his) {
  //
- // Fit a two dimensional histogram
+ // Fit a curve to a two dimensional histogram
  //
  fPoints = new TMatrixD((Long64_t)his->GetEntries(), 2);
  Long64_t entry = 0;
@@ -173,7 +234,7 @@ void AliTMinuitToolkit::FitHistogram(TH2F * his) {
 }
 
 
-void AliTMinuitToolkit::SetWeightFunction(Char_t * name, Float_t param1, Float_t param2) {
+void AliTMinuitToolkit::SetWeightFunction(const Char_t *name, Float_t param1, Float_t param2) {
  //
  // Set the weight function which must be defined on the interval [0,1].
  //
@@ -182,7 +243,7 @@ void AliTMinuitToolkit::SetWeightFunction(Char_t * name, Float_t param1, Float_t
  
  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));// !!!!!!!!!!!!!!!!!
+ if (FuncType == "ERRORFUNCTION") fWeightFunction = new TFormula("err", Form("TMath::Erfc((x-%f)/%f)", param1, param2));
  
 }
 
@@ -202,52 +263,49 @@ void AliTMinuitToolkit::FitterFCN(int &npar, double *dummy, double &fchisq, doub
   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++){
+  // calculate mean deviation for normalization or use user-defined sigma
+  Double_t dev = 0.;
+  if (fitter->GetExpectedSigma() == 0 && fitter->GetStatus() == true) {
+   for (Int_t ipoint=0; ipoint<npoints; ipoint++){
     Double_t x[100];
-    for (Int_t ivar=0; ivar<nvar; ivar++){
+     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);
+     }
+    Float_t funx = fitter->GetFormula()->EvalPar(x,gin);
+    Double_t delta = (*fitter->GetPoints())(ipoint, nvar) - funx;
+    dev += TMath::Sqrt(TMath::Abs(delta));
+   }
+   dev = dev/npoints; 
+  } else {
+   dev = fitter->GetExpectedSigma();
   }
-  //
-  // calculate chisquare
+  // 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);
-
+    Float_t funx = fitter->GetFormula()->EvalPar(x,gin);   
+    Double_t delta = TMath::Abs((*fitter->GetPoints())(ipoint, nvar) - funx);
+    if (fitter->GetStatus() == true) {
+     delta = delta/dev; // normalization
+     if (delta <= 2.5) fchisq+= delta*delta; // new metric: Huber-k-estimator
+     if (delta > 2.5) fchisq+= 2*(2.5)*delta - (2.5*2.5);
+    } else {
+     Double_t weight = (*fitter->GetWeights())(ipoint);
+     fchisq+= delta*delta*weight; //old metric
+    }
   }
-  delete Weight;
-  delete sortList;
-  delete indexList;
-}
+ }
  
 
 void AliTMinuitToolkit::Fit() {
   //
   // internal function that calls the fitter
   //
-  Int_t nparam = fParam->GetNrows();
+  Int_t nparam  = fParam->GetNrows();
+  Int_t npoints = fPoints->GetNrows();
+  Int_t nvar    = fPoints->GetNcols()-1;
   
   // set all paramter limits to infinity as default
   if (fParamLimits == 0) {
@@ -259,8 +317,11 @@ void AliTMinuitToolkit::Fit() {
   }
   
   // set all weights to 1 as default
+  Bool_t weightFlag = false;
   if (fWeightFunction == 0) {
    fWeightFunction = new TFormula("constant", "1");
+  } else {
+   weightFlag = true;
   }
   
   // migrad fit algorithm as default
@@ -268,37 +329,81 @@ void AliTMinuitToolkit::Fit() {
    fFitAlgorithm = "migrad";
   }
   
+  // assign weights
+  if (fWeights == 0) {
+   fWeights = new TVectorD(npoints);
+   for (Int_t ipoint=0; ipoint<npoints; ipoint++) (*fWeights)(ipoint) = 1;
+  }
+  
   // set up the fitter
   TVirtualFitter *minuit = TVirtualFitter::Fitter(0, nparam);
   minuit->SetObjectFit(this); 
-  minuit->SetFCN((void*)(AliTMinuitToolkit::FitterFCN));
+  minuit->SetFCN(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);
+   fFormula->SetParameter(ivar, minuit->GetParameter(ivar));
+  }
 
+  // if a weight function is specified -> enter 2nd run with weights
+  if (weightFlag == true && fUseRobust == true) {
+   // sort points for weighting
+   Double_t *sortList = new Double_t[npoints];
+   Int_t  *indexList = new Int_t[npoints];   
+   for (Int_t ipoint=0; ipoint<npoints; ipoint++){
+    Double_t funx = fFormula->Eval((*fPoints)(ipoint, 0));
+    Double_t delta = TMath::Abs((*fPoints)[ipoint][nvar] - funx);
+    sortList[ipoint] = delta;
+   } 
+   TMath::Sort(npoints, sortList, indexList, false);
+   for (Int_t ip=0; ip<npoints; ip++){
+    Double_t t = ip/(Double_t)npoints;
+    (*fWeights)(indexList[ip]) = fWeightFunction->Eval(t);
+   }
+   
+   // set up the fitter
+   fUseRobust = false;
+   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));
+   }
+   // start fitting
+   if (fMaxCalls == 500 && fPrecision == 1) minuit->ExecuteCommand(fFitAlgorithm, 0, 0); 
+   if (fMaxCalls != 500 || fPrecision != 1) minuit->ExecuteCommand(fFitAlgorithm, argList, 2);
+   fUseRobust = true;
+   
+   delete sortList; 
+   delete indexList;    
+  }
+  
   // fill parameter vector
   for (Int_t ivar=0; ivar<nparam; ivar++){
    (*fParam)(ivar) = minuit->GetParameter(ivar);
+   fFormula->SetParameter(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);
    }
   }
+  
+  if (weightFlag == false) fWeightFunction = 0;
 }
 
 
@@ -308,16 +413,9 @@ 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);
+ // 1. provide some example histogram
+ TH1F * hist = new TH1F("test", "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));
@@ -330,23 +428,33 @@ void AliTMinuitToolkit::Test() {
  TCanvas * canv = new TCanvas();
  canv->cd(1);
  hist->Draw();
- // fit it with the exponential decay
- FitHistogram(hist);
+ // 2. example fit without robust option
+ AliTMinuitToolkit * tool = new AliTMinuitToolkit();
+ TFormula *aFormExp = new TFormula("formExp", "[0]*TMath::Exp(-[1]*x)");
+ tool->SetFitFunction(aFormExp);
+ TVectorD *vec1 = new TVectorD(2); // Set initial values
+ (*vec1)(0) = 1800;
+ (*vec1)(1) = 1;
+ tool->SetInitialParam(vec1);
+ tool->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->SetParameters((*tool->GetParameters())(0), (*tool->GetParameters())(1));
  func->Draw("same");
- // robust fit
+ // 3 . robust fit 
  TVectorD *vec2 = new TVectorD(2);
  (*vec2)(0) = 1800;
  (*vec2)(1) = 1;
- SetInitialParam(vec2);
- SetWeightFunction("Box", 0.7);
- FitHistogram(hist);
+ tool->SetInitialParam(vec2);
+ tool->EnableRobust(true, 10);
+ tool->SetWeightFunction("box", 0.75);
+ tool->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->SetParameter(0, (*tool->GetParameters())(0));
+ func2->SetParameter(1, (*tool->GetParameters())(1));
  func2->SetLineColor(kRed);
  func2->Draw("same");