o) Performance studies added (mainly in runMultiplicitySelector) that evaluate the...
authorjgrosseo <jgrosseo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 11 Apr 2007 13:18:30 +0000 (13:18 +0000)
committerjgrosseo <jgrosseo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 11 Apr 2007 13:18:30 +0000 (13:18 +0000)
parameters of the chi2 and bayesian method.
o) Cleaned up AliMultiplicityCorrection
o) Implemented selector to read the spectrum from the ESD

PWG0/dNdEta/AliMultiplicityCorrection.cxx
PWG0/dNdEta/AliMultiplicityCorrection.h
PWG0/dNdEta/AliMultiplicityESDSelector.cxx
PWG0/dNdEta/AliMultiplicityESDSelector.h
PWG0/dNdEta/AliMultiplicityMCSelector.cxx
PWG0/dNdEta/AliMultiplicityMCSelector.h
PWG0/dNdEta/runMultiplicitySelector.C

index 3f28db6..bdb827f 100644 (file)
@@ -46,10 +46,11 @@ TMatrixF* AliMultiplicityCorrection::fCorrelationCovarianceMatrix = 0;
 TVectorF* AliMultiplicityCorrection::fCurrentESDVector = 0;
 AliMultiplicityCorrection::RegularizationType AliMultiplicityCorrection::fRegularizationType = AliMultiplicityCorrection::kPol1;
 Float_t AliMultiplicityCorrection::fRegularizationWeight = 1e4;
+TF1* AliMultiplicityCorrection::fNBD = 0;
 
 //____________________________________________________________________
 AliMultiplicityCorrection::AliMultiplicityCorrection() :
-  TNamed(), fLastChi2MC(0), fLastChi2Residuals(0)
+  TNamed(), fLastChi2MC(0), fLastChi2MCLimit(0), fLastChi2Residuals(0)
 {
   //
   // default constructor
@@ -74,7 +75,7 @@ AliMultiplicityCorrection::AliMultiplicityCorrection() :
 
 //____________________________________________________________________
 AliMultiplicityCorrection::AliMultiplicityCorrection(const Char_t* name, const Char_t* title) :
-  TNamed(name, title), fLastChi2MC(0), fLastChi2Residuals(0)
+  TNamed(name, title), fLastChi2MC(0), fLastChi2MCLimit(0), fLastChi2Residuals(0)
 {
   //
   // named constructor
@@ -361,17 +362,14 @@ Double_t AliMultiplicityCorrection::RegularizationPol0(Double_t *params)
 
   for (Int_t i=1; i<fgMaxParams; ++i)
   {
-    if (params[i] == 0)
-      continue;
-
-    Double_t right  = params[i]; // / fCurrentESD->GetBinWidth(i+1);
-    Double_t left   = params[i-1]; // / fCurrentESD->GetBinWidth(i);
+    Double_t right  = params[i];
+    Double_t left   = params[i-1];
 
-    Double_t diff = (right - left) / right;
+    Double_t diff = (right - left);
     chi2 += diff * diff;
   }
 
-  return chi2;
+  return chi2 * 100;
 }
 
 //____________________________________________________________________
@@ -385,25 +383,48 @@ Double_t AliMultiplicityCorrection::RegularizationPol1(Double_t *params)
 
   for (Int_t i=2; i<fgMaxParams; ++i)
   {
-    if (params[i] == 0 || params[i-1] == 0)
+    if (params[i-1] == 0)
       continue;
 
-    Double_t right  = params[i]; // / fCurrentESD->GetBinWidth(i+1);
-    Double_t middle = params[i-1]; // / fCurrentESD->GetBinWidth(i);
-    Double_t left   = params[i-2]; // / fCurrentESD->GetBinWidth(i-1);
+    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;
+  }
 
-    Double_t der1 = (right - middle); // / fCurrentESD->GetBinWidth(i);
-    Double_t der2 = (middle - left); //  / fCurrentESD->GetBinWidth(i-1);
+  return chi2;
+}
 
-    Double_t diff = (der1 - der2) / middle; /// fCurrentESD->GetBinWidth(i);
+//____________________________________________________________________
+Double_t AliMultiplicityCorrection::RegularizationTest(Double_t *params)
+{
+  // homogenity term for minuit fitting
+  // pure function of the parameters
+  // prefers linear function (pol1)
 
-    // TODO give additonal weight to big bins
-    chi2 += diff * diff; // * fCurrentESD->GetBinWidth(i) * fCurrentESD->GetBinWidth(i-1);
+  Double_t chi2 = 0;
 
-    //printf("%d --> %f\n", i, diff);
+  Float_t der[fgMaxParams];
+
+  for (Int_t i=0; i<fgMaxParams-1; ++i)
+    der[i] = params[i+1] - params[i];
+
+  for (Int_t i=0; i<fgMaxParams-6; ++i)
+  {
+    for (Int_t j=1; j<=5; ++j)
+    {
+      Double_t diff = der[i+j] - der[i];
+      chi2 += diff * diff;
+    }
   }
 
-  return chi2; // 10000
+  return chi2 * 100; // 10000
 }
 
 //____________________________________________________________________
@@ -418,25 +439,20 @@ Double_t AliMultiplicityCorrection::RegularizationTotalCurvature(Double_t *param
 
   for (Int_t i=2; i<fgMaxParams; ++i)
   {
-    if (params[i] == 0 || params[i-1] == 0)
-      continue;
+    Double_t right  = params[i];
+    Double_t middle = params[i-1];
+    Double_t left   = params[i-2];
 
-    Double_t right  = params[i]; // / fCurrentESD->GetBinWidth(i+1);
-    Double_t middle = params[i-1]; // / fCurrentESD->GetBinWidth(i);
-    Double_t left   = params[i-2]; // / fCurrentESD->GetBinWidth(i-1);
+    Double_t der1 = (right - middle);
+    Double_t der2 = (middle - left);
 
-    Double_t der1 = (right - middle); // / fCurrentESD->GetBinWidth(i);
-    Double_t der2 = (middle - left); //  / fCurrentESD->GetBinWidth(i-1);
-
-    Double_t secDer = (der1 - der2); // / fCurrentESD->GetBinWidth(i);
+    Double_t secDer = (der1 - der2);
 
     // square and weight with the bin width
-    chi2 += secDer * secDer; // * fCurrentESD->GetBinWidth(i) * fCurrentESD->GetBinWidth(i-1);
-
-    //printf("%d --> %f\n", i, secDer);
+    chi2 += secDer * secDer;
   }
 
-  return chi2; // 10
+  return chi2 * 100;
 }
 
 //____________________________________________________________________
@@ -447,31 +463,41 @@ Double_t AliMultiplicityCorrection::RegularizationEntropy(Double_t *params)
   // calculates entropy, from
   // The method of reduced cross-entropy (M. Schmelling 1993)
 
-  //static Int_t callCount = 0;
-
   Double_t paramSum = 0;
   for (Int_t i=0; i<fgMaxParams; ++i)
-    paramSum += params[i]; // / fCurrentESD->GetBinWidth(i+1);
+    paramSum += params[i];
 
   Double_t chi2 = 0;
   for (Int_t i=0; i<fgMaxParams; ++i)
   {
-    Double_t tmp = params[i] / paramSum; // / fCurrentESD->GetBinWidth(i+1);
+    Double_t tmp = params[i] / paramSum;
     if (tmp > 0)
-    {
-      chi2 += tmp * log(tmp); /* * fCurrentESD->GetBinWidth(i+1);
-      //                                   /\
-      //                                    ------------------------------------
-      // TODO WARNING the absolute fitting values seem to depend on that value |
-      // this should be impossible...
-      //if ((callCount % 100) == 0)
-      //  printf("%f => %f\n", params[i], tmp * log(tmp)); */
-    }
+      chi2 += tmp * log(tmp);
   }
-  //if ((callCount++ % 100) == 0)
-  //  printf("\n");
 
-  return chi2; // 1000
+  return 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");
+  //
+
+  fNBD->SetParameters(params[0], params[1], params[2]);
+
+  Double_t params2[fgMaxParams];
+
+  for (Int_t i=0; i<fgMaxParams; ++i)
+    params2[i] = fNBD->Eval(i);
+
+  MinuitFitFunction(unused1, unused2, chi2, params2, unused3);
+
+  printf("%f %f %f --> %f\n", params[0], params[1], params[2], chi2);
 }
 
 //____________________________________________________________________
@@ -505,68 +531,16 @@ void AliMultiplicityCorrection::MinuitFitFunction(Int_t&, Double_t*, Double_t& c
   // (Ad - m) W (Ad - m)
   chi2FromFit = paramsVector * copy;
 
-  /*printf("chi2FromFit = %f\n", chi2FromFit);
-  chi2FromFit = 0;
-
-  // loop over multiplicity (x axis of fMultiplicityESD)
-  for (Int_t i=1; i<=fCurrentESD->GetNbinsX(); ++i)
-  {
-    if (i > fCurrentCorrelation->GetNbinsY() || i > fgMaxParams)
-      break;
-
-    Double_t sum = 0;
-    //Double_t error = 0;
-    // loop over generated particles (x axis of fCorrelation) that resulted in reconstruction of i tracks
-    for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsX(); ++j)
-    {
-      if (j > fgMaxParams)
-        break;
-
-      sum += fCurrentCorrelation->GetBinContent(j, i) * params[j-1];
-
-      //if (params[j-1] > 0)
-      //  error += fCurrentCorrelation->GetBinError(j, i) * fCurrentCorrelation->GetBinError(j, i) * params[j-1];
-    }
-
-    //printf("%f\n", sum);
-
-    Double_t diff = fCurrentESD->GetBinContent(i) - sum;
-
-    // include error
-    if (fCurrentESD->GetBinError(i) > 0)
-      diff /= fCurrentESD->GetBinError(i);
-
-    //if (fCurrentESD->GetBinContent(i) > 0)
-    //  diff /= fCurrentESD->GetBinContent(i);
-    //else
-    //  diff /= fCurrentESD->Integral();
-
-    // weight with bin width
-    //diff *= fCurrentESD->GetBinWidth(i);
-
-    //error = TMath::Sqrt(error) + fCurrentESD->GetBinError(i);
-    //if (error <= 0)
-    // error = 1;
-
-    //Double_t tmp = diff / error;
-    //chi2 += tmp * tmp;
-    chi2FromFit += diff * diff;
-
-    //printf("\nExpected sum = %f; Diff for bin %d is %f\n**********************************\n", fCurrentESD->GetBinContent(i), i, diff);
-    //printf("Diff for bin %d is %f\n", i, diff);
-  }
-
-  printf("chi2FromFit = %f\n", chi2FromFit);*/
-
   Double_t penaltyVal = 0;
 
   switch (fRegularizationType)
   {
-    case kNone:      break;
-    case kPol0:      penaltyVal = RegularizationPol0(params); break;
-    case kPol1:      penaltyVal = RegularizationPol1(params); break;
-    case kCurvature: penaltyVal = RegularizationTotalCurvature(params); break;
-    case kEntropy:   penaltyVal = RegularizationEntropy(params); break;
+    case kNone:       break;
+    case kPol0:       penaltyVal = RegularizationPol0(params); break;
+    case kPol1:       penaltyVal = RegularizationPol1(params); break;
+    case kCurvature:  penaltyVal = RegularizationTotalCurvature(params); break;
+    case kEntropy:    penaltyVal = RegularizationEntropy(params); break;
+    case kTest:    penaltyVal = RegularizationTest(params); break;
   }
 
   penaltyVal *= fRegularizationWeight;
@@ -617,7 +591,7 @@ void AliMultiplicityCorrection::SetupCurrentHists(Int_t inputRange, Bool_t fullP
 }
 
 //____________________________________________________________________
-void AliMultiplicityCorrection::ApplyMinuitFit(Int_t inputRange, Bool_t fullPhaseSpace, Bool_t check)
+void AliMultiplicityCorrection::ApplyMinuitFit(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t check, TH1* inputDist)
 {
   //
   // correct spectrum using minuit chi2 method
@@ -629,7 +603,7 @@ void AliMultiplicityCorrection::ApplyMinuitFit(Int_t inputRange, Bool_t fullPhas
   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
   Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4);
 
-  SetupCurrentHists(inputRange, fullPhaseSpace, kTrVtx);
+  SetupCurrentHists(inputRange, fullPhaseSpace, eventType);
 
   fCorrelationMatrix = new TMatrixF(fgMaxParams, fgMaxParams);
   fCorrelationCovarianceMatrix = new TMatrixF(fgMaxParams, fgMaxParams);
@@ -671,45 +645,51 @@ void AliMultiplicityCorrection::ApplyMinuitFit(Int_t inputRange, Bool_t fullPhas
         fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i-1, j));
   }*/
 
-  /*TCanvas* canvas1 = new TCanvas("ApplyMinuitFit", "ApplyMinuitFit", 800, 400);
-  canvas1->Divide(2, 1);
-  canvas1->cd(1);
-  fCurrentESD->DrawCopy();
-  canvas1->cd(2);
-  fCurrentCorrelation->DrawCopy("COLZ");*/
-
   // Initialize TMinuit via generic fitter interface
   TVirtualFitter *minuit = TVirtualFitter::Fitter(0, fgMaxParams);
 
   minuit->SetFCN(MinuitFitFunction);
 
   static Double_t results[fgMaxParams];
