refactoring multiplicity analysis
authorjgrosseo <jgrosseo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 25 Feb 2009 09:59:20 +0000 (09:59 +0000)
committerjgrosseo <jgrosseo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 25 Feb 2009 09:59:20 +0000 (09:59 +0000)
new class created that does the pure unfolding
adapted the analysis accordingly

PWG0/AliUnfolding.cxx [new file with mode: 0644]
PWG0/AliUnfolding.h
PWG0/PWG0baseLinkDef.h
PWG0/libPWG0base.pkg
PWG0/multiplicity/AliMultiplicityCorrection.cxx
PWG0/multiplicity/AliMultiplicityCorrection.h
PWG0/multiplicity/correct.C

diff --git a/PWG0/AliUnfolding.cxx b/PWG0/AliUnfolding.cxx
new file mode 100644 (file)
index 0000000..6915330
--- /dev/null
@@ -0,0 +1,1103 @@
+/**************************************************************************
+ * 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: AliUnfolding.cxx 31168 2009-02-23 15:18:45Z jgrosseo $ */
+
+// This class allows 1-dimensional unfolding.
+// Methods that are implemented are chi2 minimization and bayesian unfolding.
+//
+//  Author: Jan.Fiete.Grosse-Oetringhaus@cern.ch
+
+#include "AliUnfolding.h"
+#include <TH1F.h>
+#include <TH2F.h>
+#include <TVirtualFitter.h>
+#include <TMath.h>
+#include <TCanvas.h>
+#include <TF1.h>
+
+TMatrixD* AliUnfolding::fgCorrelationMatrix = 0;
+TMatrixD* AliUnfolding::fgCorrelationCovarianceMatrix = 0;
+TVectorD* AliUnfolding::fgCurrentESDVector = 0;
+TVectorD* AliUnfolding::fgEntropyAPriori = 0;
+
+TF1* AliUnfolding::fgFitFunction = 0;
+
+AliUnfolding::MethodType AliUnfolding::fgMethodType = AliUnfolding::kInvalid;
+Int_t AliUnfolding::fgMaxInput  = -1;  // bins in measured histogram
+Int_t AliUnfolding::fgMaxParams = -1;  // bins in unfolded histogram = number of fit params
+Float_t AliUnfolding::fgOverflowBinLimit = -1;
+
+AliUnfolding::RegularizationType AliUnfolding::fgRegularizationType = AliUnfolding::kPol1;
+Float_t AliUnfolding::fgRegularizationWeight = 10000;
+Int_t AliUnfolding::fgSkipBinsBegin = 0;
+Float_t AliUnfolding::fgMinuitStepSize = 0.1;                 // (usually not needed to be changed) step size in minimization
+Bool_t AliUnfolding::fgNormalizeInput = kFALSE;                  // normalize input spectrum
+
+Float_t AliUnfolding::fgBayesianSmoothing  = 1;           // smoothing parameter (0 = no smoothing)
+Int_t   AliUnfolding::fgBayesianIterations = 10;         // number of iterations in Bayesian method
+
+Bool_t AliUnfolding::fgDebug = kFALSE;
+
+ClassImp(AliUnfolding)
+
+//____________________________________________________________________
+void AliUnfolding::SetUnfoldingMethod(MethodType methodType)
+{
+  // set unfolding method
+  fgMethodType = methodType; 
+  
+  const char* name = 0;
+  switch (methodType)
+  {
+    case kInvalid: name = "INVALID"; break;
+    case kChi2Minimization: name = "Chi2 Minimization"; break;
+    case kBayesian: name = "Bayesian unfolding"; break;
+    case kFunction: name = "Functional fit"; break;
+  }
+  Printf("AliUnfolding::SetUnfoldingMethod: %s enabled.", name);
+}
+
+//____________________________________________________________________
+void AliUnfolding::SetCreateOverflowBin(Float_t overflowBinLimit) 
+{ 
+  // enable the creation of a overflow bin that includes all statistics below the given limit
+  
+  fgOverflowBinLimit = overflowBinLimit; 
+  
+  Printf("AliUnfolding::SetCreateOverflowBin: overflow bin limit set to %f", overflowBinLimit);
+}
+
+//____________________________________________________________________
+void AliUnfolding::SetSkipBinsBegin(Int_t nBins)
+{
+  // set number of skipped bins in regularization
+  
+  fgSkipBinsBegin = nBins;
+  
+  Printf("AliUnfolding::SetSkipBinsBegin: skipping %d bins at the beginning of the spectrum in the regularization.", fgSkipBinsBegin);
+}
+
+//____________________________________________________________________
+void AliUnfolding::SetNbins(Int_t nMeasured, Int_t nUnfolded) 
+{ 
+  // set number of bins in the input (measured) distribution and in the unfolded distribution
+  fgMaxInput = nMeasured; 
+  fgMaxParams = nUnfolded; 
+  
+  Printf("AliUnfolding::SetNbins: Set %d measured bins and %d unfolded bins", nMeasured, nUnfolded);
+}
+
+//____________________________________________________________________
+void AliUnfolding::SetChi2Regularization(RegularizationType type, Float_t weight)
+{
+  //
+  // sets the parameters for chi2 minimization
+  //
+
+  fgRegularizationType = type;
+  fgRegularizationWeight = weight;
+
+  Printf("AliUnfolding::SetChi2Regularization --> Regularization set to %d with weight %f", (Int_t) type, weight);
+}
+
+//____________________________________________________________________
+void AliUnfolding::SetBayesianParameters(Float_t smoothing, Int_t nIterations)
+{
+  //
+  // sets the parameters for Bayesian unfolding
+  //
+
+  fgBayesianSmoothing = smoothing;
+  fgBayesianIterations = nIterations;
+
+  Printf("AliUnfolding::SetBayesianParameters --> Paramaters set to %d iterations with smoothing %f", fgBayesianIterations, fgBayesianSmoothing);
+}
+
+//____________________________________________________________________
+void AliUnfolding::SetFunction(TF1* function)
+{
+  // set function for unfolding with a fit function
+  
+  fgFitFunction = function;
+  
+  Printf("AliUnfolding::SetFunction: Set fit function with %d parameters.", function->GetNpar());
+}
+
+//____________________________________________________________________
+Int_t AliUnfolding::Unfold(TH2* correlation, TH1* efficiency, TH1* measured, TH1* initialConditions, TH1* result, Bool_t check)
+{
+  // unfolds with unfolding method fgMethodType
+  //
+  // parameters:
+  //  correlation: response matrix as measured vs. generated
+  //  efficiency:  (optional) efficiency that is applied on the unfolded spectrum, i.e. it has to be in unfolded variables. If 0 no efficiency is applied.
+  //  measured:    the measured spectrum
+  //  initialConditions: (optional) initial conditions for the unfolding. if 0 the measured spectrum is used as initial conditions.
+  //  result:      target for the unfolded result
+  //  check:       depends on the unfolding method, see comments in specific functions
+
+  if (fgMaxInput == -1)
+  {
+    Printf("AliUnfolding::Unfold: WARNING. Number of measured bins not set with SetNbins. Using number of bins in measured distribution");
+    fgMaxInput = measured->GetNbinsX();
+  }
+  if (fgMaxParams == -1)
+  {
+    Printf("AliUnfolding::Unfold: WARNING. Number of unfolded bins not set with SetNbins. Using number of bins in measured distribution");
+    fgMaxParams = measured->GetNbinsX();
+  }
+
+  if (fgOverflowBinLimit > 0)
+    CreateOverflowBin(correlation, measured);
+    
+  switch (fgMethodType)
+  {
+    case kInvalid:
+    {
+      Printf("AliUnfolding::Unfold: ERROR: Unfolding method not set. Use SetUnfoldingMethod. Exiting...");
+      return -1;
+    }
+    case kChi2Minimization:
+      return UnfoldWithMinuit(correlation, efficiency, measured, initialConditions, result, check);
+    case kBayesian:
+      return UnfoldWithBayesian(correlation, efficiency, measured, initialConditions, result);
+    case kFunction:
+      return UnfoldWithFunction(correlation, efficiency, measured, initialConditions, result);
+  }
+  return -1;
+}
+
+//____________________________________________________________________
+void AliUnfolding::SetStaticVariables(TH2* correlation, TH1* measured)
+{
+  // fill static variables needed for minuit fit
+
+  if (!fgCorrelationMatrix)
+    fgCorrelationMatrix = new TMatrixD(fgMaxInput, fgMaxParams);
+  if (!fgCorrelationCovarianceMatrix)
+    fgCorrelationCovarianceMatrix = new TMatrixD(fgMaxInput, fgMaxInput);
+  if (!fgCurrentESDVector)
+    fgCurrentESDVector = new TVectorD(fgMaxInput);
+  if (!fgEntropyAPriori)
+    fgEntropyAPriori = new TVectorD(fgMaxParams);
+
+  fgCorrelationMatrix->Zero();
+  fgCorrelationCovarianceMatrix->Zero();
+  fgCurrentESDVector->Zero();
+  fgEntropyAPriori->Zero();
+
+  // normalize correction for given nPart
+  for (Int_t i=1; i<=correlation->GetNbinsX(); ++i)
+  {
+    Double_t sum = correlation->Integral(i, i, 1, correlation->GetNbinsY());
+    if (sum <= 0)
+      continue;
+    Float_t maxValue = 0;
+    Int_t maxBin = -1;
+    for (Int_t j=1; j<=correlation->GetNbinsY(); ++j)
+    {
+      // find most probably value
+      if (maxValue < correlation->GetBinContent(i, j))
+      {
+        maxValue = correlation->GetBinContent(i, j);
+        maxBin = j;
+      }
+
+      // npart sum to 1
+      correlation->SetBinContent(i, j, correlation->GetBinContent(i, j) / sum);
+      correlation->SetBinError(i, j, correlation->GetBinError(i, j) / sum);
+
+      if (i <= fgMaxParams && j <= fgMaxInput)
+        (*fgCorrelationMatrix)(j-1, i-1) = correlation->GetBinContent(i, j);
+    }
+
+    //printf("MPV for Ntrue = %f is %f\n", fCurrentCorrelation->GetXaxis()->GetBinCenter(i), fCurrentCorrelation->GetYaxis()->GetBinCenter(maxBin));
+  }
+    
+  //normalize measured
+  if (fgNormalizeInput)
+    measured->Scale(1.0 / measured->Integral());
+  
+  for (Int_t i=0; i<fgMaxInput; ++i)
+  {
+    (*fgCurrentESDVector)[i] = measured->GetBinContent(i+1);
+    if (measured->GetBinError(i+1) > 0)
+      (*fgCorrelationCovarianceMatrix)(i, i) = (Double_t) 1e-6 / measured->GetBinError(i+1) / measured->GetBinError(i+1);
+
+    if ((*fgCorrelationCovarianceMatrix)(i, i) > 1e7)
+      (*fgCorrelationCovarianceMatrix)(i, i) = 0;
+    //Printf("%d, %e", i, (*fgCorrelationCovarianceMatrix)(i, i));
+  }
+}
+
+//____________________________________________________________________
+Int_t AliUnfolding::UnfoldWithMinuit(TH2* correlation, TH1* efficiency, TH1* measured, TH1* initialConditions, TH1* result, Bool_t check)
+{
+  //
+  // implementation of unfolding (internal function)
+  //
+  // unfolds <measured> using response from <correlation> and effiency <efficiency>
+  // output is in <result>
+  // <initialConditions> set the initial values for the minimization, if 0 <measured> is used
+  // if <check> is true no unfolding is made, instead only the chi2 without unfolding is printed
+  //
+  // returns minuit status (0 = success), (-1 when check was set)
+  //
+
+  SetStaticVariables(correlation, measured);
+  
+  // Initialize TMinuit via generic fitter interface
+  TVirtualFitter *minuit = TVirtualFitter::Fitter(0, fgMaxParams);
+  Double_t arglist[100];
+
+  // disable any output (-1), unfortuantly we do not see warnings anymore then. Have to find another way...
+  arglist[0] = 0;
+  minuit->ExecuteCommand("SET PRINT", arglist, 1);
+
+  // however, enable warnings
+  //minuit->ExecuteCommand("SET WAR", arglist, 0);
+
+  // set minimization function
+  minuit->SetFCN(Chi2Function);
+
+  for (Int_t i=0; i<fgMaxParams; i++)
+    (*fgEntropyAPriori)[i] = 1;
+
+  if (!initialConditions) {
+    initialConditions = measured;
+  } else {
+    Printf("AliUnfolding::UnfoldWithMinuit: Using different initial conditions...");
+    //new TCanvas; initialConditions->DrawCopy();
+    if (fgNormalizeInput)
+      initialConditions->Scale(1.0 / initialConditions->Integral());
+  }
+
+  // set initial conditions as a-priori distribution for MRX regularization
+  for (Int_t i=0; i<fgMaxParams; i++)
+    if (initialConditions->GetBinContent(i+1) > 0)
+      (*fgEntropyAPriori)[i] = initialConditions->GetBinContent(i+1);
+
+  Double_t* results = new Double_t[fgMaxParams+1];
+  for (Int_t i=0; i<fgMaxParams; ++i)
+  {
+    results[i] = initialConditions->GetBinContent(i+1);
+
+    // minuit sees squared values to prevent it from going negative...
+    results[i] = TMath::Sqrt(results[i]);
+
+    minuit->SetParameter(i, Form("param%d", i), results[i], fgMinuitStepSize, 0, 0);
+  }
+  
+  Int_t dummy = 0;
+  Double_t chi2 = 0;
+  Chi2Function(dummy, 0, chi2, results, 0);
+  printf("AliUnfolding::UnfoldWithMinuit: Chi2 of initial parameters is = %f\n", chi2);
+
+  if (check)
+  {
+    DrawGuess(results);
+    delete[] results;
+    return -1;
+  }
+
+  // first param is number of iterations, second is precision....
+  arglist[0] = 1e6;
+  //arglist[1] = 1e-5;
+  //minuit->ExecuteCommand("SCAN", arglist, 0);
+  Int_t status = minuit->ExecuteCommand("MIGRAD", arglist, 1);
+  Printf("AliUnfolding::UnfoldWithMinuit: MINUIT status is %d", status);
+  //printf("!!!!!!!!!!!!!! MIGRAD finished: Starting MINOS !!!!!!!!!!!!!!");
+  //minuit->ExecuteCommand("MINOS", arglist, 0);
+
+  for (Int_t i=0; i<fgMaxParams; ++i)
+  {
+    results[i] = minuit->GetParameter(i);
+    Double_t value = results[i] * results[i];
+    // error is : (relError) * (value) = (minuit->GetParError(i) / minuit->GetParameter(i)) * (minuit->GetParameter(i) * minuit->GetParameter(i))
+    Double_t error = minuit->GetParError(i) * results[i];
+    
+    if (efficiency)
+    {
+      if (efficiency->GetBinContent(i+1) > 0)
+      {
+        value /= efficiency->GetBinContent(i+1);
+        error /= efficiency->GetBinContent(i+1);
+      }
+      else
+      {
+        value = 0;
+        error = 0;
+      }
+    }
+    
+    result->SetBinContent(i+1, value);
+    result->SetBinError(i+1, error);
+  }
+  
+  if (fgDebug)
+    DrawGuess(results);
+
+  delete[] results;
+
+  return status;
+}
+
+//____________________________________________________________________
+Int_t AliUnfolding::UnfoldWithBayesian(TH2* correlation, TH1* aEfficiency, TH1* measured, TH1* initialConditions, TH1* aResult)
+{
+  //
+  // unfolds a spectrum using the Bayesian method
+  //
+  
+  if (measured->Integral() <= 0)
+  {
+    Printf("AliUnfolding::UnfoldWithBayesian: ERROR: The measured spectrum is empty");
+    return -1;
+  }
+
+  const Int_t kStartBin = 0;
+
+  Int_t kMaxM = fgMaxInput;  //<= fCurrentCorrelation->GetNbinsY(); // max measured axis
+  Int_t kMaxT = fgMaxParams; //<= fCurrentCorrelation->GetNbinsX(); // max true axis
+
+  // convergence limit: kMaxT * 0.001^2 = kMaxT * 1e-6 (e.g. 250 bins --> 2.5 e-4)
+  const Double_t kConvergenceLimit = kMaxT * 1e-6;
+
+  // store information in arrays, to increase processing speed (~ factor 5)
+  Double_t* measuredCopy = new Double_t[kMaxM];
+  Double_t* measuredError = new Double_t[kMaxM];
+  Double_t* prior = new Double_t[kMaxT];
+  Double_t* result = new Double_t[kMaxT];
+  Double_t* efficiency = new Double_t[kMaxT];
+
+  Double_t** response = new Double_t*[kMaxT];
+  Double_t** inverseResponse = new Double_t*[kMaxT];
+  for (Int_t i=0; i<kMaxT; i++)
+  {
+    response[i] = new Double_t[kMaxM];
+    inverseResponse[i] = new Double_t[kMaxM];
+  }
+
+  // for normalization
+  Float_t measuredIntegral = measured->Integral();
+  for (Int_t m=0; m<kMaxM; m++)
+  {
+    measuredCopy[m] = measured->GetBinContent(m+1) / measuredIntegral;
+    measuredError[m] = measured->GetBinError(m+1) / measuredIntegral;
+
+    for (Int_t t=0; t<kMaxT; t++)
+    {
+      response[t][m] = correlation->GetBinContent(t+1, m+1);
+      inverseResponse[t][m] = 0;
+    }
+  }
+
+  for (Int_t t=0; t<kMaxT; t++)
+  {
+    if (efficiency)
+    {
+      efficiency[t] = aEfficiency->GetBinContent(t+1);
+    }
+    else
+      efficiency[t] = 1;
+      
+    prior[t] = measuredCopy[t];
+    result[t] = 0;
+  }
+
+  // pick prior distribution
+  if (initialConditions)
+  {
+    printf("Using different starting conditions...\n");
+    // for normalization
+    Float_t inputDistIntegral = initialConditions->Integral();
+    for (Int_t i=0; i<kMaxT; i++)
+      prior[i] = initialConditions->GetBinContent(i+1) / inputDistIntegral;
+  }
+
+  //TH1F* convergence = new TH1F("convergence", "convergence", 200, 0.5, 200.5);
+  
+  // unfold...
+  for (Int_t i=0; i<fgBayesianIterations || fgBayesianIterations< 0; i++)
+  {
+    if (fgDebug)
+      Printf("AliUnfolding::UnfoldWithBayesian: iteration %i", i);
+
+    // calculate IR from Beyes theorem
+    // IR_ji = R_ij * prior_i / sum_k(R_kj * prior_k)
+
+    Double_t chi2Measured = 0;
+    for (Int_t m=0; m<kMaxM; m++)
+    {
+      Float_t norm = 0;
+      for (Int_t t = kStartBin; t<kMaxT; t++)
+        norm += response[t][m] * prior[t];
+
+      // calc. chi2: (measured - response * prior) / error
+      if (measuredError[m] > 0)
+      {
+        Double_t value = (measuredCopy[m] - norm) / measuredError[m];
+        chi2Measured += value * value;
+      }
+
+      if (norm > 0)
+      {
+        for (Int_t t = kStartBin; t<kMaxT; t++)
+          inverseResponse[t][m] = response[t][m] * prior[t] / norm;
+      }
+      else
+      {
+        for (Int_t t = kStartBin; t<kMaxT; t++)
+          inverseResponse[t][m] = 0;
+      }
+    }
+    //Printf("chi2Measured of the last prior is %e", chi2Measured);
+
+    for (Int_t t = kStartBin; t<kMaxT; t++)
+    {
+      Float_t value = 0;
+      for (Int_t m=0; m<kMaxM; m++)
+        value += inverseResponse[t][m] * measuredCopy[m];
+
+      if (efficiency[t] > 0)
+        result[t] = value / efficiency[t];
+      else
+        result[t] = 0;
+    }
+    
+    // draw intermediate result
+    /*
+    for (Int_t t=0; t<kMaxT; t++)
+      aResult->SetBinContent(t+1, result[t]);
+    aResult->SetMarkerStyle(20+i);
+    aResult->SetMarkerColor(2);
+    aResult->DrawCopy("P SAME HIST");
+    */
+
+    Double_t chi2LastIter = 0;
+    // regularization (simple smoothing)
+    for (Int_t t=kStartBin; t<kMaxT; t++)
+    {
+      Float_t newValue = 0;
+      
+      // 0 bin excluded from smoothing
+      if (t > kStartBin+2 && t<kMaxT-1)
+      {
+        Float_t average = (result[t-1] + result[t] + result[t+1]) / 3;
+
+        // weight the average with the regularization parameter
+        newValue = (1 - fgBayesianSmoothing) * result[t] + fgBayesianSmoothing * average;
+      }
+      else
+        newValue = result[t];
+
+      // calculate chi2 (change from last iteration)
+      if (prior[t] > 1e-5)
+      {
+        Double_t diff = (prior[t] - newValue) / prior[t];
+        chi2LastIter += diff * diff;
+      }
+
+      prior[t] = newValue;
+    }
+    //printf("Chi2 of %d iteration = %e\n", i, chi2LastIter);
+    //convergence->Fill(i+1, chi2LastIter);
+
+    if (fgBayesianIterations < 0 && chi2LastIter < kConvergenceLimit)
+    {
+      Printf("AliUnfolding::UnfoldWithBayesian: Stopped Bayesian unfolding after %d iterations at chi2(change since last iteration) of %e; chi2Measured of the last prior is %e", i, chi2LastIter, chi2Measured);
+      break;
+    }
+  } // end of iterations
+
+  //new TCanvas; convergence->DrawCopy(); gPad->SetLogy();
+  //delete convergence;
+
+  for (Int_t t=0; t<kMaxT; t++)
+    aResult->SetBinContent(t+1, result[t]);
+
+  delete[] measuredCopy;
+  delete[] measuredError;
+  delete[] prior;
+  delete[] result;
+  delete[] efficiency;
+
+  for (Int_t i=0; i<kMaxT; i++)
+  {
+    delete[] response[i];
+    delete[] inverseResponse[i];
+  }
+  delete[] response;
+  delete[] inverseResponse;
+  
+  return 0;
+
+  // ********
+  // Calculate the covariance matrix, all arguments are taken from NIM,A362,487-498,1995
+
+  /*printf("Calculating covariance matrix. This may take some time...\n");
+
+  // check if this is the right one...
+  TH1* sumHist = GetMultiplicityMC(inputRange, eventType)->ProjectionY("sumHist", 1, GetMultiplicityMC(inputRange, eventType)->GetNbinsX());
+
+  Int_t xBins = hInverseResponseBayes->GetNbinsX();
+  Int_t yBins = hInverseResponseBayes->GetNbinsY();
+
+  // calculate "unfolding matrix" Mij
+  Float_t matrixM[251][251];
+  for (Int_t i=1; i<=xBins; i++)
+  {
+    for (Int_t j=1; j<=yBins; j++)
+    {
+      if (fCurrentEfficiency->GetBinContent(i) > 0)
+        matrixM[i-1][j-1] = hInverseResponseBayes->GetBinContent(i, j) / fCurrentEfficiency->GetBinContent(i);
+      else
+        matrixM[i-1][j-1] = 0;
+    }
+  }
+
+  Float_t* vectorn = new Float_t[yBins];
+  for (Int_t j=1; j<=yBins; j++)
+    vectorn[j-1] = fCurrentESD->GetBinContent(j);
+
+  // first part of covariance matrix, depends on input distribution n(E)
+  Float_t cov1[251][251];
+
+  Float_t nEvents = fCurrentESD->Integral(); // N
+
+  xBins = 20;
+  yBins = 20;
+
+  for (Int_t k=0; k<xBins; k++)
+  {
+    printf("In Cov1: %d\n", k);
+    for (Int_t l=0; l<yBins; l++)
+    {
+      cov1[k][l] = 0;
+
+      // sum_j Mkj Mlj n(Ej) * (1 - n(Ej) / N)
+      for (Int_t j=0; j<yBins; j++)
+        cov1[k][l] += matrixM[k][j] * matrixM[l][j] * vectorn[j]
+          * (1.0 - vectorn[j] / nEvents);
+
+      // - sum_i,j (i != j) Mki Mlj n(Ei) n(Ej) / N
+      for (Int_t i=0; i<yBins; i++)
+        for (Int_t j=0; j<yBins; j++)
+        {
+          if (i == j)
+            continue;
+          cov1[k][l] -= matrixM[k][i] * matrixM[l][j] * vectorn[i]
+            * vectorn[j] / nEvents;
+         }
+    }
+  }
+
+  printf("Cov1 finished\n");
+
+  TH2F* cov = (TH2F*) hInverseResponseBayes->Clone("cov");
+  cov->Reset();
+
+  for (Int_t i=1; i<=xBins; i++)
+    for (Int_t j=1; j<=yBins; j++)
+      cov->SetBinContent(i, j, cov1[i-1][j-1]);
+
+  new TCanvas;
+  cov->Draw("COLZ");
+
+  // second part of covariance matrix, depends on response matrix
+  Float_t cov2[251][251];
+
+  // Cov[P(Er|Cu), P(Es|Cu)] term
+  Float_t covTerm[100][100][100];
+  for (Int_t r=0; r<yBins; r++)
+    for (Int_t u=0; u<xBins; u++)
+      for (Int_t s=0; s<yBins; s++)
+      {
+        if (r == s)
+          covTerm[r][u][s] = 1.0 / sumHist->GetBinContent(u+1) * hResponse->GetBinContent(u+1, r+1)
+            * (1.0 - hResponse->GetBinContent(u+1, r+1));
+        else
+          covTerm[r][u][s] = - 1.0 / sumHist->GetBinContent(u+1) * hResponse->GetBinContent(u+1, r+1)
+            * hResponse->GetBinContent(u+1, s+1);
+      }
+
+  for (Int_t k=0; k<xBins; k++)
+    for (Int_t l=0; l<yBins; l++)
+    {
+      cov2[k][l] = 0;
+      printf("In Cov2: %d %d\n", k, l);
+      for (Int_t i=0; i<yBins; i++)
+        for (Int_t j=0; j<yBins; j++)
+        {
+          //printf("In Cov2: %d %d %d %d\n", k, l, i, j);
+          // calculate Cov(Mki, Mlj) = sum{ru},{su} ...
+          Float_t tmpCov = 0;
+          for (Int_t r=0; r<yBins; r++)
+            for (Int_t u=0; u<xBins; u++)
+              for (Int_t s=0; s<yBins; s++)
+              {
+                if (hResponse->GetBinContent(u+1, r+1) == 0 || hResponse->GetBinContent(u+1, s+1) == 0
+                  || hResponse->GetBinContent(u+1, i+1) == 0)
+                  continue;
+
+                tmpCov += BayesCovarianceDerivate(matrixM, hResponse, fCurrentEfficiency, k, i, r, u)
+                  * BayesCovarianceDerivate(matrixM, hResponse, fCurrentEfficiency, l, j, s, u)
+                  * covTerm[r][u][s];
+              }
+
+          cov2[k][l] += fCurrentESD->GetBinContent(i+1) * fCurrentESD->GetBinContent(j+1) * tmpCov;
+        }
+    }
+
+  printf("Cov2 finished\n");
+
+  for (Int_t i=1; i<=xBins; i++)
+    for (Int_t j=1; j<=yBins; j++)
+      cov->SetBinContent(i, j, cov1[i-1][j-1] + cov2[i-1][j-1]);
+
+  new TCanvas;
+  cov->Draw("COLZ");*/
+}
+
+//____________________________________________________________________
+Double_t AliUnfolding::RegularizationPol0(TVectorD& params)
+{
+  // homogenity term for minuit fitting
+  // pure function of the parameters
+  // prefers constant function (pol0)
+
+  Double_t chi2 = 0;
+
+  for (Int_t i=1+fgSkipBinsBegin; i<fgMaxParams; ++i)
+  {
+    Double_t right  = params[i];
+    Double_t left   = params[i-1];
+
+    if (right != 0)
+    {
+      Double_t diff = 1 - left / right;
+      chi2 += diff * diff;
+    }
+  }
+
+  return chi2 / 100.0;
+}
+
+//____________________________________________________________________
+Double_t AliUnfolding::RegularizationPol1(TVectorD& params)
+{
+  // homogenity term for minuit fitting
+  // pure function of the parameters
+  // prefers linear function (pol1)
+
+  Double_t chi2 = 0;
+
+  for (Int_t i=2+fgSkipBinsBegin; i<fgMaxParams; ++i)
+  {
+    if (params[i-1] == 0)
+      continue;
+
+    Double_t right  = params[i];
+    Double_t middle = params[i-1];
+    Double_t left   = params[i-2];
+
+    Double_t der1 = (right - middle);
+    Double_t der2 = (middle - left);
+
+    Double_t diff = (der1 - der2) / middle;
+
+    chi2 += diff * diff;
+  }
+
+  return chi2;
+}
+
+//____________________________________________________________________
+Double_t AliUnfolding::RegularizationLog(TVectorD& params)
+{
+  // homogenity term for minuit fitting
+  // pure function of the parameters
+  // prefers linear function (pol1)
+
+  Double_t chi2 = 0;
+
+  for (Int_t i=2+fgSkipBinsBegin; i<fgMaxParams; ++i)
+  {
+    if (params[i-1] == 0)
+      continue;
+
+    Double_t right  = log(params[i]);
+    Double_t middle = log(params[i-1]);
+    Double_t left   = log(params[i-2]);
+
+    Double_t der1 = (right - middle);
+    Double_t der2 = (middle - left);
+
+    Double_t diff = (der1 - der2) / middle;
+
+    chi2 += diff * diff;
+  }
+
+  return chi2 * 100;
+}
+
+//____________________________________________________________________
+Double_t AliUnfolding::RegularizationTotalCurvature(TVectorD& params)
+{
+  // homogenity term for minuit fitting
+  // pure function of the parameters
+  // minimizes the total curvature (from Unfolding Methods In High-Energy Physics Experiments,
+  // V. Blobel (Hamburg U.) . DESY 84/118, Dec 1984. 40pp.
+
+  Double_t chi2 = 0;
+
+  for (Int_t i=2+fgSkipBinsBegin; i<fgMaxParams; ++i)
+  {
+    Double_t right  = params[i];
+    Double_t middle = params[i-1];
+    Double_t left   = params[i-2];
+
+    Double_t der1 = (right - middle);
+    Double_t der2 = (middle - left);
+
+    Double_t diff = (der1 - der2);
+
+    chi2 += diff * diff;
+  }
+
+  return chi2 * 1e4;
+}
+
+//____________________________________________________________________
+Double_t AliUnfolding::RegularizationEntropy(TVectorD& params)
+{
+  // homogenity term for minuit fitting
+  // pure function of the parameters
+  // calculates entropy, from
+  // The method of reduced cross-entropy (M. Schmelling 1993)
+
+  Double_t paramSum = 0;
+  
+  for (Int_t i=fgSkipBinsBegin; i<fgMaxParams; ++i)
+    paramSum += params[i];
+
+  Double_t chi2 = 0;
+  for (Int_t i=fgSkipBinsBegin; i<fgMaxParams; ++i)
+  {
+    //Double_t tmp = params[i] / paramSum;
+    Double_t tmp = params[i];
+    if (tmp > 0 && (*fgEntropyAPriori)[i] > 0)
+    {
+      chi2 += tmp * TMath::Log(tmp / (*fgEntropyAPriori)[i]);
+    }
+  }
+
+  return 100.0 + chi2;
+}
+
+//____________________________________________________________________
+void AliUnfolding::Chi2Function(Int_t&, Double_t*, Double_t& chi2, Double_t *params, Int_t)
+{
+  //
+  // fit function for minuit
+  // does: (m - Ad)W(m - Ad) where m = measured, A correlation matrix, d = guess, W = covariance matrix
+  //
+
+  // d
+  static TVectorD paramsVector(fgMaxParams);
+  for (Int_t i=0; i<fgMaxParams; ++i)
+    paramsVector[i] = params[i] * params[i];
+
+  // calculate penalty factor
+  Double_t penaltyVal = 0;
+  switch (fgRegularizationType)
+  {
+    case kNone:       break;
+    case kPol0:       penaltyVal = RegularizationPol0(paramsVector); break;
+    case kPol1:       penaltyVal = RegularizationPol1(paramsVector); break;
+    case kCurvature:  penaltyVal = RegularizationTotalCurvature(paramsVector); break;
+    case kEntropy:    penaltyVal = RegularizationEntropy(paramsVector); break;
+    case kLog:        penaltyVal = RegularizationLog(paramsVector); break;
+  }
+
+  //if (penaltyVal > 1e10)
+  //  paramsVector2.Print();
+
+  penaltyVal *= fgRegularizationWeight;
+
+  // Ad
+  TVectorD measGuessVector(fgMaxInput);
+  measGuessVector = (*fgCorrelationMatrix) * paramsVector;
+
+  // Ad - m
+  measGuessVector -= (*fgCurrentESDVector);
+
+  //measGuessVector.Print();
+
+  TVectorD copy(measGuessVector);
+
+  // (Ad - m) W
+  // this step can be optimized because currently only the diagonal elements of fgCorrelationCovarianceMatrix are used
+  // normal way is like this:
+  // measGuessVector *= (*fgCorrelationCovarianceMatrix);
+  // optimized way like this:
+  for (Int_t i=0; i<fgMaxInput; ++i)
+    measGuessVector[i] *= (*fgCorrelationCovarianceMatrix)(i, i);
+
+  //measGuessVector.Print();
+
+  // (Ad - m) W (Ad - m)
+  // the factor 1e6 prevents that a very small number (measGuessVector[i]) is multiplied with a very
+  // big number ((*fgCorrelationCovarianceMatrix)(i, i)) (see ApplyMinuitFit)
+  Double_t chi2FromFit = measGuessVector * copy * 1e6;
+
+  chi2 = chi2FromFit + penaltyVal;
+
+  static Int_t callCount = 0;
+  if ((callCount++ % 10000) == 0)
+    Printf("AliUnfolding::Chi2Function: Iteration %d: %f %f --> %f", callCount, chi2FromFit, penaltyVal, chi2);
+}
+
+//____________________________________________________________________
+void AliUnfolding::DrawGuess(Double_t *params)
+{
+  //
+  // draws residuals of solution suggested by params and effect of regularization
+  //
+
+  // d
+  static TVectorD paramsVector(fgMaxParams);
+  for (Int_t i=0; i<fgMaxParams; ++i)
+    paramsVector[i] = params[i] * params[i];
+
+  // Ad
+  TVectorD measGuessVector(fgMaxInput);
+  measGuessVector = (*fgCorrelationMatrix) * paramsVector;
+
+  // Ad - m
+  measGuessVector -= (*fgCurrentESDVector);
+
+  TVectorD copy(measGuessVector);
+  //copy.Print();
+
+  // (Ad - m) W
+  // this step can be optimized because currently only the diagonal elements of fgCorrelationCovarianceMatrix are used
+  // normal way is like this:
+  // measGuessVector *= (*fgCorrelationCovarianceMatrix);
+  // optimized way like this:
+  for (Int_t i=0; i<fgMaxInput; ++i)
+    measGuessVector[i] *= (*fgCorrelationCovarianceMatrix)(i, i);
+
+  // (Ad - m) W (Ad - m)
+  // the factor 1e6 prevents that a very small number (measGuessVector[i]) is multiplied with a very
+  // big number ((*fgCorrelationCovarianceMatrix)(i, i)) (see ApplyMinuitFit)
+  //Double_t chi2FromFit = measGuessVector * copy * 1e6;
+
+  // draw residuals
+  TH1* residuals = new TH1F("residuals", "residuals", fgMaxInput, -0.5, fgMaxInput - 0.5);
+  for (Int_t i=0; i<fgMaxInput; ++i)
+    residuals->SetBinContent(i+1, measGuessVector(i) * copy(i) * 1e6);
+  new TCanvas; residuals->DrawCopy(); gPad->SetLogy();
+
+  // draw penalty
+  if (fgRegularizationType != kPol1) {
+    Printf("Drawing not supported");
+    return;
+  }
+
+  TH1* penalty = new TH1F("penalty", "penalty", fgMaxParams, -0.5, fgMaxParams - 0.5);
+
+  for (Int_t i=2+1; i<fgMaxParams; ++i)
+  {
+    if (params[i-1] == 0)
+      continue;
+
+    Double_t right  = params[i];
+    Double_t middle = params[i-1];
+    Double_t left   = params[i-2];
+
+    Double_t der1 = (right - middle);
+    Double_t der2 = (middle - left);
+
+    Double_t diff = (der1 - der2) / middle;
+
+    penalty->SetBinContent(i-1, diff * diff);
+    
+    //Printf("%d %f %f %f %f", i-1, left, middle, right, diff);
+  }
+  new TCanvas; penalty->DrawCopy(); gPad->SetLogy();
+}
+
+//____________________________________________________________________
+void AliUnfolding::TF1Function(Int_t& unused1, Double_t* unused2, Double_t& chi2, Double_t *params, Int_t unused3)
+{
+  //
+  // fit function for minuit
+  // uses the TF1 stored in fgFitFunction
+  //
+
+  for (Int_t i=0; i<fgFitFunction->GetNpar(); i++)
+    fgFitFunction->SetParameter(i, params[i]);
+
+  Double_t* params2 = new Double_t[fgMaxParams];
+
+  for (Int_t i=0; i<fgMaxParams; ++i)
+    params2[i] = fgFitFunction->Eval(i);
+
+  Chi2Function(unused1, unused2, chi2, params2, unused3);
+  
+  delete[] params2;
+
+  if (fgDebug)
+    Printf("%f", chi2);
+}
+
+//____________________________________________________________________
+Int_t AliUnfolding::UnfoldWithFunction(TH2* correlation, TH1* efficiency, TH1* measured, TH1* /* initialConditions */, TH1* aResult)
+{
+  //
+  // correct spectrum using minuit chi2 method applying a functional fit
+  //
+  
+  if (!fgFitFunction)
+  {
+    Printf("AliUnfolding::UnfoldWithFunction: ERROR fit function not set. Exiting.");
+    return -1;
+  }    
+  
+  SetChi2Regularization(kNone, 0);
+  
+  SetStaticVariables(correlation, measured);
+  
+  // Initialize TMinuit via generic fitter interface
+  TVirtualFitter *minuit = TVirtualFitter::Fitter(0, fgFitFunction->GetNpar());
+
+  minuit->SetFCN(TF1Function);
+  for (Int_t i=0; i<fgFitFunction->GetNpar(); i++)
+  {
+    Double_t lower, upper;
+    fgFitFunction->GetParLimits(i, lower, upper);
+    minuit->SetParameter(i, Form("param%d",i), fgFitFunction->GetParameter(i), fgMinuitStepSize, lower, upper);
+  }
+
+  Double_t arglist[100];
+  arglist[0] = 0;
+  minuit->ExecuteCommand("SET PRINT", arglist, 1);
+  minuit->ExecuteCommand("MIGRAD", arglist, 0);
+  //minuit->ExecuteCommand("MINOS", arglist, 0);
+
+  for (Int_t i=0; i<fgFitFunction->GetNpar(); i++)
+    fgFitFunction->SetParameter(i, minuit->GetParameter(i));
+
+  for (Int_t i=0; i<fgMaxParams; ++i)
+  {
+    Double_t value = fgFitFunction->Eval(i);
+    if (fgDebug)
+      Printf("%d : %f", i, value);
+    
+    if (efficiency)
+    {
+      if (efficiency->GetBinContent(i+1) > 0)
+      {
+        value /= efficiency->GetBinContent(i+1);
+      }
+      else
+        value = 0;
+    }
+    aResult->SetBinContent(i+1, value);
+    aResult->SetBinError(i+1, 0);
+  }
+  
+  return 0;
+}
+
+//____________________________________________________________________
+void AliUnfolding::CreateOverflowBin(TH2* correlation, TH1* measured)
+{
+  // Finds the first bin where the content is below fgStatLimit and combines all values for this bin and larger bins
+  // The same limit is applied to the correlation  
+  
+  Int_t lastBin = 0;
+  for (Int_t i=1; i<=measured->GetNbinsX(); ++i)
+  {
+    if (measured->GetBinContent(i) <= fgOverflowBinLimit)
+    {
+      lastBin = i;
+      break;
+    }
+  }
+  
+  if (lastBin == 0)
+  {
+    Printf("AliUnfolding::CreateOverflowBin: WARNING: First bin is already below limit of %f", fgOverflowBinLimit);
+    return;
+  }
+  
+  Printf("AliUnfolding::CreateOverflowBin: Bin limit in measured spectrum determined to be %d", lastBin);
+  
+  TCanvas* canvas = 0;
+
+  if (fgDebug)
+  {
+    canvas = new TCanvas("StatSolution", "StatSolution", 1000, 800);
+    canvas->Divide(2, 2);
+
+    canvas->cd(1);
+    measured->SetStats(kFALSE);
+    measured->DrawCopy();
+    gPad->SetLogy();
+
+    canvas->cd(2);
+    correlation->SetStats(kFALSE);
+    correlation->DrawCopy("COLZ");
+  }
+
+  measured->SetBinContent(lastBin, measured->Integral(lastBin, measured->GetNbinsX()));
+  for (Int_t i=lastBin+1; i<=measured->GetNbinsX(); ++i)
+  {
+    measured->SetBinContent(i, 0);
+    measured->SetBinError(i, 0);
+  }
+  // the error is set to sqrt(N), better solution possible?, sum of relative errors of all contributions???
+  measured->SetBinError(lastBin, TMath::Sqrt(measured->GetBinContent(lastBin)));
+
+  Printf("AliUnfolding::CreateOverflowBin: This bin has now %f +- %f entries", measured->GetBinContent(lastBin), measured->GetBinError(lastBin));
+
+  for (Int_t i=1; i<=correlation->GetNbinsX(); ++i)
+  {
+    correlation->SetBinContent(i, lastBin, correlation->Integral(i, i, lastBin, correlation->GetNbinsY()));
+    // the error is set to sqrt(N), better solution possible?, sum of relative errors of all contributions???
+    correlation->SetBinError(i, lastBin, TMath::Sqrt(correlation->GetBinContent(i, lastBin)));
+
+    for (Int_t j=lastBin+1; j<=correlation->GetNbinsY(); ++j)
+    {
+      correlation->SetBinContent(i, j, 0);
+      correlation->SetBinError(i, j, 0);
+    }
+  }
+
+  Printf("AliUnfolding::CreateOverflowBin: Adjusted correlation matrix!");
+
+  if (canvas)
+  {
+    canvas->cd(3);
+    measured->DrawCopy();
+    gPad->SetLogy();
+
+    canvas->cd(4);
+    correlation->DrawCopy("COLZ");
+  }
+}
index f07e27b..d11de22 100644 (file)
@@ -5,7 +5,8 @@
 
 //
 // class that implements several unfolding methods
