bug fixes naming

[u/mrichter/AliRoot.git] / PWGCF / Correlations / Base / AliUEHist.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/* $Id: AliUEHist.cxx 20164 2007-08-14 15:31:50Z morsch $ */

//
//
// encapsulate histogram and corrections for one underlying event histogram
//
//
// Author: Jan Fiete Grosse-Oetringhaus, Sara Vallero

#include "AliUEHist.h"
#include "AliCFContainer.h"
#include "THnSparse.h"
#include "TMath.h"
#include "TList.h"
#include "TCollection.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TH3D.h"
#include "AliLog.h"
#include "TCanvas.h"
#include "TF1.h"
#include "AliTHn.h"

ClassImp(AliUEHist)

const Int_t AliUEHist::fgkCFSteps = 10;

AliUEHist::AliUEHist(const char* reqHist) : 
  TObject(),
  fkRegions(4),
  fEventHist(0),
  fTrackHistEfficiency(0),
  fEtaMin(0),
  fEtaMax(0),
  fPtMin(0),
  fPtMax(0),
  fCentralityMin(0),
  fCentralityMax(0),
  fZVtxMin(0),
  fZVtxMax(0),
  fContaminationEnhancement(0),
  fCombineMinMax(0),
  fCache(0),
  fHistogramType(reqHist)
{
  // Constructor
    
  for (UInt_t i=0; i<fkRegions; i++)
    fTrackHist[i] = 0;
    
  if (strlen(reqHist) == 0)
    return;
    
  AliLog::SetClassDebugLevel("AliCFContainer", -1);
  AliLog::SetClassDebugLevel("AliCFGridSparse", -3);
    
  const char* title = "";
    
  // track level
  Int_t nTrackVars = 4; // eta vs pT vs pT,lead (vs delta phi) vs multiplicity
  Int_t iTrackBin[6];
  Double_t* trackBins[6];
  const char* trackAxisTitle[6];
  
  // eta
  iTrackBin[0] = 20;
  Double_t etaBins[20+1];
  for (Int_t i=0; i<=iTrackBin[0]; i++)
    etaBins[i] = -1.0 + 0.1 * i;
  trackBins[0] = etaBins;
  trackAxisTitle[0] = "#eta";
  
  // delta eta
  const Int_t kNDeltaEtaBins = 40;
  Double_t deltaEtaBins[kNDeltaEtaBins+1];
  for (Int_t i=0; i<=kNDeltaEtaBins; i++)
    deltaEtaBins[i] = -2.0 + 0.1 * i;
  
  // pT
  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)";
  
  // pT,lead binning 1
  const Int_t kNLeadingpTBins = 100;
  Double_t leadingpTBins[kNLeadingpTBins+1];
  for (Int_t i=0; i<=kNLeadingpTBins; i++)
    leadingpTBins[i] = 0.5 * i;
  
  // pT,lead binning 2
  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; 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] = -TMath::Pi() / 2 + 1.0 / kNLeadingPhiBins * i * TMath::TwoPi();
    
  // multiplicity
  const Int_t kNMultiplicityBins = 15;
  Double_t multiplicityBins[kNMultiplicityBins+1];
  for (Int_t i=0; i<=kNMultiplicityBins; i++)
    multiplicityBins[i] = -0.5 + i;
  multiplicityBins[kNMultiplicityBins] = 200;
  
  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 };
  
  // vtx-z axis
  const Int_t kNVertexBins = 7;
  Double_t vertexBins[] = { -7, -5, -3, -1, 1, 3, 5, 7 };
  
  Bool_t useVtxAxis = kFALSE;
  
  // selection depending on requested histogram
  Int_t axis = -1; // 0 = pT,lead, 1 = phi,lead
  if (strcmp(reqHist, "NumberDensitypT") == 0)
  {
    axis = 0;
    title = "d^{2}N_{ch}/d#varphid#eta";
  }
  else if (strcmp(reqHist, "NumberDensityPhi") == 0)
  {
    axis = 1;
    title = "d^{2}N_{ch}/d#varphid#eta";
  }
  else if (strcmp(reqHist, "NumberDensityPhiCentrality") == 0 || strcmp(reqHist, "NumberDensityPhiCentralityVtx") == 0)
  {
    if (strcmp(reqHist, "NumberDensityPhiCentralityVtx") == 0)
    {
      reqHist = "NumberDensityPhiCentrality";
      fHistogramType = reqHist;
      useVtxAxis = kTRUE;
    }
    axis = 2;
    title = "d^{2}N_{ch}/d#varphid#eta";
  }
  else if (strcmp(reqHist, "SumpT") == 0)
  {
    axis = 0;
    title = "d^{2}#Sigma p_{T}/d#varphid#eta";
  }
  else
    AliFatal(Form("Invalid histogram requested: %s", reqHist));
  
  UInt_t initRegions = fkRegions;
  
  if (axis == 0)
  {
    trackBins[2] = leadingpTBins;
    iTrackBin[2] = kNLeadingpTBins;
    trackAxisTitle[2] = "leading p_{T} (GeV/c)";
    
  }
  else if (axis == 1)
  {
    nTrackVars = 5;
    initRegions = 1;
  
    iTrackBin[2] = kNLeadingpTBins2;
    trackBins[2] = leadingpTBins2;
    trackAxisTitle[2] = "leading p_{T} (GeV/c)";
    
    iTrackBin[4] = kNLeadingPhiBins;
    trackBins[4] = leadingPhiBins;
    trackAxisTitle[4] = "#Delta #varphi w.r.t. leading track";
  }
  else if (axis == 2)
  {
    nTrackVars = 5;
    initRegions = 1;
  
    iTrackBin[0] = kNDeltaEtaBins;
    trackBins[0] = deltaEtaBins;
    trackAxisTitle[0] = "#Delta#eta";
  
    iTrackBin[2] = kNLeadingpTBins2;
    trackBins[2] = leadingpTBins2;
    trackAxisTitle[2] = "leading p_{T} (GeV/c)";
    
    trackBins[3] = centralityBins;
    iTrackBin[3] = kNCentralityBins;
    trackAxisTitle[3] = "centrality";
  
    iTrackBin[4] = kNLeadingPhiBins;
    trackBins[4] = leadingPhiBins;
    trackAxisTitle[4] = "#Delta#varphi (rad.)";

    if (useVtxAxis)
    {
      nTrackVars = 6;
      iTrackBin[5] = kNVertexBins;
      trackBins[5] = vertexBins;
      trackAxisTitle[5] = "z-vtx (cm)";
    }
  }
    
  for (UInt_t i=0; i<initRegions; i++)
  {
    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++)
    {
      fTrackHist[i]->SetBinLimits(j, trackBins[j]);
      fTrackHist[i]->SetVarTitle(j, trackAxisTitle[j]);
    }
    
    SetStepNames(fTrackHist[i]);
  }
  
  // event level
  Int_t nEventVars = 2;
  Int_t iEventBin[3];

  // track 3rd and 4th axis --> event 1st and 2nd axis
  iEventBin[0] = iTrackBin[2];
  iEventBin[1] = iTrackBin[3];
  
  // plus track 5th axis (in certain cases)
  if (axis == 2 && useVtxAxis)
  {
    nEventVars = 3;
    iEventBin[2] = iTrackBin[5];
  }
  
  fEventHist = new AliCFContainer("fEventHist", title, fgkCFSteps, nEventVars, iEventBin);
  
  fEventHist->SetBinLimits(0, trackBins[2]);
  fEventHist->SetVarTitle(0, trackAxisTitle[2]);
  
  fEventHist->SetBinLimits(1, trackBins[3]);
  fEventHist->SetVarTitle(1, trackAxisTitle[3]);
  
  if (axis == 2 && useVtxAxis)
  {
    fEventHist->SetBinLimits(2, trackBins[5]);
    fEventHist->SetVarTitle(2, trackAxisTitle[5]);
  }

  SetStepNames(fEventHist);
  
  iTrackBin[2] = kNSpeciesBins;

  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]);
}

