]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG4/JetTasks/AliUEHist.cxx
updated dphi correlation analysis
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliUEHist.cxx
index 860d2640a22623ce8efe91ba7ec42e1adf38b1e4..f74935d8b84ed93311e3250592fb2b5b3f91ae90 100644 (file)
@@ -34,6 +34,7 @@
 #include "AliLog.h"
 #include "TCanvas.h"
 #include "TF1.h"
+#include "AliTHn.h"
 
 ClassImp(AliUEHist)
 
@@ -48,13 +49,18 @@ AliUEHist::AliUEHist(const char* reqHist) :
   fEtaMax(0),
   fPtMin(0),
   fPtMax(0),
+  fCentralityMin(0),
+  fCentralityMax(0),
+  fZVtxMin(0),
+  fZVtxMax(0),
   fContaminationEnhancement(0),
   fCombineMinMax(0),
-  fCache(0)
+  fCache(0),
+  fHistogramType(reqHist)
 {
   // Constructor
     
-  for (Int_t i=0; i<fkRegions; i++)
+  for (UInt_t i=0; i<fkRegions; i++)
     fTrackHist[i] = 0;
     
   if (strlen(reqHist) == 0)
@@ -67,9 +73,9 @@ AliUEHist::AliUEHist(const char* reqHist) :
     
   // track level
   Int_t nTrackVars = 4; // eta vs pT vs pT,lead (vs delta phi) vs multiplicity
-  Int_t iTrackBin[5];
-  Double_t* trackBins[5];
-  const char* trackAxisTitle[5];
+  Int_t iTrackBin[6];
+  Double_t* trackBins[6];
+  const char* trackAxisTitle[6];
   
   // eta
   iTrackBin[0] = 20;
@@ -80,14 +86,14 @@ AliUEHist::AliUEHist(const char* reqHist) :
   trackAxisTitle[0] = "#eta";
   
   // delta eta
-  const Int_t kNDeltaEtaBins = 20;
+  const Int_t kNDeltaEtaBins = 40;
   Double_t deltaEtaBins[kNDeltaEtaBins+1];
   for (Int_t i=0; i<=kNDeltaEtaBins; i++)
-    deltaEtaBins[i] = -2.0 + 0.2 * i;
+    deltaEtaBins[i] = -2.0 + 0.1 * i;
   
   // pT
-  iTrackBin[1] = 39;
-  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};
+  iTrackBin[1] = 22;
+  Double_t pTBins[] = {0.15, 0.2, 0.3, 0.4, 0.5, 0.75, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 12.0, 15.0, 20.0};
   trackBins[1] = pTBins;
   trackAxisTitle[1] = "p_{T} (GeV/c)";
   
@@ -98,14 +104,14 @@ AliUEHist::AliUEHist(const char* reqHist) :
     leadingpTBins[i] = 0.5 * i;
   
   // pT,lead binning 2
-  const Int_t kNLeadingpTBins2 = 13;
-  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 };
+  const Int_t kNLeadingpTBins2 = 9;
+  Double_t leadingpTBins2[] = { 0.5, 1.0, 2.0, 3.0, 4.0, 6.0, 8.0, 10.0, 15.0, 20.0 };
   
-  // phi,lead
-  const Int_t kNLeadingPhiBins = 40;
+  // phi,lead; this binning starts at -pi/2 and is modulo 3
+  const Int_t kNLeadingPhiBins = 72;
   Double_t leadingPhiBins[kNLeadingPhiBins+1];
   for (Int_t i=0; i<=kNLeadingPhiBins; i++)
-    leadingPhiBins[i] = -0.5 * TMath::Pi() + 1.0 / 40 * i * TMath::TwoPi();
+    leadingPhiBins[i] = -TMath::Pi() / 2 + 1.0 / kNLeadingPhiBins * i * TMath::TwoPi();
     
   // multiplicity
   const Int_t kNMultiplicityBins = 15;
@@ -114,19 +120,21 @@ AliUEHist::AliUEHist(const char* reqHist) :
     multiplicityBins[i] = -0.5 + i;
   multiplicityBins[kNMultiplicityBins] = 200;
   
-  // centrality
-  const Int_t kNCentralityBins = 15;
-  Double_t centralityBins[kNCentralityBins+1];
-  for (Int_t i=0; i<=kNCentralityBins; i++)
-    centralityBins[i] = -0.5 + i * 500;
+  trackBins[3] = multiplicityBins;
+  iTrackBin[3] = kNMultiplicityBins;
+  trackAxisTitle[3] = "multiplicity";
+  
+  // centrality (in %)
+  const Int_t kNCentralityBins = 5 + 1 + 9;
+  Double_t centralityBins[] = { 0, 1, 2, 3, 4, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100.1 };
   
   // particle species
   const Int_t kNSpeciesBins = 4; // pi, K, p, rest
   Double_t speciesBins[] = { -0.5, 0.5, 1.5, 2.5, 3.5 };
   
-  trackBins[3] = multiplicityBins;
-  iTrackBin[3] = kNMultiplicityBins;
-  trackAxisTitle[3] = "multiplicity";
+  // vtx-z axis
+  const Int_t kNVertexBins = 7;
+  Double_t vertexBins[] = { -7, -5, -3, -1, 1, 3, 5, 7 };
   
   // selection depending on requested histogram
   Int_t axis = -1; // 0 = pT,lead, 1 = phi,lead
@@ -153,7 +161,7 @@ AliUEHist::AliUEHist(const char* reqHist) :
   else
     AliFatal(Form("Invalid histogram requested: %s", reqHist));
   
-  Int_t initRegions = fkRegions;
+  UInt_t initRegions = fkRegions;
   
   if (axis == 0)
   {
@@ -177,7 +185,7 @@ AliUEHist::AliUEHist(const char* reqHist) :
   }
   else if (axis == 2)
   {
-    nTrackVars = 5;
+    nTrackVars = 6;
     initRegions = 1;
   
     iTrackBin[0] = kNDeltaEtaBins;
@@ -194,12 +202,19 @@ AliUEHist::AliUEHist(const char* reqHist) :
   
     iTrackBin[4] = kNLeadingPhiBins;
     trackBins[4] = leadingPhiBins;
-    trackAxisTitle[4] = "#Delta#phi";
+    trackAxisTitle[4] = "#Delta#phi (rad.)";
+
+    iTrackBin[5] = kNVertexBins;
+    trackBins[5] = vertexBins;
+    trackAxisTitle[5] = "z-vtx (cm)";
   }
     
-  for (Int_t i=0; i<initRegions; i++)
+  for (UInt_t i=0; i<initRegions; i++)
   {
-    fTrackHist[i] = new AliCFContainer(Form("fTrackHist_%d", i), title, fgkCFSteps, nTrackVars, iTrackBin);
+    if (strcmp(reqHist, "NumberDensityPhiCentrality") == 0)
+      fTrackHist[i] = new AliTHn(Form("fTrackHist_%d", i), title, fgkCFSteps, nTrackVars, iTrackBin);
+    else
+      fTrackHist[i] = new AliCFContainer(Form("fTrackHist_%d", i), title, fgkCFSteps, nTrackVars, iTrackBin);
     
     for (Int_t j=0; j<nTrackVars; j++)
     {
@@ -210,8 +225,22 @@ AliUEHist::AliUEHist(const char* reqHist) :
     SetStepNames(fTrackHist[i]);
   }
   
+  // event level
+  Int_t nEventVars = 2;
+  Int_t iEventBin[3];
+
   // track 3rd and 4th axis --> event 1st and 2nd axis
-  fEventHist = new AliCFContainer("fEventHist", title, fgkCFSteps, 2, iTrackBin+2);
+  iEventBin[0] = iTrackBin[2];
+  iEventBin[1] = iTrackBin[3];
+  
+  // plus track 5th axis (in certain cases)
+  if (axis == 2)
+  {
+    nEventVars = 3;
+    iEventBin[2] = iTrackBin[5];
+  }
+  
+  fEventHist = new AliCFContainer("fEventHist", title, fgkCFSteps, nEventVars, iEventBin);
   
   fEventHist->SetBinLimits(0, trackBins[2]);
   fEventHist->SetVarTitle(0, trackAxisTitle[2]);
@@ -219,17 +248,25 @@ AliUEHist::AliUEHist(const char* reqHist) :
   fEventHist->SetBinLimits(1, trackBins[3]);
   fEventHist->SetVarTitle(1, trackAxisTitle[3]);
   
+  if (axis == 2)
+  {
+    fEventHist->SetBinLimits(2, trackBins[5]);
+    fEventHist->SetVarTitle(2, trackAxisTitle[5]);
+  }
+
   SetStepNames(fEventHist);
   
   iTrackBin[2] = kNSpeciesBins;
 
-  fTrackHistEfficiency = new AliCFContainer("fTrackHistEfficiency", "Tracking efficiency", 3, 3, iTrackBin);
+  fTrackHistEfficiency = new AliCFContainer("fTrackHistEfficiency", "Tracking efficiency", 3, 4, iTrackBin);
   fTrackHistEfficiency->SetBinLimits(0, trackBins[0]);
   fTrackHistEfficiency->SetVarTitle(0, trackAxisTitle[0]);
   fTrackHistEfficiency->SetBinLimits(1, trackBins[1]);
   fTrackHistEfficiency->SetVarTitle(1, trackAxisTitle[1]);
   fTrackHistEfficiency->SetBinLimits(2, speciesBins);
   fTrackHistEfficiency->SetVarTitle(2, "particle species");
+  fTrackHistEfficiency->SetBinLimits(3, trackBins[3]);
+  fTrackHistEfficiency->SetVarTitle(3, trackAxisTitle[3]);
 }
 
 //_____________________________________________________________________________
@@ -242,9 +279,14 @@ AliUEHist::AliUEHist(const AliUEHist &c) :
   fEtaMax(0),
   fPtMin(0),
   fPtMax(0),
+  fCentralityMin(0),
+  fCentralityMax(0),
+  fZVtxMin(0),
+  fZVtxMax(0),
   fContaminationEnhancement(0),
   fCombineMinMax(0),