-// E.g. chi2 minimization and bayesian unfolding
+// I.e. chi2 minimization and bayesian unfolding
+// The whole class is static and not thread-safe (due to the fact that MINUIT unfolding is not thread-safe)
 //
 
 // TMatrixD, TVectorD defined here, because it does not seem possible to predeclare these (or i do not know how)
 
 class TH1;
 class TH2;
-class TH1F;
-class TH2F;
-class TH3F;
 class TF1;
-class TCollection;
 
 class AliUnfolding : public TObject
 {
   public:
     enum RegularizationType { kNone = 0, kPol0, kPol1, kLog, kEntropy, kCurvature };
-    enum MethodType { kChi2Minimization = 0, kBayesian = 1 };
+    enum MethodType { kInvalid = -1, kChi2Minimization = 0, kBayesian = 1, kFunction = 2};
 
-    AliUnfolding();
-    virtual ~AliUnfolding();
+    virtual ~AliUnfolding() {};
 
-    void SetInput(TH2* correlationMatrix, TH1* efficiency, TH1* measured) { fCurrentCorrelation = correlationMatrix; fCurrentEfficiency = efficiency; fCurrentESD = measured; }
-    void SetInitialConditions(TH1* initialConditions) { fInitialConditions = initialConditions; }
-    const TH1* GetResult() const { return fResult; }
+    static void SetUnfoldingMethod(MethodType methodType);
+    static void SetCreateOverflowBin(Float_t overflowBinLimit);
+    static void SetSkipBinsBegin(Int_t nBins);
+    static void SetNbins(Int_t nMeasured, Int_t nUnfolded);
+    static void SetChi2Regularization(RegularizationType type, Float_t weight);
+    static void SetMinuitStepSize(Float_t stepSize) { fgMinuitStepSize = stepSize; }
+    static void SetNormalizeInput(Bool_t flag) { fgNormalizeInput = flag; }
+    static void SetBayesianParameters(Float_t smoothing, Int_t nIterations);
+    static void SetFunction(TF1* function);
+    static void SetDebug(Bool_t flag) { fgDebug = flag; }
 
-    static void SetParameters(Int_t measuredBins, Int_t unfoldedBins, Bool_t bigbin) { fMaxInput = measuredBins; fMaxParams = unfoldedBins; fgCreateBigBin = bigbin; }
-    static void SetChi2MinimizationParameters(RegularizationType type, Float_t weight) { fgRegularizationType = type; fgRegularizationWeight = weight; }
-    static void SetRegularizationRange(Int_t start, Int_t end) { fgRegularizationRangeStart = start; fgRegularizationRangeEnd = end; }
-    static void SetBayesianParameters(Float_t smoothing, Int_t nIterations) { fgBayesianSmoothing = smoothing; fgBayesianIterations = nIterations; }
+    static Int_t Unfold(TH2* correlation, TH1* efficiency, TH1* measured, TH1* initialConditions, TH1* result, Bool_t check = kFALSE);
 
-    Int_t ApplyMinuitFit(Bool_t check = kFALSE);
-    Int_t ApplyBayesianMethod(Bool_t determineError = kTRUE);
-    Int_t ApplyNBDFit();
-    Int_t ApplyLaszloMethod();
+  protected:
+    AliUnfolding() {};
 
-    TH1* StatisticalUncertainty(MethodType methodType, Bool_t randomizeMeasured, Bool_t randomizeResponse, TH1* compareTo = 0);
+    static Int_t UnfoldWithMinuit(TH2* correlation, TH1* efficiency, TH1* measured, TH1* initialConditions, TH1* result, Bool_t check);
+    static Int_t UnfoldWithBayesian(TH2* correlation, TH1* aEfficiency, TH1* measured, TH1* initialConditions, TH1* aResult);
+    static Int_t UnfoldWithFunction(TH2* correlation, TH1* efficiency, TH1* measured, TH1* /* initialConditions */, TH1* aResult);
+
+    static void DrawGuess(Double_t *params);
+    static void CreateOverflowBin(TH2* correlation, TH1* measured); 
+    static void SetStaticVariables(TH2* correlation, TH1* measured);
 
-  protected:
     static Double_t RegularizationPol0(TVectorD& params);
     static Double_t RegularizationPol1(TVectorD& params);
     static Double_t RegularizationTotalCurvature(TVectorD& params);
     static Double_t RegularizationEntropy(TVectorD& params);
     static Double_t RegularizationLog(TVectorD& params);
 
-    static void MinuitFitFunction(Int_t&, Double_t*, Double_t& chi2, Double_t *params, Int_t);
-    static void MinuitNBD(Int_t& unused1, Double_t* unused2, Double_t& chi2, Double_t *params, Int_t unused3);
-
-    void SetupCurrentHists();
-
-    Int_t UnfoldWithBayesian(* aEfficiency, TH1* measured, TH1* initialConditions, TH1* aResult, Float_t regPar, Int_t nIterations);
-    Int_t UnfoldWithMinuit(TH1* correlation, TH1* aEfficiency, TH1* measured, TH1* initialConditions, TH1* result, Bool_t check);
-
-    Float_t BayesCovarianceDerivate(Float_t matrixM[251][251], TH2* hResponse, Int_t k, Int_t i, Int_t r, Int_t u);
-
-    TH1* fCurrentESD;         //! measured spectrum
-    TH2* fCurrentCorrelation; //! correlation matrix
-    TH1* fCurrentEfficiency;  //! efficiency
-    TH1* fInitialConditions;  //! initial conditions
-    TH1* fResult;             //! unfolding result
+    static void Chi2Function(Int_t&, Double_t*, Double_t& chi2, Double_t *params, Int_t);
+    static void TF1Function(Int_t& unused1, Double_t* unused2, Double_t& chi2, Double_t *params, Int_t unused3);
 
     // static variable to be accessed by MINUIT
-    static TMatrixD* fgCorrelationMatrix;            //! contains fCurrentCorrelation in matrix form
-    static TMatrixD* fgCorrelationCovarianceMatrix;  //! contains the errors of fCurrentESD
-    static TVectorD* fgCurrentESDVector;             //! contains fCurrentESD
-    static TVectorD* fgEntropyAPriori;               //! a-priori distribution for entropy regularization
+    static TMatrixD* fgCorrelationMatrix;            // contains fCurrentCorrelation in matrix form
+    static TMatrixD* fgCorrelationCovarianceMatrix;  // contains the errors of fCurrentESD
+    static TVectorD* fgCurrentESDVector;             // contains fCurrentESD
+    static TVectorD* fgEntropyAPriori;               // a-priori distribution for entropy regularization
+
+    static TF1* fgFitFunction;                       // fit function
 
-    static TF1* fgNBD;   //! negative binomial distribution
+    // --- configuration params follow ---
+    static MethodType fgMethodType;                  // unfolding method to be used
+    static Int_t fgMaxParams;                        // bins in unfolded histogram = number of fit params
+    static Int_t fgMaxInput;                         // bins in measured histogram
+    static Float_t fgOverflowBinLimit;               // to fix fluctuations at high multiplicities, all entries above the limit are summarized in one bin
 
-    static Int_t fgMaxParams;  //! bins in unfolded histogram = number of fit params
-    static Int_t fgMaxInput;   //! bins in measured histogram
+    static RegularizationType fgRegularizationType;  // regularization that is used during Chi2 method
+    static Float_t fgRegularizationWeight;           // factor for regularization term
+    static Int_t fgSkipBinsBegin;                    // (optional) skip the given number of bins in the regularization
+    static Float_t fgMinuitStepSize;                 // (usually not needed to be changed) step size in minimization
+    static Bool_t fgNormalizeInput;                  // normalize input spectrum
 
-    // configuration params follow
-    static RegularizationType fgRegularizationType; //! regularization that is used during Chi2 method
-    static Float_t fgRegularizationWeight;          //! factor for regularization term
-    static Int_t fgRegularizationRangeStart;        //! first bin where regularization is applied
-    static Int_t fgRegularizationRangeEnd;          //! last bin + 1 where regularization is applied
-    static Bool_t  fgCreateBigBin;                  //! to fix fluctuations at high multiplicities, all entries above a certain limit are summarized in one bin
+    static Float_t fgBayesianSmoothing;              // smoothing parameter (0 = no smoothing)
+    static Int_t   fgBayesianIterations;             // number of iterations in Bayesian method
 
-    static Float_t fgBayesianSmoothing;             //! smoothing parameter (0 = no smoothing)
-    static Int_t   fgBayesianIterations;            //! number of iterations in Bayesian method
-    // end of configuration
+    static Bool_t fgDebug;                           // debug flag
+    // --- end of configuration ---
 
- private:
+private:
     AliUnfolding(const AliUnfolding&);
     AliUnfolding& operator=(const AliUnfolding&);
 
index ef2767c..9a61f4b 100644 (file)
@@ -18,6 +18,7 @@
 #pragma link C++ class AliPWG0Helper+;
 
 #pragma link C++ class AliCorrection+;
+#pragma link C++ class AliUnfolding+;
 
 #pragma link C++ class AliMultiplicityCorrection+;
 
index a23aedb..7a17ec1 100644 (file)
@@ -9,6 +9,7 @@ HDRS = dNdEta/dNdEtaAnalysis.h \
       AliCorrectionMatrix3D.h \
       AliCorrection.h \
       AliPWG0Helper.h \
+      AliUnfolding.h \
       multiplicity/AliMultiplicityCorrection.h
 
 SRCS = $(HDRS:.h=.cxx)
index cc6f34a..45d4217 100644 (file)
@@ -29,7 +29,6 @@
 #include <TH2F.h>
 #include <TH3F.h>
 #include <TDirectory.h>
-#include <TVirtualFitter.h>
 #include <TCanvas.h>
 #include <TString.h>
 #include <TF1.h>
 
 ClassImp(AliMultiplicityCorrection)
 
-const Int_t AliMultiplicityCorrection::fgkMaxInput  = 120;  // bins in measured histogram
-const Int_t AliMultiplicityCorrection::fgkMaxParams = 120;  // bins in unfolded histogram = number of fit params
-
-TMatrixD* AliMultiplicityCorrection::fgCorrelationMatrix = 0;
-TMatrixD* AliMultiplicityCorrection::fgCorrelationCovarianceMatrix = 0;
-TVectorD* AliMultiplicityCorrection::fgCurrentESDVector = 0;
-TVectorD* AliMultiplicityCorrection::fgEntropyAPriori = 0;
-
-Bool_t AliMultiplicityCorrection::fgCreateBigBin = kTRUE;
-AliMultiplicityCorrection::RegularizationType AliMultiplicityCorrection::fgRegularizationType = AliMultiplicityCorrection::kPol1;
-Float_t AliMultiplicityCorrection::fgRegularizationWeight = 10000;
-Int_t AliMultiplicityCorrection::fgNParamsMinuit = AliMultiplicityCorrection::fgkMaxParams;
-
-TF1* AliMultiplicityCorrection::fgNBD = 0;
-
-Float_t AliMultiplicityCorrection::fgBayesianSmoothing  = 1;           // smoothing parameter (0 = no smoothing)
-Int_t   AliMultiplicityCorrection::fgBayesianIterations = 10;         // number of iterations in Bayesian method
-
 // These are the areas where the quality of the unfolding results are evaluated
 // Default defined here, call SetQualityRegions to change them
 // unit is in multiplicity (not in bin!)
@@ -221,6 +202,10 @@ AliMultiplicityCorrection::AliMultiplicityCorrection(const Char_t* name, const C
   }
 
   TH1::AddDirectory(oldStatus);
+
+  AliUnfolding::SetNbins(120, 120);
+  AliUnfolding::SetSkipBinsBegin(1);
+  AliUnfolding::SetNormalizeInput(kTRUE);
 }
 
 //____________________________________________________________________
