Updated AliUnfolding class. The regularization schemas have changed taking into accou...
[u/mrichter/AliRoot.git] / PWG0 / AliUnfolding.cxx
index f716551..b3dd89f 100644 (file)
 #include <TF1.h>
 
 TMatrixD* AliUnfolding::fgCorrelationMatrix = 0;
+TMatrixD* AliUnfolding::fgCorrelationMatrixSquared = 0;
 TMatrixD* AliUnfolding::fgCorrelationCovarianceMatrix = 0;
 TVectorD* AliUnfolding::fgCurrentESDVector = 0;
 TVectorD* AliUnfolding::fgEntropyAPriori = 0;
+TVectorD* AliUnfolding::fgEfficiency = 0;
+TVectorD* AliUnfolding::fgBinWidths = 0;
 
 TF1* AliUnfolding::fgFitFunction = 0;
 
@@ -44,13 +47,19 @@ 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::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;
+
 ClassImp(AliUnfolding)
 
 //____________________________________________________________________
@@ -97,6 +106,42 @@ 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 (fgBinWidths)
+  {
+    delete fgBinWidths;
+    fgBinWidths = 0;
+  }
+  
   Printf("AliUnfolding::SetNbins: Set %d measured bins and %d unfolded bins", nMeasured, nUnfolded);
 }
 
@@ -148,6 +193,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)
   {
@@ -181,19 +228,25 @@ Int_t AliUnfolding::Unfold(TH2* correlation, TH1* efficiency, TH1* measured, TH1
 }
 
 //____________________________________________________________________
-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 (!fgBinWidths)
+    fgBinWidths = new TVectorD(fgMaxParams);
+    
   fgCorrelationMatrix->Zero();
   fgCorrelationCovarianceMatrix->Zero();
   fgCurrentESDVector->Zero();
@@ -217,30 +270,49 @@ 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);
+    (*fgBinWidths)(i) = efficiency->GetXaxis()->GetBinWidth(i+1);
+  }
 }
 
 //____________________________________________________________________
@@ -252,15 +324,20 @@ 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];
 
   // disable any output (-1), unfortuantly we do not see warnings anymore then. Have to find another way...
@@ -276,6 +353,13 @@ Int_t AliUnfolding::UnfoldWithMinuit(TH2* correlation, TH1* efficiency, TH1* mea
   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 +369,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 = 1e100;
+  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;
@@ -322,6 +431,13 @@ Int_t AliUnfolding::UnfoldWithMinuit(TH2* correlation, TH1* efficiency, TH1* mea
   //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);
@@ -330,16 +446,16 @@ Int_t AliUnfolding::UnfoldWithMinuit(TH2* correlation, TH1* efficiency, TH1* mea
     Double_t error = minuit->GetParError(i) * results[i];
     
     if (efficiency)
-    {
+    {  
       if (efficiency->GetBinContent(i+1) > 0)
       {
-        value /= efficiency->GetBinContent(i+1);
-        error /= efficiency->GetBinContent(i+1);
+       value /= efficiency->GetBinContent(i+1);
+       error /= efficiency->GetBinContent(i+1);
       }
       else
       {
-        value = 0;
-        error = 0;
+       value = 0;
+       error = 0;
       }
     }
     
@@ -347,6 +463,10 @@ Int_t AliUnfolding::UnfoldWithMinuit(TH2* correlation, TH1* efficiency, TH1* mea
     result->SetBinError(i+1, error);
   }
   
+  fgCallCount = 0;
+  Chi2Function(dummy, 0, chi2, results, 0);
+  printf("AliUnfolding::UnfoldWithMinuit: Chi2 of final parameters is = %f\n", chi2);
+  
   if (fgDebug)
     DrawGuess(results);
 
@@ -382,6 +502,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 +537,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 +552,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 +601,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 +621,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 +651,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++)
   {
@@ -683,13 +812,13 @@ Double_t AliUnfolding::RegularizationPol0(TVectorD& params)
 
   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] / (*fgBinWidths)(i);
+    Double_t left   = params[i-1] / (*fgBinWidths)(i-1);
 
-    if (right != 0)
+    if (left != 0)
     {
-      Double_t diff = 1 - left / right;
-      chi2 += diff * diff;
+      Double_t diff = (right - left);
+      chi2 += diff * diff / left / (((*fgBinWidths)(i) + (*fgBinWidths)(i-1)) / 2);
     }
   }
 
@@ -710,16 +839,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] / (*fgBinWidths)(i);
+    Double_t middle = params[i-1] / (*fgBinWidths)(i-1);
+    Double_t left   = params[i-2] / (*fgBinWidths)(i-2);
 
-    Double_t diff = (der1 - der2) / middle;
+    Double_t der1 = (right - middle) / (((*fgBinWidths)(i) + (*fgBinWidths)(i-1)) / 2);
+    Double_t der2 = (middle - left) / (((*fgBinWidths)(i-1) + (*fgBinWidths)(i-2)) / 2);
 
-    chi2 += diff * diff;
+    //Double_t diff = (der1 - der2) / middle;
+    //chi2 += diff * diff;
+    chi2 += (der1 - der2) * (der1 - der2) / middle * (*fgBinWidths)(i-1);
   }
 
   return chi2;
