refactoring multiplicity analysis
[u/mrichter/AliRoot.git] / PWG0 / multiplicity / AliMultiplicityCorrection.cxx
index cc6f34aeaf385fa2c4e33a6ba1dd33be36f2a6ad..45d42174fc594fecc4e3e7270b74123aad104752 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());
@@ -2371,300 +1732,6 @@ void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t ful
   delete error;
 }
 
-//____________________________________________________________________
-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());