@@ -542,330 +527,13 @@ void AliMultiplicityCorrection::FillCorrection(Float_t vtx, Int_t generated05, I
 }
 
 //____________________________________________________________________
-Double_t AliMultiplicityCorrection::RegularizationPol0(TVectorD& params)
-{
-  // homogenity term for minuit fitting
-  // pure function of the parameters
-  // prefers constant function (pol0)
-
-  Double_t chi2 = 0;
-
-  // ignore the first bin here. on purpose...
-  for (Int_t i=2; i<fgNParamsMinuit; ++i)
-  {
-    Double_t right  = params[i];
-    Double_t left   = params[i-1];
-
-    if (right != 0)
-    {
-      Double_t diff = 1 - left / right;
-      chi2 += diff * diff;
-    }
-  }
-
-  return chi2 / 100.0;
-}
-
-//____________________________________________________________________
-Double_t AliMultiplicityCorrection::RegularizationPol1(TVectorD& params)
-{
-  // homogenity term for minuit fitting
-  // pure function of the parameters
-  // prefers linear function (pol1)
-
-  Double_t chi2 = 0;
-
-  // ignore the first bin here. on purpose...
-  for (Int_t i=2+1; i<fgNParamsMinuit; ++i)
-  {
-    if (params[i-1] == 0)
-      continue;
-
-    Double_t right  = params[i];
-    Double_t middle = params[i-1];
-    Double_t left   = params[i-2];
-
-    Double_t der1 = (right - middle);
-    Double_t der2 = (middle - left);
-
-    Double_t diff = (der1 - der2) / middle;
-
-    chi2 += diff * diff;
-  }
-
-  return chi2;
-}
-
-//____________________________________________________________________
-Double_t AliMultiplicityCorrection::RegularizationLog(TVectorD& params)
-{
-  // homogenity term for minuit fitting
-  // pure function of the parameters
-  // prefers linear function (pol1)
-
-  Double_t chi2 = 0;
-
-  /*Float_t der[fgkMaxParams];
-
-  for (Int_t i=0; i<fgkMaxParams-1; ++i)
-    der[i] = params[i+1] - params[i];
-
-  for (Int_t i=0; i<fgkMaxParams-6; ++i)
-  {
-    for (Int_t j=1; j<=5; ++j)
-    {
-      Double_t diff = der[i+j] - der[i];
-      chi2 += diff * diff;
-    }
-  }*/
-
-  // ignore the first bin here. on purpose...
-  for (Int_t i=2+1; i<fgNParamsMinuit; ++i)
-  {
-    if (params[i-1] == 0)
-      continue;
-
-    Double_t right  = log(params[i]);
-    Double_t middle = log(params[i-1]);
-    Double_t left   = log(params[i-2]);
-
-    Double_t der1 = (right - middle);
-    Double_t der2 = (middle - left);
-
-    Double_t diff = (der1 - der2) / middle;
-
-    chi2 += diff * diff;
-  }
-
-  return chi2 * 100;
-}
-
-//____________________________________________________________________
-Double_t AliMultiplicityCorrection::RegularizationTotalCurvature(TVectorD& params)
-{
-  // homogenity term for minuit fitting
-  // pure function of the parameters
-  // minimizes the total curvature (from Unfolding Methods In High-Energy Physics Experiments,
-  // V. Blobel (Hamburg U.) . DESY 84/118, Dec 1984. 40pp.
-
-  Double_t chi2 = 0;
-
-  // ignore the first bin here. on purpose...
-  for (Int_t i=3; i<fgNParamsMinuit; ++i)
-  {
-    Double_t right  = params[i];
-    Double_t middle = params[i-1];
-    Double_t left   = params[i-2];
-
-    Double_t der1 = (right - middle);
-    Double_t der2 = (middle - left);
-
-    Double_t diff = (der1 - der2);
-
-    chi2 += diff * diff;
-  }
-
-  return chi2 * 1e4;
-}
-
-//____________________________________________________________________
-Double_t AliMultiplicityCorrection::RegularizationEntropy(TVectorD& params)
-{
-  // homogenity term for minuit fitting
-  // pure function of the parameters
-  // calculates entropy, from
-  // The method of reduced cross-entropy (M. Schmelling 1993)
-
-  Double_t paramSum = 0;
-  // ignore the first bin here. on purpose...
-  for (Int_t i=1; i<fgNParamsMinuit; ++i)
-    paramSum += params[i];
-
-  Double_t chi2 = 0;
-  for (Int_t i=1; i<fgNParamsMinuit; ++i)
-  {
-    //Double_t tmp = params[i] / paramSum;
-    Double_t tmp = params[i];
-    if (tmp > 0 && (*fgEntropyAPriori)[i] > 0)
-    {
-      chi2 += tmp * TMath::Log(tmp / (*fgEntropyAPriori)[i]);
-    }
-  }
-
-  return 100.0 + chi2;
-}
-
-//____________________________________________________________________
-void AliMultiplicityCorrection::MinuitNBD(Int_t& unused1, Double_t* unused2, Double_t& chi2, Double_t *params, Int_t unused3)
-{
-  //
-  // fit function for minuit
-  // does: nbd
-  // func = new TF1("nbd", "[0] * TMath::Binomial([2]+TMath::Nint(x)-1, [2]-1) * pow([1] / ([1]+[2]), TMath::Nint(x)) * pow(1 + [1]/[2], -[2])", 0, 50);
-  // func->SetParNames("scaling", "averagen", "k");
-  // func->SetParLimits(0, 0.001, fCurrentESD->GetMaximum() * 1000);
-  // func->SetParLimits(1, 0.001, 1000);
-  // func->SetParLimits(2, 0.001, 1000);
-  //
-
-  fgNBD->SetParameters(params[0], params[1], params[2]);
-
-  Double_t params2[fgkMaxParams];
-
-  for (Int_t i=0; i<fgkMaxParams; ++i)
-    params2[i] = fgNBD->Eval(i);
-
-  MinuitFitFunction(unused1, unused2, chi2, params2, unused3);
-
-  printf("%f %f %f --> %f\n", params[0], params[1], params[2], chi2);
-}
-
-//____________________________________________________________________
-void AliMultiplicityCorrection::MinuitFitFunction(Int_t&, Double_t*, Double_t& chi2, Double_t *params, Int_t)
-{
-  //
-  // fit function for minuit
-  // does: (m - Ad)W(m - Ad) where m = measured, A correlation matrix, d = guess, W = covariance matrix
-  //
-
-  // d
-  static TVectorD paramsVector(fgNParamsMinuit);
-  for (Int_t i=0; i<fgNParamsMinuit; ++i)
-    paramsVector[i] = params[i] * params[i];
-
-  // calculate penalty factor
-  Double_t penaltyVal = 0;
-  switch (fgRegularizationType)
-  {
-    case kNone:       break;
-    case kPol0:       penaltyVal = RegularizationPol0(paramsVector); break;
-    case kPol1:       penaltyVal = RegularizationPol1(paramsVector); break;
-    case kCurvature:  penaltyVal = RegularizationTotalCurvature(paramsVector); break;
-    case kEntropy:    penaltyVal = RegularizationEntropy(paramsVector); break;
-    case kLog:        penaltyVal = RegularizationLog(paramsVector); break;
-  }
-
-  //if (penaltyVal > 1e10)
-  //  paramsVector2.Print();
-
-  penaltyVal *= fgRegularizationWeight;
-
-  // Ad
-  TVectorD measGuessVector(fgkMaxInput);
-  measGuessVector = (*fgCorrelationMatrix) * paramsVector;
-
-  // Ad - m
-  measGuessVector -= (*fgCurrentESDVector);
-
-  //measGuessVector.Print();
-
-  TVectorD copy(measGuessVector);
-
-  // (Ad - m) W
-  // this step can be optimized because currently only the diagonal elements of fgCorrelationCovarianceMatrix are used
-  // normal way is like this:
-  // measGuessVector *= (*fgCorrelationCovarianceMatrix);
-  // optimized way like this:
-  for (Int_t i=0; i<fgkMaxInput; ++i)
-    measGuessVector[i] *= (*fgCorrelationCovarianceMatrix)(i, i);
-
-  //measGuessVector.Print();
-
-  // reduce influence of first bin
-  //copy(0) *= 0.01;
-
-  // (Ad - m) W (Ad - m)
-  // the factor 1e6 prevents that a very small number (measGuessVector[i]) is multiplied with a very
-  // big number ((*fgCorrelationCovarianceMatrix)(i, i)) (see ApplyMinuitFit)
-  Double_t chi2FromFit = measGuessVector * copy * 1e6;
-
-  chi2 = chi2FromFit + penaltyVal;
-
-  static Int_t callCount = 0;
-  if ((callCount++ % 10000) == 0)
-    printf("%d) %f %f --> %f\n", callCount, chi2FromFit, penaltyVal, chi2);
-}
-
-//____________________________________________________________________
-void AliMultiplicityCorrection::DrawGuess(Double_t *params)
-{
-  //
-  // draws residuals of solution suggested by params and effect of regularization
-  //
-
-  // d
-  static TVectorD paramsVector(fgNParamsMinuit);
-  for (Int_t i=0; i<fgNParamsMinuit; ++i)
-    paramsVector[i] = params[i] * params[i];
-
-  // Ad
-  TVectorD measGuessVector(fgkMaxInput);
-  measGuessVector = (*fgCorrelationMatrix) * paramsVector;
-
-  // Ad - m
-  measGuessVector -= (*fgCurrentESDVector);
-
-  TVectorD copy(measGuessVector);
-  //copy.Print();
-
-  // (Ad - m) W
-  // this step can be optimized because currently only the diagonal elements of fgCorrelationCovarianceMatrix are used
-  // normal way is like this:
-  // measGuessVector *= (*fgCorrelationCovarianceMatrix);
-  // optimized way like this:
-  for (Int_t i=0; i<fgkMaxInput; ++i)
-    measGuessVector[i] *= (*fgCorrelationCovarianceMatrix)(i, i);
-
-  // (Ad - m) W (Ad - m)
-  // the factor 1e6 prevents that a very small number (measGuessVector[i]) is multiplied with a very
-  // big number ((*fgCorrelationCovarianceMatrix)(i, i)) (see ApplyMinuitFit)
-  //Double_t chi2FromFit = measGuessVector * copy * 1e6;
-
-  // draw residuals
-  TH1* residuals = new TH1F("residuals", "residuals", fgkMaxInput, -0.5, fgkMaxInput - 0.5);
-  for (Int_t i=0; i<fgkMaxInput; ++i)
-    residuals->SetBinContent(i+1, measGuessVector(i) * copy(i) * 1e6);
-  new TCanvas; residuals->DrawCopy(); gPad->SetLogy();
-
-  // draw penalty
-  if (fgRegularizationType != kPol1) {
-    Printf("Drawing not supported");
-    return;
-  }
-
-  TH1* penalty = new TH1F("penalty", "penalty", fgNParamsMinuit, -0.5, fgNParamsMinuit - 0.5);
-
-  for (Int_t i=2+1; i<fgNParamsMinuit; ++i)
-  {
-    if (params[i-1] == 0)
-      continue;
-
-    Double_t right  = params[i];
-    Double_t middle = params[i-1];
-    Double_t left   = params[i-2];
-
-    Double_t der1 = (right - middle);
-    Double_t der2 = (middle - left);
-
-    Double_t diff = (der1 - der2) / middle;
-
-    penalty->SetBinContent(i-1, diff * diff);
-    
-    //Printf("%d %f %f %f %f", i-1, left, middle, right, diff);
-  }
-  new TCanvas; penalty->DrawCopy(); gPad->SetLogy();
-}
-
-//____________________________________________________________________
-void AliMultiplicityCorrection::SetupCurrentHists(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t createBigBin)
+void AliMultiplicityCorrection::SetupCurrentHists(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType)
 {
   //
   // fills fCurrentESD, fCurrentCorrelation
   // resets fMultiplicityESDCorrected
   //
 
-  Bool_t display = kFALSE;
-
   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
 
   fMultiplicityESDCorrected[correlationID]->Reset();
@@ -886,88 +554,11 @@ void AliMultiplicityCorrection::SetupCurrentHists(Int_t inputRange, Bool_t fullP
     }
   }
 
