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