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