PWGCF/Correlations/Base/AliUEHist.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: AliUEHist.cxx 20164 2007-08-14 15:31:50Z morsch $ */
  17
  18 //
  19 //
  20 // encapsulate histogram and corrections for one underlying event histogram
  21 //
  22 //
  23 // Author: Jan Fiete Grosse-Oetringhaus, Sara Vallero
  24
  25 #include "AliUEHist.h"
  26 #include "AliCFContainer.h"
  27 #include "THnSparse.h"
  28 #include "TMath.h"
  29 #include "TList.h"
  30 #include "TCollection.h"
  31 #include "TH1D.h"
  32 #include "TH2D.h"
  33 #include "TH3D.h"
  34 #include "AliLog.h"
  35 #include "TCanvas.h"
  36 #include "TF1.h"
  37 #include "AliTHn.h"
  38 #include "THn.h"
  39
  40 ClassImp(AliUEHist)
  41
  42 const Int_t AliUEHist::fgkCFSteps = 11;
  43
  44 AliUEHist::AliUEHist(const char* reqHist, const char* binning) :
  45   TObject(),
  46   fkRegions(4),
  47   fEventHist(0),
  48   fTrackHistEfficiency(0),
  49   fFakePt(0),
  50   fEtaMin(0),
  51   fEtaMax(0),
  52   fPtMin(0),
  53   fPtMax(0),
  54   fPartSpecies(-1),
  55   fCentralityMin(0),
  56   fCentralityMax(0),
  57   fZVtxMin(0),
  58   fZVtxMax(0),
  59   fContaminationEnhancement(0),
  60   fCombineMinMax(0),
  61   fTrackEtaCut(0),
  62   fWeightPerEvent(0),
  63   fCache(0),
  64   fGetMultCacheOn(kFALSE),
  65   fGetMultCache(0),
  66   fHistogramType(reqHist)
  67 {
  68   // Constructor
  69
  70   for (UInt_t i=0; i<fkRegions; i++)
  71     fTrackHist[i] = 0;
  72
  73   if (strlen(reqHist) == 0)
  74     return;
  75
  76   Printf("Creating AliUEHist with %s (binning: %s)", reqHist, binning);
  77
  78   AliLog::SetClassDebugLevel("AliCFContainer", -1);
  79   AliLog::SetClassDebugLevel("AliCFGridSparse", -3);
  80
  81   const char* title = "";
  82
  83   // track level
  84   Int_t nTrackVars = 4; // eta vs pT vs pT,lead (vs delta phi) vs multiplicity
  85   Int_t iTrackBin[6];
  86   Double_t* trackBins[6];
  87   const char* trackAxisTitle[6];
  88
  89   // eta
  90   Int_t nEtaBins = -1;
  91   Double_t* etaBins = GetBinning(binning, "eta", nEtaBins);
  92   const char* etaTitle = "#eta";
  93
  94   iTrackBin[0] = nEtaBins;
  95   trackBins[0] = etaBins;
  96   trackAxisTitle[0] = etaTitle;
  97
  98   // delta eta
  99   Int_t nDeltaEtaBins = -1;
 100   Double_t* deltaEtaBins = GetBinning(binning, "delta_eta", nDeltaEtaBins);
 101
 102   // pT
 103   trackBins[1] = GetBinning(binning, "p_t_assoc", iTrackBin[1]);
 104   trackAxisTitle[1] = "p_{T} (GeV/c)";
 105
 106   // pT, fine
 107   Int_t npTBinsFine = -1;
 108   Double_t* pTBinsFine = GetBinning(binning, "p_t_eff", npTBinsFine);
 109
 110   // pT,lead binning 1
 111   Int_t nLeadingpTBins = -1;
 112   Double_t* leadingpTBins = GetBinning(binning, "p_t_leading", nLeadingpTBins);
 113
 114   // pT,lead binning 2
 115   Int_t nLeadingpTBins2 = -1;
 116   Double_t* leadingpTBins2 = GetBinning(binning, "p_t_leading_course", nLeadingpTBins2);
 117
 118   // phi,lead
 119   Int_t nLeadingPhiBins = -1;
 120   Double_t* leadingPhiBins = GetBinning(binning, "delta_phi", nLeadingPhiBins);
 121
 122   trackBins[3] = GetBinning(binning, "multiplicity", iTrackBin[3]);
 123   trackAxisTitle[3] = "multiplicity";
 124
 125   // particle species
 126   const Int_t kNSpeciesBins = 4; // pi, K, p, rest
 127   Double_t speciesBins[] = { -0.5, 0.5, 1.5, 2.5, 3.5 };
 128
 129   // vtx-z axis
 130   const char* vertexTitle = "z-vtx (cm)";
 131   Int_t nVertexBins = -1;
 132   Double_t* vertexBins = GetBinning(binning, "vertex", nVertexBins);
 133   Int_t nVertexBinsEff = -1;
 134   Double_t* vertexBinsEff = GetBinning(binning, "vertex_eff", nVertexBinsEff);
 135
 136   Int_t useVtxAxis = 0;
 137
 138   // selection depending on requested histogram
 139   Int_t axis = -1; // 0 = pT,lead, 1 = phi,lead
 140   if (strcmp(reqHist, "NumberDensitypT") == 0)
 141   {
 142     axis = 0;
 143     title = "d^{2}N_{ch}/d#varphid#eta";
 144   }
 145   else if (strcmp(reqHist, "NumberDensityPhi") == 0)
 146   {
 147     axis = 1;
 148     title = "d^{2}N_{ch}/d#varphid#eta";
 149   }
 150   else if (TString(reqHist).BeginsWith("NumberDensityPhiCentrality"))
 151   {
 152     if (TString(reqHist).Contains("Vtx"))
 153       useVtxAxis = 1;
 154
 155     reqHist = "NumberDensityPhiCentrality";
 156     fHistogramType = reqHist;
 157     axis = 2;
 158     title = "d^{2}N_{ch}/d#varphid#eta";
 159   }
 160   else if (strcmp(reqHist, "SumpT") == 0)
 161   {
 162     axis = 0;
 163     title = "d^{2}#Sigma p_{T}/d#varphid#eta";
 164   }
 165   else
 166     AliFatal(Form("Invalid histogram requested: %s", reqHist));
 167
 168   UInt_t initRegions = fkRegions;
 169
 170   if (axis == 0)
 171   {
 172     trackBins[2] = leadingpTBins;
 173     iTrackBin[2] = nLeadingpTBins;
 174     trackAxisTitle[2] = "leading p_{T} (GeV/c)";
 175
 176   }
 177   else if (axis == 1)
 178   {
 179     nTrackVars = 5;
 180     initRegions = 1;
 181
 182     iTrackBin[2] = nLeadingpTBins2;
 183     trackBins[2] = leadingpTBins2;
 184     trackAxisTitle[2] = "leading p_{T} (GeV/c)";
 185
 186     iTrackBin[4] = nLeadingPhiBins;
 187     trackBins[4] = leadingPhiBins;
 188     trackAxisTitle[4] = "#Delta #varphi w.r.t. leading track";
 189   }
 190   else if (axis == 2)
 191   {
 192     nTrackVars = 5;
 193     initRegions = 1;
 194
 195     iTrackBin[0] = nDeltaEtaBins;
 196     trackBins[0] = deltaEtaBins;
 197     trackAxisTitle[0] = "#Delta#eta";
 198
 199     iTrackBin[2] = nLeadingpTBins2;
 200     trackBins[2] = leadingpTBins2;
 201     trackAxisTitle[2] = "leading p_{T} (GeV/c)";
 202
 203     trackAxisTitle[3] = "centrality";
 204
 205     iTrackBin[4] = nLeadingPhiBins;
 206     trackBins[4] = leadingPhiBins;
 207     trackAxisTitle[4] = "#Delta#varphi (rad)";
 208
 209     if (useVtxAxis > 0)
 210     {
 211       nTrackVars = 6;
 212       iTrackBin[5] = nVertexBins;
 213       trackBins[5] = vertexBins;
 214       trackAxisTitle[5] = vertexTitle;
 215     }
 216   }
 217
 218   for (UInt_t i=0; i<initRegions; i++)
 219   {
 220     if (strcmp(reqHist, "NumberDensityPhiCentrality") == 0)
 221       fTrackHist[i] = new AliTHn(Form("fTrackHist_%d", i), title, fgkCFSteps, nTrackVars, iTrackBin);
 222     else
 223       fTrackHist[i] = new AliCFContainer(Form("fTrackHist_%d", i), title, fgkCFSteps, nTrackVars, iTrackBin);
 224
 225     for (Int_t j=0; j<nTrackVars; j++)
 226     {
 227       fTrackHist[i]->SetBinLimits(j, trackBins[j]);
 228       fTrackHist[i]->SetVarTitle(j, trackAxisTitle[j]);
 229     }
 230
 231     SetStepNames(fTrackHist[i]);
 232   }
 233
 234   // event level
 235   Int_t nEventVars = 2;
 236   Int_t iEventBin[3];
 237
 238   // track 3rd and 4th axis --> event 1st and 2nd axis
 239   iEventBin[0] = iTrackBin[2];
 240   iEventBin[1] = iTrackBin[3];
 241
 242   // plus track 5th axis (in certain cases)
 243   if (axis == 2 && useVtxAxis)
 244   {
 245     nEventVars = 3;
 246     iEventBin[2] = iTrackBin[5];
 247   }
 248
 249   fEventHist = new AliCFContainer("fEventHist", title, fgkCFSteps, nEventVars, iEventBin);
 250
 251   fEventHist->SetBinLimits(0, trackBins[2]);
 252   fEventHist->SetVarTitle(0, trackAxisTitle[2]);
 253
 254   fEventHist->SetBinLimits(1, trackBins[3]);
 255   fEventHist->SetVarTitle(1, trackAxisTitle[3]);
 256
 257   if (axis == 2 && useVtxAxis)
 258   {
 259     fEventHist->SetBinLimits(2, trackBins[5]);
 260     fEventHist->SetVarTitle(2, trackAxisTitle[5]);
 261   }
 262
 263   SetStepNames(fEventHist);
 264
 265   iTrackBin[0] = nEtaBins;
 266   iTrackBin[1] = npTBinsFine;
 267   iTrackBin[2] = kNSpeciesBins;
 268   iTrackBin[4] = nVertexBinsEff;
 269
 270   fTrackHistEfficiency = new AliCFContainer("fTrackHistEfficiency", "Tracking efficiency", 6, 5, iTrackBin);
 271   fTrackHistEfficiency->SetBinLimits(0, etaBins);
 272   fTrackHistEfficiency->SetVarTitle(0, etaTitle);
 273   fTrackHistEfficiency->SetBinLimits(1, pTBinsFine);
 274   fTrackHistEfficiency->SetVarTitle(1, trackAxisTitle[1]);
 275   fTrackHistEfficiency->SetBinLimits(2, speciesBins);
 276   fTrackHistEfficiency->SetVarTitle(2, "particle species");
 277   fTrackHistEfficiency->SetBinLimits(3, trackBins[3]);
 278   fTrackHistEfficiency->SetVarTitle(3, trackAxisTitle[3]);
 279   fTrackHistEfficiency->SetBinLimits(4, vertexBinsEff);
 280   fTrackHistEfficiency->SetVarTitle(4, vertexTitle);
 281
 282   fFakePt = new TH3F("fFakePt","fFakePt;p_{T,rec};p_{T};centrality", 200, 0, 20, 200, 0, 20, 20, 0, 100);
 283
 284   delete[] deltaEtaBins;
 285   delete[] pTBinsFine;
 286   delete[] leadingpTBins;
 287   delete[] leadingpTBins2;
 288   delete[] leadingPhiBins;
 289   delete[] vertexBins;
 290   delete[] vertexBinsEff;
 291 }
 292
 293 Double_t* AliUEHist::GetBinning(const char* configuration, const char* tag, Int_t& nBins)
 294 {
 295   // takes the binning from <configuration> identified by <tag>
 296   // configuration syntax example:
 297   // eta: 2.4, -2.3, -2.2, -2.1, -2.0, -1.9, -1.8, -1.7, -1.6, -1.5, -1.4, -1.3, -1.2, -1.1, -1.0, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.1, 2.2, 2.3, 2.4
 298   // phi: .....
 299   //
 300   // returns bin edges which have to be deleted by the caller
 301
 302   TString config(configuration);
 303   TObjArray* lines = config.Tokenize("\n");
 304   for (Int_t i=0; i<lines->GetEntriesFast(); i++)
 305   {
 306     TString line(lines->At(i)->GetName());
 307     if (line.BeginsWith(TString(tag) + ":"))
 308     {
 309       line.Remove(0, strlen(tag) + 1);
 310       line.ReplaceAll(" ", "");
 311       TObjArray* binning = line.Tokenize(",");
 312       Double_t* bins = new Double_t[binning->GetEntriesFast()];
 313       for (Int_t j=0; j<binning->GetEntriesFast(); j++)
 314         bins[j] = TString(binning->At(j)->GetName()).Atof();
 315
 316       nBins = binning->GetEntriesFast() - 1;
 317
 318       delete binning;
 319       delete lines;
 320       return bins;
 321     }
 322   }
 323
 324   delete lines;
 325   AliFatal(Form("Tag %s not found in %s", tag, configuration));
 326   return 0;
 327 }
 328
 329 //_____________________________________________________________________________
 330 AliUEHist::AliUEHist(const AliUEHist &c) :
 331   TObject(),
 332   fkRegions(4),
 333   fEventHist(0),
 334   fTrackHistEfficiency(0),
 335   fFakePt(0),
 336   fEtaMin(0),
 337   fEtaMax(0),
 338   fPtMin(0),
 339   fPtMax(0),
 340   fPartSpecies(-1),
 341   fCentralityMin(0),
 342   fCentralityMax(0),
 343   fZVtxMin(0),
 344   fZVtxMax(0),
 345   fContaminationEnhancement(0),
 346   fCombineMinMax(0),
 347   fTrackEtaCut(0),
 348   fWeightPerEvent(0),
 349   fCache(0),
 350   fGetMultCacheOn(kFALSE),
 351   fGetMultCache(0),
 352   fHistogramType()
 353 {
 354   //
 355   // AliUEHist copy constructor
 356   //
 357
 358   ((AliUEHist &) c).Copy(*this);
 359 }
 360
 361 //____________________________________________________________________
 362 void AliUEHist::SetStepNames(AliCFContainer* container)
 363 {
 364   // sets the names of the correction steps
 365
 366   for (Int_t i=0; i<fgkCFSteps; i++)
 367     container->SetStepTitle(i, GetStepTitle((CFStep) i));
 368 }
 369
 370 //____________________________________________________________________
 371 AliUEHist::~AliUEHist()
 372 {
 373   // Destructor
 374
 375   for (UInt_t i=0; i<fkRegions; i++)
 376   {
 377     if (fTrackHist[i])
 378     {
 379       delete fTrackHist[i];
 380       fTrackHist[i] = 0;
 381     }
 382   }
 383
 384   if (fEventHist)
 385   {
 386     delete fEventHist;
 387     fEventHist = 0;
 388   }
 389
 390   if (fTrackHistEfficiency)
 391   {
 392     delete fTrackHistEfficiency;
 393     fTrackHistEfficiency = 0;
 394   }
 395
 396   if (fFakePt)
 397   {
 398     delete fFakePt;
 399     fFakePt = 0;
 400   }
 401
 402   if (fCache)
 403   {
 404     delete fCache;
 405     fCache = 0;
 406   }
 407 }
 408
 409 //____________________________________________________________________
 410 AliUEHist &AliUEHist::operator=(const AliUEHist &c)
 411 {
 412   // assigment operator
 413
 414   if (this != &c)
 415     ((AliUEHist &) c).Copy(*this);
 416
 417   return *this;
 418 }
 419
 420 //____________________________________________________________________
 421 void AliUEHist::Copy(TObject& c) const
 422 {
 423   // copy function
 424
 425   AliUEHist& target = (AliUEHist &) c;
 426
 427   for (UInt_t i=0; i<fkRegions; i++)
 428     if (fTrackHist[i])
 429       target.fTrackHist[i] = dynamic_cast<AliCFContainer*> (fTrackHist[i]->Clone());
 430
 431   if (fEventHist)
 432     target.fEventHist = dynamic_cast<AliCFContainer*> (fEventHist->Clone());
 433
 434   if (fTrackHistEfficiency)
 435     target.fTrackHistEfficiency = dynamic_cast<AliCFContainer*> (fTrackHistEfficiency->Clone());
 436
 437   if (fFakePt)
 438     target.fFakePt = dynamic_cast<TH3F*> (fFakePt->Clone());
 439
 440   target.fEtaMin = fEtaMin;
 441   target.fEtaMax = fEtaMax;
 442   target.fPtMin = fPtMin;
 443   target.fPtMax = fPtMax;
 444   target.fPartSpecies = fPartSpecies;
 445   target.fCentralityMin = fCentralityMin;
 446   target.fCentralityMax = fCentralityMax;
 447   target.fZVtxMin = fZVtxMin;
 448   target.fZVtxMax = fZVtxMax;
 449
 450   if (fContaminationEnhancement)
 451     target.fContaminationEnhancement = dynamic_cast<TH1F*> (fContaminationEnhancement->Clone());
 452
 453   target.fCombineMinMax = fCombineMinMax;
 454   target.fTrackEtaCut = fTrackEtaCut;
 455   target.fWeightPerEvent = fWeightPerEvent;
 456   target.fHistogramType = fHistogramType;
 457 }
 458
 459 //____________________________________________________________________
 460 Long64_t AliUEHist::Merge(TCollection* list)
 461 {
 462   // Merge a list of AliUEHist objects with this (needed for
 463   // PROOF).
 464   // Returns the number of merged objects (including this).
 465
 466   if (!list)
 467     return 0;
 468
 469   if (list->IsEmpty())
 470     return 1;
 471
 472   TIterator* iter = list->MakeIterator();
 473   TObject* obj;
 474
 475   // collections of objects
 476   const UInt_t kMaxLists = fkRegions+3;
 477   TList** lists = new TList*[kMaxLists];
 478
 479   for (UInt_t i=0; i<kMaxLists; i++)
 480     lists[i] = new TList;
 481
 482   Int_t count = 0;
 483   while ((obj = iter->Next())) {
 484
 485     AliUEHist* entry = dynamic_cast<AliUEHist*> (obj);
 486     if (entry == 0)
 487       continue;
 488
 489     for (UInt_t i=0; i<fkRegions; i++)
 490       if (entry->fTrackHist[i])
 491         lists[i]->Add(entry->fTrackHist[i]);
 492
 493     lists[fkRegions]->Add(entry->fEventHist);
 494     lists[fkRegions+1]->Add(entry->fTrackHistEfficiency);
 495     if (entry->fFakePt)
 496       lists[fkRegions+2]->Add(entry->fFakePt);
 497
 498     count++;
 499   }
 500   for (UInt_t i=0; i<fkRegions; i++)
 501     if (fTrackHist[i])
 502       fTrackHist[i]->Merge(lists[i]);
 503
 504   fEventHist->Merge(lists[fkRegions]);
 505   fTrackHistEfficiency->Merge(lists[fkRegions+1]);
 506   if (fFakePt)
 507     fFakePt->Merge(lists[fkRegions+2]);
 508
 509   for (UInt_t i=0; i<kMaxLists; i++)
 510     delete lists[i];
 511
 512   delete[] lists;
 513
 514   return count+1;
 515 }
 516
 517 //____________________________________________________________________
 518 void AliUEHist::SetBinLimits(THnBase* grid)
 519 {
 520   // sets the bin limits in eta and pT defined by fEtaMin/Max, fPtMin/Max
 521
 522   if (fEtaMax > fEtaMin)
 523     grid->GetAxis(0)->SetRangeUser(fEtaMin, fEtaMax);
 524   if (fPtMax > fPtMin)
 525     grid->GetAxis(1)->SetRangeUser(fPtMin, fPtMax);
 526 }
 527
 528 //____________________________________________________________________
 529 void AliUEHist::SetBinLimits(AliCFGridSparse* grid)
 530 {
 531   // sets the bin limits in eta and pT defined by fEtaMin/Max, fPtMin/Max
 532
 533   SetBinLimits(grid->GetGrid());
 534 }
 535
 536 //____________________________________________________________________
 537 void AliUEHist::ResetBinLimits(THnBase* grid)
 538 {
 539   // resets all bin limits
 540
 541   for (Int_t i=0; i<grid->GetNdimensions(); i++)
 542     if (grid->GetAxis(i)->TestBit(TAxis::kAxisRange))
 543       grid->GetAxis(i)->SetRangeUser(0, -1);
 544 }
 545
 546 //____________________________________________________________________
 547 void AliUEHist::ResetBinLimits(AliCFGridSparse* grid)
 548 {
 549   // resets all bin limits
 550
 551   ResetBinLimits(grid->GetGrid());
 552 }
 553
 554 //____________________________________________________________________
 555 void AliUEHist::CountEmptyBins(AliUEHist::CFStep step, Float_t ptLeadMin, Float_t ptLeadMax)
 556 {
 557   // prints the number of empty bins in the track end event histograms in the given step
 558
 559   Int_t binBegin[4];
 560   Int_t binEnd[4];
 561
 562   for (Int_t i=0; i<4; i++)
 563   {
 564     binBegin[i] = 1;
 565     binEnd[i]   = fTrackHist[0]->GetNBins(i);
 566   }
 567
 568   if (fEtaMax > fEtaMin)
 569   {
 570     binBegin[0] = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(0)->FindBin(fEtaMin);
 571     binEnd[0]   = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(0)->FindBin(fEtaMax);
 572   }
 573
 574   if (fPtMax > fPtMin)
 575   {
 576     binBegin[1] = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(1)->FindBin(fPtMin);
 577     binEnd[1]   = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(1)->FindBin(fPtMax);
 578   }
 579
 580   if (ptLeadMax > ptLeadMin)
 581   {
 582     binBegin[2] = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(2)->FindBin(ptLeadMin);
 583     binEnd[2]   = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(2)->FindBin(ptLeadMax);
 584   }
 585
 586   // start from multiplicity 1
 587   binBegin[3] = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(3)->FindBin(1);
 588
 589   for (UInt_t region=0; region<fkRegions; region++)
 590   {
 591     Int_t total = 0;
 592     Int_t count = 0;
 593     Int_t vars[4];
 594
 595     for (Int_t i=0; i<4; i++)
 596       vars[i] = binBegin[i];
 597
 598     AliCFGridSparse* grid = fTrackHist[region]->GetGrid(step);
 599     while (1)
 600     {
 601       if (grid->GetElement(vars) == 0)
 602       {
 603         Printf("Empty bin at eta=%.2f pt=%.2f pt_lead=%.2f mult=%.1f",
 604           grid->GetBinCenter(0, vars[0]),
 605           grid->GetBinCenter(1, vars[1]),
 606           grid->GetBinCenter(2, vars[2]),
 607           grid->GetBinCenter(3, vars[3])
 608         );
 609         count++;
 610       }
 611
 612       vars[3]++;
 613       for (Int_t i=3; i>0; i--)
 614       {
 615         if (vars[i] == binEnd[i]+1)
 616         {
 617           vars[i] = binBegin[i];
 618           vars[i-1]++;
 619         }
 620       }
 621
 622       if (vars[0] == binEnd[0]+1)
 623         break;
 624       total++;
 625     }
 626
 627     Printf("Region %s has %d empty bins (out of %d bins)", GetRegionTitle((Region) region), count, total);
 628   }
 629 }
 630
 631 //____________________________________________________________________
 632 TH1* AliUEHist::GetUEHist(AliUEHist::CFStep step, AliUEHist::Region region, Float_t ptLeadMin, Float_t ptLeadMax, Int_t multBinBegin, Int_t multBinEnd, Int_t twoD, Bool_t etaNorm, Long64_t* normEvents)
 633 {
 634   // Extracts the UE histogram at the given step and in the given region by projection and dividing tracks by events
 635   //
 636   // ptLeadMin, ptLeadMax: Only meaningful for vs. delta phi plot (third axis is ptLead)
 637   // Histogram has to be deleted by the caller of the function
 638   //
 639   // twoD: 0: 1D histogram as function of phi
 640   //       1: 2D histogram, deltaphi, deltaeta
 641   //       10: 1D histogram, within |deltaeta| < 0.8
 642   //       11: 1D histogram, within 0.8 < |deltaeta| < 1.6
 643   //
 644   // etaNorm: if kTRUE (default), the distributions are divided by the area in delta eta
 645   //
 646   // normEvents: if non-0 the number of events/trigger particles for the normalization is filled
 647
 648   // unzoom all axes
 649   ResetBinLimits(fTrackHist[region]->GetGrid(step));
 650   ResetBinLimits(fEventHist->GetGrid(step));
 651
 652   SetBinLimits(fTrackHist[region]->GetGrid(step));
 653
 654   // z vtx
 655   if (fZVtxMax > fZVtxMin)
 656   {
 657     Printf("Using z-vtx range %f --> %f", fZVtxMin, fZVtxMax);
 658     fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(5)->SetRangeUser(fZVtxMin, fZVtxMax);
 659     fEventHist->GetGrid(step)->GetGrid()->GetAxis(2)->SetRangeUser(fZVtxMin, fZVtxMax);
 660   }
 661
 662   TH1D* tracks = 0;
 663
 664   if (ptLeadMin < 0)
 665   {
 666     tracks = fTrackHist[region]->ShowProjection(2, step);
 667     tracks->GetYaxis()->SetTitle(fTrackHist[region]->GetTitle());
 668     if (fCombineMinMax && region == kMin)
 669     {
 670       ResetBinLimits(fTrackHist[kMax]->GetGrid(step));
 671       SetBinLimits(fTrackHist[kMax]->GetGrid(step));
 672
 673       TH1D* tracks2 = fTrackHist[kMax]->ShowProjection(2, step);
 674       tracks->Add(tracks2);
 675
 676       ResetBinLimits(fTrackHist[kMax]->GetGrid(step));
 677     }
 678
 679     // normalize to get a density (deta dphi)
 680     TAxis* axis = fTrackHist[region]->GetGrid(step)->GetAxis(0);
 681     Float_t phiRegion = TMath::TwoPi() / 3;
 682     if (!fCombineMinMax && region == kMin)
 683       phiRegion /= 2;
 684     Float_t normalization = phiRegion;
 685     if (etaNorm)
 686       normalization *= (axis->GetBinUpEdge(axis->GetLast()) - axis->GetBinLowEdge(axis->GetFirst()));
 687     //Printf("Normalization: %f", normalization);
 688     tracks->Scale(1.0 / normalization);
 689
 690     TH1D* events = fEventHist->ShowProjection(0, step);
 691     tracks->Divide(events);
 692   }
 693   else
 694   {
 695     if (multBinEnd >= multBinBegin)
 696     {
 697       Printf("Using multiplicity range %d --> %d", multBinBegin, multBinEnd);
 698       fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(3)->SetRange(multBinBegin, multBinEnd);
 699       fEventHist->GetGrid(step)->GetGrid()->GetAxis(1)->SetRange(multBinBegin, multBinEnd);
 700     }
 701
 702     Int_t firstBin = fTrackHist[region]->GetAxis(2, step)->FindBin(ptLeadMin);
 703     Int_t lastBin = fTrackHist[region]->GetAxis(2, step)->FindBin(ptLeadMax);
 704
 705     Printf("Using leading pT range %d --> %d", firstBin, lastBin);
 706
 707     fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(2)->SetRange(firstBin, lastBin);
 708
 709     if (twoD == 0)
 710       tracks = (TH1D*) fTrackHist[region]->GetGrid(step)->Project(4);
 711     else
 712       tracks = (TH1D*) fTrackHist[region]->GetGrid(step)->Project(4, 0);
 713
 714     Printf("Calculated histogram --> %f tracks", tracks->Integral());
 715     fTrackHist[region]->GetGrid(step)->SetRangeUser(2, 0, -1);
 716
 717     if (twoD == 10 || twoD == 11)
 718     {
 719       Float_t etaLimit = 1.0;
 720       if (twoD == 10)
 721       {
 722         tracks = (TH1D*) ((TH2*) tracks)->ProjectionX("proj", tracks->GetYaxis()->FindBin(-etaLimit + 0.01), tracks->GetYaxis()->FindBin(etaLimit - 0.01))->Clone();
 723
 724         // TODO calculate acc with 2 * (deta - 0.5 * deta*deta / 1.6)
 725         tracks->Scale(1. / 0.75);
 726       }
 727       else if (twoD == 11)
 728       {
 729         TH1D* tracksTmp1 = (TH1D*) ((TH2*) tracks)->ProjectionX("proj1", tracks->GetYaxis()->FindBin(-etaLimit * 2 + 0.01), tracks->GetYaxis()->FindBin(-etaLimit - 0.01))->Clone();
 730         TH1D* tracksTmp2 = ((TH2*) tracks)->ProjectionX("proj2", tracks->GetYaxis()->FindBin(etaLimit + 0.01), tracks->GetYaxis()->FindBin(2 * etaLimit - 0.01));
 731
 732         tracksTmp1->Add(tracksTmp2);
 733
 734         delete tracks;
 735         tracks = tracksTmp1;
 736         delete tracksTmp2;
 737
 738         tracks->Scale(1. / 0.25);
 739       }
 740     }
 741
 742     // normalize to get a density (deta dphi)
 743     // TODO normalization may be off for 2d histogram
 744     Float_t normalization = fTrackHist[region]->GetGrid(step)->GetAxis(4)->GetBinWidth(1);
 745
 746     if (etaNorm)
 747     {
 748       TAxis* axis = fTrackHist[region]->GetGrid(step)->GetAxis(0);
 749       if (strcmp(axis->GetTitle(), "#eta") == 0)
 750       {
 751         Printf("Normalizing using eta-axis range");
 752         normalization *= axis->GetBinUpEdge(axis->GetLast()) - axis->GetBinLowEdge(axis->GetFirst());
 753       }
 754       else
 755       {
 756         Printf("Normalizing assuming accepted range of |eta| < 0.8");
 757         normalization *= 0.8 * 2;
 758       }
 759     }
 760
 761     //Printf("Normalization: %f", normalization);
 762     tracks->Scale(1.0 / normalization);
 763
 764     // NOTE fEventHist contains the number of events for the underlying event analysis and the number of trigger particles for the azimuthal correlation analysis. In the latter case the naming is therefore somewhat misleading!
 765     TH1D* events = fEventHist->ShowProjection(0, step);
 766     Long64_t nEvents = (Long64_t) events->Integral(firstBin, lastBin);
 767     Printf("Calculated histogram --> %lld events", nEvents);
 768     if (normEvents)
 769       *normEvents = nEvents;
 770
 771     if (nEvents > 0)
 772       tracks->Scale(1.0 / nEvents);
 773
 774     delete events;
 775   }
 776
 777   ResetBinLimits(fTrackHist[region]->GetGrid(step));
 778   ResetBinLimits(fEventHist->GetGrid(step));
 779
 780   return tracks;
 781 }
 782
 783 //____________________________________________________________________
 784 void AliUEHist::GetHistsZVtx(AliUEHist::CFStep step, AliUEHist::Region region, Float_t ptLeadMin, Float_t ptLeadMax, Int_t multBinBegin, Int_t multBinEnd, TH3** trackHist, TH1** eventHist)
 785 {
 786   // Calculates a 3d histogram with deltaphi, deltaeta, zvtx on track level and a 1d histogram on event level (as fct of zvtx)
 787   // Histogram has to be deleted by the caller of the function
 788
 789   if (!fTrackHist[kToward]->GetGrid(0)->GetGrid()->GetAxis(5))
 790     AliFatal("Histogram without vertex axis provided");
 791
 792   // unzoom all axes
 793   ResetBinLimits(fTrackHist[region]->GetGrid(step));
 794   ResetBinLimits(fEventHist->GetGrid(step));
 795
 796   SetBinLimits(fTrackHist[region]->GetGrid(step));
 797
 798   if (multBinEnd >= multBinBegin)
 799   {
 800     Printf("Using multiplicity range %d --> %d", multBinBegin, multBinEnd);
 801     fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(3)->SetRange(multBinBegin, multBinEnd);
 802     fEventHist->GetGrid(step)->GetGrid()->GetAxis(1)->SetRange(multBinBegin, multBinEnd);
 803   }
 804
 805   Int_t firstBin = fTrackHist[region]->GetAxis(2, step)->FindBin(ptLeadMin);
 806   Int_t lastBin = fTrackHist[region]->GetAxis(2, step)->FindBin(ptLeadMax);
 807   Printf("Using leading pT range %d --> %d", firstBin, lastBin);
 808   fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(2)->SetRange(firstBin, lastBin);
 809   fEventHist->GetGrid(step)->GetGrid()->GetAxis(0)->SetRange(firstBin, lastBin);
 810
 811   *trackHist = (TH3*) fTrackHist[region]->GetGrid(step)->Project(4, 0, 5);
 812   *eventHist = (TH1*) fEventHist->GetGrid(step)->Project(2);
 813
 814   ResetBinLimits(fTrackHist[region]->GetGrid(step));
 815   ResetBinLimits(fEventHist->GetGrid(step));
 816 }
 817
 818 //____________________________________________________________________
 819 void AliUEHist::GetHistsZVtxMult(AliUEHist::CFStep step, AliUEHist::Region region, Float_t ptLeadMin, Float_t ptLeadMax, THnBase** trackHist, TH2** eventHist)
 820 {
 821   // Calculates a 4d histogram with deltaphi, deltaeta, zvtx, multiplicity on track level and
 822   // a 2d histogram on event level (as fct of zvtx, multiplicity)
 823   // Histograms has to be deleted by the caller of the function
 824
 825   if (!fTrackHist[kToward]->GetGrid(0)->GetGrid()->GetAxis(5))
 826     AliFatal("Histogram without vertex axis provided");
 827
 828   THnBase* sparse = fTrackHist[region]->GetGrid(step)->GetGrid();
 829   if (fGetMultCacheOn)
 830   {
 831     if (!fGetMultCache)
 832     {
 833       fGetMultCache = ChangeToThn(sparse);
 834       // should work but causes SEGV in ProjectionND below
 835     }
 836     sparse = fGetMultCache;
 837   }
 838
 839   // unzoom all axes
 840   ResetBinLimits(sparse);
 841   ResetBinLimits(fEventHist->GetGrid(step));
 842
 843   SetBinLimits(sparse);
 844
 845   Int_t firstBin = sparse->GetAxis(2)->FindBin(ptLeadMin);
 846   Int_t lastBin = sparse->GetAxis(2)->FindBin(ptLeadMax);
 847   Printf("Using leading pT range %d --> %d", firstBin, lastBin);
 848   sparse->GetAxis(2)->SetRange(firstBin, lastBin);
 849   fEventHist->GetGrid(step)->GetGrid()->GetAxis(0)->SetRange(firstBin, lastBin);
 850
 851   Int_t dimensions[] = { 4, 0, 5, 3 };
 852   THnBase* tmpTrackHist = sparse->ProjectionND(4, dimensions, "E");
 853   *eventHist = (TH2*) fEventHist->GetGrid(step)->Project(2, 1);
 854
 855   ResetBinLimits(sparse);
 856   ResetBinLimits(fEventHist->GetGrid(step));
 857
 858   // convert to THn
 859   *trackHist = ChangeToThn(tmpTrackHist);
 860   delete tmpTrackHist;
 861 }
 862
 863 //____________________________________________________________________
 864 TH2* AliUEHist::GetSumOfRatios2(AliUEHist* mixed, AliUEHist::CFStep step, AliUEHist::Region region, Float_t ptLeadMin, Float_t ptLeadMax, Int_t multBinBegin, Int_t multBinEnd, Bool_t normalizePerTrigger)
 865 {
 866   // Calls GetUEHist(...) for *each* vertex bin and multiplicity bin and performs a sum of ratios:
 867   // 1_N [ (same/mixed)_1 + (same/mixed)_2 + (same/mixed)_3 + ... ]
 868   // where N is the total number of events/trigger particles and the subscript is the vertex/multiplicity bin
 869   // where mixed is normalized such that the information about the number of pairs in same is kept
 870   //
 871   // returns a 2D histogram: deltaphi, deltaeta
 872   //
 873   // Parameters:
 874   //   mixed: AliUEHist containing mixed event corresponding to this object
 875   //   <other parameters> : check documentation of AliUEHist::GetUEHist
 876   //  normalizePerTrigger: divide through number of triggers
 877
 878   // do not add this hists to the directory
 879   Bool_t oldStatus = TH1::AddDirectoryStatus();
 880   TH1::AddDirectory(kFALSE);
 881
 882   TH2* totalTracks = 0;
 883
 884   THnBase* trackSameAll = 0;
 885   THnBase* trackMixedAll = 0;
 886   TH2* eventSameAll = 0;
 887   TH2* eventMixedAll = 0;
 888
 889   Int_t totalEvents = 0;
 890   Int_t nCorrelationFunctions = 0;
 891
 892   GetHistsZVtxMult(step, region, ptLeadMin, ptLeadMax, &trackSameAll, &eventSameAll);
 893   mixed->GetHistsZVtxMult(step, region, ptLeadMin, ptLeadMax, &trackMixedAll, &eventMixedAll);
 894
 895 //   Printf("%f %f %f %f", trackSameAll->GetEntries(), eventSameAll->GetEntries(), trackMixedAll->GetEntries(), eventMixedAll->GetEntries());
 896
 897 //   TH1* normParameters = new TH1F("normParameters", "", 100, 0, 2);
 898
 899 //   trackSameAll->Dump();
 900
 901   TAxis* multAxis = trackSameAll->GetAxis(3);
 902
 903   if (multBinEnd < multBinBegin)
 904   {
 905     multBinBegin = 1;
 906     multBinEnd = multAxis->GetNbins();
 907   }
 908
 909   for (Int_t multBin = TMath::Max(1, multBinBegin); multBin <= TMath::Min(multAxis->GetNbins(), multBinEnd); multBin++)
 910   {
 911     trackSameAll->GetAxis(3)->SetRange(multBin, multBin);
 912     trackMixedAll->GetAxis(3)->SetRange(multBin, multBin);
 913
 914     // get mixed normalization correction factor: is independent of vertex bin if scaled with number of triggers
 915     trackMixedAll->GetAxis(2)->SetRange(0, -1);
 916     TH2* tracksMixed = trackMixedAll->Projection(1, 0, "E");
 917 //     Printf("%f", tracksMixed->Integral());
 918     Float_t binWidthEta = tracksMixed->GetYaxis()->GetBinWidth(1);
 919
 920     // get mixed event normalization by assuming full acceptance at deta at 0 (integrate over dphi), excluding (0, 0)
 921     Double_t mixedNormError = 0;
 922     Double_t mixedNorm = tracksMixed->IntegralAndError(1, tracksMixed->GetYaxis()->FindBin(-0.01)-1, tracksMixed->GetYaxis()->FindBin(-0.01), tracksMixed->GetYaxis()->FindBin(0.01), mixedNormError);
 923     Double_t mixedNormError2 = 0;
 924     Double_t mixedNorm2 = tracksMixed->IntegralAndError(tracksMixed->GetYaxis()->FindBin(0.01)+1, tracksMixed->GetNbinsX(), tracksMixed->GetYaxis()->FindBin(-0.01), tracksMixed->GetYaxis()->FindBin(0.01), mixedNormError2);
 925
 926     if (mixedNormError == 0 || mixedNormError2 == 0)
 927     {
 928       Printf("ERROR: Skipping multiplicity %d because mixed event is empty %f %f %f %f", multBin, mixedNorm, mixedNormError, mixedNorm2, mixedNormError2);
 929       continue;
 930     }
 931
 932     Int_t nBinsMixedNorm = (tracksMixed->GetYaxis()->FindBin(-0.01) - 1 - 1 + 1) * (tracksMixed->GetYaxis()->FindBin(0.01) - tracksMixed->GetYaxis()->FindBin(-0.01) + 1);
 933     mixedNorm /= nBinsMixedNorm;
 934     mixedNormError /= nBinsMixedNorm;
 935
 936     Int_t nBinsMixedNorm2 = (tracksMixed->GetNbinsX() - tracksMixed->GetYaxis()->FindBin(0.01) - 1 + 1) * (tracksMixed->GetYaxis()->FindBin(0.01) - tracksMixed->GetYaxis()->FindBin(-0.01) + 1);
 937     mixedNorm2 /= nBinsMixedNorm2;
 938     mixedNormError2 /= nBinsMixedNorm2;
 939
 940     mixedNorm = mixedNorm / mixedNormError / mixedNormError + mixedNorm2 / mixedNormError2 / mixedNormError2;
 941     mixedNormError = TMath::Sqrt(1.0 / (1.0 / mixedNormError / mixedNormError + 1.0 / mixedNormError2 / mixedNormError2));
 942     mixedNorm *= mixedNormError * mixedNormError;
 943
 944 /*    Double_t mixedNorm = tracksMixed->IntegralAndError(tracksMixed->GetXaxis()->FindBin(-0.01), tracksMixed->GetXaxis()->FindBin(0.01), tracksMixed->GetYaxis()->FindBin(-0.01), tracksMixed->GetYaxis()->FindBin(0.01), mixedNormError);
 945     Int_t nBinsMixedNorm = (tracksMixed->GetXaxis()->FindBin(0.01) - tracksMixed->GetXaxis()->FindBin(-0.01) + 1) * (tracksMixed->GetYaxis()->FindBin(0.01) - tracksMixed->GetYaxis()->FindBin(-0.01) + 1);*/
 946
 947     delete tracksMixed;
 948
 949     Float_t triggers = eventMixedAll->Integral(1, eventMixedAll->GetNbinsX(), multBin, multBin);
 950 //     Printf("%f +- %f | %f | %f", mixedNorm, mixedNormError, triggers, mixedNorm / triggers);
 951     if (triggers <= 0)
 952     {
 953       Printf("ERROR: Skipping multiplicity %d because mixed event is empty", multBin);
 954       continue;
 955     }
 956
 957     mixedNorm /= triggers;
 958     mixedNormError /= triggers;
 959
 960     if (mixedNorm <= 0)
 961     {
 962       Printf("ERROR: Skipping multiplicity %d because mixed event is empty at (0,0)", multBin);
 963       continue;
 964     }
 965
 966     // finite bin correction
 967     if (fTrackEtaCut > 0)
 968     {
 969       Double_t finiteBinCorrection = -1.0 / (2*fTrackEtaCut) * binWidthEta / 2 + 1;
 970       Printf("Finite bin correction: %f", finiteBinCorrection);
 971       mixedNorm *= finiteBinCorrection;
 972       mixedNormError *= finiteBinCorrection;
 973     }
 974     else
 975     {
 976       Printf("ERROR: fTrackEtaCut not set. Finite bin correction cannot be applied. Continuing anyway...");
 977     }
 978
 979 //     Printf("Norm: %f +- %f", mixedNorm, mixedNormError);
 980
 981 //     normParameters->Fill(mixedNorm);
 982
 983     TAxis* vertexAxis = trackSameAll->GetAxis(2);
 984     for (Int_t vertexBin = 1; vertexBin <= vertexAxis->GetNbins(); vertexBin++)
 985 //     for (Int_t vertexBin = 5; vertexBin <= 6; vertexBin++)
 986     {
 987       trackSameAll->GetAxis(2)->SetRange(vertexBin, vertexBin);
 988       trackMixedAll->GetAxis(2)->SetRange(vertexBin, vertexBin);
 989
 990       TH2* tracksSame = trackSameAll->Projection(1, 0, "E");
 991       tracksMixed = trackMixedAll->Projection(1, 0, "E");
 992
 993       // asssume flat in dphi, gain in statistics
 994       //     TH1* histMixedproj = mixedTwoD->ProjectionY();
 995       //     histMixedproj->Scale(1.0 / mixedTwoD->GetNbinsX());
 996       //
 997       //     for (Int_t x=1; x<=mixedTwoD->GetNbinsX(); x++)
 998       //       for (Int_t y=1; y<=mixedTwoD->GetNbinsY(); y++)
 999       //        mixedTwoD->SetBinContent(x, y, histMixedproj->GetBinContent(y));
1000
1001       //       delete histMixedproj;
1002
1003       Float_t triggers2 = eventMixedAll->Integral(vertexBin, vertexBin, multBin, multBin);
1004       if (triggers2 <= 0)
1005       {
1006         Printf("ERROR: Skipping multiplicity %d vertex bin %d because mixed event is empty", multBin, vertexBin);
1007       }
1008       else
1009       {
1010         tracksMixed->Scale(1.0 / triggers2 / mixedNorm);
1011         // tracksSame->Scale(tracksMixed->Integral() / tracksSame->Integral());
1012
1013 //      new TCanvas; tracksSame->DrawClone("SURF1");
1014 //      new TCanvas; tracksMixed->DrawClone("SURF1");
1015
1016         // some code to judge the relative contribution of the different correlation functions to the overall uncertainty
1017         Double_t sums[] = { 0, 0, 0 };
1018         Double_t errors[] = { 0, 0, 0 };
1019
1020         for (Int_t x=1; x<=tracksSame->GetNbinsX(); x++)
1021           for (Int_t y=1; y<=tracksSame->GetNbinsY(); y++)
1022           {
1023             sums[0] += tracksSame->GetBinContent(x, y);
1024             errors[0] += tracksSame->GetBinError(x, y);
1025             sums[1] += tracksMixed->GetBinContent(x, y);
1026             errors[1] += tracksMixed->GetBinError(x, y);
1027           }
1028
1029         tracksSame->Divide(tracksMixed);
1030
1031         for (Int_t x=1; x<=tracksSame->GetNbinsX(); x++)
1032           for (Int_t y=1; y<=tracksSame->GetNbinsY(); y++)
1033           {
1034             sums[2] += tracksSame->GetBinContent(x, y);
1035             errors[2] += tracksSame->GetBinError(x, y);
1036           }
1037
1038         for (Int_t x=0; x<3; x++)
1039           if (sums[x] > 0)
1040             errors[x] /= sums[x];
1041
1042         Printf("The correlation function %d %d has uncertainties %f %f %f (Ratio S/M %f)", multBin, vertexBin, errors[0], errors[1], errors[2], (errors[1] > 0) ? errors[0] / errors[1] : -1);
1043         // code to draw contributions
1044         /*
1045         TH1* proj = tracksSame->ProjectionX("projx", tracksSame->GetYaxis()->FindBin(-1.59), tracksSame->GetYaxis()->FindBin(1.59));
1046         proj->SetTitle(Form("Bin %d", vertexBin));
1047         proj->SetLineColor(vertexBin);
1048         proj->DrawCopy((vertexBin > 1) ? "SAME" : "");
1049         */
1050
1051         if (!totalTracks)
1052           totalTracks = (TH2*) tracksSame->Clone("totalTracks");
1053         else
1054           totalTracks->Add(tracksSame);
1055
1056         totalEvents += eventSameAll->GetBinContent(vertexBin, multBin);
1057
1058 //      new TCanvas; tracksMixed->DrawCopy("SURF1");
1059       }
1060
1061       delete tracksSame;
1062       delete tracksMixed;
1063
1064       nCorrelationFunctions++;
1065     }
1066   }
1067
1068   if (totalTracks) {
1069     Double_t sums[] = { 0, 0, 0 };
1070     Double_t errors[] = { 0, 0, 0 };
1071
1072     for (Int_t x=1; x<=totalTracks->GetNbinsX(); x++)
1073       for (Int_t y=1; y<=totalTracks->GetNbinsY(); y++)
1074       {
1075         sums[0] += totalTracks->GetBinContent(x, y);
1076         errors[0] += totalTracks->GetBinError(x, y);
1077       }
1078     if (sums[0] > 0)
1079       errors[0] /= sums[0];
1080
1081     if (normalizePerTrigger)
1082     {
1083       Printf("Dividing %f tracks by %d events (%d correlation function(s)) (error %f)", totalTracks->Integral(), totalEvents, nCorrelationFunctions, errors[0]);
1084       if (totalEvents > 0)
1085         totalTracks->Scale(1.0 / totalEvents);
1086     }
1087
1088     // normalizate to dphi width
1089     Float_t normalization = totalTracks->GetXaxis()->GetBinWidth(1);
1090     totalTracks->Scale(1.0 / normalization);
1091   }
1092
1093   delete trackSameAll;
1094   delete trackMixedAll;
1095   delete eventSameAll;
1096   delete eventMixedAll;
1097
1098 //   new TCanvas; normParameters->Draw();
1099
1100   TH1::AddDirectory(oldStatus);
1101
1102   return totalTracks;
1103 }
1104
1105 TH1* AliUEHist::GetTriggersAsFunctionOfMultiplicity(AliUEHist::CFStep step, Float_t ptLeadMin, Float_t ptLeadMax)
1106 {
1107   // returns the distribution of triggers as function of centrality/multiplicity
1108
1109   ResetBinLimits(fEventHist->GetGrid(step));
1110
1111   Int_t firstBin = fEventHist->GetGrid(step)->GetGrid()->GetAxis(0)->FindBin(ptLeadMin);
1112   Int_t lastBin = fEventHist->GetGrid(step)->GetGrid()->GetAxis(0)->FindBin(ptLeadMax);
1113   Printf("Using pT range %d --> %d", firstBin, lastBin);
1114   fEventHist->GetGrid(step)->GetGrid()->GetAxis(0)->SetRange(firstBin, lastBin);
1115
1116   TH1* eventHist = fEventHist->GetGrid(step)->Project(1);
1117
1118   ResetBinLimits(fEventHist->GetGrid(step));
1119
1120   return eventHist;
1121 }
1122
1123 /*
1124 TH2* AliUEHist::GetPtDistInPhiRegion(AliUEHist* mixed, AliUEHist::CFStep step, AliUEHist::Region region, Float_t ptLeadMin, Float_t ptLeadMax, Int_t multBinBegin, Int_t multBinEnd, Float_t phiBegin, Float_t phiEnd)
1125 {
1126   // Returns pT,assoc distribution in the given pt,trig, multiplicity, phi region
1127   // Does not use sum of ratios for mixed event correction (TODO to be improved)
1128   // returns a 2D histogram: deltaphi, deltaeta
1129   //
1130   // Parameters:
1131   //   mixed: AliUEHist containing mixed event corresponding to this object
1132 */
1133
1134 //____________________________________________________________________
1135 TH2* AliUEHist::GetSumOfRatios(AliUEHist* mixed, AliUEHist::CFStep step, AliUEHist::Region region, Float_t ptLeadMin, Float_t ptLeadMax, Int_t multBinBegin, Int_t multBinEnd, Bool_t etaNorm, Bool_t useVertexBins)
1136 {
1137   // Calls GetUEHist(...) for *each* multiplicity bin and performs a sum of ratios:
1138   // 1_N [ (same/mixed)_1 + (same/mixed)_2 + (same/mixed)_3 + ... ]
1139   // where N is the total number of events/trigger particles and the subscript is the multiplicity bin
1140   //
1141   // Can only be used for the 2D histogram at present
1142   //
1143   // Parameters:
1144   //   mixed: AliUEHist containing mixed event corresponding to this object
1145   //   <other parameters> : check documentation of AliUEHist::GetUEHist
1146
1147   TH2* totalTracks = 0;
1148   Int_t totalEvents = 0;
1149
1150   Int_t vertexBin = 1;
1151   TAxis* vertexAxis = fTrackHist[kToward]->GetGrid(0)->GetGrid()->GetAxis(5);
1152   if (useVertexBins && !vertexAxis)
1153   {
1154     Printf("Vertex axis requested but not available");
1155     return 0;
1156   }
1157
1158   // vertex bin loop
1159   while (1)
1160   {
1161     if (useVertexBins)
1162     {
1163       SetZVtxRange(vertexAxis->GetBinLowEdge(vertexBin) + 0.01, vertexAxis->GetBinUpEdge(vertexBin) - 0.01);
1164       mixed->SetZVtxRange(vertexAxis->GetBinLowEdge(vertexBin) + 0.01, vertexAxis->GetBinUpEdge(vertexBin) - 0.01);
1165       vertexBin++;
1166     }
1167
1168     // multiplicity loop
1169     Int_t multIter = multBinBegin;
1170     while (1)
1171     {
1172       Int_t multBinBeginLocal = multBinBegin;
1173       Int_t multBinEndLocal = multBinEnd;
1174
1175       if (multBinEnd >= multBinBegin)
1176       {
1177         multBinBeginLocal = multIter;
1178         multBinEndLocal = multIter;
1179         multIter++;
1180       }
1181
1182       Long64_t nEvents = 0;
1183       TH2* tracks = (TH2*) GetUEHist(step, region, ptLeadMin, ptLeadMax, multBinBeginLocal, multBinEndLocal, 1, etaNorm, &nEvents);
1184       // undo normalization
1185       tracks->Scale(nEvents);
1186       totalEvents += nEvents;
1187
1188       TH2* mixedTwoD = (TH2*) mixed->GetUEHist(step, region, ptLeadMin, ptLeadMax, multBinBeginLocal, multBinEndLocal, 1, etaNorm);
1189
1190       // asssume flat in dphi, gain in statistics
1191   //     TH1* histMixedproj = mixedTwoD->ProjectionY();
1192   //     histMixedproj->Scale(1.0 / mixedTwoD->GetNbinsX());
1193   //
1194   //     for (Int_t x=1; x<=mixedTwoD->GetNbinsX(); x++)
1195   //       for (Int_t y=1; y<=mixedTwoD->GetNbinsY(); y++)
1196   //    mixedTwoD->SetBinContent(x, y, histMixedproj->GetBinContent(y));
1197
1198   //       delete histMixedproj;
1199
1200       // get mixed event normalization by assuming full acceptance at deta of 0 (only works for flat dphi)
1201 /*      Double_t mixedNorm = mixedTwoD->Integral(1, mixedTwoD->GetNbinsX(), mixedTwoD->GetYaxis()->FindBin(-0.01), mixedTwoD->GetYaxis()->FindBin(0.01));
1202       mixedNorm /= mixedTwoD->GetNbinsX() * (mixedTwoD->GetYaxis()->FindBin(0.01) - mixedTwoD->GetYaxis()->FindBin(-0.01) + 1);
1203       tracks->Scale(mixedNorm);*/
1204
1205       tracks->Scale(mixedTwoD->Integral() / tracks->Integral());
1206
1207       tracks->Divide(mixedTwoD);
1208
1209       delete mixedTwoD;
1210
1211       if (!totalTracks)
1212         totalTracks = tracks;
1213       else
1214       {
1215         totalTracks->Add(tracks);
1216         delete tracks;
1217       }
1218
1219       if (multIter > multBinEnd)
1220         break;
1221     }
1222
1223     if (!useVertexBins || vertexBin > vertexAxis->GetNbins())
1224       break;
1225   }
1226
1227   if (useVertexBins)
1228     totalEvents = vertexAxis->GetNbins();
1229
1230   Printf("Dividing %f tracks by %d events", totalTracks->Integral(), totalEvents);
1231   if (totalEvents > 0)
1232     totalTracks->Scale(1.0 / totalEvents);
1233
1234   return totalTracks;
1235 }
1236
1237 //____________________________________________________________________
1238 TH1* AliUEHist::GetPtHist(CFStep step, Region region, Float_t ptLeadMin, Float_t ptLeadMax, Int_t multBinBegin, Int_t multBinEnd, Float_t phiMin, Float_t phiMax, Float_t etaMin, Float_t etaMax, Bool_t skipPhiNormalization)
1239 {
1240   // returns a histogram projected to pT,assoc with the provided cuts
1241
1242    // unzoom all axes
1243   ResetBinLimits(fTrackHist[region]->GetGrid(step));
1244   ResetBinLimits(fEventHist->GetGrid(step));
1245
1246   TH1D* tracks = 0;
1247
1248     // the efficiency to have find an event depends on leading pT and this is not corrected for because anyway per bin we calculate tracks over events
1249   // therefore the number density must be calculated here per leading pT bin and then summed
1250
1251   if (multBinEnd > multBinBegin)
1252     Printf("Using multiplicity range %d --> %d", multBinBegin, multBinEnd);
1253   fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(3)->SetRange(multBinBegin, multBinEnd);
1254   fEventHist->GetGrid(step)->GetGrid()->GetAxis(1)->SetRange(multBinBegin, multBinEnd);
1255
1256   fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(0)->SetRangeUser(etaMin, etaMax);
1257   fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(4)->SetRangeUser(phiMin, phiMax);
1258
1259   // get real phi cuts due to binning
1260   Float_t phiMinReal = fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(4)->GetBinLowEdge(fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(4)->GetFirst());
1261   Float_t phiMaxReal = fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(4)->GetBinUpEdge(fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(4)->GetLast());
1262   Printf("phi min = %.2f and max = %.2f requested; due to binning range will be from %.2f to %.2f, i.e. a %.0f%% larger window", phiMin, phiMax, phiMinReal, phiMaxReal, (phiMaxReal - phiMinReal) / (phiMax - phiMin) * 100 - 100);
1263
1264   Int_t firstBin = fTrackHist[region]->GetAxis(2, step)->FindBin(ptLeadMin);
1265   Int_t lastBin = fTrackHist[region]->GetAxis(2, step)->FindBin(ptLeadMax);
1266
1267   TH1D* events = fEventHist->ShowProjection(0, step);
1268
1269   for (Int_t bin=firstBin; bin<=lastBin; bin++)
1270   {
1271     fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(2)->SetRange(bin, bin);
1272
1273     // project to pT,assoc
1274     TH1D* tracksTmp = (TH1D*) fTrackHist[region]->GetGrid(step)->Project(1);
1275
1276     Printf("Calculated histogram in bin %d --> %f tracks", bin, tracksTmp->Integral());
1277     fTrackHist[region]->GetGrid(step)->SetRangeUser(2, 0, -1);
1278
1279     // normalize to get a density (deta dphi)
1280     Float_t normalization = 1;
1281     TAxis* axis = fTrackHist[region]->GetGrid(step)->GetAxis(0);
1282     if (strcmp(axis->GetTitle(), "#eta") == 0)
1283     {
1284       Printf("Normalizing using eta-axis range");
1285       normalization *= axis->GetBinUpEdge(axis->GetLast()) - axis->GetBinLowEdge(axis->GetFirst());
1286     }
1287     else
1288     {
1289       Printf("Normalizating assuming accepted range of |eta| < 0.8");
1290       normalization *= 0.8 * 2;
1291     }
1292
1293     // dphi
1294     if (!skipPhiNormalization)
1295       normalization *= phiMaxReal - phiMinReal;
1296
1297     //Printf("Normalization: %f", normalization);
1298     tracksTmp->Scale(1.0 / normalization);
1299
1300     // and now dpT (bins have different width)
1301     for (Int_t i=1; i<=tracksTmp->GetNbinsX(); i++)
1302     {
1303       tracksTmp->SetBinContent(i, tracksTmp->GetBinContent(i) / tracksTmp->GetXaxis()->GetBinWidth(i));
1304       tracksTmp->SetBinError  (i, tracksTmp->GetBinError(i) / tracksTmp->GetXaxis()->GetBinWidth(i));
1305     }
1306
1307     Int_t nEvents = (Int_t) events->GetBinContent(bin);
1308     if (nEvents > 0)
1309       tracksTmp->Scale(1.0 / nEvents);
1310     Printf("Calculated histogram in bin %d --> %d events", bin, nEvents);
1311
1312     if (!tracks)
1313       tracks = tracksTmp;
1314     else
1315     {
1316       tracks->Add(tracksTmp);
1317       delete tracksTmp;
1318     }
1319   }
1320
1321   delete events;
1322
1323   ResetBinLimits(fTrackHist[region]->GetGrid(step));
1324   ResetBinLimits(fEventHist->GetGrid(step));
1325
1326   return tracks;
1327 }
1328
1329 void AliUEHist::MultiplyHistograms(THnSparse* grid, THnSparse* target, TH1* histogram, Int_t var1, Int_t var2)
1330 {
1331   // multiplies <grid> with <histogram> and put the result in <target>
1332   // <grid> has usually more dimensions than histogram. The axis which are used to choose the value
1333   // from <histogram> are given in <var1> and <var2>
1334   //
1335   // if <histogram> is 0, just copies content from step1 to step2
1336
1337   // clear target histogram
1338   if (grid != target)
1339     target->Reset();
1340
1341   if (histogram != 0)
1342   {
1343     if (grid->GetAxis(var1)->GetNbins() != histogram->GetNbinsX())
1344       AliFatal(Form("Invalid binning (var1): %d %d", grid->GetAxis(var1)->GetNbins(), histogram->GetNbinsX()));
1345
1346     if (var2 >= 0 && grid->GetAxis(var2)->GetNbins() != histogram->GetNbinsY())
1347       AliFatal(Form("Invalid binning (var2): %d %d", grid->GetAxis(var2)->GetNbins(), histogram->GetNbinsY()));
1348   }
1349
1350   if (grid->GetNdimensions() > 6)
1351     AliFatal("Too many dimensions in THnSparse");
1352
1353   Int_t bins[6];
1354
1355   // optimized implementation
1356   for (Int_t binIdx = 0; binIdx < grid->GetNbins(); binIdx++)
1357   {
1358     Double_t value = grid->GetBinContent(binIdx, bins);
1359     Double_t error = grid->GetBinError(binIdx);
1360
1361     if (histogram != 0)
1362     {
1363       if (var2 < 0)
1364       {
1365         value *= histogram->GetBinContent(bins[var1]);
1366         error *= histogram->GetBinContent(bins[var1]);
1367       }
1368       else
1369       {
1370         value *= histogram->GetBinContent(bins[var1], bins[var2]);
1371         error *= histogram->GetBinContent(bins[var1], bins[var2]);
1372       }
1373     }
1374
1375     target->SetBinContent(bins, value);
1376     target->SetBinError(bins, error);
1377   }
1378 }
1379
1380 //____________________________________________________________________
1381 void AliUEHist::CorrectTracks(CFStep step1, CFStep step2, TH1* trackCorrection, Int_t var1, Int_t var2)
1382 {
1383   // corrects from step1 to step2 by multiplying the tracks with trackCorrection
1384   // trackCorrection can be a function of eta (var1 == 0), pT (var1 == 1), leading pT (var1 == 2), multiplicity (var1 == 3), delta phi (var1 == 4)
1385   // if var2 >= 0 a two dimension correction is assumed in trackCorrection
1386   //
1387   // if trackCorrection is 0, just copies content from step1 to step2
1388
1389   for (UInt_t region=0; region<fkRegions; region++)
1390     CorrectTracks(step1, step2, region, trackCorrection, var1, var2);
1391 }
1392
1393 //____________________________________________________________________
1394 void AliUEHist::CorrectTracks(CFStep step1, CFStep step2, Int_t region, TH1* trackCorrection, Int_t var1, Int_t var2)
1395 {
1396   //
1397   // see documentation of CorrectTracks above
1398   //
1399
1400   if (!fTrackHist[region])
1401     return;
1402
1403   THnSparse* grid = fTrackHist[region]->GetGrid(step1)->GetGrid();
1404   THnSparse* target = fTrackHist[region]->GetGrid(step2)->GetGrid();
1405
1406   Float_t entriesBefore = grid->GetEntries();
1407
1408   MultiplyHistograms(grid, target, trackCorrection, var1, var2);
1409
1410   Printf("AliUEHist::CorrectTracks: Corrected from %f (step %d) to %f (step %d) entries. Correction histogram: %f entries (integral: %f)", entriesBefore, step1, target->GetEntries(), step2, (trackCorrection) ? trackCorrection->GetEntries() : -1.0, (trackCorrection) ? trackCorrection->Integral() : -1.0);
1411 }
1412
1413 //____________________________________________________________________
1414 void AliUEHist::CorrectEvents(CFStep step1, CFStep step2, TH1* eventCorrection, Int_t var1, Int_t var2)
1415 {
1416   // corrects from step1 to step2 by multiplying the events with eventCorrection
1417   // eventCorrection is as function of leading pT (var == 0) or multiplicity (var == 1)
1418   //
1419   // if eventCorrection is 0, just copies content from step1 to step2
1420
1421   AliCFGridSparse* grid = fEventHist->GetGrid(step1);
1422   AliCFGridSparse* target = fEventHist->GetGrid(step2);
1423
1424   Float_t entriesBefore = grid->GetEntries();
1425
1426   MultiplyHistograms(grid->GetGrid(), target->GetGrid(), eventCorrection, var1, var2);
1427
1428   Printf("AliUEHist::CorrectEvents: Corrected from %f (step %d) to %f (step %d) entries. Correction histogram: %f entries (integral: %f)", entriesBefore, step1, target->GetEntries(), step2, (eventCorrection) ? eventCorrection->GetEntries() : -1.0, (eventCorrection) ? eventCorrection->Integral() : -1.0);
1429 }
1430
1431 //____________________________________________________________________
1432 void AliUEHist::Correct(AliUEHist* corrections)
1433 {
1434   // applies the given corrections to extract from the step kCFStepReconstructed all previous steps
1435   //
1436   // in this object the data is expected in the step kCFStepReconstructed
1437
1438   if (fHistogramType.Length() == 0)
1439   {
1440     Printf("W-AliUEHist::Correct: fHistogramType not defined. Guessing histogram type...");
1441     if (fTrackHist[kToward]->GetNVar() < 5)
1442     {
1443       if (strcmp(fTrackHist[kToward]->GetTitle(), "d^{2}N_{ch}/d#varphid#eta") == 0)
1444         fHistogramType = "NumberDensitypT";
1445       else if (strcmp(fTrackHist[kToward]->GetTitle(), "d^{2}#Sigma p_{T}/d#varphid#eta") == 0)
1446         fHistogramType = "SumpT";
1447     }
1448     else if (fTrackHist[kToward]->GetNVar() == 5)
1449     {
1450       if (strcmp(fTrackHist[kToward]->GetGrid(0)->GetGrid()->GetAxis(3)->GetTitle(), "multiplicity") == 0)
1451         fHistogramType = "NumberDensityPhi";
1452       else if (strcmp(fTrackHist[kToward]->GetGrid(0)->GetGrid()->GetAxis(3)->GetTitle(), "centrality") == 0)
1453         fHistogramType = "NumberDensityPhiCentrality";
1454     }
1455
1456     if (fHistogramType.Length() == 0)
1457       AliFatal("Cannot figure out histogram type. Try setting it manually...");
1458   }
1459
1460   Printf("AliUEHist::Correct: Correcting %s...", fHistogramType.Data());
1461
1462   if (strcmp(fHistogramType, "NumberDensitypT") == 0 || strcmp(fHistogramType, "NumberDensityPhi") == 0 || strcmp(fHistogramType, "SumpT") == 0)
1463   {
1464     // ---- track level
1465
1466     // bias due to migration in leading pT (because the leading particle is not reconstructed, and the subleading is used)
1467     // extracted as function of leading pT
1468     Bool_t biasFromMC = kFALSE;
1469     const char* projAxis = "z";
1470     Int_t secondBin = -1;
1471
1472     if (strcmp(fHistogramType, "NumberDensityPhi") == 0)
1473     {
1474       projAxis = "yz";
1475       secondBin = 4;
1476     }
1477
1478     for (UInt_t region = 0; region < fkRegions; region++)
1479     {
1480       if (!fTrackHist[region])
1481         continue;
1482
1483       TH1* leadingBiasTracks =  0;
1484       if (biasFromMC)
1485       {
1486         leadingBiasTracks = (TH1*) corrections->GetBias(kCFStepReconstructed, kCFStepTracked, region, projAxis, 0, -1, 1); // from MC
1487         Printf("WARNING: Using MC bias correction");
1488       }
1489       else
1490         leadingBiasTracks = (TH1*) GetBias(kCFStepBiasStudy, kCFStepReconstructed, region, projAxis, 0, -1, 1);          // from data
1491
1492       CorrectTracks(kCFStepReconstructed, kCFStepTracked, region, leadingBiasTracks, 2, secondBin);
1493       if (region == kMin && fCombineMinMax)
1494       {
1495         CorrectTracks(kCFStepReconstructed, kCFStepTracked, kMax, leadingBiasTracks, 2, secondBin);
1496         delete leadingBiasTracks;
1497         break;
1498       }
1499       delete leadingBiasTracks;
1500     }
1501
1502     TH1* leadingBiasEvents = 0;
1503     if (biasFromMC)
1504       leadingBiasEvents = (TH1*) corrections->GetBias(kCFStepReconstructed, kCFStepTracked, kToward, projAxis, 0, -1, 2); // from MC
1505     else
1506       leadingBiasEvents = (TH1*) GetBias(kCFStepBiasStudy, kCFStepReconstructed, kToward, projAxis, 0, -1, 2);          // from data
1507
1508     //new TCanvas; leadingBiasEvents->DrawCopy();
1509     CorrectEvents(kCFStepReconstructed, kCFStepTracked, leadingBiasEvents, 0);
1510     delete leadingBiasEvents;
1511
1512     // --- contamination correction ---
1513
1514     TH2D* contamination = corrections->GetTrackingContamination();
1515     if (corrections->fContaminationEnhancement)
1516     {
1517       Printf("Applying contamination enhancement");
1518
1519       for (Int_t x=1; x<=contamination->GetNbinsX(); x++)
1520         for (Int_t y=1; y<=contamination->GetNbinsX(); y++)
1521           contamination->SetBinContent(x, y, contamination->GetBinContent(x, y) * corrections->fContaminationEnhancement->GetBinContent(corrections->fContaminationEnhancement->GetXaxis()->FindBin(contamination->GetYaxis()->GetBinCenter(y))));
1522     }
1523     CorrectTracks(kCFStepTracked, kCFStepTrackedOnlyPrim, contamination, 0, 1);
1524     CorrectEvents(kCFStepTracked, kCFStepTrackedOnlyPrim, 0, 0);
1525     delete contamination;
1526
1527     // --- efficiency correction ---
1528     // correct single-particle efficiency for associated particles
1529     // in addition correct for efficiency on trigger particles (tracks AND events)
1530
1531     TH1* efficiencyCorrection = corrections->GetTrackingEfficiencyCorrection();
1532     // use kCFStepVertex as a temporary step (will be overwritten later anyway)
1533     CorrectTracks(kCFStepTrackedOnlyPrim, kCFStepVertex, efficiencyCorrection, 0, 1);
1534     delete efficiencyCorrection;
1535
1536     // correct pT,T
1537     efficiencyCorrection = corrections->GetTrackEfficiency(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, 1, -1, 2);
1538     CorrectEvents(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, efficiencyCorrection, 0);
1539     CorrectTracks(kCFStepVertex, kCFStepAnaTopology, efficiencyCorrection, 2);
1540     delete efficiencyCorrection;
1541
1542     // copy
1543     CorrectTracks(kCFStepAnaTopology, kCFStepVertex, 0, -1);
1544     CorrectEvents(kCFStepAnaTopology, kCFStepVertex, 0, 0);
1545
1546     // vertex correction on the level of events as function of multiplicity, weighting tracks and events with the same factor
1547     // practically independent of low pT cut
1548     TH1D* vertexCorrection = (TH1D*) corrections->GetEventEfficiency(kCFStepVertex, kCFStepTriggered, 1);
1549
1550     // convert stage from true multiplicity to observed multiplicity by simple conversion factor
1551     TH1D* vertexCorrectionObs = (TH1D*) vertexCorrection->Clone("vertexCorrection2");
1552     vertexCorrectionObs->Reset();
1553
1554     TF1* func = new TF1("func", "[1]+[0]/(x-[2])");
1555     // some defaults
1556     func->SetParameters(0.1, 1, -0.7);
1557     vertexCorrection->Fit(func, "0I", "", 0, 3);
1558     for (Int_t i=1; i<=vertexCorrectionObs->GetNbinsX(); i++)
1559     {
1560       Float_t xPos = 1.0 / 0.77 * vertexCorrectionObs->GetXaxis()->GetBinCenter(i);
1561       if (xPos < 3)
1562         vertexCorrectionObs->SetBinContent(i, func->Eval(xPos));
1563       else
1564         vertexCorrectionObs->SetBinContent(i, vertexCorrection->Interpolate(xPos));
1565     }
1566
1567     #if 0
1568     new TCanvas;
1569     vertexCorrection->DrawCopy();
1570     vertexCorrectionObs->SetLineColor(2);
1571     vertexCorrectionObs->DrawCopy("same");
1572     func->SetRange(0, 4);
1573     func->DrawClone("same");
1574     #endif
1575
1576     CorrectTracks(kCFStepVertex, kCFStepTriggered, vertexCorrectionObs, 3);
1577     CorrectEvents(kCFStepVertex, kCFStepTriggered, vertexCorrectionObs, 1);
1578     delete vertexCorrectionObs;
1579     delete vertexCorrection;
1580     delete func;
1581
1582     // copy
1583     CorrectTracks(kCFStepTriggered, kCFStepAll, 0, -1);
1584     CorrectEvents(kCFStepTriggered, kCFStepAll, 0, 0);
1585   }
1586   else if (strcmp(fHistogramType, "NumberDensityPhiCentrality") == 0)
1587   {
1588     if (fTrackHist[0]->GetNVar() <= 5)
1589     {
1590       // do corrections copying between steps
1591 //       CFStep step = kCFStepReconstructed;
1592       CFStep step = kCFStepBiasStudy;
1593
1594       // copy
1595       CorrectTracks(step, kCFStepTracked, 0, -1);
1596       CorrectEvents(step, kCFStepTracked, 0, -1);
1597
1598       // Dont use eta in the following, because it is a Delta-eta axis
1599
1600       // contamination correction
1601       // correct single-particle contamination for associated particles
1602
1603       TH1* contamination = corrections->GetTrackingContamination(1);
1604
1605       if (0)
1606       {
1607         Printf("Applying contamination enhancement");
1608
1609         for (Int_t bin = 1; bin <= contamination->GetNbinsX(); bin++)
1610         {
1611           printf("%f", contamination->GetBinContent(bin));
1612           if (contamination->GetBinContent(bin) > 0)
1613             contamination->SetBinContent(bin, 1.0 + 1.1 * (contamination->GetBinContent(bin) - 1.0));
1614           printf(" --> %f\n", contamination->GetBinContent(bin));
1615         }
1616       }
1617
1618       CorrectTracks(kCFStepTracked, kCFStepTrackedOnlyPrim, contamination, 1);
1619       delete contamination;
1620
1621       // correct for additional contamination due to trigger particle around phi ~ 0
1622       TH2* correlatedContamination = corrections->GetCorrelatedContamination();
1623       if (0)
1624       {
1625         Printf("Applying contamination enhancement");
1626
1627         for (Int_t bin = 1; bin <= correlatedContamination->GetNbinsX(); bin++)
1628           for (Int_t bin2 = 1; bin2 <= correlatedContamination->GetNbinsY(); bin2++)
1629           {
1630             printf("%f", correlatedContamination->GetBinContent(bin, bin2));
1631             if (correlatedContamination->GetBinContent(bin, bin2) > 0)
1632               correlatedContamination->SetBinContent(bin, bin2, 1.0 + 1.1 * (correlatedContamination->GetBinContent(bin, bin2) - 1.0));
1633             printf(" --> %f\n", correlatedContamination->GetBinContent(bin, bin2));
1634           }
1635       }
1636
1637 //       new TCanvas; correlatedContamination->DrawCopy("COLZ");
1638       CorrectCorrelatedContamination(kCFStepTrackedOnlyPrim, 0, correlatedContamination);
1639 //       Printf("\n\n\nWARNING ---> SKIPPING CorrectCorrelatedContamination\n\n\n");
1640
1641       delete correlatedContamination;
1642
1643       // TODO correct for contamination of trigger particles (for tracks AND events)
1644       CorrectEvents(kCFStepTracked, kCFStepTrackedOnlyPrim, 0, 0);
1645
1646       // --- efficiency correction ---
1647       // correct single-particle efficiency for associated particles
1648       // in addition correct for efficiency on trigger particles (tracks AND events)
1649
1650       // in bins of pT and centrality
1651       TH1* efficiencyCorrection = corrections->GetTrackingEfficiencyCorrectionCentrality();
1652 //       new TCanvas; efficiencyCorrection->DrawCopy("COLZ");
1653       // use kCFStepAnaTopology as a temporary step
1654       CorrectTracks(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, efficiencyCorrection, 1, 3);
1655       delete efficiencyCorrection;
1656
1657       // correct pT,T in bins of pT and centrality
1658       efficiencyCorrection = corrections->GetTrackEfficiency(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, 1, 3, 2);
1659       CorrectEvents(kCFStepTrackedOnlyPrim, kCFStepVertex, efficiencyCorrection, 0, 1);
1660       CorrectTracks(kCFStepAnaTopology, kCFStepVertex, efficiencyCorrection, 2, 3);
1661       delete efficiencyCorrection;
1662
1663       // no correction for vertex finding efficiency and trigger efficiency needed in PbPb
1664       // copy
1665       CorrectTracks(kCFStepVertex, kCFStepAll, 0, -1);
1666       CorrectEvents(kCFStepVertex, kCFStepAll, 0, -1);
1667     }
1668     else
1669     {
1670       // with 6 axes there is not enough memory, do the corrections in-place
1671       Printf(">>>>>>>> Applying corrections in place to reduce memory consumption");
1672       CFStep step = kCFStepBiasStudy;
1673 //       CFStep step = kCFStepReconstructed;
1674
1675       // Dont use eta in the following, because it is a Delta-eta axis
1676
1677       // --- contamination correction ---
1678       // correct single-particle contamination for associated particles
1679       // correct contamination for trigger particles (tracks AND events)
1680
1681       // in bins of p,T
1682       TH1* contamination = corrections->GetTrackingContamination(1);
1683
1684       if (0)
1685       {
1686         Printf("Applying contamination enhancement");
1687
1688         for (Int_t bin = 1; bin <= contamination->GetNbinsX(); bin++)
1689         {
1690           printf("%f", contamination->GetBinContent(bin));
1691           if (contamination->GetBinContent(bin) > 0)
1692             contamination->SetBinContent(bin, 1.0 + 1.1 * (contamination->GetBinContent(bin) - 1.0));
1693           printf(" --> %f\n", contamination->GetBinContent(bin));
1694         }
1695       }
1696
1697       // correct pT,A in bins of pT
1698       CorrectTracks(step, step, contamination, 1);
1699       delete contamination;
1700
1701       // correct pT,T in bins of pT (for tracks AND events)
1702       contamination = corrections->GetTrackEfficiency(kCFStepTracked, kCFStepTrackedOnlyPrim, 1, -1, 2);
1703       new TCanvas; contamination->DrawCopy();
1704       CorrectEvents(step, step, contamination, 0);
1705       CorrectTracks(step, step, contamination, 2);
1706       delete contamination;
1707
1708       // correct for additional contamination due to trigger particle around phi ~ 0
1709       if (0)
1710       {
1711         TH2* correlatedContamination = corrections->GetCorrelatedContamination();
1712         if (1)
1713         {
1714           Printf("Applying contamination enhancement");
1715
1716           for (Int_t bin = 1; bin <= correlatedContamination->GetNbinsX(); bin++)
1717             for (Int_t bin2 = 1; bin2 <= correlatedContamination->GetNbinsY(); bin2++)
1718             {
1719               printf("%f", correlatedContamination->GetBinContent(bin, bin2));
1720               if (correlatedContamination->GetBinContent(bin, bin2) > 0)
1721                 correlatedContamination->SetBinContent(bin, bin2, 1.0 + 1.1 * (correlatedContamination->GetBinContent(bin, bin2) - 1.0));
1722               printf(" --> %f\n", correlatedContamination->GetBinContent(bin, bin2));
1723             }
1724         }
1725
1726         // new TCanvas; correlatedContamination->DrawCopy("COLZ");
1727         CorrectCorrelatedContamination(step, 0, correlatedContamination);
1728
1729         delete correlatedContamination;
1730       }
1731       else
1732         Printf("\n\n\nWARNING ---> SKIPPING CorrectCorrelatedContamination\n\n\n");
1733
1734       // --- tracking efficiency correction ---
1735       // correct single-particle efficiency for associated particles
1736       // correct for efficiency on trigger particles (tracks AND events)
1737
1738       // in bins of pT and centrality
1739       TH1* efficiencyCorrection = corrections->GetTrackingEfficiencyCorrectionCentrality();
1740       // new TCanvas; efficiencyCorrection->DrawCopy("COLZ");
1741
1742       // correct pT,A in bins of pT and centrality
1743       CorrectTracks(step, step, efficiencyCorrection, 1, 3);
1744       delete efficiencyCorrection;
1745
1746       // correct pT,T in bins of pT and centrality
1747       efficiencyCorrection = corrections->GetTrackEfficiency(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, 1, 3, 2);
1748       CorrectEvents(step, step, efficiencyCorrection, 0, 1);
1749       CorrectTracks(step, step, efficiencyCorrection, 2, 3);
1750
1751       delete efficiencyCorrection;
1752     }
1753   }
1754   else
1755     AliFatal(Form("Unknown histogram for correction: %s", GetTitle()));
1756 }
1757
1758 //____________________________________________________________________
1759 THnBase* AliUEHist::GetTrackEfficiencyND(CFStep step1, CFStep step2)
1760 {
1761   // creates a track-level efficiency by dividing step2 by step1
1762   // in all dimensions but the particle species one
1763
1764   AliCFContainer* sourceContainer = fTrackHistEfficiency;
1765   // step offset because we start with kCFStepAnaTopology
1766   step1 = (CFStep) ((Int_t) step1 - (Int_t) kCFStepAnaTopology);
1767   step2 = (CFStep) ((Int_t) step2 - (Int_t) kCFStepAnaTopology);
1768
1769   ResetBinLimits(sourceContainer->GetGrid(step1));
1770   ResetBinLimits(sourceContainer->GetGrid(step2));
1771
1772   if (fEtaMax > fEtaMin)
1773   {
1774     Printf("Restricted eta-range to %f %f", fEtaMin, fEtaMax);
1775     sourceContainer->GetGrid(step1)->SetRangeUser(0, fEtaMin, fEtaMax);
1776     sourceContainer->GetGrid(step2)->SetRangeUser(0, fEtaMin, fEtaMax);
1777   }
1778
1779   Int_t dimensions[] = { 0, 1, 3, 4 };
1780   THnBase* generated = sourceContainer->GetGrid(step1)->GetGrid()->ProjectionND(4, dimensions);
1781   THnBase* measured = sourceContainer->GetGrid(step2)->GetGrid()->ProjectionND(4, dimensions);
1782
1783 //   Printf("%d %d %f %f", step1, step2, generated->GetEntries(), measured->GetEntries());
1784
1785   ResetBinLimits(sourceContainer->GetGrid(step1));
1786   ResetBinLimits(sourceContainer->GetGrid(step2));
1787
1788   measured->Divide(measured, generated, 1, 1, "B");
1789
1790   delete generated;
1791
1792   return measured;
1793 }
1794
1795 //____________________________________________________________________
1796 TH1* AliUEHist::GetTrackEfficiency(CFStep step1, CFStep step2, Int_t axis1, Int_t axis2, Int_t source, Int_t axis3)
1797 {
1798   // creates a track-level efficiency by dividing step2 by step1
1799   // projected to axis1 and axis2 (optional if >= 0)
1800   //
1801   // source: 0 = fTrackHist; 1 = fTrackHistEfficiency; 2 = fTrackHistEfficiency rebinned for pT,T / pT,lead binning
1802
1803   // integrate over regions
1804   // cache it for efficiency (usually more than one efficiency is requested)
1805
1806   AliCFContainer* sourceContainer = 0;
1807
1808   if (source == 0)
1809   {
1810     if (!fCache)
1811     {
1812       fCache = (AliCFContainer*) fTrackHist[0]->Clone();
1813       for (UInt_t i = 1; i < fkRegions; i++)
1814         if (fTrackHist[i])
1815           fCache->Add(fTrackHist[i]);
1816     }
1817     sourceContainer = fCache;
1818   }
1819   else if (source == 1 || source == 2)
1820   {
1821     sourceContainer = fTrackHistEfficiency;
1822     // step offset because we start with kCFStepAnaTopology
1823     step1 = (CFStep) ((Int_t) step1 - (Int_t) kCFStepAnaTopology);
1824     step2 = (CFStep) ((Int_t) step2 - (Int_t) kCFStepAnaTopology);
1825   }
1826   else
1827     return 0;
1828
1829   // reset all limits and set the right ones except those in axis1, axis2 and axis3
1830   ResetBinLimits(sourceContainer->GetGrid(step1));
1831   ResetBinLimits(sourceContainer->GetGrid(step2));
1832   if (fEtaMax > fEtaMin && axis1 != 0 && axis2 != 0 && axis3 != 0)
1833   {
1834     Printf("Restricted eta-range to %f %f", fEtaMin, fEtaMax);
1835     sourceContainer->GetGrid(step1)->SetRangeUser(0, fEtaMin, fEtaMax);
1836     sourceContainer->GetGrid(step2)->SetRangeUser(0, fEtaMin, fEtaMax);
1837   }
1838   if (fPtMax > fPtMin && axis1 != 1 && axis2 != 1 && axis3 != 1)
1839   {
1840     Printf("Restricted pt-range to %f %f", fPtMin, fPtMax);
1841     sourceContainer->GetGrid(step1)->SetRangeUser(1, fPtMin, fPtMax);
1842     sourceContainer->GetGrid(step2)->SetRangeUser(1, fPtMin, fPtMax);
1843   }
1844   if (fPartSpecies != -1 && axis1 != 2 && axis2 != 2 && axis3 != 2)
1845   {
1846     Printf("Restricted to particle species %d", fPartSpecies);
1847     sourceContainer->GetGrid(step1)->SetRangeUser(2, fPartSpecies, fPartSpecies);
1848     sourceContainer->GetGrid(step2)->SetRangeUser(2, fPartSpecies, fPartSpecies);
1849   }
1850   if (fCentralityMax > fCentralityMin && axis1 != 3 && axis2 != 3 && axis3 != 3)
1851   {
1852     Printf("Restricted centrality range to %f %f", fCentralityMin, fCentralityMax);
1853     sourceContainer->GetGrid(step1)->SetRangeUser(3, fCentralityMin, fCentralityMax);
1854     sourceContainer->GetGrid(step2)->SetRangeUser(3, fCentralityMin, fCentralityMax);
1855   }
1856   if (fZVtxMax > fZVtxMin && axis1 != 4 && axis2 != 4 && axis3 != 4)
1857   {
1858     Printf("Restricted z-vtx range to %f %f", fZVtxMin, fZVtxMax);
1859     sourceContainer->GetGrid(step1)->SetRangeUser(4, fZVtxMin, fZVtxMax);
1860     sourceContainer->GetGrid(step2)->SetRangeUser(4, fZVtxMin, fZVtxMax);
1861   }
1862
1863   TH1* measured = 0;
1864   TH1* generated = 0;
1865
1866   if (axis3 >= 0)
1867   {
1868     generated = sourceContainer->Project(step1, axis1, axis2, axis3);
1869     measured = sourceContainer->Project(step2, axis1, axis2, axis3);
1870   }
1871   else if (axis2 >= 0)
1872   {
1873     generated = sourceContainer->Project(step1, axis1, axis2);
1874     measured = sourceContainer->Project(step2, axis1, axis2);
1875   }
1876   else
1877   {
1878     generated = sourceContainer->Project(step1, axis1);
1879     measured = sourceContainer->Project(step2, axis1);
1880   }
1881
1882   // check for bins with less than 50 entries, print warning
1883   Int_t binBegin[3];
1884   Int_t binEnd[3];
1885
1886   binBegin[0] = 1;
1887   binBegin[1] = 1;
1888   binBegin[2] = 1;
1889
1890   binEnd[0] = generated->GetNbinsX();
1891   binEnd[1] = generated->GetNbinsY();
1892   binEnd[2] = generated->GetNbinsZ();
1893
1894   if (fEtaMax > fEtaMin)
1895   {
1896     if (axis1 == 0)
1897     {
1898       binBegin[0] = generated->GetXaxis()->FindBin(fEtaMin);
1899       binEnd[0]   = generated->GetXaxis()->FindBin(fEtaMax);
1900     }
1901     if (axis2 == 0)
1902     {
1903       binBegin[1] = generated->GetYaxis()->FindBin(fEtaMin);
1904       binEnd[1]   = generated->GetYaxis()->FindBin(fEtaMax);
1905     }
1906     if (axis3 == 0)
1907     {
1908       binBegin[2] = generated->GetZaxis()->FindBin(fEtaMin);
1909       binEnd[2]   = generated->GetZaxis()->FindBin(fEtaMax);
1910     }
1911   }
1912
1913   if (fPtMax > fPtMin)
1914   {
1915     // TODO this is just checking up to 15 for now
1916     Float_t ptMax = TMath::Min((Float_t) 15., fPtMax);
1917     if (axis1 == 1)
1918     {
1919       binBegin[0] = generated->GetXaxis()->FindBin(fPtMin);
1920       binEnd[0]   = generated->GetXaxis()->FindBin(ptMax);
1921     }
1922     if (axis2 == 1)
1923     {
1924       binBegin[1] = generated->GetYaxis()->FindBin(fPtMin);
1925       binEnd[1]   = generated->GetYaxis()->FindBin(ptMax);
1926     }
1927     if (axis3 == 1)
1928     {
1929       binBegin[2] = generated->GetZaxis()->FindBin(fPtMin);
1930       binEnd[2]   = generated->GetZaxis()->FindBin(ptMax);
1931     }
1932   }
1933
1934   Int_t total = 0;
1935   Int_t count = 0;
1936   Int_t vars[3];
1937
1938   for (Int_t i=0; i<3; i++)
1939     vars[i] = binBegin[i];
1940
1941   const Int_t limit = 50;
1942   while (1)
1943   {
1944     if (generated->GetDimension() == 1 && generated->GetBinContent(vars[0]) < limit)
1945     {
1946       Printf("Empty bin at %s=%.2f (%.2f entries)", generated->GetXaxis()->GetTitle(), generated->GetXaxis()->GetBinCenter(vars[0]), generated->GetBinContent(vars[0]));
1947       count++;
1948     }
1949     else if (generated->GetDimension() == 2 && generated->GetBinContent(vars[0], vars[1]) < limit)
1950     {
1951       Printf("Empty bin at %s=%.2f %s=%.2f (%.2f entries)",
1952         generated->GetXaxis()->GetTitle(), generated->GetXaxis()->GetBinCenter(vars[0]),
1953         generated->GetYaxis()->GetTitle(), generated->GetYaxis()->GetBinCenter(vars[1]),
1954         generated->GetBinContent(vars[0], vars[1]));
1955       count++;
1956     }
1957     else if (generated->GetDimension() == 3 && generated->GetBinContent(vars[0], vars[1], vars[2]) < limit)
1958     {
1959       Printf("Empty bin at %s=%.2f %s=%.2f %s=%.2f (%.2f entries)",
1960         generated->GetXaxis()->GetTitle(), generated->GetXaxis()->GetBinCenter(vars[0]),
1961         generated->GetYaxis()->GetTitle(), generated->GetYaxis()->GetBinCenter(vars[1]),
1962         generated->GetZaxis()->GetTitle(), generated->GetZaxis()->GetBinCenter(vars[2]),
1963         generated->GetBinContent(vars[0], vars[1], vars[2]));
1964       count++;
1965     }
1966
1967     vars[2]++;
1968     if (vars[2] == binEnd[2]+1)
1969     {
1970       vars[2] = binBegin[2];
1971       vars[1]++;
1972     }
1973
1974     if (vars[1] == binEnd[1]+1)
1975     {
1976       vars[1] = binBegin[1];
1977       vars[0]++;
1978     }
1979
1980     if (vars[0] == binEnd[0]+1)
1981       break;
1982     total++;
1983   }
1984
1985   Printf("Correction has %d empty bins (out of %d bins)", count, total);
1986
1987   // rebin if required
1988   if (source == 2)
1989   {
1990     TAxis* axis = fEventHist->GetGrid(0)->GetGrid()->GetAxis(0);
1991
1992     if (axis->GetNbins() < measured->GetNbinsX())
1993     {
1994       if (axis2 != -1)
1995       {
1996         // 2d rebin with variable axis does not exist in root
1997
1998         TH1* tmp = measured;
1999         measured = new TH2D(Form("%s_rebinned", tmp->GetName()), tmp->GetTitle(), axis->GetNbins(), axis->GetXbins()->GetArray(), tmp->GetNbinsY(), tmp->GetYaxis()->GetXbins()->GetArray());
2000         for (Int_t x=1; x<=tmp->GetNbinsX(); x++)
2001           for (Int_t y=1; y<=tmp->GetNbinsY(); y++)
2002           {
2003             ((TH2*) measured)->Fill(tmp->GetXaxis()->GetBinCenter(x), tmp->GetYaxis()->GetBinCenter(y), tmp->GetBinContent(x, y));
2004             measured->SetBinError(x, y, 0); // cannot trust bin error, set to 0
2005           }
2006         delete tmp;
2007
2008         tmp = generated;
2009         generated = new TH2D(Form("%s_rebinned", tmp->GetName()), tmp->GetTitle(), axis->GetNbins(), axis->GetXbins()->GetArray(), tmp->GetNbinsY(), tmp->GetYaxis()->GetXbins()->GetArray());
2010         for (Int_t x=1; x<=tmp->GetNbinsX(); x++)
2011           for (Int_t y=1; y<=tmp->GetNbinsY(); y++)
2012           {
2013             ((TH2*) generated)->Fill(tmp->GetXaxis()->GetBinCenter(x), tmp->GetYaxis()->GetBinCenter(y), tmp->GetBinContent(x, y));
2014             generated->SetBinError(x, y, 0); // cannot trust bin error, set to 0
2015           }
2016         delete tmp;
2017       }
2018       else
2019       {
2020         TH1* tmp = measured;
2021         measured = tmp->Rebin(axis->GetNbins(), Form("%s_rebinned", tmp->GetName()), axis->GetXbins()->GetArray());
2022         delete tmp;
2023
2024         tmp = generated;
2025         generated = tmp->Rebin(axis->GetNbins(), Form("%s_rebinned", tmp->GetName()), axis->GetXbins()->GetArray());
2026         delete tmp;
2027       }
2028     }
2029     else if (axis->GetNbins() > measured->GetNbinsX())
2030     {
2031       if (axis2 != -1)
2032         AliFatal("Rebinning only works for 1d at present");
2033
2034       // this is an unfortunate case where the number of bins has to be increased in principle
2035       // there is a region where the binning is finner in one histogram and a region where it is the other way round
2036       // this cannot be resolved in principle, but as we only calculate the ratio the bin in the second region get the same entries
2037       // only a certain binning is supported here
2038       if (axis->GetNbins() != 100 || measured->GetNbinsX() != 39)
2039         AliFatal(Form("Invalid binning --> %d %d", axis->GetNbins(), measured->GetNbinsX()));
2040
2041       Double_t newBins[] = {0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 12.0, 14.0, 16.0, 18.0, 20.0, 25.0, 30.0, 35.0, 40.0, 45.0, 50.0, 100.0};
2042
2043       // reduce binning below 5 GeV/c
2044       TH1* tmp = measured;
2045       measured = tmp->Rebin(27, Form("%s_rebinned", tmp->GetName()), newBins);
2046       delete tmp;
2047
2048       // increase binning above 5 GeV/c
2049       tmp = measured;
2050       measured = new TH1F(Form("%s_rebinned2", tmp->GetName()), tmp->GetTitle(), axis->GetNbins(), axis->GetBinLowEdge(1), axis->GetBinUpEdge(axis->GetNbins()));
2051       for (Int_t bin=1; bin<=measured->GetNbinsX(); bin++)
2052       {
2053         measured->SetBinContent(bin, tmp->GetBinContent(tmp->FindBin(measured->GetBinCenter(bin))));
2054         measured->SetBinError(bin, tmp->GetBinError(tmp->FindBin(measured->GetBinCenter(bin))));
2055       }
2056       delete tmp;
2057
2058       // reduce binning below 5 GeV/c
2059       tmp = generated;
2060       generated = tmp->Rebin(27, Form("%s_rebinned", tmp->GetName()), newBins);
2061       delete tmp;
2062
2063       // increase binning above 5 GeV/c
2064       tmp = generated;
2065       generated = new TH1F(Form("%s_rebinned2", tmp->GetName()), tmp->GetTitle(), axis->GetNbins(), axis->GetBinLowEdge(1), axis->GetBinUpEdge(axis->GetNbins()));
2066       for (Int_t bin=1; bin<=generated->GetNbinsX(); bin++)
2067       {
2068         generated->SetBinContent(bin, tmp->GetBinContent(tmp->FindBin(generated->GetBinCenter(bin))));
2069         generated->SetBinError(bin, tmp->GetBinError(tmp->FindBin(generated->GetBinCenter(bin))));
2070       }
2071       delete tmp;
2072     }
2073   }
2074
2075   measured->Divide(measured, generated, 1, 1, "B");
2076
2077   delete generated;
2078
2079   ResetBinLimits(sourceContainer->GetGrid(step1));
2080   ResetBinLimits(sourceContainer->GetGrid(step2));
2081
2082   return measured;
2083 }
2084
2085 //____________________________________________________________________
2086 TH1* AliUEHist::GetEventEfficiency(CFStep step1, CFStep step2, Int_t axis1, Int_t axis2, Float_t ptLeadMin, Float_t ptLeadMax)
2087 {
2088   // creates a event-level efficiency by dividing step2 by step1
2089   // projected to axis1 and axis2 (optional if >= 0)
2090
2091   if (ptLeadMax > ptLeadMin)
2092   {
2093     fEventHist->GetGrid(step1)->SetRangeUser(0, ptLeadMin, ptLeadMax);
2094     fEventHist->GetGrid(step2)->SetRangeUser(0, ptLeadMin, ptLeadMax);
2095   }
2096
2097   TH1* measured = 0;
2098   TH1* generated = 0;
2099
2100   if (axis2 >= 0)
2101   {
2102     generated = fEventHist->Project(step1, axis1, axis2);
2103     measured = fEventHist->Project(step2, axis1, axis2);
2104   }
2105   else
2106   {
2107     generated = fEventHist->Project(step1, axis1);
2108     measured = fEventHist->Project(step2, axis1);
2109   }
2110
2111   measured->Divide(measured, generated, 1, 1, "B");
2112
2113   delete generated;
2114
2115   if (ptLeadMax > ptLeadMin)
2116   {
2117     fEventHist->GetGrid(step1)->SetRangeUser(0, 0, -1);
2118     fEventHist->GetGrid(step2)->SetRangeUser(0, 0, -1);
2119   }
2120
2121   return measured;
2122 }
2123
2124 //____________________________________________________________________
2125 void AliUEHist::WeightHistogram(TH3* hist1, TH1* hist2)
2126 {
2127   // weights each entry of the 3d histogram hist1 with the 1d histogram hist2
2128   // where the matching is done of the z axis of hist1 with the x axis of hist2
2129
2130   if (hist1->GetNbinsZ() != hist2->GetNbinsX())
2131     AliFatal(Form("Inconsistent binning %d %d", hist1->GetNbinsZ(), hist2->GetNbinsX()));
2132
2133   for (Int_t x=1; x<=hist1->GetNbinsX(); x++)
2134   {
2135     for (Int_t y=1; y<=hist1->GetNbinsY(); y++)
2136     {
2137       for (Int_t z=1; z<=hist1->GetNbinsZ(); z++)
2138       {
2139         if (hist2->GetBinContent(z) > 0)
2140         {
2141           hist1->SetBinContent(x, y, z, hist1->GetBinContent(x, y, z) / hist2->GetBinContent(z));
2142           hist1->SetBinError(x, y, z, hist1->GetBinError(x, y, z) / hist2->GetBinContent(z));
2143         }
2144         else
2145         {
2146           hist1->SetBinContent(x, y, z, 0);
2147           hist1->SetBinError(x, y, z, 0);
2148         }
2149       }
2150     }
2151   }
2152 }
2153
2154 //____________________________________________________________________
2155 TH1* AliUEHist::GetBias(CFStep step1, CFStep step2, Int_t region, const char* axis, Float_t leadPtMin, Float_t leadPtMax, Int_t weighting)
2156 {
2157   // extracts the track-level bias (integrating out the multiplicity) between two steps (dividing step2 by step1)
2158   // in the given region (sum over all regions is calculated if region == -1)
2159   // done by weighting the track-level distribution with the number of events as function of leading pT
2160   // and then calculating the ratio between the distributions
2161   // projected to axis which is a TH3::Project3D string, e.g. "x", or "yx"
2162   //   no projection is done if axis == 0
2163   // weighting: 0 = tracks weighted with events (as discussed above)
2164   //            1 = only track bias is returned
2165   //            2 = only event bias is returned
2166
2167   AliCFContainer* tmp = 0;
2168
2169   if (region == -1)
2170   {
2171     tmp = (AliCFContainer*) fTrackHist[0]->Clone();
2172     for (UInt_t i = 1; i < fkRegions; i++)
2173       if (fTrackHist[i])
2174         tmp->Add(fTrackHist[i]);
2175   }
2176   else if (region == kMin && fCombineMinMax)
2177   {
2178     tmp = (AliCFContainer*) fTrackHist[kMin]->Clone();
2179     tmp->Add(fTrackHist[kMax]);
2180   }
2181   else
2182     tmp = fTrackHist[region];
2183
2184   ResetBinLimits(tmp->GetGrid(step1));
2185   ResetBinLimits(fEventHist->GetGrid(step1));
2186   SetBinLimits(tmp->GetGrid(step1));
2187
2188   ResetBinLimits(tmp->GetGrid(step2));
2189   ResetBinLimits(fEventHist->GetGrid(step2));
2190   SetBinLimits(tmp->GetGrid(step2));
2191
2192   TH1D* events1 = (TH1D*)fEventHist->Project(step1, 0);
2193   TH3D* hist1 = (TH3D*)tmp->Project(step1, 0, tmp->GetNVar()-1, 2);
2194   if (weighting == 0)
2195     WeightHistogram(hist1, events1);
2196
2197   TH1D* events2 = (TH1D*)fEventHist->Project(step2, 0);
2198   TH3D* hist2 = (TH3D*)tmp->Project(step2, 0, tmp->GetNVar()-1, 2);
2199   if (weighting == 0)
2200     WeightHistogram(hist2, events2);
2201
2202   TH1* generated = hist1;
2203   TH1* measured = hist2;
2204
2205   if (weighting == 0 || weighting == 1)
2206   {
2207     if (axis)
2208     {
2209       if (leadPtMax > leadPtMin)
2210       {
2211         hist1->GetZaxis()->SetRangeUser(leadPtMin, leadPtMax);
2212         hist2->GetZaxis()->SetRangeUser(leadPtMin, leadPtMax);
2213       }
2214
2215       if (fEtaMax > fEtaMin && !TString(axis).Contains("x"))
2216       {
2217         hist1->GetXaxis()->SetRangeUser(fEtaMin, fEtaMax);
2218         hist2->GetXaxis()->SetRangeUser(fEtaMin, fEtaMax);
2219       }
2220
2221       generated = hist1->Project3D(axis);
2222       measured  = hist2->Project3D(axis);
2223
2224       // delete hists here if projection has been done
2225       delete hist1;
2226       delete hist2;
2227     }
2228     delete events1;
2229     delete events2;
2230   }
2231   else if (weighting == 2)
2232   {
2233     delete hist1;
2234     delete hist2;
2235     generated = events1;
2236     measured = events2;
2237   }
2238
2239   measured->Divide(generated);
2240
2241   delete generated;
2242
2243   ResetBinLimits(tmp->GetGrid(step1));
2244   ResetBinLimits(tmp->GetGrid(step2));
2245
2246   if ((region == -1) || (region == kMin && fCombineMinMax))
2247     delete tmp;
2248
2249   return measured;
2250 }
2251
2252 //____________________________________________________________________
2253 void AliUEHist::CorrectCorrelatedContamination(CFStep step, Int_t region, TH1* trackCorrection)
2254 {
2255   // corrects for the given factor in a small delta-eta and delta-phi window as function of pT,A and pT,T
2256
2257   if (!fTrackHist[region])
2258     return;
2259
2260   THnSparse* grid = fTrackHist[region]->GetGrid(step)->GetGrid();
2261
2262   Int_t var1 = 1;
2263   Int_t var2 = 2;
2264
2265   if (grid->GetAxis(var1)->GetNbins() != trackCorrection->GetNbinsX())
2266     AliFatal(Form("Invalid binning (var1): %d %d", grid->GetAxis(var1)->GetNbins(), trackCorrection->GetNbinsX()));
2267
2268   if (grid->GetAxis(var2)->GetNbins() != trackCorrection->GetNbinsY())
2269     AliFatal(Form("Invalid binning (var2): %d %d", grid->GetAxis(var2)->GetNbins(), trackCorrection->GetNbinsY()));
2270
2271   // optimized implementation
2272   for (Int_t binIdx = 0; binIdx < grid->GetNbins(); binIdx++)
2273   {
2274     Int_t bins[6];
2275
2276     Double_t value = grid->GetBinContent(binIdx, bins);
2277     Double_t error = grid->GetBinError(binIdx);
2278
2279     // check eta and phi axes
2280     if (TMath::Abs(grid->GetAxis(0)->GetBinCenter(bins[0])) > 0.1)
2281       continue;
2282     if (TMath::Abs(grid->GetAxis(4)->GetBinCenter(bins[4])) > 0.1)
2283       continue;
2284
2285     value *= trackCorrection->GetBinContent(bins[var1], bins[var2]);
2286     error *= trackCorrection->GetBinContent(bins[var1], bins[var2]);
2287
2288     grid->SetBinContent(bins, value);
2289     grid->SetBinError(bins, error);
2290   }
2291
2292   Printf("AliUEHist::CorrectCorrelatedContamination: Corrected.");
2293 }
2294
2295 //____________________________________________________________________
2296 TH2* AliUEHist::GetCorrelatedContamination()
2297 {
2298   // contamination correlated with the trigger particle is evaluated between step kCFStepTracked and kCFStepTrackedOnlyPrim in the region of delta eta and delta phi < 0.1 (smallest bin!)
2299
2300   Int_t step1 = kCFStepTrackedOnlyPrim;
2301   Int_t step2 = kCFStepTracked;
2302
2303   fTrackHist[0]->GetGrid(step1)->SetRangeUser(0, -0.01, 0.01); // delta eta
2304   fTrackHist[0]->GetGrid(step1)->SetRangeUser(4, -0.01, 0.01); // delta phi
2305   TH2* tracksStep1 = (TH2*) fTrackHist[0]->Project(step1, 1, 2);
2306
2307   fTrackHist[0]->GetGrid(step2)->SetRangeUser(0, -0.01, 0.01); // delta eta
2308   fTrackHist[0]->GetGrid(step2)->SetRangeUser(4, -0.01, 0.01); // delta phi
2309   TH2* tracksStep2 = (TH2*) fTrackHist[0]->Project(step2, 1, 2);
2310
2311   tracksStep1->Divide(tracksStep2);
2312
2313   TH1* triggersStep1 = fEventHist->Project(step1, 0);
2314   TH1* triggersStep2 = fEventHist->Project(step2, 0);
2315
2316   TH1* singleParticle = GetTrackingContamination(1);
2317
2318   for (Int_t x=1; x<=tracksStep1->GetNbinsX(); x++)
2319     for (Int_t y=1; y<=tracksStep1->GetNbinsY(); y++)
2320       if (singleParticle->GetBinContent(x) > 0 && triggersStep1->GetBinContent(y) > 0)
2321         tracksStep1->SetBinContent(x, y, tracksStep1->GetBinContent(x, y) / triggersStep1->GetBinContent(y) * triggersStep2->GetBinContent(y) / singleParticle->GetBinContent(x));
2322       else
2323         tracksStep1->SetBinContent(x, y, 0);
2324
2325   delete singleParticle;
2326   delete tracksStep2;
2327   delete triggersStep1;
2328   delete triggersStep2;
2329
2330   return tracksStep1;
2331 }
2332
2333 //____________________________________________________________________
2334 TH2D* AliUEHist::GetTrackingEfficiency()
2335 {
2336   // extracts the tracking efficiency by calculating the efficiency from step kCFStepAnaTopology to kCFStepTrackedOnlyPrim
2337   // integrates over the regions and all other variables than pT and eta to increase the statistics
2338   //
2339   // returned histogram has to be deleted by the user
2340
2341   return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepAnaTopology, kCFStepTrackedOnlyPrim, 0, 1));
2342 }
2343
2344 //____________________________________________________________________
2345 TH2D* AliUEHist::GetFakeRate()
2346 {
2347   return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepTracked, (CFStep) (kCFStepTracked+3), 0, 1));
2348 }
2349
2350 //____________________________________________________________________
2351 TH2D* AliUEHist::GetTrackingEfficiencyCentrality()
2352 {
2353   // extracts the tracking efficiency by calculating the efficiency from step kCFStepAnaTopology to kCFStepTrackedOnlyPrim
2354   // integrates over the regions and all other variables than pT, centrality to increase the statistics
2355   //
2356   // returned histogram has to be deleted by the user
2357
2358   return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepAnaTopology, kCFStepTrackedOnlyPrim, 1, 3));
2359 }
2360
2361 //____________________________________________________________________
2362 TH1D* AliUEHist::GetTrackingEfficiency(Int_t axis)
2363 {
2364   // extracts the tracking efficiency by calculating the efficiency from step kCFStepAnaTopology to kCFStepTrackedOnlyPrim
2365   // integrates over the regions and all other variables than pT (axis == 0) and eta (axis == 1) to increase the statistics
2366
2367   return dynamic_cast<TH1D*> (GetTrackEfficiency(kCFStepAnaTopology, kCFStepTrackedOnlyPrim, axis));
2368 }
2369
2370 //____________________________________________________________________
2371 TH1D* AliUEHist::GetFakeRate(Int_t axis)
2372 {
2373   return dynamic_cast<TH1D*> (GetTrackEfficiency(kCFStepTracked, (CFStep) (kCFStepTracked+3), axis));
2374 }
2375 //____________________________________________________________________
2376 TH2D* AliUEHist::GetTrackingCorrection()
2377 {
2378   // extracts the tracking correction by calculating the efficiency from step kCFStepAnaTopology to kCFStepTracked
2379   // integrates over the regions and all other variables than pT and eta to increase the statistics
2380   //
2381   // returned histogram has to be deleted by the user
2382
2383   return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepTracked, kCFStepAnaTopology, 0, 1));
2384 }
2385
2386 //____________________________________________________________________
2387 TH1D* AliUEHist::GetTrackingCorrection(Int_t axis)
2388 {
2389   // extracts the tracking correction by calculating the efficiency from step kCFStepAnaTopology to kCFStepTracked
2390   // integrates over the regions and all other variables than pT (axis == 0) and eta (axis == 1) to increase the statistics
2391
2392   return dynamic_cast<TH1D*> (GetTrackEfficiency(kCFStepTracked, kCFStepAnaTopology, axis));
2393 }
2394
2395 //____________________________________________________________________
2396 TH2D* AliUEHist::GetTrackingEfficiencyCorrection()
2397 {
2398   // extracts the tracking correction by calculating the efficiency from step kCFStepAnaTopology to kCFStepTracked
2399   // integrates over the regions and all other variables than pT and eta to increase the statistics
2400   //
2401   // returned histogram has to be deleted by the user
2402
2403   return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, 0, 1));
2404 }
2405
2406 //____________________________________________________________________
2407 TH2D* AliUEHist::GetTrackingEfficiencyCorrectionCentrality()
2408 {
2409   // extracts the tracking correction by calculating the efficiency from step kCFStepAnaTopology to kCFStepTracked
2410   // integrates over the regions and all other variables than pT and centrality to increase the statistics
2411   //
2412   // returned histogram has to be deleted by the user
2413
2414   return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, 1, 3));
2415 }
2416
2417 //____________________________________________________________________
2418 TH1D* AliUEHist::GetTrackingEfficiencyCorrection(Int_t axis)
2419 {
2420   // extracts the tracking correction by calculating the efficiency from step kCFStepAnaTopology to kCFStepTracked
2421   // integrates over the regions and all other variables than pT (axis == 0) and eta (axis == 1) to increase the statistics
2422
2423   return dynamic_cast<TH1D*> (GetTrackEfficiency(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, axis));
2424 }
2425
2426 //____________________________________________________________________
2427 TH2D* AliUEHist::GetTrackingContamination()
2428 {
2429   // extracts the tracking contamination by secondaries by calculating the efficiency from step kCFStepTrackedOnlyPrim to kCFStepTracked
2430   // integrates over the regions and all other variables than pT and eta to increase the statistics
2431   //
2432   // returned histogram has to be deleted by the user
2433
2434   return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepTracked, kCFStepTrackedOnlyPrim, 0, 1));
2435 }
2436
2437 //____________________________________________________________________
2438 TH2D* AliUEHist::GetTrackingContaminationCentrality()
2439 {
2440   // extracts the tracking contamination by secondaries by calculating the efficiency from step kCFStepTrackedOnlyPrim to kCFStepTracked
2441   // integrates over the regions and all other variables than pT and centrality to increase the statistics
2442   //
2443   // returned histogram has to be deleted by the user
2444
2445   return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepTracked, kCFStepTrackedOnlyPrim, 1, 3));
2446 }
2447
2448 //____________________________________________________________________
2449 TH1D* AliUEHist::GetTrackingContamination(Int_t axis)
2450 {
2451   // extracts the tracking contamination by secondaries by calculating the efficiency from step kCFStepTrackedOnlyPrim to kCFStepTracked
2452   // integrates over the regions and all other variables than pT (axis == 0) and eta (axis == 1) to increase the statistics
2453
2454   return dynamic_cast<TH1D*> (GetTrackEfficiency(kCFStepTracked, kCFStepTrackedOnlyPrim, axis));
2455 }
2456
2457 //____________________________________________________________________
2458 const char* AliUEHist::GetRegionTitle(Region region)
2459 {
2460   // returns the name of the given region
2461
2462   switch (region)
2463   {
2464     case kToward:
2465       return "Towards";
2466     case kAway:
2467       return "Away";
2468     case kMin:
2469       return (fCombineMinMax) ? "Transverse" : "Min";
2470     case kMax:
2471       return "Max";
2472   }
2473
2474   return 0;
2475 }
2476
2477 //____________________________________________________________________
2478 const char* AliUEHist::GetStepTitle(CFStep step)
2479 {
2480   // returns the name of the given step
2481
2482   switch (step)
2483   {
2484     case kCFStepAll:
2485       return "All events";
2486     case kCFStepTriggered:
2487       return "Triggered";
2488     case kCFStepVertex:
2489       return "Primary Vertex";
2490     case kCFStepAnaTopology:
2491       return "Required analysis topology";
2492     case kCFStepTrackedOnlyPrim:
2493       return "Tracked (matched MC, only primaries)";
2494     case kCFStepTracked:
2495       return "Tracked (matched MC, all)";
2496     case kCFStepReconstructed:
2497       return "Reconstructed";
2498     case kCFStepRealLeading:
2499       return "Correct leading particle identified";
2500     case kCFStepBiasStudy:
2501       return "Bias study applying tracking efficiency";
2502     case kCFStepBiasStudy2:
2503       return "Bias study applying tracking efficiency in two steps";
2504     case kCFStepCorrected:
2505       return "Corrected for efficiency on-the-fly";
2506   }
2507
2508   return 0;
2509 }
2510
2511 //____________________________________________________________________
2512 void AliUEHist::CopyReconstructedData(AliUEHist* from)
2513 {
2514   // copies those histograms extracted from ESD to this object
2515
2516   // TODO at present only the pointers are copied
2517
2518   for (Int_t region=0; region<4; region++)
2519   {
2520     if (!fTrackHist[region])
2521       continue;
2522
2523     fTrackHist[region]->SetGrid(AliUEHist::kCFStepReconstructed, from->fTrackHist[region]->GetGrid(AliUEHist::kCFStepReconstructed));
2524     //fTrackHist[region]->SetGrid(AliUEHist::kCFStepTrackedOnlyPrim, from->fTrackHist[region]->GetGrid(AliUEHist::kCFStepTrackedOnlyPrim));
2525     fTrackHist[region]->SetGrid(AliUEHist::kCFStepBiasStudy,     from->fTrackHist[region]->GetGrid(AliUEHist::kCFStepBiasStudy));
2526   }
2527
2528   fEventHist->SetGrid(AliUEHist::kCFStepReconstructed, from->fEventHist->GetGrid(AliUEHist::kCFStepReconstructed));
2529   //fEventHist->SetGrid(AliUEHist::kCFStepTrackedOnlyPrim, from->fEventHist->GetGrid(AliUEHist::kCFStepTrackedOnlyPrim));
2530   fEventHist->SetGrid(AliUEHist::kCFStepBiasStudy,     from->fEventHist->GetGrid(AliUEHist::kCFStepBiasStudy));
2531 }
2532
2533 void AliUEHist::DeepCopy(AliUEHist* from)
2534 {
2535   // copies the entries of this object's members from the object <from> to this object
2536   // fills using the fill function and thus allows that the objects have different binning
2537
2538   for (Int_t region=0; region<4; region++)
2539   {
2540     if (!fTrackHist[region] || !from->fTrackHist[region])
2541       continue;
2542
2543     for (Int_t step=0; step<fTrackHist[region]->GetNStep(); step++)
2544     {
2545       Printf("Copying region %d step %d", region, step);
2546       THnSparse* target = fTrackHist[region]->GetGrid(step)->GetGrid();
2547       THnSparse* source = from->fTrackHist[region]->GetGrid(step)->GetGrid();
2548
2549       target->Reset();
2550       target->RebinnedAdd(source);
2551     }
2552   }
2553
2554   for (Int_t step=0; step<fEventHist->GetNStep(); step++)
2555   {
2556     Printf("Ev: Copying step %d", step);
2557     THnSparse* target = fEventHist->GetGrid(step)->GetGrid();
2558     THnSparse* source = from->fEventHist->GetGrid(step)->GetGrid();
2559
2560     target->Reset();
2561     target->RebinnedAdd(source);
2562   }
2563
2564   for (Int_t step=0; step<TMath::Min(fTrackHistEfficiency->GetNStep(), from->fTrackHistEfficiency->GetNStep()); step++)
2565   {
2566     if (!from->fTrackHistEfficiency->GetGrid(step))
2567       continue;
2568
2569     Printf("Eff: Copying step %d", step);
2570     THnSparse* target = fTrackHistEfficiency->GetGrid(step)->GetGrid();
2571     THnSparse* source = from->fTrackHistEfficiency->GetGrid(step)->GetGrid();
2572
2573     target->Reset();
2574     target->RebinnedAdd(source);
2575   }
2576 }
2577
2578 void AliUEHist::SymmetrizepTBins()
2579 {
2580   // copy pt,a < pt,t bins to pt,a > pt,t (inverting deltaphi and delta eta as it should be) including symmetric bins
2581
2582   for (Int_t region=0; region<4; region++)
2583   {
2584     if (!fTrackHist[region])
2585       continue;
2586
2587     for (Int_t step=0; step<fTrackHist[region]->GetNStep(); step++)
2588     {
2589       Printf("Copying region %d step %d", region, step);
2590       THnSparse* target = fTrackHist[region]->GetGrid(step)->GetGrid();
2591       if (target->GetEntries() == 0)
2592         continue;
2593
2594       // for symmetric bins
2595       THnSparse* source = (THnSparse*) target->Clone();
2596
2597       // axes: 0 delta eta; 1 pT,a; 2 pT,t; 3 centrality; 4 delta phi; 5 vtx-z
2598       for (Int_t i3 = 1; i3 <= target->GetAxis(3)->GetNbins(); i3++)
2599         for (Int_t i5 = 1; i5 <= target->GetAxis(5)->GetNbins(); i5++)
2600         {
2601           for (Int_t i1 = 1; i1 <= target->GetAxis(1)->GetNbins(); i1++)
2602             for (Int_t i2 = 1; i2 <= target->GetAxis(2)->GetNbins(); i2++)
2603             {
2604               // find source bin
2605               Int_t binA = target->GetAxis(1)->FindBin(target->GetAxis(2)->GetBinCenter(i2));
2606               Int_t binT = target->GetAxis(2)->FindBin(target->GetAxis(1)->GetBinCenter(i1));
2607
2608               Printf("(%d %d) Copying from %d %d to %d %d", i3, i5, binA, binT, i1, i2);
2609
2610               for (Int_t i0 = 1; i0 <= target->GetAxis(0)->GetNbins(); i0++)
2611                 for (Int_t i4 = 1; i4 <= target->GetAxis(4)->GetNbins(); i4++)
2612                 {
2613                   Int_t binEta = target->GetAxis(0)->FindBin(-target->GetAxis(0)->GetBinCenter(i0));
2614                   Double_t phi = -target->GetAxis(4)->GetBinCenter(i4);
2615                   if (phi < -TMath::Pi()/2)
2616                     phi += TMath::TwoPi();
2617                   Int_t binPhi = target->GetAxis(4)->FindBin(phi);
2618
2619                   Int_t binSource[] = { binEta, binA, binT, i3, binPhi, i5 };
2620                   Int_t binTarget[] = { i0, i1, i2, i3, i4, i5 };
2621
2622                   Double_t value = source->GetBinContent(binSource);
2623                   Double_t error = source->GetBinError(binSource);
2624
2625                   if (error == 0)
2626                     continue;
2627
2628                   Double_t value2 = target->GetBinContent(binTarget);
2629                   Double_t error2 = target->GetBinError(binTarget);
2630
2631                   Double_t sum = value;
2632                   Double_t err = error;
2633
2634                   if (error2 > 0)
2635                   {
2636                     sum = value + value2;
2637                     err = TMath::Sqrt(error * error + error2 * error2);
2638                   }
2639
2640                   // Printf("  Values: %f +- %f; %f +- %f --> %f +- %f", value, error, value2, error2, sum, err);
2641
2642                   target->SetBinContent(binTarget, sum);
2643                   target->SetBinError(binTarget, err);
2644                 }
2645             }
2646         }
2647
2648         delete source;
2649     }
2650   }
2651 }
2652
2653 //____________________________________________________________________
2654 void AliUEHist::ExtendTrackingEfficiency(Bool_t verbose)
2655 {
2656   // fits the tracking efficiency at high pT with a constant and fills all bins with this tracking efficiency
2657
2658   Float_t fitRangeBegin = 5.01;
2659   Float_t fitRangeEnd = 14.99;
2660   Float_t extendRangeBegin = 10.01;
2661
2662   if (fTrackHistEfficiency->GetNVar() == 3)
2663   {
2664     TH1* obj = GetTrackingEfficiency(1);
2665
2666     if (verbose)
2667     {
2668       new TCanvas;
2669       obj->Draw();
2670     }
2671
2672     obj->Fit("pol0", (verbose) ? "+" : "0+", "SAME", fitRangeBegin, fitRangeEnd);
2673
2674     Float_t trackingEff = obj->GetFunction("pol0")->GetParameter(0);
2675
2676     obj = GetTrackingContamination(1);
2677
2678     if (verbose)
2679     {
2680       new TCanvas;
2681       obj->Draw();
2682     }
2683
2684     obj->Fit("pol0", (verbose) ? "+" : "0+", "SAME", fitRangeBegin, fitRangeEnd);
2685
2686     Float_t trackingCont = obj->GetFunction("pol0")->GetParameter(0);
2687
2688     Printf("AliUEHist::ExtendTrackingEfficiency: Fitted efficiency between %f and %f and got %f tracking efficiency and %f tracking contamination correction. Extending from %f onwards (within %f < eta < %f)", fitRangeBegin, fitRangeEnd, trackingEff, trackingCont, extendRangeBegin, fEtaMin, fEtaMax);
2689
2690     // extend for full pT range
2691     for (Int_t x = fTrackHistEfficiency->GetAxis(0, 0)->FindBin(fEtaMin); x <= fTrackHistEfficiency->GetAxis(0, 0)->FindBin(fEtaMax); x++)
2692       for (Int_t y = fTrackHistEfficiency->GetAxis(1, 0)->FindBin(extendRangeBegin); y <= fTrackHistEfficiency->GetNBins(1); y++)
2693         for (Int_t z = 1; z <= fTrackHistEfficiency->GetNBins(2); z++) // particle type axis
2694         {
2695
2696           Int_t bins[3];
2697           bins[0] = x;
2698           bins[1] = y;
2699           bins[2] = z;
2700
2701           fTrackHistEfficiency->GetGrid(0)->SetElement(bins, 100);
2702           fTrackHistEfficiency->GetGrid(1)->SetElement(bins, 100.0 * trackingEff);
2703           fTrackHistEfficiency->GetGrid(2)->SetElement(bins, 100.0 * trackingEff / trackingCont);
2704         }
2705   }
2706   else if (fTrackHistEfficiency->GetNVar() == 4)
2707   {
2708     // fit in centrality intervals of 20% for efficiency, one bin for contamination
2709     Float_t* trackingEff = 0;
2710     Float_t* trackingCont = 0;
2711     Float_t centralityBins[] = { 0, 10, 20, 40, 60, 100 };
2712     Int_t nCentralityBins = 5;
2713
2714     Printf("AliUEHist::ExtendTrackingEfficiency: Fitting efficiencies between %f and %f. Extending from %f onwards (within %f < eta < %f)", fitRangeBegin, fitRangeEnd, extendRangeBegin, fEtaMin, fEtaMax);
2715
2716     // 0 = eff; 1 = cont
2717     for (Int_t caseNo = 0; caseNo < 2; caseNo++)
2718     {
2719       Float_t* target = 0;
2720       Int_t centralityBinsLocal = nCentralityBins;
2721
2722       if (caseNo == 0)
2723       {
2724         trackingEff = new Float_t[centralityBinsLocal];
2725         target = trackingEff;
2726       }
2727       else
2728       {
2729         centralityBinsLocal = 1;
2730         trackingCont = new Float_t[centralityBinsLocal];
2731         target = trackingCont;
2732       }
2733
2734       for (Int_t i=0; i<centralityBinsLocal; i++)
2735       {
2736         if (centralityBinsLocal == 1)
2737           SetCentralityRange(centralityBins[0] + 0.1, centralityBins[nCentralityBins] - 0.1);
2738         else
2739           SetCentralityRange(centralityBins[i] + 0.1, centralityBins[i+1] - 0.1);
2740         TH1* proj = (caseNo == 0) ? GetTrackingEfficiency(1) : GetTrackingContamination(1);
2741         if (verbose)
2742         {
2743           new TCanvas;
2744           proj->DrawCopy();
2745         }
2746         if ((Int_t) proj->Fit("pol0", (verbose) ? "+" : "Q0+", "SAME", fitRangeBegin, fitRangeEnd) == 0)
2747           target[i] = proj->GetFunction("pol0")->GetParameter(0);
2748         else
2749           target[i] = 0;
2750         Printf("AliUEHist::ExtendTrackingEfficiency: case %d, centrality %d, eff %f", caseNo, i, target[i]);
2751       }
2752     }
2753
2754     // extend for full pT range
2755     for (Int_t x = fTrackHistEfficiency->GetAxis(0, 0)->FindBin(fEtaMin); x <= fTrackHistEfficiency->GetAxis(0, 0)->FindBin(fEtaMax); x++)
2756       for (Int_t y = fTrackHistEfficiency->GetAxis(1, 0)->FindBin(extendRangeBegin); y <= fTrackHistEfficiency->GetNBins(1); y++)
2757         for (Int_t z = 1; z <= fTrackHistEfficiency->GetNBins(2); z++) // particle type axis
2758         {
2759           for (Int_t z2 = 1; z2 <= fTrackHistEfficiency->GetNBins(3); z2++) // centrality axis
2760           {
2761
2762             Int_t bins[4];
2763             bins[0] = x;
2764             bins[1] = y;
2765             bins[2] = z;
2766             bins[3] = z2;
2767
2768             Int_t z2Bin = 0;
2769             while (centralityBins[z2Bin+1] < fTrackHistEfficiency->GetAxis(3, 0)->GetBinCenter(z2))
2770               z2Bin++;
2771
2772             //Printf("%d %d", z2, z2Bin);
2773
2774             fTrackHistEfficiency->GetGrid(0)->SetElement(bins, 100);
2775             fTrackHistEfficiency->GetGrid(1)->SetElement(bins, 100.0 * trackingEff[z2Bin]);
2776             if (trackingCont[0] > 0)
2777               fTrackHistEfficiency->GetGrid(2)->SetElement(bins, 100.0 * trackingEff[z2Bin] / trackingCont[0]);
2778             else
2779               fTrackHistEfficiency->GetGrid(2)->SetElement(bins, 0);
2780           }
2781        }
2782
2783      delete[] trackingEff;
2784      delete[] trackingCont;
2785    }
2786
2787    SetCentralityRange(0, 0);
2788 }
2789
2790 void AliUEHist::AdditionalDPhiCorrection(Int_t step)
2791 {
2792   // corrects the dphi distribution with an extra factor close to dphi ~ 0
2793
2794   Printf("WARNING: In AliUEHist::AdditionalDPhiCorrection.");
2795
2796   THnSparse* grid = fTrackHist[0]->GetGrid(step)->GetGrid();
2797
2798   // optimized implementation
2799   for (Int_t binIdx = 0; binIdx < grid->GetNbins(); binIdx++)
2800   {
2801     Int_t bins[6];
2802     Double_t value = grid->GetBinContent(binIdx, bins);
2803     Double_t error = grid->GetBinError(binIdx);
2804
2805     Float_t binCenter = grid->GetAxis(4)->GetBinCenter(bins[4]);
2806     if (TMath::Abs(binCenter) < 0.2)
2807     {
2808       value *= 0.985;
2809       error *= 0.985;
2810     }
2811     else if (TMath::Abs(binCenter) < 0.3)
2812     {
2813       value *= 0.9925;
2814       error *= 0.9925;
2815     }
2816
2817     grid->SetBinContent(bins, value);
2818     grid->SetBinError(bins, error);
2819   }
2820 }
2821
2822 void AliUEHist::Scale(Double_t factor)
2823 {
2824   // scales all contained histograms by the given factor
2825
2826   for (Int_t i=0; i<4; i++)
2827     if (fTrackHist[i])
2828       fTrackHist[i]->Scale(factor);
2829
2830   fEventHist->Scale(factor);
2831   fTrackHistEfficiency->Scale(factor);
2832 }
2833
2834 void AliUEHist::Reset()
2835 {
2836   // resets all contained histograms
2837
2838   for (Int_t i=0; i<4; i++)
2839     if (fTrackHist[i])
2840       for (Int_t step=0; step<fTrackHist[i]->GetNStep(); step++)
2841         fTrackHist[i]->GetGrid(step)->GetGrid()->Reset();
2842
2843   for (Int_t step=0; step<fEventHist->GetNStep(); step++)
2844     fEventHist->GetGrid(step)->GetGrid()->Reset();
2845
2846   for (Int_t step=0; step<fTrackHistEfficiency->GetNStep(); step++)
2847     fTrackHistEfficiency->GetGrid(step)->GetGrid()->Reset();
2848 }
2849
2850 THnBase* AliUEHist::ChangeToThn(THnBase* sparse)
2851 {
2852   // change the object to THn for faster processing
2853
2854         // convert to THn (SEGV's for some strange reason...)
2855         // x = THn::CreateHn("a", "a", sparse);
2856
2857   // own implementation
2858   Int_t nBins[10];
2859   for (Int_t i=0; i<sparse->GetNdimensions(); i++)
2860     nBins[i] = sparse->GetAxis(i)->GetNbins();
2861   THn* tmpTHn = new THnF(Form("%s_thn", sparse->GetName()), sparse->GetTitle(), sparse->GetNdimensions(), nBins, 0, 0);
2862   for (Int_t i=0; i<sparse->GetNdimensions(); i++)
2863   {
2864     tmpTHn->SetBinEdges(i, sparse->GetAxis(i)->GetXbins()->GetArray());
2865     tmpTHn->GetAxis(i)->SetTitle(sparse->GetAxis(i)->GetTitle());
2866   }
2867   tmpTHn->RebinnedAdd(sparse);
2868
2869   return tmpTHn;
2870 }
2871
2872 void AliUEHist::CondenseBin(THnSparse* grid, THnSparse* target, Int_t axis, Float_t targetValue, Float_t from, Float_t to)
2873 {
2874   //
2875   // loops through the histogram and moves all entries to a single point <targetValue> on the axis <axis>
2876   // if <from> and <to> are set, then moving only occurs if value on <axis> is betweem <from> and <to>
2877   //
2878
2879   if (grid->GetNdimensions() > 6)
2880     AliFatal("Too many dimensions in THnSparse");
2881
2882   Int_t targetBin = grid->GetAxis(axis)->FindBin(targetValue);
2883   AliInfo(Form("Target bin on axis %d with value %f is %d", axis, targetValue, targetBin));
2884
2885   Int_t fromBin = 1;
2886   Int_t toBin = grid->GetAxis(axis)->GetNbins();
2887   if (to > from)
2888   {
2889     fromBin = grid->GetAxis(axis)->FindBin(from);
2890     toBin = grid->GetAxis(axis)->FindBin(to);
2891     AliInfo(Form("Only condensing between bin %d and %d", fromBin, toBin));
2892   }
2893
2894   Int_t bins[6];
2895   for (Int_t binIdx = 0; binIdx < grid->GetNbins(); binIdx++)
2896   {
2897     Double_t value = grid->GetBinContent(binIdx, bins);
2898     Double_t error = grid->GetBinError(binIdx);
2899
2900     if (bins[axis] >= fromBin && bins[axis] <= toBin)
2901       bins[axis] = targetBin;
2902
2903     value += target->GetBinContent(bins);
2904     error = TMath::Sqrt(error * error + target->GetBinError(bins) * target->GetBinError(bins));
2905
2906     target->SetBinContent(bins, value);
2907     target->SetBinError(bins, error);
2908   }
2909 }
2910
2911 void AliUEHist::CondenseBin(CFStep step, Int_t trackAxis, Int_t eventAxis, Float_t targetValue, Float_t from, Float_t to, CFStep tmpStep)
2912 {
2913   // loops through the histogram at <step> and moves all entries to a single point <targetValue> on the axes
2914   // <trackAxis> and <eventAxis>. <tmpStep> is used to temporary store the data
2915   // This is useful e.g. to move bin content around for MC productions where the centrality selection did
2916   // not yield the desired result
2917
2918   // reset tmpStep
2919   fEventHist->GetGrid(tmpStep)->GetGrid()->Reset();
2920   for (UInt_t i=0; i<fkRegions; i++)
2921     if (fTrackHist[i])
2922       fTrackHist[i]->GetGrid(tmpStep)->GetGrid()->Reset();
2923
2924   // copy to tmpStep
2925   CorrectTracks(step, tmpStep, 0, -1);
2926   CorrectEvents(step, tmpStep, 0, -1);
2927
2928   // reset step
2929   fEventHist->GetGrid(step)->GetGrid()->Reset();
2930   for (UInt_t i=0; i<fkRegions; i++)
2931     if (fTrackHist[i])
2932       fTrackHist[i]->GetGrid(step)->GetGrid()->Reset();
2933
2934   // rewriting
2935   for (UInt_t i=0; i<fkRegions; i++)
2936   {
2937     if (!fTrackHist[i])
2938       continue;
2939
2940     THnSparse* grid = fTrackHist[i]->GetGrid(tmpStep)->GetGrid();
2941     THnSparse* target = fTrackHist[i]->GetGrid(step)->GetGrid();
2942
2943     CondenseBin(grid, target, trackAxis, targetValue, from, to);
2944   }
2945   CondenseBin(fEventHist->GetGrid(tmpStep)->GetGrid(), fEventHist->GetGrid(step)->GetGrid(), eventAxis, targetValue, from, to);
2946
2947   // reset tmpStep
2948   fEventHist->GetGrid(tmpStep)->GetGrid()->Reset();
2949   for (UInt_t i=0; i<fkRegions; i++)
2950     if (fTrackHist[i])
2951       fTrackHist[i]->GetGrid(tmpStep)->GetGrid()->Reset();
2952 }