-  //printf("%x\n", (void*) results);
 
-  TH1* proj = fMultiplicityVtx[mcTarget]->ProjectionY();
+  if (inputDist)
+  {
+    printf("Using different starting conditions...\n");
+    new TCanvas;
+    inputDist->DrawCopy();
+  }
+  else
+    inputDist = fCurrentESD;
+
+  // normalize ESD
+  fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
+  inputDist->Scale(1.0 / inputDist->Integral());
+
+  Float_t minStartValue = 1e-3;
+
+  TH1* proj = fMultiplicityVtx[mcTarget]->ProjectionY("check-proj");
+  proj->Scale(1.0 / proj->Integral());
   for (Int_t i=0; i<fgMaxParams; ++i)
   {
-    //results[i] = 1;
-    results[i] = fCurrentESD->GetBinContent(i+1);
-    if (results[i] < 1e-2)
-      results[i] = 1e-2;
+    results[i] = inputDist->GetBinContent(i+1);
     if (check)
       results[i] = proj->GetBinContent(i+1);
-    minuit->SetParameter(i, Form("param%d", i), results[i], results[i] / 10, 0.001, fCurrentESD->GetMaximum() * 10);
+
+    printf("%f %f %f\n", inputDist->GetBinContent(i+1), TMath::Sqrt(inputDist->GetBinContent(i+1)), inputDist->GetBinError(i+1));
 
     (*fCurrentESDVector)[i] = fCurrentESD->GetBinContent(i+1);
     if (fCurrentESD->GetBinError(i+1) > 0)
-      (*fCorrelationCovarianceMatrix)(i, i) = 1.0 / (fCurrentESD->GetBinError(i+1) * fCurrentESD->GetBinError(i+1));
+      (*fCorrelationCovarianceMatrix)(i, i) = 1.0 / fCurrentESD->GetBinError(i+1);
 
-    //minuit->SetParameter(i, Form("param%d", i), fCurrentESD->GetBinWidth(i+1), 0.01, 0.001, 50);
+    minuit->SetParameter(i, Form("param%d", i), ((results[i] > minStartValue) ? results[i] : minStartValue), 0.1, 0, 1);
   }
-  minuit->SetParameter(0, "param0", results[1], ((results[1] > 1) ? (results[1] / 10) : 0), 0.001, fCurrentESD->GetMaximum() * 10);
-  //minuit->SetParameter(0, "param0", 1, 0, 1, 1);
+  // 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, 1);
 
   Int_t dummy;
   Double_t chi2;
   MinuitFitFunction(dummy, 0, chi2, results, 0);
-  printf("Chi2 of right solution is = %f\n", chi2);
+  printf("Chi2 of initial parameters is = %f\n", chi2);
 
   if (check)
     return;
@@ -718,16 +698,95 @@ void AliMultiplicityCorrection::ApplyMinuitFit(Int_t inputRange, Bool_t fullPhas
   arglist[0] = 0;
   minuit->ExecuteCommand("SET PRINT", arglist, 1);
   minuit->ExecuteCommand("MIGRAD", arglist, 1);
+  //printf("!!!!!!!!!!!!!! MIGRAD finished: Starting MINOS !!!!!!!!!!!!!!");
   //minuit->ExecuteCommand("MINOS", arglist, 0);
 
   for (Int_t i=0; i<fgMaxParams; ++i)
   {
-    results[i] = minuit->GetParameter(i);
-    fMultiplicityESDCorrected[correlationID]->SetBinContent(i+1, results[i]);
-    fMultiplicityESDCorrected[correlationID]->SetBinError(i+1, 0);
+    if (fCurrentEfficiency->GetBinContent(i+1) > 0)
+    {
+      fMultiplicityESDCorrected[correlationID]->SetBinContent(i+1, minuit->GetParameter(i) / fCurrentEfficiency->GetBinContent(i+1));
+      fMultiplicityESDCorrected[correlationID]->SetBinError(i+1, minuit->GetParError(i) / fCurrentEfficiency->GetBinContent(i+1));
+    }
+    else
+    {
+      fMultiplicityESDCorrected[correlationID]->SetBinContent(i+1, 0);
+      fMultiplicityESDCorrected[correlationID]->SetBinError(i+1, 0);
+    }
+  }
+}
+
+//____________________________________________________________________
+void AliMultiplicityCorrection::ApplyNBDFit(Int_t inputRange, Bool_t fullPhaseSpace)
+{
+  //
+  // correct spectrum using minuit chi2 method applying a NBD fit
+  //
+
+  Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
+
+  SetupCurrentHists(inputRange, fullPhaseSpace, kTrVtx);
+
+  fCorrelationMatrix = new TMatrixF(fgMaxParams, fgMaxParams);
+  fCorrelationCovarianceMatrix = new TMatrixF(fgMaxParams, fgMaxParams);
+  fCurrentESDVector = new TVectorF(fgMaxParams);
+
+  // 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 <= fgMaxParams && j <= fgMaxParams)
+        (*fCorrelationMatrix)(j-1, i-1) = fCurrentCorrelation->GetBinContent(i, j);
+    }
+  }
+
+  for (Int_t i=0; i<fgMaxParams; ++i)
+  {
+    (*fCurrentESDVector)[i] = fCurrentESD->GetBinContent(i+1);
+    if (fCurrentESD->GetBinError(i+1) > 0)
+      (*fCorrelationCovarianceMatrix)(i, i) = 1.0 / (fCurrentESD->GetBinError(i+1) * fCurrentESD->GetBinError(i+1));
+  }
+
+  // Create NBD function
+  if (!fNBD)
+  {
+    fNBD = 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);
+    fNBD->SetParNames("scaling", "averagen", "k");
   }
 
-  //printf("Penalty is %f\n", RegularizationTotalCurvature(results));
+  // 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, 1);
+  //minuit->ExecuteCommand("MINOS", arglist, 0);
+
+  fNBD->SetParameters(minuit->GetParameter(0), minuit->GetParameter(1), minuit->GetParameter(2));
+
+  for (Int_t i=0; i<fgMaxParams; ++i)
+  {
+    printf("%d : %f\n", i, fNBD->Eval(i));
+    if (fNBD->Eval(i) > 0)
+    {
+      fMultiplicityESDCorrected[correlationID]->SetBinContent(i+1, fNBD->Eval(i));
+      fMultiplicityESDCorrected[correlationID]->SetBinError(i+1, 0);
+    }
+  }
 }
 
 //____________________________________________________________________
@@ -761,20 +820,6 @@ void AliMultiplicityCorrection::NormalizeToBinWidth(TH2* hist)
 }
 
 //____________________________________________________________________
-void AliMultiplicityCorrection::ApplyMinuitFitAll()
-{
-  //
-  // fit all eta ranges
-  //
-
-  for (Int_t i=0; i<kESDHists; ++i)
-  {
-    ApplyMinuitFit(i, kFALSE);
-    ApplyMinuitFit(i, kTRUE);
-  }
-}
-
-//____________________________________________________________________
 void AliMultiplicityCorrection::DrawHistograms()
 {
   printf("ESD:\n");
@@ -830,11 +875,14 @@ void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRang
 
   // for that we convolute the response matrix with the gathered result
   // if normalizeESD is not set, the histogram is already normalized, this needs to be passed to CalculateMultiplicityESD
-  TH2* convoluted = CalculateMultiplicityESD(fMultiplicityESDCorrected[esdCorrId], esdCorrId, !normalizeESD);
+  TH1* tmpESDEfficiencyRecorrected = (TH1*) fMultiplicityESDCorrected[esdCorrId]->Clone("tmpESDEfficiencyRecorrected");
+  tmpESDEfficiencyRecorrected->Multiply(fCurrentEfficiency);
+  TH2* convoluted = CalculateMultiplicityESD(tmpESDEfficiencyRecorrected, esdCorrId, !normalizeESD);
   TH1* proj2 = convoluted->ProjectionY("proj2", -1, -1, "e");
   NormalizeToBinWidth(proj2);
-  TH1* residual = convoluted->ProjectionY("residual", -1, -1, "e");
-  residual->SetTitle("Residuals");
+
+  TH1* residual = (TH1*) proj2->Clone("residual");
+  residual->SetTitle("Residuals;Ntracks;(folded unfolded measured - measured) / measured");
 
   residual->Add(fCurrentESD, -1);
   residual->Divide(residual, fCurrentESD, 1, 1, "B");
@@ -846,6 +894,16 @@ void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRang
     residual->SetBinError(i, 0);
   }
 
+  // normalize mc to 1
+  //mcHist->Scale(1.0 / mcHist->Integral(1, 200));
+
+  // normalize unfolded esd to 1
+  //fMultiplicityESDCorrected[esdCorrId]->Scale(1.0 / fMultiplicityESDCorrected[esdCorrId]->Integral(1, 200));
+
+  // normalize unfolded result to mc hist
+  fMultiplicityESDCorrected[esdCorrId]->Scale(1.0 / fMultiplicityESDCorrected[esdCorrId]->Integral());
+  fMultiplicityESDCorrected[esdCorrId]->Scale(mcHist->Integral(1, 200));
+
   TCanvas* canvas1 = new TCanvas(tmpStr, tmpStr, 1200, 800);
   canvas1->Divide(2, 2);
 
@@ -856,7 +914,7 @@ void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRang
   if (normalizeESD)
     NormalizeToBinWidth(fMultiplicityESDCorrected[esdCorrId]);
 
-  proj->GetXaxis()->SetRangeUser(0, 150);
+  proj->GetXaxis()->SetRangeUser(0, 200);
   proj->DrawCopy("HIST");
   gPad->SetLogy();
 
@@ -865,23 +923,29 @@ void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRang
   fMultiplicityESDCorrected[esdCorrId]->DrawCopy("SAME HIST");
 
   Float_t chi2 = 0;
-  for (Int_t i=1; i<=100; ++i)
+  fLastChi2MCLimit = 0;
+  for (Int_t i=2; i<=200; ++i)
   {
-    if (fMultiplicityESDCorrected[esdCorrId]->GetBinContent(i) != 0)
+    if (proj->GetBinContent(i) != 0)
     {
-      Float_t value = (proj->GetBinContent(i) - fMultiplicityESDCorrected[esdCorrId]->GetBinContent(i)) / fMultiplicityESDCorrected[esdCorrId]->GetBinContent(i);
+      Float_t value = (proj->GetBinContent(i) - fMultiplicityESDCorrected[esdCorrId]->GetBinContent(i)) / proj->GetBinContent(i);
       chi2 += value * value;
+
+      if (chi2 < 0.1)
+        fLastChi2MCLimit = i;
+
+      if (i == 100)
+        fLastChi2MC = chi2;
     }
   }
 
-  printf("Difference (from MC) is %f for bin 1-100\n", chi2);
-  fLastChi2MC = chi2;
+  printf("Difference (from MC) is %f for bin 2-100. Limit is %d.\n", fLastChi2MC, fLastChi2MCLimit);
 
   canvas1->cd(2);
 
   NormalizeToBinWidth(fCurrentESD);
   gPad->SetLogy();
-  fCurrentESD->GetXaxis()->SetRangeUser(0, 150);
+  fCurrentESD->GetXaxis()->SetRangeUser(0, 200);
   //fCurrentESD->SetLineColor(2);
   fCurrentESD->DrawCopy("HIST");
 
@@ -891,7 +955,7 @@ void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRang
   proj2->DrawCopy("HIST SAME");
 
   chi2 = 0;
-  for (Int_t i=1; i<=100; ++i)
+  for (Int_t i=2; i<=100; ++i)
   {
     if (fCurrentESD->GetBinContent(i) != 0)
     {
@@ -900,15 +964,15 @@ void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRang
     }
   }
 
-  printf("Difference (Residuals) is %f for bin 1-100\n", chi2);
+  printf("Difference (Residuals) is %f for bin 2-100\n", chi2);
   fLastChi2Residuals = chi2;
 
   canvas1->cd(3);
   TH1* clone = dynamic_cast<TH1*> (proj->Clone("clone"));
   clone->Divide(fMultiplicityESDCorrected[esdCorrId]);
-  clone->SetTitle("Ratio;Npart;MC/ESD");
+  clone->SetTitle("Ratio;Npart;MC/Unfolded Measured");
   clone->GetYaxis()->SetRangeUser(0.8, 1.2);
-  clone->GetXaxis()->SetRangeUser(0, 150);
+  clone->GetXaxis()->SetRangeUser(0, 200);
   clone->DrawCopy("HIST");
 
   /*TLegend* legend = new TLegend(0.6, 0.7, 0.85, 0.85);
@@ -918,23 +982,26 @@ void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRang
   legend->Draw();*/
 
   canvas1->cd(4);
-
   residual->GetYaxis()->SetRangeUser(-0.2, 0.2);