-  fCurrentCorrelation = hist->Project3D("zy");
+  fCurrentCorrelation = (TH2*) hist->Project3D("zy");
   fCurrentCorrelation->Sumw2();
   
   Printf("AliMultiplicityCorrection::SetupCurrentHists: Statistics information: %.f entries in correlation map; %.f entries in measured spectrum", fCurrentCorrelation->Integral(), fCurrentESD->Integral());
 
-  fLastBinLimit = 0;
-  for (Int_t i=1; i<=fCurrentESD->GetNbinsX(); ++i)
-  {
-    if (fCurrentESD->GetBinContent(i) <= 5)
-    {
-      fLastBinLimit = i;
-      break;
-    }
-  }
-  
-  Printf("AliMultiplicityCorrection::SetupCurrentHists: Bin limit in measured spectrum determined to be %d", fLastBinLimit);
-
-  if (createBigBin)
-  {
-    if (fLastBinLimit > 0)
-    {
-      TCanvas* canvas = 0;
-
-      if (display)
-      {
-        canvas = new TCanvas("StatSolution", "StatSolution", 1000, 800);
-        canvas->Divide(2, 2);
-
-        canvas->cd(1);
-        fCurrentESD->GetXaxis()->SetRangeUser(0, 200);
-        fCurrentESD->SetStats(kFALSE);
-        fCurrentESD->SetTitle(";measured multiplicity");
-        fCurrentESD->DrawCopy();
-        gPad->SetLogy();
-
-        canvas->cd(2);
-        fCurrentCorrelation->GetXaxis()->SetRangeUser(0, 250);
-        fCurrentCorrelation->GetYaxis()->SetRangeUser(0, 250);
-        fCurrentCorrelation->SetTitle(";true multiplicity;measured multiplicity");
-        fCurrentCorrelation->SetStats(kFALSE);
-        fCurrentCorrelation->DrawCopy("COLZ");
-      }
-
-      printf("Bin limit in measured spectrum is %d.\n", fLastBinLimit);
-      fCurrentESD->SetBinContent(fLastBinLimit, fCurrentESD->Integral(fLastBinLimit, fCurrentESD->GetNbinsX()));
-      for (Int_t i=fLastBinLimit+1; i<=fCurrentESD->GetNbinsX(); ++i)
-      {
-        fCurrentESD->SetBinContent(i, 0);
-        fCurrentESD->SetBinError(i, 0);
-      }
-      // the error is set to sqrt(N), better solution possible?, sum of relative errors of all contributions???
-      fCurrentESD->SetBinError(fLastBinLimit, TMath::Sqrt(fCurrentESD->GetBinContent(fLastBinLimit)));
-
-      printf("This bin has now %f +- %f entries\n", fCurrentESD->GetBinContent(fLastBinLimit), fCurrentESD->GetBinError(fLastBinLimit));
-
-      for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
-      {
-        fCurrentCorrelation->SetBinContent(i, fLastBinLimit, fCurrentCorrelation->Integral(i, i, fLastBinLimit, fCurrentCorrelation->GetNbinsY()));
-        // the error is set to sqrt(N), better solution possible?, sum of relative errors of all contributions???
-        fCurrentCorrelation->SetBinError(i, fLastBinLimit, TMath::Sqrt(fCurrentCorrelation->GetBinContent(i, fLastBinLimit)));
-
-        for (Int_t j=fLastBinLimit+1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
-        {
-          fCurrentCorrelation->SetBinContent(i, j, 0);
-          fCurrentCorrelation->SetBinError(i, j, 0);
-        }
-      }
-
-      printf("Adjusted correlation matrix!\n");
-
-      if (canvas)
-      {
-        canvas->cd(3);
-        fCurrentESD->DrawCopy();
-        gPad->SetLogy();
-
-        canvas->cd(4);
-        fCurrentCorrelation->DrawCopy("COLZ");
-      }
-    }
-  }
-
 #if 0 // does not help
   // null bins with one entry
   Int_t nNulledBins = 0;
@@ -1032,294 +623,58 @@ TH1* AliMultiplicityCorrection::GetTriggerEfficiency(Int_t inputRange)
 }
 
 //____________________________________________________________________
