extra bit for TPC and Global constrained flagging
[u/mrichter/AliRoot.git] / PWG0 / AliUnfolding.cxx
index f716551..945ae06 100644 (file)
 #include <TMath.h>
 #include <TCanvas.h>
 #include <TF1.h>
+#include <TExec.h>
+#include "Riostream.h"
+#include "TROOT.h"
+
+using namespace std; //required for resolving the 'cout' symbol
 
 TMatrixD* AliUnfolding::fgCorrelationMatrix = 0;
+TMatrixD* AliUnfolding::fgCorrelationMatrixSquared = 0;
 TMatrixD* AliUnfolding::fgCorrelationCovarianceMatrix = 0;
 TVectorD* AliUnfolding::fgCurrentESDVector = 0;
 TVectorD* AliUnfolding::fgEntropyAPriori = 0;
+TVectorD* AliUnfolding::fgEfficiency = 0;
+
+TAxis* AliUnfolding::fgUnfoldedAxis = 0;
+TAxis* AliUnfolding::fgMeasuredAxis = 0;
 
 TF1* AliUnfolding::fgFitFunction = 0;
 
@@ -44,13 +54,34 @@ AliUnfolding::RegularizationType AliUnfolding::fgRegularizationType = AliUnfoldi
 Float_t AliUnfolding::fgRegularizationWeight = 10000;
 Int_t AliUnfolding::fgSkipBinsBegin = 0;
 Float_t AliUnfolding::fgMinuitStepSize = 0.1;                 // (usually not needed to be changed) step size in minimization
-Bool_t AliUnfolding::fgNormalizeInput = kFALSE;                  // normalize input spectrum
+Float_t AliUnfolding::fgMinuitPrecision = 1e-6;               // minuit precision
+Bool_t AliUnfolding::fgMinimumInitialValue = kFALSE;          // set all initial values at least to the smallest value among the initial values
+Float_t AliUnfolding::fgMinimumInitialValueFix = -1;
+Bool_t AliUnfolding::fgNormalizeInput = kFALSE;               // normalize input spectrum
+Float_t AliUnfolding::fgNotFoundEvents = 0;
+Bool_t AliUnfolding::fgSkipBin0InChi2 = kFALSE;
 
 Float_t AliUnfolding::fgBayesianSmoothing  = 1;           // smoothing parameter (0 = no smoothing)
 Int_t   AliUnfolding::fgBayesianIterations = 10;          // number of iterations in Bayesian method
 
 Bool_t AliUnfolding::fgDebug = kFALSE;
 
+Int_t AliUnfolding::fgCallCount = 0;
+
+Int_t AliUnfolding::fgPowern = 5;
+
+Double_t AliUnfolding::fChi2FromFit = 0.;
+Double_t AliUnfolding::fPenaltyVal  = 0.;
+Double_t AliUnfolding::fAvgResidual = 0.;
+
+Int_t AliUnfolding::fgPrintChi2Details = 0;
+
+TCanvas *AliUnfolding::fgCanvas = 0;
+TH1 *AliUnfolding::fghUnfolded = 0;     
+TH2 *AliUnfolding::fghCorrelation = 0;  
+TH1 *AliUnfolding::fghEfficiency = 0;   
+TH1 *AliUnfolding::fghMeasured = 0;     
+
 ClassImp(AliUnfolding)
 
 //____________________________________________________________________
@@ -97,6 +128,47 @@ void AliUnfolding::SetNbins(Int_t nMeasured, Int_t nUnfolded)
   fgMaxInput = nMeasured; 
   fgMaxParams = nUnfolded; 
   
+  if (fgCorrelationMatrix)
+  {
+    delete fgCorrelationMatrix;
+    fgCorrelationMatrix = 0;
+  }
+  if (fgCorrelationMatrixSquared)
+  {
+    fgCorrelationMatrixSquared = 0;
+    delete fgCorrelationMatrixSquared;
+  }
+  if (fgCorrelationCovarianceMatrix)
+  {
+    delete fgCorrelationCovarianceMatrix;
+    fgCorrelationCovarianceMatrix = 0;
+  }
+  if (fgCurrentESDVector)
+  {
+    delete fgCurrentESDVector;
+    fgCurrentESDVector = 0;
+  }
+  if (fgEntropyAPriori)
+  {
+    delete fgEntropyAPriori;
+    fgEntropyAPriori = 0;
+  }
+  if (fgEfficiency)
+  {
+    delete fgEfficiency;
+    fgEfficiency = 0;
+  }
+  if (fgUnfoldedAxis)
+  {
+    delete fgUnfoldedAxis;
+    fgUnfoldedAxis = 0;
+  }
+  if (fgMeasuredAxis)
+  {
+    delete fgMeasuredAxis;
+    fgMeasuredAxis = 0;
+  }
+  
   Printf("AliUnfolding::SetNbins: Set %d measured bins and %d unfolded bins", nMeasured, nUnfolded);
 }
 
@@ -148,6 +220,8 @@ Int_t AliUnfolding::Unfold(TH2* correlation, TH1* efficiency, TH1* measured, TH1
   //  initialConditions: (optional) initial conditions for the unfolding. if 0 the measured spectrum is used as initial conditions.
   //  result:      target for the unfolded result
   //  check:       depends on the unfolding method, see comments in specific functions
+  //
+  //  return code: see UnfoldWithMinuit/UnfoldWithBayesian/UnfoldWithFunction
 
   if (fgMaxInput == -1)
   {
@@ -177,22 +251,35 @@ Int_t AliUnfolding::Unfold(TH2* correlation, TH1* efficiency, TH1* measured, TH1
     case kFunction:
       return UnfoldWithFunction(correlation, efficiency, measured, initialConditions, result);
   }
+
+
+
   return -1;
 }
 
 //____________________________________________________________________
-void AliUnfolding::SetStaticVariables(TH2* correlation, TH1* measured)
+void AliUnfolding::SetStaticVariables(TH2* correlation, TH1* measured, TH1* efficiency)
 {
   // fill static variables needed for minuit fit
 
   if (!fgCorrelationMatrix)
     fgCorrelationMatrix = new TMatrixD(fgMaxInput, fgMaxParams);
+  if (!fgCorrelationMatrixSquared)
+    fgCorrelationMatrixSquared = new TMatrixD(fgMaxInput, fgMaxParams);
   if (!fgCorrelationCovarianceMatrix)
     fgCorrelationCovarianceMatrix = new TMatrixD(fgMaxInput, fgMaxInput);
   if (!fgCurrentESDVector)
     fgCurrentESDVector = new TVectorD(fgMaxInput);
   if (!fgEntropyAPriori)
     fgEntropyAPriori = new TVectorD(fgMaxParams);
+  if (!fgEfficiency)
+    fgEfficiency = new TVectorD(fgMaxParams);
+  if (!fgUnfoldedAxis)
+    delete fgUnfoldedAxis;
+  fgUnfoldedAxis = new TAxis(*(correlation->GetXaxis()));
+  if (!fgMeasuredAxis)
+    delete fgMeasuredAxis;
+  fgMeasuredAxis = new TAxis(*(correlation->GetYaxis()));    
 
   fgCorrelationMatrix->Zero();
   fgCorrelationCovarianceMatrix->Zero();
@@ -217,30 +304,52 @@ void AliUnfolding::SetStaticVariables(TH2* correlation, TH1* measured)
       }
 
       // npart sum to 1
-      correlation->SetBinContent(i, j, correlation->GetBinContent(i, j) / sum);
+      correlation->SetBinContent(i, j, correlation->GetBinContent(i, j) / sum);// * correlation->GetXaxis()->GetBinWidth(i));
       correlation->SetBinError(i, j, correlation->GetBinError(i, j) / sum);
 
       if (i <= fgMaxParams && j <= fgMaxInput)
+      {
         (*fgCorrelationMatrix)(j-1, i-1) = correlation->GetBinContent(i, j);
+        (*fgCorrelationMatrixSquared)(j-1, i-1) = correlation->GetBinContent(i, j) * correlation->GetBinContent(i, j);
+      }
     }
 
     //printf("MPV for Ntrue = %f is %f\n", fCurrentCorrelation->GetXaxis()->GetBinCenter(i), fCurrentCorrelation->GetYaxis()->GetBinCenter(maxBin));
   }
     
   //normalize measured
+  Float_t smallestError = 1;
   if (fgNormalizeInput)