//_____________________________________________________________________________
AliUEHist::AliUEHist(const AliUEHist &c) :
  TObject(),
  fkRegions(4),
  fEventHist(0),
  fTrackHistEfficiency(0),
  fEtaMin(0),
  fEtaMax(0),
  fPtMin(0),
  fPtMax(0),
  fCentralityMin(0),
  fCentralityMax(0),
  fZVtxMin(0),
  fZVtxMax(0),
  fContaminationEnhancement(0),
  fCombineMinMax(0),
  fCache(0),
  fHistogramType()
{
  //
  // AliUEHist copy constructor
  //

  ((AliUEHist &) c).Copy(*this);
}

//____________________________________________________________________
void AliUEHist::SetStepNames(AliCFContainer* container)
{
  // sets the names of the correction steps
  
  for (Int_t i=0; i<fgkCFSteps; i++)
    container->SetStepTitle(i, GetStepTitle((CFStep) i));
}

//____________________________________________________________________
AliUEHist::~AliUEHist()
{
  // Destructor
  
  for (UInt_t i=0; i<fkRegions; i++)
  {
    if (fTrackHist[i])
    {
      delete fTrackHist[i];
      fTrackHist[i] = 0;
    }
  }
     
  if (fEventHist)
  {
    delete fEventHist;
    fEventHist = 0;
  }
  
  if (fTrackHistEfficiency)
  {
    delete fTrackHistEfficiency;
    fTrackHistEfficiency = 0;
  }

  if (fCache)
  {
    delete fCache;
    fCache = 0;
  }
}

//____________________________________________________________________
AliUEHist &AliUEHist::operator=(const AliUEHist &c)
{
  // assigment operator

  if (this != &c)
    ((AliUEHist &) c).Copy(*this);

  return *this;
}

//____________________________________________________________________
void AliUEHist::Copy(TObject& c) const
{
  // copy function

  AliUEHist& target = (AliUEHist &) c;

  for (UInt_t i=0; i<fkRegions; i++)
    if (fTrackHist[i])
      target.fTrackHist[i] = dynamic_cast<AliCFContainer*> (fTrackHist[i]->Clone());

  if (fEventHist)
    target.fEventHist = dynamic_cast<AliCFContainer*> (fEventHist->Clone());
  
  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;
}

//____________________________________________________________________
Long64_t AliUEHist::Merge(TCollection* list)
{
  // Merge a list of AliUEHist objects with this (needed for
  // PROOF). 
  // Returns the number of merged objects (including this).

  if (!list)
    return 0;
  
  if (list->IsEmpty())
    return 1;

  TIterator* iter = list->MakeIterator();
  TObject* obj;

  // collections of objects
  const UInt_t kMaxLists = fkRegions+2;
  TList** lists = new TList*[kMaxLists];
  
  for (UInt_t i=0; i<kMaxLists; i++)
    lists[i] = new TList;
  
  Int_t count = 0;
  while ((obj = iter->Next())) {
    
    AliUEHist* entry = dynamic_cast<AliUEHist*> (obj);
    if (entry == 0) 
      continue;

    for (UInt_t i=0; i<fkRegions; i++)
      if (entry->fTrackHist[i])
        lists[i]->Add(entry->fTrackHist[i]);
    
    lists[fkRegions]->Add(entry->fEventHist);
    lists[fkRegions+1]->Add(entry->fTrackHistEfficiency);

    count++;
  }
  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 (UInt_t i=0; i<kMaxLists; i++)
    delete lists[i];
    
  delete[] lists;

  return count+1;
}

//____________________________________________________________________
void AliUEHist::SetBinLimits(AliCFGridSparse* grid)
{
  // sets the bin limits in eta and pT defined by fEtaMin/Max, fPtMin/Max
  
  if (fEtaMax > fEtaMin)
    grid->SetRangeUser(0, fEtaMin, fEtaMax);
  if (fPtMax > fPtMin)
    grid->SetRangeUser(1, fPtMin, fPtMax);
}  

