PWG0/AliUnfolding.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id: AliUnfolding.cxx 31168 2009-02-23 15:18:45Z jgrosseo $ */
  17
  18 // This class allows 1-dimensional unfolding.
  19 // Methods that are implemented are chi2 minimization and bayesian unfolding.
  20 //
  21 //  Author: Jan.Fiete.Grosse-Oetringhaus@cern.ch
  22
  23 #include "AliUnfolding.h"
  24 #include <TH1F.h>
  25 #include <TH2F.h>
  26 #include <TVirtualFitter.h>
  27 #include <TMath.h>
  28 #include <TCanvas.h>
  29 #include <TF1.h>
  30
  31 TMatrixD* AliUnfolding::fgCorrelationMatrix = 0;
  32 TMatrixD* AliUnfolding::fgCorrelationCovarianceMatrix = 0;
  33 TVectorD* AliUnfolding::fgCurrentESDVector = 0;
  34 TVectorD* AliUnfolding::fgEntropyAPriori = 0;
  35
  36 TF1* AliUnfolding::fgFitFunction = 0;
  37
  38 AliUnfolding::MethodType AliUnfolding::fgMethodType = AliUnfolding::kInvalid;
  39 Int_t AliUnfolding::fgMaxInput  = -1;  // bins in measured histogram
  40 Int_t AliUnfolding::fgMaxParams = -1;  // bins in unfolded histogram = number of fit params
  41 Float_t AliUnfolding::fgOverflowBinLimit = -1;
  42
  43 AliUnfolding::RegularizationType AliUnfolding::fgRegularizationType = AliUnfolding::kPol1;
  44 Float_t AliUnfolding::fgRegularizationWeight = 10000;
  45 Int_t AliUnfolding::fgSkipBinsBegin = 0;
  46 Float_t AliUnfolding::fgMinuitStepSize = 0.1;                 // (usually not needed to be changed) step size in minimization
  47 Bool_t AliUnfolding::fgNormalizeInput = kFALSE;                  // normalize input spectrum
  48
  49 Float_t AliUnfolding::fgBayesianSmoothing  = 1;           // smoothing parameter (0 = no smoothing)
  50 Int_t   AliUnfolding::fgBayesianIterations = 10;          // number of iterations in Bayesian method
  51
  52 Bool_t AliUnfolding::fgDebug = kFALSE;
  53
  54 ClassImp(AliUnfolding)
  55
  56 //____________________________________________________________________
  57 void AliUnfolding::SetUnfoldingMethod(MethodType methodType)
  58 {
  59   // set unfolding method
  60   fgMethodType = methodType;
  61
  62   const char* name = 0;
  63   switch (methodType)
  64   {
  65     case kInvalid: name = "INVALID"; break;
  66     case kChi2Minimization: name = "Chi2 Minimization"; break;
  67     case kBayesian: name = "Bayesian unfolding"; break;
  68     case kFunction: name = "Functional fit"; break;
  69   }
  70   Printf("AliUnfolding::SetUnfoldingMethod: %s enabled.", name);
  71 }
  72
  73 //____________________________________________________________________
  74 void AliUnfolding::SetCreateOverflowBin(Float_t overflowBinLimit)
  75 {
  76   // enable the creation of a overflow bin that includes all statistics below the given limit
  77
  78   fgOverflowBinLimit = overflowBinLimit;
  79
  80   Printf("AliUnfolding::SetCreateOverflowBin: overflow bin limit set to %f", overflowBinLimit);
  81 }
  82
  83 //____________________________________________________________________
  84 void AliUnfolding::SetSkipBinsBegin(Int_t nBins)
  85 {
  86   // set number of skipped bins in regularization
  87
  88   fgSkipBinsBegin = nBins;
  89
  90   Printf("AliUnfolding::SetSkipBinsBegin: skipping %d bins at the beginning of the spectrum in the regularization.", fgSkipBinsBegin);
  91 }
  92
  93 //____________________________________________________________________
  94 void AliUnfolding::SetNbins(Int_t nMeasured, Int_t nUnfolded)
  95 {
  96   // set number of bins in the input (measured) distribution and in the unfolded distribution
  97   fgMaxInput = nMeasured;
  98   fgMaxParams = nUnfolded;
  99
 100   Printf("AliUnfolding::SetNbins: Set %d measured bins and %d unfolded bins", nMeasured, nUnfolded);
 101 }
 102
 103 //____________________________________________________________________
 104 void AliUnfolding::SetChi2Regularization(RegularizationType type, Float_t weight)
 105 {
 106   //
 107   // sets the parameters for chi2 minimization
 108   //
 109
 110   fgRegularizationType = type;
 111   fgRegularizationWeight = weight;
 112
 113   Printf("AliUnfolding::SetChi2Regularization --> Regularization set to %d with weight %f", (Int_t) type, weight);
 114 }
 115
 116 //____________________________________________________________________
 117 void AliUnfolding::SetBayesianParameters(Float_t smoothing, Int_t nIterations)
 118 {
 119   //
 120   // sets the parameters for Bayesian unfolding
 121   //
 122
 123   fgBayesianSmoothing = smoothing;
 124   fgBayesianIterations = nIterations;
 125
 126   Printf("AliUnfolding::SetBayesianParameters --> Paramaters set to %d iterations with smoothing %f", fgBayesianIterations, fgBayesianSmoothing);
 127 }
 128
 129 //____________________________________________________________________
 130 void AliUnfolding::SetFunction(TF1* function)
 131 {
 132   // set function for unfolding with a fit function
 133
 134   fgFitFunction = function;
 135
 136   Printf("AliUnfolding::SetFunction: Set fit function with %d parameters.", function->GetNpar());
 137 }
 138
 139 //____________________________________________________________________
 140 Int_t AliUnfolding::Unfold(TH2* correlation, TH1* efficiency, TH1* measured, TH1* initialConditions, TH1* result, Bool_t check)
 141 {
 142   // unfolds with unfolding method fgMethodType
 143   //
 144   // parameters:
 145   //  correlation: response matrix as measured vs. generated
 146   //  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.
 147   //  measured:    the measured spectrum
 148   //  initialConditions: (optional) initial conditions for the unfolding. if 0 the measured spectrum is used as initial conditions.
 149   //  result:      target for the unfolded result
 150   //  check:       depends on the unfolding method, see comments in specific functions
 151
 152   if (fgMaxInput == -1)
 153   {
 154     Printf("AliUnfolding::Unfold: WARNING. Number of measured bins not set with SetNbins. Using number of bins in measured distribution");
 155     fgMaxInput = measured->GetNbinsX();
 156   }
 157   if (fgMaxParams == -1)
 158   {
 159     Printf("AliUnfolding::Unfold: WARNING. Number of unfolded bins not set with SetNbins. Using number of bins in measured distribution");
 160     fgMaxParams = measured->GetNbinsX();
 161   }
 162
 163   if (fgOverflowBinLimit > 0)
 164     CreateOverflowBin(correlation, measured);
 165
 166   switch (fgMethodType)
 167   {
 168     case kInvalid:
 169     {
 170       Printf("AliUnfolding::Unfold: ERROR: Unfolding method not set. Use SetUnfoldingMethod. Exiting...");
 171       return -1;
 172     }
 173     case kChi2Minimization:
 174       return UnfoldWithMinuit(correlation, efficiency, measured, initialConditions, result, check);
 175     case kBayesian:
 176       return UnfoldWithBayesian(correlation, efficiency, measured, initialConditions, result);
 177     case kFunction:
 178       return UnfoldWithFunction(correlation, efficiency, measured, initialConditions, result);
 179   }
 180   return -1;
 181 }
 182
 183 //____________________________________________________________________
 184 void AliUnfolding::SetStaticVariables(TH2* correlation, TH1* measured)
 185 {
 186   // fill static variables needed for minuit fit
 187
 188   if (!fgCorrelationMatrix)
 189     fgCorrelationMatrix = new TMatrixD(fgMaxInput, fgMaxParams);
 190   if (!fgCorrelationCovarianceMatrix)
 191     fgCorrelationCovarianceMatrix = new TMatrixD(fgMaxInput, fgMaxInput);
 192   if (!fgCurrentESDVector)
 193     fgCurrentESDVector = new TVectorD(fgMaxInput);
 194   if (!fgEntropyAPriori)
 195     fgEntropyAPriori = new TVectorD(fgMaxParams);
 196
 197   fgCorrelationMatrix->Zero();
 198   fgCorrelationCovarianceMatrix->Zero();
 199   fgCurrentESDVector->Zero();
 200   fgEntropyAPriori->Zero();
 201
 202   // normalize correction for given nPart
 203   for (Int_t i=1; i<=correlation->GetNbinsX(); ++i)
 204   {
 205     Double_t sum = correlation->Integral(i, i, 1, correlation->GetNbinsY());
 206     if (sum <= 0)
 207       continue;
 208     Float_t maxValue = 0;
 209     Int_t maxBin = -1;
 210     for (Int_t j=1; j<=correlation->GetNbinsY(); ++j)
 211     {
 212       // find most probably value
 213       if (maxValue < correlation->GetBinContent(i, j))
 214       {
 215         maxValue = correlation->GetBinContent(i, j);
 216         maxBin = j;
 217       }
 218
 219       // npart sum to 1
 220       correlation->SetBinContent(i, j, correlation->GetBinContent(i, j) / sum);
 221       correlation->SetBinError(i, j, correlation->GetBinError(i, j) / sum);
 222
 223       if (i <= fgMaxParams && j <= fgMaxInput)
 224         (*fgCorrelationMatrix)(j-1, i-1) = correlation->GetBinContent(i, j);
 225     }
 226
 227     //printf("MPV for Ntrue = %f is %f\n", fCurrentCorrelation->GetXaxis()->GetBinCenter(i), fCurrentCorrelation->GetYaxis()->GetBinCenter(maxBin));
 228   }
 229
 230   //normalize measured
 231   if (fgNormalizeInput)
 232     measured->Scale(1.0 / measured->Integral());
 233
 234   for (Int_t i=0; i<fgMaxInput; ++i)
 235   {
 236     (*fgCurrentESDVector)[i] = measured->GetBinContent(i+1);
 237     if (measured->GetBinError(i+1) > 0)
 238       (*fgCorrelationCovarianceMatrix)(i, i) = (Double_t) 1e-6 / measured->GetBinError(i+1) / measured->GetBinError(i+1);
 239
 240     if ((*fgCorrelationCovarianceMatrix)(i, i) > 1e7)
 241       (*fgCorrelationCovarianceMatrix)(i, i) = 0;
 242     //Printf("%d, %e", i, (*fgCorrelationCovarianceMatrix)(i, i));
 243   }
 244 }
 245
 246 //____________________________________________________________________
 247 Int_t AliUnfolding::UnfoldWithMinuit(TH2* correlation, TH1* efficiency, TH1* measured, TH1* initialConditions, TH1* result, Bool_t check)
 248 {
 249   //
 250   // implementation of unfolding (internal function)
 251   //
 252   // unfolds <measured> using response from <correlation> and effiency <efficiency>
 253   // output is in <result>
 254   // <initialConditions> set the initial values for the minimization, if 0 <measured> is used
 255   // if <check> is true no unfolding is made, instead only the chi2 without unfolding is printed
 256   //
 257   // returns minuit status (0 = success), (-1 when check was set)
 258   //
 259
 260   SetStaticVariables(correlation, measured);
 261
 262   // Initialize TMinuit via generic fitter interface
 263   TVirtualFitter *minuit = TVirtualFitter::Fitter(0, fgMaxParams);
 264   Double_t arglist[100];
 265
 266   // disable any output (-1), unfortuantly we do not see warnings anymore then. Have to find another way...
 267   arglist[0] = 0;
 268   minuit->ExecuteCommand("SET PRINT", arglist, 1);
 269
 270   // however, enable warnings
 271   //minuit->ExecuteCommand("SET WAR", arglist, 0);
 272
 273   // set minimization function
 274   minuit->SetFCN(Chi2Function);
 275
 276   for (Int_t i=0; i<fgMaxParams; i++)
 277     (*fgEntropyAPriori)[i] = 1;
 278
 279   if (!initialConditions) {
 280     initialConditions = measured;
 281   } else {
 282     Printf("AliUnfolding::UnfoldWithMinuit: Using different initial conditions...");
 283     //new TCanvas; initialConditions->DrawCopy();
 284     if (fgNormalizeInput)
 285       initialConditions->Scale(1.0 / initialConditions->Integral());
 286   }
 287
 288   // set initial conditions as a-priori distribution for MRX regularization
 289   for (Int_t i=0; i<fgMaxParams; i++)
 290     if (initialConditions->GetBinContent(i+1) > 0)
 291       (*fgEntropyAPriori)[i] = initialConditions->GetBinContent(i+1);
 292
 293   Double_t* results = new Double_t[fgMaxParams+1];
 294   for (Int_t i=0; i<fgMaxParams; ++i)
 295   {
 296     results[i] = initialConditions->GetBinContent(i+1);
 297
 298     // minuit sees squared values to prevent it from going negative...
 299     results[i] = TMath::Sqrt(results[i]);
 300
 301     minuit->SetParameter(i, Form("param%d", i), results[i], fgMinuitStepSize, 0, 0);
 302   }
 303
 304   Int_t dummy = 0;
 305   Double_t chi2 = 0;
 306   Chi2Function(dummy, 0, chi2, results, 0);
 307   printf("AliUnfolding::UnfoldWithMinuit: Chi2 of initial parameters is = %f\n", chi2);
 308
 309   if (check)
 310   {
 311     DrawGuess(results);
 312     delete[] results;
 313     return -1;
 314   }
 315
 316   // first param is number of iterations, second is precision....
 317   arglist[0] = 1e6;
 318   //arglist[1] = 1e-5;
 319   //minuit->ExecuteCommand("SCAN", arglist, 0);
 320   Int_t status = minuit->ExecuteCommand("MIGRAD", arglist, 1);
 321   Printf("AliUnfolding::UnfoldWithMinuit: MINUIT status is %d", status);
 322   //printf("!!!!!!!!!!!!!! MIGRAD finished: Starting MINOS !!!!!!!!!!!!!!");
 323   //minuit->ExecuteCommand("MINOS", arglist, 0);
 324
 325   for (Int_t i=0; i<fgMaxParams; ++i)
 326   {
 327     results[i] = minuit->GetParameter(i);
 328     Double_t value = results[i] * results[i];
 329     // error is : (relError) * (value) = (minuit->GetParError(i) / minuit->GetParameter(i)) * (minuit->GetParameter(i) * minuit->GetParameter(i))
 330     Double_t error = minuit->GetParError(i) * results[i];
 331
 332     if (efficiency)
 333     {
 334       if (efficiency->GetBinContent(i+1) > 0)
 335       {
 336         value /= efficiency->GetBinContent(i+1);
 337         error /= efficiency->GetBinContent(i+1);
 338       }
 339       else
 340       {
 341         value = 0;
 342         error = 0;
 343       }
 344     }
 345
 346     result->SetBinContent(i+1, value);
 347     result->SetBinError(i+1, error);
 348   }
 349
 350   if (fgDebug)
 351     DrawGuess(results);
 352
 353   delete[] results;
 354
 355   return status;
 356 }
 357
 358 //____________________________________________________________________
 359 Int_t AliUnfolding::UnfoldWithBayesian(TH2* correlation, TH1* aEfficiency, TH1* measured, TH1* initialConditions, TH1* aResult)
 360 {
 361   //
 362   // unfolds a spectrum using the Bayesian method
 363   //
 364
 365   if (measured->Integral() <= 0)
 366   {
 367     Printf("AliUnfolding::UnfoldWithBayesian: ERROR: The measured spectrum is empty");
 368     return -1;
 369   }
 370
 371   const Int_t kStartBin = 0;
 372
 373   Int_t kMaxM = fgMaxInput;  //<= fCurrentCorrelation->GetNbinsY(); // max measured axis
 374   Int_t kMaxT = fgMaxParams; //<= fCurrentCorrelation->GetNbinsX(); // max true axis
 375
 376   // convergence limit: kMaxT * 0.001^2 = kMaxT * 1e-6 (e.g. 250 bins --> 2.5 e-4)
 377   const Double_t kConvergenceLimit = kMaxT * 1e-6;
 378
 379   // store information in arrays, to increase processing speed (~ factor 5)
 380   Double_t* measuredCopy = new Double_t[kMaxM];
 381   Double_t* measuredError = new Double_t[kMaxM];
 382   Double_t* prior = new Double_t[kMaxT];
 383   Double_t* result = new Double_t[kMaxT];
 384   Double_t* efficiency = new Double_t[kMaxT];
 385
 386   Double_t** response = new Double_t*[kMaxT];
 387   Double_t** inverseResponse = new Double_t*[kMaxT];
 388   for (Int_t i=0; i<kMaxT; i++)
 389   {
 390     response[i] = new Double_t[kMaxM];
 391     inverseResponse[i] = new Double_t[kMaxM];
 392   }
 393
 394   // for normalization
 395   Float_t measuredIntegral = measured->Integral();
 396   for (Int_t m=0; m<kMaxM; m++)
 397   {
 398     measuredCopy[m] = measured->GetBinContent(m+1) / measuredIntegral;
 399     measuredError[m] = measured->GetBinError(m+1) / measuredIntegral;
 400
 401     for (Int_t t=0; t<kMaxT; t++)
 402     {
 403       response[t][m] = correlation->GetBinContent(t+1, m+1);
 404       inverseResponse[t][m] = 0;
 405     }
 406   }
 407
 408   for (Int_t t=0; t<kMaxT; t++)
 409   {
 410     if (aEfficiency)
 411     {
 412       efficiency[t] = aEfficiency->GetBinContent(t+1);
 413     }
 414     else
 415       efficiency[t] = 1;
 416
 417     prior[t] = measuredCopy[t];
 418     result[t] = 0;
 419   }
 420
 421   // pick prior distribution
 422   if (initialConditions)
 423   {
 424     printf("Using different starting conditions...\n");
 425     // for normalization
 426     Float_t inputDistIntegral = initialConditions->Integral();
 427     for (Int_t i=0; i<kMaxT; i++)
 428       prior[i] = initialConditions->GetBinContent(i+1) / inputDistIntegral;
 429   }
 430
 431   //TH1F* convergence = new TH1F("convergence", "convergence", 200, 0.5, 200.5);
 432
 433   // unfold...
 434   for (Int_t i=0; i<fgBayesianIterations || fgBayesianIterations< 0; i++)
 435   {
 436     if (fgDebug)
 437       Printf("AliUnfolding::UnfoldWithBayesian: iteration %i", i);
 438
 439     // calculate IR from Beyes theorem
 440     // IR_ji = R_ij * prior_i / sum_k(R_kj * prior_k)
 441
 442     Double_t chi2Measured = 0;
 443     for (Int_t m=0; m<kMaxM; m++)
 444     {
 445       Float_t norm = 0;
 446       for (Int_t t = kStartBin; t<kMaxT; t++)
 447         norm += response[t][m] * prior[t];
 448
 449       // calc. chi2: (measured - response * prior) / error
 450       if (measuredError[m] > 0)
 451       {
 452         Double_t value = (measuredCopy[m] - norm) / measuredError[m];
 453         chi2Measured += value * value;
 454       }
 455
 456       if (norm > 0)
 457       {
 458         for (Int_t t = kStartBin; t<kMaxT; t++)
 459           inverseResponse[t][m] = response[t][m] * prior[t] / norm;
 460       }
 461       else
 462       {
 463         for (Int_t t = kStartBin; t<kMaxT; t++)
 464           inverseResponse[t][m] = 0;
 465       }
 466     }
 467     //Printf("chi2Measured of the last prior is %e", chi2Measured);
 468
 469     for (Int_t t = kStartBin; t<kMaxT; t++)
 470     {
 471       Float_t value = 0;
 472       for (Int_t m=0; m<kMaxM; m++)
 473         value += inverseResponse[t][m] * measuredCopy[m];
 474
 475       if (efficiency[t] > 0)
 476         result[t] = value / efficiency[t];
 477       else
 478         result[t] = 0;
 479     }
 480
 481     // draw intermediate result
 482     /*
 483     for (Int_t t=0; t<kMaxT; t++)
 484       aResult->SetBinContent(t+1, result[t]);
 485     aResult->SetMarkerStyle(20+i);
 486     aResult->SetMarkerColor(2);
 487     aResult->DrawCopy("P SAME HIST");
 488     */
 489
 490     Double_t chi2LastIter = 0;
 491     // regularization (simple smoothing)
 492     for (Int_t t=kStartBin; t<kMaxT; t++)
 493     {
 494       Float_t newValue = 0;
 495
 496       // 0 bin excluded from smoothing
 497       if (t > kStartBin+2 && t<kMaxT-1)
 498       {
 499         Float_t average = (result[t-1] + result[t] + result[t+1]) / 3;
 500
 501         // weight the average with the regularization parameter
 502         newValue = (1 - fgBayesianSmoothing) * result[t] + fgBayesianSmoothing * average;
 503       }
 504       else
 505         newValue = result[t];
 506
 507       // calculate chi2 (change from last iteration)
 508       if (prior[t] > 1e-5)
 509       {
 510         Double_t diff = (prior[t] - newValue) / prior[t];
 511         chi2LastIter += diff * diff;
 512       }
 513
 514       prior[t] = newValue;
 515     }
 516     //printf("Chi2 of %d iteration = %e\n", i, chi2LastIter);
 517     //convergence->Fill(i+1, chi2LastIter);
 518
 519     if (fgBayesianIterations < 0 && chi2LastIter < kConvergenceLimit)
 520     {
 521       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);
 522       break;
 523     }
 524   } // end of iterations
 525
 526   //new TCanvas; convergence->DrawCopy(); gPad->SetLogy();
 527   //delete convergence;
 528
 529   for (Int_t t=0; t<kMaxT; t++)
 530     aResult->SetBinContent(t+1, result[t]);
 531
 532   delete[] measuredCopy;
 533   delete[] measuredError;
 534   delete[] prior;
 535   delete[] result;
 536   delete[] efficiency;
 537
 538   for (Int_t i=0; i<kMaxT; i++)
 539   {
 540     delete[] response[i];
 541     delete[] inverseResponse[i];
 542   }
 543   delete[] response;
 544   delete[] inverseResponse;
 545
 546   return 0;
 547
 548   // ********
 549   // Calculate the covariance matrix, all arguments are taken from NIM,A362,487-498,1995
 550
 551   /*printf("Calculating covariance matrix. This may take some time...\n");
 552
 553   // check if this is the right one...
 554   TH1* sumHist = GetMultiplicityMC(inputRange, eventType)->ProjectionY("sumHist", 1, GetMultiplicityMC(inputRange, eventType)->GetNbinsX());
 555
 556   Int_t xBins = hInverseResponseBayes->GetNbinsX();
 557   Int_t yBins = hInverseResponseBayes->GetNbinsY();
 558
 559   // calculate "unfolding matrix" Mij
 560   Float_t matrixM[251][251];
 561   for (Int_t i=1; i<=xBins; i++)
 562   {
 563     for (Int_t j=1; j<=yBins; j++)
 564     {
 565       if (fCurrentEfficiency->GetBinContent(i) > 0)
 566         matrixM[i-1][j-1] = hInverseResponseBayes->GetBinContent(i, j) / fCurrentEfficiency->GetBinContent(i);
 567       else
 568         matrixM[i-1][j-1] = 0;
 569     }
 570   }
 571
 572   Float_t* vectorn = new Float_t[yBins];
 573   for (Int_t j=1; j<=yBins; j++)
 574     vectorn[j-1] = fCurrentESD->GetBinContent(j);
 575
 576   // first part of covariance matrix, depends on input distribution n(E)
 577   Float_t cov1[251][251];
 578
 579   Float_t nEvents = fCurrentESD->Integral(); // N
 580
 581   xBins = 20;
 582   yBins = 20;
 583
 584   for (Int_t k=0; k<xBins; k++)
 585   {
 586     printf("In Cov1: %d\n", k);
 587     for (Int_t l=0; l<yBins; l++)
 588     {
 589       cov1[k][l] = 0;
 590
 591       // sum_j Mkj Mlj n(Ej) * (1 - n(Ej) / N)
 592       for (Int_t j=0; j<yBins; j++)
 593         cov1[k][l] += matrixM[k][j] * matrixM[l][j] * vectorn[j]
 594           * (1.0 - vectorn[j] / nEvents);
 595
 596       // - sum_i,j (i != j) Mki Mlj n(Ei) n(Ej) / N
 597       for (Int_t i=0; i<yBins; i++)
 598         for (Int_t j=0; j<yBins; j++)
 599         {
 600           if (i == j)
 601             continue;
 602           cov1[k][l] -= matrixM[k][i] * matrixM[l][j] * vectorn[i]
 603             * vectorn[j] / nEvents;
 604          }
 605     }
 606   }
 607
 608   printf("Cov1 finished\n");
 609
 610   TH2F* cov = (TH2F*) hInverseResponseBayes->Clone("cov");
 611   cov->Reset();
 612
 613   for (Int_t i=1; i<=xBins; i++)
 614     for (Int_t j=1; j<=yBins; j++)
 615       cov->SetBinContent(i, j, cov1[i-1][j-1]);
 616
 617   new TCanvas;
 618   cov->Draw("COLZ");
 619
 620   // second part of covariance matrix, depends on response matrix
 621   Float_t cov2[251][251];
 622
 623   // Cov[P(Er|Cu), P(Es|Cu)] term
 624   Float_t covTerm[100][100][100];
 625   for (Int_t r=0; r<yBins; r++)
 626     for (Int_t u=0; u<xBins; u++)
 627       for (Int_t s=0; s<yBins; s++)
 628       {
 629         if (r == s)
 630           covTerm[r][u][s] = 1.0 / sumHist->GetBinContent(u+1) * hResponse->GetBinContent(u+1, r+1)
 631             * (1.0 - hResponse->GetBinContent(u+1, r+1));
 632         else
 633           covTerm[r][u][s] = - 1.0 / sumHist->GetBinContent(u+1) * hResponse->GetBinContent(u+1, r+1)
 634             * hResponse->GetBinContent(u+1, s+1);
 635       }
 636
 637   for (Int_t k=0; k<xBins; k++)
 638     for (Int_t l=0; l<yBins; l++)
 639     {
 640       cov2[k][l] = 0;
 641       printf("In Cov2: %d %d\n", k, l);
 642       for (Int_t i=0; i<yBins; i++)
 643         for (Int_t j=0; j<yBins; j++)
 644         {
 645           //printf("In Cov2: %d %d %d %d\n", k, l, i, j);
 646           // calculate Cov(Mki, Mlj) = sum{ru},{su} ...
 647           Float_t tmpCov = 0;
 648           for (Int_t r=0; r<yBins; r++)
 649             for (Int_t u=0; u<xBins; u++)
 650               for (Int_t s=0; s<yBins; s++)
 651               {
 652                 if (hResponse->GetBinContent(u+1, r+1) == 0 || hResponse->GetBinContent(u+1, s+1) == 0
 653                   || hResponse->GetBinContent(u+1, i+1) == 0)
 654                   continue;
 655
 656                 tmpCov += BayesCovarianceDerivate(matrixM, hResponse, fCurrentEfficiency, k, i, r, u)
 657                   * BayesCovarianceDerivate(matrixM, hResponse, fCurrentEfficiency, l, j, s, u)
 658                   * covTerm[r][u][s];
 659               }
 660
 661           cov2[k][l] += fCurrentESD->GetBinContent(i+1) * fCurrentESD->GetBinContent(j+1) * tmpCov;
 662         }
 663     }
 664
 665   printf("Cov2 finished\n");
 666
 667   for (Int_t i=1; i<=xBins; i++)
 668     for (Int_t j=1; j<=yBins; j++)
 669       cov->SetBinContent(i, j, cov1[i-1][j-1] + cov2[i-1][j-1]);
 670
 671   new TCanvas;
 672   cov->Draw("COLZ");*/
 673 }
 674
 675 //____________________________________________________________________
 676 Double_t AliUnfolding::RegularizationPol0(TVectorD& params)
 677 {
 678   // homogenity term for minuit fitting
 679   // pure function of the parameters
 680   // prefers constant function (pol0)
 681
 682   Double_t chi2 = 0;
 683
 684   for (Int_t i=1+fgSkipBinsBegin; i<fgMaxParams; ++i)
 685   {
 686     Double_t right  = params[i];
 687     Double_t left   = params[i-1];
 688
 689     if (right != 0)
 690     {
 691       Double_t diff = 1 - left / right;
 692       chi2 += diff * diff;
 693     }
 694   }
 695
 696   return chi2 / 100.0;
 697 }
 698
 699 //____________________________________________________________________
 700 Double_t AliUnfolding::RegularizationPol1(TVectorD& params)
 701 {
 702   // homogenity term for minuit fitting
 703   // pure function of the parameters
 704   // prefers linear function (pol1)
 705
 706   Double_t chi2 = 0;
 707
 708   for (Int_t i=2+fgSkipBinsBegin; i<fgMaxParams; ++i)
 709   {
 710     if (params[i-1] == 0)
 711       continue;
 712
 713     Double_t right  = params[i];
 714     Double_t middle = params[i-1];
 715     Double_t left   = params[i-2];
 716
 717     Double_t der1 = (right - middle);
 718     Double_t der2 = (middle - left);
 719
 720     Double_t diff = (der1 - der2) / middle;
 721
 722     chi2 += diff * diff;
 723   }
 724
 725   return chi2;
 726 }
 727
 728 //____________________________________________________________________
 729 Double_t AliUnfolding::RegularizationLog(TVectorD& params)
 730 {
 731   // homogenity term for minuit fitting
 732   // pure function of the parameters
 733   // prefers linear function (pol1)
 734
 735   Double_t chi2 = 0;
 736
 737   for (Int_t i=2+fgSkipBinsBegin; i<fgMaxParams; ++i)
 738   {
 739     if (params[i-1] == 0)
 740       continue;
 741
 742     Double_t right  = log(params[i]);
 743     Double_t middle = log(params[i-1]);
 744     Double_t left   = log(params[i-2]);
 745
 746     Double_t der1 = (right - middle);
 747     Double_t der2 = (middle - left);
 748
 749     Double_t diff = (der1 - der2) / middle;
 750
 751     chi2 += diff * diff;
 752   }
 753
 754   return chi2 * 100;
 755 }
 756
 757 //____________________________________________________________________
 758 Double_t AliUnfolding::RegularizationTotalCurvature(TVectorD& params)
 759 {
 760   // homogenity term for minuit fitting
 761   // pure function of the parameters
 762   // minimizes the total curvature (from Unfolding Methods In High-Energy Physics Experiments,
 763   // V. Blobel (Hamburg U.) . DESY 84/118, Dec 1984. 40pp.
 764
 765   Double_t chi2 = 0;
 766
 767   for (Int_t i=2+fgSkipBinsBegin; i<fgMaxParams; ++i)
 768   {
 769     Double_t right  = params[i];
 770     Double_t middle = params[i-1];
 771     Double_t left   = params[i-2];
 772
 773     Double_t der1 = (right - middle);
 774     Double_t der2 = (middle - left);
 775
 776     Double_t diff = (der1 - der2);
 777
 778     chi2 += diff * diff;
 779   }
 780
 781   return chi2 * 1e4;
 782 }
 783
 784 //____________________________________________________________________
 785 Double_t AliUnfolding::RegularizationEntropy(TVectorD& params)
 786 {
 787   // homogenity term for minuit fitting
 788   // pure function of the parameters
 789   // calculates entropy, from
 790   // The method of reduced cross-entropy (M. Schmelling 1993)
 791
 792   Double_t paramSum = 0;
 793
 794   for (Int_t i=fgSkipBinsBegin; i<fgMaxParams; ++i)
 795     paramSum += params[i];
 796
 797   Double_t chi2 = 0;
 798   for (Int_t i=fgSkipBinsBegin; i<fgMaxParams; ++i)
 799   {
 800     //Double_t tmp = params[i] / paramSum;
 801     Double_t tmp = params[i];
 802     if (tmp > 0 && (*fgEntropyAPriori)[i] > 0)
 803     {
 804       chi2 += tmp * TMath::Log(tmp / (*fgEntropyAPriori)[i]);
 805     }
 806   }
 807
 808   return 100.0 + chi2;
 809 }
 810
 811 //____________________________________________________________________
 812 void AliUnfolding::Chi2Function(Int_t&, Double_t*, Double_t& chi2, Double_t *params, Int_t)
 813 {
 814   //
 815   // fit function for minuit
 816   // does: (m - Ad)W(m - Ad) where m = measured, A correlation matrix, d = guess, W = covariance matrix
 817   //
 818
 819   // d
 820   static TVectorD paramsVector(fgMaxParams);
 821   for (Int_t i=0; i<fgMaxParams; ++i)
 822     paramsVector[i] = params[i] * params[i];
 823
 824   // calculate penalty factor
 825   Double_t penaltyVal = 0;
 826   switch (fgRegularizationType)
 827   {
 828     case kNone:       break;
 829     case kPol0:       penaltyVal = RegularizationPol0(paramsVector); break;
 830     case kPol1:       penaltyVal = RegularizationPol1(paramsVector); break;
 831     case kCurvature:  penaltyVal = RegularizationTotalCurvature(paramsVector); break;
 832     case kEntropy:    penaltyVal = RegularizationEntropy(paramsVector); break;
 833     case kLog:        penaltyVal = RegularizationLog(paramsVector); break;
 834   }
 835
 836   //if (penaltyVal > 1e10)
 837   //  paramsVector2.Print();
 838
 839   penaltyVal *= fgRegularizationWeight;
 840
 841   // Ad
 842   TVectorD measGuessVector(fgMaxInput);
 843   measGuessVector = (*fgCorrelationMatrix) * paramsVector;
 844
 845   // Ad - m
 846   measGuessVector -= (*fgCurrentESDVector);
 847
 848   //measGuessVector.Print();
 849
 850   TVectorD copy(measGuessVector);
 851
 852   // (Ad - m) W
 853   // this step can be optimized because currently only the diagonal elements of fgCorrelationCovarianceMatrix are used
 854   // normal way is like this:
 855   // measGuessVector *= (*fgCorrelationCovarianceMatrix);
 856   // optimized way like this:
 857   for (Int_t i=0; i<fgMaxInput; ++i)
 858     measGuessVector[i] *= (*fgCorrelationCovarianceMatrix)(i, i);
 859
 860   //measGuessVector.Print();
 861
 862   // (Ad - m) W (Ad - m)
 863   // the factor 1e6 prevents that a very small number (measGuessVector[i]) is multiplied with a very
 864   // big number ((*fgCorrelationCovarianceMatrix)(i, i)) (see ApplyMinuitFit)
 865   Double_t chi2FromFit = measGuessVector * copy * 1e6;
 866
 867   chi2 = chi2FromFit + penaltyVal;
 868
 869   static Int_t callCount = 0;
 870   if ((callCount++ % 10000) == 0)
 871     Printf("AliUnfolding::Chi2Function: Iteration %d: %f %f --> %f", callCount, chi2FromFit, penaltyVal, chi2);
 872 }
 873
 874 //____________________________________________________________________
 875 void AliUnfolding::DrawGuess(Double_t *params)
 876 {
 877   //
 878   // draws residuals of solution suggested by params and effect of regularization
 879   //
 880
 881   // d
 882   static TVectorD paramsVector(fgMaxParams);
 883   for (Int_t i=0; i<fgMaxParams; ++i)
 884     paramsVector[i] = params[i] * params[i];
 885
 886   // Ad
 887   TVectorD measGuessVector(fgMaxInput);
 888   measGuessVector = (*fgCorrelationMatrix) * paramsVector;
 889
 890   // Ad - m
 891   measGuessVector -= (*fgCurrentESDVector);
 892
 893   TVectorD copy(measGuessVector);
 894   //copy.Print();
 895
 896   // (Ad - m) W
 897   // this step can be optimized because currently only the diagonal elements of fgCorrelationCovarianceMatrix are used
 898   // normal way is like this:
 899   // measGuessVector *= (*fgCorrelationCovarianceMatrix);
 900   // optimized way like this:
 901   for (Int_t i=0; i<fgMaxInput; ++i)
 902     measGuessVector[i] *= (*fgCorrelationCovarianceMatrix)(i, i);
 903
 904   // (Ad - m) W (Ad - m)
 905   // the factor 1e6 prevents that a very small number (measGuessVector[i]) is multiplied with a very
 906   // big number ((*fgCorrelationCovarianceMatrix)(i, i)) (see ApplyMinuitFit)
 907   //Double_t chi2FromFit = measGuessVector * copy * 1e6;
 908
 909   // draw residuals
 910   TH1* residuals = new TH1F("residuals", "residuals", fgMaxInput, -0.5, fgMaxInput - 0.5);
 911   for (Int_t i=0; i<fgMaxInput; ++i)
 912     residuals->SetBinContent(i+1, measGuessVector(i) * copy(i) * 1e6);
 913   new TCanvas; residuals->DrawCopy(); gPad->SetLogy();
 914
 915   // draw penalty
 916   if (fgRegularizationType != kPol1) {
 917     Printf("Drawing not supported");
 918     return;
 919   }
 920
 921   TH1* penalty = new TH1F("penalty", "penalty", fgMaxParams, -0.5, fgMaxParams - 0.5);
 922
 923   for (Int_t i=2+1; i<fgMaxParams; ++i)
 924   {
 925     if (params[i-1] == 0)
 926       continue;
 927
 928     Double_t right  = params[i];
 929     Double_t middle = params[i-1];
 930     Double_t left   = params[i-2];
 931
 932     Double_t der1 = (right - middle);
 933     Double_t der2 = (middle - left);
 934
 935     Double_t diff = (der1 - der2) / middle;
 936
 937     penalty->SetBinContent(i-1, diff * diff);
 938
 939     //Printf("%d %f %f %f %f", i-1, left, middle, right, diff);
 940   }
 941   new TCanvas; penalty->DrawCopy(); gPad->SetLogy();
 942 }
 943
 944 //____________________________________________________________________
 945 void AliUnfolding::TF1Function(Int_t& unused1, Double_t* unused2, Double_t& chi2, Double_t *params, Int_t unused3)
 946 {
 947   //
 948   // fit function for minuit
 949   // uses the TF1 stored in fgFitFunction
 950   //
 951
 952   for (Int_t i=0; i<fgFitFunction->GetNpar(); i++)
 953     fgFitFunction->SetParameter(i, params[i]);
 954
 955   Double_t* params2 = new Double_t[fgMaxParams];
 956
 957   for (Int_t i=0; i<fgMaxParams; ++i)
 958     params2[i] = fgFitFunction->Eval(i);
 959
 960   Chi2Function(unused1, unused2, chi2, params2, unused3);
 961
 962   delete[] params2;
 963
 964   if (fgDebug)
 965     Printf("%f", chi2);
 966 }
 967
 968 //____________________________________________________________________
 969 Int_t AliUnfolding::UnfoldWithFunction(TH2* correlation, TH1* efficiency, TH1* measured, TH1* /* initialConditions */, TH1* aResult)
 970 {
 971   //
 972   // correct spectrum using minuit chi2 method applying a functional fit
 973   //
 974
 975   if (!fgFitFunction)
 976   {
 977     Printf("AliUnfolding::UnfoldWithFunction: ERROR fit function not set. Exiting.");
 978     return -1;
 979   }
 980
 981   SetChi2Regularization(kNone, 0);
 982
 983   SetStaticVariables(correlation, measured);
 984
 985   // Initialize TMinuit via generic fitter interface
 986   TVirtualFitter *minuit = TVirtualFitter::Fitter(0, fgFitFunction->GetNpar());
 987
 988   minuit->SetFCN(TF1Function);
 989   for (Int_t i=0; i<fgFitFunction->GetNpar(); i++)
 990   {
 991     Double_t lower, upper;
 992     fgFitFunction->GetParLimits(i, lower, upper);
 993     minuit->SetParameter(i, Form("param%d",i), fgFitFunction->GetParameter(i), fgMinuitStepSize, lower, upper);
 994   }
 995
 996   Double_t arglist[100];
 997   arglist[0] = 0;
 998   minuit->ExecuteCommand("SET PRINT", arglist, 1);
 999   minuit->ExecuteCommand("MIGRAD", arglist, 0);
1000   //minuit->ExecuteCommand("MINOS", arglist, 0);
1001
1002   for (Int_t i=0; i<fgFitFunction->GetNpar(); i++)
1003     fgFitFunction->SetParameter(i, minuit->GetParameter(i));
1004
1005   for (Int_t i=0; i<fgMaxParams; ++i)
1006   {
1007     Double_t value = fgFitFunction->Eval(i);
1008     if (fgDebug)
1009       Printf("%d : %f", i, value);
1010
1011     if (efficiency)
1012     {
1013       if (efficiency->GetBinContent(i+1) > 0)
1014       {
1015         value /= efficiency->GetBinContent(i+1);
1016       }
1017       else
1018         value = 0;
1019     }
1020     aResult->SetBinContent(i+1, value);
1021     aResult->SetBinError(i+1, 0);
1022   }
1023
1024   return 0;
1025 }
1026
1027 //____________________________________________________________________
1028 void AliUnfolding::CreateOverflowBin(TH2* correlation, TH1* measured)
1029 {
1030   // Finds the first bin where the content is below fgStatLimit and combines all values for this bin and larger bins
1031   // The same limit is applied to the correlation
1032
1033   Int_t lastBin = 0;
1034   for (Int_t i=1; i<=measured->GetNbinsX(); ++i)
1035   {
1036     if (measured->GetBinContent(i) <= fgOverflowBinLimit)
1037     {
1038       lastBin = i;
1039       break;
1040     }
1041   }
1042
1043   if (lastBin == 0)
1044   {
1045     Printf("AliUnfolding::CreateOverflowBin: WARNING: First bin is already below limit of %f", fgOverflowBinLimit);
1046     return;
1047   }
1048
1049   Printf("AliUnfolding::CreateOverflowBin: Bin limit in measured spectrum determined to be %d", lastBin);
1050
1051   TCanvas* canvas = 0;
1052
1053   if (fgDebug)
1054   {
1055     canvas = new TCanvas("StatSolution", "StatSolution", 1000, 800);
1056     canvas->Divide(2, 2);
1057
1058     canvas->cd(1);
1059     measured->SetStats(kFALSE);
1060     measured->DrawCopy();
1061     gPad->SetLogy();
1062
1063     canvas->cd(2);
1064     correlation->SetStats(kFALSE);
1065     correlation->DrawCopy("COLZ");
1066   }
1067
1068   measured->SetBinContent(lastBin, measured->Integral(lastBin, measured->GetNbinsX()));
1069   for (Int_t i=lastBin+1; i<=measured->GetNbinsX(); ++i)
1070   {
1071     measured->SetBinContent(i, 0);
1072     measured->SetBinError(i, 0);
1073   }
1074   // the error is set to sqrt(N), better solution possible?, sum of relative errors of all contributions???
1075   measured->SetBinError(lastBin, TMath::Sqrt(measured->GetBinContent(lastBin)));
1076
1077   Printf("AliUnfolding::CreateOverflowBin: This bin has now %f +- %f entries", measured->GetBinContent(lastBin), measured->GetBinError(lastBin));
1078
1079   for (Int_t i=1; i<=correlation->GetNbinsX(); ++i)
1080   {
1081     correlation->SetBinContent(i, lastBin, correlation->Integral(i, i, lastBin, correlation->GetNbinsY()));
1082     // the error is set to sqrt(N), better solution possible?, sum of relative errors of all contributions???
1083     correlation->SetBinError(i, lastBin, TMath::Sqrt(correlation->GetBinContent(i, lastBin)));
1084
1085     for (Int_t j=lastBin+1; j<=correlation->GetNbinsY(); ++j)
1086     {
1087       correlation->SetBinContent(i, j, 0);
1088       correlation->SetBinError(i, j, 0);
1089     }
1090   }
1091
1092   Printf("AliUnfolding::CreateOverflowBin: Adjusted correlation matrix!");
1093
1094   if (canvas)
1095   {
1096     canvas->cd(3);
1097     measured->DrawCopy();
1098     gPad->SetLogy();
1099
1100     canvas->cd(4);
1101     correlation->DrawCopy("COLZ");
1102   }
1103 }