-    measured->Scale(1.0 / measured->Integral());
+  {
+    Float_t sumMeasured = measured->Integral();
+    measured->Scale(1.0 / sumMeasured);
+    smallestError /= sumMeasured;
+  }
   
   for (Int_t i=0; i<fgMaxInput; ++i)
   {
     (*fgCurrentESDVector)[i] = measured->GetBinContent(i+1);
     if (measured->GetBinError(i+1) > 0)
+    {
       (*fgCorrelationCovarianceMatrix)(i, i) = (Double_t) 1e-6 / measured->GetBinError(i+1) / measured->GetBinError(i+1);
+    }
+    else // in this case put error of 1, otherwise 0 bins are not added to the chi2...
+      (*fgCorrelationCovarianceMatrix)(i, i) = (Double_t) 1e-6 / smallestError / smallestError;
 
     if ((*fgCorrelationCovarianceMatrix)(i, i) > 1e7)
       (*fgCorrelationCovarianceMatrix)(i, i) = 0;
     //Printf("%d, %e", i, (*fgCorrelationCovarianceMatrix)(i, i));
   }
+
+  // efficiency is expected to match bin width of result
+  for (Int_t i=0; i<fgMaxParams; ++i)
+  {
+    (*fgEfficiency)(i) = efficiency->GetBinContent(i+1);
+  }
+
+  if (correlation->GetNbinsX() != fgMaxParams || correlation->GetNbinsY() != fgMaxInput)
+    cout << "Response histo has incorrect dimensions; expect (" << fgMaxParams << ", " << fgMaxInput << "), got (" << correlation->GetNbinsX() << ", " << correlation->GetNbinsY() << ")" << endl;
+
 }
 
 //____________________________________________________________________
