PWG0/multiplicity/AliMultiplicityCorrection.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$ */
  17
  18 // This class is used to store correction maps, raw input and results of the multiplicity
  19 // measurement with the ITS or TPC
  20 // It also contains functions to correct the spectrum using different methods.
  21 // e.g. chi2 minimization and bayesian unfolding
  22 //
  23 //  Author: Jan.Fiete.Grosse-Oetringhaus@cern.ch
  24
  25 #include "AliMultiplicityCorrection.h"
  26
  27 #include <TFile.h>
  28 #include <TH1F.h>
  29 #include <TH2F.h>
  30 #include <TH3F.h>
  31 #include <TDirectory.h>
  32 #include <TVirtualFitter.h>
  33 #include <TCanvas.h>
  34 #include <TString.h>
  35 #include <TF1.h>
  36 #include <TMath.h>
  37 #include <TCollection.h>
  38 #include <TLegend.h>
  39 #include <TLine.h>
  40 #include <TRandom.h>
  41 #include <TProfile.h>
  42 #include <TProfile2D.h>
  43
  44 ClassImp(AliMultiplicityCorrection)
  45
  46 const Int_t AliMultiplicityCorrection::fgkMaxInput  = 250;  // bins in measured histogram
  47 const Int_t AliMultiplicityCorrection::fgkMaxParams = 250;  // bins in unfolded histogram = number of fit params
  48
  49 TMatrixD* AliMultiplicityCorrection::fgCorrelationMatrix = 0;
  50 TMatrixD* AliMultiplicityCorrection::fgCorrelationCovarianceMatrix = 0;
  51 TVectorD* AliMultiplicityCorrection::fgCurrentESDVector = 0;
  52 TVectorD* AliMultiplicityCorrection::fgEntropyAPriori = 0;
  53
  54 Bool_t AliMultiplicityCorrection::fgCreateBigBin = kTRUE;
  55 AliMultiplicityCorrection::RegularizationType AliMultiplicityCorrection::fgRegularizationType = AliMultiplicityCorrection::kPol1;
  56 Float_t AliMultiplicityCorrection::fgRegularizationWeight = 10000;
  57 Int_t AliMultiplicityCorrection::fgNParamsMinuit = AliMultiplicityCorrection::fgkMaxParams;
  58
  59 TF1* AliMultiplicityCorrection::fgNBD = 0;
  60
  61 Float_t AliMultiplicityCorrection::fgBayesianSmoothing  = 1;           // smoothing parameter (0 = no smoothing)
  62 Int_t   AliMultiplicityCorrection::fgBayesianIterations = 100;         // number of iterations in Bayesian method
  63
  64 // These are the areas where the quality of the unfolding results are evaluated
  65 // Default defined here, call SetQualityRegions to change them
  66 // unit is in multiplicity (not in bin!)
  67
  68 // SPD:   peak area - flat area - low stat area
  69 Int_t AliMultiplicityCorrection::fgQualityRegionsB[kQualityRegions] = {4,  60,  180};
  70 Int_t AliMultiplicityCorrection::fgQualityRegionsE[kQualityRegions] = {20, 140, 210};
  71
  72 //____________________________________________________________________
  73 void AliMultiplicityCorrection::SetQualityRegions(Bool_t SPDStudy)
  74 {
  75   //
  76   // sets the quality region definition to TPC or SPD
  77   //
  78
  79   if (SPDStudy)
  80   {
  81     // SPD:   peak area - flat area - low stat area
  82     fgQualityRegionsB[0] = 4;
  83     fgQualityRegionsE[0] = 20;
  84
  85     fgQualityRegionsB[1] = 60;
  86     fgQualityRegionsE[1] = 140;
  87
  88     fgQualityRegionsB[2] = 180;
  89     fgQualityRegionsE[2] = 210;
  90
  91     Printf("AliMultiplicityCorrection::SetQualityRegions --> Enabled quality regions for SPD");
  92   }
  93   else
  94   {
  95     // TPC:   peak area - flat area - low stat area
  96     fgQualityRegionsB[0] = 4;
  97     fgQualityRegionsE[0] = 12;
  98
  99     fgQualityRegionsB[1] = 25;
 100     fgQualityRegionsE[1] = 55;
 101
 102     fgQualityRegionsB[2] = 88;
 103     fgQualityRegionsE[2] = 108;
 104
 105     Printf("AliMultiplicityCorrection::SetQualityRegions --> Enabled quality regions for TPC");
 106   }
 107 }
 108
 109 //____________________________________________________________________
 110 AliMultiplicityCorrection::AliMultiplicityCorrection() :
 111   TNamed(), fCurrentESD(0), fCurrentCorrelation(0), fCurrentEfficiency(0), fLastBinLimit(0), fLastChi2MC(0), fLastChi2MCLimit(0), fLastChi2Residuals(0), fRatioAverage(0)
 112 {
 113   //
 114   // default constructor
 115   //
 116
 117   for (Int_t i = 0; i < kESDHists; ++i)
 118     fMultiplicityESD[i] = 0;
 119
 120   for (Int_t i = 0; i < kMCHists; ++i)
 121   {
 122     fMultiplicityVtx[i] = 0;
 123     fMultiplicityMB[i] = 0;
 124     fMultiplicityINEL[i] = 0;
 125   }
 126
 127   for (Int_t i = 0; i < kCorrHists; ++i)
 128   {
 129     fCorrelation[i] = 0;
 130     fMultiplicityESDCorrected[i] = 0;
 131   }
 132
 133   for (Int_t i = 0; i < kQualityRegions; ++i)
 134     fQuality[i] = 0;
 135 }
 136
 137 //____________________________________________________________________
 138 AliMultiplicityCorrection::AliMultiplicityCorrection(const Char_t* name, const Char_t* title) :
 139   TNamed(name, title),
 140   fCurrentESD(0),
 141   fCurrentCorrelation(0),
 142   fCurrentEfficiency(0),
 143   fLastBinLimit(0),
 144   fLastChi2MC(0),
 145   fLastChi2MCLimit(0),
 146   fLastChi2Residuals(0),
 147   fRatioAverage(0)
 148 {
 149   //
 150   // named constructor
 151   //
 152
 153   // do not add this hists to the directory
 154   Bool_t oldStatus = TH1::AddDirectoryStatus();
 155   TH1::AddDirectory(kFALSE);
 156
 157   /*Float_t binLimitsVtx[] = {-10,-9,-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8,9,10};
 158   Float_t binLimitsN[] = {-0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5,
 159                           10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5,
 160                           20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5,
 161                           30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5,
 162                           40.5, 41.5, 42.5, 43.5, 44.5, 45.5, 46.5, 47.5, 48.5, 49.5,
 163                           50.5, 55.5, 60.5, 65.5, 70.5, 75.5, 80.5, 85.5, 90.5, 95.5,
 164                           100.5, 105.5, 110.5, 115.5, 120.5, 125.5, 130.5, 135.5, 140.5, 145.5,
 165                           150.5, 160.5, 170.5, 180.5, 190.5, 200.5, 210.5, 220.5, 230.5, 240.5,
 166                           250.5, 275.5, 300.5, 325.5, 350.5, 375.5, 400.5, 425.5, 450.5, 475.5,
 167                           500.5 };
 168                                 //525.5, 550.5, 575.5, 600.5, 625.5, 650.5, 675.5, 700.5, 725.5,
 169                           //750.5, 775.5, 800.5, 825.5, 850.5, 875.5, 900.5, 925.5, 950.5, 975.5,
 170                           //1000.5 };
 171
 172   #define VTXBINNING 10, binLimitsVtx
 173   #define NBINNING fgkMaxParams, binLimitsN*/
 174
 175   #define NBINNING 501, -0.5, 500.5
 176   #define VTXBINNING 1, -10, 10
 177
 178   for (Int_t i = 0; i < kESDHists; ++i)
 179     fMultiplicityESD[i] = new TH2F(Form("fMultiplicityESD%d", i), "fMultiplicityESD;vtx-z;Ntracks;Count", VTXBINNING, NBINNING);
 180
 181   for (Int_t i = 0; i < kMCHists; ++i)
 182   {
 183     fMultiplicityVtx[i] = dynamic_cast<TH2F*> (fMultiplicityESD[0]->Clone(Form("fMultiplicityVtx%d", i)));
 184     fMultiplicityVtx[i]->SetTitle("fMultiplicityVtx;vtx-z;Npart");
 185
 186     fMultiplicityMB[i] = dynamic_cast<TH2F*> (fMultiplicityVtx[0]->Clone(Form("fMultiplicityMB%d", i)));
 187     fMultiplicityMB[i]->SetTitle("fMultiplicityMB");
 188
 189     fMultiplicityINEL[i] = dynamic_cast<TH2F*> (fMultiplicityVtx[0]->Clone(Form("fMultiplicityINEL%d", i)));
 190     fMultiplicityINEL[i]->SetTitle("fMultiplicityINEL");
 191   }
 192
 193   for (Int_t i = 0; i < kCorrHists; ++i)
 194   {
 195     fCorrelation[i] = new TH3F(Form("fCorrelation%d", i), "fCorrelation;vtx-z;Npart;Ntracks", VTXBINNING, NBINNING, NBINNING);
 196     fMultiplicityESDCorrected[i] = new TH1F(Form("fMultiplicityESDCorrected%d", i), "fMultiplicityESDCorrected;Npart;dN/dN", NBINNING);
 197   }
 198
 199   TH1::AddDirectory(oldStatus);
 200 }
 201
 202 //____________________________________________________________________
 203 AliMultiplicityCorrection::~AliMultiplicityCorrection()
 204 {
 205   //
 206   // Destructor
 207   //
 208
 209   Printf("AliMultiplicityCorrection::~AliMultiplicityCorrection called");
 210
 211   for (Int_t i = 0; i < kESDHists; ++i)
 212   {
 213     if (fMultiplicityESD[i])
 214       delete fMultiplicityESD[i];
 215     fMultiplicityESD[i] = 0;
 216   }
 217
 218   for (Int_t i = 0; i < kMCHists; ++i)
 219   {
 220     if (fMultiplicityVtx[i])
 221       delete fMultiplicityVtx[i];
 222     fMultiplicityVtx[i] = 0;
 223
 224     if (fMultiplicityMB[i])
 225       delete fMultiplicityMB[i];
 226     fMultiplicityMB[i] = 0;
 227
 228     if (fMultiplicityINEL[i])
 229       delete fMultiplicityINEL[i];
 230     fMultiplicityINEL[i] = 0;
 231   }
 232
 233   for (Int_t i = 0; i < kCorrHists; ++i)
 234   {
 235     if (fCorrelation[i])
 236       delete fCorrelation[i];
 237     fCorrelation[i] = 0;
 238
 239     if (fMultiplicityESDCorrected[i])
 240       delete fMultiplicityESDCorrected[i];
 241     fMultiplicityESDCorrected[i] = 0;
 242   }
 243 }
 244
 245 //____________________________________________________________________
 246 Long64_t AliMultiplicityCorrection::Merge(TCollection* list)
 247 {
 248   // Merge a list of AliMultiplicityCorrection objects with this (needed for
 249   // PROOF).
 250   // Returns the number of merged objects (including this).
 251
 252   if (!list)
 253     return 0;
 254
 255   if (list->IsEmpty())
 256     return 1;
 257
 258   TIterator* iter = list->MakeIterator();
 259   TObject* obj;
 260
 261   // collections of all histograms
 262   TList collections[kESDHists+kMCHists*3+kCorrHists*2];
 263
 264   Int_t count = 0;
 265   while ((obj = iter->Next())) {
 266
 267     AliMultiplicityCorrection* entry = dynamic_cast<AliMultiplicityCorrection*> (obj);
 268     if (entry == 0)
 269       continue;
 270
 271     for (Int_t i = 0; i < kESDHists; ++i)
 272       collections[i].Add(entry->fMultiplicityESD[i]);
 273
 274     for (Int_t i = 0; i < kMCHists; ++i)
 275     {
 276       collections[kESDHists+i].Add(entry->fMultiplicityVtx[i]);
 277       collections[kESDHists+kMCHists+i].Add(entry->fMultiplicityMB[i]);
 278       collections[kESDHists+kMCHists*2+i].Add(entry->fMultiplicityINEL[i]);
 279     }
 280
 281     for (Int_t i = 0; i < kCorrHists; ++i)
 282       collections[kESDHists+kMCHists*3+i].Add(entry->fCorrelation[i]);
 283
 284     for (Int_t i = 0; i < kCorrHists; ++i)
 285       collections[kESDHists+kMCHists*3+kCorrHists+i].Add(entry->fMultiplicityESDCorrected[i]);
 286
 287     count++;
 288   }
 289
 290   for (Int_t i = 0; i < kESDHists; ++i)
 291     fMultiplicityESD[i]->Merge(&collections[i]);
 292
 293   for (Int_t i = 0; i < kMCHists; ++i)
 294   {
 295     fMultiplicityVtx[i]->Merge(&collections[kESDHists+i]);
 296     fMultiplicityMB[i]->Merge(&collections[kESDHists+kMCHists+i]);
 297     fMultiplicityINEL[i]->Merge(&collections[kESDHists+kMCHists*2+i]);
 298   }
 299
 300   for (Int_t i = 0; i < kCorrHists; ++i)
 301     fCorrelation[i]->Merge(&collections[kESDHists+kMCHists*3+i]);
 302
 303   for (Int_t i = 0; i < kCorrHists; ++i)
 304     fMultiplicityESDCorrected[i]->Merge(&collections[kESDHists+kMCHists*3+kCorrHists+i]);
 305
 306   delete iter;
 307
 308   return count+1;
 309 }
 310
 311 //____________________________________________________________________
 312 Bool_t AliMultiplicityCorrection::LoadHistograms(const Char_t* dir)
 313 {
 314   //
 315   // loads the histograms from a file
 316   // if dir is empty a directory with the name of this object is taken (like in SaveHistogram)
 317   //
 318
 319   if (!dir)
 320     dir = GetName();
 321
 322   if (!gDirectory->cd(dir))
 323     return kFALSE;
 324
 325   // store old hists to delete them later
 326   TList oldObjects;
 327   oldObjects.SetOwner(1);
 328   for (Int_t i = 0; i < kESDHists; ++i)
 329     if (fMultiplicityESD[i])
 330       oldObjects.Add(fMultiplicityESD[i]);
 331
 332   for (Int_t i = 0; i < kMCHists; ++i)
 333   {
 334     if (fMultiplicityVtx[i])
 335       oldObjects.Add(fMultiplicityVtx[i]);
 336     if (fMultiplicityMB[i])
 337       oldObjects.Add(fMultiplicityMB[i]);
 338     if (fMultiplicityINEL[i])
 339       oldObjects.Add(fMultiplicityINEL[i]);
 340   }
 341
 342   for (Int_t i = 0; i < kCorrHists; ++i)
 343     if (fCorrelation[i])
 344       oldObjects.Add(fCorrelation[i]);
 345
 346   // load histograms
 347
 348   Bool_t success = kTRUE;
 349
 350   for (Int_t i = 0; i < kESDHists; ++i)
 351   {
 352     fMultiplicityESD[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityESD[i]->GetName()));
 353     if (!fMultiplicityESD[i])
 354       success = kFALSE;
 355   }
 356
 357   for (Int_t i = 0; i < kMCHists; ++i)
 358   {
 359     fMultiplicityVtx[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityVtx[i]->GetName()));
 360     fMultiplicityMB[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityMB[i]->GetName()));
 361     fMultiplicityINEL[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityINEL[i]->GetName()));
 362     if (!fMultiplicityVtx[i] || !fMultiplicityMB[i] || !fMultiplicityINEL[i])
 363       success = kFALSE;
 364   }
 365
 366   for (Int_t i = 0; i < kCorrHists; ++i)
 367   {
 368     fCorrelation[i] = dynamic_cast<TH3F*> (gDirectory->Get(fCorrelation[i]->GetName()));
 369     if (!fCorrelation[i])
 370       success = kFALSE;
 371     fMultiplicityESDCorrected[i] = dynamic_cast<TH1F*> (gDirectory->Get(fMultiplicityESDCorrected[i]->GetName()));
 372     if (!fMultiplicityESDCorrected[i])
 373       success = kFALSE;
 374   }
 375
 376   gDirectory->cd("..");
 377
 378   // delete old hists
 379   oldObjects.Delete();
 380
 381   return success;
 382 }
 383
 384 //____________________________________________________________________
 385 void AliMultiplicityCorrection::SaveHistograms(const char* dir)
 386 {
 387   //
 388   // saves the histograms
 389   //
 390
 391   if (!dir)
 392     dir = GetName();
 393
 394   gDirectory->mkdir(dir);
 395   gDirectory->cd(dir);
 396
 397   for (Int_t i = 0; i < kESDHists; ++i)
 398     if (fMultiplicityESD[i])
 399       fMultiplicityESD[i]->Write();
 400
 401   for (Int_t i = 0; i < kMCHists; ++i)
 402   {
 403     if (fMultiplicityVtx[i])
 404       fMultiplicityVtx[i]->Write();
 405     if (fMultiplicityMB[i])
 406       fMultiplicityMB[i]->Write();
 407     if (fMultiplicityINEL[i])
 408       fMultiplicityINEL[i]->Write();
 409   }
 410
 411   for (Int_t i = 0; i < kCorrHists; ++i)
 412   {
 413     if (fCorrelation[i])
 414       fCorrelation[i]->Write();
 415     if (fMultiplicityESDCorrected[i])
 416       fMultiplicityESDCorrected[i]->Write();
 417   }
 418
 419   gDirectory->cd("..");
 420 }
 421
 422 //____________________________________________________________________
 423 void AliMultiplicityCorrection::FillGenerated(Float_t vtx, Bool_t triggered, Bool_t vertex, Int_t generated05, Int_t generated10, Int_t generated15, Int_t generated20, Int_t generatedAll)
 424 {
 425   //
 426   // Fills an event from MC
 427   //
 428
 429   if (triggered)
 430   {
 431     fMultiplicityMB[0]->Fill(vtx, generated05);
 432     fMultiplicityMB[1]->Fill(vtx, generated10);
 433     fMultiplicityMB[2]->Fill(vtx, generated15);
 434     fMultiplicityMB[3]->Fill(vtx, generated20);
 435     fMultiplicityMB[4]->Fill(vtx, generatedAll);
 436
 437     if (vertex)
 438     {
 439       fMultiplicityVtx[0]->Fill(vtx, generated05);
 440       fMultiplicityVtx[1]->Fill(vtx, generated10);
 441       fMultiplicityVtx[2]->Fill(vtx, generated15);
 442       fMultiplicityVtx[3]->Fill(vtx, generated20);
 443       fMultiplicityVtx[4]->Fill(vtx, generatedAll);
 444     }
 445   }
 446
 447   fMultiplicityINEL[0]->Fill(vtx, generated05);
 448   fMultiplicityINEL[1]->Fill(vtx, generated10);
 449   fMultiplicityINEL[2]->Fill(vtx, generated15);
 450   fMultiplicityINEL[3]->Fill(vtx, generated20);
 451   fMultiplicityINEL[4]->Fill(vtx, generatedAll);
 452 }
 453
 454 //____________________________________________________________________
 455 void AliMultiplicityCorrection::FillMeasured(Float_t vtx, Int_t measured05, Int_t measured10, Int_t measured15, Int_t measured20)
 456 {
 457   //
 458   // Fills an event from ESD
 459   //
 460
 461   fMultiplicityESD[0]->Fill(vtx, measured05);
 462   fMultiplicityESD[1]->Fill(vtx, measured10);
 463   fMultiplicityESD[2]->Fill(vtx, measured15);
 464   fMultiplicityESD[3]->Fill(vtx, measured20);
 465 }
 466
 467 //____________________________________________________________________
 468 void AliMultiplicityCorrection::FillCorrection(Float_t vtx, Int_t generated05, Int_t generated10, Int_t generated15, Int_t generated20, Int_t generatedAll, Int_t measured05, Int_t measured10, Int_t measured15, Int_t measured20)
 469 {
 470   //
 471   // Fills an event into the correlation map with the information from MC and ESD
 472   //
 473
 474   fCorrelation[0]->Fill(vtx, generated05, measured05);
 475   fCorrelation[1]->Fill(vtx, generated10, measured10);
 476   fCorrelation[2]->Fill(vtx, generated15, measured15);
 477   fCorrelation[3]->Fill(vtx, generated20, measured20);
 478
 479   fCorrelation[4]->Fill(vtx, generatedAll, measured05);
 480   fCorrelation[5]->Fill(vtx, generatedAll, measured10);
 481   fCorrelation[6]->Fill(vtx, generatedAll, measured15);
 482   fCorrelation[7]->Fill(vtx, generatedAll, measured20);
 483 }
 484
 485 //____________________________________________________________________
 486 Double_t AliMultiplicityCorrection::RegularizationPol0(TVectorD& params)
 487 {
 488   // homogenity term for minuit fitting
 489   // pure function of the parameters
 490   // prefers constant function (pol0)
 491
 492   Double_t chi2 = 0;
 493
 494   // ignore the first bin here. on purpose...
 495   for (Int_t i=2; i<fgNParamsMinuit; ++i)
 496   {
 497     Double_t right  = params[i];
 498     Double_t left   = params[i-1];
 499
 500     if (right != 0)
 501     {
 502       Double_t diff = 1 - left / right;
 503       chi2 += diff * diff;
 504     }
 505   }
 506
 507   return chi2 / 100.0;
 508 }
 509
 510 //____________________________________________________________________
 511 Double_t AliMultiplicityCorrection::RegularizationPol1(TVectorD& params)
 512 {
 513   // homogenity term for minuit fitting
 514   // pure function of the parameters
 515   // prefers linear function (pol1)
 516
 517   Double_t chi2 = 0;
 518
 519   // ignore the first bin here. on purpose...
 520   for (Int_t i=2+1; i<fgNParamsMinuit; ++i)
 521   {
 522     if (params[i-1] == 0)
 523       continue;
 524
 525     Double_t right  = params[i];
 526     Double_t middle = params[i-1];
 527     Double_t left   = params[i-2];
 528
 529     Double_t der1 = (right - middle);
 530     Double_t der2 = (middle - left);
 531
 532     Double_t diff = (der1 - der2) / middle;
 533
 534     chi2 += diff * diff;
 535   }
 536
 537   return chi2;
 538 }
 539
 540 //____________________________________________________________________
 541 Double_t AliMultiplicityCorrection::RegularizationLog(TVectorD& params)
 542 {
 543   // homogenity term for minuit fitting
 544   // pure function of the parameters
 545   // prefers linear function (pol1)
 546
 547   Double_t chi2 = 0;
 548
 549   /*Float_t der[fgkMaxParams];
 550
 551   for (Int_t i=0; i<fgkMaxParams-1; ++i)
 552     der[i] = params[i+1] - params[i];
 553
 554   for (Int_t i=0; i<fgkMaxParams-6; ++i)
 555   {
 556     for (Int_t j=1; j<=5; ++j)
 557     {
 558       Double_t diff = der[i+j] - der[i];
 559       chi2 += diff * diff;
 560     }
 561   }*/
 562
 563   // ignore the first bin here. on purpose...
 564   for (Int_t i=2+1; i<fgNParamsMinuit; ++i)
 565   {
 566     if (params[i-1] == 0)
 567       continue;
 568
 569     Double_t right  = log(params[i]);
 570     Double_t middle = log(params[i-1]);
 571     Double_t left   = log(params[i-2]);
 572
 573     Double_t der1 = (right - middle);
 574     Double_t der2 = (middle - left);
 575
 576     Double_t diff = (der1 - der2) / middle;
 577
 578     chi2 += diff * diff;
 579   }
 580
 581   return chi2 * 100;
 582 }
 583
 584 //____________________________________________________________________
 585 Double_t AliMultiplicityCorrection::RegularizationTotalCurvature(TVectorD& params)
 586 {
 587   // homogenity term for minuit fitting
 588   // pure function of the parameters
 589   // minimizes the total curvature (from Unfolding Methods In High-Energy Physics Experiments,
 590   // V. Blobel (Hamburg U.) . DESY 84/118, Dec 1984. 40pp.
 591
 592   Double_t chi2 = 0;
 593
 594   // ignore the first bin here. on purpose...
 595   for (Int_t i=3; i<fgNParamsMinuit; ++i)
 596   {
 597     Double_t right  = params[i];
 598     Double_t middle = params[i-1];
 599     Double_t left   = params[i-2];
 600
 601     Double_t der1 = (right - middle);
 602     Double_t der2 = (middle - left);
 603
 604     Double_t diff = (der1 - der2);
 605
 606     chi2 += diff * diff;
 607   }
 608
 609   return chi2 * 1e4;
 610 }
 611
 612 //____________________________________________________________________
 613 Double_t AliMultiplicityCorrection::RegularizationEntropy(TVectorD& params)
 614 {
 615   // homogenity term for minuit fitting
 616   // pure function of the parameters
 617   // calculates entropy, from
 618   // The method of reduced cross-entropy (M. Schmelling 1993)
 619
 620   Double_t paramSum = 0;
 621   // ignore the first bin here. on purpose...
 622   for (Int_t i=1; i<fgNParamsMinuit; ++i)
 623     paramSum += params[i];
 624
 625   Double_t chi2 = 0;
 626   for (Int_t i=1; i<fgNParamsMinuit; ++i)
 627   {
 628     //Double_t tmp = params[i] / paramSum;
 629     Double_t tmp = params[i];
 630     if (tmp > 0 && (*fgEntropyAPriori)[i] > 0)
 631       chi2 += tmp * TMath::Log(tmp / (*fgEntropyAPriori)[i]);
 632   }
 633
 634   return 100.0 + chi2;
 635 }
 636
 637 //____________________________________________________________________
 638 void AliMultiplicityCorrection::MinuitNBD(Int_t& unused1, Double_t* unused2, Double_t& chi2, Double_t *params, Int_t unused3)
 639 {
 640   //
 641   // fit function for minuit
 642   // does: nbd
 643   // func = new TF1("nbd", "[0] * TMath::Binomial([2]+TMath::Nint(x)-1, [2]-1) * pow([1] / ([1]+[2]), TMath::Nint(x)) * pow(1 + [1]/[2], -[2])", 0, 50);
 644   // func->SetParNames("scaling", "averagen", "k");
 645   // func->SetParLimits(0, 0.001, fCurrentESD->GetMaximum() * 1000);
 646   // func->SetParLimits(1, 0.001, 1000);
 647   // func->SetParLimits(2, 0.001, 1000);
 648   //
 649
 650   fgNBD->SetParameters(params[0], params[1], params[2]);
 651
 652   Double_t params2[fgkMaxParams];
 653
 654   for (Int_t i=0; i<fgkMaxParams; ++i)
 655     params2[i] = fgNBD->Eval(i);
 656
 657   MinuitFitFunction(unused1, unused2, chi2, params2, unused3);
 658
 659   printf("%f %f %f --> %f\n", params[0], params[1], params[2], chi2);
 660 }
 661
 662 //____________________________________________________________________
 663 void AliMultiplicityCorrection::MinuitFitFunction(Int_t&, Double_t*, Double_t& chi2, Double_t *params, Int_t)
 664 {
 665   //
 666   // fit function for minuit
 667   // does: (m - Ad)W(m - Ad) where m = measured, A correlation matrix, d = guess, W = covariance matrix
 668   //
 669
 670   // d
 671   static TVectorD paramsVector(fgNParamsMinuit);
 672   for (Int_t i=0; i<fgNParamsMinuit; ++i)
 673     paramsVector[i] = params[i] * params[i];
 674
 675   // calculate penalty factor
 676   Double_t penaltyVal = 0;
 677   switch (fgRegularizationType)
 678   {
 679     case kNone:       break;
 680     case kPol0:       penaltyVal = RegularizationPol0(paramsVector); break;
 681     case kPol1:       penaltyVal = RegularizationPol1(paramsVector); break;
 682     case kCurvature:  penaltyVal = RegularizationTotalCurvature(paramsVector); break;
 683     case kEntropy:    penaltyVal = RegularizationEntropy(paramsVector); break;
 684     case kLog:        penaltyVal = RegularizationLog(paramsVector); break;
 685   }
 686
 687   //if (penaltyVal > 1e10)
 688   //  paramsVector2.Print();
 689
 690   penaltyVal *= fgRegularizationWeight;
 691
 692   // Ad
 693   TVectorD measGuessVector(fgkMaxInput);
 694   measGuessVector = (*fgCorrelationMatrix) * paramsVector;
 695
 696   // Ad - m
 697   measGuessVector -= (*fgCurrentESDVector);
 698
 699   //measGuessVector.Print();
 700
 701   TVectorD copy(measGuessVector);
 702
 703   // (Ad - m) W
 704   // this step can be optimized because currently only the diagonal elements of fgCorrelationCovarianceMatrix are used
 705   // normal way is like this:
 706   // measGuessVector *= (*fgCorrelationCovarianceMatrix);
 707   // optimized way like this:
 708   for (Int_t i=0; i<fgkMaxInput; ++i)
 709     measGuessVector[i] *= (*fgCorrelationCovarianceMatrix)(i, i);
 710
 711   //measGuessVector.Print();
 712
 713   // reduce influence of first bin
 714   copy(0) *= 0.01;
 715
 716   // (Ad - m) W (Ad - m)
 717   // the factor 1e6 prevents that a very small number (measGuessVector[i]) is multiplied with a very
 718   // big number ((*fgCorrelationCovarianceMatrix)(i, i)) (see ApplyMinuitFit)
 719   Double_t chi2FromFit = measGuessVector * copy * 1e6;
 720
 721   chi2 = chi2FromFit + penaltyVal;
 722
 723   static Int_t callCount = 0;
 724   if ((callCount++ % 10000) == 0)
 725     printf("%d) %f %f --> %f\n", callCount, chi2FromFit, penaltyVal, chi2);
 726 }
 727
 728 //____________________________________________________________________
 729 void AliMultiplicityCorrection::DrawGuess(Double_t *params)
 730 {
 731   //
 732   // draws residuals of solution suggested by params and effect of regularization
 733   //
 734
 735   // d
 736   static TVectorD paramsVector(fgNParamsMinuit);
 737   for (Int_t i=0; i<fgNParamsMinuit; ++i)
 738     paramsVector[i] = params[i] * params[i];
 739
 740   // Ad
 741   TVectorD measGuessVector(fgkMaxInput);
 742   measGuessVector = (*fgCorrelationMatrix) * paramsVector;
 743
 744   // Ad - m
 745   measGuessVector -= (*fgCurrentESDVector);
 746
 747   TVectorD copy(measGuessVector);
 748
 749   // (Ad - m) W
 750   // this step can be optimized because currently only the diagonal elements of fgCorrelationCovarianceMatrix are used
 751   // normal way is like this:
 752   // measGuessVector *= (*fgCorrelationCovarianceMatrix);
 753   // optimized way like this:
 754   for (Int_t i=0; i<fgkMaxInput; ++i)
 755     measGuessVector[i] *= (*fgCorrelationCovarianceMatrix)(i, i);
 756
 757   // (Ad - m) W (Ad - m)
 758   // the factor 1e6 prevents that a very small number (measGuessVector[i]) is multiplied with a very
 759   // big number ((*fgCorrelationCovarianceMatrix)(i, i)) (see ApplyMinuitFit)
 760   //Double_t chi2FromFit = measGuessVector * copy * 1e6;
 761
 762   // draw residuals
 763   TH1* residuals = new TH1F("residuals", "residuals", fgkMaxInput, -0.5, fgkMaxInput - 0.5);
 764   for (Int_t i=0; i<fgkMaxInput; ++i)
 765     residuals->SetBinContent(i+1, measGuessVector(i) * copy(i) * 1e6);
 766   new TCanvas; residuals->DrawCopy(); gPad->SetLogy();
 767
 768   // draw penalty
 769   if (fgRegularizationType != kPol1) {
 770     Printf("Drawing not supported");
 771     return;
 772   }
 773
 774   TH1* penalty = new TH1F("penalty", "penalty", fgNParamsMinuit, -0.5, fgNParamsMinuit - 0.5);
 775
 776   for (Int_t i=2; i<fgNParamsMinuit; ++i)
 777   {
 778     if (params[i-1] == 0)
 779       continue;
 780
 781     Double_t right  = params[i];
 782     Double_t middle = params[i-1];
 783     Double_t left   = params[i-2];
 784
 785     Double_t der1 = (right - middle);
 786     Double_t der2 = (middle - left);
 787
 788     Double_t diff = (der1 - der2) / middle;
 789
 790     penalty->SetBinContent(i, diff * diff);
 791   }
 792   new TCanvas; penalty->DrawCopy(); gPad->SetLogy();
 793 }
 794
 795 //____________________________________________________________________
 796 void AliMultiplicityCorrection::SetupCurrentHists(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t createBigBin)
 797 {
 798   //
 799   // fills fCurrentESD, fCurrentCorrelation
 800   // resets fMultiplicityESDCorrected
 801   //
 802
 803   Bool_t display = kFALSE;
 804
 805   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
 806
 807   fMultiplicityESDCorrected[correlationID]->Reset();
 808   fMultiplicityESDCorrected[correlationID]->Sumw2();
 809
 810   // project without under/overflow bins
 811   fCurrentESD = fMultiplicityESD[inputRange]->ProjectionY("fCurrentESD", 1, fMultiplicityESD[inputRange]->GetXaxis()->GetNbins());
 812   fCurrentESD->Sumw2();
 813
 814   // empty under/overflow bins in x, otherwise Project3D takes them into account
 815   TH3* hist = fCorrelation[correlationID];
 816   for (Int_t y=0; y<=hist->GetYaxis()->GetNbins()+1; ++y)
 817   {
 818     for (Int_t z=0; z<=hist->GetZaxis()->GetNbins()+1; ++z)
 819     {
 820       hist->SetBinContent(0, y, z, 0);
 821       hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0);
 822     }
 823   }
 824
 825   fCurrentCorrelation = hist->Project3D("zy");
 826   //((TH2*) fCurrentCorrelation)->Rebin2D(2, 1);
 827   //fMultiplicityESDCorrected[correlationID]->Rebin(2);
 828   fCurrentCorrelation->Sumw2();
 829
 830   Printf("AliMultiplicityCorrection::SetupCurrentHists: Statistics information: %.f entries in correlation map; %.f entries in measured spectrum", fCurrentCorrelation->Integral(), fCurrentESD->Integral());
 831
 832   if (createBigBin)
 833   {
 834     fLastBinLimit = 0;
 835     for (Int_t i=1; i<=fCurrentESD->GetNbinsX(); ++i)
 836     {
 837       if (fCurrentESD->GetBinContent(i) <= 5)
 838       {
 839         fLastBinLimit = i;
 840         break;
 841       }
 842     }
 843
 844     if (fLastBinLimit > 0)
 845     {
 846       TCanvas* canvas = 0;
 847
 848       if (display)
 849       {
 850         canvas = new TCanvas("StatSolution", "StatSolution", 1000, 800);
 851         canvas->Divide(2, 2);
 852
 853         canvas->cd(1);
 854         fCurrentESD->GetXaxis()->SetRangeUser(0, 200);
 855         fCurrentESD->SetStats(kFALSE);
 856         fCurrentESD->SetTitle(";measured multiplicity");
 857         fCurrentESD->DrawCopy();
 858         gPad->SetLogy();
 859
 860         canvas->cd(2);
 861         fCurrentCorrelation->GetXaxis()->SetRangeUser(0, 250);
 862         fCurrentCorrelation->GetYaxis()->SetRangeUser(0, 250);
 863         fCurrentCorrelation->SetTitle(";true multiplicity;measured multiplicity");
 864         fCurrentCorrelation->SetStats(kFALSE);
 865         fCurrentCorrelation->DrawCopy("COLZ");
 866       }
 867
 868       printf("Bin limit in measured spectrum is %d.\n", fLastBinLimit);
 869       fCurrentESD->SetBinContent(fLastBinLimit, fCurrentESD->Integral(fLastBinLimit, fCurrentESD->GetNbinsX()));
 870       for (Int_t i=fLastBinLimit+1; i<=fCurrentESD->GetNbinsX(); ++i)
 871       {
 872         fCurrentESD->SetBinContent(i, 0);
 873         fCurrentESD->SetBinError(i, 0);
 874       }
 875       // the error is set to sqrt(N), better solution possible?, sum of relative errors of all contributions???
 876       fCurrentESD->SetBinError(fLastBinLimit, TMath::Sqrt(fCurrentESD->GetBinContent(fLastBinLimit)));
 877
 878       printf("This bin has now %f +- %f entries\n", fCurrentESD->GetBinContent(fLastBinLimit), fCurrentESD->GetBinError(fLastBinLimit));
 879
 880       for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
 881       {
 882         fCurrentCorrelation->SetBinContent(i, fLastBinLimit, fCurrentCorrelation->Integral(i, i, fLastBinLimit, fCurrentCorrelation->GetNbinsY()));
 883         // the error is set to sqrt(N), better solution possible?, sum of relative errors of all contributions???
 884         fCurrentCorrelation->SetBinError(i, fLastBinLimit, TMath::Sqrt(fCurrentCorrelation->GetBinContent(i, fLastBinLimit)));
 885
 886         for (Int_t j=fLastBinLimit+1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
 887         {
 888           fCurrentCorrelation->SetBinContent(i, j, 0);
 889           fCurrentCorrelation->SetBinError(i, j, 0);
 890         }
 891       }
 892
 893       printf("Adjusted correlation matrix!\n");
 894
 895       if (canvas)
 896       {
 897         canvas->cd(3);
 898         fCurrentESD->DrawCopy();
 899         gPad->SetLogy();
 900
 901         canvas->cd(4);
 902         fCurrentCorrelation->DrawCopy("COLZ");
 903       }
 904     }
 905   }
 906
 907 #if 0 // does not help
 908   // null bins with one entry
 909   Int_t nNulledBins = 0;
 910   for (Int_t x=1; x<=fCurrentCorrelation->GetXaxis()->GetNbins(); ++x)
 911     for (Int_t y=1; y<=fCurrentCorrelation->GetYaxis()->GetNbins(); ++y)
 912     {
 913       if (fCurrentCorrelation->GetBinContent(x, y) == 1)
 914       {
 915         fCurrentCorrelation->SetBinContent(x, y, 0);
 916         fCurrentCorrelation->SetBinError(x, y, 0);
 917
 918         ++nNulledBins;
 919       }
 920     }
 921   Printf("Nulled %d bins", nNulledBins);
 922 #endif
 923
 924   fCurrentEfficiency = GetEfficiency(inputRange, eventType);
 925   //fCurrentEfficiency->Rebin(2);
 926   //fCurrentEfficiency->Scale(0.5);
 927 }
 928
 929 //____________________________________________________________________
 930 TH1* AliMultiplicityCorrection::GetEfficiency(Int_t inputRange, EventType eventType)
 931 {
 932   //
 933   // calculates efficiency for given event type
 934   //
 935
 936   TH1* divisor = 0;
 937   switch (eventType)
 938   {
 939     case kTrVtx : divisor = fMultiplicityVtx[inputRange]->ProjectionY("divisor", -1, -1, "e"); break;
 940     case kMB: divisor = fMultiplicityMB[inputRange]->ProjectionY("divisor", -1, -1, "e"); break;
 941     case kINEL: divisor = fMultiplicityINEL[inputRange]->ProjectionY("divisor", -1, -1, "e"); break;
 942   }
 943   TH1* eff = fMultiplicityVtx[inputRange]->ProjectionY("CurrentEfficiency", -1, -1, "e");
 944   eff->Divide(eff, divisor, 1, 1, "B");
 945   return eff;
 946 }
 947
 948 //____________________________________________________________________
 949 void AliMultiplicityCorrection::SetRegularizationParameters(RegularizationType type, Float_t weight, Int_t minuitParams)
 950 {
 951   //
 952   // sets the parameters for chi2 minimization
 953   //
 954
 955   fgRegularizationType = type;
 956   fgRegularizationWeight = weight;
 957
 958   printf("AliMultiplicityCorrection::SetRegularizationParameters --> Regularization set to %d with weight %f\n", (Int_t) type, weight);
 959
 960   if (minuitParams > 0)
 961   {
 962     fgNParamsMinuit = minuitParams;
 963     printf("AliMultiplicityCorrection::SetRegularizationParameters --> Number of MINUIT iterations set to %d\n", minuitParams);
 964   }
 965 }
 966
 967 //____________________________________________________________________
 968 void AliMultiplicityCorrection::SetBayesianParameters(Float_t smoothing, Int_t nIterations)
 969 {
 970   //
 971   // sets the parameters for Bayesian unfolding
 972   //
 973
 974   fgBayesianSmoothing = smoothing;
 975   fgBayesianIterations = nIterations;
 976
 977   printf("AliMultiplicityCorrection::SetBayesianParameters --> Paramaters set to %d iterations with smoothing %f\n", fgBayesianIterations, fgBayesianSmoothing);
 978 }
 979
 980 //____________________________________________________________________
 981 Int_t AliMultiplicityCorrection::UnfoldWithMinuit(TH1* correlation, TH1* aEfficiency, TH1* measured, TH1* initialConditions, TH1* result, Bool_t check)
 982 {
 983   //
 984   // implemenation of unfolding (internal function)
 985   //
 986   // unfolds <measured> using response from <correlation> and effiency <aEfficiency>
 987   // output is in <result>
 988   // <initialConditions> set the initial values for the minimization, if 0 <measured> is used
 989   // if <check> is true no unfolding is made, instead only the chi2 without unfolding is printed
 990   //
 991   // returns minuit status (0 = success), (-1 when check was set)
 992   //
 993
 994   //normalize measured
 995   measured->Scale(1.0 / measured->Integral());
 996
 997   if (!fgCorrelationMatrix)
 998     fgCorrelationMatrix = new TMatrixD(fgkMaxInput, fgNParamsMinuit);
 999   if (!fgCorrelationCovarianceMatrix)
1000     fgCorrelationCovarianceMatrix = new TMatrixD(fgkMaxInput, fgkMaxInput);
1001   if (!fgCurrentESDVector)
1002     fgCurrentESDVector = new TVectorD(fgkMaxInput);
1003   if (!fgEntropyAPriori)
1004     fgEntropyAPriori = new TVectorD(fgNParamsMinuit);
1005
1006   fgCorrelationMatrix->Zero();
1007   fgCorrelationCovarianceMatrix->Zero();
1008   fgCurrentESDVector->Zero();
1009   fgEntropyAPriori->Zero();
1010
1011   // normalize correction for given nPart
1012   for (Int_t i=1; i<=correlation->GetNbinsX(); ++i)
1013   {
1014     Double_t sum = correlation->Integral(i, i, 1, correlation->GetNbinsY());
1015     if (sum <= 0)
1016       continue;
1017     Float_t maxValue = 0;
1018     Int_t maxBin = -1;
1019     for (Int_t j=1; j<=correlation->GetNbinsY(); ++j)
1020     {
1021       // find most probably value
1022       if (maxValue < correlation->GetBinContent(i, j))
1023       {
1024         maxValue = correlation->GetBinContent(i, j);
1025         maxBin = j;
1026       }
1027
1028       // npart sum to 1
1029       correlation->SetBinContent(i, j, correlation->GetBinContent(i, j) / sum);
1030       correlation->SetBinError(i, j, correlation->GetBinError(i, j) / sum);
1031
1032       if (i <= fgNParamsMinuit && j <= fgkMaxInput)
1033         (*fgCorrelationMatrix)(j-1, i-1) = correlation->GetBinContent(i, j);
1034     }
1035
1036     //printf("MPV for Ntrue = %f is %f\n", fCurrentCorrelation->GetXaxis()->GetBinCenter(i), fCurrentCorrelation->GetYaxis()->GetBinCenter(maxBin));
1037   }
1038
1039   // Initialize TMinuit via generic fitter interface
1040   TVirtualFitter *minuit = TVirtualFitter::Fitter(0, fgNParamsMinuit);
1041   Double_t arglist[100];
1042
1043   // disable any output (-1), unfortuantly we do not see warnings anymore then. Have to find another way...
1044   arglist[0] = 0;
1045   minuit->ExecuteCommand("SET PRINT", arglist, 1);
1046
1047   // however, enable warnings
1048   //minuit->ExecuteCommand("SET WAR", arglist, 0);
1049
1050   // set minimization function
1051   minuit->SetFCN(MinuitFitFunction);
1052
1053   for (Int_t i=0; i<fgNParamsMinuit; i++)
1054     (*fgEntropyAPriori)[i] = 1;
1055
1056   if (!initialConditions) {
1057     initialConditions = measured;
1058   } else {
1059     printf("Using different starting conditions...\n");
1060     //new TCanvas; initialConditions->DrawCopy();
1061     initialConditions->Scale(1.0 / initialConditions->Integral());
1062   }
1063
1064   for (Int_t i=0; i<fgNParamsMinuit; i++)
1065     if (initialConditions->GetBinContent(i+1) > 0)
1066       (*fgEntropyAPriori)[i] = initialConditions->GetBinContent(i+1);
1067
1068   Double_t* results = new Double_t[fgNParamsMinuit+1];
1069   for (Int_t i=0; i<fgNParamsMinuit; ++i)
1070   {
1071     results[i] = initialConditions->GetBinContent(i+1);
1072
1073     // minimum value
1074     if (!check)
1075       results[i] = TMath::Max(results[i], 1e-3);
1076
1077     Float_t stepSize = 0.1;
1078
1079     // minuit sees squared values to prevent it from going negative...
1080     results[i] = TMath::Sqrt(results[i]);
1081     //stepSize /= results[i]; // keep relative step size
1082
1083     minuit->SetParameter(i, Form("param%d", i), results[i], stepSize, 0, 0);
1084   }
1085   //minuit->SetParameter(fgkMaxParams, "alpha", 1e4, 1, 0, 0);
1086   // bin 0 is filled with value from bin 1 (otherwise it's 0)
1087   //minuit->SetParameter(0, "param0", (results[1] > minStartValue) ? results[1] : minStartValue, 0.1, 0, 0);
1088   //results[0] = 0;
1089   //minuit->SetParameter(0, "param0", 0, 0, 0, 0);
1090
1091   for (Int_t i=0; i<fgkMaxInput; ++i)
1092   {
1093     (*fgCurrentESDVector)[i] = measured->GetBinContent(i+1);
1094     if (measured->GetBinError(i+1) > 0)
1095       (*fgCorrelationCovarianceMatrix)(i, i) = (Double_t) 1e-6 / measured->GetBinError(i+1) / measured->GetBinError(i+1);
1096
1097     if ((*fgCorrelationCovarianceMatrix)(i, i) > 1e7)
1098       (*fgCorrelationCovarianceMatrix)(i, i) = 0;
1099   }
1100
1101   Int_t dummy = 0;
1102   Double_t chi2 = 0;
1103   MinuitFitFunction(dummy, 0, chi2, results, 0);
1104   DrawGuess(results);
1105   printf("Chi2 of initial parameters is = %f\n", chi2);
1106
1107   if (check)
1108   {
1109     delete[] results;
1110     return -1;
1111   }
1112
1113   // first param is number of iterations, second is precision....
1114   arglist[0] = 1e6;
1115   //arglist[1] = 1e-5;
1116   //minuit->ExecuteCommand("SCAN", arglist, 0);
1117   Int_t status = minuit->ExecuteCommand("MIGRAD", arglist, 1);
1118   printf("MINUIT status is %d\n", status);
1119   //minuit->ExecuteCommand("MIGRAD", arglist, 1);
1120   //minuit->ExecuteCommand("MIGRAD", arglist, 1);
1121   //printf("!!!!!!!!!!!!!! MIGRAD finished: Starting MINOS !!!!!!!!!!!!!!");
1122   //minuit->ExecuteCommand("MINOS", arglist, 0);
1123
1124   for (Int_t i=0; i<fgNParamsMinuit; ++i)
1125   {
1126     if (aEfficiency->GetBinContent(i+1) > 0)
1127     {
1128       results[i] = minuit->GetParameter(i);
1129       result->SetBinContent(i+1, results[i] * results[i] / aEfficiency->GetBinContent(i+1));
1130       // error is : (relError) * (value) = (minuit->GetParError(i) / minuit->GetParameter(i)) * (minuit->GetParameter(i) * minuit->GetParameter(i))
1131       result->SetBinError(i+1, minuit->GetParError(i) * results[i] /  aEfficiency->GetBinContent(i+1));
1132     }
1133     else
1134     {
1135       result->SetBinContent(i+1, 0);
1136       result->SetBinError(i+1, 0);
1137     }
1138   }
1139   DrawGuess(results);
1140
1141   delete[] results;
1142
1143   return status;
1144 }
1145
1146 //____________________________________________________________________
1147 Int_t AliMultiplicityCorrection::ApplyMinuitFit(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t check, TH1* initialConditions)
1148 {
1149   //
1150   // correct spectrum using minuit chi2 method
1151   //
1152   // for description of parameters, see UnfoldWithMinuit
1153   //
1154
1155   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
1156   SetupCurrentHists(inputRange, fullPhaseSpace, eventType, fgCreateBigBin);
1157
1158   return UnfoldWithMinuit(fCurrentCorrelation, fCurrentEfficiency, fCurrentESD, initialConditions, fMultiplicityESDCorrected[correlationID], check);
1159 }
1160
1161 //____________________________________________________________________
1162 void AliMultiplicityCorrection::ApplyNBDFit(Int_t inputRange, Bool_t fullPhaseSpace)
1163 {
1164   //
1165   // correct spectrum using minuit chi2 method applying a NBD fit
1166   //
1167
1168   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
1169
1170   SetupCurrentHists(inputRange, fullPhaseSpace, kTrVtx, kFALSE);
1171   //normalize ESD
1172   fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
1173
1174   fgCorrelationMatrix = new TMatrixD(fgkMaxParams, fgkMaxParams);
1175   fgCorrelationCovarianceMatrix = new TMatrixD(fgkMaxParams, fgkMaxParams);
1176   fgCurrentESDVector = new TVectorD(fgkMaxParams);
1177
1178   // normalize correction for given nPart
1179   for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
1180   {
1181     Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
1182     if (sum <= 0)
1183       continue;
1184     for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
1185     {
1186       // npart sum to 1
1187       fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
1188       fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
1189
1190       if (i <= fgkMaxParams && j <= fgkMaxParams)
1191         (*fgCorrelationMatrix)(j-1, i-1) = fCurrentCorrelation->GetBinContent(i, j);
1192     }
1193   }
1194
1195   for (Int_t i=0; i<fgkMaxParams; ++i)
1196   {
1197     (*fgCurrentESDVector)[i] = fCurrentESD->GetBinContent(i+1);
1198     if (fCurrentESD->GetBinError(i+1) > 0)
1199       (*fgCorrelationCovarianceMatrix)(i, i) = 1.0 / (fCurrentESD->GetBinError(i+1) * fCurrentESD->GetBinError(i+1));
1200   }
1201
1202   // Create NBD function
1203   if (!fgNBD)
1204   {
1205     fgNBD = new TF1("nbd", "[0] * TMath::Binomial([2]+TMath::Nint(x)-1, [2]-1) * pow([1] / ([1]+[2]), TMath::Nint(x)) * pow(1 + [1]/[2], -[2])", 0, 250);
1206     fgNBD->SetParNames("scaling", "averagen", "k");
1207   }
1208
1209   // Initialize TMinuit via generic fitter interface
1210   TVirtualFitter *minuit = TVirtualFitter::Fitter(0, 3);
1211
1212   minuit->SetFCN(MinuitNBD);
1213   minuit->SetParameter(0, "param0", 1, 0.1, 0.001, fCurrentESD->GetMaximum() * 1000);
1214   minuit->SetParameter(1, "param1", 10, 1, 0.001, 1000);
1215   minuit->SetParameter(2, "param2", 2, 1, 0.001, 1000);
1216
1217   Double_t arglist[100];
1218   arglist[0] = 0;
1219   minuit->ExecuteCommand("SET PRINT", arglist, 1);
1220   minuit->ExecuteCommand("MIGRAD", arglist, 0);
1221   //minuit->ExecuteCommand("MINOS", arglist, 0);
1222
1223   fgNBD->SetParameters(minuit->GetParameter(0), minuit->GetParameter(1), minuit->GetParameter(2));
1224
1225   for (Int_t i=0; i<fgkMaxParams; ++i)
1226   {
1227     printf("%d : %f\n", i, fgNBD->Eval(i));
1228     if (fgNBD->Eval(i) > 0)
1229     {
1230       fMultiplicityESDCorrected[correlationID]->SetBinContent(i+1, fgNBD->Eval(i));
1231       fMultiplicityESDCorrected[correlationID]->SetBinError(i+1, 0);
1232     }
1233   }
1234 }
1235
1236 //____________________________________________________________________
1237 void AliMultiplicityCorrection::DrawHistograms()
1238 {
1239   //
1240   // draws the histograms of this class
1241   //
1242
1243   printf("ESD:\n");
1244
1245   TCanvas* canvas1 = new TCanvas("fMultiplicityESD", "fMultiplicityESD", 900, 600);
1246   canvas1->Divide(3, 2);
1247   for (Int_t i = 0; i < kESDHists; ++i)
1248   {
1249     canvas1->cd(i+1);
1250     fMultiplicityESD[i]->DrawCopy("COLZ");
1251     printf("%d --> %f\n", i, (Float_t) fMultiplicityESD[i]->ProjectionY()->GetMean());
1252   }
1253
1254   printf("Vtx:\n");
1255
1256   TCanvas* canvas2 = new TCanvas("fMultiplicityMC", "fMultiplicityMC", 900, 600);
1257   canvas2->Divide(3, 2);
1258   for (Int_t i = 0; i < kMCHists; ++i)
1259   {
1260     canvas2->cd(i+1);
1261     fMultiplicityVtx[i]->DrawCopy("COLZ");
1262     printf("%d --> %f\n", i, fMultiplicityVtx[i]->ProjectionY()->GetMean());
1263   }
1264
1265   TCanvas* canvas3 = new TCanvas("fCorrelation", "fCorrelation", 900, 900);
1266   canvas3->Divide(3, 3);
1267   for (Int_t i = 0; i < kCorrHists; ++i)
1268   {
1269     canvas3->cd(i+1);
1270     TH3* hist = dynamic_cast<TH3*> (fCorrelation[i]->Clone());
1271     for (Int_t y=1; y<=hist->GetYaxis()->GetNbins(); ++y)
1272     {
1273       for (Int_t z=1; z<=hist->GetZaxis()->GetNbins(); ++z)
1274       {
1275         hist->SetBinContent(0, y, z, 0);
1276         hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0);
1277       }
1278     }
1279     TH1* proj = hist->Project3D("zy");
1280     proj->DrawCopy("COLZ");
1281   }
1282 }
1283
1284 //____________________________________________________________________
1285 void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRange, Bool_t fullPhaseSpace, Bool_t /*normalizeESD*/, TH1* mcHist, Bool_t simple)
1286 {
1287   // draw comparison plots
1288
1289
1290   //mcHist->Rebin(2);
1291
1292   Int_t esdCorrId = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
1293
1294   TString tmpStr;
1295   tmpStr.Form("%s_DrawComparison_%d", name, esdCorrId);
1296
1297   if (fMultiplicityESDCorrected[esdCorrId]->Integral() == 0)
1298   {
1299     printf("ERROR. Unfolded histogram is empty\n");
1300     return;
1301   }
1302
1303   //regain measured distribution used for unfolding, because the bins at high mult. were modified in SetupCurrentHists
1304   fCurrentESD = fMultiplicityESD[esdCorrId]->ProjectionY("fCurrentESD", 1, fMultiplicityESD[inputRange]->GetXaxis()->GetNbins());
1305   fCurrentESD->Sumw2();
1306   fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
1307
1308   // normalize unfolded result to 1
1309   fMultiplicityESDCorrected[esdCorrId]->Scale(1.0 / fMultiplicityESDCorrected[esdCorrId]->Integral());
1310
1311   //fCurrentESD->Scale(mcHist->Integral(2, 200));
1312
1313   //new TCanvas;
1314   /*TH1* ratio = (TH1*) fMultiplicityESDCorrected[esdCorrId]->Clone("ratio");
1315   ratio->Divide(mcHist);
1316   ratio->Draw("HIST");
1317   ratio->Fit("pol0", "W0", "", 20, 120);
1318   Float_t scalingFactor = ratio->GetFunction("pol0")->GetParameter(0);
1319   delete ratio;
1320   mcHist->Scale(scalingFactor);*/
1321
1322   // find bin with <= 50 entries. this is used as limit for the chi2 calculation
1323   Int_t mcBinLimit = 0;
1324   for (Int_t i=20; i<mcHist->GetNbinsX(); ++i)
1325   {
1326     if (mcHist->GetBinContent(i) > 50)
1327     {
1328       mcBinLimit = i;
1329     }
1330     else
1331       break;
1332   }
1333   Printf("AliMultiplicityCorrection::DrawComparison: MC bin limit is %d", mcBinLimit);
1334
1335   // scale to 1
1336   mcHist->Sumw2();
1337   if (mcHist->Integral() > 0)
1338     mcHist->Scale(1.0 / mcHist->Integral());
1339
1340   // calculate residual
1341
1342   // for that we convolute the response matrix with the gathered result
1343   // if normalizeESD is not set, the histogram is already normalized, this needs to be passed to CalculateMultiplicityESD
1344   TH1* tmpESDEfficiencyRecorrected = (TH1*) fMultiplicityESDCorrected[esdCorrId]->Clone("tmpESDEfficiencyRecorrected");
1345   tmpESDEfficiencyRecorrected->Multiply(fCurrentEfficiency);
1346   TH2* convoluted = CalculateMultiplicityESD(tmpESDEfficiencyRecorrected, esdCorrId);
1347   TH1* convolutedProj = convoluted->ProjectionY("convolutedProj", -1, -1, "e");
1348   if (convolutedProj->Integral() > 0)
1349   {
1350     convolutedProj->Scale(1.0 / convolutedProj->Integral());
1351   }
1352   else
1353     printf("ERROR: convolutedProj is empty. Something went wrong calculating the convoluted histogram.\n");
1354
1355   TH1* residual = (TH1*) convolutedProj->Clone("residual");
1356   residual->SetTitle("Residuals;Ntracks;(folded unfolded measured - measured) / e");
1357
1358   residual->Add(fCurrentESD, -1);
1359   //residual->Divide(residual, fCurrentESD, 1, 1, "B");
1360
1361   TH1* residualHist = new TH1F("residualHist", "residualHist", 51, -5, 5);
1362
1363   // TODO fix errors
1364   Float_t chi2 = 0;
1365   for (Int_t i=1; i<=residual->GetNbinsX(); ++i)
1366   {
1367     if (fCurrentESD->GetBinError(i) > 0)
1368     {
1369       Float_t value = residual->GetBinContent(i) / fCurrentESD->GetBinError(i);
1370       if (i > 1 && i <= fLastBinLimit)
1371         chi2 += value * value;
1372       Printf("%d --> %f (%f)", i, value * value, chi2);
1373       residual->SetBinContent(i, value);
1374       residualHist->Fill(value);
1375     }
1376     else
1377     {
1378       //printf("Residual bin %d set to 0\n", i);
1379       residual->SetBinContent(i, 0);
1380     }
1381     convolutedProj->SetBinError(i, 0);
1382     residual->SetBinError(i, 0);
1383   }
1384   fLastChi2Residuals = chi2;
1385
1386   new TCanvas; residualHist->DrawCopy();
1387
1388   //residualHist->Fit("gaus", "N");
1389   //delete residualHist;
1390
1391   printf("Difference (Residuals) is %f for bin 2-%d\n", fLastChi2Residuals, fLastBinLimit);
1392
1393   TCanvas* canvas1 = 0;
1394   if (simple)
1395   {
1396     canvas1 = new TCanvas(tmpStr, tmpStr, 900, 400);
1397     canvas1->Divide(2, 1);
1398   }
1399   else
1400   {
1401     canvas1 = new TCanvas(tmpStr, tmpStr, 1200, 1200);
1402     canvas1->Divide(2, 3);
1403   }
1404
1405   canvas1->cd(1);
1406   TH1* proj = (TH1*) mcHist->Clone("proj");
1407
1408   // normalize without 0 bin
1409   proj->Scale(1.0 / proj->Integral(2, proj->GetNbinsX()));
1410   Printf("Normalized without 0 bin!");
1411   proj->GetXaxis()->SetRangeUser(0, 200);
1412   proj->SetTitle(";true multiplicity;Entries");
1413   proj->SetStats(kFALSE);
1414
1415   fMultiplicityESDCorrected[esdCorrId]->SetLineColor(2);
1416   // normalize without 0 bin
1417   fMultiplicityESDCorrected[esdCorrId]->Scale(1.0 / fMultiplicityESDCorrected[esdCorrId]->Integral(2, fMultiplicityESDCorrected[esdCorrId]->GetNbinsX()));
1418   Printf("Normalized without 0 bin!");
1419
1420   if (proj->GetEntries() > 0) {
1421     proj->DrawCopy("HIST");
1422     fMultiplicityESDCorrected[esdCorrId]->DrawCopy("SAME HIST E");
1423   }
1424   else
1425     fMultiplicityESDCorrected[esdCorrId]->DrawCopy("HIST E");
1426
1427   gPad->SetLogy();
1428
1429   TLegend* legend = new TLegend(0.3, 0.8, 0.93, 0.93);
1430   legend->AddEntry(proj, "true distribution");
1431   legend->AddEntry(fMultiplicityESDCorrected[esdCorrId], "unfolded distribution");
1432   legend->SetFillColor(0);
1433   legend->Draw();
1434   // unfortunately does not work. maybe a bug? --> legend->SetTextSizePixels(14);
1435
1436   canvas1->cd(2);
1437
1438   gPad->SetLogy();
1439   fCurrentESD->GetXaxis()->SetRangeUser(0, 200);
1440   //fCurrentESD->SetLineColor(2);
1441   fCurrentESD->SetTitle(";measured multiplicity;Entries");
1442   fCurrentESD->SetStats(kFALSE);
1443   fCurrentESD->DrawCopy("HIST E");
1444
1445   convolutedProj->SetLineColor(2);
1446   //proj2->SetMarkerColor(2);
1447   //proj2->SetMarkerStyle(5);
1448   convolutedProj->DrawCopy("HIST SAME");
1449
1450   legend = new TLegend(0.3, 0.8, 0.93, 0.93);
1451   legend->AddEntry(fCurrentESD, "measured distribution");
1452   legend->AddEntry(convolutedProj, "R #otimes unfolded distribution");
1453   legend->SetFillColor(0);
1454   legend->Draw();
1455
1456   //TH1* diffMCUnfolded = dynamic_cast<TH1*> (proj->Clone("diffMCUnfolded"));
1457   //diffMCUnfolded->Add(fMultiplicityESDCorrected[esdCorrId], -1);
1458
1459   /*Float_t chi2 = 0;
1460   Float_t chi = 0;
1461   fLastChi2MCLimit = 0;
1462   Int_t limit = 0;
1463   for (Int_t i=2; i<=200; ++i)
1464   {
1465     if (proj->GetBinContent(i) != 0)
1466     {
1467       Float_t value = (proj->GetBinContent(i) - fMultiplicityESDCorrected[esdCorrId]->GetBinContent(i)) / proj->GetBinContent(i);
1468       chi2 += value * value;
1469       chi += TMath::Abs(value);
1470
1471       //printf("%d %f\n", i, chi);
1472
1473       if (chi2 < 0.2)
1474         fLastChi2MCLimit = i;
1475
1476       if (chi < 3)
1477         limit = i;
1478
1479     }
1480   }*/
1481
1482   /*chi2 = 0;
1483   Float_t chi = 0;
1484   Int_t limit = 0;
1485   for (Int_t i=1; i<=diffMCUnfolded->GetNbinsX(); ++i)
1486   {
1487     if (fMultiplicityESDCorrected[esdCorrId]->GetBinError(i) > 0)
1488     {
1489       Double_t value = diffMCUnfolded->GetBinContent(i) / fMultiplicityESDCorrected[esdCorrId]->GetBinError(i);
1490       if (value > 1e8)
1491         value = 1e8; //prevent arithmetic exception
1492       else if (value < -1e8)
1493         value = -1e8;
1494       if (i > 1)
1495       {
1496         chi2 += value * value;
1497         chi += TMath::Abs(value);
1498       }
1499       diffMCUnfolded->SetBinContent(i, value);
1500     }
1501     else
1502     {
1503       //printf("diffMCUnfolded bin %d set to 0\n", i);
1504       diffMCUnfolded->SetBinContent(i, 0);
1505     }
1506     if (chi2 < 1000)
1507       fLastChi2MCLimit = i;
1508     if (chi < 1000)
1509       limit = i;
1510     if (i == 150)
1511       fLastChi2MC = chi2;
1512   }
1513
1514   printf("limits %d %d\n", fLastChi2MCLimit, limit);
1515   fLastChi2MCLimit = limit;
1516
1517   printf("Difference (from MC) is %f for bin 2-150. Limit is %d.\n", fLastChi2MC, fLastChi2MCLimit);*/
1518
1519   if (!simple)
1520   {
1521     /*canvas1->cd(3);
1522
1523     diffMCUnfolded->SetTitle("#chi^{2};Npart;(MC - Unfolded) / e(unfolded)");
1524     //diffMCUnfolded->GetYaxis()->SetRangeUser(-20, 20);
1525     diffMCUnfolded->GetXaxis()->SetRangeUser(0, 200);
1526     diffMCUnfolded->DrawCopy("HIST");
1527
1528     TH1F* fluctuation = new TH1F("fluctuation", "fluctuation", 20, -5, 5);
1529     for (Int_t i=20; i<=diffMCUnfolded->GetNbinsX(); ++i)
1530       fluctuation->Fill(diffMCUnfolded->GetBinContent(i));
1531
1532     //new TCanvas; fluctuation->DrawCopy();
1533     delete fluctuation;*/
1534
1535     /*TLegend* legend = new TLegend(0.6, 0.7, 0.85, 0.85);
1536     legend->AddEntry(fMultiplicityESDCorrected, "ESD corrected");
1537     legend->AddEntry(fMultiplicityMC, "MC");
1538     legend->AddEntry(fMultiplicityESD, "ESD");
1539     legend->Draw();*/
1540
1541     canvas1->cd(4);
1542     residual->GetYaxis()->SetRangeUser(-5, 5);
1543     residual->GetXaxis()->SetRangeUser(0, 200);
1544     residual->DrawCopy();
1545
1546     canvas1->cd(5);
1547
1548     TH1* ratio = (TH1*) fMultiplicityESDCorrected[esdCorrId]->Clone("ratio");
1549     ratio->Divide(mcHist);
1550     ratio->SetTitle("Ratio;true multiplicity;Unfolded / MC");
1551     ratio->GetYaxis()->SetRangeUser(0.5, 1.5);
1552     ratio->GetXaxis()->SetRangeUser(0, 200);
1553     ratio->SetStats(kFALSE);
1554     ratio->Draw("HIST");
1555
1556     Double_t ratioChi2 = 0;
1557     fRatioAverage = 0;
1558     fLastChi2MCLimit = 0;
1559     for (Int_t i=2; i<=150; ++i)
1560     {
1561       Float_t value = ratio->GetBinContent(i) - 1;
1562       if (value > 1e8)
1563         value = 1e8; //prevent arithmetic exception
1564       else if (value < -1e8)
1565         value = -1e8;
1566
1567       ratioChi2 += value * value;
1568       fRatioAverage += TMath::Abs(value);
1569
1570       if (ratioChi2 < 0.1)
1571         fLastChi2MCLimit = i;
1572     }
1573     fRatioAverage /= 149;
1574
1575     printf("Sum over (ratio-1)^2 (2..150) is %f; average of |ratio-1| is %f\n", ratioChi2, fRatioAverage);
1576     // TODO FAKE
1577     fLastChi2MC = ratioChi2;
1578
1579     // FFT of ratio
1580     canvas1->cd(6);
1581     const Int_t kFFT = 128;
1582     Double_t fftReal[kFFT];
1583     Double_t fftImag[kFFT];
1584
1585     for (Int_t i=0; i<kFFT; ++i)
1586     {
1587       fftReal[i] = ratio->GetBinContent(i+1+10);
1588       // test: ;-)
1589       //fftReal[i] = cos(TMath::Pi() * 5 * 2 * i / 128);
1590       fftImag[i] = 0;
1591     }
1592
1593     FFT(-1, TMath::Nint(TMath::Log(kFFT) / TMath::Log(2)), fftReal, fftImag);
1594
1595     TH1* fftResult = (TH1*) ratio->Clone("fftResult");
1596     fftResult->SetTitle("FFT;true multiplicity;coeff. (10...137)");
1597     TH1* fftResult2 = (TH1*) ratio->Clone("fftResult2");
1598     TH1* fftResult3 = (TH1*) ratio->Clone("fftResult3");
1599     fftResult->Reset();
1600     fftResult2->Reset();
1601     fftResult3->Reset();
1602     fftResult->GetYaxis()->UnZoom();
1603     fftResult2->GetYaxis()->UnZoom();
1604     fftResult3->GetYaxis()->UnZoom();
1605
1606     //Printf("AFTER FFT");
1607     for (Int_t i=0; i<kFFT; ++i)
1608     {
1609       //Printf("%d: %f", i, fftReal[i]);
1610       fftResult->SetBinContent(i+1, fftReal[i]);
1611       /*if (i != 0 && TMath::Abs(fftReal[i]) > 0.5)
1612       {
1613         Printf("Nulled %d", i);
1614         fftReal[i] = 0;
1615       }*/
1616     }
1617
1618     fftResult->SetLineColor(1);
1619     fftResult->DrawCopy();
1620
1621
1622     //Printf("IMAG");
1623     for (Int_t i=0; i<kFFT; ++i)
1624     {
1625       //Printf("%d: %f", i, fftImag[i]);
1626       fftResult2->SetBinContent(i+1, fftImag[i]);
1627
1628       fftResult3->SetBinContent(i+1, TMath::Sqrt(fftReal[i] * fftReal[i] + fftImag[i] * fftImag[i]));
1629     }
1630
1631     fftResult2->SetLineColor(2);
1632     fftResult2->DrawCopy("SAME");
1633
1634     fftResult3->SetLineColor(4);
1635     fftResult3->DrawCopy("SAME");
1636
1637     for (Int_t i=1; i<kFFT - 1; ++i)
1638     {
1639       if (TMath::Sqrt(fftReal[i] * fftReal[i] + fftImag[i] * fftImag[i]) > 3)
1640       {
1641         fftReal[i-1] = 0;
1642         fftReal[i] = 0;
1643         fftReal[i+1] = 0;
1644         fftImag[i-1] = 0;
1645         fftImag[i] = 0;
1646         fftImag[i+1] = 0;
1647         //fftReal[i] = (fftReal[i-1] + fftReal[i+1]) / 2;
1648         //fftImag[i]  = (fftImag[i-1] + fftImag[i+1]) / 2;
1649         //Printf("Nulled %d to %f %f", i, fftReal[i], fftImag[i]);
1650       }
1651     }
1652
1653     FFT(1, TMath::Nint(TMath::Log(kFFT) / TMath::Log(2)), fftReal, fftImag);
1654
1655     TH1* fftResult4 = (TH1*) fftResult3->Clone("fftResult4");
1656     fftResult4->Reset();
1657
1658     //Printf("AFTER BACK-TRAFO");
1659     for (Int_t i=0; i<kFFT; ++i)
1660     {
1661       //Printf("%d: %f", i, fftReal[i]);
1662       fftResult4->SetBinContent(i+1+10, fftReal[i]);
1663     }
1664
1665     canvas1->cd(5);
1666     fftResult4->SetLineColor(4);
1667     fftResult4->DrawCopy("SAME");
1668
1669     // plot (MC - Unfolded) / error (MC)
1670     canvas1->cd(3);
1671
1672     TH1* diffMCUnfolded2 = dynamic_cast<TH1*> (proj->Clone("diffMCUnfolded2"));
1673     diffMCUnfolded2->Add(fMultiplicityESDCorrected[esdCorrId], -1);
1674
1675     Int_t ndfQual[kQualityRegions];
1676     for (Int_t region=0; region<kQualityRegions; ++region)
1677     {
1678       fQuality[region] = 0;
1679       ndfQual[region] = 0;
1680     }
1681
1682
1683     Double_t newChi2 = 0;
1684     Double_t newChi2Limit150 = 0;
1685     Int_t ndf = 0;
1686     for (Int_t i=1; i<=TMath::Min(diffMCUnfolded2->GetNbinsX(), fgkMaxParams+1); ++i)
1687     {
1688       Double_t value = 0;
1689       if (proj->GetBinError(i) > 0)
1690       {
1691         value = diffMCUnfolded2->GetBinContent(i) / proj->GetBinError(i);
1692         newChi2 += value * value;
1693         if (i > 1 && i <= mcBinLimit)
1694           newChi2Limit150 += value * value;
1695         ++ndf;
1696
1697         for (Int_t region=0; region<kQualityRegions; ++region)
1698           if (diffMCUnfolded2->GetXaxis()->GetBinCenter(i) >= fgQualityRegionsB[region] - 0.1 && diffMCUnfolded2->GetXaxis()->GetBinCenter(i) <= fgQualityRegionsE[region] + 0.1) // 0.1 to avoid e.g. 3.9999 < 4 problem
1699           {
1700             fQuality[region] += TMath::Abs(value);
1701             ++ndfQual[region];
1702           }
1703       }
1704
1705       diffMCUnfolded2->SetBinContent(i, value);
1706     }
1707
1708     // normalize region to the number of entries
1709     for (Int_t region=0; region<kQualityRegions; ++region)
1710     {
1711       if (ndfQual[region] > 0)
1712         fQuality[region] /= ndfQual[region];
1713       Printf("Quality parameter %d (%d <= mult <= %d) is %f with %d df", region, fgQualityRegionsB[region], fgQualityRegionsE[region], fQuality[region], ndfQual[region]);
1714     }
1715
1716     if (mcBinLimit > 1)
1717     {
1718       // TODO another FAKE
1719       fLastChi2MC = newChi2Limit150 / (mcBinLimit - 1);
1720       Printf("Chi2 (2..%d) from (MC - Unfolded) / e(MC) is: %.2f ndf is %d --> chi2 / ndf = %.2f", mcBinLimit, newChi2Limit150, mcBinLimit - 1, fLastChi2MC);
1721     }
1722     else
1723       fLastChi2MC = -1;
1724
1725     Printf("Chi2 (full range) from (MC - Unfolded) / e(MC) is: %.2f ndf is %d --> chi2 / ndf = %.2f", newChi2, ndf, ((ndf > 0) ? newChi2 / ndf : -1));
1726
1727     diffMCUnfolded2->SetTitle("#chi^{2};Npart;(MC - Unfolded) / e(MC)");
1728     diffMCUnfolded2->GetYaxis()->SetRangeUser(-5, 5);
1729     diffMCUnfolded2->GetXaxis()->SetRangeUser(0, 200);
1730     diffMCUnfolded2->DrawCopy("HIST");
1731   }
1732
1733   canvas1->SaveAs(Form("%s.gif", canvas1->GetName()));
1734 }
1735
1736 //____________________________________________________________________
1737 void AliMultiplicityCorrection::FFT(Int_t dir, Int_t m, Double_t *x, Double_t *y)
1738 {
1739 /*-------------------------------------------------------------------------
1740    This computes an in-place complex-to-complex FFT
1741    x and y are the real and imaginary arrays of 2^m points.
1742    dir =  1 gives forward transform
1743    dir = -1 gives reverse transform
1744
1745      Formula: forward
1746                   N-1
1747                   ---
1748               1   \          - j k 2 pi n / N
1749       X(n) = ---   >   x(k) e                    = forward transform
1750               N   /                                n=0..N-1
1751                   ---
1752                   k=0
1753
1754       Formula: reverse
1755                   N-1
1756                   ---
1757                   \          j k 2 pi n / N
1758       X(n) =       >   x(k) e                    = forward transform
1759                   /                                n=0..N-1
1760                   ---
1761                   k=0
1762 */
1763
1764    Long_t   nn, i, i1, j, k, i2, l, l1, l2;
1765    Double_t c1, c2, tx, ty, t1, t2, u1, u2, z;
1766
1767    /* Calculate the number of points */
1768    nn = 1;
1769    for (i = 0; i < m; i++)
1770        nn *= 2;
1771
1772    /* Do the bit reversal */
1773    i2 = nn >> 1;
1774    j = 0;
1775    for (i= 0; i < nn - 1; i++) {
1776        if (i < j) {
1777            tx = x[i];
1778            ty = y[i];
1779            x[i] = x[j];
1780            y[i] = y[j];
1781            x[j] = tx;
1782            y[j] = ty;
1783        }
1784        k = i2;
1785        while (k <= j) {
1786            j -= k;
1787            k >>= 1;
1788        }
1789        j += k;
1790    }
1791
1792    /* Compute the FFT */
1793    c1 = -1.0;
1794    c2 = 0.0;
1795    l2 = 1;
1796    for (l = 0; l < m; l++) {
1797        l1 = l2;
1798        l2 <<= 1;
1799        u1 = 1.0;
1800        u2 = 0.0;
1801        for (j = 0;j < l1; j++) {
1802            for (i = j; i < nn; i += l2) {
1803                i1 = i + l1;
1804                t1 = u1 * x[i1] - u2 * y[i1];
1805                t2 = u1 * y[i1] + u2 * x[i1];
1806                x[i1] = x[i] - t1;
1807                y[i1] = y[i] - t2;
1808                x[i] += t1;
1809                y[i] += t2;
1810            }
1811            z =  u1 * c1 - u2 * c2;
1812            u2 = u1 * c2 + u2 * c1;
1813            u1 = z;
1814        }
1815        c2 = TMath::Sqrt((1.0 - c1) / 2.0);
1816        if (dir == 1)
1817            c2 = -c2;
1818        c1 = TMath::Sqrt((1.0 + c1) / 2.0);
1819    }
1820
1821    /* Scaling for forward transform */
1822    if (dir == 1) {
1823        for (i=0;i<nn;i++) {
1824            x[i] /= (Double_t)nn;
1825            y[i] /= (Double_t)nn;
1826        }
1827    }
1828 }
1829
1830 //____________________________________________________________________
1831 void AliMultiplicityCorrection::GetComparisonResults(Float_t* mc, Int_t* mcLimit, Float_t* residuals, Float_t* ratioAverage) const
1832 {
1833   // Returns the chi2 between the MC and the unfolded ESD as well as between the ESD and the folded unfolded ESD
1834   // These values are computed during DrawComparison, thus this function picks up the
1835   // last calculation
1836
1837   if (mc)
1838     *mc = fLastChi2MC;
1839   if (mcLimit)
1840     *mcLimit = fLastChi2MCLimit;
1841   if (residuals)
1842     *residuals = fLastChi2Residuals;
1843   if (ratioAverage)
1844     *ratioAverage = fRatioAverage;
1845 }
1846
1847 //____________________________________________________________________
1848 TH2F* AliMultiplicityCorrection::GetMultiplicityMC(Int_t i, EventType eventType)
1849 {
1850   //
1851   // returns the corresponding MC spectrum
1852   //
1853
1854   switch (eventType)
1855   {
1856     case kTrVtx : return fMultiplicityVtx[i]; break;
1857     case kMB : return fMultiplicityMB[i]; break;
1858     case kINEL : return fMultiplicityINEL[i]; break;
1859   }
1860
1861   return 0;
1862 }
1863
1864 //____________________________________________________________________
1865 TH1* AliMultiplicityCorrection::CalculateStdDev(TH1** results, Int_t max)
1866 {
1867   // calculate standard deviation of (results[0] - results[k]) k=1...max-1
1868   // per bin one gets: sigma(r[0] - r[n]) / r[0]
1869
1870   TH1* standardDeviation = (TH1*) results[0]->Clone("standardDeviation");
1871   standardDeviation->Reset();
1872
1873   for (Int_t x=1; x<=results[0]->GetNbinsX(); x++)
1874   {
1875     if (results[0]->GetBinContent(x) > 0)
1876     {
1877       Double_t average = 0;
1878       for (Int_t n=1; n<max; ++n)
1879         average += results[n]->GetBinContent(x) - results[0]->GetBinContent(x);
1880       average /= max-1;
1881
1882       Double_t variance = 0;
1883       for (Int_t n=1; n<max; ++n)
1884       {
1885         Double_t value = results[n]->GetBinContent(x) - results[0]->GetBinContent(x) - average;
1886         variance += value * value;
1887       }
1888       variance /= max-1;
1889
1890       Double_t standardDev = TMath::Sqrt(variance);
1891       standardDeviation->SetBinContent(x, standardDev / results[0]->GetBinContent(x));
1892       //Printf("sigma_%d is %f value %f --> error %f", x, standardDev, results[0]->GetBinContent(x), standardDev / results[0]->GetBinContent(x));
1893     }
1894   }
1895
1896   return standardDeviation;
1897 }
1898
1899 //____________________________________________________________________
1900 TH1* AliMultiplicityCorrection::StatisticalUncertainty(MethodType methodType, Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t randomizeMeasured, Bool_t randomizeResponse, TH1* compareTo)
1901 {
1902   //
1903   // evaluates the uncertainty that arises from the non-infinite statistics in the response matrix
1904   // the function unfolds the spectrum using the default response matrix and several modified ones
1905   // the modified ones are created by randomizing each cell using poisson statistics with the mean = bin value
1906   // these unfolded results are compared to the first result gained with the default response OR to the histogram given
1907   // in <compareTo> (optional)
1908   //
1909   // returns the error assigned to the measurement
1910   //
1911
1912   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
1913
1914   // initialize seed with current time
1915   gRandom->SetSeed(0);
1916
1917   const Int_t kErrorIterations = 150;
1918
1919   TH1* maxError = 0;
1920   TH1* firstResult = 0;
1921
1922   TH1** results = new TH1*[kErrorIterations];
1923
1924   for (Int_t n=0; n<kErrorIterations; ++n)
1925   {
1926     Printf("Iteration %d of %d...", n, kErrorIterations);
1927
1928     // only done for chi2 minimization
1929     Bool_t createBigBin = kFALSE;
1930     if (methodType == kChi2Minimization)
1931       createBigBin = kTRUE;
1932
1933     SetupCurrentHists(inputRange, fullPhaseSpace, eventType, createBigBin);
1934
1935     TH1* measured = (TH1*) fCurrentESD->Clone("measured");
1936
1937     if (n > 0)
1938     {
1939       if (randomizeResponse)
1940       {
1941         // randomize response matrix
1942         for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
1943           for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
1944             fCurrentCorrelation->SetBinContent(i, j, gRandom->Poisson(fCurrentCorrelation->GetBinContent(i, j)));
1945       }
1946
1947       if (randomizeMeasured)
1948       {
1949         // randomize measured spectrum
1950         for (Int_t x=1; x<=measured->GetNbinsX(); x++) // mult. axis
1951         {
1952           Int_t randomValue = gRandom->Poisson(fCurrentESD->GetBinContent(x));
1953           measured->SetBinContent(x, randomValue);
1954           measured->SetBinError(x, TMath::Sqrt(randomValue));
1955         }
1956       }
1957     }
1958
1959     // only for bayesian method we have to do it before the call to Unfold...
1960     if (methodType == kBayesian)
1961     {
1962       for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
1963       {
1964         // with this it is normalized to 1
1965         Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
1966
1967         // with this normalized to the given efficiency
1968         if (fCurrentEfficiency->GetBinContent(i) > 0)
1969           sum /= fCurrentEfficiency->GetBinContent(i);
1970         else
1971           sum = 0;
1972
1973         for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
1974         {
1975           if (sum > 0)
1976           {
1977             fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
1978             fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
1979           }
1980           else
1981           {
1982             fCurrentCorrelation->SetBinContent(i, j, 0);
1983             fCurrentCorrelation->SetBinError(i, j, 0);
1984           }
1985         }
1986       }
1987     }
1988
1989     TH1* result = 0;
1990     if (n == 0 && compareTo)
1991     {
1992       // in this case we just store the histogram we want to compare to
1993       result = (TH1*) compareTo->Clone("compareTo");
1994       result->Sumw2();
1995     }
1996     else
1997     {
1998       result = (TH1*) fMultiplicityESDCorrected[correlationID]->Clone(Form("result_%d", n));
1999
2000       Int_t returnCode = -1;
2001       if (methodType == kChi2Minimization)
2002       {
2003         returnCode = UnfoldWithMinuit(fCurrentCorrelation, fCurrentEfficiency, measured, 0, result, kFALSE);
2004       }
2005       else if (methodType == kBayesian)
2006         returnCode = UnfoldWithBayesian(fCurrentCorrelation, fCurrentEfficiency, measured, 0, result, fgBayesianSmoothing, fgBayesianIterations);
2007
2008       if (returnCode != 0)
2009         return 0;
2010     }
2011
2012     // normalize
2013     result->Scale(1.0 / result->Integral());
2014
2015     if (n == 0)
2016     {
2017       firstResult = (TH1*) result->Clone("firstResult");
2018
2019       maxError = (TH1*) result->Clone("maxError");
2020       maxError->Reset();
2021     }
2022     else
2023     {
2024       // calculate ratio
2025       TH1* ratio = (TH1*) firstResult->Clone("ratio");
2026       ratio->Divide(result);
2027
2028       // find max. deviation
2029       for (Int_t x=1; x<=ratio->GetNbinsX(); x++)
2030         maxError->SetBinContent(x, TMath::Max(maxError->GetBinContent(x), TMath::Abs(1 - ratio->GetBinContent(x))));
2031
2032       delete ratio;
2033     }
2034
2035     results[n] = result;
2036   }
2037
2038   // find covariance matrix
2039   // results[n] is X_x
2040   // cov. matrix is M_xy = E ( (X_x - E(X_x)) * (X_y - E(X_y))), with E() = expectation value
2041
2042   Int_t nBins = results[0]->GetNbinsX();
2043   Float_t lowEdge = results[0]->GetXaxis()->GetBinLowEdge(1);
2044   Float_t upEdge = results[0]->GetXaxis()->GetBinUpEdge(nBins);
2045
2046   // find average, E(X)
2047   TProfile* average = new TProfile("average", "average", nBins, lowEdge, upEdge);
2048   for (Int_t n=1; n<kErrorIterations; ++n)
2049     for (Int_t x=1; x<=results[n]->GetNbinsX(); x++)
2050       average->Fill(results[n]->GetXaxis()->GetBinCenter(x), results[n]->GetBinContent(x));
2051   //new TCanvas; average->DrawClone();
2052
2053   // find cov. matrix
2054   TProfile2D* covMatrix = new TProfile2D("covMatrix", "covMatrix", nBins, lowEdge, upEdge, nBins, lowEdge, upEdge);
2055
2056   for (Int_t n=1; n<kErrorIterations; ++n)
2057     for (Int_t x=1; x<=results[n]->GetNbinsX(); x++)
2058       for (Int_t y=1; y<=results[n]->GetNbinsX(); y++)
2059       {
2060         // (X_x - E(X_x)) * (X_y - E(X_y)
2061         Float_t cov = (results[n]->GetBinContent(x) - average->GetBinContent(x)) * (results[n]->GetBinContent(y) - average->GetBinContent(y));
2062         covMatrix->Fill(results[n]->GetXaxis()->GetBinCenter(x), results[n]->GetXaxis()->GetBinCenter(y), cov);
2063       }
2064   new TCanvas; covMatrix->DrawClone("COLZ");
2065
2066 //   // fill 2D histogram that contains deviation from first
2067 //   TH2F* deviations = new TH2F("deviations", "deviations", nBins, lowEdge, upEdge, 1000, -0.01, 0.01);
2068 //   for (Int_t n=1; n<kErrorIterations; ++n)
2069 //     for (Int_t x=1; x<=results[n]->GetNbinsX(); x++)
2070 //       deviations->Fill(results[n]->GetXaxis()->GetBinCenter(x), results[n]->GetBinContent(x) - results[0]->GetBinContent(x));
2071 //   //new TCanvas; deviations->DrawCopy("COLZ");
2072 //
2073 //   // get standard deviation "by hand"
2074 //   for (Int_t x=1; x<=nBins; x++)
2075 //   {
2076 //     if (results[0]->GetBinContent(x) > 0)
2077 //     {
2078 //       TH1* proj = deviations->ProjectionY("projRMS", x, x, "e");
2079 //       Float_t standardDev = proj->GetRMS(); // this is standard deviation in fact
2080 //       //standardDeviation->SetBinContent(x, standardDev / results[0]->GetBinContent(x));
2081 //       Printf("sigma_%d is %f value %f --> error %f", x, standardDev, results[0]->GetBinContent(x), standardDev / results[0]->GetBinContent(x));
2082 //     }
2083 //   }
2084
2085   TH1* standardDeviation = CalculateStdDev(results, kErrorIterations);
2086
2087   // compare maxError to RMS of profile (variable name: average)
2088   // first: calculate rms in percent of value
2089   TH1* rmsError = (TH1*) maxError->Clone("rmsError");
2090   rmsError->Reset();
2091
2092   // enable error to be standard deviation (see http://root.cern.ch/root/html/TProfile.html#TProfile:SetErrorOption)
2093   average->SetErrorOption("s");
2094   for (Int_t x=1; x<=rmsError->GetNbinsX(); x++)
2095     if (average->GetBinContent(x) > 0)
2096       rmsError->SetBinContent(x, average->GetBinError(x) / average->GetBinContent(x));
2097
2098   // find maxError deviation from average (not from "first result"/mc as above)
2099   TH1* maxError2 = (TH1*) maxError->Clone("maxError2");
2100   maxError2->Reset();
2101   for (Int_t n=1; n<kErrorIterations; ++n)
2102     for (Int_t x=1; x<=results[n]->GetNbinsX(); x++)
2103       if (average->GetBinContent(x) > 0)
2104         maxError2->SetBinContent(x, TMath::Max(maxError2->GetBinContent(x), TMath::Abs((results[n]->GetBinContent(x) - average->GetBinContent(x)) / average->GetBinContent(x))));
2105
2106   //new TCanvas; maxError2->DrawCopy(); rmsError->SetLineColor(2); rmsError->DrawCopy("SAME"); standardDeviation->SetLineColor(3); standardDeviation->DrawCopy("SAME");
2107
2108   // plot difference between average and MC/first result
2109   TH1* averageFirstRatio = (TH1*) results[0]->Clone("averageFirstRatio");
2110   averageFirstRatio->Reset();
2111   averageFirstRatio->Divide(results[0], average);
2112
2113   //new TCanvas; results[0]->DrawCopy(); average->SetLineColor(2); average->DrawClone("SAME");
2114   //new TCanvas; averageFirstRatio->DrawCopy();
2115
2116   // clean up
2117   for (Int_t n=0; n<kErrorIterations; ++n)
2118     delete results[n];
2119   delete[] results;
2120
2121   // fill into result histogram
2122   for (Int_t i=1; i<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); ++i)
2123     fMultiplicityESDCorrected[correlationID]->SetBinContent(i, firstResult->GetBinContent(i));
2124
2125   for (Int_t i=1; i<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); ++i)
2126     fMultiplicityESDCorrected[correlationID]->SetBinError(i, maxError->GetBinContent(i) * fMultiplicityESDCorrected[correlationID]->GetBinContent(i));
2127
2128   return standardDeviation;
2129 }
2130
2131 //____________________________________________________________________
2132 void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Float_t regPar, Int_t nIterations, TH1* initialConditions, Bool_t determineError)
2133 {
2134   //
2135   // correct spectrum using bayesian method
2136   //
2137
2138   // initialize seed with current time
2139   gRandom->SetSeed(0);
2140
2141   SetupCurrentHists(inputRange, fullPhaseSpace, eventType, kFALSE);
2142
2143   // normalize correction for given nPart
2144   for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
2145   {
2146     // with this it is normalized to 1
2147     Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
2148
2149     // with this normalized to the given efficiency
2150     if (fCurrentEfficiency->GetBinContent(i) > 0)
2151       sum /= fCurrentEfficiency->GetBinContent(i);
2152     else
2153       sum = 0;
2154
2155     for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
2156     {
2157       if (sum > 0)
2158       {
2159         fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
2160         fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
2161       }
2162       else
2163       {
2164         fCurrentCorrelation->SetBinContent(i, j, 0);
2165         fCurrentCorrelation->SetBinError(i, j, 0);
2166       }
2167     }
2168   }
2169
2170   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
2171
2172   if (UnfoldWithBayesian(fCurrentCorrelation, fCurrentEfficiency, fCurrentESD, initialConditions, fMultiplicityESDCorrected[correlationID], regPar, nIterations) != 0)
2173     return;
2174
2175   if (!determineError)
2176   {
2177     Printf("AliMultiplicityCorrection::ApplyBayesianMethod: WARNING: No errors calculated.");
2178     return;
2179   }
2180
2181   // evaluate errors, this is done by randomizing the measured spectrum following Poission statistics
2182   // this (new) measured spectrum is then unfolded and the different to the result from the "real" measured
2183   // spectrum calculated. This is performed N times and the sigma is taken as the statistical
2184   // error of the unfolding method itself.
2185
2186   const Int_t kErrorIterations = 20;
2187
2188   printf("Spectrum unfolded. Determining error (%d iterations)...\n", kErrorIterations);
2189
2190   TH1* randomized = (TH1*) fCurrentESD->Clone("randomized");
2191   TH1* resultArray[kErrorIterations+1];
2192   for (Int_t n=0; n<kErrorIterations; ++n)
2193   {
2194     // randomize the content of clone following a poisson with the mean = the value of that bin
2195     for (Int_t x=1; x<=randomized->GetNbinsX(); x++) // mult. axis
2196     {
2197       Int_t randomValue = gRandom->Poisson(fCurrentESD->GetBinContent(x));
2198       //printf("%d --> %d\n", fCurrentESD->GetBinContent(x), randomValue);
2199       randomized->SetBinContent(x, randomValue);
2200       randomized->SetBinError(x, TMath::Sqrt(randomValue));
2201     }
2202
2203     TH1* result2 = (TH1*) fMultiplicityESDCorrected[correlationID]->Clone("result2");
2204     result2->Reset();
2205     if (UnfoldWithBayesian(fCurrentCorrelation, fCurrentEfficiency, randomized, initialConditions, result2, regPar, nIterations) != 0)
2206       return;
2207
2208     result2->Scale(1.0 / result2->Integral());
2209
2210     resultArray[n+1] = result2;
2211   }
2212   delete randomized;
2213
2214   resultArray[0] = fMultiplicityESDCorrected[correlationID];
2215   TH1* error = CalculateStdDev(resultArray, kErrorIterations+1);
2216
2217   for (Int_t n=0; n<kErrorIterations; ++n)
2218     delete resultArray[n+1];
2219
2220   for (Int_t i=1; i<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); ++i)
2221     fMultiplicityESDCorrected[correlationID]->SetBinError(i, error->GetBinContent(i) * fMultiplicityESDCorrected[correlationID]->GetBinContent(i));
2222
2223   delete error;
2224 }
2225
2226 //____________________________________________________________________
2227 Int_t AliMultiplicityCorrection::UnfoldWithBayesian(TH1* correlation, TH1* aEfficiency, TH1* measured, TH1* initialConditions, TH1* aResult, Float_t regPar, Int_t nIterations)
2228 {
2229   //
2230   // unfolds a spectrum
2231   //
2232
2233   if (measured->Integral() <= 0)
2234   {
2235     Printf("AliMultiplicityCorrection::UnfoldWithBayesian: ERROR: The measured spectrum is empty");
2236     return -1;
2237   }
2238
2239   const Int_t kStartBin = 0;
2240
2241   const Int_t kMaxM = fgkMaxInput;  //<= fCurrentCorrelation->GetNbinsY(); // max measured axis
2242   const Int_t kMaxT = fgkMaxParams; //<= fCurrentCorrelation->GetNbinsX(); // max true axis
2243
2244   // convergence limit: kMaxT * 0.001^2 = kMaxT * 1e-6 (e.g. 250 bins --> 2.5 e-4)
2245   const Double_t kConvergenceLimit = kMaxT * 1e-6;
2246
2247   // store information in arrays, to increase processing speed (~ factor 5)
2248   Double_t measuredCopy[kMaxM];
2249   Double_t measuredError[kMaxM];
2250   Double_t prior[kMaxT];
2251   Double_t result[kMaxT];
2252   Double_t efficiency[kMaxT];
2253
2254   Double_t response[kMaxT][kMaxM];
2255   Double_t inverseResponse[kMaxT][kMaxM];
2256
2257   // for normalization
2258   Float_t measuredIntegral = measured->Integral();
2259   for (Int_t m=0; m<kMaxM; m++)
2260   {
2261     measuredCopy[m] = measured->GetBinContent(m+1) / measuredIntegral;
2262     measuredError[m] = measured->GetBinError(m+1) / measuredIntegral;
2263
2264     for (Int_t t=0; t<kMaxT; t++)
2265     {
2266       response[t][m] = correlation->GetBinContent(t+1, m+1);
2267       inverseResponse[t][m] = 0;
2268     }
2269   }
2270
2271   for (Int_t t=0; t<kMaxT; t++)
2272   {
2273     efficiency[t] = aEfficiency->GetBinContent(t+1);
2274     prior[t] = measuredCopy[t];
2275     result[t] = 0;
2276   }
2277
2278   // pick prior distribution
2279   if (initialConditions)
2280   {
2281     printf("Using different starting conditions...\n");
2282     // for normalization
2283     Float_t inputDistIntegral = initialConditions->Integral();
2284     for (Int_t i=0; i<kMaxT; i++)
2285       prior[i] = initialConditions->GetBinContent(i+1) / inputDistIntegral;
2286   }
2287
2288   //TH1F* convergence = new TH1F("convergence", "convergence", 200, 0.5, 200.5);
2289
2290   // unfold...
2291   for (Int_t i=0; i<nIterations || nIterations < 0; i++)
2292   {
2293     //printf(" iteration %i \n", i);
2294
2295     // calculate IR from Beyes theorem
2296     // IR_ji = R_ij * prior_i / sum_k(R_kj * prior_k)
2297
2298     Double_t chi2Measured = 0;
2299     for (Int_t m=0; m<kMaxM; m++)
2300     {
2301       Float_t norm = 0;
2302       for (Int_t t = kStartBin; t<kMaxT; t++)
2303         norm += response[t][m] * prior[t];
2304
2305       // calc. chi2: (measured - response * prior) / error
2306       if (measuredError[m] > 0)
2307       {
2308         Double_t value = (measuredCopy[m] - norm) / measuredError[m];
2309         chi2Measured += value * value;
2310       }
2311
2312       if (norm > 0)
2313       {
2314         for (Int_t t = kStartBin; t<kMaxT; t++)
2315           inverseResponse[t][m] = response[t][m] * prior[t] / norm;
2316       }
2317       else
2318       {
2319         for (Int_t t = kStartBin; t<kMaxT; t++)
2320           inverseResponse[t][m] = 0;
2321       }
2322     }
2323     //Printf("chi2Measured of the last prior is %e", chi2Measured);
2324
2325     for (Int_t t = kStartBin; t<kMaxT; t++)
2326     {
2327       Float_t value = 0;
2328       for (Int_t m=0; m<kMaxM; m++)
2329         value += inverseResponse[t][m] * measuredCopy[m];
2330
2331       if (efficiency[t] > 0)
2332         result[t] = value / efficiency[t];
2333       else
2334         result[t] = 0;
2335     }
2336
2337     Double_t chi2LastIter = 0;
2338     // regularization (simple smoothing)
2339     for (Int_t t=kStartBin; t<kMaxT; t++)
2340     {
2341       Float_t newValue = 0;
2342       // 0 bin excluded from smoothing
2343       if (t > kStartBin+1 && t<kMaxT-1)
2344       {
2345         Float_t average = (result[t-1] + result[t] + result[t+1]) / 3;
2346
2347         // weight the average with the regularization parameter
2348         newValue = (1 - regPar) * result[t] + regPar * average;
2349       }
2350       else
2351         newValue = result[t];
2352
2353       // calculate chi2 (change from last iteration)
2354       if (prior[t] > 1e-5)
2355       {
2356         Double_t diff = (prior[t] - newValue) / prior[t];
2357         chi2LastIter += diff * diff;
2358       }
2359
2360       prior[t] = newValue;
2361     }
2362     //printf("Chi2 of %d iteration = %e\n", i, chi2LastIter);
2363     //convergence->Fill(i+1, chi2LastIter);
2364
2365     if (nIterations < 0 && chi2LastIter < kConvergenceLimit)
2366     {
2367       Printf("Stopped Bayesian unfolding after %d iterations at chi2(change since last iteration) of %e; chi2Measured of the last prior is %e", i, chi2LastIter, chi2Measured);
2368       break;
2369     }
2370
2371     //if (i % 5 == 0)
2372     //  DrawComparison(Form("Bayesian_%d", i), mcTarget, correlationID, kTRUE, eventType);
2373   } // end of iterations
2374
2375   //new TCanvas; convergence->DrawCopy(); gPad->SetLogy();
2376   //delete convergence;
2377
2378   for (Int_t t=0; t<kMaxT; t++)
2379     aResult->SetBinContent(t+1, result[t]);
2380
2381   return 0;
2382
2383   // ********
2384   // Calculate the covariance matrix, all arguments are taken from NIM,A362,487-498,1995
2385
2386   /*printf("Calculating covariance matrix. This may take some time...\n");
2387
2388   // TODO check if this is the right one...
2389   TH1* sumHist = GetMultiplicityMC(inputRange, eventType)->ProjectionY("sumHist", 1, GetMultiplicityMC(inputRange, eventType)->GetNbinsX());
2390
2391   Int_t xBins = hInverseResponseBayes->GetNbinsX();
2392   Int_t yBins = hInverseResponseBayes->GetNbinsY();
2393
2394   // calculate "unfolding matrix" Mij
2395   Float_t matrixM[251][251];
2396   for (Int_t i=1; i<=xBins; i++)
2397   {
2398     for (Int_t j=1; j<=yBins; j++)
2399     {
2400       if (fCurrentEfficiency->GetBinContent(i) > 0)
2401         matrixM[i-1][j-1] = hInverseResponseBayes->GetBinContent(i, j) / fCurrentEfficiency->GetBinContent(i);
2402       else
2403         matrixM[i-1][j-1] = 0;
2404     }
2405   }
2406
2407   Float_t* vectorn = new Float_t[yBins];
2408   for (Int_t j=1; j<=yBins; j++)
2409     vectorn[j-1] = fCurrentESD->GetBinContent(j);
2410
2411   // first part of covariance matrix, depends on input distribution n(E)
2412   Float_t cov1[251][251];
2413
2414   Float_t nEvents = fCurrentESD->Integral(); // N
2415
2416   xBins = 20;
2417   yBins = 20;
2418
2419   for (Int_t k=0; k<xBins; k++)
2420   {
2421     printf("In Cov1: %d\n", k);
2422     for (Int_t l=0; l<yBins; l++)
2423     {
2424       cov1[k][l] = 0;
2425
2426       // sum_j Mkj Mlj n(Ej) * (1 - n(Ej) / N)
2427       for (Int_t j=0; j<yBins; j++)
2428         cov1[k][l] += matrixM[k][j] * matrixM[l][j] * vectorn[j]
2429           * (1.0 - vectorn[j] / nEvents);
2430
2431       // - sum_i,j (i != j) Mki Mlj n(Ei) n(Ej) / N
2432       for (Int_t i=0; i<yBins; i++)
2433         for (Int_t j=0; j<yBins; j++)
2434         {
2435           if (i == j)
2436             continue;
2437           cov1[k][l] -= matrixM[k][i] * matrixM[l][j] * vectorn[i]
2438             * vectorn[j] / nEvents;
2439          }
2440     }
2441   }
2442
2443   printf("Cov1 finished\n");
2444
2445   TH2F* cov = (TH2F*) hInverseResponseBayes->Clone("cov");
2446   cov->Reset();
2447
2448   for (Int_t i=1; i<=xBins; i++)
2449     for (Int_t j=1; j<=yBins; j++)
2450       cov->SetBinContent(i, j, cov1[i-1][j-1]);
2451
2452   new TCanvas;
2453   cov->Draw("COLZ");
2454
2455   // second part of covariance matrix, depends on response matrix
2456   Float_t cov2[251][251];
2457
2458   // Cov[P(Er|Cu), P(Es|Cu)] term
2459   Float_t covTerm[100][100][100];
2460   for (Int_t r=0; r<yBins; r++)
2461     for (Int_t u=0; u<xBins; u++)
2462       for (Int_t s=0; s<yBins; s++)
2463       {
2464         if (r == s)
2465           covTerm[r][u][s] = 1.0 / sumHist->GetBinContent(u+1) * hResponse->GetBinContent(u+1, r+1)
2466             * (1.0 - hResponse->GetBinContent(u+1, r+1));
2467         else
2468           covTerm[r][u][s] = - 1.0 / sumHist->GetBinContent(u+1) * hResponse->GetBinContent(u+1, r+1)
2469             * hResponse->GetBinContent(u+1, s+1);
2470       }
2471
2472   for (Int_t k=0; k<xBins; k++)
2473     for (Int_t l=0; l<yBins; l++)
2474     {
2475       cov2[k][l] = 0;
2476       printf("In Cov2: %d %d\n", k, l);
2477       for (Int_t i=0; i<yBins; i++)
2478         for (Int_t j=0; j<yBins; j++)
2479         {
2480           //printf("In Cov2: %d %d %d %d\n", k, l, i, j);
2481           // calculate Cov(Mki, Mlj) = sum{ru},{su} ...
2482           Float_t tmpCov = 0;
2483           for (Int_t r=0; r<yBins; r++)
2484             for (Int_t u=0; u<xBins; u++)
2485               for (Int_t s=0; s<yBins; s++)
2486               {
2487                 if (hResponse->GetBinContent(u+1, r+1) == 0 || hResponse->GetBinContent(u+1, s+1) == 0
2488                   || hResponse->GetBinContent(u+1, i+1) == 0)
2489                   continue;
2490
2491                 tmpCov += BayesCovarianceDerivate(matrixM, hResponse, fCurrentEfficiency, k, i, r, u)
2492                   * BayesCovarianceDerivate(matrixM, hResponse, fCurrentEfficiency, l, j, s, u)
2493                   * covTerm[r][u][s];
2494               }
2495
2496           cov2[k][l] += fCurrentESD->GetBinContent(i+1) * fCurrentESD->GetBinContent(j+1) * tmpCov;
2497         }
2498     }
2499
2500   printf("Cov2 finished\n");
2501
2502   for (Int_t i=1; i<=xBins; i++)
2503     for (Int_t j=1; j<=yBins; j++)
2504       cov->SetBinContent(i, j, cov1[i-1][j-1] + cov2[i-1][j-1]);
2505
2506   new TCanvas;
2507   cov->Draw("COLZ");*/
2508 }
2509
2510 //____________________________________________________________________
2511 Float_t AliMultiplicityCorrection::BayesCovarianceDerivate(Float_t matrixM[251][251], TH2* hResponse, Int_t k, Int_t i, Int_t r, Int_t u)
2512 {
2513   //
2514   // helper function for the covariance matrix of the bayesian method
2515   //
2516
2517   Float_t result = 0;
2518
2519   if (k == u && r == i)
2520     result += 1.0 / hResponse->GetBinContent(u+1, r+1);
2521
2522   if (k == u)
2523     result -= 1.0 / fCurrentEfficiency->GetBinContent(u+1);
2524
2525   if (r == i)
2526     result -= matrixM[u][i] * fCurrentEfficiency->GetBinContent(u+1) / hResponse->GetBinContent(u+1, i+1);
2527
2528   result *= matrixM[k][i];
2529
2530   return result;
2531 }
2532
2533 //____________________________________________________________________
2534 void AliMultiplicityCorrection::ApplyLaszloMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType)
2535 {
2536   //
2537   // correct spectrum using bayesian method
2538   //
2539
2540   Float_t regPar = 0;
2541
2542   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
2543   Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4);
2544
2545   SetupCurrentHists(inputRange, fullPhaseSpace, eventType, kFALSE);
2546   //normalize ESD
2547   fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
2548
2549   // TODO should be taken from correlation map
2550   //TH1* sumHist = GetMultiplicityMC(inputRange, eventType)->ProjectionY("sumHist", 1, GetMultiplicityMC(inputRange, eventType)->GetNbinsX());
2551
2552   // normalize correction for given nPart
2553   for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
2554   {
2555     Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
2556     //Double_t sum = sumHist->GetBinContent(i);
2557     if (sum <= 0)
2558       continue;
2559     for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
2560     {
2561       // npart sum to 1
2562       fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
2563       fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
2564     }
2565   }
2566
2567   new TCanvas;
2568   fCurrentCorrelation->Draw("COLZ");
2569
2570   // FAKE
2571   fCurrentEfficiency = ((TH2*) fCurrentCorrelation)->ProjectionX("eff");
2572
2573   // pick prior distribution
2574   TH1F* hPrior = (TH1F*)fCurrentESD->Clone("prior");
2575   Float_t norm = 1; //hPrior->Integral();
2576   for (Int_t t=1; t<=hPrior->GetNbinsX(); t++)
2577     hPrior->SetBinContent(t, hPrior->GetBinContent(t)/norm);
2578
2579   // zero distribution
2580   TH1F* zero =  (TH1F*)hPrior->Clone("zero");
2581
2582   // define temp hist
2583   TH1F* hTemp = (TH1F*)fCurrentESD->Clone("temp");
2584   hTemp->Reset();
2585
2586   // just a shortcut
2587   TH2F* hResponse = (TH2F*) fCurrentCorrelation;
2588
2589   // unfold...
2590   Int_t iterations = 25;
2591   for (Int_t i=0; i<iterations; i++)
2592   {
2593     //printf(" iteration %i \n", i);
2594
2595     for (Int_t m=1; m<=hResponse->GetNbinsY(); m++)
2596     {
2597       Float_t value = 0;
2598       for (Int_t t = 1; t<=hResponse->GetNbinsX(); t++)
2599         value += hResponse->GetBinContent(t, m) * hPrior->GetBinContent(t);
2600       hTemp->SetBinContent(m, value);
2601       //printf("%d %f %f %f\n", m, zero->GetBinContent(m), hPrior->GetBinContent(m), value);
2602     }
2603
2604     // regularization (simple smoothing)
2605     TH1F* hTrueSmoothed = (TH1F*) hTemp->Clone("truesmoothed");
2606
2607     for (Int_t t=2; t<hTrueSmoothed->GetNbinsX(); t++)
2608     {
2609       Float_t average = (hTemp->GetBinContent(t-1) / hTemp->GetBinWidth(t-1)
2610                          + hTemp->GetBinContent(t) / hTemp->GetBinWidth(t)
2611                          + hTemp->GetBinContent(t+1) / hTemp->GetBinWidth(t+1)) / 3.;
2612       average *= hTrueSmoothed->GetBinWidth(t);
2613
2614       // weight the average with the regularization parameter
2615       hTrueSmoothed->SetBinContent(t, (1-regPar) * hTemp->GetBinContent(t) + regPar * average);
2616     }
2617
2618     for (Int_t m=1; m<=hResponse->GetNbinsY(); m++)
2619       hTemp->SetBinContent(m, zero->GetBinContent(m) + hPrior->GetBinContent(m) - hTrueSmoothed->GetBinContent(m));
2620
2621     // fill guess
2622     for (Int_t t=1; t<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); t++)
2623     {
2624       fMultiplicityESDCorrected[correlationID]->SetBinContent(t, hTemp->GetBinContent(t));
2625       fMultiplicityESDCorrected[correlationID]->SetBinError(t, 0.05 * hTemp->GetBinContent(t)); // TODO
2626
2627       //printf(" bin %d content %f \n", t, fMultiplicityESDCorrected[correlationID]->GetBinContent(t));
2628     }
2629
2630
2631     // calculate chi2 (change from last iteration)
2632     Double_t chi2 = 0;
2633
2634     // use smoothed true (normalized) as new prior
2635     norm = 1; //hTrueSmoothed->Integral();
2636
2637     for (Int_t t=1; t<hTrueSmoothed->GetNbinsX(); t++)
2638     {
2639       Float_t newValue = hTemp->GetBinContent(t)/norm;
2640       Float_t diff = hPrior->GetBinContent(t) - newValue;
2641       chi2 += (Double_t) diff * diff;
2642
2643       hPrior->SetBinContent(t, newValue);
2644     }
2645
2646     printf("Chi2 of %d iteration = %.10f\n", i, chi2);
2647
2648     //if (i % 5 == 0)
2649       DrawComparison(Form("Laszlo_%d", i), inputRange, fullPhaseSpace, kTRUE, GetMultiplicityMC(mcTarget, eventType)->ProjectionY());
2650
2651     delete hTrueSmoothed;
2652   } // end of iterations
2653
2654   DrawComparison("Laszlo", inputRange, fullPhaseSpace, kTRUE, GetMultiplicityMC(mcTarget, eventType)->ProjectionY());
2655 }
2656
2657 //____________________________________________________________________
2658 void AliMultiplicityCorrection::ApplyGaussianMethod(Int_t inputRange, Bool_t fullPhaseSpace)
2659 {
2660   //
2661   // correct spectrum using a simple Gaussian approach, that is model-dependent
2662   //
2663
2664   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
2665   Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4);
2666
2667   SetupCurrentHists(inputRange, fullPhaseSpace, kTrVtx, kFALSE);
2668   //normalize ESD
2669   fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
2670
2671   TH1D* correction = dynamic_cast<TH1D*> (fCurrentESD->Clone("GaussianMean"));
2672   correction->SetTitle("GaussianMean");
2673
2674   TH1D* correctionWidth = dynamic_cast<TH1D*> (fCurrentESD->Clone("GaussianWidth"));
2675   correction->SetTitle("GaussianWidth");
2676
2677   for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
2678   {
2679     TH1D* proj = (dynamic_cast<TH2*> (fCurrentCorrelation))->ProjectionX("_px", i, i+1);
2680     proj->Fit("gaus", "0Q");
2681     correction->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(1));
2682     correctionWidth->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(2));
2683
2684     continue;
2685
2686     // draw for debugging
2687     new TCanvas;
2688     proj->DrawCopy();
2689     proj->GetFunction("gaus")->DrawCopy("SAME");
2690   }
2691
2692   TH1* target = fMultiplicityESDCorrected[correlationID];
2693
2694   Int_t nBins = target->GetNbinsX()*10+1;
2695   Float_t* binning = new Float_t[nBins];
2696   for (Int_t i=1; i<=target->GetNbinsX(); ++i)
2697     for (Int_t j=0; j<10; ++j)
2698       binning[(i-1)*10 + j] = target->GetXaxis()->GetBinLowEdge(i) + target->GetXaxis()->GetBinWidth(i) / 10 * j;
2699
2700   binning[nBins-1] = target->GetXaxis()->GetBinUpEdge(target->GetNbinsX());
2701
2702   TH1F* fineBinned = new TH1F("targetFineBinned", "targetFineBinned", nBins-1, binning);
2703
2704   for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
2705   {
2706     Float_t mean = correction->GetBinContent(i);
2707     Float_t width = correctionWidth->GetBinContent(i);
2708
2709     Int_t fillBegin = fineBinned->FindBin(mean - width * 5);
2710     Int_t fillEnd   = fineBinned->FindBin(mean + width * 5);
2711     //printf("bin %d mean %f width %f, filling from %d to %d\n", i, mean, width, fillBegin, fillEnd);
2712
2713     for (Int_t j=fillBegin; j <= fillEnd; ++j)
2714     {
2715       fineBinned->AddBinContent(j, TMath::Gaus(fineBinned->GetXaxis()->GetBinCenter(j), mean, width, kTRUE) * fCurrentESD->GetBinContent(i));
2716     }
2717   }
2718
2719   for (Int_t i=1; i<=target->GetNbinsX(); ++i)
2720   {
2721     Float_t sum = 0;
2722     for (Int_t j=1; j<=10; ++j)
2723       sum += fineBinned->GetBinContent((i-1)*10 + j);
2724     target->SetBinContent(i, sum / 10);
2725   }
2726
2727   delete[] binning;
2728
2729   DrawComparison("Gaussian", inputRange, fullPhaseSpace, kFALSE, GetMultiplicityMC(mcTarget, kTrVtx)->ProjectionY());
2730 }
2731
2732 //____________________________________________________________________
2733 TH2F* AliMultiplicityCorrection::CalculateMultiplicityESD(TH1* inputMC, Int_t correlationMap)
2734 {
2735   // runs the distribution given in inputMC through the response matrix identified by
2736   // correlationMap and produces a measured distribution
2737   // although it is a TH2F the vertex axis is not used at the moment and all entries are filled in mid-vertex
2738   // if normalized is set, inputMC is expected to be normalized to the bin width
2739
2740   if (!inputMC)
2741     return 0;
2742
2743   if (correlationMap < 0 || correlationMap >= kCorrHists)
2744     return 0;
2745
2746   // empty under/overflow bins in x, otherwise Project3D takes them into account
2747   TH3* hist = fCorrelation[correlationMap];
2748   for (Int_t y=0; y<=hist->GetYaxis()->GetNbins()+1; ++y)
2749   {
2750     for (Int_t z=0; z<=hist->GetZaxis()->GetNbins()+1; ++z)
2751     {
2752       hist->SetBinContent(0, y, z, 0);
2753       hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0);
2754     }
2755   }
2756
2757   TH2* corr = (TH2*) hist->Project3D("zy");
2758   //corr->Rebin2D(2, 1);
2759   corr->Sumw2();
2760   Printf("Correction histogram used for convolution has %f entries", corr->Integral());
2761
2762   // normalize correction for given nPart
2763   for (Int_t i=1; i<=corr->GetNbinsX(); ++i)
2764   {
2765     Double_t sum = corr->Integral(i, i, 1, corr->GetNbinsY());
2766     if (sum <= 0)
2767       continue;
2768
2769     for (Int_t j=1; j<=corr->GetNbinsY(); ++j)
2770     {
2771       // npart sum to 1
2772       corr->SetBinContent(i, j, corr->GetBinContent(i, j) / sum);
2773       corr->SetBinError(i, j, corr->GetBinError(i, j) / sum);
2774     }
2775   }
2776
2777   TH2F* target = dynamic_cast<TH2F*> (fMultiplicityESD[0]->Clone(Form("%s_measured", inputMC->GetName())));
2778   target->Reset();
2779
2780   for (Int_t meas=1; meas<=corr->GetNbinsY(); ++meas)
2781   {
2782     Float_t measured = 0;
2783     Float_t error = 0;
2784
2785     for (Int_t gen=1; gen<=corr->GetNbinsX(); ++gen)
2786     {
2787       Int_t mcGenBin = inputMC->GetXaxis()->FindBin(corr->GetXaxis()->GetBinCenter(gen));
2788
2789       measured += inputMC->GetBinContent(mcGenBin) * corr->GetBinContent(gen, meas);
2790       error += inputMC->GetBinError(mcGenBin) * corr->GetBinContent(gen, meas);
2791     }
2792
2793     //printf("%f +- %f ; %f +- %f \n", inputMC->GetBinContent(meas), inputMC->GetBinError(meas), measured, error);
2794
2795     target->SetBinContent(1 + target->GetNbinsX() / 2, meas, measured);
2796     target->SetBinError(1 + target->GetNbinsX() / 2, meas, error);
2797   }
2798
2799   return target;
2800 }
2801
2802 //____________________________________________________________________
2803 void AliMultiplicityCorrection::SetGenMeasFromFunc(TF1* inputMC, Int_t id)
2804 {
2805   // uses the given function to fill the input MC histogram and generates from that
2806   // the measured histogram by applying the response matrix
2807   // this can be used to evaluate if the methods work indepedently of the input
2808   // distribution
2809   // WARNING does not respect the vertex distribution, just fills central vertex bin
2810
2811   if (!inputMC)
2812     return;
2813
2814   if (id < 0 || id >= kESDHists)
2815     return;
2816
2817   TH2F* mc = fMultiplicityVtx[id];
2818
2819   mc->Reset();
2820
2821   for (Int_t i=1; i<=mc->GetNbinsY(); ++i)
2822   {
2823     mc->SetBinContent(mc->GetNbinsX() / 2 + 1, i, inputMC->Eval(mc->GetYaxis()->GetBinCenter(i)) * mc->GetYaxis()->GetBinWidth(i));
2824     mc->SetBinError(mc->GetNbinsX() / 2 + 1, i, 0);
2825   }
2826
2827   //new TCanvas;
2828   //mc->Draw("COLZ");
2829
2830   TH1* proj = mc->ProjectionY();
2831
2832   TString tmp(fMultiplicityESD[id]->GetName());
2833   delete fMultiplicityESD[id];
2834   fMultiplicityESD[id] = CalculateMultiplicityESD(proj, id);
2835   fMultiplicityESD[id]->SetName(tmp);
2836 }