PWG0/dNdEta/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   fCurrentESD = fMultiplicityESD[inputRange]->ProjectionY("fCurrentESD"); //"fCurrentESD", -1, -1, "e");
 811   fCurrentESD->Sumw2();
 812
 813   // empty under/overflow bins in x, otherwise Project3D takes them into account
 814   TH3* hist = fCorrelation[correlationID];
 815   for (Int_t y=1; y<=hist->GetYaxis()->GetNbins(); ++y)
 816   {
 817     for (Int_t z=1; z<=hist->GetZaxis()->GetNbins(); ++z)
 818     {
 819       hist->SetBinContent(0, y, z, 0);
 820       hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0);
 821     }
 822   }
 823
 824   fCurrentCorrelation = hist->Project3D("zy");
 825   //((TH2*) fCurrentCorrelation)->Rebin2D(2, 1);
 826   //fMultiplicityESDCorrected[correlationID]->Rebin(2);
 827   fCurrentCorrelation->Sumw2();
 828
 829   Printf("AliMultiplicityCorrection::SetupCurrentHists: Statistics information: %.f entries in correlation map; %.f entries in measured spectrum", fCurrentCorrelation->Integral(), fCurrentESD->Integral());
 830
 831   if (createBigBin)
 832   {
 833     fLastBinLimit = 0;
 834     for (Int_t i=1; i<=fCurrentESD->GetNbinsX(); ++i)
 835     {
 836       if (fCurrentESD->GetBinContent(i) <= 5)
 837       {
 838         fLastBinLimit = i;
 839         break;
 840       }
 841     }
 842
 843     if (fLastBinLimit > 0)
 844     {
 845       TCanvas* canvas = 0;
 846
 847       if (display)
 848       {
 849         canvas = new TCanvas("StatSolution", "StatSolution", 1000, 800);
 850         canvas->Divide(2, 2);
 851
 852         canvas->cd(1);
 853         fCurrentESD->GetXaxis()->SetRangeUser(0, 200);
 854         fCurrentESD->SetStats(kFALSE);
 855         fCurrentESD->SetTitle(";measured multiplicity");
 856         fCurrentESD->DrawCopy();
 857         gPad->SetLogy();
 858
 859         canvas->cd(2);
 860         fCurrentCorrelation->GetXaxis()->SetRangeUser(0, 250);
 861         fCurrentCorrelation->GetYaxis()->SetRangeUser(0, 250);
 862         fCurrentCorrelation->SetTitle(";true multiplicity;measured multiplicity");
 863         fCurrentCorrelation->SetStats(kFALSE);
 864         fCurrentCorrelation->DrawCopy("COLZ");
 865       }
 866
 867       printf("Bin limit in measured spectrum is %d.\n", fLastBinLimit);
 868       fCurrentESD->SetBinContent(fLastBinLimit, fCurrentESD->Integral(fLastBinLimit, fCurrentESD->GetNbinsX()));
 869       for (Int_t i=fLastBinLimit+1; i<=fCurrentESD->GetNbinsX(); ++i)
 870       {
 871         fCurrentESD->SetBinContent(i, 0);
 872         fCurrentESD->SetBinError(i, 0);
 873       }
 874       // the error is set to sqrt(N), better solution possible?, sum of relative errors of all contributions???
 875       fCurrentESD->SetBinError(fLastBinLimit, TMath::Sqrt(fCurrentESD->GetBinContent(fLastBinLimit)));
 876
 877       printf("This bin has now %f +- %f entries\n", fCurrentESD->GetBinContent(fLastBinLimit), fCurrentESD->GetBinError(fLastBinLimit));
 878
 879       for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
 880       {
 881         fCurrentCorrelation->SetBinContent(i, fLastBinLimit, fCurrentCorrelation->Integral(i, i, fLastBinLimit, fCurrentCorrelation->GetNbinsY()));
 882         // the error is set to sqrt(N), better solution possible?, sum of relative errors of all contributions???
 883         fCurrentCorrelation->SetBinError(i, fLastBinLimit, TMath::Sqrt(fCurrentCorrelation->GetBinContent(i, fLastBinLimit)));
 884
 885         for (Int_t j=fLastBinLimit+1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
 886         {
 887           fCurrentCorrelation->SetBinContent(i, j, 0);
 888           fCurrentCorrelation->SetBinError(i, j, 0);
 889         }
 890       }
 891
 892       printf("Adjusted correlation matrix!\n");
 893
 894       if (canvas)
 895       {
 896         canvas->cd(3);
 897         fCurrentESD->DrawCopy();
 898         gPad->SetLogy();
 899
 900         canvas->cd(4);
 901         fCurrentCorrelation->DrawCopy("COLZ");
 902       }
 903     }
 904   }
 905
 906 #if 0 // does not help
 907   // null bins with one entry
 908   Int_t nNulledBins = 0;
 909   for (Int_t x=1; x<=fCurrentCorrelation->GetXaxis()->GetNbins(); ++x)
 910     for (Int_t y=1; y<=fCurrentCorrelation->GetYaxis()->GetNbins(); ++y)
 911     {
 912       if (fCurrentCorrelation->GetBinContent(x, y) == 1)
 913       {
 914         fCurrentCorrelation->SetBinContent(x, y, 0);
 915         fCurrentCorrelation->SetBinError(x, y, 0);
 916
 917         ++nNulledBins;
 918       }
 919     }
 920   Printf("Nulled %d bins", nNulledBins);
 921 #endif
 922
 923   fCurrentEfficiency = GetEfficiency(inputRange, eventType);
 924   //fCurrentEfficiency->Rebin(2);
 925   //fCurrentEfficiency->Scale(0.5);
 926 }
 927
 928 //____________________________________________________________________
 929 TH1* AliMultiplicityCorrection::GetEfficiency(Int_t inputRange, EventType eventType)
 930 {
 931   //
 932   // calculates efficiency for given event type
 933   //
 934
 935   TH1* divisor = 0;
 936   switch (eventType)
 937   {
 938     case kTrVtx : divisor = fMultiplicityVtx[inputRange]->ProjectionY("divisor", -1, -1, "e"); break;
 939     case kMB: divisor = fMultiplicityMB[inputRange]->ProjectionY("divisor", -1, -1, "e"); break;
 940     case kINEL: divisor = fMultiplicityINEL[inputRange]->ProjectionY("divisor", -1, -1, "e"); break;
 941   }
 942   TH1* eff = fMultiplicityVtx[inputRange]->ProjectionY("CurrentEfficiency", -1, -1, "e");
 943   eff->Divide(eff, divisor, 1, 1, "B");
 944   return eff;
 945 }
 946
 947 //____________________________________________________________________
 948 void AliMultiplicityCorrection::SetRegularizationParameters(RegularizationType type, Float_t weight, Int_t minuitParams)
 949 {
 950   //
 951   // sets the parameters for chi2 minimization
 952   //
 953
 954   fgRegularizationType = type;
 955   fgRegularizationWeight = weight;
 956
 957   printf("AliMultiplicityCorrection::SetRegularizationParameters --> Regularization set to %d with weight %f\n", (Int_t) type, weight);
 958
 959   if (minuitParams > 0)
 960   {
 961     fgNParamsMinuit = minuitParams;
 962     printf("AliMultiplicityCorrection::SetRegularizationParameters --> Number of MINUIT iterations set to %d\n", minuitParams);
 963   }
 964 }
 965
 966 //____________________________________________________________________
 967 void AliMultiplicityCorrection::SetBayesianParameters(Float_t smoothing, Int_t nIterations)
 968 {
 969   //
 970   // sets the parameters for Bayesian unfolding
 971   //
 972
 973   fgBayesianSmoothing = smoothing;
 974   fgBayesianIterations = nIterations;
 975
 976   printf("AliMultiplicityCorrection::SetBayesianParameters --> Paramaters set to %d iterations with smoothing %f\n", fgBayesianIterations, fgBayesianSmoothing);
 977 }
 978
 979 //____________________________________________________________________
 980 Int_t AliMultiplicityCorrection::UnfoldWithMinuit(TH1* correlation, TH1* aEfficiency, TH1* measured, TH1* initialConditions, TH1* result, Bool_t check)
 981 {
 982   //
 983   // implemenation of unfolding (internal function)
 984   //
 985   // unfolds <measured> using response from <correlation> and effiency <aEfficiency>
 986   // output is in <result>
 987   // <initialConditions> set the initial values for the minimization, if 0 <measured> is used
 988   // if <check> is true no unfolding is made, instead only the chi2 without unfolding is printed
 989   //
 990   // returns minuit status (0 = success), (-1 when check was set)
 991   //
 992
 993   //normalize measured
 994   measured->Scale(1.0 / measured->Integral());
 995
 996   if (!fgCorrelationMatrix)
 997     fgCorrelationMatrix = new TMatrixD(fgkMaxInput, fgNParamsMinuit);
 998   if (!fgCorrelationCovarianceMatrix)
 999     fgCorrelationCovarianceMatrix = new TMatrixD(fgkMaxInput, fgkMaxInput);
1000   if (!fgCurrentESDVector)
1001     fgCurrentESDVector = new TVectorD(fgkMaxInput);
1002   if (!fgEntropyAPriori)
1003     fgEntropyAPriori = new TVectorD(fgNParamsMinuit);
1004
1005   fgCorrelationMatrix->Zero();
1006   fgCorrelationCovarianceMatrix->Zero();
1007   fgCurrentESDVector->Zero();
1008   fgEntropyAPriori->Zero();
1009
1010   // normalize correction for given nPart
1011   for (Int_t i=1; i<=correlation->GetNbinsX(); ++i)
1012   {
1013     Double_t sum = correlation->Integral(i, i, 1, correlation->GetNbinsY());
1014     if (sum <= 0)
1015       continue;
1016     Float_t maxValue = 0;
1017     Int_t maxBin = -1;
1018     for (Int_t j=1; j<=correlation->GetNbinsY(); ++j)
1019     {
1020       // find most probably value
1021       if (maxValue < correlation->GetBinContent(i, j))
1022       {
1023         maxValue = correlation->GetBinContent(i, j);
1024         maxBin = j;
1025       }
1026
1027       // npart sum to 1
1028       correlation->SetBinContent(i, j, correlation->GetBinContent(i, j) / sum);
1029       correlation->SetBinError(i, j, correlation->GetBinError(i, j) / sum);
1030
1031       if (i <= fgNParamsMinuit && j <= fgkMaxInput)
1032         (*fgCorrelationMatrix)(j-1, i-1) = correlation->GetBinContent(i, j);
1033     }
1034
1035     //printf("MPV for Ntrue = %f is %f\n", fCurrentCorrelation->GetXaxis()->GetBinCenter(i), fCurrentCorrelation->GetYaxis()->GetBinCenter(maxBin));
1036   }
1037
1038   // Initialize TMinuit via generic fitter interface
1039   TVirtualFitter *minuit = TVirtualFitter::Fitter(0, fgNParamsMinuit);
1040   Double_t arglist[100];
1041
1042   // disable any output (-1), unfortuantly we do not see warnings anymore then. Have to find another way...
1043   arglist[0] = 0;
1044   minuit->ExecuteCommand("SET PRINT", arglist, 1);
1045
1046   // however, enable warnings
1047   //minuit->ExecuteCommand("SET WAR", arglist, 0);
1048
1049   // set minimization function
1050   minuit->SetFCN(MinuitFitFunction);
1051
1052   for (Int_t i=0; i<fgNParamsMinuit; i++)
1053     (*fgEntropyAPriori)[i] = 1;
1054
1055   if (!initialConditions) {
1056     initialConditions = measured;
1057   } else {
1058     printf("Using different starting conditions...\n");
1059     //new TCanvas; initialConditions->DrawCopy();
1060     initialConditions->Scale(1.0 / initialConditions->Integral());
1061   }
1062
1063   for (Int_t i=0; i<fgNParamsMinuit; i++)
1064     if (initialConditions->GetBinContent(i+1) > 0)
1065       (*fgEntropyAPriori)[i] = initialConditions->GetBinContent(i+1);
1066
1067   Double_t* results = new Double_t[fgNParamsMinuit+1];
1068   for (Int_t i=0; i<fgNParamsMinuit; ++i)
1069   {
1070     results[i] = initialConditions->GetBinContent(i+1);
1071
1072     // minimum value
1073     if (!check)
1074       results[i] = TMath::Max(results[i], 1e-3);
1075
1076     Float_t stepSize = 0.1;
1077
1078     // minuit sees squared values to prevent it from going negative...
1079     results[i] = TMath::Sqrt(results[i]);
1080     //stepSize /= results[i]; // keep relative step size
1081
1082     minuit->SetParameter(i, Form("param%d", i), results[i], stepSize, 0, 0);
1083   }
1084   //minuit->SetParameter(fgkMaxParams, "alpha", 1e4, 1, 0, 0);
1085   // bin 0 is filled with value from bin 1 (otherwise it's 0)
1086   //minuit->SetParameter(0, "param0", (results[1] > minStartValue) ? results[1] : minStartValue, 0.1, 0, 0);
1087   //results[0] = 0;
1088   //minuit->SetParameter(0, "param0", 0, 0, 0, 0);
1089
1090   for (Int_t i=0; i<fgkMaxInput; ++i)
1091   {
1092     (*fgCurrentESDVector)[i] = measured->GetBinContent(i+1);
1093     if (measured->GetBinError(i+1) > 0)
1094       (*fgCorrelationCovarianceMatrix)(i, i) = (Double_t) 1e-6 / measured->GetBinError(i+1) / measured->GetBinError(i+1);
1095
1096     if ((*fgCorrelationCovarianceMatrix)(i, i) > 1e7)
1097       (*fgCorrelationCovarianceMatrix)(i, i) = 0;
1098   }
1099
1100   Int_t dummy = 0;
1101   Double_t chi2 = 0;
1102   MinuitFitFunction(dummy, 0, chi2, results, 0);
1103   DrawGuess(results);
1104   printf("Chi2 of initial parameters is = %f\n", chi2);
1105
1106   if (check)
1107   {
1108     delete[] results;
1109     return -1;
1110   }
1111
1112   // first param is number of iterations, second is precision....
1113   arglist[0] = 1e6;
1114   //arglist[1] = 1e-5;
1115   //minuit->ExecuteCommand("SCAN", arglist, 0);
1116   Int_t status = minuit->ExecuteCommand("MIGRAD", arglist, 1);
1117   printf("MINUIT status is %d\n", status);
1118   //minuit->ExecuteCommand("MIGRAD", arglist, 1);
1119   //minuit->ExecuteCommand("MIGRAD", arglist, 1);
1120   //printf("!!!!!!!!!!!!!! MIGRAD finished: Starting MINOS !!!!!!!!!!!!!!");
1121   //minuit->ExecuteCommand("MINOS", arglist, 0);
1122
1123   for (Int_t i=0; i<fgNParamsMinuit; ++i)
1124   {
1125     if (aEfficiency->GetBinContent(i+1) > 0)
1126     {
1127       results[i] = minuit->GetParameter(i);
1128       result->SetBinContent(i+1, results[i] * results[i] / aEfficiency->GetBinContent(i+1));
1129       // error is : (relError) * (value) = (minuit->GetParError(i) / minuit->GetParameter(i)) * (minuit->GetParameter(i) * minuit->GetParameter(i))
1130       result->SetBinError(i+1, minuit->GetParError(i) * results[i] /  aEfficiency->GetBinContent(i+1));
1131     }
1132     else
1133     {
1134       result->SetBinContent(i+1, 0);
1135       result->SetBinError(i+1, 0);
1136     }
1137   }
1138   DrawGuess(results);
1139
1140   delete[] results;
1141
1142   return status;
1143 }
1144
1145 //____________________________________________________________________
1146 Int_t AliMultiplicityCorrection::ApplyMinuitFit(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t check, TH1* initialConditions)
1147 {
1148   //
1149   // correct spectrum using minuit chi2 method
1150   //
1151   // for description of parameters, see UnfoldWithMinuit
1152   //
1153
1154   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
1155   SetupCurrentHists(inputRange, fullPhaseSpace, eventType, fgCreateBigBin);
1156
1157   return UnfoldWithMinuit(fCurrentCorrelation, fCurrentEfficiency, fCurrentESD, initialConditions, fMultiplicityESDCorrected[correlationID], check);
1158 }
1159
1160 //____________________________________________________________________
1161 void AliMultiplicityCorrection::ApplyNBDFit(Int_t inputRange, Bool_t fullPhaseSpace)
1162 {
1163   //
1164   // correct spectrum using minuit chi2 method applying a NBD fit
1165   //
1166
1167   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
1168
1169   SetupCurrentHists(inputRange, fullPhaseSpace, kTrVtx, kFALSE);
1170   //normalize ESD
1171   fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
1172
1173   fgCorrelationMatrix = new TMatrixD(fgkMaxParams, fgkMaxParams);
1174   fgCorrelationCovarianceMatrix = new TMatrixD(fgkMaxParams, fgkMaxParams);
1175   fgCurrentESDVector = new TVectorD(fgkMaxParams);
1176
1177   // normalize correction for given nPart
1178   for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
1179   {
1180     Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
1181     if (sum <= 0)
1182       continue;
1183     for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
1184     {
1185       // npart sum to 1
1186       fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
1187       fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
1188
1189       if (i <= fgkMaxParams && j <= fgkMaxParams)
1190         (*fgCorrelationMatrix)(j-1, i-1) = fCurrentCorrelation->GetBinContent(i, j);
1191     }
1192   }
1193
1194   for (Int_t i=0; i<fgkMaxParams; ++i)
1195   {
1196     (*fgCurrentESDVector)[i] = fCurrentESD->GetBinContent(i+1);
1197     if (fCurrentESD->GetBinError(i+1) > 0)
1198       (*fgCorrelationCovarianceMatrix)(i, i) = 1.0 / (fCurrentESD->GetBinError(i+1) * fCurrentESD->GetBinError(i+1));
1199   }
1200
1201   // Create NBD function
1202   if (!fgNBD)
1203   {
1204     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);
1205     fgNBD->SetParNames("scaling", "averagen", "k");
1206   }
1207
1208   // Initialize TMinuit via generic fitter interface
1209   TVirtualFitter *minuit = TVirtualFitter::Fitter(0, 3);
1210
1211   minuit->SetFCN(MinuitNBD);
1212   minuit->SetParameter(0, "param0", 1, 0.1, 0.001, fCurrentESD->GetMaximum() * 1000);
1213   minuit->SetParameter(1, "param1", 10, 1, 0.001, 1000);
1214   minuit->SetParameter(2, "param2", 2, 1, 0.001, 1000);
1215
1216   Double_t arglist[100];
1217   arglist[0] = 0;
1218   minuit->ExecuteCommand("SET PRINT", arglist, 1);
1219   minuit->ExecuteCommand("MIGRAD", arglist, 0);
1220   //minuit->ExecuteCommand("MINOS", arglist, 0);
1221
1222   fgNBD->SetParameters(minuit->GetParameter(0), minuit->GetParameter(1), minuit->GetParameter(2));
1223
1224   for (Int_t i=0; i<fgkMaxParams; ++i)
1225   {
1226     printf("%d : %f\n", i, fgNBD->Eval(i));
1227     if (fgNBD->Eval(i) > 0)
1228     {
1229       fMultiplicityESDCorrected[correlationID]->SetBinContent(i+1, fgNBD->Eval(i));
1230       fMultiplicityESDCorrected[correlationID]->SetBinError(i+1, 0);
1231     }
1232   }
1233 }
1234
1235 //____________________________________________________________________
1236 void AliMultiplicityCorrection::DrawHistograms()
1237 {
1238   //
1239   // draws the histograms of this class
1240   //
1241
1242   printf("ESD:\n");
1243
1244   TCanvas* canvas1 = new TCanvas("fMultiplicityESD", "fMultiplicityESD", 900, 600);
1245   canvas1->Divide(3, 2);
1246   for (Int_t i = 0; i < kESDHists; ++i)
1247   {
1248     canvas1->cd(i+1);
1249     fMultiplicityESD[i]->DrawCopy("COLZ");
1250     printf("%d --> %f\n", i, (Float_t) fMultiplicityESD[i]->ProjectionY()->GetMean());
1251   }
1252
1253   printf("Vtx:\n");
1254
1255   TCanvas* canvas2 = new TCanvas("fMultiplicityMC", "fMultiplicityMC", 900, 600);
1256   canvas2->Divide(3, 2);
1257   for (Int_t i = 0; i < kMCHists; ++i)
1258   {
1259     canvas2->cd(i+1);
1260     fMultiplicityVtx[i]->DrawCopy("COLZ");
1261     printf("%d --> %f\n", i, fMultiplicityVtx[i]->ProjectionY()->GetMean());
1262   }
1263
1264   TCanvas* canvas3 = new TCanvas("fCorrelation", "fCorrelation", 900, 900);
1265   canvas3->Divide(3, 3);
1266   for (Int_t i = 0; i < kCorrHists; ++i)
1267   {
1268     canvas3->cd(i+1);
1269     TH3* hist = dynamic_cast<TH3*> (fCorrelation[i]->Clone());
1270     for (Int_t y=1; y<=hist->GetYaxis()->GetNbins(); ++y)
1271     {
1272       for (Int_t z=1; z<=hist->GetZaxis()->GetNbins(); ++z)
1273       {
1274         hist->SetBinContent(0, y, z, 0);
1275         hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0);
1276       }
1277     }
1278     TH1* proj = hist->Project3D("zy");
1279     proj->DrawCopy("COLZ");
1280   }
1281 }
1282
1283 //____________________________________________________________________
1284 void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRange, Bool_t fullPhaseSpace, Bool_t /*normalizeESD*/, TH1* mcHist, Bool_t simple)
1285 {
1286   //mcHist->Rebin(2);
1287
1288   Int_t esdCorrId = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
1289
1290   TString tmpStr;
1291   tmpStr.Form("%s_DrawComparison_%d", name, esdCorrId);
1292
1293   if (fMultiplicityESDCorrected[esdCorrId]->Integral() == 0)
1294   {
1295     printf("ERROR. Unfolded histogram is empty\n");
1296     return;
1297   }
1298
1299   //regain measured distribution used for unfolding, because the bins at high mult. were modified in SetupCurrentHists
1300   fCurrentESD = fMultiplicityESD[esdCorrId]->ProjectionY();
1301   fCurrentESD->Sumw2();
1302   fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
1303
1304   // normalize unfolded result to 1
1305   fMultiplicityESDCorrected[esdCorrId]->Scale(1.0 / fMultiplicityESDCorrected[esdCorrId]->Integral());
1306
1307   //fCurrentESD->Scale(mcHist->Integral(2, 200));
1308
1309   //new TCanvas;
1310   /*TH1* ratio = (TH1*) fMultiplicityESDCorrected[esdCorrId]->Clone("ratio");
1311   ratio->Divide(mcHist);
1312   ratio->Draw("HIST");
1313   ratio->Fit("pol0", "W0", "", 20, 120);
1314   Float_t scalingFactor = ratio->GetFunction("pol0")->GetParameter(0);
1315   delete ratio;
1316   mcHist->Scale(scalingFactor);*/
1317
1318   // find bin with <= 50 entries. this is used as limit for the chi2 calculation
1319   Int_t mcBinLimit = 0;
1320   for (Int_t i=20; i<mcHist->GetNbinsX(); ++i)
1321   {
1322     if (mcHist->GetBinContent(i) > 50)
1323     {
1324       mcBinLimit = i;
1325     }
1326     else
1327       break;
1328   }
1329   Printf("AliMultiplicityCorrection::DrawComparison: MC bin limit is %d", mcBinLimit);
1330
1331   // scale to 1
1332   mcHist->Sumw2();
1333   if (mcHist->Integral() > 0)
1334     mcHist->Scale(1.0 / mcHist->Integral());
1335
1336   // calculate residual
1337
1338   // for that we convolute the response matrix with the gathered result
1339   // if normalizeESD is not set, the histogram is already normalized, this needs to be passed to CalculateMultiplicityESD
1340   TH1* tmpESDEfficiencyRecorrected = (TH1*) fMultiplicityESDCorrected[esdCorrId]->Clone("tmpESDEfficiencyRecorrected");
1341   tmpESDEfficiencyRecorrected->Multiply(fCurrentEfficiency);
1342   TH2* convoluted = CalculateMultiplicityESD(tmpESDEfficiencyRecorrected, esdCorrId);
1343   TH1* convolutedProj = convoluted->ProjectionY("convolutedProj", -1, -1, "e");
1344   if (convolutedProj->Integral() > 0)
1345   {
1346     convolutedProj->Scale(1.0 / convolutedProj->Integral());
1347   }
1348   else
1349     printf("ERROR: convolutedProj is empty. Something went wrong calculating the convoluted histogram.\n");
1350
1351   TH1* residual = (TH1*) convolutedProj->Clone("residual");
1352   residual->SetTitle("Residuals;Ntracks;(folded unfolded measured - measured) / e");
1353
1354   residual->Add(fCurrentESD, -1);
1355   //residual->Divide(residual, fCurrentESD, 1, 1, "B");
1356
1357   TH1* residualHist = new TH1F("residualHist", "residualHist", 51, -5, 5);
1358
1359   // TODO fix errors
1360   Float_t chi2 = 0;
1361   for (Int_t i=1; i<=residual->GetNbinsX(); ++i)
1362   {
1363     if (fCurrentESD->GetBinError(i) > 0)
1364     {
1365       Float_t value = residual->GetBinContent(i) / fCurrentESD->GetBinError(i);
1366       if (i > 1 && i <= fLastBinLimit)
1367         chi2 += value * value;
1368       Printf("%d --> %f (%f)", i, value * value, chi2);
1369       residual->SetBinContent(i, value);
1370       residualHist->Fill(value);
1371     }
1372     else
1373     {
1374       //printf("Residual bin %d set to 0\n", i);
1375       residual->SetBinContent(i, 0);
1376     }
1377     convolutedProj->SetBinError(i, 0);
1378     residual->SetBinError(i, 0);
1379   }
1380   fLastChi2Residuals = chi2;
1381
1382   new TCanvas; residualHist->DrawCopy();
1383
1384   //residualHist->Fit("gaus", "N");
1385   //delete residualHist;
1386
1387   printf("Difference (Residuals) is %f for bin 2-%d\n", fLastChi2Residuals, fLastBinLimit);
1388
1389   TCanvas* canvas1 = 0;
1390   if (simple)
1391   {
1392     canvas1 = new TCanvas(tmpStr, tmpStr, 900, 400);
1393     canvas1->Divide(2, 1);
1394   }
1395   else
1396   {
1397     canvas1 = new TCanvas(tmpStr, tmpStr, 1200, 1200);
1398     canvas1->Divide(2, 3);
1399   }
1400
1401   canvas1->cd(1);
1402   TH1* proj = (TH1*) mcHist->Clone("proj");
1403
1404   proj->GetXaxis()->SetRangeUser(0, 200);
1405   proj->SetTitle(";true multiplicity;Entries");
1406   proj->SetStats(kFALSE);
1407
1408   fMultiplicityESDCorrected[esdCorrId]->SetLineColor(2);
1409
1410   if (proj->GetEntries() > 0) {
1411     proj->DrawCopy("HIST");
1412     fMultiplicityESDCorrected[esdCorrId]->DrawCopy("SAME HIST E");
1413   }
1414   else
1415     fMultiplicityESDCorrected[esdCorrId]->DrawCopy("HIST E");
1416
1417   gPad->SetLogy();
1418
1419   TLegend* legend = new TLegend(0.3, 0.8, 0.93, 0.93);
1420   legend->AddEntry(proj, "true distribution");
1421   legend->AddEntry(fMultiplicityESDCorrected[esdCorrId], "unfolded distribution");
1422   legend->SetFillColor(0);
1423   legend->Draw();
1424   // unfortunately does not work. maybe a bug? --> legend->SetTextSizePixels(14);
1425
1426   canvas1->cd(2);
1427
1428   gPad->SetLogy();
1429   fCurrentESD->GetXaxis()->SetRangeUser(0, 200);
1430   //fCurrentESD->SetLineColor(2);
1431   fCurrentESD->SetTitle(";measured multiplicity;Entries");
1432   fCurrentESD->SetStats(kFALSE);
1433   fCurrentESD->DrawCopy("HIST E");
1434
1435   convolutedProj->SetLineColor(2);
1436   //proj2->SetMarkerColor(2);
1437   //proj2->SetMarkerStyle(5);
1438   convolutedProj->DrawCopy("HIST SAME");
1439
1440   legend = new TLegend(0.3, 0.8, 0.93, 0.93);
1441   legend->AddEntry(fCurrentESD, "measured distribution");
1442   legend->AddEntry(convolutedProj, "R #otimes unfolded distribution");
1443   legend->SetFillColor(0);
1444   legend->Draw();
1445
1446   //TH1* diffMCUnfolded = dynamic_cast<TH1*> (proj->Clone("diffMCUnfolded"));
1447   //diffMCUnfolded->Add(fMultiplicityESDCorrected[esdCorrId], -1);
1448
1449   /*Float_t chi2 = 0;
1450   Float_t chi = 0;
1451   fLastChi2MCLimit = 0;
1452   Int_t limit = 0;
1453   for (Int_t i=2; i<=200; ++i)
1454   {
1455     if (proj->GetBinContent(i) != 0)
1456     {
1457       Float_t value = (proj->GetBinContent(i) - fMultiplicityESDCorrected[esdCorrId]->GetBinContent(i)) / proj->GetBinContent(i);
1458       chi2 += value * value;
1459       chi += TMath::Abs(value);
1460
1461       //printf("%d %f\n", i, chi);
1462
1463       if (chi2 < 0.2)
1464         fLastChi2MCLimit = i;
1465
1466       if (chi < 3)
1467         limit = i;
1468
1469     }
1470   }*/
1471
1472   /*chi2 = 0;
1473   Float_t chi = 0;
1474   Int_t limit = 0;
1475   for (Int_t i=1; i<=diffMCUnfolded->GetNbinsX(); ++i)
1476   {
1477     if (fMultiplicityESDCorrected[esdCorrId]->GetBinError(i) > 0)
1478     {
1479       Double_t value = diffMCUnfolded->GetBinContent(i) / fMultiplicityESDCorrected[esdCorrId]->GetBinError(i);
1480       if (value > 1e8)
1481         value = 1e8; //prevent arithmetic exception
1482       else if (value < -1e8)
1483         value = -1e8;
1484       if (i > 1)
1485       {
1486         chi2 += value * value;
1487         chi += TMath::Abs(value);
1488       }
1489       diffMCUnfolded->SetBinContent(i, value);
1490     }
1491     else
1492     {
1493       //printf("diffMCUnfolded bin %d set to 0\n", i);
1494       diffMCUnfolded->SetBinContent(i, 0);
1495     }
1496     if (chi2 < 1000)
1497       fLastChi2MCLimit = i;
1498     if (chi < 1000)
1499       limit = i;
1500     if (i == 150)
1501       fLastChi2MC = chi2;
1502   }
1503
1504   printf("limits %d %d\n", fLastChi2MCLimit, limit);
1505   fLastChi2MCLimit = limit;
1506
1507   printf("Difference (from MC) is %f for bin 2-150. Limit is %d.\n", fLastChi2MC, fLastChi2MCLimit);*/
1508
1509   if (!simple)
1510   {
1511     /*canvas1->cd(3);
1512
1513     diffMCUnfolded->SetTitle("#chi^{2};Npart;(MC - Unfolded) / e(unfolded)");
1514     //diffMCUnfolded->GetYaxis()->SetRangeUser(-20, 20);
1515     diffMCUnfolded->GetXaxis()->SetRangeUser(0, 200);
1516     diffMCUnfolded->DrawCopy("HIST");
1517
1518     TH1F* fluctuation = new TH1F("fluctuation", "fluctuation", 20, -5, 5);
1519     for (Int_t i=20; i<=diffMCUnfolded->GetNbinsX(); ++i)
1520       fluctuation->Fill(diffMCUnfolded->GetBinContent(i));
1521
1522     //new TCanvas; fluctuation->DrawCopy();
1523     delete fluctuation;*/
1524
1525     /*TLegend* legend = new TLegend(0.6, 0.7, 0.85, 0.85);
1526     legend->AddEntry(fMultiplicityESDCorrected, "ESD corrected");
1527     legend->AddEntry(fMultiplicityMC, "MC");
1528     legend->AddEntry(fMultiplicityESD, "ESD");
1529     legend->Draw();*/
1530
1531     canvas1->cd(4);
1532     //residual->GetYaxis()->SetRangeUser(-0.2, 0.2);
1533     residual->GetXaxis()->SetRangeUser(0, 200);
1534     residual->DrawCopy();
1535
1536     canvas1->cd(5);
1537
1538     TH1* ratio = (TH1*) fMultiplicityESDCorrected[esdCorrId]->Clone("ratio");
1539     ratio->Divide(mcHist);
1540     ratio->SetTitle("Ratio;true multiplicity;Unfolded / MC");
1541     ratio->GetYaxis()->SetRangeUser(0.5, 1.5);
1542     ratio->GetXaxis()->SetRangeUser(0, 200);
1543     ratio->SetStats(kFALSE);
1544     ratio->Draw("HIST");
1545
1546     Double_t ratioChi2 = 0;
1547     fRatioAverage = 0;
1548     fLastChi2MCLimit = 0;
1549     for (Int_t i=2; i<=150; ++i)
1550     {
1551       Float_t value = ratio->GetBinContent(i) - 1;
1552       if (value > 1e8)
1553         value = 1e8; //prevent arithmetic exception
1554       else if (value < -1e8)
1555         value = -1e8;
1556
1557       ratioChi2 += value * value;
1558       fRatioAverage += TMath::Abs(value);
1559
1560       if (ratioChi2 < 0.1)
1561         fLastChi2MCLimit = i;
1562     }
1563     fRatioAverage /= 149;
1564
1565     printf("Sum over (ratio-1)^2 (2..150) is %f; average of |ratio-1| is %f\n", ratioChi2, fRatioAverage);
1566     // TODO FAKE
1567     fLastChi2MC = ratioChi2;
1568
1569     // FFT of ratio
1570     canvas1->cd(6);
1571     const Int_t kFFT = 128;
1572     Double_t fftReal[kFFT];
1573     Double_t fftImag[kFFT];
1574
1575     for (Int_t i=0; i<kFFT; ++i)
1576     {
1577       fftReal[i] = ratio->GetBinContent(i+1+10);
1578       // test: ;-)
1579       //fftReal[i] = cos(TMath::Pi() * 5 * 2 * i / 128);
1580       fftImag[i] = 0;
1581     }
1582
1583     FFT(-1, TMath::Nint(log(kFFT) / log(2)), fftReal, fftImag);
1584
1585     TH1* fftResult = (TH1*) ratio->Clone("fftResult");
1586     fftResult->SetTitle("FFT;true multiplicity;coeff. (10...137)");
1587     TH1* fftResult2 = (TH1*) ratio->Clone("fftResult2");
1588     TH1* fftResult3 = (TH1*) ratio->Clone("fftResult3");
1589     fftResult->Reset();
1590     fftResult2->Reset();
1591     fftResult3->Reset();
1592     fftResult->GetYaxis()->UnZoom();
1593     fftResult2->GetYaxis()->UnZoom();
1594     fftResult3->GetYaxis()->UnZoom();
1595
1596     //Printf("AFTER FFT");
1597     for (Int_t i=0; i<kFFT; ++i)
1598     {
1599       //Printf("%d: %f", i, fftReal[i]);
1600       fftResult->SetBinContent(i+1, fftReal[i]);
1601       /*if (i != 0 && TMath::Abs(fftReal[i]) > 0.5)
1602       {
1603         Printf("Nulled %d", i);
1604         fftReal[i] = 0;
1605       }*/
1606     }
1607
1608     fftResult->SetLineColor(1);
1609     fftResult->DrawCopy();
1610
1611
1612     //Printf("IMAG");
1613     for (Int_t i=0; i<kFFT; ++i)
1614     {
1615       //Printf("%d: %f", i, fftImag[i]);
1616       fftResult2->SetBinContent(i+1, fftImag[i]);
1617
1618       fftResult3->SetBinContent(i+1, TMath::Sqrt(fftReal[i] * fftReal[i] + fftImag[i] * fftImag[i]));
1619     }
1620
1621     fftResult2->SetLineColor(2);
1622     fftResult2->DrawCopy("SAME");
1623
1624     fftResult3->SetLineColor(4);
1625     fftResult3->DrawCopy("SAME");
1626
1627     for (Int_t i=1; i<kFFT - 1; ++i)
1628     {
1629       if (TMath::Sqrt(fftReal[i] * fftReal[i] + fftImag[i] * fftImag[i]) > 3)
1630       {
1631         fftReal[i-1] = 0;
1632         fftReal[i] = 0;
1633         fftReal[i+1] = 0;
1634         fftImag[i-1] = 0;
1635         fftImag[i] = 0;
1636         fftImag[i+1] = 0;
1637         //fftReal[i] = (fftReal[i-1] + fftReal[i+1]) / 2;
1638         //fftImag[i]  = (fftImag[i-1] + fftImag[i+1]) / 2;
1639         //Printf("Nulled %d to %f %f", i, fftReal[i], fftImag[i]);
1640       }
1641     }
1642
1643     FFT(1, TMath::Nint(log(kFFT) / log(2)), fftReal, fftImag);
1644
1645     TH1* fftResult4 = (TH1*) fftResult3->Clone("fftResult4");
1646     fftResult4->Reset();
1647
1648     //Printf("AFTER BACK-TRAFO");
1649     for (Int_t i=0; i<kFFT; ++i)
1650     {
1651       //Printf("%d: %f", i, fftReal[i]);
1652       fftResult4->SetBinContent(i+1+10, fftReal[i]);
1653     }
1654
1655     canvas1->cd(5);
1656     fftResult4->SetLineColor(4);
1657     fftResult4->DrawCopy("SAME");
1658
1659     // plot (MC - Unfolded) / error (MC)
1660     canvas1->cd(3);
1661
1662     TH1* diffMCUnfolded2 = dynamic_cast<TH1*> (proj->Clone("diffMCUnfolded2"));
1663     diffMCUnfolded2->Add(fMultiplicityESDCorrected[esdCorrId], -1);
1664
1665     Int_t ndfQual[kQualityRegions];
1666     for (Int_t region=0; region<kQualityRegions; ++region)
1667     {
1668       fQuality[region] = 0;
1669       ndfQual[region] = 0;
1670     }
1671
1672
1673     Double_t newChi2 = 0;
1674     Double_t newChi2Limit150 = 0;
1675     Int_t ndf = 0;
1676     for (Int_t i=1; i<=TMath::Min(diffMCUnfolded2->GetNbinsX(), fgkMaxParams+1); ++i)
1677     {
1678       Double_t value = 0;
1679       if (proj->GetBinError(i) > 0)
1680       {
1681         value = diffMCUnfolded2->GetBinContent(i) / proj->GetBinError(i);
1682         newChi2 += value * value;
1683         if (i > 1 && i <= mcBinLimit)
1684           newChi2Limit150 += value * value;
1685         ++ndf;
1686
1687         for (Int_t region=0; region<kQualityRegions; ++region)
1688           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
1689           {
1690             fQuality[region] += TMath::Abs(value);
1691             ++ndfQual[region];
1692           }
1693       }
1694
1695       diffMCUnfolded2->SetBinContent(i, value);
1696     }
1697
1698     // normalize region to the number of entries
1699     for (Int_t region=0; region<kQualityRegions; ++region)
1700     {
1701       if (ndfQual[region] > 0)
1702         fQuality[region] /= ndfQual[region];
1703       Printf("Quality parameter %d (%d <= mult <= %d) is %f with %d df", region, fgQualityRegionsB[region], fgQualityRegionsE[region], fQuality[region], ndfQual[region]);
1704     }
1705
1706     if (mcBinLimit > 1)
1707     {
1708       // TODO another FAKE
1709       fLastChi2MC = newChi2Limit150 / (mcBinLimit - 1);
1710       Printf("Chi2 (2..%d) from (MC - Unfolded) / e(MC) is: %.2f ndf is %d --> chi2 / ndf = %.2f", mcBinLimit, newChi2Limit150, mcBinLimit - 1, fLastChi2MC);
1711     }
1712     else
1713       fLastChi2MC = -1;
1714
1715     Printf("Chi2 (full range) from (MC - Unfolded) / e(MC) is: %.2f ndf is %d --> chi2 / ndf = %.2f", newChi2, ndf, ((ndf > 0) ? newChi2 / ndf : -1));
1716
1717     diffMCUnfolded2->SetTitle("#chi^{2};Npart;(MC - Unfolded) / e(MC)");
1718     //diffMCUnfolded->GetYaxis()->SetRangeUser(-20, 20);
1719     diffMCUnfolded2->GetXaxis()->SetRangeUser(0, 200);
1720     diffMCUnfolded2->DrawCopy("HIST");
1721   }
1722
1723   canvas1->SaveAs(Form("%s.gif", canvas1->GetName()));
1724 }
1725
1726 //____________________________________________________________________
1727 void AliMultiplicityCorrection::FFT(Int_t dir, Int_t m, Double_t *x, Double_t *y)
1728 {
1729 /*-------------------------------------------------------------------------
1730    This computes an in-place complex-to-complex FFT
1731    x and y are the real and imaginary arrays of 2^m points.
1732    dir =  1 gives forward transform
1733    dir = -1 gives reverse transform
1734
1735      Formula: forward
1736                   N-1
1737                   ---
1738               1   \          - j k 2 pi n / N
1739       X(n) = ---   >   x(k) e                    = forward transform
1740               N   /                                n=0..N-1
1741                   ---
1742                   k=0
1743
1744       Formula: reverse
1745                   N-1
1746                   ---
1747                   \          j k 2 pi n / N
1748       X(n) =       >   x(k) e                    = forward transform
1749                   /                                n=0..N-1
1750                   ---
1751                   k=0
1752 */
1753
1754    Long_t   nn, i, i1, j, k, i2, l, l1, l2;
1755    Double_t c1, c2, tx, ty, t1, t2, u1, u2, z;
1756
1757    /* Calculate the number of points */
1758    nn = 1;
1759    for (i = 0; i < m; i++)
1760        nn *= 2;
1761
1762    /* Do the bit reversal */
1763    i2 = nn >> 1;
1764    j = 0;
1765    for (i= 0; i < nn - 1; i++) {
1766        if (i < j) {
1767            tx = x[i];
1768            ty = y[i];
1769            x[i] = x[j];
1770            y[i] = y[j];
1771            x[j] = tx;
1772            y[j] = ty;
1773        }
1774        k = i2;
1775        while (k <= j) {
1776            j -= k;
1777            k >>= 1;
1778        }
1779        j += k;
1780    }
1781
1782    /* Compute the FFT */
1783    c1 = -1.0;
1784    c2 = 0.0;
1785    l2 = 1;
1786    for (l = 0; l < m; l++) {
1787        l1 = l2;
1788        l2 <<= 1;
1789        u1 = 1.0;
1790        u2 = 0.0;
1791        for (j = 0;j < l1; j++) {
1792            for (i = j; i < nn; i += l2) {
1793                i1 = i + l1;
1794                t1 = u1 * x[i1] - u2 * y[i1];
1795                t2 = u1 * y[i1] + u2 * x[i1];
1796                x[i1] = x[i] - t1;
1797                y[i1] = y[i] - t2;
1798                x[i] += t1;
1799                y[i] += t2;
1800            }
1801            z =  u1 * c1 - u2 * c2;
1802            u2 = u1 * c2 + u2 * c1;
1803            u1 = z;
1804        }
1805        c2 = TMath::Sqrt((1.0 - c1) / 2.0);
1806        if (dir == 1)
1807            c2 = -c2;
1808        c1 = TMath::Sqrt((1.0 + c1) / 2.0);
1809    }
1810
1811    /* Scaling for forward transform */
1812    if (dir == 1) {
1813        for (i=0;i<nn;i++) {
1814            x[i] /= (Double_t)nn;
1815            y[i] /= (Double_t)nn;
1816        }
1817    }
1818 }
1819
1820 //____________________________________________________________________
1821 void AliMultiplicityCorrection::GetComparisonResults(Float_t* mc, Int_t* mcLimit, Float_t* residuals, Float_t* ratioAverage) const
1822 {
1823   // Returns the chi2 between the MC and the unfolded ESD as well as between the ESD and the folded unfolded ESD
1824   // These values are computed during DrawComparison, thus this function picks up the
1825   // last calculation
1826
1827   if (mc)
1828     *mc = fLastChi2MC;
1829   if (mcLimit)
1830     *mcLimit = fLastChi2MCLimit;
1831   if (residuals)
1832     *residuals = fLastChi2Residuals;
1833   if (ratioAverage)
1834     *ratioAverage = fRatioAverage;
1835 }
1836
1837 //____________________________________________________________________
1838 TH2F* AliMultiplicityCorrection::GetMultiplicityMC(Int_t i, EventType eventType)
1839 {
1840   //
1841   // returns the corresponding MC spectrum
1842   //
1843
1844   switch (eventType)
1845   {
1846     case kTrVtx : return fMultiplicityVtx[i]; break;
1847     case kMB : return fMultiplicityMB[i]; break;
1848     case kINEL : return fMultiplicityINEL[i]; break;
1849   }
1850
1851   return 0;
1852 }
1853
1854 //____________________________________________________________________
1855 TH1* AliMultiplicityCorrection::CalculateStdDev(TH1** results, Int_t max)
1856 {
1857   // calculate standard deviation of (results[0] - results[k]) k=1...max-1
1858   // per bin one gets: sigma(r[0] - r[n]) / r[0]
1859
1860   TH1* standardDeviation = (TH1*) results[0]->Clone("standardDeviation");
1861   standardDeviation->Reset();
1862
1863   for (Int_t x=1; x<=results[0]->GetNbinsX(); x++)
1864   {
1865     if (results[0]->GetBinContent(x) > 0)
1866     {
1867       Double_t average = 0;
1868       for (Int_t n=1; n<max; ++n)
1869         average += results[n]->GetBinContent(x) - results[0]->GetBinContent(x);
1870       average /= max-1;
1871
1872       Double_t variance = 0;
1873       for (Int_t n=1; n<max; ++n)
1874       {
1875         Double_t value = results[n]->GetBinContent(x) - results[0]->GetBinContent(x) - average;
1876         variance += value * value;
1877       }
1878       variance /= max-1;
1879
1880       Double_t standardDev = TMath::Sqrt(variance);
1881       standardDeviation->SetBinContent(x, standardDev / results[0]->GetBinContent(x));
1882       //Printf("sigma_%d is %f value %f --> error %f", x, standardDev, results[0]->GetBinContent(x), standardDev / results[0]->GetBinContent(x));
1883     }
1884   }
1885
1886   return standardDeviation;
1887 }
1888
1889 //____________________________________________________________________
1890 TH1* AliMultiplicityCorrection::StatisticalUncertainty(MethodType methodType, Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t randomizeMeasured, Bool_t randomizeResponse, TH1* compareTo)
1891 {
1892   //
1893   // evaluates the uncertainty that arises from the non-infinite statistics in the response matrix
1894   // the function unfolds the spectrum using the default response matrix and several modified ones
1895   // the modified ones are created by randomizing each cell using poisson statistics with the mean = bin value
1896   // these unfolded results are compared to the first result gained with the default response OR to the histogram given
1897   // in <compareTo> (optional)
1898   //
1899   // returns the error assigned to the measurement
1900   //
1901
1902   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
1903
1904   // initialize seed with current time
1905   gRandom->SetSeed(0);
1906
1907   const Int_t kErrorIterations = 150;
1908
1909   TH1* maxError = 0;
1910   TH1* firstResult = 0;
1911
1912   TH1** results = new TH1*[kErrorIterations];
1913
1914   for (Int_t n=0; n<kErrorIterations; ++n)
1915   {
1916     Printf("Iteration %d of %d...", n, kErrorIterations);
1917
1918     // only done for chi2 minimization
1919     Bool_t createBigBin = kFALSE;
1920     if (methodType == kChi2Minimization)
1921       createBigBin = kTRUE;
1922
1923     SetupCurrentHists(inputRange, fullPhaseSpace, eventType, createBigBin);
1924
1925     TH1* measured = (TH1*) fCurrentESD->Clone("measured");
1926
1927     if (n > 0)
1928     {
1929       if (randomizeResponse)
1930       {
1931         // randomize response matrix
1932         for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
1933           for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
1934             fCurrentCorrelation->SetBinContent(i, j, gRandom->Poisson(fCurrentCorrelation->GetBinContent(i, j)));
1935       }
1936
1937       if (randomizeMeasured)
1938       {
1939         // randomize measured spectrum
1940         for (Int_t x=1; x<=measured->GetNbinsX(); x++) // mult. axis
1941         {
1942           Int_t randomValue = gRandom->Poisson(fCurrentESD->GetBinContent(x));
1943           measured->SetBinContent(x, randomValue);
1944           measured->SetBinError(x, TMath::Sqrt(randomValue));
1945         }
1946       }
1947     }
1948
1949     // only for bayesian method we have to do it before the call to Unfold...
1950     if (methodType == kBayesian)
1951     {
1952       for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
1953       {
1954         // with this it is normalized to 1
1955         Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
1956
1957         // with this normalized to the given efficiency
1958         if (fCurrentEfficiency->GetBinContent(i) > 0)
1959           sum /= fCurrentEfficiency->GetBinContent(i);
1960         else
1961           sum = 0;
1962
1963         for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
1964         {
1965           if (sum > 0)
1966           {
1967             fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
1968             fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
1969           }
1970           else
1971           {
1972             fCurrentCorrelation->SetBinContent(i, j, 0);
1973             fCurrentCorrelation->SetBinError(i, j, 0);
1974           }
1975         }
1976       }
1977     }
1978
1979     TH1* result = 0;
1980     if (n == 0 && compareTo)
1981     {
1982       // in this case we just store the histogram we want to compare to
1983       result = (TH1*) compareTo->Clone("compareTo");
1984       result->Sumw2();
1985     }
1986     else
1987     {
1988       result = (TH1*) fMultiplicityESDCorrected[correlationID]->Clone(Form("result_%d", n));
1989
1990       Int_t returnCode = -1;
1991       if (methodType == kChi2Minimization)
1992       {
1993         returnCode = UnfoldWithMinuit(fCurrentCorrelation, fCurrentEfficiency, measured, 0, result, kFALSE);
1994       }
1995       else if (methodType == kBayesian)
1996         returnCode = UnfoldWithBayesian(fCurrentCorrelation, fCurrentEfficiency, measured, 0, result, fgBayesianSmoothing, fgBayesianIterations);
1997
1998       if (returnCode != 0)
1999         return 0;
2000     }
2001
2002     // normalize
2003     result->Scale(1.0 / result->Integral());
2004
2005     if (n == 0)
2006     {
2007       firstResult = (TH1*) result->Clone("firstResult");
2008
2009       maxError = (TH1*) result->Clone("maxError");
2010       maxError->Reset();
2011     }
2012     else
2013     {
2014       // calculate ratio
2015       TH1* ratio = (TH1*) firstResult->Clone("ratio");
2016       ratio->Divide(result);
2017
2018       // find max. deviation
2019       for (Int_t x=1; x<=ratio->GetNbinsX(); x++)
2020         maxError->SetBinContent(x, TMath::Max(maxError->GetBinContent(x), TMath::Abs(1 - ratio->GetBinContent(x))));
2021
2022       delete ratio;
2023     }
2024
2025     results[n] = result;
2026   }
2027
2028   // find covariance matrix
2029   // results[n] is X_x
2030   // cov. matrix is M_xy = E ( (X_x - E(X_x)) * (X_y - E(X_y))), with E() = expectation value
2031
2032   Int_t nBins = results[0]->GetNbinsX();
2033   Float_t lowEdge = results[0]->GetXaxis()->GetBinLowEdge(1);
2034   Float_t upEdge = results[0]->GetXaxis()->GetBinUpEdge(nBins);
2035
2036   // find average, E(X)
2037   TProfile* average = new TProfile("average", "average", nBins, lowEdge, upEdge);
2038   for (Int_t n=1; n<kErrorIterations; ++n)
2039     for (Int_t x=1; x<=results[n]->GetNbinsX(); x++)
2040       average->Fill(results[n]->GetXaxis()->GetBinCenter(x), results[n]->GetBinContent(x));
2041   //new TCanvas; average->DrawClone();
2042
2043   // find cov. matrix
2044   TProfile2D* covMatrix = new TProfile2D("covMatrix", "covMatrix", nBins, lowEdge, upEdge, nBins, lowEdge, upEdge);
2045
2046   for (Int_t n=1; n<kErrorIterations; ++n)
2047     for (Int_t x=1; x<=results[n]->GetNbinsX(); x++)
2048       for (Int_t y=1; y<=results[n]->GetNbinsX(); y++)
2049       {
2050         // (X_x - E(X_x)) * (X_y - E(X_y)
2051         Float_t cov = (results[n]->GetBinContent(x) - average->GetBinContent(x)) * (results[n]->GetBinContent(y) - average->GetBinContent(y));
2052         covMatrix->Fill(results[n]->GetXaxis()->GetBinCenter(x), results[n]->GetXaxis()->GetBinCenter(y), cov);
2053       }
2054   new TCanvas; covMatrix->DrawClone("COLZ");
2055
2056 //   // fill 2D histogram that contains deviation from first
2057 //   TH2F* deviations = new TH2F("deviations", "deviations", nBins, lowEdge, upEdge, 1000, -0.01, 0.01);
2058 //   for (Int_t n=1; n<kErrorIterations; ++n)
2059 //     for (Int_t x=1; x<=results[n]->GetNbinsX(); x++)
2060 //       deviations->Fill(results[n]->GetXaxis()->GetBinCenter(x), results[n]->GetBinContent(x) - results[0]->GetBinContent(x));
2061 //   //new TCanvas; deviations->DrawCopy("COLZ");
2062 //
2063 //   // get standard deviation "by hand"
2064 //   for (Int_t x=1; x<=nBins; x++)
2065 //   {
2066 //     if (results[0]->GetBinContent(x) > 0)
2067 //     {
2068 //       TH1* proj = deviations->ProjectionY("projRMS", x, x, "e");
2069 //       Float_t standardDev = proj->GetRMS(); // this is standard deviation in fact
2070 //       //standardDeviation->SetBinContent(x, standardDev / results[0]->GetBinContent(x));
2071 //       Printf("sigma_%d is %f value %f --> error %f", x, standardDev, results[0]->GetBinContent(x), standardDev / results[0]->GetBinContent(x));
2072 //     }
2073 //   }
2074
2075   TH1* standardDeviation = CalculateStdDev(results, kErrorIterations);
2076
2077   // compare maxError to RMS of profile (variable name: average)
2078   // first: calculate rms in percent of value
2079   TH1* rmsError = (TH1*) maxError->Clone("rmsError");
2080   rmsError->Reset();
2081
2082   // enable error to be standard deviation (see http://root.cern.ch/root/html/TProfile.html#TProfile:SetErrorOption)
2083   average->SetErrorOption("s");
2084   for (Int_t x=1; x<=rmsError->GetNbinsX(); x++)
2085     if (average->GetBinContent(x) > 0)
2086       rmsError->SetBinContent(x, average->GetBinError(x) / average->GetBinContent(x));
2087
2088   // find maxError deviation from average (not from "first result"/mc as above)
2089   TH1* maxError2 = (TH1*) maxError->Clone("maxError2");
2090   maxError2->Reset();
2091   for (Int_t n=1; n<kErrorIterations; ++n)
2092     for (Int_t x=1; x<=results[n]->GetNbinsX(); x++)
2093       if (average->GetBinContent(x) > 0)
2094         maxError2->SetBinContent(x, TMath::Max(maxError2->GetBinContent(x), TMath::Abs((results[n]->GetBinContent(x) - average->GetBinContent(x)) / average->GetBinContent(x))));
2095
2096   //new TCanvas; maxError2->DrawCopy(); rmsError->SetLineColor(2); rmsError->DrawCopy("SAME"); standardDeviation->SetLineColor(3); standardDeviation->DrawCopy("SAME");
2097
2098   // plot difference between average and MC/first result
2099   TH1* averageFirstRatio = (TH1*) results[0]->Clone("averageFirstRatio");
2100   averageFirstRatio->Reset();
2101   averageFirstRatio->Divide(results[0], average);
2102
2103   //new TCanvas; results[0]->DrawCopy(); average->SetLineColor(2); average->DrawClone("SAME");
2104   //new TCanvas; averageFirstRatio->DrawCopy();
2105
2106   // clean up
2107   for (Int_t n=0; n<kErrorIterations; ++n)
2108     delete results[n];
2109   delete[] results;
2110
2111   // fill into result histogram
2112   for (Int_t i=1; i<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); ++i)
2113     fMultiplicityESDCorrected[correlationID]->SetBinContent(i, firstResult->GetBinContent(i));
2114
2115   for (Int_t i=1; i<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); ++i)
2116     fMultiplicityESDCorrected[correlationID]->SetBinError(i, maxError->GetBinContent(i) * fMultiplicityESDCorrected[correlationID]->GetBinContent(i));
2117
2118   return standardDeviation;
2119 }
2120
2121 //____________________________________________________________________
2122 void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Float_t regPar, Int_t nIterations, TH1* initialConditions, Bool_t determineError)
2123 {
2124   //
2125   // correct spectrum using bayesian method
2126   //
2127
2128   // initialize seed with current time
2129   gRandom->SetSeed(0);
2130
2131   SetupCurrentHists(inputRange, fullPhaseSpace, eventType, kFALSE);
2132
2133   // normalize correction for given nPart
2134   for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
2135   {
2136     // with this it is normalized to 1
2137     Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
2138
2139     // with this normalized to the given efficiency
2140     if (fCurrentEfficiency->GetBinContent(i) > 0)
2141       sum /= fCurrentEfficiency->GetBinContent(i);
2142     else
2143       sum = 0;
2144
2145     for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
2146     {
2147       if (sum > 0)
2148       {
2149         fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
2150         fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
2151       }
2152       else
2153       {
2154         fCurrentCorrelation->SetBinContent(i, j, 0);
2155         fCurrentCorrelation->SetBinError(i, j, 0);
2156       }
2157     }
2158   }
2159
2160   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
2161
2162   if (UnfoldWithBayesian(fCurrentCorrelation, fCurrentEfficiency, fCurrentESD, initialConditions, fMultiplicityESDCorrected[correlationID], regPar, nIterations) != 0)
2163     return;
2164
2165   if (!determineError)
2166   {
2167     Printf("AliMultiplicityCorrection::ApplyBayesianMethod: WARNING: No errors calculated.");
2168     return;
2169   }
2170
2171   // evaluate errors, this is done by randomizing the measured spectrum following Poission statistics
2172   // this (new) measured spectrum is then unfolded and the different to the result from the "real" measured
2173   // spectrum calculated. This is performed N times and the sigma is taken as the statistical
2174   // error of the unfolding method itself.
2175
2176   const Int_t kErrorIterations = 20;
2177
2178   printf("Spectrum unfolded. Determining error (%d iterations)...\n", kErrorIterations);
2179
2180   TH1* randomized = (TH1*) fCurrentESD->Clone("randomized");
2181   TH1* resultArray[kErrorIterations+1];
2182   for (Int_t n=0; n<kErrorIterations; ++n)
2183   {
2184     // randomize the content of clone following a poisson with the mean = the value of that bin
2185     for (Int_t x=1; x<=randomized->GetNbinsX(); x++) // mult. axis
2186     {
2187       Int_t randomValue = gRandom->Poisson(fCurrentESD->GetBinContent(x));
2188       //printf("%d --> %d\n", fCurrentESD->GetBinContent(x), randomValue);
2189       randomized->SetBinContent(x, randomValue);
2190       randomized->SetBinError(x, TMath::Sqrt(randomValue));
2191     }
2192
2193     TH1* result2 = (TH1*) fMultiplicityESDCorrected[correlationID]->Clone("result2");
2194     result2->Reset();
2195     if (UnfoldWithBayesian(fCurrentCorrelation, fCurrentEfficiency, randomized, initialConditions, result2, regPar, nIterations) != 0)
2196       return;
2197
2198     result2->Scale(1.0 / result2->Integral());
2199
2200     resultArray[n+1] = result2;
2201   }
2202   delete randomized;
2203
2204   resultArray[0] = fMultiplicityESDCorrected[correlationID];
2205   TH1* error = CalculateStdDev(resultArray, kErrorIterations+1);
2206
2207   for (Int_t n=0; n<kErrorIterations; ++n)
2208     delete resultArray[n+1];
2209
2210   for (Int_t i=1; i<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); ++i)
2211     fMultiplicityESDCorrected[correlationID]->SetBinError(i, error->GetBinContent(i) * fMultiplicityESDCorrected[correlationID]->GetBinContent(i));
2212
2213   delete error;
2214 }
2215
2216 //____________________________________________________________________
2217 Int_t AliMultiplicityCorrection::UnfoldWithBayesian(TH1* correlation, TH1* aEfficiency, TH1* measured, TH1* initialConditions, TH1* aResult, Float_t regPar, Int_t nIterations)
2218 {
2219   //
2220   // unfolds a spectrum
2221   //
2222
2223   if (measured->Integral() <= 0)
2224   {
2225     Printf("AliMultiplicityCorrection::UnfoldWithBayesian: ERROR: The measured spectrum is empty");
2226     return -1;
2227   }
2228
2229   const Int_t kStartBin = 0;
2230
2231   const Int_t kMaxM = fgkMaxInput;  //<= fCurrentCorrelation->GetNbinsY(); // max measured axis
2232   const Int_t kMaxT = fgkMaxParams; //<= fCurrentCorrelation->GetNbinsX(); // max true axis
2233
2234   // convergence limit: kMaxT * 0.001^2 = kMaxT * 1e-6 (e.g. 250 bins --> 2.5 e-4)
2235   const Double_t kConvergenceLimit = kMaxT * 1e-6;
2236
2237   // store information in arrays, to increase processing speed (~ factor 5)
2238   Double_t measuredCopy[kMaxM];
2239   Double_t measuredError[kMaxM];
2240   Double_t prior[kMaxT];
2241   Double_t result[kMaxT];
2242   Double_t efficiency[kMaxT];
2243
2244   Double_t response[kMaxT][kMaxM];
2245   Double_t inverseResponse[kMaxT][kMaxM];
2246
2247   // for normalization
2248   Float_t measuredIntegral = measured->Integral();
2249   for (Int_t m=0; m<kMaxM; m++)
2250   {
2251     measuredCopy[m] = measured->GetBinContent(m+1) / measuredIntegral;
2252     measuredError[m] = measured->GetBinError(m+1) / measuredIntegral;
2253
2254     for (Int_t t=0; t<kMaxT; t++)
2255     {
2256       response[t][m] = correlation->GetBinContent(t+1, m+1);
2257       inverseResponse[t][m] = 0;
2258     }
2259   }
2260
2261   for (Int_t t=0; t<kMaxT; t++)
2262   {
2263     efficiency[t] = aEfficiency->GetBinContent(t+1);
2264     prior[t] = measuredCopy[t];
2265     result[t] = 0;
2266   }
2267
2268   // pick prior distribution
2269   if (initialConditions)
2270   {
2271     printf("Using different starting conditions...\n");
2272     // for normalization
2273     Float_t inputDistIntegral = initialConditions->Integral();
2274     for (Int_t i=0; i<kMaxT; i++)
2275       prior[i] = initialConditions->GetBinContent(i+1) / inputDistIntegral;
2276   }
2277
2278   //TH1F* convergence = new TH1F("convergence", "convergence", 200, 0.5, 200.5);
2279
2280   // unfold...
2281   for (Int_t i=0; i<nIterations || nIterations < 0; i++)
2282   {
2283     //printf(" iteration %i \n", i);
2284
2285     // calculate IR from Beyes theorem
2286     // IR_ji = R_ij * prior_i / sum_k(R_kj * prior_k)
2287
2288     Double_t chi2Measured = 0;
2289     for (Int_t m=0; m<kMaxM; m++)
2290     {
2291       Float_t norm = 0;
2292       for (Int_t t = kStartBin; t<kMaxT; t++)
2293         norm += response[t][m] * prior[t];
2294
2295       // calc. chi2: (measured - response * prior) / error
2296       if (measuredError[m] > 0)
2297       {
2298         Double_t value = (measuredCopy[m] - norm) / measuredError[m];
2299         chi2Measured += value * value;
2300       }
2301
2302       if (norm > 0)
2303       {
2304         for (Int_t t = kStartBin; t<kMaxT; t++)
2305           inverseResponse[t][m] = response[t][m] * prior[t] / norm;
2306       }
2307       else
2308       {
2309         for (Int_t t = kStartBin; t<kMaxT; t++)
2310           inverseResponse[t][m] = 0;
2311       }
2312     }
2313     //Printf("chi2Measured of the last prior is %e", chi2Measured);
2314
2315     for (Int_t t = kStartBin; t<kMaxT; t++)
2316     {
2317       Float_t value = 0;
2318       for (Int_t m=0; m<kMaxM; m++)
2319         value += inverseResponse[t][m] * measuredCopy[m];
2320
2321       if (efficiency[t] > 0)
2322         result[t] = value / efficiency[t];
2323       else
2324         result[t] = 0;
2325     }
2326
2327     Double_t chi2LastIter = 0;
2328     // regularization (simple smoothing)
2329     for (Int_t t=kStartBin; t<kMaxT; t++)
2330     {
2331       Float_t newValue = 0;
2332       // 0 bin excluded from smoothing
2333       if (t > kStartBin+1 && t<kMaxT-1)
2334       {
2335         Float_t average = (result[t-1] + result[t] + result[t+1]) / 3;
2336
2337         // weight the average with the regularization parameter
2338         newValue = (1 - regPar) * result[t] + regPar * average;
2339       }
2340       else
2341         newValue = result[t];
2342
2343       // calculate chi2 (change from last iteration)
2344       if (prior[t] > 1e-5)
2345       {
2346         Double_t diff = (prior[t] - newValue) / prior[t];
2347         chi2LastIter += diff * diff;
2348       }
2349
2350       prior[t] = newValue;
2351     }
2352     //printf("Chi2 of %d iteration = %e\n", i, chi2LastIter);
2353     //convergence->Fill(i+1, chi2LastIter);
2354
2355     if (nIterations < 0 && chi2LastIter < kConvergenceLimit)
2356     {
2357       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);
2358       break;
2359     }
2360
2361     //if (i % 5 == 0)
2362     //  DrawComparison(Form("Bayesian_%d", i), mcTarget, correlationID, kTRUE, eventType);
2363   } // end of iterations
2364
2365   //new TCanvas; convergence->DrawCopy(); gPad->SetLogy();
2366   //delete convergence;
2367
2368   for (Int_t t=0; t<kMaxT; t++)
2369     aResult->SetBinContent(t+1, result[t]);
2370
2371   return 0;
2372
2373   // ********
2374   // Calculate the covariance matrix, all arguments are taken from NIM,A362,487-498,1995
2375
2376   /*printf("Calculating covariance matrix. This may take some time...\n");
2377
2378   // TODO check if this is the right one...
2379   TH1* sumHist = GetMultiplicityMC(inputRange, eventType)->ProjectionY("sumHist", 1, GetMultiplicityMC(inputRange, eventType)->GetNbinsX());
2380
2381   Int_t xBins = hInverseResponseBayes->GetNbinsX();
2382   Int_t yBins = hInverseResponseBayes->GetNbinsY();
2383
2384   // calculate "unfolding matrix" Mij
2385   Float_t matrixM[251][251];
2386   for (Int_t i=1; i<=xBins; i++)
2387   {
2388     for (Int_t j=1; j<=yBins; j++)
2389     {
2390       if (fCurrentEfficiency->GetBinContent(i) > 0)
2391         matrixM[i-1][j-1] = hInverseResponseBayes->GetBinContent(i, j) / fCurrentEfficiency->GetBinContent(i);
2392       else
2393         matrixM[i-1][j-1] = 0;
2394     }
2395   }
2396
2397   Float_t* vectorn = new Float_t[yBins];
2398   for (Int_t j=1; j<=yBins; j++)
2399     vectorn[j-1] = fCurrentESD->GetBinContent(j);
2400
2401   // first part of covariance matrix, depends on input distribution n(E)
2402   Float_t cov1[251][251];
2403
2404   Float_t nEvents = fCurrentESD->Integral(); // N
2405
2406   xBins = 20;
2407   yBins = 20;
2408
2409   for (Int_t k=0; k<xBins; k++)
2410   {
2411     printf("In Cov1: %d\n", k);
2412     for (Int_t l=0; l<yBins; l++)
2413     {
2414       cov1[k][l] = 0;
2415
2416       // sum_j Mkj Mlj n(Ej) * (1 - n(Ej) / N)
2417       for (Int_t j=0; j<yBins; j++)
2418         cov1[k][l] += matrixM[k][j] * matrixM[l][j] * vectorn[j]
2419           * (1.0 - vectorn[j] / nEvents);
2420
2421       // - sum_i,j (i != j) Mki Mlj n(Ei) n(Ej) / N
2422       for (Int_t i=0; i<yBins; i++)
2423         for (Int_t j=0; j<yBins; j++)
2424         {
2425           if (i == j)
2426             continue;
2427           cov1[k][l] -= matrixM[k][i] * matrixM[l][j] * vectorn[i]
2428             * vectorn[j] / nEvents;
2429          }
2430     }
2431   }
2432
2433   printf("Cov1 finished\n");
2434
2435   TH2F* cov = (TH2F*) hInverseResponseBayes->Clone("cov");
2436   cov->Reset();
2437
2438   for (Int_t i=1; i<=xBins; i++)
2439     for (Int_t j=1; j<=yBins; j++)
2440       cov->SetBinContent(i, j, cov1[i-1][j-1]);
2441
2442   new TCanvas;
2443   cov->Draw("COLZ");
2444
2445   // second part of covariance matrix, depends on response matrix
2446   Float_t cov2[251][251];
2447
2448   // Cov[P(Er|Cu), P(Es|Cu)] term
2449   Float_t covTerm[100][100][100];
2450   for (Int_t r=0; r<yBins; r++)
2451     for (Int_t u=0; u<xBins; u++)
2452       for (Int_t s=0; s<yBins; s++)
2453       {
2454         if (r == s)
2455           covTerm[r][u][s] = 1.0 / sumHist->GetBinContent(u+1) * hResponse->GetBinContent(u+1, r+1)
2456             * (1.0 - hResponse->GetBinContent(u+1, r+1));
2457         else
2458           covTerm[r][u][s] = - 1.0 / sumHist->GetBinContent(u+1) * hResponse->GetBinContent(u+1, r+1)
2459             * hResponse->GetBinContent(u+1, s+1);
2460       }
2461
2462   for (Int_t k=0; k<xBins; k++)
2463     for (Int_t l=0; l<yBins; l++)
2464     {
2465       cov2[k][l] = 0;
2466       printf("In Cov2: %d %d\n", k, l);
2467       for (Int_t i=0; i<yBins; i++)
2468         for (Int_t j=0; j<yBins; j++)
2469         {
2470           //printf("In Cov2: %d %d %d %d\n", k, l, i, j);
2471           // calculate Cov(Mki, Mlj) = sum{ru},{su} ...
2472           Float_t tmpCov = 0;
2473           for (Int_t r=0; r<yBins; r++)
2474             for (Int_t u=0; u<xBins; u++)
2475               for (Int_t s=0; s<yBins; s++)
2476               {
2477                 if (hResponse->GetBinContent(u+1, r+1) == 0 || hResponse->GetBinContent(u+1, s+1) == 0
2478                   || hResponse->GetBinContent(u+1, i+1) == 0)
2479                   continue;
2480
2481                 tmpCov += BayesCovarianceDerivate(matrixM, hResponse, fCurrentEfficiency, k, i, r, u)
2482                   * BayesCovarianceDerivate(matrixM, hResponse, fCurrentEfficiency, l, j, s, u)
2483                   * covTerm[r][u][s];
2484               }
2485
2486           cov2[k][l] += fCurrentESD->GetBinContent(i+1) * fCurrentESD->GetBinContent(j+1) * tmpCov;
2487         }
2488     }
2489
2490   printf("Cov2 finished\n");
2491
2492   for (Int_t i=1; i<=xBins; i++)
2493     for (Int_t j=1; j<=yBins; j++)
2494       cov->SetBinContent(i, j, cov1[i-1][j-1] + cov2[i-1][j-1]);
2495
2496   new TCanvas;
2497   cov->Draw("COLZ");*/
2498 }
2499
2500 //____________________________________________________________________
2501 Float_t AliMultiplicityCorrection::BayesCovarianceDerivate(Float_t matrixM[251][251], TH2* hResponse, Int_t k, Int_t i, Int_t r, Int_t u)
2502 {
2503   //
2504   // helper function for the covariance matrix of the bayesian method
2505   //
2506
2507   Float_t result = 0;
2508
2509   if (k == u && r == i)
2510     result += 1.0 / hResponse->GetBinContent(u+1, r+1);
2511
2512   if (k == u)
2513     result -= 1.0 / fCurrentEfficiency->GetBinContent(u+1);
2514
2515   if (r == i)
2516     result -= matrixM[u][i] * fCurrentEfficiency->GetBinContent(u+1) / hResponse->GetBinContent(u+1, i+1);
2517
2518   result *= matrixM[k][i];
2519
2520   return result;
2521 }
2522
2523 //____________________________________________________________________
2524 void AliMultiplicityCorrection::ApplyLaszloMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType)
2525 {
2526   //
2527   // correct spectrum using bayesian method
2528   //
2529
2530   Float_t regPar = 0;
2531
2532   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
2533   Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4);
2534
2535   SetupCurrentHists(inputRange, fullPhaseSpace, eventType, kFALSE);
2536   //normalize ESD
2537   fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
2538
2539   // TODO should be taken from correlation map
2540   //TH1* sumHist = GetMultiplicityMC(inputRange, eventType)->ProjectionY("sumHist", 1, GetMultiplicityMC(inputRange, eventType)->GetNbinsX());
2541
2542   // normalize correction for given nPart
2543   for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
2544   {
2545     Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
2546     //Double_t sum = sumHist->GetBinContent(i);
2547     if (sum <= 0)
2548       continue;
2549     for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
2550     {
2551       // npart sum to 1
2552       fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
2553       fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
2554     }
2555   }
2556
2557   new TCanvas;
2558   fCurrentCorrelation->Draw("COLZ");
2559
2560   // FAKE
2561   fCurrentEfficiency = ((TH2*) fCurrentCorrelation)->ProjectionX("eff");
2562
2563   // pick prior distribution
2564   TH1F* hPrior = (TH1F*)fCurrentESD->Clone("prior");
2565   Float_t norm = 1; //hPrior->Integral();
2566   for (Int_t t=1; t<=hPrior->GetNbinsX(); t++)
2567     hPrior->SetBinContent(t, hPrior->GetBinContent(t)/norm);
2568
2569   // zero distribution
2570   TH1F* zero =  (TH1F*)hPrior->Clone("zero");
2571
2572   // define temp hist
2573   TH1F* hTemp = (TH1F*)fCurrentESD->Clone("temp");
2574   hTemp->Reset();
2575
2576   // just a shortcut
2577   TH2F* hResponse = (TH2F*) fCurrentCorrelation;
2578
2579   // unfold...
2580   Int_t iterations = 25;
2581   for (Int_t i=0; i<iterations; i++)
2582   {
2583     //printf(" iteration %i \n", i);
2584
2585     for (Int_t m=1; m<=hResponse->GetNbinsY(); m++)
2586     {
2587       Float_t value = 0;
2588       for (Int_t t = 1; t<=hResponse->GetNbinsX(); t++)
2589         value += hResponse->GetBinContent(t, m) * hPrior->GetBinContent(t);
2590       hTemp->SetBinContent(m, value);
2591       //printf("%d %f %f %f\n", m, zero->GetBinContent(m), hPrior->GetBinContent(m), value);
2592     }
2593
2594     // regularization (simple smoothing)
2595     TH1F* hTrueSmoothed = (TH1F*) hTemp->Clone("truesmoothed");
2596
2597     for (Int_t t=2; t<hTrueSmoothed->GetNbinsX(); t++)
2598     {
2599       Float_t average = (hTemp->GetBinContent(t-1) / hTemp->GetBinWidth(t-1)
2600                          + hTemp->GetBinContent(t) / hTemp->GetBinWidth(t)
2601                          + hTemp->GetBinContent(t+1) / hTemp->GetBinWidth(t+1)) / 3.;
2602       average *= hTrueSmoothed->GetBinWidth(t);
2603
2604       // weight the average with the regularization parameter
2605       hTrueSmoothed->SetBinContent(t, (1-regPar) * hTemp->GetBinContent(t) + regPar * average);
2606     }
2607
2608     for (Int_t m=1; m<=hResponse->GetNbinsY(); m++)
2609       hTemp->SetBinContent(m, zero->GetBinContent(m) + hPrior->GetBinContent(m) - hTrueSmoothed->GetBinContent(m));
2610
2611     // fill guess
2612     for (Int_t t=1; t<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); t++)
2613     {
2614       fMultiplicityESDCorrected[correlationID]->SetBinContent(t, hTemp->GetBinContent(t));
2615       fMultiplicityESDCorrected[correlationID]->SetBinError(t, 0.05 * hTemp->GetBinContent(t)); // TODO
2616
2617       //printf(" bin %d content %f \n", t, fMultiplicityESDCorrected[correlationID]->GetBinContent(t));
2618     }
2619
2620
2621     // calculate chi2 (change from last iteration)
2622     Double_t chi2 = 0;
2623
2624     // use smoothed true (normalized) as new prior
2625     norm = 1; //hTrueSmoothed->Integral();
2626
2627     for (Int_t t=1; t<hTrueSmoothed->GetNbinsX(); t++)
2628     {
2629       Float_t newValue = hTemp->GetBinContent(t)/norm;
2630       Float_t diff = hPrior->GetBinContent(t) - newValue;
2631       chi2 += (Double_t) diff * diff;
2632
2633       hPrior->SetBinContent(t, newValue);
2634     }
2635
2636     printf("Chi2 of %d iteration = %.10f\n", i, chi2);
2637
2638     //if (i % 5 == 0)
2639       DrawComparison(Form("Laszlo_%d", i), inputRange, fullPhaseSpace, kTRUE, GetMultiplicityMC(mcTarget, eventType)->ProjectionY());
2640
2641     delete hTrueSmoothed;
2642   } // end of iterations
2643
2644   DrawComparison("Laszlo", inputRange, fullPhaseSpace, kTRUE, GetMultiplicityMC(mcTarget, eventType)->ProjectionY());
2645 }
2646
2647 //____________________________________________________________________
2648 void AliMultiplicityCorrection::ApplyGaussianMethod(Int_t inputRange, Bool_t fullPhaseSpace)
2649 {
2650   //
2651   // correct spectrum using a simple Gaussian approach, that is model-dependent
2652   //
2653
2654   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
2655   Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4);
2656
2657   SetupCurrentHists(inputRange, fullPhaseSpace, kTrVtx, kFALSE);
2658   //normalize ESD
2659   fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
2660
2661   TH1D* correction = dynamic_cast<TH1D*> (fCurrentESD->Clone("GaussianMean"));
2662   correction->SetTitle("GaussianMean");
2663
2664   TH1D* correctionWidth = dynamic_cast<TH1D*> (fCurrentESD->Clone("GaussianWidth"));
2665   correction->SetTitle("GaussianWidth");
2666
2667   for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
2668   {
2669     TH1D* proj = (dynamic_cast<TH2*> (fCurrentCorrelation))->ProjectionX("_px", i, i+1);
2670     proj->Fit("gaus", "0Q");
2671     correction->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(1));
2672     correctionWidth->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(2));
2673
2674     continue;
2675
2676     // draw for debugging
2677     new TCanvas;
2678     proj->DrawCopy();
2679     proj->GetFunction("gaus")->DrawCopy("SAME");
2680   }
2681
2682   TH1* target = fMultiplicityESDCorrected[correlationID];
2683
2684   Int_t nBins = target->GetNbinsX()*10+1;
2685   Float_t* binning = new Float_t[nBins];
2686   for (Int_t i=1; i<=target->GetNbinsX(); ++i)
2687     for (Int_t j=0; j<10; ++j)
2688       binning[(i-1)*10 + j] = target->GetXaxis()->GetBinLowEdge(i) + target->GetXaxis()->GetBinWidth(i) / 10 * j;
2689
2690   binning[nBins-1] = target->GetXaxis()->GetBinUpEdge(target->GetNbinsX());
2691
2692   TH1F* fineBinned = new TH1F("targetFineBinned", "targetFineBinned", nBins-1, binning);
2693
2694   for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
2695   {
2696     Float_t mean = correction->GetBinContent(i);
2697     Float_t width = correctionWidth->GetBinContent(i);
2698
2699     Int_t fillBegin = fineBinned->FindBin(mean - width * 5);
2700     Int_t fillEnd   = fineBinned->FindBin(mean + width * 5);
2701     //printf("bin %d mean %f width %f, filling from %d to %d\n", i, mean, width, fillBegin, fillEnd);
2702
2703     for (Int_t j=fillBegin; j <= fillEnd; ++j)
2704     {
2705       fineBinned->AddBinContent(j, TMath::Gaus(fineBinned->GetXaxis()->GetBinCenter(j), mean, width, kTRUE) * fCurrentESD->GetBinContent(i));
2706     }
2707   }
2708
2709   for (Int_t i=1; i<=target->GetNbinsX(); ++i)
2710   {
2711     Float_t sum = 0;
2712     for (Int_t j=1; j<=10; ++j)
2713       sum += fineBinned->GetBinContent((i-1)*10 + j);
2714     target->SetBinContent(i, sum / 10);
2715   }
2716
2717   delete[] binning;
2718
2719   DrawComparison("Gaussian", inputRange, fullPhaseSpace, kFALSE, GetMultiplicityMC(mcTarget, kTrVtx)->ProjectionY());
2720 }
2721
2722 //____________________________________________________________________
2723 TH2F* AliMultiplicityCorrection::CalculateMultiplicityESD(TH1* inputMC, Int_t correlationMap)
2724 {
2725   // runs the distribution given in inputMC through the response matrix identified by
2726   // correlationMap and produces a measured distribution
2727   // although it is a TH2F the vertex axis is not used at the moment and all entries are filled in mid-vertex
2728   // if normalized is set, inputMC is expected to be normalized to the bin width
2729
2730   if (!inputMC)
2731     return 0;
2732
2733   if (correlationMap < 0 || correlationMap >= kCorrHists)
2734     return 0;
2735
2736   // empty under/overflow bins in x, otherwise Project3D takes them into account
2737   TH3* hist = fCorrelation[correlationMap];
2738   for (Int_t y=1; y<=hist->GetYaxis()->GetNbins(); ++y)
2739   {
2740     for (Int_t z=1; z<=hist->GetZaxis()->GetNbins(); ++z)
2741     {
2742       hist->SetBinContent(0, y, z, 0);
2743       hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0);
2744     }
2745   }
2746
2747   TH2* corr = (TH2*) hist->Project3D("zy");
2748   //corr->Rebin2D(2, 1);
2749   corr->Sumw2();
2750
2751   // normalize correction for given nPart
2752   for (Int_t i=1; i<=corr->GetNbinsX(); ++i)
2753   {
2754     Double_t sum = corr->Integral(i, i, 1, corr->GetNbinsY());
2755     if (sum <= 0)
2756       continue;
2757
2758     for (Int_t j=1; j<=corr->GetNbinsY(); ++j)
2759     {
2760       // npart sum to 1
2761       corr->SetBinContent(i, j, corr->GetBinContent(i, j) / sum);
2762       corr->SetBinError(i, j, corr->GetBinError(i, j) / sum);
2763     }
2764   }
2765
2766   TH2F* target = dynamic_cast<TH2F*> (fMultiplicityESD[0]->Clone(Form("%s_measured", inputMC->GetName())));
2767   target->Reset();
2768
2769   for (Int_t meas=1; meas<=corr->GetNbinsY(); ++meas)
2770   {
2771     Float_t measured = 0;
2772     Float_t error = 0;
2773
2774     for (Int_t gen=1; gen<=corr->GetNbinsX(); ++gen)
2775     {
2776       Int_t mcGenBin = inputMC->GetXaxis()->FindBin(corr->GetXaxis()->GetBinCenter(gen));
2777
2778       measured += inputMC->GetBinContent(mcGenBin) * corr->GetBinContent(gen, meas);
2779       error += inputMC->GetBinError(mcGenBin) * corr->GetBinContent(gen, meas);
2780     }
2781
2782     //printf("%f +- %f ; %f +- %f \n", inputMC->GetBinContent(meas), inputMC->GetBinError(meas), measured, error);
2783
2784     target->SetBinContent(1 + target->GetNbinsX() / 2, meas, measured);
2785     target->SetBinError(1 + target->GetNbinsX() / 2, meas, error);
2786   }
2787
2788   return target;
2789 }
2790
2791 //____________________________________________________________________
2792 void AliMultiplicityCorrection::SetGenMeasFromFunc(TF1* inputMC, Int_t id)
2793 {
2794   // uses the given function to fill the input MC histogram and generates from that
2795   // the measured histogram by applying the response matrix
2796   // this can be used to evaluate if the methods work indepedently of the input
2797   // distribution
2798   // WARNING does not respect the vertex distribution, just fills central vertex bin
2799
2800   if (!inputMC)
2801     return;
2802
2803   if (id < 0 || id >= kESDHists)
2804     return;
2805
2806   TH2F* mc = fMultiplicityVtx[id];
2807
2808   mc->Reset();
2809
2810   for (Int_t i=1; i<=mc->GetNbinsY(); ++i)
2811   {
2812     mc->SetBinContent(mc->GetNbinsX() / 2 + 1, i, inputMC->Eval(mc->GetYaxis()->GetBinCenter(i)) * mc->GetYaxis()->GetBinWidth(i));
2813     mc->SetBinError(mc->GetNbinsX() / 2 + 1, i, 0);
2814   }
2815
2816   //new TCanvas;
2817   //mc->Draw("COLZ");
2818
2819   TH1* proj = mc->ProjectionY();
2820
2821   TString tmp(fMultiplicityESD[id]->GetName());
2822   delete fMultiplicityESD[id];
2823   fMultiplicityESD[id] = CalculateMultiplicityESD(proj, id);
2824   fMultiplicityESD[id]->SetName(tmp);
2825 }