@@ -252,16 +361,22 @@ Int_t AliUnfolding::UnfoldWithMinuit(TH2* correlation, TH1* efficiency, TH1* mea
   // unfolds <measured> using response from <correlation> and effiency <efficiency>
   // output is in <result>
   // <initialConditions> set the initial values for the minimization, if 0 <measured> is used
+  //   negative values in initialConditions mean that the given parameter is fixed to the absolute of the value
   // if <check> is true no unfolding is made, instead only the chi2 without unfolding is printed
   //
   // returns minuit status (0 = success), (-1 when check was set)
   //
 
-  SetStaticVariables(correlation, measured);
+  SetStaticVariables(correlation, measured, efficiency);
   
   // Initialize TMinuit via generic fitter interface
-  TVirtualFitter *minuit = TVirtualFitter::Fitter(0, fgMaxParams);
+  Int_t params = fgMaxParams;
+  if (fgNotFoundEvents > 0)
+    params++;
+    
+  TVirtualFitter *minuit = TVirtualFitter::Fitter(0, params);
   Double_t arglist[100];
+  //  minuit->SetDefaultFitter("Minuit2");
 
   // disable any output (-1), unfortuantly we do not see warnings anymore then. Have to find another way...
   arglist[0] = 0;
@@ -273,9 +388,19 @@ Int_t AliUnfolding::UnfoldWithMinuit(TH2* correlation, TH1* efficiency, TH1* mea
   // set minimization function
   minuit->SetFCN(Chi2Function);
 
+  // set precision
+  minuit->SetPrecision(fgMinuitPrecision);
+
   for (Int_t i=0; i<fgMaxParams; i++)
     (*fgEntropyAPriori)[i] = 1;
 
+  // set initial conditions as a-priori distribution for MRX regularization
+  /*
+  for (Int_t i=0; i<fgMaxParams; i++)
+    if (initialConditions && initialConditions->GetBinContent(i+1) > 0)
+      (*fgEntropyAPriori)[i] = initialConditions->GetBinContent(i+1);
+  */
+
   if (!initialConditions) {
     initialConditions = measured;
   } else {
@@ -285,20 +410,45 @@ Int_t AliUnfolding::UnfoldWithMinuit(TH2* correlation, TH1* efficiency, TH1* mea
       initialConditions->Scale(1.0 / initialConditions->Integral());
   }
 
-  // set initial conditions as a-priori distribution for MRX regularization
-  for (Int_t i=0; i<fgMaxParams; i++)
-    if (initialConditions->GetBinContent(i+1) > 0)
-      (*fgEntropyAPriori)[i] = initialConditions->GetBinContent(i+1);
-
+  // extract minimum value from initial conditions (if we set a value to 0 it will stay 0)
+  Float_t minValue = 1e35;
+  if (fgMinimumInitialValueFix < 0)
+  {
+    for (Int_t i=0; i<fgMaxParams; ++i)
+    {
+      Int_t bin = initialConditions->GetXaxis()->FindBin(result->GetXaxis()->GetBinCenter(i+1));
+      if (initialConditions->GetBinContent(bin) > 0)
+       minValue = TMath::Min(minValue, (Float_t) initialConditions->GetBinContent(bin));
+    }
+  }
+  else
+    minValue = fgMinimumInitialValueFix;
+  
   Double_t* results = new Double_t[fgMaxParams+1];
   for (Int_t i=0; i<fgMaxParams; ++i)
   {
-    results[i] = initialConditions->GetBinContent(i+1);
+    Int_t bin = initialConditions->GetXaxis()->FindBin(result->GetXaxis()->GetBinCenter(i+1));
+    results[i] = initialConditions->GetBinContent(bin);
 
+    Bool_t fix = kFALSE;
+    if (results[i] < 0)
+    {
+      fix = kTRUE;
+      results[i] = -results[i];
+    }
+    if (!fix && fgMinimumInitialValue && results[i] < minValue)
+      results[i] = minValue;
+      
     // minuit sees squared values to prevent it from going negative...
     results[i] = TMath::Sqrt(results[i]);
 
-    minuit->SetParameter(i, Form("param%d", i), results[i], fgMinuitStepSize, 0, 0);
+    minuit->SetParameter(i, Form("param%d", i), results[i], (fix) ? 0 : fgMinuitStepSize, 0, 0);
+  }
+  if (fgNotFoundEvents > 0)
+  {
+    results[fgMaxParams] = efficiency->GetBinContent(1);
+    minuit->SetParameter(fgMaxParams, "vtx0", results[fgMaxParams], fgMinuitStepSize / 100, 0.01, 0.80);
   }
   
   Int_t dummy = 0;
@@ -316,40 +466,54 @@ Int_t AliUnfolding::UnfoldWithMinuit(TH2* correlation, TH1* efficiency, TH1* mea
   // first param is number of iterations, second is precision....
   arglist[0] = 1e6;
   //arglist[1] = 1e-5;
-  //minuit->ExecuteCommand("SCAN", arglist, 0);
+  //  minuit->ExecuteCommand("SET PRINT", arglist, 3);
+  //  minuit->ExecuteCommand("SCAN", arglist, 0);
   Int_t status = minuit->ExecuteCommand("MIGRAD", arglist, 1);
   Printf("AliUnfolding::UnfoldWithMinuit: MINUIT status is %d", status);
   //printf("!!!!!!!!!!!!!! MIGRAD finished: Starting MINOS !!!!!!!!!!!!!!");
   //minuit->ExecuteCommand("MINOS", arglist, 0);
 
+  if (fgNotFoundEvents > 0)
+  {
+    results[fgMaxParams] = minuit->GetParameter(fgMaxParams);
+    Printf("Efficiency for bin 0 changed from %f to %f", efficiency->GetBinContent(1), results[fgMaxParams]);
+    efficiency->SetBinContent(1, results[fgMaxParams]);
+  }
+  
   for (Int_t i=0; i<fgMaxParams; ++i)
   {
     results[i] = minuit->GetParameter(i);
     Double_t value = results[i] * results[i];
     // error is : (relError) * (value) = (minuit->GetParError(i) / minuit->GetParameter(i)) * (minuit->GetParameter(i) * minuit->GetParameter(i))
-    Double_t error = minuit->GetParError(i) * results[i];
+    Double_t error = 0;
+    if (TMath::IsNaN(minuit->GetParError(i)))
+      Printf("WARNING: Parameter %d error is nan", i);
+    else 
+      error = minuit->GetParError(i) * results[i];
     
     if (efficiency)
-    {
+    {  
+      //printf("value before efficiency correction: %f\n",value);
       if (efficiency->GetBinContent(i+1) > 0)
       {
-        value /= efficiency->GetBinContent(i+1);
-        error /= efficiency->GetBinContent(i+1);
+       value /= efficiency->GetBinContent(i+1);
+       error /= efficiency->GetBinContent(i+1);
       }
       else
       {
-        value = 0;
-        error = 0;
+       value = 0;
+       error = 0;
       }
     }
-    
+    //printf("value after efficiency correction: %f +/- %f\n",value,error);
     result->SetBinContent(i+1, value);
     result->SetBinError(i+1, error);
   }
   
-  if (fgDebug)
-    DrawGuess(results);
-
+  fgCallCount = 0;
+  Chi2Function(dummy, 0, chi2, results, 0);
+  Printf("AliUnfolding::UnfoldWithMinuit: Chi2 of final parameters is = %f", chi2);
+  
   delete[] results;
 
   return status;
@@ -382,6 +546,7 @@ Int_t AliUnfolding::UnfoldWithBayesian(TH2* correlation, TH1* aEfficiency, TH1*
   Double_t* prior = new Double_t[kMaxT];
   Double_t* result = new Double_t[kMaxT];
   Double_t* efficiency = new Double_t[kMaxT];
+  Double_t* binWidths = new Double_t[kMaxT];
 
   Double_t** response = new Double_t*[kMaxT];
   Double_t** inverseResponse = new Double_t*[kMaxT];
@@ -416,6 +581,7 @@ Int_t AliUnfolding::UnfoldWithBayesian(TH2* correlation, TH1* aEfficiency, TH1*
       
     prior[t] = measuredCopy[t];
     result[t] = 0;
+    binWidths[t] = aResult->GetXaxis()->GetBinWidth(t+1);
   }
 
   // pick prior distribution
@@ -430,13 +596,14 @@ Int_t AliUnfolding::UnfoldWithBayesian(TH2* correlation, TH1* aEfficiency, TH1*
 
   //TH1F* convergence = new TH1F("convergence", "convergence", 200, 0.5, 200.5);
   
+  //new TCanvas;
   // unfold...
-  for (Int_t i=0; i<fgBayesianIterations || fgBayesianIterations< 0; i++)
+  for (Int_t i=0; i<fgBayesianIterations || fgBayesianIterations < 0; i++)
   {
     if (fgDebug)
       Printf("AliUnfolding::UnfoldWithBayesian: iteration %i", i);
 
-    // calculate IR from Beyes theorem
+    // calculate IR from Bayes theorem
     // IR_ji = R_ij * prior_i / sum_k(R_kj * prior_k)
 
     Double_t chi2Measured = 0;
@@ -478,15 +645,17 @@ Int_t AliUnfolding::UnfoldWithBayesian(TH2* correlation, TH1* aEfficiency, TH1*
         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->SetMarkerStyle(24+i);
     aResult->SetMarkerColor(2);
-    aResult->DrawCopy("P SAME HIST");
+    aResult->DrawCopy((i == 0) ? "P" : "PSAME");
     */
-
     Double_t chi2LastIter = 0;
     // regularization (simple smoothing)
     for (Int_t t=kStartBin; t<kMaxT; t++)
@@ -496,7 +665,7 @@ Int_t AliUnfolding::UnfoldWithBayesian(TH2* correlation, TH1* aEfficiency, TH1*
       // 0 bin excluded from smoothing
       if (t > kStartBin+2 && t<kMaxT-1)
       {
-        Float_t average = (result[t-1] + result[t] + result[t+1]) / 3;
+        Float_t average = (result[t-1] / binWidths[t-1] + result[t] / binWidths[t] + result[t+1] / binWidths[t+1]) / 3 * binWidths[t];
 
         // weight the average with the regularization parameter
         newValue = (1 - fgBayesianSmoothing) * result[t] + fgBayesianSmoothing * average;
@@ -526,14 +695,18 @@ Int_t AliUnfolding::UnfoldWithBayesian(TH2* correlation, TH1* aEfficiency, TH1*
   //new TCanvas; convergence->DrawCopy(); gPad->SetLogy();
   //delete convergence;
 
+  Float_t factor = 1;
+  if (!fgNormalizeInput)
+    factor = measuredIntegral;
   for (Int_t t=0; t<kMaxT; t++)
-    aResult->SetBinContent(t+1, result[t]);
+    aResult->SetBinContent(t+1, result[t] * factor);
 
   delete[] measuredCopy;
   delete[] measuredError;
   delete[] prior;
   delete[] result;
   delete[] efficiency;
+  delete[] binWidths;
 
   for (Int_t i=0; i<kMaxT; i++)
   {
@@ -678,18 +851,19 @@ Double_t AliUnfolding::RegularizationPol0(TVectorD& params)
   // homogenity term for minuit fitting
   // pure function of the parameters
   // prefers constant function (pol0)
-
+  //
+  // Does not take into account efficiency
   Double_t chi2 = 0;
 
   for (Int_t i=1+fgSkipBinsBegin; i<fgMaxParams; ++i)
   {
-    Double_t right  = params[i];
-    Double_t left   = params[i-1];
+    Double_t right  = params[i] / fgUnfoldedAxis->GetBinWidth(i+1);
+    Double_t left   = params[i-1] / fgUnfoldedAxis->GetBinWidth(i);
 
-    if (right != 0)
+    if (left != 0)
     {
-      Double_t diff = 1 - left / right;
-      chi2 += diff * diff;
+      Double_t diff = (right - left);
+      chi2 += diff * diff / left / ((fgUnfoldedAxis->GetBinWidth(i+1) + fgUnfoldedAxis->GetBinWidth(i)) / 2);
     }
   }
 
@@ -702,7 +876,8 @@ Double_t AliUnfolding::RegularizationPol1(TVectorD& params)
   // homogenity term for minuit fitting
   // pure function of the parameters
   // prefers linear function (pol1)
-
+  //
+  // Does not take into account efficiency
   Double_t chi2 = 0;
 
   for (Int_t i=2+fgSkipBinsBegin; i<fgMaxParams; ++i)
@@ -710,16 +885,16 @@ Double_t AliUnfolding::RegularizationPol1(TVectorD& params)
     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 right  = params[i] / fgUnfoldedAxis->GetBinWidth(i+1);
+    Double_t middle = params[i-1] / fgUnfoldedAxis->GetBinWidth(i);
+    Double_t left   = params[i-2] / fgUnfoldedAxis->GetBinWidth(i-1);
 
-    Double_t diff = (der1 - der2) / middle;
+    Double_t der1 = (right - middle) / ((fgUnfoldedAxis->GetBinWidth(i+1) + fgUnfoldedAxis->GetBinWidth(i)) / 2);
+    Double_t der2 = (middle - left) / ((fgUnfoldedAxis->GetBinWidth(i) + fgUnfoldedAxis->GetBinWidth(i-1)) / 2);
 
-    chi2 += diff * diff;
+    //Double_t diff = (der1 - der2) / middle;
+    //chi2 += diff * diff;
+    chi2 += (der1 - der2) * (der1 - der2) / middle * fgUnfoldedAxis->GetBinWidth(i);
   }
 
   return chi2;
@@ -730,28 +905,32 @@ Double_t AliUnfolding::RegularizationLog(TVectorD& params)
 {
   // homogenity term for minuit fitting
   // pure function of the parameters
-  // prefers linear function (pol1)
+  // prefers logarithmic function (log)
+  //
+  // Does not take into account efficiency
 
   Double_t chi2 = 0;
 
   for (Int_t i=2+fgSkipBinsBegin; i<fgMaxParams; ++i)
   {
-    if (params[i-1] == 0)
-      continue;
-
-    Double_t right  = log(params[i]);
-    Double_t middle = log(params[i-1]);
-    Double_t left   = log(params[i-2]);
-
-    Double_t der1 = (right - middle);
-    Double_t der2 = (middle - left);
+    if (params[i-1] == 0 || params[i] == 0 || params[i-2] == 0)
+     continue;
 
-    Double_t diff = (der1 - der2) / middle;
+    Double_t right  = log(params[i] / fgUnfoldedAxis->GetBinWidth(i+1));
+    Double_t middle = log(params[i-1] / fgUnfoldedAxis->GetBinWidth(i));
+    Double_t left   = log(params[i-2] / fgUnfoldedAxis->GetBinWidth(i-1));
+    
+    Double_t der1 = (right - middle) / ((fgUnfoldedAxis->GetBinWidth(i+1) + fgUnfoldedAxis->GetBinWidth(i)) / 2);
+    Double_t der2 = (middle - left) / ((fgUnfoldedAxis->GetBinWidth(i) + fgUnfoldedAxis->GetBinWidth(i-1)) / 2);
+    
+    //Double_t error = 1. / params[i] + 4. / params[i-1] + 1. / params[i-2];
 
-    chi2 += diff * diff;
+    //if (fgCallCount == 0)
+    //  Printf("%d %f %f", i, (der1 - der2) * (der1 - der2), error);
+    chi2 += (der1 - der2) * (der1 - der2);// / error;
   }
 
-  return chi2 * 100;
+  return chi2;
 }
 
 //____________________________________________________________________
@@ -761,6 +940,8 @@ Double_t AliUnfolding::RegularizationTotalCurvature(TVectorD& params)
   // 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.
+  //
+  // Does not take into account efficiency
 
   Double_t chi2 = 0;
 
@@ -788,6 +969,8 @@ Double_t AliUnfolding::RegularizationEntropy(TVectorD& params)
   // pure function of the parameters
   // calculates entropy, from
   // The method of reduced cross-entropy (M. Schmelling 1993)
+  //
+  // Does not take into account efficiency
 
   Double_t paramSum = 0;
   
@@ -797,17 +980,132 @@ Double_t AliUnfolding::RegularizationEntropy(TVectorD& params)
   Double_t chi2 = 0;
   for (Int_t i=fgSkipBinsBegin; i<fgMaxParams; ++i)
   {
-    //Double_t tmp = params[i] / paramSum;
-    Double_t tmp = params[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]);
     }
+    else
+      chi2 += 100;
+  }
+
+  return -chi2;
+}
+
+//____________________________________________________________________
+Double_t AliUnfolding::RegularizationRatio(TVectorD& params)
+{
+  // homogenity term for minuit fitting
+  // pure function of the parameters
+  //
+  // Does not take into account efficiency
+
+  Double_t chi2 = 0;
+
+  for (Int_t i=5+fgSkipBinsBegin; i<fgMaxParams; ++i)
+  {
+    if (params[i-1] == 0 || params[i] == 0)
+      continue;
+
+    Double_t right  = params[i] / fgUnfoldedAxis->GetBinWidth(i+1);
+    Double_t middle = params[i-1] / fgUnfoldedAxis->GetBinWidth(i);
+    Double_t left   = params[i-2] / fgUnfoldedAxis->GetBinWidth(i-1);
+    Double_t left2   = params[i-3] / fgUnfoldedAxis->GetBinWidth(i-2);
+    Double_t left3   = params[i-4] / fgUnfoldedAxis->GetBinWidth(i-3);
+    Double_t left4   = params[i-5] / fgUnfoldedAxis->GetBinWidth(i-4);
+
+    //Double_t diff = left / middle - middle / right;
+    //Double_t diff = 2 * left / middle - middle / right - left2 / left;
+    Double_t diff = 4 * left2 / left - middle / right - left / middle - left3 / left2 - left4 / left3;
+    
+    chi2 += diff * diff;// / middle;
+  }
+
+  return chi2;
+}
+
+//____________________________________________________________________
+Double_t AliUnfolding::RegularizationPowerLaw(TVectorD& params)
+{
+  // homogenity term for minuit fitting
+  // pure function of the parameters
+  // prefers power law with n = -5
+  //
+  // Does not take into account efficiency
+
+  Double_t chi2 = 0;
+
+  Double_t right = 0.;
+  Double_t middle = 0.;
+  Double_t left = 0.;
+
+  for (Int_t i=2+fgSkipBinsBegin; i<fgMaxParams; ++i)
+  {
+    if (params[i] < 1e-8 || params[i-1] < 1e-8 || params[i-2] < 1e-8)
+      continue;
+
+    if (fgUnfoldedAxis->GetBinWidth(i+1) < 1e-8 || fgUnfoldedAxis->GetBinWidth(i) < 1e-8 || fgUnfoldedAxis->GetBinWidth(i-1) < 1e-8)
+      continue;
+    
+    middle = TMath::Power(params[i-1] / fgUnfoldedAxis->GetBinWidth(i),fgPowern);
+
+    if(middle>0) {
+      right  = TMath::Power(params[i] / fgUnfoldedAxis->GetBinWidth(i),fgPowern)/middle;
+
+      left   = TMath::Power(params[i-2] / fgUnfoldedAxis->GetBinWidth(i-1),fgPowern)/middle;
+      
+      middle = 1.;
+      
+      Double_t der1 = (right - middle) / ((fgUnfoldedAxis->GetBinWidth(i+1) + fgUnfoldedAxis->GetBinWidth(i)) / 2);
+      Double_t der2 = (middle - left) / ((fgUnfoldedAxis->GetBinWidth(i-1) + fgUnfoldedAxis->GetBinWidth(i-2)) / 2);
+    
+      chi2 += (der1 - der2) * (der1 - der2)/ ( fgUnfoldedAxis->GetBinWidth(i)/2. + fgUnfoldedAxis->GetBinWidth(i-1) + fgUnfoldedAxis->GetBinWidth(i-2)/2.)/( fgUnfoldedAxis->GetBinWidth(i)/2. + fgUnfoldedAxis->GetBinWidth(i-1) + fgUnfoldedAxis->GetBinWidth(i-2)/2.);// / error;
+      //   printf("i: %d chi2 = %f\n",i,chi2);
+    }
+
+  }
+
+  return chi2;
+}
+
+//____________________________________________________________________
+Double_t AliUnfolding::RegularizationLogLog(TVectorD& params)
+{
+  // homogenity term for minuit fitting
+  // pure function of the parameters
+  // prefers a powerlaw (linear on a log-log scale)
+  //
+  // The calculation takes into account the efficiencies
+
+  Double_t chi2 = 0;
+
+  for (Int_t i=2+fgSkipBinsBegin; i<fgMaxParams; ++i)
+  {
+    if (params[i-1] == 0 || params[i] == 0 || params[i-2] == 0)
+     continue;
+    if ((*fgEfficiency)(i-1) == 0 || (*fgEfficiency)(i) == 0 || (*fgEfficiency)(i-2) == 0)
+     continue;
+
+
+    Double_t right  = log(params[i] / (*fgEfficiency)(i) / fgUnfoldedAxis->GetBinWidth(i));
+    Double_t middle = log(params[i-1] / (*fgEfficiency)(i-1) / fgUnfoldedAxis->GetBinWidth(i-1));
+    Double_t left   = log(params[i-2] / (*fgEfficiency)(i-2) / fgUnfoldedAxis->GetBinWidth(i-2));
+    
+    Double_t der1 = (right - middle) / ( log(fgUnfoldedAxis->GetBinCenter(i+1)) - log(fgUnfoldedAxis->GetBinCenter(i)) );
+    Double_t der2 = (middle - left) /( log(fgUnfoldedAxis->GetBinCenter(i)) - log(fgUnfoldedAxis->GetBinCenter(i-1)) );
+
+    double tmp = (log(fgUnfoldedAxis->GetBinCenter(i+1)) - log(fgUnfoldedAxis->GetBinCenter(i-1)))/2.;
+    Double_t dder = (der1-der2) / tmp;
+
+    chi2 += dder * dder;
   }
 
-  return 100.0 + chi2;
+  return chi2;
 }
 
+
+
 //____________________________________________________________________
 void AliUnfolding::Chi2Function(Int_t&, Double_t*, Double_t& chi2, Double_t *params, Int_t)
 {
@@ -815,14 +1113,17 @@ void AliUnfolding::Chi2Function(Int_t&, Double_t*, Double_t& chi2, Double_t *par
   // fit function for minuit
   // does: (m - Ad)W(m - Ad) where m = measured, A correlation matrix, d = guess, W = covariance matrix
   //
+  
+  // TODO use static members for the variables here to speed up processing (no construction/deconstruction)
 
-  // d
-  static TVectorD paramsVector(fgMaxParams);
+  // d = guess
+  TVectorD paramsVector(fgMaxParams);
   for (Int_t i=0; i<fgMaxParams; ++i)
     paramsVector[i] = params[i] * params[i];
 
   // calculate penalty factor
   Double_t penaltyVal = 0;
+
   switch (fgRegularizationType)
   {
     case kNone:       break;
@@ -831,12 +1132,12 @@ void AliUnfolding::Chi2Function(Int_t&, Double_t*, Double_t& chi2, Double_t *par
     case kCurvature:  penaltyVal = RegularizationTotalCurvature(paramsVector); break;
     case kEntropy:    penaltyVal = RegularizationEntropy(paramsVector); break;
     case kLog:        penaltyVal = RegularizationLog(paramsVector); break;
+    case kRatio:      penaltyVal = RegularizationRatio(paramsVector); break;
+    case kPowerLaw:   penaltyVal = RegularizationPowerLaw(paramsVector); break;
+    case kLogLog:     penaltyVal = RegularizationLogLog(paramsVector); break;
   }
 
-  //if (penaltyVal > 1e10)
-  //  paramsVector2.Print();
-
-  penaltyVal *= fgRegularizationWeight;
+  penaltyVal *= fgRegularizationWeight; //beta*PU
 
   // Ad
   TVectorD measGuessVector(fgMaxInput);
@@ -844,9 +1145,26 @@ void AliUnfolding::Chi2Function(Int_t&, Double_t*, Double_t& chi2, Double_t *par
 
   // Ad - m
   measGuessVector -= (*fgCurrentESDVector);
+  
+#if 0
+  // new error calcuation using error on the guess
+  
+  // error from guess
+  TVectorD errorGuessVector(fgMaxInput);
+  //errorGuessVector = (*fgCorrelationMatrixSquared) * paramsVector;
+  errorGuessVector = (*fgCorrelationMatrix) * paramsVector;
 
-  //measGuessVector.Print();
+  Double_t chi2FromFit = 0;
+  for (Int_t i=0; i<fgMaxInput; ++i)
+  {
+    Float_t error = 1;
+    if (errorGuessVector(i) > 0)
+      error = errorGuessVector(i);
+    chi2FromFit += measGuessVector(i) * measGuessVector(i) / error;
+  }
 
+#else
+  // old error calcuation using the error on the measured
   TVectorD copy(measGuessVector);
 
   // (Ad - m) W
@@ -857,29 +1175,240 @@ void AliUnfolding::Chi2Function(Int_t&, Double_t*, Double_t& chi2, Double_t *par
   for (Int_t i=0; i<fgMaxInput; ++i)
     measGuessVector[i] *= (*fgCorrelationCovarianceMatrix)(i, i);
 
-  //measGuessVector.Print();
+
+  if (fgSkipBin0InChi2)
+    measGuessVector[0] = 0;
 
   // (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)
+  // big number ((*fgCorrelationCovarianceMatrix)(i, i)) (see UnfoldWithMinuit)
   Double_t chi2FromFit = measGuessVector * copy * 1e6;
+#endif
+
+  Double_t notFoundEventsConstraint = 0;
+  Double_t currentNotFoundEvents = 0;
+  Double_t errorNotFoundEvents = 0;
+  
+  if (fgNotFoundEvents > 0)
+  {
+    for (Int_t i=0; i<fgMaxParams; ++i)
+    {
+      Double_t eff = (1.0 / (*fgEfficiency)(i) - 1);
+      if (i == 0)
+       eff = (1.0 / params[fgMaxParams] - 1);
+      currentNotFoundEvents += eff * paramsVector(i);
+      errorNotFoundEvents += eff * eff * paramsVector(i); // error due to guess (paramsVector)
+      if (i <= 3)
+       errorNotFoundEvents += (eff * 0.03) * (eff * 0.03) * paramsVector(i) * paramsVector(i); // error on eff
+      //      if ((fgCallCount % 10000) == 0)
+       //Printf("%d %f %f %f", i, (*fgEfficiency)(i), paramsVector(i), currentNotFoundEvents);
+    }
+    errorNotFoundEvents += fgNotFoundEvents;
+    // TODO add error on background, background estimate
+    
+    notFoundEventsConstraint = (currentNotFoundEvents - fgNotFoundEvents) * (currentNotFoundEvents - fgNotFoundEvents) / errorNotFoundEvents;
+  }
+  
+  Float_t avg = 0;
+  Float_t sum = 0;
+  Float_t currentMult = 0;
+  // try to match dn/deta
+  for (Int_t i=0; i<fgMaxParams; ++i)
+  {
+    avg += paramsVector(i) * currentMult;
+    sum += paramsVector(i);
+    currentMult += fgUnfoldedAxis->GetBinWidth(i);
+  }
+  avg /= sum;
+  Float_t chi2Avg = 0; //(avg - 3.73) * (avg - 3.73) * 100;
+
+  chi2 = chi2FromFit + penaltyVal + notFoundEventsConstraint + chi2Avg;
+
+  if ((fgCallCount++ % 1000) == 0)
+  {
+
+    Printf("AliUnfolding::Chi2Function: Iteration %d (ev %d %d +- %f) (%f) (%f): %f %f %f %f --> %f", fgCallCount-1, (Int_t) fgNotFoundEvents, (Int_t) currentNotFoundEvents, TMath::Sqrt(errorNotFoundEvents), params[fgMaxParams-1], avg, chi2FromFit, penaltyVal, notFoundEventsConstraint, chi2Avg, chi2);
+
+    //for (Int_t i=0; i<fgMaxInput; ++i)
+    //  Printf("%d: %f", i, measGuessVector(i) * copy(i) * 1e6);
+  }
+
+  fChi2FromFit = chi2FromFit;
+  fPenaltyVal = penaltyVal;
+}
+
+//____________________________________________________________________
+void AliUnfolding::MakePads() {
+    TPad *presult = new TPad("presult","result",0,0.4,1,1);
+    presult->SetNumber(1);
+    presult->Draw();
+    presult->SetLogy();
+    TPad *pres = new TPad("pres","residuals",0,0.2,1,0.4);
+    pres->SetNumber(2);
+    pres->Draw();
+    TPad *ppen = new TPad("ppen","penalty",0,0.0,1,0.2);
+    ppen->SetNumber(3);
+    ppen->Draw();
+}
+//____________________________________________________________________
+void AliUnfolding::DrawResults(TH2* correlation, TH1* efficiency, TH1* measured, TH1* initialConditions, TCanvas *canv, Int_t reuseHists,TH1 *unfolded)
+{
+  // Draw histograms of
+  //   - Result folded with response
+  //   - Penalty factors
+  //   - Chisquare contributions
+  // (Useful for debugging/sanity checks and the interactive unfolder)
+  //
+  // If a canvas pointer is given (canv != 0), it will be used for all
+  // plots; 3 pads are made if needed.
+
+
+  Int_t blankCanvas = 0;
+  TVirtualPad *presult = 0;
+  TVirtualPad *pres = 0; 
+  TVirtualPad *ppen = 0;
+  
+  if (canv) {
+    canv->cd();
+    presult = canv->GetPad(1);
+    pres = canv->GetPad(2);
+    ppen = canv->GetPad(3); 
+    if (presult == 0 || pres == 0 || ppen == 0) {
+      canv->Clear();
+      MakePads();
+      blankCanvas = 1;
+      presult = canv->GetPad(1);
+      pres = canv->GetPad(2);
+      ppen = canv->GetPad(3); 
+    } 
+  }
+  else {
+    presult = new TCanvas;
+    pres = new TCanvas;
+    ppen = new TCanvas;
+  }
+
+
+  if (fgMaxInput == -1)
+  {
+    Printf("AliUnfolding::Unfold: WARNING. Number of measured bins not set with SetNbins. Using number of bins in measured distribution");
+    fgMaxInput = measured->GetNbinsX();
+  }
+  if (fgMaxParams == -1)
+  {
+    Printf("AliUnfolding::Unfold: WARNING. Number of unfolded bins not set with SetNbins. Using number of bins in measured distribution");
+    fgMaxParams = initialConditions->GetNbinsX();
+  }
+
+  if (fgOverflowBinLimit > 0)
+    CreateOverflowBin(correlation, measured);
+
+  // Uses Minuit methods
 
-  chi2 = chi2FromFit + penaltyVal;
+  SetStaticVariables(correlation, measured, efficiency);
 
-  static Int_t callCount = 0;
-  if ((callCount++ % 10000) == 0)
-    Printf("AliUnfolding::Chi2Function: Iteration %d: %f %f --> %f", callCount, chi2FromFit, penaltyVal, chi2);
+  Double_t* params = new Double_t[fgMaxParams+1];
+  for (Int_t i=0; i<fgMaxParams; ++i)
+  {
+    params[i] = initialConditions->GetBinContent(i+1) * efficiency->GetBinContent(i+1);
+
+    Bool_t fix = kFALSE;
+    if (params[i] < 0)
+    {
+      fix = kTRUE;
+      params[i] = -params[i];
+    }
+    params[i]=TMath::Sqrt(params[i]);
+
+    //cout << "params[" << i << "] " << params[i] << endl;
+
+  } 
+
+  Double_t chi2;
+  Int_t dummy;
+
+  //fgPrintChi2Details = kTRUE;
+  fgCallCount = 0; // To make sure that Chi2Function prints the components
+  Chi2Function(dummy, 0, chi2, params, 0);
+
+  presult->cd();
+  TH1 *meas2 = (TH1*)measured->Clone();
+  meas2->SetLineColor(1);
+  meas2->SetLineWidth(2);
+  meas2->SetMarkerColor(meas2->GetLineColor());
+  meas2->SetMarkerStyle(20);
+  Float_t scale = unfolded->GetXaxis()->GetBinWidth(1)/meas2->GetXaxis()->GetBinWidth(1);
+  meas2->Scale(scale);
+  if (blankCanvas)
+    meas2->DrawCopy();
+  else
+    meas2->DrawCopy("same");
+  //meas2->GetXaxis()->SetLimits(0,fgMaxInput);
+  meas2->SetBit(kCannotPick);
+  DrawGuess(params, presult, pres, ppen, reuseHists,unfolded);
 }
+//____________________________________________________________________
+void AliUnfolding::RedrawInteractive() {
+  // 
+  // Helper function for interactive unfolding
+  //
+  DrawResults(fghCorrelation,fghEfficiency,fghMeasured,fghUnfolded,fgCanvas,1,fghUnfolded);
+}
+//____________________________________________________________________
+void AliUnfolding::InteractiveUnfold(TH2* correlation, TH1* efficiency, TH1* measured, TH1* initialConditions) 
+{
+  //
+  // Function to do interactive unfolding
+  // A canvas is drawn with the unfolding result
+  // Change the histogram with your mouse and all histograms 
+  // will be updated automatically
+
+  fgCanvas = new TCanvas("UnfoldingCanvas","Interactive unfolding",500,800);
+  fgCanvas->cd();
+
+  MakePads();
 
+  if (fghUnfolded) {
+    delete fghUnfolded; fghUnfolded = 0;
+  }
+  
+  gROOT->SetEditHistograms(kTRUE);
+
+  fghUnfolded = new TH1F("AluUnfoldingfghUnfolded","Unfolded result (Interactive unfolder",efficiency->GetNbinsX(),efficiency->GetXaxis()->GetXmin(),efficiency->GetXaxis()->GetXmax());
+  fghUnfolded->SetLineColor(2);
+  fghUnfolded->SetMarkerColor(2);
+  fghUnfolded->SetLineWidth(2);
+
+  fghCorrelation = correlation;
+  fghEfficiency = efficiency;
+  fghMeasured = measured;
+
+  UnfoldWithMinuit(correlation,efficiency,measured,initialConditions,fghUnfolded,kFALSE);
+
+  fgCanvas->cd(1);
+  fghUnfolded->Draw();
+  DrawResults(correlation,efficiency,measured,fghUnfolded,fgCanvas,kFALSE,fghUnfolded);
+
+  TExec *execRedraw = new TExec("redraw","AliUnfolding::RedrawInteractive()");
+  fghUnfolded->GetListOfFunctions()->Add(execRedraw);
+}
 //____________________________________________________________________
-void AliUnfolding::DrawGuess(Double_t *params)
+void AliUnfolding::DrawGuess(Double_t *params, TVirtualPad *pfolded, TVirtualPad *pres, TVirtualPad *ppen, Int_t reuseHists,TH1* unfolded)
 {
   //
   // draws residuals of solution suggested by params and effect of regularization
   //
 
+  if (pfolded == 0)
+    pfolded = new TCanvas;
+  if (ppen == 0)
+    ppen = new TCanvas;
+  if (pres == 0)
+    pres = new TCanvas;
+  
   // d
-  static TVectorD paramsVector(fgMaxParams);
+  TVectorD paramsVector(fgMaxParams);
   for (Int_t i=0; i<fgMaxParams; ++i)
     paramsVector[i] = params[i] * params[i];
 
@@ -887,11 +1416,34 @@ void AliUnfolding::DrawGuess(Double_t *params)
   TVectorD measGuessVector(fgMaxInput);
   measGuessVector = (*fgCorrelationMatrix) * paramsVector;
 
+  TH1* folded = 0;
+  if (reuseHists) 
+    folded = dynamic_cast<TH1*>(gROOT->FindObject("__hfolded"));
+  if (!reuseHists || folded == 0) {
+    if (fgMeasuredAxis->GetXbins()->GetArray()) // variable bins
+      folded = new TH1F("__hfolded","Folded histo from AliUnfolding",fgMeasuredAxis->GetNbins(),fgMeasuredAxis->GetXbins()->GetArray());
+    else
+      folded = new TH1F("__hfolded","Folded histo from AliUnfolding",fgMeasuredAxis->GetNbins(),fgMeasuredAxis->GetXmin(),fgMeasuredAxis->GetXmax());
+  }
+
+  folded->SetBit(kCannotPick);
+  folded->SetLineColor(4);
+  folded->SetLineWidth(2);
+
+  for (Int_t ibin =0; ibin < fgMaxInput; ibin++) 
+    folded->SetBinContent(ibin+1, measGuessVector[ibin]);
+
+  Float_t scale = unfolded->GetXaxis()->GetBinWidth(1)/folded->GetXaxis()->GetBinWidth(1);
+  folded->Scale(scale);
+
+  pfolded->cd();
+
+  folded->Draw("same");
+
   // 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
@@ -907,38 +1459,246 @@ void AliUnfolding::DrawGuess(Double_t *params)
   //Double_t chi2FromFit = measGuessVector * copy * 1e6;
 
   // draw residuals
-  TH1* residuals = new TH1F("residuals", "residuals", fgMaxInput, -0.5, fgMaxInput - 0.5);
-  for (Int_t i=0; i<fgMaxInput; ++i)
-    residuals->SetBinContent(i+1, measGuessVector(i) * copy(i) * 1e6);
-  new TCanvas; residuals->DrawCopy(); gPad->SetLogy();
+  // Double_t pTarray[fgMaxParams+1];
+  // for(int i=0; i<fgMaxParams; i++) {
+  //   pTarray[i] = fgUnfoldedAxis->GetBinCenter(i)-0.5*fgUnfoldedAxis->GetBinWidth(i);
+  // }
+  //  pTarray[fgMaxParams] = fgUnfoldedAxis->GetBinCenter(fgMaxParams-1)+0.5*fgUnfoldedAxis->GetBinWidth(fgMaxParams-1);
+  //  TH1* residuals = new TH1F("residuals", "residuals", fgMaxParams,pTarray);
+  //  TH1* residuals = new TH1F("residuals", "residuals", fgMaxInput, -0.5, fgMaxInput - 0.5);
+  // for (Int_t i=0; i<fgMaxInput; ++i)
+  //   residuals->SetBinContent(i+1, measGuessVector(i) * copy(i) * 1e6);7
+  TH1* residuals = GetResidualsPlot(params);
+  residuals->GetXaxis()->SetTitleSize(0.06);
+  residuals->GetXaxis()->SetTitleOffset(0.7);
+  residuals->GetXaxis()->SetLabelSize(0.07);
+  residuals->GetYaxis()->SetTitleSize(0.08);
+  residuals->GetYaxis()->SetTitleOffset(0.5);
+  residuals->GetYaxis()->SetLabelSize(0.06);
+  pres->cd(); residuals->DrawCopy(); gPad->SetLogy();
+    
 
   // draw penalty
-  if (fgRegularizationType != kPol1) {
-    Printf("Drawing not supported");
-    return;
+  TH1* penalty = GetPenaltyPlot(params);
+  penalty->GetXaxis()->SetTitleSize(0.06);
+  penalty->GetXaxis()->SetTitleOffset(0.7);
+  penalty->GetXaxis()->SetLabelSize(0.07);
+  penalty->GetYaxis()->SetTitleSize(0.08);
+  penalty->GetYaxis()->SetTitleOffset(0.5);
+  penalty->GetYaxis()->SetLabelSize(0.06);
+  ppen->cd(); penalty->DrawCopy(); gPad->SetLogy();
+}
+
+//____________________________________________________________________
+TH1* AliUnfolding::GetResidualsPlot(TH1* corrected)
+{
+  //
+  // MvL: THIS MUST BE INCORRECT. 
+  // Need heff to calculate params from TH1 'corrected'
+  //
+
+  //
+  // fill residuals histogram of solution suggested by params and effect of regularization
+  //
+
+  Double_t* params = new Double_t[fgMaxParams];
+  for (Int_t i=0; i<TMath::Min(fgMaxParams, corrected->GetNbinsX()); i++)
+    params[i] = TMath::Sqrt(TMath::Abs(corrected->GetBinContent(i+1)*(*fgEfficiency)(i)));
+
+
+  return GetResidualsPlot(params);
+}
+
+//____________________________________________________________________
+TH1* AliUnfolding::GetResidualsPlot(Double_t* params)
+{
+  //
+  // fill residuals histogram of solution suggested by params and effect of regularization
+  //
+
+  // d
+  TVectorD paramsVector(fgMaxParams);
+  for (Int_t i=0; i<fgMaxParams; ++i)
+    paramsVector[i] = params[i] * params[i];
+
+  // Ad
+  TVectorD measGuessVector(fgMaxInput);
+  measGuessVector = (*fgCorrelationMatrix) * paramsVector;
+
+  // Ad - m
+  measGuessVector -= (*fgCurrentESDVector);
+
+  TVectorD copy(measGuessVector);
+
+  // (Ad - m) W
+  // this step can be optimized because currently only the diagonal elements of fgCorrelationCovarianceMatrix are used
+  // normal way is like this:
+  // measGuessVector *= (*fgCorrelationCovarianceMatrix);
+  // optimized way like this:
+  for (Int_t i=0; i<fgMaxInput; ++i)
+    measGuessVector[i] *= (*fgCorrelationCovarianceMatrix)(i, i);
+
+  // (Ad - m) W (Ad - m)
+  // the factor 1e6 prevents that a very small number (measGuessVector[i]) is multiplied with a very
+  // big number ((*fgCorrelationCovarianceMatrix)(i, i)) (see ApplyMinuitFit)
+  //Double_t chi2FromFit = measGuessVector * copy * 1e6;
+
+  // draw residuals
+  TH1* residuals = 0;
+  if (fgMeasuredAxis->GetXbins()->GetArray()) // variable bins
+    residuals = new TH1F("residuals", "residuals;unfolded pos;residual",fgMeasuredAxis->GetNbins(),fgMeasuredAxis->GetXbins()->GetArray());
+  else
+    residuals = new TH1F("residuals", "residuals;unfolded pos;residual",fgMeasuredAxis->GetNbins(),fgMeasuredAxis->GetXmin(), fgMeasuredAxis->GetXmax());
+  //  TH1* residuals = new TH1F("residuals", "residuals", fgMaxInput, -0.5, fgMaxInput - 0.5);
+
+  Double_t sumResiduals = 0.; 
+  for (Int_t i=0; i<fgMaxInput; ++i) {
+    residuals->SetBinContent(i+1, measGuessVector(i) * copy(i) * 1e6);
+    sumResiduals += measGuessVector(i) * copy(i) * 1e6;
   }
+  fAvgResidual = sumResiduals/(double)fgMaxInput;
+  return residuals;
+}
 
-  TH1* penalty = new TH1F("penalty", "penalty", fgMaxParams, -0.5, fgMaxParams - 0.5);
+//____________________________________________________________________
+TH1* AliUnfolding::GetPenaltyPlot(TH1* corrected)
+{
+  // draws the penalty factors as function of multiplicity of the current selected regularization
 
-  for (Int_t i=2+1; i<fgMaxParams; ++i)
+  Double_t* params = new Double_t[fgMaxParams];
+  for (Int_t i=0; i<TMath::Min(fgMaxParams, corrected->GetNbinsX()); i++)
+    params[i] = (*fgEfficiency)(i)*corrected->GetBinContent(i+1);
+  
+  TH1* penalty = GetPenaltyPlot(params);
+  
+  delete[] params;
+  
+  return penalty;
+}
+
+//____________________________________________________________________
+TH1* AliUnfolding::GetPenaltyPlot(Double_t* params)
+{
+  // draws the penalty factors as function of multiplicity of the current selected regularization
+  
+  //TH1* penalty = new TH1F("penalty", ";unfolded multiplicity;penalty factor", fgMaxParams, -0.5, fgMaxParams - 0.5);
+  //  TH1* penalty = new TH1F("penalty", ";unfolded pos;penalty factor", fgMaxParams, fgUnfoldedAxis->GetBinCenter(0)-0.5*fgUnfoldedAxis->GetBinWidth(0),fgUnfoldedAxis->GetBinCenter(fgMaxParams)+0.5*fgUnfoldedAxis->GetBinWidth(fgMaxParams) );
+
+  TH1* penalty = 0;
+  if (fgUnfoldedAxis->GetXbins()->GetArray())
+    penalty = new TH1F("penalty", ";unfolded pos;penalty factor", fgUnfoldedAxis->GetNbins(),fgUnfoldedAxis->GetXbins()->GetArray());
+  else
+    penalty = new TH1F("penalty", ";unfolded pos;penalty factor", fgUnfoldedAxis->GetNbins(),fgUnfoldedAxis->GetXmin(),fgUnfoldedAxis->GetXmax());
+
+  for (Int_t i=1+fgSkipBinsBegin; i<fgMaxParams; ++i)
   {
-    if (params[i-1] == 0)
-      continue;
+    Double_t diff = 0;
+    if (fgRegularizationType == kPol0)
+    {
+      Double_t right  = params[i] / fgUnfoldedAxis->GetBinWidth(i+1);
+      Double_t left   = params[i-1] / fgUnfoldedAxis->GetBinWidth(i);
 
-    Double_t right  = params[i];
-    Double_t middle = params[i-1];
-    Double_t left   = params[i-2];
+      if (left != 0)
+      {
+        Double_t diffTmp = (right - left);
+        diff = diffTmp * diffTmp / left / ((fgUnfoldedAxis->GetBinWidth(i+1) + fgUnfoldedAxis->GetBinWidth(i)) / 2) / 100;
+      }
+    } 
+    if (fgRegularizationType == kPol1 && i > 1+fgSkipBinsBegin) 
+    {
+      if (params[i-1] == 0)
+        continue;
 
-    Double_t der1 = (right - middle);
-    Double_t der2 = (middle - left);
+      Double_t right  = params[i] / fgUnfoldedAxis->GetBinWidth(i+1);
+      Double_t middle = params[i-1] / fgUnfoldedAxis->GetBinWidth(i);
+      Double_t left   = params[i-2] / fgUnfoldedAxis->GetBinWidth(i-1);
+
+      Double_t der1 = (right - middle) / ((fgUnfoldedAxis->GetBinWidth(i+1) + fgUnfoldedAxis->GetBinWidth(i)) / 2);
+      Double_t der2 = (middle - left) / ((fgUnfoldedAxis->GetBinWidth(i) + fgUnfoldedAxis->GetBinWidth(i-1)) / 2);
+
+      diff = (der1 - der2) * (der1 - der2) / middle;
+    }
+
+    if (fgRegularizationType == kLog && i > 1+fgSkipBinsBegin) 
+    {
+      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 error = 1. / params[i] + 4. / params[i-1] + 1. / params[i-2];
+      //Printf("%d %f %f", i, (der1 - der2) * (der1 - der2), error);
+       
+      diff = (der1 - der2) * (der1 - der2);// / error;
+    }
+    if (fgRegularizationType == kCurvature && i > 1+fgSkipBinsBegin)
+    {
+      Double_t right  = params[i];    // params are sqrt
+      Double_t middle = params[i-1];
+      Double_t left   = params[i-2];
+      
+      Double_t der1 = (right - middle)/0.5/(fgUnfoldedAxis->GetBinWidth(i-1) + fgUnfoldedAxis->GetBinWidth(i));
+      Double_t der2 = (middle - left)/0.5/(fgUnfoldedAxis->GetBinWidth(i) + fgUnfoldedAxis->GetBinWidth(i+1));
+      
+      diff = (der1 - der2)/(fgUnfoldedAxis->GetBinWidth(i-1) + fgUnfoldedAxis->GetBinWidth(i) + fgUnfoldedAxis->GetBinWidth(i-1))*3.0;
+      diff = 1e4*diff*diff;
+    }
+    if (fgRegularizationType == kPowerLaw && i > 1+fgSkipBinsBegin) 
+    {
+
+      if (params[i] < 1e-8 || params[i-1] < 1e-8 || params[i-2] < 1e-8)
+       continue;
+      
+      if (fgUnfoldedAxis->GetBinWidth(i+1) < 1e-8 || fgUnfoldedAxis->GetBinWidth(i) < 1e-8 || fgUnfoldedAxis->GetBinWidth(i) < 1e-8)
+       continue;
+      
+      double middle = TMath::Power(params[i-1] / fgUnfoldedAxis->GetBinWidth(i),fgPowern);
+
+      if(middle>0) {
+       double right  = TMath::Power(params[i] / fgUnfoldedAxis->GetBinWidth(i+1),fgPowern)/middle;
+       
+       double left   = TMath::Power(params[i-2] / fgUnfoldedAxis->GetBinWidth(i-1),fgPowern)/middle;
+       
+       middle = 1.;
+       
+       Double_t der1 = (right - middle) / ((fgUnfoldedAxis->GetBinWidth(i+1) + fgUnfoldedAxis->GetBinWidth(i)) / 2);
+       Double_t der2 = (middle - left) / ((fgUnfoldedAxis->GetBinWidth(i) + fgUnfoldedAxis->GetBinWidth(i-1)) / 2);
+       
+       diff = (der1 - der2) * (der1 - der2);// / error;
+      }
+    }
 
-    Double_t diff = (der1 - der2) / middle;
+    if (fgRegularizationType == kLogLog && i > 1+fgSkipBinsBegin) 
+    {
+
+      if (params[i] < 1e-8 || params[i-1] < 1e-8 || params[i-2] < 1e-8)
+       continue;
 
-    penalty->SetBinContent(i-1, diff * diff);
+      Double_t right  = log(params[i] / (*fgEfficiency)(i) / fgUnfoldedAxis->GetBinWidth(i+1));
+      Double_t middle = log(params[i-1] / (*fgEfficiency)(i-1) / fgUnfoldedAxis->GetBinWidth(i));
+      Double_t left   = log(params[i-2] / (*fgEfficiency)(i-2) / fgUnfoldedAxis->GetBinWidth(i-1));
+      
+      Double_t der1 = (right - middle) / ( log(fgUnfoldedAxis->GetBinCenter(i+1)) - log(fgUnfoldedAxis->GetBinCenter(i)) );
+      Double_t der2 = (middle - left) /( log(fgUnfoldedAxis->GetBinCenter(i)) - log(fgUnfoldedAxis->GetBinCenter(i-1)) );
+      
+      double tmp = (log(fgUnfoldedAxis->GetBinCenter(i+1)) - log(fgUnfoldedAxis->GetBinCenter(i-1)))/2.;
+      Double_t dder = (der1-der2) / tmp;
+      
+      diff = dder * dder;
+    }
+    
+    penalty->SetBinContent(i, diff*fgRegularizationWeight);
     
     //Printf("%d %f %f %f %f", i-1, left, middle, right, diff);
   }
-  new TCanvas; penalty->DrawCopy(); gPad->SetLogy();
+  
+  return penalty;
 }
 
 //____________________________________________________________________
@@ -980,7 +1740,7 @@ Int_t AliUnfolding::UnfoldWithFunction(TH2* correlation, TH1* efficiency, TH1* m
   
   SetChi2Regularization(kNone, 0);
   
-  SetStaticVariables(correlation, measured);
+  SetStaticVariables(correlation, measured, efficiency);
   
   // Initialize TMinuit via generic fitter interface
   TVirtualFitter *minuit = TVirtualFitter::Fitter(0, fgFitFunction->GetNpar());
@@ -996,6 +1756,7 @@ Int_t AliUnfolding::UnfoldWithFunction(TH2* correlation, TH1* efficiency, TH1* m
   Double_t arglist[100];
   arglist[0] = 0;
   minuit->ExecuteCommand("SET PRINT", arglist, 1);
+  minuit->ExecuteCommand("SCAN", arglist, 0);
   minuit->ExecuteCommand("MIGRAD", arglist, 0);
   //minuit->ExecuteCommand("MINOS", arglist, 0);
 
@@ -1101,3 +1862,118 @@ void AliUnfolding::CreateOverflowBin(TH2* correlation, TH1* measured)
     correlation->DrawCopy("COLZ");
   }
 }
+
+Int_t AliUnfolding::UnfoldGetBias(TH2* correlation, TH1* efficiency, TH1* measured, TH1* initialConditions, TH1* result)
+{
+  // unfolds and assigns bias as errors with Eq. 19 of Cowan, "a survey of unfolding methods for particle physics"
+  // b_i = sum_j dmu_i/dn_j (nu_j - n_j) with nu_j as folded guess, n_j as data
+  // dmu_i/dn_j is found numerically by changing the bin content and re-unfolding
+  //
+  // return code: 0 (success) -1 (error: from Unfold(...) )
+
+  if (Unfold(correlation, efficiency, measured, initialConditions, result) != 0)
+    return -1;
+
+  TMatrixD matrix(fgMaxInput, fgMaxParams);
+  
+  TH1* newResult[4];
+  for (Int_t loop=0; loop<4; loop++)
+    newResult[loop] = (TH1*) result->Clone(Form("newresult_%d", loop));
+    
+  // change bin-by-bin and built matrix of effects
+  for (Int_t m=0; m<fgMaxInput; m++)
+  {
+    if (measured->GetBinContent(m+1) < 1)
+      continue;
+      
+    for (Int_t loop=0; loop<4; loop++)
+    {
+      //Printf("%d %d", i, loop);
+      Float_t factor = 1;
+      switch (loop)
+      {
+        case 0: factor = 0.5; break;
+        case 1: factor = -0.5; break;
+        case 2: factor = 1; break;
+        case 3: factor = -1; break;
+        default: return -1;
+      }
+      
+      TH1* measuredClone = (TH1*) measured->Clone("measuredClone");
+  
+      measuredClone->SetBinContent(m+1, measured->GetBinContent(m+1) + factor * TMath::Sqrt(measured->GetBinContent(m+1)));
+      //new TCanvas; measuredClone->Draw("COLZ");
+    
+      newResult[loop]->Reset();
+      Unfold(correlation, efficiency, measuredClone, measuredClone, newResult[loop]);
+      // WARNING if we leave here, then nothing is calculated
+        //return -1;
+      
+      delete measuredClone;
+    }
+    
+    for (Int_t t=0; t<fgMaxParams; t++)
+    {
+      // one value
+      //matrix(i, j) = (result->GetBinContent(j+1) - newResult->GetBinContent(j+1)) / TMath::Sqrt(changedHist->GetBinContent(1, i+1));
+      
+      // four values from the derivate (procedure taken from ROOT -- suggestion by Ruben)
+      // = 1/2D * [ 8 (f(D/2) - f(-D/2)) - (f(D)-f(-D)) ]/3
+      
+      /*
+      // check formula
+      measured->SetBinContent(1, m+1, 100);
+      newResult[0]->SetBinContent(t+1, 5 * 0.5 + 10);
+      newResult[1]->SetBinContent(t+1, 5 * -0.5 + 10);
+      newResult[2]->SetBinContent(t+1, 5 * 1 + 10);
+      newResult[3]->SetBinContent(t+1, 5 * -1 + 10);
+      */
+      
+      matrix(m, t) = 0.5 / TMath::Sqrt(measured->GetBinContent(m+1)) * 
+        ( 8. * (newResult[0]->GetBinContent(t+1) - newResult[1]->GetBinContent(t+1)) - 
+              (newResult[2]->GetBinContent(t+1) - newResult[3]->GetBinContent(t+1))
+        ) / 3;
+    }
+  }
+  
+  for (Int_t loop=0; loop<4; loop++)
+    delete newResult[loop];
+  
+  // ...calculate folded guess  
+  TH1* convoluted = (TH1*) measured->Clone("convoluted");
+  convoluted->Reset();
+  for (Int_t m=0; m<fgMaxInput; m++)
+  {
+    Float_t value = 0;
+    for (Int_t t = 0; t<fgMaxParams; t++)
+    {
+      Float_t tmp = correlation->GetBinContent(t+1, m+1) * result->GetBinContent(t+1);
+      if (efficiency)
+        tmp *= efficiency->GetBinContent(t+1);
+      value += tmp;
+    }
+    convoluted->SetBinContent(m+1, value);
+  }
+  
+  // rescale
+  convoluted->Scale(measured->Integral() / convoluted->Integral());
+  
+  //new TCanvas; convoluted->Draw(); measured->Draw("SAME"); measured->SetLineColor(2);
+  //return;
+  
+  // difference
+  convoluted->Add(measured, -1);
+  
+  // ...and the bias
+  for (Int_t t = 0; t<fgMaxParams; t++)
+  {
+    Double_t error = 0;
+    for (Int_t m=0; m<fgMaxInput; m++)
+      error += matrix(m, t) * convoluted->GetBinContent(m+1); 
+    result->SetBinError(t+1, error);
+  }
+  
+  //new TCanvas; result->Draw(); gPad->SetLogy();
+  
+  return 0;
+}