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