PWG0/multiplicity/AliMultiplicityCorrection.cxx

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