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