@@ -736,22 +865,24 @@ Double_t AliUnfolding::RegularizationLog(TVectorD& params)
 
   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] / (*fgBinWidths)(i));
+    Double_t middle = log(params[i-1] / (*fgBinWidths)(i-1));
+    Double_t left   = log(params[i-2] / (*fgBinWidths)(i-2));
+    
+    Double_t der1 = (right - middle) / (((*fgBinWidths)(i) + (*fgBinWidths)(i-1)) / 2);
+    Double_t der2 = (middle - left) / (((*fgBinWidths)(i-1) + (*fgBinWidths)(i-2)) / 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;
 }
 
 //____________________________________________________________________
@@ -797,15 +928,47 @@ 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
+
+  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] / (*fgBinWidths)(i);
+    Double_t middle = params[i-1] / (*fgBinWidths)(i-1);
+    Double_t left   = params[i-2] / (*fgBinWidths)(i-2);
+    Double_t left2   = params[i-3] / (*fgBinWidths)(i-3);
+    Double_t left3   = params[i-4] / (*fgBinWidths)(i-4);
+    Double_t left4   = params[i-5] / (*fgBinWidths)(i-5);
+
+    //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 100.0 + chi2;
+  return chi2;
 }
 
 //____________________________________________________________________
@@ -815,9 +978,11 @@ 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);
+  TVectorD paramsVector(fgMaxParams);
   for (Int_t i=0; i<fgMaxParams; ++i)
     paramsVector[i] = params[i] * params[i];
 
@@ -831,6 +996,7 @@ 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;
   }
 
   //if (penaltyVal > 1e10)
@@ -844,9 +1010,27 @@ 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;
 
+  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;
+  }
   //measGuessVector.Print();
 
+#else
+  // old error calcuation using the error on the measured
   TVectorD copy(measGuessVector);
 
   // (Ad - m) W
@@ -859,16 +1043,60 @@ void AliUnfolding::Chi2Function(Int_t&, Double_t*, Double_t& chi2, Double_t *par
 
   //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
 
-  chi2 = chi2FromFit + penaltyVal;
+  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 += (*fgBinWidths)(i);
+  }
+  avg /= sum;
+  Float_t chi2Avg = 0; //(avg - 3.73) * (avg - 3.73) * 100;
 
-  static Int_t callCount = 0;
-  if ((callCount++ % 10000) == 0)
-    Printf("AliUnfolding::Chi2Function: Iteration %d: %f %f --> %f", callCount, chi2FromFit, penaltyVal, chi2);
+  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], avg, chi2FromFit, penaltyVal, notFoundEventsConstraint, chi2Avg, chi2);
+    //for (Int_t i=0; i<fgMaxInput; ++i)
+    //  Printf("%d: %f", i, measGuessVector(i) * copy(i) * 1e6);
+  }
 }
 
 //____________________________________________________________________
@@ -879,7 +1107,7 @@ void AliUnfolding::DrawGuess(Double_t *params)
   //
 
   // d
-  static TVectorD paramsVector(fgMaxParams);
+  TVectorD paramsVector(fgMaxParams);
   for (Int_t i=0; i<fgMaxParams; ++i)
     paramsVector[i] = params[i] * params[i];
 
@@ -913,32 +1141,86 @@ void AliUnfolding::DrawGuess(Double_t *params)
   new TCanvas; residuals->DrawCopy(); gPad->SetLogy();
 
   // draw penalty
-  if (fgRegularizationType != kPol1) {
-    Printf("Drawing not supported");
-    return;
-  }
+  TH1* penalty = GetPenaltyPlot(params);
+  
+  new TCanvas; penalty->DrawCopy(); gPad->SetLogy();
+}
 
-  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] = 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);
+
+  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] / (*fgBinWidths)(i);
+      Double_t left   = params[i-1] / (*fgBinWidths)(i-1);
 
-    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 / (((*fgBinWidths)(i) + (*fgBinWidths)(i-1)) / 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] / (*fgBinWidths)(i);
+      Double_t middle = params[i-1] / (*fgBinWidths)(i-1);
+      Double_t left   = params[i-2] / (*fgBinWidths)(i-2);
 
-    Double_t diff = (der1 - der2) / middle;
+      Double_t der1 = (right - middle) / (((*fgBinWidths)(i) + (*fgBinWidths)(i-1)) / 2);
+      Double_t der2 = (middle - left) / (((*fgBinWidths)(i-1) + (*fgBinWidths)(i-2)) / 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;
+    }
 
-    penalty->SetBinContent(i-1, diff * diff);
+    penalty->SetBinContent(i, diff);
     
     //Printf("%d %f %f %f %f", i-1, left, middle, right, diff);
   }
-  new TCanvas; penalty->DrawCopy(); gPad->SetLogy();
+  
+  return penalty;
 }
 
 //____________________________________________________________________
@@ -980,7 +1262,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());
@@ -1101,3 +1383,121 @@ 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];
+  
+  //matrix.Print();
+  //TH2* matrixHist = new TH2D(matrix); new TCanvas; matrixHist->Draw("BOX");
+  
+  // ...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;
+}