-  residual->GetXaxis()->SetRangeUser(0, 150);
+  residual->GetXaxis()->SetRangeUser(0, 200);
   residual->DrawCopy();
 
   canvas1->SaveAs(Form("%s.gif", canvas1->GetName()));
 }
 
 //____________________________________________________________________
-void AliMultiplicityCorrection::GetComparisonResults(Float_t& mc, Float_t& residuals)
+void AliMultiplicityCorrection::GetComparisonResults(Float_t* mc, Int_t* mcLimit, Float_t* residuals)
 {
   // Returns the chi2 between the MC and the unfolded ESD as well as between the ESD and the folded unfolded ESD
   // These values are computed during DrawComparison, thus this function picks up the
   // last calculation
 
-  mc = fLastChi2MC;
-  residuals = fLastChi2Residuals;
+  if (mc)
+    *mc = fLastChi2MC;
+  if (mcLimit)
+    *mcLimit = fLastChi2MCLimit;
+  if (residuals)
+    *residuals = fLastChi2Residuals;
 }
 
 
@@ -956,7 +1023,7 @@ TH2F* AliMultiplicityCorrection::GetMultiplicityMC(Int_t i, EventType eventType)
 }
 
 //____________________________________________________________________
-void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Float_t regPar)
+void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Float_t regPar, Int_t nIterations, TH1* inputDist)
 {
   //
   // correct spectrum using bayesian method
@@ -966,52 +1033,59 @@ void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t ful
 
   SetupCurrentHists(inputRange, fullPhaseSpace, eventType);
 
-  // TODO should be taken from correlation map
-  TH1* sumHist = GetMultiplicityMC(inputRange, eventType)->ProjectionY("sumHist", 1, GetMultiplicityMC(inputRange, eventType)->GetNbinsX());
-
-  //new TCanvas;
-  //sumHist->Draw();
+  // normalize ESD
+  fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
 
   // normalize correction for given nPart
   for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
   {
-    //Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
-    Double_t sum = sumHist->GetBinContent(i);
-    if (sum <= 0)
-      continue;
+    // with this it is normalized to 1
+    Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
+
+    // with this normalized to the given efficiency
+    if (fCurrentEfficiency->GetBinContent(i) > 0)
+      sum /= fCurrentEfficiency->GetBinContent(i);
+    else
+      sum = 0;
+
     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 (sum > 0)
+      {
+        fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
+        fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
+      }
+      else
+      {
+        fCurrentCorrelation->SetBinContent(i, j, 0);
+        fCurrentCorrelation->SetBinError(i, j, 0);
+      }
     }
   }
 
-  //new TCanvas;
-  //fCurrentCorrelation->Draw("COLZ");
-  //return;
+  // smooth input spectrum
+  /*
+  TH1* esdClone = (TH1*) fCurrentESD->Clone("esdClone");
 
-  // normalize correction for given nTrack
-  /*for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
-  {
-    Double_t sum = fCurrentCorrelation->Integral(1, fCurrentCorrelation->GetNbinsX(), j, j);
-    if (sum <= 0)
-      continue;
-    for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsY(); ++i)
-    {
-      // ntrack sum to 1
-      fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
-      fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
-    }
-  }*/
+  for (Int_t i=2; i<fCurrentESD->GetNbinsX(); ++i)
+    if (esdClone->GetBinContent(i) < 1e-3)
+      fCurrentESD->SetBinContent(i, (esdClone->GetBinContent(i-1) + esdClone->GetBinContent(i) + esdClone->GetBinContent(i+1)) / 3);
 
-  // FAKE
-  fCurrentEfficiency = ((TH2*) fCurrentCorrelation)->ProjectionX("eff");
-  //new TCanvas;
-  //fCurrentEfficiency->Draw();
+  delete esdClone;
+
+  // rescale to 1
+  fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
+  */
 
   // pick prior distribution
-  TH1F* hPrior = (TH1F*)fCurrentESD->Clone("prior");
+  TH1* hPrior = 0;
+  if (inputDist)
+  {
+    printf("Using different starting conditions...\n");
+    hPrior = (TH1*)inputDist->Clone("prior");
+  }
+  else
+    hPrior = (TH1*)fCurrentESD->Clone("prior");
   Float_t norm = hPrior->Integral();
   for (Int_t t=1; t<=hPrior->GetNbinsX(); t++)
     hPrior->SetBinContent(t, hPrior->GetBinContent(t)/norm);
@@ -1027,9 +1101,10 @@ void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t ful
   TH2F* hInverseResponseBayes = (TH2F*)hResponse->Clone("pji");
   hInverseResponseBayes->Reset();
 
+  TH1F* convergence = new TH1F("convergence", "convergence", 50, 0.5, 50.5);
+
   // unfold...
-  Int_t iterations = 20;
-  for (Int_t i=0; i<iterations; i++)
+  for (Int_t i=0; i<nIterations; i++)
   {
     //printf(" iteration %i \n", i);
 
@@ -1056,17 +1131,10 @@ void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t ful
       for (Int_t m=1; m<=hResponse->GetNbinsY(); m++)
         value += fCurrentESD->GetBinContent(m) * hInverseResponseBayes->GetBinContent(t,m);
 
-      /*if (eventType == kTrVtx)
-      {
-        hTemp->SetBinContent(t, value);
-      }
-      else*/
-      {
-        if (fCurrentEfficiency->GetBinContent(t) > 0)
-          hTemp->SetBinContent(t, value / fCurrentEfficiency->GetBinContent(t));
-        else
-          hTemp->SetBinContent(t, 0);
-      }
+      if (fCurrentEfficiency->GetBinContent(t) > 0)
+        hTemp->SetBinContent(t, value / fCurrentEfficiency->GetBinContent(t));
+      else
+        hTemp->SetBinContent(t, 0);
     }
 
     // this is the last guess, fill before (!) smoothing
@@ -1083,10 +1151,10 @@ void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t ful
 
     for (Int_t t=2; t<hTrueSmoothed->GetNbinsX(); t++)
     {
-      Float_t average = (hTemp->GetBinContent(t-1) / hTemp->GetBinWidth(t-1)
-                         + hTemp->GetBinContent(t) / hTemp->GetBinWidth(t)
-                         + hTemp->GetBinContent(t+1) / hTemp->GetBinWidth(t+1)) / 3.;
-      average *= hTrueSmoothed->GetBinWidth(t);
+      Float_t average = (hTemp->GetBinContent(t-1)
+                         + hTemp->GetBinContent(t)
+                         + hTemp->GetBinContent(t+1)
+                         ) / 3.;
 
       // weight the average with the regularization parameter
       hTrueSmoothed->SetBinContent(t, (1-regPar) * hTemp->GetBinContent(t) + regPar * average);
@@ -1111,19 +1179,29 @@ void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t ful
 
     //printf("Chi2 of %d iteration = %.10f\n", i, chi2);
 
+    convergence->Fill(i+1, chi2);
+
     //if (i % 5 == 0)
     //  DrawComparison(Form("Bayesian_%d", i), mcTarget, correlationID, kTRUE, eventType);
 
     delete hTrueSmoothed;
   } // end of iterations
 
+  //new TCanvas;
+  //convergence->DrawCopy();
+  //gPad->SetLogy();
+
+  delete convergence;
 
   return;
 
   // ********
   // Calculate the covariance matrix, all arguments are taken from NIM,A362,487-498,1995
 
-  printf("Calculating covariance matrix. This may take some time...\n");
+  /*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();
@@ -1241,7 +1319,31 @@ void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t ful
       cov->SetBinContent(i, j, cov1[i-1][j-1] + cov2[i-1][j-1]);
 
   new TCanvas;
-  cov->Draw("COLZ");
+  cov->Draw("COLZ");*/
+}
+
+//____________________________________________________________________
+Float_t AliMultiplicityCorrection::BayesCovarianceDerivate(Float_t matrixM[251][251], TH2* hResponse,
+  TH1* fCurrentEfficiency, Int_t k, Int_t i, Int_t r, Int_t u)
+{
+  //
+  // helper function for the covariance matrix of the bayesian method
+  //
+
+  Float_t result = 0;
+
+  if (k == u && r == i)
+    result += 1.0 / hResponse->GetBinContent(u+1, r+1);
+
+  if (k == u)
+    result -= 1.0 / fCurrentEfficiency->GetBinContent(u+1);
+
+  if (r == i)
+    result -= matrixM[u][i] * fCurrentEfficiency->GetBinContent(u+1) / hResponse->GetBinContent(u+1, i+1);
+
+  result *= matrixM[k][i];
+
+  return result;
 }
 
 //____________________________________________________________________
@@ -1251,7 +1353,7 @@ void AliMultiplicityCorrection::ApplyLaszloMethod(Int_t inputRange, Bool_t fullP
   // correct spectrum using bayesian method
   //
 
-  Float_t regPar = 0.05;
+  Float_t regPar = 0;
 
   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
   Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4);