//____________________________________________________________________
void AliUEHist::ResetBinLimits(AliCFGridSparse* grid)
{
  // resets all bin limits 
  
  for (Int_t i=0; i<grid->GetNVar(); i++)
    if (grid->GetGrid()->GetAxis(i)->TestBit(TAxis::kAxisRange))
      grid->SetRangeUser(i, 0, -1);
}
  
//____________________________________________________________________
void AliUEHist::CountEmptyBins(AliUEHist::CFStep step, Float_t ptLeadMin, Float_t ptLeadMax)
{
  // prints the number of empty bins in the track end event histograms in the given step
  
  Int_t binBegin[4];
  Int_t binEnd[4];
  
  for (Int_t i=0; i<4; i++)
  {
    binBegin[i] = 1;
    binEnd[i]   = fTrackHist[0]->GetNBins(i);
  }
  
  if (fEtaMax > fEtaMin)
  {
    binBegin[0] = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(0)->FindBin(fEtaMin);
    binEnd[0]   = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(0)->FindBin(fEtaMax);
  }
  
  if (fPtMax > fPtMin)
  {
    binBegin[1] = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(1)->FindBin(fPtMin);
    binEnd[1]   = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(1)->FindBin(fPtMax);
  }
  
  if (ptLeadMax > ptLeadMin)
  {
    binBegin[2] = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(2)->FindBin(ptLeadMin);
    binEnd[2]   = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(2)->FindBin(ptLeadMax);
  }
  
  // start from multiplicity 1
  binBegin[3] = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(3)->FindBin(1);
  
  for (UInt_t region=0; region<fkRegions; region++)
  {
    Int_t total = 0;
    Int_t count = 0;
    Int_t vars[4];
    
    for (Int_t i=0; i<4; i++)
      vars[i] = binBegin[i];
      
    AliCFGridSparse* grid = fTrackHist[region]->GetGrid(step);
    while (1)
    {
      if (grid->GetElement(vars) == 0)
      {
        Printf("Empty bin at eta=%.2f pt=%.2f pt_lead=%.2f mult=%.1f", 
          grid->GetBinCenter(0, vars[0]), 
          grid->GetBinCenter(1, vars[1]), 
          grid->GetBinCenter(2, vars[2]), 
          grid->GetBinCenter(3, vars[3])
        );
        count++;
      }
      
      vars[3]++;
      for (Int_t i=3; i>0; i--)
      {
        if (vars[i] == binEnd[i]+1)
        {
          vars[i] = binBegin[i];
          vars[i-1]++;
        }
      }
      
      if (vars[0] == binEnd[0]+1)
        break;
      total++;
    }
  
    Printf("Region %s has %d empty bins (out of %d bins)", GetRegionTitle((Region) region), count, total);
  }
}  

//____________________________________________________________________
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, Int_t* normEvents)
{
  // 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
  //
  // normEvents: if non-0 the number of events/trigger particles for the normalization is filled
  
  // 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;
  
  if (ptLeadMin < 0)
  {
    tracks = fTrackHist[region]->ShowProjection(2, step);
    tracks->GetYaxis()->SetTitle(fTrackHist[region]->GetTitle());
    if (fCombineMinMax && region == kMin)
    {
      ResetBinLimits(fTrackHist[kMax]->GetGrid(step));
      SetBinLimits(fTrackHist[kMax]->GetGrid(step));
      
      TH1D* tracks2 = fTrackHist[kMax]->ShowProjection(2, step);
      tracks->Add(tracks2);
      
      ResetBinLimits(fTrackHist[kMax]->GetGrid(step));
    }
      
    // normalize to get a density (deta dphi)
    TAxis* axis = fTrackHist[region]->GetGrid(step)->GetAxis(0);
    Float_t phiRegion = TMath::TwoPi() / 3;
    if (!fCombineMinMax && region == kMin)
      phiRegion /= 2;
    Float_t normalization = phiRegion;
    if (etaNorm)
      normalization *= (axis->GetBinUpEdge(axis->GetLast()) - axis->GetBinLowEdge(axis->GetFirst()));
    //Printf("Normalization: %f", normalization);
    tracks->Scale(1.0 / normalization);
    
    TH1D* events = fEventHist->ShowProjection(0, step);
    tracks->Divide(events);
  }
  else
  {
    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);
    }
    
    Int_t firstBin = fTrackHist[region]->GetAxis(2, step)->FindBin(ptLeadMin);
    Int_t lastBin = fTrackHist[region]->GetAxis(2, step)->FindBin(ptLeadMax);
    
    Printf("Using leading pT range %d --> %d", firstBin, lastBin);
    
    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);
      
    Printf("Calculated histogram --> %f tracks", tracks->Integral());
    fTrackHist[region]->GetGrid(step)->SetRangeUser(2, 0, -1);
    
    if (twoD == 10 || twoD == 11)
    {
      Float_t etaLimit = 1.0;
      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);
      if (strcmp(axis->GetTitle(), "#eta") == 0)
      {
	Printf("Normalizing using eta-axis range");
	normalization *= axis->GetBinUpEdge(axis->GetLast()) - axis->GetBinLowEdge(axis->GetFirst());
      }
      else
      {
	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 (normEvents)
      *normEvents = nEvents;
      
    if (nEvents > 0)
      tracks->Scale(1.0 / nEvents);
    
    delete events;
  }

  ResetBinLimits(fTrackHist[region]->GetGrid(step));
  ResetBinLimits(fEventHist->GetGrid(step));

  return tracks;
}

