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 <TCanvas.h>
  33 #include <TString.h>
  34 #include <TF1.h>
  35 #include <TMath.h>
  36 #include <TCollection.h>
  37 #include <TLegend.h>
  38 #include <TLine.h>
  39 #include <TRandom.h>
  40 #include <TProfile.h>
  41 #include <TProfile2D.h>
  42
  43 ClassImp(AliMultiplicityCorrection)
  44
  45 // These are the areas where the quality of the unfolding results are evaluated
  46 // Default defined here, call SetQualityRegions to change them
  47 // unit is in multiplicity (not in bin!)
  48
  49 // SPD:   peak area - flat area - low stat area
  50 Int_t AliMultiplicityCorrection::fgQualityRegionsB[kQualityRegions] = {1,  20, 70};
  51 Int_t AliMultiplicityCorrection::fgQualityRegionsE[kQualityRegions] = {10, 65, 80};
  52
  53 //____________________________________________________________________
  54 void AliMultiplicityCorrection::SetQualityRegions(Bool_t SPDStudy)
  55 {
  56   //
  57   // sets the quality region definition to TPC or SPD
  58   //
  59
  60   if (SPDStudy)
  61   {
  62     // SPD:   peak area - flat area - low stat area
  63     fgQualityRegionsB[0] = 1;
  64     fgQualityRegionsE[0] = 10;
  65
  66     fgQualityRegionsB[1] = 20;
  67     fgQualityRegionsE[1] = 65;
  68
  69     fgQualityRegionsB[2] = 70;
  70     fgQualityRegionsE[2] = 80;
  71
  72     Printf("AliMultiplicityCorrection::SetQualityRegions --> Enabled quality regions for SPD");
  73   }
  74   else
  75   {
  76     // TPC:   peak area - flat area - low stat area
  77     fgQualityRegionsB[0] = 4;
  78     fgQualityRegionsE[0] = 12;
  79
  80     fgQualityRegionsB[1] = 25;
  81     fgQualityRegionsE[1] = 55;
  82
  83     fgQualityRegionsB[2] = 88;
  84     fgQualityRegionsE[2] = 108;
  85
  86     Printf("AliMultiplicityCorrection::SetQualityRegions --> Enabled quality regions for TPC");
  87   }
  88 }
  89
  90 //____________________________________________________________________
  91 AliMultiplicityCorrection::AliMultiplicityCorrection() :
  92   TNamed(), fCurrentESD(0), fCurrentCorrelation(0), fCurrentEfficiency(0), fLastBinLimit(0), fLastChi2MC(0), fLastChi2MCLimit(0), fLastChi2Residuals(0), fRatioAverage(0)
  93 {
  94   //
  95   // default constructor
  96   //
  97
  98   for (Int_t i = 0; i < kESDHists; ++i)
  99     fMultiplicityESD[i] = 0;
 100
 101   for (Int_t i = 0; i < kMCHists; ++i)
 102   {
 103     fMultiplicityVtx[i] = 0;
 104     fMultiplicityMB[i] = 0;
 105     fMultiplicityINEL[i] = 0;
 106     fMultiplicityNSD[i] = 0;
 107   }
 108
 109   for (Int_t i = 0; i < kCorrHists; ++i)
 110   {
 111     fCorrelation[i] = 0;
 112     fMultiplicityESDCorrected[i] = 0;
 113   }
 114
 115   for (Int_t i = 0; i < kQualityRegions; ++i)
 116     fQuality[i] = 0;
 117 }
 118
 119 //____________________________________________________________________
 120 AliMultiplicityCorrection* AliMultiplicityCorrection::Open(const char* fileName, const char* folderName)
 121 {
 122   // opens the given file, reads the multiplicity from the given folder and returns the object
 123
 124   TFile* file = TFile::Open(fileName);
 125   if (!file)
 126   {
 127     Printf("ERROR: Could not open %s", fileName);
 128     return 0;
 129   }
 130
 131   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection(folderName, folderName);
 132   mult->LoadHistograms();
 133
 134   // TODO closing the file does not work here, because the histograms cannot be read anymore. LoadHistograms need to be adapted
 135
 136   return mult;
 137 }
 138
 139 //____________________________________________________________________
 140 AliMultiplicityCorrection::AliMultiplicityCorrection(const Char_t* name, const Char_t* title) :
 141   TNamed(name, title),
 142   fCurrentESD(0),
 143   fCurrentCorrelation(0),
 144   fCurrentEfficiency(0),
 145   fLastBinLimit(0),
 146   fLastChi2MC(0),
 147   fLastChi2MCLimit(0),
 148   fLastChi2Residuals(0),
 149   fRatioAverage(0)
 150 {
 151   //
 152   // named constructor
 153   //
 154
 155   // do not add this hists to the directory
 156   Bool_t oldStatus = TH1::AddDirectoryStatus();
 157   TH1::AddDirectory(kFALSE);
 158
 159   /*Float_t binLimitsVtx[] = {-10,-9,-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8,9,10};
 160   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,
 161                           10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5,
 162                           20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5,
 163                           30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5,
 164                           40.5, 41.5, 42.5, 43.5, 44.5, 45.5, 46.5, 47.5, 48.5, 49.5,
 165                           50.5, 55.5, 60.5, 65.5, 70.5, 75.5, 80.5, 85.5, 90.5, 95.5,
 166                           100.5, 105.5, 110.5, 115.5, 120.5, 125.5, 130.5, 135.5, 140.5, 145.5,
 167                           150.5, 160.5, 170.5, 180.5, 190.5, 200.5, 210.5, 220.5, 230.5, 240.5,
 168                           250.5, 275.5, 300.5, 325.5, 350.5, 375.5, 400.5, 425.5, 450.5, 475.5,
 169                           500.5 };
 170                                 //525.5, 550.5, 575.5, 600.5, 625.5, 650.5, 675.5, 700.5, 725.5,
 171                           //750.5, 775.5, 800.5, 825.5, 850.5, 875.5, 900.5, 925.5, 950.5, 975.5,
 172                           //1000.5 };
 173
 174   #define VTXBINNING 10, binLimitsVtx
 175   #define NBINNING fgkMaxParams, binLimitsN*/
 176
 177   #define NBINNING 201, -0.5, 200.5
 178   #define VTXBINNING 1, -6, 6
 179
 180   for (Int_t i = 0; i < kESDHists; ++i)
 181     fMultiplicityESD[i] = new TH2F(Form("fMultiplicityESD%d", i), "fMultiplicityESD;vtx-z;Ntracks;Count", VTXBINNING, NBINNING);
 182
 183   for (Int_t i = 0; i < kMCHists; ++i)
 184   {
 185     fMultiplicityVtx[i] = dynamic_cast<TH2F*> (fMultiplicityESD[0]->Clone(Form("fMultiplicityVtx%d", i)));
 186     fMultiplicityVtx[i]->SetTitle("fMultiplicityVtx;vtx-z;Npart");
 187
 188     fMultiplicityMB[i] = dynamic_cast<TH2F*> (fMultiplicityVtx[0]->Clone(Form("fMultiplicityMB%d", i)));
 189     fMultiplicityMB[i]->SetTitle("fMultiplicityMB");
 190
 191     fMultiplicityINEL[i] = dynamic_cast<TH2F*> (fMultiplicityVtx[0]->Clone(Form("fMultiplicityINEL%d", i)));
 192     fMultiplicityINEL[i]->SetTitle("fMultiplicityINEL");
 193
 194     fMultiplicityNSD[i] = dynamic_cast<TH2F*> (fMultiplicityVtx[0]->Clone(Form("fMultiplicityNSD%d", i)));
 195     fMultiplicityNSD[i]->SetTitle("fMultiplicityNSD");
 196   }
 197
 198   for (Int_t i = 0; i < kCorrHists; ++i)
 199   {
 200     fCorrelation[i] = new TH3F(Form("fCorrelation%d", i), "fCorrelation;vtx-z;Npart;Ntracks", VTXBINNING, NBINNING, NBINNING);
 201     fMultiplicityESDCorrected[i] = new TH1F(Form("fMultiplicityESDCorrected%d", i), "fMultiplicityESDCorrected;Npart;dN/dN", NBINNING);
 202   }
 203
 204   TH1::AddDirectory(oldStatus);
 205
 206   AliUnfolding::SetNbins(120, 120);
 207   AliUnfolding::SetSkipBinsBegin(1);
 208   AliUnfolding::SetNormalizeInput(kTRUE);
 209 }
 210
 211 //____________________________________________________________________
 212 AliMultiplicityCorrection::~AliMultiplicityCorrection()
 213 {
 214   //
 215   // Destructor
 216   //
 217
 218   Printf("AliMultiplicityCorrection::~AliMultiplicityCorrection called");
 219
 220   for (Int_t i = 0; i < kESDHists; ++i)
 221   {
 222     if (fMultiplicityESD[i])
 223       delete fMultiplicityESD[i];
 224     fMultiplicityESD[i] = 0;
 225   }
 226
 227   for (Int_t i = 0; i < kMCHists; ++i)
 228   {
 229     if (fMultiplicityVtx[i])
 230       delete fMultiplicityVtx[i];
 231     fMultiplicityVtx[i] = 0;
 232
 233     if (fMultiplicityMB[i])
 234       delete fMultiplicityMB[i];
 235     fMultiplicityMB[i] = 0;
 236
 237     if (fMultiplicityINEL[i])
 238       delete fMultiplicityINEL[i];
 239     fMultiplicityINEL[i] = 0;
 240
 241     if (fMultiplicityNSD[i])
 242       delete fMultiplicityNSD[i];
 243     fMultiplicityNSD[i] = 0;
 244 }
 245
 246   for (Int_t i = 0; i < kCorrHists; ++i)
 247   {
 248     if (fCorrelation[i])
 249       delete fCorrelation[i];
 250     fCorrelation[i] = 0;
 251
 252     if (fMultiplicityESDCorrected[i])
 253       delete fMultiplicityESDCorrected[i];
 254     fMultiplicityESDCorrected[i] = 0;
 255   }
 256 }
 257
 258 //____________________________________________________________________
 259 Long64_t AliMultiplicityCorrection::Merge(TCollection* list)
 260 {
 261   // Merge a list of AliMultiplicityCorrection objects with this (needed for
 262   // PROOF).
 263   // Returns the number of merged objects (including this).
 264
 265   if (!list)
 266     return 0;
 267
 268   if (list->IsEmpty())
 269     return 1;
 270
 271   TIterator* iter = list->MakeIterator();
 272   TObject* obj;
 273
 274   // collections of all histograms
 275   TList collections[kESDHists+kMCHists*4+kCorrHists*2];
 276
 277   Int_t count = 0;
 278   while ((obj = iter->Next())) {
 279
 280     AliMultiplicityCorrection* entry = dynamic_cast<AliMultiplicityCorrection*> (obj);
 281     if (entry == 0)
 282       continue;
 283
 284     for (Int_t i = 0; i < kESDHists; ++i)
 285       collections[i].Add(entry->fMultiplicityESD[i]);
 286
 287     for (Int_t i = 0; i < kMCHists; ++i)
 288     {
 289       collections[kESDHists+i].Add(entry->fMultiplicityVtx[i]);
 290       collections[kESDHists+kMCHists+i].Add(entry->fMultiplicityMB[i]);
 291       collections[kESDHists+kMCHists*2+i].Add(entry->fMultiplicityINEL[i]);
 292       collections[kESDHists+kMCHists*3+i].Add(entry->fMultiplicityNSD[i]);
 293     }
 294
 295     for (Int_t i = 0; i < kCorrHists; ++i)
 296       collections[kESDHists+kMCHists*4+i].Add(entry->fCorrelation[i]);
 297
 298     for (Int_t i = 0; i < kCorrHists; ++i)
 299       collections[kESDHists+kMCHists*4+kCorrHists+i].Add(entry->fMultiplicityESDCorrected[i]);
 300
 301     count++;
 302   }
 303
 304   for (Int_t i = 0; i < kESDHists; ++i)
 305     fMultiplicityESD[i]->Merge(&collections[i]);
 306
 307   for (Int_t i = 0; i < kMCHists; ++i)
 308   {
 309     fMultiplicityVtx[i]->Merge(&collections[kESDHists+i]);
 310     fMultiplicityMB[i]->Merge(&collections[kESDHists+kMCHists+i]);
 311     fMultiplicityINEL[i]->Merge(&collections[kESDHists+kMCHists*2+i]);
 312     fMultiplicityNSD[i]->Merge(&collections[kESDHists+kMCHists*3+i]);
 313   }
 314
 315   for (Int_t i = 0; i < kCorrHists; ++i)
 316     fCorrelation[i]->Merge(&collections[kESDHists+kMCHists*4+i]);
 317
 318   for (Int_t i = 0; i < kCorrHists; ++i)
 319     fMultiplicityESDCorrected[i]->Merge(&collections[kESDHists+kMCHists*4+kCorrHists+i]);
 320
 321   delete iter;
 322
 323   return count+1;
 324 }
 325
 326 //____________________________________________________________________
 327 Bool_t AliMultiplicityCorrection::LoadHistograms(const Char_t* dir)
 328 {
 329   //
 330   // loads the histograms from a file
 331   // if dir is empty a directory with the name of this object is taken (like in SaveHistogram)
 332   //
 333
 334   if (!dir)
 335     dir = GetName();
 336
 337   if (!gDirectory->cd(dir))
 338     return kFALSE;
 339
 340   // store old hists to delete them later
 341   TList oldObjects;
 342   oldObjects.SetOwner(1);
 343   for (Int_t i = 0; i < kESDHists; ++i)
 344     if (fMultiplicityESD[i])
 345       oldObjects.Add(fMultiplicityESD[i]);
 346
 347   for (Int_t i = 0; i < kMCHists; ++i)
 348   {
 349     if (fMultiplicityVtx[i])
 350       oldObjects.Add(fMultiplicityVtx[i]);
 351     if (fMultiplicityMB[i])
 352       oldObjects.Add(fMultiplicityMB[i]);
 353     if (fMultiplicityINEL[i])
 354       oldObjects.Add(fMultiplicityINEL[i]);
 355     if (fMultiplicityNSD[i])
 356       oldObjects.Add(fMultiplicityNSD[i]);
 357   }
 358
 359   for (Int_t i = 0; i < kCorrHists; ++i)
 360     if (fCorrelation[i])
 361       oldObjects.Add(fCorrelation[i]);
 362
 363   // load histograms
 364
 365   Bool_t success = kTRUE;
 366
 367   for (Int_t i = 0; i < kESDHists; ++i)
 368   {
 369     fMultiplicityESD[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityESD[i]->GetName()));
 370     if (!fMultiplicityESD[i])
 371       success = kFALSE;
 372   }
 373
 374   for (Int_t i = 0; i < kMCHists; ++i)
 375   {
 376     fMultiplicityVtx[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityVtx[i]->GetName()));
 377     fMultiplicityMB[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityMB[i]->GetName()));
 378     fMultiplicityINEL[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityINEL[i]->GetName()));
 379     fMultiplicityNSD[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityNSD[i]->GetName()));
 380     if (!fMultiplicityVtx[i] || !fMultiplicityMB[i] || !fMultiplicityINEL[i])
 381       success = kFALSE;
 382   }
 383
 384   for (Int_t i = 0; i < kCorrHists; ++i)
 385   {
 386     fCorrelation[i] = dynamic_cast<TH3F*> (gDirectory->Get(fCorrelation[i]->GetName()));
 387     if (!fCorrelation[i])
 388       success = kFALSE;
 389     fMultiplicityESDCorrected[i] = dynamic_cast<TH1F*> (gDirectory->Get(fMultiplicityESDCorrected[i]->GetName()));
 390     if (!fMultiplicityESDCorrected[i])
 391       success = kFALSE;
 392   }
 393
 394   gDirectory->cd("..");
 395
 396   // delete old hists
 397   oldObjects.Delete();
 398
 399   return success;
 400 }
 401
 402 //____________________________________________________________________
 403 void AliMultiplicityCorrection::SaveHistograms(const char* dir)
 404 {
 405   //
 406   // saves the histograms
 407   //
 408
 409   if (!dir)
 410     dir = GetName();
 411
 412   gDirectory->mkdir(dir);
 413   gDirectory->cd(dir);
 414
 415   for (Int_t i = 0; i < kESDHists; ++i)
 416     if (fMultiplicityESD[i])
 417     {
 418       fMultiplicityESD[i]->Write();
 419       fMultiplicityESD[i]->ProjectionY(Form("%s_px", fMultiplicityESD[i]->GetName()), 1, fMultiplicityESD[i]->GetNbinsX())->Write();
 420     }
 421
 422   for (Int_t i = 0; i < kMCHists; ++i)
 423   {
 424     if (fMultiplicityVtx[i])
 425     {
 426       fMultiplicityVtx[i]->Write();
 427       fMultiplicityVtx[i]->ProjectionY(Form("%s_px", fMultiplicityVtx[i]->GetName()), 1, fMultiplicityVtx[i]->GetNbinsX())->Write();
 428     }
 429     if (fMultiplicityMB[i])
 430     {
 431       fMultiplicityMB[i]->Write();
 432       fMultiplicityMB[i]->ProjectionY(Form("%s_px", fMultiplicityMB[i]->GetName()), 1, fMultiplicityMB[i]->GetNbinsX())->Write();
 433     }
 434     if (fMultiplicityINEL[i])
 435     {
 436       fMultiplicityINEL[i]->Write();
 437       fMultiplicityINEL[i]->ProjectionY(Form("%s_px", fMultiplicityINEL[i]->GetName()), 1, fMultiplicityINEL[i]->GetNbinsX())->Write();
 438     }
 439     if (fMultiplicityNSD[i])
 440     {
 441       fMultiplicityNSD[i]->Write();
 442       fMultiplicityNSD[i]->ProjectionY(Form("%s_px", fMultiplicityNSD[i]->GetName()), 1, fMultiplicityNSD[i]->GetNbinsX())->Write();
 443     }
 444   }
 445
 446   for (Int_t i = 0; i < kCorrHists; ++i)
 447   {
 448     if (fCorrelation[i])
 449       fCorrelation[i]->Write();
 450     if (fMultiplicityESDCorrected[i])
 451       fMultiplicityESDCorrected[i]->Write();
 452   }
 453
 454   gDirectory->cd("..");
 455 }
 456
 457 //____________________________________________________________________
 458 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)
 459 {
 460   //
 461   // Fills an event from MC
 462   //
 463
 464   if (triggered)
 465   {
 466     fMultiplicityMB[0]->Fill(vtx, generated05);
 467     fMultiplicityMB[1]->Fill(vtx, generated10);
 468     fMultiplicityMB[2]->Fill(vtx, generated15);
 469     fMultiplicityMB[3]->Fill(vtx, generated20);
 470     fMultiplicityMB[4]->Fill(vtx, generatedAll);
 471
 472     if (vertex)
 473     {
 474       fMultiplicityVtx[0]->Fill(vtx, generated05);
 475       fMultiplicityVtx[1]->Fill(vtx, generated10);
 476       fMultiplicityVtx[2]->Fill(vtx, generated15);
 477       fMultiplicityVtx[3]->Fill(vtx, generated20);
 478       fMultiplicityVtx[4]->Fill(vtx, generatedAll);
 479     }
 480   }
 481
 482   fMultiplicityINEL[0]->Fill(vtx, generated05);
 483   fMultiplicityINEL[1]->Fill(vtx, generated10);
 484   fMultiplicityINEL[2]->Fill(vtx, generated15);
 485   fMultiplicityINEL[3]->Fill(vtx, generated20);
 486   fMultiplicityINEL[4]->Fill(vtx, generatedAll);
 487
 488   if (processType != AliPWG0Helper::kSD)
 489   {
 490     fMultiplicityNSD[0]->Fill(vtx, generated05);
 491     fMultiplicityNSD[1]->Fill(vtx, generated10);
 492     fMultiplicityNSD[2]->Fill(vtx, generated15);
 493     fMultiplicityNSD[3]->Fill(vtx, generated20);
 494     fMultiplicityNSD[4]->Fill(vtx, generatedAll);
 495   }
 496 }
 497
 498 //____________________________________________________________________
 499 void AliMultiplicityCorrection::FillMeasured(Float_t vtx, Int_t measured05, Int_t measured10, Int_t measured15, Int_t measured20)
 500 {
 501   //
 502   // Fills an event from ESD
 503   //
 504
 505   fMultiplicityESD[0]->Fill(vtx, measured05);
 506   fMultiplicityESD[1]->Fill(vtx, measured10);
 507   fMultiplicityESD[2]->Fill(vtx, measured15);
 508   fMultiplicityESD[3]->Fill(vtx, measured20);
 509 }
 510
 511 //____________________________________________________________________
 512 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)
 513 {
 514   //
 515   // Fills an event into the correlation map with the information from MC and ESD
 516   //
 517
 518   fCorrelation[0]->Fill(vtx, generated05, measured05);
 519   fCorrelation[1]->Fill(vtx, generated10, measured10);
 520   fCorrelation[2]->Fill(vtx, generated15, measured15);
 521   fCorrelation[3]->Fill(vtx, generated20, measured20);
 522
 523   fCorrelation[4]->Fill(vtx, generatedAll, measured05);
 524   fCorrelation[5]->Fill(vtx, generatedAll, measured10);
 525   fCorrelation[6]->Fill(vtx, generatedAll, measured15);
 526   fCorrelation[7]->Fill(vtx, generatedAll, measured20);
 527 }
 528
 529 //____________________________________________________________________
 530 void AliMultiplicityCorrection::SetupCurrentHists(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType)
 531 {
 532   //
 533   // fills fCurrentESD, fCurrentCorrelation
 534   // resets fMultiplicityESDCorrected
 535   //
 536
 537   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
 538
 539   fMultiplicityESDCorrected[correlationID]->Reset();
 540   fMultiplicityESDCorrected[correlationID]->Sumw2();
 541
 542   // project without under/overflow bins
 543   fCurrentESD = fMultiplicityESD[inputRange]->ProjectionY("fCurrentESD", 1, fMultiplicityESD[inputRange]->GetXaxis()->GetNbins());
 544   fCurrentESD->Sumw2();
 545
 546   // empty under/overflow bins in x, otherwise Project3D takes them into account
 547   TH3* hist = fCorrelation[correlationID];
 548   for (Int_t y=0; y<=hist->GetYaxis()->GetNbins()+1; ++y)
 549   {
 550     for (Int_t z=0; z<=hist->GetZaxis()->GetNbins()+1; ++z)
 551     {
 552       hist->SetBinContent(0, y, z, 0);
 553       hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0);
 554     }
 555   }
 556
 557   fCurrentCorrelation = (TH2*) hist->Project3D("zy");
 558   fCurrentCorrelation->Sumw2();
 559
 560   Printf("AliMultiplicityCorrection::SetupCurrentHists: Statistics information: %.f entries in correlation map; %.f entries in measured spectrum", fCurrentCorrelation->Integral(), fCurrentESD->Integral());
 561
 562 #if 0 // does not help
 563   // null bins with one entry
 564   Int_t nNulledBins = 0;
 565   for (Int_t x=1; x<=fCurrentCorrelation->GetXaxis()->GetNbins(); ++x)
 566     for (Int_t y=1; y<=fCurrentCorrelation->GetYaxis()->GetNbins(); ++y)
 567     {
 568       if (fCurrentCorrelation->GetBinContent(x, y) == 1)
 569       {
 570         fCurrentCorrelation->SetBinContent(x, y, 0);
 571         fCurrentCorrelation->SetBinError(x, y, 0);
 572
 573         ++nNulledBins;
 574       }
 575     }
 576   Printf("Nulled %d bins", nNulledBins);
 577 #endif
 578
 579   fCurrentEfficiency = GetEfficiency(inputRange, eventType);
 580   //fCurrentEfficiency->Rebin(2);
 581   //fCurrentEfficiency->Scale(0.5);
 582 }
 583
 584 //____________________________________________________________________
 585 TH1* AliMultiplicityCorrection::GetEfficiency(Int_t inputRange, EventType eventType)
 586 {
 587   //
 588   // calculates efficiency for given event type
 589   //
 590
 591   TH1* divisor = 0;
 592   switch (eventType)
 593   {
 594     case kTrVtx : break;
 595     case kMB: divisor = fMultiplicityMB[inputRange]->ProjectionY("divisor", 1, fMultiplicityMB[inputRange]->GetNbinsX(), "e"); break;
 596     case kINEL: divisor = fMultiplicityINEL[inputRange]->ProjectionY("divisor", 1, fMultiplicityINEL[inputRange]->GetNbinsX(), "e"); break;
 597     case kNSD: divisor = fMultiplicityNSD[inputRange]->ProjectionY("divisor", 1, fMultiplicityNSD[inputRange]->GetNbinsX(), "e"); break;
 598   }
 599   TH1* eff = fMultiplicityVtx[inputRange]->ProjectionY("CurrentEfficiency", 1, fMultiplicityVtx[inputRange]->GetNbinsX(), "e");
 600
 601   if (eventType == kTrVtx)
 602   {
 603     for (Int_t i=0; i<= eff->GetNbinsX()+1; i++)
 604       eff->SetBinContent(i, 1);
 605   }
 606   else
 607     eff->Divide(eff, divisor, 1, 1, "B");
 608
 609   return eff;
 610 }
 611
 612 //____________________________________________________________________
 613 TH1* AliMultiplicityCorrection::GetTriggerEfficiency(Int_t inputRange)
 614 {
 615   //
 616   // calculates efficiency for given event type
 617   //
 618
 619   TH1* divisor = fMultiplicityINEL[inputRange]->ProjectionY("divisor", 1, fMultiplicityINEL[inputRange]->GetNbinsX(), "e");
 620   TH1* eff = fMultiplicityMB[inputRange]->ProjectionY("CurrentEfficiency", 1, fMultiplicityMB[inputRange]->GetNbinsX(), "e");
 621   eff->Divide(eff, divisor, 1, 1, "B");
 622   return eff;
 623 }
 624
 625 //____________________________________________________________________
 626 Int_t AliMultiplicityCorrection::ApplyMinuitFit(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t check, TH1* initialConditions)
 627 {
 628   //
 629   // correct spectrum using minuit chi2 method
 630   //
 631   // for description of parameters, see AliUnfolding::Unfold
 632   //
 633
 634   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
 635
 636   AliUnfolding::SetCreateOverflowBin(5);
 637   AliUnfolding::SetUnfoldingMethod(AliUnfolding::kChi2Minimization);
 638   SetupCurrentHists(inputRange, fullPhaseSpace, eventType);
 639
 640   if (!initialConditions)
 641   {
 642     initialConditions = (TH1*) fCurrentESD->Clone("initialConditions");
 643     initialConditions->Scale(1.0 / initialConditions->Integral());
 644     if (!check)
 645     {
 646       // set minimum value to prevent MINUIT just staying in the small value
 647       for (Int_t i=1; i<=initialConditions->GetNbinsX(); i++)
 648         initialConditions->SetBinContent(i, TMath::Max(initialConditions->GetBinContent(i), 1e-3));
 649     }
 650   }
 651
 652   return AliUnfolding::Unfold(fCurrentCorrelation, fCurrentEfficiency, fCurrentESD, initialConditions, fMultiplicityESDCorrected[correlationID], check);
 653 }
 654
 655 //____________________________________________________________________
 656 Int_t AliMultiplicityCorrection::ApplyNBDFit(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType)
 657 {
 658   //
 659   // correct spectrum using minuit chi2 method with a NBD function
 660   //
 661   // for description of parameters, see AliUnfolding::Unfold
 662   //
 663
 664   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
 665
 666   AliUnfolding::SetUnfoldingMethod(AliUnfolding::kFunction);
 667   SetupCurrentHists(inputRange, fullPhaseSpace, eventType);
 668
 669         TF1* func = new TF1("nbd", "[0] * TMath::Gamma([2]+x) / TMath::Gamma([2]) / TMath::Gamma(x+1) * pow([1] / ([1]+[2]), x) * pow(1.0 + [1]/[2], -[2])");
 670         func->SetParNames("scaling", "averagen", "k");
 671         func->SetParLimits(0, 0, 1000);
 672         func->SetParLimits(1, 1, 50);
 673         func->SetParLimits(2, 1, 10);
 674         func->SetParameters(1, 10, 2);
 675   AliUnfolding::SetFunction(func);
 676
 677   return AliUnfolding::Unfold(fCurrentCorrelation, fCurrentEfficiency, fCurrentESD, 0, fMultiplicityESDCorrected[correlationID]);
 678 }
 679
 680 //____________________________________________________________________
 681 void AliMultiplicityCorrection::DrawHistograms()
 682 {
 683   //
 684   // draws the histograms of this class
 685   //
 686
 687   printf("ESD:\n");
 688
 689   TCanvas* canvas1 = new TCanvas("fMultiplicityESD", "fMultiplicityESD", 900, 600);
 690   canvas1->Divide(3, 2);
 691   for (Int_t i = 0; i < kESDHists; ++i)
 692   {
 693     canvas1->cd(i+1);
 694     fMultiplicityESD[i]->DrawCopy("COLZ");
 695     printf("%d --> %f\n", i, (Float_t) fMultiplicityESD[i]->ProjectionY()->GetMean());
 696   }
 697
 698   printf("Vtx:\n");
 699
 700   TCanvas* canvas2 = new TCanvas("fMultiplicityMC", "fMultiplicityMC", 900, 600);
 701   canvas2->Divide(3, 2);
 702   for (Int_t i = 0; i < kMCHists; ++i)
 703   {
 704     canvas2->cd(i+1);
 705     fMultiplicityVtx[i]->DrawCopy("COLZ");
 706     printf("%d --> %f\n", i, fMultiplicityVtx[i]->ProjectionY()->GetMean());
 707   }
 708
 709   TCanvas* canvas3 = new TCanvas("fCorrelation", "fCorrelation", 900, 900);
 710   canvas3->Divide(3, 3);
 711   for (Int_t i = 0; i < kCorrHists; ++i)
 712   {
 713     canvas3->cd(i+1);
 714     TH3* hist = dynamic_cast<TH3*> (fCorrelation[i]->Clone());
 715     for (Int_t y=1; y<=hist->GetYaxis()->GetNbins(); ++y)
 716     {
 717       for (Int_t z=1; z<=hist->GetZaxis()->GetNbins(); ++z)
 718       {
 719         hist->SetBinContent(0, y, z, 0);
 720         hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0);
 721       }
 722     }
 723     TH1* proj = hist->Project3D("zy");
 724     proj->DrawCopy("COLZ");
 725   }
 726 }
 727
 728 //____________________________________________________________________
 729 void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRange, Bool_t fullPhaseSpace, Bool_t /*normalizeESD*/, TH1* mcHist, Bool_t simple)
 730 {
 731   // draw comparison plots
 732
 733
 734   //mcHist->Rebin(2);
 735
 736   Int_t esdCorrId = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
 737
 738   TString tmpStr;
 739   tmpStr.Form("%s_DrawComparison_%d", name, esdCorrId);
 740
 741   if (fMultiplicityESDCorrected[esdCorrId]->Integral() == 0)
 742   {
 743     printf("ERROR. Unfolded histogram is empty\n");
 744     return;
 745   }
 746
 747   //regain measured distribution used for unfolding, because the bins at high mult. were modified in SetupCurrentHists
 748   fCurrentESD = fMultiplicityESD[esdCorrId]->ProjectionY("fCurrentESD", 1, fMultiplicityESD[inputRange]->GetXaxis()->GetNbins());
 749   fCurrentESD->Sumw2();
 750   fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
 751
 752   // normalize unfolded result to 1
 753   fMultiplicityESDCorrected[esdCorrId]->Scale(1.0 / fMultiplicityESDCorrected[esdCorrId]->Integral());
 754
 755   //fCurrentESD->Scale(mcHist->Integral(2, 200));
 756
 757   //new TCanvas;
 758   /*TH1* ratio = (TH1*) fMultiplicityESDCorrected[esdCorrId]->Clone("ratio");
 759   ratio->Divide(mcHist);
 760   ratio->Draw("HIST");
 761   ratio->Fit("pol0", "W0", "", 20, 120);
 762   Float_t scalingFactor = ratio->GetFunction("pol0")->GetParameter(0);
 763   delete ratio;
 764   mcHist->Scale(scalingFactor);*/
 765
 766   // find bin with <= 50 entries. this is used as limit for the chi2 calculation
 767   Int_t mcBinLimit = 0;
 768   for (Int_t i=20; i<mcHist->GetNbinsX(); ++i)
 769   {
 770     if (mcHist->GetBinContent(i) > 50)
 771     {
 772       mcBinLimit = i;
 773     }
 774     else
 775       break;
 776   }
 777   Printf("AliMultiplicityCorrection::DrawComparison: MC bin limit is %d", mcBinLimit);
 778
 779   // scale to 1
 780   mcHist->Sumw2();
 781   if (mcHist->Integral() > 0)
 782     mcHist->Scale(1.0 / mcHist->Integral());
 783
 784   // calculate residual
 785
 786   // for that we convolute the response matrix with the gathered result
 787   // if normalizeESD is not set, the histogram is already normalized, this needs to be passed to CalculateMultiplicityESD
 788   TH1* tmpESDEfficiencyRecorrected = (TH1*) fMultiplicityESDCorrected[esdCorrId]->Clone("tmpESDEfficiencyRecorrected");
 789   tmpESDEfficiencyRecorrected->Multiply(fCurrentEfficiency);
 790   TH2* convoluted = CalculateMultiplicityESD(tmpESDEfficiencyRecorrected, esdCorrId);
 791   TH1* convolutedProj = convoluted->ProjectionY("convolutedProj", -1, -1, "e");
 792   if (convolutedProj->Integral() > 0)
 793   {
 794     convolutedProj->Scale(1.0 / convolutedProj->Integral());
 795   }
 796   else
 797     printf("ERROR: convolutedProj is empty. Something went wrong calculating the convoluted histogram.\n");
 798
 799   TH1* residual = (TH1*) convolutedProj->Clone("residual");
 800   residual->SetTitle("Residuals;Ntracks;(folded unfolded measured - measured) / e");
 801
 802   residual->Add(fCurrentESD, -1);
 803   //residual->Divide(residual, fCurrentESD, 1, 1, "B");
 804
 805   TH1* residualHist = new TH1F("residualHist", "residualHist", 51, -5, 5);
 806
 807   // find bin limit
 808   Int_t lastBin = 0;
 809   for (Int_t i=1; i<=fCurrentESD->GetNbinsX(); ++i)
 810   {
 811     if (fCurrentESD->GetBinContent(i) <= 5)
 812     {
 813       lastBin = i;
 814       break;
 815     }
 816   }
 817
 818   // TODO fix errors
 819   Float_t chi2 = 0;
 820   for (Int_t i=1; i<=residual->GetNbinsX(); ++i)
 821   {
 822     if (fCurrentESD->GetBinError(i) > 0)
 823     {
 824       Float_t value = residual->GetBinContent(i) / fCurrentESD->GetBinError(i);
 825       if (i > 1 && i <= lastBin)
 826         chi2 += value * value;
 827       Printf("%d --> %f (%f)", i, value * value, chi2);
 828       residual->SetBinContent(i, value);
 829       residualHist->Fill(value);
 830     }
 831     else
 832     {
 833       //printf("Residual bin %d set to 0\n", i);
 834       residual->SetBinContent(i, 0);
 835     }
 836     convolutedProj->SetBinError(i, 0);
 837     residual->SetBinError(i, 0);
 838   }
 839   fLastChi2Residuals = chi2;
 840
 841   new TCanvas; residualHist->DrawCopy();
 842
 843   //residualHist->Fit("gaus", "N");
 844   //delete residualHist;
 845
 846   printf("Difference (Residuals) is %f for bin 2-%d\n", fLastChi2Residuals, lastBin);
 847
 848   TCanvas* canvas1 = 0;
 849   if (simple)
 850   {
 851     canvas1 = new TCanvas(tmpStr, tmpStr, 1200, 600);
 852     canvas1->Divide(2, 1);
 853   }
 854   else
 855   {
 856     canvas1 = new TCanvas(tmpStr, tmpStr, 1200, 1200);
 857     canvas1->Divide(2, 3);
 858   }
 859
 860   canvas1->cd(1);
 861   canvas1->cd(1)->SetGridx();
 862   canvas1->cd(1)->SetGridy();
 863   canvas1->cd(1)->SetTopMargin(0.05);
 864   canvas1->cd(1)->SetRightMargin(0.05);
 865   canvas1->cd(1)->SetLeftMargin(0.12);
 866   canvas1->cd(1)->SetBottomMargin(0.12);
 867   TH1* proj = (TH1*) mcHist->Clone("proj");
 868
 869   // normalize without 0 bin
 870   proj->Scale(1.0 / proj->Integral(2, proj->GetNbinsX()));
 871   Printf("Normalized without 0 bin!");
 872   proj->GetXaxis()->SetRangeUser(0, 200);
 873   proj->GetYaxis()->SetTitleOffset(1.4);
 874   //proj->SetLabelSize(0.05, "xy");
 875   //proj->SetTitleSize(0.05, "xy");
 876   proj->SetTitle(Form(";True multiplicity in |#eta| < %.1f;Entries", (inputRange+1)*0.5));
 877   proj->SetStats(kFALSE);
 878
 879   fMultiplicityESDCorrected[esdCorrId]->SetLineColor(2);
 880   fMultiplicityESDCorrected[esdCorrId]->SetMarkerColor(2);
 881   //fMultiplicityESDCorrected[esdCorrId]->SetMarkerStyle(5);
 882   // normalize without 0 bin
 883   fMultiplicityESDCorrected[esdCorrId]->Scale(1.0 / fMultiplicityESDCorrected[esdCorrId]->Integral(2, fMultiplicityESDCorrected[esdCorrId]->GetNbinsX()));
 884   Printf("Normalized without 0 bin!");
 885
 886   if (proj->GetEntries() > 0) {
 887     proj->DrawCopy("HIST");
 888     fMultiplicityESDCorrected[esdCorrId]->DrawCopy("SAME HIST E");
 889   }
 890   else
 891     fMultiplicityESDCorrected[esdCorrId]->DrawCopy("HIST E");
 892
 893   gPad->SetLogy();
 894
 895   TLegend* legend = new TLegend(0.3, 0.8, 0.93, 0.93);
 896   legend->AddEntry(proj, "True distribution");
 897   legend->AddEntry(fMultiplicityESDCorrected[esdCorrId], "Unfolded distribution");
 898   legend->SetFillColor(0);
 899   legend->SetTextSize(0.04);
 900   legend->Draw();
 901   // unfortunately does not work. maybe a bug? --> legend->SetTextSizePixels(14);
 902
 903   canvas1->cd(2);
 904   canvas1->cd(2)->SetGridx();
 905   canvas1->cd(2)->SetGridy();
 906   canvas1->cd(2)->SetTopMargin(0.05);
 907   canvas1->cd(2)->SetRightMargin(0.05);
 908   canvas1->cd(2)->SetLeftMargin(0.12);
 909   canvas1->cd(2)->SetBottomMargin(0.12);
 910
 911   gPad->SetLogy();
 912   fCurrentESD->GetXaxis()->SetRangeUser(0, 200);
 913   //fCurrentESD->SetLineColor(2);
 914   fCurrentESD->SetTitle(Form(";Measured multiplicity in |#eta| < %.1f;Entries", (inputRange+1)*0.5));
 915   fCurrentESD->SetStats(kFALSE);
 916   fCurrentESD->GetYaxis()->SetTitleOffset(1.4);
 917   //fCurrentESD->SetLabelSize(0.05, "xy");
 918   //fCurrentESD->SetTitleSize(0.05, "xy");
 919   fCurrentESD->DrawCopy("HIST E");
 920
 921   convolutedProj->SetLineColor(2);
 922   convolutedProj->SetMarkerColor(2);
 923   convolutedProj->SetMarkerStyle(5);
 924   //proj2->SetMarkerColor(2);
 925   //proj2->SetMarkerStyle(5);
 926   convolutedProj->DrawCopy("HIST SAME P");
 927
 928   legend = new TLegend(0.3, 0.8, 0.93, 0.93);
 929   legend->AddEntry(fCurrentESD, "Measured distribution");
 930   legend->AddEntry(convolutedProj, "R #otimes unfolded distribution", "P");
 931   legend->SetFillColor(0);
 932   legend->SetTextSize(0.04);
 933   legend->Draw();
 934
 935   //TH1* diffMCUnfolded = dynamic_cast<TH1*> (proj->Clone("diffMCUnfolded"));
 936   //diffMCUnfolded->Add(fMultiplicityESDCorrected[esdCorrId], -1);
 937
 938   /*Float_t chi2 = 0;
 939   Float_t chi = 0;
 940   fLastChi2MCLimit = 0;
 941   Int_t limit = 0;
 942   for (Int_t i=2; i<=200; ++i)
 943   {
 944     if (proj->GetBinContent(i) != 0)
 945     {
 946       Float_t value = (proj->GetBinContent(i) - fMultiplicityESDCorrected[esdCorrId]->GetBinContent(i)) / proj->GetBinContent(i);
 947       chi2 += value * value;
 948       chi += TMath::Abs(value);
 949
 950       //printf("%d %f\n", i, chi);
 951
 952       if (chi2 < 0.2)
 953         fLastChi2MCLimit = i;
 954
 955       if (chi < 3)
 956         limit = i;
 957
 958     }
 959   }*/
 960
 961   /*chi2 = 0;
 962   Float_t chi = 0;
 963   Int_t limit = 0;
 964   for (Int_t i=1; i<=diffMCUnfolded->GetNbinsX(); ++i)
 965   {
 966     if (fMultiplicityESDCorrected[esdCorrId]->GetBinError(i) > 0)
 967     {
 968       Double_t value = diffMCUnfolded->GetBinContent(i) / fMultiplicityESDCorrected[esdCorrId]->GetBinError(i);
 969       if (value > 1e8)
 970         value = 1e8; //prevent arithmetic exception
 971       else if (value < -1e8)
 972         value = -1e8;
 973       if (i > 1)
 974       {
 975         chi2 += value * value;
 976         chi += TMath::Abs(value);
 977       }
 978       diffMCUnfolded->SetBinContent(i, value);
 979     }
 980     else
 981     {
 982       //printf("diffMCUnfolded bin %d set to 0\n", i);
 983       diffMCUnfolded->SetBinContent(i, 0);
 984     }
 985     if (chi2 < 1000)
 986       fLastChi2MCLimit = i;
 987     if (chi < 1000)
 988       limit = i;
 989     if (i == 150)
 990       fLastChi2MC = chi2;
 991   }
 992
 993   printf("limits %d %d\n", fLastChi2MCLimit, limit);
 994   fLastChi2MCLimit = limit;
 995
 996   printf("Difference (from MC) is %f for bin 2-150. Limit is %d.\n", fLastChi2MC, fLastChi2MCLimit);*/
 997
 998   if (!simple)
 999   {
1000     /*canvas1->cd(3);
1001
1002     diffMCUnfolded->SetTitle("#chi^{2};Npart;(MC - Unfolded) / e(unfolded)");
1003     //diffMCUnfolded->GetYaxis()->SetRangeUser(-20, 20);
1004     diffMCUnfolded->GetXaxis()->SetRangeUser(0, 200);
1005     diffMCUnfolded->DrawCopy("HIST");
1006
1007     TH1F* fluctuation = new TH1F("fluctuation", "fluctuation", 20, -5, 5);
1008     for (Int_t i=20; i<=diffMCUnfolded->GetNbinsX(); ++i)
1009       fluctuation->Fill(diffMCUnfolded->GetBinContent(i));
1010
1011     //new TCanvas; fluctuation->DrawCopy();
1012     delete fluctuation;*/
1013
1014     /*TLegend* legend = new TLegend(0.6, 0.7, 0.85, 0.85);
1015     legend->AddEntry(fMultiplicityESDCorrected, "ESD corrected");
1016     legend->AddEntry(fMultiplicityMC, "MC");
1017     legend->AddEntry(fMultiplicityESD, "ESD");
1018     legend->Draw();*/
1019
1020     canvas1->cd(4);
1021     residual->GetYaxis()->SetRangeUser(-5, 5);
1022     residual->GetXaxis()->SetRangeUser(0, 200);
1023     residual->DrawCopy();
1024
1025     canvas1->cd(5);
1026
1027     TH1* ratio = (TH1*) fMultiplicityESDCorrected[esdCorrId]->Clone("ratio");
1028     ratio->Divide(mcHist);
1029     ratio->SetTitle("Ratio;true multiplicity;Unfolded / MC");
1030     ratio->GetYaxis()->SetRangeUser(0.5, 1.5);
1031     ratio->GetXaxis()->SetRangeUser(0, 200);
1032     ratio->SetStats(kFALSE);
1033     ratio->Draw("HIST");
1034
1035     Double_t ratioChi2 = 0;
1036     fRatioAverage = 0;
1037     fLastChi2MCLimit = 0;
1038     for (Int_t i=2; i<=150; ++i)
1039     {
1040       Float_t value = ratio->GetBinContent(i) - 1;
1041       if (value > 1e8)
1042         value = 1e8; //prevent arithmetic exception
1043       else if (value < -1e8)
1044         value = -1e8;
1045
1046       ratioChi2 += value * value;
1047       fRatioAverage += TMath::Abs(value);
1048
1049       if (ratioChi2 < 0.1)
1050         fLastChi2MCLimit = i;
1051     }
1052     fRatioAverage /= 149;
1053
1054     printf("Sum over (ratio-1)^2 (2..150) is %f; average of |ratio-1| is %f\n", ratioChi2, fRatioAverage);
1055     // TODO FAKE
1056     fLastChi2MC = ratioChi2;
1057
1058     // FFT of ratio
1059     canvas1->cd(6);
1060     const Int_t kFFT = 128;
1061     Double_t fftReal[kFFT];
1062     Double_t fftImag[kFFT];
1063
1064     for (Int_t i=0; i<kFFT; ++i)
1065     {
1066       fftReal[i] = ratio->GetBinContent(i+1+10);
1067       // test: ;-)
1068       //fftReal[i] = cos(TMath::Pi() * 5 * 2 * i / 128);
1069       fftImag[i] = 0;
1070     }
1071
1072     FFT(-1, TMath::Nint(TMath::Log(kFFT) / TMath::Log(2)), fftReal, fftImag);
1073
1074     TH1* fftResult = (TH1*) ratio->Clone("fftResult");
1075     fftResult->SetTitle("FFT;true multiplicity;coeff. (10...137)");
1076     TH1* fftResult2 = (TH1*) ratio->Clone("fftResult2");
1077     TH1* fftResult3 = (TH1*) ratio->Clone("fftResult3");
1078     fftResult->Reset();
1079     fftResult2->Reset();
1080     fftResult3->Reset();
1081     fftResult->GetYaxis()->UnZoom();
1082     fftResult2->GetYaxis()->UnZoom();
1083     fftResult3->GetYaxis()->UnZoom();
1084
1085     //Printf("AFTER FFT");
1086     for (Int_t i=0; i<kFFT; ++i)
1087     {
1088       //Printf("%d: %f", i, fftReal[i]);
1089       fftResult->SetBinContent(i+1, fftReal[i]);
1090       /*if (i != 0 && TMath::Abs(fftReal[i]) > 0.5)
1091       {
1092         Printf("Nulled %d", i);
1093         fftReal[i] = 0;
1094       }*/
1095     }
1096
1097     fftResult->SetLineColor(1);
1098     fftResult->DrawCopy();
1099
1100
1101     //Printf("IMAG");
1102     for (Int_t i=0; i<kFFT; ++i)
1103     {
1104       //Printf("%d: %f", i, fftImag[i]);
1105       fftResult2->SetBinContent(i+1, fftImag[i]);
1106
1107       fftResult3->SetBinContent(i+1, TMath::Sqrt(fftReal[i] * fftReal[i] + fftImag[i] * fftImag[i]));
1108     }
1109
1110     fftResult2->SetLineColor(2);
1111     fftResult2->DrawCopy("SAME");
1112
1113     fftResult3->SetLineColor(4);
1114     fftResult3->DrawCopy("SAME");
1115
1116     for (Int_t i=1; i<kFFT - 1; ++i)
1117     {
1118       if (TMath::Sqrt(fftReal[i] * fftReal[i] + fftImag[i] * fftImag[i]) > 3)
1119       {
1120         fftReal[i-1] = 0;
1121         fftReal[i] = 0;
1122         fftReal[i+1] = 0;
1123         fftImag[i-1] = 0;
1124         fftImag[i] = 0;
1125         fftImag[i+1] = 0;
1126         //fftReal[i] = (fftReal[i-1] + fftReal[i+1]) / 2;
1127         //fftImag[i]  = (fftImag[i-1] + fftImag[i+1]) / 2;
1128         //Printf("Nulled %d to %f %f", i, fftReal[i], fftImag[i]);
1129       }
1130     }
1131
1132     FFT(1, TMath::Nint(TMath::Log(kFFT) / TMath::Log(2)), fftReal, fftImag);
1133
1134     TH1* fftResult4 = (TH1*) fftResult3->Clone("fftResult4");
1135     fftResult4->Reset();
1136
1137     //Printf("AFTER BACK-TRAFO");
1138     for (Int_t i=0; i<kFFT; ++i)
1139     {
1140       //Printf("%d: %f", i, fftReal[i]);
1141       fftResult4->SetBinContent(i+1+10, fftReal[i]);
1142     }
1143
1144     canvas1->cd(5);
1145     fftResult4->SetLineColor(4);
1146     fftResult4->DrawCopy("SAME");
1147
1148     // plot (MC - Unfolded) / error (MC)
1149     canvas1->cd(3);
1150
1151     TH1* diffMCUnfolded2 = dynamic_cast<TH1*> (proj->Clone("diffMCUnfolded2"));
1152     diffMCUnfolded2->Add(fMultiplicityESDCorrected[esdCorrId], -1);
1153
1154     Int_t ndfQual[kQualityRegions];
1155     for (Int_t region=0; region<kQualityRegions; ++region)
1156     {
1157       fQuality[region] = 0;
1158       ndfQual[region] = 0;
1159     }
1160
1161
1162     Double_t newChi2 = 0;
1163     Double_t newChi2Limit150 = 0;
1164     Int_t ndf = 0;
1165     for (Int_t i=1; i<=diffMCUnfolded2->GetNbinsX(); ++i)
1166     {
1167       Double_t value = 0;
1168       if (proj->GetBinError(i) > 0)
1169       {
1170         value = diffMCUnfolded2->GetBinContent(i) / proj->GetBinError(i);
1171         newChi2 += value * value;
1172         if (i > 1 && i <= mcBinLimit)
1173           newChi2Limit150 += value * value;
1174         ++ndf;
1175
1176         for (Int_t region=0; region<kQualityRegions; ++region)
1177           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
1178           {
1179             fQuality[region] += TMath::Abs(value);
1180             ++ndfQual[region];
1181           }
1182       }
1183
1184       diffMCUnfolded2->SetBinContent(i, value);
1185     }
1186
1187     // normalize region to the number of entries
1188     for (Int_t region=0; region<kQualityRegions; ++region)
1189     {
1190       if (ndfQual[region] > 0)
1191         fQuality[region] /= ndfQual[region];
1192       Printf("Quality parameter %d (%d <= mult <= %d) is %f with %d df", region, fgQualityRegionsB[region], fgQualityRegionsE[region], fQuality[region], ndfQual[region]);
1193     }
1194
1195     if (mcBinLimit > 1)
1196     {
1197       // TODO another FAKE
1198       fLastChi2MC = newChi2Limit150 / (mcBinLimit - 1);
1199       Printf("Chi2 (2..%d) from (MC - Unfolded) / e(MC) is: %.2f ndf is %d --> chi2 / ndf = %.2f", mcBinLimit, newChi2Limit150, mcBinLimit - 1, fLastChi2MC);
1200     }
1201     else
1202       fLastChi2MC = -1;
1203
1204     Printf("Chi2 (full range) from (MC - Unfolded) / e(MC) is: %.2f ndf is %d --> chi2 / ndf = %.2f", newChi2, ndf, ((ndf > 0) ? newChi2 / ndf : -1));
1205
1206     diffMCUnfolded2->SetTitle("#chi^{2};Npart;(MC - Unfolded) / e(MC)");
1207     diffMCUnfolded2->GetYaxis()->SetRangeUser(-5, 5);
1208     diffMCUnfolded2->GetXaxis()->SetRangeUser(0, 200);
1209     diffMCUnfolded2->DrawCopy("HIST");
1210   }
1211
1212   canvas1->SaveAs(Form("%s.gif", canvas1->GetName()));
1213 }
1214
1215 //____________________________________________________________________
1216 void AliMultiplicityCorrection::FFT(Int_t dir, Int_t m, Double_t *x, Double_t *y)
1217 {
1218 /*-------------------------------------------------------------------------
1219    This computes an in-place complex-to-complex FFT
1220    x and y are the real and imaginary arrays of 2^m points.
1221    dir =  1 gives forward transform
1222    dir = -1 gives reverse transform
1223
1224      Formula: forward
1225                   N-1
1226                   ---
1227               1   \          - j k 2 pi n / N
1228       X(n) = ---   >   x(k) e                    = forward transform
1229               N   /                                n=0..N-1
1230                   ---
1231                   k=0
1232
1233       Formula: reverse
1234                   N-1
1235                   ---
1236                   \          j k 2 pi n / N
1237       X(n) =       >   x(k) e                    = forward transform
1238                   /                                n=0..N-1
1239                   ---
1240                   k=0
1241 */
1242
1243    Long_t   nn, i, i1, j, k, i2, l, l1, l2;
1244    Double_t c1, c2, tx, ty, t1, t2, u1, u2, z;
1245
1246    /* Calculate the number of points */
1247    nn = 1;
1248    for (i = 0; i < m; i++)
1249        nn *= 2;
1250
1251    /* Do the bit reversal */
1252    i2 = nn >> 1;
1253    j = 0;
1254    for (i= 0; i < nn - 1; i++) {
1255        if (i < j) {
1256            tx = x[i];
1257            ty = y[i];
1258            x[i] = x[j];
1259            y[i] = y[j];
1260            x[j] = tx;
1261            y[j] = ty;
1262        }
1263        k = i2;
1264        while (k <= j) {
1265            j -= k;
1266            k >>= 1;
1267        }
1268        j += k;
1269    }
1270
1271    /* Compute the FFT */
1272    c1 = -1.0;
1273    c2 = 0.0;
1274    l2 = 1;
1275    for (l = 0; l < m; l++) {
1276        l1 = l2;
1277        l2 <<= 1;
1278        u1 = 1.0;
1279        u2 = 0.0;
1280        for (j = 0;j < l1; j++) {
1281            for (i = j; i < nn; i += l2) {
1282                i1 = i + l1;
1283                t1 = u1 * x[i1] - u2 * y[i1];
1284                t2 = u1 * y[i1] + u2 * x[i1];
1285                x[i1] = x[i] - t1;
1286                y[i1] = y[i] - t2;
1287                x[i] += t1;
1288                y[i] += t2;
1289            }
1290            z =  u1 * c1 - u2 * c2;
1291            u2 = u1 * c2 + u2 * c1;
1292            u1 = z;
1293        }
1294        c2 = TMath::Sqrt((1.0 - c1) / 2.0);
1295        if (dir == 1)
1296            c2 = -c2;
1297        c1 = TMath::Sqrt((1.0 + c1) / 2.0);
1298    }
1299
1300    /* Scaling for forward transform */
1301    if (dir == 1) {
1302        for (i=0;i<nn;i++) {
1303            x[i] /= (Double_t)nn;
1304            y[i] /= (Double_t)nn;
1305        }
1306    }
1307 }
1308
1309 //____________________________________________________________________
1310 void AliMultiplicityCorrection::GetComparisonResults(Float_t* mc, Int_t* mcLimit, Float_t* residuals, Float_t* ratioAverage) const
1311 {
1312   // Returns the chi2 between the MC and the unfolded ESD as well as between the ESD and the folded unfolded ESD
1313   // These values are computed during DrawComparison, thus this function picks up the
1314   // last calculation
1315
1316   if (mc)
1317     *mc = fLastChi2MC;
1318   if (mcLimit)
1319     *mcLimit = fLastChi2MCLimit;
1320   if (residuals)
1321     *residuals = fLastChi2Residuals;
1322   if (ratioAverage)
1323     *ratioAverage = fRatioAverage;
1324 }
1325
1326 //____________________________________________________________________
1327 TH2F* AliMultiplicityCorrection::GetMultiplicityMC(Int_t i, EventType eventType)
1328 {
1329   //
1330   // returns the corresponding MC spectrum
1331   //
1332
1333   switch (eventType)
1334   {
1335     case kTrVtx : return fMultiplicityVtx[i]; break;
1336     case kMB : return fMultiplicityMB[i]; break;
1337     case kINEL : return fMultiplicityINEL[i]; break;
1338     case kNSD : return fMultiplicityNSD[i]; break;
1339   }
1340
1341   return 0;
1342 }
1343
1344 //____________________________________________________________________
1345 void AliMultiplicityCorrection::SetMultiplicityMC(Int_t i, EventType eventType, TH2F* hist)
1346 {
1347   //
1348   // returns the corresponding MC spectrum
1349   //
1350
1351   switch (eventType)
1352   {
1353     case kTrVtx : fMultiplicityVtx[i] = hist; break;
1354     case kMB : fMultiplicityMB[i] = hist; break;
1355     case kINEL : fMultiplicityINEL[i] = hist; break;
1356     case kNSD : fMultiplicityNSD[i] = hist; break;
1357   }
1358 }
1359
1360 //____________________________________________________________________
1361 TH1* AliMultiplicityCorrection::CalculateStdDev(TH1** results, Int_t max)
1362 {
1363   // calculate standard deviation of (results[0] - results[k]) k=1...max-1
1364   // per bin one gets: sigma(r[0] - r[n]) / r[0]
1365
1366   TH1* standardDeviation = (TH1*) results[0]->Clone("standardDeviation");
1367   standardDeviation->Reset();
1368
1369   for (Int_t x=1; x<=results[0]->GetNbinsX(); x++)
1370   {
1371     if (results[0]->GetBinContent(x) > 0)
1372     {
1373       Double_t average = 0;
1374       for (Int_t n=1; n<max; ++n)
1375         average += results[n]->GetBinContent(x) - results[0]->GetBinContent(x);
1376       average /= max-1;
1377
1378       Double_t variance = 0;
1379       for (Int_t n=1; n<max; ++n)
1380       {
1381         Double_t value = results[n]->GetBinContent(x) - results[0]->GetBinContent(x) - average;
1382         variance += value * value;
1383       }
1384       variance /= max-1;
1385
1386       Double_t standardDev = TMath::Sqrt(variance);
1387       standardDeviation->SetBinContent(x, standardDev / results[0]->GetBinContent(x));
1388       //Printf("sigma_%d is %f value %f --> error %f", x, standardDev, results[0]->GetBinContent(x), standardDev / results[0]->GetBinContent(x));
1389     }
1390   }
1391
1392   return standardDeviation;
1393 }
1394
1395 //____________________________________________________________________
1396 TH1* AliMultiplicityCorrection::StatisticalUncertainty(AliUnfolding::MethodType methodType, Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t randomizeMeasured, Bool_t randomizeResponse, TH1* compareTo)
1397 {
1398   //
1399   // evaluates the uncertainty that arises from the non-infinite statistics in the response matrix
1400   // the function unfolds the spectrum using the default response matrix and several modified ones
1401   // the modified ones are created by randomizing each cell using poisson statistics with the mean = bin value
1402   // these unfolded results are compared to the first result gained with the default response OR to the histogram given
1403   // in <compareTo> (optional)
1404   //
1405   // returns the error assigned to the measurement
1406   //
1407
1408   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
1409
1410   // initialize seed with current time
1411   gRandom->SetSeed(0);
1412
1413   if (methodType == AliUnfolding::kChi2Minimization)
1414     AliUnfolding::SetCreateOverflowBin(5);
1415   AliUnfolding::SetUnfoldingMethod(methodType);
1416
1417   const Int_t kErrorIterations = 150;
1418
1419   TH1* maxError = 0;
1420   TH1* firstResult = 0;
1421
1422   TH1** results = new TH1*[kErrorIterations];
1423
1424   for (Int_t n=0; n<kErrorIterations; ++n)
1425   {
1426     Printf("Iteration %d of %d...", n, kErrorIterations);
1427
1428     SetupCurrentHists(inputRange, fullPhaseSpace, eventType);
1429
1430     TH1* measured = (TH1*) fCurrentESD->Clone("measured");
1431
1432     if (n > 0)
1433     {
1434       if (randomizeResponse)
1435       {
1436         // randomize response matrix
1437         for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
1438           for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
1439             fCurrentCorrelation->SetBinContent(i, j, gRandom->Poisson(fCurrentCorrelation->GetBinContent(i, j)));
1440       }
1441
1442       if (randomizeMeasured)
1443       {
1444         // randomize measured spectrum
1445         for (Int_t x=1; x<=measured->GetNbinsX(); x++) // mult. axis
1446         {
1447           Int_t randomValue = gRandom->Poisson(fCurrentESD->GetBinContent(x));
1448           measured->SetBinContent(x, randomValue);
1449           measured->SetBinError(x, TMath::Sqrt(randomValue));
1450         }
1451       }
1452     }
1453
1454     // only for bayesian method we have to do it before the call to Unfold...
1455     if (methodType == AliUnfolding::kBayesian)
1456     {
1457       for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
1458       {
1459         // with this it is normalized to 1
1460         Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
1461
1462         // with this normalized to the given efficiency
1463         if (fCurrentEfficiency->GetBinContent(i) > 0)
1464           sum /= fCurrentEfficiency->GetBinContent(i);
1465         else
1466           sum = 0;
1467
1468         for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
1469         {
1470           if (sum > 0)
1471           {
1472             fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
1473             fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
1474           }
1475           else
1476           {
1477             fCurrentCorrelation->SetBinContent(i, j, 0);
1478             fCurrentCorrelation->SetBinError(i, j, 0);
1479           }
1480         }
1481       }
1482     }
1483
1484     TH1* result = 0;
1485     if (n == 0 && compareTo)
1486     {
1487       // in this case we just store the histogram we want to compare to
1488       result = (TH1*) compareTo->Clone("compareTo");
1489       result->Sumw2();
1490     }
1491     else
1492     {
1493       result = (TH1*) fMultiplicityESDCorrected[correlationID]->Clone(Form("result_%d", n));
1494
1495       Int_t returnCode = AliUnfolding::Unfold(fCurrentCorrelation, fCurrentEfficiency, measured, 0, result);
1496
1497       if (returnCode != 0)
1498         return 0;
1499     }
1500
1501     // normalize
1502     result->Scale(1.0 / result->Integral());
1503
1504     if (n == 0)
1505     {
1506       firstResult = (TH1*) result->Clone("firstResult");
1507
1508       maxError = (TH1*) result->Clone("maxError");
1509       maxError->Reset();
1510     }
1511     else
1512     {
1513       // calculate ratio
1514       TH1* ratio = (TH1*) firstResult->Clone("ratio");
1515       ratio->Divide(result);
1516
1517       // find max. deviation
1518       for (Int_t x=1; x<=ratio->GetNbinsX(); x++)
1519         maxError->SetBinContent(x, TMath::Max(maxError->GetBinContent(x), TMath::Abs(1 - ratio->GetBinContent(x))));
1520
1521       delete ratio;
1522     }
1523
1524     results[n] = result;
1525   }
1526
1527   // find covariance matrix
1528   // results[n] is X_x
1529   // cov. matrix is M_xy = E ( (X_x - E(X_x)) * (X_y - E(X_y))), with E() = expectation value
1530
1531   Int_t nBins = results[0]->GetNbinsX();
1532   Float_t lowEdge = results[0]->GetXaxis()->GetBinLowEdge(1);
1533   Float_t upEdge = results[0]->GetXaxis()->GetBinUpEdge(nBins);
1534
1535   // find average, E(X)
1536   TProfile* average = new TProfile("average", "average", nBins, lowEdge, upEdge);
1537   for (Int_t n=1; n<kErrorIterations; ++n)
1538     for (Int_t x=1; x<=results[n]->GetNbinsX(); x++)
1539       average->Fill(results[n]->GetXaxis()->GetBinCenter(x), results[n]->GetBinContent(x));
1540   //new TCanvas; average->DrawClone();
1541
1542   // find cov. matrix
1543   TProfile2D* covMatrix = new TProfile2D("covMatrix", "covMatrix", nBins, lowEdge, upEdge, nBins, lowEdge, upEdge);
1544
1545   for (Int_t n=1; n<kErrorIterations; ++n)
1546     for (Int_t x=1; x<=results[n]->GetNbinsX(); x++)
1547       for (Int_t y=1; y<=results[n]->GetNbinsX(); y++)
1548       {
1549         // (X_x - E(X_x)) * (X_y - E(X_y)
1550         Float_t cov = (results[n]->GetBinContent(x) - average->GetBinContent(x)) * (results[n]->GetBinContent(y) - average->GetBinContent(y));
1551         covMatrix->Fill(results[n]->GetXaxis()->GetBinCenter(x), results[n]->GetXaxis()->GetBinCenter(y), cov);
1552       }
1553   TCanvas* c = new TCanvas; c->cd(); covMatrix->DrawCopy("COLZ");
1554
1555 //   // fill 2D histogram that contains deviation from first
1556 //   TH2F* deviations = new TH2F("deviations", "deviations", nBins, lowEdge, upEdge, 1000, -0.01, 0.01);
1557 //   for (Int_t n=1; n<kErrorIterations; ++n)
1558 //     for (Int_t x=1; x<=results[n]->GetNbinsX(); x++)
1559 //       deviations->Fill(results[n]->GetXaxis()->GetBinCenter(x), results[n]->GetBinContent(x) - results[0]->GetBinContent(x));
1560 //   //new TCanvas; deviations->DrawCopy("COLZ");
1561 //
1562 //   // get standard deviation "by hand"
1563 //   for (Int_t x=1; x<=nBins; x++)
1564 //   {
1565 //     if (results[0]->GetBinContent(x) > 0)
1566 //     {
1567 //       TH1* proj = deviations->ProjectionY("projRMS", x, x, "e");
1568 //       Float_t standardDev = proj->GetRMS(); // this is standard deviation in fact
1569 //       //standardDeviation->SetBinContent(x, standardDev / results[0]->GetBinContent(x));
1570 //       Printf("sigma_%d is %f value %f --> error %f", x, standardDev, results[0]->GetBinContent(x), standardDev / results[0]->GetBinContent(x));
1571 //     }
1572 //   }
1573
1574   TH1* standardDeviation = CalculateStdDev(results, kErrorIterations);
1575
1576   // compare maxError to RMS of profile (variable name: average)
1577   // first: calculate rms in percent of value
1578   TH1* rmsError = (TH1*) maxError->Clone("rmsError");
1579   rmsError->Reset();
1580
1581   // enable error to be standard deviation (see http://root.cern.ch/root/html/TProfile.html#TProfile:SetErrorOption)
1582   average->SetErrorOption("s");
1583   for (Int_t x=1; x<=rmsError->GetNbinsX(); x++)
1584     if (average->GetBinContent(x) > 0)
1585       rmsError->SetBinContent(x, average->GetBinError(x) / average->GetBinContent(x));
1586
1587   // find maxError deviation from average (not from "first result"/mc as above)
1588   TH1* maxError2 = (TH1*) maxError->Clone("maxError2");
1589   maxError2->Reset();
1590   for (Int_t n=1; n<kErrorIterations; ++n)
1591     for (Int_t x=1; x<=results[n]->GetNbinsX(); x++)
1592       if (average->GetBinContent(x) > 0)
1593         maxError2->SetBinContent(x, TMath::Max(maxError2->GetBinContent(x), TMath::Abs((results[n]->GetBinContent(x) - average->GetBinContent(x)) / average->GetBinContent(x))));
1594
1595   //new TCanvas; maxError2->DrawCopy(); rmsError->SetLineColor(2); rmsError->DrawCopy("SAME"); standardDeviation->SetLineColor(3); standardDeviation->DrawCopy("SAME");
1596
1597   // plot difference between average and MC/first result
1598   TH1* averageFirstRatio = (TH1*) results[0]->Clone("averageFirstRatio");
1599   averageFirstRatio->Reset();
1600   averageFirstRatio->Divide(results[0], average);
1601
1602   //new TCanvas; results[0]->DrawCopy(); average->SetLineColor(2); average->DrawClone("SAME");
1603   //new TCanvas; averageFirstRatio->DrawCopy();
1604
1605   static TH1* temp = 0;
1606   if (!temp)
1607   {
1608     temp = (TH1*) standardDeviation->Clone("temp");
1609     for (Int_t x=1; x<=results[0]->GetNbinsX(); x++)
1610       temp->SetBinContent(x, temp->GetBinContent(x) * results[0]->GetBinContent(x));
1611   }
1612   else
1613   {
1614     // find difference from result[0] as TH2
1615     TH2F* pulls = new TH2F("pulls", "pulls;multiplicity;difference", nBins, lowEdge, upEdge, 1000, -10, 10);
1616     for (Int_t n=1; n<kErrorIterations; ++n)
1617       for (Int_t x=1; x<=results[n]->GetNbinsX(); x++)
1618         if (temp->GetBinContent(x) > 0)
1619           pulls->Fill(results[n]->GetXaxis()->GetBinCenter(x), (results[0]->GetBinContent(x) - results[n]->GetBinContent(x)) / temp->GetBinContent(x));
1620     new TCanvas("pulls", "pulls", 800, 600); pulls->DrawCopy(); pulls->FitSlicesY();
1621   }
1622
1623   // clean up
1624   for (Int_t n=0; n<kErrorIterations; ++n)
1625     delete results[n];
1626   delete[] results;
1627
1628   // fill into result histogram
1629   for (Int_t i=1; i<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); ++i)
1630     fMultiplicityESDCorrected[correlationID]->SetBinContent(i, firstResult->GetBinContent(i));
1631
1632   for (Int_t i=1; i<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); ++i)
1633     fMultiplicityESDCorrected[correlationID]->SetBinError(i, maxError->GetBinContent(i) * fMultiplicityESDCorrected[correlationID]->GetBinContent(i));
1634
1635   return standardDeviation;
1636 }
1637
1638 //____________________________________________________________________
1639 void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Float_t regPar, Int_t nIterations, TH1* initialConditions, Bool_t determineError)
1640 {
1641   //
1642   // correct spectrum using bayesian method
1643   //
1644
1645   // initialize seed with current time
1646   gRandom->SetSeed(0);
1647
1648   SetupCurrentHists(inputRange, fullPhaseSpace, eventType);
1649
1650   // normalize correction for given nPart
1651   for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
1652   {
1653     // with this it is normalized to 1
1654     Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
1655
1656     // with this normalized to the given efficiency
1657     if (fCurrentEfficiency->GetBinContent(i) > 0)
1658       sum /= fCurrentEfficiency->GetBinContent(i);
1659     else
1660       sum = 0;
1661
1662     for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
1663     {
1664       if (sum > 0)
1665       {
1666         fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
1667         fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
1668       }
1669       else
1670       {
1671         fCurrentCorrelation->SetBinContent(i, j, 0);
1672         fCurrentCorrelation->SetBinError(i, j, 0);
1673       }
1674     }
1675   }
1676
1677   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
1678
1679   AliUnfolding::SetBayesianParameters(regPar, nIterations);
1680   AliUnfolding::SetUnfoldingMethod(AliUnfolding::kBayesian);
1681   if (AliUnfolding::Unfold(fCurrentCorrelation, fCurrentEfficiency, fCurrentESD, initialConditions, fMultiplicityESDCorrected[correlationID]) != 0)
1682     return;
1683
1684   if (!determineError)
1685   {
1686     Printf("AliMultiplicityCorrection::ApplyBayesianMethod: WARNING: No errors calculated.");
1687     return;
1688   }
1689
1690   // evaluate errors, this is done by randomizing the measured spectrum following Poission statistics
1691   // this (new) measured spectrum is then unfolded and the different to the result from the "real" measured
1692   // spectrum calculated. This is performed N times and the sigma is taken as the statistical
1693   // error of the unfolding method itself.
1694
1695   const Int_t kErrorIterations = 20;
1696
1697   printf("Spectrum unfolded. Determining error (%d iterations)...\n", kErrorIterations);
1698
1699   TH1* randomized = (TH1*) fCurrentESD->Clone("randomized");
1700   TH1* resultArray[kErrorIterations+1];
1701   for (Int_t n=0; n<kErrorIterations; ++n)
1702   {
1703     // randomize the content of clone following a poisson with the mean = the value of that bin
1704     for (Int_t x=1; x<=randomized->GetNbinsX(); x++) // mult. axis
1705     {
1706       Int_t randomValue = gRandom->Poisson(fCurrentESD->GetBinContent(x));
1707       //printf("%d --> %d\n", fCurrentESD->GetBinContent(x), randomValue);
1708       randomized->SetBinContent(x, randomValue);
1709       randomized->SetBinError(x, TMath::Sqrt(randomValue));
1710     }
1711
1712     TH1* result2 = (TH1*) fMultiplicityESDCorrected[correlationID]->Clone("result2");
1713     result2->Reset();
1714     if (AliUnfolding::Unfold(fCurrentCorrelation, fCurrentEfficiency, randomized, initialConditions, result2) != 0)
1715       return;
1716
1717     result2->Scale(1.0 / result2->Integral());
1718
1719     resultArray[n+1] = result2;
1720   }
1721   delete randomized;
1722
1723   resultArray[0] = fMultiplicityESDCorrected[correlationID];
1724   TH1* error = CalculateStdDev(resultArray, kErrorIterations+1);
1725
1726   for (Int_t n=0; n<kErrorIterations; ++n)
1727     delete resultArray[n+1];
1728
1729   for (Int_t i=1; i<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); ++i)
1730     fMultiplicityESDCorrected[correlationID]->SetBinError(i, error->GetBinContent(i) * fMultiplicityESDCorrected[correlationID]->GetBinContent(i));
1731
1732   delete error;
1733 }
1734
1735 //____________________________________________________________________
1736 Float_t AliMultiplicityCorrection::BayesCovarianceDerivate(Float_t matrixM[251][251], TH2* hResponse, Int_t k, Int_t i, Int_t r, Int_t u)
1737 {
1738   //
1739   // helper function for the covariance matrix of the bayesian method
1740   //
1741
1742   Float_t result = 0;
1743
1744   if (k == u && r == i)
1745     result += 1.0 / hResponse->GetBinContent(u+1, r+1);
1746
1747   if (k == u)
1748     result -= 1.0 / fCurrentEfficiency->GetBinContent(u+1);
1749
1750   if (r == i)
1751     result -= matrixM[u][i] * fCurrentEfficiency->GetBinContent(u+1) / hResponse->GetBinContent(u+1, i+1);
1752
1753   result *= matrixM[k][i];
1754
1755   return result;
1756 }
1757
1758 //____________________________________________________________________
1759 void AliMultiplicityCorrection::ApplyLaszloMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType)
1760 {
1761   //
1762   // correct spectrum using bayesian method
1763   //
1764
1765   Float_t regPar = 0;
1766
1767   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
1768   Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4);
1769
1770   SetupCurrentHists(inputRange, fullPhaseSpace, eventType);
1771   //normalize ESD
1772   fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
1773
1774   // TODO should be taken from correlation map
1775   //TH1* sumHist = GetMultiplicityMC(inputRange, eventType)->ProjectionY("sumHist", 1, GetMultiplicityMC(inputRange, eventType)->GetNbinsX());
1776
1777   // normalize correction for given nPart
1778   for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
1779   {
1780     Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
1781     //Double_t sum = sumHist->GetBinContent(i);
1782     if (sum <= 0)
1783       continue;
1784     for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
1785     {
1786       // npart sum to 1
1787       fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
1788       fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
1789     }
1790   }
1791
1792   new TCanvas;
1793   fCurrentCorrelation->Draw("COLZ");
1794
1795   // FAKE
1796   fCurrentEfficiency = ((TH2*) fCurrentCorrelation)->ProjectionX("eff");
1797
1798   // pick prior distribution
1799   TH1F* hPrior = (TH1F*)fCurrentESD->Clone("prior");
1800   Float_t norm = 1; //hPrior->Integral();
1801   for (Int_t t=1; t<=hPrior->GetNbinsX(); t++)
1802     hPrior->SetBinContent(t, hPrior->GetBinContent(t)/norm);
1803
1804   // zero distribution
1805   TH1F* zero =  (TH1F*)hPrior->Clone("zero");
1806
1807   // define temp hist
1808   TH1F* hTemp = (TH1F*)fCurrentESD->Clone("temp");
1809   hTemp->Reset();
1810
1811   // just a shortcut
1812   TH2F* hResponse = (TH2F*) fCurrentCorrelation;
1813
1814   // unfold...
1815   Int_t iterations = 25;
1816   for (Int_t i=0; i<iterations; i++)
1817   {
1818     //printf(" iteration %i \n", i);
1819
1820     for (Int_t m=1; m<=hResponse->GetNbinsY(); m++)
1821     {
1822       Float_t value = 0;
1823       for (Int_t t = 1; t<=hResponse->GetNbinsX(); t++)
1824         value += hResponse->GetBinContent(t, m) * hPrior->GetBinContent(t);
1825       hTemp->SetBinContent(m, value);
1826       //printf("%d %f %f %f\n", m, zero->GetBinContent(m), hPrior->GetBinContent(m), value);
1827     }
1828
1829     // regularization (simple smoothing)
1830     TH1F* hTrueSmoothed = (TH1F*) hTemp->Clone("truesmoothed");
1831
1832     for (Int_t t=2; t<hTrueSmoothed->GetNbinsX(); t++)
1833     {
1834       Float_t average = (hTemp->GetBinContent(t-1) / hTemp->GetBinWidth(t-1)
1835                          + hTemp->GetBinContent(t) / hTemp->GetBinWidth(t)
1836                          + hTemp->GetBinContent(t+1) / hTemp->GetBinWidth(t+1)) / 3.;
1837       average *= hTrueSmoothed->GetBinWidth(t);
1838
1839       // weight the average with the regularization parameter
1840       hTrueSmoothed->SetBinContent(t, (1-regPar) * hTemp->GetBinContent(t) + regPar * average);
1841     }
1842
1843     for (Int_t m=1; m<=hResponse->GetNbinsY(); m++)
1844       hTemp->SetBinContent(m, zero->GetBinContent(m) + hPrior->GetBinContent(m) - hTrueSmoothed->GetBinContent(m));
1845
1846     // fill guess
1847     for (Int_t t=1; t<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); t++)
1848     {
1849       fMultiplicityESDCorrected[correlationID]->SetBinContent(t, hTemp->GetBinContent(t));
1850       fMultiplicityESDCorrected[correlationID]->SetBinError(t, 0.05 * hTemp->GetBinContent(t)); // TODO
1851
1852       //printf(" bin %d content %f \n", t, fMultiplicityESDCorrected[correlationID]->GetBinContent(t));
1853     }
1854
1855
1856     // calculate chi2 (change from last iteration)
1857     Double_t chi2 = 0;
1858
1859     // use smoothed true (normalized) as new prior
1860     norm = 1; //hTrueSmoothed->Integral();
1861
1862     for (Int_t t=1; t<hTrueSmoothed->GetNbinsX(); t++)
1863     {
1864       Float_t newValue = hTemp->GetBinContent(t)/norm;
1865       Float_t diff = hPrior->GetBinContent(t) - newValue;
1866       chi2 += (Double_t) diff * diff;
1867
1868       hPrior->SetBinContent(t, newValue);
1869     }
1870
1871     printf("Chi2 of %d iteration = %.10f\n", i, chi2);
1872
1873     //if (i % 5 == 0)
1874       DrawComparison(Form("Laszlo_%d", i), inputRange, fullPhaseSpace, kTRUE, GetMultiplicityMC(mcTarget, eventType)->ProjectionY());
1875
1876     delete hTrueSmoothed;
1877   } // end of iterations
1878
1879   DrawComparison("Laszlo", inputRange, fullPhaseSpace, kTRUE, GetMultiplicityMC(mcTarget, eventType)->ProjectionY());
1880 }
1881
1882 //____________________________________________________________________
1883 void AliMultiplicityCorrection::ApplyGaussianMethod(Int_t inputRange, Bool_t fullPhaseSpace)
1884 {
1885   //
1886   // correct spectrum using a simple Gaussian approach, that is model-dependent
1887   //
1888
1889   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
1890   Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4);
1891
1892   SetupCurrentHists(inputRange, fullPhaseSpace, kTrVtx);
1893   //normalize ESD
1894   fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
1895
1896   TH1D* correction = dynamic_cast<TH1D*> (fCurrentESD->Clone("GaussianMean"));
1897   correction->SetTitle("GaussianMean");
1898
1899   TH1D* correctionWidth = dynamic_cast<TH1D*> (fCurrentESD->Clone("GaussianWidth"));
1900   correction->SetTitle("GaussianWidth");
1901
1902   for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
1903   {
1904     TH1D* proj = (dynamic_cast<TH2*> (fCurrentCorrelation))->ProjectionX("_px", i, i+1);
1905     proj->Fit("gaus", "0Q");
1906     correction->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(1));
1907     correctionWidth->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(2));
1908
1909     continue;
1910
1911     // draw for debugging
1912     new TCanvas;
1913     proj->DrawCopy();
1914     proj->GetFunction("gaus")->DrawCopy("SAME");
1915   }
1916
1917   TH1* target = fMultiplicityESDCorrected[correlationID];
1918
1919   Int_t nBins = target->GetNbinsX()*10+1;
1920   Float_t* binning = new Float_t[nBins];
1921   for (Int_t i=1; i<=target->GetNbinsX(); ++i)
1922     for (Int_t j=0; j<10; ++j)
1923       binning[(i-1)*10 + j] = target->GetXaxis()->GetBinLowEdge(i) + target->GetXaxis()->GetBinWidth(i) / 10 * j;
1924
1925   binning[nBins-1] = target->GetXaxis()->GetBinUpEdge(target->GetNbinsX());
1926
1927   TH1F* fineBinned = new TH1F("targetFineBinned", "targetFineBinned", nBins-1, binning);
1928
1929   for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
1930   {
1931     Float_t mean = correction->GetBinContent(i);
1932     Float_t width = correctionWidth->GetBinContent(i);
1933
1934     Int_t fillBegin = fineBinned->FindBin(mean - width * 5);
1935     Int_t fillEnd   = fineBinned->FindBin(mean + width * 5);
1936     //printf("bin %d mean %f width %f, filling from %d to %d\n", i, mean, width, fillBegin, fillEnd);
1937
1938     for (Int_t j=fillBegin; j <= fillEnd; ++j)
1939     {
1940       fineBinned->AddBinContent(j, TMath::Gaus(fineBinned->GetXaxis()->GetBinCenter(j), mean, width, kTRUE) * fCurrentESD->GetBinContent(i));
1941     }
1942   }
1943
1944   for (Int_t i=1; i<=target->GetNbinsX(); ++i)
1945   {
1946     Float_t sum = 0;
1947     for (Int_t j=1; j<=10; ++j)
1948       sum += fineBinned->GetBinContent((i-1)*10 + j);
1949     target->SetBinContent(i, sum / 10);
1950   }
1951
1952   delete[] binning;
1953
1954   DrawComparison("Gaussian", inputRange, fullPhaseSpace, kFALSE, GetMultiplicityMC(mcTarget, kTrVtx)->ProjectionY());
1955 }
1956
1957 //____________________________________________________________________
1958 TH2F* AliMultiplicityCorrection::CalculateMultiplicityESD(TH1* inputMC, Int_t correlationMap)
1959 {
1960   // runs the distribution given in inputMC through the response matrix identified by
1961   // correlationMap and produces a measured distribution
1962   // although it is a TH2F the vertex axis is not used at the moment and all entries are filled in mid-vertex
1963   // if normalized is set, inputMC is expected to be normalized to the bin width
1964
1965   if (!inputMC)
1966     return 0;
1967
1968   if (correlationMap < 0 || correlationMap >= kCorrHists)
1969     return 0;
1970
1971   // empty under/overflow bins in x, otherwise Project3D takes them into account
1972   TH3* hist = fCorrelation[correlationMap];
1973   for (Int_t y=0; y<=hist->GetYaxis()->GetNbins()+1; ++y)
1974   {
1975     for (Int_t z=0; z<=hist->GetZaxis()->GetNbins()+1; ++z)
1976     {
1977       hist->SetBinContent(0, y, z, 0);
1978       hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0);
1979     }
1980   }
1981
1982   TH2* corr = (TH2*) hist->Project3D("zy");
1983   //corr->Rebin2D(2, 1);
1984   corr->Sumw2();
1985   Printf("Correction histogram used for convolution has %f entries", corr->Integral());
1986
1987   // normalize correction for given nPart
1988   for (Int_t i=1; i<=corr->GetNbinsX(); ++i)
1989   {
1990     Double_t sum = corr->Integral(i, i, 1, corr->GetNbinsY());
1991     if (sum <= 0)
1992       continue;
1993
1994     for (Int_t j=1; j<=corr->GetNbinsY(); ++j)
1995     {
1996       // npart sum to 1
1997       corr->SetBinContent(i, j, corr->GetBinContent(i, j) / sum);
1998       corr->SetBinError(i, j, corr->GetBinError(i, j) / sum);
1999     }
2000   }
2001
2002   TH2F* target = dynamic_cast<TH2F*> (fMultiplicityESD[0]->Clone(Form("%s_measured", inputMC->GetName())));
2003   target->Reset();
2004
2005   for (Int_t meas=1; meas<=corr->GetNbinsY(); ++meas)
2006   {
2007     Float_t measured = 0;
2008     Float_t error = 0;
2009
2010     for (Int_t gen=1; gen<=corr->GetNbinsX(); ++gen)
2011     {
2012       Int_t mcGenBin = inputMC->GetXaxis()->FindBin(corr->GetXaxis()->GetBinCenter(gen));
2013
2014       measured += inputMC->GetBinContent(mcGenBin) * corr->GetBinContent(gen, meas);
2015       error += inputMC->GetBinError(mcGenBin) * corr->GetBinContent(gen, meas);
2016     }
2017
2018     //printf("%f +- %f ; %f +- %f \n", inputMC->GetBinContent(meas), inputMC->GetBinError(meas), measured, error);
2019
2020     target->SetBinContent(1 + target->GetNbinsX() / 2, meas, measured);
2021     target->SetBinError(1 + target->GetNbinsX() / 2, meas, error);
2022   }
2023
2024   return target;
2025 }
2026
2027 //____________________________________________________________________
2028 void AliMultiplicityCorrection::SetGenMeasFromFunc(TF1* inputMC, Int_t id)
2029 {
2030   // uses the given function to fill the input MC histogram and generates from that
2031   // the measured histogram by applying the response matrix
2032   // this can be used to evaluate if the methods work indepedently of the input
2033   // distribution
2034   // WARNING does not respect the vertex distribution, just fills central vertex bin
2035
2036   if (!inputMC)
2037     return;
2038
2039   if (id < 0 || id >= kESDHists)
2040     return;
2041
2042   // fill histogram used for random generation
2043   TH1* tmp = fMultiplicityVtx[id]->ProjectionY("tmp");
2044   tmp->Reset();
2045
2046   for (Int_t i=1; i<=tmp->GetNbinsX(); ++i)
2047     tmp->SetBinContent(i, inputMC->Eval(tmp->GetXaxis()->GetBinCenter(i)) * tmp->GetXaxis()->GetBinWidth(i));
2048
2049   TH1* mcRnd = fMultiplicityVtx[id]->ProjectionY("mcRnd");
2050   mcRnd->Reset();
2051   mcRnd->FillRandom(tmp, tmp->Integral());
2052
2053   //new TCanvas; tmp->Draw();
2054   //new TCanvas; mcRnd->Draw();
2055
2056   // and move into 2d histogram
2057   TH1* mc = fMultiplicityVtx[id];
2058   mc->Reset();
2059   for (Int_t i=1; i<=mc->GetNbinsY(); ++i)
2060   {
2061     mc->SetBinContent(mc->GetNbinsX() / 2 + 1, i, mcRnd->GetBinContent(i));
2062     mc->SetBinError(mc->GetNbinsX() / 2 + 1, i, TMath::Sqrt(mcRnd->GetBinContent(i)));
2063   }
2064
2065   //new TCanvas; mc->Draw("COLZ");
2066
2067   // now randomize the measured histogram; funcMeasured is used as pilot function to generated the measured entries
2068   TH1* funcMeasured = CalculateMultiplicityESD(tmp, id)->ProjectionY("funcMeasured");
2069
2070   //new TCanvas; funcMeasured->Draw();
2071
2072   fMultiplicityESD[id]->Reset();
2073
2074   TH1* measRnd = fMultiplicityESD[id]->ProjectionY("measRnd");
2075   measRnd->FillRandom(funcMeasured, tmp->Integral());
2076
2077   //new TCanvas; measRnd->Draw();
2078
2079   fMultiplicityESD[id]->Reset();
2080   for (Int_t i=1; i<=fMultiplicityESD[id]->GetNbinsY(); ++i)
2081   {
2082     fMultiplicityESD[id]->SetBinContent(fMultiplicityESD[id]->GetNbinsX() / 2 + 1, i, measRnd->GetBinContent(i));
2083     fMultiplicityESD[id]->SetBinError(fMultiplicityESD[id]->GetNbinsX() / 2 + 1, i, TMath::Sqrt(measRnd->GetBinContent(i)));
2084   }
2085 }