-  fCache(0)
+  fCache(0),
+  fHistogramType()
 {
   //
   // AliUEHist copy constructor
@@ -267,7 +309,7 @@ AliUEHist::~AliUEHist()
 {
   // Destructor
   
-  for (Int_t i=0; i<fkRegions; i++)
+  for (UInt_t i=0; i<fkRegions; i++)
   {
     if (fTrackHist[i])
     {
@@ -313,7 +355,7 @@ void AliUEHist::Copy(TObject& c) const
 
   AliUEHist& target = (AliUEHist &) c;
 
-  for (Int_t i=0; i<fkRegions; i++)
+  for (UInt_t i=0; i<fkRegions; i++)
     if (fTrackHist[i])
       target.fTrackHist[i] = dynamic_cast<AliCFContainer*> (fTrackHist[i]->Clone());
 
@@ -322,6 +364,21 @@ void AliUEHist::Copy(TObject& c) const
   
   if (fTrackHistEfficiency)
     target.fTrackHistEfficiency = dynamic_cast<AliCFContainer*> (fTrackHistEfficiency->Clone());
+    
+  target.fEtaMin = fEtaMin;
+  target.fEtaMax = fEtaMax;
+  target.fPtMin = fPtMin;
+  target.fPtMax = fPtMax;
+  target.fCentralityMin = fCentralityMin;
+  target.fCentralityMax = fCentralityMax;
+  target.fZVtxMin = fZVtxMin;
+  target.fZVtxMax = fZVtxMax;
+  
+  if (fContaminationEnhancement)
+    target.fContaminationEnhancement = dynamic_cast<TH1F*> (fContaminationEnhancement->Clone());
+    
+  target.fCombineMinMax = fCombineMinMax;
+  target.fHistogramType = fHistogramType;
 }
 
 //____________________________________________________________________
@@ -341,10 +398,10 @@ Long64_t AliUEHist::Merge(TCollection* list)
   TObject* obj;
 
   // collections of objects
-  const Int_t kMaxLists = fkRegions+2;
+  const UInt_t kMaxLists = fkRegions+2;
   TList** lists = new TList*[kMaxLists];
   
-  for (Int_t i=0; i<kMaxLists; i++)
+  for (UInt_t i=0; i<kMaxLists; i++)
     lists[i] = new TList;
   
   Int_t count = 0;
@@ -354,7 +411,7 @@ Long64_t AliUEHist::Merge(TCollection* list)
     if (entry == 0) 
       continue;
 
-    for (Int_t i=0; i<fkRegions; i++)
+    for (UInt_t i=0; i<fkRegions; i++)
       if (entry->fTrackHist[i])
         lists[i]->Add(entry->fTrackHist[i]);
     
@@ -363,14 +420,14 @@ Long64_t AliUEHist::Merge(TCollection* list)
 
     count++;
   }
-  for (Int_t i=0; i<fkRegions; i++)
+  for (UInt_t i=0; i<fkRegions; i++)
     if (fTrackHist[i])
       fTrackHist[i]->Merge(lists[i]);
   
   fEventHist->Merge(lists[fkRegions]);
   fTrackHistEfficiency->Merge(lists[fkRegions+1]);
 
-  for (Int_t i=0; i<kMaxLists; i++)
+  for (UInt_t i=0; i<kMaxLists; i++)
     delete lists[i];
     
   delete[] lists;
@@ -434,7 +491,7 @@ void AliUEHist::CountEmptyBins(AliUEHist::CFStep step, Float_t ptLeadMin, Float_
   // start from multiplicity 1
   binBegin[3] = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(3)->FindBin(1);
   
-  for (Int_t region=0; region<fkRegions; region++)
+  for (UInt_t region=0; region<fkRegions; region++)
   {
     Int_t total = 0;
     Int_t count = 0;
@@ -477,18 +534,33 @@ void AliUEHist::CountEmptyBins(AliUEHist::CFStep step, Float_t ptLeadMin, Float_
 }  
 
 //____________________________________________________________________
-TH1* AliUEHist::GetUEHist(AliUEHist::CFStep step, AliUEHist::Region region, Float_t ptLeadMin, Float_t ptLeadMax, Int_t multBinBegin, Int_t multBinEnd, Bool_t twoD)
+TH1* AliUEHist::GetUEHist(AliUEHist::CFStep step, AliUEHist::Region region, Float_t ptLeadMin, Float_t ptLeadMax, Int_t multBinBegin, Int_t multBinEnd, Int_t twoD, Bool_t etaNorm)
 {
   // Extracts the UE histogram at the given step and in the given region by projection and dividing tracks by events
   //
   // ptLeadMin, ptLeadMax: Only meaningful for vs. delta phi plot (third axis is ptLead)
   // Histogram has to be deleted by the caller of the function
+  //
+  // twoD: 0: 1D histogram as function of phi
+  //       1: 2D histogram, deltaphi, deltaeta
+  //       10: 1D histogram, within |deltaeta| < 0.8
+  //       11: 1D histogram, within 0.8 < |deltaeta| < 1.6
+  //
+  // etaNorm: if kTRUE (default), the distributions are divided by the area in delta eta
   
   // unzoom all axes
   ResetBinLimits(fTrackHist[region]->GetGrid(step));
   ResetBinLimits(fEventHist->GetGrid(step));
   
   SetBinLimits(fTrackHist[region]->GetGrid(step));
+  
+  // z vtx
+  if (fZVtxMax > fZVtxMin)
+  {
+    Printf("Using z-vtx range %f --> %f", fZVtxMin, fZVtxMax);
+    fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(5)->SetRangeUser(fZVtxMin, fZVtxMax);
+    fEventHist->GetGrid(step)->GetGrid()->GetAxis(2)->SetRangeUser(fZVtxMin, fZVtxMax);
+  }
     
   TH1D* tracks = 0;
   
@@ -512,7 +584,9 @@ TH1* AliUEHist::GetUEHist(AliUEHist::CFStep step, AliUEHist::Region region, Floa
     Float_t phiRegion = TMath::TwoPi() / 3;
     if (!fCombineMinMax && region == kMin)
       phiRegion /= 2;
-    Float_t normalization = phiRegion * (axis->GetBinUpEdge(axis->GetLast()) - axis->GetBinLowEdge(axis->GetFirst()));
+    Float_t normalization = phiRegion;
+    if (etaNorm)
+      normalization *= (axis->GetBinUpEdge(axis->GetLast()) - axis->GetBinLowEdge(axis->GetFirst()));
     //Printf("Normalization: %f", normalization);
     tracks->Scale(1.0 / normalization);
     
@@ -521,50 +595,83 @@ TH1* AliUEHist::GetUEHist(AliUEHist::CFStep step, AliUEHist::Region region, Floa
   }
   else
   {
-    // 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
-    // therefore the number density must be calculated here per leading pT bin and then summed
-  
-    if (multBinEnd > multBinBegin)
+    if (multBinEnd >= multBinBegin)
+    {
       Printf("Using multiplicity range %d --> %d", multBinBegin, multBinEnd);
-    fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(3)->SetRange(multBinBegin, multBinEnd);
-    fEventHist->GetGrid(step)->GetGrid()->GetAxis(1)->SetRange(multBinBegin, multBinEnd);
+      fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(3)->SetRange(multBinBegin, multBinEnd);
+      fEventHist->GetGrid(step)->GetGrid()->GetAxis(1)->SetRange(multBinBegin, multBinEnd);
+    }
     
     Int_t firstBin = fTrackHist[region]->GetAxis(2, step)->FindBin(ptLeadMin);
     Int_t lastBin = fTrackHist[region]->GetAxis(2, step)->FindBin(ptLeadMax);
     
-    TH1D* events = fEventHist->ShowProjection(0, step);
+    Printf("Using leading pT range %d --> %d", firstBin, lastBin);
     
-    for (Int_t bin=firstBin; bin<=lastBin; bin++)
-    {
-      fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(2)->SetRange(bin, bin);
-      TH1D* tracksTmp = 0;
-      if (twoD)
-        tracksTmp = (TH1D*) fTrackHist[region]->GetGrid(step)->Project(4, 0);
-      else
-        tracksTmp = (TH1D*) fTrackHist[region]->GetGrid(step)->Project(4);
-      Printf("Calculated histogram in bin %d --> %f tracks", bin, tracksTmp->Integral());
-      fTrackHist[region]->GetGrid(step)->SetRangeUser(2, 0, -1);
+    fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(2)->SetRange(firstBin, lastBin);
+    
+    if (twoD == 0)
+      tracks = (TH1D*) fTrackHist[region]->GetGrid(step)->Project(4);
+    else
+      tracks = (TH1D*) fTrackHist[region]->GetGrid(step)->Project(4, 0);
       
-      // normalize to get a density (deta dphi)
+    Printf("Calculated histogram --> %f tracks", tracks->Integral());
+    fTrackHist[region]->GetGrid(step)->SetRangeUser(2, 0, -1);
+    
+    if (twoD == 10 || twoD == 11)
+    {
+      Float_t etaLimit = 0.8;
+      if (twoD == 10)
+      {
+        tracks = (TH1D*) ((TH2*) tracks)->ProjectionX("proj", tracks->GetYaxis()->FindBin(-etaLimit + 0.01), tracks->GetYaxis()->FindBin(etaLimit - 0.01))->Clone();
+        
+        // TODO calculate acc with 2 * (deta - 0.5 * deta*deta / 1.6)
+        tracks->Scale(1. / 0.75);
+      }
+      else if (twoD == 11)
+      {
+        TH1D* tracksTmp1 = (TH1D*) ((TH2*) tracks)->ProjectionX("proj1", tracks->GetYaxis()->FindBin(-etaLimit * 2 + 0.01), tracks->GetYaxis()->FindBin(-etaLimit - 0.01))->Clone();
+        TH1D* tracksTmp2 = ((TH2*) tracks)->ProjectionX("proj2", tracks->GetYaxis()->FindBin(etaLimit + 0.01), tracks->GetYaxis()->FindBin(2 * etaLimit - 0.01));
+        
+        tracksTmp1->Add(tracksTmp2);
+        
+        delete tracks;
+        tracks = tracksTmp1;
+        delete tracksTmp2;
+        
+        tracks->Scale(1. / 0.25);
+      }
+    }
+    
+    // normalize to get a density (deta dphi)
+    // TODO normalization may be off for 2d histogram
+    Float_t normalization = fTrackHist[region]->GetGrid(step)->GetAxis(4)->GetBinWidth(1);
+    
+    if (etaNorm)
+    {
       TAxis* axis = fTrackHist[region]->GetGrid(step)->GetAxis(0);
-      Float_t normalization = fTrackHist[region]->GetGrid(step)->GetAxis(4)->GetBinWidth(1) * (axis->GetBinUpEdge(axis->GetLast()) - axis->GetBinLowEdge(axis->GetFirst()));
-      //Printf("Normalization: %f", normalization);
-      tracksTmp->Scale(1.0 / normalization);
-      
-      Int_t nEvents = (Int_t) events->GetBinContent(bin);
-      if (nEvents > 0)
-        tracksTmp->Scale(1.0 / nEvents);
-      Printf("Calculated histogram in bin %d --> %d events", bin, nEvents);
-      
-      if (!tracks)
-        tracks = tracksTmp;
+      if (strcmp(axis->GetTitle(), "#eta") == 0)
+      {
+       Printf("Normalizing using eta-axis range");
+       normalization *= axis->GetBinUpEdge(axis->GetLast()) - axis->GetBinLowEdge(axis->GetFirst());
+      }
       else
       {
-        tracks->Add(tracksTmp);
-        delete tracksTmp;
+       Printf("Normalizing assuming accepted range of |eta| < 0.8");
+       normalization *= 0.8 * 2;
       }
     }
     
+    //Printf("Normalization: %f", normalization);
+    tracks->Scale(1.0 / normalization);
+    
+    // NOTE fEventHist contains the number of events for the underlying event analysis and the number of trigger particles for the azimuthal correlation analysis. In the latter case the naming is therefore somewhat misleading!
+    TH1D* events = fEventHist->ShowProjection(0, step);
+    Int_t nEvents = (Int_t) events->Integral(firstBin, lastBin);
+    Printf("Calculated histogram --> %d events", nEvents);
+      
+    if (nEvents > 0)
+      tracks->Scale(1.0 / nEvents);
+    
     delete events;
   }
 
@@ -574,6 +681,98 @@ TH1* AliUEHist::GetUEHist(AliUEHist::CFStep step, AliUEHist::Region region, Floa
   return tracks;
 }
 
+//____________________________________________________________________
+TH1* AliUEHist::GetPtHist(CFStep step, Region region, Float_t ptLeadMin, Float_t ptLeadMax, Int_t multBinBegin, Int_t multBinEnd, Float_t phiMin, Float_t phiMax, Float_t etaMin, Float_t etaMax, Bool_t skipPhiNormalization)
+{
+  // returns a histogram projected to pT,assoc with the provided cuts
+  
+   // unzoom all axes
+  ResetBinLimits(fTrackHist[region]->GetGrid(step));
+  ResetBinLimits(fEventHist->GetGrid(step));
+  
+  TH1D* tracks = 0;
+  
+    // 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
+  // therefore the number density must be calculated here per leading pT bin and then summed
+
+  if (multBinEnd > multBinBegin)
+    Printf("Using multiplicity range %d --> %d", multBinBegin, multBinEnd);
+  fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(3)->SetRange(multBinBegin, multBinEnd);
+  fEventHist->GetGrid(step)->GetGrid()->GetAxis(1)->SetRange(multBinBegin, multBinEnd);
+  
+  fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(0)->SetRangeUser(etaMin, etaMax);
+  fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(4)->SetRangeUser(phiMin, phiMax);
+  
+  // get real phi cuts due to binning
+  Float_t phiMinReal = fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(4)->GetBinLowEdge(fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(4)->GetFirst());
+  Float_t phiMaxReal = fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(4)->GetBinUpEdge(fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(4)->GetLast());
+  Printf("phi min = %.2f and max = %.2f requested; due to binning range will be from %.2f to %.2f, i.e. a %.0f%% larger window", phiMin, phiMax, phiMinReal, phiMaxReal, (phiMaxReal - phiMinReal) / (phiMax - phiMin) * 100 - 100);
+  
+  Int_t firstBin = fTrackHist[region]->GetAxis(2, step)->FindBin(ptLeadMin);
+  Int_t lastBin = fTrackHist[region]->GetAxis(2, step)->FindBin(ptLeadMax);
+  
+  TH1D* events = fEventHist->ShowProjection(0, step);
+  
+  for (Int_t bin=firstBin; bin<=lastBin; bin++)
+  {
+    fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(2)->SetRange(bin, bin);
+    
+    // project to pT,assoc
+    TH1D* tracksTmp = (TH1D*) fTrackHist[region]->GetGrid(step)->Project(1);
+    
+    Printf("Calculated histogram in bin %d --> %f tracks", bin, tracksTmp->Integral());
+    fTrackHist[region]->GetGrid(step)->SetRangeUser(2, 0, -1);
+    
+    // normalize to get a density (deta dphi)
+    Float_t normalization = 1;
+    TAxis* axis = fTrackHist[region]->GetGrid(step)->GetAxis(0);
+    if (strcmp(axis->GetTitle(), "#eta") == 0)
+    {
+      Printf("Normalizing using eta-axis range");
+      normalization *= axis->GetBinUpEdge(axis->GetLast()) - axis->GetBinLowEdge(axis->GetFirst());
+    }
+    else
+    {
+      Printf("Normalizating assuming accepted range of |eta| < 0.8");
+      normalization *= 0.8 * 2;
+    }
+    
+    // dphi
+    if (!skipPhiNormalization)
+      normalization *= phiMaxReal - phiMinReal;
+    
+    //Printf("Normalization: %f", normalization);
+    tracksTmp->Scale(1.0 / normalization);
+    
+    // and now dpT (bins have different width)
+    for (Int_t i=1; i<=tracksTmp->GetNbinsX(); i++)
+    {
+      tracksTmp->SetBinContent(i, tracksTmp->GetBinContent(i) / tracksTmp->GetXaxis()->GetBinWidth(i));
+      tracksTmp->SetBinError  (i, tracksTmp->GetBinError(i) / tracksTmp->GetXaxis()->GetBinWidth(i));
+    }
+     
+    Int_t nEvents = (Int_t) events->GetBinContent(bin);
+    if (nEvents > 0)
+      tracksTmp->Scale(1.0 / nEvents);
+    Printf("Calculated histogram in bin %d --> %d events", bin, nEvents);
+    
+    if (!tracks)
+      tracks = tracksTmp;
+    else
+    {
+      tracks->Add(tracksTmp);
+      delete tracksTmp;
+    }
+  }
+  
+  delete events;
+
+  ResetBinLimits(fTrackHist[region]->GetGrid(step));
+  ResetBinLimits(fEventHist->GetGrid(step));
+
+  return tracks;
+}
+
 //____________________________________________________________________
 void AliUEHist::CorrectTracks(CFStep step1, CFStep step2, TH1* trackCorrection, Int_t var1, Int_t var2)
 {
@@ -583,7 +782,7 @@ void AliUEHist::CorrectTracks(CFStep step1, CFStep step2, TH1* trackCorrection,
   //
   // if trackCorrection is 0, just copies content from step1 to step2
   
-  for (Int_t region=0; region<fkRegions; region++)
+  for (UInt_t region=0; region<fkRegions; region++)
     CorrectTracks(step1, step2, region, trackCorrection, var1, var2);
 }
 
@@ -641,7 +840,7 @@ void AliUEHist::CorrectTracks(CFStep step1, CFStep step2, Int_t region, TH1* tra
 }
 
 //____________________________________________________________________
-void AliUEHist::CorrectEvents(CFStep step1, CFStep step2, TH1D* eventCorrection, Int_t var)
+void AliUEHist::CorrectEvents(CFStep step1, CFStep step2, TH1* eventCorrection, Int_t var1, Int_t var2)
 {
   // corrects from step1 to step2 by multiplying the events with eventCorrection
   // eventCorrection is as function of leading pT (var == 0) or multiplicity (var == 1)
@@ -654,8 +853,11 @@ void AliUEHist::CorrectEvents(CFStep step1, CFStep step2, TH1D* eventCorrection,
   // clear target histogram
   target->GetGrid()->Reset();
   
-  if (eventCorrection != 0 && grid->GetNBins(var) != eventCorrection->GetNbinsX())
-    AliFatal(Form("Invalid binning: %d %d", grid->GetNBins(var), eventCorrection->GetNbinsX()));
+  if (eventCorrection != 0 && grid->GetNBins(var1) != eventCorrection->GetNbinsX())
+    AliFatal(Form("Invalid binning: %d %d", grid->GetNBins(var1), eventCorrection->GetNbinsX()));
+  
+  if (eventCorrection != 0 && var2 != -1 && grid->GetNBins(var2) != eventCorrection->GetNbinsY())
+    AliFatal(Form("Invalid binning: %d %d", grid->GetNBins(var2), eventCorrection->GetNbinsY()));
   
   Int_t bins[2];
   for (Int_t x = 1; x <= grid->GetNBins(0); x++)
@@ -672,8 +874,16 @@ void AliUEHist::CorrectEvents(CFStep step1, CFStep step2, TH1D* eventCorrection,
         
         if (eventCorrection != 0)
         {
-          value *= eventCorrection->GetBinContent(bins[var]);
-          error *= eventCorrection->GetBinContent(bins[var]);
+          if (var2 == -1)
+          {
+            value *= eventCorrection->GetBinContent(bins[var1]);
+            error *= eventCorrection->GetBinContent(bins[var1]);
+          }
+          else
+          {
+            value *= eventCorrection->GetBinContent(bins[var1], bins[var2]);
+            error *= eventCorrection->GetBinContent(bins[var1], bins[var2]);
+          }
         }
         
         target->SetElement(bins, value);
@@ -692,113 +902,239 @@ void AliUEHist::Correct(AliUEHist* corrections)
   //
   // in this object the data is expected in the step kCFStepReconstructed
   
-  // ---- track level
+  if (fHistogramType.Length() == 0)
+  {
+    Printf("W-AliUEHist::Correct: fHistogramType not defined. Guessing histogram type...");
+    if (fTrackHist[kToward]->GetNVar() < 5)
+    {
+      if (strcmp(fTrackHist[kToward]->GetTitle(), "d^{2}N_{ch}/d#phid#eta") == 0)
+        fHistogramType = "NumberDensitypT";
+      else if (strcmp(fTrackHist[kToward]->GetTitle(), "d^{2}#Sigma p_{T}/d#phid#eta") == 0)
+        fHistogramType = "SumpT";
+    }
+    else if (fTrackHist[kToward]->GetNVar() == 5)
+    {
+      if (strcmp(fTrackHist[kToward]->GetGrid(0)->GetGrid()->GetAxis(3)->GetTitle(), "multiplicity") == 0)
+        fHistogramType = "NumberDensityPhi";
+      else if (strcmp(fTrackHist[kToward]->GetGrid(0)->GetGrid()->GetAxis(3)->GetTitle(), "centrality") == 0)
+        fHistogramType = "NumberDensityPhiCentrality";
+    }
+      
+    if (fHistogramType.Length() == 0)
+      AliFatal("Cannot figure out histogram type. Try setting it manually...");
+  }
+  
+  Printf("AliUEHist::Correct: Correcting %s...", fHistogramType.Data());
   
-  // bias due to migration in leading pT (because the leading particle is not reconstructed, and the subleading is used)
-  // extracted as function of leading pT
-  for (Int_t region = 0; region < fkRegions; region++)
+  if (strcmp(fHistogramType, "NumberDensitypT") == 0 || strcmp(fHistogramType, "NumberDensityPhi") == 0 || strcmp(fHistogramType, "SumpT") == 0) // the last is for backward compatibilty
   {
-    if (!fTrackHist[region])
-      continue;
-
+    // ---- track level
+    
+    // bias due to migration in leading pT (because the leading particle is not reconstructed, and the subleading is used)
+    // extracted as function of leading pT
+    Bool_t biasFromMC = kFALSE;
     const char* projAxis = "z";
     Int_t secondBin = -1;
 
-    if (fTrackHist[region]->GetNVar() == 5)
+    if (strcmp(fHistogramType, "NumberDensityPhi") == 0)
     {
       projAxis = "yz";
       secondBin = 4;
     }
     
-    #if 0
-      TH1* leadingBias = (TH1*) corrections->GetBias(kCFStepReconstructed, kCFStepTracked, region, projAxis); // from MC
-      Printf("WARNING: Using MC bias correction");
-    #else
-      TH1* leadingBias = (TH1*) GetBias(kCFStepBiasStudy, kCFStepReconstructed, region, projAxis);          // from data
-    #endif
-    CorrectTracks(kCFStepReconstructed, kCFStepTracked, region, leadingBias, 2, secondBin);
-    if (region == kMin && fCombineMinMax)
+    for (UInt_t region = 0; region < fkRegions; region++)
     {
-      CorrectTracks(kCFStepReconstructed, kCFStepTracked, kMax, leadingBias, 2, secondBin);
-      delete leadingBias;
-      break;
+      if (!fTrackHist[region])
+        continue;
+   
+      TH1* leadingBiasTracks =  0;
+      if (biasFromMC)
+      {
+        leadingBiasTracks = (TH1*) corrections->GetBias(kCFStepReconstructed, kCFStepTracked, region, projAxis, 0, -1, 1); // from MC
+        Printf("WARNING: Using MC bias correction");
+      }
+      else
+        leadingBiasTracks = (TH1*) GetBias(kCFStepBiasStudy, kCFStepReconstructed, region, projAxis, 0, -1, 1);          // from data
+        
+      CorrectTracks(kCFStepReconstructed, kCFStepTracked, region, leadingBiasTracks, 2, secondBin);
+      if (region == kMin && fCombineMinMax)
+      {
+        CorrectTracks(kCFStepReconstructed, kCFStepTracked, kMax, leadingBiasTracks, 2, secondBin);
+        delete leadingBiasTracks;
+        break;
+      }
+      delete leadingBiasTracks;
+    }
+    
+    TH1* leadingBiasEvents = 0;
+    if (biasFromMC)
+      leadingBiasEvents = (TH1*) corrections->GetBias(kCFStepReconstructed, kCFStepTracked, kToward, projAxis, 0, -1, 2); // from MC
+    else
+      leadingBiasEvents = (TH1*) GetBias(kCFStepBiasStudy, kCFStepReconstructed, kToward, projAxis, 0, -1, 2);          // from data
+      
+    //new TCanvas; leadingBiasEvents->DrawCopy();
+    CorrectEvents(kCFStepReconstructed, kCFStepTracked, leadingBiasEvents, 0);
+    delete leadingBiasEvents;
+    
+    // --- contamination correction ---
+    
+    TH2D* contamination = corrections->GetTrackingContamination();
+    if (corrections->fContaminationEnhancement)
+    {
+      Printf("Applying contamination enhancement");
+      
+      for (Int_t x=1; x<=contamination->GetNbinsX(); x++)
+        for (Int_t y=1; y<=contamination->GetNbinsX(); y++)
+          contamination->SetBinContent(x, y, contamination->GetBinContent(x, y) * corrections->fContaminationEnhancement->GetBinContent(corrections->fContaminationEnhancement->GetXaxis()->FindBin(contamination->GetYaxis()->GetBinCenter(y))));
+    }
+    CorrectTracks(kCFStepTracked, kCFStepTrackedOnlyPrim, contamination, 0, 1);
+    CorrectEvents(kCFStepTracked, kCFStepTrackedOnlyPrim, 0, 0);
+    delete contamination;
+    
+    // --- efficiency correction ---
+    // correct single-particle efficiency for associated particles
+    // in addition correct for efficiency on trigger particles (tracks AND events)
+    
+    TH1* efficiencyCorrection = corrections->GetTrackingEfficiencyCorrection();
+    // use kCFStepVertex as a temporary step (will be overwritten later anyway)
+    CorrectTracks(kCFStepTrackedOnlyPrim, kCFStepVertex, efficiencyCorrection, 0, 1);
+    delete efficiencyCorrection;
+    
+    // correct pT,T
+    efficiencyCorrection = corrections->GetTrackEfficiency(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, 1, -1, 2);
+    CorrectEvents(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, efficiencyCorrection, 0);
+    CorrectTracks(kCFStepVertex, kCFStepAnaTopology, efficiencyCorrection, 2);
+    delete efficiencyCorrection;
+    
+    // copy 
+    CorrectTracks(kCFStepAnaTopology, kCFStepVertex, 0, -1);
+    CorrectEvents(kCFStepAnaTopology, kCFStepVertex, 0, 0);
+    
+    // vertex correction on the level of events as function of multiplicity, weighting tracks and events with the same factor
+    // practically independent of low pT cut 
+    TH1D* vertexCorrection = (TH1D*) corrections->GetEventEfficiency(kCFStepVertex, kCFStepTriggered, 1);
+    
+    // convert stage from true multiplicity to observed multiplicity by simple conversion factor
+    TH1D* vertexCorrectionObs = (TH1D*) vertexCorrection->Clone("vertexCorrection2");
+    vertexCorrectionObs->Reset();
+    
+    TF1* func = new TF1("func", "[1]+[0]/(x-[2])");
+    // some defaults
+    func->SetParameters(0.1, 1, -0.7);
+    vertexCorrection->Fit(func, "0I", "", 0, 3);
+    for (Int_t i=1; i<=vertexCorrectionObs->GetNbinsX(); i++)
+    {
+      Float_t xPos = 1.0 / 0.77 * vertexCorrectionObs->GetXaxis()->GetBinCenter(i);
+      if (xPos < 3)
+        vertexCorrectionObs->SetBinContent(i, func->Eval(xPos));
+      else
+        vertexCorrectionObs->SetBinContent(i, vertexCorrection->Interpolate(xPos));
     }
-    delete leadingBias;
-  }
-  CorrectEvents(kCFStepReconstructed, kCFStepTracked, 0, 0);
   
-  // correct with kCFStepTracked --> kCFStepTrackedOnlyPrim
-  TH2D* contamination = corrections->GetTrackingContamination();
-  if (corrections->fContaminationEnhancement)
-  {
-    Printf("Applying contamination enhancement");
+    #if 0
+    new TCanvas;
+    vertexCorrection->DrawCopy();
+    vertexCorrectionObs->SetLineColor(2);
+    vertexCorrectionObs->DrawCopy("same");
+    func->SetRange(0, 4);
+    func->DrawClone("same");
+    #endif
     
-    for (Int_t x=1; x<=contamination->GetNbinsX(); x++)
-      for (Int_t y=1; y<=contamination->GetNbinsX(); y++)
-        contamination->SetBinContent(x, y, contamination->GetBinContent(x, y) * corrections->fContaminationEnhancement->GetBinContent(corrections->fContaminationEnhancement->GetXaxis()->FindBin(contamination->GetYaxis()->GetBinCenter(y))));
+    CorrectTracks(kCFStepVertex, kCFStepTriggered, vertexCorrectionObs, 3);
+    CorrectEvents(kCFStepVertex, kCFStepTriggered, vertexCorrectionObs, 1);
+    delete vertexCorrectionObs;
+    delete vertexCorrection;
+    delete func;
+    
+    // copy 
+    CorrectTracks(kCFStepTriggered, kCFStepAll, 0, -1);
+    CorrectEvents(kCFStepTriggered, kCFStepAll, 0, 0);
   }
-  CorrectTracks(kCFStepTracked, kCFStepTrackedOnlyPrim, contamination, 0, 1);
-  CorrectEvents(kCFStepTracked, kCFStepTrackedOnlyPrim, 0, 0);
-  delete contamination;
-  
-  // correct with kCFStepTrackedOnlyPrim --> kCFStepAnaTopology
-  TH2D* efficiencyCorrection = corrections->GetTrackingEfficiencyCorrection();
-  CorrectTracks(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, efficiencyCorrection, 0, 1);
-  CorrectEvents(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, 0, 0);
-  delete efficiencyCorrection;
-  
-  // copy 
-  CorrectTracks(kCFStepAnaTopology, kCFStepVertex, 0, -1);
-  CorrectEvents(kCFStepAnaTopology, kCFStepVertex, 0, 0);
-  
-  // vertex correction on the level of events as function of multiplicity, weighting tracks and events with the same factor
-  // practically independent of low pT cut 
-  TH1D* vertexCorrection = (TH1D*) corrections->GetEventEfficiency(kCFStepVertex, kCFStepTriggered, 1);
-  
-  // convert stage from true multiplicity to observed multiplicity by simple conversion factor
-  TH1D* vertexCorrectionObs = (TH1D*) vertexCorrection->Clone("vertexCorrection2");
-  vertexCorrectionObs->Reset();
-  
-  TF1* func = new TF1("func", "[1]+[0]/(x-[2])");
-  // some defaults
-  func->SetParameters(0.1, 1, -0.7);
-  vertexCorrection->Fit(func, "0I", "", 0, 3);
-  for (Int_t i=1; i<=vertexCorrectionObs->GetNbinsX(); i++)
+  else if (strcmp(fHistogramType, "NumberDensityPhiCentrality") == 0)
   {
-    Float_t xPos = 1.0 / 0.77 * vertexCorrectionObs->GetXaxis()->GetBinCenter(i);
-    if (xPos < 3)
-      vertexCorrectionObs->SetBinContent(i, func->Eval(xPos));
-    else
-      vertexCorrectionObs->SetBinContent(i, vertexCorrection->Interpolate(xPos));
+    // copy 
+    CorrectTracks(kCFStepReconstructed, kCFStepTracked, 0, -1);
+    CorrectEvents(kCFStepReconstructed, kCFStepTracked, 0, -1);
+    
+    // Dont use eta in the following, because it is a Delta-eta axis
+    
+    // contamination correction
+    // correct single-particle contamination for associated particles
+    
+    TH1* contamination = corrections->GetTrackingContamination(1);
+    
+    if (1)
+    {
+      Printf("Applying contamination enhancement");
+      
+      for (Int_t bin = 1; bin <= contamination->GetNbinsX(); bin++)
+      {
+        printf("%f", contamination->GetBinContent(bin));
+        if (contamination->GetBinContent(bin) > 0)
+          contamination->SetBinContent(bin, 1.0 + 1.1 * (contamination->GetBinContent(bin) - 1.0));
+        printf(" --> %f\n", contamination->GetBinContent(bin));
+      }
+    }
+      
+    CorrectTracks(kCFStepTracked, kCFStepTrackedOnlyPrim, contamination, 1);
+    delete contamination;    
+    
+    // correct for additional contamination due to trigger particle around phi ~ 0
+    TH2* correlatedContamination = corrections->GetCorrelatedContamination();
+    if (1)
+    {
+      Printf("Applying contamination enhancement");
+      
+      for (Int_t bin = 1; bin <= correlatedContamination->GetNbinsX(); bin++)
+        for (Int_t bin2 = 1; bin2 <= correlatedContamination->GetNbinsY(); bin2++)
+        {
+          printf("%f", correlatedContamination->GetBinContent(bin, bin2));
+          if (correlatedContamination->GetBinContent(bin, bin2) > 0)
+            correlatedContamination->SetBinContent(bin, bin2, 1.0 + 1.1 * (correlatedContamination->GetBinContent(bin, bin2) - 1.0));
+          printf(" --> %f\n", correlatedContamination->GetBinContent(bin, bin2));
+        }
+    }
+    
+    //new TCanvas; correlatedContamination->DrawCopy("COLZ");
+    CorrectCorrelatedContamination(kCFStepTrackedOnlyPrim, 0, correlatedContamination);
+    
+    delete correlatedContamination;
+    
+    // TODO correct for contamination of trigger particles (for tracks AND events)
+    CorrectEvents(kCFStepTracked, kCFStepTrackedOnlyPrim, 0, 0);
+    
+    // --- efficiency correction ---
+    // correct single-particle efficiency for associated particles
+    // in addition correct for efficiency on trigger particles (tracks AND events)
+    
+    // in bins of pT and centrality
+    TH1* efficiencyCorrection = corrections->GetTrackingEfficiencyCorrectionCentrality();
+    // use kCFStepAnaTopology as a temporary step 
+    CorrectTracks(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, efficiencyCorrection, 1, 3);
+    delete efficiencyCorrection;
+    
+    // correct pT,T in bins of pT and centrality
+    efficiencyCorrection = corrections->GetTrackEfficiency(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, 1, 3, 2);
+    CorrectEvents(kCFStepTrackedOnlyPrim, kCFStepVertex, efficiencyCorrection, 0, 1);
+    CorrectTracks(kCFStepAnaTopology, kCFStepVertex, efficiencyCorrection, 2, 3);
+    delete efficiencyCorrection;
+    
+    // no correction for vertex finding efficiency and trigger efficiency needed in PbPb
+    // copy 
+    CorrectTracks(kCFStepVertex, kCFStepAll, 0, -1);
+    CorrectEvents(kCFStepVertex, kCFStepAll, 0, -1);
   }
-  #if 1
-  new TCanvas;
-  vertexCorrection->DrawCopy();
-  vertexCorrectionObs->SetLineColor(2);
-  vertexCorrectionObs->DrawCopy("same");
-  func->SetRange(0, 4);
-  func->DrawClone("same");
-  #endif
-  
-  CorrectTracks(kCFStepVertex, kCFStepTriggered, vertexCorrectionObs, 3);
-  CorrectEvents(kCFStepVertex, kCFStepTriggered, vertexCorrectionObs, 1);
-  delete vertexCorrectionObs;
-  delete vertexCorrection;
-  delete func;
-  
-  // copy 
-  CorrectTracks(kCFStepTriggered, kCFStepAll, 0, -1);
-  CorrectEvents(kCFStepTriggered, kCFStepAll, 0, 0);
+  else
+    AliFatal(Form("Unknown histogram for correction: %s", GetTitle()));
 }
 
 //____________________________________________________________________
-TH1* AliUEHist::GetTrackEfficiency(CFStep step1, CFStep step2, Int_t axis1, Int_t axis2, Int_t source)
+TH1* AliUEHist::GetTrackEfficiency(CFStep step1, CFStep step2, Int_t axis1, Int_t axis2, Int_t source, Int_t axis3)
 {
   // creates a track-level efficiency by dividing step2 by step1
   // projected to axis1 and axis2 (optional if >= 0)
   //
-  // source: 0 = fTrackHist; 1 = fTrackHistEfficiency
+  // source: 0 = fTrackHist; 1 = fTrackHistEfficiency; 2 = fTrackHistEfficiency rebinned for pT,T / pT,lead binning
   
   // integrate over regions
   // cache it for efficiency (usually more than one efficiency is requested)
@@ -810,13 +1146,13 @@ TH1* AliUEHist::GetTrackEfficiency(CFStep step1, CFStep step2, Int_t axis1, Int_
     if (!fCache)
     {
       fCache = (AliCFContainer*) fTrackHist[0]->Clone();
-      for (Int_t i = 1; i < fkRegions; i++)
+      for (UInt_t i = 1; i < fkRegions; i++)
         if (fTrackHist[i])
           fCache->Add(fTrackHist[i]);
     }
     sourceContainer = fCache;
   }
-  else if (source == 1)
+  else if (source == 1 || source == 2)
   {
     sourceContainer = fTrackHistEfficiency;
     // step offset because we start with kCFStepAnaTopology
@@ -826,43 +1162,55 @@ TH1* AliUEHist::GetTrackEfficiency(CFStep step1, CFStep step2, Int_t axis1, Int_
   else
     return 0;
         
-  // reset all limits and set the right ones except those in axis1 and axis2
+  // reset all limits and set the right ones except those in axis1, axis2 and axis3
   ResetBinLimits(sourceContainer->GetGrid(step1));
   ResetBinLimits(sourceContainer->GetGrid(step2));
-  if (fEtaMax > fEtaMin && axis1 != 0 && axis2 != 0)
+  if (fEtaMax > fEtaMin && axis1 != 0 && axis2 != 0 && axis3 != 0)
   {
     sourceContainer->GetGrid(step1)->SetRangeUser(0, fEtaMin, fEtaMax);
     sourceContainer->GetGrid(step2)->SetRangeUser(0, fEtaMin, fEtaMax);
   }
-  if (fPtMax > fPtMin && axis1 != 1 && axis2 != 1)
+  if (fPtMax > fPtMin && axis1 != 1 && axis2 != 1 && axis3 != 1)
   {
     sourceContainer->GetGrid(step1)->SetRangeUser(1, fPtMin, fPtMax);
     sourceContainer->GetGrid(step2)->SetRangeUser(1, fPtMin, fPtMax);
   }
+  if (fCentralityMax > fCentralityMin && axis1 != 3 && axis2 != 3 && axis3 != 3)
+  {
+    sourceContainer->GetGrid(step1)->SetRangeUser(3, fCentralityMin, fCentralityMax);
+    sourceContainer->GetGrid(step2)->SetRangeUser(3, fCentralityMin, fCentralityMax);
+  }
   
   TH1* measured = 0;
   TH1* generated = 0;
     
-  if (axis2 >= 0)
+  if (axis3 >= 0)
   {
-    generated = sourceContainer->Project(axis1, axis2, step1);
-    measured = sourceContainer->Project(axis1, axis2, step2);
+    generated = sourceContainer->Project(step1, axis1, axis2, axis3);
+    measured = sourceContainer->Project(step2, axis1, axis2, axis3);
+  }
+  else if (axis2 >= 0)
+  {
+    generated = sourceContainer->Project(step1, axis1, axis2);
+    measured = sourceContainer->Project(step2, axis1, axis2);
   }
   else
   {
-    generated = sourceContainer->Project(axis1, step1);
-    measured = sourceContainer->Project(axis1, step2);
+    generated = sourceContainer->Project(step1, axis1);
+    measured = sourceContainer->Project(step2, axis1);
   }
   
-  // check for bins with less than 100 entries, print warning
-  Int_t binBegin[2];
-  Int_t binEnd[2];
+  // check for bins with less than 50 entries, print warning
+  Int_t binBegin[3];
+  Int_t binEnd[3];
   
   binBegin[0] = 1;
   binBegin[1] = 1;
+  binBegin[2] = 1;
   
   binEnd[0] = generated->GetNbinsX();
   binEnd[1] = generated->GetNbinsY();
+  binEnd[2] = generated->GetNbinsZ();
   
   if (fEtaMax > fEtaMin)
   {
@@ -876,6 +1224,11 @@ TH1* AliUEHist::GetTrackEfficiency(CFStep step1, CFStep step2, Int_t axis1, Int_
       binBegin[1] = generated->GetYaxis()->FindBin(fEtaMin);
       binEnd[1]   = generated->GetYaxis()->FindBin(fEtaMax);
     }
+    if (axis3 == 0)
+    {
+      binBegin[2] = generated->GetZaxis()->FindBin(fEtaMin);
+      binEnd[2]   = generated->GetZaxis()->FindBin(fEtaMax);
+    }
   }
   
   if (fPtMax > fPtMin)
@@ -892,13 +1245,18 @@ TH1* AliUEHist::GetTrackEfficiency(CFStep step1, CFStep step2, Int_t axis1, Int_
       binBegin[1] = generated->GetYaxis()->FindBin(fPtMin);
       binEnd[1]   = generated->GetYaxis()->FindBin(ptMax);
     }
+    if (axis3 == 1)
+    {
+      binBegin[2] = generated->GetZaxis()->FindBin(fPtMin);
+      binEnd[2]   = generated->GetZaxis()->FindBin(ptMax);
+    }
   }
   
   Int_t total = 0;
   Int_t count = 0;
-  Int_t vars[2];
+  Int_t vars[3];
   
-  for (Int_t i=0; i<2; i++)
+  for (Int_t i=0; i<3; i++)
     vars[i] = binBegin[i];
     
   const Int_t limit = 50;
@@ -917,8 +1275,23 @@ TH1* AliUEHist::GetTrackEfficiency(CFStep step1, CFStep step2, Int_t axis1, Int_
         generated->GetBinContent(vars[0], vars[1]));
       count++;
     }
+    else if (generated->GetDimension() == 3 && generated->GetBinContent(vars[0], vars[1], vars[2]) < limit)
+    {
+      Printf("Empty bin at %s=%.2f %s=%.2f %s=%.2f (%.2f entries)", 
+        generated->GetXaxis()->GetTitle(), generated->GetXaxis()->GetBinCenter(vars[0]),
+        generated->GetYaxis()->GetTitle(), generated->GetYaxis()->GetBinCenter(vars[1]),
+        generated->GetZaxis()->GetTitle(), generated->GetZaxis()->GetBinCenter(vars[2]),
+        generated->GetBinContent(vars[0], vars[1], vars[2]));
+      count++;
+    }
+    
+    vars[2]++;
+    if (vars[2] == binEnd[2]+1)
+    {
+      vars[2] = binBegin[2];
+      vars[1]++;
+    }
     
-    vars[1]++;
     if (vars[1] == binEnd[1]+1)
     {
       vars[1] = binBegin[1];
@@ -932,6 +1305,94 @@ TH1* AliUEHist::GetTrackEfficiency(CFStep step1, CFStep step2, Int_t axis1, Int_
 
   Printf("Correction has %d empty bins (out of %d bins)", count, total);
   
+  // rebin if required
+  if (source == 2)
+  {
+    TAxis* axis = fEventHist->GetGrid(0)->GetGrid()->GetAxis(0);
+    
+    if (axis->GetNbins() < measured->GetNbinsX())
+    {
+      if (axis2 != -1)
+      {
+        // 2d rebin with variable axis does not exist in root
+        
+        TH1* tmp = measured;
+        measured = new TH2D(Form("%s_rebinned", tmp->GetName()), tmp->GetTitle(), axis->GetNbins(), axis->GetXbins()->GetArray(), tmp->GetNbinsY(), tmp->GetYaxis()->GetXbins()->GetArray());
+        for (Int_t x=1; x<=tmp->GetNbinsX(); x++)
+          for (Int_t y=1; y<=tmp->GetNbinsY(); y++)
+          {
+            ((TH2*) measured)->Fill(tmp->GetXaxis()->GetBinCenter(x), tmp->GetYaxis()->GetBinCenter(y), tmp->GetBinContent(x, y));
+            measured->SetBinError(x, y, 0); // cannot trust bin error, set to 0
+          }
+        delete tmp;
+        
+        tmp = generated;
+        generated = new TH2D(Form("%s_rebinned", tmp->GetName()), tmp->GetTitle(), axis->GetNbins(), axis->GetXbins()->GetArray(), tmp->GetNbinsY(), tmp->GetYaxis()->GetXbins()->GetArray());
+        for (Int_t x=1; x<=tmp->GetNbinsX(); x++)
+          for (Int_t y=1; y<=tmp->GetNbinsY(); y++)
+          {
+            ((TH2*) generated)->Fill(tmp->GetXaxis()->GetBinCenter(x), tmp->GetYaxis()->GetBinCenter(y), tmp->GetBinContent(x, y));
+            generated->SetBinError(x, y, 0); // cannot trust bin error, set to 0
+          }
+        delete tmp;
+      }
+      else
+      {
+        TH1* tmp = measured;
+        measured = tmp->Rebin(axis->GetNbins(), Form("%s_rebinned", tmp->GetName()), axis->GetXbins()->GetArray());
+        delete tmp;
+        
+        tmp = generated;
+        generated = tmp->Rebin(axis->GetNbins(), Form("%s_rebinned", tmp->GetName()), axis->GetXbins()->GetArray());
+        delete tmp;
+      }
+    }
+    else if (axis->GetNbins() > measured->GetNbinsX())
+    {
+      if (axis2 != -1)
+        AliFatal("Rebinning only works for 1d at present");
+  
+      // this is an unfortunate case where the number of bins has to be increased in principle
+      // there is a region where the binning is finner in one histogram and a region where it is the other way round
+      // this cannot be resolved in principle, but as we only calculate the ratio the bin in the second region get the same entries
+      // only a certain binning is supported here
+      if (axis->GetNbins() != 100 || measured->GetNbinsX() != 39)
+        AliFatal(Form("Invalid binning --> %d %d", axis->GetNbins(), measured->GetNbinsX()));
+      
+      Double_t newBins[] = {0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 12.0, 14.0, 16.0, 18.0, 20.0, 25.0, 30.0, 35.0, 40.0, 45.0, 50.0, 100.0};
+      
+      // reduce binning below 5 GeV/c
+      TH1* tmp = measured;
+      measured = tmp->Rebin(27, Form("%s_rebinned", tmp->GetName()), newBins);
+      delete tmp;
+      
+      // increase binning above 5 GeV/c
+      tmp = measured;
+      measured = new TH1F(Form("%s_rebinned2", tmp->GetName()), tmp->GetTitle(), axis->GetNbins(), axis->GetBinLowEdge(1), axis->GetBinUpEdge(axis->GetNbins()));
+      for (Int_t bin=1; bin<=measured->GetNbinsX(); bin++)
+      {
+        measured->SetBinContent(bin, tmp->GetBinContent(tmp->FindBin(measured->GetBinCenter(bin))));
+        measured->SetBinError(bin, tmp->GetBinError(tmp->FindBin(measured->GetBinCenter(bin))));
+      }
+      delete tmp;
+      
+      // reduce binning below 5 GeV/c
+      tmp = generated;
+      generated = tmp->Rebin(27, Form("%s_rebinned", tmp->GetName()), newBins);
+      delete tmp;
+      
+      // increase binning above 5 GeV/c
+      tmp = generated;
+      generated = new TH1F(Form("%s_rebinned2", tmp->GetName()), tmp->GetTitle(), axis->GetNbins(), axis->GetBinLowEdge(1), axis->GetBinUpEdge(axis->GetNbins()));
+      for (Int_t bin=1; bin<=generated->GetNbinsX(); bin++)
+      {
+        generated->SetBinContent(bin, tmp->GetBinContent(tmp->FindBin(generated->GetBinCenter(bin))));
+        generated->SetBinError(bin, tmp->GetBinError(tmp->FindBin(generated->GetBinCenter(bin))));
+      }
+      delete tmp;
+    }
+  }
+  
   measured->Divide(measured, generated, 1, 1, "B");
   
   delete generated;
@@ -959,13 +1420,13 @@ TH1* AliUEHist::GetEventEfficiency(CFStep step1, CFStep step2, Int_t axis1, Int_
     
   if (axis2 >= 0)
   {
-    generated = fEventHist->Project(axis1, axis2, step1);
-    measured = fEventHist->Project(axis1, axis2, step2);
+    generated = fEventHist->Project(step1, axis1, axis2);
+    measured = fEventHist->Project(step2, axis1, axis2);
   }
   else
   {
-    generated = fEventHist->Project(axis1, step1);
-    measured = fEventHist->Project(axis1, step2);
+    generated = fEventHist->Project(step1, axis1);
+    measured = fEventHist->Project(step2, axis1);
   }
   
   measured->Divide(measured, generated, 1, 1, "B");
@@ -1012,7 +1473,7 @@ void AliUEHist::WeightHistogram(TH3* hist1, TH1* hist2)
 }  
 
 //____________________________________________________________________
-TH1* AliUEHist::GetBias(CFStep step1, CFStep step2, Int_t region, const char* axis, Float_t leadPtMin, Float_t leadPtMax)
+TH1* AliUEHist::GetBias(CFStep step1, CFStep step2, Int_t region, const char* axis, Float_t leadPtMin, Float_t leadPtMax, Int_t weighting)
 {
   // extracts the track-level bias (integrating out the multiplicity) between two steps (dividing step2 by step1)
   // in the given region (sum over all regions is calculated if region == -1)
@@ -1020,13 +1481,16 @@ TH1* AliUEHist::GetBias(CFStep step1, CFStep step2, Int_t region, const char* ax
   // and then calculating the ratio between the distributions
   // projected to axis which is a TH3::Project3D string, e.g. "x", or "yx"
   //   no projection is done if axis == 0
+  // weighting: 0 = tracks weighted with events (as discussed above)
+  //            1 = only track bias is returned
+  //            2 = only event bias is returned
   
   AliCFContainer* tmp = 0;
   
   if (region == -1)
   {
     tmp = (AliCFContainer*) fTrackHist[0]->Clone();
-    for (Int_t i = 1; i < fkRegions; i++)
+    for (UInt_t i = 1; i < fkRegions; i++)
       if (fTrackHist[i])
        tmp->Add(fTrackHist[i]);
   }
@@ -1046,43 +1510,55 @@ TH1* AliUEHist::GetBias(CFStep step1, CFStep step2, Int_t region, const char* ax
   ResetBinLimits(fEventHist->GetGrid(step2));
   SetBinLimits(tmp->GetGrid(step2));
   
-  TH1D* events1 = (TH1D*)fEventHist->Project(0, step1);
-  TH3D* hist1 = (TH3D*)tmp->Project(0, tmp->GetNVar()-1, 2, step1);
-  WeightHistogram(hist1, events1);
+  TH1D* events1 = (TH1D*)fEventHist->Project(step1, 0);
+  TH3D* hist1 = (TH3D*)tmp->Project(step1, 0, tmp->GetNVar()-1, 2);
+  if (weighting == 0)
+    WeightHistogram(hist1, events1);
   
-  TH1D* events2 = (TH1D*)fEventHist->Project(0, step2);
-  TH3D* hist2 = (TH3D*)tmp->Project(0, tmp->GetNVar()-1, 2, step2);
-  WeightHistogram(hist2, events2);
+  TH1D* events2 = (TH1D*)fEventHist->Project(step2, 0);
+  TH3D* hist2 = (TH3D*)tmp->Project(step2, 0, tmp->GetNVar()-1, 2);
+  if (weighting == 0)
+    WeightHistogram(hist2, events2);
   
   TH1* generated = hist1;
   TH1* measured = hist2;
   
-  if (axis)
+  if (weighting == 0 || weighting == 1)
   {
-    if (leadPtMax > leadPtMin)
+    if (axis)
     {
-      hist1->GetZaxis()->SetRangeUser(leadPtMin, leadPtMax);
-      hist2->GetZaxis()->SetRangeUser(leadPtMin, leadPtMax);
-    }
+      if (leadPtMax > leadPtMin)
+      {
+        hist1->GetZaxis()->SetRangeUser(leadPtMin, leadPtMax);
+        hist2->GetZaxis()->SetRangeUser(leadPtMin, leadPtMax);
+      }
+      
+      if (fEtaMax > fEtaMin && !TString(axis).Contains("x"))
+      {
+        hist1->GetXaxis()->SetRangeUser(fEtaMin, fEtaMax);
+        hist2->GetXaxis()->SetRangeUser(fEtaMin, fEtaMax);
+      }
     
-    if (fEtaMax > fEtaMin && !TString(axis).Contains("x"))
-    {
-      hist1->GetXaxis()->SetRangeUser(fEtaMin, fEtaMax);
-      hist2->GetXaxis()->SetRangeUser(fEtaMin, fEtaMax);
+      generated = hist1->Project3D(axis);
+      measured  = hist2->Project3D(axis);
+      
+      // delete hists here if projection has been done
+      delete hist1;
+      delete hist2;
     }
-   
-    generated = hist1->Project3D(axis);
-    measured  = hist2->Project3D(axis);
-    
-    // delete hists here if projection has been done
+    delete events1;
+    delete events2;
+  }
+  else if (weighting == 2)
+  {
     delete hist1;
     delete hist2;
+    generated = events1;
+    measured = events2;
   }
   
   measured->Divide(generated);
   
-  delete events1;
-  delete events2;
   delete generated;
   
   ResetBinLimits(tmp->GetGrid(step1));
@@ -1094,6 +1570,87 @@ TH1* AliUEHist::GetBias(CFStep step1, CFStep step2, Int_t region, const char* ax
   return measured;
 }
 
+//____________________________________________________________________
+void AliUEHist::CorrectCorrelatedContamination(CFStep step, Int_t region, TH1* trackCorrection)
+{
+  // corrects for the given factor in a small delta-eta and delta-phi window as function of pT,A and pT,T
+  
+  if (!fTrackHist[region])
+    return;
+   
+  THnSparse* grid = fTrackHist[region]->GetGrid(step)->GetGrid();
+  
+  Int_t var1 = 1;
+  Int_t var2 = 2;
+  
+  if (grid->GetAxis(var1)->GetNbins() != trackCorrection->GetNbinsX())
+    AliFatal(Form("Invalid binning (var1): %d %d", grid->GetAxis(var1)->GetNbins(), trackCorrection->GetNbinsX()));
+    
+  if (grid->GetAxis(var2)->GetNbins() != trackCorrection->GetNbinsY())
+    AliFatal(Form("Invalid binning (var2): %d %d", grid->GetAxis(var2)->GetNbins(), trackCorrection->GetNbinsY()));
+  
+  // optimized implementation
+  for (Int_t binIdx = 0; binIdx < grid->GetNbins(); binIdx++)
+  {
+    Int_t bins[5];
+    
+    Double_t value = grid->GetBinContent(binIdx, bins);
+    Double_t error = grid->GetBinError(binIdx);
+    
+    // check eta and phi axes
+    if (TMath::Abs(grid->GetAxis(0)->GetBinCenter(bins[0])) > 0.1)
+      continue;
+    if (TMath::Abs(grid->GetAxis(4)->GetBinCenter(bins[4])) > 0.1)
+      continue;
+    
+    value *= trackCorrection->GetBinContent(bins[var1], bins[var2]);
+    error *= trackCorrection->GetBinContent(bins[var1], bins[var2]);
+    
+    grid->SetBinContent(bins, value);
+    grid->SetBinError(bins, error);
+  }
+  Printf("AliUEHist::CorrectCorrelatedContamination: Corrected.");
+}
+
+//____________________________________________________________________
+TH2* AliUEHist::GetCorrelatedContamination()
+{
+  // contamination correlated with the trigger particle is evaluated between step kCFStepTracked and kCFStepTrackedOnlyPrim in the region of delta eta and delta phi < 0.1 (smallest bin!)
+  
+  Int_t step1 = kCFStepTrackedOnlyPrim;
+  Int_t step2 = kCFStepTracked;
+  
+  fTrackHist[0]->GetGrid(step1)->SetRangeUser(0, -0.01, 0.01); // delta eta
+  fTrackHist[0]->GetGrid(step1)->SetRangeUser(4, -0.01, 0.01); // delta phi
+  TH2* tracksStep1 = (TH2*) fTrackHist[0]->Project(step1, 1, 2);
+  
+  fTrackHist[0]->GetGrid(step2)->SetRangeUser(0, -0.01, 0.01); // delta eta
+  fTrackHist[0]->GetGrid(step2)->SetRangeUser(4, -0.01, 0.01); // delta phi
+  TH2* tracksStep2 = (TH2*) fTrackHist[0]->Project(step2, 1, 2);
+  
+  tracksStep1->Divide(tracksStep2);
+  
+  TH1* triggersStep1 = fEventHist->Project(step1, 0);
+  TH1* triggersStep2 = fEventHist->Project(step2, 0);
+  
+  TH1* singleParticle = GetTrackingContamination(1);
+  
+  for (Int_t x=1; x<=tracksStep1->GetNbinsX(); x++)
+    for (Int_t y=1; y<=tracksStep1->GetNbinsY(); y++)
+      if (singleParticle->GetBinContent(x) > 0)
+        tracksStep1->SetBinContent(x, y, tracksStep1->GetBinContent(x, y) / triggersStep1->GetBinContent(y) * triggersStep2->GetBinContent(y) / singleParticle->GetBinContent(x));
+      else
+        tracksStep1->SetBinContent(x, y, 0);
+        
+  delete singleParticle;
+  delete tracksStep2;
+  delete triggersStep1;
+  delete triggersStep2;
+        
+  return tracksStep1;
+}
+
 //____________________________________________________________________
 TH2D* AliUEHist::GetTrackingEfficiency()
 {
@@ -1105,6 +1662,17 @@ TH2D* AliUEHist::GetTrackingEfficiency()
   return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepAnaTopology, kCFStepTrackedOnlyPrim, 0, 1));
 }
   
+//____________________________________________________________________
+TH2D* AliUEHist::GetTrackingEfficiencyCentrality()
+{
+  // extracts the tracking efficiency by calculating the efficiency from step kCFStepAnaTopology to kCFStepTrackedOnlyPrim
+  // integrates over the regions and all other variables than pT, centrality to increase the statistics
+  //
+  // returned histogram has to be deleted by the user
+
+  return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepAnaTopology, kCFStepTrackedOnlyPrim, 1, 3));
+}
+
 //____________________________________________________________________
 TH1D* AliUEHist::GetTrackingEfficiency(Int_t axis)
 {
@@ -1145,6 +1713,17 @@ TH2D* AliUEHist::GetTrackingEfficiencyCorrection()
   return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, 0, 1));
 }
   
+//____________________________________________________________________
+TH2D* AliUEHist::GetTrackingEfficiencyCorrectionCentrality()
+{
+  // extracts the tracking correction by calculating the efficiency from step kCFStepAnaTopology to kCFStepTracked
+  // integrates over the regions and all other variables than pT and centrality to increase the statistics
+  //
+  // returned histogram has to be deleted by the user
+
+  return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, 1, 3));
+}
+  
 //____________________________________________________________________
 TH1D* AliUEHist::GetTrackingEfficiencyCorrection(Int_t axis)
 {
@@ -1165,6 +1744,17 @@ TH2D* AliUEHist::GetTrackingContamination()
   return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepTracked, kCFStepTrackedOnlyPrim, 0, 1));
 }
   
+//____________________________________________________________________
+TH2D* AliUEHist::GetTrackingContaminationCentrality()
+{
+  // extracts the tracking contamination by secondaries by calculating the efficiency from step kCFStepTrackedOnlyPrim to kCFStepTracked
+  // integrates over the regions and all other variables than pT and centrality to increase the statistics
+  //
+  // returned histogram has to be deleted by the user
+
+  return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepTracked, kCFStepTrackedOnlyPrim, 1, 3));
+}
+  
 //____________________________________________________________________
 TH1D* AliUEHist::GetTrackingContamination(Int_t axis)
 {
@@ -1257,38 +1847,125 @@ void AliUEHist::ExtendTrackingEfficiency(Bool_t verbose)
   Float_t fitRangeEnd = 14.99;
   Float_t extendRangeBegin = 10.01;
 
-  TH1* obj = GetTrackingEfficiency(1);
-
-  if (verbose)
-    new TCanvas; obj->Draw();
-  obj->Fit("pol0", (verbose) ? "+" : "+0", "SAME", fitRangeBegin, fitRangeEnd);
-
-  Float_t trackingEff = obj->GetFunction("pol0")->GetParameter(0);
-
-  obj = GetTrackingContamination(1);
-
-  if (verbose)
-    new TCanvas; obj->Draw();
-  obj->Fit("pol0", (verbose) ? "+" : "+0", "SAME", fitRangeBegin, fitRangeEnd);
-
-  Float_t trackingCont = obj->GetFunction("pol0")->GetParameter(0);
-
-  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);
-  // extend up to pT 100
-  for (Int_t x = fTrackHistEfficiency->GetAxis(0, 0)->FindBin(fEtaMin); x <= fTrackHistEfficiency->GetAxis(0, 0)->FindBin(fEtaMax); x++)
-    for (Int_t y = fTrackHistEfficiency->GetAxis(1, 0)->FindBin(extendRangeBegin); y <= fTrackHistEfficiency->GetNBins(1); y++)
-      for (Int_t z = 1; z <= fTrackHistEfficiency->GetNBins(2); z++) // particle type axis
+  if (fTrackHistEfficiency->GetNVar() == 3)
+  {
+    TH1* obj = GetTrackingEfficiency(1);
+  
+    if (verbose)
+    {
+      new TCanvas; 
+      obj->Draw();
+    }
+    
+    obj->Fit("pol0", (verbose) ? "+" : "0+", "SAME", fitRangeBegin, fitRangeEnd);
+  
+    Float_t trackingEff = obj->GetFunction("pol0")->GetParameter(0);
+  
+    obj = GetTrackingContamination(1);
+  
+    if (verbose)
+    {
+      new TCanvas; 
+      obj->Draw();
+    }
+    
+    obj->Fit("pol0", (verbose) ? "+" : "0+", "SAME", fitRangeBegin, fitRangeEnd);
+  
+    Float_t trackingCont = obj->GetFunction("pol0")->GetParameter(0);
+  
+    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);
+  
+    // extend for full pT range
+    for (Int_t x = fTrackHistEfficiency->GetAxis(0, 0)->FindBin(fEtaMin); x <= fTrackHistEfficiency->GetAxis(0, 0)->FindBin(fEtaMax); x++)
+      for (Int_t y = fTrackHistEfficiency->GetAxis(1, 0)->FindBin(extendRangeBegin); y <= fTrackHistEfficiency->GetNBins(1); y++)
+        for (Int_t z = 1; z <= fTrackHistEfficiency->GetNBins(2); z++) // particle type axis
+        {
+          
+          Int_t bins[3];
+          bins[0] = x;
+          bins[1] = y;
+          bins[2] = z;
+          
+          fTrackHistEfficiency->GetGrid(0)->SetElement(bins, 100);
+          fTrackHistEfficiency->GetGrid(1)->SetElement(bins, 100.0 * trackingEff);
+          fTrackHistEfficiency->GetGrid(2)->SetElement(bins, 100.0 * trackingEff / trackingCont);
+        }
+  }
+  else if (fTrackHistEfficiency->GetNVar() == 4)
+  {
+    // fit in centrality intervals of 20% for efficiency, one bin for contamination
+    Float_t* trackingEff = 0;
+    Float_t* trackingCont = 0;
+    Int_t centralityBins = 5;
+    
+    Printf("AliUEHist::ExtendTrackingEfficiency: Fitting efficiencies between %f and %f. Extending from %f onwards (within %f < eta < %f)", fitRangeBegin, fitRangeEnd, extendRangeBegin, fEtaMin, fEtaMax);
+    
+    // 0 = eff; 1 = cont
+    for (Int_t caseNo = 0; caseNo < 2; caseNo++)
+    {
+      Float_t* target = 0;
+      Int_t centralityBinsLocal = centralityBins;
+      
+      if (caseNo == 0)
       {
-        Int_t bins[3];
-        bins[0] = x;
-        bins[1] = y;
-        bins[2] = z;
-        
-        fTrackHistEfficiency->GetGrid(0)->SetElement(bins, 100);
-        fTrackHistEfficiency->GetGrid(1)->SetElement(bins, 100.0 * trackingEff);
-        fTrackHistEfficiency->GetGrid(2)->SetElement(bins, 100.0 * trackingEff / trackingCont);
+        trackingEff = new Float_t[centralityBinsLocal];
+        target = trackingEff;
+      }
+      else
+      {
+        centralityBinsLocal = 1;
+        trackingCont = new Float_t[centralityBinsLocal];
+        target = trackingCont;
+      }
+    
+      for (Int_t i=0; i<centralityBinsLocal; i++)
+      {
+        SetCentralityRange(100.0/centralityBinsLocal*i + 0.1, 100.0/centralityBinsLocal*(i+1) - 0.1);
+        TH1* proj = (caseNo == 0) ? GetTrackingEfficiency(1) : GetTrackingContamination(1);
+        if (verbose)
+        {
+          new TCanvas;
+          proj->DrawCopy();
+        }
+        if ((Int_t) proj->Fit("pol0", (verbose) ? "+" : "Q0+", "SAME", fitRangeBegin, fitRangeEnd) == 0)
+          target[i] = proj->GetFunction("pol0")->GetParameter(0);
+        else
+          target[i] = 0;
+        Printf("AliUEHist::ExtendTrackingEfficiency: case %d, centrality %d, eff %f", caseNo, i, target[i]);
       }
+    }
+  
+    // extend for full pT range
+    for (Int_t x = fTrackHistEfficiency->GetAxis(0, 0)->FindBin(fEtaMin); x <= fTrackHistEfficiency->GetAxis(0, 0)->FindBin(fEtaMax); x++)
+      for (Int_t y = fTrackHistEfficiency->GetAxis(1, 0)->FindBin(extendRangeBegin); y <= fTrackHistEfficiency->GetNBins(1); y++)
+        for (Int_t z = 1; z <= fTrackHistEfficiency->GetNBins(2); z++) // particle type axis
+        {
+          for (Int_t z2 = 1; z2 <= fTrackHistEfficiency->GetNBins(3); z2++) // centrality axis
+          {
+            
+            Int_t bins[4];
+            bins[0] = x;
+            bins[1] = y;
+            bins[2] = z;
+            bins[3] = z2;
+            
+            Int_t z2Bin = (Int_t) (fTrackHistEfficiency->GetAxis(3, 0)->GetBinCenter(z2) / (100.0 / centralityBins));
+            //Printf("%d %d", z2, z2Bin);
+            
+            fTrackHistEfficiency->GetGrid(0)->SetElement(bins, 100);
+            fTrackHistEfficiency->GetGrid(1)->SetElement(bins, 100.0 * trackingEff[z2Bin]);
+            if (trackingCont[0] > 0)
+              fTrackHistEfficiency->GetGrid(2)->SetElement(bins, 100.0 * trackingEff[z2Bin] / trackingCont[0]);
+            else  
+              fTrackHistEfficiency->GetGrid(2)->SetElement(bins, 0);
+          }
+       }
+       
+     delete[] trackingEff;
+     delete[] trackingCont;
+   }
+   
+   SetCentralityRange(0, 0);
 }
 
 void AliUEHist::AdditionalDPhiCorrection(Int_t step)
@@ -1322,3 +1999,31 @@ void AliUEHist::AdditionalDPhiCorrection(Int_t step)
     grid->SetBinError(bins, error);
   }
 }
+
+void AliUEHist::Scale(Double_t factor)
+{
+  // scales all contained histograms by the given factor
+  
+  for (Int_t i=0; i<4; i++)
+    if (fTrackHist[i])
+      fTrackHist[i]->Scale(factor);
+  
+  fEventHist->Scale(factor);
+  fTrackHistEfficiency->Scale(factor);
+}
+
+void AliUEHist::Reset()
+{
+  // resets all contained histograms
+  
+  for (Int_t i=0; i<4; i++)
+    if (fTrackHist[i])
+      for (Int_t step=0; step<fTrackHist[i]->GetNStep(); step++)
+        fTrackHist[i]->GetGrid(step)->GetGrid()->Reset();
+  
+  for (Int_t step=0; step<fEventHist->GetNStep(); step++)
+    fEventHist->GetGrid(step)->GetGrid()->Reset();
+    
+  for (Int_t step=0; step<fTrackHistEfficiency->GetNStep(); step++)
+    fTrackHistEfficiency->GetGrid(step)->GetGrid()->Reset();
+}