//____________________________________________________________________
void AliUEHist::GetHistsZVtx(AliUEHist::CFStep step, AliUEHist::Region region, Float_t ptLeadMin, Float_t ptLeadMax, Int_t multBinBegin, Int_t multBinEnd, TH3** trackHist, TH1** eventHist)
{
  // Calculates a 3d histogram with deltaphi, deltaeta, zvtx on track level and a 1d histogram on event level (as fct of zvtx)
  // Histogram has to be deleted by the caller of the function
  
  if (!fTrackHist[kToward]->GetGrid(0)->GetGrid()->GetAxis(5))
    AliFatal("Histogram without vertex axis provided");

  // unzoom all axes
  ResetBinLimits(fTrackHist[region]->GetGrid(step));
  ResetBinLimits(fEventHist->GetGrid(step));
  
  SetBinLimits(fTrackHist[region]->GetGrid(step));
  
  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);
  }
    
  Int_t firstBin = fTrackHist[region]->GetAxis(2, step)->FindBin(ptLeadMin);
  Int_t lastBin = fTrackHist[region]->GetAxis(2, step)->FindBin(ptLeadMax);
  Printf("Using leading pT range %d --> %d", firstBin, lastBin);
  fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(2)->SetRange(firstBin, lastBin);
  fEventHist->GetGrid(step)->GetGrid()->GetAxis(0)->SetRange(firstBin, lastBin);
    
  *trackHist = (TH3*) fTrackHist[region]->GetGrid(step)->Project(4, 0, 5);
  *eventHist = (TH1*) fEventHist->GetGrid(step)->Project(2);

  ResetBinLimits(fTrackHist[region]->GetGrid(step));
  ResetBinLimits(fEventHist->GetGrid(step));
}

//____________________________________________________________________
TH2* AliUEHist::GetSumOfRatios2(AliUEHist* mixed, AliUEHist::CFStep step, AliUEHist::Region region, Float_t ptLeadMin, Float_t ptLeadMax, Int_t multBinBegin, Int_t multBinEnd)
{
  // Calls GetUEHist(...) for *each* vertex bin and performs a sum of ratios:
  // 1_N [ (same/mixed)_1 + (same/mixed)_2 + (same/mixed)_3 + ... ]
  // where N is the total number of events/trigger particles and the subscript is the vertex bin
  // where mixed is normalized such that the information about the number of pairs in same is kept
  //
  // returns a 2D histogram: deltaphi, deltaeta
  //
  // Parameters:
  //   mixed: AliUEHist containing mixed event corresponding to this object
  //   <other parameters> : check documentation of AliUEHist::GetUEHist
  
  TH2* totalTracks = 0;
  
  TH3* trackSameAll = 0;
  TH3* trackMixedAll = 0;
  TH1* eventSameAll = 0;
  TH1* eventMixedAll = 0;
  
  GetHistsZVtx(step, region, ptLeadMin, ptLeadMax, multBinBegin, multBinEnd, &trackSameAll, &eventSameAll);
  mixed->GetHistsZVtx(step, region, ptLeadMin, ptLeadMax, multBinBegin, multBinEnd, &trackMixedAll, &eventMixedAll);
  
  TAxis* vertexAxis = trackSameAll->GetZaxis();
  for (Int_t vertexBin = 1; vertexBin <= vertexAxis->GetNbins(); vertexBin++)
  {
    trackSameAll->GetZaxis()->SetRange(vertexBin, vertexBin);
    trackMixedAll->GetZaxis()->SetRange(vertexBin, vertexBin);

    TH2* tracksSame = (TH2*) trackSameAll->Project3D("yx1");
    TH2* tracksMixed = (TH2*) trackMixedAll->Project3D("yx2");
    
    // asssume flat in dphi, gain in statistics
//     TH1* histMixedproj = mixedTwoD->ProjectionY();
//     histMixedproj->Scale(1.0 / mixedTwoD->GetNbinsX());
//     
//     for (Int_t x=1; x<=mixedTwoD->GetNbinsX(); x++)
//       for (Int_t y=1; y<=mixedTwoD->GetNbinsY(); y++)
// 	mixedTwoD->SetBinContent(x, y, histMixedproj->GetBinContent(y));

//       delete histMixedproj;

    // get mixed event normalization by assuming full acceptance at deta of 0
    Double_t mixedNorm = tracksMixed->Integral(tracksMixed->GetXaxis()->FindBin(-0.01), tracksMixed->GetXaxis()->FindBin(0.01), tracksMixed->GetYaxis()->FindBin(-0.01), tracksMixed->GetYaxis()->FindBin(0.01));
    mixedNorm /= (tracksMixed->GetXaxis()->FindBin(0.01) - tracksMixed->GetXaxis()->FindBin(-0.01) + 1) * (tracksMixed->GetYaxis()->FindBin(0.01) - tracksMixed->GetYaxis()->FindBin(-0.01) + 1);
    tracksMixed->Scale(1.0 / mixedNorm);

//     tracksSame->Scale(tracksMixed->Integral() / tracksSame->Integral());
    
    tracksSame->Divide(tracksMixed);
    
    // code to draw contributions
    /*
    TH1* proj = tracksSame->ProjectionX("projx", tracksSame->GetYaxis()->FindBin(-1.59), tracksSame->GetYaxis()->FindBin(1.59));
    proj->SetTitle(Form("Bin %d", vertexBin));
    proj->SetLineColor(vertexBin);
    proj->DrawCopy((vertexBin > 1) ? "SAME" : "");
    */
    
    if (!totalTracks)
      totalTracks = (TH2*) tracksSame->Clone("totalTracks");
    else
      totalTracks->Add(tracksSame);

    delete tracksSame;
    delete tracksMixed;
  }

  Int_t totalEvents = eventSameAll->Integral();
  if (totalTracks) {
    Printf("Dividing %f tracks by %d events", totalTracks->Integral(), totalEvents);
    if (totalEvents > 0)
      totalTracks->Scale(1.0 / totalEvents);
  }
  
  // normalizate to dphi width
  Float_t normalization = totalTracks->GetXaxis()->GetBinWidth(1);
  totalTracks->Scale(1.0 / normalization);
  
  delete trackSameAll;
  delete trackMixedAll;
  delete eventSameAll;
  delete eventMixedAll;
  
  return totalTracks;
}

