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