PWG4/JetTasks/AliUEHist.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id: AliUEHist.cxx 20164 2007-08-14 15:31:50Z morsch $ */
  17
  18 //
  19 //
  20 // encapsulate histogram and corrections for one underlying event histogram
  21 //
  22 //
  23 // Author: Jan Fiete Grosse-Oetringhaus, Sara Vallero
  24
  25 #include "AliUEHist.h"
  26 #include "AliCFContainer.h"
  27 #include "THnSparse.h"
  28 #include "TMath.h"
  29 #include "TList.h"
  30 #include "TCollection.h"
  31 #include "TH1D.h"
  32 #include "TH2D.h"
  33 #include "TH3D.h"
  34 #include "AliLog.h"
  35 #include "TCanvas.h"
  36 #include "TF1.h"
  37
  38 ClassImp(AliUEHist)
  39
  40 const Int_t AliUEHist::fgkCFSteps = 10;
  41
  42 AliUEHist::AliUEHist(const char* reqHist) :
  43   TObject(),
  44   fkRegions(4),
  45   fEventHist(0),
  46   fTrackHistEfficiency(0),
  47   fEtaMin(0),
  48   fEtaMax(0),
  49   fPtMin(0),
  50   fPtMax(0),
  51   fContaminationEnhancement(0),
  52   fCombineMinMax(0),
  53   fCache(0)
  54 {
  55   // Constructor
  56
  57   for (Int_t i=0; i<fkRegions; i++)
  58     fTrackHist[i] = 0;
  59
  60   if (strlen(reqHist) == 0)
  61     return;
  62
  63   AliLog::SetClassDebugLevel("AliCFContainer", -1);
  64   AliLog::SetClassDebugLevel("AliCFGridSparse", -3);
  65
  66   const char* title = "";
  67
  68   // track level
  69   Int_t nTrackVars = 4; // eta vs pT vs pT,lead (vs delta phi) vs multiplicity
  70   Int_t iTrackBin[5];
  71   Double_t* trackBins[5];
  72   const char* trackAxisTitle[5];
  73
  74   // eta
  75   iTrackBin[0] = 20;
  76   Double_t etaBins[20+1];
  77   for (Int_t i=0; i<=iTrackBin[0]; i++)
  78     etaBins[i] = -1.0 + 0.1 * i;
  79   trackBins[0] = etaBins;
  80   trackAxisTitle[0] = "#eta";
  81
  82   // pT
  83   iTrackBin[1] = 39;
  84   Double_t pTBins[] = {0.0, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 12.0, 14.0, 16.0, 18.0, 20.0, 25.0, 30.0, 35.0, 40.0, 45.0, 50.0, 100.0};
  85   trackBins[1] = pTBins;
  86   trackAxisTitle[1] = "p_{T} (GeV/c)";
  87
  88   // pT,lead binning 1
  89   const Int_t kNLeadingpTBins = 100;
  90   Double_t leadingpTBins[kNLeadingpTBins+1];
  91   for (Int_t i=0; i<=kNLeadingpTBins; i++)
  92     leadingpTBins[i] = 0.5 * i;
  93
  94   // pT,lead binning 2
  95   const Int_t kNLeadingpTBins2 = 13;
  96   Double_t leadingpTBins2[] = { 0.0, 0.5, 1.0, 2.0, 4.0, 6.0, 8.0, 10.0, 15.0, 20.0, 30.0, 40.0, 50.0, 100.0 };
  97
  98   // phi,lead
  99   const Int_t kNLeadingPhiBins = 40;
 100   Double_t leadingPhiBins[kNLeadingPhiBins+1];
 101   for (Int_t i=0; i<=kNLeadingPhiBins; i++)
 102     leadingPhiBins[i] = -0.5 * TMath::Pi() + 1.0 / 40 * i * TMath::TwoPi();
 103
 104   // multiplicity
 105   const Int_t kNMultiplicityBins = 15;
 106   Double_t multiplicityBins[kNMultiplicityBins+1];
 107   for (Int_t i=0; i<=kNMultiplicityBins; i++)
 108     multiplicityBins[i] = -0.5 + i;
 109   multiplicityBins[kNMultiplicityBins] = 200;
 110
 111   // particle species
 112   const Int_t kNSpeciesBins = 4; // pi, K, p, rest
 113   Double_t speciesBins[] = { -0.5, 0.5, 1.5, 2.5, 3.5 };
 114
 115   trackBins[3] = multiplicityBins;
 116   iTrackBin[3] = kNMultiplicityBins;
 117   trackAxisTitle[3] = "multiplicity";
 118
 119   // selection depending on requested histogram
 120   Int_t axis = -1; // 0 = pT,lead, 1 = phi,lead
 121   if (strcmp(reqHist, "NumberDensitypT") == 0)
 122   {
 123     axis = 0;
 124     title = "d^{2}N_{ch}/d#phid#eta";
 125   }
 126   else if (strcmp(reqHist, "NumberDensityPhi") == 0)
 127   {
 128     axis = 1;
 129     title = "d^{2}N_{ch}/d#phid#eta";
 130   }
 131   else if (strcmp(reqHist, "SumpT") == 0)
 132   {
 133     axis = 0;
 134     title = "d^{2}#Sigma p_{T}/d#phid#eta";
 135   }
 136   else
 137     AliFatal(Form("Invalid histogram requested: %s", reqHist));
 138
 139   Int_t initRegions = fkRegions;
 140
 141   if (axis == 0)
 142   {
 143     trackBins[2] = leadingpTBins;
 144     iTrackBin[2] = kNLeadingpTBins;
 145     trackAxisTitle[2] = "leading p_{T} (GeV/c)";
 146
 147   }
 148   else if (axis == 1)
 149   {
 150     nTrackVars = 5;
 151     initRegions = 1;
 152
 153     iTrackBin[2] = kNLeadingpTBins2;
 154     trackBins[2] = leadingpTBins2;
 155     trackAxisTitle[2] = "leading p_{T} (GeV/c)";
 156
 157     iTrackBin[4] = kNLeadingPhiBins;
 158     trackBins[4] = leadingPhiBins;
 159     trackAxisTitle[4] = "#phi w.r.t leading track";
 160   }
 161
 162   for (Int_t i=0; i<initRegions; i++)
 163   {
 164     fTrackHist[i] = new AliCFContainer(Form("fTrackHist_%d", i), title, fgkCFSteps, nTrackVars, iTrackBin);
 165
 166     for (Int_t j=0; j<nTrackVars; j++)
 167     {
 168       fTrackHist[i]->SetBinLimits(j, trackBins[j]);
 169       fTrackHist[i]->SetVarTitle(j, trackAxisTitle[j]);
 170     }
 171
 172     SetStepNames(fTrackHist[i]);
 173   }
 174
 175   // track 3rd and 4th axis --> event 1st and 2nd axis
 176   fEventHist = new AliCFContainer("fEventHist", title, fgkCFSteps, 2, iTrackBin+2);
 177
 178   fEventHist->SetBinLimits(0, trackBins[2]);
 179   fEventHist->SetVarTitle(0, trackAxisTitle[2]);
 180
 181   fEventHist->SetBinLimits(1, trackBins[3]);
 182   fEventHist->SetVarTitle(1, trackAxisTitle[3]);
 183
 184   SetStepNames(fEventHist);
 185
 186   iTrackBin[2] = kNSpeciesBins;
 187
 188   fTrackHistEfficiency = new AliCFContainer("fTrackHistEfficiency", "Tracking efficiency", 3, 3, iTrackBin);
 189   fTrackHistEfficiency->SetBinLimits(0, trackBins[0]);
 190   fTrackHistEfficiency->SetVarTitle(0, trackAxisTitle[0]);
 191   fTrackHistEfficiency->SetBinLimits(1, trackBins[1]);
 192   fTrackHistEfficiency->SetVarTitle(1, trackAxisTitle[1]);
 193   fTrackHistEfficiency->SetBinLimits(2, speciesBins);
 194   fTrackHistEfficiency->SetVarTitle(2, "particle species");
 195 }
 196
 197 //_____________________________________________________________________________
 198 AliUEHist::AliUEHist(const AliUEHist &c) :
 199   TObject(),
 200   fkRegions(4),
 201   fEventHist(0),
 202   fTrackHistEfficiency(0),
 203   fEtaMin(0),
 204   fEtaMax(0),
 205   fPtMin(0),
 206   fPtMax(0),
 207   fContaminationEnhancement(0),
 208   fCombineMinMax(0),
 209   fCache(0)
 210 {
 211   //
 212   // AliUEHist copy constructor
 213   //
 214
 215   ((AliUEHist &) c).Copy(*this);
 216 }
 217
 218 //____________________________________________________________________
 219 void AliUEHist::SetStepNames(AliCFContainer* container)
 220 {
 221   // sets the names of the correction steps
 222
 223   for (Int_t i=0; i<fgkCFSteps; i++)
 224     container->SetStepTitle(i, GetStepTitle((CFStep) i));
 225 }
 226
 227 //____________________________________________________________________
 228 AliUEHist::~AliUEHist()
 229 {
 230   // Destructor
 231
 232   for (Int_t i=0; i<fkRegions; i++)
 233   {
 234     if (fTrackHist[i])
 235     {
 236       delete fTrackHist[i];
 237       fTrackHist[i] = 0;
 238     }
 239   }
 240
 241   if (fEventHist)
 242   {
 243     delete fEventHist;
 244     fEventHist = 0;
 245   }
 246
 247   if (fTrackHistEfficiency)
 248   {
 249     delete fTrackHistEfficiency;
 250     fTrackHistEfficiency = 0;
 251   }
 252
 253   if (fCache)
 254   {
 255     delete fCache;
 256     fCache = 0;
 257   }
 258 }
 259
 260 //____________________________________________________________________
 261 AliUEHist &AliUEHist::operator=(const AliUEHist &c)
 262 {
 263   // assigment operator
 264
 265   if (this != &c)
 266     ((AliUEHist &) c).Copy(*this);
 267
 268   return *this;
 269 }
 270
 271 //____________________________________________________________________
 272 void AliUEHist::Copy(TObject& c) const
 273 {
 274   // copy function
 275
 276   AliUEHist& target = (AliUEHist &) c;
 277
 278   for (Int_t i=0; i<fkRegions; i++)
 279     if (fTrackHist[i])
 280       target.fTrackHist[i] = dynamic_cast<AliCFContainer*> (fTrackHist[i]->Clone());
 281
 282   if (fEventHist)
 283     target.fEventHist = dynamic_cast<AliCFContainer*> (fEventHist->Clone());
 284
 285   if (fTrackHistEfficiency)
 286     target.fTrackHistEfficiency = dynamic_cast<AliCFContainer*> (fTrackHistEfficiency->Clone());
 287 }
 288
 289 //____________________________________________________________________
 290 Long64_t AliUEHist::Merge(TCollection* list)
 291 {
 292   // Merge a list of AliUEHist objects with this (needed for
 293   // PROOF).
 294   // Returns the number of merged objects (including this).
 295
 296   if (!list)
 297     return 0;
 298
 299   if (list->IsEmpty())
 300     return 1;
 301
 302   TIterator* iter = list->MakeIterator();
 303   TObject* obj;
 304
 305   // collections of objects
 306   const Int_t kMaxLists = fkRegions+2;
 307   TList** lists = new TList*[kMaxLists];
 308
 309   for (Int_t i=0; i<kMaxLists; i++)
 310     lists[i] = new TList;
 311
 312   Int_t count = 0;
 313   while ((obj = iter->Next())) {
 314
 315     AliUEHist* entry = dynamic_cast<AliUEHist*> (obj);
 316     if (entry == 0)
 317       continue;
 318
 319     for (Int_t i=0; i<fkRegions; i++)
 320       if (entry->fTrackHist[i])
 321         lists[i]->Add(entry->fTrackHist[i]);
 322
 323     lists[fkRegions]->Add(entry->fEventHist);
 324     lists[fkRegions+1]->Add(entry->fTrackHistEfficiency);
 325
 326     count++;
 327   }
 328   for (Int_t i=0; i<fkRegions; i++)
 329     if (fTrackHist[i])
 330       fTrackHist[i]->Merge(lists[i]);
 331
 332   fEventHist->Merge(lists[fkRegions]);
 333   fTrackHistEfficiency->Merge(lists[fkRegions+1]);
 334
 335   for (Int_t i=0; i<kMaxLists; i++)
 336     delete lists[i];
 337
 338   delete[] lists;
 339
 340   return count+1;
 341 }
 342
 343 //____________________________________________________________________
 344 void AliUEHist::SetBinLimits(AliCFGridSparse* grid)
 345 {
 346   // sets the bin limits in eta and pT defined by fEtaMin/Max, fPtMin/Max
 347
 348   if (fEtaMax > fEtaMin)
 349     grid->SetRangeUser(0, fEtaMin, fEtaMax);
 350   if (fPtMax > fPtMin)
 351     grid->SetRangeUser(1, fPtMin, fPtMax);
 352 }
 353
 354 //____________________________________________________________________
 355 void AliUEHist::ResetBinLimits(AliCFGridSparse* grid)
 356 {
 357   // resets all bin limits
 358
 359   for (Int_t i=0; i<grid->GetNVar(); i++)
 360     if (grid->GetGrid()->GetAxis(i)->TestBit(TAxis::kAxisRange))
 361       grid->SetRangeUser(i, 0, -1);
 362 }
 363
 364 //____________________________________________________________________
 365 void AliUEHist::CountEmptyBins(AliUEHist::CFStep step, Float_t ptLeadMin, Float_t ptLeadMax)
 366 {
 367   // prints the number of empty bins in the track end event histograms in the given step
 368
 369   Int_t binBegin[4];
 370   Int_t binEnd[4];
 371
 372   for (Int_t i=0; i<4; i++)
 373   {
 374     binBegin[i] = 1;
 375     binEnd[i]   = fTrackHist[0]->GetNBins(i);
 376   }
 377
 378   if (fEtaMax > fEtaMin)
 379   {
 380     binBegin[0] = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(0)->FindBin(fEtaMin);
 381     binEnd[0]   = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(0)->FindBin(fEtaMax);
 382   }
 383
 384   if (fPtMax > fPtMin)
 385   {
 386     binBegin[1] = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(1)->FindBin(fPtMin);
 387     binEnd[1]   = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(1)->FindBin(fPtMax);
 388   }
 389
 390   if (ptLeadMax > ptLeadMin)
 391   {
 392     binBegin[2] = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(2)->FindBin(ptLeadMin);
 393     binEnd[2]   = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(2)->FindBin(ptLeadMax);
 394   }
 395
 396   // start from multiplicity 1
 397   binBegin[3] = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(3)->FindBin(1);
 398
 399   for (Int_t region=0; region<fkRegions; region++)
 400   {
 401     Int_t total = 0;
 402     Int_t count = 0;
 403     Int_t vars[4];
 404
 405     for (Int_t i=0; i<4; i++)
 406       vars[i] = binBegin[i];
 407
 408     AliCFGridSparse* grid = fTrackHist[region]->GetGrid(step);
 409     while (1)
 410     {
 411       if (grid->GetElement(vars) == 0)
 412       {
 413         Printf("Empty bin at eta=%.2f pt=%.2f pt_lead=%.2f mult=%.1f",
 414           grid->GetBinCenter(0, vars[0]),
 415           grid->GetBinCenter(1, vars[1]),
 416           grid->GetBinCenter(2, vars[2]),
 417           grid->GetBinCenter(3, vars[3])
 418         );
 419         count++;
 420       }
 421
 422       vars[3]++;
 423       for (Int_t i=3; i>0; i--)
 424       {
 425         if (vars[i] == binEnd[i]+1)
 426         {
 427           vars[i] = binBegin[i];
 428           vars[i-1]++;
 429         }
 430       }
 431
 432       if (vars[0] == binEnd[0]+1)
 433         break;
 434       total++;
 435     }
 436
 437     Printf("Region %s has %d empty bins (out of %d bins)", GetRegionTitle((Region) region), count, total);
 438   }
 439 }
 440
 441 //____________________________________________________________________
 442 TH1D* AliUEHist::GetUEHist(AliUEHist::CFStep step, AliUEHist::Region region, Float_t ptLeadMin, Float_t ptLeadMax)
 443 {
 444   // Extracts the UE histogram at the given step and in the given region by projection and dividing tracks by events
 445   //
 446   // ptLeadMin, ptLeadMax: Only meaningful for vs. delta phi plot (third axis is ptLead)
 447   // Histogram has to be deleted by the caller of the function
 448
 449   // unzoom all axes
 450   ResetBinLimits(fTrackHist[region]->GetGrid(step));
 451   ResetBinLimits(fEventHist->GetGrid(step));
 452
 453   SetBinLimits(fTrackHist[region]->GetGrid(step));
 454
 455   TH1D* tracks = 0;
 456
 457   if (ptLeadMin < 0)
 458   {
 459     tracks = fTrackHist[region]->ShowProjection(2, step);
 460     tracks->GetYaxis()->SetTitle(fTrackHist[region]->GetTitle());
 461     if (fCombineMinMax && region == kMin)
 462     {
 463       ResetBinLimits(fTrackHist[kMax]->GetGrid(step));
 464       SetBinLimits(fTrackHist[kMax]->GetGrid(step));
 465
 466       TH1D* tracks2 = fTrackHist[kMax]->ShowProjection(2, step);
 467       tracks->Add(tracks2);
 468
 469       ResetBinLimits(fTrackHist[kMax]->GetGrid(step));
 470     }
 471
 472     // normalize to get a density (deta dphi)
 473     TAxis* axis = fTrackHist[region]->GetGrid(step)->GetAxis(0);
 474     Float_t phiRegion = TMath::TwoPi() / 3;
 475     if (!fCombineMinMax && region == kMin)
 476       phiRegion /= 2;
 477     Float_t normalization = phiRegion * (axis->GetBinUpEdge(axis->GetLast()) - axis->GetBinLowEdge(axis->GetFirst()));
 478     //Printf("Normalization: %f", normalization);
 479     tracks->Scale(1.0 / normalization);
 480
 481     TH1D* events = fEventHist->ShowProjection(0, step);
 482     tracks->Divide(events);
 483   }
 484   else
 485   {
 486     // 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
 487     // therefore the number density must be calculated here per leading pT bin and then summed
 488
 489     Int_t firstBin = fTrackHist[region]->GetAxis(2, step)->FindBin(ptLeadMin);
 490     Int_t lastBin = fTrackHist[region]->GetAxis(2, step)->FindBin(ptLeadMax);
 491
 492     for (Int_t bin=firstBin; bin<=lastBin; bin++)
 493     {
 494       fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(2)->SetRange(bin, bin);
 495       TH1D* tracksTmp = (TH1D*) fTrackHist[region]->GetGrid(step)->Project(4);
 496       Printf("Calculated histogram in bin %d --> %f tracks", bin, tracksTmp->Integral());
 497       fTrackHist[region]->GetGrid(step)->SetRangeUser(2, 0, -1);
 498
 499       // normalize to get a density (deta dphi)
 500       TAxis* axis = fTrackHist[region]->GetGrid(step)->GetAxis(0);
 501       Float_t normalization = fTrackHist[region]->GetGrid(step)->GetAxis(4)->GetBinWidth(1) * (axis->GetBinUpEdge(axis->GetLast()) - axis->GetBinLowEdge(axis->GetFirst()));
 502       //Printf("Normalization: %f", normalization);
 503       tracksTmp->Scale(1.0 / normalization);
 504
 505       TH1D* events = fEventHist->ShowProjection(0, step);
 506       Int_t nEvents = (Int_t) events->GetBinContent(bin);
 507       if (nEvents > 0)
 508         tracksTmp->Scale(1.0 / nEvents);
 509       Printf("Calculated histogram in bin %d --> %d events", bin, nEvents);
 510
 511       if (!tracks)
 512         tracks = tracksTmp;
 513       else
 514       {
 515         tracks->Add(tracksTmp);
 516         delete tracksTmp;
 517       }
 518     }
 519   }
 520
 521   ResetBinLimits(fTrackHist[region]->GetGrid(step));
 522
 523   return tracks;
 524 }
 525
 526 //____________________________________________________________________
 527 void AliUEHist::CorrectTracks(CFStep step1, CFStep step2, TH1* trackCorrection, Int_t var1, Int_t var2)
 528 {
 529   // corrects from step1 to step2 by multiplying the tracks with trackCorrection
 530   // trackCorrection can be a function of eta (var1 == 0), pT (var1 == 1), leading pT (var1 == 2), multiplicity (var1 == 3), delta phi (var1 == 4)
 531   // if var2 >= 0 a two dimension correction is assumed in trackCorrection
 532   //
 533   // if trackCorrection is 0, just copies content from step1 to step2
 534
 535   for (Int_t region=0; region<fkRegions; region++)
 536     CorrectTracks(step1, step2, region, trackCorrection, var1, var2);
 537 }
 538
 539 //____________________________________________________________________
 540 void AliUEHist::CorrectTracks(CFStep step1, CFStep step2, Int_t region, TH1* trackCorrection, Int_t var1, Int_t var2)
 541 {
 542   //
 543   // see documentation of CorrectTracks above
 544   //
 545
 546   if (!fTrackHist[region])
 547     return;
 548
 549   THnSparse* grid = fTrackHist[region]->GetGrid(step1)->GetGrid();
 550   THnSparse* target = fTrackHist[region]->GetGrid(step2)->GetGrid();
 551
 552   // clear target histogram
 553   target->Reset();
 554
 555   if (trackCorrection != 0)
 556   {
 557     if (grid->GetAxis(var1)->GetNbins() != trackCorrection->GetNbinsX())
 558       AliFatal(Form("Invalid binning (var1): %d %d", grid->GetAxis(var1)->GetNbins(), trackCorrection->GetNbinsX()));
 559
 560     if (var2 >= 0 && grid->GetAxis(var2)->GetNbins() != trackCorrection->GetNbinsY())
 561       AliFatal(Form("Invalid binning (var2): %d %d", grid->GetAxis(var2)->GetNbins(), trackCorrection->GetNbinsY()));
 562   }
 563
 564   // optimized implementation
 565   for (Int_t binIdx = 0; binIdx < grid->GetNbins(); binIdx++)
 566   {
 567     Int_t bins[5];
 568     Double_t value = grid->GetBinContent(binIdx, bins);
 569     Double_t error = grid->GetBinError(binIdx);
 570
 571     if (trackCorrection != 0)
 572     {
 573       if (var2 < 0)
 574       {
 575         value *= trackCorrection->GetBinContent(bins[var1]);
 576         error *= trackCorrection->GetBinContent(bins[var1]);
 577       }
 578       else
 579       {
 580         value *= trackCorrection->GetBinContent(bins[var1], bins[var2]);
 581         error *= trackCorrection->GetBinContent(bins[var1], bins[var2]);
 582       }
 583     }
 584
 585     target->SetBinContent(bins, value);
 586     target->SetBinError(bins, error);
 587   }
 588
 589   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);
 590 }
 591
 592 //____________________________________________________________________
 593 void AliUEHist::CorrectEvents(CFStep step1, CFStep step2, TH1D* eventCorrection, Int_t var)
 594 {
 595   // corrects from step1 to step2 by multiplying the events with eventCorrection
 596   // eventCorrection is as function of leading pT (var == 0) or multiplicity (var == 1)
 597   //
 598   // if eventCorrection is 0, just copies content from step1 to step2
 599
 600   AliCFGridSparse* grid = fEventHist->GetGrid(step1);
 601   AliCFGridSparse* target = fEventHist->GetGrid(step2);
 602
 603   // clear target histogram
 604   target->GetGrid()->Reset();
 605
 606   if (eventCorrection != 0 && grid->GetNBins(var) != eventCorrection->GetNbinsX())
 607     AliFatal(Form("Invalid binning: %d %d", grid->GetNBins(var), eventCorrection->GetNbinsX()));
 608
 609   Int_t bins[2];
 610   for (Int_t x = 1; x <= grid->GetNBins(0); x++)
 611   {
 612     bins[0] = x;
 613     for (Int_t y = 1; y <= grid->GetNBins(1); y++)
 614     {
 615       bins[1] = y;
 616
 617       Double_t value = grid->GetElement(bins);
 618       if (value != 0)
 619       {
 620         Double_t error = grid->GetElementError(bins);
 621
 622         if (eventCorrection != 0)
 623         {
 624           value *= eventCorrection->GetBinContent(bins[var]);
 625           error *= eventCorrection->GetBinContent(bins[var]);
 626         }
 627
 628         target->SetElement(bins, value);
 629         target->SetElementError(bins, error);
 630       }
 631     }
 632   }
 633
 634   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);
 635 }
 636
 637 //____________________________________________________________________
 638 void AliUEHist::Correct(AliUEHist* corrections)
 639 {
 640   // applies the given corrections to extract from the step kCFStepReconstructed all previous steps
 641   //
 642   // in this object the data is expected in the step kCFStepReconstructed
 643
 644   // ---- track level
 645
 646   // bias due to migration in leading pT (because the leading particle is not reconstructed, and the subleading is used)
 647   // extracted as function of leading pT
 648   for (Int_t region = 0; region < fkRegions; region++)
 649   {
 650     if (!fTrackHist[region])
 651       continue;
 652
 653     const char* projAxis = "z";
 654     Int_t secondBin = -1;
 655
 656     if (fTrackHist[region]->GetNVar() == 5)
 657     {
 658       projAxis = "yz";
 659       secondBin = 4;
 660     }
 661
 662     #if 0
 663       TH1* leadingBias = (TH1*) corrections->GetBias(kCFStepReconstructed, kCFStepTracked, region, projAxis); // from MC
 664       Printf("WARNING: Using MC bias correction");
 665     #else
 666       TH1* leadingBias = (TH1*) GetBias(kCFStepBiasStudy, kCFStepReconstructed, region, projAxis);          // from data
 667     #endif
 668     CorrectTracks(kCFStepReconstructed, kCFStepTracked, region, leadingBias, 2, secondBin);
 669     if (region == kMin && fCombineMinMax)
 670     {
 671       CorrectTracks(kCFStepReconstructed, kCFStepTracked, kMax, leadingBias, 2, secondBin);
 672       delete leadingBias;
 673       break;
 674     }
 675     delete leadingBias;
 676   }
 677   CorrectEvents(kCFStepReconstructed, kCFStepTracked, 0, 0);
 678
 679   // correct with kCFStepTracked --> kCFStepTrackedOnlyPrim
 680   TH2D* contamination = corrections->GetTrackingContamination();
 681   if (corrections->fContaminationEnhancement)
 682   {
 683     Printf("Applying contamination enhancement");
 684
 685     for (Int_t x=1; x<=contamination->GetNbinsX(); x++)
 686       for (Int_t y=1; y<=contamination->GetNbinsX(); y++)
 687         contamination->SetBinContent(x, y, contamination->GetBinContent(x, y) * corrections->fContaminationEnhancement->GetBinContent(corrections->fContaminationEnhancement->GetXaxis()->FindBin(contamination->GetYaxis()->GetBinCenter(y))));
 688   }
 689   CorrectTracks(kCFStepTracked, kCFStepTrackedOnlyPrim, contamination, 0, 1);
 690   CorrectEvents(kCFStepTracked, kCFStepTrackedOnlyPrim, 0, 0);
 691   delete contamination;
 692
 693   // correct with kCFStepTrackedOnlyPrim --> kCFStepAnaTopology
 694   TH2D* efficiencyCorrection = corrections->GetTrackingEfficiencyCorrection();
 695   CorrectTracks(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, efficiencyCorrection, 0, 1);
 696   CorrectEvents(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, 0, 0);
 697   delete efficiencyCorrection;
 698
 699   // copy
 700   CorrectTracks(kCFStepAnaTopology, kCFStepVertex, 0, -1);
 701   CorrectEvents(kCFStepAnaTopology, kCFStepVertex, 0, 0);
 702
 703   // vertex correction on the level of events as function of multiplicity, weighting tracks and events with the same factor
 704   // practically independent of low pT cut
 705   TH1D* vertexCorrection = (TH1D*) corrections->GetEventEfficiency(kCFStepVertex, kCFStepTriggered, 1);
 706
 707   // convert stage from true multiplicity to observed multiplicity by simple conversion factor
 708   TH1D* vertexCorrectionObs = (TH1D*) vertexCorrection->Clone("vertexCorrection2");
 709   vertexCorrectionObs->Reset();
 710
 711   TF1* func = new TF1("func", "[1]+[0]/(x-[2])");
 712   // some defaults
 713   func->SetParameters(0.1, 1, -0.7);
 714   vertexCorrection->Fit(func, "0I", "", 0, 3);
 715   for (Int_t i=1; i<=vertexCorrectionObs->GetNbinsX(); i++)
 716   {
 717     Float_t xPos = 1.0 / 0.77 * vertexCorrectionObs->GetXaxis()->GetBinCenter(i);
 718     if (xPos < 3)
 719       vertexCorrectionObs->SetBinContent(i, func->Eval(xPos));
 720     else
 721       vertexCorrectionObs->SetBinContent(i, vertexCorrection->Interpolate(xPos));
 722   }
 723
 724   #if 1
 725   new TCanvas;
 726   vertexCorrection->DrawCopy();
 727   vertexCorrectionObs->SetLineColor(2);
 728   vertexCorrectionObs->DrawCopy("same");
 729   func->SetRange(0, 4);
 730   func->DrawClone("same");
 731   #endif
 732
 733   CorrectTracks(kCFStepVertex, kCFStepTriggered, vertexCorrectionObs, 3);
 734   CorrectEvents(kCFStepVertex, kCFStepTriggered, vertexCorrectionObs, 1);
 735   delete vertexCorrectionObs;
 736   delete vertexCorrection;
 737   delete func;
 738
 739   // copy
 740   CorrectTracks(kCFStepTriggered, kCFStepAll, 0, -1);
 741   CorrectEvents(kCFStepTriggered, kCFStepAll, 0, 0);
 742 }
 743
 744 //____________________________________________________________________
 745 TH1* AliUEHist::GetTrackEfficiency(CFStep step1, CFStep step2, Int_t axis1, Int_t axis2, Int_t source)
 746 {
 747   // creates a track-level efficiency by dividing step2 by step1
 748   // projected to axis1 and axis2 (optional if >= 0)
 749   //
 750   // source: 0 = fTrackHist; 1 = fTrackHistEfficiency
 751
 752   // integrate over regions
 753   // cache it for efficiency (usually more than one efficiency is requested)
 754
 755   AliCFContainer* sourceContainer = 0;
 756
 757   if (source == 0)
 758   {
 759     if (!fCache)
 760     {
 761       fCache = (AliCFContainer*) fTrackHist[0]->Clone();
 762       for (Int_t i = 1; i < fkRegions; i++)
 763         if (fTrackHist[i])
 764           fCache->Add(fTrackHist[i]);
 765     }
 766     sourceContainer = fCache;
 767   }
 768   else if (source == 1)
 769   {
 770     sourceContainer = fTrackHistEfficiency;
 771     // step offset because we start with kCFStepAnaTopology
 772     step1 = (CFStep) ((Int_t) step1 - (Int_t) kCFStepAnaTopology);
 773     step2 = (CFStep) ((Int_t) step2 - (Int_t) kCFStepAnaTopology);
 774   }
 775   else
 776     return 0;
 777
 778   // reset all limits and set the right ones except those in axis1 and axis2
 779   ResetBinLimits(sourceContainer->GetGrid(step1));
 780   ResetBinLimits(sourceContainer->GetGrid(step2));
 781   if (fEtaMax > fEtaMin && axis1 != 0 && axis2 != 0)
 782   {
 783     sourceContainer->GetGrid(step1)->SetRangeUser(0, fEtaMin, fEtaMax);
 784     sourceContainer->GetGrid(step2)->SetRangeUser(0, fEtaMin, fEtaMax);
 785   }
 786   if (fPtMax > fPtMin && axis1 != 1 && axis2 != 1)
 787   {
 788     sourceContainer->GetGrid(step1)->SetRangeUser(1, fPtMin, fPtMax);
 789     sourceContainer->GetGrid(step2)->SetRangeUser(1, fPtMin, fPtMax);
 790   }
 791
 792   TH1* measured = 0;
 793   TH1* generated = 0;
 794
 795   if (axis2 >= 0)
 796   {
 797     generated = sourceContainer->Project(axis1, axis2, step1);
 798     measured = sourceContainer->Project(axis1, axis2, step2);
 799   }
 800   else
 801   {
 802     generated = sourceContainer->Project(axis1, step1);
 803     measured = sourceContainer->Project(axis1, step2);
 804   }
 805
 806   // check for bins with less than 100 entries, print warning
 807   Int_t binBegin[2];
 808   Int_t binEnd[2];
 809
 810   binBegin[0] = 1;
 811   binBegin[1] = 1;
 812
 813   binEnd[0] = generated->GetNbinsX();
 814   binEnd[1] = generated->GetNbinsY();
 815
 816   if (fEtaMax > fEtaMin)
 817   {
 818     if (axis1 == 0)
 819     {
 820       binBegin[0] = generated->GetXaxis()->FindBin(fEtaMin);
 821       binEnd[0]   = generated->GetXaxis()->FindBin(fEtaMax);
 822     }
 823     if (axis2 == 0)
 824     {
 825       binBegin[1] = generated->GetYaxis()->FindBin(fEtaMin);
 826       binEnd[1]   = generated->GetYaxis()->FindBin(fEtaMax);
 827     }
 828   }
 829
 830   if (fPtMax > fPtMin)
 831   {
 832     // TODO this is just checking up to 15 for now
 833     Float_t ptMax = TMath::Min((Float_t) 15., fPtMax);
 834     if (axis1 == 1)
 835     {
 836       binBegin[0] = generated->GetXaxis()->FindBin(fPtMin);
 837       binEnd[0]   = generated->GetXaxis()->FindBin(ptMax);
 838     }
 839     if (axis2 == 1)
 840     {
 841       binBegin[1] = generated->GetYaxis()->FindBin(fPtMin);
 842       binEnd[1]   = generated->GetYaxis()->FindBin(ptMax);
 843     }
 844   }
 845
 846   Int_t total = 0;
 847   Int_t count = 0;
 848   Int_t vars[2];
 849
 850   for (Int_t i=0; i<2; i++)
 851     vars[i] = binBegin[i];
 852
 853   const Int_t limit = 50;
 854   while (1)
 855   {
 856     if (generated->GetDimension() == 1 && generated->GetBinContent(vars[0]) < limit)
 857     {
 858       Printf("Empty bin at %s=%.2f (%.2f entries)", generated->GetXaxis()->GetTitle(), generated->GetXaxis()->GetBinCenter(vars[0]), generated->GetBinContent(vars[0]));
 859       count++;
 860     }
 861     else if (generated->GetDimension() == 2 && generated->GetBinContent(vars[0], vars[1]) < limit)
 862     {
 863       Printf("Empty bin at %s=%.2f %s=%.2f (%.2f entries)",
 864         generated->GetXaxis()->GetTitle(), generated->GetXaxis()->GetBinCenter(vars[0]),
 865         generated->GetYaxis()->GetTitle(), generated->GetYaxis()->GetBinCenter(vars[1]),
 866         generated->GetBinContent(vars[0], vars[1]));
 867       count++;
 868     }
 869
 870     vars[1]++;
 871     if (vars[1] == binEnd[1]+1)
 872     {
 873       vars[1] = binBegin[1];
 874       vars[0]++;
 875     }
 876
 877     if (vars[0] == binEnd[0]+1)
 878       break;
 879     total++;
 880   }
 881
 882   Printf("Correction has %d empty bins (out of %d bins)", count, total);
 883
 884   measured->Divide(measured, generated, 1, 1, "B");
 885
 886   delete generated;
 887
 888   ResetBinLimits(sourceContainer->GetGrid(step1));
 889   ResetBinLimits(sourceContainer->GetGrid(step2));
 890
 891   return measured;
 892 }
 893
 894 //____________________________________________________________________
 895 TH1* AliUEHist::GetEventEfficiency(CFStep step1, CFStep step2, Int_t axis1, Int_t axis2, Float_t ptLeadMin, Float_t ptLeadMax)
 896 {
 897   // creates a event-level efficiency by dividing step2 by step1
 898   // projected to axis1 and axis2 (optional if >= 0)
 899
 900   if (ptLeadMax > ptLeadMin)
 901   {
 902     fEventHist->GetGrid(step1)->SetRangeUser(0, ptLeadMin, ptLeadMax);
 903     fEventHist->GetGrid(step2)->SetRangeUser(0, ptLeadMin, ptLeadMax);
 904   }
 905
 906   TH1* measured = 0;
 907   TH1* generated = 0;
 908
 909   if (axis2 >= 0)
 910   {
 911     generated = fEventHist->Project(axis1, axis2, step1);
 912     measured = fEventHist->Project(axis1, axis2, step2);
 913   }
 914   else
 915   {
 916     generated = fEventHist->Project(axis1, step1);
 917     measured = fEventHist->Project(axis1, step2);
 918   }
 919
 920   measured->Divide(measured, generated, 1, 1, "B");
 921
 922   delete generated;
 923
 924   if (ptLeadMax > ptLeadMin)
 925   {
 926     fEventHist->GetGrid(step1)->SetRangeUser(0, 0, -1);
 927     fEventHist->GetGrid(step2)->SetRangeUser(0, 0, -1);
 928   }
 929
 930   return measured;
 931 }
 932
 933 //____________________________________________________________________
 934 void AliUEHist::WeightHistogram(TH3* hist1, TH1* hist2)
 935 {
 936   // weights each entry of the 3d histogram hist1 with the 1d histogram hist2
 937   // where the matching is done of the z axis of hist1 with the x axis of hist2
 938
 939   if (hist1->GetNbinsZ() != hist2->GetNbinsX())
 940     AliFatal(Form("Inconsistent binning %d %d", hist1->GetNbinsZ(), hist2->GetNbinsX()));
 941
 942   for (Int_t x=1; x<=hist1->GetNbinsX(); x++)
 943   {
 944     for (Int_t y=1; y<=hist1->GetNbinsY(); y++)
 945     {
 946       for (Int_t z=1; z<=hist1->GetNbinsZ(); z++)
 947       {
 948         if (hist2->GetBinContent(z) > 0)
 949         {
 950           hist1->SetBinContent(x, y, z, hist1->GetBinContent(x, y, z) / hist2->GetBinContent(z));
 951           hist1->SetBinError(x, y, z, hist1->GetBinError(x, y, z) / hist2->GetBinContent(z));
 952         }
 953         else
 954         {
 955           hist1->SetBinContent(x, y, z, 0);
 956           hist1->SetBinError(x, y, z, 0);
 957         }
 958       }
 959     }
 960   }
 961 }
 962
 963 //____________________________________________________________________
 964 TH1* AliUEHist::GetBias(CFStep step1, CFStep step2, Int_t region, const char* axis, Float_t leadPtMin, Float_t leadPtMax)
 965 {
 966   // extracts the track-level bias (integrating out the multiplicity) between two steps (dividing step2 by step1)
 967   // in the given region (sum over all regions is calculated if region == -1)
 968   // done by weighting the track-level distribution with the number of events as function of leading pT
 969   // and then calculating the ratio between the distributions
 970   // projected to axis which is a TH3::Project3D string, e.g. "x", or "yx"
 971   //   no projection is done if axis == 0
 972
 973   AliCFContainer* tmp = 0;
 974
 975   if (region == -1)
 976   {
 977     tmp = (AliCFContainer*) fTrackHist[0]->Clone();
 978     for (Int_t i = 1; i < fkRegions; i++)
 979       if (fTrackHist[i])
 980         tmp->Add(fTrackHist[i]);
 981   }
 982   else if (region == kMin && fCombineMinMax)
 983   {
 984     tmp = (AliCFContainer*) fTrackHist[kMin]->Clone();
 985     tmp->Add(fTrackHist[kMax]);
 986   }
 987   else
 988     tmp = fTrackHist[region];
 989
 990   ResetBinLimits(tmp->GetGrid(step1));
 991   ResetBinLimits(fEventHist->GetGrid(step1));
 992   SetBinLimits(tmp->GetGrid(step1));
 993
 994   ResetBinLimits(tmp->GetGrid(step2));
 995   ResetBinLimits(fEventHist->GetGrid(step2));
 996   SetBinLimits(tmp->GetGrid(step2));
 997
 998   TH1D* events1 = (TH1D*)fEventHist->Project(0, step1);
 999   TH3D* hist1 = (TH3D*)tmp->Project(0, tmp->GetNVar()-1, 2, step1);
1000   WeightHistogram(hist1, events1);
1001
1002   TH1D* events2 = (TH1D*)fEventHist->Project(0, step2);
1003   TH3D* hist2 = (TH3D*)tmp->Project(0, tmp->GetNVar()-1, 2, step2);
1004   WeightHistogram(hist2, events2);
1005
1006   TH1* generated = hist1;
1007   TH1* measured = hist2;
1008
1009   if (axis)
1010   {
1011     if (leadPtMax > leadPtMin)
1012     {
1013       hist1->GetZaxis()->SetRangeUser(leadPtMin, leadPtMax);
1014       hist2->GetZaxis()->SetRangeUser(leadPtMin, leadPtMax);
1015     }
1016
1017     if (fEtaMax > fEtaMin && !TString(axis).Contains("x"))
1018     {
1019       hist1->GetXaxis()->SetRangeUser(fEtaMin, fEtaMax);
1020       hist2->GetXaxis()->SetRangeUser(fEtaMin, fEtaMax);
1021     }
1022
1023     generated = hist1->Project3D(axis);
1024     measured  = hist2->Project3D(axis);
1025
1026     // delete hists here if projection has been done
1027     delete hist1;
1028     delete hist2;
1029   }
1030
1031   measured->Divide(generated);
1032
1033   delete events1;
1034   delete events2;
1035   delete generated;
1036
1037   ResetBinLimits(tmp->GetGrid(step1));
1038   ResetBinLimits(tmp->GetGrid(step2));
1039
1040   if ((region == -1) || (region == kMin && fCombineMinMax))
1041     delete tmp;
1042
1043   return measured;
1044 }
1045
1046 //____________________________________________________________________
1047 TH2D* AliUEHist::GetTrackingEfficiency()
1048 {
1049   // extracts the tracking efficiency by calculating the efficiency from step kCFStepAnaTopology to kCFStepTrackedOnlyPrim
1050   // integrates over the regions and all other variables than pT and eta to increase the statistics
1051   //
1052   // returned histogram has to be deleted by the user
1053
1054   return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepAnaTopology, kCFStepTrackedOnlyPrim, 0, 1));
1055 }
1056
1057 //____________________________________________________________________
1058 TH1D* AliUEHist::GetTrackingEfficiency(Int_t axis)
1059 {
1060   // extracts the tracking efficiency by calculating the efficiency from step kCFStepAnaTopology to kCFStepTrackedOnlyPrim
1061   // integrates over the regions and all other variables than pT (axis == 0) and eta (axis == 1) to increase the statistics
1062
1063   return dynamic_cast<TH1D*> (GetTrackEfficiency(kCFStepAnaTopology, kCFStepTrackedOnlyPrim, axis));
1064 }
1065
1066 //____________________________________________________________________
1067 TH2D* AliUEHist::GetTrackingCorrection()
1068 {
1069   // extracts the tracking correction by calculating the efficiency from step kCFStepAnaTopology to kCFStepTracked
1070   // integrates over the regions and all other variables than pT and eta to increase the statistics
1071   //
1072   // returned histogram has to be deleted by the user
1073
1074   return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepTracked, kCFStepAnaTopology, 0, 1));
1075 }
1076
1077 //____________________________________________________________________
1078 TH1D* AliUEHist::GetTrackingCorrection(Int_t axis)
1079 {
1080   // extracts the tracking correction by calculating the efficiency from step kCFStepAnaTopology to kCFStepTracked
1081   // integrates over the regions and all other variables than pT (axis == 0) and eta (axis == 1) to increase the statistics
1082
1083   return dynamic_cast<TH1D*> (GetTrackEfficiency(kCFStepTracked, kCFStepAnaTopology, axis));
1084 }
1085
1086 //____________________________________________________________________
1087 TH2D* AliUEHist::GetTrackingEfficiencyCorrection()
1088 {
1089   // extracts the tracking correction by calculating the efficiency from step kCFStepAnaTopology to kCFStepTracked
1090   // integrates over the regions and all other variables than pT and eta to increase the statistics
1091   //
1092   // returned histogram has to be deleted by the user
1093
1094   return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, 0, 1));
1095 }
1096
1097 //____________________________________________________________________
1098 TH1D* AliUEHist::GetTrackingEfficiencyCorrection(Int_t axis)
1099 {
1100   // extracts the tracking correction by calculating the efficiency from step kCFStepAnaTopology to kCFStepTracked
1101   // integrates over the regions and all other variables than pT (axis == 0) and eta (axis == 1) to increase the statistics
1102
1103   return dynamic_cast<TH1D*> (GetTrackEfficiency(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, axis));
1104 }
1105
1106 //____________________________________________________________________
1107 TH2D* AliUEHist::GetTrackingContamination()
1108 {
1109   // extracts the tracking contamination by secondaries by calculating the efficiency from step kCFStepTrackedOnlyPrim to kCFStepTracked
1110   // integrates over the regions and all other variables than pT and eta to increase the statistics
1111   //
1112   // returned histogram has to be deleted by the user
1113
1114   return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepTracked, kCFStepTrackedOnlyPrim, 0, 1));
1115 }
1116
1117 //____________________________________________________________________
1118 TH1D* AliUEHist::GetTrackingContamination(Int_t axis)
1119 {
1120   // extracts the tracking contamination by secondaries by calculating the efficiency from step kCFStepTrackedOnlyPrim to kCFStepTracked
1121   // integrates over the regions and all other variables than pT (axis == 0) and eta (axis == 1) to increase the statistics
1122
1123   return dynamic_cast<TH1D*> (GetTrackEfficiency(kCFStepTracked, kCFStepTrackedOnlyPrim, axis));
1124 }
1125
1126 //____________________________________________________________________
1127 const char* AliUEHist::GetRegionTitle(Region region)
1128 {
1129   // returns the name of the given region
1130
1131   switch (region)
1132   {
1133     case kToward:
1134       return "Towards";
1135     case kAway:
1136       return "Away";
1137     case kMin:
1138       return (fCombineMinMax) ? "Transverse" : "Min";
1139     case kMax:
1140       return "Max";
1141   }
1142
1143   return 0;
1144 }
1145
1146 //____________________________________________________________________
1147 const char* AliUEHist::GetStepTitle(CFStep step)
1148 {
1149   // returns the name of the given step
1150
1151   switch (step)
1152   {
1153     case kCFStepAll:
1154       return "All events";
1155     case kCFStepTriggered:
1156       return "Triggered";
1157     case kCFStepVertex:
1158       return "Primary Vertex";
1159     case kCFStepAnaTopology:
1160       return "Required analysis topology";
1161     case kCFStepTrackedOnlyPrim:
1162       return "Tracked (matched MC, only primaries)";
1163     case kCFStepTracked:
1164       return "Tracked (matched MC, all)";
1165     case kCFStepReconstructed:
1166       return "Reconstructed";
1167     case kCFStepRealLeading:
1168       return "Correct leading particle identified";
1169     case kCFStepBiasStudy:
1170       return "Bias study applying tracking efficiency";
1171     case kCFStepBiasStudy2:
1172       return "Bias study applying tracking efficiency in two steps";
1173   }
1174
1175   return 0;
1176 }
1177
1178 //____________________________________________________________________
1179 void AliUEHist::CopyReconstructedData(AliUEHist* from)
1180 {
1181   // copies those histograms extracted from ESD to this object
1182
1183   // TODO at present only the pointers are copied
1184
1185   for (Int_t region=0; region<4; region++)
1186   {
1187     if (!fTrackHist[region])
1188       continue;
1189
1190     fTrackHist[region]->SetGrid(AliUEHist::kCFStepReconstructed, from->fTrackHist[region]->GetGrid(AliUEHist::kCFStepReconstructed));
1191     //fTrackHist[region]->SetGrid(AliUEHist::kCFStepTrackedOnlyPrim, from->fTrackHist[region]->GetGrid(AliUEHist::kCFStepTrackedOnlyPrim));
1192     fTrackHist[region]->SetGrid(AliUEHist::kCFStepBiasStudy,     from->fTrackHist[region]->GetGrid(AliUEHist::kCFStepBiasStudy));
1193   }
1194
1195   fEventHist->SetGrid(AliUEHist::kCFStepReconstructed, from->fEventHist->GetGrid(AliUEHist::kCFStepReconstructed));
1196   //fEventHist->SetGrid(AliUEHist::kCFStepTrackedOnlyPrim, from->fEventHist->GetGrid(AliUEHist::kCFStepTrackedOnlyPrim));
1197   fEventHist->SetGrid(AliUEHist::kCFStepBiasStudy,     from->fEventHist->GetGrid(AliUEHist::kCFStepBiasStudy));
1198 }
1199
1200 //____________________________________________________________________
1201 void AliUEHist::ExtendTrackingEfficiency(Bool_t verbose)
1202 {
1203   // fits the tracking efficiency at high pT with a constant and fills all bins with this tracking efficiency
1204
1205   Float_t fitRangeBegin = 5.01;
1206   Float_t fitRangeEnd = 14.99;
1207   Float_t extendRangeBegin = 10.01;
1208
1209   TH1* obj = GetTrackingEfficiency(1);
1210
1211   if (verbose)
1212     new TCanvas; obj->Draw();
1213   obj->Fit("pol0", (verbose) ? "+" : "+0", "SAME", fitRangeBegin, fitRangeEnd);
1214
1215   Float_t trackingEff = obj->GetFunction("pol0")->GetParameter(0);
1216
1217   obj = GetTrackingContamination(1);
1218
1219   if (verbose)
1220     new TCanvas; obj->Draw();
1221   obj->Fit("pol0", (verbose) ? "+" : "+0", "SAME", fitRangeBegin, fitRangeEnd);
1222
1223   Float_t trackingCont = obj->GetFunction("pol0")->GetParameter(0);
1224
1225   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);
1226
1227   // extend up to pT 100
1228   for (Int_t x = fTrackHistEfficiency->GetAxis(0, 0)->FindBin(fEtaMin); x <= fTrackHistEfficiency->GetAxis(0, 0)->FindBin(fEtaMax); x++)
1229     for (Int_t y = fTrackHistEfficiency->GetAxis(1, 0)->FindBin(extendRangeBegin); y <= fTrackHistEfficiency->GetNBins(1); y++)
1230       for (Int_t z = 1; z <= fTrackHistEfficiency->GetNBins(2); z++) // particle type axis
1231       {
1232         Int_t bins[3];
1233         bins[0] = x;
1234         bins[1] = y;
1235         bins[2] = z;
1236
1237         fTrackHistEfficiency->GetGrid(0)->SetElement(bins, 100);
1238         fTrackHistEfficiency->GetGrid(1)->SetElement(bins, 100.0 * trackingEff);
1239         fTrackHistEfficiency->GetGrid(2)->SetElement(bins, 100.0 * trackingEff / trackingCont);
1240       }
1241 }
1242
1243 void AliUEHist::AdditionalDPhiCorrection(Int_t step)
1244 {
1245   // corrects the dphi distribution with an extra factor close to dphi ~ 0
1246
1247   Printf("WARNING: In AliUEHist::AdditionalDPhiCorrection.");
1248
1249   THnSparse* grid = fTrackHist[0]->GetGrid(step)->GetGrid();
1250
1251   // optimized implementation
1252   for (Int_t binIdx = 0; binIdx < grid->GetNbins(); binIdx++)
1253   {
1254     Int_t bins[5];
1255     Double_t value = grid->GetBinContent(binIdx, bins);
1256     Double_t error = grid->GetBinError(binIdx);
1257
1258     Float_t binCenter = grid->GetAxis(4)->GetBinCenter(bins[4]);
1259     if (TMath::Abs(binCenter) < 0.2)
1260     {
1261       value *= 0.985;
1262       error *= 0.985;
1263     }
1264     else if (TMath::Abs(binCenter) < 0.3)
1265     {
1266       value *= 0.9925;
1267       error *= 0.9925;
1268     }
1269
1270     grid->SetBinContent(bins, value);
1271     grid->SetBinError(bins, error);
1272   }
1273 }