@@ -1299,7 +1401,7 @@ void AliMultiplicityCorrection::ApplyLaszloMethod(Int_t inputRange, Bool_t fullP
   TH2F* hResponse = (TH2F*) fCurrentCorrelation;
 
   // unfold...
-  Int_t iterations = 20;
+  Int_t iterations = 25;
   for (Int_t i=0; i<iterations; i++)
   {
     //printf(" iteration %i \n", i);
@@ -1348,7 +1450,7 @@ void AliMultiplicityCorrection::ApplyLaszloMethod(Int_t inputRange, Bool_t fullP
 
     for (Int_t t=1; t<hTrueSmoothed->GetNbinsX(); t++)
     {
-      Float_t newValue = hTrueSmoothed->GetBinContent(t)/norm;
+      Float_t newValue = hTemp->GetBinContent(t)/norm;
       Float_t diff = hPrior->GetBinContent(t) - newValue;
       chi2 += (Double_t) diff * diff;
 
@@ -1357,7 +1459,7 @@ void AliMultiplicityCorrection::ApplyLaszloMethod(Int_t inputRange, Bool_t fullP
 
     printf("Chi2 of %d iteration = %.10f\n", i, chi2);
 
-    if (i % 5 == 0)
+    //if (i % 5 == 0)
       DrawComparison(Form("Laszlo_%d", i), inputRange, fullPhaseSpace, kTRUE, GetMultiplicityMC(mcTarget, eventType)->ProjectionY());
 
     delete hTrueSmoothed;
@@ -1367,30 +1469,6 @@ void AliMultiplicityCorrection::ApplyLaszloMethod(Int_t inputRange, Bool_t fullP
 }
 
 //____________________________________________________________________
-Float_t AliMultiplicityCorrection::BayesCovarianceDerivate(Float_t matrixM[251][251], TH2* hResponse,
-  TH1* fCurrentEfficiency, Int_t k, Int_t i, Int_t r, Int_t u)
-{
-  //
-  //
-  //
-
-  Float_t result = 0;
-
-  if (k == u && r == i)
-    result += 1.0 / hResponse->GetBinContent(u+1, r+1);
-
-  if (k == u)
-    result -= 1.0 / fCurrentEfficiency->GetBinContent(u+1);
-
-  if (r == i)
-    result -= matrixM[u][i] * fCurrentEfficiency->GetBinContent(u+1) / hResponse->GetBinContent(u+1, i+1);
-
-  result *= matrixM[k][i];
-
-  return result;
-}
-
-//____________________________________________________________________
 void AliMultiplicityCorrection::ApplyGaussianMethod(Int_t inputRange, Bool_t fullPhaseSpace)
 {
   //
@@ -1556,7 +1634,7 @@ void AliMultiplicityCorrection::SetGenMeasFromFunc(TF1* inputMC, Int_t id)
   if (id < 0 || id >= kESDHists)
     return;
 
-  TH2F* mc = fMultiplicityINEL[id];
+  TH2F* mc = fMultiplicityVtx[id];
 
   mc->Reset();
 
@@ -1571,5 +1649,8 @@ void AliMultiplicityCorrection::SetGenMeasFromFunc(TF1* inputMC, Int_t id)
 
   TH1* proj = mc->ProjectionY();
 
+  TString tmp(fMultiplicityESD[id]->GetName());
+  delete fMultiplicityESD[id];
   fMultiplicityESD[id] = CalculateMultiplicityESD(proj, id);
+  fMultiplicityESD[id]->SetName(tmp);
 }
index 5ce8747..45c94f5 100644 (file)
@@ -24,7 +24,7 @@ class TCollection;
 class AliMultiplicityCorrection : public TNamed {
   public:
     enum EventType { kTrVtx = 0, kMB, kINEL };
-    enum RegularizationType { kNone = 0, kPol0, kPol1, kCurvature, kEntropy };
+    enum RegularizationType { kNone = 0, kPol0, kPol1, kCurvature, kEntropy, kTest };
 
     AliMultiplicityCorrection();
     AliMultiplicityCorrection(const Char_t* name, const Char_t* title);
@@ -42,11 +42,12 @@ class AliMultiplicityCorrection : public TNamed {
     void DrawHistograms();
     void DrawComparison(const char* name, Int_t inputRange, Bool_t fullPhaseSpace, Bool_t normalizeESD, TH1* mcHist);
 
-    void ApplyMinuitFit(Int_t inputRange, Bool_t fullPhaseSpace, Bool_t check = kFALSE);
-    void ApplyMinuitFitAll();
+    void ApplyMinuitFit(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t check = kFALSE, TH1* inputDist = 0);
     void SetRegularizationParameters(RegularizationType type, Float_t weight) { fRegularizationType = type; fRegularizationWeight = weight; };
 
-    void ApplyBayesianMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Float_t regPar = 0.07);
+    void ApplyNBDFit(Int_t inputRange, Bool_t fullPhaseSpace);
+
+    void ApplyBayesianMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Float_t regPar = 0.07, Int_t nIterations = 30, TH1* inputDist = 0);
 
     void ApplyGaussianMethod(Int_t inputRange, Bool_t fullPhaseSpace);
 
@@ -58,22 +59,22 @@ class AliMultiplicityCorrection : public TNamed {
     TH2F* GetMultiplicityINEL(Int_t i) { return fMultiplicityINEL[i]; }
     TH2F* GetMultiplicityMC(Int_t i, EventType eventType);
     TH3F* GetCorrelation(Int_t i) { return fCorrelation[i]; }
+    TH1F* GetMultiplicityESDCorrected(Int_t i) { return fMultiplicityESDCorrected[i]; }
 
     void SetMultiplicityESD(Int_t i, TH2F* hist)  { fMultiplicityESD[i] = hist; }
     void SetMultiplicityVtx(Int_t i, TH2F* hist)  { fMultiplicityVtx[i] = hist; }
     void SetMultiplicityMB(Int_t i, TH2F* hist)   { fMultiplicityMB[i] = hist; }
     void SetMultiplicityINEL(Int_t i, TH2F* hist) { fMultiplicityINEL[i] = hist; }
     void SetCorrelation(Int_t i, TH3F* hist) { fCorrelation[i] = hist; }
+    void SetMultiplicityESDCorrected(Int_t i, TH1F* hist) { fMultiplicityESDCorrected[i] = hist; }
 
     void SetGenMeasFromFunc(TF1* inputMC, Int_t id);
     TH2F* CalculateMultiplicityESD(TH1* inputMC, Int_t correlationMap, Bool_t normalized = kFALSE);
 
-    TH1F* GetMultiplicityESDCorrected(Int_t i) { return fMultiplicityESDCorrected[i]; }
-
     static void NormalizeToBinWidth(TH1* hist);
     static void NormalizeToBinWidth(TH2* hist);
 
-    void GetComparisonResults(Float_t& mc, Float_t& residuals);
+    void GetComparisonResults(Float_t* mc, Int_t* mcLimit, Float_t* residuals);
 
   protected:
     enum { kESDHists = 4, kMCHists = 5, kCorrHists = 8 };
@@ -84,8 +85,10 @@ class AliMultiplicityCorrection : public TNamed {
     static Double_t RegularizationPol1(Double_t *params);
     static Double_t RegularizationTotalCurvature(Double_t *params);
     static Double_t RegularizationEntropy(Double_t *params);
+    static Double_t RegularizationTest(Double_t *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 inputRange, Bool_t fullPhaseSpace, EventType eventType);
 
@@ -99,6 +102,8 @@ class AliMultiplicityCorrection : public TNamed {
     static TMatrixF* fCorrelationCovarianceMatrix;  // contains the errors of fCurrentESD
     static TVectorF* fCurrentESDVector;             // contains fCurrentESD
 
+    static TF1* fNBD;   // negative binomial distribution
+
     static RegularizationType fRegularizationType; // regularization that is used during Chi2 method
     static Float_t fRegularizationWeight;          // factor for regularization term
 
@@ -111,6 +116,7 @@ class AliMultiplicityCorrection : public TNamed {
     TH1F* fMultiplicityESDCorrected[kCorrHists]; // corrected histograms
 
     Float_t fLastChi2MC;        // last Chi2 between MC and unfolded ESD (calculated in DrawComparison)
+    Int_t   fLastChi2MCLimit;   // bin where the last chi2 breached a certain threshold, used to evaluate the multiplicity reach (calc. in DrawComparison)
     Float_t fLastChi2Residuals; // last Chi2 of the ESD and the folded unfolded ESD (calculated in DrawComparison)
 
  private:
index 4a57794..7f26a66 100644 (file)
@@ -2,23 +2,22 @@
 
 #include "AliMultiplicityESDSelector.h"
 
-#include <TStyle.h>
-#include <TSystem.h>
-#include <TCanvas.h>
 #include <TVector3.h>
-#include <TChain.h>
 #include <TFile.h>
-#include <TH1F.h>
+#include <TH2F.h>
+#include <TH3F.h>
+#include <TTree.h>
 
 #include <AliLog.h>
 #include <AliESD.h>
+#include <AliMultiplicity.h>
 
 #include "esdTrackCuts/AliESDtrackCuts.h"
 #include "AliPWG0Helper.h"
+#include "dNdEta/AliMultiplicityCorrection.h"
 
-#ifdef ALISELECTOR_USEMONALISA
-  #include <TMonaLisaWriter.h>
-#endif
+//#define TPCMEASUREMENT
+#define ITSMEASUREMENT
 
 ClassImp(AliMultiplicityESDSelector)
 
@@ -26,9 +25,6 @@ AliMultiplicityESDSelector::AliMultiplicityESDSelector() :
   AliSelector(),
   fMultiplicity(0),
   fEsdTrackCuts(0)
-#ifdef ALISELECTOR_USEMONALISA
-  ,fMonaLisaWriter(0)
-#endif
 {
   //
   // Constructor. Initialization of pointers
@@ -78,25 +74,30 @@ void AliMultiplicityESDSelector::SlaveBegin(TTree* tree)
 
   ReadUserObjects(tree);
 
-  fMultiplicity = new TH1F("fMultiplicity", "multiplicity", 201, 0.5, 200.5);
+  fMultiplicity = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
+}
 
-  #ifdef ALISELECTOR_USEMONALISA
-    TNamed *nm = 0;
-    if (fInput)
-      nm = dynamic_cast<TNamed*> (fInput->FindObject("PROOF_QueryTag"));
-    if (!nm)
-    {
-      AliDebug(AliLog::kError, "Query tag not found. Cannot enable monitoring");
-      return;
-    }
+void AliMultiplicityESDSelector::Init(TTree* tree)
+{
+  // read the user objects
 
-    TString option = GetOption();
-    option.ReplaceAll("#+", "");
+  AliSelector::Init(tree);
 
-    TString id;
-    id.Form("%s_%s%d", gSystem->HostName(), nm->GetTitle(), gSystem->GetPid());
-    fMonaLisaWriter = new TMonaLisaWriter(option, id, "CAF", "aliendb6.cern.ch");
-  #endif
+  // enable only the needed branches
+  if (tree)
+  {
+    tree->SetBranchStatus("*", 0);
+    tree->SetBranchStatus("fTriggerMask", 1);
+    tree->SetBranchStatus("fSPDVertex*", 1);
+
+    #ifdef ITSMEASUREMENT
+      tree->SetBranchStatus("fSPDMult*", 1);
+    #endif
+
+    #ifdef TPCMEASUREMENT
+      AliESDtrackCuts::EnableNeededBranches(tree);
+    #endif
+  }
 }
 
 Bool_t AliMultiplicityESDSelector::Process(Long64_t entry)
@@ -129,22 +130,103 @@ Bool_t AliMultiplicityESDSelector::Process(Long64_t entry)
     return kFALSE;
   }
 
+  Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD);
+  Bool_t eventVertex = AliPWG0Helper::IsVertexReconstructed(fESD);
+
+  if (!eventTriggered || !eventVertex)
+    return kTRUE;
+
+  // get the ESD vertex
+  const AliESDVertex* vtxESD = fESD->GetVertex();
+  Double_t vtx[3];
+  vtxESD->GetXYZ(vtx);
+
+  Int_t nESDTracks05 = 0;
+  Int_t nESDTracks10 = 0;
+  Int_t nESDTracks15 = 0;
+  Int_t nESDTracks20 = 0;
+
+#ifdef ITSMEASUREMENT
+  // get tracklets
+  const AliMultiplicity* mult = fESD->GetMultiplicity();
+  if (!mult)
+  {
+    AliDebug(AliLog::kError, "AliMultiplicity not available");
+    return kFALSE;
+  }
+
+  // get multiplicity from ITS tracklets
+  for (Int_t i=0; i<mult->GetNumberOfTracklets(); ++i)
+  {
+    //printf("%d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), mult->GetDeltaPhi(i));
+
+    // this removes non-tracklets. Very bad solution. SPD guys are working on better solution...
+    if (mult->GetDeltaPhi(i) < -1000)
+      continue;
+
+    Float_t theta = mult->GetTheta(i);
+    Float_t eta   = -TMath::Log(TMath::Tan(theta/2.));
+
+    if (TMath::Abs(eta) < 0.5)
+      nESDTracks05++;
+
+    if (TMath::Abs(eta) < 1.0)
+      nESDTracks10++;
+
+    if (TMath::Abs(eta) < 1.5)
+      nESDTracks15++;
+
+    if (TMath::Abs(eta) < 2.0)
+      nESDTracks20++;
+  }
+#endif
+
+#ifdef TPCMEASUREMENT
   if (!fEsdTrackCuts)
   {
     AliDebug(AliLog::kError, "fESDTrackCuts not available");
     return kFALSE;
   }
 
-  if (AliPWG0Helper::IsEventTriggered(fESD) == kFALSE)
-    return kTRUE;
+  // get multiplicity from ESD tracks
+  TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD);
+  Int_t nGoodTracks = list->GetEntries();
+  // loop over esd tracks
+  for (Int_t i=0; i<nGoodTracks; i++)
+  {
+    AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (list->At(i));
+    if (!esdTrack)
+    {
+      AliDebug(AliLog::kError, Form("ERROR: Could not retrieve track %d.", i));
+      continue;
+    }
 
-  if (AliPWG0Helper::IsVertexReconstructed(fESD) == kFALSE)
-    return kTRUE;
+    Double_t p[3];
+    esdTrack->GetConstrainedPxPyPz(p); // ### TODO should be okay because we have a vertex, however GetInnerPxPyPy / GetOuterPxPyPy also exist
+    TVector3 vector(p);
+
+    Float_t theta = vector.Theta();
+    Float_t eta   = -TMath::Log(TMath::Tan(theta/2.));
+    Float_t pt = vector.Pt();
+
+    //if (pt < kPtCut)
+    //  continue;
 
-  // get number of "good" tracks
-  Int_t nGoodTracks = fEsdTrackCuts->CountAcceptedTracks(fESD);
+    if (TMath::Abs(eta) < 0.5)
+      nESDTracks05++;
 
-  fMultiplicity->Fill(nGoodTracks);
+    if (TMath::Abs(eta) < 1.0)
+      nESDTracks10++;
+
+    if (TMath::Abs(eta) < 1.5)
+      nESDTracks15++;
+
+    if (TMath::Abs(eta) < 2.0)
+      nESDTracks20++;
+  }
+#endif
+
+  fMultiplicity->FillMeasured(vtx[2], nESDTracks05, nESDTracks10, nESDTracks15, nESDTracks20);
 
   return kTRUE;
 }
@@ -157,14 +239,6 @@ void AliMultiplicityESDSelector::SlaveTerminate()
 
   AliSelector::SlaveTerminate();
 
-  #ifdef ALISELECTOR_USEMONALISA
-    if (fMonaLisaWriter)
-    {
-      delete fMonaLisaWriter;
-      fMonaLisaWriter = 0;
-    }
-  #endif
-
   // Add the histograms to the output on each slave server
   if (!fOutput)
   {
@@ -183,15 +257,19 @@ void AliMultiplicityESDSelector::Terminate()
 
   AliSelector::Terminate();
 
-  fMultiplicity = dynamic_cast<TH1F*> (fOutput->FindObject("fMultiplicity"));
+  fMultiplicity = dynamic_cast<AliMultiplicityCorrection*> (fOutput->FindObject("Multiplicity"));
 
   if (!fMultiplicity)
   {
-    AliDebug(AliLog::kError, Form("ERROR: Histogram not available %p", (void*) fMultiplicity));
+    AliDebug(AliLog::kError, Form("ERROR: Histograms not available %p", (void*) fMultiplicity));
     return;
   }
 
   TFile* file = TFile::Open("multiplicityESD.root", "RECREATE");
-  fMultiplicity->Write();
+
+  fMultiplicity->SaveHistograms();
+
   file->Close();
+
+  fMultiplicity->DrawHistograms();
 }
index f818ad3..6487fef 100644 (file)
@@ -5,15 +5,8 @@
 
 #include "AliSelector.h"
 
-// uncomment this to enable mona lisa monitoring
-//#define ALISELECTOR_USEMONALISA
-
 class AliESDtrackCuts;
-class TH1F;
-
-#ifdef ALISELECTOR_USEMONALISA
-  class TMonaLisaWriter;
-#endif
+class AliMultiplicityCorrection;
 
 class AliMultiplicityESDSelector : public AliSelector {
   public:
@@ -22,6 +15,7 @@ class AliMultiplicityESDSelector : public AliSelector {
 
     virtual void    Begin(TTree* tree);
     virtual void    SlaveBegin(TTree *tree);
+    virtual void    Init(TTree *tree);
     virtual Bool_t  Process(Long64_t entry);
     virtual void    SlaveTerminate();
     virtual void    Terminate();
@@ -29,18 +23,13 @@ class AliMultiplicityESDSelector : public AliSelector {
  protected:
     void ReadUserObjects(TTree* tree);
 
-    TH1F* fMultiplicity; // multiplicity histogram
-
-    AliESDtrackCuts*  fEsdTrackCuts;     // Object containing the parameters of the esd track cuts
+    AliMultiplicityCorrection* fMultiplicity; // object containing the extracted data
+    AliESDtrackCuts* fEsdTrackCuts;           // Object containing the parameters of the esd track cuts
 
  private:
     AliMultiplicityESDSelector(const AliMultiplicityESDSelector&);
     AliMultiplicityESDSelector& operator=(const AliMultiplicityESDSelector&);
 
-#ifdef ALISELECTOR_USEMONALISA
-    TMonaLisaWriter* fMonaLisaWriter; //! ML instance for monitoring
-#endif
-
   ClassDef(AliMultiplicityESDSelector, 0);
 };
 
index aae710d..dbecc09 100644 (file)
@@ -96,7 +96,10 @@ void AliMultiplicityMCSelector::Init(TTree* tree)
     tree->SetBranchStatus("*", 0);
     tree->SetBranchStatus("fTriggerMask", 1);
     tree->SetBranchStatus("fSPDVertex*", 1);
-    tree->SetBranchStatus("fSPDMult*", 1);
+
+    #ifdef ITSMEASUREMENT
+      tree->SetBranchStatus("fSPDMult*", 1);
+    #endif
 
     #ifdef TPCMEASUREMENT
       AliESDtrackCuts::EnableNeededBranches(tree);
@@ -298,7 +301,7 @@ Bool_t AliMultiplicityMCSelector::Process(Long64_t entry)
   }
 #endif
 
-  fMultiplicity->FillMeasured(vtxMC[2], nESDTracks05, nESDTracks10, nESDTracks15, nESDTracks20);
+  fMultiplicity->FillMeasured(vtx[2], nESDTracks05, nESDTracks10, nESDTracks15, nESDTracks20);
 
   // TODO should this be vtx[2] or vtxMC[2] ?
   fMultiplicity->FillCorrection(vtxMC[2], nMCTracks05, nMCTracks10, nMCTracks15, nMCTracks20, nMCTracksAll, nESDTracks05, nESDTracks10, nESDTracks15, nESDTracks20);
index 675643b..3dfec89 100644 (file)
@@ -6,8 +6,6 @@
 #include "AliSelectorRL.h"
 
 class AliESDtrackCuts;
-class TH2F;
-class TH3F;
 class AliMultiplicityCorrection;
 
 class AliMultiplicityMCSelector : public AliSelectorRL {
index 33109a4..9dd5e60 100644 (file)
@@ -24,10 +24,6 @@ void runMultiplicitySelector(Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool_
   if (!prepareQuery(libraries, packages, kTRUE))
     return;
 
-  gProof->SetParallel(0);
-  gProof->SetLogLevel(4);
-  gProof->SetParallel(9999);
-
   gROOT->ProcessLine(".L CreateCuts.C");
   gROOT->ProcessLine(".L drawPlots.C");
 
@@ -82,18 +78,39 @@ void* fit(const char* fileName = "multiplicityMC.root", Int_t hist = 2, Int_t ev
   TFile::Open(fileName);
   mult->LoadHistograms("Multiplicity");
 
-  mult->ApplyBayesianMethod(hist, kFALSE, eventType);
-  mult->DrawComparison("Bayesian", hist, kFALSE, kTRUE, mult->GetMultiplicityMC(hist, eventType)->ProjectionY());
+  //mult->ApplyLaszloMethod(hist, kFALSE, AliMultiplicityCorrection::kTrVtx);
+
+  //return;
+
+
+  mult->ApplyBayesianMethod(hist, kFALSE, AliMultiplicityCorrection::kTrVtx);
+  mult->DrawComparison("Bayesian", hist, kFALSE, kTRUE, mult->GetMultiplicityMC(hist, AliMultiplicityCorrection::kTrVtx)->ProjectionY());
+
+  return;
+
+  TStopwatch timer;
 
-  //mult->ApplyMinuitFit(hist, kFALSE);
-  //mult->DrawComparison("MinuitChi2", hist, kFALSE, kTRUE, mult->GetMultiplicityMC(hist, AliMultiplicityCorrection::kTrVtx)->ProjectionY());
+  timer.Start();
+
+  mult->SetRegularizationParameters(AliMultiplicityCorrection::kTest, 1);
+  mult->ApplyMinuitFit(hist, kFALSE, AliMultiplicityCorrection::kTrVtx);
+  mult->DrawComparison("MinuitChi2", hist, kFALSE, kFALSE, mult->GetMultiplicityMC(hist, AliMultiplicityCorrection::kTrVtx)->ProjectionY());
+
+  timer.Stop();
+  timer.Print();
+
+  return 0;
 
   //mult->ApplyGaussianMethod(hist, kFALSE);
 
+  mult->SetRegularizationParameters(AliMultiplicityCorrection::kNone, 0);
+  mult->ApplyNBDFit(hist, kFALSE);
+  mult->DrawComparison("NBDChi2Fit", hist, kFALSE, kTRUE, mult->GetMultiplicityMC(hist, eventType)->ProjectionY());
+
   return mult;
 }
 
-void* fitOther(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityESD.root", Int_t histID = 2)
+void* fitOther(const char* fileNameMC = "multiplicityMC_3M.root", const char* fileNameESD = "multiplicityMC_3M.root", Int_t histID = 2)
 {
   gSystem->Load("libPWG0base");
 
@@ -105,16 +122,32 @@ void* fitOther(const char* fileNameMC = "multiplicityMC.root", const char* fileN
   TFile::Open(fileNameESD);
   TH2F* hist = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityESD%d", histID));
   TH2F* hist2 = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityVtx%d", histID));
+  //hist2 = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityMB%d", histID));
 
   mult->SetMultiplicityESD(histID, hist);
 
+  mult->SetMultiplicityVtx(histID, hist2);
+  mult->ApplyLaszloMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx);
+  return;
+
+  mult->SetRegularizationParameters(AliMultiplicityCorrection::kTest, 1.1);
+  mult->ApplyMinuitFit(histID, kFALSE, AliMultiplicityCorrection::kTrVtx);
+  mult->DrawComparison("MinuitChi2", histID, kFALSE, kTRUE, hist2->ProjectionY());
+
+  return;
+
   //mult->ApplyGaussianMethod(histID, kFALSE);
-  //for (Float_t f=0; f<0.1; f+=0.01)
-  //mult->ApplyBayesianMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx);
-  mult->ApplyMinuitFit(histID, kFALSE);
-  mult->DrawComparison("MinuitChi2", hist, kFALSE, kTRUE, hist2->ProjectionY());
 
-  //mult->ApplyLaszloMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx);
+  for (Float_t f=0.1; f<=0.11; f+=0.05)
+  {
+    mult->ApplyBayesianMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, f);
+    mult->DrawComparison(Form("Bayesian_%f", f), histID, kFALSE, kTRUE, hist2->ProjectionY());
+  }
+
+  //mult->SetRegularizationParameters(AliMultiplicityCorrection::kEntropy, 1e7);
+  //mult->ApplyMinuitFit(histID, kFALSE);
+  //mult->DrawComparison("MinuitChi2", histID, kFALSE, kTRUE, hist2->ProjectionY());
+
 
   return mult;
 }
@@ -128,12 +161,188 @@ const char* GetRegName(Int_t type)
     case AliMultiplicityCorrection::kPol1:      return "Pol1"; break;
     case AliMultiplicityCorrection::kCurvature: return "TotalCurvature"; break;
     case AliMultiplicityCorrection::kEntropy:   return "Reduced cross-entropy"; break;
+    case AliMultiplicityCorrection::kTest   :   return "Test"; break;
+  }
+  return 0;
+}
+
+const char* GetEventTypeName(Int_t type)
+{
+  switch (type)
+  {
+    case AliMultiplicityCorrection::kTrVtx:   return "trigger, vertex"; break;
+    case AliMultiplicityCorrection::kMB:      return "minimum bias"; break;
+    case AliMultiplicityCorrection::kINEL:    return "inelastic"; break;
   }
   return 0;
 }
 
-void EvaluateChi2Method(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityMC.root", Int_t histID = 2)
+void EvaluateBayesianMethod(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityMC.root", const char* targetDir, Int_t histID = 2)
 {
+  gSystem->mkdir(targetDir);
+
+  gSystem->Load("libPWG0base");
+
+  AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
+  TFile::Open(fileNameMC);
+  mult->LoadHistograms("Multiplicity");
+
+  AliMultiplicityCorrection* multESD = new AliMultiplicityCorrection("MultiplicityESD", "MultiplicityESD");
+  TFile::Open(fileNameESD);
+  multESD->LoadHistograms("Multiplicity");
+  mult->SetMultiplicityESD(histID, multESD->GetMultiplicityESD(histID));
+
+  TCanvas* canvas = new TCanvas("EvaluateBayesianMethod", "EvaluateBayesianMethod", 800, 600);
+  TLegend* legend = new TLegend(0.2, 0.7, 0.58, 0.98);
+  legend->SetFillColor(0);
+
+  Float_t min = 1e10;
+  Float_t max = 0;
+
+  TGraph* first = 0;
+  Int_t count = 0; // just to order the saved images...
+
+  for (AliMultiplicityCorrection::EventType type = AliMultiplicityCorrection::kTrVtx; type <= AliMultiplicityCorrection::kINEL; ++type)
+  {
+    TGraph* fitResultsMC = new TGraph;
+    fitResultsMC->SetTitle(";Weight Parameter");
+    TGraph* fitResultsRes = new TGraph;
+    fitResultsRes->SetTitle(";Weight Parameter");
+
+    fitResultsMC->SetFillColor(0);
+    fitResultsRes->SetFillColor(0);
+    fitResultsMC->SetMarkerStyle(20+type);
+    fitResultsRes->SetMarkerStyle(24+type);
+    fitResultsRes->SetMarkerColor(kRed);
+    fitResultsRes->SetLineColor(kRed);
+
+    legend->AddEntry(fitResultsMC, Form("%s MC chi2", GetEventTypeName(type)));
+    legend->AddEntry(fitResultsRes, Form("%s residual chi2", GetEventTypeName(type)));
+
+    for (Float_t weight = 0; weight < 0.301; weight += 0.02)
+    {
+      Float_t chi2MC = 0;
+      Float_t residuals = 0;
+
+      mult->ApplyBayesianMethod(histID, kFALSE, type, weight);
+      mult->DrawComparison(Form("%s/Bayesian_%02d_%d_%f", targetDir, count++, type, weight), histID, kFALSE, kTRUE, multESD->GetMultiplicityMC(histID, type)->ProjectionY());
+      mult->GetComparisonResults(&chi2MC, 0, &residuals);
+
+      fitResultsMC->SetPoint(fitResultsMC->GetN(), weight, chi2MC);
+      fitResultsRes->SetPoint(fitResultsRes->GetN(), weight, residuals);
+
+      min = TMath::Min(min, TMath::Min(chi2MC, residuals));
+      max = TMath::Max(max, TMath::Max(chi2MC, residuals));
+    }
+
+    fitResultsMC->Print();
+    fitResultsRes->Print();
+
+    canvas->cd();
+    fitResultsMC->Draw(Form("%s CP", (first == 0) ? "A" : "SAME"));
+    fitResultsRes->Draw("SAME CP");
+
+    if (first == 0)
+      first = fitResultsMC;
+  }
+
+  gPad->SetLogy();
+  printf("min = %f, max = %f\n", min, max);
+  if (min <= 0)
+    min = 1e-5;
+  first->GetYaxis()->SetRangeUser(min * 0.5, max * 1.5);
+
+  legend->Draw();
+
+  canvas->SaveAs(Form("%s/%s.gif", targetDir, canvas->GetName()));
+  canvas->SaveAs(Form("%s/%s.C", targetDir, canvas->GetName()));
+}
+
+void EvaluateBayesianMethodIterations(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityMC.root", const char* targetDir, Int_t histID = 2)
+{
+  gSystem->mkdir(targetDir);
+
+  gSystem->Load("libPWG0base");
+
+  AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
+  TFile::Open(fileNameMC);
+  mult->LoadHistograms("Multiplicity");
+
+  AliMultiplicityCorrection* multESD = new AliMultiplicityCorrection("MultiplicityESD", "MultiplicityESD");
+  TFile::Open(fileNameESD);
+  multESD->LoadHistograms("Multiplicity");
+  mult->SetMultiplicityESD(histID, multESD->GetMultiplicityESD(histID));
+
+  TCanvas* canvas = new TCanvas("EvaluateBayesianMethodIterations", "EvaluateBayesianMethodIterations", 800, 600);
+  TLegend* legend = new TLegend(0.2, 0.7, 0.58, 0.98);
+  legend->SetFillColor(0);
+
+  Float_t min = 1e10;
+  Float_t max = 0;
+
+  TGraph* first = 0;
+  Int_t count = 0; // just to order the saved images...
+
+  for (AliMultiplicityCorrection::EventType type = AliMultiplicityCorrection::kTrVtx; type <= AliMultiplicityCorrection::kINEL; ++type)
+  {
+    TGraph* fitResultsMC = new TGraph;
+    fitResultsMC->SetTitle(";Iterations");
+    TGraph* fitResultsRes = new TGraph;
+    fitResultsRes->SetTitle(";Iterations");
+
+    fitResultsMC->SetFillColor(0);
+    fitResultsRes->SetFillColor(0);
+    fitResultsMC->SetMarkerStyle(20+type);
+    fitResultsRes->SetMarkerStyle(24+type);
+    fitResultsRes->SetMarkerColor(kRed);
+    fitResultsRes->SetLineColor(kRed);
+
+    legend->AddEntry(fitResultsMC, Form("%s MC chi2", GetEventTypeName(type)));
+    legend->AddEntry(fitResultsRes, Form("%s residual chi2", GetEventTypeName(type)));
+
+    for (Int_t iter = 5; iter <= 50; iter += 5)
+    {
+      Float_t chi2MC = 0;
+      Float_t residuals = 0;
+
+      mult->ApplyBayesianMethod(histID, kFALSE, type, 0.1, iter);
+      mult->DrawComparison(Form("%s/BayesianIter_%02d_%d_%d", targetDir, count++, type, iter), histID, kFALSE, kTRUE, multESD->GetMultiplicityMC(histID, type)->ProjectionY());
+      mult->GetComparisonResults(&chi2MC, 0, &residuals);
+
+      fitResultsMC->SetPoint(fitResultsMC->GetN(), iter, chi2MC);
+      fitResultsRes->SetPoint(fitResultsRes->GetN(), iter, residuals);
+
+      min = TMath::Min(min, TMath::Min(chi2MC, residuals));
+      max = TMath::Max(max, TMath::Max(chi2MC, residuals));
+    }
+
+    fitResultsMC->Print();
+    fitResultsRes->Print();
+
+    canvas->cd();
+    fitResultsMC->Draw(Form("%s CP", (first == 0) ? "A" : "SAME"));
+    fitResultsRes->Draw("SAME CP");
+
+    if (first == 0)
+      first = fitResultsMC;
+  }
+
+  gPad->SetLogy();
+  printf("min = %f, max = %f\n", min, max);
+  if (min <= 0)
+    min = 1e-5;
+  first->GetYaxis()->SetRangeUser(min * 0.5, max * 1.5);
+
+  legend->Draw();
+
+  canvas->SaveAs(Form("%s/%s.gif", targetDir, canvas->GetName()));
+  canvas->SaveAs(Form("%s/%s.C", targetDir, canvas->GetName()));
+}
+
+void EvaluateChi2Method(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityMC.root", const char* targetDir, Int_t histID = 2)
+{
+  gSystem->mkdir(targetDir);
+
   gSystem->Load("libPWG0base");
 
   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
@@ -148,15 +357,19 @@ void EvaluateChi2Method(const char* fileNameMC = "multiplicityMC.root", const ch
   mult->SetMultiplicityESD(histID, hist);
 
   TCanvas* canvas = new TCanvas("EvaluateChi2Method", "EvaluateChi2Method", 800, 600);
-  TLegend* legend = new TLegend(0.6, 0.1, 0.98, 0.4);
+  TLegend* legend = new TLegend(0.2, 0.7, 0.58, 0.98);
   legend->SetFillColor(0);
 
   Float_t min = 1e10;
   Float_t max = 0;
 
   TGraph* first = 0;
+  Int_t count = 0; // just to order the saved images...
+
+  Bool_t firstLoop = kTRUE;
 
   for (AliMultiplicityCorrection::RegularizationType type = AliMultiplicityCorrection::kPol0; type <= AliMultiplicityCorrection::kEntropy; ++type)
+  //for (AliMultiplicityCorrection::RegularizationType type = AliMultiplicityCorrection::kPol1; type <= AliMultiplicityCorrection::kPol1; ++type)
   {
     TGraph* fitResultsMC = new TGraph;
     fitResultsMC->SetTitle(";Weight Parameter");
@@ -166,7 +379,7 @@ void EvaluateChi2Method(const char* fileNameMC = "multiplicityMC.root", const ch
     fitResultsMC->SetFillColor(0);
     fitResultsRes->SetFillColor(0);
     fitResultsMC->SetMarkerStyle(19+type);
-    fitResultsRes->SetMarkerStyle(19+type);
+    fitResultsRes->SetMarkerStyle(23+type);
     fitResultsRes->SetMarkerColor(kRed);
     fitResultsRes->SetLineColor(kRed);
 
@@ -176,15 +389,16 @@ void EvaluateChi2Method(const char* fileNameMC = "multiplicityMC.root", const ch
     if (first == 0)
       first = fitResultsMC;
 
-    for (Float_t weight = 1e-2; weight < 2e4; weight *= 1e2)
+    for (Float_t weight = 1e-4; weight < 2e4; weight *= 10)
+    //for (Float_t weight = 0.1; weight < 10; weight *= TMath::Sqrt(TMath::Sqrt(10)))
     {
       Float_t chi2MC = 0;
       Float_t residuals = 0;
 
       mult->SetRegularizationParameters(type, weight);
-      mult->ApplyMinuitFit(histID, kFALSE);
-      mult->DrawComparison(Form("MinuitChi2_%d_%f", type, weight), histID, kFALSE, kTRUE, hist2->ProjectionY());
-      mult->GetComparisonResults(chi2MC, residuals);
+      mult->ApplyMinuitFit(histID, kFALSE, AliMultiplicityCorrection::kTrVtx);
+      mult->DrawComparison(Form("%s/MinuitChi2_%02d_%d_%f", targetDir, count++, type, weight), histID, kFALSE, kTRUE, hist2->ProjectionY());
+      mult->GetComparisonResults(&chi2MC, 0, &residuals);
 
       fitResultsMC->SetPoint(fitResultsMC->GetN(), weight, chi2MC);
       fitResultsRes->SetPoint(fitResultsRes->GetN(), weight, residuals);
@@ -197,8 +411,10 @@ void EvaluateChi2Method(const char* fileNameMC = "multiplicityMC.root", const ch
     fitResultsRes->Print();
 
     canvas->cd();
-    fitResultsMC->Draw(Form("%s CP", (type == AliMultiplicityCorrection::kPol0) ? "A" : "SAME"));
+    fitResultsMC->Draw(Form("%s CP", (firstLoop) ? "A" : "SAME"));
     fitResultsRes->Draw("SAME CP");
+
+    firstLoop = kFALSE;
   }
 
   gPad->SetLogx();
@@ -210,7 +426,416 @@ void EvaluateChi2Method(const char* fileNameMC = "multiplicityMC.root", const ch
 
   legend->Draw();
 
+  canvas->SaveAs(Form("%s/%s.gif", targetDir, canvas->GetName()));
+  canvas->SaveAs(Form("%s/%s.C", targetDir, canvas->GetName()));
+}
+
+void EvaluateChi2MethodAll()
+{
+  EvaluateChi2Method("multiplicityMC_3M.root", "multiplicityMC_3M.root", "eval-3M-3M");
+  EvaluateChi2Method("multiplicityMC_2M.root", "multiplicityMC_1M_3.root", "eval-2M-1M");
+  EvaluateChi2Method("multiplicityMC_3M.root", "multiplicityMC_3M_NBD.root", "eval-3M-NBD");
+  EvaluateChi2Method("multiplicityMC_2M_smoothed.root", "multiplicityMC_1M_3.root", "eval-2MS-1M");
+  EvaluateChi2Method("multiplicityMC_2M_smoothed.root", "multiplicityMC_3M_NBD.root", "eval-2MS-NBD");
+}
+
+void EvaluateBayesianMethodAll()
+{
+  EvaluateBayesianMethod("multiplicityMC_3M.root", "multiplicityMC_3M.root", "eval-3M-3M");
+  EvaluateBayesianMethod("multiplicityMC_2M.root", "multiplicityMC_1M_3.root", "eval-2M-1M");
+  EvaluateBayesianMethod("multiplicityMC_3M.root", "multiplicityMC_3M_NBD.root", "eval-3M-NBD");
+  EvaluateBayesianMethod("multiplicityMC_2M_smoothed.root", "multiplicityMC_1M_3.root", "eval-2MS-1M");
+  EvaluateBayesianMethod("multiplicityMC_2M_smoothed.root", "multiplicityMC_3M_NBD.root", "eval-2MS-NBD");
+}
+
+void CompareMethods(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityMC.root", const char* targetDir, Int_t histID = 2)
+{
+  gSystem->mkdir(targetDir);
+
+  gSystem->Load("libPWG0base");
+
+  AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
+
+  TFile::Open(fileNameMC);
+  mult->LoadHistograms("Multiplicity");
+
+  TFile::Open(fileNameESD);
+  AliMultiplicityCorrection* multESD = new AliMultiplicityCorrection("MultiplicityESD", "MultiplicityESD");
+  multESD->LoadHistograms("Multiplicity");
+
+  mult->SetMultiplicityESD(histID, multESD->GetMultiplicityESD(histID));
+
+  TCanvas* canvas = new TCanvas("CompareMethods", "CompareMethods", 1200, 800);
+  canvas->Divide(3, 2);
+
+  Int_t count = 0;
+
+  for (AliMultiplicityCorrection::EventType type = AliMultiplicityCorrection::kTrVtx; type <= AliMultiplicityCorrection::kINEL; ++type)
+  {
+    TH1* mc = multESD->GetMultiplicityMC(histID, type)->ProjectionY();
+
+    mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 0.3);
+    mult->ApplyMinuitFit(histID, kFALSE, type);
+    mult->DrawComparison(Form("%s/CompareMethods_%d_MinuitChi2", targetDir, type), histID, kFALSE, kTRUE, mc);
+    TH1* chi2Result = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone("chi2Result");
+
+    mult->ApplyBayesianMethod(histID, kFALSE, type, 0.1);
+    mult->DrawComparison(Form("%s/CompareMethods_%d_Bayesian", targetDir, type), histID, kFALSE, kTRUE, mc);
+    TH1* bayesResult = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone("bayesResult");
+
+    mc->GetXaxis()->SetRangeUser(0, 150);
+    chi2Result->GetXaxis()->SetRangeUser(0, 150);
+
+    // skip errors for now
+    for (Int_t i=1; i<=chi2Result->GetNbinsX(); ++i)
+    {
+      chi2Result->SetBinError(i, 0);
+      bayesResult->SetBinError(i, 0);
+    }
+
+    canvas->cd(++count);
+    mc->SetFillColor(kYellow);
+    mc->DrawCopy();
+    chi2Result->SetLineColor(kRed);
+    chi2Result->DrawCopy("SAME");
+    bayesResult->SetLineColor(kBlue);
+    bayesResult->DrawCopy("SAME");
+    gPad->SetLogy();
+
+    canvas->cd(count + 3);
+    chi2Result->Divide(chi2Result, mc, 1, 1, "B");
+    bayesResult->Divide(bayesResult, mc, 1, 1, "B");
+
+    // skip errors for now
+    for (Int_t i=1; i<=chi2Result->GetNbinsX(); ++i)
+    {
+      chi2Result->SetBinError(i, 0);
+      bayesResult->SetBinError(i, 0);
+    }
+
+    chi2Result->GetYaxis()->SetRangeUser(0.5, 1.5);
+
+    chi2Result->DrawCopy("");
+    bayesResult->DrawCopy("SAME");
+  }
+
+  canvas->SaveAs(Form("%s/%s.gif", targetDir, canvas->GetName()));
+  canvas->SaveAs(Form("%s/%s.C", targetDir, canvas->GetName()));
+}
+
+void StatisticsPlot(const char* fileNameMC = "multiplicityMC.root", Int_t histID = 2)
+{
+  gSystem->Load("libPWG0base");
+
+  AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
+
+  TFile::Open(fileNameMC);
+  mult->LoadHistograms("Multiplicity");
+
+  const char* files[] = { "multiplicityMC_0.5M.root", "multiplicityMC_0.75M.root", "multiplicityMC_1M_3.root", "multiplicityMC_1.25M.root", "multiplicityMC_1.5M.root" };
+
+  TGraph* fitResultsChi2 = new TGraph;
+  fitResultsChi2->SetTitle(";Nevents;Chi2");
+  TGraph* fitResultsBayes = new TGraph;
+  fitResultsBayes->SetTitle(";Nevents;Chi2");
+  TGraph* fitResultsChi2Limit = new TGraph;
+  fitResultsChi2Limit->SetTitle(";Nevents;Multiplicity reach");
+  TGraph* fitResultsBayesLimit = new TGraph;
+  fitResultsBayesLimit->SetTitle(";Nevents;Multiplicity reach");
+
+  TCanvas* canvas = new TCanvas("StatisticsPlot", "StatisticsPlot", 1200, 600);
+  canvas->Divide(5, 2);
+
+  Float_t min = 1e10;
+  Float_t max = 0;
+
+  for (Int_t i=0; i<5; ++i)
+  {
+    TFile::Open(files[i]);
+    AliMultiplicityCorrection* multESD = new AliMultiplicityCorrection("MultiplicityESD", "MultiplicityESD");
+    multESD->LoadHistograms("Multiplicity");
+
+    mult->SetMultiplicityESD(histID, multESD->GetMultiplicityESD(histID));
+    Int_t nEntries = multESD->GetMultiplicityESD(histID)->GetEntries();
+    TH1* mc = multESD->GetMultiplicityMC(histID, AliMultiplicityCorrection::kTrVtx)->ProjectionY();
+
+    mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 0.3);
+    mult->ApplyMinuitFit(histID, kFALSE, AliMultiplicityCorrection::kTrVtx);
+    mult->DrawComparison(Form("StatisticsPlot_%d_MinuitChi2", i), histID, kFALSE, kTRUE, mc);
+
+    Float_t chi2MC = 0;
+    Int_t chi2MCLimit = 0;
+    mult->GetComparisonResults(&chi2MC, &chi2MCLimit, 0);
+    fitResultsChi2->SetPoint(fitResultsChi2->GetN(), nEntries, chi2MC);
+    fitResultsChi2Limit->SetPoint(fitResultsChi2Limit->GetN(), nEntries, chi2MCLimit);
+    min = TMath::Min(min, chi2MC);
+    max = TMath::Max(max, chi2MC);
+
+    TH1* chi2Result = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone("chi2Result");
+
+    mult->ApplyBayesianMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, 0.1);
+    mult->DrawComparison(Form("StatisticsPlot_%d_Bayesian", i), histID, kFALSE, kTRUE, mc);
+    TH1* bayesResult = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone("bayesResult");
+    mult->GetComparisonResults(&chi2MC, &chi2MCLimit, 0);
+    fitResultsBayes->SetPoint(fitResultsBayes->GetN(), nEntries, chi2MC);
+    fitResultsBayesLimit->SetPoint(fitResultsBayesLimit->GetN(), nEntries, chi2MCLimit);
+
+    min = TMath::Min(min, chi2MC);
+    max = TMath::Max(max, chi2MC);
+    mc->GetXaxis()->SetRangeUser(0, 150);
+    chi2Result->GetXaxis()->SetRangeUser(0, 150);
+
+    // skip errors for now
+    for (Int_t j=0; j<=chi2Result->GetNbinsX(); ++j)
+    {
+      chi2Result->SetBinError(j, 0);
+      bayesResult->SetBinError(j, 0);
+    }
+
+    canvas->cd(i+1);
+    mc->SetFillColor(kYellow);
+    mc->DrawCopy();
+    chi2Result->SetLineColor(kRed);
+    chi2Result->DrawCopy("SAME");
+    bayesResult->SetLineColor(kBlue);
+    bayesResult->DrawCopy("SAME");
+    gPad->SetLogy();
+
+    canvas->cd(i+6);
+    chi2Result->Divide(chi2Result, mc, 1, 1, "B");
+    bayesResult->Divide(bayesResult, mc, 1, 1, "B");
+
+    // skip errors for now
+    for (Int_t j=0; j<=chi2Result->GetNbinsX(); ++j)
+    {
+      chi2Result->SetBinError(j, 0);
+      bayesResult->SetBinError(j, 0);
+    }
+
+    chi2Result->SetTitle("Ratios;Npart;unfolded measured/MC");
+    chi2Result->GetYaxis()->SetRangeUser(0.5, 1.5);
+
+    chi2Result->DrawCopy("");
+    bayesResult->DrawCopy("SAME");
+  }
+
+  canvas->SaveAs(Form("%s.gif", canvas->GetName()));
+  canvas->SaveAs(Form("%s.C", canvas->GetName()));
+
+  TCanvas* canvas2 = new TCanvas("StatisticsPlot2", "StatisticsPlot2", 800, 400);
+  canvas2->Divide(2, 1);
+
+  canvas2->cd(1);
+  fitResultsChi2->SetMarkerStyle(20);
+  fitResultsChi2->GetYaxis()->SetRangeUser(0.5 * min, 1.5 * max);
+  fitResultsChi2->Draw("AP");
+
+  fitResultsBayes->SetMarkerStyle(3);
+  fitResultsBayes->SetMarkerColor(2);
+  fitResultsBayes->Draw("P SAME");
+
+  gPad->SetLogy();
+
+  canvas2->cd(2);
+  fitResultsChi2Limit->SetMarkerStyle(20);
+  fitResultsChi2Limit->GetYaxis()->SetRangeUser(0.9 * TMath::Min(fitResultsChi2Limit->GetYaxis()->GetXmin(), fitResultsBayesLimit->GetYaxis()->GetXmin()), 1.1 * TMath::Max(fitResultsChi2Limit->GetYaxis()->GetXmax(), fitResultsBayesLimit->GetYaxis()->GetXmax()));
+  fitResultsChi2Limit->Draw("AP");
+
+  fitResultsBayesLimit->SetMarkerStyle(3);
+  fitResultsBayesLimit->SetMarkerColor(2);
+  fitResultsBayesLimit->Draw("P SAME");
+
+  canvas2->SaveAs(Form("%s.gif", canvas2->GetName()));
+  canvas2->SaveAs(Form("%s.C", canvas2->GetName()));
+}
+
+void StartingConditions(const char* fileNameMC = "multiplicityMC.root", Int_t histID = 2)
+{
+  gSystem->Load("libPWG0base");
+
+  AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
+
+  TFile::Open(fileNameMC);
+  mult->LoadHistograms("Multiplicity");
+
+  const char* files[] = { "multiplicityMC_0.5M.root", "multiplicityMC_0.75M.root", "multiplicityMC_1M_3.root", "multiplicityMC_1.25M.root", "multiplicityMC_1.5M.root" };
+
+  // this one we try to unfold
+  TFile::Open(files[0]);
+  AliMultiplicityCorrection* multESD = new AliMultiplicityCorrection("MultiplicityESD", "MultiplicityESD");
+  multESD->LoadHistograms("Multiplicity");
+  mult->SetMultiplicityESD(histID, multESD->GetMultiplicityESD(histID));
+  TH1* mc = multESD->GetMultiplicityMC(histID, AliMultiplicityCorrection::kTrVtx)->ProjectionY();
+
+  TGraph* fitResultsChi2 = new TGraph;
+  fitResultsChi2->SetTitle(";Input Dist ID;Chi2");
+  TGraph* fitResultsBayes = new TGraph;
+  fitResultsBayes->SetTitle(";Input Dist ID;Chi2");
+  TGraph* fitResultsChi2Limit = new TGraph;
+  fitResultsChi2Limit->SetTitle(";Input Dist ID;Multiplicity reach");
+  TGraph* fitResultsBayesLimit = new TGraph;
+  fitResultsBayesLimit->SetTitle(";Input Dist ID;Multiplicity reach");
+
+  TCanvas* canvas = new TCanvas("StartingConditions", "StartingConditions", 1200, 600);
+  canvas->Divide(8, 2);
+
+  TCanvas* canvas3 = new TCanvas("StartingConditions3", "StartingConditions3", 1000, 400);
+  canvas3->Divide(2, 1);
+
+  Float_t min = 1e10;
+  Float_t max = 0;
+
+  TH1* firstChi = 0;
+  TH1* firstBayesian = 0;
+  TH1* startCond = multESD->GetMultiplicityESD(histID)->ProjectionY("startCond");
+
+  TLegend* legend = new TLegend(0.7, 0.7, 1, 1);
+
+  for (Int_t i=0; i<8; ++i)
+  {
+    if (i == 0)
+    {
+      startCond = (TH1*) mc->Clone("startCond2");
+    }
+    else if (i < 6)
+    {
+      TFile::Open(files[i-1]);
+      AliMultiplicityCorrection* multESD2 = new AliMultiplicityCorrection("MultiplicityESD2", "MultiplicityESD2");
+      multESD2->LoadHistograms("Multiplicity");
+      startCond = multESD2->GetMultiplicityESD(histID)->ProjectionY("startCond");
+    }
+    else if (i == 6)
+    {
+      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, 1e-3, 1e10);
+      func->SetParLimits(1, 0.001, 1000);
+      func->SetParLimits(2, 0.001, 1000);
+
+      func->SetParameters(1, 10, 2);
+      for (Int_t j=2; j<=startCond->GetNbinsX(); j++)
+        startCond->SetBinContent(j, func->Eval(j-1));
+    }
+    else if (i == 7)
+    {
+      for (Int_t j=1; j<=startCond->GetNbinsX(); j++)
+        startCond->SetBinContent(j, 1);
+    }
+
+    mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 0.3);
+    mult->ApplyMinuitFit(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE, startCond);
+    mult->DrawComparison(Form("StartingConditions_%d_MinuitChi2", i), histID, kFALSE, kTRUE, mc);
+
+    Float_t chi2MC = 0;
+    Int_t chi2MCLimit = 0;
+    mult->GetComparisonResults(&chi2MC, &chi2MCLimit, 0);
+    fitResultsChi2->SetPoint(fitResultsChi2->GetN(), i, chi2MC);
+    fitResultsChi2Limit->SetPoint(fitResultsChi2Limit->GetN(), i, chi2MCLimit);
+    min = TMath::Min(min, chi2MC);
+    max = TMath::Max(max, chi2MC);
+
+    TH1* chi2Result = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone("chi2Result");
+    if (!firstChi)
+      firstChi = (TH1*) chi2Result->Clone("firstChi");
+
+    mult->ApplyBayesianMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, 0.1, 30, startCond);
+    mult->DrawComparison(Form("StartingConditions_%d_Bayesian", i), histID, kFALSE, kTRUE, mc);
+    TH1* bayesResult = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone("bayesResult");
+    if (!firstBayesian)
+      firstBayesian = (TH1*) bayesResult->Clone("firstBayesian");
+
+    mult->GetComparisonResults(&chi2MC, &chi2MCLimit, 0);
+    fitResultsBayes->SetPoint(fitResultsBayes->GetN(), i, chi2MC);
+    fitResultsBayesLimit->SetPoint(fitResultsBayesLimit->GetN(), i, chi2MCLimit);
+
+    min = TMath::Min(min, chi2MC);
+    max = TMath::Max(max, chi2MC);
+    mc->GetXaxis()->SetRangeUser(0, 150);
+    chi2Result->GetXaxis()->SetRangeUser(0, 150);
+
+    // skip errors for now
+    for (Int_t j=0; j<=chi2Result->GetNbinsX(); ++j)
+    {
+      chi2Result->SetBinError(j, 0);
+      bayesResult->SetBinError(j, 0);
+    }
+
+    canvas3->cd(1);
+    TH1* tmp = (TH1*) chi2Result->Clone("tmp");
+    tmp->SetTitle("Difference to best initial conditions;Npart;Ratio");
+    tmp->Divide(firstChi);
+    tmp->GetYaxis()->SetRangeUser(0.5, 1.5);
+    tmp->GetXaxis()->SetRangeUser(0, 200);
+    tmp->SetLineColor(i+1);
+    legend->AddEntry(tmp, Form("%d", i));
+    tmp->DrawCopy((i > 0) ? "SAME HIST" : "HIST");
+
+    canvas3->cd(2);
+    tmp = (TH1*) bayesResult->Clone("tmp");
+    tmp->SetTitle("Difference to best initial conditions;Npart;Ratio");
+    tmp->Divide(firstBayesian);
+    tmp->SetLineColor(i+1);
+    tmp->GetYaxis()->SetRangeUser(0.5, 1.5);
+    tmp->GetXaxis()->SetRangeUser(0, 200);
+    tmp->DrawCopy((i > 0) ? "SAME HIST" : "HIST");
+
+    canvas->cd(i+1);
+    mc->SetFillColor(kYellow);
+    mc->DrawCopy();
+    chi2Result->SetLineColor(kRed);
+    chi2Result->DrawCopy("SAME");
+    bayesResult->SetLineColor(kBlue);
+    bayesResult->DrawCopy("SAME");
+    gPad->SetLogy();
+
+    canvas->cd(i+9);
+    chi2Result->Divide(chi2Result, mc, 1, 1, "B");
+    bayesResult->Divide(bayesResult, mc, 1, 1, "B");
+
+    // skip errors for now
+    for (Int_t j=0; j<=chi2Result->GetNbinsX(); ++j)
+    {
+      chi2Result->SetBinError(j, 0);
+      bayesResult->SetBinError(j, 0);
+    }
+
+    chi2Result->SetTitle("Ratios;Npart;unfolded measured/MC");
+    chi2Result->GetYaxis()->SetRangeUser(0.5, 1.5);
+
+    chi2Result->DrawCopy("");
+    bayesResult->DrawCopy("SAME");
+  }
+
+  canvas3->cd(1);
+  legend->Draw();
+
   canvas->SaveAs(Form("%s.gif", canvas->GetName()));
+
+  TCanvas* canvas2 = new TCanvas("StartingConditions2", "StartingConditions2", 800, 400);
+  canvas2->Divide(2, 1);
+
+  canvas2->cd(1);
+  fitResultsChi2->SetMarkerStyle(20);
+  fitResultsChi2->GetYaxis()->SetRangeUser(0.5 * min, 1.5 * max);
+  fitResultsChi2->Draw("AP");
+
+  fitResultsBayes->SetMarkerStyle(3);
+  fitResultsBayes->SetMarkerColor(2);
+  fitResultsBayes->Draw("P SAME");
+
+  gPad->SetLogy();
+
+  canvas2->cd(2);
+  fitResultsChi2Limit->SetMarkerStyle(20);
+  fitResultsChi2Limit->GetYaxis()->SetRangeUser(0.9 * TMath::Min(fitResultsChi2Limit->GetYaxis()->GetXmin(), fitResultsBayesLimit->GetYaxis()->GetXmin()), 1.1 * TMath::Max(fitResultsChi2Limit->GetYaxis()->GetXmax(), fitResultsBayesLimit->GetYaxis()->GetXmax()));
+  fitResultsChi2Limit->Draw("AP");
+
+  fitResultsBayesLimit->SetMarkerStyle(3);
+  fitResultsBayesLimit->SetMarkerColor(2);
+  fitResultsBayesLimit->Draw("P SAME");
+
+  canvas2->SaveAs(Form("%s.gif", canvas2->GetName()));
+  canvas3->SaveAs(Form("%s.gif", canvas3->GetName()));
 }
 
 void Merge(Int_t n, const char** files, const char* output)
@@ -264,17 +889,251 @@ void testMethod(Int_t caseNo, const char* fileName = "multiplicityMC.root")
     case 1: func = new TF1("flat", "501-x"); break;
     case 2: func = new TF1("flat", "1000 * 1/(x+1)"); break;
     case 3: func = new TF1("flat", "1000 * TMath::Landau(x, 10, 5)"); break;
-    case 4: func->SetParameters(1000, 10, 2); break;
-    case 5: func->SetParameters(1000, 20, 3); break;
-    case 6: func->SetParameters(1000, 30, 4); break;
-    case 7: func->SetParameters(1000, 40, 5); break;
+    case 4: func->SetParameters(1e7, 10, 2); break;
+    case 5: func->SetParameters(1e7, 20, 3); break;
+    case 6: func->SetParameters(1e7, 30, 4); break;
+    case 7: func->SetParameters(1e7, 70, 2); break;
+    case 8: func = new TF1("testlaszlo", "10*1000*x*exp(-0.1*x)"); break;
 
     default: return;
   }
 
   mult->SetGenMeasFromFunc(func, 2);
 
-  mult->ApplyBayesianMethod(2, kFALSE);
-  mult->ApplyMinuitFit(2, kFALSE);
+  TFile::Open("out.root", "RECREATE");
+  mult->SaveHistograms();
+
+  //mult->ApplyBayesianMethod(2, kFALSE);
+  //mult->ApplyMinuitFit(2, kFALSE);
   //mult->ApplyGaussianMethod(2, kFALSE);
+  mult->ApplyLaszloMethod(2, kFALSE, AliMultiplicityCorrection::kTrVtx);
+}
+
+void smoothCorrelationMap(const char* fileName = "multiplicityMC.root", Int_t corrMatrix = 2)
+{
+  gSystem->Load("libPWG0base");
+
+  AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
+
+  TFile::Open(fileName);
+  mult->LoadHistograms("Multiplicity");
+
+  // empty under/overflow bins in x, otherwise Project3D takes them into account
+  TH3* corr = mult->GetCorrelation(corrMatrix);
+  for (Int_t j=1; j<=corr->GetYaxis()->GetNbins(); ++j)
+  {
+    for (Int_t k=1; k<=corr->GetZaxis()->GetNbins(); ++k)
+    {
+      corr->SetBinContent(0, j, k, 0);
+      corr->SetBinContent(corr->GetXaxis()->GetNbins()+1, j, k, 0);
+    }
+  }
+
+  TH2* proj = (TH2*) corr->Project3D("zy");
+
+  // normalize correction for given nPart
+  for (Int_t i=1; i<=proj->GetNbinsX(); ++i)
+  {
+    Double_t sum = proj->Integral(i, i, 1, proj->GetNbinsY());
+    if (sum <= 0)
+      continue;
+
+    for (Int_t j=1; j<=proj->GetNbinsY(); ++j)
+    {
+      // npart sum to 1
+      proj->SetBinContent(i, j, proj->GetBinContent(i, j) / sum);
+      proj->SetBinError(i, j, proj->GetBinError(i, j) / sum);
+    }
+  }
+
+  new TCanvas;
+  proj->Draw("COLZ");
+
+  TH1* scaling = proj->ProjectionY("scaling", 1, 1);
+  scaling->Reset();
+  scaling->SetMarkerStyle(3);
+  //scaling->GetXaxis()->SetRangeUser(0, 50);
+  TH1* mean = (TH1F*) scaling->Clone("mean");
+  TH1* width = (TH1F*) scaling->Clone("width");
+
+  TF1* lognormal = new TF1("lognormal", "[0]*exp(-(log(x)-[1])^2/(2*[2]^2))/(x*[2]*TMath::Sqrt(2*TMath::Pi()))", 0.01, 500);
+  lognormal->SetParNames("scaling", "mean", "sigma");
+  lognormal->SetParameters(1, 1, 1);
+  lognormal->SetParLimits(0, 1, 1);
+  lognormal->SetParLimits(1, 0, 100);
+  lognormal->SetParLimits(2, 1e-3, 1);
+
+  TF1* nbd = 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);
+  nbd->SetParNames("scaling", "averagen", "k");
+  nbd->SetParameters(1, 13, 5);
+  nbd->SetParLimits(0, 1, 1);
+  nbd->SetParLimits(1, 1, 100);
+  nbd->SetParLimits(2, 1, 1e8);
+
+  TF1* poisson = new TF1("poisson", "[0] * exp(-(x+[2])) * (x+[2])**[1] / TMath::Factorial([1])", 0.01, 50);
+  poisson->SetParNames("scaling", "k", "deltax");
+  poisson->SetParameters(1, 1, 0);
+  poisson->SetParLimits(0, 0, 10);
+  poisson->SetParLimits(1, 0.01, 100);
+  poisson->SetParLimits(2, 0, 10);
+
+  TF1* mygaus = new TF1("mygaus", "[0] * exp(-(x-[1])**2 / 2 / [2] - [3] * log(x + [4]) / [5])", 0.01, 50);
+  mygaus->SetParNames("scaling", "mean", "width", "scale2log", "logmean", "logwidth");
+  mygaus->SetParameters(1, 0, 1, 1, 0, 1);
+  mygaus->SetParLimits(2, 1e-5, 10);
+  mygaus->SetParLimits(4, 1, 1);
+  mygaus->SetParLimits(5, 1e-5, 10);
+
+  //TF1* sqrt = new TF1("sqrt", "[0] + [1] * sqrt((x + [3]) * [2])", 0, 50);
+  TF1* sqrt = new TF1("sqrt", "[0] + (x + [1])**[2]", 0, 50);
+  sqrt->SetParNames("ydelta", "exp", "xdelta");
+  sqrt->SetParameters(0, 0, 1);
+  sqrt->SetParLimits(1, 0, 10);
+
+  const char* fitWith = "gaus";
+
+  for (Int_t i=1; i<=150; ++i)
+  {
+    printf("Fitting %d...\n", i);
+
+    TH1* hist = proj->ProjectionY(Form("proj%d", i), i, i, "e");
+    //hist->GetXaxis()->SetRangeUser(0, 50);
+    //lognormal->SetParameter(0, hist->GetMaximum());
+    hist->Fit(fitWith, "0 M", "");
+
+    TF1* func = hist->GetListOfFunctions()->FindObject(fitWith);
+
+    if (((i-1) % 15 == 0) || ((i % 5 == 0) && i < 30))
+    {
+      new TCanvas;
+      hist->Draw();
+      func->Clone()->Draw("SAME");
+      gPad->SetLogy();
+    }
+
+    scaling->Fill(i, func->GetParameter(0));
+    mean->Fill(i, func->GetParameter(1));
+    width->Fill(i, func->GetParameter(2));
+  }
+
+  TF1* log = new TF1("log", "[0] + [1] * log([2] * x)", 0.01, 500);
+  log->SetParameters(0, 1, 1);
+  log->SetParLimits(1, 0, 100);
+  log->SetParLimits(2, 1e-3, 10);
+
+  TF1* over = new TF1("over", "[0] + [1] / (x+[2])", 0.01, 500);
+  over->SetParameters(0, 1, 0);
+  //over->SetParLimits(0, 0, 100);
+  over->SetParLimits(1, 1e-3, 10);
+  over->SetParLimits(2, 0, 100);
+
+  c1 = new TCanvas("fitparams", "fitparams", 1200, 400);
+  c1->Divide(3, 1);
+
+  c1->cd(1);
+  scaling->Draw("P");
+
+  //TF1* scalingFit = new TF1("mypol0", "[0]");
+  TF1* scalingFit = over;
+  scaling->Fit(scalingFit, "", "", 3, 100);
+
+  c1->cd(2);
+  mean->Draw("P");
+
+  //TF1* meanFit = log;
+  TF1* meanFit = new TF1("mypol1", "[0]+[1]*x");
+  mean->Fit(meanFit, "", "", 3, 100);
+
+  c1->cd(3);
+  width->Draw("P");
+
+  //TF1* widthFit = over;
+  TF1* widthFit = new TF1("mypol2", "[0]+[1]*x+[2]*x*x");
+  width->Fit(widthFit, "", "", 5, 100);
+
+  // build new correction matrix
+  TH2* new = (TH2*) proj->Clone("new");
+  new->Reset();
+  Float_t x, y;
+  for (Int_t i=1; i<=new->GetXaxis()->GetNbins(); i+=1)
+  {
+    TF1* func = (TF1*) gROOT->FindObject(fitWith);
+    x = new->GetXaxis()->GetBinCenter(i);
+    //if (i == 1)
+    //  x = 0.1;
+    x++;
+    func->SetParameters(scalingFit->Eval(x), meanFit->Eval(x), widthFit->Eval(x));
+    printf("%f %f %f %f\n", x, scalingFit->Eval(x), meanFit->Eval(x), widthFit->Eval(x));
+
+    for (Int_t j=1; j<=new->GetYaxis()->GetNbins(); j+=1)
+    {
+      if (i < 21)
+      {
+        // leave bins 1..20 untouched
+        new->SetBinContent(i, j, corr->Integral(1, corr->GetNbinsX(), i, i, j, j));
+      }
+      else
+      {
+        y = new->GetYaxis()->GetBinCenter(j);
+        if (j == 1)
+          y = 0.1;
+        if (func->Eval(y) > 1e-4)
+          new->SetBinContent(i, j, func->Eval(y));
+      }
+    }
+  }
+
+  // fill 0 multiplicity bins, this cannot be done with the function because it does not accept 0
+  // we take the values from the old response matrix
+  //for (Int_t i=1; i<=new->GetXaxis()->GetNbins(); i+=1)
+  //  new->SetBinContent(i, 1, proj->GetBinContent(i, 1));
+
+  //for (Int_t j=1; j<=new->GetYaxis()->GetNbins(); j+=1)
+  //  new->SetBinContent(1, j, proj->GetBinContent(1, j));
+
+  // normalize correction for given nPart
+  for (Int_t i=1; i<=new->GetNbinsX(); ++i)
+  {
+    Double_t sum = new->Integral(i, i, 1, proj->GetNbinsY());
+    if (sum <= 0)
+      continue;
+
+    for (Int_t j=1; j<=new->GetNbinsY(); ++j)
+    {
+      // npart sum to 1
+      new->SetBinContent(i, j, new->GetBinContent(i, j) / sum);
+      new->SetBinError(i, j, new->GetBinError(i, j) / sum);
+    }
+  }
+
+  new TCanvas;
+  new->Draw("COLZ");
+
+  TH2* diff = (TH2*) new->Clone("diff");
+  diff->Add(proj, -1);
+
+  new TCanvas;
+  diff->Draw("COLZ");
+  diff->SetMinimum(-0.05);
+  diff->SetMaximum(0.05);
+
+  corr->Reset();
+
+  for (Int_t i=1; i<=new->GetNbinsX(); ++i)
+    for (Int_t j=1; j<=new->GetNbinsY(); ++j)
+      corr->SetBinContent(corr->GetXaxis()->GetNbins() / 2, i, j, new->GetBinContent(i, j));
+
+  new TCanvas;
+  corr->Project3D("zy")->Draw("COLZ");
+
+  TFile::Open("out.root", "RECREATE");
+  mult->SaveHistograms();
+
+  TH1* proj1 = proj->ProjectionY("proj1", 36, 36);
+  TH1* proj2 = new->ProjectionY("proj2", 36, 36);
+  proj2->SetLineColor(2);
+
+  new TCanvas;
+  proj1->Draw();
+  proj2->Draw("SAME");
 }