//____________________________________________________________________
TH2* AliUEHist::GetSumOfRatios(AliUEHist* mixed, AliUEHist::CFStep step, AliUEHist::Region region, Float_t ptLeadMin, Float_t ptLeadMax, Int_t multBinBegin, Int_t multBinEnd, Bool_t etaNorm, Bool_t useVertexBins)
{
  // Calls GetUEHist(...) for *each* multiplicity bin and performs a sum of ratios:
  // 1_N [ (same/mixed)_1 + (same/mixed)_2 + (same/mixed)_3 + ... ]
  // where N is the total number of events/trigger particles and the subscript is the multiplicity bin
  //
  // Can only be used for the 2D histogram at present
  //
  // Parameters:
  //   mixed: AliUEHist containing mixed event corresponding to this object
  //   <other parameters> : check documentation of AliUEHist::GetUEHist
  
  TH2* totalTracks = 0;
  Int_t totalEvents = 0;
  
  Int_t vertexBin = 1;
  TAxis* vertexAxis = fTrackHist[kToward]->GetGrid(0)->GetGrid()->GetAxis(5);
  if (useVertexBins && !vertexAxis)
  {
    Printf("Vertex axis requested but not available");
    return 0;
  }
  
  // vertex bin loop
  while (1)
  {
    if (useVertexBins)
    {
      SetZVtxRange(vertexAxis->GetBinLowEdge(vertexBin) + 0.01, vertexAxis->GetBinUpEdge(vertexBin) - 0.01);
      mixed->SetZVtxRange(vertexAxis->GetBinLowEdge(vertexBin) + 0.01, vertexAxis->GetBinUpEdge(vertexBin) - 0.01);
      vertexBin++;
    }
    
    // multiplicity loop
    Int_t multIter = multBinBegin;
    while (1)
    {
      Int_t multBinBeginLocal = multBinBegin;
      Int_t multBinEndLocal = multBinEnd;
      
      if (multBinEnd >= multBinBegin)
      {
	multBinBeginLocal = multIter;
	multBinEndLocal = multIter;
	multIter++;
      }
	
      Int_t nEvents = 0;
      TH2* tracks = (TH2*) GetUEHist(step, region, ptLeadMin, ptLeadMax, multBinBeginLocal, multBinEndLocal, 1, etaNorm, &nEvents);
      // undo normalization
      tracks->Scale(nEvents);
      totalEvents += nEvents;
      
      TH2* mixedTwoD = (TH2*) mixed->GetUEHist(step, region, ptLeadMin, ptLeadMax, multBinBeginLocal, multBinEndLocal, 1, etaNorm);
      
      // asssume flat in dphi, gain in statistics
  //     TH1* histMixedproj = mixedTwoD->ProjectionY();
  //     histMixedproj->Scale(1.0 / mixedTwoD->GetNbinsX());
  //     
  //     for (Int_t x=1; x<=mixedTwoD->GetNbinsX(); x++)
  //       for (Int_t y=1; y<=mixedTwoD->GetNbinsY(); y++)
  // 	mixedTwoD->SetBinContent(x, y, histMixedproj->GetBinContent(y));

  //       delete histMixedproj;
  
      // get mixed event normalization by assuming full acceptance at deta of 0 (only works for flat dphi)
/*      Double_t mixedNorm = mixedTwoD->Integral(1, mixedTwoD->GetNbinsX(), mixedTwoD->GetYaxis()->FindBin(-0.01), mixedTwoD->GetYaxis()->FindBin(0.01));
      mixedNorm /= mixedTwoD->GetNbinsX() * (mixedTwoD->GetYaxis()->FindBin(0.01) - mixedTwoD->GetYaxis()->FindBin(-0.01) + 1);
      tracks->Scale(mixedNorm);*/
      
      tracks->Scale(mixedTwoD->Integral() / tracks->Integral());

      tracks->Divide(mixedTwoD);
      
      delete mixedTwoD;
      
      if (!totalTracks)
	totalTracks = tracks;
      else
      {
	totalTracks->Add(tracks);
	delete tracks;
      }

      if (multIter > multBinEnd)
	break;
    }
    
    if (!useVertexBins || vertexBin > vertexAxis->GetNbins())
      break;
  }

  if (useVertexBins)
    totalEvents = vertexAxis->GetNbins();
  if (totalTracks) {
    Printf("Dividing %f tracks by %d events", totalTracks->Integral(), totalEvents);
    if (totalEvents > 0)
      totalTracks->Scale(1.0 / totalEvents);
  }
  
  return totalTracks;
}

//____________________________________________________________________
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::MultiplyHistograms(THnSparse* grid, THnSparse* target, TH1* histogram, Int_t var1, Int_t var2)
{
  // multiplies <grid> with <histogram> and put the result in <target>
  // <grid> has usually more dimensions than histogram. The axis which are used to choose the value 
  // from <histogram> are given in <var1> and <var2>
  //
  // if <histogram> is 0, just copies content from step1 to step2
  
  // clear target histogram
  target->Reset();
  
  if (histogram != 0)
  {
    if (grid->GetAxis(var1)->GetNbins() != histogram->GetNbinsX())
      AliFatal(Form("Invalid binning (var1): %d %d", grid->GetAxis(var1)->GetNbins(), histogram->GetNbinsX()));
      
    if (var2 >= 0 && grid->GetAxis(var2)->GetNbins() != histogram->GetNbinsY())
      AliFatal(Form("Invalid binning (var2): %d %d", grid->GetAxis(var2)->GetNbins(), histogram->GetNbinsY()));
  }

  if (grid->GetNdimensions() > 6)
    AliFatal("Too many dimensions in THnSparse");
  
  Int_t bins[6];
  
  // optimized implementation
  for (Int_t binIdx = 0; binIdx < grid->GetNbins(); binIdx++)
  {
    Double_t value = grid->GetBinContent(binIdx, bins);
    Double_t error = grid->GetBinError(binIdx);
    
    if (histogram != 0)
    {
      if (var2 < 0)
      {
        value *= histogram->GetBinContent(bins[var1]);
        error *= histogram->GetBinContent(bins[var1]);
      }
      else
      {
        value *= histogram->GetBinContent(bins[var1], bins[var2]);
        error *= histogram->GetBinContent(bins[var1], bins[var2]);
      }
    }
    
    target->SetBinContent(bins, value);
    target->SetBinError(bins, error);
  }
}