-void AliMultiplicityCorrection::SetRegularizationParameters(RegularizationType type, Float_t weight, Int_t minuitParams)
-{
-  //
-  // sets the parameters for chi2 minimization
-  //
-
-  fgRegularizationType = type;
-  fgRegularizationWeight = weight;
-
-  printf("AliMultiplicityCorrection::SetRegularizationParameters --> Regularization set to %d with weight %f\n", (Int_t) type, weight);
-
-  if (minuitParams > 0)
-  {
-    fgNParamsMinuit = minuitParams;
-    printf("AliMultiplicityCorrection::SetRegularizationParameters --> Number of MINUIT iterations set to %d\n", minuitParams);
-  }
-}
-
-//____________________________________________________________________
-void AliMultiplicityCorrection::SetBayesianParameters(Float_t smoothing, Int_t nIterations)
-{
-  //
-  // sets the parameters for Bayesian unfolding
-  //
-
-  fgBayesianSmoothing = smoothing;
-  fgBayesianIterations = nIterations;
-
-  printf("AliMultiplicityCorrection::SetBayesianParameters --> Paramaters set to %d iterations with smoothing %f\n", fgBayesianIterations, fgBayesianSmoothing);
-}
-
-//____________________________________________________________________
-Int_t AliMultiplicityCorrection::UnfoldWithMinuit(TH1* correlation, TH1* aEfficiency, TH1* measured, TH1* initialConditions, TH1* result, Bool_t check)
+Int_t AliMultiplicityCorrection::ApplyMinuitFit(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t check, TH1* initialConditions)
 {
   //
-  // implementation of unfolding (internal function)
-  //
-  // unfolds <measured> using response from <correlation> and effiency <aEfficiency>
-  // output is in <result>
-  // <initialConditions> set the initial values for the minimization, if 0 <measured> is used
-  // if <check> is true no unfolding is made, instead only the chi2 without unfolding is printed
+  // correct spectrum using minuit chi2 method
   //
-  // returns minuit status (0 = success), (-1 when check was set)
+  // for description of parameters, see AliUnfolding::Unfold
   //
 
-  //normalize measured
-  measured->Scale(1.0 / measured->Integral());
-
-  if (!fgCorrelationMatrix)
-    fgCorrelationMatrix = new TMatrixD(fgkMaxInput, fgNParamsMinuit);
-  if (!fgCorrelationCovarianceMatrix)
-    fgCorrelationCovarianceMatrix = new TMatrixD(fgkMaxInput, fgkMaxInput);
-  if (!fgCurrentESDVector)
-    fgCurrentESDVector = new TVectorD(fgkMaxInput);
-  if (!fgEntropyAPriori)
-    fgEntropyAPriori = new TVectorD(fgNParamsMinuit);
-
-  fgCorrelationMatrix->Zero();
-  fgCorrelationCovarianceMatrix->Zero();
-  fgCurrentESDVector->Zero();
-  fgEntropyAPriori->Zero();
-
-  // normalize correction for given nPart
-  for (Int_t i=1; i<=correlation->GetNbinsX(); ++i)
+  Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
+  
+  AliUnfolding::SetCreateOverflowBin(5);
+  AliUnfolding::SetUnfoldingMethod(AliUnfolding::kChi2Minimization);
+  SetupCurrentHists(inputRange, fullPhaseSpace, eventType);
+  
+  if (!initialConditions)
   {
-    Double_t sum = correlation->Integral(i, i, 1, correlation->GetNbinsY());
-    if (sum <= 0)
-      continue;
-    Float_t maxValue = 0;
-    Int_t maxBin = -1;
-    for (Int_t j=1; j<=correlation->GetNbinsY(); ++j)
-    {
-      // find most probably value
-      if (maxValue < correlation->GetBinContent(i, j))
-      {
-        maxValue = correlation->GetBinContent(i, j);
-        maxBin = j;
-      }
-
-      // npart sum to 1
-      correlation->SetBinContent(i, j, correlation->GetBinContent(i, j) / sum);
-      correlation->SetBinError(i, j, correlation->GetBinError(i, j) / sum);
-
-      if (i <= fgNParamsMinuit && j <= fgkMaxInput)
-        (*fgCorrelationMatrix)(j-1, i-1) = correlation->GetBinContent(i, j);
-    }
-
-    //printf("MPV for Ntrue = %f is %f\n", fCurrentCorrelation->GetXaxis()->GetBinCenter(i), fCurrentCorrelation->GetYaxis()->GetBinCenter(maxBin));
-  }
-
-  // Initialize TMinuit via generic fitter interface
-  TVirtualFitter *minuit = TVirtualFitter::Fitter(0, fgNParamsMinuit);
-  Double_t arglist[100];
-
-  // disable any output (-1), unfortuantly we do not see warnings anymore then. Have to find another way...
-  arglist[0] = 0;
-  minuit->ExecuteCommand("SET PRINT", arglist, 1);
-
-  // however, enable warnings
-  //minuit->ExecuteCommand("SET WAR", arglist, 0);
-
-  // set minimization function
-  minuit->SetFCN(MinuitFitFunction);
-
-  for (Int_t i=0; i<fgNParamsMinuit; i++)
-    (*fgEntropyAPriori)[i] = 1;
-
-  if (!initialConditions) {
-    initialConditions = measured;
-  } else {
-    printf("Using different starting conditions...\n");
-    //new TCanvas; initialConditions->DrawCopy();
+    initialConditions = (TH1*) fCurrentESD->Clone("initialConditions");
     initialConditions->Scale(1.0 / initialConditions->Integral());
-  }
-
-  for (Int_t i=0; i<fgNParamsMinuit; i++)
-    if (initialConditions->GetBinContent(i+1) > 0)
-      (*fgEntropyAPriori)[i] = initialConditions->GetBinContent(i+1);
-
-  Double_t* results = new Double_t[fgNParamsMinuit+1];
-  for (Int_t i=0; i<fgNParamsMinuit; ++i)
-  {
-    results[i] = initialConditions->GetBinContent(i+1);
-
-    // minimum value
     if (!check)
-      results[i] = TMath::Max(results[i], 1e-3);
-
-    Float_t stepSize = 0.1;
-
-    // minuit sees squared values to prevent it from going negative...
-    results[i] = TMath::Sqrt(results[i]);
-    //stepSize /= results[i]; // keep relative step size
-
-    minuit->SetParameter(i, Form("param%d", i), results[i], stepSize, 0, 0);
-  }
-  //minuit->SetParameter(fgkMaxParams, "alpha", 1e4, 1, 0, 0);
-  // bin 0 is filled with value from bin 1 (otherwise it's 0)
-  //minuit->SetParameter(0, "param0", (results[1] > minStartValue) ? results[1] : minStartValue, 0.1, 0, 0);
-  //results[0] = 0;
-  //minuit->SetParameter(0, "param0", 0, 0, 0, 0);
-
-  for (Int_t i=0; i<fgkMaxInput; ++i)
-  {
-    (*fgCurrentESDVector)[i] = measured->GetBinContent(i+1);
-    if (measured->GetBinError(i+1) > 0)
-      (*fgCorrelationCovarianceMatrix)(i, i) = (Double_t) 1e-6 / measured->GetBinError(i+1) / measured->GetBinError(i+1);
-
-    if ((*fgCorrelationCovarianceMatrix)(i, i) > 1e7)
-      (*fgCorrelationCovarianceMatrix)(i, i) = 0;
-    //Printf("%d, %e", i, (*fgCorrelationCovarianceMatrix)(i, i));
-  }
-
-  Int_t dummy = 0;
-  Double_t chi2 = 0;
-  MinuitFitFunction(dummy, 0, chi2, results, 0);
-  printf("Chi2 of initial parameters is = %f\n", chi2);
-
-  if (check)
-  {
-    DrawGuess(results);
-    delete[] results;
-    return -1;
-  }
-
-  // first param is number of iterations, second is precision....
-  arglist[0] = 1e6;
-  //arglist[1] = 1e-5;
-  //minuit->ExecuteCommand("SCAN", arglist, 0);
-  Int_t status = minuit->ExecuteCommand("MIGRAD", arglist, 1);
-  printf("MINUIT status is %d\n", status);
-  //minuit->ExecuteCommand("MIGRAD", arglist, 1);
-  //minuit->ExecuteCommand("MIGRAD", arglist, 1);
-  //printf("!!!!!!!!!!!!!! MIGRAD finished: Starting MINOS !!!!!!!!!!!!!!");
-  //minuit->ExecuteCommand("MINOS", arglist, 0);
-
-  //new TCanvas; aEfficiency->Draw();
-
-  for (Int_t i=0; i<fgNParamsMinuit; ++i)
-  {
-    results[i] = minuit->GetParameter(i);
-    if (aEfficiency->GetBinContent(i+1) > 0)
-    {
-      result->SetBinContent(i+1, results[i] * results[i] / aEfficiency->GetBinContent(i+1));
-      // error is : (relError) * (value) = (minuit->GetParError(i) / minuit->GetParameter(i)) * (minuit->GetParameter(i) * minuit->GetParameter(i))
-      result->SetBinError(i+1, minuit->GetParError(i) * results[i] /  aEfficiency->GetBinContent(i+1));
-    }
-    else
     {
-      result->SetBinContent(i+1, 0);
-      result->SetBinError(i+1, 0);
+      // set minimum value to prevent MINUIT just staying in the small value
+      for (Int_t i=1; i<=initialConditions->GetNbinsX(); i++)
+        initialConditions->SetBinContent(i, TMath::Max(initialConditions->GetBinContent(i), 1e-3));
     }
   }
-  //DrawGuess(results);
 
-  delete[] results;
-
-  return status;
+  return AliUnfolding::Unfold(fCurrentCorrelation, fCurrentEfficiency, fCurrentESD, initialConditions, fMultiplicityESDCorrected[correlationID], check);
 }
 
 //____________________________________________________________________
