PWG0/AliUnfolding.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id: AliUnfolding.cxx 31168 2009-02-23 15:18:45Z jgrosseo $ */
  17
  18 // This class allows 1-dimensional unfolding.
  19 // Methods that are implemented are chi2 minimization and bayesian unfolding.
  20 //
  21 //  Author: Jan.Fiete.Grosse-Oetringhaus@cern.ch
  22
  23 #include "AliUnfolding.h"
  24 #include <TH1F.h>
  25 #include <TH2F.h>
  26 #include <TVirtualFitter.h>
  27 #include <TMath.h>
  28 #include <TCanvas.h>
  29 #include <TF1.h>
  30
  31 TMatrixD* AliUnfolding::fgCorrelationMatrix = 0;
  32 TMatrixD* AliUnfolding::fgCorrelationMatrixSquared = 0;
  33 TMatrixD* AliUnfolding::fgCorrelationCovarianceMatrix = 0;
  34 TVectorD* AliUnfolding::fgCurrentESDVector = 0;
  35 TVectorD* AliUnfolding::fgEntropyAPriori = 0;
  36 TVectorD* AliUnfolding::fgEfficiency = 0;
  37 TVectorD* AliUnfolding::fgBinWidths = 0;
  38
  39 TF1* AliUnfolding::fgFitFunction = 0;
  40
  41 AliUnfolding::MethodType AliUnfolding::fgMethodType = AliUnfolding::kInvalid;
  42 Int_t AliUnfolding::fgMaxInput  = -1;  // bins in measured histogram
  43 Int_t AliUnfolding::fgMaxParams = -1;  // bins in unfolded histogram = number of fit params
  44 Float_t AliUnfolding::fgOverflowBinLimit = -1;
  45
  46 AliUnfolding::RegularizationType AliUnfolding::fgRegularizationType = AliUnfolding::kPol1;
  47 Float_t AliUnfolding::fgRegularizationWeight = 10000;
  48 Int_t AliUnfolding::fgSkipBinsBegin = 0;
  49 Float_t AliUnfolding::fgMinuitStepSize = 0.1;                 // (usually not needed to be changed) step size in minimization
  50 Bool_t AliUnfolding::fgMinimumInitialValue = kFALSE;          // set all initial values at least to the smallest value among the initial values
  51 Float_t AliUnfolding::fgMinimumInitialValueFix = -1;
  52 Bool_t AliUnfolding::fgNormalizeInput = kFALSE;                  // normalize input spectrum
  53 Float_t AliUnfolding::fgNotFoundEvents = 0;
  54 Bool_t AliUnfolding::fgSkipBin0InChi2 = kFALSE;
  55
  56 Float_t AliUnfolding::fgBayesianSmoothing  = 1;           // smoothing parameter (0 = no smoothing)
  57 Int_t   AliUnfolding::fgBayesianIterations = 10;          // number of iterations in Bayesian method
  58
  59 Bool_t AliUnfolding::fgDebug = kFALSE;
  60
  61 Int_t AliUnfolding::fgCallCount = 0;
  62
  63 ClassImp(AliUnfolding)
  64
  65 //____________________________________________________________________
  66 void AliUnfolding::SetUnfoldingMethod(MethodType methodType)
  67 {
  68   // set unfolding method
  69   fgMethodType = methodType;
  70
  71   const char* name = 0;
  72   switch (methodType)
  73   {
  74     case kInvalid: name = "INVALID"; break;
  75     case kChi2Minimization: name = "Chi2 Minimization"; break;
  76     case kBayesian: name = "Bayesian unfolding"; break;
  77     case kFunction: name = "Functional fit"; break;
  78   }
  79   Printf("AliUnfolding::SetUnfoldingMethod: %s enabled.", name);
  80 }
  81
  82 //____________________________________________________________________
  83 void AliUnfolding::SetCreateOverflowBin(Float_t overflowBinLimit)
  84 {
  85   // enable the creation of a overflow bin that includes all statistics below the given limit
  86
  87   fgOverflowBinLimit = overflowBinLimit;
  88
  89   Printf("AliUnfolding::SetCreateOverflowBin: overflow bin limit set to %f", overflowBinLimit);
  90 }
  91
  92 //____________________________________________________________________
  93 void AliUnfolding::SetSkipBinsBegin(Int_t nBins)
  94 {
  95   // set number of skipped bins in regularization
  96
  97   fgSkipBinsBegin = nBins;
  98
  99   Printf("AliUnfolding::SetSkipBinsBegin: skipping %d bins at the beginning of the spectrum in the regularization.", fgSkipBinsBegin);
 100 }
 101
 102 //____________________________________________________________________
 103 void AliUnfolding::SetNbins(Int_t nMeasured, Int_t nUnfolded)
 104 {
 105   // set number of bins in the input (measured) distribution and in the unfolded distribution
 106   fgMaxInput = nMeasured;
 107   fgMaxParams = nUnfolded;
 108
 109   if (fgCorrelationMatrix)
 110   {
 111     delete fgCorrelationMatrix;
 112     fgCorrelationMatrix = 0;
 113   }
 114   if (fgCorrelationMatrixSquared)
 115   {
 116     fgCorrelationMatrixSquared = 0;
 117     delete fgCorrelationMatrixSquared;
 118   }
 119   if (fgCorrelationCovarianceMatrix)
 120   {
 121     delete fgCorrelationCovarianceMatrix;
 122     fgCorrelationCovarianceMatrix = 0;
 123   }
 124   if (fgCurrentESDVector)
 125   {
 126     delete fgCurrentESDVector;
 127     fgCurrentESDVector = 0;
 128   }
 129   if (fgEntropyAPriori)
 130   {
 131     delete fgEntropyAPriori;
 132     fgEntropyAPriori = 0;
 133   }
 134   if (fgEfficiency)
 135   {
 136     delete fgEfficiency;
 137     fgEfficiency = 0;
 138   }
 139   if (fgBinWidths)
 140   {
 141     delete fgBinWidths;
 142     fgBinWidths = 0;
 143   }
 144
 145   Printf("AliUnfolding::SetNbins: Set %d measured bins and %d unfolded bins", nMeasured, nUnfolded);
 146 }
 147
 148 //____________________________________________________________________
 149 void AliUnfolding::SetChi2Regularization(RegularizationType type, Float_t weight)
 150 {
 151   //
 152   // sets the parameters for chi2 minimization
 153   //
 154
 155   fgRegularizationType = type;
 156   fgRegularizationWeight = weight;
 157
 158   Printf("AliUnfolding::SetChi2Regularization --> Regularization set to %d with weight %f", (Int_t) type, weight);
 159 }
 160
 161 //____________________________________________________________________
 162 void AliUnfolding::SetBayesianParameters(Float_t smoothing, Int_t nIterations)
 163 {
 164   //
 165   // sets the parameters for Bayesian unfolding
 166   //
 167
 168   fgBayesianSmoothing = smoothing;
 169   fgBayesianIterations = nIterations;
 170
 171   Printf("AliUnfolding::SetBayesianParameters --> Paramaters set to %d iterations with smoothing %f", fgBayesianIterations, fgBayesianSmoothing);
 172 }
 173
 174 //____________________________________________________________________
 175 void AliUnfolding::SetFunction(TF1* function)
 176 {
 177   // set function for unfolding with a fit function
 178
 179   fgFitFunction = function;
 180
 181   Printf("AliUnfolding::SetFunction: Set fit function with %d parameters.", function->GetNpar());
 182 }
 183
 184 //____________________________________________________________________
 185 Int_t AliUnfolding::Unfold(TH2* correlation, TH1* efficiency, TH1* measured, TH1* initialConditions, TH1* result, Bool_t check)
 186 {
 187   // unfolds with unfolding method fgMethodType
 188   //
 189   // parameters:
 190   //  correlation: response matrix as measured vs. generated
 191   //  efficiency:  (optional) efficiency that is applied on the unfolded spectrum, i.e. it has to be in unfolded variables. If 0 no efficiency is applied.
 192   //  measured:    the measured spectrum
 193   //  initialConditions: (optional) initial conditions for the unfolding. if 0 the measured spectrum is used as initial conditions.
 194   //  result:      target for the unfolded result
 195   //  check:       depends on the unfolding method, see comments in specific functions
 196   //
 197   //  return code: see UnfoldWithMinuit/UnfoldWithBayesian/UnfoldWithFunction
 198
 199   if (fgMaxInput == -1)
 200   {
 201     Printf("AliUnfolding::Unfold: WARNING. Number of measured bins not set with SetNbins. Using number of bins in measured distribution");
 202     fgMaxInput = measured->GetNbinsX();
 203   }
 204   if (fgMaxParams == -1)
 205   {
 206     Printf("AliUnfolding::Unfold: WARNING. Number of unfolded bins not set with SetNbins. Using number of bins in measured distribution");
 207     fgMaxParams = measured->GetNbinsX();
 208   }
 209
 210   if (fgOverflowBinLimit > 0)
 211     CreateOverflowBin(correlation, measured);
 212
 213   switch (fgMethodType)
 214   {
 215     case kInvalid:
 216     {
 217       Printf("AliUnfolding::Unfold: ERROR: Unfolding method not set. Use SetUnfoldingMethod. Exiting...");
 218       return -1;
 219     }
 220     case kChi2Minimization:
 221       return UnfoldWithMinuit(correlation, efficiency, measured, initialConditions, result, check);
 222     case kBayesian:
 223       return UnfoldWithBayesian(correlation, efficiency, measured, initialConditions, result);
 224     case kFunction:
 225       return UnfoldWithFunction(correlation, efficiency, measured, initialConditions, result);
 226   }
 227   return -1;
 228 }
 229
 230 //____________________________________________________________________
 231 void AliUnfolding::SetStaticVariables(TH2* correlation, TH1* measured, TH1* efficiency)
 232 {
 233   // fill static variables needed for minuit fit
 234
 235   if (!fgCorrelationMatrix)
 236     fgCorrelationMatrix = new TMatrixD(fgMaxInput, fgMaxParams);
 237   if (!fgCorrelationMatrixSquared)
 238     fgCorrelationMatrixSquared = new TMatrixD(fgMaxInput, fgMaxParams);
 239   if (!fgCorrelationCovarianceMatrix)
 240     fgCorrelationCovarianceMatrix = new TMatrixD(fgMaxInput, fgMaxInput);
 241   if (!fgCurrentESDVector)
 242     fgCurrentESDVector = new TVectorD(fgMaxInput);
 243   if (!fgEntropyAPriori)
 244     fgEntropyAPriori = new TVectorD(fgMaxParams);
 245   if (!fgEfficiency)
 246     fgEfficiency = new TVectorD(fgMaxParams);
 247   if (!fgBinWidths)
 248     fgBinWidths = new TVectorD(fgMaxParams);
 249
 250   fgCorrelationMatrix->Zero();
 251   fgCorrelationCovarianceMatrix->Zero();
 252   fgCurrentESDVector->Zero();
 253   fgEntropyAPriori->Zero();
 254
 255   // normalize correction for given nPart
 256   for (Int_t i=1; i<=correlation->GetNbinsX(); ++i)
 257   {
 258     Double_t sum = correlation->Integral(i, i, 1, correlation->GetNbinsY());
 259     if (sum <= 0)
 260       continue;
 261     Float_t maxValue = 0;
 262     Int_t maxBin = -1;
 263     for (Int_t j=1; j<=correlation->GetNbinsY(); ++j)
 264     {
 265       // find most probably value
 266       if (maxValue < correlation->GetBinContent(i, j))
 267       {
 268         maxValue = correlation->GetBinContent(i, j);
 269         maxBin = j;
 270       }
 271
 272       // npart sum to 1
 273       correlation->SetBinContent(i, j, correlation->GetBinContent(i, j) / sum);// * correlation->GetXaxis()->GetBinWidth(i));
 274       correlation->SetBinError(i, j, correlation->GetBinError(i, j) / sum);
 275
 276       if (i <= fgMaxParams && j <= fgMaxInput)
 277       {
 278         (*fgCorrelationMatrix)(j-1, i-1) = correlation->GetBinContent(i, j);
 279         (*fgCorrelationMatrixSquared)(j-1, i-1) = correlation->GetBinContent(i, j) * correlation->GetBinContent(i, j);
 280       }
 281     }
 282
 283     //printf("MPV for Ntrue = %f is %f\n", fCurrentCorrelation->GetXaxis()->GetBinCenter(i), fCurrentCorrelation->GetYaxis()->GetBinCenter(maxBin));
 284   }
 285
 286   //normalize measured
 287   Float_t smallestError = 1;
 288   if (fgNormalizeInput)
 289   {
 290     Float_t sumMeasured = measured->Integral();
 291     measured->Scale(1.0 / sumMeasured);
 292     smallestError /= sumMeasured;
 293   }
 294
 295   for (Int_t i=0; i<fgMaxInput; ++i)
 296   {
 297     (*fgCurrentESDVector)[i] = measured->GetBinContent(i+1);
 298     if (measured->GetBinError(i+1) > 0)
 299     {
 300       (*fgCorrelationCovarianceMatrix)(i, i) = (Double_t) 1e-6 / measured->GetBinError(i+1) / measured->GetBinError(i+1);
 301     }
 302     else // in this case put error of 1, otherwise 0 bins are not added to the chi2...
 303       (*fgCorrelationCovarianceMatrix)(i, i) = (Double_t) 1e-6 / smallestError / smallestError;
 304
 305     if ((*fgCorrelationCovarianceMatrix)(i, i) > 1e7)
 306       (*fgCorrelationCovarianceMatrix)(i, i) = 0;
 307     //Printf("%d, %e", i, (*fgCorrelationCovarianceMatrix)(i, i));
 308   }
 309
 310   // efficiency is expected to match bin width of result
 311   for (Int_t i=0; i<fgMaxParams; ++i)
 312   {
 313     (*fgEfficiency)(i) = efficiency->GetBinContent(i+1);
 314     (*fgBinWidths)(i) = efficiency->GetXaxis()->GetBinWidth(i+1);
 315   }
 316 }
 317
 318 //____________________________________________________________________
 319 Int_t AliUnfolding::UnfoldWithMinuit(TH2* correlation, TH1* efficiency, TH1* measured, TH1* initialConditions, TH1* result, Bool_t check)
 320 {
 321   //
 322   // implementation of unfolding (internal function)
 323   //
 324   // unfolds <measured> using response from <correlation> and effiency <efficiency>
 325   // output is in <result>
 326   // <initialConditions> set the initial values for the minimization, if 0 <measured> is used
 327   //   negative values in initialConditions mean that the given parameter is fixed to the absolute of the value
 328   // if <check> is true no unfolding is made, instead only the chi2 without unfolding is printed
 329   //
 330   // returns minuit status (0 = success), (-1 when check was set)
 331   //
 332
 333   SetStaticVariables(correlation, measured, efficiency);
 334
 335   // Initialize TMinuit via generic fitter interface
 336   Int_t params = fgMaxParams;
 337   if (fgNotFoundEvents > 0)
 338     params++;
 339
 340   TVirtualFitter *minuit = TVirtualFitter::Fitter(0, params);
 341   Double_t arglist[100];
 342
 343   // disable any output (-1), unfortuantly we do not see warnings anymore then. Have to find another way...
 344   arglist[0] = 0;
 345   minuit->ExecuteCommand("SET PRINT", arglist, 1);
 346
 347   // however, enable warnings
 348   //minuit->ExecuteCommand("SET WAR", arglist, 0);
 349
 350   // set minimization function
 351   minuit->SetFCN(Chi2Function);
 352
 353   for (Int_t i=0; i<fgMaxParams; i++)
 354     (*fgEntropyAPriori)[i] = 1;
 355
 356   // set initial conditions as a-priori distribution for MRX regularization
 357   /*
 358   for (Int_t i=0; i<fgMaxParams; i++)
 359     if (initialConditions && initialConditions->GetBinContent(i+1) > 0)
 360       (*fgEntropyAPriori)[i] = initialConditions->GetBinContent(i+1);
 361   */
 362
 363   if (!initialConditions) {
 364     initialConditions = measured;
 365   } else {
 366     Printf("AliUnfolding::UnfoldWithMinuit: Using different initial conditions...");
 367     //new TCanvas; initialConditions->DrawCopy();
 368     if (fgNormalizeInput)
 369       initialConditions->Scale(1.0 / initialConditions->Integral());
 370   }
 371
 372   // extract minimum value from initial conditions (if we set a value to 0 it will stay 0)
 373   Float_t minValue = 1e35;
 374   if (fgMinimumInitialValueFix < 0)
 375   {
 376     for (Int_t i=0; i<fgMaxParams; ++i)
 377     {
 378       Int_t bin = initialConditions->GetXaxis()->FindBin(result->GetXaxis()->GetBinCenter(i+1));
 379       if (initialConditions->GetBinContent(bin) > 0)
 380         minValue = TMath::Min(minValue, (Float_t) initialConditions->GetBinContent(bin));
 381     }
 382   }
 383   else
 384     minValue = fgMinimumInitialValueFix;
 385
 386   Double_t* results = new Double_t[fgMaxParams+1];
 387   for (Int_t i=0; i<fgMaxParams; ++i)
 388   {
 389     Int_t bin = initialConditions->GetXaxis()->FindBin(result->GetXaxis()->GetBinCenter(i+1));
 390     results[i] = initialConditions->GetBinContent(bin);
 391
 392     Bool_t fix = kFALSE;
 393     if (results[i] < 0)
 394     {
 395       fix = kTRUE;
 396       results[i] = -results[i];
 397     }
 398
 399     if (!fix && fgMinimumInitialValue && results[i] < minValue)
 400       results[i] = minValue;
 401
 402     // minuit sees squared values to prevent it from going negative...
 403     results[i] = TMath::Sqrt(results[i]);
 404
 405     minuit->SetParameter(i, Form("param%d", i), results[i], (fix) ? 0 : fgMinuitStepSize, 0, 0);
 406   }
 407   if (fgNotFoundEvents > 0)
 408   {
 409     results[fgMaxParams] = efficiency->GetBinContent(1);
 410     minuit->SetParameter(fgMaxParams, "vtx0", results[fgMaxParams], fgMinuitStepSize / 100, 0.01, 0.80);
 411   }
 412
 413   Int_t dummy = 0;
 414   Double_t chi2 = 0;
 415   Chi2Function(dummy, 0, chi2, results, 0);
 416   printf("AliUnfolding::UnfoldWithMinuit: Chi2 of initial parameters is = %f\n", chi2);
 417
 418   if (check)
 419   {
 420     DrawGuess(results);
 421     delete[] results;
 422     return -1;
 423   }
 424
 425   // first param is number of iterations, second is precision....
 426   arglist[0] = 1e6;
 427   //arglist[1] = 1e-5;
 428   //minuit->ExecuteCommand("SCAN", arglist, 0);
 429   Int_t status = minuit->ExecuteCommand("MIGRAD", arglist, 1);
 430   Printf("AliUnfolding::UnfoldWithMinuit: MINUIT status is %d", status);
 431   //printf("!!!!!!!!!!!!!! MIGRAD finished: Starting MINOS !!!!!!!!!!!!!!");
 432   //minuit->ExecuteCommand("MINOS", arglist, 0);
 433
 434   if (fgNotFoundEvents > 0)
 435   {
 436     results[fgMaxParams] = minuit->GetParameter(fgMaxParams);
 437     Printf("Efficiency for bin 0 changed from %f to %f", efficiency->GetBinContent(1), results[fgMaxParams]);
 438     efficiency->SetBinContent(1, results[fgMaxParams]);
 439   }
 440
 441   for (Int_t i=0; i<fgMaxParams; ++i)
 442   {
 443     results[i] = minuit->GetParameter(i);
 444     Double_t value = results[i] * results[i];
 445     // error is : (relError) * (value) = (minuit->GetParError(i) / minuit->GetParameter(i)) * (minuit->GetParameter(i) * minuit->GetParameter(i))
 446     Double_t error = 0;
 447     if (TMath::IsNaN(minuit->GetParError(i)))
 448       Printf("WARNING: Parameter %d error is nan", i);
 449     else
 450       error = minuit->GetParError(i) * results[i];
 451
 452     if (efficiency)
 453     {
 454       if (efficiency->GetBinContent(i+1) > 0)
 455       {
 456         value /= efficiency->GetBinContent(i+1);
 457         error /= efficiency->GetBinContent(i+1);
 458       }
 459       else
 460       {
 461         value = 0;
 462         error = 0;
 463       }
 464     }
 465
 466     result->SetBinContent(i+1, value);
 467     result->SetBinError(i+1, error);
 468   }
 469
 470   fgCallCount = 0;
 471   Chi2Function(dummy, 0, chi2, results, 0);
 472   Printf("AliUnfolding::UnfoldWithMinuit: Chi2 of final parameters is = %f", chi2);
 473
 474   if (fgDebug)
 475     DrawGuess(results);
 476
 477   delete[] results;
 478
 479   return status;
 480 }
 481
 482 //____________________________________________________________________
 483 Int_t AliUnfolding::UnfoldWithBayesian(TH2* correlation, TH1* aEfficiency, TH1* measured, TH1* initialConditions, TH1* aResult)
 484 {
 485   //
 486   // unfolds a spectrum using the Bayesian method
 487   //
 488
 489   if (measured->Integral() <= 0)
 490   {
 491     Printf("AliUnfolding::UnfoldWithBayesian: ERROR: The measured spectrum is empty");
 492     return -1;
 493   }
 494
 495   const Int_t kStartBin = 0;
 496
 497   Int_t kMaxM = fgMaxInput;  //<= fCurrentCorrelation->GetNbinsY(); // max measured axis
 498   Int_t kMaxT = fgMaxParams; //<= fCurrentCorrelation->GetNbinsX(); // max true axis
 499
 500   // convergence limit: kMaxT * 0.001^2 = kMaxT * 1e-6 (e.g. 250 bins --> 2.5 e-4)
 501   const Double_t kConvergenceLimit = kMaxT * 1e-6;
 502
 503   // store information in arrays, to increase processing speed (~ factor 5)
 504   Double_t* measuredCopy = new Double_t[kMaxM];
 505   Double_t* measuredError = new Double_t[kMaxM];
 506   Double_t* prior = new Double_t[kMaxT];
 507   Double_t* result = new Double_t[kMaxT];
 508   Double_t* efficiency = new Double_t[kMaxT];
 509   Double_t* binWidths = new Double_t[kMaxT];
 510
 511   Double_t** response = new Double_t*[kMaxT];
 512   Double_t** inverseResponse = new Double_t*[kMaxT];
 513   for (Int_t i=0; i<kMaxT; i++)
 514   {
 515     response[i] = new Double_t[kMaxM];
 516     inverseResponse[i] = new Double_t[kMaxM];
 517   }
 518
 519   // for normalization
 520   Float_t measuredIntegral = measured->Integral();
 521   for (Int_t m=0; m<kMaxM; m++)
 522   {
 523     measuredCopy[m] = measured->GetBinContent(m+1) / measuredIntegral;
 524     measuredError[m] = measured->GetBinError(m+1) / measuredIntegral;
 525
 526     for (Int_t t=0; t<kMaxT; t++)
 527     {
 528       response[t][m] = correlation->GetBinContent(t+1, m+1);
 529       inverseResponse[t][m] = 0;
 530     }
 531   }
 532
 533   for (Int_t t=0; t<kMaxT; t++)
 534   {
 535     if (aEfficiency)
 536     {
 537       efficiency[t] = aEfficiency->GetBinContent(t+1);
 538     }
 539     else
 540       efficiency[t] = 1;
 541
 542     prior[t] = measuredCopy[t];
 543     result[t] = 0;
 544     binWidths[t] = aResult->GetXaxis()->GetBinWidth(t+1);
 545   }
 546
 547   // pick prior distribution
 548   if (initialConditions)
 549   {
 550     printf("Using different starting conditions...\n");
 551     // for normalization
 552     Float_t inputDistIntegral = initialConditions->Integral();
 553     for (Int_t i=0; i<kMaxT; i++)
 554       prior[i] = initialConditions->GetBinContent(i+1) / inputDistIntegral;
 555   }
 556
 557   //TH1F* convergence = new TH1F("convergence", "convergence", 200, 0.5, 200.5);
 558
 559   //new TCanvas;
 560   // unfold...
 561   for (Int_t i=0; i<fgBayesianIterations || fgBayesianIterations < 0; i++)
 562   {
 563     if (fgDebug)
 564       Printf("AliUnfolding::UnfoldWithBayesian: iteration %i", i);
 565
 566     // calculate IR from Bayes theorem
 567     // IR_ji = R_ij * prior_i / sum_k(R_kj * prior_k)
 568
 569     Double_t chi2Measured = 0;
 570     for (Int_t m=0; m<kMaxM; m++)
 571     {
 572       Float_t norm = 0;
 573       for (Int_t t = kStartBin; t<kMaxT; t++)
 574         norm += response[t][m] * prior[t];
 575
 576       // calc. chi2: (measured - response * prior) / error
 577       if (measuredError[m] > 0)
 578       {
 579         Double_t value = (measuredCopy[m] - norm) / measuredError[m];
 580         chi2Measured += value * value;
 581       }
 582
 583       if (norm > 0)
 584       {
 585         for (Int_t t = kStartBin; t<kMaxT; t++)
 586           inverseResponse[t][m] = response[t][m] * prior[t] / norm;
 587       }
 588       else
 589       {
 590         for (Int_t t = kStartBin; t<kMaxT; t++)
 591           inverseResponse[t][m] = 0;
 592       }
 593     }
 594     //Printf("chi2Measured of the last prior is %e", chi2Measured);
 595
 596     for (Int_t t = kStartBin; t<kMaxT; t++)
 597     {
 598       Float_t value = 0;
 599       for (Int_t m=0; m<kMaxM; m++)
 600         value += inverseResponse[t][m] * measuredCopy[m];
 601
 602       if (efficiency[t] > 0)
 603         result[t] = value / efficiency[t];
 604       else
 605         result[t] = 0;
 606     }
 607
 608     /*
 609     // draw intermediate result
 610     for (Int_t t=0; t<kMaxT; t++)
 611     {
 612       aResult->SetBinContent(t+1, result[t]);
 613     }
 614     aResult->SetMarkerStyle(24+i);
 615     aResult->SetMarkerColor(2);
 616     aResult->DrawCopy((i == 0) ? "P" : "PSAME");
 617     */
 618
 619     Double_t chi2LastIter = 0;
 620     // regularization (simple smoothing)
 621     for (Int_t t=kStartBin; t<kMaxT; t++)
 622     {
 623       Float_t newValue = 0;
 624
 625       // 0 bin excluded from smoothing
 626       if (t > kStartBin+2 && t<kMaxT-1)
 627       {
 628         Float_t average = (result[t-1] / binWidths[t-1] + result[t] / binWidths[t] + result[t+1] / binWidths[t+1]) / 3 * binWidths[t];
 629
 630         // weight the average with the regularization parameter
 631         newValue = (1 - fgBayesianSmoothing) * result[t] + fgBayesianSmoothing * average;
 632       }
 633       else
 634         newValue = result[t];
 635
 636       // calculate chi2 (change from last iteration)
 637       if (prior[t] > 1e-5)
 638       {
 639         Double_t diff = (prior[t] - newValue) / prior[t];
 640         chi2LastIter += diff * diff;
 641       }
 642
 643       prior[t] = newValue;
 644     }
 645     //printf("Chi2 of %d iteration = %e\n", i, chi2LastIter);
 646     //convergence->Fill(i+1, chi2LastIter);
 647
 648     if (fgBayesianIterations < 0 && chi2LastIter < kConvergenceLimit)
 649     {
 650       Printf("AliUnfolding::UnfoldWithBayesian: Stopped Bayesian unfolding after %d iterations at chi2(change since last iteration) of %e; chi2Measured of the last prior is %e", i, chi2LastIter, chi2Measured);
 651       break;
 652     }
 653   } // end of iterations
 654
 655   //new TCanvas; convergence->DrawCopy(); gPad->SetLogy();
 656   //delete convergence;
 657
 658   Float_t factor = 1;
 659   if (!fgNormalizeInput)
 660     factor = measuredIntegral;
 661   for (Int_t t=0; t<kMaxT; t++)
 662     aResult->SetBinContent(t+1, result[t] * factor);
 663
 664   delete[] measuredCopy;
 665   delete[] measuredError;
 666   delete[] prior;
 667   delete[] result;
 668   delete[] efficiency;
 669   delete[] binWidths;
 670
 671   for (Int_t i=0; i<kMaxT; i++)
 672   {
 673     delete[] response[i];
 674     delete[] inverseResponse[i];
 675   }
 676   delete[] response;
 677   delete[] inverseResponse;
 678
 679   return 0;
 680
 681   // ********
 682   // Calculate the covariance matrix, all arguments are taken from NIM,A362,487-498,1995
 683
 684   /*printf("Calculating covariance matrix. This may take some time...\n");
 685
 686   // check if this is the right one...
 687   TH1* sumHist = GetMultiplicityMC(inputRange, eventType)->ProjectionY("sumHist", 1, GetMultiplicityMC(inputRange, eventType)->GetNbinsX());
 688
 689   Int_t xBins = hInverseResponseBayes->GetNbinsX();
 690   Int_t yBins = hInverseResponseBayes->GetNbinsY();
 691
 692   // calculate "unfolding matrix" Mij
 693   Float_t matrixM[251][251];
 694   for (Int_t i=1; i<=xBins; i++)
 695   {
 696     for (Int_t j=1; j<=yBins; j++)
 697     {
 698       if (fCurrentEfficiency->GetBinContent(i) > 0)
 699         matrixM[i-1][j-1] = hInverseResponseBayes->GetBinContent(i, j) / fCurrentEfficiency->GetBinContent(i);
 700       else
 701         matrixM[i-1][j-1] = 0;
 702     }
 703   }
 704
 705   Float_t* vectorn = new Float_t[yBins];
 706   for (Int_t j=1; j<=yBins; j++)
 707     vectorn[j-1] = fCurrentESD->GetBinContent(j);
 708
 709   // first part of covariance matrix, depends on input distribution n(E)
 710   Float_t cov1[251][251];
 711
 712   Float_t nEvents = fCurrentESD->Integral(); // N
 713
 714   xBins = 20;
 715   yBins = 20;
 716
 717   for (Int_t k=0; k<xBins; k++)
 718   {
 719     printf("In Cov1: %d\n", k);
 720     for (Int_t l=0; l<yBins; l++)
 721     {
 722       cov1[k][l] = 0;
 723
 724       // sum_j Mkj Mlj n(Ej) * (1 - n(Ej) / N)
 725       for (Int_t j=0; j<yBins; j++)
 726         cov1[k][l] += matrixM[k][j] * matrixM[l][j] * vectorn[j]
 727           * (1.0 - vectorn[j] / nEvents);
 728
 729       // - sum_i,j (i != j) Mki Mlj n(Ei) n(Ej) / N
 730       for (Int_t i=0; i<yBins; i++)
 731         for (Int_t j=0; j<yBins; j++)
 732         {
 733           if (i == j)
 734             continue;
 735           cov1[k][l] -= matrixM[k][i] * matrixM[l][j] * vectorn[i]
 736             * vectorn[j] / nEvents;
 737          }
 738     }
 739   }
 740
 741   printf("Cov1 finished\n");
 742
 743   TH2F* cov = (TH2F*) hInverseResponseBayes->Clone("cov");
 744   cov->Reset();
 745
 746   for (Int_t i=1; i<=xBins; i++)
 747     for (Int_t j=1; j<=yBins; j++)
 748       cov->SetBinContent(i, j, cov1[i-1][j-1]);
 749
 750   new TCanvas;
 751   cov->Draw("COLZ");
 752
 753   // second part of covariance matrix, depends on response matrix
 754   Float_t cov2[251][251];
 755
 756   // Cov[P(Er|Cu), P(Es|Cu)] term
 757   Float_t covTerm[100][100][100];
 758   for (Int_t r=0; r<yBins; r++)
 759     for (Int_t u=0; u<xBins; u++)
 760       for (Int_t s=0; s<yBins; s++)
 761       {
 762         if (r == s)
 763           covTerm[r][u][s] = 1.0 / sumHist->GetBinContent(u+1) * hResponse->GetBinContent(u+1, r+1)
 764             * (1.0 - hResponse->GetBinContent(u+1, r+1));
 765         else
 766           covTerm[r][u][s] = - 1.0 / sumHist->GetBinContent(u+1) * hResponse->GetBinContent(u+1, r+1)
 767             * hResponse->GetBinContent(u+1, s+1);
 768       }
 769
 770   for (Int_t k=0; k<xBins; k++)
 771     for (Int_t l=0; l<yBins; l++)
 772     {
 773       cov2[k][l] = 0;
 774       printf("In Cov2: %d %d\n", k, l);
 775       for (Int_t i=0; i<yBins; i++)
 776         for (Int_t j=0; j<yBins; j++)
 777         {
 778           //printf("In Cov2: %d %d %d %d\n", k, l, i, j);
 779           // calculate Cov(Mki, Mlj) = sum{ru},{su} ...
 780           Float_t tmpCov = 0;
 781           for (Int_t r=0; r<yBins; r++)
 782             for (Int_t u=0; u<xBins; u++)
 783               for (Int_t s=0; s<yBins; s++)
 784               {
 785                 if (hResponse->GetBinContent(u+1, r+1) == 0 || hResponse->GetBinContent(u+1, s+1) == 0
 786                   || hResponse->GetBinContent(u+1, i+1) == 0)
 787                   continue;
 788
 789                 tmpCov += BayesCovarianceDerivate(matrixM, hResponse, fCurrentEfficiency, k, i, r, u)
 790                   * BayesCovarianceDerivate(matrixM, hResponse, fCurrentEfficiency, l, j, s, u)
 791                   * covTerm[r][u][s];
 792               }
 793
 794           cov2[k][l] += fCurrentESD->GetBinContent(i+1) * fCurrentESD->GetBinContent(j+1) * tmpCov;
 795         }
 796     }
 797
 798   printf("Cov2 finished\n");
 799
 800   for (Int_t i=1; i<=xBins; i++)
 801     for (Int_t j=1; j<=yBins; j++)
 802       cov->SetBinContent(i, j, cov1[i-1][j-1] + cov2[i-1][j-1]);
 803
 804   new TCanvas;
 805   cov->Draw("COLZ");*/
 806 }
 807
 808 //____________________________________________________________________
 809 Double_t AliUnfolding::RegularizationPol0(TVectorD& params)
 810 {
 811   // homogenity term for minuit fitting
 812   // pure function of the parameters
 813   // prefers constant function (pol0)
 814
 815   Double_t chi2 = 0;
 816
 817   for (Int_t i=1+fgSkipBinsBegin; i<fgMaxParams; ++i)
 818   {
 819     Double_t right  = params[i] / (*fgBinWidths)(i);
 820     Double_t left   = params[i-1] / (*fgBinWidths)(i-1);
 821
 822     if (left != 0)
 823     {
 824       Double_t diff = (right - left);
 825       chi2 += diff * diff / left / (((*fgBinWidths)(i) + (*fgBinWidths)(i-1)) / 2);
 826     }
 827   }
 828
 829   return chi2 / 100.0;
 830 }
 831
 832 //____________________________________________________________________
 833 Double_t AliUnfolding::RegularizationPol1(TVectorD& params)
 834 {
 835   // homogenity term for minuit fitting
 836   // pure function of the parameters
 837   // prefers linear function (pol1)
 838
 839   Double_t chi2 = 0;
 840
 841   for (Int_t i=2+fgSkipBinsBegin; i<fgMaxParams; ++i)
 842   {
 843     if (params[i-1] == 0)
 844       continue;
 845
 846     Double_t right  = params[i] / (*fgBinWidths)(i);
 847     Double_t middle = params[i-1] / (*fgBinWidths)(i-1);
 848     Double_t left   = params[i-2] / (*fgBinWidths)(i-2);
 849
 850     Double_t der1 = (right - middle) / (((*fgBinWidths)(i) + (*fgBinWidths)(i-1)) / 2);
 851     Double_t der2 = (middle - left) / (((*fgBinWidths)(i-1) + (*fgBinWidths)(i-2)) / 2);
 852
 853     //Double_t diff = (der1 - der2) / middle;
 854     //chi2 += diff * diff;
 855     chi2 += (der1 - der2) * (der1 - der2) / middle * (*fgBinWidths)(i-1);
 856   }
 857
 858   return chi2;
 859 }
 860
 861 //____________________________________________________________________
 862 Double_t AliUnfolding::RegularizationLog(TVectorD& params)
 863 {
 864   // homogenity term for minuit fitting
 865   // pure function of the parameters
 866   // prefers linear function (pol1)
 867
 868   Double_t chi2 = 0;
 869
 870   for (Int_t i=2+fgSkipBinsBegin; i<fgMaxParams; ++i)
 871   {
 872     if (params[i-1] == 0 || params[i] == 0 || params[i-2] == 0)
 873      continue;
 874
 875     Double_t right  = log(params[i] / (*fgBinWidths)(i));
 876     Double_t middle = log(params[i-1] / (*fgBinWidths)(i-1));
 877     Double_t left   = log(params[i-2] / (*fgBinWidths)(i-2));
 878
 879     Double_t der1 = (right - middle) / (((*fgBinWidths)(i) + (*fgBinWidths)(i-1)) / 2);
 880     Double_t der2 = (middle - left) / (((*fgBinWidths)(i-1) + (*fgBinWidths)(i-2)) / 2);
 881
 882     //Double_t error = 1. / params[i] + 4. / params[i-1] + 1. / params[i-2];
 883
 884     //if (fgCallCount == 0)
 885     //  Printf("%d %f %f", i, (der1 - der2) * (der1 - der2), error);
 886     chi2 += (der1 - der2) * (der1 - der2);// / error;
 887   }
 888
 889   return chi2;
 890 }
 891
 892 //____________________________________________________________________
 893 Double_t AliUnfolding::RegularizationTotalCurvature(TVectorD& params)
 894 {
 895   // homogenity term for minuit fitting
 896   // pure function of the parameters
 897   // minimizes the total curvature (from Unfolding Methods In High-Energy Physics Experiments,
 898   // V. Blobel (Hamburg U.) . DESY 84/118, Dec 1984. 40pp.
 899
 900   Double_t chi2 = 0;
 901
 902   for (Int_t i=2+fgSkipBinsBegin; i<fgMaxParams; ++i)
 903   {
 904     Double_t right  = params[i];
 905     Double_t middle = params[i-1];
 906     Double_t left   = params[i-2];
 907
 908     Double_t der1 = (right - middle);
 909     Double_t der2 = (middle - left);
 910
 911     Double_t diff = (der1 - der2);
 912
 913     chi2 += diff * diff;
 914   }
 915
 916   return chi2 * 1e4;
 917 }
 918
 919 //____________________________________________________________________
 920 Double_t AliUnfolding::RegularizationEntropy(TVectorD& params)
 921 {
 922   // homogenity term for minuit fitting
 923   // pure function of the parameters
 924   // calculates entropy, from
 925   // The method of reduced cross-entropy (M. Schmelling 1993)
 926
 927   Double_t paramSum = 0;
 928
 929   for (Int_t i=fgSkipBinsBegin; i<fgMaxParams; ++i)
 930     paramSum += params[i];
 931
 932   Double_t chi2 = 0;
 933   for (Int_t i=fgSkipBinsBegin; i<fgMaxParams; ++i)
 934   {
 935     Double_t tmp = params[i] / paramSum;
 936     //Double_t tmp = params[i];
 937     if (tmp > 0 && (*fgEntropyAPriori)[i] > 0)
 938     {
 939       chi2 += tmp * TMath::Log(tmp / (*fgEntropyAPriori)[i]);
 940     }
 941     else
 942       chi2 += 100;
 943   }
 944
 945   return -chi2;
 946 }
 947
 948 //____________________________________________________________________
 949 Double_t AliUnfolding::RegularizationRatio(TVectorD& params)
 950 {
 951   // homogenity term for minuit fitting
 952   // pure function of the parameters
 953
 954   Double_t chi2 = 0;
 955
 956   for (Int_t i=5+fgSkipBinsBegin; i<fgMaxParams; ++i)
 957   {
 958     if (params[i-1] == 0 || params[i] == 0)
 959       continue;
 960
 961     Double_t right  = params[i] / (*fgBinWidths)(i);
 962     Double_t middle = params[i-1] / (*fgBinWidths)(i-1);
 963     Double_t left   = params[i-2] / (*fgBinWidths)(i-2);
 964     Double_t left2   = params[i-3] / (*fgBinWidths)(i-3);
 965     Double_t left3   = params[i-4] / (*fgBinWidths)(i-4);
 966     Double_t left4   = params[i-5] / (*fgBinWidths)(i-5);
 967
 968     //Double_t diff = left / middle - middle / right;
 969     //Double_t diff = 2 * left / middle - middle / right - left2 / left;
 970     Double_t diff = 4 * left2 / left - middle / right - left / middle - left3 / left2 - left4 / left3;
 971
 972     chi2 += diff * diff;// / middle;
 973   }
 974
 975   return chi2;
 976 }
 977
 978 //____________________________________________________________________
 979 void AliUnfolding::Chi2Function(Int_t&, Double_t*, Double_t& chi2, Double_t *params, Int_t)
 980 {
 981   //
 982   // fit function for minuit
 983   // does: (m - Ad)W(m - Ad) where m = measured, A correlation matrix, d = guess, W = covariance matrix
 984   //
 985
 986   // TODO use static members for the variables here to speed up processing (no construction/deconstruction)
 987
 988   // d
 989   TVectorD paramsVector(fgMaxParams);
 990   for (Int_t i=0; i<fgMaxParams; ++i)
 991     paramsVector[i] = params[i] * params[i];
 992
 993   // calculate penalty factor
 994   Double_t penaltyVal = 0;
 995   switch (fgRegularizationType)
 996   {
 997     case kNone:       break;
 998     case kPol0:       penaltyVal = RegularizationPol0(paramsVector); break;
 999     case kPol1:       penaltyVal = RegularizationPol1(paramsVector); break;
1000     case kCurvature:  penaltyVal = RegularizationTotalCurvature(paramsVector); break;
1001     case kEntropy:    penaltyVal = RegularizationEntropy(paramsVector); break;
1002     case kLog:        penaltyVal = RegularizationLog(paramsVector); break;
1003     case kRatio:      penaltyVal = RegularizationRatio(paramsVector); break;
1004   }
1005
1006   //if (penaltyVal > 1e10)
1007   //  paramsVector2.Print();
1008
1009   penaltyVal *= fgRegularizationWeight;
1010
1011   // Ad
1012   TVectorD measGuessVector(fgMaxInput);
1013   measGuessVector = (*fgCorrelationMatrix) * paramsVector;
1014
1015   // Ad - m
1016   measGuessVector -= (*fgCurrentESDVector);
1017
1018 #if 0
1019   // new error calcuation using error on the guess
1020
1021   // error from guess
1022   TVectorD errorGuessVector(fgMaxInput);
1023   //errorGuessVector = (*fgCorrelationMatrixSquared) * paramsVector;
1024   errorGuessVector = (*fgCorrelationMatrix) * paramsVector;
1025
1026   Double_t chi2FromFit = 0;
1027   for (Int_t i=0; i<fgMaxInput; ++i)
1028   {
1029     Float_t error = 1;
1030     if (errorGuessVector(i) > 0)
1031       error = errorGuessVector(i);
1032     chi2FromFit += measGuessVector(i) * measGuessVector(i) / error;
1033   }
1034   //measGuessVector.Print();
1035
1036 #else
1037   // old error calcuation using the error on the measured
1038   TVectorD copy(measGuessVector);
1039
1040   // (Ad - m) W
1041   // this step can be optimized because currently only the diagonal elements of fgCorrelationCovarianceMatrix are used
1042   // normal way is like this:
1043   // measGuessVector *= (*fgCorrelationCovarianceMatrix);
1044   // optimized way like this:
1045   for (Int_t i=0; i<fgMaxInput; ++i)
1046     measGuessVector[i] *= (*fgCorrelationCovarianceMatrix)(i, i);
1047
1048   //measGuessVector.Print();
1049
1050   if (fgSkipBin0InChi2)
1051     measGuessVector[0] = 0;
1052
1053   // (Ad - m) W (Ad - m)
1054   // the factor 1e6 prevents that a very small number (measGuessVector[i]) is multiplied with a very
1055   // big number ((*fgCorrelationCovarianceMatrix)(i, i)) (see UnfoldWithMinuit)
1056   Double_t chi2FromFit = measGuessVector * copy * 1e6;
1057 #endif
1058
1059   Double_t notFoundEventsConstraint = 0;
1060   Double_t currentNotFoundEvents = 0;
1061   Double_t errorNotFoundEvents = 0;
1062
1063   if (fgNotFoundEvents > 0)
1064   {
1065     for (Int_t i=0; i<fgMaxParams; ++i)
1066     {
1067       Double_t eff = (1.0 / (*fgEfficiency)(i) - 1);
1068       if (i == 0)
1069         eff = (1.0 / params[fgMaxParams] - 1);
1070       currentNotFoundEvents += eff * paramsVector(i);
1071       errorNotFoundEvents += eff * eff * paramsVector(i); // error due to guess (paramsVector)
1072       if (i <= 3)
1073         errorNotFoundEvents += (eff * 0.03) * (eff * 0.03) * paramsVector(i) * paramsVector(i); // error on eff
1074       //if ((fgCallCount % 10000) == 0)
1075       //  Printf("%d %f %f %f", i, (*fgEfficiency)(i), paramsVector(i), currentNotFoundEvents);
1076     }
1077     errorNotFoundEvents += fgNotFoundEvents;
1078     // TODO add error on background, background estimate
1079
1080     notFoundEventsConstraint = (currentNotFoundEvents - fgNotFoundEvents) * (currentNotFoundEvents - fgNotFoundEvents) / errorNotFoundEvents;
1081   }
1082
1083   Float_t avg = 0;
1084   Float_t sum = 0;
1085   Float_t currentMult = 0;
1086   // try to match dn/deta
1087   for (Int_t i=0; i<fgMaxParams; ++i)
1088   {
1089     avg += paramsVector(i) * currentMult;
1090     sum += paramsVector(i);
1091     currentMult += (*fgBinWidths)(i);
1092   }
1093   avg /= sum;
1094   Float_t chi2Avg = 0; //(avg - 3.73) * (avg - 3.73) * 100;
1095
1096   chi2 = chi2FromFit + penaltyVal + notFoundEventsConstraint + chi2Avg;
1097
1098   if ((fgCallCount++ % 1000) == 0)
1099   {
1100     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);
1101     //for (Int_t i=0; i<fgMaxInput; ++i)
1102     //  Printf("%d: %f", i, measGuessVector(i) * copy(i) * 1e6);
1103   }
1104 }
1105
1106 //____________________________________________________________________
1107 void AliUnfolding::DrawGuess(Double_t *params)
1108 {
1109   //
1110   // draws residuals of solution suggested by params and effect of regularization
1111   //
1112
1113   // d
1114   TVectorD paramsVector(fgMaxParams);
1115   for (Int_t i=0; i<fgMaxParams; ++i)
1116     paramsVector[i] = params[i] * params[i];
1117
1118   // Ad
1119   TVectorD measGuessVector(fgMaxInput);
1120   measGuessVector = (*fgCorrelationMatrix) * paramsVector;
1121
1122   // Ad - m
1123   measGuessVector -= (*fgCurrentESDVector);
1124
1125   TVectorD copy(measGuessVector);
1126   //copy.Print();
1127
1128   // (Ad - m) W
1129   // this step can be optimized because currently only the diagonal elements of fgCorrelationCovarianceMatrix are used
1130   // normal way is like this:
1131   // measGuessVector *= (*fgCorrelationCovarianceMatrix);
1132   // optimized way like this:
1133   for (Int_t i=0; i<fgMaxInput; ++i)
1134     measGuessVector[i] *= (*fgCorrelationCovarianceMatrix)(i, i);
1135
1136   // (Ad - m) W (Ad - m)
1137   // the factor 1e6 prevents that a very small number (measGuessVector[i]) is multiplied with a very
1138   // big number ((*fgCorrelationCovarianceMatrix)(i, i)) (see ApplyMinuitFit)
1139   //Double_t chi2FromFit = measGuessVector * copy * 1e6;
1140
1141   // draw residuals
1142   TH1* residuals = new TH1F("residuals", "residuals", fgMaxInput, -0.5, fgMaxInput - 0.5);
1143   for (Int_t i=0; i<fgMaxInput; ++i)
1144     residuals->SetBinContent(i+1, measGuessVector(i) * copy(i) * 1e6);
1145   new TCanvas; residuals->DrawCopy(); gPad->SetLogy();
1146
1147   // draw penalty
1148   TH1* penalty = GetPenaltyPlot(params);
1149
1150   new TCanvas; penalty->DrawCopy(); gPad->SetLogy();
1151 }
1152
1153 //____________________________________________________________________
1154 TH1* AliUnfolding::GetPenaltyPlot(TH1* corrected)
1155 {
1156   // draws the penalty factors as function of multiplicity of the current selected regularization
1157
1158   Double_t* params = new Double_t[fgMaxParams];
1159   for (Int_t i=0; i<TMath::Min(fgMaxParams, corrected->GetNbinsX()); i++)
1160     params[i] = corrected->GetBinContent(i+1);
1161
1162   TH1* penalty = GetPenaltyPlot(params);
1163
1164   delete[] params;
1165
1166   return penalty;
1167 }
1168
1169 //____________________________________________________________________
1170 TH1* AliUnfolding::GetPenaltyPlot(Double_t* params)
1171 {
1172   // draws the penalty factors as function of multiplicity of the current selected regularization
1173
1174   TH1* penalty = new TH1F("penalty", ";unfolded multiplicity;penalty factor", fgMaxParams, -0.5, fgMaxParams - 0.5);
1175
1176   for (Int_t i=1+fgSkipBinsBegin; i<fgMaxParams; ++i)
1177   {
1178     Double_t diff = 0;
1179     if (fgRegularizationType == kPol0)
1180     {
1181       Double_t right  = params[i] / (*fgBinWidths)(i);
1182       Double_t left   = params[i-1] / (*fgBinWidths)(i-1);
1183
1184       if (left != 0)
1185       {
1186         Double_t diffTmp = (right - left);
1187         diff = diffTmp * diffTmp / left / (((*fgBinWidths)(i) + (*fgBinWidths)(i-1)) / 2) / 100;
1188       }
1189     }
1190     if (fgRegularizationType == kPol1 && i > 1+fgSkipBinsBegin)
1191     {
1192       if (params[i-1] == 0)
1193         continue;
1194
1195       Double_t right  = params[i] / (*fgBinWidths)(i);
1196       Double_t middle = params[i-1] / (*fgBinWidths)(i-1);
1197       Double_t left   = params[i-2] / (*fgBinWidths)(i-2);
1198
1199       Double_t der1 = (right - middle) / (((*fgBinWidths)(i) + (*fgBinWidths)(i-1)) / 2);
1200       Double_t der2 = (middle - left) / (((*fgBinWidths)(i-1) + (*fgBinWidths)(i-2)) / 2);
1201
1202       diff = (der1 - der2) * (der1 - der2) / middle;
1203     }
1204     if (fgRegularizationType == kLog && i > 1+fgSkipBinsBegin)
1205     {
1206       if (params[i-1] == 0)
1207         continue;
1208
1209       Double_t right  = log(params[i]);
1210       Double_t middle = log(params[i-1]);
1211       Double_t left   = log(params[i-2]);
1212
1213       Double_t der1 = (right - middle);
1214       Double_t der2 = (middle - left);
1215
1216       //Double_t error = 1. / params[i] + 4. / params[i-1] + 1. / params[i-2];
1217       //Printf("%d %f %f", i, (der1 - der2) * (der1 - der2), error);
1218
1219       diff = (der1 - der2) * (der1 - der2);// / error;
1220     }
1221
1222     penalty->SetBinContent(i, diff);
1223
1224     //Printf("%d %f %f %f %f", i-1, left, middle, right, diff);
1225   }
1226
1227   return penalty;
1228 }
1229
1230 //____________________________________________________________________
1231 void AliUnfolding::TF1Function(Int_t& unused1, Double_t* unused2, Double_t& chi2, Double_t *params, Int_t unused3)
1232 {
1233   //
1234   // fit function for minuit
1235   // uses the TF1 stored in fgFitFunction
1236   //
1237
1238   for (Int_t i=0; i<fgFitFunction->GetNpar(); i++)
1239     fgFitFunction->SetParameter(i, params[i]);
1240
1241   Double_t* params2 = new Double_t[fgMaxParams];
1242
1243   for (Int_t i=0; i<fgMaxParams; ++i)
1244     params2[i] = fgFitFunction->Eval(i);
1245
1246   Chi2Function(unused1, unused2, chi2, params2, unused3);
1247
1248   delete[] params2;
1249
1250   if (fgDebug)
1251     Printf("%f", chi2);
1252 }
1253
1254 //____________________________________________________________________
1255 Int_t AliUnfolding::UnfoldWithFunction(TH2* correlation, TH1* efficiency, TH1* measured, TH1* /* initialConditions */, TH1* aResult)
1256 {
1257   //
1258   // correct spectrum using minuit chi2 method applying a functional fit
1259   //
1260
1261   if (!fgFitFunction)
1262   {
1263     Printf("AliUnfolding::UnfoldWithFunction: ERROR fit function not set. Exiting.");
1264     return -1;
1265   }
1266
1267   SetChi2Regularization(kNone, 0);
1268
1269   SetStaticVariables(correlation, measured, efficiency);
1270
1271   // Initialize TMinuit via generic fitter interface
1272   TVirtualFitter *minuit = TVirtualFitter::Fitter(0, fgFitFunction->GetNpar());
1273
1274   minuit->SetFCN(TF1Function);
1275   for (Int_t i=0; i<fgFitFunction->GetNpar(); i++)
1276   {
1277     Double_t lower, upper;
1278     fgFitFunction->GetParLimits(i, lower, upper);
1279     minuit->SetParameter(i, Form("param%d",i), fgFitFunction->GetParameter(i), fgMinuitStepSize, lower, upper);
1280   }
1281
1282   Double_t arglist[100];
1283   arglist[0] = 0;
1284   minuit->ExecuteCommand("SET PRINT", arglist, 1);
1285   minuit->ExecuteCommand("MIGRAD", arglist, 0);
1286   //minuit->ExecuteCommand("MINOS", arglist, 0);
1287
1288   for (Int_t i=0; i<fgFitFunction->GetNpar(); i++)
1289     fgFitFunction->SetParameter(i, minuit->GetParameter(i));
1290
1291   for (Int_t i=0; i<fgMaxParams; ++i)
1292   {
1293     Double_t value = fgFitFunction->Eval(i);
1294     if (fgDebug)
1295       Printf("%d : %f", i, value);
1296
1297     if (efficiency)
1298     {
1299       if (efficiency->GetBinContent(i+1) > 0)
1300       {
1301         value /= efficiency->GetBinContent(i+1);
1302       }
1303       else
1304         value = 0;
1305     }
1306     aResult->SetBinContent(i+1, value);
1307     aResult->SetBinError(i+1, 0);
1308   }
1309
1310   return 0;
1311 }
1312
1313 //____________________________________________________________________
1314 void AliUnfolding::CreateOverflowBin(TH2* correlation, TH1* measured)
1315 {
1316   // Finds the first bin where the content is below fgStatLimit and combines all values for this bin and larger bins
1317   // The same limit is applied to the correlation
1318
1319   Int_t lastBin = 0;
1320   for (Int_t i=1; i<=measured->GetNbinsX(); ++i)
1321   {
1322     if (measured->GetBinContent(i) <= fgOverflowBinLimit)
1323     {
1324       lastBin = i;
1325       break;
1326     }
1327   }
1328
1329   if (lastBin == 0)
1330   {
1331     Printf("AliUnfolding::CreateOverflowBin: WARNING: First bin is already below limit of %f", fgOverflowBinLimit);
1332     return;
1333   }
1334
1335   Printf("AliUnfolding::CreateOverflowBin: Bin limit in measured spectrum determined to be %d", lastBin);
1336
1337   TCanvas* canvas = 0;
1338
1339   if (fgDebug)
1340   {
1341     canvas = new TCanvas("StatSolution", "StatSolution", 1000, 800);
1342     canvas->Divide(2, 2);
1343
1344     canvas->cd(1);
1345     measured->SetStats(kFALSE);
1346     measured->DrawCopy();
1347     gPad->SetLogy();
1348
1349     canvas->cd(2);
1350     correlation->SetStats(kFALSE);
1351     correlation->DrawCopy("COLZ");
1352   }
1353
1354   measured->SetBinContent(lastBin, measured->Integral(lastBin, measured->GetNbinsX()));
1355   for (Int_t i=lastBin+1; i<=measured->GetNbinsX(); ++i)
1356   {
1357     measured->SetBinContent(i, 0);
1358     measured->SetBinError(i, 0);
1359   }
1360   // the error is set to sqrt(N), better solution possible?, sum of relative errors of all contributions???
1361   measured->SetBinError(lastBin, TMath::Sqrt(measured->GetBinContent(lastBin)));
1362
1363   Printf("AliUnfolding::CreateOverflowBin: This bin has now %f +- %f entries", measured->GetBinContent(lastBin), measured->GetBinError(lastBin));
1364
1365   for (Int_t i=1; i<=correlation->GetNbinsX(); ++i)
1366   {
1367     correlation->SetBinContent(i, lastBin, correlation->Integral(i, i, lastBin, correlation->GetNbinsY()));
1368     // the error is set to sqrt(N), better solution possible?, sum of relative errors of all contributions???
1369     correlation->SetBinError(i, lastBin, TMath::Sqrt(correlation->GetBinContent(i, lastBin)));
1370
1371     for (Int_t j=lastBin+1; j<=correlation->GetNbinsY(); ++j)
1372     {
1373       correlation->SetBinContent(i, j, 0);
1374       correlation->SetBinError(i, j, 0);
1375     }
1376   }
1377
1378   Printf("AliUnfolding::CreateOverflowBin: Adjusted correlation matrix!");
1379
1380   if (canvas)
1381   {
1382     canvas->cd(3);
1383     measured->DrawCopy();
1384     gPad->SetLogy();
1385
1386     canvas->cd(4);
1387     correlation->DrawCopy("COLZ");
1388   }
1389 }
1390
1391 Int_t AliUnfolding::UnfoldGetBias(TH2* correlation, TH1* efficiency, TH1* measured, TH1* initialConditions, TH1* result)
1392 {
1393   // unfolds and assigns bias as errors with Eq. 19 of Cowan, "a survey of unfolding methods for particle physics"
1394   // b_i = sum_j dmu_i/dn_j (nu_j - n_j) with nu_j as folded guess, n_j as data
1395   // dmu_i/dn_j is found numerically by changing the bin content and re-unfolding
1396   //
1397   // return code: 0 (success) -1 (error: from Unfold(...) )
1398
1399   if (Unfold(correlation, efficiency, measured, initialConditions, result) != 0)
1400     return -1;
1401
1402   TMatrixD matrix(fgMaxInput, fgMaxParams);
1403
1404   TH1* newResult[4];
1405   for (Int_t loop=0; loop<4; loop++)
1406     newResult[loop] = (TH1*) result->Clone(Form("newresult_%d", loop));
1407
1408   // change bin-by-bin and built matrix of effects
1409   for (Int_t m=0; m<fgMaxInput; m++)
1410   {
1411     if (measured->GetBinContent(m+1) < 1)
1412       continue;
1413
1414     for (Int_t loop=0; loop<4; loop++)
1415     {
1416       //Printf("%d %d", i, loop);
1417       Float_t factor = 1;
1418       switch (loop)
1419       {
1420         case 0: factor = 0.5; break;
1421         case 1: factor = -0.5; break;
1422         case 2: factor = 1; break;
1423         case 3: factor = -1; break;
1424         default: return -1;
1425       }
1426
1427       TH1* measuredClone = (TH1*) measured->Clone("measuredClone");
1428
1429       measuredClone->SetBinContent(m+1, measured->GetBinContent(m+1) + factor * TMath::Sqrt(measured->GetBinContent(m+1)));
1430       //new TCanvas; measuredClone->Draw("COLZ");
1431
1432       newResult[loop]->Reset();
1433       Unfold(correlation, efficiency, measuredClone, measuredClone, newResult[loop]);
1434       // WARNING if we leave here, then nothing is calculated
1435         //return -1;
1436
1437       delete measuredClone;
1438     }
1439
1440     for (Int_t t=0; t<fgMaxParams; t++)
1441     {
1442       // one value
1443       //matrix(i, j) = (result->GetBinContent(j+1) - newResult->GetBinContent(j+1)) / TMath::Sqrt(changedHist->GetBinContent(1, i+1));
1444
1445       // four values from the derivate (procedure taken from ROOT -- suggestion by Ruben)
1446       // = 1/2D * [ 8 (f(D/2) - f(-D/2)) - (f(D)-f(-D)) ]/3
1447
1448       /*
1449       // check formula
1450       measured->SetBinContent(1, m+1, 100);
1451       newResult[0]->SetBinContent(t+1, 5 * 0.5 + 10);
1452       newResult[1]->SetBinContent(t+1, 5 * -0.5 + 10);
1453       newResult[2]->SetBinContent(t+1, 5 * 1 + 10);
1454       newResult[3]->SetBinContent(t+1, 5 * -1 + 10);
1455       */
1456
1457       matrix(m, t) = 0.5 / TMath::Sqrt(measured->GetBinContent(m+1)) *
1458         ( 8. * (newResult[0]->GetBinContent(t+1) - newResult[1]->GetBinContent(t+1)) -
1459               (newResult[2]->GetBinContent(t+1) - newResult[3]->GetBinContent(t+1))
1460         ) / 3;
1461     }
1462   }
1463
1464   for (Int_t loop=0; loop<4; loop++)
1465     delete newResult[loop];
1466
1467   //matrix.Print();
1468   //TH2* matrixHist = new TH2D(matrix); new TCanvas; matrixHist->Draw("BOX");
1469
1470   // ...calculate folded guess
1471   TH1* convoluted = (TH1*) measured->Clone("convoluted");
1472   convoluted->Reset();
1473   for (Int_t m=0; m<fgMaxInput; m++)
1474   {
1475     Float_t value = 0;
1476     for (Int_t t = 0; t<fgMaxParams; t++)
1477     {
1478       Float_t tmp = correlation->GetBinContent(t+1, m+1) * result->GetBinContent(t+1);
1479       if (efficiency)
1480         tmp *= efficiency->GetBinContent(t+1);
1481       value += tmp;
1482     }
1483     convoluted->SetBinContent(m+1, value);
1484   }
1485
1486   // rescale
1487   convoluted->Scale(measured->Integral() / convoluted->Integral());
1488
1489   //new TCanvas; convoluted->Draw(); measured->Draw("SAME"); measured->SetLineColor(2);
1490   //return;
1491
1492   // difference
1493   convoluted->Add(measured, -1);
1494
1495   // ...and the bias
1496   for (Int_t t = 0; t<fgMaxParams; t++)
1497   {
1498     Double_t error = 0;
1499     for (Int_t m=0; m<fgMaxInput; m++)
1500       error += matrix(m, t) * convoluted->GetBinContent(m+1);
1501     result->SetBinError(t+1, error);
1502   }
1503
1504   //new TCanvas; result->Draw(); gPad->SetLogy();
1505
1506   return 0;
1507 }