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
  37 ClassImp(AliUEHist)
  38
  39 const Int_t AliUEHist::fgkCFSteps = 10;
  40
  41 AliUEHist::AliUEHist(const char* reqHist) :
  42   TObject(),
  43   fkRegions(4),
  44   fEventHist(0),
  45   fEtaMin(0),
  46   fEtaMax(0),
  47   fPtMin(0),
  48   fPtMax(0),
  49   fCombineMinMax(0),
  50   fCache(0)
  51 {
  52   // Constructor
  53
  54   for (Int_t i=0; i<fkRegions; i++)
  55     fTrackHist[i] = 0;
  56
  57   if (strlen(reqHist) == 0)
  58     return;
  59
  60   const char* title = "";
  61
  62   // track level
  63   Int_t nTrackVars = 4; // eta vs pT vs pT,lead (vs delta phi) vs multiplicity
  64   Int_t iTrackBin[5];
  65   Double_t* trackBins[5];
  66   const char* trackAxisTitle[5];
  67
  68   // eta
  69   iTrackBin[0] = 20;
  70   Double_t etaBins[20+1];
  71   for (Int_t i=0; i<=iTrackBin[0]; i++)
  72     etaBins[i] = -1.0 + 0.1 * i;
  73   trackBins[0] = etaBins;
  74   trackAxisTitle[0] = "#eta";
  75
  76   // pT
  77   iTrackBin[1] = 39;
  78   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};
  79   trackBins[1] = pTBins;
  80   trackAxisTitle[1] = "p_{T} (GeV/c)";
  81
  82   // pT,lead binning 1
  83   const Int_t kNLeadingpTBins = 100;
  84   Double_t leadingpTBins[kNLeadingpTBins+1];
  85   for (Int_t i=0; i<=kNLeadingpTBins; i++)
  86     leadingpTBins[i] = 0.5 * i;
  87
  88   // pT,lead binning 2
  89   const Int_t kNLeadingpTBins2 = 10;
  90   Double_t leadingpTBins2[kNLeadingpTBins2+1];
  91   for (Int_t i=0; i<=kNLeadingpTBins2; i++)
  92     leadingpTBins2[i] = 5.0 * i;
  93
  94   // phi,lead
  95   const Int_t kNLeadingPhiBins = 40;
  96   Double_t leadingPhiBins[kNLeadingPhiBins+1];
  97   for (Int_t i=0; i<=kNLeadingPhiBins; i++)
  98     leadingPhiBins[i] = -1.5 * TMath::Pi() + 1.0 / 40 * i * TMath::TwoPi();
  99
 100   // multiplicity
 101   const Int_t kNMultiplicityBins = 15;
 102   Double_t multiplicityBins[kNMultiplicityBins+1];
 103   for (Int_t i=0; i<=kNMultiplicityBins; i++)
 104     multiplicityBins[i] = -0.5 + i;
 105   multiplicityBins[kNMultiplicityBins] = 200;
 106
 107   trackBins[3] = multiplicityBins;
 108   iTrackBin[3] = kNMultiplicityBins;
 109   trackAxisTitle[3] = "multiplicity";
 110
 111   // selection depending on requested histogram
 112   Int_t axis = -1; // 0 = pT,lead, 1 = phi,lead
 113   if (strcmp(reqHist, "NumberDensitypT") == 0)
 114   {
 115     axis = 0;
 116     title = "d^{2}N_{ch}/d#phid#eta";
 117   }
 118   else if (strcmp(reqHist, "NumberDensityPhi") == 0)
 119   {
 120     axis = 1;
 121     title = "d^{2}N_{ch}/d#phid#eta";
 122   }
 123   else if (strcmp(reqHist, "SumpT") == 0)
 124   {
 125     axis = 0;
 126     title = "d^{2}#Sigma p_{T}/d#phid#eta";
 127   }
 128   else
 129     AliFatal(Form("Invalid histogram requested: %s", reqHist));
 130
 131   Int_t initRegions = fkRegions;
 132
 133   if (axis == 0)
 134   {
 135     trackBins[2] = leadingpTBins;
 136     iTrackBin[2] = kNLeadingpTBins;
 137     trackAxisTitle[2] = "leading p_{T} (GeV/c)";
 138
 139   }
 140   else if (axis == 1)
 141   {
 142     nTrackVars = 5;
 143     initRegions = 1;
 144
 145     iTrackBin[2] = kNLeadingpTBins2;
 146     trackBins[2] = leadingpTBins2;
 147     trackAxisTitle[2] = "leading p_{T} (GeV/c)";
 148
 149     iTrackBin[4] = kNLeadingPhiBins;
 150     trackBins[4] = leadingPhiBins;
 151     trackAxisTitle[4] = "#phi w.r.t leading track";
 152   }
 153
 154   for (Int_t i=0; i<initRegions; i++)
 155   {
 156     fTrackHist[i] = new AliCFContainer(Form("fTrackHist_%d", i), title, fgkCFSteps, nTrackVars, iTrackBin);
 157
 158     for (Int_t j=0; j<nTrackVars; j++)
 159     {
 160       fTrackHist[i]->SetBinLimits(j, trackBins[j]);
 161       fTrackHist[i]->SetVarTitle(j, trackAxisTitle[j]);
 162     }
 163
 164     SetStepNames(fTrackHist[i]);
 165   }
 166
 167   // track 3rd and 4th axis --> event 1st and 2nd axis
 168   fEventHist = new AliCFContainer("fEventHist", title, fgkCFSteps, 2, iTrackBin+2);
 169
 170   fEventHist->SetBinLimits(0, trackBins[2]);
 171   fEventHist->SetVarTitle(0, trackAxisTitle[2]);
 172
 173   fEventHist->SetBinLimits(1, trackBins[3]);
 174   fEventHist->SetVarTitle(1, trackAxisTitle[3]);
 175
 176   SetStepNames(fEventHist);
 177 }
 178
 179 //_____________________________________________________________________________
 180 AliUEHist::AliUEHist(const AliUEHist &c) :
 181   TObject(),
 182   fkRegions(4),
 183   fEventHist(0),
 184   fEtaMin(0),
 185   fEtaMax(0),
 186   fPtMin(0),
 187   fPtMax(0),
 188   fCombineMinMax(0),
 189   fCache(0)
 190 {
 191   //
 192   // AliUEHist copy constructor
 193   //
 194
 195   ((AliUEHist &) c).Copy(*this);
 196 }
 197
 198 //____________________________________________________________________
 199 void AliUEHist::SetStepNames(AliCFContainer* container)
 200 {
 201   // sets the names of the correction steps
 202
 203   for (Int_t i=0; i<fgkCFSteps; i++)
 204     container->SetStepTitle(i, GetStepTitle((CFStep) i));
 205 }
 206
 207 //____________________________________________________________________
 208 AliUEHist::~AliUEHist()
 209 {
 210   // Destructor
 211
 212   for (Int_t i=0; i<fkRegions; i++)
 213   {
 214     if (fTrackHist[i])
 215     {
 216       delete fTrackHist[i];
 217       fTrackHist[i] = 0;
 218     }
 219   }
 220
 221   if (fEventHist)
 222   {
 223     delete fEventHist;
 224     fEventHist = 0;
 225   }
 226
 227   if (fCache)
 228   {
 229     delete fCache;
 230     fCache = 0;
 231   }
 232 }
 233
 234 //____________________________________________________________________
 235 AliUEHist &AliUEHist::operator=(const AliUEHist &c)
 236 {
 237   // assigment operator
 238
 239   if (this != &c)
 240     ((AliUEHist &) c).Copy(*this);
 241
 242   return *this;
 243 }
 244
 245 //____________________________________________________________________
 246 void AliUEHist::Copy(TObject& c) const
 247 {
 248   // copy function
 249
 250   AliUEHist& target = (AliUEHist &) c;
 251
 252   for (Int_t i=0; i<fkRegions; i++)
 253     if (fTrackHist[i])
 254       target.fTrackHist[i] = dynamic_cast<AliCFContainer*> (fTrackHist[i]->Clone());
 255
 256   if (fEventHist)
 257     target.fEventHist = dynamic_cast<AliCFContainer*> (fEventHist->Clone());
 258 }
 259
 260 //____________________________________________________________________
 261 Long64_t AliUEHist::Merge(TCollection* list)
 262 {
 263   // Merge a list of AliUEHist objects with this (needed for
 264   // PROOF).
 265   // Returns the number of merged objects (including this).
 266
 267   if (!list)
 268     return 0;
 269
 270   if (list->IsEmpty())
 271     return 1;
 272
 273   TIterator* iter = list->MakeIterator();
 274   TObject* obj;
 275
 276   // collections of objects
 277   const Int_t kMaxLists = fkRegions+1;
 278   TList** lists = new TList*[kMaxLists];
 279
 280   for (Int_t i=0; i<kMaxLists; i++)
 281     lists[i] = new TList;
 282
 283   Int_t count = 0;
 284   while ((obj = iter->Next())) {
 285
 286     AliUEHist* entry = dynamic_cast<AliUEHist*> (obj);
 287     if (entry == 0)
 288       continue;
 289
 290     for (Int_t i=0; i<fkRegions; i++)
 291       if (entry->fTrackHist[i])
 292         lists[i]->Add(entry->fTrackHist[i]);
 293
 294     lists[fkRegions]->Add(entry->fEventHist);
 295
 296     count++;
 297   }
 298   for (Int_t i=0; i<fkRegions; i++)
 299     if (fTrackHist[i])
 300       fTrackHist[i]->Merge(lists[i]);
 301
 302   fEventHist->Merge(lists[fkRegions]);
 303
 304   for (Int_t i=0; i<kMaxLists; i++)
 305     delete lists[i];
 306
 307   delete[] lists;
 308
 309   return count+1;
 310 }
 311
 312 //____________________________________________________________________
 313 void AliUEHist::SetBinLimits(AliCFGridSparse* grid)
 314 {
 315   // sets the bin limits in eta and pT defined by fEtaMin/Max, fPtMin/Max
 316
 317   if (fEtaMax > fEtaMin)
 318     grid->SetRangeUser(0, fEtaMin, fEtaMax);
 319   if (fPtMax > fPtMin)
 320     grid->SetRangeUser(1, fPtMin, fPtMax);
 321 }
 322
 323 //____________________________________________________________________
 324 void AliUEHist::ResetBinLimits(AliCFGridSparse* grid)
 325 {
 326   // resets all bin limits
 327
 328   for (Int_t i=0; i<grid->GetNVar(); i++)
 329     if (grid->GetGrid()->GetAxis(i)->TestBit(TH1::kIsZoomed))
 330       grid->SetRangeUser(i, 0, -1);
 331 }
 332
 333 //____________________________________________________________________
 334 TH1D* AliUEHist::GetUEHist(AliUEHist::CFStep step, AliUEHist::Region region, Float_t ptLeadMin, Float_t ptLeadMax)
 335 {
 336   // Extracts the UE histogram at the given step and in the given region by projection and dividing tracks by events
 337   //
 338   // ptLeadMin, ptLeadMax: Only meaningful for vs. delta phi plot (third axis is ptLead)
 339   // Histogram has to be deleted by the caller of the function
 340
 341   // unzoom all axes
 342   ResetBinLimits(fTrackHist[region]->GetGrid(step));
 343   ResetBinLimits(fEventHist->GetGrid(step));
 344
 345   SetBinLimits(fTrackHist[region]->GetGrid(step));
 346
 347   TH1D* tracks = 0;
 348
 349   if (ptLeadMin < 0)
 350   {
 351     tracks = fTrackHist[region]->ShowProjection(2, step);
 352     tracks->GetYaxis()->SetTitle(fTrackHist[region]->GetTitle());
 353     if (fCombineMinMax && region == kMin)
 354     {
 355       ResetBinLimits(fTrackHist[kMax]->GetGrid(step));
 356       SetBinLimits(fTrackHist[kMax]->GetGrid(step));
 357
 358       TH1D* tracks2 = fTrackHist[kMax]->ShowProjection(2, step);
 359       tracks->Add(tracks2);
 360
 361       ResetBinLimits(fTrackHist[kMax]->GetGrid(step));
 362     }
 363
 364     // normalize to get a density (deta dphi)
 365     TAxis* axis = fTrackHist[region]->GetGrid(step)->GetAxis(0);
 366     Float_t phiRegion = TMath::TwoPi() / 3;
 367     if (!fCombineMinMax && region == kMin)
 368       phiRegion /= 2;
 369     Float_t normalization = phiRegion * (axis->GetBinUpEdge(axis->GetLast()) - axis->GetBinLowEdge(axis->GetFirst()));
 370     //Printf("Normalization: %f", normalization);
 371     tracks->Scale(1.0 / normalization);
 372
 373     TH1D* events = fEventHist->ShowProjection(0, step);
 374     tracks->Divide(events);
 375   }
 376   else
 377   {
 378     fTrackHist[region]->GetGrid(step)->SetRangeUser(2, ptLeadMin, ptLeadMax);
 379     tracks = fTrackHist[region]->GetGrid(step)->Project(4);
 380     fTrackHist[region]->GetGrid(step)->SetRangeUser(2, 0, -1);
 381
 382     // normalize to get a density (deta dphi)
 383     TAxis* axis = fTrackHist[region]->GetGrid(step)->GetAxis(0);
 384     Float_t normalization = fTrackHist[region]->GetGrid(step)->GetAxis(4)->GetBinWidth(1) * (axis->GetBinUpEdge(axis->GetLast()) - axis->GetBinLowEdge(axis->GetFirst()));
 385     //Printf("Normalization: %f", normalization);
 386     tracks->Scale(1.0 / normalization);
 387
 388     TH1D* events = fEventHist->ShowProjection(0, step);
 389     Int_t nEvents = (Int_t) events->Integral(events->FindBin(ptLeadMin), events->FindBin(ptLeadMax));
 390     if (nEvents > 0)
 391       tracks->Scale(1.0 / nEvents);
 392   }
 393
 394   ResetBinLimits(fTrackHist[region]->GetGrid(step));
 395
 396   return tracks;
 397 }
 398
 399 //____________________________________________________________________
 400 void AliUEHist::CorrectTracks(CFStep step1, CFStep step2, TH1* trackCorrection, Int_t var1, Int_t var2)
 401 {
 402   // corrects from step1 to step2 by multiplying the tracks with trackCorrection
 403   // trackCorrection can be a function of eta (var1 == 0), pT (var1 == 1), leading pT (var1 == 2), multiplicity (var1 == 3), delta phi (var1 == 4)
 404   // if var2 >= 0 a two dimension correction is assumed in trackCorrection
 405   //
 406   // if trackCorrection is 0, just copies content from step1 to step2
 407
 408   for (Int_t region=0; region<fkRegions; region++)
 409   {
 410     if (!fTrackHist[region])
 411       continue;
 412
 413     THnSparse* grid = fTrackHist[region]->GetGrid(step1)->GetGrid();
 414     THnSparse* target = fTrackHist[region]->GetGrid(step2)->GetGrid();
 415
 416     if (trackCorrection != 0)
 417     {
 418       if (grid->GetAxis(var1)->GetNbins() != trackCorrection->GetNbinsX())
 419         AliFatal(Form("Invalid binning (var1): %d %d", grid->GetAxis(var1)->GetNbins(), trackCorrection->GetNbinsX()));
 420
 421       if (var2 >= 0 && grid->GetAxis(var2)->GetNbins() != trackCorrection->GetNbinsY())
 422         AliFatal(Form("Invalid binning (var2): %d %d", grid->GetAxis(var2)->GetNbins(), trackCorrection->GetNbinsY()));
 423     }
 424
 425     // optimized implementation
 426     for (Int_t binIdx = 0; binIdx < grid->GetNbins(); binIdx++)
 427     {
 428       Int_t bins[5];
 429       Double_t value = grid->GetBinContent(binIdx, bins);
 430       Double_t error = grid->GetBinError(binIdx);
 431
 432       if (trackCorrection != 0)
 433       {
 434         if (var2 < 0)
 435         {
 436           value *= trackCorrection->GetBinContent(bins[var1]);
 437           error *= trackCorrection->GetBinContent(bins[var1]);
 438         }
 439         else
 440         {
 441           value *= trackCorrection->GetBinContent(bins[var1], bins[var2]);
 442           error *= trackCorrection->GetBinContent(bins[var1], bins[var2]);
 443         }
 444       }
 445
 446       target->SetBinContent(bins, value);
 447       target->SetBinError(bins, error);
 448     }
 449   }
 450 }
 451
 452 //____________________________________________________________________
 453 void AliUEHist::CorrectEvents(CFStep step1, CFStep step2, TH1D* eventCorrection, Int_t var)
 454 {
 455   // corrects from step1 to step2 by multiplying the events with eventCorrection
 456   // eventCorrection is as function of leading pT (var == 0) or multiplicity (var == 1)
 457   //
 458   // if eventCorrection is 0, just copies content from step1 to step2
 459
 460   AliCFGridSparse* grid = fEventHist->GetGrid(step1);
 461   AliCFGridSparse* target = fEventHist->GetGrid(step2);
 462
 463   if (eventCorrection != 0 && grid->GetNBins(var) != eventCorrection->GetNbinsX())
 464     AliFatal(Form("Invalid binning: %d %d", grid->GetNBins(var), eventCorrection->GetNbinsX()));
 465
 466   Int_t bins[2];
 467   for (Int_t x = 1; x <= grid->GetNBins(0); x++)
 468   {
 469     bins[0] = x;
 470     for (Int_t y = 1; y <= grid->GetNBins(1); y++)
 471     {
 472       bins[1] = y;
 473
 474       Double_t value = grid->GetElement(bins);
 475       if (value != 0)
 476       {
 477         Double_t error = grid->GetElementError(bins);
 478
 479         if (eventCorrection != 0)
 480         {
 481           value *= eventCorrection->GetBinContent(bins[var]);
 482           error *= eventCorrection->GetBinContent(bins[var]);
 483         }
 484
 485         target->SetElement(bins, value);
 486         target->SetElementError(bins, error);
 487       }
 488     }
 489   }
 490 }
 491
 492 //____________________________________________________________________
 493 void AliUEHist::Correct(AliUEHist* corrections)
 494 {
 495   // applies the given corrections to extract from the step kCFStepReconstructed all previous steps
 496   //
 497   // in this object the data is expected in the step kCFStepReconstructed
 498
 499   // ---- track level
 500
 501   // bias due to migration in leading pT (because the leading particle is not reconstructed, and the subleading is used)
 502   // extracted as function of leading pT
 503   //TH1D* leadingBias = (TH1D*) corrections->GetBias(kCFStepReconstructed, kCFStepTracked, "z"); // from MC
 504   TH1D* leadingBias = (TH1D*) GetBias(kCFStepBiasStudy, kCFStepReconstructed, "z");              // from data
 505   CorrectTracks(kCFStepReconstructed, kCFStepTracked, leadingBias, 2);
 506   CorrectEvents(kCFStepReconstructed, kCFStepTracked, 0, 0);
 507   delete leadingBias;
 508
 509   // correct with kCFStepTracked --> kCFStepTrackedOnlyPrim
 510   TH2D* contamination = corrections->GetTrackingContamination();
 511   CorrectTracks(kCFStepTracked, kCFStepTrackedOnlyPrim, contamination, 0, 1);
 512   CorrectEvents(kCFStepTracked, kCFStepTrackedOnlyPrim, 0, 0);
 513   delete contamination;
 514
 515   // correct with kCFStepTracked --> kCFStepAnaTopology
 516   TH2D* trackingCorrection = corrections->GetTrackingCorrection();
 517   CorrectTracks(kCFStepTracked, kCFStepAnaTopology, trackingCorrection, 0, 1);
 518   CorrectEvents(kCFStepTracked, kCFStepAnaTopology, 0, 0);
 519   delete trackingCorrection;
 520
 521   // copy
 522   CorrectTracks(kCFStepAnaTopology, kCFStepVertex, 0, -1);
 523   CorrectEvents(kCFStepAnaTopology, kCFStepVertex, 0, 0);
 524
 525   // vertex correction on the level of events as function of multiplicity, weighting tracks and events with the same factor
 526   // practically independent of low pT cut
 527   TH1D* vertexCorrection = (TH1D*) corrections->GetEventEfficiency(kCFStepVertex, kCFStepTriggered, 1);
 528
 529   // convert stage from true multiplicity to observed multiplicity by simple conversion factor
 530   TH1D* vertexCorrectionObs = (TH1D*) vertexCorrection->Clone("vertexCorrection2");
 531   vertexCorrectionObs->Reset();
 532
 533   for (Int_t i=1; i<=vertexCorrectionObs->GetNbinsX(); i++)
 534     vertexCorrectionObs->SetBinContent(i, vertexCorrection->Interpolate(1.0 / 0.77 * vertexCorrectionObs->GetXaxis()->GetBinCenter(i)));
 535
 536 /*  new TCanvas;
 537   vertexCorrection->DrawCopy();
 538   vertexCorrectionObs->SetLineColor(2);
 539   vertexCorrectionObs->DrawCopy("same");*/
 540
 541   CorrectTracks(kCFStepVertex, kCFStepTriggered, vertexCorrectionObs, 3);
 542   CorrectEvents(kCFStepVertex, kCFStepTriggered, vertexCorrectionObs, 1);
 543   delete vertexCorrectionObs;
 544   delete vertexCorrection;
 545
 546   // copy
 547   CorrectTracks(kCFStepTriggered, kCFStepAll, 0, -1);
 548   CorrectEvents(kCFStepTriggered, kCFStepAll, 0, 0);
 549 }
 550
 551 //____________________________________________________________________
 552 TH1* AliUEHist::GetTrackEfficiency(CFStep step1, CFStep step2, Int_t axis1, Int_t axis2)
 553 {
 554   // creates a track-level efficiency by dividing step2 by step1
 555   // projected to axis1 and axis2 (optional if >= 0)
 556
 557   // integrate over regions
 558   // cache it for efficiency (usually more than one efficiency is requested)
 559   if (!fCache)
 560   {
 561     fCache = (AliCFContainer*) fTrackHist[0]->Clone();
 562     for (Int_t i = 1; i < fkRegions; i++)
 563       if (fTrackHist[i])
 564         fCache->Add(fTrackHist[i]);
 565   }
 566
 567   // reset all limits and set the right ones except those in axis1 and axis2
 568   ResetBinLimits(fCache->GetGrid(step1));
 569   ResetBinLimits(fCache->GetGrid(step2));
 570   if (fEtaMax > fEtaMin && axis1 != 0 && axis2 != 0)
 571   {
 572     fCache->GetGrid(step1)->SetRangeUser(0, fEtaMin, fEtaMax);
 573     fCache->GetGrid(step2)->SetRangeUser(0, fEtaMin, fEtaMax);
 574   }
 575   if (fPtMax > fPtMin && axis1 != 1 && axis2 != 1)
 576   {
 577     fCache->GetGrid(step1)->SetRangeUser(1, fPtMin, fPtMax);
 578     fCache->GetGrid(step2)->SetRangeUser(1, fPtMin, fPtMax);
 579   }
 580
 581   TH1* measured = 0;
 582   TH1* generated = 0;
 583
 584   if (axis2 > 0)
 585   {
 586     generated = fCache->Project(axis1, axis2, step1);
 587     measured = fCache->Project(axis1, axis2, step2);
 588   }
 589   else
 590   {
 591     generated = fCache->Project(axis1, step1);
 592     measured = fCache->Project(axis1, step2);
 593   }
 594
 595   measured->Divide(generated);
 596
 597   delete generated;
 598
 599   return measured;
 600 }
 601
 602 //____________________________________________________________________
 603 TH1* AliUEHist::GetEventEfficiency(CFStep step1, CFStep step2, Int_t axis1, Int_t axis2, Float_t ptLeadMin, Float_t ptLeadMax)
 604 {
 605   // creates a event-level efficiency by dividing step2 by step1
 606   // projected to axis1 and axis2 (optional if >= 0)
 607
 608   if (ptLeadMax > ptLeadMin)
 609   {
 610     fEventHist->GetGrid(step1)->SetRangeUser(0, ptLeadMin, ptLeadMax);
 611     fEventHist->GetGrid(step2)->SetRangeUser(0, ptLeadMin, ptLeadMax);
 612   }
 613
 614   TH1* measured = 0;
 615   TH1* generated = 0;
 616
 617   if (axis2 > 0)
 618   {
 619     generated = fEventHist->Project(axis1, axis2, step1);
 620     measured = fEventHist->Project(axis1, axis2, step2);
 621   }
 622   else
 623   {
 624     generated = fEventHist->Project(axis1, step1);
 625     measured = fEventHist->Project(axis1, step2);
 626   }
 627
 628   measured->Divide(generated);
 629
 630   delete generated;
 631
 632   if (ptLeadMax > ptLeadMin)
 633   {
 634     fEventHist->GetGrid(step1)->SetRangeUser(0, 0, -1);
 635     fEventHist->GetGrid(step2)->SetRangeUser(0, 0, -1);
 636   }
 637
 638   return measured;
 639 }
 640
 641 //____________________________________________________________________
 642 void AliUEHist::WeightHistogram(TH3* hist1, TH1* hist2)
 643 {
 644   // weights each entry of the 3d histogram hist1 with the 1d histogram hist2
 645   // where the matching is done of the z axis of hist1 with the x axis of hist2
 646
 647   if (hist1->GetNbinsZ() != hist2->GetNbinsX())
 648     AliFatal(Form("Inconsistent binning %d %d", hist1->GetNbinsZ(), hist2->GetNbinsX()));
 649
 650   for (Int_t x=1; x<=hist1->GetNbinsX(); x++)
 651   {
 652     for (Int_t y=1; y<=hist1->GetNbinsY(); y++)
 653     {
 654       for (Int_t z=1; z<=hist1->GetNbinsZ(); z++)
 655       {
 656         if (hist2->GetBinContent(z) > 0)
 657         {
 658           hist1->SetBinContent(x, y, z, hist1->GetBinContent(x, y, z) / hist2->GetBinContent(z));
 659           hist1->SetBinError(x, y, z, hist1->GetBinError(x, y, z) / hist2->GetBinContent(z));
 660         }
 661         else
 662         {
 663           hist1->SetBinContent(x, y, z, 0);
 664           hist1->SetBinError(x, y, z, 0);
 665         }
 666       }
 667     }
 668   }
 669 }
 670
 671 //____________________________________________________________________
 672 TH1* AliUEHist::GetBias(CFStep step1, CFStep step2, const char* axis, Float_t leadPtMin, Float_t leadPtMax)
 673 {
 674   // extracts the track-level bias (integrating out the multiplicity) between two steps (dividing step2 by step1)
 675   // done by weighting the track-level distribution with the number of events as function of leading pT
 676   // and then calculating the ratio between the distributions
 677   // projected to axis which is a TH3::Project3D string, e.g. "x", or "yx"
 678   //   no projection is done if axis == 0
 679
 680   // integrate over regions
 681   AliCFContainer* tmp = (AliCFContainer*) fTrackHist[0]->Clone();
 682   for (Int_t i = 1; i < fkRegions; i++)
 683     if (fTrackHist[i])
 684       tmp->Add(fTrackHist[i]);
 685
 686   TH1D* events1 = fEventHist->Project(0, step1);
 687   TH3D* hist1 = tmp->Project(0, 1, 2, step1);
 688   WeightHistogram(hist1, events1);
 689
 690   TH1D* events2 = fEventHist->Project(0, step2);
 691   TH3D* hist2 = tmp->Project(0, 1, 2, step2);
 692   WeightHistogram(hist2, events2);
 693
 694   TH1* generated = hist1;
 695   TH1* measured = hist2;
 696
 697   if (axis)
 698   {
 699     if (leadPtMax > leadPtMin)
 700     {
 701       hist1->GetZaxis()->SetRangeUser(leadPtMin, leadPtMax);
 702       hist2->GetZaxis()->SetRangeUser(leadPtMin, leadPtMax);
 703     }
 704
 705     if (fEtaMax > fEtaMin && !TString(axis).Contains("x"))
 706     {
 707       hist1->GetXaxis()->SetRangeUser(fEtaMin, fEtaMax);
 708       hist2->GetXaxis()->SetRangeUser(fEtaMin, fEtaMax);
 709     }
 710     if (fPtMax > fPtMin && !TString(axis).Contains("y"))
 711     {
 712       hist1->GetYaxis()->SetRangeUser(fPtMin, fPtMax);
 713       hist2->GetYaxis()->SetRangeUser(fPtMin, fPtMax);
 714     }
 715
 716     generated = hist1->Project3D(axis);
 717     measured  = hist2->Project3D(axis);
 718
 719     // delete hists here if projection has been done
 720     delete hist1;
 721     delete hist2;
 722   }
 723
 724   measured->Divide(generated);
 725
 726   delete events1;
 727   delete events2;
 728   delete generated;
 729   delete tmp;
 730
 731   return measured;
 732 }
 733
 734 //____________________________________________________________________
 735 TH2D* AliUEHist::GetTrackingEfficiency()
 736 {
 737   // extracts the tracking efficiency by calculating the efficiency from step kCFStepAnaTopology to kCFStepTrackedOnlyPrim
 738   // integrates over the regions and all other variables than pT and eta to increase the statistics
 739   //
 740   // returned histogram has to be deleted by the user
 741
 742   return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepAnaTopology, kCFStepTrackedOnlyPrim, 0, 1));
 743 }
 744
 745 //____________________________________________________________________
 746 TH1D* AliUEHist::GetTrackingEfficiency(Int_t axis)
 747 {
 748   // extracts the tracking efficiency by calculating the efficiency from step kCFStepAnaTopology to kCFStepTrackedOnlyPrim
 749   // integrates over the regions and all other variables than pT (axis == 0) and eta (axis == 1) to increase the statistics
 750
 751   return dynamic_cast<TH1D*> (GetTrackEfficiency(kCFStepAnaTopology, kCFStepTrackedOnlyPrim, axis));
 752 }
 753
 754 //____________________________________________________________________
 755 TH2D* AliUEHist::GetTrackingCorrection()
 756 {
 757   // extracts the tracking correction by calculating the efficiency from step kCFStepAnaTopology to kCFStepTracked
 758   // integrates over the regions and all other variables than pT and eta to increase the statistics
 759   //
 760   // returned histogram has to be deleted by the user
 761
 762   return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepTracked, kCFStepAnaTopology, 0, 1));
 763 }
 764
 765 //____________________________________________________________________
 766 TH1D* AliUEHist::GetTrackingCorrection(Int_t axis)
 767 {
 768   // extracts the tracking correction by calculating the efficiency from step kCFStepAnaTopology to kCFStepTracked
 769   // integrates over the regions and all other variables than pT (axis == 0) and eta (axis == 1) to increase the statistics
 770
 771   return dynamic_cast<TH1D*> (GetTrackEfficiency(kCFStepTracked, kCFStepAnaTopology, axis));
 772 }
 773
 774 //____________________________________________________________________
 775 TH2D* AliUEHist::GetTrackingContamination()
 776 {
 777   // extracts the tracking contamination by secondaries by calculating the efficiency from step kCFStepTrackedOnlyPrim to kCFStepTracked
 778   // integrates over the regions and all other variables than pT and eta to increase the statistics
 779   //
 780   // returned histogram has to be deleted by the user
 781
 782   return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepTracked, kCFStepTrackedOnlyPrim, 0, 1));
 783 }
 784
 785 //____________________________________________________________________
 786 TH1D* AliUEHist::GetTrackingContamination(Int_t axis)
 787 {
 788   // extracts the tracking contamination by secondaries by calculating the efficiency from step kCFStepTrackedOnlyPrim to kCFStepTracked
 789   // integrates over the regions and all other variables than pT (axis == 0) and eta (axis == 1) to increase the statistics
 790
 791   return dynamic_cast<TH1D*> (GetTrackEfficiency(kCFStepTracked, kCFStepTrackedOnlyPrim, axis));
 792 }
 793
 794 //____________________________________________________________________
 795 const char* AliUEHist::GetRegionTitle(Region region)
 796 {
 797   // returns the name of the given region
 798
 799   switch (region)
 800   {
 801     case kToward:
 802       return "Towards";
 803     case kAway:
 804       return "Away";
 805     case kMin:
 806       return (fCombineMinMax) ? "Transverse" : "Min";
 807     case kMax:
 808       return "Max";
 809   }
 810
 811   return 0;
 812 }
 813
 814 //____________________________________________________________________
 815 const char* AliUEHist::GetStepTitle(CFStep step)
 816 {
 817   // returns the name of the given step
 818
 819   switch (step)
 820   {
 821     case kCFStepAll:
 822       return "All events";
 823     case kCFStepTriggered:
 824       return "Triggered";
 825     case kCFStepVertex:
 826       return "Primary Vertex";
 827     case kCFStepAnaTopology:
 828       return "Required analysis topology";
 829     case kCFStepTrackedOnlyPrim:
 830       return "Tracked (matched MC, only primaries)";
 831     case kCFStepTracked:
 832       return "Tracked (matched MC, all)";
 833     case kCFStepReconstructed:
 834       return "Reconstructed";
 835     case kCFStepRealLeading:
 836       return "Correct leading particle identified";
 837     case kCFStepBiasStudy:
 838       return "Bias study applying tracking efficiency";
 839     case kCFStepBiasStudy2:
 840       return "Bias study applying tracking efficiency in two steps";
 841   }
 842
 843   return 0;
 844 }