-Int_t AliMultiplicityCorrection::ApplyMinuitFit(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t check, TH1* initialConditions)
+Int_t AliMultiplicityCorrection::ApplyNBDFit(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType)
 {
   //
-  // correct spectrum using minuit chi2 method
-  //
-  // for description of parameters, see UnfoldWithMinuit
-  //
-
-  Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
-  SetupCurrentHists(inputRange, fullPhaseSpace, eventType, fgCreateBigBin);
-
-  return UnfoldWithMinuit(fCurrentCorrelation, fCurrentEfficiency, fCurrentESD, initialConditions, fMultiplicityESDCorrected[correlationID], check);
-}
-
-//____________________________________________________________________
-void AliMultiplicityCorrection::ApplyNBDFit(Int_t inputRange, Bool_t fullPhaseSpace)
-{
+  // correct spectrum using minuit chi2 method with a NBD function
   //
-  // correct spectrum using minuit chi2 method applying a NBD fit
+  // for description of parameters, see AliUnfolding::Unfold
   //
 
   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
-
-  SetupCurrentHists(inputRange, fullPhaseSpace, kTrVtx, kFALSE);
-  //normalize ESD
-  fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
-
-  fgCorrelationMatrix = new TMatrixD(fgkMaxParams, fgkMaxParams);
-  fgCorrelationCovarianceMatrix = new TMatrixD(fgkMaxParams, fgkMaxParams);
-  fgCurrentESDVector = new TVectorD(fgkMaxParams);
-
-  // normalize correction for given nPart
-  for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
-  {
-    Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
-    if (sum <= 0)
-      continue;
-    for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
-    {
-      // npart sum to 1
-      fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
-      fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
-
-      if (i <= fgkMaxParams && j <= fgkMaxParams)
-        (*fgCorrelationMatrix)(j-1, i-1) = fCurrentCorrelation->GetBinContent(i, j);
-    }
-  }
-
-  for (Int_t i=0; i<fgkMaxParams; ++i)
-  {
-    (*fgCurrentESDVector)[i] = fCurrentESD->GetBinContent(i+1);
-    if (fCurrentESD->GetBinError(i+1) > 0)
-      (*fgCorrelationCovarianceMatrix)(i, i) = 1.0 / (fCurrentESD->GetBinError(i+1) * fCurrentESD->GetBinError(i+1));
-  }
-
-  // Create NBD function
-  if (!fgNBD)
-  {
-    fgNBD = new TF1("nbd", "[0] * TMath::Binomial([2]+TMath::Nint(x)-1, [2]-1) * pow([1] / ([1]+[2]), TMath::Nint(x)) * pow(1 + [1]/[2], -[2])", 0, 250);
-    fgNBD->SetParNames("scaling", "averagen", "k");
-  }
-
-  // Initialize TMinuit via generic fitter interface
-  TVirtualFitter *minuit = TVirtualFitter::Fitter(0, 3);
-
-  minuit->SetFCN(MinuitNBD);
-  minuit->SetParameter(0, "param0", 1, 0.1, 0.001, fCurrentESD->GetMaximum() * 1000);
-  minuit->SetParameter(1, "param1", 10, 1, 0.001, 1000);
-  minuit->SetParameter(2, "param2", 2, 1, 0.001, 1000);
-
-  Double_t arglist[100];
-  arglist[0] = 0;
-  minuit->ExecuteCommand("SET PRINT", arglist, 1);
-  minuit->ExecuteCommand("MIGRAD", arglist, 0);
-  //minuit->ExecuteCommand("MINOS", arglist, 0);
-
-  fgNBD->SetParameters(minuit->GetParameter(0), minuit->GetParameter(1), minuit->GetParameter(2));
-
-  for (Int_t i=0; i<fgkMaxParams; ++i)
-  {
-    printf("%d : %f\n", i, fgNBD->Eval(i));
-    if (fgNBD->Eval(i) > 0)
-    {
-      fMultiplicityESDCorrected[correlationID]->SetBinContent(i+1, fgNBD->Eval(i));
-      fMultiplicityESDCorrected[correlationID]->SetBinError(i+1, 0);
-    }
-  }
+  
+  AliUnfolding::SetUnfoldingMethod(AliUnfolding::kFunction);
+  SetupCurrentHists(inputRange, fullPhaseSpace, eventType);
+  
+       TF1* func = new TF1("nbd", "[0] * TMath::Gamma([2]+x) / TMath::Gamma([2]) / TMath::Gamma(x+1) * pow([1] / ([1]+[2]), x) * pow(1.0 + [1]/[2], -[2])");
+       func->SetParNames("scaling", "averagen", "k");
+       func->SetParLimits(0, 0, 1000);
+       func->SetParLimits(1, 1, 50);
+       func->SetParLimits(2, 1, 10);
+       func->SetParameters(1, 10, 2);
+  AliUnfolding::SetFunction(func);
+  
+  return AliUnfolding::Unfold(fCurrentCorrelation, fCurrentEfficiency, fCurrentESD, 0, fMultiplicityESDCorrected[correlationID]);
 }
 
 //____________________________________________________________________
