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