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