@@ -1449,6 +804,17 @@ void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRang
 
   TH1* residualHist = new TH1F("residualHist", "residualHist", 51, -5, 5);
 
+  // find bin limit
+  Int_t lastBin = 0;
+  for (Int_t i=1; i<=fCurrentESD->GetNbinsX(); ++i)
+  {
+    if (fCurrentESD->GetBinContent(i) <= 5)
+    {
+      lastBin = i;
+      break;
+    }
+  }
+  
   // TODO fix errors
   Float_t chi2 = 0;
   for (Int_t i=1; i<=residual->GetNbinsX(); ++i)
@@ -1456,7 +822,7 @@ void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRang
     if (fCurrentESD->GetBinError(i) > 0)
     {
       Float_t value = residual->GetBinContent(i) / fCurrentESD->GetBinError(i);
-      if (i > 1 && i <= fLastBinLimit)
+      if (i > 1 && i <= lastBin)
         chi2 += value * value;
       Printf("%d --> %f (%f)", i, value * value, chi2);
       residual->SetBinContent(i, value);
@@ -1477,7 +843,7 @@ void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRang
   //residualHist->Fit("gaus", "N");
   //delete residualHist;
 
-  printf("Difference (Residuals) is %f for bin 2-%d\n", fLastChi2Residuals, fLastBinLimit);
+  printf("Difference (Residuals) is %f for bin 2-%d\n", fLastChi2Residuals, lastBin);
 
   TCanvas* canvas1 = 0;
   if (simple)
@@ -1796,7 +1162,7 @@ void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRang
     Double_t newChi2 = 0;
     Double_t newChi2Limit150 = 0;
     Int_t ndf = 0;
-    for (Int_t i=1; i<=TMath::Min(diffMCUnfolded2->GetNbinsX(), fgkMaxParams+1); ++i)
+    for (Int_t i=1; i<=diffMCUnfolded2->GetNbinsX(); ++i)
     {
       Double_t value = 0;
       if (proj->GetBinError(i) > 0)
@@ -2027,7 +1393,7 @@ TH1* AliMultiplicityCorrection::CalculateStdDev(TH1** results, Int_t max)
 }
 
 //____________________________________________________________________
-TH1* AliMultiplicityCorrection::StatisticalUncertainty(MethodType methodType, Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t randomizeMeasured, Bool_t randomizeResponse, TH1* compareTo)
+TH1* AliMultiplicityCorrection::StatisticalUncertainty(AliUnfolding::MethodType methodType, Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t randomizeMeasured, Bool_t randomizeResponse, TH1* compareTo)
 {
   //
   // evaluates the uncertainty that arises from the non-infinite statistics in the response matrix
@@ -2043,6 +1409,10 @@ TH1* AliMultiplicityCorrection::StatisticalUncertainty(MethodType methodType, In
 
   // initialize seed with current time
   gRandom->SetSeed(0);
+  
+  if (methodType == AliUnfolding::kChi2Minimization)
+    AliUnfolding::SetCreateOverflowBin(5);
+  AliUnfolding::SetUnfoldingMethod(methodType);
 
   const Int_t kErrorIterations = 150;
 
@@ -2055,12 +1425,7 @@ TH1* AliMultiplicityCorrection::StatisticalUncertainty(MethodType methodType, In
   {
     Printf("Iteration %d of %d...", n, kErrorIterations);
 
-    // only done for chi2 minimization
-    Bool_t createBigBin = kFALSE;
-    if (methodType == kChi2Minimization)
-      createBigBin = kTRUE;
-
-    SetupCurrentHists(inputRange, fullPhaseSpace, eventType, createBigBin);
+    SetupCurrentHists(inputRange, fullPhaseSpace, eventType);
 
     TH1* measured = (TH1*) fCurrentESD->Clone("measured");
 
@@ -2087,7 +1452,7 @@ TH1* AliMultiplicityCorrection::StatisticalUncertainty(MethodType methodType, In
     }
 
     // only for bayesian method we have to do it before the call to Unfold...
-    if (methodType == kBayesian)
+    if (methodType == AliUnfolding::kBayesian)
     {
       for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
       {
@@ -2127,13 +1492,7 @@ TH1* AliMultiplicityCorrection::StatisticalUncertainty(MethodType methodType, In
     {
       result = (TH1*) fMultiplicityESDCorrected[correlationID]->Clone(Form("result_%d", n));
 
-      Int_t returnCode = -1;
-      if (methodType == kChi2Minimization)
-      {
-        returnCode = UnfoldWithMinuit(fCurrentCorrelation, fCurrentEfficiency, measured, 0, result, kFALSE);
-      }
-      else if (methodType == kBayesian)
-        returnCode = UnfoldWithBayesian(fCurrentCorrelation, fCurrentEfficiency, measured, 0, result, fgBayesianSmoothing, fgBayesianIterations);
+      Int_t returnCode = AliUnfolding::Unfold(fCurrentCorrelation, fCurrentEfficiency, measured, 0, result);
 
       if (returnCode != 0)
         return 0;
@@ -2286,7 +1645,7 @@ void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t ful
   // initialize seed with current time
   gRandom->SetSeed(0);
 
-  SetupCurrentHists(inputRange, fullPhaseSpace, eventType, kFALSE);
+  SetupCurrentHists(inputRange, fullPhaseSpace, eventType);
 
   // normalize correction for given nPart
   for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
@@ -2317,7 +1676,9 @@ void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t ful
 
   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
 
-  if (UnfoldWithBayesian(fCurrentCorrelation, fCurrentEfficiency, fCurrentESD, initialConditions, fMultiplicityESDCorrected[correlationID], regPar, nIterations) != 0)
+  AliUnfolding::SetBayesianParameters(regPar, nIterations);
+  AliUnfolding::SetUnfoldingMethod(AliUnfolding::kBayesian);
+  if (AliUnfolding::Unfold(fCurrentCorrelation, fCurrentEfficiency, fCurrentESD, initialConditions, fMultiplicityESDCorrected[correlationID]) != 0)
     return;
 
   if (!determineError)
@@ -2350,7 +1711,7 @@ void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t ful
 
     TH1* result2 = (TH1*) fMultiplicityESDCorrected[correlationID]->Clone("result2");
     result2->Reset();
-    if (UnfoldWithBayesian(fCurrentCorrelation, fCurrentEfficiency, randomized, initialConditions, result2, regPar, nIterations) != 0)
+    if (AliUnfolding::Unfold(fCurrentCorrelation, fCurrentEfficiency, randomized, initialConditions, result2) != 0)
       return;
 
     result2->Scale(1.0 / result2->Integral());
@@ -2372,300 +1733,6 @@ void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t ful
 }
 
 //____________________________________________________________________
-Int_t AliMultiplicityCorrection::UnfoldWithBayesian(TH1* correlation, TH1* aEfficiency, TH1* measured, TH1* initialConditions, TH1* aResult, Float_t regPar, Int_t nIterations)
-{
-  //
-  // unfolds a spectrum
-  //
-  
-  if (measured->Integral() <= 0)
-  {
-    Printf("AliMultiplicityCorrection::UnfoldWithBayesian: ERROR: The measured spectrum is empty");
-    return -1;
-  }
-
-  const Int_t kStartBin = 0;
-
-  const Int_t kMaxM = fgkMaxInput;  //<= fCurrentCorrelation->GetNbinsY(); // max measured axis
-  const Int_t kMaxT = fgkMaxParams; //<= fCurrentCorrelation->GetNbinsX(); // max true axis
-
-  // convergence limit: kMaxT * 0.001^2 = kMaxT * 1e-6 (e.g. 250 bins --> 2.5 e-4)
-  const Double_t kConvergenceLimit = kMaxT * 1e-6;
-
-  // store information in arrays, to increase processing speed (~ factor 5)
-  Double_t measuredCopy[kMaxM];
-  Double_t measuredError[kMaxM];
-  Double_t prior[kMaxT];
-  Double_t result[kMaxT];
-  Double_t efficiency[kMaxT];
-
-  Double_t response[kMaxT][kMaxM];
-  Double_t inverseResponse[kMaxT][kMaxM];
-
-  // for normalization
-  Float_t measuredIntegral = measured->Integral();
-  for (Int_t m=0; m<kMaxM; m++)
-  {
-    measuredCopy[m] = measured->GetBinContent(m+1) / measuredIntegral;
-    measuredError[m] = measured->GetBinError(m+1) / measuredIntegral;
-
-    for (Int_t t=0; t<kMaxT; t++)
-    {
-      response[t][m] = correlation->GetBinContent(t+1, m+1);
-      inverseResponse[t][m] = 0;
-    }
-  }
-
-  for (Int_t t=0; t<kMaxT; t++)
-  {
-    efficiency[t] = aEfficiency->GetBinContent(t+1);
-    prior[t] = measuredCopy[t];
-    result[t] = 0;
-  }
-
-  // pick prior distribution
-  if (initialConditions)
-  {
-    printf("Using different starting conditions...\n");
-    // for normalization
-    Float_t inputDistIntegral = initialConditions->Integral();
-    for (Int_t i=0; i<kMaxT; i++)
-      prior[i] = initialConditions->GetBinContent(i+1) / inputDistIntegral;
-  }
-
-  //TH1F* convergence = new TH1F("convergence", "convergence", 200, 0.5, 200.5);
-  
-  // unfold...
-  for (Int_t i=0; i<nIterations || nIterations < 0; i++)
-  {
-    //printf(" iteration %i \n", i);
-
-    // calculate IR from Beyes theorem
-    // IR_ji = R_ij * prior_i / sum_k(R_kj * prior_k)
-
-    Double_t chi2Measured = 0;
-    for (Int_t m=0; m<kMaxM; m++)
-    {
-      Float_t norm = 0;
-      for (Int_t t = kStartBin; t<kMaxT; t++)
-        norm += response[t][m] * prior[t];
-
-      // calc. chi2: (measured - response * prior) / error
-      if (measuredError[m] > 0)
-      {
-        Double_t value = (measuredCopy[m] - norm) / measuredError[m];
-        chi2Measured += value * value;
-      }
-
-      if (norm > 0)
-      {
-        for (Int_t t = kStartBin; t<kMaxT; t++)
-          inverseResponse[t][m] = response[t][m] * prior[t] / norm;
-      }
-      else
-      {
-        for (Int_t t = kStartBin; t<kMaxT; t++)
-          inverseResponse[t][m] = 0;
-      }
-    }
-    //Printf("chi2Measured of the last prior is %e", chi2Measured);
-
-    for (Int_t t = kStartBin; t<kMaxT; t++)
-    {
-      Float_t value = 0;
-      for (Int_t m=0; m<kMaxM; m++)
-        value += inverseResponse[t][m] * measuredCopy[m];
-
-      if (efficiency[t] > 0)
-        result[t] = value / efficiency[t];
-      else
-        result[t] = 0;
-    }
-    
-    // draw intermediate result
-    /*
-    for (Int_t t=0; t<kMaxT; t++)
-      aResult->SetBinContent(t+1, result[t]);
-    aResult->SetMarkerStyle(20+i);
-    aResult->SetMarkerColor(2);
-    aResult->DrawCopy("P SAME HIST");
-    */
-
-    Double_t chi2LastIter = 0;
-    // regularization (simple smoothing)
-    for (Int_t t=kStartBin; t<kMaxT; t++)
-    {
-      Float_t newValue = 0;
-      
-      // 0 bin excluded from smoothing
-      if (t > kStartBin+2 && t<kMaxT-1)
-      {
-        Float_t average = (result[t-1] + result[t] + result[t+1]) / 3;
-
-        // weight the average with the regularization parameter
-        newValue = (1 - regPar) * result[t] + regPar * average;
-      }
-      else
-        newValue = result[t];
-
-      // calculate chi2 (change from last iteration)
-      if (prior[t] > 1e-5)
-      {
-        Double_t diff = (prior[t] - newValue) / prior[t];
-        chi2LastIter += diff * diff;
-      }
-
-      prior[t] = newValue;
-    }
-    //printf("Chi2 of %d iteration = %e\n", i, chi2LastIter);
-    //convergence->Fill(i+1, chi2LastIter);
-
-    if (nIterations < 0 && chi2LastIter < kConvergenceLimit)
-    {
-      Printf("Stopped Bayesian unfolding after %d iterations at chi2(change since last iteration) of %e; chi2Measured of the last prior is %e", i, chi2LastIter, chi2Measured);
-      break;
-    }
-
-    //if (i % 5 == 0)
-    //  DrawComparison(Form("Bayesian_%d", i), mcTarget, correlationID, kTRUE, eventType);
-  } // end of iterations
-
-  //new TCanvas; convergence->DrawCopy(); gPad->SetLogy();
-  //delete convergence;
-
-  for (Int_t t=0; t<kMaxT; t++)
-    aResult->SetBinContent(t+1, result[t]);
-
-  return 0;
-
-  // ********
-  // Calculate the covariance matrix, all arguments are taken from NIM,A362,487-498,1995
-
-  /*printf("Calculating covariance matrix. This may take some time...\n");
-
-  // TODO check if this is the right one...
-  TH1* sumHist = GetMultiplicityMC(inputRange, eventType)->ProjectionY("sumHist", 1, GetMultiplicityMC(inputRange, eventType)->GetNbinsX());
-
-  Int_t xBins = hInverseResponseBayes->GetNbinsX();
-  Int_t yBins = hInverseResponseBayes->GetNbinsY();
-
-  // calculate "unfolding matrix" Mij
-  Float_t matrixM[251][251];
-  for (Int_t i=1; i<=xBins; i++)
-  {
-    for (Int_t j=1; j<=yBins; j++)
-    {
-      if (fCurrentEfficiency->GetBinContent(i) > 0)
-        matrixM[i-1][j-1] = hInverseResponseBayes->GetBinContent(i, j) / fCurrentEfficiency->GetBinContent(i);
-      else
-        matrixM[i-1][j-1] = 0;
-    }
-  }
-
-  Float_t* vectorn = new Float_t[yBins];
-  for (Int_t j=1; j<=yBins; j++)
-    vectorn[j-1] = fCurrentESD->GetBinContent(j);
-
-  // first part of covariance matrix, depends on input distribution n(E)
-  Float_t cov1[251][251];
-
-  Float_t nEvents = fCurrentESD->Integral(); // N
-
-  xBins = 20;
-  yBins = 20;
-
-  for (Int_t k=0; k<xBins; k++)
-  {
-    printf("In Cov1: %d\n", k);
-    for (Int_t l=0; l<yBins; l++)
-    {
-      cov1[k][l] = 0;
-
-      // sum_j Mkj Mlj n(Ej) * (1 - n(Ej) / N)
-      for (Int_t j=0; j<yBins; j++)
-        cov1[k][l] += matrixM[k][j] * matrixM[l][j] * vectorn[j]
-          * (1.0 - vectorn[j] / nEvents);
-
-      // - sum_i,j (i != j) Mki Mlj n(Ei) n(Ej) / N
-      for (Int_t i=0; i<yBins; i++)
-        for (Int_t j=0; j<yBins; j++)
-        {
-          if (i == j)
-            continue;
-          cov1[k][l] -= matrixM[k][i] * matrixM[l][j] * vectorn[i]
-            * vectorn[j] / nEvents;
-         }
-    }
-  }
-
-  printf("Cov1 finished\n");
-
-  TH2F* cov = (TH2F*) hInverseResponseBayes->Clone("cov");
-  cov->Reset();
-
-  for (Int_t i=1; i<=xBins; i++)
-    for (Int_t j=1; j<=yBins; j++)
-      cov->SetBinContent(i, j, cov1[i-1][j-1]);
-
-  new TCanvas;
-  cov->Draw("COLZ");
-
-  // second part of covariance matrix, depends on response matrix
-  Float_t cov2[251][251];
-
-  // Cov[P(Er|Cu), P(Es|Cu)] term
-  Float_t covTerm[100][100][100];
-  for (Int_t r=0; r<yBins; r++)
-    for (Int_t u=0; u<xBins; u++)
-      for (Int_t s=0; s<yBins; s++)
-      {
-        if (r == s)
-          covTerm[r][u][s] = 1.0 / sumHist->GetBinContent(u+1) * hResponse->GetBinContent(u+1, r+1)
-            * (1.0 - hResponse->GetBinContent(u+1, r+1));
-        else
-          covTerm[r][u][s] = - 1.0 / sumHist->GetBinContent(u+1) * hResponse->GetBinContent(u+1, r+1)
-            * hResponse->GetBinContent(u+1, s+1);
-      }
-
-  for (Int_t k=0; k<xBins; k++)
-    for (Int_t l=0; l<yBins; l++)
-    {
-      cov2[k][l] = 0;
-      printf("In Cov2: %d %d\n", k, l);
-      for (Int_t i=0; i<yBins; i++)
-        for (Int_t j=0; j<yBins; j++)
-        {
-          //printf("In Cov2: %d %d %d %d\n", k, l, i, j);
-          // calculate Cov(Mki, Mlj) = sum{ru},{su} ...
-          Float_t tmpCov = 0;
-          for (Int_t r=0; r<yBins; r++)
-            for (Int_t u=0; u<xBins; u++)
-              for (Int_t s=0; s<yBins; s++)
-              {
-                if (hResponse->GetBinContent(u+1, r+1) == 0 || hResponse->GetBinContent(u+1, s+1) == 0
-                  || hResponse->GetBinContent(u+1, i+1) == 0)
-                  continue;
-
-                tmpCov += BayesCovarianceDerivate(matrixM, hResponse, fCurrentEfficiency, k, i, r, u)
-                  * BayesCovarianceDerivate(matrixM, hResponse, fCurrentEfficiency, l, j, s, u)
-                  * covTerm[r][u][s];
-              }
-
-          cov2[k][l] += fCurrentESD->GetBinContent(i+1) * fCurrentESD->GetBinContent(j+1) * tmpCov;
-        }
-    }
-
-  printf("Cov2 finished\n");
-
-  for (Int_t i=1; i<=xBins; i++)
-    for (Int_t j=1; j<=yBins; j++)
-      cov->SetBinContent(i, j, cov1[i-1][j-1] + cov2[i-1][j-1]);
-
-  new TCanvas;
-  cov->Draw("COLZ");*/
-}
-
-//____________________________________________________________________
 Float_t AliMultiplicityCorrection::BayesCovarianceDerivate(Float_t matrixM[251][251], TH2* hResponse, Int_t k, Int_t i, Int_t r, Int_t u)
 {
   //
@@ -2700,7 +1767,7 @@ void AliMultiplicityCorrection::ApplyLaszloMethod(Int_t inputRange, Bool_t fullP
   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
   Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4);
 
-  SetupCurrentHists(inputRange, fullPhaseSpace, eventType, kFALSE);
+  SetupCurrentHists(inputRange, fullPhaseSpace, eventType);
   //normalize ESD
   fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
 
@@ -2822,7 +1889,7 @@ void AliMultiplicityCorrection::ApplyGaussianMethod(Int_t inputRange, Bool_t ful
   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
   Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4);
 
-  SetupCurrentHists(inputRange, fullPhaseSpace, kTrVtx, kFALSE);
+  SetupCurrentHists(inputRange, fullPhaseSpace, kTrVtx);
   //normalize ESD
   fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
 
index fede05b..63204d3 100644 (file)
@@ -27,12 +27,11 @@ class TCollection;
 #include <TMatrixD.h>
 #include <TVectorD.h>
 #include <AliPWG0Helper.h>
+#include <AliUnfolding.h>
 
 class AliMultiplicityCorrection : public TNamed {
   public:
     enum EventType { kTrVtx = 0, kMB, kINEL, kNSD };
-    enum RegularizationType { kNone = 0, kPol0, kPol1, kLog, kEntropy, kCurvature };
-    enum MethodType { kChi2Minimization = 0, kBayesian = 1 };
     enum { kESDHists = 4, kMCHists = 5, kCorrHists = 8, kQualityRegions = 3 };
 
     AliMultiplicityCorrection();
@@ -55,17 +54,12 @@ class AliMultiplicityCorrection : public TNamed {
 
     Int_t ApplyMinuitFit(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t check = kFALSE, TH1* initialConditions = 0);
 
-    static void SetRegularizationParameters(RegularizationType type, Float_t weight, Int_t minuitParams = -1);
-    static void SetBayesianParameters(Float_t smoothing, Int_t nIterations);
-    static void SetCreateBigBin(Bool_t flag) { fgCreateBigBin = flag; }
-
-    void ApplyNBDFit(Int_t inputRange, Bool_t fullPhaseSpace);
-
     void ApplyBayesianMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Float_t regPar = 1, Int_t nIterations = 100, TH1* initialConditions = 0, Bool_t determineError = kTRUE);
 
     static TH1* CalculateStdDev(TH1** results, Int_t max);
-    TH1* StatisticalUncertainty(MethodType methodType, Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t randomizeMeasured, Bool_t randomizeResponse, TH1* compareTo = 0);
+    TH1* StatisticalUncertainty(AliUnfolding::MethodType methodType, Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t randomizeMeasured, Bool_t randomizeResponse, TH1* compareTo = 0);
 
+    Int_t ApplyNBDFit(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType);
     void ApplyGaussianMethod(Int_t inputRange, Bool_t fullPhaseSpace);
 
     void ApplyLaszloMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType);
@@ -102,48 +96,14 @@ class AliMultiplicityCorrection : public TNamed {
     void FFT(Int_t dir, Int_t m, Double_t *x, Double_t *y);
 
   protected:
-    static const Int_t fgkMaxParams;  //! bins in unfolded histogram = number of fit params
-    static const Int_t fgkMaxInput;   //! bins in measured histogram
-
-    static Double_t RegularizationPol0(TVectorD& params);
-    static Double_t RegularizationPol1(TVectorD& params);
-    static Double_t RegularizationTotalCurvature(TVectorD& params);
-    static Double_t RegularizationEntropy(TVectorD& params);
-    static Double_t RegularizationLog(TVectorD& params);
-
-    static void MinuitFitFunction(Int_t&, Double_t*, Double_t& chi2, Double_t *params, Int_t);
-    static void MinuitNBD(Int_t& unused1, Double_t* unused2, Double_t& chi2, Double_t *params, Int_t unused3);
-
-    static void DrawGuess(Double_t *params);
-
-    void SetupCurrentHists(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t createBigBin);
+    void SetupCurrentHists(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType);
 
     Float_t BayesCovarianceDerivate(Float_t matrixM[251][251], TH2* hResponse, Int_t k, Int_t i, Int_t r, Int_t u);
-    static Int_t UnfoldWithBayesian(TH1* correlation, TH1* aEfficiency, TH1* measured, TH1* initialConditions, TH1* aResult, Float_t regPar, Int_t nIterations);
-    static Int_t UnfoldWithMinuit(TH1* correlation, TH1* aEfficiency, TH1* measured, TH1* initialConditions, TH1* result, Bool_t check);
-
+    
     TH1* fCurrentESD;         //! current input esd
-    TH1* fCurrentCorrelation; //! current correlation
+    TH2* fCurrentCorrelation; //! current correlation
     TH1* fCurrentEfficiency;  //! current efficiency
 
-    // static variable to be accessed by MINUIT
-    static TMatrixD* fgCorrelationMatrix;            //! contains fCurrentCorrelation in matrix form
-    static TMatrixD* fgCorrelationCovarianceMatrix;  //! contains the errors of fCurrentESD
-    static TVectorD* fgCurrentESDVector;             //! contains fCurrentESD
-    static TVectorD* fgEntropyAPriori;               //! a-priori distribution for entropy regularization
-
-    static TF1* fgNBD;   //! negative binomial distribution
-
-    // configuration params follow
-    static RegularizationType fgRegularizationType; //! regularization that is used during Chi2 method
-    static Float_t fgRegularizationWeight;          //! factor for regularization term
-    static Bool_t  fgCreateBigBin;                  //! to fix fluctuations at high multiplicities, all entries above a certain limit are summarized in one bin
-    static Int_t   fgNParamsMinuit;                 //! number of parameters minuit uses for unfolding (todo: to be merged w/ fgkMaxParams that has to be const. for the moment)
-
-    static Float_t fgBayesianSmoothing;             //! smoothing parameter (0 = no smoothing)
-    static Int_t   fgBayesianIterations;            //! number of iterations in Bayesian method
-    // end of configuration
-
     TH2F* fMultiplicityESD[kESDHists]; // multiplicity histogram: vtx vs multiplicity; array: |eta| < 0.5, 1.0, 1.5, 2 (0..3)
 
     TH2F* fMultiplicityVtx[kMCHists];  // multiplicity histogram of events that have a reconstructed vertex : vtx vs multiplicity; array: |eta| < 0.5, 1.0, 1.5, 2, inf (0..4)
index c692553..25fd74f 100644 (file)
@@ -88,7 +88,7 @@ void correct(const char* fileNameMC = "multiplicityMC.root", const char* folder
 
   if (chi2)
   {
-    mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol0, beta);
+    AliUnfolding::SetChi2Regularization(AliUnfolding::kPol0, beta);
     //mult->SetCreateBigBin(kFALSE);
     //mult->SetRegularizationParameters(AliMultiplicityCorrection::kNone, 0); //mult->SetCreateBigBin(kFALSE);
     //mult->SetRegularizationParameters(AliMultiplicityCorrection::kNone, 0, 125); mult->SetCreateBigBin(kFALSE);
@@ -99,6 +99,7 @@ void correct(const char* fileNameMC = "multiplicityMC.root", const char* folder
     //mult->SetMultiplicityESDCorrected(histID, (TH1F*) mcCompare);
     
     mult->ApplyMinuitFit(histID, fullPhaseSpace, eventType, kFALSE); //hist2->ProjectionY("mymchist"));
+    //mult->ApplyNBDFit(histID, fullPhaseSpace, eventType);
   }
   else
   {