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