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