//____________________________________________________________________
void AliUEHist::CorrectTracks(CFStep step1, CFStep step2, TH1* trackCorrection, Int_t var1, Int_t var2)
{
  // corrects from step1 to step2 by multiplying the tracks with trackCorrection
  // trackCorrection can be a function of eta (var1 == 0), pT (var1 == 1), leading pT (var1 == 2), multiplicity (var1 == 3), delta phi (var1 == 4)
  // if var2 >= 0 a two dimension correction is assumed in trackCorrection
  //
  // if trackCorrection is 0, just copies content from step1 to step2
  
  for (UInt_t region=0; region<fkRegions; region++)
    CorrectTracks(step1, step2, region, trackCorrection, var1, var2);
}

//____________________________________________________________________
void AliUEHist::CorrectTracks(CFStep step1, CFStep step2, Int_t region, TH1* trackCorrection, Int_t var1, Int_t var2)
{
  //
  // see documentation of CorrectTracks above
  //
  
  if (!fTrackHist[region])
    return;
   
  THnSparse* grid = fTrackHist[region]->GetGrid(step1)->GetGrid();
  THnSparse* target = fTrackHist[region]->GetGrid(step2)->GetGrid();
  
  MultiplyHistograms(grid, target, trackCorrection, var1, var2);
  
  Printf("AliUEHist::CorrectTracks: Corrected from %f to %f entries. Correction histogram: %f entries (integral: %f)", grid->GetEntries(), target->GetEntries(), (trackCorrection) ? trackCorrection->GetEntries() : -1.0, (trackCorrection) ? trackCorrection->Integral() : -1.0); 
}

//____________________________________________________________________
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)
  //
  // if eventCorrection is 0, just copies content from step1 to step2
  
  AliCFGridSparse* grid = fEventHist->GetGrid(step1);
  AliCFGridSparse* target = fEventHist->GetGrid(step2);
  
  MultiplyHistograms(grid->GetGrid(), target->GetGrid(), eventCorrection, var1, var2);

  Printf("AliUEHist::CorrectEvents: Corrected from %f to %f entries. Correction histogram: %f entries (integral: %f)", grid->GetEntries(), target->GetEntries(), (eventCorrection) ? eventCorrection->GetEntries() : -1.0, (eventCorrection) ? eventCorrection->Integral() : -1.0); 
}

