ANALYSIS/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 #include <TExec.h>
  31 #include "Riostream.h"
  32 #include "TROOT.h"
  33
  34 using namespace std; //required for resolving the 'cout' symbol
  35
  36 TMatrixD* AliUnfolding::fgCorrelationMatrix = 0;
  37 TMatrixD* AliUnfolding::fgCorrelationMatrixSquared = 0;
  38 TMatrixD* AliUnfolding::fgCorrelationCovarianceMatrix = 0;
  39 TVectorD* AliUnfolding::fgCurrentESDVector = 0;
  40 TVectorD* AliUnfolding::fgEntropyAPriori = 0;
  41 TVectorD* AliUnfolding::fgEfficiency = 0;
  42
  43 TAxis* AliUnfolding::fgUnfoldedAxis = 0;
  44 TAxis* AliUnfolding::fgMeasuredAxis = 0;
  45
  46 TF1* AliUnfolding::fgFitFunction = 0;
  47
  48 AliUnfolding::MethodType AliUnfolding::fgMethodType = AliUnfolding::kInvalid;
  49 Int_t AliUnfolding::fgMaxInput  = -1;  // bins in measured histogram
  50 Int_t AliUnfolding::fgMaxParams = -1;  // bins in unfolded histogram = number of fit params
  51 Float_t AliUnfolding::fgOverflowBinLimit = -1;
  52
  53 AliUnfolding::RegularizationType AliUnfolding::fgRegularizationType = AliUnfolding::kPol1;
  54 Float_t AliUnfolding::fgRegularizationWeight = 10000;
  55 Int_t AliUnfolding::fgSkipBinsBegin = 0;
  56 Float_t AliUnfolding::fgMinuitStepSize = 0.1;                 // (usually not needed to be changed) step size in minimization
  57 Float_t AliUnfolding::fgMinuitPrecision = 1e-6;               // minuit precision
  58 Int_t   AliUnfolding::fgMinuitMaxIterations = 1000000;            // minuit maximum number of iterations
  59 Double_t AliUnfolding::fgMinuitStrategy = 1.;                 // minuit strategy
  60 Bool_t AliUnfolding::fgMinimumInitialValue = kFALSE;          // set all initial values at least to the smallest value among the initial values
  61 Float_t AliUnfolding::fgMinimumInitialValueFix = -1;
  62 Bool_t AliUnfolding::fgNormalizeInput = kFALSE;               // normalize input spectrum
  63 Float_t AliUnfolding::fgNotFoundEvents = 0;
  64 Bool_t AliUnfolding::fgSkipBin0InChi2 = kFALSE;
  65
  66 Float_t AliUnfolding::fgBayesianSmoothing  = 1;           // smoothing parameter (0 = no smoothing)
  67 Int_t   AliUnfolding::fgBayesianIterations = 10;          // number of iterations in Bayesian method
  68
  69 Bool_t AliUnfolding::fgDebug = kFALSE;
  70
  71 Int_t AliUnfolding::fgCallCount = 0;
  72
  73 Int_t AliUnfolding::fgPowern = 5;
  74
  75 Double_t AliUnfolding::fChi2FromFit = 0.;
  76 Double_t AliUnfolding::fPenaltyVal  = 0.;
  77 Double_t AliUnfolding::fAvgResidual = 0.;
  78
  79 Int_t AliUnfolding::fgPrintChi2Details = 0;
  80
  81 TCanvas *AliUnfolding::fgCanvas = 0;
  82 TH1 *AliUnfolding::fghUnfolded = 0;
  83 TH2 *AliUnfolding::fghCorrelation = 0;
  84 TH1 *AliUnfolding::fghEfficiency = 0;
  85 TH1 *AliUnfolding::fghMeasured = 0;
  86
  87 ClassImp(AliUnfolding)
  88
  89 //____________________________________________________________________
  90 void AliUnfolding::SetUnfoldingMethod(MethodType methodType)
  91 {
  92   // set unfolding method
  93   fgMethodType = methodType;
  94
  95   const char* name = 0;
  96   switch (methodType)
  97   {
  98     case kInvalid: name = "INVALID"; break;
  99     case kChi2Minimization: name = "Chi2 Minimization"; break;
 100     case kBayesian: name = "Bayesian unfolding"; break;
 101     case kFunction: name = "Functional fit"; break;
 102   }
 103   Printf("AliUnfolding::SetUnfoldingMethod: %s enabled.", name);
 104 }
 105
 106 //____________________________________________________________________
 107 void AliUnfolding::SetCreateOverflowBin(Float_t overflowBinLimit)
 108 {
 109   // enable the creation of a overflow bin that includes all statistics below the given limit
 110
 111   fgOverflowBinLimit = overflowBinLimit;
 112
 113   Printf("AliUnfolding::SetCreateOverflowBin: overflow bin limit set to %f", overflowBinLimit);
 114 }
 115
 116 //____________________________________________________________________
 117 void AliUnfolding::SetSkipBinsBegin(Int_t nBins)
 118 {
 119   // set number of skipped bins in regularization
 120
 121   fgSkipBinsBegin = nBins;
 122
 123   Printf("AliUnfolding::SetSkipBinsBegin: skipping %d bins at the beginning of the spectrum in the regularization.", fgSkipBinsBegin);
 124 }
 125
 126 //____________________________________________________________________
 127 void AliUnfolding::SetNbins(Int_t nMeasured, Int_t nUnfolded)
 128 {
 129   // set number of bins in the input (measured) distribution and in the unfolded distribution
 130   fgMaxInput = nMeasured;
 131   fgMaxParams = nUnfolded;
 132
 133   if (fgCorrelationMatrix)
 134   {
 135     delete fgCorrelationMatrix;
 136     fgCorrelationMatrix = 0;
 137   }
 138   if (fgCorrelationMatrixSquared)
 139   {
 140     fgCorrelationMatrixSquared = 0;
 141     delete fgCorrelationMatrixSquared;
 142   }
 143   if (fgCorrelationCovarianceMatrix)
 144   {
 145     delete fgCorrelationCovarianceMatrix;
 146     fgCorrelationCovarianceMatrix = 0;
 147   }
 148   if (fgCurrentESDVector)
 149   {
 150     delete fgCurrentESDVector;
 151     fgCurrentESDVector = 0;
 152   }
 153   if (fgEntropyAPriori)
 154   {
 155     delete fgEntropyAPriori;
 156     fgEntropyAPriori = 0;
 157   }
 158   if (fgEfficiency)
 159   {
 160     delete fgEfficiency;
 161     fgEfficiency = 0;
 162   }
 163   if (fgUnfoldedAxis)
 164   {
 165     delete fgUnfoldedAxis;
 166     fgUnfoldedAxis = 0;
 167   }
 168   if (fgMeasuredAxis)
 169   {
 170     delete fgMeasuredAxis;
 171     fgMeasuredAxis = 0;
 172   }
 173
 174   Printf("AliUnfolding::SetNbins: Set %d measured bins and %d unfolded bins", nMeasured, nUnfolded);
 175 }
 176
 177 //____________________________________________________________________
 178 void AliUnfolding::SetChi2Regularization(RegularizationType type, Float_t weight)
 179 {
 180   //
 181   // sets the parameters for chi2 minimization
 182   //
 183
 184   fgRegularizationType = type;
 185   fgRegularizationWeight = weight;
 186
 187   Printf("AliUnfolding::SetChi2Regularization --> Regularization set to %d with weight %f", (Int_t) type, weight);
 188 }
 189
 190 //____________________________________________________________________
 191 void AliUnfolding::SetBayesianParameters(Float_t smoothing, Int_t nIterations)
 192 {
 193   //
 194   // sets the parameters for Bayesian unfolding
 195   //
 196
 197   fgBayesianSmoothing = smoothing;
 198   fgBayesianIterations = nIterations;
 199
 200   Printf("AliUnfolding::SetBayesianParameters --> Parameters set to %d iterations with smoothing %f", fgBayesianIterations, fgBayesianSmoothing);
 201 }
 202
 203 //____________________________________________________________________
 204 void AliUnfolding::SetFunction(TF1* function)
 205 {
 206   // set function for unfolding with a fit function
 207
 208   fgFitFunction = function;
 209
 210   Printf("AliUnfolding::SetFunction: Set fit function with %d parameters.", function->GetNpar());
 211 }
 212
 213 //____________________________________________________________________
 214 Int_t AliUnfolding::Unfold(TH2* correlation, TH1* efficiency, TH1* measured, TH1* initialConditions, TH1* result, Bool_t check)
 215 {
 216   // unfolds with unfolding method fgMethodType
 217   //
 218   // parameters:
 219   //  correlation: response matrix as measured vs. generated
 220   //  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.
 221   //  measured:    the measured spectrum
 222   //  initialConditions: (optional) initial conditions for the unfolding. if 0 the measured spectrum is used as initial conditions.
 223   //  result:      target for the unfolded result
 224   //  check:       depends on the unfolding method, see comments in specific functions
 225   //
 226   //  return code: see UnfoldWithMinuit/UnfoldWithBayesian/UnfoldWithFunction
 227
 228   if (fgMaxInput == -1)
 229   {
 230     Printf("AliUnfolding::Unfold: WARNING. Number of measured bins not set with SetNbins. Using number of bins in measured distribution");
 231     fgMaxInput = measured->GetNbinsX();
 232   }
 233   if (fgMaxParams == -1)
 234   {
 235     Printf("AliUnfolding::Unfold: WARNING. Number of unfolded bins not set with SetNbins. Using number of bins in measured distribution");
 236     fgMaxParams = measured->GetNbinsX();
 237   }
 238
 239   if (fgOverflowBinLimit > 0)
 240     CreateOverflowBin(correlation, measured);
 241
 242   switch (fgMethodType)
 243   {
 244     case kInvalid:
 245     {
 246       Printf("AliUnfolding::Unfold: ERROR: Unfolding method not set. Use SetUnfoldingMethod. Exiting...");
 247       return -1;
 248     }
 249     case kChi2Minimization:
 250       return UnfoldWithMinuit(correlation, efficiency, measured, initialConditions, result, check);
 251     case kBayesian:
 252       return UnfoldWithBayesian(correlation, efficiency, measured, initialConditions, result);
 253     case kFunction:
 254       return UnfoldWithFunction(correlation, efficiency, measured, initialConditions, result);
 255   }
 256
 257
 258
 259   return -1;
 260 }
 261
 262 //____________________________________________________________________
 263 void AliUnfolding::SetStaticVariables(TH2* correlation, TH1* measured, TH1* efficiency)
 264 {
 265   // fill static variables needed for minuit fit
 266
 267   if (!fgCorrelationMatrix)
 268     fgCorrelationMatrix = new TMatrixD(fgMaxInput, fgMaxParams);
 269   if (!fgCorrelationMatrixSquared)
 270     fgCorrelationMatrixSquared = new TMatrixD(fgMaxInput, fgMaxParams);
 271   if (!fgCorrelationCovarianceMatrix)
 272     fgCorrelationCovarianceMatrix = new TMatrixD(fgMaxInput, fgMaxInput);
 273   if (!fgCurrentESDVector)
 274     fgCurrentESDVector = new TVectorD(fgMaxInput);
 275   if (!fgEntropyAPriori)
 276     fgEntropyAPriori = new TVectorD(fgMaxParams);
 277   if (!fgEfficiency)
 278     fgEfficiency = new TVectorD(fgMaxParams);
 279   if (fgUnfoldedAxis)
 280   {
 281     delete fgUnfoldedAxis;
 282     fgUnfoldedAxis = 0;
 283   }
 284   if (!fgUnfoldedAxis)
 285     fgUnfoldedAxis = new TAxis(*(correlation->GetXaxis()));
 286   if (fgMeasuredAxis)
 287   {
 288     delete fgMeasuredAxis;
 289     fgMeasuredAxis = 0;
 290   }
 291   if (!fgMeasuredAxis)
 292     fgMeasuredAxis = new TAxis(*(correlation->GetYaxis()));
 293
 294   fgCorrelationMatrix->Zero();
 295   fgCorrelationCovarianceMatrix->Zero();
 296   fgCurrentESDVector->Zero();
 297   fgEntropyAPriori->Zero();
 298
 299   // normalize correction for given nPart
 300   for (Int_t i=1; i<=correlation->GetNbinsX(); ++i)
 301   {
 302     Double_t sum = correlation->Integral(i, i, 1, correlation->GetNbinsY());
 303     if (sum <= 0)
 304       continue;
 305     Float_t maxValue = 0;
 306     Int_t maxBin = -1;
 307     for (Int_t j=1; j<=correlation->GetNbinsY(); ++j)
 308     {
 309       // find most probably value
 310       if (maxValue < correlation->GetBinContent(i, j))
 311       {
 312         maxValue = correlation->GetBinContent(i, j);
 313         maxBin = j;
 314       }
 315
 316       // npart sum to 1
 317       correlation->SetBinContent(i, j, correlation->GetBinContent(i, j) / sum);// * correlation->GetXaxis()->GetBinWidth(i));
 318       correlation->SetBinError(i, j, correlation->GetBinError(i, j) / sum);
 319
 320       if (i <= fgMaxParams && j <= fgMaxInput)
 321       {
 322         (*fgCorrelationMatrix)(j-1, i-1) = correlation->GetBinContent(i, j);
 323         (*fgCorrelationMatrixSquared)(j-1, i-1) = correlation->GetBinContent(i, j) * correlation->GetBinContent(i, j);
 324       }
 325     }
 326
 327     //printf("MPV for Ntrue = %f is %f\n", fCurrentCorrelation->GetXaxis()->GetBinCenter(i), fCurrentCorrelation->GetYaxis()->GetBinCenter(maxBin));
 328   }
 329
 330   //normalize measured
 331   Float_t smallestError = 1;
 332   if (fgNormalizeInput)
 333   {
 334     Float_t sumMeasured = measured->Integral();
 335     measured->Scale(1.0 / sumMeasured);
 336     smallestError /= sumMeasured;
 337   }
 338
 339   for (Int_t i=0; i<fgMaxInput; ++i)
 340   {
 341     (*fgCurrentESDVector)[i] = measured->GetBinContent(i+1);
 342     if (measured->GetBinError(i+1) > 0)
 343     {
 344       (*fgCorrelationCovarianceMatrix)(i, i) = (Double_t) 1e-6 / measured->GetBinError(i+1) / measured->GetBinError(i+1);
 345     }
 346     else // in this case put error of 1, otherwise 0 bins are not added to the chi2...
 347       (*fgCorrelationCovarianceMatrix)(i, i) = (Double_t) 1e-6 / smallestError / smallestError;
 348
 349     if ((*fgCorrelationCovarianceMatrix)(i, i) > 1e7)
 350       (*fgCorrelationCovarianceMatrix)(i, i) = 0;
 351     //Printf("%d, %e", i, (*fgCorrelationCovarianceMatrix)(i, i));
 352   }
 353
 354   // efficiency is expected to match bin width of result
 355   for (Int_t i=0; i<fgMaxParams; ++i)
 356   {
 357     (*fgEfficiency)(i) = efficiency->GetBinContent(i+1);
 358   }
 359
 360   if (correlation->GetNbinsX() != fgMaxParams || correlation->GetNbinsY() != fgMaxInput)
 361     cout << "Response histo has incorrect dimensions; expect (" << fgMaxParams << ", " << fgMaxInput << "), got (" << correlation->GetNbinsX() << ", " << correlation->GetNbinsY() << ")" << endl;
 362
 363 }
 364
 365 //____________________________________________________________________
 366 Int_t AliUnfolding::UnfoldWithMinuit(TH2* correlation, TH1* efficiency, TH1* measured, TH1* initialConditions, TH1* result, Bool_t check)
 367 {
 368   //
 369   // implementation of unfolding (internal function)
 370   //
 371   // unfolds <measured> using response from <correlation> and effiency <efficiency>
 372   // output is in <result>
 373   // <initialConditions> set the initial values for the minimization, if 0 <measured> is used
 374   //   negative values in initialConditions mean that the given parameter is fixed to the absolute of the value
 375   // if <check> is true no unfolding is made, instead only the chi2 without unfolding is printed
 376   //
 377   // returns minuit status (0 = success), (-1 when check was set)
 378   //
 379
 380   SetStaticVariables(correlation, measured, efficiency);
 381
 382   // Initialize TMinuit via generic fitter interface
 383   Int_t params = fgMaxParams;
 384   if (fgNotFoundEvents > 0)
 385     params++;
 386
 387   TVirtualFitter *minuit = TVirtualFitter::Fitter(0, params);
 388   Double_t arglist[100];
 389   //  minuit->SetDefaultFitter("Minuit2");
 390
 391   // disable any output (-1), unfortuantly we do not see warnings anymore then. Have to find another way...
 392   arglist[0] = 0;
 393   minuit->ExecuteCommand("SET PRINT", arglist, 1);
 394
 395   // however, enable warnings
 396   //minuit->ExecuteCommand("SET WAR", arglist, 0);
 397
 398   // set minimization function
 399   minuit->SetFCN(Chi2Function);
 400
 401   // set precision
 402   minuit->SetPrecision(fgMinuitPrecision);
 403
 404   minuit->SetMaxIterations(fgMinuitMaxIterations);
 405
 406   minuit->ExecuteCommand("SET STRATEGY",&fgMinuitStrategy,1);
 407
 408   for (Int_t i=0; i<fgMaxParams; i++)
 409     (*fgEntropyAPriori)[i] = 1;
 410
 411   // set initial conditions as a-priori distribution for MRX regularization
 412   /*
 413   for (Int_t i=0; i<fgMaxParams; i++)
 414     if (initialConditions && initialConditions->GetBinContent(i+1) > 0)
 415       (*fgEntropyAPriori)[i] = initialConditions->GetBinContent(i+1);
 416   */
 417
 418   if (!initialConditions) {
 419     initialConditions = measured;
 420   } else {
 421     Printf("AliUnfolding::UnfoldWithMinuit: Using different initial conditions...");
 422     //new TCanvas; initialConditions->DrawCopy();
 423     if (fgNormalizeInput)
 424       initialConditions->Scale(1.0 / initialConditions->Integral());
 425   }
 426
 427   // extract minimum value from initial conditions (if we set a value to 0 it will stay 0)
 428   Float_t minValue = 1e35;
 429   if (fgMinimumInitialValueFix < 0)
 430   {
 431     for (Int_t i=0; i<fgMaxParams; ++i)
 432     {
 433       Int_t bin = initialConditions->GetXaxis()->FindBin(result->GetXaxis()->GetBinCenter(i+1));
 434       if (initialConditions->GetBinContent(bin) > 0)
 435         minValue = TMath::Min(minValue, (Float_t) initialConditions->GetBinContent(bin));
 436     }
 437   }
 438   else
 439     minValue = fgMinimumInitialValueFix;
 440
 441   Double_t* results = new Double_t[fgMaxParams+1];
 442   for (Int_t i=0; i<fgMaxParams; ++i)
 443   {
 444     Int_t bin = initialConditions->GetXaxis()->FindBin(result->GetXaxis()->GetBinCenter(i+1));
 445     results[i] = initialConditions->GetBinContent(bin);
 446
 447     Bool_t fix = kFALSE;
 448     if (results[i] < 0)
 449     {
 450       fix = kTRUE;
 451       results[i] = -results[i];
 452     }
 453
 454     if (!fix && fgMinimumInitialValue && results[i] < minValue)
 455       results[i] = minValue;
 456
 457     // minuit sees squared values to prevent it from going negative...
 458     results[i] = TMath::Sqrt(results[i]);
 459
 460     minuit->SetParameter(i, Form("param%d", i), results[i], (fix) ? 0 : fgMinuitStepSize, 0, 0);
 461   }
 462   if (fgNotFoundEvents > 0)
 463   {
 464     results[fgMaxParams] = efficiency->GetBinContent(1);
 465     minuit->SetParameter(fgMaxParams, "vtx0", results[fgMaxParams], fgMinuitStepSize / 100, 0.01, 0.80);
 466   }
 467
 468   Int_t dummy = 0;
 469   Double_t chi2 = 0;
 470   Chi2Function(dummy, 0, chi2, results, 0);
 471   printf("AliUnfolding::UnfoldWithMinuit: Chi2 of initial parameters is = %f\n", chi2);
 472
 473   if (check)
 474   {
 475     DrawGuess(results);
 476     delete[] results;
 477     return -1;
 478   }
 479
 480   // first param is number of iterations, second is precision....
 481   arglist[0] = (float)fgMinuitMaxIterations;
 482   //  arglist[1] = 1e-5;
 483   //  minuit->ExecuteCommand("SET PRINT", arglist, 3);
 484   //  minuit->ExecuteCommand("SCAN", arglist, 0);
 485   Int_t status = minuit->ExecuteCommand("MIGRAD", arglist, 1);
 486   Printf("AliUnfolding::UnfoldWithMinuit: MINUIT status is %d", status);
 487   //printf("!!!!!!!!!!!!!! MIGRAD finished: Starting MINOS !!!!!!!!!!!!!!");
 488   //minuit->ExecuteCommand("MINOS", arglist, 0);
 489
 490   if (fgNotFoundEvents > 0)
 491   {
 492     results[fgMaxParams] = minuit->GetParameter(fgMaxParams);
 493     Printf("Efficiency for bin 0 changed from %f to %f", efficiency->GetBinContent(1), results[fgMaxParams]);
 494     efficiency->SetBinContent(1, results[fgMaxParams]);
 495   }
 496
 497   for (Int_t i=0; i<fgMaxParams; ++i)
 498   {
 499     results[i] = minuit->GetParameter(i);
 500     Double_t value = results[i] * results[i];
 501    // error is : 2 * (relError on results[i]) * (value) = 2 * (minuit->GetParError(i) / minuit->GetParameter(i)) * (minuit->GetParameter(i) * minuit->GetParameter(i))
 502     Double_t error = 0;
 503     if (TMath::IsNaN(minuit->GetParError(i)))
 504       Printf("WARNING: Parameter %d error is nan", i);
 505     else
 506       error = 2 * minuit->GetParError(i) * results[i];
 507
 508     if (efficiency)
 509     {
 510       //printf("value before efficiency correction: %f\n",value);
 511       if (efficiency->GetBinContent(i+1) > 0)
 512       {
 513         value /= efficiency->GetBinContent(i+1);
 514         error /= efficiency->GetBinContent(i+1);
 515       }
 516       else
 517       {
 518         value = 0;
 519         error = 0;
 520       }
 521     }
 522     //printf("value after efficiency correction: %f +/- %f\n",value,error);
 523     result->SetBinContent(i+1, value);
 524     result->SetBinError(i+1, error);
 525   }
 526
 527   Int_t tmpCallCount = fgCallCount;
 528   fgCallCount = 0; // needs to be 0 so that the Chi2Function prints its output
 529   Chi2Function(dummy, 0, chi2, results, 0);
 530
 531   Printf("AliUnfolding::UnfoldWithMinuit: iterations %d. Chi2 of final parameters is = %f", tmpCallCount, chi2);
 532
 533   delete[] results;
 534
 535   return status;
 536 }
 537
 538 //____________________________________________________________________
 539 Int_t AliUnfolding::UnfoldWithBayesian(TH2* correlation, TH1* aEfficiency, TH1* measured, TH1* initialConditions, TH1* aResult)
 540 {
 541   //
 542   // unfolds a spectrum using the Bayesian method
 543   //
 544
 545   if (measured->Integral() <= 0)
 546   {
 547     Printf("AliUnfolding::UnfoldWithBayesian: ERROR: The measured spectrum is empty");
 548     return -1;
 549   }
 550
 551   const Int_t kStartBin = 0;
 552
 553   Int_t kMaxM = fgMaxInput;  //<= fCurrentCorrelation->GetNbinsY(); // max measured axis
 554   Int_t kMaxT = fgMaxParams; //<= fCurrentCorrelation->GetNbinsX(); // max true axis
 555
 556   // convergence limit: kMaxT * 0.001^2 = kMaxT * 1e-6 (e.g. 250 bins --> 2.5 e-4)
 557   const Double_t kConvergenceLimit = kMaxT * 1e-6;
 558
 559   // store information in arrays, to increase processing speed (~ factor 5)
 560   Double_t* measuredCopy = new Double_t[kMaxM];
 561   Double_t* measuredError = new Double_t[kMaxM];
 562   Double_t* prior = new Double_t[kMaxT];
 563   Double_t* result = new Double_t[kMaxT];
 564   Double_t* efficiency = new Double_t[kMaxT];
 565   Double_t* binWidths = new Double_t[kMaxT];
 566   Double_t* normResponse = new Double_t[kMaxT];
 567
 568   Double_t** response = new Double_t*[kMaxT];
 569   Double_t** inverseResponse = new Double_t*[kMaxT];
 570   for (Int_t i=0; i<kMaxT; i++)
 571   {
 572     response[i] = new Double_t[kMaxM];
 573     inverseResponse[i] = new Double_t[kMaxM];
 574     normResponse[i] = 0;
 575   }
 576
 577   // for normalization
 578   Float_t measuredIntegral = measured->Integral();
 579   for (Int_t m=0; m<kMaxM; m++)
 580   {
 581     measuredCopy[m] = measured->GetBinContent(m+1) / measuredIntegral;
 582     measuredError[m] = measured->GetBinError(m+1) / measuredIntegral;
 583
 584     for (Int_t t=0; t<kMaxT; t++)
 585     {
 586       response[t][m] = correlation->GetBinContent(t+1, m+1);
 587       inverseResponse[t][m] = 0;
 588       normResponse[t] += correlation->GetBinContent(t+1, m+1);
 589     }
 590   }
 591
 592   for (Int_t t=0; t<kMaxT; t++)
 593   {
 594     if (aEfficiency)
 595     {
 596       efficiency[t] = aEfficiency->GetBinContent(t+1);
 597     }
 598     else
 599       efficiency[t] = 1;
 600
 601     prior[t] = measuredCopy[t];
 602     result[t] = 0;
 603     binWidths[t] = aResult->GetXaxis()->GetBinWidth(t+1);
 604
 605     for (Int_t m=0; m<kMaxM; m++) { // Normalise response matrix
 606       if (normResponse[t] != 0)
 607         response[t][m] /= normResponse[t];
 608       else
 609         Printf("AliUnfolding::UnfoldWithBayesian: Empty row,column in response matrix, for truth bin %d",t);
 610     }
 611   }
 612
 613   // pick prior distribution
 614   if (initialConditions)
 615   {
 616     printf("Using different starting conditions...\n");
 617     // for normalization
 618     Float_t inputDistIntegral = initialConditions->Integral();
 619     for (Int_t i=0; i<kMaxT; i++)
 620       prior[i] = initialConditions->GetBinContent(i+1) / inputDistIntegral;
 621   }
 622
 623   //TH1F* convergence = new TH1F("convergence", "convergence", 200, 0.5, 200.5);
 624
 625   //new TCanvas;
 626   // unfold...
 627   for (Int_t i=0; i<fgBayesianIterations || fgBayesianIterations < 0; i++)
 628   {
 629     if (fgDebug)
 630       Printf("AliUnfolding::UnfoldWithBayesian: iteration %i", i);
 631
 632     // calculate IR from Bayes theorem
 633     // IR_ji = R_ij * prior_i / sum_k(R_kj * prior_k)
 634
 635     Double_t chi2Measured = 0;
 636     for (Int_t m=0; m<kMaxM; m++)
 637     {
 638       Float_t norm = 0;
 639       for (Int_t t = kStartBin; t<kMaxT; t++)
 640         norm += response[t][m] * prior[t] * efficiency[t];
 641
 642       // calc. chi2: (measured - response * prior) / error
 643       if (measuredError[m] > 0)
 644       {
 645         Double_t value = (measuredCopy[m] - norm) / measuredError[m];
 646         chi2Measured += value * value;
 647       }
 648
 649       if (norm > 0)
 650       {
 651         for (Int_t t = kStartBin; t<kMaxT; t++)
 652           inverseResponse[t][m] = response[t][m] * prior[t] * efficiency[t] / norm;
 653       }
 654       else
 655       {
 656         for (Int_t t = kStartBin; t<kMaxT; t++)
 657           inverseResponse[t][m] = 0;
 658       }
 659     }
 660     //Printf("chi2Measured of the last prior is %e", chi2Measured);
 661
 662     for (Int_t t = kStartBin; t<kMaxT; t++)
 663     {
 664       Float_t value = 0;
 665       for (Int_t m=0; m<kMaxM; m++)
 666         value += inverseResponse[t][m] * measuredCopy[m];
 667
 668       if (efficiency[t] > 0)
 669         result[t] = value / efficiency[t];
 670       else
 671         result[t] = 0;
 672     }
 673
 674     /*
 675     // draw intermediate result
 676     for (Int_t t=0; t<kMaxT; t++)
 677     {
 678       aResult->SetBinContent(t+1, result[t]);
 679     }
 680     aResult->SetMarkerStyle(24+i);
 681     aResult->SetMarkerColor(2);
 682     aResult->DrawCopy((i == 0) ? "P" : "PSAME");
 683     */
 684
 685     Double_t chi2LastIter = 0;
 686     // regularization (simple smoothing)
 687     for (Int_t t=kStartBin; t<kMaxT; t++)
 688     {
 689       Float_t newValue = 0;
 690
 691       // 0 bin excluded from smoothing
 692       if (t > kStartBin+2 && t<kMaxT-1)
 693       {
 694         Float_t average = (result[t-1] / binWidths[t-1] + result[t] / binWidths[t] + result[t+1] / binWidths[t+1]) / 3 * binWidths[t];
 695
 696         // weight the average with the regularization parameter
 697         newValue = (1 - fgBayesianSmoothing) * result[t] + fgBayesianSmoothing * average;
 698       }
 699       else
 700         newValue = result[t];
 701
 702       // calculate chi2 (change from last iteration)
 703       if (prior[t] > 1e-5)
 704       {
 705         Double_t diff = (prior[t] - newValue) / prior[t];
 706         chi2LastIter += diff * diff;
 707       }
 708
 709       prior[t] = newValue;
 710     }
 711     //printf("Chi2 of %d iteration = %e\n", i, chi2LastIter);
 712     //convergence->Fill(i+1, chi2LastIter);
 713
 714     if (fgBayesianIterations < 0 && chi2LastIter < kConvergenceLimit)
 715     {
 716       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);
 717       break;
 718     }
 719   } // end of iterations
 720
 721   //new TCanvas; convergence->DrawCopy(); gPad->SetLogy();
 722   //delete convergence;
 723
 724   Float_t factor = 1;
 725   if (!fgNormalizeInput)
 726     factor = measuredIntegral;
 727   for (Int_t t=0; t<kMaxT; t++)
 728     aResult->SetBinContent(t+1, result[t] * factor);
 729
 730   delete[] measuredCopy;
 731   delete[] measuredError;
 732   delete[] prior;
 733   delete[] result;
 734   delete[] efficiency;
 735   delete[] binWidths;
 736   delete[] normResponse;
 737
 738   for (Int_t i=0; i<kMaxT; i++)
 739   {
 740     delete[] response[i];
 741     delete[] inverseResponse[i];
 742   }
 743   delete[] response;
 744   delete[] inverseResponse;
 745
 746   return 0;
 747
 748   // ********
 749   // Calculate the covariance matrix, all arguments are taken from NIM,A362,487-498,1995
 750
 751   /*printf("Calculating covariance matrix. This may take some time...\n");
 752
 753   // check if this is the right one...
 754   TH1* sumHist = GetMultiplicityMC(inputRange, eventType)->ProjectionY("sumHist", 1, GetMultiplicityMC(inputRange, eventType)->GetNbinsX());
 755
 756   Int_t xBins = hInverseResponseBayes->GetNbinsX();
 757   Int_t yBins = hInverseResponseBayes->GetNbinsY();
 758
 759   // calculate "unfolding matrix" Mij
 760   Float_t matrixM[251][251];
 761   for (Int_t i=1; i<=xBins; i++)
 762   {
 763     for (Int_t j=1; j<=yBins; j++)
 764     {
 765       if (fCurrentEfficiency->GetBinContent(i) > 0)
 766         matrixM[i-1][j-1] = hInverseResponseBayes->GetBinContent(i, j) / fCurrentEfficiency->GetBinContent(i);
 767       else
 768         matrixM[i-1][j-1] = 0;
 769     }
 770   }
 771
 772   Float_t* vectorn = new Float_t[yBins];
 773   for (Int_t j=1; j<=yBins; j++)
 774     vectorn[j-1] = fCurrentESD->GetBinContent(j);
 775
 776   // first part of covariance matrix, depends on input distribution n(E)
 777   Float_t cov1[251][251];
 778
 779   Float_t nEvents = fCurrentESD->Integral(); // N
 780
 781   xBins = 20;
 782   yBins = 20;
 783
 784   for (Int_t k=0; k<xBins; k++)
 785   {
 786     printf("In Cov1: %d\n", k);
 787     for (Int_t l=0; l<yBins; l++)
 788     {
 789       cov1[k][l] = 0;
 790
 791       // sum_j Mkj Mlj n(Ej) * (1 - n(Ej) / N)
 792       for (Int_t j=0; j<yBins; j++)
 793         cov1[k][l] += matrixM[k][j] * matrixM[l][j] * vectorn[j]
 794           * (1.0 - vectorn[j] / nEvents);
 795
 796       // - sum_i,j (i != j) Mki Mlj n(Ei) n(Ej) / N
 797       for (Int_t i=0; i<yBins; i++)
 798         for (Int_t j=0; j<yBins; j++)
 799         {
 800           if (i == j)
 801             continue;
 802           cov1[k][l] -= matrixM[k][i] * matrixM[l][j] * vectorn[i]
 803             * vectorn[j] / nEvents;
 804          }
 805     }
 806   }
 807
 808   printf("Cov1 finished\n");
 809
 810   TH2F* cov = (TH2F*) hInverseResponseBayes->Clone("cov");
 811   cov->Reset();
 812
 813   for (Int_t i=1; i<=xBins; i++)
 814     for (Int_t j=1; j<=yBins; j++)
 815       cov->SetBinContent(i, j, cov1[i-1][j-1]);
 816
 817   new TCanvas;
 818   cov->Draw("COLZ");
 819
 820   // second part of covariance matrix, depends on response matrix
 821   Float_t cov2[251][251];
 822
 823   // Cov[P(Er|Cu), P(Es|Cu)] term
 824   Float_t covTerm[100][100][100];
 825   for (Int_t r=0; r<yBins; r++)
 826     for (Int_t u=0; u<xBins; u++)
 827       for (Int_t s=0; s<yBins; s++)
 828       {
 829         if (r == s)
 830           covTerm[r][u][s] = 1.0 / sumHist->GetBinContent(u+1) * hResponse->GetBinContent(u+1, r+1)
 831             * (1.0 - hResponse->GetBinContent(u+1, r+1));
 832         else
 833           covTerm[r][u][s] = - 1.0 / sumHist->GetBinContent(u+1) * hResponse->GetBinContent(u+1, r+1)
 834             * hResponse->GetBinContent(u+1, s+1);
 835       }
 836
 837   for (Int_t k=0; k<xBins; k++)
 838     for (Int_t l=0; l<yBins; l++)
 839     {
 840       cov2[k][l] = 0;
 841       printf("In Cov2: %d %d\n", k, l);
 842       for (Int_t i=0; i<yBins; i++)
 843         for (Int_t j=0; j<yBins; j++)
 844         {
 845           //printf("In Cov2: %d %d %d %d\n", k, l, i, j);
 846           // calculate Cov(Mki, Mlj) = sum{ru},{su} ...
 847           Float_t tmpCov = 0;
 848           for (Int_t r=0; r<yBins; r++)
 849             for (Int_t u=0; u<xBins; u++)
 850               for (Int_t s=0; s<yBins; s++)
 851               {
 852                 if (hResponse->GetBinContent(u+1, r+1) == 0 || hResponse->GetBinContent(u+1, s+1) == 0
 853                   || hResponse->GetBinContent(u+1, i+1) == 0)
 854                   continue;
 855
 856                 tmpCov += BayesCovarianceDerivate(matrixM, hResponse, fCurrentEfficiency, k, i, r, u)
 857                   * BayesCovarianceDerivate(matrixM, hResponse, fCurrentEfficiency, l, j, s, u)
 858                   * covTerm[r][u][s];
 859               }
 860
 861           cov2[k][l] += fCurrentESD->GetBinContent(i+1) * fCurrentESD->GetBinContent(j+1) * tmpCov;
 862         }
 863     }
 864
 865   printf("Cov2 finished\n");
 866
 867   for (Int_t i=1; i<=xBins; i++)
 868     for (Int_t j=1; j<=yBins; j++)
 869       cov->SetBinContent(i, j, cov1[i-1][j-1] + cov2[i-1][j-1]);
 870
 871   new TCanvas;
 872   cov->Draw("COLZ");*/
 873 }
 874
 875 //____________________________________________________________________
 876 Double_t AliUnfolding::RegularizationPol0(TVectorD& params)
 877 {
 878   // homogenity term for minuit fitting
 879   // pure function of the parameters
 880   // prefers constant function (pol0)
 881   //
 882   // Does not take into account efficiency
 883   Double_t chi2 = 0;
 884
 885   for (Int_t i=1+fgSkipBinsBegin; i<fgMaxParams; ++i)
 886   {
 887     Double_t right  = params[i] / fgUnfoldedAxis->GetBinWidth(i+1);
 888     Double_t left   = params[i-1] / fgUnfoldedAxis->GetBinWidth(i);
 889
 890     if (left != 0)
 891     {
 892       Double_t diff = (right - left);
 893       chi2 += diff * diff / left / ((fgUnfoldedAxis->GetBinWidth(i+1) + fgUnfoldedAxis->GetBinWidth(i)) / 2);
 894     }
 895   }
 896
 897   return chi2 / 100.0;
 898 }
 899
 900 //____________________________________________________________________
 901 Double_t AliUnfolding::RegularizationPol1(TVectorD& params)
 902 {
 903   // homogenity term for minuit fitting
 904   // pure function of the parameters
 905   // prefers linear function (pol1)
 906   //
 907   // Does not take into account efficiency
 908   Double_t chi2 = 0;
 909
 910   for (Int_t i=2+fgSkipBinsBegin; i<fgMaxParams; ++i)
 911   {
 912     if (params[i-1] == 0)
 913       continue;
 914
 915     Double_t right  = params[i] / fgUnfoldedAxis->GetBinWidth(i+1);
 916     Double_t middle = params[i-1] / fgUnfoldedAxis->GetBinWidth(i);
 917     Double_t left   = params[i-2] / fgUnfoldedAxis->GetBinWidth(i-1);
 918
 919     Double_t der1 = (right - middle) / ((fgUnfoldedAxis->GetBinWidth(i+1) + fgUnfoldedAxis->GetBinWidth(i)) / 2);
 920     Double_t der2 = (middle - left) / ((fgUnfoldedAxis->GetBinWidth(i) + fgUnfoldedAxis->GetBinWidth(i-1)) / 2);
 921
 922     //Double_t diff = (der1 - der2) / middle;
 923     //chi2 += diff * diff;
 924     chi2 += (der1 - der2) * (der1 - der2) / middle * fgUnfoldedAxis->GetBinWidth(i);
 925   }
 926
 927   return chi2;
 928 }
 929
 930 //____________________________________________________________________
 931 Double_t AliUnfolding::RegularizationLog(TVectorD& params)
 932 {
 933   // homogenity term for minuit fitting
 934   // pure function of the parameters
 935   // prefers logarithmic function (log)
 936   //
 937   // Does not take into account efficiency
 938
 939   Double_t chi2 = 0;
 940
 941   for (Int_t i=2+fgSkipBinsBegin; i<fgMaxParams; ++i)
 942   {
 943     if (params[i-1] == 0 || params[i] == 0 || params[i-2] == 0)
 944      continue;
 945
 946     Double_t right  = log(params[i] / fgUnfoldedAxis->GetBinWidth(i+1));
 947     Double_t middle = log(params[i-1] / fgUnfoldedAxis->GetBinWidth(i));
 948     Double_t left   = log(params[i-2] / fgUnfoldedAxis->GetBinWidth(i-1));
 949
 950     Double_t der1 = (right - middle) / ((fgUnfoldedAxis->GetBinWidth(i+1) + fgUnfoldedAxis->GetBinWidth(i)) / 2);
 951     Double_t der2 = (middle - left) / ((fgUnfoldedAxis->GetBinWidth(i) + fgUnfoldedAxis->GetBinWidth(i-1)) / 2);
 952
 953     //Double_t error = 1. / params[i] + 4. / params[i-1] + 1. / params[i-2];
 954
 955     //if (fgCallCount == 0)
 956     //  Printf("%d %f %f", i, (der1 - der2) * (der1 - der2), error);
 957     chi2 += (der1 - der2) * (der1 - der2);// / error;
 958   }
 959
 960   return chi2;
 961 }
 962
 963 //____________________________________________________________________
 964 Double_t AliUnfolding::RegularizationTotalCurvature(TVectorD& params)
 965 {
 966   // homogenity term for minuit fitting
 967   // pure function of the parameters
 968   // minimizes the total curvature (from Unfolding Methods In High-Energy Physics Experiments,
 969   // V. Blobel (Hamburg U.) . DESY 84/118, Dec 1984. 40pp.
 970   //
 971   // Does not take into account efficiency
 972
 973   Double_t chi2 = 0;
 974
 975   for (Int_t i=2+fgSkipBinsBegin; i<fgMaxParams; ++i)
 976   {
 977     Double_t right  = params[i];
 978     Double_t middle = params[i-1];
 979     Double_t left   = params[i-2];
 980
 981     Double_t der1 = (right - middle);
 982     Double_t der2 = (middle - left);
 983
 984     Double_t diff = (der1 - der2);
 985
 986     chi2 += diff * diff;
 987   }
 988
 989   return chi2 * 1e4;
 990 }
 991
 992 //____________________________________________________________________
 993 Double_t AliUnfolding::RegularizationEntropy(TVectorD& params)
 994 {
 995   // homogenity term for minuit fitting
 996   // pure function of the parameters
 997   // calculates entropy, from
 998   // The method of reduced cross-entropy (M. Schmelling 1993)
 999   //
1000   // Does not take into account efficiency
1001
1002   Double_t paramSum = 0;
1003
1004   for (Int_t i=fgSkipBinsBegin; i<fgMaxParams; ++i)
1005     paramSum += params[i];
1006
1007   Double_t chi2 = 0;
1008   for (Int_t i=fgSkipBinsBegin; i<fgMaxParams; ++i)
1009   {
1010     Double_t tmp = params[i] / paramSum;
1011     //Double_t tmp = params[i];
1012     if (tmp > 0 && (*fgEntropyAPriori)[i] > 0)
1013     {
1014       chi2 += tmp * TMath::Log(tmp / (*fgEntropyAPriori)[i]);
1015     }
1016     else
1017       chi2 += 100;
1018   }
1019
1020   return -chi2;
1021 }
1022
1023 //____________________________________________________________________
1024 Double_t AliUnfolding::RegularizationRatio(TVectorD& params)
1025 {
1026   // homogenity term for minuit fitting
1027   // pure function of the parameters
1028   //
1029   // Does not take into account efficiency
1030
1031   Double_t chi2 = 0;
1032
1033   for (Int_t i=5+fgSkipBinsBegin; i<fgMaxParams; ++i)
1034   {
1035     if (params[i-1] == 0 || params[i] == 0)
1036       continue;
1037
1038     Double_t right  = params[i] / fgUnfoldedAxis->GetBinWidth(i+1);
1039     Double_t middle = params[i-1] / fgUnfoldedAxis->GetBinWidth(i);
1040     Double_t left   = params[i-2] / fgUnfoldedAxis->GetBinWidth(i-1);
1041     Double_t left2   = params[i-3] / fgUnfoldedAxis->GetBinWidth(i-2);
1042     Double_t left3   = params[i-4] / fgUnfoldedAxis->GetBinWidth(i-3);
1043     Double_t left4   = params[i-5] / fgUnfoldedAxis->GetBinWidth(i-4);
1044
1045     //Double_t diff = left / middle - middle / right;
1046     //Double_t diff = 2 * left / middle - middle / right - left2 / left;
1047     Double_t diff = 4 * left2 / left - middle / right - left / middle - left3 / left2 - left4 / left3;
1048
1049     chi2 += diff * diff;// / middle;
1050   }
1051
1052   return chi2;
1053 }
1054
1055 //____________________________________________________________________
1056 Double_t AliUnfolding::RegularizationPowerLaw(TVectorD& params)
1057 {
1058   // homogenity term for minuit fitting
1059   // pure function of the parameters
1060   // prefers power law with n = -5
1061   //
1062   // Does not take into account efficiency
1063
1064   Double_t chi2 = 0;
1065
1066   Double_t right = 0.;
1067   Double_t middle = 0.;
1068   Double_t left = 0.;
1069
1070   for (Int_t i=2+fgSkipBinsBegin; i<fgMaxParams; ++i)
1071   {
1072     if (params[i] < 1e-8 || params[i-1] < 1e-8 || params[i-2] < 1e-8)
1073       continue;
1074
1075     if (fgUnfoldedAxis->GetBinWidth(i+1) < 1e-8 || fgUnfoldedAxis->GetBinWidth(i) < 1e-8 || fgUnfoldedAxis->GetBinWidth(i-1) < 1e-8)
1076       continue;
1077
1078     middle = TMath::Power(params[i-1] / fgUnfoldedAxis->GetBinWidth(i),fgPowern);
1079
1080     if(middle>0) {
1081       right  = TMath::Power(params[i] / fgUnfoldedAxis->GetBinWidth(i),fgPowern)/middle;
1082
1083       left   = TMath::Power(params[i-2] / fgUnfoldedAxis->GetBinWidth(i-1),fgPowern)/middle;
1084
1085       middle = 1.;
1086
1087       Double_t der1 = (right - middle) / ((fgUnfoldedAxis->GetBinWidth(i+1) + fgUnfoldedAxis->GetBinWidth(i)) / 2);
1088       Double_t der2 = (middle - left) / ((fgUnfoldedAxis->GetBinWidth(i-1) + fgUnfoldedAxis->GetBinWidth(i-2)) / 2);
1089
1090       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;
1091       //   printf("i: %d chi2 = %f\n",i,chi2);
1092     }
1093
1094   }
1095
1096   return chi2;
1097 }
1098
1099 //____________________________________________________________________
1100 Double_t AliUnfolding::RegularizationLogLog(TVectorD& params)
1101 {
1102   // homogenity term for minuit fitting
1103   // pure function of the parameters
1104   // prefers a powerlaw (linear on a log-log scale)
1105   //
1106   // The calculation takes into account the efficiencies
1107
1108   Double_t chi2 = 0;
1109
1110   for (Int_t i=2+fgSkipBinsBegin; i<fgMaxParams; ++i)
1111   {
1112     if (params[i-1] == 0 || params[i] == 0 || params[i-2] == 0)
1113      continue;
1114     if ((*fgEfficiency)(i-1) == 0 || (*fgEfficiency)(i) == 0 || (*fgEfficiency)(i-2) == 0)
1115      continue;
1116
1117
1118     Double_t right  = log(params[i] / (*fgEfficiency)(i) / fgUnfoldedAxis->GetBinWidth(i));
1119     Double_t middle = log(params[i-1] / (*fgEfficiency)(i-1) / fgUnfoldedAxis->GetBinWidth(i-1));
1120     Double_t left   = log(params[i-2] / (*fgEfficiency)(i-2) / fgUnfoldedAxis->GetBinWidth(i-2));
1121
1122     Double_t der1 = (right - middle) / ( log(fgUnfoldedAxis->GetBinCenter(i+1)) - log(fgUnfoldedAxis->GetBinCenter(i)) );
1123     Double_t der2 = (middle - left) /( log(fgUnfoldedAxis->GetBinCenter(i)) - log(fgUnfoldedAxis->GetBinCenter(i-1)) );
1124
1125     double tmp = (log(fgUnfoldedAxis->GetBinCenter(i+1)) - log(fgUnfoldedAxis->GetBinCenter(i-1)))/2.;
1126     Double_t dder = (der1-der2) / tmp;
1127
1128     chi2 += dder * dder;
1129   }
1130
1131   return chi2;
1132 }
1133
1134
1135
1136 //____________________________________________________________________
1137 void AliUnfolding::Chi2Function(Int_t&, Double_t*, Double_t& chi2, Double_t *params, Int_t)
1138 {
1139   //
1140   // fit function for minuit
1141   // does: (m - Ad)W(m - Ad) where m = measured, A correlation matrix, d = guess, W = covariance matrix
1142   //
1143
1144   // TODO use static members for the variables here to speed up processing (no construction/deconstruction)
1145
1146   // d = guess
1147   TVectorD paramsVector(fgMaxParams);
1148   for (Int_t i=0; i<fgMaxParams; ++i)
1149     paramsVector[i] = params[i] * params[i];
1150
1151   // calculate penalty factor
1152   Double_t penaltyVal = 0;
1153
1154   switch (fgRegularizationType)
1155   {
1156     case kNone:       break;
1157     case kPol0:       penaltyVal = RegularizationPol0(paramsVector); break;
1158     case kPol1:       penaltyVal = RegularizationPol1(paramsVector); break;
1159     case kCurvature:  penaltyVal = RegularizationTotalCurvature(paramsVector); break;
1160     case kEntropy:    penaltyVal = RegularizationEntropy(paramsVector); break;
1161     case kLog:        penaltyVal = RegularizationLog(paramsVector); break;
1162     case kRatio:      penaltyVal = RegularizationRatio(paramsVector); break;
1163     case kPowerLaw:   penaltyVal = RegularizationPowerLaw(paramsVector); break;
1164     case kLogLog:     penaltyVal = RegularizationLogLog(paramsVector); break;
1165   }
1166
1167   penaltyVal *= fgRegularizationWeight; //beta*PU
1168
1169   // Ad
1170   TVectorD measGuessVector(fgMaxInput);
1171   measGuessVector = (*fgCorrelationMatrix) * paramsVector;
1172
1173   // Ad - m
1174   measGuessVector -= (*fgCurrentESDVector);
1175
1176 #if 0
1177   // new error calcuation using error on the guess
1178
1179   // error from guess
1180   TVectorD errorGuessVector(fgMaxInput);
1181   //errorGuessVector = (*fgCorrelationMatrixSquared) * paramsVector;
1182   errorGuessVector = (*fgCorrelationMatrix) * paramsVector;
1183
1184   Double_t chi2FromFit = 0;
1185   for (Int_t i=0; i<fgMaxInput; ++i)
1186   {
1187     Float_t error = 1;
1188     if (errorGuessVector(i) > 0)
1189       error = errorGuessVector(i);
1190     chi2FromFit += measGuessVector(i) * measGuessVector(i) / error;
1191   }
1192
1193 #else
1194   // old error calcuation using the error on the measured
1195   TVectorD copy(measGuessVector);
1196
1197   // (Ad - m) W
1198   // this step can be optimized because currently only the diagonal elements of fgCorrelationCovarianceMatrix are used
1199   // normal way is like this:
1200   // measGuessVector *= (*fgCorrelationCovarianceMatrix);
1201   // optimized way like this:
1202   for (Int_t i=0; i<fgMaxInput; ++i)
1203     measGuessVector[i] *= (*fgCorrelationCovarianceMatrix)(i, i);
1204
1205
1206   if (fgSkipBin0InChi2)
1207     measGuessVector[0] = 0;
1208
1209   // (Ad - m) W (Ad - m)
1210   // the factor 1e6 prevents that a very small number (measGuessVector[i]) is multiplied with a very
1211   // big number ((*fgCorrelationCovarianceMatrix)(i, i)) (see UnfoldWithMinuit)
1212   Double_t chi2FromFit = measGuessVector * copy * 1e6;
1213 #endif
1214
1215   Double_t notFoundEventsConstraint = 0;
1216   Double_t currentNotFoundEvents = 0;
1217   Double_t errorNotFoundEvents = 0;
1218
1219   if (fgNotFoundEvents > 0)
1220   {
1221     for (Int_t i=0; i<fgMaxParams; ++i)
1222     {
1223       Double_t eff = (1.0 / (*fgEfficiency)(i) - 1);
1224       if (i == 0)
1225         eff = (1.0 / params[fgMaxParams] - 1);
1226       currentNotFoundEvents += eff * paramsVector(i);
1227       errorNotFoundEvents += eff * eff * paramsVector(i); // error due to guess (paramsVector)
1228       if (i <= 3)
1229         errorNotFoundEvents += (eff * 0.03) * (eff * 0.03) * paramsVector(i) * paramsVector(i); // error on eff
1230       //      if ((fgCallCount % 10000) == 0)
1231         //Printf("%d %f %f %f", i, (*fgEfficiency)(i), paramsVector(i), currentNotFoundEvents);
1232     }
1233     errorNotFoundEvents += fgNotFoundEvents;
1234     // TODO add error on background, background estimate
1235
1236     notFoundEventsConstraint = (currentNotFoundEvents - fgNotFoundEvents) * (currentNotFoundEvents - fgNotFoundEvents) / errorNotFoundEvents;
1237   }
1238
1239   Float_t avg = 0;
1240   Float_t sum = 0;
1241   Float_t currentMult = 0;
1242   // try to match dn/deta
1243   for (Int_t i=0; i<fgMaxParams; ++i)
1244   {
1245     avg += paramsVector(i) * currentMult;
1246     sum += paramsVector(i);
1247     currentMult += fgUnfoldedAxis->GetBinWidth(i);
1248   }
1249   avg /= sum;
1250   Float_t chi2Avg = 0; //(avg - 3.73) * (avg - 3.73) * 100;
1251
1252   chi2 = chi2FromFit + penaltyVal + notFoundEventsConstraint + chi2Avg;
1253
1254   if ((fgCallCount++ % 1000) == 0)
1255   {
1256
1257     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);
1258
1259     //for (Int_t i=0; i<fgMaxInput; ++i)
1260     //  Printf("%d: %f", i, measGuessVector(i) * copy(i) * 1e6);
1261   }
1262
1263   fChi2FromFit = chi2FromFit;
1264   fPenaltyVal = penaltyVal;
1265 }
1266
1267 //____________________________________________________________________
1268 void AliUnfolding::MakePads() {
1269     TPad *presult = new TPad("presult","result",0,0.4,1,1);
1270     presult->SetNumber(1);
1271     presult->Draw();
1272     presult->SetLogy();
1273     TPad *pres = new TPad("pres","residuals",0,0.2,1,0.4);
1274     pres->SetNumber(2);
1275     pres->Draw();
1276     TPad *ppen = new TPad("ppen","penalty",0,0.0,1,0.2);
1277     ppen->SetNumber(3);
1278     ppen->Draw();
1279
1280 }
1281 //____________________________________________________________________
1282 void AliUnfolding::DrawResults(TH2* correlation, TH1* efficiency, TH1* measured, TH1* initialConditions, TCanvas *canv, Int_t reuseHists,TH1 *unfolded)
1283 {
1284   // Draw histograms of
1285   //   - Result folded with response
1286   //   - Penalty factors
1287   //   - Chisquare contributions
1288   // (Useful for debugging/sanity checks and the interactive unfolder)
1289   //
1290   // If a canvas pointer is given (canv != 0), it will be used for all
1291   // plots; 3 pads are made if needed.
1292
1293
1294   Int_t blankCanvas = 0;
1295   TVirtualPad *presult = 0;
1296   TVirtualPad *pres = 0;
1297   TVirtualPad *ppen = 0;
1298
1299   if (canv) {
1300     canv->cd();
1301     presult = canv->GetPad(1);
1302     pres = canv->GetPad(2);
1303     ppen = canv->GetPad(3);
1304     if (presult == 0 || pres == 0 || ppen == 0) {
1305       canv->Clear();
1306       MakePads();
1307       blankCanvas = 1;
1308       presult = canv->GetPad(1);
1309       pres = canv->GetPad(2);
1310       ppen = canv->GetPad(3);
1311     }
1312   }
1313   else {
1314     presult = new TCanvas;
1315     pres = new TCanvas;
1316     ppen = new TCanvas;
1317   }
1318
1319
1320   if (fgMaxInput == -1)
1321   {
1322     Printf("AliUnfolding::Unfold: WARNING. Number of measured bins not set with SetNbins. Using number of bins in measured distribution");
1323     fgMaxInput = measured->GetNbinsX();
1324   }
1325   if (fgMaxParams == -1)
1326   {
1327     Printf("AliUnfolding::Unfold: WARNING. Number of unfolded bins not set with SetNbins. Using number of bins in measured distribution");
1328     fgMaxParams = initialConditions->GetNbinsX();
1329   }
1330
1331   if (fgOverflowBinLimit > 0)
1332     CreateOverflowBin(correlation, measured);
1333
1334   // Uses Minuit methods
1335
1336   SetStaticVariables(correlation, measured, efficiency);
1337
1338   Double_t* params = new Double_t[fgMaxParams+1];
1339   for (Int_t i=0; i<fgMaxParams; ++i)
1340   {
1341     params[i] = initialConditions->GetBinContent(i+1) * efficiency->GetBinContent(i+1);
1342
1343     Bool_t fix = kFALSE;
1344     if (params[i] < 0)
1345     {
1346       fix = kTRUE;
1347       params[i] = -params[i];
1348     }
1349     params[i]=TMath::Sqrt(params[i]);
1350
1351     //cout << "params[" << i << "] " << params[i] << endl;
1352
1353   }
1354
1355   Double_t chi2;
1356   Int_t dummy;
1357
1358   //fgPrintChi2Details = kTRUE;
1359   fgCallCount = 0; // To make sure that Chi2Function prints the components
1360   Chi2Function(dummy, 0, chi2, params, 0);
1361
1362   presult->cd();
1363   TH1 *meas2 = (TH1*)measured->Clone("meas2");
1364   meas2->SetTitle("meas2");
1365   meas2->SetName("meas2");
1366   meas2->SetLineColor(1);
1367   meas2->SetLineWidth(2);
1368   meas2->SetMarkerColor(meas2->GetLineColor());
1369   meas2->SetMarkerStyle(20);
1370   Float_t scale = unfolded->GetXaxis()->GetBinWidth(1)/meas2->GetXaxis()->GetBinWidth(1);
1371   meas2->Scale(scale);
1372   if (blankCanvas)
1373     meas2->DrawCopy();
1374   else
1375     meas2->DrawCopy("same");
1376   //meas2->GetXaxis()->SetLimits(0,fgMaxInput);
1377   meas2->SetBit(kCannotPick);
1378   DrawGuess(params, presult, pres, ppen, reuseHists,unfolded);
1379   delete [] params;
1380 }
1381 //____________________________________________________________________
1382 void AliUnfolding::RedrawInteractive() {
1383   //
1384   // Helper function for interactive unfolding
1385   //
1386   DrawResults(fghCorrelation,fghEfficiency,fghMeasured,fghUnfolded,fgCanvas,1,fghUnfolded);
1387 }
1388 //____________________________________________________________________
1389 void AliUnfolding::InteractiveUnfold(TH2* correlation, TH1* efficiency, TH1* measured, TH1* initialConditions)
1390 {
1391   //
1392   // Function to do interactive unfolding
1393   // A canvas is drawn with the unfolding result
1394   // Change the histogram with your mouse and all histograms
1395   // will be updated automatically
1396
1397   fgCanvas = new TCanvas("UnfoldingCanvas","Interactive unfolding",500,800);
1398   fgCanvas->cd();
1399
1400   MakePads();
1401
1402   if (fghUnfolded) {
1403     delete fghUnfolded; fghUnfolded = 0;
1404   }
1405
1406   gROOT->SetEditHistograms(kTRUE);
1407
1408   fghUnfolded = new TH1F("AliUnfoldingfghUnfolded","Unfolded result (Interactive unfolder",efficiency->GetNbinsX(),efficiency->GetXaxis()->GetXmin(),efficiency->GetXaxis()->GetXmax());
1409
1410   fghCorrelation = correlation;
1411   fghEfficiency = efficiency;
1412   fghMeasured = measured;
1413
1414   if(fgMinuitMaxIterations>0)
1415     UnfoldWithMinuit(correlation,efficiency,measured,initialConditions,fghUnfolded,kFALSE);
1416   else
1417     fghUnfolded = initialConditions;
1418
1419   fghUnfolded->SetLineColor(2);
1420   fghUnfolded->SetMarkerColor(2);
1421   fghUnfolded->SetLineWidth(2);
1422
1423
1424   fgCanvas->cd(1);
1425   fghUnfolded->Draw();
1426   DrawResults(correlation,efficiency,measured,fghUnfolded,fgCanvas,kFALSE,fghUnfolded);
1427
1428   TExec *execRedraw = new TExec("redraw","AliUnfolding::RedrawInteractive()");
1429   fghUnfolded->GetListOfFunctions()->Add(execRedraw);
1430 }
1431 //____________________________________________________________________
1432 void AliUnfolding::DrawGuess(Double_t *params, TVirtualPad *pfolded, TVirtualPad *pres, TVirtualPad *ppen, Int_t reuseHists,TH1* unfolded)
1433 {
1434   //
1435   // draws residuals of solution suggested by params and effect of regularization
1436   //
1437
1438   if (pfolded == 0)
1439     pfolded = new TCanvas;
1440   if (ppen == 0)
1441     ppen = new TCanvas;
1442   if (pres == 0)
1443     pres = new TCanvas;
1444
1445   // d
1446   TVectorD paramsVector(fgMaxParams);
1447   for (Int_t i=0; i<fgMaxParams; ++i)
1448     paramsVector[i] = params[i] * params[i];
1449
1450   // Ad
1451   TVectorD measGuessVector(fgMaxInput);
1452   measGuessVector = (*fgCorrelationMatrix) * paramsVector;
1453
1454   TH1* folded = 0;
1455   if (reuseHists)
1456     folded = dynamic_cast<TH1*>(gROOT->FindObject("__hfolded"));
1457   if (!reuseHists || folded == 0) {
1458     if (fgMeasuredAxis->GetXbins()->GetArray()) // variable bins
1459       folded = new TH1F("__hfolded","Folded histo from AliUnfolding",fgMeasuredAxis->GetNbins(),fgMeasuredAxis->GetXbins()->GetArray());
1460     else
1461       folded = new TH1F("__hfolded","Folded histo from AliUnfolding",fgMeasuredAxis->GetNbins(),fgMeasuredAxis->GetXmin(),fgMeasuredAxis->GetXmax());
1462   }
1463
1464   folded->SetBit(kCannotPick);
1465   folded->SetLineColor(4);
1466   folded->SetLineWidth(2);
1467
1468   for (Int_t ibin =0; ibin < fgMaxInput; ibin++)
1469     folded->SetBinContent(ibin+1, measGuessVector[ibin]);
1470
1471   Float_t scale = unfolded->GetXaxis()->GetBinWidth(1)/folded->GetXaxis()->GetBinWidth(1);
1472   folded->Scale(scale);
1473
1474   pfolded->cd();
1475
1476   folded->Draw("same");
1477
1478   // Ad - m
1479   measGuessVector -= (*fgCurrentESDVector);
1480
1481   TVectorD copy(measGuessVector);
1482
1483   // (Ad - m) W
1484   // this step can be optimized because currently only the diagonal elements of fgCorrelationCovarianceMatrix are used
1485   // normal way is like this:
1486   // measGuessVector *= (*fgCorrelationCovarianceMatrix);
1487   // optimized way like this:
1488   for (Int_t i=0; i<fgMaxInput; ++i)
1489     measGuessVector[i] *= (*fgCorrelationCovarianceMatrix)(i, i);
1490
1491   // (Ad - m) W (Ad - m)
1492   // the factor 1e6 prevents that a very small number (measGuessVector[i]) is multiplied with a very
1493   // big number ((*fgCorrelationCovarianceMatrix)(i, i)) (see ApplyMinuitFit)
1494   //Double_t chi2FromFit = measGuessVector * copy * 1e6;
1495
1496   // draw residuals
1497   // Double_t pTarray[fgMaxParams+1];
1498   // for(int i=0; i<fgMaxParams; i++) {
1499   //   pTarray[i] = fgUnfoldedAxis->GetBinCenter(i)-0.5*fgUnfoldedAxis->GetBinWidth(i);
1500   // }
1501   //  pTarray[fgMaxParams] = fgUnfoldedAxis->GetBinCenter(fgMaxParams-1)+0.5*fgUnfoldedAxis->GetBinWidth(fgMaxParams-1);
1502   //  TH1* residuals = new TH1F("residuals", "residuals", fgMaxParams,pTarray);
1503   //  TH1* residuals = new TH1F("residuals", "residuals", fgMaxInput, -0.5, fgMaxInput - 0.5);
1504   // for (Int_t i=0; i<fgMaxInput; ++i)
1505   //   residuals->SetBinContent(i+1, measGuessVector(i) * copy(i) * 1e6);7
1506   TH1* residuals = GetResidualsPlot(params);
1507   residuals->GetXaxis()->SetTitleSize(0.06);
1508   residuals->GetXaxis()->SetTitleOffset(0.7);
1509   residuals->GetXaxis()->SetLabelSize(0.07);
1510   residuals->GetYaxis()->SetTitleSize(0.08);
1511   residuals->GetYaxis()->SetTitleOffset(0.5);
1512   residuals->GetYaxis()->SetLabelSize(0.06);
1513   pres->cd(); residuals->DrawCopy(); gPad->SetLogy();
1514
1515
1516   // draw penalty
1517   TH1* penalty = GetPenaltyPlot(params);
1518   penalty->GetXaxis()->SetTitleSize(0.06);
1519   penalty->GetXaxis()->SetTitleOffset(0.7);
1520   penalty->GetXaxis()->SetLabelSize(0.07);
1521   penalty->GetYaxis()->SetTitleSize(0.08);
1522   penalty->GetYaxis()->SetTitleOffset(0.5);
1523   penalty->GetYaxis()->SetLabelSize(0.06);
1524   ppen->cd(); penalty->DrawCopy(); gPad->SetLogy();
1525 }
1526
1527 //____________________________________________________________________
1528 TH1* AliUnfolding::GetResidualsPlot(TH1* corrected)
1529 {
1530   //
1531   // MvL: THIS MUST BE INCORRECT.
1532   // Need heff to calculate params from TH1 'corrected'
1533   //
1534
1535   //
1536   // fill residuals histogram of solution suggested by params and effect of regularization
1537   //
1538
1539   Double_t* params = new Double_t[fgMaxParams];
1540   for (Int_t i=0; i<fgMaxParams; i++)
1541     params[i] = 0;
1542
1543   for (Int_t i=0; i<TMath::Min(fgMaxParams, corrected->GetNbinsX()); i++)
1544     params[i] = TMath::Sqrt(TMath::Abs(corrected->GetBinContent(i+1)*(*fgEfficiency)(i)));
1545
1546   TH1 * plot = GetResidualsPlot(params);
1547   delete [] params;
1548   return plot;
1549 }
1550
1551 //____________________________________________________________________
1552 TH1* AliUnfolding::GetResidualsPlot(Double_t* params)
1553 {
1554   //
1555   // fill residuals histogram of solution suggested by params and effect of regularization
1556   //
1557
1558   // d
1559   TVectorD paramsVector(fgMaxParams);
1560   for (Int_t i=0; i<fgMaxParams; ++i)
1561     paramsVector[i] = params[i] * params[i];
1562
1563   // Ad
1564   TVectorD measGuessVector(fgMaxInput);
1565   measGuessVector = (*fgCorrelationMatrix) * paramsVector;
1566
1567   // Ad - m
1568   measGuessVector -= (*fgCurrentESDVector);
1569
1570   TVectorD copy(measGuessVector);
1571
1572   // (Ad - m) W
1573   // this step can be optimized because currently only the diagonal elements of fgCorrelationCovarianceMatrix are used
1574   // normal way is like this:
1575   // measGuessVector *= (*fgCorrelationCovarianceMatrix);
1576   // optimized way like this:
1577   for (Int_t i=0; i<fgMaxInput; ++i)
1578     measGuessVector[i] *= (*fgCorrelationCovarianceMatrix)(i, i);
1579
1580   // (Ad - m) W (Ad - m)
1581   // the factor 1e6 prevents that a very small number (measGuessVector[i]) is multiplied with a very
1582   // big number ((*fgCorrelationCovarianceMatrix)(i, i)) (see ApplyMinuitFit)
1583   //Double_t chi2FromFit = measGuessVector * copy * 1e6;
1584
1585   // draw residuals
1586   TH1* residuals = 0;
1587   if (fgMeasuredAxis->GetXbins()->GetArray()) // variable bins
1588     residuals = new TH1F("residuals", "residuals;unfolded pos;residual",fgMeasuredAxis->GetNbins(),fgMeasuredAxis->GetXbins()->GetArray());
1589   else
1590     residuals = new TH1F("residuals", "residuals;unfolded pos;residual",fgMeasuredAxis->GetNbins(),fgMeasuredAxis->GetXmin(), fgMeasuredAxis->GetXmax());
1591   //  TH1* residuals = new TH1F("residuals", "residuals", fgMaxInput, -0.5, fgMaxInput - 0.5);
1592
1593   Double_t sumResiduals = 0.;
1594   for (Int_t i=0; i<fgMaxInput; ++i) {
1595     residuals->SetBinContent(i+1, measGuessVector(i) * copy(i) * 1e6);
1596     sumResiduals += measGuessVector(i) * copy(i) * 1e6;
1597   }
1598   fAvgResidual = sumResiduals/(double)fgMaxInput;
1599
1600   return residuals;
1601 }
1602
1603 //____________________________________________________________________
1604 TH1* AliUnfolding::GetPenaltyPlot(TH1* corrected)
1605 {
1606   // draws the penalty factors as function of multiplicity of the current selected regularization
1607
1608   Double_t* params = new Double_t[fgMaxParams];
1609   for (Int_t i=0; i<fgMaxParams; i++)
1610     params[i] = 0;
1611
1612   for (Int_t i=0; i<TMath::Min(fgMaxParams, corrected->GetNbinsX()); i++)
1613     params[i] = (*fgEfficiency)(i)*corrected->GetBinContent(i+1);
1614
1615   TH1* penalty = GetPenaltyPlot(params);
1616
1617   delete[] params;
1618
1619   return penalty;
1620 }
1621
1622 //____________________________________________________________________
1623 TH1* AliUnfolding::GetPenaltyPlot(Double_t* params)
1624 {
1625   // draws the penalty factors as function of multiplicity of the current selected regularization
1626
1627   //TH1* penalty = new TH1F("penalty", ";unfolded multiplicity;penalty factor", fgMaxParams, -0.5, fgMaxParams - 0.5);
1628   //  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) );
1629
1630   TH1* penalty = 0;
1631   if (fgUnfoldedAxis->GetXbins()->GetArray())
1632     penalty = new TH1F("penalty", ";unfolded pos;penalty factor", fgUnfoldedAxis->GetNbins(),fgUnfoldedAxis->GetXbins()->GetArray());
1633   else
1634     penalty = new TH1F("penalty", ";unfolded pos;penalty factor", fgUnfoldedAxis->GetNbins(),fgUnfoldedAxis->GetXmin(),fgUnfoldedAxis->GetXmax());
1635
1636   for (Int_t i=1+fgSkipBinsBegin; i<fgMaxParams; ++i)
1637   {
1638     Double_t diff = 0;
1639     if (fgRegularizationType == kPol0)
1640     {
1641       Double_t right  = params[i] / fgUnfoldedAxis->GetBinWidth(i+1);
1642       Double_t left   = params[i-1] / fgUnfoldedAxis->GetBinWidth(i);
1643
1644       if (left != 0)
1645       {
1646         Double_t diffTmp = (right - left);
1647         diff = diffTmp * diffTmp / left / ((fgUnfoldedAxis->GetBinWidth(i+1) + fgUnfoldedAxis->GetBinWidth(i)) / 2) / 100;
1648       }
1649     }
1650     if (fgRegularizationType == kPol1 && i > 1+fgSkipBinsBegin)
1651     {
1652       if (params[i-1] == 0)
1653         continue;
1654
1655       Double_t right  = params[i] / fgUnfoldedAxis->GetBinWidth(i+1);
1656       Double_t middle = params[i-1] / fgUnfoldedAxis->GetBinWidth(i);
1657       Double_t left   = params[i-2] / fgUnfoldedAxis->GetBinWidth(i-1);
1658
1659       Double_t der1 = (right - middle) / ((fgUnfoldedAxis->GetBinWidth(i+1) + fgUnfoldedAxis->GetBinWidth(i)) / 2);
1660       Double_t der2 = (middle - left) / ((fgUnfoldedAxis->GetBinWidth(i) + fgUnfoldedAxis->GetBinWidth(i-1)) / 2);
1661
1662       diff = (der1 - der2) * (der1 - der2) / middle;
1663     }
1664
1665     if (fgRegularizationType == kLog && i > 1+fgSkipBinsBegin)
1666     {
1667       if (params[i-1] == 0)
1668         continue;
1669
1670       Double_t right  = log(params[i]);
1671       Double_t middle = log(params[i-1]);
1672       Double_t left   = log(params[i-2]);
1673
1674       Double_t der1 = (right - middle);
1675       Double_t der2 = (middle - left);
1676
1677       //Double_t error = 1. / params[i] + 4. / params[i-1] + 1. / params[i-2];
1678       //Printf("%d %f %f", i, (der1 - der2) * (der1 - der2), error);
1679
1680       diff = (der1 - der2) * (der1 - der2);// / error;
1681     }
1682     if (fgRegularizationType == kCurvature && i > 1+fgSkipBinsBegin)
1683     {
1684       Double_t right  = params[i];    // params are sqrt
1685       Double_t middle = params[i-1];
1686       Double_t left   = params[i-2];
1687
1688       Double_t der1 = (right - middle)/0.5/(fgUnfoldedAxis->GetBinWidth(i-1) + fgUnfoldedAxis->GetBinWidth(i));
1689       Double_t der2 = (middle - left)/0.5/(fgUnfoldedAxis->GetBinWidth(i) + fgUnfoldedAxis->GetBinWidth(i+1));
1690
1691       diff = (der1 - der2)/(fgUnfoldedAxis->GetBinWidth(i-1) + fgUnfoldedAxis->GetBinWidth(i) + fgUnfoldedAxis->GetBinWidth(i-1))*3.0;
1692       diff = 1e4*diff*diff;
1693     }
1694     if (fgRegularizationType == kPowerLaw && i > 1+fgSkipBinsBegin)
1695     {
1696
1697       if (params[i] < 1e-8 || params[i-1] < 1e-8 || params[i-2] < 1e-8)
1698         continue;
1699
1700       if (fgUnfoldedAxis->GetBinWidth(i+1) < 1e-8 || fgUnfoldedAxis->GetBinWidth(i) < 1e-8 || fgUnfoldedAxis->GetBinWidth(i) < 1e-8)
1701         continue;
1702
1703       double middle = TMath::Power(params[i-1] / fgUnfoldedAxis->GetBinWidth(i),fgPowern);
1704
1705       if(middle>0) {
1706         double right  = TMath::Power(params[i] / fgUnfoldedAxis->GetBinWidth(i+1),fgPowern)/middle;
1707
1708         double left   = TMath::Power(params[i-2] / fgUnfoldedAxis->GetBinWidth(i-1),fgPowern)/middle;
1709
1710         middle = 1.;
1711
1712         Double_t der1 = (right - middle) / ((fgUnfoldedAxis->GetBinWidth(i+1) + fgUnfoldedAxis->GetBinWidth(i)) / 2);
1713         Double_t der2 = (middle - left) / ((fgUnfoldedAxis->GetBinWidth(i) + fgUnfoldedAxis->GetBinWidth(i-1)) / 2);
1714
1715         diff = (der1 - der2) * (der1 - der2);// / error;
1716       }
1717     }
1718
1719     if (fgRegularizationType == kLogLog && i > 1+fgSkipBinsBegin)
1720     {
1721
1722       if (params[i] < 1e-8 || params[i-1] < 1e-8 || params[i-2] < 1e-8)
1723         continue;
1724
1725       Double_t right  = log(params[i] / (*fgEfficiency)(i) / fgUnfoldedAxis->GetBinWidth(i+1));
1726       Double_t middle = log(params[i-1] / (*fgEfficiency)(i-1) / fgUnfoldedAxis->GetBinWidth(i));
1727       Double_t left   = log(params[i-2] / (*fgEfficiency)(i-2) / fgUnfoldedAxis->GetBinWidth(i-1));
1728
1729       Double_t der1 = (right - middle) / ( log(fgUnfoldedAxis->GetBinCenter(i+1)) - log(fgUnfoldedAxis->GetBinCenter(i)) );
1730       Double_t der2 = (middle - left) /( log(fgUnfoldedAxis->GetBinCenter(i)) - log(fgUnfoldedAxis->GetBinCenter(i-1)) );
1731
1732       double tmp = (log(fgUnfoldedAxis->GetBinCenter(i+1)) - log(fgUnfoldedAxis->GetBinCenter(i-1)))/2.;
1733       Double_t dder = (der1-der2) / tmp;
1734
1735       diff = dder * dder;
1736     }
1737
1738     penalty->SetBinContent(i, diff*fgRegularizationWeight);
1739
1740     //Printf("%d %f %f %f %f", i-1, left, middle, right, diff);
1741   }
1742
1743   return penalty;
1744 }
1745
1746 //____________________________________________________________________
1747 void AliUnfolding::TF1Function(Int_t& unused1, Double_t* unused2, Double_t& chi2, Double_t *params, Int_t unused3)
1748 {
1749   //
1750   // fit function for minuit
1751   // uses the TF1 stored in fgFitFunction
1752   //
1753
1754   for (Int_t i=0; i<fgFitFunction->GetNpar(); i++)
1755     fgFitFunction->SetParameter(i, params[i]);
1756
1757   Double_t* params2 = new Double_t[fgMaxParams];
1758
1759   for (Int_t i=0; i<fgMaxParams; ++i)
1760     params2[i] = fgFitFunction->Eval(i);
1761
1762   Chi2Function(unused1, unused2, chi2, params2, unused3);
1763
1764   delete[] params2;
1765
1766   if (fgDebug)
1767     Printf("%f", chi2);
1768 }
1769
1770 //____________________________________________________________________
1771 Int_t AliUnfolding::UnfoldWithFunction(TH2* correlation, TH1* efficiency, TH1* measured, TH1* /* initialConditions */, TH1* aResult)
1772 {
1773   //
1774   // correct spectrum using minuit chi2 method applying a functional fit
1775   //
1776
1777   if (!fgFitFunction)
1778   {
1779     Printf("AliUnfolding::UnfoldWithFunction: ERROR fit function not set. Exiting.");
1780     return -1;
1781   }
1782
1783   SetChi2Regularization(kNone, 0);
1784
1785   SetStaticVariables(correlation, measured, efficiency);
1786
1787   // Initialize TMinuit via generic fitter interface
1788   TVirtualFitter *minuit = TVirtualFitter::Fitter(0, fgFitFunction->GetNpar());
1789
1790   minuit->SetFCN(TF1Function);
1791   for (Int_t i=0; i<fgFitFunction->GetNpar(); i++)
1792   {
1793     Double_t lower, upper;
1794     fgFitFunction->GetParLimits(i, lower, upper);
1795     minuit->SetParameter(i, Form("param%d",i), fgFitFunction->GetParameter(i), fgMinuitStepSize, lower, upper);
1796   }
1797
1798   Double_t arglist[100];
1799   arglist[0] = 0;
1800   minuit->ExecuteCommand("SET PRINT", arglist, 1);
1801   minuit->ExecuteCommand("SCAN", arglist, 0);
1802   minuit->ExecuteCommand("MIGRAD", arglist, 0);
1803   //minuit->ExecuteCommand("MINOS", arglist, 0);
1804
1805   for (Int_t i=0; i<fgFitFunction->GetNpar(); i++)
1806     fgFitFunction->SetParameter(i, minuit->GetParameter(i));
1807
1808   for (Int_t i=0; i<fgMaxParams; ++i)
1809   {
1810     Double_t value = fgFitFunction->Eval(i);
1811     if (fgDebug)
1812       Printf("%d : %f", i, value);
1813
1814     if (efficiency)
1815     {
1816       if (efficiency->GetBinContent(i+1) > 0)
1817       {
1818         value /= efficiency->GetBinContent(i+1);
1819       }
1820       else
1821         value = 0;
1822     }
1823     aResult->SetBinContent(i+1, value);
1824     aResult->SetBinError(i+1, 0);
1825   }
1826
1827   return 0;
1828 }
1829
1830 //____________________________________________________________________
1831 void AliUnfolding::CreateOverflowBin(TH2* correlation, TH1* measured)
1832 {
1833   // Finds the first bin where the content is below fgStatLimit and combines all values for this bin and larger bins
1834   // The same limit is applied to the correlation
1835
1836   Int_t lastBin = 0;
1837   for (Int_t i=1; i<=measured->GetNbinsX(); ++i)
1838   {
1839     if (measured->GetBinContent(i) <= fgOverflowBinLimit)
1840     {
1841       lastBin = i;
1842       break;
1843     }
1844   }
1845
1846   if (lastBin == 0)
1847   {
1848     Printf("AliUnfolding::CreateOverflowBin: WARNING: First bin is already below limit of %f", fgOverflowBinLimit);
1849     return;
1850   }
1851
1852   Printf("AliUnfolding::CreateOverflowBin: Bin limit in measured spectrum determined to be %d", lastBin);
1853
1854   TCanvas* canvas = 0;
1855
1856   if (fgDebug)
1857   {
1858     canvas = new TCanvas("StatSolution", "StatSolution", 1000, 800);
1859     canvas->Divide(2, 2);
1860
1861     canvas->cd(1);
1862     measured->SetStats(kFALSE);
1863     measured->DrawCopy();
1864     gPad->SetLogy();
1865
1866     canvas->cd(2);
1867     correlation->SetStats(kFALSE);
1868     correlation->DrawCopy("COLZ");
1869   }
1870
1871   measured->SetBinContent(lastBin, measured->Integral(lastBin, measured->GetNbinsX()));
1872   for (Int_t i=lastBin+1; i<=measured->GetNbinsX(); ++i)
1873   {
1874     measured->SetBinContent(i, 0);
1875     measured->SetBinError(i, 0);
1876   }
1877   // the error is set to sqrt(N), better solution possible?, sum of relative errors of all contributions???
1878   measured->SetBinError(lastBin, TMath::Sqrt(measured->GetBinContent(lastBin)));
1879
1880   Printf("AliUnfolding::CreateOverflowBin: This bin has now %f +- %f entries", measured->GetBinContent(lastBin), measured->GetBinError(lastBin));
1881
1882   for (Int_t i=1; i<=correlation->GetNbinsX(); ++i)
1883   {
1884     correlation->SetBinContent(i, lastBin, correlation->Integral(i, i, lastBin, correlation->GetNbinsY()));
1885     // the error is set to sqrt(N), better solution possible?, sum of relative errors of all contributions???
1886     correlation->SetBinError(i, lastBin, TMath::Sqrt(correlation->GetBinContent(i, lastBin)));
1887
1888     for (Int_t j=lastBin+1; j<=correlation->GetNbinsY(); ++j)
1889     {
1890       correlation->SetBinContent(i, j, 0);
1891       correlation->SetBinError(i, j, 0);
1892     }
1893   }
1894
1895   Printf("AliUnfolding::CreateOverflowBin: Adjusted correlation matrix!");
1896
1897   if (canvas)
1898   {
1899     canvas->cd(3);
1900     measured->DrawCopy();
1901     gPad->SetLogy();
1902
1903     canvas->cd(4);
1904     correlation->DrawCopy("COLZ");
1905   }
1906 }
1907
1908 Int_t AliUnfolding::UnfoldGetBias(TH2* correlation, TH1* efficiency, TH1* measured, TH1* initialConditions, TH1* result)
1909 {
1910   // unfolds and assigns bias as errors with Eq. 19 of Cowan, "a survey of unfolding methods for particle physics"
1911   // b_i = sum_j dmu_i/dn_j (nu_j - n_j) with nu_j as folded guess, n_j as data
1912   // dmu_i/dn_j is found numerically by changing the bin content and re-unfolding
1913   //
1914   // return code: 0 (success) -1 (error: from Unfold(...) )
1915
1916   if (Unfold(correlation, efficiency, measured, initialConditions, result) != 0)
1917     return -1;
1918
1919   TMatrixD matrix(fgMaxInput, fgMaxParams);
1920
1921   TH1* newResult[4];
1922   for (Int_t loop=0; loop<4; loop++)
1923     newResult[loop] = (TH1*) result->Clone(Form("newresult_%d", loop));
1924
1925   // change bin-by-bin and built matrix of effects
1926   for (Int_t m=0; m<fgMaxInput; m++)
1927   {
1928     if (measured->GetBinContent(m+1) < 1)
1929       continue;
1930
1931     for (Int_t loop=0; loop<4; loop++)
1932     {
1933       //Printf("%d %d", i, loop);
1934       Float_t factor = 1;
1935       switch (loop)
1936       {
1937         case 0: factor = 0.5; break;
1938         case 1: factor = -0.5; break;
1939         case 2: factor = 1; break;
1940         case 3: factor = -1; break;
1941         default: return -1;
1942       }
1943
1944       TH1* measuredClone = (TH1*) measured->Clone("measuredClone");
1945
1946       measuredClone->SetBinContent(m+1, measured->GetBinContent(m+1) + factor * TMath::Sqrt(measured->GetBinContent(m+1)));
1947       //new TCanvas; measuredClone->Draw("COLZ");
1948
1949       newResult[loop]->Reset();
1950       Unfold(correlation, efficiency, measuredClone, measuredClone, newResult[loop]);
1951       // WARNING if we leave here, then nothing is calculated
1952         //return -1;
1953
1954       delete measuredClone;
1955     }
1956
1957     for (Int_t t=0; t<fgMaxParams; t++)
1958     {
1959       // one value
1960       //matrix(i, j) = (result->GetBinContent(j+1) - newResult->GetBinContent(j+1)) / TMath::Sqrt(changedHist->GetBinContent(1, i+1));
1961
1962       // four values from the derivate (procedure taken from ROOT -- suggestion by Ruben)
1963       // = 1/2D * [ 8 (f(D/2) - f(-D/2)) - (f(D)-f(-D)) ]/3
1964
1965       /*
1966       // check formula
1967       measured->SetBinContent(1, m+1, 100);
1968       newResult[0]->SetBinContent(t+1, 5 * 0.5 + 10);
1969       newResult[1]->SetBinContent(t+1, 5 * -0.5 + 10);
1970       newResult[2]->SetBinContent(t+1, 5 * 1 + 10);
1971       newResult[3]->SetBinContent(t+1, 5 * -1 + 10);
1972       */
1973
1974       matrix(m, t) = 0.5 / TMath::Sqrt(measured->GetBinContent(m+1)) *
1975         ( 8. * (newResult[0]->GetBinContent(t+1) - newResult[1]->GetBinContent(t+1)) -
1976               (newResult[2]->GetBinContent(t+1) - newResult[3]->GetBinContent(t+1))
1977         ) / 3;
1978     }
1979   }
1980
1981   for (Int_t loop=0; loop<4; loop++)
1982     delete newResult[loop];
1983
1984   // ...calculate folded guess
1985   TH1* convoluted = (TH1*) measured->Clone("convoluted");
1986   convoluted->Reset();
1987   for (Int_t m=0; m<fgMaxInput; m++)
1988   {
1989     Float_t value = 0;
1990     for (Int_t t = 0; t<fgMaxParams; t++)
1991     {
1992       Float_t tmp = correlation->GetBinContent(t+1, m+1) * result->GetBinContent(t+1);
1993       if (efficiency)
1994         tmp *= efficiency->GetBinContent(t+1);
1995       value += tmp;
1996     }
1997     convoluted->SetBinContent(m+1, value);
1998   }
1999
2000   // rescale
2001   convoluted->Scale(measured->Integral() / convoluted->Integral());
2002
2003   //new TCanvas; convoluted->Draw(); measured->Draw("SAME"); measured->SetLineColor(2);
2004   //return;
2005
2006   // difference
2007   convoluted->Add(measured, -1);
2008
2009   // ...and the bias
2010   for (Int_t t = 0; t<fgMaxParams; t++)
2011   {
2012     Double_t error = 0;
2013     for (Int_t m=0; m<fgMaxInput; m++)
2014       error += matrix(m, t) * convoluted->GetBinContent(m+1);
2015     result->SetBinError(t+1, error);
2016   }
2017
2018   //new TCanvas; result->Draw(); gPad->SetLogy();
2019
2020   return 0;
2021 }