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