//____________________________________________________________________
void AliUEHist::Correct(AliUEHist* corrections)
{
  // applies the given corrections to extract from the step kCFStepReconstructed all previous steps
  //
  // in this object the data is expected in the step kCFStepReconstructed
  
  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#varphid#eta") == 0)
        fHistogramType = "NumberDensitypT";
      else if (strcmp(fTrackHist[kToward]->GetTitle(), "d^{2}#Sigma p_{T}/d#varphid#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());
  
  if (strcmp(fHistogramType, "NumberDensitypT") == 0 || strcmp(fHistogramType, "NumberDensityPhi") == 0 || strcmp(fHistogramType, "SumpT") == 0)
  {
    // ---- 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 (strcmp(fHistogramType, "NumberDensityPhi") == 0)
    {
      projAxis = "yz";
      secondBin = 4;
    }
    
    for (UInt_t region = 0; region < fkRegions; region++)
    {
      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));
    }
  
    #if 0
    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 if (strcmp(fHistogramType, "NumberDensityPhiCentrality") == 0)
  {
    // 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);
//     Printf("\n\n\nWARNING ---> SKIPPING CorrectCorrelatedContamination\n\n\n");
    
    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);
  }
  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, 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; 2 = fTrackHistEfficiency rebinned for pT,T / pT,lead binning
  
  // integrate over regions
  // cache it for efficiency (usually more than one efficiency is requested)
  
  AliCFContainer* sourceContainer = 0;
  
  if (source == 0)
  {
    if (!fCache)
    {
      fCache = (AliCFContainer*) fTrackHist[0]->Clone();
      for (UInt_t i = 1; i < fkRegions; i++)
        if (fTrackHist[i])
          fCache->Add(fTrackHist[i]);
    }
    sourceContainer = fCache;
  }
  else if (source == 1 || source == 2)
  {
    sourceContainer = fTrackHistEfficiency;
    // step offset because we start with kCFStepAnaTopology
    step1 = (CFStep) ((Int_t) step1 - (Int_t) kCFStepAnaTopology);
    step2 = (CFStep) ((Int_t) step2 - (Int_t) kCFStepAnaTopology);
  }
  else
    return 0;
        
  // 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 && axis3 != 0)
  {
    Printf("Restricted eta-range to %f %f", fEtaMin, fEtaMax);
    sourceContainer->GetGrid(step1)->SetRangeUser(0, fEtaMin, fEtaMax);
    sourceContainer->GetGrid(step2)->SetRangeUser(0, fEtaMin, fEtaMax);
  }
  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 (axis3 >= 0)
  {
    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(step1, axis1);
    measured = sourceContainer->Project(step2, axis1);
  }
  
  // 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)
  {
    if (axis1 == 0)
    {
      binBegin[0] = generated->GetXaxis()->FindBin(fEtaMin);
      binEnd[0]   = generated->GetXaxis()->FindBin(fEtaMax);
    }
    if (axis2 == 0)
    {
      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)
  {
    // TODO this is just checking up to 15 for now
    Float_t ptMax = TMath::Min((Float_t) 15., fPtMax);
    if (axis1 == 1)
    {
      binBegin[0] = generated->GetXaxis()->FindBin(fPtMin);
      binEnd[0]   = generated->GetXaxis()->FindBin(ptMax);
    }
    if (axis2 == 1)
    {
      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[3];
  
  for (Int_t i=0; i<3; i++)
    vars[i] = binBegin[i];
    
  const Int_t limit = 50;
  while (1)
  {
    if (generated->GetDimension() == 1 && generated->GetBinContent(vars[0]) < limit)
    {
      Printf("Empty bin at %s=%.2f (%.2f entries)", generated->GetXaxis()->GetTitle(), generated->GetXaxis()->GetBinCenter(vars[0]), generated->GetBinContent(vars[0]));
      count++;
    } 
    else if (generated->GetDimension() == 2 && generated->GetBinContent(vars[0], vars[1]) < limit)
    {
      Printf("Empty bin at %s=%.2f %s=%.2f (%.2f entries)", 
        generated->GetXaxis()->GetTitle(), generated->GetXaxis()->GetBinCenter(vars[0]),
        generated->GetYaxis()->GetTitle(), generated->GetYaxis()->GetBinCenter(vars[1]),
        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]++;
    }
    
    if (vars[1] == binEnd[1]+1)
    {
      vars[1] = binBegin[1];
      vars[0]++;
    }
    
    if (vars[0] == binEnd[0]+1)
      break;
    total++;
  }

  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;
  
  ResetBinLimits(sourceContainer->GetGrid(step1));
  ResetBinLimits(sourceContainer->GetGrid(step2));
  
  return measured;
}

//____________________________________________________________________
TH1* AliUEHist::GetEventEfficiency(CFStep step1, CFStep step2, Int_t axis1, Int_t axis2, Float_t ptLeadMin, Float_t ptLeadMax)
{
  // creates a event-level efficiency by dividing step2 by step1
  // projected to axis1 and axis2 (optional if >= 0)
  
  if (ptLeadMax > ptLeadMin)
  {
    fEventHist->GetGrid(step1)->SetRangeUser(0, ptLeadMin, ptLeadMax);
    fEventHist->GetGrid(step2)->SetRangeUser(0, ptLeadMin, ptLeadMax);
  }
  
  TH1* measured = 0;
  TH1* generated = 0;
    
  if (axis2 >= 0)
  {
    generated = fEventHist->Project(step1, axis1, axis2);
    measured = fEventHist->Project(step2, axis1, axis2);
  }
  else
  {
    generated = fEventHist->Project(step1, axis1);
    measured = fEventHist->Project(step2, axis1);
  }
  
  measured->Divide(measured, generated, 1, 1, "B");
  
  delete generated;
  
  if (ptLeadMax > ptLeadMin)
  {
    fEventHist->GetGrid(step1)->SetRangeUser(0, 0, -1);
    fEventHist->GetGrid(step2)->SetRangeUser(0, 0, -1);
  }
  
  return measured;
}

//____________________________________________________________________
void AliUEHist::WeightHistogram(TH3* hist1, TH1* hist2)
{
  // weights each entry of the 3d histogram hist1 with the 1d histogram hist2 
  // where the matching is done of the z axis of hist1 with the x axis of hist2
  
  if (hist1->GetNbinsZ() != hist2->GetNbinsX())
    AliFatal(Form("Inconsistent binning %d %d", hist1->GetNbinsZ(), hist2->GetNbinsX()));
  
  for (Int_t x=1; x<=hist1->GetNbinsX(); x++)
  {
    for (Int_t y=1; y<=hist1->GetNbinsY(); y++)
    {
      for (Int_t z=1; z<=hist1->GetNbinsZ(); z++)
      {
        if (hist2->GetBinContent(z) > 0)
        {
          hist1->SetBinContent(x, y, z, hist1->GetBinContent(x, y, z) / hist2->GetBinContent(z));
          hist1->SetBinError(x, y, z, hist1->GetBinError(x, y, z) / hist2->GetBinContent(z));
        }
        else
        {
          hist1->SetBinContent(x, y, z, 0);
          hist1->SetBinError(x, y, z, 0);
        }
      }
    }
  }
}  

//____________________________________________________________________
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)
  // done by weighting the track-level distribution with the number of events as function of leading pT
  // 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 (UInt_t i = 1; i < fkRegions; i++)
      if (fTrackHist[i])
	tmp->Add(fTrackHist[i]);
  }
  else if (region == kMin && fCombineMinMax)
  {
    tmp = (AliCFContainer*) fTrackHist[kMin]->Clone();
    tmp->Add(fTrackHist[kMax]);
  }
  else
    tmp = fTrackHist[region];
  
  ResetBinLimits(tmp->GetGrid(step1));
  ResetBinLimits(fEventHist->GetGrid(step1));
  SetBinLimits(tmp->GetGrid(step1));
  
  ResetBinLimits(tmp->GetGrid(step2));
  ResetBinLimits(fEventHist->GetGrid(step2));
  SetBinLimits(tmp->GetGrid(step2));
  
  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(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 (weighting == 0 || weighting == 1)
  {
    if (axis)
    {
      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);
      }
    
      generated = hist1->Project3D(axis);
      measured  = hist2->Project3D(axis);
      
      // delete hists here if projection has been done
      delete hist1;
      delete hist2;
    }
    delete events1;
    delete events2;
  }
  else if (weighting == 2)
  {
    delete hist1;
    delete hist2;
    generated = events1;
    measured = events2;
  }
  
  measured->Divide(generated);
  
  delete generated;
  
  ResetBinLimits(tmp->GetGrid(step1));
  ResetBinLimits(tmp->GetGrid(step2));
  
  if ((region == -1) || (region == kMin && fCombineMinMax))
    delete tmp;
  
  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[6];
    
    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 && triggersStep1->GetBinContent(y) > 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()
{
  // extracts the tracking efficiency by calculating the efficiency from step kCFStepAnaTopology to kCFStepTrackedOnlyPrim
  // integrates over the regions and all other variables than pT and eta to increase the statistics
  //
  // returned histogram has to be deleted by the user

  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)
{
  // extracts the tracking efficiency by calculating the efficiency from step kCFStepAnaTopology to kCFStepTrackedOnlyPrim
  // integrates over the regions and all other variables than pT (axis == 0) and eta (axis == 1) to increase the statistics

  return dynamic_cast<TH1D*> (GetTrackEfficiency(kCFStepAnaTopology, kCFStepTrackedOnlyPrim, axis));
}

//____________________________________________________________________
TH2D* AliUEHist::GetTrackingCorrection()
{
  // extracts the tracking correction by calculating the efficiency from step kCFStepAnaTopology to kCFStepTracked
  // integrates over the regions and all other variables than pT and eta to increase the statistics
  //
  // returned histogram has to be deleted by the user

  return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepTracked, kCFStepAnaTopology, 0, 1));
}
  
