PWG4/JetTasks/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
  38 ClassImp(AliUEHist)
  39
  40 const Int_t AliUEHist::fgkCFSteps = 10;
  41
  42 AliUEHist::AliUEHist(const char* reqHist) :
  43   TObject(),
  44   fkRegions(4),
  45   fEventHist(0),
  46   fTrackHistEfficiency(0),
  47   fEtaMin(0),
  48   fEtaMax(0),
  49   fPtMin(0),
  50   fPtMax(0),
  51   fCentralityMin(0),
  52   fCentralityMax(0),
  53   fContaminationEnhancement(0),
  54   fCombineMinMax(0),
  55   fCache(0),
  56   fHistogramType(reqHist)
  57 {
  58   // Constructor
  59
  60   for (Int_t i=0; i<fkRegions; i++)
  61     fTrackHist[i] = 0;
  62
  63   if (strlen(reqHist) == 0)
  64     return;
  65
  66   AliLog::SetClassDebugLevel("AliCFContainer", -1);
  67   AliLog::SetClassDebugLevel("AliCFGridSparse", -3);
  68
  69   const char* title = "";
  70
  71   // track level
  72   Int_t nTrackVars = 4; // eta vs pT vs pT,lead (vs delta phi) vs multiplicity
  73   Int_t iTrackBin[5];
  74   Double_t* trackBins[5];
  75   const char* trackAxisTitle[5];
  76
  77   // eta
  78   iTrackBin[0] = 20;
  79   Double_t etaBins[20+1];
  80   for (Int_t i=0; i<=iTrackBin[0]; i++)
  81     etaBins[i] = -1.0 + 0.1 * i;
  82   trackBins[0] = etaBins;
  83   trackAxisTitle[0] = "#eta";
  84
  85   // delta eta
  86   const Int_t kNDeltaEtaBins = 32;
  87   Double_t deltaEtaBins[kNDeltaEtaBins+1];
  88   for (Int_t i=0; i<=kNDeltaEtaBins; i++)
  89     deltaEtaBins[i] = -1.6 + 0.1 * i;
  90
  91   // pT
  92   iTrackBin[1] = 39;
  93   Double_t pTBins[] = {0.0, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.6, 0.7, 0.8, 0.9, 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};
  94   trackBins[1] = pTBins;
  95   trackAxisTitle[1] = "p_{T} (GeV/c)";
  96
  97   // pT,lead binning 1
  98   const Int_t kNLeadingpTBins = 100;
  99   Double_t leadingpTBins[kNLeadingpTBins+1];
 100   for (Int_t i=0; i<=kNLeadingpTBins; i++)
 101     leadingpTBins[i] = 0.5 * i;
 102
 103   // pT,lead binning 2
 104   const Int_t kNLeadingpTBins2 = 14;
 105   Double_t leadingpTBins2[] = { 0.0, 0.5, 1.0, 2.0, 3.0, 4.0, 6.0, 8.0, 10.0, 15.0, 20.0, 30.0, 40.0, 50.0, 100.0 };
 106
 107   // phi,lead; this binning starts at -pi/2 and is modulo 3
 108   const Int_t kNLeadingPhiBins = 36;
 109   Double_t leadingPhiBins[kNLeadingPhiBins+1];
 110   for (Int_t i=0; i<=kNLeadingPhiBins; i++)
 111     leadingPhiBins[i] = -TMath::Pi() / 2 + 1.0 / kNLeadingPhiBins * i * TMath::TwoPi();
 112
 113   // multiplicity
 114   const Int_t kNMultiplicityBins = 15;
 115   Double_t multiplicityBins[kNMultiplicityBins+1];
 116   for (Int_t i=0; i<=kNMultiplicityBins; i++)
 117     multiplicityBins[i] = -0.5 + i;
 118   multiplicityBins[kNMultiplicityBins] = 200;
 119
 120   // centrality (in %)
 121   const Int_t kNCentralityBins = 5 + 3 + 8;
 122   Double_t centralityBins[] = { 0, 1, 2, 3, 4, 5, 10, 15, 20, 30, 40, 50, 60, 70, 80, 90, 100.1 };
 123
 124   // particle species
 125   const Int_t kNSpeciesBins = 4; // pi, K, p, rest
 126   Double_t speciesBins[] = { -0.5, 0.5, 1.5, 2.5, 3.5 };
 127
 128   trackBins[3] = multiplicityBins;
 129   iTrackBin[3] = kNMultiplicityBins;
 130   trackAxisTitle[3] = "multiplicity";
 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#phid#eta";
 138   }
 139   else if (strcmp(reqHist, "NumberDensityPhi") == 0)
 140   {
 141     axis = 1;
 142     title = "d^{2}N_{ch}/d#phid#eta";
 143   }
 144   else if (strcmp(reqHist, "NumberDensityPhiCentrality") == 0)
 145   {
 146     axis = 2;
 147     title = "d^{2}N_{ch}/d#phid#eta";
 148   }
 149   else if (strcmp(reqHist, "SumpT") == 0)
 150   {
 151     axis = 0;
 152     title = "d^{2}#Sigma p_{T}/d#phid#eta";
 153   }
 154   else
 155     AliFatal(Form("Invalid histogram requested: %s", reqHist));
 156
 157   Int_t initRegions = fkRegions;
 158
 159   if (axis == 0)
 160   {
 161     trackBins[2] = leadingpTBins;
 162     iTrackBin[2] = kNLeadingpTBins;
 163     trackAxisTitle[2] = "leading p_{T} (GeV/c)";
 164
 165   }
 166   else if (axis == 1)
 167   {
 168     nTrackVars = 5;
 169     initRegions = 1;
 170
 171     iTrackBin[2] = kNLeadingpTBins2;
 172     trackBins[2] = leadingpTBins2;
 173     trackAxisTitle[2] = "leading p_{T} (GeV/c)";
 174
 175     iTrackBin[4] = kNLeadingPhiBins;
 176     trackBins[4] = leadingPhiBins;
 177     trackAxisTitle[4] = "#Delta #phi w.r.t. leading track";
 178   }
 179   else if (axis == 2)
 180   {
 181     nTrackVars = 5;
 182     initRegions = 1;
 183
 184     iTrackBin[0] = kNDeltaEtaBins;
 185     trackBins[0] = deltaEtaBins;
 186     trackAxisTitle[0] = "#Delta#eta";
 187
 188     iTrackBin[2] = kNLeadingpTBins2;
 189     trackBins[2] = leadingpTBins2;
 190     trackAxisTitle[2] = "leading p_{T} (GeV/c)";
 191
 192     trackBins[3] = centralityBins;
 193     iTrackBin[3] = kNCentralityBins;
 194     trackAxisTitle[3] = "centrality";
 195
 196     iTrackBin[4] = kNLeadingPhiBins;
 197     trackBins[4] = leadingPhiBins;
 198     trackAxisTitle[4] = "#Delta#phi (rad.)";
 199   }
 200
 201   for (Int_t i=0; i<initRegions; i++)
 202   {
 203     fTrackHist[i] = new AliCFContainer(Form("fTrackHist_%d", i), title, fgkCFSteps, nTrackVars, iTrackBin);
 204
 205     for (Int_t j=0; j<nTrackVars; j++)
 206     {
 207       fTrackHist[i]->SetBinLimits(j, trackBins[j]);
 208       fTrackHist[i]->SetVarTitle(j, trackAxisTitle[j]);
 209     }
 210
 211     SetStepNames(fTrackHist[i]);
 212   }
 213
 214   // track 3rd and 4th axis --> event 1st and 2nd axis
 215   fEventHist = new AliCFContainer("fEventHist", title, fgkCFSteps, 2, iTrackBin+2);
 216
 217   fEventHist->SetBinLimits(0, trackBins[2]);
 218   fEventHist->SetVarTitle(0, trackAxisTitle[2]);
 219
 220   fEventHist->SetBinLimits(1, trackBins[3]);
 221   fEventHist->SetVarTitle(1, trackAxisTitle[3]);
 222
 223   SetStepNames(fEventHist);
 224
 225   iTrackBin[2] = kNSpeciesBins;
 226
 227   fTrackHistEfficiency = new AliCFContainer("fTrackHistEfficiency", "Tracking efficiency", 3, 4, iTrackBin);
 228   fTrackHistEfficiency->SetBinLimits(0, trackBins[0]);
 229   fTrackHistEfficiency->SetVarTitle(0, trackAxisTitle[0]);
 230   fTrackHistEfficiency->SetBinLimits(1, trackBins[1]);
 231   fTrackHistEfficiency->SetVarTitle(1, trackAxisTitle[1]);
 232   fTrackHistEfficiency->SetBinLimits(2, speciesBins);
 233   fTrackHistEfficiency->SetVarTitle(2, "particle species");
 234   fTrackHistEfficiency->SetBinLimits(3, trackBins[3]);
 235   fTrackHistEfficiency->SetVarTitle(3, trackAxisTitle[3]);
 236 }
 237
 238 //_____________________________________________________________________________
 239 AliUEHist::AliUEHist(const AliUEHist &c) :
 240   TObject(),
 241   fkRegions(4),
 242   fEventHist(0),
 243   fTrackHistEfficiency(0),
 244   fEtaMin(0),
 245   fEtaMax(0),
 246   fPtMin(0),
 247   fPtMax(0),
 248   fCentralityMin(0),
 249   fCentralityMax(0),
 250   fContaminationEnhancement(0),
 251   fCombineMinMax(0),
 252   fCache(0),
 253   fHistogramType()
 254 {
 255   //
 256   // AliUEHist copy constructor
 257   //
 258
 259   ((AliUEHist &) c).Copy(*this);
 260 }
 261
 262 //____________________________________________________________________
 263 void AliUEHist::SetStepNames(AliCFContainer* container)
 264 {
 265   // sets the names of the correction steps
 266
 267   for (Int_t i=0; i<fgkCFSteps; i++)
 268     container->SetStepTitle(i, GetStepTitle((CFStep) i));
 269 }
 270
 271 //____________________________________________________________________
 272 AliUEHist::~AliUEHist()
 273 {
 274   // Destructor
 275
 276   for (Int_t i=0; i<fkRegions; i++)
 277   {
 278     if (fTrackHist[i])
 279     {
 280       delete fTrackHist[i];
 281       fTrackHist[i] = 0;
 282     }
 283   }
 284
 285   if (fEventHist)
 286   {
 287     delete fEventHist;
 288     fEventHist = 0;
 289   }
 290
 291   if (fTrackHistEfficiency)
 292   {
 293     delete fTrackHistEfficiency;
 294     fTrackHistEfficiency = 0;
 295   }
 296
 297   if (fCache)
 298   {
 299     delete fCache;
 300     fCache = 0;
 301   }
 302 }
 303
 304 //____________________________________________________________________
 305 AliUEHist &AliUEHist::operator=(const AliUEHist &c)
 306 {
 307   // assigment operator
 308
 309   if (this != &c)
 310     ((AliUEHist &) c).Copy(*this);
 311
 312   return *this;
 313 }
 314
 315 //____________________________________________________________________
 316 void AliUEHist::Copy(TObject& c) const
 317 {
 318   // copy function
 319
 320   AliUEHist& target = (AliUEHist &) c;
 321
 322   for (Int_t i=0; i<fkRegions; i++)
 323     if (fTrackHist[i])
 324       target.fTrackHist[i] = dynamic_cast<AliCFContainer*> (fTrackHist[i]->Clone());
 325
 326   if (fEventHist)
 327     target.fEventHist = dynamic_cast<AliCFContainer*> (fEventHist->Clone());
 328
 329   if (fTrackHistEfficiency)
 330     target.fTrackHistEfficiency = dynamic_cast<AliCFContainer*> (fTrackHistEfficiency->Clone());
 331
 332   target.fEtaMin = fEtaMin;
 333   target.fEtaMax = fEtaMax;
 334   target.fPtMin = fPtMin;
 335   target.fPtMax = fPtMax;
 336   target.fCentralityMin = fCentralityMin;
 337   target.fCentralityMax = fCentralityMax;
 338
 339   if (fContaminationEnhancement)
 340     target.fContaminationEnhancement = dynamic_cast<TH1F*> (fContaminationEnhancement->Clone());
 341
 342   target.fCombineMinMax = fCombineMinMax;
 343   target.fHistogramType = fHistogramType;
 344 }
 345
 346 //____________________________________________________________________
 347 Long64_t AliUEHist::Merge(TCollection* list)
 348 {
 349   // Merge a list of AliUEHist objects with this (needed for
 350   // PROOF).
 351   // Returns the number of merged objects (including this).
 352
 353   if (!list)
 354     return 0;
 355
 356   if (list->IsEmpty())
 357     return 1;
 358
 359   TIterator* iter = list->MakeIterator();
 360   TObject* obj;
 361
 362   // collections of objects
 363   const Int_t kMaxLists = fkRegions+2;
 364   TList** lists = new TList*[kMaxLists];
 365
 366   for (Int_t i=0; i<kMaxLists; i++)
 367     lists[i] = new TList;
 368
 369   Int_t count = 0;
 370   while ((obj = iter->Next())) {
 371
 372     AliUEHist* entry = dynamic_cast<AliUEHist*> (obj);
 373     if (entry == 0)
 374       continue;
 375
 376     for (Int_t i=0; i<fkRegions; i++)
 377       if (entry->fTrackHist[i])
 378         lists[i]->Add(entry->fTrackHist[i]);
 379
 380     lists[fkRegions]->Add(entry->fEventHist);
 381     lists[fkRegions+1]->Add(entry->fTrackHistEfficiency);
 382
 383     count++;
 384   }
 385   for (Int_t i=0; i<fkRegions; i++)
 386     if (fTrackHist[i])
 387       fTrackHist[i]->Merge(lists[i]);
 388
 389   fEventHist->Merge(lists[fkRegions]);
 390   fTrackHistEfficiency->Merge(lists[fkRegions+1]);
 391
 392   for (Int_t i=0; i<kMaxLists; i++)
 393     delete lists[i];
 394
 395   delete[] lists;
 396
 397   return count+1;
 398 }
 399
 400 //____________________________________________________________________
 401 void AliUEHist::SetBinLimits(AliCFGridSparse* grid)
 402 {
 403   // sets the bin limits in eta and pT defined by fEtaMin/Max, fPtMin/Max
 404
 405   if (fEtaMax > fEtaMin)
 406     grid->SetRangeUser(0, fEtaMin, fEtaMax);
 407   if (fPtMax > fPtMin)
 408     grid->SetRangeUser(1, fPtMin, fPtMax);
 409 }
 410
 411 //____________________________________________________________________
 412 void AliUEHist::ResetBinLimits(AliCFGridSparse* grid)
 413 {
 414   // resets all bin limits
 415
 416   for (Int_t i=0; i<grid->GetNVar(); i++)
 417     if (grid->GetGrid()->GetAxis(i)->TestBit(TAxis::kAxisRange))
 418       grid->SetRangeUser(i, 0, -1);
 419 }
 420
 421 //____________________________________________________________________
 422 void AliUEHist::CountEmptyBins(AliUEHist::CFStep step, Float_t ptLeadMin, Float_t ptLeadMax)
 423 {
 424   // prints the number of empty bins in the track end event histograms in the given step
 425
 426   Int_t binBegin[4];
 427   Int_t binEnd[4];
 428
 429   for (Int_t i=0; i<4; i++)
 430   {
 431     binBegin[i] = 1;
 432     binEnd[i]   = fTrackHist[0]->GetNBins(i);
 433   }
 434
 435   if (fEtaMax > fEtaMin)
 436   {
 437     binBegin[0] = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(0)->FindBin(fEtaMin);
 438     binEnd[0]   = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(0)->FindBin(fEtaMax);
 439   }
 440
 441   if (fPtMax > fPtMin)
 442   {
 443     binBegin[1] = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(1)->FindBin(fPtMin);
 444     binEnd[1]   = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(1)->FindBin(fPtMax);
 445   }
 446
 447   if (ptLeadMax > ptLeadMin)
 448   {
 449     binBegin[2] = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(2)->FindBin(ptLeadMin);
 450     binEnd[2]   = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(2)->FindBin(ptLeadMax);
 451   }
 452
 453   // start from multiplicity 1
 454   binBegin[3] = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(3)->FindBin(1);
 455
 456   for (Int_t region=0; region<fkRegions; region++)
 457   {
 458     Int_t total = 0;
 459     Int_t count = 0;
 460     Int_t vars[4];
 461
 462     for (Int_t i=0; i<4; i++)
 463       vars[i] = binBegin[i];
 464
 465     AliCFGridSparse* grid = fTrackHist[region]->GetGrid(step);
 466     while (1)
 467     {
 468       if (grid->GetElement(vars) == 0)
 469       {
 470         Printf("Empty bin at eta=%.2f pt=%.2f pt_lead=%.2f mult=%.1f",
 471           grid->GetBinCenter(0, vars[0]),
 472           grid->GetBinCenter(1, vars[1]),
 473           grid->GetBinCenter(2, vars[2]),
 474           grid->GetBinCenter(3, vars[3])
 475         );
 476         count++;
 477       }
 478
 479       vars[3]++;
 480       for (Int_t i=3; i>0; i--)
 481       {
 482         if (vars[i] == binEnd[i]+1)
 483         {
 484           vars[i] = binBegin[i];
 485           vars[i-1]++;
 486         }
 487       }
 488
 489       if (vars[0] == binEnd[0]+1)
 490         break;
 491       total++;
 492     }
 493
 494     Printf("Region %s has %d empty bins (out of %d bins)", GetRegionTitle((Region) region), count, total);
 495   }
 496 }
 497
 498 //____________________________________________________________________
 499 TH1* AliUEHist::GetUEHist(AliUEHist::CFStep step, AliUEHist::Region region, Float_t ptLeadMin, Float_t ptLeadMax, Int_t multBinBegin, Int_t multBinEnd, Int_t twoD)
 500 {
 501   // Extracts the UE histogram at the given step and in the given region by projection and dividing tracks by events
 502   //
 503   // ptLeadMin, ptLeadMax: Only meaningful for vs. delta phi plot (third axis is ptLead)
 504   // Histogram has to be deleted by the caller of the function
 505   //
 506   // twoD: 0: 1D histogram as function of phi
 507   //       1: 2D histogram, deltaphi, deltaeta
 508   //       10: 1D histogram, within |deltaeta| < 0.8
 509   //       11: 1D histogram, within 0.8 < |deltaeta| < 1.6
 510
 511   // unzoom all axes
 512   ResetBinLimits(fTrackHist[region]->GetGrid(step));
 513   ResetBinLimits(fEventHist->GetGrid(step));
 514
 515   SetBinLimits(fTrackHist[region]->GetGrid(step));
 516
 517   TH1D* tracks = 0;
 518
 519   if (ptLeadMin < 0)
 520   {
 521     tracks = fTrackHist[region]->ShowProjection(2, step);
 522     tracks->GetYaxis()->SetTitle(fTrackHist[region]->GetTitle());
 523     if (fCombineMinMax && region == kMin)
 524     {
 525       ResetBinLimits(fTrackHist[kMax]->GetGrid(step));
 526       SetBinLimits(fTrackHist[kMax]->GetGrid(step));
 527
 528       TH1D* tracks2 = fTrackHist[kMax]->ShowProjection(2, step);
 529       tracks->Add(tracks2);
 530
 531       ResetBinLimits(fTrackHist[kMax]->GetGrid(step));
 532     }
 533
 534     // normalize to get a density (deta dphi)
 535     TAxis* axis = fTrackHist[region]->GetGrid(step)->GetAxis(0);
 536     Float_t phiRegion = TMath::TwoPi() / 3;
 537     if (!fCombineMinMax && region == kMin)
 538       phiRegion /= 2;
 539     Float_t normalization = phiRegion * (axis->GetBinUpEdge(axis->GetLast()) - axis->GetBinLowEdge(axis->GetFirst()));
 540     //Printf("Normalization: %f", normalization);
 541     tracks->Scale(1.0 / normalization);
 542
 543     TH1D* events = fEventHist->ShowProjection(0, step);
 544     tracks->Divide(events);
 545   }
 546   else
 547   {
 548     if (multBinEnd >= multBinBegin)
 549       Printf("Using multiplicity range %d --> %d", multBinBegin, multBinEnd);
 550     fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(3)->SetRange(multBinBegin, multBinEnd);
 551     fEventHist->GetGrid(step)->GetGrid()->GetAxis(1)->SetRange(multBinBegin, multBinEnd);
 552
 553     Int_t firstBin = fTrackHist[region]->GetAxis(2, step)->FindBin(ptLeadMin);
 554     Int_t lastBin = fTrackHist[region]->GetAxis(2, step)->FindBin(ptLeadMax);
 555
 556     Printf("Using leading pT range %d --> %d", firstBin, lastBin);
 557
 558     fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(2)->SetRange(firstBin, lastBin);
 559
 560     if (twoD == 0)
 561       tracks = (TH1D*) fTrackHist[region]->GetGrid(step)->Project(4);
 562     else
 563       tracks = (TH1D*) fTrackHist[region]->GetGrid(step)->Project(4, 0);
 564
 565     Printf("Calculated histogram --> %f tracks", tracks->Integral());
 566     fTrackHist[region]->GetGrid(step)->SetRangeUser(2, 0, -1);
 567
 568     if (twoD == 10 || twoD == 11)
 569     {
 570       Float_t etaLimit = 0.8;
 571       if (twoD == 10)
 572       {
 573         tracks = (TH1D*) ((TH2*) tracks)->ProjectionX("proj", tracks->GetYaxis()->FindBin(-etaLimit + 0.01), tracks->GetYaxis()->FindBin(etaLimit - 0.01))->Clone();
 574
 575         // TODO calculate acc with 2 * (deta - 0.5 * deta*deta / 1.6)
 576         tracks->Scale(1. / 0.75);
 577       }
 578       else if (twoD == 11)
 579       {
 580         TH1D* tracksTmp1 = (TH1D*) ((TH2*) tracks)->ProjectionX("proj1", tracks->GetYaxis()->FindBin(-etaLimit * 2 + 0.01), tracks->GetYaxis()->FindBin(-etaLimit - 0.01))->Clone();
 581         TH1D* tracksTmp2 = ((TH2*) tracks)->ProjectionX("proj2", tracks->GetYaxis()->FindBin(etaLimit + 0.01), tracks->GetYaxis()->FindBin(2 * etaLimit - 0.01));
 582
 583         tracksTmp1->Add(tracksTmp2);
 584
 585         delete tracks;
 586         tracks = tracksTmp1;
 587         delete tracksTmp2;
 588
 589         tracks->Scale(1. / 0.25);
 590       }
 591     }
 592
 593     // normalize to get a density (deta dphi)
 594     // TODO normalization may be off for 2d histogram
 595     Float_t normalization = fTrackHist[region]->GetGrid(step)->GetAxis(4)->GetBinWidth(1);
 596     TAxis* axis = fTrackHist[region]->GetGrid(step)->GetAxis(0);
 597     if (strcmp(axis->GetTitle(), "#eta") == 0)
 598     {
 599       Printf("Normalizing using eta-axis range");
 600       normalization *= axis->GetBinUpEdge(axis->GetLast()) - axis->GetBinLowEdge(axis->GetFirst());
 601     }
 602     else
 603     {
 604       Printf("Normalizing assuming accepted range of |eta| < 0.8");
 605       normalization *= 0.8 * 2;
 606     }
 607
 608     //Printf("Normalization: %f", normalization);
 609     tracks->Scale(1.0 / normalization);
 610
 611     // 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!
 612     TH1D* events = fEventHist->ShowProjection(0, step);
 613     Int_t nEvents = (Int_t) events->Integral(firstBin, lastBin);
 614     Printf("Calculated histogram --> %d events", nEvents);
 615
 616     if (nEvents > 0)
 617       tracks->Scale(1.0 / nEvents);
 618
 619     delete events;
 620   }
 621
 622   ResetBinLimits(fTrackHist[region]->GetGrid(step));
 623   ResetBinLimits(fEventHist->GetGrid(step));
 624
 625   return tracks;
 626 }
 627
 628 //____________________________________________________________________
 629 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)
 630 {
 631   // returns a histogram projected to pT,assoc with the provided cuts
 632
 633    // unzoom all axes
 634   ResetBinLimits(fTrackHist[region]->GetGrid(step));
 635   ResetBinLimits(fEventHist->GetGrid(step));
 636
 637   TH1D* tracks = 0;
 638
 639     // 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
 640   // therefore the number density must be calculated here per leading pT bin and then summed
 641
 642   if (multBinEnd > multBinBegin)
 643     Printf("Using multiplicity range %d --> %d", multBinBegin, multBinEnd);
 644   fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(3)->SetRange(multBinBegin, multBinEnd);
 645   fEventHist->GetGrid(step)->GetGrid()->GetAxis(1)->SetRange(multBinBegin, multBinEnd);
 646
 647   fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(0)->SetRangeUser(etaMin, etaMax);
 648   fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(4)->SetRangeUser(phiMin, phiMax);
 649
 650   // get real phi cuts due to binning
 651   Float_t phiMinReal = fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(4)->GetBinLowEdge(fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(4)->GetFirst());
 652   Float_t phiMaxReal = fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(4)->GetBinUpEdge(fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(4)->GetLast());
 653   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);
 654
 655   Int_t firstBin = fTrackHist[region]->GetAxis(2, step)->FindBin(ptLeadMin);
 656   Int_t lastBin = fTrackHist[region]->GetAxis(2, step)->FindBin(ptLeadMax);
 657
 658   TH1D* events = fEventHist->ShowProjection(0, step);
 659
 660   for (Int_t bin=firstBin; bin<=lastBin; bin++)
 661   {
 662     fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(2)->SetRange(bin, bin);
 663
 664     // project to pT,assoc
 665     TH1D* tracksTmp = (TH1D*) fTrackHist[region]->GetGrid(step)->Project(1);
 666
 667     Printf("Calculated histogram in bin %d --> %f tracks", bin, tracksTmp->Integral());
 668     fTrackHist[region]->GetGrid(step)->SetRangeUser(2, 0, -1);
 669
 670     // normalize to get a density (deta dphi)
 671     Float_t normalization = 1;
 672     TAxis* axis = fTrackHist[region]->GetGrid(step)->GetAxis(0);
 673     if (strcmp(axis->GetTitle(), "#eta") == 0)
 674     {
 675       Printf("Normalizing using eta-axis range");
 676       normalization *= axis->GetBinUpEdge(axis->GetLast()) - axis->GetBinLowEdge(axis->GetFirst());
 677     }
 678     else
 679     {
 680       Printf("Normalizating assuming accepted range of |eta| < 0.8");
 681       normalization *= 0.8 * 2;
 682     }
 683
 684     // dphi
 685     if (!skipPhiNormalization)
 686       normalization *= phiMaxReal - phiMinReal;
 687
 688     //Printf("Normalization: %f", normalization);
 689     tracksTmp->Scale(1.0 / normalization);
 690
 691     // and now dpT (bins have different width)
 692     for (Int_t i=1; i<=tracksTmp->GetNbinsX(); i++)
 693     {
 694       tracksTmp->SetBinContent(i, tracksTmp->GetBinContent(i) / tracksTmp->GetXaxis()->GetBinWidth(i));
 695       tracksTmp->SetBinError  (i, tracksTmp->GetBinError(i) / tracksTmp->GetXaxis()->GetBinWidth(i));
 696     }
 697
 698     Int_t nEvents = (Int_t) events->GetBinContent(bin);
 699     if (nEvents > 0)
 700       tracksTmp->Scale(1.0 / nEvents);
 701     Printf("Calculated histogram in bin %d --> %d events", bin, nEvents);
 702
 703     if (!tracks)
 704       tracks = tracksTmp;
 705     else
 706     {
 707       tracks->Add(tracksTmp);
 708       delete tracksTmp;
 709     }
 710   }
 711
 712   delete events;
 713
 714   ResetBinLimits(fTrackHist[region]->GetGrid(step));
 715   ResetBinLimits(fEventHist->GetGrid(step));
 716
 717   return tracks;
 718 }
 719
 720 //____________________________________________________________________
 721 void AliUEHist::CorrectTracks(CFStep step1, CFStep step2, TH1* trackCorrection, Int_t var1, Int_t var2)
 722 {
 723   // corrects from step1 to step2 by multiplying the tracks with trackCorrection
 724   // trackCorrection can be a function of eta (var1 == 0), pT (var1 == 1), leading pT (var1 == 2), multiplicity (var1 == 3), delta phi (var1 == 4)
 725   // if var2 >= 0 a two dimension correction is assumed in trackCorrection
 726   //
 727   // if trackCorrection is 0, just copies content from step1 to step2
 728
 729   for (Int_t region=0; region<fkRegions; region++)
 730     CorrectTracks(step1, step2, region, trackCorrection, var1, var2);
 731 }
 732
 733 //____________________________________________________________________
 734 void AliUEHist::CorrectTracks(CFStep step1, CFStep step2, Int_t region, TH1* trackCorrection, Int_t var1, Int_t var2)
 735 {
 736   //
 737   // see documentation of CorrectTracks above
 738   //
 739
 740   if (!fTrackHist[region])
 741     return;
 742
 743   THnSparse* grid = fTrackHist[region]->GetGrid(step1)->GetGrid();
 744   THnSparse* target = fTrackHist[region]->GetGrid(step2)->GetGrid();
 745
 746   // clear target histogram
 747   target->Reset();
 748
 749   if (trackCorrection != 0)
 750   {
 751     if (grid->GetAxis(var1)->GetNbins() != trackCorrection->GetNbinsX())
 752       AliFatal(Form("Invalid binning (var1): %d %d", grid->GetAxis(var1)->GetNbins(), trackCorrection->GetNbinsX()));
 753
 754     if (var2 >= 0 && grid->GetAxis(var2)->GetNbins() != trackCorrection->GetNbinsY())
 755       AliFatal(Form("Invalid binning (var2): %d %d", grid->GetAxis(var2)->GetNbins(), trackCorrection->GetNbinsY()));
 756   }
 757
 758   // optimized implementation
 759   for (Int_t binIdx = 0; binIdx < grid->GetNbins(); binIdx++)
 760   {
 761     Int_t bins[5];
 762     Double_t value = grid->GetBinContent(binIdx, bins);
 763     Double_t error = grid->GetBinError(binIdx);
 764
 765     if (trackCorrection != 0)
 766     {
 767       if (var2 < 0)
 768       {
 769         value *= trackCorrection->GetBinContent(bins[var1]);
 770         error *= trackCorrection->GetBinContent(bins[var1]);
 771       }
 772       else
 773       {
 774         value *= trackCorrection->GetBinContent(bins[var1], bins[var2]);
 775         error *= trackCorrection->GetBinContent(bins[var1], bins[var2]);
 776       }
 777     }
 778
 779     target->SetBinContent(bins, value);
 780     target->SetBinError(bins, error);
 781   }
 782
 783   Printf("AliUEHist::CorrectTracks: Corrected from %f to %f entries. Correction histogram: %f entries (integral: %f)", grid->GetEntries(), target->GetEntries(), (trackCorrection) ? trackCorrection->GetEntries() : -1.0, (trackCorrection) ? trackCorrection->Integral() : -1.0);
 784 }
 785
 786 //____________________________________________________________________
 787 void AliUEHist::CorrectEvents(CFStep step1, CFStep step2, TH1* eventCorrection, Int_t var1, Int_t var2)
 788 {
 789   // corrects from step1 to step2 by multiplying the events with eventCorrection
 790   // eventCorrection is as function of leading pT (var == 0) or multiplicity (var == 1)
 791   //
 792   // if eventCorrection is 0, just copies content from step1 to step2
 793
 794   AliCFGridSparse* grid = fEventHist->GetGrid(step1);
 795   AliCFGridSparse* target = fEventHist->GetGrid(step2);
 796
 797   // clear target histogram
 798   target->GetGrid()->Reset();
 799
 800   if (eventCorrection != 0 && grid->GetNBins(var1) != eventCorrection->GetNbinsX())
 801     AliFatal(Form("Invalid binning: %d %d", grid->GetNBins(var1), eventCorrection->GetNbinsX()));
 802
 803   if (eventCorrection != 0 && var2 != -1 && grid->GetNBins(var2) != eventCorrection->GetNbinsY())
 804     AliFatal(Form("Invalid binning: %d %d", grid->GetNBins(var2), eventCorrection->GetNbinsY()));
 805
 806   Int_t bins[2];
 807   for (Int_t x = 1; x <= grid->GetNBins(0); x++)
 808   {
 809     bins[0] = x;
 810     for (Int_t y = 1; y <= grid->GetNBins(1); y++)
 811     {
 812       bins[1] = y;
 813
 814       Double_t value = grid->GetElement(bins);
 815       if (value != 0)
 816       {
 817         Double_t error = grid->GetElementError(bins);
 818
 819         if (eventCorrection != 0)
 820         {
 821           if (var2 == -1)
 822           {
 823             value *= eventCorrection->GetBinContent(bins[var1]);
 824             error *= eventCorrection->GetBinContent(bins[var1]);
 825           }
 826           else
 827           {
 828             value *= eventCorrection->GetBinContent(bins[var1], bins[var2]);
 829             error *= eventCorrection->GetBinContent(bins[var1], bins[var2]);
 830           }
 831         }
 832
 833         target->SetElement(bins, value);
 834         target->SetElementError(bins, error);
 835       }
 836     }
 837   }
 838
 839   Printf("AliUEHist::CorrectEvents: Corrected from %f to %f entries. Correction histogram: %f entries (integral: %f)", grid->GetEntries(), target->GetEntries(), (eventCorrection) ? eventCorrection->GetEntries() : -1.0, (eventCorrection) ? eventCorrection->Integral() : -1.0);
 840 }
 841
 842 //____________________________________________________________________
 843 void AliUEHist::Correct(AliUEHist* corrections)
 844 {
 845   // applies the given corrections to extract from the step kCFStepReconstructed all previous steps
 846   //
 847   // in this object the data is expected in the step kCFStepReconstructed
 848
 849   if (fHistogramType.Length() == 0)
 850   {
 851     Printf("W-AliUEHist::Correct: fHistogramType not defined. Guessing histogram type...");
 852     if (fTrackHist[kToward]->GetNVar() < 5)
 853     {
 854       if (strcmp(fTrackHist[kToward]->GetTitle(), "d^{2}N_{ch}/d#phid#eta") == 0)
 855         fHistogramType = "NumberDensitypT";
 856       else if (strcmp(fTrackHist[kToward]->GetTitle(), "d^{2}#Sigma p_{T}/d#phid#eta") == 0)
 857         fHistogramType = "SumpT";
 858     }
 859     else if (fTrackHist[kToward]->GetNVar() == 5)
 860     {
 861       if (strcmp(fTrackHist[kToward]->GetGrid(0)->GetGrid()->GetAxis(3)->GetTitle(), "multiplicity") == 0)
 862         fHistogramType = "NumberDensityPhi";
 863       else if (strcmp(fTrackHist[kToward]->GetGrid(0)->GetGrid()->GetAxis(3)->GetTitle(), "centrality") == 0)
 864         fHistogramType = "NumberDensityPhiCentrality";
 865     }
 866
 867     if (fHistogramType.Length() == 0)
 868       AliFatal("Cannot figure out histogram type. Try setting it manually...");
 869   }
 870
 871   Printf("AliUEHist::Correct: Correcting %s...", fHistogramType.Data());
 872
 873   if (strcmp(fHistogramType, "NumberDensitypT") == 0 || strcmp(fHistogramType, "NumberDensityPhi") == 0 || strcmp(fHistogramType, "SumpT") == 0) // the last is for backward compatibilty
 874   {
 875     // ---- track level
 876
 877     // bias due to migration in leading pT (because the leading particle is not reconstructed, and the subleading is used)
 878     // extracted as function of leading pT
 879     Bool_t biasFromMC = kFALSE;
 880     const char* projAxis = "z";
 881     Int_t secondBin = -1;
 882
 883     if (strcmp(fHistogramType, "NumberDensityPhi") == 0)
 884     {
 885       projAxis = "yz";
 886       secondBin = 4;
 887     }
 888
 889     for (Int_t region = 0; region < fkRegions; region++)
 890     {
 891       if (!fTrackHist[region])
 892         continue;
 893
 894       TH1* leadingBiasTracks =  0;
 895       if (biasFromMC)
 896       {
 897         leadingBiasTracks = (TH1*) corrections->GetBias(kCFStepReconstructed, kCFStepTracked, region, projAxis, 0, -1, 1); // from MC
 898         Printf("WARNING: Using MC bias correction");
 899       }
 900       else
 901         leadingBiasTracks = (TH1*) GetBias(kCFStepBiasStudy, kCFStepReconstructed, region, projAxis, 0, -1, 1);          // from data
 902
 903       CorrectTracks(kCFStepReconstructed, kCFStepTracked, region, leadingBiasTracks, 2, secondBin);
 904       if (region == kMin && fCombineMinMax)
 905       {
 906         CorrectTracks(kCFStepReconstructed, kCFStepTracked, kMax, leadingBiasTracks, 2, secondBin);
 907         delete leadingBiasTracks;
 908         break;
 909       }
 910       delete leadingBiasTracks;
 911     }
 912
 913     TH1* leadingBiasEvents = 0;
 914     if (biasFromMC)
 915       leadingBiasEvents = (TH1*) corrections->GetBias(kCFStepReconstructed, kCFStepTracked, kToward, projAxis, 0, -1, 2); // from MC
 916     else
 917       leadingBiasEvents = (TH1*) GetBias(kCFStepBiasStudy, kCFStepReconstructed, kToward, projAxis, 0, -1, 2);          // from data
 918
 919     //new TCanvas; leadingBiasEvents->DrawCopy();
 920     CorrectEvents(kCFStepReconstructed, kCFStepTracked, leadingBiasEvents, 0);
 921     delete leadingBiasEvents;
 922
 923     // --- contamination correction ---
 924
 925     TH2D* contamination = corrections->GetTrackingContamination();
 926     if (corrections->fContaminationEnhancement)
 927     {
 928       Printf("Applying contamination enhancement");
 929
 930       for (Int_t x=1; x<=contamination->GetNbinsX(); x++)
 931         for (Int_t y=1; y<=contamination->GetNbinsX(); y++)
 932           contamination->SetBinContent(x, y, contamination->GetBinContent(x, y) * corrections->fContaminationEnhancement->GetBinContent(corrections->fContaminationEnhancement->GetXaxis()->FindBin(contamination->GetYaxis()->GetBinCenter(y))));
 933     }
 934     CorrectTracks(kCFStepTracked, kCFStepTrackedOnlyPrim, contamination, 0, 1);
 935     CorrectEvents(kCFStepTracked, kCFStepTrackedOnlyPrim, 0, 0);
 936     delete contamination;
 937
 938     // --- efficiency correction ---
 939     // correct single-particle efficiency for associated particles
 940     // in addition correct for efficiency on trigger particles (tracks AND events)
 941
 942     TH1* efficiencyCorrection = corrections->GetTrackingEfficiencyCorrection();
 943     // use kCFStepVertex as a temporary step (will be overwritten later anyway)
 944     CorrectTracks(kCFStepTrackedOnlyPrim, kCFStepVertex, efficiencyCorrection, 0, 1);
 945     delete efficiencyCorrection;
 946
 947     // correct pT,T
 948     efficiencyCorrection = corrections->GetTrackEfficiency(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, 1, -1, 2);
 949     CorrectEvents(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, efficiencyCorrection, 0);
 950     CorrectTracks(kCFStepVertex, kCFStepAnaTopology, efficiencyCorrection, 2);
 951     delete efficiencyCorrection;
 952
 953     // copy
 954     CorrectTracks(kCFStepAnaTopology, kCFStepVertex, 0, -1);
 955     CorrectEvents(kCFStepAnaTopology, kCFStepVertex, 0, 0);
 956
 957     // vertex correction on the level of events as function of multiplicity, weighting tracks and events with the same factor
 958     // practically independent of low pT cut
 959     TH1D* vertexCorrection = (TH1D*) corrections->GetEventEfficiency(kCFStepVertex, kCFStepTriggered, 1);
 960
 961     // convert stage from true multiplicity to observed multiplicity by simple conversion factor
 962     TH1D* vertexCorrectionObs = (TH1D*) vertexCorrection->Clone("vertexCorrection2");
 963     vertexCorrectionObs->Reset();
 964
 965     TF1* func = new TF1("func", "[1]+[0]/(x-[2])");
 966     // some defaults
 967     func->SetParameters(0.1, 1, -0.7);
 968     vertexCorrection->Fit(func, "0I", "", 0, 3);
 969     for (Int_t i=1; i<=vertexCorrectionObs->GetNbinsX(); i++)
 970     {
 971       Float_t xPos = 1.0 / 0.77 * vertexCorrectionObs->GetXaxis()->GetBinCenter(i);
 972       if (xPos < 3)
 973         vertexCorrectionObs->SetBinContent(i, func->Eval(xPos));
 974       else
 975         vertexCorrectionObs->SetBinContent(i, vertexCorrection->Interpolate(xPos));
 976     }
 977
 978     #if 0
 979     new TCanvas;
 980     vertexCorrection->DrawCopy();
 981     vertexCorrectionObs->SetLineColor(2);
 982     vertexCorrectionObs->DrawCopy("same");
 983     func->SetRange(0, 4);
 984     func->DrawClone("same");
 985     #endif
 986
 987     CorrectTracks(kCFStepVertex, kCFStepTriggered, vertexCorrectionObs, 3);
 988     CorrectEvents(kCFStepVertex, kCFStepTriggered, vertexCorrectionObs, 1);
 989     delete vertexCorrectionObs;
 990     delete vertexCorrection;
 991     delete func;
 992
 993     // copy
 994     CorrectTracks(kCFStepTriggered, kCFStepAll, 0, -1);
 995     CorrectEvents(kCFStepTriggered, kCFStepAll, 0, 0);
 996   }
 997   else if (strcmp(fHistogramType, "NumberDensityPhiCentrality") == 0)
 998   {
 999     // copy
1000     CorrectTracks(kCFStepReconstructed, kCFStepTracked, 0, -1);
1001     CorrectEvents(kCFStepReconstructed, kCFStepTracked, 0, -1);
1002
1003     // Dont use eta in the following, because it is a Delta-eta axis
1004
1005     // contamination correction
1006     // correct single-particle contamination for associated particles
1007
1008     TH1* contamination = corrections->GetTrackingContamination(1);
1009
1010     if (1)
1011     {
1012       Printf("Applying contamination enhancement");
1013
1014       for (Int_t bin = 1; bin <= contamination->GetNbinsX(); bin++)
1015       {
1016         printf("%f", contamination->GetBinContent(bin));
1017         if (contamination->GetBinContent(bin) > 0)
1018           contamination->SetBinContent(bin, 1.0 + 1.1 * (contamination->GetBinContent(bin) - 1.0));
1019         printf(" --> %f\n", contamination->GetBinContent(bin));
1020       }
1021     }
1022
1023     CorrectTracks(kCFStepTracked, kCFStepTrackedOnlyPrim, contamination, 1);
1024     delete contamination;
1025
1026     // correct for additional contamination due to trigger particle around phi ~ 0
1027     TH2* correlatedContamination = corrections->GetCorrelatedContamination();
1028     if (1)
1029     {
1030       Printf("Applying contamination enhancement");
1031
1032       for (Int_t bin = 1; bin <= correlatedContamination->GetNbinsX(); bin++)
1033         for (Int_t bin2 = 1; bin2 <= correlatedContamination->GetNbinsY(); bin2++)
1034         {
1035           printf("%f", correlatedContamination->GetBinContent(bin, bin2));
1036           if (correlatedContamination->GetBinContent(bin, bin2) > 0)
1037             correlatedContamination->SetBinContent(bin, bin2, 1.0 + 1.1 * (correlatedContamination->GetBinContent(bin, bin2) - 1.0));
1038           printf(" --> %f\n", correlatedContamination->GetBinContent(bin, bin2));
1039         }
1040     }
1041
1042     //new TCanvas; correlatedContamination->DrawCopy("COLZ");
1043     CorrectCorrelatedContamination(kCFStepTrackedOnlyPrim, 0, correlatedContamination);
1044
1045     delete correlatedContamination;
1046
1047     // TODO correct for contamination of trigger particles (for tracks AND events)
1048     CorrectEvents(kCFStepTracked, kCFStepTrackedOnlyPrim, 0, 0);
1049
1050     // --- efficiency correction ---
1051     // correct single-particle efficiency for associated particles
1052     // in addition correct for efficiency on trigger particles (tracks AND events)
1053
1054     // in bins of pT and centrality
1055     TH1* efficiencyCorrection = corrections->GetTrackingEfficiencyCorrectionCentrality();
1056     // use kCFStepAnaTopology as a temporary step
1057     CorrectTracks(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, efficiencyCorrection, 1, 3);
1058     delete efficiencyCorrection;
1059
1060     // correct pT,T in bins of pT and centrality
1061     efficiencyCorrection = corrections->GetTrackEfficiency(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, 1, 3, 2);
1062     CorrectEvents(kCFStepTrackedOnlyPrim, kCFStepVertex, efficiencyCorrection, 0, 1);
1063     CorrectTracks(kCFStepAnaTopology, kCFStepVertex, efficiencyCorrection, 2, 3);
1064     delete efficiencyCorrection;
1065
1066     // no correction for vertex finding efficiency and trigger efficiency needed in PbPb
1067     // copy
1068     CorrectTracks(kCFStepVertex, kCFStepAll, 0, -1);
1069     CorrectEvents(kCFStepVertex, kCFStepAll, 0, -1);
1070   }
1071   else
1072     AliFatal(Form("Unknown histogram for correction: %s", GetTitle()));
1073 }
1074
1075 //____________________________________________________________________
1076 TH1* AliUEHist::GetTrackEfficiency(CFStep step1, CFStep step2, Int_t axis1, Int_t axis2, Int_t source, Int_t axis3)
1077 {
1078   // creates a track-level efficiency by dividing step2 by step1
1079   // projected to axis1 and axis2 (optional if >= 0)
1080   //
1081   // source: 0 = fTrackHist; 1 = fTrackHistEfficiency; 2 = fTrackHistEfficiency rebinned for pT,T / pT,lead binning
1082
1083   // integrate over regions
1084   // cache it for efficiency (usually more than one efficiency is requested)
1085
1086   AliCFContainer* sourceContainer = 0;
1087
1088   if (source == 0)
1089   {
1090     if (!fCache)
1091     {
1092       fCache = (AliCFContainer*) fTrackHist[0]->Clone();
1093       for (Int_t i = 1; i < fkRegions; i++)
1094         if (fTrackHist[i])
1095           fCache->Add(fTrackHist[i]);
1096     }
1097     sourceContainer = fCache;
1098   }
1099   else if (source == 1 || source == 2)
1100   {
1101     sourceContainer = fTrackHistEfficiency;
1102     // step offset because we start with kCFStepAnaTopology
1103     step1 = (CFStep) ((Int_t) step1 - (Int_t) kCFStepAnaTopology);
1104     step2 = (CFStep) ((Int_t) step2 - (Int_t) kCFStepAnaTopology);
1105   }
1106   else
1107     return 0;
1108
1109   // reset all limits and set the right ones except those in axis1, axis2 and axis3
1110   ResetBinLimits(sourceContainer->GetGrid(step1));
1111   ResetBinLimits(sourceContainer->GetGrid(step2));
1112   if (fEtaMax > fEtaMin && axis1 != 0 && axis2 != 0 && axis3 != 0)
1113   {
1114     sourceContainer->GetGrid(step1)->SetRangeUser(0, fEtaMin, fEtaMax);
1115     sourceContainer->GetGrid(step2)->SetRangeUser(0, fEtaMin, fEtaMax);
1116   }
1117   if (fPtMax > fPtMin && axis1 != 1 && axis2 != 1 && axis3 != 1)
1118   {
1119     sourceContainer->GetGrid(step1)->SetRangeUser(1, fPtMin, fPtMax);
1120     sourceContainer->GetGrid(step2)->SetRangeUser(1, fPtMin, fPtMax);
1121   }
1122   if (fCentralityMax > fCentralityMin && axis1 != 3 && axis2 != 3 && axis3 != 3)
1123   {
1124     sourceContainer->GetGrid(step1)->SetRangeUser(3, fCentralityMin, fCentralityMax);
1125     sourceContainer->GetGrid(step2)->SetRangeUser(3, fCentralityMin, fCentralityMax);
1126   }
1127
1128   TH1* measured = 0;
1129   TH1* generated = 0;
1130
1131   if (axis3 >= 0)
1132   {
1133     generated = sourceContainer->Project(step1, axis1, axis2, axis3);
1134     measured = sourceContainer->Project(step2, axis1, axis2, axis3);
1135   }
1136   else if (axis2 >= 0)
1137   {
1138     generated = sourceContainer->Project(step1, axis1, axis2);
1139     measured = sourceContainer->Project(step2, axis1, axis2);
1140   }
1141   else
1142   {
1143     generated = sourceContainer->Project(step1, axis1);
1144     measured = sourceContainer->Project(step2, axis1);
1145   }
1146
1147   // check for bins with less than 50 entries, print warning
1148   Int_t binBegin[3];
1149   Int_t binEnd[3];
1150
1151   binBegin[0] = 1;
1152   binBegin[1] = 1;
1153   binBegin[2] = 1;
1154
1155   binEnd[0] = generated->GetNbinsX();
1156   binEnd[1] = generated->GetNbinsY();
1157   binEnd[2] = generated->GetNbinsZ();
1158
1159   if (fEtaMax > fEtaMin)
1160   {
1161     if (axis1 == 0)
1162     {
1163       binBegin[0] = generated->GetXaxis()->FindBin(fEtaMin);
1164       binEnd[0]   = generated->GetXaxis()->FindBin(fEtaMax);
1165     }
1166     if (axis2 == 0)
1167     {
1168       binBegin[1] = generated->GetYaxis()->FindBin(fEtaMin);
1169       binEnd[1]   = generated->GetYaxis()->FindBin(fEtaMax);
1170     }
1171     if (axis3 == 0)
1172     {
1173       binBegin[2] = generated->GetZaxis()->FindBin(fEtaMin);
1174       binEnd[2]   = generated->GetZaxis()->FindBin(fEtaMax);
1175     }
1176   }
1177
1178   if (fPtMax > fPtMin)
1179   {
1180     // TODO this is just checking up to 15 for now
1181     Float_t ptMax = TMath::Min((Float_t) 15., fPtMax);
1182     if (axis1 == 1)
1183     {
1184       binBegin[0] = generated->GetXaxis()->FindBin(fPtMin);
1185       binEnd[0]   = generated->GetXaxis()->FindBin(ptMax);
1186     }
1187     if (axis2 == 1)
1188     {
1189       binBegin[1] = generated->GetYaxis()->FindBin(fPtMin);
1190       binEnd[1]   = generated->GetYaxis()->FindBin(ptMax);
1191     }
1192     if (axis3 == 1)
1193     {
1194       binBegin[2] = generated->GetZaxis()->FindBin(fPtMin);
1195       binEnd[2]   = generated->GetZaxis()->FindBin(ptMax);
1196     }
1197   }
1198
1199   Int_t total = 0;
1200   Int_t count = 0;
1201   Int_t vars[3];
1202
1203   for (Int_t i=0; i<3; i++)
1204     vars[i] = binBegin[i];
1205
1206   const Int_t limit = 50;
1207   while (1)
1208   {
1209     if (generated->GetDimension() == 1 && generated->GetBinContent(vars[0]) < limit)
1210     {
1211       Printf("Empty bin at %s=%.2f (%.2f entries)", generated->GetXaxis()->GetTitle(), generated->GetXaxis()->GetBinCenter(vars[0]), generated->GetBinContent(vars[0]));
1212       count++;
1213     }
1214     else if (generated->GetDimension() == 2 && generated->GetBinContent(vars[0], vars[1]) < limit)
1215     {
1216       Printf("Empty bin at %s=%.2f %s=%.2f (%.2f entries)",
1217         generated->GetXaxis()->GetTitle(), generated->GetXaxis()->GetBinCenter(vars[0]),
1218         generated->GetYaxis()->GetTitle(), generated->GetYaxis()->GetBinCenter(vars[1]),
1219         generated->GetBinContent(vars[0], vars[1]));
1220       count++;
1221     }
1222     else if (generated->GetDimension() == 3 && generated->GetBinContent(vars[0], vars[1], vars[2]) < limit)
1223     {
1224       Printf("Empty bin at %s=%.2f %s=%.2f %s=%.2f (%.2f entries)",
1225         generated->GetXaxis()->GetTitle(), generated->GetXaxis()->GetBinCenter(vars[0]),
1226         generated->GetYaxis()->GetTitle(), generated->GetYaxis()->GetBinCenter(vars[1]),
1227         generated->GetZaxis()->GetTitle(), generated->GetZaxis()->GetBinCenter(vars[2]),
1228         generated->GetBinContent(vars[0], vars[1], vars[2]));
1229       count++;
1230     }
1231
1232     vars[2]++;
1233     if (vars[2] == binEnd[2]+1)
1234     {
1235       vars[2] = binBegin[2];
1236       vars[1]++;
1237     }
1238
1239     if (vars[1] == binEnd[1]+1)
1240     {
1241       vars[1] = binBegin[1];
1242       vars[0]++;
1243     }
1244
1245     if (vars[0] == binEnd[0]+1)
1246       break;
1247     total++;
1248   }
1249
1250   Printf("Correction has %d empty bins (out of %d bins)", count, total);
1251
1252   // rebin if required
1253   if (source == 2)
1254   {
1255     TAxis* axis = fEventHist->GetGrid(0)->GetGrid()->GetAxis(0);
1256
1257     if (axis->GetNbins() < measured->GetNbinsX())
1258     {
1259       if (axis2 != -1)
1260       {
1261         // 2d rebin with variable axis does not exist in root
1262
1263         TH1* tmp = measured;
1264         measured = new TH2D(Form("%s_rebinned", tmp->GetName()), tmp->GetTitle(), axis->GetNbins(), axis->GetXbins()->GetArray(), tmp->GetNbinsY(), tmp->GetYaxis()->GetXbins()->GetArray());
1265         for (Int_t x=1; x<=tmp->GetNbinsX(); x++)
1266           for (Int_t y=1; y<=tmp->GetNbinsY(); y++)
1267           {
1268             ((TH2*) measured)->Fill(tmp->GetXaxis()->GetBinCenter(x), tmp->GetYaxis()->GetBinCenter(y), tmp->GetBinContent(x, y));
1269             measured->SetBinError(x, y, 0); // cannot trust bin error, set to 0
1270           }
1271         delete tmp;
1272
1273         tmp = generated;
1274         generated = new TH2D(Form("%s_rebinned", tmp->GetName()), tmp->GetTitle(), axis->GetNbins(), axis->GetXbins()->GetArray(), tmp->GetNbinsY(), tmp->GetYaxis()->GetXbins()->GetArray());
1275         for (Int_t x=1; x<=tmp->GetNbinsX(); x++)
1276           for (Int_t y=1; y<=tmp->GetNbinsY(); y++)
1277           {
1278             ((TH2*) generated)->Fill(tmp->GetXaxis()->GetBinCenter(x), tmp->GetYaxis()->GetBinCenter(y), tmp->GetBinContent(x, y));
1279             generated->SetBinError(x, y, 0); // cannot trust bin error, set to 0
1280           }
1281         delete tmp;
1282       }
1283       else
1284       {
1285         TH1* tmp = measured;
1286         measured = tmp->Rebin(axis->GetNbins(), Form("%s_rebinned", tmp->GetName()), axis->GetXbins()->GetArray());
1287         delete tmp;
1288
1289         tmp = generated;
1290         generated = tmp->Rebin(axis->GetNbins(), Form("%s_rebinned", tmp->GetName()), axis->GetXbins()->GetArray());
1291         delete tmp;
1292       }
1293     }
1294     else if (axis->GetNbins() > measured->GetNbinsX())
1295     {
1296       if (axis2 != -1)
1297         AliFatal("Rebinning only works for 1d at present");
1298
1299       // this is an unfortunate case where the number of bins has to be increased in principle
1300       // there is a region where the binning is finner in one histogram and a region where it is the other way round
1301       // this cannot be resolved in principle, but as we only calculate the ratio the bin in the second region get the same entries
1302       // only a certain binning is supported here
1303       if (axis->GetNbins() != 100 || measured->GetNbinsX() != 39)
1304         AliFatal(Form("Invalid binning --> %d %d", axis->GetNbins(), measured->GetNbinsX()));
1305
1306       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};
1307
1308       // reduce binning below 5 GeV/c
1309       TH1* tmp = measured;
1310       measured = tmp->Rebin(27, Form("%s_rebinned", tmp->GetName()), newBins);
1311       delete tmp;
1312
1313       // increase binning above 5 GeV/c
1314       tmp = measured;
1315       measured = new TH1F(Form("%s_rebinned2", tmp->GetName()), tmp->GetTitle(), axis->GetNbins(), axis->GetBinLowEdge(1), axis->GetBinUpEdge(axis->GetNbins()));
1316       for (Int_t bin=1; bin<=measured->GetNbinsX(); bin++)
1317       {
1318         measured->SetBinContent(bin, tmp->GetBinContent(tmp->FindBin(measured->GetBinCenter(bin))));
1319         measured->SetBinError(bin, tmp->GetBinError(tmp->FindBin(measured->GetBinCenter(bin))));
1320       }
1321       delete tmp;
1322
1323       // reduce binning below 5 GeV/c
1324       tmp = generated;
1325       generated = tmp->Rebin(27, Form("%s_rebinned", tmp->GetName()), newBins);
1326       delete tmp;
1327
1328       // increase binning above 5 GeV/c
1329       tmp = generated;
1330       generated = new TH1F(Form("%s_rebinned2", tmp->GetName()), tmp->GetTitle(), axis->GetNbins(), axis->GetBinLowEdge(1), axis->GetBinUpEdge(axis->GetNbins()));
1331       for (Int_t bin=1; bin<=generated->GetNbinsX(); bin++)
1332       {
1333         generated->SetBinContent(bin, tmp->GetBinContent(tmp->FindBin(generated->GetBinCenter(bin))));
1334         generated->SetBinError(bin, tmp->GetBinError(tmp->FindBin(generated->GetBinCenter(bin))));
1335       }
1336       delete tmp;
1337     }
1338   }
1339
1340   measured->Divide(measured, generated, 1, 1, "B");
1341
1342   delete generated;
1343
1344   ResetBinLimits(sourceContainer->GetGrid(step1));
1345   ResetBinLimits(sourceContainer->GetGrid(step2));
1346
1347   return measured;
1348 }
1349
1350 //____________________________________________________________________
1351 TH1* AliUEHist::GetEventEfficiency(CFStep step1, CFStep step2, Int_t axis1, Int_t axis2, Float_t ptLeadMin, Float_t ptLeadMax)
1352 {
1353   // creates a event-level efficiency by dividing step2 by step1
1354   // projected to axis1 and axis2 (optional if >= 0)
1355
1356   if (ptLeadMax > ptLeadMin)
1357   {
1358     fEventHist->GetGrid(step1)->SetRangeUser(0, ptLeadMin, ptLeadMax);
1359     fEventHist->GetGrid(step2)->SetRangeUser(0, ptLeadMin, ptLeadMax);
1360   }
1361
1362   TH1* measured = 0;
1363   TH1* generated = 0;
1364
1365   if (axis2 >= 0)
1366   {
1367     generated = fEventHist->Project(step1, axis1, axis2);
1368     measured = fEventHist->Project(step2, axis1, axis2);
1369   }
1370   else
1371   {
1372     generated = fEventHist->Project(step1, axis1);
1373     measured = fEventHist->Project(step2, axis1);
1374   }
1375
1376   measured->Divide(measured, generated, 1, 1, "B");
1377
1378   delete generated;
1379
1380   if (ptLeadMax > ptLeadMin)
1381   {
1382     fEventHist->GetGrid(step1)->SetRangeUser(0, 0, -1);
1383     fEventHist->GetGrid(step2)->SetRangeUser(0, 0, -1);
1384   }
1385
1386   return measured;
1387 }
1388
1389 //____________________________________________________________________
1390 void AliUEHist::WeightHistogram(TH3* hist1, TH1* hist2)
1391 {
1392   // weights each entry of the 3d histogram hist1 with the 1d histogram hist2
1393   // where the matching is done of the z axis of hist1 with the x axis of hist2
1394
1395   if (hist1->GetNbinsZ() != hist2->GetNbinsX())
1396     AliFatal(Form("Inconsistent binning %d %d", hist1->GetNbinsZ(), hist2->GetNbinsX()));
1397
1398   for (Int_t x=1; x<=hist1->GetNbinsX(); x++)
1399   {
1400     for (Int_t y=1; y<=hist1->GetNbinsY(); y++)
1401     {
1402       for (Int_t z=1; z<=hist1->GetNbinsZ(); z++)
1403       {
1404         if (hist2->GetBinContent(z) > 0)
1405         {
1406           hist1->SetBinContent(x, y, z, hist1->GetBinContent(x, y, z) / hist2->GetBinContent(z));
1407           hist1->SetBinError(x, y, z, hist1->GetBinError(x, y, z) / hist2->GetBinContent(z));
1408         }
1409         else
1410         {
1411           hist1->SetBinContent(x, y, z, 0);
1412           hist1->SetBinError(x, y, z, 0);
1413         }
1414       }
1415     }
1416   }
1417 }
1418
1419 //____________________________________________________________________
1420 TH1* AliUEHist::GetBias(CFStep step1, CFStep step2, Int_t region, const char* axis, Float_t leadPtMin, Float_t leadPtMax, Int_t weighting)
1421 {
1422   // extracts the track-level bias (integrating out the multiplicity) between two steps (dividing step2 by step1)
1423   // in the given region (sum over all regions is calculated if region == -1)
1424   // done by weighting the track-level distribution with the number of events as function of leading pT
1425   // and then calculating the ratio between the distributions
1426   // projected to axis which is a TH3::Project3D string, e.g. "x", or "yx"
1427   //   no projection is done if axis == 0
1428   // weighting: 0 = tracks weighted with events (as discussed above)
1429   //            1 = only track bias is returned
1430   //            2 = only event bias is returned
1431
1432   AliCFContainer* tmp = 0;
1433
1434   if (region == -1)
1435   {
1436     tmp = (AliCFContainer*) fTrackHist[0]->Clone();
1437     for (Int_t i = 1; i < fkRegions; i++)
1438       if (fTrackHist[i])
1439         tmp->Add(fTrackHist[i]);
1440   }
1441   else if (region == kMin && fCombineMinMax)
1442   {
1443     tmp = (AliCFContainer*) fTrackHist[kMin]->Clone();
1444     tmp->Add(fTrackHist[kMax]);
1445   }
1446   else
1447     tmp = fTrackHist[region];
1448
1449   ResetBinLimits(tmp->GetGrid(step1));
1450   ResetBinLimits(fEventHist->GetGrid(step1));
1451   SetBinLimits(tmp->GetGrid(step1));
1452
1453   ResetBinLimits(tmp->GetGrid(step2));
1454   ResetBinLimits(fEventHist->GetGrid(step2));
1455   SetBinLimits(tmp->GetGrid(step2));
1456
1457   TH1D* events1 = (TH1D*)fEventHist->Project(step1, 0);
1458   TH3D* hist1 = (TH3D*)tmp->Project(step1, 0, tmp->GetNVar()-1, 2);
1459   if (weighting == 0)
1460     WeightHistogram(hist1, events1);
1461
1462   TH1D* events2 = (TH1D*)fEventHist->Project(step2, 0);
1463   TH3D* hist2 = (TH3D*)tmp->Project(step2, 0, tmp->GetNVar()-1, 2);
1464   if (weighting == 0)
1465     WeightHistogram(hist2, events2);
1466
1467   TH1* generated = hist1;
1468   TH1* measured = hist2;
1469
1470   if (weighting == 0 || weighting == 1)
1471   {
1472     if (axis)
1473     {
1474       if (leadPtMax > leadPtMin)
1475       {
1476         hist1->GetZaxis()->SetRangeUser(leadPtMin, leadPtMax);
1477         hist2->GetZaxis()->SetRangeUser(leadPtMin, leadPtMax);
1478       }
1479
1480       if (fEtaMax > fEtaMin && !TString(axis).Contains("x"))
1481       {
1482         hist1->GetXaxis()->SetRangeUser(fEtaMin, fEtaMax);
1483         hist2->GetXaxis()->SetRangeUser(fEtaMin, fEtaMax);
1484       }
1485
1486       generated = hist1->Project3D(axis);
1487       measured  = hist2->Project3D(axis);
1488
1489       // delete hists here if projection has been done
1490       delete hist1;
1491       delete hist2;
1492     }
1493     delete events1;
1494     delete events2;
1495   }
1496   else if (weighting == 2)
1497   {
1498     delete hist1;
1499     delete hist2;
1500     generated = events1;
1501     measured = events2;
1502   }
1503
1504   measured->Divide(generated);
1505
1506   delete generated;
1507
1508   ResetBinLimits(tmp->GetGrid(step1));
1509   ResetBinLimits(tmp->GetGrid(step2));
1510
1511   if ((region == -1) || (region == kMin && fCombineMinMax))
1512     delete tmp;
1513
1514   return measured;
1515 }
1516
1517 //____________________________________________________________________
1518 void AliUEHist::CorrectCorrelatedContamination(CFStep step, Int_t region, TH1* trackCorrection)
1519 {
1520   // corrects for the given factor in a small delta-eta and delta-phi window as function of pT,A and pT,T
1521
1522   if (!fTrackHist[region])
1523     return;
1524
1525   THnSparse* grid = fTrackHist[region]->GetGrid(step)->GetGrid();
1526
1527   Int_t var1 = 1;
1528   Int_t var2 = 2;
1529
1530   if (grid->GetAxis(var1)->GetNbins() != trackCorrection->GetNbinsX())
1531     AliFatal(Form("Invalid binning (var1): %d %d", grid->GetAxis(var1)->GetNbins(), trackCorrection->GetNbinsX()));
1532
1533   if (grid->GetAxis(var2)->GetNbins() != trackCorrection->GetNbinsY())
1534     AliFatal(Form("Invalid binning (var2): %d %d", grid->GetAxis(var2)->GetNbins(), trackCorrection->GetNbinsY()));
1535
1536   // optimized implementation
1537   for (Int_t binIdx = 0; binIdx < grid->GetNbins(); binIdx++)
1538   {
1539     Int_t bins[5];
1540
1541     Double_t value = grid->GetBinContent(binIdx, bins);
1542     Double_t error = grid->GetBinError(binIdx);
1543
1544     // check eta and phi axes
1545     if (TMath::Abs(grid->GetAxis(0)->GetBinCenter(bins[0])) > 0.1)
1546       continue;
1547     if (TMath::Abs(grid->GetAxis(4)->GetBinCenter(bins[4])) > 0.1)
1548       continue;
1549
1550     value *= trackCorrection->GetBinContent(bins[var1], bins[var2]);
1551     error *= trackCorrection->GetBinContent(bins[var1], bins[var2]);
1552
1553     grid->SetBinContent(bins, value);
1554     grid->SetBinError(bins, error);
1555   }
1556
1557   Printf("AliUEHist::CorrectCorrelatedContamination: Corrected.");
1558 }
1559
1560 //____________________________________________________________________
1561 TH2* AliUEHist::GetCorrelatedContamination()
1562 {
1563   // 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!)
1564
1565   Int_t step1 = kCFStepTrackedOnlyPrim;
1566   Int_t step2 = kCFStepTracked;
1567
1568   fTrackHist[0]->GetGrid(step1)->SetRangeUser(0, -0.01, 0.01); // delta eta
1569   fTrackHist[0]->GetGrid(step1)->SetRangeUser(4, -0.01, 0.01); // delta phi
1570   TH2* tracksStep1 = (TH2*) fTrackHist[0]->Project(step1, 1, 2);
1571
1572   fTrackHist[0]->GetGrid(step2)->SetRangeUser(0, -0.01, 0.01); // delta eta
1573   fTrackHist[0]->GetGrid(step2)->SetRangeUser(4, -0.01, 0.01); // delta phi
1574   TH2* tracksStep2 = (TH2*) fTrackHist[0]->Project(step2, 1, 2);
1575
1576   tracksStep1->Divide(tracksStep2);
1577
1578   TH1* triggersStep1 = fEventHist->Project(step1, 0);
1579   TH1* triggersStep2 = fEventHist->Project(step2, 0);
1580
1581   TH1* singleParticle = GetTrackingContamination(1);
1582
1583   for (Int_t x=1; x<=tracksStep1->GetNbinsX(); x++)
1584     for (Int_t y=1; y<=tracksStep1->GetNbinsY(); y++)
1585       if (singleParticle->GetBinContent(x) > 0)
1586         tracksStep1->SetBinContent(x, y, tracksStep1->GetBinContent(x, y) / triggersStep1->GetBinContent(y) * triggersStep2->GetBinContent(y) / singleParticle->GetBinContent(x));
1587       else
1588         tracksStep1->SetBinContent(x, y, 0);
1589
1590   delete singleParticle;
1591   delete tracksStep2;
1592   delete triggersStep1;
1593   delete triggersStep2;
1594
1595   return tracksStep1;
1596 }
1597
1598 //____________________________________________________________________
1599 TH2D* AliUEHist::GetTrackingEfficiency()
1600 {
1601   // extracts the tracking efficiency by calculating the efficiency from step kCFStepAnaTopology to kCFStepTrackedOnlyPrim
1602   // integrates over the regions and all other variables than pT and eta to increase the statistics
1603   //
1604   // returned histogram has to be deleted by the user
1605
1606   return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepAnaTopology, kCFStepTrackedOnlyPrim, 0, 1));
1607 }
1608
1609 //____________________________________________________________________
1610 TH2D* AliUEHist::GetTrackingEfficiencyCentrality()
1611 {
1612   // extracts the tracking efficiency by calculating the efficiency from step kCFStepAnaTopology to kCFStepTrackedOnlyPrim
1613   // integrates over the regions and all other variables than pT, centrality to increase the statistics
1614   //
1615   // returned histogram has to be deleted by the user
1616
1617   return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepAnaTopology, kCFStepTrackedOnlyPrim, 1, 3));
1618 }
1619
1620 //____________________________________________________________________
1621 TH1D* AliUEHist::GetTrackingEfficiency(Int_t axis)
1622 {
1623   // extracts the tracking efficiency by calculating the efficiency from step kCFStepAnaTopology to kCFStepTrackedOnlyPrim
1624   // integrates over the regions and all other variables than pT (axis == 0) and eta (axis == 1) to increase the statistics
1625
1626   return dynamic_cast<TH1D*> (GetTrackEfficiency(kCFStepAnaTopology, kCFStepTrackedOnlyPrim, axis));
1627 }
1628
1629 //____________________________________________________________________
1630 TH2D* AliUEHist::GetTrackingCorrection()
1631 {
1632   // extracts the tracking correction by calculating the efficiency from step kCFStepAnaTopology to kCFStepTracked
1633   // integrates over the regions and all other variables than pT and eta to increase the statistics
1634   //
1635   // returned histogram has to be deleted by the user
1636
1637   return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepTracked, kCFStepAnaTopology, 0, 1));
1638 }
1639
1640 //____________________________________________________________________
1641 TH1D* AliUEHist::GetTrackingCorrection(Int_t axis)
1642 {
1643   // extracts the tracking correction by calculating the efficiency from step kCFStepAnaTopology to kCFStepTracked
1644   // integrates over the regions and all other variables than pT (axis == 0) and eta (axis == 1) to increase the statistics
1645
1646   return dynamic_cast<TH1D*> (GetTrackEfficiency(kCFStepTracked, kCFStepAnaTopology, axis));
1647 }
1648
1649 //____________________________________________________________________
1650 TH2D* AliUEHist::GetTrackingEfficiencyCorrection()
1651 {
1652   // extracts the tracking correction by calculating the efficiency from step kCFStepAnaTopology to kCFStepTracked
1653   // integrates over the regions and all other variables than pT and eta to increase the statistics
1654   //
1655   // returned histogram has to be deleted by the user
1656
1657   return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, 0, 1));
1658 }
1659
1660 //____________________________________________________________________
1661 TH2D* AliUEHist::GetTrackingEfficiencyCorrectionCentrality()
1662 {
1663   // extracts the tracking correction by calculating the efficiency from step kCFStepAnaTopology to kCFStepTracked
1664   // integrates over the regions and all other variables than pT and centrality to increase the statistics
1665   //
1666   // returned histogram has to be deleted by the user
1667
1668   return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, 1, 3));
1669 }
1670
1671 //____________________________________________________________________
1672 TH1D* AliUEHist::GetTrackingEfficiencyCorrection(Int_t axis)
1673 {
1674   // extracts the tracking correction by calculating the efficiency from step kCFStepAnaTopology to kCFStepTracked
1675   // integrates over the regions and all other variables than pT (axis == 0) and eta (axis == 1) to increase the statistics
1676
1677   return dynamic_cast<TH1D*> (GetTrackEfficiency(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, axis));
1678 }
1679
1680 //____________________________________________________________________
1681 TH2D* AliUEHist::GetTrackingContamination()
1682 {
1683   // extracts the tracking contamination by secondaries by calculating the efficiency from step kCFStepTrackedOnlyPrim to kCFStepTracked
1684   // integrates over the regions and all other variables than pT and eta to increase the statistics
1685   //
1686   // returned histogram has to be deleted by the user
1687
1688   return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepTracked, kCFStepTrackedOnlyPrim, 0, 1));
1689 }
1690
1691 //____________________________________________________________________
1692 TH2D* AliUEHist::GetTrackingContaminationCentrality()
1693 {
1694   // extracts the tracking contamination by secondaries by calculating the efficiency from step kCFStepTrackedOnlyPrim to kCFStepTracked
1695   // integrates over the regions and all other variables than pT and centrality to increase the statistics
1696   //
1697   // returned histogram has to be deleted by the user
1698
1699   return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepTracked, kCFStepTrackedOnlyPrim, 1, 3));
1700 }
1701
1702 //____________________________________________________________________
1703 TH1D* AliUEHist::GetTrackingContamination(Int_t axis)
1704 {
1705   // extracts the tracking contamination by secondaries by calculating the efficiency from step kCFStepTrackedOnlyPrim to kCFStepTracked
1706   // integrates over the regions and all other variables than pT (axis == 0) and eta (axis == 1) to increase the statistics
1707
1708   return dynamic_cast<TH1D*> (GetTrackEfficiency(kCFStepTracked, kCFStepTrackedOnlyPrim, axis));
1709 }
1710
1711 //____________________________________________________________________
1712 const char* AliUEHist::GetRegionTitle(Region region)
1713 {
1714   // returns the name of the given region
1715
1716   switch (region)
1717   {
1718     case kToward:
1719       return "Towards";
1720     case kAway:
1721       return "Away";
1722     case kMin:
1723       return (fCombineMinMax) ? "Transverse" : "Min";
1724     case kMax:
1725       return "Max";
1726   }
1727
1728   return 0;
1729 }
1730
1731 //____________________________________________________________________
1732 const char* AliUEHist::GetStepTitle(CFStep step)
1733 {
1734   // returns the name of the given step
1735
1736   switch (step)
1737   {
1738     case kCFStepAll:
1739       return "All events";
1740     case kCFStepTriggered:
1741       return "Triggered";
1742     case kCFStepVertex:
1743       return "Primary Vertex";
1744     case kCFStepAnaTopology:
1745       return "Required analysis topology";
1746     case kCFStepTrackedOnlyPrim:
1747       return "Tracked (matched MC, only primaries)";
1748     case kCFStepTracked:
1749       return "Tracked (matched MC, all)";
1750     case kCFStepReconstructed:
1751       return "Reconstructed";
1752     case kCFStepRealLeading:
1753       return "Correct leading particle identified";
1754     case kCFStepBiasStudy:
1755       return "Bias study applying tracking efficiency";
1756     case kCFStepBiasStudy2:
1757       return "Bias study applying tracking efficiency in two steps";
1758   }
1759
1760   return 0;
1761 }
1762
1763 //____________________________________________________________________
1764 void AliUEHist::CopyReconstructedData(AliUEHist* from)
1765 {
1766   // copies those histograms extracted from ESD to this object
1767
1768   // TODO at present only the pointers are copied
1769
1770   for (Int_t region=0; region<4; region++)
1771   {
1772     if (!fTrackHist[region])
1773       continue;
1774
1775     fTrackHist[region]->SetGrid(AliUEHist::kCFStepReconstructed, from->fTrackHist[region]->GetGrid(AliUEHist::kCFStepReconstructed));
1776     //fTrackHist[region]->SetGrid(AliUEHist::kCFStepTrackedOnlyPrim, from->fTrackHist[region]->GetGrid(AliUEHist::kCFStepTrackedOnlyPrim));
1777     fTrackHist[region]->SetGrid(AliUEHist::kCFStepBiasStudy,     from->fTrackHist[region]->GetGrid(AliUEHist::kCFStepBiasStudy));
1778   }
1779
1780   fEventHist->SetGrid(AliUEHist::kCFStepReconstructed, from->fEventHist->GetGrid(AliUEHist::kCFStepReconstructed));
1781   //fEventHist->SetGrid(AliUEHist::kCFStepTrackedOnlyPrim, from->fEventHist->GetGrid(AliUEHist::kCFStepTrackedOnlyPrim));
1782   fEventHist->SetGrid(AliUEHist::kCFStepBiasStudy,     from->fEventHist->GetGrid(AliUEHist::kCFStepBiasStudy));
1783 }
1784
1785 //____________________________________________________________________
1786 void AliUEHist::ExtendTrackingEfficiency(Bool_t verbose)
1787 {
1788   // fits the tracking efficiency at high pT with a constant and fills all bins with this tracking efficiency
1789
1790   Float_t fitRangeBegin = 5.01;
1791   Float_t fitRangeEnd = 14.99;
1792   Float_t extendRangeBegin = 10.01;
1793
1794   if (fTrackHistEfficiency->GetNVar() == 3)
1795   {
1796     TH1* obj = GetTrackingEfficiency(1);
1797
1798     if (verbose)
1799     {
1800       new TCanvas;
1801       obj->Draw();
1802     }
1803
1804     obj->Fit("pol0", (verbose) ? "+" : "0+", "SAME", fitRangeBegin, fitRangeEnd);
1805
1806     Float_t trackingEff = obj->GetFunction("pol0")->GetParameter(0);
1807
1808     obj = GetTrackingContamination(1);
1809
1810     if (verbose)
1811     {
1812       new TCanvas;
1813       obj->Draw();
1814     }
1815
1816     obj->Fit("pol0", (verbose) ? "+" : "0+", "SAME", fitRangeBegin, fitRangeEnd);
1817
1818     Float_t trackingCont = obj->GetFunction("pol0")->GetParameter(0);
1819
1820     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);
1821
1822     // extend for full pT range
1823     for (Int_t x = fTrackHistEfficiency->GetAxis(0, 0)->FindBin(fEtaMin); x <= fTrackHistEfficiency->GetAxis(0, 0)->FindBin(fEtaMax); x++)
1824       for (Int_t y = fTrackHistEfficiency->GetAxis(1, 0)->FindBin(extendRangeBegin); y <= fTrackHistEfficiency->GetNBins(1); y++)
1825         for (Int_t z = 1; z <= fTrackHistEfficiency->GetNBins(2); z++) // particle type axis
1826         {
1827
1828           Int_t bins[3];
1829           bins[0] = x;
1830           bins[1] = y;
1831           bins[2] = z;
1832
1833           fTrackHistEfficiency->GetGrid(0)->SetElement(bins, 100);
1834           fTrackHistEfficiency->GetGrid(1)->SetElement(bins, 100.0 * trackingEff);
1835           fTrackHistEfficiency->GetGrid(2)->SetElement(bins, 100.0 * trackingEff / trackingCont);
1836         }
1837   }
1838   else if (fTrackHistEfficiency->GetNVar() == 4)
1839   {
1840     // fit in centrality intervals of 20% for efficiency, one bin for contamination
1841     Float_t* trackingEff = 0;
1842     Float_t* trackingCont = 0;
1843     Int_t centralityBins = 5;
1844
1845     Printf("AliUEHist::ExtendTrackingEfficiency: Fitting efficiencies between %f and %f. Extending from %f onwards (within %f < eta < %f)", fitRangeBegin, fitRangeEnd, extendRangeBegin, fEtaMin, fEtaMax);
1846
1847     // 0 = eff; 1 = cont
1848     for (Int_t caseNo = 0; caseNo < 2; caseNo++)
1849     {
1850       Float_t* target = 0;
1851       Int_t centralityBinsLocal = centralityBins;
1852
1853       if (caseNo == 0)
1854       {
1855         trackingEff = new Float_t[centralityBinsLocal];
1856         target = trackingEff;
1857       }
1858       else
1859       {
1860         centralityBinsLocal = 1;
1861         trackingCont = new Float_t[centralityBinsLocal];
1862         target = trackingCont;
1863       }
1864
1865       for (Int_t i=0; i<centralityBinsLocal; i++)
1866       {
1867         SetCentralityRange(100.0/centralityBinsLocal*i + 0.1, 100.0/centralityBinsLocal*(i+1) - 0.1);
1868         TH1* proj = (caseNo == 0) ? GetTrackingEfficiency(1) : GetTrackingContamination(1);
1869         if (verbose)
1870         {
1871           new TCanvas;
1872           proj->DrawCopy();
1873         }
1874         if ((Int_t) proj->Fit("pol0", (verbose) ? "+" : "Q0+", "SAME", fitRangeBegin, fitRangeEnd) == 0)
1875           target[i] = proj->GetFunction("pol0")->GetParameter(0);
1876         else
1877           target[i] = 0;
1878         Printf("AliUEHist::ExtendTrackingEfficiency: case %d, centrality %d, eff %f", caseNo, i, target[i]);
1879       }
1880     }
1881
1882     // extend for full pT range
1883     for (Int_t x = fTrackHistEfficiency->GetAxis(0, 0)->FindBin(fEtaMin); x <= fTrackHistEfficiency->GetAxis(0, 0)->FindBin(fEtaMax); x++)
1884       for (Int_t y = fTrackHistEfficiency->GetAxis(1, 0)->FindBin(extendRangeBegin); y <= fTrackHistEfficiency->GetNBins(1); y++)
1885         for (Int_t z = 1; z <= fTrackHistEfficiency->GetNBins(2); z++) // particle type axis
1886         {
1887           for (Int_t z2 = 1; z2 <= fTrackHistEfficiency->GetNBins(3); z2++) // centrality axis
1888           {
1889
1890             Int_t bins[4];
1891             bins[0] = x;
1892             bins[1] = y;
1893             bins[2] = z;
1894             bins[3] = z2;
1895
1896             Int_t z2Bin = (Int_t) (fTrackHistEfficiency->GetAxis(3, 0)->GetBinCenter(z2) / (100.0 / centralityBins));
1897             //Printf("%d %d", z2, z2Bin);
1898
1899             fTrackHistEfficiency->GetGrid(0)->SetElement(bins, 100);
1900             fTrackHistEfficiency->GetGrid(1)->SetElement(bins, 100.0 * trackingEff[z2Bin]);
1901             if (trackingCont[0] > 0)
1902               fTrackHistEfficiency->GetGrid(2)->SetElement(bins, 100.0 * trackingEff[z2Bin] / trackingCont[0]);
1903             else
1904               fTrackHistEfficiency->GetGrid(2)->SetElement(bins, 0);
1905           }
1906        }
1907
1908      delete[] trackingEff;
1909      delete[] trackingCont;
1910    }
1911
1912    SetCentralityRange(0, 0);
1913 }
1914
1915 void AliUEHist::AdditionalDPhiCorrection(Int_t step)
1916 {
1917   // corrects the dphi distribution with an extra factor close to dphi ~ 0
1918
1919   Printf("WARNING: In AliUEHist::AdditionalDPhiCorrection.");
1920
1921   THnSparse* grid = fTrackHist[0]->GetGrid(step)->GetGrid();
1922
1923   // optimized implementation
1924   for (Int_t binIdx = 0; binIdx < grid->GetNbins(); binIdx++)
1925   {
1926     Int_t bins[5];
1927     Double_t value = grid->GetBinContent(binIdx, bins);
1928     Double_t error = grid->GetBinError(binIdx);
1929
1930     Float_t binCenter = grid->GetAxis(4)->GetBinCenter(bins[4]);
1931     if (TMath::Abs(binCenter) < 0.2)
1932     {
1933       value *= 0.985;
1934       error *= 0.985;
1935     }
1936     else if (TMath::Abs(binCenter) < 0.3)
1937     {
1938       value *= 0.9925;
1939       error *= 0.9925;
1940     }
1941
1942     grid->SetBinContent(bins, value);
1943     grid->SetBinError(bins, error);
1944   }
1945 }
1946
1947 void AliUEHist::Scale(Double_t factor)
1948 {
1949   // scales all contained histograms by the given factor
1950
1951   for (Int_t i=0; i<4; i++)
1952     if (fTrackHist[i])
1953       fTrackHist[i]->Scale(factor);
1954
1955   fEventHist->Scale(factor);
1956   fTrackHistEfficiency->Scale(factor);
1957 }
1958
1959 void AliUEHist::Reset()
1960 {
1961   // resets all contained histograms
1962
1963   for (Int_t i=0; i<4; i++)
1964     if (fTrackHist[i])
1965       for (Int_t step=0; step<fTrackHist[i]->GetNStep(); step++)
1966         fTrackHist[i]->GetGrid(step)->GetGrid()->Reset();
1967
1968   for (Int_t step=0; step<fEventHist->GetNStep(); step++)
1969     fEventHist->GetGrid(step)->GetGrid()->Reset();
1970
1971   for (Int_t step=0; step<fTrackHistEfficiency->GetNStep(); step++)
1972     fTrackHistEfficiency->GetGrid(step)->GetGrid()->Reset();
1973 }