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 //____________________________________________________________________
 783 void AliUEHist::CorrectTracks(CFStep step1, CFStep step2, TH1* trackCorrection, Int_t var1, Int_t var2)
 784 {
 785   // corrects from step1 to step2 by multiplying the tracks with trackCorrection
 786   // trackCorrection can be a function of eta (var1 == 0), pT (var1 == 1), leading pT (var1 == 2), multiplicity (var1 == 3), delta phi (var1 == 4)
 787   // if var2 >= 0 a two dimension correction is assumed in trackCorrection
 788   //
 789   // if trackCorrection is 0, just copies content from step1 to step2
 790
 791   for (UInt_t region=0; region<fkRegions; region++)
 792     CorrectTracks(step1, step2, region, trackCorrection, var1, var2);
 793 }
 794
 795 //____________________________________________________________________
 796 void AliUEHist::CorrectTracks(CFStep step1, CFStep step2, Int_t region, TH1* trackCorrection, Int_t var1, Int_t var2)
 797 {
 798   //
 799   // see documentation of CorrectTracks above
 800   //
 801
 802   if (!fTrackHist[region])
 803     return;
 804
 805   THnSparse* grid = fTrackHist[region]->GetGrid(step1)->GetGrid();
 806   THnSparse* target = fTrackHist[region]->GetGrid(step2)->GetGrid();
 807
 808   // clear target histogram
 809   target->Reset();
 810
 811   if (trackCorrection != 0)
 812   {
 813     if (grid->GetAxis(var1)->GetNbins() != trackCorrection->GetNbinsX())
 814       AliFatal(Form("Invalid binning (var1): %d %d", grid->GetAxis(var1)->GetNbins(), trackCorrection->GetNbinsX()));
 815
 816     if (var2 >= 0 && grid->GetAxis(var2)->GetNbins() != trackCorrection->GetNbinsY())
 817       AliFatal(Form("Invalid binning (var2): %d %d", grid->GetAxis(var2)->GetNbins(), trackCorrection->GetNbinsY()));
 818   }
 819
 820   // optimized implementation
 821   for (Int_t binIdx = 0; binIdx < grid->GetNbins(); binIdx++)
 822   {
 823     Int_t bins[5];
 824     Double_t value = grid->GetBinContent(binIdx, bins);
 825     Double_t error = grid->GetBinError(binIdx);
 826
 827     if (trackCorrection != 0)
 828     {
 829       if (var2 < 0)
 830       {
 831         value *= trackCorrection->GetBinContent(bins[var1]);
 832         error *= trackCorrection->GetBinContent(bins[var1]);
 833       }
 834       else
 835       {
 836         value *= trackCorrection->GetBinContent(bins[var1], bins[var2]);
 837         error *= trackCorrection->GetBinContent(bins[var1], bins[var2]);
 838       }
 839     }
 840
 841     target->SetBinContent(bins, value);
 842     target->SetBinError(bins, error);
 843   }
 844
 845   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);
 846 }
 847
 848 //____________________________________________________________________
 849 void AliUEHist::CorrectEvents(CFStep step1, CFStep step2, TH1* eventCorrection, Int_t var1, Int_t var2)
 850 {
 851   // corrects from step1 to step2 by multiplying the events with eventCorrection
 852   // eventCorrection is as function of leading pT (var == 0) or multiplicity (var == 1)
 853   //
 854   // if eventCorrection is 0, just copies content from step1 to step2
 855
 856   AliCFGridSparse* grid = fEventHist->GetGrid(step1);
 857   AliCFGridSparse* target = fEventHist->GetGrid(step2);
 858
 859   // clear target histogram
 860   target->GetGrid()->Reset();
 861
 862   if (eventCorrection != 0 && grid->GetNBins(var1) != eventCorrection->GetNbinsX())
 863     AliFatal(Form("Invalid binning: %d %d", grid->GetNBins(var1), eventCorrection->GetNbinsX()));
 864
 865   if (eventCorrection != 0 && var2 != -1 && grid->GetNBins(var2) != eventCorrection->GetNbinsY())
 866     AliFatal(Form("Invalid binning: %d %d", grid->GetNBins(var2), eventCorrection->GetNbinsY()));
 867
 868   Int_t bins[2];
 869   for (Int_t x = 1; x <= grid->GetNBins(0); x++)
 870   {
 871     bins[0] = x;
 872     for (Int_t y = 1; y <= grid->GetNBins(1); y++)
 873     {
 874       bins[1] = y;
 875
 876       Double_t value = grid->GetElement(bins);
 877       if (value != 0)
 878       {
 879         Double_t error = grid->GetElementError(bins);
 880
 881         if (eventCorrection != 0)
 882         {
 883           if (var2 == -1)
 884           {
 885             value *= eventCorrection->GetBinContent(bins[var1]);
 886             error *= eventCorrection->GetBinContent(bins[var1]);
 887           }
 888           else
 889           {
 890             value *= eventCorrection->GetBinContent(bins[var1], bins[var2]);
 891             error *= eventCorrection->GetBinContent(bins[var1], bins[var2]);
 892           }
 893         }
 894
 895         target->SetElement(bins, value);
 896         target->SetElementError(bins, error);
 897       }
 898     }
 899   }
 900
 901   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);
 902 }
 903
 904 //____________________________________________________________________
 905 void AliUEHist::Correct(AliUEHist* corrections)
 906 {
 907   // applies the given corrections to extract from the step kCFStepReconstructed all previous steps
 908   //
 909   // in this object the data is expected in the step kCFStepReconstructed
 910
 911   if (fHistogramType.Length() == 0)
 912   {
 913     Printf("W-AliUEHist::Correct: fHistogramType not defined. Guessing histogram type...");
 914     if (fTrackHist[kToward]->GetNVar() < 5)
 915     {
 916       if (strcmp(fTrackHist[kToward]->GetTitle(), "d^{2}N_{ch}/d#phid#eta") == 0)
 917         fHistogramType = "NumberDensitypT";
 918       else if (strcmp(fTrackHist[kToward]->GetTitle(), "d^{2}#Sigma p_{T}/d#phid#eta") == 0)
 919         fHistogramType = "SumpT";
 920     }
 921     else if (fTrackHist[kToward]->GetNVar() == 5)
 922     {
 923       if (strcmp(fTrackHist[kToward]->GetGrid(0)->GetGrid()->GetAxis(3)->GetTitle(), "multiplicity") == 0)
 924         fHistogramType = "NumberDensityPhi";
 925       else if (strcmp(fTrackHist[kToward]->GetGrid(0)->GetGrid()->GetAxis(3)->GetTitle(), "centrality") == 0)
 926         fHistogramType = "NumberDensityPhiCentrality";
 927     }
 928
 929     if (fHistogramType.Length() == 0)
 930       AliFatal("Cannot figure out histogram type. Try setting it manually...");
 931   }
 932
 933   Printf("AliUEHist::Correct: Correcting %s...", fHistogramType.Data());
 934
 935   if (strcmp(fHistogramType, "NumberDensitypT") == 0 || strcmp(fHistogramType, "NumberDensityPhi") == 0 || strcmp(fHistogramType, "SumpT") == 0) // the last is for backward compatibilty
 936   {
 937     // ---- track level
 938
 939     // bias due to migration in leading pT (because the leading particle is not reconstructed, and the subleading is used)
 940     // extracted as function of leading pT
 941     Bool_t biasFromMC = kFALSE;
 942     const char* projAxis = "z";
 943     Int_t secondBin = -1;
 944
 945     if (strcmp(fHistogramType, "NumberDensityPhi") == 0)
 946     {
 947       projAxis = "yz";
 948       secondBin = 4;
 949     }
 950
 951     for (UInt_t region = 0; region < fkRegions; region++)
 952     {
 953       if (!fTrackHist[region])
 954         continue;
 955
 956       TH1* leadingBiasTracks =  0;
 957       if (biasFromMC)
 958       {
 959         leadingBiasTracks = (TH1*) corrections->GetBias(kCFStepReconstructed, kCFStepTracked, region, projAxis, 0, -1, 1); // from MC
 960         Printf("WARNING: Using MC bias correction");
 961       }
 962       else
 963         leadingBiasTracks = (TH1*) GetBias(kCFStepBiasStudy, kCFStepReconstructed, region, projAxis, 0, -1, 1);          // from data
 964
 965       CorrectTracks(kCFStepReconstructed, kCFStepTracked, region, leadingBiasTracks, 2, secondBin);
 966       if (region == kMin && fCombineMinMax)
 967       {
 968         CorrectTracks(kCFStepReconstructed, kCFStepTracked, kMax, leadingBiasTracks, 2, secondBin);
 969         delete leadingBiasTracks;
 970         break;
 971       }
 972       delete leadingBiasTracks;
 973     }
 974
 975     TH1* leadingBiasEvents = 0;
 976     if (biasFromMC)
 977       leadingBiasEvents = (TH1*) corrections->GetBias(kCFStepReconstructed, kCFStepTracked, kToward, projAxis, 0, -1, 2); // from MC
 978     else
 979       leadingBiasEvents = (TH1*) GetBias(kCFStepBiasStudy, kCFStepReconstructed, kToward, projAxis, 0, -1, 2);          // from data
 980
 981     //new TCanvas; leadingBiasEvents->DrawCopy();
 982     CorrectEvents(kCFStepReconstructed, kCFStepTracked, leadingBiasEvents, 0);
 983     delete leadingBiasEvents;
 984
 985     // --- contamination correction ---
 986
 987     TH2D* contamination = corrections->GetTrackingContamination();
 988     if (corrections->fContaminationEnhancement)
 989     {
 990       Printf("Applying contamination enhancement");
 991
 992       for (Int_t x=1; x<=contamination->GetNbinsX(); x++)
 993         for (Int_t y=1; y<=contamination->GetNbinsX(); y++)
 994           contamination->SetBinContent(x, y, contamination->GetBinContent(x, y) * corrections->fContaminationEnhancement->GetBinContent(corrections->fContaminationEnhancement->GetXaxis()->FindBin(contamination->GetYaxis()->GetBinCenter(y))));
 995     }
 996     CorrectTracks(kCFStepTracked, kCFStepTrackedOnlyPrim, contamination, 0, 1);
 997     CorrectEvents(kCFStepTracked, kCFStepTrackedOnlyPrim, 0, 0);
 998     delete contamination;
 999
1000     // --- efficiency correction ---
1001     // correct single-particle efficiency for associated particles
1002     // in addition correct for efficiency on trigger particles (tracks AND events)
1003
1004     TH1* efficiencyCorrection = corrections->GetTrackingEfficiencyCorrection();
1005     // use kCFStepVertex as a temporary step (will be overwritten later anyway)
1006     CorrectTracks(kCFStepTrackedOnlyPrim, kCFStepVertex, efficiencyCorrection, 0, 1);
1007     delete efficiencyCorrection;
1008
1009     // correct pT,T
1010     efficiencyCorrection = corrections->GetTrackEfficiency(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, 1, -1, 2);
1011     CorrectEvents(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, efficiencyCorrection, 0);
1012     CorrectTracks(kCFStepVertex, kCFStepAnaTopology, efficiencyCorrection, 2);
1013     delete efficiencyCorrection;
1014
1015     // copy
1016     CorrectTracks(kCFStepAnaTopology, kCFStepVertex, 0, -1);
1017     CorrectEvents(kCFStepAnaTopology, kCFStepVertex, 0, 0);
1018
1019     // vertex correction on the level of events as function of multiplicity, weighting tracks and events with the same factor
1020     // practically independent of low pT cut
1021     TH1D* vertexCorrection = (TH1D*) corrections->GetEventEfficiency(kCFStepVertex, kCFStepTriggered, 1);
1022
1023     // convert stage from true multiplicity to observed multiplicity by simple conversion factor
1024     TH1D* vertexCorrectionObs = (TH1D*) vertexCorrection->Clone("vertexCorrection2");
1025     vertexCorrectionObs->Reset();
1026
1027     TF1* func = new TF1("func", "[1]+[0]/(x-[2])");
1028     // some defaults
1029     func->SetParameters(0.1, 1, -0.7);
1030     vertexCorrection->Fit(func, "0I", "", 0, 3);
1031     for (Int_t i=1; i<=vertexCorrectionObs->GetNbinsX(); i++)
1032     {
1033       Float_t xPos = 1.0 / 0.77 * vertexCorrectionObs->GetXaxis()->GetBinCenter(i);
1034       if (xPos < 3)
1035         vertexCorrectionObs->SetBinContent(i, func->Eval(xPos));
1036       else
1037         vertexCorrectionObs->SetBinContent(i, vertexCorrection->Interpolate(xPos));
1038     }
1039
1040     #if 0
1041     new TCanvas;
1042     vertexCorrection->DrawCopy();
1043     vertexCorrectionObs->SetLineColor(2);
1044     vertexCorrectionObs->DrawCopy("same");
1045     func->SetRange(0, 4);
1046     func->DrawClone("same");
1047     #endif
1048
1049     CorrectTracks(kCFStepVertex, kCFStepTriggered, vertexCorrectionObs, 3);
1050     CorrectEvents(kCFStepVertex, kCFStepTriggered, vertexCorrectionObs, 1);
1051     delete vertexCorrectionObs;
1052     delete vertexCorrection;
1053     delete func;
1054
1055     // copy
1056     CorrectTracks(kCFStepTriggered, kCFStepAll, 0, -1);
1057     CorrectEvents(kCFStepTriggered, kCFStepAll, 0, 0);
1058   }
1059   else if (strcmp(fHistogramType, "NumberDensityPhiCentrality") == 0)
1060   {
1061     // copy
1062     CorrectTracks(kCFStepReconstructed, kCFStepTracked, 0, -1);
1063     CorrectEvents(kCFStepReconstructed, kCFStepTracked, 0, -1);
1064
1065     // Dont use eta in the following, because it is a Delta-eta axis
1066
1067     // contamination correction
1068     // correct single-particle contamination for associated particles
1069
1070     TH1* contamination = corrections->GetTrackingContamination(1);
1071
1072     if (1)
1073     {
1074       Printf("Applying contamination enhancement");
1075
1076       for (Int_t bin = 1; bin <= contamination->GetNbinsX(); bin++)
1077       {
1078         printf("%f", contamination->GetBinContent(bin));
1079         if (contamination->GetBinContent(bin) > 0)
1080           contamination->SetBinContent(bin, 1.0 + 1.1 * (contamination->GetBinContent(bin) - 1.0));
1081         printf(" --> %f\n", contamination->GetBinContent(bin));
1082       }
1083     }
1084
1085     CorrectTracks(kCFStepTracked, kCFStepTrackedOnlyPrim, contamination, 1);
1086     delete contamination;
1087
1088     // correct for additional contamination due to trigger particle around phi ~ 0
1089     TH2* correlatedContamination = corrections->GetCorrelatedContamination();
1090     if (1)
1091     {
1092       Printf("Applying contamination enhancement");
1093
1094       for (Int_t bin = 1; bin <= correlatedContamination->GetNbinsX(); bin++)
1095         for (Int_t bin2 = 1; bin2 <= correlatedContamination->GetNbinsY(); bin2++)
1096         {
1097           printf("%f", correlatedContamination->GetBinContent(bin, bin2));
1098           if (correlatedContamination->GetBinContent(bin, bin2) > 0)
1099             correlatedContamination->SetBinContent(bin, bin2, 1.0 + 1.1 * (correlatedContamination->GetBinContent(bin, bin2) - 1.0));
1100           printf(" --> %f\n", correlatedContamination->GetBinContent(bin, bin2));
1101         }
1102     }
1103
1104     //new TCanvas; correlatedContamination->DrawCopy("COLZ");
1105     CorrectCorrelatedContamination(kCFStepTrackedOnlyPrim, 0, correlatedContamination);
1106
1107     delete correlatedContamination;
1108
1109     // TODO correct for contamination of trigger particles (for tracks AND events)
1110     CorrectEvents(kCFStepTracked, kCFStepTrackedOnlyPrim, 0, 0);
1111
1112     // --- efficiency correction ---
1113     // correct single-particle efficiency for associated particles
1114     // in addition correct for efficiency on trigger particles (tracks AND events)
1115
1116     // in bins of pT and centrality
1117     TH1* efficiencyCorrection = corrections->GetTrackingEfficiencyCorrectionCentrality();
1118     // use kCFStepAnaTopology as a temporary step
1119     CorrectTracks(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, efficiencyCorrection, 1, 3);
1120     delete efficiencyCorrection;
1121
1122     // correct pT,T in bins of pT and centrality
1123     efficiencyCorrection = corrections->GetTrackEfficiency(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, 1, 3, 2);
1124     CorrectEvents(kCFStepTrackedOnlyPrim, kCFStepVertex, efficiencyCorrection, 0, 1);
1125     CorrectTracks(kCFStepAnaTopology, kCFStepVertex, efficiencyCorrection, 2, 3);
1126     delete efficiencyCorrection;
1127
1128     // no correction for vertex finding efficiency and trigger efficiency needed in PbPb
1129     // copy
1130     CorrectTracks(kCFStepVertex, kCFStepAll, 0, -1);
1131     CorrectEvents(kCFStepVertex, kCFStepAll, 0, -1);
1132   }
1133   else
1134     AliFatal(Form("Unknown histogram for correction: %s", GetTitle()));
1135 }
1136
1137 //____________________________________________________________________
1138 TH1* AliUEHist::GetTrackEfficiency(CFStep step1, CFStep step2, Int_t axis1, Int_t axis2, Int_t source, Int_t axis3)
1139 {
1140   // creates a track-level efficiency by dividing step2 by step1
1141   // projected to axis1 and axis2 (optional if >= 0)
1142   //
1143   // source: 0 = fTrackHist; 1 = fTrackHistEfficiency; 2 = fTrackHistEfficiency rebinned for pT,T / pT,lead binning
1144
1145   // integrate over regions
1146   // cache it for efficiency (usually more than one efficiency is requested)
1147
1148   AliCFContainer* sourceContainer = 0;
1149
1150   if (source == 0)
1151   {
1152     if (!fCache)
1153     {
1154       fCache = (AliCFContainer*) fTrackHist[0]->Clone();
1155       for (UInt_t i = 1; i < fkRegions; i++)
1156         if (fTrackHist[i])
1157           fCache->Add(fTrackHist[i]);
1158     }
1159     sourceContainer = fCache;
1160   }
1161   else if (source == 1 || source == 2)
1162   {
1163     sourceContainer = fTrackHistEfficiency;
1164     // step offset because we start with kCFStepAnaTopology
1165     step1 = (CFStep) ((Int_t) step1 - (Int_t) kCFStepAnaTopology);
1166     step2 = (CFStep) ((Int_t) step2 - (Int_t) kCFStepAnaTopology);
1167   }
1168   else
1169     return 0;
1170
1171   // reset all limits and set the right ones except those in axis1, axis2 and axis3
1172   ResetBinLimits(sourceContainer->GetGrid(step1));
1173   ResetBinLimits(sourceContainer->GetGrid(step2));
1174   if (fEtaMax > fEtaMin && axis1 != 0 && axis2 != 0 && axis3 != 0)
1175   {
1176     sourceContainer->GetGrid(step1)->SetRangeUser(0, fEtaMin, fEtaMax);
1177     sourceContainer->GetGrid(step2)->SetRangeUser(0, fEtaMin, fEtaMax);
1178   }
1179   if (fPtMax > fPtMin && axis1 != 1 && axis2 != 1 && axis3 != 1)
1180   {
1181     sourceContainer->GetGrid(step1)->SetRangeUser(1, fPtMin, fPtMax);
1182     sourceContainer->GetGrid(step2)->SetRangeUser(1, fPtMin, fPtMax);
1183   }
1184   if (fCentralityMax > fCentralityMin && axis1 != 3 && axis2 != 3 && axis3 != 3)
1185   {
1186     sourceContainer->GetGrid(step1)->SetRangeUser(3, fCentralityMin, fCentralityMax);
1187     sourceContainer->GetGrid(step2)->SetRangeUser(3, fCentralityMin, fCentralityMax);
1188   }
1189
1190   TH1* measured = 0;
1191   TH1* generated = 0;
1192
1193   if (axis3 >= 0)
1194   {
1195     generated = sourceContainer->Project(step1, axis1, axis2, axis3);
1196     measured = sourceContainer->Project(step2, axis1, axis2, axis3);
1197   }
1198   else if (axis2 >= 0)
1199   {
1200     generated = sourceContainer->Project(step1, axis1, axis2);
1201     measured = sourceContainer->Project(step2, axis1, axis2);
1202   }
1203   else
1204   {
1205     generated = sourceContainer->Project(step1, axis1);
1206     measured = sourceContainer->Project(step2, axis1);
1207   }
1208
1209   // check for bins with less than 50 entries, print warning
1210   Int_t binBegin[3];
1211   Int_t binEnd[3];
1212
1213   binBegin[0] = 1;
1214   binBegin[1] = 1;
1215   binBegin[2] = 1;
1216
1217   binEnd[0] = generated->GetNbinsX();
1218   binEnd[1] = generated->GetNbinsY();
1219   binEnd[2] = generated->GetNbinsZ();
1220
1221   if (fEtaMax > fEtaMin)
1222   {
1223     if (axis1 == 0)
1224     {
1225       binBegin[0] = generated->GetXaxis()->FindBin(fEtaMin);
1226       binEnd[0]   = generated->GetXaxis()->FindBin(fEtaMax);
1227     }
1228     if (axis2 == 0)
1229     {
1230       binBegin[1] = generated->GetYaxis()->FindBin(fEtaMin);
1231       binEnd[1]   = generated->GetYaxis()->FindBin(fEtaMax);
1232     }
1233     if (axis3 == 0)
1234     {
1235       binBegin[2] = generated->GetZaxis()->FindBin(fEtaMin);
1236       binEnd[2]   = generated->GetZaxis()->FindBin(fEtaMax);
1237     }
1238   }
1239
1240   if (fPtMax > fPtMin)
1241   {
1242     // TODO this is just checking up to 15 for now
1243     Float_t ptMax = TMath::Min((Float_t) 15., fPtMax);
1244     if (axis1 == 1)
1245     {
1246       binBegin[0] = generated->GetXaxis()->FindBin(fPtMin);
1247       binEnd[0]   = generated->GetXaxis()->FindBin(ptMax);
1248     }
1249     if (axis2 == 1)
1250     {
1251       binBegin[1] = generated->GetYaxis()->FindBin(fPtMin);
1252       binEnd[1]   = generated->GetYaxis()->FindBin(ptMax);
1253     }
1254     if (axis3 == 1)
1255     {
1256       binBegin[2] = generated->GetZaxis()->FindBin(fPtMin);
1257       binEnd[2]   = generated->GetZaxis()->FindBin(ptMax);
1258     }
1259   }
1260
1261   Int_t total = 0;
1262   Int_t count = 0;
1263   Int_t vars[3];
1264
1265   for (Int_t i=0; i<3; i++)
1266     vars[i] = binBegin[i];
1267
1268   const Int_t limit = 50;
1269   while (1)
1270   {
1271     if (generated->GetDimension() == 1 && generated->GetBinContent(vars[0]) < limit)
1272     {
1273       Printf("Empty bin at %s=%.2f (%.2f entries)", generated->GetXaxis()->GetTitle(), generated->GetXaxis()->GetBinCenter(vars[0]), generated->GetBinContent(vars[0]));
1274       count++;
1275     }
1276     else if (generated->GetDimension() == 2 && generated->GetBinContent(vars[0], vars[1]) < limit)
1277     {
1278       Printf("Empty bin at %s=%.2f %s=%.2f (%.2f entries)",
1279         generated->GetXaxis()->GetTitle(), generated->GetXaxis()->GetBinCenter(vars[0]),
1280         generated->GetYaxis()->GetTitle(), generated->GetYaxis()->GetBinCenter(vars[1]),
1281         generated->GetBinContent(vars[0], vars[1]));
1282       count++;
1283     }
1284     else if (generated->GetDimension() == 3 && generated->GetBinContent(vars[0], vars[1], vars[2]) < limit)
1285     {
1286       Printf("Empty bin at %s=%.2f %s=%.2f %s=%.2f (%.2f entries)",
1287         generated->GetXaxis()->GetTitle(), generated->GetXaxis()->GetBinCenter(vars[0]),
1288         generated->GetYaxis()->GetTitle(), generated->GetYaxis()->GetBinCenter(vars[1]),
1289         generated->GetZaxis()->GetTitle(), generated->GetZaxis()->GetBinCenter(vars[2]),
1290         generated->GetBinContent(vars[0], vars[1], vars[2]));
1291       count++;
1292     }
1293
1294     vars[2]++;
1295     if (vars[2] == binEnd[2]+1)
1296     {
1297       vars[2] = binBegin[2];
1298       vars[1]++;
1299     }
1300
1301     if (vars[1] == binEnd[1]+1)
1302     {
1303       vars[1] = binBegin[1];
1304       vars[0]++;
1305     }
1306
1307     if (vars[0] == binEnd[0]+1)
1308       break;
1309     total++;
1310   }
1311
1312   Printf("Correction has %d empty bins (out of %d bins)", count, total);
1313
1314   // rebin if required
1315   if (source == 2)
1316   {
1317     TAxis* axis = fEventHist->GetGrid(0)->GetGrid()->GetAxis(0);
1318
1319     if (axis->GetNbins() < measured->GetNbinsX())
1320     {
1321       if (axis2 != -1)
1322       {
1323         // 2d rebin with variable axis does not exist in root
1324
1325         TH1* tmp = measured;
1326         measured = new TH2D(Form("%s_rebinned", tmp->GetName()), tmp->GetTitle(), axis->GetNbins(), axis->GetXbins()->GetArray(), tmp->GetNbinsY(), tmp->GetYaxis()->GetXbins()->GetArray());
1327         for (Int_t x=1; x<=tmp->GetNbinsX(); x++)
1328           for (Int_t y=1; y<=tmp->GetNbinsY(); y++)
1329           {
1330             ((TH2*) measured)->Fill(tmp->GetXaxis()->GetBinCenter(x), tmp->GetYaxis()->GetBinCenter(y), tmp->GetBinContent(x, y));
1331             measured->SetBinError(x, y, 0); // cannot trust bin error, set to 0
1332           }
1333         delete tmp;
1334
1335         tmp = generated;
1336         generated = new TH2D(Form("%s_rebinned", tmp->GetName()), tmp->GetTitle(), axis->GetNbins(), axis->GetXbins()->GetArray(), tmp->GetNbinsY(), tmp->GetYaxis()->GetXbins()->GetArray());
1337         for (Int_t x=1; x<=tmp->GetNbinsX(); x++)
1338           for (Int_t y=1; y<=tmp->GetNbinsY(); y++)
1339           {
1340             ((TH2*) generated)->Fill(tmp->GetXaxis()->GetBinCenter(x), tmp->GetYaxis()->GetBinCenter(y), tmp->GetBinContent(x, y));
1341             generated->SetBinError(x, y, 0); // cannot trust bin error, set to 0
1342           }
1343         delete tmp;
1344       }
1345       else
1346       {
1347         TH1* tmp = measured;
1348         measured = tmp->Rebin(axis->GetNbins(), Form("%s_rebinned", tmp->GetName()), axis->GetXbins()->GetArray());
1349         delete tmp;
1350
1351         tmp = generated;
1352         generated = tmp->Rebin(axis->GetNbins(), Form("%s_rebinned", tmp->GetName()), axis->GetXbins()->GetArray());
1353         delete tmp;
1354       }
1355     }
1356     else if (axis->GetNbins() > measured->GetNbinsX())
1357     {
1358       if (axis2 != -1)
1359         AliFatal("Rebinning only works for 1d at present");
1360
1361       // this is an unfortunate case where the number of bins has to be increased in principle
1362       // there is a region where the binning is finner in one histogram and a region where it is the other way round
1363       // this cannot be resolved in principle, but as we only calculate the ratio the bin in the second region get the same entries
1364       // only a certain binning is supported here
1365       if (axis->GetNbins() != 100 || measured->GetNbinsX() != 39)
1366         AliFatal(Form("Invalid binning --> %d %d", axis->GetNbins(), measured->GetNbinsX()));
1367
1368       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};
1369
1370       // reduce binning below 5 GeV/c
1371       TH1* tmp = measured;
1372       measured = tmp->Rebin(27, Form("%s_rebinned", tmp->GetName()), newBins);
1373       delete tmp;
1374
1375       // increase binning above 5 GeV/c
1376       tmp = measured;
1377       measured = new TH1F(Form("%s_rebinned2", tmp->GetName()), tmp->GetTitle(), axis->GetNbins(), axis->GetBinLowEdge(1), axis->GetBinUpEdge(axis->GetNbins()));
1378       for (Int_t bin=1; bin<=measured->GetNbinsX(); bin++)
1379       {
1380         measured->SetBinContent(bin, tmp->GetBinContent(tmp->FindBin(measured->GetBinCenter(bin))));
1381         measured->SetBinError(bin, tmp->GetBinError(tmp->FindBin(measured->GetBinCenter(bin))));
1382       }
1383       delete tmp;
1384
1385       // reduce binning below 5 GeV/c
1386       tmp = generated;
1387       generated = tmp->Rebin(27, Form("%s_rebinned", tmp->GetName()), newBins);
1388       delete tmp;
1389
1390       // increase binning above 5 GeV/c
1391       tmp = generated;
1392       generated = new TH1F(Form("%s_rebinned2", tmp->GetName()), tmp->GetTitle(), axis->GetNbins(), axis->GetBinLowEdge(1), axis->GetBinUpEdge(axis->GetNbins()));
1393       for (Int_t bin=1; bin<=generated->GetNbinsX(); bin++)
1394       {
1395         generated->SetBinContent(bin, tmp->GetBinContent(tmp->FindBin(generated->GetBinCenter(bin))));
1396         generated->SetBinError(bin, tmp->GetBinError(tmp->FindBin(generated->GetBinCenter(bin))));
1397       }
1398       delete tmp;
1399     }
1400   }
1401
1402   measured->Divide(measured, generated, 1, 1, "B");
1403
1404   delete generated;
1405
1406   ResetBinLimits(sourceContainer->GetGrid(step1));
1407   ResetBinLimits(sourceContainer->GetGrid(step2));
1408
1409   return measured;
1410 }
1411
1412 //____________________________________________________________________
1413 TH1* AliUEHist::GetEventEfficiency(CFStep step1, CFStep step2, Int_t axis1, Int_t axis2, Float_t ptLeadMin, Float_t ptLeadMax)
1414 {
1415   // creates a event-level efficiency by dividing step2 by step1
1416   // projected to axis1 and axis2 (optional if >= 0)
1417
1418   if (ptLeadMax > ptLeadMin)
1419   {
1420     fEventHist->GetGrid(step1)->SetRangeUser(0, ptLeadMin, ptLeadMax);
1421     fEventHist->GetGrid(step2)->SetRangeUser(0, ptLeadMin, ptLeadMax);
1422   }
1423
1424   TH1* measured = 0;
1425   TH1* generated = 0;
1426
1427   if (axis2 >= 0)
1428   {
1429     generated = fEventHist->Project(step1, axis1, axis2);
1430     measured = fEventHist->Project(step2, axis1, axis2);
1431   }
1432   else
1433   {
1434     generated = fEventHist->Project(step1, axis1);
1435     measured = fEventHist->Project(step2, axis1);
1436   }
1437
1438   measured->Divide(measured, generated, 1, 1, "B");
1439
1440   delete generated;
1441
1442   if (ptLeadMax > ptLeadMin)
1443   {
1444     fEventHist->GetGrid(step1)->SetRangeUser(0, 0, -1);
1445     fEventHist->GetGrid(step2)->SetRangeUser(0, 0, -1);
1446   }
1447
1448   return measured;
1449 }
1450
1451 //____________________________________________________________________
1452 void AliUEHist::WeightHistogram(TH3* hist1, TH1* hist2)
1453 {
1454   // weights each entry of the 3d histogram hist1 with the 1d histogram hist2
1455   // where the matching is done of the z axis of hist1 with the x axis of hist2
1456
1457   if (hist1->GetNbinsZ() != hist2->GetNbinsX())
1458     AliFatal(Form("Inconsistent binning %d %d", hist1->GetNbinsZ(), hist2->GetNbinsX()));
1459
1460   for (Int_t x=1; x<=hist1->GetNbinsX(); x++)
1461   {
1462     for (Int_t y=1; y<=hist1->GetNbinsY(); y++)
1463     {
1464       for (Int_t z=1; z<=hist1->GetNbinsZ(); z++)
1465       {
1466         if (hist2->GetBinContent(z) > 0)
1467         {
1468           hist1->SetBinContent(x, y, z, hist1->GetBinContent(x, y, z) / hist2->GetBinContent(z));
1469           hist1->SetBinError(x, y, z, hist1->GetBinError(x, y, z) / hist2->GetBinContent(z));
1470         }
1471         else
1472         {
1473           hist1->SetBinContent(x, y, z, 0);
1474           hist1->SetBinError(x, y, z, 0);
1475         }
1476       }
1477     }
1478   }
1479 }
1480
1481 //____________________________________________________________________
1482 TH1* AliUEHist::GetBias(CFStep step1, CFStep step2, Int_t region, const char* axis, Float_t leadPtMin, Float_t leadPtMax, Int_t weighting)
1483 {
1484   // extracts the track-level bias (integrating out the multiplicity) between two steps (dividing step2 by step1)
1485   // in the given region (sum over all regions is calculated if region == -1)
1486   // done by weighting the track-level distribution with the number of events as function of leading pT
1487   // and then calculating the ratio between the distributions
1488   // projected to axis which is a TH3::Project3D string, e.g. "x", or "yx"
1489   //   no projection is done if axis == 0
1490   // weighting: 0 = tracks weighted with events (as discussed above)
1491   //            1 = only track bias is returned
1492   //            2 = only event bias is returned
1493
1494   AliCFContainer* tmp = 0;
1495
1496   if (region == -1)
1497   {
1498     tmp = (AliCFContainer*) fTrackHist[0]->Clone();
1499     for (UInt_t i = 1; i < fkRegions; i++)
1500       if (fTrackHist[i])
1501         tmp->Add(fTrackHist[i]);
1502   }
1503   else if (region == kMin && fCombineMinMax)
1504   {
1505     tmp = (AliCFContainer*) fTrackHist[kMin]->Clone();
1506     tmp->Add(fTrackHist[kMax]);
1507   }
1508   else
1509     tmp = fTrackHist[region];
1510
1511   ResetBinLimits(tmp->GetGrid(step1));
1512   ResetBinLimits(fEventHist->GetGrid(step1));
1513   SetBinLimits(tmp->GetGrid(step1));
1514
1515   ResetBinLimits(tmp->GetGrid(step2));
1516   ResetBinLimits(fEventHist->GetGrid(step2));
1517   SetBinLimits(tmp->GetGrid(step2));
1518
1519   TH1D* events1 = (TH1D*)fEventHist->Project(step1, 0);
1520   TH3D* hist1 = (TH3D*)tmp->Project(step1, 0, tmp->GetNVar()-1, 2);
1521   if (weighting == 0)
1522     WeightHistogram(hist1, events1);
1523
1524   TH1D* events2 = (TH1D*)fEventHist->Project(step2, 0);
1525   TH3D* hist2 = (TH3D*)tmp->Project(step2, 0, tmp->GetNVar()-1, 2);
1526   if (weighting == 0)
1527     WeightHistogram(hist2, events2);
1528
1529   TH1* generated = hist1;
1530   TH1* measured = hist2;
1531
1532   if (weighting == 0 || weighting == 1)
1533   {
1534     if (axis)
1535     {
1536       if (leadPtMax > leadPtMin)
1537       {
1538         hist1->GetZaxis()->SetRangeUser(leadPtMin, leadPtMax);
1539         hist2->GetZaxis()->SetRangeUser(leadPtMin, leadPtMax);
1540       }
1541
1542       if (fEtaMax > fEtaMin && !TString(axis).Contains("x"))
1543       {
1544         hist1->GetXaxis()->SetRangeUser(fEtaMin, fEtaMax);
1545         hist2->GetXaxis()->SetRangeUser(fEtaMin, fEtaMax);
1546       }
1547
1548       generated = hist1->Project3D(axis);
1549       measured  = hist2->Project3D(axis);
1550
1551       // delete hists here if projection has been done
1552       delete hist1;
1553       delete hist2;
1554     }
1555     delete events1;
1556     delete events2;
1557   }
1558   else if (weighting == 2)
1559   {
1560     delete hist1;
1561     delete hist2;
1562     generated = events1;
1563     measured = events2;
1564   }
1565
1566   measured->Divide(generated);
1567
1568   delete generated;
1569
1570   ResetBinLimits(tmp->GetGrid(step1));
1571   ResetBinLimits(tmp->GetGrid(step2));
1572
1573   if ((region == -1) || (region == kMin && fCombineMinMax))
1574     delete tmp;
1575
1576   return measured;
1577 }
1578
1579 //____________________________________________________________________
1580 void AliUEHist::CorrectCorrelatedContamination(CFStep step, Int_t region, TH1* trackCorrection)
1581 {
1582   // corrects for the given factor in a small delta-eta and delta-phi window as function of pT,A and pT,T
1583
1584   if (!fTrackHist[region])
1585     return;
1586
1587   THnSparse* grid = fTrackHist[region]->GetGrid(step)->GetGrid();
1588
1589   Int_t var1 = 1;
1590   Int_t var2 = 2;
1591
1592   if (grid->GetAxis(var1)->GetNbins() != trackCorrection->GetNbinsX())
1593     AliFatal(Form("Invalid binning (var1): %d %d", grid->GetAxis(var1)->GetNbins(), trackCorrection->GetNbinsX()));
1594
1595   if (grid->GetAxis(var2)->GetNbins() != trackCorrection->GetNbinsY())
1596     AliFatal(Form("Invalid binning (var2): %d %d", grid->GetAxis(var2)->GetNbins(), trackCorrection->GetNbinsY()));
1597
1598   // optimized implementation
1599   for (Int_t binIdx = 0; binIdx < grid->GetNbins(); binIdx++)
1600   {
1601     Int_t bins[5];
1602
1603     Double_t value = grid->GetBinContent(binIdx, bins);
1604     Double_t error = grid->GetBinError(binIdx);
1605
1606     // check eta and phi axes
1607     if (TMath::Abs(grid->GetAxis(0)->GetBinCenter(bins[0])) > 0.1)
1608       continue;
1609     if (TMath::Abs(grid->GetAxis(4)->GetBinCenter(bins[4])) > 0.1)
1610       continue;
1611
1612     value *= trackCorrection->GetBinContent(bins[var1], bins[var2]);
1613     error *= trackCorrection->GetBinContent(bins[var1], bins[var2]);
1614
1615     grid->SetBinContent(bins, value);
1616     grid->SetBinError(bins, error);
1617   }
1618
1619   Printf("AliUEHist::CorrectCorrelatedContamination: Corrected.");
1620 }
1621
1622 //____________________________________________________________________
1623 TH2* AliUEHist::GetCorrelatedContamination()
1624 {
1625   // 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!)
1626
1627   Int_t step1 = kCFStepTrackedOnlyPrim;
1628   Int_t step2 = kCFStepTracked;
1629
1630   fTrackHist[0]->GetGrid(step1)->SetRangeUser(0, -0.01, 0.01); // delta eta
1631   fTrackHist[0]->GetGrid(step1)->SetRangeUser(4, -0.01, 0.01); // delta phi
1632   TH2* tracksStep1 = (TH2*) fTrackHist[0]->Project(step1, 1, 2);
1633
1634   fTrackHist[0]->GetGrid(step2)->SetRangeUser(0, -0.01, 0.01); // delta eta
1635   fTrackHist[0]->GetGrid(step2)->SetRangeUser(4, -0.01, 0.01); // delta phi
1636   TH2* tracksStep2 = (TH2*) fTrackHist[0]->Project(step2, 1, 2);
1637
1638   tracksStep1->Divide(tracksStep2);
1639
1640   TH1* triggersStep1 = fEventHist->Project(step1, 0);
1641   TH1* triggersStep2 = fEventHist->Project(step2, 0);
1642
1643   TH1* singleParticle = GetTrackingContamination(1);
1644
1645   for (Int_t x=1; x<=tracksStep1->GetNbinsX(); x++)
1646     for (Int_t y=1; y<=tracksStep1->GetNbinsY(); y++)
1647       if (singleParticle->GetBinContent(x) > 0)
1648         tracksStep1->SetBinContent(x, y, tracksStep1->GetBinContent(x, y) / triggersStep1->GetBinContent(y) * triggersStep2->GetBinContent(y) / singleParticle->GetBinContent(x));
1649       else
1650         tracksStep1->SetBinContent(x, y, 0);
1651
1652   delete singleParticle;
1653   delete tracksStep2;
1654   delete triggersStep1;
1655   delete triggersStep2;
1656
1657   return tracksStep1;
1658 }
1659
1660 //____________________________________________________________________
1661 TH2D* AliUEHist::GetTrackingEfficiency()
1662 {
1663   // extracts the tracking efficiency by calculating the efficiency from step kCFStepAnaTopology to kCFStepTrackedOnlyPrim
1664   // integrates over the regions and all other variables than pT and eta to increase the statistics
1665   //
1666   // returned histogram has to be deleted by the user
1667
1668   return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepAnaTopology, kCFStepTrackedOnlyPrim, 0, 1));
1669 }
1670
1671 //____________________________________________________________________
1672 TH2D* AliUEHist::GetTrackingEfficiencyCentrality()
1673 {
1674   // extracts the tracking efficiency by calculating the efficiency from step kCFStepAnaTopology to kCFStepTrackedOnlyPrim
1675   // integrates over the regions and all other variables than pT, centrality to increase the statistics
1676   //
1677   // returned histogram has to be deleted by the user
1678
1679   return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepAnaTopology, kCFStepTrackedOnlyPrim, 1, 3));
1680 }
1681
1682 //____________________________________________________________________
1683 TH1D* AliUEHist::GetTrackingEfficiency(Int_t axis)
1684 {
1685   // extracts the tracking efficiency by calculating the efficiency from step kCFStepAnaTopology to kCFStepTrackedOnlyPrim
1686   // integrates over the regions and all other variables than pT (axis == 0) and eta (axis == 1) to increase the statistics
1687
1688   return dynamic_cast<TH1D*> (GetTrackEfficiency(kCFStepAnaTopology, kCFStepTrackedOnlyPrim, axis));
1689 }
1690
1691 //____________________________________________________________________
1692 TH2D* AliUEHist::GetTrackingCorrection()
1693 {
1694   // extracts the tracking correction by calculating the efficiency from step kCFStepAnaTopology to kCFStepTracked
1695   // integrates over the regions and all other variables than pT and eta to increase the statistics
1696   //
1697   // returned histogram has to be deleted by the user
1698
1699   return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepTracked, kCFStepAnaTopology, 0, 1));
1700 }
1701
1702 //____________________________________________________________________
1703 TH1D* AliUEHist::GetTrackingCorrection(Int_t axis)
1704 {
1705   // extracts the tracking correction by calculating the efficiency from step kCFStepAnaTopology 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, kCFStepAnaTopology, axis));
1709 }
1710
1711 //____________________________________________________________________
1712 TH2D* AliUEHist::GetTrackingEfficiencyCorrection()
1713 {
1714   // extracts the tracking correction by calculating the efficiency from step kCFStepAnaTopology to kCFStepTracked
1715   // integrates over the regions and all other variables than pT and eta to increase the statistics
1716   //
1717   // returned histogram has to be deleted by the user
1718
1719   return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, 0, 1));
1720 }
1721
1722 //____________________________________________________________________
1723 TH2D* AliUEHist::GetTrackingEfficiencyCorrectionCentrality()
1724 {
1725   // extracts the tracking correction by calculating the efficiency from step kCFStepAnaTopology to kCFStepTracked
1726   // integrates over the regions and all other variables than pT and centrality to increase the statistics
1727   //
1728   // returned histogram has to be deleted by the user
1729
1730   return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, 1, 3));
1731 }
1732
1733 //____________________________________________________________________
1734 TH1D* AliUEHist::GetTrackingEfficiencyCorrection(Int_t axis)
1735 {
1736   // extracts the tracking correction by calculating the efficiency from step kCFStepAnaTopology to kCFStepTracked
1737   // integrates over the regions and all other variables than pT (axis == 0) and eta (axis == 1) to increase the statistics
1738
1739   return dynamic_cast<TH1D*> (GetTrackEfficiency(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, axis));
1740 }
1741
1742 //____________________________________________________________________
1743 TH2D* AliUEHist::GetTrackingContamination()
1744 {
1745   // extracts the tracking contamination by secondaries by calculating the efficiency from step kCFStepTrackedOnlyPrim to kCFStepTracked
1746   // integrates over the regions and all other variables than pT and eta to increase the statistics
1747   //
1748   // returned histogram has to be deleted by the user
1749
1750   return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepTracked, kCFStepTrackedOnlyPrim, 0, 1));
1751 }
1752
1753 //____________________________________________________________________
1754 TH2D* AliUEHist::GetTrackingContaminationCentrality()
1755 {
1756   // extracts the tracking contamination by secondaries by calculating the efficiency from step kCFStepTrackedOnlyPrim to kCFStepTracked
1757   // integrates over the regions and all other variables than pT and centrality to increase the statistics
1758   //
1759   // returned histogram has to be deleted by the user
1760
1761   return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepTracked, kCFStepTrackedOnlyPrim, 1, 3));
1762 }
1763
1764 //____________________________________________________________________
1765 TH1D* AliUEHist::GetTrackingContamination(Int_t axis)
1766 {
1767   // extracts the tracking contamination by secondaries by calculating the efficiency from step kCFStepTrackedOnlyPrim to kCFStepTracked
1768   // integrates over the regions and all other variables than pT (axis == 0) and eta (axis == 1) to increase the statistics
1769
1770   return dynamic_cast<TH1D*> (GetTrackEfficiency(kCFStepTracked, kCFStepTrackedOnlyPrim, axis));
1771 }
1772
1773 //____________________________________________________________________
1774 const char* AliUEHist::GetRegionTitle(Region region)
1775 {
1776   // returns the name of the given region
1777
1778   switch (region)
1779   {
1780     case kToward:
1781       return "Towards";
1782     case kAway:
1783       return "Away";
1784     case kMin:
1785       return (fCombineMinMax) ? "Transverse" : "Min";
1786     case kMax:
1787       return "Max";
1788   }
1789
1790   return 0;
1791 }
1792
1793 //____________________________________________________________________
1794 const char* AliUEHist::GetStepTitle(CFStep step)
1795 {
1796   // returns the name of the given step
1797
1798   switch (step)
1799   {
1800     case kCFStepAll:
1801       return "All events";
1802     case kCFStepTriggered:
1803       return "Triggered";
1804     case kCFStepVertex:
1805       return "Primary Vertex";
1806     case kCFStepAnaTopology:
1807       return "Required analysis topology";
1808     case kCFStepTrackedOnlyPrim:
1809       return "Tracked (matched MC, only primaries)";
1810     case kCFStepTracked:
1811       return "Tracked (matched MC, all)";
1812     case kCFStepReconstructed:
1813       return "Reconstructed";
1814     case kCFStepRealLeading:
1815       return "Correct leading particle identified";
1816     case kCFStepBiasStudy:
1817       return "Bias study applying tracking efficiency";
1818     case kCFStepBiasStudy2:
1819       return "Bias study applying tracking efficiency in two steps";
1820   }
1821
1822   return 0;
1823 }
1824
1825 //____________________________________________________________________
1826 void AliUEHist::CopyReconstructedData(AliUEHist* from)
1827 {
1828   // copies those histograms extracted from ESD to this object
1829
1830   // TODO at present only the pointers are copied
1831
1832   for (Int_t region=0; region<4; region++)
1833   {
1834     if (!fTrackHist[region])
1835       continue;
1836
1837     fTrackHist[region]->SetGrid(AliUEHist::kCFStepReconstructed, from->fTrackHist[region]->GetGrid(AliUEHist::kCFStepReconstructed));
1838     //fTrackHist[region]->SetGrid(AliUEHist::kCFStepTrackedOnlyPrim, from->fTrackHist[region]->GetGrid(AliUEHist::kCFStepTrackedOnlyPrim));
1839     fTrackHist[region]->SetGrid(AliUEHist::kCFStepBiasStudy,     from->fTrackHist[region]->GetGrid(AliUEHist::kCFStepBiasStudy));
1840   }
1841
1842   fEventHist->SetGrid(AliUEHist::kCFStepReconstructed, from->fEventHist->GetGrid(AliUEHist::kCFStepReconstructed));
1843   //fEventHist->SetGrid(AliUEHist::kCFStepTrackedOnlyPrim, from->fEventHist->GetGrid(AliUEHist::kCFStepTrackedOnlyPrim));
1844   fEventHist->SetGrid(AliUEHist::kCFStepBiasStudy,     from->fEventHist->GetGrid(AliUEHist::kCFStepBiasStudy));
1845 }
1846
1847 //____________________________________________________________________
1848 void AliUEHist::ExtendTrackingEfficiency(Bool_t verbose)
1849 {
1850   // fits the tracking efficiency at high pT with a constant and fills all bins with this tracking efficiency
1851
1852   Float_t fitRangeBegin = 5.01;
1853   Float_t fitRangeEnd = 14.99;
1854   Float_t extendRangeBegin = 10.01;
1855
1856   if (fTrackHistEfficiency->GetNVar() == 3)
1857   {
1858     TH1* obj = GetTrackingEfficiency(1);
1859
1860     if (verbose)
1861     {
1862       new TCanvas;
1863       obj->Draw();
1864     }
1865
1866     obj->Fit("pol0", (verbose) ? "+" : "0+", "SAME", fitRangeBegin, fitRangeEnd);
1867
1868     Float_t trackingEff = obj->GetFunction("pol0")->GetParameter(0);
1869
1870     obj = GetTrackingContamination(1);
1871
1872     if (verbose)
1873     {
1874       new TCanvas;
1875       obj->Draw();
1876     }
1877
1878     obj->Fit("pol0", (verbose) ? "+" : "0+", "SAME", fitRangeBegin, fitRangeEnd);
1879
1880     Float_t trackingCont = obj->GetFunction("pol0")->GetParameter(0);
1881
1882     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);
1883
1884     // extend for full pT range
1885     for (Int_t x = fTrackHistEfficiency->GetAxis(0, 0)->FindBin(fEtaMin); x <= fTrackHistEfficiency->GetAxis(0, 0)->FindBin(fEtaMax); x++)
1886       for (Int_t y = fTrackHistEfficiency->GetAxis(1, 0)->FindBin(extendRangeBegin); y <= fTrackHistEfficiency->GetNBins(1); y++)
1887         for (Int_t z = 1; z <= fTrackHistEfficiency->GetNBins(2); z++) // particle type axis
1888         {
1889
1890           Int_t bins[3];
1891           bins[0] = x;
1892           bins[1] = y;
1893           bins[2] = z;
1894
1895           fTrackHistEfficiency->GetGrid(0)->SetElement(bins, 100);
1896           fTrackHistEfficiency->GetGrid(1)->SetElement(bins, 100.0 * trackingEff);
1897           fTrackHistEfficiency->GetGrid(2)->SetElement(bins, 100.0 * trackingEff / trackingCont);
1898         }
1899   }
1900   else if (fTrackHistEfficiency->GetNVar() == 4)
1901   {
1902     // fit in centrality intervals of 20% for efficiency, one bin for contamination
1903     Float_t* trackingEff = 0;
1904     Float_t* trackingCont = 0;
1905     Int_t centralityBins = 5;
1906
1907     Printf("AliUEHist::ExtendTrackingEfficiency: Fitting efficiencies between %f and %f. Extending from %f onwards (within %f < eta < %f)", fitRangeBegin, fitRangeEnd, extendRangeBegin, fEtaMin, fEtaMax);
1908
1909     // 0 = eff; 1 = cont
1910     for (Int_t caseNo = 0; caseNo < 2; caseNo++)
1911     {
1912       Float_t* target = 0;
1913       Int_t centralityBinsLocal = centralityBins;
1914
1915       if (caseNo == 0)
1916       {
1917         trackingEff = new Float_t[centralityBinsLocal];
1918         target = trackingEff;
1919       }
1920       else
1921       {
1922         centralityBinsLocal = 1;
1923         trackingCont = new Float_t[centralityBinsLocal];
1924         target = trackingCont;
1925       }
1926
1927       for (Int_t i=0; i<centralityBinsLocal; i++)
1928       {
1929         SetCentralityRange(100.0/centralityBinsLocal*i + 0.1, 100.0/centralityBinsLocal*(i+1) - 0.1);
1930         TH1* proj = (caseNo == 0) ? GetTrackingEfficiency(1) : GetTrackingContamination(1);
1931         if (verbose)
1932         {
1933           new TCanvas;
1934           proj->DrawCopy();
1935         }
1936         if ((Int_t) proj->Fit("pol0", (verbose) ? "+" : "Q0+", "SAME", fitRangeBegin, fitRangeEnd) == 0)
1937           target[i] = proj->GetFunction("pol0")->GetParameter(0);
1938         else
1939           target[i] = 0;
1940         Printf("AliUEHist::ExtendTrackingEfficiency: case %d, centrality %d, eff %f", caseNo, i, target[i]);
1941       }
1942     }
1943
1944     // extend for full pT range
1945     for (Int_t x = fTrackHistEfficiency->GetAxis(0, 0)->FindBin(fEtaMin); x <= fTrackHistEfficiency->GetAxis(0, 0)->FindBin(fEtaMax); x++)
1946       for (Int_t y = fTrackHistEfficiency->GetAxis(1, 0)->FindBin(extendRangeBegin); y <= fTrackHistEfficiency->GetNBins(1); y++)
1947         for (Int_t z = 1; z <= fTrackHistEfficiency->GetNBins(2); z++) // particle type axis
1948         {
1949           for (Int_t z2 = 1; z2 <= fTrackHistEfficiency->GetNBins(3); z2++) // centrality axis
1950           {
1951
1952             Int_t bins[4];
1953             bins[0] = x;
1954             bins[1] = y;
1955             bins[2] = z;
1956             bins[3] = z2;
1957
1958             Int_t z2Bin = (Int_t) (fTrackHistEfficiency->GetAxis(3, 0)->GetBinCenter(z2) / (100.0 / centralityBins));
1959             //Printf("%d %d", z2, z2Bin);
1960
1961             fTrackHistEfficiency->GetGrid(0)->SetElement(bins, 100);
1962             fTrackHistEfficiency->GetGrid(1)->SetElement(bins, 100.0 * trackingEff[z2Bin]);
1963             if (trackingCont[0] > 0)
1964               fTrackHistEfficiency->GetGrid(2)->SetElement(bins, 100.0 * trackingEff[z2Bin] / trackingCont[0]);
1965             else
1966               fTrackHistEfficiency->GetGrid(2)->SetElement(bins, 0);
1967           }
1968        }
1969
1970      delete[] trackingEff;
1971      delete[] trackingCont;
1972    }
1973
1974    SetCentralityRange(0, 0);
1975 }
1976
1977 void AliUEHist::AdditionalDPhiCorrection(Int_t step)
1978 {
1979   // corrects the dphi distribution with an extra factor close to dphi ~ 0
1980
1981   Printf("WARNING: In AliUEHist::AdditionalDPhiCorrection.");
1982
1983   THnSparse* grid = fTrackHist[0]->GetGrid(step)->GetGrid();
1984
1985   // optimized implementation
1986   for (Int_t binIdx = 0; binIdx < grid->GetNbins(); binIdx++)
1987   {
1988     Int_t bins[5];
1989     Double_t value = grid->GetBinContent(binIdx, bins);
1990     Double_t error = grid->GetBinError(binIdx);
1991
1992     Float_t binCenter = grid->GetAxis(4)->GetBinCenter(bins[4]);
1993     if (TMath::Abs(binCenter) < 0.2)
1994     {
1995       value *= 0.985;
1996       error *= 0.985;
1997     }
1998     else if (TMath::Abs(binCenter) < 0.3)
1999     {
2000       value *= 0.9925;
2001       error *= 0.9925;
2002     }
2003
2004     grid->SetBinContent(bins, value);
2005     grid->SetBinError(bins, error);
2006   }
2007 }
2008
2009 void AliUEHist::Scale(Double_t factor)
2010 {
2011   // scales all contained histograms by the given factor
2012
2013   for (Int_t i=0; i<4; i++)
2014     if (fTrackHist[i])
2015       fTrackHist[i]->Scale(factor);
2016
2017   fEventHist->Scale(factor);
2018   fTrackHistEfficiency->Scale(factor);
2019 }
2020
2021 void AliUEHist::Reset()
2022 {
2023   // resets all contained histograms
2024
2025   for (Int_t i=0; i<4; i++)
2026     if (fTrackHist[i])
2027       for (Int_t step=0; step<fTrackHist[i]->GetNStep(); step++)
2028         fTrackHist[i]->GetGrid(step)->GetGrid()->Reset();
2029
2030   for (Int_t step=0; step<fEventHist->GetNStep(); step++)
2031     fEventHist->GetGrid(step)->GetGrid()->Reset();
2032
2033   for (Int_t step=0; step<fTrackHistEfficiency->GetNStep(); step++)
2034     fTrackHistEfficiency->GetGrid(step)->GetGrid()->Reset();
2035 }