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