//____________________________________________________________________
TH1D* AliUEHist::GetTrackingCorrection(Int_t axis)
{
  // extracts the tracking correction by calculating the efficiency from step kCFStepAnaTopology to kCFStepTracked
  // integrates over the regions and all other variables than pT (axis == 0) and eta (axis == 1) to increase the statistics

  return dynamic_cast<TH1D*> (GetTrackEfficiency(kCFStepTracked, kCFStepAnaTopology, axis));
}

//____________________________________________________________________
TH2D* AliUEHist::GetTrackingEfficiencyCorrection()
{
  // extracts the tracking correction by calculating the efficiency from step kCFStepAnaTopology to kCFStepTracked
  // integrates over the regions and all other variables than pT and eta to increase the statistics
  //
  // returned histogram has to be deleted by the user

  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)
{
  // extracts the tracking correction by calculating the efficiency from step kCFStepAnaTopology to kCFStepTracked
  // integrates over the regions and all other variables than pT (axis == 0) and eta (axis == 1) to increase the statistics

  return dynamic_cast<TH1D*> (GetTrackEfficiency(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, axis));
}

//____________________________________________________________________
TH2D* AliUEHist::GetTrackingContamination()
{
  // 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 eta to increase the statistics
  //
  // returned histogram has to be deleted by the user

  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)
{
  // 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 (axis == 0) and eta (axis == 1) to increase the statistics

  return dynamic_cast<TH1D*> (GetTrackEfficiency(kCFStepTracked, kCFStepTrackedOnlyPrim, axis));
}

//____________________________________________________________________
const char* AliUEHist::GetRegionTitle(Region region)
{
  // returns the name of the given region
  
  switch (region)
  {
    case kToward:
      return "Towards";
    case kAway:
      return "Away";
    case kMin:
      return (fCombineMinMax) ? "Transverse" : "Min";
    case kMax:
      return "Max";
  }
  
  return 0;
}
  
//____________________________________________________________________
const char* AliUEHist::GetStepTitle(CFStep step)
{
  // returns the name of the given step
  
  switch (step)
  {
    case kCFStepAll:
      return "All events";
    case kCFStepTriggered:
      return "Triggered";
    case kCFStepVertex:
      return "Primary Vertex";
    case kCFStepAnaTopology:
      return "Required analysis topology";
    case kCFStepTrackedOnlyPrim:
      return "Tracked (matched MC, only primaries)";
    case kCFStepTracked:
      return "Tracked (matched MC, all)";
    case kCFStepReconstructed:
      return "Reconstructed";
    case kCFStepRealLeading:
      return "Correct leading particle identified";
    case kCFStepBiasStudy:
      return "Bias study applying tracking efficiency";
    case kCFStepBiasStudy2:
      return "Bias study applying tracking efficiency in two steps";
  }
  
  return 0;
}

//____________________________________________________________________
void AliUEHist::CopyReconstructedData(AliUEHist* from)
{
  // copies those histograms extracted from ESD to this object
  
  // TODO at present only the pointers are copied
  
  for (Int_t region=0; region<4; region++)
  {
    if (!fTrackHist[region])
      continue;
  
    fTrackHist[region]->SetGrid(AliUEHist::kCFStepReconstructed, from->fTrackHist[region]->GetGrid(AliUEHist::kCFStepReconstructed));
    //fTrackHist[region]->SetGrid(AliUEHist::kCFStepTrackedOnlyPrim, from->fTrackHist[region]->GetGrid(AliUEHist::kCFStepTrackedOnlyPrim));
    fTrackHist[region]->SetGrid(AliUEHist::kCFStepBiasStudy,     from->fTrackHist[region]->GetGrid(AliUEHist::kCFStepBiasStudy));
  }
    
  fEventHist->SetGrid(AliUEHist::kCFStepReconstructed, from->fEventHist->GetGrid(AliUEHist::kCFStepReconstructed));
  //fEventHist->SetGrid(AliUEHist::kCFStepTrackedOnlyPrim, from->fEventHist->GetGrid(AliUEHist::kCFStepTrackedOnlyPrim));
  fEventHist->SetGrid(AliUEHist::kCFStepBiasStudy,     from->fEventHist->GetGrid(AliUEHist::kCFStepBiasStudy));
}

//____________________________________________________________________
void AliUEHist::ExtendTrackingEfficiency(Bool_t verbose)
{
  // fits the tracking efficiency at high pT with a constant and fills all bins with this tracking efficiency

  Float_t fitRangeBegin = 5.01;
  Float_t fitRangeEnd = 14.99;
  Float_t extendRangeBegin = 10.01;

  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;
    Float_t centralityBins[] = { 0, 10, 20, 40, 60, 100 };
    Int_t nCentralityBins = 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 = nCentralityBins;
      
      if (caseNo == 0)
      {
        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++)
      {
	if (centralityBinsLocal == 1)
	  SetCentralityRange(centralityBins[0] + 0.1, centralityBins[nCentralityBins] - 0.1);
	else
	  SetCentralityRange(centralityBins[i] + 0.1, centralityBins[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 = 0;
	    while (centralityBins[z2Bin+1] < fTrackHistEfficiency->GetAxis(3, 0)->GetBinCenter(z2))
	      z2Bin++;
	    
            //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)
{
  // corrects the dphi distribution with an extra factor close to dphi ~ 0

  Printf("WARNING: In AliUEHist::AdditionalDPhiCorrection.");

  THnSparse* grid = fTrackHist[0]->GetGrid(step)->GetGrid();
  
  // optimized implementation
  for (Int_t binIdx = 0; binIdx < grid->GetNbins(); binIdx++)
  {
    Int_t bins[6];
    Double_t value = grid->GetBinContent(binIdx, bins);
    Double_t error = grid->GetBinError(binIdx);
    
    Float_t binCenter = grid->GetAxis(4)->GetBinCenter(bins[4]);
    if (TMath::Abs(binCenter) < 0.2)
    {
      value *= 0.985;
      error *= 0.985;
    }
    else if (TMath::Abs(binCenter) < 0.3)
    {
      value *= 0.9925;
      error *= 0.9925;
    }
    
    grid->SetBinContent(bins, value);
    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();
}